From 100b3d63fb0d27116ccc92e572c3284a1eef3011 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Thu, 16 May 2024 15:24:15 -0700 Subject: [PATCH] update for new state policies --- docker/Dockerfile.llvm | 2 + include/micm/process/process.hpp | 24 +++--- include/micm/process/process_set.hpp | 75 ++++++++-------- include/micm/solver/rosenbrock.hpp | 8 +- include/micm/solver/rosenbrock.inl | 12 +-- include/micm/solver/solver_builder.hpp | 20 ++--- include/micm/solver/solver_builder.inl | 85 +++++++++++-------- include/micm/solver/state.hpp | 8 +- include/micm/solver/state.inl | 44 +++++----- test/integration/cmake/test_micm.cpp | 6 +- test/integration/e5.hpp | 4 +- test/integration/hires.hpp | 4 +- test/integration/oregonator.hpp | 4 +- test/tutorial/test_README_example.cpp | 51 ++++++++--- test/unit/process/test_jit_process_set.cpp | 15 ++-- test/unit/process/test_process.cpp | 16 ++-- test/unit/process/test_process_set.cpp | 26 +++--- test/unit/process/test_process_set_policy.hpp | 22 ++--- test/unit/solver/test_chapman_ode_solver.cpp | 8 +- 19 files changed, 235 insertions(+), 199 deletions(-) diff --git a/docker/Dockerfile.llvm b/docker/Dockerfile.llvm index 33b3272fa..005f24b1d 100644 --- a/docker/Dockerfile.llvm +++ b/docker/Dockerfile.llvm @@ -10,6 +10,7 @@ RUN dnf -y update \ make \ zlib-devel \ llvm-devel \ + openmpi-devel \ valgrind \ && dnf clean all @@ -24,6 +25,7 @@ RUN mkdir /build \ -D MICM_ENABLE_CLANG_TIDY:BOOL=FALSE \ -D MICM_ENABLE_LLVM:BOOL=TRUE \ -D MICM_ENABLE_MEMCHECK:BOOL=TRUE \ + -D MICM_ENABLE_OPENMP:BOOL=TRUE \ ../micm \ && make install -j 8 diff --git a/include/micm/process/process.hpp b/include/micm/process/process.hpp index 5d0379b06..5bc86cc5e 100644 --- a/include/micm/process/process.hpp +++ b/include/micm/process/process.hpp @@ -85,14 +85,14 @@ namespace micm /// @brief Update the solver state rate constants /// @param processes The set of processes being solved /// @param state The solver state to update - template class MatrixPolicy, class SparseMatrixPolicy> - requires(!VectorizableDense>) static void UpdateState( + template + requires(!VectorizableDense) static void UpdateState( const std::vector& processes, - State& state); - template class MatrixPolicy, class SparseMatrixPolicy> - requires(VectorizableDense>) static void UpdateState( + State& state); + template + requires(VectorizableDense) static void UpdateState( const std::vector& processes, - State& state); + State& state); friend class ProcessBuilder; static ProcessBuilder Create(); @@ -148,10 +148,10 @@ namespace micm ProcessBuilder& SetPhase(const Phase& phase); }; - template class MatrixPolicy, class SparseMatrixPolicy> - requires(!VectorizableDense>) void Process::UpdateState( + template + requires(!VectorizableDense) void Process::UpdateState( const std::vector& processes, - State& state) + State& state) { MICM_PROFILE_FUNCTION(); @@ -173,10 +173,10 @@ namespace micm } } - template class MatrixPolicy, class SparseMatrixPolicy> - requires(VectorizableDense>) void Process::UpdateState( + template + requires(VectorizableDense) void Process::UpdateState( const std::vector& processes, - State& state) + State& state) { MICM_PROFILE_FUNCTION(); const auto& v_custom_parameters = state.custom_rate_parameters_.AsVector(); diff --git a/include/micm/process/process_set.hpp b/include/micm/process/process_set.hpp index b495c69a7..31bfd1005 100644 --- a/include/micm/process/process_set.hpp +++ b/include/micm/process/process_set.hpp @@ -91,33 +91,32 @@ namespace micm /// @param rate_constants Current values for the process rate constants (grid cell, process) /// @param state_variables Current state variable values (grid cell, state variable) /// @param forcing Forcing terms for each state variable (grid cell, state variable) - template typename MatrixPolicy> - requires(!VectorizableDense>) void AddForcingTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - MatrixPolicy& forcing) const; - template typename MatrixPolicy> - requires VectorizableDense> + template + requires(!VectorizableDense) void AddForcingTerms( + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, + DenseMatrixPolicy& forcing) const; + template + requires VectorizableDense void AddForcingTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - MatrixPolicy& forcing) const; + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, + DenseMatrixPolicy& forcing) const; /// @brief Subtract Jacobian terms for the set of processes for the current conditions /// @param rate_constants Current values for the process rate constants (grid cell, process) /// @param state_variables Current state variable values (grid cell, state variable) /// @param jacobian Jacobian matrix for the system (grid cell, dependent variable, independent variable) - // template class MatrixPolicy, template class SparseMatrixPolicy> - template - requires(!VectorizableDense || !VectorizableSparse) void SubtractJacobianTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, + template + requires(!VectorizableDense || !VectorizableSparse) void SubtractJacobianTerms( + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, SparseMatrixPolicy& jacobian) const; // template class MatrixPolicy, template class SparseMatrixPolicy> - template - requires(VectorizableDense&& VectorizableSparse) void SubtractJacobianTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, + template + requires(VectorizableDense&& VectorizableSparse) void SubtractJacobianTerms( + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, SparseMatrixPolicy& jacobian) const; /// @brief Returns the set of species used in a set of processes @@ -217,11 +216,11 @@ namespace micm } } - template typename MatrixPolicy> - requires(!VectorizableDense>) inline void ProcessSet::AddForcingTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - MatrixPolicy& forcing) const + template + requires(!VectorizableDense) inline void ProcessSet::AddForcingTerms( + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, + DenseMatrixPolicy& forcing) const { MICM_PROFILE_FUNCTION(); @@ -255,12 +254,12 @@ namespace micm } }; - template typename MatrixPolicy> - requires VectorizableDense> + template + requires VectorizableDense inline void ProcessSet::AddForcingTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - MatrixPolicy& forcing) const + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, + DenseMatrixPolicy& forcing) const { MICM_PROFILE_FUNCTION(); @@ -300,12 +299,11 @@ namespace micm } // Forming the Jacobian matrix "J" and returning "-J" to be consistent with the CUDA implementation - // template class MatrixPolicy, template class SparseMatrixPolicy> - template - requires(!VectorizableDense || !VectorizableSparse) inline void ProcessSet:: + template + requires(!VectorizableDense || !VectorizableSparse) inline void ProcessSet:: SubtractJacobianTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, SparseMatrixPolicy& jacobian) const { MICM_PROFILE_FUNCTION(); @@ -350,12 +348,11 @@ namespace micm } // Forming the Jacobian matrix "J" and returning "-J" to be consistent with the CUDA implementation - // template class MatrixPolicy, template class SparseMatrixPolicy> - template - requires(VectorizableDense&& VectorizableSparse) inline void ProcessSet:: + template + requires(VectorizableDense&& VectorizableSparse) inline void ProcessSet:: SubtractJacobianTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, SparseMatrixPolicy& jacobian) const { MICM_PROFILE_FUNCTION(); diff --git a/include/micm/solver/rosenbrock.hpp b/include/micm/solver/rosenbrock.hpp index c50c03b72..12bda0117 100644 --- a/include/micm/solver/rosenbrock.hpp +++ b/include/micm/solver/rosenbrock.hpp @@ -149,12 +149,12 @@ namespace micm /// @brief Returns a state object for use with the solver /// @return A object that can hold the full state of the chemical system - virtual State GetState() const; + virtual State, SparseMatrixPolicy> GetState() const; /// @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, State& state) noexcept; + SolverResult Solve(double time_step, State, SparseMatrixPolicy>& state) noexcept; /// @brief Calculate a chemical forcing /// @param rate_constants List of rate constants for each needed species @@ -175,7 +175,7 @@ namespace micm /// @brief Update the rate constants for the environment state /// @param state The current state of the chemical system - void UpdateState(State& state); + void UpdateState(State, SparseMatrixPolicy>& state); /// @brief Compute the derivative of the forcing w.r.t. each chemical, and return the negative jacobian /// @param rate_constants List of rate constants for each needed species @@ -198,7 +198,7 @@ namespace micm bool& singular, const MatrixPolicy& number_densities, SolverStats& stats, - State& state); + State, SparseMatrixPolicy>& state); /// @brief Computes the scaled norm of the vector errors /// @param y the original vector diff --git a/include/micm/solver/rosenbrock.inl b/include/micm/solver/rosenbrock.inl index 351647174..93a7c3b56 100644 --- a/include/micm/solver/rosenbrock.inl +++ b/include/micm/solver/rosenbrock.inl @@ -237,10 +237,10 @@ namespace micm class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> - inline State + inline State, SparseMatrixPolicy> RosenbrockSolver::GetState() const { - return State{ state_parameters_ }; + return State, SparseMatrixPolicy>{ state_parameters_ }; } template< @@ -252,7 +252,7 @@ namespace micm inline typename RosenbrockSolver::SolverResult RosenbrockSolver::Solve( double time_step, - State& state) noexcept + State, SparseMatrixPolicy>& state) noexcept { MICM_PROFILE_FUNCTION(); @@ -460,7 +460,7 @@ namespace micm { MICM_PROFILE_FUNCTION(); std::fill(forcing.AsVector().begin(), forcing.AsVector().end(), 0.0); - process_set_.template AddForcingTerms(rate_constants, number_densities, forcing); + process_set_.template AddForcingTerms>(rate_constants, number_densities, forcing); } template< @@ -530,7 +530,7 @@ namespace micm class LinearSolverPolicy, class ProcessSetPolicy> inline void RosenbrockSolver::UpdateState( - State& state) + State, SparseMatrixPolicy>& state) { Process::UpdateState(processes_, state); } @@ -547,7 +547,7 @@ namespace micm bool& singular, const MatrixPolicy& number_densities, SolverStats& stats, - State& state) + State, SparseMatrixPolicy>& state) { MICM_PROFILE_FUNCTION(); diff --git a/include/micm/solver/solver_builder.hpp b/include/micm/solver/solver_builder.hpp index 1e887a892..bf28b6f40 100644 --- a/include/micm/solver/solver_builder.hpp +++ b/include/micm/solver/solver_builder.hpp @@ -20,6 +20,7 @@ namespace micm { + template class SolverBuilder { protected: @@ -58,13 +59,12 @@ namespace micm /// @brief /// @return - template - Solver> Build(); + Solver> Build(); protected: /// @brief /// @return - virtual Solver BuildBackwardEulerSolver() = 0; + virtual Solver> BuildBackwardEulerSolver() = 0; virtual ~SolverBuilder() = default; /// @brief @@ -75,7 +75,7 @@ namespace micm /// @brief Get a species map properly ordered /// @return - template + template std::map GetSpeciesMap() const; /// @brief Set the absolute tolerances per species @@ -91,18 +91,18 @@ namespace micm std::vector GetJacobianDiagonalElements(auto jacobian) const; }; - template - class CpuSolverBuilder : public SolverBuilder + template + class CpuSolverBuilder : public SolverBuilder { public: - Solver BuildBackwardEulerSolver() override; + Solver> BuildBackwardEulerSolver() override; }; - template - class GpuSolverBuilder : public SolverBuilder + template + class GpuSolverBuilder : public SolverBuilder { public: - Solver BuildBackwardEulerSolver() override; + Solver> BuildBackwardEulerSolver() override; }; } // namespace micm diff --git a/include/micm/solver/solver_builder.inl b/include/micm/solver/solver_builder.inl index cf830e2a7..457520670 100644 --- a/include/micm/solver/solver_builder.inl +++ b/include/micm/solver/solver_builder.inl @@ -47,25 +47,29 @@ inline std::error_code make_error_code(MicmSolverBuilderErrc e) namespace micm { - inline SolverBuilder& SolverBuilder::SetSystem(System system) + template + inline SolverBuilder& SolverBuilder::SetSystem(System system) { system_ = system; return *this; } - inline SolverBuilder& SolverBuilder::SetReactions(std::vector reactions) + template + inline SolverBuilder& SolverBuilder::SetReactions(std::vector reactions) { reactions_ = reactions; return *this; } - inline SolverBuilder& SolverBuilder::SetNumberOfGridCells(int number_of_grid_cells) + template + inline SolverBuilder& SolverBuilder::SetNumberOfGridCells(int number_of_grid_cells) { number_of_grid_cells_ = number_of_grid_cells; return *this; } - inline SolverBuilder& SolverBuilder::SetSolverParameters(const RosenbrockSolverParameters& options) + template + inline SolverBuilder& SolverBuilder::SetSolverParameters(const RosenbrockSolverParameters& options) { if (!std::holds_alternative(options_)) { @@ -75,7 +79,8 @@ namespace micm return *this; } - inline SolverBuilder& SolverBuilder::SetSolverParameters(const BackwardEulerSolverParameters& options) + template + inline SolverBuilder& SolverBuilder::SetSolverParameters(const BackwardEulerSolverParameters& options) { if (!std::holds_alternative(options_)) { @@ -85,7 +90,8 @@ namespace micm return *this; } - inline Solver SolverBuilder::Build() + template + inline Solver> SolverBuilder::Build() { if (std::holds_alternative(options_)) { @@ -99,8 +105,9 @@ namespace micm throw std::runtime_error("No solver type set"); } + template template - inline void SolverBuilder::UnusedSpeciesCheck() + inline void SolverBuilder::UnusedSpeciesCheck() { if (ignore_unused_species_) { @@ -128,8 +135,9 @@ namespace micm } } - template - inline std::map SolverBuilder::GetSpeciesMap() const + template + template + inline std::map SolverBuilder::GetSpeciesMap() const { std::map species_map; std::function& variables, const std::size_t i)> state_reordering; @@ -143,7 +151,7 @@ namespace micm auto unsorted_process_set = ProcessSetPolicy(reactions_, species_map); auto unsorted_jac_elements = unsorted_process_set.NonZeroJacobianElements(); - using Matrix = typename MatrixPolicy::IntMatrix; + using Matrix = typename DenseMatrixPolicy::IntMatrix; Matrix unsorted_jac_non_zeros(system_.StateSize(), system_.StateSize(), 0); for (auto& elem : unsorted_jac_elements) unsorted_jac_non_zeros[elem.first][elem.second] = 1; @@ -160,7 +168,8 @@ namespace micm return species_map; } - inline void SolverBuilder::SetAbsoluteTolerances( + template + inline void SolverBuilder::SetAbsoluteTolerances( std::vector& tolerances, const std::map& species_map) const { @@ -188,7 +197,8 @@ namespace micm } } - inline std::vector SolverBuilder::GetCustomParameterLabels() const + template + inline std::vector SolverBuilder::GetCustomParameterLabels() const { std::vector param_labels{}; for (const auto& reaction : reactions_) @@ -198,7 +208,8 @@ namespace micm return param_labels; } - inline std::vector SolverBuilder::GetJacobianDiagonalElements(auto jacobian) const + template + inline std::vector SolverBuilder::GetJacobianDiagonalElements(auto jacobian) const { std::vector jacobian_diagonal_elements; @@ -212,41 +223,45 @@ namespace micm return jacobian_diagonal_elements; } - template - inline Solver CpuSolverBuilder::BuildBackwardEulerSolver() + template + inline Solver> CpuSolverBuilder::BuildBackwardEulerSolver() { using ProcessSetPolicy = ProcessSet; + using LinearSolverPolicy = LinearSolver; - auto parameters = std::get(options_); - auto species_map = GetSpeciesMap(); - auto labels = GetCustomParameterLabels(); - std::size_t number_of_species = system_.StateSize(); + auto parameters = std::get(this->options_); + auto species_map = this->template GetSpeciesMap(); + auto labels = this->GetCustomParameterLabels(); + std::size_t number_of_species = this->system_.StateSize(); - UnusedSpeciesCheck(); - SetAbsoluteTolerances(parameters.absolute_tolerance_, species_map); + this->template UnusedSpeciesCheck(); + this->SetAbsoluteTolerances(parameters.absolute_tolerance_, species_map); - ProcessSetPolicy process_set(reactions_, species_map); + ProcessSetPolicy process_set(this->reactions_, species_map); auto jacobian = - BuildJacobian(process_set.NonZeroJacobianElements(), number_of_grid_cells_, number_of_species); - auto diagonal_elements = GetJacobianDiagonalElements(jacobian); + BuildJacobian(process_set.NonZeroJacobianElements(), this->number_of_grid_cells_, number_of_species); process_set.SetJacobianFlatIds(jacobian); - LinearSolver linear_solver(jacobian, 1e-30); + LinearSolverPolicy linear_solver(jacobian, 1e-30); + + std::vector variable_names{number_of_species}; + for (auto& species_pair : species_map) + variable_names[species_pair.second] = species_pair.first; - StateParameters state_parameters_ = { .number_of_grid_cells_ = parameters_.number_of_grid_cells_, - .number_of_species_ = variable_map.size(), - .number_of_rate_constants_ = processes_.size(), - .variable_names_ = system.UniqueNames(state_reordering), - .custom_rate_parameter_labels_ = param_labels, - .jacobian_diagonal_elements_ = jacobian_diagonal_elements, - .nonzero_jacobian_elements_ = process_set_.NonZeroJacobianElements() }; + StateParameters state_parameters_ = { .number_of_grid_cells_ = this->number_of_grid_cells_, + .number_of_species_ = number_of_species, + .number_of_rate_constants_ = this->reactions_.size(), + .variable_names_ = variable_names, + .custom_rate_parameter_labels_ = this->GetCustomParameterLabels(), + .jacobian_diagonal_elements_ = this->GetJacobianDiagonalElements(jacobian), + .nonzero_jacobian_elements_ = process_set.NonZeroJacobianElements() }; - return Solver>( + return Solver>( new SolverImpl(), parameters, - number_of_grid_cells_, + this->number_of_grid_cells_, number_of_species, - reactions_.size()); + this->reactions_.size()); } } // namespace micm \ No newline at end of file diff --git a/include/micm/solver/state.hpp b/include/micm/solver/state.hpp index 3aaad6753..4ff66d41d 100644 --- a/include/micm/solver/state.hpp +++ b/include/micm/solver/state.hpp @@ -35,15 +35,15 @@ namespace micm std::set> nonzero_jacobian_elements_{}; }; - template + template struct State { /// @brief The concentration of chemicals, varies through time - MatrixPolicy variables_; + DenseMatrixPolicy variables_; /// @brief Rate paramters particular to user-defined rate constants, may vary in time - MatrixPolicy custom_rate_parameters_; + DenseMatrixPolicy custom_rate_parameters_; /// @brief The reaction rates, may vary in time - MatrixPolicy rate_constants_; + DenseMatrixPolicy rate_constants_; /// @brief Atmospheric conditions, varies in time std::vector conditions_; /// @brief The jacobian structure, varies for each solve diff --git a/include/micm/solver/state.inl b/include/micm/solver/state.inl index 10779a18f..186b7d0f8 100644 --- a/include/micm/solver/state.inl +++ b/include/micm/solver/state.inl @@ -58,8 +58,8 @@ inline std::error_code make_error_code(MicmStateErrc e) namespace micm { - template - inline State::State() + template + inline State::State() : conditions_(), variables_(), custom_rate_parameters_(), @@ -68,8 +68,8 @@ namespace micm { } - template - inline State::State(const StateParameters& parameters) + template + inline State::State(const StateParameters& parameters) : conditions_(parameters.number_of_grid_cells_), variables_(parameters.number_of_grid_cells_, parameters.variable_names_.size(), 0.0), custom_rate_parameters_(parameters.number_of_grid_cells_, parameters.custom_rate_parameter_labels_.size(), 0.0), @@ -100,8 +100,8 @@ namespace micm upper_matrix_ = upper_matrix; } - template - inline void State::SetConcentrations( + template + inline void State::SetConcentrations( const std::unordered_map>& species_to_concentration) { const std::size_t num_grid_cells = conditions_.size(); @@ -109,8 +109,8 @@ namespace micm SetConcentration({ pair.first }, pair.second); } - template - inline void State::SetConcentration(const Species& species, double concentration) + template + inline void State::SetConcentration(const Species& species, double concentration) { auto var = variable_map_.find(species.name_); if (var == variable_map_.end()) @@ -120,8 +120,8 @@ namespace micm variables_[0][variable_map_[species.name_]] = concentration; } - template - inline void State::SetConcentration( + template + inline void State::SetConcentration( const Species& species, const std::vector& concentration) { @@ -135,8 +135,8 @@ namespace micm variables_[i][i_species] = concentration[i]; } - template - inline void State::UnsafelySetCustomRateParameters( + template + inline void State::UnsafelySetCustomRateParameters( const std::vector>& parameters) { if (parameters.size() != variables_.NumRows()) @@ -152,16 +152,16 @@ namespace micm } } - template - inline void State::SetCustomRateParameters( + template + inline void State::SetCustomRateParameters( const std::unordered_map>& parameters) { for (auto& pair : parameters) SetCustomRateParameter(pair.first, pair.second); } - template - inline void State::SetCustomRateParameter(const std::string& label, double value) + template + inline void State::SetCustomRateParameter(const std::string& label, double value) { auto param = custom_rate_parameter_map_.find(label); if (param == custom_rate_parameter_map_.end()) @@ -172,8 +172,8 @@ namespace micm custom_rate_parameters_[0][param->second] = value; } - template - inline void State::SetCustomRateParameter( + template + inline void State::SetCustomRateParameter( const std::string& label, const std::vector& values) { @@ -187,8 +187,8 @@ namespace micm custom_rate_parameters_[i][param->second] = values[i]; } - template - inline void State::PrintHeader() + template + inline void State::PrintHeader() { auto largest_str_iter = std::max_element( variable_names_.begin(), variable_names_.end(), [](const auto& a, const auto& b) { return a.size() < b.size(); }); @@ -208,8 +208,8 @@ namespace micm std::cout << std::endl; } - template class MatrixPolicy, class SparseMatrixPolicy> - inline void State::PrintState(double time) + template + inline void State::PrintState(double time) { std::ios oldState(nullptr); oldState.copyfmt(std::cout); diff --git a/test/integration/cmake/test_micm.cpp b/test/integration/cmake/test_micm.cpp index ccd02e835..0310bbc00 100644 --- a/test/integration/cmake/test_micm.cpp +++ b/test/integration/cmake/test_micm.cpp @@ -11,8 +11,8 @@ void print_header() << "," << std::setw(10) << "C" << std::endl; } -template class T> -void print_state(double time, State& state) +template +void print_state(double time, StatePolicy& state) { std::ios oldState(nullptr); oldState.copyfmt(std::cout); @@ -42,7 +42,7 @@ int main() auto reactions = solver_params.processes_; RosenbrockSolver<> solver{ chemical_system, reactions, params }; - State state = solver.GetState(); + auto state = solver.GetState(); // mol m-3 state.variables_[0] = { 1, 0, 0 }; diff --git a/test/integration/e5.hpp b/test/integration/e5.hpp index 931bad40b..09e97a2b9 100644 --- a/test/integration/e5.hpp +++ b/test/integration/e5.hpp @@ -61,9 +61,9 @@ class E5 : public micm::RosenbrockSolver GetState() const override + micm::State, SparseMatrixPolicy> GetState() const override { - auto state = micm::State{ this->state_parameters_ }; + auto state = micm::State, SparseMatrixPolicy>{ this->state_parameters_ }; state.jacobian_ = micm::BuildJacobian( nonzero_jacobian_elements_, diff --git a/test/integration/hires.hpp b/test/integration/hires.hpp index 7f97aa3f5..66c036706 100644 --- a/test/integration/hires.hpp +++ b/test/integration/hires.hpp @@ -61,9 +61,9 @@ class HIRES : public micm::RosenbrockSolver GetState() const override + micm::State, SparseMatrixPolicy> GetState() const override { - auto state = micm::State{ this->state_parameters_ }; + auto state = micm::State, SparseMatrixPolicy>{ this->state_parameters_ }; state.jacobian_ = micm::BuildJacobian( nonzero_jacobian_elements_, diff --git a/test/integration/oregonator.hpp b/test/integration/oregonator.hpp index c64627a67..707351974 100644 --- a/test/integration/oregonator.hpp +++ b/test/integration/oregonator.hpp @@ -62,9 +62,9 @@ class Oregonator : public micm::RosenbrockSolver GetState() const override + micm::State, SparseMatrixPolicy> GetState() const override { - auto state = micm::State{ this->state_parameters_ }; + auto state = micm::State, SparseMatrixPolicy>{ this->state_parameters_ }; state.jacobian_ = micm::BuildJacobian( nonzero_jacobian_elements_, diff --git a/test/tutorial/test_README_example.cpp b/test/tutorial/test_README_example.cpp index 099a7f830..4be9fbb66 100644 --- a/test/tutorial/test_README_example.cpp +++ b/test/tutorial/test_README_example.cpp @@ -30,21 +30,52 @@ int main(const int argc, const char *argv[]) std::vector reactions{ r1, r2 }; - RosenbrockSolver<> solver{ chemical_system, reactions, RosenbrockSolverParameters::ThreeStageRosenbrockParameters() }; +#if 0 + auto rosenbrock_solver = Solver::Create() + .SetSystem(chemical_system) + .SetReactions(reactions) + .SetNumberOfGridCells(19) + .Rosenbrock(RosenbrockSolverParameters::ThreeStageRosenbrockParameters()); - State state = solver.GetState(); + auto backward_euler_solver = Solver::Create() + .SetSystem(chemical_system) + .SetReactions(reactions) + .SetNumberOfGridCells(19) + .BackwardEuler(BackwardEulerParameters::CamChemParameters()); - state.conditions_[0].temperature_ = 287.45; // K - state.conditions_[0].pressure_ = 101319.9; // Pa - state.SetConcentration(foo, 20.0); // mol m-3 + auto cuda_solver = Solver::Create() + .SetSystem(chemical_system) + .SetReactions(reactions) + .SetNumberOfGridCells(19) + .Rosenbrock(RosenbrockSolverParameters::ThreeStageRosenbrockParameter()); + - state.PrintHeader(); + auto rosenbrock_state = rosenbrock_solver.GetState(); + auto backward_euler_state = backward_euler_solver.GetState(); + auto cuda_state = cuda_solver.GetState(); + + rosenbrock_state.conditions_[0].temperature_ = 287.45; // K + rosenbrock_state.conditions_[0].pressure_ = 101319.9; // Pa + rosenbrock_state.SetConcentration(foo, 20.0); // mol m-3 + + backward_euler_state.conditions_[0].temperature_ = 287.45; // K + backward_euler_state.conditions_[0].pressure_ = 101319.9; // Pa + backward_euler_state.SetConcentration(foo, 20.0); // mol m-3 + + cuda_state.conditions_[0].temperature_ = 287.45; // K + cuda_state.conditions_[0].pressure_ = 101319.9; // Pa + cuda_state.SetConcentration(foo, 20.0); // mol m-3 + + cuda_state.PrintHeader(); for (int i = 0; i < 10; ++i) { - auto result = solver.Solve(500.0, state); - state.variables_ = result.result_; - state.PrintState(i * 500); + cuda_state.CalculateRateConstants(); // calculate rate constants on the host + cuda_state.SyncToDevice(); + auto result = cuda_solver.Solve(500.0, cuda_state); // solve on the device + cuda_state.variables_ = result.result_; // device to device copy with overloaded operator= + cuda_state.SyncToHost(); + cuda_state.PrintState(i * 500); } - +#endif return 0; } \ No newline at end of file diff --git a/test/unit/process/test_jit_process_set.cpp b/test/unit/process/test_jit_process_set.cpp index 3f2e43864..5128e402a 100644 --- a/test/unit/process/test_jit_process_set.cpp +++ b/test/unit/process/test_jit_process_set.cpp @@ -9,17 +9,12 @@ #include -template -using Group2VectorMatrix = micm::VectorMatrix; -template -using Group200VectorMatrix = micm::VectorMatrix; -template -using Group300VectorMatrix = micm::VectorMatrix; -template -using Group400VectorMatrix = micm::VectorMatrix; +using Group2VectorMatrix = micm::VectorMatrix; +using Group200VectorMatrix = micm::VectorMatrix; +using Group300VectorMatrix = micm::VectorMatrix; +using Group400VectorMatrix = micm::VectorMatrix; -template -using Group2SparseVectorMatrix = micm::SparseMatrix>; +using Group2SparseVectorMatrix = micm::SparseMatrix>; TEST(JitProcessSet, VectorMatrix) { diff --git a/test/unit/process/test_process.cpp b/test/unit/process/test_process.cpp index 1061846a6..b4a7e4f01 100644 --- a/test/unit/process/test_process.cpp +++ b/test/unit/process/test_process.cpp @@ -9,7 +9,7 @@ #include -template class MatrixPolicy> +template void testProcessUpdateState(const std::size_t number_of_grid_cells) { micm::Species foo("foo", { { "molecular weight [kg mol-1]", 0.025 }, { "diffusion coefficient [m2 s-1]", 2.3e2 } }); @@ -29,14 +29,14 @@ void testProcessUpdateState(const std::size_t number_of_grid_cells) for (const auto& process : processes) for (auto& label : process.rate_constant_->CustomParameters()) param_labels.push_back(label); - micm::State state{ micm::StateParameters{ + micm::State state{ micm::StateParameters{ .number_of_grid_cells_ = number_of_grid_cells, .number_of_rate_constants_ = processes.size(), .variable_names_ = { "foo", "bar'" }, .custom_rate_parameter_labels_ = param_labels, } }; - MatrixPolicy expected_rate_constants(number_of_grid_cells, 3, 0.0); + DenseMatrixPolicy expected_rate_constants(number_of_grid_cells, 3, 0.0); std::vector params = { 0.0, 0.0, 0.0 }; auto get_double = std::bind(std::lognormal_distribution(0.0, 0.01), std::default_random_engine()); @@ -82,15 +82,15 @@ using Group4VectorMatrix = micm::VectorMatrix; TEST(Process, Matrix) { - testProcessUpdateState(5); + testProcessUpdateState>(5); } TEST(Process, VectorMatrix) { - testProcessUpdateState(5); - testProcessUpdateState(5); - testProcessUpdateState(5); - testProcessUpdateState(5); + testProcessUpdateState>(5); + testProcessUpdateState>(5); + testProcessUpdateState>(5); + testProcessUpdateState>(5); } TEST(Process, SurfaceRateConstantOnlyHasOneReactant) diff --git a/test/unit/process/test_process_set.cpp b/test/unit/process/test_process_set.cpp index ab866148b..a1c0cbd81 100644 --- a/test/unit/process/test_process_set.cpp +++ b/test/unit/process/test_process_set.cpp @@ -13,14 +13,10 @@ using SparseMatrixTest = micm::SparseMatrix; -template -using Group1VectorMatrix = micm::VectorMatrix; -template -using Group2VectorMatrix = micm::VectorMatrix; -template -using Group3VectorMatrix = micm::VectorMatrix; -template -using Group4VectorMatrix = micm::VectorMatrix; +using Group1VectorMatrix = micm::VectorMatrix; +using Group2VectorMatrix = micm::VectorMatrix; +using Group3VectorMatrix = micm::VectorMatrix; +using Group4VectorMatrix = micm::VectorMatrix; using Group1SparseVectorMatrix = micm::SparseMatrix>; using Group2SparseVectorMatrix = micm::SparseMatrix>; @@ -29,7 +25,7 @@ using Group4SparseVectorMatrix = micm::SparseMatrix(); + testProcessSet, SparseMatrixTest, micm::ProcessSet>(); } TEST(ProcessSet, VectorMatrix) @@ -42,28 +38,28 @@ TEST(ProcessSet, VectorMatrix) TEST(RandomProcessSet, Matrix) { - testRandomSystem( + testRandomSystem, SparseMatrixTest, micm::ProcessSet>( 200, 50, 40, [](const std::vector& processes, - const micm::State& state) -> micm::ProcessSet { + const micm::State, SparseMatrixTest>& state) -> micm::ProcessSet { return micm::ProcessSet{ processes, state.variable_map_ }; }); - testRandomSystem( + testRandomSystem, SparseMatrixTest, micm::ProcessSet>( 300, 30, 20, [](const std::vector& processes, - const micm::State& state) -> micm::ProcessSet { + const micm::State, SparseMatrixTest>& state) -> micm::ProcessSet { return micm::ProcessSet{ processes, state.variable_map_ }; }); - testRandomSystem( + testRandomSystem, SparseMatrixTest, micm::ProcessSet>( 400, 100, 80, [](const std::vector& processes, - const micm::State& state) -> micm::ProcessSet { + const micm::State, SparseMatrixTest>& state) -> micm::ProcessSet { return micm::ProcessSet{ processes, state.variable_map_ }; }); } \ No newline at end of file diff --git a/test/unit/process/test_process_set_policy.hpp b/test/unit/process/test_process_set_policy.hpp index f000cf195..f32d89cf8 100644 --- a/test/unit/process/test_process_set_policy.hpp +++ b/test/unit/process/test_process_set_policy.hpp @@ -13,7 +13,7 @@ void compare_pair(const index_pair& a, const index_pair& b) EXPECT_EQ(a.second, b.second); } -template class MatrixPolicy, class SparseMatrixPolicy, class ProcessSetPolicy> +template void testProcessSet() { auto foo = micm::Species("foo"); @@ -27,7 +27,7 @@ void testProcessSet() micm::Phase gas_phase{ std::vector{ foo, bar, qux, baz, quz, quuz, corge } }; - micm::State state( + micm::State state( micm::StateParameters{ .number_of_grid_cells_ = 2, .number_of_rate_constants_ = 3, .variable_names_{ "foo", "bar", "baz", "quz", "quuz", "corge" } }); @@ -61,13 +61,13 @@ void testProcessSet() EXPECT_EQ(state.variables_.NumColumns(), 6); state.variables_[0] = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.0 }; state.variables_[1] = { 1.1, 1.2, 1.3, 1.4, 1.5, 0.0 }; - MatrixPolicy rate_constants{ 2, 3 }; + DenseMatrixPolicy rate_constants{ 2, 3 }; rate_constants[0] = { 10.0, 20.0, 30.0 }; rate_constants[1] = { 110.0, 120.0, 130.0 }; - MatrixPolicy forcing{ 2, 5, 1000.0 }; + DenseMatrixPolicy forcing{ 2, 5, 1000.0 }; - set.template AddForcingTerms(rate_constants, state.variables_, forcing); + set.template AddForcingTerms(rate_constants, state.variables_, forcing); EXPECT_EQ(forcing[0][0], 1000.0 - 10.0 * 0.1 * 0.3 + 20.0 * 0.2); EXPECT_EQ(forcing[1][0], 1000.0 - 110.0 * 1.1 * 1.3 + 120.0 * 1.2); EXPECT_EQ(forcing[0][1], 1000.0 + 10.0 * 0.1 * 0.3 - 20.0 * 0.2); @@ -133,13 +133,13 @@ void testProcessSet() EXPECT_DOUBLE_EQ(jacobian[1][4][2], 100.0 - 2.4 * 110.0 * 1.1); } -template class MatrixPolicy, class SparseMatrixPolicy, class ProcessSetPolicy> +template void testRandomSystem( std::size_t n_cells, std::size_t n_reactions, std::size_t n_species, const std::function< - ProcessSetPolicy(const std::vector&, const micm::State&)> + ProcessSetPolicy(const std::vector&, const micm::State&)> create_set) { auto get_n_react = std::bind(std::uniform_int_distribution<>(0, 3), std::default_random_engine()); @@ -155,7 +155,7 @@ void testRandomSystem( species_names.push_back(std::to_string(i)); } micm::Phase gas_phase{ species }; - micm::State state{ micm::StateParameters{ + micm::State state{ micm::StateParameters{ .number_of_grid_cells_ = n_cells, .number_of_rate_constants_ = n_reactions, .variable_names_{ species_names }, @@ -183,10 +183,10 @@ void testRandomSystem( for (auto& elem : state.variables_.AsVector()) elem = get_double(); - MatrixPolicy rate_constants{ n_cells, n_reactions }; + DenseMatrixPolicy rate_constants{ n_cells, n_reactions }; for (auto& elem : rate_constants.AsVector()) elem = get_double(); - MatrixPolicy forcing{ n_cells, n_species, 1000.0 }; + DenseMatrixPolicy forcing{ n_cells, n_species, 1000.0 }; - set.template AddForcingTerms(rate_constants, state.variables_, forcing); + set.template AddForcingTerms(rate_constants, state.variables_, forcing); } diff --git a/test/unit/solver/test_chapman_ode_solver.cpp b/test/unit/solver/test_chapman_ode_solver.cpp index 07740598e..c27cc0437 100644 --- a/test/unit/solver/test_chapman_ode_solver.cpp +++ b/test/unit/solver/test_chapman_ode_solver.cpp @@ -210,7 +210,7 @@ TEST(ChapmanMechanismHardCodedAndGeneral, dforce_dy_times_vector) void TestSolve(micm::ChapmanODESolver& solver) { - micm::State state = solver.GetState(); + micm::State> state = solver.GetState(); state.variables_[0] = { 1, 3.92e-1, 1.69e-2, 0, 3.29e1, 0, 0, 8.84, 0 }; //"M" "Ar" "CO2", "H2O", "N2", "O1D", "O", "O2", "O3", state.conditions_[0].temperature_ = 273.15; @@ -245,7 +245,7 @@ TEST(ChapmanMechanismHardCodedAndGeneral, Solve) void TestSolve10TimesLarger(micm::ChapmanODESolver& solver) { - micm::State state = solver.GetState(); + micm::State> state = solver.GetState(); state.variables_[0] = { 1, 3.92e-1, 1.69e-2, 0, 3.29e1, 0, 0, 8.84, 0 }; //"M" "Ar" "CO2", "H2O", "N2", "O1D", "O", "O2", "O3", state.conditions_[0].temperature_ = 273.15; @@ -283,7 +283,7 @@ TEST(ChapmanMechanismHardCodedAndGeneral, solve_10_times_larger) void TestSolve10TimesSmaller(micm::ChapmanODESolver& solver) { - micm::State state = solver.GetState(); + micm::State> state = solver.GetState(); state.variables_[0] = { 1, 3.92e-1, 1.69e-2, 0, 3.29e1, 0, 0, 8.84, 0 }; //"M" "Ar" "CO2", "H2O", "N2", "O1D", "O", "O2", "O3", state.conditions_[0].temperature_ = 273.15; @@ -321,7 +321,7 @@ TEST(RegressionChapmanODESolver, solve_10_times_smaller) void TestSolveWithRandomNumberDensities(micm::ChapmanODESolver& solver) { - micm::State state = solver.GetState(); + micm::State> state = solver.GetState(); state.conditions_[0].temperature_ = 273.15; state.conditions_[0].pressure_ = 1000 * 100; // 1000 hPa state.conditions_[0].air_density_ = 2.7e19;