Simulation of decay photons through the D1S method

[paulromano]

Description

This PR introduces a new capability to carry out direct 1-step (D1S) shutdown dose rate calculations, which are based on a coupled neutron–photon transport calculation where the results are scaled by "time correction factors". For those interested in more of the details, I'd recommend reading through a recent paper on the topic, such as that by Sauvan et al. The implementation here has a few pieces that come together in a relatively user-friendly workflow, consisting of functions in a new openmc.deplete.d1s module. How it works is:

1. The user tells OpenMC to use decay photons instead of prompt photons from neutron reactions by setting settings.use_decay_photons = True. Internally, this will cause OpenMC to read decay source information from a chain file and use it to construct reaction products representing decay photons.
2. A d1s.prepare_tallies function modifies tallies by adding a ParentNuclideFilter, a new filter introduced here that allows us to filter photon tallies by which radionuclide produced the photon in the first place. This is necessary for correctly applying the time correction factors, which are different for each radionuclide.
3. Time correction factors are calculation with the d1s.time_correction_factors function, which takes an irradiation schedule like any other depletion calculation (timesteps + source rates).
4. A simulation is carried out as normal. When it is done, the d1s.apply_time_correction function takes any tally with a ParentNuclideFilter and applies the time correction factors to get properly normalized results, which it returns in a new tally object.

Altogether, this looks like the following in code:

from openmc.deplete import d1s

model = openmc.Model()
…
model.settings.use_decay_photons = True

nuclides = d1s.prepare_tallies(model)
factors = d1s.time_correction_factors(nuclides, timesteps, source_rates)

output_path = model.run()

with openmc.StatePoint(output_path) as sp:
    tally = sp.get_tally(name=’dose tally’)

tally = d1s.apply_time_correction(tally, factors, time_index)

V&V

I presented some results using this implementation at TOFE earlier this year. On the FNG dose rate benchmark from SINBAD, we get the following results with D1S relative to the R2S method:

Remaining tasks

• Settle on method for sampling decay photons. There are a some tradeoffs for how different decay-photon-producing reactions are sampled because it's influenced not only by the cross section of the original activation reaction but also the decay constant of the radionuclide.
• Provide some cutoff for ignoring radionuclides with very small half-lives
• Add test for problem with analytical solution (if possible)

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 15) 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)