diff --git a/Cantera/src/thermo/HMWSoln.cpp b/Cantera/src/thermo/HMWSoln.cpp index f0fdc7fb0..38dc43f5a 100644 --- a/Cantera/src/thermo/HMWSoln.cpp +++ b/Cantera/src/thermo/HMWSoln.cpp @@ -277,6 +277,7 @@ namespace Cantera { m_Psi_ijk_L = b.m_Psi_ijk_L; m_Psi_ijk_LL = b.m_Psi_ijk_LL; m_Psi_ijk_P = b.m_Psi_ijk_P; + m_Psi_ijk_coeff = b.m_Psi_ijk_coeff; m_Lambda_ij = b.m_Lambda_ij; m_Lambda_ij_L = b.m_Lambda_ij_L; m_Lambda_ij_LL = b.m_Lambda_ij_LL; @@ -1589,10 +1590,12 @@ namespace Cantera { m_Theta_ij_P.resize(maxCounterIJlen, 0.0); m_Theta_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0); + int n = m_kk*m_kk*m_kk; m_Psi_ijk.resize(m_kk*m_kk*m_kk, 0.0); m_Psi_ijk_L.resize(m_kk*m_kk*m_kk, 0.0); m_Psi_ijk_LL.resize(m_kk*m_kk*m_kk, 0.0); m_Psi_ijk_P.resize(m_kk*m_kk*m_kk, 0.0); + m_Psi_ijk_coeff.resize(TCoeffLength, n, 0.0); m_Lambda_ij.resize(leng, leng, 0.0); m_Lambda_ij_L.resize(leng, leng, 0.0); @@ -2033,6 +2036,41 @@ namespace Cantera { } } + + for (i = 0; i < m_kk; i++) { + for (j = 0; j < m_kk; j++) { + for (int k = 0; k < m_kk; k++) { + n = i * m_kk *m_kk + j * m_kk + k ; + const double *Psi_coeff = m_Psi_ijk_coeff.ptrColumn(n); + switch (m_formPitzerTemp) { + case PITZER_TEMP_CONSTANT: + m_Psi_ijk[n] = Psi_coeff[n]; + break; + case PITZER_TEMP_LINEAR: + m_Psi_ijk[n] = Psi_coeff[0] + Psi_coeff[1]*tlin; + m_Psi_ijk_L[n] = Psi_coeff[1]; + m_Psi_ijk_LL[n] = 0.0; + case PITZER_TEMP_COMPLEX1: + m_Psi_ijk[n] = Psi_coeff[0] + + Psi_coeff[1]*tlin + + Psi_coeff[2]*tquad + + Psi_coeff[3]*tinv + + Psi_coeff[4]*tln; + + m_Psi_ijk_L[n] = Psi_coeff[1] + + Psi_coeff[2]*2.0*T + - Psi_coeff[3]/(T*T) + + Psi_coeff[4]/T; + + m_Psi_ijk_LL[n] = + Psi_coeff[2]*2.0 + + 2.0*Psi_coeff[3]/(T*T*T) + - Psi_coeff[4]/(T*T); + } + } + } + } + } /* diff --git a/Cantera/src/thermo/HMWSoln.h b/Cantera/src/thermo/HMWSoln.h index 4f6c97246..6975458d9 100644 --- a/Cantera/src/thermo/HMWSoln.h +++ b/Cantera/src/thermo/HMWSoln.h @@ -2644,25 +2644,42 @@ namespace Cantera { * The first two coordinates are symmetric wrt cations, * and the last two coordinates are symmetric wrt anions. */ - vector_fp m_Psi_ijk; + mutable vector_fp m_Psi_ijk; //! Derivitive of Psi_ijk[n] wrt T /*! * see m_Psi_ijk for reference on the indexing into this variable. */ - vector_fp m_Psi_ijk_L; + mutable vector_fp m_Psi_ijk_L; //! Derivitive of Psi_ijk[n] wrt TT /*! * see m_Psi_ijk for reference on the indexing into this variable. */ - vector_fp m_Psi_ijk_LL; + mutable vector_fp m_Psi_ijk_LL; //! Derivitive of Psi_ijk[n] wrt P /*! * see m_Psi_ijk for reference on the indexing into this variable. */ - vector_fp m_Psi_ijk_P; + mutable vector_fp m_Psi_ijk_P; + + //! Array of coefficients for Psi_ijk[n] in the Pitzer/HMW formulation. + /*! + * Psi_ijk[n] is the value of the psi coefficient for the + * ijk interaction where + * + * n = k + j * m_kk + i * m_kk * m_kk; + * + * It is potentially nonzero everywhere. + * The first two coordinates are symmetric wrt cations, + * and the last two coordinates are symmetric wrt anions. + * + * + * m_Psi_ijk_coeff.ptrColumn(n) is a double* containing + * the vector of coefficients for the n interaction. + */ + Array2D m_Psi_ijk_coeff; //! Lambda coefficient for the ij interaction /*! diff --git a/Cantera/src/thermo/HMWSoln_input.cpp b/Cantera/src/thermo/HMWSoln_input.cpp index 00b720ece..d5bba130b 100644 --- a/Cantera/src/thermo/HMWSoln_input.cpp +++ b/Cantera/src/thermo/HMWSoln_input.cpp @@ -452,6 +452,8 @@ namespace Cantera { } double *charge = DATA_PTR(m_speciesCharge); string stemp; + vector_fp vParams; + int nParamsFound = 0; string kName = BinSalt.attrib("cation"); if (kName == "") { throw CanteraError("HMWSoln::readXMLPsiCommonCation", "no cation attrib"); @@ -512,25 +514,76 @@ namespace Cantera { } } if (nodeName == "psi") { - stemp = xmlChild.value(); - double param = atofCheck(stemp.c_str()); + getFloatArray(xmlChild, vParams, false, "", "Psi"); + nParamsFound = vParams.size(); n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies ; - m_Psi_ijk[n] = param; + + if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) { + if (nParamsFound != 1) { + throw CanteraError("HMWSoln::readXMLPsiCommonCation::Psi for " + + kName + "::" + iName + "::" + jName, + "wrong number of params found"); + } + m_Psi_ijk_coeff(0,n) = vParams[0]; + m_Psi_ijk[n] = vParams[0]; + } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) { + if (nParamsFound != 2) { + throw CanteraError("HMWSoln::readXMLPsiCation::Psi for " + + kName + "::" + iName + "::" + jName, + "wrong number of params found"); + } + m_Psi_ijk_coeff(0,n) = vParams[0]; + m_Psi_ijk_coeff(1,n) = vParams[1]; + m_Psi_ijk[n] = vParams[0]; + } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) { + if (nParamsFound == 1) { + vParams.resize(5, 0.0); + nParamsFound = 5; + } else if (nParamsFound != 5) { + throw CanteraError("HMWSoln::readXMLPsiCation::Psi for " + + kName + "::" + iName + "::" + jName, + "wrong number of params found"); + } + for (i = 0; i < nParamsFound; i++) { + m_Psi_ijk_coeff(i, n) = vParams[i]; + } + m_Psi_ijk[n] = vParams[0]; + } + + + // fill in the duplicate entries n = iSpecies * m_kk *m_kk + kSpecies * m_kk + jSpecies ; - m_Psi_ijk[n] = param; + for (i = 0; i < nParamsFound; i++) { + m_Psi_ijk_coeff(i, n) = vParams[i]; + } + m_Psi_ijk[n] = vParams[0]; + n = jSpecies * m_kk *m_kk + iSpecies * m_kk + kSpecies ; - m_Psi_ijk[n] = param; + for (i = 0; i < nParamsFound; i++) { + m_Psi_ijk_coeff(i, n) = vParams[i]; + } + m_Psi_ijk[n] = vParams[0]; + n = jSpecies * m_kk *m_kk + kSpecies * m_kk + iSpecies ; - m_Psi_ijk[n] = param; + for (i = 0; i < nParamsFound; i++) { + m_Psi_ijk_coeff(i, n) = vParams[i]; + } + m_Psi_ijk[n] = vParams[0]; + n = kSpecies * m_kk *m_kk + jSpecies * m_kk + iSpecies ; - m_Psi_ijk[n] = param; + for (i = 0; i < nParamsFound; i++) { + m_Psi_ijk_coeff(i, n) = vParams[i]; + } + m_Psi_ijk[n] = vParams[0]; + n = kSpecies * m_kk *m_kk + iSpecies * m_kk + jSpecies ; - m_Psi_ijk[n] = param; + for (i = 0; i < nParamsFound; i++) { + m_Psi_ijk_coeff(i, n) = vParams[i]; + } + m_Psi_ijk[n] = vParams[0]; } } } - - /** * Process an XML node called "PsiCommonAnion". @@ -545,6 +598,8 @@ namespace Cantera { } double *charge = DATA_PTR(m_speciesCharge); string stemp; + vector_fp vParams; + int nParamsFound = 0; string kName = BinSalt.attrib("anion"); if (kName == "") { throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "no anion attrib"); @@ -604,20 +659,75 @@ namespace Cantera { } } if (nodeName == "psi") { - stemp = xmlChild.value(); - double param = atofCheck(stemp.c_str()); + + getFloatArray(xmlChild, vParams, false, "", "Psi"); + nParamsFound = vParams.size(); n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies ; - m_Psi_ijk[n] = param; + + if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) { + if (nParamsFound != 1) { + throw CanteraError("HMWSoln::readXMLPsiCommonAnion::Psi for " + + kName + "::" + iName + "::" + jName, + "wrong number of params found"); + } + m_Psi_ijk_coeff(0,n) = vParams[0]; + m_Psi_ijk[n] = vParams[0]; + } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) { + if (nParamsFound != 2) { + throw CanteraError("HMWSoln::readXMLPsiAnion::Psi for " + + kName + "::" + iName + "::" + jName, + "wrong number of params found"); + } + m_Psi_ijk_coeff(0,n) = vParams[0]; + m_Psi_ijk_coeff(1,n) = vParams[1]; + m_Psi_ijk[n] = vParams[0]; + } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) { + if (nParamsFound == 1) { + vParams.resize(5, 0.0); + nParamsFound = 5; + } else if (nParamsFound != 5) { + throw CanteraError("HMWSoln::readXMLPsiAnion::Psi for " + + kName + "::" + iName + "::" + jName, + "wrong number of params found"); + } + for (i = 0; i < nParamsFound; i++) { + m_Psi_ijk_coeff(i, n) = vParams[i]; + } + m_Psi_ijk[n] = vParams[0]; + } + + + // fill in the duplicate entries n = iSpecies * m_kk *m_kk + kSpecies * m_kk + jSpecies ; - m_Psi_ijk[n] = param; + for (i = 0; i < nParamsFound; i++) { + m_Psi_ijk_coeff(i, n) = vParams[i]; + } + m_Psi_ijk[n] = vParams[0]; + n = jSpecies * m_kk *m_kk + iSpecies * m_kk + kSpecies ; - m_Psi_ijk[n] = param; + for (i = 0; i < nParamsFound; i++) { + m_Psi_ijk_coeff(i, n) = vParams[i]; + } + m_Psi_ijk[n] = vParams[0]; + n = jSpecies * m_kk *m_kk + kSpecies * m_kk + iSpecies ; - m_Psi_ijk[n] = param; + for (i = 0; i < nParamsFound; i++) { + m_Psi_ijk_coeff(i, n) = vParams[i]; + } + m_Psi_ijk[n] = vParams[0]; + n = kSpecies * m_kk *m_kk + jSpecies * m_kk + iSpecies ; - m_Psi_ijk[n] = param; + for (i = 0; i < nParamsFound; i++) { + m_Psi_ijk_coeff(i, n) = vParams[i]; + } + m_Psi_ijk[n] = vParams[0]; + n = kSpecies * m_kk *m_kk + iSpecies * m_kk + jSpecies ; - m_Psi_ijk[n] = param; + for (i = 0; i < nParamsFound; i++) { + m_Psi_ijk_coeff(i, n) = vParams[i]; + } + m_Psi_ijk[n] = vParams[0]; + } } }