feat(GridFire): Much more robust starting network

GraphEngine now can initialize with a much more robust set of reactions (including the entire reaction set). The jacobian can still be efficiently evaluated using CppAD's sparse jacobian feature. Further, the primeing network has been signifiganty enhanced to handle much hotter termperatures

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

@@ -5,8 +5,12 @@
 #include "gridfire/screening/screening_abstract.h"
 #include "gridfire/screening/screening_types.h"
 
+#include "gridfire/engine/types/reporting.h"
+#include "gridfire/engine/types/building.h"
+
 #include <vector>
 #include <unordered_map>
+#include <utility>
 
 /**
  * @file engine_abstract.h
@@ -55,6 +59,8 @@ namespace gridfire {
         T nuclearEnergyGenerationRate = T(0.0); ///< Specific energy generation rate (e.g., erg/g/s).
     };
 
+    using SparsityPattern = std::vector<std::pair<size_t, size_t>>;
+
     /**
      * @brief Abstract base class for a reaction network engine.
      *
@@ -132,9 +138,19 @@ namespace gridfire {
          */
         virtual void generateJacobianMatrix(
             const std::vector<double>& Y_dynamic,
-            double T9, double rho
+            double T9,
+            double rho
         ) = 0;
 
+        virtual void generateJacobianMatrix(
+            const std::vector<double>& Y_dynamic,
+            double T9,
+            double rho,
+            const SparsityPattern& sparsityPattern
+        ) {
+            throw std::logic_error("Sparsity pattern not supported by this engine.");
+        }
+
         /**
          * @brief Get an entry from the previously generated Jacobian matrix.
          *
@@ -149,6 +165,7 @@ namespace gridfire {
             int j
         ) const = 0;
 
+
         /**
          * @brief Generate the stoichiometry matrix for the network.
          *
@@ -269,5 +286,15 @@ namespace gridfire {
         [[nodiscard]] virtual int getSpeciesIndex(const fourdst::atomic::Species &species) const = 0;
 
         [[nodiscard]] virtual std::vector<double> mapNetInToMolarAbundanceVector(const NetIn &netIn) const = 0;
+
+        [[nodiscard]] virtual PrimingReport primeEngine(const NetIn &netIn) = 0;
+
+        [[nodiscard]] virtual BuildDepthType getDepth() const {
+            throw std::logic_error("Network depth not supported by this engine.");
+        }
+
+        virtual void rebuild(const fourdst::composition::Composition& comp, BuildDepthType depth) {
+            throw std::logic_error("Setting network depth not supported by this engine.");
+        }
     };
 }

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

@@ -11,7 +11,7 @@
 #include "gridfire/screening/screening_abstract.h"
 #include "gridfire/screening/screening_types.h"
 #include "gridfire/partition/partition_abstract.h"
-#include "gridfire/partition/partition_ground.h"
+#include "gridfire/engine/procedures/construction.h"
 
 #include <string>
 #include <unordered_map>
@@ -22,6 +22,12 @@
 #include <boost/numeric/ublas/matrix_sparse.hpp>
 
 #include "cppad/cppad.hpp"
+#include "cppad/utility/sparse_rc.hpp"
+#include "cppad/speed/sparse_jac_fun.hpp"
+
+#include "procedures/priming.h"
+
+
 #include "quill/LogMacros.h"
 
 // PERF: The function getNetReactionStoichiometry returns a map of species to their stoichiometric coefficients for a given reaction.
@@ -30,6 +36,7 @@
 //       REACLIBReactions are quite large data structures, so this could be a performance bottleneck.
 
 namespace gridfire {
+    static bool s_debug = false; // Global debug flag for the GraphEngine
     /**
      * @brief Alias for CppAD AD type for double precision.
      *
@@ -103,10 +110,15 @@ namespace gridfire {
          *
          * @see build_reaclib_nuclear_network
          */
-        explicit GraphEngine(const fourdst::composition::Composition &composition);
+        explicit GraphEngine(
+            const fourdst::composition::Composition &composition,
+            const BuildDepthType = NetworkBuildDepth::Full
+        );
 
