feat(KINSOL): Switch from Eigen to KINSOL

Previously QSE solving was done using Eigen. While this worked we were
limited in the ability to use previous iterations to speed up later
steps. We have switched to KINSOL, from SUNDIALS, for linear solving.
This has drastically speed up the process of solving for QSE abundances,
primarily because the jacobian matrix does not need to be generated
every single time time a QSE abundance is requested.

.gitignore

@@ -79,6 +79,7 @@ subprojects/fourdst/
 subprojects/libplugin/
 subprojects/minizip-ng-4.0.10/
 subprojects/cvode-*/
+subprojects/kinsol-*/
 *.fbundle
 *.wraplock
 

build-config/meson.build

@@ -3,7 +3,7 @@ cmake = import('cmake')
 subdir('fourdst')
 subdir('libplugin')
 
-subdir('cvode')
+subdir('sundials')
 
 subdir('boost')
 subdir('cppad')

build-config/sundials/cvode/meson.build

@@ -34,7 +34,7 @@ sundials_sunmatrixdense_dep = cvode_sp.dependency('sundials_sunmatrixdense_share
 # For the dense linear solver library
 sundials_sunlinsoldense_dep = cvode_sp.dependency('sundials_sunlinsoldense_shared')
 
-sundials_dep = declare_dependency(
+cvode_dep = declare_dependency(
     dependencies: [
         sundials_core_dep,
         sundials_cvode_dep,

build-config/sundials/kinsol/meson.build

@@ -0,0 +1,29 @@
+cmake = import('cmake')
+
+kinsol_cmake_options = cmake.subproject_options()
+
+kinsol_cmake_options.add_cmake_defines({
+    'CMAKE_CXX_FLAGS'   : '-Wno-deprecated-declarations',
+    'CMAKE_C_FLAGS'     : '-Wno-deprecated-declarations',
+    'BUILD_SHARED_LIBS' : 'ON',
+    'BUILD_STATIC_LIBS' : 'OFF',
+})
+
+kinsol_cmake_options.add_cmake_defines({
+    'CMAKE_INSTALL_LIBDIR': get_option('libdir'),
+    'CMAKE_INSTALL_INCLUDEDIR': get_option('includedir')
+})
+
+kinsol_sp = cmake.subproject(
+    'kinsol',
+    options: kinsol_cmake_options,
+)
+
+sundials_kinsol_shared = kinsol_sp.dependency('sundials_kinsol_shared')
+
+kinsol_dep = declare_dependency(
+    dependencies: [
+        sundials_kinsol_shared,
+    ]
+)
+

build-config/sundials/meson.build

@@ -0,0 +1,9 @@
+subdir('cvode')
+subdir('kinsol')
+
+sundials_dep = declare_dependency(
+    dependencies: [
+        cvode_dep,
+        kinsol_dep,
+    ],
+)

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

@@ -3,6 +3,10 @@
 #include "gridfire/engine/engine_abstract.h"
 #include "gridfire/engine/views/engine_view_abstract.h"
 #include "gridfire/engine/engine_graph.h"
+#include "sundials/sundials_linearsolver.h"
+#include "sundials/sundials_matrix.h"
+#include "sundials/sundials_nvector.h"
+#include "sundials/sundials_types.h"
 
 #include "unsupported/Eigen/NonLinearOptimization"
 
@@ -78,6 +82,8 @@ namespace gridfire {
          */
         explicit MultiscalePartitioningEngineView(DynamicEngine& baseEngine);
 
+        ~MultiscalePartitioningEngineView() override;
+
         /**
          * @brief Gets the list of species in the network.
          * @return A const reference to the vector of `Species` objects representing all species
@@ -611,120 +617,81 @@ namespace gridfire {
             [[nodiscard]] [[maybe_unused]] std::string toString() const;
         };
 
-        /**
+        class QSESolver {
-         * @brief Functor for solving QSE abundances using Eigen's nonlinear optimization.
+        public:
-         *
+            QSESolver(
-         * @par Purpose
+                const std::vector<fourdst::atomic::Species>& species,
-         * This struct provides the objective function (`operator()`) and its Jacobian
+                const DynamicEngine& engine,
-         * (`df`) to Eigen's Levenberg-Marquardt solver. The goal is to find the abundances
+                SUNContext sun_ctx
-         * of algebraic species that make their time derivatives (`dY/dt`) equal to zero.
+            );
-         *
+
-         * @how
+            QSESolver(
-         * - **`operator()`**: Takes a vector `v_qse` (scaled abundances of algebraic species) as input.
+                const QSESolver& other
-         *   It constructs a full trial abundance vector `y_trial`, calls the base engine's
+            ) = delete;
-         *   `calculateRHSAndEnergy`, and returns the `dY/dt` values for the algebraic species.
+            QSESolver& operator=(
-         *   The solver attempts to drive this return vector to zero.
+                const QSESolver& other
-         * - **`df`**: Computes the Jacobian of the objective function. It calls the base engine's
+            ) = delete;
-         *   `generateJacobianMatrix` and extracts the sub-matrix corresponding to the algebraic
+
-         *   species. It applies the chain rule to account for the `asinh` scaling used on the
+            ~QSESolver();
-         *   abundances.
+
-         *
+            fourdst::composition::Composition solve(
-         * The abundances are scaled using `asinh` to handle the large dynamic range and ensure positivity.
+                const fourdst::composition::Composition& comp,
-         */
+                double T9,
-        struct EigenFunctor {
+                double rho
-            using InputType = Eigen::Matrix<double, Eigen::Dynamic, 1>;
+            ) const;
-            using OutputType = Eigen::Matrix<double, Eigen::Dynamic, 1>;
+
-            using JacobianType = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>;
+            size_t solves() const;
-            enum {
+
-                InputsAtCompileTime = Eigen::Dynamic,
+            void log_diagnostics() const;
-                ValuesAtCompileTime = Eigen::Dynamic
+        private:
+
+            static int sys_func(
+                N_Vector y,
+                N_Vector f,
+                void *user_data
+            );
+            static int sys_jac(
+                N_Vector y,
+                N_Vector fy,
+                SUNMatrix J,
+                void *user_data,
+                N_Vector tmp1,
+                N_Vector tmp2
+            );
+
+            struct UserData {
+                const DynamicEngine& engine;
+                double T9;
+                double rho;
+                fourdst::composition::Composition& comp;
+                const std::unordered_map<fourdst::atomic::Species, size_t>& qse_solve_species_index_map;
+                const std::vector<fourdst::atomic::Species>& qse_solve_species;
             };
 
-            /**
+        private:
-             * @brief Pointer to the MultiscalePartitioningEngineView instance.
+            mutable size_t m_solves = 0;
-             */
+            size_t m_N;
-            const MultiscalePartitioningEngineView& m_view;
+            const DynamicEngine& m_engine;
-            /**
+            std::vector<fourdst::atomic::Species> m_species;
-             * @brief The set of species to solve for in the QSE group.
+            std::unordered_map<fourdst::atomic::Species, size_t> m_speciesMap;
-             */
-            const std::set<fourdst::atomic::Species>& m_qse_solve_species;
-            /**
-             * @brief Initial abundances of all species in the full network.
-             */
-            const fourdst::composition::CompositionAbstract& m_initial_comp;
-            /**
-             * @brief Temperature in units of 10^9 K.
-             */
-            const double m_T9;
-            /**
-             * @brief Density in g/cm^3.
-             */
-            const double m_rho;
-            /**
-             * @brief Scaling factors for the species abundances, used to improve solver stability.
-             */
-            const Eigen::VectorXd& m_Y_scale;
 
-            /**
+            // --- SUNDIALS / KINSOL Persistent resources
-             * @brief Mapping from species to their indices in the QSE solve vector.
+            SUNContext m_sun_ctx   = nullptr;
-             */
+            void* m_kinsol_mem     = nullptr;
-            const std::unordered_map<fourdst::atomic::Species, size_t> m_qse_solve_species_index_map;
 
-            mutable std::optional<JacobianType> m_cached_jacobian = std::nullopt;
+            N_Vector m_Y           = nullptr;
+            N_Vector m_scale       = nullptr;
+            N_Vector m_f_scale     = nullptr;
+            N_Vector m_constraints = nullptr;
+            N_Vector m_func_tmpl   = nullptr;
 
-            /**
+            SUNMatrix m_J          = nullptr;
-             * @brief Constructs an EigenFunctor.
+            SUNLinearSolver m_LS   = nullptr;
-             *
-             * @param view The MultiscalePartitioningEngineView instance.
-             * @param qse_solve_species Species to solve for in the QSE group.
-             * @param initial_comp Initial abundances of all species in the full network.
-             * @param T9 Temperature in units of 10^9 K.
-             * @param rho Density in g/cm^3.
-             * @param Y_scale Scaling factors for the species abundances.
-             * @param qse_solve_species_index_map Mapping from species to their indices in the QSE solve vector.
-             */
-            EigenFunctor(
-                const MultiscalePartitioningEngineView& view,
-                const std::set<fourdst::atomic::Species>& qse_solve_species,
-                const fourdst::composition::CompositionAbstract& initial_comp,
-                const double T9,
-                const double rho,
-                const Eigen::VectorXd& Y_scale,
-                const std::unordered_map<fourdst::atomic::Species, size_t>& qse_solve_species_index_map
-            ) :
-            m_view(view),
-            m_qse_solve_species(qse_solve_species),
-            m_initial_comp(initial_comp),
-            m_T9(T9),
-            m_rho(rho),
-            m_Y_scale(Y_scale),
-            m_qse_solve_species_index_map(qse_solve_species_index_map) {}
 
-            /**
+            static quill::Logger* getLogger() {
-             * @brief Gets the number of output values from the functor (size of the residual vector).
+                return LogManager::getInstance().getLogger("log");
-             * @return The number of algebraic species being solved.
+            }
-             */
-            [[nodiscard]] size_t values() const { return m_qse_solve_species.size(); }
-            /**
-             * @brief Gets the number of input values to the functor (size of the variable vector).
-             * @return The number of algebraic species being solved.
-             */
-            [[nodiscard]] size_t inputs() const { return m_qse_solve_species.size(); }
 
-            /**
-             * @brief Evaluates the functor's residual vector `f_qse = dY_alg/dt`.
-             * @param v_qse The input vector of scaled algebraic abundances.
-             * @param f_qse The output residual vector.
-             * @return 0 on success.
-             */
-            int operator()(const InputType& v_qse, OutputType& f_qse) const;
-            /**
-             * @brief Evaluates the Jacobian of the functor, `J_qse = d(f_qse)/d(v_qse)`.
-             * @param v_qse The input vector of scaled algebraic abundances.
-             * @param J_qse The output Jacobian matrix.
-             * @return 0 on success.
-             */
-            int df(const InputType& v_qse, JacobianType& J_qse) const;
         };
 
         struct FluxValidationResult {
@@ -746,6 +713,11 @@ namespace gridfire {
          * @brief The list of identified equilibrium groups.
          */
         std::vector<QSEGroup> m_qse_groups;
+
+        /**
+         * @brief A set of solvers, one for each QSE group
+         */
+        std::vector<std::unique_ptr<QSESolver>> m_qse_solvers;
         /**
          * @brief The simplified set of species presented to the solver (the "slow" species).
          */
@@ -771,6 +743,8 @@ namespace gridfire {
 
         mutable std::unordered_map<uint64_t, fourdst::composition::Composition> m_composition_cache;
 
+        SUNContext m_sun_ctx = nullptr; // TODO: initialize and safely destroy this
+
     private:
         /**
          * @brief Partitions the network by timescale.

src/include/gridfire/utils/sundials.h

@@ -0,0 +1,66 @@
+#pragma once
+#include <unordered_map>
+#include <nvector/nvector_serial.h>
+
+#include "gridfire/exceptions/error_solver.h"
+#include "sundials/sundials_nvector.h"
+
+namespace gridfire::utils {
+    inline 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."},
+        {2, "CV_ROOT_RETURN: A root was found."},
+        {-99, "CV_WARNING: CVODE succeeded but in an unusual manner"},
+        {-1, "CV_TOO_MUCH_WORK: The solver took too many internal steps."},
+        {-2, "CV_TOO_MUCH_ACC: The solver could not satisfy the accuracy requested."},
+        {-3, "CV_ERR_FAILURE: The solver encountered a non-recoverable error."},
+        {-4, "CV_CONV_FAILURE: The solver failed to converge."},
+        {-5, "CV_LINIT_FAIL: The linear solver's initialization function failed."},
+        {-6, "CV_LSETUP_FAIL: The linear solver's setup function failed."},
+        {-7, "CV_LSOLVE_FAIL: The linear solver's solve function failed."},
+        {-8, "CV_RHSFUNC_FAIL: The right-hand side function failed in an unrecoverable manner."},
+        {-9, "CV_FIRST_RHSFUNC_ERR: The right-hand side function failed at the first call."},
+        {-10, "CV_REPTD_RHSFUNC_ERR: The right-hand side function repeatedly failed recoverable."},
+        {-11, "CV_UNREC_RHSFUNC_ERR: The right-hand side function failed unrecoverably."},
+        {-12, "CV_RTFUNC_FAIL: The rootfinding function failed in an unrecoverable manner."},
+        {-13, "CV_NLS_INIT_FAIL: The nonlinear solver's initialization function failed."},
+        {-14, "CV_NLS_SETUP_FAIL: The nonlinear solver's setup function failed."},
+        {-15, "CV_CONSTR_FAIL : The inequality constraint was violated and the solver was unable to recover."},
+        {-16, "CV_NLS_FAIL: The nonlinear solver's solve function failed."},
+        {-20, "CV_MEM_FAIL: Memory allocation failed."},
+        {-21, "CV_MEM_NULL: The CVODE memory structure is NULL."},
+        {-22, "CV_ILL_INPUT: An illegal input was detected."},
+        {-23, "CV_NO_MALLOC: The CVODE memory structure has not been allocated."},
+        {-24, "CV_BAD_K: The value of k is invalid."},
+        {-25, "CV_BAD_T: The value of t is invalid."},
+        {-26, "CV_BAD_DKY: The value of dky is invalid."},
+        {-27, "CV_TOO_CLOSE: The time points are too close together."},
+        {-28, "CV_VECTOROP_ERR: A vector operation failed."},
+        {-29, "CV_PROJ_MEM_NULL: The projection memory structure is NULL."},
+        {-30, "CV_PROJFUNC_FAIL: The projection function failed in an unrecoverable manner."},
+        {-31, "CV_REPTD_PROJFUNC_ERR: The projection function has repeated recoverable errors."}
+    };
+
+    inline void check_cvode_flag(const int flag, const std::string& func_name) {
+        if (flag < 0) {
+            if (!cvode_ret_code_map.contains(flag)) {
+                throw exceptions::CVODESolverFailureError("CVODE error in " + func_name + ": Unknown error code: " + std::to_string(flag));
+            }
+            throw exceptions::CVODESolverFailureError("CVODE error in " + func_name + ": " + cvode_ret_code_map.at(flag));
+        }
+    }
+
+    inline N_Vector init_sun_vector(uint64_t size, SUNContext sun_ctx) {
+        #ifdef SUNDIALS_HAVE_OPENMP
+            const N_Vector vec = N_VNew_OpenMP(size, 0, sun_ctx);
+        #elif SUNDIALS_HAVE_PTHREADS
+            const N_Vector vec = N_VNew_Pthreads(size, sun_ctx);
+        #else
+            const N_Vector vec = N_VNew_Serial(static_cast<long long>(size), sun_ctx);
+        #endif
+
+        check_cvode_flag(vec == nullptr ? -1 : 0, "N_VNew");
+        return vec;
+    }
+
+}

src/lib/engine/views/engine_multiscale.cpp

@@ -1,6 +1,7 @@
 #include "gridfire/engine/views/engine_multiscale.h"
 #include "gridfire/exceptions/error_engine.h"
 #include "gridfire/engine/procedures/priming.h"
+#include "gridfire/utils/sundials.h"
 
 #include <stdexcept>
 #include <vector>
@@ -17,12 +18,15 @@
 #include "quill/LogMacros.h"
 #include "quill/Logger.h"
 
+#include "kinsol/kinsol.h"
+#include "sundials/sundials_context.h"
+#include "sunmatrix/sunmatrix_dense.h"
+#include "sunlinsol/sunlinsol_dense.h"
+
 #include  "xxhash64.h"
 #include "fourdst/composition/utils/composition_hash.h"
 
 namespace {
-    constexpr double MIN_ABS_NORM_ALWAYS_CONVERGED = 1e-30;
-
     using namespace fourdst::atomic;
 
 
@@ -159,6 +163,18 @@ namespace {
         {Eigen::LevenbergMarquardtSpace::Status::GtolTooSmall, "GtolTooSmall"},
         {Eigen::LevenbergMarquardtSpace::Status::UserAsked, "UserAsked"}
     };
+
+    void QuietErrorRouter(int line, const char *func, const char *file, const char *msg,
+                      SUNErrCode err_code, void *err_user_data, SUNContext sunctx) {
+
+        // LIST OF ERRORS TO IGNORE
+        if (err_code == KIN_LINESEARCH_NONCONV) {
+            return;
+        }
+
+        // For everything else, use the default SUNDIALS logger (or your own)
+        SUNLogErrHandlerFn(line, func, file, msg, err_code, err_user_data, sunctx);
+    }
 }
 
 namespace gridfire {
@@ -166,7 +182,23 @@ namespace gridfire {
 
     MultiscalePartitioningEngineView::MultiscalePartitioningEngineView(
         DynamicEngine& baseEngine
-    ) : m_baseEngine(baseEngine) {}
+    ) : m_baseEngine(baseEngine) {
+        const int flag = SUNContext_Create(SUN_COMM_NULL, &m_sun_ctx);
+        if (flag != 0) {
+            LOG_CRITICAL(m_logger, "Error while creating SUNContext in MultiscalePartitioningEngineView");
+            throw std::runtime_error("Error creating SUNContext in MultiscalePartitioningEngineView");
+        }
+        SUNContext_PushErrHandler(m_sun_ctx, QuietErrorRouter, nullptr);
+    }
+
+    MultiscalePartitioningEngineView::~MultiscalePartitioningEngineView() {
+        LOG_TRACE_L1(m_logger, "Cleaning up MultiscalePartitioningEngineView...");
+        m_qse_solvers.clear();
+        if (m_sun_ctx) {
+            SUNContext_Free(&m_sun_ctx);
+            m_sun_ctx = nullptr;
+        }
+    }
 
     const std::vector<Species> & MultiscalePartitioningEngineView::getNetworkSpecies() const {
         return m_baseEngine.getNetworkSpecies();
@@ -598,6 +630,7 @@ namespace gridfire {
         LOG_TRACE_L1(m_logger, "Partitioning network...");
         LOG_TRACE_L1(m_logger, "Clearing previous state...");
         m_qse_groups.clear();
+        m_qse_solvers.clear();
         m_dynamic_species.clear();
         m_algebraic_species.clear();
         m_composition_cache.clear(); // We need to clear the cache now cause the same comp, temp, and density may result in a different value
@@ -790,7 +823,18 @@ namespace gridfire {
                 m_dynamic_species.push_back(species);
             }
         }
-        return getNormalizedEquilibratedComposition(comp, T9, rho);
+
+        for (const auto& group : m_qse_groups) {
+            std::vector<Species> groupAlgebraicSpecies;
+            for (const auto& species : group.algebraic_species) {
+                groupAlgebraicSpecies.push_back(species);
+            }
+            m_qse_solvers.push_back(std::make_unique<QSESolver>(groupAlgebraicSpecies, m_baseEngine, m_sun_ctx));
+        }
+
+        fourdst::composition::Composition result = getNormalizedEquilibratedComposition(comp, T9, rho);
+
+        return result;
     }
 
     void MultiscalePartitioningEngineView::exportToDot(
@@ -1502,98 +1546,18 @@ namespace gridfire {
         const double T9,
         const double rho
     ) const {
-        LOG_TRACE_L1(m_logger, "Solving for QSE abundances...");
+        LOG_TRACE_L2(m_logger, "Solving for QSE abundances...");
-        LOG_TRACE_L2(m_logger, "Composition before QSE solving: {}", [&comp]() -> std::string {
-            std::stringstream ss;
-            size_t i = 0;
-            for (const auto& [sp, y] : comp) {
-                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) {
+        for (const auto& [group, solver]: std::views::zip(m_qse_groups, m_qse_solvers)) {
-            LOG_TRACE_L2(m_logger, "Working on QSE group with algebraic species: {}",
+            const fourdst::composition::Composition groupResult = solver->solve(outputComposition, T9, rho);
-                [&]() -> std::string {
+            for (const auto& [sp, y] : groupResult) {
-                    std::stringstream ss;
+                outputComposition.setMolarAbundance(sp, y);
-                    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;
-            }
-
-            Eigen::VectorXd Y_scale(algebraic_species.size());
-            Eigen::VectorXd v_initial(algebraic_species.size());
-            long i = 0;
-            std::unordered_map<Species, size_t> species_to_index_map;
-            for (const auto& species : algebraic_species) {
-                constexpr double abundance_floor = 1.0e-100;
-                const double initial_abundance = comp.getMolarAbundance(species);
-                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_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;
-                msg << "While working on QSE group with algebraic species: ";
-                size_t count = 0;
-                for (const auto& species: algebraic_species) {
-                    msg << species;
-                    if (count < algebraic_species.size() - 1) {
-                        msg << ", ";
-                    }
-                    count++;
-                }
-                msg << " the QSE solver failed to converge with status: " << lm_status_map.at(status) << " (No. " << status << ").";
-                LOG_ERROR(m_logger, "{}", msg.str());
-                throw std::runtime_error(msg.str());
-            }
-            LOG_TRACE_L1(m_logger, "QSE Group minimization succeeded with status: {}", lm_status_map.at(status));
-            Eigen::VectorXd Y_final_qse = v_initial.array().exp(); // Convert back to physical abundances using exponential scaling
-            i = 0;
-            for (const auto& species: algebraic_species) {
-                LOG_TRACE_L1(
-                    m_logger,
-                    "During QSE solving species {} started with a molar abundance of {} and ended with an abundance of {}.",
-                    species.name(),
-                    comp.getMolarAbundance(species),
-                    Y_final_qse(i)
-                );
-                // double Xi = Y_final_qse(i) * species.mass(); // Convert from molar abundance to mass fraction
-                if (!outputComposition.contains (species)) {
-                    outputComposition.registerSpecies(species);
-                }
-                outputComposition.setMolarAbundance(species, Y_final_qse(i));
-                i++;
             }
+            solver->log_diagnostics();
         }
+        LOG_TRACE_L2(m_logger, "Done solving for QSE abundances!");
         return outputComposition;
     }
 
@@ -1795,121 +1759,249 @@ namespace gridfire {
         return candidate_groups;
     }
 
-
+     //////////////////////////////////
-     //////////////////////////////////////
+    /// QSESolver Member Functions ///
-    /// Eigen Functor Member Functions ///
+    //////////////////////////////////
-   /////////////////////////////////////
 
 
-    int MultiscalePartitioningEngineView::EigenFunctor::operator()(const InputType &v_qse, OutputType &f_qse) const {
+    MultiscalePartitioningEngineView::QSESolver::QSESolver(
-        fourdst::composition::Composition comp_trial(m_initial_comp.getRegisteredSpecies());
+        const std::vector<fourdst::atomic::Species>& species,
-        for (const auto& [sp, y] : m_initial_comp) {
+        const DynamicEngine& engine,
-            comp_trial.setMolarAbundance(sp, y);
+        const SUNContext sun_ctx
-        }
+    ) :
-        Eigen::VectorXd y_qse = v_qse.array().exp(); // Convert to physical abundances using exponential scaling
+    m_N(species.size()),
+    m_engine(engine),
+    m_species(species),
+    m_sun_ctx(sun_ctx) {
+        m_Y = utils::init_sun_vector(m_N, m_sun_ctx);
+        m_scale = N_VClone(m_Y);
+        m_f_scale = N_VClone(m_Y);
+        m_constraints = N_VClone(m_Y);
+        m_func_tmpl = N_VClone(m_Y);
 
-        for (const auto& species: m_qse_solve_species) {
+        if (!m_Y || !m_scale || !m_constraints || !m_func_tmpl) {
-            auto index = static_cast<long>(m_qse_solve_species_index_map.at(species));
+            LOG_CRITICAL(getLogger(), "Failed to allocate SUNVectors for QSE solver.");
-            comp_trial.setMolarAbundance(species, y_qse[index]);
+            throw std::runtime_error("Failed to allocate SUNVectors for QSE solver.");
         }
 
-        const auto result = m_view.getBaseEngine().calculateRHSAndEnergy(comp_trial, m_T9, m_rho);
+        for (size_t i = 0; i < m_N; ++i) {
-        if (!result) {
+            m_speciesMap[m_species[i]] = i;
-            throw exceptions::StaleEngineError("Failed to calculate RHS and energy due to stale engine state");
         }
-        const auto&[dydt, nuclearEnergyGenerationRate, _] = result.value();
+
-        f_qse.resize(static_cast<long>(m_qse_solve_species.size()));
+        N_VConst(1.0, m_constraints);
-        long i = 0;
+        m_kinsol_mem = KINCreate(m_sun_ctx);
-        // TODO: make sure that just counting up i is a valid approach, this is a possible place an indexing bug may have crept in
+
-        for (const auto& species: m_qse_solve_species) {
+        utils::check_cvode_flag(m_kinsol_mem ? 0 : -1, "KINCreate");
-            const double dydt_i = dydt.at(species);
+        utils::check_cvode_flag(KINInit(m_kinsol_mem, sys_func, m_func_tmpl), "KINInit");
-            f_qse(i) = dydt_i/y_qse(i); // We square the residuals to improve numerical stability in the solver
+        utils::check_cvode_flag(KINSetConstraints(m_kinsol_mem, m_constraints), "KINSetConstraints");
-            i++;
+
+        m_J = SUNDenseMatrix(static_cast<sunindextype>(m_N), static_cast<sunindextype>(m_N), m_sun_ctx);
+        utils::check_cvode_flag(m_J ? 0 : -1, "SUNDenseMatrix");
+
+        m_LS = SUNLinSol_Dense(m_Y, m_J, m_sun_ctx);
+        utils::check_cvode_flag(m_LS ? 0 : -1, "SUNLinSol_Dense");
+
+        utils::check_cvode_flag(KINSetLinearSolver(m_kinsol_mem, m_LS, m_J), "KINSetLinearSolver");
+        utils::check_cvode_flag(KINSetJacFn(m_kinsol_mem, sys_jac), "KINSetJacFn");
+
+        utils::check_cvode_flag(KINSetMaxSetupCalls(m_kinsol_mem, 20), "KINSetMaxSetupCalls");
+
+        utils::check_cvode_flag(KINSetFuncNormTol(m_kinsol_mem, 1e-6), "KINSetFuncNormTol");
+        utils::check_cvode_flag(KINSetNumMaxIters(m_kinsol_mem, 200), "KINSetNumMaxIters");
+
+        // We want to effectively disable this since enormous changes in order of magnitude are realistic for this problem.
+        utils::check_cvode_flag(KINSetMaxNewtonStep(m_kinsol_mem, 200), "KINSetMaxNewtonStep");
+
     }
-        LOG_TRACE_L2(
+
-            m_view.m_logger,
+    MultiscalePartitioningEngineView::QSESolver::~QSESolver() {
-            "Functor evaluation at T9 = {}, rho = {}, y_qse (v_qse) = <{}> complete. ||f|| = {}",
+        if (m_Y) {
-            m_T9,
+            N_VDestroy(m_Y);
-            m_rho,
+            m_Y = nullptr;
-            [&]() -> std::string {
+        }
-                std::stringstream ss;
+        if (m_scale) {
-                for (long j = 0; j < y_qse.size(); ++j) {
+            N_VDestroy(m_scale);
-                    ss << y_qse(j);
+            m_scale = nullptr;
-                    ss << "(" << v_qse(j) << ")";
+        }
-                    if (j < y_qse.size() - 1) {
+        if (m_f_scale) {
-                        ss << ", ";
+            N_VDestroy(m_f_scale);
+            m_f_scale = nullptr;
+        }
+        if (m_constraints) {
+            N_VDestroy(m_constraints);
+            m_constraints = nullptr;
+        }
+        if (m_func_tmpl) {
+            N_VDestroy(m_func_tmpl);
+            m_func_tmpl = nullptr;
+        }
+        if (m_kinsol_mem) {
+            KINFree(&m_kinsol_mem);
+            m_kinsol_mem = nullptr;
+        }
+        if (m_J) {
+            SUNMatDestroy(m_J);
+            m_J = nullptr;
+        }
+        if (m_LS) {
+            SUNLinSolFree(m_LS);
+            m_LS = nullptr;
         }
     }
-                return ss.str();
+
-            }(),
+    fourdst::composition::Composition MultiscalePartitioningEngineView::QSESolver::solve(
-            f_qse.norm()
+        const fourdst::composition::Composition &comp,
+        const double T9,
+        const double rho
+    ) const {
+
+        fourdst::composition::Composition result = comp;
+        UserData data {
+            m_engine,
+            T9,
+            rho,
+            result,
+            m_speciesMap,
+            m_species
+        };
+
+        utils::check_cvode_flag(KINSetUserData(m_kinsol_mem, &data), "KINSetUserData");
+
+        sunrealtype* y_data = N_VGetArrayPointer(m_Y);
+        sunrealtype* scale_data = N_VGetArrayPointer(m_scale);
+
+        // It is more cache optimized to do a standard as opposed to range based for-loop here
+        for (size_t i = 0; i < m_N; ++i) {
+            const auto& species = m_species[i];
+            double Y = result.getMolarAbundance(species);
+
+            constexpr double abundance_floor = 1.0e-100;
+            Y = std::max(abundance_floor, Y);
+            y_data[i] = Y;
+            scale_data[i] = 1.0;
+        }
+
+        auto initial_rhs = m_engine.calculateRHSAndEnergy(result, T9, rho);
+        if (!initial_rhs) {
+            throw std::runtime_error("In QSE solver failed to calculate initial RHS");
+        }
+
+        sunrealtype* f_scale_data = N_VGetArrayPointer(m_f_scale);
+        for (size_t i = 0; i < m_N; ++i) {
+            const auto& species = m_species[i];
+            double dydt = std::abs(initial_rhs.value().dydt.at(species));
+            f_scale_data[i] = 1.0 / (dydt + 1e-15);
+        }
+
+        if (m_solves > 0) {
+            // After the initial solution we want to allow kinsol to reuse its state
+            utils::check_cvode_flag(KINSetNoInitSetup(m_kinsol_mem, SUNTRUE), "KINSetNoInitSetup");
+        } else {
+            utils::check_cvode_flag(KINSetNoInitSetup(m_kinsol_mem, SUNFALSE), "KINSetNoInitSetup");
+        }
+
+        const int flag = KINSol(m_kinsol_mem, m_Y, KIN_LINESEARCH, m_scale, m_f_scale);
+        if (flag < 0) {
+            LOG_WARNING(getLogger(), "KINSol failed to converge while solving QSE abundances with flag {}.", utils::cvode_ret_code_map.at(flag));
+            return comp;
+        }
+
+        for (size_t i = 0; i < m_N; ++i) {
+            const auto& species = m_species[i];
+            result.setMolarAbundance(species, y_data[i]);
+        }
+
+        m_solves++;
+        return result;
+    }
+
+    size_t MultiscalePartitioningEngineView::QSESolver::solves() const {
+        return m_solves;
+    }
+
+    void MultiscalePartitioningEngineView::QSESolver::log_diagnostics() const {
+        long int nni, nfe, nje;
+
+        int flag = KINGetNumNonlinSolvIters(m_kinsol_mem, &nni);
+
+        flag = KINGetNumFuncEvals(m_kinsol_mem, &nfe);
+
+        flag = KINGetNumJacEvals(m_kinsol_mem, &nje);
+
+        LOG_INFO(getLogger(),
+            "QSE Stats | Iters: {} | RHS Evals: {} | Jac Evals: {} | Ratio (J/I): {:.2f}",
+            nni, nfe, nje, static_cast<double>(nje) / static_cast<double>(nni)
         );
-        LOG_TRACE_L3(
+        getLogger()->flush_log(true);
-            m_view.m_logger,
-            "{}",
-            [&]() -> std::string {
-                std::stringstream ss;
-                const std::vector species(m_qse_solve_species.begin(), m_qse_solve_species.end());
-                for (long j = 0; j < f_qse.size(); ++j) {
-                    ss << "Residual for species " << species.at(j).name() << " f(" << j << ") = " << f_qse(j) << "\n";
-                }
-                return ss.str();
-            }()
-        );
-        return 0;
     }
 
-    int MultiscalePartitioningEngineView::EigenFunctor::df(const InputType &v_qse, JacobianType &J_qse) const {
-        fourdst::composition::Composition comp_trial(m_initial_comp.getRegisteredSpecies());
-        for (const auto& [sp, y] : m_initial_comp) {
-            comp_trial.setMolarAbundance(sp, y);
-        }
-        Eigen::VectorXd y_qse = v_qse.array().exp(); // Convert to physical abundances using exponential scaling
 
-        for (const auto& species: m_qse_solve_species) {
+    int MultiscalePartitioningEngineView::QSESolver::sys_func(
-            const double molarAbundance = y_qse[static_cast<long>(m_qse_solve_species_index_map.at(species))];
+        const N_Vector y,
-            comp_trial.setMolarAbundance(species, molarAbundance);
+        const N_Vector f,
+        void *user_data
+    ) {
+        const auto* data = static_cast<UserData*>(user_data);
+
+        const sunrealtype* y_data = N_VGetArrayPointer(y);
+        sunrealtype* f_data = N_VGetArrayPointer(f);
+
+        const auto& map = data->qse_solve_species_index_map;
+        for (size_t index = 0; index < data->qse_solve_species.size(); ++index) {
+            const auto& species = data->qse_solve_species[index];
+            data->comp.setMolarAbundance(species, y_data[index]);
         }
 
-        std::vector<Species> qse_species_vector(m_qse_solve_species.begin(), m_qse_solve_species.end());
+        const auto result = data->engine.calculateRHSAndEnergy(data->comp, data->T9, data->rho);
-        NetworkJacobian jac = m_view.getBaseEngine().generateJacobianMatrix(comp_trial, m_T9, m_rho, qse_species_vector);
+
-        const auto result = m_view.getBaseEngine().calculateRHSAndEnergy(comp_trial, m_T9, m_rho);
         if (!result) {
-            throw exceptions::StaleEngineError("Failed to calculate RHS and energy due to stale engine state");
+            return 1; // Potentially recoverable error
-        }
-        const auto&[dydt, nuclearEnergyGenerationRate, _] = result.value();
-
-        const long N = static_cast<long>(m_qse_solve_species.size());
-        J_qse.resize(N, N);
-        long rowID = 0;
-        for (const auto& rowSpecies : m_qse_solve_species) {
-            long colID = 0;
-            for (const auto& colSpecies: m_qse_solve_species) {
-                J_qse(rowID, colID) = jac(rowSpecies, colSpecies);
-                colID += 1;
-                LOG_TRACE_L3(m_view.m_logger, "Jacobian[{}, {}] (d(dY({}))/dY({})) = {}", rowID, colID - 1, rowSpecies.name(), colSpecies.name(), J_qse(rowID, colID - 1));
-            }
-            rowID += 1;
         }
 
-        for (long i = 0; i < J_qse.rows(); ++i) {
+        const auto& dydt = result.value().dydt;
-            for (long j = 0; j < J_qse.cols(); ++j) {
-                double on_diag_correction = 0.0;
-                if (i == j) {
-                    auto rowSpecies = *(std::next(m_qse_solve_species.begin(), i));
-                    const double Fi = dydt.at(rowSpecies);
-                    on_diag_correction = Fi / y_qse(i);
-                }
-                J_qse(i, j) = y_qse(j) * (J_qse(i, j) - on_diag_correction) / y_qse(i); // Apply chain rule J'(i,j) = y_j * (J(i,j) - δ_ij(F_i/Y_i)) / Y_i
-            }
-        }
 
-        m_cached_jacobian = J_qse; // Cache the computed Jacobian for future use
+        for (const auto& [species, index] : map) {
+            f_data[index] = dydt.at(species);
+        }
 
         return 0; // Success
     }
 
 
+    int MultiscalePartitioningEngineView::QSESolver::sys_jac(
+        const N_Vector y,
+        N_Vector fy,
+        SUNMatrix J,
+        void *user_data,
+        N_Vector tmp1,
+        N_Vector tmp2
+    ) {
+        const auto* data = static_cast<UserData*>(user_data);
+
+        const sunrealtype* y_data = N_VGetArrayPointer(y);
+        const auto& map = data->qse_solve_species_index_map;
+        for (const auto& [species, index] : map) {
+            data->comp.setMolarAbundance(species, y_data[index]);
+        }
+
+        const NetworkJacobian jac = data->engine.generateJacobianMatrix(
+            data->comp,
+            data->T9,
+            data->rho,
+            data->qse_solve_species
+        );
+
+        sunrealtype* J_data = SUNDenseMatrix_Data(J);
+        const sunindextype N = SUNDenseMatrix_Columns(J);
+
+        for (const auto& [col_species, col_idx] : map) {
+            for (const auto& [row_species, row_idx] : map) {
+                J_data[col_idx * N + row_idx] = jac(row_species, col_species);
+            }
+        }
+
+        return 0;
+    }
 
       /////////////////////////////////
      /// QSEGroup Member Functions ///

src/lib/solver/strategies/CVODE_solver_strategy.cpp

@@ -22,65 +22,8 @@
 #include "gridfire/solver/strategies/triggers/engine_partitioning_trigger.h"
 #include "gridfire/trigger/procedures/trigger_pprint.h"
 #include "gridfire/exceptions/error_solver.h"
-#include "gridfire/utils/logging.h"
+#include "gridfire/utils/sundials.h"
 
-namespace {
-
-    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."},
-        {2, "CV_ROOT_RETURN: A root was found."},
-        {-99, "CV_WARNING: CVODE succeeded but in an unusual manner"},
-        {-1, "CV_TOO_MUCH_WORK: The solver took too many internal steps."},
-        {-2, "CV_TOO_MUCH_ACC: The solver could not satisfy the accuracy requested."},
-        {-3, "CV_ERR_FAILURE: The solver encountered a non-recoverable error."},
-        {-4, "CV_CONV_FAILURE: The solver failed to converge."},
-        {-5, "CV_LINIT_FAIL: The linear solver's initialization function failed."},
-        {-6, "CV_LSETUP_FAIL: The linear solver's setup function failed."},
-        {-7, "CV_LSOLVE_FAIL: The linear solver's solve function failed."},
-        {-8, "CV_RHSFUNC_FAIL: The right-hand side function failed in an unrecoverable manner."},
-        {-9, "CV_FIRST_RHSFUNC_ERR: The right-hand side function failed at the first call."},
-        {-10, "CV_REPTD_RHSFUNC_ERR: The right-hand side function repeatedly failed recoverable."},
-        {-11, "CV_UNREC_RHSFUNC_ERR: The right-hand side function failed unrecoverably."},
-        {-12, "CV_RTFUNC_FAIL: The rootfinding function failed in an unrecoverable manner."},
-        {-13, "CV_NLS_INIT_FAIL: The nonlinear solver's initialization function failed."},
-        {-14, "CV_NLS_SETUP_FAIL: The nonlinear solver's setup function failed."},
-        {-15, "CV_CONSTR_FAIL : The inequality constraint was violated and the solver was unable to recover."},
-        {-16, "CV_NLS_FAIL: The nonlinear solver's solve function failed."},
-        {-20, "CV_MEM_FAIL: Memory allocation failed."},
-        {-21, "CV_MEM_NULL: The CVODE memory structure is NULL."},
-        {-22, "CV_ILL_INPUT: An illegal input was detected."},
-        {-23, "CV_NO_MALLOC: The CVODE memory structure has not been allocated."},
-        {-24, "CV_BAD_K: The value of k is invalid."},
-        {-25, "CV_BAD_T: The value of t is invalid."},
-        {-26, "CV_BAD_DKY: The value of dky is invalid."},
-        {-27, "CV_TOO_CLOSE: The time points are too close together."},
-        {-28, "CV_VECTOROP_ERR: A vector operation failed."},
-        {-29, "CV_PROJ_MEM_NULL: The projection memory structure is NULL."},
-        {-30, "CV_PROJFUNC_FAIL: The projection function failed in an unrecoverable manner."},
-        {-31, "CV_REPTD_PROJFUNC_ERR: The projection function has repeated recoverable errors."}
-    };
-    void check_cvode_flag(const int flag, const std::string& func_name) {
-        if (flag < 0) {
-            if (!cvode_ret_code_map.contains(flag)) {
-                throw gridfire::exceptions::CVODESolverFailureError("CVODE error in " + func_name + ": Unknown error code: " + std::to_string(flag));
-            }
-            throw gridfire::exceptions::CVODESolverFailureError("CVODE error in " + func_name + ": " + cvode_ret_code_map.at(flag));
-        }
-    }
-
-    N_Vector init_sun_vector(uint64_t size, SUNContext sun_ctx) {
-#ifdef SUNDIALS_HAVE_OPENMP
-        const N_Vector vec = N_VNew_OpenMP(size, 0, sun_ctx);
-#elif SUNDIALS_HAVE_PTHREADS
-        const N_Vector vec = N_VNew_Pthreads(size, sun_ctx);
-#else
-        const N_Vector vec = N_VNew_Serial(static_cast<long long>(size), sun_ctx);
-#endif
-        check_cvode_flag(vec == nullptr ? -1 : 0, "N_VNew");
-        return vec;
-    }
-}
 
 namespace gridfire::solver {
 
@@ -175,7 +118,7 @@ namespace gridfire::solver {
         LOG_TRACE_L1(m_logger, "Number of species: {} ({} independent variables)", numSpecies, N);
         LOG_TRACE_L1(m_logger, "Initializing CVODE resources");
         m_cvode_mem = CVodeCreate(CV_BDF, m_sun_ctx);
-        check_cvode_flag(m_cvode_mem == nullptr ? -1 : 0, "CVodeCreate");
+        utils::check_cvode_flag(m_cvode_mem == nullptr ? -1 : 0, "CVodeCreate");
 
         initialize_cvode_integration_resources(N, numSpecies, 0.0, equilibratedComposition, absTol, relTol, 0.0);
 
@@ -206,11 +149,11 @@ namespace gridfire::solver {
             user_data.networkSpecies = &m_engine.getNetworkSpecies();
             user_data.captured_exception.reset();
 
-            check_cvode_flag(CVodeSetUserData(m_cvode_mem, &user_data), "CVodeSetUserData");
+            utils::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));
+            LOG_TRACE_L2(m_logger, "CVODE step complete. Current time: {}, step status: {}", current_time, utils::cvode_ret_code_map.at(flag));
 
             if (user_data.captured_exception){
                 std::rethrow_exception(std::make_exception_ptr(*user_data.captured_exception));
@@ -220,7 +163,7 @@ namespace gridfire::solver {
             // TODO: Come up with some way to save these to a file rather than spamming stdout. JSON maybe? OPAT?
             // log_step_diagnostics(user_data, true, false, false);
             // exit(0);
-            check_cvode_flag(flag, "CVode");
+            utils::check_cvode_flag(flag, "CVode");
 
             long int n_steps;
             double last_step_size;
@@ -449,11 +392,11 @@ namespace gridfire::solver {
                 cleanup_cvode_resources(true);
 
                 m_cvode_mem = CVodeCreate(CV_BDF, m_sun_ctx);
-                check_cvode_flag(m_cvode_mem == nullptr ? -1 : 0, "CVodeCreate");
+                utils::check_cvode_flag(m_cvode_mem == nullptr ? -1 : 0, "CVodeCreate");
 
                 initialize_cvode_integration_resources(N, numSpecies, current_time, currentComposition, absTol, relTol, accumulated_energy);
 
-                check_cvode_flag(CVodeReInit(m_cvode_mem, current_time, m_Y), "CVodeReInit");
+                utils::check_cvode_flag(CVodeReInit(m_cvode_mem, current_time, m_Y), "CVodeReInit");
             }
 
         }
@@ -487,7 +430,7 @@ namespace gridfire::solver {
         NetOut netOut;
         netOut.composition = outputComposition;
         netOut.energy = accumulated_energy;
-        check_cvode_flag(CVodeGetNumSteps(m_cvode_mem, reinterpret_cast<long int *>(&netOut.num_steps)), "CVodeGetNumSteps");
+        utils::check_cvode_flag(CVodeGetNumSteps(m_cvode_mem, reinterpret_cast<long int *>(&netOut.num_steps)), "CVodeGetNumSteps");
 
         LOG_TRACE_L2(m_logger, "generating final nuclear energy generation rate derivatives...");
         auto [dEps_dT, dEps_dRho] = m_engine.calculateEpsDerivatives(
@@ -688,7 +631,7 @@ 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());
+        LOG_TRACE_L2(m_logger, "Calculating RHS at time {} with {} species in composition", t, composition.size());
         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);
@@ -732,7 +675,7 @@ namespace gridfire::solver {
         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);
+        m_Y = utils::init_sun_vector(N, m_sun_ctx);
         m_YErr = N_VClone(m_Y);
 
         sunrealtype *y_data = N_VGetArrayPointer(m_Y);
@@ -747,8 +690,8 @@ namespace gridfire::solver {
         y_data[numSpecies] = accumulatedEnergy; // Specific energy rate, initialized to zero
 
 
-        check_cvode_flag(CVodeInit(m_cvode_mem, cvode_rhs_wrapper, current_time, m_Y), "CVodeInit");
+        utils::check_cvode_flag(CVodeInit(m_cvode_mem, cvode_rhs_wrapper, current_time, m_Y), "CVodeInit");
-        check_cvode_flag(CVodeSStolerances(m_cvode_mem, relTol, absTol), "CVodeSStolerances");
+        utils::check_cvode_flag(CVodeSStolerances(m_cvode_mem, relTol, absTol), "CVodeSStolerances");
 
         // Constraints
         // We constrain the solution vector using CVODE's built in constraint flags as outlines on page 53 of the CVODE manual
@@ -768,17 +711,17 @@ namespace gridfire::solver {
         }
         N_VConst(1.0, m_constraints); // Set all constraints to >= 0 (note this is where the flag values are set)
 
-        check_cvode_flag(CVodeSetConstraints(m_cvode_mem, m_constraints), "CVodeSetConstraints");
+        utils::check_cvode_flag(CVodeSetConstraints(m_cvode_mem, m_constraints), "CVodeSetConstraints");
 
-        check_cvode_flag(CVodeSetMaxStep(m_cvode_mem, 1.0e20), "CVodeSetMaxStep");
+        utils::check_cvode_flag(CVodeSetMaxStep(m_cvode_mem, 1.0e20), "CVodeSetMaxStep");
 
         m_J = SUNDenseMatrix(static_cast<sunindextype>(N), static_cast<sunindextype>(N), m_sun_ctx);
-        check_cvode_flag(m_J == nullptr ? -1 : 0, "SUNDenseMatrix");
+        utils::check_cvode_flag(m_J == nullptr ? -1 : 0, "SUNDenseMatrix");
         m_LS = SUNLinSol_Dense(m_Y, m_J, m_sun_ctx);
-        check_cvode_flag(m_LS == nullptr ? -1 : 0, "SUNLinSol_Dense");
+        utils::check_cvode_flag(m_LS == nullptr ? -1 : 0, "SUNLinSol_Dense");
 
-        check_cvode_flag(CVodeSetLinearSolver(m_cvode_mem, m_LS, m_J), "CVodeSetLinearSolver");
+        utils::check_cvode_flag(CVodeSetLinearSolver(m_cvode_mem, m_LS, m_J), "CVodeSetLinearSolver");
-        check_cvode_flag(CVodeSetJacFn(m_cvode_mem, cvode_jac_wrapper), "CVodeSetJacFn");
+        utils::check_cvode_flag(CVodeSetJacFn(m_cvode_mem, cvode_jac_wrapper), "CVodeSetJacFn");
         LOG_TRACE_L2(m_logger, "CVODE solver initialized");
     }
 
@@ -814,10 +757,10 @@ namespace gridfire::solver {
         sunrealtype hlast, hcur, tcur;
         int qlast;
 
-        check_cvode_flag(CVodeGetLastStep(m_cvode_mem, &hlast), "CVodeGetLastStep");
+        utils::check_cvode_flag(CVodeGetLastStep(m_cvode_mem, &hlast), "CVodeGetLastStep");
-        check_cvode_flag(CVodeGetCurrentStep(m_cvode_mem, &hcur), "CVodeGetCurrentStep");
+        utils::check_cvode_flag(CVodeGetCurrentStep(m_cvode_mem, &hcur), "CVodeGetCurrentStep");
-        check_cvode_flag(CVodeGetLastOrder(m_cvode_mem, &qlast), "CVodeGetLastOrder");
+        utils::check_cvode_flag(CVodeGetLastOrder(m_cvode_mem, &qlast), "CVodeGetLastOrder");
-        check_cvode_flag(CVodeGetCurrentTime(m_cvode_mem, &tcur), "CVodeGetCurrentTime");
+        utils::check_cvode_flag(CVodeGetCurrentTime(m_cvode_mem, &tcur), "CVodeGetCurrentTime");
 
         {
             std::vector<std::string> labels = {"Current Time (tcur)", "Last Step (hlast)", "Current Step (hcur)", "Last Order (qlast)"};
@@ -834,13 +777,13 @@ namespace gridfire::solver {
         // These are the CRITICAL counters for diagnosing your problem
         long int nsteps, nfevals, nlinsetups, netfails, nniters, nconvfails, nsetfails;
 
-        check_cvode_flag(CVodeGetNumSteps(m_cvode_mem, &nsteps), "CVodeGetNumSteps");
+        utils::check_cvode_flag(CVodeGetNumSteps(m_cvode_mem, &nsteps), "CVodeGetNumSteps");
-        check_cvode_flag(CVodeGetNumRhsEvals(m_cvode_mem, &nfevals), "CVodeGetNumRhsEvals");
+        utils::check_cvode_flag(CVodeGetNumRhsEvals(m_cvode_mem, &nfevals), "CVodeGetNumRhsEvals");
-        check_cvode_flag(CVodeGetNumLinSolvSetups(m_cvode_mem, &nlinsetups), "CVodeGetNumLinSolvSetups");
+        utils::check_cvode_flag(CVodeGetNumLinSolvSetups(m_cvode_mem, &nlinsetups), "CVodeGetNumLinSolvSetups");
-        check_cvode_flag(CVodeGetNumErrTestFails(m_cvode_mem, &netfails), "CVodeGetNumErrTestFails");
+        utils::check_cvode_flag(CVodeGetNumErrTestFails(m_cvode_mem, &netfails), "CVodeGetNumErrTestFails");
-        check_cvode_flag(CVodeGetNumNonlinSolvIters(m_cvode_mem, &nniters), "CVodeGetNumNonlinSolvIters");
+        utils::check_cvode_flag(CVodeGetNumNonlinSolvIters(m_cvode_mem, &nniters), "CVodeGetNumNonlinSolvIters");
-        check_cvode_flag(CVodeGetNumNonlinSolvConvFails(m_cvode_mem, &nconvfails), "CVodeGetNumNonlinSolvConvFails");
+        utils::check_cvode_flag(CVodeGetNumNonlinSolvConvFails(m_cvode_mem, &nconvfails), "CVodeGetNumNonlinSolvConvFails");
-        check_cvode_flag(CVodeGetNumLinConvFails(m_cvode_mem, &nsetfails), "CVodeGetNumLinConvFails");
+        utils::check_cvode_flag(CVodeGetNumLinConvFails(m_cvode_mem, &nsetfails), "CVodeGetNumLinConvFails");
 
 
         {
@@ -864,7 +807,7 @@ namespace gridfire::solver {
         }
 
         // --- 3. Get Estimated Local Errors (Your Original Logic) ---
-        check_cvode_flag(CVodeGetEstLocalErrors(m_cvode_mem, m_YErr), "CVodeGetEstLocalErrors");
+        utils::check_cvode_flag(CVodeGetEstLocalErrors(m_cvode_mem, m_YErr), "CVodeGetEstLocalErrors");
 
         sunrealtype *y_data = N_VGetArrayPointer(m_Y);
         sunrealtype *y_err_data = N_VGetArrayPointer(m_YErr);

subprojects/cvode.wrap

@@ -1,6 +1,6 @@
 [wrap-file]
-directory = cvode-7.3.0
+directory = cvode-7.5.0
 
-source_url = https://github.com/LLNL/sundials/releases/download/v7.3.0/cvode-7.3.0.tar.gz
+source_url = https://github.com/LLNL/sundials/releases/download/v7.5.0/cvode-7.5.0.tar.gz
-source_filename = cvode-7.3.0.tar.gz
+source_filename = cvode-7.5.0.tar.gz
-source_hash = 8b15a646882f2414b1915cad4d53136717a077539e7cfc480f2002c5898ae568
+source_hash = e588d1c569c2fe92403c6b1a8bcce973624bf6ead4a6a302ab7cc3c4f39523bd

subprojects/kinsol.wrap

@@ -0,0 +1,6 @@
+[wrap-file]
+directory = kinsol-7.5.0
+
+source_url = https://github.com/LLNL/sundials/releases/download/v7.5.0/kinsol-7.5.0.tar.gz
+source_filename = kinsol-7.5.0.tar.gz
+source_hash = ce712ea29aaec537a6cd234d1a6bea024592ef54a26a5d457dda2320f0c49d07

tests/graphnet_sandbox/main.cpp

@@ -66,7 +66,7 @@ gridfire::NetIn init(const double temp) {
     netIn.density = 1.5e2;
     netIn.energy = 0;
 
-    netIn.tMax = 3e17;
+    netIn.tMax = 1e17;
     netIn.dt0 = 1e-12;
 
     return netIn;