Hybrid delta tracking

[nuclearkevin]

Description

This PR adds hybrid tracking on top of #3971 to mitigate performance deficiencies of delta tracking in problems with localized absorbers. The first hybrid scheme is the method first implemented in Serpent, which uses the local ratio of the total to majorant cross section to determine whether delta or surface tracking should be used. This is done with a user-specified threshold which bounds the minimum rejection sampling efficiency of delta tracking. This is referred to as hybrid-in-cross-section tracking in OpenMC.

The second method added is the hybrid-in-energy scheme, which compares the particle energy to a given cutoff. Particles with an energy greater than the cutoff use delta tracking, while particles with an energy less than the cutoff use surface tracking. I chose to implement two hybrid schemes as the hybrid-in-energy approach has been found to yield better performance in OpenMC at the cost of requiring more user intuition, while the hybrid-in-cross-section approach is more familiar due to its use in other Monte Carlo codes.

Implementation Details

Hybrid tracking replaces the history and event based delta tracking implementation. If the particle is flagged for delta tracking the distance to the nearest external boundary is computed, the distance to the next collision is sampled from the majorant, and a call to global point containment is performed after advancing the particle. If the particle collisions before crossing an external boundary cross sections are computed and the collision is rejection sampled. If the particle is flagged for surface tracking the distance to the nearest surface is computed and the distance to next collision is sampled from the local total cross section. Then, if the surface distance is less than the collision distance the surface is crossed and the process repeats. Otherwise, a collision is processed. After every real or null collision the delta tracking flag is set to either true or false depending on the behaviour of the chosen hybrid scheme. Even if the particle is running surface tracking in a given iteration, the majorant is kept in sync as it may be used to change to delta tracking.

To better facilitate the implementation of event-based hybrid tracking, a new coarse-grained particle event for computing the majorant was introduced. Event-based hybrid tracking was then implemented such that it yields identical particle tracks to history-based transport. Event-based mode is implemented for the sake of feature-parity, and is likely not optimal.

Enabling hybrid tracking and specifying a hybrid scheme (and the associated parameters) is done by setting the delta_tracking dictionary in the openmc.Settings object. For the hybrid-in-cross-section scheme, this looks like the following:

settings.delta_tracking = {
  'enable' : True,
  'hybrid_type' : 'cross_section',
  'xs_threshold' : 0.9
}

Here, xs_threshold is the hybrid parameter which has a default value of 0.9 based on parameter studies (shown below). Setting xs_threshold to 0.0 will result in surface tracking, setting it to 1.0 will result in delta tracking. For the hybrid-in-energy scheme, the delta_tracking dictionary looks like:

settings.delta_tracking = {
  'enable' : True,
  'hybrid_type' : 'energy',
  'neutron_energy_threshold' : 10,
  'photon_energy_threshold' : 1e5
}

neutron_energy_threshold and photon_energy_threshold are the neutron and photon threshold energies for delta tracking (in units of eV). Separate energies are required for different particles as the performance of this hybrid scheme is very sensitive to the variation of the majorant and total cross sections with respect to energy. Defaults are set to 10 eV (neutrons) and 100 keV (photons) based on parameter studies (shown below).

Testing

All existing delta tracking tests become tests for hybrid tracking when using the hybrid-in-cross-section scheme and a hybrid parameter of 1.0 (a regold was not required as particle histories are the same). Eight new tests have been added, four for each hybrid scheme. Each test is a permutation of (history vs event) and (neutron vs neutron+photon).

Verification and Performance Results

Verification of hybrid tracking was performed in an identical manner to standard delta tracking. This included comparisons of integral metrics, spectra, and spatial fields (flux and total reaction rates). I've included whole-core integral comparisons here as they compliment the performance results. If people are interested in the other comparisons I can add them as requested.

To determine default values of the different hybrid parameters, parametric studies are performed on four different large-scale models:

1. The Gas-Cooled Microreactor (GCMR). This is a TRISO-fueled (40% packing fraction) prismatic HTGR with gadolinium poison particles (25% packing fraction), boron carbide control drums, and boron carbide control rods. The model is run with 10000 particles per batch, 100 inactive batches, and 1000 active batches.
2. BEAVRS. This is a prototypic whole-core PWR using oxide fuel of varying levels of enrichment. Absorbers include borosilicate glass, boron carbide control rods, and cadmium control rods. The model is run with 1000000 particles per batch, 100 inactive batches, and 1000 active batches.
3. ARC-100. This is a prototypic small SFR using metallic fuel of various levels of enrichment. Boron carbide control rods are the only absorbers in the geometry. The model is run with 1000000 particles per batch, 100 inactive batches, and 1000 active batches.
4. Octomak. This is a simple tokamak constructed entirely with CSG to test the performance of hybrid tracking when universes are flat as opposed to hierarchical/lattice-based. The model is run with 1000000 particles per batch and 100 active batches.

All simulations are performed on a single node of Improv @ ANL, where 8 MPI ranks and 16 threads per rank are used (ranks are explicitly bound to NUMA domains). The total runtimes of each hybrid tracking simulation (using collision estimators) is normalized to the runtime of an equivalent surface tracking simulation which uses tracklength estimators. The total runtimes are used for comparisons as it includes the time spent calculating the majorant, which can be significant. We start by varying the hybrid-in-cross-section threshold between 0.0 and 1.0 for both neutron and neutron+photon transport simulations, the results of which can be found below. We see maximum speedups of ~3.5x for the ARC-100, ~1.8x for the GCMR and ~1.6x for BEAVRS. Photon transport does not benefit much from hybrid-in-cross-section tracking as photons do not travel over multiple lattice elements (unlike neutrons). The Octomak problem sees a severe performance decrease when hybrid tracking is used, which is caused by the increase in relative cost of point location calls compared to next surface searches (which are accelerated with neighbourhood lists). In general, the behaviour found here mirrors the behaviour reported in literature.

Similarly, the hybrid-in-energy thresholds for neutron transport and neutron+photon transport are varied between 1 eV and 5 MeV. For the coupled transport problems, the neutron energy threshold which yielded the largest speedup is used while the photon transport energy threshold is varied. We see maximum speedup of ~3.5x for the ARC-100, ~3x for the GCMR, and ~1.6x for BEAVRS. Photon transport does not benefit much from hybrid-in-energy tracking outside of the GCMR, which sees a ~2.3x speedup (likely due to photons above the energy threshold being able to cross over multiple TRISO particles before colliding). As expected, the Octomak also does not perform well with the hybrid-in-energy for the same reason as the hybrid-in-cross-section method.

We compare the true runtimes, particle rates, and integral metrics below for the fission reactor problems at the optimal values of the hybrid parameters (as determined by the parameter study above). Integral quantities agree across the board to $3\sigma$.

GCMR

QoI	Hybrid-in-Cross-Section	Hybrid-in-Energy	Surface Tracking
Total Runtime (min)	83.6	51.7	154.8
Inactive Rate (p/s)	2272	3673	1200
Active Rate (p/s)	2279	3783	1186
k-eff (combined estimator)	1.01751 +/- 0.00030	1.01722 +/- 0.00031	1.01748 +/- 0.00029
Leakage Fraction	0.08880 +/- 0.00010	0.08908 +/- 0.00010	0.08863 +/- 0.00010

BEAVRS

QoI	Hybrid-in-Cross-Section	Hybrid-in-Energy	Surface Tracking
Total Runtime (min)	27.6	28.0	46.4
Inactive Rate (p/s)	731539	721782	520551
Active Rate (p/s)	676290	666265	387293
k-eff (combined estimator)	1.00020 +/- 0.00003	1.00018 +/- 0.00003	1.00014 +/- 0.00003
Leakage Fraction	0.00003 +/- 0.00000	0.00003 +/- 0.00000	0.00003 +/- 0.00000

ARC-100

QoI	Hybrid-in-Cross-Section	Hybrid-in-Energy	Surface Tracking
Total Runtime (min)	27.1	27.2	96.8
Inactive Rate (p/s)	738603	736850	252234
Active Rate (p/s)	679519	676248	185196
k-eff (combined estimator)	1.02695 +/- 0.00002	1.02693 +/- 0.00002	1.02698 +/- 0.00002
Leakage Fraction	0.24883 +/- 0.00002	0.24884 +/- 0.00002	0.24879 +/- 0.00001

Closes #3975

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)