From e33fe6904d573ddd141ce5a619a24fb873137746 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 2 Aug 2017 22:39:40 -0400 Subject: [PATCH] [Thermo] Solvent is always the first species No useful capabilities are provided by allowing the solvent species to vary, and there are many places where the solvent was already implicitly assumed to be the first species. --- include/cantera/thermo/MolalityVPSSTP.h | 9 +- src/thermo/DebyeHuckel.cpp | 252 +++++++++--------------- src/thermo/HMWSoln.cpp | 82 +++----- src/thermo/HMWSoln_input.cpp | 41 +--- src/thermo/IdealMolalSoln.cpp | 100 +++------- src/thermo/MolalityVPSSTP.cpp | 32 ++- 6 files changed, 160 insertions(+), 356 deletions(-) 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); }