Remove solver_.parameters_

examples/profile_example.cpp

@@ -65,14 +65,15 @@ int Run(const char* filepath, const char* initial_conditions, const std::string&
   auto reactions = solver_params.processes_;
 
   auto params = RosenbrockSolverParameters::ThreeStageRosenbrockParameters();
-  params.relative_tolerance_ = 0.1;
+
 
   auto solver = VectorBuilder<L>(params)
                     .SetSystem(chemical_system)
                     .SetReactions(reactions)
                     .Build();
   auto state = solver.GetState();
-
+  state.SetRelativeTolerances(0.1);
+
   state.conditions_[0].temperature_ = dataMap.environments["temperature"];  // K
   state.conditions_[0].pressure_ = dataMap.environments["pressure"];        // Pa
   state.conditions_[0].air_density_ = dataMap.environments["air_density"];  // mol m-3

include/micm/configure/solver_config.hpp

@@ -102,6 +102,7 @@ namespace micm
     System system_;
     std::vector<Process> processes_;
     RosenbrockSolverParameters parameters_;
+    double relative_tolerance_;    
 
     SolverParameters(const System& system, std::vector<Process>&& processes, const RosenbrockSolverParameters&& parameters)
         : system_(system),
@@ -116,6 +117,13 @@ namespace micm
           parameters_(parameters)
     {
     }
+    SolverParameters(System&& system, std::vector<Process>&& processes, RosenbrockSolverParameters&& parameters, double relative_tolerance)
+        : system_(system),
+          processes_(processes),
+          parameters_(parameters), 
+          relative_tolerance_(relative_tolerance)
+    {
+    }
   };
 
   class ConfigReaderPolicy
@@ -138,6 +146,7 @@ namespace micm
     std::unordered_map<std::string, Phase> phases_;
     std::vector<Process> processes_;
     RosenbrockSolverParameters parameters_;
+    double relative_tolerance_;
 
     // Common YAML
     inline static const std::string DEFAULT_CONFIG_FILE_JSON = "config.json";
@@ -273,7 +282,7 @@ namespace micm
       processes_.clear();
 
       // Parse species object array
-      ParseSpeciesArray(species_objects);
+      ParseSpeciesArray(species_objects);     
 
       // Assign the parsed 'Species' to 'Phase'
       std::vector<Species> species_arr;
@@ -426,7 +435,7 @@ namespace micm
     void ParseRelativeTolerance(const objectType& object)
     {
       ValidateSchema(object, { "value", "type" }, {});
-      this->parameters_.relative_tolerance_ = object["value"].as<double>();
+      this->relative_tolerance_ = object["value"].as<double>();
     }
 
     void ParseMechanism(const objectType& object)
@@ -959,7 +968,7 @@ namespace micm
     SolverParameters GetSolverParams()
     {
       return SolverParameters(
-          std::move(System(this->gas_phase_, this->phases_)), std::move(this->processes_), std::move(this->parameters_));
+          std::move(System(this->gas_phase_, this->phases_)), std::move(this->processes_), std::move(this->parameters_), std::move(this->relative_tolerance_));
     }
   };
 }  // namespace micm

include/micm/cuda/solver/cuda_rosenbrock.hpp

@@ -115,7 +115,7 @@ namespace micm
     /// @param errors The computed errors
     /// @return The scaled norm of the errors
     template<class DenseMatrixPolicy>
