merged with upstream

.gitignore

@@ -14,7 +14,9 @@ __pycache__
 *.swp
 
 # Output
+*output.h5
 output*.h5
+*output.h5
 *.png
 *.mp4
 *.gif
@@ -35,3 +37,6 @@ tmp_*
 # test cache
 .pytest_cache
 pytestdebug.log
+
+# profiler
+*.prof

README.md

@@ -78,7 +78,7 @@ python input.py
 
 A more advanced input example that includes setting up multigroup (in energy and 
 delayed precursor) materials, lattice geometry, and continuous movements of 
-control rods is provided in `MCDC/example/fixed_source/C5G7-TDX`.
+control rods is provided in `MCDC/example/c5g7/3d/TDX`.
 
 ### Output
 

examples/fixed_source/C5G7-TDX/input.py → examples/c5g7/3d/TDX/input.py

@@ -9,7 +9,7 @@
 # =============================================================================
 
 # Load material data
-lib = h5py.File("c5g7.h5", "r")
+lib = h5py.File("../../c5g7_xs.h5", "r")
 
 
 # Setter
@@ -48,14 +48,18 @@ def set_mat(mat):
 # Control rod banks fractions
 #   All out: 0.0
 #   All in : 1.0
-cr1 = np.array([1.0, 1.0, 1.0, 0.84, 1.0])
+cr1 = np.array([1.0, 1.0, 1.0, 0.889, 1.0])
 cr1_t = np.array([0.0, 5.0, 10.0, 15.0, 15.0 + 1.0 - cr1[-2]])
-cr2 = np.array([1.0, 1.0, 0.0, 0.0, 0.25])
-cr2_t = np.array([0.0, 5.0, 10.0, 15.0, 15.25])
-cr3 = np.array([1.0, 1.0, 0.0, 0.0, 1.0])
-cr3_t = np.array([0.0, 5.0, 10.0, 15.0, 16.0])
-cr4 = np.array([0.5, 0.5])
-cr4_t = np.array([0.0, 20.0])
+
+cr2 = np.array([1.0, 1.0, 0.0, 0.0, 0.8])
+cr2_t = np.array([0.0, 5.0, 10.0, 15.0, 15.8])
+
+cr3 = np.array([0.75, 0.75, 1.0])
+cr3_t = np.array([0.0, 15.0, 15.25])
+
+cr4 = np.array([1.0, 1.0, 0.5, 0.5, 1.0])
+cr4_t = np.array([0.0, 5.0, 7.5, 15.0, 15.5])
+
 # Control rod banks interfaces
 cr1 = core_height * (0.5 - cr1)
 cr2 = core_height * (0.5 - cr2)
@@ -279,7 +283,7 @@ def set_mat(mat):
 # At highest energy
 
 energy = np.zeros(7)
-energy[-1] = 1.0
+energy[0] = 1.0
 
 source = mcdc.source(
     point=[pitch * 17 / 2, -pitch * 17 / 2, 0.0], time=[0.0, 15.0], energy=energy
@@ -290,19 +294,11 @@ def set_mat(mat):
 # =============================================================================
 
 # Tally
-x_grid = np.linspace(0.0, pitch * 17 * 3, 17 * 3 + 1)
-y_grid = np.linspace(-pitch * 17 * 3, 0.0, 17 * 3 + 1)
-z_grid = np.linspace(
-    -(core_height / 2 + refl_thick), (core_height / 2 + refl_thick), 102 + 17 * 2 + 1
-)
-x_grid = np.linspace(0.0, pitch * 17 * 3, 2)
-y_grid = np.linspace(-pitch * 17 * 3, 0.0, 2)
-z_grid = np.linspace(-(core_height / 2 + refl_thick), (core_height / 2 + refl_thick), 2)
 t_grid = np.linspace(0.0, 20.0, 201)
 mcdc.tally(scores=["fission"], t=t_grid)
 
 # Setting
-mcdc.setting(N_particle=1e5, active_bank_buff=10000)
+mcdc.setting(N_particle=1e2, active_bank_buff=1000)
 
 # Run
 mcdc.run()

examples/fixed_source/C5G7-TDX/process.py → examples/c5g7/3d/TDX/process.py

@@ -4,62 +4,38 @@
 from matplotlib import cm
 from matplotlib import colors
 
+
 # =============================================================================
 # Plot results
 # =============================================================================
 
-# Results
+# Get results
 with h5py.File("output.h5", "r") as f:
     fis_avg = f["tally/fission/mean"][:]
     fis_sd = f["tally/fission/sdev"][:]
     t = f["tally/grid/t"][:]
-    x = f["tally/grid/x"][:]
-    z = f["tally/grid/z"][:]
-dx = x[1] - x[0]
-dz = z[1] - z[0]
-dV = dx * dx * dz
-
 t_mid = 0.5 * (t[:-1] + t[1:])
 dt = t[1:] - t[:-1]
 
-# Normalize
-norm = np.sum(fis_avg[:, 0]) / dt[0]
-fis_avg /= norm
-fis_sd /= norm
-fis_ = np.zeros_like(fis_avg[0])
-fis__sd = np.zeros_like(fis_avg[0])
-for i in range(7):
-    fis_ += fis_avg[i]
-    fis__sd += np.square(fis_sd[i])
-fis__sd = np.sqrt(fis__sd)
-fis_ /= dt
-fis__sd /= dt
-
-# Reference
+# Get reference
 with h5py.File("reference.h5", "r") as f:
-    fis_avg_ref = f["tally/fission/mean"][:]
-    fis_sd_ref = f["tally/fission/sdev"][:]
+    fis_ref = f["tally/fission/mean"][:]
+    fis_ref_sd = f["tally/fission/sdev"][:]
 
 # Normalize
-norm = np.sum(fis_avg_ref[:, 0]) / dt[0]
-fis_avg_ref /= norm
-fis_sd_ref /= norm
-fis__ref = np.zeros_like(fis_avg_ref[0])
-fis__sd_ref = np.zeros_like(fis_avg_ref[0])
-for i in range(7):
-    fis__ref += fis_avg_ref[i]
-    fis__sd_ref += np.square(fis_sd_ref[i])
-fis__sd_ref = np.sqrt(fis__sd_ref)
-fis__ref /= dt
-fis__sd_ref /= dt
-
+norm = fis_ref[0] / dt[0]
+fis_avg /= norm * dt
+fis_sd /= norm * dt
+fis_ref /= norm * dt
+fis_ref_sd /= norm * dt
 
-plt.plot(t_mid, fis__ref, "-m", fillstyle="none", label="Ref. (MC)")
+# Plot
+plt.plot(t_mid, fis_ref, "-m", fillstyle="none", label="Ref. (MC)")
 plt.fill_between(
-    t_mid, fis__ref - fis__sd_ref, fis__ref + fis__sd_ref, alpha=0.2, color="m"
+    t_mid, fis_ref - fis_ref_sd, fis_ref + fis_ref_sd, alpha=0.2, color="m"
 )
-plt.plot(t_mid, fis_, "ok", fillstyle="none", label="MC")
-plt.fill_between(t_mid, fis_ - fis__sd, fis_ + fis__sd, alpha=0.2, color="k")
+plt.plot(t_mid, fis_avg, "ok", fillstyle="none", label="MC")
+plt.fill_between(t_mid, fis_avg - fis_sd, fis_avg + fis_sd, alpha=0.2, color="k")
 plt.yscale("log")
 plt.ylabel("Normalized fission rate")
 plt.xlabel("time [s]")