Cleanup examples folder

examples/eigenvalue/2d_c5g7/input.py

@@ -8,7 +8,7 @@
 # =============================================================================
 
 # Load material data
-lib = h5py.File("c5g7.h5", "r")
+lib = h5py.File("2d_c5g7_xs.h5", "r")
 
 
 # Materials
@@ -17,22 +17,17 @@ def set_mat(mat):
         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"],
+        nu_p=mat["nu"][:],
+        chi_p=mat["chi"][:],
     )
 
 
 mat_uo2 = set_mat(lib["uo2"])
 mat_mox43 = set_mat(lib["mox43"])
-mat_mox7 = set_mat(lib["mox7"])
+mat_mox7 = set_mat(lib["mox70"])
 mat_mox87 = set_mat(lib["mox87"])
 mat_gt = set_mat(lib["gt"])
 mat_fc = set_mat(lib["fc"])
-mat_cr = set_mat(lib["cr"])
 mat_mod = set_mat(lib["mod"])
 
 # =============================================================================
@@ -52,7 +47,6 @@ def set_mat(mat):
 mox8 = mcdc.cell([-cy], mat_mox87)
 gt = mcdc.cell([-cy], mat_gt)
 fc = mcdc.cell([-cy], mat_fc)
-cr = mcdc.cell([-cy], mat_cr)
 mod = mcdc.cell([+cy], mat_mod)
 modi = mcdc.cell([-cy], mat_mod)  # For all-water lattice
 
@@ -63,7 +57,6 @@ def set_mat(mat):
 n = mcdc.universe([mox8, mod])["ID"]
 g = mcdc.universe([gt, mod])["ID"]
 f = mcdc.universe([fc, mod])["ID"]
-c = mcdc.universe([cr, mod])["ID"]
 w = mcdc.universe([modi, mod])["ID"]
 
 # =============================================================================

examples/eigenvalue/2d_c5g7_iqmc/input.py

@@ -164,7 +164,7 @@ def set_mat(mat):
 # =============================================================================
 # iQMC Parameters
 # =============================================================================
-N = 1e3
+N = 1e4
 maxit = 10
 tol = 1e-3
 Nx = 17 * 3 * 2

examples/eigenvalue/2d_c5g7_iqmc/process.py

@@ -8,10 +8,10 @@
 # =============================================================================
 
 # Load iqmc result
-with h5py.File("PI_output.h5", "r") as f:
+with h5py.File("output.h5", "r") as f:
     x = f["iqmc/grid/x"][:]
     y = f["iqmc/grid/y"][:]
-    phi_avg = f["iqmc/flux"][:]
+    phi_avg = f["iqmc/tally/flux"][:]
     f.close()
 
 

examples/eigenvalue/slab_2gu_iqmc/input.py

@@ -9,17 +9,6 @@
 # "Analytical Benchmark Test Set For Criticality Code Verification"
 
 # Set materials
-
-# 1G-PU Slab data
-# m1 = mcdc.material(
-#         capture=np.array([0.019584]),
-#         scatter=np.array([[0.225216]]),
-#         fission=np.array([0.081600]),
-#         nu_p=np.array([2.84]),
-#         chi_p=np.array([1.0]),
-#     )
-# R = [0.0, 4.513502]
-
 # 2G-U Slab data
 m1 = mcdc.material(
     capture=np.array([0.01344, 0.00384]),

examples/eigenvalue/slab_kornreich_iqmc/input.py

@@ -35,14 +35,13 @@
 # =============================================================================
 # iQMC Parameters
 # =============================================================================
-N = 1000
+N = 5000
 maxit = 25
 tol = 1e-3
 x = np.arange(0.0, 2.6, 0.1)
 Nx = len(x) - 1
 generator = "halton"
 solver = "power_iteration"
-inner_solver = "source_iteration"
 fixed_source = np.zeros(Nx)
 phi0 = np.ones((Nx))
 
@@ -58,8 +57,7 @@
     tol=tol,
     generator=generator,
     eigenmode_solver=solver,
-    fixed_source_solver=inner_solver,
-    # score=['tilt-x']
+    score=["tilt-x"],
 )
 # Setting
 mcdc.setting(N_particle=N)

examples/eigenvalue/slab_kornreich_iqmc/process.py

@@ -3,6 +3,26 @@
 import h5py
 import sys
 
