Forced collisions variance reduction technique

docs/source/pythonapi/index.rst

@@ -87,6 +87,9 @@ Defining techniques
 Techniques are enabled by calling methods on the ``mcdc.simulation`` singleton:
 
 - ``mcdc.simulation.implicit_capture(active=True)``
+- ``mcdc.simulation.forced_collisions(cells=[], weight_thresholds=[], weight_targets=[])``
+- ``mcdc.simulation.weighted_emission(active=True, weight_target=1.0)``
+- ``mcdc.simulation.weight_windows(weight_windows, mesh=None, energy=None)``
 - ``mcdc.simulation.global_weight_roulette(weight_threshold=0.0, weight_target=1.0)``
 - ``mcdc.simulation.population_control(active=True)``
 - ``mcdc.simulation.weighted_emission(active=True, weight_target=1.0)``

docs/source/theory/variance_reduction.rst

@@ -42,6 +42,59 @@ effective in highly absorbing media.
    to eliminate very-low-weight particles, a memory overhead can build
    up over time.
 
+
+Forced Collisions
+------------------
+
+.. warning::
+
+   Forced collisions are currently valid only for neutral particles.
+
+Forced collisions force a neutral particle to undergo a collision within selected material cells before it reaches the next surface.
+This is done by splitting the incident particle into two components:
+
+- a transmitted component that travels to the next surface without collision, with weight
+
+  .. math::
+
+     w_{\text{trans}} = w_0 e^{-\Sigma_t d}
+
+- a collided component that undergoes a forced collision in the cell, with weight
+
+  .. math::
+
+     w_{\text{coll}} = w_0 \left(1 - e^{-\Sigma_t d}\right)
+
+where :math:`d` is the distance to the next surface and :math:`\Sigma_t` is the total macroscopic cross section.
+The collision distance is sampled from the following distribution:
+
+.. math::
+
+   s = -\frac{1}{\Sigma_t} \ln \left(1 - \xi \left(1 - e^{-\Sigma_t d}\right)\right)
+
+Additionally, collided particles are continually forced to undergo collisions in the cell of interest.
+To prevent tracking particles with extremely low weight, weight roulette is used.
+
+**Usage:**
+
+.. code-block:: python3
+
+  mcdc.simulation.forced_collisions(cells=[cell])
+
+Optional roulette parameters can be supplied for each forced-collision cell:
+
+.. code-block:: python3
+
+  mcdc.simulation.forced_collisions(
+      cells=[cell_1, cell_2],
+      threshold_weights=[0.25, 0.25],
+      target_weights=[1.0, 1.0]
+  )
+
+If ``threshold_weights`` or ``target_weights`` are omitted, default values of ``0.5`` and ``1.0`` are assumed for each cell.
+Forced collisions may only be enabled on cells with material fills.
+
+
 Weight Roulette
 ----------------
 

mcdc/mcdc_get/__init__.py

@@ -86,6 +86,8 @@
 
 import mcdc.mcdc_get.settings as settings
 
+import mcdc.mcdc_get.forced_collisions as forced_collisions
+
 import mcdc.mcdc_get.global_weight_roulette as global_weight_roulette
 
 import mcdc.mcdc_get.implicit_capture as implicit_capture

mcdc/mcdc_get/forced_collisions.py

