correctly using tolerances

include/micm/cuda/solver/cuda_rosenbrock.cuh

@@ -42,7 +42,7 @@ namespace micm
         const CudaMatrixParam& y_new_param,
         const CudaMatrixParam& errors_param,
         const CudaMatrixParam& absolute_tolerance_param,
-        const double* relative_tolerance,
+        const double relative_tolerance,
         CudaRosenbrockSolverParam devstruct);
   }  // namespace cuda
 }  // namespace micm

include/micm/cuda/solver/cuda_rosenbrock.hpp

@@ -121,8 +121,8 @@ namespace micm
           y_old.AsDeviceParam(),
           y_new.AsDeviceParam(),
           errors.AsDeviceParam(),
-          state.absolute_tolerance_.AsDeviceParam(), 
-          state.cuda_relative_tolerance_,
+          state.absolute_tolerance_param_, 
+          state.relative_tolerance_,
           this->devstruct_);
     }
   };  // end CudaRosenbrockSolver

include/micm/cuda/solver/cuda_state.hpp

@@ -21,11 +21,12 @@ namespace micm
     CudaState(CudaState&&) = default;
     CudaState& operator=(CudaState&&) = default;
 
-    double* cuda_relative_tolerance_{ nullptr };
+    CudaMatrixParam absolute_tolerance_param_;
 
     ~CudaState()
     {
-      CHECK_CUDA_ERROR(cudaFree(cuda_relative_tolerance_), "cudaFree");
+      CHECK_CUDA_ERROR(micm::cuda::FreeVector(absolute_tolerance_param_), "cudaFree");
+      absolute_tolerance_param_.d_data_ = nullptr;
     }
 
     /// @brief Constructor which takes the state dimension information as input
@@ -38,14 +39,12 @@ namespace micm
     {
       this->variables_.CopyToDevice();
       this->rate_constants_.CopyToDevice();
-      this->absolute_tolerance_.CopyToDevice();
 
-      size_t relative_tolerance_bytes = sizeof(double);
+      absolute_tolerance_param_.number_of_elements_ = this->absolute_tolerance_.size();
+      absolute_tolerance_param_.number_of_grid_cells_ = 1;
 
-      CHECK_CUDA_ERROR(
-          cudaMallocAsync(
-              &(cuda_relative_tolerance_), relative_tolerance_bytes, micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)),
-          "cudaMalloc");
+      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_, this->absolute_tolerance_), "cudaMemcpyHostToDevice");
     }
 
     /// @brief Copy output variables to the host

include/micm/solver/backward_euler.inl

@@ -133,7 +133,7 @@ namespace micm
           continue;
 
         // check for convergence