+# =============================================================================
+# Import data
+# =============================================================================
+
+with h5py.File("output.h5", "r") as f:
+    # Note the spatial (dx) and source strength (100+1) normalization
+    keff = f["k_eff"][()]
+    phi_avg = f["iqmc/tally/flux"][:]
+    sweeps = f["iqmc/sweep_count"][()]
+    x = f["iqmc/grid/x"][:]
+    dx = x[1] - x[0]
+    x_mid = 0.5 * (x[:-1] + x[1:])
+    f.close()
+
+tmp = 0.5 * (phi_avg[1:] + phi_avg[:-1])
+norm = np.sum(tmp * dx)
+phi_avg /= norm
+print("Keff = ", keff)
+print("Number of QMC Transport Sweeps = ", sweeps)
+
 # =============================================================================
 # Reference solution
 # =============================================================================
@@ -69,27 +89,7 @@
 plt.figure(dpi=300, figsize=(8, 5))
 plt.plot(x_exact, phi_exact, label="sol")
 
-with h5py.File("output.h5", "r") as f:
-    # Note the spatial (dx) and source strength (100+1) normalization
-    keff = f["k_eff"][()]
-    phi_avg = f["iqmc/flux"][:]
-    sweeps = f["iqmc/sweep_count"][()]
-    x = f["iqmc/grid/x"][:]
-    dx = x[1] - x[0]
-    x_mid = 0.5 * (x[:-1] + x[1:])
-    f.close()
-
-tmp = 0.5 * (phi_avg[1:] + phi_avg[:-1])
-norm = np.sum(tmp * dx)
-phi_avg /= norm
-print("Keff = ", keff)
-print(sweeps)
-plt.plot(x_mid, phi_avg, label="PI")
-
-
-# =============================================================================
-# Finish Plot
-# =============================================================================
+plt.plot(x_mid, phi_avg)
 plt.title("Kornreich et al. Slab")
 plt.ylabel(r"$\phi(x)$")
 plt.xlabel(r"$x$")

examples/fixed_source/cooper2_iqmc/input.py

@@ -35,13 +35,12 @@
 # =============================================================================
 N = 1e4
 Nx = Ny = 40
-maxit = 15
-tol = 1e-9
+maxit = 30
+tol = 1e-6
 x = np.linspace(0, 4, num=Nx + 1)
 y = np.linspace(0, 4, num=Ny + 1)
 generator = "halton"
 solver = "gmres"
-tilt = 1
 
 # fixed source in lower left corner
 fixed_source = np.zeros((Nx, Ny))
@@ -57,8 +56,8 @@
     maxitt=maxit,
     tol=tol,
     generator=generator,
-    source_tilt=tilt,
     fixed_source_solver=solver,
+    score=["tilt-x", "tilt-y"],
 )
 
 # =============================================================================

examples/fixed_source/cooper2_iqmc/process.py

@@ -1,36 +1,14 @@
 import matplotlib.pyplot as plt
-
-plt.rcParams.update({"font.size": 14})
-
-import matplotlib as mpl
 import h5py
 import numpy as np
 
-# from mpl_toolkits.mplot3d import axes3d, Axes3D
-
-
-# Load iqmc result
-# with h5py.File("output.h5", "r") as f:
-#     x = f["iqmc/grid/x"][:]
-#     dx = [x[1:] - x[:-1]][-1]
-#     x_mid = 0.5 * (x[:-1] + x[1:])
-
-#     phi = f["iqmc/flux"][:]
-#     f.close()
-
 # Load iqmc result
-with h5py.File("output2.h5", "r") as f:
+with h5py.File("output.h5", "r") as f:
     meshx = f["iqmc/grid/x"][:]
     meshy = f["iqmc/grid/y"][:]
     dx = [meshx[1:] - meshx[:-1]][-1]
     x_mid = 0.5 * (meshx[:-1] + meshx[1:])
-    lowX = meshx[:-1]
-    highX = meshx[1:]
-
-    phi = f["iqmc/flux"][:]
-    q = f["iqmc/source"][:][0, 0, :, :, 0]
-    qdotx = f["iqmc/Q11/x"][:][0, 0, :, :, 0]
-    qdoty = f["iqmc/Q11/x"][:][0, 0, :, :, 0]
+    phi = f["iqmc/tally/flux"][:]
 
     f.close()
 
@@ -46,115 +24,4 @@
 ax.set_zlabel(r"log($\phi$)", rotation=180)
 
 ax.view_init(elev=15, azim=20)