-        explicit GraphEngine(const fourdst::composition::Composition &composition,
+        explicit GraphEngine(
-                             const partition::PartitionFunction& partitionFunction
+            const fourdst::composition::Composition &composition,
+            const partition::PartitionFunction& partitionFunction,
+            const BuildDepthType buildDepth = NetworkBuildDepth::Full
         );
 
         /**
@@ -157,6 +169,13 @@ namespace gridfire {
             const double rho
         ) override;
 
+        void generateJacobianMatrix(
+            const std::vector<double> &Y_dynamic,
+            double T9,
+            double rho,
+            const SparsityPattern &sparsityPattern
+        ) override;
+
         /**
          * @brief Generates the stoichiometry matrix for the network.
          *
@@ -322,28 +341,34 @@ namespace gridfire {
             double T9
         ) const;
 
-        double calculateReverseRateTwoBody(
+        [[nodiscard]] double calculateReverseRateTwoBody(
             const reaction::Reaction &reaction,
             const double T9,
             const double forwardRate,
             const double expFactor
         ) const;
 
-        double calculateReverseRateTwoBodyDerivative(
+        [[nodiscard]] double calculateReverseRateTwoBodyDerivative(
             const reaction::Reaction &reaction,
             const double T9,
             const double reverseRate
         ) const;
 
-        bool isUsingReverseReactions() const;
+        [[nodiscard]] bool isUsingReverseReactions() const;
 
         void setUseReverseReactions(bool useReverse);
 
-        int getSpeciesIndex(
+        [[nodiscard]] int getSpeciesIndex(
             const fourdst::atomic::Species& species
         ) const override;
 
-        std::vector<double> mapNetInToMolarAbundanceVector(const NetIn &netIn) const override;
+        [[nodiscard]] std::vector<double> mapNetInToMolarAbundanceVector(const NetIn &netIn) const override;
+
+        [[nodiscard]] PrimingReport primeEngine(const NetIn &netIn) override;
+
+        [[nodiscard]] BuildDepthType getDepth() const override;
+
+        void rebuild(const fourdst::composition::Composition& comp, const BuildDepthType depth) override;
 
 
 
@@ -427,6 +452,8 @@ 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.
+        CppAD::sparse_jac_work m_jac_work; ///< Work object for sparse Jacobian calculations.
+        CppAD::sparse_rc<std::vector<size_t>> m_full_jacobian_sparsity_pattern; ///< Full sparsity pattern for the Jacobian matrix.
         std::vector<std::unique_ptr<AtomicReverseRate>> m_atomicReverseRates;
 
         screening::ScreeningType m_screeningType = screening::ScreeningType::BARE; ///< Screening type for the reaction network. Default to no screening.
@@ -436,6 +463,8 @@ namespace gridfire {
 
         bool m_useReverseReactions = true; ///< Flag to enable or disable reverse reactions. If false, only forward reactions are considered.
 
+        BuildDepthType m_depth;
+
         std::vector<PrecomputedReaction> m_precomputedReactions; ///< Precomputed reactions for efficiency.
         std::unique_ptr<partition::PartitionFunction> m_partitionFunction; ///< Partition function for the network.
 
@@ -636,6 +665,14 @@ namespace gridfire {
         if (!m_useReverseReactions) {
             return static_cast<T>(0.0); // If reverse reactions are not used, return zero
         }
+        s_debug = false;
+        if (reaction.peName() == "p(p,e+)d" || reaction.peName() =="d(d,n)he3" || reaction.peName() == "c12(p,g)n13") {
+            if constexpr (!std::is_same_v<T, ADDouble>) {
+                s_debug = true;
+                std::cout << "Calculating reverse molar flow for reaction: " << reaction.peName() << std::endl;
+                std::cout << "\tT9: " << T9 << ", rho: " << rho << std::endl;
+            }
+        }
         T reverseMolarFlow = static_cast<T>(0.0);
 
         if (reaction.qValue() != 0.0) {
@@ -655,6 +692,9 @@ namespace gridfire {
                 }
             } else {
                 // A,B If not calling with an AD type, calculate the reverse rate directly
+                if (s_debug) {
+                    std::cout << "\tUsing double overload\n";
+                }
                 reverseRateConstant = calculateReverseRate(reaction, T9);
             }
 
@@ -673,7 +713,7 @@ namespace gridfire {
             // 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);
+                const unsigned long speciesIndex = m_speciesToIndexMap.at(species);
                 productAbundanceTerm *= CppAD::pow(Y[speciesIndex], count);
             }
 
@@ -757,7 +797,7 @@ namespace gridfire {
             ss << "Forward: " << forwardMolarReactionFlow
                << ", Reverse: " << reverseMolarFlow
                << ", Net: " << molarReactionFlow;
-            LOG_DEBUG(
+            LOG_TRACE_L2(
                 m_logger,
                 "Reaction: {}, {}",
                 reaction.peName(),

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

@@ -1,8 +1,17 @@
-//
+#pragma once
-// Created by Emily Boudreaux on 7/14/25.
-//
 
-#ifndef CONSTRUCTION_H
+#include "gridfire/reaction/reaction.h"
-#define CONSTRUCTION_H
+#include "gridfire/engine/types/building.h"
 
-#endif //CONSTRUCTION_H
+#include "fourdst/composition/composition.h"
+
+#include <variant>
+
+namespace gridfire {
+
+    reaction::LogicalReactionSet build_reaclib_nuclear_network(
+        const fourdst::composition::Composition &composition,
+        BuildDepthType maxLayers = NetworkBuildDepth::Full,
+        bool reverse = false
+    );
+}

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

@@ -6,9 +6,15 @@
 #include "fourdst/composition/composition.h"
 #include "fourdst/composition/atomicSpecies.h"
 
+#include <map>
+#include <ranges>
+#include <sstream>
+
 
 namespace gridfire {
-    fourdst::composition::Composition primeNetwork(
+
+
+    PrimingReport primeNetwork(
         const NetIn&,
         DynamicEngine& engine
     );

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

@@ -0,0 +1,16 @@
+#pragma once
+
+#include <variant>
+
+namespace gridfire {
+    enum class NetworkBuildDepth {
+        Full = -1,
+        Shallow = 1,
+        SecondOrder = 2,
+        ThirdOrder = 3,
+        FourthOrder = 4,
+        FifthOrder = 5
+    };
+
+    using BuildDepthType = std::variant<NetworkBuildDepth, int>;
+}

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

@@ -1,8 +1,42 @@
-//
+#pragma once
-// Created by Emily Boudreaux on 7/14/25.
-//
 
-#ifndef REPORTING_H
+#include <map>
-#define REPORTING_H
+#include <string>
+#include <ranges>
 
-#endif //REPORTING_H
+namespace gridfire {
+    enum class PrimingReportStatus {
+        FULL_SUCCESS = 0,
+        NO_SPECIES_TO_PRIME = 1,
+        MAX_ITERATIONS_REACHED = 2,
+        FAILED_TO_FINALIZE_COMPOSITION = 3,
+        FAILED_TO_FIND_CREATION_CHANNEL = 4,
+        FAILED_TO_FIND_PRIMING_REACTIONS = 5,
+        BASE_NETWORK_TOO_SHALLOW = 6
+    };
+
+    inline std::map<PrimingReportStatus, std::string> PrimingReportStatusStrings = {
+        {PrimingReportStatus::FULL_SUCCESS, "Full Success"},
+        {PrimingReportStatus::NO_SPECIES_TO_PRIME, "No Species to Prime"},
+        {PrimingReportStatus::MAX_ITERATIONS_REACHED, "Max Iterations Reached"},
+        {PrimingReportStatus::FAILED_TO_FINALIZE_COMPOSITION, "Failed to Finalize Composition"},
+        {PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL, "Failed to Find Creation Channel"},
+        {PrimingReportStatus::FAILED_TO_FIND_PRIMING_REACTIONS, "Failed to Find Priming Reactions"},
+        {PrimingReportStatus::BASE_NETWORK_TOO_SHALLOW, "Base Network Too Shallow"}
+    };
+
+    struct PrimingReport {
+        fourdst::composition::Composition primedComposition;
+        std::vector<std::pair<fourdst::atomic::Species, double>> massFractionChanges; ///< Species and their destruction/creation rates
+        bool success;
+        PrimingReportStatus status;
+
+        friend std::ostream& operator<<(std::ostream& os, const PrimingReport& report) {
+            std::stringstream ss;
+            std::string successStr = report.success ? "true" : "false";
+            ss << "PrimingReport(success=" << successStr
+               << ", status=" << PrimingReportStatusStrings[report.status] << ")";
+            return os << ss.str();
+        }
+    };
+}

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

@@ -9,6 +9,9 @@
 #include "fourdst/config/config.h"
 #include "fourdst/logging/logging.h"
 
+#include "gridfire/engine/procedures/priming.h"
+#include "gridfire/engine/procedures/construction.h"
+
 #include "quill/Logger.h"
 
 namespace gridfire {
@@ -260,6 +263,8 @@ namespace gridfire {
         [[nodiscard]] int getSpeciesIndex(const fourdst::atomic::Species &species) const override;
 
         [[nodiscard]] std::vector<double> mapNetInToMolarAbundanceVector(const NetIn &netIn) const override;
+
+        [[nodiscard]] PrimingReport primeEngine(const NetIn &netIn) override;
     private:
         using Config = fourdst::config::Config;
         using LogManager = fourdst::logging::LogManager;

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

@@ -161,6 +161,8 @@ namespace gridfire{
         [[nodiscard]] int getSpeciesIndex(const fourdst::atomic::Species &species) const override;
 
         [[nodiscard]] std::vector<double> mapNetInToMolarAbundanceVector(const NetIn &netIn) const override;
+
+        [[nodiscard]] PrimingReport primeEngine(const NetIn &netIn) override;
     protected:
         bool m_isStale = true;
         DynamicEngine& m_baseEngine;

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

@@ -83,18 +83,20 @@ namespace gridfire {
 
         [[nodiscard]] std::vector<double> mapNetInToMolarAbundanceVector(const NetIn &netIn) const override;
 
-        std::vector<fourdst::atomic::Species> getFastSpecies() const;
+        [[nodiscard]] PrimingReport primeEngine(const NetIn &netIn) override;
-        const std::vector<fourdst::atomic::Species>& getDynamicSpecies() const;
 
-        void equilibrateNetwork(
+        [[nodiscard]] std::vector<fourdst::atomic::Species> getFastSpecies() const;
-            const std::vector<double>& Y,
+        [[nodiscard]] const std::vector<fourdst::atomic::Species>& getDynamicSpecies() const;
+
+        fourdst::composition::Composition equilibrateNetwork(
+            const std::vector<double> &Y,
             double T9,
             double rho,
             double dt_control
         );
 
-        void equilibrateNetwork(
+        fourdst::composition::Composition equilibrateNetwork(
-            const NetIn& netIn,
+            const NetIn &netIn,
             const double dt_control
         );
 
@@ -137,8 +139,8 @@ namespace gridfire {
             m_rho(rho),
             m_Y_scale(Y_scale) {}
 
-            int values() const { return m_qse_solve_indices.size(); }
+            [[nodiscard]] int values() const { return m_qse_solve_indices.size(); }
-            int inputs() const { return m_qse_solve_indices.size(); }
+            [[nodiscard]] 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;

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

@@ -20,13 +20,10 @@ namespace gridfire {
     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(

src/network/include/gridfire/network.h

@@ -104,7 +104,5 @@ namespace gridfire {
     };
 
 
-    reaction::LogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition, bool reverse);
-
 
 } // namespace nuclearNetwork

src/network/lib/engine/engine_graph.cpp

@@ -2,6 +2,9 @@
 #include "gridfire/reaction/reaction.h"
 #include "gridfire/network.h"
 #include "gridfire/screening/screening_types.h"
+#include "gridfire/engine/procedures/priming.h"
+#include "gridfire/partition/partition_ground.h"
+#include "gridfire/engine/procedures/construction.h"
 
 #include "fourdst/composition/species.h"
 #include "fourdst/composition/atomicSpecies.h"
@@ -22,17 +25,24 @@
 
 #include <boost/numeric/odeint.hpp>
 
+#include "cppad/cppad.hpp"
+#include "cppad/utility/sparse_rc.hpp"
+#include "cppad/utility/sparse_rcv.hpp"
+
 
 namespace gridfire {
     GraphEngine::GraphEngine(
-        const fourdst::composition::Composition &composition
+        const fourdst::composition::Composition &composition,
-    ): GraphEngine(composition, partition::GroundStatePartitionFunction()) {}
+        const BuildDepthType buildDepth
+    ): GraphEngine(composition, partition::GroundStatePartitionFunction(), buildDepth) {}
 
     GraphEngine::GraphEngine(
         const fourdst::composition::Composition &composition,
-        const partition::PartitionFunction& partitionFunction) :
+        const partition::PartitionFunction& partitionFunction,
-    m_reactions(build_reaclib_nuclear_network(composition, false)),
+        const BuildDepthType buildDepth) :
-    m_partitionFunction(partitionFunction.clone())  // Clone the partition function to ensure ownership
+    m_reactions(build_reaclib_nuclear_network(composition, buildDepth, false)),
+    m_partitionFunction(partitionFunction.clone()),
+    m_depth(buildDepth)
     {
         syncInternalMaps();
     }
@@ -77,6 +87,16 @@ namespace gridfire {
         generateStoichiometryMatrix();
         reserveJacobianMatrix();
         recordADTape();
+
+        const size_t n = m_rhsADFun.Domain();
+        const size_t m = m_rhsADFun.Range();
+
+        std::vector<bool> select_domain(n, true);
+        std::vector<bool> select_range(m, true);
+
+        m_rhsADFun.subgraph_sparsity(select_domain, select_range, false, m_full_jacobian_sparsity_pattern);
+        m_jac_work.clear();
+
         precomputeNetwork();
         LOG_INFO(m_logger, "Internal maps synchronized. Network contains {} species and {} reactions.",
                  m_networkSpecies.size(), m_reactions.size());
@@ -243,6 +263,14 @@ namespace gridfire {
             return 0.0; // If reverse reactions are not used, return 0.0
         }
         const double expFactor = std::exp(-reaction.qValue() / (m_constants.kB * T9));
+        if (s_debug) {
+            std::cout << "\texpFactor = exp(-qValue/(kB * T9))\n";
+            std::cout << "\texpFactor: " << expFactor << " for reaction: " << reaction.peName() << std::endl;
+            std::cout << "\tQ-value: " << reaction.qValue() << " at T9: " << T9 << std::endl;
+            std::cout << "\tT9: " << T9 << std::endl;
+            std::cout << "\tkB * T9: " << m_constants.kB * T9 << std::endl;
+            std::cout << "\tqValue/(kB * T9): " << reaction.qValue() / (m_constants.kB * T9) << std::endl;
+        }
         double reverseRate = 0.0;
         const double forwardRate = reaction.calculate_rate(T9);
 
@@ -393,6 +421,47 @@ namespace gridfire {
         return Y; // Return the vector of molar abundances
     }
 
+    PrimingReport GraphEngine::primeEngine(const NetIn &netIn) {
+        NetIn fullNetIn;
+        fourdst::composition::Composition composition;
+
+        std::vector<std::string> symbols;
+        symbols.reserve(m_networkSpecies.size());
+        for (const auto &symbol: m_networkSpecies) {
+            symbols.emplace_back(symbol.name());
+        }
+        composition.registerSymbol(symbols);
+        for (const auto& [symbol, entry] : netIn.composition) {
+            if (m_networkSpeciesMap.contains(symbol)) {
+                composition.setMassFraction(symbol, entry.mass_fraction());
+            } else {
+                composition.setMassFraction(symbol, 0.0);
+            }
+        }
+        composition.finalize(true);
+        fullNetIn.composition = composition;
+        fullNetIn.temperature = netIn.temperature;
+        fullNetIn.density = netIn.density;
+
+        auto primingReport = primeNetwork(fullNetIn, *this);
+
+        return primingReport;
+    }
+
+    BuildDepthType GraphEngine::getDepth() const {
+        return m_depth;
+    }
+
+    void GraphEngine::rebuild(const fourdst::composition::Composition& comp, const BuildDepthType depth) {
+        if (depth != m_depth) {
+            m_depth = depth;
+            m_reactions = build_reaclib_nuclear_network(comp, m_depth, false);
+            syncInternalMaps(); // Resync internal maps after changing the depth
+        } else {
+            LOG_DEBUG(m_logger, "Rebuild requested with the same depth. No changes made to the network.");
+        }
+    }
+
     StepDerivatives<double> GraphEngine::calculateAllDerivativesUsingPrecomputation(
         const std::vector<double> &Y_in,
         const std::vector<double> &bare_rates,
@@ -414,23 +483,23 @@ namespace gridfire {
         molarReactionFlows.reserve(m_precomputedReactions.size());
 
         for (const auto& precomp : m_precomputedReactions) {
-            double abundanceProduct = 1.0;
+            double forwardAbundanceProduct = 1.0;
-            bool below_threshold = false;
+            // bool below_threshold = false;
             for (size_t i = 0; i < precomp.unique_reactant_indices.size(); ++i) {
                 const size_t reactantIndex = precomp.unique_reactant_indices[i];
                 const int power = precomp.reactant_powers[i];
-                const double abundance = Y_in[reactantIndex];
+                // const double abundance = Y_in[reactantIndex];
-                if (abundance < MIN_ABUNDANCE_THRESHOLD) {
+                // if (abundance < MIN_ABUNDANCE_THRESHOLD) {
-                    below_threshold = true;
+                //     below_threshold = true;
-                    break;
+                //     break;
-                }
+                // }
 
-                abundanceProduct *= std::pow(Y_in[reactantIndex], power);
+                forwardAbundanceProduct *= std::pow(Y_in[reactantIndex], power);
-            }
-            if (below_threshold) {
-                molarReactionFlows.push_back(0.0);
-                continue; // Skip this reaction if any reactant is below the abundance threshold
             }
+            // if (below_threshold) {
+            //     molarReactionFlows.push_back(0.0);
+            //     continue; // Skip this reaction if any reactant is below the abundance threshold
+            // }
 
             const double bare_rate = bare_rates[precomp.reaction_index];
             const double screeningFactor = screeningFactors[precomp.reaction_index];
@@ -441,8 +510,8 @@ namespace gridfire {
                     screeningFactor *
                     bare_rate *
                     precomp.symmetry_factor *
-                    abundanceProduct *
+                    forwardAbundanceProduct *
-                    std::pow(rho, numReactants >  1 ? numReactants - 1 : 0);
+                    std::pow(rho, numReactants >  1 ? numReactants - 1 : 0.0);
 
             double reverseMolarReactionFlow = 0.0;
             if (precomp.reverse_symmetry_factor != 0.0 and m_useReverseReactions) {
@@ -455,7 +524,7 @@ namespace gridfire {
                     bare_reverse_rate *
                     precomp.reverse_symmetry_factor *
                     reverseAbundanceProduct *
-                    std::pow(rho, numProducts > 1 ? numProducts - 1 : 0);
+                    std::pow(rho, numProducts > 1 ? numProducts - 1 : 0.0);
             }
 
             molarReactionFlows.push_back(forwardMolarReactionFlow - reverseMolarReactionFlow);
@@ -474,7 +543,7 @@ namespace gridfire {
                 const int stoichiometricCoefficient = precomp.stoichiometric_coefficients[i];
 
                 // Update the derivative for this species
-                result.dydt[speciesIndex] += static_cast<double>(stoichiometricCoefficient) * R_j / rho;
+                result.dydt[speciesIndex] += static_cast<double>(stoichiometricCoefficient) * R_j;
             }
         }
 
@@ -641,6 +710,63 @@ namespace gridfire {
         LOG_TRACE_L1(m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2());
     }
 
+    void GraphEngine::generateJacobianMatrix(
+        const std::vector<double> &Y_dynamic,
+        const double T9,
+        const double rho,
+        const SparsityPattern &sparsityPattern
+    ) {
+        // --- Pack the input variables into a vector for CppAD ---
+        const size_t numSpecies = m_networkSpecies.size();
+        std::vector<double> x(numSpecies + 2, 0.0);
+        for (size_t i = 0; i < numSpecies; ++i) {
+           x[i] = Y_dynamic[i];
+        }
+        x[numSpecies] = T9;
+        x[numSpecies + 1] = rho;
+
+        // --- Convert into CppAD Sparsity pattern ---
+        const size_t nnz = sparsityPattern.size(); // Number of non-zero entries in the sparsity pattern
+        std::vector<size_t> row_indices(nnz);
+        std::vector<size_t> col_indices(nnz);
+
+        for (size_t k = 0; k < nnz; ++k) {
+            row_indices[k] = sparsityPattern[k].first;
+            col_indices[k] = sparsityPattern[k].second;
+        }
+
+        std::vector<double> values(nnz);
+        const size_t num_rows_jac = numSpecies;
+        const size_t num_cols_jac = numSpecies + 2; // +2 for T9 and rho
+
+        CppAD::sparse_rc<std::vector<size_t>> CppAD_sparsity_pattern(num_rows_jac, num_cols_jac, nnz);
+        for (size_t k = 0; k < nnz; ++k) {
+            CppAD_sparsity_pattern.set(k, sparsityPattern[k].first, sparsityPattern[k].second);
+        }
+
+        CppAD::sparse_rcv<std::vector<size_t>, std::vector<double>> jac_subset(CppAD_sparsity_pattern);
+
+        m_rhsADFun.sparse_jac_rev(
+            x,
+            jac_subset, // Sparse Jacobian output
+            m_full_jacobian_sparsity_pattern,
+            "cppad",
+            m_jac_work // Work vector for CppAD
+        );
+
+        // --- Convert the sparse Jacobian back to the Boost uBLAS format ---
+        m_jacobianMatrix.clear();
+        for (size_t k = 0; k < nnz; ++k) {
+            const size_t row = jac_subset.row()[k];
+            const size_t col = jac_subset.col()[k];
+            const double value = jac_subset.val()[k];
+
+            if (std::abs(value) > MIN_JACOBIAN_THRESHOLD) {
+                m_jacobianMatrix(row, col) = value; // Insert into the sparse matrix
+            }
+        }
+    }
+
     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);

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

@@ -1,3 +1,113 @@
-//
+#include "gridfire/engine/procedures/construction.h"
-// Created by Emily Boudreaux on 7/14/25.
+
-//
+#include <ranges>
+#include <stdexcept>
+
+#include "gridfire/reaction/reaction.h"
+#include "gridfire/reaction/reaclib.h"
+
+#include "fourdst/composition/composition.h"
+
+#include "fourdst/logging/logging.h"
+
+#include "quill/Logger.h"
+#include "quill/LogMacros.h"
+
+namespace gridfire {
+    using reaction::LogicalReactionSet;
+    using reaction::ReactionSet;
+    using reaction::Reaction;
+    using fourdst::composition::Composition;
+    using fourdst::atomic::Species;
+
+
+    LogicalReactionSet build_reaclib_nuclear_network(
+        const Composition &composition,
+        BuildDepthType maxLayers,
+        bool reverse
+    ) {
+        int depth;
+        if (std::holds_alternative<NetworkBuildDepth>(maxLayers)) {
+            depth = static_cast<int>(std::get<NetworkBuildDepth>(maxLayers));
+        } else {
+            depth = std::get<int>(maxLayers);
+        }
+        auto logger = fourdst::logging::LogManager::getInstance().getLogger("log");
+        if (depth == 0) {
+            LOG_ERROR(logger, "Network build depth is set to 0. No reactions will be collected.");
+            throw std::logic_error("Network build depth is set to 0. No reactions will be collected.");
+        }
+
+        const auto allReactions = reaclib::get_all_reactions();
+        std::vector<Reaction> remainingReactions;
+        for (const auto& reaction : allReactions) {
+            if (reaction.is_reverse() == reverse) {
+                remainingReactions.push_back(reaction);
+            }
+        }
+
+        if (depth == static_cast<int>(NetworkBuildDepth::Full)) {
+            LOG_INFO(logger, "Building full nuclear network with a total of {} reactions.", allReactions.size());
+            const ReactionSet reactionSet(remainingReactions);
+            return reaction::packReactionSetToLogicalReactionSet(reactionSet);
+        }
+
+        std::unordered_set<Species> availableSpecies;
+        for (const auto &entry: composition | std::views::values) {
+            if (entry.mass_fraction() > 0.0) {
+                availableSpecies.insert(entry.isotope());
+            }
+        }
+
+
+        std::vector<Reaction> collectedReactions;
+
+        LOG_INFO(logger, "Starting network construction with {} available species.", availableSpecies.size());
+        for (int layer = 0; layer < depth && !remainingReactions.empty(); ++layer) {
+            LOG_TRACE_L1(logger, "Collecting reactions for layer {} with {} remaining reactions. Currently there are {} available species", layer, remainingReactions.size(), availableSpecies.size());
+            std::vector<Reaction> reactionsForNextPass;
+            std::unordered_set<Species> newProductsThisLayer;
+            bool newReactionsAdded = false;
+
+            reactionsForNextPass.reserve(remainingReactions.size());
+
+            for (const auto &reaction : remainingReactions) {
+                bool allReactantsAvailable = true;
+                for (const auto& reactant : reaction.reactants()) {
+                    if (!availableSpecies.contains(reactant)) {
+                        allReactantsAvailable = false;
+                        break;
+                    }
+                }
+
+                if (allReactantsAvailable) {
+                    collectedReactions.push_back(reaction);
+                    newReactionsAdded = true;
+
+                    for (const auto& product : reaction.products()) {
+                        newProductsThisLayer.insert(product);
+                    }
+                } else {
+                    reactionsForNextPass.push_back(reaction);
+                }
+            }
+
+            if (!newReactionsAdded) {
+                LOG_INFO(logger, "No new reactions added in layer {}. Stopping network construction with {} reactions collected.", layer, collectedReactions.size());
+                break;
+            }
+
+            LOG_TRACE_L1(logger, "Layer {}: Collected {} reactions. New products this layer: {}", layer, collectedReactions.size(), newProductsThisLayer.size());
+            availableSpecies.insert(newProductsThisLayer.begin(), newProductsThisLayer.end());
+
+            remainingReactions = std::move(reactionsForNextPass);
+        }
+
+        LOG_INFO(logger, "Network construction completed with {} reactions collected.", collectedReactions.size());
+        const ReactionSet reactionSet(collectedReactions);
+        return reaction::packReactionSetToLogicalReactionSet(reactionSet);
+    }
+
+
+
+}

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

@@ -1,48 +1,189 @@
 #include "gridfire/engine/procedures/priming.h"
 #include "gridfire/engine/views/engine_priming.h"
+#include "gridfire/engine/procedures/construction.h"
 #include "gridfire/solver/solver.h"
 
 #include "gridfire/engine/engine_abstract.h"
 #include "gridfire/network.h"
 
+#include "fourdst/logging/logging.h"
+#include "quill/Logger.h"
+#include "quill/LogMacros.h"
+
 namespace gridfire {
-    fourdst::composition::Composition primeNetwork(const NetIn& netIn, DynamicEngine& engine) {
     using fourdst::composition::Composition;
     using fourdst::atomic::Species;
 
+    const reaction::Reaction* findDominantCreationChannel (
+        const DynamicEngine& engine,
+        const Species& species,
+        const std::vector<double>& Y,
+        const double T9,
+        const double rho
+    ) {
+        const reaction::Reaction* dominateReaction = nullptr;
+        double maxFlow = -1.0;
+        for (const auto& reaction : engine.getNetworkReactions()) {
+            if (reaction.contains(species) && reaction.stoichiometry(species) > 0) {
+                const double flow = engine.calculateMolarReactionFlow(reaction, Y, T9, rho);
+                if (flow > maxFlow) {
+                    maxFlow = flow;
+                    dominateReaction = &reaction;
+                }
+            }
+        }
+        return dominateReaction;
+    }
+
+
+    PrimingReport primeNetwork(const NetIn& netIn, DynamicEngine& engine) {
+        auto logger = LogManager::getInstance().getLogger("log");
 
-        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());
             }
         }
+        LOG_INFO(logger, "Priming {} species in the network.", speciesToPrime.size());
 
-
+        PrimingReport report;
         if (speciesToPrime.empty()) {
-            return netIn.composition; // No priming needed
+            report.primedComposition = netIn.composition;
+            report.success = true;
+            report.status = PrimingReportStatus::NO_SPECIES_TO_PRIME;
+            return report;
         }
 
-        const double T9 = netIn.temperature / 1e9; // Convert temperature to T9
+        const double T9 = netIn.temperature / 1e9;
-        const double rho = netIn.density; // Density in g/cm^3
+        const double rho = netIn.density;
+        const auto initialDepth = engine.getDepth();
+
+        report.status = PrimingReportStatus::FULL_SUCCESS;
+        report.success = true;
+
+        // --- 1: pack composition into internal map ---
+        std::unordered_map<Species, double> currentMassFractions;
+        for (const auto& entry : netIn.composition | std::views::values) {
+            currentMassFractions[entry.isotope()] = entry.mass_fraction();
+        }
+        for (const auto& entry : speciesToPrime) {
+            currentMassFractions[entry] = 0.0; // Initialize priming species with 0 mass fraction
+        }
+
+        std::unordered_map<Species, double> totalMassFractionChanges;
+
+        engine.rebuild(netIn.composition, NetworkBuildDepth::Full);
+
+        for (const auto& primingSpecies : speciesToPrime) {
+            LOG_DEBUG(logger, "Priming species: {}", primingSpecies.name());
+
+            // Create a temporary composition from the current internal state for the primer
+            Composition tempComp;
+            for(const auto& [sp, mf] : currentMassFractions) {
+                tempComp.registerSymbol(std::string(sp.name()));
+                if (mf < 0.0 && std::abs(mf) < 1e-16) {
+                    tempComp.setMassFraction(sp, 0.0); // Avoid negative mass fractions
+                } else {
+                    tempComp.setMassFraction(sp, mf);
+                }
+            }
+            tempComp.finalize(true); // Finalize with normalization
+
+            NetIn tempNetIn = netIn;
+            tempNetIn.composition = tempComp;
 
-        for (const auto& primingSpecies :speciesToPrime) {
             NetworkPrimingEngineView primer(primingSpecies, engine);
 
-            const auto Y = primer.mapNetInToMolarAbundanceVector(currentConditions);
+            if (primer.getNetworkReactions().size() == 0) {
-            const double destructionRateConstant = calculateDestructionRateConstant(primer, primingSpecies, Y, T9, rho);
+                LOG_ERROR(logger, "No priming reactions found for species {}.", primingSpecies.name());
-            const double creationRate = calculateCreationRate(primer, primingSpecies, Y, T9, rho);
+                report.success = false;
-            const double equilibriumMolarAbundance = creationRate / destructionRateConstant;
+                report.status = PrimingReportStatus::FAILED_TO_FIND_PRIMING_REACTIONS;
-            const double equilibriumMassFraction = equilibriumMolarAbundance * primingSpecies.mass();
+                continue;
-
-            currentConditions.composition.setMassFraction(
-                std::string(primingSpecies.name()),
-                equilibriumMassFraction
-            );
             }
 
-        return currentConditions.composition;
+            const auto Y = primer.mapNetInToMolarAbundanceVector(tempNetIn);
+            const double destructionRateConstant = calculateDestructionRateConstant(primer, primingSpecies, Y, T9, rho);
+
+            double equilibriumMassFraction = 0.0;
+
+            if (destructionRateConstant > 1e-99) {
+                const double creationRate = calculateCreationRate(primer, primingSpecies, Y, T9, rho);
+                equilibriumMassFraction = (creationRate / destructionRateConstant) * primingSpecies.mass();
+                LOG_INFO(logger, "Found equilibrium for {}: X_eq = {:.4e}", primingSpecies.name(), equilibriumMassFraction);
+
+                const reaction::Reaction* dominantChannel = findDominantCreationChannel(primer, primingSpecies, Y, T9, rho);
+
+                if (dominantChannel) {
+                    LOG_DEBUG(logger, "Dominant creation channel for {}: {}", primingSpecies.name(), dominantChannel->peName());
+
+                    double totalReactantMass = 0.0;
+                    for (const auto& reactant : dominantChannel->reactants()) {
+                        totalReactantMass += reactant.mass();
+                    }
+
+                    double scalingFactor = 1.0;
+                    for (const auto& reactant : dominantChannel->reactants()) {
+                        const double massToSubtract = equilibriumMassFraction * (reactant.mass() / totalReactantMass);
+                        double availableMass = 0.0;
+                        if (currentMassFractions.contains(reactant)) {
+                            availableMass = currentMassFractions.at(reactant);
+                        }
+                        if (massToSubtract > availableMass && availableMass > 0) {
+                            scalingFactor = std::min(scalingFactor, availableMass / massToSubtract);
+                        }
+                    }
+
+                    if (scalingFactor < 1.0) {
+                        LOG_WARNING(logger, "Priming for {} was limited by reactant availability. Scaling transfer by {:.4e}", primingSpecies.name(), scalingFactor);
+                        equilibriumMassFraction *= scalingFactor;
+                    }
+
+                    // Update the internal mass fraction map and accumulate total changes
+                    totalMassFractionChanges[primingSpecies] += equilibriumMassFraction;
+                    currentMassFractions[primingSpecies] += equilibriumMassFraction;
+
+                    for (const auto& reactant : dominantChannel->reactants()) {
+                        const double massToSubtract = equilibriumMassFraction * (reactant.mass() / totalReactantMass);
+                        totalMassFractionChanges[reactant] -= massToSubtract;
+                        currentMassFractions[reactant] -= massToSubtract;
+                    }
+                } else {
+                     LOG_ERROR(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name());
+                     report.status = PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL;
+                     totalMassFractionChanges[primingSpecies] += 1e-40;
+                     currentMassFractions[primingSpecies] += 1e-40;
+                }
+            } else {
+                LOG_WARNING(logger, "No destruction channel found for {}. Using fallback abundance.", primingSpecies.name());
+                totalMassFractionChanges[primingSpecies] += 1e-40;
+                currentMassFractions[primingSpecies] += 1e-40;
+                report.status = PrimingReportStatus::BASE_NETWORK_TOO_SHALLOW;
+            }
+        }
+
+        // --- Final Composition Construction ---
+        std::vector<std::string> final_symbols;
+        std::vector<double> final_mass_fractions;
+        for(const auto& [species, mass_fraction] : currentMassFractions) {
+            final_symbols.push_back(std::string(species.name()));
+            if (mass_fraction < 0.0 && std::abs(mass_fraction) < 1e-16) {
+                final_mass_fractions.push_back(0.0); // Avoid negative mass fractions
+            } else {
+                final_mass_fractions.push_back(mass_fraction);
+            }
+        }
+
+        // Create the final composition object from the pre-normalized mass fractions
+        Composition primedComposition(final_symbols, final_mass_fractions, true);
+
+        report.primedComposition = primedComposition;
+        for (const auto& [species, change] : totalMassFractionChanges) {
+            report.massFractionChanges.emplace_back(species, change);
+        }
+
+        engine.rebuild(netIn.composition, initialDepth);
+        return report;
     }
 
   double calculateDestructionRateConstant(

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

@@ -229,6 +229,10 @@ namespace gridfire {
         return Y; // Return the vector of molar abundances
     }
 
+    PrimingReport AdaptiveEngineView::primeEngine(const NetIn &netIn) {
+        return m_baseEngine.primeEngine(netIn);
+    }
+
     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()) {

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

@@ -76,8 +76,6 @@ namespace gridfire {
         return m_baseEngine;
     }
 
-
-
     const std::vector<Species> & DefinedEngineView::getNetworkSpecies() const {
         return m_activeSpecies;
     }
@@ -217,6 +215,10 @@ namespace gridfire {
         return Y; // Return the vector of molar abundances
     }
 
+    PrimingReport DefinedEngineView::primeEngine(const NetIn &netIn) {
+        return m_baseEngine.primeEngine(netIn);
+    }
+
     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;

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

@@ -6,6 +6,7 @@
 #include <queue>
 
 namespace gridfire {
+    static int s_operator_parens_called = 0;
     using fourdst::atomic::Species;
 
     MultiscalePartitioningEngineView::MultiscalePartitioningEngineView(
@@ -58,8 +59,8 @@ namespace gridfire {
     double MultiscalePartitioningEngineView::calculateMolarReactionFlow(
         const reaction::Reaction &reaction,
         const std::vector<double> &Y,
-        double T9,
+        const double T9,
-        double rho
+        const double rho
     ) const {
         return m_baseEngine.calculateMolarReactionFlow(reaction, Y, T9, rho);
     }
@@ -230,18 +231,45 @@ namespace gridfire {
         return m_dynamic_species;
     }
 
-    void MultiscalePartitioningEngineView::equilibrateNetwork(
+    PrimingReport MultiscalePartitioningEngineView::primeEngine(const NetIn &netIn) {
+        return m_baseEngine.primeEngine(netIn);
+    }
+
+    fourdst::composition::Composition 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);
+        const std::vector<double> Y_equilibrated = solveQSEAbundances(Y, T9, rho);
+        fourdst::composition::Composition composition;
+
+        std::vector<std::string> symbols;
+        symbols.reserve(m_baseEngine.getNetworkSpecies().size());
+        for (const auto& species : m_baseEngine.getNetworkSpecies()) {
+            symbols.emplace_back(species.name());
+        }
+        composition.registerSymbol(symbols);
+
+        std::vector<double> X;
+        X.reserve(Y_equilibrated.size());
+        for (size_t i = 0; i < Y_equilibrated.size(); ++i) {
+            const double molarMass = m_baseEngine.getNetworkSpecies()[i].mass();
+            X.push_back(Y_equilibrated[i] * molarMass); // Convert from molar abundance to mass fraction
         }
 
-    void MultiscalePartitioningEngineView::equilibrateNetwork(
+        for (size_t i = 0; i < Y_equilibrated.size(); ++i) {
-        const NetIn& netIn,
+            const auto& species = m_baseEngine.getNetworkSpecies()[i];
+            composition.setMassFraction(std::string(species.name()), X[i]);
+        }
+
+        composition.finalize(true);
+        return composition;
+    }
+
+    fourdst::composition::Composition MultiscalePartitioningEngineView::equilibrateNetwork(
+        const NetIn &netIn,
         const double dt_control
     ) {
         std::vector<double> Y(m_baseEngine.getNetworkSpecies().size(), 0.0);
@@ -257,10 +285,7 @@ namespace gridfire {
         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);
+        return equilibrateNetwork(Y, T9, rho, dt_control);
-        std::vector<double> Y_equilibrated = solveQSEAbundances(Y, T9, rho);
-
-
     }
 
     int MultiscalePartitioningEngineView::getSpeciesIndex(const fourdst::atomic::Species &species) const {
@@ -566,14 +591,13 @@ namespace gridfire {
             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 {
+        s_operator_parens_called++;
         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
 
@@ -586,6 +610,7 @@ namespace gridfire {
         for (size_t i = 0; i < m_qse_solve_indices.size(); ++i) {
             f_qse(i) = dydt[m_qse_solve_indices[i]];
         }
+
         return 0; // Success
     }
 

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

@@ -44,7 +44,8 @@ namespace gridfire {
         ),
         baseEngine
     ),
-    m_primingSpecies(primingSpecies) {}
+    m_primingSpecies(primingSpecies) {
+    }
 
 
     std::vector<std::string> NetworkPrimingEngineView::constructPrimingReactionSet(

src/network/lib/network.cpp

@@ -61,32 +61,6 @@ namespace gridfire {
         return oldFormat;
     }
 
-    reaction::LogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition, bool reverse) {
-        using namespace reaction;
-        std::vector<Reaction> reaclibReactions;
-        auto logger = fourdst::logging::LogManager::getInstance().getLogger("log");
-
-        for (const auto &reaction: reaclib::get_all_reactions()) {
-            if (reaction.is_reverse() != reverse) {
-                continue; // Skip reactions that do not match the requested direction
-            }
-            bool gotReaction = true;
-            const auto& reactants = reaction.reactants();
-            for (const auto& reactant : reactants) {
-                if (!composition.contains(reactant)) {
-                    gotReaction = false;
-                    break; // If any reactant is not in the composition, skip this reaction
-                }
-            }
-            if (gotReaction) {
-                LOG_TRACE_L3(logger, "Adding reaction {} to REACLIB reaction set.", reaction.peName());
-                reaclibReactions.push_back(reaction);
-            }
-        }
-        const ReactionSet reactionSet(reaclibReactions);
-        return packReactionSetToLogicalReactionSet(reactionSet);
-    }
-
     // Trim whitespace from both ends of a string
     std::string trim_whitespace(const std::string& str) {
         auto startIt = str.begin();

src/network/lib/solver/solver.cpp

@@ -4,7 +4,6 @@
 
 
 #include "gridfire/utils/logging.h"
-#include "gridfire/utils/qse_rules.h"
 
 #include "fourdst/composition/atomicSpecies.h"
 #include "fourdst/composition/composition.h"
@@ -176,11 +175,11 @@ namespace gridfire::solver {
 
             const double final_timescale = std::min(network_timescale, decay_timescale);
 
-            const bool isQSE = is_species_in_qse(
+            const bool isQSE = false;
-                network_timescale,
+            //     network_timescale,
-                decay_timescale,
+            //     decay_timescale,
-                abundance
+            //     abundance
-            );
+            // );
 
             if (isQSE) {
                 LOG_TRACE_L2(m_logger, "{} is in QSE based on rules in qse_rules.h", species.name());

src/network/meson.build

@@ -8,6 +8,7 @@ network_sources = files(
     'lib/engine/views/engine_multiscale.cpp',
     'lib/engine/views/engine_priming.cpp',
     'lib/engine/procedures/priming.cpp',
+    'lib/engine/procedures/construction.cpp',
     'lib/reaction/reaction.cpp',
     'lib/reaction/reaclib.cpp',
     'lib/io/network_file.cpp',
@@ -59,6 +60,7 @@ network_headers = files(
     'include/gridfire/engine/views/engine_multiscale.h',
     'include/gridfire/engine/views/engine_priming.h',
     'include/gridfire/engine/procedures/priming.h',
+    'include/gridfire/engine/procedures/construction.h',
     'include/gridfire/reaction/reaction.h',
     'include/gridfire/reaction/reaclib.h',
     'include/gridfire/io/network_file.h',
@@ -72,6 +74,5 @@ 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')