[Kinetics] Always use temperature of "reacting" phase

This is always the lowest-dimensional phase, e.g. surface or edge for
heterogeneous systems.

src/kinetics/ImplicitSurfChem.cpp

@@ -180,7 +180,7 @@ void ImplicitSurfChem::solvePseudoSteadyStateProblem(int ifuncOverride,
     //  2) Temperature and pressure
     getConcSpecies(m_concSpecies.data());
     InterfaceKinetics* ik = m_vecKinPtrs[0];
-    ThermoPhase& tp = ik->thermo(0);
+    ThermoPhase& tp = ik->thermo(ik->reactionPhaseIndex());
     doublereal TKelvin = tp.temperature();
     doublereal PGas = tp.pressure();
 

src/kinetics/InterfaceKinetics.cpp

@@ -129,7 +129,7 @@ void InterfaceKinetics::updateKc()
          * and m_mu0_Kc[]
          */
         updateMu0();
-        doublereal rrt = 1.0 / thermo(0).RT();
+        doublereal rrt = 1.0 / thermo(reactionPhaseIndex()).RT();
 
         // compute Delta mu^0 for all reversible reactions
         getRevReactionDelta(m_mu0_Kc.data(), m_rkcn.data());
@@ -159,7 +159,8 @@ void InterfaceKinetics::updateMu0()
         thermo(n).getStandardChemPotentials(m_mu0.data() + m_start[n]);
         for (size_t k = 0; k < thermo(n).nSpecies(); k++) {
             m_mu0_Kc[ik] = m_mu0[ik] + Faraday * m_phi[n] * thermo(n).charge(k);
-            m_mu0_Kc[ik] -= thermo(0).RT() * thermo(n).logStandardConc(k);
+            m_mu0_Kc[ik] -= thermo(reactionPhaseIndex()).RT()
+                            * thermo(n).logStandardConc(k);
             ik++;
         }
     }
@@ -168,7 +169,7 @@ void InterfaceKinetics::updateMu0()
 void InterfaceKinetics::getEquilibriumConstants(doublereal* kc)
 {
     updateMu0();
-    doublereal rrt = 1.0 / thermo(0).RT();
+    doublereal rrt = 1.0 / thermo(reactionPhaseIndex()).RT();
     std::fill(kc, kc + nReactions(), 0.0);
     getReactionDelta(m_mu0_Kc.data(), kc);
     for (size_t i = 0; i < nReactions(); i++) {
@@ -239,7 +240,7 @@ void InterfaceKinetics::applyVoltageKfwdCorrection(doublereal* const kf)
         if (m_ctrxn_BVform[i] == 0) {
             double eamod = m_beta[i] * deltaElectricEnergy_[irxn];
             if (eamod != 0.0) {
-                kf[irxn] *= exp(-eamod/thermo(0).RT());
+                kf[irxn] *= exp(-eamod/thermo(reactionPhaseIndex()).RT());
             }
         }
     }
@@ -265,7 +266,8 @@ void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* con
             // come out of this calculation.
             if (m_ctrxn_BVform[i] == 0) {
                 //  Calculate the term and modify the forward reaction
-                double tmp = exp(- m_beta[i] * m_deltaG0[irxn] / thermo(0).RT());
+                double tmp = exp(- m_beta[i] * m_deltaG0[irxn]
+                                 / thermo(reactionPhaseIndex()).RT());
                 tmp *= 1.0 / m_ProdStanConcReac[irxn] / Faraday;
                 kfwd[irxn] *= tmp;
             }
@@ -280,7 +282,8 @@ void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* con
                 // Calculate the term and modify the forward reaction rate
                 // constant so that it's in the exchange current density
                 // formulation format
-                double tmp = exp(m_beta[i] * m_deltaG0[irxn] * thermo(0).RT());
+                double tmp = exp(m_beta[i] * m_deltaG0[irxn]
+                                 * thermo(reactionPhaseIndex()).RT());
                 tmp *= Faraday * m_ProdStanConcReac[irxn];
                 kfwd[irxn] *= tmp;
             }
