feat(GridFire): stabalized network, increased performance, evolving over 10Gyr timescales now with ~correct results

src/network/include/gridfire/engine/engine_abstract.h

@@ -106,7 +106,7 @@ namespace gridfire {
          * time derivatives of all species and the specific nuclear energy generation
          * rate for the current state.
          */
-        [[nodiscard]] virtual StepDerivatives<double> calculateRHSAndEnergy(
+        [[nodiscard]] virtual std::expected<StepDerivatives<double>, expectations::StaleEngineError> calculateRHSAndEnergy(
             const std::vector<double>& Y,
             double T9,
             double rho
@@ -217,6 +217,8 @@ namespace gridfire {
          */
         [[nodiscard]] virtual const reaction::LogicalReactionSet& getNetworkReactions() const = 0;
 
+        virtual void setNetworkReactions(const reaction::LogicalReactionSet& reactions) = 0;
+
         /**
          * @brief Compute timescales for all species in the network.
          *
@@ -228,7 +230,13 @@ namespace gridfire {
          * This method estimates the timescale for abundance change of each species,
          * which can be used for timestep control, diagnostics, and reaction network culling.
          */
-        [[nodiscard]] virtual std::unordered_map<fourdst::atomic::Species, double> getSpeciesTimescales(
+        [[nodiscard]] virtual std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesTimescales(
+            const std::vector<double>& Y,
+            double T9,
+            double rho
+        ) const = 0;
+
+        [[nodiscard]] virtual std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesDestructionTimescales(
             const std::vector<double>& Y,
             double T9,
             double rho

src/network/include/gridfire/engine/engine_graph.h

@@ -144,7 +144,7 @@ namespace gridfire {
          *
          * @see StepDerivatives
          */
-        [[nodiscard]] StepDerivatives<double> calculateRHSAndEnergy(
+        [[nodiscard]] std::expected<StepDerivatives<double>, expectations::StaleEngineError> calculateRHSAndEnergy(
             const std::vector<double>& Y,
             const double T9,
             const double rho
@@ -215,6 +215,8 @@ namespace gridfire {
          */
         [[nodiscard]] const reaction::LogicalReactionSet& getNetworkReactions() const override;
 
+        void setNetworkReactions(const reaction::LogicalReactionSet& reactions) override;
+
         /**
          * @brief Gets an entry from the previously generated Jacobian matrix.
          *
@@ -268,7 +270,13 @@ namespace gridfire {
          * This method estimates the timescale for abundance change of each species,
          * which can be used for timestep control or diagnostics.
          */
-        [[nodiscard]] std::unordered_map<fourdst::atomic::Species, double> getSpeciesTimescales(
+        [[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesTimescales(
+            const std::vector<double>& Y,
+            double T9,
+            double rho
+        ) const override;
+
+        [[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesDestructionTimescales(
             const std::vector<double>& Y,
             double T9,
             double rho
@@ -511,7 +519,7 @@ namespace gridfire {
          * to store the partial derivatives of the right-hand side of the ODE
          * with respect to the species abundances.
          */
-        void reserveJacobianMatrix();
+        void reserveJacobianMatrix() const;
 
         /**
          * @brief Records the AD tape for the right-hand side of the ODE.
@@ -795,43 +803,12 @@ namespace gridfire {
             );
 
             const T molarReactionFlow = forwardMolarReactionFlow - reverseMolarFlow; // Net molar reaction flow
-            // std::stringstream ss;
-            // ss << "Forward: " << forwardMolarReactionFlow
-            //    << ", Reverse: " << reverseMolarFlow
-            //    << ", Net: " << molarReactionFlow;
-            // LOG_TRACE_L3(
-            //     m_logger,
-            //     "Reaction: {}, {}",
-            //     reaction.peName(),
-            //     ss.str()
-            // );
-            // std::cout << "Forward molar flow for reaction " << reaction.peName() << ": " << forwardMolarReactionFlow << std::endl;
-            // std::cout << "Reverse molar flow for reaction " << reaction.peName() << ": " << reverseMolarFlow << std::endl;
-            // std::cout << "Net molar flow for reaction " << reaction.peName() << ": " << molarReactionFlow << std::endl;
 
             // 3. Use the rate to update all relevant species derivatives (dY/dt)
             for (size_t speciesIndex = 0; speciesIndex < m_networkSpecies.size(); ++speciesIndex) {
-                const auto& species = m_networkSpecies[speciesIndex];
                 const T nu_ij = static_cast<T>(m_stoichiometryMatrix(speciesIndex, reactionIndex));
-                // bool taping = false;
-                // if constexpr (std::is_same_v<T, CppAD::AD<double>>) {
-                //     taping = true;
-                // }
-                // if (species.name() == "F-17" && nu_ij != static_cast<T>(0.0)) {
-                //     T f17dydt = result.dydt[speciesIndex] + threshold_flag * nu_ij * molarReactionFlow;
-                //     std::string tstring = taping ? "During Taping" : "During Evaluation";
-                //     std::cout << "(" << tstring << ") F-17 Debugging. For Reaction " << reaction.id() << " (" << reaction.peName() << "): "
-                //               << "nu_ij: " << nu_ij << ", molarReactionFlow: " << molarReactionFlow
-                //               << ", threshold_flag: " << threshold_flag << ", dY/dt: " << f17dydt << ", Y: " << Y[speciesIndex]
-                //               << std::endl;
-                //     isF17 = true;
-                // }
                 result.dydt[speciesIndex] += threshold_flag * nu_ij * molarReactionFlow;
             }
-            // if (isF17) {
-            //     std::cout << "=========================================================" << std::endl;
-            //     isF17 = false;
-            // }
         }
 
         T massProductionRate = static_cast<T>(0.0); // [mol][s^-1]

src/network/include/gridfire/engine/types/reporting.h

@@ -33,7 +33,7 @@ namespace gridfire {
 
         friend std::ostream& operator<<(std::ostream& os, const PrimingReport& report) {
             std::stringstream ss;
-            std::string successStr = report.success ? "true" : "false";
+            const std::string successStr = report.success ? "true" : "false";
             ss << "PrimingReport(success=" << successStr
                << ", status=" << PrimingReportStatusStrings[report.status] << ")";
             return os << ss.str();

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

@@ -102,7 +102,7 @@ namespace gridfire {
          * @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called).
          * @see AdaptiveEngineView::update()
          */
-        [[nodiscard]] StepDerivatives<double> calculateRHSAndEnergy(
+        [[nodiscard]] std::expected<StepDerivatives<double>, expectations::StaleEngineError> calculateRHSAndEnergy(
             const std::vector<double> &Y_culled,
             const double T9,
             const double rho
@@ -205,6 +205,8 @@ namespace gridfire {
          */
         [[nodiscard]] const reaction::LogicalReactionSet& getNetworkReactions() const override;
 
+        void setNetworkReactions(const reaction::LogicalReactionSet& reactions) override;
+
         /**
          * @brief Computes timescales for all active species in the network.
          *
@@ -218,12 +220,18 @@ namespace gridfire {
          *
          * @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called).
          */
-        [[nodiscard]] std::unordered_map<fourdst::atomic::Species, double> getSpeciesTimescales(
+        [[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesTimescales(
             const std::vector<double> &Y_culled,
             double T9,
             double rho
         ) const override;
 
+        [[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesDestructionTimescales(
+            const std::vector<double> &Y,
+            double T9,
+            double rho
+        ) const override;
+
         /**
          * @brief Gets the base engine.
          * @return A const reference to the base engine.
@@ -443,7 +451,6 @@ namespace gridfire {
 
         typedef std::pair<std::unordered_set<const reaction::LogicalReaction*>, std::unordered_set<fourdst::atomic::Species>> RescueSet;
         [[nodiscard]] RescueSet rescueEdgeSpeciesDestructionChannel(
-            const std::vector<ReactionFlow>& allFlows,
             const std::vector<double>& Y_full,
             const double T9,
             const double rho,

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

@@ -37,7 +37,7 @@ namespace gridfire{
          *
          * @throws std::runtime_error If the view is stale (i.e., `update()` has not been called after `setNetworkFile()`).
          */
-        StepDerivatives<double> calculateRHSAndEnergy(
+        std::expected<StepDerivatives<double>, expectations::StaleEngineError> calculateRHSAndEnergy(
             const std::vector<double>& Y_defined,
             const double T9,
             const double rho
@@ -115,6 +115,8 @@ namespace gridfire{
          * @throws std::runtime_error If the view is stale.
          */
         const reaction::LogicalReactionSet& getNetworkReactions() const override;
+
+        void setNetworkReactions(const reaction::LogicalReactionSet& reactions) override;
         /**
          * @brief Computes timescales for all active species in the network.
          *
@@ -125,7 +127,13 @@ namespace gridfire{
          *
          * @throws std::runtime_error If the view is stale.
          */
-        std::unordered_map<fourdst::atomic::Species, double> getSpeciesTimescales(
+        [[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesTimescales(
+            const std::vector<double>& Y_defined,
+            const double T9,
+            const double rho
+        ) const override;
+
+        [[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesDestructionTimescales(
             const std::vector<double>& Y_defined,
             const double T9,
             const double rho
@@ -243,6 +251,9 @@ namespace gridfire{
         size_t mapViewToFullReactionIndex(size_t definedReactionIndex) const;
 
         void validateNetworkState() const;
+
+        void collect(const std::vector<std::string>& peNames);
+
     };
 
     class FileDefinedEngineView final: public DefinedEngineView {

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

@@ -66,7 +66,7 @@ namespace gridfire {
 
         [[nodiscard]] const std::vector<fourdst::atomic::Species> & getNetworkSpecies() const override;
 
-        [[nodiscard]] StepDerivatives<double> calculateRHSAndEnergy(
+        [[nodiscard]] std::expected<StepDerivatives<double>, expectations::StaleEngineError> calculateRHSAndEnergy(
             const std::vector<double> &Y_full,
             double T9,
             double rho
@@ -99,7 +99,17 @@ namespace gridfire {
 
         [[nodiscard]] const reaction::LogicalReactionSet & getNetworkReactions() const override;
 
-        [[nodiscard]] std::unordered_map<fourdst::atomic::Species, double> getSpeciesTimescales(
+        void setNetworkReactions(
+            const reaction::LogicalReactionSet &reactions
+        ) override;
+
+        [[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesTimescales(
+            const std::vector<double> &Y,
+            double T9,
+            double rho
+        ) const override;
+
+        [[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesDestructionTimescales(
             const std::vector<double> &Y,
             double T9,
             double rho
@@ -119,6 +129,13 @@ namespace gridfire {
 
         const DynamicEngine & getBaseEngine() const override;
 
+        std::vector<std::vector<size_t>> analyzeTimescalePoolConnectivity(
+            const std::vector<std::vector<size_t>> &timescale_pools,
+            const std::vector<double> &Y,
+            double T9,
+            double rho
+        ) const;
+
         void partitionNetwork(
             const std::vector<double>& Y,
             double T9,
@@ -162,6 +179,12 @@ namespace gridfire {
             bool is_in_equilibrium = false;      ///< Flag set by flux analysis.
             std::set<size_t> algebraic_indices; ///< Indices of algebraic species in this group.
             std::set<size_t> seed_indices; ///< Indices of dynamic species in this group.
+            double mean_timescale; ///< Mean timescale of the group.
+
+            bool operator<(const QSEGroup& other) const;
+            bool operator>(const QSEGroup& other) const;
+            bool operator==(const QSEGroup& other) const;
+            bool operator!=(const QSEGroup& other) const;
 
             friend std::ostream& operator<<(std::ostream& os, const QSEGroup& group) {
                 os << "QSEGroup(species_indices: [";
@@ -175,7 +198,8 @@ namespace gridfire {
                 }
                 count = 0;
                 os << "], is_in_equilibrium: " << group.is_in_equilibrium
-                   << ", algebraic_indices: [";
+                   << ", mean_timescale: " << group.mean_timescale
+                   << " s, algebraic_indices: [";
                 for (const auto& idx : group.algebraic_indices) {
                     os << idx;
                     if (count < group.algebraic_indices.size() - 1) {
@@ -242,6 +266,7 @@ namespace gridfire {
                 GenerateJacobianMatrix,
                 CalculateMolarReactionFlow,
                 GetSpeciesTimescales,
+                GetSpeciesDestructionTimescales,
                 Other,
                 All
             };
@@ -251,6 +276,7 @@ namespace gridfire {
                 {operators::GenerateJacobianMatrix, "generateJacobianMatrix"},
                 {operators::CalculateMolarReactionFlow, "calculateMolarReactionFlow"},
                 {operators::GetSpeciesTimescales, "getSpeciesTimescales"},
+                {operators::GetSpeciesDestructionTimescales, "getSpeciesDestructionTimescales"},
                 {operators::Other, "other"}
             };
 
@@ -262,6 +288,7 @@ namespace gridfire {
                 {operators::GenerateJacobianMatrix, 0},
                 {operators::CalculateMolarReactionFlow, 0},
                 {operators::GetSpeciesTimescales, 0},
+                {operators::GetSpeciesDestructionTimescales, 0},
                 {operators::Other, 0}
             };
 
@@ -351,11 +378,14 @@ namespace gridfire {
             double rho
         ) const;
 
+        std::unordered_map<size_t, std::vector<size_t>> buildConnectivityGraph(
+            const std::vector<size_t>& species_pool
+        ) const;
+
         std::vector<QSEGroup> constructCandidateGroups(
-            const std::vector<std::vector<size_t>>& timescale_pools,
+            const std::vector<std::vector<size_t>>& candidate_pools,
             const std::vector<double>& Y,
-            double T9,
+            double T9, double rho
-            double rho
         ) const;
     };
 }

src/network/include/gridfire/exceptions/error_engine.h

@@ -7,6 +7,58 @@
 namespace gridfire::exceptions {
     class EngineError : public std::exception {};
 
+    class StaleEngineTrigger final : public EngineError {
+    public:
+        struct state {
+            double m_T9;
+            double m_rho;
+            std::vector<double> m_Y;
+            double m_t;
+            int m_total_steps;
+            double m_eps_nuc;
+        };
+        explicit StaleEngineTrigger(const state &s)
+            :  m_state(s) {}
+
+        const char* what() const noexcept override{
+            return "Engine reports stale state. This means that the caller should trigger a update of the engine state before continuing with the integration. If you as an end user are seeing this error, it is likely a bug in the code that should be reported. Please provide the input parameters and the context in which this error occurred. Thank you for your help!";
+        }
+
+        state getState() const {
+            return m_state;
+        }
+
+        size_t numSpecies() const {
+            return m_state.m_Y.size();
+        }
+
+        size_t totalSteps() const {
+            return m_state.m_total_steps;
+        }
+
+        double energy() const {
+            return m_state.m_eps_nuc;
+        }
+
+        double getMolarAbundance(const size_t index) const {
+            if (index > m_state.m_Y.size() - 1) {
+                throw std::out_of_range("Index out of bounds for molar abundance vector.");
+            }
+            return m_state.m_Y[index];
+        }
+
+        double temperature() const {
+            return m_state.m_T9 * 1e9; // Convert T9 back to Kelvin
+        }
+
+        double density() const {
+            return m_state.m_rho;
+        }
+    private:
+        state m_state;
+
+    };
+
     class StaleEngineError final : public EngineError {
     public:
         explicit StaleEngineError(const std::string& message)
@@ -15,6 +67,7 @@ namespace gridfire::exceptions {
         const char* what() const noexcept override {
             return m_message.c_str();
         }
+
     private:
         std::string m_message;
     };
@@ -31,4 +84,29 @@ namespace gridfire::exceptions {
         std::string m_message;
     };
 
+    class NetworkResizedError final : public EngineError {
+    public:
+        explicit NetworkResizedError(const std::string& message)
+            : m_message(message) {}
+
+        const char* what() const noexcept override {
+            return m_message.c_str();
+        }
+    private:
+        std::string m_message;
+    };
+
+    class UnableToSetNetworkReactionsError final : public EngineError {
+    public:
+        explicit UnableToSetNetworkReactionsError(const std::string& message)
+            : m_message(message) {}
+
+        const char* what() const noexcept override {
+            return m_message.c_str();
+        }
+
+    private:
+        std::string m_message;
+    };
+
 }

src/network/include/gridfire/expectations/expected_engine.h

@@ -10,6 +10,10 @@ namespace gridfire::expectations {
         STALE
     };
 
+    enum class StaleEngineErrorTypes {
+        SYSTEM_RESIZED
+    };
+
     struct EngineError {
         std::string m_message;
         EngineErrorTypes type = EngineErrorTypes::FAILURE;
@@ -30,5 +34,17 @@ namespace gridfire::expectations {
 
     struct StaleEngineError : EngineError {
         EngineErrorTypes type = EngineErrorTypes::STALE;
+        StaleEngineErrorTypes staleType;
+
+        explicit StaleEngineError(StaleEngineErrorTypes staleType) : staleType(staleType) {}
+
+        explicit operator std::string() const {
+            switch (staleType) {
+                case (StaleEngineErrorTypes::SYSTEM_RESIZED):
+                    return "StaleEngineError: System resized, please update the engine.";
+                default:
+                    return "StaleEngineError: Unknown stale error type.";
+            }
+        }
     };
 }

src/network/include/gridfire/solver/solver.h

@@ -10,25 +10,9 @@
 
 #include "quill/Logger.h"
 
-#include "unsupported/Eigen/NonLinearOptimization" // Required for LevenbergMarquardt
-
 #include <vector>
 
 namespace gridfire::solver {
-
-    /**
-     * @struct dynamicQSESpeciesIndices
-     * @brief Structure to hold indices of dynamic and QSE species.
-     *
-     * This structure is used by the QSENetworkSolver to store the indices of species
-     * that are treated dynamically and those that are assumed to be in Quasi-Steady-State
-     * Equilibrium (QSE).
-     */
-    struct dynamicQSESpeciesIndices {
-        std::vector<size_t> dynamicSpeciesIndices; ///< Indices of slow species that are not in QSE.
-        std::vector<size_t> QSESpeciesIndices;  ///< Indices of fast species that are in QSE.
-    };
-
     /**
      * @class NetworkSolverStrategy
      * @brief Abstract base class for network solver strategies.
@@ -68,316 +52,6 @@ namespace gridfire::solver {
      */
     using DynamicNetworkSolverStrategy = NetworkSolverStrategy<DynamicEngine>;
 
-    /**
-     * @brief Type alias for a network solver strategy that uses an AdaptiveEngineView.
-     */
-    using AdaptiveNetworkSolverStrategy = NetworkSolverStrategy<AdaptiveEngineView>;
-
-    /**
-     * @brief Type alias for a network solver strategy that uses a static Engine.
-     */
-    using StaticNetworkSolverStrategy = NetworkSolverStrategy<Engine>;
-
-    /**
-     * @class QSENetworkSolver
-     * @brief A network solver that uses a Quasi-Steady-State Equilibrium (QSE) approach.
-     *
-     * This solver partitions the network into "fast" species in QSE and "slow" (dynamic) species.
-     * The abundances of the fast species are determined by solving a system of algebraic
-     * equations, while the abundances of the slow species are integrated using an ODE solver.
-     * This hybrid approach is highly effective for stiff networks with disparate timescales.
-     *
-     * The QSE solver uses an AdaptiveEngineView to dynamically cull unimportant species and
-     * reactions, which significantly improves performance for large networks.
-     *
-     * @implements AdaptiveNetworkSolverStrategy
-     *
-     * @see AdaptiveEngineView
-     * @see DynamicEngine::getSpeciesTimescales()
-     */
-    class QSENetworkSolver final : public DynamicNetworkSolverStrategy {
-    public:
-        /**
-         * @brief Constructor for the QSENetworkSolver.
-         * @param engine The adaptive engine view to use for evaluating the network.
-         */
-        using DynamicNetworkSolverStrategy::DynamicNetworkSolverStrategy;
-
-        /**
-         * @brief Evaluates the network for a given timestep using the QSE approach.
-         * @param netIn The input conditions for the network.
-         * @return The output conditions after the timestep.
-         *
-         * This method performs the following steps:
-         *   1. Updates the adaptive engine view (if necessary).
-         *   2. Partitions the species into dynamic and QSE species based on their timescales.
-         *   3. Calculates the steady-state abundances of the QSE species.
-         *   4. Integrates the ODEs for the dynamic species using a Runge-Kutta solver.
-         *   5. Marshals the output variables into a NetOut struct.
-         *
-         * @throws std::runtime_error If the steady-state abundances cannot be calculated.
-         *
-         * @see AdaptiveEngineView::update()
-         * @see packSpeciesTypeIndexVectors()
-         * @see calculateSteadyStateAbundances()
-         */
-        NetOut evaluate(const NetIn& netIn) override;
-    private: // methods
-        /**
-         * @brief Packs the species indices into vectors based on their type (dynamic or QSE).
-         * @param Y Vector of current abundances for all species.
-         * @param T9 Temperature in units of 10^9 K.
-         * @param rho Density in g/cm^3.
-         * @return A dynamicQSESpeciesIndices struct containing the indices of the dynamic and QSE species.
-         *
-         * This method determines whether each species should be treated dynamically or as
-         * being in QSE based on its timescale and abundance.  Species with short timescales
-         * or low abundances are assumed to be in QSE.
-         *
-         * @see DynamicEngine::getSpeciesTimescales()
-         */
-        dynamicQSESpeciesIndices packSpeciesTypeIndexVectors(
-            const std::vector<double>& Y,
-            const double T9,
-            const double rho
-        ) const;
-
-        /**
-         * @brief Calculates the steady-state abundances of the QSE species.
-         * @param Y Vector of current abundances for all species.
-         * @param T9 Temperature in units of 10^9 K.
-         * @param rho Density in g/cm^3.
-         * @param indices A dynamicQSESpeciesIndices struct containing the indices of the dynamic and QSE species.
-         * @return An Eigen::VectorXd containing the steady-state abundances of the QSE species.
-         *
-         * This method solves a system of algebraic equations to determine the steady-state
-         * abundances of the QSE species.
-         *
-         * @throws std::runtime_error If the steady-state abundances cannot be calculated.
-         */
-        Eigen::VectorXd calculateSteadyStateAbundances(
-            const std::vector<double>& Y,
-            const double T9,
-            const double rho,
-            const dynamicQSESpeciesIndices& indices
-        ) const;
-
-        /**
-         * @brief Initializes the network with a short ignition phase.
-         * @param netIn The input conditions for the network.
-         * @return The output conditions after the ignition phase.
-         *
-         * This method performs a short integration of the network at a high temperature and
-         * density to ignite the network and bring it closer to equilibrium.  This can improve
-         * the convergence of the QSE solver.
-         *
-         * @see DirectNetworkSolver::evaluate()
-         */
-        NetOut initializeNetworkWithShortIgnition(
-            const NetIn& netIn
-        ) const;
-
-        /**
-         * @brief Determines whether the adaptive engine view should be updated.
-         * @param conditions The current input conditions.
-         * @return True if the view should be updated, false otherwise.
-         *
-         * This method implements a policy for determining when the adaptive engine view
-         * should be updated.  The view is updated if the temperature or density has changed
-         * significantly, or if a primary fuel source has been depleted.
-         *
-         * @see AdaptiveEngineView::update()
-         */
-        bool shouldUpdateView(const NetIn& conditions) const;
-    private: // Nested functors for ODE integration
-        /**
-         * @struct RHSFunctor
-         * @brief Functor for calculating the right-hand side of the ODEs for the dynamic species.
-         *
-         * This functor is used by the ODE solver to calculate the time derivatives of the
-         * dynamic species.  It takes the current abundances of the dynamic species as input
-         * and returns the time derivatives of those abundances.
-         */
-        struct RHSFunctor {
-            DynamicEngine& m_engine; ///< The engine used to evaluate the network.
-            const std::vector<size_t>& m_dynamicSpeciesIndices; ///< Indices of the dynamic species.
-            const std::vector<size_t>& m_QSESpeciesIndices; ///< Indices of the QSE species.
-            const Eigen::VectorXd& m_Y_QSE; ///< Steady-state abundances of the QSE species.
-            const double m_T9; ///< Temperature in units of 10^9 K.
-            const double m_rho; ///< Density in g/cm^3.
-            quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); ///< Logger instance for trace and debug information.
-
-            bool m_isViewInitialized = false;
-
-            /**
-             * @brief Constructor for the RHSFunctor.
-             * @param engine The engine used to evaluate the network.
-             * @param dynamicSpeciesIndices Indices of the dynamic species.
-             * @param QSESpeciesIndices Indices of the QSE species.
-             * @param Y_QSE Steady-state abundances of the QSE species.
-             * @param T9 Temperature in units of 10^9 K.
-             * @param rho Density in g/cm^3.
-             */
-            RHSFunctor(
-                DynamicEngine& engine,
-                const std::vector<size_t>& dynamicSpeciesIndices,
-                const std::vector<size_t>& QSESpeciesIndices,
-                const Eigen::VectorXd& Y_QSE,
-                const double T9,
-                const double rho
-            ) :
-            m_engine(engine),
-            m_dynamicSpeciesIndices(dynamicSpeciesIndices),
-            m_QSESpeciesIndices(QSESpeciesIndices),
-            m_Y_QSE(Y_QSE),
-            m_T9(T9),
-            m_rho(rho) {}
-
-            /**
-             * @brief Calculates the time derivatives of the dynamic species.
-             * @param YDynamic Vector of current abundances for the dynamic species.
-             * @param dYdtDynamic Vector to store the time derivatives of the dynamic species.
-             * @param t Current time.
-             */
-            void operator()(
-                const boost::numeric::ublas::vector<double>& YDynamic,
-                boost::numeric::ublas::vector<double>& dYdtDynamic,
-                double t
-            ) const;
-
-        };
-
-        /**
-         * @struct JacobianFunctor
-         * @brief Functor for calculating the Jacobian matrix of the ODEs for the dynamic species.
-         *
-         * This functor is used by the ODE solver to calculate the Jacobian matrix of the
-         * ODEs for the dynamic species.  It takes the current abundances of the dynamic
-         * species as input and returns the Jacobian matrix.
-         *
-         * @todo Implement the Jacobian functor.
-         */
-        struct JacobianFunctor {
-            DynamicEngine& m_engine; ///< The engine used to evaluate the network.
-            const std::vector<size_t>& m_dynamicSpeciesIndices; ///< Indices of the dynamic species.
-            const std::vector<size_t>& m_QSESpeciesIndices; ///< Indices of the QSE species.
-            const double m_T9; ///< Temperature in units of 10^9 K.
-            const double m_rho; ///< Density in g/cm^3.
-
-            /**
-             * @brief Constructor for the JacobianFunctor.
-             * @param engine The engine used to evaluate the network.
-             * @param dynamicSpeciesIndices Indices of the dynamic species.
-             * @param QSESpeciesIndices Indices of the QSE species.
-             * @param T9 Temperature in units of 10^9 K.
-             * @param rho Density in g/cm^3.
-             */
-            JacobianFunctor(
-                DynamicEngine& engine,
-                const std::vector<size_t>& dynamicSpeciesIndices,
-                const std::vector<size_t>& QSESpeciesIndices,
-                const double T9,
-                const double rho
-            ) :
-            m_engine(engine),
-            m_dynamicSpeciesIndices(dynamicSpeciesIndices),
-            m_QSESpeciesIndices(QSESpeciesIndices),
-            m_T9(T9),
-            m_rho(rho) {}
-
-            /**
-             * @brief Calculates the Jacobian matrix of the ODEs for the dynamic species.
-             * @param YDynamic Vector of current abundances for the dynamic species.
-             * @param JDynamic Matrix to store the Jacobian matrix.
-             * @param t Current time.
-             * @param dfdt Vector to store the time derivatives of the dynamic species (not used).
-             */
-            void operator()(
-                const boost::numeric::ublas::vector<double>& YDynamic,
-                boost::numeric::ublas::matrix<double>& JDynamic,
-                double t,
-                boost::numeric::ublas::vector<double>& dfdt
-            ) const;
-        };
-
-        /**
-         * @struct EigenFunctor
-         * @brief Functor for calculating the residual and Jacobian for the QSE species using Eigen.
-         *
-         * @tparam T The numeric type to use for the calculation (double or ADDouble).
-         */
-        template<typename T>
-        struct EigenFunctor {
-            using InputType = Eigen::Matrix<T, Eigen::Dynamic, 1>;
-            using OutputType = Eigen::Matrix<T, Eigen::Dynamic, 1>;
-            using JacobianType = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;
-            enum {
-                InputsAtCompileTime = Eigen::Dynamic,
-                ValuesAtCompileTime = Eigen::Dynamic
-            };
-
-            DynamicEngine& m_engine; ///< The engine used to evaluate the network.
-            quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); ///< Logger instance for trace and debug information.
-            const std::vector<double>& m_YFull; ///< The full, initial abundance vector
-            const std::vector<size_t>& m_dynamicSpeciesIndices; ///< Indices of the dynamic species.
-            const std::vector<size_t>& m_QSESpeciesIndices; ///< Indices of the QSE species.
-            const double m_T9; ///< Temperature in units of 10^9 K.
-            const double m_rho; ///< Density in g/cm^3.
-            const Eigen::VectorXd& m_YScale; ///< Scaling vector for the QSE species for asinh scaling.
-
-            /**
-             * @brief Constructor for the EigenFunctor.
-             * @param engine The engine used to evaluate the network.
-             * @param YFull Abundances of the dynamic species.
-             * @param dynamicSpeciesIndices Indices of the dynamic species.
-             * @param QSESpeciesIndices Indices of the QSE species.
-             * @param T9 Temperature in units of 10^9 K.
-             * @param rho Density in g/cm^3.
-             */
-            EigenFunctor(
-                DynamicEngine& engine,
-                const std::vector<double>& YFull,
-                const std::vector<size_t>& dynamicSpeciesIndices,
-                const std::vector<size_t>& QSESpeciesIndices,
-                const double T9,
-                const double rho,
-                const Eigen::VectorXd& YScale
-            ) :
-            m_engine(engine),
-            m_YFull(YFull),
-            m_dynamicSpeciesIndices(dynamicSpeciesIndices),
-            m_QSESpeciesIndices(QSESpeciesIndices),
-            m_T9(T9),
-            m_rho(rho),
-            m_YScale(YScale) {}
-
-            int values() const { return m_QSESpeciesIndices.size(); }
-            int inputs() const { return m_QSESpeciesIndices.size(); }
-
-            /**
-             * @brief Calculates the residual vector for the QSE species.
-             * @param v_QSE_log Input vector of QSE species abundances (logarithmic).
-             * @param f_QSE Output vector of residuals.
-             * @return 0 for success.
-             */
-            int operator()(const InputType& v_QSE_log, OutputType& f_QSE) const;
-
-            /**
-             * @brief Calculates the Jacobian matrix for the QSE species.
-             * @param v_QSE_log Input vector of QSE species abundances (logarithmic).
-             * @param J_QSE Output Jacobian matrix.
-             * @return 0 for success.
-             */
-            int df(const InputType& v_QSE_log, JacobianType& J_QSE) const;
-        };
-    private:
-        quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); ///< Logger instance.
-        fourdst::config::Config& m_config = fourdst::config::Config::getInstance(); ///< Configuration instance.
-
-        bool m_isViewInitialized = false; ///< Flag indicating whether the adaptive engine view has been initialized.
-        NetIn m_lastSeenConditions; ///< The last seen input conditions.
-    };
-
     /**
      * @class DirectNetworkSolver
      * @brief A network solver that directly integrates the reaction network ODEs.
@@ -411,12 +85,20 @@ namespace gridfire::solver {
          * species abundances.  It takes the current abundances as input and returns the
          * time derivatives.
          */
-        struct RHSFunctor {
+        struct RHSManager {
             DynamicEngine& m_engine; ///< The engine used to evaluate the network.
             const double m_T9; ///< Temperature in units of 10^9 K.
             const double m_rho; ///< Density in g/cm^3.
-            const size_t m_numSpecies; ///< The number of species in the network.
+
-            quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); ///< Logger instance.
+            mutable double m_cached_time;
+            mutable std::optional<StepDerivatives<double>> m_cached_result;
+
+            mutable double m_last_observed_time = 0.0; ///< Last time the state was observed.
+
+
+            quill::Logger* m_logger = LogManager::getInstance().newFileLogger("integration.log", "GridFire"); ///< Logger instance.
+            mutable int m_num_steps = 0;
+            mutable double m_last_step_time = 1e-20;
 
             /**
              * @brief Constructor for the RHSFunctor.
@@ -424,7 +106,7 @@ namespace gridfire::solver {
              * @param T9 Temperature in units of 10^9 K.
              * @param rho Density in g/cm^3.
              */
-            RHSFunctor(
+            RHSManager(
                 DynamicEngine& engine,
                 const double T9,
                 const double rho
@@ -432,7 +114,7 @@ namespace gridfire::solver {
             m_engine(engine),
             m_T9(T9),
             m_rho(rho),
-            m_numSpecies(engine.getNetworkSpecies().size()) {}
+            m_cached_time(0) {}
 
             /**
              * @brief Calculates the time derivatives of the species abundances.
@@ -445,6 +127,10 @@ namespace gridfire::solver {
                 boost::numeric::ublas::vector<double>& dYdt,
                 double t
             ) const;
+
+            void observe(const boost::numeric::ublas::vector<double>& state, double t) const;
+            void compute_and_cache(const boost::numeric::ublas::vector<double>& state, double t) const;
+
         };
 
         /**
@@ -458,7 +144,6 @@ namespace gridfire::solver {
             DynamicEngine& m_engine; ///< The engine used to evaluate the network.
             const double m_T9; ///< Temperature in units of 10^9 K.
             const double m_rho; ///< Density in g/cm^3.
-            const size_t m_numSpecies; ///< The number of species in the network.
 
             /**
              * @brief Constructor for the JacobianFunctor.
@@ -473,8 +158,7 @@ namespace gridfire::solver {
             ) :
             m_engine(engine),
             m_T9(T9),
-            m_rho(rho),
+            m_rho(rho) {}
-            m_numSpecies(engine.getNetworkSpecies().size()) {}
 
             /**
              * @brief Calculates the Jacobian matrix.
@@ -493,97 +177,7 @@ namespace gridfire::solver {
         };
 
     private:
-        quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); ///< Logger instance.
+        quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); ///< Logger instance.
-        fourdst::config::Config& m_config = fourdst::config::Config::getInstance(); ///< Configuration instance.
+        Config& m_config = Config::getInstance(); ///< Configuration instance.
     };
-
-    template<typename T>
-    int QSENetworkSolver::EigenFunctor<T>::operator()(const InputType &v_QSE, OutputType &f_QSE) const {
-        std::vector<double> y = m_YFull; // Full vector of species abundances
-        Eigen::VectorXd Y_QSE = m_YScale.array() * v_QSE.array().sinh(); // Convert to physical abundances using asinh scaling
-        LOG_DEBUG(
-            m_logger,
-            "Calling QSENetworkSolver::EigenFunctor::operator() with QSE abundances {}",
-            [&]() -> std::string {
-                std::stringstream ss;
-                ss << std::scientific << std::setprecision(5);
-                for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
-                    ss << std::string(m_engine.getNetworkSpecies()[m_QSESpeciesIndices[i]].name()) << ": " << Y_QSE(i) << ", ";
-                }
-                return ss.str();
-            }());
-
-        for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
-            y[m_QSESpeciesIndices[i]] = Y_QSE(i);
-        }
-
-
-        const auto [dydt, specificEnergyRate] = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho);
-        LOG_DEBUG(
-            m_logger,
-            "Calling QSENetworkSolver::EigenFunctor::operator() found QSE dydt {}",
-            [&]() -> std::string {
-                std::stringstream ss;
-                ss << std::scientific << std::setprecision(5);
-                for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
-                    ss << std::string(m_engine.getNetworkSpecies()[m_QSESpeciesIndices[i]].name()) << ": " << dydt[m_QSESpeciesIndices[i]] << ", ";
-                }
-                return ss.str();
-            }());
-        f_QSE.resize(m_QSESpeciesIndices.size());
-        for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
-            f_QSE(i) = dydt[m_QSESpeciesIndices[i]];
-        }
-        return 0; // Success
-    }
-
-    template <typename T>
-    int QSENetworkSolver::EigenFunctor<T>::df(const InputType& v_QSE, JacobianType& J_QSE) const {
-        std::vector<double> y = m_YFull;
-        Eigen::VectorXd Y_QSE = m_YScale.array() * v_QSE.array().sinh();
-
-        for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
-            y[m_QSESpeciesIndices[i]] = Y_QSE(i);
-        }
-
-        m_engine.generateJacobianMatrix(y, m_T9, m_rho);
-
-        J_QSE.resize(m_QSESpeciesIndices.size(), m_QSESpeciesIndices.size());
-        for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
-            for (size_t j = 0; j < m_QSESpeciesIndices.size(); ++j) {
-                J_QSE(i, j) = m_engine.getJacobianMatrixEntry(m_QSESpeciesIndices[i], m_QSESpeciesIndices[j]);
-            }
-        }
-
-        std::string format_jacobian = [&]() -> std::string {
-            std::stringstream ss;
-            ss << std::scientific << std::setprecision(5);
-            for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
-                for (size_t j = 0; j < m_QSESpeciesIndices.size(); ++j) {
-                    ss << J_QSE(i, j) << " ";
-                }
-                ss << "\n";
-            }
-            return ss.str();
-        }();
-
-
-        for (long j = 0; j < J_QSE.cols(); ++j) {
-            LOG_DEBUG(
-                m_logger,
-                "Jacobian ({},{}) before chain rule: {}",
-                j, j, J_QSE(j, j)
-            );
-            const double dY_dv = m_YScale(j) * CppAD::cosh(v_QSE(j));
-            J_QSE.col(j) *= dY_dv; // chain rule for log space transformation
-            LOG_DEBUG(
-                m_logger,
-                "Jacobian ({},{}) after chain rule: {} (dY/dv = {})",
-                j, j, J_QSE(j, j), dY_dv
-            );
-        }
-        return 0;
-    }
-
-
 }

src/network/lib/engine/engine_graph.cpp

@@ -54,7 +54,7 @@ namespace gridfire {
         syncInternalMaps();
     }
 
-    StepDerivatives<double> GraphEngine::calculateRHSAndEnergy(
+    std::expected<StepDerivatives<double>, expectations::StaleEngineError> GraphEngine::calculateRHSAndEnergy(
         const std::vector<double> &Y,
         const double T9,
         const double rho
@@ -91,8 +91,8 @@ namespace gridfire {
         const size_t n = m_rhsADFun.Domain();
         const size_t m = m_rhsADFun.Range();
 
-        std::vector<bool> select_domain(n, true);
+        const std::vector<bool> select_domain(n, true);
-        std::vector<bool> select_range(m, true);
+        const std::vector<bool> select_range(m, true);
 
         m_rhsADFun.subgraph_sparsity(select_domain, select_range, false, m_full_jacobian_sparsity_pattern);
         m_jac_work.clear();
@@ -148,10 +148,10 @@ namespace gridfire {
         }
     }
 
-    void GraphEngine::reserveJacobianMatrix() {
+    void GraphEngine::reserveJacobianMatrix() const {
         // The implementation of this function (and others) constrains this nuclear network to a constant temperature and density during
         // each evaluation.
-        size_t numSpecies = m_networkSpecies.size();
+        const size_t numSpecies = m_networkSpecies.size();
         m_jacobianMatrix.clear();
         m_jacobianMatrix.resize(numSpecies, numSpecies, false); // Sparse matrix, no initial values
         LOG_TRACE_L2(m_logger, "Jacobian matrix resized to {} rows and {} columns.",
@@ -171,6 +171,11 @@ namespace gridfire {
         return m_reactions;
     }
 
+    void GraphEngine::setNetworkReactions(const reaction::LogicalReactionSet &reactions) {
+        m_reactions = reactions;
+        syncInternalMaps();
+    }
+
     bool GraphEngine::involvesSpecies(const fourdst::atomic::Species& species) const {
         // Checks if a given species is present in the network's species map for efficient lookup.
         const bool found = m_networkSpeciesMap.contains(species.name());
@@ -259,7 +264,7 @@ namespace gridfire {
         const double T9
     ) const {
         if (!m_useReverseReactions) {
-            LOG_TRACE_L2(m_logger, "Reverse reactions are disabled. Returning 0.0 for reverse rate of reaction '{}'.", reaction.id());
+            LOG_TRACE_L3_LIMIT_EVERY_N(std::numeric_limits<int>::max(), m_logger, "Reverse reactions are disabled. Returning 0.0 for reverse rate of reaction '{}'.", reaction.id());
             return 0.0; // If reverse reactions are not used, return 0.0
         }
         const double expFactor = std::exp(-reaction.qValue() / (m_constants.kB * T9));
@@ -279,7 +284,7 @@ namespace gridfire {
         } else {
             LOG_WARNING_LIMIT_EVERY_N(1000000, m_logger, "Reverse rate calculation for reactions with more than two reactants or products is not implemented (reaction id {}).", reaction.peName());
         }
-        LOG_TRACE_L2(m_logger, "Calculated reverse rate for reaction '{}': {:.3E} at T9={:.3E}.", reaction.id(), reverseRate, T9);
+        LOG_TRACE_L2_LIMIT_EVERY_N(1000, m_logger, "Calculated reverse rate for reaction '{}': {:.3E} at T9={:.3E}.", reaction.id(), reverseRate, T9);
         return reverseRate;
     }
 
@@ -352,7 +357,7 @@ namespace gridfire {
         const double CT = std::pow(massNumerator/massDenominator, 1.5) *
             (partitionFunctionNumerator/partitionFunctionDenominator);
 
-        double reverseRate = forwardRate * symmetryFactor * CT * expFactor;
+        const double reverseRate = forwardRate * symmetryFactor * CT * expFactor;
         if (!std::isfinite(reverseRate)) {
             return 0.0; // If the reverse rate is not finite, return 0.0
         }
@@ -366,7 +371,7 @@ namespace gridfire {
         const double reverseRate
     ) const {
         if (!m_useReverseReactions) {
-            LOG_TRACE_L2(m_logger, "Reverse reactions are disabled. Returning 0.0 for reverse rate derivative of reaction '{}'.", reaction.id());
+            LOG_TRACE_L3_LIMIT_EVERY_N(std::numeric_limits<int>::max(), m_logger, "Reverse reactions are disabled. Returning 0.0 for reverse rate of reaction '{}'.", reaction.id());
             return 0.0; // If reverse reactions are not used, return 0.0
         }
         const double d_log_kFwd = reaction.calculate_forward_rate_log_derivative(T9);
@@ -652,7 +657,7 @@ namespace gridfire {
         const double rho
     ) const {
 
-        LOG_TRACE_L1(m_logger, "Generating jacobian matrix for T9={}, rho={}..", T9, rho);
+        LOG_TRACE_L1_LIMIT_EVERY_N(1000, m_logger, "Generating jacobian matrix for T9={}, rho={}..", T9, rho);
         const size_t numSpecies = m_networkSpecies.size();
 
         // 1. Pack the input variables into a vector for CppAD
@@ -662,20 +667,20 @@ namespace gridfire {
         }
         adInput[numSpecies]     = T9;  // T9
         adInput[numSpecies + 1] = rho; // rho
-        LOG_DEBUG(
+        // LOG_DEBUG(
-            m_logger,
+        //     m_logger,
-            "AD Input to jacobian {}",
+        //     "AD Input to jacobian {}",
-            [&]() -> std::string {
+        //     [&]() -> std::string {
-                std::stringstream ss;
+        //         std::stringstream ss;
-                ss << std::scientific << std::setprecision(5);
+        //         ss << std::scientific << std::setprecision(5);
-                for (size_t i = 0; i < adInput.size(); ++i) {
+        //         for (size_t i = 0; i < adInput.size(); ++i) {
-                    ss << adInput[i];
+        //             ss << adInput[i];
-                    if (i < adInput.size() - 1) {
+        //             if (i < adInput.size() - 1) {
-                        ss << ", ";
+        //                 ss << ", ";
-                    }
+        //             }
-                }
+        //         }
-                return ss.str();
+        //         return ss.str();
-            }());
+        //     }());
 
         // 2. Calculate the full jacobian
         const std::vector<double> dotY = m_rhsADFun.Jacobian(adInput);
@@ -690,32 +695,32 @@ namespace gridfire {
                 }
             }
         }
-        LOG_DEBUG(
+        // LOG_DEBUG(
-            m_logger,
+        //     m_logger,
-            "Final Jacobian is:\n{}",
+        //     "Final Jacobian is:\n{}",
-            [&]() -> std::string {
+        //     [&]() -> std::string {
-                std::stringstream ss;
+        //         std::stringstream ss;
-                ss << std::scientific << std::setprecision(5);
+        //         ss << std::scientific << std::setprecision(5);
-                for (size_t i = 0; i < m_jacobianMatrix.size1(); ++i) {
+        //         for (size_t i = 0; i < m_jacobianMatrix.size1(); ++i) {
-                    ss << getNetworkSpecies()[i].name();
+        //             ss << getNetworkSpecies()[i].name();
-                    if (i < m_jacobianMatrix.size1() - 1) {
+        //             if (i < m_jacobianMatrix.size1() - 1) {
-                        ss << ", ";
+        //                 ss << ", ";
-                    }
+        //             }
-                }
+        //         }
-                ss << "\n";
+        //         ss << "\n";
-                for (size_t i = 0; i < m_jacobianMatrix.size1(); ++i) {
+        //         for (size_t i = 0; i < m_jacobianMatrix.size1(); ++i) {
-                    ss << getNetworkSpecies()[i].name() << ": ";
+        //             ss << getNetworkSpecies()[i].name() << ": ";
-                    for (size_t j = 0; j < m_jacobianMatrix.size2(); ++j) {
+        //             for (size_t j = 0; j < m_jacobianMatrix.size2(); ++j) {
-                        ss << m_jacobianMatrix(i, j);
+        //                 ss << m_jacobianMatrix(i, j);
-                        if (j < m_jacobianMatrix.size2() - 1) {
+        //                 if (j < m_jacobianMatrix.size2() - 1) {
-                            ss << ", ";
+        //                     ss << ", ";
-                        }
+        //                 }
-                    }
+        //             }
-                    ss << "\n";
+        //             ss << "\n";
-                }
+        //         }
-                return ss.str();
+        //         return ss.str();
-            }());
+        //     }());
-        LOG_TRACE_L1(m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2());
+        LOG_TRACE_L1_LIMIT_EVERY_N(1000, m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2());
     }
 
     void GraphEngine::generateJacobianMatrix(
@@ -895,7 +900,7 @@ namespace gridfire {
         LOG_TRACE_L1(m_logger, "Successfully exported network graph to {}", filename);
     }
 
-    std::unordered_map<fourdst::atomic::Species, double> GraphEngine::getSpeciesTimescales(
+    std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> GraphEngine::getSpeciesTimescales(
         const std::vector<double> &Y,
         const double T9,
         const double rho
@@ -914,8 +919,41 @@ namespace gridfire {
         return speciesTimescales;
     }
 
+    std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> GraphEngine::getSpeciesDestructionTimescales(
+        const std::vector<double> &Y,
+        const double T9,
+        const double rho
+    ) const {
+        auto [dydt, _] = calculateAllDerivatives<double>(Y, T9, rho);
+        std::unordered_map<fourdst::atomic::Species, double> speciesDestructionTimescales;
+        speciesDestructionTimescales.reserve(m_networkSpecies.size());
+        for (const auto& species : m_networkSpecies) {
+            double netDestructionFlow = 0.0;
+            for (const auto& reaction : m_reactions) {
+                if (reaction.stoichiometry(species) < 0) {
+                    const double flow = calculateMolarReactionFlow<double>(reaction, Y, T9, rho);
+                    netDestructionFlow += flow;
+                }
+            }
+            double timescale = std::numeric_limits<double>::infinity();
+            if (netDestructionFlow != 0.0) {
+                timescale = Y[getSpeciesIndex(species)] / netDestructionFlow;
+            }
+            speciesDestructionTimescales.emplace(species, timescale);
+        }
+        return speciesDestructionTimescales;
+    }
+
     fourdst::composition::Composition GraphEngine::update(const NetIn &netIn) {
-        return netIn.composition;
+        fourdst::composition::Composition baseUpdatedComposition = netIn.composition;
+        for (const auto& species : m_networkSpecies) {
+            if (!netIn.composition.contains(species)) {
+                baseUpdatedComposition.registerSpecies(species);
+                baseUpdatedComposition.setMassFraction(species, 0.0);
+            }
+        }
+        baseUpdatedComposition.finalize(false);
+        return baseUpdatedComposition;
     }
 
     bool GraphEngine::isStale(const NetIn &netIn) {

src/network/lib/engine/procedures/priming.cpp

@@ -45,7 +45,7 @@ namespace gridfire {
                 speciesToPrime.push_back(entry.isotope());
             }
         }
-        LOG_INFO(logger, "Priming {} species in the network.", speciesToPrime.size());
+        LOG_DEBUG(logger, "Priming {} species in the network.", speciesToPrime.size());
 
         PrimingReport report;
         if (speciesToPrime.empty()) {
@@ -57,7 +57,7 @@ namespace gridfire {
 
         const double T9 = netIn.temperature / 1e9;
         const double rho = netIn.density;
-        const auto initialDepth = engine.getDepth();
+        const auto initialReactionSet = engine.getNetworkReactions();
 
         report.status = PrimingReportStatus::FULL_SUCCESS;
         report.success = true;
@@ -76,7 +76,7 @@ namespace gridfire {
         engine.rebuild(netIn.composition, NetworkBuildDepth::Full);
 
         for (const auto& primingSpecies : speciesToPrime) {
-            LOG_DEBUG(logger, "Priming species: {}", primingSpecies.name());
+            LOG_TRACE_L3(logger, "Priming species: {}", primingSpecies.name());
 
             // Create a temporary composition from the current internal state for the primer
             Composition tempComp;
@@ -114,12 +114,12 @@ namespace gridfire {
                     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_INFO(logger, "Found equilibrium for {}: X_eq = {:.4e}", primingSpecies.name(), equilibriumMassFraction);
+                LOG_TRACE_L3(logger, "Found equilibrium for {}: X_eq = {:.4e}", primingSpecies.name(), equilibriumMassFraction);
 
                 const reaction::Reaction* dominantChannel = findDominantCreationChannel(primer, primingSpecies, Y, T9, rho);
 
                 if (dominantChannel) {
-                    LOG_DEBUG(logger, "Dominant creation channel for {}: {}", primingSpecies.name(), dominantChannel->peName());
+                    LOG_TRACE_L3(logger, "Dominant creation channel for {}: {}", primingSpecies.name(), dominantChannel->peName());
 
                     double totalReactantMass = 0.0;
                     for (const auto& reactant : dominantChannel->reactants()) {
@@ -186,7 +186,7 @@ namespace gridfire {
             report.massFractionChanges.emplace_back(species, change);
         }
 
-        engine.rebuild(netIn.composition, initialDepth);
+        engine.setNetworkReactions(initialReactionSet);
         return report;
     }
 

src/network/lib/engine/views/engine_adaptive.cpp

@@ -6,6 +6,7 @@
 
 
 #include "gridfire/network.h"
+#include "gridfire/exceptions/error_engine.h"
 
 #include "quill/LogMacros.h"
 #include "quill/Logger.h"
@@ -86,11 +87,21 @@ namespace gridfire {
 
     fourdst::composition::Composition AdaptiveEngineView::update(const NetIn &netIn) {
         fourdst::composition::Composition baseUpdatedComposition = m_baseEngine.update(netIn);
+        NetIn updatedNetIn = netIn;
+
+        // for (const auto &entry: netIn.composition | std::views::values) {
+        //     if (baseUpdatedComposition.contains(entry.isotope())) {
+        //         updatedNetIn.composition.setMassFraction(entry.isotope(), baseUpdatedComposition.getMassFraction(entry.isotope()));
+        //     }
+        // }
+        updatedNetIn.composition = baseUpdatedComposition;
+
+        updatedNetIn.composition.finalize(false);
 
         LOG_TRACE_L1(m_logger, "Updating AdaptiveEngineView with new network input...");
 
         std::vector<double> Y_Full;
-        std::vector<ReactionFlow> allFlows = calculateAllReactionFlows(netIn, Y_Full);
+        std::vector<ReactionFlow> allFlows = calculateAllReactionFlows(updatedNetIn, Y_Full);
 
         double maxFlow = 0.0;
 
@@ -101,24 +112,24 @@ namespace gridfire {
         }
         LOG_DEBUG(m_logger, "Maximum reaction flow rate in adaptive engine view: {:0.3E} [mol/s]", maxFlow);
 
-        const std::unordered_set<Species> reachableSpecies = findReachableSpecies(netIn);
+        const std::unordered_set<Species> reachableSpecies = findReachableSpecies(updatedNetIn);
         LOG_DEBUG(m_logger, "Found {} reachable species in adaptive engine view.", reachableSpecies.size());
 
         const std::vector<const reaction::LogicalReaction*> finalReactions = cullReactionsByFlow(allFlows, reachableSpecies, Y_Full, maxFlow);
 
         finalizeActiveSet(finalReactions);
 
-        // auto [rescuedReactions, rescuedSpecies] = rescueEdgeSpeciesDestructionChannel(allFlows, Y_Full, netIn.temperature/1e9, netIn.density, m_activeSpecies, m_activeReactions);
+        auto [rescuedReactions, rescuedSpecies] = rescueEdgeSpeciesDestructionChannel(Y_Full, netIn.temperature/1e9, netIn.density, m_activeSpecies, m_activeReactions);
-        //
+
-        // for (const auto& reactionPtr : rescuedReactions) {
+        for (const auto& reactionPtr : rescuedReactions) {
-        //     m_activeReactions.add_reaction(*reactionPtr);
+            m_activeReactions.add_reaction(*reactionPtr);
-        // }
+        }
-        //
+
-        // for (const auto& species : rescuedSpecies) {
+        for (const auto& species : rescuedSpecies) {
-        //     if (!std::ranges::contains(m_activeSpecies, species)) {
+            if (!std::ranges::contains(m_activeSpecies, species)) {
-        //         m_activeSpecies.push_back(species);
+                m_activeSpecies.push_back(species);
-        //     }
+            }
-        // }
+        }
 
         m_speciesIndexMap = constructSpeciesIndexMap();
         m_reactionIndexMap = constructReactionIndexMap();
@@ -127,18 +138,18 @@ namespace gridfire {
 
         LOG_INFO(m_logger, "AdaptiveEngineView updated successfully with {} active species and {} active reactions.", m_activeSpecies.size(), m_activeReactions.size());
 
-        return baseUpdatedComposition;
+        return updatedNetIn.composition;
     }
 
     bool AdaptiveEngineView::isStale(const NetIn &netIn) {
-        return m_isStale;
+        return m_isStale || m_baseEngine.isStale(netIn);
     }
 
     const std::vector<Species> & AdaptiveEngineView::getNetworkSpecies() const {
         return m_activeSpecies;
     }
 
-    StepDerivatives<double> AdaptiveEngineView::calculateRHSAndEnergy(
+    std::expected<StepDerivatives<double>, expectations::StaleEngineError> AdaptiveEngineView::calculateRHSAndEnergy(
         const std::vector<double> &Y_culled,
         const double T9,
         const double rho
@@ -147,8 +158,13 @@ namespace gridfire {
 
         const auto Y_full = mapCulledToFull(Y_culled);
 
-        const auto [dydt, nuclearEnergyGenerationRate] = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho);
+        auto result = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho);
 
+        if (!result) {
+            return std::unexpected{result.error()};
+        }
+
+        const auto [dydt, nuclearEnergyGenerationRate] = result.value();
         StepDerivatives<double> culledResults;
         culledResults.nuclearEnergyGenerationRate = nuclearEnergyGenerationRate;
         culledResults.dydt = mapFullToCulled(dydt);
@@ -214,14 +230,26 @@ namespace gridfire {
         return m_activeReactions;
     }
 
-    std::unordered_map<Species, double> AdaptiveEngineView::getSpeciesTimescales(
+    void AdaptiveEngineView::setNetworkReactions(const reaction::LogicalReactionSet &reactions) {
+        LOG_CRITICAL(m_logger, "AdaptiveEngineView does not support setting network reactions directly. Use update() with NetIn instead. Perhaps you meant to call this on the base engine?");
+        throw exceptions::UnableToSetNetworkReactionsError("AdaptiveEngineView does not support setting network reactions directly. Use update() with NetIn instead. Perhaps you meant to call this on the base engine?");
+    }
+
+    std::expected<std::unordered_map<Species, double>, expectations::StaleEngineError> AdaptiveEngineView::getSpeciesTimescales(
         const std::vector<double> &Y_culled,
         const double T9,
         const double rho
     ) const {
         validateState();
         const auto Y_full = mapCulledToFull(Y_culled);
-        const auto fullTimescales = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho);
+        const auto result = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho);
+
+        if (!result) {
+            return std::unexpected{result.error()};
+        }
+
+        const std::unordered_map<Species, double> fullTimescales = result.value();
+
 
         std::unordered_map<Species, double> culledTimescales;
         culledTimescales.reserve(m_activeSpecies.size());
@@ -234,6 +262,33 @@ namespace gridfire {
 
     }
 
+    std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError>
+    AdaptiveEngineView::getSpeciesDestructionTimescales(
+        const std::vector<double> &Y,
+        double T9,
+        double rho
+    ) const {
+        validateState();
+
+        const auto Y_full = mapCulledToFull(Y);
+        const auto result = m_baseEngine.getSpeciesDestructionTimescales(Y_full, T9, rho);
+
+        if (!result) {
+            return std::unexpected{result.error()};
+        }
+
+        const std::unordered_map<Species, double> destructionTimescales = result.value();
+
+        std::unordered_map<Species, double> culledTimescales;
+        culledTimescales.reserve(m_activeSpecies.size());
+        for (const auto& active_species : m_activeSpecies) {
+            if (destructionTimescales.contains(active_species)) {
+                culledTimescales[active_species] = destructionTimescales.at(active_species);
+            }
+        }
+        return culledTimescales;
+    }
+
     void AdaptiveEngineView::setScreeningModel(const screening::ScreeningType model) {
         m_baseEngine.setScreeningModel(model);
     }
@@ -255,7 +310,7 @@ namespace gridfire {
     }
 
     int AdaptiveEngineView::getSpeciesIndex(const fourdst::atomic::Species &species) const {
-        auto it = std::find(m_activeSpecies.begin(), m_activeSpecies.end(), species);
+        const auto it = std::ranges::find(m_activeSpecies, species);
         if (it != m_activeSpecies.end()) {
             return static_cast<int>(std::distance(m_activeSpecies.begin(), it));
         } else {
@@ -428,14 +483,18 @@ namespace gridfire {
     }
 
     AdaptiveEngineView::RescueSet AdaptiveEngineView::rescueEdgeSpeciesDestructionChannel(
-        const std::vector<ReactionFlow> &allFlows,
         const std::vector<double> &Y_full,
         const double T9,
         const double rho,
         const std::vector<Species> &activeSpecies,
         const reaction::LogicalReactionSet &activeReactions
     ) const {
-        const auto timescales = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho);
+        const auto result = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho);
+        if (!result) {
+            LOG_ERROR(m_logger, "Failed to get species timescales due to stale engine state.");
+            throw exceptions::StaleEngineError("Failed to get species timescales");
+        }
+        std::unordered_map<Species, double> timescales = result.value();
         std::set<Species> onlyProducedSpecies;
         for (const auto& reaction : activeReactions) {
             const std::vector<Species> products = reaction.products();
@@ -521,7 +580,7 @@ namespace gridfire {
                 }
                 int count = 0;
                 ss << "(";
-                for (const auto& [species, reaction] : reactionsToRescue) {
+                for (const auto &reaction : reactionsToRescue | std::views::values) {
                     ss << reaction->id();
                     if (count < reactionsToRescue.size() - 1) {
                         ss << ", ";

src/network/lib/engine/views/engine_defined.cpp

@@ -16,60 +16,7 @@ namespace gridfire {
 
     DefinedEngineView::DefinedEngineView(const std::vector<std::string>& peNames, DynamicEngine& baseEngine) :
     m_baseEngine(baseEngine) {
-        m_activeReactions.clear();
+        collect(peNames);
-        m_activeSpecies.clear();
-
-        std::unordered_set<Species> seenSpecies;
-
-        const auto& fullNetworkReactionSet = m_baseEngine.getNetworkReactions();
-        for (const auto& peName : peNames) {
-            if (!fullNetworkReactionSet.contains(peName)) {
-                LOG_ERROR(m_logger, "Reaction with name '{}' not found in the base engine's network reactions. Aborting...", peName);
-                m_logger->flush_log();
-                throw std::runtime_error("Reaction with name '" + std::string(peName) + "' not found in the base engine's network reactions.");
-            }
-            auto reaction = fullNetworkReactionSet[peName];
-            for (const auto& reactant : reaction.reactants()) {
-                if (!seenSpecies.contains(reactant)) {
-                    seenSpecies.insert(reactant);
-                    m_activeSpecies.push_back(reactant);
-                }
-            }
-            for (const auto& product : reaction.products()) {
-                if (!seenSpecies.contains(product)) {
-                    seenSpecies.insert(product);
-                    m_activeSpecies.push_back(product);
-                }
-            }
-            m_activeReactions.add_reaction(reaction);
-        }
-        LOG_TRACE_L1(m_logger, "DefinedEngineView built with {} active species and {} active reactions.", m_activeSpecies.size(), m_activeReactions.size());
-        LOG_DEBUG(m_logger, "Active species: {}", [this]() -> std::string {
-            std::string result;
-            for (const auto& species : m_activeSpecies) {
-                result += std::string(species.name()) + ", ";
-            }
-            if (!result.empty()) {
-                result.pop_back(); // Remove last space
-                result.pop_back(); // Remove last comma
-            }
-            return result;
-        }());
-        LOG_DEBUG(m_logger, "Active reactions: {}", [this]() -> std::string {
-            std::string result;
-            for (const auto& reaction : m_activeReactions) {
-                result += std::string(reaction.id()) + ", ";
-            }
-            if (!result.empty()) {
-                result.pop_back(); // Remove last space
-                result.pop_back(); // Remove last comma
-            }
-            return result;
-        }());
-        m_speciesIndexMap = constructSpeciesIndexMap();
-        m_reactionIndexMap = constructReactionIndexMap();
-        m_isStale = false;
-
     }
 
     const DynamicEngine & DefinedEngineView::getBaseEngine() const {
@@ -80,7 +27,7 @@ namespace gridfire {
         return m_activeSpecies;
     }
 
-    StepDerivatives<double> DefinedEngineView::calculateRHSAndEnergy(
+    std::expected<StepDerivatives<double>, expectations::StaleEngineError> DefinedEngineView::calculateRHSAndEnergy(
         const std::vector<double> &Y_defined,
         const double T9,
         const double rho
@@ -88,8 +35,13 @@ namespace gridfire {
         validateNetworkState();
 
         const auto Y_full = mapViewToFull(Y_defined);
-        const auto [dydt, nuclearEnergyGenerationRate] = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho);
+        const auto result = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho);
 
+        if (!result) {
+            return std::unexpected{result.error()};
+        }
+
+        const auto [dydt, nuclearEnergyGenerationRate] = result.value();
         StepDerivatives<double> definedResults;
         definedResults.nuclearEnergyGenerationRate = nuclearEnergyGenerationRate;
         definedResults.dydt = mapFullToView(dydt);
@@ -159,7 +111,15 @@ namespace gridfire {
         return m_activeReactions;
     }
 
-    std::unordered_map<Species, double> DefinedEngineView::getSpeciesTimescales(
+    void DefinedEngineView::setNetworkReactions(const reaction::LogicalReactionSet &reactions) {
+        std::vector<std::string> peNames;
+        for (const auto& reaction : reactions) {
+            peNames.push_back(std::string(reaction.id()));
+        }
+        collect(peNames);
+    }
+
+    std::expected<std::unordered_map<Species, double>, expectations::StaleEngineError> DefinedEngineView::getSpeciesTimescales(
         const std::vector<double> &Y_defined,
         const double T9,
         const double rho
@@ -167,7 +127,11 @@ namespace gridfire {
         validateNetworkState();
 
         const auto Y_full = mapViewToFull(Y_defined);
-        const auto fullTimescales = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho);
+        const auto result = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho);
+        if (!result) {
+            return std::unexpected{result.error()};
+        }
+        const auto& fullTimescales = result.value();
 
         std::unordered_map<Species, double> definedTimescales;
         for (const auto& active_species : m_activeSpecies) {
@@ -178,12 +142,38 @@ namespace gridfire {
         return definedTimescales;
     }
 
+    std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError>
+    DefinedEngineView::getSpeciesDestructionTimescales(
+        const std::vector<double> &Y_defined,
+        const double T9,
+        const double rho
+    ) const {
+        validateNetworkState();
+
+        const auto Y_full = mapViewToFull(Y_defined);
+        const auto result = m_baseEngine.getSpeciesDestructionTimescales(Y_full, T9, rho);
+
+        if (!result) {
+            return std::unexpected{result.error()};
+        }
+
+        const auto& destructionTimescales = result.value();
+
+        std::unordered_map<Species, double> definedTimescales;
+        for (const auto& active_species : m_activeSpecies) {
+            if (destructionTimescales.contains(active_species)) {
+                definedTimescales[active_species] = destructionTimescales.at(active_species);
+            }
+        }
+        return definedTimescales;
+    }
+
     fourdst::composition::Composition DefinedEngineView::update(const NetIn &netIn) {
         return m_baseEngine.update(netIn);
     }
 
     bool DefinedEngineView::isStale(const NetIn &netIn) {
-        return false;
+        return m_baseEngine.isStale(netIn);
     }
 
 
@@ -198,7 +188,7 @@ namespace gridfire {
     int DefinedEngineView::getSpeciesIndex(const Species &species) const {
         validateNetworkState();
 
-        auto it = std::find(m_activeSpecies.begin(), m_activeSpecies.end(), species);
+        const auto it = std::ranges::find(m_activeSpecies, species);
         if (it != m_activeSpecies.end()) {
             return static_cast<int>(std::distance(m_activeSpecies.begin(), it));
         } else {
@@ -224,7 +214,7 @@ namespace gridfire {
     }
 
     std::vector<size_t> DefinedEngineView::constructSpeciesIndexMap() const {
-        LOG_TRACE_L1(m_logger, "Constructing species index map for DefinedEngineView...");
+        LOG_TRACE_L3(m_logger, "Constructing species index map for DefinedEngineView...");
         std::unordered_map<Species, size_t> fullSpeciesReverseMap;
         const auto& fullSpeciesList = m_baseEngine.getNetworkSpecies();
 
@@ -247,13 +237,13 @@ namespace gridfire {
                 throw std::runtime_error("Species not found in full species map: " + std::string(active_species.name()));
             }
         }
-        LOG_TRACE_L1(m_logger, "Species index map constructed with {} entries.", speciesIndexMap.size());
+        LOG_TRACE_L3(m_logger, "Species index map constructed with {} entries.", speciesIndexMap.size());
         return speciesIndexMap;
 
     }
 
     std::vector<size_t> DefinedEngineView::constructReactionIndexMap() const {
-        LOG_TRACE_L1(m_logger, "Constructing reaction index map for DefinedEngineView...");
+        LOG_TRACE_L3(m_logger, "Constructing reaction index map for DefinedEngineView...");
 
         // --- Step 1: Create a reverse map using the reaction's unique ID as the key. ---
         std::unordered_map<std::string_view, size_t> fullReactionReverseMap;
@@ -280,7 +270,7 @@ namespace gridfire {
             }
         }
 
-        LOG_TRACE_L1(m_logger, "Reaction index map constructed with {} entries.", reactionIndexMap.size());
+        LOG_TRACE_L3(m_logger, "Reaction index map constructed with {} entries.", reactionIndexMap.size());
         return reactionIndexMap;
     }
 
@@ -328,6 +318,59 @@ namespace gridfire {
         }
     }
 
+    void DefinedEngineView::collect(const std::vector<std::string> &peNames) {
+        std::unordered_set<Species> seenSpecies;
+
+        const auto& fullNetworkReactionSet = m_baseEngine.getNetworkReactions();
+        for (const auto& peName : peNames) {
+            if (!fullNetworkReactionSet.contains(peName)) {
+                LOG_ERROR(m_logger, "Reaction with name '{}' not found in the base engine's network reactions. Aborting...", peName);
+                m_logger->flush_log();
+                throw std::runtime_error("Reaction with name '" + std::string(peName) + "' not found in the base engine's network reactions.");
+            }
+            auto reaction = fullNetworkReactionSet[peName];
+            for (const auto& reactant : reaction.reactants()) {
+                if (!seenSpecies.contains(reactant)) {
+                    seenSpecies.insert(reactant);
+                    m_activeSpecies.push_back(reactant);
+                }
+            }
+            for (const auto& product : reaction.products()) {
+                if (!seenSpecies.contains(product)) {
+                    seenSpecies.insert(product);
+                    m_activeSpecies.push_back(product);
+                }
+            }
+            m_activeReactions.add_reaction(reaction);
+        }
+        LOG_TRACE_L3(m_logger, "DefinedEngineView built with {} active species and {} active reactions.", m_activeSpecies.size(), m_activeReactions.size());
+        LOG_TRACE_L3(m_logger, "Active species: {}", [this]() -> std::string {
+            std::string result;
+            for (const auto& species : m_activeSpecies) {
+                result += std::string(species.name()) + ", ";
+            }
+            if (!result.empty()) {
+                result.pop_back(); // Remove last space
+                result.pop_back(); // Remove last comma
+            }
+            return result;
+        }());
+        LOG_TRACE_L3(m_logger, "Active reactions: {}", [this]() -> std::string {
+            std::string result;
+            for (const auto& reaction : m_activeReactions) {
+                result += std::string(reaction.id()) + ", ";
+            }
+            if (!result.empty()) {
+                result.pop_back(); // Remove last space
+                result.pop_back(); // Remove last comma
+            }
+            return result;
+        }());
+        m_speciesIndexMap = constructSpeciesIndexMap();
+        m_reactionIndexMap = constructReactionIndexMap();
+        m_isStale = false;
+    }
+
 
     ////////////////////////////////////////////
     /// FileDefinedEngineView Implementation ///

src/network/lib/network.cpp

@@ -64,7 +64,7 @@ namespace gridfire {
     // Trim whitespace from both ends of a string
     std::string trim_whitespace(const std::string& str) {
         auto startIt = str.begin();
-        auto endIt   = str.end();
+        const auto endIt   = str.end();
 
         while (startIt != endIt && std::isspace(static_cast<unsigned char>(*startIt))) {
             ++startIt;
@@ -72,8 +72,8 @@ namespace gridfire {
         if (startIt == endIt) {
             return "";
         }
-        auto ritr = std::find_if(str.rbegin(), std::string::const_reverse_iterator(startIt),
+        const auto ritr = std::find_if(str.rbegin(), std::string::const_reverse_iterator(startIt),
-                                 [](unsigned char ch){ return !std::isspace(ch); });
+                                 [](const unsigned char ch){ return !std::isspace(ch); });
         return std::string(startIt, ritr.base());
     }
 

src/network/lib/partition/partition_rauscher_thielemann.cpp

@@ -22,17 +22,17 @@ namespace gridfire::partition {
         m_partitionData.reserve(numRecords);
         const auto* records = reinterpret_cast<const record::RauscherThielemannPartitionDataRecord*>(rauscher_thielemann_partition_data);
         for (size_t i = 0; i < numRecords; ++i) {
-            const auto& record = records[i];
+            const auto&[z, a, ground_state_spin, normalized_g_values] = records[i];
             IsotopeData data;
-            data.ground_state_spin = record.ground_state_spin;
+            data.ground_state_spin = ground_state_spin;
-            std::ranges::copy(record.normalized_g_values, data.normalized_g_values.begin());
+            std::ranges::copy(normalized_g_values, data.normalized_g_values.begin());
-            const int key = make_key(record.z, record.a);
+            const int key = make_key(z, a);
             LOG_TRACE_L3_LIMIT_EVERY_N(
                 100,
                 m_logger,
                 "(EVERY 100) Adding Rauscher-Thielemann partition data for Z={} A={} (key={})",
-                record.z,
+                z,
-                record.a,
+                a,
                 key
             );
 

src/network/lib/reaction/reaclib.cpp

@@ -12,7 +12,7 @@
 
 std::string trim_whitespace(const std::string& str) {
     auto startIt = str.begin();
-    auto endIt   = str.end();
+    const auto endIt   = str.end();
 
     while (startIt != endIt && std::isspace(static_cast<unsigned char>(*startIt))) {
         ++startIt;
@@ -20,8 +20,8 @@ std::string trim_whitespace(const std::string& str) {
     if (startIt == endIt) {
         return "";
     }
-    auto ritr = std::find_if(str.rbegin(), std::string::const_reverse_iterator(startIt),
+    const auto ritr = std::find_if(str.rbegin(), std::string::const_reverse_iterator(startIt),
-                             [](unsigned char ch){ return !std::isspace(ch); });
+                             [](const unsigned char ch){ return !std::isspace(ch); });
     return std::string(startIt, ritr.base());
 }
 
@@ -83,41 +83,41 @@ namespace gridfire::reaclib {
 
         // Cast the raw byte data to our structured record format.
         const auto* records = reinterpret_cast<const ReactionRecord*>(raw_reactions_data);
-        const size_t num_reactions = raw_reactions_data_len / sizeof(ReactionRecord);
+        constexpr size_t num_reactions = raw_reactions_data_len / sizeof(ReactionRecord);
 
         std::vector<reaction::Reaction> reaction_list;
         reaction_list.reserve(num_reactions);
 
         for (size_t i = 0; i < num_reactions; ++i) {
-            const auto& record = records[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(record.label, strnlen(record.label, sizeof(record.label)));
+            const std::string_view label_sv(label, strnlen(label, sizeof(label)));
-            const std::string_view rpName_sv(record.rpName, strnlen(record.rpName, sizeof(record.rpName)));
+            const std::string_view rpName_sv(rpName, strnlen(rpName, sizeof(rpName)));
-            const std::string_view reactants_sv(record.reactants_str, strnlen(record.reactants_str, sizeof(record.reactants_str)));
+            const std::string_view reactants_sv(reactants_str, strnlen(reactants_str, sizeof(reactants_str)));
-            const std::string_view products_sv(record.products_str, strnlen(record.products_str, sizeof(record.products_str)));
+            const std::string_view products_sv(products_str, strnlen(products_str, sizeof(products_str)));
 
             auto reactants = parseSpeciesString(reactants_sv);
             auto products = parseSpeciesString(products_sv);
 
             const reaction::RateCoefficientSet rate_coeffs = {
-                record.coeffs[0], record.coeffs[1], record.coeffs[2],
+                coeffs[0], coeffs[1], coeffs[2],
-                record.coeffs[3], record.coeffs[4], record.coeffs[5],
+                coeffs[3], coeffs[4], coeffs[5],
-                record.coeffs[6]
+                coeffs[6]
             };
 
             // Construct the Reaction object. We use rpName for both the unique ID and the human-readable name.
             reaction_list.emplace_back(
                 rpName_sv,
                 rpName_sv,
-                record.chapter,
+                chapter,
                 reactants,
                 products,
-                record.qValue,
+                qValue,
                 label_sv,
                 rate_coeffs,
-                record.reverse
+                reverse
             );
         }
 

src/network/lib/reaction/reaction.cpp

@@ -154,7 +154,7 @@ namespace gridfire::reaction {
         return (reactantMass - productMass) * AMU2MeV;
     }
 
-    uint64_t Reaction::hash(uint64_t seed) const {
+    uint64_t Reaction::hash(const uint64_t seed) const {
         return XXHash64::hash(m_id.data(), m_id.size(), seed);
     }
 

src/network/lib/solver/solver.cpp

@@ -3,518 +3,22 @@
 #include "gridfire/network.h"
 #include "gridfire/exceptions/error_engine.h"
 
-
-#include "gridfire/utils/logging.h"
-
 #include "fourdst/composition/atomicSpecies.h"
 #include "fourdst/composition/composition.h"
 #include "fourdst/config/config.h"
 
-#include "Eigen/Dense"
 #include "unsupported/Eigen/NonLinearOptimization"
 
 #include <boost/numeric/odeint.hpp>
 
 #include <vector>
-#include <unordered_map>
 #include <string>
 #include <stdexcept>
 #include <iomanip>
 
-#include "gridfire/engine/views/engine_multiscale.h"
-
 #include "quill/LogMacros.h"
 
-struct adaptive_integrator_observer {
-    double m_last_t = 0.0;
-    int m_total_steps = 0;
-
-    void operator()(boost::numeric::ublas::vector<double> &state, double t) {
-        // std::cout << "(Step: " << m_total_steps << ") Current Time: " << t << " (dt: " << t - m_last_t << ")\n";
-        m_last_t = t;
-        m_total_steps++;
-
-    }
-};
-
 namespace gridfire::solver {
-    double s_prevTimestep = 0.0;
-
-    NetOut QSENetworkSolver::evaluate(const NetIn &netIn) {
-        // --- Use the policy to decide whether to update the view ---
-        if (shouldUpdateView(netIn)) {
-            LOG_DEBUG(m_logger, "Solver update policy triggered, network view updating...");
-            m_engine.update(netIn);
-            LOG_DEBUG(m_logger, "Network view updated!");
-
-            m_lastSeenConditions = netIn;
-            m_isViewInitialized = true;
-        }
-        m_engine.generateJacobianMatrix(netIn.MolarAbundance(), netIn.temperature / 1e9, netIn.density);
-        using state_type = boost::numeric::ublas::vector<double>;
-        using namespace boost::numeric::odeint;
-
-        NetOut postIgnition = initializeNetworkWithShortIgnition(netIn);
-
-        constexpr double abundance_floor = 1.0e-30;
-        std::vector<double>Y_sanitized_initial;
-        Y_sanitized_initial.reserve(m_engine.getNetworkSpecies().size());
-
-        LOG_DEBUG(m_logger, "Sanitizing initial abundances with a floor of {:0.3E}...", abundance_floor);
-        for (const auto& species : m_engine.getNetworkSpecies()) {
-            double molar_abundance = 0.0;
-            if (postIgnition.composition.contains(species)) {
-                molar_abundance = postIgnition.composition.getMolarAbundance(std::string(species.name()));
-            }
-
-            if (molar_abundance < abundance_floor) {
-                molar_abundance = abundance_floor;
-            }
-            Y_sanitized_initial.push_back(molar_abundance);
-        }
-        const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9)
-        const double rho = netIn.density; // Density in g/cm^3
-
-        const auto indices = packSpeciesTypeIndexVectors(Y_sanitized_initial, T9, rho);
-        Eigen::VectorXd Y_QSE;
-        try {
-            Y_QSE = calculateSteadyStateAbundances(Y_sanitized_initial, T9, rho, indices);
-            LOG_DEBUG(m_logger, "QSE Abundances: {}", [*this](const dynamicQSESpeciesIndices& indices, const Eigen::VectorXd& Y_QSE) -> std::string {
-                std::stringstream ss;
-                ss << std::scientific << std::setprecision(5);
-                for (size_t i = 0; i < indices.QSESpeciesIndices.size(); ++i) {
-                    ss << std::string(m_engine.getNetworkSpecies()[indices.QSESpeciesIndices[i]].name()) + ": ";
-                    ss << Y_QSE(i);
-                    if (i < indices.QSESpeciesIndices.size() - 2) {
-                        ss << ", ";
-                    } else if (i == indices.QSESpeciesIndices.size() - 2) {
-                        ss << ", and ";
-                    }
-
-                }
-                return ss.str();
-            }(indices, Y_QSE));
-        } catch (const std::runtime_error& e) {
-            LOG_ERROR(m_logger, "Failed to calculate steady state abundances. Aborting QSE evaluation.");
-            m_logger->flush_log();
-            throw std::runtime_error("Failed to calculate steady state abundances: " + std::string(e.what()));
-        }
-
-        state_type YDynamic_ublas(indices.dynamicSpeciesIndices.size() + 1);
-        for (size_t i = 0; i < indices.dynamicSpeciesIndices.size(); ++i) {
-            YDynamic_ublas(i) = Y_sanitized_initial[indices.dynamicSpeciesIndices[i]];
-        }
-        YDynamic_ublas(indices.dynamicSpeciesIndices.size()) = 0.0; // Placeholder for specific energy rate
-
-        const RHSFunctor rhs_functor(m_engine, indices.dynamicSpeciesIndices, indices.QSESpeciesIndices, Y_QSE, T9, rho);
-        const auto stepper = make_controlled<runge_kutta_dopri5<state_type>>(1.0e-8, 1.0e-8);
-
-        size_t stepCount = integrate_adaptive(
-            stepper,
-            rhs_functor,
-            YDynamic_ublas,
-            0.0, // Start time
-            netIn.tMax,
-            netIn.dt0
-        );
-
-        std::vector<double> YFinal = Y_sanitized_initial;
-        for (size_t i = 0; i <indices.dynamicSpeciesIndices.size(); ++i) {
-            YFinal[indices.dynamicSpeciesIndices[i]] = YDynamic_ublas(i);
-        }
-        for (size_t i = 0; i < indices.QSESpeciesIndices.size(); ++i) {
-            YFinal[indices.QSESpeciesIndices[i]] = Y_QSE(i);
-        }
-
-        const double finalSpecificEnergyRate = YDynamic_ublas(indices.dynamicSpeciesIndices.size());
-
-        // --- Marshal output variables ---
-        std::vector<std::string> speciesNames(m_engine.getNetworkSpecies().size());
-        std::vector<double> finalMassFractions(m_engine.getNetworkSpecies().size());
-        double massFractionSum = 0.0;
-
-        for (size_t i = 0; i < speciesNames.size(); ++i) {
-            const auto& species = m_engine.getNetworkSpecies()[i];
-            speciesNames[i] = species.name();
-            finalMassFractions[i] = YFinal[i] * species.mass(); // Convert from molar abundance to mass fraction
-            massFractionSum += finalMassFractions[i];
-        }
-        for (auto& mf : finalMassFractions) {
-            mf /= massFractionSum; // Normalize to get mass fractions
-        }
-
-        fourdst::composition::Composition outputComposition(speciesNames, finalMassFractions);
-        NetOut netOut;
-        netOut.composition = outputComposition;
-        netOut.energy = finalSpecificEnergyRate; // Specific energy rate
-        netOut.num_steps = stepCount;
-        return netOut;
-    }
-
-    dynamicQSESpeciesIndices QSENetworkSolver::packSpeciesTypeIndexVectors(
-        const std::vector<double>& Y,
-        const double T9,
-        const double rho
-    ) const {
-        constexpr double timescaleCutoff = 1.0e-5;
-        constexpr double abundanceCutoff = 1.0e-15;
-
-        LOG_INFO(m_logger, "Partitioning species using T9={:0.2f} and ρ={:0.2e}", T9, rho);
-        LOG_INFO(m_logger, "Timescale Cutoff: {:.1e} s, Abundance Cutoff: {:.1e}", timescaleCutoff, abundanceCutoff);
-
-        std::vector<size_t>dynamicSpeciesIndices; // Slow species that are not in QSE
-        std::vector<size_t>QSESpeciesIndices;  // Fast species that are in QSE
-
-        std::unordered_map<fourdst::atomic::Species, double> speciesTimescale = m_engine.getSpeciesTimescales(Y, T9, rho);
-        LOG_DEBUG(
-            m_logger,
-            "{}",
-            gridfire::utils::formatNuclearTimescaleLogString(
-                m_engine,
-                Y,
-                T9,
-                rho
-            ));
-
-        for (size_t i = 0; i < m_engine.getNetworkSpecies().size(); ++i) {
-            const auto& species = m_engine.getNetworkSpecies()[i];
-            const double network_timescale = speciesTimescale.at(species);
-            const double abundance = Y[i];
-
-            double decay_timescale = std::numeric_limits<double>::infinity();
-            const double half_life = species.halfLife();
-            if (half_life > 0 && !std::isinf(half_life)) {
-                constexpr double LN2 = 0.69314718056;
-                decay_timescale = half_life / LN2;
-            }
-
-            const double final_timescale = std::min(network_timescale, decay_timescale);
-
-            const bool isQSE = false;
-            //     network_timescale,
-            //     decay_timescale,
-            //     abundance
-            // );
-
-            if (isQSE) {
-                LOG_TRACE_L2(m_logger, "{} is in QSE based on rules in qse_rules.h", species.name());
-                QSESpeciesIndices.push_back(i);
-            } else {
-                LOG_TRACE_L2(m_logger, "{} is dynamic based on rules in qse_rules.h", species.name());
-                dynamicSpeciesIndices.push_back(i);
-            }
-        }
-        LOG_INFO(m_logger, "Partitioning complete. Dynamical species: {}, QSE species: {}.", dynamicSpeciesIndices.size(), QSESpeciesIndices.size());
-        LOG_INFO(m_logger, "Dynamic species: {}", [dynamicSpeciesIndices](const DynamicEngine& engine_wrapper) -> std::string {
-            std::string result;
-            int count = 0;
-            for (const auto& i : dynamicSpeciesIndices) {
-                result += std::string(engine_wrapper.getNetworkSpecies()[i].name());
-                if (count < dynamicSpeciesIndices.size() - 2) {
-                    result += ", ";
-                } else if (count == dynamicSpeciesIndices.size() - 2) {
-                    result += " and ";
-                }
-                count++;
-            }
-            return result;
-        }(m_engine));
-        LOG_INFO(m_logger, "QSE species: {}", [QSESpeciesIndices](const DynamicEngine& engine_wrapper) -> std::string {
-            std::string result;
-            int count = 0;
-            for (const auto& i : QSESpeciesIndices) {
-                result += std::string(engine_wrapper.getNetworkSpecies()[i].name());
-                if (count < QSESpeciesIndices.size() - 2) {
-                    result += ", ";
-                } else if (count == QSESpeciesIndices.size() - 2) {
-                    result += " and ";
-                }
-                count++;
-            }
-            return result;
-        }(m_engine));
-        return {dynamicSpeciesIndices, QSESpeciesIndices};
-    }
-
-    Eigen::VectorXd QSENetworkSolver::calculateSteadyStateAbundances(
-    const std::vector<double> &Y,
-    const double T9,
-    const double rho,
-    const dynamicQSESpeciesIndices &indices
-    ) const {
-        LOG_TRACE_L1(m_logger, "Calculating steady state abundances for QSE species...");
-        LOG_TRACE_L1(
-            m_logger,
-            "Initial QSE species abundances: {}",
-            [&]() -> std::string {
-                std::stringstream ss;
-                ss << std::scientific << std::setprecision(5);
-                int count = 0;
-                for (const auto& i : indices.QSESpeciesIndices) {
-                    ss << std::string(m_engine.getNetworkSpecies()[i].name()) << ": " << Y[i] << ", ";
-                    if (count < indices.QSESpeciesIndices.size() - 2) {
-                        ss << ", ";
-                    } else if (count == indices.QSESpeciesIndices.size() - 2) {
-                        ss << " and ";
-                    }
-                    count++;
-                }
-                return ss.str();
-            }());
-
-        if (indices.QSESpeciesIndices.empty()) {
-            LOG_DEBUG(m_logger, "No QSE species to solve for.");
-            return Eigen::VectorXd(0);
-        }
-
-        Eigen::VectorXd YScale(indices.QSESpeciesIndices.size());
-        const double abundanceFloor = 1e-15;
-        for (size_t i = 0; i < indices.QSESpeciesIndices.size(); i++) {
-            const double initial_abundance = Y[indices.QSESpeciesIndices[i]];
-            YScale(i) = std::max(initial_abundance, abundanceFloor);
-        }
-
-        // Use the EigenFunctor with Eigen's nonlinear solver
-        EigenFunctor<double> functor(
-            m_engine,
-            Y,
-            indices.dynamicSpeciesIndices,
-            indices.QSESpeciesIndices,
-            T9,
-            rho,
-            YScale
-        );
-
-        Eigen::VectorXd v_qse_initial(indices.QSESpeciesIndices.size());
-        for (size_t i = 0; i < indices.QSESpeciesIndices.size(); ++i) {
-            const double initial_abundance = Y[indices.QSESpeciesIndices[i]];
-            v_qse_initial(i) = std::asinh(initial_abundance/ YScale(i)); // Use asinh for better numerical stability compared to log
-        }
-        LOG_TRACE_L1(
-            m_logger,
-            "v_qse_log_initial: {}",
-            [&]() -> std::string {
-                std::stringstream ss;
-                ss << std::scientific << std::setprecision(5);
-                for (size_t i = 0; i < v_qse_initial.size(); ++i) {
-                    ss << "log(X_" << std::string(m_engine.getNetworkSpecies()[indices.QSESpeciesIndices[i]].name()) << ") = " << v_qse_initial(i);
-                    if (i < v_qse_initial.size() - 2) {
-                        ss << ", ";
-                    } else if (i == v_qse_initial.size() - 2) {
-                        ss << " and ";
-                    }
-                }
-                return ss.str();
-            }());
-
-        const static std::unordered_map<Eigen::LevenbergMarquardtSpace::Status, const char*> statusMessages = {
-            {Eigen::LevenbergMarquardtSpace::NotStarted, "Not started"},
-            {Eigen::LevenbergMarquardtSpace::Running, "Running"},
-            {Eigen::LevenbergMarquardtSpace::ImproperInputParameters, "Improper input parameters"},
-            {Eigen::LevenbergMarquardtSpace::RelativeReductionTooSmall, "Relative reduction too small"},
-            {Eigen::LevenbergMarquardtSpace::RelativeErrorTooSmall, "Relative error too small"},
-            {Eigen::LevenbergMarquardtSpace::RelativeErrorAndReductionTooSmall, "Relative error and reduction too small"},
-            {Eigen::LevenbergMarquardtSpace::CosinusTooSmall, "Cosine too small"},
-            {Eigen::LevenbergMarquardtSpace::TooManyFunctionEvaluation, "Too many function evaluations"},
-            {Eigen::LevenbergMarquardtSpace::FtolTooSmall, "Function tolerance too small"},
-            {Eigen::LevenbergMarquardtSpace::XtolTooSmall, "X tolerance too small"},
-            {Eigen::LevenbergMarquardtSpace::GtolTooSmall, "Gradient tolerance too small"},
-            {Eigen::LevenbergMarquardtSpace::UserAsked, "User asked to stop"}
-        };
-
-        Eigen::LevenbergMarquardt lm(functor);
-        lm.parameters.xtol = 1e-24;
-        lm.parameters.ftol = 1e-24;
-        const Eigen::LevenbergMarquardtSpace::Status info = lm.minimize(v_qse_initial);
-
-        if (info <= 0 || info >= 4) {
-            LOG_ERROR(m_logger, "QSE species minimization failed with status: {} ({})",
-                      static_cast<int>(info), statusMessages.at(info));
-            throw std::runtime_error(
-                "QSE species minimization failed with status: " + std::to_string(static_cast<int>(info)) +
-                " (" + std::string(statusMessages.at(info)) + ")"
-            );
-        }
-        LOG_DEBUG(m_logger, "QSE species minimization completed successfully with status: {} ({})",
-                  static_cast<int>(info), statusMessages.at(info));
-        return YScale.array() * v_qse_initial.array().sinh(); // Convert back to molar abundances
-
-    }
-
-    NetOut QSENetworkSolver::initializeNetworkWithShortIgnition(const NetIn &netIn) const {
-        const auto ignitionTemperature = m_config.get<double>(
-            "gridfire:solver:QSE:ignition:temperature",
-            2e8
-        ); // 0.2 GK
-        const auto ignitionDensity = m_config.get<double>(
-            "gridfire:solver:QSE:ignition:density",
-            1e6
-        ); // 1e6 g/cm^3
-        const auto ignitionTime = m_config.get<double>(
-            "gridfire:solver:QSE:ignition:tMax",
-            1e-7
-        ); // 0.1 μs
-        const auto ignitionStepSize = m_config.get<double>(
-            "gridfire:solver:QSE:ignition:dt0",
-            1e-15
-        ); // 1e-15 seconds
-
-        LOG_INFO(
-            m_logger,
-            "Igniting network with T={:<5.3E}, ρ={:<5.3E}, tMax={:<5.3E}, dt0={:<5.3E}...",
-            ignitionTemperature,
-            ignitionDensity,
-            ignitionTime,
-            ignitionStepSize
-        );
-        LOG_INFO(
-            m_logger,
-            "Pre-ignition composition: {}",
-            [&]() -> std::string {
-                std::stringstream ss;
-                ss << std::scientific << std::setprecision(3);
-                int count = 0;
-                for (const auto& entry : netIn.composition | std::views::values) {
-                    ss << "X_" << std::string(entry.isotope().name()) << ": " << entry.mass_fraction() << " ";
-                    if (count < netIn.composition.getRegisteredSymbols().size() - 2) {
-                        ss << ", ";
-                    } else if (count == netIn.composition.getRegisteredSymbols().size() - 2) {
-                        ss << " and ";
-                    }
-                    count++;
-                }
-                return ss.str();
-            }());
-
-        NetIn preIgnition = netIn;
-        preIgnition.temperature = ignitionTemperature;
-        preIgnition.density = ignitionDensity;
-        preIgnition.tMax = ignitionTime;
-        preIgnition.dt0 = ignitionStepSize;
-
-        const auto prevScreeningModel = m_engine.getScreeningModel();
-        LOG_DEBUG(m_logger, "Setting screening model to BARE for high temperature and density ignition.");
-        m_engine.setScreeningModel(screening::ScreeningType::BARE);
-        DirectNetworkSolver ignitionSolver(m_engine);
-        NetOut postIgnition = ignitionSolver.evaluate(preIgnition);
-        LOG_INFO(m_logger, "Network ignition completed in {} steps.", postIgnition.num_steps);
-        LOG_INFO(
-            m_logger,
-            "Post-ignition composition: {}",
-            [&]() -> std::string {
-                std::stringstream ss;
-                ss << std::scientific << std::setprecision(3);
-                int count = 0;
-                for (const auto& entry : postIgnition.composition | std::views::values) {
-                    ss << "X_" << std::string(entry.isotope().name()) << ": " << entry.mass_fraction() << " ";
-                    if (count < postIgnition.composition.getRegisteredSymbols().size() - 2) {
-                        ss << ", ";
-                    } else if (count == postIgnition.composition.getRegisteredSymbols().size() - 2) {
-                        ss << " and ";
-                    }
-                    count++;
-                }
-                return ss.str();
-            }());
-        LOG_DEBUG(
-            m_logger,
-            "Average change in composition during ignition: {}",
-            [&]() -> std::string {
-                std::stringstream ss;
-                ss << std::scientific << std::setprecision(3);
-                for (const auto& entry : postIgnition.composition | std::views::values) {
-                    if (!postIgnition.composition.contains(entry.isotope()) || !netIn.composition.contains(entry.isotope())) {
-                        ss << std::string(entry.isotope().name()) << ": inf %, ";
-                        continue;
-                    }
-                    const double initialAbundance = netIn.composition.getMolarAbundance(std::string(entry.isotope().name()));
-                    const double finalAbundance = postIgnition.composition.getMolarAbundance(std::string(entry.isotope().name()));
-                    const double change = (finalAbundance - initialAbundance) / initialAbundance * 100.0; // Percentage change
-                    ss << std::string(entry.isotope().name()) << ": " << change << " %, ";
-                }
-                return ss.str();
-            }());
-        m_engine.setScreeningModel(prevScreeningModel);
-        LOG_DEBUG(m_logger, "Restoring previous screening model: {}", static_cast<int>(prevScreeningModel));
-        return postIgnition;
-    }
-
-    bool QSENetworkSolver::shouldUpdateView(const NetIn &conditions) const {
-        // Policy 1: If the view has never been initialized, we must update.
-        if (!m_isViewInitialized) {
-            return true;
-        }
-
-        // Policy 2: Check for significant relative change in temperature.
-        // Reaction rates are exponentially sensitive to temperature, so we use a tight threshold.
-        const double temp_threshold = m_config.get<double>("gridfire:solver:policy:temp_threshold", 0.05); // 5%
-        const double temp_relative_change = std::abs(conditions.temperature - m_lastSeenConditions.temperature) / m_lastSeenConditions.temperature;
-        if (temp_relative_change > temp_threshold) {
-            LOG_DEBUG(m_logger, "Temperature changed by {:.1f}%, triggering view update.", temp_relative_change * 100);
-            return true;
-        }
-
-        // Policy 3: Check for significant relative change in density.
-        const double rho_threshold = m_config.get<double>("gridfire:solver:policy:rho_threshold", 0.10); // 10%
-        const double rho_relative_change = std::abs(conditions.density - m_lastSeenConditions.density) / m_lastSeenConditions.density;
-        if (rho_relative_change > rho_threshold) {
-            LOG_DEBUG(m_logger, "Density changed by {:.1f}%, triggering view update.", rho_relative_change * 100);
-            return true;
-        }
-
-        // Policy 4: Check for fuel depletion.
-        // If a primary fuel source changes significantly, the network structure may change.
-        const double fuel_threshold = m_config.get<double>("gridfire:solver:policy:fuel_threshold", 0.15); // 15%
-
-        // Example: Check hydrogen abundance
-        const double h1_old = m_lastSeenConditions.composition.getMassFraction("H-1");
-        const double h1_new = conditions.composition.getMassFraction("H-1");
-        if (h1_old > 1e-12) { // Avoid division by zero
-            const double h1_relative_change = std::abs(h1_new - h1_old) / h1_old;
-            if (h1_relative_change > fuel_threshold) {
-                LOG_DEBUG(m_logger, "H-1 mass fraction changed by {:.1f}%, triggering view update.", h1_relative_change * 100);
-                return true;
-            }
-        }
-
-        // If none of the above conditions are met, the current view is still good enough.
-        return false;
-    }
-
-    void QSENetworkSolver::RHSFunctor::operator()(
-        const boost::numeric::ublas::vector<double> &YDynamic,
-        boost::numeric::ublas::vector<double> &dYdtDynamic,
-        double t
-    ) const {
-        LOG_TRACE_L1(
-            m_logger,
-            "Timestepping at t={:0.3E} (dt={:0.3E}) with T9={:0.3E}, ρ={:0.3E}",
-            t,
-            t - s_prevTimestep,
-            m_T9,
-            m_rho
-        );
-        s_prevTimestep = t;
-        // --- Populate the slow / dynamic species vector ---
-        std::vector<double> YFull(m_engine.getNetworkSpecies().size());
-        for (size_t i = 0; i < m_dynamicSpeciesIndices.size(); ++i) {
-            YFull[m_dynamicSpeciesIndices[i]] = YDynamic(i);
-        }
-
-        // --- Populate the QSE species vector ---
-        for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
-            YFull[m_QSESpeciesIndices[i]] = m_Y_QSE(i);
-        }
-
-        auto [full_dYdt, specificEnergyRate] = m_engine.calculateRHSAndEnergy(YFull, m_T9, m_rho);
-
-        dYdtDynamic.resize(m_dynamicSpeciesIndices.size() + 1);
-        for (size_t i = 0; i < m_dynamicSpeciesIndices.size(); ++i) {
-            dYdtDynamic(i) = full_dYdt[m_dynamicSpeciesIndices[i]];
-        }
-        dYdtDynamic[m_dynamicSpeciesIndices.size()] = specificEnergyRate;
-    }
-
     NetOut DirectNetworkSolver::evaluate(const NetIn &netIn) {
         namespace ublas = boost::numeric::ublas;
         namespace odeint = boost::numeric::odeint;
@@ -526,44 +30,97 @@ namespace gridfire::solver {
         const auto absTol = m_config.get<double>("gridfire:solver:DirectNetworkSolver:absTol", 1.0e-8);
         const auto relTol = m_config.get<double>("gridfire:solver:DirectNetworkSolver:relTol", 1.0e-8);
 
-        size_t stepCount = 0;
-
         Composition equilibratedComposition = m_engine.update(netIn);
-        const unsigned long numSpecies = m_engine.getNetworkSpecies().size();
+        size_t numSpecies = m_engine.getNetworkSpecies().size();
         ublas::vector<double> Y(numSpecies + 1);
 
-        RHSFunctor rhsFunctor(m_engine, T9, netIn.density);
+        RHSManager manager(m_engine, T9, netIn.density);
         JacobianFunctor jacobianFunctor(m_engine, T9, netIn.density);
 
+        auto populateY = [&](const Composition& comp) {
+            const size_t numSpeciesInternal = m_engine.getNetworkSpecies().size();
+            Y.resize(numSpeciesInternal + 1);
+            for (size_t i = 0; i < numSpeciesInternal; i++) {
+                const auto& species = m_engine.getNetworkSpecies()[i];
+                if (!comp.contains(species)) {
+                    double lim = std::numeric_limits<double>::min();
+                    LOG_DEBUG(m_logger, "Species '{}' not found in composition. Setting abundance to {:0.3E}.", species.name(), lim);
+                    Y(i) = lim; // Species not in the composition, set to zero
+                } else {
+                    Y(i) = comp.getMolarAbundance(species);
+                }
+            }
+
+            // TODO: a good starting point to make the temperature, density, and energy self consistent would be to turn this into an accumulator
+            Y(numSpeciesInternal) = 0.0; // Specific energy rate, initialized to zero
+        };
 
         // This is a quick debug that can be turned on. For solar code input parameters (T~1.5e7K, ρ~1.5e3 g/cm^3) this should be near 8e-17
         // std::cout << "D/H: " << equilibratedComposition.getMolarAbundance("H-2") / equilibratedComposition.getMolarAbundance("H-1") << std::endl;
 
-        for (size_t i = 0; i < numSpecies; ++i) {
+        populateY(equilibratedComposition);
-            const auto& species = m_engine.getNetworkSpecies()[i];
-            if (!equilibratedComposition.contains(species)) {
-                LOG_DEBUG(m_logger, "Species '{}' not found in composition. Setting abundance to 1e-99.", species.name());
-                Y(i) = 1e-99;
-                continue;
-            }
-            double abun = equilibratedComposition.getMolarAbundance(species);
-            // if (abun == 0.0) {
-            //     abun = 1e-99; // Avoid zero abundances
-            // }
-            Y(i) = abun;
-        }
-        Y(numSpecies) = 0.0;
         const auto stepper = odeint::make_controlled<odeint::rosenbrock4<double>>(absTol, relTol);
-        adaptive_integrator_observer observer;
+
-        stepCount = odeint::integrate_adaptive(
+        double current_time = 0.0;
+        double current_initial_timestep = netIn.dt0;
+        double accumulated_energy = 0.0;
+        // size_t total_update_stages_triggered = 0;
+
+        while (current_time < netIn.tMax) {
+            try {
+                odeint::integrate_adaptive(
                     stepper,
-            std::make_pair(rhsFunctor, jacobianFunctor),
+                    std::make_pair(manager, jacobianFunctor),
                     Y,
-            0.0,
+                    current_time,
                     netIn.tMax,
-            netIn.dt0,
+                    current_initial_timestep,
-            observer
+                    [&](const auto& state, double t) {
+                        current_time = t;
+                        manager.observe(state, t);
+                    }
                 );
+                current_time = netIn.tMax;
+            } catch (const exceptions::StaleEngineTrigger &e) {
+                LOG_INFO(m_logger, "Catching StaleEngineTrigger at t = {:0.3E} with T9 = {:0.3E}, rho = {:0.3E}. Triggering update stage (last stage took {} steps).", current_time, T9, netIn.density, e.totalSteps());
+                exceptions::StaleEngineTrigger::state staleState = e.getState();
+                accumulated_energy += e.energy(); // Add the specific energy rate to the accumulated energy
+                // total_update_stages_triggered++;
+
+                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) {
+                    const auto& species = m_engine.getNetworkSpecies()[i];
+                    temp_comp.registerSpecies(species);
+                    mass_fractions.push_back(e.getMolarAbundance(i) * species.mass()); // Convert from molar abundance to mass fraction
+                }
+                temp_comp.setMassFraction(m_engine.getNetworkSpecies(), mass_fractions);
+                temp_comp.finalize(true);
+
+                NetIn netInTemp = netIn;
+                netInTemp.temperature = e.temperature();
+                netInTemp.density = e.density();
+                netInTemp.composition = std::move(temp_comp);
+
+                Composition currentComposition = m_engine.update(netInTemp);
+                populateY(currentComposition);
+                Y(Y.size() - 1) = e.energy(); // Set the specific energy rate from the stale state
+                numSpecies = m_engine.getNetworkSpecies().size();
+
+                // current_initial_timestep = 0.001 * manager.m_last_step_time; // set the new timestep to the last successful timestep before repartitioning
+            }
+        }
+
+        accumulated_energy += Y(Y.size() - 1); // Add the specific energy rate to the accumulated energy
 
         std::vector<double> finalMassFractions(numSpecies);
         for (size_t i = 0; i < numSpecies; ++i) {
@@ -586,72 +143,73 @@ namespace gridfire::solver {
 
         NetOut netOut;
         netOut.composition = std::move(outputComposition);
-        netOut.energy = Y(numSpecies); // Specific energy rate
+        netOut.energy = accumulated_energy; // Specific energy rate
-        netOut.num_steps = stepCount;
+        netOut.num_steps = manager.m_num_steps;
 
         return netOut;
     }
 
-    void DirectNetworkSolver::RHSFunctor::operator()(
+    void DirectNetworkSolver::RHSManager::operator()(
         const boost::numeric::ublas::vector<double> &Y,
         boost::numeric::ublas::vector<double> &dYdt,
         const double t
     ) const {
-        const std::vector<double> y(Y.begin(), m_numSpecies + Y.begin());
+        const size_t numSpecies = m_engine.getNetworkSpecies().size();
-        // const auto timescales = m_engine.getSpeciesTimescales(y, m_T9, m_rho);
+        if (t != m_cached_time || !m_cached_result.has_value() || m_cached_result.value().dydt.size() != numSpecies + 1) {
-        StepDerivatives<double> deriv;
+            compute_and_cache(Y, t);
-
+        }
-        // TODO: Refactor this to use std::expected
+        const auto&[dydt, nuclearEnergyGenerationRate] = m_cached_result.value();
-        try {
+        dYdt.resize(numSpecies + 1);
-            deriv = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho);
+        std::ranges::copy(dydt, dYdt.begin());
-        } catch (const exceptions::StaleEngineError& e) { // If the engine reports that it is stale, we need to update it
+        dYdt(numSpecies) = nuclearEnergyGenerationRate; // Set the last element to the specific energy rate
-            LOG_INFO(m_logger, "Engine reports stale state. Calling engine update method...");
-
-            NetIn netInTemp;
-            fourdst::composition::Composition compositionTemp;
-
-            for (const auto& species : m_engine.getNetworkSpecies()) {
-                compositionTemp.registerSpecies(species);
     }
 
-            std::vector<double> X(y.size(), 0.0);
+    void DirectNetworkSolver::RHSManager::observe(
-            for (size_t i = 0; i < y.size(); ++i) {
+        const boost::numeric::ublas::vector<double> &state,
-                double Xi = y[i] * m_engine.getNetworkSpecies()[i].mass(); // Convert from molar abundance to mass fraction
+        const double t
-                if (Xi < 0 && std::abs(Xi) < 1e-20) {
+    ) const {
-                    Xi = 0.0; // Avoid negative mass fractions
+        double dt = t - m_last_observed_time;
-                }
+        compute_and_cache(state, t);
-                X[i] = Xi;
+        LOG_INFO(
+            m_logger,
+            "(Step {}) Observed state at t = {:0.3E} (dt = {:0.3E})",
+            m_num_steps,
+            t,
+            dt
+        );
+        std::ostringstream oss;
+        oss << std::scientific << std::setprecision(3);
+        oss << "(Step: " << std::setw(10) << m_num_steps << ") t = " << t << " (dt = " << dt << ", eps_nuc: " << state(state.size() - 1) << " [erg])\n";
+        std::cout << oss.str();
+        m_last_observed_time = t;
+        m_last_step_time = dt;
+
     }
 
-            compositionTemp.setMassFraction(m_engine.getNetworkSpecies(), X);
+    void DirectNetworkSolver::RHSManager::compute_and_cache(
-            compositionTemp.finalize(true);
+        const boost::numeric::ublas::vector<double> &state,
+        double t
+    ) const {
+        std::vector<double> y_vec(state.begin(), state.end() - 1);
+        std::ranges::replace_if(
+            y_vec,
+            [](const double yi){
+                return yi < 0.0;
+            },
+            0.0 // Avoid negative abundances
+        );
 
-            netInTemp.composition = std::move(compositionTemp);
+        const auto result = m_engine.calculateRHSAndEnergy(y_vec, m_T9, m_rho);
-            netInTemp.temperature = m_T9 * 1e9; // Convert T9 back to Kelvin
+        if (!result) {
-            netInTemp.density = m_rho; // Density in g/cm^3
+            LOG_INFO(m_logger,
-
+                "Triggering update stage due to stale engine detected at t = {:0.3E} with T9 = {:0.3E}, rho = {:0.3E}, y_vec (size: {})",
-            m_engine.update(netInTemp);
+                t, m_T9, m_rho, y_vec.size());
-
+            throw exceptions::StaleEngineTrigger({m_T9, m_rho, y_vec, t, m_num_steps, state(state.size() - 1)});
-            LOG_INFO(m_logger, "Engine update complete. Recalculating RHS and energy...");
-            deriv = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho);
-            LOG_INFO(m_logger, "RHS and energy recalculated successfully!");
         }
+        m_cached_result = result.value();
+        m_cached_time = t;
 
-        dYdt.resize(m_numSpecies + 1);
+        m_num_steps++;
-        // for (size_t i = 0; i < m_numSpecies; ++i) {
-        //     std::cout << "Species: " << m_engine.getNetworkSpecies()[i].name() << ", dYdt: " << deriv.dydt[i] << '\n';
-        // }
-
-        // std::vector<double> S;
-        // S.reserve(m_numSpecies);
-        // for (size_t i = 0; i < m_numSpecies; ++i) {
-        //     S.push_back(std::abs(deriv.dydt[i]) / (1e-8 + 1e-8 * y[i]));
-        // }
-        // for (size_t i = 0; i < m_numSpecies; ++i) {
-        //     std::cout << "Species: " << m_engine.getNetworkSpecies()[i].name() << ", S: " << S[i] << '\n';
-        // }
-        std::ranges::copy(deriv.dydt, dYdt.begin());
-        dYdt(m_numSpecies) = deriv.nuclearEnergyGenerationRate;
     }
 
     void DirectNetworkSolver::JacobianFunctor::operator()(
@@ -660,13 +218,14 @@ namespace gridfire::solver {
         double t,
         boost::numeric::ublas::vector<double> &dfdt
     ) const {
-        J.resize(m_numSpecies+1, m_numSpecies+1);
+        size_t numSpecies = m_engine.getNetworkSpecies().size();
+        J.resize(numSpecies+1, numSpecies+1);
         J.clear();
-        for (int i = 0; i < m_numSpecies; ++i) {
+        for (int i = 0; i < numSpecies; ++i) {
-            for (int j = 0; j < m_numSpecies; ++j) {
+            for (int j = 0; j < numSpecies; ++j) {
                 J(i, j) = m_engine.getJacobianMatrixEntry(i, j);
-                // std::cout << "J(" << m_engine.getNetworkSpecies()[i].name() << ", " << m_engine.getNetworkSpecies()[j].name() << ") = " << J(i, j) << '\n';
             }
         }
     }
+
 }

src/network/lib/utils/logging.cpp

@@ -16,7 +16,13 @@ std::string gridfire::utils::formatNuclearTimescaleLogString(
     const double T9,
     const double rho
 ) {
-    auto const& timescales = engine.getSpeciesTimescales(Y, T9, rho);
+    auto const& result = engine.getSpeciesTimescales(Y, T9, rho);
+    if (!result) {
+        std::ostringstream ss;
+        ss << "Failed to get species timescales: " << result.error();
+        return ss.str();
+    }
+    const std::unordered_map<fourdst::atomic::Species, double>& timescales = result.value();
 
     // Figure out how wide the "Species" column needs to be:
     std::size_t maxNameLen = std::string_view("Species").size();

subprojects/libcomposition.wrap

@@ -1,4 +1,4 @@
 [wrap-git]
 url = https://github.com/4D-STAR/libcomposition.git
-revision = v1.4.0
+revision = v1.4.1
 depth = 1

tests/graphnet_sandbox/main.cpp

@@ -72,24 +72,30 @@ int main() {
     NetIn netIn;
     netIn.composition = composition;
     netIn.temperature = 1.5e7;
-    netIn.density = 1.5e3;
+    netIn.density = 1.6e2;
     netIn.energy = 0;
-    netIn.tMax = 3.546e17;
+    netIn.tMax = 3.1536e17; // ~ 10Gyr
     netIn.dt0 = 1e-12;
-
+    for (const auto& [symbol, entry] : netIn.composition) {
-    GraphEngine ReaclibEngine(composition, partitionFunction, NetworkBuildDepth::SecondOrder);
-    ReaclibEngine.setPrecomputation(true);
-    ReaclibEngine.setUseReverseReactions(false);
-    MultiscalePartitioningEngineView partitioningView(ReaclibEngine);
-    AdaptiveEngineView adaptiveView(partitioningView);
-
-    solver::DirectNetworkSolver solver(adaptiveView);
-    NetOut netOut;
-    measure_execution_time([&](){netOut = solver.evaluate(netIn);}, "DirectNetworkSolver Evaluation");
-    std::cout << "DirectNetworkSolver completed in " << netOut.num_steps << " steps.\n";
-    std::cout << "Final composition:\n";
-    for (const auto& [symbol, entry] : netOut.composition) {
         std::cout << symbol << ": " << entry.mass_fraction() << "\n";
     }
 
+    // GraphEngine ReaclibEngine(composition, partitionFunction, NetworkBuildDepth::SecondOrder);
+    //
+    // ReaclibEngine.setPrecomputation(true);
+    // ReaclibEngine.setUseReverseReactions(false);
+    // ReaclibEngine.setScreeningModel(screening::ScreeningType::WEAK);
+    //
+    // MultiscalePartitioningEngineView partitioningView(ReaclibEngine);
+    // AdaptiveEngineView adaptiveView(partitioningView);
+    //
+    // solver::DirectNetworkSolver solver(adaptiveView);
+    // NetOut netOut;
+    // measure_execution_time([&](){netOut = solver.evaluate(netIn);}, "DirectNetworkSolver Evaluation");
+    // std::cout << "DirectNetworkSolver completed in " << netOut.num_steps << " steps.\n";
+    // std::cout << "Final composition:\n";
+    // for (const auto& [symbol, entry] : netOut.composition) {
+    //     std::cout << symbol << ": " << entry.mass_fraction() << "\n";
+    // }
+
 }