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.SetRelativeTolerance(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.cuh

@@ -41,7 +41,8 @@ namespace micm
         const CudaMatrixParam& y_old_param,
         const CudaMatrixParam& y_new_param,
         const CudaMatrixParam& errors_param,
-        const RosenbrockSolverParameters& ros_param,
+        const CudaMatrixParam& absolute_tolerance_param,
+        const double relative_tolerance,
         CudaRosenbrockSolverParam devstruct);
   }  // namespace cuda
 }  // namespace micm

include/micm/cuda/solver/cuda_rosenbrock.hpp

@@ -38,7 +38,6 @@ namespace micm
     {
       other.devstruct_.errors_input_ = nullptr;
       other.devstruct_.errors_output_ = nullptr;
-      other.devstruct_.absolute_tolerance_ = nullptr;
       other.devstruct_.jacobian_diagonal_elements_ = nullptr;
     };
 
@@ -69,21 +68,21 @@ namespace micm
         RosenbrockSolverParameters parameters,
         LinearSolverPolicy&& linear_solver,
         RatesPolicy&& rates,
-        auto& jacobian)
+        auto& jacobian,
+        const size_t number_of_species)
         : AbstractRosenbrockSolver<RatesPolicy, LinearSolverPolicy, CudaRosenbrockSolver<RatesPolicy, LinearSolverPolicy>>(
               parameters,
               std::move(linear_solver),
               std::move(rates),
-              jacobian)
+              jacobian,
+              number_of_species)
     {
       CudaRosenbrockSolverParam hoststruct;
       // jacobian.GroupVectorSize() is the same as the number of grid cells for the CUDA implementation
       // the absolute tolerance size is the same as the number of solved variables in one grid cell
-      hoststruct.errors_size_ = jacobian.GroupVectorSize() * this->parameters_.absolute_tolerance_.size();
+      hoststruct.errors_size_ = jacobian.GroupVectorSize() * number_of_species;
       hoststruct.jacobian_diagonal_elements_ = this->jacobian_diagonal_elements_.data();
       hoststruct.jacobian_diagonal_elements_size_ = this->jacobian_diagonal_elements_.size();
-      hoststruct.absolute_tolerance_ = this->parameters_.absolute_tolerance_.data();
-      hoststruct.absolute_tolerance_size_ = this->parameters_.absolute_tolerance_.size();
       // Copy the data from host struct to device struct
       this->devstruct_ = micm::cuda::CopyConstData(hoststruct);
     };
@@ -115,11 +114,16 @@ 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(
-          y_old.AsDeviceParam(), y_new.AsDeviceParam(), errors.AsDeviceParam(), this->parameters_, this->devstruct_);
+          y_old.AsDeviceParam(),
+          y_new.AsDeviceParam(),
+          errors.AsDeviceParam(),
+          state.absolute_tolerance_param_, 
+          state.GetRelativeTolerance(),
+          this->devstruct_);
     }
   };  // end CudaRosenbrockSolver
 }  // namespace micm

include/micm/cuda/solver/cuda_state.hpp

@@ -16,13 +16,53 @@ namespace micm
    public:
     CudaState(const CudaState&) = delete;
     CudaState& operator=(const CudaState&) = delete;
-    CudaState(CudaState&&) = default;
-    CudaState& operator=(CudaState&&) = default;
+
+    CudaMatrixParam absolute_tolerance_param_;
+
+    ~CudaState()
+    {
+      CHECK_CUDA_ERROR(micm::cuda::FreeVector(absolute_tolerance_param_), "cudaFree");
+    }
 
     /// @brief Constructor which takes the state dimension information as input
     /// @param parameters State dimension information
     CudaState(const StateParameters& parameters)
-        : State<DenseMatrixPolicy, SparseMatrixPolicy>(parameters){};
+        : State<DenseMatrixPolicy, SparseMatrixPolicy>(parameters)
+    {
+      const auto& atol = this->GetAbsoluteTolerances();
+
+      absolute_tolerance_param_.number_of_elements_ = atol.size();
+      absolute_tolerance_param_.number_of_grid_cells_ = 1;
+
+      CHECK_CUDA_ERROR(micm::cuda::MallocVector<double>(absolute_tolerance_param_, absolute_tolerance_param_.number_of_elements_), "cudaMalloc");
+      CHECK_CUDA_ERROR(micm::cuda::CopyToDevice<double>(absolute_tolerance_param_, atol), "cudaMemcpyHostToDevice");
+    };
+
+    /// @brief Move constructor
+    CudaState(CudaState&& other)
+        : State<DenseMatrixPolicy, SparseMatrixPolicy>(std::move(other))
+    {
+      absolute_tolerance_param_ = other.absolute_tolerance_param_;
+      other.absolute_tolerance_param_.d_data_ = nullptr;
+    }
+
+    /// @brief Move assignment operator
+    CudaState& operator=(CudaState&& other)
+    {
+      if (this != &other)
+      {
+        State<DenseMatrixPolicy, SparseMatrixPolicy>::operator=(std::move(other));
+        absolute_tolerance_param_ = other.absolute_tolerance_param_;
+        other.absolute_tolerance_param_.d_data_ = nullptr;
+      }
+      return *this;
+    }
+
+    void SetAbsoluteTolerances(const std::vector<double>& absoluteTolerance) override
+    {
+      State<DenseMatrixPolicy, SparseMatrixPolicy>::SetAbsoluteTolerances(absoluteTolerance);
+      CHECK_CUDA_ERROR(micm::cuda::CopyToDevice<double>(absolute_tolerance_param_, absoluteTolerance), "cudaMemcpyHostToDevice");
+    }
 
     /// @brief Copy input variables to the device
     void SyncInputsToDevice() requires(CudaMatrix<DenseMatrixPolicy>&& VectorizableDense<DenseMatrixPolicy>)

