diff --git a/CMakeLists.txt b/CMakeLists.txt index 48f3bfc398b..abed2307066 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index 2a9d0876cd2..a5a8de6d900 100644 --- a/docs/source/pythonapi/base.rst +++ b/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 diff --git a/docs/source/pythonapi/capi.rst b/docs/source/pythonapi/capi.rst index 67eca009471..0e9836c89d6 100644 --- a/docs/source/pythonapi/capi.rst +++ b/docs/source/pythonapi/capi.rst @@ -69,6 +69,7 @@ Classes EnergyFunctionFilter EnergyoutFilter Filter + FissionYieldsFilter LegendreFilter Material MaterialFilter diff --git a/include/openmc/chain.h b/include/openmc/chain.h index a3bc6f3a364..5b6a7f61ccb 100644 --- a/include/openmc/chain.h +++ b/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 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>& 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> 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 chain_nuclide_map; extern vector> chain_nuclides; +extern vector> fission_yields; } // namespace data diff --git a/include/openmc/endf.h b/include/openmc/endf.h index 4a737eb8816..ae79ee315da 100644 --- a/include/openmc/endf.h +++ b/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& x, vector& y); + //! Evaluate the tabulated function //! \param[in] x independent variable //! \return Function evaluated at x diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index 65098597a5d..713f45cc60a 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -28,6 +28,7 @@ enum class FilterType { ENERGY_FUNCTION, ENERGY, ENERGY_OUT, + FISSION_YIELDS, LEGENDRE, MATERIAL, MATERIALFROM, diff --git a/include/openmc/tallies/filter_fission_yields.h b/include/openmc/tallies/filter_fission_yields.h new file mode 100644 index 00000000000..ac58f6650fe --- /dev/null +++ b/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& bins() const { return bins_; } + +private: + //---------------------------------------------------------------------------- + // Data members + + vector bins_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_FISSION_ENERGY_H diff --git a/openmc/filter.py b/openmc/filter.py index f29c540857c..d1bc7b528a1 100644 --- a/openmc/filter.py +++ b/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. diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index 55b681d89ad..1cba6f812f4 100644 --- a/openmc/lib/filter.py +++ b/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, diff --git a/src/chain.cpp b/src/chain.cpp index e279d1f5916..418a9fdf931 100644 --- a/src/chain.cpp +++ b/src/chain.cpp @@ -6,6 +6,8 @@ #include // for getenv #include // for make_unique #include // for stod +#include +#include #include #include @@ -17,6 +19,34 @@ namespace openmc { +//============================================================================== +// FissionYields implementation +//============================================================================== + +FissionYields::FissionYields(pugi::xml_node node) +{ + std::unordered_map>> temp; + auto energies = get_node_array(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(fy, "products"); + auto data = get_node_array(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 x_; + vector 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(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 chain_nuclide_map; vector> chain_nuclides; +vector> fission_yields; } // namespace data diff --git a/src/endf.cpp b/src/endf.cpp index c0c1d2e7e8b..b01d160640e 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -164,6 +164,14 @@ Tabulated1D::Tabulated1D(hid_t dset) n_pairs_ = x_.size(); } +Tabulated1D::Tabulated1D(vector& x, vector& 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 diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index 79817981db0..81da354fe15 100644 --- a/src/tallies/filter.cpp +++ b/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(id); } else if (type == "energyout") { return Filter::create(id); + } else if (type == "fissionyields") { + return Filter::create(id); } else if (type == "legendre") { return Filter::create(id); } else if (type == "material") { diff --git a/src/tallies/filter_fission_yields.cpp b/src/tallies/filter_fission_yields.cpp new file mode 100644 index 00000000000..2c7b03bd5ae --- /dev/null +++ b/src/tallies/filter_fission_yields.cpp @@ -0,0 +1,56 @@ +#include "openmc/tallies/filter_fission_yields.h" + +#include + +#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(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 diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 54840b5188e..3022bcdd556 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -527,6 +527,7 @@ void Tally::set_scores(const vector& scores) // Check for the presence of certain restrictive filters. bool energyout_present = energyout_filter_ != C_NONE; + bool fission_yields_present = false; bool legendre_present = false; bool cell_present = false; bool cellfrom_present = false; @@ -550,6 +551,8 @@ void Tally::set_scores(const vector& scores) surface_present = true; } else if (filt->type() == FilterType::MESH_SURFACE) { meshsurface_present = true; + } else if (filt->type() == FilterType::FISSION_YIELDS) { + fission_yields_present = true; } } @@ -575,11 +578,18 @@ void Tally::set_scores(const vector& scores) case SCORE_TOTAL: case SCORE_ABSORPTION: + if (energyout_present) + fatal_error("Cannot tally " + score_str + + " reaction rate with an " + "outgoing energy filter"); + break; case SCORE_FISSION: if (energyout_present) fatal_error("Cannot tally " + score_str + " reaction rate with an " "outgoing energy filter"); + if (fission_yields_present && estimator_ == TallyEstimator::TRACKLENGTH) + estimator_ = TallyEstimator::COLLISION; break; case SCORE_SCATTER: diff --git a/tests/unit_tests/test_filter_fission_yields.py b/tests/unit_tests/test_filter_fission_yields.py new file mode 100644 index 00000000000..3e3dd177715 --- /dev/null +++ b/tests/unit_tests/test_filter_fission_yields.py @@ -0,0 +1,37 @@ +import openmc +import numpy as np + + +def test_fissionyieldsfilter(run_in_tmpdir): + steel = openmc.Material(name='U') + steel.set_density('g/cm3', 8.00) + steel.add_nuclide('U235', 1.0) + + sphere = openmc.Sphere(r=50.0, boundary_type='vacuum') + cell = openmc.Cell(region=-sphere, fill=steel) + model = openmc.Model() + model.geometry = openmc.Geometry([cell]) + model.settings.particles = 100 + model.settings.batches = 10 + + model.settings.source = openmc.IndependentSource( + energy=openmc.stats.delta_function(1e6), + ) + model.settings.run_mode = "eigenvalue" + + fissionyields_filter = openmc.FissionYieldsFilter(['Xe135','I135']) + + tally = openmc.Tally() + tally.filters = [fissionyields_filter] + tally.scores = ['fission'] + model.tallies = openmc.Tallies([tally]) + + model.export_to_model_xml() + + # Run OpenMC + model.run(apply_tally_results=True) + + assert np.all(tally.mean.reshape(2)>0) + + +