Implementation of P-log rate expressions

Also includes skeleton for implementing Chebyshev rates.

include/cantera/kinetics/GasKinetics.h

@@ -386,6 +386,11 @@ public:
 
 
     void _update_rates_T();
+
+    //! Update properties that depend on concentrations.
+    //! Currently the enhanced collision partner concentrations are updated
+    //! here, as well as the pressure-dependent portion of P-log and Chebyshev
+    //! reactions.
     void _update_rates_C();
 
     //@}
@@ -409,6 +414,9 @@ protected:
 
     std::vector<size_t> m_irrev;
 
+    Rate1<Plog> m_plog_rates;
+    Rate1<ChebyshevRate> m_cheb_rates;
+
     ReactionStoichMgr*                   m_rxnstoich;
 
     std::vector<size_t> m_fwdOrder;

include/cantera/kinetics/RateCoeffMgr.h

@@ -104,6 +104,10 @@ public:
         output1 << key;
     }
 
+    size_t nReactions() const {
+        return m_rates.size();
+    }
+
 protected:
     std::vector<R>             m_rates;
     std::vector<size_t>           m_rxn;

include/cantera/kinetics/ReactionData.h

@@ -74,7 +74,9 @@ public:
 
     //! Arrhenius parameters for P-log reactions.
     //! The keys are the pressures corresponding to each Arrhenius expression.
-    std::map<double, vector_fp> plogParameters;
+    //! Multiple sets of Arrhenius parameters may be specified at a given
+    //! pressure.
+    std::multimap<double, vector_fp> plogParameters;
 
     double chebTmin; //!< Minimum temperature for Chebyshev fit
     double chebTmax; //!< Maximum temperature for Chebyshev fit

include/cantera/kinetics/RxnRates.h

@@ -468,6 +468,239 @@ protected:
 };
 
 