include/micm/cuda/util/cuda_matrix.cuh

@@ -29,7 +29,7 @@ namespace micm
     /// @param h_data Host data to copy from
     /// @returns Error code from copying to device from the host, if any
     template<typename T>
-    cudaError_t CopyToDevice(CudaMatrixParam& vectorMatrix, std::vector<T>& h_data);
+    cudaError_t CopyToDevice(CudaMatrixParam& vectorMatrix, const std::vector<T>& h_data);
 
     /// @brief Copies data from the device to the host
     /// @param vectorMatrix Struct containing allocated device memory

include/micm/cuda/util/cuda_param.hpp

@@ -87,8 +87,6 @@ struct CudaRosenbrockSolverParam
   // for NormalizedError function
   double* errors_input_ = nullptr;
   double* errors_output_ = nullptr;
-  double* absolute_tolerance_ = nullptr;
-  std::size_t absolute_tolerance_size_;
   std::size_t errors_size_;
   // for AlphaMinusJacobian function
   std::size_t* jacobian_diagonal_elements_ = nullptr;

include/micm/jit/solver/jit_rosenbrock.hpp

@@ -73,12 +73,13 @@ namespace micm
         RosenbrockSolverParameters parameters,
         LinearSolverPolicy linear_solver,
         RatesPolicy rates,
-        auto& jacobian)
+        auto& jacobian,
+        const size_t number_of_species)
         : AbstractRosenbrockSolver<RatesPolicy, LinearSolverPolicy, JitRosenbrockSolver<RatesPolicy, LinearSolverPolicy>>(
               parameters,
               std::move(linear_solver),
               std::move(rates),
-              jacobian)
+              jacobian, number_of_species)
     {
       this->GenerateAlphaMinusJacobian(jacobian);
     }

include/micm/solver/backward_euler.hpp

@@ -49,7 +49,8 @@ namespace micm
         const BackwardEulerSolverParameters& parameters,
         LinearSolverPolicy&& linear_solver,
         RatesPolicy&& rates,
-        auto& jacobian)
+        auto& jacobian,
+        const size_t number_of_species)
         : parameters_(parameters),
           linear_solver_(std::move(linear_solver)),
           rates_(std::move(rates)),
@@ -78,14 +79,14 @@ namespace micm
     /// @return true if the residual is small enough to stop the iteration
     template<class DenseMatrixPolicy>
     static bool IsConverged(
-        const BackwardEulerSolverParameters& parameters,
-        const DenseMatrixPolicy& residual,
-        const DenseMatrixPolicy& state) requires(!VectorizableDense<DenseMatrixPolicy>);
+      const BackwardEulerSolverParameters& parameters,
+      const DenseMatrixPolicy& residual,
+      const DenseMatrixPolicy& Yn1, const std::vector<double>& absolute_tolerance, double relative_tolerance) requires(!VectorizableDense<DenseMatrixPolicy>);
     template<class DenseMatrixPolicy>
     static bool IsConverged(
-        const BackwardEulerSolverParameters& parameters,
-        const DenseMatrixPolicy& residual,
-        const DenseMatrixPolicy& state) requires(VectorizableDense<DenseMatrixPolicy>);
+      const BackwardEulerSolverParameters& parameters,
+      const DenseMatrixPolicy& residual,
+      const DenseMatrixPolicy& Yn1, const std::vector<double>& absolute_tolerance, double relative_tolerance) 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.GetAbsoluteTolerances(), state.GetRelativeTolerance());
       } while (!converged && iterations < max_iter);
 
       if (!converged)
