feat(priming): Priming now uses solver to initialize species abundances

Instead of a complex system of identifying dominate channels and
approximating abundances based on that, priming now simply ignites a
basic network at 1e7K for 1e-15s. This is an effective approach to prime
relevant species while being short enough to not change the abundance of
any fuel species.

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

@@ -27,13 +27,8 @@ namespace gridfire {
      * @see PrimingReport for data associated with each status.
      */
     enum class PrimingReportStatus {
-        FULL_SUCCESS = 0,
+        SUCCESS,
-        NO_SPECIES_TO_PRIME = 1,
+        SOLVER_FAILURE,
-        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
     };
 
     /**
@@ -43,13 +38,8 @@ namespace gridfire {
      * The map contains entries for all PrimingReportStatus values.
      */
     inline std::map<PrimingReportStatus, std::string> PrimingReportStatusStrings = {
-        {PrimingReportStatus::FULL_SUCCESS, "Full Success"},
+        {PrimingReportStatus::SUCCESS, "SUCCESS"},
-        {PrimingReportStatus::NO_SPECIES_TO_PRIME, "No Species to Prime"},
+        {PrimingReportStatus::SOLVER_FAILURE, "SOLVER_FAILURE"},
-        {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"}
     };
 
     /**
@@ -65,11 +55,6 @@ namespace gridfire {
     struct PrimingReport {
         /** Finalized composition after priming. */
         fourdst::composition::Composition primedComposition;
-        /**
-         * List of pairs (species, rate change) representing destruction (<0)
-         * or creation (>0) rates of species during priming.
-         */
-        std::vector<std::pair<fourdst::atomic::Species, double>> massFractionChanges;
         /** True if priming completed without error. */
         bool success;
         /** Detailed status code indicating the result. */
@@ -93,4 +78,11 @@ namespace gridfire {
         }
     };
 
+    enum class SpeciesStatus {
+        ACTIVE,
+        EQUILIBRIUM,
+        INACTIVE_FLOW,
+        NOT_PRESENT
+    };
+
 }

src/lib/engine/engine_graph.cpp

@@ -25,6 +25,7 @@
 #include <ranges>
 
 #include <boost/numeric/odeint.hpp>
+#include <boost/numeric/ublas/matrix_sparse.hpp>
 
 #include "cppad/cppad.hpp"
 #include "cppad/utility/sparse_rc.hpp"
@@ -180,7 +181,6 @@ namespace gridfire {
         populateSpeciesToIndexMap();
         collectAtomicReverseRateAtomicBases();
         generateStoichiometryMatrix();
-        reserveJacobianMatrix();
 
         recordADTape(); // Record the AD tape for the RHS of the ODE (dY/di and dEps/di) for all independent variables i
 
@@ -264,18 +264,6 @@ namespace gridfire {
         }
     }
 