-        converged = IsConverged(parameters_, forcing, Yn1, state.absolute_tolerance_.AsVector(), state.relative_tolerance_);
+        converged = IsConverged(parameters_, forcing, Yn1, state.absolute_tolerance_, state.relative_tolerance_);
       } while (!converged && iterations < max_iter);
 
       if (!converged)

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_.AsVector().size();
+    const std::size_t n_vars = state.absolute_tolerance_.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_.AsVector()[i % n_vars] + state.relative_tolerance_ * ymax);
+          _errors[i] / (state.absolute_tolerance_[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_.AsVector().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();
 
@@ -302,7 +302,7 @@ namespace micm
     for (std::size_t i = 0; i < whole_blocks; ++i)
     {
       errors_over_scale =
-          *errors_iter / (state.absolute_tolerance_.AsVector()[(i / L) % n_vars] +
+          *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;
@@ -321,7 +321,7 @@ namespace micm
         {
           const std::size_t idx = y * L + x;
           errors_over_scale = errors_iter[idx] /
-                              (state.absolute_tolerance_.AsVector()[y] +
+                              (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/state.hpp

@@ -61,7 +61,7 @@ namespace micm
     std::size_t number_of_grid_cells_;
     std::unique_ptr<TemporaryVariables> temporary_variables_;
     double relative_tolerance_;
-    DenseMatrixPolicy absolute_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

include/micm/solver/state.inl

@@ -82,10 +82,8 @@ namespace micm
         state_size_(parameters.variable_names_.size()),
         number_of_grid_cells_(parameters.number_of_grid_cells_),
         relative_tolerance_(parameters.relative_tolerance_),
-        absolute_tolerance_()
+        absolute_tolerance_(parameters.absolute_tolerance_)
   {
-    absolute_tolerance_.AsVector() = parameters.absolute_tolerance_;
-
     std::size_t index = 0;
     for (auto& name : variable_names_)
       variable_map_[name] = index++;

src/solver/rosenbrock.cu

@@ -117,7 +117,7 @@ namespace micm
         const CudaMatrixParam y_old_param,
         const CudaMatrixParam y_new_param,
         const CudaMatrixParam absolute_tolerance_param,
-        const double* relative_tolerance,
+        const double relative_tolerance,
         CudaRosenbrockSolverParam devstruct,
         const size_t n,
         bool is_first_call)
@@ -127,7 +127,7 @@ namespace micm
       double* const d_errors_input = devstruct.errors_input_;
       double* const d_errors_output = devstruct.errors_output_;
       const double* const atol = absolute_tolerance_param.d_data_;
-      const double rtol = *relative_tolerance;
+      const double rtol = relative_tolerance;
       const size_t number_of_grid_cells = y_old_param.number_of_grid_cells_;
 
       // Declares a dynamically-sized shared memory array.
@@ -226,7 +226,7 @@ namespace micm
         const CudaMatrixParam y_old_param,
         const CudaMatrixParam y_new_param,
         const CudaMatrixParam absolute_tolerance_param,
-        const double* relative_tolerance,
+        const double relative_tolerance,
         CudaRosenbrockSolverParam devstruct)
     {
       // Local device variables
@@ -235,7 +235,7 @@ namespace micm
       const double* const d_y_new = y_new_param.d_data_;
       double* const d_errors = devstruct.errors_input_;
       const double* const atol = absolute_tolerance_param.d_data_;
-      const double rtol = *relative_tolerance;
+      const double rtol = relative_tolerance;
       const size_t num_elements = devstruct.errors_size_;
       const size_t number_of_grid_cells = y_old_param.number_of_grid_cells_;
 
@@ -270,7 +270,7 @@ namespace micm
         const CudaMatrixParam& y_new_param,
         const CudaMatrixParam& errors_param,
         const CudaMatrixParam& absolute_tolerance_param,
-        const double* relative_tolerance,
+        const double relative_tolerance,
         CudaRosenbrockSolverParam devstruct)
     {
       double normalized_error;

test/integration/analytical_policy.hpp

@@ -1332,7 +1332,7 @@ void test_analytical_robertson(
 
   auto state = solver.GetState();
   state.relative_tolerance_ = 1e-10;
-  state.absolute_tolerance_.AsVector() = std::vector<double>(3, state.relative_tolerance_ * 1e-2);
+  state.absolute_tolerance_ = std::vector<double>(3, state.relative_tolerance_ * 1e-2);
 
   double k1 = 0.04;
   double k2 = 3e7;
@@ -1548,7 +1548,7 @@ void test_analytical_oregonator(
   auto state = solver.GetState();
 
   state.relative_tolerance_ = 1e-8;
-  state.absolute_tolerance_.AsVector() = std::vector<double>(5, state.relative_tolerance_ * 1e-6);
+  state.absolute_tolerance_ = std::vector<double>(5, state.relative_tolerance_ * 1e-6);
 
   state.SetCustomRateParameter("r1", 1.34 * 0.06);
   state.SetCustomRateParameter("r2", 1.6e9);
@@ -1726,7 +1726,7 @@ void test_analytical_hires(
 
   auto state = solver.GetState();
   state.relative_tolerance_ = 1e-6;
-  state.absolute_tolerance_.AsVector() = std::vector<double>(8, state.relative_tolerance_ * 1e-2);
+  state.absolute_tolerance_ = std::vector<double>(8, state.relative_tolerance_ * 1e-2);
 
   state.SetCustomRateParameter("r1", 1.71);
   state.SetCustomRateParameter("r2", 8.75);
@@ -1886,10 +1886,10 @@ void test_analytical_e5(
   auto state = solver.GetState();
 
   state.relative_tolerance_ = 1e-13;
-  state.absolute_tolerance_.AsVector() = std::vector<double>(6, 1e-17);
-  state.absolute_tolerance_[0][0] = 1e-7;
-  state.absolute_tolerance_[0][4] = 1e-7;
-  state.absolute_tolerance_[0][5] = 1e-7;
+  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);

test/integration/test_chapman_integration.cpp

@@ -43,7 +43,7 @@ TEST(ChapmanIntegration, CanBuildChapmanSystemUsingConfig)
   micm::State state = solver.GetState();
   state.SetRelativeTolerances(solver_params.relative_tolerance_); 
   EXPECT_EQ(state.relative_tolerance_, 1.0e-4);
-  auto& abs_tol = state.absolute_tolerance_.AsVector();
+  auto& abs_tol = state.absolute_tolerance_;
 
   for (size_t n_grid_cell = 0; n_grid_cell < state.number_of_grid_cells_; ++n_grid_cell)
   {

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_.AsVector();
+  auto& atol = state.absolute_tolerance_;
   double rtol = state.relative_tolerance_;
 
   auto y_old = micm::CudaDenseMatrix<double, L>(number_of_grid_cells, state.state_size_, 1.0);
@@ -133,6 +133,8 @@ void testNormalizedErrorConst()
   y_new.CopyToDevice();
   errors.CopyToDevice();
 
+  state.SyncInputsToDevice();
+
   double error = gpu_solver.solver_.NormalizedError(y_old, y_new, errors, state);
 
   auto expected_error = 0.0;
@@ -166,12 +168,14 @@ 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_.AsVector();  
+  auto& atol = state.absolute_tolerance_;  
   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);
   auto errors = micm::CudaDenseMatrix<double, L>(number_of_grid_cells, state.state_size_, 81.57);
 
+  state.SyncInputsToDevice();
+
   double expected_error = 0.0;
   for (size_t i = 0; i < number_of_grid_cells; ++i)
   {

test/unit/solver/test_backward_euler.cpp

@@ -48,7 +48,7 @@ TEST(BackwardEuler, CanCallSolve)
   double time_step = 1.0;
 
   auto state = be.GetState();
-  state.absolute_tolerance_.AsVector() = { 1e-6, 1e-6, 1e-6 };
+  state.absolute_tolerance_ = { 1e-6, 1e-6, 1e-6 };
 
   state.variables_[0] = { 1.0, 0.0, 0.0 };
   state.conditions_[0].temperature_ = 272.5;

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_.AsVector();  
+  std::vector<double> atol = state.absolute_tolerance_;  
   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_.AsVector().size(), 2);    
-    EXPECT_EQ(state.absolute_tolerance_.AsVector()[0], 1.0e-07);
-    EXPECT_EQ(state.absolute_tolerance_.AsVector()[1], 1.0e-08);
+    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);
   }
 }