diff --git a/include/cantera/thermo/MolalityVPSSTP.h b/include/cantera/thermo/MolalityVPSSTP.h
index 9bd80ad31..ab5642f3f 100644
--- a/include/cantera/thermo/MolalityVPSSTP.h
+++ b/include/cantera/thermo/MolalityVPSSTP.h
@@ -225,10 +225,14 @@ public:
* molality.
*
* @param k the solvent index number
+ * @deprecated The solvent is always the first species in the phase. To be
+ * removed after Cantera 2.4.
*/
void setSolvent(size_t k);
//! Returns the solvent index.
+ //! @deprecated The solvent is always the first species in the phase. To be
+ //! removed after Cantera 2.4.
size_t solventIndex() const;
/**
@@ -563,11 +567,6 @@ private:
virtual size_t findCLMIndex() const;
protected:
-
- //! Index of the solvent. Currently the index of the solvent is hard-coded
- //! to the value 0
- size_t m_indexSolvent;
-
//! Scaling to be used for output of single-ion species activity
//! coefficients.
/*!
diff --git a/src/thermo/DebyeHuckel.cpp b/src/thermo/DebyeHuckel.cpp
index dc30bddf1..6945b9bc0 100644
--- a/src/thermo/DebyeHuckel.cpp
+++ b/src/thermo/DebyeHuckel.cpp
@@ -146,7 +146,7 @@ void DebyeHuckel::getActivityConcentrations(doublereal* c) const
doublereal DebyeHuckel::standardConcentration(size_t k) const
{
- double mvSolvent = m_speciesSize[m_indexSolvent];
+ double mvSolvent = m_speciesSize[0];
return 1.0 / mvSolvent;
}
@@ -157,14 +157,11 @@ void DebyeHuckel::getActivities(doublereal* ac) const
// Update the molality array, m_molalities(). This requires an update due to
// mole fractions
s_update_lnMolalityActCoeff();
- for (size_t k = 0; k < m_kk; k++) {
- if (k != m_indexSolvent) {
- ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal[k]);
- }
+ for (size_t k = 1; k < m_kk; k++) {
+ ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal[k]);
}
- double xmolSolvent = moleFraction(m_indexSolvent);
- ac[m_indexSolvent] =
- exp(m_lnActCoeffMolal[m_indexSolvent]) * xmolSolvent;
+ double xmolSolvent = moleFraction(0);
+ ac[0] = exp(m_lnActCoeffMolal[0]) * xmolSolvent;
}
void DebyeHuckel::getMolalityActivityCoefficients(doublereal* acMolality) const
@@ -191,16 +188,13 @@ void DebyeHuckel::getChemPotentials(doublereal* mu) const
// Update the activity coefficients. This also updates the internal molality
// array.
s_update_lnMolalityActCoeff();
- 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]);
- }
+ double xmolSolvent = moleFraction(0);
+ for (size_t k = 1; k < m_kk; k++) {
+ xx = std::max(m_molalities[k], SmallNumber);
+ mu[k] += RT() * (log(xx) + m_lnActCoeffMolal[k]);
}
xx = std::max(xmolSolvent, SmallNumber);
- mu[m_indexSolvent] +=
- RT() * (log(xx) + m_lnActCoeffMolal[m_indexSolvent]);
+ mu[0] += RT() * (log(xx) + m_lnActCoeffMolal[0]);
}
void DebyeHuckel::getPartialMolarEnthalpies(doublereal* hbar) const
@@ -246,15 +240,13 @@ void DebyeHuckel::getPartialMolarEntropies(doublereal* sbar) const
// First we will add in the obvious dependence on the T term out front of
// the log activity term
doublereal mm;
- for (size_t k = 0; k < m_kk; k++) {
- if (k != m_indexSolvent) {
- mm = std::max(SmallNumber, m_molalities[k]);
- sbar[k] -= GasConstant * (log(mm) + m_lnActCoeffMolal[k]);
- }
+ for (size_t k = 1; k < m_kk; k++) {
+ mm = std::max(SmallNumber, m_molalities[k]);
+ sbar[k] -= GasConstant * (log(mm) + m_lnActCoeffMolal[k]);
}
- double xmolSolvent = moleFraction(m_indexSolvent);
+ double xmolSolvent = moleFraction(0);
mm = std::max(SmallNumber, xmolSolvent);
- sbar[m_indexSolvent] -= GasConstant *(log(mm) + m_lnActCoeffMolal[m_indexSolvent]);
+ sbar[0] -= GasConstant *(log(mm) + m_lnActCoeffMolal[0]);
// Check to see whether activity coefficients are temperature dependent. If
// they are, then calculate the their temperature derivatives and add them
@@ -428,40 +420,6 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
setDebyeHuckelModel("Dilute_limit");
}
- // Reconcile the solvent name and index.
-
- // Get the Name of the Solvent:
- // solventName
- std::string solventName = "";
- if (thermoNode.hasChild("solvent")) {
- XML_Node& scNode = thermoNode.child("solvent");
- vector nameSolventa;
- getStringArray(scNode, nameSolventa);
- if (nameSolventa.size() != 1) {
- throw CanteraError("DebyeHuckel::initThermoXML",
- "badly formed solvent XML node");
- }
- solventName = nameSolventa[0];
- }
- for (size_t k = 0; k < m_kk; k++) {
- std::string sname = speciesName(k);
- if (solventName == sname) {
- m_indexSolvent = k;
- break;
- }
- }
- if (m_indexSolvent == npos) {
- cout << "DebyeHuckel::initThermoXML: Solvent Name not found"
- << endl;
- throw CanteraError("DebyeHuckel::initThermoXML",
- "Solvent name not found");
- }
- if (m_indexSolvent != 0) {
- throw CanteraError("DebyeHuckel::initThermoXML",
- "Solvent " + solventName +
- " should be first species");
- }
-
// Go get all of the coefficients and factors in the activityCoefficients
// XML block
XML_Node* acNodePtr = 0;
@@ -822,10 +780,8 @@ double DebyeHuckel::_lnactivityWaterHelgesonFixedForm() const
calcMolalities();
double oc = _osmoticCoeffHelgesonFixedForm();
double sum = 0.0;
- for (size_t k = 0; k < m_kk; k++) {
- if (k != m_indexSolvent) {
- sum += std::max(m_molalities[k], 0.0);
- }
+ for (size_t k = 1; k < m_kk; k++) {
+ sum += std::max(m_molalities[k], 0.0);
}
if (sum > 2.0 * m_maxIionicStrength) {
sum = 2.0 * m_maxIionicStrength;
@@ -876,7 +832,7 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
m_A_Debye = A_Debye_TP();
// Calculate a safe value for the mole fraction of the solvent
- double xmolSolvent = moleFraction(m_indexSolvent);
+ double xmolSolvent = moleFraction(0);
xmolSolvent = std::max(8.689E-3, xmolSolvent);
int est;
@@ -920,7 +876,7 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
tmp = 0.0;
if (denomTmp > 0.0) {
for (size_t k = 0; k < m_kk; k++) {
- if (k != m_indexSolvent || m_Aionic[k] != 0.0) {
+ if (k != 0 || m_Aionic[k] != 0.0) {
y = denomTmp * m_Aionic[k];
yp1 = y + 1.0;
sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
@@ -931,9 +887,9 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
}
lnActivitySolvent += coeff * tmp;
tmp = 0.0;
- for (size_t k = 0; k < m_kk; k++) {
+ for (size_t k = 1; k < m_kk; k++) {
z_k = m_speciesCharge[k];
- if ((k != m_indexSolvent) && (z_k != 0.0)) {
+ if (z_k != 0.0) {
tmp += m_B_Dot[k] * m_molalities[k];
}
}
@@ -967,9 +923,9 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
2.0 /3.0 * m_A_Debye * m_Mnaught *
m_IionicMolality * sqrt(m_IionicMolality) * sigma;
tmp = 0.0;
- for (size_t k = 0; k < m_kk; k++) {
+ for (size_t k = 1; k < m_kk; k++) {
z_k = m_speciesCharge[k];
- if ((k != m_indexSolvent) && (z_k != 0.0)) {
+ if (z_k != 0.0) {
tmp += m_B_Dot[k] * m_molalities[k];
}
}
@@ -983,15 +939,13 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
lnActivitySolvent =
(xmolSolvent - 1.0)/xmolSolvent;
- for (size_t k = 0; k < m_kk; k++) {
- if (k != m_indexSolvent) {
- z_k = m_speciesCharge[k];
- m_lnActCoeffMolal[k] =
- - z_k * z_k * numTmp / (1.0 + denomTmp);
- for (size_t j = 0; j < m_kk; j++) {
- double beta = m_Beta_ij.value(k, j);
- m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] * beta;
- }
+ for (size_t k = 1; k < m_kk; k++) {
+ z_k = m_speciesCharge[k];
+ m_lnActCoeffMolal[k] =
+ - z_k * z_k * numTmp / (1.0 + denomTmp);
+ for (size_t j = 0; j < m_kk; j++) {
+ double beta = m_Beta_ij.value(k, j);
+ m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] * beta;
}
}
if (denomTmp > 0.0) {
@@ -1020,18 +974,16 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
denomTmp *= m_Aionic[0];
numTmp = m_A_Debye * sqrt(m_IionicMolality);
tmpLn = log(1.0 + denomTmp);
- for (size_t k = 0; k < m_kk; k++) {
- if (k != m_indexSolvent) {
- z_k = m_speciesCharge[k];
- m_lnActCoeffMolal[k] =
- - z_k * z_k * numTmp / 3.0 / (1.0 + denomTmp);
- m_lnActCoeffMolal[k] +=
- - 2.0 * z_k * z_k * m_A_Debye * tmpLn /
- (3.0 * m_B_Debye * m_Aionic[0]);
- for (size_t j = 0; j < m_kk; j++) {
- m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] *
- m_Beta_ij.value(k, j);
- }
+ for (size_t k = 1; k < m_kk; k++) {
+ z_k = m_speciesCharge[k];
+ m_lnActCoeffMolal[k] =
+ - z_k * z_k * numTmp / 3.0 / (1.0 + denomTmp);
+ m_lnActCoeffMolal[k] +=
+ - 2.0 * z_k * z_k * m_A_Debye * tmpLn /
+ (3.0 * m_B_Debye * m_Aionic[0]);
+ for (size_t j = 0; j < m_kk; j++) {
+ m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] *
+ m_Beta_ij.value(k, j);
}
}
sigma = 1.0 / (1.0 + denomTmp);
@@ -1056,9 +1008,8 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
// Above, we calculated the ln(activitySolvent). Translate that into the
// molar-based activity coefficient by dividing by the solvent mole
// fraction. Solvents are not on the molality scale.
- xmolSolvent = moleFraction(m_indexSolvent);
- m_lnActCoeffMolal[m_indexSolvent] =
- lnActivitySolvent - log(xmolSolvent);
+ xmolSolvent = moleFraction(0);
+ m_lnActCoeffMolal[0] = lnActivitySolvent - log(xmolSolvent);
}
void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
@@ -1074,7 +1025,7 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
}
// Calculate a safe value for the mole fraction of the solvent
- double xmolSolvent = moleFraction(m_indexSolvent);
+ double xmolSolvent = moleFraction(0);
xmolSolvent = std::max(8.689E-3, xmolSolvent);
double sqrtI = sqrt(m_IionicMolality);
double numdAdTTmp = dAdT * sqrtI;
@@ -1089,7 +1040,7 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
}
d_lnActivitySolvent_dT = 2.0 / 3.0 * dAdT * m_Mnaught *
m_IionicMolality * sqrt(m_IionicMolality);
- m_dlnActCoeffMolaldT[m_indexSolvent] = d_lnActivitySolvent_dT;
+ m_dlnActCoeffMolaldT[0] = d_lnActivitySolvent_dT;
break;
case DHFORM_BDOT_AK:
@@ -1099,7 +1050,7 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
- z_k * z_k * numdAdTTmp / (1.0 + denomTmp * m_Aionic[k]);
}
- m_dlnActCoeffMolaldT[m_indexSolvent] = 0.0;
+ m_dlnActCoeffMolaldT[0] = 0.0;
coeff = 2.0 / 3.0 * dAdT * m_Mnaught * sqrtI;
tmp = 0.0;
if (denomTmp > 0.0) {
@@ -1111,7 +1062,7 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
}
}
- m_dlnActCoeffMolaldT[m_indexSolvent] += coeff * tmp;
+ m_dlnActCoeffMolaldT[0] += coeff * tmp;
break;
case DHFORM_BDOT_ACOMMON:
@@ -1128,19 +1079,15 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
} else {
sigma = 0.0;
}
- m_dlnActCoeffMolaldT[m_indexSolvent] =
- 2.0 /3.0 * dAdT * m_Mnaught *
+ m_dlnActCoeffMolaldT[0] = 2.0 /3.0 * dAdT * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
case DHFORM_BETAIJ:
denomTmp *= m_Aionic[0];
- for (size_t k = 0; k < m_kk; k++) {
- if (k != m_indexSolvent) {
- z_k = m_speciesCharge[k];
- m_dlnActCoeffMolaldT[k] =
- - z_k * z_k * numdAdTTmp / (1.0 + denomTmp);
- }
+ for (size_t k = 1; k < m_kk; k++) {
+ z_k = m_speciesCharge[k];
+ m_dlnActCoeffMolaldT[k] = -z_k*z_k * numdAdTTmp / (1.0 + denomTmp);
}
if (denomTmp > 0.0) {
y = denomTmp;
@@ -1149,28 +1096,23 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
} else {
sigma = 0.0;
}
- m_dlnActCoeffMolaldT[m_indexSolvent] =
- 2.0 /3.0 * dAdT * m_Mnaught *
+ m_dlnActCoeffMolaldT[0] = 2.0 /3.0 * dAdT * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
case DHFORM_PITZER_BETAIJ:
denomTmp *= m_Aionic[0];
tmpLn = log(1.0 + denomTmp);
- for (size_t k = 0; k < m_kk; k++) {
- if (k != m_indexSolvent) {
- z_k = m_speciesCharge[k];
- m_dlnActCoeffMolaldT[k] =
- - z_k * z_k * numdAdTTmp / (1.0 + denomTmp)
- - 2.0 * z_k * z_k * dAdT * tmpLn
- / (m_B_Debye * m_Aionic[0]);
- m_dlnActCoeffMolaldT[k] /= 3.0;
- }
+ for (size_t k = 1; k < m_kk; k++) {
+ z_k = m_speciesCharge[k];
+ m_dlnActCoeffMolaldT[k] =
+ - z_k * z_k * numdAdTTmp / (1.0 + denomTmp)
+ - 2.0 * z_k * z_k * dAdT * tmpLn / (m_B_Debye * m_Aionic[0]);
+ m_dlnActCoeffMolaldT[k] /= 3.0;
}
sigma = 1.0 / (1.0 + denomTmp);
- m_dlnActCoeffMolaldT[m_indexSolvent] =
- 2.0 /3.0 * dAdT * m_Mnaught *
+ m_dlnActCoeffMolaldT[0] = 2.0 /3.0 * dAdT * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
@@ -1193,7 +1135,7 @@ void DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2() const
}
// Calculate a safe value for the mole fraction of the solvent
- double xmolSolvent = moleFraction(m_indexSolvent);
+ double xmolSolvent = moleFraction(0);
xmolSolvent = std::max(8.689E-3, xmolSolvent);
double sqrtI = sqrt(m_IionicMolality);
double numd2AdT2Tmp = d2AdT2 * sqrtI;
@@ -1214,7 +1156,7 @@ void DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2() const
- z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp * m_Aionic[k]);
}
- m_d2lnActCoeffMolaldT2[m_indexSolvent] = 0.0;
+ m_d2lnActCoeffMolaldT2[0] = 0.0;
coeff = 2.0 / 3.0 * d2AdT2 * m_Mnaught * sqrtI;
tmp = 0.0;
if (denomTmp > 0.0) {
@@ -1226,7 +1168,7 @@ void DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2() const
tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
}
}
- m_d2lnActCoeffMolaldT2[m_indexSolvent] += coeff * tmp;
+ m_d2lnActCoeffMolaldT2[0] += coeff * tmp;
break;
case DHFORM_BDOT_ACOMMON:
@@ -1243,19 +1185,15 @@ void DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2() const
} else {
sigma = 0.0;
}
- m_d2lnActCoeffMolaldT2[m_indexSolvent] =
- 2.0 /3.0 * d2AdT2 * m_Mnaught *
+ m_d2lnActCoeffMolaldT2[0] = 2.0 /3.0 * d2AdT2 * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
case DHFORM_BETAIJ:
denomTmp *= m_Aionic[0];
- for (size_t k = 0; k < m_kk; k++) {
- if (k != m_indexSolvent) {
- z_k = m_speciesCharge[k];
- m_d2lnActCoeffMolaldT2[k] =
- - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp);
- }
+ for (size_t k = 1; k < m_kk; k++) {
+ z_k = m_speciesCharge[k];
+ m_d2lnActCoeffMolaldT2[k] = -z_k*z_k * numd2AdT2Tmp / (1.0 + denomTmp);
}
if (denomTmp > 0.0) {
y = denomTmp;
@@ -1264,28 +1202,23 @@ void DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2() const
} else {
sigma = 0.0;
}
- m_d2lnActCoeffMolaldT2[m_indexSolvent] =
- 2.0 /3.0 * d2AdT2 * m_Mnaught *
+ m_d2lnActCoeffMolaldT2[0] = 2.0 /3.0 * d2AdT2 * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
case DHFORM_PITZER_BETAIJ:
denomTmp *= m_Aionic[0];
tmpLn = log(1.0 + denomTmp);
- for (size_t k = 0; k < m_kk; k++) {
- if (k != m_indexSolvent) {
- z_k = m_speciesCharge[k];
- m_d2lnActCoeffMolaldT2[k] =
- - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp)
- - 2.0 * z_k * z_k * d2AdT2 * tmpLn
- / (m_B_Debye * m_Aionic[0]);
- m_d2lnActCoeffMolaldT2[k] /= 3.0;
- }
+ for (size_t k = 1; k < m_kk; k++) {
+ z_k = m_speciesCharge[k];
+ m_d2lnActCoeffMolaldT2[k] =
+ - z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp)
+ - 2.0 * z_k * z_k * d2AdT2 * tmpLn / (m_B_Debye * m_Aionic[0]);
+ m_d2lnActCoeffMolaldT2[k] /= 3.0;
}
sigma = 1.0 / (1.0 + denomTmp);
- m_d2lnActCoeffMolaldT2[m_indexSolvent] =
- 2.0 /3.0 * d2AdT2 * m_Mnaught *
+ m_d2lnActCoeffMolaldT2[0] = 2.0 /3.0 * d2AdT2 * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
@@ -1308,7 +1241,7 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dP() const
}
// Calculate a safe value for the mole fraction of the solvent
- double xmolSolvent = moleFraction(m_indexSolvent);
+ double xmolSolvent = moleFraction(0);
xmolSolvent = std::max(8.689E-3, xmolSolvent);
double sqrtI = sqrt(m_IionicMolality);
double numdAdPTmp = dAdP * sqrtI;
@@ -1334,7 +1267,7 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dP() const
}
}
- m_dlnActCoeffMolaldP[m_indexSolvent] = 0.0;
+ m_dlnActCoeffMolaldP[0] = 0.0;
coeff = 2.0 / 3.0 * dAdP * m_Mnaught * sqrtI;
tmp = 0.0;
if (denomTmp > 0.0) {
@@ -1346,7 +1279,7 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dP() const
tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
}
}
- m_dlnActCoeffMolaldP[m_indexSolvent] += coeff * tmp;
+ m_dlnActCoeffMolaldP[0] += coeff * tmp;
break;
case DHFORM_BDOT_ACOMMON:
@@ -1363,19 +1296,16 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dP() const
} else {
sigma = 0.0;
}
- m_dlnActCoeffMolaldP[m_indexSolvent] =
+ m_dlnActCoeffMolaldP[0] =
2.0 /3.0 * dAdP * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
case DHFORM_BETAIJ:
denomTmp *= m_Aionic[0];
- for (size_t k = 0; k < m_kk; k++) {
- if (k != m_indexSolvent) {
- z_k = m_speciesCharge[k];
- m_dlnActCoeffMolaldP[k] =
- - z_k * z_k * numdAdPTmp / (1.0 + denomTmp);
- }
+ for (size_t k = 1; k < m_kk; k++) {
+ z_k = m_speciesCharge[k];
+ m_dlnActCoeffMolaldP[k] = - z_k*z_k * numdAdPTmp / (1.0 + denomTmp);
}
if (denomTmp > 0.0) {
y = denomTmp;
@@ -1384,28 +1314,24 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dP() const
} else {
sigma = 0.0;
}
- m_dlnActCoeffMolaldP[m_indexSolvent] =
- 2.0 /3.0 * dAdP * m_Mnaught *
+ m_dlnActCoeffMolaldP[0] = 2.0 /3.0 * dAdP * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
case DHFORM_PITZER_BETAIJ:
denomTmp *= m_Aionic[0];
tmpLn = log(1.0 + denomTmp);
- for (size_t k = 0; k < m_kk; k++) {
- if (k != m_indexSolvent) {
- z_k = m_speciesCharge[k];
- m_dlnActCoeffMolaldP[k] =
- - z_k * z_k * numdAdPTmp / (1.0 + denomTmp)
- - 2.0 * z_k * z_k * dAdP * tmpLn
- / (m_B_Debye * m_Aionic[0]);
- m_dlnActCoeffMolaldP[k] /= 3.0;
- }
+ for (size_t k = 1; k < m_kk; k++) {
+ z_k = m_speciesCharge[k];
+ m_dlnActCoeffMolaldP[k] =
+ - z_k * z_k * numdAdPTmp / (1.0 + denomTmp)
+ - 2.0 * z_k * z_k * dAdP * tmpLn
+ / (m_B_Debye * m_Aionic[0]);
+ m_dlnActCoeffMolaldP[k] /= 3.0;
}
sigma = 1.0 / (1.0 + denomTmp);
- m_dlnActCoeffMolaldP[m_indexSolvent] =
- 2.0 /3.0 * dAdP * m_Mnaught *
+ m_dlnActCoeffMolaldP[0] = 2.0 /3.0 * dAdP * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
diff --git a/src/thermo/HMWSoln.cpp b/src/thermo/HMWSoln.cpp
index 99116c5cb..1ed47315a 100644
--- a/src/thermo/HMWSoln.cpp
+++ b/src/thermo/HMWSoln.cpp
@@ -307,7 +307,7 @@ void HMWSoln::getActivityConcentrations(doublereal* c) const
doublereal HMWSoln::standardConcentration(size_t k) const
{
getStandardVolumes(m_tmpV.data());
- double mvSolvent = m_tmpV[m_indexSolvent];
+ double mvSolvent = m_tmpV[0];
if (k > 0) {
return m_Mnaught / mvSolvent;
}
@@ -323,14 +323,11 @@ void HMWSoln::getActivities(doublereal* ac) const
s_update_lnMolalityActCoeff();
// Now calculate the array of activities.
- for (size_t k = 0; k < m_kk; k++) {
- if (k != m_indexSolvent) {
- ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal_Scaled[k]);
- }
+ for (size_t k = 1; k < m_kk; k++) {
+ ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal_Scaled[k]);
}
- double xmolSolvent = moleFraction(m_indexSolvent);
- ac[m_indexSolvent] =
- exp(m_lnActCoeffMolal_Scaled[m_indexSolvent]) * xmolSolvent;
+ double xmolSolvent = moleFraction(0);
+ ac[0] = exp(m_lnActCoeffMolal_Scaled[0]) * xmolSolvent;
}
void HMWSoln::getUnscaledMolalityActivityCoefficients(doublereal* acMolality) const
@@ -357,16 +354,13 @@ void HMWSoln::getChemPotentials(doublereal* mu) const
// Update the activity coefficients. This also updates the internal molality
// array.
s_update_lnMolalityActCoeff();
- 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]);
- }
+ double xmolSolvent = moleFraction(0);
+ for (size_t k = 1; k < m_kk; k++) {
+ xx = std::max(m_molalities[k], SmallNumber);
+ 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]);
+ mu[0] += RT() * (log(xx) + m_lnActCoeffMolal_Scaled[0]);
}
void HMWSoln::getPartialMolarEnthalpies(doublereal* hbar) const
@@ -406,15 +400,13 @@ void HMWSoln::getPartialMolarEntropies(doublereal* sbar) const
// First we will add in the obvious dependence on the T term out front of
// the log activity term
doublereal mm;
- for (size_t k = 0; k < m_kk; k++) {
- if (k != m_indexSolvent) {
- mm = std::max(SmallNumber, m_molalities[k]);
- sbar[k] -= GasConstant * (log(mm) + m_lnActCoeffMolal_Scaled[k]);
- }
+ for (size_t k = 1; k < m_kk; k++) {
+ mm = std::max(SmallNumber, m_molalities[k]);
+ sbar[k] -= GasConstant * (log(mm) + m_lnActCoeffMolal_Scaled[k]);
}
- double xmolSolvent = moleFraction(m_indexSolvent);
+ double xmolSolvent = moleFraction(0);
mm = std::max(SmallNumber, xmolSolvent);
- sbar[m_indexSolvent] -= GasConstant *(log(mm) + m_lnActCoeffMolal_Scaled[m_indexSolvent]);
+ sbar[0] -= GasConstant *(log(mm) + m_lnActCoeffMolal_Scaled[0]);
// Check to see whether activity coefficients are temperature dependent. If
// they are, then calculate the their temperature derivatives and add them
@@ -781,7 +773,7 @@ void HMWSoln::s_update_lnMolalityActCoeff() const
// Now do the main calculation.
s_updatePitzer_lnMolalityActCoeff();
- double xmolSolvent = moleFraction(m_indexSolvent);
+ double xmolSolvent = moleFraction(0);
double xx = std::max(m_xmolSolventMIN, xmolSolvent);
double lnActCoeffMolal0 = - log(xx) + (xx - 1.0)/xx;
double lnxs = log(xx);
@@ -929,7 +921,7 @@ void HMWSoln::calcMolalitiesCropped() const
if (cropMethod == 1) {
double* molF = m_gamma_tmp.data();
getMoleFractions(molF);
- double xmolSolvent = molF[m_indexSolvent];
+ double xmolSolvent = molF[0];
if (xmolSolvent >= MC_X_o_cutoff_) {
return;
}
@@ -1231,12 +1223,6 @@ void HMWSoln::s_updatePitzer_CoeffWRTemp(int doDerivs) const
void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
{
- // HKM -> Assumption is made that the solvent is species 0.
- if (m_indexSolvent != 0) {
- throw CanteraError("HMWSoln::s_updatePitzer_lnMolalityActCoeff",
- "Wrong index solvent value!");
- }
-
// Use the CROPPED molality of the species in solution.
const vector_fp& molality = m_molalitiesCropped;
@@ -1913,7 +1899,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
//
// We have just computed act_0. However, this routine returns
// ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
- double xmolSolvent = moleFraction(m_indexSolvent);
+ double xmolSolvent = moleFraction(0);
double xx = std::max(m_xmolSolventMIN, xmolSolvent);
m_lnActCoeffMolal_Unscaled[0] = lnwateract - log(xx);
if (m_debugCalc) {
@@ -1959,12 +1945,6 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
// immediately preceding the calling of this routine. Therefore, some
// quantities do not need to be recalculated in this routine.
- // HKM -> Assumption is made that the solvent is species 0.
- if (m_indexSolvent != 0) {
- throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
- "Wrong index solvent value!");
- }
-
const vector_fp& molality = m_molalitiesCropped;
double* d_gamma_dT_Unscaled = m_gamma_tmp.data();
@@ -2569,12 +2549,6 @@ void HMWSoln::s_update_d2lnMolalityActCoeff_dT2() const
void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
{
- // HKM -> Assumption is made that the solvent is species 0.
- if (m_indexSolvent != 0) {
- throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
- "Wrong index solvent value!");
- }
-
const double* molality = m_molalitiesCropped.data();
// Local variables defined by Coltrin
@@ -3171,12 +3145,6 @@ void HMWSoln::s_update_dlnMolalityActCoeff_dP() const
void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
{
- // HKM -> Assumption is made that the solvent is species 0.
- if (m_indexSolvent != 0) {
- throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
- "Wrong index solvent value!");
- }
-
const double* molality = m_molalitiesCropped.data();
// Local variables defined by Coltrin
@@ -3830,27 +3798,27 @@ void HMWSoln::s_updateIMS_lnMolalityActCoeff() const
// Calculate the molalities. Currently, the molalities may not be current
// with respect to the contents of the State objects' data.
calcMolalities();
- double xmolSolvent = moleFraction(m_indexSolvent);
+ double xmolSolvent = moleFraction(0);
double xx = std::max(m_xmolSolventMIN, xmolSolvent);
if (IMS_typeCutoff_ == 0) {
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
- IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
+ IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
return;
} else if (IMS_typeCutoff_ == 1) {
if (xmolSolvent > 3.0 * IMS_X_o_cutoff_/2.0) {
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
- IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
+ IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
return;
} else if (xmolSolvent < IMS_X_o_cutoff_/2.0) {
double tmp = log(xx * IMS_gamma_k_min_);
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
- IMS_lnActCoeffMolal_[m_indexSolvent] = log(IMS_gamma_o_min_);
+ IMS_lnActCoeffMolal_[0] = log(IMS_gamma_o_min_);
return;
} else {
// If we are in the middle region, calculate the connecting polynomials
@@ -3887,7 +3855,7 @@ void HMWSoln::s_updateIMS_lnMolalityActCoeff() const
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
- IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
+ IMS_lnActCoeffMolal_[0] = lngammao;
}
} else if (IMS_typeCutoff_ == 2) {
// Exponentials - trial 2
@@ -3895,7 +3863,7 @@ void HMWSoln::s_updateIMS_lnMolalityActCoeff() const
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
- IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
+ IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
return;
} else {
double xoverc = xmolSolvent/IMS_cCut_;
@@ -3921,7 +3889,7 @@ void HMWSoln::s_updateIMS_lnMolalityActCoeff() const
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
- IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
+ IMS_lnActCoeffMolal_[0] = lngammao;
}
}
return;
diff --git a/src/thermo/HMWSoln_input.cpp b/src/thermo/HMWSoln_input.cpp
index d1cc3d730..dd376a299 100644
--- a/src/thermo/HMWSoln_input.cpp
+++ b/src/thermo/HMWSoln_input.cpp
@@ -993,49 +993,10 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
- // Get the Name of the Solvent:
- // solventName
- string solventName = "";
- if (thermoNode.hasChild("solvent")) {
- XML_Node& scNode = thermoNode.child("solvent");
- vector nameSolventa;
- getStringArray(scNode, nameSolventa);
- if (nameSolventa.size() != 1) {
- throw CanteraError("HMWSoln::initThermoXML",
- "badly formed solvent XML node");
- }
- solventName = nameSolventa[0];
- }
-
// Initialize all of the lengths of arrays in the object
// now that we know what species are in the phase.
initLengths();
- // Reconcile the solvent name and index.
- for (size_t k = 0; k < m_kk; k++) {
- string sname = speciesName(k);
- if (solventName == sname) {
- setSolvent(k);
- if (k != 0) {
- throw CanteraError("HMWSoln::initThermoXML",
- "Solvent must be species 0 atm");
- }
- m_indexSolvent = k;
- break;
- }
- }
- if (m_indexSolvent == npos) {
- std::cout << "HMWSoln::initThermo: Solvent Name not found"
- << std::endl;
- throw CanteraError("HMWSoln::initThermoXML",
- "Solvent name not found");
- }
- if (m_indexSolvent != 0) {
- throw CanteraError("HMWSoln::initThermoXML",
- "Solvent " + solventName +
- " should be first species");
- }
-
// Now go get the specification of the standard states for species in the
// solution. This includes the molar volumes data blocks for incompressible
// species.
@@ -1239,7 +1200,7 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
m_electrolyteSpeciesType[k] = cEST_nonpolarNeutral;
}
}
- m_electrolyteSpeciesType[m_indexSolvent] = cEST_solvent;
+ m_electrolyteSpeciesType[0] = cEST_solvent;
// First look at the species database. Look for the subelement
// "stoichIsMods" in each of the species SS databases.
diff --git a/src/thermo/IdealMolalSoln.cpp b/src/thermo/IdealMolalSoln.cpp
index 48a9e4cbf..7ef3b7fb4 100644
--- a/src/thermo/IdealMolalSoln.cpp
+++ b/src/thermo/IdealMolalSoln.cpp
@@ -180,10 +180,10 @@ doublereal IdealMolalSoln::standardConcentration(size_t k) const
case 0:
break;
case 1:
- return c0 = 1.0 /m_speciesMolarVolume[m_indexSolvent];
+ return c0 = 1.0 /m_speciesMolarVolume[0];
break;
case 2:
- c0 = 1.0 / m_speciesMolarVolume[m_indexSolvent];
+ c0 = 1.0 / m_speciesMolarVolume[0];
break;
}
return c0;
@@ -200,12 +200,11 @@ void IdealMolalSoln::getActivities(doublereal* ac) const
for (size_t k = 0; k < m_kk; k++) {
ac[k] = m_molalities[k];
}
- double xmolSolvent = moleFraction(m_indexSolvent);
+ double xmolSolvent = moleFraction(0);
// Limit the activity coefficient to be finite as the solvent mole
// fraction goes to zero.
xmolSolvent = std::max(m_xmolSolventMIN, xmolSolvent);
- ac[m_indexSolvent] =
- exp((xmolSolvent - 1.0)/xmolSolvent);
+ ac[0] = exp((xmolSolvent - 1.0)/xmolSolvent);
} else {
s_updateIMS_lnMolalityActCoeff();
@@ -214,9 +213,8 @@ void IdealMolalSoln::getActivities(doublereal* ac) const
for (size_t k = 1; k < m_kk; k++) {
ac[k] = m_molalities[k] * exp(IMS_lnActCoeffMolal_[k]);
}
- double xmolSolvent = moleFraction(m_indexSolvent);
- ac[m_indexSolvent] =
- exp(IMS_lnActCoeffMolal_[m_indexSolvent]) * xmolSolvent;
+ double xmolSolvent = moleFraction(0);
+ ac[0] = exp(IMS_lnActCoeffMolal_[0]) * xmolSolvent;
}
}
@@ -226,12 +224,11 @@ void IdealMolalSoln::getMolalityActivityCoefficients(doublereal* acMolality) con
for (size_t k = 0; k < m_kk; k++) {
acMolality[k] = 1.0;
}
- double xmolSolvent = moleFraction(m_indexSolvent);
+ double xmolSolvent = moleFraction(0);
// Limit the activity coefficient to be finite as the solvent mole
// fraction goes to zero.
xmolSolvent = std::max(m_xmolSolventMIN, xmolSolvent);
- acMolality[m_indexSolvent] =
- exp((xmolSolvent - 1.0)/xmolSolvent) / xmolSolvent;
+ acMolality[0] = exp((xmolSolvent - 1.0)/xmolSolvent) / xmolSolvent;
} else {
s_updateIMS_lnMolalityActCoeff();
std::copy(IMS_lnActCoeffMolal_.begin(), IMS_lnActCoeffMolal_.end(), acMolality);
@@ -245,9 +242,6 @@ void IdealMolalSoln::getMolalityActivityCoefficients(doublereal* acMolality) con
void IdealMolalSoln::getChemPotentials(doublereal* mu) const
{
- // Assertion is made for speed
- AssertThrow(m_indexSolvent == 0, "solvent not the first species");
-
// First get the standard chemical potentials. This requires updates of
// standard state as a function of T and P These are defined at unit
// molality.
@@ -258,7 +252,7 @@ void IdealMolalSoln::getChemPotentials(doublereal* mu) const
calcMolalities();
// get the solvent mole fraction
- double xmolSolvent = moleFraction(m_indexSolvent);
+ double xmolSolvent = moleFraction(0);
if (IMS_typeCutoff_ == 0 || xmolSolvent > 3.* IMS_X_o_cutoff_/2.0) {
for (size_t k = 1; k < m_kk; k++) {
@@ -269,8 +263,7 @@ void IdealMolalSoln::getChemPotentials(doublereal* mu) const
// Do the solvent
// -> see my notes
double xx = std::max(xmolSolvent, SmallNumber);
- mu[m_indexSolvent] +=
- (RT() * (xmolSolvent - 1.0) / xx);
+ mu[0] += (RT() * (xmolSolvent - 1.0) / xx);
} else {
// Update the activity coefficients. This also updates the internal
// molality array.
@@ -281,8 +274,7 @@ void IdealMolalSoln::getChemPotentials(doublereal* mu) const
mu[k] += RT() * (log(xx) + IMS_lnActCoeffMolal_[k]);
}
double xx = std::max(xmolSolvent, SmallNumber);
- mu[m_indexSolvent] +=
- RT() * (log(xx) + IMS_lnActCoeffMolal_[m_indexSolvent]);
+ mu[0] += RT() * (log(xx) + IMS_lnActCoeffMolal_[0]);
}
}
@@ -299,14 +291,12 @@ void IdealMolalSoln::getPartialMolarEntropies(doublereal* sbar) const
getEntropy_R(sbar);
calcMolalities();
if (IMS_typeCutoff_ == 0) {
- for (size_t k = 0; k < m_kk; k++) {
- if (k != m_indexSolvent) {
- doublereal mm = std::max(SmallNumber, m_molalities[k]);
- sbar[k] -= GasConstant * log(mm);
- }
+ for (size_t k = 1; k < m_kk; k++) {
+ doublereal mm = std::max(SmallNumber, m_molalities[k]);
+ sbar[k] -= GasConstant * log(mm);
}
- double xmolSolvent = moleFraction(m_indexSolvent);
- sbar[m_indexSolvent] -= (GasConstant * (xmolSolvent - 1.0) / xmolSolvent);
+ double xmolSolvent = moleFraction(0);
+ sbar[0] -= (GasConstant * (xmolSolvent - 1.0) / xmolSolvent);
} else {
// Update the activity coefficients, This also update the internally
// stored molalities.
@@ -315,15 +305,13 @@ void IdealMolalSoln::getPartialMolarEntropies(doublereal* sbar) const
// First we will add in the obvious dependence on the T term out front
// of the log activity term
doublereal mm;
- for (size_t k = 0; k < m_kk; k++) {
- if (k != m_indexSolvent) {
- mm = std::max(SmallNumber, m_molalities[k]);
- sbar[k] -= GasConstant * (log(mm) + IMS_lnActCoeffMolal_[k]);
- }
+ for (size_t k = 1; k < m_kk; k++) {
+ mm = std::max(SmallNumber, m_molalities[k]);
+ sbar[k] -= GasConstant * (log(mm) + IMS_lnActCoeffMolal_[k]);
}
- double xmolSolvent = moleFraction(m_indexSolvent);
+ double xmolSolvent = moleFraction(0);
mm = std::max(SmallNumber, xmolSolvent);
- sbar[m_indexSolvent] -= GasConstant *(log(mm) + IMS_lnActCoeffMolal_[m_indexSolvent]);
+ sbar[0] -= GasConstant *(log(mm) + IMS_lnActCoeffMolal_[0]);
}
}
@@ -377,19 +365,6 @@ void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
setStandardConcentrationModel(scNode["model"]);
}
- // Get the Name of the Solvent:
- // solventName
- std::string solventName = "";
- if (thermoNode.hasChild("solvent")) {
- std::vector nameSolventa;
- getStringArray(thermoNode.child("solvent"), nameSolventa);
- if (nameSolventa.size() != 1) {
- throw CanteraError("IdealMolalSoln::initThermoXML",
- "badly formed solvent XML node");
- }
- solventName = nameSolventa[0];
- }
-
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& acNode = thermoNode.child("activityCoefficients");
std::string modelString = acNode.attrib("model");
@@ -425,25 +400,6 @@ void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
setCutoffModel("none");
}
}
-
- // Reconcile the solvent name and index.
- for (size_t k = 0; k < m_kk; k++) {
- if (solventName == speciesName(k)) {
- m_indexSolvent = k;
- break;
- }
- }
- if (m_indexSolvent == npos) {
- std::cout << "IdealMolalSoln::initThermo: Solvent Name not found"
- << std::endl;
- throw CanteraError("IdealMolalSoln::initThermo",
- "Solvent name not found");
- }
- if (m_indexSolvent != 0) {
- throw CanteraError("IdealMolalSoln::initThermo",
- "Solvent " + solventName +
- " should be first species");
- }
}
void IdealMolalSoln::initThermo()
@@ -494,28 +450,28 @@ void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
// with respect to the contents of the State objects' data.
calcMolalities();
- double xmolSolvent = moleFraction(m_indexSolvent);
+ double xmolSolvent = moleFraction(0);
double xx = std::max(m_xmolSolventMIN, xmolSolvent);
if (IMS_typeCutoff_ == 0) {
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
- IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
+ IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
return;
} else if (IMS_typeCutoff_ == 1) {
if (xmolSolvent > 3.0 * IMS_X_o_cutoff_/2.0) {
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
- IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
+ IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
return;
} else if (xmolSolvent < IMS_X_o_cutoff_/2.0) {
double tmp = log(xx * IMS_gamma_k_min_);
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
- IMS_lnActCoeffMolal_[m_indexSolvent] = log(IMS_gamma_o_min_);
+ IMS_lnActCoeffMolal_[0] = log(IMS_gamma_o_min_);
return;
} else {
// If we are in the middle region, calculate the connecting polynomials
@@ -552,7 +508,7 @@ void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
- IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
+ IMS_lnActCoeffMolal_[0] = lngammao;
}
} else if (IMS_typeCutoff_ == 2) {
// Exponentials - trial 2
@@ -560,7 +516,7 @@ void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
- IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
+ IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
return;
} else {
double xoverc = xmolSolvent/IMS_cCut_;
@@ -584,7 +540,7 @@ void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
- IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
+ IMS_lnActCoeffMolal_[0] = lngammao;
}
}
}
diff --git a/src/thermo/MolalityVPSSTP.cpp b/src/thermo/MolalityVPSSTP.cpp
index 18c8865f7..6e2e2febf 100644
--- a/src/thermo/MolalityVPSSTP.cpp
+++ b/src/thermo/MolalityVPSSTP.cpp
@@ -22,7 +22,6 @@ namespace Cantera
{
MolalityVPSSTP::MolalityVPSSTP() :
- m_indexSolvent(0),
m_pHScalingType(PHSCALE_PITZER),
m_indexCLM(npos),
m_weightSolvent(18.01528),
@@ -53,20 +52,15 @@ int MolalityVPSSTP::pHScale() const
void MolalityVPSSTP::setSolvent(size_t k)
{
- if (k >= m_kk) {
- throw CanteraError("MolalityVPSSTP::setSolute ",
- "bad value");
- }
- m_indexSolvent = k;
- AssertThrowMsg(m_indexSolvent==0, "MolalityVPSSTP::setSolvent",
- "Molality-based methods limit solvent id to being 0");
- m_weightSolvent = molecularWeight(k);
- m_Mnaught = m_weightSolvent / 1000.;
+ warn_deprecated("MolalityVPSSTP::setSolvent", "Solvent is always the first"
+ " species. To be removed after Cantera 2.4.");
}
size_t MolalityVPSSTP::solventIndex() const
{
- return m_indexSolvent;
+ warn_deprecated("MolalityVPSSTP::solventIndex", "Solvent is always the"
+ " first species. To be removed after Cantera 2.4.");
+ return 0;
}
void MolalityVPSSTP::setMoleFSolventMin(doublereal xmolSolventMIN)
@@ -87,7 +81,7 @@ doublereal MolalityVPSSTP::moleFSolventMin() const
void MolalityVPSSTP::calcMolalities() const
{
getMoleFractions(m_molalities.data());
- double xmolSolvent = std::max(m_molalities[m_indexSolvent], m_xmolSolventMIN);
+ double xmolSolvent = std::max(m_molalities[0], m_xmolSolventMIN);
double denomInv = 1.0/ (m_Mnaught * xmolSolvent);
for (size_t k = 0; k < m_kk; k++) {
m_molalities[k] *= denomInv;
@@ -110,8 +104,8 @@ void MolalityVPSSTP::setMolalities(const doublereal* const molal)
Lsum += molal[k];
}
double tmp = 1.0 / Lsum;
- m_molalities[m_indexSolvent] = tmp / m_Mnaught;
- double sum = m_molalities[m_indexSolvent];
+ m_molalities[0] = tmp / m_Mnaught;
+ double sum = m_molalities[0];
for (size_t k = 1; k < m_kk; k++) {
m_molalities[k] = tmp * molal[k];
sum += m_molalities[k];
@@ -137,7 +131,7 @@ void MolalityVPSSTP::setMolalitiesByName(const compositionMap& mMap)
// Get a vector of mole fractions
vector_fp mf(m_kk, 0.0);
getMoleFractions(mf.data());
- double xmolSmin = std::max(mf[m_indexSolvent], m_xmolSolventMIN);
+ double xmolSmin = std::max(mf[0], m_xmolSolventMIN);
for (size_t k = 0; k < m_kk; k++) {
double mol_k = getValue(mMap, speciesName(k), 0.0);
if (mol_k > 0) {
@@ -228,8 +222,7 @@ void MolalityVPSSTP::getActivities(doublereal* ac) const
void MolalityVPSSTP::getActivityCoefficients(doublereal* ac) const
{
getMolalityActivityCoefficients(ac);
- AssertThrow(m_indexSolvent==0, "MolalityVPSSTP::getActivityCoefficients");
- double xmolSolvent = std::max(moleFraction(m_indexSolvent), m_xmolSolventMIN);
+ double xmolSolvent = std::max(moleFraction(0), m_xmolSolventMIN);
for (size_t k = 1; k < m_kk; k++) {
ac[k] /= xmolSolvent;
}
@@ -254,7 +247,7 @@ doublereal MolalityVPSSTP::osmoticCoefficient() const
}
double oc = 1.0;
if (sum > 1.0E-200) {
- oc = - log(act[m_indexSolvent]) / (m_Mnaught * sum);
+ oc = - log(act[0]) / (m_Mnaught * sum);
}
return oc;
}
@@ -371,7 +364,8 @@ bool MolalityVPSSTP::addSpecies(shared_ptr spec)
if (added) {
if (m_kk == 1) {
// The solvent defaults to species 0
- setSolvent(0);
+ m_weightSolvent = molecularWeight(0);
+ m_Mnaught = m_weightSolvent / 1000.;
}
m_molalities.push_back(0.0);
}