Moving old stuff out of the way.

Author: ctkelley

NDA_Notebook.ipynb

@@ -11,6 +11,7 @@ using PyPlot
 using Printf
 using Sobol
 
+export fprintTeX
 export testgauss
 export sn_angles
 export source_iteration
@@ -26,6 +27,7 @@ export writetab
 export readtab
 export qmc_vs_sn
 export makeqmctab
+export makeqmctex
 export readdata
 export sn_validate
 export qmc_sweep
@@ -59,4 +61,5 @@ 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")
 end

src/Data_Orig/Siewert:s=1.dat ctk_nda_src/Data_Orig/Siewert:s=1.dat

@@ -0,0 +1,76 @@
+#using LinearAlgebra
+#include("qmc_sweep.jl")
+"""
+ 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
+
+ Radiative transfer in finite inhomogeneous plane-parallel atmospheres
+ by Garcia and Siewert
+ JQSRT (27), 1982 pp 141-148.
+
+ Inputs:
+ -------
+     N:
+         Number of internal source particles
+     Nx:
+         Number of scalar flux and current tallying cells
+     na2:
+         Number of exit angular flux tally bins
+     s:
+         Scattering cross section varies with exp(-x/s)
+ Outputs:
+ --------
+     phi_avg:
+         Array of length Nx, average scalar flux across each cell
+     phi_edge:
+         Array of length Nx+1, scalar flux at cell boundaries
+     J_avg:
+         Array of length Nx, average cell current
+     J_edge:
+         Array of length Nx+1, current at cell boundaries
+     psi_right:
+         Array of length nbins, angular flux distribution at right boundary
+     psi_left:
+         Array of length nbins, angular flux distribution at left boundary
+     history:
+         Infinity norm of change in phi_avg with each iteration
+     itt:
+         Number of iterations
+"""
+
+
+function qmc_source_iteration(s, qmc_data, tol=1.e-4)
+    #
+    # precomputed data
+    #
+    N = qmc_data.N
+    Nx = qmc_data.Nx
+    itt = 0
+    delflux = 1
+    phi_avg = zeros(Nx)
+    phi_avg_old = zeros(Nx)
+    reshist = []
+    #globalize variables
+    phi_edge = J_avg = J_edge = exit_right_bins = exit_left_bins = 0
+
+    while itt < 200 && delflux > tol
+        phi_avg, phi_edge, J_avg, J_edge, exit_right_bins, exit_left_bins = qmc_sweep(phi_avg,qmc_data)
+        delflux = norm(phi_avg - phi_avg_old, Inf)
+        itt += 1
+        push!(reshist, delflux)
+        phi_avg_old .= phi_avg
+        #println("**********************")
+        #println("Iteration: ", itt," change = ",delflux)
+    end
+
+    return (phi_avg = phi_avg,
+            phi_edge = phi_edge,
+            J_avg = J_avg,
+            J_edge = J_edge,
+            psi_right = exit_right_bins,
+            psi_left = exit_left_bins,
+            history = reshist,
+            itt = itt)
+end

src/Data_Orig/Siewert:s=Inf.dat ctk_nda_src/Data_Orig/Siewert:s=Inf.dat

