From 263af0741638f6d61b90c7a4dcee7e697410e05a Mon Sep 17 00:00:00 2001 From: Montek Thind Date: Tue, 24 Sep 2024 15:37:39 -0700 Subject: [PATCH 01/31] draft of 621 --- include/micm/solver/solver.hpp | 7 +++++++ test/unit/cuda/solver/test_cuda_rosenbrock.cpp | 4 ++-- test/unit/jit/solver/test_jit_rosenbrock.cpp | 4 ++-- test/unit/solver/test_rosenbrock.cpp | 11 +++++------ 4 files changed, 16 insertions(+), 10 deletions(-) diff --git a/include/micm/solver/solver.hpp b/include/micm/solver/solver.hpp index f8fc816d0..a387866c0 100644 --- a/include/micm/solver/solver.hpp +++ b/include/micm/solver/solver.hpp @@ -63,6 +63,13 @@ namespace micm return solver_.Solve(time_step, state); } + // Overloaded Solve function to change parameters + SolverResult Solve(double time_step, StatePolicy& state, RosenbrockSolverParameters params) + { + solver_.parameters_.h_start_ = params.h_start_; + return solver_.Solve(time_step, state); + } + /// @brief Returns the number of grid cells /// @return std::size_t GetNumberOfGridCells() const diff --git a/test/unit/cuda/solver/test_cuda_rosenbrock.cpp b/test/unit/cuda/solver/test_cuda_rosenbrock.cpp index f40d11ed5..825d111c1 100644 --- a/test/unit/cuda/solver/test_cuda_rosenbrock.cpp +++ b/test/unit/cuda/solver/test_cuda_rosenbrock.cpp @@ -284,12 +284,12 @@ TEST(RosenbrockSolver, SingularSystemZeroInBottomRightOfU) // so H needs to be 1 / ( (-k1 - k2) * gamma) // since H is positive we need -k1 -k2 to be positive, hence the smaller, negative value for k1 double H = 1 / ((-k1 - k2) * params.gamma_[0]); - vector_solver.solver_.parameters_.h_start_ = H; + params.h_start_ = H; vector_solver.CalculateRateConstants(vector_state); vector_state.SyncInputsToDevice(); - auto vector_result = vector_solver.Solve(2 * H, vector_state); + auto vector_result = vector_solver.Solve(2 * H, vector_state, params); vector_state.SyncOutputsToHost(); EXPECT_NE(vector_result.stats_.singular_, 0); } diff --git a/test/unit/jit/solver/test_jit_rosenbrock.cpp b/test/unit/jit/solver/test_jit_rosenbrock.cpp index 542bde6f3..4a1d8eaff 100644 --- a/test/unit/jit/solver/test_jit_rosenbrock.cpp +++ b/test/unit/jit/solver/test_jit_rosenbrock.cpp @@ -105,11 +105,11 @@ TEST(JitRosenbrockSolver, SingularSystemZeroInBottomRightOfU) // so H needs to be 1 / ( (-k1 - k2) * gamma) // since H is positive we need -k1 -k2 to be positive, hence the smaller, negative value for k1 double H = 1 / ((-k1 - k2) * params.gamma_[0]); - vector_solver.solver_.parameters_.h_start_ = H; + params.h_start_ = H; vector_solver.CalculateRateConstants(vector_state); - auto vector_result = vector_solver.Solve(2 * H, vector_state); + auto vector_result = vector_solver.Solve(2 * H, vector_state, params); EXPECT_NE(vector_result.stats_.singular_, 0); } diff --git a/test/unit/solver/test_rosenbrock.cpp b/test/unit/solver/test_rosenbrock.cpp index 06f3cd4bf..568b8e6b9 100644 --- a/test/unit/solver/test_rosenbrock.cpp +++ b/test/unit/solver/test_rosenbrock.cpp @@ -17,7 +17,7 @@ void testNormalizedErrorDiff(SolverBuilderPolicy builder, std::size_t number_of_ { builder = getSolver(builder); auto solver = builder.SetNumberOfGridCells(number_of_grid_cells).Build(); - std::vector atol = solver.solver_.parameters_.absolute_tolerance_; + std::vector atol = solver.solver_.parameters_.absolute_tolerance_; double rtol = solver.solver_.parameters_.relative_tolerance_; auto state = solver.GetState(); @@ -172,18 +172,17 @@ TEST(RosenbrockSolver, SingularSystemZeroInBottomRightOfU) // alpha is 1 / (H * gamma), where H is the time step and gamma is the gamma value from // the rosenbrock paramters // so H needs to be 1 / ( (-k1 - k2) * gamma) - // since H is positive we need -k1 -k2 to be positive, hence the smaller, negative value for k1 + // since H is positive we need -k1 -k2 to be positive, hence the smaller, negative value for k1 double H = 1 / ((-k1 - k2) * params.gamma_[0]); - standard_solver.solver_.parameters_.h_start_ = H; - vector_solver.solver_.parameters_.h_start_ = H; + params.h_start_ = H; standard_solver.CalculateRateConstants(standard_state); vector_solver.CalculateRateConstants(vector_state); - auto standard_result = standard_solver.Solve(2 * H, standard_state); + auto standard_result = standard_solver.Solve(2 * H, standard_state, params); EXPECT_NE(standard_result.stats_.singular_, 0); - auto vector_result = vector_solver.Solve(2 * H, vector_state); + auto vector_result = vector_solver.Solve(2 * H, vector_state, params); EXPECT_NE(vector_result.stats_.singular_, 0); } From da7cd22ff494a4f81af5bb58fe84bbaf582af5b5 Mon Sep 17 00:00:00 2001 From: Montek Thind Date: Mon, 30 Sep 2024 13:11:51 -0700 Subject: [PATCH 02/31] second draft --- include/micm/configure/solver_config.hpp | 2 +- include/micm/solver/backward_euler.hpp | 4 +-- include/micm/solver/backward_euler.inl | 14 +++++----- .../backward_euler_solver_parameters.hpp | 3 +-- include/micm/solver/rosenbrock.hpp | 4 +-- include/micm/solver/rosenbrock.inl | 22 +++++++++------- .../solver/rosenbrock_solver_parameters.hpp | 11 ++------ include/micm/solver/solver.hpp | 6 ++--- include/micm/solver/solver_builder.inl | 6 +++-- include/micm/solver/state.hpp | 8 ++++++ include/micm/solver/state.inl | 4 ++- test/integration/analytical_policy.hpp | 13 ++++++++++ .../cuda/test_cuda_analytical_rosenbrock.cpp | 12 ++++----- .../jit/test_jit_analytical_rosenbrock.cpp | 12 ++++----- test/integration/terminator.hpp | 1 + .../test_analytical_rosenbrock.cpp | 26 +++++++++---------- test/integration/test_chapman_integration.cpp | 22 ++++++++-------- test/integration/test_terminator.cpp | 4 +-- .../unit/cuda/solver/test_cuda_rosenbrock.cpp | 17 +++++------- test/unit/solver/test_backward_euler.cpp | 5 ++-- test/unit/solver/test_rosenbrock.cpp | 17 ++++++------ test/unit/solver/test_solver_builder.cpp | 4 ++- 22 files changed, 119 insertions(+), 98 deletions(-) diff --git a/include/micm/configure/solver_config.hpp b/include/micm/configure/solver_config.hpp index c23c5a795..0358e8cf6 100644 --- a/include/micm/configure/solver_config.hpp +++ b/include/micm/configure/solver_config.hpp @@ -426,7 +426,7 @@ namespace micm void ParseRelativeTolerance(const objectType& object) { ValidateSchema(object, { "value", "type" }, {}); - this->parameters_.relative_tolerance_ = object["value"].as(); + //this->parameters_.relative_tolerance_ = object["value"].as(); } void ParseMechanism(const objectType& object) diff --git a/include/micm/solver/backward_euler.hpp b/include/micm/solver/backward_euler.hpp index 808ce8d05..a0bfdf861 100644 --- a/include/micm/solver/backward_euler.hpp +++ b/include/micm/solver/backward_euler.hpp @@ -80,12 +80,12 @@ namespace micm static bool IsConverged( const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& state) requires(!VectorizableDense); + const DenseMatrixPolicy& state, auto& stateParams) requires(!VectorizableDense); template static bool IsConverged( const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& state) requires(VectorizableDense); + const DenseMatrixPolicy& state, auto& stateParams) requires(VectorizableDense); }; } // namespace micm diff --git a/include/micm/solver/backward_euler.inl b/include/micm/solver/backward_euler.inl index ec7e4b9aa..a46d2ba1b 100644 --- a/include/micm/solver/backward_euler.inl +++ b/include/micm/solver/backward_euler.inl @@ -135,7 +135,7 @@ namespace micm continue; // check for convergence - converged = IsConverged(parameters_, forcing, Yn1); + converged = IsConverged(parameters_, forcing, Yn1, state); } while (!converged && iterations < max_iter); if (!converged) @@ -182,11 +182,11 @@ namespace micm inline bool BackwardEuler::IsConverged( const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& state) requires(!VectorizableDense) + const DenseMatrixPolicy& state, auto& stateParams) requires(!VectorizableDense) { double small = parameters.small_; - double rel_tol = parameters.relative_tolerance_; - auto& abs_tol = parameters.absolute_tolerance_; + double rel_tol = stateParams.relative_tolerance_; + auto& abs_tol = stateParams.absolute_tolerance_; auto residual_iter = residual.AsVector().begin(); auto state_iter = state.AsVector().begin(); const std::size_t n_elem = residual.NumRows() * residual.NumColumns(); @@ -208,11 +208,11 @@ namespace micm inline bool BackwardEuler::IsConverged( const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& state) requires(VectorizableDense) + const DenseMatrixPolicy& state, auto& stateParams) requires(VectorizableDense) { double small = parameters.small_; - double rel_tol = parameters.relative_tolerance_; - auto& abs_tol = parameters.absolute_tolerance_; + double rel_tol = stateParams.relative_tolerance_; + auto& abs_tol = stateParams.absolute_tolerance_; auto residual_iter = residual.AsVector().begin(); auto state_iter = state.AsVector().begin(); const std::size_t n_elem = residual.NumRows() * residual.NumColumns(); diff --git a/include/micm/solver/backward_euler_solver_parameters.hpp b/include/micm/solver/backward_euler_solver_parameters.hpp index 59b6ea03e..f9706e9b7 100644 --- a/include/micm/solver/backward_euler_solver_parameters.hpp +++ b/include/micm/solver/backward_euler_solver_parameters.hpp @@ -19,8 +19,7 @@ namespace micm template using SolverType = BackwardEuler; - std::vector absolute_tolerance_; - double relative_tolerance_{ 1.0e-8 }; + //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 diff --git a/include/micm/solver/rosenbrock.hpp b/include/micm/solver/rosenbrock.hpp index 714dce978..d91dbce1d 100644 --- a/include/micm/solver/rosenbrock.hpp +++ b/include/micm/solver/rosenbrock.hpp @@ -123,10 +123,10 @@ namespace micm /// @param errors The computed errors /// @return template - 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); template - 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); }; // end of Abstract Rosenbrock Solver diff --git a/include/micm/solver/rosenbrock.inl b/include/micm/solver/rosenbrock.inl index 8f463cdc2..37782b6c7 100644 --- a/include/micm/solver/rosenbrock.inl +++ b/include/micm/solver/rosenbrock.inl @@ -122,7 +122,7 @@ namespace micm Yerror.Axpy(parameters_.e_[stage], K[stage]); // Compute the normalized error - auto error = static_cast(this)->NormalizedError(Y, Ynew, Yerror); + auto error = static_cast(this)->NormalizedError(Y, Ynew, Yerror, state); // New step size is bounded by FacMin <= Hnew/H <= FacMax double fac = std::min( @@ -267,7 +267,8 @@ namespace micm inline double AbstractRosenbrockSolver::NormalizedError( const DenseMatrixPolicy& Y, const DenseMatrixPolicy& Ynew, - const DenseMatrixPolicy& errors) const requires(!VectorizableDense) + const DenseMatrixPolicy& errors, + auto& state) const requires(!VectorizableDense) { // Solving Ordinary Differential Equations II, page 123 // https://link-springer-com.cuucar.idm.oclc.org/book/10.1007/978-3-642-05221-7 @@ -278,7 +279,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 = parameters_.absolute_tolerance_.size(); + const std::size_t n_vars = state.absolute_tolerance_.size(); double ymax = 0; double errors_over_scale = 0; @@ -288,7 +289,7 @@ namespace micm { ymax = std::max(std::abs(_y[i]), std::abs(_ynew[i])); errors_over_scale = - _errors[i] / (parameters_.absolute_tolerance_[i % n_vars] + parameters_.relative_tolerance_ * ymax); + _errors[i] / (state.absolute_tolerance_[i % n_vars] + state.relative_tolerance_ * ymax); error += errors_over_scale * errors_over_scale; } @@ -302,7 +303,8 @@ namespace micm inline double AbstractRosenbrockSolver::NormalizedError( const DenseMatrixPolicy& Y, const DenseMatrixPolicy& Ynew, - const DenseMatrixPolicy& errors) const requires(VectorizableDense) + const DenseMatrixPolicy& errors, + auto& state) const requires(VectorizableDense) { // Solving Ordinary Differential Equations II, page 123 // https://link-springer-com.cuucar.idm.oclc.org/book/10.1007/978-3-642-05221-7 @@ -314,7 +316,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 = parameters_.absolute_tolerance_.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(); @@ -325,8 +327,8 @@ namespace micm for (std::size_t i = 0; i < whole_blocks; ++i) { errors_over_scale = - *errors_iter / (parameters_.absolute_tolerance_[(i / L) % n_vars] + - parameters_.relative_tolerance_ * std::max(std::abs(*y_iter), std::abs(*ynew_iter))); + *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; ++ynew_iter; @@ -344,8 +346,8 @@ namespace micm { const std::size_t idx = y * L + x; errors_over_scale = errors_iter[idx] / - (parameters_.absolute_tolerance_[y] + - parameters_.relative_tolerance_ * std::max(std::abs(y_iter[idx]), std::abs(ynew_iter[idx]))); + (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/rosenbrock_solver_parameters.hpp b/include/micm/solver/rosenbrock_solver_parameters.hpp index 6ae0a6b6f..9f7338ec1 100644 --- a/include/micm/solver/rosenbrock_solver_parameters.hpp +++ b/include/micm/solver/rosenbrock_solver_parameters.hpp @@ -9,6 +9,8 @@ #include #include +//#include "../include/micm/solver/solver.hpp" + namespace micm { @@ -60,9 +62,6 @@ namespace micm // Gamma_i = \sum_j gamma_{i,j} std::array gamma_{}; - std::vector absolute_tolerance_; - double relative_tolerance_{ 1e-6 }; - bool check_singularity_{ false }; // Check for singular A matrix in linear solve of A x = b // Print RosenbrockSolverParameters to console @@ -132,12 +131,6 @@ namespace micm std::cout << val << " "; std::cout << std::endl; std::cout << "absolute_tolerance: "; - for (auto& val : absolute_tolerance_) - { - std::cout << val << " "; - } - std::cout << std::endl; - std::cout << "relative_tolerance: " << relative_tolerance_ << std::endl; } inline RosenbrockSolverParameters RosenbrockSolverParameters::TwoStageRosenbrockParameters() diff --git a/include/micm/solver/solver.hpp b/include/micm/solver/solver.hpp index a387866c0..72b975073 100644 --- a/include/micm/solver/solver.hpp +++ b/include/micm/solver/solver.hpp @@ -64,9 +64,9 @@ namespace micm } // Overloaded Solve function to change parameters - SolverResult Solve(double time_step, StatePolicy& state, RosenbrockSolverParameters params) + SolverResult Solve(double time_step, StatePolicy& state, SolverParametersType& params) { - solver_.parameters_.h_start_ = params.h_start_; + solver_.parameters_ = params; return solver_.Solve(time_step, state); } @@ -88,7 +88,7 @@ namespace micm { return state_parameters_.number_of_rate_constants_; } - + StatePolicy GetState() const { auto state = std::move(StatePolicy(state_parameters_)); diff --git a/include/micm/solver/solver_builder.inl b/include/micm/solver/solver_builder.inl index ee188134e..a3889e56b 100644 --- a/include/micm/solver/solver_builder.inl +++ b/include/micm/solver/solver_builder.inl @@ -379,8 +379,7 @@ namespace micm } this->UnusedSpeciesCheck(); - this->SetAbsoluteTolerances(options.absolute_tolerance_, species_map); - + RatesPolicy rates(this->reactions_, species_map); auto nonzero_elements = rates.NonZeroJacobianElements(); auto jacobian = BuildJacobian(nonzero_elements, this->number_of_grid_cells_, number_of_species); @@ -398,6 +397,9 @@ namespace micm .variable_names_ = variable_names, .custom_rate_parameter_labels_ = labels, .nonzero_jacobian_elements_ = nonzero_elements }; + + + this->SetAbsoluteTolerances(state_parameters.absolute_tolerance_, species_map); //WHERE TO ADD THIS?????? return Solver( SolverPolicy(options, std::move(linear_solver), std::move(rates), jacobian), diff --git a/include/micm/solver/state.hpp b/include/micm/solver/state.hpp index 29bd50d37..6e1773d63 100644 --- a/include/micm/solver/state.hpp +++ b/include/micm/solver/state.hpp @@ -31,6 +31,8 @@ namespace micm std::vector variable_names_{}; std::vector custom_rate_parameter_labels_{}; std::set> nonzero_jacobian_elements_{}; + double relative_tolerance_; + std::vector absolute_tolerance_ {}; }; template @@ -58,6 +60,8 @@ namespace micm std::size_t state_size_; std::size_t number_of_grid_cells_; std::unique_ptr temporary_variables_; + double relative_tolerance_; + std::vector absolute_tolerance_; /// @brief Copy constructor /// @param other The state object to be copied @@ -76,6 +80,8 @@ namespace micm state_size_ = other.state_size_; number_of_grid_cells_ = other.number_of_grid_cells_; temporary_variables_ = std::make_unique(*other.temporary_variables_); + relative_tolerance_ = other.relative_tolerance_; + absolute_tolerance_ = other.absolute_tolerance_; } /// @brief Assignment operator @@ -98,6 +104,8 @@ namespace micm state_size_ = other.state_size_; number_of_grid_cells_ = other.number_of_grid_cells_; temporary_variables_ = std::make_unique(*other.temporary_variables_); + relative_tolerance_ = other.relative_tolerance_; + absolute_tolerance_ = other.absolute_tolerance_; } return *this; } diff --git a/include/micm/solver/state.inl b/include/micm/solver/state.inl index 20a3d725b..654befaed 100644 --- a/include/micm/solver/state.inl +++ b/include/micm/solver/state.inl @@ -80,7 +80,9 @@ namespace micm lower_matrix_(), upper_matrix_(), state_size_(parameters.variable_names_.size()), - number_of_grid_cells_(parameters.number_of_grid_cells_) + number_of_grid_cells_(parameters.number_of_grid_cells_), + relative_tolerance_(1e-06), + absolute_tolerance_(parameters.absolute_tolerance_) { std::size_t index = 0; for (auto& name : variable_names_) diff --git a/test/integration/analytical_policy.hpp b/test/integration/analytical_policy.hpp index 67b154b22..0ff45068e 100644 --- a/test/integration/analytical_policy.hpp +++ b/test/integration/analytical_policy.hpp @@ -1375,6 +1375,8 @@ void test_analytical_robertson( double air_density = 1e6; auto state = solver.GetState(); + state.relative_tolerance_ = 1e-10; + state.absolute_tolerance_ = std::vector(3, state.relative_tolerance_ * 1e-2); double k1 = 0.04; double k2 = 3e7; @@ -1583,6 +1585,9 @@ void test_analytical_oregonator( auto state = solver.GetState(); + state.relative_tolerance_ = 1e-6; + state.absolute_tolerance_ = std::vector(5, state.relative_tolerance_ * 1e-2); + state.SetCustomRateParameter("r1", 1.34 * 0.06); state.SetCustomRateParameter("r2", 1.6e9); state.SetCustomRateParameter("r3", 8e3 * 0.06); @@ -1758,6 +1763,8 @@ void test_analytical_hires( }; auto state = solver.GetState(); + state.relative_tolerance_ = 1e-6; + state.absolute_tolerance_ = std::vector(8, state.relative_tolerance_ * 1e-2); state.SetCustomRateParameter("r1", 1.71); state.SetCustomRateParameter("r2", 8.75); @@ -1916,6 +1923,12 @@ void test_analytical_e5( auto state = solver.GetState(); + state.relative_tolerance_ = 1e-13; + 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); state.SetCustomRateParameter("r3", 1.1e7); diff --git a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp index d2543561c..d109ff8ca 100644 --- a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp +++ b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp @@ -180,14 +180,14 @@ TEST(AnalyticalExamplesCudaRosenbrock, E5) { auto rosenbrock_solver = [](auto params) { - params.relative_tolerance_ = 1e-13; + /*params.relative_tolerance_ = 1e-13; params.absolute_tolerance_ = std::vector(6, 1e-17); // this paper https://archimede.uniba.it/~testset/report/e5.pdf // says that the first variable should have a much looser tolerance than the other species params.absolute_tolerance_[0] = 1e-7; // these last two aren't actually provided values and we don't care how they behave params.absolute_tolerance_[4] = 1e-7; - params.absolute_tolerance_[5] = 1e-7; + params.absolute_tolerance_[5] = 1e-7;*/ return builderType1Cell(params); }; @@ -212,8 +212,8 @@ TEST(AnalyticalExamplesCudaRosenbrock, Oregonator) auto rosenbrock_solver = [](auto params) { // anything below 1e-6 is too strict for the Oregonator - params.relative_tolerance_ = 1e-6; - params.absolute_tolerance_ = std::vector(5, params.relative_tolerance_ * 1e-2); + //params.relative_tolerance_ = 1e-6; + //params.absolute_tolerance_ = std::vector(5, params.relative_tolerance_ * 1e-2); return builderType1Cell(params); }; @@ -237,8 +237,8 @@ TEST(AnalyticalExamplesCudaRosenbrock, HIRES) { auto rosenbrock_solver = [](auto params) { - params.relative_tolerance_ = 1e-6; - params.absolute_tolerance_ = std::vector(8, params.relative_tolerance_ * 1e-2); + //params.relative_tolerance_ = 1e-6; + //params.absolute_tolerance_ = std::vector(8, params.relative_tolerance_ * 1e-2); return builderType1Cell(params); }; diff --git a/test/integration/jit/test_jit_analytical_rosenbrock.cpp b/test/integration/jit/test_jit_analytical_rosenbrock.cpp index e9e91ac04..18dbc891e 100644 --- a/test/integration/jit/test_jit_analytical_rosenbrock.cpp +++ b/test/integration/jit/test_jit_analytical_rosenbrock.cpp @@ -183,14 +183,14 @@ TEST(AnalyticalExamplesJitRosenbrock, E5) { auto rosenbrock_solver = [](auto params) { - params.relative_tolerance_ = 1e-13; + /*params.relative_tolerance_ = 1e-13; params.absolute_tolerance_ = std::vector(6, 1e-17); // this paper https://archimede.uniba.it/~testset/report/e5.pdf // says that the first variable should have a much looser tolerance than the other species params.absolute_tolerance_[0] = 1e-7; // these last two aren't actually provided values and we don't care how they behave params.absolute_tolerance_[4] = 1e-7; - params.absolute_tolerance_[5] = 1e-7; + params.absolute_tolerance_[5] = 1e-7;*/ return BuilderType<1>(params); }; @@ -215,8 +215,8 @@ TEST(AnalyticalExamples, Oregonator) auto rosenbrock_solver = [](auto params) { // anything below 1e-6 is too strict for the Oregonator - params.relative_tolerance_ = 1e-6; - params.absolute_tolerance_ = std::vector(5, params.relative_tolerance_ * 1e-2); + //params.relative_tolerance_ = 1e-6; + //params.absolute_tolerance_ = std::vector(5, params.relative_tolerance_ * 1e-2); return BuilderType<1>(params); }; @@ -240,8 +240,8 @@ TEST(AnalyticalExamples, HIRES) { auto rosenbrock_solver = [](auto params) { - params.relative_tolerance_ = 1e-6; - params.absolute_tolerance_ = std::vector(8, params.relative_tolerance_ * 1e-2); + //params.relative_tolerance_ = 1e-6; + //params.absolute_tolerance_ = std::vector(8, params.relative_tolerance_ * 1e-2); return BuilderType<1>(params); }; diff --git a/test/integration/terminator.hpp b/test/integration/terminator.hpp index e78e77ba9..fe9832c6b 100644 --- a/test/integration/terminator.hpp +++ b/test/integration/terminator.hpp @@ -48,6 +48,7 @@ void TestTerminator(BuilderPolicy& builder, std::size_t number_of_grid_cells) .SetNumberOfGridCells(number_of_grid_cells) .Build(); auto state = solver.GetState(); + state.relative_tolerance_ = 1.0e-8; auto get_double = std::bind(std::lognormal_distribution(-2.0, 2.0), std::default_random_engine()); std::unordered_map> concentrations{ { "Cl2", {} }, { "Cl", {} } }; diff --git a/test/integration/test_analytical_rosenbrock.cpp b/test/integration/test_analytical_rosenbrock.cpp index 3f6211706..127b537f7 100644 --- a/test/integration/test_analytical_rosenbrock.cpp +++ b/test/integration/test_analytical_rosenbrock.cpp @@ -9,7 +9,7 @@ #include #include - +/* using BuilderType = micm::CpuSolverBuilder; using StateType = micm::State; @@ -198,8 +198,8 @@ TEST(AnalyticalExamples, Robertson) { auto rosenbrock_solver = [](auto params) { - params.relative_tolerance_ = 1e-10; - params.absolute_tolerance_ = std::vector(3, params.relative_tolerance_ * 1e-2); + // params.relative_tolerance_ = 1e-10; + // params.absolute_tolerance_ = std::vector(3, params.relative_tolerance_ * 1e-2); return BuilderType(params); }; @@ -223,14 +223,14 @@ TEST(AnalyticalExamples, E5) { auto rosenbrock_solver = [](auto params) { - params.relative_tolerance_ = 1e-13; - params.absolute_tolerance_ = std::vector(6, 1e-17); + //params.relative_tolerance_ = 1e-13; + //params.absolute_tolerance_ = std::vector(6, 1e-17); // this paper https://archimede.uniba.it/~testset/report/e5.pdf // says that the first variable should have a much looser tolerance than the other species - params.absolute_tolerance_[0] = 1e-7; + //params.absolute_tolerance_[0] = 1e-7; // these last two aren't actually provided values and we don't care how they behave - params.absolute_tolerance_[4] = 1e-7; - params.absolute_tolerance_[5] = 1e-7; + //params.absolute_tolerance_[4] = 1e-7; + //params.absolute_tolerance_[5] = 1e-7; return BuilderType(params); }; @@ -255,8 +255,8 @@ TEST(AnalyticalExamples, Oregonator) auto rosenbrock_solver = [](auto params) { // anything below 1e-6 is too strict for the Oregonator - params.relative_tolerance_ = 1e-6; - params.absolute_tolerance_ = std::vector(5, params.relative_tolerance_ * 1e-2); + //params.relative_tolerance_ = 1e-6; + //params.absolute_tolerance_ = std::vector(5, params.relative_tolerance_ * 1e-2); return micm::CpuSolverBuilder(params); }; @@ -289,8 +289,8 @@ TEST(AnalyticalExamples, HIRES) { auto rosenbrock_solver = [](auto params) { - params.relative_tolerance_ = 1e-6; - params.absolute_tolerance_ = std::vector(8, params.relative_tolerance_ * 1e-2); + //params.relative_tolerance_ = 1e-6; + //params.absolute_tolerance_ = std::vector(8, params.relative_tolerance_ * 1e-2); return micm::CpuSolverBuilder(params); }; @@ -308,4 +308,4 @@ TEST(AnalyticalExamples, HIRES) solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_hires(solver, 1e-3); -} +}*/ diff --git a/test/integration/test_chapman_integration.cpp b/test/integration/test_chapman_integration.cpp index 0793e2d12..5cea1e877 100644 --- a/test/integration/test_chapman_integration.cpp +++ b/test/integration/test_chapman_integration.cpp @@ -32,9 +32,7 @@ TEST(ChapmanIntegration, CanBuildChapmanSystemUsingConfig) // Get solver parameters ('System', the collection of 'Process') micm::SolverParameters solver_params = solverConfig.GetSolverParams(); - auto options = solver_params.parameters_; - - EXPECT_EQ(options.relative_tolerance_, 1.0e-4); + auto options = solver_params.parameters_; auto solver = micm::CpuSolverBuilder(options) .SetSystem(solver_params.system_) @@ -44,16 +42,18 @@ TEST(ChapmanIntegration, CanBuildChapmanSystemUsingConfig) micm::State state = solver.GetState(); + EXPECT_EQ(state.relative_tolerance_, 1.0e-6); + for (size_t n_grid_cell = 0; n_grid_cell < state.number_of_grid_cells_; ++n_grid_cell) { - EXPECT_EQ(solver.solver_.parameters_.absolute_tolerance_[state.variable_map_["Ar"]], 1.0e-12); - EXPECT_EQ(solver.solver_.parameters_.absolute_tolerance_[state.variable_map_["CO2"]], 1.0e-12); - EXPECT_EQ(solver.solver_.parameters_.absolute_tolerance_[state.variable_map_["H2O"]], 1.0e-12); - EXPECT_EQ(solver.solver_.parameters_.absolute_tolerance_[state.variable_map_["N2"]], 1.0e-12); - EXPECT_EQ(solver.solver_.parameters_.absolute_tolerance_[state.variable_map_["O1D"]], 1.0e-12); - EXPECT_EQ(solver.solver_.parameters_.absolute_tolerance_[state.variable_map_["O"]], 1.0e-12); - EXPECT_EQ(solver.solver_.parameters_.absolute_tolerance_[state.variable_map_["O2"]], 1.0e-12); - EXPECT_EQ(solver.solver_.parameters_.absolute_tolerance_[state.variable_map_["O3"]], 1.0e-12); + EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["Ar"]], 1.0e-12); + EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["CO2"]], 1.0e-12); + EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["H2O"]], 1.0e-12); + EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["N2"]], 1.0e-12); + EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["O1D"]], 1.0e-12); + EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["O"]], 1.0e-12); + EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["O2"]], 1.0e-12); + EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["O3"]], 1.0e-12); } // User gives an input of concentrations diff --git a/test/integration/test_terminator.cpp b/test/integration/test_terminator.cpp index b86d35cb5..7fc947e61 100644 --- a/test/integration/test_terminator.cpp +++ b/test/integration/test_terminator.cpp @@ -14,7 +14,7 @@ TEST(RosenbrockSolver, Terminator) { auto parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - parameters.relative_tolerance_ = 1.0e-8; + //parameters.relative_tolerance_ = 1.0e-8; parameters.max_number_of_steps_ = 100000; { auto builder = micm::CpuSolverBuilder(parameters).SetIgnoreUnusedSpecies(true); @@ -60,7 +60,7 @@ using VectorBuilder = micm::CpuSolverBuilder< TEST(RosenbrockSolver, VectorTerminator) { auto parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - parameters.relative_tolerance_ = 1.0e-8; + //parameters.relative_tolerance_ = 1.0e-8; parameters.max_number_of_steps_ = 100000; { auto builder = VectorBuilder<1>(parameters).SetIgnoreUnusedSpecies(true); diff --git a/test/unit/cuda/solver/test_cuda_rosenbrock.cpp b/test/unit/cuda/solver/test_cuda_rosenbrock.cpp index 825d111c1..3d5d8791b 100644 --- a/test/unit/cuda/solver/test_cuda_rosenbrock.cpp +++ b/test/unit/cuda/solver/test_cuda_rosenbrock.cpp @@ -121,11 +121,10 @@ void testNormalizedErrorConst() auto gpu_builder = GpuBuilder(micm::CudaRosenbrockSolverParameters::ThreeStageRosenbrockParameters()); 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_; + double rtol = state.relative_tolerance_; - std::vector atol = gpu_solver.solver_.parameters_.absolute_tolerance_; - double rtol = gpu_solver.solver_.parameters_.relative_tolerance_; - - auto state = gpu_solver.GetState(); auto y_old = micm::CudaDenseMatrix(number_of_grid_cells, state.state_size_, 1.0); auto y_new = micm::CudaDenseMatrix(number_of_grid_cells, state.state_size_, 2.0); auto errors = micm::CudaDenseMatrix(number_of_grid_cells, state.state_size_, 3.0); @@ -134,7 +133,7 @@ void testNormalizedErrorConst() y_new.CopyToDevice(); errors.CopyToDevice(); - double error = gpu_solver.solver_.NormalizedError(y_old, y_new, errors); + double error = gpu_solver.solver_.NormalizedError(y_old, y_new, errors, state); auto expected_error = 0.0; for (size_t i = 0; i < state.state_size_; ++i) @@ -166,11 +165,9 @@ void testNormalizedErrorDiff() auto gpu_builder = GpuBuilder(micm::CudaRosenbrockSolverParameters::ThreeStageRosenbrockParameters()); gpu_builder = getSolver(gpu_builder); auto gpu_solver = gpu_builder.SetNumberOfGridCells(number_of_grid_cells).Build(); - - std::vector atol = gpu_solver.solver_.parameters_.absolute_tolerance_; - double rtol = gpu_solver.solver_.parameters_.relative_tolerance_; - auto state = gpu_solver.GetState(); + std::vector 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); @@ -195,7 +192,7 @@ void testNormalizedErrorDiff() y_new.CopyToDevice(); errors.CopyToDevice(); - double computed_error = gpu_solver.solver_.NormalizedError(y_old, y_new, errors); + double computed_error = gpu_solver.solver_.NormalizedError(y_old, y_new, errors, state); auto relative_error = std::abs(computed_error - expected_error) / std::max(std::abs(computed_error), std::abs(expected_error)); diff --git a/test/unit/solver/test_backward_euler.cpp b/test/unit/solver/test_backward_euler.cpp index df10aa3f3..36c062f65 100644 --- a/test/unit/solver/test_backward_euler.cpp +++ b/test/unit/solver/test_backward_euler.cpp @@ -39,7 +39,6 @@ namespace TEST(BackwardEuler, CanCallSolve) { auto params = micm::BackwardEulerSolverParameters(); - params.absolute_tolerance_ = { 1e-6, 1e-6, 1e-6 }; auto be = micm::CpuSolverBuilder(params) .SetSystem(the_system) @@ -49,6 +48,7 @@ TEST(BackwardEuler, CanCallSolve) double time_step = 1.0; auto state = be.GetState(); + state.absolute_tolerance_ = { 1e-6, 1e-6, 1e-6 }; state.variables_[0] = { 1.0, 0.0, 0.0 }; state.conditions_[0].temperature_ = 272.5; @@ -71,7 +71,7 @@ void CheckIsConverged() DenseMatrixPolicy state{ 4, 3, 0.0 }; parameters.small_ = 1e-6; - parameters.relative_tolerance_ = 1e-3; + /*parameters.relative_tolerance_ = 1e-3; parameters.absolute_tolerance_ = { 1e-6, 1e-6, 1e-6 }; ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state)); @@ -87,6 +87,7 @@ void CheckIsConverged() ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, state)); parameters.absolute_tolerance_[2] = 1.0; ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state)); + */ } TEST(BackwardEuler, IsConverged) diff --git a/test/unit/solver/test_rosenbrock.cpp b/test/unit/solver/test_rosenbrock.cpp index 568b8e6b9..5d3a992a7 100644 --- a/test/unit/solver/test_rosenbrock.cpp +++ b/test/unit/solver/test_rosenbrock.cpp @@ -16,11 +16,11 @@ template void testNormalizedErrorDiff(SolverBuilderPolicy builder, std::size_t number_of_grid_cells) { builder = getSolver(builder); - auto solver = builder.SetNumberOfGridCells(number_of_grid_cells).Build(); - std::vector atol = solver.solver_.parameters_.absolute_tolerance_; - double rtol = solver.solver_.parameters_.relative_tolerance_; - + auto solver = builder.SetNumberOfGridCells(number_of_grid_cells).Build(); auto state = solver.GetState(); + std::vector atol = state.absolute_tolerance_; + double rtol = state.relative_tolerance_; + using MatrixPolicy = decltype(state.variables_); auto y_old = MatrixPolicy(number_of_grid_cells, state.state_size_, 7.7); auto y_new = MatrixPolicy(number_of_grid_cells, state.state_size_, -13.9); @@ -42,7 +42,7 @@ void testNormalizedErrorDiff(SolverBuilderPolicy builder, std::size_t number_of_ double error_min_ = 1.0e-10; expected_error = std::max(std::sqrt(expected_error / (number_of_grid_cells * state.state_size_)), error_min_); - double computed_error = solver.solver_.NormalizedError(y_old, y_new, errors); + double computed_error = solver.solver_.NormalizedError(y_old, y_new, errors, state); auto relative_error = std::abs(computed_error - expected_error) / std::max(std::abs(computed_error), std::abs(expected_error)); @@ -106,9 +106,10 @@ TEST(RosenbrockSolver, CanSetTolerances) .SetReactions(std::vector{ r1 }) .SetNumberOfGridCells(number_of_grid_cells) .Build(); - EXPECT_EQ(solver.solver_.parameters_.absolute_tolerance_.size(), 2); - EXPECT_EQ(solver.solver_.parameters_.absolute_tolerance_[0], 1.0e-07); - EXPECT_EQ(solver.solver_.parameters_.absolute_tolerance_[1], 1.0e-08); + 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); } } diff --git a/test/unit/solver/test_solver_builder.cpp b/test/unit/solver/test_solver_builder.cpp index 28c53bb5f..5f045cd25 100644 --- a/test/unit/solver/test_solver_builder.cpp +++ b/test/unit/solver/test_solver_builder.cpp @@ -106,6 +106,7 @@ TEST(SolverBuilder, CanBuildRosenbrock) .Build(); } +/* TEST(SolverBuilder, MismatchedToleranceSizeIsCaught) { auto params = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); @@ -127,4 +128,5 @@ TEST(SolverBuilder, MismatchedToleranceSizeIsCaught) micm::SparseMatrix>>(params); EXPECT_ANY_THROW(builder.SetSystem(the_system).SetReactions(reactions).SetNumberOfGridCells(1).Build();); -} \ No newline at end of file +} +*/ \ No newline at end of file From 32a16e2e906e082dd8391d3836f44d4892fc481a Mon Sep 17 00:00:00 2001 From: Montek Thind Date: Tue, 1 Oct 2024 13:53:19 -0700 Subject: [PATCH 03/31] third draft --- include/micm/configure/solver_config.hpp | 15 ++++++++++++--- .../micm/solver/rosenbrock_solver_parameters.hpp | 2 -- include/micm/solver/solver_builder.hpp | 1 + include/micm/solver/solver_builder.inl | 4 ++-- include/micm/solver/state.hpp | 3 ++- include/micm/solver/state.inl | 6 ++++++ test/integration/test_chapman_integration.cpp | 10 +++++----- test/unit/solver/test_solver_builder.cpp | 4 +--- 8 files changed, 29 insertions(+), 16 deletions(-) diff --git a/include/micm/configure/solver_config.hpp b/include/micm/configure/solver_config.hpp index 0358e8cf6..27af08b7e 100644 --- a/include/micm/configure/solver_config.hpp +++ b/include/micm/configure/solver_config.hpp @@ -102,6 +102,7 @@ namespace micm System system_; std::vector processes_; RosenbrockSolverParameters parameters_; + double relative_tolerance_; SolverParameters(const System& system, std::vector&& processes, const RosenbrockSolverParameters&& parameters) : system_(system), @@ -116,6 +117,13 @@ namespace micm parameters_(parameters) { } + SolverParameters(System&& system, std::vector&& 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 phases_; std::vector 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_arr; @@ -426,7 +435,7 @@ namespace micm void ParseRelativeTolerance(const objectType& object) { ValidateSchema(object, { "value", "type" }, {}); - //this->parameters_.relative_tolerance_ = object["value"].as(); + this->relative_tolerance_ = object["value"].as(); } 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 diff --git a/include/micm/solver/rosenbrock_solver_parameters.hpp b/include/micm/solver/rosenbrock_solver_parameters.hpp index 9f7338ec1..8ae4c67b0 100644 --- a/include/micm/solver/rosenbrock_solver_parameters.hpp +++ b/include/micm/solver/rosenbrock_solver_parameters.hpp @@ -9,8 +9,6 @@ #include #include -//#include "../include/micm/solver/solver.hpp" - namespace micm { diff --git a/include/micm/solver/solver_builder.hpp b/include/micm/solver/solver_builder.hpp index c3a8d7a3a..fc5d741fe 100644 --- a/include/micm/solver/solver_builder.hpp +++ b/include/micm/solver/solver_builder.hpp @@ -50,6 +50,7 @@ namespace micm bool reorder_state_ = true; bool valid_system_ = false; bool valid_reactions_ = false; + double relative_tolerance_ = 1e-06; public: SolverBuilder() = delete; diff --git a/include/micm/solver/solver_builder.inl b/include/micm/solver/solver_builder.inl index a3889e56b..e98e43538 100644 --- a/include/micm/solver/solver_builder.inl +++ b/include/micm/solver/solver_builder.inl @@ -316,7 +316,7 @@ namespace micm } } } - } + } template< class SolverParametersPolicy, @@ -399,7 +399,7 @@ namespace micm .nonzero_jacobian_elements_ = nonzero_elements }; - this->SetAbsoluteTolerances(state_parameters.absolute_tolerance_, species_map); //WHERE TO ADD THIS?????? + this->SetAbsoluteTolerances(state_parameters.absolute_tolerance_, species_map); return Solver( SolverPolicy(options, std::move(linear_solver), std::move(rates), jacobian), diff --git a/include/micm/solver/state.hpp b/include/micm/solver/state.hpp index 6e1773d63..cb38eeddc 100644 --- a/include/micm/solver/state.hpp +++ b/include/micm/solver/state.hpp @@ -31,7 +31,7 @@ namespace micm std::vector variable_names_{}; std::vector custom_rate_parameter_labels_{}; std::set> nonzero_jacobian_elements_{}; - double relative_tolerance_; + double relative_tolerance_{ 1e-06 }; std::vector absolute_tolerance_ {}; }; @@ -140,6 +140,7 @@ namespace micm /// @param value new parameter value void SetCustomRateParameter(const std::string& label, double value); void SetCustomRateParameter(const std::string& label, const std::vector& values); + void SetRelativeTolerances(double relativeTolerance); /// @brief Print a header of species to display concentrations with respect to time void PrintHeader(); diff --git a/include/micm/solver/state.inl b/include/micm/solver/state.inl index 654befaed..37f3962b5 100644 --- a/include/micm/solver/state.inl +++ b/include/micm/solver/state.inl @@ -187,6 +187,12 @@ namespace micm for (std::size_t i = 0; i < custom_rate_parameters_.NumRows(); ++i) custom_rate_parameters_[i][param->second] = values[i]; } + + template + inline void State::SetRelativeTolerances(double relativeTolerance) + { + this->relative_tolerance_ = relativeTolerance; + } template inline void State::PrintHeader() diff --git a/test/integration/test_chapman_integration.cpp b/test/integration/test_chapman_integration.cpp index 5cea1e877..8ebdc001f 100644 --- a/test/integration/test_chapman_integration.cpp +++ b/test/integration/test_chapman_integration.cpp @@ -32,17 +32,17 @@ TEST(ChapmanIntegration, CanBuildChapmanSystemUsingConfig) // Get solver parameters ('System', the collection of 'Process') micm::SolverParameters solver_params = solverConfig.GetSolverParams(); - auto options = solver_params.parameters_; + auto options = solver_params.parameters_; auto solver = micm::CpuSolverBuilder(options) .SetSystem(solver_params.system_) .SetReactions(solver_params.processes_) - .SetIgnoreUnusedSpecies(true) - .Build(); + .SetIgnoreUnusedSpecies(true) + .Build(); micm::State state = solver.GetState(); - - EXPECT_EQ(state.relative_tolerance_, 1.0e-6); + state.SetRelativeTolerances(solver_params.relative_tolerance_); + EXPECT_EQ(state.relative_tolerance_, 1.0e-4); for (size_t n_grid_cell = 0; n_grid_cell < state.number_of_grid_cells_; ++n_grid_cell) { diff --git a/test/unit/solver/test_solver_builder.cpp b/test/unit/solver/test_solver_builder.cpp index 5f045cd25..fc9209ef2 100644 --- a/test/unit/solver/test_solver_builder.cpp +++ b/test/unit/solver/test_solver_builder.cpp @@ -105,7 +105,6 @@ TEST(SolverBuilder, CanBuildRosenbrock) .SetNumberOfGridCells(1) .Build(); } - /* TEST(SolverBuilder, MismatchedToleranceSizeIsCaught) { @@ -128,5 +127,4 @@ TEST(SolverBuilder, MismatchedToleranceSizeIsCaught) micm::SparseMatrix>>(params); EXPECT_ANY_THROW(builder.SetSystem(the_system).SetReactions(reactions).SetNumberOfGridCells(1).Build();); -} -*/ \ No newline at end of file +}*/ \ No newline at end of file From f125b1615bceea89f12379ada77b0e4049f9ec8f Mon Sep 17 00:00:00 2001 From: Montek Thind Date: Tue, 1 Oct 2024 16:20:55 -0700 Subject: [PATCH 04/31] darft fourth --- examples/profile_example.cpp | 5 +++-- include/micm/solver/backward_euler.inl | 1 - include/micm/version.hpp | 2 -- .../cuda/test_cuda_analytical_rosenbrock.cpp | 13 ------------- .../jit/test_jit_analytical_rosenbrock.cpp | 13 ------------- .../test_analytical_rosenbrock.cpp | 19 ++----------------- test/integration/test_terminator.cpp | 2 -- test/unit/solver/test_backward_euler.cpp | 4 ++-- 8 files changed, 7 insertions(+), 52 deletions(-) diff --git a/examples/profile_example.cpp b/examples/profile_example.cpp index cb206e0d2..852e0f9fd 100644 --- 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(params) .SetSystem(chemical_system) .SetReactions(reactions) .Build(); auto state = solver.GetState(); - + state.SetRelativeTolerances(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/solver/backward_euler.inl b/include/micm/solver/backward_euler.inl index a46d2ba1b..b5805bdea 100644 --- a/include/micm/solver/backward_euler.inl +++ b/include/micm/solver/backward_euler.inl @@ -219,7 +219,6 @@ namespace micm 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) { diff --git a/include/micm/version.hpp b/include/micm/version.hpp index c0ae8410d..fd5398ebf 100644 --- a/include/micm/version.hpp +++ b/include/micm/version.hpp @@ -1,5 +1,3 @@ -// Copyright (C) 2023-2024 National Center for Atmospheric Research -// SPDX-License-Identifier: Apache-2.0 // clang-format off #pragma once diff --git a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp index d109ff8ca..d1551e98c 100644 --- a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp +++ b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp @@ -180,14 +180,6 @@ TEST(AnalyticalExamplesCudaRosenbrock, E5) { auto rosenbrock_solver = [](auto params) { - /*params.relative_tolerance_ = 1e-13; - params.absolute_tolerance_ = std::vector(6, 1e-17); - // this paper https://archimede.uniba.it/~testset/report/e5.pdf - // says that the first variable should have a much looser tolerance than the other species - params.absolute_tolerance_[0] = 1e-7; - // these last two aren't actually provided values and we don't care how they behave - params.absolute_tolerance_[4] = 1e-7; - params.absolute_tolerance_[5] = 1e-7;*/ return builderType1Cell(params); }; @@ -211,9 +203,6 @@ TEST(AnalyticalExamplesCudaRosenbrock, Oregonator) { auto rosenbrock_solver = [](auto params) { - // anything below 1e-6 is too strict for the Oregonator - //params.relative_tolerance_ = 1e-6; - //params.absolute_tolerance_ = std::vector(5, params.relative_tolerance_ * 1e-2); return builderType1Cell(params); }; @@ -237,8 +226,6 @@ TEST(AnalyticalExamplesCudaRosenbrock, HIRES) { auto rosenbrock_solver = [](auto params) { - //params.relative_tolerance_ = 1e-6; - //params.absolute_tolerance_ = std::vector(8, params.relative_tolerance_ * 1e-2); return builderType1Cell(params); }; diff --git a/test/integration/jit/test_jit_analytical_rosenbrock.cpp b/test/integration/jit/test_jit_analytical_rosenbrock.cpp index 18dbc891e..1688c0118 100644 --- a/test/integration/jit/test_jit_analytical_rosenbrock.cpp +++ b/test/integration/jit/test_jit_analytical_rosenbrock.cpp @@ -183,14 +183,6 @@ TEST(AnalyticalExamplesJitRosenbrock, E5) { auto rosenbrock_solver = [](auto params) { - /*params.relative_tolerance_ = 1e-13; - params.absolute_tolerance_ = std::vector(6, 1e-17); - // this paper https://archimede.uniba.it/~testset/report/e5.pdf - // says that the first variable should have a much looser tolerance than the other species - params.absolute_tolerance_[0] = 1e-7; - // these last two aren't actually provided values and we don't care how they behave - params.absolute_tolerance_[4] = 1e-7; - params.absolute_tolerance_[5] = 1e-7;*/ return BuilderType<1>(params); }; @@ -214,9 +206,6 @@ TEST(AnalyticalExamples, Oregonator) { auto rosenbrock_solver = [](auto params) { - // anything below 1e-6 is too strict for the Oregonator - //params.relative_tolerance_ = 1e-6; - //params.absolute_tolerance_ = std::vector(5, params.relative_tolerance_ * 1e-2); return BuilderType<1>(params); }; @@ -240,8 +229,6 @@ TEST(AnalyticalExamples, HIRES) { auto rosenbrock_solver = [](auto params) { - //params.relative_tolerance_ = 1e-6; - //params.absolute_tolerance_ = std::vector(8, params.relative_tolerance_ * 1e-2); return BuilderType<1>(params); }; diff --git a/test/integration/test_analytical_rosenbrock.cpp b/test/integration/test_analytical_rosenbrock.cpp index 127b537f7..3a9d54792 100644 --- a/test/integration/test_analytical_rosenbrock.cpp +++ b/test/integration/test_analytical_rosenbrock.cpp @@ -9,7 +9,7 @@ #include #include -/* + using BuilderType = micm::CpuSolverBuilder; using StateType = micm::State; @@ -198,8 +198,6 @@ TEST(AnalyticalExamples, Robertson) { auto rosenbrock_solver = [](auto params) { - // params.relative_tolerance_ = 1e-10; - // params.absolute_tolerance_ = std::vector(3, params.relative_tolerance_ * 1e-2); return BuilderType(params); }; @@ -223,14 +221,6 @@ TEST(AnalyticalExamples, E5) { auto rosenbrock_solver = [](auto params) { - //params.relative_tolerance_ = 1e-13; - //params.absolute_tolerance_ = std::vector(6, 1e-17); - // this paper https://archimede.uniba.it/~testset/report/e5.pdf - // says that the first variable should have a much looser tolerance than the other species - //params.absolute_tolerance_[0] = 1e-7; - // these last two aren't actually provided values and we don't care how they behave - //params.absolute_tolerance_[4] = 1e-7; - //params.absolute_tolerance_[5] = 1e-7; return BuilderType(params); }; @@ -254,9 +244,6 @@ TEST(AnalyticalExamples, Oregonator) { auto rosenbrock_solver = [](auto params) { - // anything below 1e-6 is too strict for the Oregonator - //params.relative_tolerance_ = 1e-6; - //params.absolute_tolerance_ = std::vector(5, params.relative_tolerance_ * 1e-2); return micm::CpuSolverBuilder(params); }; @@ -289,8 +276,6 @@ TEST(AnalyticalExamples, HIRES) { auto rosenbrock_solver = [](auto params) { - //params.relative_tolerance_ = 1e-6; - //params.absolute_tolerance_ = std::vector(8, params.relative_tolerance_ * 1e-2); return micm::CpuSolverBuilder(params); }; @@ -308,4 +293,4 @@ TEST(AnalyticalExamples, HIRES) solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_hires(solver, 1e-3); -}*/ +} diff --git a/test/integration/test_terminator.cpp b/test/integration/test_terminator.cpp index 7fc947e61..804d12972 100644 --- a/test/integration/test_terminator.cpp +++ b/test/integration/test_terminator.cpp @@ -14,7 +14,6 @@ TEST(RosenbrockSolver, Terminator) { auto parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - //parameters.relative_tolerance_ = 1.0e-8; parameters.max_number_of_steps_ = 100000; { auto builder = micm::CpuSolverBuilder(parameters).SetIgnoreUnusedSpecies(true); @@ -60,7 +59,6 @@ using VectorBuilder = micm::CpuSolverBuilder< TEST(RosenbrockSolver, VectorTerminator) { auto parameters = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - //parameters.relative_tolerance_ = 1.0e-8; parameters.max_number_of_steps_ = 100000; { auto builder = VectorBuilder<1>(parameters).SetIgnoreUnusedSpecies(true); diff --git a/test/unit/solver/test_backward_euler.cpp b/test/unit/solver/test_backward_euler.cpp index 36c062f65..4be954f4c 100644 --- a/test/unit/solver/test_backward_euler.cpp +++ b/test/unit/solver/test_backward_euler.cpp @@ -86,8 +86,8 @@ void CheckIsConverged() residual[3][2] = 1e-1; ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, state)); parameters.absolute_tolerance_[2] = 1.0; - ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state)); - */ + ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state));*/ + } TEST(BackwardEuler, IsConverged) From 03a8105611a0a88816cb93f55b5e9e0fa6692a9d Mon Sep 17 00:00:00 2001 From: Montek Thind Date: Wed, 2 Oct 2024 08:58:33 -0700 Subject: [PATCH 05/31] test fixed --- test/unit/solver/test_backward_euler.cpp | 22 ++++++++++----------- test/unit/solver/test_solver_builder.cpp | 25 +----------------------- 2 files changed, 12 insertions(+), 35 deletions(-) diff --git a/test/unit/solver/test_backward_euler.cpp b/test/unit/solver/test_backward_euler.cpp index 4be954f4c..b48141585 100644 --- a/test/unit/solver/test_backward_euler.cpp +++ b/test/unit/solver/test_backward_euler.cpp @@ -69,25 +69,25 @@ void CheckIsConverged() micm::BackwardEulerSolverParameters parameters; DenseMatrixPolicy residual{ 4, 3, 0.0 }; DenseMatrixPolicy state{ 4, 3, 0.0 }; + micm::StateParameters state_params; parameters.small_ = 1e-6; - /*parameters.relative_tolerance_ = 1e-3; - parameters.absolute_tolerance_ = { 1e-6, 1e-6, 1e-6 }; + state_params.relative_tolerance_ = 1e-3; + state_params.absolute_tolerance_ = { 1e-6, 1e-6, 1e-6 }; - ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state)); + ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state, state_params)); residual[0][1] = 1e-5; - ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, state)); + ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, state, state_params)); parameters.small_ = 1e-4; - ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state)); + ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state, state_params)); residual[3][2] = 1e-3; - ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, state)); + ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, state, state_params)); state[3][2] = 10.0; - ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state)); + ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state, state_params)); residual[3][2] = 1e-1; - ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, state)); - parameters.absolute_tolerance_[2] = 1.0; - ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state));*/ - + ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, state, state_params)); + state_params.absolute_tolerance_[2] = 1.0; + ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state, state_params)); } TEST(BackwardEuler, IsConverged) diff --git a/test/unit/solver/test_solver_builder.cpp b/test/unit/solver/test_solver_builder.cpp index fc9209ef2..047c0972c 100644 --- a/test/unit/solver/test_solver_builder.cpp +++ b/test/unit/solver/test_solver_builder.cpp @@ -104,27 +104,4 @@ TEST(SolverBuilder, CanBuildRosenbrock) .SetReactions(reactions) .SetNumberOfGridCells(1) .Build(); -} -/* -TEST(SolverBuilder, MismatchedToleranceSizeIsCaught) -{ - auto params = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - // too many - params.absolute_tolerance_ = { 1e-6, 1e-6, 1e-6, 1e-6, 1e-6 }; - EXPECT_ANY_THROW(micm::CpuSolverBuilder(params) - .SetSystem(the_system) - .SetReactions(reactions) - .SetNumberOfGridCells(1) - .Build();); - - constexpr std::size_t L = 4; - // too few - params.absolute_tolerance_ = { 1e-6, 1e-6 }; - - auto builder = micm::CpuSolverBuilder< - micm::RosenbrockSolverParameters, - micm::VectorMatrix, - micm::SparseMatrix>>(params); - - EXPECT_ANY_THROW(builder.SetSystem(the_system).SetReactions(reactions).SetNumberOfGridCells(1).Build();); -}*/ \ No newline at end of file +} \ No newline at end of file From 0904e97f003a6bb571a37775fdf5cb84eba4cf16 Mon Sep 17 00:00:00 2001 From: Montek Thind Date: Wed, 2 Oct 2024 10:32:16 -0700 Subject: [PATCH 06/31] fix cuda test --- test/unit/cuda/solver/test_cuda_rosenbrock.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/unit/cuda/solver/test_cuda_rosenbrock.cpp b/test/unit/cuda/solver/test_cuda_rosenbrock.cpp index 412af8745..2492aed2b 100644 --- a/test/unit/cuda/solver/test_cuda_rosenbrock.cpp +++ b/test/unit/cuda/solver/test_cuda_rosenbrock.cpp @@ -133,7 +133,7 @@ void testNormalizedErrorConst() y_new.CopyToDevice(); errors.CopyToDevice(); - double error = gpu_solver.solver_.NormalizedError(y_old, y_new, errors, state); + double error = gpu_solver.solver_.NormalizedError(y_old, y_new, errors); auto expected_error = 0.0; for (size_t i = 0; i < state.state_size_; ++i) @@ -192,7 +192,7 @@ void testNormalizedErrorDiff() y_new.CopyToDevice(); errors.CopyToDevice(); - double computed_error = gpu_solver.solver_.NormalizedError(y_old, y_new, errors, state); + double computed_error = gpu_solver.solver_.NormalizedError(y_old, y_new, errors); auto relative_error = std::abs(computed_error - expected_error) / std::max(std::abs(computed_error), std::abs(expected_error)); From f4cf73afc6bcc784141b590ccd81d11e3026ec7a Mon Sep 17 00:00:00 2001 From: Montek Thind Date: Wed, 2 Oct 2024 11:01:16 -0700 Subject: [PATCH 07/31] fix cuda tests --- include/micm/cuda/solver/cuda_rosenbrock.hpp | 2 +- test/integration/cuda/test_cuda_analytical_rosenbrock.cpp | 2 -- test/unit/cuda/solver/test_cuda_rosenbrock.cpp | 4 ++-- 3 files changed, 3 insertions(+), 5 deletions(-) diff --git a/include/micm/cuda/solver/cuda_rosenbrock.hpp b/include/micm/cuda/solver/cuda_rosenbrock.hpp index 2044f97c9..e74a8c0f1 100644 --- a/include/micm/cuda/solver/cuda_rosenbrock.hpp +++ b/include/micm/cuda/solver/cuda_rosenbrock.hpp @@ -116,7 +116,7 @@ namespace micm /// @return The scaled norm of the errors template double NormalizedError(const DenseMatrixPolicy& y_old, const DenseMatrixPolicy& y_new, const DenseMatrixPolicy& errors) - const requires(CudaMatrix&& VectorizableDense) + const requires(CudaMatrix&& VectorizableDense, auto& state) { return micm::cuda::NormalizedErrorDriver( y_old.AsDeviceParam(), y_new.AsDeviceParam(), errors.AsDeviceParam(), this->parameters_, this->devstruct_); diff --git a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp index d1551e98c..09c366e98 100644 --- a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp +++ b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp @@ -146,8 +146,6 @@ TEST(AnalyticalExamplesCudaRosenbrock, Robertson) { auto rosenbrock_solver = [](auto params) { - params.relative_tolerance_ = 1e-10; - params.absolute_tolerance_ = std::vector(3, params.relative_tolerance_ * 1e-2); return builderType1Cell(params); }; diff --git a/test/unit/cuda/solver/test_cuda_rosenbrock.cpp b/test/unit/cuda/solver/test_cuda_rosenbrock.cpp index 2492aed2b..412af8745 100644 --- a/test/unit/cuda/solver/test_cuda_rosenbrock.cpp +++ b/test/unit/cuda/solver/test_cuda_rosenbrock.cpp @@ -133,7 +133,7 @@ void testNormalizedErrorConst() y_new.CopyToDevice(); errors.CopyToDevice(); - double error = gpu_solver.solver_.NormalizedError(y_old, y_new, errors); + double error = gpu_solver.solver_.NormalizedError(y_old, y_new, errors, state); auto expected_error = 0.0; for (size_t i = 0; i < state.state_size_; ++i) @@ -192,7 +192,7 @@ void testNormalizedErrorDiff() y_new.CopyToDevice(); errors.CopyToDevice(); - double computed_error = gpu_solver.solver_.NormalizedError(y_old, y_new, errors); + double computed_error = gpu_solver.solver_.NormalizedError(y_old, y_new, errors, state); auto relative_error = std::abs(computed_error - expected_error) / std::max(std::abs(computed_error), std::abs(expected_error)); From d67add64fdbc9a48fc4feb85a2cec287b5ae07e1 Mon Sep 17 00:00:00 2001 From: Montek Thind Date: Wed, 2 Oct 2024 11:08:29 -0700 Subject: [PATCH 08/31] fix again --- include/micm/cuda/solver/cuda_rosenbrock.hpp | 4 ++-- test/integration/cuda/test_cuda_analytical_rosenbrock.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/include/micm/cuda/solver/cuda_rosenbrock.hpp b/include/micm/cuda/solver/cuda_rosenbrock.hpp index e74a8c0f1..cae46e76d 100644 --- a/include/micm/cuda/solver/cuda_rosenbrock.hpp +++ b/include/micm/cuda/solver/cuda_rosenbrock.hpp @@ -115,8 +115,8 @@ namespace micm /// @param errors The computed errors /// @return The scaled norm of the errors template - double NormalizedError(const DenseMatrixPolicy& y_old, const DenseMatrixPolicy& y_new, const DenseMatrixPolicy& errors) - const requires(CudaMatrix&& VectorizableDense, auto& state) + double NormalizedError(const DenseMatrixPolicy& y_old, const DenseMatrixPolicy& y_new, const DenseMatrixPolicy& errors, auto& state) + const requires(CudaMatrix&& VectorizableDense) { return micm::cuda::NormalizedErrorDriver( y_old.AsDeviceParam(), y_new.AsDeviceParam(), errors.AsDeviceParam(), this->parameters_, this->devstruct_); diff --git a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp index 09c366e98..69e93dab4 100644 --- a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp +++ b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp @@ -33,7 +33,7 @@ auto six_da_1_cell = builderType1Cell(micm::RosenbrockSolverParameters::SixStage auto copy_to_device = [](auto& state) -> void { state.SyncInputsToDevice(); }; auto copy_to_host = [](auto& state) -> void { state.SyncOutputsToHost(); }; - +/* TEST(AnalyticalExamplesCudaRosenbrock, Troe) { test_analytical_troe(two, 2e-3, copy_to_device, copy_to_host); @@ -241,4 +241,4 @@ TEST(AnalyticalExamplesCudaRosenbrock, HIRES) solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_hires(solver, 1e-3, copy_to_device, copy_to_host); -} \ No newline at end of file +}*/ \ No newline at end of file From ecb857a080920240e73192dec6cd266068fae5b7 Mon Sep 17 00:00:00 2001 From: Montek Thind Date: Wed, 2 Oct 2024 18:57:34 -0700 Subject: [PATCH 09/31] cuda test fix --- .../cuda/solver/test_cuda_solver_builder.cpp | 16 ---------------- test/unit/jit/solver/test_jit_solver_builder.cpp | 16 ---------------- 2 files changed, 32 deletions(-) diff --git a/test/unit/cuda/solver/test_cuda_solver_builder.cpp b/test/unit/cuda/solver/test_cuda_solver_builder.cpp index 9958e3b87..1b0ffcfec 100644 --- a/test/unit/cuda/solver/test_cuda_solver_builder.cpp +++ b/test/unit/cuda/solver/test_cuda_solver_builder.cpp @@ -47,20 +47,4 @@ TEST(SolverBuilder, CanBuildCudaSolvers) .SetReactions(reactions) .SetNumberOfGridCells(L) .Build(); -} - -TEST(SolverBuilder, MismatchedToleranceSizeIsCaught) -{ - auto params = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - // too many - params.absolute_tolerance_ = { 1e-6, 1e-6, 1e-6, 1e-6, 1e-6 }; - - constexpr std::size_t L = 4; - using cuda_builder = micm::CudaSolverBuilder; - - EXPECT_ANY_THROW(cuda_builder(params).SetSystem(the_system).SetReactions(reactions).SetNumberOfGridCells(L).Build();); - - // too few - params.absolute_tolerance_ = { 1e-6, 1e-6 }; - EXPECT_ANY_THROW(cuda_builder(params).SetSystem(the_system).SetReactions(reactions).SetNumberOfGridCells(L).Build();); } \ No newline at end of file diff --git a/test/unit/jit/solver/test_jit_solver_builder.cpp b/test/unit/jit/solver/test_jit_solver_builder.cpp index 87e0e0a02..7cf1f02c1 100644 --- a/test/unit/jit/solver/test_jit_solver_builder.cpp +++ b/test/unit/jit/solver/test_jit_solver_builder.cpp @@ -47,20 +47,4 @@ TEST(SolverBuilder, CanBuildJitRosenbrock) .SetReactions(reactions) .SetNumberOfGridCells(L) .Build(); -} - -TEST(SolverBuilder, MismatchedToleranceSizeIsCaught) -{ - auto params = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - // too many - params.absolute_tolerance_ = { 1e-6, 1e-6, 1e-6, 1e-6, 1e-6 }; - - constexpr std::size_t L = 4; - using jit_builder = micm::JitSolverBuilder; - - EXPECT_ANY_THROW(jit_builder(params).SetSystem(the_system).SetReactions(reactions).SetNumberOfGridCells(L).Build();); - - // too few - params.absolute_tolerance_ = { 1e-6, 1e-6 }; - EXPECT_ANY_THROW(jit_builder(params).SetSystem(the_system).SetReactions(reactions).SetNumberOfGridCells(L).Build();); } \ No newline at end of file From fefab051fb754b5bbda9769cbeeee0c953d2015e Mon Sep 17 00:00:00 2001 From: Montek Thind Date: Mon, 14 Oct 2024 08:59:42 -0700 Subject: [PATCH 10/31] removed params --- include/micm/solver/rosenbrock_solver_parameters.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/include/micm/solver/rosenbrock_solver_parameters.hpp b/include/micm/solver/rosenbrock_solver_parameters.hpp index 923e022b9..8ee8a5d18 100644 --- a/include/micm/solver/rosenbrock_solver_parameters.hpp +++ b/include/micm/solver/rosenbrock_solver_parameters.hpp @@ -59,8 +59,9 @@ namespace micm std::array alpha_{}; // Gamma_i = \sum_j gamma_{i,j} std::array gamma_{}; - std::vector absolute_tolerance_; - double relative_tolerance_{ 1e-6 }; + + //std::vector absolute_tolerance_; + //double relative_tolerance_{ 1e-6 }; // Print RosenbrockSolverParameters to console void Print() const; From c481a0bf887021d7d53f8f5671fa365f75c5183e Mon Sep 17 00:00:00 2001 From: Montek Thind Date: Mon, 14 Oct 2024 12:43:04 -0700 Subject: [PATCH 11/31] added tolerance back --- include/micm/solver/backward_euler_solver_parameters.hpp | 3 ++- include/micm/solver/rosenbrock_solver_parameters.hpp | 3 --- test/integration/cuda/test_cuda_analytical_rosenbrock.cpp | 4 ++-- test/regression/RosenbrockChapman/chapman_ode_solver.hpp | 1 - 4 files changed, 4 insertions(+), 7 deletions(-) diff --git a/include/micm/solver/backward_euler_solver_parameters.hpp b/include/micm/solver/backward_euler_solver_parameters.hpp index f9706e9b7..59b6ea03e 100644 --- a/include/micm/solver/backward_euler_solver_parameters.hpp +++ b/include/micm/solver/backward_euler_solver_parameters.hpp @@ -19,7 +19,8 @@ namespace micm template using SolverType = BackwardEuler; - //double relative_tolerance_{ 1.0e-8 }; + std::vector 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 diff --git a/include/micm/solver/rosenbrock_solver_parameters.hpp b/include/micm/solver/rosenbrock_solver_parameters.hpp index 8ee8a5d18..456daa21b 100644 --- a/include/micm/solver/rosenbrock_solver_parameters.hpp +++ b/include/micm/solver/rosenbrock_solver_parameters.hpp @@ -59,9 +59,6 @@ namespace micm std::array alpha_{}; // Gamma_i = \sum_j gamma_{i,j} std::array gamma_{}; - - //std::vector absolute_tolerance_; - //double relative_tolerance_{ 1e-6 }; // Print RosenbrockSolverParameters to console void Print() const; diff --git a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp index 69e93dab4..09c366e98 100644 --- a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp +++ b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp @@ -33,7 +33,7 @@ auto six_da_1_cell = builderType1Cell(micm::RosenbrockSolverParameters::SixStage auto copy_to_device = [](auto& state) -> void { state.SyncInputsToDevice(); }; auto copy_to_host = [](auto& state) -> void { state.SyncOutputsToHost(); }; -/* + TEST(AnalyticalExamplesCudaRosenbrock, Troe) { test_analytical_troe(two, 2e-3, copy_to_device, copy_to_host); @@ -241,4 +241,4 @@ TEST(AnalyticalExamplesCudaRosenbrock, HIRES) solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_hires(solver, 1e-3, copy_to_device, copy_to_host); -}*/ \ No newline at end of file +} \ No newline at end of file diff --git a/test/regression/RosenbrockChapman/chapman_ode_solver.hpp b/test/regression/RosenbrockChapman/chapman_ode_solver.hpp index 16d4c5d16..37f0155f6 100644 --- a/test/regression/RosenbrockChapman/chapman_ode_solver.hpp +++ b/test/regression/RosenbrockChapman/chapman_ode_solver.hpp @@ -789,7 +789,6 @@ namespace micm std::function)> is_successful = [](const std::vector& jacobian) { return true; }; std::vector ode_jacobian; - uint64_t n_consecutive = 0; double alpha = 1 / (H * gamma); // compute jacobian decomposition of alpha*I - dforce_dy From fdb1fb6da2bd10c1879efd3f1de095dafbd75410 Mon Sep 17 00:00:00 2001 From: Montek Thind Date: Wed, 16 Oct 2024 13:23:52 -0700 Subject: [PATCH 12/31] fix cuda tests --- include/micm/solver/rosenbrock_solver_parameters.hpp | 4 ++++ include/micm/solver/solver_builder.inl | 1 + test/integration/cuda/test_cuda_analytical_rosenbrock.cpp | 6 +++++- 3 files changed, 10 insertions(+), 1 deletion(-) diff --git a/include/micm/solver/rosenbrock_solver_parameters.hpp b/include/micm/solver/rosenbrock_solver_parameters.hpp index 456daa21b..e2ebfc675 100644 --- a/include/micm/solver/rosenbrock_solver_parameters.hpp +++ b/include/micm/solver/rosenbrock_solver_parameters.hpp @@ -60,6 +60,10 @@ namespace micm // Gamma_i = \sum_j gamma_{i,j} std::array gamma_{}; + //We have this to pass CUDA Tests. once we figure out a path forward we can get rid of these. + std::vector absolute_tolerance_; + double relative_tolerance_{ 1e-6 }; + // Print RosenbrockSolverParameters to console void Print() const; diff --git a/include/micm/solver/solver_builder.inl b/include/micm/solver/solver_builder.inl index 3bd429350..9623711c0 100644 --- a/include/micm/solver/solver_builder.inl +++ b/include/micm/solver/solver_builder.inl @@ -400,6 +400,7 @@ namespace micm this->SetAbsoluteTolerances(state_parameters.absolute_tolerance_, species_map); + options.absolute_tolerance_ = state_parameters.absolute_tolerance_; return Solver( SolverPolicy(options, std::move(linear_solver), std::move(rates), jacobian), diff --git a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp index 09c366e98..b9a9ed052 100644 --- a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp +++ b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp @@ -174,6 +174,9 @@ TEST(AnalyticalExamplesCudaRosenbrock, SurfaceRxn) test_analytical_surface_rxn(six_da_1_cell, 1e-7, copy_to_device, copy_to_host); } +//Commented out the following tests as it is failing to copy State to Device +//Will uncomment them once we have a plan how to deal with state (relative tolearnce and absolute tolerance) on GPU +/* TEST(AnalyticalExamplesCudaRosenbrock, E5) { auto rosenbrock_solver = [](auto params) @@ -241,4 +244,5 @@ TEST(AnalyticalExamplesCudaRosenbrock, HIRES) solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_hires(solver, 1e-3, copy_to_device, copy_to_host); -} \ No newline at end of file +} +*/ \ No newline at end of file From 9f4b7b39a56216ff1385c8b00cd014ee60a2463d Mon Sep 17 00:00:00 2001 From: Kyle Shores Date: Wed, 30 Oct 2024 08:57:59 -0600 Subject: [PATCH 13/31] trying to copy tolerances from state --- include/micm/cuda/solver/cuda_rosenbrock.cuh | 1 + include/micm/cuda/solver/cuda_rosenbrock.hpp | 4 +--- include/micm/cuda/solver/cuda_state.hpp | 1 + .../micm/solver/rosenbrock_solver_parameters.hpp | 2 -- include/micm/solver/solver_builder.inl | 1 - include/micm/solver/state.hpp | 2 +- include/micm/solver/state.inl | 5 ++++- include/micm/version.hpp | 2 ++ src/solver/rosenbrock.cu | 13 ++++++++----- 9 files changed, 18 insertions(+), 13 deletions(-) diff --git a/include/micm/cuda/solver/cuda_rosenbrock.cuh b/include/micm/cuda/solver/cuda_rosenbrock.cuh index 3a03fb8b9..8ab81c554 100644 --- a/include/micm/cuda/solver/cuda_rosenbrock.cuh +++ b/include/micm/cuda/solver/cuda_rosenbrock.cuh @@ -41,6 +41,7 @@ namespace micm const CudaMatrixParam& y_old_param, const CudaMatrixParam& y_new_param, const CudaMatrixParam& errors_param, + const CudaMatrixParam& absolute_tolerance_param, const RosenbrockSolverParameters& ros_param, CudaRosenbrockSolverParam devstruct); } // namespace cuda diff --git a/include/micm/cuda/solver/cuda_rosenbrock.hpp b/include/micm/cuda/solver/cuda_rosenbrock.hpp index cae46e76d..2c8317280 100644 --- a/include/micm/cuda/solver/cuda_rosenbrock.hpp +++ b/include/micm/cuda/solver/cuda_rosenbrock.hpp @@ -82,8 +82,6 @@ namespace micm hoststruct.errors_size_ = jacobian.GroupVectorSize() * this->parameters_.absolute_tolerance_.size(); 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); }; @@ -119,7 +117,7 @@ namespace micm const requires(CudaMatrix&& VectorizableDense) { return micm::cuda::NormalizedErrorDriver( - y_old.AsDeviceParam(), y_new.AsDeviceParam(), errors.AsDeviceParam(), this->parameters_, this->devstruct_); + y_old.AsDeviceParam(), y_new.AsDeviceParam(), errors.AsDeviceParam(), this->parameters_, state.absolute_tolerance_.AsDeviceParam(), this->devstruct_); } }; // end CudaRosenbrockSolver } // namespace micm \ No newline at end of file diff --git a/include/micm/cuda/solver/cuda_state.hpp b/include/micm/cuda/solver/cuda_state.hpp index 13ab3ac36..8fc0c1f21 100644 --- a/include/micm/cuda/solver/cuda_state.hpp +++ b/include/micm/cuda/solver/cuda_state.hpp @@ -29,6 +29,7 @@ namespace micm { this->variables_.CopyToDevice(); this->rate_constants_.CopyToDevice(); + this->absolute_tolerance_.CopyToDevice(); } /// @brief Copy output variables to the host diff --git a/include/micm/solver/rosenbrock_solver_parameters.hpp b/include/micm/solver/rosenbrock_solver_parameters.hpp index e2ebfc675..92b9e97d6 100644 --- a/include/micm/solver/rosenbrock_solver_parameters.hpp +++ b/include/micm/solver/rosenbrock_solver_parameters.hpp @@ -60,8 +60,6 @@ namespace micm // Gamma_i = \sum_j gamma_{i,j} std::array gamma_{}; - //We have this to pass CUDA Tests. once we figure out a path forward we can get rid of these. - std::vector absolute_tolerance_; double relative_tolerance_{ 1e-6 }; // Print RosenbrockSolverParameters to console diff --git a/include/micm/solver/solver_builder.inl b/include/micm/solver/solver_builder.inl index 9623711c0..3bd429350 100644 --- a/include/micm/solver/solver_builder.inl +++ b/include/micm/solver/solver_builder.inl @@ -400,7 +400,6 @@ namespace micm this->SetAbsoluteTolerances(state_parameters.absolute_tolerance_, species_map); - options.absolute_tolerance_ = state_parameters.absolute_tolerance_; return Solver( SolverPolicy(options, std::move(linear_solver), std::move(rates), jacobian), diff --git a/include/micm/solver/state.hpp b/include/micm/solver/state.hpp index 91cfae3fd..5bbf06aa3 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_; - std::vector absolute_tolerance_; + DenseMatrixPolicy 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 37f3962b5..2576c6c1f 100644 --- a/include/micm/solver/state.inl +++ b/include/micm/solver/state.inl @@ -82,8 +82,11 @@ namespace micm state_size_(parameters.variable_names_.size()), number_of_grid_cells_(parameters.number_of_grid_cells_), relative_tolerance_(1e-06), - absolute_tolerance_(parameters.absolute_tolerance_) + absolute_tolerance_() { + absolute_tolerance_ = DenseMatrixPolicy(1, parameters.variable_names_.size()); + absolute_tolerance_.AsVector() = parameters.absolute_tolerance_; + std::size_t index = 0; for (auto& name : variable_names_) variable_map_[name] = index++; diff --git a/include/micm/version.hpp b/include/micm/version.hpp index ee99b8093..ab33a2d99 100644 --- a/include/micm/version.hpp +++ b/include/micm/version.hpp @@ -1,3 +1,5 @@ +// Copyright (C) 2023-2024 National Center for Atmospheric Research +// SPDX-License-Identifier: Apache-2.0 // clang-format off #pragma once diff --git a/src/solver/rosenbrock.cu b/src/solver/rosenbrock.cu index 6a6ab6961..e7053da1a 100644 --- a/src/solver/rosenbrock.cu +++ b/src/solver/rosenbrock.cu @@ -137,6 +137,7 @@ namespace micm __global__ void NormalizedErrorKernel( const CudaMatrixParam y_old_param, const CudaMatrixParam y_new_param, + const CudaMatrixParam absolute_tolerance_param, const RosenbrockSolverParameters ros_param, CudaRosenbrockSolverParam devstruct, const size_t n, @@ -245,6 +246,7 @@ namespace micm __global__ void ScaledErrorKernel( const CudaMatrixParam y_old_param, const CudaMatrixParam y_new_param, + const CudaMatrixParam absolute_tolerance_param, const RosenbrockSolverParameters ros_param, CudaRosenbrockSolverParam devstruct) { @@ -253,7 +255,7 @@ namespace micm const double* const d_y_old = y_old_param.d_data_; const double* const d_y_new = y_new_param.d_data_; double* const d_errors = devstruct.errors_input_; - const double* const atol = devstruct.absolute_tolerance_; + const double* const atol = absolute_tolerance_param.d_data_; const double rtol = ros_param.relative_tolerance_; const size_t num_elements = devstruct.errors_size_; const size_t number_of_grid_cells = y_old_param.number_of_grid_cells_; @@ -288,6 +290,7 @@ namespace micm const CudaMatrixParam& y_old_param, const CudaMatrixParam& y_new_param, const CudaMatrixParam& errors_param, + const CudaMatrixParam& absolute_tolerance_param, const RosenbrockSolverParameters& ros_param, CudaRosenbrockSolverParam devstruct) { @@ -318,7 +321,7 @@ namespace micm BLOCK_SIZE, 0, micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)>>>( - y_old_param, y_new_param, ros_param, devstruct); + y_old_param, y_new_param, absolute_tolerance_param, ros_param, devstruct); // call cublas function to perform the norm: // https://docs.nvidia.com/cuda/cublas/index.html?highlight=dnrm2#cublas-t-nrm2 CHECK_CUBLAS_ERROR( @@ -341,7 +344,7 @@ namespace micm BLOCK_SIZE, BLOCK_SIZE * sizeof(double), micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)>>>( - y_old_param, y_new_param, ros_param, devstruct, number_of_elements, is_first_call); + y_old_param, y_new_param, absolute_tolerance_param, ros_param, devstruct, number_of_elements, is_first_call); is_first_call = false; while (number_of_blocks > 1) { @@ -355,7 +358,7 @@ namespace micm BLOCK_SIZE, BLOCK_SIZE * sizeof(double), micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)>>>( - y_old_param, y_new_param, ros_param, devstruct, number_of_blocks, is_first_call); + y_old_param, y_new_param, absolute_tolerance_param, ros_param, devstruct, number_of_blocks, is_first_call); break; } NormalizedErrorKernel<<< @@ -363,7 +366,7 @@ namespace micm BLOCK_SIZE, BLOCK_SIZE * sizeof(double), micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)>>>( - y_old_param, y_new_param, ros_param, devstruct, number_of_blocks, is_first_call); + y_old_param, y_new_param, absolute_tolerance_param, ros_param, devstruct, number_of_blocks, is_first_call); number_of_blocks = new_number_of_blocks; } From 4c84691b99d43187d2605c3ec43ebbb9061a19a6 Mon Sep 17 00:00:00 2001 From: Kyle Shores Date: Wed, 6 Nov 2024 11:44:56 -0700 Subject: [PATCH 14/31] updating usage of absolute tolerance --- include/micm/cuda/solver/cuda_rosenbrock.hpp | 15 +++++++++++---- include/micm/solver/rosenbrock.hpp | 9 ++++++--- include/micm/solver/rosenbrock.inl | 10 +++++----- include/micm/solver/solver_builder.inl | 2 +- test/unit/cuda/solver/test_cuda_rosenbrock.cpp | 4 ++-- test/unit/solver/test_rosenbrock.cpp | 8 ++++---- 6 files changed, 29 insertions(+), 19 deletions(-) diff --git a/include/micm/cuda/solver/cuda_rosenbrock.hpp b/include/micm/cuda/solver/cuda_rosenbrock.hpp index 2c8317280..084d8c205 100644 --- a/include/micm/cuda/solver/cuda_rosenbrock.hpp +++ b/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>( 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&& VectorizableDense) { 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 \ No newline at end of file diff --git a/include/micm/solver/rosenbrock.hpp b/include/micm/solver/rosenbrock.hpp index 3939315ad..da359706a 100644 --- a/include/micm/solver/rosenbrock.hpp +++ b/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>( parameters, std::move(linear_solver), std::move(rates), - jacobian) + jacobian, + number_of_species) { } RosenbrockSolver(const RosenbrockSolver&) = delete; diff --git a/include/micm/solver/rosenbrock.inl b/include/micm/solver/rosenbrock.inl index 2e6484b03..939c9521c 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_.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; } diff --git a/include/micm/solver/solver_builder.inl b/include/micm/solver/solver_builder.inl index 3bd429350..4f5bc76d4 100644 --- a/include/micm/solver/solver_builder.inl +++ b/include/micm/solver/solver_builder.inl @@ -402,7 +402,7 @@ namespace micm this->SetAbsoluteTolerances(state_parameters.absolute_tolerance_, species_map); return Solver( - 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_); diff --git a/test/unit/cuda/solver/test_cuda_rosenbrock.cpp b/test/unit/cuda/solver/test_cuda_rosenbrock.cpp index 412af8745..b3bcaa7bb 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_; + std::vector atol = state.absolute_tolerance_.AsVector(); double rtol = state.relative_tolerance_; auto y_old = micm::CudaDenseMatrix(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 atol = state.absolute_tolerance_; + std::vector atol = state.absolute_tolerance_.AsVector(); 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); diff --git a/test/unit/solver/test_rosenbrock.cpp b/test/unit/solver/test_rosenbrock.cpp index 01aa8a510..ee760b54c 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_; + std::vector 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); } } From 1b67f7982539eae96113e451c6be3c8f0ff56d16 Mon Sep 17 00:00:00 2001 From: Kyle Shores Date: Thu, 7 Nov 2024 12:38:25 -0700 Subject: [PATCH 15/31] removed compile errors --- include/micm/solver/backward_euler.hpp | 15 +++++----- include/micm/solver/backward_euler.inl | 28 +++++++++---------- .../backward_euler_solver_parameters.hpp | 2 -- test/integration/analytical_policy.hpp | 14 +++++----- test/integration/test_chapman_integration.cpp | 17 +++++------ test/unit/solver/test_backward_euler.cpp | 25 ++++++++--------- 6 files changed, 50 insertions(+), 51 deletions(-) diff --git a/include/micm/solver/backward_euler.hpp b/include/micm/solver/backward_euler.hpp index a0bfdf861..c3a6c65fa 100644 --- a/include/micm/solver/backward_euler.hpp +++ b/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 static bool IsConverged( - const BackwardEulerSolverParameters& parameters, - const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& state, auto& stateParams) requires(!VectorizableDense); + const BackwardEulerSolverParameters& parameters, + const DenseMatrixPolicy& residual, + const DenseMatrixPolicy& Yn1, std::vector& absolute_tolerance, double relative_tolerance) requires(!VectorizableDense); template static bool IsConverged( - const BackwardEulerSolverParameters& parameters, - const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& state, auto& stateParams) requires(VectorizableDense); + const BackwardEulerSolverParameters& parameters, + const DenseMatrixPolicy& residual, + const DenseMatrixPolicy& Yn1, std::vector& absolute_tolerance, double relative_tolerance) requires(VectorizableDense); }; } // namespace micm diff --git a/include/micm/solver/backward_euler.inl b/include/micm/solver/backward_euler.inl index 5bb48ea95..427ceea12 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); + converged = IsConverged(parameters_, forcing, Yn1, state.absolute_tolerance_.AsVector(), state.relative_tolerance_); } while (!converged && iterations < max_iter); if (!converged) @@ -180,23 +180,23 @@ namespace micm inline bool BackwardEuler::IsConverged( const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& state, auto& stateParams) requires(!VectorizableDense) + const DenseMatrixPolicy& Yn1, std::vector& absolute_tolerance, double relative_tolerance) requires(!VectorizableDense) { double small = parameters.small_; - double rel_tol = stateParams.relative_tolerance_; - auto& abs_tol = stateParams.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,13 +206,13 @@ namespace micm inline bool BackwardEuler::IsConverged( const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& state, auto& stateParams) requires(VectorizableDense) + const DenseMatrixPolicy& Yn1, std::vector& absolute_tolerance, double relative_tolerance) requires(VectorizableDense) { double small = parameters.small_; - double rel_tol = stateParams.relative_tolerance_; - auto& abs_tol = stateParams.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(); @@ -221,11 +221,11 @@ namespace micm 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 @@ -238,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; } diff --git a/include/micm/solver/backward_euler_solver_parameters.hpp b/include/micm/solver/backward_euler_solver_parameters.hpp index 59b6ea03e..3226c2198 100644 --- a/include/micm/solver/backward_euler_solver_parameters.hpp +++ b/include/micm/solver/backward_euler_solver_parameters.hpp @@ -19,8 +19,6 @@ namespace micm template using SolverType = BackwardEuler; - std::vector 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 diff --git a/test/integration/analytical_policy.hpp b/test/integration/analytical_policy.hpp index 13800603f..556545f78 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_ = std::vector(3, state.relative_tolerance_ * 1e-2); + state.absolute_tolerance_.AsVector() = 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-6; - state.absolute_tolerance_ = std::vector(5, state.relative_tolerance_ * 1e-2); + state.absolute_tolerance_.AsVector() = std::vector(5, state.relative_tolerance_ * 1e-2); 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_ = std::vector(8, state.relative_tolerance_ * 1e-2); + state.absolute_tolerance_.AsVector() = 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_ = std::vector(6, 1e-17); - state.absolute_tolerance_[0] = 1e-7; - state.absolute_tolerance_[4] = 1e-7; - state.absolute_tolerance_[5] = 1e-7; + 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.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 8ebdc001f..cbd876f41 100644 --- a/test/integration/test_chapman_integration.cpp +++ b/test/integration/test_chapman_integration.cpp @@ -43,17 +43,18 @@ 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(); for (size_t n_grid_cell = 0; n_grid_cell < state.number_of_grid_cells_; ++n_grid_cell) { - EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["Ar"]], 1.0e-12); - EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["CO2"]], 1.0e-12); - EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["H2O"]], 1.0e-12); - EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["N2"]], 1.0e-12); - EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["O1D"]], 1.0e-12); - EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["O"]], 1.0e-12); - EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["O2"]], 1.0e-12); - EXPECT_EQ(state.absolute_tolerance_[state.variable_map_["O3"]], 1.0e-12); + EXPECT_EQ(abs_tol[state.variable_map_["Ar"]], 1.0e-12); + EXPECT_EQ(abs_tol[state.variable_map_["CO2"]], 1.0e-12); + EXPECT_EQ(abs_tol[state.variable_map_["H2O"]], 1.0e-12); + EXPECT_EQ(abs_tol[state.variable_map_["N2"]], 1.0e-12); + EXPECT_EQ(abs_tol[state.variable_map_["O1D"]], 1.0e-12); + EXPECT_EQ(abs_tol[state.variable_map_["O"]], 1.0e-12); + EXPECT_EQ(abs_tol[state.variable_map_["O2"]], 1.0e-12); + EXPECT_EQ(abs_tol[state.variable_map_["O3"]], 1.0e-12); } // User gives an input of concentrations diff --git a/test/unit/solver/test_backward_euler.cpp b/test/unit/solver/test_backward_euler.cpp index b48141585..eb976b350 100644 --- a/test/unit/solver/test_backward_euler.cpp +++ b/test/unit/solver/test_backward_euler.cpp @@ -68,26 +68,25 @@ void CheckIsConverged() micm::BackwardEulerSolverParameters parameters; DenseMatrixPolicy residual{ 4, 3, 0.0 }; - DenseMatrixPolicy state{ 4, 3, 0.0 }; - micm::StateParameters state_params; + DenseMatrixPolicy Yn1{ 4, 3, 0.0 }; parameters.small_ = 1e-6; - state_params.relative_tolerance_ = 1e-3; - state_params.absolute_tolerance_ = { 1e-6, 1e-6, 1e-6 }; + double relative_tolerance = 1e-3; + std::vector absolute_tolerance = { 1e-6, 1e-6, 1e-6 }; - ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state, state_params)); + ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance)); residual[0][1] = 1e-5; - ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, state, state_params)); + ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance)); parameters.small_ = 1e-4; - ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state, state_params)); + ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance)); residual[3][2] = 1e-3; - ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, state, state_params)); - state[3][2] = 10.0; - ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state, state_params)); + ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance)); + Yn1[3][2] = 10.0; + ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance)); residual[3][2] = 1e-1; - ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, state, state_params)); - state_params.absolute_tolerance_[2] = 1.0; - ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state, state_params)); + ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance)); + absolute_tolerance[2] = 1.0; + ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance)); } TEST(BackwardEuler, IsConverged) From 3e0d09fef6a26e44e57cc9d6ce8da334978b1551 Mon Sep 17 00:00:00 2001 From: Kyle Shores Date: Thu, 7 Nov 2024 14:14:31 -0700 Subject: [PATCH 16/31] moving relative tolerance --- include/micm/cuda/solver/cuda_rosenbrock.cuh | 2 +- include/micm/cuda/solver/cuda_rosenbrock.hpp | 3 +- include/micm/cuda/solver/cuda_state.hpp | 16 ++++++++ include/micm/cuda/util/cuda_param.hpp | 2 - .../solver/rosenbrock_solver_parameters.hpp | 2 - include/micm/solver/state.hpp | 2 + src/solver/rosenbrock.cu | 41 +++++-------------- 7 files changed, 30 insertions(+), 38 deletions(-) diff --git a/include/micm/cuda/solver/cuda_rosenbrock.cuh b/include/micm/cuda/solver/cuda_rosenbrock.cuh index 8ab81c554..95cfd1e9a 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 RosenbrockSolverParameters& ros_param, + 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 084d8c205..15cbfea6d 100644 --- 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; }; @@ -123,7 +122,7 @@ namespace micm y_new.AsDeviceParam(), errors.AsDeviceParam(), state.absolute_tolerance_.AsDeviceParam(), - this->parameters_, + state.cuda_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 8fc0c1f21..38428a318 100644 --- a/include/micm/cuda/solver/cuda_state.hpp +++ b/include/micm/cuda/solver/cuda_state.hpp @@ -7,6 +7,8 @@ #include #include +#include + namespace micm { /// @brief Construct a state variable for CUDA tests @@ -19,6 +21,13 @@ namespace micm CudaState(CudaState&&) = default; CudaState& operator=(CudaState&&) = default; + double* cuda_relative_tolerance_{ nullptr }; + + ~CudaState() + { + CHECK_CUDA_ERROR(cudaFree(cuda_relative_tolerance_), "cudaFree"); + } + /// @brief Constructor which takes the state dimension information as input /// @param parameters State dimension information CudaState(const StateParameters& parameters) @@ -30,6 +39,13 @@ namespace micm this->variables_.CopyToDevice(); this->rate_constants_.CopyToDevice(); this->absolute_tolerance_.CopyToDevice(); + + size_t relative_tolerance_bytes = sizeof(double); + + CHECK_CUDA_ERROR( + cudaMallocAsync( + &(cuda_relative_tolerance_), relative_tolerance_bytes, micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)), + "cudaMalloc"); } /// @brief Copy output variables to the host diff --git a/include/micm/cuda/util/cuda_param.hpp b/include/micm/cuda/util/cuda_param.hpp index 766587c14..b16548b33 100644 --- a/include/micm/cuda/util/cuda_param.hpp +++ b/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; diff --git a/include/micm/solver/rosenbrock_solver_parameters.hpp b/include/micm/solver/rosenbrock_solver_parameters.hpp index 92b9e97d6..456daa21b 100644 --- a/include/micm/solver/rosenbrock_solver_parameters.hpp +++ b/include/micm/solver/rosenbrock_solver_parameters.hpp @@ -60,8 +60,6 @@ namespace micm // Gamma_i = \sum_j gamma_{i,j} std::array gamma_{}; - double relative_tolerance_{ 1e-6 }; - // Print RosenbrockSolverParameters to console void Print() const; diff --git a/include/micm/solver/state.hpp b/include/micm/solver/state.hpp index 5bbf06aa3..14dbd6406 100644 --- a/include/micm/solver/state.hpp +++ b/include/micm/solver/state.hpp @@ -164,6 +164,8 @@ namespace micm return *this; } + virtual ~State() = default; + /// @brief Set species' concentrations /// @param species_to_concentration void SetConcentrations(const std::unordered_map>& species_to_concentration); diff --git a/src/solver/rosenbrock.cu b/src/solver/rosenbrock.cu index e7053da1a..dc725b0a3 100644 --- a/src/solver/rosenbrock.cu +++ b/src/solver/rosenbrock.cu @@ -41,7 +41,6 @@ namespace micm /// Calculate the memory space of each temporary variable size_t errors_bytes = sizeof(double) * hoststruct.errors_size_; - size_t tolerance_bytes = sizeof(double) * hoststruct.absolute_tolerance_size_; /// Create a struct whose members contain the addresses in the device memory. CudaRosenbrockSolverParam devstruct; @@ -59,12 +58,6 @@ namespace micm jacobian_diagonal_elements_bytes, micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)), "cudaMalloc"); - CHECK_CUDA_ERROR( - cudaMallocAsync( - &(devstruct.absolute_tolerance_), - tolerance_bytes, - micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)), - "cudaMalloc"); /// Copy the data from host to device CHECK_CUDA_ERROR( @@ -76,18 +69,8 @@ namespace micm micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)), "cudaMemcpy"); - CHECK_CUDA_ERROR( - cudaMemcpyAsync( - devstruct.absolute_tolerance_, - hoststruct.absolute_tolerance_, - tolerance_bytes, - cudaMemcpyHostToDevice, - micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)), - "cudaMemcpy"); - devstruct.errors_size_ = hoststruct.errors_size_; devstruct.jacobian_diagonal_elements_size_ = hoststruct.jacobian_diagonal_elements_size_; - devstruct.absolute_tolerance_size_ = hoststruct.absolute_tolerance_size_; return devstruct; } @@ -109,10 +92,6 @@ namespace micm cudaFreeAsync( devstruct.jacobian_diagonal_elements_, micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)), "cudaFree"); - if (devstruct.absolute_tolerance_ != nullptr) - CHECK_CUDA_ERROR( - cudaFreeAsync(devstruct.absolute_tolerance_, micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)), - "cudaFree"); } // Specific CUDA device function to do reduction within a warp @@ -138,7 +117,7 @@ namespace micm const CudaMatrixParam y_old_param, const CudaMatrixParam y_new_param, const CudaMatrixParam absolute_tolerance_param, - const RosenbrockSolverParameters ros_param, + const double* relative_tolerance, CudaRosenbrockSolverParam devstruct, const size_t n, bool is_first_call) @@ -147,8 +126,8 @@ namespace micm const double* const d_y_new = y_new_param.d_data_; double* const d_errors_input = devstruct.errors_input_; double* const d_errors_output = devstruct.errors_output_; - const double* const atol = devstruct.absolute_tolerance_; - const double rtol = ros_param.relative_tolerance_; + const double* const atol = absolute_tolerance_param.d_data_; + 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. @@ -247,7 +226,7 @@ namespace micm const CudaMatrixParam y_old_param, const CudaMatrixParam y_new_param, const CudaMatrixParam absolute_tolerance_param, - const RosenbrockSolverParameters ros_param, + const double* relative_tolerance, CudaRosenbrockSolverParam devstruct) { // Local device variables @@ -256,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 = ros_param.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_; @@ -291,7 +270,7 @@ namespace micm const CudaMatrixParam& y_new_param, const CudaMatrixParam& errors_param, const CudaMatrixParam& absolute_tolerance_param, - const RosenbrockSolverParameters& ros_param, + const double* relative_tolerance, CudaRosenbrockSolverParam devstruct) { double normalized_error; @@ -321,7 +300,7 @@ namespace micm BLOCK_SIZE, 0, micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)>>>( - y_old_param, y_new_param, absolute_tolerance_param, ros_param, devstruct); + y_old_param, y_new_param, absolute_tolerance_param, relative_tolerance, devstruct); // call cublas function to perform the norm: // https://docs.nvidia.com/cuda/cublas/index.html?highlight=dnrm2#cublas-t-nrm2 CHECK_CUBLAS_ERROR( @@ -344,7 +323,7 @@ namespace micm BLOCK_SIZE, BLOCK_SIZE * sizeof(double), micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)>>>( - y_old_param, y_new_param, absolute_tolerance_param, ros_param, devstruct, number_of_elements, is_first_call); + y_old_param, y_new_param, absolute_tolerance_param, relative_tolerance, devstruct, number_of_elements, is_first_call); is_first_call = false; while (number_of_blocks > 1) { @@ -358,7 +337,7 @@ namespace micm BLOCK_SIZE, BLOCK_SIZE * sizeof(double), micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)>>>( - y_old_param, y_new_param, absolute_tolerance_param, ros_param, devstruct, number_of_blocks, is_first_call); + y_old_param, y_new_param, absolute_tolerance_param, relative_tolerance, devstruct, number_of_blocks, is_first_call); break; } NormalizedErrorKernel<<< @@ -366,7 +345,7 @@ namespace micm BLOCK_SIZE, BLOCK_SIZE * sizeof(double), micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)>>>( - y_old_param, y_new_param, absolute_tolerance_param, ros_param, devstruct, number_of_blocks, is_first_call); + y_old_param, y_new_param, absolute_tolerance_param, relative_tolerance, devstruct, number_of_blocks, is_first_call); number_of_blocks = new_number_of_blocks; } From a5796b125d9f2208a1e0129b0d56d9b4b3d702c4 Mon Sep 17 00:00:00 2001 From: Kyle Shores Date: Thu, 7 Nov 2024 15:55:13 -0700 Subject: [PATCH 17/31] fixed some backward euler stuff --- include/micm/solver/backward_euler.inl | 5 +++++ include/micm/solver/state.hpp | 4 +++- include/micm/solver/state.inl | 5 +++-- test/unit/solver/test_backward_euler.cpp | 2 +- 4 files changed, 12 insertions(+), 4 deletions(-) diff --git a/include/micm/solver/backward_euler.inl b/include/micm/solver/backward_euler.inl index 427ceea12..c11f7f759 100644 --- a/include/micm/solver/backward_euler.inl +++ b/include/micm/solver/backward_euler.inl @@ -182,6 +182,11 @@ namespace micm const DenseMatrixPolicy& residual, const DenseMatrixPolicy& Yn1, std::vector& absolute_tolerance, double relative_tolerance) requires(!VectorizableDense) { + std::cout << "rel_tol: " << relative_tolerance << std::endl; + for(auto& val : absolute_tolerance) + { + std::cout << "abs_tol: " << val << std::endl; + } double small = parameters.small_; double rel_tol = relative_tolerance; auto& abs_tol = absolute_tolerance; diff --git a/include/micm/solver/state.hpp b/include/micm/solver/state.hpp index 14dbd6406..6fe8945c6 100644 --- a/include/micm/solver/state.hpp +++ b/include/micm/solver/state.hpp @@ -133,7 +133,9 @@ namespace micm upper_matrix_(std::move(other.upper_matrix_)), state_size_(other.state_size_), number_of_grid_cells_(other.number_of_grid_cells_), - temporary_variables_(std::move(other.temporary_variables_)) + temporary_variables_(std::move(other.temporary_variables_)), + relative_tolerance_(other.relative_tolerance_), + absolute_tolerance_(std::move(other.absolute_tolerance_)) { } diff --git a/include/micm/solver/state.inl b/include/micm/solver/state.inl index 2576c6c1f..b9aa211c4 100644 --- a/include/micm/solver/state.inl +++ b/include/micm/solver/state.inl @@ -81,11 +81,12 @@ namespace micm upper_matrix_(), state_size_(parameters.variable_names_.size()), number_of_grid_cells_(parameters.number_of_grid_cells_), - relative_tolerance_(1e-06), + relative_tolerance_(parameters.relative_tolerance_), absolute_tolerance_() { - absolute_tolerance_ = DenseMatrixPolicy(1, parameters.variable_names_.size()); + std::cout << "absolute tolerance size" << absolute_tolerance_.AsVector().size() << std::endl; absolute_tolerance_.AsVector() = parameters.absolute_tolerance_; + std::cout << "absolute tolerance size" << absolute_tolerance_.AsVector().size() << std::endl; std::size_t index = 0; for (auto& name : variable_names_) diff --git a/test/unit/solver/test_backward_euler.cpp b/test/unit/solver/test_backward_euler.cpp index eb976b350..dd6d4f1b6 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_ = { 1e-6, 1e-6, 1e-6 }; + state.absolute_tolerance_.AsVector() = { 1e-6, 1e-6, 1e-6 }; state.variables_[0] = { 1.0, 0.0, 0.0 }; state.conditions_[0].temperature_ = 272.5; From 752519e6c92f2eecaf9dc77a6c3c2295251fba42 Mon Sep 17 00:00:00 2001 From: Kyle Shores Date: Fri, 8 Nov 2024 12:35:12 -0700 Subject: [PATCH 18/31] removing debug print statements --- include/micm/solver/state.inl | 2 -- 1 file changed, 2 deletions(-) diff --git a/include/micm/solver/state.inl b/include/micm/solver/state.inl index b9aa211c4..cd6a7e59b 100644 --- a/include/micm/solver/state.inl +++ b/include/micm/solver/state.inl @@ -84,9 +84,7 @@ namespace micm relative_tolerance_(parameters.relative_tolerance_), absolute_tolerance_() { - std::cout << "absolute tolerance size" << absolute_tolerance_.AsVector().size() << std::endl; absolute_tolerance_.AsVector() = parameters.absolute_tolerance_; - std::cout << "absolute tolerance size" << absolute_tolerance_.AsVector().size() << std::endl; std::size_t index = 0; for (auto& name : variable_names_) From 601b79e6e8464a1e375f0d4b7774fc02a5f2d1c5 Mon Sep 17 00:00:00 2001 From: Kyle Shores Date: Fri, 8 Nov 2024 15:00:07 -0600 Subject: [PATCH 19/31] correctly setting oregonator tolerances --- test/integration/analytical_policy.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/integration/analytical_policy.hpp b/test/integration/analytical_policy.hpp index 556545f78..8554db6cf 100644 --- a/test/integration/analytical_policy.hpp +++ b/test/integration/analytical_policy.hpp @@ -1547,8 +1547,8 @@ void test_analytical_oregonator( auto state = solver.GetState(); - state.relative_tolerance_ = 1e-6; - state.absolute_tolerance_.AsVector() = std::vector(5, state.relative_tolerance_ * 1e-2); + state.relative_tolerance_ = 1e-8; + state.absolute_tolerance_.AsVector() = std::vector(5, state.relative_tolerance_ * 1e-6); state.SetCustomRateParameter("r1", 1.34 * 0.06); state.SetCustomRateParameter("r2", 1.6e9); From aee4f4fbd4d1df8d52a69321a63e5f89b3a85dc1 Mon Sep 17 00:00:00 2001 From: Kyle Shores Date: Fri, 8 Nov 2024 15:43:09 -0700 Subject: [PATCH 20/31] correctly using tolerances --- include/micm/cuda/solver/cuda_rosenbrock.cuh | 2 +- include/micm/cuda/solver/cuda_rosenbrock.hpp | 4 ++-- include/micm/cuda/solver/cuda_state.hpp | 15 +++++++-------- include/micm/solver/backward_euler.inl | 2 +- include/micm/solver/rosenbrock.inl | 10 +++++----- include/micm/solver/state.hpp | 2 +- include/micm/solver/state.inl | 4 +--- src/solver/rosenbrock.cu | 10 +++++----- test/integration/analytical_policy.hpp | 14 +++++++------- test/integration/test_chapman_integration.cpp | 2 +- test/unit/cuda/solver/test_cuda_rosenbrock.cpp | 8 ++++++-- test/unit/solver/test_backward_euler.cpp | 2 +- test/unit/solver/test_rosenbrock.cpp | 8 ++++---- 13 files changed, 42 insertions(+), 41 deletions(-) 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); } } From f42bfca4c7d66adc07d8231928b3d287965046e0 Mon Sep 17 00:00:00 2001 From: Kyle Shores Date: Fri, 8 Nov 2024 16:48:23 -0600 Subject: [PATCH 21/31] updating jit constructor --- include/micm/jit/solver/jit_rosenbrock.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/micm/jit/solver/jit_rosenbrock.hpp b/include/micm/jit/solver/jit_rosenbrock.hpp index 36d7e3ad0..8bef434e3 100644 --- a/include/micm/jit/solver/jit_rosenbrock.hpp +++ b/include/micm/jit/solver/jit_rosenbrock.hpp @@ -73,7 +73,8 @@ namespace micm RosenbrockSolverParameters parameters, LinearSolverPolicy linear_solver, RatesPolicy rates, - auto& jacobian) + auto& jacobian, + const size_t number_of_species) : AbstractRosenbrockSolver>( parameters, std::move(linear_solver), From 96ddc42647c7ebc292864bc7cf9e0add9cd4369b Mon Sep 17 00:00:00 2001 From: Kyle Shores Date: Fri, 8 Nov 2024 16:53:03 -0600 Subject: [PATCH 22/31] passing the number of species to the constructor --- include/micm/jit/solver/jit_rosenbrock.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/micm/jit/solver/jit_rosenbrock.hpp b/include/micm/jit/solver/jit_rosenbrock.hpp index 8bef434e3..cc6ba45ee 100644 --- a/include/micm/jit/solver/jit_rosenbrock.hpp +++ b/include/micm/jit/solver/jit_rosenbrock.hpp @@ -79,7 +79,7 @@ namespace micm parameters, std::move(linear_solver), std::move(rates), - jacobian) + jacobian, number_of_species) { this->GenerateAlphaMinusJacobian(jacobian); } From 96fbfd013bb3720198f9a3a3fd4ec46b7bd56d18 Mon Sep 17 00:00:00 2001 From: Kyle Shores Date: Mon, 11 Nov 2024 08:33:41 -0600 Subject: [PATCH 23/31] correcting jit tests --- test/integration/jit/test_jit_analytical_rosenbrock.cpp | 2 -- test/integration/jit/test_jit_terminator.cpp | 1 - 2 files changed, 3 deletions(-) diff --git a/test/integration/jit/test_jit_analytical_rosenbrock.cpp b/test/integration/jit/test_jit_analytical_rosenbrock.cpp index ac33be82f..3792563da 100644 --- a/test/integration/jit/test_jit_analytical_rosenbrock.cpp +++ b/test/integration/jit/test_jit_analytical_rosenbrock.cpp @@ -148,8 +148,6 @@ TEST(AnalyticalExamplesJitRosenbrock, Robertson) { auto rosenbrock_solver = [](auto params) { - params.relative_tolerance_ = 1e-10; - params.absolute_tolerance_ = std::vector(3, params.relative_tolerance_ * 1e-2); return BuilderType<1>(params); }; diff --git a/test/integration/jit/test_jit_terminator.cpp b/test/integration/jit/test_jit_terminator.cpp index 0e633e265..2e80f1f00 100644 --- a/test/integration/jit/test_jit_terminator.cpp +++ b/test/integration/jit/test_jit_terminator.cpp @@ -27,7 +27,6 @@ using Group4SparseVectorMatrix = micm::SparseMatrix(parameters).SetIgnoreUnusedSpecies(true); From 09c4df23c4042202bb6d87634e4cb5b7720644a4 Mon Sep 17 00:00:00 2001 From: Kyle Shores Date: Mon, 11 Nov 2024 08:40:44 -0600 Subject: [PATCH 24/31] uncommenting tests --- .../cuda/test_cuda_analytical_rosenbrock.cpp | 64 ++++++------------- .../jit/test_jit_analytical_rosenbrock.cpp | 60 ++++++----------- .../test_analytical_rosenbrock.cpp | 60 ++++++----------- 3 files changed, 60 insertions(+), 124 deletions(-) diff --git a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp index c64eaf443..28fbdc1e0 100644 --- a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp +++ b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp @@ -144,24 +144,19 @@ TEST(AnalyticalExamplesCudaRosenbrock, BranchedSuperStiffButAnalytical) TEST(AnalyticalExamplesCudaRosenbrock, Robertson) { - auto rosenbrock_solver = [](auto params) - { - return builderType1Cell(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); + auto solver = builderType1Cell(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); test_analytical_robertson(solver, 2e-1, copy_to_device, copy_to_host); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + solver = builderType1Cell(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); test_analytical_robertson(solver, 2e-1, copy_to_device, copy_to_host); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); + solver = builderType1Cell(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); test_analytical_robertson(solver, 2e-1, copy_to_device, copy_to_host); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); + solver = builderType1Cell(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_robertson(solver, 2e-1, copy_to_device, copy_to_host); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); + solver = builderType1Cell(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_robertson(solver, 2e-1, copy_to_device, copy_to_host); } @@ -174,75 +169,56 @@ TEST(AnalyticalExamplesCudaRosenbrock, SurfaceRxn) test_analytical_surface_rxn(six_da_1_cell, 1e-7, copy_to_device, copy_to_host); } -//Commented out the following tests as it is failing to copy State to Device -//Will uncomment them once we have a plan how to deal with state (relative tolearnce and absolute tolerance) on GPU -/* TEST(AnalyticalExamplesCudaRosenbrock, E5) { - auto rosenbrock_solver = [](auto params) - { - return builderType1Cell(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); + auto solver = builderType1Cell(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); test_analytical_e5(solver, 1e-3, copy_to_device, copy_to_host); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + solver = builderType1Cell(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); test_analytical_e5(solver, 1e-3, copy_to_device, copy_to_host); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); + solver = builderType1Cell(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); test_analytical_e5(solver, 1e-3, copy_to_device, copy_to_host); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); + solver = builderType1Cell(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_e5(solver, 1e-3, copy_to_device, copy_to_host); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); + solver = builderType1Cell(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_e5(solver, 1e-3, copy_to_device, copy_to_host); } TEST(AnalyticalExamplesCudaRosenbrock, Oregonator) { - auto rosenbrock_solver = [](auto params) - { - return builderType1Cell(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); + auto solver = builderType1Cell(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); test_analytical_oregonator(solver, 2e-3, copy_to_device, copy_to_host); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + solver = builderType1Cell(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); test_analytical_oregonator(solver, 2e-3, copy_to_device, copy_to_host); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); + solver = builderType1Cell(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); test_analytical_oregonator(solver, 2e-3, copy_to_device, copy_to_host); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); + solver = builderType1Cell(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_oregonator(solver, 2e-3, copy_to_device, copy_to_host); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); + solver = builderType1Cell(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_oregonator(solver, 2e-3, copy_to_device, copy_to_host); } TEST(AnalyticalExamplesCudaRosenbrock, HIRES) { - auto rosenbrock_solver = [](auto params) - { - return builderType1Cell(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); + auto solver = builderType1Cell(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); test_analytical_hires(solver, 1e-6, copy_to_device, copy_to_host); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + solver = builderType1Cell(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); test_analytical_hires(solver, 1e-7, copy_to_device, copy_to_host); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); + solver = builderType1Cell(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); test_analytical_hires(solver, 1e-7, copy_to_device, copy_to_host); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); + solver = builderType1Cell(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_hires(solver, 1e-6, copy_to_device, copy_to_host); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); + solver = builderType1Cell(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_hires(solver, 1e-6, copy_to_device, copy_to_host); } -*/ diff --git a/test/integration/jit/test_jit_analytical_rosenbrock.cpp b/test/integration/jit/test_jit_analytical_rosenbrock.cpp index 3792563da..a3dfb903c 100644 --- a/test/integration/jit/test_jit_analytical_rosenbrock.cpp +++ b/test/integration/jit/test_jit_analytical_rosenbrock.cpp @@ -146,24 +146,19 @@ TEST(AnalyticalExamplesJitRosenbrock, BranchedSuperStiffButAnalytical) TEST(AnalyticalExamplesJitRosenbrock, Robertson) { - auto rosenbrock_solver = [](auto params) - { - return BuilderType<1>(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); + auto solver = BuilderType<1>(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); test_analytical_robertson, StateType<1>>(solver, 1e-1); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + solver = BuilderType<1>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); test_analytical_robertson, StateType<1>>(solver, 1e-1); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); + solver = BuilderType<1>(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); test_analytical_robertson, StateType<1>>(solver, 1e-1); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); + solver = BuilderType<1>(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_robertson, StateType<1>>(solver, 1e-1); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); + solver = BuilderType<1>(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_robertson, StateType<1>>(solver, 1e-1); } @@ -178,69 +173,54 @@ TEST(AnalyticalExamplesJitRosenbrock, SurfaceRxn) TEST(AnalyticalExamplesJitRosenbrock, E5) { - auto rosenbrock_solver = [](auto params) - { - return BuilderType<1>(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); + auto solver = BuilderType<1>(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); test_analytical_e5, StateType<1>>(solver, 1e-3); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + solver = BuilderType<1>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); test_analytical_e5, StateType<1>>(solver, 1e-3); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); + solver = BuilderType<1>(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); test_analytical_e5, StateType<1>>(solver, 1e-3); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); + solver = BuilderType<1>(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_e5, StateType<1>>(solver, 1e-3); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); + solver = BuilderType<1>(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_e5, StateType<1>>(solver, 1e-3); } TEST(AnalyticalExamplesJitRosenbrock, Oregonator) { - auto rosenbrock_solver = [](auto params) - { - return BuilderType<1>(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); + auto solver = BuilderType<1>(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); test_analytical_oregonator, StateType<1>>(solver, 4e-4); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + solver = BuilderType<1>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); test_analytical_oregonator, StateType<1>>(solver, 4e-4); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); + solver = BuilderType<1>(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); test_analytical_oregonator, StateType<1>>(solver, 4e-4); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); + solver = BuilderType<1>(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_oregonator, StateType<1>>(solver, 4e-4); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); + solver = BuilderType<1>(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_oregonator, StateType<1>>(solver, 4e-4); } TEST(AnalyticalExamplesJitRosenbrock, HIRES) { - auto rosenbrock_solver = [](auto params) - { - return BuilderType<1>(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); + auto solver = BuilderType<1>(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); test_analytical_hires, StateType<1>>(solver, 1e-6); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + solver = BuilderType<1>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); test_analytical_hires, StateType<1>>(solver, 1e-7); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); + solver = BuilderType<1>(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); test_analytical_hires, StateType<1>>(solver, 1e-7); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); + solver = BuilderType<1>(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_hires, StateType<1>>(solver, 1e-6); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); + solver = BuilderType<1>(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_hires, StateType<1>>(solver, 1e-6); } \ No newline at end of file diff --git a/test/integration/test_analytical_rosenbrock.cpp b/test/integration/test_analytical_rosenbrock.cpp index cba463545..cd6a1bd91 100644 --- a/test/integration/test_analytical_rosenbrock.cpp +++ b/test/integration/test_analytical_rosenbrock.cpp @@ -196,70 +196,55 @@ TEST(AnalyticalExamples, BranchedSuperStiffButAnalytical) TEST(AnalyticalExamples, Robertson) { - auto rosenbrock_solver = [](auto params) - { - return BuilderType(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); + auto solver = BuilderType(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); test_analytical_robertson(solver, 1e-1); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + solver = BuilderType(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); test_analytical_robertson(solver, 1e-1); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); + solver = BuilderType(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); test_analytical_robertson(solver, 1e-1); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); + solver = BuilderType(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_robertson(solver, 1e-1); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); + solver = BuilderType(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_robertson(solver, 1e-1); } TEST(AnalyticalExamples, E5) { - auto rosenbrock_solver = [](auto params) - { - return BuilderType(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); + auto solver = BuilderType(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); test_analytical_e5(solver, 1e-3); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + solver = BuilderType(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); test_analytical_e5(solver, 1e-3); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); + solver = BuilderType(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); test_analytical_e5(solver, 1e-3); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); + solver = BuilderType(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_e5(solver, 1e-3); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); + solver = BuilderType(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_e5(solver, 1e-3); } TEST(AnalyticalExamples, Oregonator) { - auto rosenbrock_solver = [](auto params) - { - return micm::CpuSolverBuilder(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); + auto solver = BuilderType(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); test_analytical_oregonator(solver, 4e-6); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + solver = BuilderType(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); test_analytical_oregonator(solver, 4e-6); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); + solver = BuilderType(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); test_analytical_oregonator(solver, 4e-6); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); + solver = BuilderType(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_oregonator(solver, 4e-6); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); + solver = BuilderType(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_oregonator(solver, 4e-6); } @@ -274,23 +259,18 @@ TEST(AnalyticalExamples, SurfaceRxn) TEST(AnalyticalExamples, HIRES) { - auto rosenbrock_solver = [](auto params) - { - return micm::CpuSolverBuilder(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); + auto solver = BuilderType(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); test_analytical_hires(solver, 1e-6); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); + solver = BuilderType(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); test_analytical_hires(solver, 1e-7); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); + solver = BuilderType(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); test_analytical_hires(solver, 1e-7); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); + solver = BuilderType(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_hires(solver, 1e-6); - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); + solver = BuilderType(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); test_analytical_hires(solver, 1e-6); } From f3c71f79809d8f55ec5050e34ea401171a5db60a Mon Sep 17 00:00:00 2001 From: Kyle Shores Date: Mon, 11 Nov 2024 07:55:27 -0700 Subject: [PATCH 25/31] removing comment --- include/micm/solver/backward_euler.inl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/include/micm/solver/backward_euler.inl b/include/micm/solver/backward_euler.inl index fc6fa479d..f39b0deaf 100644 --- a/include/micm/solver/backward_euler.inl +++ b/include/micm/solver/backward_euler.inl @@ -182,11 +182,6 @@ namespace micm const DenseMatrixPolicy& residual, const DenseMatrixPolicy& Yn1, std::vector& absolute_tolerance, double relative_tolerance) requires(!VectorizableDense) { - std::cout << "rel_tol: " << relative_tolerance << std::endl; - for(auto& val : absolute_tolerance) - { - std::cout << "abs_tol: " << val << std::endl; - } double small = parameters.small_; double rel_tol = relative_tolerance; auto& abs_tol = absolute_tolerance; From 59460c6302ec7e1343d4eb9629770003e9db378c Mon Sep 17 00:00:00 2001 From: Kyle Shores Date: Mon, 11 Nov 2024 09:16:13 -0600 Subject: [PATCH 26/31] using same calling convention for analytical tests --- .../jit/test_jit_analytical_rosenbrock.cpp | 86 +++++------------ .../test_analytical_rosenbrock.cpp | 94 ++++++------------- 2 files changed, 54 insertions(+), 126 deletions(-) diff --git a/test/integration/jit/test_jit_analytical_rosenbrock.cpp b/test/integration/jit/test_jit_analytical_rosenbrock.cpp index a3dfb903c..a099ee351 100644 --- a/test/integration/jit/test_jit_analytical_rosenbrock.cpp +++ b/test/integration/jit/test_jit_analytical_rosenbrock.cpp @@ -144,24 +144,6 @@ TEST(AnalyticalExamplesJitRosenbrock, BranchedSuperStiffButAnalytical) test_analytical_stiff_branched, StateType>(BuilderType(param_six_da), 2e-3); } -TEST(AnalyticalExamplesJitRosenbrock, Robertson) -{ - auto solver = BuilderType<1>(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_robertson, StateType<1>>(solver, 1e-1); - - solver = BuilderType<1>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_robertson, StateType<1>>(solver, 1e-1); - - solver = BuilderType<1>(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_robertson, StateType<1>>(solver, 1e-1); - - solver = BuilderType<1>(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_robertson, StateType<1>>(solver, 1e-1); - - solver = BuilderType<1>(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_robertson, StateType<1>>(solver, 1e-1); -} - TEST(AnalyticalExamplesJitRosenbrock, SurfaceRxn) { test_analytical_surface_rxn, StateType<1>>(two, 1e-4); @@ -171,56 +153,38 @@ TEST(AnalyticalExamplesJitRosenbrock, SurfaceRxn) test_analytical_surface_rxn, StateType<1>>(six_da, 1e-5); } -TEST(AnalyticalExamplesJitRosenbrock, E5) +TEST(AnalyticalExamplesJitRosenbrock, Robertson) { - auto solver = BuilderType<1>(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_e5, StateType<1>>(solver, 1e-3); - - solver = BuilderType<1>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_e5, StateType<1>>(solver, 1e-3); - - solver = BuilderType<1>(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_e5, StateType<1>>(solver, 1e-3); - - solver = BuilderType<1>(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_e5, StateType<1>>(solver, 1e-3); + test_analytical_robertson, StateType<1>>(BuilderType<1>(param_two), 1e-1); + test_analytical_robertson, StateType<1>>(BuilderType<1>(param_three), 1e-1); + test_analytical_robertson, StateType<1>>(BuilderType<1>(param_four), 1e-1); + test_analytical_robertson, StateType<1>>(BuilderType<1>(param_four_da), 1e-1); + test_analytical_robertson, StateType<1>>(BuilderType<1>(param_six_da), 1e-1); +} - solver = BuilderType<1>(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_e5, StateType<1>>(solver, 1e-3); +TEST(AnalyticalExamplesJitRosenbrock, E5) +{ + test_analytical_e5, StateType<1>>(BuilderType<1>(param_two), 1e-3); + test_analytical_e5, StateType<1>>(BuilderType<1>(param_three), 1e-3); + test_analytical_e5, StateType<1>>(BuilderType<1>(param_four), 1e-3); + test_analytical_e5, StateType<1>>(BuilderType<1>(param_four_da), 1e-3); + test_analytical_e5, StateType<1>>(BuilderType<1>(param_six_da), 1e-3); } TEST(AnalyticalExamplesJitRosenbrock, Oregonator) { - auto solver = BuilderType<1>(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_oregonator, StateType<1>>(solver, 4e-4); - - solver = BuilderType<1>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_oregonator, StateType<1>>(solver, 4e-4); - - solver = BuilderType<1>(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_oregonator, StateType<1>>(solver, 4e-4); - - solver = BuilderType<1>(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_oregonator, StateType<1>>(solver, 4e-4); - - solver = BuilderType<1>(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_oregonator, StateType<1>>(solver, 4e-4); + test_analytical_oregonator, StateType<1>>(BuilderType<1>(param_two), 4e-4); + test_analytical_oregonator, StateType<1>>(BuilderType<1>(param_three), 4e-4); + test_analytical_oregonator, StateType<1>>(BuilderType<1>(param_four), 4e-4); + test_analytical_oregonator, StateType<1>>(BuilderType<1>(param_four_da), 4e-4); + test_analytical_oregonator, StateType<1>>(BuilderType<1>(param_six_da), 4e-4); } TEST(AnalyticalExamplesJitRosenbrock, HIRES) { - auto solver = BuilderType<1>(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_hires, StateType<1>>(solver, 1e-6); - - solver = BuilderType<1>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_hires, StateType<1>>(solver, 1e-7); - - solver = BuilderType<1>(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_hires, StateType<1>>(solver, 1e-7); - - solver = BuilderType<1>(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_hires, StateType<1>>(solver, 1e-6); - - solver = BuilderType<1>(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_hires, StateType<1>>(solver, 1e-6); -} \ No newline at end of file + test_analytical_hires, StateType<1>>(BuilderType<1>(param_two), 1e-6); + test_analytical_hires, StateType<1>>(BuilderType<1>(param_three), 1e-7); + test_analytical_hires, StateType<1>>(BuilderType<1>(param_four), 1e-7); + test_analytical_hires, StateType<1>>(BuilderType<1>(param_four_da), 1e-6); + test_analytical_hires, StateType<1>>(BuilderType<1>(param_six_da), 1e-6); +} diff --git a/test/integration/test_analytical_rosenbrock.cpp b/test/integration/test_analytical_rosenbrock.cpp index cd6a1bd91..2bd3ca1fa 100644 --- a/test/integration/test_analytical_rosenbrock.cpp +++ b/test/integration/test_analytical_rosenbrock.cpp @@ -194,83 +194,47 @@ TEST(AnalyticalExamples, BranchedSuperStiffButAnalytical) test_analytical_stiff_branched, VectorStateType<4>>(rosenbrock_vector_4, 2e-3); } -TEST(AnalyticalExamples, Robertson) +TEST(AnalyticalExamples, SurfaceRxn) { - auto solver = BuilderType(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_robertson(solver, 1e-1); - - solver = BuilderType(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_robertson(solver, 1e-1); - - solver = BuilderType(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_robertson(solver, 1e-1); - - solver = BuilderType(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_robertson(solver, 1e-1); - - solver = BuilderType(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_robertson(solver, 1e-1); + test_analytical_surface_rxn(rosenbrock_2stage, 1e-2); + test_analytical_surface_rxn(rosenbrock_3stage, 1e-5); + test_analytical_surface_rxn(rosenbrock_4stage, 1e-6); + test_analytical_surface_rxn(rosenbrock_4stage_da, 1e-5); + test_analytical_surface_rxn(rosenbrock_6stage_da, 1e-7); } -TEST(AnalyticalExamples, E5) +TEST(AnalyticalExamples, Robertson) { - auto solver = BuilderType(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_e5(solver, 1e-3); - - solver = BuilderType(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_e5(solver, 1e-3); - - solver = BuilderType(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_e5(solver, 1e-3); - - solver = BuilderType(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_e5(solver, 1e-3); - - solver = BuilderType(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_e5(solver, 1e-3); + test_analytical_robertson(rosenbrock_2stage, 1e-1); + test_analytical_robertson(rosenbrock_3stage, 1e-1); + test_analytical_robertson(rosenbrock_4stage, 1e-1); + test_analytical_robertson(rosenbrock_4stage_da, 1e-1); + test_analytical_robertson(rosenbrock_6stage_da, 1e-1); } -TEST(AnalyticalExamples, Oregonator) +TEST(AnalyticalExamples, E5) { - auto solver = BuilderType(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_oregonator(solver, 4e-6); - - solver = BuilderType(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_oregonator(solver, 4e-6); - - solver = BuilderType(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_oregonator(solver, 4e-6); - - solver = BuilderType(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_oregonator(solver, 4e-6); - - solver = BuilderType(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_oregonator(solver, 4e-6); + test_analytical_e5(rosenbrock_2stage, 1e-3); + test_analytical_e5(rosenbrock_3stage, 1e-3); + test_analytical_e5(rosenbrock_4stage, 1e-3); + test_analytical_e5(rosenbrock_4stage_da, 1e-3); + test_analytical_e5(rosenbrock_6stage_da, 1e-3); } -TEST(AnalyticalExamples, SurfaceRxn) +TEST(AnalyticalExamples, Oregonator) { - test_analytical_surface_rxn(rosenbrock_2stage, 1e-2); - test_analytical_surface_rxn(rosenbrock_3stage, 1e-5); - test_analytical_surface_rxn(rosenbrock_4stage, 1e-6); - test_analytical_surface_rxn(rosenbrock_4stage_da, 1e-5); - test_analytical_surface_rxn(rosenbrock_6stage_da, 1e-7); + test_analytical_oregonator(rosenbrock_2stage, 4e-6); + test_analytical_oregonator(rosenbrock_3stage, 4e-6); + test_analytical_oregonator(rosenbrock_4stage, 4e-6); + test_analytical_oregonator(rosenbrock_4stage_da, 4e-6); + test_analytical_oregonator(rosenbrock_6stage_da, 4e-6); } TEST(AnalyticalExamples, HIRES) { - auto solver = BuilderType(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_hires(solver, 1e-6); - - solver = BuilderType(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_hires(solver, 1e-7); - - solver = BuilderType(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_hires(solver, 1e-7); - - solver = BuilderType(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_hires(solver, 1e-6); - - solver = BuilderType(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_hires(solver, 1e-6); + test_analytical_hires(rosenbrock_2stage, 1e-6); + test_analytical_hires(rosenbrock_3stage, 1e-7); + test_analytical_hires(rosenbrock_4stage, 1e-7); + test_analytical_hires(rosenbrock_4stage_da, 1e-6); + test_analytical_hires(rosenbrock_6stage_da, 1e-6); } From 4b81426b43b44aaf5f23f5ec2e9f1d2e418ae1e5 Mon Sep 17 00:00:00 2001 From: Kyle Shores Date: Mon, 11 Nov 2024 10:00:57 -0700 Subject: [PATCH 27/31] simplifying cuda analytical test interface --- .../cuda/test_cuda_analytical_rosenbrock.cpp | 84 ++++++------------- 1 file changed, 24 insertions(+), 60 deletions(-) diff --git a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp index 28fbdc1e0..d4fb6e266 100644 --- a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp +++ b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp @@ -142,24 +142,6 @@ TEST(AnalyticalExamplesCudaRosenbrock, BranchedSuperStiffButAnalytical) test_analytical_stiff_branched(six_da, 2e-3, copy_to_device, copy_to_host); } -TEST(AnalyticalExamplesCudaRosenbrock, Robertson) -{ - auto solver = builderType1Cell(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_robertson(solver, 2e-1, copy_to_device, copy_to_host); - - solver = builderType1Cell(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_robertson(solver, 2e-1, copy_to_device, copy_to_host); - - solver = builderType1Cell(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_robertson(solver, 2e-1, copy_to_device, copy_to_host); - - solver = builderType1Cell(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_robertson(solver, 2e-1, copy_to_device, copy_to_host); - - solver = builderType1Cell(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_robertson(solver, 2e-1, copy_to_device, copy_to_host); -} - TEST(AnalyticalExamplesCudaRosenbrock, SurfaceRxn) { test_analytical_surface_rxn(two_1_cell, 1e-2, copy_to_device, copy_to_host); @@ -169,56 +151,38 @@ TEST(AnalyticalExamplesCudaRosenbrock, SurfaceRxn) test_analytical_surface_rxn(six_da_1_cell, 1e-7, copy_to_device, copy_to_host); } -TEST(AnalyticalExamplesCudaRosenbrock, E5) +TEST(AnalyticalExamplesCudaRosenbrock, Robertson) { - auto solver = builderType1Cell(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_e5(solver, 1e-3, copy_to_device, copy_to_host); - - solver = builderType1Cell(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_e5(solver, 1e-3, copy_to_device, copy_to_host); - - solver = builderType1Cell(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_e5(solver, 1e-3, copy_to_device, copy_to_host); - - solver = builderType1Cell(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_e5(solver, 1e-3, copy_to_device, copy_to_host); + test_analytical_robertson(two_1_cell, 2e-1, copy_to_device, copy_to_host); + test_analytical_robertson(three_1_cell, 2e-1, copy_to_device, copy_to_host); + test_analytical_robertson(four_1_cell, 2e-1, copy_to_device, copy_to_host); + test_analytical_robertson(four_da_1_cell, 2e-1, copy_to_device, copy_to_host); + test_analytical_robertson(six_da_1_cell, 2e-1, copy_to_device, copy_to_host); +} - solver = builderType1Cell(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_e5(solver, 1e-3, copy_to_device, copy_to_host); +TEST(AnalyticalExamplesCudaRosenbrock, E5) +{ + test_analytical_e5(two_1_cell, 1e-3, copy_to_device, copy_to_host); + test_analytical_e5(three_1_cell, 1e-3, copy_to_device, copy_to_host); + test_analytical_e5(four_1_cell, 1e-3, copy_to_device, copy_to_host); + test_analytical_e5(four_da_1_cell, 1e-3, copy_to_device, copy_to_host); + test_analytical_e5(six_da_1_cell, 1e-3, copy_to_device, copy_to_host); } TEST(AnalyticalExamplesCudaRosenbrock, Oregonator) { - auto solver = builderType1Cell(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_oregonator(solver, 2e-3, copy_to_device, copy_to_host); - - solver = builderType1Cell(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_oregonator(solver, 2e-3, copy_to_device, copy_to_host); - - solver = builderType1Cell(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_oregonator(solver, 2e-3, copy_to_device, copy_to_host); - - solver = builderType1Cell(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_oregonator(solver, 2e-3, copy_to_device, copy_to_host); - - solver = builderType1Cell(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_oregonator(solver, 2e-3, copy_to_device, copy_to_host); + test_analytical_oregonator(two_1_cell, 2e-3, copy_to_device, copy_to_host); + test_analytical_oregonator(three_1_cell, 2e-3, copy_to_device, copy_to_host); + test_analytical_oregonator(four_1_cell, 2e-3, copy_to_device, copy_to_host); + test_analytical_oregonator(four_da_1_cell, 2e-3, copy_to_device, copy_to_host); + test_analytical_oregonator(six_da_1_cell, 2e-3, copy_to_device, copy_to_host); } TEST(AnalyticalExamplesCudaRosenbrock, HIRES) { - auto solver = builderType1Cell(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_hires(solver, 1e-6, copy_to_device, copy_to_host); - - solver = builderType1Cell(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_hires(solver, 1e-7, copy_to_device, copy_to_host); - - solver = builderType1Cell(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_hires(solver, 1e-7, copy_to_device, copy_to_host); - - solver = builderType1Cell(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_hires(solver, 1e-6, copy_to_device, copy_to_host); - - solver = builderType1Cell(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_hires(solver, 1e-6, copy_to_device, copy_to_host); + test_analytical_hires(two_1_cell, 1e-6, copy_to_device, copy_to_host); + test_analytical_hires(three_1_cell, 1e-7, copy_to_device, copy_to_host); + test_analytical_hires(four_1_cell, 1e-7, copy_to_device, copy_to_host); + test_analytical_hires(four_da_1_cell, 1e-6, copy_to_device, copy_to_host); + test_analytical_hires(six_da_1_cell, 1e-6, copy_to_device, copy_to_host); } From 6112c2cae1cd2fd45d04d1c23620dd18d9cfebf1 Mon Sep 17 00:00:00 2001 From: Montek Thind Date: Wed, 13 Nov 2024 11:30:57 -0800 Subject: [PATCH 28/31] added tests --- include/micm/cuda/solver/cuda_state.hpp | 2 - test/unit/solver/test_solver_builder.cpp | 58 ++++++++++++++++++++++++ 2 files changed, 58 insertions(+), 2 deletions(-) diff --git a/include/micm/cuda/solver/cuda_state.hpp b/include/micm/cuda/solver/cuda_state.hpp index 367dde531..8c4ea2156 100644 --- a/include/micm/cuda/solver/cuda_state.hpp +++ b/include/micm/cuda/solver/cuda_state.hpp @@ -7,8 +7,6 @@ #include #include -#include - namespace micm { /// @brief Construct a state variable for CUDA tests diff --git a/test/unit/solver/test_solver_builder.cpp b/test/unit/solver/test_solver_builder.cpp index 047c0972c..04b55c23c 100644 --- a/test/unit/solver/test_solver_builder.cpp +++ b/test/unit/solver/test_solver_builder.cpp @@ -104,4 +104,62 @@ TEST(SolverBuilder, CanBuildRosenbrock) .SetReactions(reactions) .SetNumberOfGridCells(1) .Build(); +} + +TEST(SolverBuilder, CanBuildBackwardEulerOverloadedSolverMethod) +{ + auto solver = micm::CpuSolverBuilder(micm::BackwardEulerSolverParameters{}) + .SetSystem(the_system) + .SetReactions(reactions) + .SetNumberOfGridCells(1) + .Build(); + auto state = solver.GetState(); + auto options = micm::BackwardEulerSolverParameters(); + auto solve = solver.Solve(5, state, options); + + ASSERT_EQ(solve.final_time_, 5); + ASSERT_EQ(solve.stats_.function_calls_, 2); + ASSERT_EQ(solve.stats_.jacobian_updates_, 2); + ASSERT_EQ(solve.stats_.number_of_steps_, 2); + ASSERT_EQ(solve.stats_.solves_, 2); + + options.small_ = 1.0; + options.max_number_of_steps_ = 1.0; + + solve = solver.Solve(5, state, options); + + ASSERT_EQ(solve.final_time_, 0.03125); + ASSERT_EQ(solve.stats_.function_calls_, 6); + ASSERT_EQ(solve.stats_.jacobian_updates_, 6); + ASSERT_EQ(solve.stats_.number_of_steps_, 6); + ASSERT_EQ(solve.stats_.solves_, 6); +} + +TEST(SolverBuilder, CanBuildRosenbrockOverloadedSolveMethod) +{ + auto solver = micm::CpuSolverBuilder(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) + .SetSystem(the_system) + .SetReactions(reactions) + .SetNumberOfGridCells(1) + .Build(); + auto state = solver.GetState(); + auto options = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); + auto solve = solver.Solve(5, state, options); + + ASSERT_EQ(solve.final_time_, 5); + ASSERT_EQ(solve.stats_.function_calls_, 20); + ASSERT_EQ(solve.stats_.jacobian_updates_, 10); + ASSERT_EQ(solve.stats_.number_of_steps_, 10); + ASSERT_EQ(solve.stats_.solves_, 30); + + options.h_min_ = 15.0; + options.max_number_of_steps_ = 6.0; + + solve = solver.Solve(5, state, options); + + ASSERT_EQ(solve.final_time_, 5); + ASSERT_EQ(solve.stats_.function_calls_, 2); + ASSERT_EQ(solve.stats_.jacobian_updates_, 1); + ASSERT_EQ(solve.stats_.number_of_steps_, 1); + ASSERT_EQ(solve.stats_.solves_, 3); } \ No newline at end of file From bd1829e0e5b8972da51a659806cb04fa56308efc Mon Sep 17 00:00:00 2001 From: Jian Sun Date: Fri, 15 Nov 2024 15:21:50 -0700 Subject: [PATCH 29/31] move the copy of absolute tolerance into the cuda state constructor --- include/micm/cuda/solver/cuda_rosenbrock.hpp | 2 +- include/micm/cuda/solver/cuda_state.hpp | 46 +++++++++++++++---- include/micm/cuda/util/cuda_matrix.cuh | 2 +- include/micm/solver/backward_euler.hpp | 4 +- include/micm/solver/backward_euler.inl | 6 +-- include/micm/solver/rosenbrock.inl | 18 +++++--- include/micm/solver/state.hpp | 20 ++++++-- include/micm/solver/state.inl | 20 +++++++- src/util/cuda_matrix.cu | 4 +- test/integration/analytical_policy.hpp | 24 +++++----- .../cuda/test_cuda_analytical_rosenbrock.cpp | 16 +++---- test/integration/terminator.hpp | 2 +- test/integration/test_chapman_integration.cpp | 6 +-- .../unit/cuda/solver/test_cuda_rosenbrock.cpp | 8 ++-- test/unit/solver/test_backward_euler.cpp | 2 +- test/unit/solver/test_rosenbrock.cpp | 11 +++-- 16 files changed, 128 insertions(+), 63 deletions(-) diff --git a/include/micm/cuda/solver/cuda_rosenbrock.hpp b/include/micm/cuda/solver/cuda_rosenbrock.hpp index 36b96d27e..138e5b4e6 100644 --- a/include/micm/cuda/solver/cuda_rosenbrock.hpp +++ b/include/micm/cuda/solver/cuda_rosenbrock.hpp @@ -122,7 +122,7 @@ namespace micm y_new.AsDeviceParam(), errors.AsDeviceParam(), state.absolute_tolerance_param_, - state.relative_tolerance_, + state.GetRelativeTolerance(), this->devstruct_); } }; // end CudaRosenbrockSolver diff --git a/include/micm/cuda/solver/cuda_state.hpp b/include/micm/cuda/solver/cuda_state.hpp index 367dde531..aa218c6a4 100644 --- a/include/micm/cuda/solver/cuda_state.hpp +++ b/include/micm/cuda/solver/cuda_state.hpp @@ -18,33 +18,59 @@ 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"); - absolute_tolerance_param_.d_data_ = nullptr; } /// @brief Constructor which takes the state dimension information as input /// @param parameters State dimension information CudaState(const StateParameters& parameters) - : State(parameters){}; + : State(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(absolute_tolerance_param_, absolute_tolerance_param_.number_of_elements_), "cudaMalloc"); + CHECK_CUDA_ERROR(micm::cuda::CopyToDevice(absolute_tolerance_param_, atol), "cudaMemcpyHostToDevice"); + }; + + /// @brief Move constructor + CudaState(CudaState&& other) + : State(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::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& absoluteTolerance) override + { + State::SetAbsoluteTolerances(absoluteTolerance); + CHECK_CUDA_ERROR(micm::cuda::CopyToDevice(absolute_tolerance_param_, absoluteTolerance), "cudaMemcpyHostToDevice"); + } /// @brief Copy input variables to the device void SyncInputsToDevice() requires(CudaMatrix&& VectorizableDense) { this->variables_.CopyToDevice(); this->rate_constants_.CopyToDevice(); - - absolute_tolerance_param_.number_of_elements_ = this->absolute_tolerance_.size(); - absolute_tolerance_param_.number_of_grid_cells_ = 1; - - 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/cuda/util/cuda_matrix.cuh b/include/micm/cuda/util/cuda_matrix.cuh index 61bfff6e5..6040fc682 100644 --- a/include/micm/cuda/util/cuda_matrix.cuh +++ b/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 - cudaError_t CopyToDevice(CudaMatrixParam& vectorMatrix, std::vector& h_data); + cudaError_t CopyToDevice(CudaMatrixParam& vectorMatrix, const std::vector& h_data); /// @brief Copies data from the device to the host /// @param vectorMatrix Struct containing allocated device memory diff --git a/include/micm/solver/backward_euler.hpp b/include/micm/solver/backward_euler.hpp index c3a6c65fa..6e135cfeb 100644 --- a/include/micm/solver/backward_euler.hpp +++ b/include/micm/solver/backward_euler.hpp @@ -81,12 +81,12 @@ namespace micm static bool IsConverged( const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& Yn1, std::vector& absolute_tolerance, double relative_tolerance) requires(!VectorizableDense); + const DenseMatrixPolicy& Yn1, const std::vector& absolute_tolerance, double relative_tolerance) requires(!VectorizableDense); template static bool IsConverged( const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& Yn1, std::vector& absolute_tolerance, double relative_tolerance) requires(VectorizableDense); + const DenseMatrixPolicy& Yn1, const std::vector& absolute_tolerance, double relative_tolerance) requires(VectorizableDense); }; } // namespace micm diff --git a/include/micm/solver/backward_euler.inl b/include/micm/solver/backward_euler.inl index f39b0deaf..dc41ab340 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_, state.relative_tolerance_); + converged = IsConverged(parameters_, forcing, Yn1, state.GetAbsoluteTolerances(), state.GetRelativeTolerance()); } while (!converged && iterations < max_iter); if (!converged) @@ -180,7 +180,7 @@ namespace micm inline bool BackwardEuler::IsConverged( const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& Yn1, std::vector& absolute_tolerance, double relative_tolerance) requires(!VectorizableDense) + const DenseMatrixPolicy& Yn1, const std::vector& absolute_tolerance, double relative_tolerance) requires(!VectorizableDense) { double small = parameters.small_; double rel_tol = relative_tolerance; @@ -206,7 +206,7 @@ namespace micm inline bool BackwardEuler::IsConverged( const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& Yn1, std::vector& absolute_tolerance, double relative_tolerance) requires(VectorizableDense) + const DenseMatrixPolicy& Yn1, const std::vector& absolute_tolerance, double relative_tolerance) requires(VectorizableDense) { double small = parameters.small_; double rel_tol = relative_tolerance; diff --git a/include/micm/solver/rosenbrock.inl b/include/micm/solver/rosenbrock.inl index 2e6484b03..78818c752 100644 --- a/include/micm/solver/rosenbrock.inl +++ b/include/micm/solver/rosenbrock.inl @@ -253,8 +253,10 @@ namespace micm auto& _y = Y.AsVector(); auto& _ynew = Ynew.AsVector(); auto& _errors = errors.AsVector(); + const auto& atol = state.GetAbsoluteTolerances(); + const auto& rtol = state.GetRelativeTolerance(); const std::size_t N = Y.AsVector().size(); - const std::size_t n_vars = state.absolute_tolerance_.size(); + const std::size_t n_vars = atol.size(); double ymax = 0; double errors_over_scale = 0; @@ -264,7 +266,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] / (atol[i % n_vars] + rtol * ymax); error += errors_over_scale * errors_over_scale; } @@ -289,9 +291,11 @@ namespace micm auto y_iter = Y.AsVector().begin(); auto ynew_iter = Ynew.AsVector().begin(); auto errors_iter = errors.AsVector().begin(); + const auto& atol = state.GetAbsoluteTolerances(); + auto rtol = state.GetRelativeTolerance(); 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 = atol.size(); const std::size_t whole_blocks = std::floor(Y.NumRows() / Y.GroupVectorSize()) * Y.GroupSize(); @@ -302,8 +306,8 @@ namespace micm for (std::size_t i = 0; i < whole_blocks; ++i) { errors_over_scale = - *errors_iter / (state.absolute_tolerance_[(i / L) % n_vars] + - state.relative_tolerance_ * std::max(std::abs(*y_iter), std::abs(*ynew_iter))); + *errors_iter / (atol[(i / L) % n_vars] + + rtol * std::max(std::abs(*y_iter), std::abs(*ynew_iter))); error += errors_over_scale * errors_over_scale; ++y_iter; ++ynew_iter; @@ -321,8 +325,8 @@ namespace micm { const std::size_t idx = y * L + x; errors_over_scale = errors_iter[idx] / - (state.absolute_tolerance_[y] + - state.relative_tolerance_ * std::max(std::abs(y_iter[idx]), std::abs(ynew_iter[idx]))); + (atol[y] + + rtol * 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 082edde4b..ac7914819 100644 --- a/include/micm/solver/state.hpp +++ b/include/micm/solver/state.hpp @@ -60,9 +60,12 @@ namespace micm std::size_t state_size_; std::size_t number_of_grid_cells_; std::unique_ptr temporary_variables_; - double relative_tolerance_; - std::vector absolute_tolerance_; + private: + double relative_tolerance_; + std::vector absolute_tolerance_; + + public: /// @brief Default constructor /// Only defined to be used to create default values in types, but a default constructed state is not useable State(); @@ -191,8 +194,19 @@ namespace micm /// @param value new parameter value void SetCustomRateParameter(const std::string& label, double value); void SetCustomRateParameter(const std::string& label, const std::vector& values); - void SetRelativeTolerances(double relativeTolerance); + /// @brief Set the relative tolerances + /// @param relativeTolerance relative tolerance + void SetRelativeTolerance(double relativeTolerance); + + /// @brief Set the absolute tolerances per species + /// @param absoluteTolerance absolute tolerance + virtual void SetAbsoluteTolerances(const std::vector& absoluteTolerance); + + double GetRelativeTolerance() const; + + const std::vector& GetAbsoluteTolerances() const; + /// @brief Print a header of species to display concentrations with respect to time void PrintHeader(); diff --git a/include/micm/solver/state.inl b/include/micm/solver/state.inl index 6bc63897e..5d76ee50c 100644 --- a/include/micm/solver/state.inl +++ b/include/micm/solver/state.inl @@ -189,11 +189,29 @@ namespace micm } template - inline void State::SetRelativeTolerances(double relativeTolerance) + inline void State::SetRelativeTolerance(double relativeTolerance) { this->relative_tolerance_ = relativeTolerance; } + template + inline void State::SetAbsoluteTolerances(const std::vector& absoluteTolerance) + { + this->absolute_tolerance_ = absoluteTolerance; + } + + template + inline double State::GetRelativeTolerance() const + { + return this->relative_tolerance_; + } + + template + inline const std::vector& State::GetAbsoluteTolerances() const + { + return this->absolute_tolerance_; + } + template inline void State::PrintHeader() { diff --git a/src/util/cuda_matrix.cu b/src/util/cuda_matrix.cu index 08a3946bf..fd5ac1948 100644 --- a/src/util/cuda_matrix.cu +++ b/src/util/cuda_matrix.cu @@ -36,7 +36,7 @@ namespace micm } template - cudaError_t CopyToDevice(CudaMatrixParam& param, std::vector& h_data) + cudaError_t CopyToDevice(CudaMatrixParam& param, const std::vector& h_data) { cudaError_t err = cudaMemcpyAsync( param.d_data_, @@ -98,7 +98,7 @@ namespace micm // source code needs the instantiation of the template template cudaError_t MallocVector(CudaMatrixParam& param, std::size_t number_of_elements); template cudaError_t MallocVector(CudaMatrixParam& param, std::size_t number_of_elements); - template cudaError_t CopyToDevice(CudaMatrixParam& param, std::vector& h_data); + template cudaError_t CopyToDevice(CudaMatrixParam& param, const std::vector& h_data); template cudaError_t CopyToHost(CudaMatrixParam& param, std::vector& h_data); template cudaError_t CopyToDeviceFromDevice( CudaMatrixParam& vectorMatrixDest, diff --git a/test/integration/analytical_policy.hpp b/test/integration/analytical_policy.hpp index 44e50c690..07dbb2014 100644 --- a/test/integration/analytical_policy.hpp +++ b/test/integration/analytical_policy.hpp @@ -1331,8 +1331,8 @@ void test_analytical_robertson( double air_density = 1e6; auto state = solver.GetState(); - state.relative_tolerance_ = 1e-10; - state.absolute_tolerance_ = std::vector(3, state.relative_tolerance_ * 1e-2); + state.SetRelativeTolerance(1e-10); + state.SetAbsoluteTolerances(std::vector(3, state.GetRelativeTolerance() * 1e-2)); double k1 = 0.04; double k2 = 3e7; @@ -1547,8 +1547,8 @@ void test_analytical_oregonator( auto state = solver.GetState(); - state.relative_tolerance_ = 1e-8; - state.absolute_tolerance_ = std::vector(5, state.relative_tolerance_ * 1e-6); + state.SetRelativeTolerance(1e-8); + state.SetAbsoluteTolerances(std::vector(5, state.GetRelativeTolerance() * 1e-6)); state.SetCustomRateParameter("r1", 1.34 * 0.06); state.SetCustomRateParameter("r2", 1.6e9); @@ -1725,8 +1725,8 @@ void test_analytical_hires( }; auto state = solver.GetState(); - state.relative_tolerance_ = 1e-6; - state.absolute_tolerance_ = std::vector(8, state.relative_tolerance_ * 1e-2); + state.SetRelativeTolerance(1e-6); + state.SetAbsoluteTolerances(std::vector(8, state.GetRelativeTolerance() * 1e-2)); state.SetCustomRateParameter("r1", 1.71); state.SetCustomRateParameter("r2", 8.75); @@ -1885,11 +1885,13 @@ void test_analytical_e5( auto state = solver.GetState(); - state.relative_tolerance_ = 1e-13; - 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.SetRelativeTolerance(1e-13); + state.SetAbsoluteTolerances(std::vector(6, 1e-17)); + auto atol = state.GetAbsoluteTolerances(); + atol[0] = 1e-7; + atol[4] = 1e-7; + atol[5] = 1e-7; + state.SetAbsoluteTolerances(atol); state.SetCustomRateParameter("r1", 7.89e-10); state.SetCustomRateParameter("r2", 1.13e9); diff --git a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp index d4fb6e266..a3b6f66c0 100644 --- a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp +++ b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp @@ -178,11 +178,11 @@ TEST(AnalyticalExamplesCudaRosenbrock, Oregonator) test_analytical_oregonator(six_da_1_cell, 2e-3, copy_to_device, copy_to_host); } -TEST(AnalyticalExamplesCudaRosenbrock, HIRES) -{ - test_analytical_hires(two_1_cell, 1e-6, copy_to_device, copy_to_host); - test_analytical_hires(three_1_cell, 1e-7, copy_to_device, copy_to_host); - test_analytical_hires(four_1_cell, 1e-7, copy_to_device, copy_to_host); - test_analytical_hires(four_da_1_cell, 1e-6, copy_to_device, copy_to_host); - test_analytical_hires(six_da_1_cell, 1e-6, copy_to_device, copy_to_host); -} +// TEST(AnalyticalExamplesCudaRosenbrock, HIRES) +// { +// test_analytical_hires(two_1_cell, 1e-6, copy_to_device, copy_to_host); +// test_analytical_hires(three_1_cell, 1e-7, copy_to_device, copy_to_host); +// test_analytical_hires(four_1_cell, 1e-7, copy_to_device, copy_to_host); +// test_analytical_hires(four_da_1_cell, 1e-6, copy_to_device, copy_to_host); +// test_analytical_hires(six_da_1_cell, 1e-6, copy_to_device, copy_to_host); +// } diff --git a/test/integration/terminator.hpp b/test/integration/terminator.hpp index fe9832c6b..7ffc1329c 100644 --- a/test/integration/terminator.hpp +++ b/test/integration/terminator.hpp @@ -48,7 +48,7 @@ void TestTerminator(BuilderPolicy& builder, std::size_t number_of_grid_cells) .SetNumberOfGridCells(number_of_grid_cells) .Build(); auto state = solver.GetState(); - state.relative_tolerance_ = 1.0e-8; + state.SetRelativeTolerance(1.0e-8); auto get_double = std::bind(std::lognormal_distribution(-2.0, 2.0), std::default_random_engine()); std::unordered_map> concentrations{ { "Cl2", {} }, { "Cl", {} } }; diff --git a/test/integration/test_chapman_integration.cpp b/test/integration/test_chapman_integration.cpp index 545571a17..a9d265dca 100644 --- a/test/integration/test_chapman_integration.cpp +++ b/test/integration/test_chapman_integration.cpp @@ -41,9 +41,9 @@ TEST(ChapmanIntegration, CanBuildChapmanSystemUsingConfig) .Build(); 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_; + state.SetRelativeTolerance(solver_params.relative_tolerance_); + EXPECT_EQ(state.GetRelativeTolerance(), 1.0e-4); + auto& abs_tol = state.GetAbsoluteTolerances(); 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 eceb9d4a6..e4b838567 100644 --- a/test/unit/cuda/solver/test_cuda_rosenbrock.cpp +++ b/test/unit/cuda/solver/test_cuda_rosenbrock.cpp @@ -122,8 +122,8 @@ void testNormalizedErrorConst() gpu_builder = getSolver(gpu_builder); auto gpu_solver = gpu_builder.SetNumberOfGridCells(number_of_grid_cells).Build(); auto state = gpu_solver.GetState(); - auto& atol = state.absolute_tolerance_; - double rtol = state.relative_tolerance_; + auto& atol = state.GetAbsoluteTolerances() ; + double rtol = state.GetRelativeTolerance(); auto y_old = micm::CudaDenseMatrix(number_of_grid_cells, state.state_size_, 1.0); auto y_new = micm::CudaDenseMatrix(number_of_grid_cells, state.state_size_, 2.0); @@ -168,8 +168,8 @@ void testNormalizedErrorDiff() gpu_builder = getSolver(gpu_builder); auto gpu_solver = gpu_builder.SetNumberOfGridCells(number_of_grid_cells).Build(); auto state = gpu_solver.GetState(); - auto& atol = state.absolute_tolerance_; - double rtol = state.relative_tolerance_; + auto& atol = state.GetAbsoluteTolerances(); + double rtol = state.GetRelativeTolerance(); 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); diff --git a/test/unit/solver/test_backward_euler.cpp b/test/unit/solver/test_backward_euler.cpp index eb976b350..b507b6045 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_ = { 1e-6, 1e-6, 1e-6 }; + state.SetAbsoluteTolerances({ 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 01aa8a510..812093efb 100644 --- a/test/unit/solver/test_rosenbrock.cpp +++ b/test/unit/solver/test_rosenbrock.cpp @@ -18,8 +18,8 @@ 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_; - double rtol = state.relative_tolerance_; + const std::vector& atol = state.GetAbsoluteTolerances(); + double rtol = state.GetRelativeTolerance(); using MatrixPolicy = decltype(state.variables_); auto y_old = MatrixPolicy(number_of_grid_cells, state.state_size_, 7.7); @@ -107,9 +107,10 @@ 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); + auto absolute_tolerances = state.GetAbsoluteTolerances(); + EXPECT_EQ(absolute_tolerances.size(), 2); + EXPECT_EQ(absolute_tolerances[0], 1.0e-07); + EXPECT_EQ(absolute_tolerances[1], 1.0e-08); } } From 5a612d62c5fc9db90556768b5edf59262c29a332 Mon Sep 17 00:00:00 2001 From: Jian Sun Date: Fri, 15 Nov 2024 15:28:01 -0700 Subject: [PATCH 30/31] fix broken CI test --- examples/profile_example.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/profile_example.cpp b/examples/profile_example.cpp index 852e0f9fd..f8214e40a 100644 --- a/examples/profile_example.cpp +++ b/examples/profile_example.cpp @@ -72,7 +72,7 @@ int Run(const char* filepath, const char* initial_conditions, const std::string& .SetReactions(reactions) .Build(); auto state = solver.GetState(); - state.SetRelativeTolerances(0.1); + state.SetRelativeTolerance(0.1); state.conditions_[0].temperature_ = dataMap.environments["temperature"]; // K state.conditions_[0].pressure_ = dataMap.environments["pressure"]; // Pa From bacb5587a8be6bf90d93bf96b30840bda9847eb0 Mon Sep 17 00:00:00 2001 From: Jian Sun Date: Mon, 18 Nov 2024 14:40:49 -0700 Subject: [PATCH 31/31] uncomment a GPU test --- .../cuda/test_cuda_analytical_rosenbrock.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp index a3b6f66c0..d4fb6e266 100644 --- a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp +++ b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp @@ -178,11 +178,11 @@ TEST(AnalyticalExamplesCudaRosenbrock, Oregonator) test_analytical_oregonator(six_da_1_cell, 2e-3, copy_to_device, copy_to_host); } -// TEST(AnalyticalExamplesCudaRosenbrock, HIRES) -// { -// test_analytical_hires(two_1_cell, 1e-6, copy_to_device, copy_to_host); -// test_analytical_hires(three_1_cell, 1e-7, copy_to_device, copy_to_host); -// test_analytical_hires(four_1_cell, 1e-7, copy_to_device, copy_to_host); -// test_analytical_hires(four_da_1_cell, 1e-6, copy_to_device, copy_to_host); -// test_analytical_hires(six_da_1_cell, 1e-6, copy_to_device, copy_to_host); -// } +TEST(AnalyticalExamplesCudaRosenbrock, HIRES) +{ + test_analytical_hires(two_1_cell, 1e-6, copy_to_device, copy_to_host); + test_analytical_hires(three_1_cell, 1e-7, copy_to_device, copy_to_host); + test_analytical_hires(four_1_cell, 1e-7, copy_to_device, copy_to_host); + test_analytical_hires(four_da_1_cell, 1e-6, copy_to_device, copy_to_host); + test_analytical_hires(six_da_1_cell, 1e-6, copy_to_device, copy_to_host); +}