A new FissionYieldsFilter for tallying fission yields.

CMakeLists.txt

@@ -431,6 +431,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_distribcell.cpp
   src/tallies/filter_energy.cpp
   src/tallies/filter_energyfunc.cpp
+  src/tallies/filter_fission_yields.cpp
   src/tallies/filter_legendre.cpp
   src/tallies/filter_material.cpp
   src/tallies/filter_materialfrom.cpp

docs/source/pythonapi/base.rst

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

docs/source/pythonapi/capi.rst

@@ -69,6 +69,7 @@ Classes
    EnergyFunctionFilter
    EnergyoutFilter
    Filter
+   FissionYieldsFilter
    LegendreFilter
    Material
    MaterialFilter

include/openmc/chain.h

@@ -12,11 +12,25 @@
 
 #include "openmc/angle_energy.h" // for AngleEnergy
 #include "openmc/distribution.h" // for UPtrDist
-#include "openmc/memory.h"       // for unique_ptr
+#include "openmc/endf.h"
+#include "openmc/memory.h" // for unique_ptr
 #include "openmc/vector.h"
 
 namespace openmc {
 
+//==============================================================================
+// Fission Yield Data
+//==============================================================================
+
+class FissionYields {
+public:
+  // Constructors, destructors
+  FissionYields(pugi::xml_node node);
+
+  // Data members
+  std::unordered_map<std::string, Tabulated1D> yields_;
+};
+
 //==============================================================================
 // Data for a nuclide in the depletion chain
 //==============================================================================
@@ -36,6 +50,7 @@ class ChainNuclide {
   //! Compute the decay constant for the nuclide
   //! \return Decay constant in [1/s]
   double decay_constant() const { return std::log(2.0) / half_life_; }
+  FissionYields* fission_yields();
 
   const Distribution* photon_energy() const { return photon_energy_.get(); }
   const std::unordered_map<int, vector<Product>>& reaction_products() const
@@ -45,9 +60,11 @@ class ChainNuclide {
 
 private:
   // Data members
-  std::string name_;          //!< Name of nuclide
-  double half_life_ {0.0};    //!< Half-life in [s]
-  double decay_energy_ {0.0}; //!< Decay energy in [eV]
+  std::string name_;                        //!< Name of nuclide
+  double half_life_ {0.0};                  //!< Half-life in [s]
+  double decay_energy_ {0.0};               //!< Decay energy in [eV]
+  FissionYields* fission_yields_ = nullptr; //!< Fission yields
+  std::string fission_yields_parent_ {""};  //!< Fission yields parent
   std::unordered_map<int, vector<Product>>
     reaction_products_;    //!< Map of MT to reaction products
   UPtrDist photon_energy_; //!< Decay photon energy distribution
@@ -83,6 +100,7 @@ namespace data {
 
 extern std::unordered_map<std::string, int> chain_nuclide_map;
 extern vector<unique_ptr<ChainNuclide>> chain_nuclides;
+extern vector<unique_ptr<FissionYields>> fission_yields;
 
 } // namespace data
 

include/openmc/endf.h

@@ -78,6 +78,8 @@ class Tabulated1D : public Function1D {
   //! \param[in] dset Dataset containing tabulated data
   explicit Tabulated1D(hid_t dset);
 
+  explicit Tabulated1D(vector<double>& x, vector<double>& y);
+
   //! Evaluate the tabulated function
   //! \param[in] x independent variable
   //! \return Function evaluated at x

include/openmc/tallies/filter.h

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

include/openmc/tallies/filter_fission_yields.h

@@ -0,0 +1,47 @@
+#ifndef OPENMC_TALLIES_FILTER_FISSION_ENERGY_H
+#define OPENMC_TALLIES_FILTER_FISSION_ENERGY_H
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Multiplies tally scores by fission yields.
+//==============================================================================
+
+class FissionYieldsFilter : public Filter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  ~FissionYieldsFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "fissionyields"; }
+  FilterType type() const override { return FilterType::FISSION_YIELDS; }
+
+  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<std::string>& bins() const { return bins_; }
+
+private:
+  //----------------------------------------------------------------------------
+  // Data members
+
+  vector<std::string> bins_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_FISSION_ENERGY_H

openmc/filter.py

@@ -23,7 +23,7 @@
 _FILTER_TYPES = (
     'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy',
     'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
-    'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre',
+    'energyfunction', 'fissionyields', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre',
     'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance',
     'collision', 'time', 'parentnuclide', 'weight', 'meshborn', 'meshsurface',
     'meshmaterial',
@@ -2509,6 +2509,31 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs):
 
         return df
 
+class FissionYieldsFilter(ParticleFilter):
+    """Bins tally fission events based on fission yields
+
+    Parameters
+    ----------
+    bins : str, or iterable of str
+        Names of nuclides (e.g., 'Ni65')
+    filter_id : int
+        Unique identifier for the filter
+
+    Attributes
+    ----------
+    bins : iterable of str
+        Names of nuclides
+    id : int
+        Unique identifier for the filter
+    num_bins : Integral
+        The number of filter bins
+
+    """
+    @Filter.bins.setter
+    def bins(self, bins):
+        bins = np.atleast_1d(bins)
+        cv.check_iterable_type('filter bins', bins, str)
+        self._bins = bins
 
 class WeightFilter(RealFilter):
     """Bins tally events based on the incoming particle weight.

openmc/lib/filter.py

@@ -319,6 +319,8 @@ def _get_attr(self, cfunc):
         cfunc(self._index, n, array_p)
         return as_array(array_p, (n.value, ))
 
+class FissionYieldsFilter(Filter):
+    filter_type = 'fissionyields'
 
 class LegendreFilter(Filter):
     filter_type = 'legendre'
@@ -657,6 +659,7 @@ class ZernikeRadialFilter(ZernikeFilter):
     'energy': EnergyFilter,
     'energyout': EnergyoutFilter,
     'energyfunction': EnergyFunctionFilter,
+    'fissionyields': FissionYieldsFilter,
     'legendre': LegendreFilter,
     'material': MaterialFilter,
     'materialfrom': MaterialFromFilter,

src/chain.cpp

@@ -6,6 +6,8 @@
 #include <cstdlib> // for getenv
 #include <memory>  // for make_unique
 #include <string>  // for stod
+#include <unordered_map>
+#include <utility>
 
 #include <fmt/core.h>
 #include <pugixml.hpp>
@@ -17,6 +19,34 @@
 
 namespace openmc {
 
+//==============================================================================
+// FissionYields implementation
+//==============================================================================
+
+FissionYields::FissionYields(pugi::xml_node node)
+{
+  std::unordered_map<std::string, vector<std::pair<double, double>>> temp;
+  auto energies = get_node_array<double>(node, "energies");
+  for (pugi::xml_node fy : node.children("fission_yields")) {
+    double energy = std::stod(get_node_value(fy, "energy"));
+    auto products = get_node_array<std::string>(fy, "products");
+    auto data = get_node_array<double>(fy, "data");
+    for (int32_t i = 0; i < products.size(); ++i) {
+      temp.insert({products[i], {}});
+      temp[products[i]].push_back(std::make_pair(energy, data[i]));
+    }
+  }
+  for (const auto& pair : temp) {
+    vector<double> x_;
+    vector<double> y_;
+    for (const auto& v : pair.second) {
+      x_.push_back(v.first);
+      y_.push_back(v.second);
+    }
+    yields_[pair.first] = Tabulated1D(x_, y_);
+  }
+}
+
 //==============================================================================
 // ChainNuclide implementation
 //==============================================================================
@@ -54,6 +84,16 @@ ChainNuclide::ChainNuclide(pugi::xml_node node)
     }
   }
 
+  if (check_for_node(node, "neutron_fission_yields")) {
+    pugi::xml_node nfy = node.child("neutron_fission_yields");
+    if (check_for_node(nfy, "parent")) {
+      fission_yields_parent_ = get_node_value(nfy, "parent");
+    } else {
+      data::fission_yields.push_back(std::make_unique<FissionYields>(nfy));
+      fission_yields_ = data::fission_yields.back().get();
+    }
+  }
+
   // Set entry in mapping
   data::chain_nuclide_map[name_] = data::chain_nuclides.size();
 }
@@ -63,6 +103,15 @@ ChainNuclide::~ChainNuclide()
   data::chain_nuclide_map.erase(name_);
 }
 
+FissionYields* ChainNuclide::fission_yields()
+{
+  if (fission_yields_parent_.size() > 0) {
+    return data::chain_nuclides[data::chain_nuclide_map[fission_yields_parent_]]
+      ->fission_yields();
+  } else {
+    return fission_yields_;
+  }
+}
 //==============================================================================
 // DecayPhotonAngleEnergy implementation
 //==============================================================================
@@ -82,6 +131,7 @@ namespace data {
 
 std::unordered_map<std::string, int> chain_nuclide_map;
 vector<unique_ptr<ChainNuclide>> chain_nuclides;
+vector<unique_ptr<FissionYields>> fission_yields;
 
 } // namespace data
 

src/endf.cpp

@@ -164,6 +164,14 @@ Tabulated1D::Tabulated1D(hid_t dset)
   n_pairs_ = x_.size();
 }
 
+Tabulated1D::Tabulated1D(vector<double>& x, vector<double>& y)
+{
+  n_regions_ = 0;
+  std::copy(x.begin(), x.end(), std::back_inserter(x_));
+  std::copy(y.begin(), y.end(), std::back_inserter(y_));
+  n_pairs_ = x_.size();
+}
+
 double Tabulated1D::operator()(double x) const
 {
   // find which bin the abscissa is in -- if the abscissa is outside the

src/tallies/filter.cpp

@@ -20,6 +20,7 @@
 #include "openmc/tallies/filter_distribcell.h"
 #include "openmc/tallies/filter_energy.h"
 #include "openmc/tallies/filter_energyfunc.h"
+#include "openmc/tallies/filter_fission_yields.h"
 #include "openmc/tallies/filter_legendre.h"
 #include "openmc/tallies/filter_material.h"
 #include "openmc/tallies/filter_materialfrom.h"
@@ -124,6 +125,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 == "fissionyields") {
+    return Filter::create<FissionYieldsFilter>(id);
   } else if (type == "legendre") {
     return Filter::create<LegendreFilter>(id);
   } else if (type == "material") {

src/tallies/filter_fission_yields.cpp

@@ -0,0 +1,56 @@
+#include "openmc/tallies/filter_fission_yields.h"
+
+#include <pugixml.hpp>
+
+#include "openmc/chain.h"
+#include "openmc/endf.h"
+#include "openmc/error.h"
+#include "openmc/nuclide.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void FissionYieldsFilter::from_xml(pugi::xml_node node)
+{
+  if (!settings::run_CE)
+    fatal_error("FissionYieldsFilter filters are only supported for "
+                "continuous-energy transport calculations");
+
+  if (!check_for_node(node, "bins"))
+    fatal_error("Bins not specified for FissionYieldsFilter.");
+
+  bins_ = get_node_array<std::string>(node, "bins");
+  n_bins_ = bins_.size();
+}
+
+void FissionYieldsFilter::get_all_bins(
+  const Particle& p, TallyEstimator estimator, FilterMatch& match) const
+{
+  if (p.neutron_xs(p.event_nuclide()).fission > 0) {
+    auto nuc = data::nuclides[p.event_nuclide()]->name_;
+    if (data::chain_nuclide_map.find(nuc) != data::chain_nuclide_map.end()) {
+      auto fy = data::chain_nuclides[data::chain_nuclide_map[nuc]]
+                  ->fission_yields()
+                  ->yields_;
+      for (int i = 0; i < bins_.size(); ++i) {
+        if (fy.find(bins_[i]) != fy.end()) {
+          match.bins_.push_back(i);
+          match.weights_.push_back(fy[bins_[i]](p.E_last()));
+        }
+      }
+    }
+  }
+}
+
+void FissionYieldsFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "bins", bins_);
+}
+
+std::string FissionYieldsFilter::text_label(int bin) const
+{
+  return fmt::format("Fission Yield [{}]", bins_[bin]);
+}
+
+} // namespace openmc