Adding Bi-CGSTAB to the linear solver list.

Author: ctkelley

NDA_Notebook.ipynb

@@ -3,4 +3,5 @@ using LinearAlgebra
 using BenchmarkTools
 using PyPlot
 push!(LOAD_PATH,"./src")
+push!(LOAD_PATH,".")
 using NDA_QMC

nda_note.jl

@@ -74,100 +74,128 @@
 
 <tr valign="top">
 <td align="right" class="bibtexnumber">
-[<a name="knoll08">KPS08</a>]
+[<a name="ctk:fajulia">Kel20c</a>]
 </td>
 <td class="bibtexitem">
-D.&nbsp;A. Knoll, H.&nbsp;Park, and K.&nbsp;Smith.
- A new look at nonlinear acceleration.
- <em>Nuclear Science and Engineering</em>, 99:332-334, 2008.
+C.&nbsp;T. Kelley.
+ Solving Nonlinear Equations with Iterative Methods: Solvers and
+  Examples in Julia, 2020.
+ Unpublished book ms, under contract with SIAM.
 
 </td>
 </tr>
 
 
 <tr valign="top">
 <td align="right" class="bibtexnumber">
-[<a name="knollnda">KPS11</a>]
+[<a name="ctk:notebooknl">Kel20a</a>]
+</td>
+<td class="bibtexitem">
+C.&nbsp;T. Kelley.
+ Notebook for Solving Nonlinear Equations with Iterative Methods:
+  Solvers and Examples in Julia.
+ https://github.com/ctkelley/NotebookSIAMFANL, 2020.
+ IJulia Notebook.
+[&nbsp;<a href="http://dx.doi.org/10.5281/zenodo.4284687">DOI</a>&nbsp;| 
+<a href="https://github.com/ctkelley/NotebookSIAMFANL">http</a>&nbsp;]
+
+</td>
+</tr>
+
+
+<tr valign="top">
+<td align="right" class="bibtexnumber">
+[<a name="ctk:siamfanl">Kel20b</a>]
+</td>
+<td class="bibtexitem">
+C.&nbsp;T. Kelley.
+ SIAMFANLEquations.jl.
+ https://github.com/ctkelley/SIAMFANLEquations.jl, 2020.
+ Julia Package.
+[&nbsp;<a href="http://dx.doi.org/10.5281/zenodo.4284807">DOI</a>&nbsp;| 
+<a href="https://github.com/ctkelley/SIAMFANLEquations.jl">http</a>&nbsp;]
+
+</td>
+</tr>
+
+
+<tr valign="top">
+<td align="right" class="bibtexnumber">
+[<a name="knoll08">KPS08</a>]
 </td>
 <td class="bibtexitem">
 D.&nbsp;A. Knoll, H.&nbsp;Park, and K.&nbsp;Smith.
- Application of the Jacobian-free Newton-Krylov method to
-  nonlinear acceleration of transport source iteration in slab geometry.
- <em>Nuclear Science and Engineering</em>, 167:122-132, 2011.
+ A new look at nonlinear acceleration.
+ <em>Nuclear Science and Engineering</em>, 99:332-334, 2008.
 
 </td>
 </tr>
 
 
 <tr valign="top">
 <td align="right" class="bibtexnumber">
-[<a name="ctk:jeff1">WKKP13b</a>]
+[<a name="knollnda">KPS11</a>]
 </td>
 <td class="bibtexitem">
-Jeff Willert, C.&nbsp;T. Kelley, D.&nbsp;A. Knoll, and H.&nbsp;K. Park.
- Hybrid deterministic/Monte Carlo neutronics.
- 35:S62-S83, 2013.
+D.&nbsp;A. Knoll, H.&nbsp;Park, and K.&nbsp;Smith.
+ Application of the Jacobian-free Newton-Krylov method to
+  nonlinear acceleration of transport source iteration in slab geometry.
+ <em>Nuclear Science and Engineering</em>, 167:122-132, 2011.
 
 </td>
 </tr>
 
 
 <tr valign="top">
 <td align="right" class="bibtexnumber">
