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..498f06f63 100644 --- a/include/micm/configure/solver_config.hpp +++ b/include/micm/configure/solver_config.hpp @@ -102,18 +102,28 @@ namespace micm System system_; std::vector processes_; RosenbrockSolverParameters parameters_; + double relative_tolerance_; SolverParameters(const System& system, std::vector&& processes, const RosenbrockSolverParameters&& parameters) : system_(system), processes_(processes), - parameters_(parameters) + parameters_(parameters), + relative_tolerance_(0) { } SolverParameters(System&& system, std::vector&& processes, RosenbrockSolverParameters&& parameters) : system_(system), processes_(processes), - parameters_(parameters) + parameters_(parameters), + relative_tolerance_(0) + { + } + SolverParameters(System&& system, std::vector&& processes, RosenbrockSolverParameters&& parameters, double relative_tolerance) + : system_(system), + processes_(processes), + parameters_(parameters), + relative_tolerance_(relative_tolerance) { } }; @@ -138,6 +148,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 +284,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 +437,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 +970,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_), 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 836533d7c..451722d5a 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; }; @@ -66,24 +65,22 @@ namespace micm /// @param rates Rates calculator /// @param jacobian Jacobian matrix CudaRosenbrockSolver( - 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,12 +112,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) - const - requires(CudaMatrix && VectorizableDense) + 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.relative_tolerance_, + 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 a20ad1a29..8ffb8da41 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) + { + auto& atol = this->absolute_tolerance_; + + 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() diff --git a/include/micm/cuda/util/cuda_matrix.cuh b/include/micm/cuda/util/cuda_matrix.cuh index 713f5844c..9e58328d9 100644 --- a/include/micm/cuda/util/cuda_matrix.cuh +++ b/include/micm/cuda/util/cuda_matrix.cuh @@ -41,7 +41,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 79868bfa1..685ec1247 100644 --- a/include/micm/cuda/util/cuda_param.hpp +++ b/include/micm/cuda/util/cuda_param.hpp @@ -90,8 +90,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..2d85b8974 100644 --- a/include/micm/jit/solver/jit_rosenbrock.hpp +++ b/include/micm/jit/solver/jit_rosenbrock.hpp @@ -65,20 +65,18 @@ namespace micm } /// @brief Builds a Rosenbrock solver for the given system and solver parameters - /// @param parameters Solver parameters /// @param linear_solver Linear solver /// @param rates Rates calculator /// @param jacobian Jacobian matrix JitRosenbrockSolver( - 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 30424a6a3..288449137 100644 --- a/include/micm/solver/backward_euler.hpp +++ b/include/micm/solver/backward_euler.hpp @@ -32,7 +32,6 @@ namespace micm class BackwardEuler { public: - BackwardEulerSolverParameters parameters_; LinearSolverPolicy linear_solver_; RatesPolicy rates_; std::vector jacobian_diagonal_elements_; @@ -41,17 +40,15 @@ namespace micm using ParametersType = BackwardEulerSolverParameters; /// @brief Default constructor - /// @param parameters Solver parameters /// @param linear_solver Linear solver /// @param rates Rates calculator /// @param jacobian Jacobian matrix BackwardEuler( - const BackwardEulerSolverParameters& parameters, LinearSolverPolicy&& linear_solver, RatesPolicy&& rates, - auto& jacobian) - : parameters_(parameters), - linear_solver_(std::move(linear_solver)), + auto& jacobian, + const size_t number_of_species) + : linear_solver_(std::move(linear_solver)), rates_(std::move(rates)), jacobian_diagonal_elements_(jacobian.DiagonalIndices(0)) { @@ -68,26 +65,23 @@ namespace micm /// @param time_step Time [s] to advance the state by /// @param state The state to advance /// @return result of the solver (success or failure, and statistics) - SolverResult Solve(double time_step, auto& state) const; + SolverResult Solve(double time_step, auto& state, const BackwardEulerSolverParameters& parameters) const; /// @brief Determines whether the residual is small enough to stop the /// internal solver iteration - /// @param parameters Solver parameters /// @param residual The residual to check /// @param state The current state being solved for /// @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 b810a7c8d..1e7ac42b4 100644 --- a/include/micm/solver/backward_euler.inl +++ b/include/micm/solver/backward_euler.inl @@ -43,7 +43,7 @@ inline std::error_code make_error_code(MicmBackwardEulerErrc e) namespace micm { template - inline SolverResult BackwardEuler::Solve(double time_step, auto& state) const + inline SolverResult BackwardEuler::Solve(double time_step, auto& state, const BackwardEulerSolverParameters& parameters) const { // A fully implicit euler implementation is given by the following equation: // y_{n+1} = y_n + H * f(t_{n+1}, y_{n+1}) @@ -59,10 +59,10 @@ namespace micm SolverResult result; - std::size_t max_iter = parameters_.max_number_of_steps_; - const auto time_step_reductions = parameters_.time_step_reductions_; + std::size_t max_iter = parameters.max_number_of_steps_; + const auto time_step_reductions = parameters.time_step_reductions_; - double H = parameters_.h_start_ == 0.0 ? time_step : parameters_.h_start_; + double H = parameters.h_start_ == 0.0 ? time_step : parameters.h_start_; double t = 0.0; std::size_t n_successful_integrations = 0; std::size_t n_convergence_failures = 0; @@ -148,7 +148,7 @@ namespace micm continue; // check for convergence - converged = IsConverged(parameters_, forcing, Yn1); + converged = IsConverged(parameters, forcing, Yn1, state.absolute_tolerance_, state.relative_tolerance_); } while (!converged && iterations < max_iter); if (!converged) @@ -195,24 +195,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; } @@ -222,28 +221,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(); constexpr std::size_t L = DenseMatrixPolicy::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 @@ -256,7 +253,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 71ed89a94..9b8bb5398 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 }; double h_start_{ 0.0 }; size_t max_number_of_steps_{ 11 }; diff --git a/include/micm/solver/lu_decomposition_mozart.hpp b/include/micm/solver/lu_decomposition_mozart.hpp index 39aad214c..47132cfec 100644 --- a/include/micm/solver/lu_decomposition_mozart.hpp +++ b/include/micm/solver/lu_decomposition_mozart.hpp @@ -93,6 +93,12 @@ namespace micm LuDecompositionMozart(LuDecompositionMozart&& other) = default; LuDecompositionMozart& operator=(LuDecompositionMozart&&) = default; + /// @brief Construct an LU decomposition algorithm for a given sparse matrix + /// @param matrix Sparse matrix + template + requires(SparseMatrixConcept) + LuDecompositionMozart(const SparseMatrixPolicy& matrix); + /// @brief Construct an LU decomposition algorithm for a given sparse matrix /// @param matrix Sparse matrix template diff --git a/include/micm/solver/rosenbrock.hpp b/include/micm/solver/rosenbrock.hpp index 72dba2f3d..27a6e0e36 100644 --- a/include/micm/solver/rosenbrock.hpp +++ b/include/micm/solver/rosenbrock.hpp @@ -51,7 +51,6 @@ namespace micm class AbstractRosenbrockSolver { public: - RosenbrockSolverParameters parameters_; LinearSolverPolicy linear_solver_; RatesPolicy rates_; std::vector jacobian_diagonal_elements_; @@ -62,19 +61,17 @@ namespace micm using ParametersType = RosenbrockSolverParameters; /// @brief Default constructor - /// @param parameters Solver parameters /// @param linear_solver Linear solver /// @param rates Rates calculator /// @param jacobian Jacobian matrix /// /// Note: This constructor is not intended to be used directly. Instead, use the SolverBuilder to create a solver AbstractRosenbrockSolver( - const RosenbrockSolverParameters& parameters, LinearSolverPolicy&& linear_solver, RatesPolicy&& rates, - auto& jacobian) - : parameters_(parameters), - linear_solver_(std::move(linear_solver)), + auto& jacobian, + const size_t number_of_species) + : linear_solver_(std::move(linear_solver)), rates_(std::move(rates)), jacobian_diagonal_elements_(jacobian.DiagonalIndices(0)) { @@ -90,7 +87,7 @@ namespace micm /// @brief Advances the given step over the specified time step /// @param time_step Time [s] to advance the state by /// @return A struct containing results and a status code - SolverResult Solve(double time_step, auto& state) const noexcept; + SolverResult Solve(double time_step, auto& state, const RosenbrockSolverParameters& parameters) const noexcept; /// @brief compute [alpha * I - dforce_dy] /// @param jacobian Jacobian matrix (dforce_dy) @@ -115,11 +112,11 @@ namespace micm /// @param errors The computed errors /// @return template - double NormalizedError(const DenseMatrixPolicy& y, const DenseMatrixPolicy& y_new, const DenseMatrixPolicy& errors) const - requires(!VectorizableDense); + 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 - requires(VectorizableDense); + double NormalizedError(const DenseMatrixPolicy& y, const DenseMatrixPolicy& y_new, const DenseMatrixPolicy& errors, auto& state) const + requires(VectorizableDense); }; // end of Abstract Rosenbrock Solver template @@ -128,22 +125,21 @@ namespace micm { public: /// @brief Default constructor - /// @param parameters Solver parameters /// @param linear_solver Linear solver /// @param rates Rates calculator /// @param jacobian Jacobian matrix /// /// Note: This constructor is not intended to be used directly. Instead, use the SolverBuilder to create a solver RosenbrockSolver( - 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 cf3af5e70..2960b1d7a 100644 --- a/include/micm/solver/rosenbrock.inl +++ b/include/micm/solver/rosenbrock.inl @@ -6,7 +6,7 @@ namespace micm template inline SolverResult AbstractRosenbrockSolver::Solve( double time_step, - auto& state) const noexcept + auto& state, const RosenbrockSolverParameters& parameters) const noexcept { MICM_PROFILE_FUNCTION(); using MatrixPolicy = decltype(state.variables_); @@ -20,16 +20,16 @@ namespace micm auto& initial_forcing = derived_class_temporary_variables->initial_forcing_; auto& K = derived_class_temporary_variables->K_; auto& Yerror = derived_class_temporary_variables->Yerror_; - const double h_max = parameters_.h_max_ == 0.0 ? time_step : std::min(time_step, parameters_.h_max_); + const double h_max = parameters.h_max_ == 0.0 ? time_step : std::min(time_step, parameters.h_max_); const double h_start = - parameters_.h_start_ == 0.0 ? std::max(parameters_.h_min_, DELTA_MIN) : std::min(h_max, parameters_.h_start_); + parameters.h_start_ == 0.0 ? std::max(parameters.h_min_, DELTA_MIN) : std::min(h_max, parameters.h_start_); SolverStats stats; double present_time = 0.0; - double H = std::min(std::max(std::abs(parameters_.h_min_), std::abs(h_start)), std::abs(h_max)); + double H = std::min(std::max(std::abs(parameters.h_min_), std::abs(h_start)), std::abs(h_max)); - if (std::abs(H) <= 10 * parameters_.round_off_) + if (std::abs(H) <= 10 * parameters.round_off_) H = DELTA_MIN; // TODO: the logic above this point should be moved to the constructor and should return an error @@ -38,15 +38,15 @@ namespace micm bool reject_last_h = false; bool reject_more_h = false; - while ((present_time - time_step + parameters_.round_off_) <= 0 && (result.state_ == SolverState::Running)) + while ((present_time - time_step + parameters.round_off_) <= 0 && (result.state_ == SolverState::Running)) { - if (stats.number_of_steps_ > parameters_.max_number_of_steps_) + if (stats.number_of_steps_ > parameters.max_number_of_steps_) { result.state_ = SolverState::ConvergenceExceededMaxSteps; break; } - if (((present_time + 0.1 * H) == present_time) || (H <= parameters_.round_off_)) + if (((present_time + 0.1 * H) == present_time) || (H <= parameters.round_off_)) { result.state_ = SolverState::StepSizeTooSmall; break; @@ -72,14 +72,14 @@ namespace micm { // Compute alpha accounting for the last alpha value // This is necessary to avoid the need to re-factor the jacobian - double alpha = 1.0 / (H * parameters_.gamma_[0]) - last_alpha; + double alpha = 1.0 / (H * parameters.gamma_[0]) - last_alpha; // Form and factor the rosenbrock ode jacobian LinearFactor(alpha, Y, stats, state); last_alpha = alpha; // Compute the stages - for (uint64_t stage = 0; stage < parameters_.stages_; ++stage) + for (uint64_t stage = 0; stage < parameters.stages_; ++stage) { double stage_combinations = ((stage + 1) - 1) * ((stage + 1) - 2) / 2; if (stage == 0) @@ -88,25 +88,25 @@ namespace micm } else { - if (parameters_.new_function_evaluation_[stage]) + if (parameters.new_function_evaluation_[stage]) { Ynew.Copy(Y); for (uint64_t j = 0; j < stage; ++j) { - Ynew.Axpy(parameters_.a_[stage_combinations + j], K[j]); + Ynew.Axpy(parameters.a_[stage_combinations + j], K[j]); } K[stage].Fill(0); rates_.AddForcingTerms(state.rate_constants_, Ynew, K[stage]); stats.function_calls_ += 1; } } - if (stage + 1 < parameters_.stages_ && !parameters_.new_function_evaluation_[stage + 1]) + if (stage + 1 < parameters.stages_ && !parameters.new_function_evaluation_[stage + 1]) { K[stage + 1].Copy(K[stage]); } for (uint64_t j = 0; j < stage; ++j) { - K[stage].Axpy(parameters_.c_[stage_combinations + j] / H, K[j]); + K[stage].Axpy(parameters.c_[stage_combinations + j] / H, K[j]); } linear_solver_.Solve(K[stage], state.lower_matrix_, state.upper_matrix_); stats.solves_ += 1; @@ -114,22 +114,22 @@ namespace micm // Compute the new solution Ynew.Copy(Y); - for (uint64_t stage = 0; stage < parameters_.stages_; ++stage) - Ynew.Axpy(parameters_.m_[stage], K[stage]); + for (uint64_t stage = 0; stage < parameters.stages_; ++stage) + Ynew.Axpy(parameters.m_[stage], K[stage]); Yerror.Fill(0); - for (uint64_t stage = 0; stage < parameters_.stages_; ++stage) - Yerror.Axpy(parameters_.e_[stage], K[stage]); + for (uint64_t stage = 0; stage < parameters.stages_; ++stage) + 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( - parameters_.factor_max_, + parameters.factor_max_, std::max( - parameters_.factor_min_, - parameters_.safety_factor_ / std::pow(error, 1 / parameters_.estimator_of_local_order_))); + parameters.factor_min_, + parameters.safety_factor_ / std::pow(error, 1 / parameters.estimator_of_local_order_))); double Hnew = H * fac; stats.number_of_steps_ += 1; @@ -147,12 +147,12 @@ namespace micm result.state_ = SolverState::InfDetected; break; } - else if ((error < 1) || (H < parameters_.h_min_)) + else if ((error < 1) || (H < parameters.h_min_)) { stats.accepted_ += 1; present_time = present_time + H; Y.Copy(Ynew); - Hnew = std::max(parameters_.h_min_, std::min(Hnew, h_max)); + Hnew = std::max(parameters.h_min_, std::min(Hnew, h_max)); if (reject_last_h) { // No step size increase after a rejected step @@ -168,7 +168,7 @@ namespace micm // Reject step if (reject_more_h) { - Hnew = H * parameters_.rejection_factor_decrease_; + Hnew = H * parameters.rejection_factor_decrease_; } reject_more_h = reject_last_h; reject_last_h = true; @@ -248,8 +248,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 @@ -259,8 +259,10 @@ namespace micm auto& _y = Y.AsVector(); auto& _ynew = Ynew.AsVector(); auto& _errors = errors.AsVector(); + const auto& atol = state.absolute_tolerance_; + const auto& rtol = state.relative_tolerance_; 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; @@ -270,7 +272,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; } @@ -284,8 +286,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 @@ -295,11 +297,12 @@ namespace micm auto y_iter = Y.AsVector().begin(); auto ynew_iter = Ynew.AsVector().begin(); auto errors_iter = errors.AsVector().begin(); + const auto& atol = state.absolute_tolerance_; + auto rtol = state.relative_tolerance_; const std::size_t N = Y.NumRows() * Y.NumColumns(); constexpr std::size_t L = DenseMatrixPolicy::GroupVectorSize(); - const std::size_t n_vars = parameters_.absolute_tolerance_.size(); - const std::size_t whole_blocks = std::floor(Y.NumRows() / Y.GroupVectorSize()) * Y.GroupSize(); + const std::size_t n_vars = atol.size(); double errors_over_scale = 0; double error = 0; @@ -308,8 +311,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; @@ -327,8 +330,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 187d6e46d..38e2a5ed4 100644 --- a/include/micm/solver/solver.hpp +++ b/include/micm/solver/solver.hpp @@ -21,11 +21,11 @@ namespace micm using DenseMatrixType = typename StatePolicy::DenseMatrixPolicyType; StateParameters state_parameters_; - SolverParametersType solver_parameters_; std::vector processes_; public: SolverPolicy solver_; + SolverParametersType solver_parameters_; Solver( SolverPolicy&& solver, @@ -60,11 +60,18 @@ namespace micm SolverResult Solve(double time_step, StatePolicy& state) { - auto result = solver_.Solve(time_step, state); + auto result = solver_.Solve(time_step, state, solver_parameters_); state.variables_.Max(0.0); return result; } + // Overloaded Solve function to change parameters + SolverResult Solve(double time_step, StatePolicy& state, const SolverParametersType& params) + { + solver_parameters_ = params; + return solver_.Solve(time_step, state, params); + } + /// @brief Returns the number of grid cells /// @return std::size_t GetNumberOfGridCells() const @@ -83,7 +90,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 5b15e9dd1..5b51d65ef 100644 --- a/include/micm/solver/solver_builder.inl +++ b/include/micm/solver/solver_builder.inl @@ -337,7 +337,7 @@ namespace micm } } } - } + } template< class SolverParametersPolicy, @@ -404,8 +404,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); @@ -423,9 +422,11 @@ 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(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 469faabcf..d9a80bfdd 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< @@ -63,7 +65,10 @@ 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_; + public: /// @brief Default constructor /// Only defined to be used to create default values in types, but a default constructed state is not useable State(); @@ -89,6 +94,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 @@ -111,6 +118,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; } @@ -130,7 +139,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_)) { } @@ -161,6 +172,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); @@ -185,6 +198,14 @@ 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); + /// @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 6b6e9852b..b810f8ac9 100644 --- a/include/micm/solver/state.inl +++ b/include/micm/solver/state.inl @@ -91,7 +91,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_) @@ -250,6 +252,28 @@ namespace micm for (std::size_t i = 0; i < custom_rate_parameters_.NumRows(); ++i) custom_rate_parameters_[i][param->second] = values[i]; } + + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> + inline void State::SetRelativeTolerance(double relativeTolerance) + { + this->relative_tolerance_ = relativeTolerance; + } + + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> + inline void State::SetAbsoluteTolerances(const std::vector& absoluteTolerance) + { + this->absolute_tolerance_ = absoluteTolerance; + } template< class DenseMatrixPolicy, diff --git a/src/solver/rosenbrock.cu b/src/solver/rosenbrock.cu index e2a980e82..036def798 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 97fa1624a..e49c48208 100644 --- a/src/util/cuda_matrix.cu +++ b/src/util/cuda_matrix.cu @@ -76,7 +76,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_, @@ -140,7 +140,7 @@ namespace micm template cudaError_t MallocVector(CudaMatrixParam& param, std::size_t number_of_elements); template cudaError_t MatrixMax(CudaMatrixParam& param, double val); template cudaError_t MatrixMin(CudaMatrixParam& param, double val); - 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..8ac6b9c7f 100644 --- a/test/integration/analytical_policy.hpp +++ b/test/integration/analytical_policy.hpp @@ -140,15 +140,15 @@ using yields = std::pair; using SparseMatrixTest = micm::SparseMatrix; // Test the analytical solution for a simple A -k1-> B -k2-> C system -template +template void test_simple_system( const std::string& test_label, BuilderPolicy builder, double absolute_tolerances, std::function calculate_k1, std::function calculate_k2, - std::function prepare_for_solve, - std::function postpare_for_solve, + std::function prepare_for_solve, + std::function postpare_for_solve, std::unordered_map> custom_parameters = {}) { auto solver = builder.SetNumberOfGridCells(NUM_CELLS).Build(); @@ -256,15 +256,15 @@ void test_simple_system( } // Test the analytical solution for a simple stiff A1<-fast->A2 -k1-> B -k2-> C system -template +template void test_simple_stiff_system( const std::string& test_label, BuilderPolicy builder, double absolute_tolerances, std::function calculate_k1, std::function calculate_k2, - std::function prepare_for_solve, - std::function postpare_for_solve, + std::function prepare_for_solve, + std::function postpare_for_solve, std::unordered_map> custom_parameters = {}) { auto solver = builder.SetNumberOfGridCells(NUM_CELLS).Build(); @@ -374,12 +374,12 @@ void test_simple_stiff_system( // Specific analytical tests // /////////////////////////////// -template> +template void test_analytical_troe( BuilderPolicy builder, double absolute_tolerances = 1e-10, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { /* * A -> B, k1 @@ -416,7 +416,7 @@ void test_analytical_troe( auto processes = std::vector{ r1, r2 }; builder.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })).SetReactions(processes); - test_simple_system( + test_simple_system( "troe", builder, absolute_tolerances, @@ -440,12 +440,12 @@ void test_analytical_troe( postpare_for_solve); } -template> +template void test_analytical_stiff_troe( BuilderPolicy builder, double absolute_tolerances = 1e-5, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { /* * A1 -> B, k1 @@ -503,7 +503,7 @@ void test_analytical_stiff_troe( auto processes = std::vector{ r1, r2, r3, r4, r5 }; builder.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })).SetReactions(processes); - test_simple_stiff_system( + test_simple_stiff_system( "troe", builder, absolute_tolerances, @@ -527,12 +527,12 @@ void test_analytical_stiff_troe( postpare_for_solve); } -template> +template void test_analytical_photolysis( BuilderPolicy builder, double absolute_tolerances = 2e-6, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { /* * A -> B, k1 @@ -566,7 +566,7 @@ void test_analytical_photolysis( { "photoA", std::vector(NUM_CELLS, 2e-3) }, { "photoB", std::vector(NUM_CELLS, 3e-3) } }; - test_simple_system( + test_simple_system( "photolysis", builder, absolute_tolerances, @@ -585,12 +585,12 @@ void test_analytical_photolysis( custom_parameters); } -template> +template void test_analytical_stiff_photolysis( BuilderPolicy builder, double absolute_tolerances = 2e-5, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { /* * A1 -> B, k1 @@ -647,7 +647,7 @@ void test_analytical_stiff_photolysis( { "photoB", std::vector(NUM_CELLS, 3e-3) } }; - test_simple_stiff_system( + test_simple_stiff_system( "photolysis", builder, absolute_tolerances, @@ -666,12 +666,12 @@ void test_analytical_stiff_photolysis( custom_parameters); } -template> +template void test_analytical_ternary_chemical_activation( BuilderPolicy builder, double absolute_tolerances = 1e-08, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { /* * A -> B, k1 @@ -709,7 +709,7 @@ void test_analytical_ternary_chemical_activation( auto processes = std::vector{ r1, r2 }; builder.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })).SetReactions(processes); - test_simple_system( + test_simple_system( "ternary_chemical_activation", builder, absolute_tolerances, @@ -733,12 +733,12 @@ void test_analytical_ternary_chemical_activation( postpare_for_solve); } -template> +template void test_analytical_stiff_ternary_chemical_activation( BuilderPolicy builder, double absolute_tolerances = 1e-6, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { /* * A1 -> B, k1 @@ -796,7 +796,7 @@ void test_analytical_stiff_ternary_chemical_activation( auto processes = std::vector{ r1, r2, r3, r4, r5 }; builder.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })).SetReactions(processes); - test_simple_stiff_system( + test_simple_stiff_system( "ternary_chemical_activation", builder, absolute_tolerances, @@ -820,12 +820,12 @@ void test_analytical_stiff_ternary_chemical_activation( postpare_for_solve); } -template> +template void test_analytical_tunneling( BuilderPolicy builder, double absolute_tolerances = 1e-8, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { /* * A -> B, k1 @@ -856,7 +856,7 @@ void test_analytical_tunneling( auto processes = std::vector{ r1, r2 }; builder.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })).SetReactions(processes); - test_simple_system( + test_simple_system( "tunneling", builder, absolute_tolerances, @@ -876,12 +876,12 @@ void test_analytical_tunneling( postpare_for_solve); } -template> +template void test_analytical_stiff_tunneling( BuilderPolicy builder, double absolute_tolerances = 1e-6, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { /* * A1 -> B, k1 @@ -932,7 +932,7 @@ void test_analytical_stiff_tunneling( auto processes = std::vector{ r1, r2, r3, r4, r5 }; builder.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })).SetReactions(processes); - test_simple_stiff_system( + test_simple_stiff_system( "tunneling", builder, absolute_tolerances, @@ -952,12 +952,12 @@ void test_analytical_stiff_tunneling( postpare_for_solve); } -template> +template void test_analytical_arrhenius( BuilderPolicy builder, double absolute_tolerances = 1e-9, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { /* * A -> B, k1 @@ -987,7 +987,7 @@ void test_analytical_arrhenius( auto processes = std::vector{ r1, r2 }; builder.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })).SetReactions(processes); - test_simple_system( + test_simple_system( "arrhenius", builder, absolute_tolerances, @@ -1007,12 +1007,12 @@ void test_analytical_arrhenius( postpare_for_solve); } -template> +template void test_analytical_stiff_arrhenius( BuilderPolicy builder, double absolute_tolerances = 1e-6, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { /* * A1 -> B, k1 @@ -1064,7 +1064,7 @@ void test_analytical_stiff_arrhenius( auto processes = std::vector{ r1, r2, r3, r4, r5 }; builder.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })).SetReactions(processes); - test_simple_stiff_system( + test_simple_stiff_system( "arrhenius", builder, absolute_tolerances, @@ -1084,12 +1084,12 @@ void test_analytical_stiff_arrhenius( postpare_for_solve); } -template> +template void test_analytical_branched( BuilderPolicy builder, double absolute_tolerances = 1e-13, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { /* * A -> B, k1 @@ -1129,7 +1129,7 @@ void test_analytical_branched( auto processes = std::vector{ r1, r2 }; builder.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })).SetReactions(processes); - test_simple_system( + test_simple_system( "branched", builder, absolute_tolerances, @@ -1167,12 +1167,12 @@ void test_analytical_branched( postpare_for_solve); } -template> +template void test_analytical_stiff_branched( BuilderPolicy builder, double absolute_tolerances = 1e-6, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { /* * A1 -> B, k1 @@ -1238,7 +1238,7 @@ void test_analytical_stiff_branched( auto processes = std::vector{ r1, r2, r3, r4, r5 }; builder.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })).SetReactions(processes); - test_simple_stiff_system( + test_simple_stiff_system( "branched", builder, absolute_tolerances, @@ -1276,12 +1276,12 @@ void test_analytical_stiff_branched( postpare_for_solve); } -template> +template void test_analytical_robertson( BuilderPolicy builder, double relative_tolerance = 1e-8, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { /* * A -> B, k1 = 0.04 @@ -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.relative_tolerance_ * 1e-2)); double k1 = 0.04; double k2 = 3e7; @@ -1421,12 +1423,12 @@ void test_analytical_robertson( } } -template> +template void test_analytical_oregonator( BuilderPolicy builder, double absolute_tolerance = 1e-8, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { /* * This problem is described in @@ -1545,6 +1547,9 @@ void test_analytical_oregonator( auto state = solver.GetState(); + state.SetRelativeTolerance(1e-6); + state.SetAbsoluteTolerances(std::vector(5, state.relative_tolerance_ * 1e-6)); + state.SetCustomRateParameter("r1", 1.34 * 0.06); state.SetCustomRateParameter("r2", 1.6e9); state.SetCustomRateParameter("r3", 8e3 * 0.06); @@ -1596,12 +1601,12 @@ void test_analytical_oregonator( } } -template> +template void test_analytical_hires( BuilderPolicy builder, double absolute_tolerance = 1e-8, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { /* * This problem is described in @@ -1720,6 +1725,8 @@ void test_analytical_hires( }; auto state = solver.GetState(); + state.SetRelativeTolerance(1e-6); + state.SetAbsoluteTolerances(std::vector(8, state.relative_tolerance_ * 1e-2)); state.SetCustomRateParameter("r1", 1.71); state.SetCustomRateParameter("r2", 8.75); @@ -1790,12 +1797,12 @@ void test_analytical_hires( } } -template> +template void test_analytical_e5( BuilderPolicy builder, double relative_tolerance = 1e-8, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { /* * A1 -> A2 + A3, k1 = 7.89e-10 @@ -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.absolute_tolerance_; + 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/analytical_surface_rxn_policy.hpp b/test/integration/analytical_surface_rxn_policy.hpp index 5a1710ee3..3eb1ca8d3 100644 --- a/test/integration/analytical_surface_rxn_policy.hpp +++ b/test/integration/analytical_surface_rxn_policy.hpp @@ -7,12 +7,12 @@ #include -template> +template void test_analytical_surface_rxn( BuilderPolicy& builder, double tolerance = 1e-8, - std::function prepare_for_solve = [](StateType& state) {}, - std::function postpare_for_solve = [](StateType& state) {}) + std::function prepare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}, + std::function postpare_for_solve = [](typename BuilderPolicy::StatePolicyType& state) {}) { // parameters, from CAMP/test/unit_rxn_data/test_rxn_surface.F90 const double mode_GMD = 1.0e-6; // mode geometric mean diameter [m] diff --git a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp index c39a11507..28c302cc4 100644 --- a/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp +++ b/test/integration/cuda/test_cuda_analytical_rosenbrock.cpp @@ -42,224 +42,153 @@ auto copy_to_host = [](auto& state) -> void { state.SyncOutputsToHost(); }; TEST(AnalyticalExamplesCudaRosenbrock, Troe) { - test_analytical_troe(two, 1e-10, copy_to_device, copy_to_host); - test_analytical_troe(three, 1e-10, copy_to_device, copy_to_host); - test_analytical_troe(four, 1e-10, copy_to_device, copy_to_host); - test_analytical_troe(four_da, 1e-10, copy_to_device, copy_to_host); - test_analytical_troe(six_da, 1e-10, copy_to_device, copy_to_host); + test_analytical_troe(two, 1e-10, copy_to_device, copy_to_host); + test_analytical_troe(three, 1e-10, copy_to_device, copy_to_host); + test_analytical_troe(four, 1e-10, copy_to_device, copy_to_host); + test_analytical_troe(four_da, 1e-10, copy_to_device, copy_to_host); + test_analytical_troe(six_da, 1e-10, copy_to_device, copy_to_host); } TEST(AnalyticalExamplesCudaRosenbrock, TroeSuperStiffButAnalytical) { - test_analytical_stiff_troe(two, 1e-5, copy_to_device, copy_to_host); - test_analytical_stiff_troe(three, 1e-5, copy_to_device, copy_to_host); - test_analytical_stiff_troe(four, 1e-5, copy_to_device, copy_to_host); - test_analytical_stiff_troe(four_da, 1e-5, copy_to_device, copy_to_host); - test_analytical_stiff_troe(six_da, 1e-5, copy_to_device, copy_to_host); + test_analytical_stiff_troe(two, 1e-5, copy_to_device, copy_to_host); + test_analytical_stiff_troe(three, 1e-5, copy_to_device, copy_to_host); + test_analytical_stiff_troe(four, 1e-5, copy_to_device, copy_to_host); + test_analytical_stiff_troe(four_da, 1e-5, copy_to_device, copy_to_host); + test_analytical_stiff_troe(six_da, 1e-5, copy_to_device, copy_to_host); } TEST(AnalyticalExamplesCudaRosenbrock, Photolysis) { - test_analytical_photolysis(two, 2e-6, copy_to_device, copy_to_host); - test_analytical_photolysis(three, 2e-6, copy_to_device, copy_to_host); - test_analytical_photolysis(four, 2e-8, copy_to_device, copy_to_host); - test_analytical_photolysis(four_da, 2e-6, copy_to_device, copy_to_host); - test_analytical_photolysis(six_da, 2e-6, copy_to_device, copy_to_host); + test_analytical_photolysis(two, 2e-6, copy_to_device, copy_to_host); + test_analytical_photolysis(three, 2e-6, copy_to_device, copy_to_host); + test_analytical_photolysis(four, 2e-8, copy_to_device, copy_to_host); + test_analytical_photolysis(four_da, 2e-6, copy_to_device, copy_to_host); + test_analytical_photolysis(six_da, 2e-6, copy_to_device, copy_to_host); } TEST(AnalyticalExamplesCudaRosenbrock, PhotolysisSuperStiffButAnalytical) { - test_analytical_stiff_photolysis(two, 2e-5, copy_to_device, copy_to_host); - test_analytical_stiff_photolysis(three, 2e-5, copy_to_device, copy_to_host); - test_analytical_stiff_photolysis(four, 2e-5, copy_to_device, copy_to_host); - test_analytical_stiff_photolysis(four_da, 2e-5, copy_to_device, copy_to_host); - test_analytical_stiff_photolysis(six_da, 2e-5, copy_to_device, copy_to_host); + test_analytical_stiff_photolysis(two, 2e-5, copy_to_device, copy_to_host); + test_analytical_stiff_photolysis(three, 2e-5, copy_to_device, copy_to_host); + test_analytical_stiff_photolysis(four, 2e-5, copy_to_device, copy_to_host); + test_analytical_stiff_photolysis(four_da, 2e-5, copy_to_device, copy_to_host); + test_analytical_stiff_photolysis(six_da, 2e-5, copy_to_device, copy_to_host); } TEST(AnalyticalExamplesCudaRosenbrock, TernaryChemicalActivation) { - test_analytical_ternary_chemical_activation(two, 1e-8, copy_to_device, copy_to_host); - test_analytical_ternary_chemical_activation(three, 1e-8, copy_to_device, copy_to_host); - test_analytical_ternary_chemical_activation(four, 1e-8, copy_to_device, copy_to_host); - test_analytical_ternary_chemical_activation(four_da, 1e-8, copy_to_device, copy_to_host); - test_analytical_ternary_chemical_activation(six_da, 1e-8, copy_to_device, copy_to_host); + test_analytical_ternary_chemical_activation(two, 1e-8, copy_to_device, copy_to_host); + test_analytical_ternary_chemical_activation(three, 1e-8, copy_to_device, copy_to_host); + test_analytical_ternary_chemical_activation(four, 1e-8, copy_to_device, copy_to_host); + test_analytical_ternary_chemical_activation(four_da, 1e-8, copy_to_device, copy_to_host); + test_analytical_ternary_chemical_activation(six_da, 1e-8, copy_to_device, copy_to_host); } TEST(AnalyticalExamplesCudaRosenbrock, TernaryChemicalActivationSuperStiffButAnalytical) { - test_analytical_stiff_ternary_chemical_activation(two, 2e-3, copy_to_device, copy_to_host); - test_analytical_stiff_ternary_chemical_activation(three, 2e-3, copy_to_device, copy_to_host); - test_analytical_stiff_ternary_chemical_activation(four, 2e-3, copy_to_device, copy_to_host); - test_analytical_stiff_ternary_chemical_activation(four_da, 2e-3, copy_to_device, copy_to_host); - test_analytical_stiff_ternary_chemical_activation(six_da, 2e-3, copy_to_device, copy_to_host); + test_analytical_stiff_ternary_chemical_activation(two, 2e-3, copy_to_device, copy_to_host); + test_analytical_stiff_ternary_chemical_activation(three, 2e-3, copy_to_device, copy_to_host); + test_analytical_stiff_ternary_chemical_activation(four, 2e-3, copy_to_device, copy_to_host); + test_analytical_stiff_ternary_chemical_activation(four_da, 2e-3, copy_to_device, copy_to_host); + test_analytical_stiff_ternary_chemical_activation(six_da, 2e-3, copy_to_device, copy_to_host); } TEST(AnalyticalExamplesCudaRosenbrock, Tunneling) { - test_analytical_tunneling(two, 2e-5, copy_to_device, copy_to_host); - test_analytical_tunneling(three, 1e-6, copy_to_device, copy_to_host); - test_analytical_tunneling(four, 1e-6, copy_to_device, copy_to_host); - test_analytical_tunneling(four_da, 1e-6, copy_to_device, copy_to_host); - test_analytical_tunneling(six_da, 1e-6, copy_to_device, copy_to_host); + test_analytical_tunneling(two, 2e-5, copy_to_device, copy_to_host); + test_analytical_tunneling(three, 1e-6, copy_to_device, copy_to_host); + test_analytical_tunneling(four, 1e-6, copy_to_device, copy_to_host); + test_analytical_tunneling(four_da, 1e-6, copy_to_device, copy_to_host); + test_analytical_tunneling(six_da, 1e-6, copy_to_device, copy_to_host); } TEST(AnalyticalExamplesCudaRosenbrock, TunnelingSuperStiffButAnalytical) { - test_analytical_stiff_tunneling(two, 1e-4, copy_to_device, copy_to_host); - test_analytical_stiff_tunneling(three, 1e-4, copy_to_device, copy_to_host); - test_analytical_stiff_tunneling(four, 1e-4, copy_to_device, copy_to_host); - test_analytical_stiff_tunneling(four_da, 1e-4, copy_to_device, copy_to_host); - test_analytical_stiff_tunneling(six_da, 1e-4, copy_to_device, copy_to_host); + test_analytical_stiff_tunneling(two, 1e-4, copy_to_device, copy_to_host); + test_analytical_stiff_tunneling(three, 1e-4, copy_to_device, copy_to_host); + test_analytical_stiff_tunneling(four, 1e-4, copy_to_device, copy_to_host); + test_analytical_stiff_tunneling(four_da, 1e-4, copy_to_device, copy_to_host); + test_analytical_stiff_tunneling(six_da, 1e-4, copy_to_device, copy_to_host); } TEST(AnalyticalExamplesCudaRosenbrock, Arrhenius) { - test_analytical_arrhenius(two, 4e-6, copy_to_device, copy_to_host); - test_analytical_arrhenius(three, 1e-9, copy_to_device, copy_to_host); - test_analytical_arrhenius(four, 1e-9, copy_to_device, copy_to_host); - test_analytical_arrhenius(four_da, 1e-9, copy_to_device, copy_to_host); - test_analytical_arrhenius(six_da, 1e-9, copy_to_device, copy_to_host); + test_analytical_arrhenius(two, 4e-6, copy_to_device, copy_to_host); + test_analytical_arrhenius(three, 1e-9, copy_to_device, copy_to_host); + test_analytical_arrhenius(four, 1e-9, copy_to_device, copy_to_host); + test_analytical_arrhenius(four_da, 1e-9, copy_to_device, copy_to_host); + test_analytical_arrhenius(six_da, 1e-9, copy_to_device, copy_to_host); } TEST(AnalyticalExamplesCudaRosenbrock, ArrheniusSuperStiffButAnalytical) { - test_analytical_stiff_arrhenius(two, 1e-4, copy_to_device, copy_to_host); - test_analytical_stiff_arrhenius(three, 2e-5, copy_to_device, copy_to_host); - test_analytical_stiff_arrhenius(four, 2e-5, copy_to_device, copy_to_host); - test_analytical_stiff_arrhenius(four_da, 2e-5, copy_to_device, copy_to_host); - test_analytical_stiff_arrhenius(six_da, 1e-5, copy_to_device, copy_to_host); + test_analytical_stiff_arrhenius(two, 1e-4, copy_to_device, copy_to_host); + test_analytical_stiff_arrhenius(three, 2e-5, copy_to_device, copy_to_host); + test_analytical_stiff_arrhenius(four, 2e-5, copy_to_device, copy_to_host); + test_analytical_stiff_arrhenius(four_da, 2e-5, copy_to_device, copy_to_host); + test_analytical_stiff_arrhenius(six_da, 1e-5, copy_to_device, copy_to_host); } TEST(AnalyticalExamplesCudaRosenbrock, Branched) { - test_analytical_branched(two, 1e-10, copy_to_device, copy_to_host); - test_analytical_branched(three, 1e-13, copy_to_device, copy_to_host); - test_analytical_branched(four, 1e-13, copy_to_device, copy_to_host); - test_analytical_branched(four_da, 1e-13, copy_to_device, copy_to_host); - test_analytical_branched(six_da, 1e-13, copy_to_device, copy_to_host); + test_analytical_branched(two, 1e-10, copy_to_device, copy_to_host); + test_analytical_branched(three, 1e-13, copy_to_device, copy_to_host); + test_analytical_branched(four, 1e-13, copy_to_device, copy_to_host); + test_analytical_branched(four_da, 1e-13, copy_to_device, copy_to_host); + test_analytical_branched(six_da, 1e-13, copy_to_device, copy_to_host); } TEST(AnalyticalExamplesCudaRosenbrock, BranchedSuperStiffButAnalytical) { - test_analytical_stiff_branched(two, 2e-3, copy_to_device, copy_to_host); - test_analytical_stiff_branched(three, 2e-3, copy_to_device, copy_to_host); - test_analytical_stiff_branched(four, 2e-3, copy_to_device, copy_to_host); - test_analytical_stiff_branched(four_da, 2e-3, copy_to_device, copy_to_host); - test_analytical_stiff_branched(six_da, 2e-3, copy_to_device, copy_to_host); + test_analytical_stiff_branched(two, 2e-3, copy_to_device, copy_to_host); + test_analytical_stiff_branched(three, 2e-3, copy_to_device, copy_to_host); + test_analytical_stiff_branched(four, 2e-3, copy_to_device, copy_to_host); + test_analytical_stiff_branched(four_da, 2e-3, copy_to_device, copy_to_host); + test_analytical_stiff_branched(six_da, 2e-3, copy_to_device, copy_to_host); } -TEST(AnalyticalExamplesCudaRosenbrock, Robertson) +TEST(AnalyticalExamplesCudaRosenbrock, SurfaceRxn) { - 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_analytical_surface_rxn(two_1_cell, 1e-2, copy_to_device, copy_to_host); + test_analytical_surface_rxn(three_1_cell, 1e-5, copy_to_device, copy_to_host); + test_analytical_surface_rxn(four_1_cell, 1e-6, copy_to_device, copy_to_host); + test_analytical_surface_rxn(four_da_1_cell, 1e-5, copy_to_device, copy_to_host); + test_analytical_surface_rxn(six_da_1_cell, 1e-7, copy_to_device, copy_to_host); } -TEST(AnalyticalExamplesCudaRosenbrock, SurfaceRxn) +TEST(AnalyticalExamplesCudaRosenbrock, Robertson) { - test_analytical_surface_rxn(two_1_cell, 1e-2, copy_to_device, copy_to_host); - test_analytical_surface_rxn(three_1_cell, 1e-5, copy_to_device, copy_to_host); - test_analytical_surface_rxn(four_1_cell, 1e-6, copy_to_device, copy_to_host); - test_analytical_surface_rxn(four_da_1_cell, 1e-5, copy_to_device, copy_to_host); - test_analytical_surface_rxn(six_da_1_cell, 1e-7, 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); } 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); - }; - - 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); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_e5(solver, 1e-3, copy_to_device, copy_to_host); + 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 7729138d7..2008e0a64 100644 --- a/test/integration/jit/test_jit_analytical_rosenbrock.cpp +++ b/test/integration/jit/test_jit_analytical_rosenbrock.cpp @@ -30,234 +30,153 @@ auto param_six_da = micm::RosenbrockSolverParameters::SixStageDifferentialAlgebr TEST(AnalyticalExamplesJitRosenbrock, Troe) { - test_analytical_troe, StateType>(BuilderType(param_two)); - test_analytical_troe, StateType>(BuilderType(param_three)); - test_analytical_troe, StateType>(BuilderType(param_four)); - test_analytical_troe, StateType>(BuilderType(param_four_da)); - test_analytical_troe, StateType>(BuilderType(param_six_da)); + test_analytical_troe(BuilderType(param_two)); + test_analytical_troe(BuilderType(param_three)); + test_analytical_troe(BuilderType(param_four)); + test_analytical_troe(BuilderType(param_four_da)); + test_analytical_troe(BuilderType(param_six_da)); } TEST(AnalyticalExamplesJitRosenbrock, TroeSuperStiffButAnalytical) { - test_analytical_stiff_troe, StateType>(BuilderType(param_two)); - test_analytical_stiff_troe, StateType>(BuilderType(param_three)); - test_analytical_stiff_troe, StateType>(BuilderType(param_four)); - test_analytical_stiff_troe, StateType>(BuilderType(param_four_da)); - test_analytical_stiff_troe, StateType>(BuilderType(param_six_da)); + test_analytical_stiff_troe(BuilderType(param_two)); + test_analytical_stiff_troe(BuilderType(param_three)); + test_analytical_stiff_troe(BuilderType(param_four)); + test_analytical_stiff_troe(BuilderType(param_four_da)); + test_analytical_stiff_troe(BuilderType(param_six_da)); } TEST(AnalyticalExamplesJitRosenbrock, Photolysis) { - test_analytical_photolysis, StateType>(BuilderType(param_two)); - test_analytical_photolysis, StateType>(BuilderType(param_three)); - test_analytical_photolysis, StateType>(BuilderType(param_four)); - test_analytical_photolysis, StateType>(BuilderType(param_four_da)); - test_analytical_photolysis, StateType>(BuilderType(param_six_da)); + test_analytical_photolysis(BuilderType(param_two)); + test_analytical_photolysis(BuilderType(param_three)); + test_analytical_photolysis(BuilderType(param_four)); + test_analytical_photolysis(BuilderType(param_four_da)); + test_analytical_photolysis(BuilderType(param_six_da)); } TEST(AnalyticalExamplesJitRosenbrock, PhotolysisSuperStiffButAnalytical) { - test_analytical_stiff_photolysis, StateType>(BuilderType(param_two)); - test_analytical_stiff_photolysis, StateType>(BuilderType(param_three)); - test_analytical_stiff_photolysis, StateType>(BuilderType(param_four)); - test_analytical_stiff_photolysis, StateType>(BuilderType(param_four_da)); - test_analytical_stiff_photolysis, StateType>(BuilderType(param_six_da)); + test_analytical_stiff_photolysis(BuilderType(param_two)); + test_analytical_stiff_photolysis(BuilderType(param_three)); + test_analytical_stiff_photolysis(BuilderType(param_four)); + test_analytical_stiff_photolysis(BuilderType(param_four_da)); + test_analytical_stiff_photolysis(BuilderType(param_six_da)); } TEST(AnalyticalExamplesJitRosenbrock, TernaryChemicalActivation) { - test_analytical_ternary_chemical_activation, StateType>( - BuilderType(param_two)); - test_analytical_ternary_chemical_activation, StateType>( - BuilderType(param_three)); - test_analytical_ternary_chemical_activation, StateType>( - BuilderType(param_four)); - test_analytical_ternary_chemical_activation, StateType>( - BuilderType(param_four_da)); - test_analytical_ternary_chemical_activation, StateType>( - BuilderType(param_six_da)); + test_analytical_ternary_chemical_activation(BuilderType(param_two)); + test_analytical_ternary_chemical_activation(BuilderType(param_three)); + test_analytical_ternary_chemical_activation(BuilderType(param_four)); + test_analytical_ternary_chemical_activation(BuilderType(param_four_da)); + test_analytical_ternary_chemical_activation(BuilderType(param_six_da)); } TEST(AnalyticalExamplesJitRosenbrock, TernaryChemicalActivationSuperStiffButAnalytical) { - test_analytical_stiff_ternary_chemical_activation, StateType>( - BuilderType(param_two), 2e-3); - test_analytical_stiff_ternary_chemical_activation, StateType>( - BuilderType(param_three), 2e-3); - test_analytical_stiff_ternary_chemical_activation, StateType>( - BuilderType(param_four), 2e-3); - test_analytical_stiff_ternary_chemical_activation, StateType>( - BuilderType(param_four_da), 2e-3); - test_analytical_stiff_ternary_chemical_activation, StateType>( - BuilderType(param_six_da), 2e-3); + test_analytical_stiff_ternary_chemical_activation(BuilderType(param_two), 2e-3); + test_analytical_stiff_ternary_chemical_activation(BuilderType(param_three), 2e-3); + test_analytical_stiff_ternary_chemical_activation(BuilderType(param_four), 2e-3); + test_analytical_stiff_ternary_chemical_activation(BuilderType(param_four_da), 2e-3); + test_analytical_stiff_ternary_chemical_activation(BuilderType(param_six_da), 2e-3); } TEST(AnalyticalExamplesJitRosenbrock, Tunneling) { - test_analytical_tunneling, StateType>(BuilderType(param_two), 2e-5); - test_analytical_tunneling, StateType>(BuilderType(param_three)); - test_analytical_tunneling, StateType>(BuilderType(param_four)); - test_analytical_tunneling, StateType>(BuilderType(param_four_da)); - test_analytical_tunneling, StateType>(BuilderType(param_six_da)); + test_analytical_tunneling(BuilderType(param_two), 2e-5); + test_analytical_tunneling(BuilderType(param_three)); + test_analytical_tunneling(BuilderType(param_four)); + test_analytical_tunneling(BuilderType(param_four_da)); + test_analytical_tunneling(BuilderType(param_six_da)); } TEST(AnalyticalExamplesJitRosenbrock, TunnelingSuperStiffButAnalytical) { - test_analytical_stiff_tunneling, StateType>(BuilderType(param_two), 1e-4); - test_analytical_stiff_tunneling, StateType>(BuilderType(param_three), 1e-4); - test_analytical_stiff_tunneling, StateType>(BuilderType(param_four), 1e-4); - test_analytical_stiff_tunneling, StateType>(BuilderType(param_four_da), 1e-4); - test_analytical_stiff_tunneling, StateType>(BuilderType(param_six_da), 1e-4); + test_analytical_stiff_tunneling(BuilderType(param_two), 1e-4); + test_analytical_stiff_tunneling(BuilderType(param_three), 1e-4); + test_analytical_stiff_tunneling(BuilderType(param_four), 1e-4); + test_analytical_stiff_tunneling(BuilderType(param_four_da), 1e-4); + test_analytical_stiff_tunneling(BuilderType(param_six_da), 1e-4); } TEST(AnalyticalExamplesJitRosenbrock, Arrhenius) { - test_analytical_arrhenius, StateType>(BuilderType(param_two), 4e-6); - test_analytical_arrhenius, StateType>(BuilderType(param_three)); - test_analytical_arrhenius, StateType>(BuilderType(param_four)); - test_analytical_arrhenius, StateType>(BuilderType(param_four_da)); - test_analytical_arrhenius, StateType>(BuilderType(param_six_da)); + test_analytical_arrhenius(BuilderType(param_two), 4e-6); + test_analytical_arrhenius(BuilderType(param_three)); + test_analytical_arrhenius(BuilderType(param_four)); + test_analytical_arrhenius(BuilderType(param_four_da)); + test_analytical_arrhenius(BuilderType(param_six_da)); } TEST(AnalyticalExamplesJitRosenbrock, ArrheniusSuperStiffButAnalytical) { - test_analytical_stiff_arrhenius, StateType>(BuilderType(param_two), 1e-4); - test_analytical_stiff_arrhenius, StateType>(BuilderType(param_three), 2e-5); - test_analytical_stiff_arrhenius, StateType>(BuilderType(param_four), 2e-5); - test_analytical_stiff_arrhenius, StateType>(BuilderType(param_four_da), 2e-5); - test_analytical_stiff_arrhenius, StateType>(BuilderType(param_six_da), 1e-5); + test_analytical_stiff_arrhenius(BuilderType(param_two), 1e-4); + test_analytical_stiff_arrhenius(BuilderType(param_three), 2e-5); + test_analytical_stiff_arrhenius(BuilderType(param_four), 2e-5); + test_analytical_stiff_arrhenius(BuilderType(param_four_da), 2e-5); + test_analytical_stiff_arrhenius(BuilderType(param_six_da), 1e-5); } TEST(AnalyticalExamplesJitRosenbrock, Branched) { - test_analytical_branched, StateType>(BuilderType(param_two), 1e-10); - test_analytical_branched, StateType>(BuilderType(param_three)); - test_analytical_branched, StateType>(BuilderType(param_four)); - test_analytical_branched, StateType>(BuilderType(param_four_da)); - test_analytical_branched, StateType>(BuilderType(param_six_da)); + test_analytical_branched(BuilderType(param_two), 1e-10); + test_analytical_branched(BuilderType(param_three)); + test_analytical_branched(BuilderType(param_four)); + test_analytical_branched(BuilderType(param_four_da)); + test_analytical_branched(BuilderType(param_six_da)); } TEST(AnalyticalExamplesJitRosenbrock, BranchedSuperStiffButAnalytical) { - test_analytical_stiff_branched, StateType>(BuilderType(param_two), 2e-3); - test_analytical_stiff_branched, StateType>(BuilderType(param_three), 2e-3); - test_analytical_stiff_branched, StateType>(BuilderType(param_four), 2e-3); - test_analytical_stiff_branched, StateType>(BuilderType(param_four_da), 2e-3); - test_analytical_stiff_branched, StateType>(BuilderType(param_six_da), 2e-3); + test_analytical_stiff_branched(BuilderType(param_two), 2e-3); + test_analytical_stiff_branched(BuilderType(param_three), 2e-3); + test_analytical_stiff_branched(BuilderType(param_four), 2e-3); + test_analytical_stiff_branched(BuilderType(param_four_da), 2e-3); + test_analytical_stiff_branched(BuilderType(param_six_da), 2e-3); } -TEST(AnalyticalExamplesJitRosenbrock, Robertson) +TEST(AnalyticalExamplesJitRosenbrock, SurfaceRxn) { - 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_analytical_surface_rxn(two, 1e-4); + test_analytical_surface_rxn(three, 1e-5); + test_analytical_surface_rxn(four, 1e-5); + test_analytical_surface_rxn(four_da, 1e-5); + test_analytical_surface_rxn(six_da, 1e-5); } -TEST(AnalyticalExamplesJitRosenbrock, SurfaceRxn) +TEST(AnalyticalExamplesJitRosenbrock, Robertson) { - test_analytical_surface_rxn, StateType<1>>(two, 1e-4); - test_analytical_surface_rxn, StateType<1>>(three, 1e-5); - test_analytical_surface_rxn, StateType<1>>(four, 1e-5); - test_analytical_surface_rxn, StateType<1>>(four_da, 1e-5); - test_analytical_surface_rxn, StateType<1>>(six_da, 1e-5); + test_analytical_robertson(BuilderType<1>(param_two), 1e-1); + test_analytical_robertson(BuilderType<1>(param_three), 1e-1); + test_analytical_robertson(BuilderType<1>(param_four), 1e-1); + test_analytical_robertson(BuilderType<1>(param_four_da), 1e-1); + test_analytical_robertson(BuilderType<1>(param_six_da), 1e-1); } 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); - }; - - 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); - - solver = rosenbrock_solver(micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters()); - test_analytical_e5, StateType<1>>(solver, 1e-3); + test_analytical_e5(BuilderType<1>(param_two), 1e-3); + test_analytical_e5(BuilderType<1>(param_three), 1e-3); + test_analytical_e5(BuilderType<1>(param_four), 1e-3); + test_analytical_e5(BuilderType<1>(param_four_da), 1e-3); + test_analytical_e5(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(BuilderType<1>(param_two), 4e-4); + test_analytical_oregonator(BuilderType<1>(param_three), 4e-4); + test_analytical_oregonator(BuilderType<1>(param_four), 4e-4); + test_analytical_oregonator(BuilderType<1>(param_four_da), 4e-4); + test_analytical_oregonator(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(BuilderType<1>(param_two), 1e-6); + test_analytical_hires(BuilderType<1>(param_three), 1e-7); + test_analytical_hires(BuilderType<1>(param_four), 1e-7); + test_analytical_hires(BuilderType<1>(param_four_da), 1e-6); + test_analytical_hires(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_backward_euler.cpp b/test/integration/test_analytical_backward_euler.cpp index 9ed9653ad..0c0d9cc52 100644 --- a/test/integration/test_analytical_backward_euler.cpp +++ b/test/integration/test_analytical_backward_euler.cpp @@ -115,157 +115,129 @@ auto backward_euler_vector_mozart_in_place_4 = VectorBackwardEulerMozartInPlace< TEST(AnalyticalExamples, Troe) { test_analytical_troe(backward_euler, 1e-6); - test_analytical_troe, VectorStateType<1>>(backard_euler_vector_1, 1e-6); - test_analytical_troe, VectorStateType<2>>(backard_euler_vector_2, 1e-6); - test_analytical_troe, VectorStateType<3>>(backard_euler_vector_3, 1e-6); - test_analytical_troe, VectorStateType<4>>(backard_euler_vector_4, 1e-6); - test_analytical_troe, VectorStateTypeDoolittle<1>>( - backward_euler_vector_doolittle_1, 1e-6); - test_analytical_troe, VectorStateTypeDoolittle<2>>( - backward_euler_vector_doolittle_2, 1e-6); - test_analytical_troe, VectorStateTypeDoolittle<3>>( - backward_euler_vector_doolittle_3, 1e-6); - test_analytical_troe, VectorStateTypeDoolittle<4>>( - backward_euler_vector_doolittle_4, 1e-6); - test_analytical_troe, VectorStateTypeDoolittleCSC<1>>( - backward_euler_vector_doolittle_csc_1, 1e-6); - test_analytical_troe, VectorStateTypeDoolittleCSC<2>>( - backward_euler_vector_doolittle_csc_2, 1e-6); - test_analytical_troe, VectorStateTypeDoolittleCSC<3>>( - backward_euler_vector_doolittle_csc_3, 1e-6); - test_analytical_troe, VectorStateTypeDoolittleCSC<4>>( - backward_euler_vector_doolittle_csc_4, 1e-6); - test_analytical_troe, VectorStateTypeMozart<1>>(backward_euler_vector_mozart_1, 1e-6); - test_analytical_troe, VectorStateTypeMozart<2>>(backward_euler_vector_mozart_2, 1e-6); - test_analytical_troe, VectorStateTypeMozart<3>>(backward_euler_vector_mozart_3, 1e-6); - test_analytical_troe, VectorStateTypeMozart<4>>(backward_euler_vector_mozart_4, 1e-6); - test_analytical_troe, VectorStateTypeMozartCSC<1>>( - backward_euler_vector_mozart_csc_1, 1e-6); - test_analytical_troe, VectorStateTypeMozartCSC<2>>( - backward_euler_vector_mozart_csc_2, 1e-6); - test_analytical_troe, VectorStateTypeMozartCSC<3>>( - backward_euler_vector_mozart_csc_3, 1e-6); - test_analytical_troe, VectorStateTypeMozartCSC<4>>( - backward_euler_vector_mozart_csc_4, 1e-6); - test_analytical_troe, VectorStateTypeDoolittleInPlace<1>>( - backward_euler_vector_doolittle_in_place_1, 1e-6); - test_analytical_troe, VectorStateTypeDoolittleInPlace<2>>( - backward_euler_vector_doolittle_in_place_2, 1e-6); - test_analytical_troe, VectorStateTypeDoolittleInPlace<3>>( - backward_euler_vector_doolittle_in_place_3, 1e-6); - test_analytical_troe, VectorStateTypeDoolittleInPlace<4>>( - backward_euler_vector_doolittle_in_place_4, 1e-6); - test_analytical_troe, VectorStateTypeMozartInPlace<1>>( - backward_euler_vector_mozart_in_place_1, 1e-6); - test_analytical_troe, VectorStateTypeMozartInPlace<2>>( - backward_euler_vector_mozart_in_place_2, 1e-6); - test_analytical_troe, VectorStateTypeMozartInPlace<3>>( - backward_euler_vector_mozart_in_place_3, 1e-6); - test_analytical_troe, VectorStateTypeMozartInPlace<4>>( - backward_euler_vector_mozart_in_place_4, 1e-6); + test_analytical_troe(backard_euler_vector_1, 1e-6); + test_analytical_troe(backard_euler_vector_2, 1e-6); + test_analytical_troe(backard_euler_vector_3, 1e-6); + test_analytical_troe(backard_euler_vector_4, 1e-6); + test_analytical_troe(backward_euler_vector_doolittle_1, 1e-6); + test_analytical_troe(backward_euler_vector_doolittle_2, 1e-6); + test_analytical_troe(backward_euler_vector_doolittle_3, 1e-6); + test_analytical_troe(backward_euler_vector_doolittle_4, 1e-6); + test_analytical_troe(backward_euler_vector_doolittle_csc_1, 1e-6); + test_analytical_troe(backward_euler_vector_doolittle_csc_2, 1e-6); + test_analytical_troe(backward_euler_vector_doolittle_csc_3, 1e-6); + test_analytical_troe(backward_euler_vector_doolittle_csc_4, 1e-6); + test_analytical_troe(backward_euler_vector_mozart_1, 1e-6); + test_analytical_troe(backward_euler_vector_mozart_2, 1e-6); + test_analytical_troe(backward_euler_vector_mozart_3, 1e-6); + test_analytical_troe(backward_euler_vector_mozart_4, 1e-6); + test_analytical_troe(backward_euler_vector_mozart_csc_1, 1e-6); + test_analytical_troe(backward_euler_vector_mozart_csc_2, 1e-6); + test_analytical_troe(backward_euler_vector_mozart_csc_3, 1e-6); + test_analytical_troe(backward_euler_vector_mozart_csc_4, 1e-6); + test_analytical_troe(backward_euler_vector_mozart_in_place_1, 1e-6); + test_analytical_troe(backward_euler_vector_mozart_in_place_2, 1e-6); + test_analytical_troe(backward_euler_vector_mozart_in_place_3, 1e-6); + test_analytical_troe(backward_euler_vector_mozart_in_place_4, 1e-6); } TEST(AnalyticalExamples, TroeSuperStiffButAnalytical) { test_analytical_stiff_troe(backward_euler); - test_analytical_stiff_troe, VectorStateType<1>>(backard_euler_vector_1); - test_analytical_stiff_troe, VectorStateType<2>>(backard_euler_vector_2); - test_analytical_stiff_troe, VectorStateType<3>>(backard_euler_vector_3); - test_analytical_stiff_troe, VectorStateType<4>>(backard_euler_vector_4); + test_analytical_stiff_troe(backard_euler_vector_1); + test_analytical_stiff_troe(backard_euler_vector_2); + test_analytical_stiff_troe(backard_euler_vector_3); + test_analytical_stiff_troe(backard_euler_vector_4); } TEST(AnalyticalExamples, Photolysis) { test_analytical_photolysis(backward_euler, 1e-3); - test_analytical_photolysis, VectorStateType<1>>(backard_euler_vector_1, 1e-3); - test_analytical_photolysis, VectorStateType<2>>(backard_euler_vector_2, 1e-3); - test_analytical_photolysis, VectorStateType<3>>(backard_euler_vector_3, 1e-3); - test_analytical_photolysis, VectorStateType<4>>(backard_euler_vector_4, 1e-3); + test_analytical_photolysis(backard_euler_vector_1, 1e-3); + test_analytical_photolysis(backard_euler_vector_2, 1e-3); + test_analytical_photolysis(backard_euler_vector_3, 1e-3); + test_analytical_photolysis(backard_euler_vector_4, 1e-3); } TEST(AnalyticalExamples, PhotolysisSuperStiffButAnalytical) { test_analytical_stiff_photolysis(backward_euler, 1e-3); - test_analytical_stiff_photolysis, VectorStateType<1>>(backard_euler_vector_1, 1e-3); - test_analytical_stiff_photolysis, VectorStateType<2>>(backard_euler_vector_2, 1e-3); - test_analytical_stiff_photolysis, VectorStateType<3>>(backard_euler_vector_3, 1e-3); - test_analytical_stiff_photolysis, VectorStateType<4>>(backard_euler_vector_4, 1e-3); + test_analytical_stiff_photolysis(backard_euler_vector_1, 1e-3); + test_analytical_stiff_photolysis(backard_euler_vector_2, 1e-3); + test_analytical_stiff_photolysis(backard_euler_vector_3, 1e-3); + test_analytical_stiff_photolysis(backard_euler_vector_4, 1e-3); } TEST(AnalyticalExamples, TernaryChemicalActivation) { test_analytical_ternary_chemical_activation(backward_euler, 1e-5); - test_analytical_ternary_chemical_activation, VectorStateType<1>>(backard_euler_vector_1, 1e-5); - test_analytical_ternary_chemical_activation, VectorStateType<2>>(backard_euler_vector_2, 1e-5); - test_analytical_ternary_chemical_activation, VectorStateType<3>>(backard_euler_vector_3, 1e-5); - test_analytical_ternary_chemical_activation, VectorStateType<4>>(backard_euler_vector_4, 1e-5); + test_analytical_ternary_chemical_activation(backard_euler_vector_1, 1e-5); + test_analytical_ternary_chemical_activation(backard_euler_vector_2, 1e-5); + test_analytical_ternary_chemical_activation(backard_euler_vector_3, 1e-5); + test_analytical_ternary_chemical_activation(backard_euler_vector_4, 1e-5); } TEST(AnalyticalExamples, TernaryChemicalActivationSuperStiffButAnalytical) { test_analytical_stiff_ternary_chemical_activation(backward_euler, 1e-2); - test_analytical_stiff_ternary_chemical_activation, VectorStateType<1>>( - backard_euler_vector_1, 1e-2); - test_analytical_stiff_ternary_chemical_activation, VectorStateType<2>>( - backard_euler_vector_2, 1e-2); - test_analytical_stiff_ternary_chemical_activation, VectorStateType<3>>( - backard_euler_vector_3, 1e-2); - test_analytical_stiff_ternary_chemical_activation, VectorStateType<4>>( - backard_euler_vector_4, 1e-2); + test_analytical_stiff_ternary_chemical_activation(backard_euler_vector_1, 1e-2); + test_analytical_stiff_ternary_chemical_activation(backard_euler_vector_2, 1e-2); + test_analytical_stiff_ternary_chemical_activation(backard_euler_vector_3, 1e-2); + test_analytical_stiff_ternary_chemical_activation(backard_euler_vector_4, 1e-2); } TEST(AnalyticalExamples, Tunneling) { test_analytical_tunneling(backward_euler, 1e-3); - test_analytical_tunneling, VectorStateType<1>>(backard_euler_vector_1, 1e-3); - test_analytical_tunneling, VectorStateType<2>>(backard_euler_vector_2, 1e-3); - test_analytical_tunneling, VectorStateType<3>>(backard_euler_vector_3, 1e-3); - test_analytical_tunneling, VectorStateType<4>>(backard_euler_vector_4, 1e-3); + test_analytical_tunneling(backard_euler_vector_1, 1e-3); + test_analytical_tunneling(backard_euler_vector_2, 1e-3); + test_analytical_tunneling(backard_euler_vector_3, 1e-3); + test_analytical_tunneling(backard_euler_vector_4, 1e-3); } TEST(AnalyticalExamples, TunnelingSuperStiffButAnalytical) { test_analytical_stiff_tunneling(backward_euler, 1e-3); - test_analytical_stiff_tunneling, VectorStateType<1>>(backard_euler_vector_1, 1e-3); - test_analytical_stiff_tunneling, VectorStateType<2>>(backard_euler_vector_2, 1e-3); - test_analytical_stiff_tunneling, VectorStateType<3>>(backard_euler_vector_3, 1e-3); - test_analytical_stiff_tunneling, VectorStateType<4>>(backard_euler_vector_4, 1e-3); + test_analytical_stiff_tunneling(backard_euler_vector_1, 1e-3); + test_analytical_stiff_tunneling(backard_euler_vector_2, 1e-3); + test_analytical_stiff_tunneling(backard_euler_vector_3, 1e-3); + test_analytical_stiff_tunneling(backard_euler_vector_4, 1e-3); } TEST(AnalyticalExamples, Arrhenius) { test_analytical_arrhenius(backward_euler, 1e-3); - test_analytical_arrhenius, VectorStateType<1>>(backard_euler_vector_1, 1e-3); - test_analytical_arrhenius, VectorStateType<2>>(backard_euler_vector_2, 1e-3); - test_analytical_arrhenius, VectorStateType<3>>(backard_euler_vector_3, 1e-3); - test_analytical_arrhenius, VectorStateType<4>>(backard_euler_vector_4, 1e-3); + test_analytical_arrhenius(backard_euler_vector_1, 1e-3); + test_analytical_arrhenius(backard_euler_vector_2, 1e-3); + test_analytical_arrhenius(backard_euler_vector_3, 1e-3); + test_analytical_arrhenius(backard_euler_vector_4, 1e-3); } TEST(AnalyticalExamples, ArrheniusSuperStiffButAnalytical) { test_analytical_stiff_arrhenius(backward_euler, 1e-3); - test_analytical_stiff_arrhenius, VectorStateType<1>>(backard_euler_vector_1, 1e-3); - test_analytical_stiff_arrhenius, VectorStateType<2>>(backard_euler_vector_2, 1e-3); - test_analytical_stiff_arrhenius, VectorStateType<3>>(backard_euler_vector_3, 1e-3); - test_analytical_stiff_arrhenius, VectorStateType<4>>(backard_euler_vector_4, 1e-3); + test_analytical_stiff_arrhenius(backard_euler_vector_1, 1e-3); + test_analytical_stiff_arrhenius(backard_euler_vector_2, 1e-3); + test_analytical_stiff_arrhenius(backard_euler_vector_3, 1e-3); + test_analytical_stiff_arrhenius(backard_euler_vector_4, 1e-3); } TEST(AnalyticalExamples, Branched) { test_analytical_branched(backward_euler, 1e-5); - test_analytical_branched, VectorStateType<1>>(backard_euler_vector_1, 1e-5); - test_analytical_branched, VectorStateType<2>>(backard_euler_vector_2, 1e-5); - test_analytical_branched, VectorStateType<3>>(backard_euler_vector_3, 1e-5); - test_analytical_branched, VectorStateType<4>>(backard_euler_vector_4, 1e-5); + test_analytical_branched(backard_euler_vector_1, 1e-5); + test_analytical_branched(backard_euler_vector_2, 1e-5); + test_analytical_branched(backard_euler_vector_3, 1e-5); + test_analytical_branched(backard_euler_vector_4, 1e-5); } TEST(AnalyticalExamples, BranchedSuperStiffButAnalytical) { test_analytical_stiff_branched(backward_euler, 1e-2); - test_analytical_stiff_branched, VectorStateType<1>>(backard_euler_vector_1, 1e-2); - test_analytical_stiff_branched, VectorStateType<2>>(backard_euler_vector_2, 1e-2); - test_analytical_stiff_branched, VectorStateType<3>>(backard_euler_vector_3, 1e-2); - test_analytical_stiff_branched, VectorStateType<4>>(backard_euler_vector_4, 1e-2); + test_analytical_stiff_branched(backard_euler_vector_1, 1e-2); + test_analytical_stiff_branched(backard_euler_vector_2, 1e-2); + test_analytical_stiff_branched(backard_euler_vector_3, 1e-2); + test_analytical_stiff_branched(backard_euler_vector_4, 1e-2); } TEST(AnalyticalExamples, SurfaceRxn) @@ -276,41 +248,33 @@ TEST(AnalyticalExamples, SurfaceRxn) TEST(AnalyticalExamples, HIRES) { test_analytical_hires(backward_euler, 1e-1); - test_analytical_hires, VectorStateType<1>>(backard_euler_vector_1, 1e-1); - test_analytical_hires, VectorStateType<2>>(backard_euler_vector_2, 1e-1); - test_analytical_hires, VectorStateType<3>>(backard_euler_vector_3, 1e-1); - test_analytical_hires, VectorStateType<4>>(backard_euler_vector_4, 1e-1); - test_analytical_hires, VectorStateTypeDoolittle<1>>( - backward_euler_vector_doolittle_1, 1e-1); - test_analytical_hires, VectorStateTypeDoolittle<2>>( - backward_euler_vector_doolittle_2, 1e-1); - test_analytical_hires, VectorStateTypeDoolittle<3>>( - backward_euler_vector_doolittle_3, 1e-1); - test_analytical_hires, VectorStateTypeDoolittle<4>>( - backward_euler_vector_doolittle_4, 1e-1); - test_analytical_hires, VectorStateTypeMozart<1>>(backward_euler_vector_mozart_1, 1e-1); - test_analytical_hires, VectorStateTypeMozart<2>>(backward_euler_vector_mozart_2, 1e-1); - test_analytical_hires, VectorStateTypeMozart<3>>(backward_euler_vector_mozart_3, 1e-1); - test_analytical_hires, VectorStateTypeMozart<4>>(backward_euler_vector_mozart_4, 1e-1); + test_analytical_hires(backard_euler_vector_1, 1e-1); + test_analytical_hires(backard_euler_vector_2, 1e-1); + test_analytical_hires(backard_euler_vector_3, 1e-1); + test_analytical_hires(backard_euler_vector_4, 1e-1); + test_analytical_hires(backward_euler_vector_doolittle_1, 1e-1); + test_analytical_hires(backward_euler_vector_doolittle_2, 1e-1); + test_analytical_hires(backward_euler_vector_doolittle_3, 1e-1); + test_analytical_hires(backward_euler_vector_doolittle_4, 1e-1); + test_analytical_hires(backward_euler_vector_mozart_1, 1e-1); + test_analytical_hires(backward_euler_vector_mozart_2, 1e-1); + test_analytical_hires(backward_euler_vector_mozart_3, 1e-1); + test_analytical_hires(backward_euler_vector_mozart_4, 1e-1); } TEST(AnalyticalExamples, Oregonator) { test_analytical_oregonator(backward_euler, 1e-3); - test_analytical_oregonator, VectorStateType<1>>(backard_euler_vector_1, 1e-3); - test_analytical_oregonator, VectorStateType<2>>(backard_euler_vector_2, 1e-3); - test_analytical_oregonator, VectorStateType<3>>(backard_euler_vector_3, 1e-3); - test_analytical_oregonator, VectorStateType<4>>(backard_euler_vector_4, 1e-3); - test_analytical_oregonator, VectorStateTypeDoolittle<1>>( - backward_euler_vector_doolittle_1, 1e-3); - test_analytical_oregonator, VectorStateTypeDoolittle<2>>( - backward_euler_vector_doolittle_2, 1e-3); - test_analytical_oregonator, VectorStateTypeDoolittle<3>>( - backward_euler_vector_doolittle_3, 1e-3); - test_analytical_oregonator, VectorStateTypeDoolittle<4>>( - backward_euler_vector_doolittle_4, 1e-3); - test_analytical_oregonator, VectorStateTypeMozart<1>>(backward_euler_vector_mozart_1, 1e-3); - test_analytical_oregonator, VectorStateTypeMozart<2>>(backward_euler_vector_mozart_2, 1e-3); - test_analytical_oregonator, VectorStateTypeMozart<3>>(backward_euler_vector_mozart_3, 1e-3); - test_analytical_oregonator, VectorStateTypeMozart<4>>(backward_euler_vector_mozart_4, 1e-3); + test_analytical_oregonator(backard_euler_vector_1, 1e-3); + test_analytical_oregonator(backard_euler_vector_2, 1e-3); + test_analytical_oregonator(backard_euler_vector_3, 1e-3); + test_analytical_oregonator(backard_euler_vector_4, 1e-3); + test_analytical_oregonator(backward_euler_vector_doolittle_1, 1e-3); + test_analytical_oregonator(backward_euler_vector_doolittle_2, 1e-3); + test_analytical_oregonator(backward_euler_vector_doolittle_3, 1e-3); + test_analytical_oregonator(backward_euler_vector_doolittle_4, 1e-3); + test_analytical_oregonator(backward_euler_vector_mozart_1, 1e-3); + test_analytical_oregonator(backward_euler_vector_mozart_2, 1e-3); + test_analytical_oregonator(backward_euler_vector_mozart_3, 1e-3); + test_analytical_oregonator(backward_euler_vector_mozart_4, 1e-3); } diff --git a/test/integration/test_analytical_rosenbrock.cpp b/test/integration/test_analytical_rosenbrock.cpp index 3cb0167eb..9b4cc33bc 100644 --- a/test/integration/test_analytical_rosenbrock.cpp +++ b/test/integration/test_analytical_rosenbrock.cpp @@ -112,20 +112,20 @@ TEST(AnalyticalExamples, Troe) test_analytical_troe(rosenbrock_4stage); test_analytical_troe(rosenbrock_4stage_da); test_analytical_troe(rosenbrock_6stage_da); - test_analytical_troe, VectorStateType<1>>(rosenbrock_vector_1); - test_analytical_troe, VectorStateType<2>>(rosenbrock_vector_2); - test_analytical_troe, VectorStateType<3>>(rosenbrock_vector_3); - test_analytical_troe, VectorStateType<4>>(rosenbrock_vector_4); - test_analytical_troe(rosenbrock_standard_doolittle); - test_analytical_troe, VectorStateTypeDoolittle<1>>(rosenbrock_vector_doolittle_1); - test_analytical_troe, VectorStateTypeDoolittle<2>>(rosenbrock_vector_doolittle_2); - test_analytical_troe, VectorStateTypeDoolittle<3>>(rosenbrock_vector_doolittle_3); - test_analytical_troe, VectorStateTypeDoolittle<4>>(rosenbrock_vector_doolittle_4); - test_analytical_troe(rosenbrock_standard_mozart); - test_analytical_troe, VectorStateTypeMozart<1>>(rosenbrock_vector_mozart_1); - test_analytical_troe, VectorStateTypeMozart<2>>(rosenbrock_vector_mozart_2); - test_analytical_troe, VectorStateTypeMozart<3>>(rosenbrock_vector_mozart_3); - test_analytical_troe, VectorStateTypeMozart<4>>(rosenbrock_vector_mozart_4); + test_analytical_troe(rosenbrock_vector_1); + test_analytical_troe(rosenbrock_vector_2); + test_analytical_troe(rosenbrock_vector_3); + test_analytical_troe(rosenbrock_vector_4); + test_analytical_troe(rosenbrock_standard_doolittle); + test_analytical_troe(rosenbrock_vector_doolittle_1); + test_analytical_troe(rosenbrock_vector_doolittle_2); + test_analytical_troe(rosenbrock_vector_doolittle_3); + test_analytical_troe(rosenbrock_vector_doolittle_4); + test_analytical_troe(rosenbrock_standard_mozart); + test_analytical_troe(rosenbrock_vector_mozart_1); + test_analytical_troe(rosenbrock_vector_mozart_2); + test_analytical_troe(rosenbrock_vector_mozart_3); + test_analytical_troe(rosenbrock_vector_mozart_4); } TEST(AnalyticalExamples, TroeSuperStiffButAnalytical) @@ -135,10 +135,10 @@ TEST(AnalyticalExamples, TroeSuperStiffButAnalytical) test_analytical_stiff_troe(rosenbrock_4stage); test_analytical_stiff_troe(rosenbrock_4stage_da); test_analytical_stiff_troe(rosenbrock_6stage_da); - test_analytical_stiff_troe, VectorStateType<1>>(rosenbrock_vector_1); - test_analytical_stiff_troe, VectorStateType<2>>(rosenbrock_vector_2); - test_analytical_stiff_troe, VectorStateType<3>>(rosenbrock_vector_3); - test_analytical_stiff_troe, VectorStateType<4>>(rosenbrock_vector_4); + test_analytical_stiff_troe(rosenbrock_vector_1); + test_analytical_stiff_troe(rosenbrock_vector_2); + test_analytical_stiff_troe(rosenbrock_vector_3); + test_analytical_stiff_troe(rosenbrock_vector_4); } TEST(AnalyticalExamples, Photolysis) @@ -148,10 +148,10 @@ TEST(AnalyticalExamples, Photolysis) test_analytical_photolysis(rosenbrock_4stage); test_analytical_photolysis(rosenbrock_4stage_da); test_analytical_photolysis(rosenbrock_6stage_da); - test_analytical_photolysis, VectorStateType<1>>(rosenbrock_vector_1); - test_analytical_photolysis, VectorStateType<2>>(rosenbrock_vector_2); - test_analytical_photolysis, VectorStateType<3>>(rosenbrock_vector_3); - test_analytical_photolysis, VectorStateType<4>>(rosenbrock_vector_4); + test_analytical_photolysis(rosenbrock_vector_1); + test_analytical_photolysis(rosenbrock_vector_2); + test_analytical_photolysis(rosenbrock_vector_3); + test_analytical_photolysis(rosenbrock_vector_4); } TEST(AnalyticalExamples, PhotolysisSuperStiffButAnalytical) @@ -161,10 +161,10 @@ TEST(AnalyticalExamples, PhotolysisSuperStiffButAnalytical) test_analytical_stiff_photolysis(rosenbrock_4stage); test_analytical_stiff_photolysis(rosenbrock_4stage_da); test_analytical_stiff_photolysis(rosenbrock_6stage_da); - test_analytical_stiff_photolysis, VectorStateType<1>>(rosenbrock_vector_1); - test_analytical_stiff_photolysis, VectorStateType<2>>(rosenbrock_vector_2); - test_analytical_stiff_photolysis, VectorStateType<3>>(rosenbrock_vector_3); - test_analytical_stiff_photolysis, VectorStateType<4>>(rosenbrock_vector_4); + test_analytical_stiff_photolysis(rosenbrock_vector_1); + test_analytical_stiff_photolysis(rosenbrock_vector_2); + test_analytical_stiff_photolysis(rosenbrock_vector_3); + test_analytical_stiff_photolysis(rosenbrock_vector_4); } TEST(AnalyticalExamples, TernaryChemicalActivation) @@ -174,10 +174,10 @@ TEST(AnalyticalExamples, TernaryChemicalActivation) test_analytical_ternary_chemical_activation(rosenbrock_4stage); test_analytical_ternary_chemical_activation(rosenbrock_4stage_da); test_analytical_ternary_chemical_activation(rosenbrock_6stage_da); - test_analytical_ternary_chemical_activation, VectorStateType<1>>(rosenbrock_vector_1); - test_analytical_ternary_chemical_activation, VectorStateType<2>>(rosenbrock_vector_2); - test_analytical_ternary_chemical_activation, VectorStateType<3>>(rosenbrock_vector_3); - test_analytical_ternary_chemical_activation, VectorStateType<4>>(rosenbrock_vector_4); + test_analytical_ternary_chemical_activation(rosenbrock_vector_1); + test_analytical_ternary_chemical_activation(rosenbrock_vector_2); + test_analytical_ternary_chemical_activation(rosenbrock_vector_3); + test_analytical_ternary_chemical_activation(rosenbrock_vector_4); } TEST(AnalyticalExamples, TernaryChemicalActivationSuperStiffButAnalytical) @@ -187,10 +187,10 @@ TEST(AnalyticalExamples, TernaryChemicalActivationSuperStiffButAnalytical) test_analytical_stiff_ternary_chemical_activation(rosenbrock_4stage, 2e-3); test_analytical_stiff_ternary_chemical_activation(rosenbrock_4stage_da, 2e-3); test_analytical_stiff_ternary_chemical_activation(rosenbrock_6stage_da, 2e-3); - test_analytical_stiff_ternary_chemical_activation, VectorStateType<1>>(rosenbrock_vector_1, 2e-3); - test_analytical_stiff_ternary_chemical_activation, VectorStateType<2>>(rosenbrock_vector_2, 2e-3); - test_analytical_stiff_ternary_chemical_activation, VectorStateType<3>>(rosenbrock_vector_3, 2e-3); - test_analytical_stiff_ternary_chemical_activation, VectorStateType<4>>(rosenbrock_vector_4, 2e-3); + test_analytical_stiff_ternary_chemical_activation(rosenbrock_vector_1, 2e-3); + test_analytical_stiff_ternary_chemical_activation(rosenbrock_vector_2, 2e-3); + test_analytical_stiff_ternary_chemical_activation(rosenbrock_vector_3, 2e-3); + test_analytical_stiff_ternary_chemical_activation(rosenbrock_vector_4, 2e-3); } TEST(AnalyticalExamples, Tunneling) @@ -200,10 +200,10 @@ TEST(AnalyticalExamples, Tunneling) test_analytical_tunneling(rosenbrock_4stage); test_analytical_tunneling(rosenbrock_4stage_da); test_analytical_tunneling(rosenbrock_6stage_da); - test_analytical_tunneling, VectorStateType<1>>(rosenbrock_vector_1); - test_analytical_tunneling, VectorStateType<2>>(rosenbrock_vector_2); - test_analytical_tunneling, VectorStateType<3>>(rosenbrock_vector_3); - test_analytical_tunneling, VectorStateType<4>>(rosenbrock_vector_4); + test_analytical_tunneling(rosenbrock_vector_1); + test_analytical_tunneling(rosenbrock_vector_2); + test_analytical_tunneling(rosenbrock_vector_3); + test_analytical_tunneling(rosenbrock_vector_4); } TEST(AnalyticalExamples, TunnelingSuperStiffButAnalytical) @@ -213,10 +213,10 @@ TEST(AnalyticalExamples, TunnelingSuperStiffButAnalytical) test_analytical_stiff_tunneling(rosenbrock_4stage, 1e-4); test_analytical_stiff_tunneling(rosenbrock_4stage_da, 1e-4); test_analytical_stiff_tunneling(rosenbrock_6stage_da, 1e-4); - test_analytical_stiff_tunneling, VectorStateType<1>>(rosenbrock_vector_1, 1e-4); - test_analytical_stiff_tunneling, VectorStateType<2>>(rosenbrock_vector_2, 1e-4); - test_analytical_stiff_tunneling, VectorStateType<3>>(rosenbrock_vector_3, 1e-4); - test_analytical_stiff_tunneling, VectorStateType<4>>(rosenbrock_vector_4, 1e-4); + test_analytical_stiff_tunneling(rosenbrock_vector_1, 1e-4); + test_analytical_stiff_tunneling(rosenbrock_vector_2, 1e-4); + test_analytical_stiff_tunneling(rosenbrock_vector_3, 1e-4); + test_analytical_stiff_tunneling(rosenbrock_vector_4, 1e-4); } TEST(AnalyticalExamples, Arrhenius) @@ -226,10 +226,10 @@ TEST(AnalyticalExamples, Arrhenius) test_analytical_arrhenius(rosenbrock_4stage); test_analytical_arrhenius(rosenbrock_4stage_da); test_analytical_arrhenius(rosenbrock_6stage_da); - test_analytical_arrhenius, VectorStateType<1>>(rosenbrock_vector_1); - test_analytical_arrhenius, VectorStateType<2>>(rosenbrock_vector_2); - test_analytical_arrhenius, VectorStateType<3>>(rosenbrock_vector_3); - test_analytical_arrhenius, VectorStateType<4>>(rosenbrock_vector_4); + test_analytical_arrhenius(rosenbrock_vector_1); + test_analytical_arrhenius(rosenbrock_vector_2); + test_analytical_arrhenius(rosenbrock_vector_3); + test_analytical_arrhenius(rosenbrock_vector_4); } TEST(AnalyticalExamples, ArrheniusSuperStiffButAnalytical) @@ -239,10 +239,10 @@ TEST(AnalyticalExamples, ArrheniusSuperStiffButAnalytical) test_analytical_stiff_arrhenius(rosenbrock_4stage, 2e-5); test_analytical_stiff_arrhenius(rosenbrock_4stage_da, 2e-5); test_analytical_stiff_arrhenius(rosenbrock_6stage_da, 1e-5); - test_analytical_stiff_arrhenius, VectorStateType<1>>(rosenbrock_vector_1, 2e-5); - test_analytical_stiff_arrhenius, VectorStateType<2>>(rosenbrock_vector_2, 2e-5); - test_analytical_stiff_arrhenius, VectorStateType<3>>(rosenbrock_vector_3, 2e-5); - test_analytical_stiff_arrhenius, VectorStateType<4>>(rosenbrock_vector_4, 2e-5); + test_analytical_stiff_arrhenius(rosenbrock_vector_1, 2e-5); + test_analytical_stiff_arrhenius(rosenbrock_vector_2, 2e-5); + test_analytical_stiff_arrhenius(rosenbrock_vector_3, 2e-5); + test_analytical_stiff_arrhenius(rosenbrock_vector_4, 2e-5); } TEST(AnalyticalExamples, Branched) @@ -252,10 +252,10 @@ TEST(AnalyticalExamples, Branched) test_analytical_branched(rosenbrock_4stage); test_analytical_branched(rosenbrock_4stage_da); test_analytical_branched(rosenbrock_6stage_da); - test_analytical_branched, VectorStateType<1>>(rosenbrock_vector_1); - test_analytical_branched, VectorStateType<2>>(rosenbrock_vector_2); - test_analytical_branched, VectorStateType<3>>(rosenbrock_vector_3); - test_analytical_branched, VectorStateType<4>>(rosenbrock_vector_4); + test_analytical_branched(rosenbrock_vector_1); + test_analytical_branched(rosenbrock_vector_2); + test_analytical_branched(rosenbrock_vector_3); + test_analytical_branched(rosenbrock_vector_4); } TEST(AnalyticalExamples, BranchedSuperStiffButAnalytical) @@ -265,227 +265,89 @@ TEST(AnalyticalExamples, BranchedSuperStiffButAnalytical) test_analytical_stiff_branched(rosenbrock_4stage, 2e-3); test_analytical_stiff_branched(rosenbrock_4stage_da, 2e-3); test_analytical_stiff_branched(rosenbrock_6stage_da, 2e-3); - test_analytical_stiff_branched, VectorStateType<1>>(rosenbrock_vector_1, 2e-3); - test_analytical_stiff_branched, VectorStateType<2>>(rosenbrock_vector_2, 2e-3); - test_analytical_stiff_branched, VectorStateType<3>>(rosenbrock_vector_3, 2e-3); - test_analytical_stiff_branched, VectorStateType<4>>(rosenbrock_vector_4, 2e-3); + test_analytical_stiff_branched(rosenbrock_vector_1, 2e-3); + test_analytical_stiff_branched(rosenbrock_vector_2, 2e-3); + test_analytical_stiff_branched(rosenbrock_vector_3, 2e-3); + test_analytical_stiff_branched(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); - - auto params = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - params.relative_tolerance_ = 1e-10; - params.absolute_tolerance_ = std::vector(3, params.relative_tolerance_ * 1e-2); - - auto standard_dootlittle = StandardRosenbrockDoolittle(params); - test_analytical_robertson(standard_dootlittle, 1e-1); - auto vector_1_dootlittle = VectorRosenbrockDoolittle<1>(params); - test_analytical_robertson, VectorStateTypeDoolittle<1>>(vector_1_dootlittle, 1e-1); - auto vector_2_dootlittle = VectorRosenbrockDoolittle<2>(params); - test_analytical_robertson, VectorStateTypeDoolittle<2>>(vector_2_dootlittle, 1e-1); - auto vector_3_dootlittle = VectorRosenbrockDoolittle<3>(params); - test_analytical_robertson, VectorStateTypeDoolittle<3>>(vector_3_dootlittle, 1e-1); - auto vector_4_dootlittle = VectorRosenbrockDoolittle<4>(params); - test_analytical_robertson, VectorStateTypeDoolittle<4>>(vector_4_dootlittle, 1e-1); - auto standard_mozart = StandardRosenbrockMozart(params); - test_analytical_robertson(standard_mozart, 1e-1); - auto vector_1_mozart = VectorRosenbrockMozart<1>(params); - test_analytical_robertson, VectorStateTypeMozart<1>>(vector_1_mozart, 1e-1); - auto vector_2_mozart = VectorRosenbrockMozart<2>(params); - test_analytical_robertson, VectorStateTypeMozart<2>>(vector_2_mozart, 1e-1); - auto vector_3_mozart = VectorRosenbrockMozart<3>(params); - test_analytical_robertson, VectorStateTypeMozart<3>>(vector_3_mozart, 1e-1); - auto vector_4_mozart = VectorRosenbrockMozart<4>(params); - test_analytical_robertson, VectorStateTypeMozart<4>>(vector_4_mozart, 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); - - auto params = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - params.relative_tolerance_ = 1e-13; - params.absolute_tolerance_ = std::vector(6, 1e-17); - params.absolute_tolerance_[0] = 1e-7; - params.absolute_tolerance_[4] = 1e-7; - params.absolute_tolerance_[5] = 1e-7; - - auto standard_dootlittle = StandardRosenbrockDoolittle(params); - test_analytical_e5(standard_dootlittle, 1e-3); - auto vector_1_dootlittle = VectorRosenbrockDoolittle<1>(params); - test_analytical_e5, VectorStateTypeDoolittle<1>>(vector_1_dootlittle, 1e-3); - auto vector_2_dootlittle = VectorRosenbrockDoolittle<2>(params); - test_analytical_e5, VectorStateTypeDoolittle<2>>(vector_2_dootlittle, 1e-3); - auto vector_3_dootlittle = VectorRosenbrockDoolittle<3>(params); - test_analytical_e5, VectorStateTypeDoolittle<3>>(vector_3_dootlittle, 1e-3); - auto vector_4_dootlittle = VectorRosenbrockDoolittle<4>(params); - test_analytical_e5, VectorStateTypeDoolittle<4>>(vector_4_dootlittle, 1e-3); - auto standard_mozart = StandardRosenbrockMozart(params); - test_analytical_e5(standard_mozart, 1e-3); - auto vector_1_mozart = VectorRosenbrockMozart<1>(params); - test_analytical_e5, VectorStateTypeMozart<1>>(vector_1_mozart, 1e-3); - auto vector_2_mozart = VectorRosenbrockMozart<2>(params); - test_analytical_e5, VectorStateTypeMozart<2>>(vector_2_mozart, 1e-3); - auto vector_3_mozart = VectorRosenbrockMozart<3>(params); - test_analytical_e5, VectorStateTypeMozart<3>>(vector_3_mozart, 1e-3); - auto vector_4_mozart = VectorRosenbrockMozart<4>(params); - test_analytical_e5, VectorStateTypeMozart<4>>(vector_4_mozart, 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_analytical_robertson(rosenbrock_vector_doolittle_1, 1e-1); + test_analytical_robertson(rosenbrock_vector_doolittle_2, 1e-1); + test_analytical_robertson(rosenbrock_vector_doolittle_3, 1e-1); + test_analytical_robertson(rosenbrock_vector_doolittle_4, 1e-1); + test_analytical_robertson(rosenbrock_standard_mozart, 1e-1); + test_analytical_robertson(rosenbrock_vector_mozart_1, 1e-1); + test_analytical_robertson(rosenbrock_vector_mozart_2, 1e-1); + test_analytical_robertson(rosenbrock_vector_mozart_3, 1e-1); + test_analytical_robertson(rosenbrock_vector_mozart_4, 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); - - auto params = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - params.relative_tolerance_ = 1e-8; - params.absolute_tolerance_ = std::vector(5, params.relative_tolerance_ * 1e-6); - - auto standard_dootlittle = StandardRosenbrockDoolittle(params); - test_analytical_oregonator(standard_dootlittle, 4e-6); - auto vector_1_dootlittle = VectorRosenbrockDoolittle<1>(params); - test_analytical_oregonator, VectorStateTypeDoolittle<1>>(vector_1_dootlittle, 4e-6); - auto vector_2_dootlittle = VectorRosenbrockDoolittle<2>(params); - test_analytical_oregonator, VectorStateTypeDoolittle<2>>(vector_2_dootlittle, 4e-6); - auto vector_3_dootlittle = VectorRosenbrockDoolittle<3>(params); - test_analytical_oregonator, VectorStateTypeDoolittle<3>>(vector_3_dootlittle, 4e-6); - auto vector_4_dootlittle = VectorRosenbrockDoolittle<4>(params); - test_analytical_oregonator, VectorStateTypeDoolittle<4>>(vector_4_dootlittle, 4e-6); - auto standard_mozart = StandardRosenbrockMozart(params); - test_analytical_oregonator(standard_mozart, 4e-6); - auto vector_1_mozart = VectorRosenbrockMozart<1>(params); - test_analytical_oregonator, VectorStateTypeMozart<1>>(vector_1_mozart, 4e-6); - auto vector_2_mozart = VectorRosenbrockMozart<2>(params); - test_analytical_oregonator, VectorStateTypeMozart<2>>(vector_2_mozart, 4e-6); - auto vector_3_mozart = VectorRosenbrockMozart<3>(params); - test_analytical_oregonator, VectorStateTypeMozart<3>>(vector_3_mozart, 4e-6); - auto vector_4_mozart = VectorRosenbrockMozart<4>(params); - test_analytical_oregonator, VectorStateTypeMozart<4>>(vector_4_mozart, 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_analytical_e5(rosenbrock_vector_doolittle_1, 1e-3); + test_analytical_e5(rosenbrock_vector_doolittle_2, 1e-3); + test_analytical_e5(rosenbrock_vector_doolittle_3, 1e-3); + test_analytical_e5(rosenbrock_vector_doolittle_4, 1e-3); + test_analytical_e5(rosenbrock_standard_mozart, 1e-3); + test_analytical_e5(rosenbrock_vector_mozart_1, 1e-3); + test_analytical_e5(rosenbrock_vector_mozart_2, 1e-3); + test_analytical_e5(rosenbrock_vector_mozart_3, 1e-3); + test_analytical_e5(rosenbrock_vector_mozart_4, 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_analytical_oregonator(rosenbrock_vector_doolittle_1, 4e-6); + test_analytical_oregonator(rosenbrock_vector_doolittle_2, 4e-6); + test_analytical_oregonator(rosenbrock_vector_doolittle_3, 4e-6); + test_analytical_oregonator(rosenbrock_vector_doolittle_4, 4e-6); + test_analytical_oregonator(rosenbrock_standard_mozart, 4e-6); + test_analytical_oregonator(rosenbrock_vector_mozart_1, 4e-6); + test_analytical_oregonator(rosenbrock_vector_mozart_2, 4e-6); + test_analytical_oregonator(rosenbrock_vector_mozart_3, 4e-6); + test_analytical_oregonator(rosenbrock_vector_mozart_4, 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); - - auto params = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - params.relative_tolerance_ = 1e-6; - params.absolute_tolerance_ = std::vector(8, params.relative_tolerance_ * 1e-2); - - auto standard_dootlittle = StandardRosenbrockDoolittle(params); - test_analytical_hires(standard_dootlittle, 1e-7); - auto vector_1_dootlittle = VectorRosenbrockDoolittle<1>(params); - test_analytical_hires, VectorStateTypeDoolittle<1>>(vector_1_dootlittle, 1e-7); - auto vector_2_dootlittle = VectorRosenbrockDoolittle<2>(params); - test_analytical_hires, VectorStateTypeDoolittle<2>>(vector_2_dootlittle, 1e-7); - auto vector_3_dootlittle = VectorRosenbrockDoolittle<3>(params); - test_analytical_hires, VectorStateTypeDoolittle<3>>(vector_3_dootlittle, 1e-7); - auto vector_4_dootlittle = VectorRosenbrockDoolittle<4>(params); - test_analytical_hires, VectorStateTypeDoolittle<4>>(vector_4_dootlittle, 1e-7); - auto standard_mozart = StandardRosenbrockMozart(params); - test_analytical_hires(standard_mozart, 1e-7); - auto vector_1_mozart = VectorRosenbrockMozart<1>(params); - test_analytical_hires, VectorStateTypeMozart<1>>(vector_1_mozart, 1e-7); - auto vector_2_mozart = VectorRosenbrockMozart<2>(params); - test_analytical_hires, VectorStateTypeMozart<2>>(vector_2_mozart, 1e-7); - auto vector_3_mozart = VectorRosenbrockMozart<3>(params); - test_analytical_hires, VectorStateTypeMozart<3>>(vector_3_mozart, 1e-7); - auto vector_4_mozart = VectorRosenbrockMozart<4>(params); - test_analytical_hires, VectorStateTypeMozart<4>>(vector_4_mozart, 1e-7); + 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); + test_analytical_hires(rosenbrock_vector_doolittle_1, 1e-7); + test_analytical_hires(rosenbrock_vector_doolittle_2, 1e-7); + test_analytical_hires(rosenbrock_vector_doolittle_3, 1e-7); + test_analytical_hires(rosenbrock_vector_doolittle_4, 1e-7); + test_analytical_hires(rosenbrock_standard_mozart, 1e-7); + test_analytical_hires(rosenbrock_vector_mozart_1, 1e-7); + test_analytical_hires(rosenbrock_vector_mozart_2, 1e-7); + test_analytical_hires(rosenbrock_vector_mozart_3, 1e-7); + test_analytical_hires(rosenbrock_vector_mozart_4, 1e-7); } diff --git a/test/integration/test_chapman_integration.cpp b/test/integration/test_chapman_integration.cpp index 0793e2d12..6bd8f0668 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.relative_tolerance_, 1.0e-4); + auto& abs_tol = state.absolute_tolerance_; 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..8a0d8035d 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.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,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.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) { @@ -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_lu_decomposition_policy.hpp b/test/unit/solver/test_lu_decomposition_policy.hpp index afea2f18f..2741b834d 100644 --- a/test/unit/solver/test_lu_decomposition_policy.hpp +++ b/test/unit/solver/test_lu_decomposition_policy.hpp @@ -113,8 +113,7 @@ void testDenseMatrix() A[0][2][1] = -2; A[0][2][2] = 8; - LuDecompositionPolicy lud = - LuDecompositionPolicy::template Create(A); + LuDecompositionPolicy lud = LuDecompositionPolicy::template Create(A); auto LU = LuDecompositionPolicy::template GetLUMatrices(A, 0); lud.template Decompose(A, LU.first, LU.second); check_results( diff --git a/test/unit/solver/test_rosenbrock.cpp b/test/unit/solver/test_rosenbrock.cpp index 357f2b683..c9f9170b1 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.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,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.absolute_tolerance_; + 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..80fdc148f 100644 --- a/test/unit/solver/test_solver_builder.cpp +++ b/test/unit/solver/test_solver_builder.cpp @@ -106,25 +106,66 @@ 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 options = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); + auto solver = micm::CpuSolverBuilder(options) + .SetSystem(the_system) + .SetReactions(reactions) + .SetNumberOfGridCells(1) + .Build(); + auto state = solver.GetState(); + state.variables_[0] = { 1.0, 0.0, 0.0 }; + + auto solve = solver.Solve(5, state); + + 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; + + state.variables_[0] = { 1.0, 0.0, 0.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); + + ASSERT_EQ(solver.solver_parameters_.h_min_, 15.0); + ASSERT_EQ(solver.solver_parameters_.max_number_of_steps_, 6.0); } \ No newline at end of file