Compute material volumes in mesh elements based on raytracing

docs/source/pythonapi/base.rst

@@ -150,6 +150,7 @@ Constructing Tallies
    openmc.CylindricalMesh
    openmc.SphericalMesh
    openmc.UnstructuredMesh
+   openmc.MeshMaterialVolumes
    openmc.Trigger
    openmc.TallyDerivative
    openmc.Tally

include/openmc/bounding_box.h

@@ -4,6 +4,7 @@
 #include <algorithm> // for min, max
 
 #include "openmc/constants.h"
+#include "openmc/position.h"
 
 namespace openmc {
 
@@ -54,6 +55,9 @@ struct BoundingBox {
     zmax = std::max(zmax, other.zmax);
     return *this;
   }
+
+  inline Position min() const { return {xmin, ymin, zmin}; }
+  inline Position max() const { return {xmax, ymax, zmax}; }
 };
 
 } // namespace openmc

include/openmc/capi.h

@@ -107,8 +107,8 @@ int openmc_mesh_get_id(int32_t index, int32_t* id);
 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 n_sample, int bin,
-  int result_size, void* result, int* hits, uint64_t* seed);
+int openmc_mesh_material_volumes(int32_t index, int nx, int ny, int nz,
+  int max_mats, int32_t* materials, double* volumes);
 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

@@ -69,13 +69,57 @@ extern const libMesh::Parallel::Communicator* libmesh_comm;
 } // namespace settings
 #endif
 
+//==============================================================================
+//! Helper class for keeping track of volume for each material in a mesh element
+//==============================================================================
+
+namespace detail {
+
+class MaterialVolumes {
+public:
+  MaterialVolumes(int32_t* mats, double* vols, int max_materials)
+    : materials_(mats), volumes_(vols), n_mats_(max_materials)
+  {}
+
+  //! 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
+  //! \param[in] volume Volume to add
+  void add_volume(int index_elem, int index_material, double volume);
+
+  // Accessors
+  int32_t& materials(int i, int j) { return materials_[i * n_mats_ + j]; }
+  const int32_t& materials(int i, int j) const
+  {
+    return materials_[i * n_mats_ + j];
+  }
+
+  double& volumes(int i, int j) { return volumes_[i * n_mats_ + j]; }
+  const double& volumes(int i, int j) const
+  {
+    return volumes_[i * n_mats_ + j];
+  }
+
+  bool too_many_mats() const { return too_many_mats_; }
+
+private:
+  int32_t* materials_; //!< material index (bins, max_mats)
+  double* volumes_;    //!< volume in [cm^3] (bins, max_mats)
+  int n_mats_;         //!< Maximum number of materials in a single mesh element
+  bool too_many_mats_ = false; //!< Whether the maximum number of materials has
+                               //!< been exceeded
+};
+
+} // namespace detail
+
+//==============================================================================
+//! Base mesh class
+//==============================================================================
+
 class Mesh {
 public:
   // Types, aliases
-  struct MaterialVolume {
-    int32_t material; //!< material index
-    double volume;    //!< volume in [cm^3]
-  };
 
   // Constructors and destructor
   Mesh() = default;
@@ -167,24 +211,17 @@ class Mesh {
 
   virtual std::string get_mesh_type() const = 0;
 
-  //! Determine volume of materials within a single mesh elemenet
-  //
-  //! \param[in] n_sample Number of samples within each element
-  //! \param[in] bin Index of mesh element
-  //! \param[out] Array of (material index, volume) for desired element
-  //! \param[inout] seed Pseudorandom number seed
-  //! \return Number of materials within element
-  int material_volumes(int n_sample, int bin, gsl::span<MaterialVolume> volumes,
-    uint64_t* seed) const;
-
-  //! Determine volume of materials within a single mesh elemenet
+  //! Determine volume of materials within each mesh element
   //
-  //! \param[in] n_sample Number of samples within each element
-  //! \param[in] bin Index of mesh element
-  //! \param[inout] seed Pseudorandom number seed
-  //! \return Vector of (material index, volume) for desired element
-  vector<MaterialVolume> material_volumes(
-    int n_sample, int bin, uint64_t* seed) const;
+  //! \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
+  void material_volumes(int nx, int ny, int nz, int max_materials,
+    int32_t* materials, double* volumes) const;
 
   //! Determine bounding box of mesh
   //

openmc/lib/mesh.py

@@ -1,7 +1,7 @@
 from collections.abc import Mapping, Sequence
-from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER, Structure,
-                    create_string_buffer, c_uint64, c_size_t)
-from random import getrandbits
+from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER,
+                    create_string_buffer, c_size_t)
+from math import sqrt
 import sys
 from weakref import WeakValueDictionary
 
@@ -12,22 +12,18 @@
 from . import _dll
 from .core import _FortranObjectWithID
 from .error import _error_handler
-from .material import Material
 from .plot import _Position
 from ..bounding_box import BoundingBox
+from ..mesh import MeshMaterialVolumes
 
 __all__ = [
     'Mesh', 'RegularMesh', 'RectilinearMesh', 'CylindricalMesh',
-    'SphericalMesh', 'UnstructuredMesh', 'meshes'
+    'SphericalMesh', 'UnstructuredMesh', 'meshes', 'MeshMaterialVolumes'
 ]
 
 