-    double NormalizedError(const DenseMatrixPolicy& y_old, const DenseMatrixPolicy& y_new, const DenseMatrixPolicy& errors)
+    double NormalizedError(const DenseMatrixPolicy& y_old, const DenseMatrixPolicy& y_new, const DenseMatrixPolicy& errors, auto& state)
         const requires(CudaMatrix<DenseMatrixPolicy>&& VectorizableDense<DenseMatrixPolicy>)
     {
       return micm::cuda::NormalizedErrorDriver(

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

@@ -133,7 +133,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)
@@ -180,11 +180,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();
@@ -206,18 +206,17 @@ 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();
     const std::size_t L = residual.GroupVectorSize();
     const std::size_t n_vars = abs_tol.size();
     const std::size_t whole_blocks = std::floor(residual.NumRows() / L) * residual.GroupSize();
-
     // evaluate the rows that fit exactly into the vectorizable dimension (L)
     for (std::size_t i = 0; i < whole_blocks; ++i)
     {

include/micm/solver/rosenbrock.hpp

@@ -116,10 +116,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

@@ -116,7 +116,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(
@@ -242,7 +242,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
@@ -253,7 +254,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;
@@ -263,7 +264,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;
     }
 
@@ -277,7 +278,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
@@ -289,7 +291,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();
 
@@ -300,8 +302,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;
@@ -319,8 +321,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

@@ -60,6 +60,7 @@ namespace micm
     // Gamma_i = \sum_j  gamma_{i,j}
     std::array<double, 6> gamma_{};
 
+    //We have this to pass CUDA Tests. once we figure out a path forward we can get rid of these. 
     std::vector<double> absolute_tolerance_;
     double relative_tolerance_{ 1e-6 };
 
@@ -130,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

@@ -63,6 +63,13 @@ namespace micm
       return solver_.Solve(time_step, state);
     }
 
+    // Overloaded Solve function to change parameters
+    SolverResult Solve(double time_step, StatePolicy& state, SolverParametersType& params)
+    {
+      solver_.parameters_ = params;
+      return solver_.Solve(time_step, state);
+    }
+
     /// @brief Returns the number of grid cells
     /// @return
     std::size_t GetNumberOfGridCells() const
@@ -81,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.hpp

@@ -50,6 +50,7 @@ namespace micm
     bool reorder_state_ = true;
     bool valid_system_ = false;
     bool valid_reactions_ = false;
+    double relative_tolerance_ = 1e-06;
 
    public:
     SolverBuilder() = delete;

include/micm/solver/solver_builder.inl

@@ -316,7 +316,7 @@ namespace micm
         }
       }
     }
-  }
+  }  
 
   template<
       class SolverParametersPolicy,
@@ -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,10 @@ namespace micm
                                          .variable_names_ = variable_names,
                                          .custom_rate_parameter_labels_ = labels,
                                          .nonzero_jacobian_elements_ = nonzero_elements };
+
+
+    this->SetAbsoluteTolerances(state_parameters.absolute_tolerance_, species_map);
+    options.absolute_tolerance_ = state_parameters.absolute_tolerance_;
 
     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_{ 1e-06 };
+    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 Default constructor
     /// Only defined to be used to create default values in types, but a default constructed state is not useable
@@ -84,6 +88,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
@@ -106,6 +112,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;
     }
@@ -179,6 +187,7 @@ namespace micm
     /// @param value new parameter value
     void SetCustomRateParameter(const std::string& label, double value);
     void SetCustomRateParameter(const std::string& label, const std::vector<double>& values);
+    void SetRelativeTolerances(double relativeTolerance);
 
     /// @brief Print a header of species to display concentrations with respect to time
     void PrintHeader();

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_)
@@ -185,6 +187,12 @@ namespace micm
     for (std::size_t i = 0; i < custom_rate_parameters_.NumRows(); ++i)
       custom_rate_parameters_[i][param->second] = values[i];
   }
+
+  template<class DenseMatrixPolicy, class SparseMatrixPolicy>
+  inline void State<DenseMatrixPolicy, SparseMatrixPolicy>::SetRelativeTolerances(double relativeTolerance)
+  {
+    this->relative_tolerance_ = relativeTolerance;
+  }
 
   template<class DenseMatrixPolicy, class SparseMatrixPolicy>
   inline void State<DenseMatrixPolicy, SparseMatrixPolicy>::PrintHeader()

include/micm/version.hpp

@@ -1,5 +1,3 @@
-// Copyright (C) 2023-2024 National Center for Atmospheric Research
-// SPDX-License-Identifier: Apache-2.0
 // clang-format off
 #pragma once