feat(jacobian): Added regularization

There are times when the jacobian matrix has infinities or nans. If
these cases correspond to species (rows or columns) which have
effectivley zero abundance (i.e. if Y(Cl-32) ~ 1e-310 and
(dY(H-2)/dt)/dY(Cl-32) is inf) then it is safe to regularize these
entries to 0. If this is not done then the solver will end up finding
NaN values for the molar abundances on subsequent steps. This has been
implimented through a small regularization function in the
CVODE_solver_strategy file.

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

@@ -9,6 +9,8 @@
 #include <unordered_map>
 
 namespace gridfire {
+    using JacobianEntry = std::pair<std::pair<fourdst::atomic::Species, fourdst::atomic::Species>, double>;
+
     class NetworkJacobian {
     public:
         explicit NetworkJacobian(
@@ -16,6 +18,10 @@ namespace gridfire {
             const std::function<fourdst::atomic::Species(size_t)> &indexToSpeciesFunc
         );
 
+        NetworkJacobian(const NetworkJacobian& jacobian);
+        NetworkJacobian(NetworkJacobian&& jacobian) noexcept;
+        NetworkJacobian& operator=(NetworkJacobian&& jacobian) noexcept;
+
         double operator()(
             const fourdst::atomic::Species& row,
             const fourdst::atomic::Species& col
@@ -26,12 +32,34 @@ namespace gridfire {
             size_t j
         ) const;
 
+        void set(
+            const fourdst::atomic::Species& row,
+            const fourdst::atomic::Species& col,
+            double value
+        );
+
+        void set(
+            size_t i,
+            size_t j,
+            double value
+        );
+
+        void set(
+            const JacobianEntry &entry
+        );
+
         std::tuple<size_t, size_t> shape() const;
         size_t rank() const;
         size_t nnz() const;
 
         bool singular() const;
 
+        [[nodiscard]] std::vector<JacobianEntry> infs() const;
+        [[nodiscard]] std::vector<JacobianEntry> nans() const;
+
+        [[nodiscard]] Eigen::SparseMatrix<double> data() const;
+        [[nodiscard]] const std::unordered_map<fourdst::atomic::Species, size_t>& mapping() const;
+
     private:
         Eigen::SparseMatrix<double> m_jacobianMatrix;
         std::unordered_map<fourdst::atomic::Species, size_t> m_speciesToIndexMap;

src/include/gridfire/exceptions/error_solver.h

@@ -22,4 +22,8 @@ namespace gridfire::exceptions {
     class SingularJacobianError final : public SolverError {
         using SolverError::SolverError;
     };
+
+    class IllConditionedJacobianError final : public SolverError {
+        using SolverError::SolverError;
+    };
 }

src/lib/engine/engine_graph.cpp

@@ -80,9 +80,11 @@ namespace gridfire {
         const double rho,
         const reaction::ReactionSet &activeReactions
     ) const {
+        LOG_TRACE_L2(m_logger, "Calculating RHS and Energy in GraphEngine at T9 = {}, rho = {}.", T9, rho);
         const double Ye = comp.getElectronAbundance();
-        const double mue = 5.0e-3 * std::pow(rho * Ye, 1.0 / 3.0) + 0.5 * T9;
+        const double mue = 0.0; // TODO: Remove
         if (m_usePrecomputation) {
+            LOG_TRACE_L2(m_logger, "Using precomputation for reaction rates in GraphEngine calculateRHSAndEnergy.");
             std::vector<double> bare_rates;
             std::vector<double> bare_reverse_rates;
             bare_rates.reserve(activeReactions.size());
@@ -96,9 +98,12 @@ namespace gridfire {
                 }
             }
 
+            LOG_TRACE_L2(m_logger, "Precomputed {} forward and {} reverse reaction rates for active reactions.", bare_rates.size(), bare_reverse_rates.size());
+
             // --- The public facing interface can always use the precomputed version since taping is done internally ---
             return calculateAllDerivativesUsingPrecomputation(comp, bare_rates, bare_reverse_rates, T9, rho, activeReactions);
         } else {
+            LOG_TRACE_L2(m_logger, "Not using precomputation for reaction rates in GraphEngine calculateRHSAndEnergy.");
             StepDerivatives<double> result = calculateAllDerivatives<double>(
                 comp.getMolarAbundanceVector(),
                 T9,
@@ -598,6 +603,7 @@ namespace gridfire {
         const double rho,
         const reaction::ReactionSet &activeReactions
     ) const {
+        LOG_TRACE_L2(m_logger, "Computing screening factors for {} active reactions.", activeReactions.size());
         // --- Calculate screening factors ---
         const std::vector<double> screeningFactors = m_screeningModel->calculateScreeningFactors(
             activeReactions,
@@ -672,6 +678,8 @@ namespace gridfire {
             reactionCounter++;
         }
 
+        LOG_TRACE_L2(m_logger, "Computed {} molar reaction flows for active reactions. Assembling these into RHS", molarReactionFlows.size());
+
         // --- Assemble molar abundance derivatives ---
         StepDerivatives<double> result;
         for (const auto& species: m_networkSpecies) {

src/lib/engine/types/jacobian.cpp

@@ -21,6 +21,31 @@ namespace gridfire {
         }
     }
 
+    NetworkJacobian::NetworkJacobian(
+        const NetworkJacobian &jacobian
+    ) : m_jacobianMatrix(jacobian.m_jacobianMatrix),
+        m_speciesToIndexMap(jacobian.m_speciesToIndexMap),
+        m_rank(jacobian.m_rank)
+    {}
+
+    NetworkJacobian::NetworkJacobian(
+        NetworkJacobian &&jacobian
+    ) noexcept : m_jacobianMatrix(std::move(jacobian.m_jacobianMatrix)),
+                 m_speciesToIndexMap(std::move(jacobian.m_speciesToIndexMap)),
+                 m_rank(jacobian.m_rank)
+    {}
+
+    NetworkJacobian & NetworkJacobian::operator=(
+        NetworkJacobian &&jacobian
+    ) noexcept {
+        if (this != &jacobian) {
+            m_jacobianMatrix = std::move(jacobian.m_jacobianMatrix);
+            m_speciesToIndexMap = std::move(jacobian.m_speciesToIndexMap);
+            m_rank = jacobian.m_rank;
+        }
+        return *this;
+    }
+
     double NetworkJacobian::operator()(const fourdst::atomic::Species &row, const fourdst::atomic::Species &col) const {
         if (!m_speciesToIndexMap.contains(row) || !m_speciesToIndexMap.contains(col)) {
             throw std::out_of_range("Species not found in NetworkJacobian operator().");
@@ -37,6 +62,26 @@ namespace gridfire {
         return m_jacobianMatrix.coeff(i, j);
     }
 
+    void NetworkJacobian::set(const fourdst::atomic::Species &row, const fourdst::atomic::Species &col, const double value) {
+        if (!m_speciesToIndexMap.contains(row) || !m_speciesToIndexMap.contains(col)) {
+            throw std::out_of_range("Species not found in NetworkJacobian set().");
+        }
+        const size_t i = m_speciesToIndexMap.at(row);
+        const size_t j = m_speciesToIndexMap.at(col);
+        set(i, j, value);
+    }
+
+    void NetworkJacobian::set(const size_t i, const size_t j, const double value) {
+        if (i >= m_jacobianMatrix.rows() || j >= m_jacobianMatrix.cols()) {
+            throw std::out_of_range("Index out of bounds in NetworkJacobian set().");
+        }
+        m_jacobianMatrix.coeffRef(i, j) = value;
+    }
+
+    void NetworkJacobian::set(const JacobianEntry &entry) {
+        set(entry.first.first, entry.first.second, entry.second);
+    }
+
     std::tuple<size_t, size_t> NetworkJacobian::shape() const {
         return {m_jacobianMatrix.rows(), m_jacobianMatrix.cols()};
     }
@@ -56,5 +101,60 @@ namespace gridfire {
         return m_rank < minDim;
     }
 
+    std::vector<JacobianEntry> NetworkJacobian::infs() const {
+        std::vector<JacobianEntry> infs;
+        for (int k=0; k<m_jacobianMatrix.outerSize(); ++k) {
+            for (Eigen::SparseMatrix<double>::InnerIterator it(m_jacobianMatrix,k); it; ++it) {
+                if (std::isinf(it.value())) {
+                    fourdst::atomic::Species rowSpecies = std::ranges::find_if(
+                        m_speciesToIndexMap,
+                        [it](const auto& pair) {
+                            return pair.second == static_cast<size_t>(it.row());
+                        })->first;
+
+                    fourdst::atomic::Species colSpecies = std::ranges::find_if(
+                        m_speciesToIndexMap,
+                        [it](const auto& pair) {
+                            return pair.second == static_cast<size_t>(it.col());
+                        })->first;
+
+                    infs.emplace_back(std::make_pair(rowSpecies, colSpecies), it.value());
+                }
+            }
+        }
+        return infs;
+    }
+
+    std::vector<JacobianEntry> NetworkJacobian::nans() const {
+        std::vector<JacobianEntry> nans;
+        for (int k=0; k<m_jacobianMatrix.outerSize(); ++k) {
+            for (Eigen::SparseMatrix<double>::InnerIterator it(m_jacobianMatrix,k); it; ++it) {
+                if (std::isnan(it.value())) {
+                    fourdst::atomic::Species rowSpecies = std::ranges::find_if(
+                        m_speciesToIndexMap,
+                        [it](const auto& pair) {
+                            return pair.second == static_cast<size_t>(it.row());
+                        })->first;
+
+                    fourdst::atomic::Species colSpecies = std::ranges::find_if(
+                        m_speciesToIndexMap,
+                        [it](const auto& pair) {
+                            return pair.second == static_cast<size_t>(it.col());
+                        })->first;
+
+                    nans.emplace_back(std::make_pair(rowSpecies, colSpecies), it.value());
+                }
+            }
+        }
+        return nans;
 
     }
+
+    Eigen::SparseMatrix<double> NetworkJacobian::data() const {
+        return m_jacobianMatrix;
+    }
+
+    const std::unordered_map<fourdst::atomic::Species, size_t> & NetworkJacobian::mapping() const {
+        return m_speciesToIndexMap;
+    }
+}

src/lib/engine/views/engine_adaptive.cpp

@@ -82,11 +82,47 @@ namespace gridfire {
         const double T9,
         const double rho
     ) const {
+        LOG_TRACE_L2(m_logger, "Calculating RHS and Energy in AdaptiveEngineView at T9 = {}, rho = {}.", T9, rho);
         validateState();
+        LOG_TRACE_L2(
+            m_logger,
+            "Adaptive engine view state validated prior to composition collection. Input Composition: {}",
+            [&comp]() -> std::string {
+                std::stringstream ss;
+                size_t i = 0;
+                for (const auto& [species, abundance] : comp) {
+                    ss << species.name() << ": " << abundance;
+                    if (i < comp.size() - 1) {
+                        ss << ", ";
+                    }
+                    i++;
+                }
+                return ss.str();
+            }());
         fourdst::composition::Composition collectedComp = collectComposition(comp, T9, rho);
+        LOG_TRACE_L2(
+            m_logger,
+            "Composition Collected prior to passing to base engine. Collected Composition: {}",
+            [&comp, &collectedComp]() -> std::string {
+                std::stringstream ss;
+                size_t i = 0;
+                for (const auto& [species, abundance] : collectedComp) {
+                    ss << species.name() << ": " << abundance;
+                    if (comp.contains(species)) {
+                        ss << " (input: " << comp.getMolarAbundance(species) << ")";
+                    }
+                    if (i < collectedComp.size() - 1) {
+                        ss << ", ";
+                    }
+                    i++;
+                }
+                return ss.str();
+            }());
         auto result = m_baseEngine.calculateRHSAndEnergy(collectedComp, T9, rho);
+        LOG_TRACE_L2(m_logger, "Base engine calculation of RHS and Energy complete.");
 
         if (!result) {
+            LOG_TRACE_L2(m_logger, "Base engine returned stale error during RHS and Energy calculation.");
             return std::unexpected{result.error()};
         }
 

src/lib/engine/views/engine_multiscale.cpp

@@ -163,19 +163,51 @@ namespace gridfire {
         const double T9,
         const double rho
     ) const {
-        // const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho);
+        LOG_TRACE_L2(m_logger, "Calculating RHS and Energy in MultiscalePartitioningEngineView at T9 = {}, rho = {}.", T9, rho);
+        LOG_TRACE_L2(m_logger, "Input composition is {}", [&comp]() -> std::string {
+            std::stringstream ss;
+            size_t i = 0;
+            for (const auto& [species, abundance] : comp) {
+                ss << species.name() << ": " << abundance;
+                if (i < comp.size() - 1) {
+                    ss << ", ";
+                }
+                i++;
+            }
+            return ss.str();
+        }());
+        const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho);
+        LOG_TRACE_L2(m_logger, "Equilibrated composition prior to calling base engine is {}", [&qseComposition, &comp]() -> std::string {
+            std::stringstream ss;
+            size_t i = 0;
+            for (const auto& [species, abundance] : qseComposition) {
+                ss << species.name() << ": " << abundance;
+                if (comp.contains(species)) {
+                    ss << " (input: " << comp.getMolarAbundance(species) << ")";
+                }
+                if (i < qseComposition.size() - 1) {
+                    ss << ", ";
+                }
+                i++;
+            }
+            return ss.str();
+        }());
 
         const auto result = m_baseEngine.calculateRHSAndEnergy(comp, T9, rho);
+        LOG_TRACE_L2(m_logger, "Base engine calculation of RHS and Energy complete.");
 
         if (!result) {
+            LOG_TRACE_L2(m_logger, "Base engine returned stale error during RHS and Energy calculation.");
             return std::unexpected{result.error()};
         }
 
         auto deriv = result.value();
 
+        LOG_TRACE_L2(m_logger, "Zeroing out algebraic species derivatives.");
         for (const auto& species : m_algebraic_species) {
             deriv.dydt[species] = 0.0; // Fix the algebraic species to the equilibrium abundances we calculate.
         }
+        LOG_TRACE_L2(m_logger, "Done Zeroing out algebraic species derivatives.");
         return deriv;
     }
 
@@ -1157,13 +1189,35 @@ namespace gridfire {
         const double rho
     ) const {
         LOG_TRACE_L1(m_logger, "Solving for QSE abundances...");
-
+        LOG_TRACE_L2(m_logger, "Composition before QSE solving: {}", [&comp]() -> std::string {
-        fourdst::composition::Composition outputComposition(comp.getRegisteredSpecies());
+            std::stringstream ss;
+            size_t i = 0;
             for (const auto& [sp, y] : comp) {
-            outputComposition.setMolarAbundance(sp, y);
+                ss << std::format("{}: {}", sp.name(), y);
+                if (i < comp.size() - 1) {
+                    ss << ", ";
                 }
+                i++;
+            }
+            return ss.str();
+        }());
+
+        fourdst::composition::Composition outputComposition(comp);
 
         for (const auto&[is_in_equilibrium, algebraic_species, seed_species, mean_timescale] : m_qse_groups) {
+            LOG_TRACE_L2(m_logger, "Working on QSE group with algebraic species: {}",
+                [&]() -> std::string {
+                    std::stringstream ss;
+                    size_t count = 0;
+                    for (const auto& species: algebraic_species) {
+                        ss << species.name();
+                        if (count < algebraic_species.size() - 1) {
+                            ss << ", ";
+                        }
+                        count++;
+                    }
+                    return ss.str();
+                }());
             if (!is_in_equilibrium || (algebraic_species.empty() && seed_species.empty())) {
                 continue;
             }
@@ -1178,16 +1232,19 @@ namespace gridfire {
                 const double Y = std::max(initial_abundance, abundance_floor);
                 v_initial(i) = std::log(Y);
                 species_to_index_map.emplace(species, i);
+                LOG_TRACE_L2(m_logger, "For species {} initial molar abundance is {}, log scaled to {}. Species placed at index {}.", species.name(), Y, v_initial(i), i);
                 i++;
             }
 
+            LOG_TRACE_L2(m_logger, "Setting up Eigen Levenberg-Marquardt solver for QSE group...");
             EigenFunctor functor(*this, algebraic_species, comp, T9, rho, Y_scale, species_to_index_map);
             Eigen::LevenbergMarquardt lm(functor);
             lm.parameters.ftol = 1.0e-10;
             lm.parameters.xtol = 1.0e-10;
 
-            LOG_TRACE_L1(m_logger, "Minimizing functor...");
+            LOG_TRACE_L2(m_logger, "Minimizing functor...");
             Eigen::LevenbergMarquardtSpace::Status status = lm.minimize(v_initial);
+            LOG_TRACE_L2(m_logger, "Minimizing functor status: {}", lm_status_map.at(status));
 
             if (status <= 0 || status >= 4) {
                 std::stringstream msg;
@@ -1459,13 +1516,14 @@ namespace gridfire {
         }
         LOG_TRACE_L2(
             m_view.m_logger,
-            "Functor evaluation at T9 = {}, rho = {}, y_qse = <{}> complete. ||f|| = {}",
+            "Functor evaluation at T9 = {}, rho = {}, y_qse (v_qse) = <{}> complete. ||f|| = {}",
             m_T9,
             m_rho,
             [&]() -> std::string {
                 std::stringstream ss;
                 for (long j = 0; j < y_qse.size(); ++j) {
                     ss << y_qse(j);
+                    ss << "(" << v_qse(j) << ")";
                     if (j < y_qse.size() - 1) {
                         ss << ", ";
                     }

src/lib/solver/strategies/CVODE_solver_strategy.cpp

@@ -24,6 +24,8 @@
 #include "gridfire/exceptions/error_solver.h"
 
 namespace {
+    constexpr double MIN_ABUNDANCE_TO_CONTRIBUTE_TO_JACOBIAN = 1e-100;
+
     std::unordered_map<int, std::string> cvode_ret_code_map {
         {0, "CV_SUCCESS: The solver succeeded."},
         {1, "CV_TSTOP_RETURN: The solver reached the specified stopping time."},
@@ -78,6 +80,34 @@ namespace {
         check_cvode_flag(vec == nullptr ? -1 : 0, "N_VNew");
         return vec;
     }
+
+    gridfire::NetworkJacobian regularize_jacobian(const gridfire::NetworkJacobian& jacobian, const fourdst::composition::CompositionAbstract& comp, std::optional<quill::Logger*> logger = std::nullopt) {
+        const std::vector<gridfire::JacobianEntry> infs = jacobian.infs();
+        const std::vector<gridfire::JacobianEntry> nans = jacobian.nans();
+        if (infs.size() == 0 && nans.size() == 0) {
+            return jacobian;
+        }
+
+        gridfire::NetworkJacobian newJacobian = jacobian;
+        for (const auto& [iSp, dSp] : infs | std::views::keys) {
+            if (comp.getMolarAbundance(iSp) < MIN_ABUNDANCE_TO_CONTRIBUTE_TO_JACOBIAN || comp.getMolarAbundance(dSp) < MIN_ABUNDANCE_TO_CONTRIBUTE_TO_JACOBIAN) {
+                newJacobian.set(iSp, dSp, 0.0);
+                if (logger) {
+                    LOG_TRACE_L1(logger.value(), "Regularized Jacobian entry ({}, {}) from inf to 0.0 due to low abundance.", iSp.name(), dSp.name());
+                }
+            }
+        }
+        for (const auto& [iSp, dSp] : nans | std::views::keys) {
+            if (comp.getMolarAbundance(iSp) < MIN_ABUNDANCE_TO_CONTRIBUTE_TO_JACOBIAN || comp.getMolarAbundance(dSp) < MIN_ABUNDANCE_TO_CONTRIBUTE_TO_JACOBIAN) {
+                newJacobian.set(iSp, dSp, 0.0);
+                if (logger) {
+                    LOG_TRACE_L1(logger.value(), "Regularized Jacobian entry ({}, {}) from inf to 0.0 due to low abundance.", iSp.name(), dSp.name());
+                }
+            }
+        }
+
+        return newJacobian;
+    }
 }
 
 namespace gridfire::solver {
@@ -206,7 +236,9 @@ namespace gridfire::solver {
 
             check_cvode_flag(CVodeSetUserData(m_cvode_mem, &user_data), "CVodeSetUserData");
 
+            LOG_TRACE_L2(m_logger, "Taking one CVODE step...");
             int flag = CVode(m_cvode_mem, netIn.tMax, m_Y, &current_time, CV_ONE_STEP);
+            LOG_TRACE_L2(m_logger, "CVODE step complete. Current time: {}, step status: {}", current_time, cvode_ret_code_map.at(flag));
 
             if (user_data.captured_exception){
                 std::rethrow_exception(std::make_exception_ptr(*user_data.captured_exception));
@@ -539,13 +571,17 @@ namespace gridfire::solver {
         const auto* instance = data->solver_instance;
 
         try {
+            LOG_TRACE_L2(instance->m_logger, "CVODE RHS wrapper called at time {}", t);
             const CVODERHSOutputData out = instance->calculate_rhs(t, y, ydot, data);
             data->reaction_contribution_map = out.reaction_contribution_map;
+            LOG_TRACE_L2(instance->m_logger, "CVODE RHS wrapper completed successfully at time {}", t);
             return 0;
         } catch (const exceptions::StaleEngineTrigger& e) {
+            LOG_ERROR(instance->m_logger, "StaleEngineTrigger caught in CVODE RHS wrapper at time {}: {}", t, e.what());
             data->captured_exception = std::make_unique<exceptions::StaleEngineTrigger>(e);
             return 1; // 1 Indicates a recoverable error, CVODE will retry the step
         } catch (...) {
+            LOG_CRITICAL(instance->m_logger, "Unrecoverable and Unknown exception caught in CVODE RHS wrapper at time {}", t);
             return -1; // Some unrecoverable error
         }
     }
@@ -562,10 +598,13 @@ namespace gridfire::solver {
     ) {
         const auto* data = static_cast<CVODEUserData*>(user_data);
         const auto* engine = data->engine;
+        const auto* solver_instance = data->solver_instance;
 
+        LOG_TRACE_L2(solver_instance->m_logger, "CVODE Jacobian wrapper starting");
         const size_t numSpecies = engine->getNetworkSpecies().size();
         sunrealtype* y_data = N_VGetArrayPointer(y);
 
+
         // Solver constraints should keep these values very close to 0 but floating point noise can still result in very
         // small negative numbers which can result in NaN's and more immediate crashes in the composition
         // finalization stage
@@ -576,12 +615,75 @@ namespace gridfire::solver {
         }
         std::vector<double> y_vec(y_data, y_data + numSpecies);
         fourdst::composition::Composition composition(engine->getNetworkSpecies(), y_vec);
+        LOG_TRACE_L2(solver_instance->m_logger, "Generating Jacobian matrix at time {} with {} species in composition (mean molecular mass: {})", t, composition.size(), composition.getMeanParticleMass());
+        LOG_TRACE_L2(solver_instance->m_logger, "Composition is {}", [&composition]() -> std::string {
+            std::stringstream ss;
+            size_t i = 0;
+            for (const auto& [species, abundance] : composition) {
+                ss << species.name() << ": " << abundance;
+                if (i < composition.size() - 1) {
+                    ss << ", ";
+                }
+                i++;
+            }
+            return ss.str();
+        }());
 
+        LOG_TRACE_L2(solver_instance->m_logger, "Generating Jacobian matrix at time {}", t);
         NetworkJacobian jac = engine->generateJacobianMatrix(composition, data->T9, data->rho);
+        LOG_TRACE_L2(solver_instance->m_logger, "Regularizing Jacobian matrix at time {}", t);
+        jac = regularize_jacobian(jac, composition, solver_instance->m_logger);
+        LOG_TRACE_L2(solver_instance->m_logger, "Done regularizing Jacobian matrix at time {}", t);
+        if (jac.infs().size() != 0 || jac.nans().size() != 0) {
+            auto infString = [&jac]() -> std::string {
+                std::stringstream ss;
+                size_t i = 0;
+                std::vector<JacobianEntry> entries = jac.infs();
+                for (const auto &[fst, snd]: entries | std::views::keys) {
+                    ss << "J(" << fst.name() << ", " << snd.name() << ")";
+                    if (i < entries.size() - 1) {
+                        ss << ", ";
+                    }
+                    i++;
+                }
+                if (entries.size() == 0) {
+                    ss << "None";
+                }
+                return ss.str();
+            };
+            auto nanString = [&jac]() -> std::string {
+                std::stringstream ss;
+                size_t i = 0;
+                std::vector<JacobianEntry> entries = jac.nans();
+                for (const auto &[fst, snd]: entries | std::views::keys) {
+                    ss << "J(" << fst.name() << ", " << snd.name() << ")";
+                    if (i < entries.size() - 1) {
+                        ss << ", ";
+                    }
+                    i++;
+                }
+                if (entries.size() == 0) {
+                    ss << "None";
+                }
+                return ss.str();
+            };
+            LOG_ERROR(
+                solver_instance->m_logger,
+                "Jacobian matrix generated at time {} contains {} infinite entries ({}) and {} NaN entries ({}). This will lead to a solver failure. GridFire will now halt.",
+                t,
+                jac.infs().size(),
+                infString(),
+                jac.nans().size(),
+                nanString()
+            );
+            throw exceptions::IllConditionedJacobianError(std::format("Jacobian matrix generated at time {} contains {} infinite entries ({}) and {} NaN entries ({}). This will lead to a solver failure. In order to ensure tractability GridFire will not proceed. Focus on improving conditioning of the Jacobian matrix. If you believe this is an error please contact the GridFire developers.", t, jac.infs().size(), infString(), jac.nans().size(), nanString()));
+        }
+        LOG_TRACE_L2(solver_instance->m_logger, "Jacobian matrix created at time {} of shape ({} x {}) and rank {}", t, std::get<0>(jac.shape()), std::get<1>(jac.shape()), jac.rank());
 
         sunrealtype* J_data = SUNDenseMatrix_Data(J);
         const long int N = SUNDenseMatrix_Columns(J);
 
+        LOG_TRACE_L2(solver_instance->m_logger, "Transferring Jacobian matrix data to SUNDenseMatrix format at time {}", t);
         for (size_t j = 0; j < numSpecies; ++j) {
             const fourdst::atomic::Species& species_j = engine->getNetworkSpecies()[j];
             for (size_t i = 0; i < numSpecies; ++i) {
@@ -589,16 +691,10 @@ namespace gridfire::solver {
                 // J(i,j) = d(f_i)/d(y_j)
                 // Column-major order format for SUNDenseMatrix: J_data[j*N + i] indexes J(i,j)
                 const double dYi_dt = jac(species_i, species_j);
-                // if (i == j && dYi_dt == 0 && engine->getSpeciesStatus(species_i) == SpeciesStatus::ACTIVE) {
-                //     std::cerr << "Warning: Jacobian matrix has a zero on the diagonal for species " << species_i.name() << ". This may lead to solver failure or pathological stiffness.\n";
-                //     // throw exceptions::SingularJacobianError(
-                //     //     "Jacobian matrix has a zero on the diagonal for species " + std::string(species_i.name()) +
-                //     //     ". This will either lead to solver failure or pathological stiffness. In order to ensure tractability GridFire will not proceed. Focus on improving conditioning of the Jacobian matrix. If you believe this is an error please contact the GridFire developers."
-                //     // );
-                // }
                 J_data[j * N + i] = dYi_dt;
             }
         }
+        LOG_TRACE_L2(solver_instance->m_logger, "Done transferring Jacobian matrix data to SUNDenseMatrix format at time {}", t);
 
         // For now assume that the energy derivatives wrt. abundances are zero
         // TODO: Need a better way to build this part of the output jacobian so it properly pushes the solver
@@ -631,13 +727,28 @@ namespace gridfire::solver {
         std::vector<double> y_vec(y_data, y_data + numSpecies);
         fourdst::composition::Composition composition(m_engine.getNetworkSpecies(), y_vec);
 
+        LOG_TRACE_L2(m_logger, "Calculating RHS at time {} with {} species in composition (mean molecular mass: {})", t, composition.size(), composition.getMeanParticleMass());
         const auto result = m_engine.calculateRHSAndEnergy(composition, data->T9, data->rho);
         if (!result) {
+            LOG_WARNING(m_logger, "StaleEngineTrigger thrown during RHS calculation at time {}", t);
             throw exceptions::StaleEngineTrigger({data->T9, data->rho, y_vec, t, m_num_steps, y_data[numSpecies]});
         }
 
         sunrealtype* ydot_data = N_VGetArrayPointer(ydot);
         const auto& [dydt, nuclearEnergyGenerationRate, reactionContributions] = result.value();
+        LOG_TRACE_L2(m_logger, "Done calculating RHS at time {}, specific nuclear energy generation rate: {}", t, nuclearEnergyGenerationRate);
+        LOG_TRACE_L2(m_logger, "RHS at time {} is {}", t, [&dydt]() -> std::string {
+            std::stringstream ss;
+            size_t i = 0;
+            for (const auto& [species, rate] : dydt) {
+                ss << "dY(" << species.name() << ")/dt" << ": " << rate;
+                if (i < dydt.size() - 1) {
+                    ss << ", ";
+                }
+                i++;
+            }
+            return ss.str();
+        }());
 
         for (size_t i = 0; i < numSpecies; ++i) {
             fourdst::atomic::Species species = m_engine.getNetworkSpecies()[i];
@@ -657,6 +768,7 @@ namespace gridfire::solver {
         const double relTol,
         const double accumulatedEnergy
     ) {
+        LOG_TRACE_L2(m_logger, "Initializing CVODE integration resources with N: {}, current_time: {}, absTol: {}, relTol: {}", N, current_time, absTol, relTol);
         cleanup_cvode_resources(false); // Cleanup any existing resources before initializing new ones
 
         m_Y = init_sun_vector(N, m_sun_ctx);
@@ -706,9 +818,11 @@ namespace gridfire::solver {
 
         check_cvode_flag(CVodeSetLinearSolver(m_cvode_mem, m_LS, m_J), "CVodeSetLinearSolver");
         check_cvode_flag(CVodeSetJacFn(m_cvode_mem, cvode_jac_wrapper), "CVodeSetJacFn");
+        LOG_TRACE_L2(m_logger, "CVODE solver initialized");
     }
 
     void CVODESolverStrategy::cleanup_cvode_resources(const bool memFree) {
+        LOG_TRACE_L2(m_logger, "Cleaning up cvode resources");
         if (m_LS) SUNLinSolFree(m_LS);
         if (m_J) SUNMatDestroy(m_J);
         if (m_Y) N_VDestroy(m_Y);
@@ -725,6 +839,7 @@ namespace gridfire::solver {
             if (m_cvode_mem) CVodeFree(&m_cvode_mem);
             m_cvode_mem = nullptr;
         }
+        LOG_TRACE_L2(m_logger, "Done Cleaning up cvode resources");
     }
 
     void CVODESolverStrategy::log_step_diagnostics(const CVODEUserData &user_data, bool displayJacobianStiffness) const {