second draft

include/micm/configure/solver_config.hpp

@@ -426,7 +426,7 @@ namespace micm
     void ParseRelativeTolerance(const objectType& object)
     {
       ValidateSchema(object, { "value", "type" }, {});
-      this->parameters_.relative_tolerance_ = object["value"].as<double>();
+      //this->parameters_.relative_tolerance_ = object["value"].as<double>();
     }
 
     void ParseMechanism(const objectType& object)

include/micm/solver/backward_euler.hpp

@@ -80,12 +80,12 @@ namespace micm
     static bool IsConverged(
         const BackwardEulerSolverParameters& parameters,
         const DenseMatrixPolicy& residual,
-        const DenseMatrixPolicy& state) requires(!VectorizableDense<DenseMatrixPolicy>);
+        const DenseMatrixPolicy& state, auto& stateParams) requires(!VectorizableDense<DenseMatrixPolicy>);
     template<class DenseMatrixPolicy>
     static bool IsConverged(
         const BackwardEulerSolverParameters& parameters,
         const DenseMatrixPolicy& residual,
-        const DenseMatrixPolicy& state) requires(VectorizableDense<DenseMatrixPolicy>);
+        const DenseMatrixPolicy& state, auto& stateParams) requires(VectorizableDense<DenseMatrixPolicy>);
   };
 
 }  // namespace micm

include/micm/solver/backward_euler.inl

@@ -135,7 +135,7 @@ namespace micm
           continue;
 
         // check for convergence
-        converged = IsConverged(parameters_, forcing, Yn1);
+        converged = IsConverged(parameters_, forcing, Yn1, state);
       } while (!converged && iterations < max_iter);
 
       if (!converged)
@@ -182,11 +182,11 @@ namespace micm
   inline bool BackwardEuler<RatesPolicy, LinearSolverPolicy>::IsConverged(
       const BackwardEulerSolverParameters& parameters,
       const DenseMatrixPolicy& residual,
-      const DenseMatrixPolicy& state) requires(!VectorizableDense<DenseMatrixPolicy>)
+      const DenseMatrixPolicy& state, auto& stateParams) requires(!VectorizableDense<DenseMatrixPolicy>)
   {
     double small = parameters.small_;
-    double rel_tol = parameters.relative_tolerance_;
-    auto& abs_tol = parameters.absolute_tolerance_;
+    double rel_tol = stateParams.relative_tolerance_;
+    auto& abs_tol = stateParams.absolute_tolerance_;
     auto residual_iter = residual.AsVector().begin();
     auto state_iter = state.AsVector().begin();
     const std::size_t n_elem = residual.NumRows() * residual.NumColumns();
@@ -208,11 +208,11 @@ namespace micm
   inline bool BackwardEuler<RatesPolicy, LinearSolverPolicy>::IsConverged(
       const BackwardEulerSolverParameters& parameters,
       const DenseMatrixPolicy& residual,
-      const DenseMatrixPolicy& state) requires(VectorizableDense<DenseMatrixPolicy>)
+      const DenseMatrixPolicy& state, auto& stateParams) requires(VectorizableDense<DenseMatrixPolicy>)
   {
     double small = parameters.small_;
-    double rel_tol = parameters.relative_tolerance_;
-    auto& abs_tol = parameters.absolute_tolerance_;
+    double rel_tol = stateParams.relative_tolerance_;
+    auto& abs_tol = stateParams.absolute_tolerance_;
     auto residual_iter = residual.AsVector().begin();
     auto state_iter = state.AsVector().begin();
     const std::size_t n_elem = residual.NumRows() * residual.NumColumns();

include/micm/solver/backward_euler_solver_parameters.hpp

@@ -19,8 +19,7 @@ namespace micm
     template<class RatesPolicy, class LinearSolverPolicy>
     using SolverType = BackwardEuler<RatesPolicy, LinearSolverPolicy>;
 
-    std::vector<double> absolute_tolerance_;
-    double relative_tolerance_{ 1.0e-8 };
+    //double relative_tolerance_{ 1.0e-8 };
     double small_{ 1.0e-40 };
     size_t max_number_of_steps_{ 11 };
     // The time step reductions are used to determine the time step after a failed solve

include/micm/solver/rosenbrock.hpp

@@ -123,10 +123,10 @@ namespace micm
     /// @param errors The computed errors
     /// @return
     template<class DenseMatrixPolicy>
-    double NormalizedError(const DenseMatrixPolicy& y, const DenseMatrixPolicy& y_new, const DenseMatrixPolicy& errors) const
+    double NormalizedError(const DenseMatrixPolicy& y, const DenseMatrixPolicy& y_new, const DenseMatrixPolicy& errors, auto& state) const
         requires(!VectorizableDense<DenseMatrixPolicy>);
     template<class DenseMatrixPolicy>
-    double NormalizedError(const DenseMatrixPolicy& y, const DenseMatrixPolicy& y_new, const DenseMatrixPolicy& errors) const
+    double NormalizedError(const DenseMatrixPolicy& y, const DenseMatrixPolicy& y_new, const DenseMatrixPolicy& errors, auto& state) const
         requires(VectorizableDense<DenseMatrixPolicy>);
   };  // end of Abstract Rosenbrock Solver
 

