New stuff for validation

Author: ctkelley

Siewert:s=1.dat

@@ -1,2 +1,2 @@
-�_ĺ�?ݦ�'i�?�����?|�<�F�?0��
-q��?�&To65�?�������?T��*�?��d[��?xG�'�$�?r�����?x�D��)�?	���Ƕ�?�2�����?�Q�.���?#܉�;�?�ˌB&+�?��<�A�?6�5T�?l�b�?��`�l�?(��ʸ�?
+x�����?f�E'i�?
'����?p?e	F�?mT4&q��?,�y�65�?u��ą��?����*�?$�>u[��?�,�6�$�?B	
+����?`�%b�)�?��k�Ŷ�?0&!?���?��7����?�Lg;�?�+��%+�?4$��A�?�ʟ�4T�?���*�b�?�i#$�l�?���ʸ�?

Siewert:s=Inf.dat

@@ -28,6 +28,7 @@ export ctk_qmc_test
 export qmc_vs_sn
 export makeqmctab
 export readdata
+export sn_validate
 
 include("validate.jl")
 include("sn_angles.jl")

src/NDA_QMC.jl

@@ -1,2 +1,2 @@
-�_ĺ�?ݦ�'i�?�����?|�<�F�?0��
-q��?�&To65�?�������?T��*�?��d[��?xG�'�$�?r�����?x�D��)�?	���Ƕ�?�2�����?�Q�.���?#܉�;�?�ˌB&+�?��<�A�?6�5T�?l�b�?��`�l�?(��ʸ�?
+x�����?f�E'i�?
'����?p?e	F�?mT4&q��?,�y�65�?u��ą��?����*�?$�>u[��?�,�6�$�?B	
+����?`�%b�)�?��k�Ŷ�?0&!?���?��7����?�Lg;�?�+��%+�?4$��A�?�ʟ�4T�?���*�b�?�i#$�l�?���ʸ�?

src/Siewert:s=1.dat

@@ -1,7 +1,7 @@
 #using LinearAlgebra
 #include("qmc_sweep.jl")
 """
- qmc_source_iteration(s, qmc_data,tol=1.e-8)
+ qmc_source_iteration(s, qmc_data,tol=1.e-4)
  quasi-monte carlo source iteration example script for transport equation.
 
  This is one of the test cases from

src/Siewert:s=Inf.dat

@@ -2,7 +2,7 @@
 #include("ckt_qmc_init.jl")
 #include("ctk_qmc_source_iteration.jl")
 
-function ctk_qmc_test(N=10^4, Nx=50, na2=11; s=1, plotme=false)
+function ctk_qmc_test(N=10^4, Nx=50, na2=11; s=1, plotme=false, tol=1.e-4)
 
 ###############################################################################
 #### Inputs
@@ -19,7 +19,7 @@ function ctk_qmc_test(N=10^4, Nx=50, na2=11; s=1, plotme=false)
 # initialize qmc
 qmc_data = qmc_init(N, Nx, na2, s)
 # call qmc SI
-phi_avg, phi_edge, J_avg, J_edge, psi_right, psi_left, history, itt = qmc_source_iteration(s,qmc_data)
+phi_avg, phi_edge, J_avg, J_edge, psi_right, psi_left, history, itt = qmc_source_iteration(s,qmc_data,tol)
 #
 if plotme 
 plotqmd(phi_avg, phi_edge, J_avg, J_edge, 

src/Tim_QMC/ctk_qmc_source_iteration.jl

@@ -3,9 +3,9 @@ siewert(s=1.0)
 
 Gets the Garcia-Siewert tabular data and parks it where it's easy to find.
 """
-function siewert(s=1.0; filedump=false)
-    na2 = 20
-    nx = 2001
+function siewert(s=1.0; filedump=false, na2=80, nx=16001)
+#    na2 = 20; nx = 2001;
+#    na2 = 80; nx = 16001;
     #
     # precomputed data
     #
@@ -15,7 +15,7 @@ function siewert(s=1.0; filedump=false)
     #
     gout=gmres_iteration(sn_data,s)
     fgm = gout.sol
-    tabout=sn_tabulate(s, nx, fgm; maketab=false)
+    tabout=sn_tabulate(s, nx, fgm; maketab=false, phiedge=true)
     Dout=[tabout.left tabout.right]
     if filedump
     s == 1.0 ? sstr=string(1) : sstr="Inf"

src/Tim_QMC/ctk_qmc_test.jl

@@ -1,5 +1,5 @@
 """
-sn_init(nx, na2, s)
+sn_init(nx, na2, s; siewert=false, phiedge=true)
 
 I pass a named tuple of precomputed and preallocated data to
 all the functions and solvers. 
@@ -15,7 +15,7 @@ s = parameter in Garcia/Siewert
 
 angles, weights = nodes and weights for the angular quadrature
 """
-function sn_init(nx, na2, s; siewert=false)
+function sn_init(nx, na2, s; siewert=false, phiedge=true)
     if siewert
     angles = [-.05; collect(-.1:-.1:-1.0); 0.05; collect(0.1:0.1:1.0)]
     weights = angles
@@ -31,7 +31,8 @@ function sn_init(nx, na2, s; siewert=false)
     psi_left = ones(na)
     source = zeros(nx)
     ptmp = zeros(na)
-    psi = zeros(na2, nx)
+    phiedge ? np=nx : np=nx+1
+    psi = zeros(na2, np)
     return sn_data = (
         c = c,
         dx = dx,

src/siewert.jl

@@ -13,7 +13,7 @@ function sn_tabulate(s, nx, flux; maketab=true, phiedge=true)
     #
     na2 = length(angleout)
     na = floor(Int, na2 / 2)
-    sn_data = sn_init(nx, na2, s; siewert=true)
+    sn_data = sn_init(nx, na2, s; siewert=true, phiedge=phiedge)
     psi = sn_data.psi
     psi = transport_sweep!(psi, flux, sn_data; phiedge=phiedge)
     if maketab

src/sn_init.jl

@@ -17,11 +17,12 @@ function transport_sweep!(psi, phi, sn_data; phiedge=true)
     na2 = length(angles)
     na = floor(Int, na2 / 2)
     nx = length(phi)
+    phiedge ? np=nx : np=nx+1
     #
     # Initialize the angular flux and set the boundary conditions.
     #
     psi *= 0.0
-    @views psi[1:na, nx] = psi_right
+    @views psi[1:na, np] = psi_right
     @views psi[na+1:na2, 1] = psi_left
     #
     #
@@ -41,7 +42,7 @@ function transport_sweep!(psi, phi, sn_data; phiedge=true)
     #
     # Forward sweep
     #
-    @views for ix = 2:nx
+    @views for ix = 2:np
         copy!(psi[na+1:na2, ix], psi[na+1:na2, ix-1])
         psi[na+1:na2, ix] .*= vfr
         psi[na+1:na2, ix] .+= source_average[ix-1]
@@ -50,7 +51,7 @@ function transport_sweep!(psi, phi, sn_data; phiedge=true)
     #
     # Backward sweep
     #
-    @views for ix = nx-1:-1:1
+    @views for ix = np-1:-1:1
         copy!(psi[1:na, ix], psi[1:na, ix+1])
         psi[1:na, ix] .*= vfr
         psi[1:na, ix] .+= source_average[ix]

src/sn_tabulate.jl

@@ -2,11 +2,11 @@
 This makes the table in the notebook showing the errors in QMC as a 
 function of Nx and N
 
-qmc_vs_sn(levels=6, Nrange=4; s=1.0)
+qmc_vs_sn(levels=6, Nrange=3; s=1.0)
 """
-function qmc_vs_sn(levels=6, Nrange=4; s=1.0)
+function qmc_vs_sn(levels=6, Nrange=3; s=1.0)
 NVals=zeros(Int64,Nrange)
-NVals[1]=1000
+NVals[1]=2000
 for i=2:Nrange
 NVals[i]=2*NVals[i-1]
 end
@@ -25,6 +25,23 @@ makeqmctab(Dcomp,NVals,NxVals)
 return Dcomp
 end
 
+function sn_validate(Nrange=3; s=1.0)
+NVals=zeros(Int64,Nrange)
+NVals[1]=100
+for i=2:Nrange
+NVals[i]=2*NVals[i-1]
+end
+Dcomp=zeros(Nrange);
+#
+DataGS=readdata(s)
+for idv=1:Nrange
+Dval=siewert(s; nx=NVals[idv], na2=80)
+Dcomp[idv]=norm(Dval-DataGS,Inf)
+end
+return Dcomp
+end
+
+
 function makeqmctab(Dcomp,NVals,NxVals)
 printf(fmt::String, args...) = @eval @printf($fmt, $(args...))
 levels=length(NxVals)
@@ -49,9 +66,11 @@ function validate(N, levels, DataGS, s=1.0)
 Nx = 50
 Diffs=zeros(levels,)
 for il=0:levels-1
-DataQMC=ctk_qmc_test(N, Nx; s=s, plotme=false)
+DataQMC=ctk_qmc_test(N, Nx; s=s, plotme=false, tol=1.e-4)
 EDiff=(DataGS-DataQMC)./DataGS;
 Diffs[il+1]=norm(EDiff,Inf);
+#EDiff=(DataGS-DataQMC)
+#Diffs[il+1]=norm(EDiff)/norm(DataGS)
 Nx *= 2
 end
 return Diffs