Implement a new MeshMaterialFilter (#3406)

Co-authored-by: Jonathan Shimwell <drshimwell@gmail.com>
Co-authored-by: Patrick Shriwise <pshriwise@gmail.com>

Author: 3 people

CMakeLists.txt

@@ -414,6 +414,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_materialfrom.cpp
   src/tallies/filter_mesh.cpp
   src/tallies/filter_meshborn.cpp
+  src/tallies/filter_meshmaterial.cpp
   src/tallies/filter_meshsurface.cpp
   src/tallies/filter_mu.cpp
   src/tallies/filter_musurface.cpp

docs/source/pythonapi/base.rst

@@ -129,6 +129,7 @@ Constructing Tallies
    openmc.SurfaceFilter
    openmc.MeshFilter
    openmc.MeshBornFilter
+   openmc.MeshMaterialFilter
    openmc.MeshSurfaceFilter
    openmc.EnergyFilter
    openmc.EnergyoutFilter

include/openmc/tallies/filter.h

@@ -33,6 +33,7 @@ enum class FilterType {
   MATERIALFROM,
   MESH,
   MESHBORN,
+  MESH_MATERIAL,
   MESH_SURFACE,
   MU,
   MUSURFACE,

include/openmc/tallies/filter_mesh.h

@@ -9,9 +9,9 @@
 namespace openmc {
 
 //==============================================================================
-//! Indexes the location of particle events to a regular mesh.  For tracklength
-//! tallies, it will produce multiple valid bins and the bin weight will
-//! correspond to the fraction of the track length that lies in that bin.
+//! Indexes the location of particle events to a mesh. For tracklength tallies,
+//! it will produce multiple valid bins and the bin weight will correspond to
+//! the fraction of the track length that lies in that bin.
 //==============================================================================
 
 class MeshFilter : public Filter {

include/openmc/tallies/filter_meshmaterial.h

@@ -0,0 +1,114 @@
+#ifndef OPENMC_TALLIES_FILTER_MESHMATERIAL_H
+#define OPENMC_TALLIES_FILTER_MESHMATERIAL_H
+
+#include <cstdint>
+#include <string>
+#include <unordered_map>
+#include <unordered_set>
+
+#include "openmc/position.h"
+#include "openmc/random_ray/source_region.h"
+#include "openmc/span.h"
+#include "openmc/tallies/filter.h"
+#include "openmc/vector.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Helper structs that define a combination of a mesh element index and a
+//! material index and a functor for hashing to place in an unordered_map
+//==============================================================================
+
+struct ElementMat {
+  //! Check for equality
+  bool operator==(const ElementMat& other) const
+  {
+    return index_element == other.index_element && index_mat == other.index_mat;
+  }
+
+  int32_t index_element;
+  int32_t index_mat;
+};
+
+struct ElementMatHash {
+  std::size_t operator()(const ElementMat& k) const
+  {
+    size_t seed = 0;
+    hash_combine(seed, k.index_element);
+    hash_combine(seed, k.index_mat);
+    return seed;
+  }
+};
+
+//==============================================================================
+//! Indexes the location of particle events to combinations of mesh element
+//! index and material.  For tracklength tallies, it will produce multiple valid
+//! bins and the bin weight will correspond to the fraction of the track length
+//! that lies in that bin.
+//==============================================================================
+
+class MeshMaterialFilter : public Filter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  ~MeshMaterialFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "meshmaterial"; }
+  FilterType type() const override { return FilterType::MESH_MATERIAL; }
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator,
+    FilterMatch& match) const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //----------------------------------------------------------------------------
+  // Accessors
+
+  int32_t mesh() const { return mesh_; }
+
+  void set_mesh(int32_t mesh);
+
+  //! Set the bins based on a flat vector of alternating element index and
+  //! material IDs
+  void set_bins(span<int32_t> bins);
+
+  //! Set the bins based on a vector of (element, material index) pairs
+  void set_bins(vector<ElementMat>&& bins);
+
+  virtual void set_translation(const Position& translation);
+
+  virtual void set_translation(const double translation[3]);
+
+  virtual const Position& translation() const { return translation_; }
+
+  virtual bool translated() const { return translated_; }
+
+private:
+  //----------------------------------------------------------------------------
+  // Data members
+
+  int32_t mesh_;            //!< Index of the mesh
+  bool translated_ {false}; //!< Whether or not the filter is translated
+  Position translation_ {0.0, 0.0, 0.0}; //!< Filter translation
+
+  //! The indices of the mesh element-material combinations binned by this
+  //! filter.
+  vector<ElementMat> bins_;
+
+  //! The set of materials used in this filter
+  std::unordered_set<int32_t> materials_;
+
+  //! A map from mesh element-material indices to filter bin indices.
+  std::unordered_map<ElementMat, int32_t, ElementMatHash> map_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_MESHMATERIAL_H

openmc/deplete/microxs.py

@@ -56,8 +56,8 @@ def get_microxs_and_flux(
     ----------
     model : openmc.Model
         OpenMC model object. Must contain geometry, materials, and settings.
-    domains : list of openmc.Material or openmc.Cell or openmc.Universe, or openmc.MeshBase
-        Domains in which to tally reaction rates.
+    domains : list of openmc.Material or openmc.Cell or openmc.Universe, or openmc.MeshBase, or openmc.Filter
+        Domains in which to tally reaction rates, or a spatial tally filter.
     nuclides : list of str
         Nuclides to get cross sections for. If not specified, all burnable
         nuclides from the depletion chain file are used.
@@ -103,7 +103,10 @@ def get_microxs_and_flux(
         energy_filter = openmc.EnergyFilter.from_group_structure(energies)
     else:
         energy_filter = openmc.EnergyFilter(energies)
-    if isinstance(domains, openmc.MeshBase):
+
+    if isinstance(domains, openmc.Filter):
+        domain_filter = domains
+    elif isinstance(domains, openmc.MeshBase):
         domain_filter = openmc.MeshFilter(domains)
     elif isinstance(domains[0], openmc.Material):
         domain_filter = openmc.MaterialFilter(domains)

openmc/filter.py

@@ -25,7 +25,8 @@
     'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
     'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre',
     'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance',
-    'collision', 'time', 'parentnuclide', 'weight'
+    'collision', 'time', 'parentnuclide', 'weight', 'meshborn', 'meshsurface',
+    'meshmaterial',
 )
 
 _CURRENT_NAMES = (
@@ -900,8 +901,6 @@ def mesh(self, mesh):
 
     @property
     def shape(self):
-        if isinstance(self, MeshSurfaceFilter):
-            return (self.num_bins,)
         return self.mesh.dimension
 
     @property
@@ -992,8 +991,11 @@ def to_xml_element(self):
             XML element containing filter data
 
         """
-        element = super().to_xml_element()
-        element[0].text = str(self.mesh.id)
+        element = ET.Element('filter')
+        element.set('id', str(self.id))
+        element.set('type', self.short_name.lower())
+        subelement = ET.SubElement(element, 'bins')
+        subelement.text = str(self.mesh.id)
         if self.translation is not None:
             element.set('translation', ' '.join(map(str, self.translation)))
         return element
@@ -1039,6 +1041,181 @@ class MeshBornFilter(MeshFilter):
     """
 
 
+class MeshMaterialFilter(MeshFilter):
+    """Filter events by combinations of mesh elements and materials.
+
+    .. versionadded:: 0.15.3
+
+    Parameters
+    ----------
+    mesh : openmc.MeshBase
+        The mesh object that events will be tallied onto
+    bins : iterable of 2-tuples or numpy.ndarray
+        Combinations of (mesh element, material) to tally, given as 2-tuples.
+        The first value in the tuple represents the index of the mesh element,
+        and the second value indicates the material (either a
+        :class:`openmc.Material` instance of the ID).
+    filter_id : int
+        Unique identifier for the filter
+
+    """
+    def __init__(self, mesh: openmc.MeshBase, bins, filter_id=None):
+        self.mesh = mesh
+        self.bins = bins
+        self.id = filter_id
+        self._translation = None
+
+    @classmethod
+    def from_volumes(cls, mesh: openmc.MeshBase, volumes: openmc.MeshMaterialVolumes):
+        """Construct a MeshMaterialFilter from a MeshMaterialVolumes object.
+
+        Parameters
+        ----------
+        mesh : openmc.MeshBase
+            The mesh object that events will be tallied onto
+        volumes : openmc.MeshMaterialVolumes
+            The mesh material volumes to use for the filter
+
+        Returns
+        -------
+        MeshMaterialFilter
+            A new MeshMaterialFilter instance
+
+        """
+        bins = list(zip(*np.where(volumes._materials > -1)))
+        return cls(mesh, bins)
+
+    def __hash__(self):
+        data = (type(self).__name__, self.mesh.id, tuple(self.bins.ravel()))
+        return hash(data)
+
+    def __repr__(self):
+        string = type(self).__name__ + '\n'
+        string += '{: <16}=\t{}\n'.format('\tID', self.id)
+        string += '{: <16}=\t{}\n'.format('\tMesh ID', self.mesh.id)
+        string += '{: <16}=\n{}\n'.format('\tBins', self.bins)
+        string += '{: <16}=\t{}\n'.format('\tTranslation', self.translation)
+        return string
+
+    @property
+    def shape(self):
+        return (self.num_bins,)
+
+    @property
+    def mesh(self):
+        return self._mesh
+
+    @mesh.setter
+    def mesh(self, mesh):
+        cv.check_type('filter mesh', mesh, openmc.MeshBase)
+        self._mesh = mesh
+
+    @Filter.bins.setter
+    def bins(self, bins):
+        pairs = np.empty((len(bins), 2), dtype=int)
+        for i, (elem, mat) in enumerate(bins):
+            cv.check_type('element', elem, Integral)
+            cv.check_type('material', mat, (Integral, openmc.Material))
+            pairs[i, 0] = elem
+            pairs[i, 1] = mat if isinstance(mat, Integral) else mat.id
+        self._bins = pairs
+
+    def to_xml_element(self):
+        """Return XML element representing the filter.
+
+        Returns
+        -------
+        element : lxml.etree._Element
+            XML element containing filter data
+
+        """
+        element = ET.Element('filter')
+        element.set('id', str(self.id))
+        element.set('type', self.short_name.lower())
+        element.set('mesh', str(self.mesh.id))
+
+        if self.translation is not None:
+            element.set('translation', ' '.join(map(str, self.translation)))
+
+        subelement = ET.SubElement(element, 'bins')
+        subelement.text = ' '.join(str(i) for i in self.bins.ravel())
+
+        return element
+
+    @classmethod
+    def from_xml_element(cls, elem: ET.Element, **kwargs) -> MeshMaterialFilter:
+        filter_id = int(elem.get('id'))
+        mesh_id = int(elem.get('mesh'))
+        mesh_obj = kwargs['meshes'][mesh_id]
+        bins = [int(x) for x in get_text(elem, 'bins').split()]
+        bins = list(zip(bins[::2], bins[1::2]))
+        out = cls(mesh_obj, bins, filter_id=filter_id)
+
+        translation = elem.get('translation')
+        if translation:
+            out.translation = [float(x) for x in translation.split()]
+        return out
+
+    @classmethod
+    def from_hdf5(cls, group, **kwargs):
+        if group['type'][()].decode() != cls.short_name.lower():
+            raise ValueError("Expected HDF5 data for filter type '"
+                             + cls.short_name.lower() + "' but got '"
+                             + group['type'][()].decode() + " instead")
+
+        if 'meshes' not in kwargs:
+            raise ValueError(cls.__name__ + " requires a 'meshes' keyword "
+                             "argument.")
+
+        mesh_id = group['mesh'][()]
+        mesh_obj = kwargs['meshes'][mesh_id]
+        bins = group['bins'][()]
+        filter_id = int(group.name.split('/')[-1].lstrip('filter '))
+        out = cls(mesh_obj, bins, filter_id=filter_id)
+
+        translation = group.get('translation')
+        if translation:
+            out.translation = translation[()]
+
+        return out
+
+    def get_pandas_dataframe(self, data_size, stride, **kwargs):
+        """Builds a Pandas DataFrame for the Filter's bins.
+
+        This method constructs a Pandas DataFrame object for the filter with
+        columns annotated by filter bin information. This is a helper method for
+        :meth:`Tally.get_pandas_dataframe`.
+
+        Parameters
+        ----------
+        data_size : int
+            The total number of bins in the tally corresponding to this filter
+        stride : int
+            Stride in memory for the filter
+
+        Returns
+        -------
+        pandas.DataFrame
+            A Pandas DataFrame with a multi-index column for the cell instance.
+            The number of rows in the DataFrame is the same as the total number
+            of bins in the corresponding tally, with the filter bin appropriately
+            tiled to map to the corresponding tally bins.
+
+        See also
+        --------
+        Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
+
+        """
+        # Repeat and tile bins as necessary to account for other filters.
+        bins = np.repeat(self.bins, stride, axis=0)
+        tile_factor = data_size // len(bins)
+        bins = np.tile(bins, (tile_factor, 1))
+
+        columns = pd.MultiIndex.from_product([[self.short_name.lower()],
+                                              ['element', 'material']])
+        return pd.DataFrame(bins, columns=columns)
+
+
 class MeshSurfaceFilter(MeshFilter):
     """Filter events by surface crossings on a mesh.
 
@@ -1065,6 +1242,9 @@ class MeshSurfaceFilter(MeshFilter):
         The number of filter bins
 
     """
+    @property
+    def shape(self):
+        return (self.num_bins,)
 
     @MeshFilter.mesh.setter
     def mesh(self, mesh):