include/micm/solver/rosenbrock.inl

@@ -122,7 +122,7 @@ namespace micm
           Yerror.Axpy(parameters_.e_[stage], K[stage]);
 
         // Compute the normalized error
-        auto error = static_cast<const Derived*>(this)->NormalizedError(Y, Ynew, Yerror);
+        auto error = static_cast<const Derived*>(this)->NormalizedError(Y, Ynew, Yerror, state);
 
         // New step size is bounded by FacMin <= Hnew/H <= FacMax
         double fac = std::min(
@@ -267,7 +267,8 @@ namespace micm
   inline double AbstractRosenbrockSolver<RatesPolicy, LinearSolverPolicy, Derived>::NormalizedError(
       const DenseMatrixPolicy& Y,
       const DenseMatrixPolicy& Ynew,
-      const DenseMatrixPolicy& errors) const requires(!VectorizableDense<DenseMatrixPolicy>)
+      const DenseMatrixPolicy& errors, 
+      auto& state) const requires(!VectorizableDense<DenseMatrixPolicy>)
   {
     // Solving Ordinary Differential Equations II, page 123
     // https://link-springer-com.cuucar.idm.oclc.org/book/10.1007/978-3-642-05221-7
@@ -278,7 +279,7 @@ namespace micm
     auto& _ynew = Ynew.AsVector();
     auto& _errors = errors.AsVector();
     const std::size_t N = Y.AsVector().size();
-    const std::size_t n_vars = parameters_.absolute_tolerance_.size();
+    const std::size_t n_vars = state.absolute_tolerance_.size();
 
     double ymax = 0;
     double errors_over_scale = 0;
@@ -288,7 +289,7 @@ namespace micm
     {
       ymax = std::max(std::abs(_y[i]), std::abs(_ynew[i]));
       errors_over_scale =
-          _errors[i] / (parameters_.absolute_tolerance_[i % n_vars] + parameters_.relative_tolerance_ * ymax);
+          _errors[i] / (state.absolute_tolerance_[i % n_vars] + state.relative_tolerance_ * ymax);
       error += errors_over_scale * errors_over_scale;
     }
 
@@ -302,7 +303,8 @@ namespace micm
   inline double AbstractRosenbrockSolver<RatesPolicy, LinearSolverPolicy, Derived>::NormalizedError(
       const DenseMatrixPolicy& Y,
       const DenseMatrixPolicy& Ynew,
-      const DenseMatrixPolicy& errors) const requires(VectorizableDense<DenseMatrixPolicy>)
+      const DenseMatrixPolicy& errors, 
+      auto& state) const requires(VectorizableDense<DenseMatrixPolicy>)
   {
     // Solving Ordinary Differential Equations II, page 123
     // https://link-springer-com.cuucar.idm.oclc.org/book/10.1007/978-3-642-05221-7
@@ -314,7 +316,7 @@ namespace micm
     auto errors_iter = errors.AsVector().begin();
     const std::size_t N = Y.NumRows() * Y.NumColumns();
     const std::size_t L = Y.GroupVectorSize();
-    const std::size_t n_vars = parameters_.absolute_tolerance_.size();
+    const std::size_t n_vars = state.absolute_tolerance_.size();
 
     const std::size_t whole_blocks = std::floor(Y.NumRows() / Y.GroupVectorSize()) * Y.GroupSize();
 
@@ -325,8 +327,8 @@ namespace micm
     for (std::size_t i = 0; i < whole_blocks; ++i)
     {
       errors_over_scale =
-          *errors_iter / (parameters_.absolute_tolerance_[(i / L) % n_vars] +
-                          parameters_.relative_tolerance_ * std::max(std::abs(*y_iter), std::abs(*ynew_iter)));
+          *errors_iter / (state.absolute_tolerance_[(i / L) % n_vars] +
+                          state.relative_tolerance_ * std::max(std::abs(*y_iter), std::abs(*ynew_iter)));
       error += errors_over_scale * errors_over_scale;
       ++y_iter;
       ++ynew_iter;
@@ -344,8 +346,8 @@ namespace micm
         {
           const std::size_t idx = y * L + x;
           errors_over_scale = errors_iter[idx] /
-                              (parameters_.absolute_tolerance_[y] +
-                               parameters_.relative_tolerance_ * std::max(std::abs(y_iter[idx]), std::abs(ynew_iter[idx])));
+                              (state.absolute_tolerance_[y] +
+                               state.relative_tolerance_ * std::max(std::abs(y_iter[idx]), std::abs(ynew_iter[idx])));
           error += errors_over_scale * errors_over_scale;
         }
       }

include/micm/solver/rosenbrock_solver_parameters.hpp

@@ -9,6 +9,8 @@
 #include <limits>
 #include <vector>
 
+//#include "../include/micm/solver/solver.hpp"
+
 namespace micm
 {
 
@@ -60,9 +62,6 @@ namespace micm
     // Gamma_i = \sum_j  gamma_{i,j}
     std::array<double, 6> gamma_{};
 
-    std::vector<double> absolute_tolerance_;
-    double relative_tolerance_{ 1e-6 };
-
     bool check_singularity_{ false };  // Check for singular A matrix in linear solve of A x = b
 
     // Print RosenbrockSolverParameters to console
@@ -132,12 +131,6 @@ namespace micm
       std::cout << val << " ";
     std::cout << std::endl;
     std::cout << "absolute_tolerance: ";
-    for (auto& val : absolute_tolerance_)
-    {
-      std::cout << val << " ";
-    }
-    std::cout << std::endl;
-    std::cout << "relative_tolerance: " << relative_tolerance_ << std::endl;
   }
 
   inline RosenbrockSolverParameters RosenbrockSolverParameters::TwoStageRosenbrockParameters()

include/micm/solver/solver.hpp

@@ -64,9 +64,9 @@ namespace micm
     }
 
     // Overloaded Solve function to change parameters
