Random Ray Base Source Region Refactor

[jtramm]

Description

This PR lays the groundwork for an upcoming PR by Max Kraus that will introduce dynamic domain decomposition for random ray with MPI.

The main problem this PR solves has to do with how source regions are initialized. Currently, we can determine the exact number of "base" source regions analytically by analyzing the geometry. These "base" source regions are synonymous with the idea of a material filled cell instance in Monte Carlo.

In the original OpenMC implementation, these "base" source regions were used directly to index into all data needed to represent source regions. However, in PR #3333, we added the ability to subdivide source regions by overlaying a mesh so as to improve the accuracy of random ray without having to require the user to manually alter the geometry definition to make smaller cells. As we don't know apriori which mesh bins will actually intersect with a particular material filled cell instance, we no longer were able to allocate all data storage ahead of time. Instead, PR #3333 introduced the ability to dynamically generate new source regions when they were encountered.

Even though source regions can be dynamically generated, they will still pull from the "base source region" data that is computed up front. For instance, the starting scalar flux guess, the external source terms in each cell, and the overall starting source term are all computed once up front and stored in memory. This worked fine on a single rank, but would greatly complicate full decomposition of data in the context of MPI, as ranks might need to access any of the (potentially very many) base source regions when initializing local source regions they are claiming ownership of.

This PR removes any data storage for "base source regions". Source regions that are dynamically generated on-the-fly when a ray first visits them now have the required logic to figure out what their starting source terms etc should be. This has the immediate benefit of ensuring that we never have to use a lot of wasted data storing base source region data. The logic is now in place to do everything on the fly when needed. This is key for introducing MPI in the upcoming PR.

Changes

1. Removal of the base_source_regions_ object within the FlatSourceDomain class. New logic is added in a number of areas to dynamically initialize source regions when they are actually encountered by a ray.

2. To support the above item, two new map objects were added external_volumetric_source_map_ and mesh_map_ that are used to allow us to determine which external sources and meshes should be applied to which source regions. These are built once during initialization, and then used whenever a new source region is found to help populate it with the right source terms etc.

3. The scheme for figuring out what the source region index is given a particle's location is now contained in a single set of functions, rather than being repeated many times throughout the code.

4. The top level driving logic is simplified to no longer have a different object for the forward and adjoint solves. Rather, the same object is reused. This should significantly reduce the memory footprint of adjoint solves, and reduces a lot of the complexity of the top level driver (openmc_run_random_ray()) function.

5. Point sources can now stack. I.e., it is now valid to place multiple point sources that end up within the same source region. Previously, this would result in undefined behavior. This is a very useful feature for fusion simulations that often define the plasma source in terms of a large number of discrete point sources.

Testing Change

The "point source location" test needed to be updated. This was due to a flaw in the current logic of the code, where point sources are applied when the source region is discovered, but new source regions were still inheriting from the base source region scalar fluxes and sources. In effect, this meant that the first iteration when a source region is discovered will have a source that (erroneously) does not include the point source. In subsequent iterations, the point source will contribute to the source regions source, but not in the very first one. This does not really induce any significant/detectable error as this almost always happens in one of the first few inactive batches, such that subsequent inactive batches wash out that bias. However, for the small/short CI tests, this was causing a very slight ~0.01% difference in tallies.

The new PR fixes this logic implicitly as source terms are now computed on-the-fly, so are inclusive of external point sources. Thus, the change in the test results for this one are due to a small bug in the existing code that is fixed by this PR.

Performance

The PR should generally make no performance difference on problems that use the mesh overlay feature introduced in PR #3333. However, if the user has manually subdivided their source regions by adding more cells/surfaces to their code, and is not overlaying a mesh, then the new strategy in this PR may incur some added computational cost. This is because in the old version, some of the logic required to support dynamic generation of source regions (e.g., checking maps to determine where a source region's data is stored) could be skipped when no mesh overlay was being used. However, we now are always using that more expensive logic. For 2D C5G7 with a manually defined flat source mesh with ~142k FSRs, we go from a runtime of 145 seconds to 179 seconds (about a 23% increase in runtime). While this is a pretty significant performance loss, there are the following items to consider:

• Almost all of the cost stems from read operations in a std::unordered_map object. The C++ standard library implements this to minimize memory usage. In my testing, using one of the other available template implementations online (e.g., google has a good one) should mitigate most of this expense.

• It's useful from a reduced program complexity standpoint, as the overall indexing scheme only follows one set of logic, so is a little easier to follow and maintain.

• The people who have been using random ray almost entirely rely on the mesh overlay feature, meaning that maintaining code paths that optimize for manually subdivided geometries of less value to end users.

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)

[Copilot]

Pull Request Overview

This PR refactors the random ray implementation to remove base source region storage and simplify the handling of source regions, particularly in preparation for future MPI domain decomposition support.

• Removes the base_source_regions_ object and associated mesh subdivision flag, replacing with dynamic source region initialization
• Consolidates source region lookup and initialization logic into centralized functions
• Simplifies the simulation driver to reuse the same object for both forward and adjoint solves

Reviewed Changes

Copilot reviewed 11 out of 11 changed files in this pull request and generated no comments.

Show a summary per file

File	Description
tests/regression_tests/random_ray_point_source_locator/results_true.dat	Updates expected test results due to fixed point source initialization
src/settings.cpp	Removes mesh subdivision flag setting
src/random_ray/source_region.cpp	Changes source initialization and adds conditional flux reset logic
src/random_ray/random_ray_simulation.cpp	Major refactor removing separate forward/adjoint objects and simplifying driver
src/random_ray/random_ray.cpp	Removes mesh subdivision checks and consolidates source region lookup
src/random_ray/linear_source_domain.cpp	Refactors to use single source region updates instead of batch processing
src/random_ray/flat_source_domain.cpp	Major changes removing base source regions, adding lookup functions, and restructuring external source handling
include/openmc/random_ray/random_ray_simulation.h	Removes k_eff member and simplifies adjoint preparation method signature
include/openmc/random_ray/random_ray.h	Removes mesh subdivision enabled flag
include/openmc/random_ray/linear_source_domain.h	Updates method signature for single source region updates
include/openmc/random_ray/flat_source_domain.h	Major header changes adding lookup methods, removing base source regions, and restructuring data members

---

_{Tip: Customize your code reviews with copilot-instructions.md. Create the file or learn how to get started.}

[jtramm]

There was a race condition with multiple threads that affected results by a few pcm, but I fixed that. PR should be good to go now.

[jtramm]

Forgot to mention: this PR closes #3470.

[jtramm]

I've tacked on an additional minor change to this PR. Materials with macroscopic cross sections below 1e-6 are now treated as being pure void, so they follow the true void code paths. This prevents fluxes from blowing up/going negative due to issues with 1/Sigma_t when Sigma_t is extremely low. This issue was encountered by @shimwell as some fusion models will define a plasma material with density around 1e-9. The new check will make sure these materials can still be used without the results being corrupted. I also added a new test to make sure the check is working.

I had originally just wanted to make a fresh PR to add the low density check, but the logic for the check is heavily dependent on the original refactoring changes of this PR, so it seemed much easier just to add this in directly, given that it's just a few added lines of code.

[paulromano]

Really nice work, thanks a lot @jtramm!