fix(cache-+-tolerance): added cache clearning and dynamic runtime tolerances

Previously caches were unbound and this leand to linear memrory growth with the number of steps taken. Further, the numeric tolerance was distributed equally amoung all species. This lead to cases where trace species could dominate the error budget leading to far more steps than truley needed. Now the RHS cache is cleared after ever accepted CVODE step and numeric tolerances are distributed based on the relative abundance of species. This has lead to a dramatic speed up (~4x) in certain cases while also preventing catastrophic memory growth

src/include/gridfire/solver/strategies/PointSolver.h

@@ -263,9 +263,16 @@ namespace gridfire::solver {
             PointSolverContext* sctx;    // Pointer to the solver context
             engine::scratch::StateBlob& ctx;
             const engine::DynamicEngine* engine{};
+            double abs_tol{};
+            double rel_tol{};
             double T9{};
             double rho{};
-            double energy{};
+            double energy_generation_rate{};
+            double last_step_size{};
+            double previous_accepted_energy{};
+            double last_energy_change{};
+            std::vector<double> previous_accepted_abundances;
+            std::vector<double> last_mass_fraction_changes;
             const std::vector<fourdst::atomic::Species>* networkSpecies{};
             std::unique_ptr<exceptions::EngineError> captured_exception = nullptr;
             std::optional<std::map<fourdst::atomic::Species, std::unordered_map<std::string, double>>> reaction_contribution_map;
@@ -287,6 +294,11 @@ namespace gridfire::solver {
          * @return 0 on success; 1 on recoverable StaleEngineTrigger; -1 on other failures.
          */
         static int cvode_rhs_wrapper(sunrealtype t, N_Vector y, N_Vector ydot, void *user_data);
+
+        /**
+         * @brief Builds CVODE error weights from runtime species and energy scales.
+         */
+        static int cvode_error_weight_wrapper(N_Vector y, N_Vector ewt, void* user_data);
         /**
          * @brief CVODE dense Jacobian C-wrapper that fills SUNDenseMatrix using the engine.
          *
@@ -320,7 +332,7 @@ namespace gridfire::solver {
          *
          * State vector m_Y is sized to N (numSpecies + 1). Species slots are initialized from Composition
          * molar abundances when present, otherwise a tiny positive value; the last slot is set to
-         * accumulatedEnergy. Sets scalar tolerances, non-negativity constraints for species, maximum
+         * accumulatedEnergy. Sets runtime-derived error weights, non-negativity constraints for species, maximum
          * step size, creates a dense matrix and dense linear solver, and registers the Jacobian.
          */
         void initialize_cvode_integration_resources(

src/lib/solver/strategies/PointSolver.cpp

@@ -3,6 +3,7 @@
 #include "gridfire/types/types.h"
 #include "gridfire/utils/table_format.h"
 #include "gridfire/engine/diagnostics/dynamic_engine_diagnostics.h"
+#include "gridfire/engine/scratchpads/engine_graph_scratchpad.h"
 
 #include "quill/LogMacros.h"
 
@@ -10,6 +11,7 @@
 
 // ReSharper disable once CppUnusedIncludeDirective
 #include <cstdint>
+#include <cmath>
 #include <limits>
 #include <string>
 #include <unordered_map>
@@ -291,7 +293,6 @@ namespace gridfire::solver {
                 }
             }
             y_data[numSpecies] = 0.0; // Reset energy accumulator
-            utils::check_cvode_flag(CVodeSStolerances(sctx_p->cvode_mem, sctx_p->rel_tol.value(), sctx_p->abs_tol.value()), "CVodeSStolerances");
             utils::check_cvode_flag(CVodeReInit(sctx_p->cvode_mem, 0.0, sctx_p->Y), "CVodeReInit");
 
             equilibratedComposition = netIn.composition; // Use the provided composition as-is if we already have validated CVODE resources and that the composition is consistent with the previous state
@@ -303,7 +304,19 @@ namespace gridfire::solver {
             .sctx = sctx_p,
             .ctx = *sctx_p->engine_ctx,
             .engine = &m_engine,
+            .abs_tol = sctx_p->abs_tol.value(),
+            .rel_tol = sctx_p->rel_tol.value(),
         };
+        auto reset_error_weight_history = [&]() {
+            const auto* y_data = N_VGetArrayPointer(sctx_p->Y);
+            user_data.previous_accepted_abundances.assign(y_data, y_data + numSpecies);
+            user_data.last_mass_fraction_changes.assign(numSpecies, 0.0);
+            user_data.previous_accepted_energy = y_data[numSpecies];
+            user_data.last_energy_change = 0.0;
+            user_data.last_step_size = 0.0;
+            user_data.energy_generation_rate = 0.0;
+        };
+        reset_error_weight_history();
         LOG_TRACE_L1(m_logger, "CVODE resources successfully initialized!");
 
         double current_time = 0;
@@ -345,6 +358,10 @@ namespace gridfire::solver {
 
             utils::check_cvode_flag(flag, "CVode");
 
+            if (auto graph_state = sctx_p->engine_ctx->get<scratch::GraphEngineScratchPad>(); graph_state) {
+                graph_state.value()->stepDerivativesCache.clear();
+            }
+
             long int n_steps;
             double last_step_size;
             CVodeGetNumSteps(sctx_p->cvode_mem, &n_steps);
@@ -356,6 +373,19 @@ namespace gridfire::solver {
             sunrealtype* y_data = N_VGetArrayPointer(sctx_p->Y);
             const double current_energy = y_data[numSpecies]; // Specific energy rate
 
+            if (user_data.previous_accepted_abundances.size() == numSpecies) {
+                for (size_t i = 0; i < numSpecies; ++i) {
+                    const double mass_number = static_cast<double>(user_data.networkSpecies->at(i).a());
+                    user_data.last_mass_fraction_changes[i] = mass_number * std::abs(
+                        y_data[i] - user_data.previous_accepted_abundances[i]
+                    );
+                    user_data.previous_accepted_abundances[i] = y_data[i];
+                }
+            }
+            user_data.last_energy_change = std::abs(current_energy - user_data.previous_accepted_energy);
+            user_data.previous_accepted_energy = current_energy;
+            user_data.last_step_size = std::abs(last_step_size);
+
             // TODO: Accumulate neutrino loss through the state vector directly which will allow CVODE to properly integrate it
             accumulated_neutrino_energy_loss += user_data.neutrino_energy_loss_rate * last_step_size;
             accumulated_total_neutrino_flux += user_data.total_neutrino_flux * last_step_size;
@@ -607,6 +637,8 @@ namespace gridfire::solver {
                 initialize_cvode_integration_resources(sctx_p, N, numSpecies, current_time, currentComposition, sctx_p->abs_tol.value(), sctx_p->rel_tol.value(), accumulated_energy);
 
                 utils::check_cvode_flag(CVodeReInit(sctx_p->cvode_mem, current_time, sctx_p->Y), "CVodeReInit");
+                user_data.networkSpecies = &m_engine.getNetworkSpecies(*sctx_p->engine_ctx);
+                reset_error_weight_history();
                 LOG_INFO(m_logger, "Done reinitializing CVODE after engine update. The next log messages will be from the first step after reinitialization...");
             }
 
@@ -706,6 +738,7 @@ namespace gridfire::solver {
             data->reaction_contribution_map = result.reaction_contribution_map;
             data->neutrino_energy_loss_rate = result.neutrino_energy_loss_rate;
             data->total_neutrino_flux = result.total_neutrino_flux;
+            data->energy_generation_rate = N_VGetArrayPointer(ydot)[data->networkSpecies->size()];
             LOG_TRACE_L2(instance->m_logger, "CVODE RHS wrapper completed successfully at time {}", t);
             return 0;
         } catch (const exceptions::EngineError& e) {
@@ -718,6 +751,81 @@ namespace gridfire::solver {
         }
     }
 
+    int PointSolver::cvode_error_weight_wrapper(
+        const N_Vector y,
+        const N_Vector ewt,
+        void* user_data
+    ) {
+        const auto* data = static_cast<CVODEUserData*>(user_data);
+        if (data == nullptr || data->networkSpecies == nullptr) {
+            return -1;
+        }
+        if (!std::isfinite(data->abs_tol) || !std::isfinite(data->rel_tol) ||
+            data->abs_tol < 0.0 || data->rel_tol < 0.0) {
+            return -1;
+        }
+
+        const size_t num_species = data->networkSpecies->size();
+        const auto expected_length = static_cast<sunindextype>(num_species + 1);
+        if (N_VGetLength(y) != expected_length || N_VGetLength(ewt) != expected_length) {
+            return -1;
+        }
+
+        const auto* y_data = N_VGetArrayPointer(y);
+        auto* weight_data = N_VGetArrayPointer(ewt);
+
+        double total_mass_fraction = 0.0;
+        for (size_t i = 0; i < num_species; ++i) {
+            const double mass_number = static_cast<double>(data->networkSpecies->at(i).a());
+            total_mass_fraction += mass_number * std::abs(y_data[i]);
+        }
+        if (!std::isfinite(total_mass_fraction) || total_mass_fraction <= 0.0) {
+            return -1;
+        }
+
+        const double numerical_floor = std::numeric_limits<double>::epsilon() * total_mass_fraction;
+        const bool have_change_history = data->last_mass_fraction_changes.size() == num_species;
+        const auto importance = [&](const size_t i) {
+            const double mass_number = static_cast<double>(data->networkSpecies->at(i).a());
+            const double current_mass_fraction = mass_number * std::abs(y_data[i]);
+            const double recent_change = have_change_history ? data->last_mass_fraction_changes[i] : 0.0;
+            return std::max({current_mass_fraction, recent_change, numerical_floor});
+        };
+
+        double total_importance = 0.0;
+        for (size_t i = 0; i < num_species; ++i) {
+            total_importance += importance(i);
+        }
+        if (!std::isfinite(total_importance) || total_importance <= 0.0) {
+            return -1;
+        }
+
+        for (size_t i = 0; i < num_species; ++i) {
+            const double mass_number = static_cast<double>(data->networkSpecies->at(i).a());
+            const double component_importance = importance(i);
+            const double mass_fraction_tolerance =
+                data->rel_tol * component_importance +
+                data->abs_tol * component_importance / total_importance;
+            const double molar_abundance_tolerance = mass_fraction_tolerance / mass_number;
+            if (!std::isfinite(molar_abundance_tolerance) || molar_abundance_tolerance <= 0.0) {
+                return -1;
+            }
+            weight_data[i] = 1.0 / molar_abundance_tolerance;
+        }
+
+        const double energy_scale = std::max({
+            std::abs(y_data[num_species]),
+            data->last_energy_change,
+            std::abs(data->last_step_size * data->energy_generation_rate)
+        });
+        const double energy_tolerance = data->rel_tol * energy_scale + data->abs_tol;
+        if (!std::isfinite(energy_tolerance) || energy_tolerance <= 0.0) {
+            return -1;
+        }
+        weight_data[num_species] = 1.0 / energy_tolerance;
+        return 0;
+    }
+
     int PointSolver::cvode_jac_wrapper(
         sunrealtype t,
         N_Vector y,
@@ -942,7 +1050,7 @@ namespace gridfire::solver {
 
 
         utils::check_cvode_flag(CVodeInit(sctx_p->cvode_mem, cvode_rhs_wrapper, current_time, sctx_p->Y), "CVodeInit");
-        utils::check_cvode_flag(CVodeSStolerances(sctx_p->cvode_mem, relTol, absTol), "CVodeSStolerances");
+        utils::check_cvode_flag(CVodeWFtolerances(sctx_p->cvode_mem, cvode_error_weight_wrapper), "CVodeWFtolerances");
         utils::check_cvode_flag(CVodeSetInitStep(sctx_p->cvode_mem, 1.0e-8), "CVodeSetInitStep");
 
         // Constraints