From ceefc5ecb05dd220db1081f326f9c016a38b4794 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sat, 1 Aug 2015 23:23:47 -0400 Subject: [PATCH] Consistently use ThermoPhase::RT() --- include/cantera/thermo/ConstDensityThermo.h | 2 +- include/cantera/thermo/ThermoPhase.h | 14 ++++- src/equil/ChemEquil.cpp | 6 +-- src/kinetics/AqueousKinetics.cpp | 7 +-- src/kinetics/BulkKinetics.cpp | 6 +-- src/kinetics/InterfaceKinetics.cpp | 25 +++------ src/thermo/ConstDensityThermo.cpp | 3 +- src/thermo/DebyeHuckel.cpp | 33 ++++-------- src/thermo/FixedChemPotSSTP.cpp | 12 ++--- src/thermo/HMWSoln.cpp | 36 +++++-------- src/thermo/IdealGasPhase.cpp | 20 +++----- src/thermo/IdealMolalSoln.cpp | 12 ++--- src/thermo/IdealSolidSolnPhase.cpp | 12 ++--- src/thermo/IdealSolnGasVPSS.cpp | 6 +-- src/thermo/IonsFromNeutralVPSSTP.cpp | 10 ++-- src/thermo/LatticePhase.cpp | 12 ++--- src/thermo/MargulesVPSSTP.cpp | 35 +++++-------- src/thermo/MaskellSolidSolnPhase.cpp | 8 ++- src/thermo/MetalSHEelectrons.cpp | 3 +- src/thermo/MineralEQ3.cpp | 9 ++-- src/thermo/MixedSolventElectrolyte.cpp | 35 +++++-------- src/thermo/MixtureFugacityTP.cpp | 24 ++++----- src/thermo/MolarityIonicVPSSTP.cpp | 3 +- src/thermo/PhaseCombo_Interaction.cpp | 3 +- src/thermo/PureFluidPhase.cpp | 4 +- src/thermo/RedlichKisterVPSSTP.cpp | 3 +- src/thermo/RedlichKwongMFTP.cpp | 57 +++++++++------------ src/thermo/SurfPhase.cpp | 3 +- src/thermo/ThermoPhase.cpp | 3 +- src/thermo/VPStandardStateTP.cpp | 6 +-- src/thermo/WaterSSTP.cpp | 3 +- test/kinetics/rates.cpp | 2 +- 32 files changed, 160 insertions(+), 257 deletions(-) diff --git a/include/cantera/thermo/ConstDensityThermo.h b/include/cantera/thermo/ConstDensityThermo.h index 37c86aec6..808dcc006 100644 --- a/include/cantera/thermo/ConstDensityThermo.h +++ b/include/cantera/thermo/ConstDensityThermo.h @@ -153,7 +153,7 @@ public: */ virtual void getPureGibbs(doublereal* gpure) const { const vector_fp& gibbsrt = gibbs_RT(); - scale(gibbsrt.begin(), gibbsrt.end(), gpure, _RT()); + scale(gibbsrt.begin(), gibbsrt.end(), gpure, RT()); } //! Get the nondimensional Enthalpy functions for the species diff --git a/include/cantera/thermo/ThermoPhase.h b/include/cantera/thermo/ThermoPhase.h index 618f93afc..eba0eb6cf 100644 --- a/include/cantera/thermo/ThermoPhase.h +++ b/include/cantera/thermo/ThermoPhase.h @@ -12,6 +12,7 @@ #include "Phase.h" #include "SpeciesThermo.h" +#include "cantera/base/global.h" namespace Cantera { @@ -865,9 +866,20 @@ public: //! Return the Gas Constant multiplied by the current temperature /*! - * The units are Joules kmol-1 + * The units are Joules kmol-1. + * @deprecated use RT() instead. To be removed after Cantera 2.3. */ doublereal _RT() const { + warn_deprecated("ThermoPhase::_RT()", + "use RT() instead. To be removed after Cantera 2.3."); + return temperature() * GasConstant; + } + + //! Return the Gas Constant multiplied by the current temperature + /*! + * The units are Joules kmol-1 + */ + doublereal RT() const { return temperature() * GasConstant; } diff --git a/src/equil/ChemEquil.cpp b/src/equil/ChemEquil.cpp index b79ebcf9c..575199011 100644 --- a/src/equil/ChemEquil.cpp +++ b/src/equil/ChemEquil.cpp @@ -552,12 +552,11 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, if (useThermoPhaseElementPotentials) { bool haveEm = s.getElementPotentials(DATA_PTR(x)); if (haveEm) { - doublereal rt = GasConstant * s.temperature(); if (s.temperature() < 100.) { printf("we are here %g\n", s.temperature()); } for (m = 0; m < m_mm; m++) { - x[m] /= rt; + x[m] *= 1.0 / s.RT(); } } else { estimateElementPotentials(s, x, elMolesGoal); @@ -663,9 +662,8 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, if (iter > 0 && passThis && fabs(deltax) < options.relTolerance && fabs(deltay) < options.relTolerance) { options.iterations = iter; - doublereal rt = GasConstant* s.temperature(); for (m = 0; m < m_mm; m++) { - m_lambda[m] = x[m]*rt; + m_lambda[m] = x[m]* s.RT(); } if (m_eloc != npos) { diff --git a/src/kinetics/AqueousKinetics.cpp b/src/kinetics/AqueousKinetics.cpp index a3a6e21a2..b4ce9ec14 100644 --- a/src/kinetics/AqueousKinetics.cpp +++ b/src/kinetics/AqueousKinetics.cpp @@ -48,13 +48,11 @@ void AqueousKinetics::_update_rates_C() void AqueousKinetics::updateKc() { - doublereal rt = GasConstant * m_temp; - thermo().getStandardChemPotentials(&m_grt[0]); fill(m_rkcn.begin(), m_rkcn.end(), 0.0); for (size_t k = 0; k < thermo().nSpecies(); k++) { doublereal logStandConc_k = thermo().logStandardConc(k); - m_grt[k] -= rt * logStandConc_k; + m_grt[k] -= GasConstant * m_temp * logStandConc_k; } // compute Delta G^0 for all reversible reactions @@ -77,10 +75,9 @@ void AqueousKinetics::getEquilibriumConstants(doublereal* kc) thermo().getStandardChemPotentials(&m_grt[0]); fill(m_rkcn.begin(), m_rkcn.end(), 0.0); - doublereal rt = GasConstant * m_temp; for (size_t k = 0; k < thermo().nSpecies(); k++) { doublereal logStandConc_k = thermo().logStandardConc(k); - m_grt[k] -= rt * logStandConc_k; + m_grt[k] -= GasConstant * m_temp * logStandConc_k; } // compute Delta G^0 for all reactions diff --git a/src/kinetics/BulkKinetics.cpp b/src/kinetics/BulkKinetics.cpp index ecccf47e6..838999713 100644 --- a/src/kinetics/BulkKinetics.cpp +++ b/src/kinetics/BulkKinetics.cpp @@ -64,9 +64,8 @@ void BulkKinetics::getDeltaSSEnthalpy(doublereal* deltaH) { // Get the standard state enthalpies of the species. thermo().getEnthalpy_RT(&m_grt[0]); - doublereal RT = thermo().temperature() * GasConstant; for (size_t k = 0; k < m_kk; k++) { - m_grt[k] *= RT; + m_grt[k] *= thermo().RT(); } // Use the stoichiometric manager to find deltaH for each reaction. getReactionDelta(&m_grt[0], deltaH); @@ -78,9 +77,8 @@ void BulkKinetics::getDeltaSSEntropy(doublereal* deltaS) // the entropies of the pure species at the temperature and pressure of the // solution. thermo().getEntropy_R(&m_grt[0]); - doublereal R = GasConstant; for (size_t k = 0; k < m_kk; k++) { - m_grt[k] *= R; + m_grt[k] *= GasConstant; } // Use the stoichiometric manager to find deltaS for each reaction. getReactionDelta(&m_grt[0], deltaS); diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index b33470a3d..40bda7bbd 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -240,14 +240,13 @@ void InterfaceKinetics::updateMu0() * kinetics object and store it in m_mu0[] and in m_mu0_Kc[] */ size_t nsp, ik = 0; - doublereal rt = GasConstant * thermo(0).temperature(); size_t np = nPhases(); for (size_t n = 0; n < np; n++) { thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]); nsp = thermo(n).nSpecies(); for (size_t k = 0; k < nsp; k++) { m_mu0_Kc[ik] = m_mu0[ik] + Faraday * m_phi[n] * thermo(n).charge(k); - m_mu0_Kc[ik] -= rt * thermo(n).logStandardConc(k); + m_mu0_Kc[ik] -= thermo(0).RT() * thermo(n).logStandardConc(k); ik++; } } @@ -261,8 +260,7 @@ void InterfaceKinetics::checkPartialEquil() vector_fp dmu(nTotalSpecies(), 0.0); vector_fp rmu(std::max(nReactions(), 1), 0.0); if (m_nrev > 0) { - doublereal rt = GasConstant*thermo(0).temperature(); - cout << "T = " << thermo(0).temperature() << " " << rt << endl; + cout << "T = " << thermo(0).temperature() << " " << thermo(0).RT() << endl; size_t nsp, ik=0; doublereal delta; for (size_t n = 0; n < nPhases(); n++) { @@ -281,7 +279,7 @@ void InterfaceKinetics::checkPartialEquil() for (size_t i = 0; i < m_nrev; i++) { size_t irxn = m_revindex[i]; cout << "Reaction " << reactionString(irxn) - << " " << rmu[irxn]/rt << endl; + << " " << rmu[irxn]/thermo(0).RT() << endl; printf("%12.6e %12.6e %12.6e %12.6e \n", m_ropf[irxn], m_ropr[irxn], m_ropnet[irxn], m_ropnet[irxn]/(m_ropf[irxn] + m_ropr[irxn])); @@ -372,9 +370,7 @@ void InterfaceKinetics::applyVoltageKfwdCorrection(doublereal* const kf) if (m_ctrxn_BVform[i] == 0) { eamod = m_beta[i] * deltaElectricEnergy_[irxn]; if (eamod != 0.0) { - doublereal rt = GasConstant*thermo(0).temperature(); - doublereal rrt = 1.0/rt; - kf[irxn] *= exp(-eamod*rrt); + kf[irxn] *= exp(-eamod/thermo(0).RT()); } } } @@ -383,9 +379,6 @@ void InterfaceKinetics::applyVoltageKfwdCorrection(doublereal* const kf) void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* const kfwd) { updateExchangeCurrentQuantities(); - doublereal rt = GasConstant * thermo(0).temperature(); - doublereal rrt = 1.0/rt; - // Loop over all reactions which are defined to have a voltage transfer coefficient that // affects the activity energy for the reaction for (size_t i = 0; i < m_ctrxn.size(); i++) { @@ -400,7 +393,7 @@ void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* con // We need to have the straight chemical reaction rate constant to 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] * rrt); + double tmp = exp(- m_beta[i] * m_deltaG0[irxn] / thermo(0).RT()); double tmp2 = m_ProdStanConcReac[irxn]; tmp *= 1.0 / tmp2 / Faraday; kfwd[irxn] *= tmp; @@ -415,7 +408,7 @@ void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* con if (m_ctrxn_BVform[i] != 0) { // 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] * rrt); + double tmp = exp(m_beta[i] * m_deltaG0[irxn] * thermo(0).RT()); double tmp2 = m_ProdStanConcReac[irxn]; tmp *= Faraday * tmp2; kfwd[irxn] *= tmp; @@ -647,9 +640,8 @@ void InterfaceKinetics::getDeltaSSEnthalpy(doublereal* deltaH) for (size_t n = 0; n < nPhases(); n++) { thermo(n).getEnthalpy_RT(DATA_PTR(m_grt) + m_start[n]); } - doublereal RT = thermo(0).temperature() * GasConstant; for (size_t k = 0; k < m_kk; k++) { - m_grt[k] *= RT; + m_grt[k] *= thermo(0).RT(); } /* * Use the stoichiometric manager to find deltaG for each @@ -668,9 +660,8 @@ void InterfaceKinetics::getDeltaSSEntropy(doublereal* deltaS) for (size_t n = 0; n < nPhases(); n++) { thermo(n).getEntropy_R(DATA_PTR(m_grt) + m_start[n]); } - doublereal R = GasConstant; for (size_t k = 0; k < m_kk; k++) { - m_grt[k] *= R; + m_grt[k] *= GasConstant; } /* * Use the stoichiometric manager to find deltaS for each diff --git a/src/thermo/ConstDensityThermo.cpp b/src/thermo/ConstDensityThermo.cpp index 00379e5a0..085ce3ab7 100644 --- a/src/thermo/ConstDensityThermo.cpp +++ b/src/thermo/ConstDensityThermo.cpp @@ -97,11 +97,10 @@ void ConstDensityThermo::getChemPotentials(doublereal* mu) const { doublereal vdp = (pressure() - m_spthermo->refPressure())/ molarDensity(); - doublereal rt = temperature() * GasConstant; const vector_fp& g_RT = gibbs_RT(); for (size_t k = 0; k < m_kk; k++) { double xx = std::max(SmallNumber, moleFraction(k)); - mu[k] = rt*(g_RT[k] + log(xx)) + vdp; + mu[k] = RT()*(g_RT[k] + log(xx)) + vdp; } } diff --git a/src/thermo/DebyeHuckel.cpp b/src/thermo/DebyeHuckel.cpp index c59b96c16..072128dbc 100644 --- a/src/thermo/DebyeHuckel.cpp +++ b/src/thermo/DebyeHuckel.cpp @@ -362,17 +362,16 @@ void DebyeHuckel::getChemPotentials(doublereal* mu) const * This also updates the internal molality array. */ s_update_lnMolalityActCoeff(); - doublereal RT = GasConstant * temperature(); double xmolSolvent = moleFraction(m_indexSolvent); for (size_t k = 0; k < m_kk; k++) { if (m_indexSolvent != k) { xx = std::max(m_molalities[k], SmallNumber); - mu[k] += RT * (log(xx) + m_lnActCoeffMolal[k]); + mu[k] += RT() * (log(xx) + m_lnActCoeffMolal[k]); } } xx = std::max(xmolSolvent, SmallNumber); mu[m_indexSolvent] += - RT * (log(xx) + m_lnActCoeffMolal[m_indexSolvent]); + RT() * (log(xx) + m_lnActCoeffMolal[m_indexSolvent]); } void DebyeHuckel::getPartialMolarEnthalpies(doublereal* hbar) const @@ -384,10 +383,8 @@ void DebyeHuckel::getPartialMolarEnthalpies(doublereal* hbar) const /* * Dimensionalize it. */ - double T = temperature(); - double RT = GasConstant * T; for (size_t k = 0; k < m_kk; k++) { - hbar[k] *= RT; + hbar[k] *= RT(); } /* * Check to see whether activity coefficients are temperature @@ -402,9 +399,8 @@ void DebyeHuckel::getPartialMolarEnthalpies(doublereal* hbar) const */ s_update_lnMolalityActCoeff(); s_update_dlnMolalityActCoeff_dT(); - double RTT = GasConstant * T * T; for (size_t k = 0; k < m_kk; k++) { - hbar[k] -= RTT * m_dlnActCoeffMolaldT[k]; + hbar[k] -= RT() * temperature() * m_dlnActCoeffMolaldT[k]; } } } @@ -419,9 +415,8 @@ void DebyeHuckel::getPartialMolarEntropies(doublereal* sbar) const /* * Dimensionalize the entropies */ - doublereal R = GasConstant; for (size_t k = 0; k < m_kk; k++) { - sbar[k] *= R; + sbar[k] *= GasConstant; } /* * Update the activity coefficients, This also update the @@ -436,12 +431,12 @@ void DebyeHuckel::getPartialMolarEntropies(doublereal* sbar) const for (size_t k = 0; k < m_kk; k++) { if (k != m_indexSolvent) { mm = std::max(SmallNumber, m_molalities[k]); - sbar[k] -= R * (log(mm) + m_lnActCoeffMolal[k]); + sbar[k] -= GasConstant * (log(mm) + m_lnActCoeffMolal[k]); } } double xmolSolvent = moleFraction(m_indexSolvent); mm = std::max(SmallNumber, xmolSolvent); - sbar[m_indexSolvent] -= R *(log(mm) + m_lnActCoeffMolal[m_indexSolvent]); + sbar[m_indexSolvent] -= GasConstant *(log(mm) + m_lnActCoeffMolal[m_indexSolvent]); /* * Check to see whether activity coefficients are temperature * dependent. If they are, then calculate the their temperature @@ -450,9 +445,8 @@ void DebyeHuckel::getPartialMolarEntropies(doublereal* sbar) const double dAdT = dA_DebyedT_TP(); if (dAdT != 0.0) { s_update_dlnMolalityActCoeff_dT(); - double RT = R * temperature(); for (size_t k = 0; k < m_kk; k++) { - sbar[k] -= RT * m_dlnActCoeffMolaldT[k]; + sbar[k] -= RT() * m_dlnActCoeffMolaldT[k]; } } } @@ -465,10 +459,8 @@ void DebyeHuckel::getPartialMolarVolumes(doublereal* vbar) const */ s_update_lnMolalityActCoeff(); s_update_dlnMolalityActCoeff_dP(); - double T = temperature(); - double RT = GasConstant * T; for (size_t k = 0; k < m_kk; k++) { - vbar[k] += RT * m_dlnActCoeffMolaldP[k]; + vbar[k] += RT() * m_dlnActCoeffMolaldP[k]; } } @@ -498,12 +490,9 @@ void DebyeHuckel::getPartialMolarCp(doublereal* cpbar) const s_update_lnMolalityActCoeff(); s_update_dlnMolalityActCoeff_dT(); s_update_d2lnMolalityActCoeff_dT2(); - double T = temperature(); - double RT = GasConstant * T; - double RTT = RT * T; for (size_t k = 0; k < m_kk; k++) { - cpbar[k] -= (2.0 * RT * m_dlnActCoeffMolaldT[k] + - RTT * m_d2lnActCoeffMolaldT2[k]); + cpbar[k] -= (2.0 * RT() * m_dlnActCoeffMolaldT[k] + + RT() * temperature() * m_d2lnActCoeffMolaldT2[k]); } } } diff --git a/src/thermo/FixedChemPotSSTP.cpp b/src/thermo/FixedChemPotSSTP.cpp index 93aa33264..e021abb4f 100644 --- a/src/thermo/FixedChemPotSSTP.cpp +++ b/src/thermo/FixedChemPotSSTP.cpp @@ -204,8 +204,7 @@ void FixedChemPotSSTP::getStandardChemPotentials(doublereal* mu0) const void FixedChemPotSSTP::getEnthalpy_RT(doublereal* hrt) const { - double rt = _RT(); - hrt[0] = chemPot_ / rt; + hrt[0] = chemPot_ / RT(); } void FixedChemPotSSTP::getEntropy_R(doublereal* sr) const @@ -215,8 +214,7 @@ void FixedChemPotSSTP::getEntropy_R(doublereal* sr) const void FixedChemPotSSTP::getGibbs_RT(doublereal* grt) const { - double rt = _RT(); - grt[0] = chemPot_ / rt; + grt[0] = chemPot_ / RT(); } void FixedChemPotSSTP::getCp_R(doublereal* cpr) const @@ -245,8 +243,7 @@ void FixedChemPotSSTP::getIntEnergy_RT_ref(doublereal* urt) const void FixedChemPotSSTP::getEnthalpy_RT_ref(doublereal* hrt) const { - double rt = _RT(); - hrt[0] = chemPot_ / rt; + hrt[0] = chemPot_ / RT(); } void FixedChemPotSSTP::getEntropy_R_ref(doublereal* sr) const @@ -256,8 +253,7 @@ void FixedChemPotSSTP::getEntropy_R_ref(doublereal* sr) const void FixedChemPotSSTP::getGibbs_RT_ref(doublereal* grt) const { - double rt = _RT(); - grt[0] = chemPot_ / rt; + grt[0] = chemPot_ / RT(); } void FixedChemPotSSTP::getGibbs_ref(doublereal* g) const diff --git a/src/thermo/HMWSoln.cpp b/src/thermo/HMWSoln.cpp index 3b103eeb0..97b5d497e 100644 --- a/src/thermo/HMWSoln.cpp +++ b/src/thermo/HMWSoln.cpp @@ -429,9 +429,8 @@ doublereal HMWSoln::relative_enthalpy() const getPartialMolarEnthalpies(DATA_PTR(m_tmpV)); double hbar = mean_X(m_tmpV); getEnthalpy_RT(DATA_PTR(m_gamma_tmp)); - double RT = GasConstant * temperature(); for (size_t k = 0; k < m_kk; k++) { - m_gamma_tmp[k] *= RT; + m_gamma_tmp[k] *= RT(); } double h0bar = mean_X(m_gamma_tmp); return hbar - h0bar; @@ -669,17 +668,16 @@ void HMWSoln::getChemPotentials(doublereal* mu) const * This also updates the internal molality array. */ s_update_lnMolalityActCoeff(); - doublereal RT = GasConstant * temperature(); double xmolSolvent = moleFraction(m_indexSolvent); for (size_t k = 0; k < m_kk; k++) { if (m_indexSolvent != k) { xx = std::max(m_molalities[k], SmallNumber); - mu[k] += RT * (log(xx) + m_lnActCoeffMolal_Scaled[k]); + mu[k] += RT() * (log(xx) + m_lnActCoeffMolal_Scaled[k]); } } xx = std::max(xmolSolvent, SmallNumber); mu[m_indexSolvent] += - RT * (log(xx) + m_lnActCoeffMolal_Scaled[m_indexSolvent]); + RT() * (log(xx) + m_lnActCoeffMolal_Scaled[m_indexSolvent]); } void HMWSoln::getPartialMolarEnthalpies(doublereal* hbar) const @@ -691,10 +689,8 @@ void HMWSoln::getPartialMolarEnthalpies(doublereal* hbar) const /* * dimensionalize it. */ - double T = temperature(); - double RT = GasConstant * T; for (size_t k = 0; k < m_kk; k++) { - hbar[k] *= RT; + hbar[k] *= RT(); } /* * Update the activity coefficients, This also update the @@ -702,9 +698,8 @@ void HMWSoln::getPartialMolarEnthalpies(doublereal* hbar) const */ s_update_lnMolalityActCoeff(); s_update_dlnMolalityActCoeff_dT(); - double RTT = RT * T; for (size_t k = 0; k < m_kk; k++) { - hbar[k] -= RTT * m_dlnActCoeffMolaldT_Scaled[k]; + hbar[k] -= RT() * temperature() * m_dlnActCoeffMolaldT_Scaled[k]; } } @@ -718,9 +713,8 @@ void HMWSoln::getPartialMolarEntropies(doublereal* sbar) const /* * Dimensionalize the entropies */ - doublereal R = GasConstant; for (size_t k = 0; k < m_kk; k++) { - sbar[k] *= R; + sbar[k] *= GasConstant; } /* * Update the activity coefficients, This also update the @@ -735,21 +729,20 @@ void HMWSoln::getPartialMolarEntropies(doublereal* sbar) const for (size_t k = 0; k < m_kk; k++) { if (k != m_indexSolvent) { mm = std::max(SmallNumber, m_molalities[k]); - sbar[k] -= R * (log(mm) + m_lnActCoeffMolal_Scaled[k]); + sbar[k] -= GasConstant * (log(mm) + m_lnActCoeffMolal_Scaled[k]); } } double xmolSolvent = moleFraction(m_indexSolvent); mm = std::max(SmallNumber, xmolSolvent); - sbar[m_indexSolvent] -= R *(log(mm) + m_lnActCoeffMolal_Scaled[m_indexSolvent]); + sbar[m_indexSolvent] -= GasConstant *(log(mm) + m_lnActCoeffMolal_Scaled[m_indexSolvent]); /* * Check to see whether activity coefficients are temperature * dependent. If they are, then calculate the their temperature * derivatives and add them into the result. */ s_update_dlnMolalityActCoeff_dT(); - double RT = R * temperature(); for (size_t k = 0; k < m_kk; k++) { - sbar[k] -= RT * m_dlnActCoeffMolaldT_Scaled[k]; + sbar[k] -= RT() * m_dlnActCoeffMolaldT_Scaled[k]; } } @@ -764,10 +757,8 @@ void HMWSoln::getPartialMolarVolumes(doublereal* vbar) const */ s_update_lnMolalityActCoeff(); s_update_dlnMolalityActCoeff_dP(); - double T = temperature(); - double RT = GasConstant * T; for (size_t k = 0; k < m_kk; k++) { - vbar[k] += RT * m_dlnActCoeffMolaldP_Scaled[k]; + vbar[k] += RT() * m_dlnActCoeffMolaldP_Scaled[k]; } } @@ -788,12 +779,9 @@ void HMWSoln::getPartialMolarCp(doublereal* cpbar) const s_update_lnMolalityActCoeff(); s_update_dlnMolalityActCoeff_dT(); s_update_d2lnMolalityActCoeff_dT2(); - double T = temperature(); - double RT = GasConstant * T; - double RTT = RT * T; for (size_t k = 0; k < m_kk; k++) { - cpbar[k] -= (2.0 * RT * m_dlnActCoeffMolaldT_Scaled[k] + - RTT * m_d2lnActCoeffMolaldT2_Scaled[k]); + cpbar[k] -= (2.0 * RT() * m_dlnActCoeffMolaldT_Scaled[k] + + RT() * temperature() * m_d2lnActCoeffMolaldT2_Scaled[k]); } } diff --git a/src/thermo/IdealGasPhase.cpp b/src/thermo/IdealGasPhase.cpp index 32e58ad2d..074ee6bb5 100644 --- a/src/thermo/IdealGasPhase.cpp +++ b/src/thermo/IdealGasPhase.cpp @@ -97,7 +97,7 @@ void IdealGasPhase::getActivityCoefficients(doublereal* ac) const void IdealGasPhase::getStandardChemPotentials(doublereal* muStar) const { const vector_fp& gibbsrt = gibbs_RT_ref(); - scale(gibbsrt.begin(), gibbsrt.end(), muStar, _RT()); + scale(gibbsrt.begin(), gibbsrt.end(), muStar, RT()); double tmp = log(pressure() / m_spthermo->refPressure()); tmp *= GasConstant * temperature(); for (size_t k = 0; k < m_kk; k++) { @@ -110,18 +110,16 @@ void IdealGasPhase::getStandardChemPotentials(doublereal* muStar) const void IdealGasPhase::getChemPotentials(doublereal* mu) const { getStandardChemPotentials(mu); - doublereal rt = temperature() * GasConstant; for (size_t k = 0; k < m_kk; k++) { double xx = std::max(SmallNumber, moleFraction(k)); - mu[k] += rt * log(xx); + mu[k] += RT() * log(xx); } } void IdealGasPhase::getPartialMolarEnthalpies(doublereal* hbar) const { const vector_fp& _h = enthalpy_RT_ref(); - doublereal rt = GasConstant * temperature(); - scale(_h.begin(), _h.end(), hbar, rt); + scale(_h.begin(), _h.end(), hbar, RT()); } void IdealGasPhase::getPartialMolarEntropies(doublereal* sbar) const @@ -138,9 +136,8 @@ void IdealGasPhase::getPartialMolarEntropies(doublereal* sbar) const void IdealGasPhase::getPartialMolarIntEnergies(doublereal* ubar) const { const vector_fp& _h = enthalpy_RT_ref(); - doublereal rt = GasConstant * temperature(); for (size_t k = 0; k < m_kk; k++) { - ubar[k] = rt * (_h[k] - 1.0); + ubar[k] = RT() * (_h[k] - 1.0); } } @@ -189,9 +186,8 @@ void IdealGasPhase::getGibbs_RT(doublereal* grt) const void IdealGasPhase::getPureGibbs(doublereal* gpure) const { const vector_fp& gibbsrt = gibbs_RT_ref(); - scale(gibbsrt.begin(), gibbsrt.end(), gpure, _RT()); - double tmp = log(pressure() / m_spthermo->refPressure()); - tmp *= _RT(); + scale(gibbsrt.begin(), gibbsrt.end(), gpure, RT()); + double tmp = log(pressure() / m_spthermo->refPressure()) * RT(); for (size_t k = 0; k < m_kk; k++) { gpure[k] += tmp; } @@ -236,7 +232,7 @@ void IdealGasPhase::getGibbs_RT_ref(doublereal* grt) const void IdealGasPhase::getGibbs_ref(doublereal* g) const { const vector_fp& gibbsrt = gibbs_RT_ref(); - scale(gibbsrt.begin(), gibbsrt.end(), g, _RT()); + scale(gibbsrt.begin(), gibbsrt.end(), g, RT()); } void IdealGasPhase::getEntropy_R_ref(doublereal* er) const @@ -261,7 +257,7 @@ void IdealGasPhase::getCp_R_ref(doublereal* cprt) const void IdealGasPhase::getStandardVolumes_ref(doublereal* vol) const { - doublereal tmp = _RT() / m_p0; + doublereal tmp = RT() / m_p0; for (size_t k = 0; k < m_kk; k++) { vol[k] = tmp; } diff --git a/src/thermo/IdealMolalSoln.cpp b/src/thermo/IdealMolalSoln.cpp index ce66b1f9e..4fdad07fb 100644 --- a/src/thermo/IdealMolalSoln.cpp +++ b/src/thermo/IdealMolalSoln.cpp @@ -342,12 +342,11 @@ void IdealMolalSoln::getChemPotentials(doublereal* mu) const * get the solvent mole fraction */ double xmolSolvent = moleFraction(m_indexSolvent); - doublereal RT = GasConstant * temperature(); if (IMS_typeCutoff_ == 0 || xmolSolvent > 3.* IMS_X_o_cutoff_/2.0) { for (size_t k = 1; k < m_kk; k++) { double xx = std::max(m_molalities[k], SmallNumber); - mu[k] += RT * log(xx); + mu[k] += RT() * log(xx); } /* * Do the solvent @@ -355,7 +354,7 @@ void IdealMolalSoln::getChemPotentials(doublereal* mu) const */ double xx = std::max(xmolSolvent, SmallNumber); mu[m_indexSolvent] += - (RT * (xmolSolvent - 1.0) / xx); + (RT() * (xmolSolvent - 1.0) / xx); } else { /* * Update the activity coefficients @@ -365,20 +364,19 @@ void IdealMolalSoln::getChemPotentials(doublereal* mu) const for (size_t k = 1; k < m_kk; k++) { double xx = std::max(m_molalities[k], SmallNumber); - mu[k] += RT * (log(xx) + IMS_lnActCoeffMolal_[k]); + mu[k] += RT() * (log(xx) + IMS_lnActCoeffMolal_[k]); } double xx = std::max(xmolSolvent, SmallNumber); mu[m_indexSolvent] += - RT * (log(xx) + IMS_lnActCoeffMolal_[m_indexSolvent]); + RT() * (log(xx) + IMS_lnActCoeffMolal_[m_indexSolvent]); } } void IdealMolalSoln::getPartialMolarEnthalpies(doublereal* hbar) const { getEnthalpy_RT(hbar); - doublereal RT = _RT(); for (size_t k = 0; k < m_kk; k++) { - hbar[k] *= RT; + hbar[k] *= RT(); } } diff --git a/src/thermo/IdealSolidSolnPhase.cpp b/src/thermo/IdealSolidSolnPhase.cpp index bce4887b1..6ba032378 100644 --- a/src/thermo/IdealSolidSolnPhase.cpp +++ b/src/thermo/IdealSolidSolnPhase.cpp @@ -292,11 +292,10 @@ void IdealSolidSolnPhase::getActivityCoefficients(doublereal* ac) const void IdealSolidSolnPhase::getChemPotentials(doublereal* mu) const { doublereal delta_p = m_Pcurrent - m_Pref; - doublereal RT = temperature() * GasConstant; const vector_fp& g_RT = gibbs_RT_ref(); for (size_t k = 0; k < m_kk; k++) { double xx = std::max(SmallNumber, moleFraction(k)); - mu[k] = RT * (g_RT[k] + log(xx)) + mu[k] = RT() * (g_RT[k] + log(xx)) + delta_p * m_speciesMolarVolume[k]; } } @@ -351,20 +350,18 @@ void IdealSolidSolnPhase::getPartialMolarVolumes(doublereal* vbar) const void IdealSolidSolnPhase::getPureGibbs(doublereal* gpure) const { const vector_fp& gibbsrt = gibbs_RT_ref(); - doublereal RT = _RT(); const doublereal* const gk = DATA_PTR(gibbsrt); doublereal delta_p = (m_Pcurrent - m_Pref); for (size_t k = 0; k < m_kk; k++) { - gpure[k] = RT * gk[k] + delta_p * m_speciesMolarVolume[k]; + gpure[k] = RT() * gk[k] + delta_p * m_speciesMolarVolume[k]; } } void IdealSolidSolnPhase::getGibbs_RT(doublereal* grt) const { const vector_fp& gibbsrt = gibbs_RT_ref(); - doublereal RT = _RT(); const doublereal* const gk = DATA_PTR(gibbsrt); - doublereal delta_prt = (m_Pcurrent - m_Pref)/ RT; + doublereal delta_prt = (m_Pcurrent - m_Pref)/ RT(); for (size_t k = 0; k < m_kk; k++) { grt[k] = gk[k] + delta_prt * m_speciesMolarVolume[k]; } @@ -373,8 +370,7 @@ void IdealSolidSolnPhase::getGibbs_RT(doublereal* grt) const void IdealSolidSolnPhase::getEnthalpy_RT(doublereal* hrt) const { const vector_fp& _h = enthalpy_RT_ref(); - doublereal delta_prt = ((m_Pcurrent - m_Pref) / - (GasConstant * temperature())); + doublereal delta_prt = (m_Pcurrent - m_Pref) / RT(); for (size_t k = 0; k < m_kk; k++) { hrt[k] = _h[k] + delta_prt * m_speciesMolarVolume[k]; } diff --git a/src/thermo/IdealSolnGasVPSS.cpp b/src/thermo/IdealSolnGasVPSS.cpp index dbe7ae29f..1b27b12bc 100644 --- a/src/thermo/IdealSolnGasVPSS.cpp +++ b/src/thermo/IdealSolnGasVPSS.cpp @@ -207,19 +207,17 @@ void IdealSolnGasVPSS::getActivityCoefficients(doublereal* ac) const void IdealSolnGasVPSS::getChemPotentials_RT(doublereal* muRT) const { getChemPotentials(muRT); - doublereal invRT = 1.0 / _RT(); for (size_t k = 0; k < m_kk; k++) { - muRT[k] *= invRT; + muRT[k] *= 1.0 / RT(); } } void IdealSolnGasVPSS::getChemPotentials(doublereal* mu) const { getStandardChemPotentials(mu); - doublereal rt = temperature() * GasConstant; for (size_t k = 0; k < m_kk; k++) { double xx = std::max(SmallNumber, moleFraction(k)); - mu[k] += rt*(log(xx)); + mu[k] += RT() * log(xx); } } diff --git a/src/thermo/IonsFromNeutralVPSSTP.cpp b/src/thermo/IonsFromNeutralVPSSTP.cpp index ae1970b92..707466efa 100644 --- a/src/thermo/IonsFromNeutralVPSSTP.cpp +++ b/src/thermo/IonsFromNeutralVPSSTP.cpp @@ -377,10 +377,8 @@ void IonsFromNeutralVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const /* * dimensionalize it. */ - double T = temperature(); - double RT = GasConstant * T; for (size_t k = 0; k < m_kk; k++) { - hbar[k] *= RT; + hbar[k] *= RT(); } /* * Update the activity coefficients, This also update the @@ -388,9 +386,8 @@ void IonsFromNeutralVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const */ s_update_lnActCoeff(); s_update_dlnActCoeffdT(); - double RTT = RT * T; for (size_t k = 0; k < m_kk; k++) { - hbar[k] -= RTT * dlnActCoeffdT_Scaled_[k]; + hbar[k] -= RT() * temperature() * dlnActCoeffdT_Scaled_[k]; } } @@ -400,7 +397,6 @@ void IonsFromNeutralVPSSTP::getPartialMolarEntropies(doublereal* sbar) const * Get the nondimensional standard state entropies */ getEntropy_R(sbar); - double T = temperature(); /* * Update the activity coefficients, This also update the * internally stored molalities. @@ -410,7 +406,7 @@ void IonsFromNeutralVPSSTP::getPartialMolarEntropies(doublereal* sbar) const for (size_t k = 0; k < m_kk; k++) { double xx = std::max(moleFractions_[k], SmallNumber); - sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k]; + sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - temperature() * dlnActCoeffdT_Scaled_[k]; } /* * dimensionalize it. diff --git a/src/thermo/LatticePhase.cpp b/src/thermo/LatticePhase.cpp index 2513b4979..438e5d60e 100644 --- a/src/thermo/LatticePhase.cpp +++ b/src/thermo/LatticePhase.cpp @@ -151,11 +151,10 @@ doublereal LatticePhase::logStandardConc(size_t k) const void LatticePhase::getChemPotentials(doublereal* mu) const { doublereal delta_p = m_Pcurrent - m_Pref; - doublereal RT = temperature() * GasConstant; const vector_fp& g_RT = gibbs_RT_ref(); for (size_t k = 0; k < m_kk; k++) { double xx = std::max(SmallNumber, moleFraction(k)); - mu[k] = RT * (g_RT[k] + log(xx)) + mu[k] = RT() * (g_RT[k] + log(xx)) + delta_p * m_speciesMolarVolume[k]; } } @@ -191,23 +190,22 @@ void LatticePhase::getPartialMolarVolumes(doublereal* vbar) const void LatticePhase::getStandardChemPotentials(doublereal* mu0) const { const vector_fp& gibbsrt = gibbs_RT_ref(); - scale(gibbsrt.begin(), gibbsrt.end(), mu0, _RT()); + scale(gibbsrt.begin(), gibbsrt.end(), mu0, RT()); } void LatticePhase::getPureGibbs(doublereal* gpure) const { const vector_fp& gibbsrt = gibbs_RT_ref(); doublereal delta_p = (m_Pcurrent - m_Pref); - double RT = GasConstant * temperature(); for (size_t k = 0; k < m_kk; k++) { - gpure[k] = RT * gibbsrt[k] + delta_p * m_speciesMolarVolume[k]; + gpure[k] = RT() * gibbsrt[k] + delta_p * m_speciesMolarVolume[k]; } } void LatticePhase::getEnthalpy_RT(doublereal* hrt) const { const vector_fp& _h = enthalpy_RT_ref(); - doublereal delta_prt = ((m_Pcurrent - m_Pref) / (GasConstant * temperature())); + doublereal delta_prt = (m_Pcurrent - m_Pref) / RT(); for (size_t k = 0; k < m_kk; k++) { hrt[k] = _h[k] + delta_prt * m_speciesMolarVolume[k]; } @@ -222,7 +220,7 @@ void LatticePhase::getEntropy_R(doublereal* sr) const void LatticePhase::getGibbs_RT(doublereal* grt) const { const vector_fp& gibbsrt = gibbs_RT_ref(); - doublereal delta_prt = (m_Pcurrent - m_Pref) / _RT(); + doublereal delta_prt = (m_Pcurrent - m_Pref) / RT(); for (size_t k = 0; k < m_kk; k++) { grt[k] = gibbsrt[k] + delta_prt * m_speciesMolarVolume[k]; } diff --git a/src/thermo/MargulesVPSSTP.cpp b/src/thermo/MargulesVPSSTP.cpp index 9165dc2f0..996d608b5 100644 --- a/src/thermo/MargulesVPSSTP.cpp +++ b/src/thermo/MargulesVPSSTP.cpp @@ -126,10 +126,9 @@ void MargulesVPSSTP::getChemPotentials(doublereal* mu) const * Update the activity coefficients */ s_update_lnActCoeff(); - doublereal RT = GasConstant * temperature(); for (size_t k = 0; k < m_kk; k++) { double xx = std::max(moleFractions_[k], SmallNumber); - mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]); + mu[k] += RT() * (log(xx) + lnActCoeff_Scaled_[k]); } } @@ -183,10 +182,8 @@ void MargulesVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const /* * dimensionalize it. */ - double T = temperature(); - double RT = GasConstant * T; for (size_t k = 0; k < m_kk; k++) { - hbar[k] *= RT; + hbar[k] *= RT(); } /* * Update the activity coefficients, This also update the @@ -194,9 +191,8 @@ void MargulesVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const */ s_update_lnActCoeff(); s_update_dlnActCoeff_dT(); - double RTT = RT * T; for (size_t k = 0; k < m_kk; k++) { - hbar[k] -= RTT * dlnActCoeffdT_Scaled_[k]; + hbar[k] -= RT() * temperature() * dlnActCoeffdT_Scaled_[k]; } } @@ -349,13 +345,12 @@ void MargulesVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_) void MargulesVPSSTP::s_update_lnActCoeff() const { double T = temperature(); - double invRT = 1.0 / (GasConstant*T); lnActCoeff_Scaled_.assign(m_kk, 0.0); for (size_t i = 0; i < numBinaryInteractions_; i++) { size_t iA = m_pSpecies_A_ij[i]; size_t iB = m_pSpecies_B_ij[i]; - double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) * invRT; - double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) * invRT; + double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT(); + double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT(); double XA = moleFractions_[iA]; double XB = moleFractions_[iB]; const doublereal XAXB = XA * XB; @@ -418,7 +413,6 @@ void MargulesVPSSTP::getdlnActCoeffds(const doublereal dTds, const doublereal* c doublereal* dlnActCoeffds) const { double T = temperature(); - double RT = GasConstant*T; s_update_dlnActCoeff_dT(); for (size_t iK = 0; iK < m_kk; iK++) { dlnActCoeffds[iK] = 0.0; @@ -431,8 +425,8 @@ void MargulesVPSSTP::getdlnActCoeffds(const doublereal dTds, const doublereal* c double XB = moleFractions_[iB]; double dXA = dXds[iA]; double dXB = dXds[iB]; - double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; - double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT(); + double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT(); const doublereal g02g1XB = g0 + 2*g1*XB; const doublereal g2XAdXB = 2*g1*XA*dXB; const doublereal all = (-XB * dXA - XA *dXB) * g02g1XB - XB *g2XAdXB; @@ -447,7 +441,6 @@ void MargulesVPSSTP::getdlnActCoeffds(const doublereal dTds, const doublereal* c void MargulesVPSSTP::s_update_dlnActCoeff_dlnN_diag() const { double T = temperature(); - double RT = GasConstant*T; dlnActCoeffdlnN_diag_.assign(m_kk, 0.0); for (size_t iK = 0; iK < m_kk; iK++) { @@ -468,8 +461,8 @@ void MargulesVPSSTP::s_update_dlnActCoeff_dlnN_diag() const double XA = moleFractions_[iA]; double XB = moleFractions_[iB]; - double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; - double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT(); + double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT(); dlnActCoeffdlnN_diag_[iK] += 2*(delBK-XB)*(g0*(delAK-XA)+g1*(2*(delAK-XA)*XB+XA*(delBK-XB))); } @@ -480,7 +473,6 @@ void MargulesVPSSTP::s_update_dlnActCoeff_dlnN_diag() const void MargulesVPSSTP::s_update_dlnActCoeff_dlnN() const { double T = temperature(); - double RT = GasConstant*T; dlnActCoeffdlnN_.zero(); /* @@ -509,8 +501,8 @@ void MargulesVPSSTP::s_update_dlnActCoeff_dlnN() const double XA = moleFractions_[iA]; double XB = moleFractions_[iB]; - double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; - double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT(); + double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT(); dlnActCoeffdlnN_(iK,iM) += g0*((delAM-XA)*(delBK-XB)+(delAK-XA)*(delBM-XB)); dlnActCoeffdlnN_(iK,iM) += 2*g1*((delAM-XA)*(delBK-XB)*XB+(delAK-XA)*(delBM-XB)*XB+(delBM-XB)*(delBK-XB)*XA); } @@ -523,7 +515,6 @@ void MargulesVPSSTP::s_update_dlnActCoeff_dlnX_diag() const { doublereal T = temperature(); dlnActCoeffdlnX_diag_.assign(m_kk, 0.0); - doublereal RT = GasConstant * T; for (size_t i = 0; i < numBinaryInteractions_; i++) { size_t iA = m_pSpecies_A_ij[i]; @@ -532,8 +523,8 @@ void MargulesVPSSTP::s_update_dlnActCoeff_dlnX_diag() const doublereal XA = moleFractions_[iA]; doublereal XB = moleFractions_[iB]; - doublereal g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; - doublereal g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + doublereal g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT(); + doublereal g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT(); dlnActCoeffdlnX_diag_[iA] += XA*XB*(2*g1*-2*g0-6*g1*XB); dlnActCoeffdlnX_diag_[iB] += XA*XB*(2*g1*-2*g0-6*g1*XB); diff --git a/src/thermo/MaskellSolidSolnPhase.cpp b/src/thermo/MaskellSolidSolnPhase.cpp index 996f14f7d..485ad2bb5 100644 --- a/src/thermo/MaskellSolidSolnPhase.cpp +++ b/src/thermo/MaskellSolidSolnPhase.cpp @@ -166,21 +166,19 @@ void MaskellSolidSolnPhase::getChemPotentials(doublereal* mu) const const doublereal r = moleFraction(product_species_index); const doublereal pval = p(r); const doublereal rfm = r * fm(r); - const doublereal RT = GasConstant * temperature(); const doublereal DgbarDr = pval * h_mixing + GasConstant * temperature() * std::log( (std::pow(1 - rfm, pval) * std::pow(rfm, pval) * std::pow(r - rfm, 1 - pval) * r) / (std::pow(1 - r - rfm, 1 + pval) * (1 - r)) ); - mu[product_species_index] = RT * m_g0_RT[product_species_index] + DgbarDr; - mu[reactant_species_index] = RT * m_g0_RT[reactant_species_index] - DgbarDr; + mu[product_species_index] = RT() * m_g0_RT[product_species_index] + DgbarDr; + mu[reactant_species_index] = RT() * m_g0_RT[reactant_species_index] - DgbarDr; } void MaskellSolidSolnPhase::getChemPotentials_RT(doublereal* mu) const { - const doublereal invRT = 1.0 / (GasConstant * temperature()); getChemPotentials(mu); for (size_t sp=0; sp < m_kk; ++sp) { - mu[sp] *= invRT; + mu[sp] *= 1.0 / RT(); } } diff --git a/src/thermo/MetalSHEelectrons.cpp b/src/thermo/MetalSHEelectrons.cpp index fd33ddbd0..64e1d2d62 100644 --- a/src/thermo/MetalSHEelectrons.cpp +++ b/src/thermo/MetalSHEelectrons.cpp @@ -195,8 +195,7 @@ void MetalSHEelectrons::getIntEnergy_RT(doublereal* urt) const void MetalSHEelectrons::getIntEnergy_RT_ref(doublereal* urt) const { _updateThermo(); - doublereal RT = GasConstant * temperature(); - urt[0] = m_h0_RT[0] - m_p0 / molarDensity() / RT; + urt[0] = m_h0_RT[0] - m_p0 / molarDensity() / RT(); } /* diff --git a/src/thermo/MineralEQ3.cpp b/src/thermo/MineralEQ3.cpp index 8889bbdb7..ebef9986a 100644 --- a/src/thermo/MineralEQ3.cpp +++ b/src/thermo/MineralEQ3.cpp @@ -155,9 +155,8 @@ void MineralEQ3::getStandardChemPotentials(doublereal* mu0) const void MineralEQ3::getEnthalpy_RT(doublereal* hrt) const { getEnthalpy_RT_ref(hrt); - doublereal RT = GasConstant * temperature(); doublereal presCorrect = (m_press - m_p0) / molarDensity(); - hrt[0] += presCorrect / RT; + hrt[0] += presCorrect / RT(); } void MineralEQ3::getEntropy_R(doublereal* sr) const @@ -180,8 +179,7 @@ void MineralEQ3::getCp_R(doublereal* cpr) const void MineralEQ3::getIntEnergy_RT(doublereal* urt) const { _updateThermo(); - doublereal RT = GasConstant * temperature(); - urt[0] = m_h0_RT[0] - m_p0 / molarDensity() / RT; + urt[0] = m_h0_RT[0] - m_p0 / molarDensity() / RT(); } /* @@ -191,8 +189,7 @@ void MineralEQ3::getIntEnergy_RT(doublereal* urt) const void MineralEQ3::getIntEnergy_RT_ref(doublereal* urt) const { _updateThermo(); - doublereal RT = GasConstant * temperature(); - urt[0] = m_h0_RT[0] - m_p0 / molarDensity() / RT; + urt[0] = m_h0_RT[0] - m_p0 / molarDensity() / RT(); } /* diff --git a/src/thermo/MixedSolventElectrolyte.cpp b/src/thermo/MixedSolventElectrolyte.cpp index ecf66ef69..ba415c2e3 100644 --- a/src/thermo/MixedSolventElectrolyte.cpp +++ b/src/thermo/MixedSolventElectrolyte.cpp @@ -130,10 +130,9 @@ void MixedSolventElectrolyte::getChemPotentials(doublereal* mu) const * Update the activity coefficients */ s_update_lnActCoeff(); - doublereal RT = GasConstant * temperature(); for (size_t k = 0; k < m_kk; k++) { double xx = std::max(moleFractions_[k], SmallNumber); - mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]); + mu[k] += RT() * (log(xx) + lnActCoeff_Scaled_[k]); } } @@ -184,10 +183,8 @@ void MixedSolventElectrolyte::getPartialMolarEnthalpies(doublereal* hbar) const /* * dimensionalize it. */ - double T = temperature(); - double RT = GasConstant * T; for (size_t k = 0; k < m_kk; k++) { - hbar[k] *= RT; + hbar[k] *= RT(); } /* * Update the activity coefficients, This also update the @@ -195,9 +192,8 @@ void MixedSolventElectrolyte::getPartialMolarEnthalpies(doublereal* hbar) const */ s_update_lnActCoeff(); s_update_dlnActCoeff_dT(); - double RTT = RT * T; for (size_t k = 0; k < m_kk; k++) { - hbar[k] -= RTT * dlnActCoeffdT_Scaled_[k]; + hbar[k] -= RT() * temperature() * dlnActCoeffdT_Scaled_[k]; } } @@ -351,7 +347,6 @@ void MixedSolventElectrolyte::initThermoXML(XML_Node& phaseNode, const std::stri void MixedSolventElectrolyte::s_update_lnActCoeff() const { double T = temperature(); - double RT = GasConstant*T; lnActCoeff_Scaled_.assign(m_kk, 0.0); for (size_t iK = 0; iK < m_kk; iK++) { for (size_t i = 0; i < numBinaryInteractions_; i++) { @@ -366,8 +361,8 @@ void MixedSolventElectrolyte::s_update_lnActCoeff() const } double XA = moleFractions_[iA]; double XB = moleFractions_[iB]; - double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; - double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT(); + double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT(); lnActCoeff_Scaled_[iK] += (delAK * XB + XA * delBK - XA * XB) * (g0 + g1 * XB) + XA * XB * (delBK - XB) * g1; } } @@ -421,7 +416,6 @@ void MixedSolventElectrolyte::getdlnActCoeffds(const doublereal dTds, const doub doublereal* dlnActCoeffds) const { double T = temperature(); - double RT = GasConstant*T; s_update_dlnActCoeff_dT(); for (size_t iK = 0; iK < m_kk; iK++) { @@ -442,8 +436,8 @@ void MixedSolventElectrolyte::getdlnActCoeffds(const doublereal dTds, const doub double XB = moleFractions_[iB]; double dXA = dXds[iA]; double dXB = dXds[iB]; - double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; - double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT(); + double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT(); dlnActCoeffds[iK] += ((delBK-XB)*dXA + (delAK-XA)*dXB)*(g0+2*g1*XB) + (delBK-XB)*2*g1*XA*dXB + dlnActCoeffdT_Scaled_[iK]*dTds; } @@ -453,7 +447,6 @@ void MixedSolventElectrolyte::getdlnActCoeffds(const doublereal dTds, const doub void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN_diag() const { double T = temperature(); - double RT = GasConstant*T; dlnActCoeffdlnN_diag_.assign(m_kk, 0); for (size_t iK = 0; iK < m_kk; iK++) { @@ -472,8 +465,8 @@ void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN_diag() const double XA = moleFractions_[iA]; double XB = moleFractions_[iB]; - double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; - double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT(); + double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT(); dlnActCoeffdlnN_diag_[iK] += 2*(delBK-XB)*(g0*(delAK-XA)+g1*(2*(delAK-XA)*XB+XA*(delBK-XB))); } @@ -484,7 +477,6 @@ void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN_diag() const void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN() const { double T = temperature(); - double RT = GasConstant*T; dlnActCoeffdlnN_.zero(); /* @@ -513,8 +505,8 @@ void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN() const double XA = moleFractions_[iA]; double XB = moleFractions_[iB]; - double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; - double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT(); + double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT(); dlnActCoeffdlnN_(iK,iM) += g0*((delAM-XA)*(delBK-XB)+(delAK-XA)*(delBM-XB)); dlnActCoeffdlnN_(iK,iM) += 2*g1*((delAM-XA)*(delBK-XB)*XB+(delAK-XA)*(delBM-XB)*XB+(delBM-XB)*(delBK-XB)*XA); } @@ -527,15 +519,14 @@ void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnX_diag() const { doublereal T = temperature(); dlnActCoeffdlnX_diag_.assign(m_kk, 0); - doublereal RT = GasConstant * T; for (size_t i = 0; i < numBinaryInteractions_; i++) { size_t iA = m_pSpecies_A_ij[i]; size_t iB = m_pSpecies_B_ij[i]; double XA = moleFractions_[iA]; double XB = moleFractions_[iB]; - double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT; - double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; + double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT(); + double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT(); dlnActCoeffdlnX_diag_[iA] += XA*XB*(2*g1*-2*g0-6*g1*XB); dlnActCoeffdlnX_diag_[iB] += XA*XB*(2*g1*-2*g0-6*g1*XB); } diff --git a/src/thermo/MixtureFugacityTP.cpp b/src/thermo/MixtureFugacityTP.cpp index 19428e2d2..42a50f151 100644 --- a/src/thermo/MixtureFugacityTP.cpp +++ b/src/thermo/MixtureFugacityTP.cpp @@ -95,9 +95,8 @@ int MixtureFugacityTP::reportSolnBranchActual() const void MixtureFugacityTP::getChemPotentials_RT(doublereal* muRT) const { getChemPotentials(muRT); - doublereal invRT = 1.0 / _RT(); for (size_t k = 0; k < m_kk; k++) { - muRT[k] *= invRT; + muRT[k] *= 1.0 / RT(); } } @@ -109,10 +108,9 @@ void MixtureFugacityTP::getStandardChemPotentials(doublereal* g) const { _updateReferenceStateThermo(); copy(m_g0_RT.begin(), m_g0_RT.end(), g); - doublereal RT = _RT(); double tmp = log(pressure() /m_spthermo->refPressure()); for (size_t k = 0; k < m_kk; k++) { - g[k] = RT * (g[k] + tmp); + g[k] = RT() * (g[k] + tmp); } } @@ -150,8 +148,8 @@ void MixtureFugacityTP::getGibbs_RT(doublereal* grt) const void MixtureFugacityTP::getPureGibbs(doublereal* g) const { _updateReferenceStateThermo(); - scale(m_g0_RT.begin(), m_g0_RT.end(), g, _RT()); - double tmp = log(pressure() /m_spthermo->refPressure()) * _RT(); + scale(m_g0_RT.begin(), m_g0_RT.end(), g, RT()); + double tmp = log(pressure() /m_spthermo->refPressure()) * RT(); for (size_t k = 0; k < m_kk; k++) { g[k] += tmp; } @@ -175,9 +173,8 @@ void MixtureFugacityTP::getCp_R(doublereal* cpr) const void MixtureFugacityTP::getStandardVolumes(doublereal* vol) const { _updateReferenceStateThermo(); - doublereal v0 = _RT() / pressure(); for (size_t i = 0; i < m_kk; i++) { - vol[i]= v0; + vol[i] = RT() / pressure(); } } @@ -201,7 +198,7 @@ void MixtureFugacityTP::getGibbs_RT_ref(doublereal* grt) const void MixtureFugacityTP::getGibbs_ref(doublereal* g) const { const vector_fp& gibbsrt = gibbs_RT_ref(); - scale(gibbsrt.begin(), gibbsrt.end(), g, _RT()); + scale(gibbsrt.begin(), gibbsrt.end(), g, RT()); } const vector_fp& MixtureFugacityTP::gibbs_RT_ref() const @@ -225,9 +222,8 @@ void MixtureFugacityTP::getCp_R_ref(doublereal* cpr) const void MixtureFugacityTP::getStandardVolumes_ref(doublereal* vol) const { _updateReferenceStateThermo(); - doublereal v0 = _RT() / refPressure(); for (size_t i = 0; i < m_kk; i++) { - vol[i]= v0; + vol[i]= RT() / refPressure(); } } @@ -436,7 +432,7 @@ void MixtureFugacityTP::initThermoXML(XML_Node& phaseNode, const std::string& id doublereal MixtureFugacityTP::z() const { - return pressure() * meanMolecularWeight() / (density() * _RT()); + return pressure() * meanMolecularWeight() / (density() * RT()); } doublereal MixtureFugacityTP::sresid() const @@ -657,7 +653,7 @@ int MixtureFugacityTP::corr0(doublereal TKelvin, doublereal pres, doublereal& de } else { densLiqGuess = densLiq; setState_TR(TKelvin, densLiq); - liqGRT = gibbs_mole() / _RT(); + liqGRT = gibbs_mole() / RT(); } doublereal densGas = densityCalc(TKelvin, pres, FLUID_GAS, densGasGuess); @@ -671,7 +667,7 @@ int MixtureFugacityTP::corr0(doublereal TKelvin, doublereal pres, doublereal& de } else { densGasGuess = densGas; setState_TR(TKelvin, densGas); - gasGRT = gibbs_mole() / _RT(); + gasGRT = gibbs_mole() / RT(); } return retn; } diff --git a/src/thermo/MolarityIonicVPSSTP.cpp b/src/thermo/MolarityIonicVPSSTP.cpp index 01549169c..754c60cb3 100644 --- a/src/thermo/MolarityIonicVPSSTP.cpp +++ b/src/thermo/MolarityIonicVPSSTP.cpp @@ -120,10 +120,9 @@ void MolarityIonicVPSSTP::getChemPotentials(doublereal* mu) const * Update the activity coefficients */ s_update_lnActCoeff(); - doublereal RT = GasConstant * temperature(); for (size_t k = 0; k < m_kk; k++) { double xx = std::max(moleFractions_[k], SmallNumber); - mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]); + mu[k] += RT() * (log(xx) + lnActCoeff_Scaled_[k]); } } diff --git a/src/thermo/PhaseCombo_Interaction.cpp b/src/thermo/PhaseCombo_Interaction.cpp index 90644a8da..bd5234eea 100644 --- a/src/thermo/PhaseCombo_Interaction.cpp +++ b/src/thermo/PhaseCombo_Interaction.cpp @@ -135,10 +135,9 @@ void PhaseCombo_Interaction::getChemPotentials(doublereal* mu) const */ s_update_lnActCoeff(); - doublereal RT = GasConstant * temperature(); for (size_t k = 0; k < m_kk; k++) { double xx = std::max(moleFractions_[k], SmallNumber); - mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]); + mu[k] += RT() * (log(xx) + lnActCoeff_Scaled_[k]); } } diff --git a/src/thermo/PureFluidPhase.cpp b/src/thermo/PureFluidPhase.cpp index 014fff7c7..bd326de18 100644 --- a/src/thermo/PureFluidPhase.cpp +++ b/src/thermo/PureFluidPhase.cpp @@ -219,7 +219,7 @@ void PureFluidPhase::getStandardChemPotentials(doublereal* mu) const void PureFluidPhase::getEnthalpy_RT(doublereal* hrt) const { - hrt[0] = enthalpy_mole() / _RT(); + hrt[0] = enthalpy_mole() / RT(); } void PureFluidPhase::getEntropy_R(doublereal* sr) const @@ -229,7 +229,7 @@ void PureFluidPhase::getEntropy_R(doublereal* sr) const void PureFluidPhase::getGibbs_RT(doublereal* grt) const { - grt[0] = gibbs_mole() / _RT(); + grt[0] = gibbs_mole() / RT(); } void PureFluidPhase::getEnthalpy_RT_ref(doublereal* hrt) const diff --git a/src/thermo/RedlichKisterVPSSTP.cpp b/src/thermo/RedlichKisterVPSSTP.cpp index b880d6102..73408ecbe 100644 --- a/src/thermo/RedlichKisterVPSSTP.cpp +++ b/src/thermo/RedlichKisterVPSSTP.cpp @@ -125,10 +125,9 @@ void RedlichKisterVPSSTP::getChemPotentials(doublereal* mu) const */ s_update_lnActCoeff(); - doublereal RT = GasConstant * temperature(); for (size_t k = 0; k < m_kk; k++) { double xx = std::max(moleFractions_[k], SmallNumber); - mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]); + mu[k] += RT() * (log(xx) + lnActCoeff_Scaled_[k]); } } diff --git a/src/thermo/RedlichKwongMFTP.cpp b/src/thermo/RedlichKwongMFTP.cpp index 136babef9..2768731a1 100644 --- a/src/thermo/RedlichKwongMFTP.cpp +++ b/src/thermo/RedlichKwongMFTP.cpp @@ -150,7 +150,7 @@ int RedlichKwongMFTP::eosType() const doublereal RedlichKwongMFTP::enthalpy_mole() const { _updateReferenceStateThermo(); - doublereal h_ideal = _RT() * mean_X(m_h0_RT); + doublereal h_ideal = RT() * mean_X(m_h0_RT); doublereal h_nonideal = hresid(); return h_ideal + h_nonideal; } @@ -278,10 +278,8 @@ doublereal RedlichKwongMFTP::standardConcentration(size_t k) const void RedlichKwongMFTP::getActivityCoefficients(doublereal* ac) const { - doublereal TKelvin = temperature(); - doublereal rt = TKelvin * GasConstant; doublereal mv = molarVolume(); - doublereal sqt = sqrt(TKelvin); + doublereal sqt = sqrt(temperature()); doublereal vpb = mv + m_b_current; doublereal vmb = mv - m_b_current; @@ -295,16 +293,16 @@ void RedlichKwongMFTP::getActivityCoefficients(doublereal* ac) const doublereal pres = pressure(); for (size_t k = 0; k < m_kk; k++) { - ac[k] = (- rt * log(pres * mv / rt) - + rt * log(mv / vmb) - + rt * b_vec_Curr_[k] / vmb + ac[k] = (- RT() * log(pres * mv / RT()) + + RT() * log(mv / vmb) + + RT() * b_vec_Curr_[k] / vmb - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv) + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb) ); } for (size_t k = 0; k < m_kk; k++) { - ac[k] = exp(ac[k]/rt); + ac[k] = exp(ac[k]/RT()); } } @@ -315,24 +313,21 @@ void RedlichKwongMFTP::getActivityCoefficients(doublereal* ac) const void RedlichKwongMFTP::getChemPotentials_RT(doublereal* muRT) const { getChemPotentials(muRT); - doublereal invRT = 1.0 / _RT(); for (size_t k = 0; k < m_kk; k++) { - muRT[k] *= invRT; + muRT[k] *= 1.0 / RT(); } } void RedlichKwongMFTP::getChemPotentials(doublereal* mu) const { getGibbs_ref(mu); - doublereal rt = temperature() * GasConstant; for (size_t k = 0; k < m_kk; k++) { double xx = std::max(SmallNumber, moleFraction(k)); - mu[k] += rt*(log(xx)); + mu[k] += RT()*(log(xx)); } - doublereal TKelvin = temperature(); doublereal mv = molarVolume(); - doublereal sqt = sqrt(TKelvin); + doublereal sqt = sqrt(temperature()); doublereal vpb = mv + m_b_current; doublereal vmb = mv - m_b_current; @@ -347,9 +342,9 @@ void RedlichKwongMFTP::getChemPotentials(doublereal* mu) const doublereal refP = refPressure(); for (size_t k = 0; k < m_kk; k++) { - mu[k] += (rt * log(pres/refP) - rt * log(pres * mv / rt) - + rt * log(mv / vmb) - + rt * b_vec_Curr_[k] / vmb + mu[k] += (RT() * log(pres/refP) - RT() * log(pres * mv / RT()) + + RT() * log(mv / vmb) + + RT() * b_vec_Curr_[k] / vmb - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv) + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb) @@ -363,8 +358,7 @@ void RedlichKwongMFTP::getPartialMolarEnthalpies(doublereal* hbar) const * First we get the reference state contributions */ getEnthalpy_RT_ref(hbar); - doublereal rt = GasConstant * temperature(); - scale(hbar, hbar+m_kk, hbar, rt); + scale(hbar, hbar+m_kk, hbar, RT()); /* * We calculate dpdni_ @@ -382,7 +376,7 @@ void RedlichKwongMFTP::getPartialMolarEnthalpies(doublereal* hbar) const } } for (size_t k = 0; k < m_kk; k++) { - dpdni_[k] = rt/vmb + rt * b_vec_Curr_[k] / (vmb * vmb) - 2.0 * m_pp[k] / (sqt * mv * vpb) + dpdni_[k] = RT()/vmb + RT() * b_vec_Curr_[k] / (vmb * vmb) - 2.0 * m_pp[k] / (sqt * mv * vpb) + m_a_current * b_vec_Curr_[k]/(sqt * mv * vpb * vpb); } doublereal dadt = da_dt(); @@ -399,7 +393,7 @@ void RedlichKwongMFTP::getPartialMolarEnthalpies(doublereal* hbar) const pressureDerivatives(); doublereal fac2 = mv + TKelvin * dpdT_ / dpdV_; for (size_t k = 0; k < m_kk; k++) { - double hE_v = (mv * dpdni_[k] - rt - b_vec_Curr_[k]/ (m_b_current * m_b_current * sqt) * log(vpb/mv)*fac + double hE_v = (mv * dpdni_[k] - RT() - b_vec_Curr_[k]/ (m_b_current * m_b_current * sqt) * log(vpb/mv)*fac + 1.0 / (m_b_current * sqt) * log(vpb/mv) * m_tmpV[k] + b_vec_Curr_[k] / vpb / (m_b_current * sqt) * fac); hbar[k] = hbar[k] + hE_v; @@ -410,8 +404,7 @@ void RedlichKwongMFTP::getPartialMolarEnthalpies(doublereal* hbar) const void RedlichKwongMFTP::getPartialMolarEntropies(doublereal* sbar) const { getEntropy_R_ref(sbar); - doublereal r = GasConstant; - scale(sbar, sbar+m_kk, sbar, r); + scale(sbar, sbar+m_kk, sbar, GasConstant); doublereal TKelvin = temperature(); doublereal sqt = sqrt(TKelvin); doublereal mv = molarVolume(); @@ -419,7 +412,7 @@ void RedlichKwongMFTP::getPartialMolarEntropies(doublereal* sbar) const for (size_t k = 0; k < m_kk; k++) { doublereal xx = std::max(SmallNumber, moleFraction(k)); - sbar[k] += r * (- log(xx)); + sbar[k] += GasConstant * (- log(xx)); } for (size_t k = 0; k < m_kk; k++) { m_pp[k] = 0.0; @@ -462,15 +455,13 @@ void RedlichKwongMFTP::getPartialMolarEntropies(doublereal* sbar) const void RedlichKwongMFTP::getPartialMolarIntEnergies(doublereal* ubar) const { getIntEnergy_RT(ubar); - doublereal rt = GasConstant * temperature(); - scale(ubar, ubar+m_kk, ubar, rt); + scale(ubar, ubar+m_kk, ubar, RT()); } void RedlichKwongMFTP::getPartialMolarCp(doublereal* cpbar) const { getCp_R(cpbar); - doublereal r = GasConstant; - scale(cpbar, cpbar+m_kk, cpbar, r); + scale(cpbar, cpbar+m_kk, cpbar, GasConstant); } void RedlichKwongMFTP::getPartialMolarVolumes(doublereal* vbar) const @@ -490,19 +481,17 @@ void RedlichKwongMFTP::getPartialMolarVolumes(doublereal* vbar) const } } - doublereal TKelvin = temperature(); - doublereal sqt = sqrt(TKelvin); + doublereal sqt = sqrt(temperature()); doublereal mv = molarVolume(); - doublereal rt = GasConstant * TKelvin; doublereal vmb = mv - m_b_current; doublereal vpb = mv + m_b_current; for (size_t k = 0; k < m_kk; k++) { - doublereal num = (rt + rt * m_b_current/ vmb + rt * b_vec_Curr_[k] / vmb - + rt * m_b_current * b_vec_Curr_[k] /(vmb * vmb) + doublereal num = (RT() + RT() * m_b_current/ vmb + RT() * b_vec_Curr_[k] / vmb + + RT() * m_b_current * b_vec_Curr_[k] /(vmb * vmb) - 2.0 * m_pp[k] / (sqt * vpb) + m_a_current * b_vec_Curr_[k] / (sqt * vpb * vpb) ); - doublereal denom = (m_Pcurrent + rt * m_b_current/(vmb * vmb) - m_a_current / (sqt * vpb * vpb) + doublereal denom = (m_Pcurrent + RT() * m_b_current/(vmb * vmb) - m_a_current / (sqt * vpb * vpb) ); vbar[k] = num / denom; } diff --git a/src/thermo/SurfPhase.cpp b/src/thermo/SurfPhase.cpp index 6df6468d5..18533e37f 100644 --- a/src/thermo/SurfPhase.cpp +++ b/src/thermo/SurfPhase.cpp @@ -126,9 +126,8 @@ doublereal SurfPhase::cv_mole() const void SurfPhase::getPartialMolarEnthalpies(doublereal* hbar) const { getEnthalpy_RT(hbar); - doublereal rt = GasConstant * temperature(); for (size_t k = 0; k < m_kk; k++) { - hbar[k] *= rt; + hbar[k] *= RT(); } } diff --git a/src/thermo/ThermoPhase.cpp b/src/thermo/ThermoPhase.cpp index ee596c8ce..1ae9dbcdb 100644 --- a/src/thermo/ThermoPhase.cpp +++ b/src/thermo/ThermoPhase.cpp @@ -963,7 +963,6 @@ std::string ThermoPhase::report(bool show_thermo, doublereal threshold) const getMoleFractions(&x[0]); getMassFractions(&y[0]); getChemPotentials(&mu[0]); - doublereal rt = GasConstant * temperature(); int nMinor = 0; doublereal xMinor = 0.0; doublereal yMinor = 0.0; @@ -978,7 +977,7 @@ std::string ThermoPhase::report(bool show_thermo, doublereal threshold) const if (x[k] >= threshold) { if (x[k] > SmallNumber) { sprintf(p, "%18s %12.6g %12.6g %12.6g\n", - speciesName(k).c_str(), x[k], y[k], mu[k]/rt); + speciesName(k).c_str(), x[k], y[k], mu[k]/RT()); } else { sprintf(p, "%18s %12.6g %12.6g \n", speciesName(k).c_str(), x[k], y[k]); diff --git a/src/thermo/VPStandardStateTP.cpp b/src/thermo/VPStandardStateTP.cpp index 4c293285a..b1916ce29 100644 --- a/src/thermo/VPStandardStateTP.cpp +++ b/src/thermo/VPStandardStateTP.cpp @@ -122,9 +122,8 @@ int VPStandardStateTP::standardStateConvention() const void VPStandardStateTP::getChemPotentials_RT(doublereal* muRT) const { getChemPotentials(muRT); - doublereal invRT = 1.0 / _RT(); for (size_t k = 0; k < m_kk; k++) { - muRT[k] *= invRT; + muRT[k] *= 1.0 / RT(); } } @@ -134,9 +133,8 @@ void VPStandardStateTP::getChemPotentials_RT(doublereal* muRT) const void VPStandardStateTP::getStandardChemPotentials(doublereal* g) const { getGibbs_RT(g); - doublereal RT = _RT(); for (size_t k = 0; k < m_kk; k++) { - g[k] *= RT; + g[k] *= RT(); } } diff --git a/src/thermo/WaterSSTP.cpp b/src/thermo/WaterSSTP.cpp index f6fc66c74..3388799fb 100644 --- a/src/thermo/WaterSSTP.cpp +++ b/src/thermo/WaterSSTP.cpp @@ -278,9 +278,8 @@ void WaterSSTP::getGibbs_RT_ref(doublereal* grt) const void WaterSSTP::getGibbs_ref(doublereal* g) const { getGibbs_RT_ref(g); - doublereal rt = _RT(); for (size_t k = 0; k < m_kk; k++) { - g[k] *= rt; + g[k] *= RT(); } } diff --git a/test/kinetics/rates.cpp b/test/kinetics/rates.cpp index 0ed4c04f5..cb944aa53 100644 --- a/test/kinetics/rates.cpp +++ b/test/kinetics/rates.cpp @@ -137,7 +137,7 @@ TEST_F(FracCoeffTest, EquilibriumConstants) double deltaG0_1 = mu0[kH2O] - 0.7 * mu0[kH2] - 0.6 * mu0[kOH] - 0.2 * mu0[kO2]; double pRef = therm.refPressure(); - double RT = GasConstant * therm.temperature(); + double RT = therm.RT(); // Net stoichiometric coefficients are 1.2 and -0.5 EXPECT_NEAR(exp(-deltaG0_0/RT) * pow(pRef/RT, 1.2), Kc[0], 1e-13 * Kc[0]);