Add domain decomposition for random ray solver

[Suark94]

Description

This PR adds a domain decomposition capability to the random ray solver, allowing users to run random ray calculations in parallel across multiple MPI ranks. This lets random ray problems that are too large to fit in the memory of a single computational node be spread across multiple nodes and enables greater scalability on HPC clusters.

This feature is enabled automatically when the random ray solver is run with more than one MPI rank. No additional input settings are required. The geometry is decomposed with a capacity-constrained Voronoi tessellation, where each Voronoi region corresponds to one MPI rank's subdomain. A source region is assigned to a Voronoi region/ MPI rank based on the squared distance between the source region centroid and the Voronoi region centroid, adjusted with an additive weight. Using the weighted squared distance rather than the simple Euclidean distance yields compact, centroidal subdomains and prevents elongated or scattered subdomains. During the first 5 iterations, the weights are adjusted according to the measured rank load to balance the work across subdomains.

Both CSG and CAD geometries are supported. CAD models require MOAB version 5.2.0 or later.

The full method documentation is included in this PR under docs/source/methods/random_ray.rst (see the new "Domain Decomposition" section).

Fixes #3009

Verification

The implemented domain decomposition scheme preserves the original source region mesh and ray sampling routines regardless of the number of parallel MPI processes. Simulation results are therefore fully reproducible.

For verification, the random ray solver with and without domain decomposition was tested on the 2D C5G7 benchmark problem. This is a multi-group benchmark problem with 7 energy groups, which describes four simplified 2D reactor assemblies surrounded by a reflector region.

For the verification, a simulation mesh with 142,964 flat source regions was used. The simulation was run twice with the same random number seed: once with the original version of OpenMC with a single MPI rank and 128 OpenMP threads (no domain decomposition), and once with the implemented domain decomposition scheme on 128 MPI ranks with 1 thread per rank. Both runs produced an identical eigenvalue of 1.18626 ± 0.00038. Pin powers were also compared against the reference multi-group Monte Carlo benchmark solution, with average absolute pin power errors of 0.5534%, and a maximum pin power error of 2.0456% in both cases.

C5G7 geometry	Rank subdomains (128 MPI ranks)

Performance

The performance of the domain decomposition scheme was tested with a strong scaling study on the UKAEA “Simple Tokamak” benchmark, a large fixed source fusion neutronics simulation problem.

This problem is a nearly “worst case” challenge problem for domain decomposition as it features extreme spatial mesh resolution differences, with fine detail and small cells in areas like the divertor, and coarse mesh resolution in void areas of the problem, giving 1,626,536 source regions in total. Furthermore, the high complexity of the wall of the fusion device results in some cells being bounded by hundreds of CSG surfaces (resulting in extremely high ray tracing costs), while other areas (like the shield wall) are constrained with just a few planar surfaces. All results were obtained from simulations on the Improv cluster at Argonne National Laboratory. It was found that OpenMC performed best on the Improv architecture when 8 MPI ranks (each with 16 OpenMP threads) were used on each node.

Simulations were scaled from one node out to 40 nodes (5120 cores) while the mesh size was kept fixed. The simulations were run with 4 million rays per batch.

Runtime	Times by code component

The domain decomposition feature enables a continuous reduction in runtime all the way out to 40 nodes. As the problem is scaled to higher node counts, an increasing deviation from ideal scaling behavior can be observed due to the progressive loss of parallel efficiency. At the maximum of 40 nodes, only 38 % parallel efficiency was achieved compared to single node operation. At this scale, the cost of the domain decomposition operations begin to dominate the runtime of the code.

Summary of changes

New classes / files

• DecompositionMap (decomposition_map): Contains the subdomain_map_, which maps source regions to MPI ranks. Provides methods for initializing the Voronoi volumes, estimating and balancing load, exchanging source regions and determining the owner rank of newly discovered source regions.
• RayBank (ray_bank): Contains the list of rays to be sent to their new owner ranks and the methods for exchanging ray data between MPI ranks.

Random ray transport (random_ray)

• Added two structs: RayBufferContainer, which holds the scalar members and vector fields of a ray while it is stored in the buffer, and RayExchangeData, which collects the scalar members for communication between MPI ranks.
• Rays are now stopped once they leave their current MPI rank's subdomain, and their state is buffered in a RayBufferContainer. restart_ray() reinitializes a ray that has been received by its new owner rank. Additional methods check whether a ray has left its subdomain and pack it into the buffer.
• The number of ray-trace operations is now counted for use in load balancing.

Random ray simulation (random_ray_simulation)

• Added a separate transport_sweep_decomp() that buffers rays as they leave a subdomain
• Additional domain-decomposition timers and output data are printed in the simulation summary.

Source regions (source_region)

• Added ScalarSourceRegionFields, which holds all scalar fields of a source region.
• Added SourceRegion::merge(), which merges two source regions that were discovered and contested by separate ranks in the same transport sweep.
• Flat source regions now also store centroids, which are used for MPI rank assignment during load optimization.

Flat source domain (flat_source_domain)

• Added is_geometry_3D() to probe the geometry type
• Added output_to_vtk_decomp() to compose the VTK output data from the separate MPI processes.
• External source strength and tally volumes are reduced across MPI ranks.

Supporting changes

• parallel_map: added size() and erase(). operator[] now raises an error when the requested key is not present to help with debugging.
• cell: Cell can now return its number of surfaces.
• constants: added MAX_N_HANDLES (max DAGMC entity handles sent when exchanging rays, set to 5) and ITER_LOAD_BALANCE (max batches over which load balancing is performed, set to 5).
• source: satisfies_spatial_constraints moved from protected to public so it can be used to confine Voronoi regions to the real geometry at the start of a simulation.
• New timers quantify the cost of the domain-decomposition-specific operations.

Testing

No new regression tests were added, as the existing test suite already covers this feature. Domain decomposition activates automatically whenever the random ray solver runs on more than one MPI rank, so every existing random_ray_* regression test doubles as a domain decomposition test when the test suite is run with the --mpi flag. Since the reference results for these tests are generated from serial runs, the MPI runs directly verify that the decomposed solver reproduces the single-rank solution for the eigenvalue, fixed source, linear source, adjoint and mesh-tally variants.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) 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)

[GuySten]

LGTM.
The only thing left to do is to rerun the benchmark after the fixes I introduced.

[GuySten]

I've checked that the results are identical when using mpi ranks 16,32,64 on a C5G7-like model.

I am planning to merge in a few days if no one objects.