feat(AdaptiveEngine): system much more stable

dramatically increased stability of jacobian. System is burning again with much more robust physics

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

@@ -8,9 +8,12 @@
 #include "gridfire/engine/types/reporting.h"
 #include "gridfire/engine/types/building.h"
 
+#include "gridfire/expectations/expected_engine.h"
+
 #include <vector>
 #include <unordered_map>
 #include <utility>
+#include <expected>
 
 /**
  * @file engine_abstract.h
@@ -140,14 +143,14 @@ namespace gridfire {
             const std::vector<double>& Y_dynamic,
             double T9,
             double rho
-        ) = 0;
+        ) const = 0;
 
         virtual void generateJacobianMatrix(
             const std::vector<double>& Y_dynamic,
             double T9,
             double rho,
             const SparsityPattern& sparsityPattern
-        ) {
+        ) const {
             throw std::logic_error("Sparsity pattern not supported by this engine.");
         }
 
@@ -251,7 +254,9 @@ namespace gridfire {
          *
          * @post The internal state of the engine is updated to reflect the new conditions.
          */
-        virtual void update(const NetIn& netIn) = 0;
+        virtual fourdst::composition::Composition update(const NetIn &netIn) = 0;
+
+        virtual bool isStale(const NetIn& netIn) = 0;
 
         /**
          * @brief Set the electron screening model.
@@ -296,5 +301,6 @@ namespace gridfire {
         virtual void rebuild(const fourdst::composition::Composition& comp, BuildDepthType depth) {
             throw std::logic_error("Setting network depth not supported by this engine.");
         }
+
     };
 }

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

@@ -34,7 +34,7 @@
 //       this makes extra copies of the species, which is not ideal and could be optimized further.
 //       Even more relevant is the member m_reactionIDMap which makes copies of a REACLIBReaction for each reaction ID.
 //       REACLIBReactions are quite large data structures, so this could be a performance bottleneck.
-
+// static bool isF17 = false;
 namespace gridfire {
     static bool s_debug = false; // Global debug flag for the GraphEngine
     /**
@@ -167,14 +167,14 @@ namespace gridfire {
             const std::vector<double>& Y_dynamic,
             const double T9,
             const double rho
-        ) override;
+        ) const override;
 
         void generateJacobianMatrix(
             const std::vector<double> &Y_dynamic,
             double T9,
             double rho,
             const SparsityPattern &sparsityPattern
-        ) override;
+        ) const override;
 
         /**
          * @brief Generates the stoichiometry matrix for the network.
@@ -274,7 +274,9 @@ namespace gridfire {
             double rho
         ) const override;
 
-        void update(const NetIn& netIn) override;
+        fourdst::composition::Composition update(const NetIn &netIn) override;
+
+        bool isStale(const NetIn &netIn) override;
 
         /**
          * @brief Checks if a given species is involved in the network.
@@ -371,7 +373,6 @@ namespace gridfire {
         void rebuild(const fourdst::composition::Composition& comp, const BuildDepthType depth) override;
 
 
-
     private:
         struct PrecomputedReaction {
             // Forward cacheing
@@ -449,11 +450,12 @@ namespace gridfire {
         std::unordered_map<fourdst::atomic::Species, size_t> m_speciesToIndexMap; ///< Map from species to their index in the stoichiometry matrix.
 
         boost::numeric::ublas::compressed_matrix<int> m_stoichiometryMatrix; ///< Stoichiometry matrix (species x reactions).
-        boost::numeric::ublas::compressed_matrix<double> m_jacobianMatrix; ///< Jacobian matrix (species x species).
 
-        CppAD::ADFun<double> m_rhsADFun; ///< CppAD function for the right-hand side of the ODE.
+        mutable boost::numeric::ublas::compressed_matrix<double> m_jacobianMatrix; ///< Jacobian matrix (species x species).
-        CppAD::sparse_jac_work m_jac_work; ///< Work object for sparse Jacobian calculations.
+        mutable CppAD::ADFun<double> m_rhsADFun; ///< CppAD function for the right-hand side of the ODE.
+        mutable CppAD::sparse_jac_work m_jac_work; ///< Work object for sparse Jacobian calculations.
         CppAD::sparse_rc<std::vector<size_t>> m_full_jacobian_sparsity_pattern; ///< Full sparsity pattern for the Jacobian matrix.
+
         std::vector<std::unique_ptr<AtomicReverseRate>> m_atomicReverseRates;
 
         screening::ScreeningType m_screeningType = screening::ScreeningType::BARE; ///< Screening type for the reaction network. Default to no screening.
@@ -793,25 +795,43 @@ namespace gridfire {
             );
 
             const T molarReactionFlow = forwardMolarReactionFlow - reverseMolarFlow; // Net molar reaction flow
-            std::stringstream ss;
+            // std::stringstream ss;
-            ss << "Forward: " << forwardMolarReactionFlow
+            // ss << "Forward: " << forwardMolarReactionFlow
-               << ", Reverse: " << reverseMolarFlow
+            //    << ", Reverse: " << reverseMolarFlow
-               << ", Net: " << molarReactionFlow;
+            //    << ", Net: " << molarReactionFlow;
-            LOG_TRACE_L2(
+            // LOG_TRACE_L3(
-                m_logger,
+            //     m_logger,
-                "Reaction: {}, {}",
+            //     "Reaction: {}, {}",
-                reaction.peName(),
+            //     reaction.peName(),
-                ss.str()
+            //     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/views/engine_adaptive.h

@@ -76,7 +76,9 @@ namespace gridfire {
          * @see AdaptiveEngineView::constructSpeciesIndexMap()
          * @see AdaptiveEngineView::constructReactionIndexMap()
          */
-        void update(const NetIn& netIn) override;
+        fourdst::composition::Composition update(const NetIn &netIn) override;
+
+        bool isStale(const NetIn& netIn) override;
 
         /**
          * @brief Gets the list of active species in the network.
@@ -123,7 +125,7 @@ namespace gridfire {
             const std::vector<double> &Y_dynamic,
             const double T9,
             const double rho
-        ) override;
+        ) const override;
 
         /**
          * @brief Gets an entry from the Jacobian matrix for the active species.
@@ -438,6 +440,16 @@ namespace gridfire {
             const std::vector<double>& Y_full,
             double maxFlow
         ) const;
+
+        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,
+            const std::vector<fourdst::atomic::Species>& activeSpecies,
+            const reaction::LogicalReactionSet& activeReactions
+        ) const;
         /**
          * @brief Finalizes the set of active species and reactions.
          *

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

@@ -55,7 +55,7 @@ namespace gridfire{
             const std::vector<double>& Y_dynamic,
             const double T9,
             const double rho
-        ) override;
+        ) const override;
         /**
          * @brief Gets an entry from the Jacobian matrix for the active species.
          *
@@ -142,7 +142,9 @@ namespace gridfire{
          *
          * @post If the view was stale, it is rebuilt and is no longer stale.
          */
