Added uq_inf_lat test problem

test/verification/analytic/process.py

@@ -8,8 +8,13 @@
 
 for name in task.keys():
     os.chdir(name)
-    N_min = task[name]["N_lim"][0]
-    N_max = task[name]["N_lim"][1]
     N = task[name]["N"]
-    os.system("python process.py %i %i %i" % (N_min, N_max, N))
+    if "N_batch" in task[name]:
+        N_h = task[name]["N_lim"]
+        N_b = task[name]["N_batch"]
+        os.system("python process.py %i %i %i %i %i" % (N_h[0], N_h[1], N_b[0], N_b[1], N))
+    else:
+        N_min = task[name]["N_lim"][0]
+        N_max = task[name]["N_lim"][1]
+        os.system("python process.py %i %i %i" % (N_min, N_max, N))
     os.chdir(r"..")

test/verification/analytic/run.py

@@ -20,8 +20,12 @@ def run(N_hist, name, N_batch=1):
     """
     N_hist: number of histories
     name: name for the job
+    N_batch: number of batches
     """
-    output = "output_%i" % N_hist
+    if N_batch == 1:
+        output = "output_%i" % N_hist
+    else:
+        output = "output_%i_%i" % (N_hist, N_batch)
 
     if srun > 1:
         os.system(
@@ -46,14 +50,28 @@ def run(N_hist, name, N_batch=1):
 
 for name in task.keys():
     os.chdir(name)
+    N = task[name]["N"]
     if "N_batch" in task[name]:
-            x = 1
+        N_hist = int(10**int(task[name]["N_lim"][0]))
+        N_min = task[name]["N_batch"][0]
+        N_max = task[name]["N_batch"][1]
+        for N_batch in np.logspace(N_min, N_max, N):
+            N_batch = int(N_batch)
+            print("N_batch %i" % N_batch)
+            run(N_hist, name, N_batch)
+        N_batch = int(10**int(task[name]["N_batch"][0]))
+        N_min = task[name]["N_lim"][0]
+        N_max = task[name]["N_lim"][1]
+        hist_list = np.logspace(N_min, N_max, N)
+        for i in range(1, N):
+            N_hist = int(hist_list[i])
+            print("N_hist %i" % N_hist)
+            run(N_hist, name, N_batch)
     else:
         N_min = task[name]["N_lim"][0]
         N_max = task[name]["N_lim"][1]
-        N = task[name]["N"]
         for N_hist in np.logspace(N_min, N_max, N):
             N_hist = int(N_hist)
             print(name, N_hist)
             run(N_hist, name)
-        os.chdir(r"..")
+    os.chdir(r"..")

test/verification/analytic/task.py

@@ -20,7 +20,13 @@
 task["slab_absorbium"]["N_lim"] = [5, 9]
 task["slab_absorbium"]["N"] = 9
 
-# Slab isotropic beam
-task["slab_isobeam_td"] = {}
-task["slab_isobeam_td"]["N_lim"] = [3, 7]
-task["slab_isobeam_td"]["N"] = 9
+ # Slab isotropic beam
+ task["slab_isobeam_td"] = {}
+ task["slab_isobeam_td"]["N_lim"] = [3, 7]
+ task["slab_isobeam_td"]["N"] = 9
+
+# UQ - infinite 7-group
+task["uq_inf_lat"] = {}
+task["uq_inf_lat"]["N_lim"] = [3, 6]
+task["uq_inf_lat"]["N_batch"] = [1, 3]
+task["uq_inf_lat"]["N"] = 3

test/verification/analytic/uq_inf_lat/input.py

@@ -0,0 +1,65 @@
+import numpy as np
+import h5py
+
+import mcdc
+
+
+# =========================================================================
+# Set model and run
+# =========================================================================
+
+lib = h5py.File("c5g7.h5", "r")
+
+def set_mat(mat):
+    return mcdc.material(
+        capture=mat["capture"][:],
+        scatter=mat["scatter"][:],
+        fission=mat["fission"][:],
+        nu_p=mat["nu_p"][:],
+        nu_d=mat["nu_d"][:],
+        chi_p=mat["chi_p"][:],
+        chi_d=mat["chi_d"][:],
+        speed=mat["speed"],
+        decay=mat["decay"],
+    )
+
+
+def uq_mat(mat, param, c):
+    return c*mat[param][:]
+
+
+mat_uo2 = set_mat(lib["uo2"])  # Fuel: UO2
+mat_mod = set_mat(lib["mod"])  # Moderator
+mat_cr = set_mat(lib["cr"])  # Control rod
+
+s1 = mcdc.surface("plane-x", x=0.0, bc="reflective")
+s2 = mcdc.surface("plane-x", x=0.5)
+s3 = mcdc.surface("plane-x", x=1.5)
+s4 = mcdc.surface("plane-x", x=2.0, bc="reflective")
+
+mcdc.cell([+s1, -s2], mat_uo2)
+mcdc.cell([+s2, -s3], mat_mod)
+mcdc.cell([+s3, -s4], mat_cr)
+
+mcdc.source(point=[1.0, 0.0, 0.0], energy=[1, 0, 0, 0, 0, 0, 0], isotropic=True)
+
+scores = ["flux"]
+grid = 200
+mcdc.tally(scores=scores, x=np.linspace(0.0, 2.0, grid+1), g=[-.5, 3.5, 6.5])
+
+mcdc.setting(N_particle=1E2, N_batch=1E2)
+
+mcdc.uq(material=mat_mod, distribution='uniform',
+         capture=uq_mat(lib["mod"], 'capture', 0.8),
+         scatter=uq_mat(lib["mod"], 'scatter', 0.8),
+         )
+
+mcdc.uq(material=mat_cr, distribution='uniform',
+        capture=uq_mat(lib["cr"], 'capture', 0.8),
+        scatter=uq_mat(lib["cr"], 'scatter', 0.8),
+        )
+
+mcdc.run()
+
+
+

test/verification/analytic/uq_inf_lat/process.py

@@ -0,0 +1,45 @@
+import numpy as np
+import h5py
+import sys
+import ast
+
+sys.path.append("../")
+import tool
+
+# Reference solution
+with h5py.File("reference.h5","r") as f:
+    x = f["tally/grid/x"][:]
+    dx = (x[1:] - x[:-1])[0]
+    xmid = 0.5 * (x[:-1] + x[1:])
+    var_ref = f["tally/flux/varp"][:]
+
+for i in range(1,6):
+    sys.argv[i] = int(sys.argv[i])
+
+hist_lim = [sys.argv[1], sys.argv[2]]
+batch_lim = [sys.argv[3], sys.argv[4]]
+N = int(sys.argv[5])
+
+# Batch cases run
+N_hist = 10**int(hist_lim[0])
+N_batch_list = np.logspace(batch_lim[0], batch_lim[1], N)
+error = np.zeros(N)
+error_max = np.zeros(N)
+for i, N_batch in enumerate(N_batch_list):
+    # Get results
+    with h5py.File("output_%i_%i.h5" % (int(N_hist), int(N_batch)),"r") as f:
+        varp = f["tally/flux/uq_var"][:]
+    error[i] = tool.error(varp, var_ref)
+    error_max[i] = tool.error_max(varp, var_ref)
+tool.plot_convergence("uq_inf_lat_batches", N_batch_list, error, error_max)
+
+# History cases run
+N_batch = 10**int(batch_lim[0])
+N_hist_list = np.logspace(hist_lim[0], hist_lim[1], N)
+for i, N_hist in enumerate(N_hist_list):
+    # Get results
+    with h5py.File("output_%i_%i.h5" % (int(N_hist), int(N_batch)),"r") as f:
+        varp = f["tally/flux/uq_var"][:]
+    error[i] = tool.error(varp, var_ref)
+    error_max[i] = tool.error_max(varp, var_ref)
+tool.plot_convergence("uq_inf_lat_hist", N_hist_list, error, error_max)

test/verification/analytic/uq_inf_lat/reference.py

@@ -0,0 +1,47 @@
+import numpy as np
+import h5py
+
+import h5py
+import numpy as np
+import matplotlib.pyplot as plt
+
+n_runs = 100
+p = np.array([[2, 0], [5, 0], [1, 1], [5, 1], [1, 2], [1, 3], [2, 3]])
+b = np.array([[5, 3], [2, 3], [1, 3], [2, 2], [1, 2], [1, 1], [5, 0]])
+
+assert len(p) == len(b), "Number of runs not equal."
+with h5py.File('p16b11/combined_results.h5','r') as g:
+        x = g['tally/grid/x'][:]
+        dx = (x[1:] - x[:-1])[0]
+        xmid = 0.5 * (x[:-1] + x[1:])
+        phi_ref = g['tally/flux/mean'][:]
+        var_ref = g['tally/flux/varp'][:] / (dx**2)
+
+n_eta = np.zeros(len(p))
+n_xi = n_eta.copy()
+run_tot = n_eta.copy()
+run_prep = n_eta.copy()
+run_out = n_eta.copy()
+varp_norm = n_eta.copy()
+vart_norm = varp_norm.copy()
+for i in range(len(p)):
+    n_eta[i] = p[i][0] * (10**p[i][1])
+    n_xi[i] = n_runs * b[i][0] * (10**b[i][1])
+    outdir = 'p' + str(p[i][0]) + str(p[i][1]) + 'b' + str(b[i][0]) + str(b[i][1])
+    with h5py.File(outdir+'/combined_results.h5','r') as f:
+        phi = f['tally/flux/mean'][:]
+        varp = f['tally/flux/varp'][:] / (dx**2)
+        vart = f['tally/flux/vart'][:] / (dx**2)
+        run_tot[i] = f['runtime/total'][:]
+        run_prep[i] = f['runtime/preparation'][:]
+        run_out[i] = f['runtime/output'][:]
+
+        varp_err = (varp - var_ref)/var_ref
+        varp_norm[i] = np.sqrt(np.sum(np.square(varp_err)))
+        vart_err = (vart - var_ref)/var_ref
+        vart_norm[i] = np.sqrt(np.sum(np.square(vart_err)))
+
+cost_tot = n_eta * n_xi
+run_time = run_tot - (run_prep+run_out)/100*99
+varp_fom = 1/(varp_norm**2)/run_time
+vart_fom = 1/(vart_norm**2)/run_time