Ability to compute material volume fractions over mesh elements (open…

…mc-dev#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 @@ 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
+  //
+  //! \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 @@ def __set__(self, instance, value):
 
 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 @@ def __init__(self, uid=None, new=True, index=None):
         # 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 . 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_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 @@ def id(self):
     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.