research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Consistently use ThermoPhase::RT()

Browse source
This commit is contained in:
Ray Speth 2015-08-01 23:23:47 -04:00
parent 9cdfd12172
commit ceefc5ecb0
32 changed files with 160 additions and 257 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
include/cantera/thermo/ConstDensityThermo.h
View file

@ -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

14
include/cantera/thermo/ThermoPhase.h
View file

@ -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;
}

6
src/equil/ChemEquil.cpp
View file

@ -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) {

7
src/kinetics/AqueousKinetics.cpp
View file

@ -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

6
src/kinetics/BulkKinetics.cpp
View file

@ -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);

25
src/kinetics/InterfaceKinetics.cpp
View file

@ -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<size_t>(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

3
src/thermo/ConstDensityThermo.cpp
View file

@ -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;
}
}

33
src/thermo/DebyeHuckel.cpp
View file

@ -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]);
}
}
}

12
src/thermo/FixedChemPotSSTP.cpp
View file

@ -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

36
src/thermo/HMWSoln.cpp
View file

@ -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]);
}
}

20
src/thermo/IdealGasPhase.cpp
View file

@ -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;
}

12
src/thermo/IdealMolalSoln.cpp
View file

@ -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();
}
}

12
src/thermo/IdealSolidSolnPhase.cpp
View file

@ -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];
}

6
src/thermo/IdealSolnGasVPSS.cpp
View file

@ -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);
}
}

10
src/thermo/IonsFromNeutralVPSSTP.cpp
View file

@ -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.

12
src/thermo/LatticePhase.cpp
View file

@ -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];
}

35
src/thermo/MargulesVPSSTP.cpp
View file

@ -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);

8
src/thermo/MaskellSolidSolnPhase.cpp
View file

@ -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();
}
}

3
src/thermo/MetalSHEelectrons.cpp
View file

@ -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();
}
/*

9
src/thermo/MineralEQ3.cpp
View file

@ -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();
}
/*

35
src/thermo/MixedSolventElectrolyte.cpp
View file

@ -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);
}

24
src/thermo/MixtureFugacityTP.cpp
View file

@ -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;
}

3
src/thermo/MolarityIonicVPSSTP.cpp
View file

@ -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]);
}
}

3
src/thermo/PhaseCombo_Interaction.cpp
View file

@ -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]);
}
}

4
src/thermo/PureFluidPhase.cpp
View file

@ -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

3
src/thermo/RedlichKisterVPSSTP.cpp
View file

@ -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]);
}
}

57
src/thermo/RedlichKwongMFTP.cpp
View file

@ -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;
}

3
src/thermo/SurfPhase.cpp
View file

@ -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();
}
}

3
src/thermo/ThermoPhase.cpp
View file

@ -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]);

6
src/thermo/VPStandardStateTP.cpp
View file

@ -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();
}
}

3
src/thermo/WaterSSTP.cpp
View file

@ -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();
}
}

2
test/kinetics/rates.cpp
View file

@ -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]);