From 78b2985cd4d8cc62b0eebbaf854d54d92cbb598f Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sun, 13 Mar 2022 20:42:09 -0400 Subject: [PATCH 1/3] Remove unused return value from MultiRate::replace --- include/cantera/kinetics/MultiRate.h | 4 +--- include/cantera/kinetics/MultiRateBase.h | 2 +- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/include/cantera/kinetics/MultiRate.h b/include/cantera/kinetics/MultiRate.h index d8e73e5ba36..811c84c5bf5 100644 --- a/include/cantera/kinetics/MultiRate.h +++ b/include/cantera/kinetics/MultiRate.h @@ -39,7 +39,7 @@ class MultiRate final : public MultiRateBase m_shared.invalidateCache(); } - virtual bool replace(const size_t rxn_index, ReactionRate& rate) override { + virtual void replace(const size_t rxn_index, ReactionRate& rate) override { if (!m_rxn_rates.size()) { throw CanteraError("MultiRate::replace", "Invalid operation: cannot replace rate object " @@ -54,9 +54,7 @@ class MultiRate final : public MultiRateBase if (m_indices.find(rxn_index) != m_indices.end()) { size_t j = m_indices[rxn_index]; m_rxn_rates.at(j).second = dynamic_cast(rate); - return true; } - return false; } virtual void resize(size_t n_species, size_t n_reactions) override { diff --git a/include/cantera/kinetics/MultiRateBase.h b/include/cantera/kinetics/MultiRateBase.h index 7df318626c0..0c182d1effc 100644 --- a/include/cantera/kinetics/MultiRateBase.h +++ b/include/cantera/kinetics/MultiRateBase.h @@ -39,7 +39,7 @@ class MultiRateBase //! Replace reaction rate object handled by the evaluator //! @param rxn_index index of reaction //! @param rate reaction rate object - virtual bool replace(size_t rxn_index, ReactionRate& rate) = 0; + virtual void replace(size_t rxn_index, ReactionRate& rate) = 0; //! Update number of species and reactions //! @param n_species number of species From acebd5ca56bb69d8253eae17286f7d7ad752005c Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sun, 13 Mar 2022 15:00:57 -0400 Subject: [PATCH 2/3] Improve performance of ArrheniusRate evaluation - Use Eigen to vectorize rate evaluations - Only do calculations for rates that are actually temperature-dependent --- include/cantera/kinetics/Arrhenius.h | 55 ++++++++++++++++++- include/cantera/kinetics/MultiRate.h | 2 +- src/kinetics/Arrhenius.cpp | 79 ++++++++++++++++++++++++++++ src/kinetics/BulkKinetics.cpp | 3 ++ src/kinetics/InterfaceKinetics.cpp | 3 ++ 5 files changed, 140 insertions(+), 2 deletions(-) diff --git a/include/cantera/kinetics/Arrhenius.h b/include/cantera/kinetics/Arrhenius.h index 683dbe779a3..5cb8d472fcc 100644 --- a/include/cantera/kinetics/Arrhenius.h +++ b/include/cantera/kinetics/Arrhenius.h @@ -9,6 +9,8 @@ #include "cantera/kinetics/ReactionData.h" #include "ReactionRate.h" #include "MultiRate.h" +#include "cantera/numerics/eigen_dense.h" +#include "cantera/base/global.h" namespace Cantera { @@ -450,10 +452,61 @@ class BulkRate : public RateType } }; -typedef BulkRate ArrheniusRate; typedef BulkRate TwoTempPlasmaRate; typedef BulkRate BlowersMaselRate; +class ArrheniusRate : public BulkRate +{ + using BulkRate::BulkRate; + unique_ptr newMultiRate() const override; +}; + +//! A specialized rate evaluator for ArrheniusRate that avoids unnecessary rate +//! calculations for reactions where the rate is simply a constant +class MultiArrheniusRate : public MultiRate +{ +public: + virtual void resize(size_t n_species, size_t n_reactions) override; + virtual void replace(const size_t rxn_index, ReactionRate& rate) override; + virtual void getRateConstants(double* kf) override; + +protected: + //! Partition rates based on whether they are temperature dependent or not + void partitionRates(); + + //! Pre-exponential factors for reactions where the rate is temperature dependent. + //! The corresponding element of #m_variable gives the reaction index. + Eigen::ArrayXd m_A; + + //! Temperature exponents for reactions where the rate is temperature dependent + //! The corresponding element of #m_variable gives the reaction index. + Eigen::ArrayXd m_b; + + //! Activation energies for reactions where the rate is temperature dependent + //! The corresponding element of #m_variable gives the reaction index. + Eigen::ArrayXd m_Ea_R; + + //! Forward rate constants for reactions where the rate is temperature dependent + //! The corresponding element of #m_variable gives the reaction index. + Eigen::ArrayXd m_kf_var; + + //! Forward rate constants for reactions where the rate is not temperature dependent + //! The corresponding element of #m_const gives the reaction index. + Eigen::ArrayXd m_kf_const; + + //! Reaction indices for reactions where the rate is not temperature dependent + std::vector m_const; + + //! Reaction indices for reactions where the rate is temperature dependent + std::vector m_variable; + + //! Map from reaction index to index in #m_const + std::map m_const_indices; + + //! Map from reaction index to index in #m_variable + std::map m_variable_indices; +}; + } #endif diff --git a/include/cantera/kinetics/MultiRate.h b/include/cantera/kinetics/MultiRate.h index 811c84c5bf5..43d9432084b 100644 --- a/include/cantera/kinetics/MultiRate.h +++ b/include/cantera/kinetics/MultiRate.h @@ -17,7 +17,7 @@ namespace Cantera //! A class template handling ReactionRate specializations. template -class MultiRate final : public MultiRateBase +class MultiRate : public MultiRateBase { CT_DEFINE_HAS_MEMBER(has_update, updateFromStruct) CT_DEFINE_HAS_MEMBER(has_ddT, ddTScaledFromStruct) diff --git a/src/kinetics/Arrhenius.cpp b/src/kinetics/Arrhenius.cpp index 469bba1171b..3384ba8d813 100644 --- a/src/kinetics/Arrhenius.cpp +++ b/src/kinetics/Arrhenius.cpp @@ -194,4 +194,83 @@ void BlowersMasel::setContext(const Reaction& rxn, const Kinetics& kin) } } +unique_ptr ArrheniusRate::newMultiRate() const +{ + return unique_ptr(new MultiArrheniusRate()); +} + +void MultiArrheniusRate::resize(size_t n_species, size_t n_reactions) +{ + MultiRate::resize(n_species, n_reactions); + + if (m_kf_var.size() + m_kf_const.size() != m_rxn_rates.size()) { + partitionRates(); + } +} + +void MultiArrheniusRate::partitionRates() +{ + m_const.clear(); + m_variable.clear(); + m_const_indices.clear(); + m_variable_indices.clear(); + + vector_fp A_const, A_var, b, Ea_R; + for (size_t i = 0; i < m_rxn_rates.size(); i++) { + auto& rate = m_rxn_rates[i].second; + if (rate.temperatureExponent() == 0 && rate.activationEnergy() == 0) { + m_const.push_back(m_rxn_rates[i].first); + m_const_indices[m_rxn_rates[i].first] = m_const.size() - 1; + A_const.push_back(rate.preExponentialFactor()); + } else { + m_variable.push_back(m_rxn_rates[i].first); + m_variable_indices[m_rxn_rates[i].first] = m_variable.size() - 1; + A_var.push_back(rate.preExponentialFactor()); + b.push_back(rate.temperatureExponent()); + Ea_R.push_back(rate.activationEnergy() / GasConstant); + } + } + + m_kf_const = Eigen::Map(A_const.data(), A_const.size()); + m_A = Eigen::Map(A_var.data(), A_var.size()); + m_b = Eigen::Map(b.data(), b.size()); + m_Ea_R = Eigen::Map(Ea_R.data(), Ea_R.size()); + m_kf_var.resize(A_var.size()); +} + +void MultiArrheniusRate::replace(const size_t rxn_index, ReactionRate& rate) +{ + auto& old_rate = m_rxn_rates[m_indices.at(rxn_index)].second; + bool old_const = (old_rate.temperatureExponent() == 0 && + old_rate.activationEnergy() == 0); + MultiRate::replace(rxn_index, rate); + auto& new_rate = m_rxn_rates[m_indices.at(rxn_index)].second; + bool new_const = (new_rate.temperatureExponent() == 0 && + new_rate.activationEnergy() == 0); + if (old_const && new_const) { + size_t j = m_const_indices.at(rxn_index); + m_kf_const[j] = new_rate.preExponentialFactor(); + } else if (!old_const && !new_const) { + size_t j = m_variable_indices.at(rxn_index); + m_A[j] = new_rate.preExponentialFactor(); + m_b[j] = new_rate.temperatureExponent(); + m_Ea_R[j] = new_rate.activationEnergy() / GasConstant; + } else { + partitionRates(); + } +} + +void MultiArrheniusRate::getRateConstants(double* kf) +{ + AssertThrow(m_kf_var.size() + m_kf_const.size() == m_rxn_rates.size(), + "MultiArrheniusRate::getRateConstants"); + m_kf_var = m_A * (m_b * m_shared.logT - m_Ea_R * m_shared.recipT).exp(); + for (size_t i = 0; i < m_kf_var.size(); i++) { + kf[m_variable[i]] = m_kf_var[i]; + } + for (size_t i = 0; i < m_kf_const.size(); i++) { + kf[m_const[i]] = m_kf_const[i]; + } +} + } diff --git a/src/kinetics/BulkKinetics.cpp b/src/kinetics/BulkKinetics.cpp index 1136a26e889..4cb772f9a74 100644 --- a/src/kinetics/BulkKinetics.cpp +++ b/src/kinetics/BulkKinetics.cpp @@ -151,6 +151,9 @@ bool BulkKinetics::addReaction(shared_ptr r, bool resize) // Add reaction rate to evaluator size_t index = m_bulk_types[rate->type()]; m_bulk_rates[index]->add(nReactions() - 1, *rate); + if (resize) { + m_bulk_rates[index]->resize(nReactions(), nTotalSpecies()); + } // Add reaction to third-body evaluator if (r->thirdBody() != nullptr) { diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 760f8e045d7..8556460c6c0 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -548,6 +548,9 @@ bool InterfaceKinetics::addReaction(shared_ptr r_base, bool resize) // Add reaction rate to evaluator size_t index = m_interface_types[rate->type()]; m_interface_rates[index]->add(nReactions() - 1, *rate); + if (resize) { + m_interface_rates[index]->resize(nReactions(), nTotalSpecies()); + } } else if (r_base->reaction_type == BMINTERFACE_RXN) { throw NotImplementedError("InterfaceKinetics::addReaction"); From 46b3f2790ff4e21ee2e93fe00d12a6d683031e66 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 14 Mar 2022 11:38:39 -0400 Subject: [PATCH 3/3] [Test] Eliminate some exact floating point equality checks --- test/kinetics/kineticsFromScratch.cpp | 23 ++++++++++++++--------- test/kinetics/kineticsFromScratch3.cpp | 17 ++++++++++------- 2 files changed, 24 insertions(+), 16 deletions(-) diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index 3867e9d59f4..c905036439d 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -40,11 +40,11 @@ class KineticsFromScratch : public testing::Test kin.getFwdRateConstants(&k[0]); kin_ref->getFwdRateConstants(&k_ref[0]); - EXPECT_DOUBLE_EQ(k_ref[iRef], k[0]); + EXPECT_NEAR(k_ref[iRef], k[0], 1e-14 * k_ref[iRef]); kin.getRevRateConstants(&k[0]); kin_ref->getRevRateConstants(&k_ref[0]); - EXPECT_DOUBLE_EQ(k_ref[iRef], k[0]); + EXPECT_NEAR(k_ref[iRef], k[0], 1e-14 * k_ref[iRef]); } }; @@ -344,11 +344,11 @@ class InterfaceKineticsFromScratch : public testing::Test kin.getFwdRateConstants(&k[0]); kin_ref->getFwdRateConstants(&k_ref[0]); - EXPECT_DOUBLE_EQ(k_ref[iRef], k[0]); + EXPECT_NEAR(k_ref[iRef], k[0], 1e-14 * k_ref[iRef]); kin.getRevRateConstants(&k[0]); kin_ref->getRevRateConstants(&k_ref[0]); - EXPECT_DOUBLE_EQ(k_ref[iRef], k[0]); + EXPECT_NEAR(k_ref[iRef], k[0], 1e-14 * k_ref[iRef]); } }; @@ -433,32 +433,37 @@ class KineticsAddSpecies : public testing::Test kin.getFwdRateConstants(k.data()); kin_ref->getFwdRateConstants(k_ref.data()); for (size_t i = 0; i < kin.nReactions(); i++) { - EXPECT_DOUBLE_EQ(k_ref[i], k[i]) << "i = " << i << "; N = " << N; + EXPECT_NEAR(k_ref[i], k[i], 1e-14 * k_ref[i]) + << "i = " << i << "; N = " << N; } kin.getFwdRatesOfProgress(k.data()); kin_ref->getFwdRatesOfProgress(k_ref.data()); for (size_t i = 0; i < kin.nReactions(); i++) { - EXPECT_DOUBLE_EQ(k_ref[i], k[i]) << "i = " << i << "; N = " << N; + EXPECT_NEAR(k_ref[i], k[i], 1e-14 * k_ref[i]) + << "i = " << i << "; N = " << N; } kin.getRevRateConstants(k.data()); kin_ref->getRevRateConstants(k_ref.data()); for (size_t i = 0; i < kin.nReactions(); i++) { - EXPECT_DOUBLE_EQ(k_ref[i], k[i]) << "i = " << i << "; N = " << N; + EXPECT_NEAR(k_ref[i], k[i], 1e-14 * k_ref[i]) + << "i = " << i << "; N = " << N; } kin.getRevRatesOfProgress(k.data()); kin_ref->getRevRatesOfProgress(k_ref.data()); for (size_t i = 0; i < kin.nReactions(); i++) { - EXPECT_DOUBLE_EQ(k_ref[i], k[i]) << "i = " << i << "; N = " << N; + EXPECT_NEAR(k_ref[i], k[i], 1e-14 * k_ref[i]) + << "i = " << i << "; N = " << N; } kin.getCreationRates(w.data()); kin_ref->getCreationRates(w_ref.data()); for (size_t i = 0; i < kin.nTotalSpecies(); i++) { size_t iref = p_ref.speciesIndex(p.speciesName(i)); - EXPECT_DOUBLE_EQ(w_ref[iref], w[i]) << "sp = " << p.speciesName(i) << "; N = " << N; + EXPECT_NEAR(w_ref[iref], w[i], 1e-14 * w_ref[iref]) + << "sp = " << p.speciesName(i) << "; N = " << N; } } }; diff --git a/test/kinetics/kineticsFromScratch3.cpp b/test/kinetics/kineticsFromScratch3.cpp index 7c912d564f1..089640afe89 100644 --- a/test/kinetics/kineticsFromScratch3.cpp +++ b/test/kinetics/kineticsFromScratch3.cpp @@ -42,11 +42,11 @@ class KineticsFromScratch3 : public testing::Test kin.getFwdRateConstants(&k[0]); kin_ref->getFwdRateConstants(&k_ref[0]); - EXPECT_DOUBLE_EQ(k_ref[iRef], k[0]); + EXPECT_NEAR(k_ref[iRef], k[0], 1e-14 * k_ref[iRef]); kin.getRevRateConstants(&k[0]); kin_ref->getRevRateConstants(&k_ref[0]); - EXPECT_DOUBLE_EQ(k_ref[iRef], k[0]); + EXPECT_NEAR(k_ref[iRef], k[0], 1e-14 * k_ref[iRef]); } }; @@ -354,11 +354,11 @@ class InterfaceKineticsFromScratch3 : public testing::Test kin.getFwdRateConstants(&k[0]); kin_ref->getFwdRateConstants(&k_ref[0]); - EXPECT_DOUBLE_EQ(k_ref[iRef], k[0]); + EXPECT_NEAR(k_ref[iRef], k[0], 1e-14 * k_ref[iRef]); kin.getRevRateConstants(&k[0]); kin_ref->getRevRateConstants(&k_ref[0]); - EXPECT_DOUBLE_EQ(k_ref[iRef], k[0]); + EXPECT_NEAR(k_ref[iRef], k[0], 1e-14 * k_ref[iRef]); } }; @@ -447,13 +447,15 @@ class KineticsAddSpecies3 : public testing::Test kin.getFwdRateConstants(k.data()); kin_ref->getFwdRateConstants(k_ref.data()); for (size_t i = 0; i < kin.nReactions(); i++) { - EXPECT_DOUBLE_EQ(k_ref[i], k[i]) << "i = " << i << "; N = " << N; + EXPECT_NEAR(k_ref[i], k[i], 1e-14 * k_ref[i]) + << "i = " << i << "; N = " << N; } kin.getFwdRatesOfProgress(k.data()); kin_ref->getFwdRatesOfProgress(k_ref.data()); for (size_t i = 0; i < kin.nReactions(); i++) { - EXPECT_DOUBLE_EQ(k_ref[i], k[i]) << "i = " << i << "; N = " << N; + EXPECT_NEAR(k_ref[i], k[i], 1e-14 * k_ref[i]) + << "i = " << i << "; N = " << N; } kin.getRevRateConstants(k.data()); @@ -472,7 +474,8 @@ class KineticsAddSpecies3 : public testing::Test kin_ref->getCreationRates(w_ref.data()); for (size_t i = 0; i < kin.nTotalSpecies(); i++) { size_t iref = p_ref.speciesIndex(p.speciesName(i)); - EXPECT_NEAR(w_ref[iref], w[i], w_ref[iref]*1e-12) << "sp = " << p.speciesName(i) << "; N = " << N; + EXPECT_NEAR(w_ref[iref], w[i], w_ref[iref]*1e-12) + << "sp = " << p.speciesName(i) << "; N = " << N; } } };