-        void update(const NetIn &netIn) override;
+        fourdst::composition::Composition update(const NetIn &netIn) override;
+
+        bool isStale(const NetIn& netIn) override;
 
         /**
          * @brief Sets the screening model for the base engine.

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

@@ -6,6 +6,58 @@
 
 #include "unsupported/Eigen/NonLinearOptimization"
 
+namespace gridfire {
+    struct QSECacheConfig {
+        double T9_tol; ///< Absolute tolerance to produce the same hash for T9.
+        double rho_tol; ///< Absolute tolerance to produce the same hash for rho.
+        double Yi_tol; ///< Absolute tolerance to produce the same hash for species abundances.
+    };
+
+    struct QSECacheKey {
+        double m_T9;
+        double m_rho;
+        std::vector<double> m_Y; ///< Note that the ordering of Y must match the dynamic species indices in the view.
+
+        std::size_t m_hash = 0; ///< Precomputed hash value for this key.
+
+        // TODO: We should probably sort out how to adjust these from absolute to relative tolerances.
+        QSECacheConfig m_cacheConfig = {
+            1e-3, // Default tolerance for T9
+            1e-1, // Default tolerance for rho
+            1e-3  // Default tolerance for species abundances
+        };
+
+        QSECacheKey(
+            const double T9,
+            const double rho,
+            const std::vector<double>& Y
+        );
+
+        size_t hash() const;
+
+        static long bin(double value, double tol);
+
+        bool operator==(const QSECacheKey& other) const;
+
+    };
+}
+
+// Needs to be in this order (splitting gridfire namespace up) to avoid some issues with forward declarations and the () operator.
+namespace std {
+    template <>
+    struct hash<gridfire::QSECacheKey> {
+        /**
+         * @brief Computes the hash of a QSECacheKey.
+         * @param key The QSECacheKey to hash.
+         * @return The pre-computed hash value of the key.
+         */
+        size_t operator()(const gridfire::QSECacheKey& key) const noexcept {
+            // The hash is pre-computed, so we just return it.
+            return key.m_hash;
+        }
+    };
+} // namespace std
+
 namespace gridfire {
     class MultiscalePartitioningEngineView final: public DynamicEngine, public EngineView<DynamicEngine> {
         typedef std::tuple<std::vector<fourdst::atomic::Species>, std::vector<size_t>, std::vector<fourdst::atomic::Species>, std::vector<size_t>> QSEPartition;
@@ -15,20 +67,20 @@ namespace gridfire {
         [[nodiscard]] const std::vector<fourdst::atomic::Species> & getNetworkSpecies() const override;
 
         [[nodiscard]] StepDerivatives<double> calculateRHSAndEnergy(
-            const std::vector<double> &Y,
+            const std::vector<double> &Y_full,
             double T9,
             double rho
         ) const override;
 
         void generateJacobianMatrix(
-            const std::vector<double> &Y_dynamic,
+            const std::vector<double> &Y_full,
             double T9,
             double rho
-        ) override;
+        ) const override;
 
         [[nodiscard]] double getJacobianMatrixEntry(
-            int i,
+            int i_full,
-            int j
+            int j_full
         ) const override;
 
         void generateStoichiometryMatrix() override;
@@ -40,7 +92,7 @@ namespace gridfire {
 
         [[nodiscard]] double calculateMolarReactionFlow(
             const reaction::Reaction &reaction,
-            const std::vector<double> &Y,
+            const std::vector<double> &Y_full,
             double T9,
             double rho
         ) const override;
@@ -53,10 +105,12 @@ namespace gridfire {
             double rho
         ) const override;
 
-        void update(
+        fourdst::composition::Composition update(
             const NetIn &netIn
         ) override;
 
+        bool isStale(const NetIn& netIn) override;
+
         void setScreeningModel(
             screening::ScreeningType model
         ) override;
@@ -181,6 +235,72 @@ namespace gridfire {
             int df(const InputType& v_qse, JacobianType& J_qse) const;
         };
 
+
+        struct CacheStats {
+            enum class operators {
+                CalculateRHSAndEnergy,
+                GenerateJacobianMatrix,
+                CalculateMolarReactionFlow,
+                GetSpeciesTimescales,
+                Other,
+                All
+            };
+
+            std::map<operators, std::string> operatorsNameMap = {
+                {operators::CalculateRHSAndEnergy, "calculateRHSAndEnergy"},
+                {operators::GenerateJacobianMatrix, "generateJacobianMatrix"},
+                {operators::CalculateMolarReactionFlow, "calculateMolarReactionFlow"},
+                {operators::GetSpeciesTimescales, "getSpeciesTimescales"},
+                {operators::Other, "other"}
+            };
+
+            size_t m_hit = 0;
+            size_t m_miss = 0;
+
+            std::map<operators, size_t> m_operatorHits = {
+                {operators::CalculateRHSAndEnergy, 0},
+                {operators::GenerateJacobianMatrix, 0},
+                {operators::CalculateMolarReactionFlow, 0},
+                {operators::GetSpeciesTimescales, 0},
+                {operators::Other, 0}
+            };
+
+
+            friend std::ostream& operator<<(std::ostream& os, const CacheStats& stats) {
+                os << "CacheStats(hit: " << stats.m_hit << ", miss: " << stats.m_miss << ")";
+                return os;
+            }
+
+            void hit(const operators op=operators::Other) {
+                if (op==operators::All) {
+                    throw std::invalid_argument("Cannot use 'All' as an operator for hit/miss.");
+                }
+                m_hit++;
+                m_operatorHits[op]++;
+            }
+            void miss(const operators op=operators::Other) {
+                if (op==operators::All) {
+                    throw std::invalid_argument("Cannot use 'All' as an operator for hit/miss.");
+                }
+                m_miss++;
+                m_operatorHits[op]++;
+            }
+
+            [[nodiscard]] size_t hits(const operators op=operators::All) const {
+                if (op==operators::All) {
+                    return m_hit;
+                }
+                return m_operatorHits.at(op);
+            }
+            [[nodiscard]] size_t misses(const operators op=operators::All) const {
+                if (op==operators::All) {
+                    return m_miss;
+                }
+                return m_operatorHits.at(op);
+            }
+        };
+
+
     private:
         quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
         GraphEngine& m_baseEngine; ///< The base engine to which this view delegates calculations.
@@ -189,7 +309,16 @@ namespace gridfire {
         std::vector<size_t> m_dynamic_species_indices; ///< Indices mapping the dynamic species back to the master engine's list.
         std::vector<fourdst::atomic::Species> m_algebraic_species; ///< Species that are algebraic in the QSE groups.
         std::vector<size_t> m_algebraic_species_indices; ///< Indices of algebraic species in the full network.
-        std::unordered_map<size_t, std::vector<size_t>> m_connectivity_graph;
+
+        std::vector<size_t> m_activeSpeciesIndices; ///< Indices of active species in the full network.
+        std::vector<size_t> m_activeReactionIndices; ///< Indices of active reactions in the full network.
+
+        // TODO: Enhance the hashing for the cache to consider not just T and rho but also the current abundance in some careful way that automatically ignores small changes (i.e. network should only be repartitioned sometimes)
+        std::unordered_map<QSECacheKey, std::vector<double>> m_qse_abundance_cache; ///< Cache for QSE abundances based on T9 and rho.
+        mutable CacheStats m_cacheStats; ///< Statistics for the QSE abundance cache.
+
+
+
 
     private:
         std::vector<std::vector<size_t>> partitionByTimescale(
@@ -230,3 +359,4 @@ namespace gridfire {
         ) const;
     };
 }
+

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

@@ -0,0 +1,34 @@
+#pragma once
+
+#include <exception>
+#include <string>
+#include <iostream>
+
+namespace gridfire::exceptions {
+    class EngineError : public std::exception {};
+
+    class StaleEngineError final : public EngineError {
+    public:
+        explicit StaleEngineError(const std::string& message)
+            : m_message(message) {}
+
+        const char* what() const noexcept override{
+            return m_message.c_str();
+        }
+    private:
+        std::string m_message;
+    };
+
+    class FailedToPartitionEngineError final : public EngineError {
+    public:
+        explicit FailedToPartitionEngineError(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

@@ -0,0 +1,34 @@
+#pragma once
+
+#include <string>
+#include <iostream>
+
+namespace gridfire::expectations {
+    enum class EngineErrorTypes {
+        FAILURE,
+        INDEX,
+        STALE
+    };
+
+    struct EngineError {
+        std::string m_message;
+        EngineErrorTypes type = EngineErrorTypes::FAILURE;
+        friend std::ostream& operator<<(std::ostream& os, const EngineError& e) {
+            os << "EngineError: " << e.m_message;
+            return os;
+        }
+    };
+
+    struct EngineIndexError : EngineError {
+        int m_index;
+        EngineErrorTypes type = EngineErrorTypes::INDEX;
+        friend std::ostream& operator<<(std::ostream& os, const EngineIndexError& e) {
+            os << "EngineIndexError: " << e.m_message << " at index " << e.m_index;
+            return os;
+        }
+    };
+
+    struct StaleEngineError : EngineError {
+        EngineErrorTypes type = EngineErrorTypes::STALE;
+    };
+}

src/network/lib/engine/engine_graph.cpp

@@ -65,6 +65,7 @@ namespace gridfire {
             bare_rates.reserve(m_reactions.size());
             bare_reverse_rates.reserve(m_reactions.size());
 
+            // TODO: Add cache to this
             for (const auto& reaction: m_reactions) {
                 bare_rates.push_back(reaction.calculate_rate(T9));
                 bare_reverse_rates.push_back(calculateReverseRate(reaction, T9));
@@ -77,7 +78,6 @@ namespace gridfire {
         }
     }
 
-
     void GraphEngine::syncInternalMaps() {
         LOG_INFO(m_logger, "Synchronizing internal maps for REACLIB graph network (serif::network::GraphNetwork)...");
         collectNetworkSpecies();
@@ -650,7 +650,7 @@ namespace gridfire {
         const std::vector<double> &Y_dynamic,
         const double T9,
         const double rho
-    ) {
+    ) const {
 
         LOG_TRACE_L1(m_logger, "Generating jacobian matrix for T9={}, rho={}..", T9, rho);
         const size_t numSpecies = m_networkSpecies.size();
@@ -658,7 +658,7 @@ namespace gridfire {
         // 1. Pack the input variables into a vector for CppAD
         std::vector<double> adInput(numSpecies + 2, 0.0); // +2 for T9 and rho
         for (size_t i = 0; i < numSpecies; ++i) {
-            adInput[i] = Y_dynamic[i];
+            adInput[i] = std::max(Y_dynamic[i], 1e-99); // regularize the jacobian...
         }
         adInput[numSpecies]     = T9;  // T9
         adInput[numSpecies + 1] = rho; // rho
@@ -685,28 +685,36 @@ namespace gridfire {
         for (size_t i = 0; i < numSpecies; ++i) {
             for (size_t j = 0; j < numSpecies; ++j) {
                 const double value = dotY[i * (numSpecies + 2) + j];
-                if (std::abs(value) > MIN_JACOBIAN_THRESHOLD) {
+                if (std::abs(value) > MIN_JACOBIAN_THRESHOLD || i == j) { // Always keep diagonal elements to avoid pathological stiffness
                     m_jacobianMatrix(i, j) = value;
                 }
             }
         }
-        // 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) {
-        //             for (size_t j = 0; j < m_jacobianMatrix.size2(); ++j) {
+                    ss << getNetworkSpecies()[i].name();
-        //                 ss << m_jacobianMatrix(i, j);
+                    if (i < m_jacobianMatrix.size1() - 1) {
-        //                 if (j < m_jacobianMatrix.size2() - 1) {
+                        ss << ", ";
-        //                     ss << ", ";
+                    }
-        //                 }
+                }
-        //             }
+                ss << "\n";
-        //             ss << "\n";
+                for (size_t i = 0; i < m_jacobianMatrix.size1(); ++i) {
-        //         }
+                    ss << getNetworkSpecies()[i].name() << ": ";
-        //         return ss.str();
+                    for (size_t j = 0; j < m_jacobianMatrix.size2(); ++j) {
-        //     }());
+                        ss << m_jacobianMatrix(i, j);
+                        if (j < m_jacobianMatrix.size2() - 1) {
+                            ss << ", ";
+                        }
+                    }
+                    ss << "\n";
+                }
+                return ss.str();
+            }());
         LOG_TRACE_L1(m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2());
     }
 
@@ -715,7 +723,7 @@ namespace gridfire {
         const double T9,
         const double rho,
         const SparsityPattern &sparsityPattern
-    ) {
+    ) const {
         // --- Pack the input variables into a vector for CppAD ---
         const size_t numSpecies = m_networkSpecies.size();
         std::vector<double> x(numSpecies + 2, 0.0);
@@ -768,7 +776,7 @@ namespace gridfire {
     }
 
     double GraphEngine::getJacobianMatrixEntry(const int i, const int j) const {
-        LOG_TRACE_L3(m_logger, "Getting jacobian matrix entry for {},{} = {}", i, j, m_jacobianMatrix(i, j));
+        // LOG_TRACE_L3(m_logger, "Getting jacobian matrix entry for {},{} = {}", i, j, m_jacobianMatrix(i, j));
         return m_jacobianMatrix(i, j);
     }
 
@@ -906,8 +914,12 @@ namespace gridfire {
         return speciesTimescales;
     }
 
-    void GraphEngine::update(const NetIn &netIn) {
+    fourdst::composition::Composition GraphEngine::update(const NetIn &netIn) {
-        // No-op for GraphEngine, as it does not support manually triggering updates.
+        return netIn.composition;
+    }
+
+    bool GraphEngine::isStale(const NetIn &netIn) {
+        return false;
     }
 
     void GraphEngine::recordADTape() {

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

@@ -1,7 +1,9 @@
-#include "../../../include/gridfire/engine/views/engine_adaptive.h"
+#include "gridfire/engine/views/engine_adaptive.h"
 
 #include <ranges>
 #include <queue>
+#include <algorithm>
+
 
 #include "gridfire/network.h"
 
@@ -82,7 +84,9 @@ namespace gridfire {
         return reactionIndexMap;
     }
 
-    void AdaptiveEngineView::update(const NetIn& netIn) {
+    fourdst::composition::Composition AdaptiveEngineView::update(const NetIn &netIn) {
+        fourdst::composition::Composition baseUpdatedComposition = m_baseEngine.update(netIn);
+
         LOG_TRACE_L1(m_logger, "Updating AdaptiveEngineView with new network input...");
 
         std::vector<double> Y_Full;
@@ -104,6 +108,17 @@ namespace gridfire {
 
         finalizeActiveSet(finalReactions);
 
+        // auto [rescuedReactions, rescuedSpecies] = rescueEdgeSpeciesDestructionChannel(allFlows, Y_Full, netIn.temperature/1e9, netIn.density, m_activeSpecies, m_activeReactions);
+        //
+        // for (const auto& reactionPtr : rescuedReactions) {
+        //     m_activeReactions.add_reaction(*reactionPtr);
+        // }
+        //
+        // for (const auto& species : rescuedSpecies) {
+        //     if (!std::ranges::contains(m_activeSpecies, species)) {
+        //         m_activeSpecies.push_back(species);
+        //     }
+        // }
 
         m_speciesIndexMap = constructSpeciesIndexMap();
         m_reactionIndexMap = constructReactionIndexMap();
@@ -111,6 +126,12 @@ namespace gridfire {
         m_isStale = false;
 
         LOG_INFO(m_logger, "AdaptiveEngineView updated successfully with {} active species and {} active reactions.", m_activeSpecies.size(), m_activeReactions.size());
+
+        return baseUpdatedComposition;
+    }
+
+    bool AdaptiveEngineView::isStale(const NetIn &netIn) {
+        return m_isStale;
     }
 
     const std::vector<Species> & AdaptiveEngineView::getNetworkSpecies() const {
@@ -139,7 +160,7 @@ namespace gridfire {
         const std::vector<double> &Y_dynamic,
         const double T9,
         const double rho
-    ) {
+    ) const {
         validateState();
         const auto Y_full = mapCulledToFull(Y_dynamic);
 
@@ -288,6 +309,7 @@ namespace gridfire {
         }
     }
 
+    // TODO: Change this to use a return value instead of an output parameter.
     std::vector<AdaptiveEngineView::ReactionFlow> AdaptiveEngineView::calculateAllReactionFlows(
         const NetIn &netIn,
         std::vector<double> &out_Y_Full
@@ -314,7 +336,7 @@ namespace gridfire {
         for (const auto& reaction : fullReactionSet) {
             const double flow = m_baseEngine.calculateMolarReactionFlow(reaction, out_Y_Full, T9, rho);
             reactionFlows.push_back({&reaction, flow});
-            LOG_TRACE_L2(m_logger, "Reaction '{}' has flow rate: {:0.3E} [mol/s]", reaction.id(), flow);
+            LOG_TRACE_L1(m_logger, "Reaction '{}' has flow rate: {:0.3E} [mol/s/g]", reaction.id(), flow);
         }
         return reactionFlows;
     }
@@ -380,12 +402,12 @@ namespace gridfire {
                 keepReaction = true;
             } else {
                 bool zero_flow_due_to_reachable_reactants = false;
-                if (flowRate < 1e-99) {
+                if (flowRate < 1e-99 && flowRate > 0.0) {
                     for (const auto& reactant: reactionPtr->reactants()) {
                         const auto it = std::ranges::find(m_baseEngine.getNetworkSpecies(), reactant);
                         const size_t index = std::distance(m_baseEngine.getNetworkSpecies().begin(), it);
                         if (Y_full[index] < 1e-99 && reachableSpecies.contains(reactant)) {
-                            LOG_TRACE_L2(m_logger, "Maintaining reaction '{}' with zero flow due to reachable reactant '{}'.", reactionPtr->id(), reactant.name());
+                            LOG_TRACE_L1(m_logger, "Maintaining reaction '{}' with low flow ({:0.3E} [mol/s/g]) due to reachable reactant '{}'.", reactionPtr->id(), flowRate, reactant.name());
                             zero_flow_due_to_reachable_reactants = true;
                             break;
                         }
@@ -405,6 +427,158 @@ namespace gridfire {
         return culledReactions;
     }
 
+    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);
+        std::set<Species> onlyProducedSpecies;
+        for (const auto& reaction : activeReactions) {
+            const std::vector<Species> products = reaction.products();
+            onlyProducedSpecies.insert(products.begin(), products.end());
+        }
+
+        // Remove species that are consumed by any one of the active reactions.
+        std::erase_if(
+            onlyProducedSpecies,
+            [&](const Species &species) {
+                for (const auto& reaction : activeReactions) {
+                    if (reaction.contains_reactant(species)) {
+                        return true; // If any active reaction consumes the species then erase it from the set.
+                    }
+                }
+                return false;
+            }
+        );
+
+        // Remove species that have a non-zero timescale (these are expected to be +inf as they should be the equilibrium species if using with a MultiscalePartitioningEngineView)
+        std::erase_if(
+            onlyProducedSpecies,
+            [&](const Species &species) {
+                return std::isinf(timescales.at(species));
+            }
+        );
+
+        LOG_TRACE_L1(
+            m_logger,
+            "Found {} species {}that are produced but not consumed in any reaction and do not have an inf timescale (those are expected to be equilibrium species and do not contribute to the stiffness of the network).",
+            onlyProducedSpecies.size(),
+            [&]() -> std::string {
+                std::ostringstream ss;
+                if (onlyProducedSpecies.empty()) {
+                    return "";
+                }
+                int count = 0;
+                ss << "(";
+                for (const auto& species : onlyProducedSpecies) {
+                    ss << species.name();
+                    if (count < onlyProducedSpecies.size() - 1) {
+                        ss << ", ";
+                    }
+                    count++;
+                }
+                ss << ") ";
+                return ss.str();
+            }()
+        );
+
+        std::unordered_map<Species, const reaction::LogicalReaction*> reactionsToRescue;
+        for (const auto& species : onlyProducedSpecies) {
+            double maxSpeciesConsumptionRate = 0.0;
+            for (const auto& reaction : m_baseEngine.getNetworkReactions()) {
+                const bool speciesToCheckIsConsumed = reaction.contains_reactant(species);
+                if (!speciesToCheckIsConsumed) {
+                    continue; // If the species is not consumed by this reaction, skip it.
+                }
+                bool allOtherReactantsAreAvailable = true;
+                for (const auto& reactant : reaction.reactants()) {
+                    const bool reactantIsAvailable = std::ranges::contains(activeSpecies, reactant);
+                    if (!reactantIsAvailable && reactant != species) {
+                        allOtherReactantsAreAvailable = false;
+                    }
+                }
+                if (allOtherReactantsAreAvailable && speciesToCheckIsConsumed) {
+                    double rate = reaction.calculate_rate(T9);
+                    if (rate > maxSpeciesConsumptionRate) {
+                        maxSpeciesConsumptionRate = rate;
+                        reactionsToRescue[species] = &reaction;
+                    }
+                }
+            }
+        }
+        LOG_TRACE_L1(
+            m_logger,
+            "Rescuing {} {}reactions",
+            reactionsToRescue.size(),
+            [&]() -> std::string {
+                std::ostringstream ss;
+                if (reactionsToRescue.empty()) {
+                    return "";
+                }
+                int count = 0;
+                ss << "(";
+                for (const auto& [species, reaction] : reactionsToRescue) {
+                    ss << reaction->id();
+                    if (count < reactionsToRescue.size() - 1) {
+                        ss << ", ";
+                    }
+                    count++;
+                }
+                ss << ") ";
+                return ss.str();
+            }()
+        );
+
+        LOG_TRACE_L1(
+            m_logger,
+            "Timescale adjustments due to reaction rescue: {}",
+            [&]() -> std::string {
+                std::stringstream ss;
+                if (reactionsToRescue.empty()) {
+                    return "No reactions rescued...";
+                }
+                int count = 0;
+                for (const auto& [species, reaction] : reactionsToRescue) {
+                    ss << "(Species: " << species.name() << " started with a timescale of " << timescales.at(species);
+                    ss << ", rescued by reaction: " << reaction->id();
+                    ss << " whose product timescales are --- [";
+                    int icount = 0;
+                    for (const auto& product : reaction->products()) {
+                        ss << product.name() << ": " << timescales.at(product);
+                        if (icount < reaction->products().size() - 1) {
+                            ss << ", ";
+                        }
+                        icount++;
+                    }
+                    ss << "])";
+
+                    if (count < reactionsToRescue.size() - 1) {
+                        ss << ", ";
+                    }
+                    count++;
+                }
+
+                return ss.str();
+            }()
+        );
+
+        RescueSet rescueSet;
+        std::unordered_set<const reaction::LogicalReaction*> newReactions;
+        std::unordered_set<Species> newSpecies;
+
+        for (const auto &reactionPtr: reactionsToRescue | std::views::values) {
+            newReactions.insert(reactionPtr);
+            for (const auto& product : reactionPtr->products()) {
+                newSpecies.insert(product);
+            }
+        }
+        return {std::move(newReactions), std::move(newSpecies)};
+    }
+
     void AdaptiveEngineView::finalizeActiveSet(
         const std::vector<const reaction::LogicalReaction *> &finalReactions
     ) {
@@ -413,13 +587,24 @@ namespace gridfire {
         for (const auto* reactionPtr: finalReactions) {
             m_activeReactions.add_reaction(*reactionPtr);
             for (const auto& reactant : reactionPtr->reactants()) {
+                if (!finalSpeciesSet.contains(reactant)) {
+                    LOG_TRACE_L1(m_logger, "Adding reactant '{}' to active species set through reaction {}.", reactant.name(), reactionPtr->id());
+                } else {
+                    LOG_TRACE_L1(m_logger, "Reactant '{}' already in active species set through another reaction.", reactant.name());
+                }
                 finalSpeciesSet.insert(reactant);
             }
             for (const auto& product : reactionPtr->products()) {
+                if (!finalSpeciesSet.contains(product)) {
+                    LOG_TRACE_L1(m_logger, "Adding product '{}' to active species set through reaction {}.", product.name(), reactionPtr->id());
+                } else {
+                    LOG_TRACE_L1(m_logger, "Product '{}' already in active species set through another reaction.", product.name());
+                }
                 finalSpeciesSet.insert(product);
             }
         }
-        m_activeSpecies.assign(finalSpeciesSet.begin(), finalSpeciesSet.end());
+        m_activeSpecies.clear();
+        m_activeSpecies = std::vector<Species>(finalSpeciesSet.begin(), finalSpeciesSet.end());
         std::ranges::sort(
             m_activeSpecies,
             [](const Species &a, const Species &b) { return a.mass() < b.mass(); }

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

@@ -100,7 +100,7 @@ namespace gridfire {
         const std::vector<double> &Y_dynamic,
         const double T9,
         const double rho
-    ) {
+    ) const {
         validateNetworkState();
 
         const auto Y_full = mapViewToFull(Y_dynamic);
@@ -178,8 +178,12 @@ namespace gridfire {
         return definedTimescales;
     }
 
-    void DefinedEngineView::update(const NetIn &netIn) {
+    fourdst::composition::Composition DefinedEngineView::update(const NetIn &netIn) {
-        return;
+        return m_baseEngine.update(netIn);
+    }
+
+    bool DefinedEngineView::isStale(const NetIn &netIn) {
+        return false;
     }
 
 

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

@@ -1,4 +1,5 @@
 #include "gridfire/engine/views/engine_multiscale.h"
+#include "gridfire/exceptions/error_engine.h"
 
 #include <numeric>
 
@@ -8,9 +9,33 @@
 #include <unordered_map>
 #include <unordered_set>
 
+#include <ranges>
+
 #include "quill/LogMacros.h"
 #include "quill/Logger.h"
 
+namespace {
+    std::vector<double> packCompositionToVector(const fourdst::composition::Composition& composition, const gridfire::GraphEngine& engine) {
+        std::vector<double> Y(engine.getNetworkSpecies().size(), 0.0);
+        const auto& allSpecies = engine.getNetworkSpecies();
+        for (size_t i = 0; i < allSpecies.size(); ++i) {
+            const auto& species = allSpecies[i];
+            if (!composition.contains(species)) {
+                Y[i] = 0.0; // Species not in the composition, set to zero
+            } else {
+                Y[i] = composition.getMolarAbundance(species);
+            }
+        }
+        return Y;
+    }
+
+    template <class T>
+    void hash_combine(std::size_t& seed, const T& v) {
+        std::hash<T> hashed;
+        seed ^= hashed(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
+    }
+}
+
 namespace gridfire {
     static int s_operator_parens_called = 0;
     using fourdst::atomic::Species;
@@ -20,35 +45,84 @@ namespace gridfire {
     ) : m_baseEngine(baseEngine) {}
 
     const std::vector<Species> & MultiscalePartitioningEngineView::getNetworkSpecies() const {
-        return m_dynamic_species;
+        return m_baseEngine.getNetworkSpecies();
     }
 
     StepDerivatives<double> MultiscalePartitioningEngineView::calculateRHSAndEnergy(
-        const std::vector<double> &Y,
+        const std::vector<double> &Y_full,
         const double T9,
         const double rho
     ) const {
-        return m_baseEngine.calculateRHSAndEnergy(Y, T9, rho);
+        if (Y_full.size() != getNetworkSpecies().size()) {
+            LOG_ERROR(
+                m_logger,
+                "Y_full size ({}) does not match the number of species in the network ({})",
+                Y_full.size(),
+                getNetworkSpecies().size()
+            );
+            throw std::runtime_error("Y_full size does not match the number of species in the network. See logs for more details...");
+        }
+
+        // Check the cache to see if the network needs to be repartitioned. Note that the QSECacheKey manages binning of T9, rho, and Y_full to ensure that small changes (which would likely not result in a repartitioning) do not trigger a cache miss.
+        const QSECacheKey key(T9, rho, Y_full);
+        if (! m_qse_abundance_cache.contains(key)) {
+            m_cacheStats.miss(CacheStats::operators::CalculateRHSAndEnergy);
+            LOG_ERROR(
+                m_logger,
+                "QSE abundance cache miss for T9 = {}, rho = {} (misses: {}, hits: {}). calculateRHSAndEnergy does not receive sufficient context to partition and stabilize the network. Throwing an error which should be caught by the caller and trigger a re-partition stage.",
+                T9,
+                rho,
+                m_cacheStats.misses(),
+                m_cacheStats.hits()
+            );
+            throw exceptions::StaleEngineError("QSE Cache Miss while lacking context for partitioning. This should be caught by the caller and trigger a re-partition stage.");
+        }
+        m_cacheStats.hit(CacheStats::operators::CalculateRHSAndEnergy);
+        StepDerivatives<double> deriv = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho);
+
+        for (size_t i = 0; i < m_algebraic_species_indices.size(); ++i) {
+            const size_t species_index = m_algebraic_species_indices[i];
+            deriv.dydt[species_index] = 0.0; // Fix the algebraic species to the equilibrium abundances we calculate.
+        }
+        return deriv;
     }
 
     void MultiscalePartitioningEngineView::generateJacobianMatrix(
-        const std::vector<double> &Y_dynamic,
+        const std::vector<double> &Y_full,
         const double T9,
         const double rho
-    ) {
+    ) const {
-        std::vector<double> Y_full(m_baseEngine.getNetworkSpecies().size(), 0.0);
+        const QSECacheKey key(T9, rho, Y_full);
-        for (size_t i = 0; i < m_dynamic_species_indices.size(); ++i) {
+        if (!m_qse_abundance_cache.contains(key)) {
-            Y_full[m_dynamic_species_indices[i]] = Y_dynamic[i];
+            m_cacheStats.miss(CacheStats::operators::GenerateJacobianMatrix);
+            LOG_ERROR(
+                m_logger,
+                "QSE abundance cache miss for T9 = {}, rho = {} (misses: {}, hits: {}). generateJacobianMatrix does not receive sufficient context to partition and stabilize the network. Throwing an error which should be caught by the caller and trigger a re-partition stage.",
+                T9,
+                rho,
+                m_cacheStats.misses(),
+                m_cacheStats.hits()
+            );
+            throw exceptions::StaleEngineError("QSE Cache Miss while lacking context for partitioning. This should be caught by the caller and trigger a re-partition stage.");
         }
-        // solveQSEAbundances(Y_full, T9, rho);
+        m_cacheStats.hit(CacheStats::operators::GenerateJacobianMatrix);
-        m_baseEngine.generateJacobianMatrix(Y_dynamic, T9, rho);
+
+        // TODO: Add sparsity pattern to this to prevent base engine from doing unnecessary work.
+        m_baseEngine.generateJacobianMatrix(Y_full, T9, rho);
     }
 
     double MultiscalePartitioningEngineView::getJacobianMatrixEntry(
-        const int i,
+        const int i_full,
-        const int j
+        const int j_full
     ) const {
-        return m_baseEngine.getJacobianMatrixEntry(i, j);
+        // // Check if the species we are differentiating with respect to is algebraic or dynamic. If it is algebraic we can reduce the work significantly...
+        // if (std::ranges::contains(m_algebraic_species_indices, j_full)) {
+        //     const auto& species = m_baseEngine.getNetworkSpecies()[j_full];
+        //     // If j is algebraic, we can return 0.0 since the Jacobian entry for algebraic species is always zero.
+        //     return 0.0;
+        // }
+        // Otherwise we need to query the full jacobian
+        return m_baseEngine.getJacobianMatrixEntry(i_full, j_full);
     }
 
     void MultiscalePartitioningEngineView::generateStoichiometryMatrix() {
@@ -64,11 +138,35 @@ namespace gridfire {
 
     double MultiscalePartitioningEngineView::calculateMolarReactionFlow(
         const reaction::Reaction &reaction,
-        const std::vector<double> &Y,
+        const std::vector<double> &Y_full,
         const double T9,
         const double rho
     ) const {
-        return m_baseEngine.calculateMolarReactionFlow(reaction, Y, T9, rho);
+        const auto key = QSECacheKey(T9, rho, Y_full);
+        if (!m_qse_abundance_cache.contains(key)) {
+            m_cacheStats.miss(CacheStats::operators::CalculateMolarReactionFlow);
+            LOG_ERROR(
+                m_logger,
+                "QSE abundance cache miss for T9 = {}, rho = {} (misses: {}, hits: {}). calculateMolarReactionFlow does not receive sufficient context to partition and stabilize the network. Throwing an error which should be caught by the caller and trigger a re-partition stage.",
+                T9,
+                rho,
+                m_cacheStats.misses(),
+                m_cacheStats.hits()
+            );
+            throw exceptions::StaleEngineError("QSE Cache Miss while lacking context for partitioning. This should be caught by the caller and trigger a re-partition stage.");
+        }
+        m_cacheStats.hit(CacheStats::operators::CalculateMolarReactionFlow);
+        std::vector<double> Y_algebraic = m_qse_abundance_cache.at(key);
+
+        assert(Y_algebraic.size() == m_algebraic_species_indices.size());
+
+        // Fix the algebraic species to the equilibrium abundances we calculate.
+        std::vector<double> Y_mutable = Y_full;
+        for (const auto& [index, Yi] : std::views::zip(m_algebraic_species_indices, Y_algebraic)) {
+            Y_mutable[index] = Yi;
+
+        }
+        return m_baseEngine.calculateMolarReactionFlow(reaction, Y_mutable, T9, rho);
     }
 
     const reaction::LogicalReactionSet & MultiscalePartitioningEngineView::getNetworkReactions() const {
@@ -80,11 +178,60 @@ namespace gridfire {
         const double T9,
         const double rho
     ) const {
-        return m_baseEngine.getSpeciesTimescales(Y, T9, rho);
+        const auto key = QSECacheKey(T9, rho, Y);
+        if (!m_qse_abundance_cache.contains(key)) {
+            m_cacheStats.miss(CacheStats::operators::GetSpeciesTimescales);
+            LOG_ERROR(
+                m_logger,
+                "QSE abundance cache miss for T9 = {}, rho = {} (misses: {}, hits: {}). getSpeciesTimescales does not receive sufficient context to partition and stabilize the network. Throwing an error which should be caught by the caller and trigger a re-partition stage.",
+                T9,
+                rho,
+                m_cacheStats.misses(),
+                m_cacheStats.hits()
+            );
+            throw exceptions::StaleEngineError("QSE Cache Miss while lacking context for partitioning. This should be caught by the caller and trigger a re-partition stage.");
+        }
+        m_cacheStats.hit(CacheStats::operators::GetSpeciesTimescales);
+        std::unordered_map<Species, double> speciesTimescales = m_baseEngine.getSpeciesTimescales(Y, T9, rho);
+        for (const auto& algebraicSpecies : m_algebraic_species) {
+            speciesTimescales[algebraicSpecies] = std::numeric_limits<double>::infinity(); // Algebraic species have infinite timescales.
+        }
+        return speciesTimescales;
     }
 
-    void MultiscalePartitioningEngineView::update(const NetIn &netIn) {
+    fourdst::composition::Composition MultiscalePartitioningEngineView::update(const NetIn &netIn) {
-        return; // call partition eventually
+        const fourdst::composition::Composition baseUpdatedComposition = m_baseEngine.update(netIn);
+
+        const auto key = QSECacheKey(
+            netIn.temperature / 1.0e9, // Convert temperature from Kelvin to T9 (T9 = T / 1e9)
+            netIn.density,
+            packCompositionToVector(baseUpdatedComposition, m_baseEngine)
+        );
+        if (m_qse_abundance_cache.contains(key)) {
+            return baseUpdatedComposition;
+        }
+        NetIn baseUpdatedNetIn = netIn;
+        baseUpdatedNetIn.composition = baseUpdatedComposition;
+        const fourdst::composition::Composition equilibratedComposition = equilibrateNetwork(baseUpdatedNetIn);
+        std::vector<double> Y_algebraic(m_algebraic_species_indices.size(), 0.0);
+        for (size_t i = 0; i < m_algebraic_species_indices.size(); ++i) {
+            const size_t species_index = m_algebraic_species_indices[i];
+            Y_algebraic[i] = equilibratedComposition.getMolarAbundance(m_baseEngine.getNetworkSpecies()[species_index]);
+        }
+        m_qse_abundance_cache[key] = std::move(Y_algebraic);
+        return equilibratedComposition;
+    }
+
+    bool MultiscalePartitioningEngineView::isStale(const NetIn &netIn) {
+        const auto key = QSECacheKey(
+            netIn.temperature,
+            netIn.density,
+            packCompositionToVector(netIn.composition, m_baseEngine)
+        );
+        if (m_qse_abundance_cache.contains(key)) {
+            return false; // The cache hit indicates the engine is not stale for the given conditions.
+        }
+        return true;
     }
 
     void MultiscalePartitioningEngineView::setScreeningModel(
@@ -112,7 +259,6 @@ namespace gridfire {
         m_qse_groups.clear();
         m_dynamic_species.clear();
         m_dynamic_species_indices.clear();
-        m_connectivity_graph.clear();
         m_algebraic_species.clear();
         m_algebraic_species_indices.clear();
 
@@ -223,16 +369,7 @@ namespace gridfire {
     void MultiscalePartitioningEngineView::partitionNetwork(
         const NetIn &netIn
     ) {
-        std::vector<double> Y(m_baseEngine.getNetworkSpecies().size(), 0.0);
+        const std::vector<double> Y = packCompositionToVector(netIn.composition, m_baseEngine);
-        for (const auto& [symbol, entry]: netIn.composition ) {
-            if (!m_baseEngine.involvesSpecies(entry.isotope())) {
-                LOG_WARNING(m_logger, "Species {} is not part of the network. This is likely due to the base network not being deep enough. For rare species this may not be an issue. Skipping...", entry.isotope().name());
-                continue; // Skip species not in the network
-            }
-            const size_t species_index = m_baseEngine.getSpeciesIndex(entry.isotope());
-            Y[species_index] = netIn.composition.getMolarAbundance(symbol);
-        }
-
         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
 
@@ -300,9 +437,9 @@ namespace gridfire {
 
             // Determine color based on category. A species can be a seed and also in the core dynamic group.
             // The more specific category (algebraic, then seed) takes precedence.
-            if (algebraic_indices.count(i)) {
+            if (algebraic_indices.contains(i)) {
                 fillcolor = "#e0f2fe"; // Light Blue: Algebraic (in QSE)
-            } else if (seed_indices.count(i)) {
+            } else if (seed_indices.contains(i)) {
                 fillcolor = "#a7f3d0"; // Light Green: Seed (Dynamic, feeds a QSE group)
             } else if (std::ranges::contains(m_dynamic_species_indices, i)) {
                 fillcolor = "#dcfce7"; // Pale Green: Core Dynamic
@@ -318,7 +455,7 @@ namespace gridfire {
         // --- Layout and Ranking ---
         // Enforce a top-down layout based on mass number.
         dotFile << "    // --- Layout using Ranks ---\n";
-        for (const auto& [mass, species_list] : species_by_mass) {
+        for (const auto &species_list: species_by_mass | std::views::values) {
             dotFile << "    { rank=same; ";
             for (const auto& name : species_list) {
                 dotFile << "\"" << name << "\"; ";
@@ -332,7 +469,7 @@ namespace gridfire {
         for (const auto& [mass, species_list] : species_by_mass) {
             // Find the next largest mass in the species list
             int minLargestMass = std::numeric_limits<int>::max();
-            for (const auto& [next_mass, _] : species_by_mass) {
+            for (const auto &next_mass: species_by_mass | std::views::keys) {
                 if (next_mass > mass && next_mass < minLargestMass) {
                     minLargestMass = next_mass;
                 }
@@ -387,7 +524,7 @@ namespace gridfire {
                 dotFile << "            " << ss.str() << "    [label=\"" << all_species[species_idx].name() << "\"];\n";
                 algebraic_node_ids.push_back(ss.str());
             }
-            // Make invisible edges between algebraic indicies to keep them in top down order
+            // Make invisible edges between algebraic indices to keep them in top-down order
             for (size_t i = 0; i < algebraic_node_ids.size() - 1; ++i) {
                 dotFile << "            " << algebraic_node_ids[i] << " -> " << algebraic_node_ids[i + 1] << " [style=invis];\n";
             }
@@ -511,21 +648,15 @@ namespace gridfire {
         }
 
         composition.finalize(true);
+
         return composition;
     }
 
     fourdst::composition::Composition MultiscalePartitioningEngineView::equilibrateNetwork(
         const NetIn &netIn
     ) {
-        std::vector<double> Y(m_baseEngine.getNetworkSpecies().size(), 0.0);
+        const PrimingReport primingReport = m_baseEngine.primeEngine(netIn);
-        for (const auto& [symbol, entry]: netIn.composition ) {
+        const std::vector<double> Y = packCompositionToVector(primingReport.primedComposition, m_baseEngine);
-            if (!m_baseEngine.involvesSpecies(entry.isotope())) {
-                LOG_WARNING(m_logger, "Species {} is not part of the network. This is likely due to the base network not being deep enough. For rare species this may not be an issue. Skipping...", entry.isotope().name());
-                continue; // Skip species not in the network
-            }
-            const size_t species_index = m_baseEngine.getSpeciesIndex(entry.isotope());
-            Y[species_index] = netIn.composition.getMolarAbundance(symbol);
-        }
 
         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
@@ -1128,6 +1259,7 @@ namespace gridfire {
         return candidate_groups;
     }
 
+
     int MultiscalePartitioningEngineView::EigenFunctor::operator()(const InputType &v_qse, OutputType &f_qse) const {
         s_operator_parens_called++;
         std::vector<double> y_trial = m_Y_full_initial;
@@ -1137,7 +1269,7 @@ namespace gridfire {
             y_trial[m_qse_solve_indices[i]] = y_qse(i);
         }
 
-        const auto [dydt, energy] = m_view->calculateRHSAndEnergy(y_trial, m_T9, m_rho);
+        const auto [dydt, energy] = m_view->getBaseEngine().calculateRHSAndEnergy(y_trial, m_T9, m_rho);
         f_qse.resize(m_qse_solve_indices.size());
         for (size_t i = 0; i < m_qse_solve_indices.size(); ++i) {
             f_qse(i) = dydt[m_qse_solve_indices[i]];
@@ -1155,13 +1287,13 @@ namespace gridfire {
         }
 
         // TODO: Think about if the jacobian matrix should be mutable so that generateJacobianMatrix can be const
-        m_view->generateJacobianMatrix(y_trial, m_T9, m_rho);
+        m_view->getBaseEngine().generateJacobianMatrix(y_trial, m_T9, m_rho);
 
         // TODO: Think very carefully about the indices here.
         J_qse.resize(m_qse_solve_indices.size(), m_qse_solve_indices.size());
         for (size_t i = 0; i < m_qse_solve_indices.size(); ++i) {
             for (size_t j = 0; j < m_qse_solve_indices.size(); ++j) {
-                J_qse(i, j) = m_view->getJacobianMatrixEntry(
+                J_qse(i, j) = m_view->getBaseEngine().getJacobianMatrixEntry(
                     m_qse_solve_indices[i],
                     m_qse_solve_indices[j]
                 );
@@ -1175,4 +1307,45 @@ namespace gridfire {
         }
         return 0; // Success
     }
+
+
+    QSECacheKey::QSECacheKey(
+        const double T9,
+        const double rho,
+        const std::vector<double> &Y
+    ) :
+    m_T9(T9),
+    m_rho(rho),
+    m_Y(Y) {
+        m_hash = hash();
+    }
+
+    size_t QSECacheKey::hash() const {
+        std::size_t seed = 0;
+
+        hash_combine(seed, m_Y.size());
+
+        hash_combine(seed, bin(m_T9, m_cacheConfig.T9_tol));
+        hash_combine(seed, bin(m_rho, m_cacheConfig.rho_tol));
+
+        for (double Yi : m_Y) {
+            if (Yi < 0.0 && std::abs(Yi) < 1e-20) {
+                Yi = 0.0; // Avoid negative abundances
+            } else if (Yi < 0.0 && std::abs(Yi) >= 1e-20) {
+                throw std::invalid_argument("Yi should be positive for valid hashing (expected Yi > 0, received " + std::to_string(Yi) + ")");
+            }
+            hash_combine(seed, bin(Yi, m_cacheConfig.Yi_tol));
+        }
+
+        return seed;
+
+    }
+
+    long QSECacheKey::bin(const double value, const double tol) {
+        return static_cast<long>(std::floor(value / tol));
+    }
+
+    bool QSECacheKey::operator==(const QSECacheKey &other) const {
+        return m_hash == other.m_hash;
+    }
 }

src/network/lib/solver/solver.cpp

@@ -1,6 +1,7 @@
 #include "gridfire/solver/solver.h"
 #include "gridfire/engine/engine_graph.h"
 #include "gridfire/network.h"
+#include "gridfire/exceptions/error_engine.h"
 
 
 #include "gridfire/utils/logging.h"
@@ -20,8 +21,22 @@
 #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;
 
@@ -507,37 +522,47 @@ namespace gridfire::solver {
 
 
         const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9)
-        const unsigned long numSpecies = m_engine.getNetworkSpecies().size();
 
         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();
+        ublas::vector<double> Y(numSpecies + 1);
+
         RHSFunctor rhsFunctor(m_engine, T9, netIn.density);
         JacobianFunctor jacobianFunctor(m_engine, T9, netIn.density);
 
-        ublas::vector<double> Y(numSpecies + 1);
+
+        // 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) {
             const auto& species = m_engine.getNetworkSpecies()[i];
-            try {
+            if (!equilibratedComposition.contains(species)) {
-                Y(i) = netIn.composition.getMolarAbundance(std::string(species.name()));
+                LOG_DEBUG(m_logger, "Species '{}' not found in composition. Setting abundance to 1e-99.", species.name());
-            } catch (const std::runtime_error) {
+                Y(i) = 1e-99;
-                LOG_DEBUG(m_logger, "Species '{}' not found in composition. Setting abundance to 0.0.", species.name());
+                continue;
-                Y(i) = 0.0;
             }
+            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(
             stepper,
             std::make_pair(rhsFunctor, jacobianFunctor),
             Y,
             0.0,
             netIn.tMax,
-            netIn.dt0
+            netIn.dt0,
+            observer
         );
 
         std::vector<double> finalMassFractions(numSpecies);
@@ -570,13 +595,63 @@ namespace gridfire::solver {
     void DirectNetworkSolver::RHSFunctor::operator()(
         const boost::numeric::ublas::vector<double> &Y,
         boost::numeric::ublas::vector<double> &dYdt,
-        double t
+        const double t
     ) const {
         const std::vector<double> y(Y.begin(), m_numSpecies + Y.begin());
-        auto [dydt, eps] = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho);
+        // const auto timescales = m_engine.getSpeciesTimescales(y, m_T9, m_rho);
+        StepDerivatives<double> deriv;
+
+        // TODO: Refactor this to use std::expected
+        try {
+            deriv = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho);
+        } catch (const exceptions::StaleEngineError& e) { // If the engine reports that it is stale, we need to update it
+            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);
+            for (size_t i = 0; i < y.size(); ++i) {
+                double Xi = y[i] * m_engine.getNetworkSpecies()[i].mass(); // Convert from molar abundance to mass fraction
+                if (Xi < 0 && std::abs(Xi) < 1e-20) {
+                    Xi = 0.0; // Avoid negative mass fractions
+                }
+                X[i] = Xi;
+            }
+
+            compositionTemp.setMassFraction(m_engine.getNetworkSpecies(), X);
+            compositionTemp.finalize(true);
+
+            netInTemp.composition = std::move(compositionTemp);
+            netInTemp.temperature = m_T9 * 1e9; // Convert T9 back to Kelvin
+            netInTemp.density = m_rho; // Density in g/cm^3
+
+            m_engine.update(netInTemp);
+
+            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!");
+        }
+
         dYdt.resize(m_numSpecies + 1);
-        std::ranges::copy(dydt, dYdt.begin());
+        // for (size_t i = 0; i < m_numSpecies; ++i) {
-        dYdt(m_numSpecies) = eps;
+        //     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()(
@@ -590,6 +665,7 @@ namespace gridfire::solver {
         for (int i = 0; i < m_numSpecies; ++i) {
             for (int j = 0; j < m_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';
             }
         }
     }

tests/graphnet_sandbox/main.cpp

@@ -51,7 +51,7 @@ void quill_terminate_handler()
 int main() {
     g_previousHandler = std::set_terminate(quill_terminate_handler);
     quill::Logger* logger = fourdst::logging::LogManager::getInstance().getLogger("log");
-    logger->set_log_level(quill::LogLevel::TraceL1);
+    logger->set_log_level(quill::LogLevel::Info);
     LOG_DEBUG(logger, "Starting Adaptive Engine View Example...");
 
     using namespace gridfire;
@@ -74,24 +74,22 @@ int main() {
     netIn.temperature = 1.5e7;
     netIn.density = 1.5e3;
     netIn.energy = 0;
-    netIn.tMax = 3e17;
+    netIn.tMax = 3.546e17;
     netIn.dt0 = 1e-12;
 
     GraphEngine ReaclibEngine(composition, partitionFunction, NetworkBuildDepth::SecondOrder);
-    ReaclibEngine.exportToDot("GraphEngine.dot");
+    ReaclibEngine.setPrecomputation(true);
-    auto primedReport = ReaclibEngine.primeEngine(netIn);
+    ReaclibEngine.setUseReverseReactions(false);
-    if (!primedReport.success) {
+    MultiscalePartitioningEngineView partitioningView(ReaclibEngine);
-        LOG_CRITICAL(logger, "Failed to prime the network!");
+    AdaptiveEngineView adaptiveView(partitioningView);
-        return 1;
+
+    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";
     }
 
-    NetIn primedNetIn = netIn;
-    primedNetIn.composition = primedReport.primedComposition;
-
-    std::cout << primedReport.primedComposition << std::endl;
-
-    MultiscalePartitioningEngineView partitioningView(ReaclibEngine);
-    fourdst::composition::Composition qseComp = partitioningView.equilibrateNetwork(primedNetIn);
-    std::cout << qseComp.getMolarAbundance("H-2") / qseComp.getMolarAbundance("H-1") << std::endl;
-
 }