Secondary energy filter

docs/source/pythonapi/base.rst

@@ -133,6 +133,7 @@ Constructing Tallies
    openmc.MeshSurfaceFilter
    openmc.EnergyFilter
    openmc.EnergyoutFilter
+   openmc.SecondaryEnergyFilter
    openmc.MuFilter
    openmc.MuSurfaceFilter
    openmc.PolarFilter

include/openmc/constants.h

@@ -290,7 +290,7 @@ enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION };
 enum class TallyEvent { SURFACE, LATTICE, KILL, SCATTER, ABSORB };
 
 // Tally score type -- if you change these, make sure you also update the
-// _SCORES dictionary in openmc/capi/tally.py
+// _SCORES dictionary in openmc/lib/tally.py
 //
 // These are kept as a normal enum and made negative, since variables which
 // store one of these enum values usually also may be responsible for storing

include/openmc/particle_data.h

@@ -466,6 +466,11 @@ class ParticleData : public GeometryState {
 
   vector<SourceSite> secondary_bank_;
 
+  // Keeps track of how many secondary particles were created
+  // in the collision the particle most recently experienced.
+  // This is used to tally secondary photon energies.
+  int secondaries_this_collision_ {0};
+
   int64_t current_work_;
 
   vector<double> flux_derivs_;
@@ -615,7 +620,17 @@ class ParticleData : public GeometryState {
 
   // secondary particle bank
   SourceSite& secondary_bank(int i) { return secondary_bank_[i]; }
+  const SourceSite& secondary_bank(int i) const { return secondary_bank_[i]; }
   decltype(secondary_bank_)& secondary_bank() { return secondary_bank_; }
+  decltype(secondary_bank_) const& secondary_bank() const
+  {
+    return secondary_bank_;
+  }
+  const int& secondaries_this_collision() const
+  {
+    return secondaries_this_collision_;
+  }
+  int& secondaries_this_collision() { return secondaries_this_collision_; }
 
   // Current simulation work index
   int64_t& current_work() { return current_work_; }

include/openmc/tallies/filter.h

@@ -28,6 +28,7 @@ enum class FilterType {
   ENERGY_FUNCTION,
   ENERGY,
   ENERGY_OUT,
+  SECONDARY_ENERGY,
   LEGENDRE,
   MATERIAL,
   MATERIALFROM,

include/openmc/tallies/filter_energy.h

@@ -1,6 +1,7 @@
 #ifndef OPENMC_TALLIES_FILTER_ENERGY_H
 #define OPENMC_TALLIES_FILTER_ENERGY_H
 
+#include "openmc/particle.h"
 #include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
@@ -72,5 +73,38 @@ class EnergyoutFilter : public EnergyFilter {
   std::string text_label(int bin) const override;
 };
 
+//==============================================================================
+//! Bins the outgoing energy of secondary particles
+//!
+//! This is used to get the photon production matrix for multigroup photon
+//! transport. It could also be used to find the energy distribution of
+//! neutron secondaries or others, for example.
+//!
+//! Using anything other than analog estimators here would be complicated
+//==============================================================================
+
+class SecondaryEnergyFilter : public EnergyFilter {
+public:
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "secondaryenergy"; }
+  FilterType type() const override { return FilterType::SECONDARY_ENERGY; }
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator,
+    FilterMatch& match) const override;
+
+  std::string text_label(int bin) const override;
+
+  void from_xml(pugi::xml_node node) override;
+
+protected:
+  // This filter could simultaneously use different particle types, for
+  // example creating one set of energy bins for photons and another for
+  // electrons. However, for typical use only photons are of interest.
+  // Covering multiple particles can be done by defining separate tallies.
+  ParticleType secondary_type_ {ParticleType::photon};
+};
+
 } // namespace openmc
 #endif // OPENMC_TALLIES_FILTER_ENERGY_H

openmc/filter.py

@@ -127,9 +127,9 @@ def __eq__(self, other):
     def __gt__(self, other):
         if type(self) is not type(other):
             if self.short_name in _FILTER_TYPES and \
-                other.short_name in _FILTER_TYPES:
+                    other.short_name in _FILTER_TYPES:
                 delta = _FILTER_TYPES.index(self.short_name) - \
-                        _FILTER_TYPES.index(other.short_name)
+                    _FILTER_TYPES.index(other.short_name)
                 return delta > 0
             else:
                 return False
@@ -278,7 +278,6 @@ def from_xml_element(cls, elem, **kwargs):
             if filter_type == subclass.short_name.lower():
                 return subclass.from_xml_element(elem, **kwargs)
 
-
     def can_merge(self, other):
         """Determine if filter can be merged with another.
 
@@ -429,6 +428,7 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs):
 
 class WithIDFilter(Filter):
     """Abstract parent for filters of types with IDs (Cell, Material, etc.)."""
+
     def __init__(self, bins, filter_id=None):
         bins = np.atleast_1d(bins)
 
@@ -632,6 +632,7 @@ class CellInstanceFilter(Filter):
     DistribcellFilter
 
     """
