feat(MultiscalePartitioningEngineView): added *much* more robust qse group identifiction and solving

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

@@ -123,7 +123,7 @@ namespace gridfire {
         /**
          * @brief Generate the Jacobian matrix for the current state.
          *
-         * @param Y Vector of current abundances.
+         * @param Y_dynamic Vector of current abundances.
          * @param T9 Temperature in units of 10^9 K.
          * @param rho Density in g/cm^3.
          *
@@ -131,7 +131,7 @@ namespace gridfire {
          * for the current state. The matrix can then be accessed via getJacobianMatrixEntry().
          */
         virtual void generateJacobianMatrix(
-            const std::vector<double>& Y,
+            const std::vector<double>& Y_dynamic,
             double T9, double rho
         ) = 0;
 
@@ -265,5 +265,9 @@ namespace gridfire {
          * @endcode
          */
         [[nodiscard]] virtual screening::ScreeningType getScreeningModel() const = 0;
+
+        [[nodiscard]] virtual int getSpeciesIndex(const fourdst::atomic::Species &species) const = 0;
+
+        [[nodiscard]] virtual std::vector<double> mapNetInToMolarAbundanceVector(const NetIn &netIn) const = 0;
     };
 }

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

@@ -17,10 +17,12 @@
 #include <unordered_map>
 #include <vector>
 #include <memory>
+#include <ranges>
 
 #include <boost/numeric/ublas/matrix_sparse.hpp>
 
 #include "cppad/cppad.hpp"
+#include "quill/LogMacros.h"
 
 // PERF: The function getNetReactionStoichiometry returns a map of species to their stoichiometric coefficients for a given reaction.
 //       this makes extra copies of the species, which is not ideal and could be optimized further.