-    SolverResult Solve(double time_step, StatePolicy& state, RosenbrockSolverParameters params)
+    SolverResult Solve(double time_step, StatePolicy& state, SolverParametersType& params)
     {
-      solver_.parameters_.h_start_ = params.h_start_;
+      solver_.parameters_ = params;
       return solver_.Solve(time_step, state);
     }
 
@@ -88,7 +88,7 @@ namespace micm
     {
       return state_parameters_.number_of_rate_constants_;
     }
-
+    
     StatePolicy GetState() const
     {
       auto state = std::move(StatePolicy(state_parameters_));

include/micm/solver/solver_builder.inl

@@ -379,8 +379,7 @@ namespace micm
     }
 
     this->UnusedSpeciesCheck();
-    this->SetAbsoluteTolerances(options.absolute_tolerance_, species_map);
-
+
     RatesPolicy rates(this->reactions_, species_map);
     auto nonzero_elements = rates.NonZeroJacobianElements();
     auto jacobian = BuildJacobian<SparseMatrixPolicy>(nonzero_elements, this->number_of_grid_cells_, number_of_species);
@@ -398,6 +397,9 @@ namespace micm
                                          .variable_names_ = variable_names,
                                          .custom_rate_parameter_labels_ = labels,
                                          .nonzero_jacobian_elements_ = nonzero_elements };
+
+
+    this->SetAbsoluteTolerances(state_parameters.absolute_tolerance_, species_map); //WHERE TO ADD THIS??????
 
     return Solver<SolverPolicy, StatePolicy>(
         SolverPolicy(options, std::move(linear_solver), std::move(rates), jacobian),

include/micm/solver/state.hpp

@@ -31,6 +31,8 @@ namespace micm
     std::vector<std::string> variable_names_{};
     std::vector<std::string> custom_rate_parameter_labels_{};
     std::set<std::pair<std::size_t, std::size_t>> nonzero_jacobian_elements_{};
+    double relative_tolerance_;
+    std::vector<double> absolute_tolerance_ {};
   };
 
   template<class DenseMatrixPolicy = StandardDenseMatrix, class SparseMatrixPolicy = StandardSparseMatrix>
@@ -58,6 +60,8 @@ namespace micm
     std::size_t state_size_;
     std::size_t number_of_grid_cells_;
     std::unique_ptr<TemporaryVariables> temporary_variables_;
+    double relative_tolerance_;
+    std::vector<double> absolute_tolerance_;
 
     /// @brief Copy constructor
     /// @param other The state object to be copied
@@ -76,6 +80,8 @@ namespace micm
       state_size_ = other.state_size_;
       number_of_grid_cells_ = other.number_of_grid_cells_;
       temporary_variables_ = std::make_unique<TemporaryVariables>(*other.temporary_variables_);
+      relative_tolerance_ = other.relative_tolerance_;
+      absolute_tolerance_ = other.absolute_tolerance_;
     }
 
     /// @brief Assignment operator
@@ -98,6 +104,8 @@ namespace micm
         state_size_ = other.state_size_;
         number_of_grid_cells_ = other.number_of_grid_cells_;
         temporary_variables_ = std::make_unique<TemporaryVariables>(*other.temporary_variables_);