+
     def __init__(self, bins, filter_id=None):
         self.bins = bins
         self.id = filter_id
@@ -753,6 +754,7 @@ class ParticleFilter(Filter):
         The number of filter bins
 
     """
+
     def __eq__(self, other):
         if type(self) is not type(other):
             return False
@@ -1059,6 +1061,7 @@ class MeshMaterialFilter(MeshFilter):
         Unique identifier for the filter
 
     """
+
     def __init__(self, mesh: openmc.MeshBase, bins, filter_id=None):
         self.mesh = mesh
         self.bins = bins
@@ -1388,6 +1391,7 @@ class RealFilter(Filter):
         The number of filter bins
 
     """
+
     def __init__(self, values, filter_id=None):
         self.values = np.asarray(values)
         self.bins = np.vstack((self.values[:-1], self.values[1:])).T
@@ -1663,7 +1667,8 @@ def from_group_structure(cls, group_structure):
 
         """
 
-        cv.check_value('group_structure', group_structure, openmc.mgxs.GROUP_STRUCTURES.keys())
+        cv.check_value('group_structure', group_structure,
+                       openmc.mgxs.GROUP_STRUCTURES.keys())
         return cls(openmc.mgxs.GROUP_STRUCTURES[group_structure.upper()])
 
 
@@ -1693,6 +1698,66 @@ class EnergyoutFilter(EnergyFilter):
 
     """
 
+
+class SecondaryEnergyFilter(EnergyFilter):
+    """Bins tally events based on energy of secondary particles.
+
+    .. versionadded:: 0.15.3
+
+    This filter collects the energies of secondary particles (e.g., photons
+    or electrons) produced in a reaction. This is useful for constructing
+    production matrices or analyzing secondary particle spectra.
+
+    The primary particle type should be filtered by the usual ParticleFilter.
+
+    Parameters
+    ----------
+    values : Iterable of Real
+        A list of energy boundaries in [eV]; each successive pair defines a bin.
+    filter_id : int, optional
+        Unique identifier for the filter
+    particle : str
+        Type of secondary particle to tally ('photon', 'neutron', etc.)
+
+    Attributes
+    ----------
+    values : numpy.ndarray
+        Energy boundaries in [eV]
+    bins : numpy.ndarray
+        Array of (low, high) energy bin pairs
+    num_bins : int
+        Number of filter bins
+    particle : str
+        The secondary particle type this filter applies to
+    """
+
+    def __init__(self, values, filter_id=None, particle='photon'):
+        self.particle = particle
+        super().__init__(values, filter_id)
+
+    @property
+    def particle(self):
+        return self._particle
+
+    @particle.setter
+    def particle(self, particle):
+        cv.check_value('particle', particle, _PARTICLES)
+        self._particle = particle
+
+    def to_xml_element(self):
+        element = super().to_xml_element()
+        subelement = ET.SubElement(element, 'particle')
+        subelement.text = self.particle
+        return element
+
+    @classmethod
+    def from_xml_element(cls, elem, **kwargs):
+        filter_id = int(elem.get('id'))
+        values = [float(x) for x in get_text(elem, 'bins').split()]
+        particle = get_text(elem, 'particle')
+        return cls(values, filter_id=filter_id, particle=particle)
+
+
 class TimeFilter(RealFilter):
     """Bins tally events based on the particle's time.
 
@@ -1851,7 +1916,7 @@ def bins(self, bins):
         # Make sure there is only 1 bin.
         if not len(bins) == 1:
             msg = (f'Unable to add bins "{bins}" to a DistribcellFilter since '
-                  'only a single distribcell can be used per tally')
+                   'only a single distribcell can be used per tally')
             raise ValueError(msg)
 
         # Check the type and extract the id, if necessary.
@@ -2004,7 +2069,7 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs):
         # requests Summary geometric information
         filter_bins = _repeat_and_tile(
             np.arange(self.num_bins), stride, data_size)
-        df = pd.DataFrame({self.short_name.lower() : filter_bins})
+        df = pd.DataFrame({self.short_name.lower(): filter_bins})
 
         # Concatenate with DataFrame of distribcell instance IDs
         if level_df is not None:
@@ -2043,6 +2108,7 @@ class MuFilter(RealFilter):
         The number of filter bins
 
     """
+
     def __init__(self, values, filter_id=None):
         if isinstance(values, Integral):
             values = np.linspace(-1., 1., values + 1)
@@ -2208,6 +2274,7 @@ class DelayedGroupFilter(Filter):
         The number of filter bins
 
     """
+
     def check_bins(self, bins):
         # Check the bin values.
         for g in bins:
@@ -2276,9 +2343,9 @@ def __eq__(self, other):
     def __gt__(self, other):
         if type(self) is not type(other):
             if self.short_name in _FILTER_TYPES and \
-                other.short_name in _FILTER_TYPES:
+                    other.short_name in _FILTER_TYPES:
                 delta = _FILTER_TYPES.index(self.short_name) - \
