fix(LogicalReaclibReaction): fixed bug where reaclib reverse reactions were being included in forward reaction sets erroneously

Previously, due to a indexing issue, reverse rates were sometimes included in the forward rates for strong reactions. This lead to erroneously high molar reaction flows for some reactions. This has been resolved so now each logical reaction set is guareteed to include only either forward reactions or reverse reactions (not both)

src/include/gridfire/reaction/reaction.h

@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ranges>
 #include <string_view>
 
 #include "fourdst/composition/atomicSpecies.h"
@@ -613,12 +614,22 @@ namespace gridfire::reaction {
     class LogicalReaclibReaction final : public ReaclibReaction {
     public:
         /**
-         * @brief Constructs a LogicalReaction from a vector of `Reaction` objects.
+         * @brief Constructs a LogicalReaction from a vector of `Reaction` objects. Implicitly assumes that the
+         * logical reaction is for a forward (i.e. not reverse) 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);
 
+        /**
+         * @breif Constructs a LogicalReaction from a vector of `Reaction` objects and allows the user
+         * to specify if the logical set is for a reverse reaction explicitly
+         * @param reactions A vector of reactions that represent the same logical process
+         * @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);
+
         /**
          * @brief Adds another `Reaction` source to this logical reaction.
          * @param reaction The reaction to add.
@@ -719,7 +730,7 @@ namespace gridfire::reaction {
             const T T953 = CppAD::pow(T9, 5.0/3.0);
             const T logT9 = CppAD::log(T9);
             // ReSharper disable once CppUseStructuredBinding
-            for (const auto& rate : m_rates) {
+            for (const auto& [rate, source] : std::views::zip(m_rates, m_sources)) {
                 const T exponent = rate.a0 +
                        rate.a1 / T9 +
                        rate.a2 / T913 +

src/lib/engine/engine_graph.cpp

@@ -854,9 +854,18 @@ namespace gridfire {
         // --- Pack the input variables into a vector for CppAD ---
         const size_t numSpecies = m_networkSpecies.size();
         std::vector<double> x(numSpecies + 2, 0.0);
-        const std::vector<double>& Y_dynamic = comp.getMolarAbundanceVector();
+        // const std::vector<double>& Y_dynamic = comp.getMolarAbundanceVector();
-        for (size_t i = 0; i < numSpecies; ++i) {
+        // for (size_t i = 0; i < numSpecies; ++i) {
-           x[i] = Y_dynamic[i];
+        //    x[i] = Y_dynamic[i];
+        // }
+        size_t i = 0;
+        for (const auto& species: m_networkSpecies) {
+            double Yi = 0.0;
+            if (comp.contains(species)) {
+                Yi = comp.getMolarAbundance(species);
+            }
+            x[i] = Yi;
+            i++;
         }
         x[numSpecies] = T9;
         x[numSpecies + 1] = rho;

src/lib/engine/views/engine_multiscale.cpp

@@ -225,9 +225,21 @@ namespace gridfire {
         const bool activeSpeciesIsSubset = std::ranges::any_of(activeSpecies, [&](const auto& species) -> bool {
             return !involvesSpecies(species);
         });
-        if (activeSpeciesIsSubset) {
+        if (!activeSpeciesIsSubset) {
-            LOG_CRITICAL(m_logger, "Active species contains species not in the network partition. Cannot generate Jacobian matrix for active species.");
+            std::string msg = std::format(
-            throw std::runtime_error("Active species contains species not in the network partition. Cannot generate Jacobian matrix for active species.");
+                "Active species set contains species ({}) not present in network partition. Cannot generate jacobian matrix due to this.",
+                [&]() -> std::string {
+                    std::stringstream ss;
+                    for (const auto& species : activeSpecies) {
+                        if (!this->involvesSpecies(species)) {
+                            ss << species << " ";
+                        }
+                    }
+                    return ss.str();
+                }()
+            );
+            LOG_CRITICAL(m_logger, "{}", msg);
+            throw std::runtime_error(msg);
         }
 
         std::vector<Species> dynamicActiveSpeciesIntersection;

src/lib/reaction/reaclib.cpp

@@ -90,7 +90,6 @@ namespace gridfire::reaclib {
 
         for (size_t i = 0; i < num_reactions; ++i) {
             const auto&[chapter, qValue, coeffs, reverse, label, rpName, reactants_str, products_str] = records[i];
-
             // The char arrays from the binary are not guaranteed to be null-terminated
             // if the string fills the entire buffer. We create null-terminated string_views.
             const std::string_view label_sv(label, strnlen(label, sizeof(label)));

src/lib/reaction/reaction.cpp

@@ -182,22 +182,29 @@ namespace gridfire::reaction {
     }
 
 
-    LogicalReaclibReaction::LogicalReaclibReaction(const std::vector<ReaclibReaction>& reactants) :
+    LogicalReaclibReaction::LogicalReaclibReaction(
-        ReaclibReaction(
+        const std::vector<ReaclibReaction>& reactions
-                 safe_check_reactant_id(reactants), // Use this first to check if the reactants array is empty and safely exit if so
+    ) : LogicalReaclibReaction(reactions, false) {}
-                 reactants.front().peName(),
-                 reactants.front().chapter(),
-                 reactants.front().reactants(),
-                 reactants.front().products(),
-                 reactants.front().qValue(),
-                 reactants.front().sourceLabel(),
-                 reactants.front().rateCoefficients(),
-                 reactants.front().is_reverse()) {
 
-        m_sources.reserve(reactants.size());
+    LogicalReaclibReaction::LogicalReaclibReaction(
-        m_rates.reserve(reactants.size());
+        const std::vector<ReaclibReaction> &reactions,
-        for (const auto& reaction : reactants) {
+        const bool reverse
-            if (std::abs(std::abs(reaction.qValue()) - std::abs(m_qValue)) > 1e-6) {
+    ) :
+    ReaclibReaction(
+        safe_check_reactant_id(reactions),
+        reactions.front().peName(),
+        reactions.front().chapter(),
+        reactions.front().reactants(),
+        reactions.front().products(),
+        reactions.front().qValue(),
+        reactions.front().sourceLabel(),
+        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) {
                 LOG_ERROR(
                     m_logger,
                     "LogicalReaclibReaction constructed with reactions having different Q-values. Expected {} got {}.",
@@ -212,6 +219,7 @@ namespace gridfire::reaction {
         }
     }
 
+
     void LogicalReaclibReaction::add_reaction(const ReaclibReaction& reaction) {
         if (reaction.peName() != m_id) {
             LOG_ERROR(m_logger, "Cannot add reaction with different peName to LogicalReaclibReaction. Expected {} got {}.", m_id, reaction.peName());
@@ -521,7 +529,8 @@ namespace gridfire::reaction {
             switch (reaction_ptr->type()) {
                 case ReactionType::REACLIB: {
                     const auto& reaclib_cast_reaction = static_cast<const ReaclibReaction&>(*reaction_ptr); // NOLINT(*-pro-type-static-cast-downcast)
-                    groupedReaclibReactions[std::string(reaclib_cast_reaction.peName())].push_back(reaclib_cast_reaction);
+                    std::string rid = std::format("{}{}", reaclib_cast_reaction.peName(), (reaclib_cast_reaction.is_reverse() ? "_r" : ""));
+                    groupedReaclibReactions[rid].push_back(reaclib_cast_reaction);
                     break;
                 }
                 case ReactionType::LOGICAL_REACLIB: {

src/lib/reaction/weak/weak.cpp

@@ -21,6 +21,13 @@ namespace {
         {fourdst::atomic::SpeciesErrorType::SPECIES_SYMBOL_NOT_FOUND, "Species symbol not found ((A,Z) out of range)"}
     };
 
+    std::string normalize_species_id(const fourdst::atomic::Species& species) {
+        auto result = std::string(species.name());
+        std::ranges::transform(result, result.begin(), ::tolower);
+        std::erase(result, '-');
+        return result;
+    }
+
     fourdst::atomic::Species resolve_weak_product(
         const gridfire::rates::weak::WeakReactionType type,
         const fourdst::atomic::Species& reactant
@@ -68,18 +75,20 @@ namespace {
         using namespace gridfire::rates::weak;
 
         std::string id;
+        std::string reactant_id = normalize_species_id(reactant);
+        std::string product_id = normalize_species_id(product);
         switch (type) {
             case WeakReactionType::BETA_MINUS_DECAY:
-                id = std::format("{}(,ν|)e-,{}", reactant.name(), product.name());
+                id = std::format("{}(,ν|)e-,{}", reactant_id, product_id);
                 break;
             case WeakReactionType::BETA_PLUS_DECAY:
-                id = std::format("{}(,ν)e+,{}", reactant.name(), product.name());
+                id = std::format("{}(,ν)e+,{}", reactant_id, product_id);
                 break;
             case WeakReactionType::ELECTRON_CAPTURE:
-                id = std::format("{}(e-,ν){}", reactant.name(), product.name());
+                id = std::format("{}(e-,ν){}", reactant_id, product_id);
                 break;
             case WeakReactionType::POSITRON_CAPTURE:
-                id = std::format("{}(e+,ν|){}", reactant.name(), product.name());
+                id = std::format("{}(e+,ν|){}", reactant_id, product_id);
                 break;
         }
         return id;

src/lib/solver/strategies/CVODE_solver_strategy.cpp

@@ -22,8 +22,6 @@
 #include "gridfire/trigger/procedures/trigger_pprint.h"
 #include "gridfire/utils/general_composition.h"
 
-static size_t s_call_counter = 0;
-
 namespace {
     std::unordered_map<int, std::string> cvode_ret_code_map {
         {0, "CV_SUCCESS: The solver succeeded."},
@@ -57,7 +55,7 @@ namespace {
         {-28, "CV_VECTOROP_ERR: A vector operation failed."},
         {-29, "CV_PROJ_MEM_NULL: The projection memory structure is NULL."},
         {-30, "CV_PROJFUNC_FAIL: The projection function failed in an unrecoverable manner."},
-        {-31, "CV_REPTD_PROJFUNC_ERR: THe projection function has repeated recoverable errors."}
+        {-31, "CV_REPTD_PROJFUNC_ERR: The projection function has repeated recoverable errors."}
     };
     void check_cvode_flag(const int flag, const std::string& func_name) {
         if (flag < 0) {
@@ -74,7 +72,7 @@ namespace {
 #elif SUNDIALS_HAVE_PTHREADS
         N_Vector vec = N_VNew_Pthreads(size, sun_ctx);
 #else
-        N_Vector vec = N_VNew_Serial(static_cast<long long>(size), sun_ctx);
+        const N_Vector vec = N_VNew_Serial(static_cast<long long>(size), sun_ctx);
 #endif
         check_cvode_flag(vec == nullptr ? -1 : 0, "N_VNew");
         return vec;
@@ -129,7 +127,7 @@ namespace gridfire::solver {
     }
 
     CVODESolverStrategy::~CVODESolverStrategy() {
-        std::cout << "Cleaning up CVODE resources..." << std::endl;
+        LOG_TRACE_L1(m_logger, "Cleaning up CVODE resources...");
         cleanup_cvode_resources(true);
 
         if (m_sun_ctx) {
@@ -138,10 +136,10 @@ namespace gridfire::solver {
     }
 
     NetOut CVODESolverStrategy::evaluate(const NetIn& netIn) {
-        LOG_TRACE_L2(m_logger, "Starting solver evaluation with T9: {} and rho: {}", netIn.temperature/1e9, netIn.density);
+        LOG_TRACE_L1(m_logger, "Starting solver evaluation with T9: {} and rho: {}", netIn.temperature/1e9, netIn.density);
-        LOG_TRACE_L2(m_logger, "Building engine update trigger....");
+        LOG_TRACE_L1(m_logger, "Building engine update trigger....");
         auto trigger = trigger::solver::CVODE::makeEnginePartitioningTrigger(1e12, 1e10, 1, true, 10);
-        LOG_TRACE_L2(m_logger, "Engine update trigger built!");
+        LOG_TRACE_L1(m_logger, "Engine update trigger built!");
 
 
         const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9)
@@ -149,15 +147,15 @@ namespace gridfire::solver {
         const auto absTol = m_config.get<double>("gridfire:solver:CVODESolverStrategy:absTol", 1.0e-8);
         const auto relTol = m_config.get<double>("gridfire:solver:CVODESolverStrategy:relTol", 1.0e-8);
 
-        LOG_TRACE_L2(m_logger, "Starting engine update chain...");
+        LOG_TRACE_L1(m_logger, "Starting engine update chain...");
         fourdst::composition::Composition equilibratedComposition = m_engine.update(netIn);
-        LOG_TRACE_L2(m_logger, "Engine updated and equilibrated composition found!");
+        LOG_TRACE_L1(m_logger, "Engine updated and equilibrated composition found!");
 
         size_t numSpecies = m_engine.getNetworkSpecies().size();
         uint64_t N = numSpecies + 1;
 
-        LOG_TRACE_L2(m_logger, "Number of species: {}", N);
+        LOG_TRACE_L1(m_logger, "Number of species: {}", N);
-        LOG_TRACE_L2(m_logger, "Initializing CVODE resources");
+        LOG_TRACE_L1(m_logger, "Initializing CVODE resources");
         m_cvode_mem = CVodeCreate(CV_BDF, m_sun_ctx);
         check_cvode_flag(m_cvode_mem == nullptr ? -1 : 0, "CVodeCreate");
 
@@ -166,7 +164,7 @@ namespace gridfire::solver {
         CVODEUserData user_data;
         user_data.solver_instance = this;
         user_data.engine = &m_engine;
-        LOG_TRACE_L2(m_logger, "CVODE resources successfully initialized!");
+        LOG_TRACE_L1(m_logger, "CVODE resources successfully initialized!");
 
         double current_time = 0;
         // ReSharper disable once CppTooWideScope
@@ -174,7 +172,7 @@ namespace gridfire::solver {
         m_num_steps = 0;
         double accumulated_energy = 0.0;
         int total_update_stages_triggered = 0;
-        LOG_TRACE_L2(m_logger, "Starting CVODE iteration");
+        LOG_TRACE_L1(m_logger, "Starting CVODE iteration");
         while (current_time < netIn.tMax) {
             user_data.T9 = T9;
             user_data.rho = netIn.density;
@@ -194,6 +192,7 @@ namespace gridfire::solver {
                 std::rethrow_exception(std::make_exception_ptr(*user_data.captured_exception));
             }
 
+            log_step_diagnostics(user_data);
             check_cvode_flag(flag, "CVode");
 
             long int n_steps;
@@ -359,9 +358,9 @@ namespace gridfire::solver {
     }
 
     int CVODESolverStrategy::cvode_rhs_wrapper(
-        sunrealtype t,
+        const sunrealtype t,
-        N_Vector y,
+        const N_Vector y,
-        N_Vector ydot,
+        const N_Vector ydot,
         void *user_data
     ) {
         auto* data = static_cast<CVODEUserData*>(user_data);
@@ -454,17 +453,9 @@ namespace gridfire::solver {
         sunrealtype* ydot_data = N_VGetArrayPointer(ydot);
         const auto& [dydt, nuclearEnergyGenerationRate] = result.value();
 
-        std::cout << "=========================\n";
         for (size_t i = 0; i < numSpecies; ++i) {
             fourdst::atomic::Species species = m_engine.getNetworkSpecies()[i];
             ydot_data[i] = dydt.at(species);
-            std::cout << "\t" << species << " dydt = " << dydt.at(species) << "\n";
-        }
-        std::cout << "\tε = " << nuclearEnergyGenerationRate << std::endl;
-        s_call_counter++;
-        std::cout << "\tMethod called " << s_call_counter << " times" << std::endl;
-        if (s_call_counter == 31) {
-            std::cout << "Last reliable step\n";
         }
         ydot_data[numSpecies] = nuclearEnergyGenerationRate; // Set the last element to the specific energy rate
     }
@@ -575,6 +566,8 @@ namespace gridfire::solver {
             0.0
         );
 
+        std::vector<double> M;
+        M.reserve(num_components);
         for (size_t i = 0; i < num_components - 1; i++) {
             const double weight = relTol * std::abs(y_data[i]) + absTol;
             if (weight == 0.0) continue; // Skip components with zero weight
@@ -583,8 +576,10 @@ namespace gridfire::solver {
 
             err_ratios[i] = err_ratio;
             speciesNames.emplace_back(user_data.networkSpecies->at(i).name());
+            M.push_back(user_data.networkSpecies->at(i).mass());
         }
-        fourdst::composition::Composition composition(speciesNames, Y_full);
+        std::vector<double> X = utils::massFractionFromMolarAbundanceAndMolarMass(Y_full, M);
+        fourdst::composition::Composition composition(speciesNames, X);
 
         if (err_ratios.empty()) {
             return;
@@ -620,9 +615,11 @@ namespace gridfire::solver {
         columns.push_back(std::make_unique<utils::Column<double>>("Error Ratio", sorted_err_ratios));
 
         std::cout << utils::format_table("Species Error Ratios", columns) << std::endl;
+
         diagnostics::inspect_jacobian_stiffness(*user_data.engine, composition, user_data.T9, user_data.rho);
-        diagnostics::inspect_species_balance(*user_data.engine, "N-14", composition, user_data.T9, user_data.rho);
+        for (const auto& species : sorted_speciesNames) {
-        diagnostics::inspect_species_balance(*user_data.engine, "n-1", composition, user_data.T9, user_data.rho);
+            diagnostics::inspect_species_balance(*user_data.engine, species, composition, user_data.T9, user_data.rho);
+        }
 
     }
 }