Hybrid delta tracking

CMakeLists.txt

@@ -399,6 +399,7 @@ list(APPEND libopenmc_SOURCES
   src/ifp.cpp
   src/initialize.cpp
   src/lattice.cpp
+  src/majorant.cpp
   src/material.cpp
   src/math_functions.cpp
   src/mcpl_interface.cpp

docs/source/io_formats/settings.rst

@@ -236,6 +236,50 @@ from MF=1, MT=458 on an ENDF evaluation.
 
   *Default*: true
 
+----------------------------
+``<delta_tracking>``
+----------------------------
+
+The ``<delta_tracking>`` element lists the settings used for hybrid delta
+tracking. It has the following attributes/sub-elements:
+
+  :enable:
+    A boolean which dictates if OpenMC should run the simulation with delta
+    tracking or not. If false, surface tracking is used and the other
+    attributes/sub-elements of ``<delta_tracking>`` are ignored.
+
+  :hybrid_type:
+    The hybrid delta tracking scheme to use. Options are the
+    hybrid-in-cross-section scheme (``cross_section``) or the
+    hybrid-in-energy scheme (``energy``).
+
+    *Default*: ``cross_section``
+
+  :xs_threshold:
+    The hybrid-in-cross-section threshold (:math:`c`) used to determine whether
+    or not a particle should use delta tracking or surface tracking. When
+    :math:`\Sigma_{t} / \Sigma_{maj} > 1 - c`, the particle switches to delta
+    tracking. Otherwise, the particle switches to surface tracking. This parameter
+    is only used when ``hybrid_type`` is set to ``cross_section``.
+
+    *Default*: 0.9
+
+  :neutron_energy_threshold:
+    The energy threshold (eV) for neutrons in the hybrid-in-energy approach. When
+    the neutron energy is larger than ``neutron_energy_threshold``, delta
+    tracking is used. Otherwise, surface tracking is used. This parameter is
+    only valid when ``hybrid_type`` is set to ``energy``.
+
+    *Default*: 10
+
+  :photon_energy_threshold:
+    The energy threshold (eV) for photons in the hybrid-in-energy approach. When
+    the photon energy is larger than ``photon_energy_threshold``, delta
+    tracking is used. Otherwise, surface tracking is used. This parameter is
+    only valid when ``hybrid_type`` is set to ``energy``.
+
+    *Default*: 100000
+
 --------------------------------
 ``<electron_treatment>`` Element
 --------------------------------

docs/source/io_formats/statepoint.rst

@@ -28,6 +28,10 @@ The current version of the statepoint file format is 18.2.
              'continuous-energy' or 'multi-group'.
            - **run_mode** (*char[]*) -- Run mode used, either 'eigenvalue' or
              'fixed source'.
+           - **photon_transport** (*bool*) -- Whether photon transport was enabled
+             or not.
+           - **delta_tracking** (*bool*) -- Whether delta tracking was enabled
+             or not.
            - **n_particles** (*int8_t*) -- Number of particles used per generation.
            - **n_batches** (*int*) -- Number of batches to simulate.
            - **current_batch** (*int*) -- The number of batches already simulated.
@@ -178,6 +182,8 @@ All values are given in seconds and are measured on the master process.
              allocating arrays, etc.
            - **reading cross sections** (*double*) -- Time spent loading cross
              section libraries (this is a subset of initialization).
+           - **build majorant** (*double*) -- Time spent constructing the majorant
+             cross sections. This is only saved if running with delta tracking.
            - **simulation** (*double*) -- Time spent between initialization and
              finalization.
            - **transport** (*double*) -- Time spent transporting particles.

docs/source/methods/neutron_physics.rst

@@ -52,6 +52,156 @@ the formula usually used to calculate the distance to next collision is
 
     \ell = -\frac{\ln \xi}{\Sigma_t}
 
+.. _method_surface_tracking:
+
+~~~~~~~~~~~~~~~~~~~~~~~~~
+Surface Tracking
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The development of Equation :eq:`sample-distance-2` requires the assumption
+that the medium under consideration is homogeneous. To accomodate heterogeneous
+geometries, a resampling scheme is used. First, the distance to the next
+collision is sampled with Equation :eq:`sample-distance-2`. Then, the distance
+to the nearest surface from the particle position along its current trajectory
+is computed (discussed in the :ref:`methods_geometry` section). If the distance
+to the nearest surface is smaller than the distance to the next collision, the
+sampled distance is not statistically valid. The particle is moved to the
+surface and is considered to be contained by the next geometric region. Cross
+sections are recomputed, and a new distance to the next collision is sampled
+with Equation :eq:`sample-distance-2`. This process repeats until the distance
+to the next collision is smaller than the distance to the nearest surface,
+which is when a collision is accepted. This procedure is known as surface
+tracking.
+
+Surface tracking is quite efficient when used in problems with short mean free
+paths relative to the size of individual regions in the problem geometry.
+Surface tracking also admits the use of the track length estimator (discussed
+in the :ref:`methods_tallies` section). In problems with long mean free paths
+relative to the size of geometry regions, surface tracking will require a large
+number of surface distance calculations per collision. The cost of finding the
+nearest surface is also non-trivial for problems that contain many geometric
+regions at the same cell level (e.g. TRISO-fueled fission reactors).
+
+.. _method_delta_tracking:
+
+~~~~~~~~~~~~~~~~~~~~~~~~~
+Delta Tracking
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The disadvantages of surface tracking for certain classes of problems motivates
+the development of alternative approaches which do not require distance
+to surface checks. Delta tracking (also known by Woodcock tracking,
+delta scattering, and null scattering) is one approach to
+avoid surface geometry queries [Woodcock]_. In delta tracking, the domain is
+homogenized to a majorant cross section (:math:`\Sigma_{maj}(E)`) which is
+computed as the maximum total cross section over the entire problem
+
+.. math::
+    :label: majorant-xs-1
+
+    \Sigma_{maj}(E) = \max_{\mathbf{r}}\left(\Sigma_{t}(\mathbf{r}, E)\right).
+
+Particles move through this homogenized problem by sampling a distance to the
+next collision with the majorant cross section instead of the total cross section
+
+.. math::
+    :label: sample-distance-maj
+
+    \ell = -\frac{\ln \xi}{\Sigma_{maj}(E)}.
+
+To recover the spatial heterogeneity present in the original simulation, the
+delta tracking method formally defines the majorant cross section to be the
+sum of the total cross section and a spatially-varying fictitious delta
+scattering cross sections :math:`\Sigma_{\delta}(\mathbf{r}, E)`
+
+.. math::
+    :label: majorant-xs-2
+
+    \Sigma_{maj}(E) = \Sigma_{t}(\mathbf{r}, E) + \Sigma_{\delta}(\mathbf{r},
+    E).
+
+At the collision point computed from :eq:`sample-distance-maj` one of
+two reactions could occur. The first is a real collision, which is processed
+as usual. The second is known as a delta scatter event (also referred to
+as a virtual or null collision), which is the reaction type associated with
+:math:`\Sigma_{\delta}(\mathbf{r}, E)`. A rejection sampling test is used to
+determine which collision type occurs; a random number :math:`\xi` on the
+interval :math:`[0,1)` is drawn and used to check
+
+.. math::
+    :label: delta-real-collision
+
+    \xi < \frac{\Sigma_t (\mathbf{r}, E)}{\Sigma_{maj} (E)}.
+
+If the condition above is true, the collision is accepted as real. If the condition
+is false, a delta scatter event has occurred and the particle continues along
+its trajectory with the same energy and direction. Boundary conditions
+are applied by testing the distance to the nearest external boundary and
+comparing this to the distance sampled with Equation :eq:`sample-distance-maj`.
+If the distance to the nearest boundary is less than the sampled distance to
+the next collision, the particle crosses the external boundary.
+
+Delta tracking is advantageous as it only requires point location checks to
+determine the total cross section at each collision point to test
+:eq:`delta-real-collision`. This allows for the use of continuously-varying
+material properties and avoids computationally expensive distance-to-nearest-surface
+calculations. Problems which contain small regions with large total
+cross sections (such as burnable absorbers) will have majorant cross sections
+several orders of magnitude larger than the total cross section over the majority
+of the domain [Leppänen]_. This decreases the number of real collisions, and
+therefore the effectiveness of delta tracking. Material discontinuities are not
+considered in delta tracking, which prohibits the use of track length
+estimators for quantities restricted to material/geometric subdomains (such as
+reaction rates) and forces the use of the higher-variance collision
+estimator (discussed in detail in the :ref:`methods_tallies` section). When compared
+with surface tracking, delta tracking often performs better in problems where
+the particle mean free path is larger than the distance between surfaces.
+
+.. _method_hybrid_tracking:
+
+~~~~~~~~~~~~~~~~~~~~~~~~~
+Hybrid Tracking
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The complimentary performance bottlenecks of surface and delta tracking lead to the
+development of hybrid schemes that combined the tracking methodologies. The most
+common method is the approach first implemented in Serpent [Leppänen]_, which is
+referred to as the hybrid-in-cross-section approach. After every collision (real or
+virtual) the following inequality is tested:
+
+.. math::
+    :label: hybrid-cross-section-test
+
+    \frac{\Sigma_t (\mathbf{r}, E)}{\Sigma_{maj} (E)} > 1 - c,
+
+where :math:`c \in [0, 1]` is a user-provided hybrid-in-cross-section threshold.
+If equation :eq:`hybrid-cross-section-test` is true, the particle will use delta
+tracking for the next distance-to-collision calculation. Otherwise, the particle
+will use surface tracking. This ensures delta tracking is not used when the majorant
+cross section is sufficiently large and the efficiency of rejection sampling is
+reduced. The value of :math:`c` necessary to obtain optimal performance depends both
+on the specific model, and how the model is constructed. Problems constructed with
+lattices benefit more from delta tracking as point containment checks are cheap,
+and should use a value of :math:`c` between 0.8 and unity. Problems containing
+many cells per universe benefit more from surface tracking (which can be
+accelerated with cell neighborhood lists), and should use a value close to zero.
+
+An alternative hybrid tracking approach is the hybrid-in-energy scheme [Morgan]_.
+After every collision the particle energy is compared to an energy
+threshold :math:`E_{th}`. If the energy is greater than :math:`E_{th}`, the particle
+will use delta tracking for the next distance-to-collision calculation. If the energy
+is less than or equal to :math:`E_{th}` the particle will use surface tracking for
+the next distance-to-collision calculation. A sufficiently large value of
+:math:`E_{th}` ensures that particles do not use delta tracking in the low
+energy domain, where the majorant cross section is inflated by absorbers.
+Most neutron transport problems should use a value for neutrons greater than
+1 eV to avoid the thermal region where absorbers decrease the efficiency of
+delta tracking. Photons should use a value greater than 100 keV to avoid
+the region where the Compton cross section dominates, resulting in a
+decrease in the efficiency of delta tracking. Otherwise, the specific value
+of :math:`E_{th}` that maximizes performance depends on the problem and
+its construction in a similar manner to the hybrid-in-cross-section method.
+
 ----------------------------------------------------
 :math:`(n,\gamma)` and Other Disappearance Reactions
 ----------------------------------------------------
@@ -1631,6 +1781,13 @@ the unresolved range to get the actual cross sections. Lastly, the total cross
 section is calculated as the sum of the elastic, fission, capture, and inelastic
 cross sections.
 
+Unresolved resonance probability tables pose a challenge when computing a majorant
+cross section for :ref:`method_delta_tracking`. OpenMC implements a conservative approach:
+the maximum total cross section is computed over all bands, which is then
+interpolated to the corresponding energy using either linear or logarithmic
+interpolation. This ensures the majorant bounds the total cross section at the
+cost of increasing the number of delta scatters in the unresolved range.
+
 -----------------------------
 Variance Reduction Techniques
 -----------------------------
@@ -1735,13 +1892,26 @@ types.
 .. [Gelbard] Ely M. Gelbard, "Epithermal Scattering in VIM," FRA-TM-123, Argonne
    National Laboratory (1979).
 
+.. [Leppänen] J. Leppänen. "Performance of Woodcock Delta-Tracking in Lattice
+   Physics Applications using the Serpent Monte Carlo Reactor Physics Burnup
+   Calculation Code", *Annals of Nuclear Energy*, 37:715-722, 2010.
+
+.. [Morgan] J. P. Morgan, I. Variansyah, K. B. Clements, T. S. Palmer,
+   and K. E. Niemeyer. "Hybrid Delta Tracking Schemes Using a Track-Length
+   Estimator", *Journal of Computational and Theoretical Transport*, 2026.
+
 .. [Squires] G. L. Squires, *Introduction to the Theory of Thermal Neutron
    Scattering*, Cambridge University Press (1978).
 
 .. [Williams] M. M. R. Williams, *The Slowing Down and Thermalization of
    Neutrons*, North-Holland Publishing Co., Amsterdam (1966). **Note:** This
    book can be obtained for free from the OECD_.
 
+.. [Woodcock] E.R. Woodcock, T. Murphy, P.J. Hemmings, and T.C. Longworth.
+   "Techniques used in the GEM Code for Monte Carlo Neutronics Calculations
+   in Reactors and other Systems of Complex Geometry", ANL-7050,
+   Argonne National Laboratory (1965).
+
 .. |sab| replace:: S(:math:`\alpha,\beta,T`)
 
 .. _SIGMA1 method: https://doi.org/10.13182/NSE76-1

docs/source/usersguide/settings.rst

@@ -52,6 +52,79 @@ would need to instantiate a :class:`openmc.Settings` object and assign the
 
 If you don't specify a run mode, the default run mode is 'eigenvalue'.
 
+.. _usersguide_tracking:
+
+-------------------
+Tracking Algorithms
+-------------------
+
+By default, OpenMC and many other Monte Carlo transport solvers use
+:ref:`surface tracking <method_surface_tracking>` to sample distances
+between collisions as it allows for the use of the low-variance tracklength
+estimator. Surface tracking often performs poorly in problems that contain
+many repeated unit cells, as particles must determine the distance to the
+nearest surface many times between collisions. Examples of problems which
+often see this performance degradation are fission reactors, namely fast
+reactors (where neutron mean free paths are long) and reactors with very small
+dispersed heterogeneities (such as TRISO-fueled reactors).
+
+To remedy the performance penalty of using surface tracking in these systems,
+OpenMC includes an implementation of
+:ref:`hybrid tracking <method_hybrid_tracking>`. This family of techniques
+is based on :ref:`delta tracking <method_delta_tracking>`, which replaces
+distance to nearest surface calculations with "where am I" lookups that
+are extremely fast in lattice geometries. The performance of delta tracking
+degrades when problems contain small localized absorbers, and so hybrid
+tracking mixes surface tracking and delta tracking to maximize performance.
+The two approaches implemented in OpenMC are the hybrid-in-cross-section
+method (identical to the approach used by Serpent) and the hybrid-in-energy
+method. To use hybrid-in-cross-section tracking, you need to set
+:attr:`Settings.delta_tracking` to the following::
+
+    settings.delta_tracking = {
+      'enable' : True,
+      'hybrid_type' : 'cross_section'
+    }
+
+The efficiency of this hybrid scheme is determined by how often delta tracking
+is run relalative to surface tracking, which is controlled by the unitless
+hybrid tracking threshold :math:`0 \leq c \leq 1`. This parameter is roughly
+equivalent to the probability that a collision is *rejected* by delta tracking
+at a given point in space. When that probability is high, the cost of delta
+tracking becomes large and surface tracking is comparatively more efficient.
+Accordingly, a value of :math:`c = 0` will run pure surface tracking,
+:math:`c = 1` will run pure delta tracking, and a value between 0 and 1 will
+run a combination of the two to bound the rejection probability. A value of
+:math:`c` between 0.8 and 0.95 is often recommended for the vast majority
+of fission reactor problems based on parameter studies. OpenMC uses a default
+of :math:`c = 0.9`; this can be changed by setting::
+
+    settings.delta_tracking['xs_threshold'] = 0.95
+
+The hybrid-in-energy-approach choses when surface tracking or delta tracking
+should be used based on a particle-specific energy threshold :math:`E_{th}`
+(in units of eV). If the energy of a particle is greater than :math:`E_{th}`
+delta tracking is used, otherwise surface tracking is used. The use of
+particle energy as a hybrid parameter is based on the observation that
+neutrons collide less often at energies at high energies (passing through
+several lattice cells), and so delta tracking will be more efficient in
+those regimes. The hybrid-in-energy method has been found to be up to
+1.5x faster than the hybrid-in-cross-section method in OpenMC at the
+cost of requiring a particle-dependent threshold. To use hybrid-in-energy
+tracking, set :attr:`Settings.delta_tracking` to the following::
+
+    settings.delta_tracking = {
+      'enable' : True,
+      'hybrid_type' : 'energy'
+    }
+
+The default value of :math:`E_{th}` in OpenMC are 10 eV for neutrons and
+100 keV for photons, which were determined based on fission reactor
+parameter studies. These thresholds can be adjusted, in units of eV, as::
+
+    settings.delta_tracking['neutron_energy_threshold'] = 15
+    settings.delta_tracking['photon_energy_threshold'] = 1e6
+
 .. _usersguide_particles:
 
 ------------