Fixed Source Random Ray

[jtramm]

Description

This PR introduces fixed source support to the random ray solver mode, which currently only supports eigenvalue solves. This PR implements feature request #2829 as part of the Random Ray FW-CADIS project.

Unlike in MC where fixed source involves a change in how particles themselves are sampled, the sampling and movement of rays throughout the geometry in random ray mode is identical in fixed source vs. eigenvalue modes. Instead, the fixed source term is simply added on to the overall source term in each source region when the source is being updated each iteration.

From a user's perspective, definition of sources and usage of fixed source mode is mostly the same as with the Monte Carlo solver. There are a few exceptions:

• While there are a variety of source types available in MC, this PR only supports volumetric, isotropic, discrete energy source terms for random ray. This type of source is the easiest to support in random ray, as the fixed volumetric term can simply be added to the multigroup scattering source term that is already being calculated for each source region (cell) in the problem.
• The same interface is used for defining fixed sources, but the domain argument is used to specify the materials, cells, or universes that the source should be applied to. OpenMC will then generate a fixed source term for each flat source region that is within the domain specified. Cells and universes can be specified at any level of the geometry hierarchy and the appropriate source regions will be found.
• In fixed source mode, both the ray source as well as a fixed particle source are required. The ray source is also required for eigenvalue mode, and has the same restrictions (namely, it must cover the full geometry and be uniform in space and angle).

A basic example of setting up a random ray source and a fixed particle source are shown below:

    # Define geometry, etc.
    ...
    ebins = [1e-5, 1e-1, 20.0e6]
    source_cell = openmc.Cell(fill=source_mat, name='cell where fixed source will be')
    ...
    # Define physical neutron fixed source
    energy_points = [1.0e-2, 1.0e1]
    strengths = [0.25, 0.75]
    energy_distribution = openmc.stats.Discrete(x=energy_points,p=strengths)
    neutron_source = openmc.IndependentSource(energy=energy_distribution, domains=[source_cell], strength=1.0)

    # Add fixed neutron source to settings
    settings.source = [neutron_source]

    # Create an initial uniform spatial source distribution for sampling rays
    pitch = 1.26
    lower_left  = (-pitch, -pitch, -pitch)
    upper_right = ( pitch,  pitch,  pitch)
    uniform_dist = openmc.stats.Box(lower_left, upper_right)

    # Add ray source to settings
    settings.random_ray['ray_source'] = openmc.IndependentSource(space=uniform_dist)

More details and examples are given in the updates to the documentation.

Complexities

Numerically, the only real change being made is the addition of a few lines to FlatSourceDomain::update_neutron_source(double k_eff) where we add in the fixed source term:

// Add fixed source source if in fixed source mode
#pragma omp parallel for
    for (int se = 0; se < n_source_elements_; se++) {
      source_[se] += fixed_source_[se];
    }
  }

However there is some complexity in how we translate a user-specified fixed source into a total source term for each source region. In my implementation, this involves several hundred lines of rather dense code that begins with the FlatSourceDomain::convert_fixed_sources() function and calls several other functions. The reason for this complexity is that it is actually not trivial to come up with a unique list of all material filled cell instances that correspond with a particular material, universe, or cell. Given that only material filled cell instances will correspond to random ray source regions, a navigation of the geometry hierarchy needs to be done to find contained material filled cell instances for higher level cells and universes.

Thankfully, there was a nicely written utility function (Cell::get_contained_cells()) that I believe @pshriwise wrote, which I was able to use to reduce the complexity and speed this process up. The basic idea here is that we loop over all sources and call this function to get a list of material filled cell instances (which correspond 1:1 with random ray FSRs and can be used to index to specific FSRs) to apply the source to. I also tested a more naive approach where for each FSR, we checked against each source to see if the FSR was contained within the domain for the source, but the inverse approach did not scale well and became really expensive for a lot of test problems.

Tally Normalization

This PR also changes the way tallies are normalized so as to get the same results for flux tallies as in MC fixed source. While reaction rate tallies in the current eigenvalue random ray solver are correct, I realized that the original solver was not normalizing flux tallies by volume correctly. The new tally normalization method corrects that issue by independently tracking the total volume of each tally region, so flux tallies should now be correct for both eigenvalue and fixed source problems. This resulted in changing the expected outputs for the existing two random ray eigenvalue python tests.

Validation

To validate the new solver mode, I utilized the standard Kobayashi dog leg problem. This is a simple 3D problem with specified 1 group cross sections, which represents a streaming duct channel surrounded by an absorbing material. A single fixed source region is located at one end of the duct channel. Tallies along three 1D lines through the geometry are specified as part of the benchmark, which we label as lines 3A, 3B, and 3C corresponding to the benchmark specification. The benchmark also gives reference solutions along points, but for making it easier to compare MGMC and random ray, we regenerate our reference solution using spatial boxes around the tally locations with MGMC to eliminate consideration of bias from space vs. point estimator tallies. A basic visualization, showing the duct is below:

And the key results compared to OpenMC's multigroup Monte Carlo mode is below:

The results show close agreement, though with some difference appearing along tally line 3C in position 4. Notably, the uncertainty for the MGMC solve at the 3C tally position 4 has a very high variance at +/- 57%, so the difference we're seeing is likely due to the high uncertainty in the MGMC solver given the low flux in this area (comparatively, the random ray solver has an uncertainty of +/- 1.15%). High MGMC variance (+/- 6.1%) is also the cause of the difference seen at Tally line 3B position 4. It would be good to generate a more tightly converged MGMC reference solution, and to also study the impacts of random ray FSR resolution, but for now this shows the solver is clearly working well.

