Merge pull request #116 from ilhamv/main

Refactoring loops

examples/eigenvalue/2d_c5g7_iqmc/input.py

@@ -205,7 +205,7 @@ def set_mat(mat):
 
 
 # Setting
-mcdc.setting(N_particle=N, output="davidson_output")
+mcdc.setting(N_particle=N, output_name="davidson_output")
 mcdc.eigenmode()
 
 # Run

examples/fixed_source/slab_absorbium/input.py

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

mcdc/__init__.py

@@ -14,7 +14,7 @@
     weighted_emission,
     population_control,
     branchless_collision,
-    census,
+    time_census,
     weight_window,
     iQMC,
     weight_roulette,

mcdc/card.py

@@ -1,6 +1,6 @@
 import numpy as np
 
-from mcdc.constant import INF, GR_ALL, PCT_NONE, PI, SHIFT
+from mcdc.constant import INF, GYRATION_RADIUS_ALL, PCT_NONE, PI, SHIFT
 
 
 class InputDeck:
@@ -66,27 +66,31 @@ def reset(self):
         self.setting = {
             "tag": "Setting",
             "N_particle": 0,
+            "rng_seed": 1,
+            "time_boundary": INF,
+            "progress_bar": True,
+            "output_name": "output",
+            "save_input_deck": True,
+            "track_particle": False,
+            "mode_eigenvalue": False,
+            "k_init": 1.0,
             "N_inactive": 0,
             "N_active": 0,
             "N_cycle": 0,
-            "mode_eigenvalue": False,
-            "time_boundary": INF,
-            "rng_seed": 1,
-            "bank_active_buff": 100,
-            "bank_census_buff": 1.0,
-            "k_init": 1.0,
-            "output": "output",
-            "progress_bar": True,
+            "save_particle": False,
             "gyration_radius": False,
-            "gyration_radius_type": GR_ALL,
+            "gyration_radius_type": GYRATION_RADIUS_ALL,
+            "N_census": 1,
+            "census_time": np.array([INF]),
             "source_file": False,
             "source_file_name": "",
-            "track_particle": False,
-            "save_particle": False,
-            "save_input_deck": True,
             "IC_file": False,
             "IC_file_name": "",
             "N_precursor": 0,
+            # Below are parameters not copied to mcdc.setting
+            "bank_active_buff": 100,
+            "bank_census_buff": 1.0,
+            # TODO: Move to technique
             "N_sensitivity": 0,
         }
 
@@ -138,8 +142,6 @@ def reset(self):
             "IC_neutron_density_max": 0.0,
             "IC_precursor_density_max": 0.0,
             "IC_cycle_stretch": 1.0,
-            "time_census": False,
-            "census_time": np.array([INF]),
             "branchless_collision": False,
             "dsm_order": 1,
         }

mcdc/constant.py

@@ -22,13 +22,13 @@
 MESH_T = 3
 
 # Gyration raius type
-GR_ALL = 0
-GR_INFINITE_X = 1
-GR_INFINITE_Y = 2
-GR_INFINITE_Z = 3
-GR_ONLY_X = 4
-GR_ONLY_Y = 5
-GR_ONLY_Z = 6
+GYRATION_RADIUS_ALL = 0
+GYRATION_RADIUS_INFINITE_X = 1
+GYRATION_RADIUS_INFINITE_Y = 2
+GYRATION_RADIUS_INFINITE_Z = 3
+GYRATION_RADIUS_ONLY_X = 4
+GYRATION_RADIUS_ONLY_Y = 5
+GYRATION_RADIUS_ONLY_Z = 6
 
 # Population control
 PCT_NONE = 0

mcdc/input_.py

