Linear Source Random Ray

[RufusN]

Description

This PR introduces Linear Source (LS) support to both the eigenvalue and fixed source random ray solver mode. This PR implements feature request #3010.

Using LS changes the source in each source region from constant to a linear function. This allows the user to coarsen the
mesh and reduce the ray population used, allowing for significant overall computational time savings and potential improvements in accuracy.

3D LS is activated by setting the source_shape field
in the openmc.Settings.random_ray dictionary to 'linear' as:

settings.random_ray['source_shape'] = 'linear' 

If LS is used on a problem that uses reflective or unbounded z boundary conditions, with no variation, such that the problem is effectively 2D, instead use:

settings.random_ray['source_shape'] = 'linear_xy'

this zeros the z directional linear source component, reducing noise in the LS.

Complexities

To implement LS, the LinearSourceDomain class has been created which inherits from the FlatSourceDomain class. To begin we solve for the additional first order moments of source, included in the LinearSourceDomain::update_neutron_source() function, which is calculated from the first order scalar flux spatial moments. The source moments are then used to calculate the source gradient vector, this is done through the inversion of a spatial moment matrix that represents the spatial distribution of a source region.

With LS enabled, the source gradient is then used in the RandomRay::attenuate_flux_linear_source() function to calculate the angular flux alongside the first order angular flux spatial moment. The scalar flux moments are then calculated by the summation of the angular flux spatial moments and angular flux average, before being normalised in the same way as the scalar flux.

To facilitate these changes two new classes have been introduced, the MomentArray class and the MomentMatrix class. The MomentArray class allows for convenient manipulation of the scalar and source moment arrays (consisting of x, y and z components). The MomentMatrix class handles the inversion, calculation and determinant calculation of the moment matrix.

A LS random ray basic derivation has been added to the OpenMC docs, while a more detailed derivation has been documented at: https://www.ans.org/pubs/proceedings/article-55554/.

Validation

To validate the LS implementation we used the C5G7 benchmark, the 7 energy group miniature light water reactor. We analysed the error in k-eff and the maximum and average pin powers across varying ray population and source region number. Results are displayed in the table below.

The results show LS offers comparable, if not slightly improved, accuracy on the fine mesh case. On the coarse mesh, with the same ray population, LS maintains the same power errors, while flat source sees an 85% increase in average power error and 169% in maximum power error. On this case, with a LS runtime of 39 seconds, we can achieve the same accuracy 2.3x faster than flat source on the fine mesh. On the coarse mesh, with reduced population, flat source sees further increasing error, while LS shows only a 47% increase in average power error and 33% increase in maximum power error. However, these settings result in LS being 5.5x faster than flat source on the fine mesh, with only a marginal decrease in accuracy.

LS random ray on C5G7 was also tested in SCONE and the results can be found alongside the derivation.

Extrapolating from the scalar flux moments and scalar flux values for a specific region, multiple flux values can be calculated per region. This allows for a visual representation of LS, displayed below on the left, which can be compared to flat source on the same mesh, shown on the right.

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)

[jtramm]

Looks like we forgot to reduce the flux moments each iteration when doing domain replication with MPI. I'll work on a fix for this.

[jtramm]

Thanks @RufusN -- this is a fantastic additional to the random ray solver. Appears to paying even bigger dividends on 3D fixed source problems in my testing as compared with 2D k-eff problems. E.g., going from 1x1x1 cm mesh resolution with flat source to 10x10x10cm resolution on the Kobayashi dog-leg is pretty amazing (1000x fewer cells!). I've incorporated all my review comments directly as edits to the branch, so it is ready to go now I think.

[jtramm]

As a few others (e.g., @paulromano) had expressed some interest in taking a quick look before this one goes in, I will hold off a couple more days before merging with develop.

[jtramm]

One question that may come up: why not have the RandomRay class be polymorphic in the same way that the FlatSourceDomain class became, so we can avoid having the function:

void RandomRay::attenuate_flux(double distance, bool is_active)
{
  switch (source_shape_) {
  case RandomRaySourceShape::FLAT:
    attenuate_flux_flat_source(distance, is_active);
    break;
  case RandomRaySourceShape::LINEAR:
  case RandomRaySourceShape::LINEAR_XY:
    attenuate_flux_linear_source(distance, is_active);
    break;
  default:
    fatal_error("Unknown source shape for random ray transport.");
  }
}

?

The answer is that my goal for the random ray area is to make it trivially portable to GPU. The FlatSourceDomain polymorphism all occurs on the host, but RandomRay polymorphism has to occur on device, so is a problem. As only one or two switch statements had to be added to avoid polymorphism in the RandomRay class, I thought it was a minor enough sacrifice in cleanliness to keep the GPU delta to basically zero.