Random Ray has also been tested out more rigorously on this problem in SCONE and is documented at: https://www.researchgate.net/publication/380099087_The_random_ray_method_applied_to_fixed_source_transport_problems

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]

Thanks @makeclean -- let me know if you run into any issues! Besides the Kobayashi dog leg, I've also been using the UKAEA octomak benchmark for testing/development. There's a few more key additions that are in the pipeline that should be super useful for fusion models (in particular #2989, #2990, #2832).

[gridley]

This is awesome!

As far as I can tell, this would put us pretty dang close to being to do automatic WW generation if we just had some simple pre-generated multigroup neutron/photon library. Regions of interest would then just have volumetric source terms with the energy spectrum coming from the dose-response energy functional.

Because the WW mesh likely is Cartesian and some approximation error is acceptable in generating a WW map, the geometry could just be voxelized and re-defined on a lattice in order to handle the FSR discretization being fine enough. I can't remember the details on moving from an adjoint solution to weight window values, but IIRC it's not hard at all once you have an approximation to the adjoint fluxes.

Am I missing anything? Maybe it doesn't make sense to discuss it here. Will do a review here in a few days!

[jtramm]

Thanks @gridley! Yes, I agree that this brings things much closer to WW capabilities.

Because the WW mesh likely is Cartesian and some approximation error is acceptable in generating a WW map, the geometry could just be voxelized and re-defined on a lattice in order to handle the FSR discretization being fine enough. I can't remember the details on moving from an adjoint solution to weight window values, but IIRC it's not hard at all once you have an approximation to the adjoint fluxes.

This is an interesting idea -- sounds like it would work. For #2832 I have so far coded up a "cell-under-voxel" type treatment, where we basically introduce a second level of ray tracing on a mesh in order to allow for virtual subdivision of FSRs. The implementation I have cooked up so far allows for the user to just define some number of meshes, and then apply them to a list of cell, material, or universe domains in the same manner as we're doing in this PR with fixed sources. The upside is that this allows for varying the subdivision resolution quite easily, for instance

• applying a really coarse mesh to void regions (where the source is 0) and a fine mesh to all other materials.
• applying a mesh to the root universe so that it gets applied everywhere.
• subdividing only specific cells and leaving others untouched

etc. The upside here is that the original geometry is preserved so there's no homogenization or anything. The downside is that extremely small FSRs can be generated, which can create a few issues that we need to deal with.

The voxel homogenization idea you're proposing sounds like it would be nice in that we wouldn't have the small FSR problem, and for just generating WW it seems like the error may be acceptable. Could also be considerably faster as there would be fewer FSRs. Might be interesting to try the idea out at some point. Feel free to add in thoughts to the discussion of #2832, or open another issue to contain discussion on it as well.

And thanks for offering to give the PR a look!

[paulromano]

Thanks for the enhancement @jtramm! A few comments for your consideration:

[paulromano]

The three tests that you've added here are almost identical other than one line, which represents a lot of redundant code. Can you please combine them into a single file? With pytest, a clean way of doing this is to have a parametrized test.

[paulromano]

@jtramm I just made an update to consolidate the three tests by using a separate directory only for the reference inputs/results. All good now!

[paulromano]

@gridley Let me know if you'd like to review before we merge.

[gridley]

Hey Paul! Yeah, I would like to review this real quick. I'll do that now.

[gridley]

I recommend checking that is is not null here. The check may happen elsewhere too, but it's good defensive programming. I wish we used assert from cassert because it'd be perfect for this kind of thing. A simple assert(is) would be great.

And it's nice that this is a no-op if compiled with -DNDEBUG.

[gridley]

Same here!

[gridley]

Alright, reasonable enough to me. It's easy to debug even if it did cause a segfault from a failed dynamic cast.

[jtramm]

Alright, reasonable enough to me. It's easy to debug even if it did cause a segfault from a failed dynamic cast.

The opportunity for a seg fault is always a valid concern!

In this specific case though, the way the random ray solver is setup I believe the validate_random_ray_inputs() will be run no matter if the C API is used or not, so there doesn't seem to be any possibility of a seg fault. The C API doesn't even allow for sources to be added, and even if that capability were added down the line, the validate_random_ray_inputs() would still be run before the random ray solver begins. If the random ray solver were to be changed in the future to expose more of it to the C API (e.g., adding next_batch() support), we'd need to add significant structural changes to the solver to re-apply fixed source terms to FSRs between batches, among other things, so would likely want to just call validate_random_ray_inputs() again after any regions where the user might have made changes.

[yardasol]

Amazing work @jtramm. Sorry I didn't get around to reviewing this. I do have a quick comment: on lines 47-49 of random_ray.rst, the following section of text was removed from the guide:

Similar to Monte Carlo,
:ref:`Shannon entropy <usersguide_entropy>` can be used to gauge whether the
combined scattering and fission source has fully developed.

I think it would be useful to keep this in somewhere, but I missed my chance to make a review! Perhaps we can add it back in somewhere during the next PR.

EDIT: Browsing through the open PRs it seems like someone is adding Shannon entropy for random ray (#3030)

[jtramm]

Thanks @yardasol! Yes, random ray Shannon entropy for random ray is currently under a separate PR (#3030) so the discussion in the docs on this topic is already slated to be expanded there.