Random Ray Source Region Refactor

[jtramm]

Description

This PR changes the way source region data is stored in OpenMC. It is generally a refactor -- it does not change any program behavior or output, save for one bug fix (as described later on below).

Previously, source region data was stored in "Structure of Arrays" (SoA) format within FlatSourceDomain or LinearSourceDomain objects, with a serialized array constituting each variable (with values for each source region and energy group within each array). This is the most efficient storage format for this use case, and allows for easy contiguous looping over all elements of each variable type.

This has worked great to date, but the downside is that this storage technique is not very friendly to the abstraction of a single source region object. Having a single object that holds only the info for a single source region, e.g., a SourceRegion class,will be useful in upcoming work. Additionally, the current SoA format is not very amenable to dynamically appending source regions, which will be desirable for upcoming features.

In this PR, a SourceRegion class is made which contains all data for a single source region (including all energy groups - making this object also a SoA). Naively we could just store this object in a vector to represent all source regions, as an Array of Structures of Arrays (AoSoA). However, single variables (like the scalar flux) would no longer be contiguous in memory, and would also involve another layer of indirection. I tried this AoSoA storage scheme out and found that it was about 30% slower on 2D C5G7 as compared to the original SoA.

Thus, to allow for both abstraction and efficiency, I have introduced a second new class SourceRegionContainer that essentially retains the original SoA format, but emits accessors that abstract the access to the underlying data. Additionally, it can handle appending SourceRegion objects dynamically, with the contents distributed automatically to the various vectors. It can also handle reduction of itself.

Another question is why we didn't just leave the SoA fields in the FlatSourceDomain and LinearSourceDomain classes and add methods to them for appending the new SourceRegion object. This would have worked, but I like the new hierarchy as it allows for the SourceRegionContainer class to really just store source region-specific information only, and only have methods that have to do with the storage of that data. The higher level FlatSourceDomain classes contain additional data that is not specific to each source region (e.g., MGXS data, etc), as well as the physics methods for updating fluxes between iterations etc. Basically, we are now separating out the storage vs physics/simulation functionality into two different classes, which I think will be a useful abstraction.

One might also wonder why there is only one type of SourceRegion instead of having flat and linear versions. It seems wasteful to include the extra fields when just doing a flat source simulation, however, when the objects are added into the SourceRegionContainer, the SoA arrays for the linear terms will not be allocated if doing flat source, so there will be no wasted memory.

Motivation

Specifically, this PR is in support of the upcoming "cell-under-voxel" feature for random ray (#2832) that will allow the overlaying of an arbitrary mesh geometry over the simulation so as to allow for automated subdivision of source regions. As the total number of physical source regions will not be known a priori when doing automated on-the-fly source region subdivision, the algorithm becomes highly dynamic. SourceRegion objects will be created when new physical source regions are discovered, and then travel through a variety of shared memory data structures so as to minimize locking. Eventually, they will get added to the new SourceRegionContainer class for permanent residence. Base material filled cell instances will also be created (as a second SourceRegionContainer instance) where external source data will reside.

OpenMPMutex Bug Fix

Each source region has a OpenMPMutex object that is used to control memory access during transport. Due to the way I initialize the new data structure, the copy constructor of OpenMPMutex is now being called when it wasn't before, resulting in the discovery of an error in this class. I corrected the copy constructor so that it is working as expected now, and the underlying openmp locks are being handled correctly. I thought about putting this in a separate PR, but it is good to include here to avoid having to do backflips to make this PR work independent of the mutex bug fix.

In the old implementation, the copy constructor was (attempting to) call the regular class constructor, though this was actually just initializing a new object and immediately destroying it, resulting in the new object having a lock that is uninitialized. The new version appropriately initializes the lock in the copy constructor. Notably, it doesn't have to copy anything from the object it is copying, as locks need to be unique and should not actually be copied, so we just make a new one when the copy constructor is being called. The old assignment operator was working, but had an unneeded object creation/destruction, so that has been removed.

Performance Impacts

As the underlying contiguous storage representation is basically unchanged, I have not observed any performance degradation. In my testing, it was actually about 2% faster (though this could just be noise). There is no additional memory usage from the new method. In contrast, the simpler AoSoA storage that I tried would have increased memory usage and had about a 30% performance penalty.

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)