[Kinetics] Implement direct constructor for P-log reactions

include/cantera/kinetics/RxnRates.h

@@ -17,6 +17,8 @@
 namespace Cantera
 {
 
+class Plog;
+
 //! Arrhenius reaction rate type depends only on temperature
 /**
  * A reaction rate coefficient of the following form.
@@ -104,6 +106,8 @@ public:
 
 protected:
     doublereal m_logA, m_b, m_E, m_A;
+
+    friend class Plog;
 };
 
 
@@ -303,6 +307,9 @@ public:
     //! Constructor from ReactionData.
     explicit Plog(const ReactionData& rdata);
 
+    //! Constructor from Arrhenius rate expressions at a set of pressures
+    explicit Plog(const std::multimap<double, Arrhenius>& rates);
+
     //! Update concentration-dependent parts of the rate coefficient.
     //! @param c natural log of the pressure in Pa
     void update_C(const doublereal* c) {

src/kinetics/RxnRates.cpp

@@ -181,6 +181,62 @@ Plog::Plog(const ReactionData& rdata)
     }
 }
 
+Plog::Plog(const std::multimap<double, Arrhenius>& rates)
+    : logP_(-1000)
+    , logP1_(1000)
+    , logP2_(-1000)
+    , m1_(npos)
+    , m2_(npos)
+    , rDeltaP_(-1.0)
+    , maxRates_(1)
+{
+    size_t j = 0;
+    size_t rateCount = 0;
+    // Insert intermediate pressures
+    for (std::multimap<double, Arrhenius>::const_iterator iter = rates.begin();
+            iter != rates.end();
+            iter++) {
+        double logp = std::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.m_A);
+        n_.push_back(iter->second.m_b);
+        Ea_.push_back(iter->second.m_E);
+    }
+
+    // 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] = std::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_);
+}
+
 void Plog::validate(const ReactionData& rdata)
 {
     double T[] = {200.0, 500.0, 1000.0, 2000.0, 5000.0, 10000.0};

test/data/kineticsfromscratch.cti

@@ -18,3 +18,9 @@ falloff_reaction('2 OH (+ M) <=> H2O2 (+ M)',
                  kf0=[2.300000e+12, -0.9, -1700.0],
                  efficiencies='AR:0.7 H2:2.0 H2O:6.0',
                  falloff=Troe(A=0.7346, T3=94.0, T1=1756.0, T2=5182.0))
+
+pdep_arrhenius('H2 + O2 <=> 2 OH',
+               [(0.01, 'atm'), 1.212400e+16, -0.5779, 10872.7],
+               [(1.0, 'atm'), 4.910800e+31, -4.8507, 24772.8],
+               [(10.0, 'atm'), 1.286600e+47, -9.0246, 39796.5],
+               [(100.0, 'atm'), 5.963200e+56, -11.529, 52599.6])

test/kinetics/kineticsFromScratch.cpp

@@ -101,3 +101,26 @@ TEST_F(KineticsFromScratch, add_falloff_reaction)
     kin.finalize();
     check_rates(2);
 }
+
+TEST_F(KineticsFromScratch, add_plog_reaction)
+{
+    // reaction 3:
+    // pdep_arrhenius('H2 + O2 <=> 2 OH',
+    //                [(0.01, 'atm'), 1.212400e+16, -0.5779, 10872.7],
+    //                [(1.0, 'atm'), 4.910800e+31, -4.8507, 24772.8],
+    //                [(10.0, 'atm'), 1.286600e+47, -9.0246, 39796.5],
+    //                [(100.0, 'atm'), 5.963200e+56, -11.529, 52599.6])
+    Composition reac = parseCompString("H2:1, O2:1");
+    Composition prod = parseCompString("OH:2");
+    std::multimap<double, Arrhenius> rates;
+    typedef std::multimap<double, Arrhenius>::value_type item;
+    rates.insert(item(0.01*101325, Arrhenius(1.212400e+16, -0.5779, 10872.7 / GasConst_cal_mol_K)));
+    rates.insert(item(1.0*101325, Arrhenius(4.910800e+31, -4.8507, 24772.8 / GasConst_cal_mol_K)));
+    rates.insert(item(10.0*101325, Arrhenius(1.286600e+47, -9.0246, 39796.5 / GasConst_cal_mol_K)));
+    rates.insert(item(100.0*101325, Arrhenius(5.963200e+56, -11.529, 52599.6 / GasConst_cal_mol_K)));
+
+    shared_ptr<PlogReaction> R(new PlogReaction(reac, prod, Plog(rates)));
+    kin.addReaction(R);
+    kin.finalize();
+    check_rates(3);
+}