@@ -139,7 +141,7 @@ namespace gridfire {
         /**
          * @brief Generates the Jacobian matrix for the current state.
          *
-         * @param Y Vector of current abundances.
+         * @param Y_dynamic Vector of current abundances.
          * @param T9 Temperature in units of 10^9 K.
          * @param rho Density in g/cm^3.
          *
@@ -150,7 +152,7 @@ namespace gridfire {
          * @see getJacobianMatrixEntry()
          */
         void generateJacobianMatrix(
-            const std::vector<double>& Y,
+            const std::vector<double>& Y_dynamic,
             const double T9,
             const double rho
         ) override;
@@ -317,18 +319,48 @@ namespace gridfire {
 
         [[nodiscard]] double calculateReverseRate(
             const reaction::Reaction &reaction,
-            double T9,
-            double expFactor
+            double T9
         ) const;
 
+        double calculateReverseRateTwoBody(
+            const reaction::Reaction &reaction,
+            const double T9,
+            const double forwardRate,
+            const double expFactor
+        ) const;
+
+        double calculateReverseRateTwoBodyDerivative(
+            const reaction::Reaction &reaction,
+            const double T9,
+            const double reverseRate
+        ) const;
+
+        bool isUsingReverseReactions() const;
+
+        void setUseReverseReactions(bool useReverse);
+
+        int getSpeciesIndex(
+            const fourdst::atomic::Species& species
+        ) const override;
+
+        std::vector<double> mapNetInToMolarAbundanceVector(const NetIn &netIn) const override;
+
+
+
     private:
         struct PrecomputedReaction {
+            // Forward cacheing
             size_t reaction_index;
             std::vector<size_t> unique_reactant_indices;
             std::vector<int> reactant_powers;
             double symmetry_factor;
             std::vector<size_t> affected_species_indices;
             std::vector<int> stoichiometric_coefficients;
+
+            // Reverse cacheing
+            std::vector<size_t> unique_product_indices; ///< Unique product indices for reverse reactions.
+            std::vector<int> product_powers; ///< Powers of each unique product in the reverse reaction.
+            double reverse_symmetry_factor; ///< Symmetry factor for reverse reactions.
         };
 
         struct constants {
@@ -337,7 +369,47 @@ namespace gridfire {
             const double c = Constants::getInstance().get("c").value; ///< Speed of light in cm/s.
             const double kB = Constants::getInstance().get("kB").value; ///< Boltzmann constant in erg/K.
         };
+    private:
+        class AtomicReverseRate final : public CppAD::atomic_base<double> {
+        public:
+            AtomicReverseRate(
+                const reaction::Reaction& reaction,
+                const GraphEngine& engine
+            ):
+            atomic_base<double>("AtomicReverseRate"),
+            m_reaction(reaction),
+            m_engine(engine) {}
 
+            bool forward(
+                size_t p,
+                size_t q,
+                const CppAD::vector<bool>& vx,
+                CppAD::vector<bool>& vy,
+                const CppAD::vector<double>& tx,
+                CppAD::vector<double>& ty
+            ) override;
+            bool reverse(
+                size_t q,
+                const CppAD::vector<double>& tx,
+                const CppAD::vector<double>& ty,
+                CppAD::vector<double>& px,
+                const CppAD::vector<double>& py
+            ) override;
+            bool for_sparse_jac(
+                size_t q,
+                const CppAD::vector<std::set<size_t>>&r,
+                CppAD::vector<std::set<size_t>>& s
+            ) override;
+            bool rev_sparse_jac(
+                size_t q,
+                const CppAD::vector<std::set<size_t>>&rt,
+                CppAD::vector<std::set<size_t>>& st
+            ) override;
+
+        private:
+            const reaction::Reaction& m_reaction;
+            const GraphEngine& m_engine;
+        };
     private:
         Config& m_config = Config::getInstance();
         quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
@@ -355,12 +427,15 @@ namespace gridfire {
         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.
+        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.
         std::unique_ptr<screening::ScreeningModel> m_screeningModel = screening::selectScreeningModel(m_screeningType);
 
         bool m_usePrecomputation = true; ///< Flag to enable or disable using precomputed reactions for efficiency. Mathematically, this should not change the results. Generally end users should not need to change this.
 
+        bool m_useReverseReactions = true; ///< Flag to enable or disable reverse reactions. If false, only forward reactions are considered.
+
         std::vector<PrecomputedReaction> m_precomputedReactions; ///< Precomputed reactions for efficiency.
         std::unique_ptr<partition::PartitionFunction> m_partitionFunction; ///< Partition function for the network.
 
@@ -418,8 +493,11 @@ namespace gridfire {
          */
         void recordADTape();
 
+        void collectAtomicReverseRateAtomicBases();
+
         void precomputeNetwork();
 
+
         /**
          * @brief Validates mass and charge conservation across all reactions.
          *
@@ -449,24 +527,12 @@ namespace gridfire {
         );
 
 
-        double calculateReverseRateTwoBody(
-            const reaction::Reaction &reaction,
-            const double T9,
-            const double forwardRate,
-            const double expFactor
-        ) const;
-
-        double GraphEngine::calculateReverseRateTwoBodyDerivative(
-            const reaction::Reaction &reaction,
-            const double T9,
-            const double reverseRate
-        ) const;
 
         [[nodiscard]] StepDerivatives<double> calculateAllDerivativesUsingPrecomputation(
             const std::vector<double> &Y_in,
             const std::vector<double>& bare_rates,
-            double T9,
-            double rho
+            const std::vector<double> &bare_reverse_rates,
+            double T9, double rho
         ) const;
 
         /**
@@ -490,6 +556,16 @@ namespace gridfire {
             const T rho
         ) const;
 
+        template<IsArithmeticOrAD T>
+        T calculateReverseMolarReactionFlow(
+            T T9,
+            T rho,
+            std::vector<T> screeningFactors,
+            std::vector<T> Y,
+            size_t reactionIndex,
+            const reaction::LogicalReaction &reaction
+        ) const;
+
         /**
          * @brief Calculates all derivatives (dY/dt) and the energy generation rate.
          *
@@ -547,6 +623,74 @@ namespace gridfire {
     };
 
 
+
+    template <IsArithmeticOrAD T>
+    T GraphEngine::calculateReverseMolarReactionFlow(
+        T T9,
+        T rho,
+        std::vector<T> screeningFactors,
+        std::vector<T> Y,
+        size_t reactionIndex,
+        const reaction::LogicalReaction &reaction
+    ) const {
+        if (!m_useReverseReactions) {
+            return static_cast<T>(0.0); // If reverse reactions are not used, return zero
+        }
+        T reverseMolarFlow = static_cast<T>(0.0);
+
+        if (reaction.qValue() != 0.0) {
+            T reverseRateConstant = static_cast<T>(0.0);
+            if constexpr (std::is_same_v<T, ADDouble>) { // Check if T is an AD type at compile time
+                const auto& atomic_func_ptr = m_atomicReverseRates[reactionIndex];
+                if (atomic_func_ptr != nullptr) {
+                    // A. Instantiate the atomic operator for the specific reaction
+                    // B. Marshal the input vector
+                    std::vector<T> ax = { T9 };
+
+                    std::vector<T> ay(1);
+                    (*atomic_func_ptr)(ax, ay);
+                    reverseRateConstant = static_cast<T>(ay[0]);
+                } else {
+                    return reverseMolarFlow; // If no atomic function is available, return zero
+                }
+            } else {
+                // A,B If not calling with an AD type, calculate the reverse rate directly
+                reverseRateConstant = calculateReverseRate(reaction, T9);
+            }
+
+            // C. Get product multiplicities
+            std::unordered_map<fourdst::atomic::Species, int> productCounts;
+            for (const auto& product : reaction.products()) {
+                productCounts[product]++;
+            }
+
+            // D. Calculate the symmetry factor
+            T reverseSymmetryFactor = static_cast<T>(1.0);
+            for (const auto &count: productCounts | std::views::values) {
+                reverseSymmetryFactor /= static_cast<T>(std::tgamma(static_cast<double>(count + 1))); // Gamma function for factorial
+            }
+
+            // E. Calculate the abundance term
+            T productAbundanceTerm = static_cast<T>(1.0);
+            for (const auto& [species, count] : productCounts) {
+                const int speciesIndex = m_speciesToIndexMap.at(species);
+                productAbundanceTerm *= CppAD::pow(Y[speciesIndex], count);
+            }
+
+            // F. Determine the power for the density term
+            const size_t num_products = reaction.products().size();
+            const T rho_power = CppAD::pow(rho, static_cast<T>(num_products > 1 ? num_products - 1 : 0)); // Density raised to the power of (N-1) for N products
+
+            // G. Assemble the reverse molar flow rate
+            reverseMolarFlow = screeningFactors[reactionIndex] *
+                               reverseRateConstant *
+                               productAbundanceTerm *
+                               reverseSymmetryFactor *
+                               rho_power;
+        }
+        return reverseMolarFlow;
+    }
+
     template<IsArithmeticOrAD T>
     StepDerivatives<T> GraphEngine::calculateAllDerivatives(
         const std::vector<T> &Y_in, T T9, T rho) const {
@@ -594,13 +738,39 @@ namespace gridfire {
         for (size_t reactionIndex = 0; reactionIndex < m_reactions.size(); ++reactionIndex) {
             const auto& reaction = m_reactions[reactionIndex];
 
-            // 1. Calculate reaction rate
-            const T molarReactionFlow = screeningFactors[reactionIndex] * calculateMolarReactionFlow<T>(reaction, Y, T9, rho);
+            // 1. Calculate forward reaction rate
+            const T forwardMolarReactionFlow = screeningFactors[reactionIndex] *
+                calculateMolarReactionFlow<T>(reaction, Y, T9, rho);
 
-            // 2. Use the rate to update all relevant species derivatives (dY/dt)
+            // 2. Calculate reverse reaction rate
+            T reverseMolarFlow = calculateReverseMolarReactionFlow<T>(
+                T9,
+                rho,
+                screeningFactors,
+                Y,
+                reactionIndex,
+                reaction
+            );
+
+            const T molarReactionFlow = forwardMolarReactionFlow - reverseMolarFlow; // Net molar reaction flow
+            std::stringstream ss;
+            ss << "Forward: " << forwardMolarReactionFlow
+               << ", Reverse: " << reverseMolarFlow
+               << ", Net: " << molarReactionFlow;
+            LOG_DEBUG(
+                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 T nu_ij = static_cast<T>(m_stoichiometryMatrix(speciesIndex, reactionIndex));
-                result.dydt[speciesIndex] += threshold_flag * nu_ij * molarReactionFlow / rho;
+                result.dydt[speciesIndex] += threshold_flag * nu_ij * molarReactionFlow;
             }
         }
 
@@ -644,6 +814,7 @@ namespace gridfire {
         for (const auto& reactant : reaction.reactants()) {
             reactant_counts[std::string(reactant.name())]++;
         }
+        const int totalReactants = static_cast<int>(reaction.reactants().size());
 
         // --- Accumulator for the molar concentration ---
         auto molar_concentration_product = static_cast<T>(1.0);
@@ -657,56 +828,23 @@ namespace gridfire {
             const T Yi = Y[species_index];
 
             // --- Check if the species abundance is below the threshold where we ignore reactions ---
-            threshold_flag *= CppAD::CondExpLt(Yi, Y_threshold, zero, one);
-
-            // --- Convert from molar abundance to molar concentration ---
-            T molar_concentration = Yi * rho;
+            // threshold_flag *= CppAD::CondExpLt(Yi, Y_threshold, zero, one);
 
             // --- If count is > 1 , we need to raise the molar concentration to the power of count since there are really count bodies in that reaction ---
-            molar_concentration_product *= CppAD::pow(molar_concentration, static_cast<T>(count)); // ni^count
+            molar_concentration_product *= CppAD::pow(Yi, static_cast<T>(count)); // ni^count
 
             // --- Apply factorial correction for identical reactions ---
             if (count > 1) {
                 molar_concentration_product /= static_cast<T>(std::tgamma(static_cast<double>(count + 1))); // Gamma function for factorial
             }
         }
-        // --- Final reaction flow calculation [mol][s^-1][cm^-3] ---
+        // --- Final reaction flow calculation [mol][s^-1][g^-1] ---
         // Note: If the threshold flag ever gets set to zero this will return zero.
         //       This will result basically in multiple branches being written to the AD tape, which will make
         //       the tape more expensive to record, but it will also mean that we only need to record it once for
         //       the entire network.
-        return molar_concentration_product * k_reaction * threshold_flag;
+        const T densityTerm = CppAD::pow(rho, totalReactants > 1 ? static_cast<T>(totalReactants - 1) : zero); // Density raised to the power of (N-1) for N reactants
+        return molar_concentration_product * k_reaction * threshold_flag * densityTerm;
     }
 
-    class AtomicReverseRate final : public CppAD::atomic_base<double> {
-    public:
-        AtomicReverseRate(
-            const reaction::Reaction& reaction,
-            const GraphEngine& engine
-        ):
-        atomic_base<double>("AtomicReverseRate"),
-        m_reaction(reaction),
-        m_engine(engine) {}
-
-        bool forward(
-            size_t p,
-            size_t q,
-            const CppAD::vector<bool>& vx,
-            CppAD::vector<bool>& vy,
-            const CppAD::vector<double>& tx,
-            CppAD::vector<double>& ty
-        ) override;
-        bool reverse(
-            size_t id,
-            size_t an,
-            const CppAD::vector<double>& tx,
-            const CppAD::vector<double>& ty,
-            CppAD::vector<double>& px,
-            const CppAD::vector<double>& py
-        );
-    private:
-        const double m_kB = Constants::getInstance().get("k_b").value; ///< Boltzmann constant in erg/K.
-        const reaction::Reaction& m_reaction;
-        const GraphEngine& m_engine;
-    };
 };

src/network/include/gridfire/engine/procedures/priming.h

@@ -0,0 +1,31 @@
+#pragma once
+
+#include "gridfire/engine/engine_abstract.h"
+#include "gridfire/network.h"
+
+#include "fourdst/composition/composition.h"
+#include "fourdst/composition/atomicSpecies.h"
+
+
+namespace gridfire {
+    fourdst::composition::Composition primeNetwork(
+        const NetIn&,
+        DynamicEngine& engine
+    );
+
+    double calculateDestructionRateConstant(
+        const DynamicEngine& engine,
+        const fourdst::atomic::Species& species,
+        const std::vector<double>& Y,
+        double T9,
+        double rho
+    );
+
+    double calculateCreationRate(
+        const DynamicEngine& engine,
+        const fourdst::atomic::Species& species,
+        const std::vector<double>& Y,
+        double T9,
+        double rho
+    );
+}

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

@@ -106,7 +106,7 @@ namespace gridfire {
         /**
          * @brief Generates the Jacobian matrix for the active species.
          *
-         * @param Y_culled A vector of abundances for the active species.
+         * @param Y_dynamic A vector of abundances for the active species.
          * @param T9 The temperature in units of 10^9 K.
          * @param rho The density in g/cm^3.
          *
@@ -117,7 +117,7 @@ namespace gridfire {
          * @see AdaptiveEngineView::update()
          */
         void generateJacobianMatrix(
-            const std::vector<double> &Y_culled,
+            const std::vector<double> &Y_dynamic,
             const double T9,
             const double rho
         ) override;
@@ -256,6 +256,10 @@ namespace gridfire {
          * @endcode
          */
         [[nodiscard]] screening::ScreeningType getScreeningModel() const override;
+
+        [[nodiscard]] int getSpeciesIndex(const fourdst::atomic::Species &species) const override;
+
+        [[nodiscard]] std::vector<double> mapNetInToMolarAbundanceVector(const NetIn &netIn) const override;
     private:
         using Config = fourdst::config::Config;
         using LogManager = fourdst::logging::LogManager;

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

@@ -13,56 +13,9 @@
 #include <string>
 
 namespace gridfire{
-    /**
-     * @class FileDefinedEngineView
-     * @brief An engine view that uses a user-defined reaction network from a file.
-     *
-     * This class implements an EngineView that restricts the reaction network to a specific set of
-     * reactions defined in an external file. It acts as a filter or a view on a larger, more
-     * comprehensive base engine. The file provides a list of reaction identifiers, and this view
-     * will only consider those reactions and the species involved in them.
-     *
-     * This is useful for focusing on a specific sub-network for analysis, debugging, or performance
-     * reasons, without modifying the underlying full network.
-     *
-     * The view maintains mappings between the indices of its active (defined) species and reactions
-     * and the corresponding indices in the full network of the base engine. All calculations are
-     * delegated to the base engine after mapping the inputs from the view's context to the full
-     * network context, and the results are mapped back.
-     *
-     * @implements DynamicEngine
-     * @implements EngineView<DynamicEngine>
-     */
-    class FileDefinedEngineView final: public DynamicEngine, public EngineView<DynamicEngine> {
+    class DefinedEngineView : public DynamicEngine, public EngineView<DynamicEngine> {
     public:
-        /**
-         * @brief Constructs a FileDefinedEngineView.
-         *
-         * @param baseEngine The underlying DynamicEngine to which this view delegates calculations.
-         * @param fileName The path to the file that defines the reaction network for this view.
-         * @param parser A reference to a parser object capable of parsing the network file.
-         *
-         * @par Usage Example:
-         * @code
-         * MyParser parser;
-         * DynamicEngine baseEngine(...);
-         * FileDefinedEngineView view(baseEngine, "my_network.net", parser);
-         * @endcode
-         *
-         * @post The view is initialized with the reactions and species from the specified file.
-         * @throws std::runtime_error If a reaction from the file is not found in the base engine.
-         */
-        explicit FileDefinedEngineView(
-            DynamicEngine& baseEngine,
-            const std::string& fileName,
-            const io::NetworkFileParser& parser
-        );
-
-        // --- EngineView Interface ---
-        /**
-         * @brief Gets the base engine.
-         * @return A const reference to the base engine.
-         */
+        DefinedEngineView(const std::vector<std::string>& peNames, DynamicEngine& baseEngine);
         const DynamicEngine& getBaseEngine() const override;
 
         // --- Engine Interface ---
@@ -92,14 +45,14 @@ namespace gridfire{
         /**
          * @brief Generates the Jacobian matrix for the active species.
          *
-         * @param Y_defined A vector of abundances for the active species.
+         * @param Y_dynamic A vector of abundances for the active species.
          * @param T9 The temperature in units of 10^9 K.
          * @param rho The density in g/cm^3.
          *
          * @throws std::runtime_error If the view is stale.
          */
         void generateJacobianMatrix(
-            const std::vector<double>& Y_defined,
+            const std::vector<double>& Y_dynamic,
             const double T9,
             const double rho
         ) override;
@@ -191,21 +144,6 @@ namespace gridfire{
          */
         void update(const NetIn &netIn) override;
 
-        /**
-         * @brief Sets a new network file to define the active reactions.
-         *
-         * @param fileName The path to the new network definition file.
-         *
-         * @par Usage Example:
-         * @code
-         * view.setNetworkFile("another_network.net");
-         * view.update(netIn); // Must be called before using the view again
-         * @endcode
-         *
-         * @post The view is marked as stale. `update()` must be called before further use.
-         */
-        void setNetworkFile(const std::string& fileName);
-
         /**
          * @brief Sets the screening model for the base engine.
          *
@@ -219,21 +157,15 @@ namespace gridfire{
          * @return The current screening model type.
          */
         [[nodiscard]] screening::ScreeningType getScreeningModel() const override;
-    private:
-        using Config = fourdst::config::Config;
-        using LogManager = fourdst::logging::LogManager;
-        /** @brief A reference to the singleton Config instance. */
-        Config& m_config = Config::getInstance();
-        /** @brief A pointer to the logger instance. */
-        quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
 
-        /** @brief The underlying engine to which this view delegates calculations. */
+        [[nodiscard]] int getSpeciesIndex(const fourdst::atomic::Species &species) const override;
+
+        [[nodiscard]] std::vector<double> mapNetInToMolarAbundanceVector(const NetIn &netIn) const override;
+    protected:
+        bool m_isStale = true;
         DynamicEngine& m_baseEngine;
-        ///< Name of the file defining the reaction set considered by the engine view.
-        std::string m_fileName;
-        ///< Parser for the network file.
-        const io::NetworkFileParser& m_parser;
-
+    private:
+        quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); ///< Logger instance for trace and debug information.
         ///< Active species in the defined engine.
         std::vector<fourdst::atomic::Species> m_activeSpecies;
         ///< Active reactions in the defined engine.
@@ -243,30 +175,7 @@ namespace gridfire{
         std::vector<size_t> m_speciesIndexMap;
         ///< Maps indices of active reactions to indices in the full network.
         std::vector<size_t> m_reactionIndexMap;
-
-        /** @brief A flag indicating whether the view is stale and needs to be updated. */
-        bool m_isStale = true;
-
     private:
-        /**
-         * @brief Builds the active species and reaction sets from a file.
-         *
-         * This method uses the provided parser to read reaction names from the given file.
-         * It then finds these reactions in the base engine's full network and populates
-         * the `m_activeReactions` and `m_activeSpecies` members. Finally, it constructs
-         * the index maps for the active sets.
-         *
-         * @param fileName The path to the network definition file.
-         *
-         * @post
-         * - `m_activeReactions` and `m_activeSpecies` are populated.
-         * - `m_speciesIndexMap` and `m_reactionIndexMap` are constructed.
-         * - `m_isStale` is set to false.
-         *
-         * @throws std::runtime_error If a reaction from the file is not found in the base engine.
-         */
-        void buildFromFile(const std::string& fileName);
-
         /**
          * @brief Constructs the species index map.
          *
@@ -329,11 +238,25 @@ namespace gridfire{
          */
         size_t mapViewToFullReactionIndex(size_t definedReactionIndex) const;
 
-        /**
-         * @brief Validates that the FileDefinedEngineView is not stale.
-         *
-         * @throws std::runtime_error If the view is stale (i.e., `update()` has not been called after the view was made stale).
-         */
         void validateNetworkState() const;
     };
+
+    class FileDefinedEngineView final: public DefinedEngineView {
+    public:
+        explicit FileDefinedEngineView(
+            DynamicEngine& baseEngine,
+            const std::string& fileName,
+            const io::NetworkFileParser& parser
+        );
+        std::string getNetworkFile() const { return m_fileName; }
+        const io::NetworkFileParser& getParser() const { return m_parser; }
+    private:
+        using Config = fourdst::config::Config;
+        using LogManager = fourdst::logging::LogManager;
+        Config& m_config = Config::getInstance();
+        quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
+        std::string m_fileName;
+        ///< Parser for the network file.
+        const io::NetworkFileParser& m_parser;
+    };
 }

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

@@ -0,0 +1,188 @@
+#pragma once
+
+#include "gridfire/engine/engine_abstract.h"
+#include "gridfire/engine/views/engine_view_abstract.h"
+#include "gridfire/engine/engine_graph.h"
+
+#include "Eigen/Dense"
+#include "unsupported/Eigen/NonLinearOptimization"
+
+namespace gridfire {
+    class MultiscalePartitioningEngineView final: public DynamicEngine, public EngineView<DynamicEngine> {
+    public:
+        explicit MultiscalePartitioningEngineView(GraphEngine& baseEngine);
+
+        [[nodiscard]] const std::vector<fourdst::atomic::Species> & getNetworkSpecies() const override;
+
+        [[nodiscard]] StepDerivatives<double> calculateRHSAndEnergy(
+            const std::vector<double> &Y,
+            double T9,
+            double rho
+        ) const override;
+
+        void generateJacobianMatrix(
+            const std::vector<double> &Y_dynamic,
+            double T9,
+            double rho
+        ) override;
+
+        [[nodiscard]] double getJacobianMatrixEntry(
+            int i,
+            int j
+        ) const override;
+
+        void generateStoichiometryMatrix() override;
+
+        [[nodiscard]] int getStoichiometryMatrixEntry(
+            int speciesIndex,
+            int reactionIndex
+        ) const override;
+
+        [[nodiscard]] double calculateMolarReactionFlow(
+            const reaction::Reaction &reaction,
+            const std::vector<double> &Y,
+            double T9,
+            double rho
+        ) const override;
+
+        [[nodiscard]] const reaction::LogicalReactionSet & getNetworkReactions() const override;
+
+        [[nodiscard]] std::unordered_map<fourdst::atomic::Species, double> getSpeciesTimescales(
+            const std::vector<double> &Y,
+            double T9,
+            double rho
+        ) const override;
+
+        void update(
+            const NetIn &netIn
+        ) override;
+
+        void setScreeningModel(
+            screening::ScreeningType model
+        ) override;
+
+        [[nodiscard]] screening::ScreeningType getScreeningModel() const override;
+
+        const DynamicEngine & getBaseEngine() const override;
+
+        void partitionNetwork(
+            const std::vector<double>& Y,
+            double T9,
+            double rho,
+            double dt_control
+        );
+
+        void partitionNetwork(
+            const NetIn& netIn,
+            double dt_control
+        );
+
+        void exportToDot(const std::string& filename) const;
+
+        [[nodiscard]] int getSpeciesIndex(const fourdst::atomic::Species &species) const override;
+
+        [[nodiscard]] std::vector<double> mapNetInToMolarAbundanceVector(const NetIn &netIn) const override;
+
+        std::vector<fourdst::atomic::Species> getFastSpecies() const;
+        const std::vector<fourdst::atomic::Species>& getDynamicSpecies() const;
+
+        void equilibrateNetwork(
+            const std::vector<double>& Y,
+            double T9,
+            double rho,
+            double dt_control
+        );
+
+        void equilibrateNetwork(
+            const NetIn& netIn,
+            const double dt_control
+        );
+
+
+    private:
+        struct QSEGroup {
+            std::vector<size_t> species_indices; ///< Indices of all species in this group.
+            size_t seed_nucleus_index;           ///< Index of the one species that will be evolved dynamically.
+            bool is_in_equilibrium = false;      ///< Flag set by flux analysis.
+        };
+
+        struct EigenFunctor {
+            using InputType = Eigen::Matrix<double, Eigen::Dynamic, 1>;
+            using OutputType = Eigen::Matrix<double, Eigen::Dynamic, 1>;
+            using JacobianType = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>;
+            enum {
+                InputsAtCompileTime = Eigen::Dynamic,
+                ValuesAtCompileTime = Eigen::Dynamic
+            };
+
+            MultiscalePartitioningEngineView* m_view;
+            const std::vector<size_t>& m_qse_solve_indices;
+            const std::vector<double>& m_Y_full_initial;
+            const double m_T9;
+            const double m_rho;
+            const Eigen::VectorXd& m_Y_scale;
+
+            EigenFunctor(
+                MultiscalePartitioningEngineView& view,
+                const std::vector<size_t>& qse_solve_indices,
+                const std::vector<double>& Y_full_initial,
+                const double T9,
+                const double rho,
+                const Eigen::VectorXd& Y_scale
+            ) :
+            m_view(&view),
+            m_qse_solve_indices(qse_solve_indices),
+            m_Y_full_initial(Y_full_initial),
+            m_T9(T9),
+            m_rho(rho),
+            m_Y_scale(Y_scale) {}
+
+            int values() const { return m_qse_solve_indices.size(); }
+            int inputs() const { return m_qse_solve_indices.size(); }
+
+            int operator()(const InputType& v_qse, OutputType& f_qse) const;
+            int df(const InputType& v_qse, JacobianType& J_qse) const;
+        };
+
+    private:
+        quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
+        GraphEngine& m_baseEngine; ///< The base engine to which this view delegates calculations.
+        std::vector<QSEGroup> m_qse_groups; ///< The list of identified equilibrium groups.
+        std::vector<fourdst::atomic::Species> m_dynamic_species; ///< The simplified set of species presented to the solver.
+        std::vector<size_t> m_dynamic_species_indices; ///< Indices mapping the dynamic species back to the master engine's list.
+        std::unordered_map<size_t, std::vector<size_t>> m_connectivity_graph;
+
+    private:
+        std::unordered_set<size_t> identifyFastReactions(
+            const std::vector<double>& Y_full,
+            double T9,
+            double rho,
+            double dt_control
+        ) const;
+
+        void buildConnectivityGraph(
+            const std::unordered_set<size_t>& fast_reaction_indices
+        );
+
+        void findConnectedComponents();
+
+        void validateGroupsWithFluxAnalysis(
+            const std::vector<double>& Y,
+            double T9,
+            double rho
+        );
+
+        std::pair<std::vector<fourdst::atomic::Species>, std::vector<size_t>> identifyDynamicSpecies(
+            const std::vector<double>& Y,
+            const std::vector<QSEGroup>& qse_groups,
+            double T9,
+            double rho
+        ) const;
+
+        std::vector<double> solveQSEAbundances(
+            const std::vector<double> &Y_full,
+            double T9,
+            double rho
+        );
+    };
+}

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

@@ -0,0 +1,38 @@
+#pragma once
+
+#include "gridfire/engine/engine_abstract.h"
+#include "gridfire/engine/views/engine_defined.h"
+
+#include "gridfire/network.h"
+
+#include "fourdst/logging/logging.h"
+#include "fourdst/composition/atomicSpecies.h"
+#include "fourdst/composition/composition.h"
+
+#include "quill/Logger.h"
+
+#include <vector>
+#include <string>
+
+namespace gridfire {
+
+    class NetworkPrimingEngineView final : public DefinedEngineView {
+    public:
+        NetworkPrimingEngineView(const std::string& primingSymbol, DynamicEngine& baseEngine);
+        NetworkPrimingEngineView(const fourdst::atomic::Species& primingSpecies, DynamicEngine& baseEngine);
+        const std::vector<std::string>& getPrimingReactionNames() const { return m_peNames; }
+        const fourdst::atomic::Species& getPrimingSpecies() const { return m_primingSpecies; }
+
+
+    private:
+        quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
+        std::vector<std::string> m_peNames; ///< Names of the priming reactions.
+        fourdst::atomic::Species m_primingSpecies; ///< The priming species, if specified.
+    private:
+        std::vector<std::string> constructPrimingReactionSet(
+            const fourdst::atomic::Species& primingSpecies,
+            const DynamicEngine& baseEngine
+        ) const;
+    };
+
+}

src/network/include/gridfire/io/network_file.h

@@ -9,23 +9,7 @@
 #include <vector>
 
 namespace gridfire::io {
-
-    /**
-     * @struct ParsedNetworkData
-     * @brief Holds the data parsed from a network file.
-     *
-     * This struct is used to return the results of parsing a reaction network
-     * file. It contains the list of reaction names that define the network.
-     */
-    struct ParsedNetworkData {
-        /**
-         * @brief A vector of reaction names in their PEN-style format.
-         *
-         * Projectile, Ejectile style names p(p,e+)d is a common format for representing
-         * nuclear reactions as strings.
-         */
-        std::vector<std::string> reactionPENames;
-    };
+    typedef std::vector<std::string> ParsedNetworkData;
 
     /**
      * @class NetworkFileParser

src/network/include/gridfire/reaction/reaction.h

@@ -411,6 +411,8 @@ namespace gridfire::reaction {
          */
         explicit TemplatedReactionSet(std::vector<ReactionT> reactions);
 
+        TemplatedReactionSet();
+
         /**
          * @brief Copy constructor.
          * @param other The ReactionSet to copy.
@@ -577,6 +579,9 @@ namespace gridfire::reaction {
         }
     }
 
+    template<typename ReactionT>
+    TemplatedReactionSet<ReactionT>::TemplatedReactionSet() {}
+
     template <typename ReactionT>
     TemplatedReactionSet<ReactionT>::TemplatedReactionSet(const TemplatedReactionSet<ReactionT> &other) {
         m_reactions.reserve(other.m_reactions.size());

src/network/include/gridfire/screening/screening_weak.h

@@ -145,7 +145,7 @@ namespace gridfire::screening {
         const T T9,
         const T rho
     ) const {
-        LOG_TRACE_L1(
+        LOG_TRACE_L3(
             m_logger,
             "Calculating weak screening factors for {} reactions...",
             reactions.size()

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

@@ -205,6 +205,7 @@ namespace gridfire::solver {
             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;
 
@@ -316,11 +317,13 @@ namespace gridfire::solver {
             };
 
             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.
@@ -337,14 +340,16 @@ namespace gridfire::solver {
                 const std::vector<size_t>& dynamicSpeciesIndices,
                 const std::vector<size_t>& QSESpeciesIndices,
                 const double T9,
-                const double rho
+                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_rho(rho),
+            m_YScale(YScale) {}
 
             int values() const { return m_QSESpeciesIndices.size(); }
             int inputs() const { return m_QSESpeciesIndices.size(); }
@@ -493,9 +498,20 @@ namespace gridfire::solver {
     };
 
     template<typename T>
-    int QSENetworkSolver::EigenFunctor<T>::operator()(const InputType &v_QSE_log, OutputType &f_QSE) const {
+    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 = v_QSE_log.array().exp();
+        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);
@@ -503,6 +519,17 @@ namespace gridfire::solver {
 
 
         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]];
@@ -511,9 +538,9 @@ namespace gridfire::solver {
     }
 
     template <typename T>
-    int QSENetworkSolver::EigenFunctor<T>::df(const InputType& v_QSE_log, JacobianType& J_QSE) const {
+    int QSENetworkSolver::EigenFunctor<T>::df(const InputType& v_QSE, JacobianType& J_QSE) const {
         std::vector<double> y = m_YFull;
-        Eigen::VectorXd Y_QSE = v_QSE_log.array().exp();
+        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);
@@ -528,8 +555,32 @@ namespace gridfire::solver {
             }
         }
 
+        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) {
-            J_QSE(j, j) *= Y_QSE(j); // chain rule for log space transformation
+            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

@@ -26,12 +26,7 @@
 namespace gridfire {
     GraphEngine::GraphEngine(
         const fourdst::composition::Composition &composition
-    ):
-     m_reactions(build_reaclib_nuclear_network(composition, false)),
-     m_partitionFunction(std::make_unique<partition::GroundStatePartitionFunction>()){
-        syncInternalMaps();
-        precomputeNetwork();
-    }
+    ): GraphEngine(composition, partition::GroundStatePartitionFunction()) {}
 
     GraphEngine::GraphEngine(
         const fourdst::composition::Composition &composition,
@@ -40,7 +35,6 @@ namespace gridfire {
     m_partitionFunction(partitionFunction.clone())  // Clone the partition function to ensure ownership
     {
         syncInternalMaps();
-        precomputeNetwork();
     }
 
     GraphEngine::GraphEngine(
@@ -48,7 +42,6 @@ namespace gridfire {
     ) :
     m_reactions(reactions) {
         syncInternalMaps();
-        precomputeNetwork();
     }
 
     StepDerivatives<double> GraphEngine::calculateRHSAndEnergy(
@@ -58,13 +51,17 @@ namespace gridfire {
     ) const {
         if (m_usePrecomputation) {
             std::vector<double> bare_rates;
+            std::vector<double> bare_reverse_rates;
             bare_rates.reserve(m_reactions.size());
+            bare_reverse_rates.reserve(m_reactions.size());
+
             for (const auto& reaction: m_reactions) {
                 bare_rates.push_back(reaction.calculate_rate(T9));
+                bare_reverse_rates.push_back(calculateReverseRate(reaction, T9));
             }
 
             // --- The public facing interface can always use the precomputed version since taping is done internally ---
-            return calculateAllDerivativesUsingPrecomputation(Y, bare_rates, T9, rho);
+            return calculateAllDerivativesUsingPrecomputation(Y, bare_rates, bare_reverse_rates, T9, rho);
         } else {
             return calculateAllDerivatives<double>(Y, T9, rho);
         }
@@ -72,12 +69,17 @@ namespace gridfire {
 
 
     void GraphEngine::syncInternalMaps() {
+        LOG_INFO(m_logger, "Synchronizing internal maps for REACLIB graph network (serif::network::GraphNetwork)...");
         collectNetworkSpecies();
         populateReactionIDMap();
         populateSpeciesToIndexMap();
+        collectAtomicReverseRateAtomicBases();
         generateStoichiometryMatrix();
         reserveJacobianMatrix();
         recordADTape();
+        precomputeNetwork();
+        LOG_INFO(m_logger, "Internal maps synchronized. Network contains {} species and {} reactions.",
+                 m_networkSpecies.size(), m_reactions.size());
     }
 
     // --- Network Graph Construction Methods ---
@@ -234,17 +236,22 @@ namespace gridfire {
 
     double GraphEngine::calculateReverseRate(
         const reaction::Reaction &reaction,
-        const double T9,
-        const double expFactor
+        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());
+            return 0.0; // If reverse reactions are not used, return 0.0
+        }
+        const double expFactor = std::exp(-reaction.qValue() / (m_constants.kB * T9));
         double reverseRate = 0.0;
         const double forwardRate = reaction.calculate_rate(T9);
 
         if (reaction.reactants().size() == 2 && reaction.products().size() == 2) {
             reverseRate = calculateReverseRateTwoBody(reaction, T9, forwardRate, expFactor);
         } else {
-            LOG_WARNING(m_logger, "Reverse rate calculation for reactions with more than two reactants or products is not implemented.");
+            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);
         return reverseRate;
     }
 
@@ -298,7 +305,9 @@ namespace gridfire {
         );
 
         // Accumulate partition functions
-        auto pf_op = [&](double acc, const auto& species) { return acc * m_partitionFunction->evaluate(species.z(), species.a(), T9); };
+        auto pf_op = [&](double acc, const auto& species) {
+            return acc * m_partitionFunction->evaluate(species.z(), species.a(), T9);
+        };
         const double partitionFunctionNumerator = std::accumulate(
             reaction.reactants().begin(),
             reaction.reactants().end(),
@@ -312,9 +321,14 @@ namespace gridfire {
             pf_op
         );
 
-        const double CT = std::pow(massNumerator/massDenominator, 1.5) * (partitionFunctionNumerator/partitionFunctionDenominator);
+        const double CT = std::pow(massNumerator/massDenominator, 1.5) *
+            (partitionFunctionNumerator/partitionFunctionDenominator);
 
-        return forwardRate * symmetryFactor * CT * expFactor;
+        double reverseRate = forwardRate * symmetryFactor * CT * expFactor;
+        if (!std::isfinite(reverseRate)) {
+            return 0.0; // If the reverse rate is not finite, return 0.0
+        }
+        return reverseRate; // Return the calculated reverse rate
 
     }
 
@@ -323,6 +337,10 @@ namespace gridfire {
         const double T9,
         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());
+            return 0.0; // If reverse reactions are not used, return 0.0
+        }
         const double d_log_kFwd = reaction.calculate_forward_rate_log_derivative(T9);
 
         auto log_deriv_pf_op = [&](double acc, const auto& species) {
@@ -355,9 +373,30 @@ namespace gridfire {
 
     }
 
+    bool GraphEngine::isUsingReverseReactions() const {
+        return m_useReverseReactions;
+    }
+
+    void GraphEngine::setUseReverseReactions(const bool useReverse) {
+        m_useReverseReactions = useReverse;
+    }
+
+    int GraphEngine::getSpeciesIndex(const fourdst::atomic::Species &species) const {
+        return m_speciesToIndexMap.at(species); // Returns the index of the species in the stoichiometry matrix
+    }
+
+    std::vector<double> GraphEngine::mapNetInToMolarAbundanceVector(const NetIn &netIn) const {
+        std::vector<double> Y(m_networkSpecies.size(), 0.0); // Initialize with zeros
+        for (const auto& [symbol, entry] : netIn.composition) {
+            Y[getSpeciesIndex(entry.isotope())] = netIn.composition.getMolarAbundance(symbol); // Map species to their molar abundance
+        }
+        return Y; // Return the vector of molar abundances
+    }
+
     StepDerivatives<double> GraphEngine::calculateAllDerivativesUsingPrecomputation(
         const std::vector<double> &Y_in,
         const std::vector<double> &bare_rates,
+        const std::vector<double> &bare_reverse_rates,
         const double T9,
         const double rho
     ) const {
@@ -396,15 +435,31 @@ namespace gridfire {
             const double bare_rate = bare_rates[precomp.reaction_index];
             const double screeningFactor = screeningFactors[precomp.reaction_index];
             const size_t numReactants = m_reactions[precomp.reaction_index].reactants().size();
+            const size_t numProducts = m_reactions[precomp.reaction_index].products().size();
 
-            const double molarReactionFlow =
+            const double forwardMolarReactionFlow =
                     screeningFactor *
                     bare_rate *
                     precomp.symmetry_factor *
                     abundanceProduct *
-                    std::pow(rho, numReactants);
+                    std::pow(rho, numReactants >  1 ? numReactants - 1 : 0);
+
+            double reverseMolarReactionFlow = 0.0;
+            if (precomp.reverse_symmetry_factor != 0.0 and m_useReverseReactions) {
+                const double bare_reverse_rate = bare_reverse_rates[precomp.reaction_index];
+                double reverseAbundanceProduct = 1.0;
+                for (size_t i = 0; i < precomp.unique_product_indices.size(); ++i) {
+                    reverseAbundanceProduct *= std::pow(Y_in[precomp.unique_product_indices[i]], precomp.product_powers[i]);
+                }
+                reverseMolarReactionFlow = screeningFactor *
+                    bare_reverse_rate *
+                    precomp.reverse_symmetry_factor *
+                    reverseAbundanceProduct *
+                    std::pow(rho, numProducts > 1 ? numProducts - 1 : 0);
+            }
+
+            molarReactionFlows.push_back(forwardMolarReactionFlow - reverseMolarReactionFlow);
 
-            molarReactionFlows.push_back(molarReactionFlow);
         }
 
         // --- Assemble molar abundance derivatives ---
@@ -523,7 +578,7 @@ namespace gridfire {
     }
 
     void GraphEngine::generateJacobianMatrix(
-        const std::vector<double> &Y,
+        const std::vector<double> &Y_dynamic,
         const double T9,
         const double rho
     ) {
@@ -534,10 +589,24 @@ 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[i];
+            adInput[i] = Y_dynamic[i];
         }
         adInput[numSpecies]     = T9;  // T9
         adInput[numSpecies + 1] = rho; // rho
+        LOG_DEBUG(
+            m_logger,
+            "AD Input to jacobian {}",
+            [&]() -> std::string {
+                std::stringstream ss;
+                ss << std::scientific << std::setprecision(5);
+                for (size_t i = 0; i < adInput.size(); ++i) {
+                    ss << adInput[i];
+                    if (i < adInput.size() - 1) {
+                        ss << ", ";
+                    }
+                }
+                return ss.str();
+            }());
 
         // 2. Calculate the full jacobian
         const std::vector<double> dotY = m_rhsADFun.Jacobian(adInput);
@@ -552,10 +621,28 @@ namespace gridfire {
                 }
             }
         }
+        LOG_DEBUG(
+            m_logger,
+            "Final Jacobian is:\n{}",
+            [&]() -> std::string {
+                std::stringstream ss;
+                ss << std::scientific << std::setprecision(5);
+                for (size_t i = 0; i < m_jacobianMatrix.size1(); ++i) {
+                    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());
     }
 
     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));
         return m_jacobianMatrix(i, j);
     }
 
@@ -674,8 +761,11 @@ namespace gridfire {
         LOG_TRACE_L1(m_logger, "Successfully exported network graph to {}", filename);
     }
 
-    std::unordered_map<fourdst::atomic::Species, double> GraphEngine::getSpeciesTimescales(const std::vector<double> &Y, const double T9,
-        const double rho) const {
+    std::unordered_map<fourdst::atomic::Species, double> GraphEngine::getSpeciesTimescales(
+        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> speciesTimescales;
         speciesTimescales.reserve(m_networkSpecies.size());
@@ -739,6 +829,19 @@ namespace gridfire {
                  adInput.size());
     }
 
+    void GraphEngine::collectAtomicReverseRateAtomicBases() {
+        m_atomicReverseRates.clear();
+        m_atomicReverseRates.reserve(m_reactions.size());
+
+        for (const auto& reaction: m_reactions) {
+            if (reaction.qValue() != 0.0) {
+                m_atomicReverseRates.push_back(std::make_unique<AtomicReverseRate>(reaction, *this));
+            } else {
+                m_atomicReverseRates.push_back(nullptr);
+            }
+        }
+    }
+
     void GraphEngine::precomputeNetwork() {
         LOG_TRACE_L1(m_logger, "Pre-computing constant components of GraphNetwork state...");
 
@@ -756,7 +859,7 @@ namespace gridfire {
             PrecomputedReaction precomp;
             precomp.reaction_index = i;
 
-            // --- Precompute reactant information ---
+            // --- Precompute forward reaction information ---
             // Count occurrences for each reactant to determine powers and symmetry
             std::unordered_map<size_t, int> reactantCounts;
             for (const auto& reactant: reaction.reactants()) {
@@ -769,10 +872,30 @@ namespace gridfire {
                 precomp.unique_reactant_indices.push_back(index);
                 precomp.reactant_powers.push_back(count);
 
-                symmetryDenominator *= 1.0/std::tgamma(count + 1);
+                symmetryDenominator *= std::tgamma(count + 1);
             }
 
-            precomp.symmetry_factor = symmetryDenominator;
+            precomp.symmetry_factor = 1.0/symmetryDenominator;
+
+            // --- Precompute reverse reaction information ---
+            if (reaction.qValue() != 0.0) {
+                std::unordered_map<size_t, int> productCounts;
+                for (const auto& product : reaction.products()) {
+                    productCounts[speciesIndexMap.at(product)]++;
+                }
+                double reverseSymmetryDenominator = 1.0;
+                for (const auto& [index, count] : productCounts) {
+                    precomp.unique_product_indices.push_back(index);
+                    precomp.product_powers.push_back(count);
+                    reverseSymmetryDenominator *= std::tgamma(count + 1);
+                }
+
+                precomp.reverse_symmetry_factor = 1.0/reverseSymmetryDenominator;
+            } else {
+                precomp.unique_product_indices.clear();
+                precomp.product_powers.clear();
+                precomp.reverse_symmetry_factor = 0.0; // No reverse reaction for Q = 0 reactions
+            }
 
             // --- Precompute stoichiometry information ---
             const auto stoichiometryMap = reaction.stoichiometry();
@@ -788,7 +911,7 @@ namespace gridfire {
         }
     }
 
-    bool AtomicReverseRate::forward(
+    bool GraphEngine::AtomicReverseRate::forward(
         const size_t p,
         const size_t q,
         const CppAD::vector<bool> &vx,
@@ -796,12 +919,53 @@ namespace gridfire {
         const CppAD::vector<double> &tx,
         CppAD::vector<double> &ty
     ) {
+
+        if ( p != 0) { return false; }
         const double T9 = tx[0];
-        const double expFactor = std::exp(-m_reaction.qValue() / (m_kB * T9 * 1e9)); // Convert MeV to erg
-        if (p == 0) {
-            // --- Zeroth order forward sweep ---
-            const auto k_rev = m_engine.calculateReverseRate(m_reaction, T9, expFactor);
+
+        const double reverseRate = m_engine.calculateReverseRate(m_reaction, T9);
+        // std::cout << m_reaction.peName() << " reverseRate: " << reverseRate << " at T9: " << T9 << "\n";
+        ty[0] = reverseRate; // Store the reverse rate in the output vector
+
+        if (vx.size() > 0) {
+            vy[0] = vx[0];
         }
+        return true;
     }
 
+    bool GraphEngine::AtomicReverseRate::reverse(
+        size_t q,
+        const CppAD::vector<double> &tx,
+        const CppAD::vector<double> &ty,
+        CppAD::vector<double> &px,
+        const CppAD::vector<double> &py
+    ) {
+        const double T9 = tx[0];
+        const double reverseRate = ty[0];
+
+        const double derivative = m_engine.calculateReverseRateTwoBodyDerivative(m_reaction, T9, reverseRate);
+        // std::cout << m_reaction.peName() << " reverseRate Derivative: " << derivative << "\n";
+
+        px[0] = py[0] * derivative; // Return the derivative of the reverse rate with respect to T9
+
+        return true;
+    }
+
+    bool GraphEngine::AtomicReverseRate::for_sparse_jac(
+        size_t q,
+        const CppAD::vector<std::set<size_t>> &r,
+        CppAD::vector<std::set<size_t>> &s
+    ) {
+        s[0] = r[0];
+        return true;
+    }
+
+    bool GraphEngine::AtomicReverseRate::rev_sparse_jac(
+        size_t q,
+        const CppAD::vector<std::set<size_t>> &rt,
+        CppAD::vector<std::set<size_t>> &st
+    ) {
+        st[0] = rt[0];
+        return true;
+    }
 }

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

@@ -0,0 +1,84 @@
+#include "gridfire/engine/procedures/priming.h"
+#include "gridfire/engine/views/engine_priming.h"
+#include "gridfire/solver/solver.h"
+
+#include "gridfire/engine/engine_abstract.h"
+#include "gridfire/network.h"
+
+namespace gridfire {
+    fourdst::composition::Composition primeNetwork(const NetIn& netIn, DynamicEngine& engine) {
+        using fourdst::composition::Composition;
+        using fourdst::atomic::Species;
+
+
+        NetIn currentConditions = netIn;
+        std::vector<Species> speciesToPrime;
+        for (const auto &entry: netIn.composition | std::views::values) {
+            if (entry.mass_fraction() == 0.0) {
+                speciesToPrime.push_back(entry.isotope());
+            }
+        }
+
+
+        if (speciesToPrime.empty()) {
+            return netIn.composition; // No priming needed
+        }
+
+        const double T9 = netIn.temperature / 1e9; // Convert temperature to T9
+        const double rho = netIn.density; // Density in g/cm^3
+
+        for (const auto& primingSpecies :speciesToPrime) {
+            NetworkPrimingEngineView primer(primingSpecies, engine);
+
+            const auto Y = primer.mapNetInToMolarAbundanceVector(currentConditions);
+            const double destructionRateConstant = calculateDestructionRateConstant(primer, primingSpecies, Y, T9, rho);
+            const double creationRate = calculateCreationRate(primer, primingSpecies, Y, T9, rho);
+            const double equilibriumMolarAbundance = creationRate / destructionRateConstant;
+            const double equilibriumMassFraction = equilibriumMolarAbundance * primingSpecies.mass();
+
+            currentConditions.composition.setMassFraction(
+                std::string(primingSpecies.name()),
+                equilibriumMassFraction
+            );
+        }
+
+        return currentConditions.composition;
+    }
+
+    double calculateDestructionRateConstant(
+        const DynamicEngine& engine,
+        const fourdst::atomic::Species& species,
+        const std::vector<double>& Y,
+        const double T9,
+        const double rho
+    ) {
+        const int speciesIndex = engine.getSpeciesIndex(species);
+        std::vector<double> Y_scaled(Y.begin(), Y.end());
+        Y_scaled[speciesIndex] = 1.0; // Set the abundance of the species to 1.0 for rate constant calculation
+        double destructionRateConstant = 0.0;
+        for (const auto& reaction: engine.getNetworkReactions()) {
+            if (reaction.contains_reactant(species)) {
+                const int stoichiometry = reaction.stoichiometry(species);
+                destructionRateConstant += std::abs(stoichiometry) * engine.calculateMolarReactionFlow(reaction, Y_scaled, T9, rho);
+            }
+        }
+        return destructionRateConstant;
+    }
+
+    double calculateCreationRate(
+        const DynamicEngine& engine,
+        const fourdst::atomic::Species& species,
+        const std::vector<double>& Y,
+        const double T9,
+        const double rho
+    ) {
+        double creationRate = 0.0;
+        for (const auto& reaction: engine.getNetworkReactions()) {
+            const int stoichiometry = reaction.stoichiometry(species);
+            if (stoichiometry > 0) {
+                creationRate += stoichiometry * engine.calculateMolarReactionFlow(reaction, Y, T9, rho);
+            }
+        }
+        return creationRate;
+    }
+}

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

@@ -136,12 +136,12 @@ namespace gridfire {
     }
 
     void AdaptiveEngineView::generateJacobianMatrix(
-        const std::vector<double> &Y_culled,
+        const std::vector<double> &Y_dynamic,
         const double T9,
         const double rho
     ) {
         validateState();
-        const auto Y_full = mapCulledToFull(Y_culled);
+        const auto Y_full = mapCulledToFull(Y_dynamic);
 
         m_baseEngine.generateJacobianMatrix(Y_full, T9, rho);
     }
@@ -221,6 +221,25 @@ namespace gridfire {
         return m_baseEngine.getScreeningModel();
     }
 
+    std::vector<double> AdaptiveEngineView::mapNetInToMolarAbundanceVector(const NetIn &netIn) const {
+        std::vector<double> Y(m_activeSpecies.size(), 0.0); // Initialize with zeros
+        for (const auto& [symbol, entry] : netIn.composition) {
+            Y[getSpeciesIndex(entry.isotope())] = netIn.composition.getMolarAbundance(symbol); // Map species to their molar abundance
+        }
+        return Y; // Return the vector of molar abundances
+    }
+
+    int AdaptiveEngineView::getSpeciesIndex(const fourdst::atomic::Species &species) const {
+        auto it = std::find(m_activeSpecies.begin(), m_activeSpecies.end(), species);
+        if (it != m_activeSpecies.end()) {
+            return static_cast<int>(std::distance(m_activeSpecies.begin(), it));
+        } else {
+            LOG_ERROR(m_logger, "Species '{}' not found in active species list.", species.name());
+            m_logger->flush_log();
+            throw std::runtime_error("Species not found in active species list: " + std::string(species.name()));
+        }
+    }
+
     std::vector<double> AdaptiveEngineView::mapCulledToFull(const std::vector<double>& culled) const {
         std::vector<double> full(m_baseEngine.getNetworkSpecies().size(), 0.0);
         for (size_t i_culled = 0; i_culled < culled.size(); ++i_culled) {

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

@@ -1,32 +1,88 @@
 #include "gridfire/engine/views/engine_defined.h"
 
+#include <ranges>
+
 #include "quill/LogMacros.h"
 
+#include <string>
+#include <vector>
+#include <unordered_set>
+#include <stdexcept>
+#include <unordered_map>
+#include <utility>
+
 namespace gridfire {
     using fourdst::atomic::Species;
 
-    FileDefinedEngineView::FileDefinedEngineView(
-        DynamicEngine &baseEngine,
-        const std::string &fileName,
-        const io::NetworkFileParser &parser
-    ):
-    m_baseEngine(baseEngine),
-    m_fileName(fileName),
-    m_parser(parser),
-    m_activeSpecies(baseEngine.getNetworkSpecies()),
-    m_activeReactions(baseEngine.getNetworkReactions()) {
-        buildFromFile(fileName);
+    DefinedEngineView::DefinedEngineView(const std::vector<std::string>& peNames, DynamicEngine& baseEngine) :
+    m_baseEngine(baseEngine) {
+        m_activeReactions.clear();
+        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 & FileDefinedEngineView::getBaseEngine() const {
+    const DynamicEngine & DefinedEngineView::getBaseEngine() const {
         return m_baseEngine;
     }
 
-    const std::vector<Species> & FileDefinedEngineView::getNetworkSpecies() const {
+
+
+    const std::vector<Species> & DefinedEngineView::getNetworkSpecies() const {
         return m_activeSpecies;
     }
 
-    StepDerivatives<double> FileDefinedEngineView::calculateRHSAndEnergy(
+    StepDerivatives<double> DefinedEngineView::calculateRHSAndEnergy(
         const std::vector<double> &Y_defined,
         const double T9,
         const double rho
@@ -42,18 +98,18 @@ namespace gridfire {
         return definedResults;
     }
 
-    void FileDefinedEngineView::generateJacobianMatrix(
-        const std::vector<double> &Y_defined,
+    void DefinedEngineView::generateJacobianMatrix(
+        const std::vector<double> &Y_dynamic,
         const double T9,
         const double rho
     ) {
         validateNetworkState();
 
-        const auto Y_full = mapViewToFull(Y_defined);
+        const auto Y_full = mapViewToFull(Y_dynamic);
         m_baseEngine.generateJacobianMatrix(Y_full, T9, rho);
     }
 
-    double FileDefinedEngineView::getJacobianMatrixEntry(
+    double DefinedEngineView::getJacobianMatrixEntry(
         const int i_defined,
         const int j_defined
     ) const {
@@ -65,13 +121,13 @@ namespace gridfire {
         return m_baseEngine.getJacobianMatrixEntry(i_full, j_full);
     }
 
-    void FileDefinedEngineView::generateStoichiometryMatrix() {
+    void DefinedEngineView::generateStoichiometryMatrix() {
         validateNetworkState();
 
         m_baseEngine.generateStoichiometryMatrix();
     }
 
-    int FileDefinedEngineView::getStoichiometryMatrixEntry(
+    int DefinedEngineView::getStoichiometryMatrixEntry(
         const int speciesIndex_defined,
         const int reactionIndex_defined
     ) const {
@@ -82,7 +138,7 @@ namespace gridfire {
         return m_baseEngine.getStoichiometryMatrixEntry(i_full, j_full);
     }
 
-    double FileDefinedEngineView::calculateMolarReactionFlow(
+    double DefinedEngineView::calculateMolarReactionFlow(
         const reaction::Reaction &reaction,
         const std::vector<double> &Y_defined,
         const double T9,
@@ -91,7 +147,7 @@ namespace gridfire {
         validateNetworkState();
 
         if (!m_activeReactions.contains(reaction)) {
-            LOG_ERROR(m_logger, "Reaction '{}' is not part of the active reactions in the file defined engine view.", reaction.id());
+            LOG_ERROR(m_logger, "Reaction '{}' is not part of the active reactions in the DefinedEngineView.", reaction.id());
             m_logger -> flush_log();
             throw std::runtime_error("Reaction not found in active reactions: " + std::string(reaction.id()));
         }
@@ -99,13 +155,13 @@ namespace gridfire {
         return m_baseEngine.calculateMolarReactionFlow(reaction, Y_full, T9, rho);
     }
 
-    const reaction::LogicalReactionSet & FileDefinedEngineView::getNetworkReactions() const {
+    const reaction::LogicalReactionSet & DefinedEngineView::getNetworkReactions() const {
         validateNetworkState();
 
         return m_activeReactions;
     }
 
-    std::unordered_map<Species, double> FileDefinedEngineView::getSpeciesTimescales(
+    std::unordered_map<Species, double> DefinedEngineView::getSpeciesTimescales(
         const std::vector<double> &Y_defined,
         const double T9,
         const double rho
@@ -124,28 +180,45 @@ namespace gridfire {
         return definedTimescales;
     }
 
-    void FileDefinedEngineView::update(const NetIn &netIn) {
-        if (m_isStale) {
-            buildFromFile(m_fileName);
-        }
+    void DefinedEngineView::update(const NetIn &netIn) {
+        return;
     }
 
-    void FileDefinedEngineView::setNetworkFile(const std::string &fileName) {
-        m_isStale = true;
-        m_fileName = fileName;
-        LOG_DEBUG(m_logger, "File '{}' set to '{}'. FileDefinedNetworkView is now stale! You MUST call update() before you use it!", m_fileName, fileName);
-    }
 
-    void FileDefinedEngineView::setScreeningModel(const screening::ScreeningType model) {
+    void DefinedEngineView::setScreeningModel(const screening::ScreeningType model) {
         m_baseEngine.setScreeningModel(model);
     }
 
-    screening::ScreeningType FileDefinedEngineView::getScreeningModel() const {
+    screening::ScreeningType DefinedEngineView::getScreeningModel() const {
         return m_baseEngine.getScreeningModel();
     }
 
-    std::vector<size_t> FileDefinedEngineView::constructSpeciesIndexMap() const {
-        LOG_TRACE_L1(m_logger, "Constructing species index map for file defined engine view...");
+    int DefinedEngineView::getSpeciesIndex(const Species &species) const {
+        validateNetworkState();
+
+        auto it = std::find(m_activeSpecies.begin(), m_activeSpecies.end(), species);
+        if (it != m_activeSpecies.end()) {
+            return static_cast<int>(std::distance(m_activeSpecies.begin(), it));
+        } else {
+            LOG_ERROR(m_logger, "Species '{}' not found in active species list.", species.name());
+            m_logger->flush_log();
+            throw std::runtime_error("Species not found in active species list: " + std::string(species.name()));
+        }
+    }
+
+    std::vector<double> DefinedEngineView::mapNetInToMolarAbundanceVector(const NetIn &netIn) const {
+        std::vector<double> Y(m_activeSpecies.size(), 0.0); // Initialize with zeros
+        for (const auto& [symbol, entry] : netIn.composition) {
+            auto it = std::ranges::find(m_activeSpecies, entry.isotope());
+            if (it != m_activeSpecies.end()) {
+                Y[getSpeciesIndex(entry.isotope())] = netIn.composition.getMolarAbundance(symbol); // Map species to their molar abundance
+            }
+        }
+        return Y; // Return the vector of molar abundances
+    }
+
+    std::vector<size_t> DefinedEngineView::constructSpeciesIndexMap() const {
+        LOG_TRACE_L1(m_logger, "Constructing species index map for DefinedEngineView...");
         std::unordered_map<Species, size_t> fullSpeciesReverseMap;
         const auto& fullSpeciesList = m_baseEngine.getNetworkSpecies();
 
@@ -173,8 +246,8 @@ namespace gridfire {
 
     }
 
-    std::vector<size_t> FileDefinedEngineView::constructReactionIndexMap() const {
-        LOG_TRACE_L1(m_logger, "Constructing reaction index map for file defined engine view...");
+    std::vector<size_t> DefinedEngineView::constructReactionIndexMap() const {
+        LOG_TRACE_L1(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;
@@ -205,66 +278,7 @@ namespace gridfire {
         return reactionIndexMap;
     }
 
-    void FileDefinedEngineView::buildFromFile(const std::string &fileName) {
-        LOG_TRACE_L1(m_logger, "Building file defined engine view from {}...", fileName);
-        auto [reactionPENames] = m_parser.parse(fileName);
-
-        m_activeReactions.clear();
-        m_activeSpecies.clear();
-
-        std::unordered_set<Species> seenSpecies;
-
-        const auto& fullNetworkReactionSet = m_baseEngine.getNetworkReactions();
-        for (const auto& peName : reactionPENames) {
-            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, "File defined engine view 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;
-    }
-
-    std::vector<double> FileDefinedEngineView::mapViewToFull(const std::vector<double>& culled) const {
+    std::vector<double> DefinedEngineView::mapViewToFull(const std::vector<double>& culled) const {
         std::vector<double> full(m_baseEngine.getNetworkSpecies().size(), 0.0);
         for (size_t i_culled = 0; i_culled < culled.size(); ++i_culled) {
             const size_t i_full = m_speciesIndexMap[i_culled];
@@ -273,7 +287,7 @@ namespace gridfire {
         return full;
     }
 
-    std::vector<double> FileDefinedEngineView::mapFullToView(const std::vector<double>& full) const {
+    std::vector<double> DefinedEngineView::mapFullToView(const std::vector<double>& full) const {
         std::vector<double> culled(m_activeSpecies.size(), 0.0);
         for (size_t i_culled = 0; i_culled < m_activeSpecies.size(); ++i_culled) {
             const size_t i_full = m_speciesIndexMap[i_culled];
@@ -282,7 +296,7 @@ namespace gridfire {
         return culled;
     }
 
-    size_t FileDefinedEngineView::mapViewToFullSpeciesIndex(size_t culledSpeciesIndex) const {
+    size_t DefinedEngineView::mapViewToFullSpeciesIndex(size_t culledSpeciesIndex) const {
         if (culledSpeciesIndex < 0 || culledSpeciesIndex >= static_cast<int>(m_speciesIndexMap.size())) {
             LOG_ERROR(m_logger, "Defined index {} is out of bounds for species index map of size {}.", culledSpeciesIndex, m_speciesIndexMap.size());
             m_logger->flush_log();
@@ -291,7 +305,7 @@ namespace gridfire {
         return m_speciesIndexMap[culledSpeciesIndex];
     }
 
-    size_t FileDefinedEngineView::mapViewToFullReactionIndex(size_t culledReactionIndex) const {
+    size_t DefinedEngineView::mapViewToFullReactionIndex(size_t culledReactionIndex) const {
         if (culledReactionIndex < 0 || culledReactionIndex >= static_cast<int>(m_reactionIndexMap.size())) {
             LOG_ERROR(m_logger, "Defined index {} is out of bounds for reaction index map of size {}.", culledReactionIndex, m_reactionIndexMap.size());
             m_logger->flush_log();
@@ -300,11 +314,25 @@ namespace gridfire {
         return m_reactionIndexMap[culledReactionIndex];
     }
 
-    void FileDefinedEngineView::validateNetworkState() const {
+    void DefinedEngineView::validateNetworkState() const {
         if (m_isStale) {
-            LOG_ERROR(m_logger, "File defined engine view is stale. Please call update() with a valid NetIn object.");
+            LOG_ERROR(m_logger, "DefinedEngineView is stale. Please call update() with a valid NetIn object.");
             m_logger->flush_log();
-            throw std::runtime_error("File defined engine view is stale. Please call update() with a valid NetIn object.");
+            throw std::runtime_error("DefinedEngineView is stale. Please call update() with a valid NetIn object.");
         }
     }
+
+
+    ////////////////////////////////////////////
+    /// FileDefinedEngineView Implementation ///
+    /////////////////////////////////////////////
+
+    FileDefinedEngineView::FileDefinedEngineView(
+        DynamicEngine &baseEngine,
+        const std::string &fileName,
+        const io::NetworkFileParser &parser
+    ):
+    DefinedEngineView(parser.parse(fileName), baseEngine),
+    m_fileName(fileName),
+    m_parser(parser) {}
 }

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

@@ -0,0 +1,623 @@
+#include "gridfire/engine/views/engine_multiscale.h"
+
+#include <vector>
+#include <unordered_map>
+#include <unordered_set>
+#include <queue>
+
+namespace gridfire {
+    using fourdst::atomic::Species;
+
+    MultiscalePartitioningEngineView::MultiscalePartitioningEngineView(
+        GraphEngine& baseEngine
+    ) : m_baseEngine(baseEngine) {}
+
+    const std::vector<Species> & MultiscalePartitioningEngineView::getNetworkSpecies() const {
+        return m_baseEngine.getNetworkSpecies();
+    }
+
+    StepDerivatives<double> MultiscalePartitioningEngineView::calculateRHSAndEnergy(
+        const std::vector<double> &Y,
+        const double T9,
+        const double rho
+    ) const {
+        return m_baseEngine.calculateRHSAndEnergy(Y, T9, rho);
+    }
+
+    void MultiscalePartitioningEngineView::generateJacobianMatrix(
+        const std::vector<double> &Y_dynamic,
+        const double T9,
+        const double rho
+    ) {
+        std::vector<double> Y_full(m_baseEngine.getNetworkSpecies().size(), 0.0);
+        for (size_t i = 0; i < m_dynamic_species_indices.size(); ++i) {
+            Y_full[m_dynamic_species_indices[i]] = Y_dynamic[i];
+        }
+        // solveQSEAbundances(Y_full, T9, rho);
+        m_baseEngine.generateJacobianMatrix(Y_dynamic, T9, rho);
+    }
+
+    double MultiscalePartitioningEngineView::getJacobianMatrixEntry(
+        const int i,
+        const int j
+    ) const {
+        return m_baseEngine.getJacobianMatrixEntry(i, j);
+    }
+
+    void MultiscalePartitioningEngineView::generateStoichiometryMatrix() {
+        m_baseEngine.generateStoichiometryMatrix();
+    }
+
+    int MultiscalePartitioningEngineView::getStoichiometryMatrixEntry(
+        const int speciesIndex,
+        const int reactionIndex
+    ) const {
+        return m_baseEngine.getStoichiometryMatrixEntry(speciesIndex, reactionIndex);
+    }
+
+    double MultiscalePartitioningEngineView::calculateMolarReactionFlow(
+        const reaction::Reaction &reaction,
+        const std::vector<double> &Y,
+        double T9,
+        double rho
+    ) const {
+        return m_baseEngine.calculateMolarReactionFlow(reaction, Y, T9, rho);
+    }
+
+    const reaction::LogicalReactionSet & MultiscalePartitioningEngineView::getNetworkReactions() const {
+        return m_baseEngine.getNetworkReactions();
+    }
+
+    std::unordered_map<Species, double> MultiscalePartitioningEngineView::getSpeciesTimescales(
+        const std::vector<double> &Y,
+        const double T9,
+        const double rho
+    ) const {
+        return m_baseEngine.getSpeciesTimescales(Y, T9, rho);
+    }
+
+    void MultiscalePartitioningEngineView::update(const NetIn &netIn) {
+        return; // call partition eventually
+    }
+
+    void MultiscalePartitioningEngineView::setScreeningModel(
+        const screening::ScreeningType model
+    ) {
+        m_baseEngine.setScreeningModel(model);
+    }
+
+    screening::ScreeningType MultiscalePartitioningEngineView::getScreeningModel() const {
+        return m_baseEngine.getScreeningModel();
+    }
+
+    const DynamicEngine & MultiscalePartitioningEngineView::getBaseEngine() const {
+        return m_baseEngine;
+    }
+
+    void MultiscalePartitioningEngineView::partitionNetwork(
+        const std::vector<double> &Y,
+        const double T9,
+        const double rho,
+        const double dt_control
+    ) {
+        // --- Step 0. Clear previous state ---
+        m_qse_groups.clear();
+        m_dynamic_species.clear();
+        m_dynamic_species_indices.clear();
+        m_connectivity_graph.clear();
+
+        // --- Step 1. Identify fast reactions ---
+        const std::unordered_set<size_t> fast_reaction_indices = identifyFastReactions(Y, T9, rho, dt_control);
+
+        // --- Step 2. Build connectivity graph based on fast reactions ---
+        buildConnectivityGraph(fast_reaction_indices);
+
+        // --- Step 3. Find connected components (candidate QSE clusters) ---
+        findConnectedComponents();
+
+        // --- Step 4. Validate clusters via flux analysis ---
+        validateGroupsWithFluxAnalysis(Y, T9, rho);
+
+        // --- Step 5. Identify dynamic species ---
+        const auto [fst, snd] = identifyDynamicSpecies(Y, m_qse_groups, T9, rho);
+        m_dynamic_species = fst;
+        m_dynamic_species_indices = snd;
+
+    }
+
+    void MultiscalePartitioningEngineView::partitionNetwork(
+        const NetIn &netIn,
+        const double dt_control
+    ) {
+        std::vector<double> Y(m_baseEngine.getNetworkSpecies().size(), 0.0);
+        for (const auto& [symbol, entry]: netIn.composition ) {
+            if (!m_baseEngine.involvesSpecies(entry.isotope())) {
+                LOG_ERROR(m_logger, "Species {} is not part of the network. Exiting...", entry.isotope().name());
+                throw std::runtime_error("Species " + std::string(entry.isotope().name()) + " is not part of 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
+
+        partitionNetwork(Y, T9, rho, dt_control);
+    }
+
+    void MultiscalePartitioningEngineView::exportToDot(const std::string &filename) const {
+        std::ofstream dotFile(filename);
+        if (!dotFile.is_open()) {
+            LOG_ERROR(m_logger, "Failed to open file for writing: {}", filename);
+            throw std::runtime_error("Failed to open file for writing: " + filename);
+        }
+
+        dotFile << "digraph PartitionedNetwork {\n";
+        dotFile << "    graph [rankdir=LR, splines=true, overlap=false, bgcolor=\"#fdfdfd\", label=\"Multiscale Partitioned Network View\"];\n";
+        dotFile << "    node [shape=circle, style=filled, fillcolor=\"#eafaf1\", fontname=\"Helvetica\"];\n";
+        dotFile << "    edge [fontname=\"Helvetica\", fontsize=10, color=\"#5d6d7e\"];\n\n";
+
+        // 1. Define all species from the base network as nodes
+        const auto& all_species = m_baseEngine.getNetworkSpecies();
+        dotFile << "    // --- Species Nodes ---\n";
+        for (const auto& species : all_species) {
+            dotFile << "    \"" << species.name() << "\" [label=\"" << species.name() << "\"];\n";
+        }
+        dotFile << "\n";
+
+        // 2. Define all reactions and their connections from the base network
+        dotFile << "    // --- Reaction Edges ---\n";
+        for (const auto& reaction : m_baseEngine.getNetworkReactions()) {
+            std::string reactionNodeId = "reaction_" + std::string(reaction.id());
+            dotFile << "    \"" << reactionNodeId << "\" [shape=point, fillcolor=black, width=0.08, height=0.08];\n";
+            for (const auto& reactant : reaction.reactants()) {
+                dotFile << "    \"" << reactant.name() << "\" -> \"" << reactionNodeId << "\";\n";
+            }
+            for (const auto& product : reaction.products()) {
+                dotFile << "    \"" << reactionNodeId << "\" -> \"" << product.name() << "\";\n";
+            }
+            dotFile << "\n";
+        }
+
+        // 3. Draw clusters around the identified QSE groups
+        dotFile << "    // --- QSE Group Clusters ---\n";
+        int group_counter = 0;
+        for (const auto& group : m_qse_groups) {
+            dotFile << "    subgraph cluster_" << group_counter << " {\n";
+            dotFile << "        style = \"filled,rounded\";\n";
+            if (group.is_in_equilibrium) {
+                dotFile << "        label = \"Valid QSE Group" << group_counter + 1 << "\";\n";
+                dotFile << "        color = \"#fdebd0\";\n";
+            } else {
+                dotFile << "        label = \"Non valid QSE Group" << group_counter + 1 << "\";\n";
+                dotFile << "        color = \"#d5f5e3\";\n"; // Different color for non-equilibrium groups
+            }
+
+            for (const size_t species_idx : group.species_indices) {
+                dotFile << "        \"" << all_species[species_idx].name() << "\";\n";
+            }
+
+            dotFile << "    }\n";
+            group_counter++;
+        }
+
+        dotFile << "}\n";
+        dotFile.close();
+    }
+
+    std::vector<double> MultiscalePartitioningEngineView::mapNetInToMolarAbundanceVector(const NetIn &netIn) const {
+        std::vector<double> Y(m_dynamic_species.size(), 0.0); // Initialize with zeros
+        for (const auto& [symbol, entry] : netIn.composition) {
+            Y[getSpeciesIndex(entry.isotope())] = netIn.composition.getMolarAbundance(symbol); // Map species to their molar abundance
+        }
+        return Y; // Return the vector of molar abundances
+    }
+
+    std::vector<Species> MultiscalePartitioningEngineView::getFastSpecies() const {
+        const auto& all_species = m_baseEngine.getNetworkSpecies();
+        std::vector<Species> fast_species;
+        fast_species.reserve(all_species.size() - m_dynamic_species.size());
+        for (const auto& species : all_species) {
+            auto it = std::ranges::find(m_dynamic_species, species);
+            if (it == m_dynamic_species.end()) {
+                fast_species.push_back(species);
+            }
+        }
+        return fast_species;
+    }
+
+    const std::vector<Species> & MultiscalePartitioningEngineView::getDynamicSpecies() const {
+        return m_dynamic_species;
+    }
+
+    void MultiscalePartitioningEngineView::equilibrateNetwork(
+        const std::vector<double> &Y,
+        const double T9,
+        const double rho,
+        const double dt_control
+    ) {
+        partitionNetwork(Y, T9, rho, dt_control);
+        std::vector<double> Y_equilibrated = solveQSEAbundances(Y, T9, rho);
+    }
+
+    void MultiscalePartitioningEngineView::equilibrateNetwork(
+        const NetIn& netIn,
+        const double dt_control
+    ) {
+        std::vector<double> Y(m_baseEngine.getNetworkSpecies().size(), 0.0);
+        for (const auto& [symbol, entry]: netIn.composition ) {
+            if (!m_baseEngine.involvesSpecies(entry.isotope())) {
+                LOG_ERROR(m_logger, "Species {} is not part of the network. Exiting...", entry.isotope().name());
+                throw std::runtime_error("Species " + std::string(entry.isotope().name()) + " is not part of 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
+
+        partitionNetwork(Y, T9, rho, dt_control);
+        std::vector<double> Y_equilibrated = solveQSEAbundances(Y, T9, rho);
+
+
+    }
+
+    int MultiscalePartitioningEngineView::getSpeciesIndex(const fourdst::atomic::Species &species) const {
+        return m_baseEngine.getSpeciesIndex(species);
+    }
+
+    std::unordered_set<size_t> MultiscalePartitioningEngineView::identifyFastReactions(
+        const std::vector<double> &Y_full,
+        const double T9,
+        const double rho,
+        const double dt_control
+    ) const {
+        std::unordered_set<size_t> fast_reaction_indices;
+        std::unordered_set<size_t> fast_species_indices;
+
+        const double timescale_threshold = 0.01 * dt_control;
+        constexpr int FAST_BOOTSTRAPPING_ITERATIONS = 2;
+
+        const auto& all_reactions = m_baseEngine.getNetworkReactions();
+
+        // Loop a few times to refine the fast species and reactions.
+        for (int iteration = 0; iteration < FAST_BOOTSTRAPPING_ITERATIONS; ++iteration) {
+            // --- Step A: Identify fast species based on their individual reaction timescales ---
+            for (size_t i = 0; i < all_reactions.size(); ++i) {
+                const auto& reaction = all_reactions[i];
+
+                // Calculate the forward molar flow, which represents the rate of destruction.
+                // Note: We use the base engine's direct calculation method.
+                const double forward_flow = m_baseEngine.calculateMolarReactionFlow(reaction, Y_full, T9, rho);
+
+                if (forward_flow == 0.0) continue;
+
+                // Determine the timescale for each reactant against this specific reaction.
+                for (const auto& reactant : reaction.reactants()) {
+                    // ~ Timescale = Abundance / Destruction Rate
+
+                    // TODO: Should this call base engine getSpeciesIndex or this view's getSpeciesIndex?
+                    const size_t reactant_idx = m_baseEngine.getSpeciesIndex(reactant);
+                    const double abundance = Y_full[reactant_idx];
+
+                    if (abundance > 0) {
+                        const double timescale = abundance / forward_flow;
+                        if (timescale < timescale_threshold) {
+                            fast_species_indices.insert(reactant_idx);
+                        }
+                    }
+                }
+            }
+
+            // --- Step B: Identify fast reactions based on reactants ---
+            for (size_t i = 0; i < all_reactions.size(); ++i) {
+                const auto& reaction = all_reactions[i];
+                bool all_reactants_are_fast = true;
+                if (reaction.reactants().empty()) {
+                    all_reactants_are_fast = false;
+                }
+                for (const auto& reactant : reaction.reactants()) {
+                    if (!fast_species_indices.contains(m_baseEngine.getSpeciesIndex(reactant))) {
+                        all_reactants_are_fast = false;
+                        break; // No need to check further if one reactant is not fast
+                    }
+                }
+                if (all_reactants_are_fast) {
+                    fast_reaction_indices.insert(i);
+                }
+            }
+
+            // --- Step C: Propagate "fast" status to products of the fast reactions (this handles things like n-1) ---
+            for (size_t reaction_idx : fast_reaction_indices) {
+                for (const auto& product : all_reactions[reaction_idx].products()) {
+                    fast_species_indices.insert(m_baseEngine.getSpeciesIndex(product));
+                }
+            }
+        }
+        return fast_reaction_indices;
+    }
+
+    void MultiscalePartitioningEngineView::buildConnectivityGraph(
+        const std::unordered_set<size_t> &fast_reaction_indices
+    ) {
+        const auto& all_reactions = m_baseEngine.getNetworkReactions();
+        for (const size_t reaction_idx : fast_reaction_indices) {
+            const auto& reaction = all_reactions[reaction_idx];
+            const auto& reactants = reaction.reactants();
+            const auto& products = reaction.products();
+
+            // For each fast reaction, create edges between all reactants and all products.
+            // This represents that nucleons can flow quickly between these species.
+            for (const auto& reactant : reactants) {
+                const size_t reactant_idx = m_baseEngine.getSpeciesIndex(reactant);
+                for (const auto& product : products) {
+                    const size_t product_idx = m_baseEngine.getSpeciesIndex(product);
+
+                    // Add a two-way edge to the adjacency list.
+                    m_connectivity_graph[reactant_idx].push_back(product_idx);
+                    m_connectivity_graph[product_idx].push_back(reactant_idx);
+                }
+            }
+        }
+
+
+    }
+
+    void MultiscalePartitioningEngineView::findConnectedComponents() {
+        const size_t num_total_species = m_baseEngine.getNetworkSpecies().size();
+        std::vector<bool> visited(num_total_species, false);
+
+        for (size_t i = 0; i < num_total_species; ++i) {
+            if (!visited[i]) {
+                // Start a new BFS traversal from this unvisited node.
+                // This will find one complete connected component.
+                QSEGroup new_group;
+                std::queue<size_t> q;
+
+                q.push(i);
+                visited[i] = true;
+
+                while (!q.empty()) {
+                    const size_t current_species_idx = q.front();
+                    q.pop();
+                    new_group.species_indices.push_back(current_species_idx);
+
+                    // Check if the node exists in the graph (it might be an isolated species)
+                    if (m_connectivity_graph.contains(current_species_idx)) {
+                        for (const size_t neighbor_idx : m_connectivity_graph.at(current_species_idx)) {
+                            if (!visited[neighbor_idx]) {
+                                visited[neighbor_idx] = true;
+                                q.push(neighbor_idx);
+                            }
+                        }
+                    }
+                }
+
+                // A "group" must contain more than one species to be a cluster.
+                // Isolated species are treated as dynamic by default.
+                if (new_group.species_indices.size() > 1) {
+                    m_qse_groups.push_back(new_group);
+                }
+            }
+        }
+
+    }
+
+    void MultiscalePartitioningEngineView::validateGroupsWithFluxAnalysis(
+        const std::vector<double> &Y,
+        const double T9,
+        const double rho
+    ) {
+        constexpr double FLUX_RATIO_THRESHOLD = 100;
+        for (auto& group : m_qse_groups) {
+            double internal_flux = 0.0;
+            double external_flux = 0.0;
+
+            const std::unordered_set<size_t> group_members(
+                group.species_indices.begin(),
+                group.species_indices.end()
+            );
+
+            for (const auto& reaction: m_baseEngine.getNetworkReactions()) {
+                const double flow = std::abs(m_baseEngine.calculateMolarReactionFlow(reaction, Y, T9, rho));
+                if (flow == 0.0) {
+                    continue; // Skip reactions with zero flow
+                }
+                bool has_internal_reactant = false;
+                bool has_external_reactant = false;
+
+                for (const auto& reactant : reaction.reactants()) {
+                    if (group_members.contains(m_baseEngine.getSpeciesIndex(reactant))) {
+                        has_internal_reactant = true;
+                    } else {
+                        has_external_reactant = true;
+                    }
+                }
+
+                bool has_internal_product = false;
+                bool has_external_product = false;
+
+                for (const auto& product : reaction.products()) {
+                    if (group_members.contains(m_baseEngine.getSpeciesIndex(product))) {
+                        has_internal_product = true;
+                    } else {
+                        has_external_product = true;
+                    }
+                }
+
+                // Classify the reaction based on its participants
+                if ((has_internal_reactant && has_internal_product) && !(has_external_reactant || has_external_product)) {
+                    internal_flux += flow; // Internal flux within the group
+                } else if ((has_internal_reactant || has_internal_product) && (has_external_reactant || has_external_product)) {
+                    external_flux += flow; // External flux to/from the group
+                }
+                // otherwise the reaction is fully decoupled from the QSE group and can be ignored.
+            }
+            if (internal_flux > FLUX_RATIO_THRESHOLD * external_flux) {
+                group.is_in_equilibrium = true; // This group is in equilibrium if internal flux is significantly larger than external flux.
+            } else {
+                group.is_in_equilibrium = false;
+            }
+        }
+    }
+
+    std::pair<std::vector<Species>, std::vector<size_t>> MultiscalePartitioningEngineView::identifyDynamicSpecies(
+        const std::vector<double>& Y,
+        const std::vector<QSEGroup>& qse_groups,
+        const double T9,
+        const double rho
+    ) const {
+        // This set will contain the indices of species that are truly in qse (as opposed to those which might be in a qse group but are not in equilibrium).
+        std::unordered_set<size_t> algebraic_qse_species_indices;
+        auto species_timescales = m_baseEngine.getSpeciesTimescales(Y, T9, rho);
+
+        for (auto& group : qse_groups) {
+            if (!group.is_in_equilibrium || group.species_indices.empty()) continue;
+
+            for (size_t species_idx : group.species_indices) {
+                const double timescale = species_timescales.at(m_baseEngine.getNetworkSpecies()[species_idx]);
+
+                if (timescale < 1) { // If the timescale is less than 1 second, we consider it algebraic.
+                    algebraic_qse_species_indices.insert(species_idx);
+                }
+            }
+        }
+
+        std::vector<fourdst::atomic::Species> dynamic_species;
+        std::unordered_set<size_t> qse_species_indices_set;
+        std::vector<size_t> dynamic_species_indices;
+        const auto& all_base_species = m_baseEngine.getNetworkSpecies();
+
+        for (size_t i = 0; i < all_base_species.size(); ++i) {
+            if (!algebraic_qse_species_indices.contains(i)) {
+                dynamic_species.push_back(all_base_species[i]);
+                dynamic_species_indices.push_back(i);
+            }
+        }
+        return {dynamic_species, dynamic_species_indices};
+    }
+
+    std::vector<double> MultiscalePartitioningEngineView::solveQSEAbundances(
+        const std::vector<double> &Y_full,
+        const double T9,
+        const double rho
+    ) {
+        auto Y_out = Y_full;
+        auto species_timescales = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho);
+
+        for (const auto& qse_group : m_qse_groups) {
+            if (!qse_group.is_in_equilibrium || qse_group.species_indices.empty()) continue;
+            std::vector<size_t> qse_solve_indices;
+            std::vector<size_t> seed_indices;
+
+            for (size_t species_idx : qse_group.species_indices) {
+                const double timescale = species_timescales.at(m_baseEngine.getNetworkSpecies()[species_idx]);
+                if (timescale < 1.0) { // If the timescale is less than 1 second, we consider it algebraic.
+                    qse_solve_indices.push_back(species_idx);
+                } else {
+                    seed_indices.push_back(species_idx);
+                }
+            }
+
+            if (qse_solve_indices.empty()) continue;
+
+            Eigen::VectorXd Y_scale(qse_solve_indices.size());
+            Eigen::VectorXd v_initial(qse_solve_indices.size());
+            for (size_t i = 0; i < qse_solve_indices.size(); ++i) {
+                constexpr double abundance_floor = 1.0e-15;
+                const double initial_abundance = Y_full[qse_solve_indices[i]];
+                Y_scale(i) = std::max(initial_abundance, abundance_floor);
+                v_initial(i) = std::asinh(initial_abundance / Y_scale(i)); // Scale the initial abundances using asinh
+            }
+
+            EigenFunctor functor(*this, qse_solve_indices, Y_full, T9, rho, Y_scale);
+            Eigen::LevenbergMarquardt lm(functor);
+            lm.parameters.ftol = 1.0e-10;
+            lm.parameters.xtol = 1.0e-10;
+
+            Eigen::LevenbergMarquardtSpace::Status status = lm.minimize(v_initial);
+
+            if (status <= 0 || status >= 4) {
+                std::stringstream msg;
+                msg << "QSE solver failed with status: " << status << " for group with seed ";
+                if (seed_indices.size() == 1) {
+                    msg << "nucleus " << m_baseEngine.getNetworkSpecies()[seed_indices[0]].name();
+                } else {
+                    msg << "nuclei ";
+                    // TODO: Refactor nice list printing into utility function somewhere
+                    int counter = 0;
+                    for (const auto& seed_idx : seed_indices) {
+                        msg << m_baseEngine.getNetworkSpecies()[seed_idx].name();
+                        if (counter < seed_indices.size() - 2) {
+                            msg << ", ";
+                        } else if (counter == seed_indices.size() - 2) {
+                            if (seed_indices.size() < 2) {
+                                msg << " and ";
+                            } else {
+                                msg << ", and ";
+                            }
+                        }
+                        ++counter;
+                    }
+                }
+                throw std::runtime_error(msg.str());
+            }
+            Eigen::VectorXd Y_final_qse = Y_scale.array() * v_initial.array().sinh(); // Convert back to physical abundances using asinh scaling
+            for (size_t i = 0; i < qse_solve_indices.size(); ++i) {
+                Y_out[qse_solve_indices[i]] = Y_final_qse(i);
+                std::cout << "Updating species " << m_baseEngine.getNetworkSpecies()[qse_solve_indices[i]].name()
+                          << " to QSE abundance: " << Y_final_qse(i) << std::endl;
+            }
+        }
+        return Y_out;
+    }
+
+    int MultiscalePartitioningEngineView::EigenFunctor::operator()(const InputType &v_qse, OutputType &f_qse) const {
+        std::vector<double> y_trial = m_Y_full_initial;
+        Eigen::VectorXd y_qse = m_Y_scale.array() * v_qse.array().sinh(); // Convert to physical abundances using asinh scaling
+
+        for (size_t i = 0; i < m_qse_solve_indices.size(); ++i) {
+            y_trial[m_qse_solve_indices[i]] = y_qse(i);
+        }
+
+        const auto [dydt, energy] = m_view->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]];
+        }
+        return 0; // Success
+    }
+
+    int MultiscalePartitioningEngineView::EigenFunctor::df(const InputType &v_qse, JacobianType &J_qse) const {
+        std::vector<double> y_trial = m_Y_full_initial;
+        Eigen::VectorXd y_qse = m_Y_scale.array() * v_qse.array().sinh(); // Convert to physical abundances using asinh scaling
+
+        for (size_t i = 0; i < m_qse_solve_indices.size(); ++i) {
+            y_trial[m_qse_solve_indices[i]] = y_qse(i);
+        }
+
+        // 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);
+
+        // 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(
+                    m_qse_solve_indices[i],
+                    m_qse_solve_indices[j]
+                );
+            }
+        }
+
+        // Chain rule for asinh scaling:
+        for (long j = 0; j < J_qse.cols(); ++j) {
+            const double dY_dv = m_Y_scale(j) * std::cosh(v_qse(j));
+            J_qse.col(j) *= dY_dv; // Scale the column by the derivative of the asinh scaling
+        }
+        return 0; // Success
+    }
+
+
+}

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

@@ -0,0 +1,72 @@
+#include "gridfire/engine/views/engine_priming.h"
+#include "gridfire/solver/solver.h"
+
+#include "fourdst/composition/species.h"
+#include "fourdst/logging/logging.h"
+
+#include "quill/LogMacros.h"
+#include "quill/Logger.h"
+
+#include <vector>
+#include <string>
+#include <unordered_set>
+#include <stdexcept>
+#include <unordered_map>
+#include <utility>
+#include <ranges>
+#include <cmath>
+
+
+namespace gridfire {
+    using fourdst::atomic::species;
+
+    NetworkPrimingEngineView::NetworkPrimingEngineView(
+        const std::string &primingSymbol,
+        DynamicEngine &baseEngine
+    ) :
+    DefinedEngineView(
+        constructPrimingReactionSet(
+            species.at(primingSymbol),
+            baseEngine
+        ),
+        baseEngine
+    ),
+    m_primingSpecies(species.at(primingSymbol)) {}
+
+    NetworkPrimingEngineView::NetworkPrimingEngineView(
+        const fourdst::atomic::Species &primingSpecies,
+        DynamicEngine &baseEngine
+    ) :
+    DefinedEngineView(
+        constructPrimingReactionSet(
+            primingSpecies,
+            baseEngine
+        ),
+        baseEngine
+    ),
+    m_primingSpecies(primingSpecies) {}
+
+
+    std::vector<std::string> NetworkPrimingEngineView::constructPrimingReactionSet(
+        const fourdst::atomic::Species &primingSpecies,
+        const DynamicEngine &baseEngine
+    ) const {
+        std::unordered_set<std::string> primeReactions;
+        for (const auto &reaction : baseEngine.getNetworkReactions()) {
+            if (reaction.contains(primingSpecies)) {
+                primeReactions.insert(std::string(reaction.peName()));
+            }
+        }
+        if (primeReactions.empty()) {
+            LOG_ERROR(m_logger, "No priming reactions found for species '{}'.", primingSpecies.name());
+            m_logger->flush_log();
+            throw std::runtime_error("No priming reactions found for species '" + std::string(primingSpecies.name()) + "'.");
+        }
+        std::vector<std::string> primingReactionSet(primeReactions.begin(), primeReactions.end());
+        // LOG_INFO(m_logger, "Constructed priming reaction set with {} reactions for species '{}'.", primingReactionSet.size(), primingSpecies.name());
+        return primingReactionSet;
+    }
+
+
+
+}

src/network/lib/io/network_file.cpp

@@ -68,9 +68,9 @@ namespace gridfire::io {
             if (line.empty()) {
                 continue; // Skip empty lines
             }
-            parsed.reactionPENames.push_back(line);
+            parsed.push_back(line);
         }
-        LOG_TRACE_L1(m_logger, "Parsed {} reactions from file: {}", parsed.reactionPENames.size(), filename);
+        LOG_TRACE_L1(m_logger, "Parsed {} reactions from file: {}", parsed.size(), filename);
         return parsed;
     }
 

src/network/lib/partition/composite/partition_composite.cpp

@@ -25,13 +25,13 @@ namespace gridfire::partition {
     }
 
     double CompositePartitionFunction::evaluate(int z, int a, double T9) const {
-        LOG_TRACE_L1(m_logger, "Evaluating partition function for Z={} A={} T9={}", z, a, T9);
+        LOG_TRACE_L3(m_logger, "Evaluating partition function for Z={} A={} T9={}", z, a, T9);
         for (const auto& partitionFunction : m_partitionFunctions) {
             if (partitionFunction->supports(z, a)) {
-                LOG_TRACE_L2(m_logger, "Partition function of type {} supports Z={} A={}", partitionFunction->type(), z, a);
+                LOG_TRACE_L3(m_logger, "Partition function of type {} supports Z={} A={}", partitionFunction->type(), z, a);
                 return partitionFunction->evaluate(z, a, T9);
             } else {
-                LOG_TRACE_L2(m_logger, "Partition function of type {} does not support Z={} A={}", partitionFunction->type(), z, a);
+                LOG_TRACE_L3(m_logger, "Partition function of type {} does not support Z={} A={}", partitionFunction->type(), z, a);
             }
         }
         LOG_ERROR(
@@ -48,7 +48,7 @@ namespace gridfire::partition {
     double CompositePartitionFunction::evaluateDerivative(int z, int a, double T9) const {
         for (const auto& partitionFunction : m_partitionFunctions) {
             if (partitionFunction->supports(z, a)) {
-                LOG_TRACE_L2(m_logger, "Evaluating derivative of partition function for Z={} A={} T9={}", z, a, T9);
+                LOG_TRACE_L3(m_logger, "Evaluating derivative of partition function for Z={} A={} T9={}", z, a, T9);
                 return partitionFunction->evaluateDerivative(z, a, T9);
             }
         }

src/network/lib/solver/solver.cpp

@@ -2,7 +2,9 @@
 #include "gridfire/engine/engine_graph.h"
 #include "gridfire/network.h"
 
+
 #include "gridfire/utils/logging.h"
+#include "gridfire/utils/qse_rules.h"
 
 #include "fourdst/composition/atomicSpecies.h"
 #include "fourdst/composition/composition.h"
@@ -22,6 +24,7 @@
 #include "quill/LogMacros.h"
 
 namespace gridfire::solver {
+    double s_prevTimestep = 0.0;
 
     NetOut QSENetworkSolver::evaluate(const NetIn &netIn) {
         // --- Use the policy to decide whether to update the view ---
@@ -149,6 +152,15 @@ namespace gridfire::solver {
         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];
@@ -164,9 +176,17 @@ namespace gridfire::solver {
 
             const double final_timescale = std::min(network_timescale, decay_timescale);
 
-            if (std::isinf(final_timescale) || abundance < abundanceCutoff || final_timescale <= timescaleCutoff) {
+            const bool isQSE = is_species_in_qse(
+                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);
             }
         }
@@ -209,11 +229,37 @@ namespace gridfire::solver {
     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,
@@ -221,13 +267,31 @@ namespace gridfire::solver {
             indices.dynamicSpeciesIndices,
             indices.QSESpeciesIndices,
             T9,
-            rho
+            rho,
+            YScale
         );
 
-        Eigen::VectorXd v_qse_log_initial(indices.QSESpeciesIndices.size());
+        Eigen::VectorXd v_qse_initial(indices.QSESpeciesIndices.size());
         for (size_t i = 0; i < indices.QSESpeciesIndices.size(); ++i) {
-            v_qse_log_initial(i) = std::log(std::max(Y[indices.QSESpeciesIndices[i]], 1e-99));
+            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"},
@@ -245,7 +309,9 @@ namespace gridfire::solver {
         };
 
         Eigen::LevenbergMarquardt lm(functor);
-        const Eigen::LevenbergMarquardtSpace::Status info = lm.minimize(v_qse_log_initial);
+        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: {} ({})",
@@ -257,7 +323,7 @@ namespace gridfire::solver {
         }
         LOG_DEBUG(m_logger, "QSE species minimization completed successfully with status: {} ({})",
                   static_cast<int>(info), statusMessages.at(info));
-        return v_qse_log_initial.array().exp();
+        return YScale.array() * v_qse_initial.array().sinh(); // Convert back to molar abundances
 
     }
 
@@ -287,6 +353,24 @@ namespace gridfire::solver {
             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;
@@ -300,6 +384,42 @@ namespace gridfire::solver {
         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;
@@ -352,6 +472,15 @@ namespace gridfire::solver {
         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) {
@@ -445,15 +574,6 @@ namespace gridfire::solver {
         double t
     ) const {
         const std::vector<double> y(Y.begin(), m_numSpecies + Y.begin());
-
-        // std::string timescales = utils::formatNuclearTimescaleLogString(
-        //     m_engine,
-        //     y,
-        //     m_T9,
-        //     m_rho
-        // );
-        // LOG_TRACE_L2(m_logger, "{}", timescales);
-
         auto [dydt, eps] = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho);
         dYdt.resize(m_numSpecies + 1);
         std::ranges::copy(dydt, dYdt.begin());

src/network/meson.build

@@ -5,6 +5,9 @@ network_sources = files(
     'lib/engine/engine_graph.cpp',
     'lib/engine/views/engine_adaptive.cpp',
     'lib/engine/views/engine_defined.cpp',
+    'lib/engine/views/engine_multiscale.cpp',
+    'lib/engine/views/engine_priming.cpp',
+    'lib/engine/procedures/priming.cpp',
     'lib/reaction/reaction.cpp',
     'lib/reaction/reaclib.cpp',
     'lib/io/network_file.cpp',
@@ -53,6 +56,9 @@ network_headers = files(
     'include/gridfire/engine/engine_graph.h',
     'include/gridfire/engine/views/engine_adaptive.h',
     'include/gridfire/engine/views/engine_defined.h',
+    'include/gridfire/engine/views/engine_multiscale.h',
+    'include/gridfire/engine/views/engine_priming.h',
+    'include/gridfire/engine/procedures/priming.h',
     'include/gridfire/reaction/reaction.h',
     'include/gridfire/reaction/reaclib.h',
     'include/gridfire/io/network_file.h',
@@ -66,5 +72,6 @@ network_headers = files(
     'include/gridfire/partition/partition_ground.h',
     'include/gridfire/partition/composite/partition_composite.h',
     'include/gridfire/utils/logging.h',
+    'include/gridfire/utils/qse_rules.h',
 )
 install_headers(network_headers, subdir : 'gridfire')

subprojects/libcomposition.wrap

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

tests/graphnet_sandbox/main.cpp

@@ -6,6 +6,8 @@
 #include "gridfire/engine/views/engine_adaptive.h"
 #include "gridfire/partition/partition_types.h"
 #include "gridfire/engine/views/engine_defined.h"
+#include "gridfire/engine/views/engine_multiscale.h"
+#include "gridfire/engine/procedures/priming.h"
 #include "gridfire/io/network_file.h"
 
 #include "gridfire/solver/solver.h"
@@ -58,17 +60,25 @@ 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::TraceL2);
+    logger->set_log_level(quill::LogLevel::TraceL1);
     LOG_DEBUG(logger, "Starting Adaptive Engine View Example...");
 
     using namespace gridfire;
     const std::vector<double> comp = {0.708, 0.0, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4};
     const std::vector<std::string> symbols = {"H-1", "H-2", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"};
 
+
     fourdst::composition::Composition composition;
     composition.registerSymbol(symbols, true);
     composition.setMassFraction(symbols, comp);
     composition.finalize(true);
+    using partition::BasePartitionType;
+    const auto partitionFunction = partition::CompositePartitionFunction({
+        BasePartitionType::RauscherThielemann,
+        BasePartitionType::GroundState
+    });
+    GraphEngine ReaclibEngine(composition, partitionFunction);
+
 
 
     NetIn netIn;
@@ -77,38 +87,30 @@ int main() {
     netIn.density = 1e2;
     netIn.energy = 0.0;
 
-    netIn.tMax = 3.1536e17;
+    netIn.tMax = 1; // 1 year in seconds
     netIn.dt0 = 1e12;
 
     NetOut netOut;
-
-    // netIn.dt0 = 1e12;
-    // approx8::Approx8Network approx8Network;
-    // measure_execution_time([&]() {
-    //     netOut = approx8Network.evaluate(netIn);
-    // }, "Approx8 Network Initialization");
-    // std::cout << "Approx8 Network H-1: " << netOut.composition.getMassFraction("H-1") << " in " << netOut.num_steps << " steps." << std::endl;
-
-    using partition::BasePartitionType;
-    const auto partitionFunction = partition::CompositePartitionFunction({
-        BasePartitionType::RauscherThielemann,
-        BasePartitionType::GroundState
-    });
-    std::cout << "Partition Function for Mg-24: " << partitionFunction.evaluate(12, 24, 8) << std::endl;
-    std::cout << "Partition Function for F-23: " << partitionFunction.evaluate(9, 23, 8) << std::endl;
-    std::cout << "Partition Function for O-13: " << partitionFunction.evaluate(8, 13, 8) << std::endl;
-
     netIn.dt0 = 1e-15;
 
-    GraphEngine ReaclibEngine(composition, partitionFunction);
-    std::cout << ReaclibEngine.getPartitionFunction().type() << std::endl;
+    const auto primedComposition = primeNetwork(netIn, ReaclibEngine);
+    netIn.composition = primedComposition;
+    netIn.composition.finalize(true);
+    MultiscalePartitioningEngineView multiscaleEngine(ReaclibEngine);
+    multiscaleEngine.equilibrateNetwork(netIn, 1e12);
+
+
+    // // std::cout << ReaclibEngine.getPartitionFunction().type() << std::endl;
     // ReaclibEngine.setPrecomputation(true);
-    // // AdaptiveEngineView adaptiveEngine(ReaclibEngine);
+    // ReaclibEngine.setUseReverseReactions(false);
+    // // // AdaptiveEngineView adaptiveEngine(ReaclibEngine);
     // io::SimpleReactionListFileParser parser{};
     // FileDefinedEngineView approx8EngineView(ReaclibEngine, "approx8.net", parser);
     // approx8EngineView.setScreeningModel(screening::ScreeningType::WEAK);
-    // solver::QSENetworkSolver solver(approx8EngineView);
+    // solver::DirectNetworkSolver solver(approx8EngineView);
     // netOut = solver.evaluate(netIn);
+    // std::cout << netOut.composition << std::endl;
+
 
     // measure_execution_time([&]() {
     //     netOut = solver.evaluate(netIn);