From ecca4c55ed0dda38846fb28a348fb364078f9f0b Mon Sep 17 00:00:00 2001 From: GuySten Date: Sun, 13 Jul 2025 22:03:43 +0300 Subject: [PATCH 01/11] wip --- include/openmc/chain.h | 213 ++++++++++++++++++++----------------- src/chain.cpp | 236 +++++++++++++++++++++-------------------- 2 files changed, 236 insertions(+), 213 deletions(-) diff --git a/include/openmc/chain.h b/include/openmc/chain.h index a3bc6f3a364..bd827cb5939 100644 --- a/include/openmc/chain.h +++ b/include/openmc/chain.h @@ -1,97 +1,116 @@ -//! \file chain.h -//! \brief Depletion chain and associated information - -#ifndef OPENMC_CHAIN_H -#define OPENMC_CHAIN_H - -#include -#include -#include - -#include "pugixml.hpp" - -#include "openmc/angle_energy.h" // for AngleEnergy -#include "openmc/distribution.h" // for UPtrDist -#include "openmc/memory.h" // for unique_ptr -#include "openmc/vector.h" - -namespace openmc { - -//============================================================================== -// Data for a nuclide in the depletion chain -//============================================================================== - -class ChainNuclide { -public: - // Types - struct Product { - std::string name; //!< Reaction product name - double branching_ratio; //!< Branching ratio - }; - - // Constructors, destructors - ChainNuclide(pugi::xml_node node); - ~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_; } - - const Distribution* photon_energy() const { return photon_energy_.get(); } - const std::unordered_map>& reaction_products() const - { - return reaction_products_; - } - -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::unordered_map> - reaction_products_; //!< Map of MT to reaction products - UPtrDist photon_energy_; //!< Decay photon energy distribution -}; - -//============================================================================== -// Angle-energy distribution for decay photon -//============================================================================== - -class DecayPhotonAngleEnergy : public AngleEnergy { -public: - explicit DecayPhotonAngleEnergy(const Distribution* dist) - : photon_energy_(dist) - {} - - //! Sample distribution for an angle and energy - //! \param[in] E_in Incoming energy in [eV] - //! \param[out] E_out Outgoing energy in [eV] - //! \param[out] mu Outgoing cosine with respect to current direction - //! \param[inout] seed Pseudorandom seed pointer - void sample( - double E_in, double& E_out, double& mu, uint64_t* seed) const override; - -private: - const Distribution* photon_energy_; -}; - -//============================================================================== -// Global variables -//============================================================================== - -namespace data { - -extern std::unordered_map chain_nuclide_map; -extern vector> chain_nuclides; - -} // namespace data - -//============================================================================== -// Non-member functions -//============================================================================== - -void read_chain_file_xml(); - -} // namespace openmc - -#endif // OPENMC_CHAIN_H +//! \file chain.h +//! \brief Depletion chain and associated information + +#ifndef OPENMC_CHAIN_H +#define OPENMC_CHAIN_H + +#include +#include +#include + +#include "pugixml.hpp" + +#include "openmc/angle_energy.h" // for AngleEnergy +#include "openmc/distribution.h" // for UPtrDist +#include "openmc/memory.h" // for unique_ptr +#include "openmc/vector.h" + +namespace openmc { + +//============================================================================== +// Data for a nuclide in the depletion chain +//============================================================================== + +class ChainNuclide { +public: + // Types + struct Product { + std::string name; //!< Reaction product name + double branching_ratio; //!< Branching ratio + }; + + // Constructors, destructors + ChainNuclide(pugi::xml_node node); + ~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_; } + + const Distribution* photon_energy() const { return photon_energy_.get(); } + const std::unordered_map>& reaction_products() const + { + return reaction_products_; + } + +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] + bool fissionable_ {false}; //!< Can do fission + double fission_energy_ {0.0}; //!< Fission energy in [eV] + FissionYields* fission_yields_; //!< Fission yields + std::unordered_map> + reaction_products_; //!< Map of MT to reaction products + UPtrDist photon_energy_; //!< Decay photon energy distribution +}; + +//============================================================================== +// Angle-energy distribution for decay photon +//============================================================================== + +class DecayPhotonAngleEnergy : public AngleEnergy { +public: + explicit DecayPhotonAngleEnergy(const Distribution* dist) + : photon_energy_(dist) + {} + + //! Sample distribution for an angle and energy + //! \param[in] E_in Incoming energy in [eV] + //! \param[out] E_out Outgoing energy in [eV] + //! \param[out] mu Outgoing cosine with respect to current direction + //! \param[inout] seed Pseudorandom seed pointer + void sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const override; + +private: + const Distribution* photon_energy_; +}; + +//============================================================================== +// Fission Yield Data +//============================================================================== + +class FissionYields { +public: + // Constructors, destructors + FissionYields(pugi::xml_node node); + ~FissionYields(); + +private: + // Data members + std::unordered_map> yields_; +}; + +//============================================================================== +// Global variables +//============================================================================== + +namespace data { + +extern std::unordered_map chain_nuclide_map; +extern vector> chain_nuclides; +extern vector> fission_yields; + +} // namespace data + +//============================================================================== +// Non-member functions +//============================================================================== + +void read_chain_file_xml(); + +} // namespace openmc + +#endif // OPENMC_CHAIN_H diff --git a/src/chain.cpp b/src/chain.cpp index e279d1f5916..27639c4fb77 100644 --- a/src/chain.cpp +++ b/src/chain.cpp @@ -1,116 +1,120 @@ -//! \file chain.cpp -//! \brief Depletion chain and associated information - -#include "openmc/chain.h" - -#include // for getenv -#include // for make_unique -#include // for stod - -#include -#include - -#include "openmc/distribution.h" // for distribution_from_xml -#include "openmc/error.h" -#include "openmc/reaction.h" -#include "openmc/xml_interface.h" // for get_node_value - -namespace openmc { - -//============================================================================== -// ChainNuclide implementation -//============================================================================== - -ChainNuclide::ChainNuclide(pugi::xml_node node) -{ - name_ = get_node_value(node, "name"); - if (check_for_node(node, "half_life")) { - half_life_ = std::stod(get_node_value(node, "half_life")); - } - if (check_for_node(node, "decay_energy")) { - decay_energy_ = std::stod(get_node_value(node, "decay_energy")); - } - - // Read reactions to store MT -> product map - for (pugi::xml_node reaction_node : node.children("reaction")) { - std::string rx_name = get_node_value(reaction_node, "type"); - if (!reaction_node.attribute("target")) - continue; - std::string rx_target = get_node_value(reaction_node, "target"); - double branching_ratio = 1.0; - if (reaction_node.attribute("branching_ratio")) { - branching_ratio = - std::stod(get_node_value(reaction_node, "branching_ratio")); - } - int mt = reaction_type(rx_name); - reaction_products_[mt].push_back({rx_target, branching_ratio}); - } - - for (pugi::xml_node source_node : node.children("source")) { - auto particle = get_node_value(source_node, "particle"); - if (particle == "photon") { - photon_energy_ = distribution_from_xml(source_node); - break; - } - } - - // Set entry in mapping - data::chain_nuclide_map[name_] = data::chain_nuclides.size(); -} - -ChainNuclide::~ChainNuclide() -{ - data::chain_nuclide_map.erase(name_); -} - -//============================================================================== -// DecayPhotonAngleEnergy implementation -//============================================================================== - -void DecayPhotonAngleEnergy::sample( - double E_in, double& E_out, double& mu, uint64_t* seed) const -{ - E_out = photon_energy_->sample(seed); - mu = Uniform(-1., 1.).sample(seed); -} - -//============================================================================== -// Global variables -//============================================================================== - -namespace data { - -std::unordered_map chain_nuclide_map; -vector> chain_nuclides; - -} // namespace data - -//============================================================================== -// Non-member functions -//============================================================================== - -void read_chain_file_xml() -{ - char* chain_file_path = std::getenv("OPENMC_CHAIN_FILE"); - if (!chain_file_path) { - return; - } - - write_message(5, "Reading chain file: {}...", chain_file_path); - - pugi::xml_document doc; - auto result = doc.load_file(chain_file_path); - if (!result) { - fatal_error( - fmt::format("Error processing chain file: {}", chain_file_path)); - } - - // Get root element - pugi::xml_node root = doc.document_element(); - - for (auto node : root.children("nuclide")) { - data::chain_nuclides.push_back(std::make_unique(node)); - } -} - -} // namespace openmc +//! \file chain.cpp +//! \brief Depletion chain and associated information + +#include "openmc/chain.h" + +#include // for getenv +#include // for make_unique +#include // for stod + +#include +#include + +#include "openmc/distribution.h" // for distribution_from_xml +#include "openmc/error.h" +#include "openmc/reaction.h" +#include "openmc/xml_interface.h" // for get_node_value + +namespace openmc { + +//============================================================================== +// ChainNuclide implementation +//============================================================================== + +ChainNuclide::ChainNuclide(pugi::xml_node node) +{ + name_ = get_node_value(node, "name"); + if (check_for_node(node, "half_life")) { + half_life_ = std::stod(get_node_value(node, "half_life")); + } + if (check_for_node(node, "decay_energy")) { + decay_energy_ = std::stod(get_node_value(node, "decay_energy")); + } + + // Read reactions to store MT -> product map + for (pugi::xml_node reaction_node : node.children("reaction")) { + std::string rx_name = get_node_value(reaction_node, "type"); + if (rx_name == "fission") { + fission_energy_ = std::stod(get_node_value(reaction_node, "Q")); + fissionable_ = true; + } + if (!reaction_node.attribute("target")) + continue; + std::string rx_target = get_node_value(reaction_node, "target"); + double branching_ratio = 1.0; + if (reaction_node.attribute("branching_ratio")) { + branching_ratio = + std::stod(get_node_value(reaction_node, "branching_ratio")); + } + int mt = reaction_type(rx_name); + reaction_products_[mt].push_back({rx_target, branching_ratio}); + } + + for (pugi::xml_node source_node : node.children("source")) { + auto particle = get_node_value(source_node, "particle"); + if (particle == "photon") { + photon_energy_ = distribution_from_xml(source_node); + break; + } + } + + // Set entry in mapping + data::chain_nuclide_map[name_] = data::chain_nuclides.size(); +} + +ChainNuclide::~ChainNuclide() +{ + data::chain_nuclide_map.erase(name_); +} + +//============================================================================== +// DecayPhotonAngleEnergy implementation +//============================================================================== + +void DecayPhotonAngleEnergy::sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const +{ + E_out = photon_energy_->sample(seed); + mu = Uniform(-1., 1.).sample(seed); +} + +//============================================================================== +// Global variables +//============================================================================== + +namespace data { + +std::unordered_map chain_nuclide_map; +vector> chain_nuclides; + +} // namespace data + +//============================================================================== +// Non-member functions +//============================================================================== + +void read_chain_file_xml() +{ + char* chain_file_path = std::getenv("OPENMC_CHAIN_FILE"); + if (!chain_file_path) { + return; + } + + write_message(5, "Reading chain file: {}...", chain_file_path); + + pugi::xml_document doc; + auto result = doc.load_file(chain_file_path); + if (!result) { + fatal_error( + fmt::format("Error processing chain file: {}", chain_file_path)); + } + + // Get root element + pugi::xml_node root = doc.document_element(); + + for (auto node : root.children("nuclide")) { + data::chain_nuclides.push_back(std::make_unique(node)); + } +} + +} // namespace openmc From 0484799bd296f6b63aa94d6562349306fdd50824 Mon Sep 17 00:00:00 2001 From: GuySten Date: Mon, 14 Jul 2025 22:13:48 +0300 Subject: [PATCH 02/11] added new fission yields filter --- CMakeLists.txt | 1 + include/openmc/chain.h | 45 +- include/openmc/endf.h | 324 +++++----- include/openmc/tallies/filter.h | 1 + .../openmc/tallies/filter_fission_yields.h | 48 ++ openmc/filter.py | 27 +- src/chain.cpp | 50 ++ src/endf.cpp | 612 +++++++++--------- src/tallies/filter.cpp | 3 + src/tallies/filter_fission_yields.cpp | 56 ++ src/tallies/tally.cpp | 10 + 11 files changed, 691 insertions(+), 486 deletions(-) create mode 100644 include/openmc/tallies/filter_fission_yields.h create mode 100644 src/tallies/filter_fission_yields.cpp 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/include/openmc/chain.h b/include/openmc/chain.h index bd827cb5939..201fd88d1f1 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,12 +60,13 @@ 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] - bool fissionable_ {false}; //!< Can do fission - double fission_energy_ {0.0}; //!< Fission energy in [eV] - FissionYields* fission_yields_; //!< Fission yields + std::string name_; //!< Name of nuclide + double half_life_ {0.0}; //!< Half-life in [s] + double decay_energy_ {0.0}; //!< Decay energy in [eV] + bool fissionable_ {false}; //!< Can do fission + double fission_energy_ {0.0}; //!< Fission 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 @@ -78,21 +94,6 @@ class DecayPhotonAngleEnergy : public AngleEnergy { const Distribution* photon_energy_; }; -//============================================================================== -// Fission Yield Data -//============================================================================== - -class FissionYields { -public: - // Constructors, destructors - FissionYields(pugi::xml_node node); - ~FissionYields(); - -private: - // Data members - std::unordered_map> yields_; -}; - //============================================================================== // Global variables //============================================================================== diff --git a/include/openmc/endf.h b/include/openmc/endf.h index 4a737eb8816..ced89cff493 100644 --- a/include/openmc/endf.h +++ b/include/openmc/endf.h @@ -1,161 +1,163 @@ -//! \file endf.h -//! Classes and functions related to the ENDF-6 format - -#ifndef OPENMC_ENDF_H -#define OPENMC_ENDF_H - -#include "hdf5.h" - -#include "openmc/constants.h" -#include "openmc/memory.h" -#include "openmc/vector.h" - -namespace openmc { - -//! Convert integer representing interpolation law to enum -//! \param[in] i Intereger (e.g. 1=histogram, 2=lin-lin) -//! \return Corresponding enum value -Interpolation int2interp(int i); - -//! Determine whether MT number corresponds to a fission reaction -//! \param[in] MT ENDF MT value -//! \return Whether corresponding reaction is a fission reaction -bool is_fission(int MT); - -//! Determine if a given MT number is that of a disappearance reaction, i.e., a -//! reaction with no neutron in the exit channel -//! \param[in] MT ENDF MT value -//! \return Whether corresponding reaction is a disappearance reaction -bool is_disappearance(int MT); - -//! Determine if a given MT number is that of an inelastic scattering reaction -//! \param[in] MT ENDF MT value -//! \return Whether corresponding reaction is an inelastic scattering reaction -bool is_inelastic_scatter(int MT); - -//============================================================================== -//! Abstract one-dimensional function -//============================================================================== - -class Function1D { -public: - virtual double operator()(double x) const = 0; - virtual ~Function1D() = default; -}; - -//============================================================================== -//! One-dimensional function expressed as a polynomial -//============================================================================== - -class Polynomial : public Function1D { -public: - //! Construct polynomial from HDF5 data - //! \param[in] dset Dataset containing coefficients - explicit Polynomial(hid_t dset); - - //! Construct polynomial from coefficients - //! \param[in] coef Polynomial coefficients - explicit Polynomial(vector coef) : coef_(coef) {} - - //! Evaluate the polynomials - //! \param[in] x independent variable - //! \return Polynomial evaluated at x - double operator()(double x) const override; - -private: - vector coef_; //!< Polynomial coefficients -}; - -//============================================================================== -//! One-dimensional interpolable function -//============================================================================== - -class Tabulated1D : public Function1D { -public: - Tabulated1D() = default; - - //! Construct function from HDF5 data - //! \param[in] dset Dataset containing tabulated data - explicit Tabulated1D(hid_t dset); - - //! Evaluate the tabulated function - //! \param[in] x independent variable - //! \return Function evaluated at x - double operator()(double x) const override; - - // Accessors - const vector& x() const { return x_; } - const vector& y() const { return y_; } - -private: - std::size_t n_regions_ {0}; //!< number of interpolation regions - vector nbt_; //!< values separating interpolation regions - vector int_; //!< interpolation schemes - std::size_t n_pairs_; //!< number of (x,y) pairs - vector x_; //!< values of abscissa - vector y_; //!< values of ordinate -}; - -//============================================================================== -//! Coherent elastic scattering data from a crystalline material -//============================================================================== - -class CoherentElasticXS : public Function1D { -public: - explicit CoherentElasticXS(hid_t dset); - - double operator()(double E) const override; - - const vector& bragg_edges() const { return bragg_edges_; } - const vector& factors() const { return factors_; } - -private: - vector bragg_edges_; //!< Bragg edges in [eV] - vector factors_; //!< Partial sums of structure factors [eV-b] -}; - -//============================================================================== -//! Incoherent elastic scattering cross section -//============================================================================== - -class IncoherentElasticXS : public Function1D { -public: - explicit IncoherentElasticXS(hid_t dset); - - double operator()(double E) const override; - -private: - double bound_xs_; //!< Characteristic bound xs in [b] - double - debye_waller_; //!< Debye-Waller integral divided by atomic mass in [eV^-1] -}; - -//============================================================================== -//! Sum of multiple 1D functions -//============================================================================== - -class Sum1D : public Function1D { -public: - // Constructors - explicit Sum1D(hid_t group); - - //! Evaluate each function and sum results - //! \param[in] x independent variable - //! \return Function evaluated at x - double operator()(double E) const override; - - const unique_ptr& functions(int i) const { return functions_[i]; } - -private: - vector> functions_; //!< individual functions -}; - -//! Read 1D function from HDF5 dataset -//! \param[in] group HDF5 group containing dataset -//! \param[in] name Name of dataset -//! \return Unique pointer to 1D function -unique_ptr read_function(hid_t group, const char* name); - -} // namespace openmc - -#endif // OPENMC_ENDF_H +//! \file endf.h +//! Classes and functions related to the ENDF-6 format + +#ifndef OPENMC_ENDF_H +#define OPENMC_ENDF_H + +#include "hdf5.h" + +#include "openmc/constants.h" +#include "openmc/memory.h" +#include "openmc/vector.h" + +namespace openmc { + +//! Convert integer representing interpolation law to enum +//! \param[in] i Intereger (e.g. 1=histogram, 2=lin-lin) +//! \return Corresponding enum value +Interpolation int2interp(int i); + +//! Determine whether MT number corresponds to a fission reaction +//! \param[in] MT ENDF MT value +//! \return Whether corresponding reaction is a fission reaction +bool is_fission(int MT); + +//! Determine if a given MT number is that of a disappearance reaction, i.e., a +//! reaction with no neutron in the exit channel +//! \param[in] MT ENDF MT value +//! \return Whether corresponding reaction is a disappearance reaction +bool is_disappearance(int MT); + +//! Determine if a given MT number is that of an inelastic scattering reaction +//! \param[in] MT ENDF MT value +//! \return Whether corresponding reaction is an inelastic scattering reaction +bool is_inelastic_scatter(int MT); + +//============================================================================== +//! Abstract one-dimensional function +//============================================================================== + +class Function1D { +public: + virtual double operator()(double x) const = 0; + virtual ~Function1D() = default; +}; + +//============================================================================== +//! One-dimensional function expressed as a polynomial +//============================================================================== + +class Polynomial : public Function1D { +public: + //! Construct polynomial from HDF5 data + //! \param[in] dset Dataset containing coefficients + explicit Polynomial(hid_t dset); + + //! Construct polynomial from coefficients + //! \param[in] coef Polynomial coefficients + explicit Polynomial(vector coef) : coef_(coef) {} + + //! Evaluate the polynomials + //! \param[in] x independent variable + //! \return Polynomial evaluated at x + double operator()(double x) const override; + +private: + vector coef_; //!< Polynomial coefficients +}; + +//============================================================================== +//! One-dimensional interpolable function +//============================================================================== + +class Tabulated1D : public Function1D { +public: + Tabulated1D() = default; + + //! Construct function from HDF5 data + //! \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 + double operator()(double x) const override; + + // Accessors + const vector& x() const { return x_; } + const vector& y() const { return y_; } + +private: + std::size_t n_regions_ {0}; //!< number of interpolation regions + vector nbt_; //!< values separating interpolation regions + vector int_; //!< interpolation schemes + std::size_t n_pairs_; //!< number of (x,y) pairs + vector x_; //!< values of abscissa + vector y_; //!< values of ordinate +}; + +//============================================================================== +//! Coherent elastic scattering data from a crystalline material +//============================================================================== + +class CoherentElasticXS : public Function1D { +public: + explicit CoherentElasticXS(hid_t dset); + + double operator()(double E) const override; + + const vector& bragg_edges() const { return bragg_edges_; } + const vector& factors() const { return factors_; } + +private: + vector bragg_edges_; //!< Bragg edges in [eV] + vector factors_; //!< Partial sums of structure factors [eV-b] +}; + +//============================================================================== +//! Incoherent elastic scattering cross section +//============================================================================== + +class IncoherentElasticXS : public Function1D { +public: + explicit IncoherentElasticXS(hid_t dset); + + double operator()(double E) const override; + +private: + double bound_xs_; //!< Characteristic bound xs in [b] + double + debye_waller_; //!< Debye-Waller integral divided by atomic mass in [eV^-1] +}; + +//============================================================================== +//! Sum of multiple 1D functions +//============================================================================== + +class Sum1D : public Function1D { +public: + // Constructors + explicit Sum1D(hid_t group); + + //! Evaluate each function and sum results + //! \param[in] x independent variable + //! \return Function evaluated at x + double operator()(double E) const override; + + const unique_ptr& functions(int i) const { return functions_[i]; } + +private: + vector> functions_; //!< individual functions +}; + +//! Read 1D function from HDF5 dataset +//! \param[in] group HDF5 group containing dataset +//! \param[in] name Name of dataset +//! \return Unique pointer to 1D function +unique_ptr read_function(hid_t group, const char* name); + +} // namespace openmc + +#endif // OPENMC_ENDF_H 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..cddef070e38 --- /dev/null +++ b/include/openmc/tallies/filter_fission_yields.h @@ -0,0 +1,48 @@ +#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 an arbitrary function of incident energy +//! described by a piecewise linear-linear interpolation. +//============================================================================== + +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& nuclides() const { return nuclides_; } + +private: + //---------------------------------------------------------------------------- + // Data members + + vector nuclides_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_FISSION_ENERGY_H diff --git a/openmc/filter.py b/openmc/filter.py index f29c540857c..11aeed9be1a 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', 'fission_yields', '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/src/chain.cpp b/src/chain.cpp index 27639c4fb77..e892c4593ca 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 //============================================================================== @@ -58,6 +88,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(); } @@ -67,6 +107,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 //============================================================================== @@ -86,6 +135,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..5608e195cbf 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -1,302 +1,310 @@ -#include "openmc/endf.h" - -#include // for copy -#include // for log, exp -#include // for back_inserter -#include // for runtime_error - -#include "xtensor/xarray.hpp" -#include "xtensor/xview.hpp" - -#include "openmc/array.h" -#include "openmc/constants.h" -#include "openmc/hdf5_interface.h" -#include "openmc/search.h" - -namespace openmc { - -//============================================================================== -// Functions -//============================================================================== - -Interpolation int2interp(int i) -{ - // TODO: We are ignoring specification of two-dimensional interpolation - // schemes (method of corresponding points and unit base interpolation). Those - // should be accounted for in the distribution classes somehow. - - switch (i) { - case 1: - case 11: - case 21: - return Interpolation::histogram; - case 2: - case 12: - case 22: - return Interpolation::lin_lin; - case 3: - case 13: - case 23: - return Interpolation::lin_log; - case 4: - case 14: - case 24: - return Interpolation::log_lin; - case 5: - case 15: - case 25: - return Interpolation::log_log; - default: - throw std::runtime_error {"Invalid interpolation code."}; - } -} - -bool is_fission(int mt) -{ - return mt == N_FISSION || mt == N_F || mt == N_NF || mt == N_2NF || - mt == N_3NF; -} - -bool is_disappearance(int mt) -{ - if (mt >= N_DISAPPEAR && mt <= N_DA) { - return true; - } else if (mt >= N_P0 && mt <= N_AC) { - return true; - } else if (mt == N_TA || mt == N_DT || mt == N_P3HE || mt == N_D3HE || - mt == N_3HEA || mt == N_3P) { - return true; - } else { - return false; - } -} - -bool is_inelastic_scatter(int mt) -{ - if (mt < 100) { - if (is_fission(mt)) { - return false; - } else { - return mt >= MISC && mt != 27; - } - } else if (mt <= 200) { - return !is_disappearance(mt); - } else if (mt >= N_2N0 && mt <= N_2NC) { - return true; - } else { - return false; - } -} - -unique_ptr read_function(hid_t group, const char* name) -{ - hid_t obj_id = open_object(group, name); - std::string func_type; - read_attribute(obj_id, "type", func_type); - unique_ptr func; - if (func_type == "Tabulated1D") { - func = make_unique(obj_id); - } else if (func_type == "Polynomial") { - func = make_unique(obj_id); - } else if (func_type == "CoherentElastic") { - func = make_unique(obj_id); - } else if (func_type == "IncoherentElastic") { - func = make_unique(obj_id); - } else if (func_type == "Sum") { - func = make_unique(obj_id); - } else { - throw std::runtime_error {"Unknown function type " + func_type + - " for dataset " + object_name(obj_id)}; - } - close_object(obj_id); - return func; -} - -//============================================================================== -// Polynomial implementation -//============================================================================== - -Polynomial::Polynomial(hid_t dset) -{ - // Read coefficients into a vector - read_dataset(dset, coef_); -} - -double Polynomial::operator()(double x) const -{ - // Use Horner's rule to evaluate polynomial. Note that coefficients are - // ordered in increasing powers of x. - double y = 0.0; - for (auto c = coef_.crbegin(); c != coef_.crend(); ++c) { - y = y * x + *c; - } - return y; -} - -//============================================================================== -// Tabulated1D implementation -//============================================================================== - -Tabulated1D::Tabulated1D(hid_t dset) -{ - read_attribute(dset, "breakpoints", nbt_); - n_regions_ = nbt_.size(); - - // Change 1-indexing to 0-indexing - for (auto& b : nbt_) - --b; - - vector int_temp; - read_attribute(dset, "interpolation", int_temp); - - // Convert vector of ints into Interpolation - for (const auto i : int_temp) - int_.push_back(int2interp(i)); - - xt::xarray arr; - read_dataset(dset, arr); - - auto xs = xt::view(arr, 0); - auto ys = xt::view(arr, 1); - - std::copy(xs.begin(), xs.end(), std::back_inserter(x_)); - std::copy(ys.begin(), ys.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 - // tabulated range, the first or last point is chosen, i.e. no interpolation - // is done outside the energy range - int i; - if (x < x_[0]) { - return y_[0]; - } else if (x > x_[n_pairs_ - 1]) { - return y_[n_pairs_ - 1]; - } else { - i = lower_bound_index(x_.begin(), x_.end(), x); - } - - // determine interpolation scheme - Interpolation interp; - if (n_regions_ == 0) { - interp = Interpolation::lin_lin; - } else { - interp = int_[0]; - for (int j = 0; j < n_regions_; ++j) { - if (i < nbt_[j]) { - interp = int_[j]; - break; - } - } - } - - // handle special case of histogram interpolation - if (interp == Interpolation::histogram) - return y_[i]; - - // determine bounding values - double x0 = x_[i]; - double x1 = x_[i + 1]; - double y0 = y_[i]; - double y1 = y_[i + 1]; - - // determine interpolation factor and interpolated value - double r; - switch (interp) { - case Interpolation::lin_lin: - r = (x - x0) / (x1 - x0); - return y0 + r * (y1 - y0); - case Interpolation::lin_log: - r = log(x / x0) / log(x1 / x0); - return y0 + r * (y1 - y0); - case Interpolation::log_lin: - r = (x - x0) / (x1 - x0); - return y0 * exp(r * log(y1 / y0)); - case Interpolation::log_log: - r = log(x / x0) / log(x1 / x0); - return y0 * exp(r * log(y1 / y0)); - default: - throw std::runtime_error {"Invalid interpolation scheme."}; - } -} - -//============================================================================== -// CoherentElasticXS implementation -//============================================================================== - -CoherentElasticXS::CoherentElasticXS(hid_t dset) -{ - // Read 2D array from dataset - xt::xarray arr; - read_dataset(dset, arr); - - // Get views for Bragg edges and structure factors - auto E = xt::view(arr, 0); - auto s = xt::view(arr, 1); - - // Copy Bragg edges and partial sums of structure factors - std::copy(E.begin(), E.end(), std::back_inserter(bragg_edges_)); - std::copy(s.begin(), s.end(), std::back_inserter(factors_)); -} - -double CoherentElasticXS::operator()(double E) const -{ - if (E < bragg_edges_[0]) { - // If energy is below that of the lowest Bragg peak, the elastic cross - // section will be zero - return 0.0; - } else { - auto i_grid = - lower_bound_index(bragg_edges_.begin(), bragg_edges_.end(), E); - return factors_[i_grid] / E; - } -} - -//============================================================================== -// IncoherentElasticXS implementation -//============================================================================== - -IncoherentElasticXS::IncoherentElasticXS(hid_t dset) -{ - array tmp; - read_dataset(dset, nullptr, tmp); - bound_xs_ = tmp[0]; - debye_waller_ = tmp[1]; -} - -double IncoherentElasticXS::operator()(double E) const -{ - // Determine cross section using ENDF-102, Eq. (7.5) - double W = debye_waller_; - return bound_xs_ / 2.0 * ((1 - std::exp(-4.0 * E * W)) / (2.0 * E * W)); -} - -//============================================================================== -// Sum1D implementation -//============================================================================== - -Sum1D::Sum1D(hid_t group) -{ - // Get number of functions - int n; - read_attribute(group, "n", n); - - // Get each function - for (int i = 0; i < n; ++i) { - auto dset_name = fmt::format("func_{}", i + 1); - functions_.push_back(read_function(group, dset_name.c_str())); - } -} - -double Sum1D::operator()(double x) const -{ - double result = 0.0; - for (auto& func : functions_) { - result += (*func)(x); - } - return result; -} - -} // namespace openmc +#include "openmc/endf.h" + +#include // for copy +#include // for log, exp +#include // for back_inserter +#include // for runtime_error + +#include "xtensor/xarray.hpp" +#include "xtensor/xview.hpp" + +#include "openmc/array.h" +#include "openmc/constants.h" +#include "openmc/hdf5_interface.h" +#include "openmc/search.h" + +namespace openmc { + +//============================================================================== +// Functions +//============================================================================== + +Interpolation int2interp(int i) +{ + // TODO: We are ignoring specification of two-dimensional interpolation + // schemes (method of corresponding points and unit base interpolation). Those + // should be accounted for in the distribution classes somehow. + + switch (i) { + case 1: + case 11: + case 21: + return Interpolation::histogram; + case 2: + case 12: + case 22: + return Interpolation::lin_lin; + case 3: + case 13: + case 23: + return Interpolation::lin_log; + case 4: + case 14: + case 24: + return Interpolation::log_lin; + case 5: + case 15: + case 25: + return Interpolation::log_log; + default: + throw std::runtime_error {"Invalid interpolation code."}; + } +} + +bool is_fission(int mt) +{ + return mt == N_FISSION || mt == N_F || mt == N_NF || mt == N_2NF || + mt == N_3NF; +} + +bool is_disappearance(int mt) +{ + if (mt >= N_DISAPPEAR && mt <= N_DA) { + return true; + } else if (mt >= N_P0 && mt <= N_AC) { + return true; + } else if (mt == N_TA || mt == N_DT || mt == N_P3HE || mt == N_D3HE || + mt == N_3HEA || mt == N_3P) { + return true; + } else { + return false; + } +} + +bool is_inelastic_scatter(int mt) +{ + if (mt < 100) { + if (is_fission(mt)) { + return false; + } else { + return mt >= MISC && mt != 27; + } + } else if (mt <= 200) { + return !is_disappearance(mt); + } else if (mt >= N_2N0 && mt <= N_2NC) { + return true; + } else { + return false; + } +} + +unique_ptr read_function(hid_t group, const char* name) +{ + hid_t obj_id = open_object(group, name); + std::string func_type; + read_attribute(obj_id, "type", func_type); + unique_ptr func; + if (func_type == "Tabulated1D") { + func = make_unique(obj_id); + } else if (func_type == "Polynomial") { + func = make_unique(obj_id); + } else if (func_type == "CoherentElastic") { + func = make_unique(obj_id); + } else if (func_type == "IncoherentElastic") { + func = make_unique(obj_id); + } else if (func_type == "Sum") { + func = make_unique(obj_id); + } else { + throw std::runtime_error {"Unknown function type " + func_type + + " for dataset " + object_name(obj_id)}; + } + close_object(obj_id); + return func; +} + +//============================================================================== +// Polynomial implementation +//============================================================================== + +Polynomial::Polynomial(hid_t dset) +{ + // Read coefficients into a vector + read_dataset(dset, coef_); +} + +double Polynomial::operator()(double x) const +{ + // Use Horner's rule to evaluate polynomial. Note that coefficients are + // ordered in increasing powers of x. + double y = 0.0; + for (auto c = coef_.crbegin(); c != coef_.crend(); ++c) { + y = y * x + *c; + } + return y; +} + +//============================================================================== +// Tabulated1D implementation +//============================================================================== + +Tabulated1D::Tabulated1D(hid_t dset) +{ + read_attribute(dset, "breakpoints", nbt_); + n_regions_ = nbt_.size(); + + // Change 1-indexing to 0-indexing + for (auto& b : nbt_) + --b; + + vector int_temp; + read_attribute(dset, "interpolation", int_temp); + + // Convert vector of ints into Interpolation + for (const auto i : int_temp) + int_.push_back(int2interp(i)); + + xt::xarray arr; + read_dataset(dset, arr); + + auto xs = xt::view(arr, 0); + auto ys = xt::view(arr, 1); + + std::copy(xs.begin(), xs.end(), std::back_inserter(x_)); + std::copy(ys.begin(), ys.end(), std::back_inserter(y_)); + 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 + // tabulated range, the first or last point is chosen, i.e. no interpolation + // is done outside the energy range + int i; + if (x < x_[0]) { + return y_[0]; + } else if (x > x_[n_pairs_ - 1]) { + return y_[n_pairs_ - 1]; + } else { + i = lower_bound_index(x_.begin(), x_.end(), x); + } + + // determine interpolation scheme + Interpolation interp; + if (n_regions_ == 0) { + interp = Interpolation::lin_lin; + } else { + interp = int_[0]; + for (int j = 0; j < n_regions_; ++j) { + if (i < nbt_[j]) { + interp = int_[j]; + break; + } + } + } + + // handle special case of histogram interpolation + if (interp == Interpolation::histogram) + return y_[i]; + + // determine bounding values + double x0 = x_[i]; + double x1 = x_[i + 1]; + double y0 = y_[i]; + double y1 = y_[i + 1]; + + // determine interpolation factor and interpolated value + double r; + switch (interp) { + case Interpolation::lin_lin: + r = (x - x0) / (x1 - x0); + return y0 + r * (y1 - y0); + case Interpolation::lin_log: + r = log(x / x0) / log(x1 / x0); + return y0 + r * (y1 - y0); + case Interpolation::log_lin: + r = (x - x0) / (x1 - x0); + return y0 * exp(r * log(y1 / y0)); + case Interpolation::log_log: + r = log(x / x0) / log(x1 / x0); + return y0 * exp(r * log(y1 / y0)); + default: + throw std::runtime_error {"Invalid interpolation scheme."}; + } +} + +//============================================================================== +// CoherentElasticXS implementation +//============================================================================== + +CoherentElasticXS::CoherentElasticXS(hid_t dset) +{ + // Read 2D array from dataset + xt::xarray arr; + read_dataset(dset, arr); + + // Get views for Bragg edges and structure factors + auto E = xt::view(arr, 0); + auto s = xt::view(arr, 1); + + // Copy Bragg edges and partial sums of structure factors + std::copy(E.begin(), E.end(), std::back_inserter(bragg_edges_)); + std::copy(s.begin(), s.end(), std::back_inserter(factors_)); +} + +double CoherentElasticXS::operator()(double E) const +{ + if (E < bragg_edges_[0]) { + // If energy is below that of the lowest Bragg peak, the elastic cross + // section will be zero + return 0.0; + } else { + auto i_grid = + lower_bound_index(bragg_edges_.begin(), bragg_edges_.end(), E); + return factors_[i_grid] / E; + } +} + +//============================================================================== +// IncoherentElasticXS implementation +//============================================================================== + +IncoherentElasticXS::IncoherentElasticXS(hid_t dset) +{ + array tmp; + read_dataset(dset, nullptr, tmp); + bound_xs_ = tmp[0]; + debye_waller_ = tmp[1]; +} + +double IncoherentElasticXS::operator()(double E) const +{ + // Determine cross section using ENDF-102, Eq. (7.5) + double W = debye_waller_; + return bound_xs_ / 2.0 * ((1 - std::exp(-4.0 * E * W)) / (2.0 * E * W)); +} + +//============================================================================== +// Sum1D implementation +//============================================================================== + +Sum1D::Sum1D(hid_t group) +{ + // Get number of functions + int n; + read_attribute(group, "n", n); + + // Get each function + for (int i = 0; i < n; ++i) { + auto dset_name = fmt::format("func_{}", i + 1); + functions_.push_back(read_function(group, dset_name.c_str())); + } +} + +double Sum1D::operator()(double x) const +{ + double result = 0.0; + for (auto& func : functions_) { + result += (*func)(x); + } + return result; +} + +} // namespace openmc 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..cbbff09af91 --- /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("Nuclides not specified for FissionYieldsFilter."); + + nuclides_ = get_node_array(node, "bins"); + n_bins_ = nuclides_.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 < nuclides_.size(); ++i) { + if (fy.find(nuclides_[i]) != fy.end()) { + match.bins_.push_back(i); + match.weights_.push_back(fy[nuclides_[i]](p.E_last())); + } + } + } + } +} + +void FissionYieldsFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + write_dataset(filter_group, "nuclides", nuclides_); +} + +std::string FissionYieldsFilter::text_label(int bin) const +{ + return fmt::format("Fission Yield [{}]", nuclides_[bin]); +} + +} // namespace openmc diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 54840b5188e..94a1ebfa924 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::ANALOG; break; case SCORE_SCATTER: From 3370a2cb8be0cf70d47627f85cbd9987894a779a Mon Sep 17 00:00:00 2001 From: GuySten Date: Mon, 14 Jul 2025 22:30:33 +0300 Subject: [PATCH 03/11] added test --- .../unit_tests/test_filter_fission_yields.py | 38 +++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 tests/unit_tests/test_filter_fission_yields.py 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..2bc4b41f534 --- /dev/null +++ b/tests/unit_tests/test_filter_fission_yields.py @@ -0,0 +1,38 @@ +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','Te135','I135']) + + tally = openmc.Tally() + tally.filters = [fissionyields_filter] + tally.estimator = 'analog' + 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(3)>0) + + + From ff597b7fec0f9c6fa346af276c9507269cf42958 Mon Sep 17 00:00:00 2001 From: GuySten Date: Mon, 14 Jul 2025 22:31:50 +0300 Subject: [PATCH 04/11] dos2unix --- include/openmc/chain.h | 234 ++++++++-------- include/openmc/endf.h | 326 +++++++++++----------- src/chain.cpp | 340 +++++++++++----------- src/endf.cpp | 620 ++++++++++++++++++++--------------------- 4 files changed, 760 insertions(+), 760 deletions(-) diff --git a/include/openmc/chain.h b/include/openmc/chain.h index 201fd88d1f1..12353fcb6d0 100644 --- a/include/openmc/chain.h +++ b/include/openmc/chain.h @@ -1,117 +1,117 @@ -//! \file chain.h -//! \brief Depletion chain and associated information - -#ifndef OPENMC_CHAIN_H -#define OPENMC_CHAIN_H - -#include -#include -#include - -#include "pugixml.hpp" - -#include "openmc/angle_energy.h" // for AngleEnergy -#include "openmc/distribution.h" // for UPtrDist -#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 -//============================================================================== - -class ChainNuclide { -public: - // Types - struct Product { - std::string name; //!< Reaction product name - double branching_ratio; //!< Branching ratio - }; - - // Constructors, destructors - ChainNuclide(pugi::xml_node node); - ~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 - { - return reaction_products_; - } - -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] - bool fissionable_ {false}; //!< Can do fission - double fission_energy_ {0.0}; //!< Fission 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 -}; - -//============================================================================== -// Angle-energy distribution for decay photon -//============================================================================== - -class DecayPhotonAngleEnergy : public AngleEnergy { -public: - explicit DecayPhotonAngleEnergy(const Distribution* dist) - : photon_energy_(dist) - {} - - //! Sample distribution for an angle and energy - //! \param[in] E_in Incoming energy in [eV] - //! \param[out] E_out Outgoing energy in [eV] - //! \param[out] mu Outgoing cosine with respect to current direction - //! \param[inout] seed Pseudorandom seed pointer - void sample( - double E_in, double& E_out, double& mu, uint64_t* seed) const override; - -private: - const Distribution* photon_energy_; -}; - -//============================================================================== -// Global variables -//============================================================================== - -namespace data { - -extern std::unordered_map chain_nuclide_map; -extern vector> chain_nuclides; -extern vector> fission_yields; - -} // namespace data - -//============================================================================== -// Non-member functions -//============================================================================== - -void read_chain_file_xml(); - -} // namespace openmc - -#endif // OPENMC_CHAIN_H +//! \file chain.h +//! \brief Depletion chain and associated information + +#ifndef OPENMC_CHAIN_H +#define OPENMC_CHAIN_H + +#include +#include +#include + +#include "pugixml.hpp" + +#include "openmc/angle_energy.h" // for AngleEnergy +#include "openmc/distribution.h" // for UPtrDist +#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 +//============================================================================== + +class ChainNuclide { +public: + // Types + struct Product { + std::string name; //!< Reaction product name + double branching_ratio; //!< Branching ratio + }; + + // Constructors, destructors + ChainNuclide(pugi::xml_node node); + ~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 + { + return reaction_products_; + } + +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] + bool fissionable_ {false}; //!< Can do fission + double fission_energy_ {0.0}; //!< Fission 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 +}; + +//============================================================================== +// Angle-energy distribution for decay photon +//============================================================================== + +class DecayPhotonAngleEnergy : public AngleEnergy { +public: + explicit DecayPhotonAngleEnergy(const Distribution* dist) + : photon_energy_(dist) + {} + + //! Sample distribution for an angle and energy + //! \param[in] E_in Incoming energy in [eV] + //! \param[out] E_out Outgoing energy in [eV] + //! \param[out] mu Outgoing cosine with respect to current direction + //! \param[inout] seed Pseudorandom seed pointer + void sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const override; + +private: + const Distribution* photon_energy_; +}; + +//============================================================================== +// Global variables +//============================================================================== + +namespace data { + +extern std::unordered_map chain_nuclide_map; +extern vector> chain_nuclides; +extern vector> fission_yields; + +} // namespace data + +//============================================================================== +// Non-member functions +//============================================================================== + +void read_chain_file_xml(); + +} // namespace openmc + +#endif // OPENMC_CHAIN_H diff --git a/include/openmc/endf.h b/include/openmc/endf.h index ced89cff493..ae79ee315da 100644 --- a/include/openmc/endf.h +++ b/include/openmc/endf.h @@ -1,163 +1,163 @@ -//! \file endf.h -//! Classes and functions related to the ENDF-6 format - -#ifndef OPENMC_ENDF_H -#define OPENMC_ENDF_H - -#include "hdf5.h" - -#include "openmc/constants.h" -#include "openmc/memory.h" -#include "openmc/vector.h" - -namespace openmc { - -//! Convert integer representing interpolation law to enum -//! \param[in] i Intereger (e.g. 1=histogram, 2=lin-lin) -//! \return Corresponding enum value -Interpolation int2interp(int i); - -//! Determine whether MT number corresponds to a fission reaction -//! \param[in] MT ENDF MT value -//! \return Whether corresponding reaction is a fission reaction -bool is_fission(int MT); - -//! Determine if a given MT number is that of a disappearance reaction, i.e., a -//! reaction with no neutron in the exit channel -//! \param[in] MT ENDF MT value -//! \return Whether corresponding reaction is a disappearance reaction -bool is_disappearance(int MT); - -//! Determine if a given MT number is that of an inelastic scattering reaction -//! \param[in] MT ENDF MT value -//! \return Whether corresponding reaction is an inelastic scattering reaction -bool is_inelastic_scatter(int MT); - -//============================================================================== -//! Abstract one-dimensional function -//============================================================================== - -class Function1D { -public: - virtual double operator()(double x) const = 0; - virtual ~Function1D() = default; -}; - -//============================================================================== -//! One-dimensional function expressed as a polynomial -//============================================================================== - -class Polynomial : public Function1D { -public: - //! Construct polynomial from HDF5 data - //! \param[in] dset Dataset containing coefficients - explicit Polynomial(hid_t dset); - - //! Construct polynomial from coefficients - //! \param[in] coef Polynomial coefficients - explicit Polynomial(vector coef) : coef_(coef) {} - - //! Evaluate the polynomials - //! \param[in] x independent variable - //! \return Polynomial evaluated at x - double operator()(double x) const override; - -private: - vector coef_; //!< Polynomial coefficients -}; - -//============================================================================== -//! One-dimensional interpolable function -//============================================================================== - -class Tabulated1D : public Function1D { -public: - Tabulated1D() = default; - - //! Construct function from HDF5 data - //! \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 - double operator()(double x) const override; - - // Accessors - const vector& x() const { return x_; } - const vector& y() const { return y_; } - -private: - std::size_t n_regions_ {0}; //!< number of interpolation regions - vector nbt_; //!< values separating interpolation regions - vector int_; //!< interpolation schemes - std::size_t n_pairs_; //!< number of (x,y) pairs - vector x_; //!< values of abscissa - vector y_; //!< values of ordinate -}; - -//============================================================================== -//! Coherent elastic scattering data from a crystalline material -//============================================================================== - -class CoherentElasticXS : public Function1D { -public: - explicit CoherentElasticXS(hid_t dset); - - double operator()(double E) const override; - - const vector& bragg_edges() const { return bragg_edges_; } - const vector& factors() const { return factors_; } - -private: - vector bragg_edges_; //!< Bragg edges in [eV] - vector factors_; //!< Partial sums of structure factors [eV-b] -}; - -//============================================================================== -//! Incoherent elastic scattering cross section -//============================================================================== - -class IncoherentElasticXS : public Function1D { -public: - explicit IncoherentElasticXS(hid_t dset); - - double operator()(double E) const override; - -private: - double bound_xs_; //!< Characteristic bound xs in [b] - double - debye_waller_; //!< Debye-Waller integral divided by atomic mass in [eV^-1] -}; - -//============================================================================== -//! Sum of multiple 1D functions -//============================================================================== - -class Sum1D : public Function1D { -public: - // Constructors - explicit Sum1D(hid_t group); - - //! Evaluate each function and sum results - //! \param[in] x independent variable - //! \return Function evaluated at x - double operator()(double E) const override; - - const unique_ptr& functions(int i) const { return functions_[i]; } - -private: - vector> functions_; //!< individual functions -}; - -//! Read 1D function from HDF5 dataset -//! \param[in] group HDF5 group containing dataset -//! \param[in] name Name of dataset -//! \return Unique pointer to 1D function -unique_ptr read_function(hid_t group, const char* name); - -} // namespace openmc - -#endif // OPENMC_ENDF_H +//! \file endf.h +//! Classes and functions related to the ENDF-6 format + +#ifndef OPENMC_ENDF_H +#define OPENMC_ENDF_H + +#include "hdf5.h" + +#include "openmc/constants.h" +#include "openmc/memory.h" +#include "openmc/vector.h" + +namespace openmc { + +//! Convert integer representing interpolation law to enum +//! \param[in] i Intereger (e.g. 1=histogram, 2=lin-lin) +//! \return Corresponding enum value +Interpolation int2interp(int i); + +//! Determine whether MT number corresponds to a fission reaction +//! \param[in] MT ENDF MT value +//! \return Whether corresponding reaction is a fission reaction +bool is_fission(int MT); + +//! Determine if a given MT number is that of a disappearance reaction, i.e., a +//! reaction with no neutron in the exit channel +//! \param[in] MT ENDF MT value +//! \return Whether corresponding reaction is a disappearance reaction +bool is_disappearance(int MT); + +//! Determine if a given MT number is that of an inelastic scattering reaction +//! \param[in] MT ENDF MT value +//! \return Whether corresponding reaction is an inelastic scattering reaction +bool is_inelastic_scatter(int MT); + +//============================================================================== +//! Abstract one-dimensional function +//============================================================================== + +class Function1D { +public: + virtual double operator()(double x) const = 0; + virtual ~Function1D() = default; +}; + +//============================================================================== +//! One-dimensional function expressed as a polynomial +//============================================================================== + +class Polynomial : public Function1D { +public: + //! Construct polynomial from HDF5 data + //! \param[in] dset Dataset containing coefficients + explicit Polynomial(hid_t dset); + + //! Construct polynomial from coefficients + //! \param[in] coef Polynomial coefficients + explicit Polynomial(vector coef) : coef_(coef) {} + + //! Evaluate the polynomials + //! \param[in] x independent variable + //! \return Polynomial evaluated at x + double operator()(double x) const override; + +private: + vector coef_; //!< Polynomial coefficients +}; + +//============================================================================== +//! One-dimensional interpolable function +//============================================================================== + +class Tabulated1D : public Function1D { +public: + Tabulated1D() = default; + + //! Construct function from HDF5 data + //! \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 + double operator()(double x) const override; + + // Accessors + const vector& x() const { return x_; } + const vector& y() const { return y_; } + +private: + std::size_t n_regions_ {0}; //!< number of interpolation regions + vector nbt_; //!< values separating interpolation regions + vector int_; //!< interpolation schemes + std::size_t n_pairs_; //!< number of (x,y) pairs + vector x_; //!< values of abscissa + vector y_; //!< values of ordinate +}; + +//============================================================================== +//! Coherent elastic scattering data from a crystalline material +//============================================================================== + +class CoherentElasticXS : public Function1D { +public: + explicit CoherentElasticXS(hid_t dset); + + double operator()(double E) const override; + + const vector& bragg_edges() const { return bragg_edges_; } + const vector& factors() const { return factors_; } + +private: + vector bragg_edges_; //!< Bragg edges in [eV] + vector factors_; //!< Partial sums of structure factors [eV-b] +}; + +//============================================================================== +//! Incoherent elastic scattering cross section +//============================================================================== + +class IncoherentElasticXS : public Function1D { +public: + explicit IncoherentElasticXS(hid_t dset); + + double operator()(double E) const override; + +private: + double bound_xs_; //!< Characteristic bound xs in [b] + double + debye_waller_; //!< Debye-Waller integral divided by atomic mass in [eV^-1] +}; + +//============================================================================== +//! Sum of multiple 1D functions +//============================================================================== + +class Sum1D : public Function1D { +public: + // Constructors + explicit Sum1D(hid_t group); + + //! Evaluate each function and sum results + //! \param[in] x independent variable + //! \return Function evaluated at x + double operator()(double E) const override; + + const unique_ptr& functions(int i) const { return functions_[i]; } + +private: + vector> functions_; //!< individual functions +}; + +//! Read 1D function from HDF5 dataset +//! \param[in] group HDF5 group containing dataset +//! \param[in] name Name of dataset +//! \return Unique pointer to 1D function +unique_ptr read_function(hid_t group, const char* name); + +} // namespace openmc + +#endif // OPENMC_ENDF_H diff --git a/src/chain.cpp b/src/chain.cpp index e892c4593ca..b1c3c24656b 100644 --- a/src/chain.cpp +++ b/src/chain.cpp @@ -1,170 +1,170 @@ -//! \file chain.cpp -//! \brief Depletion chain and associated information - -#include "openmc/chain.h" - -#include // for getenv -#include // for make_unique -#include // for stod -#include -#include - -#include -#include - -#include "openmc/distribution.h" // for distribution_from_xml -#include "openmc/error.h" -#include "openmc/reaction.h" -#include "openmc/xml_interface.h" // for get_node_value - -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 -//============================================================================== - -ChainNuclide::ChainNuclide(pugi::xml_node node) -{ - name_ = get_node_value(node, "name"); - if (check_for_node(node, "half_life")) { - half_life_ = std::stod(get_node_value(node, "half_life")); - } - if (check_for_node(node, "decay_energy")) { - decay_energy_ = std::stod(get_node_value(node, "decay_energy")); - } - - // Read reactions to store MT -> product map - for (pugi::xml_node reaction_node : node.children("reaction")) { - std::string rx_name = get_node_value(reaction_node, "type"); - if (rx_name == "fission") { - fission_energy_ = std::stod(get_node_value(reaction_node, "Q")); - fissionable_ = true; - } - if (!reaction_node.attribute("target")) - continue; - std::string rx_target = get_node_value(reaction_node, "target"); - double branching_ratio = 1.0; - if (reaction_node.attribute("branching_ratio")) { - branching_ratio = - std::stod(get_node_value(reaction_node, "branching_ratio")); - } - int mt = reaction_type(rx_name); - reaction_products_[mt].push_back({rx_target, branching_ratio}); - } - - for (pugi::xml_node source_node : node.children("source")) { - auto particle = get_node_value(source_node, "particle"); - if (particle == "photon") { - photon_energy_ = distribution_from_xml(source_node); - break; - } - } - - 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(); -} - -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 -//============================================================================== - -void DecayPhotonAngleEnergy::sample( - double E_in, double& E_out, double& mu, uint64_t* seed) const -{ - E_out = photon_energy_->sample(seed); - mu = Uniform(-1., 1.).sample(seed); -} - -//============================================================================== -// Global variables -//============================================================================== - -namespace data { - -std::unordered_map chain_nuclide_map; -vector> chain_nuclides; -vector> fission_yields; - -} // namespace data - -//============================================================================== -// Non-member functions -//============================================================================== - -void read_chain_file_xml() -{ - char* chain_file_path = std::getenv("OPENMC_CHAIN_FILE"); - if (!chain_file_path) { - return; - } - - write_message(5, "Reading chain file: {}...", chain_file_path); - - pugi::xml_document doc; - auto result = doc.load_file(chain_file_path); - if (!result) { - fatal_error( - fmt::format("Error processing chain file: {}", chain_file_path)); - } - - // Get root element - pugi::xml_node root = doc.document_element(); - - for (auto node : root.children("nuclide")) { - data::chain_nuclides.push_back(std::make_unique(node)); - } -} - -} // namespace openmc +//! \file chain.cpp +//! \brief Depletion chain and associated information + +#include "openmc/chain.h" + +#include // for getenv +#include // for make_unique +#include // for stod +#include +#include + +#include +#include + +#include "openmc/distribution.h" // for distribution_from_xml +#include "openmc/error.h" +#include "openmc/reaction.h" +#include "openmc/xml_interface.h" // for get_node_value + +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 +//============================================================================== + +ChainNuclide::ChainNuclide(pugi::xml_node node) +{ + name_ = get_node_value(node, "name"); + if (check_for_node(node, "half_life")) { + half_life_ = std::stod(get_node_value(node, "half_life")); + } + if (check_for_node(node, "decay_energy")) { + decay_energy_ = std::stod(get_node_value(node, "decay_energy")); + } + + // Read reactions to store MT -> product map + for (pugi::xml_node reaction_node : node.children("reaction")) { + std::string rx_name = get_node_value(reaction_node, "type"); + if (rx_name == "fission") { + fission_energy_ = std::stod(get_node_value(reaction_node, "Q")); + fissionable_ = true; + } + if (!reaction_node.attribute("target")) + continue; + std::string rx_target = get_node_value(reaction_node, "target"); + double branching_ratio = 1.0; + if (reaction_node.attribute("branching_ratio")) { + branching_ratio = + std::stod(get_node_value(reaction_node, "branching_ratio")); + } + int mt = reaction_type(rx_name); + reaction_products_[mt].push_back({rx_target, branching_ratio}); + } + + for (pugi::xml_node source_node : node.children("source")) { + auto particle = get_node_value(source_node, "particle"); + if (particle == "photon") { + photon_energy_ = distribution_from_xml(source_node); + break; + } + } + + 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(); +} + +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 +//============================================================================== + +void DecayPhotonAngleEnergy::sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const +{ + E_out = photon_energy_->sample(seed); + mu = Uniform(-1., 1.).sample(seed); +} + +//============================================================================== +// Global variables +//============================================================================== + +namespace data { + +std::unordered_map chain_nuclide_map; +vector> chain_nuclides; +vector> fission_yields; + +} // namespace data + +//============================================================================== +// Non-member functions +//============================================================================== + +void read_chain_file_xml() +{ + char* chain_file_path = std::getenv("OPENMC_CHAIN_FILE"); + if (!chain_file_path) { + return; + } + + write_message(5, "Reading chain file: {}...", chain_file_path); + + pugi::xml_document doc; + auto result = doc.load_file(chain_file_path); + if (!result) { + fatal_error( + fmt::format("Error processing chain file: {}", chain_file_path)); + } + + // Get root element + pugi::xml_node root = doc.document_element(); + + for (auto node : root.children("nuclide")) { + data::chain_nuclides.push_back(std::make_unique(node)); + } +} + +} // namespace openmc diff --git a/src/endf.cpp b/src/endf.cpp index 5608e195cbf..b01d160640e 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -1,310 +1,310 @@ -#include "openmc/endf.h" - -#include // for copy -#include // for log, exp -#include // for back_inserter -#include // for runtime_error - -#include "xtensor/xarray.hpp" -#include "xtensor/xview.hpp" - -#include "openmc/array.h" -#include "openmc/constants.h" -#include "openmc/hdf5_interface.h" -#include "openmc/search.h" - -namespace openmc { - -//============================================================================== -// Functions -//============================================================================== - -Interpolation int2interp(int i) -{ - // TODO: We are ignoring specification of two-dimensional interpolation - // schemes (method of corresponding points and unit base interpolation). Those - // should be accounted for in the distribution classes somehow. - - switch (i) { - case 1: - case 11: - case 21: - return Interpolation::histogram; - case 2: - case 12: - case 22: - return Interpolation::lin_lin; - case 3: - case 13: - case 23: - return Interpolation::lin_log; - case 4: - case 14: - case 24: - return Interpolation::log_lin; - case 5: - case 15: - case 25: - return Interpolation::log_log; - default: - throw std::runtime_error {"Invalid interpolation code."}; - } -} - -bool is_fission(int mt) -{ - return mt == N_FISSION || mt == N_F || mt == N_NF || mt == N_2NF || - mt == N_3NF; -} - -bool is_disappearance(int mt) -{ - if (mt >= N_DISAPPEAR && mt <= N_DA) { - return true; - } else if (mt >= N_P0 && mt <= N_AC) { - return true; - } else if (mt == N_TA || mt == N_DT || mt == N_P3HE || mt == N_D3HE || - mt == N_3HEA || mt == N_3P) { - return true; - } else { - return false; - } -} - -bool is_inelastic_scatter(int mt) -{ - if (mt < 100) { - if (is_fission(mt)) { - return false; - } else { - return mt >= MISC && mt != 27; - } - } else if (mt <= 200) { - return !is_disappearance(mt); - } else if (mt >= N_2N0 && mt <= N_2NC) { - return true; - } else { - return false; - } -} - -unique_ptr read_function(hid_t group, const char* name) -{ - hid_t obj_id = open_object(group, name); - std::string func_type; - read_attribute(obj_id, "type", func_type); - unique_ptr func; - if (func_type == "Tabulated1D") { - func = make_unique(obj_id); - } else if (func_type == "Polynomial") { - func = make_unique(obj_id); - } else if (func_type == "CoherentElastic") { - func = make_unique(obj_id); - } else if (func_type == "IncoherentElastic") { - func = make_unique(obj_id); - } else if (func_type == "Sum") { - func = make_unique(obj_id); - } else { - throw std::runtime_error {"Unknown function type " + func_type + - " for dataset " + object_name(obj_id)}; - } - close_object(obj_id); - return func; -} - -//============================================================================== -// Polynomial implementation -//============================================================================== - -Polynomial::Polynomial(hid_t dset) -{ - // Read coefficients into a vector - read_dataset(dset, coef_); -} - -double Polynomial::operator()(double x) const -{ - // Use Horner's rule to evaluate polynomial. Note that coefficients are - // ordered in increasing powers of x. - double y = 0.0; - for (auto c = coef_.crbegin(); c != coef_.crend(); ++c) { - y = y * x + *c; - } - return y; -} - -//============================================================================== -// Tabulated1D implementation -//============================================================================== - -Tabulated1D::Tabulated1D(hid_t dset) -{ - read_attribute(dset, "breakpoints", nbt_); - n_regions_ = nbt_.size(); - - // Change 1-indexing to 0-indexing - for (auto& b : nbt_) - --b; - - vector int_temp; - read_attribute(dset, "interpolation", int_temp); - - // Convert vector of ints into Interpolation - for (const auto i : int_temp) - int_.push_back(int2interp(i)); - - xt::xarray arr; - read_dataset(dset, arr); - - auto xs = xt::view(arr, 0); - auto ys = xt::view(arr, 1); - - std::copy(xs.begin(), xs.end(), std::back_inserter(x_)); - std::copy(ys.begin(), ys.end(), std::back_inserter(y_)); - 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 - // tabulated range, the first or last point is chosen, i.e. no interpolation - // is done outside the energy range - int i; - if (x < x_[0]) { - return y_[0]; - } else if (x > x_[n_pairs_ - 1]) { - return y_[n_pairs_ - 1]; - } else { - i = lower_bound_index(x_.begin(), x_.end(), x); - } - - // determine interpolation scheme - Interpolation interp; - if (n_regions_ == 0) { - interp = Interpolation::lin_lin; - } else { - interp = int_[0]; - for (int j = 0; j < n_regions_; ++j) { - if (i < nbt_[j]) { - interp = int_[j]; - break; - } - } - } - - // handle special case of histogram interpolation - if (interp == Interpolation::histogram) - return y_[i]; - - // determine bounding values - double x0 = x_[i]; - double x1 = x_[i + 1]; - double y0 = y_[i]; - double y1 = y_[i + 1]; - - // determine interpolation factor and interpolated value - double r; - switch (interp) { - case Interpolation::lin_lin: - r = (x - x0) / (x1 - x0); - return y0 + r * (y1 - y0); - case Interpolation::lin_log: - r = log(x / x0) / log(x1 / x0); - return y0 + r * (y1 - y0); - case Interpolation::log_lin: - r = (x - x0) / (x1 - x0); - return y0 * exp(r * log(y1 / y0)); - case Interpolation::log_log: - r = log(x / x0) / log(x1 / x0); - return y0 * exp(r * log(y1 / y0)); - default: - throw std::runtime_error {"Invalid interpolation scheme."}; - } -} - -//============================================================================== -// CoherentElasticXS implementation -//============================================================================== - -CoherentElasticXS::CoherentElasticXS(hid_t dset) -{ - // Read 2D array from dataset - xt::xarray arr; - read_dataset(dset, arr); - - // Get views for Bragg edges and structure factors - auto E = xt::view(arr, 0); - auto s = xt::view(arr, 1); - - // Copy Bragg edges and partial sums of structure factors - std::copy(E.begin(), E.end(), std::back_inserter(bragg_edges_)); - std::copy(s.begin(), s.end(), std::back_inserter(factors_)); -} - -double CoherentElasticXS::operator()(double E) const -{ - if (E < bragg_edges_[0]) { - // If energy is below that of the lowest Bragg peak, the elastic cross - // section will be zero - return 0.0; - } else { - auto i_grid = - lower_bound_index(bragg_edges_.begin(), bragg_edges_.end(), E); - return factors_[i_grid] / E; - } -} - -//============================================================================== -// IncoherentElasticXS implementation -//============================================================================== - -IncoherentElasticXS::IncoherentElasticXS(hid_t dset) -{ - array tmp; - read_dataset(dset, nullptr, tmp); - bound_xs_ = tmp[0]; - debye_waller_ = tmp[1]; -} - -double IncoherentElasticXS::operator()(double E) const -{ - // Determine cross section using ENDF-102, Eq. (7.5) - double W = debye_waller_; - return bound_xs_ / 2.0 * ((1 - std::exp(-4.0 * E * W)) / (2.0 * E * W)); -} - -//============================================================================== -// Sum1D implementation -//============================================================================== - -Sum1D::Sum1D(hid_t group) -{ - // Get number of functions - int n; - read_attribute(group, "n", n); - - // Get each function - for (int i = 0; i < n; ++i) { - auto dset_name = fmt::format("func_{}", i + 1); - functions_.push_back(read_function(group, dset_name.c_str())); - } -} - -double Sum1D::operator()(double x) const -{ - double result = 0.0; - for (auto& func : functions_) { - result += (*func)(x); - } - return result; -} - -} // namespace openmc +#include "openmc/endf.h" + +#include // for copy +#include // for log, exp +#include // for back_inserter +#include // for runtime_error + +#include "xtensor/xarray.hpp" +#include "xtensor/xview.hpp" + +#include "openmc/array.h" +#include "openmc/constants.h" +#include "openmc/hdf5_interface.h" +#include "openmc/search.h" + +namespace openmc { + +//============================================================================== +// Functions +//============================================================================== + +Interpolation int2interp(int i) +{ + // TODO: We are ignoring specification of two-dimensional interpolation + // schemes (method of corresponding points and unit base interpolation). Those + // should be accounted for in the distribution classes somehow. + + switch (i) { + case 1: + case 11: + case 21: + return Interpolation::histogram; + case 2: + case 12: + case 22: + return Interpolation::lin_lin; + case 3: + case 13: + case 23: + return Interpolation::lin_log; + case 4: + case 14: + case 24: + return Interpolation::log_lin; + case 5: + case 15: + case 25: + return Interpolation::log_log; + default: + throw std::runtime_error {"Invalid interpolation code."}; + } +} + +bool is_fission(int mt) +{ + return mt == N_FISSION || mt == N_F || mt == N_NF || mt == N_2NF || + mt == N_3NF; +} + +bool is_disappearance(int mt) +{ + if (mt >= N_DISAPPEAR && mt <= N_DA) { + return true; + } else if (mt >= N_P0 && mt <= N_AC) { + return true; + } else if (mt == N_TA || mt == N_DT || mt == N_P3HE || mt == N_D3HE || + mt == N_3HEA || mt == N_3P) { + return true; + } else { + return false; + } +} + +bool is_inelastic_scatter(int mt) +{ + if (mt < 100) { + if (is_fission(mt)) { + return false; + } else { + return mt >= MISC && mt != 27; + } + } else if (mt <= 200) { + return !is_disappearance(mt); + } else if (mt >= N_2N0 && mt <= N_2NC) { + return true; + } else { + return false; + } +} + +unique_ptr read_function(hid_t group, const char* name) +{ + hid_t obj_id = open_object(group, name); + std::string func_type; + read_attribute(obj_id, "type", func_type); + unique_ptr func; + if (func_type == "Tabulated1D") { + func = make_unique(obj_id); + } else if (func_type == "Polynomial") { + func = make_unique(obj_id); + } else if (func_type == "CoherentElastic") { + func = make_unique(obj_id); + } else if (func_type == "IncoherentElastic") { + func = make_unique(obj_id); + } else if (func_type == "Sum") { + func = make_unique(obj_id); + } else { + throw std::runtime_error {"Unknown function type " + func_type + + " for dataset " + object_name(obj_id)}; + } + close_object(obj_id); + return func; +} + +//============================================================================== +// Polynomial implementation +//============================================================================== + +Polynomial::Polynomial(hid_t dset) +{ + // Read coefficients into a vector + read_dataset(dset, coef_); +} + +double Polynomial::operator()(double x) const +{ + // Use Horner's rule to evaluate polynomial. Note that coefficients are + // ordered in increasing powers of x. + double y = 0.0; + for (auto c = coef_.crbegin(); c != coef_.crend(); ++c) { + y = y * x + *c; + } + return y; +} + +//============================================================================== +// Tabulated1D implementation +//============================================================================== + +Tabulated1D::Tabulated1D(hid_t dset) +{ + read_attribute(dset, "breakpoints", nbt_); + n_regions_ = nbt_.size(); + + // Change 1-indexing to 0-indexing + for (auto& b : nbt_) + --b; + + vector int_temp; + read_attribute(dset, "interpolation", int_temp); + + // Convert vector of ints into Interpolation + for (const auto i : int_temp) + int_.push_back(int2interp(i)); + + xt::xarray arr; + read_dataset(dset, arr); + + auto xs = xt::view(arr, 0); + auto ys = xt::view(arr, 1); + + std::copy(xs.begin(), xs.end(), std::back_inserter(x_)); + std::copy(ys.begin(), ys.end(), std::back_inserter(y_)); + 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 + // tabulated range, the first or last point is chosen, i.e. no interpolation + // is done outside the energy range + int i; + if (x < x_[0]) { + return y_[0]; + } else if (x > x_[n_pairs_ - 1]) { + return y_[n_pairs_ - 1]; + } else { + i = lower_bound_index(x_.begin(), x_.end(), x); + } + + // determine interpolation scheme + Interpolation interp; + if (n_regions_ == 0) { + interp = Interpolation::lin_lin; + } else { + interp = int_[0]; + for (int j = 0; j < n_regions_; ++j) { + if (i < nbt_[j]) { + interp = int_[j]; + break; + } + } + } + + // handle special case of histogram interpolation + if (interp == Interpolation::histogram) + return y_[i]; + + // determine bounding values + double x0 = x_[i]; + double x1 = x_[i + 1]; + double y0 = y_[i]; + double y1 = y_[i + 1]; + + // determine interpolation factor and interpolated value + double r; + switch (interp) { + case Interpolation::lin_lin: + r = (x - x0) / (x1 - x0); + return y0 + r * (y1 - y0); + case Interpolation::lin_log: + r = log(x / x0) / log(x1 / x0); + return y0 + r * (y1 - y0); + case Interpolation::log_lin: + r = (x - x0) / (x1 - x0); + return y0 * exp(r * log(y1 / y0)); + case Interpolation::log_log: + r = log(x / x0) / log(x1 / x0); + return y0 * exp(r * log(y1 / y0)); + default: + throw std::runtime_error {"Invalid interpolation scheme."}; + } +} + +//============================================================================== +// CoherentElasticXS implementation +//============================================================================== + +CoherentElasticXS::CoherentElasticXS(hid_t dset) +{ + // Read 2D array from dataset + xt::xarray arr; + read_dataset(dset, arr); + + // Get views for Bragg edges and structure factors + auto E = xt::view(arr, 0); + auto s = xt::view(arr, 1); + + // Copy Bragg edges and partial sums of structure factors + std::copy(E.begin(), E.end(), std::back_inserter(bragg_edges_)); + std::copy(s.begin(), s.end(), std::back_inserter(factors_)); +} + +double CoherentElasticXS::operator()(double E) const +{ + if (E < bragg_edges_[0]) { + // If energy is below that of the lowest Bragg peak, the elastic cross + // section will be zero + return 0.0; + } else { + auto i_grid = + lower_bound_index(bragg_edges_.begin(), bragg_edges_.end(), E); + return factors_[i_grid] / E; + } +} + +//============================================================================== +// IncoherentElasticXS implementation +//============================================================================== + +IncoherentElasticXS::IncoherentElasticXS(hid_t dset) +{ + array tmp; + read_dataset(dset, nullptr, tmp); + bound_xs_ = tmp[0]; + debye_waller_ = tmp[1]; +} + +double IncoherentElasticXS::operator()(double E) const +{ + // Determine cross section using ENDF-102, Eq. (7.5) + double W = debye_waller_; + return bound_xs_ / 2.0 * ((1 - std::exp(-4.0 * E * W)) / (2.0 * E * W)); +} + +//============================================================================== +// Sum1D implementation +//============================================================================== + +Sum1D::Sum1D(hid_t group) +{ + // Get number of functions + int n; + read_attribute(group, "n", n); + + // Get each function + for (int i = 0; i < n; ++i) { + auto dset_name = fmt::format("func_{}", i + 1); + functions_.push_back(read_function(group, dset_name.c_str())); + } +} + +double Sum1D::operator()(double x) const +{ + double result = 0.0; + for (auto& func : functions_) { + result += (*func)(x); + } + return result; +} + +} // namespace openmc From 8db2f029584526f1e67b9471476e5545e182769a Mon Sep 17 00:00:00 2001 From: GuySten Date: Mon, 14 Jul 2025 22:38:12 +0300 Subject: [PATCH 05/11] fixed estimator --- src/tallies/tally.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 94a1ebfa924..3022bcdd556 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -588,8 +588,8 @@ void Tally::set_scores(const vector& scores) fatal_error("Cannot tally " + score_str + " reaction rate with an " "outgoing energy filter"); - if (fission_yields_present) - estimator_ = TallyEstimator::ANALOG; + if (fission_yields_present && estimator_ == TallyEstimator::TRACKLENGTH) + estimator_ = TallyEstimator::COLLISION; break; case SCORE_SCATTER: From ac5553fd6d0ffb2b577d72d8966b1645437c9e95 Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 15 Jul 2025 17:23:20 +0300 Subject: [PATCH 06/11] fixed test --- tests/unit_tests/test_filter_fission_yields.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_filter_fission_yields.py b/tests/unit_tests/test_filter_fission_yields.py index 2bc4b41f534..a8e5760227f 100644 --- a/tests/unit_tests/test_filter_fission_yields.py +++ b/tests/unit_tests/test_filter_fission_yields.py @@ -19,7 +19,7 @@ def test_fissionyieldsfilter(run_in_tmpdir): ) model.settings.run_mode = "eigenvalue" - fissionyields_filter = openmc.FissionYieldsFilter(['Xe135','Te135','I135']) + fissionyields_filter = openmc.FissionYieldsFilter(['Xe135','I135']) tally = openmc.Tally() tally.filters = [fissionyields_filter] @@ -32,7 +32,7 @@ def test_fissionyieldsfilter(run_in_tmpdir): # Run OpenMC model.run(apply_tally_results=True) - assert np.all(tally.mean.reshape(3)>0) + assert np.all(tally.mean.reshape(2)>0) From ba0538b6a33a176c5ee4562ceba1504251f51248 Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 15 Jul 2025 23:15:10 +0300 Subject: [PATCH 07/11] added docs --- docs/source/pythonapi/base.rst | 1 + docs/source/pythonapi/capi.rst | 1 + openmc/filter.py | 2 +- openmc/lib/filter.py | 3 +++ 4 files changed, 6 insertions(+), 1 deletion(-) 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/openmc/filter.py b/openmc/filter.py index 11aeed9be1a..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', 'fission_yields', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre', + 'energyfunction', 'fissionyields', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre', 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance', 'collision', 'time', 'parentnuclide', 'weight', 'meshborn', 'meshsurface', 'meshmaterial', 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, From ee7af4e744af543decf0cc8755e129ad5d39605f Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 15 Jul 2025 23:34:17 +0300 Subject: [PATCH 08/11] nuclides->bins --- include/openmc/tallies/filter_fission_yields.h | 4 ++-- src/tallies/filter_fission_yields.cpp | 16 ++++++++-------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/include/openmc/tallies/filter_fission_yields.h b/include/openmc/tallies/filter_fission_yields.h index cddef070e38..b08ca699dde 100644 --- a/include/openmc/tallies/filter_fission_yields.h +++ b/include/openmc/tallies/filter_fission_yields.h @@ -35,13 +35,13 @@ class FissionYieldsFilter : public Filter { //---------------------------------------------------------------------------- // Accessors - const vector& nuclides() const { return nuclides_; } + const vector& bins() const { return bins_; } private: //---------------------------------------------------------------------------- // Data members - vector nuclides_; + vector bins_; }; } // namespace openmc diff --git a/src/tallies/filter_fission_yields.cpp b/src/tallies/filter_fission_yields.cpp index cbbff09af91..2c7b03bd5ae 100644 --- a/src/tallies/filter_fission_yields.cpp +++ b/src/tallies/filter_fission_yields.cpp @@ -17,10 +17,10 @@ void FissionYieldsFilter::from_xml(pugi::xml_node node) "continuous-energy transport calculations"); if (!check_for_node(node, "bins")) - fatal_error("Nuclides not specified for FissionYieldsFilter."); + fatal_error("Bins not specified for FissionYieldsFilter."); - nuclides_ = get_node_array(node, "bins"); - n_bins_ = nuclides_.size(); + bins_ = get_node_array(node, "bins"); + n_bins_ = bins_.size(); } void FissionYieldsFilter::get_all_bins( @@ -32,10 +32,10 @@ void FissionYieldsFilter::get_all_bins( auto fy = data::chain_nuclides[data::chain_nuclide_map[nuc]] ->fission_yields() ->yields_; - for (int i = 0; i < nuclides_.size(); ++i) { - if (fy.find(nuclides_[i]) != fy.end()) { + 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[nuclides_[i]](p.E_last())); + match.weights_.push_back(fy[bins_[i]](p.E_last())); } } } @@ -45,12 +45,12 @@ void FissionYieldsFilter::get_all_bins( void FissionYieldsFilter::to_statepoint(hid_t filter_group) const { Filter::to_statepoint(filter_group); - write_dataset(filter_group, "nuclides", nuclides_); + write_dataset(filter_group, "bins", bins_); } std::string FissionYieldsFilter::text_label(int bin) const { - return fmt::format("Fission Yield [{}]", nuclides_[bin]); + return fmt::format("Fission Yield [{}]", bins_[bin]); } } // namespace openmc From 62192e5fba67772dd19640724cb5fc75c9d47768 Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Tue, 15 Jul 2025 23:35:49 +0300 Subject: [PATCH 09/11] Update filter_fission_yields.h --- include/openmc/tallies/filter_fission_yields.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/include/openmc/tallies/filter_fission_yields.h b/include/openmc/tallies/filter_fission_yields.h index b08ca699dde..ac58f6650fe 100644 --- a/include/openmc/tallies/filter_fission_yields.h +++ b/include/openmc/tallies/filter_fission_yields.h @@ -6,8 +6,7 @@ namespace openmc { //============================================================================== -//! Multiplies tally scores by an arbitrary function of incident energy -//! described by a piecewise linear-linear interpolation. +//! Multiplies tally scores by fission yields. //============================================================================== class FissionYieldsFilter : public Filter { From d4aa3de61409546e7eff88a46e1e2ff03ef0d986 Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 15 Jul 2025 23:38:43 +0300 Subject: [PATCH 10/11] simplified hain --- include/openmc/chain.h | 2 -- src/chain.cpp | 4 ---- 2 files changed, 6 deletions(-) diff --git a/include/openmc/chain.h b/include/openmc/chain.h index 12353fcb6d0..5b6a7f61ccb 100644 --- a/include/openmc/chain.h +++ b/include/openmc/chain.h @@ -63,8 +63,6 @@ class ChainNuclide { std::string name_; //!< Name of nuclide double half_life_ {0.0}; //!< Half-life in [s] double decay_energy_ {0.0}; //!< Decay energy in [eV] - bool fissionable_ {false}; //!< Can do fission - double fission_energy_ {0.0}; //!< Fission energy in [eV] FissionYields* fission_yields_ = nullptr; //!< Fission yields std::string fission_yields_parent_ {""}; //!< Fission yields parent std::unordered_map> diff --git a/src/chain.cpp b/src/chain.cpp index b1c3c24656b..418a9fdf931 100644 --- a/src/chain.cpp +++ b/src/chain.cpp @@ -64,10 +64,6 @@ ChainNuclide::ChainNuclide(pugi::xml_node node) // Read reactions to store MT -> product map for (pugi::xml_node reaction_node : node.children("reaction")) { std::string rx_name = get_node_value(reaction_node, "type"); - if (rx_name == "fission") { - fission_energy_ = std::stod(get_node_value(reaction_node, "Q")); - fissionable_ = true; - } if (!reaction_node.attribute("target")) continue; std::string rx_target = get_node_value(reaction_node, "target"); From d100d055f4d858f498e7956c820e07e7170bea90 Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Tue, 15 Jul 2025 23:42:43 +0300 Subject: [PATCH 11/11] Update test_filter_fission_yields.py --- tests/unit_tests/test_filter_fission_yields.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/unit_tests/test_filter_fission_yields.py b/tests/unit_tests/test_filter_fission_yields.py index a8e5760227f..3e3dd177715 100644 --- a/tests/unit_tests/test_filter_fission_yields.py +++ b/tests/unit_tests/test_filter_fission_yields.py @@ -23,7 +23,6 @@ def test_fissionyieldsfilter(run_in_tmpdir): tally = openmc.Tally() tally.filters = [fissionyields_filter] - tally.estimator = 'analog' tally.scores = ['fission'] model.tallies = openmc.Tallies([tally])