@@ -0,0 +1,92 @@
+#using PyPlot
+#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, tol=1.e-4)
+
+###############################################################################
+#### Inputs
+###############################################################################
+#N = 10^4 #number of QMC particles
+#Nx = 50 #number of tally cells
+#na2 = 11 #number of angles for angular mesh
+#s = 1 #parameter in Garcia/Siewert
+
+#phi_edge = qmc_data.phi_edge
+###############################################################################
+#### Function Call
+###############################################################################
+# 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,tol)
+#
+if plotme 
+plotqmd(phi_avg, phi_edge, J_avg, J_edge, 
+               psi_right, psi_left, history, itt, qmc_data)
+end
+qmctab=sn_tabulate(s, Nx, phi_avg; maketab=false, phiedge=false)
+return [qmctab.left qmctab.right]
+#
+end
+
+function plotqmd(phi_avg, phi_edge, J_avg, J_edge, 
+               psi_right, psi_left, history, itt, qmc_data)
+###############################################################################
+#### Plotting
+###############################################################################
+####
+iteration = 1:itt
+
+midpoints = qmc_data.midpoints
+edges = qmc_data.edges
+
+####
+figure(figsize = (6,12))
+
+subplot(311)
+plot(midpoints,phi_avg)
+ylabel("phi midpoints")
+xlabel("midpoints")
+title("Cell Averaged Scalar Flux")
+
+# cell Averaged Spatial Derivative of Scalar Flux
+subplot(312)
+plot(edges,phi_edge)
+ylabel("phi edges")
+xlabel("cell edges")
+title("Cell edge Scalar Flux")
+
+# cell averaged current
+subplot(313)
+plot(midpoints,J_avg)
+ylabel("current")
+xlabel("midpoints")
+title("Cell Averaged Current")
+
+#left exit bins
+figure(figsize = (5,10))
+subplot(211)
+gsSol = zeros((11,2))
+gsSol[:,1] = -1*[0.05,.1,.2,.3,.4,.5,.6,.7,.8,.9,1]
+gsSol[:,2] = [0.58966,0.53112,0.44328,0.38031,0.33296,0.29609,0.26656,0.24239,
+                                                        0.22223,0.20517,0.19055]
+plot(gsSol[:,1], gsSol[:,2],"+",label="Garcia et al.")
+plot(psi_left[:,1], psi_left[:,2],label="QMC_sweep")
+xlabel("-mu")
+ylabel("Left Boundary - Angular Flux Dist.")
+title("Angular Flux Exit Distributions from Garcia et al., c = 1")
+legend()
+
+#right exit bins
+subplot(212)
+gsSol = zeros((11,2))
+gsSol[:,1] = 1*[0.05,.1,.2,.3,.4,.5,.6,.7,.8,.9,1]
+gsSol[:,2] = [6.08E-06,6.93E-06,9.64E-06,1.62E-05,4.39E-05,1.69E-04,5.73E-04,
+                                            1.51E-03,3.24E-03,5.96E-03,9.77E-03]
+plot(gsSol[:,1], gsSol[:,2],"+", label="Garcia et al.")
+plot(psi_right[:,1], psi_right[:,2], label="QMC_sweep")
+xlabel("mu")
+ylabel("Right Boundary - Angular Flux Dist.")
+legend()
+end

src/NDA_QMC.jl ctk_nda_src/NDA_QMC.jl

@@ -0,0 +1,16 @@
+module NDA_QMC
+using LinearAlgebra
+using LinearAlgebra.BLAS
+using SIAMFANLEquations
+using FastGaussQuadrature
+
+export testgauss
+export sn_angles
+export source_iteration
+
+include("testgauss.jl")
+include("sn_angles.jl")
+include("source_iteration.jl")
+include("transport_sweep.jl")
+end
+

ctk_nda_src/Old_Stuff/Old_Tim_QMC/ctk_qmc_source_iteration.jl

@@ -0,0 +1,18 @@
+"""
+sn_angles(na=20)
+
+Get double Gauss nodes and weights for SN
+"""
+function sn_angles(na=20)
+na2=floor(Int,na/2)
+2*na2 == na || error("odd number of angles")
+baseangles, baseweights = gausslegendre(na2)
+posweights=baseweights * .5
+negweights=copy(posweights)
+posangles=(baseangles .+ 1.0) *.5
+negangles=-copy(posangles)
+weights=[negweights; posweights]
+angles=[negangles; posangles]
+angles, weights
+end
+

ctk_nda_src/Old_Stuff/Old_Tim_QMC/ctk_qmc_test.jl

