feat(dynamic-engine): added derivitves for energy generation rate. dε/dT and dε/dρ have been added to NetOut and computed with auto diff

meson.build

@@ -79,9 +79,19 @@ subdir('build-config')
 
 subdir('src')
 
-subdir('build-python')
+if get_option('build-python')
+    message('Configuring Python bindings...')
+    subdir('src-pybind')
+else
+    message('Skipping Python bindings...')
+endif
 
-subdir('tests')
+if get_option('build-tests')
+    message('Setting up tests for GridFire...')
+    subdir('tests')
+else
+    message('Skipping tests for GridFire...')
+endif
 
 
 if get_option('pkg-config')

meson_options.txt

@@ -1,2 +1,4 @@
 option('log_level', type: 'combo', choices: ['traceL3', 'traceL2', 'traceL1', 'debug', 'info', 'warning', 'error', 'critial'], value: 'info', description: 'Set the log level for the GridFire library')
 option('pkg-config', type: 'boolean', value: true, description: 'generate pkg-config file for GridFire (gridfire.pc)')
+option('build-python', type: 'boolean', value: true, description: 'build the python bindings so you can use GridFire from python')
+option('build-tests', type: 'boolean', value: true, description: 'build the test suite')

src/include/gridfire/engine/diagnostics/dynamic_engine_diagnostics.h

@@ -0,0 +1,34 @@
+#pragma once
+
+#include "gridfire/engine/engine_abstract.h"
+
+#include <vector>
+#include <string>
+
+namespace gridfire::diagnostics {
+    void report_limiting_species(
+        const DynamicEngine& engine,
+        const std::vector<double>& Y_full,
+        const std::vector<double>& E_full,
+        const std::vector<double>& dydt_full,
+        double relTol,
+        double absTol,
+        size_t top_n = 10
+    );
+
+    void inspect_species_balance(
+        const DynamicEngine& engine,
+        const std::string& species_name,
+        const std::vector<double>& Y_full,
+        double T9,
+        double rho
+    );
+
+    void inspect_jacobian_stiffness(
+        const DynamicEngine& engine,
+        const std::vector<double>& Y_full,
+        double T9,
+        double rho
+    );
+
+}

src/include/gridfire/engine/engine_abstract.h

@@ -64,6 +64,16 @@ namespace gridfire {
 
     using SparsityPattern = std::vector<std::pair<size_t, size_t>>;
 
+    struct EnergyDerivatives {
+        double dEps_dT = 0.0;  ///< Partial derivative of energy generation rate with respect to temperature.
+        double dEps_dRho = 0.0;///< Partial derivative of energy generation rate with respect to density.
+
+        friend std::ostream& operator<<(std::ostream& os, const EnergyDerivatives& ed) {
+            os << "<dε/dT: " << ed.dEps_dT << ", dε/dρ: " << ed.dEps_dRho << ">";
+            return os;
+        }
+    };
+
     /**
      * @brief Abstract base class for a reaction network engine.
      *
@@ -210,6 +220,23 @@ namespace gridfire {
             double rho
         ) const = 0;
 
+        /**
+         * @brief Calculate the derivatives of the energy generation rate with respect to T and rho.
+         *
+         * @param Y Vector of current abundances.
+         * @param T9 Temperature in units of 10^9 K.
+         * @param rho Density in g/cm^3.
+         * @return EnergyDerivatives containing dEps/dT and dEps/dRho.
+         *
+         * This method computes the partial derivatives of the specific nuclear energy
+         * generation rate with respect to temperature and density for the current state.
+         */
+        [[nodiscard]] virtual EnergyDerivatives calculateEpsDerivatives(
+            const std::vector<double>& Y,
+            const double T9,
+            const double rho
+        ) const = 0;
+
         /**
          * @brief Get the set of logical reactions in the network.
          *

src/include/gridfire/engine/engine_graph.h

@@ -68,6 +68,7 @@ namespace gridfire {
     static constexpr double MIN_JACOBIAN_THRESHOLD = 1e-24;
 
 
+
     /**
      * @class GraphEngine
      * @brief A reaction network engine that uses a graph-based representation.
@@ -143,6 +144,12 @@ namespace gridfire {
             const double rho
         ) const override;
 
+        [[nodiscard]] EnergyDerivatives calculateEpsDerivatives(
+            const std::vector<double>& Y,
+            const double T9,
+            const double rho
+        ) const override;
+
         /**
          * @brief Generates the Jacobian matrix for the current state.
          *
@@ -581,6 +588,7 @@ namespace gridfire {
 
         mutable boost::numeric::ublas::compressed_matrix<double> m_jacobianMatrix; ///< Jacobian matrix (species x species).
         mutable CppAD::ADFun<double> m_rhsADFun; ///< CppAD function for the right-hand side of the ODE.
+        mutable CppAD::ADFun<double> m_epsADFun; ///< CppAD function for the energy generation rate.
         mutable CppAD::sparse_jac_work m_jac_work; ///< Work object for sparse Jacobian calculations.
         CppAD::sparse_rc<std::vector<size_t>> m_full_jacobian_sparsity_pattern; ///< Full sparsity pattern for the Jacobian matrix.
 
@@ -652,6 +660,8 @@ namespace gridfire {
          */
         void recordADTape() const;
 
+        void recordEpsADTape() const;
+
         void collectAtomicReverseRateAtomicBases();
 
         void precomputeNetwork();

src/include/gridfire/engine/views/engine_adaptive.h

@@ -107,6 +107,12 @@ namespace gridfire {
             const double rho
         ) const override;
 
+        [[nodiscard]] EnergyDerivatives calculateEpsDerivatives(
+            const std::vector<double> &Y_culled,
+            const double T9,
+            const double rho
+        ) const override;
+
         /**
          * @brief Generates the Jacobian matrix for the active species.
          *

src/include/gridfire/engine/views/engine_defined.h

@@ -42,6 +42,13 @@ namespace gridfire{
             const double T9,
             const double rho
         ) const override;
+
+        [[nodiscard]] EnergyDerivatives calculateEpsDerivatives(
+            const std::vector<double> &Y,
+            const double T9,
+            const double rho
+        ) const override;
+
         /**
          * @brief Generates the Jacobian matrix for the active species.
          *

src/include/gridfire/engine/views/engine_multiscale.h

@@ -236,6 +236,12 @@ namespace gridfire {
             double rho
         ) const override;
 
+        [[nodiscard]] EnergyDerivatives calculateEpsDerivatives(
+            const std::vector<double> &Y,
+            const double T9,
+            const double rho
+        ) const override;
+
         /**
          * @brief Generates the Jacobian matrix for the current state.
          *
@@ -754,6 +760,10 @@ namespace gridfire {
              * @return True if the sets of species indices are not identical.
              */
             bool operator!=(const QSEGroup& other) const;
+
+            std::string toString() const;
+
+            std::string toString(DynamicEngine &engine) const;
         };
 
         /**
@@ -1056,9 +1066,11 @@ namespace gridfire {
          *      flux exceeds a configurable threshold, the group is considered valid and is added
          *      to the returned vector.
          */
