Updates and bug fixes

Author: ctkelley

ctk_nda_src/Tim_CMFD/cmfd.jl

@@ -22,13 +22,14 @@ end
 
 
 """
-cmfd_krylov(N=4096*4, Nx=512, Nc=64, na2=11; s=1.0, tol=1.e-10)
+cmfd_krylov(N=4096*4, Nx=512, Nc=64, na2=11; s=1.0, tol=1.e-10, tau=5.0)
 Solves the linear system you get from CMFD with GMRES.
 """
-function cmfd_krylov(N=4096*4, Nx=512, Nc=64, na2=11; s=1.0, tol=1.e-10)
+function cmfd_krylov(N=4096*4, Nx=512, Nc=64, na2=11; s=1.0, 
+                  tol=1.e-10, tau=5.0)
 phi_c = zeros(Nc)
 phi_c0 = zeros(Nc)
-cmfd_data=cmfd_init(N, Nx, Nc, na2, s)
+cmfd_data=cmfd_init(N, Nx, Nc, na2, s; tau=tau)
 linop_data=linop_init(Nc, cmfd_data)
 b=linop_data.frhs
 V=zeros(Nc,20)
@@ -91,14 +92,14 @@ return Mprod
 end
 
 
-function cmfd_init(N, Nx, Nc, na2, s)
+function cmfd_init(N, Nx, Nc, na2, s; tau=5.0)
 phi_f = zeros(Nx)
 J_c = zeros(Nc)
-sn_data=sn_init(Nx+1, 2*na2, s)
-qmc_data = qmc_init(N, Nx, na2, s);
-nda_data=nda_cmfd_init(Nx, Nc, 2*na2, s)
+sn_data=sn_init(Nx+1, 2*na2, s; tau=tau)
+qmc_data = qmc_init(N, Nx, na2, s; Lx = tau);
+nda_data=nda_cmfd_init(Nx, Nc, 2*na2, s; tau=tau)
 return(phi_f=phi_f, J_c = J_c, qmc_data=qmc_data, sn_data=nda_data.sn_data, 
-       nda_data=nda_data)
+       nda_data=nda_data, tau=tau)
 end
 
 function linop_init(Nc, cmfd_data)

ctk_nda_src/Tim_CMFD/compare_cmfd.jl

@@ -1,17 +1,24 @@
-function compare_cmfd(N=4096*4, Nx=1024, Nc=64, na2=11; s=1.0, tol=1.e-10)
+function compare_cmfd(N=4096*4, Nx=1024, Nc=64, na2=11; 
+                    tau=5.0, s=1.0, tol=1.e-10, plotme=true)
 phi_c = zeros(Nc)
 phi_c0 = zeros(Nc)
-cmfd_data=cmfd_init(N, Nx, Nc, na2, s)
+dx=tau/Nc
+x=dx*.5:dx:tau-dx*.5
+cmfd_data=cmfd_init(N, Nx, Nc, na2, s; tau=tau)
 #
 # QMC-GMRES is the baseline 
 #
-goutgm=qmc_krylov(N, Nx, na2; s=s)
+goutgm=qmc_krylov(N, Nx, na2; s=s, tau=tau)
+solgmf=goutgm.sol
+solgm=zeros(Nc)
+solgm.=ftoc(solgmf,solgm)
 reshistg=goutgm.reshistg
 itgmax=length(reshistg)
 #
 # Newton-GMRES + NDA
 #
-nout=cmfd_nsoli(Nc; s=s)
+nout=cmfd_nsoli(N, Nx, Nc; s=s, tau=tau)
+soln=nout.solution
 fout=nout.stats.ifun+nout.stats.ijac+nout.stats.iarm
 #
 # count the warm-up qmc-si sweep
@@ -22,18 +29,32 @@ reshistn=nout.history
 #
 # NDA-Picard
 #
-pout=test_cmfd(Nc; s=s)
+pout=test_cmfd(N, Nx, Nc; s=s, tau=tau)
+solp=pout.solution
+println(norm(soln-solp,Inf),"  ",norm(soln-solgm,Inf))
 reshistp=pout.history
 itpic=length(reshistp)
 #
 # Plot the results
 #
