Delta tracking

[nuclearkevin]

Description

This pull request implements standard delta tracking in OpenMC for the continuous-energy Monte Carlo solver. This includes history-based mode, event-based mode, and photon transport. Both the history-based and event-based implementations support the shared secondary bank recently implemented in #3863.

Neutron majorants are generated in a two-step process. The energy grids for all nuclides (including both smooth cross sections and unresolved resonance probability tables) in the problem are unionized. Afterwards, the total cross section for each nuclide is interpolated to the unionized grid to facilitate the computation of a per-material maximum cross section (the maximum value of which is used for the global majorant cross section). Thermal scattering laws are applied for materials which specify them. Unresolved resonance probability tables are treated in a conservative manner, where the maximum value over all bands is computed and then interpolated onto the unionized grid. Temperatures are handled in a conservative manner, where the majorant is the maximum over the entire temperature grid. Distributed cell densities are applied by finding the maximum density multiplier over the domain for each material, caching the result, and using the corresponding maximum density multiplier to adjust the per-material maximum cross section. Photon majorants are generated in a similar manner, but their computation is more straight-forward. The smooth cross section grids of all required photon interactions are unionized, and cross sections are then computed on the new unionized grid to compute the majorant.

The delta tracking history-based algorithm itself is quite simple. Particles are advanced by the smallest of either the distance to next collision (sampled from the majorant cross section) or the distance to the nearest external boundary surface. At that point, the particles are located by a call to global point containment. If the particle is crossing a boundary surface, the BC is applied and the particle goes back to advancing (if its alive). Otherwise, cross sections at the real/virtual collision point are computed and used to rejection sample the majorant. If the collision is virtual, the particle goes back to advancing. Otherwise, the collision is processed.

The event-based algorithm essentially implements this process verbatim with the current event queues and different logic to shuffle particles between queues. This is great in that it yields identical particle histories between approaches (which can be seen in the testing results), but is almost certainly leaving performance on the table. The event-based implementation is provided for feature parity only, and we leave performance improvements in this part of delta tracking to future work.

Limitations

This PR does not support multigroup mode (under development by folks at LLNL, see #3902). Additionally, windowed multi-pole cross sections are not currently supported, nor are NCrystal materials. The usual delta tracking limitation of collision estimators is also enforced in this PR.

Testing

Regression tests are provided for history-based delta tracking, event-based delta tracking, distributed cell densities, and all boundary conditions. The first three tests focus on the majorant calculation (to ensure it hits features that impact cross section generation), and the last test ensures boundary conditions function as intended. I intended to add a temperature interpolation test, but none of the nuclides in the NNDC library appear to have temperature interpolation tables. I was also considering adding a test for distributed materials and the shared secondary bank algorithm, but decided against it for the initial PR due to the large size. I am more than happy to add more tests if people are fine with reviewing them. 😃

Verification & Initial Performance Results

Initial verification and performance results can be found in this delta tracking benchmark repo, alongside the models used. The following systems are modelled:

• Jezebel;
• A LWR pincell;
• A SFR pincell;
• A HTGR compact (15% packing fraction).

All simulations are run in history-based mode in both neutron-only and neutron+photon coupled transport. The cases use 1000, 10000, or 100000 particles per batch with 100 inactive batches and 1100 total batches. Several quantities are tallied. These include neutron & photon spectra (101 energy bins), total reaction rates and particle fluxes (structured mesh tallies with 51x51 bins in the x-y plane) that filter by particle type, and global quantities (leakage fraction and $k_{eff}$ with the collision estimator). It should be noted that the LWR, SFR, and Jezebel problems do not contain surface subdivisions (the best-case scenario for surface tracking). All problems are run on a single node of Improv @ ANL (128 available OpenMP threads).

Integral Quantities

We start by comparing the convergence of integral quantities. Surface tracking and delta tracking agree within $3\sigma$ (the error bounds in the plots) for both leakage fractions and $k_{eff}$, which is exactly what we want to see! This gives confidence that the implementation of both reflective and vacuum boundary conditions is correct, and the majorant is conservative. We do not hit the rejection sampling guard for any of these problems even after simulating ~ $10^8$ particle histories.

Spectra

The neutron spectra (computed with single-physics runs to improve statistics) and photon spectra can be found below for the 100000 particle/batch case. The agreement for all cases is within $3\sigma$.

In addition to the raw spectra, the $L_2$ difference between the tracking methods is plotted below. We'd expect to see that the difference follows a $1/\sqrt{N}$ trend with respect to the number of particles per batch, which would indicate that delta tracking and surface tracking are converging to the same values. We're close to the expected $1/\sqrt{N}$ convergence rate, with some noise due to the low number of samples in the ultra thermal range.

Spatial Distributions

It is difficult to visualize differences between the spatial distributions, and so only the $L_{2}$ differences in the total reaction rates are plotted. Once again, the expected $1/\sqrt{N}$ rate of convergence is observed for both the neutron and photon reaction rates!

Performance

We compare the performance of delta tracking to surface tracking by defining the following speedup: $S = R_{delta}/R_{surface}$. $R_{delta}$ is the active tracking rate with delta tracking, and $R_{surface}$ is the active tracking rate with surface tracking. We see modest speedups when running delta tracking in Jezebel and the pincell problems (SFR & LWR). The expected performance improvement materializes when running the HTGR compact problem, where delta tracking is close to 5x faster than surface tracking!

Future Work

This is only an initial implementation of delta tracking. Next steps will include the implementation of hybrid tracking (in a manner similar to Serpent) to mitigate performance degradations that are commonly observed when using delta tracking in problems with localized absorbers. I'm also looking at improving cell containment routines using novel data structures. Performance benchmarking is ongoing with production models (whole core problems). I am also quite interested in optimizing the event-based delta tracking approach to pave the way for a GPU implementation of the method.

Closes #24, closes #3342, closes #3343

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)

[sallustius]

Looping in @realmisch on this PR for awareness.

[nuclearkevin]

Additional Production Testing

Many thanks to the folks who provided models for testing! All models are run with 128 threads of a single node of Improv @ ANL.

Gas-Cooled Microreactor (model provided by Lewis Gross at UW-Madison)

This model 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. These absorbers inflate the majorant in the thermal portion of the spectrum, resulting in reduced delta tracking performance. We run the model with 10000 particles per batch, 100 inactive batches, and 1000 active batches. I believe the reduced inactive rate (compared to the active rate) is caused by the building of cell neighbourhood lists in inactive batches, which is substantially slower than tallying in this problem due to the TRISO/poison particle surfaces.

QoI	Delta Tracking	Surface Tracking
Inactive Rate	375.992	567.406
Active Rate	427.563	595.406
Inactive Speedup	~0.66x	N/A
Active Speedup	~0.71x	N/A
k-eff (collision estimator)	1.01790 +/- 0.00035	1.01773 +/- 0.00034
Leakage Fraction	0.08870 +/- 0.00010	0.08863 +/- 0.00010

BEAVRS

This model 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. In an identical manner to the GCMR, these absorbers inflate the majorant in the thermal region (just not to the same degree as the poison particles). We run the model with 1000000 particles per batch, 100 inactive batches, and 1000 active batches.

QoI	Delta Tracking	Surface Tracking
Inactive Rate	238087	548034
Active Rate	232218	510892
Inactive Speedup	~0.43x	N/A
Active Speedup	~0.45x	N/A
k-eff (collision estimator)	1.00021 +/- 0.00003	1.00016 +/- 0.00004
Leakage Fraction	0.00003 +/- 0.00000	0.00003 +/- 0.00000

ARC-100 (model provided by Nathan Glaser at OSU)

This model 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 vast majority of neutrons in this SFR transport above energies where the inflated majorant degrades performance, and so we see a fairly significant speedup when running delta tracking. We run the model with 1000000 particles per batch, 100 inactive batches, and 1000 active batches.

QoI	Delta Tracking	Surface Tracking
Inactive Rate	736634	248838
Active Rate	678635	237658
Inactive Speedup	~2.96x	N/A
Active Speedup	~2.85x	N/A
k-eff (collision estimator)	1.02694 +/- 0.00003	1.02697 +/- 0.00003
Leakage Fraction	0.24884 +/- 0.00002	0.24879 +/- 0.00001

Takeaways

Delta tracking (as implemented in this PR) can improve performance for certain whole-core systems which are not susceptible to inflated majorant cross sections, namely fast reactors. We anticipate that hybrid delta tracking (as implemented in Serpent) will mitigate the observed performance degradations in thermal systems. A first cut of hybrid tracking in OpenMC based on this PR shows promise, where we see decent speedups for thermal problems with strong absorbers (results to come in the hybrid tracking PR). 😄

[pshriwise]

Great stuff @nuclearkevin!! Really looking forward to finally getting this feature in. Here's a first round of review to consider. I'll start working through the tests in the next pass.

[nuclearkevin]

The most recent set of commits (a913f1a -> 7f497ff) fix a bug I found when attempting to combined the collision and absorption $k_{eff}$ estimates. The issue I was running into was standard deviations of the two-component combined estimate were several orders of magnitude smaller than the individual estimators. We expect the combined estimator to have a lower standard deviation as its merging two correlated estimates of the mean, but multiple orders of magnitude is non-physical. To demonstrate the problem, I varied the number of active batches in both an LWR and SFR simulation with delta tracking and plotted the standard deviations of the different estimates. Blue and red are the collision and absorption estimators. Purple is the old two-component combined estimator, which is clearly too small and doesn't follow the expected convergence rate predicted by the central limit theorem. Green is the fixed two-component combined estimator, which follows the central limit theorem.

The fix implements the two-component estimator as its written in Urbatsch's original report (LA-12658-MS, Eq. 36 and Eq. 40). The main issue was the use of the matrix $\mathbf{\Sigma}$ instead of the matrix $\mathbf{S}$ (where $\mathbf{\Sigma} = \mathbf{S} / (n-1)$ ), which works for the mean value of $k_{eff}$ as the factors of $n-1$ cancel out. This doesn't work when estimating the combined standard deviation.

The two-component combined estimator is used by a fallback in surface tracking to catch edge cases where the three-component combined estimator becomes numerically unstable, and doesn't see much use outside of the pathological case (multi-group Monte Carlo with survival biasing). This explains why it took so long to find this bug.

I had to regold the multi-group survival biasing test and all of the delta tracking tests to account for the fix. Only the standard deviation required a regold, the mean values are unchanged.

[pshriwise]

Another round of review for your consideration @nuclearkevin.

Some higher level questions:

• Thank you for a very comprehensive set of tests to go with this new feature! The test directory layout, particularly for photon transport using True/False aren't very clear at first glance. If we can flatten this layout a bit or use more descriptive subdirectory names I think it would help.
• This is relatively new ground, so kind of an open question: The inputs for the regression models all seem to be based on the newly added example of a 2x2 PWR lattice. What do you of reducing that repitition by importing and modifying the example model as needed in the test files?
• It looks like there are places where we might be able to use the openmc::vector instead of the std:: types. We're trying to stick to these types when possible to ease a future integration of GPU code. Can you sub these types in where appropriate?
• There are some features we may also need to put some safeguards/warnings in for like NCrystal and surface tallies.
• Have we tested this for DAGMC geometries as well?

[pshriwise]

We'll want to include this parameter in the "settings.xml" documentation.

[pshriwise]

Let's keep the docstrings a litle higher level, leaving out the details on the std:: functions being called withing. If the implementation or work done within this method ever changes they will become stale and we're likley to miss that.

[pshriwise]

I'm a little confused as to why we pass the grid in when it's already a member of the class here.

[pshriwise]

Suggested change

	//!< A vector of materials contained in maj_universe_.
	//!< A vector of materials contained in the majorant's universe.

Can you also specify whether these are material IDs or indices?

[pshriwise]

Naming suggestion: calculate_max_smooth_micro_xs.

It's called out appropriately in the docstrings, but best this provides the best clarity for where it's called in the code I'd say.

[pshriwise]

Let's place this in majorant.h and rename to MAJORANT_URR_TOL.

[nuclearkevin]

I've also gone through the other static member constexpr and renamed them for consistency (since they're effectively global variables).

[pshriwise]

Can you clarify why smooth_xs is included here?

[nuclearkevin]

Ah, this shouldn't be there as we take the max with the smooth cross section again outside of the function. I've removed it

[pshriwise]

I don't think we're worried about negative XS here, but maybe > 0.0 is a little safer?

[nuclearkevin]

From my understanding, the != 0 exists as photon cross sections get incredibly large at low energies and therefore have the potential to overflow (which is caught and clamped to zero on initialization). This guard is also in photon.h (which I copied to the majorant calculation) to catch those zeros. Setting it to > will be fine, though it is different from the true photon interaction behaviour.

[nuclearkevin]

Ah, I just remembered that these cross sections are stored and interpolated in log space. Negatives therefore should be allowed. This must be a != 0.

[pshriwise]

Let's add the include for unordered_set here.

[pshriwise]

If I'm understanding correctly, it looks like an explicit setting of the estimator to TRACKLENGTH will still be accepted below. We should probably add one more after the XML value has been parsed.

[nuclearkevin]

Explicitly setting the estimator to tracklength will not work, as tracklength estimators are not allowed to override collision or analog default estimators. We set the estimator to collision before this block of code checks for the user-set estimator:

else if (est == "tracklength" || est == "track-length" ||
           est == "pathlength" || est == "path-length") {
  // If the estimator was set to an analog estimator, this means the
  // tally needs post-collision information
  if (estimator_ == TallyEstimator::ANALOG ||
      estimator_ == TallyEstimator::COLLISION) {
    throw std::runtime_error {fmt::format("Cannot use track-length "
                                          "estimator for tally {}",
      id_)};
  }

  // Set estimator to track-length estimator
  estimator_ = TallyEstimator::TRACKLENGTH;

}

I can make the error message more clear for delta tracking specifically, but that would run counter to how these estimators are restricted in OpenMC.

[nuclearkevin]

Thanks for the review @pshriwise! I've addressed all of the code comments, here's my responses to the high-level comments:

Thank you for a very comprehensive set of tests to go with this new feature! The test directory layout, particularly for photon transport using True/False aren't very clear at first glance. If we can flatten this layout a bit or use more descriptive subdirectory names I think it would help.

I completely agree! I've modified the directory naming so the subdirectories are photon_True / photon_False instead (I wrote these tests when I was still getting used to pytest's parameterization and didn't want to break anything).

This is relatively new ground, so kind of an open question: The inputs for the regression models all seem to be based on the newly added example of a 2x2 PWR lattice. What do you of reducing that repitition by importing and modifying the example model as needed in the test files?

I'm a little confused on what you're asking here. I believe code duplication is minimized as the individual test files already import the model and modify it in-place.

It looks like there are places where we might be able to use the openmc::vector instead of the std:: types. We're trying to stick to these types when possible to ease a future integration of GPU code. Can you sub these types in where appropriate?

Good point! I've subbed these in for the majorant calculation.

There are some features we may also need to put some safeguards/warnings in for like NCrystal and surface tallies.

I've already added a guard for NCrystal materials (its in materials.cpp as that's where we check for NCrystal). Good catch on surface tallies, I hadn't thought of that!

Have we tested this for DAGMC geometries as well?

I actually haven't tested this with DAGMC geometry. My main reason for not doing testing with DAGMC is that the geometry kernel is functionally unmodified, and so it should work (plus building DAGMC + Moab is a pain when rapidly iterating and swapping between branches). Obviously "should" isn't the same as "does", I'll run a few DAGMC models to make sure no particles tunnel.