-class _MaterialVolume(Structure):
-    _fields_ = [
-        ("material", c_int32),
-        ("volume", c_double)
-    ]
-
+arr_2d_int32 = np.ctypeslib.ndpointer(dtype=np.int32, ndim=2, flags='CONTIGUOUS')
+arr_2d_double = np.ctypeslib.ndpointer(dtype=np.double, ndim=2, flags='CONTIGUOUS')
 
 # Mesh functions
 _dll.openmc_extend_meshes.argtypes = [c_int32, c_char_p, POINTER(c_int32),
@@ -51,8 +47,7 @@ class _MaterialVolume(Structure):
 _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, POINTER(_MaterialVolume),
-    POINTER(c_int), POINTER(c_uint64)]
+    c_int32, c_int, c_int, c_int, c_int, arr_2d_int32, arr_2d_double]
 _dll.openmc_mesh_material_volumes.restype = c_int
 _dll.openmc_mesh_material_volumes.errcheck = _error_handler
 _dll.openmc_mesh_get_plot_bins.argtypes = [
@@ -190,58 +185,71 @@ def bounding_box(self) -> BoundingBox:
 
     def material_volumes(
             self,
-            n_samples: int = 10_000,
-            prn_seed: int | None = None
-    ) -> list[list[tuple[Material, float]]]:
-        """Determine volume of materials in each mesh element
+            n_samples: int | tuple[int, int, int] = 10_000,
+            max_materials: int = 4
+    ) -> MeshMaterialVolumes:
+        """Determine volume of materials in each mesh element.
+
+        This method works by raytracing repeatedly through the mesh to count the
+        estimated volume of each material in all mesh elements. Three sets of
+        rays are used: one set parallel to the x-axis, one parallel to the
+        y-axis, and one parallel to the z-axis.
 
         .. versionadded:: 0.15.0
 
+        .. versionchanged:: 0.15.1
+            Material volumes are now determined by raytracing rather than by
+            point sampling.
+
         Parameters
         ----------
-        n_samples : int
-            Number of samples in each mesh element
-        prn_seed : int
-            Pseudorandom number generator (PRNG) seed; if None, one will be
-            generated randomly.
+        n_samples : int or 3-tuple of int
+            Total number of rays to sample. The number of rays in each direction
+            is determined by the aspect ratio of the mesh bounding box. When
+            specified as a 3-tuple, it is interpreted as the number of rays in
+            the x, y, and z dimensions.
+        max_materials : int, optional
+            Estimated maximum number of materials in any given mesh element.
 
         Returns
         -------
-        List of tuple of (material, volume) for each mesh element. Void volume
-        is represented by having a value of None in the first element of a
-        tuple.
+        Dictionary-like object that maps material IDs to an array of volumes
+        equal in size to the number of mesh elements.
 
         """
-        if n_samples <= 0:
-            raise ValueError("Number of samples must be positive")
-        if prn_seed is None:
-            prn_seed = getrandbits(63)
-        prn_seed = c_uint64(prn_seed)
-
-        # Preallocate space for MaterialVolume results
-        size = 16
-        result = (_MaterialVolume * size)()
-
-        hits = c_int()  # Number of materials hit in a given element
-        volumes = []
-        for i_element in range(self.n_elements):
-            while True:
-                try:
-                    _dll.openmc_mesh_material_volumes(
-                        self._index, n_samples, i_element, size, result, hits, prn_seed)
-                except AllocationError:
-                    # Increase size of result array and try again
-                    size *= 2
-                    result = (_MaterialVolume * size)()
-                else:
-                    # If no error, break out of loop
-                    break
+        if isinstance(n_samples, int):
+            # Determine number of rays in each direction based on aspect ratios
+            # and using the relation (nx*ny + ny*nz + nx*nz) = n_samples
+            width_x, width_y, width_z = self.bounding_box.width
+            ax = width_x / width_z
+            ay = width_y / width_z
+            f = sqrt(n_samples/(ax*ay + ax + ay))
+            nx = round(f * ax)
+            ny = round(f * ay)
+            nz = round(f)
+        else:
+            nx, ny, nz = n_samples
+
+        # Preallocate arrays for material indices and volumes
+        n = self.n_elements
+        materials = np.full((n, max_materials), -2, dtype=np.int32)
+        volumes = np.zeros((n, max_materials), dtype=np.float64)
+
+        # Run material volume calculation
+        while True:
+            try:
+                _dll.openmc_mesh_material_volumes(
+                    self._index, nx, ny, nz, max_materials, materials, volumes)
+            except AllocationError:
+                # Increase size of result array and try again
+                max_materials *= 2
+                materials = np.full((n, max_materials), -2, dtype=np.int32)
+                volumes = np.zeros((n, max_materials), dtype=np.float64)
+            else:
+                # If no error, break out of loop
+                break
 
-            volumes.append([
-                (Material(index=r.material), r.volume)
-                for r in result[:hits.value]
-            ])
-        return volumes
+        return MeshMaterialVolumes(materials, volumes)
 
     def get_plot_bins(
             self,
@@ -719,7 +727,6 @@ def __getitem__(self, key):
             raise KeyError(str(e))
         return _get_mesh(index.value)
 
-
     def __iter__(self):
         for i in range(len(self)):
             yield _get_mesh(i).id