@@ -468,7 +471,7 @@ void InterfaceKinetics::getDeltaSSEnthalpy(doublereal* deltaH)
         thermo(n).getEnthalpy_RT(m_grt.data() + m_start[n]);
     }
     for (size_t k = 0; k < m_kk; k++) {
-        m_grt[k] *= thermo(0).RT();
+        m_grt[k] *= thermo(reactionPhaseIndex()).RT();
     }
 
     // Use the stoichiometric manager to find deltaH for each reaction.

test_problems/surfSolverTest/haca2.cti

@@ -15,7 +15,7 @@ units(length = 'cm', quantity = 'mol', act_energy = 'kcal/mol')
 ideal_gas(name = 'gas',
           elements = 'O H C N Ar',
           species = 'gri30: H N2 CH3 CH4 C2H2 H2 OH H2O CO O2',
-          initial_state = state(temperature = 1400.0,
+          initial_state = state(temperature = 1100.0,
                                 pressure = OneAtm,
                     mole_fractions = 'H:0.01, N2:0.8899, H2:0.04, CH4:0.01 C2H2:0.01 \
                                       OH:0.0001 H2O:0.04 O2:0.001'))
@@ -53,7 +53,7 @@ ideal_interface(name = 'soot_interface',
                 reactions = 'all',
                 phases = 'gas soot',
                 site_density = (3.8E-9, 'mol/cm2'),
-                initial_state = state(temperature= 1000.0,
+                initial_state = state(temperature= 1400.0,
                                       coverages = 'Csoot-*:0.1, Csoot-H:0.9'))
 
 # HKM -> Note, thermo from the following source:

test_problems/surfSolverTest/surfaceSolver.cpp

@@ -371,19 +371,10 @@ int main(int argc, char** argv)
         /*  Now Tweak the inputs and do a quick calculation */
         /****************************************************************************/
 
-        /*
-         * Set the Gas State:
-         * -> note that the states are set in the XML files too
-         */
-
-        /*
-         * Set the Gas State:
-         * -> note that the states are set in the XML files too
-         */
-        pres = gasTP->pressure();
-        double temp = gasTP->temperature();
+        pres = surfPhaseTP->pressure();
+        double temp = surfPhaseTP->temperature();
         temp += 95;
-        gasTP->setState_TP(temp, pres);
+        surfPhaseTP->setState_TP(temp, pres);
 
         iKin_ptr->solvePseudoSteadyStateProblem();
         iKin_ptr->getNetProductionRates(src);
@@ -399,7 +390,7 @@ int main(int argc, char** argv)
         /*****************************************************************************/
         /*  Now Don't Tweak the inputs at all */
         /****************************************************************************/
-        gasTP->setState_TP(temp, pres);
+        surfPhaseTP->setState_TP(temp, pres);
 
         iKin_ptr->solvePseudoSteadyStateProblem();
         iKin_ptr->getNetProductionRates(src);

test_problems/surfSolverTest/surfaceSolver2.cpp

@@ -399,19 +399,10 @@ int main(int argc, char** argv)
         /*  Now Tweak the inputs and do a quick calculation */
         /****************************************************************************/
 
-        /*
-         * Set the Gas State:
-         * -> note that the states are set in the XML files too
-         */
-
-        /*
-         * Set the Gas State:
-         * -> note that the states are set in the XML files too
-         */
-        pres = gasTP->pressure();
-        double temp = gasTP->temperature();
+        pres = surfPhaseTP->pressure();
+        double temp = surfPhaseTP->temperature();
         temp += 95;
-        gasTP->setState_TP(temp, pres);
+        surfPhaseTP->setState_TP(temp, pres);
 
         surfaceProb->solvePseudoSteadyStateProblem();
         iKin_ptr->getNetProductionRates(src);