Random Ray Source Region Mesh Subdivision (Cell-Under-Voxel Geometry) (#3333)

Co-authored-by: Paul Romano <[email protected]>

Author: jtramm

docs/source/_images/2x2_sr_mesh.png

@@ -470,6 +470,24 @@ found in the :ref:`random ray user guide <random_ray>`.
 
     *Default*: prng
 
+  :source_region_meshes:
+    Relates meshes to spatial domains for subdividing source regions with each domain.
+
+    :mesh:
+      Contains an ``id`` attribute and one or more ``<domain>`` sub-elements.
+
+      :id:
+        The unique identifier for the mesh.
+
+      :domain:
+        Each domain element has an ``id`` attribute and a ``type`` attribute.
+
+        :id:
+          The unique identifier for the domain.
+
+        :type:
+          The type of the domain. Can be ``material``, ``cell``, or ``universe``.
+
 ----------------------------------
 ``<resonance_scattering>`` Element
 ----------------------------------

docs/source/io_formats/settings.rst

@@ -319,20 +319,24 @@ Default behavior using OpenMC's native PRNG can be manually specified as::
 
 .. _subdivision_fsr:
 
-----------------------------------
-Subdivision of Flat Source Regions
-----------------------------------
-
-While the scattering and fission sources in Monte Carlo
-are treated continuously, they are assumed to be invariant (flat) within a
-MOC or random ray flat source region (FSR). This introduces bias into the
-simulation, which can be remedied by reducing the physical size of the FSR
-to dimensions below that of typical mean free paths of particles.
+-----------------------------
+Subdivision of Source Regions
+-----------------------------
 
-In OpenMC, this subdivision currently must be done manually. The level of
+While the scattering and fission sources in Monte Carlo are treated
+continuously, they are assumed to have a shape (flat or linear) within a MOC or
+random ray source region (SR). This introduces bias into the simulation that can
+be remedied by reducing the physical size of the SR to be smaller than the
+typical mean free paths of particles. While use of linear sources in OpenMC
+greatly reduces the error stemming from this approximation, subdivision is still
+typically required.
+
+In OpenMC, this subdivision can be done either manually by the user (by defining
+additional surfaces and cells in the geometry) or automatically by assigning a
+mesh to one or more cells, universes, or material types. The level of
 subdivision needed will be dependent on the fidelity the user requires. For
-typical light water reactor analysis, consider the following example subdivision
-of a two-dimensional 2x2 reflective pincell lattice:
+typical light water reactor analysis, consider the following example of manual
+subdivision of a two-dimensional 2x2 reflective pincell lattice:
 
 .. figure:: ../_images/2x2_materials.jpeg
     :class: with-border
@@ -344,9 +348,79 @@ of a two-dimensional 2x2 reflective pincell lattice:
     :class: with-border
     :width: 400
 
-    FSR decomposition for an asymmetrical 2x2 lattice (1.26 cm pitch)
+    Manual decomposition for an asymmetrical 2x2 lattice (1.26 cm pitch)
+
+Geometry cells can also be subdivided into small source regions by assigning a
+mesh to a list of domains, with each domain being of type
+:class:`openmc.Material`, :class:`openmc.Cell`, or :class:`openmc.Universe`. The
+idea of defining a source region as a combination of a base geometry cell and a
+mesh element is known as "cell-under-voxel" style geometry, although in OpenMC
+the mesh can be any kind and is not restricted to 3D regular voxels. An example
+of overlaying a simple 2D mesh over a geometry is given as::
+
+  sr_mesh = openmc.RegularMesh()
+  sr_mesh.dimension = (n, n)
+  sr_mesh.lower_left = (0.0, 0.0)
+  sr_mesh.upper_right = (x, y)
+  domain = geometry.root_universe
+  settings.random_ray['source_region_meshes'] = [(sr_mesh, [domain])]
+
+In the above example, we apply a single :math:`n \times n` uniform mesh over the
+entire domain by assigning it to the root universe of the geometry.
+Alternatively, we might want to apply a finer or coarser mesh to different
+regions of a 3D problem, for instance, as::
+
+  fuel = openmc.Material(name='UO2 fuel')
+  ...
+  water = openmc.Material(name='hot borated water')
+  ...
+  clad = openmc.Material(name='Zr cladding')
+  ...
+
+  coarse_mesh = openmc.RegularMesh()
+  coarse_mesh.dimension = (n, n, n)
+  coarse_mesh.lower_left = (0.0, 0.0, 0.0)
+  coarse_mesh.upper_right = (x, y, z)
+
+  fine_mesh = openmc.RegularMesh()
+  fine_mesh.dimension = (2*n, 2*n, 2*n)
+  fine_mesh.lower_left = (0.0, 0.0, 0.0)
+  fine_mesh.upper_right = (x, y, z)
+
+  settings.random_ray['source_region_meshes'] = [(fine_mesh, [fuel, clad]), (coarse_mesh, [water])]
+
+Note that we don't need to adjust the outer bounds of the mesh to tightly wrap
+the domain we assign the mesh to. Rather, OpenMC will dynamically generate
+source regions based on the mesh bins rays actually visit, such that no
+additional memory is wasted even if a domain only intersects a few mesh bins.
+Going back to our 2x2 lattice example, if using a mesh-based subdivision, this
+might look as below:
+
+.. figure:: ../_images/2x2_sr_mesh.png
+    :class: with-border
+    :width: 400
 
-In the future, automated subdivision of FSRs via mesh overlay may be supported.
+    20x20 overlaid "cell-under-voxel" mesh decomposition for an asymmetrical 2x2 lattice (1.26 cm pitch)
+
+Note that mesh-bashed subdivision is much easier for a user to implement but
+does have a few downsides compared to manual subdivision. Manual subdivision can
+be done with the specifics of the geometry in mind. As in the pincell example,
+it is more efficient to subdivide the fuel region into azimuthal sectors and
+radial rings as opposed to a Cartesian mesh. This is more efficient because the
+regions are a more uniform size and follow the material boundaries closer,
+resulting in the need for fewer source regions. Fewer source regions tends to
+equate to a faster computational speed and/or the need for fewer rays per batch
+to achieve good statistics. Additionally, applying a mesh often tends to create
+a few very small source regions, as shown in the above picture where corners of
+the mesh happen to intersect close to the actual fuel-moderator interface. These
+small regions are rarely visited by rays, which can result in inaccurate
+estimates of the source within those small regions and, thereby, numerical
+instability. However, OpenMC utilizes several techniques to detect these small
+source regions and mitigate instabilities that are associated with them. In
+conclusion, mesh overlay is a great way to subdivide any geometry into smaller
+source regions. It can be used while retaining stability, though typically at
+the cost of generating more source regions relative to an optimal manual
+subdivision.
 
 .. _usersguide_flux_norm:
 

docs/source/usersguide/random_ray.rst

@@ -130,9 +130,15 @@ random ray mode can be found in the :ref:`Random Ray User Guide <random_ray>`.
 
 .. warning::
     If using FW-CADIS weight window generation, ensure that the selected weight
-    window mesh does not subdivide any cells in the problem. In the future, this
-    restriction is intended to be relaxed, but for now subdivision of cells by a
-    mesh tally will result in undefined behavior.
+    window mesh does not subdivide any source regions in the problem. This can
+    be ensured by assigning the weight window tally mesh to the root universe so
+    as to create source region boundaries that conform to the mesh, as in the
+    example below.
+
+::
+    
+    root = model.geometry.root_universe    
+    settings.random_ray['source_region_meshes'] = [(ww_mesh, [root])]
 
 6. When running your multigroup random ray input deck, OpenMC will automatically
    run a forward solve followed by an adjoint solve, with a

docs/source/usersguide/variance_reduction.rst

@@ -59,6 +59,10 @@ constexpr double RADIAL_MESH_TOL {1e-10};
 // Maximum number of random samples per history
 constexpr int MAX_SAMPLE {100000};
 
+// Avg. number of hits per batch to be defined as a "small"
+// source region in the random ray solver
+constexpr double MIN_HITS_PER_BATCH {1.5};
+
 // ============================================================================
 // MATH AND PHYSICAL CONSTANTS
 

include/openmc/constants.h

@@ -4,8 +4,10 @@
 #include "openmc/constants.h"
 #include "openmc/openmp_interface.h"
 #include "openmc/position.h"
+#include "openmc/random_ray/parallel_map.h"
 #include "openmc/random_ray/source_region.h"
 #include "openmc/source.h"
+#include <unordered_map>
 #include <unordered_set>
 
 namespace openmc {
@@ -37,7 +39,6 @@ class FlatSourceDomain {
   void random_ray_tally();
   virtual void accumulate_iteration_flux();
   void output_to_vtk() const;
-  void all_reduce_replicated_source_regions();
   void convert_external_sources();
   void count_external_source_regions();
   void set_adjoint_sources(const vector<double>& forward_flux);
@@ -47,12 +48,34 @@ class FlatSourceDomain {
   void flatten_xs();
   void transpose_scattering_matrix();
   void serialize_final_fluxes(vector<double>& flux);
+  void apply_meshes();
+  void apply_mesh_to_cell_instances(int32_t i_cell, int32_t mesh_idx,
+    int target_material_id, const vector<int32_t>& instances,
+    bool is_target_void);
+  void apply_mesh_to_cell_and_children(int32_t i_cell, int32_t mesh_idx,
+    int32_t target_material_id, bool is_target_void);
+  void prepare_base_source_regions();
+  SourceRegionHandle get_subdivided_source_region_handle(
+    int64_t sr, int mesh_bin, Position r, double dist, Direction u);
+  void finalize_discovered_source_regions();
+  int64_t n_source_regions() const
+  {
+    return source_regions_.n_source_regions();
+  }
+  int64_t n_source_elements() const
+  {
+    return source_regions_.n_source_regions() * negroups_;
+  }
 
   //----------------------------------------------------------------------------
   // Static Data members
   static bool volume_normalized_flux_tallies_;
   static bool adjoint_; // If the user wants outputs based on the adjoint flux
 
+  // Static variables to store source region meshes and domains
+  static std::unordered_map<int, vector<std::pair<Source::DomainType, int>>>
+    mesh_domain_map_;
+
   //----------------------------------------------------------------------------
   // Static data members
   static RandomRayVolumeEstimator volume_estimator_;
@@ -61,7 +84,6 @@ class FlatSourceDomain {
   // Public Data members
   bool mapped_all_tallies_ {false}; // If all source regions have been visited
 
-  int64_t n_source_regions_ {0}; // Total number of source regions in the model
   int64_t n_external_source_regions_ {0}; // Total number of source regions with
                                           // non-zero external source terms
 
@@ -84,6 +106,27 @@ class FlatSourceDomain {
   // The abstract container holding all source region-specific data
   SourceRegionContainer source_regions_;
 
+  // Base source region container. When source region subdivision via mesh
+  // is in use, this container holds the original (non-subdivided) material
+  // filled cell instance source regions. These are useful as they can be
+  // initialized with external source and mesh domain information ahead of time.
+  // Then, dynamically discovered source regions can be initialized by cloning
+  // their base region.
+  SourceRegionContainer base_source_regions_;
+
+  // Parallel hash map holding all source regions discovered during
+  // a single iteration. This is a threadsafe data structure that is cleaned
+  // out after each iteration and stored in the "source_regions_" container.
+  // It is keyed with a SourceRegionKey, which combines the base source
+  // region index and the mesh bin.
+  ParallelMap<SourceRegionKey, SourceRegion, SourceRegionKey::HashFunctor>
+    discovered_source_regions_;
+
+  // Map that relates a SourceRegionKey to the index at which the source
+  // region can be found in the "source_regions_" container.
+  std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>
+    source_region_map_;
+
 protected:
   //----------------------------------------------------------------------------
   // Methods
@@ -100,9 +143,7 @@ class FlatSourceDomain {
 
   //----------------------------------------------------------------------------
   // Private data members
-  int negroups_;                  // Number of energy groups in simulation
-  int64_t n_source_elements_ {0}; // Total number of source regions in the model
-                                  // times the number of energy groups
+  int negroups_; // Number of energy groups in simulation
 
   double
     simulation_volume_; // Total physical volume of the simulation domain, as