Ability to compute material volume fractions over mesh elements (#2802)

include/openmc/capi.h

@@ -105,6 +105,9 @@ int openmc_mesh_filter_get_translation(int32_t index, double translation[3]);
 int openmc_mesh_filter_set_translation(int32_t index, double translation[3]);
 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_material_volumes(int32_t index, int n_sample, int bin,
+  int result_size, void* result, int* hits, uint64_t* seed);
 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

@@ -9,6 +9,7 @@
 #include "hdf5.h"
 #include "pugixml.hpp"
 #include "xtensor/xtensor.hpp"
+#include <gsl/gsl-lite.hpp>
 
 #include "openmc/error.h"
 #include "openmc/memory.h" // for unique_ptr
@@ -69,6 +70,12 @@ extern const libMesh::Parallel::Communicator* libmesh_comm;
 
 class Mesh {
 public:
+  // Types, aliases
+  struct MaterialVolume {
+    int32_t material; //!< material index
+    double volume;    //!< volume in [cm^3]
+  };
+
   // Constructors and destructor
   Mesh() = default;
   Mesh(pugi::xml_node node);
@@ -159,9 +166,28 @@ public:
 
   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
+  //
+  //! \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;
+
   // Data members
-  int id_ {-1};     //!< User-specified ID
-  int n_dimension_; //!< Number of dimensions
+  int id_ {-1};          //!< User-specified ID
+  int n_dimension_ {-1}; //!< Number of dimensions
 };
 
 class StructuredMesh : public Mesh {

include/openmc/volume_calc.h

@@ -1,18 +1,23 @@
 #ifndef OPENMC_VOLUME_CALC_H
 #define OPENMC_VOLUME_CALC_H
 
+#include <algorithm> // for find
 #include <cstdint>
 #include <string>
+#include <vector>
 
 #include "openmc/array.h"
+#include "openmc/openmp_interface.h"
 #include "openmc/position.h"
 #include "openmc/tallies/trigger.h"
 #include "openmc/vector.h"
 
 #include "pugixml.hpp"
 #include "xtensor/xtensor.hpp"
-
 #include <gsl/gsl-lite.hpp>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
 
 namespace openmc {
 
@@ -89,6 +94,35 @@ extern vector<VolumeCalculation> volume_calcs;
 // Non-member functions
 //==============================================================================
 
+//! Reduce vector of indices and hits from each thread to a single copy
+//
+//! \param[in] local_indices Indices specific to each thread
+//! \param[in] local_hits Hit count specific to each thread
+//! \param[out] indices Reduced vector of indices
+//! \param[out] hits Reduced vector of hits
+template<typename T, typename T2>
+void reduce_indices_hits(const vector<T>& local_indices,
+  const vector<T2>& local_hits, vector<T>& indices, vector<T2>& hits)
+{
+  const int n_threads = num_threads();
+
+#pragma omp for ordered schedule(static, 1)
+  for (int i = 0; i < n_threads; ++i) {
+#pragma omp ordered
+    for (int j = 0; j < local_indices.size(); ++j) {
+      // Check if this material has been added to the master list and if
+      // so, accumulate the number of hits
+      auto it = std::find(indices.begin(), indices.end(), local_indices[j]);
+      if (it == indices.end()) {
+        indices.push_back(local_indices[j]);
+        hits.push_back(local_hits[j]);
+      } else {
+        hits[it - indices.begin()] += local_hits[j];
+      }
+    }
+  }
+}
+
 void free_memory_volume();
 
 } // namespace openmc

openmc/lib/core.py

@@ -611,7 +611,7 @@ class _DLLGlobal:
 
 class _FortranObject:
     def __repr__(self):
-        return "{}[{}]".format(type(self).__name__, self._index)
+        return "<{}(index={})>".format(type(self).__name__, self._index)
 
 
 class _FortranObjectWithID(_FortranObject):
@@ -622,6 +622,9 @@ class _FortranObjectWithID(_FortranObject):
         # OutOfBoundsError will be raised here by virtue of referencing self.id
         self.id
 
+    def __repr__(self):
+        return "<{}(id={})>".format(type(self).__name__, self.id)
+
 
 @contextmanager
 def quiet_dll(output=True):

openmc/lib/mesh.py

@@ -1,6 +1,8 @@
 from collections.abc import Mapping
-from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER,
-                    create_string_buffer)
+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 typing import Optional, List, Tuple
 from weakref import WeakValueDictionary
 
 import numpy as np
@@ -10,9 +12,18 @@ from ..exceptions import AllocationError, InvalidIDError
 from . import _dll
 from .core import _FortranObjectWithID
 from .error import _error_handler
+from .material import Material
 
 __all__ = ['RegularMesh', 'RectilinearMesh', 'CylindricalMesh', 'SphericalMesh', 'UnstructuredMesh', 'meshes']
 
+
+class _MaterialVolume(Structure):
+    _fields_ = [
+        ("material", c_int32),
+        ("volume", c_double)
+    ]
+
+
 # Mesh functions
 _dll.openmc_extend_meshes.argtypes = [c_int32, c_char_p, POINTER(c_int32),
                                       POINTER(c_int32)]
@@ -24,6 +35,14 @@ _dll.openmc_mesh_get_id.errcheck = _error_handler
 _dll.openmc_mesh_set_id.argtypes = [c_int32, c_int32]
 _dll.openmc_mesh_set_id.restype = c_int
 _dll.openmc_mesh_set_id.errcheck = _error_handler
+_dll.openmc_mesh_get_n_elements.argtypes = [c_int32, POINTER(c_size_t)]
+_dll.openmc_mesh_get_n_elements.restype = c_int
+_dll.openmc_mesh_get_n_elements.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)]
+_dll.openmc_mesh_material_volumes.restype = c_int
+_dll.openmc_mesh_material_volumes.errcheck = _error_handler
 _dll.openmc_get_mesh_index.argtypes = [c_int32, POINTER(c_int32)]
 _dll.openmc_get_mesh_index.restype = c_int
 _dll.openmc_get_mesh_index.errcheck = _error_handler
@@ -123,6 +142,67 @@ class Mesh(_FortranObjectWithID):
     def id(self, mesh_id):
         _dll.openmc_mesh_set_id(self._index, mesh_id)
 
+    @property
+    def n_elements(self):
+        n = c_size_t()
+        _dll.openmc_mesh_get_n_elements(self._index, n)
+        return n.value
+
+    def material_volumes(
+            self,
+            n_samples: int = 10_000,
+            prn_seed: Optional[int] = None
+    ) -> List[List[Tuple[Material, float]]]:
+        """Determine volume of materials in each mesh element
+
+        .. versionadded:: 0.14.1
+
+        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.
+
+        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.
+
+        """
+        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
+
+            volumes.append([
+                (Material(index=r.material), r.volume)
+                for r in result[:hits.value]
+            ])
+        return volumes
+
 
 class RegularMesh(Mesh):
     """RegularMesh stored internally.

src/mesh.cpp

@@ -22,15 +22,19 @@
 #include "openmc/container_util.h"
 #include "openmc/error.h"
 #include "openmc/file_utils.h"
+#include "openmc/geometry.h"
 #include "openmc/hdf5_interface.h"
+#include "openmc/material.h"
 #include "openmc/memory.h"
 #include "openmc/message_passing.h"
 #include "openmc/openmp_interface.h"
+#include "openmc/particle_data.h"
 #include "openmc/random_dist.h"
 #include "openmc/search.h"
 #include "openmc/settings.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
+#include "openmc/volume_calc.h"
 #include "openmc/xml_interface.h"
 
 #ifdef LIBMESH
@@ -141,6 +145,92 @@ vector<double> Mesh::volumes() const
   return volumes;
 }
 
+int Mesh::material_volumes(
+  int n_sample, int bin, gsl::span<MaterialVolume> result, uint64_t* seed) const
+{
+  vector<int32_t> materials;
+  vector<int64_t> hits;
+
+#pragma omp parallel
+  {
+    vector<int32_t> local_materials;
+    vector<int64_t> local_hits;
+    GeometryState geom;
+
+#pragma omp for
+    for (int i = 0; i < n_sample; ++i) {
+      // Get seed for i-th sample
+      uint64_t seed_i = future_seed(3 * i, *seed);
+
+      // Sample position and set geometry state
+      geom.r() = this->sample_element(bin, &seed_i);
+      geom.u() = {1., 0., 0.};
+      geom.n_coord() = 1;
+
+      // If this location is not in the geometry at all, move on to next block
+      if (!exhaustive_find_cell(geom))
+        continue;
+
+      int i_material = geom.material();
+
+      // Check if this material was previously hit and if so, increment count
+      auto it =
+        std::find(local_materials.begin(), local_materials.end(), i_material);
+      if (it == local_materials.end()) {
+        local_materials.push_back(i_material);
+        local_hits.push_back(1);
+      } else {
+        local_hits[it - local_materials.begin()]++;
+      }
+    } // omp for
+
+    // Reduce index/hits lists from each thread into a single copy
+    reduce_indices_hits(local_materials, local_hits, materials, hits);
+  } // omp parallel
+
+  // Advance RNG seed
+  advance_prn_seed(3 * n_sample, seed);
+
+  // Make sure span passed in is large enough
+  if (hits.size() > result.size()) {
+    return -1;
+  }
+
+  // Convert hits to fractions
+  for (int i_mat = 0; i_mat < hits.size(); ++i_mat) {
+    double fraction = double(hits[i_mat]) / n_sample;
+    result[i_mat].material = materials[i_mat];
+    result[i_mat].volume = fraction * this->volume(bin);
+  }
+  return hits.size();
+}
+
+vector<Mesh::MaterialVolume> Mesh::material_volumes(
+  int n_sample, int bin, uint64_t* seed) const
+{
+  // Create result vector with space for 8 pairs
+  vector<Mesh::MaterialVolume> result;
+  result.reserve(8);
+
+  int size = -1;
+  while (true) {
+    // Get material volumes
+    size = this->material_volumes(
+      n_sample, bin, {result.data(), result.data() + result.capacity()}, seed);
+
+    // If capacity was sufficient, resize the vector and return
+    if (size >= 0) {
+      result.resize(size);
+      break;
+    }
+
+    // Otherwise, increase capacity of the vector
+    result.reserve(2 * result.capacity());
+  }
+
+  return result;
+}
+
 //==============================================================================
 // Structured Mesh implementation
 //==============================================================================
@@ -736,7 +826,7 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
     fatal_error("Must specify either <upper_right> or <width> on a mesh.");
   }
 
-  // Set volume fraction
+  // Set material volumes
   volume_frac_ = 1.0 / xt::prod(shape)();
 
   element_volume_ = 1.0;
@@ -1035,7 +1125,7 @@ double RectilinearMesh::volume(const MeshIndex& ijk) const
   double vol {1.0};
 
   for (int i = 0; i < n_dimension_; i++) {
-    vol *= grid_[i][ijk[i] + 1] - grid_[i][ijk[i]];
+    vol *= grid_[i][ijk[i]] - grid_[i][ijk[i] - 1];
   }
   return vol;
 }
@@ -1781,6 +1871,33 @@ extern "C" int openmc_mesh_set_id(int32_t index, int32_t id)
   return 0;
 }
 
+//! Get the number of elements in a mesh
+extern "C" int openmc_mesh_get_n_elements(int32_t index, size_t* n)
+{
+  if (int err = check_mesh(index))
+    return err;
+  *n = model::meshes[index]->n_bins();
+  return 0;
+}
+
+extern "C" int openmc_mesh_material_volumes(int32_t index, int n_sample,
+  int bin, int result_size, void* result, int* hits, uint64_t* seed)
+{
+  auto result_ = reinterpret_cast<Mesh::MaterialVolume*>(result);
+  if (!result_) {
+    set_errmsg("Invalid result pointer passed to openmc_mesh_material_volumes");
+    return OPENMC_E_INVALID_ARGUMENT;
+  }
+
+  if (int err = check_mesh(index))
+    return err;
+
+  int n = model::meshes[index]->material_volumes(
+    n_sample, bin, {result_, result_ + result_size}, seed);
+  *hits = n;
+  return (n == -1) ? OPENMC_E_ALLOCATE : 0;
+}
+
 //! Get the dimension of a regular mesh
 extern "C" int openmc_regular_mesh_get_dimension(
   int32_t index, int** dims, int* n)
@@ -1835,6 +1952,11 @@ extern "C" int openmc_regular_mesh_set_params(
     return err;
   RegularMesh* m = dynamic_cast<RegularMesh*>(model::meshes[index].get());
 
+  if (m->n_dimension_ == -1) {
+    set_errmsg("Need to set mesh dimension before setting parameters.");
+    return OPENMC_E_UNASSIGNED;
+  }
+
   vector<std::size_t> shape = {static_cast<std::size_t>(n)};
   if (ll && ur) {
     m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape);
@@ -1853,6 +1975,16 @@ extern "C" int openmc_regular_mesh_set_params(
     return OPENMC_E_INVALID_ARGUMENT;
   }
 
+  // Set material volumes
+
+  // TODO: incorporate this into method in RegularMesh that can be called from
+  // here and from constructor
+  m->volume_frac_ = 1.0 / xt::prod(m->get_x_shape())();
+  m->element_volume_ = 1.0;
+  for (int i = 0; i < m->n_dimension_; i++) {
+    m->element_volume_ *= m->width_[i];
+  }
+
   return 0;
 }
 

src/volume_calc.cpp

@@ -203,32 +203,9 @@ vector<VolumeCalculation::Result> VolumeCalculation::execute() const
       // At this point, each thread has its own pair of index/hits lists and we
       // now need to reduce them. OpenMP is not nearly smart enough to do this
       // on its own, so we have to manually reduce them
-      const int n_threads = num_threads();
-
-#pragma omp for ordered schedule(static)
-      for (int i = 0; i < n_threads; ++i) {
-#pragma omp ordered
-        for (int i_domain = 0; i_domain < n; ++i_domain) {
-          for (int j = 0; j < indices[i_domain].size(); ++j) {
-            // Check if this material has been added to the master list and if
-            // so, accumulate the number of hits
-            bool already_added = false;
-            for (int k = 0; k < master_indices[i_domain].size(); k++) {
-              if (indices[i_domain][j] == master_indices[i_domain][k]) {
-                master_hits[i_domain][k] += hits[i_domain][j];
-                already_added = true;
-                break;
-              }
-            }
-            if (!already_added) {
-              // If we made it here, the material hasn't yet been added to the
-              // master list, so add entries to the master indices and master
-              // hits lists
-              master_indices[i_domain].push_back(indices[i_domain][j]);
-              master_hits[i_domain].push_back(hits[i_domain][j]);
-            }
-          }
-        }
+      for (int i_domain = 0; i_domain < n; ++i_domain) {
+        reduce_indices_hits(indices[i_domain], hits[i_domain],
+          master_indices[i_domain], master_hits[i_domain]);
       }
     } // omp parallel
 

tests/unit_tests/test_lib.py

@@ -1,4 +1,5 @@
 from collections.abc import Mapping
+from math import pi
 import os
 
 import numpy as np
@@ -126,7 +127,7 @@ def test_cell(lib_init):
     cell = openmc.lib.cells[1]
     assert isinstance(cell.fill, openmc.lib.Material)
     cell.fill = openmc.lib.materials[1]
-    assert str(cell) == 'Cell[0]'
+    assert str(cell) == '<Cell(id=1)>'
     assert cell.name == "Fuel"
     cell.name = "Not fuel"
     assert cell.name == "Not fuel"
@@ -584,6 +585,25 @@ def test_regular_mesh(lib_init):
     msf.translation = translation
     assert msf.translation == translation
 
+    # Test material volumes
+    mesh = openmc.lib.RegularMesh()
+    mesh.dimension = (2, 2, 1)
+    mesh.set_parameters(lower_left=(-0.63, -0.63, -0.5),
+                        upper_right=(0.63, 0.63, 0.5))
+    vols = mesh.material_volumes()
+    assert len(vols) == 4
+    for elem_vols in vols:
+        assert sum(f[1] for f in elem_vols) == pytest.approx(1.26 * 1.26 / 4)
+
+    # If the mesh extends beyond the boundaries of the model, the volumes should
+    # still be reported correctly
+    mesh.dimension = (1, 1, 1)
+    mesh.set_parameters(lower_left=(-1.0, -1.0, -0.5),
+                        upper_right=(1.0, 1.0, 0.5))
+    vols = mesh.material_volumes(100_000)
+    for elem_vols in vols:
+        assert sum(f[1] for f in elem_vols) == pytest.approx(1.26 * 1.26, 1e-2)
+
 
 def test_rectilinear_mesh(lib_init):
     mesh = openmc.lib.RectilinearMesh()
@@ -604,7 +624,7 @@ def test_rectilinear_mesh(lib_init):
 
     meshes = openmc.lib.meshes
     assert isinstance(meshes, Mapping)
-    assert len(meshes) == 2
+    assert len(meshes) == 3
 
     mesh = meshes[mesh.id]
     assert isinstance(mesh, openmc.lib.RectilinearMesh)
@@ -615,8 +635,21 @@ def test_rectilinear_mesh(lib_init):
     msf = openmc.lib.MeshSurfaceFilter(mesh)
     assert msf.mesh == mesh
 
+    # Test material volumes
+    mesh = openmc.lib.RectilinearMesh()
+    w = 1.26
+    mesh.set_grid([-w/2, -w/4, w/2], [-w/2, -w/4, w/2], [-0.5, 0.5])
+
+    vols = mesh.material_volumes()
+    assert len(vols) == 4
+    assert sum(f[1] for f in vols[0]) == pytest.approx(w/4 * w/4)
+    assert sum(f[1] for f in vols[1]) == pytest.approx(w/4 * 3*w/4)
+    assert sum(f[1] for f in vols[2]) == pytest.approx(3*w/4 * w/4)
+    assert sum(f[1] for f in vols[3]) == pytest.approx(3*w/4 * 3*w/4)
+
+
 def test_cylindrical_mesh(lib_init):
-    deg2rad = lambda deg: deg*np.pi/180
+    deg2rad = lambda deg: deg*pi/180
     mesh = openmc.lib.CylindricalMesh()
     r_grid = [0., 5., 10.]
     phi_grid = np.radians([0., 10., 20.])
@@ -635,7 +668,7 @@ def test_cylindrical_mesh(lib_init):
 
     meshes = openmc.lib.meshes
     assert isinstance(meshes, Mapping)
-    assert len(meshes) == 3
+    assert len(meshes) == 5
 
     mesh = meshes[mesh.id]
     assert isinstance(mesh, openmc.lib.CylindricalMesh)
@@ -646,6 +679,21 @@ def test_cylindrical_mesh(lib_init):
     msf = openmc.lib.MeshSurfaceFilter(mesh)
     assert msf.mesh == mesh
 
+    # Test material volumes
+    mesh = openmc.lib.CylindricalMesh()
+    r_grid = (0., 0.25, 0.5)
+    phi_grid = np.linspace(0., 2.0*pi, 4)
+    z_grid = (-0.5, 0.5)
+    mesh.set_grid(r_grid, phi_grid, z_grid)
+
+    vols = mesh.material_volumes()
+    assert len(vols) == 6
+    for i in range(0, 6, 2):
+        assert sum(f[1] for f in vols[i]) == pytest.approx(pi * 0.25**2 / 3)
+    for i in range(1, 6, 2):
+        assert sum(f[1] for f in vols[i]) == pytest.approx(pi * (0.5**2 - 0.25**2) / 3)
+
+
 def test_spherical_mesh(lib_init):
     deg2rad = lambda deg: deg*np.pi/180
     mesh = openmc.lib.SphericalMesh()
@@ -666,7 +714,7 @@ def test_spherical_mesh(lib_init):
 
     meshes = openmc.lib.meshes
     assert isinstance(meshes, Mapping)
-    assert len(meshes) == 4
+    assert len(meshes) == 7
 
     mesh = meshes[mesh.id]
     assert isinstance(mesh, openmc.lib.SphericalMesh)
@@ -677,6 +725,24 @@ def test_spherical_mesh(lib_init):
     msf = openmc.lib.MeshSurfaceFilter(mesh)
     assert msf.mesh == mesh
 
+    # Test material volumes
+    mesh = openmc.lib.SphericalMesh()
+    r_grid = (0., 0.25, 0.5)
+    theta_grid = np.linspace(0., pi, 3)
+    phi_grid = np.linspace(0., 2.0*pi, 4)
+    mesh.set_grid(r_grid, theta_grid, phi_grid)
+
+    vols = mesh.material_volumes()
+    assert len(vols) == 12
+    d_theta = theta_grid[1] - theta_grid[0]
+    d_phi = phi_grid[1] - phi_grid[0]
+    for i in range(0, 12, 2):
+        assert sum(f[1] for f in vols[i]) == pytest.approx(
+            0.25**3 / 3 * d_theta * d_phi * 2/pi)
+    for i in range(1, 12, 2):
+        assert sum(f[1] for f in vols[i]) == pytest.approx(
+            (0.5**3 - 0.25**3) / 3 * d_theta * d_phi * 2/pi)
+
 
 def test_restart(lib_init, mpi_intracomm):
     # Finalize and re-init to make internal state consistent with XML.