diff --git a/include/micm/process/process.hpp b/include/micm/process/process.hpp index cd7aee650..786f4bda8 100644 --- a/include/micm/process/process.hpp +++ b/include/micm/process/process.hpp @@ -36,13 +36,11 @@ namespace micm /// @param processes The set of processes being solved /// @param state The solver state to update template class MatrixPolicy> - requires(!VectorizableDense>) static void UpdateState( - const std::vector& processes, - State& state); + requires(!VectorizableDense>) + static void UpdateState(const std::vector& processes, State& state); template class MatrixPolicy> - requires(VectorizableDense>) static void UpdateState( - const std::vector& processes, - State& state); + requires(VectorizableDense>) + static void UpdateState(const std::vector& processes, State& state); friend class ProcessBuilder; static ProcessBuilder create(); @@ -92,9 +90,8 @@ namespace micm }; template class MatrixPolicy> - requires(!VectorizableDense>) void Process::UpdateState( - const std::vector& processes, - State& state) + requires(!VectorizableDense>) + void Process::UpdateState(const std::vector& processes, State& state) { for (std::size_t i{}; i < state.custom_rate_parameters_.size(); ++i) { @@ -103,17 +100,20 @@ namespace micm std::size_t i_rate_constant = 0; for (auto& process : processes) { + double fixed_reactants = 1.0; + for (auto& reactant : process.reactants_) + if (reactant.IsParameterized()) + fixed_reactants *= reactant.parameterize_(state.conditions_[i]); state.rate_constants_[i][(i_rate_constant++)] = - process.rate_constant_->calculate(state.conditions_[i], custom_parameters_iter); + process.rate_constant_->calculate(state.conditions_[i], custom_parameters_iter) * fixed_reactants; custom_parameters_iter += process.rate_constant_->SizeCustomParameters(); } } } template class MatrixPolicy> - requires(VectorizableDense>) void Process::UpdateState( - const std::vector& processes, - State& state) + requires(VectorizableDense>) + void Process::UpdateState(const std::vector& processes, State& state) { const auto& v_custom_parameters = state.custom_rate_parameters_.AsVector(); auto& v_rate_constants = state.rate_constants_.AsVector(); @@ -133,8 +133,13 @@ namespace micm params[i_param] = v_custom_parameters[offset_params + i_param * L + i_cell]; } std::vector::const_iterator custom_parameters_iter = params.begin(); + double fixed_reactants = 1.0; + for (auto& reactant : process.reactants_) + if (reactant.IsParameterized()) + fixed_reactants *= reactant.parameterize_(state.conditions_[i_group * L + i_cell]); v_rate_constants[offset_rc + i_cell] = - process.rate_constant_->calculate(state.conditions_[i_group * L + i_cell], custom_parameters_iter); + process.rate_constant_->calculate(state.conditions_[i_group * L + i_cell], custom_parameters_iter) * + fixed_reactants; } offset_params += params.size() * L; offset_rc += L; diff --git a/include/micm/process/process_set.hpp b/include/micm/process/process_set.hpp index 74e04eb64..39d1e5126 100644 --- a/include/micm/process/process_set.hpp +++ b/include/micm/process/process_set.hpp @@ -47,12 +47,13 @@ namespace micm /// @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( + requires(!VectorizableDense>) + void AddForcingTerms( const MatrixPolicy& rate_constants, const MatrixPolicy& state_variables, MatrixPolicy& forcing) const; template typename MatrixPolicy> - requires VectorizableDense> + requires VectorizableDense> void AddForcingTerms( const MatrixPolicy& rate_constants, const MatrixPolicy& state_variables, @@ -63,12 +64,14 @@ namespace micm /// @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> - requires(!VectorizableDense> || !VectorizableSparse>) void AddJacobianTerms( + requires(!VectorizableDense> || !VectorizableSparse>) + void AddJacobianTerms( const MatrixPolicy& rate_constants, const MatrixPolicy& state_variables, SparseMatrixPolicy& jacobian) const; template class MatrixPolicy, template class SparseMatrixPolicy> - requires(VectorizableDense>&& VectorizableSparse>) void AddJacobianTerms( + requires(VectorizableDense> && VectorizableSparse>) + void AddJacobianTerms( const MatrixPolicy& rate_constants, const MatrixPolicy& state_variables, SparseMatrixPolicy& jacobian) const; @@ -84,17 +87,25 @@ namespace micm { for (auto& process : processes) { - number_of_reactants_.push_back(process.reactants_.size()); - number_of_products_.push_back(process.products_.size()); + std::size_t number_of_reactants = 0; + std::size_t number_of_products = 0; for (auto& reactant : process.reactants_) { + if (reactant.IsParameterized()) + continue; // Skip reactants that are parameterizations reactant_ids_.push_back(state.variable_map_.at(reactant.name_)); + ++number_of_reactants; } for (auto& product : process.products_) { + if (product.first.IsParameterized()) + continue; // Skip products that are parameterizations product_ids_.push_back(state.variable_map_.at(product.first.name_)); yields_.push_back(product.second); + ++number_of_products; } + number_of_reactants_.push_back(number_of_reactants); + number_of_products_.push_back(number_of_products); } }; @@ -147,7 +158,8 @@ namespace micm } template typename MatrixPolicy> - requires(!VectorizableDense>) inline void ProcessSet::AddForcingTerms( + requires(!VectorizableDense>) + inline void ProcessSet::AddForcingTerms( const MatrixPolicy& rate_constants, const MatrixPolicy& state_variables, MatrixPolicy& forcing) const @@ -183,7 +195,7 @@ namespace micm }; template typename MatrixPolicy> - requires VectorizableDense> + requires VectorizableDense> inline void ProcessSet::AddForcingTerms( const MatrixPolicy& rate_constants, const MatrixPolicy& state_variables, @@ -224,12 +236,11 @@ namespace micm } template class MatrixPolicy, template class SparseMatrixPolicy> - requires( - !VectorizableDense> || !VectorizableSparse>) inline void ProcessSet:: - AddJacobianTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - SparseMatrixPolicy& jacobian) const + requires(!VectorizableDense> || !VectorizableSparse>) + inline void ProcessSet::AddJacobianTerms( + const MatrixPolicy& rate_constants, + const MatrixPolicy& state_variables, + SparseMatrixPolicy& jacobian) const { auto cell_jacobian = jacobian.AsVector().begin(); @@ -271,11 +282,11 @@ namespace micm } template class MatrixPolicy, template class SparseMatrixPolicy> - requires(VectorizableDense>&& VectorizableSparse>) inline void ProcessSet:: - AddJacobianTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - SparseMatrixPolicy& jacobian) const + requires(VectorizableDense> && VectorizableSparse>) + inline void ProcessSet::AddJacobianTerms( + const MatrixPolicy& rate_constants, + const MatrixPolicy& state_variables, + SparseMatrixPolicy& jacobian) const { const auto& v_rate_constants = rate_constants.AsVector(); const auto& v_state_variables = state_variables.AsVector(); diff --git a/include/micm/solver/rosenbrock.hpp b/include/micm/solver/rosenbrock.hpp index b692d1151..99f95af55 100644 --- a/include/micm/solver/rosenbrock.hpp +++ b/include/micm/solver/rosenbrock.hpp @@ -50,9 +50,11 @@ namespace micm double rejection_factor_decrease_{ 0.1 }; // used to decrease the step after 2 successive rejections double safety_factor_{ 0.9 }; // safety factor in new step size computation - double h_min_{ 0 }; // step size min - double h_max_{ 0.5 }; // step size max - double h_start_{ 0.005 }; // step size start + double h_min_{ 0.0 }; // step size min [s] + double h_max_{ + 0.0 + }; // step size max [s] (if zero or greater than the solver time-step, the time-step passed to the solver will be used) + double h_start_{ 0.0 }; // step size start [s] (if zero, 1.0e-6 will be used, if greater than h_max, h_max will be used) // Does the stage i require a new function evaluation (ros_NewF(i)=TRUE) // or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) @@ -173,7 +175,8 @@ namespace micm uint64_t decompositions{}; /// @brief The number of linear solvers uint64_t solves{}; - /// @brief The number of times a singular matrix is detected. For now, this will always be zero as we assume the matrix is never singular + /// @brief The number of times a singular matrix is detected. For now, this will always be zero as we assume the matrix + /// is never singular uint64_t singular{}; /// @brief The cumulative amount of time spent calculating the forcing function std::chrono::duration total_forcing_time{}; diff --git a/include/micm/solver/rosenbrock.inl b/include/micm/solver/rosenbrock.inl index 7d1cc6f9b..bdb582fb3 100644 --- a/include/micm/solver/rosenbrock.inl +++ b/include/micm/solver/rosenbrock.inl @@ -1,17 +1,20 @@ // Copyright (C) 2023 National Center for Atmospheric Research // SPDX-License-Identifier: Apache-2.0 -#define TIMED_METHOD(assigned_increment, time_it, method, ...) \ - { \ - if constexpr (time_it) { \ - auto start = std::chrono::high_resolution_clock::now(); \ - method(__VA_ARGS__); \ - auto end = std::chrono::high_resolution_clock::now(); \ - assigned_increment += std::chrono::duration_cast(end - start); \ - } else { \ - method(__VA_ARGS__); \ - } \ - } \ +#define TIMED_METHOD(assigned_increment, time_it, method, ...) \ + { \ + if constexpr (time_it) \ + { \ + auto start = std::chrono::high_resolution_clock::now(); \ + method(__VA_ARGS__); \ + auto end = std::chrono::high_resolution_clock::now(); \ + assigned_increment += std::chrono::duration_cast(end - start); \ + } \ + else \ + { \ + method(__VA_ARGS__); \ + } \ + } namespace micm { @@ -517,9 +520,9 @@ namespace micm MatrixPolicy forcing(Y.size(), Y[0].size(), 0.0); MatrixPolicy temp(Y.size(), Y[0].size(), 0.0); std::vector> K{}; - - parameters_.h_max_ = time_step; - parameters_.h_start_ = std::max(parameters_.h_min_, delta_min_); + 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_); stats_.Reset(); UpdateState(state); @@ -528,12 +531,14 @@ namespace micm K.push_back(MatrixPolicy(Y.size(), Y[0].size(), 0.0)); double present_time = 0.0; - double H = - std::min(std::max(std::abs(parameters_.h_min_), std::abs(parameters_.h_start_)), std::abs(parameters_.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_) H = delta_min_; + // TODO: the logic above this point should be moved to the constructor and should return an error + // if the parameters are invalid (e.g., h_min > h_max or h_start > h_max) + bool reject_last_h = false; bool reject_more_h = false; @@ -639,7 +644,7 @@ namespace micm stats_.accepted += 1; present_time = present_time + H; Y.AsVector().assign(Ynew.AsVector().begin(), Ynew.AsVector().end()); - Hnew = std::max(parameters_.h_min_, std::min(Hnew, parameters_.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 diff --git a/include/micm/system/phase.hpp b/include/micm/system/phase.hpp index 00ee36053..3928a5c2e 100644 --- a/include/micm/system/phase.hpp +++ b/include/micm/system/phase.hpp @@ -49,5 +49,21 @@ namespace micm species_ = other.species_; return *this; } + + /// @brief Returns the number of non-parameterized species + std::size_t StateSize() const + { + return std::count_if(species_.begin(), species_.end(), [](const Species& s) { return !s.IsParameterized(); }); + } + + /// @brief Returns a set of unique names for each non-parameterized species + std::vector UniqueNames() const + { + std::vector names{}; + for (const auto& species : species_) + if (!species.IsParameterized()) + names.push_back(species.name_); + return names; + } }; } // namespace micm \ No newline at end of file diff --git a/include/micm/system/species.hpp b/include/micm/system/species.hpp index e64adcef4..f4b4f0dc1 100644 --- a/include/micm/system/species.hpp +++ b/include/micm/system/species.hpp @@ -3,9 +3,12 @@ // SPDX-License-Identifier: Apache-2.0 #pragma once +#include #include #include +#include "conditions.hpp" + namespace micm { @@ -19,6 +22,12 @@ namespace micm /// @brief A list of properties of this species std::map properties_; + /// @brief A function that if provided will be used to parameterize + /// the concentration of this species during solving. + /// Species with this function defined will be excluded from + /// the solver state. + std::function parameterize_{ nullptr }; + /// @brief Copy assignment /// @param other species to copy Species& operator=(const Species& other); @@ -35,6 +44,12 @@ namespace micm /// @param name The name of the species /// @param properties The properties of the species Species(const std::string& name, const std::map& properties); + + /// @brief Returns whether a species is parameterized + bool IsParameterized() const; + + /// @brief Return a Species instance parameterized on air density + static Species ThirdBody(); }; inline Species& Species::operator=(const Species& other) @@ -46,13 +61,15 @@ namespace micm name_ = other.name_; properties_ = other.properties_; // This performs a shallow copy + parameterize_ = other.parameterize_; return *this; } inline Species::Species(const Species& other) : name_(other.name_), - properties_(other.properties_){}; + properties_(other.properties_), + parameterize_(other.parameterize_){}; inline Species::Species(const std::string& name) : name_(name){}; @@ -61,4 +78,15 @@ namespace micm : name_(name), properties_(properties){}; + inline bool Species::IsParameterized() const + { + return parameterize_ != nullptr; + } + + inline Species Species::ThirdBody() + { + Species third_body{ "M" }; + third_body.parameterize_ = [](const Conditions& c) { return c.air_density_; }; + return third_body; + } } // namespace micm diff --git a/include/micm/system/system.hpp b/include/micm/system/system.hpp index 843345617..c0b881b0f 100644 --- a/include/micm/system/system.hpp +++ b/include/micm/system/system.hpp @@ -92,10 +92,10 @@ namespace micm inline size_t System::StateSize() const { - size_t state_size = gas_phase_.species_.size(); + size_t state_size = gas_phase_.StateSize(); for (const auto& phase : phases_) { - state_size += phase.second.species_.size(); + state_size += phase.second.StateSize(); } return state_size; } @@ -108,17 +108,11 @@ namespace micm inline std::vector System::UniqueNames( const std::function& variables, const std::size_t i)> f) const { - std::vector names{}; - for (auto& species : gas_phase_.species_) - { - names.push_back(species.name_); - } + std::vector names = gas_phase_.UniqueNames(); for (auto& phase : phases_) { - for (auto& species : phase.second.species_) - { - names.push_back(phase.first + "." + species.name_); - } + for (auto& species_name : phase.second.UniqueNames()) + names.push_back(phase.first + "." + species_name); } if (f) { diff --git a/test/integration/CMakeLists.txt b/test/integration/CMakeLists.txt index 463e70d86..df835627c 100644 --- a/test/integration/CMakeLists.txt +++ b/test/integration/CMakeLists.txt @@ -8,7 +8,9 @@ include(test_util) create_standard_test(NAME chapman_integration SOURCES chapman.cpp) create_standard_test(NAME analytical_rosenbrock_integration SOURCES analytical_rosenbrock.cpp) +create_standard_test(NAME terminator_integration SOURCES terminator.cpp) if(ENABLE_LLVM) create_standard_test(NAME analytical_jit_rosenbrock_integration SOURCES analytical_jit_rosenbrock.cpp) + create_standard_test(NAME jit_terminator_integration SOURCES jit_terminator.cpp) endif() \ No newline at end of file diff --git a/test/integration/jit_terminator.cpp b/test/integration/jit_terminator.cpp new file mode 100644 index 000000000..1499054fc --- /dev/null +++ b/test/integration/jit_terminator.cpp @@ -0,0 +1,69 @@ +#include + +#include +#include +#include +#include +#include +#include + +#include "terminator.hpp" + +template class MatrixPolicy, template class SparseMatrixPolicy> +void RunTerminatorTest() +{ + TestTerminator< + MatrixPolicy, + micm::JitRosenbrockSolver< + MatrixPolicy, + SparseMatrixPolicy, + micm::JitLinearSolver>>( + [&](const micm::System& s, const std::vector& p) + -> micm::JitRosenbrockSolver< + MatrixPolicy, + SparseMatrixPolicy, + micm::JitLinearSolver> + { + auto jit{ micm::JitCompiler::create() }; + if (auto err = jit.takeError()) + { + llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]"); + EXPECT_TRUE(false); + } + auto solver_params = micm::RosenbrockSolverParameters::three_stage_rosenbrock_parameters(number_of_grid_cells, true); + solver_params.absolute_tolerance_ = 1.0e-20; + solver_params.relative_tolerance_ = 1.0e-8; + solver_params.max_number_of_steps_ = 100000; + return micm::JitRosenbrockSolver< + MatrixPolicy, + SparseMatrixPolicy, + micm::JitLinearSolver>{ jit.get(), s, p, solver_params }; + }, + number_of_grid_cells); +} + +template +using Group1VectorMatrix = micm::VectorMatrix; +template +using Group2VectorMatrix = micm::VectorMatrix; +template +using Group3VectorMatrix = micm::VectorMatrix; +template +using Group4VectorMatrix = micm::VectorMatrix; + +template +using Group1SparseVectorMatrix = micm::SparseMatrix>; +template +using Group2SparseVectorMatrix = micm::SparseMatrix>; +template +using Group3SparseVectorMatrix = micm::SparseMatrix>; +template +using Group4SparseVectorMatrix = micm::SparseMatrix>; + +TEST(JitRosenbrockSolver, Terminator) +{ + RunTerminatorTest<1, Group1VectorMatrix, Group1SparseVectorMatrix>(); + RunTerminatorTest<2, Group2VectorMatrix, Group2SparseVectorMatrix>(); + RunTerminatorTest<3, Group3VectorMatrix, Group3SparseVectorMatrix>(); + RunTerminatorTest<4, Group4VectorMatrix, Group4SparseVectorMatrix>(); +} \ No newline at end of file diff --git a/test/integration/terminator.cpp b/test/integration/terminator.cpp new file mode 100644 index 000000000..813a37c6c --- /dev/null +++ b/test/integration/terminator.cpp @@ -0,0 +1,61 @@ +#include "terminator.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include + +template class MatrixPolicy, template class SparseMatrixPolicy, class LinearSolverPolicy> +void RunTerminatorTest(std::size_t number_of_grid_cells) +{ + TestTerminator>( + [&](const micm::System& s, const std::vector& p) + -> micm::RosenbrockSolver + { + auto solver_params = micm::RosenbrockSolverParameters::three_stage_rosenbrock_parameters(number_of_grid_cells, true); + solver_params.absolute_tolerance_ = 1.0e-20; + solver_params.relative_tolerance_ = 1.0e-8; + solver_params.max_number_of_steps_ = 100000; + return micm::RosenbrockSolver{ s, p, solver_params }; + }, + number_of_grid_cells); +} + +TEST(RosenbrockSolver, Terminator) +{ + RunTerminatorTest>(2); + RunTerminatorTest>(2); + RunTerminatorTest>(3); + RunTerminatorTest>(4); +} + +template +using Group1VectorMatrix = micm::VectorMatrix; +template +using Group2VectorMatrix = micm::VectorMatrix; +template +using Group3VectorMatrix = micm::VectorMatrix; +template +using Group4VectorMatrix = micm::VectorMatrix; + +template +using Group1SparseVectorMatrix = micm::SparseMatrix>; +template +using Group2SparseVectorMatrix = micm::SparseMatrix>; +template +using Group3SparseVectorMatrix = micm::SparseMatrix>; +template +using Group4SparseVectorMatrix = micm::SparseMatrix>; + +TEST(RosenbrockSolver, VectorTerminator) +{ + RunTerminatorTest>(1); + RunTerminatorTest>(4); + RunTerminatorTest>(3); + RunTerminatorTest>(2); +} \ No newline at end of file diff --git a/test/integration/terminator.hpp b/test/integration/terminator.hpp new file mode 100644 index 000000000..1a57f475f --- /dev/null +++ b/test/integration/terminator.hpp @@ -0,0 +1,101 @@ +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/// @brief A test of the "Terminator" mechanism: +/// +/// Cl2 --hv--> 2 Cl +/// Cl + Cl --> Cl2 +/// +/// More details including analytical solution can be found here: +/// https://github.com/ESCOMP/CAM/blob/8cd44c50fe107c0b93ccd48b61eaa3d10a5b4e2f/src/chemistry/pp_terminator/chemistry.F90#L1-L434 +template class MatrixPolicy, class OdeSolverPolicy> +void TestTerminator( + const std::function&)> create_solver, + std::size_t number_of_grid_cells) +{ + auto cl2 = micm::Species("Cl2"); + auto cl = micm::Species("Cl"); + + micm::Phase gas_phase{ std::vector{ cl2, cl } }; + + micm::Process toy_r1 = micm::Process::create() + .reactants({ cl2 }) + .products({ micm::Yield(cl, 2.0) }) + .phase(gas_phase) + .rate_constant(micm::UserDefinedRateConstant({ .label_ = "toy_k1" })); + + constexpr double k2 = 1.0; + micm::Process toy_r2 = micm::Process::create() + .reactants({ cl, cl }) + .products({ micm::Yield(cl2, 1.0) }) + .phase(gas_phase) + .rate_constant(micm::ArrheniusRateConstant({ .A_ = k2 })); + + auto solver = create_solver( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ toy_r1, toy_r2 }); + micm::State state = solver.GetState(); + + auto get_double = std::bind(std::lognormal_distribution(-2.0, 2.0), std::default_random_engine()); + std::unordered_map> concentrations{ { "Cl2", {} }, { "Cl", {} } }; + for (std::size_t i_cell = 0; i_cell < number_of_grid_cells; ++i_cell) + { + concentrations["Cl2"].push_back(get_double() * 1.0e-6); + concentrations["Cl"].push_back(get_double() * 1.0e-10); + } + state.SetConcentrations(concentrations); + + std::unordered_map> custom_rate_constants{ + { "toy_k1", std::vector(number_of_grid_cells) } + }; + for (double lon = 0.0; lon < 2.0 * M_PI; lon += 0.3) + { + for (std::size_t i_cell = 0; i_cell < number_of_grid_cells; ++i_cell) + { + constexpr double k1_lat_center = M_PI * 20.0 / 180.0; + constexpr double k1_lon_center = M_PI * 300.0 / 180.0; + double lat = M_PI / 180.0 * (i_cell * (90.0 / number_of_grid_cells)); + + double k1 = std::max( + 0.0, + std::sin(lat) * std::sin(k1_lat_center) + std::cos(lat) * std::cos(k1_lat_center) * std::cos(lon - k1_lon_center)); + custom_rate_constants["toy_k1"][i_cell] = k1; + state.conditions_[i_cell].temperature_ = 298.0; // K + state.conditions_[i_cell].pressure_ = 101300.0; // Pa + state.conditions_[i_cell].air_density_ = 42.0; // mol m-3 + } + state.SetCustomRateParameters(custom_rate_constants); + + double dt = 30.0; + auto result = solver.Solve(dt, state); + + EXPECT_EQ(result.state_, micm::SolverState::Converged); + + for (std::size_t i_cell = 0; i_cell < number_of_grid_cells; ++i_cell) + { + double r = custom_rate_constants["toy_k1"][i_cell] / (4.0 * k2); + double cl_i = concentrations["Cl"][i_cell]; + double cl2_i = concentrations["Cl2"][i_cell]; + double cly = cl_i + 2.0 * cl2_i; + double det = std::sqrt(r * r + 2.0 * r * cly); + double e = std::exp(-4.0 * k2 * det * dt); + double l = (det * k2 * dt) > 1.0e-16 ? (1.0 - e) / det / dt : 4.0 * k2; + double cl_f = -l * (cl_i - det + r) * (cl_i + det + r) / (1.0 + e + dt * l * (cl_i + r)); + double cl2_f = -cl_f / 2.0; + EXPECT_NEAR( + result.result_[i_cell][state.variable_map_["Cl"]], cl_i + dt * cl_f, (cl_i + dt * cl_f) * 1.0e-8 + 1.0e-15); + EXPECT_NEAR( + result.result_[i_cell][state.variable_map_["Cl2"]], cl2_i + dt * cl2_f, (cl2_i + dt * cl2_f) * 1.0e-8 + 1.0e-15); + } + } +} diff --git a/test/unit/process/test_process.cpp b/test/unit/process/test_process.cpp index 8a1b88a45..8e29b1752 100644 --- a/test/unit/process/test_process.cpp +++ b/test/unit/process/test_process.cpp @@ -14,11 +14,14 @@ template class MatrixPolicy> 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 } }); + micm::Species bar("bar"); + bar.parameterize_ = [](const micm::Conditions& c) { return c.air_density_ * 0.82; }; + micm::ArrheniusRateConstant rc1({ .A_ = 12.2, .C_ = 300.0 }); micm::SurfaceRateConstant rc2({ .label_ = "foo_surf", .species_ = foo }); micm::UserDefinedRateConstant rc3({ .label_ = "bar_user" }); - micm::Process r1 = micm::Process::create().rate_constant(rc1); + micm::Process r1 = micm::Process::create().rate_constant(rc1).reactants({ foo, bar }); micm::Process r2 = micm::Process::create().rate_constant(rc2); micm::Process r3 = micm::Process::create().rate_constant(rc3); std::vector processes = { r1, r2, r3 }; @@ -39,6 +42,7 @@ void testProcessUpdateState(const std::size_t number_of_grid_cells) { state.conditions_[i_cell].temperature_ = get_double() * 285.0; state.conditions_[i_cell].pressure_ = get_double() * 101100.0; + state.conditions_[i_cell].air_density_ = get_double() * 10.0; params[0] = get_double() * 1.0e-8; params[1] = get_double() * 1.0e5; params[2] = get_double() * 1.0e-2; @@ -48,7 +52,8 @@ void testProcessUpdateState(const std::size_t number_of_grid_cells) params[1]; state.custom_rate_parameters_[i_cell][state.custom_rate_parameter_map_["bar_user"]] = params[2]; std::vector::const_iterator param_iter = params.begin(); - expected_rate_constants[i_cell][0] = rc1.calculate(state.conditions_[i_cell], param_iter); + expected_rate_constants[i_cell][0] = + rc1.calculate(state.conditions_[i_cell], param_iter) * (state.conditions_[i_cell].air_density_ * 0.82); param_iter += rc1.SizeCustomParameters(); expected_rate_constants[i_cell][1] = rc2.calculate(state.conditions_[i_cell], param_iter); param_iter += rc2.SizeCustomParameters(); diff --git a/test/unit/process/test_process_set_policy.hpp b/test/unit/process/test_process_set_policy.hpp index 1e050b006..38cb54ed4 100644 --- a/test/unit/process/test_process_set_policy.hpp +++ b/test/unit/process/test_process_set_policy.hpp @@ -21,8 +21,10 @@ void testProcessSet( auto baz = micm::Species("baz"); auto quz = micm::Species("quz"); auto quuz = micm::Species("quuz"); + auto qux = micm::Species("qux"); + qux.parameterize_ = [](const micm::Conditions& c) { return c.air_density_ * 0.72; }; - micm::Phase gas_phase{ std::vector{ foo, bar, baz, quz, quuz } }; + micm::Phase gas_phase{ std::vector{ foo, bar, qux, baz, quz, quuz } }; micm::State state{ micm::StateParameters{ .state_variable_names_{ "foo", "bar", "baz", "quz", "quuz" }, .number_of_grid_cells_ = 2, @@ -32,7 +34,7 @@ void testProcessSet( micm::Process::create().reactants({ foo, baz }).products({ yields(bar, 1), yields(quuz, 2.4) }).phase(gas_phase); micm::Process r2 = - micm::Process::create().reactants({ bar }).products({ yields(foo, 1), yields(quz, 1.4) }).phase(gas_phase); + micm::Process::create().reactants({ bar, qux }).products({ yields(foo, 1), yields(quz, 1.4) }).phase(gas_phase); micm::Process r3 = micm::Process::create().reactants({ quz }).products({}).phase(gas_phase); diff --git a/test/unit/system/test_phase.cpp b/test/unit/system/test_phase.cpp index 284433129..0565af37d 100644 --- a/test/unit/system/test_phase.cpp +++ b/test/unit/system/test_phase.cpp @@ -8,4 +8,26 @@ TEST(Phase, ConstructorWithVector) micm::Phase phase(std::vector({ micm::Species("species1"), micm::Species("species2") })); EXPECT_EQ(phase.species_.size(), 2); + EXPECT_EQ(phase.StateSize(), 2); + + auto names = phase.UniqueNames(); + EXPECT_EQ(names[0], "species1"); + EXPECT_EQ(names[1], "species2"); +} + +TEST(Phase, ConstructorWithParameterizedSpecies) +{ + auto foo = micm::Species("foo"); + auto bar = micm::Species("bar"); + auto baz = micm::Species("baz"); + + bar.parameterize_ = [](const micm::Conditions& c) { return 42.0; }; + micm::Phase phase(std::vector({ foo, bar, baz })); + + EXPECT_EQ(phase.species_.size(), 3); + EXPECT_EQ(phase.StateSize(), 2); + + auto names = phase.UniqueNames(); + EXPECT_EQ(names[0], "foo"); + EXPECT_EQ(names[1], "baz"); } \ No newline at end of file diff --git a/test/unit/system/test_species.cpp b/test/unit/system/test_species.cpp index a4f4fcc51..dcfab1f89 100644 --- a/test/unit/system/test_species.cpp +++ b/test/unit/system/test_species.cpp @@ -5,8 +5,11 @@ TEST(Species, StringConstructor) { micm::Species species("thing"); - EXPECT_EQ(species.name_, "thing"); + EXPECT_FALSE(species.IsParameterized()); + species.parameterize_ = [](const micm::Conditions& c) { return c.temperature_ * 100.0; }; + EXPECT_TRUE(species.IsParameterized()); + EXPECT_EQ(species.parameterize_({ .temperature_ = 12.5 }), 12.5 * 100.0); } TEST(Species, StringAndVectorConstructor) @@ -19,26 +22,66 @@ TEST(Species, StringAndVectorConstructor) EXPECT_EQ(species.properties_["name2 [units2]"], 2.0); } +TEST(Species, ThirdBody) +{ + micm::Species species = micm::Species::ThirdBody(); + EXPECT_EQ(species.name_, "M"); + EXPECT_TRUE(species.IsParameterized()); + EXPECT_EQ(species.parameterize_({ .air_density_ = 42.4 }), 42.4); +} + TEST(Species, CopyConstructor) { - micm::Species species("thing", { { "name [units]", 1.0 }, { "name2 [units2]", 2.0 } }); + { + micm::Species species("thing", { { "name [units]", 1.0 }, { "name2 [units2]", 2.0 } }); - micm::Species species2(species); + micm::Species species2(species); - EXPECT_EQ(species2.name_, "thing"); - EXPECT_EQ(species2.properties_.size(), 2); - EXPECT_EQ(species2.properties_["name [units]"], 1.0); - EXPECT_EQ(species2.properties_["name2 [units2]"], 2.0); + EXPECT_EQ(species2.name_, "thing"); + EXPECT_EQ(species2.properties_.size(), 2); + EXPECT_EQ(species2.properties_["name [units]"], 1.0); + EXPECT_EQ(species2.properties_["name2 [units2]"], 2.0); + EXPECT_FALSE(species2.IsParameterized()); + } + { + micm::Species species("thing", { { "name [units]", 1.0 }, { "name2 [units2]", 2.0 } }); + species.parameterize_ = [](const micm::Conditions& c) { return 15.4; }; + + micm::Species species2(species); + + EXPECT_EQ(species2.name_, "thing"); + EXPECT_EQ(species2.properties_.size(), 2); + EXPECT_EQ(species2.properties_["name [units]"], 1.0); + EXPECT_EQ(species2.properties_["name2 [units2]"], 2.0); + EXPECT_TRUE(species2.IsParameterized()); + EXPECT_EQ(species.parameterize_({}), 15.4); + } } TEST(Species, CopyAssignment) { - micm::Species species("thing", { { "name [units]", 1.0 }, { "name2 [units2]", 2.0 } }); + { + micm::Species species("thing", { { "name [units]", 1.0 }, { "name2 [units2]", 2.0 } }); + + micm::Species species2 = species; + + EXPECT_EQ(species2.name_, "thing"); + EXPECT_EQ(species2.properties_.size(), 2); + EXPECT_EQ(species2.properties_["name [units]"], 1.0); + EXPECT_EQ(species2.properties_["name2 [units2]"], 2.0); + EXPECT_FALSE(species2.IsParameterized()); + } + { + micm::Species species("thing", { { "name [units]", 1.0 }, { "name2 [units2]", 2.0 } }); + species.parameterize_ = [](const micm::Conditions& c) { return 15.4; }; - micm::Species species2 = species; + micm::Species species2 = species; - EXPECT_EQ(species2.name_, "thing"); - EXPECT_EQ(species2.properties_.size(), 2); - EXPECT_EQ(species2.properties_["name [units]"], 1.0); - EXPECT_EQ(species2.properties_["name2 [units2]"], 2.0); + EXPECT_EQ(species2.name_, "thing"); + EXPECT_EQ(species2.properties_.size(), 2); + EXPECT_EQ(species2.properties_["name [units]"], 1.0); + EXPECT_EQ(species2.properties_["name2 [units2]"], 2.0); + EXPECT_TRUE(species2.IsParameterized()); + EXPECT_EQ(species.parameterize_({}), 15.4); + } } diff --git a/test/unit/system/test_system.cpp b/test/unit/system/test_system.cpp index 750a17b38..ec60f3a20 100644 --- a/test/unit/system/test_system.cpp +++ b/test/unit/system/test_system.cpp @@ -27,6 +27,45 @@ TEST(System, ConstructorWithAllParameters) EXPECT_NE(std::find(names.begin(), names.end(), "phase2.species3"), names.end()); EXPECT_NE(std::find(names.begin(), names.end(), "phase2.species4"), names.end()); + std::vector reorder{ 3, 2, 1, 0, 5, 4 }; + auto reordered_names = system.UniqueNames([&](const std::vector variables, const std::size_t i) + { return variables[reorder[i]]; }); + EXPECT_EQ(reordered_names.size(), 6); + EXPECT_EQ(reordered_names[0], names[3]); + EXPECT_EQ(reordered_names[1], names[2]); + EXPECT_EQ(reordered_names[2], names[1]); + EXPECT_EQ(reordered_names[3], names[0]); + EXPECT_EQ(reordered_names[4], names[5]); + EXPECT_EQ(reordered_names[5], names[4]); +} + +TEST(System, ConstructorWithParameterizedSpecies) +{ + std::vector speciesA = { micm::Species("species1"), micm::Species("species2") }; + std::vector speciesB = { micm::Species("species3"), micm::Species("species4") }; + auto param_species = micm::Species("paramSpecies"); + param_species.parameterize_ = [](const micm::Conditions& c) { return 64.2; }; + speciesA.push_back(param_species); + + micm::Phase phase = speciesA; + std::unordered_map phases = { { "phase1", speciesA }, { "phase2", speciesB } }; + + micm::System system = { micm::SystemParameters{ .gas_phase_ = phase, .phases_ = phases } }; + + EXPECT_EQ(system.gas_phase_.species_.size(), 3); + EXPECT_EQ(system.phases_.size(), 2); + EXPECT_EQ(system.StateSize(), 6); + + auto names = system.UniqueNames(); + + EXPECT_EQ(names.size(), 6); + EXPECT_NE(std::find(names.begin(), names.end(), "species1"), names.end()); + EXPECT_NE(std::find(names.begin(), names.end(), "species2"), names.end()); + EXPECT_NE(std::find(names.begin(), names.end(), "phase1.species1"), names.end()); + EXPECT_NE(std::find(names.begin(), names.end(), "phase1.species2"), names.end()); + EXPECT_NE(std::find(names.begin(), names.end(), "phase2.species3"), names.end()); + EXPECT_NE(std::find(names.begin(), names.end(), "phase2.species4"), names.end()); + std::vector reorder{ 3, 2, 1, 0, 5, 4 }; auto reordered_names = system.UniqueNames([&](const std::vector variables, const std::size_t i) { return variables[reorder[i]]; });