Remove solver_.parameters_ (#668)

Moves the absolute and relative tolerances off of the solver parameters and onto the state iwth other mechanism-specific data --------- Co-authored-by: Kyle Shores <[email protected]> Co-authored-by: Kyle Shores <[email protected]> Co-authored-by: Jian Sun <[email protected]>
NCAR · Jan 14, 2025 · 71302b7 · 71302b7
1 parent faf21bd
commit 71302b7
Show file tree

Hide file tree

Showing 39 changed files with 831 additions and 1,063 deletions.
diff --git a/examples/profile_example.cpp b/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

diff --git a/include/micm/configure/solver_config.hpp b/include/micm/configure/solver_config.hpp
@@ -102,18 +102,28 @@ 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),
           processes_(processes),
-          parameters_(parameters)
+          parameters_(parameters),
+          relative_tolerance_(0)
     {
     }
 
     SolverParameters(System&& system, std::vector<Process>&& processes, RosenbrockSolverParameters&& parameters)
         : system_(system),
           processes_(processes),
-          parameters_(parameters)
+          parameters_(parameters),
+          relative_tolerance_(0)
+    {
+    }
+    SolverParameters(System&& system, std::vector<Process>&& processes, RosenbrockSolverParameters&& parameters, double relative_tolerance)
+        : system_(system),
+          processes_(processes),
+          parameters_(parameters), 
+          relative_tolerance_(relative_tolerance)
     {
     }
   };
@@ -138,6 +148,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 +284,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 +437,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 +970,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_), this->relative_tolerance_);
     }
   };
 }  // namespace micm
diff --git a/include/micm/cuda/solver/cuda_rosenbrock.cuh b/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
diff --git a/include/micm/cuda/solver/cuda_rosenbrock.hpp b/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;
     };
 
@@ -66,24 +65,22 @@ namespace micm
     /// @param rates Rates calculator
     /// @param jacobian Jacobian matrix
     CudaRosenbrockSolver(
-        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,12 +112,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)
-        const
-      requires(CudaMatrix<DenseMatrixPolicy> && VectorizableDense<DenseMatrixPolicy>)
+    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.relative_tolerance_,
+          this->devstruct_);
     }
   };  // end CudaRosenbrockSolver
 }  // namespace micm
diff --git a/include/micm/cuda/solver/cuda_state.hpp b/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, LuDecompositionPolicy>(parameters){};
+        : State<DenseMatrixPolicy, SparseMatrixPolicy, LuDecompositionPolicy>(parameters)
+    {
+      auto& atol = this->absolute_tolerance_;
+
+      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, LuDecompositionPolicy>(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, LuDecompositionPolicy>::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, LuDecompositionPolicy>::SetAbsoluteTolerances(absoluteTolerance);
+      CHECK_CUDA_ERROR(micm::cuda::CopyToDevice<double>(absolute_tolerance_param_, absoluteTolerance), "cudaMemcpyHostToDevice");
+    }
 
     /// @brief Copy input variables to the device
     void SyncInputsToDevice()

diff --git a/include/micm/cuda/util/cuda_matrix.cuh b/include/micm/cuda/util/cuda_matrix.cuh
@@ -41,7 +41,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

diff --git a/include/micm/cuda/util/cuda_param.hpp b/include/micm/cuda/util/cuda_param.hpp
@@ -90,8 +90,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;

diff --git a/include/micm/jit/solver/jit_rosenbrock.hpp b/include/micm/jit/solver/jit_rosenbrock.hpp
@@ -65,20 +65,18 @@ namespace micm
     }
 
     /// @brief Builds a Rosenbrock solver for the given system and solver parameters
-    /// @param parameters Solver parameters
     /// @param linear_solver Linear solver
     /// @param rates Rates calculator
     /// @param jacobian Jacobian matrix
     JitRosenbrockSolver(
-        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);
     }

diff --git a/include/micm/solver/backward_euler.hpp b/include/micm/solver/backward_euler.hpp
@@ -32,7 +32,6 @@ namespace micm
   class BackwardEuler
   {
    public:
-    BackwardEulerSolverParameters parameters_;
     LinearSolverPolicy linear_solver_;
     RatesPolicy rates_;
     std::vector<std::size_t> jacobian_diagonal_elements_;
@@ -41,17 +40,15 @@ namespace micm
     using ParametersType = BackwardEulerSolverParameters;
 
     /// @brief Default constructor
-    /// @param parameters Solver parameters
     /// @param linear_solver Linear solver
     /// @param rates Rates calculator
     /// @param jacobian Jacobian matrix
     BackwardEuler(
-        const BackwardEulerSolverParameters& parameters,
         LinearSolverPolicy&& linear_solver,
         RatesPolicy&& rates,
-        auto& jacobian)
-        : parameters_(parameters),
-          linear_solver_(std::move(linear_solver)),
+        auto& jacobian,
+        const size_t number_of_species)
+        : linear_solver_(std::move(linear_solver)),
           rates_(std::move(rates)),
           jacobian_diagonal_elements_(jacobian.DiagonalIndices(0))
     {
@@ -68,26 +65,23 @@ namespace micm
     /// @param time_step Time [s] to advance the state by
     /// @param state The state to advance
     /// @return result of the solver (success or failure, and statistics)
-    SolverResult Solve(double time_step, auto& state) const;
+    SolverResult Solve(double time_step, auto& state, const BackwardEulerSolverParameters& parameters) const;
 
     /// @brief Determines whether the residual is small enough to stop the
     ///        internal solver iteration
-    /// @param parameters Solver parameters
     /// @param residual The residual to check
     /// @param state The current state being solved for
     /// @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