Improve performance of Arrhenius rate evaluations

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<Arrhenius3, ArrheniusData> ArrheniusRate;
 typedef BulkRate<TwoTempPlasma, TwoTempPlasmaData> TwoTempPlasmaRate;
 typedef BulkRate<BlowersMasel, BlowersMaselData> BlowersMaselRate;
 
+class ArrheniusRate : public BulkRate<Arrhenius3, ArrheniusData>
+{
+    using BulkRate::BulkRate;
+    unique_ptr<MultiRateBase> 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<ArrheniusRate, ArrheniusData>
+{
+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<size_t> m_const;
+
+    //! Reaction indices for reactions where the rate is temperature dependent
+    std::vector<size_t> m_variable;
+
+    //! Map from reaction index to index in #m_const
+    std::map<size_t, size_t> m_const_indices;
+
+    //! Map from reaction index to index in #m_variable
+    std::map<size_t, size_t> m_variable_indices;
+};
+
 }
 
 #endif

include/cantera/kinetics/MultiRate.h

@@ -17,7 +17,7 @@ namespace Cantera
 
 //! A class template handling ReactionRate specializations.
 template <class RateType, class DataType>
-class MultiRate final : public MultiRateBase
+class MultiRate : public MultiRateBase
 {
     CT_DEFINE_HAS_MEMBER(has_update, updateFromStruct)
     CT_DEFINE_HAS_MEMBER(has_ddT, ddTScaledFromStruct)
@@ -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<RateType&>(rate);
-            return true;
         }
-        return false;
     }
 
     virtual void resize(size_t n_species, size_t n_reactions) override {

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

src/kinetics/Arrhenius.cpp

@@ -194,4 +194,83 @@ void BlowersMasel::setContext(const Reaction& rxn, const Kinetics& kin)
     }
 }
 
+unique_ptr<MultiRateBase> ArrheniusRate::newMultiRate() const
+{
+    return unique_ptr<MultiRateBase>(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<Eigen::ArrayXd>(A_const.data(), A_const.size());
+    m_A = Eigen::Map<Eigen::ArrayXd>(A_var.data(), A_var.size());
+    m_b = Eigen::Map<Eigen::ArrayXd>(b.data(), b.size());
+    m_Ea_R = Eigen::Map<Eigen::ArrayXd>(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];
+    }
+}
+
 }

src/kinetics/BulkKinetics.cpp

@@ -151,6 +151,9 @@ bool BulkKinetics::addReaction(shared_ptr<Reaction> 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) {

src/kinetics/InterfaceKinetics.cpp

@@ -548,6 +548,9 @@ bool InterfaceKinetics::addReaction(shared_ptr<Reaction> 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");

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;
         }
     }
 };

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;
         }
     }
 };