+Itau=Int(tau)
+if plotme
+figure(1)
+plot(x,soln,"k-",x,solgm,"k--",x,solp,"k-.")
+figure(2)
 semilogy(1:itgmax,reshistg/reshistg[1],"k-",
                1:itpic,reshistp/reshistp[1],"k.-",
                ftot,reshistn/reshistn[1],"k--")
-title("N=$N, Nx=$Nx, Nc=$Nc")
+stri=L"\infty"
+xxst="N=$N, Nx=$Nx, Nc=$Nc, s="*stri*", Lx = $Itau"
+if s==Inf
+#title("N=$N, Nx=$Nx, Nc=$Nc, s=, Lx = $Itau")
+title(xxst)
+else
+title("N=$N, Nx=$Nx, Nc=$Nc, s=1.0, Lx = $Itau")
+end
 xlabel("Transport sweeps")
 legend(["QMC-GMRES","NDA-Picard","NDA-Newton"])
-
+end
 end
 

ctk_nda_src/Tim_CMFD/nda_cmfd.jl

@@ -7,12 +7,13 @@ FS.=phi-cmfd_fixed(phi,cmfd_data)
 return FS 
 end
 
-function cmfd_nsoli(Nc=64; s=1.0)
+function cmfd_nsoli(N=4096*4, Nx=1024, Nc=64; s=1.0, tau=5.0)
 FS=zeros(Nc,)
 FPS=zeros(Nc,20)
 phi0=zeros(Nc,)
-N=4096*4; Nx=1024; na2=11;
-cmfd_data=cmfd_init(N, Nx, Nc, na2, s);
+#N=4096*4; Nx=1024; 
+na2=11;
+cmfd_data=cmfd_init(N, Nx, Nc, na2, s; tau=tau);
 #
 # Fix the initial iterate if you want decent results.
 #
@@ -42,7 +43,8 @@ function cmfd_fixed(phi_in, cmfd_data)
 # Solve high order problem
 #
 qmc_data=cmfd_data.qmc_data
-nc=length(phi_in); dx=5.0/nc
+tau=cmfd_data.tau
+nc=length(phi_in); dx=tau/nc
 #
 # we get the cell average flux, so convert to cell edges
 #
@@ -83,8 +85,8 @@ phiout=copy(phi)
 return phiout
 end
 
-function nda_cmfd_init(nx, nc, na2, s)
-sn_data=sn_init(nc+1, na2, s)
+function nda_cmfd_init(nx, nc, na2, s; tau=5.0)
+sn_data=sn_init(nc+1, na2, s; tau=tau)
 dx=sn_data.dx
 c=sn_data.c
 nout=cmfdD2(nc+1,dx,c)

ctk_nda_src/Tim_CMFD/result_quality.jl

@@ -1,13 +1,13 @@
-function result_quality(Nc = 64)
+function result_quality(Nc = 64; tau=5.0)
 na2=11
 Nx=512
-dx=5.0/Float64(Nc)
+dx=tau/Float64(Nc)
 x=.5*dx:dx:5.0-.5*dx
 sol=zeros(Nc,3)
 N=[1024, 2048, 4096]
 #N *= 4
 for in=1:3
-nout=cmfd_krylov(N[in], Nx, Nc, na2)
+nout=cmfd_krylov(N[in], Nx, Nc, na2; tau=tau)
 sol[:,in]=nout.sol
 end
 plot(x,sol[:,1],"k-",x,sol[:,2],"k--",x,sol[:,3],"k-.")

ctk_nda_src/Tim_CMFD/test_cmfd.jl

@@ -1,15 +1,16 @@
-function test_cmfd(Nc=64; s=1.0)
-N=4096*4; Nx=1024; na2=11; 
+function test_cmfd(N=4096*4, Nx=1024, Nc=64; s=1.0, tau=5.0)
+#N=4096*4; Nx=1024; 
+na2=11; 
 #N=4096*2; Nx=1024; na2=11;
