Initial commit for polynomial Ea coverage dependency

Add documentation

doc/sphinx/yaml/reactions.rst

@@ -388,7 +388,9 @@ Includes the fields of an :ref:`sec-yaml-elementary` reaction plus:
 
     ``E``
         Activation energy dependence on coverage, which uses the same sign convention
-        as the leading-order activation energy term
+        as the leading-order activation energy term. This can be a scalar value for
+        the linear dependency or a list of four values for the polynomial dependency
+        given in the order of 1st, 2nd, 3rd, and 4th-order coefficients
 
     or a list containing the three elements above, in the given order.
 
@@ -403,12 +405,21 @@ Examples::
       rate-constant: {A: 3.7e21 cm^2/mol/s, b: 0, Ea: 67400 J/mol}
       coverage-dependencies: {H(s): {a: 0, m: 0, E: -6000 J/mol}}
 
+    - equation: 2 O(S) => O2 + 2 Pt(S)
+      rate-constant: {A: 3.7e+21, b: 0, Ea: 213200 J/mol}
+      coverage-dependencies: {O(S): {a: 0.0, m: 0.0,
+        E: [1.0e3 J/mol, 3.0e3 J/mol , -7.0e4 J/mol , 5.0e3 J/mol]}
+
     - equation: CH4 + PT(S) + O(S) => CH3(S) + OH(S)
       rate-constant: {A: 5.0e+18, b: 0.7, Ea: 4.2e+04}
       coverage-dependencies:
         O(S): [0, 0, 8000]
         PT(S): [0, -1.0, 0]
 
+    - equation: 2 O(S) => O2 + 2 Pt(S)
+      rate-constant: {A: 3.7e+21, b: 0, Ea: 213200 J/mol}
+      coverage-dependencies:
+        O(S): [0, 0, [1.0e6, 3.0e6, -7.0e7, 5.0e6]]
 
 .. _sec-yaml-interface-Blowers-Masel:
 

include/cantera/kinetics/InterfaceRate.h

@@ -84,9 +84,20 @@ struct InterfaceData : public BlowersMaselData
  * It is evident that this expression combines a regular modified Arrhenius rate
  * expression \f$ A T^b \exp \left( - \frac{E_a}{RT} \right) \f$ with coverage-related
  * terms, where the parameters \f$ (a_k, E_k, m_k) \f$ describe the dependency on the
- * surface coverage of species \f$ k, \theta_k \f$. The InterfaceRateBase class implements
- * terms related to coverage only, which allows for combinations with arbitrary rate
- * parameterizations (for example Arrhenius and BlowersMaselRate).
+ * surface coverage of species \f$ k, \theta_k \f$. In addition to the linear coverage
+ * dependence on the activation energy modifier \f$ E_k \f$, polynomial coverage
+ * dependence is also available. When the dependence parameter \f$ E_k \f$ is given as
+ * a scalar value, the linear dependency is applied whereas if a list of four values
+ * are given as \f$ [E^{(1)}_k-E^{(4)}_k] \f$, a polynomial dependency is applied as
+ *  \f[
+ *      k_f = A T^b \exp \left( - \frac{E_a}{RT} \right)
+ *          \prod_k 10^{a_k \theta_k} \theta_k^{m_k}
+ *          \exp \left( \frac{- E^{(1)}_k \theta_k - E^{(2)}_k \theta_k^2
+ *          - E^{(3)}_k \theta_k^3 - E^{(4)}_k \theta_k^4}{RT} \right)
+ *  \f]
+ * The InterfaceRateBase class implements terms related to coverage only, which allows
+ * for combinations with arbitrary rate parameterizations (for example Arrhenius and
+ * BlowersMaselRate).
  */
 class InterfaceRateBase
 {
@@ -116,7 +127,7 @@ class InterfaceRateBase
     //! Add a coverage dependency for species *sp*, with exponential dependence
     //! *a*, power-law exponent *m*, and activation energy dependence *e*,
     //! where *e* is in Kelvin, that is, energy divided by the molar gas constant.
-    virtual void addCoverageDependence(const string& sp, double a, double m, double e);
+    virtual void addCoverageDependence(const string& sp, double a, double m, vector_fp e);
 
     //! Boolean indicating whether rate uses exchange current density formulation
     bool exchangeCurrentDensityFormulation() {
@@ -243,7 +254,11 @@ class InterfaceRateBase
     map<size_t, size_t> m_indices;
     vector<string> m_cov; //!< Vector holding names of coverage species
     vector_fp m_ac; //!< Vector holding coverage-specific exponential dependence
-    vector_fp m_ec; //!< Vector holding coverage-specific activation energy dependence
+    //! Vector holding coverage-specific activation energy dependence as a
+    //! 5-membered array of polynomial coeffcients starting from 0th-order to
+    //! 4th-order coefficients
+    std::vector<vector_fp> m_ec;
+    std::vector<bool> m_lindep; //!< Vector holding boolean for linear dependence
     vector_fp m_mc; //!< Vector holding coverage-specific power-law exponents
 
 private:
@@ -434,7 +449,7 @@ class InterfaceRate : public RateType, public InterfaceRateBase
         return RateType::activationEnergy() + m_ecov * GasConstant;
     }
 
-    void addCoverageDependence(const string& sp, double a, double m, double e) override {
+    void addCoverageDependence(const string& sp, double a, double m, vector_fp e) override {
         InterfaceRateBase::addCoverageDependence(sp, a, m, e);
         RateType::setCompositionDependence(true);
     }

src/kinetics/InterfaceRate.cpp

@@ -9,6 +9,7 @@
 #include "cantera/thermo/ThermoPhase.h"
 #include "cantera/thermo/SurfPhase.h"
 #include "cantera/base/AnyMap.h"
+#include "cantera/base/utilities.h"
 
 namespace Cantera
 {
@@ -140,18 +141,38 @@ void InterfaceRateBase::setCoverageDependencies(const AnyMap& dependencies,
     m_ac.clear();
     m_ec.clear();
     m_mc.clear();
+    m_lindep.clear();
     for (const auto& [species, coeffs] : dependencies) {
-        double a, E, m;
+        double a, m;
+        vector_fp E(5, 0.0);
         if (coeffs.is<AnyMap>()) {
             auto& cov_map = coeffs.as<AnyMap>();
             a = cov_map["a"].asDouble();
             m = cov_map["m"].asDouble();
-            E = units.convertActivationEnergy(cov_map["E"], "K");
+            if (cov_map["E"].isScalar()) {
+                m_lindep.push_back(true);
+                E[1] = units.convertActivationEnergy(cov_map["E"], "K");
+            } else {
+                m_lindep.push_back(false);
+                auto& E_temp = cov_map["E"].asVector<AnyValue>(4);
+                for (size_t i = 0; i < E_temp.size(); i++) {
+                    E[i+1] = units.convertActivationEnergy(E_temp[i], "K");
+                }
+            }
         } else {
             auto& cov_vec = coeffs.asVector<AnyValue>(3);
             a = cov_vec[0].asDouble();
             m = cov_vec[1].asDouble();
-            E = units.convertActivationEnergy(cov_vec[2], "K");
+            if (cov_vec[2].isScalar()) {
+                m_lindep.push_back(true);
+                E[1] = units.convertActivationEnergy(cov_vec[2], "K");
+            } else {
+                m_lindep.push_back(false);
+                auto& E_temp = cov_vec[2].asVector<AnyValue>(4);
+                for (size_t i = 0; i < E_temp.size(); i++) {
+                    E[i+1] = units.convertActivationEnergy(E_temp[i], "K");
+                }
+            }
         }
         addCoverageDependence(species, a, m, E);
     }
@@ -166,22 +187,39 @@ void InterfaceRateBase::getCoverageDependencies(AnyMap& dependencies,
             warn_deprecated("InterfaceRateBase::getCoverageDependencies",
                 "To be changed after Cantera 3.0: second argument will be removed.");
             vector_fp dep(3);
+            if (m_lindep[k]) {
+                dep[2] = m_ec[k][1];
+            } else {
+                std::vector<AnyValue> dep(3);
+                vector_fp E_temp(4);
+                for (size_t i = 0; i < m_ec[k].size() - 1; i++) {
+                    E_temp[i] = m_ec[k][i+1];
+                }
+                dep[2] = E_temp;
+            }
             dep[0] = m_ac[k];
             dep[1] = m_mc[k];
-            dep[2] = m_ec[k];
             dependencies[m_cov[k]] = std::move(dep);
         } else {
             AnyMap dep;
             dep["a"] = m_ac[k];
             dep["m"] = m_mc[k];
-            dep["E"].setQuantity(m_ec[k], "K", true);
+            if (m_lindep[k]) {
+                dep["E"].setQuantity(m_ec[k][1], "K", true);
+            } else {
+                std::vector<AnyValue> E_temp(4);
+                for (size_t i = 0; i < m_ec[k].size() - 1; i++) {
+                    E_temp[i].setQuantity(m_ec[k][i+1], "K", true);
+                }
+                dep["E"] = E_temp;
+            }
             dependencies[m_cov[k]] = std::move(dep);
         }
     }
 }
 
 void InterfaceRateBase::addCoverageDependence(const std::string& sp,
-                                              double a, double m, double e)
+                                              double a, double m, vector_fp e)
 {
     if (std::find(m_cov.begin(), m_cov.end(), sp) == m_cov.end()) {
         m_cov.push_back(sp);
@@ -226,7 +264,7 @@ void InterfaceRateBase::updateFromStruct(const InterfaceData& shared_data) {
     m_mcov = 0.0;
     for (auto& [iCov, iKin] : m_indices) {
         m_acov += m_ac[iCov] * shared_data.coverages[iKin];
-        m_ecov += m_ec[iCov] * shared_data.coverages[iKin];
+        m_ecov += poly4(shared_data.coverages[iKin], m_ec[iCov].data());
         m_mcov += m_mc[iCov] * shared_data.logCoverages[iKin];
     }
 

test/data/surface-phases.yaml

@@ -148,6 +148,10 @@ Pt-reactions:
 Pt-multisite-reactions:
 - equation: 2 O(s) <=> O2(s)
   rate-constant: {A: 3.1e9 cm^2/mol/s, b: 0.0, Ea: 0.0}
+  # dummy reaction for polynomial Ea coverage dependence
+  coverage-dependencies:
+    {O(s): {a: 0, m: 0, E: [1.0e3 J/mol, 3.0e3 J/mol,
+                            -7.0e4 J/mol, 5.0e3 J/mol]}}
 
 graphite-anode-species:
 - name: EC(e)

test/kinetics/kineticsFromYaml.cpp

@@ -605,8 +605,23 @@ TEST(Kinetics, InterfaceKineticsFromYaml)
     auto& cov_map = coverage_deps["H(s)"].as<Cantera::AnyMap>();
     EXPECT_DOUBLE_EQ(cov_map["E"].asDouble(), -6e6);
 
+
     auto R3 = kin->reaction(2);
     EXPECT_TRUE(std::dynamic_pointer_cast<StickingArrheniusRate>(R3->rate()));
+
+    auto soln2 = newInterface("surface-phases.yaml", "Pt-multi-sites");
+    auto kin2 = soln2->kinetics();
+    auto R4 = kin2->reaction(3);
+    vector_fp coeffs{1.0e6, 3.0e6, -7.0e7, 5.0e6};
+    const auto rate4 = std::dynamic_pointer_cast<InterfaceArrheniusRate>(R4->rate());
+    Cantera::AnyMap coverage_deps2;
+    rate4->getCoverageDependencies(coverage_deps2);
+    coverage_deps2.applyUnits();
+    auto& cov_map2 = coverage_deps2["O(s)"].as<Cantera::AnyMap>();
+    auto& poly_dep = cov_map2["E"].asVector<AnyValue>();
+    for (size_t i = 0; i < poly_dep.size(); i++) {
+        EXPECT_DOUBLE_EQ(poly_dep[i].asDouble(), coeffs[i]);
+    }
 }
 
 TEST(Kinetics, BMInterfaceKineticsFromYaml)