@@ -19,13 +19,13 @@
     make_card_source,
 )
 from mcdc.constant import (
-    GR_ALL,
-    GR_INFINITE_X,
-    GR_INFINITE_Y,
-    GR_INFINITE_Z,
-    GR_ONLY_X,
-    GR_ONLY_Y,
-    GR_ONLY_Z,
+    GYRATION_RADIUS_ALL,
+    GYRATION_RADIUS_INFINITE_X,
+    GYRATION_RADIUS_INFINITE_Y,
+    GYRATION_RADIUS_INFINITE_Z,
+    GYRATION_RADIUS_ONLY_X,
+    GYRATION_RADIUS_ONLY_Y,
+    GYRATION_RADIUS_ONLY_Z,
     PCT_NONE,
     PCT_COMBING,
     PCT_COMBING_WEIGHT,
@@ -900,21 +900,21 @@ def setting(**kw):
     # Check the suplied keyword arguments
     for key in kw.keys():
         check_support(
-            "source parameter",
+            "setting parameter",
             key,
             [
+                "progress_bar",
                 "N_particle",
-                "time_boundary",
                 "rng_seed",
-                "output",
-                "progress_bar",
-                "k_eff",
-                "active_bank_buff",
-                "census_bank_buff",
-                "source_file",
+                "time_boundary",
                 "particle_tracker",
                 "save_input_deck",
+                "output_name",
+                "source_file",
                 "IC_file",
+                "active_bank_buff",
+                "census_bank_buff",
+                "k_eff",
             ],
             False,
         )
@@ -923,7 +923,7 @@ def setting(**kw):
     N_particle = kw.get("N_particle")
     time_boundary = kw.get("time_boundary")
     rng_seed = kw.get("rng_seed")
-    output = kw.get("output")
+    output = kw.get("output_name")
     progress_bar = kw.get("progress_bar")
     k_eff = kw.get("k_eff")
     bank_active_buff = kw.get("active_bank_buff")
@@ -950,7 +950,7 @@ def setting(**kw):
 
     # Output .h5 file name
     if output is not None:
-        card["output"] = output
+        card["output_name"] = output
 
     # Progress bar
     if progress_bar is not None:
@@ -1020,19 +1020,19 @@ def eigenmode(
     if gyration_radius is not None:
         card["gyration_radius"] = True
         if gyration_radius == "all":
-            card["gyration_radius_type"] = GR_ALL
+            card["gyration_radius_type"] = GYRATION_RADIUS_ALL
         elif gyration_radius == "infinite-x":
-            card["gyration_radius_type"] = GR_INFINITE_X
+            card["gyration_radius_type"] = GYRATION_RADIUS_INFINITE_X
         elif gyration_radius == "infinite-y":
-            card["gyration_radius_type"] = GR_INFINITE_Y
+            card["gyration_radius_type"] = GYRATION_RADIUS_INFINITE_Y
         elif gyration_radius == "infinite-z":
-            card["gyration_radius_type"] = GR_INFINITE_Z
+            card["gyration_radius_type"] = GYRATION_RADIUS_INFINITE_Z
         elif gyration_radius == "only-x":
-            card["gyration_radius_type"] = GR_ONLY_X
+            card["gyration_radius_type"] = GYRATION_RADIUS_ONLY_X
         elif gyration_radius == "only-y":
-            card["gyration_radius_type"] = GR_ONLY_Y
+            card["gyration_radius_type"] = GYRATION_RADIUS_ONLY_Y
         elif gyration_radius == "only-z":
-            card["gyration_radius_type"] = GR_ONLY_Z
+            card["gyration_radius_type"] = GYRATION_RADIUS_ONLY_Z
         else:
             print_error("Unknown gyration radius type")
 
@@ -1075,12 +1075,21 @@ def branchless_collision():
     card["weighted_emission"] = False
 
 
-def census(t, pct="none"):
-    card = mcdc.input_deck.technique
-    card["time_census"] = True
+def time_census(t):
+    # Remove census beyond the final tally time grid point
+    while True:
+        if t[-1] >= mcdc.input_deck.tally["mesh"]["t"][-1]:
+            t = t[:-1]
+        else:
+            break
+
+    # Add the default, final census-at-infinity
+    t = np.append(t, INF)
+
+    # Set the time census parameters
+    card = mcdc.input_deck.setting
     card["census_time"] = t
-    if pct != "none":
-        population_control(pct)
+    card["N_census"] = len(t)
 
 
 def weight_window(x=None, y=None, z=None, t=None, window=None, width=None):

mcdc/kernel.py

@@ -1561,7 +1561,7 @@ def eigenvalue_tally(P, distance, mcdc):
 def eigenvalue_tally_closeout_history(mcdc):
     N_particle = mcdc["setting"]["N_particle"]
 
-    i_cycle = mcdc["i_cycle"]
+    idx_cycle = mcdc["idx_cycle"]
 
     # MPI Allreduce
     buff_nuSigmaF = np.zeros(1, np.float64)
@@ -1592,7 +1592,7 @@ def eigenvalue_tally_closeout_history(mcdc):
 
     # Update and store k_eff
     mcdc["k_eff"] = buff_nuSigmaF[0] / N_particle
-    mcdc["k_cycle"][i_cycle] = mcdc["k_eff"]
+    mcdc["k_cycle"][idx_cycle] = mcdc["k_eff"]
 
     # Normalize other eigenvalue/global tallies
     tally_n = buff_n[0] / N_particle
@@ -1612,7 +1612,7 @@ def eigenvalue_tally_closeout_history(mcdc):
         mcdc["C_avg"] += tally_C
         mcdc["C_sdv"] += tally_C * tally_C
 
-        N = 1 + mcdc["i_cycle"] - mcdc["setting"]["N_inactive"]
+        N = 1 + mcdc["idx_cycle"] - mcdc["setting"]["N_inactive"]
         mcdc["k_avg_running"] = mcdc["k_avg"] / N
         if N == 1:
             mcdc["k_sdv_running"] = 0.0
@@ -1658,35 +1658,35 @@ def eigenvalue_tally_closeout_history(mcdc):
         rms_local = np.zeros(1, np.float64)
         rms = np.zeros(1, np.float64)
         gr_type = mcdc["setting"]["gyration_radius_type"]
-        if gr_type == GR_ALL:
+        if gr_type == GYRATION_RADIUS_ALL:
             for i in range(N_local):
                 P = mcdc["bank_census"]["particles"][i]
                 rms_local[0] += (
                     (P["x"] - com_x) ** 2
                     + (P["y"] - com_y) ** 2
                     + (P["z"] - com_z) ** 2
                 ) * P["w"]
-        elif gr_type == GR_INFINITE_X:
+        elif gr_type == GYRATION_RADIUS_INFINITE_X:
             for i in range(N_local):
                 P = mcdc["bank_census"]["particles"][i]
                 rms_local[0] += ((P["y"] - com_y) ** 2 + (P["z"] - com_z) ** 2) * P["w"]
-        elif gr_type == GR_INFINITE_Y:
+        elif gr_type == GYRATION_RADIUS_INFINITE_Y:
             for i in range(N_local):
                 P = mcdc["bank_census"]["particles"][i]
                 rms_local[0] += ((P["x"] - com_x) ** 2 + (P["z"] - com_z) ** 2) * P["w"]
-        elif gr_type == GR_INFINITE_Z:
+        elif gr_type == GYRATION_RADIUS_INFINITE_Z:
             for i in range(N_local):
                 P = mcdc["bank_census"]["particles"][i]
                 rms_local[0] += ((P["x"] - com_x) ** 2 + (P["y"] - com_y) ** 2) * P["w"]
-        elif gr_type == GR_ONLY_X:
+        elif gr_type == GYRATION_RADIUS_ONLY_X:
             for i in range(N_local):
                 P = mcdc["bank_census"]["particles"][i]
                 rms_local[0] += ((P["x"] - com_x) ** 2) * P["w"]
-        elif gr_type == GR_ONLY_Y:
+        elif gr_type == GYRATION_RADIUS_ONLY_Y:
             for i in range(N_local):
                 P = mcdc["bank_census"]["particles"][i]
                 rms_local[0] += ((P["y"] - com_y) ** 2) * P["w"]
-        elif gr_type == GR_ONLY_Z:
+        elif gr_type == GYRATION_RADIUS_ONLY_Z:
             for i in range(N_local):
                 P = mcdc["bank_census"]["particles"][i]
                 rms_local[0] += ((P["z"] - com_z) ** 2) * P["w"]
@@ -1697,7 +1697,7 @@ def eigenvalue_tally_closeout_history(mcdc):
         rms = math.sqrt(rms[0] / W)
 
         # Gyration radius
-        mcdc["gyration_radius"][i_cycle] = rms
+        mcdc["gyration_radius"][idx_cycle] = rms
 
 
 @njit
@@ -1743,8 +1743,8 @@ def move_to_event(P, mcdc):
     d_time_boundary = speed * (mcdc["setting"]["time_boundary"] - P["t"])
 
     # Distance to census time
-    idx = mcdc["technique"]["census_idx"]
-    d_time_census = speed * (mcdc["technique"]["census_time"][idx] - P["t"])
+    idx = mcdc["idx_census"]
+    d_time_census = speed * (mcdc["setting"]["census_time"][idx] - P["t"])
 
     # Distance to collision
     if mcdc["technique"]["iQMC"]: