Implement tally filter for filtering by reaction

CMakeLists.txt

@@ -429,6 +429,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_particle.cpp
   src/tallies/filter_particle_production.cpp
   src/tallies/filter_polar.cpp
+  src/tallies/filter_reaction.cpp
   src/tallies/filter_sph_harm.cpp
   src/tallies/filter_sptl_legendre.cpp
   src/tallies/filter_surface.cpp

docs/source/pythonapi/base.rst

@@ -149,6 +149,7 @@ Constructing Tallies
    openmc.ZernikeRadialFilter
    openmc.ParentNuclideFilter
    openmc.ParticleFilter
+   openmc.ReactionFilter
    openmc.MeshMaterialVolumes
    openmc.Trigger
    openmc.TallyDerivative

docs/source/pythonapi/capi.rst

@@ -82,7 +82,9 @@ Classes
    Nuclide
    ParentNuclideFilter
    ParticleFilter
+   ParticleProductionFilter
    PolarFilter
+   ReactionFilter
    RectilinearMesh
    RegularMesh
    SpatialLegendreFilter

include/openmc/constants.h

@@ -227,6 +227,7 @@ enum ReactionType {
   N_XA = 207,
   HEATING = 301,
   DAMAGE_ENERGY = 444,
+  PHOTON_TOTAL = 501,
   COHERENT = 502,
   INCOHERENT = 504,
   PAIR_PROD_ELEC = 515,

include/openmc/endf.h

@@ -33,6 +33,13 @@ bool is_disappearance(int MT);
 //! \return Whether corresponding reaction is an inelastic scattering reaction
 bool is_inelastic_scatter(int MT);
 
+//! Determine whether an MT number matches a target MT, considering that the
+//! target may be a summation reaction.
+//! \param[in] event_mt MT number of the actual event
+//! \param[in] target_mt MT number to check against
+//! \return Whether event_mt is a component of target_mt (or equal to it)
+bool mt_matches(int event_mt, int target_mt);
+
 //==============================================================================
 //! Abstract one-dimensional function
 //==============================================================================

include/openmc/reaction.h

@@ -76,11 +76,21 @@ class Reaction {
 //! \return Name of the corresponding reaction
 std::string reaction_name(int mt);
 
-//! Return reaction type (MT value) given a reaction name
+//! Return MT value for given a reaction name (including special tally MT
+//! values)
 //
 //! \param[in] name  Reaction name
-//! \return Corresponding reaction type (MT value)
-int reaction_type(std::string name);
+//! \return Corresponding MT number or special tally score
+int reaction_tally_mt(std::string name);
+
+//! Return ENDF MT number given a reaction name
+//
+//! Unlike reaction_tally_mt(), this function always returns a positive ENDF MT
+//! number and never a special negative tally score.
+//!
+//! \param[in] name  Reaction name
+//! \return Corresponding ENDF MT number
+int reaction_mt(const std::string& name);
 
 } // namespace openmc
 

include/openmc/tallies/filter.h

@@ -41,6 +41,7 @@ enum class FilterType {
   PARTICLE,
   PARTICLE_PRODUCTION,
   POLAR,
+  REACTION,
   SPHERICAL_HARMONICS,
   SPATIAL_LEGENDRE,
   SURFACE,

include/openmc/tallies/filter_reaction.h

@@ -0,0 +1,52 @@
+#ifndef OPENMC_TALLIES_FILTER_REACTION_H
+#define OPENMC_TALLIES_FILTER_REACTION_H
+
+#include "openmc/span.h"
+#include "openmc/tallies/filter.h"
+#include "openmc/vector.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Bins tally events based on the reaction type (MT number).
+//==============================================================================
+
+class ReactionFilter : public Filter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  ~ReactionFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "reaction"; }
+  FilterType type() const override { return FilterType::REACTION; }
+
+  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
+
+  const vector<int>& bins() const { return bins_; }
+  void set_bins(span<const int> bins);
+
+protected:
+  //----------------------------------------------------------------------------
+  // Data members
+
+  //! MT numbers to match
+  vector<int> bins_;
+};
+
+} // namespace openmc
+
+#endif // OPENMC_TALLIES_FILTER_REACTION_H

openmc/data/reaction.py

@@ -21,6 +21,7 @@
 from .kalbach_mann import KalbachMann
 from .laboratory import LaboratoryAngleEnergy
 from .nbody import NBodyPhaseSpace
+from .photon import _SUBSHELLS
 from .product import Product
 from .uncorrelated import UncorrelatedAngleEnergy
 
@@ -54,6 +55,10 @@
                  198: '(n,n3p)', 199: '(n,3n2pa)', 200: '(n,5n2p)', 203: '(n,Xp)',
                  204: '(n,Xd)', 205: '(n,Xt)', 206: '(n,X3He)', 207: '(n,Xa)',
                  301: 'heating', 444: 'damage-energy',
+                 501: 'photon-total', 502: 'coherent-scatter',
+                 504: 'incoherent-scatter', 515: 'pair-production-electron',
+                 516: 'pair-production', 517: 'pair-production-nuclear',
+                 522: 'photoelectric',
                  649: '(n,pc)', 699: '(n,dc)', 749: '(n,tc)', 799: '(n,3Hec)',
                  849: '(n,ac)', 891: '(n,2nc)', 901: 'heating-local'}
 REACTION_NAME.update({i: f'(n,n{i - 50})' for i in range(51, 91)})
@@ -63,9 +68,16 @@
 REACTION_NAME.update({i: f'(n,3He{i - 750})' for i in range(750, 799)})
 REACTION_NAME.update({i: f'(n,a{i - 800})' for i in range(800, 849)})
 REACTION_NAME.update({i: f'(n,2n{i - 875})' for i in range(875, 891)})
+REACTION_NAME.update(
+    {534 + i: f'photoelectric-{shell}' for i, shell in enumerate(_SUBSHELLS[1:])}
+)
 
 REACTION_MT = {name: mt for mt, name in REACTION_NAME.items()}
+REACTION_MT['total'] = 1
+REACTION_MT['elastic'] = 2
 REACTION_MT['fission'] = 18
+REACTION_MT['absorption'] = 27
+REACTION_MT['capture'] = 102
 
 FISSION_MTS = (18, 19, 20, 21, 38)
 

openmc/filter.py

@@ -13,6 +13,7 @@
 import openmc
 import openmc.checkvalue as cv
 from .cell import Cell
+from .data.reaction import REACTION_NAME, REACTION_MT
 from .material import Material
 from .mixin import IDManagerMixin
 from .surface import Surface
