diff --git a/examples/profile_example.cpp b/examples/profile_example.cpp index cb206e0d2..f8214e40a 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.SetRelativeTolerance(0.1); + state.conditions_[0].temperature_ = dataMap.environments["temperature"]; // K state.conditions_[0].pressure_ = dataMap.environments["pressure"]; // Pa state.conditions_[0].air_density_ = dataMap.environments["air_density"]; // mol m-3 diff --git a/include/micm/configure/solver_config.hpp b/include/micm/configure/solver_config.hpp index c23c5a795..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/cuda/solver/cuda_rosenbrock.cuh b/include/micm/cuda/solver/cuda_rosenbrock.cuh index 3a03fb8b9..ea00f3da1 100644 --- a/include/micm/cuda/solver/cuda_rosenbrock.cuh +++ b/include/micm/cuda/solver/cuda_rosenbrock.cuh @@ -41,7 +41,8 @@ namespace micm const CudaMatrixParam& y_old_param, const CudaMatrixParam& y_new_param, const CudaMatrixParam& errors_param, - const RosenbrockSolverParameters& ros_param, + const CudaMatrixParam& absolute_tolerance_param, + const double relative_tolerance, CudaRosenbrockSolverParam devstruct); } // namespace cuda } // namespace micm \ 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 2044f97c9..138e5b4e6 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; }; @@ -69,21 +68,21 @@ 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(); - hoststruct.absolute_tolerance_ = this->parameters_.absolute_tolerance_.data(); - hoststruct.absolute_tolerance_size_ = this->parameters_.absolute_tolerance_.size(); // Copy the data from host struct to device struct this->devstruct_ = micm::cuda::CopyConstData(hoststruct); }; @@ -115,11 +114,16 @@ 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) + 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_); + y_old.AsDeviceParam(), + y_new.AsDeviceParam(), + errors.AsDeviceParam(), + state.absolute_tolerance_param_, + state.GetRelativeTolerance(), + 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..953859bde 100644 --- a/include/micm/cuda/solver/cuda_state.hpp +++ b/include/micm/cuda/solver/cuda_state.hpp @@ -16,13 +16,53 @@ namespace micm public: CudaState(const CudaState&) = delete; CudaState& operator=(const CudaState&) = delete; - CudaState(CudaState&&) = default; - CudaState& operator=(CudaState&&) = default; + + CudaMatrixParam absolute_tolerance_param_; + + ~CudaState() + { + CHECK_CUDA_ERROR(micm::cuda::FreeVector(absolute_tolerance_param_), "cudaFree"); + } /// @brief Constructor which takes the state dimension information as input /// @param parameters State dimension information CudaState(const StateParameters& parameters) - : State(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) 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/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/jit/solver/jit_rosenbrock.hpp b/include/micm/jit/solver/jit_rosenbrock.hpp index 36d7e3ad0..cc6ba45ee 100644 --- a/include/micm/jit/solver/jit_rosenbrock.hpp +++ b/include/micm/jit/solver/jit_rosenbrock.hpp @@ -73,12 +73,13 @@ 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) { this->GenerateAlphaMinusJacobian(jacobian); } diff --git a/include/micm/solver/backward_euler.hpp b/include/micm/solver/backward_euler.hpp index 808ce8d05..6e135cfeb 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) requires(!VectorizableDense); + const BackwardEulerSolverParameters& parameters, + const DenseMatrixPolicy& residual, + 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& state) requires(VectorizableDense); + const BackwardEulerSolverParameters& parameters, + const DenseMatrixPolicy& residual, + 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 25fe68353..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); + converged = IsConverged(parameters_, forcing, Yn1, state.GetAbsoluteTolerances(), state.GetRelativeTolerance()); } 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) requires(!VectorizableDense) + const DenseMatrixPolicy& Yn1, const std::vector& absolute_tolerance, double relative_tolerance) requires(!VectorizableDense) { double small = parameters.small_; - double rel_tol = parameters.relative_tolerance_; - auto& abs_tol = parameters.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,27 +206,26 @@ namespace micm inline bool BackwardEuler::IsConverged( const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& state) requires(VectorizableDense) + const DenseMatrixPolicy& Yn1, const std::vector& absolute_tolerance, double relative_tolerance) requires(VectorizableDense) { double small = parameters.small_; - double rel_tol = parameters.relative_tolerance_; - auto& abs_tol = parameters.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(); 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) { 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 @@ -239,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/include/micm/solver/rosenbrock.hpp b/include/micm/solver/rosenbrock.hpp index a2baa59cd..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)), @@ -116,10 +117,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 @@ -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 bb2226819..78818c752 100644 --- a/include/micm/solver/rosenbrock.inl +++ b/include/micm/solver/rosenbrock.inl @@ -116,7 +116,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( @@ -242,7 +242,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 @@ -252,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 = parameters_.absolute_tolerance_.size(); + const std::size_t n_vars = atol.size(); double ymax = 0; double errors_over_scale = 0; @@ -263,7 +266,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] / (atol[i % n_vars] + rtol * ymax); error += errors_over_scale * errors_over_scale; } @@ -277,7 +280,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 @@ -287,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 = parameters_.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(); @@ -300,8 +306,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 / (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; @@ -319,8 +325,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]))); + (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/rosenbrock_solver_parameters.hpp b/include/micm/solver/rosenbrock_solver_parameters.hpp index 1c5e6037c..456daa21b 100644 --- a/include/micm/solver/rosenbrock_solver_parameters.hpp +++ b/include/micm/solver/rosenbrock_solver_parameters.hpp @@ -60,9 +60,6 @@ namespace micm // 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; @@ -130,12 +127,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 f8fc816d0..72b975073 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, SolverParametersType& params) + { + solver_.parameters_ = params; + return solver_.Solve(time_step, state); + } + /// @brief Returns the number of grid cells /// @return std::size_t GetNumberOfGridCells() const @@ -81,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.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 10677664a..4f5bc76d4 100644 --- a/include/micm/solver/solver_builder.inl +++ b/include/micm/solver/solver_builder.inl @@ -316,7 +316,7 @@ namespace micm } } } - } + } template< class SolverParametersPolicy, @@ -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,9 +397,12 @@ namespace micm .variable_names_ = variable_names, .custom_rate_parameter_labels_ = labels, .nonzero_jacobian_elements_ = nonzero_elements }; + + + 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/include/micm/solver/state.hpp b/include/micm/solver/state.hpp index ff8f23817..ac7914819 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_{ 1e-06 }; + std::vector absolute_tolerance_ {}; }; template @@ -59,6 +61,11 @@ namespace micm std::size_t number_of_grid_cells_; std::unique_ptr temporary_variables_; + 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(); @@ -84,6 +91,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 @@ -106,6 +115,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; } @@ -125,7 +136,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_)) { } @@ -156,6 +169,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); @@ -180,6 +195,18 @@ namespace micm void SetCustomRateParameter(const std::string& label, double value); void SetCustomRateParameter(const std::string& label, const std::vector& values); + /// @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 20a3d725b..5d76ee50c 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_(parameters.relative_tolerance_), + absolute_tolerance_(parameters.absolute_tolerance_) { std::size_t index = 0; for (auto& name : variable_names_) @@ -185,6 +187,30 @@ 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::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/solver/rosenbrock.cu b/src/solver/rosenbrock.cu index 6a6ab6961..1fa0f453e 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 @@ -137,7 +116,8 @@ namespace micm __global__ void NormalizedErrorKernel( const CudaMatrixParam y_old_param, const CudaMatrixParam y_new_param, - const RosenbrockSolverParameters ros_param, + const CudaMatrixParam absolute_tolerance_param, + const double relative_tolerance, CudaRosenbrockSolverParam devstruct, const size_t n, bool is_first_call) @@ -146,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. @@ -245,7 +225,8 @@ namespace micm __global__ void ScaledErrorKernel( const CudaMatrixParam y_old_param, const CudaMatrixParam y_new_param, - const RosenbrockSolverParameters ros_param, + const CudaMatrixParam absolute_tolerance_param, + const double relative_tolerance, CudaRosenbrockSolverParam devstruct) { // Local device variables @@ -253,8 +234,8 @@ 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 rtol = ros_param.relative_tolerance_; + const double* const atol = absolute_tolerance_param.d_data_; + 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_; @@ -288,7 +269,8 @@ namespace micm const CudaMatrixParam& y_old_param, const CudaMatrixParam& y_new_param, const CudaMatrixParam& errors_param, - const RosenbrockSolverParameters& ros_param, + const CudaMatrixParam& absolute_tolerance_param, + const double relative_tolerance, CudaRosenbrockSolverParam devstruct) { double normalized_error; @@ -318,7 +300,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, 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( @@ -341,7 +323,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, relative_tolerance, devstruct, number_of_elements, is_first_call); is_first_call = false; while (number_of_blocks > 1) { @@ -355,7 +337,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, relative_tolerance, devstruct, number_of_blocks, is_first_call); break; } NormalizedErrorKernel<<< @@ -363,7 +345,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, relative_tolerance, devstruct, number_of_blocks, is_first_call); number_of_blocks = new_number_of_blocks; } 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 aeea71fdb..07dbb2014 100644 --- a/test/integration/analytical_policy.hpp +++ b/test/integration/analytical_policy.hpp @@ -1331,6 +1331,8 @@ void test_analytical_robertson( double air_density = 1e6; auto state = solver.GetState(); + state.SetRelativeTolerance(1e-10); + state.SetAbsoluteTolerances(std::vector(3, state.GetRelativeTolerance() * 1e-2)); double k1 = 0.04; double k2 = 3e7; @@ -1545,6 +1547,9 @@ void test_analytical_oregonator( auto state = solver.GetState(); + 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); state.SetCustomRateParameter("r3", 8e3 * 0.06); @@ -1720,6 +1725,8 @@ void test_analytical_hires( }; auto state = solver.GetState(); + state.SetRelativeTolerance(1e-6); + state.SetAbsoluteTolerances(std::vector(8, state.GetRelativeTolerance() * 1e-2)); state.SetCustomRateParameter("r1", 1.71); state.SetCustomRateParameter("r2", 8.75); @@ -1878,6 +1885,14 @@ void test_analytical_e5( auto state = solver.GetState(); + 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); 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 f5fa72d5d..d4fb6e266 100644 --- a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp +++ b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp @@ -142,31 +142,6 @@ TEST(AnalyticalExamplesCudaRosenbrock, BranchedSuperStiffButAnalytical) test_analytical_stiff_branched(six_da, 2e-3, copy_to_device, copy_to_host); } -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); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_robertson(solver, 2e-1, copy_to_device, copy_to_host); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_robertson(solver, 2e-1, copy_to_device, copy_to_host); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_robertson(solver, 2e-1, copy_to_device, copy_to_host); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_robertson(solver, 2e-1, copy_to_device, copy_to_host); - - solver = rosenbrock_solver(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); @@ -176,84 +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 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); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_e5(solver, 1e-3, copy_to_device, copy_to_host); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_e5(solver, 1e-3, copy_to_device, copy_to_host); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_e5(solver, 1e-3, copy_to_device, copy_to_host); - - solver = rosenbrock_solver(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 = rosenbrock_solver(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 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); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_oregonator(solver, 2e-3, copy_to_device, copy_to_host); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_oregonator(solver, 2e-3, copy_to_device, copy_to_host); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_oregonator(solver, 2e-3, copy_to_device, copy_to_host); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_oregonator(solver, 2e-3, copy_to_device, copy_to_host); - - solver = rosenbrock_solver(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 rosenbrock_solver = [](auto params) - { - params.relative_tolerance_ = 1e-6; - params.absolute_tolerance_ = std::vector(8, params.relative_tolerance_ * 1e-2); - return builderType1Cell(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_hires(solver, 1e-6, copy_to_device, copy_to_host); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_hires(solver, 1e-7, copy_to_device, copy_to_host); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_hires(solver, 1e-7, copy_to_device, copy_to_host); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_hires(solver, 1e-6, copy_to_device, copy_to_host); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_hires(solver, 1e-6, copy_to_device, copy_to_host); -} \ No newline at end of file + 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/jit/test_jit_analytical_rosenbrock.cpp b/test/integration/jit/test_jit_analytical_rosenbrock.cpp index 7307e68e1..a099ee351 100644 --- a/test/integration/jit/test_jit_analytical_rosenbrock.cpp +++ b/test/integration/jit/test_jit_analytical_rosenbrock.cpp @@ -144,31 +144,6 @@ TEST(AnalyticalExamplesJitRosenbrock, BranchedSuperStiffButAnalytical) test_analytical_stiff_branched, StateType>(BuilderType(param_six_da), 2e-3); } -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); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_robertson, StateType<1>>(solver, 1e-1); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_robertson, StateType<1>>(solver, 1e-1); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_robertson, StateType<1>>(solver, 1e-1); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_robertson, StateType<1>>(solver, 1e-1); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_robertson, StateType<1>>(solver, 1e-1); -} - TEST(AnalyticalExamplesJitRosenbrock, SurfaceRxn) { test_analytical_surface_rxn, StateType<1>>(two, 1e-4); @@ -178,84 +153,38 @@ TEST(AnalyticalExamplesJitRosenbrock, SurfaceRxn) test_analytical_surface_rxn, StateType<1>>(six_da, 1e-5); } -TEST(AnalyticalExamplesJitRosenbrock, E5) +TEST(AnalyticalExamplesJitRosenbrock, Robertson) { - 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); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_e5, StateType<1>>(solver, 1e-3); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_e5, StateType<1>>(solver, 1e-3); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_e5, StateType<1>>(solver, 1e-3); - - solver = rosenbrock_solver(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 = rosenbrock_solver(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 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); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_oregonator, StateType<1>>(solver, 4e-4); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_oregonator, StateType<1>>(solver, 4e-4); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_oregonator, StateType<1>>(solver, 4e-4); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_oregonator, StateType<1>>(solver, 4e-4); - - solver = rosenbrock_solver(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 rosenbrock_solver = [](auto params) - { - params.relative_tolerance_ = 1e-6; - params.absolute_tolerance_ = std::vector(8, params.relative_tolerance_ * 1e-2); - return BuilderType<1>(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_hires, StateType<1>>(solver, 1e-6); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_hires, StateType<1>>(solver, 1e-7); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_hires, StateType<1>>(solver, 1e-7); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_hires, StateType<1>>(solver, 1e-6); - - solver = rosenbrock_solver(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/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); diff --git a/test/integration/terminator.hpp b/test/integration/terminator.hpp index e78e77ba9..7ffc1329c 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.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_analytical_rosenbrock.cpp b/test/integration/test_analytical_rosenbrock.cpp index 0c9799cd5..2bd3ca1fa 100644 --- a/test/integration/test_analytical_rosenbrock.cpp +++ b/test/integration/test_analytical_rosenbrock.cpp @@ -194,118 +194,47 @@ TEST(AnalyticalExamples, BranchedSuperStiffButAnalytical) test_analytical_stiff_branched, VectorStateType<4>>(rosenbrock_vector_4, 2e-3); } -TEST(AnalyticalExamples, Robertson) +TEST(AnalyticalExamples, SurfaceRxn) { - auto rosenbrock_solver = [](auto params) - { - params.relative_tolerance_ = 1e-10; - params.absolute_tolerance_ = std::vector(3, params.relative_tolerance_ * 1e-2); - return BuilderType(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_robertson(solver, 1e-1); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_robertson(solver, 1e-1); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_robertson(solver, 1e-1); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_robertson(solver, 1e-1); - - solver = rosenbrock_solver(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 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); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_e5(solver, 1e-3); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_e5(solver, 1e-3); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_e5(solver, 1e-3); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_e5(solver, 1e-3); - - solver = rosenbrock_solver(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 rosenbrock_solver = [](auto params) - { - // anything below 1e-6 is too strict for the Oregonator - params.relative_tolerance_ = 1e-8; - params.absolute_tolerance_ = std::vector(5, params.relative_tolerance_ * 1e-6); - return micm::CpuSolverBuilder(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_oregonator(solver, 4e-6); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_oregonator(solver, 4e-6); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_oregonator(solver, 4e-6); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_oregonator(solver, 4e-6); - - solver = rosenbrock_solver(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 rosenbrock_solver = [](auto params) - { - params.relative_tolerance_ = 1e-6; - params.absolute_tolerance_ = std::vector(8, params.relative_tolerance_ * 1e-2); - return micm::CpuSolverBuilder(params); - }; - - auto solver = rosenbrock_solver(micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters()); - test_analytical_hires(solver, 1e-6); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - test_analytical_hires(solver, 1e-7); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageRosenbrockParameters()); - test_analytical_hires(solver, 1e-7); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_hires(solver, 1e-6); - - solver = rosenbrock_solver(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); } diff --git a/test/integration/test_chapman_integration.cpp b/test/integration/test_chapman_integration.cpp index 0793e2d12..a9d265dca 100644 --- a/test/integration/test_chapman_integration.cpp +++ b/test/integration/test_chapman_integration.cpp @@ -32,28 +32,29 @@ 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_) .SetReactions(solver_params.processes_) - .SetIgnoreUnusedSpecies(true) - .Build(); + .SetIgnoreUnusedSpecies(true) + .Build(); micm::State state = solver.GetState(); + 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) { - 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(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/integration/test_terminator.cpp b/test/integration/test_terminator.cpp index 0384301de..c9eee218e 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); @@ -59,7 +58,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/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 diff --git a/test/unit/cuda/solver/test_cuda_rosenbrock.cpp b/test/unit/cuda/solver/test_cuda_rosenbrock.cpp index 73586f7e1..e4b838567 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(); + auto& atol = state.GetAbsoluteTolerances() ; + double rtol = state.GetRelativeTolerance(); - 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,9 @@ void testNormalizedErrorConst() y_new.CopyToDevice(); errors.CopyToDevice(); - double error = gpu_solver.solver_.NormalizedError(y_old, y_new, errors); + state.SyncInputsToDevice(); + + 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,15 +167,15 @@ 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(); + 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); + state.SyncInputsToDevice(); + double expected_error = 0.0; for (size_t i = 0; i < number_of_grid_cells; ++i) { @@ -195,7 +196,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/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 diff --git a/test/unit/solver/test_backward_euler.cpp b/test/unit/solver/test_backward_euler.cpp index df10aa3f3..b507b6045 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.SetAbsoluteTolerances({ 1e-6, 1e-6, 1e-6 }); state.variables_[0] = { 1.0, 0.0, 0.0 }; state.conditions_[0].temperature_ = 272.5; @@ -68,25 +68,25 @@ void CheckIsConverged() micm::BackwardEulerSolverParameters parameters; DenseMatrixPolicy residual{ 4, 3, 0.0 }; - DenseMatrixPolicy state{ 4, 3, 0.0 }; + DenseMatrixPolicy Yn1{ 4, 3, 0.0 }; parameters.small_ = 1e-6; - parameters.relative_tolerance_ = 1e-3; - parameters.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)); + ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance)); residual[0][1] = 1e-5; - ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, state)); + ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance)); parameters.small_ = 1e-4; - ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state)); + ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, Yn1, absolute_tolerance, relative_tolerance)); residual[3][2] = 1e-3; - ASSERT_FALSE(BackwardEuler::IsConverged(parameters, residual, state)); - state[3][2] = 10.0; - ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state)); + 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)); - parameters.absolute_tolerance_[2] = 1.0; - ASSERT_TRUE(BackwardEuler::IsConverged(parameters, residual, state)); + 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) diff --git a/test/unit/solver/test_rosenbrock.cpp b/test/unit/solver/test_rosenbrock.cpp index 357f2b683..812093efb 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(); + 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); 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,11 @@ 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(); + 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); } } @@ -134,4 +136,4 @@ TEST(RosenbrockSolver, VectorNormalizedError) testNormalizedErrorDiff(VectorBuilder<4>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()), 3); testNormalizedErrorDiff(VectorBuilder<8>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()), 5); testNormalizedErrorDiff(VectorBuilder<10>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()), 3); -} \ No newline at end of file +} diff --git a/test/unit/solver/test_solver_builder.cpp b/test/unit/solver/test_solver_builder.cpp index 28c53bb5f..04b55c23c 100644 --- a/test/unit/solver/test_solver_builder.cpp +++ b/test/unit/solver/test_solver_builder.cpp @@ -106,25 +106,60 @@ TEST(SolverBuilder, CanBuildRosenbrock) .Build(); } -TEST(SolverBuilder, MismatchedToleranceSizeIsCaught) +TEST(SolverBuilder, CanBuildBackwardEulerOverloadedSolverMethod) { - 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); + 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); +} - EXPECT_ANY_THROW(builder.SetSystem(the_system).SetReactions(reactions).SetNumberOfGridCells(1).Build();); +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