-tau=5.0; dx=tau/Nc; x=.5*dx:dx:tau-.5*dx
-cmfd_data=cmfd_init(N, Nx, Nc, na2, s);
+dx=tau/Nc; x=.5*dx:dx:tau-.5*dx
+cmfd_data=cmfd_init(N, Nx, Nc, na2, s; tau=tau);
 phi0=zeros(Nc);
 phic=copy(phi0)
 reshist=Float64[]
 plotit=false
 itc=0
 resid=1.0
-while (itc < 100) && (resid > 1.e-11)
+while (itc < 200) && (resid > 1.e-11)
 itc+=1
     phip=cmfd_fixed(phic, cmfd_data)
 resid=norm(phip-phic)

ctk_nda_src/Tim_QMC/ctk_qmc_init.jl

@@ -1,8 +1,8 @@
 #using Sobol
 
-function qmc_init(N, Nx, na2, s)
+function qmc_init(N, Nx, na2, s; Lx=5.0)
 
-    Lx = 5.0
+#    Lx = 5.0
     dx = Lx/Nx
     #define tally mesh
     low_edges = range(0, stop=Lx-dx, length=Nx)

ctk_nda_src/Tim_QMC/qmc_krylov.jl

@@ -1,12 +1,13 @@
 """
-qmc_krylov(N=10^4, Nx=100, na2=11; s=1.0)
+qmc_krylov(N=10^4, Nx=100, na2=11; s=1.0, tau=5.0)
 Solves the linear system you get from QMC with GMRES.
 """
-function qmc_krylov(N=1024, Nx=128, na2=11; s=1.0, tol=1.e-10, onlygmres=true)
-mdata=mdata_init(N, Nx, na2, s)
+function qmc_krylov(N=1024, Nx=128, na2=11; 
+          s=1.0, tol=1.e-10, onlygmres=true, tau=5.0)
+mdata=mdata_init(N, Nx, na2, s; tau=tau)
 phi0=zeros(Nx,);
 b=mdata.frhs;
-V=zeros(Nx,20)
+V=zeros(Nx,80)
 eta=tol
 #
 # kl_gmres and kl_bicgstab solve the problem.
@@ -38,11 +39,12 @@ end
 
 
 """
-qmc_si(N=10^4, Nx=100, na2=11; s=1.0, tol=1.e-8, maxit=100)
+qmc_si(N=10^4, Nx=100, na2=11; tau=5.0, 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)
+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; tau=tau)
 phic=zeros(Nx,);
 delflux=1.0
 reshist=[]
@@ -86,9 +88,9 @@ 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)
+function mdata_init(N, Nx, na2, s; tau=5.0)
 # Precomputed data
-qmc_data = qmc_init(N, Nx, na2, s);
+qmc_data = qmc_init(N, Nx, na2, s; Lx=tau);
 # Sweep with zero RHS
 phizero=zeros(Nx,);
 nullout=qmc_sweep(phizero,qmc_data);

ctk_nda_src/qmc_init.jl

@@ -1,8 +1,8 @@
 using Sobol
 
-function qmc_init(N, Nx, na2, s)
+function qmc_init(N, Nx, na2, s; Lx = 5.0)
 
-    Lx = 5.0
+#    Lx = 5.0
     dx = Lx/Nx
     #define tally mesh
     low_edges = range(0, stop=Lx-dx, length=Nx)

ctk_nda_src/sn_init.jl

@@ -14,14 +14,14 @@ s = parameter in Garcia/Siewert
 
 angles, weights = nodes and weights for the angular quadrature
 """
-function sn_init(nx, na2, s; siewert=false, phiedge=true)
+function sn_init(nx, na2, s; siewert=false, phiedge=true, tau=5.0)
     if siewert
     angles = [-.05; collect(-.1:-.1:-1.0); 0.05; collect(0.1:0.1:1.0)]
     weights = angles
     else
     (angles, weights) = sn_angles(na2)
     end
-    tau = 5.0
+#    tau = 5.0
     dx = tau / (nx - 1)
     na = floor(Int, na2 / 2)
     x = collect(0:dx:tau)