updating usage of absolute tolerance

include/micm/cuda/solver/cuda_rosenbrock.hpp

@@ -69,17 +69,19 @@ 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();
       // Copy the data from host struct to device struct
@@ -117,7 +119,12 @@ namespace micm
         const requires(CudaMatrix<DenseMatrixPolicy>&& VectorizableDense<DenseMatrixPolicy>)
     {
       return micm::cuda::NormalizedErrorDriver(
-          y_old.AsDeviceParam(), y_new.AsDeviceParam(), errors.AsDeviceParam(), this->parameters_, state.absolute_tolerance_.AsDeviceParam(), this->devstruct_);
+          y_old.AsDeviceParam(),
+          y_new.AsDeviceParam(),
+          errors.AsDeviceParam(),
+          state.absolute_tolerance_.AsDeviceParam(), 
+          this->parameters_, 
+          this->devstruct_);
     }
   };  // end CudaRosenbrockSolver
 }  // namespace micm

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)),
@@ -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;

include/micm/solver/rosenbrock.inl

@@ -254,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 = state.absolute_tolerance_.size();
+    const std::size_t n_vars = state.absolute_tolerance_.AsVector().size();
 
     double ymax = 0;
     double errors_over_scale = 0;
@@ -264,7 +264,7 @@ namespace micm
     {
       ymax = std::max(std::abs(_y[i]), std::abs(_ynew[i]));
       errors_over_scale =
-          _errors[i] / (state.absolute_tolerance_[i % n_vars] + state.relative_tolerance_ * ymax);
+          _errors[i] / (state.absolute_tolerance_.AsVector()[i % n_vars] + state.relative_tolerance_ * ymax);
       error += errors_over_scale * errors_over_scale;
     }
 
@@ -291,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 = state.absolute_tolerance_.size();
+    const std::size_t n_vars = state.absolute_tolerance_.AsVector().size();
 
     const std::size_t whole_blocks = std::floor(Y.NumRows() / Y.GroupVectorSize()) * Y.GroupSize();
 
@@ -302,7 +302,7 @@ namespace micm
     for (std::size_t i = 0; i < whole_blocks; ++i)
     {
       errors_over_scale =
-          *errors_iter / (state.absolute_tolerance_[(i / L) % n_vars] +
+          *errors_iter / (state.absolute_tolerance_.AsVector()[(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;
@@ -321,7 +321,7 @@ namespace micm
         {
           const std::size_t idx = y * L + x;
           errors_over_scale = errors_iter[idx] /
-                              (state.absolute_tolerance_[y] +
+                              (state.absolute_tolerance_.AsVector()[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/solver_builder.inl

@@ -402,7 +402,7 @@ namespace micm
     this->SetAbsoluteTolerances(state_parameters.absolute_tolerance_, species_map);
 
     return Solver<SolverPolicy, StatePolicy>(
-        SolverPolicy(options, std::move(linear_solver), std::move(rates), jacobian),
+        SolverPolicy(options, std::move(linear_solver), std::move(rates), jacobian, number_of_species),
         state_parameters,
         options,
         this->reactions_);

test/unit/cuda/solver/test_cuda_rosenbrock.cpp

@@ -122,7 +122,7 @@ void testNormalizedErrorConst()
   gpu_builder = getSolver(gpu_builder);
   auto gpu_solver = gpu_builder.SetNumberOfGridCells(number_of_grid_cells).Build();
   auto state = gpu_solver.GetState();  
-  std::vector<double> atol = state.absolute_tolerance_;
+  std::vector<double> atol = state.absolute_tolerance_.AsVector();
   double rtol = state.relative_tolerance_;
 
   auto y_old = micm::CudaDenseMatrix<double, L>(number_of_grid_cells, state.state_size_, 1.0);
@@ -166,7 +166,7 @@ void testNormalizedErrorDiff()
   gpu_builder = getSolver(gpu_builder);
   auto gpu_solver = gpu_builder.SetNumberOfGridCells(number_of_grid_cells).Build();
   auto state = gpu_solver.GetState();
-  std::vector<double> atol = state.absolute_tolerance_;  
+  std::vector<double> atol = state.absolute_tolerance_.AsVector();  
   double rtol = state.relative_tolerance_;
   auto y_old = micm::CudaDenseMatrix<double, L>(number_of_grid_cells, state.state_size_, 7.7);
   auto y_new = micm::CudaDenseMatrix<double, L>(number_of_grid_cells, state.state_size_, -13.9);

test/unit/solver/test_rosenbrock.cpp

@@ -18,7 +18,7 @@ void testNormalizedErrorDiff(SolverBuilderPolicy builder, std::size_t number_of_
   builder = getSolver(builder);
   auto solver = builder.SetNumberOfGridCells(number_of_grid_cells).Build();  
   auto state = solver.GetState();
-  std::vector<double> atol = state.absolute_tolerance_;  
+  std::vector<double> atol = state.absolute_tolerance_.AsVector();  
   double rtol = state.relative_tolerance_;
 
   using MatrixPolicy = decltype(state.variables_);
@@ -107,9 +107,9 @@ TEST(RosenbrockSolver, CanSetTolerances)
                       .SetNumberOfGridCells(number_of_grid_cells)
                       .Build();
     auto state = solver.GetState();
-    EXPECT_EQ(state.absolute_tolerance_.size(), 2);    
-    EXPECT_EQ(state.absolute_tolerance_[0], 1.0e-07);
-    EXPECT_EQ(state.absolute_tolerance_[1], 1.0e-08);
+    EXPECT_EQ(state.absolute_tolerance_.AsVector().size(), 2);    
+    EXPECT_EQ(state.absolute_tolerance_.AsVector()[0], 1.0e-07);
+    EXPECT_EQ(state.absolute_tolerance_.AsVector()[1], 1.0e-08);
   }
 }