From 1500f863b6be1d2a2f2c696e3fac74fb34d04744 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 14 Nov 2025 10:45:45 -0500 Subject: [PATCH] fix(Reaction): resolved bug which prevented counting number of reactants and products properly This bug was introduced when switching the internal tracking of products and reactants from a vector to a set. Initially we had done this to improve lookup performance. However, due to the uniquness of items in a set this broke counting methods. We have resolved this so that all old code will still work while maintaing the efficency gains by using an unordered map which stored species and counts. Further we have added countReactantOccurences and countProductOffurences methods to make this process more explicit and avalible to callers. --- src/include/gridfire/engine/types/jacobian.h | 0 src/include/gridfire/reaction/reaction.h | 32 +++++++++++++---- src/include/gridfire/reaction/weak/weak.h | 5 +++ src/lib/engine/types/jacobian.cpp | 0 src/lib/reaction/reaction.cpp | 38 ++++++++++---------- src/lib/reaction/weak/weak.cpp | 20 +++++++++++ 6 files changed, 71 insertions(+), 24 deletions(-) create mode 100644 src/include/gridfire/engine/types/jacobian.h create mode 100644 src/lib/engine/types/jacobian.cpp diff --git a/src/include/gridfire/engine/types/jacobian.h b/src/include/gridfire/engine/types/jacobian.h new file mode 100644 index 00000000..e69de29b diff --git a/src/include/gridfire/reaction/reaction.h b/src/include/gridfire/reaction/reaction.h index 4b054333..997b8396 100644 --- a/src/include/gridfire/reaction/reaction.h +++ b/src/include/gridfire/reaction/reaction.h @@ -217,6 +217,10 @@ namespace gridfire::reaction { */ [[nodiscard]] virtual std::unordered_set product_species() const = 0; + [[nodiscard]] virtual size_t countReactantOccurrences(const fourdst::atomic::Species& species) const = 0; + + [[nodiscard]] virtual size_t countProductOccurrences(const fourdst::atomic::Species& species) const = 0; + /** * @brief Number of unique species involved in the reaction. * @return Count of distinct species across reactants and products. @@ -525,8 +529,10 @@ namespace gridfire::reaction { if (!m_reactantsVec) { m_reactantsVec.emplace(std::vector()); m_reactantsVec->reserve(m_reactants.size()); - for (const auto& reactant : m_reactants) { - m_reactantsVec->push_back(reactant); + for (const auto& [reactant, count] : m_reactants) { + for (size_t i = 0; i < count; ++i) { + m_reactantsVec->push_back(reactant); + } } } return m_reactantsVec.value(); @@ -540,8 +546,10 @@ namespace gridfire::reaction { if (!m_productsVec) { m_productsVec.emplace(std::vector()); m_productsVec->reserve(m_products.size()); - for (const auto& product : m_products) { - m_productsVec->push_back(product); + for (const auto& [product, count] : m_products) { + for (size_t i = 0; i < count; ++i) { + m_productsVec->push_back(product); + } } } return m_productsVec.value(); @@ -587,14 +595,24 @@ namespace gridfire::reaction { return os << "(ReaclibReaction:" << r.m_id << ")"; } + size_t countReactantOccurrences(const fourdst::atomic::Species& species) const override { + if (!m_reactants.contains(species)) return 0; + return m_reactants.at(species); + } + + size_t countProductOccurrences(const fourdst::atomic::Species& species) const override { + if (!m_products.contains(species)) return 0; + return m_products.at(species); + } + protected: quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); std::string m_id; ///< Unique identifier for the reaction (e.g., "h1+h1=>h2+e+nu"). std::string m_peName; ///< Name of the reaction in (projectile, ejectile) notation (e.g. "p(p,g)d"). int m_chapter; ///< Chapter number from the REACLIB database, defining the reaction structure. double m_qValue = 0.0; ///< Q-value of the reaction in MeV. - std::set m_reactants; ///< Reactants of the reaction. - std::set m_products; ///< Products of the reaction. + std::unordered_map m_reactants; ///< Reactants of the reaction. + std::unordered_map m_products; ///< Products of the reaction. mutable std::optional> m_reactantsVec; mutable std::optional> m_productsVec; @@ -836,6 +854,8 @@ namespace gridfire::reaction { */ [[nodiscard]] size_t size() const { return m_reactions.size(); } + [[nodiscard]] bool empty() const {return m_reactions.empty(); } + /** * @brief Removes all reactions from the set. */ diff --git a/src/include/gridfire/reaction/weak/weak.h b/src/include/gridfire/reaction/weak/weak.h index 00659fe1..0428a1e9 100644 --- a/src/include/gridfire/reaction/weak/weak.h +++ b/src/include/gridfire/reaction/weak/weak.h @@ -400,6 +400,11 @@ namespace gridfire::rates::weak { */ double get_log_neutrino_loss_from_payload(const WeakRatePayload& payload) const; + public: + [[nodiscard]] size_t countReactantOccurrences(const fourdst::atomic::Species &species) const override; + + [[nodiscard]] size_t countProductOccurrences(const fourdst::atomic::Species &species) const override; + private: /** * @brief CppAD atomic that wraps weak-rate interpolation for AD evaluation. diff --git a/src/lib/engine/types/jacobian.cpp b/src/lib/engine/types/jacobian.cpp new file mode 100644 index 00000000..e69de29b diff --git a/src/lib/reaction/reaction.cpp b/src/lib/reaction/reaction.cpp index 3a227243..6bdd8af7 100644 --- a/src/lib/reaction/reaction.cpp +++ b/src/lib/reaction/reaction.cpp @@ -40,11 +40,16 @@ namespace gridfire::reaction { m_peName(peName), m_chapter(chapter), m_qValue(qValue), - m_reactants(reactants.begin(), reactants.end()), - m_products(products.begin(), products.end()), m_sourceLabel(label), m_rateCoefficients(sets), - m_reverse(reverse) {} + m_reverse(reverse) { + for (const auto& reactant : reactants) { + m_reactants[reactant]++; + } + for (const auto& product : products) { + m_products[product]++; + } + } double ReaclibReaction::calculate_rate( const double T9, @@ -107,11 +112,11 @@ namespace gridfire::reaction { bool ReaclibReaction::contains_reactant(const Species& species) const { - return std::ranges::any_of(m_reactants, [&](const Species& r) { return r == species; }); + return m_reactants.contains(species); } bool ReaclibReaction::contains_product(const Species& species) const { - return std::ranges::any_of(m_products, [&](const Species& p) { return p == species; }); + return m_products.contains(species); } std::unordered_set ReaclibReaction::all_species() const { @@ -123,7 +128,7 @@ namespace gridfire::reaction { std::unordered_set ReaclibReaction::reactant_species() const { std::unordered_set reactantsSet; - for (const auto& reactant : m_reactants) { + for (const auto& reactant : m_reactants | std::views::keys) { reactantsSet.insert(reactant); } return reactantsSet; @@ -131,7 +136,7 @@ namespace gridfire::reaction { std::unordered_set ReaclibReaction::product_species() const { std::unordered_set productsSet; - for (const auto& product : m_products) { + for (const auto& product : m_products | std::views::keys) { productsSet.insert(product); } return productsSet; @@ -139,14 +144,14 @@ namespace gridfire::reaction { int ReaclibReaction::stoichiometry(const Species& species) const { int s = 0; - for (const auto& reactant : m_reactants) { + for (const auto& [reactant, count] : m_reactants) { if (reactant == species) { - s--; + s -= count; } } - for (const auto& product : m_products) { + for (const auto& [product, count] : m_products) { if (product == species) { - s++; + s += count; } } return s; @@ -158,11 +163,8 @@ namespace gridfire::reaction { std::unordered_map ReaclibReaction::stoichiometry() const { std::unordered_map stoichiometryMap; - for (const auto& reactant : m_reactants) { - stoichiometryMap[reactant]--; - } - for (const auto& product : m_products) { - stoichiometryMap[product]++; + for (const auto& sp : all_species()) { + stoichiometryMap[sp] = stoichiometry(sp); } return stoichiometryMap; } @@ -171,10 +173,10 @@ namespace gridfire::reaction { double reactantMass = 0.0; double productMass = 0.0; constexpr double AMU2MeV = 931.494893; // Conversion factor from atomic mass unit to MeV - for (const auto& reactant : m_reactants) { + for (const auto& reactant : m_reactants | std::views::keys) { reactantMass += reactant.mass(); } - for (const auto& product : m_products) { + for (const auto& product : m_products | std::views::keys) { productMass += product.mass(); } return (reactantMass - productMass) * AMU2MeV; diff --git a/src/lib/reaction/weak/weak.cpp b/src/lib/reaction/weak/weak.cpp index 93e75c08..34d57eae 100644 --- a/src/lib/reaction/weak/weak.cpp +++ b/src/lib/reaction/weak/weak.cpp @@ -490,6 +490,26 @@ namespace gridfire::rates::weak { return logNeutrinoLoss; } + size_t WeakReaction::countReactantOccurrences(const fourdst::atomic::Species &species) const { + size_t count = 0; + for (const auto& reactant : m_reactants) { + if (reactant == species) { + count++; + } + } + return count; + } + + size_t WeakReaction::countProductOccurrences(const fourdst::atomic::Species &species) const { + size_t count = 0; + for (const auto& product : m_products) { + if (product == species) { + count++; + } + } + return count; + } + // Note that the input vector tx is of size 2: [T9, log10(rho*Ye)] bool WeakReaction::AtomicWeakRate::forward ( const size_t p,