-                        _FILTER_TYPES.index(other.short_name)
+                    _FILTER_TYPES.index(other.short_name)
                 return delta > 0
             else:
                 return False
@@ -2288,9 +2355,9 @@ def __gt__(self, other):
     def __lt__(self, other):
         if type(self) is not type(other):
             if self.short_name in _FILTER_TYPES and \
-                other.short_name in _FILTER_TYPES:
+                    other.short_name in _FILTER_TYPES:
                 delta = _FILTER_TYPES.index(self.short_name) - \
-                        _FILTER_TYPES.index(other.short_name)
+                    _FILTER_TYPES.index(other.short_name)
                 return delta < 0
             else:
                 return False
@@ -2301,14 +2368,16 @@ def __hash__(self):
         string = type(self).__name__ + '\n'
         string += '{: <16}=\t{}\n'.format('\tEnergy', self.energy)
         string += '{: <16}=\t{}\n'.format('\tInterpolant', self.y)
-        string += '{: <16}=\t{}\n'.format('\tInterpolation', self.interpolation)
+        string += '{: <16}=\t{}\n'.format('\tInterpolation',
+                                          self.interpolation)
         return hash(string)
 
     def __repr__(self):
         string = type(self).__name__ + '\n'
         string += '{: <16}=\t{}\n'.format('\tEnergy', self.energy)
         string += '{: <16}=\t{}\n'.format('\tInterpolant', self.y)
-        string += '{: <16}=\t{}\n'.format('\tInterpolation', self.interpolation)
+        string += '{: <16}=\t{}\n'.format('\tInterpolation',
+                                          self.interpolation)
         string += '{: <16}=\t{}\n'.format('\tID', self.id)
         return string
 
@@ -2394,10 +2463,12 @@ def interpolation(self):
     @interpolation.setter
     def interpolation(self, val):
         cv.check_type('interpolation', val, str)
-        cv.check_value('interpolation', val, self.INTERPOLATION_SCHEMES.values())
+        cv.check_value('interpolation', val,
+                       self.INTERPOLATION_SCHEMES.values())
 
         if val == 'quadratic' and len(self.energy) < 3:
-            raise ValueError('Quadratic interpolation requires 3 or more values.')
+            raise ValueError(
+                'Quadratic interpolation requires 3 or more values.')
 
         if val == 'cubic' and len(self.energy) < 4:
             raise ValueError('Cubic interpolation requires 3 or more values.')

src/particle.cpp

@@ -84,6 +84,11 @@ bool Particle::create_secondary(
     return false;
   }
 
+  // This is used to count backward in the secondary source bank for tallying
+  // outgoing photon energies from a neutron in the SecondaryPhotonEnergy
+  // filter.
+  secondaries_this_collision()++;
+
   auto& bank = secondary_bank().emplace_back();
   bank.particle = type;
   bank.wgt = wgt;
@@ -320,6 +325,7 @@ void Particle::event_cross_surface()
 
 void Particle::event_collide()
 {
+
   // Score collision estimate of keff
   if (settings::run_mode == RunMode::EIGENVALUE &&
       type() == ParticleType::neutron) {
@@ -364,6 +370,11 @@ void Particle::event_collide()
   n_bank() = 0;
   bank_second_E() = 0.0;
   wgt_bank() = 0.0;
+
+  // Clear number of secondaries in this collision. This is
+  // distinct from the number of created neutrons n_bank() above!
+  secondaries_this_collision() = 0;
+
   zero_delayed_bank();
 
   // Reset fission logical

src/tallies/filter.cpp

@@ -124,6 +124,8 @@ Filter* Filter::create(const std::string& type, int32_t id)
     return Filter::create<CollisionFilter>(id);
   } else if (type == "energyout") {
     return Filter::create<EnergyoutFilter>(id);
+  } else if (type == "secondaryenergy") {
+    return Filter::create<SecondaryEnergyFilter>(id);
   } else if (type == "legendre") {
     return Filter::create<LegendreFilter>(id);
   } else if (type == "material") {

src/tallies/filter_delayedgroup.cpp

@@ -27,7 +27,7 @@ void DelayedGroupFilter::set_groups(span<int> groups)
     } else if (group > MAX_DELAYED_GROUPS) {
       throw std::invalid_argument {
         "Encountered delayedgroup bin with index " + std::to_string(group) +
-        " which is greater than MAX_DELATED_GROUPS (" +
+        " which is greater than MAX_DELAYED_GROUPS (" +
         std::to_string(MAX_DELAYED_GROUPS) + ")"};
     }
     groups_.push_back(group);
@@ -39,6 +39,10 @@ void DelayedGroupFilter::set_groups(span<int> groups)
 void DelayedGroupFilter::get_all_bins(
   const Particle& p, TallyEstimator estimator, FilterMatch& match) const
 {
+  // Note that the bin is set to zero here, but bins outside zero are
+  // tallied to regardless. This is because that logic has to be handled
+  // in the scoring code instead where looping over the delayed
+  // group takes place (tally_scoring.cpp).
   match.bins_.push_back(0);
   match.weights_.push_back(1.0);
 }