-# iqmc_label = mpl.lines.Line2D([0], [0], linestyle="none", c="b", marker="o")
-# mc_label = mpl.lines.Line2D([0], [0], linestyle="none", c="r", marker="o")
-# ax.legend([iqmc_label, mc_label], ["iQMC", "MC/DC"], numpoints=1)
-# plt.tight_layout()
-plt.show()
-
-# =============================================================================
-# Source Plots
-# =============================================================================
-X, Y = np.meshgrid(x_mid, x_mid)
-Z = np.log10(q)
-fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, dpi=300, figsize=(12, 10))
-ax.plot_surface(Y, X, Z, edgecolor="b", linewidth=0.5, cmap="viridis")
-ax.set_xlabel("x")
-ax.set_ylabel("y")
-ax.set_zlabel(r"Source", rotation=180)
-ax.view_init(elev=15, azim=20)
-
-
-num_steps = 800
-y_mid = x_mid.copy()
-Nx = num_steps
-Ny = num_steps
-x = np.linspace(lowX[0], highX[-1], num=num_steps)
-y = np.linspace(lowX[0], highX[-1], num=num_steps)
-# Nx = meshx.size
-# Ny = meshy.size
-# x = meshx
-# y = meshy
-z1 = np.zeros(shape=(Nx, Ny))
-z2 = np.zeros(shape=(Nx, Ny))
-z3 = np.zeros(shape=(Nx, Ny))
-z4 = np.zeros(shape=(Nx, Ny))
-points = np.zeros(shape=(int(Nx * Ny), 2))
-count = 0
-for i in range(Nx):
-    zonex = np.argmax((x[i] > lowX) * (x[i] <= highX))
-    midx = x_mid[zonex]
-    for j in range(Ny):
-        zoney = np.argmax((y[j] > lowX) * (y[j] <= highX))
-        midy = y_mid[zoney]
-        points[count, :] = [x[i], y[j]]
-        count += 1
-        z1[i, j] = q[zonex, zoney]
-        z2[i, j] = q[zonex, zoney] + qdotx[zonex, zoney] * (x[i] - midx)
-        z3[i, j] = q[zonex, zoney] + qdoty[zonex, zoney] * (y[j] - midy)
-        z4[i, j] = (
-            q[zonex, zoney]
-            + qdotx[zonex, zoney] * (x[i] - midx)
-            + qdoty[zonex, zoney] * (y[j] - midy)
-        )
-
-X, Y = np.meshgrid(x, y)
-Z1 = np.log10(z1)
-fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, dpi=300, figsize=(12, 10))
-ax.plot_surface(Y, X, z1, edgecolor="b", linewidth=0.5, cmap="viridis")
-ax.set_xlabel("x")
-ax.set_ylabel("y")
-ax.set_zlabel(r"Source", rotation=180)
-ax.view_init(elev=15, azim=20)
-plt.show()
-
-Z2 = np.log10(z2)
-fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, dpi=300, figsize=(12, 10))
-ax.plot_surface(Y, X, z2, edgecolor="b", linewidth=0.5, cmap="viridis")
-ax.set_xlabel("x")
-ax.set_ylabel("y")
-ax.set_zlabel(r"Source - Tilted X", rotation=180)
-ax.view_init(elev=15, azim=20)
-plt.show()
-
-Z3 = np.log10(z3)
-fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, dpi=300, figsize=(12, 10))
-ax.plot_surface(Y, X, z3, edgecolor="b", linewidth=0.5, cmap="viridis")
-ax.set_xlabel("x")
-ax.set_ylabel("y")
-ax.set_zlabel(r"Source - Tilted Y", rotation=180)
-ax.view_init(elev=15, azim=20)
-plt.show()
-
-Z4 = np.log10(z4)
-fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, dpi=300, figsize=(12, 10))
-ax.plot_surface(Y, X, z4, edgecolor="b", linewidth=0.5, cmap="viridis")
-ax.set_xlabel("x")
-ax.set_ylabel("y")
-ax.set_zlabel(r"Source - Tilted X,Y", rotation=180)
-ax.view_init(elev=15, azim=20)
-plt.show()
-
-# =============================================================================
-# Scipy Interpolate plots
-# =============================================================================
-
-from scipy.interpolate import RegularGridInterpolator
-
-X, Y = np.meshgrid(x_mid, x_mid)
-interp = RegularGridInterpolator(
-    (x_mid, x_mid), q, method="linear", bounds_error=False, fill_value=None
-)
-result = interp(points)
-result = np.reshape(result, (Nx, Ny))
-X, Y = np.meshgrid(x, y)
-fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, dpi=300, figsize=(12, 10))
-ax.plot_surface(Y, X, np.log10(result), edgecolor="b", linewidth=0.5, cmap="viridis")
-ax.set_xlabel("x")
-ax.set_ylabel("y")
-ax.set_zlabel(r"Source")
-ax.view_init(elev=15, azim=20)
-# plt.xlim((0,10))
-# plt.ylim((0,10))
-plt.title("Source - Scipy Interp. X, Y, XY")
 plt.show()