Merge pull request #118 from ilhamv/main

add batch loop for fixed source

examples/fixed_source/slab_absorbium/input.py

@@ -42,7 +42,7 @@
 )
 
 # Setting
-mcdc.setting(N_particle=1e4)
+mcdc.setting(N_particle=1e7)
 
 # Run
 mcdc.run()

mcdc/card.py

@@ -66,6 +66,7 @@ def reset(self):
         self.setting = {
             "tag": "Setting",
             "N_particle": 0,
+            "N_batch": 1,
             "rng_seed": 1,
             "time_boundary": INF,
             "progress_bar": True,

mcdc/constant.py

@@ -49,6 +49,7 @@
 RNG_MOD = nb.uint64(0x8000000000000000)
 
 # RNG splitter seeds
+SEED_SPLIT_CENSUS = nb.uint64(0x43454D654E54)
 SEED_SPLIT_SOURCE = nb.uint64(0x43616D696C6C65)
 SEED_SPLIT_SOURCE_PRECURSOR = nb.uint64(0x546F6464)
 SEED_SPLIT_BANK = nb.uint64(0x5279616E)

mcdc/input_.py

@@ -903,35 +903,37 @@ def setting(**kw):
             "setting parameter",
             key,
             [
-                "progress_bar",
                 "N_particle",
+                "N_batch",
                 "rng_seed",
                 "time_boundary",
-                "particle_tracker",
-                "save_input_deck",
+                "progress_bar",
                 "output_name",
+                "save_input_deck",
+                "particle_tracker",
+                "k_eff",
                 "source_file",
                 "IC_file",
                 "active_bank_buff",
                 "census_bank_buff",
-                "k_eff",
             ],
             False,
         )
 
     # Get keyword arguments
     N_particle = kw.get("N_particle")
-    time_boundary = kw.get("time_boundary")
+    N_batch = kw.get("N_batch")
     rng_seed = kw.get("rng_seed")
-    output = kw.get("output_name")
+    time_boundary = kw.get("time_boundary")
     progress_bar = kw.get("progress_bar")
+    output = kw.get("output_name")
+    save_input_deck = kw.get("save_input_deck")
+    particle_tracker = kw.get("particle_tracker")
     k_eff = kw.get("k_eff")
-    bank_active_buff = kw.get("active_bank_buff")
-    bank_census_buff = kw.get("census_bank_buff")
     source_file = kw.get("source_file")
-    particle_tracker = kw.get("particle_tracker")
-    save_input_deck = kw.get("save_input_deck")
     IC_file = kw.get("IC_file")
+    bank_active_buff = kw.get("active_bank_buff")
+    bank_census_buff = kw.get("census_bank_buff")
 
     # Check if setting card has been initialized
     card = mcdc.input_deck.setting
@@ -940,6 +942,10 @@ def setting(**kw):
     if N_particle is not None:
         card["N_particle"] = int(N_particle)
 
+    # Number of batches
+    if N_batch is not None:
+        card["N_batch"] = int(N_batch)
+
     # Time boundary
     if time_boundary is not None:
         card["time_boundary"] = time_boundary

mcdc/kernel.py

@@ -1465,8 +1465,12 @@ def score_closeout_history(score):
 def score_closeout(score, mcdc):
     N_history = mcdc["setting"]["N_particle"]
 
-    if mcdc["setting"]["mode_eigenvalue"]:
+    if mcdc["setting"]["N_batch"] > 1:
+        N_history = mcdc["setting"]["N_batch"]
+
+    elif mcdc["setting"]["mode_eigenvalue"]:
         N_history = mcdc["setting"]["N_active"]
+
     else:
         # MPI Reduce
         buff = np.zeros_like(score["mean"])

mcdc/loop.py

