Optionally compute bounding boxes in Mesh.material_volumes (#3731)

Author: paulromano

include/openmc/capi.h

@@ -116,7 +116,7 @@ int openmc_mesh_set_id(int32_t index, int32_t id);
 int openmc_mesh_get_n_elements(int32_t index, size_t* n);
 int openmc_mesh_get_volumes(int32_t index, double* volumes);
 int openmc_mesh_material_volumes(int32_t index, int nx, int ny, int nz,
-  int max_mats, int32_t* materials, double* volumes);
+  int max_mats, int32_t* materials, double* volumes, double* bboxes);
 int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh);
 int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh);
 int openmc_new_filter(const char* type, int32_t* index);

include/openmc/mesh.h

@@ -88,17 +88,24 @@ namespace detail {
 
 class MaterialVolumes {
 public:
+  MaterialVolumes(int32_t* mats, double* vols, double* bboxes, int table_size)
+    : materials_(mats), volumes_(vols), bboxes_(bboxes), table_size_(table_size)
+  {}
+
   MaterialVolumes(int32_t* mats, double* vols, int table_size)
-    : materials_(mats), volumes_(vols), table_size_(table_size)
+    : MaterialVolumes(mats, vols, nullptr, table_size)
   {}
 
   //! Add volume for a given material in a mesh element
   //
   //! \param[in] index_elem Index of the mesh element
   //! \param[in] index_material Index of the material within the model
   //! \param[in] volume Volume to add
-  void add_volume(int index_elem, int index_material, double volume);
-  void add_volume_unsafe(int index_elem, int index_material, double volume);
+  //! \param[in] bbox Bounding box to union into the result (optional)
+  void add_volume(int index_elem, int index_material, double volume,
+    const BoundingBox* bbox = nullptr);
+  void add_volume_unsafe(int index_elem, int index_material, double volume,
+    const BoundingBox* bbox = nullptr);
 
   // Accessors
   int32_t& materials(int i, int j) { return materials_[i * table_size_ + j]; }
@@ -113,11 +120,23 @@ class MaterialVolumes {
     return volumes_[i * table_size_ + j];
   }
 
+  double& bboxes(int i, int j, int k)
+  {
+    return bboxes_[(i * table_size_ + j) * 6 + k];
+  }
+  const double& bboxes(int i, int j, int k) const
+  {
+    return bboxes_[(i * table_size_ + j) * 6 + k];
+  }
+
+  bool has_bboxes() const { return bboxes_ != nullptr; }
+
   bool table_full() const { return table_full_; }
 
 private:
   int32_t* materials_;      //!< material index (bins, table_size)
   double* volumes_;         //!< volume in [cm^3] (bins, table_size)
+  double* bboxes_;          //!< bounding boxes (bins, table_size, 6)
   int table_size_;          //!< Size of hash table for each mesh element
   bool table_full_ {false}; //!< Whether the hash table is full
 };
@@ -240,6 +259,19 @@ class Mesh {
   void material_volumes(int nx, int ny, int nz, int max_materials,
     int32_t* materials, double* volumes) const;
 
+  //! Determine volume and bounding boxes of materials within each mesh element
+  //
+  //! \param[in] nx Number of samples in x direction
+  //! \param[in] ny Number of samples in y direction
+  //! \param[in] nz Number of samples in z direction
+  //! \param[in] max_materials Maximum number of materials in a single mesh
+  //!                          element
+  //! \param[inout] materials Array storing material indices
+  //! \param[inout] volumes Array storing volumes
+  //! \param[inout] bboxes Array storing bounding boxes (n_elems, table_size, 6)
+  void material_volumes(int nx, int ny, int nz, int max_materials,
+    int32_t* materials, double* volumes, double* bboxes) const;
+
   //! Determine bounding box of mesh
   //
   //! \return Bounding box of mesh

openmc/deplete/r2s.py

@@ -269,6 +269,7 @@ def step1_neutron_transport(
             # Compute material volume fractions on the mesh
             if mat_vol_kwargs is None:
                 mat_vol_kwargs = {}
+            mat_vol_kwargs.setdefault('bounding_boxes', True)
             self.results['mesh_material_volumes'] = mmv = comm.bcast(
                 self.domains.material_volumes(self.neutron_model, **mat_vol_kwargs))
 
@@ -539,14 +540,15 @@ def step3_photon_transport(
     def get_decay_photon_source_mesh(
         self,
         time_index: int = -1
-    ) -> list[openmc.MeshSource]:
+    ) -> list[openmc.IndependentSource]:
         """Create decay photon source for a mesh-based calculation.
 
-        This function creates N :class:`MeshSource` objects where N is the
-        maximum number of unique materials that appears in a single mesh
-        element. For each mesh element-material combination, and
-        IndependentSource instance is created with a spatial constraint limited
-        the sampled decay photons to the correct region.
+        For each mesh element-material combination, an
+        :class:`~openmc.IndependentSource` is created with a
+        :class:`~openmc.stats.Box` spatial distribution based on the bounding
+        box of the material within the mesh element. A material constraint is
+        also applied so that sampled source sites are limited to the correct
+        region.
 
         When the photon transport model is different from the neutron model, the
         photon MeshMaterialVolumes is used to determine whether an (element,
@@ -559,19 +561,15 @@ def get_decay_photon_source_mesh(
 
         Returns
         -------
-        list of openmc.MeshSource
-            A list of MeshSource objects, each containing IndependentSource
-            instances for the decay photons in the corresponding mesh element.
+        list of openmc.IndependentSource
+            A list of IndependentSource objects for the decay photons, one for
+            each mesh element-material combination with non-zero source strength.
 
         """
         mat_dict = self.neutron_model._get_all_materials()
 
-        # Some MeshSource objects will have empty positions; create a "null source"
-        # that is used for this case
-        null_source = openmc.IndependentSource(particle='photon', strength=0.0)
-
-        # List to hold sources for each MeshSource (length = N)
-        source_lists = []
+        # List to hold all sources
+        sources = []
 
         # Index in the overall list of activated materials
         index_mat = 0
@@ -594,7 +592,7 @@ def get_decay_photon_source_mesh(
                     if mat_id is not None
                 }
 
-            for j, (mat_id, _) in enumerate(mat_vols.by_element(index_elem)):
+            for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True):
                 # Skip void volume
                 if mat_id is None:
                     continue
@@ -604,30 +602,27 @@ def get_decay_photon_source_mesh(
                     index_mat += 1
                     continue
 
-                # Check whether a new MeshSource object is needed
-                if j >= len(source_lists):
-                    source_lists.append([null_source]*n_elements)
-
                 # Get activated material composition
                 original_mat = materials[index_mat]
                 activated_mat = results[time_index].get_material(str(original_mat.id))
 
-                # Create decay photon source source
+                # Create decay photon source
                 energy = activated_mat.get_decay_photon_energy()
                 if energy is not None:
                     strength = energy.integral()
-                    source_lists[j][index_elem] = openmc.IndependentSource(
+                    space = openmc.stats.Box(*bbox)
+                    sources.append(openmc.IndependentSource(
+                        space=space,
                         energy=energy,
                         particle='photon',
                         strength=strength,
                         constraints={'domains': [mat_dict[mat_id]]}
-                    )
+                    ))
 
                 # Increment index of activated material
                 index_mat += 1
 
-        # Return list of mesh sources
-        return [openmc.MeshSource(self.domains, sources) for sources in source_lists]
+        return sources
 
     def load_results(self, path: PathLike):
         """Load results from a previous R2S calculation.

openmc/lib/mesh.py

@@ -1,5 +1,5 @@
 from collections.abc import Mapping, Sequence
-from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER,
+from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER, c_void_p,
                     create_string_buffer, c_size_t)
 from math import sqrt
 import sys
@@ -47,7 +47,8 @@
 _dll.openmc_mesh_bounding_box.restype = c_int
 _dll.openmc_mesh_bounding_box.errcheck = _error_handler
 _dll.openmc_mesh_material_volumes.argtypes = [
-    c_int32, c_int, c_int, c_int, c_int, arr_2d_int32, arr_2d_double]
+    c_int32, c_int, c_int, c_int, c_int, arr_2d_int32, arr_2d_double,
+    c_void_p]
 _dll.openmc_mesh_material_volumes.restype = c_int
 _dll.openmc_mesh_material_volumes.errcheck = _error_handler
 _dll.openmc_mesh_get_plot_bins.argtypes = [
@@ -188,6 +189,7 @@ def material_volumes(
             n_samples: int | tuple[int, int, int] = 10_000,
             max_materials: int = 4,
             output: bool = True,
+            bounding_boxes: bool = False,
     ) -> MeshMaterialVolumes:
         """Determine volume of materials in each mesh element.
 
@@ -213,6 +215,11 @@ def material_volumes(
             Estimated maximum number of materials in any given mesh element.
         output : bool, optional
             Whether or not to show output.
+        bounding_boxes : bool, optional
+            Whether or not to compute an axis-aligned bounding box for each
+            (mesh element, material) combination. When enabled, the bounding
+            box encloses the ray-estimator prisms used for the volume
+            estimation.
 
         Returns
         -------
@@ -243,23 +250,36 @@ def material_volumes(
         table_size = slot_factor*max_materials
         materials = np.full((n, table_size), EMPTY_SLOT, dtype=np.int32)
         volumes = np.zeros((n, table_size), dtype=np.float64)
+        bboxes = None
+        if bounding_boxes:
+            bboxes = np.empty((n, table_size, 6), dtype=np.float64)
+            bboxes[..., 0:3] = np.inf
+            bboxes[..., 3:6] = -np.inf
 
         # Run material volume calculation
         while True:
             try:
+                bboxes_ptr = None
+                if bboxes is not None:
+                    bboxes_ptr = bboxes.ctypes.data_as(POINTER(c_double))
                 with quiet_dll(output):
                     _dll.openmc_mesh_material_volumes(
-                        self._index, nx, ny, nz, table_size, materials, volumes)
+                        self._index, nx, ny, nz, table_size, materials,
+                        volumes, bboxes_ptr)
             except AllocationError:
                 # Increase size of result array and try again
                 table_size *= 2
                 materials = np.full((n, table_size), EMPTY_SLOT, dtype=np.int32)
                 volumes = np.zeros((n, table_size), dtype=np.float64)
+                if bounding_boxes:
+                    bboxes = np.empty((n, table_size, 6), dtype=np.float64)
+                    bboxes[..., 0:3] = np.inf
+                    bboxes[..., 3:6] = -np.inf
             else:
                 # If no error, break out of loop
                 break
 
-        return MeshMaterialVolumes(materials, volumes)
+        return MeshMaterialVolumes(materials, volumes, bboxes)
 
     def get_plot_bins(
             self,