-        std::vector<QSEGroup> validateGroupsWithFluxAnalysis(
+        std::pair<std::vector<MultiscalePartitioningEngineView::QSEGroup>, std::vector<MultiscalePartitioningEngineView
+        ::
+        QSEGroup>> validateGroupsWithFluxAnalysis(
             const std::vector<QSEGroup> &candidate_groups,
-            const std::vector<double>& Y,
+            const std::vector<double> &Y,
             double T9,
             double rho
         ) const;

src/include/gridfire/network.h

@@ -57,16 +57,18 @@ namespace gridfire {
         double energy; ///< Energy in ergs
         double culling = 0.0; ///< Culling threshold for reactions (default is 0.0, meaning no culling)
 
-        std::vector<double> MolarAbundance() const;
+        [[nodiscard]] std::vector<double> MolarAbundance() const;
     };
 
     struct NetOut {
         fourdst::composition::Composition composition; ///< Composition of the network after evaluation
         int num_steps; ///< Number of steps taken in the evaluation
         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
 
         friend std::ostream& operator<<(std::ostream& os, const NetOut& netOut) {
-            os << "NetOut(composition=" << netOut.composition << ", num_steps=" << netOut.num_steps << ", energy=" << netOut.energy << ")";
+            os << "NetOut(composition=" << netOut.composition << ", num_steps=" << netOut.num_steps << ", ε=" << netOut.energy << ", dε/dT=" << netOut.dEps_dT << ", dε/dρ=" << netOut.dEps_dRho << ")";
             return os;
         }
     };

src/include/gridfire/partition/partition_abstract.h

@@ -4,7 +4,6 @@
 #include <memory>
 
 namespace gridfire::partition {
-
     /**
      * @class PartitionFunction
      * @brief Abstract interface for evaluating nuclear partition functions.

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

@@ -98,16 +98,18 @@ namespace gridfire::solver {
             DynamicEngine* engine;
             double T9;
             double rho;
+            double energy;
             const std::vector<fourdst::atomic::Species>* networkSpecies;
             std::unique_ptr<exceptions::StaleEngineTrigger> captured_exception = nullptr;
         };
 
     private:
         Config& m_config = Config::getInstance();
+        quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
         static int cvode_rhs_wrapper(sunrealtype t, N_Vector y, N_Vector ydot, void *user_data);
         static int cvode_jac_wrapper(sunrealtype t, N_Vector y, N_Vector ydot, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
 
-        int calculate_rhs(sunrealtype t, N_Vector y, N_Vector ydot, const CVODEUserData* data) const;
+        void calculate_rhs(sunrealtype t, N_Vector y, N_Vector ydot, const CVODEUserData* data) const;
 
         void initialize_cvode_integration_resources(
             uint64_t N,
@@ -120,10 +122,13 @@ namespace gridfire::solver {
         );
 
         void cleanup_cvode_resources(bool memFree);
+
+        void log_step_diagnostics(const CVODEUserData& user_data) const;
     private:
         SUNContext m_sun_ctx = nullptr;
         void* m_cvode_mem = nullptr;
         N_Vector m_Y = nullptr;
+        N_Vector m_YErr = nullptr;
         SUNMatrix m_J = nullptr;
         SUNLinearSolver m_LS = nullptr;
 

src/include/gridfire/utils/logging.h

@@ -4,6 +4,12 @@
 
 #include <string>
 #include <vector>
+#include <sstream>
+#include <iomanip>
+#include <algorithm>
+#include <numeric>
+#include <iostream>
+
 
 namespace gridfire::utils {
     /**
@@ -62,4 +68,6 @@ namespace gridfire::utils {
         const double T9,
         const double rho
     );
+
+
 }

src/include/gridfire/utils/table_format.h

@@ -0,0 +1,121 @@
+#pragma once
+
+#include <iostream>
+#include <vector>
+#include <string>
+#include <sstream>
+#include <iomanip>
+#include <algorithm>
+#include <numeric>
+#include <memory>
+
+
+namespace gridfire::utils {
+    class ColumnBase {
+    public:
+        virtual ~ColumnBase() = default;
+        // Gets the string representation of the data at a given row
+        virtual std::string getCellData(size_t rowIndex) const = 0;
+        // Gets the header text for the column
+        virtual std::string getHeader() const = 0;
+        // Gets the number of data rows in the column
+        virtual size_t getRowCount() const = 0;
+    };
+
+    template<typename T>
+    class Column final : public ColumnBase {
+    public:
+        Column(const std::string& header, const std::vector<T>& data)
+            : m_header(header), m_data(data) {}
+
+        std::string getCellData(size_t rowIndex) const override {
+            std::stringstream ss;
+            if (rowIndex < m_data.size()) {
+                ss << m_data[rowIndex];
+            }
+            return ss.str();
+        }
+
+        std::string getHeader() const override {
+            return m_header;
+        }
+
+        size_t getRowCount() const override {
+            return m_data.size();
+        }
+    private:
+        std::string m_header;
+        std::vector<T> m_data;
+
+    };
+
+    inline std::string format_table(const std::string& tableName, const std::vector<std::unique_ptr<ColumnBase>>& columns) {
+        // --- 1. Handle Empty Table ---
+        if (columns.empty()) {
+            return tableName + "\n(Table has no columns)\n";
+        }
+
+        // --- 2. Determine dimensions and calculate column widths ---
+        size_t num_cols = columns.size();
+        size_t num_rows = 0;
+        for(const auto& col : columns) {
+            num_rows = std::max(num_rows, col->getRowCount());
+        }
+
+        std::vector<size_t> col_widths(num_cols);
+        for (size_t j = 0; j < num_cols; ++j) {
+            col_widths[j] = columns[j]->getHeader().length();
+            for (size_t i = 0; i < num_rows; ++i) {
+                col_widths[j] = std::max(col_widths[j], columns[j]->getCellData(i).length());
+            }
+        }
+
+        // --- 3. Build the table string using stringstream ---
+        std::stringstream table_ss;
+
+        // --- Table Title ---
+        size_t total_width = std::accumulate(col_widths.begin(), col_widths.end(), 0) + (num_cols * 3) + 1;
+        size_t title_padding = (total_width > tableName.length()) ? (total_width - tableName.length()) / 2 : 0;
+        table_ss << std::string(title_padding, ' ') << tableName << "\n";
+
+        // --- Helper to draw horizontal border ---
+        auto draw_border = [&]() {
+            table_ss << "+";
+            for (size_t width : col_widths) {
+                table_ss << std::string(width + 2, '-'); // +2 for padding
+                table_ss << "+";
+            }
+            table_ss << "\n";
+        };
+
+        // --- Draw Top Border ---
+        draw_border();
+
+        // --- Draw Header Row ---
+        table_ss << "|";
+        for (size_t j = 0; j < num_cols; ++j) {
+            table_ss << " " << std::left << std::setw(col_widths[j]) << columns[j]->getHeader() << " |";
+        }
+        table_ss << "\n";
+
+        // --- Draw Separator ---
+        draw_border();
+
+        // --- Draw Data Rows ---
+        for (size_t i = 0; i < num_rows; ++i) {
+            table_ss << "|";
+            for (size_t j = 0; j < num_cols; ++j) {
+                table_ss << " " << std::left << std::setw(col_widths[j]) << columns[j]->getCellData(i) << " |";
+            }
+            table_ss << "\n";
+        }
+
+        // --- Draw Bottom Border ---
+        draw_border();
+
+        return table_ss.str();
+    }
+
+
+
+}

src/lib/engine/diagnostics/dynamic_engine_diagnostics.cpp

@@ -0,0 +1,168 @@
+#include "gridfire/engine/diagnostics/dynamic_engine_diagnostics.h"
+#include "gridfire/engine/engine_abstract.h"
+#include "gridfire/utils/table_format.h"
+
+#include <vector>
+#include <string>
+#include <iostream>
+#include <iomanip>
+#include <algorithm>
+
+namespace gridfire::diagnostics {
+    void report_limiting_species(
+        const DynamicEngine& engine,
+        const std::vector<double>& Y_full,
+        const std::vector<double>& E_full,
+        const std::vector<double>& dydt_full,
+        const double relTol,
+        const double absTol,
+        const size_t top_n
+    ) {
+        struct SpeciesError {
+            std::string name;
+            double ratio;
+            double abundance;
+            double dydt;
+        };
+
+        const auto& species_list = engine.getNetworkSpecies();
+        std::vector<SpeciesError> errors;
+
+        for (size_t i = 0; i < species_list.size(); ++i) {
+            const double weight = relTol * std::abs(Y_full[i]) + absTol;
+            if (weight > 1e-99) { // Avoid division by zero for zero-abundance species
+                const double ratio = std::abs(E_full[i]) / weight;
+                errors.push_back({
+                    std::string(species_list[i].name()),
+                    ratio,
+                    Y_full[i],
+                    dydt_full[i]
+                });
+            }
+        }
+
+        // Sort by error ratio in descending order
+        std::ranges::sort(
+            errors,
+            [](const auto& a, const auto& b) {
+                return a.ratio > b.ratio;
+            }
+        );
+
+        std::vector<std::string> sorted_speciesNames;
+        std::vector<double> sorted_err_ratios;
+        std::vector<double> sorted_abundances;
+        std::vector<double> sorted_dydt;
+
+        for (size_t i = 0; i < std::min(top_n, errors.size()); ++i) {
+            sorted_speciesNames.push_back(errors[i].name);
+            sorted_err_ratios.push_back(errors[i].ratio);
+            sorted_abundances.push_back(errors[i].abundance);
+            sorted_dydt.push_back(errors[i].dydt);
+        }
+
+        std::vector<std::unique_ptr<utils::ColumnBase>> columns;
+        columns.push_back(std::make_unique<utils::Column<std::string>>("Species", sorted_speciesNames));
+        columns.push_back(std::make_unique<utils::Column<double>>("Error Ratio", sorted_err_ratios));
+        columns.push_back(std::make_unique<utils::Column<double>>("Abundance", sorted_abundances));
+        columns.push_back(std::make_unique<utils::Column<double>>("dY/dt", sorted_dydt));
+
+        std::cout << utils::format_table("Timestep Limiting Species", columns) << std::endl;
+    }
+
+    void inspect_species_balance(
+        const DynamicEngine& engine,
+        const std::string& species_name,
+        const std::vector<double>& Y_full,
+        const double T9,
+        const double rho
+    ) {
+        const auto& species_obj = fourdst::atomic::species.at(species_name);
+
+        std::vector<std::string> creation_ids, destruction_ids;
+        std::vector<int> creation_stoichiometry, destruction_stoichiometry;
+        std::vector<double> creation_flows, destruction_flows;
+        double total_creation_flow = 0.0;
+        double total_destruction_flow = 0.0;
+
+        for (const auto& reaction : engine.getNetworkReactions()) {
+            const int stoichiometry = reaction->stoichiometry(species_obj);
+            if (stoichiometry == 0) continue;
+
+            const double flow = engine.calculateMolarReactionFlow(*reaction, Y_full, T9, rho);
+
+            if (stoichiometry > 0) {
+                creation_ids.push_back(std::string(reaction->id()));
+                creation_stoichiometry.push_back(stoichiometry);
+                creation_flows.push_back(flow);
+                total_creation_flow += stoichiometry * flow;
+            } else {
+                destruction_ids.push_back(std::string(reaction->id()));
+                destruction_stoichiometry.push_back(stoichiometry);
+                destruction_flows.push_back(flow);
+                total_destruction_flow += std::abs(stoichiometry) * flow;
+            }
+        }
+
+        {
+            std::vector<std::unique_ptr<utils::ColumnBase>> columns;
+            columns.push_back(std::make_unique<utils::Column<std::string>>("Reaction ID", creation_ids));
+            columns.push_back(std::make_unique<utils::Column<int>>("Stoichiometry", creation_stoichiometry));
+            columns.push_back(std::make_unique<utils::Column<double>>("Molar Flow", creation_flows));
+            std::cout << utils::format_table("Creation Reactions for " + species_name, columns) << std::endl;
+        }
+
+        {
+            std::vector<std::unique_ptr<utils::ColumnBase>> columns;
+            columns.push_back(std::make_unique<utils::Column<std::string>>("Reaction ID", destruction_ids));
+            columns.push_back(std::make_unique<utils::Column<int>>("Stoichiometry", destruction_stoichiometry));
+            columns.push_back(std::make_unique<utils::Column<double>>("Molar Flow", destruction_flows));
+            std::cout << utils::format_table("Destruction Reactions for " + species_name, columns) << std::endl;
+        }
+
+        std::cout << "--- Balance Summary for " << species_name << " ---" << std::endl;
+        std::cout << "  Total Creation Rate:  " << std::scientific << total_creation_flow << " [mol/g/s]" << std::endl;
+        std::cout << "  Total Destruction Rate: " << std::scientific << total_destruction_flow << " [mol/g/s]" << std::endl;
+        std::cout << "  Net dY/dt:              " << std::scientific << (total_creation_flow - total_destruction_flow) << std::endl;
+        std::cout << "-----------------------------------" << std::endl;
+    }
+
+    void inspect_jacobian_stiffness(
+        const DynamicEngine& engine,
+        const std::vector<double>& Y_full,
+        const double T9,
+        const double rho
+    ) {
+        engine.generateJacobianMatrix(Y_full, T9, rho);
+        const auto& species_list = engine.getNetworkSpecies();
+
+        double max_diag = 0.0;
+        double max_off_diag = 0.0;
+        int max_diag_idx = -1;
+        int max_off_diag_i = -1, max_off_diag_j = -1;
+
+        for (size_t i = 0; i < species_list.size(); ++i) {
+            for (size_t j = 0; j < species_list.size(); ++j) {
+                const double val = std::abs(engine.getJacobianMatrixEntry(i, j));
+                if (i == j) {
+                    if (val > max_diag) { max_diag = val; max_diag_idx = i; }
+                } else {
+                    if (val > max_off_diag) { max_off_diag = val; max_off_diag_i = i; max_off_diag_j = j; }
+                }
+            }
+        }
+
+        std::cout << "\n--- Jacobian Stiffness Report ---" << std::endl;
+        if (max_diag_idx != -1) {
+            std::cout << "  Largest Diagonal Element (d(dYi/dt)/dYi): " << std::scientific << max_diag
+                      << " for species " << species_list[max_diag_idx].name() << std::endl;
+        }
+        if (max_off_diag_i != -1) {
+            std::cout << "  Largest Off-Diagonal Element (d(dYi/dt)/dYj): " << std::scientific << max_off_diag
+                      << " for d(" << species_list[max_off_diag_i].name()
+                      << ")/d(" << species_list[max_off_diag_j].name() << ")" << std::endl;
+        }
+        std::cout << "---------------------------------" << std::endl;
+    }
+
+}

src/lib/engine/engine_graph.cpp

@@ -78,6 +78,49 @@ namespace gridfire {
         }
     }
 
+    EnergyDerivatives GraphEngine::calculateEpsDerivatives(
+        const std::vector<double> &Y,
+        const double T9,
+        const double rho
+    ) const {
+        const size_t numSpecies = m_networkSpecies.size();
+        const size_t numADInputs = numSpecies + 2; // +2 for T9 and rho
+
+        if (Y.size() != numSpecies) {
+            LOG_ERROR(m_logger, "Input abundance vector size ({}) does not match number of species in the network ({}).",
+                      Y.size(), numSpecies);
+            throw std::invalid_argument("Input abundance vector size does not match number of species in the network.");
+        }
+
+        std::vector<double> x(numADInputs);
+        for (size_t i = 0; i < numSpecies; ++i) {
+            x[i] = Y[i];
+        }
+
+        x[numSpecies] = T9;
+        x[numSpecies + 1] = rho;
+
+        // Use reverse mode to get the gradient. W selects which dependent variable we care about, the Eps AD tape only has eps as a dependent variable so we just select set the 0th element to 1.
+        std::vector<double> w(1);
+        w[0] = 1.0; // We want the derivative of the energy generation rate
+
+        // Sweep the tape forward to record the function value at x
+        m_epsADFun.Forward(0, x);
+
+        // Extract the gradient at the previously evaluated point x using reverse mode
+        const std::vector<double> eps_derivatives = m_epsADFun.Reverse(1, w);
+
+        const double dEps_dT9 = eps_derivatives[numSpecies];
+        const double dEps_dRho = eps_derivatives[numSpecies + 1];
+
+        // Chain rule to scale from deps/dT9 to deps/dT
+        // dT9/dT = 1e-9
+        const double dEps_dT = dEps_dT9 * 1e-9;
+
+
+        return {dEps_dT, dEps_dRho};
+    }
+
     void GraphEngine::syncInternalMaps() {
         LOG_INFO(m_logger, "Synchronizing internal maps for REACLIB graph network (serif::network::GraphNetwork)...");
         collectNetworkSpecies();
@@ -86,7 +129,8 @@ namespace gridfire {
         collectAtomicReverseRateAtomicBases();
         generateStoichiometryMatrix();
         reserveJacobianMatrix();
-        recordADTape();
+        recordADTape(); // Record the AD tape for the RHS function
+        recordEpsADTape(); // Record the AD tape for the energy generation rate function
 
         const size_t n = m_rhsADFun.Domain();
         const size_t m = m_rhsADFun.Range();
@@ -161,13 +205,13 @@ namespace gridfire {
     // --- Basic Accessors and Queries ---
     const std::vector<fourdst::atomic::Species>& GraphEngine::getNetworkSpecies() const {
         // Returns a constant reference to the vector of unique species in the network.
-        LOG_TRACE_L3(m_logger, "Providing access to network species vector. Size: {}.", m_networkSpecies.size());
+        // LOG_TRACE_L3(m_logger, "Providing access to network species vector. Size: {}.", m_networkSpecies.size());
         return m_networkSpecies;
     }
 
     const reaction::ReactionSet& GraphEngine::getNetworkReactions() const {
         // Returns a constant reference to the set of reactions in the network.
-        LOG_TRACE_L3(m_logger, "Providing access to network reactions set. Size: {}.", m_reactions.size());
+        // LOG_TRACE_L3(m_logger, "Providing access to network reactions set. Size: {}.", m_reactions.size());
         return m_reactions;
     }
 
@@ -931,6 +975,43 @@ namespace gridfire {
                  adInput.size());
     }
 
+    void GraphEngine::recordEpsADTape() const {
+        LOG_TRACE_L1(m_logger, "Recording AD tape for the EPS calculation...");
+
+        const size_t numSpecies = m_networkSpecies.size();
+        if (numSpecies == 0) {
+            LOG_ERROR(m_logger, "Cannot record EPS tape: No species in the network.");
+            throw std::runtime_error("Cannot record EPS AD tape: No species in the network.");
+        }
+
+        const size_t numADInputs = numSpecies + 2; // Y + T9 + rho
+
+        std::vector<CppAD::AD<double>> adInput(numADInputs, 0.0);
+
+        for (size_t i = 0; i < numSpecies; ++i) {
+            adInput[i] = static_cast<CppAD::AD<double>>(1.0 / static_cast<double>(numSpecies)); // Uniform distribution
+        }
+        adInput[numSpecies]     = 1.0; // Dummy T9
+        adInput[numSpecies + 1] = 1.0; // Dummy rho
+
+        CppAD::Independent(adInput);
+
+        const std::vector<CppAD::AD<double>> adY(adInput.begin(), adInput.begin() + numSpecies);
+        const CppAD::AD<double> adT9  = adInput[numSpecies];
+        const CppAD::AD<double> adRho = adInput[numSpecies + 1];
+
+        StepDerivatives<CppAD::AD<double>> derivatives = calculateAllDerivatives<CppAD::AD<double>>(adY, adT9, adRho);
+
+        std::vector<CppAD::AD<double>> adOutput(1);
+        adOutput[0] = derivatives.nuclearEnergyGenerationRate;
+        m_epsADFun.Dependent(adInput, adOutput);
+        // m_epsADFun.optimize();
+
+        LOG_TRACE_L1(m_logger, "AD tape recorded successfully for the EPS calculation. Number of independent variables: {}.", adInput.size());
+
+
+    }
+
     void GraphEngine::collectAtomicReverseRateAtomicBases() {
         m_atomicReverseRates.clear();
         m_atomicReverseRates.reserve(m_reactions.size());

src/lib/engine/procedures/priming.cpp

@@ -35,10 +35,42 @@ namespace gridfire {
         return dominateReaction;
     }
 
-
+    /**
+     * @brief Primes absent species in the network to their equilibrium abundances using a robust, two-stage approach.
+     *
+     * @details This function implements a robust network priming algorithm that avoids the pitfalls of
+     * sequential, one-by-one priming. The previous, brittle method could allow an early priming
+     * reaction to consume all of a shared reactant, starving later reactions. This new, two-stage
+     * method ensures that all priming reactions are considered collectively, competing for the
+     * same limited pool of initial reactants in a physically consistent manner.
+     *
+     * The algorithm proceeds in three main stages:
+     * 1.  **Calculation Stage:** It first loops through all species that need priming. For each one,
+     * it calculates its theoretical equilibrium mass fraction and identifies the dominant
+     * creation channel. Crucially, it *does not* modify any abundances at this stage. Instead,
+     * it stores these calculations as a list of "mass transfer requests".
+     *
+     * 2.  **Collective Scaling Stage:** It then processes the full list of requests to determine the
+     * total "debit" required from each reactant. By comparing these total debits to the
+     * initially available mass of each reactant, it calculates a single, global `scalingFactor`.
+     * If any reactant is overdrawn, this factor will be less than 1.0, ensuring that no
+     * reactant's abundance can go negative.
+     *
+     * 3.  **Application Stage:** Finally, it loops through the requests again, applying the mass
+     * transfers. Each calculated equilibrium mass fraction and corresponding reactant debit is
+     * multiplied by the global `scalingFactor` before being applied to the final composition.
+     * This ensures that if resources are limited, all primed species are scaled down proportionally.
+     *
+     * @param netIn Input network data containing initial composition, temperature, and density.
+     * @param engine DynamicEngine used to build and evaluate the reaction network.
+     * @return PrimingReport encapsulating the results of the priming operation, including the new
+     * robustly primed composition.
+     */
     PrimingReport primeNetwork(const NetIn& netIn, DynamicEngine& engine) {
         auto logger = LogManager::getInstance().getLogger("log");
 
+        // --- Initial Setup ---
+        // Identify all species with zero initial mass fraction that need to be primed.
         std::vector<Species> speciesToPrime;
         for (const auto &entry: netIn.composition | std::views::values) {
             if (entry.mass_fraction() == 0.0) {
@@ -47,6 +79,7 @@ namespace gridfire {
         }
         LOG_DEBUG(logger, "Priming {} species in the network.", speciesToPrime.size());
 
+        // If no species need priming, return immediately.
         PrimingReport report;
         if (speciesToPrime.empty()) {
             report.primedComposition = netIn.composition;
@@ -58,43 +91,46 @@ namespace gridfire {
         const double T9 = netIn.temperature / 1e9;
         const double rho = netIn.density;
         const auto initialReactionSet = engine.getNetworkReactions();
-
         report.status = PrimingReportStatus::FULL_SUCCESS;
         report.success = true;
 
-        // --- 1: pack composition into internal map ---
+        // Create a mutable map of the mass fractions that we will modify.
         std::unordered_map<Species, double> currentMassFractions;
         for (const auto& entry : netIn.composition | std::views::values) {
             currentMassFractions[entry.isotope()] = entry.mass_fraction();
         }
+        // Ensure all species to be primed exist in the map, initialized to zero.
         for (const auto& entry : speciesToPrime) {
-            currentMassFractions[entry] = 0.0; // Initialize priming species with 0 mass fraction
+            currentMassFractions[entry] = 0.0;
         }
 
-        std::unordered_map<Species, double> totalMassFractionChanges;
+        // Rebuild the engine with the full network to ensure all possible creation channels are available.
-
         engine.rebuild(netIn.composition, NetworkBuildDepth::Full);
 
-        for (const auto& primingSpecies : speciesToPrime) {
+        // --- STAGE 1: Calculation and Bookkeeping (No Modifications) ---
-            LOG_TRACE_L3(logger, "Priming species: {}", primingSpecies.name());
+        // In this stage, we calculate all required mass transfers but do not apply them yet.
 
-            // Create a temporary composition from the current internal state for the primer
+        // A struct to hold the result of each individual priming calculation.
+        struct MassTransferRequest {
+            Species species_to_prime;
+            double equilibrium_mass_fraction;
+            std::vector<Species> reactants;
+        };
+        std::vector<MassTransferRequest> requests;
+
+        for (const auto& primingSpecies : speciesToPrime) {
+            // Create a temporary composition reflecting the current state for rate calculations.
             Composition tempComp;
             for(const auto& [sp, mf] : currentMassFractions) {
                 tempComp.registerSymbol(std::string(sp.name()));
-                if (mf < 0.0 && std::abs(mf) < 1e-16) {
+                tempComp.setMassFraction(sp, std::max(0.0, mf));
-                    tempComp.setMassFraction(sp, 0.0); // Avoid negative mass fractions
-                } else {
-                    tempComp.setMassFraction(sp, mf);
-                }
             }
-            tempComp.finalize(true); // Finalize with normalization
+            tempComp.finalize(true);
 
             NetIn tempNetIn = netIn;
             tempNetIn.composition = tempComp;
 
             NetworkPrimingEngineView primer(primingSpecies, engine);
-
             if (primer.getNetworkReactions().size() == 0) {
                 LOG_ERROR(logger, "No priming reactions found for species {}.", primingSpecies.name());
                 report.success = false;
@@ -106,60 +142,87 @@ namespace gridfire {
             const double destructionRateConstant = calculateDestructionRateConstant(primer, primingSpecies, Y, T9, rho);
 
             if (destructionRateConstant > 1e-99) {
-                double equilibriumMassFraction = 0.0;
                 const double creationRate = calculateCreationRate(primer, primingSpecies, Y, T9, rho);
-                equilibriumMassFraction = (creationRate / destructionRateConstant) * primingSpecies.mass();
+                double equilibriumMassFraction = (creationRate / destructionRateConstant) * primingSpecies.mass();
-                if (std::isnan(equilibriumMassFraction)) {
+                if (std::isnan(equilibriumMassFraction)) equilibriumMassFraction = 0.0;
-                    LOG_WARNING(logger, "Equilibrium mass fraction for {} is NaN. Setting to 0.0. This is likely not an issue. It probably originates from all reactions leading to creation and destruction being frozen out. In that case 0.0 should be a good approximation. Hint: This happens often when the network temperature is very the low. ", primingSpecies.name());
-                    equilibriumMassFraction = 0.0;
-                }
-                LOG_TRACE_L3(logger, "Found equilibrium for {}: X_eq = {:.4e}", primingSpecies.name(), equilibriumMassFraction);
 
                 if (const reaction::Reaction* dominantChannel = findDominantCreationChannel(primer, primingSpecies, Y, T9, rho)) {
-                    LOG_TRACE_L3(logger, "Dominant creation channel for {}: {}", primingSpecies.name(), dominantChannel->id());
+                    // Store the request instead of applying it immediately.
-
+                    requests.push_back({primingSpecies, equilibriumMassFraction, dominantChannel->reactants()});
-                    double totalReactantMass = 0.0;
-                    for (const auto& reactant : dominantChannel->reactants()) {
-                        totalReactantMass += reactant.mass();
-                    }
-
-                    double scalingFactor = 1.0;
-                    for (const auto& reactant : dominantChannel->reactants()) {
-                        const double massToSubtract = equilibriumMassFraction * (reactant.mass() / totalReactantMass);
-                        double availableMass = 0.0;
-                        if (currentMassFractions.contains(reactant)) {
-                            availableMass = currentMassFractions.at(reactant);
-                        }
-                        if (massToSubtract > availableMass && availableMass > 0) {
-                            scalingFactor = std::min(scalingFactor, availableMass / massToSubtract);
-                        }
-                    }
-
-                    if (scalingFactor < 1.0) {
-                        LOG_WARNING(logger, "Priming for {} was limited by reactant availability. Scaling transfer by {:.4e}", primingSpecies.name(), scalingFactor);
-                        equilibriumMassFraction *= scalingFactor;
-                    }
-
-                    // Update the internal mass fraction map and accumulate total changes
-                    totalMassFractionChanges[primingSpecies] += equilibriumMassFraction;
-                    currentMassFractions[primingSpecies] += equilibriumMassFraction;
-
-                    for (const auto& reactant : dominantChannel->reactants()) {
-                        const double massToSubtract = equilibriumMassFraction * (reactant.mass() / totalReactantMass);
-                        totalMassFractionChanges[reactant] -= massToSubtract;
-                        currentMassFractions[reactant] -= massToSubtract;
-                    }
                 } else {
-                     LOG_ERROR(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name());
+                    LOG_ERROR(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name());
-                     report.status = PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL;
+                    report.status = PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL;
-                     totalMassFractionChanges[primingSpecies] += 1e-40;
-                     currentMassFractions[primingSpecies] += 1e-40;
                 }
             } else {
                 LOG_WARNING(logger, "No destruction channel found for {}. Using fallback abundance.", primingSpecies.name());
-                totalMassFractionChanges[primingSpecies] += 1e-40;
+                // For species with no destruction, we can't calculate an equilibrium.
-                currentMassFractions[primingSpecies] += 1e-40;
+                // We add a request with a tiny fallback abundance to ensure it exists in the network.
-                report.status = PrimingReportStatus::BASE_NETWORK_TOO_SHALLOW;
+                requests.push_back({primingSpecies, 1e-40, {}});
+            }
+        }
+
+        // --- STAGE 2: Collective Scaling Based on Reactant Availability ---
+        // Now, we determine the total demand for each reactant and find a global scaling factor.
+
+        std::unordered_map<Species, double> total_mass_debits;
+        for (const auto& req : requests) {
+            if (req.reactants.empty()) continue; // Skip fallbacks which don't consume reactants.
+
+            double totalReactantMass = 0.0;
+            for (const auto& reactant : req.reactants) {
+                totalReactantMass += reactant.mass();
+            }
+            if (totalReactantMass == 0.0) continue;
+
+            for (const auto& reactant : req.reactants) {
+                const double massToSubtract = req.equilibrium_mass_fraction * (reactant.mass() / totalReactantMass);
+                total_mass_debits[reactant] += massToSubtract;
+            }
+        }
+
+        double globalScalingFactor = 1.0;
+        for (const auto& [reactant, total_debit] : total_mass_debits) {
+            double availableMass;
+            if (currentMassFractions.contains(reactant)) {
+                availableMass = currentMassFractions.at(reactant);
+            } else {
+                availableMass = 0.0;
+            }
+            if (total_debit > availableMass && availableMass > 0) {
+                globalScalingFactor = std::min(globalScalingFactor, availableMass / total_debit);
+            }
+        }
+
+        if (globalScalingFactor < 1.0) {
+            LOG_WARNING(logger, "Priming was limited by reactant availability. All transfers will be scaled by {:.4e}", globalScalingFactor);
+        }
+
+        // --- STAGE 3: Application of Scaled Mass Transfers ---
+        // Finally, apply all the transfers, scaled by our global factor.
+
+        std::unordered_map<Species, double> totalMassFractionChanges;
+        for (const auto&[species_to_prime, equilibrium_mass_fraction, reactants] : requests) {
+            const double scaled_equilibrium_mf = equilibrium_mass_fraction * globalScalingFactor;
+
+            // Add the scaled mass to the primed species.
+            currentMassFractions.at(species_to_prime) += scaled_equilibrium_mf;
+            totalMassFractionChanges[species_to_prime] += scaled_equilibrium_mf;
+
+            // Subtract the scaled mass from the reactants.
+            if (!reactants.empty()) {
+                double totalReactantMass = 0.0;
+                for (const auto& reactant : reactants) {
+                    totalReactantMass += reactant.mass();
+                }
+                if (totalReactantMass == 0.0) continue;
+
+                for (const auto& reactant : reactants) {
+                    const double massToSubtract = scaled_equilibrium_mf * (reactant.mass() / totalReactantMass);
+                    if (massToSubtract != 0) {
+                        currentMassFractions.at(reactant) -= massToSubtract;
+                        totalMassFractionChanges[reactant] -= massToSubtract;
+                    }
+                }
             }
         }
 
@@ -168,14 +231,9 @@ namespace gridfire {
         std::vector<double> final_mass_fractions;
         for(const auto& [species, mass_fraction] : currentMassFractions) {
             final_symbols.emplace_back(species.name());
-            if (mass_fraction < 0.0 && std::abs(mass_fraction) < 1e-16) {
+            final_mass_fractions.push_back(std::max(0.0, mass_fraction)); // Ensure no negative mass fractions.
-                final_mass_fractions.push_back(0.0); // Avoid negative mass fractions
-            } else {
-                final_mass_fractions.push_back(mass_fraction);
-            }
         }
 
-        // Create the final composition object from the pre-normalized mass fractions
         Composition primedComposition(final_symbols, final_mass_fractions, true);
 
         report.primedComposition = primedComposition;
@@ -183,10 +241,166 @@ namespace gridfire {
             report.massFractionChanges.emplace_back(species, change);
         }
 
+        // Restore the engine to its original, smaller network state.
         engine.setNetworkReactions(initialReactionSet);
         return report;
     }
 
+
+
+    // PrimingReport primeNetwork(const NetIn& netIn, DynamicEngine& engine) {
+    //     auto logger = LogManager::getInstance().getLogger("log");
+    //
+    //     std::vector<Species> speciesToPrime;
+    //     for (const auto &entry: netIn.composition | std::views::values) {
+    //         std::cout << "Checking species: " << entry.isotope().name() << " with mass fraction: " << entry.mass_fraction() << std::endl;
+    //         if (entry.mass_fraction() == 0.0) {
+    //             speciesToPrime.push_back(entry.isotope());
+    //         }
+    //     }
+    //     LOG_DEBUG(logger, "Priming {} species in the network.", speciesToPrime.size());
+    //
+    //     PrimingReport report;
+    //     if (speciesToPrime.empty()) {
+    //         report.primedComposition = netIn.composition;
+    //         report.success = true;
+    //         report.status = PrimingReportStatus::NO_SPECIES_TO_PRIME;
+    //         return report;
+    //     }
+    //
+    //     const double T9 = netIn.temperature / 1e9;
+    //     const double rho = netIn.density;
+    //     const auto initialReactionSet = engine.getNetworkReactions();
+    //
+    //     report.status = PrimingReportStatus::FULL_SUCCESS;
+    //     report.success = true;
+    //
+    //     // --- 1: pack composition into internal map ---
+    //     std::unordered_map<Species, double> currentMassFractions;
+    //     for (const auto& entry : netIn.composition | std::views::values) {
+    //         currentMassFractions[entry.isotope()] = entry.mass_fraction();
+    //     }
+    //     for (const auto& entry : speciesToPrime) {
+    //         currentMassFractions[entry] = 0.0; // Initialize priming species with 0 mass fraction
+    //     }
+    //
+    //     std::unordered_map<Species, double> totalMassFractionChanges;
+    //
+    //     engine.rebuild(netIn.composition, NetworkBuildDepth::Full);
+    //
+    //     for (const auto& primingSpecies : speciesToPrime) {
+    //         LOG_TRACE_L3(logger, "Priming species: {}", primingSpecies.name());
+    //
+    //         // Create a temporary composition from the current internal state for the primer
+    //         Composition tempComp;
+    //         for(const auto& [sp, mf] : currentMassFractions) {
+    //             tempComp.registerSymbol(std::string(sp.name()));
+    //             if (mf < 0.0 && std::abs(mf) < 1e-16) {
+    //                 tempComp.setMassFraction(sp, 0.0); // Avoid negative mass fractions
+    //             } else {
+    //                 tempComp.setMassFraction(sp, mf);
+    //             }
+    //         }
+    //         tempComp.finalize(true); // Finalize with normalization
+    //
+    //         NetIn tempNetIn = netIn;
+    //         tempNetIn.composition = tempComp;
+    //
+    //         NetworkPrimingEngineView primer(primingSpecies, engine);
+    //
+    //         if (primer.getNetworkReactions().size() == 0) {
+    //             LOG_ERROR(logger, "No priming reactions found for species {}.", primingSpecies.name());
+    //             report.success = false;
+    //             report.status = PrimingReportStatus::FAILED_TO_FIND_PRIMING_REACTIONS;
+    //             continue;
+    //         }
+    //
+    //         const auto Y = primer.mapNetInToMolarAbundanceVector(tempNetIn);
+    //         const double destructionRateConstant = calculateDestructionRateConstant(primer, primingSpecies, Y, T9, rho);
+    //
+    //         if (destructionRateConstant > 1e-99) {
+    //             double equilibriumMassFraction = 0.0;
+    //             const double creationRate = calculateCreationRate(primer, primingSpecies, Y, T9, rho);
+    //             equilibriumMassFraction = (creationRate / destructionRateConstant) * primingSpecies.mass();
+    //             if (std::isnan(equilibriumMassFraction)) {
+    //                 LOG_WARNING(logger, "Equilibrium mass fraction for {} is NaN. Setting to 0.0. This is likely not an issue. It probably originates from all reactions leading to creation and destruction being frozen out. In that case 0.0 should be a good approximation. Hint: This happens often when the network temperature is very the low. ", primingSpecies.name());
+    //                 equilibriumMassFraction = 0.0;
+    //             }
+    //             LOG_TRACE_L3(logger, "Found equilibrium for {}: X_eq = {:.4e}", primingSpecies.name(), equilibriumMassFraction);
+    //
+    //             if (const reaction::Reaction* dominantChannel = findDominantCreationChannel(primer, primingSpecies, Y, T9, rho)) {
+    //                 LOG_TRACE_L3(logger, "Dominant creation channel for {}: {}", primingSpecies.name(), dominantChannel->id());
+    //
+    //                 double totalReactantMass = 0.0;
+    //                 for (const auto& reactant : dominantChannel->reactants()) {
+    //                     totalReactantMass += reactant.mass();
+    //                 }
+    //
+    //                 double scalingFactor = 1.0;
+    //                 for (const auto& reactant : dominantChannel->reactants()) {
+    //                     const double massToSubtract = equilibriumMassFraction * (reactant.mass() / totalReactantMass);
+    //                     double availableMass = 0.0;
+    //                     if (currentMassFractions.contains(reactant)) {
+    //                         availableMass = currentMassFractions.at(reactant);
+    //                     }
+    //                     if (massToSubtract > availableMass && availableMass > 0) {
+    //                         scalingFactor = std::min(scalingFactor, availableMass / massToSubtract);
+    //                     }
+    //                 }
+    //
+    //                 if (scalingFactor < 1.0) {
+    //                     LOG_WARNING(logger, "Priming for {} was limited by reactant availability. Scaling transfer by {:.4e}", primingSpecies.name(), scalingFactor);
+    //                     equilibriumMassFraction *= scalingFactor;
+    //                 }
+    //
+    //                 // Update the internal mass fraction map and accumulate total changes
+    //                 totalMassFractionChanges[primingSpecies] += equilibriumMassFraction;
+    //                 currentMassFractions.at(primingSpecies) += equilibriumMassFraction;
+    //
+    //                 for (const auto& reactant : dominantChannel->reactants()) {
+    //                     const double massToSubtract = equilibriumMassFraction * (reactant.mass() / totalReactantMass);
+    //                     std::cout << "[Priming: " << primingSpecies.name() << ", Channel: " << dominantChannel->id() <<  "] Subtracting " << massToSubtract << " from reactant " << reactant.name() << std::endl;
+    //                     totalMassFractionChanges[reactant] -= massToSubtract;
+    //                     currentMassFractions[reactant] -= massToSubtract;
+    //                 }
+    //             } else {
+    //                  LOG_ERROR(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name());
+    //                  report.status = PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL;
+    //                  totalMassFractionChanges[primingSpecies] += 1e-40;
+    //                  currentMassFractions.at(primingSpecies) += 1e-40;
+    //             }
+    //         } else {
+    //             LOG_WARNING(logger, "No destruction channel found for {}. Using fallback abundance.", primingSpecies.name());
+    //             totalMassFractionChanges.at(primingSpecies) += 1e-40;
+    //             currentMassFractions.at(primingSpecies) += 1e-40;
+    //             report.status = PrimingReportStatus::BASE_NETWORK_TOO_SHALLOW;
+    //         }
+    //     }
+    //
+    //     // --- Final Composition Construction ---
+    //     std::vector<std::string> final_symbols;
+    //     std::vector<double> final_mass_fractions;
+    //     for(const auto& [species, mass_fraction] : currentMassFractions) {
+    //         final_symbols.emplace_back(species.name());
+    //         if (mass_fraction < 0.0 && std::abs(mass_fraction) < 1e-16) {
+    //             final_mass_fractions.push_back(0.0); // Avoid negative mass fractions
+    //         } else {
+    //             final_mass_fractions.push_back(mass_fraction);
+    //         }
+    //     }
+    //
+    //     // Create the final composition object from the pre-normalized mass fractions
+    //     Composition primedComposition(final_symbols, final_mass_fractions, true);
+    //
+    //     report.primedComposition = primedComposition;
+    //     for (const auto& [species, change] : totalMassFractionChanges) {
+    //         report.massFractionChanges.emplace_back(species, change);
+    //     }
+    //
+    //     engine.setNetworkReactions(initialReactionSet);
+    //     return report;
+    // }
+
   double calculateDestructionRateConstant(
         const DynamicEngine& engine,
         const fourdst::atomic::Species& species,

src/lib/engine/views/engine_adaptive.cpp

@@ -167,6 +167,17 @@ namespace gridfire {
         return culledResults;
     }
 
+    EnergyDerivatives AdaptiveEngineView::calculateEpsDerivatives(
+        const std::vector<double> &Y_culled,
+        const double T9,
+        const double rho
+    ) const {
+        validateState();
+        const auto Y_full = mapCulledToFull(Y_culled);
+
+        return m_baseEngine.calculateEpsDerivatives(Y_full, T9, rho);
+    }
+
     void AdaptiveEngineView::generateJacobianMatrix(
         const std::vector<double> &Y_dynamic,
         const double T9,

src/lib/engine/views/engine_defined.cpp

@@ -48,6 +48,17 @@ namespace gridfire {
         return definedResults;
     }
 
+    EnergyDerivatives DefinedEngineView::calculateEpsDerivatives(
+        const std::vector<double> &Y_dynamic,
+        const double T9,
+        const double rho
+    ) const {
+        validateNetworkState();
+
+        const auto Y_full = mapViewToFull(Y_dynamic);
+        return m_baseEngine.calculateEpsDerivatives(Y_full, T9, rho);
+    }
+
     void DefinedEngineView::generateJacobianMatrix(
         const std::vector<double> &Y_dynamic,
         const double T9,

src/lib/engine/views/engine_multiscale.cpp

@@ -1,5 +1,6 @@
 #include "gridfire/engine/views/engine_multiscale.h"
 #include "gridfire/exceptions/error_engine.h"
+#include "gridfire/engine/procedures/priming.h"
 
 #include <stdexcept>
 #include <vector>
@@ -11,6 +12,7 @@
 #include <ranges>
 #include <algorithm>
 
+#include "gridfire/utils/logging.h"
 #include "quill/LogMacros.h"
 #include "quill/Logger.h"
 
@@ -200,6 +202,14 @@ namespace gridfire {
         return deriv;
     }
 
+    EnergyDerivatives MultiscalePartitioningEngineView::calculateEpsDerivatives(
+        const std::vector<double> &Y,
+        const double T9,
+        const double rho
+    ) const {
+        return m_baseEngine.calculateEpsDerivatives(Y, T9, rho);
+    }
+
     void MultiscalePartitioningEngineView::generateJacobianMatrix(
         const std::vector<double> &Y_full,
         const double T9,
@@ -228,12 +238,12 @@ namespace gridfire {
         const int i_full,
         const int j_full
     ) const {
-        // // Check if the species we are differentiating with respect to is algebraic or dynamic. If it is algebraic we can reduce the work significantly...
+        // Check if the species we are differentiating with respect to is algebraic or dynamic. If it is algebraic we can reduce the work significantly...
-        // if (std::ranges::contains(m_algebraic_species_indices, j_full)) {
+        if (std::ranges::contains(m_algebraic_species_indices, j_full)) {
-        //     const auto& species = m_baseEngine.getNetworkSpecies()[j_full];
+            // const auto& species = m_baseEngine.getNetworkSpecies()[j_full];
-        //     // If j is algebraic, we can return 0.0 since the Jacobian entry for algebraic species is always zero.
+            // If j is algebraic, we can return 0.0 since the Jacobian entry for algebraic species is always zero.
-        //     return 0.0;
+            return 0.0;
-        // }
+        }
         // Otherwise we need to query the full jacobian
         return m_baseEngine.getJacobianMatrixEntry(i_full, j_full);
     }
@@ -481,9 +491,42 @@ namespace gridfire {
         LOG_TRACE_L1(m_logger, "Identifying potential seed species for candidate pools...");
         const std::vector<QSEGroup> candidate_groups = constructCandidateGroups(connected_pools, Y, T9, rho);
         LOG_TRACE_L1(m_logger, "Found {} candidate QSE groups for further analysis", candidate_groups.size());
+        LOG_TRACE_L2(
+            m_logger,
+            "{}",
+             [&]() -> std::string {
+                 std::stringstream ss;
+                 int j = 0;
+                 for (const auto& group : candidate_groups) {
+                     ss << "CandidateQSEGroup(Algebraic: {";
+                     int i = 0;
+                     for (const auto& index : group.algebraic_indices) {
+                         ss << m_baseEngine.getNetworkSpecies()[index].name();
+                         if (i < group.algebraic_indices.size() - 1) {
+                             ss << ", ";
+                         }
+                     }
+                     ss << "}, Seed: {";
+                     i = 0;
+                     for (const auto& index : group.seed_indices) {
+                         ss << m_baseEngine.getNetworkSpecies()[index].name();
+                         if (i < group.seed_indices.size() - 1) {
+                             ss << ", ";
+                         }
+                         i++;
+                     }
+                     ss << "})";
+                     if (j < candidate_groups.size() - 1) {
+                         ss << ", ";
+                     }
+                     j++;
+                 }
+                 return ss.str();
+             }()
+        );
 
         LOG_TRACE_L1(m_logger, "Validating candidate groups with flux analysis...");
-        const std::vector<QSEGroup> validated_groups = validateGroupsWithFluxAnalysis(candidate_groups, Y, T9, rho);
+        const auto [validated_groups, invalidate_groups] = validateGroupsWithFluxAnalysis(candidate_groups, Y, T9, rho);
         LOG_TRACE_L1(
             m_logger,
             "Validated {} group(s) QSE groups. {}",
@@ -507,6 +550,13 @@ namespace gridfire {
             }()
         );
 
+        // Push the invalidated groups' species into the dynamic set
+        for (const auto& group : invalidate_groups) {
+            for (const auto& index : group.algebraic_indices) {
+                m_dynamic_species.push_back(m_baseEngine.getNetworkSpecies()[index]);
+            }
+        }
+
         m_qse_groups = validated_groups;
         LOG_TRACE_L1(m_logger, "Identified {} QSE groups.", m_qse_groups.size());
 
@@ -542,6 +592,10 @@ namespace gridfire {
             m_qse_groups.size() == 1 ? "" : "s"
         );
 
+        // throw std::runtime_error(
+        //     "Partitioning complete. Throwing an error to end the program during debugging. This error should not be caught by the caller. "
+        // );
+
     }
 
     void MultiscalePartitioningEngineView::partitionNetwork(
@@ -846,6 +900,8 @@ namespace gridfire {
         return m_baseEngine.getSpeciesIndex(species);
     }
 
+
+
     std::vector<std::vector<size_t>> MultiscalePartitioningEngineView::partitionByTimescale(
         const std::vector<double>& Y_full,
         const double T9,
@@ -853,19 +909,31 @@ namespace gridfire {
     ) const {
         LOG_TRACE_L1(m_logger, "Partitioning by timescale...");
         const auto result= m_baseEngine.getSpeciesDestructionTimescales(Y_full, T9, rho);
+        const auto netTimescale = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho);
+        std::cout << gridfire::utils::formatNuclearTimescaleLogString(m_baseEngine, Y_full, T9, rho) << std::endl;
         if (!result) {
-            LOG_ERROR(m_logger, "Failed to get species timescales due to stale engine state");
+            LOG_ERROR(m_logger, "Failed to get species destruction timescales due to stale engine state");
             m_logger->flush_log();
-            throw exceptions::StaleEngineError("Failed to get species timescales due to stale engine state");
+            throw exceptions::StaleEngineError("Failed to get species destruction timescales due to stale engine state");
+        }
+        if (!netTimescale) {
+            LOG_ERROR(m_logger, "Failed to get net species timescales due to stale engine state");
+            m_logger->flush_log();
+            throw exceptions::StaleEngineError("Failed to get net species timescales due to stale engine state");
         }
         const std::unordered_map<Species, double>& all_timescales = result.value();
+        const std::unordered_map<Species, double>& net_timescales = netTimescale.value();
         const auto& all_species = m_baseEngine.getNetworkSpecies();
 
         std::vector<std::pair<double, size_t>> sorted_timescales;
         for (size_t i = 0; i < all_species.size(); ++i) {
             double timescale = all_timescales.at(all_species[i]);
+            double net_timescale = net_timescales.at(all_species[i]);
             if (std::isfinite(timescale) && timescale > 0) {
+                LOG_TRACE_L3(m_logger, "Species {} has finite destruction timescale: destruction: {} s, net: {} s", all_species[i].name(), timescale, net_timescale);
                 sorted_timescales.emplace_back(timescale, i);
+            } else {
+                LOG_TRACE_L3(m_logger, "Species {} has infinite or negative destruction timescale: destruction: {} s, net: {} s", all_species[i].name(), timescale, net_timescale);
             }
         }
 
@@ -971,98 +1039,120 @@ namespace gridfire {
 
     }
 
-    std::vector<MultiscalePartitioningEngineView::QSEGroup>
+    std::pair<std::vector<MultiscalePartitioningEngineView::QSEGroup>, std::vector<MultiscalePartitioningEngineView::
+    QSEGroup>>
     MultiscalePartitioningEngineView::validateGroupsWithFluxAnalysis(
         const std::vector<QSEGroup> &candidate_groups,
         const std::vector<double> &Y,
         const double T9, const double rho
     ) const {
-        constexpr double FLUX_RATIO_THRESHOLD = 100;
+        constexpr double FLUX_RATIO_THRESHOLD = 5;
-        std::vector<QSEGroup> validated_groups = candidate_groups;
+        std::vector<QSEGroup> validated_groups;
-        for (auto& group : validated_groups) {
+        std::vector<QSEGroup> invalidated_groups;
-            double internal_flux = 0.0;
+        validated_groups.reserve(candidate_groups.size());
-            double external_flux = 0.0;
+        for (auto& group : candidate_groups) {
-
+            const std::unordered_set<size_t> algebraic_group_members(
-            const std::unordered_set<size_t> group_members(
+                group.algebraic_indices.begin(),
-                group.species_indices.begin(),
+                group.algebraic_indices.end()
-                group.species_indices.end()
             );
 
+            const std::unordered_set<size_t> seed_group_members(
+                group.seed_indices.begin(),
+                group.seed_indices.end()
+            );
+
+            double coupling_flux = 0.0;
+            double leakage_flux = 0.0;
+
             for (const auto& reaction: m_baseEngine.getNetworkReactions()) {
                 const double flow = std::abs(m_baseEngine.calculateMolarReactionFlow(*reaction, Y, T9, rho));
                 if (flow == 0.0) {
                     continue; // Skip reactions with zero flow
                 }
-                bool has_internal_reactant = false;
+                bool has_internal_algebraic_reactant = false;
-                bool has_external_reactant = false;
 
                 for (const auto& reactant : reaction->reactants()) {
-                    if (group_members.contains(m_baseEngine.getSpeciesIndex(reactant))) {
+                    if (algebraic_group_members.contains(m_baseEngine.getSpeciesIndex(reactant))) {
-                        has_internal_reactant = true;
+                        has_internal_algebraic_reactant = true;
-                    } else {
-                        has_external_reactant = true;
                     }
                 }
 
-                bool has_internal_product = false;
+                bool has_internal_algebraic_product = false;
-                bool has_external_product = false;
 
                 for (const auto& product : reaction->products()) {
-                    if (group_members.contains(m_baseEngine.getSpeciesIndex(product))) {
+                    if (algebraic_group_members.contains(m_baseEngine.getSpeciesIndex(product))) {
-                        has_internal_product = true;
+                        has_internal_algebraic_product = true;
-                    } else {
-                        has_external_product = true;
                     }
                 }
 
-                // Classify the reaction based on its participants
+                if (!has_internal_algebraic_product && !has_internal_algebraic_reactant) {
-                if ((has_internal_reactant && has_internal_product) && !(has_external_reactant || has_external_product)) {
+                    LOG_TRACE_L3(m_logger, "{}: Skipping reaction {} as it has no internal algebraic species in reactants or products.", group.toString(m_baseEngine), reaction->id());
-                    LOG_TRACE_L3(
+                    continue;
-                        m_logger,
-                        "Reaction {} is internal to the group containing {} and contributes to internal flux by {}",
-                        reaction->id(),
-                        [&]() -> std::string {
-                            std::stringstream ss;
-                            int count = 0;
-                            for (const auto& idx : group.algebraic_indices) {
-                                ss << m_baseEngine.getNetworkSpecies()[idx].name();
-                                if (count < group.species_indices.size() - 1) {
-                                    ss << ", ";
-                                }
-                                count++;
-                            }
-                            return ss.str();
-                        }(),
-                        flow
-                    );
-                    internal_flux += flow; // Internal flux within the group
-                } else if ((has_internal_reactant || has_internal_product) && (has_external_reactant || has_external_product)) {
-                    LOG_TRACE_L3(
-                        m_logger,
-                        "Reaction {} is external to the group containing {} and contributes to external flux by {}",
-                        reaction->id(),
-                        [&]() -> std::string {
-                            std::stringstream ss;
-                            int count = 0;
-                            for (const auto& idx : group.algebraic_indices) {
-                                ss << m_baseEngine.getNetworkSpecies()[idx].name();
-                                if (count < group.species_indices.size() - 1) {
-                                    ss << ", ";
-                                }
-                                count++;
-                            }
-                            return ss.str();
-                        }(),
-                        flow
-                    );
-                    external_flux += flow; // External flux to/from the group
                 }
-                // otherwise the reaction is fully decoupled from the QSE group and can be ignored.
+
+                double algebraic_participants = 0;
+                double seed_participants = 0;
+                double external_participants = 0;
+
+                std::unordered_set<Species> participants;
+                for(const auto& p : reaction->reactants()) participants.insert(p);
+                for(const auto& p : reaction->products()) participants.insert(p);
+
+                for (const auto& species : participants) {
+                    const size_t species_idx = m_baseEngine.getSpeciesIndex(species);
+                    if (algebraic_group_members.contains(species_idx)) {
+                        LOG_TRACE_L3(m_logger, "{}: Species {} is an algebraic participant in reaction {}.", group.toString(m_baseEngine), species.name(), reaction->id());
+                        algebraic_participants++;
+                    } else if (seed_group_members.contains(species_idx)) {
+                        LOG_TRACE_L3(m_logger, "{}: Species {} is a seed participant in reaction {}.", group.toString(m_baseEngine), species.name(), reaction->id());
+                        seed_participants++;
+                    } else {
+                        LOG_TRACE_L3(m_logger, "{}: Species {} is an external participant in reaction {}.", group.toString(m_baseEngine), species.name(), reaction->id());
+                        external_participants++;
+                    }
+                }
+
+                const double total_participants = algebraic_participants + seed_participants + external_participants;
+
+                if (total_participants == 0) {
+                    LOG_CRITICAL(m_logger, "Some catastrophic error has occurred. Reaction {} has no participants.", reaction->id());
+                    throw std::runtime_error("Some catastrophic error has occurred. Reaction " + std::string(reaction->id()) + " has no participants.");
+                }
+
+                const double leakage_fraction = external_participants / total_participants;
+                const double coupling_fraction = (algebraic_participants + seed_participants) / total_participants;
+
+                leakage_flux += flow * leakage_fraction;
+                coupling_flux += flow * coupling_fraction;
             }
-            if (internal_flux > FLUX_RATIO_THRESHOLD * external_flux) {
+
+            // if (leakage_flux < 1e-99) {
+            //     LOG_TRACE_L1(
+            //         m_logger,
+            //         "Group containing {} is in equilibrium due to vanishing leakage: leakage flux = {}, coupling flux = {}, ratio = {}",
+            //         [&]() -> std::string {
+            //             std::stringstream ss;
+            //             int count = 0;
+            //             for (const auto& idx : group.algebraic_indices) {
+            //                 ss << m_baseEngine.getNetworkSpecies()[idx].name();
+            //                 if (count < group.species_indices.size() - 1) {
+            //                     ss << ", ";
+            //                 }
+            //                 count++;
+            //             }
+            //             return ss.str();
+            //         }(),
+            //         leakage_flux,
+            //         coupling_flux,
+            //         coupling_flux / leakage_flux
+            //     );
+            //     validated_groups.emplace_back(group);
+            //     validated_groups.back().is_in_equilibrium = true;
+            // } else if ((coupling_flux / leakage_flux ) > FLUX_RATIO_THRESHOLD) {
+            if ((coupling_flux / leakage_flux ) > FLUX_RATIO_THRESHOLD) {
                 LOG_TRACE_L1(
                     m_logger,
-                    "Group containing {} is in equilibrium: internal flux = {}, external flux = {}, ratio = {}",
+                    "Group containing {} is in equilibrium due to high coupling flux threshold: leakage flux = {}, coupling flux = {}, ratio = {} (Threshold: {})",
                     [&]() -> std::string {
                         std::stringstream ss;
                         int count = 0;
@@ -1075,15 +1165,17 @@ namespace gridfire {
                         }
                         return ss.str();
                     }(),
-                    internal_flux,
+                    leakage_flux,
-                    external_flux,
+                    coupling_flux,
-                    internal_flux / external_flux
+                    coupling_flux / leakage_flux,
+                    FLUX_RATIO_THRESHOLD
                 );
-                group.is_in_equilibrium = true; // This group is in equilibrium if internal flux is significantly larger than external flux.
+                validated_groups.emplace_back(group);
+                validated_groups.back().is_in_equilibrium = true;
             } else {
                 LOG_TRACE_L1(
                     m_logger,
-                    "Group containing {} is NOT in equilibrium: internal flux = {}, external flux = {}, ratio = {}",
+                    "Group containing {} is NOT in equilibrium: leakage flux = {}, coupling flux = {}, ratio = {} (Threshold: {})",
                     [&]() -> std::string {
                         std::stringstream ss;
                         int count = 0;
@@ -1096,14 +1188,16 @@ namespace gridfire {
                         }
                         return ss.str();
                     }(),
-                    internal_flux,
+                    leakage_flux,
-                    external_flux,
+                    coupling_flux,
-                    internal_flux / external_flux
+                    coupling_flux / leakage_flux,
+                    FLUX_RATIO_THRESHOLD
                 );
-                group.is_in_equilibrium = false;
+                invalidated_groups.emplace_back(group);
+                invalidated_groups.back().is_in_equilibrium = false;
             }
         }
-        return validated_groups;
+        return {validated_groups, invalidated_groups};
     }
 
     std::vector<double> MultiscalePartitioningEngineView::solveQSEAbundances(
@@ -1543,10 +1637,11 @@ namespace gridfire {
         hash_combine(seed, bin(m_T9, m_cacheConfig.T9_tol));
         hash_combine(seed, bin(m_rho, m_cacheConfig.rho_tol));
 
+        double negThresh = 1e-10; // Threshold for considering a value as negative.
         for (double Yi : m_Y) {
-            if (Yi < 0.0 && std::abs(Yi) < 1e-20) {
+            if (Yi < 0.0 && std::abs(Yi) < negThresh) {
                 Yi = 0.0; // Avoid negative abundances
-            } else if (Yi < 0.0 && std::abs(Yi) >= 1e-20) {
+            } else if (Yi < 0.0 && std::abs(Yi) >= negThresh) {
                 throw std::invalid_argument("Yi should be positive for valid hashing (expected Yi > 0, received " + std::to_string(Yi) + ")");
             }
             hash_combine(seed, bin(Yi, m_cacheConfig.Yi_tol));
@@ -1580,6 +1675,55 @@ namespace gridfire {
         return !(*this == other);
     }
 
+    std::string MultiscalePartitioningEngineView::QSEGroup::toString() const {
+        std::stringstream ss;
+        ss << "QSEGroup(Algebraic: [";
+        size_t count = 0;
+        for (const auto& idx : algebraic_indices) {
+            ss << idx;
+            if (count < algebraic_indices.size() - 1) {
+                ss << ", ";
+            }
+            count++;
+        }
+        ss << "], Seed: [";
+        count = 0;
+        for (const auto& idx : seed_indices) {
+            ss << idx;
+            if (count < seed_indices.size() - 1) {
+                ss << ", ";
+            }
+            count++;
+        }
+        ss << "], Mean Timescale: " << mean_timescale << ", Is In Equilibrium: " << (is_in_equilibrium ? "True" : "False") << ")";
+        return ss.str();
+
+    }
+
+    std::string MultiscalePartitioningEngineView::QSEGroup::toString(DynamicEngine &engine) const {
+        std::stringstream ss;
+        ss << "QSEGroup(Algebraic: [";
+        size_t count = 0;
+        for (const auto& idx : algebraic_indices) {
+            ss << engine.getNetworkSpecies()[idx].name();
+            if (count < algebraic_indices.size() - 1) {
+                ss << ", ";
+            }
+            count++;
+        }
+        ss << "], Seed: [";
+        count = 0;
+        for (const auto& idx : seed_indices) {
+            ss << engine.getNetworkSpecies()[idx].name();
+            if (count < seed_indices.size() - 1) {
+                ss << ", ";
+            }
+            count++;
+        }
+        ss << "], Mean Timescale: " << mean_timescale << ", Is In Equilibrium: " << (is_in_equilibrium ? "True" : "False") << ")";
+        return ss.str();
+    }
+
     void MultiscalePartitioningEngineView::CacheStats::hit(const operators op) {
         if (op == operators::All) {
             throw std::invalid_argument("Cannot use 'ALL' as an operator for a hit");

src/lib/solver/solver.cpp

@@ -137,7 +137,7 @@ namespace gridfire::solver {
         std::vector<std::string> speciesNames;
         speciesNames.reserve(numSpecies);
         for (const auto& species : m_engine.getNetworkSpecies()) {
-            speciesNames.push_back(std::string(species.name()));
+            speciesNames.emplace_back(species.name());
         }
 
         Composition outputComposition(speciesNames);
@@ -145,10 +145,19 @@ namespace gridfire::solver {
         outputComposition.finalize(true);
 
         NetOut netOut;
-        netOut.composition = std::move(outputComposition);
+        netOut.composition = outputComposition;
         netOut.energy = accumulated_energy; // Specific energy rate
         netOut.num_steps = manager.m_num_steps;
 
+        auto [dEps_dT, dEps_dRho] = m_engine.calculateEpsDerivatives(
+            std::vector<double>(Y.begin(), Y.begin() + numSpecies), // TODO: This narrowing should probably be solved. Its possible unforeseen bugs will arise from this
+            T9,
+            netIn.density
+        );
+
+        netOut.dEps_dT = dEps_dT;
+        netOut.dEps_dRho = dEps_dRho;
+
         return netOut;
     }
 

src/lib/solver/strategies/CVODE_solver_strategy.cpp

@@ -1,6 +1,10 @@
 #include "gridfire/solver/strategies/CVODE_solver_strategy.h"
 
 #include "gridfire/network.h"
+#include "gridfire/utils/table_format.h"
+#include "gridfire/engine/diagnostics/dynamic_engine_diagnostics.h"
+
+#include "quill/LogMacros.h"
 
 #include "fourdst/composition/composition.h"
 
@@ -10,6 +14,9 @@
 #include <string>
 #include <unordered_map>
 #include <stdexcept>
+#include <algorithm>
+
+#include "fourdst/composition/exceptions/exceptions_composition.h"
 
 
 namespace {
@@ -75,11 +82,12 @@ namespace gridfire::solver {
         // TODO: In order to support MPI this function must be changed
         const int flag = SUNContext_Create(SUN_COMM_NULL, &m_sun_ctx);
         if (flag < 0) {
-            throw std::runtime_error("Failed to create SUNDIALS context (Errno: " + std::to_string(flag) + ")");
+            throw std::runtime_error("Failed to create SUNDIALS context (SUNDIALS Errno: " + std::to_string(flag) + ")");
         }
     }
 
     CVODESolverStrategy::~CVODESolverStrategy() {
+        std::cout << "Cleaning up CVODE resources..." << std::endl;
         cleanup_cvode_resources(true);
 
         if (m_sun_ctx) {
@@ -111,7 +119,7 @@ namespace gridfire::solver {
         [[maybe_unused]] double last_callback_time = 0;
         m_num_steps = 0;
         double accumulated_energy = 0.0;
-
+        int total_update_stages_triggered = 0;
         while (current_time < netIn.tMax) {
             try {
                 user_data.T9 = T9;
@@ -138,22 +146,47 @@ namespace gridfire::solver {
                 double last_step_size;
                 CVodeGetNumSteps(m_cvode_mem, &n_steps);
                 CVodeGetLastStep(m_cvode_mem, &last_step_size);
-                std::cout << std::scientific << std::setprecision(3) << "Step: " << std::setw(6) << n_steps << " | Time: " << current_time << " | Last Step Size: " << last_step_size << std::endl;
+                long int nliters, nlcfails;
+                CVodeGetNumNonlinSolvIters(m_cvode_mem, &nliters);
+                CVodeGetNumNonlinSolvConvFails(m_cvode_mem, &nlcfails);
+
+                sunrealtype* y_data = N_VGetArrayPointer(m_Y);
+                const double current_energy = y_data[numSpecies]; // Specific energy rate
+                std::cout << std::scientific << std::setprecision(3)
+                    << "Step: " << std::setw(6) << n_steps
+                    << " | Time: " << current_time << " [s]"
+                    << " | Last Step Size: " << last_step_size
+                    << " | Accumulated Energy: " << current_energy << " [erg/g]"
+                    << " | NonlinIters: " << std::setw(2) << nliters
+                    << " | ConvFails: " << std::setw(2) << nlcfails
+                    << std::endl;
+
+                // if (n_steps % 50 == 0) {
+                //     std::cout << "Logging step diagnostics at step " << n_steps << "..." << std::endl;
+                //     log_step_diagnostics(user_data);
+                // }
+                // if (n_steps == 300) {
+                //     log_step_diagnostics(user_data);
+                //     exit(0);
+                // }
+                // log_step_diagnostics(user_data);
 
             } catch (const exceptions::StaleEngineTrigger& e) {
                 exceptions::StaleEngineTrigger::state staleState = e.getState();
                 accumulated_energy += e.energy(); // Add the specific energy rate to the accumulated energy
-                // total_update_stages_triggered++;
+                LOG_INFO(
+                    m_logger,
+                    "Engine Update Triggered due to StaleEngineTrigger exception at time {} ({} update{} triggered). Current total specific energy {} [erg/g]",
+                    current_time,
+                    total_update_stages_triggered,
+                    total_update_stages_triggered == 1 ? "" : "s",
+                    accumulated_energy);
+                total_update_stages_triggered++;
 
                 fourdst::composition::Composition temp_comp;
                 std::vector<double> mass_fractions;
                 size_t num_species_at_stop = e.numSpecies();
 
-                if (num_species_at_stop != m_engine.getNetworkSpecies().size()) {
-                    throw std::runtime_error(
-                        "StaleEngineError state has a different number of species than the engine. This should not happen."
-                    );
-                }
                 mass_fractions.reserve(num_species_at_stop);
 
                 for (size_t i = 0; i < num_species_at_stop; ++i) {
@@ -170,13 +203,22 @@ namespace gridfire::solver {
                 netInTemp.composition = temp_comp;
 
                 fourdst::composition::Composition currentComposition = m_engine.update(netInTemp);
+                LOG_INFO(
+                    m_logger,
+                    "Due to a triggered stale engine the composition was updated from size {} to {} species.",
+                    num_species_at_stop,
+                    m_engine.getNetworkSpecies().size()
+                );
 
                 numSpecies = m_engine.getNetworkSpecies().size();
                 N = numSpecies + 1;
 
-                initialize_cvode_integration_resources(N, numSpecies, current_time, temp_comp, absTol, relTol, accumulated_energy);
+                initialize_cvode_integration_resources(N, numSpecies, current_time, currentComposition, absTol, relTol, accumulated_energy);
 
                 check_cvode_flag(CVodeReInit(m_cvode_mem, current_time, m_Y), "CVodeReInit");
+            } catch (fourdst::composition::exceptions::InvalidCompositionError& e) {
+                log_step_diagnostics(user_data);
+                std::rethrow_exception(std::make_exception_ptr(e));
             }
         }
 
@@ -195,7 +237,7 @@ namespace gridfire::solver {
         std::vector<std::string> speciesNames;
         speciesNames.reserve(numSpecies);
         for (const auto& species : m_engine.getNetworkSpecies()) {
-            speciesNames.push_back(std::string(species.name()));
+            speciesNames.emplace_back(species.name());
         }
 
         fourdst::composition::Composition outputComposition(speciesNames);
@@ -203,9 +245,21 @@ namespace gridfire::solver {
         outputComposition.finalize(true);
 
         NetOut netOut;
-        netOut.composition = std::move(outputComposition);
+        netOut.composition = outputComposition;
         netOut.energy = accumulated_energy;
         check_cvode_flag(CVodeGetNumSteps(m_cvode_mem, reinterpret_cast<long int *>(&netOut.num_steps)), "CVodeGetNumSteps");
+
+        outputComposition.setCompositionMode(false); // set to number fraction mode
+        std::vector<double> Y = outputComposition.getNumberFractionVector(); // TODO need to ensure that the canonical vector representation is used throughout the code to make sure tracking does not get messed up
+        auto [dEps_dT, dEps_dRho] = m_engine.calculateEpsDerivatives(
+            std::vector<double>(Y.begin(), Y.begin() + numSpecies), // TODO: This narrowing should probably be solved. Its possible unforeseen bugs will arise from this
+            T9,
+            netIn.density
+        );
+
+        netOut.dEps_dT = dEps_dT;
+        netOut.dEps_dRho = dEps_dRho;
+
         return netOut;
     }
 
@@ -232,10 +286,11 @@ namespace gridfire::solver {
         void *user_data
     ) {
         auto* data = static_cast<CVODEUserData*>(user_data);
-        auto* instance = data->solver_instance;
+        const auto* instance = data->solver_instance;
 
         try {
-            return instance->calculate_rhs(t, y, ydot, data);
+            instance->calculate_rhs(t, y, ydot, data);
+            return 0;
         } catch (const exceptions::StaleEngineTrigger& e) {
             data->captured_exception = std::make_unique<exceptions::StaleEngineTrigger>(e);
             return 1; // 1 Indicates a recoverable error, CVODE will retry the step
@@ -279,7 +334,7 @@ namespace gridfire::solver {
         return 0;
     }
 
-    int CVODESolverStrategy::calculate_rhs(
+    void CVODESolverStrategy::calculate_rhs(
         const sunrealtype t,
         const N_Vector y,
         const N_Vector ydot,
@@ -304,7 +359,6 @@ namespace gridfire::solver {
             ydot_data[i] = dydt[i];
         }
         ydot_data[numSpecies] = nuclearEnergyGenerationRate; // Set the last element to the specific energy rate
-        return 0;
     }
 
     void CVODESolverStrategy::initialize_cvode_integration_resources(
@@ -319,6 +373,7 @@ namespace gridfire::solver {
         cleanup_cvode_resources(false); // Cleanup any existing resources before initializing new ones
 
         m_Y = init_sun_vector(N, m_sun_ctx);
+        m_YErr = N_VClone(m_Y);
 
         sunrealtype *y_data = N_VGetArrayPointer(m_Y);
         for (size_t i = 0; i < numSpecies; i++) {
@@ -335,6 +390,8 @@ namespace gridfire::solver {
         check_cvode_flag(CVodeInit(m_cvode_mem, cvode_rhs_wrapper, current_time, m_Y), "CVodeInit");
         check_cvode_flag(CVodeSStolerances(m_cvode_mem, relTol, absTol), "CVodeSStolerances");
 
+        check_cvode_flag(CVodeSetMaxStep(m_cvode_mem, 1.0e20), "CVodeSetMaxStep");
+
         m_J = SUNDenseMatrix(static_cast<sunindextype>(N), static_cast<sunindextype>(N), m_sun_ctx);
         check_cvode_flag(m_J == nullptr ? -1 : 0, "SUNDenseMatrix");
         m_LS = SUNLinSol_Dense(m_Y, m_J, m_sun_ctx);
@@ -348,10 +405,12 @@ namespace gridfire::solver {
         if (m_LS) SUNLinSolFree(m_LS);
         if (m_J) SUNMatDestroy(m_J);
         if (m_Y) N_VDestroy(m_Y);
+        if (m_YErr) N_VDestroy(m_YErr);
 
         m_LS = nullptr;
         m_J = nullptr;
         m_Y = nullptr;
+        m_YErr = nullptr;
 
         if (memFree) {
             if (m_cvode_mem) CVodeFree(&m_cvode_mem);
@@ -359,4 +418,79 @@ namespace gridfire::solver {
         }
     }
 
+    void CVODESolverStrategy::log_step_diagnostics(const CVODEUserData &user_data) const {
+        check_cvode_flag(CVodeGetEstLocalErrors(m_cvode_mem, m_YErr), "CVodeGetEstLocalErrors");
+
+        sunrealtype *y_data = N_VGetArrayPointer(m_Y);
+        sunrealtype *y_err_data = N_VGetArrayPointer(m_YErr);
+
+
+        std::vector<double> err_ratios;
+        std::vector<std::string> speciesNames;
+
+        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);
+
+        const size_t num_components = N_VGetLength(m_Y);
+        err_ratios.resize(num_components - 1);
+
+        std::vector<double> Y_full(y_data, y_data + num_components - 1);
+
+        std::ranges::replace_if(
+            Y_full,
+            [](const double val) {
+                return val < 0.0 && val > -1.0e-16;
+            },
+            0.0
+        );
+
+        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
+
+            const double err_ratio = std::abs(y_err_data[i]) / weight;
+
+            err_ratios[i] = err_ratio;
+            speciesNames.push_back(std::string(user_data.networkSpecies->at(i).name()));
+        }
+
+        if (err_ratios.empty()) {
+            return;
+        }
+
+        std::vector<size_t> indices(speciesNames.size());
+        for (size_t i = 0; i < indices.size(); ++i) {
+            indices[i] = i;
+        }
+
+        std::ranges::sort(
+            indices,
+            [&err_ratios](const size_t i1, const size_t i2) {
+                return err_ratios[i1] > err_ratios[i2];
+            }
+        );
+
+        std::vector<std::string> sorted_speciesNames;
+        std::vector<double> sorted_err_ratios;
+
+        sorted_speciesNames.reserve(indices.size());
+        sorted_err_ratios.reserve(indices.size());
+
+        for (const auto idx: indices) {
+            sorted_speciesNames.push_back(speciesNames[idx]);
+            sorted_err_ratios.push_back(err_ratios[idx]);
+        }
+
+
+
+        std::vector<std::unique_ptr<utils::ColumnBase>> columns;
+        columns.push_back(std::make_unique<utils::Column<std::string>>("Species", sorted_speciesNames));
+        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, Y_full, user_data.T9, user_data.rho);
+        diagnostics::inspect_species_balance(*user_data.engine, "N-14", Y_full, user_data.T9, user_data.rho);
+        diagnostics::inspect_species_balance(*user_data.engine, "n-1", Y_full, user_data.T9, user_data.rho);
+
+    }
 }

src/meson.build

@@ -9,6 +9,7 @@ gridfire_sources = files(
     'lib/engine/views/engine_priming.cpp',
     'lib/engine/procedures/priming.cpp',
     'lib/engine/procedures/construction.cpp',
+    'lib/engine/diagnostics/dynamic_engine_diagnostics.cpp',
     'lib/reaction/reaction.cpp',
     'lib/reaction/reaclib.cpp',
     'lib/io/network_file.cpp',
@@ -53,4 +54,9 @@ gridfire_dep = declare_dependency(
 
 install_subdir('include/gridfire', install_dir: get_option('includedir'))
 
-subdir('python')
+if get_option('build-python')
+    message('Configuring Python bindings...')
+    subdir('src-pybind')
+else
+    message('Skipping Python bindings...')
+endif

subprojects/fourdst.wrap

@@ -1,4 +1,4 @@
 [wrap-git]
 url = https://github.com/4D-STAR/fourdst
-revision = v0.7.0
+revision = v0.8.0
 depth = 1

tests/graphnet_sandbox/main.cpp

@@ -99,7 +99,7 @@ int main(int argc, char* argv[]){
 
     g_previousHandler = std::set_terminate(quill_terminate_handler);
     quill::Logger* logger = fourdst::logging::LogManager::getInstance().getLogger("log");
-    logger->set_log_level(quill::LogLevel::TraceL1);
+    logger->set_log_level(quill::LogLevel::TraceL3);
     LOG_INFO(logger, "Starting Adaptive Engine View Example...");
 
     using namespace gridfire;
@@ -122,12 +122,15 @@ int main(int argc, char* argv[]){
     netIn.temperature = 1.5e7;
     netIn.density = 1.6e2;
     netIn.energy = 0;
-    netIn.tMax = 5e17;
+    netIn.tMax = 1e2;
     // netIn.tMax = 1e-14;
     netIn.dt0 = 1e-12;
 
     GraphEngine ReaclibEngine(composition, partitionFunction, NetworkBuildDepth::SecondOrder);
+
     ReaclibEngine.setUseReverseReactions(false);
+    ReaclibEngine.setPrecomputation(false);
+
     MultiscalePartitioningEngineView partitioningView(ReaclibEngine);
     AdaptiveEngineView adaptiveView(partitioningView);
 
@@ -144,5 +147,6 @@ int main(int argc, char* argv[]){
     double finalHydrogen = netOut.composition.getMassFraction("H-1");
     double fractionalConsumedHydrogen = (initialHydrogen - finalHydrogen) / initialHydrogen * 100.0;
     std::cout << "Fractional consumed hydrogen: " << fractionalConsumedHydrogen << "%" << std::endl;
+    std::cout << netOut << std::endl;
 
 }