Random Ray Explicit Void Treatment

[jtramm]

Description

This PR introduces an explicit void estimator for the random ray solver mode. It now correctly deals with void cells defined as:

cell = openmc.Cell()

or

cell = openmc.Cell(fill=None)

Previously, the solver would give errors due to attempting to access data at the -1 index. I also add a check that forbids loading of macroscopic data with 0 values for Sigma_t, so as to encourage the user to just define their cells as above. Low sigma_t values (e.g., as in a void-like region) still work fine and are handled normally as before.

New Void Treatment

Together with @ChasingNeutrons, we found that the characteristic ODE actually simplifies significantly for void regions where Sigma_t is zero. The general form is:

dψ(s,Ω,E)/ds + Σt * ψ(s,Ω,E) = Q(s,Ω,E)

but in the case of void, we can make several simplifications. Σt is zero so, so that term disappears. Additionally, as scattering and fission are no longer possible, Q (the source) is left as only containing external contributions. As the random ray solver already assumes that external sources are defined to be constant within any particular cell, this makes Q(s,Ω,E) decay into a constant Q, such that the ODE becomes:

dψ(s,Ω,E)/ds = Q

we can thus solve for ψ(s) as:

ψ(s) = ψ(0) + Q * s

The estimator for the flux in the region can be given as:

ϕ= Q/2 * sum(ℓ^2) / sum(ℓ) + sum(ψ(0) * ℓ) / sum(ℓ)

Where the sums are over each ray crossing in a given region, with starting angular flux ψ(0) and ray crossing length ℓ.

Notably the Q/2 * sum(ℓ^2) / sum(ℓ) term is a function of geometry only, so we average this value over the entire course of the simulation, even if using the naive volume estimator.

This PR adds the above treatments for void cells, in the form of some alternate logic in the update functions between each batch, and alternative void-specific versions of the flux attenuation function during transport.

Notes on Linear Source

If we place an external source in a void material, it is important to note that any source regions in this area will effectively behave as flat source regions (as the only source contribution is from the (assumed flat) external source). While the source is flat, there are often pretty steep scalar flux gradients in void regions with external sources, so a finer source region mesh resolution than typical may be required if tallying flux at the location of the external void source. However -- it is interesting to note that high flux error in void regions does not affect the solution in solid areas of the problem, as rays are not affected by the flux in void regions (they are only impacted by any flat external sources in voids).

Validation

To ensure that the new solver is working, I created a modified Kobayashi dog leg model that replaces the duct and source materials with true void. To validate, I compared against multigroup Monte Carlo (MGMC). The random ray tallies converge as expected with MGMC. As discussed above, it is notable that the errors in the void source region are much higher (~20%) than we would typically expect if the source region were solid at the 5cm mesh level. As the mesh resolution is refined though the error drops, and as I note above, high flux error in void regions does not affect the solution in solid areas of the problem.

Testing

A new test was added to the CI. It modifies the standard 3 region random ray cube problem to create some void materials, and tests both flat and linear sources, as new void logic in the code is added depending on source shape.

Closes #2989.

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)