Adding BiCGSTAB to the linear solver mix; getting ready for QMC-NDA

Author: ctkelley

NDA_Notebook.ipynb

@@ -1,6 +1,5 @@
 [deps]
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
-JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"

Project.toml

@@ -14,7 +14,7 @@ using Sobol
 export testgauss
 export sn_angles
 export source_iteration
-export gmres_iteration
+export krylov_iteration
 export sn_init
 export sn_tabulate
 export compare
@@ -37,7 +37,7 @@ export tab_test
 include("validate.jl")
 include("sn_angles.jl")
 include("source_iteration.jl")
-include("gmres_iteration.jl")
+include("krylov_iteration.jl")
 include("transport_sweep.jl")
 include("sn_init.jl")
 include("sn_tabulate.jl")

src/Siewert:s=1.dat src/Data_Orig/Siewert:s=1.datsrc/Siewert:s=1.dat renamed to src/Data_Orig/Siewert:s=1.dat

@@ -176,6 +176,7 @@ function qmc_sweep(phi_avg, qmc_data)
 
     return (phi_avg = phi_avg,
             phi_edge = phi_edge,
+            dphi = dphi,
             J_avg = J_avg,
             J_edge = J_edge,
             exit_right_bins = exit_right_bins,

src/Siewert:s=Inf.dat src/Data_Orig/Siewert:s=Inf.datsrc/Siewert:s=Inf.dat renamed to src/Data_Orig/Siewert:s=Inf.dat

@@ -1,5 +1,6 @@
 """
 qmc_gmres(N=10^4, Nx=100, na2=11; s=1.0)
+Solves the linear system you get from QMC with GMRES.
 """
 function qmc_gmres(N=10^3, Nx=100, na2=11; s=1.0, tol=1.e-10)
 mdata=mdata_init(N, Nx, na2, s)
@@ -12,13 +13,21 @@ return gmout
 end
 
 
+"""
+tab_test(N=10^3, Nx=100, na2=11; s=1.0, tol=1.e-8)
+Makes a table to compare to Garcia-Siewert
+"""
 function tab_test(N=10^3, Nx=100, na2=11; s=1.0, tol=1.e-8)
 itout=qmc_gmres(N, Nx, na2; s=s, tol=tol)
 qmctab=sn_tabulate(s, Nx, itout.sol; maketab=false, phiedge=false)
 return [qmctab.left qmctab.right]
 end
 
 
+"""
+qmc_si(N=10^4, Nx=100, na2=11; s=1.0, tol=1.e-8, maxit=100)
+Source iteration using QMC. Nothing magic here.
+"""
 function qmc_si(N=10^4, Nx=100, na2=11; s=1.0, tol=1.e-8, maxit=100)
 qmc_data=qmc_init(N, Nx, na2, s)
 phic=zeros(Nx,);
@@ -37,21 +46,40 @@ return ( sol=phic, history=reshist)
 end
 
 
+"""
+axb(phi,mdata)
+Matrix-vector product to feed to linear solvers.
+"""
 function axb(phi,mdata)
 qmc_data=mdata.qmc_data;
 frhs=mdata.frhs
 matvec=Qlin(phi,qmc_data,frhs)
 end
 
+"""
+Qlin(phi,qmc_data,frhs)
+Does a QMC transport sweep with flux phi, subtracts off the zero bc solution
+to obtain the linear integral operator. The subtract that from the identity. 
+"""
 function Qlin(phi,qmc_data,frhs)
 nullout=qmc_sweep(phi,qmc_data)
 Mprod=phi-(nullout.phi_avg-frhs)
 end
 
+
+"""
+mdata_init(N, Nx, na2, s)
+Collects the precomputed data for QMC and does a sweep with 
+zero boundary data to get the right hand side for the linear
+equation forumlation.
+"""
 function mdata_init(N, Nx, na2, s)
+# Precomputed data
 qmc_data = qmc_init(N, Nx, na2, s);
+# Sweep with zero RHS
 phizero=zeros(Nx,);
 nullout=qmc_sweep(phizero,qmc_data);
 frhs=nullout.phi_avg;
+#
 mdata=(frhs=frhs, qmc_data =qmc_data)
 end

src/NDA_QMC.jl

@@ -22,21 +22,26 @@ function compare(s = 1.0; qmc=false)
     fout = sout.flux
     end
     #
-    # and now GMRES
+    # and now GMRES and BiCGSTAB
     #
     if qmc
     gout=qmc_gmres(N, nx, na2; s=s)
     else
-    gout=gmres_iteration(sn_data,s)
+    gout=krylov_iteration(sn_data,s)
     end
     fgm = gout.sol
+    fgmb = gout.solb
     #
     # Make a nice plot
     #
     snhist = sout.history ./ sout.history[1]
-    ghist = gout.reshist ./ gout.reshist[1]
-    semilogy(ghist, "k-", snhist, "k--")
-    legend(["gmres", "source iteration"])
+    slen=collect(1:length(snhist))
+    ghist = gout.reshistg ./ gout.reshistg[1]
+    glen=collect(1:length(ghist))
+    bhist = gout.reshistb ./ gout.reshistb[1]
+    blen=2*collect(1:length(bhist))
+    semilogy(glen,ghist, "k-", slen,snhist, "k--", blen, bhist, "k-.")
+    legend(["gmres", "source iteration", "bicgstab"])
     xlabel("Transport Sweeps")
     ylabel("Relative Residual")
     s == Inf ? strs = L"\infty" : strs = string(s)
@@ -53,5 +58,6 @@ function compare(s = 1.0; qmc=false)
 #    (solverdelta < 1.e-6) || error("results not the same")
     (solverdelta < 1.e-6) || println(solverdelta)
     sn_tabulate(s,nx,fgm)
-    println("Norm of result difference = ", norm(fout - fgm, Inf))
+    println("Norm of result differences: GMRES = ", norm(fout - fgm, Inf),
+            ",    BiCGSTAB = ",norm(fout - fgmb, Inf))
 end

src/Tim_QMC/ctk_qmc_sweep.jl

@@ -1,6 +1,10 @@
-function gmres_iteration(sn_data,s,tol=1.e-8)
+"""
+krlov_iteration(sn_data,s,tol=1.e-8; onlygmres=false)
+Solves the transport equation with GMRES and BiCGSTAB
+"""
+function krylov_iteration(sn_data,s,tol=1.e-8; onlygmres=false)
     #
-    # Set up the linear system for GMRES
+    # Set up the linear system for the Krylov solvers
     #
     nx=sn_data.nx
     #
@@ -23,16 +27,23 @@ function gmres_iteration(sn_data,s,tol=1.e-8)
     f0 = zeros(size(rhs))
     V = zeros(nx, 20)
     #
-    # kl_gmres solves the problem.
+    # kl_gmres and kl_bicgstab solve the problem.
     #
     gout = kl_gmres(f0, rhs, sn_matvec, V, tol; pdata = sn_data)
+    if onlygmres
+    bout=(solb=[], reshist=[])
+    else
+    bout = kl_bicgstab(f0, rhs, sn_matvec, V, tol; pdata = sn_data)
+    end
     #
     sol = gout.sol
-    reshist=gout.reshist
+    solb = bout.sol
+    reshistg=gout.reshist
+    reshistb=bout.reshist
     #
     # Tabulate the exit distributions to check results.
     #
-    return(sol=sol , reshist=reshist)
+    return(sol=sol, solb=solb, reshistg=reshistg, reshistb=reshistb)
 end
 
 """

src/Tim_QMC/gmres_test.jl

@@ -13,7 +13,7 @@ function siewert(s=1.0; filedump=false, na2=80, nx=16001, phiedge=true)
     #
     # Solve the problem
     #
-    gout=gmres_iteration(sn_data,s)
+    gout=krylov_iteration(sn_data,s; onlygmres=true)
     if phiedge
     fgm = gout.sol
     else