@@ -180,23 +180,23 @@ namespace micm
   inline bool BackwardEuler<RatesPolicy, LinearSolverPolicy>::IsConverged(
       const BackwardEulerSolverParameters& parameters,
       const DenseMatrixPolicy& residual,
-      const DenseMatrixPolicy& state) requires(!VectorizableDense<DenseMatrixPolicy>)
+      const DenseMatrixPolicy& Yn1, const std::vector<double>& absolute_tolerance, double relative_tolerance) requires(!VectorizableDense<DenseMatrixPolicy>)
   {
     double small = parameters.small_;
-    double rel_tol = parameters.relative_tolerance_;
-    auto& abs_tol = parameters.absolute_tolerance_;
+    double rel_tol = relative_tolerance;
+    auto& abs_tol = absolute_tolerance;
     auto residual_iter = residual.AsVector().begin();
-    auto state_iter = state.AsVector().begin();
+    auto Yn1_iter = Yn1.AsVector().begin();
     const std::size_t n_elem = residual.NumRows() * residual.NumColumns();
     const std::size_t n_vars = abs_tol.size();
     for (std::size_t i = 0; i < n_elem; ++i)
     {
       if (std::abs(*residual_iter) > small && std::abs(*residual_iter) > abs_tol[i % n_vars] &&
-          std::abs(*residual_iter) > rel_tol * std::abs(*state_iter))
+          std::abs(*residual_iter) > rel_tol * std::abs(*Yn1_iter))
       {
         return false;
       }
-      ++residual_iter, ++state_iter;
+      ++residual_iter, ++Yn1_iter;
     }
     return true;
   }
@@ -206,27 +206,26 @@ namespace micm
   inline bool BackwardEuler<RatesPolicy, LinearSolverPolicy>::IsConverged(
       const BackwardEulerSolverParameters& parameters,
       const DenseMatrixPolicy& residual,
-      const DenseMatrixPolicy& state) requires(VectorizableDense<DenseMatrixPolicy>)
+      const DenseMatrixPolicy& Yn1, const std::vector<double>& absolute_tolerance, double relative_tolerance) requires(VectorizableDense<DenseMatrixPolicy>)
   {
     double small = parameters.small_;
-    double rel_tol = parameters.relative_tolerance_;
-    auto& abs_tol = parameters.absolute_tolerance_;
+    double rel_tol = relative_tolerance;
+    auto& abs_tol = absolute_tolerance;
     auto residual_iter = residual.AsVector().begin();
-    auto state_iter = state.AsVector().begin();
+    auto Yn1_iter = Yn1.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)
     {
       if (std::abs(*residual_iter) > small && std::abs(*residual_iter) > abs_tol[(i / L) % n_vars] &&
-          std::abs(*residual_iter) > rel_tol * std::abs(*state_iter))
+          std::abs(*residual_iter) > rel_tol * std::abs(*Yn1_iter))
       {
         return false;
       }
-      ++residual_iter, ++state_iter;
+      ++residual_iter, ++Yn1_iter;
     }
 
     // evaluate the remaining rows
@@ -239,7 +238,7 @@ namespace micm
         for (std::size_t i = offset; i < offset + remaining_rows; ++i)
         {
           if (std::abs(residual_iter[i]) > small && std::abs(residual_iter[i]) > abs_tol[y] &&
-              std::abs(residual_iter[i]) > rel_tol * std::abs(state_iter[i]))
+              std::abs(residual_iter[i]) > rel_tol * std::abs(Yn1_iter[i]))
           {
             return false;
           }

include/micm/solver/backward_euler_solver_parameters.hpp

@@ -19,8 +19,6 @@ namespace micm
     template<class RatesPolicy, class LinearSolverPolicy>
     using SolverType = BackwardEuler<RatesPolicy, LinearSolverPolicy>;
 
-    std::vector<double> absolute_tolerance_;
-    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

@@ -72,7 +72,8 @@ namespace micm
         const RosenbrockSolverParameters& parameters,
         LinearSolverPolicy&& linear_solver,
         RatesPolicy&& rates,
-        auto& jacobian)
+        auto& jacobian,
+        const size_t number_of_species)
         : parameters_(parameters),
           linear_solver_(std::move(linear_solver)),
           rates_(std::move(rates)),
@@ -116,10 +117,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
 
@@ -139,12 +140,14 @@ namespace micm
         const RosenbrockSolverParameters& parameters,
         LinearSolverPolicy&& linear_solver,
         RatesPolicy&& rates,
-        auto& jacobian)
+        auto& jacobian,
+        const size_t number_of_species)
         : AbstractRosenbrockSolver<RatesPolicy, LinearSolverPolicy, RosenbrockSolver<RatesPolicy, LinearSolverPolicy>>(
               parameters,
               std::move(linear_solver),
               std::move(rates),
-              jacobian)
+              jacobian,
+              number_of_species)
     {
     }
     RosenbrockSolver(const RosenbrockSolver&) = delete;