@@ -0,0 +1,90 @@
+# The following is automatically generated by code_factory.py
+
+from numba import njit
+
+
+@njit
+def cell_IDs(index, forced_collisions, data):
+    offset = forced_collisions["cell_IDs_offset"]
+    return data[offset + index]
+
+
+@njit
+def cell_IDs_all(forced_collisions, data):
+    start = forced_collisions["cell_IDs_offset"]
+    size = forced_collisions["cell_IDs_length"]
+    end = start + size
+    return data[start:end]
+
+
+@njit
+def cell_IDs_last(forced_collisions, data):
+    start = forced_collisions["cell_IDs_offset"]
+    size = forced_collisions["cell_IDs_length"]
+    end = start + size
+    return data[end - 1]
+
+
+@njit
+def cell_IDs_chunk(start, length, forced_collisions, data):
+    start += forced_collisions["cell_IDs_offset"]
+    end = start + length
+    return data[start:end]
+
+
+@njit
+def threshold_weights(index, forced_collisions, data):
+    offset = forced_collisions["threshold_weights_offset"]
+    return data[offset + index]
+
+
+@njit
+def threshold_weights_all(forced_collisions, data):
+    start = forced_collisions["threshold_weights_offset"]
+    size = forced_collisions["threshold_weights_length"]
+    end = start + size
+    return data[start:end]
+
+
+@njit
+def threshold_weights_last(forced_collisions, data):
+    start = forced_collisions["threshold_weights_offset"]
+    size = forced_collisions["threshold_weights_length"]
+    end = start + size
+    return data[end - 1]
+
+
+@njit
+def threshold_weights_chunk(start, length, forced_collisions, data):
+    start += forced_collisions["threshold_weights_offset"]
+    end = start + length
+    return data[start:end]
+
+
+@njit
+def target_weights(index, forced_collisions, data):
+    offset = forced_collisions["target_weights_offset"]
+    return data[offset + index]
+
+
+@njit
+def target_weights_all(forced_collisions, data):
+    start = forced_collisions["target_weights_offset"]
+    size = forced_collisions["target_weights_length"]
+    end = start + size
+    return data[start:end]
+
+
+@njit
+def target_weights_last(forced_collisions, data):
+    start = forced_collisions["target_weights_offset"]
+    size = forced_collisions["target_weights_length"]
+    end = start + size
+    return data[end - 1]
+
+
+@njit
+def target_weights_chunk(start, length, forced_collisions, data):
+    start += forced_collisions["target_weights_offset"]
+    end = start + length
+    return data[start:end]

mcdc/mcdc_set/__init__.py

@@ -86,6 +86,8 @@
 
 import mcdc.mcdc_set.settings as settings
 
+import mcdc.mcdc_set.forced_collisions as forced_collisions
+
 import mcdc.mcdc_set.global_weight_roulette as global_weight_roulette
 
 import mcdc.mcdc_set.implicit_capture as implicit_capture

mcdc/mcdc_set/forced_collisions.py

@@ -0,0 +1,90 @@
+# The following is automatically generated by code_factory.py
+
+from numba import njit
+
+
+@njit
+def cell_IDs(index, forced_collisions, data, value):
+    offset = forced_collisions["cell_IDs_offset"]
+    data[offset + index] = value
+
+
+@njit
+def cell_IDs_all(forced_collisions, data, value):
+    start = forced_collisions["cell_IDs_offset"]
+    size = forced_collisions["cell_IDs_length"]
+    end = start + size
+    data[start:end] = value
+
+
+@njit
+def cell_IDs_last(forced_collisions, data, value):
+    start = forced_collisions["cell_IDs_offset"]
+    size = forced_collisions["cell_IDs_length"]
+    end = start + size
+    data[end - 1] = value
+
+
+@njit
+def cell_IDs_chunk(start, length, forced_collisions, data, value):
+    start += forced_collisions["cell_IDs_offset"]
+    end = start + length
+    data[start:end] = value
+
+
+@njit
+def threshold_weights(index, forced_collisions, data, value):
+    offset = forced_collisions["threshold_weights_offset"]
+    data[offset + index] = value
+
+
+@njit
+def threshold_weights_all(forced_collisions, data, value):
+    start = forced_collisions["threshold_weights_offset"]
+    size = forced_collisions["threshold_weights_length"]
+    end = start + size
+    data[start:end] = value
+
+
+@njit
+def threshold_weights_last(forced_collisions, data, value):
+    start = forced_collisions["threshold_weights_offset"]
+    size = forced_collisions["threshold_weights_length"]
+    end = start + size
+    data[end - 1] = value
+
+
+@njit
+def threshold_weights_chunk(start, length, forced_collisions, data, value):
+    start += forced_collisions["threshold_weights_offset"]
+    end = start + length
+    data[start:end] = value
+
+
+@njit
+def target_weights(index, forced_collisions, data, value):
+    offset = forced_collisions["target_weights_offset"]
+    data[offset + index] = value
+
+
+@njit
+def target_weights_all(forced_collisions, data, value):
+    start = forced_collisions["target_weights_offset"]
+    size = forced_collisions["target_weights_length"]
+    end = start + size
+    data[start:end] = value
+
+
+@njit
+def target_weights_last(forced_collisions, data, value):
+    start = forced_collisions["target_weights_offset"]
+    size = forced_collisions["target_weights_length"]
+    end = start + size
+    data[end - 1] = value
+
+
+@njit
+def target_weights_chunk(start, length, forced_collisions, data, value):
+    start += forced_collisions["target_weights_offset"]
+    end = start + length
+    data[start:end] = value