-    void GraphEngine::reserveJacobianMatrix() const {
-        // The implementation of this function (and others) constrains this nuclear network to a constant temperature and density during
-        // each evaluation.
-        const size_t numSpecies = m_networkSpecies.size();
-        m_jacobianMatrix.clear();
-        m_jacobianMatrix.resize(numSpecies, numSpecies, false); // Sparse matrix, no initial values
-        LOG_TRACE_L2(m_logger, "Jacobian matrix resized to {} rows and {} columns.",
-                 m_jacobianMatrix.size1(), m_jacobianMatrix.size2());
-
-        m_jacobianMatrixState = JacobianMatrixState::UNINITIALIZED;
-    }
-
     // --- Basic Accessors and Queries ---
     const std::vector<fourdst::atomic::Species>& GraphEngine::getNetworkSpecies() const {
         return m_networkSpecies;
@@ -287,7 +275,6 @@ namespace gridfire {
 
     void GraphEngine::setNetworkReactions(const reaction::ReactionSet &reactions) {
         m_reactions = reactions;
-        m_jacobianMatrixState = JacobianMatrixState::STALE;
         syncInternalMaps();
     }
 
@@ -521,9 +508,6 @@ namespace gridfire {
     }
 
     void GraphEngine::setUseReverseReactions(const bool useReverse) {
-        if (useReverse != m_useReverseReactions) {
-            m_jacobianMatrixState = JacobianMatrixState::STALE;
-        }
         m_useReverseReactions = useReverse;
     }
 
@@ -572,7 +556,6 @@ namespace gridfire {
         if (depth != m_depth) {
             m_depth = depth;
             m_reactions = build_nuclear_network(comp, m_weakRateInterpolator, m_depth);
-            m_jacobianMatrixState = JacobianMatrixState::STALE;
             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.");
@@ -599,6 +582,14 @@ namespace gridfire {
         return result;
     }
 
+    SpeciesStatus GraphEngine::getSpeciesStatus(const fourdst::atomic::Species &species) const {
+        if (m_networkSpeciesMap.contains(species.name())) {
+            return SpeciesStatus::ACTIVE;
+        }
+        return SpeciesStatus::NOT_PRESENT;
+
+    }
+
     StepDerivatives<double> GraphEngine::calculateAllDerivativesUsingPrecomputation(
         const fourdst::composition::CompositionAbstract &comp,
         const std::vector<double> &bare_rates,
@@ -784,7 +775,6 @@ namespace gridfire {
     void GraphEngine::setScreeningModel(const screening::ScreeningType model) {
         m_screeningModel = screening::selectScreeningModel(model);
         m_screeningType = model;
-        m_jacobianMatrixState = JacobianMatrixState::STALE; // The screening model affects the jacobian so if its changed the jacobian must be made stale
     }
 
     screening::ScreeningType GraphEngine::getScreeningModel() const {
@@ -828,7 +818,7 @@ namespace gridfire {
         );
     }
 
-    void GraphEngine::generateJacobianMatrix(
+    NetworkJacobian GraphEngine::generateJacobianMatrix(
         const fourdst::composition::CompositionAbstract &comp,
         const double T9,
         const double rho
@@ -847,7 +837,7 @@ namespace gridfire {
         std::vector<double> adInput(numSpecies + 2, 0.0); // +2 for T9 and rho
         const std::vector<double>& Y_dynamic = mutableComp.getMolarAbundanceVector();
         for (size_t i = 0; i < numSpecies; ++i) {
-            adInput[i] = std::max(Y_dynamic[i], 1e-99); // regularize the jacobian...
+            adInput[i] = Y_dynamic[i];
         }
         adInput[numSpecies]     = T9;  // T9
         adInput[numSpecies + 1] = rho; // rho
@@ -856,26 +846,34 @@ namespace gridfire {
         const std::vector<double> dotY = m_rhsADFun.Jacobian(adInput);
 
         // 3. Pack jacobian vector into sparse matrix
-        m_jacobianMatrix.clear();
+        Eigen::SparseMatrix<double> jacobianMatrix(numSpecies, numSpecies);
-        // std::vector<std::unique_ptr<utils::ColumnBase>> columns;
+        std::vector<Eigen::Triplet<double> > triplets;
         for (size_t i = 0; i < numSpecies; ++i) {
-            // std::vector<double> colData;
             for (size_t j = 0; j < numSpecies; ++j) {
-                const double value = dotY[i * (numSpecies + 2) + j];
+                double value = dotY[i * (numSpecies + 2) + j];
                 if (std::abs(value) > MIN_JACOBIAN_THRESHOLD || i == j) { // Always keep diagonal elements to avoid pathological stiffness
-                    m_jacobianMatrix(i, j) = value;
+                    if (i == j && value == 0) {
+                        LOG_WARNING(m_logger, "While generating the Jacobian matrix, a zero diagonal element was encountered at index ({}, {}) (species: {}, abundance: {}). This may lead to numerical instability. Setting to -1 to avoid singularity", i, j, m_networkSpecies[i].name(), adInput[i]);
+                        // value = -1.0;
+                    }
+                    triplets.emplace_back(i, j, value);
                 }
-                // colData.push_back(value);
             }
-            // columns.push_back(std::make_unique<utils::Column<double>>(std::to_string(i), colData));
         }
-        // std::cout << utils::format_table("Jacobian after dense calculation", columns) << std::endl;
+        LOG_TRACE_L1_LIMIT_EVERY_N(1000, m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", jacobianMatrix.rows(), jacobianMatrix.cols());
-        // exit(0);
+
-        LOG_TRACE_L1_LIMIT_EVERY_N(1000, m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2());
+        auto index_to_species = [this](const size_t index) -> fourdst::atomic::Species {
-        m_jacobianMatrixState = JacobianMatrixState::READY_DENSE;
+            if (index < m_networkSpecies.size()) {
+                return m_networkSpecies[index];
+            }
+            throw std::out_of_range("Index out of range in index_to_species mapping.");
+        };
+        jacobianMatrix.setFromTriplets(triplets.begin(), triplets.end());
+        NetworkJacobian jac(jacobianMatrix, index_to_species);
+        return jac;
     }
 
-    void GraphEngine::generateJacobianMatrix(
+    NetworkJacobian GraphEngine::generateJacobianMatrix(
         const fourdst::composition::CompositionAbstract &comp,
         const double T9,
         const double rho,
@@ -904,10 +902,10 @@ namespace gridfire {
         }
 
         // --- 3. Call the sparse reverse-mode implementation ---
-        generateJacobianMatrix(comp, T9, rho, sparsityPattern);
+        return generateJacobianMatrix(comp, T9, rho, sparsityPattern);
     }
 
-    void GraphEngine::generateJacobianMatrix(
+    NetworkJacobian GraphEngine::generateJacobianMatrix(
         const fourdst::composition::CompositionAbstract &comp,
         const double T9,
         const double rho,
@@ -929,7 +927,7 @@ namespace gridfire {
         // }
         size_t i = 0;
         for (const auto& species: m_networkSpecies) {
-            double Yi = 0.0;
+            double Yi = 0.0; // Small floor to avoid issues with zero abundances
             if (comp.contains(species)) {
                 Yi = comp.getMolarAbundance(species);
             }
@@ -974,46 +972,26 @@ namespace gridfire {
             m_jac_work // Work vector for CppAD
         );
 
-        // --- Convert the sparse Jacobian back to the Boost uBLAS format ---
+        Eigen::SparseMatrix<double> jacobianMatrix(numSpecies, numSpecies);
-        m_jacobianMatrix.clear();
+        std::vector<Eigen::Triplet<double> > triplets;
         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 || row == col) { // Always keep diagonal elements to avoid pathological stiffness
-                m_jacobianMatrix(row, col) = value; // Insert into the sparse matrix
+                triplets.emplace_back(row, col, value);
             }
         }
-        m_jacobianMatrixState = JacobianMatrixState::READY_SPARSE;
+        jacobianMatrix.setFromTriplets(triplets.begin(), triplets.end());
-    }
+        auto index_to_species = [this](const size_t index) -> fourdst::atomic::Species {
-
+            if (index < m_networkSpecies.size()) {
-    double GraphEngine::getJacobianMatrixEntry(
+                return m_networkSpecies[index];
-        const fourdst::atomic::Species& rowSpecies,
-        const fourdst::atomic::Species& colSpecies
-    ) const {
-        switch (m_jacobianMatrixState) {
-            case JacobianMatrixState::STALE: {
-                const std::string staleMsg = std::format("Cannot retrieve jacobian entry for row {}, column {} as jacobian matrix is stale (has not been regenerated since last network topology change)", rowSpecies.name(), colSpecies.name());
-                throw exceptions::StaleJacobianError(staleMsg);
             }
-            case JacobianMatrixState::UNINITIALIZED: {
+            throw std::out_of_range("Index out of range in index_to_species mapping.");
-                const std::string unInitMsg = std::format("Cannot retrieve jacobian entry for row {}, column {} as jacobian matrix is uninitialized (will return all 0s)", rowSpecies.name(), colSpecies.name());
+        };
-                throw exceptions::UninitializedJacobianError(unInitMsg);
+        NetworkJacobian jac(jacobianMatrix, index_to_species);
-            }
+        return jac;
-            case JacobianMatrixState::READY_DENSE:
-                [[fallthrough]];
-            case JacobianMatrixState::READY_SPARSE: {
-                const size_t i = getSpeciesIndex(rowSpecies);
-                const size_t j = getSpeciesIndex(colSpecies);
-                return m_jacobianMatrix(i, j);
-            }
-            default: {
-                // Code should not be able to get into this state
-                const std::string msg = std::format("An unknown error has occurred while attempting to retrieve the jacobian element at row {}, column {}. This should be taken as a catastrophic failure and reported to GridFire developers.", rowSpecies.name(), colSpecies.name());
-                throw exceptions::UnknownJacobianError(msg);
-            }
-        }
     }
 
     std::unordered_map<fourdst::atomic::Species, int> GraphEngine::getNetReactionStoichiometry(
@@ -1287,7 +1265,6 @@ namespace gridfire {
             }
         );
 
-
         // Extract the raw vector from the associative map
         std::vector<CppAD::AD<double>> dependentVector;
         dependentVector.reserve(dydt.size() + 1);
@@ -1302,8 +1279,7 @@ namespace gridfire {
 
         m_rhsADFun.Dependent(adInput, dependentVector);
 
-        LOG_TRACE_L1(m_logger, "AD tape recorded successfully for the RHS and Eps calculation. Number of independent variables: {}.",
+        LOG_TRACE_L1(m_logger, "AD tape recorded successfully for the RHS and Eps calculation. Number of independent variables: {}.", adInput.size());
-                 adInput.size());
     }
 
     void GraphEngine::collectAtomicReverseRateAtomicBases() {

src/lib/engine/procedures/priming.cpp

@@ -3,317 +3,79 @@
 #include "fourdst/atomic/species.h"
 #include "fourdst/composition/utils.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 "gridfire/exceptions/error_solver.h"
 
 #include "fourdst/logging/logging.h"
+#include "gridfire/solver/strategies/CVODE_solver_strategy.h"
 #include "quill/Logger.h"
 #include "quill/LogMacros.h"
 
-namespace {
-    bool isReactionIgnorable(
-        const gridfire::reaction::Reaction& reaction,
-        const std::optional<std::vector<gridfire::reaction::ReactionType>>& reactionsTypesToIgnore
-    ) {
-        if (reactionsTypesToIgnore.has_value()) {
-            for (const auto& type : reactionsTypesToIgnore.value()) {
-                if (reaction.type() == type) {
-                    return true;
-                }
-            }
-        }
-        return false;
-    }
-}
 
 namespace gridfire {
     using fourdst::composition::Composition;
     using fourdst::atomic::Species;
 
-    const reaction::Reaction* findDominantCreationChannel (
-        const DynamicEngine& engine,
-        const Species& species,
-        const Composition &comp,
-        const double T9,
-        const double rho,
-        const std::optional<std::vector<reaction::ReactionType>> &reactionsTypesToIgnore
-    ) {
-        const reaction::Reaction* dominateReaction = nullptr;
-        double maxFlow = -1.0;
-        for (const auto& reaction : engine.getNetworkReactions()) {
-            if (isReactionIgnorable(*reaction, reactionsTypesToIgnore)) continue;
-
-            if (reaction->contains(species) && reaction->stoichiometry(species) > 0) {
-                const double flow = engine.calculateMolarReactionFlow(*reaction, comp, T9, rho);
-                if (flow > maxFlow) {
-                    maxFlow = flow;
-                    dominateReaction = reaction.get();
-                }
-            }
-        }
-        return dominateReaction;
-    }
-
-
-
-    /**
-     * @brief Primes absent species in the network to their equilibrium abundances using a robust, two-stage approach.
-     *
-     * @details This function implements a robust network priming algorithm that avoids the pitfalls of
-     * sequential, one-by-one priming. The previous, brittle method could allow an early priming
-     * reaction to consume all of a shared reactant, starving later reactions. This new, two-stage
-     * method ensures that all priming reactions are considered collectively, competing for the
-     * same limited pool of initial reactants in a physically consistent manner.
-     *
-     * The algorithm proceeds in three main stages:
-     * 1.  **Calculation Stage:** It first loops through all species that need priming. For each one,
-     * it calculates its theoretical equilibrium mass fraction and identifies the dominant
-     * creation channel. Crucially, it *does not* modify any abundances at this stage. Instead,
-     * it stores these calculations as a list of "mass transfer requests".
-     *
-     * 2.  **Collective Scaling Stage:** It then processes the full list of requests to determine the
-     * total "debit" required from each reactant. By comparing these total debits to the
-     * initially available mass of each reactant, it calculates a single, global `scalingFactor`.
-     * If any reactant is overdrawn, this factor will be less than 1.0, ensuring that no
-     * reactant's abundance can go negative.
-     *
-     * 3.  **Application Stage:** Finally, it loops through the requests again, applying the mass
-     * transfers. Each calculated equilibrium mass fraction and corresponding reactant debit is
-     * multiplied by the global `scalingFactor` before being applied to the final composition.
-     * This ensures that if resources are limited, all primed species are scaled down proportionally.
-     *
-     * @param netIn Input network data containing initial composition, temperature, and density.
-     * @param engine DynamicEngine used to build and evaluate the reaction network.
-     * @param ignoredReactionTypes Types of reactions to ignore during priming (e.g., weak reactions).
-     * @return PrimingReport encapsulating the results of the priming operation, including the new
-     * robustly primed composition.
-     */
     PrimingReport primeNetwork(
-        const NetIn& netIn,
+    const NetIn& netIn,
-        GraphEngine& engine,
+    GraphEngine& engine,
-        const std::optional<std::vector<reaction::ReactionType>>& ignoredReactionTypes
+    const std::optional<std::vector<reaction::ReactionType>>& ignoredReactionTypes
     ) {
-        auto logger = LogManager::getInstance().getLogger("log");
+        const auto logger = LogManager::getInstance().getLogger("log");
+        solver::CVODESolverStrategy integrator(engine);
+        integrator.set_stdout_logging_enabled(false);
+        NetIn solverInput(netIn);
 
-        // --- Initial Setup ---
+        solverInput.tMax = 1e-15;
-        // Identify all species with zero initial abundance that need to be primed.
+        solverInput.temperature = 1e7;
-        std::vector<Species> speciesToPrime;
-        for (const auto &[sp, y]: netIn.composition) {
-            if (y == 0.0) {
-                speciesToPrime.push_back(sp);
-            }
-        }
 
-        // Sort priming species by mass number, lightest to heaviest.
+        LOG_INFO(logger, "Short timescale ({}) network ignition started.", solverInput.tMax);
-        std::ranges::sort(speciesToPrime, [](const Species& a, const Species& b) {
-            return a.mass() < b.mass();
-        });
-
-        LOG_DEBUG(logger, "Priming {} species in the network.", speciesToPrime.size());
-
-        // If no species need priming, return immediately.
         PrimingReport report;
-        if (speciesToPrime.empty()) {
+        try {
-            report.primedComposition = netIn.composition;
+            const NetOut netOut = integrator.evaluate(solverInput, false);
-            report.success = true;
+            LOG_INFO(logger, "Network ignition completed.");
-            report.status = PrimingReportStatus::NO_SPECIES_TO_PRIME;
+            LOG_TRACE_L2(
-            return report;
+                logger,
-        }
+                "After ignition composition is {}",
-
+                [netOut, netIn]() -> std::string {
-        const double T9 = netIn.temperature / 1e9;
+                    std::stringstream ss;
-        const double rho = netIn.density;
+                    size_t i = 0;
-        const auto initialReactionSet = engine.getNetworkReactions();
+                    for (const auto& [species, abundance] : netOut.composition) {
-
+                        ss << species.name() << ": " << abundance << " (prior: " << netIn.composition.getMolarAbundance(species);
-        report.status = PrimingReportStatus::FULL_SUCCESS;
+                        ss << ", fractional change: " << (abundance - netIn.composition.getMolarAbundance(species)) / netIn.composition.getMolarAbundance(species) * 100.0 << "%)";
-        report.success = true;
+                        if (i < netOut.composition.size() - 1) {
-
+                            ss << ", ";
-        // Build full set of species
+                        }
-        std::set<Species> allSpecies;
+                        ++i;
-        for (const auto &sp: netIn.composition | std::views::keys) {
+                    }
-            allSpecies.insert(sp);
+                    return ss.str();
-        }
+                }()
-        for (const auto& sp : speciesToPrime) {
-            allSpecies.insert(sp);
-        }
-
-        // Rebuild the engine with the full network to ensure all possible creation channels are available.
-        engine.rebuild(netIn.composition, NetworkBuildDepth::Full);
-
-        // Initialize mutable molar abundances for all species
-        std::map<Species, double> molarAbundances;
-        for (const auto& sp : allSpecies) {
-            molarAbundances[sp] = netIn.composition.contains(sp) ? netIn.composition.getMolarAbundance(sp) : 0.0;
-        }
-
-        // --- Prime Each Species ---
-        // Since molar abundances are independent, we can directly calculate and apply changes
-        std::unordered_map<Species, double> totalMolarAbundanceChanges;
-
-        for (const auto& primingSpecies : speciesToPrime) {
-            // Create a temporary composition reflecting the current state for rate calculations.
-            Composition tempComp(allSpecies);
-            for (const auto& [sp, abundance] : molarAbundances) {
-                tempComp.setMolarAbundance(sp, abundance);
-            }
-
-            NetworkPrimingEngineView primer(primingSpecies, engine);
-            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;
-            }
-
-            const double destructionRateConstant = calculateDestructionRateConstant(
-                primer,
-                primingSpecies,
-                tempComp,
-                T9,
-                rho,
-                ignoredReactionTypes
             );
-
+            report.primedComposition = netOut.composition;
-            double equilibriumMolarAbundance = 0.0;
+            std::unordered_set<Species> unprimedSpecies;
-            std::vector<Species> reactants;
+            double minAbundance = std::numeric_limits<double>::max();
-
+            for (const auto& [sp, y] : report.primedComposition) {
-            if (destructionRateConstant > 1e-99) {
+                if (y == 0) {
-                const double creationRate = calculateCreationRate(
+                    unprimedSpecies.insert(sp);
-                    primer,
-                    primingSpecies,
-                    tempComp,
-                    T9,
-                    rho,
-                    ignoredReactionTypes
-                );
-
-                // Equilibrium: creation rate = destruction rate constant * molar abundance
-                equilibriumMolarAbundance = creationRate / destructionRateConstant;
-
-                if (std::isnan(equilibriumMolarAbundance)) {
-                    equilibriumMolarAbundance = 0.0;
                 }
-
+                if (y != 0 && y < minAbundance) {
-                if (const reaction::Reaction* dominantChannel = findDominantCreationChannel(
+                    minAbundance = y;
-                    primer,
-                    primingSpecies,
-                    tempComp,
-                    T9,
-                    rho,
-                    ignoredReactionTypes)
-                ) {
-                    reactants = dominantChannel->reactants();
-                } else {
-                    LOG_TRACE_L1(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name());
-                    report.status = PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL;
-                    reactants.clear(); // Use fallback
-                }
-            } else {
-                LOG_TRACE_L2(logger, "No destruction channel found for {}. Using fallback abundance.", primingSpecies.name());
-                // For species with no destruction, use a tiny fallback abundance
-                equilibriumMolarAbundance = 1e-40;
-            }
-
-            // Add the equilibrium molar abundance to the primed species
-            molarAbundances.at(primingSpecies) += equilibriumMolarAbundance;
-            totalMolarAbundanceChanges[primingSpecies] += equilibriumMolarAbundance;
-
-            // Subtract from reactants proportionally to their stoichiometry
-            if (!reactants.empty()) {
-                const double totalStoichiometry = static_cast<double>(reactants.size());
-                const double abundancePerReactant = equilibriumMolarAbundance / totalStoichiometry;
-
-                for (const auto& reactant : reactants) {
-                    LOG_TRACE_L1(logger, "Transferring {:.4e} molar abundance from {} to {} during priming.",
-                        abundancePerReactant, reactant.name(), primingSpecies.name());
-
-                    if (!molarAbundances.contains(reactant)) {
-                        continue;
-                    }
-                    molarAbundances.at(reactant) -= abundancePerReactant;
-                    totalMolarAbundanceChanges[reactant] -= abundancePerReactant;
-
-                    // Ensure non-negative abundances
-                    if (molarAbundances.at(reactant) < 0.0) {
-                        LOG_WARNING(logger, "Species {} went negative during priming. Clamping to zero.", reactant.name());
-                        totalMolarAbundanceChanges[reactant] += molarAbundances.at(reactant); // Adjust change to reflect clamp
-                        molarAbundances.at(reactant) = 0.0;
-                    }
                 }
             }
+            double abundanceForUnprimedSpecies = minAbundance / 1e10;
+            for (const auto& sp : unprimedSpecies) {
+                LOG_TRACE_L1(logger, "Clamping Species {}: initial abundance {}, primed abundance {} to {}", sp.name(), netIn.composition.getMolarAbundance(sp), report.primedComposition.getMolarAbundance(sp), abundanceForUnprimedSpecies);
+                report.primedComposition.setMolarAbundance(sp, abundanceForUnprimedSpecies);
+            }
+            report.success = true;
+            report.status = PrimingReportStatus::SUCCESS;
+        } catch (const exceptions::SolverError& e) {
+            LOG_ERROR(logger, "Failed to prime network: solver failure encountered: {}", e.what());
+            std::rethrow_exception(std::current_exception());
         }
-
-        // --- Final Composition Construction ---
-        std::vector<Species> final_species;
-        std::vector<double> final_molar_abundances;
-
-        for (const auto& [species, abundance] : molarAbundances) {
-            final_species.emplace_back(species);
-            final_molar_abundances.push_back(std::max(0.0, abundance));
-
-            LOG_TRACE_L1(logger, "After priming, species {} has molar abundance {:.4e} (had {:0.4e} prior to priming).",
-                species.name(),
-                std::max(0.0, abundance),
-                netIn.composition.getMolarAbundance(species));
-        }
-
-        Composition primedComposition(final_species, final_molar_abundances);
-
-        report.primedComposition = primedComposition;
-
-        // Convert molar abundance changes to mass fraction changes for reporting
-        for (const auto& [species, molarChange] : totalMolarAbundanceChanges) {
-            double massFractionChange = molarChange * species.mass();
-            report.massFractionChanges.emplace_back(species, massFractionChange);
-        }
-
-        // Restore the engine to its original, smaller network state.
-        engine.setNetworkReactions(initialReactionSet);
-
         return report;
     }
-
-  double calculateDestructionRateConstant(
-      const DynamicEngine& engine,
-      const Species& species,
-      const Composition& composition,
-      const double T9,
-      const double rho,
-      const std::optional<std::vector<reaction::ReactionType>> &reactionTypesToIgnore
-  ) {
-        Composition unrestrictedComp(composition);
-        unrestrictedComp.setMolarAbundance(species, 1.0);
-        double destructionRateConstant = 0.0;
-        for (const auto& reaction: engine.getNetworkReactions()) {
-            if (isReactionIgnorable(*reaction, reactionTypesToIgnore)) continue;
-
-            const int stoichiometry = reaction->stoichiometry(species);
-            if (stoichiometry < 0) {
-                destructionRateConstant += std::abs(stoichiometry) * engine.calculateMolarReactionFlow(*reaction, unrestrictedComp, T9, rho);
-            }
-        }
-        return destructionRateConstant;
-    }
-
-    double calculateCreationRate(
-        const DynamicEngine& engine,
-        const Species& species,
-        const Composition& composition,
-        const double T9,
-        const double rho,
-        const std::optional<std::vector<reaction::ReactionType>> &reactionTypesToIgnore
-    ) {
-        double creationRate = 0.0;
-        for (const auto& reaction: engine.getNetworkReactions()) {
-            if (isReactionIgnorable(*reaction, reactionTypesToIgnore)) continue;
-
-            const int stoichiometry = reaction->stoichiometry(species);
-            if (stoichiometry > 0) {
-                creationRate += stoichiometry * engine.calculateMolarReactionFlow(*reaction, composition, T9, rho);
-            }
-        }
-        return creationRate;
-    }
 }