@@ -10,6 +10,7 @@
 
 from mcdc.constant import *
 from mcdc.print_ import (
+    print_header_batch,
     print_progress,
     print_progress_eigenvalue,
     print_progress_iqmc,
@@ -25,29 +26,45 @@
 
 @njit
 def loop_fixed_source(mcdc):
-    # Loop over time censuses
-    for idx_census in range(mcdc["setting"]["N_census"]):
-        mcdc["idx_census"] = idx_census
-        seed_census = kernel.split_seed(idx_census, mcdc["setting"]["rng_seed"])
+    # Loop over batches
+    for idx_batch in range(mcdc["setting"]["N_batch"]):
+        mcdc["idx_batch"] = idx_batch
+        seed_batch = kernel.split_seed(idx_batch, mcdc["setting"]["rng_seed"])
 
-        # Loop over source particles
-        seed_source = kernel.split_seed(seed_census, SEED_SPLIT_SOURCE)
-        loop_source(seed_source, mcdc)
-
-        # Loop over source precursors
-        if mcdc["bank_precursor"]["size"] > 0:
-            seed_source_precursor = kernel.split_seed(
-                seed_census, SEED_SPLIT_SOURCE_PRECURSOR
-            )
-            loop_source_precursor(seed_source_precursor, mcdc)
-
-        # Time census closeout
-        if idx_census < mcdc["setting"]["N_census"] - 1:
-            # TODO: Output tally (optional)
-
-            # Manage particle banks: population control and work rebalance
-            seed_bank = kernel.split_seed(seed_census, SEED_SPLIT_BANK)
-            kernel.manage_particle_banks(seed_bank, mcdc)
+        # Print multi-batch header
+        if mcdc["setting"]["N_batch"] > 1:
+            with objmode():
+                print_header_batch(mcdc)
+
+        # Loop over time censuses
+        for idx_census in range(mcdc["setting"]["N_census"]):
+            mcdc["idx_census"] = idx_census
+            seed_census = kernel.split_seed(seed_batch, SEED_SPLIT_CENSUS)
+
+            # Loop over source particles
+            seed_source = kernel.split_seed(seed_census, SEED_SPLIT_SOURCE)
+            loop_source(seed_source, mcdc)
+
+            # Loop over source precursors
+            if mcdc["bank_precursor"]["size"] > 0:
+                seed_source_precursor = kernel.split_seed(
+                    seed_census, SEED_SPLIT_SOURCE_PRECURSOR
+                )
+                loop_source_precursor(seed_source_precursor, mcdc)
+
+            # Time census closeout
+            if idx_census < mcdc["setting"]["N_census"] - 1:
+                # TODO: Output tally (optional)
+
+                # Manage particle banks: population control and work rebalance
+                seed_bank = kernel.split_seed(seed_census, SEED_SPLIT_BANK)
+                kernel.manage_particle_banks(seed_bank, mcdc)
+
+        # Multi-batch closeout
+        if mcdc["setting"]["N_batch"] > 1:
+            # Tally history closeout
+            kernel.tally_reduce_bin(mcdc)
+            kernel.tally_closeout_history(mcdc)
 
     # Tally closeout
     kernel.tally_closeout(mcdc)
@@ -60,8 +77,6 @@ def loop_fixed_source(mcdc):
 
 @njit
 def loop_eigenvalue(mcdc):
-    simulation_end = False
-
     # Loop over power iteration cycles
     for idx_cycle in range(mcdc["setting"]["N_cycle"]):
         seed_cycle = kernel.split_seed(idx_cycle, mcdc["setting"]["rng_seed"])
@@ -105,7 +120,9 @@ def loop_source(seed, mcdc):
     N_prog = 0
 
     # Loop over particle sources
-    for idx_work in range(mcdc["mpi_work_size"]):
+    work_start = mcdc["mpi_work_start"]
+    work_end = work_start + mcdc["mpi_work_size"]
+    for idx_work in range(work_start, work_end):
         seed_work = kernel.split_seed(idx_work, seed)
 
         # Particle tracker
@@ -158,8 +175,8 @@ def loop_source(seed, mcdc):
         # Closeout
         # =====================================================================
 
-        # Tally history closeout for fixed-source simulation
-        if not mcdc["setting"]["mode_eigenvalue"]:
+        # Tally history closeout for one-batch fixed-source simulation
+        if not mcdc["setting"]["mode_eigenvalue"] and mcdc["setting"]["N_batch"] == 1:
             kernel.tally_closeout_history(mcdc)
 
         # Progress printout

mcdc/print_.py

@@ -103,6 +103,12 @@ def print_header_eigenvalue(mcdc):
             print(" ====  =======  ===================")
 
 
+def print_header_batch(mcdc):
+    idx_batch = mcdc["idx_batch"]
+    if master:
+        print("\nBatch %i/%i" % (idx_batch + 1, mcdc["setting"]["N_batch"]))
+
+
 def print_progress_eigenvalue(mcdc):
     if master:
         idx_cycle = mcdc["idx_cycle"]

mcdc/type_.py

@@ -455,6 +455,7 @@ def make_type_setting(deck):
     struct = [
         # Basic MC simulation parameters
         ("N_particle", uint64),
+        ("N_batch", uint64),
         ("rng_seed", uint64),
         ("time_boundary", float64),
         # Misc.
@@ -731,6 +732,7 @@ def make_type_global(card):
             ("eigenvalue_tally_n", float64),
             ("eigenvalue_tally_C", float64),
             ("idx_census", int64),
+            ("idx_batch", int64),
             ("mpi_size", int64),
             ("mpi_rank", int64),
             ("mpi_master", bool_),

test/verification/analytic/fixed_source/inf_shem361/input.py

@@ -60,7 +60,7 @@
 mcdc.setting(
     N_particle=N_particle,
     output_name="output_" + str(N_particle),
-    active_bank_buff=1000,
+    active_bank_buff=10000,
     progress_bar=False,
 )
 

test/verification/analytic/fixed_source/inf_shem361/plot.py

@@ -6,9 +6,6 @@
 from reference import reference
 
 
-# Reference solution
-phi_ref = reference() / dE * E_mid
-
 # Load results
 output = sys.argv[1]
 with np.load("SHEM-361.npz") as data:
@@ -20,6 +17,9 @@
     phi = f["tally/flux/mean"][:] / dE * E_mid
     phi_sd = f["tally/flux/sdev"][:] / dE * E_mid
 
+# Reference solution
+phi_ref = reference() / dE * E_mid
+
 # Flux
 plt.step(E_mid, phi, "-b", label="MC", where="mid")
 plt.fill_between(E_mid, phi - phi_sd, phi + phi_sd, alpha=0.2, color="b", step="mid")

test/verification/analytic/fixed_source/inf_shem361_td/process.py

@@ -1,4 +1,3 @@
-from plotter import plot_convergence
 import numpy as np
 import h5py
 import sys
@@ -14,8 +13,8 @@
 
 # Reference solution
 data = np.load("reference.npz")
-phi_ref = data["phi"].T
-n_ref = data["n"].T
+phi_ref = data["phi"][1:].T
+n_ref = data["n"][1:].T
 
 # Error container
 error = np.zeros(len(N_particle_list))
@@ -25,7 +24,7 @@
 error_max_n = np.zeros(len(N_particle_list))
 
 # Calculate error
-for k, N_particle in enumerate(N_particle_list):
+for i, N_particle in enumerate(N_particle_list):
     # Get results
     with h5py.File("output_%i.h5" % int(N_particle), "r") as f:
         phi = f["tally/flux-t/mean"][:]
@@ -34,14 +33,15 @@
         speed = data["v"]
 
     # Neutron density
+    K = phi.shape[1]
     n = np.zeros(K)
     for k in range(K):
         n[k] = np.sum(phi[:, k] / speed)
 
-    error[k] = tool.error(phi, phi_ref)
-    error_n[k] = tool.error(n, n_ref)
-    error_max[k] = tool.error_max(phi, phi_ref)
-    error_max_n[k] = tool.error_max(n, n_ref)
+    error[i] = tool.error(phi, phi_ref)
+    error_n[i] = tool.error(n, n_ref)
+    error_max[i] = tool.error_max(phi, phi_ref)
+    error_max_n[i] = tool.error_max(n, n_ref)
 
 tool.plot_convergence("inf_shem361_td_flux", N_particle_list, error, error_max)
 tool.plot_convergence("inf_shem361_td_n", N_particle_list, error_n, error_max_n)

test/verification/analytic/fixed_source/reed/process.py

@@ -3,10 +3,11 @@
 import h5py
 import sys
 
-sys.path.append("../../util")
-from plotter import plot_convergence
+sys.path.append("../../")
+import tool
 
 
+# Cases run
 N_min = int(sys.argv[1])
 N_max = int(sys.argv[2])
 N_particle_list = np.logspace(N_min, N_max, (N_max - N_min) * 2 + 1)
@@ -17,17 +18,25 @@
     dx = x[1:] - x[:-1]
 phi_ref, phi_x_ref = reference(x)
 
-error = []
-error_x = []
+# Error containers
+error = np.zeros(len(N_particle_list))
+error_x = np.zeros(len(N_particle_list))
 
-for N_particle in N_particle_list:
+error_max = np.zeros(len(N_particle_list))
+error_max_x = np.zeros(len(N_particle_list))
+
+# Calculate error
+for k, N_particle in enumerate(N_particle_list):
     # Get results
     with h5py.File("output_%i.h5" % int(N_particle), "r") as f:
         phi = f["tally/flux/mean"][:] / dx * 101.0
         phi_x = f["tally/flux-x/mean"][:] * 101.0
 
-    error.append(np.linalg.norm((phi - phi_ref) / phi_ref))
-    error_x.append(np.linalg.norm((phi_x - phi_x_ref) / phi_x_ref))
+    # Get error
+    error[k] = tool.error(phi, phi_ref)
+    error_x[k] = tool.error(phi_x, phi_x_ref)
+    error_max[k] = tool.error_max(phi, phi_ref)
+    error_max_x[k] = tool.error_max(phi_x, phi_x_ref)
 
-plot_convergence("reed_flux", N_particle_list, error)
-plot_convergence("reed_flux_x", N_particle_list, error_x)
+tool.plot_convergence("reed_flux", N_particle_list, error, error_max)
+tool.plot_convergence("reed_flux_x", N_particle_list, error_x, error_max_x)

test/verification/analytic/tool.py

@@ -32,6 +32,6 @@ def plot_convergence(name, N_particle, error, error_max):
     plt.legend()
     plt.grid()
     plt.title(name)
-    plt.savefig("../../results/" + name + ".png")
+    # plt.savefig("../../results/" + name + ".png")
     plt.show()
     plt.clf()