+class Plog
+{
+public:
+    //! return the rate coefficient type.
+    static int type()
+    {
+        return PLOG_REACTION_RATECOEFF_TYPE;
+    }
+
+    //! Default constructor.
+    Plog() {}
+
+    //! Constructor from ReactionData.
+    explicit Plog(const ReactionData& rdata) :
+        logP1_(1000),
+        logP2_(-1000),
+        maxRates_(1)
+    {
+        typedef std::multimap<double, vector_fp>::const_iterator iter_t;
+
+        size_t j = 0;
+        size_t rateCount = 0;
+        // Insert intermediate pressures
+        for (iter_t iter = rdata.plogParameters.begin();
+             iter != rdata.plogParameters.end();
+             iter++) {
+            double logp = log(iter->first);
+            if (pressures_.empty() || pressures_.rbegin()->first != logp) {
+                // starting a new group
+                pressures_[logp] = std::make_pair(j, j+1);
+                rateCount = 1;
+            } else {
+                // another rate expression at the same pressure
+                pressures_[logp].second = j+1;
+                rateCount++;
+            }
+            maxRates_ = std::max(rateCount, maxRates_);
+
+            j++;
+            A_.push_back(iter->second[0]);
+            n_.push_back(iter->second[1]);
+            Ea_.push_back(iter->second[2]);
+        }
+
+        // For pressures with only one Arrhenius expression, it is more
+        // efficient to work with log(A)
+        for (pressureIter iter = pressures_.begin();
+             iter != pressures_.end();
+             iter++) {
+            if (iter->second.first == iter->second.second - 1) {
+                A_[iter->second.first] = log(A_[iter->second.first]);
+            }
+        }
+
+        // Duplicate the first and last groups to handle P < P_0 and P > P_N
+        pressures_.insert(std::make_pair(-1000.0, pressures_.begin()->second));
+        pressures_.insert(std::make_pair(1000.0, pressures_.rbegin()->second));
+
+        // Resize work arrays
+        A1_.resize(maxRates_);
+        A2_.resize(maxRates_);
+        n1_.resize(maxRates_);
+        n2_.resize(maxRates_);
+        Ea1_.resize(maxRates_);
+        Ea2_.resize(maxRates_);
+    }
+
+    //! Update concentration-dependent parts of the rate coefficient.
+    //! @param c natural log of the pressure in Pa
+    void update_C(const doublereal* c)
+    {
+        logP_ = c[0];
+        if (logP_ > logP1_ && logP_ < logP2_) {
+            return;
+        }
+
+        pressureIter iter = pressures_.upper_bound(c[0]);
+        AssertThrowMsg(iter != pressures_.end(), "Plog::update_C",
+                       "Pressure out of range: " + fp2str(logP));
+        AssertThrowMsg(iter != pressures.begin(), "Plog::update_C",
+                       "Pressure out of range: " + fp2str(logP));
+
+        // upper interpolation pressure
+        logP2_ = iter->first;
+        size_t start = iter->second.first;
+        m2_ = iter->second.second - start;
+        for (size_t m = 0; m < m2_; m++) {
+            A2_[m] = A_[start+m];
+            n2_[m] = n_[start+m];
+            Ea2_[m] = Ea_[start+m];
+        }
+
+        // lower interpolation pressure
+        logP1_ = (--iter)->first;
+        start = iter->second.first;
+        m1_ = iter->second.second - start;
+        for (size_t m = 0; m < m1_; m++) {
+            A1_[m] = A_[start+m];
+            n1_[m] = n_[start+m];
+            Ea1_[m] = Ea_[start+m];
+        }
+
+        rDeltaP_ = 1.0 / (logP2_ - logP1_);
+    }
+
+    /**
+     * Update the value of the logarithm of the rate constant.
+     */
+    doublereal update(doublereal logT, doublereal recipT) const
+    {
+        double log_k1, log_k2;
+        if (m1_ == 1) {
+            log_k1 = A1_[0] + n1_[0] * logT - Ea1_[0] * recipT;
+        } else {
+            double k = 0.0;
+            for (size_t m = 0; m < m1_; m++) {
+                k += A1_[m] * exp(n1_[m] * logT - Ea1_[m] * recipT);
+            }
+            log_k1 = log(k);
+        }
+
+        if (m2_ == 1) {
+            log_k2 = A2_[0] + n2_[0] * logT - Ea2_[0] * recipT;
+        } else {
+            double k = 0.0;
+            for (size_t m = 0; m < m2_; m++) {
+                k += A2_[m] * exp(n2_[m] * logT - Ea2_[m] * recipT);
+            }
+            log_k2 = log(k);
+        }
+
+        return log_k1 + (log_k2 - log_k1) * (logP_ - logP1_) * rDeltaP_;
+    }
+
+    /**
+     * Update the value the rate constant.
+     *
+     * This function returns the actual value of the rate constant.
+     */
+    doublereal updateRC(doublereal logT, doublereal recipT) const {
+        return exp(update(logT, recipT));
+    }
+
+    doublereal activationEnergy_R() const {
+        throw CanteraError("Plog::activationEnergy_R", "Not implemented");
+    }
+
+    static bool alwaysComputeRate() {
+        return false;
+    }
+
+protected:
+    //! log(p) to (index range) in A_, n, Ea vectors
+    std::map<double, std::pair<size_t, size_t> > pressures_;
+    typedef std::map<double, std::pair<size_t, size_t> >::iterator pressureIter;
+
+    vector_fp A_; //!< Pre-exponential factor at each pressure (or log(A))
+    vector_fp n_; //!< Temperature exponent at each pressure [dimensionless]
+    vector_fp Ea_; //!< Activation energy at each pressure [K]
+
+    double logP_; //!< log(p) at the current state
+    double logP1_, logP2_; //!< log(p) at the lower / upper pressure reference
+
+    //! Pre-exponential factors at lower / upper pressure reference.
+    //! Stored as log(A) when there is only one at the corresponding pressure.
+    vector_fp A1_, A2_;
+    vector_fp n1_, n2_; //!< n at lower / upper pressure reference
+    vector_fp Ea1_, Ea2_; //!< Activation energy at lower / upper pressure reference
+
+    //! Number of Arrhenius expressions at lower / upper pressure references
+    size_t m1_, m2_;
+    double rDeltaP_; //!< reciprocal of (logP2 - logP1)
+
+    size_t maxRates_; //!< The maximum number of rates at any given pressure
+};
+
+
+class ChebyshevRate
+{
+public:
+    //! return the rate coefficient type.
+    static int type()
+    {
+        return CHEBYSHEV_REACTION_RATECOEFF_TYPE;
+    }
+
+    //! Default constructor.
+    ChebyshevRate() {}
+
+    //! Constructor from ReactionData.
+    explicit ChebyshevRate(const ReactionData& rdata)
+    {
+    }
+
+    //! Update concentration-dependent parts of the rate coefficient.
+    //! @param c natural log of the pressure in Pa
+    void update_C(const doublereal* c)
+    {
+    }
+
+    /**
+     * Update the value of the logarithm of the rate constant.
+     *
+     * Note, this function should never be called for negative A values.
+     * If it does then it will produce a negative overflow result, and
+     * a zero net forwards reaction rate, instead of a negative reaction
+     * rate constant that is the expected result.
+     */
+    doublereal update(doublereal logT, doublereal recipT) const
+    {
+        return 0.0;
+    }
+
+    /**
+     * Update the value the rate constant.
+     *
+     * This function returns the actual value of the rate constant.
+     */
+    doublereal updateRC(doublereal logT, doublereal recipT) const {
+        return exp(update(logT, recipT));
+    }
+
+    doublereal activationEnergy_R() const {
+        return 0.0;
+    }
+
+    static bool alwaysComputeRate() {
+        return false;
+    }
+
+protected:
+};
+
 //     class LandauTeller {
 
 //     public:

src/kinetics/GasKinetics.cpp

@@ -129,6 +129,8 @@ GasKinetics& GasKinetics::operator=(const GasKinetics& right)
     m_3b_concm = right.m_3b_concm;
     m_falloff_concm = right.m_falloff_concm;
     m_irrev = right.m_irrev;