@@ -0,0 +1,59 @@
+"""
+ source_iteration(na2=20, s=Inf)
+ Source iteration example script for transport equation.
+
+ This is one of the test cases from
+
+ Radiative transfer in finite inhomogeneous plane-parallel atmospheres
+ by Garcia and Siewert
+ JQSRT (27), 1982 pp 141-148.
+
+"""
+function source_iteration(na2=20, s=Inf)
+nx=2001
+(angles, weights) = sn_angles(na2)
+tau=5.0; dx=tau/(nx-1);
+na=floor(Int,na2/2)
+#
+# Problem data
+#
+x=collect(0:dx:tau);
+c=exp.(-x/s);
+psi_right=zeros(na,);
+psi_left=ones(na,);
+source=zeros(nx,);
+ptmp=zeros(na,)
+sweep_data=(c=c, dx=dx, ptmp=ptmp)
+psi=zeros(nx,na2)
+#
+# Source iteration
+#
+itt=0;
+delflux=1;
+phi=zeros(nx,)
+flux=zeros(nx,)
+ithist=[]
+while itt < 200 && delflux > 1.e-8
+   psi = transport_sweep!(psi, phi, psi_left, psi_right, 
+              angles, source, sweep_data);
+   flux.=psi*weights;
+   delflux=norm(flux-phi,Inf); itt=itt+1;
+   push!(ithist,delflux)
+   phi.=flux;
+end
+#
+# Tabulate the exit distributions to check results.
+#
+#angleout=[collect(-1.0:.1:-.1); -.05; .05; collect(.1:.1:1.0)]
+angleout=[-.05; collect(-.1:-.1:-1.0); .05; collect(.1:.1:1.0)]
+na2=length(angleout)
+na=floor(Int,na2/2)
+psi_left=ones(na,);
+psi_right=zeros(na,);
+psi = zeros(nx,na2)
+ptmp=zeros(na,)
+sweep_data=(c=c, dx=dx, ptmp=ptmp)
+psi = transport_sweep!(psi, phi, psi_left, psi_right, 
+                            angleout, source, sweep_data);
+return (left=psi[1,1:na], right=psi[nx,na+1:na2], history=ithist)
+end

ctk_nda_src/Old_Stuff/Transpose_Psi/Jan20/NDA_QMC.jl

@@ -0,0 +1,14 @@
+function testgauss(m)
+x, w = gausslegendre(m)
+# Map to [0,1]
+w .*= .5; x .+= 1.0; x .*= .5;
+errs=zeros(2*m,)
+for p=0:2*m-1
+    exact=1.0/(1.0+p)
+    y=x.^p
+#    approx=y'*w
+    approx=dot(w,y)
+    errs[p+1]=abs(exact-approx)
+end
+maximum(errs)
+end

ctk_nda_src/Old_Stuff/Transpose_Psi/Jan20/sn_angles.jl

@@ -0,0 +1,59 @@
+"""
+transport_sweep!(psi, phi, psi_left, psi_right, angles, source, 
+              sweep_data)
+
+Take a single transport sweep.
+"""
+function transport_sweep!(psi, phi, psi_left, psi_right, angles, source, 
+              sweep_data)
+
+c=sweep_data.c;
+q=source;
+dx=sweep_data.dx;
+ptmp=sweep_data.ptmp;
+na2=length(angles);
+na=floor(Int,na2/2);
+nx=length(phi);
+#
+# Initialize the intensities.
+#
+psi *= 0.0
+@views psi[nx,1:na]=psi_right;
+@views psi[1,na+1:na2]=psi_left;
+#
+#
+#
+source_total = .5*c.*phi + q;
+@views source_average=(source_total[2:nx]+source_total[1:nx-1])*.5;
+@views forward_angles=angles[na+1:na2];
+@views backward_angles=angles[1:na];
+#
+vfl=(forward_angles/dx) .+ .5;
+vfl=1.0 ./ vfl;
+vfr=(forward_angles/dx) .- .5;
+dfr=Diagonal(vfr);
+#
+# Forward sweep
+#
+ptmp.=psi[1,na+1:na2]
+@views for ix=2:nx
+    psi[ix,na+1:na2] .= ptmp.*vfr
+#    psi[ix,na+1:na2] .= psi[ix-1,na+1:na2]
+#    psi[ix,na+1:na2] .*= vfr
+    psi[ix,na+1:na2] .+= source_average[ix-1]
+    psi[ix,na+1:na2] .*= vfl
+    ptmp.=psi[ix,na+1:na2]
+end
+#
+#
+# Backward sweep
+#
+ptmp.=psi[nx,1:na]
+@views for ix=nx-1:-1:1
+    psi[ix,1:na] .= ptmp.*vfr
+    psi[ix,1:na] .+= source_average[ix]
+    psi[ix,1:na] .*= vfl
+    ptmp.=psi[ix,1:na]
+end
+return psi
+end