Merge pull request CEMeNT-PSAAP#140 from spasmann/main

iqmc source tilting and refactorization

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

@@ -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"]
 
 # =============================================================================
@@ -171,29 +164,28 @@ def set_mat(mat):
 # =============================================================================
 # iQMC Parameters
 # =============================================================================
-N = 1e5
+N = 1e4
 maxit = 10
 tol = 1e-3
-pre_sweeps = 9
-x_grid = np.linspace(0.0, pitch * 17 * 3, 17 * 3 * 2 + 1)
-y_grid = np.linspace(-pitch * 17 * 3, 0.0, 17 * 3 * 2 + 1)
+Nx = 17 * 3 * 2
+Ny = 17 * 3 * 2
+G = 7
+x_grid = np.linspace(0.0, pitch * 17 * 3, Nx + 1)
+y_grid = np.linspace(-pitch * 17 * 3, 0.0, Ny + 1)
 
 generator = "halton"
-solver = "davidson"
-
-phi0 = np.zeros((x_grid.size - 1, y_grid.size - 1))
-np.random.seed(123456)
-phi0[: int(pitch * 17 * 2), int(-pitch * 17 * 2) :] = np.random.random((42, 42))
+solver = "power_iteration"
 
+phi0 = np.ones((G, Nx, Ny))
 fixed_source = np.zeros_like(phi0)
 
 mcdc.iQMC(
     x=x_grid,
     y=y_grid,
+    g=np.ones(G),
     phi0=phi0,
     fixed_source=fixed_source,
     eigenmode_solver=solver,
-    preconditioner_sweeps=pre_sweeps,
     maxitt=maxit,
     tol=tol,
     generator=generator,
@@ -203,9 +195,8 @@ def set_mat(mat):
 # run mcdc
 # =============================================================================
 
-
 # Setting
-mcdc.setting(N_particle=N, output_name="davidson_output")
+mcdc.setting(N_particle=N)
 mcdc.eigenmode()
 
 # Run

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]),
@@ -48,9 +37,7 @@
 generator = "halton"
 solver = "davidson"
 fixed_source = np.zeros((2, Nx))
-np.random.seed(123456)
-phi0 = np.random.random((2, Nx))
-# phi0 = np.ones((2, Nx))
+phi0 = np.ones((2, Nx))
 
 # =============================================================================
 # Set tally, setting, and run mcdc

examples/eigenvalue/slab_kornreich_iqmc/input.py

@@ -35,13 +35,13 @@
 # =============================================================================
 # iQMC Parameters
 # =============================================================================
-N = 1000
-maxit = 10
+N = 5000
+maxit = 25
 tol = 1e-3
 x = np.arange(0.0, 2.6, 0.1)
 Nx = len(x) - 1
 generator = "halton"
-solver = "davidson"
+solver = "power_iteration"
 fixed_source = np.zeros(Nx)
 phi0 = np.ones((Nx))
 
@@ -57,6 +57,7 @@
     tol=tol,
     generator=generator,
     eigenmode_solver=solver,
+    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,72 +89,7 @@
 plt.figure(dpi=300, figsize=(8, 5))
 plt.plot(x_exact, phi_exact, label="sol")
 
-# =============================================================================
-# MCDC results
-# =============================================================================
-# Results
-with h5py.File("mc_output.h5", "r") as f:
-    x = f["tally/grid/x"][:]
-    dx = x[1:] - x[:-1]
-    phi_avg = f["tally/flux-x/mean"][:]
-    phi_sd = f["tally/flux-x/sdev"][:]
-    k = f["k_cycle"][:]
-    k_avg = f["k_mean"][()]
-    k_sd = f["k_sdev"][()]
-    rg = f["gyration_radius"][:]
-    f.close()
-# Note the spatial (dx) and source strength (100+1) normalization
-tmp = 0.5 * (phi_avg[1:] + phi_avg[:-1])
-norm = np.sum(tmp * dx)
-phi_avg /= norm
-
-print("MC Keff = ", k_avg)
-plt.plot(x, phi_avg, label="MC")
-
-
-# =============================================================================
-# Power Iteration Results
-# =============================================================================
-
-with h5py.File("PI_output.h5", "r") as f:
-    # Note the spatial (dx) and source strength (100+1) normalization
-    keff = f["k_eff"][()]
-    phi_avg = f["tally/iqmc_flux"][:]
-    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("PI Keff = ", keff)
-plt.plot(x_mid, phi_avg, label="PI")
-
-# =============================================================================
-# Davidson Results
-# =============================================================================
-
-with h5py.File("davidson_output.h5", "r") as f:
-    # Note the spatial (dx) and source strength (100+1) normalization
-    keff = f["k_eff"][()]
-    phi_avg = f["tally/iqmc_flux"][:]
-    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("davidson Keff = ", keff)
-plt.plot(x_mid, phi_avg, label="davidson")
-
-
-# =============================================================================
-# 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

@@ -33,13 +33,14 @@
 # =============================================================================
 # iQMC Parameters
 # =============================================================================
-N = 1e1
+N = 1e4
 Nx = Ny = 40
-maxit = 1
-tol = 1e-3
+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"
 
 # fixed source in lower left corner
 fixed_source = np.zeros((Nx, Ny))
@@ -55,6 +56,8 @@
     maxitt=maxit,
     tol=tol,
     generator=generator,
+    fixed_source_solver=solver,
+    score=["tilt-x", "tilt-y"],
 )
 
 # =============================================================================