diff --git a/src/network/include/gridfire/engine/engine_abstract.h b/src/network/include/gridfire/engine/engine_abstract.h
index b7aa6a9d..584062dd 100644
--- a/src/network/include/gridfire/engine/engine_abstract.h
+++ b/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.");
+        }
     };
 }
\ No newline at end of file
diff --git a/src/network/include/gridfire/engine/engine_graph.h b/src/network/include/gridfire/engine/engine_graph.h
index 100c49c9..ac162d42 100644
--- a/src/network/include/gridfire/engine/engine_graph.h
+++ b/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,11 +110,16 @@ 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,
-                             const partition::PartitionFunction& partitionFunction
-                             );
+        explicit GraphEngine(
+            const fourdst::composition::Composition &composition,
+            const partition::PartitionFunction& partitionFunction,
+            const BuildDepthType buildDepth = NetworkBuildDepth::Full
+        );
 
         /**
          * @brief Constructs a GraphEngine from a set of reactions.
@@ -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(),
diff --git a/src/network/include/gridfire/engine/procedures/construction.h b/src/network/include/gridfire/engine/procedures/construction.h
index 39859633..525eea53 100644
--- a/src/network/include/gridfire/engine/procedures/construction.h
+++ b/src/network/include/gridfire/engine/procedures/construction.h
@@ -1,8 +1,17 @@
-//
-// Created by Emily Boudreaux on 7/14/25.
-//
+#pragma once
 
-#ifndef CONSTRUCTION_H
-#define CONSTRUCTION_H
+#include "gridfire/reaction/reaction.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
+    );
+}
\ No newline at end of file
diff --git a/src/network/include/gridfire/engine/procedures/priming.h b/src/network/include/gridfire/engine/procedures/priming.h
index d3cc91e1..22f365a0 100644
--- a/src/network/include/gridfire/engine/procedures/priming.h
+++ b/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
     );
diff --git a/src/network/include/gridfire/engine/types/building.h b/src/network/include/gridfire/engine/types/building.h
new file mode 100644
index 00000000..1fd1a7c7
--- /dev/null
+++ b/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>;
+}
diff --git a/src/network/include/gridfire/engine/types/reporting.h b/src/network/include/gridfire/engine/types/reporting.h
index 238582df..52a55c3a 100644
--- a/src/network/include/gridfire/engine/types/reporting.h
+++ b/src/network/include/gridfire/engine/types/reporting.h
@@ -1,8 +1,42 @@
-//
-// Created by Emily Boudreaux on 7/14/25.
-//
+#pragma once
 
-#ifndef REPORTING_H
-#define REPORTING_H
+#include <map>
+#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();
+        }
+    };
+}
\ No newline at end of file
diff --git a/src/network/include/gridfire/engine/views/engine_adaptive.h b/src/network/include/gridfire/engine/views/engine_adaptive.h
index 510bc6c3..71f2e8b9 100644
--- a/src/network/include/gridfire/engine/views/engine_adaptive.h
+++ b/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;
diff --git a/src/network/include/gridfire/engine/views/engine_defined.h b/src/network/include/gridfire/engine/views/engine_defined.h
index 454199b5..76303196 100644
--- a/src/network/include/gridfire/engine/views/engine_defined.h
+++ b/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;
diff --git a/src/network/include/gridfire/engine/views/engine_multiscale.h b/src/network/include/gridfire/engine/views/engine_multiscale.h
index bf66980c..7eaeeb28 100644
--- a/src/network/include/gridfire/engine/views/engine_multiscale.h
+++ b/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;
-        const std::vector<fourdst::atomic::Species>& getDynamicSpecies() const;
+        [[nodiscard]] PrimingReport primeEngine(const NetIn &netIn) override;
 
-        void equilibrateNetwork(
-            const std::vector<double>& Y,
+        [[nodiscard]] std::vector<fourdst::atomic::Species> getFastSpecies() const;
+        [[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(
-            const NetIn& netIn,
+        fourdst::composition::Composition equilibrateNetwork(
+            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(); }
-            int inputs() const { return m_qse_solve_indices.size(); }
+            [[nodiscard]] int values() 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;
diff --git a/src/network/include/gridfire/engine/views/engine_priming.h b/src/network/include/gridfire/engine/views/engine_priming.h
index d38471df..93c88706 100644
--- a/src/network/include/gridfire/engine/views/engine_priming.h
+++ b/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(
diff --git a/src/network/include/gridfire/network.h b/src/network/include/gridfire/network.h
index 74c87e65..c5ee5f33 100644
--- a/src/network/include/gridfire/network.h
+++ b/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
diff --git a/src/network/lib/engine/engine_graph.cpp b/src/network/lib/engine/engine_graph.cpp
index 2de28e86..5e762c4b 100644
--- a/src/network/lib/engine/engine_graph.cpp
+++ b/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
-    ): GraphEngine(composition, partition::GroundStatePartitionFunction()) {}
+        const fourdst::composition::Composition &composition,
+        const BuildDepthType buildDepth
+    ): GraphEngine(composition, partition::GroundStatePartitionFunction(), buildDepth) {}
 
     GraphEngine::GraphEngine(
         const fourdst::composition::Composition &composition,
-        const partition::PartitionFunction& partitionFunction) :
-    m_reactions(build_reaclib_nuclear_network(composition, false)),
-    m_partitionFunction(partitionFunction.clone())  // Clone the partition function to ensure ownership
+        const partition::PartitionFunction& partitionFunction,
+        const BuildDepthType buildDepth) :
+    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;
-            bool below_threshold = false;
+            double forwardAbundanceProduct = 1.0;
+            // 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];
-                if (abundance < MIN_ABUNDANCE_THRESHOLD) {
-                    below_threshold = true;
-                    break;
-                }
+                // const double abundance = Y_in[reactantIndex];
+                // if (abundance < MIN_ABUNDANCE_THRESHOLD) {
+                //     below_threshold = true;
+                //     break;
+                // }
 
-                abundanceProduct *= 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
+                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
+            // }
 
             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 *
-                    std::pow(rho, numReactants >  1 ? numReactants - 1 : 0);
+                    forwardAbundanceProduct *
+                    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);
diff --git a/src/network/lib/engine/procedures/construction.cpp b/src/network/lib/engine/procedures/construction.cpp
index 93dfc79b..d891c67e 100644
--- a/src/network/lib/engine/procedures/construction.cpp
+++ b/src/network/lib/engine/procedures/construction.cpp
@@ -1,3 +1,113 @@
-//
-// Created by Emily Boudreaux on 7/14/25.
-//
+#include "gridfire/engine/procedures/construction.h"
+
+#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);
+    }
+
+
+
+}
\ No newline at end of file
diff --git a/src/network/lib/engine/procedures/priming.cpp b/src/network/lib/engine/procedures/priming.cpp
index b1850a30..c57b4f73 100644
--- a/src/network/lib/engine/procedures/priming.cpp
+++ b/src/network/lib/engine/procedures/priming.cpp
@@ -1,51 +1,192 @@
 #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;
+    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;
+    }
 
 
-        NetIn currentConditions = netIn;
+    PrimingReport primeNetwork(const NetIn& netIn, DynamicEngine& engine) {
+        auto logger = LogManager::getInstance().getLogger("log");
+
         std::vector<Species> speciesToPrime;
         for (const auto &entry: netIn.composition | std::views::values) {
             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 rho = netIn.density; // Density in g/cm^3
+        const double T9 = netIn.temperature / 1e9;
+        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);
-            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();
+            if (primer.getNetworkReactions().size() == 0) {
+                LOG_ERROR(logger, "No priming reactions found for species {}.", primingSpecies.name());
+                report.success = false;
+                report.status = PrimingReportStatus::FAILED_TO_FIND_PRIMING_REACTIONS;
+                continue;
+            }
 
-            currentConditions.composition.setMassFraction(
-                std::string(primingSpecies.name()),
-                equilibriumMassFraction
-            );
+            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;
+            }
         }
 
-        return currentConditions.composition;
+        // --- 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(
+  double calculateDestructionRateConstant(
         const DynamicEngine& engine,
         const fourdst::atomic::Species& species,
         const std::vector<double>& Y,
diff --git a/src/network/lib/engine/views/engine_adaptive.cpp b/src/network/lib/engine/views/engine_adaptive.cpp
index 9cd09ebf..d40f36cb 100644
--- a/src/network/lib/engine/views/engine_adaptive.cpp
+++ b/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()) {
diff --git a/src/network/lib/engine/views/engine_defined.cpp b/src/network/lib/engine/views/engine_defined.cpp
index d5a922de..f2818bac 100644
--- a/src/network/lib/engine/views/engine_defined.cpp
+++ b/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;
diff --git a/src/network/lib/engine/views/engine_multiscale.cpp b/src/network/lib/engine/views/engine_multiscale.cpp
index 4d8f5509..003dad8a 100644
--- a/src/network/lib/engine/views/engine_multiscale.cpp
+++ b/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,
-        double rho
+        const double T9,
+        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
+        }
+
+        for (size_t i = 0; i < Y_equilibrated.size(); ++i) {
+            const auto& species = m_baseEngine.getNetworkSpecies()[i];
+            composition.setMassFraction(std::string(species.name()), X[i]);
+        }
+
+        composition.finalize(true);
+        return composition;
     }
 
-    void MultiscalePartitioningEngineView::equilibrateNetwork(
-        const NetIn& netIn,
+    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);
-        std::vector<double> Y_equilibrated = solveQSEAbundances(Y, T9, rho);
-
-
+        return equilibrateNetwork(Y, T9, rho, dt_control);
     }
 
     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
     }
 
diff --git a/src/network/lib/engine/views/engine_priming.cpp b/src/network/lib/engine/views/engine_priming.cpp
index 1abc3a65..7132eff8 100644
--- a/src/network/lib/engine/views/engine_priming.cpp
+++ b/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(
diff --git a/src/network/lib/network.cpp b/src/network/lib/network.cpp
index 1c161069..850a3043 100644
--- a/src/network/lib/network.cpp
+++ b/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();
diff --git a/src/network/lib/solver/solver.cpp b/src/network/lib/solver/solver.cpp
index 65701d7c..52e1ca46 100644
--- a/src/network/lib/solver/solver.cpp
+++ b/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(
-                network_timescale,
-                decay_timescale,
-                abundance
-            );
+            const bool isQSE = false;
+            //     network_timescale,
+            //     decay_timescale,
+            //     abundance
+            // );
 
             if (isQSE) {
                 LOG_TRACE_L2(m_logger, "{} is in QSE based on rules in qse_rules.h", species.name());
diff --git a/src/network/meson.build b/src/network/meson.build
index f31e8a21..7016d081 100644
--- a/src/network/meson.build
+++ b/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')