@@ -22,11 +23,11 @@
 
 _FILTER_TYPES = (
     'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy',
-    'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
-    'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre',
-    'sphericalharmonics', 'zernike', 'zernikeradial', 'particle',
+    'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell',
+    'delayedgroup', 'energyfunction', 'cellfrom', 'materialfrom', 'legendre',
+    'spatiallegendre', 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle',
     'particleproduction', 'cellinstance', 'collision', 'time', 'parentnuclide',
-    'weight', 'meshborn', 'meshsurface', 'meshmaterial',
+    'weight', 'meshborn', 'meshsurface', 'meshmaterial', 'reaction',
 )
 
 _CURRENT_NAMES = (
@@ -1415,6 +1416,80 @@ def check_bins(self, bins):
             cv.check_greater_than('filter value', x, 0, equality=True)
 
 
+class ReactionFilter(Filter):
+    """Bins tally events based on the reaction type (MT number).
+
+    .. versionadded:: 0.15.4
+
+    Parameters
+    ----------
+    bins : str, int, or iterable thereof
+        The reaction types to tally. Can be reaction name strings
+        (e.g., ``'(n,elastic)'``, ``'(n,gamma)'``) or integer MT numbers
+        (e.g., 2, 102). Integer MT values are automatically converted to their
+        canonical string representation.
+    filter_id : int
+        Unique identifier for the filter
+
+    Attributes
+    ----------
+    bins : numpy.ndarray of str
+        Reaction name strings
+    id : int
+        Unique identifier for the filter
+    num_bins : int
+        The number of filter bins
+
+    """
+
+    def __init__(self, bins, filter_id=None):
+        self.bins = bins
+        self.id = filter_id
+
+    @Filter.bins.setter
+    def bins(self, bins):
+        if isinstance(bins, (str, Integral)):
+            bins = [bins]
+        elif not isinstance(bins, list):
+            bins = list(bins)
+        normalized = []
+        for b in bins:
+            if isinstance(b, Integral):
+                if int(b) not in REACTION_NAME:
+                    raise ValueError(f"No known reaction for MT={b}")
+                normalized.append(REACTION_NAME[int(b)])
+            elif isinstance(b, str):
+                if b == 'total':
+                    warnings.warn(
+                        "The reaction name 'total' is ambiguous. Use "
+                        "'(n,total)' for neutron total cross section or "
+                        "'photon-total' for photon total. Interpreting as"
+                        "'(n,total)'.")
+                if b not in REACTION_MT:
+                    raise ValueError(f"Unknown reaction name '{b}'")
+                normalized.append(REACTION_NAME[REACTION_MT[b]])
+            else:
+                raise TypeError(f"Expected str or int for reaction filter "
+                                f"bin, got {type(b)}")
+        self._bins = np.array(normalized, dtype=str)
+
+    @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")
+        bins = [b.decode() for b in group['bins'][()]]
+        filter_id = int(group.name.split('/')[-1].lstrip('filter '))
+        return cls(bins, filter_id=filter_id)
+
+    @classmethod
+    def from_xml_element(cls, elem, **kwargs):
+        filter_id = int(get_text(elem, "id"))
+        bins = get_elem_list(elem, "bins", str) or []
+        return cls(bins, filter_id=filter_id)
+
+
 class RealFilter(Filter):
     """Tally modifier that describes phase-space and other characteristics
 

openmc/lib/filter.py

@@ -22,9 +22,10 @@
     'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter',
     'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshBornFilter',
     'MeshMaterialFilter', 'MeshSurfaceFilter', 'MuFilter', 'MuSurfaceFilter',
-    'ParentNuclideFilter', 'ParticleFilter', 'PolarFilter', 'SphericalHarmonicsFilter',
-    'SpatialLegendreFilter', 'SurfaceFilter', 'TimeFilter', 'UniverseFilter',
-    'WeightFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters'
+    'ParentNuclideFilter', 'ParticleFilter', 'ParticleProductionFilter', 'PolarFilter',
+    'ReactionFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter',
+    'SurfaceFilter', 'TimeFilter', 'UniverseFilter', 'WeightFilter', 'ZernikeFilter',
+    'ZernikeRadialFilter', 'filters'
 ]
 
 # Tally functions
@@ -607,10 +608,18 @@ def bins(self):
         return [ParticleType(i) for i in particle_i]
 
 
+class ParticleProductionFilter(Filter):
+    filter_type = 'particleproduction'
+
+
 class PolarFilter(Filter):
     filter_type = 'polar'
 
 
+class ReactionFilter(Filter):
+    filter_type = 'reaction'
+
+
 class SphericalHarmonicsFilter(Filter):
     filter_type = 'sphericalharmonics'
 
@@ -711,7 +720,9 @@ class ZernikeRadialFilter(ZernikeFilter):
     'musurface': MuSurfaceFilter,
     'parentnuclide': ParentNuclideFilter,
     'particle': ParticleFilter,
+    'particleproduction': ParticleProductionFilter,
     'polar': PolarFilter,
+    'reaction': ReactionFilter,
     'sphericalharmonics': SphericalHarmonicsFilter,
     'spatiallegendre': SpatialLegendreFilter,
     'surface': SurfaceFilter,

src/chain.cpp

@@ -42,7 +42,7 @@ ChainNuclide::ChainNuclide(pugi::xml_node node)
       branching_ratio =
         std::stod(get_node_value(reaction_node, "branching_ratio"));
     }
-    int mt = reaction_type(rx_name);
+    int mt = reaction_mt(rx_name);
     reaction_products_[mt].push_back({rx_target, branching_ratio});
   }