mcdc/numba_types.py

@@ -561,6 +561,16 @@
     ('gpu_storage', int64),
 ])
 
+forced_collisions = into_dtype([
+    ('active', bool),
+    ('cell_IDs_offset', int64),
+    ('cell_IDs_length', int64),
+    ('threshold_weights_offset', int64),
+    ('threshold_weights_length', int64),
+    ('target_weights_offset', int64),
+    ('target_weights_length', int64),
+])
+
 global_weight_roulette = into_dtype([
     ('active', bool),
     ('weight_threshold', float64),
@@ -831,6 +841,7 @@ def set_simulation(N: dict):
         ('N_tracklength_tally', int64),
         ('settings', settings),
         ('implicit_capture', implicit_capture),
+        ('forced_collisions', forced_collisions),
         ('weighted_emission', weighted_emission),
         ('global_weight_roulette', global_weight_roulette),
         ('weight_windows', weight_windows),

mcdc/object_/simulation.py

@@ -3,6 +3,7 @@
 
 from mcdc.object_.technique import (
     ImplicitCapture,
+    ForcedCollisions,
     PopulationControl,
     GlobalWeightRoulette,
     WeightWindows,
@@ -81,6 +82,7 @@ class Simulation(ObjectSingleton):
 
     # Techniques
     implicit_capture: ImplicitCapture
+    forced_collisions: ForcedCollisions
     weighted_emission: WeightedEmission
     global_weight_roulette: GlobalWeightRoulette
     weight_windows: WeightWindows
@@ -167,6 +169,7 @@ def __init__(self):
 
         # Techniques
         self.implicit_capture = ImplicitCapture()
+        self.forced_collisions = ForcedCollisions()
         self.weighted_emission = WeightedEmission()
         self.global_weight_roulette = GlobalWeightRoulette()
         self.weight_windows = WeightWindows()

mcdc/object_/technique.py

@@ -1,8 +1,11 @@
+from typing import TYPE_CHECKING, Annotated
 from numpy.typing import NDArray
 import numpy as np
-from typing import Annotated
-from mcdc.constant import INF
+from mcdc.constant import FILL_MATERIAL, INF
 from mcdc.object_.base import ObjectSingleton
+
+if TYPE_CHECKING:
+    from mcdc.object_.cell import Cell
 from mcdc.object_.mesh import MeshBase, MeshUniform
 from mcdc.print_ import print_error
 
@@ -23,6 +26,47 @@ def __call__(self, active: bool = True):
         self.active = active
 
 
+# ======================================================================================
+# ForcedCollisions
+# ======================================================================================
+
+
+class ForcedCollisions(ObjectSingleton):
+    # Annotations for Numba mode
+    label: str = "forced_collisions"
+    active: bool
+
+    cell_IDs: list[np.int64]
+    threshold_weights: list[float]
+    target_weights: list[float]
+
+    def __init__(self):
+        self.active = False
+        self.cell_IDs = []
+        self.threshold_weights = []
+        self.target_weights = []
+
+    def __call__(self, cells, threshold_weights=None, target_weights=None):
+        self.active = True
+        if threshold_weights is None:
+            threshold_weights = [0.5] * len(cells)
+        if target_weights is None:
+            target_weights = [1.0] * len(cells)
+        if len(cells) != len(threshold_weights) or len(cells) != len(target_weights):
+            print_error(
+                f"Expected cells, threshold_weights, and target_weights to be the same size, but got {len(cells)}, {len(threshold_weights)}, and {len(target_weights)} instead"
+            )
+
+        for cell in cells:
+            if cell.fill_type != FILL_MATERIAL:
+                print_error(
+                    f"Invalid cell fill on cell: \n{cell}\nForced collision technique is only valid on cells with material fill"
+                )
+            self.cell_IDs.append(cell.ID)
+        self.threshold_weights = threshold_weights
+        self.target_weights = target_weights
+
+
 # ======================================================================================
 # Weighted emission
 # ======================================================================================