-[<a name="ctk:jeff2">WCK15</a>]
+[<a name="gmres">SS86</a>]
 </td>
 <td class="bibtexitem">
-Jeffrey Willert, Xiaojun Chen, and C.&nbsp;T. Kelley.
- Newton's method for Monte Carlo-based residuals.
- <em>SIAM J. Numer. Anal.</em>, 53:1738-1757, 2015.
-[&nbsp;<a href="http://dx.doi.org/10.1137/130905691">DOI</a>&nbsp;]
+Y.&nbsp;Saad and M.H. Schultz.
+ GMRES a generalized minimal residual algorithm for solving
+  nonsymmetric linear systems.
+ <em>SIAM J. Sci. Stat. Comp.</em>, 7:856-869, 1986.
 
 </td>
 </tr>
 
 
 <tr valign="top">
 <td align="right" class="bibtexnumber">
-[<a name="ctk:fajulia">Kel20c</a>]
+[<a name="bicgstab">vdV92</a>]
 </td>
 <td class="bibtexitem">
-C.&nbsp;T. Kelley.
- Solving Nonlinear Equations with Iterative Methods: Solvers and
-  Examples in Julia, 2020.
- Unpublished book ms, under contract with SIAM.
+H.&nbsp;A. van&nbsp;der Vorst.
+ Bi-CGSTAB: A fast and smoothly converging variant to Bi-CG for
+  the solution of nonsymmetric systems.
+ 13:631-644, 1992.
 
 </td>
 </tr>
 
 
 <tr valign="top">
 <td align="right" class="bibtexnumber">
-[<a name="ctk:notebooknl">Kel20a</a>]
+[<a name="ctk:jeff1">WKKP13b</a>]
 </td>
 <td class="bibtexitem">
-C.&nbsp;T. Kelley.
- Notebook for Solving Nonlinear Equations with Iterative Methods:
-  Solvers and Examples in Julia.
- https://github.com/ctkelley/NotebookSIAMFANL, 2020.
- IJulia Notebook.
-[&nbsp;<a href="http://dx.doi.org/10.5281/zenodo.4284687">DOI</a>&nbsp;| 
-<a href="https://github.com/ctkelley/NotebookSIAMFANL">http</a>&nbsp;]
+Jeff Willert, C.&nbsp;T. Kelley, D.&nbsp;A. Knoll, and H.&nbsp;K. Park.
+ Hybrid deterministic/Monte Carlo neutronics.
+ <em>SIAM J. Sci. Comp.</em>, 35:S62-S83, 2013.
 
 </td>
 </tr>
 
 
 <tr valign="top">
 <td align="right" class="bibtexnumber">
-[<a name="ctk:siamfanl">Kel20b</a>]
+[<a name="ctk:jeff2">WCK15</a>]
 </td>
 <td class="bibtexitem">
-C.&nbsp;T. Kelley.
- SIAMFANLEquations.jl.
- https://github.com/ctkelley/SIAMFANLEquations.jl, 2020.
- Julia Package.
-[&nbsp;<a href="http://dx.doi.org/10.5281/zenodo.4284807">DOI</a>&nbsp;| 
-<a href="https://github.com/ctkelley/SIAMFANLEquations.jl">http</a>&nbsp;]
+Jeffrey Willert, Xiaojun Chen, and C.&nbsp;T. Kelley.
+ Newton's method for Monte Carlo-based residuals.
+ <em>SIAM J. Numer. Anal.</em>, 53:1738-1757, 2015.
+[&nbsp;<a href="http://dx.doi.org/10.1137/130905691">DOI</a>&nbsp;]
 
 </td>
 </tr>

ndaqmc.html

@@ -30,9 +30,13 @@ export readdata
 export sn_validate
 export qmc_sweep
 export qmc_init
-export qmc_gmres
+export qmc_krylov
 export qmc_si
 export tab_test
+export nda_fixed
+export nda_init
+export dhateval!
+export ctk_qmc_nda_test
 
 include("validate.jl")
 include("sn_angles.jl")
@@ -47,6 +51,7 @@ include("nda.jl")
 include("nda_iteration.jl")
 include("siewert.jl")
 include("Tim_QMC/ctk_qmc_init.jl")
-include("Tim_QMC/gmres_test.jl")
+include("Tim_QMC/qmc_krylov.jl")
 include("Tim_QMC/ctk_qmc_sweep.jl")
+include("Tim_QMC/ctk_qmc_nda_test.jl")
 end

src/NDA_QMC.jl

@@ -0,0 +1,78 @@
+function ctk_qmc_nda_test(N,Nx)
+s=1.0
+na=20
+nx=Nx+1
+nda_data=nda_init(nx,na,s)
+sn_data=nda_data.sn_data
+angles=sn_data.angles;
+nx=sn_data.nx;
+dx=sn_data.dx;
+na=length(angles);
+psi=sn_data.psi;
+phi=zeros(nx,)
+phix=copy(phi)
+phiout=flux_map!(phix,sn_data)
+psix=copy(psi); phiy=copy(phi)
+psiout=transport_sweep!(psix, phiy, sn_data; phiedge=true)
+dhat=zeros(nx-1,)
+classout=fluxdata(psiout,sn_data)
+qmc_data=qmc_init(N, Nx, 10, 1.0)
+phiq=nda_data.AV*phi
+qout=qmc_sweep(phiq, qmc_data)
+phip=nda_data.AV*qout.phi_edge
+phix=nda_data.AV*phiout
+println("L infty absolute errors")
+println("Flux edge difference = ", norm(qout.phi_edge-phiout,Inf))
+println("Flux average difference = ", norm(phip-phix,Inf))
+println("Current edge difference = ",norm(classout.hocurrent-qout.J_edge,Inf))
+Jx=nda_data.AV*classout.hocurrent
+println("Current average difference = ",norm(Jx-qout.J_avg,Inf))
+println("dphi difference  = ", norm(3.0*classout.dflux-qout.dphi,Inf))
+# relative l infty errors
+println("L infty relative errors")
+println("Flux edge difference = ", relerr(qout.phi_edge,phiout,Inf))
+println("Flux average difference = ", relerr(phip,phix,Inf))
+println("Current edge difference = ",relerr(classout.hocurrent,qout.J_edge,Inf))
+println("Current average difference = ",relerr(Jx,qout.J_avg,Inf))
+println("dphi difference  = ", relerr(3.0*classout.dflux,qout.dphi,Inf))
+# relative l2 errors
+println("L2 relative errors")
+println("Flux edge difference = ", relerr(qout.phi_edge,phiout))
+println("Flux average difference = ", relerr(phip,phix))
+println("Current edge difference = ", relerr(classout.hocurrent,qout.J_edge))
+println("Current average difference = ",relerr(Jx,qout.J_avg))
+println("dphi difference  = ", relerr(3.0*classout.dflux,qout.dphi))
+end
+
+function relerr(xf,yf,p=2)
+df=xf-yf
+relerr=norm(df,p)/norm(yf,p)
+return relerr
+end
+
+function fluxdata(psi,sn_data)
+weights=sn_data.weights;
+angles=sn_data.angles;
+nx=sn_data.nx;
+dx=sn_data.dx;
+hoflux=zeros(nx,);
+hocurrent=zeros(nx,);
+pa=weights.*angles;
+hoflux=(weights' * psi)';
+bcl=hoflux[1];
+bcr=hoflux[nx];
+hocurrent=(pa' * psi)';
+dflux=zeros(nx-1,)
+dhat=zeros(nx-1,)
+phiav=zeros(nx-1,)
+jav=zeros(nx-1,)
+dflux.=(hoflux[2:nx]-hoflux[1:nx-1])
+dflux./= (3.0*dx)
+phiav.=.5*(hoflux[2:nx] .+ hoflux[1:nx-1])
+jav.=.5*(hocurrent[2:nx] .+ hocurrent[1:nx-1])
+dhat.=(jav + dflux)./phiav
+return (dhat=dhat, bcl=bcl, bcr=bcr, dflux=dflux, 
+       hoflux=hoflux, hocurrent=hocurrent)
+end
+
+

src/Tim_QMC/.ctk_qmc_nda_test.jl.swp

@@ -170,7 +170,7 @@ function qmc_sweep(phi_avg, qmc_data)
         dphi[i] = phi_edge[i+1] - phi_edge[i]
     end
 
-    dphi /= (dx)
+    dphi ./= (dx)
     exit_left_bins[:,2] /= dmu
     exit_right_bins[:,2] /= dmu
 

src/Tim_QMC/ctk_qmc_nda_test.jl

@@ -1,15 +1,31 @@
 """
-qmc_gmres(N=10^4, Nx=100, na2=11; s=1.0)
+qmc_krylov(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)
+function qmc_krylov(N=10^3, Nx=100, na2=11; s=1.0, tol=1.e-10, onlygmres=false)
 mdata=mdata_init(N, Nx, na2, s)
 phi0=zeros(Nx,);
 b=mdata.frhs;
 V=zeros(Nx,20)
 eta=tol
-gmout=kl_gmres(phi0,b,axb,V,eta; pdata=mdata);
-return gmout
+#
+# kl_gmres and kl_bicgstab solve the problem.
+#
+gout = kl_gmres(phi0, b, axb, V, tol; pdata = mdata)
+if onlygmres
+    bout=(solb=[], reshist=[])
+    else
+    bout = kl_bicgstab(phi0, b, axb, V, tol; pdata = mdata)
+end
+#gmout=kl_gmres(phi0,b,axb,V,eta; pdata=mdata);
+sol=gout.sol
+solb=bout.sol
+reshistg=gout.reshist
+reshistb=bout.reshist
+    #
+    # Tabulate the exit distributions to check results.
+    #
+    return(sol=sol, solb=solb, reshistg=reshistg, reshistb=reshistb)
 end
 
 
@@ -18,7 +34,7 @@ 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)
+itout=qmc_krylov(N, Nx, na2; s=s, tol=tol)
 qmctab=sn_tabulate(s, Nx, itout.sol; maketab=false, phiedge=false)
 return [qmctab.left qmctab.right]
 end

src/Tim_QMC/ctk_qmc_sweep.jl

@@ -1,3 +1,8 @@
+"""
+compare(s = 1.0; qmc=false)
+
+Makes the Krylov comparison tables/figures
+"""
 function compare(s = 1.0; qmc=false)
     na2 = 20
     N=1000
@@ -22,10 +27,10 @@ function compare(s = 1.0; qmc=false)
     fout = sout.flux
     end
     #
-    # and now GMRES and BiCGSTAB
+    # and now GMRES and Bi-CGSTAB
     #
     if qmc
-    gout=qmc_gmres(N, nx, na2; s=s)
+    gout=qmc_krylov(N, nx, na2; s=s)
     else
     gout=krylov_iteration(sn_data,s)
     end
@@ -40,6 +45,7 @@ function compare(s = 1.0; qmc=false)
     glen=collect(1:length(ghist))
     bhist = gout.reshistb ./ gout.reshistb[1]
     blen=2*collect(1:length(bhist))
+    blen[1]=1
     semilogy(glen,ghist, "k-", slen,snhist, "k--", blen, bhist, "k-.")
     legend(["gmres", "source iteration", "bicgstab"])
     xlabel("Transport Sweeps")
@@ -59,5 +65,5 @@ function compare(s = 1.0; qmc=false)
     (solverdelta < 1.e-6) || println(solverdelta)
     sn_tabulate(s,nx,fgm)
     println("Norm of result differences: GMRES = ", norm(fout - fgm, Inf),
-            ",    BiCGSTAB = ",norm(fout - fgmb, Inf))
+            ",    Bi-CGSTAB = ",norm(fout - fgmb, Inf))
 end

src/Tim_QMC/gmres_test.jl src/Tim_QMC/qmc_krylov.jl

@@ -40,9 +40,6 @@ function krylov_iteration(sn_data,s,tol=1.e-8; onlygmres=false)
     solb = bout.sol
     reshistg=gout.reshist
     reshistb=bout.reshist
-    #
-    # Tabulate the exit distributions to check results.
-    #
     return(sol=sol, solb=solb, reshistg=reshistg, reshistb=reshistb)
 end