feat(neutrino): Started framework for neutrino loss

Neutrino loss is essential for neutrino cooling. Started adding
framework to track this. Reaclib reactions use a simple heuristic where
electron capture reactions loss 100% of their energy to neutrinos
whereas beta decay reactions loose 50% of their energy to neutrinos

src/include/gridfire/engine/engine_abstract.h

@@ -54,6 +54,9 @@ namespace gridfire::engine {
         std::map<fourdst::atomic::Species, T> dydt{}; ///< Derivatives of abundances (dY/dt for each species).
         T nuclearEnergyGenerationRate = T(0.0); ///< Specific energy generation rate (e.g., erg/g/s).
         std::map<fourdst::atomic::Species, std::unordered_map<std::string, T>> reactionContributions{};
+        T neutrinoEnergyLossRate = T(0.0); // (erg/g/s)
+        T totalNeutrinoFlux = T(0.0);      // (neutrinos/g/s)
+
 
         StepDerivatives() : dydt(), nuclearEnergyGenerationRate(T(0.0)) {}
     };

src/include/gridfire/engine/engine_graph.h

@@ -778,6 +778,7 @@ namespace gridfire::engine {
             const double Na = Constants::getInstance().get("N_a").value; ///< Avogadro's number.
             const double c = Constants::getInstance().get("c").value; ///< Speed of light in cm/s.
             const double kB = Constants::getInstance().get("kB").value; ///< Boltzmann constant in erg/K.
+            const double MeV_to_erg = Constants::getInstance().get("MeV_to_erg").value; ///< Conversion factor from MeV to erg.
         };
 
         enum class JacobianMatrixState {

src/include/gridfire/reaction/reaclib.h

@@ -16,4 +16,12 @@ namespace gridfire::reaclib {
      */
     const reaction::ReactionSet &get_all_reaclib_reactions();
 
+    // Simple heuristic to check if a reaclib reaction is a strong or weak reaction
+    /*  A weak reaction is defined here as one where:
+        - The number of reactants is equal to the number of products
+        - There is only one reactant and one product
+        - The mass number (A) of the reactant is equal to the mass number (A) of the product
+    */
+    bool reaction_is_weak(const reaction::Reaction& reaction);
+
 } // namespace gridfire::reaclib

src/include/gridfire/reaction/reaction.h

@@ -27,19 +27,26 @@ namespace gridfire::reaction {
     enum class ReactionType {
         WEAK,
         REACLIB,
+        REACLIB_WEAK,
         LOGICAL_REACLIB,
+        LOGICAL_REACLIB_WEAK,
     };
 
     static std::unordered_map<ReactionType, std::string> ReactionTypeNames = {
         {ReactionType::WEAK, "weak"},
         {ReactionType::REACLIB, "reaclib"},
+        {ReactionType::REACLIB_WEAK, "reaclib_weak"},
         {ReactionType::LOGICAL_REACLIB, "logical_reaclib"},
+        {ReactionType::LOGICAL_REACLIB_WEAK, "logical_reaclib_weak"},
+
     };
 
     static std::unordered_map<ReactionType, std::string> ReactionPhysicalTypeNames = {
         {ReactionType::WEAK, "Weak"},
         {ReactionType::REACLIB, "Strong"},
         {ReactionType::LOGICAL_REACLIB, "Strong"},
+        {ReactionType::REACLIB_WEAK, "Weak"},
+        {ReactionType::LOGICAL_REACLIB_WEAK, "Weak"},
     };
     /**
      * @struct RateCoefficientSet
@@ -363,7 +370,7 @@ namespace gridfire::reaction {
          * @param reactants A vector of reactant species.
          * @param products A vector of product species.
          * @param qValue The Q-value of the reaction in MeV.
-         * @param label The source label for the rate data (e.g., "wc12", "st08").
+         * @param label The sources label for the rate data (e.g., "wc12", "st08").
          * @param sets The set of rate coefficients.
          * @param reverse True if this is a reverse reaction rate.
          */
@@ -645,6 +652,27 @@ namespace gridfire::reaction {
         }
     };
 
+    class WeakReaclibReaction final : public ReaclibReaction {
+    public:
+        using ReaclibReaction::ReaclibReaction;
+
+        [[nodiscard]] ReactionType type() const override { return ReactionType::REACLIB_WEAK; }
+
+        [[nodiscard]] std::unique_ptr<Reaction> clone() const override {
+            return std::make_unique<WeakReaclibReaction>(
+                m_id,
+                m_peName,
+                m_chapter,
+                reactants(),
+                products(),
+                m_qValue,
+                m_sourceLabel,
+                m_rateCoefficients,
+                m_reverse
+            );
+        }
+    };
+
 
     /**
      * @class LogicalReaclibReaction
@@ -664,7 +692,7 @@ namespace gridfire::reaction {
          * @param reactions A vector of reactions that represent the same logical process.
          * @throws std::runtime_error if the provided reactions have inconsistent Q-values.
          */
-        explicit LogicalReaclibReaction(const std::vector<ReaclibReaction> &reactions);
+        explicit LogicalReaclibReaction(const std::vector<std::unique_ptr<ReaclibReaction>> &reactions);
 
         /**
          * @breif Constructs a LogicalReaction from a vector of `Reaction` objects and allows the user
@@ -673,7 +701,7 @@ namespace gridfire::reaction {
          * @param reverse A flag to control if this logical reaction is reverse or not
          * @returns std::runtime_error if the provided reactions have inconsistent Q-values.
          */
-        explicit LogicalReaclibReaction(const std::vector<ReaclibReaction> &reactions, bool reverse);
+        explicit LogicalReaclibReaction(const std::vector<std::unique_ptr<ReaclibReaction>> &reactions, bool reverse);
 
         /**
          * @brief Adds another `Reaction` source to this logical reaction.
@@ -718,7 +746,12 @@ namespace gridfire::reaction {
             double Ye, double mue, const fourdst::composition::Composition& comp
         ) const override;
 
-        [[nodiscard]] ReactionType type() const override { return ReactionType::LOGICAL_REACLIB; }
+        [[nodiscard]] ReactionType type() const override {
+            if (m_weak) {
+                return ReactionType::LOGICAL_REACLIB_WEAK;
+            }
+            return ReactionType::LOGICAL_REACLIB;
+        }
 
         [[nodiscard]] std::unique_ptr<Reaction> clone() const override;
 
@@ -760,6 +793,7 @@ namespace gridfire::reaction {
     private:
         std::vector<std::string> m_sources; ///< List of source labels.
         std::vector<RateCoefficientSet> m_rates; ///< List of rate coefficient sets from each source.
+        bool m_weak = false;
 
     private:
         /**

src/include/gridfire/solver/strategies/CVODE_solver_strategy.h

@@ -232,10 +232,14 @@ namespace gridfire::solver {
             const std::vector<fourdst::atomic::Species>* networkSpecies{};
             std::unique_ptr<exceptions::EngineError> captured_exception = nullptr;
             std::optional<std::map<fourdst::atomic::Species, std::unordered_map<std::string, double>>> reaction_contribution_map;
+            double neutrino_energy_loss_rate = 0.0;
+            double total_neutrino_flux = 0.0;
         };
 
         struct CVODERHSOutputData {
             std::map<fourdst::atomic::Species, std::unordered_map<std::string, double>> reaction_contribution_map;
+            double neutrino_energy_loss_rate;
+            double total_neutrino_flux;
         };
 
     private:

src/include/gridfire/types/types.h

@@ -40,6 +40,8 @@ namespace gridfire {
         double energy; ///< Energy in ergs after evaluation
         double dEps_dT; ///< Partial derivative of energy generation rate with respect to temperature
         double dEps_dRho; ///< Partial derivative of energy generation rate with respect to density
+        double neutrino_energy_loss_rate; ///< Neutrino energy loss rate in ergs/g/s
+        double total_neutrino_flux; ///< Total neutrino flux in neutrinos/g/s
 
         friend std::ostream& operator<<(std::ostream& os, const NetOut& netOut) {
             os << "NetOut(composition=" << netOut.composition << ", num_steps=" << netOut.num_steps << ", ε=" << netOut.energy << ", dε/dT=" << netOut.dEps_dT << ", dε/dρ=" << netOut.dEps_dRho << ")";

src/lib/engine/engine_graph.cpp

@@ -32,6 +32,69 @@
 #include "cppad/utility/sparse_rcv.hpp"
 
 
+namespace {
+    enum class REACLIB_WEAK_TYPES {
+        BETA_PLUS_DECAY,
+        BETA_MINUS_DECAY,
+        ELECTRON_CAPTURE,
+        POSITRON_CAPTURE,
+        NONE
+    };
+
+    REACLIB_WEAK_TYPES get_weak_reaclib_reaction_type(const gridfire::reaction::Reaction& r) {
+        if (r.type() != gridfire::reaction::ReactionType::REACLIB_WEAK) {
+            return REACLIB_WEAK_TYPES::NONE;
+        }
+
+        // Get the () part of the id
+        const std::string_view id = r.id();
+        const size_t open_paren_pos = id.find('(');
+        const size_t close_paren_pos = id.find(')');
+        if (open_paren_pos == std::string_view::npos || close_paren_pos == std::string_view::npos || close_paren_pos <= open_paren_pos) {
+            throw gridfire::exceptions::ReactionParsingError("Invalid REACLIB weak reaction ID format.", std::string(id));
+        }
+
+        const std::string_view reaction_type_str = id.substr(open_paren_pos + 1, close_paren_pos - open_paren_pos - 1);
+
+        // Find the comma and extract the part to the left and right of it
+        const size_t reaction_type_pos = reaction_type_str.find(',');
+        if (reaction_type_pos == std::string_view::npos) {
+            throw gridfire::exceptions::ReactionParsingError("Invalid REACLIB weak reaction ID format: missing comma.", std::string(id));
+        }
+
+        const std::string_view projectiles_str = reaction_type_str.substr(0, reaction_type_pos);
+        const std::string_view ejectiles_str = reaction_type_str.substr(reaction_type_pos + 1);
+
+        // Check if the projectiles string has "e+" or "e-"
+        const bool has_captured_positron = (projectiles_str.find("e+") != std::string_view::npos);
+        const bool has_captured_electron = (projectiles_str.find("e-") != std::string_view::npos);
+
+        const bool has_ejected_electron = (ejectiles_str.find("e-") != std::string_view::npos);
+        const bool has_ejected_positron = (ejectiles_str.find("e+") != std::string_view::npos);
+
+        // Assert that only one of the four possibilities is true
+        const int true_count = static_cast<int>(has_captured_positron) +
+                         static_cast<int>(has_captured_electron) +
+                         static_cast<int>(has_ejected_electron) +
+                         static_cast<int>(has_ejected_positron);
+        if (true_count != 1) {
+            throw gridfire::exceptions::ReactionParsingError("Invalid REACLIB weak reaction ID format: must have exactly one of e+, e- in projectiles or ejectiles.", std::string(id));
+        }
+
+        if (has_ejected_positron) {
+            return REACLIB_WEAK_TYPES::BETA_PLUS_DECAY;
+        } if (has_ejected_electron) {
+            return REACLIB_WEAK_TYPES::BETA_MINUS_DECAY;
+        } if (has_captured_electron) {
+            return REACLIB_WEAK_TYPES::ELECTRON_CAPTURE;
+        } if (has_captured_positron) {
+            return REACLIB_WEAK_TYPES::POSITRON_CAPTURE;
+        }
+
+        return REACLIB_WEAK_TYPES::NONE;
+    }
+}
+
 namespace gridfire::engine {
     GraphEngine::GraphEngine(
         const fourdst::composition::Composition &composition,
@@ -82,7 +145,6 @@ namespace gridfire::engine {
     ) const {
         LOG_TRACE_L3(m_logger, "Calculating RHS and Energy in GraphEngine at T9 = {}, rho = {}.", T9, rho);
         const double Ye = comp.getElectronAbundance();
-        const double mue = 0.0; // TODO: Remove
         if (m_usePrecomputation) {
             LOG_TRACE_L3(m_logger, "Using precomputation for reaction rates in GraphEngine calculateRHSAndEnergy.");
             std::vector<double> bare_rates;
@@ -92,7 +154,7 @@ namespace gridfire::engine {
 
             for (const auto& reaction: activeReactions) {
                 assert(m_reactions.contains(*reaction)); // A bug which results in this failing indicates a serious internal inconsistency and should only be present during development.
-                bare_rates.push_back(reaction->calculate_rate(T9, rho, Ye, mue, comp.getMolarAbundanceVector(), m_indexToSpeciesMap));
+                bare_rates.push_back(reaction->calculate_rate(T9, rho, Ye, 0.0, comp.getMolarAbundanceVector(), m_indexToSpeciesMap));
                 if (reaction->type() != reaction::ReactionType::WEAK) {
                     bare_reverse_rates.push_back(calculateReverseRate(*reaction, T9, rho, comp));
                 }
@@ -109,7 +171,7 @@ namespace gridfire::engine {
                 T9,
                 rho,
                 Ye,
-                mue,
+                0.0,
                 [&comp](const fourdst::atomic::Species& species) -> std::optional<size_t> {
                     if (comp.contains(species)) {
                         return comp.getSpeciesIndex(species); // Return the index of the species in the composition
@@ -543,15 +605,6 @@ namespace gridfire::engine {
         fullNetIn.temperature = netIn.temperature;
         fullNetIn.density = netIn.density;
 
-        // Short circuit path if already primed
-        // if (m_has_been_primed) {
-        //     PrimingReport report;
-        //     report.primedComposition = composition;
-        //     report.success = true;
-        //     report.status = PrimingReportStatus::ALREADY_PRIMED;
-        //     return report;
-        // }
-
         std::optional<std::vector<reaction::ReactionType>> reactionTypesToIgnore = std::nullopt;
         if (!m_useReverseReactions) {
             reactionTypesToIgnore = {reaction::ReactionType::WEAK};
@@ -623,6 +676,7 @@ namespace gridfire::engine {
             rho
         );
 
+        StepDerivatives<double> result;
 
         // --- Optimized loop ---
         std::vector<double> molarReactionFlows;
@@ -672,6 +726,7 @@ namespace gridfire::engine {
                 throw exceptions::BadRHSEngineError("Non-finite forward molar reaction flow computed.");
             }
 
+
             // --- Reverse reaction flow ---
             // Only do this is the reaction has a non-zero reverse symmetry factor (i.e. is reversible)
             double reverseMolarReactionFlow = 0.0;
@@ -693,13 +748,37 @@ namespace gridfire::engine {
             }
 
             molarReactionFlows.push_back(forwardMolarReactionFlow - reverseMolarReactionFlow);
+
+            if (reaction->type() == reaction::ReactionType::REACLIB_WEAK) {
+                double q_abs = std::abs(reaction->qValue());
+                REACLIB_WEAK_TYPES weakType = get_weak_reaclib_reaction_type(*reaction);
+                double neutrino_loss_fraction = 0.0;
+                switch (weakType) {
+                    case REACLIB_WEAK_TYPES::BETA_PLUS_DECAY:
+                        [[fallthrough]];
+                    case REACLIB_WEAK_TYPES::BETA_MINUS_DECAY:
+                        neutrino_loss_fraction = 0.5; // Approximate 50% energy loss to neutrinos for beta decays
+                        break;
+                    case REACLIB_WEAK_TYPES::ELECTRON_CAPTURE:
+                        [[fallthrough]];
+                    case REACLIB_WEAK_TYPES::POSITRON_CAPTURE:
+                        neutrino_loss_fraction = 1.0;
+                        break;
+                    default: ;
+                }
+
+                double local_neutrino_loss =  molarReactionFlows.back() * q_abs * neutrino_loss_fraction * m_constants.Na * m_constants.MeV_to_erg;
+                double local_neutrino_flux = molarReactionFlows.back() * m_constants.Na;
+
+                result.totalNeutrinoFlux += local_neutrino_flux;
+                result.neutrinoEnergyLossRate += local_neutrino_loss;
+            }
             reactionCounter++;
         }
 
         LOG_TRACE_L3(m_logger, "Computed {} molar reaction flows for active reactions. Assembling these into RHS", molarReactionFlows.size());
 
         // --- Assemble molar abundance derivatives ---
-        StepDerivatives<double> result;
         for (const auto& species: m_networkSpecies) {
             result.dydt[species] = 0.0; // Initialize the change in abundance for each network species to 0
         }
@@ -724,34 +803,13 @@ namespace gridfire::engine {
             reactionCounter++;
         }
 
-        // std::vector<std::string> reactionIDs;
-        // for (const auto& reaction: activeReactions) {
-        //     reactionIDs.push_back(std::string(reaction->id()));
-        // }
-        //
-        // std::vector<std::unique_ptr<utils::ColumnBase>> columns;
-        // columns.push_back(std::make_unique<utils::Column<std::string>>("Reaction", reactionIDs));
-        // for (const auto& [species, contributions] : reactionContributions) {
-        //     std::vector<double> speciesData;
-        //     for (const auto& reactionID : reactionIDs) {
-        //         if (contributions.contains(reactionID)) {
-        //             speciesData.push_back(contributions.at(reactionID));
-        //         } else {
-        //             speciesData.push_back(0.0);
-        //         }
-        //     }
-        //     columns.push_back(std::make_unique<utils::Column<double>>(std::string(species.name()), speciesData));
-        // }
-
-        // utils::print_table("Contributions", columns);
-        // exit(0);
-
         // --- Calculate the nuclear energy generation rate ---
         double massProductionRate = 0.0; // [mol][s^-1]
         for (const auto & species : m_networkSpecies) {
             massProductionRate += result.dydt[species] * species.mass() * m_constants.u;
         }
         result.nuclearEnergyGenerationRate = -massProductionRate * m_constants.Na * m_constants.c * m_constants.c; // [erg][s^-1][g^-1]
+        result.nuclearEnergyGenerationRate -= result.neutrinoEnergyLossRate;
         return result;
 
     }
@@ -1130,7 +1188,7 @@ namespace gridfire::engine {
     ) const {
         const double Ye = comp.getElectronAbundance();
 
-        auto [dydt, _, __] = calculateAllDerivatives<double>(
+        auto result = calculateAllDerivatives<double>(
             comp.getMolarAbundanceVector(),
             T9,
             rho,
@@ -1146,6 +1204,8 @@ namespace gridfire::engine {
                 return activeReactions.contains(reaction);
             }
         );
+        const std::map<fourdst::atomic::Species, double>& dydt = result.dydt;
+
         std::unordered_map<fourdst::atomic::Species, double> speciesTimescales;
         speciesTimescales.reserve(m_networkSpecies.size());
         for (const auto& species : m_networkSpecies) {
@@ -1182,18 +1242,6 @@ namespace gridfire::engine {
             return std::nullopt; // Species not present
         };
 
-        auto [dydt, _, __] = calculateAllDerivatives<double>(
-            Y,
-            T9,
-            rho,
-            Ye,
-            0.0,
-            speciesLookup,
-            [&activeReactions](const reaction::Reaction& reaction) -> bool {
-                return activeReactions.contains(reaction);
-            }
-        );
-
         std::unordered_map<fourdst::atomic::Species, double> speciesDestructionTimescales;
         speciesDestructionTimescales.reserve(m_networkSpecies.size());
         for (const auto& species : m_networkSpecies) {
@@ -1273,7 +1321,7 @@ namespace gridfire::engine {
 
         // 5. Call the actual templated function
         // We let T9 and rho be constant, so we pass them as fixed values.
-        auto [dydt, nuclearEnergyGenerationRate, _] = calculateAllDerivatives<CppAD::AD<double>>(
+        auto result = calculateAllDerivatives<CppAD::AD<double>>(
             adY,
             adT9,
             adRho,
@@ -1289,15 +1337,15 @@ namespace gridfire::engine {
 
         // Extract the raw vector from the associative map
         std::vector<CppAD::AD<double>> dependentVector;
-        dependentVector.reserve(dydt.size() + 1);
+        dependentVector.reserve(result.dydt.size() + 1);
         std::ranges::transform(
-            dydt,
+            result.dydt,
             std::back_inserter(dependentVector),
             [](const auto& kv) {
                 return kv.second;
             }
         );
-        dependentVector.push_back(nuclearEnergyGenerationRate);
+        dependentVector.push_back(result.nuclearEnergyGenerationRate);
 
         m_rhsADFun.Dependent(adInput, dependentVector);
 

src/lib/engine/procedures/construction.cpp

@@ -18,32 +18,6 @@
 #include "quill/Logger.h"
 #include "quill/LogMacros.h"
 namespace {
-    // Simple heuristic to check if a reaclib reaction is a strong or weak reaction
-    /*  A weak reaction is defined here as one where:
-        - The number of reactants is equal to the number of products
-        - There is only one reactant and one product
-        - The mass number (A) of the reactant is equal to the mass number (A) of the product
-    */
-    bool reaclib_reaction_is_weak(const gridfire::reaction::Reaction& reaction) {
-        const std::vector<fourdst::atomic::Species>& reactants = reaction.reactants();
-        const std::vector<fourdst::atomic::Species>& products = reaction.products();
-
-        if (reactants.size() != products.size()) {
-            return false;
-        }
-
-        if (reactants.size() != 1 || products.size() != 1) {
-            return false;
-        }
-
-        if (std::floor(reactants[0].a()) != std::floor(products[0].a())) {
-            return false;
-        }
-
-        return true;
-    }
-
-
     gridfire::reaction::ReactionSet register_weak_reactions(
         const gridfire::rates::weak::WeakRateInterpolator &weakInterpolator,
         const gridfire::engine::NetworkConstructionFlags reactionTypes
@@ -114,7 +88,7 @@ namespace {
         if (has_flag(reaction_types, gridfire::engine::NetworkConstructionFlags::REACLIB_STRONG)) {
             const auto& allReaclibReactions = gridfire::reaclib::get_all_reaclib_reactions();
             for (const auto& reaction : allReaclibReactions) {
-                const bool isWeakReaction = reaclib_reaction_is_weak(*reaction);
+                const bool isWeakReaction = gridfire::reaclib::reaction_is_weak(*reaction);
                 const bool okayToUseReaclibWeakReaction = has_flag(reaction_types, gridfire::engine::NetworkConstructionFlags::REACLIB_WEAK);
 
                 const bool reaclibWeakOkay = !isWeakReaction || okayToUseReaclibWeakReaction;
@@ -126,10 +100,10 @@ namespace {
         return strong_reaction_pool;
     }
 
-    bool validate_unique_weak_set(gridfire::engine::NetworkConstructionFlags flag) {
+    bool validate_unique_weak_set(const gridfire::engine::NetworkConstructionFlags flag) {
         // This method ensures that weak reactions will only be fetched from either reaclib or the weak reaction library (WRL)
         // but not both
-        const std::array<gridfire::engine::NetworkConstructionFlags, 4> WRL_Flags = {
+        constexpr std::array<gridfire::engine::NetworkConstructionFlags, 4> WRL_Flags = {
             gridfire::engine::NetworkConstructionFlags::WRL_BETA_PLUS,
             gridfire::engine::NetworkConstructionFlags::WRL_ELECTRON_CAPTURE,
             gridfire::engine::NetworkConstructionFlags::WRL_POSITRON_CAPTURE,

src/lib/reaction/reaclib.cpp

@@ -223,4 +223,23 @@ namespace gridfire::reaclib {
         }
         return *s_all_reaclib_reactions_ptr;
     }
+
+    bool reaction_is_weak(const reaction::Reaction& reaction) {
+        const std::vector<fourdst::atomic::Species>& reactants = reaction.reactants();
+        const std::vector<fourdst::atomic::Species>& products = reaction.products();
+
+        if (reactants.size() != products.size()) {
+            return false;
+        }
+
+        if (reactants.size() != 1 || products.size() != 1) {
+            return false;
+        }
+
+        if (std::floor(reactants[0].a()) != std::floor(products[0].a())) {
+            return false;
+        }
+
+        return true;
+    }
 } // namespace gridfire::reaclib

src/lib/reaction/reaction.cpp

@@ -13,13 +13,15 @@
 #include "fourdst/atomic/atomicSpecies.h"
 
 #include "xxhash64.h"
+#include "gridfire/exceptions/exceptions.h"
+#include "gridfire/reaction/reaclib.h"
 
 namespace {
-    std::string_view safe_check_reactant_id(const std::vector<gridfire::reaction::ReaclibReaction>& reactions) {
+    std::string_view safe_check_reactant_id(const std::vector<std::unique_ptr<gridfire::reaction::ReaclibReaction>>& reactions) {
         if (reactions.empty()) {
             throw std::runtime_error("No reactions found in the REACLIB reaction set.");
         }
-        return reactions.front().peName();
+        return reactions.front()->peName();
     }
 }
 
@@ -192,40 +194,59 @@ namespace gridfire::reaction {
 
 
     LogicalReaclibReaction::LogicalReaclibReaction(
-        const std::vector<ReaclibReaction>& reactions
+        const std::vector<std::unique_ptr<ReaclibReaction>>& reactions
     ) : LogicalReaclibReaction(reactions, false) {}
 
     LogicalReaclibReaction::LogicalReaclibReaction(
-        const std::vector<ReaclibReaction> &reactions,
+        const std::vector<std::unique_ptr<ReaclibReaction>> &reactions,
         const bool reverse
     ) :
     ReaclibReaction(
         safe_check_reactant_id(reactions),
-        reactions.front().peName(),
+        reactions.front()->peName(),
-        reactions.front().chapter(),
+        reactions.front()->chapter(),
-        reactions.front().reactants(),
+        reactions.front()->reactants(),
-        reactions.front().products(),
+        reactions.front()->products(),
-        reactions.front().qValue(),
+        reactions.front()->qValue(),
-        reactions.front().sourceLabel(),
+        reactions.front()->sourceLabel(),
-        reactions.front().rateCoefficients(),
+        reactions.front()->rateCoefficients(),
         reverse)
     {
         m_sources.reserve(reactions.size());
         m_rates.reserve(reactions.size());
         for (const auto& reaction : reactions) {
-            if (std::abs(reaction.qValue() - m_qValue) > 1e-6) {
+            if (std::abs(reaction->qValue() - m_qValue) > 1e-6) {
                 LOG_ERROR(
                     m_logger,
                     "LogicalReaclibReaction constructed with reactions having different Q-values. Expected {} got {}.",
                     m_qValue,
-                    reaction.qValue()
+                    reaction->qValue()
                 );
                 m_logger -> flush_log();
-                throw std::runtime_error("LogicalReaclibReaction constructed with reactions having different Q-values. Expected " + std::to_string(m_qValue) + " got " + std::to_string(reaction.qValue()) + " (difference : " + std::to_string(std::abs(reaction.qValue() - m_qValue)) + ").");
+                throw std::runtime_error("LogicalReaclibReaction constructed with reactions having different Q-values. Expected " + std::to_string(m_qValue) + " got " + std::to_string(reaction->qValue()) + " (difference : " + std::to_string(std::abs(reaction->qValue() - m_qValue)) + ").");
             }
-            m_sources.emplace_back(reaction.sourceLabel());
+            m_sources.emplace_back(reaction->sourceLabel());
-            m_rates.push_back(reaction.rateCoefficients());
+            m_rates.push_back(reaction->rateCoefficients());
         }
+
+        std::unordered_set<bool> reaction_weak_types;
+        for (const auto& reaction : reactions) {
+            reaction_weak_types.insert(reaclib::reaction_is_weak(*reaction));
+        }
+
+        if (reaction_weak_types.size() != 1) {
+            LOG_ERROR(
+                m_logger,
+                "LogicalReaclibReaction constructed with reactions of mixed weak/strong types. Each LogicalReaclibReaction must contain only weak or only strong reactions."
+            );
+
+            throw exceptions::ReactionError(
+                "LogicalReaclibReaction constructed with reactions of mixed weak/strong types. Each LogicalReaclibReaction must contain only weak or only strong reactions.",
+                m_id
+            );
+        }
+
+        m_weak = *reaction_weak_types.begin();
     }
 
 
@@ -550,17 +571,21 @@ namespace gridfire::reaction {
     }
 
     ReactionSet packReactionSet(const ReactionSet& reactionSet) {
-        std::unordered_map<std::string, std::vector<ReaclibReaction>> groupedReaclibReactions;
+        std::unordered_map<std::string, std::vector<std::unique_ptr<ReaclibReaction>>> groupedReaclibReactions;
         ReactionSet finalReactionSet;
 
         for (const auto& reaction_ptr : reactionSet) {
             switch (reaction_ptr->type()) {
+                case ReactionType::REACLIB_WEAK:
+                    [[fallthrough]];
                 case ReactionType::REACLIB: {
                     const auto& reaclib_cast_reaction = static_cast<const ReaclibReaction&>(*reaction_ptr); // NOLINT(*-pro-type-static-cast-downcast)
                     std::string rid = std::format("{}{}", reaclib_cast_reaction.peName(), (reaclib_cast_reaction.is_reverse() ? "_r" : ""));
-                    groupedReaclibReactions[rid].push_back(reaclib_cast_reaction);
+                    groupedReaclibReactions[rid].push_back(std::make_unique<ReaclibReaction>(reaclib_cast_reaction));
                     break;
                 }
+                case ReactionType::LOGICAL_REACLIB_WEAK:
+                    [[fallthrough]];
                 case ReactionType::LOGICAL_REACLIB: {
                     // It doesn't make sense to pack an already-packed reaction.
                     throw std::runtime_error("packReactionSet: Cannot pack a LogicalReaclibReaction.");
@@ -573,22 +598,15 @@ namespace gridfire::reaction {
         }
 
         // Now, process the grouped REACLIB reactions
-        for (const auto &[key, reactionsGroup]: groupedReaclibReactions) {
+        for (const auto &reactionsGroup: groupedReaclibReactions | std::views::values) {
             if (reactionsGroup.empty()) {
                 continue;
             }
             if (reactionsGroup.size() == 1) {
-                finalReactionSet.add_reaction(reactionsGroup.front());
+                finalReactionSet.add_reaction(reactionsGroup.front()->clone());
             }
             else {
-                // Check that is_reverse is consistent across the group
+                const auto logicalReaction = std::make_unique<LogicalReaclibReaction>(reactionsGroup, reactionsGroup.front()->is_reverse());
-                assert(std::ranges::all_of(
-                    reactionsGroup,
-                    [&reactionsGroup](const ReaclibReaction& r) {
-                        return r.is_reverse() == reactionsGroup.front().is_reverse();
-                    }
-                ) && "Inconsistent is_reverse values in grouped REACLIB reactions.");
-                const auto logicalReaction = std::make_unique<LogicalReaclibReaction>(reactionsGroup, reactionsGroup.front().is_reverse());
                 finalReactionSet.add_reaction(logicalReaction->clone());
             }
         }

src/lib/solver/strategies/CVODE_solver_strategy.cpp

@@ -147,6 +147,9 @@ namespace gridfire::solver {
         m_num_steps = 0;
         double accumulated_energy = 0.0;
 
+        double accumulated_neutrino_energy_loss = 0.0;
+        double accumulated_total_neutrino_flux = 0.0;
+
         size_t total_convergence_failures = 0;
         size_t total_nonlinear_iterations = 0;
         size_t total_update_stages_triggered = 0;
@@ -187,6 +190,8 @@ namespace gridfire::solver {
 
             sunrealtype* y_data = N_VGetArrayPointer(m_Y);
             const double current_energy = y_data[numSpecies]; // Specific energy rate
+            accumulated_neutrino_energy_loss += user_data.neutrino_energy_loss_rate;
+            accumulated_total_neutrino_flux += user_data.total_neutrino_flux;
             size_t iter_diff = (total_nonlinear_iterations + nliters) - prev_nonlinear_iterations;
             size_t convFail_diff = (total_convergence_failures + nlcfails) - prev_convergence_failures;
             if (m_stdout_logging_enabled) {
@@ -503,6 +508,8 @@ namespace gridfire::solver {
 
         netOut.dEps_dT = dEps_dT;
         netOut.dEps_dRho = dEps_dRho;
+        netOut.neutrino_energy_loss_rate = accumulated_neutrino_energy_loss;
+        netOut.total_neutrino_flux = accumulated_total_neutrino_flux;
         LOG_TRACE_L2(m_logger, "Output data built!");
         LOG_TRACE_L2(m_logger, "Solver evaluation complete!.");
 
@@ -561,8 +568,11 @@ namespace gridfire::solver {
 
         try {
             LOG_TRACE_L2(instance->m_logger, "CVODE RHS wrapper called at time {}", t);
-            const auto [reaction_contribution_map] = instance->calculate_rhs(t, y, ydot, data);
+            // ReSharper disable once CppUseStructuredBinding
-            data->reaction_contribution_map = reaction_contribution_map;
+            const auto result = instance->calculate_rhs(t, y, ydot, data);
+            data->reaction_contribution_map = result.reaction_contribution_map;
+            data->neutrino_energy_loss_rate = result.neutrino_energy_loss_rate;
+            data->total_neutrino_flux = result.total_neutrino_flux;
             LOG_TRACE_L2(instance->m_logger, "CVODE RHS wrapper completed successfully at time {}", t);
             return 0;
         } catch (const exceptions::EngineError& e) {
@@ -724,7 +734,14 @@ namespace gridfire::solver {
         }
 
         sunrealtype* ydot_data = N_VGetArrayPointer(ydot);
-        const auto& [dydt, nuclearEnergyGenerationRate, reactionContributions] = result.value();
+        const auto& [
+            dydt,
+            nuclearEnergyGenerationRate,
+            reactionContributions,
+            neutrinoEnergyLossRate,
+            totalNeutrinoFlux
+        ] = result.value();
+
         LOG_TRACE_L2(m_logger, "Done calculating RHS at time {}, specific nuclear energy generation rate: {}", t, nuclearEnergyGenerationRate);
         LOG_TRACE_L2(m_logger, "RHS at time {} is {}", t, [&dydt]() -> std::string {
             std::stringstream ss;
@@ -745,7 +762,7 @@ namespace gridfire::solver {
         }
         ydot_data[numSpecies] = nuclearEnergyGenerationRate; // Set the last element to the specific energy rate
 
-        return {reactionContributions};
+        return {reactionContributions, result.value().neutrinoEnergyLossRate, result.value().totalNeutrinoFlux};
     }
 
     void CVODESolverStrategy::initialize_cvode_integration_resources(

src/python/reaction/bindings.cpp

@@ -183,12 +183,12 @@ void register_reaction_bindings(py::module &m) {
 
     py::class_<gridfire::reaction::LogicalReaclibReaction, gridfire::reaction::ReaclibReaction>(m, "LogicalReaclibReaction")
         .def(
-            py::init<const std::vector<gridfire::reaction::ReaclibReaction>>(),
+            py::init<const std::vector<std::unique_ptr<gridfire::reaction::ReaclibReaction>>>(),
              py::arg("reactions"),
              "Construct a LogicalReaclibReaction from a vector of ReaclibReaction objects."
         )
         .def(
-            py::init<const std::vector<gridfire::reaction::ReaclibReaction>, bool>(),
+            py::init<const std::vector<std::unique_ptr<gridfire::reaction::ReaclibReaction>>, bool>(),
              py::arg("reactions"),
              py::arg("is_reverse"),
              "Construct a LogicalReaclibReaction from a vector of ReaclibReaction objects."

tests/graphnet_sandbox/main.cpp

@@ -95,6 +95,16 @@ void log_results(const gridfire::NetOut& netOut, const gridfire::NetIn& netIn) {
     delta.push_back(netOut.dEps_dRho);
     fractional.push_back(0.0);
 
+    initial.push_back(0.0);
+    final.push_back(netOut.neutrino_energy_loss_rate);
+    delta.push_back(netOut.neutrino_energy_loss_rate);
+    fractional.push_back(0.0);
+
+    initial.push_back(0.0);
+    final.push_back(netOut.total_neutrino_flux);
+    delta.push_back(netOut.total_neutrino_flux);
+    fractional.push_back(0.0);
+
     initial.push_back(netIn.composition.getMeanParticleMass());
     final.push_back(netOut.composition.getMeanParticleMass());
     delta.push_back(final.back() - initial.back());
@@ -108,6 +118,8 @@ void log_results(const gridfire::NetOut& netOut, const gridfire::NetIn& netIn) {
         labels.push_back("ε");
         labels.push_back("dε/dT");
         labels.push_back("dε/dρ");
+        labels.push_back("Eν");
+        labels.push_back("Fν");
         labels.push_back("<μ>");
         return labels;
     }();