+    m_plog_rates = right.m_plog_rates;
+    m_cheb_rates = right.m_cheb_rates;
 
     *m_rxnstoich = *(right.m_rxnstoich);
 
@@ -195,6 +197,7 @@ _update_rates_T()
     if (!m_kdata->m_rfn.empty()) {
         m_rates.update(T, logT, &m_kdata->m_rfn[0]);
     }
+
     if (!m_kdata->m_rfn_low.empty()) {
         m_falloff_low_rates.update(T, logT, &m_kdata->m_rfn_low[0]);
         m_falloff_high_rates.update(T, logT, &m_kdata->m_rfn_high[0]);
@@ -202,29 +205,50 @@ _update_rates_T()
     if (!m_kdata->falloff_work.empty()) {
         m_falloffn.updateTemp(T, &m_kdata->falloff_work[0]);
     }
+    if (m_plog_rates.nReactions()) {
+        m_plog_rates.update(T, logT, &m_kdata->m_rfn[0]);
+    }
+
+    if (m_cheb_rates.nReactions()) {
+        m_cheb_rates.update(T, logT, &m_kdata->m_rfn[0]);
+    }
+
     m_kdata->m_temp = T;
     updateKc();
     m_kdata->m_ROP_ok = false;
-    //}
 };
 
 //====================================================================================================================
-/**
- * Update properties that depend on concentrations. Currently only
- * the enhanced collision partner concentrations are updated here.
- */
+
 void GasKinetics::
 _update_rates_C()
 {
     thermo().getActivityConcentrations(&m_conc[0]);
     doublereal ctot = thermo().molarDensity();
+
+    // 3-body reactions
     if (!m_kdata->concm_3b_values.empty()) {
         m_3b_concm.update(m_conc, ctot, &m_kdata->concm_3b_values[0]);
     }
+
+    // Falloff reactions
     if (!m_kdata->concm_falloff_values.empty()) {
         m_falloff_concm.update(m_conc, ctot,
                                &m_kdata->concm_falloff_values[0]);
     }
+
+    double logP = log(thermo().pressure());
+
+    // P-log reactions
+    if (m_plog_rates.nReactions()) {
+        m_plog_rates.update_C(&logP);
+    }
+
+    // Chebyshev reactions
+    if (m_cheb_rates.nReactions()) {
+        m_cheb_rates.update_C(&logP);
+    }
+
     m_kdata->m_ROP_ok = false;
 }
 //====================================================================================================================
@@ -523,9 +547,8 @@ void GasKinetics::processFalloffReactions()
 //====================================================================================================================
 void GasKinetics::updateROP()
 {
-
-    _update_rates_T();
     _update_rates_C();
+    _update_rates_T();
 
     if (m_kdata->m_ROP_ok) {
         return;
@@ -587,8 +610,8 @@ void GasKinetics::updateROP()
 void GasKinetics::
 getFwdRateConstants(doublereal* kfwd)
 {
-    _update_rates_T();
     _update_rates_C();
+    _update_rates_T();
 
     // copy rate coefficients into ropf
     const vector_fp& rf = m_kdata->m_rfn;
@@ -761,17 +784,30 @@ addThreeBodyReaction(ReactionData& r)
 
 void GasKinetics::addPlogReaction(ReactionData& r)
 {
-    // @todo: Not yet implemented
+    // install rate coefficient calculator
+    size_t iloc = m_plog_rates.install(reactionNumber(), r);
+
+    // add a dummy entry in m_rfn, where computed rate coeff will be put
+    m_kdata->m_rfn.push_back(0.0);
+
+    m_fwdOrder.push_back(r.reactants.size());
+    registerReaction(reactionNumber(), PLOG_RXN, iloc);
 }
 
 void GasKinetics::addChebyshevReaction(ReactionData& r)
 {
-    // @todo: Not yet implemented
+    // install rate coefficient calculator
+    size_t iloc = m_cheb_rates.install(reactionNumber(), r);
+
+    // add a dummy entry in m_rfn, where computed rate coeff will be put
+    m_kdata->m_rfn.push_back(0.0);
+
+    m_fwdOrder.push_back(r.reactants.size());
+    registerReaction(reactionNumber(), CHEBYSHEV_RXN, iloc);
 }
 
 void GasKinetics::installReagents(const ReactionData& r)
 {
-
     m_kdata->m_ropf.push_back(0.0);     // extend by one for new rxn
     m_kdata->m_ropr.push_back(0.0);
     m_kdata->m_ropnet.push_back(0.0);

src/kinetics/importKinetics.cpp

@@ -509,7 +509,8 @@ void getRateCoefficient(const XML_Node& kf, Kinetics& kin,
         for (size_t m = 0; m < kf.nChildren(); m++) {
             const XML_Node& node = kf.child(m);
             double p = getFloat(node, "P", "toSI");
-            vector_fp& rate = rdata.plogParameters[p];
+            vector_fp& rate = rdata.plogParameters.insert(
+                std::make_pair(p, vector_fp()))->second;
             rate.resize(3);
             rate[0] = getFloat(node, "A", "toSI");
             rate[1] = getFloat(node, "b");