+        relative_tolerance_ = other.relative_tolerance_;
+        absolute_tolerance_ = other.absolute_tolerance_;
       }
       return *this;
     }

include/micm/solver/state.inl

@@ -80,7 +80,9 @@ namespace micm
         lower_matrix_(),
         upper_matrix_(),
         state_size_(parameters.variable_names_.size()),
-        number_of_grid_cells_(parameters.number_of_grid_cells_)
+        number_of_grid_cells_(parameters.number_of_grid_cells_),
+        relative_tolerance_(1e-06),
+        absolute_tolerance_(parameters.absolute_tolerance_)
   {
     std::size_t index = 0;
     for (auto& name : variable_names_)

test/integration/analytical_policy.hpp

@@ -1375,6 +1375,8 @@ void test_analytical_robertson(
   double air_density = 1e6;
 
   auto state = solver.GetState();
+  state.relative_tolerance_ = 1e-10;
+  state.absolute_tolerance_ = std::vector<double>(3, state.relative_tolerance_ * 1e-2);
 
   double k1 = 0.04;
   double k2 = 3e7;
@@ -1583,6 +1585,9 @@ void test_analytical_oregonator(
 
   auto state = solver.GetState();
 
+  state.relative_tolerance_ = 1e-6;
+  state.absolute_tolerance_ = std::vector<double>(5, state.relative_tolerance_ * 1e-2);
+
   state.SetCustomRateParameter("r1", 1.34 * 0.06);
   state.SetCustomRateParameter("r2", 1.6e9);
   state.SetCustomRateParameter("r3", 8e3 * 0.06);
@@ -1758,6 +1763,8 @@ void test_analytical_hires(
   };
 
   auto state = solver.GetState();
+  state.relative_tolerance_ = 1e-6;
+  state.absolute_tolerance_ = std::vector<double>(8, state.relative_tolerance_ * 1e-2);
 
   state.SetCustomRateParameter("r1", 1.71);
   state.SetCustomRateParameter("r2", 8.75);
@@ -1916,6 +1923,12 @@ void test_analytical_e5(
 
   auto state = solver.GetState();
 
+  state.relative_tolerance_ = 1e-13;
+  state.absolute_tolerance_ = std::vector<double>(6, 1e-17);
+  state.absolute_tolerance_[0] = 1e-7;
+  state.absolute_tolerance_[4] = 1e-7;
+  state.absolute_tolerance_[5] = 1e-7;
+
   state.SetCustomRateParameter("r1", 7.89e-10);
   state.SetCustomRateParameter("r2", 1.13e9);
   state.SetCustomRateParameter("r3", 1.1e7);

test/integration/cuda/test_cuda_analytical_rosenbrock.cpp

@@ -180,14 +180,14 @@ TEST(AnalyticalExamplesCudaRosenbrock, E5)
 {
   auto rosenbrock_solver = [](auto params)
   {
-    params.relative_tolerance_ = 1e-13;
+    /*params.relative_tolerance_ = 1e-13;
     params.absolute_tolerance_ = std::vector<double>(6, 1e-17);
     // this paper https://archimede.uniba.it/~testset/report/e5.pdf
     // says that the first variable should have a much looser tolerance than the other species
     params.absolute_tolerance_[0] = 1e-7;
     // these last two aren't actually provided values and we don't care how they behave
     params.absolute_tolerance_[4] = 1e-7;
-    params.absolute_tolerance_[5] = 1e-7;
+    params.absolute_tolerance_[5] = 1e-7;*/
     return builderType1Cell(params);
   };
 
@@ -212,8 +212,8 @@ TEST(AnalyticalExamplesCudaRosenbrock, Oregonator)
   auto rosenbrock_solver = [](auto params)
   {
     // anything below 1e-6 is too strict for the Oregonator
-    params.relative_tolerance_ = 1e-6;
-    params.absolute_tolerance_ = std::vector<double>(5, params.relative_tolerance_ * 1e-2);
+    //params.relative_tolerance_ = 1e-6;
+    //params.absolute_tolerance_ = std::vector<double>(5, params.relative_tolerance_ * 1e-2);
     return builderType1Cell(params);
   };
 
@@ -237,8 +237,8 @@ TEST(AnalyticalExamplesCudaRosenbrock, HIRES)
 {
   auto rosenbrock_solver = [](auto params)
   {
-    params.relative_tolerance_ = 1e-6;
-    params.absolute_tolerance_ = std::vector<double>(8, params.relative_tolerance_ * 1e-2);
+    //params.relative_tolerance_ = 1e-6;
+    //params.absolute_tolerance_ = std::vector<double>(8, params.relative_tolerance_ * 1e-2);
     return builderType1Cell(params);
   };