Merge pull request #273 from ilhamv/upgrade_time_census

Fix bug in time census

mcdc/kernel.py

@@ -2591,14 +2591,6 @@ def move_to_event(P_arr, data, mcdc):
     elif geometry.check_coincidence(d_collision, distance):
         P["event"] += EVENT_COLLISION
 
-    # Check distance to time boundary
-    if d_time_boundary < distance - COINCIDENCE_TOLERANCE:
-        distance = d_time_boundary
-        P["event"] = EVENT_TIME_BOUNDARY
-        P["surface_ID"] = -1
-    elif geometry.check_coincidence(d_time_boundary, distance):
-        P["event"] += EVENT_TIME_BOUNDARY
-
     # Check distance to time census
     if d_time_census < distance - COINCIDENCE_TOLERANCE:
         distance = d_time_census
@@ -2607,6 +2599,12 @@ def move_to_event(P_arr, data, mcdc):
     elif geometry.check_coincidence(d_time_census, distance):
         P["event"] += EVENT_TIME_CENSUS
 
+    # Check distance to time boundary (exclusive event)
+    if d_time_boundary < distance + COINCIDENCE_TOLERANCE:
+        distance = d_time_boundary
+        P["event"] = EVENT_TIME_BOUNDARY
+        P["surface_ID"] = -1
+
     # =========================================================================
     # Move particle
     # =========================================================================

mcdc/loop.py

@@ -98,12 +98,15 @@ def loop_fixed_source(data_arr, mcdc_arr):
                 )
                 loop_source_precursor(seed_source_precursor, data, mcdc)
 
-            # Time census closeout
-            if idx_census < mcdc["setting"]["N_census"] - 1:
-                # TODO: Output tally (optional)
+            # TODO: Output tally
 
-                # Manage particle banks: population control and work rebalance
-                seed_bank = kernel.split_seed(seed_census, SEED_SPLIT_BANK)
+            # No particles censused?
+            if mcdc["bank_census"]["size"] == 0:
+                break
+
+            # 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:

mcdc/main.py

@@ -1156,7 +1156,7 @@ def prepare():
     ):
         t_limit = INF
 
-    # Check if time boundary is above the final tally mesh time grid
+    # Set appropriate time boundary
     if mcdc["setting"]["time_boundary"] > t_limit:
         mcdc["setting"]["time_boundary"] = t_limit
 

pyproject.toml

@@ -48,7 +48,7 @@ classifiers = [
     'Operating System :: MacOS',
 ]
 dependencies = [
-    "numba>=0.58.0",
+    "numba>0.60.0",
     "numpy",
     "scipy",
     "matplotlib",

test/regression/azurv1_census/input.py

@@ -0,0 +1,50 @@
+import numpy as np
+from pprint import pprint
+import mcdc
+
+# =============================================================================
+# Set model
+# =============================================================================
+# Infinite medium with isotropic plane surface at the center
+# Based on Ganapol LA-UR-01-1854 (AZURV1 benchmark)
+# Effective scattering ratio c = 1.1
+
+# Set materials
+m = mcdc.material(
+    capture=np.array([1.0 / 3.0]),
+    scatter=np.array([[1.0 / 3.0]]),
+    fission=np.array([1.0 / 3.0]),
+    nu_p=np.array([2.3]),
+)
+
+# Set surfaces
+s1 = mcdc.surface("plane-x", x=-1e10, bc="reflective")
+s2 = mcdc.surface("plane-x", x=1e10, bc="reflective")
+
+# Set cells
+mcdc.cell(+s1 & -s2, m)
+
+# =============================================================================
+# Set source
+# =============================================================================
+# Isotropic pulse at x=t=0
+
+mcdc.source(point=[0.0, 0.0, 0.0], isotropic=True, time=[1e-10, 1e-10])
+
+# =============================================================================
+# Set tally, setting, and run mcdc
+# =============================================================================
+
+mcdc.tally.mesh_tally(
+    scores=["flux"],
+    x=np.linspace(-20.5, 20.5, 202),
+    t=np.linspace(0.0, 20.0, 21),
+)
+
+# Setting
+mcdc.setting(N_particle=30, census_bank_buff=10)
+mcdc.time_census(np.linspace(0.0, 20.0, 21)[1:-1])
+mcdc.population_control()
+
+# Run
+mcdc.run()

test/regression/run.py

@@ -118,16 +118,19 @@
 
     # Run the test problem (redirect the stdout)
     if mpiexec > 1:
+        gpus_per_task = ""
+        if target == "gpu":
+            gpus_per_task = f"--gpus-per-task=1 "
         os.system(
-            "mpiexec -n %i python input.py --mode=%s --target=%s --output=output --no-progress-bar > tmp 2>&1"
-            % (mpiexec, mode, target)
+            "mpiexec -n %i %s python input.py --mode=%s --target=%s --output=output --no-progress-bar> tmp 2>&1"
+            % (mpiexec, gpus_per_task, mode, target)
         )
     elif srun > 1:
         gpus_per_task = ""
         if target == "gpu":
             gpus_per_task = f"--gpus-per-task=1 "
         os.system(
-            "srun -n %i %s python input.py --mode=%s --target=%s --output=output --no-progress-bar %s> tmp 2>&1"
+            "srun -n %i %s python input.py --mode=%s --target=%s --output=output --no-progress-bar> tmp 2>&1"
             % (srun, gpus_per_task, mode, target)
         )
     else: