Random Ray Kinetic Simulation Mode

[yardasol]

Description

This PR adds a kinetic simulation mode to OpenMC's random ray solver to allow simulation of kinetic transients (i.e. time-dependent transients without temperature/reactivity feedback)

Specifically, this PR adds the following:

• A new kinetic_simulation toggle setting.
• Material density transient capabilities.
• A new multi-group cross section type, CHI, to complement CHI_PROMPT and CHI_DELAYED
• A new score type, PRECURSORS
• A new example model, random_ray_pin_cell.
• Relevant content in the Theory Manual
• Cleaned up writing random ray metrics to the statepoint file (credits to @spasmann for getting most of the machinery off the ground for this in Add random ray info to statepoint file #3499 )
• New testing harnesses for testing kinetic simulations

I realize this is quite a large feature and may take a while to review. I'd appreciate any discussion on necessary theory/user guide additions to the doc pages.

Verification

I utilized two models to verify the kinetic simulation machinery

1. A pin-cell problem I compared against point kinetics
2. The C5G7-TD Exercise 3-1 Benchmark Problem

Pincell	C5G7 Detail

The pin cell model served as a quick sanity check for the reasonability of the my implementation’s solution. A fuel pin cell is easy to model and does not take a prohibitively long time to run. It is also relatively simple to model a pin cell using point-reactor kinetics, which serve as a robust check of my approach.

The more expensive C5G7 model was chosen due to

• Being an older OECD/NEA benchmark, cross-section data and model specification is openly available.
• An existing C5G7 model for OpenMC is available.
• Many time-dependent neutron transport publications utilize the C5G7-TD suite of benchmarks to do code-to-code comparison, so there are a plethora of reference results to compare against.

I utilize the same transient for both the pin cell model and the C5G7 model. The transient is based on Exercise 3-1 of the C5G7-TD benchmark. Moderator density is assumed to be at its nominal values at t=0 s . The moderator experiences a linear density decrease until it reaches a minimal value, $\omega=0.95$ for Exercise 3-1, at $t=1 s$ . The moderator then experiences a linear density increase until it returns to its nominal value at $t = 2s$ . The C5G7-TD reference results continue the simulation for another 8 seconds with the moderator at its nominal value. The moderator density changes imposed in the Exercise 3-1 transient, as well as the other three sub-cases for Exercise 3, are shown in the figure below:

Pincell Results

I compared the kinetic Random Ray implementation against PyRK, an open-source point-kinetics code. To generate the proper feedback coefficient, I ran several k-eigenvalue simulations along the range of densities in the Exercise 3-1 transient. Below are the PyRK and OpenMC results plotted against each other:

Flux	Precursors

There is a large deviation between the time-dependent Random Ray solution and the point kinetics solution. While the two are very relatively close for the first leg of the density transient, once the density begins to increase the scalar flux and precursor solutions diverge. We know that the radial flux solution for a cylindrical reactor is the zeroth-order Bessel function $J_0$, which decreases in magnitude as the radial point of interest moves away from the center of the reactor. The point kinetics flattens out this spatial dependence to a singular value, so it is expected that point kinetics overestimates both the scalar flux and the precursor concentration, as seen in the above figures. Similar statements can be made about the precursor solution.

I also wanted to compare the two methods (Time Isotropic (TI) and Source Derivative Propagation (SDP)). These two methods should be equivalent for this simple test case. TI is faster but will perform poorly near highly absorbing regions. SDP makes a higher order approximation for a greater execution cost.

Flux difference	Precursors difference

The SDP method produces a slightly lower total neutron flux at the start of the transient, on the order of $10^-3$ below the TI scalar flux , but begins to approach the TI scalar flux until the inflection point at $t=1 s$ when it becomes slightly larger than the SDP solution. At $t=2 s$, the SDP solution begins approaching the TI solution again before ending slightly below it for the remainder of the simulation. These errors are small enough that it is reasonable to conclude that the SDP and TI methods are equivalent for small models.

C5G7 results

I compared the results of the C5G7-TD transient run with OpenMC with reference data from Shen et al. (2019) on the same model and transient using the MOC code MPACT. The active length, inactive length, and rays per batch come from recommended parameters given by Tramm for the C5G7 model. The active and inactive batches are based on numbers used by Kraus et al. (2025). These simulation parameters are visible in the table below:

Item	Value
Inactive Length [cm]	628.0
Active Length [cm]	13.0
Batches	10000
Inactive Batches	3500
Rays per batch	650
Time step size [s]	0.01
Backwards Difference Order	3

The figures below compare the MPACT reference solution to the OpenMC implementation of TI and SDP. All three methods are virtually indistinguishable during the density decrease and increase. The absolute error reaches a negative maximum of -0.02 and a positive maximum of 0.02 at the end of the density ramp. During the flat part of the density transient ( t=2 s and onward), the normalized fission rates of the SDP and TI methods are \sim0.05 below the MPACT normalized fission rate. Notice the oscillations in the absolute error. This is due to the stochastic nature of time-dependent Random Ray. The shape of the error plot and oscillations are very similar to Kraus’ findings.

Fission Rate	Fission Rate Absolute Error

Run Time Analysis

Model	Time Derivative Method	Avg. CPU Time [s]
Pin Cell	TI	21.268 $\pm$ 0.246
Pin Cell	SDP	22.350 $\pm$ 0.234
C5G7	TI	830.348 $\pm$ 2.786
C5G7	SDP	1020.249 $\pm$ 2.759

While the TI and SDP methods are equvialent, the TI method is around $\sim 5$% and $\sim 19$% faster than the SDP method for the pin cell and C5G7 Models.

Issues/Planned improvements

• Construction of the tally_delay_tasks object in FlatSourceDomain::convert_source_regions_to_tallies copies a lot of existing code to construct tally_tasks. This shared code should ideally be parameterized in a helper function. An even better solution would be generalization of tally_tasks to include delay groups
• There is some missing functionality, namely treatment of explicit void regions in the time-dependent case, support for time-dependent external sources, and continuing a finished simulation from a final timestep. This functionality will be added in a PR coming in the next several months.
• The openmc_run_random_ray() function is getting quite long and has a lot of repeated code. I have a plan to parameterize this but this is for a future PR.

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)