diff --git a/include/micm/cuda/solver/cuda_rosenbrock.cuh b/include/micm/cuda/solver/cuda_rosenbrock.cuh index 95cfd1e9a..ea00f3da1 100644 --- a/include/micm/cuda/solver/cuda_rosenbrock.cuh +++ b/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 \ No newline at end of file diff --git a/include/micm/cuda/solver/cuda_rosenbrock.hpp b/include/micm/cuda/solver/cuda_rosenbrock.hpp index 15cbfea6d..36b96d27e 100644 --- a/include/micm/cuda/solver/cuda_rosenbrock.hpp +++ b/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 diff --git a/include/micm/cuda/solver/cuda_state.hpp b/include/micm/cuda/solver/cuda_state.hpp index 38428a318..367dde531 100644 --- a/include/micm/cuda/solver/cuda_state.hpp +++ b/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(absolute_tolerance_param_, absolute_tolerance_param_.number_of_elements_), "cudaMalloc"); + CHECK_CUDA_ERROR(micm::cuda::CopyToDevice(absolute_tolerance_param_, this->absolute_tolerance_), "cudaMemcpyHostToDevice"); } /// @brief Copy output variables to the host diff --git a/include/micm/solver/backward_euler.inl b/include/micm/solver/backward_euler.inl index c11f7f759..fc6fa479d 100644 --- a/include/micm/solver/backward_euler.inl +++ b/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) diff --git a/include/micm/solver/rosenbrock.inl b/include/micm/solver/rosenbrock.inl index 939c9521c..2e6484b03 100644 --- a/include/micm/solver/rosenbrock.inl +++ b/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; } diff --git a/include/micm/solver/state.hpp b/include/micm/solver/state.hpp index 6fe8945c6..082edde4b 100644 --- a/include/micm/solver/state.hpp +++ b/include/micm/solver/state.hpp @@ -61,7 +61,7 @@ namespace micm std::size_t number_of_grid_cells_; std::unique_ptr temporary_variables_; double relative_tolerance_; - DenseMatrixPolicy absolute_tolerance_; + std::vector absolute_tolerance_; /// @brief Default constructor /// Only defined to be used to create default values in types, but a default constructed state is not useable diff --git a/include/micm/solver/state.inl b/include/micm/solver/state.inl index cd6a7e59b..6bc63897e 100644 --- a/include/micm/solver/state.inl +++ b/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++; diff --git a/src/solver/rosenbrock.cu b/src/solver/rosenbrock.cu index dc725b0a3..1fa0f453e 100644 --- a/src/solver/rosenbrock.cu +++ b/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; diff --git a/test/integration/analytical_policy.hpp b/test/integration/analytical_policy.hpp index 8554db6cf..44e50c690 100644 --- a/test/integration/analytical_policy.hpp +++ b/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(3, state.relative_tolerance_ * 1e-2); + state.absolute_tolerance_ = std::vector(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(5, state.relative_tolerance_ * 1e-6); + state.absolute_tolerance_ = std::vector(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(8, state.relative_tolerance_ * 1e-2); + state.absolute_tolerance_ = std::vector(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(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(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); diff --git a/test/integration/test_chapman_integration.cpp b/test/integration/test_chapman_integration.cpp index cbd876f41..545571a17 100644 --- a/test/integration/test_chapman_integration.cpp +++ b/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) { diff --git a/test/unit/cuda/solver/test_cuda_rosenbrock.cpp b/test/unit/cuda/solver/test_cuda_rosenbrock.cpp index b3bcaa7bb..eceb9d4a6 100644 --- a/test/unit/cuda/solver/test_cuda_rosenbrock.cpp +++ b/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 atol = state.absolute_tolerance_.AsVector(); + auto& atol = state.absolute_tolerance_; double rtol = state.relative_tolerance_; auto y_old = micm::CudaDenseMatrix(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 atol = state.absolute_tolerance_.AsVector(); + auto& atol = state.absolute_tolerance_; double rtol = state.relative_tolerance_; auto y_old = micm::CudaDenseMatrix(number_of_grid_cells, state.state_size_, 7.7); auto y_new = micm::CudaDenseMatrix(number_of_grid_cells, state.state_size_, -13.9); auto errors = micm::CudaDenseMatrix(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) { diff --git a/test/unit/solver/test_backward_euler.cpp b/test/unit/solver/test_backward_euler.cpp index dd6d4f1b6..eb976b350 100644 --- a/test/unit/solver/test_backward_euler.cpp +++ b/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; diff --git a/test/unit/solver/test_rosenbrock.cpp b/test/unit/solver/test_rosenbrock.cpp index ee760b54c..01aa8a510 100644 --- a/test/unit/solver/test_rosenbrock.cpp +++ b/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 atol = state.absolute_tolerance_.AsVector(); + std::vector 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); } }