Updates and progress on cmfd

Author: ctkelley

.DS_Store

@@ -11,6 +11,12 @@ using PyPlot
 using Printf
 using Sobol
 
+export sn_flux
+export plot_test
+export cmfd_si
+export cmfd_krylov
+export ftoc
+export ctof
 export ndaerr
 export averop
 export fprintTeX
@@ -58,14 +64,17 @@ include("compare.jl")
 include("nda.jl")
 include("nda_iteration.jl")
 include("siewert.jl")
+include("sn_flux.jl")
 include("Tim_QMC/qmc_nda.jl")
 include("Tim_QMC/ctk_qmc_init.jl")
 include("Tim_QMC/qmc_krylov.jl")
 include("Tim_QMC/ctk_qmc_sweep.jl")
 include("Tim_QMC/ctk_qmc_nda_test.jl")
 include("Tim_QMC/sn_test.jl")
 include("fprintTex.jl")
-include("averop.jl")
+include("Tim_CMFD/averop.jl")
 include("Tim_NDA/sanity.jl")
 include("Tim_NDA/ndaerr.jl")
+include("Tim_CMFD/cmfd.jl")
+include("Tim_CMFD/plot_test.jl")
 end

QMC-NDA-Woes.ipynb

@@ -0,0 +1,31 @@
+function plot_test(N, Nx, Nc)
+tau=5.0; dx=tau/Nc; x=.5*dx:dx:tau-.5*dx;
+z=dx:dx:tau-dx;
+gm_out=cmfd_krylov(N, Nx, Nc)
+dsola=zeros(Nc)
+#
+# cell centered flux
+#
+sol=gm_out.sol
+#
+# derivative of flux at interior cell edges
+#
+@views dsol=(sol[2:Nc]-sol[1:Nc-1])/dx
+#
+# Map to cell centers
+#
+@views dsola[2:Nc-1] = .5*(dsol[1:Nc-2] + dsol[2:Nc-1])
+dsola[1]=2.0 * dsol[1]-dsol[2]
+dsola[Nc]=2.0 * dsol[Nc-1]-dsol[Nc-2]
+#
+# Now do the currents and D-hat
+#
+J_c = gm_out.J_c
+Dhat = ((dsola/(3.0*dx)) + J_c)./sol
+#
+# You can see why we need large N in the plots
+#
+plot(x,sol,"k-", x, dsola, "k--", x, J_c, "k.", x, Dhat, "k-.))
+title("Flux data: N = $N, Nx = $Nx, Nc = $Nc")
+legend([L" \phi", L"d \phi/dx", L"J_c", L"\hat D "])
+end

QMC-NDA_is_Better.ipynb

@@ -0,0 +1,61 @@
+function averop(p,n)
+d=2.0*ones(n)/4.0;
+d[1]=2.0/3.0;
+d[n]=2.0/3.0;
+du=ones(n-1)/4.0;
+du[1]=1.0/3.0;
+dl=ones(n-1)/4.0;
+dl[n-1]=1.0/3.0;
+T=Tridiagonal(dl,d,du);
+A=sparse(T^p);
+return A
+end
+
+
+"""
+ftoc(uf, uc)
+
+fine to coarse intergrid transfer. This is a cell average deal and
+the number of cells must be a power of 2
+"""
+function ftoc(uf, uc)
+nf=length(uf)
+nc=length(uc)
+lnf=Int(log2(nf))
+(2^lnf == nf) || error("number of cells must be power of 2")
+lnc=Int(log2(nc))
+(2^lnc == nc) || error("number of cells must be power of 2")
+ncells = Int(nf/nc)
+for ic=1:nc
+  cbase=(ic - 1) * ncells + 1
+@views  uc[ic]=sum(uf[cbase:cbase+ncells-1])
+end
+uc ./= ncells
+return uc
+end
+
+
+"""
+ctof(uc,uf)
+
+coarse to fine intergrid transfer. This is a cell average deal and
+the number of cells must be a power of 2
+"""
+function ctof(uc,uf)
+nf=length(uf)
+nc=length(uc)
+lnf=Int(log2(nf))
+(2^lnf == nf) || error("number of cells must be power of 2")
+lnc=Int(log2(nc))
+(2^lnc == nc) || error("number of cells must be power of 2")
+ncells = Int(nf/nc)
+for ic=1:nc
+  cbase=(ic - 1) * ncells + 1
+  uf[cbase:cbase+ncells-1] .= uc[ic]
+end
+return uf
+end
+
+
+
+

ctk_nda_src/Tim_NDA/.ndaerr.jl.swp ctk_nda_src/.nda.jl (Mac-mini.attlocal.net's conflicted copy 2023-03-26).swp

@@ -0,0 +1,8 @@
+function avplot(N=4096, Nx=512, Nc=64, na2=11; s=1.0, tol=1.e-10)
+kout=qmc_krylov(N, Nx, na2; s=s, tol=tol)
+uf=kout.sol
+uc=zeros(Nc)
+ftoc(uf,uc)
+plot(uc)
+return uc
+end

ctk_nda_src/NDA_QMC.jl

@@ -0,0 +1,102 @@
+function cmfd_si(N=4096*4, Nx=512, Nc=64, na2=11; s=1.0, tol=1.e-10)
+outputok=true
+phi_c = zeros(Nc)
+phi_c0 = zeros(Nc)
+cmfd_data=cmfd_init(N, Nx, Nc, na2, s)
+linop_data=linop_init(Nc, cmfd_data)
+reshist=Float64[]
+resid=1.0
+for iz=1:20
+   (phi_c, J_c) = cmfd_sweep(phi_c, cmfd_data)
+   resid=norm(phi_c - phi_c0, Inf)
+   push!(reshist, resid)
+   phi_c0 .= phi_c
+end
+if outputok == true
+lin_resid=cmfd_matvec(phi_c, linop_data) - linop_data.frhs
+println(norm(lin_resid,Inf))
+plot(phi_c)
+end
+return (sol=phi_c, reshist=reshist)
+end
+
+
+"""
+cmfd_krylov(N=4096*4, Nx=512, Nc=64, na2=11; s=1.0, tol=1.e-10)
+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)
+phi_c = zeros(Nc)
+phi_c0 = zeros(Nc)
+cmfd_data=cmfd_init(N, Nx, Nc, na2, s)
+linop_data=linop_init(Nc, cmfd_data)
+b=linop_data.frhs
+V=zeros(Nc,20)
+eta=tol
+#
+# kl_gmres and kl_bicgstab solve the problem.
+#
+gout = kl_gmres(phi_c, b, cmfd_matvec, V, tol; pdata = linop_data)
+sol=gout.sol
+reshistg=gout.reshist
+(phi, J_c) = cmfd_sweep(sol, cmfd_data)
+return(sol=sol, J_c=J_c, reshistg=reshistg)
+end
+
+"""
+cmfd_sweep(phi_c, cmfd_data)
+
+Returns coarse mesh cell average flux and current
+"""
+function cmfd_sweep(phi_c, cmfd_data)
+#
+# Map coarse mesh flux to fine mesh
+#
+phi_f = cmfd_data.phi_f
+ctof(phi_c,phi_f)
+#
+# Do Sam's QMC sweep
+#
+qmc_data=cmfd_data.qmc_data
+qmc_out=qmc_sweep(phi_f, qmc_data)
+phi_f = qmc_out.phi_avg
+#
+# Map output flux and current to coarse mesh
+#
+phi_c = ftoc(phi_f, phi_c)
+J_f = qmc_out.J_avg
+J_c = cmfd_data.J_c
+J_c = ftoc(J_f, J_c)
+return (phi_c, J_c)
+end
+
+"""
+cmfd_matvec(phi,linop_data)
+Matvec for CMFD-QMC
+Does a QMC transport sweep with flux phi, subtracts off the zero bc solution
+to obtain the linear integral operator. Then subtract that from the identity.
+"""
+function cmfd_matvec(phic,linop_data)
+phi_copy=linop_data.phi_copy
+phi_copy .= phic
+frhs=linop_data.frhs
+cmfd_data=linop_data.cmfd_data
+(phi_c, J_c) =cmfd_sweep(phi_copy,cmfd_data)
+Mprod=phic-(phi_c - frhs)
+return Mprod
+end
+
+
+function cmfd_init(N, Nx, Nc, na2, s)
+phi_f = zeros(Nx)
+J_c = zeros(Nc)
+qmc_data = qmc_init(N, Nx, na2, s);
+return(phi_f=phi_f, J_c = J_c, qmc_data=qmc_data)
+end
+
+function linop_init(Nc, cmfd_data)
+nullin=zeros(Nc)
+phi_copy=zeros(Nc)
+(frhs, Jrhs)=cmfd_sweep(nullin, cmfd_data)
+return (frhs = frhs, cmfd_data = cmfd_data, phi_copy = phi_copy)
+end

ctk_nda_src/Tim_CMFD/'

@@ -0,0 +1,31 @@
+function plot_test(N, Nx, Nc)
+tau=5.0; dx=tau/Nc; x=.5*dx:dx:tau-.5*dx;
+z=dx:dx:tau-dx;
+gm_out=cmfd_krylov(N, Nx, Nc)
+dsola=zeros(Nc)
+#
+# cell centered flux
+#
+sol=gm_out.sol
+#
+# derivative of flux at interior cell edges
+#
+@views dsol=(sol[2:Nc]-sol[1:Nc-1])/dx
+#
+# Map to cell centers
+#
+@views dsola[2:Nc-1] = .5*(dsol[1:Nc-2] + dsol[2:Nc-1])
+dsola[1]=2.0 * dsol[1]-dsol[2]
+dsola[Nc]=2.0 * dsol[Nc-1]-dsol[Nc-2]
+#
+# Now do the currents and D-hat
+#
+J_c = gm_out.J_c
+Dhat = ((dsola/3.0) + J_c)./sol
+#
+# You can see why we need large N in the plots
+#
+plot(x,sol,"k-", x, dsola, "k--", x, J_c, "k.", x, Dhat, "k-.")
+title("Flux data: N = $N, Nx = $Nx, Nc = $Nc")
+legend([L" \phi", L"d \phi/dx", L"J_c", L"\hat{D}"])
+end

ctk_nda_src/Tim_CMFD/averop.jl

@@ -1,4 +1,4 @@
-function sanity(Nx=2048, na=40, s=1.0)
+function sanity(Nx=2048, na=40, s=1.0, N=10^4)
 #
 # NDA and NDA-Newton
 #
@@ -23,5 +23,41 @@ snvsnda=kout.sol-ndaout.sol
 snndadiff=norm(snvsnda,Inf)
 snndaok=(snndadiff<2.e-6)
 snndaok || println("SN vs NDA consistency error")
+#
+# Compare SN and QMC
+#
+qout=qmc_krylov(N, Nx, na)
+qmdiff=qout.sol-kout.sol
+println(norm(qmdiff,2)/sqrt(Float64(Nx)))
+#
+# Now see if there is any hope for a correct current
+#
+qmc_data=qmc_init(N, Nx, na, s)
+qdout=qmc_sweep(kout.sol,qmc_data)
+qediff=qdout.phi_avg-kout.sol
+subplot(121)
+semilogy(qediff)
+println(norm(qediff,2)/sqrt(Float64(Nx)))
+(flux, current) = getJ(kout.sol, sn_data)
+cdiff=qdout.J_avg - current
+println(norm(cdiff,2)/sqrt(Float64(Nx)))
+subplot(122)
+semilogy(cdiff)
+#
 return (snndaok && snok && ndaok)
 end
+
+function getJ(phi,sn_data)
+    psi = sn_data.psi
+    psi = transport_sweep!(psi, phi, sn_data)
+    #
+    # Take the 0th moment to get the flux.
+    #
+    weights = sn_data.weights
+    angles = sn_data.angles
+    flux=copy(phi)
+    flux .= (weights' * psi)'
+    current = copy(phi)
+    current .= ((weights.*angles)' * psi)'
+    return (flux=flux, current=current)
+end

ctk_nda_src/Tim_CMFD/avplot.jl

@@ -10,7 +10,7 @@ function qmc_init(N, Nx, na2, s)
     dxs = high_edges - low_edges
     midpoints = 0.5*(high_edges + low_edges)
     edges = range(0, stop=Lx, length=Nx+1)
-    p=2
+    p=4
 #    AV=[]
     AV=averop(p,Nx+1)
     #define angular flux mesh

ctk_nda_src/Tim_CMFD/cmfd.jl

@@ -2,7 +2,7 @@
 qmc_krylov(N=10^4, Nx=100, na2=11; s=1.0)
 Solves the linear system you get from QMC with GMRES.
 """
-function qmc_krylov(N=10^3, Nx=100, na2=11; s=1.0, tol=1.e-10, onlygmres=true)
+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)
 phi0=zeros(Nx,);
 b=mdata.frhs;

ctk_nda_src/Tim_CMFD/plot_test.jl

@@ -57,7 +57,7 @@ dphi[id]=(aphi[id+1]-aphi[id])/dx
 end
 #return (phi_edge=phi_edge, dphi=dphi)
 phi_edge .= aphi
-return (phi_edge=aphi, dphi=dphi)
+return (phi_edge=phi_edge, dphi=dphi)
 end
 
 function nda_qmc_fixed(phi,qmc_nda_data)
@@ -69,9 +69,11 @@ phi_edge=zeros(Nx+1,)
 dphi=zeros(Nx,)
 dx=1.0/Nx
 qout=qmc_sweep(phi,qmc_data)
-#avg2edge!(phi_edge, dphi, qout.phi_avg,qmc_nda_data)
 avg2edge!(phi_edge, dphi, qout.phi_avg, dx,AVE)
-dhat=(qout.J_avg + (dphi./3.0))./qout.phi_avg
+JAV=AVE*qout.J_edge
+#dhat=(qout.J_avg + (dphi./3.0))./qout.phi_avg
+Javg=.5*(JAV[1:Nx] + JAV[2:Nx+1])
+dhat=(Javg + (dphi./3.0))./qout.phi_avg
 #println(norm(dhat-dhatp,Inf),"  ",norm(dphi-dphip))
 #
 # Something's busted here. My normal code has phi at edges. The
@@ -85,17 +87,21 @@ D1=nda_data.D1
 AV=nda_data.AV
 LT = L2 + (D1*(Dhat * AV))
 LT[1,1]=1.0; LT[2:Nx].=0.0;
-LT[Nx,Nx]=1.0; LT[NX,1:Nx-1].=0.0;
+LT[Nx,Nx]=1.0; LT[Nx,1:Nx-1].=0.0;
 bcl=phi_edge[1];
 bcr=phi_edge[Nx+1];
 #bcl=qout.phi_avg[1];
 #bcr=qout.phi_avg[Nx];
+ravg=AVE*qout.phi_edge;
+#bcl=ravg[1];
+#bcr=ravg[Nx];
 rhs=zeros(Nx+1,)
 rhs[1]=bcl
 rhs[Nx+1]=bcr
+println(bcl,"  ",bcr)
 phiout=LT\rhs
 phiave=.5*(phiout[1:Nx]+phiout[2:Nx+1])
-#phiout=qout.phi_avg
+return phiave
 end