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

[Thermo] Improve usage of local variables

Browse source 
Also eliminate some unnecessary temporary variables.
This commit is contained in:
Ray Speth 2015-02-12 04:04:18 +00:00
parent 004798aef4
commit 7acf0378c6
51 changed files with 1260 additions and 1990 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
src/thermo/ConstDensityThermo.cpp
View file

@ -106,11 +106,10 @@ void ConstDensityThermo::getChemPotentials(doublereal* mu) const
{
doublereal vdp = (pressure() - m_spthermo->refPressure())/
molarDensity();
doublereal xx;
doublereal rt = temperature() * GasConstant;
const vector_fp& g_RT = gibbs_RT();
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(SmallNumber, moleFraction(k));
double xx = std::max(SmallNumber, moleFraction(k));
mu[k] = rt*(g_RT[k] + log(xx)) + vdp;
}
}

30
src/thermo/DebyeHuckel.cpp
View file

@ -628,8 +628,6 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
m_formDH = DHFORM_DILUTE_LIMIT;
}
std::string stemp;
/*
* Possibly change the form of the standard concentrations
*/
@ -923,8 +921,9 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
* lack of agreement (HKM -> may be changed in the
* future).
*/
map<std::string,std::string>::const_iterator _b = m.begin();
for (; _b != m.end(); ++_b) {
for (map<std::string,std::string>::const_iterator _b = m.begin();
_b != m.end();
++_b) {
size_t kk = speciesIndex(_b->first);
m_Aionic[kk] = fpValue(_b->second) * Afactor;
@ -966,14 +965,13 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
* in each of the species SS databases.
*/
std::vector<const XML_Node*> xspecies= speciesData();
std::string kname, jname;
size_t jj = xspecies.size();
for (size_t k = 0; k < m_kk; k++) {
size_t jmap = npos;
kname = speciesName(k);
std::string kname = speciesName(k);
for (size_t j = 0; j < jj; j++) {
const XML_Node& sp = *xspecies[j];
jname = sp["name"];
std::string jname = sp["name"];
if (jname == kname) {
jmap = j;
break;
@ -997,8 +995,9 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
map<std::string, std::string> msIs;
getMap(sIsNode, msIs);
map<std::string,std::string>::const_iterator _b = msIs.begin();
for (; _b != msIs.end(); ++_b) {
for (map<std::string,std::string>::const_iterator _b = msIs.begin();
_b != msIs.end();
++_b) {
size_t kk = speciesIndex(_b->first);
double val = fpValue(_b->second);
m_speciesCharge_Stoich[kk] = val;
@ -1034,11 +1033,9 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
* in each of the species SS databases.
*/
std::vector<const XML_Node*> xspecies= speciesData();
const XML_Node* spPtr = 0;
std::string kname;
for (size_t k = 0; k < m_kk; k++) {
kname = speciesName(k);
spPtr = xspecies[k];
std::string kname = speciesName(k);
const XML_Node* spPtr = xspecies[k];
if (!spPtr) {
if (spPtr->hasChild("electrolyteSpeciesType")) {
std::string est = getChildValue(*spPtr, "electrolyteSpeciesType");
@ -1057,8 +1054,9 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType");
map<std::string, std::string> msEST;
getMap(ESTNode, msEST);
map<std::string,std::string>::const_iterator _b = msEST.begin();
for (; _b != msEST.end(); ++_b) {
for (map<std::string,std::string>::const_iterator _b = msEST.begin();
_b != msEST.end();
++_b) {
size_t kk = speciesIndex(_b->first);
std::string est = _b->second;
if ((m_electrolyteSpeciesType[kk] = interp_est(est)) == -1) {
@ -1069,8 +1067,6 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Lastly set the state
*/

895
src/thermo/HMWSoln.cpp
View file

File diff suppressed because it is too large Load diff

99
src/thermo/HMWSoln_input.cpp
View file

@ -58,8 +58,7 @@ void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
throw CanteraError("HMWSoln::readXMLBinarySalt",
"Incorrect name for processing this routine: " + xname);
}
string stemp;
size_t nParamsFound, i;
vector_fp vParams;
string iName = BinSalt.attrib("cation");
if (iName == "") {
@ -94,7 +93,7 @@ void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
int counter = m_CounterIJ[n];
for (size_t iChild = 0; iChild < BinSalt.nChildren(); iChild++) {
XML_Node& xmlChild = BinSalt.child(iChild);
stemp = xmlChild.name();
string stemp = xmlChild.name();
string nodeName = lowercase(stemp);
/*
* Process the binary salt child elements
@ -104,7 +103,7 @@ void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
* Get the string containing all of the values
*/
getFloatArray(xmlChild, vParams, false, "", "beta0");
nParamsFound = vParams.size();
size_t nParamsFound = vParams.size();
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
if (nParamsFound != 1) {
throw CanteraError("HMWSoln::readXMLBinarySalt::beta0 for " + ispName
@ -128,7 +127,7 @@ void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
+ "::" + jspName,
"wrong number of params found");
}
for (i = 0; i < nParamsFound; i++) {
for (size_t i = 0; i < nParamsFound; i++) {
m_Beta0MX_ij_coeff(i, counter) = vParams[i];
}
m_Beta0MX_ij[counter] = vParams[0];
@ -140,7 +139,7 @@ void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
* Get the string containing all of the values
*/
getFloatArray(xmlChild, vParams, false, "", "beta1");
nParamsFound = vParams.size();
size_t nParamsFound = vParams.size();
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
if (nParamsFound != 1) {
throw CanteraError("HMWSoln::readXMLBinarySalt::beta1 for " + ispName
@ -164,7 +163,7 @@ void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
+ "::" + jspName,
"wrong number of params found");
}
for (i = 0; i < nParamsFound; i++) {
for (size_t i = 0; i < nParamsFound; i++) {
m_Beta1MX_ij_coeff(i, counter) = vParams[i];
}
m_Beta1MX_ij[counter] = vParams[0];
@ -172,7 +171,7 @@ void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
}
if (nodeName == "beta2") {
getFloatArray(xmlChild, vParams, false, "", "beta2");
nParamsFound = vParams.size();
size_t nParamsFound = vParams.size();
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
if (nParamsFound != 1) {
throw CanteraError("HMWSoln::readXMLBinarySalt::beta2 for " + ispName
@ -196,7 +195,7 @@ void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
+ "::" + jspName,
"wrong number of params found");
}
for (i = 0; i < nParamsFound; i++) {
for (size_t i = 0; i < nParamsFound; i++) {
m_Beta2MX_ij_coeff(i, counter) = vParams[i];
}
m_Beta2MX_ij[counter] = vParams[0];
@ -208,7 +207,7 @@ void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
* Get the string containing all of the values
*/
getFloatArray(xmlChild, vParams, false, "", "Cphi");
nParamsFound = vParams.size();
size_t nParamsFound = vParams.size();
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
if (nParamsFound != 1) {
throw CanteraError("HMWSoln::readXMLBinarySalt::Cphi for " + ispName
@ -232,7 +231,7 @@ void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
+ "::" + jspName,
"wrong number of params found");
}
for (i = 0; i < nParamsFound; i++) {
for (size_t i = 0; i < nParamsFound; i++) {
m_CphiMX_ij_coeff(i, counter) = vParams[i];
}
m_CphiMX_ij[counter] = vParams[0];
@ -255,7 +254,6 @@ void HMWSoln::readXMLThetaAnion(XML_Node& BinSalt)
{
string xname = BinSalt.name();
vector_fp vParams;
size_t nParamsFound = 0;
if (xname != "thetaAnion") {
throw CanteraError("HMWSoln::readXMLThetaAnion",
"Incorrect name for processing this routine: " + xname);
@ -296,7 +294,7 @@ void HMWSoln::readXMLThetaAnion(XML_Node& BinSalt)
string nodeName = lowercase(stemp);
if (nodeName == "theta") {
getFloatArray(xmlChild, vParams, false, "", stemp);
nParamsFound = vParams.size();
size_t nParamsFound = vParams.size();
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
if (nParamsFound != 1) {
throw CanteraError("HMWSoln::readXMLThetaAnion::Theta for " + ispName
@ -336,12 +334,10 @@ void HMWSoln::readXMLThetaCation(XML_Node& BinSalt)
{
string xname = BinSalt.name();
vector_fp vParams;
size_t nParamsFound = 0;
if (xname != "thetaCation") {
throw CanteraError("HMWSoln::readXMLThetaCation",
"Incorrect name for processing this routine: " + xname);
}
string stemp;
string ispName = BinSalt.attrib("cation1");
if (ispName == "") {
throw CanteraError("HMWSoln::readXMLThetaCation", "no cation1 attrib");
@ -373,11 +369,11 @@ void HMWSoln::readXMLThetaCation(XML_Node& BinSalt)
int counter = m_CounterIJ[n];
for (size_t i = 0; i < BinSalt.nChildren(); i++) {
XML_Node& xmlChild = BinSalt.child(i);
stemp = xmlChild.name();
string stemp = xmlChild.name();
string nodeName = lowercase(stemp);
if (nodeName == "theta") {
getFloatArray(xmlChild, vParams, false, "", stemp);
nParamsFound = vParams.size();
size_t nParamsFound = vParams.size();
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
if (nParamsFound != 1) {
throw CanteraError("HMWSoln::readXMLThetaCation::Theta for " + ispName
@ -420,9 +416,7 @@ void HMWSoln::readXMLPsiCommonCation(XML_Node& BinSalt)
throw CanteraError("HMWSoln::readXMLPsiCommonCation",
"Incorrect name for processing this routine: " + xname);
}
string stemp;
vector_fp vParams;
size_t nParamsFound = 0;
string kName = BinSalt.attrib("cation");
if (kName == "") {
throw CanteraError("HMWSoln::readXMLPsiCommonCation", "no cation attrib");
@ -468,7 +462,7 @@ void HMWSoln::readXMLPsiCommonCation(XML_Node& BinSalt)
int counter = m_CounterIJ[n];
for (size_t i = 0; i < BinSalt.nChildren(); i++) {
XML_Node& xmlChild = BinSalt.child(i);
stemp = xmlChild.name();
string stemp = xmlChild.name();
string nodeName = lowercase(stemp);
if (nodeName == "theta") {
stemp = xmlChild.value();
@ -483,7 +477,7 @@ void HMWSoln::readXMLPsiCommonCation(XML_Node& BinSalt)
}
if (nodeName == "psi") {
getFloatArray(xmlChild, vParams, false, "", stemp);
nParamsFound = vParams.size();
size_t nParamsFound = vParams.size();
n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies ;
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
@ -560,9 +554,7 @@ void HMWSoln::readXMLPsiCommonAnion(XML_Node& BinSalt)
throw CanteraError("HMWSoln::readXMLPsiCommonAnion",
"Incorrect name for processing this routine: " + xname);
}
string stemp;
vector_fp vParams;
size_t nParamsFound = 0;
string kName = BinSalt.attrib("anion");
if (kName == "") {
throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "no anion attrib");
@ -607,7 +599,7 @@ void HMWSoln::readXMLPsiCommonAnion(XML_Node& BinSalt)
int counter = m_CounterIJ[n];
for (size_t i = 0; i < BinSalt.nChildren(); i++) {
XML_Node& xmlChild = BinSalt.child(i);
stemp = xmlChild.name();
string stemp = xmlChild.name();
string nodeName = lowercase(stemp);
if (nodeName == "theta") {
stemp = xmlChild.value();
@ -623,7 +615,7 @@ void HMWSoln::readXMLPsiCommonAnion(XML_Node& BinSalt)
if (nodeName == "psi") {
getFloatArray(xmlChild, vParams, false, "", stemp);
nParamsFound = vParams.size();
size_t nParamsFound = vParams.size();
n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies ;
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
@ -698,7 +690,6 @@ void HMWSoln::readXMLLambdaNeutral(XML_Node& BinSalt)
{
string xname = BinSalt.name();
vector_fp vParams;
size_t nParamsFound;
if (xname != "lambdaNeutral") {
throw CanteraError("HMWSoln::readXMLLanbdaNeutral",
"Incorrect name for processing this routine: " + xname);
@ -736,7 +727,7 @@ void HMWSoln::readXMLLambdaNeutral(XML_Node& BinSalt)
if (nodeName == "lambda") {
size_t nCount = iSpecies*m_kk + jSpecies;
getFloatArray(xmlChild, vParams, false, "", stemp);
nParamsFound = vParams.size();
size_t nParamsFound = vParams.size();
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
if (nParamsFound != 1) {
throw CanteraError("HMWSoln::readXMLLambdaNeutral::Lambda for " + iName
@ -778,7 +769,6 @@ void HMWSoln::readXMLMunnnNeutral(XML_Node& BinSalt)
{
string xname = BinSalt.name();
vector_fp vParams;
size_t nParamsFound;
if (xname != "MunnnNeutral") {
throw CanteraError("HMWSoln::readXMLMunnnNeutral",
"Incorrect name for processing this routine: " + xname);
@ -808,7 +798,7 @@ void HMWSoln::readXMLMunnnNeutral(XML_Node& BinSalt)
string nodeName = lowercase(stemp);
if (nodeName == "munnn") {
getFloatArray(xmlChild, vParams, false, "", "Munnn");
nParamsFound = vParams.size();
size_t nParamsFound = vParams.size();
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
if (nParamsFound != 1) {
throw CanteraError("HMWSoln::readXMLMunnnNeutral::Munnn for " + iName,
@ -850,9 +840,7 @@ void HMWSoln::readXMLZetaCation(const XML_Node& BinSalt)
throw CanteraError("HMWSoln::readXMLZetaCation",
"Incorrect name for processing this routine: " + xname);
}
string stemp;
vector_fp vParams;
size_t nParamsFound = 0;
string iName = BinSalt.attrib("neutral");
if (iName == "") {
@ -898,11 +886,11 @@ void HMWSoln::readXMLZetaCation(const XML_Node& BinSalt)
for (size_t i = 0; i < BinSalt.nChildren(); i++) {
XML_Node& xmlChild = BinSalt.child(i);
stemp = xmlChild.name();
string stemp = xmlChild.name();
string nodeName = lowercase(stemp);
if (nodeName == "zeta") {
getFloatArray(xmlChild, vParams, false, "", "zeta");
nParamsFound = vParams.size();
size_t nParamsFound = vParams.size();
size_t n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies ;
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
@ -1007,7 +995,6 @@ void HMWSoln::constructPhaseFile(std::string inputFile, std::string id_)
void HMWSoln::constructPhaseXML(XML_Node& phaseNode, std::string id_)
{
string stemp;
if (id_.size() > 0) {
string idp = phaseNode.id();
if (idp != id_) {
@ -1031,7 +1018,7 @@ void HMWSoln::constructPhaseXML(XML_Node& phaseNode, std::string id_)
if (thermoNode.hasChild("standardConc")) {
XML_Node& scNode = thermoNode.child("standardConc");
m_formGC = 2;
stemp = scNode.attrib("model");
string stemp = scNode.attrib("model");
string formString = lowercase(stemp);
if (formString != "") {
if (formString == "unity") {
@ -1072,7 +1059,7 @@ void HMWSoln::constructPhaseXML(XML_Node& phaseNode, std::string id_)
*/
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& scNode = thermoNode.child("activityCoefficients");
stemp = scNode.attrib("model");
string stemp = scNode.attrib("model");
string formString = lowercase(stemp);
if (formString != "") {
if (formString == "pitzer" || formString == "default") {
@ -1130,7 +1117,6 @@ void HMWSoln::constructPhaseXML(XML_Node& phaseNode, std::string id_)
void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
string stemp;
if (id_.size() > 0) {
string idp = phaseNode.id();
if (idp != id_) {
@ -1154,7 +1140,7 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
if (thermoNode.hasChild("standardConc")) {
XML_Node& scNode = thermoNode.child("standardConc");
m_formGC = 2;
stemp = scNode.attrib("model");
string stemp = scNode.attrib("model");
string formString = lowercase(stemp);
if (formString != "") {
if (formString == "unity") {
@ -1180,7 +1166,7 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
*/
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& scNode = thermoNode.child("activityCoefficients");
stemp = scNode.attrib("model");
string stemp = scNode.attrib("model");
string formString = lowercase(stemp);
if (formString != "") {
if (formString == "pitzer" || formString == "default") {
@ -1388,7 +1374,7 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
if (acNode.hasChild("A_Debye")) {
XML_Node& ADebye = acNode.child("A_Debye");
m_form_A_Debye = A_DEBYE_CONST;
stemp = "model";
string stemp = "model";
if (ADebye.hasAttrib(stemp)) {
string atemp = ADebye.attrib(stemp);
stemp = lowercase(atemp);
@ -1445,14 +1431,12 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
*/
std::vector<const XML_Node*> xspecies = speciesData();
string kname, jname;
size_t jj = xspecies.size();
for (size_t k = 0; k < m_kk; k++) {
size_t jmap = npos;
kname = speciesName(k);
for (size_t j = 0; j < jj; j++) {
string kname = speciesName(k);
for (size_t j = 0; j < xspecies.size(); j++) {
const XML_Node& sp = *xspecies[j];
jname = sp["name"];
string jname = sp["name"];
if (jname == kname) {
jmap = j;
break;
@ -1473,8 +1457,9 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
map<string, string> msIs;
getMap(sIsNode, msIs);
map<string,string>::const_iterator _b = msIs.begin();
for (; _b != msIs.end(); ++_b) {
for (map<string,string>::const_iterator _b = msIs.begin();
_b != msIs.end();
++_b) {
size_t kk = speciesIndex(_b->first);
if (kk != npos) {
double val = fpValue(_b->second);
@ -1492,7 +1477,7 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
if (acNodePtr) {
for (size_t i = 0; i < acNodePtr->nChildren(); i++) {
XML_Node& xmlACChild = acNodePtr->child(i);
stemp = xmlACChild.name();
string stemp = xmlACChild.name();
string nodeName = lowercase(stemp);
/*
* Process a binary salt field, or any of the other XML fields
@ -1549,11 +1534,8 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
* in each of the species SS databases.
*/
std::vector<const XML_Node*> xspecies = speciesData();
const XML_Node* spPtr = 0;
string kname;
for (size_t k = 0; k < m_kk; k++) {
kname = speciesName(k);
spPtr = xspecies[k];
const XML_Node* spPtr = xspecies[k];
if (!spPtr) {
if (spPtr->hasChild("electrolyteSpeciesType")) {
string est = getChildValue(*spPtr, "electrolyteSpeciesType");
@ -1572,8 +1554,9 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType");
map<string, string> msEST;
getMap(ESTNode, msEST);
map<string,string>::const_iterator _b = msEST.begin();
for (; _b != msEST.end(); ++_b) {
for (map<string,string>::const_iterator _b = msEST.begin();
_b != msEST.end();
++_b) {
size_t kk = speciesIndex(_b->first);
if (kk != npos) {
string est = _b->second;
@ -1670,8 +1653,7 @@ void HMWSoln::calcIMSCutoffParams_()
IMS_efCut_ = 0.0;
bool converged = false;
double oldV = 0.0;
int its;
for (its = 0; its < 100 && !converged; its++) {
for (int its = 0; its < 100 && !converged; its++) {
oldV = IMS_efCut_;
IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_) -IMS_efCut_;
IMS_bfCut_ = IMS_afCut_ / IMS_cCut_ + IMS_slopefCut_ - 1.0;
@ -1693,7 +1675,7 @@ void HMWSoln::calcIMSCutoffParams_()
double f_0 = IMS_afCut_ + IMS_efCut_;
double f_prime_0 = 1.0 - IMS_afCut_ / IMS_cCut_ + IMS_bfCut_;
IMS_egCut_ = 0.0;
for (its = 0; its < 100 && !converged; its++) {
for (int its = 0; its < 100 && !converged; its++) {
oldV = IMS_egCut_;
double lng_0 = -log(IMS_gamma_o_min_) - f_prime_0 / f_0;
IMS_agCut_ = exp(lng_0) - IMS_egCut_;
@ -1726,9 +1708,8 @@ void HMWSoln::calcMCCutoffParams_()
MC_epCut_ = 0.0;
bool converged = false;
double oldV = 0.0;
int its;
double damp = 0.5;
for (its = 0; its < 500 && !converged; its++) {
for (int its = 0; its < 500 && !converged; its++) {
oldV = MC_epCut_;
MC_apCut_ = damp *(MC_X_o_min_ - MC_epCut_) + (1-damp) * MC_apCut_;
double MC_bpCutNew = MC_apCut_ / MC_cpCut_ + MC_slopepCut_ - 1.0;

16
src/thermo/IdealGasPhase.cpp
View file

@ -139,8 +139,7 @@ doublereal IdealGasPhase::cv_vib(const int k, const doublereal T) const
doublereal IdealGasPhase::standardConcentration(size_t k) const
{
double p = pressure();
return p / (GasConstant * temperature());
return pressure() / (GasConstant * temperature());
}
void IdealGasPhase::getActivityCoefficients(doublereal* ac) const
@ -166,10 +165,9 @@ void IdealGasPhase::getStandardChemPotentials(doublereal* muStar) const
void IdealGasPhase::getChemPotentials(doublereal* mu) const
{
getStandardChemPotentials(mu);
doublereal xx;
doublereal rt = temperature() * GasConstant;
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(SmallNumber, moleFraction(k));
double xx = std::max(SmallNumber, moleFraction(k));
mu[k] += rt * (log(xx));
}
}
@ -184,12 +182,11 @@ void IdealGasPhase::getPartialMolarEnthalpies(doublereal* hbar) const
void IdealGasPhase::getPartialMolarEntropies(doublereal* sbar) const
{
const vector_fp& _s = entropy_R_ref();
doublereal r = GasConstant;
scale(_s.begin(), _s.end(), sbar, r);
scale(_s.begin(), _s.end(), sbar, GasConstant);
doublereal logp = log(pressure() / m_spthermo->refPressure());
for (size_t k = 0; k < m_kk; k++) {
doublereal xx = std::max(SmallNumber, moleFraction(k));
sbar[k] += r * (-logp - log(xx));
sbar[k] += GasConstant * (-logp - log(xx));
}
}
@ -339,7 +336,6 @@ void IdealGasPhase::initThermo()
void IdealGasPhase::setToEquilState(const doublereal* mu_RT)
{
double tmp, tmp2;
const vector_fp& grt = gibbs_RT_ref();
/*
@ -352,11 +348,11 @@ void IdealGasPhase::setToEquilState(const doublereal* mu_RT)
*/
doublereal pres = 0.0;
for (size_t k = 0; k < m_kk; k++) {
tmp = -grt[k] + mu_RT[k];
double tmp = -grt[k] + mu_RT[k];
if (tmp < -600.) {
m_pp[k] = 0.0;
} else if (tmp > 300.0) {
tmp2 = tmp / 300.;
double tmp2 = tmp / 300.;
tmp2 *= tmp2;
m_pp[k] = m_p0 * exp(300.) * tmp2;
} else {

59
src/thermo/IdealMolalSoln.cpp
View file

@ -206,8 +206,7 @@ doublereal IdealMolalSoln::thermalExpansionCoeff() const
void IdealMolalSoln::setDensity(const doublereal rho)
{
double dens = density();
if (rho != dens) {
if (rho != density()) {
throw CanteraError("Idea;MolalSoln::setDensity",
"Density is not an independent variable");
}
@ -215,8 +214,7 @@ void IdealMolalSoln::setDensity(const doublereal rho)
void IdealMolalSoln::setMolarDensity(const doublereal conc)
{
double concI = Phase::molarDensity();
if (conc != concI) {
if (conc != Phase::molarDensity()) {
throw CanteraError("IdealMolalSoln::setMolarDensity",
"molarDensity/denisty is not an independent variable");
}
@ -253,16 +251,15 @@ void IdealMolalSoln::getActivityConcentrations(doublereal* c) const
doublereal IdealMolalSoln::standardConcentration(size_t k) const
{
double c0 = 1.0, mvSolvent;
double c0 = 1.0;
switch (m_formGC) {
case 0:
break;
case 1:
c0 = 1.0 /m_speciesMolarVolume[m_indexSolvent];
return c0 = 1.0 /m_speciesMolarVolume[m_indexSolvent];
break;
case 2:
mvSolvent = m_speciesMolarVolume[m_indexSolvent];
c0 = 1.0 / mvSolvent;
c0 = 1.0 / m_speciesMolarVolume[m_indexSolvent];
break;
}
return c0;
@ -360,8 +357,6 @@ void IdealMolalSoln::getMolalityActivityCoefficients(doublereal* acMolality) con
void IdealMolalSoln::getChemPotentials(doublereal* mu) const
{
double xx;
// Assertion is made for speed
AssertThrow(m_indexSolvent == 0, "solvent not the first species");
@ -386,7 +381,7 @@ void IdealMolalSoln::getChemPotentials(doublereal* mu) const
if (IMS_typeCutoff_ == 0 || xmolSolvent > 3.* IMS_X_o_cutoff_/2.0) {
for (size_t k = 1; k < m_kk; k++) {
xx = std::max(m_molalities[k], SmallNumber);
double xx = std::max(m_molalities[k], SmallNumber);
mu[k] += RT * log(xx);
}
/*
@ -394,7 +389,7 @@ void IdealMolalSoln::getChemPotentials(doublereal* mu) const
* -> see my notes
*/
xx = std::max(xmolSolvent, SmallNumber);
double xx = std::max(xmolSolvent, SmallNumber);
mu[m_indexSolvent] +=
(RT * (xmolSolvent - 1.0) / xx);
} else {
@ -406,10 +401,10 @@ void IdealMolalSoln::getChemPotentials(doublereal* mu) const
for (size_t k = 1; k < m_kk; k++) {
xx = std::max(m_molalities[k], SmallNumber);
double xx = std::max(m_molalities[k], SmallNumber);
mu[k] += RT * (log(xx) + IMS_lnActCoeffMolal_[k]);
}
xx = std::max(xmolSolvent, SmallNumber);
double xx = std::max(xmolSolvent, SmallNumber);
mu[m_indexSolvent] +=
RT * (log(xx) + IMS_lnActCoeffMolal_[m_indexSolvent]);
}
@ -428,17 +423,16 @@ void IdealMolalSoln::getPartialMolarEnthalpies(doublereal* hbar) const
void IdealMolalSoln::getPartialMolarEntropies(doublereal* sbar) const
{
getEntropy_R(sbar);
doublereal R = GasConstant;
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] -= R * log(mm);
sbar[k] -= GasConstant * log(mm);
}
}
double xmolSolvent = moleFraction(m_indexSolvent);
sbar[m_indexSolvent] -= (R * (xmolSolvent - 1.0) / xmolSolvent);
sbar[m_indexSolvent] -= (GasConstant * (xmolSolvent - 1.0) / xmolSolvent);
} else {
/*
* Update the activity coefficients, This also update the
@ -453,12 +447,12 @@ void IdealMolalSoln::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) + IMS_lnActCoeffMolal_[k]);
sbar[k] -= GasConstant * (log(mm) + IMS_lnActCoeffMolal_[k]);
}
}
double xmolSolvent = moleFraction(m_indexSolvent);
mm = std::max(SmallNumber, xmolSolvent);
sbar[m_indexSolvent] -= R *(log(mm) + IMS_lnActCoeffMolal_[m_indexSolvent]);
sbar[m_indexSolvent] -= GasConstant *(log(mm) + IMS_lnActCoeffMolal_[m_indexSolvent]);
}
}
@ -507,8 +501,7 @@ void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
initThermo();
if (id_.size() > 0) {
std::string idp = phaseNode.id();
if (idp != id_) {
if (phaseNode.id() != id_) {
throw CanteraError("IdealMolalSoln::initThermo",
"phasenode and Id are incompatible");
}
@ -550,9 +543,8 @@ void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
*/
std::string solventName = "";
if (thermoNode.hasChild("solvent")) {
XML_Node& scNode = thermoNode.child("solvent");
std::vector<std::string> nameSolventa;
getStringArray(scNode, nameSolventa);
getStringArray(thermoNode.child("solvent"), nameSolventa);
if (nameSolventa.size() != 1) {
throw CanteraError("IdealMolalSoln::initThermoXML",
"badly formed solvent XML node");
@ -609,8 +601,7 @@ void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
* Reconcile the solvent name and index.
*/
for (size_t k = 0; k < m_kk; k++) {
std::string sname = speciesName(k);
if (solventName == sname) {
if (solventName == speciesName(k)) {
m_indexSolvent = k;
break;
}
@ -668,7 +659,6 @@ void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
{
double tmp;
/*
* Calculate the molalities. Currently, the molalities
* may not be current with respect to the contents of the
@ -693,7 +683,7 @@ void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
return;
} else if (xmolSolvent < IMS_X_o_cutoff_/2.0) {
tmp = log(xx * IMS_gamma_k_min_);
double tmp = log(xx * IMS_gamma_k_min_);
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
@ -729,7 +719,7 @@ void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
double g = h2 + h1_g;
double g_prime = h2_prime + h1_g_prime;
tmp = (xmolSolvent/ g * g_prime + (1.0-xmolSolvent) / f * f_prime);
double tmp = (xmolSolvent/ g * g_prime + (1.0-xmolSolvent) / f * f_prime);
double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
@ -764,7 +754,7 @@ void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
double g_prime = 1.0 + eterm*gptmp;
double g = xmolSolvent + IMS_egCut_ + eterm * (IMS_agCut_ + xmolSolvent * (IMS_bgCut_ + IMS_dgCut_*xmolSolvent));
tmp = (xmolSolvent / g * g_prime + (1.0 - xmolSolvent) / f * f_prime);
double tmp = (xmolSolvent / g * g_prime + (1.0 - xmolSolvent) / f * f_prime);
double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
@ -775,7 +765,6 @@ void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
}
}
return;
}
void IdealMolalSoln::initLengths()
@ -796,10 +785,8 @@ void IdealMolalSoln::calcIMSCutoffParams_()
IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_);
IMS_efCut_ = 0.0;
bool converged = false;
double oldV = 0.0;
int its;
for (its = 0; its < 100 && !converged; its++) {
oldV = IMS_efCut_;
for (int its = 0; its < 100 && !converged; its++) {
double oldV = IMS_efCut_;
IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_) - IMS_efCut_;
IMS_bfCut_ = IMS_afCut_ / IMS_cCut_ + IMS_slopefCut_ - 1.0;
IMS_dfCut_ = ((- IMS_afCut_/IMS_cCut_ + IMS_bfCut_ - IMS_bfCut_*IMS_X_o_cutoff_/IMS_cCut_)
@ -820,8 +807,8 @@ void IdealMolalSoln::calcIMSCutoffParams_()
double f_0 = IMS_afCut_ + IMS_efCut_;
double f_prime_0 = 1.0 - IMS_afCut_ / IMS_cCut_ + IMS_bfCut_;
IMS_egCut_ = 0.0;
for (its = 0; its < 100 && !converged; its++) {
oldV = IMS_egCut_;
for (int its = 0; its < 100 && !converged; its++) {
double oldV = IMS_egCut_;
double lng_0 = -log(IMS_gamma_o_min_) - f_prime_0 / f_0;
IMS_agCut_ = exp(lng_0) - IMS_egCut_;
IMS_bgCut_ = IMS_agCut_ / IMS_cCut_ + IMS_slopegCut_ - 1.0;

57
src/thermo/IdealSolidSolnPhase.cpp
View file

@ -128,8 +128,7 @@ doublereal IdealSolidSolnPhase::entropy_mole() const
doublereal IdealSolidSolnPhase::gibbs_mole() const
{
const double* dptr = DATA_PTR(gibbs_RT_ref());
doublereal g = mean_X(dptr);
return GasConstant * temperature() * (g + sum_xlogx());
return GasConstant * temperature() * (mean_X(dptr) + sum_xlogx());
}
doublereal IdealSolidSolnPhase::cp_mole() const
@ -154,8 +153,7 @@ void IdealSolidSolnPhase::calcDensity()
* Set the density in the parent State object directly,
* by calling the Phase::setDensity() function.
*/
double dens = 1.0/invDens;
Phase::setDensity(dens);
Phase::setDensity(1.0/invDens);
}
void IdealSolidSolnPhase::setDensity(const doublereal rho)
@ -166,8 +164,7 @@ void IdealSolidSolnPhase::setDensity(const doublereal rho)
* exception. This is because the density is NOT an
* independent variable.
*/
double dens = density();
if (rho != dens) {
if (rho != density()) {
throw CanteraError("IdealSolidSolnPhase::setDensity",
"Density is not an independent variable");
}
@ -330,11 +327,10 @@ void IdealSolidSolnPhase::getActivityCoefficients(doublereal* ac) const
void IdealSolidSolnPhase::getChemPotentials(doublereal* mu) const
{
doublereal delta_p = m_Pcurrent - m_Pref;
doublereal xx;
doublereal RT = temperature() * GasConstant;
const vector_fp& g_RT = gibbs_RT_ref();
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(SmallNumber, moleFraction(k));
double xx = std::max(SmallNumber, moleFraction(k));
mu[k] = RT * (g_RT[k] + log(xx))
+ delta_p * m_speciesMolarVolume[k];
}
@ -342,12 +338,10 @@ void IdealSolidSolnPhase::getChemPotentials(doublereal* mu) const
void IdealSolidSolnPhase::getChemPotentials_RT(doublereal* mu) const
{
doublereal RT = temperature() * GasConstant;
doublereal delta_pdRT = (m_Pcurrent - m_Pref) / RT;
doublereal xx;
doublereal delta_pdRT = (m_Pcurrent - m_Pref) / (temperature() * GasConstant);
const vector_fp& g_RT = gibbs_RT_ref();
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(SmallNumber, moleFraction(k));
double xx = std::max(SmallNumber, moleFraction(k));
mu[k] = (g_RT[k] + log(xx))
+ delta_pdRT * m_speciesMolarVolume[k];
}
@ -360,18 +354,15 @@ void IdealSolidSolnPhase::getChemPotentials_RT(doublereal* mu) const
void IdealSolidSolnPhase::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, GasConstant * temperature());
}
void IdealSolidSolnPhase::getPartialMolarEntropies(doublereal* sbar) const
{
const vector_fp& _s = entropy_R_ref();
doublereal r = GasConstant;
doublereal xx;
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(SmallNumber, moleFraction(k));
sbar[k] = r * (_s[k] - log(xx));
double xx = std::max(SmallNumber, moleFraction(k));
sbar[k] = GasConstant * (_s[k] - log(xx));
}
}
@ -522,11 +513,9 @@ const vector_fp& IdealSolidSolnPhase::entropy_R_ref() const
void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
string subname = "IdealSolidSolnPhase::initThermoXML";
if (id_.size() > 0) {
string idp = phaseNode.id();
if (idp != id_) {
throw CanteraError(subname.c_str(),
if (phaseNode.id() != id_) {
throw CanteraError("IdealSolidSolnPhase::initThermoXML",
"phasenode and Id are incompatible");
}
}
@ -537,14 +526,13 @@ void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string&
*/
if (phaseNode.hasChild("thermo")) {
XML_Node& thNode = phaseNode.child("thermo");
string mStringa = thNode.attrib("model");
string mString = lowercase(mStringa);
if (mString != "idealsolidsolution") {
throw CanteraError(subname.c_str(),
"Unknown thermo model: " + mStringa);
string mString = thNode.attrib("model");
if (lowercase(mString) != "idealsolidsolution") {
throw CanteraError("IdealSolidSolnPhase::initThermoXML",
"Unknown thermo model: " + mString);
}
} else {
throw CanteraError(subname.c_str(),
throw CanteraError("IdealSolidSolnPhase::initThermoXML",
"Unspecified thermo model");
}
@ -566,11 +554,11 @@ void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string&
} else if (formString == "solvent_volume") {
m_formGC = 2;
} else {
throw CanteraError(subname.c_str(),
throw CanteraError("IdealSolidSolnPhase::initThermoXML",
"Unknown standardConc model: " + formStringa);
}
} else {
throw CanteraError(subname.c_str(),
throw CanteraError("IdealSolidSolnPhase::initThermoXML",
"Unspecified standardConc model");
}
@ -585,10 +573,9 @@ void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string&
XML_Node& speciesList = phaseNode.child("speciesArray");
XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"],
&phaseNode.root());
const vector<string>&sss = speciesNames();
for (size_t k = 0; k < m_kk; k++) {
XML_Node* s = speciesDB->findByAttr("name", sss[k]);
XML_Node* s = speciesDB->findByAttr("name", speciesName(k));
XML_Node* ss = s->findByName("standardState");
m_speciesMolarVolume[k] = ctml::getFloat(*ss, "molarVolume", "toSI");
}
@ -635,8 +622,7 @@ void IdealSolidSolnPhase::setToEquilState(const doublereal* lambda_RT)
m_pp[k] = m_Pref * exp(m_pp[k]);
pres += m_pp[k];
}
doublereal* dptr = DATA_PTR(m_pp);
setState_PX(pres, dptr);
setState_PX(pres, &m_pp[0]);
}
double IdealSolidSolnPhase::speciesMolarVolume(int k) const
@ -660,9 +646,8 @@ void IdealSolidSolnPhase::_updateThermo() const
DATA_PTR(m_s0_R));
m_tlast = tnow;
doublereal rrt = 1.0 / (GasConstant * tnow);
doublereal deltaE;
for (size_t k = 0; k < m_kk; k++) {
deltaE = rrt * m_pe[k];
double deltaE = rrt * m_pe[k];
m_h0_RT[k] += deltaE;
m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
}

28
src/thermo/IdealSolnGasVPSS.cpp
View file

@ -133,12 +133,11 @@ void IdealSolnGasVPSS::calcDensity()
} else {
const doublereal* const dtmp = moleFractdivMMW();
const vector_fp& vss = m_VPSS_ptr->getStandardVolumes();
double invDens = dot(vss.begin(), vss.end(), dtmp);
double dens = 1.0 / dot(vss.begin(), vss.end(), dtmp);
/*
* Set the density in the parent State object directly,
* by calling the Phase::setDensity() function.
*/
double dens = 1.0/invDens;
Phase::setDensity(dens);
}
}
@ -202,8 +201,7 @@ doublereal IdealSolnGasVPSS::standardConcentration(size_t k) const
void IdealSolnGasVPSS::getUnitsStandardConc(double* uA, int, int sizeUA) const
{
int eos = eosType();
if (eos == cIdealSolnGasPhase0) {
if (eosType() == cIdealSolnGasPhase0) {
for (int i = 0; i < sizeUA; i++) {
uA[i] = 0.0;
}
@ -254,10 +252,9 @@ void IdealSolnGasVPSS::getChemPotentials_RT(doublereal* muRT) const
void IdealSolnGasVPSS::getChemPotentials(doublereal* mu) const
{
getStandardChemPotentials(mu);
doublereal xx;
doublereal rt = temperature() * GasConstant;
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(SmallNumber, moleFraction(k));
double xx = std::max(SmallNumber, moleFraction(k));
mu[k] += rt*(log(xx));
}
}
@ -265,33 +262,29 @@ void IdealSolnGasVPSS::getChemPotentials(doublereal* mu) const
void IdealSolnGasVPSS::getPartialMolarEnthalpies(doublereal* hbar) const
{
getEnthalpy_RT(hbar);
doublereal rt = GasConstant * temperature();
scale(hbar, hbar+m_kk, hbar, rt);
scale(hbar, hbar+m_kk, hbar, GasConstant * temperature());
}
void IdealSolnGasVPSS::getPartialMolarEntropies(doublereal* sbar) const
{
getEntropy_R(sbar);
doublereal r = GasConstant;
scale(sbar, sbar+m_kk, sbar, r);
scale(sbar, sbar+m_kk, sbar, GasConstant);
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));
}
}
void IdealSolnGasVPSS::getPartialMolarIntEnergies(doublereal* ubar) const
{
getIntEnergy_RT(ubar);
doublereal rt = GasConstant * temperature();
scale(ubar, ubar+m_kk, ubar, rt);
scale(ubar, ubar+m_kk, ubar, GasConstant * temperature());
}
void IdealSolnGasVPSS::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 IdealSolnGasVPSS::getPartialMolarVolumes(doublereal* vbar) const
@ -307,7 +300,6 @@ void IdealSolnGasVPSS::initThermo()
void IdealSolnGasVPSS::setToEquilState(const doublereal* mu_RT)
{
double tmp, tmp2;
updateStandardStateThermo();
const vector_fp& grt = m_VPSS_ptr->Gibbs_RT_ref();
@ -322,11 +314,11 @@ void IdealSolnGasVPSS::setToEquilState(const doublereal* mu_RT)
doublereal pres = 0.0;
double m_p0 = m_VPSS_ptr->refPressure();
for (size_t k = 0; k < m_kk; k++) {
tmp = -grt[k] + mu_RT[k];
double tmp = -grt[k] + mu_RT[k];
if (tmp < -600.) {
m_pp[k] = 0.0;
} else if (tmp > 500.0) {
tmp2 = tmp / 500.;
double tmp2 = tmp / 500.;
tmp2 *= tmp2;
m_pp[k] = m_p0 * exp(500.) * tmp2;
} else {

71
src/thermo/IonsFromNeutralVPSSTP.cpp
View file

@ -191,9 +191,9 @@ void IonsFromNeutralVPSSTP::constructPhaseFile(std::string inputFile, std::strin
* The phase object automatically constructs an XML object.
* Use this object to store information.
*/
XML_Node* fxml = new XML_Node();
fxml->build(fin);
XML_Node* fxml_phase = findXMLPhase(fxml, id_);
XML_Node fxml;
fxml.build(fin);
XML_Node* fxml_phase = findXMLPhase(&fxml, id_);
if (!fxml_phase) {
throw CanteraError("MargulesVPSSTP:constructPhaseFile",
"ERROR: Can not find phase named " +
@ -201,15 +201,12 @@ void IonsFromNeutralVPSSTP::constructPhaseFile(std::string inputFile, std::strin
}
setXMLdata(*fxml_phase);
constructPhaseXML(*fxml_phase, id_);
delete fxml;
}
void IonsFromNeutralVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string id_)
{
string stemp;
if (id_.size() > 0) {
string idp = phaseNode.id();
if (idp != id_) {
if (phaseNode.id() != id_) {
throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
"phasenode and Id are incompatible");
}
@ -229,8 +226,7 @@ void IonsFromNeutralVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string i
/*
* Make sure that the thermo model is IonsFromNeutralMolecule
*/
stemp = thermoNode.attrib("model");
string formString = lowercase(stemp);
string formString = lowercase(thermoNode.attrib("model"));
if (formString != "ionsfromneutralmolecule") {
throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
"model name isn't IonsFromNeutralMolecule: " + formString);
@ -245,8 +241,7 @@ void IonsFromNeutralVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string i
}
XML_Node& neutralMoleculeNode = thermoNode.child("neutralMoleculePhase");
string nsource = neutralMoleculeNode["datasrc"];
XML_Node* neut_ptr = get_XML_Node(nsource, 0);
XML_Node* neut_ptr = get_XML_Node(neutralMoleculeNode["datasrc"], 0);
if (!neut_ptr) {
throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
"neut_ptr = 0");
@ -427,7 +422,6 @@ void IonsFromNeutralVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const
void IonsFromNeutralVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
{
double xx;
/*
* Get the nondimensional standard state entropies
*/
@ -441,7 +435,7 @@ void IonsFromNeutralVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
s_update_dlnActCoeffdT();
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(moleFractions_[k], SmallNumber);
double xx = std::max(moleFractions_[k], SmallNumber);
sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
}
/*
@ -486,14 +480,12 @@ void IonsFromNeutralVPSSTP::getdlnActCoeffdlnN(const size_t ld, doublereal* dlnA
void IonsFromNeutralVPSSTP::setTemperature(const doublereal temp)
{
double p = pressure();
IonsFromNeutralVPSSTP::setState_TP(temp, p);
IonsFromNeutralVPSSTP::setState_TP(temp, pressure());
}
void IonsFromNeutralVPSSTP::setPressure(doublereal p)
{
double t = temperature();
IonsFromNeutralVPSSTP::setState_TP(t, p);
IonsFromNeutralVPSSTP::setState_TP(temperature(), p);
}
void IonsFromNeutralVPSSTP::setState_TP(doublereal t, doublereal p)
@ -508,13 +500,11 @@ void IonsFromNeutralVPSSTP::setState_TP(doublereal t, doublereal p)
/*
* Calculate the partial molar volumes, and then the density of the fluid
*/
double dd = neutralMoleculePhase_->density();
Phase::setDensity(dd);
Phase::setDensity(neutralMoleculePhase_->density());
}
void IonsFromNeutralVPSSTP::calcIonMoleFractions(doublereal* const mf) const
{
doublereal fmij;
/*
* Download the neutral mole fraction vector into the
* vector, NeutralMolecMoleFractions_[]
@ -531,7 +521,7 @@ void IonsFromNeutralVPSSTP::calcIonMoleFractions(doublereal* const mf) const
*/
for (size_t jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
for (size_t k = 0; k < m_kk; k++) {
fmij = fm_neutralMolec_ions_[k + jNeut * m_kk];
double fmij = fm_neutralMolec_ions_[k + jNeut * m_kk];
mf[k] += fmij * NeutralMolecMoleFractions_[jNeut];
}
}
@ -658,7 +648,6 @@ void IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions() const
void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal* const dx, doublereal* const dy) const
{
doublereal fmij;
doublereal sumy, sumdy;
//check sum dx = 0
@ -683,7 +672,7 @@ void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal* const
size_t icat = cationList_[k];
size_t jNeut = fm_invert_ionForNeutral[icat];
if (jNeut != npos) {
fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
double fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
AssertTrace(fmij != 0.0);
const doublereal temp = 1.0/fmij;
dy[jNeut] += dx[icat] * temp;
@ -694,7 +683,7 @@ void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal* const
for (size_t k = 0; k < numPassThroughSpecies_; k++) {
size_t icat = passThroughList_[k];
size_t jNeut = fm_invert_ionForNeutral[icat];
fmij = fm_neutralMolec_ions_[ icat + jNeut * m_kk];
double fmij = fm_neutralMolec_ions_[ icat + jNeut * m_kk];
const doublereal temp = 1.0/fmij;
dy[jNeut] += dx[icat] * temp;
y_[jNeut] += moleFractions_[icat] * temp;
@ -863,10 +852,8 @@ static double factorOverlap(const std::vector<std::string>& elnamesVN ,
}
void IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
string stemp;
if (id_.size() > 0) {
string idp = phaseNode.id();
if (idp != id_) {
if (phaseNode.id() != id_) {
throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
"phasenode and Id are incompatible");
}
@ -886,8 +873,7 @@ void IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string
/*
* Make sure that the thermo model is IonsFromNeutralMolecule
*/
stemp = thermoNode.attrib("model");
string formString = lowercase(stemp);
string formString = lowercase(thermoNode.attrib("model"));
if (formString != "ionsfromneutralmolecule") {
throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
"model name isn't IonsFromNeutralMolecule: " + formString);
@ -902,8 +888,7 @@ void IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string
}
XML_Node& neutralMoleculeNode = thermoNode.child("neutralMoleculePhase");
string nsource = neutralMoleculeNode["datasrc"];
XML_Node* neut_ptr = get_XML_Node(nsource, 0);
XML_Node* neut_ptr = get_XML_Node(neutralMoleculeNode["datasrc"], 0);
if (!neut_ptr) {
throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
"neut_ptr = 0");
@ -950,10 +935,10 @@ void IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string
}
}
PDSS_IonsFromNeutral* speciesSS = 0;
indexSpecialSpecies_ = npos;
for (size_t k = 0; k < m_kk; k++) {
speciesSS = dynamic_cast<PDSS_IonsFromNeutral*>(providePDSS(k));
PDSS_IonsFromNeutral* speciesSS =
dynamic_cast<PDSS_IonsFromNeutral*>(providePDSS(k));
if (!speciesSS) {
throw CanteraError("initThermoXML", "Dynamic cast failed");
}
@ -1063,7 +1048,6 @@ void IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string
void IonsFromNeutralVPSSTP::s_update_lnActCoeff() const
{
size_t icat, jNeut;
doublereal fmij;
/*
* Get the activity coefficiens of the neutral molecules
*/
@ -1079,7 +1063,7 @@ void IonsFromNeutralVPSSTP::s_update_lnActCoeff() const
//! Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
double fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
lnActCoeff_Scaled_[icat] = lnActCoeff_NeutralMolecule_[jNeut] / fmij;
}
@ -1113,7 +1097,6 @@ void IonsFromNeutralVPSSTP::getdlnActCoeffds(const doublereal dTds, const double
doublereal* dlnActCoeffds) const
{
size_t icat, jNeut;
doublereal fmij;
/*
* Get the activity coefficients of the neutral molecules
*/
@ -1144,7 +1127,7 @@ void IonsFromNeutralVPSSTP::getdlnActCoeffds(const doublereal dTds, const double
//! Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
double fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
dlnActCoeffds[icat] = dlnActCoeff_NeutralMolecule_[jNeut]/fmij;
}
@ -1177,7 +1160,6 @@ void IonsFromNeutralVPSSTP::getdlnActCoeffds(const doublereal dTds, const double
void IonsFromNeutralVPSSTP::s_update_dlnActCoeffdT() const
{
size_t icat, jNeut;
doublereal fmij;
/*
* Get the activity coefficients of the neutral molecules
*/
@ -1198,7 +1180,7 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeffdT() const
//! Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
double fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
dlnActCoeffdT_Scaled_[icat] = dlnActCoeffdT_NeutralMolecule_[jNeut]/fmij;
}
@ -1231,7 +1213,6 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeffdT() const
void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnX_diag() const
{
size_t icat, jNeut;
doublereal fmij;
/*
* Get the activity coefficients of the neutral molecules
*/
@ -1252,7 +1233,7 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnX_diag() const
//! Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
double fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
dlnActCoeffdlnX_diag_[icat] = dlnActCoeffdlnX_diag_NeutralMolecule_[jNeut]/fmij;
}
@ -1285,7 +1266,6 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnX_diag() const
void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN_diag() const
{
size_t icat, jNeut;
doublereal fmij;
/*
* Get the activity coefficients of the neutral molecules
*/
@ -1306,7 +1286,7 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN_diag() const
//! Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
double fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
dlnActCoeffdlnN_diag_[icat] = dlnActCoeffdlnN_diag_NeutralMolecule_[jNeut]/fmij;
}
@ -1339,7 +1319,7 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN_diag() const
void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN() const
{
size_t kcat = 0, kNeut = 0, mcat = 0, mNeut = 0;
doublereal fmij = 0.0, mfmij;
doublereal fmij = 0.0;
dlnActCoeffdlnN_.zero();
/*
* Get the activity coefficients of the neutral molecules
@ -1366,7 +1346,7 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN() const
mcat = cationList_[m];
mNeut = fm_invert_ionForNeutral[mcat];
mfmij = fm_neutralMolec_ions_[mcat + mNeut * m_kk];
double mfmij = fm_neutralMolec_ions_[mcat + mNeut * m_kk];
dlnActCoeffdlnN_(kcat,mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut,mNeut) * mfmij / fmij;
@ -1402,7 +1382,6 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN() const
for (size_t m = 0; m < cationList_.size(); m++) {
mcat = cationList_[m];
mNeut = fm_invert_ionForNeutral[mcat];
mfmij = fm_neutralMolec_ions_[mcat + mNeut * m_kk];
dlnActCoeffdlnN_(kcat, mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut,mNeut);
}

52
src/thermo/LatticePhase.cpp
View file

@ -67,16 +67,14 @@ ThermoPhase* LatticePhase::duplMyselfAsThermoPhase() const
doublereal LatticePhase::enthalpy_mole() const
{
doublereal p0 = m_spthermo->refPressure();
return GasConstant * temperature() *
mean_X(&enthalpy_RT_ref()[0])
+ (pressure() - p0)/molarDensity();
+ (pressure() - m_Pref)/molarDensity();
}
doublereal LatticePhase::entropy_mole() const
{
return GasConstant * (mean_X(&entropy_R_ref()[0]) -
sum_xlogx());
return GasConstant * (mean_X(&entropy_R_ref()[0]) - sum_xlogx());
}
doublereal LatticePhase::cp_mole() const
@ -92,9 +90,7 @@ doublereal LatticePhase::cv_mole() const
doublereal LatticePhase::calcDensity()
{
setMolarDensity(m_site_density);
doublereal mw = meanMolecularWeight();
doublereal dens = mw * m_site_density;
return dens;
return meanMolecularWeight() * m_site_density;
}
void LatticePhase::setPressure(doublereal p)
@ -158,11 +154,10 @@ doublereal LatticePhase::logStandardConc(size_t k) const
void LatticePhase::getChemPotentials(doublereal* mu) const
{
doublereal delta_p = m_Pcurrent - m_Pref;
doublereal xx;
doublereal RT = temperature() * GasConstant;
const vector_fp& g_RT = gibbs_RT_ref();
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(SmallNumber, moleFraction(k));
double xx = std::max(SmallNumber, moleFraction(k));
mu[k] = RT * (g_RT[k] + log(xx))
+ delta_p * m_speciesMolarVolume[k];
}
@ -172,18 +167,15 @@ void LatticePhase::getChemPotentials(doublereal* mu) const
void LatticePhase::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, GasConstant * temperature());
}
void LatticePhase::getPartialMolarEntropies(doublereal* sbar) const
{
const vector_fp& _s = entropy_R_ref();
doublereal r = GasConstant;
doublereal xx;
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(SmallNumber, moleFraction(k));
sbar[k] = r * (_s[k] - log(xx));
double xx = std::max(SmallNumber, moleFraction(k));
sbar[k] = GasConstant * (_s[k] - log(xx));
}
}
@ -234,8 +226,7 @@ void LatticePhase::getEntropy_R(doublereal* sr) const
void LatticePhase::getGibbs_RT(doublereal* grt) const
{
const vector_fp& gibbsrt = gibbs_RT_ref();
doublereal RT = _RT();
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];
}
@ -295,20 +286,18 @@ const vector_fp& LatticePhase::cp_R_ref() const
void LatticePhase::initThermo()
{
m_Pref = refPressure();
size_t leng = m_kk;
m_h0_RT.resize(leng);
m_g0_RT.resize(leng);
m_cp0_R.resize(leng);
m_s0_R.resize(leng);
m_speciesMolarVolume.resize(leng, 0.0);
m_h0_RT.resize(m_kk);
m_g0_RT.resize(m_kk);
m_cp0_R.resize(m_kk);
m_s0_R.resize(m_kk);
m_speciesMolarVolume.resize(m_kk, 0.0);
ThermoPhase::initThermo();
}
void LatticePhase::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
std::string idattrib = phaseNode.id();
if (!id_.empty() && id_ != idattrib) {
if (!id_.empty() && id_ != phaseNode.id()) {
throw CanteraError("LatticePhase::initThermoXML",
"ids don't match");
}
@ -320,11 +309,10 @@ void LatticePhase::initThermoXML(XML_Node& phaseNode, const std::string& id_)
*/
if (phaseNode.hasChild("thermo")) {
XML_Node& thNode = phaseNode.child("thermo");
std::string mStringa = thNode.attrib("model");
std::string mString = lowercase(mStringa);
if (mString != "lattice") {
std::string mString = thNode.attrib("model");
if (lowercase(mString) != "lattice") {
throw CanteraError(subname.c_str(),
"Unknown thermo model: " + mStringa);
"Unknown thermo model: " + mString);
}
} else {
throw CanteraError(subname.c_str(),
@ -335,11 +323,10 @@ void LatticePhase::initThermoXML(XML_Node& phaseNode, const std::string& id_)
*/
XML_Node& speciesList = phaseNode.child("speciesArray");
XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"], &phaseNode.root());
const std::vector<std::string> &sss = speciesNames();
for (size_t k = 0; k < m_kk; k++) {
m_speciesMolarVolume[k] = m_site_density;
XML_Node* s = speciesDB->findByAttr("name", sss[k]);
XML_Node* s = speciesDB->findByAttr("name", speciesName(k));
if (!s) {
throw CanteraError(" LatticePhase::initThermoXML", "database problems");
}
@ -379,8 +366,7 @@ void LatticePhase::setParameters(int n, doublereal* const c)
void LatticePhase::getParameters(int& n, doublereal* const c) const
{
double d = molarDensity();
c[0] = d;
c[0] = molarDensity();
n = 1;
}

68
src/thermo/LatticeSolidPhase.cpp
View file

@ -200,9 +200,9 @@ doublereal LatticeSolidPhase::calcDensity()
void LatticeSolidPhase::setMoleFractions(const doublereal* const x)
{
size_t nsp, strt = 0;
size_t strt = 0;
for (size_t n = 0; n < m_nlattice; n++) {
nsp = m_lattice[n]->nSpecies();
size_t nsp = m_lattice[n]->nSpecies();
m_lattice[n]->setMoleFractions(x + strt);
strt += nsp;
}
@ -215,13 +215,12 @@ void LatticeSolidPhase::setMoleFractions(const doublereal* const x)
void LatticeSolidPhase::getMoleFractions(doublereal* const x) const
{
size_t nsp, strt = 0;
size_t strt = 0;
// the ifdef block should be the way we calculate this.!!!!!
Phase::getMoleFractions(x);
doublereal sum;
for (size_t n = 0; n < m_nlattice; n++) {
nsp = m_lattice[n]->nSpecies();
sum = 0.0;
size_t nsp = m_lattice[n]->nSpecies();
double sum = 0.0;
for (size_t k = 0; k < nsp; k++) {
sum += (x + strt)[k];
}
@ -332,36 +331,27 @@ void LatticeSolidPhase::installSlavePhases(Cantera::XML_Node* phaseNode)
size_t kstart = 0;
m_speciesData.clear();
XML_Node& eosdata = phaseNode->child("thermo");
XML_Node& la = eosdata.child("LatticeArray");
XML_Node& la = phaseNode->child("thermo").child("LatticeArray");
std::vector<XML_Node*> lattices = la.getChildren("phase");
for (size_t n = 0; n < m_nlattice; n++) {
LatticePhase* lp = m_lattice[n];
size_t nsp = lp->nSpecies();
vector<doublereal> constArr(lp->nElements());
const vector_fp& aws = lp->atomicWeights();
for (size_t es = 0; es < lp->nElements(); es++) {
string esName = lp->elementName(es);
double wt = aws[es];
int an = lp->atomicNumber(es);
int e298 = lp->entropyElement298(es); //! @todo Why is this an int instead of a double?
int et = lp->elementType(es);
addElement(esName, wt, an, e298, et);
addElement(lp->elementName(es), aws[es], lp->atomicNumber(es),
lp->entropyElement298(es), lp->elementType(es));
}
const std::vector<const XML_Node*> & spNode = lp->speciesData();
kstart = kk;
for (size_t k = 0; k < nsp; k++) {
std::string sname = lp->speciesName(k);
for (size_t k = 0; k < lp->nSpecies(); k++) {
std::map<std::string, double> comp;
lp->getAtoms(k, DATA_PTR(constArr));
size_t nel = nElements();
vector_fp ecomp(nel, 0.0);
vector_fp ecomp(nElements(), 0.0);
for (size_t m = 0; m < lp->nElements(); m++) {
if (constArr[m] != 0.0) {
std::string oldEname = lp->elementName(m);
size_t newIndex = elementIndex(oldEname);
size_t newIndex = elementIndex(lp->elementName(m));
if (newIndex == npos) {
throw CanteraError("LatticeSolidPhase::installSlavePhases",
"element not found");
@ -369,9 +359,8 @@ void LatticeSolidPhase::installSlavePhases(Cantera::XML_Node* phaseNode)
ecomp[newIndex] = constArr[m];
}
}
double chrg = lp->charge(k);
double sz = lp->size(k);
addUniqueSpecies(sname, &ecomp[0], chrg, sz);
addUniqueSpecies(lp->speciesName(k), &ecomp[0], lp->charge(k),
lp->size(k));
SpeciesThermoInterpType* stit = newSpeciesThermoInterpType(*spNode[k]);
stit->setIndex(kk);
stit->validate(spNode[k]->attrib("name"));
@ -383,17 +372,14 @@ void LatticeSolidPhase::installSlavePhases(Cantera::XML_Node* phaseNode)
* Add in the lattice stoichiometry constraint
*/
if (n > 0) {
string econ = "LC_";
econ += int2str(n);
econ += "_" + id();
string econ = "LC_" + int2str(n) + "_" + id();
size_t m = addElement(econ, 0.0, 0, 0.0, CT_ELEM_TYPE_LATTICERATIO);
size_t mm = nElements();
LatticePhase* lp0 = m_lattice[0];
size_t nsp0 = lp0->nSpecies();
size_t nsp0 = m_lattice[0]->nSpecies();
for (size_t k = 0; k < nsp0; k++) {
m_speciesComp[k * mm + m] = -theta_[0];
}
for (size_t k = 0; k < nsp; k++) {
for (size_t k = 0; k < lp->nSpecies(); k++) {
size_t ks = kstart + k;
m_speciesComp[ks * mm + m] = theta_[n];
}
@ -404,9 +390,9 @@ void LatticeSolidPhase::installSlavePhases(Cantera::XML_Node* phaseNode)
void LatticeSolidPhase::initThermo()
{
initLengths();
size_t nsp, loc = 0;
size_t loc = 0;
for (size_t n = 0; n < m_nlattice; n++) {
nsp = m_lattice[n]->nSpecies();
size_t nsp = m_lattice[n]->nSpecies();
lkstart_[n] = loc;
for (size_t k = 0; k < nsp; k++) {
m_x[loc] =m_lattice[n]->moleFraction(k) / (double) m_nlattice;
@ -446,10 +432,9 @@ void LatticeSolidPhase::setLatticeMoleFractionsByName(int nn, const std::string&
{
m_lattice[nn]->setMoleFractionsByName(x);
size_t loc = 0;
doublereal ndens;
for (size_t n = 0; n < m_nlattice; n++) {
size_t nsp = m_lattice[n]->nSpecies();
ndens = m_lattice[n]->molarDensity();
double ndens = m_lattice[n]->molarDensity();
for (size_t k = 0; k < nsp; k++) {
m_x[loc] = ndens * m_lattice[n]->moleFraction(k);
loc++;
@ -463,21 +448,18 @@ void LatticeSolidPhase::setParametersFromXML(const XML_Node& eosdata)
eosdata._require("model","LatticeSolid");
XML_Node& la = eosdata.child("LatticeArray");
std::vector<XML_Node*> lattices = la.getChildren("phase");
size_t nl = lattices.size();
m_nlattice = nl;
for (size_t n = 0; n < nl; n++) {
XML_Node& i = *lattices[n];
m_lattice.push_back((LatticePhase*)newPhase(i));
m_nlattice = lattices.size();
for (size_t n = 0; n < m_nlattice; n++) {
m_lattice.push_back((LatticePhase*)newPhase(*lattices[n]));
}
std::vector<string> pnam;
std::vector<string> pval;
XML_Node& ls = eosdata.child("LatticeStoichiometry");
int np = ctml::getPairs(ls, pnam, pval);
theta_.resize(nl);
int np = ctml::getPairs(eosdata.child("LatticeStoichiometry"), pnam, pval);
theta_.resize(m_nlattice);
for (int i = 0; i < np; i++) {
double val = fpValueCheck(pval[i]);
bool found = false;
for (size_t j = 0; j < nl; j++) {
for (size_t j = 0; j < m_nlattice; j++) {
ThermoPhase& tp = *(m_lattice[j]);
string idj = tp.id();
if (idj == pnam[i]) {

202
src/thermo/MargulesVPSSTP.cpp
View file

@ -174,7 +174,6 @@ void MargulesVPSSTP::getElectrochemPotentials(doublereal* mu) const
void MargulesVPSSTP::getChemPotentials(doublereal* mu) const
{
doublereal xx;
/*
* First get the standard chemical potentials in
* molar form.
@ -188,7 +187,7 @@ void MargulesVPSSTP::getChemPotentials(doublereal* mu) const
s_update_lnActCoeff();
doublereal RT = GasConstant * temperature();
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(moleFractions_[k], SmallNumber);
double xx = std::max(moleFractions_[k], SmallNumber);
mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]);
}
}
@ -287,7 +286,6 @@ void MargulesVPSSTP::getPartialMolarCp(doublereal* cpbar) const
void MargulesVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
{
double xx;
/*
* Get the nondimensional standard state entropies
*/
@ -301,7 +299,7 @@ void MargulesVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
s_update_dlnActCoeff_dT();
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(moleFractions_[k], SmallNumber);
double xx = std::max(moleFractions_[k], SmallNumber);
sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
}
/*
@ -314,9 +312,6 @@ void MargulesVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
void MargulesVPSSTP::getPartialMolarVolumes(doublereal* vbar) const
{
size_t iA, iB;
double XA, XB, g0 , g1;
double T = temperature();
/*
@ -325,12 +320,12 @@ void MargulesVPSSTP::getPartialMolarVolumes(doublereal* vbar) const
getStandardVolumes(vbar);
for (size_t i = 0; i < numBinaryInteractions_; i++) {
iA = m_pSpecies_A_ij[i];
iB = m_pSpecies_B_ij[i];
XA = moleFractions_[iA];
XB = moleFractions_[iB];
g0 = (m_VHE_b_ij[i] - T * m_VSE_b_ij[i]);
g1 = (m_VHE_c_ij[i] - T * m_VSE_c_ij[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_VHE_b_ij[i] - T * m_VSE_b_ij[i]);
double g1 = (m_VHE_c_ij[i] - T * m_VSE_c_ij[i]);
const doublereal temp1 = g0 + g1 * XB;
const doublereal all = -1.0*XA*XB*temp1 - XA*XB*XB*g1;
@ -356,12 +351,10 @@ void MargulesVPSSTP::initLengths()
void MargulesVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
string stemp;
string subname = "MargulesVPSSTP::initThermoXML";
if ((int) id_.size() > 0) {
string idp = phaseNode.id();
if (idp != id_) {
throw CanteraError(subname, "phasenode and Id are incompatible");
throw CanteraError("MargulesVPSSTP::initThermoXML", "phasenode and Id are incompatible");
}
}
@ -369,7 +362,7 @@ void MargulesVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
* Find the Thermo XML node
*/
if (!phaseNode.hasChild("thermo")) {
throw CanteraError(subname,
throw CanteraError("MargulesVPSSTP::initThermoXML",
"no thermo XML node");
}
XML_Node& thermoNode = phaseNode.child("thermo");
@ -377,10 +370,9 @@ void MargulesVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
/*
* Make sure that the thermo model is Margules
*/
stemp = thermoNode.attrib("model");
string formString = lowercase(stemp);
string formString = lowercase(thermoNode.attrib("model"));
if (formString != "margules") {
throw CanteraError(subname,
throw CanteraError("MargulesVPSSTP::initThermoXML",
"model name isn't Margules: " + formString);
}
@ -389,28 +381,22 @@ void MargulesVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
* Go get all of the coefficients and factors in the
* activityCoefficients XML block
*/
XML_Node* acNodePtr = 0;
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& acNode = thermoNode.child("activityCoefficients");
acNodePtr = &acNode;
string mStringa = acNode.attrib("model");
string mString = lowercase(mStringa);
if (mString != "margules") {
throw CanteraError(subname.c_str(),
if (lowercase(mStringa) != "margules") {
throw CanteraError("MargulesVPSSTP::initThermoXML",
"Unknown activity coefficient model: " + mStringa);
}
for (size_t i = 0; i < acNodePtr->nChildren(); i++) {
XML_Node& xmlACChild = acNodePtr->child(i);
stemp = xmlACChild.name();
string nodeName = lowercase(stemp);
for (size_t i = 0; i < acNode.nChildren(); i++) {
XML_Node& xmlACChild = acNode.child(i);
/*
* Process a binary salt field, or any of the other XML fields
* that make up the Pitzer Database. Entries will be ignored
* if any of the species in the entry isn't in the solution.
*/
if (nodeName == "binaryneutralspeciesparameters") {
if (lowercase(xmlACChild.name()) == "binaryneutralspeciesparameters") {
readXMLBinarySpecies(xmlACChild);
}
}
}
@ -419,31 +405,24 @@ void MargulesVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
* Go down the chain
*/
GibbsExcessVPSSTP::initThermoXML(phaseNode, id_);
}
void MargulesVPSSTP::s_update_lnActCoeff() const
{
size_t iA, iB, iK;
double XA, XB, g0 , g1;
double T = temperature();
double invRT = 1.0 / (GasConstant*T);
lnActCoeff_Scaled_.resize(m_kk);
for (iK = 0; iK < m_kk; iK++) {
lnActCoeff_Scaled_[iK] = 0.0;
}
lnActCoeff_Scaled_.assign(m_kk, 0.0);
for (size_t i = 0; i < numBinaryInteractions_; i++) {
iA = m_pSpecies_A_ij[i];
iB = m_pSpecies_B_ij[i];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) * invRT;
g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) * invRT;
XA = moleFractions_[iA];
XB = moleFractions_[iB];
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 XA = moleFractions_[iA];
double XB = moleFractions_[iB];
const doublereal XAXB = XA * XB;
const doublereal g0g1XB = (g0 + g1 * XB);
const doublereal all = -1.0 * XAXB * g0g1XB - XAXB * XB * g1;
for (iK = 0; iK < m_kk; iK++) {
for (size_t iK = 0; iK < m_kk; iK++) {
lnActCoeff_Scaled_[iK] += all;
}
lnActCoeff_Scaled_[iA] += XB * g0g1XB;
@ -453,29 +432,23 @@ void MargulesVPSSTP::s_update_lnActCoeff() const
void MargulesVPSSTP::s_update_dlnActCoeff_dT() const
{
size_t iA, iB, iK;
doublereal XA, XB, g0, g1;
doublereal invT = 1.0 / temperature();
doublereal invRTT = 1.0 / (GasConstant)*invT*invT;
dlnActCoeffdT_Scaled_.resize(m_kk);
d2lnActCoeffdT2_Scaled_.resize(m_kk);
for (iK = 0; iK < m_kk; iK++) {
dlnActCoeffdT_Scaled_[iK] = 0.0;
d2lnActCoeffdT2_Scaled_[iK] = 0.0;
}
dlnActCoeffdT_Scaled_.assign(m_kk, 0.0);
d2lnActCoeffdT2_Scaled_.assign(m_kk, 0.0);
for (size_t i = 0; i < numBinaryInteractions_; i++) {
iA = m_pSpecies_A_ij[i];
iB = m_pSpecies_B_ij[i];
XA = moleFractions_[iA];
XB = moleFractions_[iB];
g0 = -m_HE_b_ij[i] * invRTT;
g1 = -m_HE_c_ij[i] * invRTT;
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] * invRTT;
double g1 = -m_HE_c_ij[i] * invRTT;
const doublereal XAXB = XA * XB;
const doublereal g0g1XB = (g0 + g1 * XB);
const doublereal all = -1.0 * XAXB * g0g1XB - XAXB * XB * g1;
const doublereal mult = 2.0 * invT;
const doublereal dT2all = mult * all;
for (iK = 0; iK < m_kk; iK++) {
for (size_t iK = 0; iK < m_kk; iK++) {
dlnActCoeffdT_Scaled_[iK] += all;
d2lnActCoeffdT2_Scaled_[iK] -= dT2all;
}
@ -505,29 +478,26 @@ void MargulesVPSSTP::getd2lnActCoeffdT2(doublereal* d2lnActCoeffdT2) const
void MargulesVPSSTP::getdlnActCoeffds(const doublereal dTds, const doublereal* const dXds,
doublereal* dlnActCoeffds) const
{
size_t iA, iB, iK;
double XA, XB, g0 , g1, dXA, dXB;
double T = temperature();
double RT = GasConstant*T;
s_update_dlnActCoeff_dT();
for (iK = 0; iK < m_kk; iK++) {
for (size_t iK = 0; iK < m_kk; iK++) {
dlnActCoeffds[iK] = 0.0;
}
for (size_t i = 0; i < numBinaryInteractions_; i++) {
iA = m_pSpecies_A_ij[i];
iB = m_pSpecies_B_ij[i];
XA = moleFractions_[iA];
XB = moleFractions_[iB];
dXA = dXds[iA];
dXB = dXds[iB];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
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 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;
const doublereal g02g1XB = g0 + 2*g1*XB;
const doublereal g2XAdXB = 2*g1*XA*dXB;
const doublereal all = (-XB * dXA - XA *dXB) * g02g1XB - XB *g2XAdXB;
for (iK = 0; iK < m_kk; iK++) {
for (size_t iK = 0; iK < m_kk; iK++) {
dlnActCoeffds[iK] += all + dlnActCoeffdT_Scaled_[iK]*dTds;
}
dlnActCoeffds[iA] += dXB * g02g1XB;
@ -537,24 +507,19 @@ void MargulesVPSSTP::getdlnActCoeffds(const doublereal dTds, const doublereal*
void MargulesVPSSTP::s_update_dlnActCoeff_dlnN_diag() const
{
size_t iA, iB, iK, delAK, delBK;
double XA, XB, XK, g0 , g1;
double T = temperature();
double RT = GasConstant*T;
dlnActCoeffdlnN_diag_.assign(m_kk, 0.0);
for (iK = 0; iK < m_kk; iK++) {
XK = moleFractions_[iK];
for (size_t iK = 0; iK < m_kk; iK++) {
double XK = moleFractions_[iK];
for (size_t i = 0; i < numBinaryInteractions_; i++) {
iA = m_pSpecies_A_ij[i];
iB = m_pSpecies_B_ij[i];
delAK = 0;
delBK = 0;
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
size_t delAK = 0;
size_t delBK = 0;
if (iA==iK) {
delAK = 1;
@ -562,11 +527,11 @@ void MargulesVPSSTP::s_update_dlnActCoeff_dlnN_diag() const
delBK = 1;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
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)));
}
@ -576,14 +541,8 @@ void MargulesVPSSTP::s_update_dlnActCoeff_dlnN_diag() const
void MargulesVPSSTP::s_update_dlnActCoeff_dlnN() const
{
size_t iA, iB;
doublereal delAK, delBK;
double XA, XB, g0, g1,XM;
double T = temperature();
double RT = GasConstant*T;
doublereal delAM, delBM;
dlnActCoeffdlnN_.zero();
/*
@ -591,16 +550,16 @@ void MargulesVPSSTP::s_update_dlnActCoeff_dlnN() const
*/
for (size_t iK = 0; iK < m_kk; iK++) {
for (size_t iM = 0; iM < m_kk; iM++) {
XM = moleFractions_[iM];
double XM = moleFractions_[iM];
for (size_t i = 0; i < numBinaryInteractions_; i++) {
iA = m_pSpecies_A_ij[i];
iB = m_pSpecies_B_ij[i];
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
delAK = 0.0;
delBK = 0.0;
delAM = 0.0;
delBM = 0.0;
double delAK = 0.0;
double delBK = 0.0;
double delAM = 0.0;
double delBM = 0.0;
if (iA==iK) {
delAK = 1.0;
} else if (iB==iK) {
@ -612,11 +571,11 @@ void MargulesVPSSTP::s_update_dlnActCoeff_dlnN() const
delBM = 1.0;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
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);
@ -701,8 +660,6 @@ void MargulesVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies",
"Incorrect name for processing this routine: " + xname);
}
string stemp;
size_t nParamsFound;
vector_fp vParams;
string aName = xmLBinarySpecies.attrib("speciesA");
if (aName == "") {
@ -726,18 +683,16 @@ void MargulesVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
// @TODO Figure out what the original reason is for putting an error condition for charged species
// Seems OK to me.
//
double chargeA = charge(aSpecies);
if (chargeA != 0.0) {
throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", "speciesA has a charge: " + fp2str(chargeA));
if (charge(aSpecies) != 0.0) {
throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", "speciesA has a charge: " + fp2str(charge(aSpecies)));
}
size_t bSpecies = speciesIndex(bName);
if (bSpecies == npos) {
return;
}
string bspName = speciesName(bSpecies);
double chargeB = charge(bSpecies);
if (chargeB != 0.0) {
throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", "speciesB has a charge: " + fp2str(chargeB));
if (charge(bSpecies) != 0.0) {
throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", "speciesB has a charge: " + fp2str(charge(bSpecies)));
}
resizeNumInteractions(numBinaryInteractions_ + 1);
@ -747,8 +702,7 @@ void MargulesVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
for (size_t iChild = 0; iChild < xmLBinarySpecies.nChildren(); iChild++) {
XML_Node& xmlChild = xmLBinarySpecies.child(iChild);
stemp = xmlChild.name();
string nodeName = lowercase(stemp);
string nodeName = lowercase(xmlChild.name());
/*
* Process the binary species interaction parameters.
* They are in subblocks labeled:
@ -764,9 +718,7 @@ void MargulesVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
* Get the string containing all of the values
*/
ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessEnthalpy");
nParamsFound = vParams.size();
if (nParamsFound != 2) {
if (vParams.size() != 2) {
throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessEnthalpy for " + aspName
+ "::" + bspName,
"wrong number of params found. Need 2");
@ -780,9 +732,7 @@ void MargulesVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
* Get the string containing all of the values
*/
ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessEntropy");
nParamsFound = vParams.size();
if (nParamsFound != 2) {
if (vParams.size() != 2) {
throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessEntropy for " + aspName
+ "::" + bspName,
"wrong number of params found. Need 2");
@ -796,9 +746,7 @@ void MargulesVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
* Get the string containing all of the values
*/
ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Enthalpy");
nParamsFound = vParams.size();
if (nParamsFound != 2) {
if (vParams.size() != 2) {
throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessVolume_Enthalpy for " + aspName
+ "::" + bspName,
"wrong number of params found. Need 2");
@ -812,9 +760,7 @@ void MargulesVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
* Get the string containing all of the values
*/
ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Entropy");
nParamsFound = vParams.size();
if (nParamsFound != 2) {
if (vParams.size() != 2) {
throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessVolume_Entropy for " + aspName
+ "::" + bspName,
"wrong number of params found. Need 2");

53
src/thermo/MaskellSolidSolnPhase.cpp
View file

@ -125,8 +125,7 @@ void MaskellSolidSolnPhase::calcDensity()
for (size_t i = 0; i < m_kk; i++) {
vtotal += vbar[i] * moleFracs[i];
}
doublereal dd = meanMolecularWeight() / vtotal;
Phase::setDensity(dd);
Phase::setDensity(meanMolecularWeight() / vtotal);
}
void MaskellSolidSolnPhase::setPressure(doublereal p)
@ -154,12 +153,10 @@ void MaskellSolidSolnPhase::getActivityCoefficients(doublereal* ac) const
const doublereal r = moleFraction(product_species_index);
const doublereal pval = p(r);
const doublereal fmval = fm(r);
const doublereal rfm = r * fmval;
const doublereal RT = GasConstant * temperature();
const doublereal rfm = r * fm(r);
const doublereal A = (std::pow(1 - rfm, pval) * std::pow(rfm, pval) * std::pow(r - rfm, 1 - pval)) /
(std::pow(1 - r - rfm, 1 + pval) * (1 - r));
const doublereal B = pval * h_mixing / RT;
const doublereal B = pval * h_mixing / (GasConstant * temperature());
cached.value[product_species_index] = A * std::exp(B);
cached.value[reactant_species_index] = 1 / (A * r * (1-r) ) * std::exp(-B);
}
@ -171,8 +168,7 @@ void MaskellSolidSolnPhase::getChemPotentials(doublereal* mu) const
_updateThermo();
const doublereal r = moleFraction(product_species_index);
const doublereal pval = p(r);
const doublereal fmval = fm(r);
const doublereal rfm = r * fmval;
const doublereal rfm = r * fm(r);
const doublereal RT = GasConstant * temperature();
const doublereal DgbarDr = pval * h_mixing +
GasConstant * temperature() *
@ -237,13 +233,9 @@ void MaskellSolidSolnPhase::getStandardChemPotentials(doublereal* mu) const
*********************************************************************/
void MaskellSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
std::string subname = "MaskellSolidSolnPhase::initThermoXML";
if (id_.size() > 0) {
std::string idp = phaseNode.id();
if (idp != id_) {
throw CanteraError(subname.c_str(),
"phasenode and Id are incompatible");
}
if (id_.size() > 0 && phaseNode.id() != id_) {
throw CanteraError("MaskellSolidSolnPhase::initThermoXML",
"phasenode and Id are incompatible");
}
/*
@ -252,32 +244,28 @@ void MaskellSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string
*/
if (phaseNode.hasChild("thermo")) {
XML_Node& thNode = phaseNode.child("thermo");
std::string mStringa = thNode.attrib("model");
std::string mString = lowercase(mStringa);
if (mString != "maskellsolidsolnphase") {
throw CanteraError(subname.c_str(),
"Unknown thermo model: " + mStringa);
std::string mString = thNode.attrib("model");
if (lowercase(mString) != "maskellsolidsolnphase") {
throw CanteraError("MaskellSolidSolnPhase::initThermoXML",
"Unknown thermo model: " + mString);
}
/*
* Parse the enthalpy of mixing constant
*/
if (thNode.hasChild("h_mix")) {
XML_Node& scNode = thNode.child("h_mix");
set_h_mix(fpValue(scNode.value()));
set_h_mix(fpValue(thNode.child("h_mix").value()));
} else {
throw CanteraError(subname.c_str(),
throw CanteraError("MaskellSolidSolnPhase::initThermoXML",
"Mixing enthalpy parameter not specified.");
}
if (thNode.hasChild("product_species")) {
XML_Node& scNode = thNode.child("product_species");
std::string product_species_name = scNode.value();
std::string product_species_name = thNode.child("product_species").value();
product_species_index = speciesIndex(product_species_name);
if( product_species_index == static_cast<int>(npos) )
{
throw CanteraError(subname.c_str(),
throw CanteraError("MaskellSolidSolnPhase::initThermoXML",
"Species " + product_species_name + " not found.");
}
if( product_species_index == 0 )
@ -289,15 +277,15 @@ void MaskellSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string
}
} else {
throw CanteraError(subname.c_str(),
throw CanteraError("MaskellSolidSolnPhase::initThermoXML",
"Unspecified thermo model");
}
// Confirm that the phase only contains 2 species
if( m_kk != 2 )
{
throw CanteraError( subname.c_str(), "MaskellSolidSolution model requires exactly 2 species.");
if (m_kk != 2) {
throw CanteraError("MaskellSolidSolnPhase::initThermoXML",
"MaskellSolidSolution model requires exactly 2 species.");
}
/*
@ -333,8 +321,7 @@ doublereal MaskellSolidSolnPhase::s() const
doublereal MaskellSolidSolnPhase::fm(const doublereal r) const
{
const doublereal sval = s();
return (1 - std::sqrt(1 - 4*r*(1-r)/sval)) / (2*r);
return (1 - std::sqrt(1 - 4*r*(1-r)/s())) / (2*r);
}
doublereal MaskellSolidSolnPhase::p(const doublereal r) const

20
src/thermo/MetalSHEelectrons.cpp
View file

@ -48,9 +48,7 @@ MetalSHEelectrons::MetalSHEelectrons(const std::string& infile, std::string id_)
"Couldn't find phase name in file:" + id_);
}
// Check the model name to ensure we have compatibility
const XML_Node& th = xphase->child("thermo");
std::string model = th["model"];
if (model != "MetalSHEelectrons") {
if (xphase->child("thermo")["model"] != "MetalSHEelectrons") {
throw CanteraError("MetalSHEelectrons::MetalSHEelectrons",
"thermo model attribute must be MetalSHEelectrons");
}
@ -61,15 +59,12 @@ MetalSHEelectrons::MetalSHEelectrons(XML_Node& xmlphase, const std::string& id_)
xdef_(0)
{
if (id_ != "") {
std::string idxml = xmlphase["id"];
if (id_ != idxml) {
if (id_ != xmlphase["id"]) {
throw CanteraError("MetalSHEelectrons::MetalSHEelectrons",
"id's don't match");
}
}
const XML_Node& th = xmlphase.child("thermo");
std::string model = th["model"];
if (model != "MetalSHEelectrons") {
if (xmlphase.child("thermo")["model"] != "MetalSHEelectrons") {
throw CanteraError("MetalSHEelectrons::MetalSHEelectrons",
"thermo model attribute must be MetalSHEelectrons");
}
@ -293,21 +288,18 @@ XML_Node* MetalSHEelectrons::makeDefaultXMLTree()
void MetalSHEelectrons::setParameters(int n, doublereal* const c)
{
doublereal rho = c[0];
setDensity(rho);
setDensity(c[0]);
}
void MetalSHEelectrons::getParameters(int& n, doublereal* const c) const
{
doublereal rho = density();
n = 1;
c[0] = rho;
c[0] = density();
}
void MetalSHEelectrons::setParametersFromXML(const XML_Node& eosdata)
{
std::string model = eosdata["model"];
if (model != "MetalSHEelectrons") {
if ( eosdata["model"] != "MetalSHEelectrons") {
throw CanteraError("MetalSHEelectrons::setParametersFromXML",
"thermo model attribute must be MetalSHEelectrons");
}

44
src/thermo/MineralEQ3.cpp
View file

@ -43,8 +43,7 @@ MineralEQ3::MineralEQ3(const std::string& infile, std::string id_)
"Couldn't find phase name in file:" + id_);
}
// Check the model name to ensure we have compatibility
const XML_Node& th = xphase->child("thermo");
std::string model = th["model"];
std::string model = xphase->child("thermo")["model"];
if (model != "StoichSubstance" && model != "MineralEQ3") {
throw CanteraError("MineralEQ3::MineralEQ3",
"thermo model attribute must be StoichSubstance");
@ -55,14 +54,12 @@ MineralEQ3::MineralEQ3(const std::string& infile, std::string id_)
MineralEQ3::MineralEQ3(XML_Node& xmlphase, const std::string& id_)
{
if (id_ != "") {
std::string idxml = xmlphase["id"];
if (id_ != idxml) {
if (id_ != xmlphase["id"]) {
throw CanteraError("MineralEQ3::MineralEQ3",
"id's don't match");
}
}
const XML_Node& th = xmlphase.child("thermo");
std::string model = th["model"];
std::string model = xmlphase.child("thermo")["model"];
if (model != "StoichSubstance" && model != "MineralEQ3") {
throw CanteraError("MineralEQ3::MineralEQ3",
"thermo model attribute must be StoichSubstance");
@ -217,15 +214,13 @@ void MineralEQ3::getIntEnergy_RT_ref(doublereal* urt) const
void MineralEQ3::setParameters(int n, doublereal* const c)
{
doublereal rho = c[0];
setDensity(rho);
setDensity(c[0]);
}
void MineralEQ3::getParameters(int& n, doublereal* const c) const
{
doublereal rho = density();
n = 1;
c[0] = rho;
c[0] = density();
}
void MineralEQ3::initThermoXML(XML_Node& phaseNode, const std::string& id_)
@ -238,8 +233,7 @@ void MineralEQ3::initThermoXML(XML_Node& phaseNode, const std::string& id_)
"no thermo XML node");
}
std::vector<const XML_Node*> xspecies = speciesData();
const XML_Node* xsp = xspecies[0];
const XML_Node* xsp = speciesData()[0];
XML_Node* aStandardState = 0;
if (xsp->hasChild("standardState")) {
@ -249,18 +243,15 @@ void MineralEQ3::initThermoXML(XML_Node& phaseNode, const std::string& id_)
"no standard state mode");
}
doublereal volVal = 0.0;
string smodel = aStandardState->attrib("model");
if (smodel != "constantVolume") {
if (aStandardState->attrib("model") != "constantVolume") {
throw CanteraError("MineralEQ3::initThermoXML",
"wrong standard state mode");
}
if (aStandardState->hasChild("V0_Pr_Tr")) {
XML_Node& aV = aStandardState->child("V0_Pr_Tr");
string Aunits = "";
double Afactor = toSI("cm3/gmol");
if (aV.hasAttrib("units")) {
Aunits = aV.attrib("units");
Afactor = toSI(Aunits);
Afactor = toSI(aV.attrib("units"));
}
volVal = ctml::getFloat(*aStandardState, "V0_Pr_Tr");
m_V0_pr_tr= volVal;
@ -270,12 +261,9 @@ void MineralEQ3::initThermoXML(XML_Node& phaseNode, const std::string& id_)
throw CanteraError("MineralEQ3::initThermoXML",
"wrong standard state mode");
}
doublereal rho = molecularWeight(0) / volVal;
setDensity(rho);
const XML_Node& sThermo = xsp->child("thermo");
const XML_Node& MinEQ3node = sThermo.child("MinEQ3");
setDensity(molecularWeight(0) / volVal);
const XML_Node& MinEQ3node = xsp->child("thermo").child("MinEQ3");
m_deltaG_formation_pr_tr =
ctml::getFloatDefaultUnits(MinEQ3node, "DG0_f_Pr_Tr", "cal/gmol", "actEnergy");
@ -291,8 +279,7 @@ void MineralEQ3::initThermoXML(XML_Node& phaseNode, const std::string& id_)
void MineralEQ3::setParametersFromXML(const XML_Node& eosdata)
{
std::string model = eosdata["model"];
if (model != "MineralEQ3") {
if (eosdata["model"] != "MineralEQ3") {
throw CanteraError("MineralEQ3::MineralEQ3",
"thermo model attribute must be MineralEQ3");
}
@ -318,17 +305,12 @@ void MineralEQ3::convertDGFormation()
/*
* Ok let's get the element compositions and conversion factors.
*/
doublereal na;
doublereal ge;
string ename;
doublereal totalSum = 0.0;
for (size_t m = 0; m < nElements(); m++) {
na = nAtoms(0, m);
double na = nAtoms(0, m);
if (na > 0.0) {
ename = elementName(m);
ge = LookupGe(ename);
totalSum += na * ge;
totalSum += na * LookupGe(elementName(m));
}
}
// Ok, now do the calculation. Convert to joules kmol-1

197
src/thermo/MixedSolventElectrolyte.cpp
View file

@ -177,7 +177,6 @@ void MixedSolventElectrolyte::getElectrochemPotentials(doublereal* mu) const
void MixedSolventElectrolyte::getChemPotentials(doublereal* mu) const
{
doublereal xx;
/*
* First get the standard chemical potentials in
* molar form.
@ -194,18 +193,17 @@ void MixedSolventElectrolyte::getChemPotentials(doublereal* mu) const
*/
doublereal RT = GasConstant * temperature();
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(moleFractions_[k], SmallNumber);
double xx = std::max(moleFractions_[k], SmallNumber);
mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]);
}
}
doublereal MixedSolventElectrolyte::enthalpy_mole() const
{
size_t kk = nSpecies();
double h = 0;
vector_fp hbar(kk);
vector_fp hbar(m_kk);
getPartialMolarEnthalpies(&hbar[0]);
for (size_t i = 0; i < kk; i++) {
for (size_t i = 0; i < m_kk; i++) {
h += moleFractions_[i]*hbar[i];
}
return h;
@ -213,11 +211,10 @@ doublereal MixedSolventElectrolyte::enthalpy_mole() const
doublereal MixedSolventElectrolyte::entropy_mole() const
{
size_t kk = nSpecies();
double s = 0;
vector_fp sbar(kk);
vector_fp sbar(m_kk);
getPartialMolarEntropies(&sbar[0]);
for (size_t i = 0; i < kk; i++) {
for (size_t i = 0; i < m_kk; i++) {
s += moleFractions_[i]*sbar[i];
}
return s;
@ -225,11 +222,10 @@ doublereal MixedSolventElectrolyte::entropy_mole() const
doublereal MixedSolventElectrolyte::cp_mole() const
{
size_t kk = nSpecies();
double cp = 0;
vector_fp cpbar(kk);
vector_fp cpbar(m_kk);
getPartialMolarCp(&cpbar[0]);
for (size_t i = 0; i < kk; i++) {
for (size_t i = 0; i < m_kk; i++) {
cp += moleFractions_[i]*cpbar[i];
}
return cp;
@ -293,7 +289,6 @@ void MixedSolventElectrolyte::getPartialMolarCp(doublereal* cpbar) const
void MixedSolventElectrolyte::getPartialMolarEntropies(doublereal* sbar) const
{
double xx;
/*
* Get the nondimensional standard state entropies
*/
@ -307,7 +302,7 @@ void MixedSolventElectrolyte::getPartialMolarEntropies(doublereal* sbar) const
s_update_dlnActCoeff_dT();
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(moleFractions_[k], SmallNumber);
double xx = std::max(moleFractions_[k], SmallNumber);
sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
}
/*
@ -320,8 +315,6 @@ void MixedSolventElectrolyte::getPartialMolarEntropies(doublereal* sbar) const
void MixedSolventElectrolyte::getPartialMolarVolumes(doublereal* vbar) const
{
int delAK, delBK;
double XA, XB, g0 , g1;
double T = temperature();
/*
@ -330,8 +323,8 @@ void MixedSolventElectrolyte::getPartialMolarVolumes(doublereal* vbar) const
getStandardVolumes(vbar);
for (size_t iK = 0; iK < m_kk; iK++) {
delAK = 0;
delBK = 0;
int delAK = 0;
int delBK = 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];
@ -342,11 +335,11 @@ void MixedSolventElectrolyte::getPartialMolarVolumes(doublereal* vbar) const
delBK = 1;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
g0 = (m_VHE_b_ij[i] - T * m_VSE_b_ij[i]);
g1 = (m_VHE_c_ij[i] - T * m_VSE_c_ij[i]);
double g0 = (m_VHE_b_ij[i] - T * m_VSE_b_ij[i]);
double g1 = (m_VHE_c_ij[i] - T * m_VSE_c_ij[i]);
vbar[iK] += XA*XB*(g0+g1*XB)+((delAK-XA)*XB+XA*(delBK-XB))*(g0+g1*XB)+XA*XB*(delBK-XB)*g1;
}
@ -367,14 +360,9 @@ void MixedSolventElectrolyte::initLengths()
void MixedSolventElectrolyte::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
string subname = "MixedSolventElectrolyte::initThermoXML";
string stemp;
if ((int) id_.size() > 0) {
string idp = phaseNode.id();
if (idp != id_) {
throw CanteraError(subname, "phasenode and Id are incompatible");
}
if ((int) id_.size() > 0 && phaseNode.id() != id_) {
throw CanteraError("MixedSolventElectrolyte::initThermoXML",
"phasenode and Id are incompatible");
}
/*
@ -382,40 +370,35 @@ void MixedSolventElectrolyte::initThermoXML(XML_Node& phaseNode, const std::stri
* <thermo model="MixedSolventElectrolyte" />
*/
if (!phaseNode.hasChild("thermo")) {
throw CanteraError(subname, "no thermo XML node");
throw CanteraError("MixedSolventElectrolyte::initThermoXML",
"no thermo XML node");
}
XML_Node& thermoNode = phaseNode.child("thermo");
string mStringa = thermoNode.attrib("model");
string mString = lowercase(mStringa);
if (mString != "mixedsolventelectrolyte") {
throw CanteraError(subname, "Unknown thermo model: " + mStringa);
string mString = thermoNode.attrib("model");
if (lowercase(mString) != "mixedsolventelectrolyte") {
throw CanteraError("MixedSolventElectrolyte::initThermoXML",
"Unknown thermo model: " + mString);
}
/*
* Go get all of the coefficients and factors in the
* activityCoefficients XML block
*/
XML_Node* acNodePtr = 0;
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& acNode = thermoNode.child("activityCoefficients");
acNodePtr = &acNode;
mStringa = acNode.attrib("model");
mString = lowercase(mStringa);
if (mString != "margules") {
throw CanteraError(subname.c_str(),
"Unknown activity coefficient model: " + mStringa);
mString = acNode.attrib("model");
if (lowercase(mString) != "margules") {
throw CanteraError("MixedSolventElectrolyte::initThermoXML",
"Unknown activity coefficient model: " + mString);
}
size_t n = acNodePtr->nChildren();
for (size_t i = 0; i < n; i++) {
XML_Node& xmlACChild = acNodePtr->child(i);
stemp = xmlACChild.name();
string nodeName = lowercase(stemp);
for (size_t i = 0; i < acNode.nChildren(); i++) {
XML_Node& xmlACChild = acNode.child(i);
/*
* Process a binary salt field, or any of the other XML fields
* that make up the Pitzer Database. Entries will be ignored
* if any of the species in the entry isn't in the solution.
*/
if (nodeName == "binaryneutralspeciesparameters") {
if (lowercase(xmlACChild.name()) == "binaryneutralspeciesparameters") {
readXMLBinarySpecies(xmlACChild);
}
@ -432,8 +415,6 @@ void MixedSolventElectrolyte::initThermoXML(XML_Node& phaseNode, const std::stri
void MixedSolventElectrolyte::s_update_lnActCoeff() const
{
int delAK, delBK;
double XA, XB, g0, g1;
double T = temperature();
double RT = GasConstant*T;
lnActCoeff_Scaled_.assign(m_kk, 0.0);
@ -441,17 +422,17 @@ void MixedSolventElectrolyte::s_update_lnActCoeff() const
for (size_t i = 0; i < numBinaryInteractions_; i++) {
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
delAK = 0;
delBK = 0;
int delAK = 0;
int delBK = 0;
if (iA==iK) {
delAK = 1;
} else if (iB==iK) {
delBK = 1;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
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;
lnActCoeff_Scaled_[iK] += (delAK * XB + XA * delBK - XA * XB) * (g0 + g1 * XB) + XA * XB * (delBK - XB) * g1;
}
}
@ -459,8 +440,6 @@ void MixedSolventElectrolyte::s_update_lnActCoeff() const
void MixedSolventElectrolyte::s_update_dlnActCoeff_dT() const
{
int delAK, delBK;
doublereal XA, XB, g0, g1;
doublereal T = temperature();
doublereal RTT = GasConstant*T*T;
dlnActCoeffdT_Scaled_.assign(m_kk, 0.0);
@ -469,17 +448,17 @@ void MixedSolventElectrolyte::s_update_dlnActCoeff_dT() const
for (size_t i = 0; i < numBinaryInteractions_; i++) {
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
delAK = 0;
delBK = 0;
int delAK = 0;
int delBK = 0;
if (iA==iK) {
delAK = 1;
} else if (iB==iK) {
delBK = 1;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
g0 = -m_HE_b_ij[i] / RTT;
g1 = -m_HE_c_ij[i] / RTT;
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
double g0 = -m_HE_b_ij[i] / RTT;
double g1 = -m_HE_c_ij[i] / RTT;
double temp = (delAK * XB + XA * delBK - XA * XB) * (g0 + g1 * XB) + XA * XB * (delBK - XB) * g1;
dlnActCoeffdT_Scaled_[iK] += temp;
d2lnActCoeffdT2_Scaled_[iK] -= 2.0 * temp / T;
@ -506,11 +485,9 @@ void MixedSolventElectrolyte::getd2lnActCoeffdT2(doublereal* d2lnActCoeffdT2) co
void MixedSolventElectrolyte::getdlnActCoeffds(const doublereal dTds, const doublereal* const dXds,
doublereal* dlnActCoeffds) const
{
int delAK, delBK;
double XA, XB, g0, g1, dXA, dXB;
double T = temperature();
double RT = GasConstant*T;
s_update_dlnActCoeff_dT();
s_update_dlnActCoeff_dT();
for (size_t iK = 0; iK < m_kk; iK++) {
dlnActCoeffds[iK] = 0.0;
@ -519,8 +496,8 @@ void MixedSolventElectrolyte::getdlnActCoeffds(const doublereal dTds, const dou
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
delAK = 0;
delBK = 0;
int delAK = 0;
int delBK = 0;
if (iA==iK) {
delAK = 1;
@ -528,14 +505,14 @@ void MixedSolventElectrolyte::getdlnActCoeffds(const doublereal dTds, const dou
delBK = 1;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
dXA = dXds[iA];
dXB = dXds[iB];
double dXA = dXds[iA];
double dXB = dXds[iB];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
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;
@ -545,24 +522,21 @@ void MixedSolventElectrolyte::getdlnActCoeffds(const doublereal dTds, const dou
void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN_diag() const
{
int delAK, delBK;
double XA, XB, XK, g0, g1;
double T = temperature();
double RT = GasConstant*T;
dlnActCoeffdlnN_diag_.assign(m_kk, 0);
for (size_t iK = 0; iK < m_kk; iK++) {
XK = moleFractions_[iK];
double XK = moleFractions_[iK];
for (size_t i = 0; i < numBinaryInteractions_; i++) {
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
delAK = 0;
delBK = 0;
int delAK = 0;
int delBK = 0;
if (iA==iK) {
delAK = 1;
@ -570,11 +544,11 @@ void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN_diag() const
delBK = 1;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
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)));
}
@ -584,12 +558,9 @@ void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN_diag() const
void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN() const
{
doublereal delAK, delBK;
double XA, XB, g0, g1,XM;
double T = temperature();
double RT = GasConstant*T;
doublereal delAM, delBM;
dlnActCoeffdlnN_.zero();
/*
@ -597,16 +568,16 @@ void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN() const
*/
for (size_t iK = 0; iK < m_kk; iK++) {
for (size_t iM = 0; iM < m_kk; iM++) {
XM = moleFractions_[iM];
double XM = moleFractions_[iM];
for (size_t i = 0; i < numBinaryInteractions_; i++) {
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
delAK = 0.0;
delBK = 0.0;
delAM = 0.0;
delBM = 0.0;
double delAK = 0.0;
double delBK = 0.0;
double delAM = 0.0;
double delBM = 0.0;
if (iA==iK) {
delAK = 1.0;
} else if (iB==iK) {
@ -618,11 +589,11 @@ void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN() const
delBM = 1.0;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
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);
@ -634,9 +605,7 @@ void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN() const
void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnX_diag() const
{
doublereal XA, XB, g0 , g1;
doublereal T = temperature();
dlnActCoeffdlnX_diag_.assign(m_kk, 0);
doublereal RT = GasConstant * T;
@ -645,11 +614,11 @@ void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnX_diag() const
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
XA = moleFractions_[iA];
XB = moleFractions_[iB];
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
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);
@ -712,8 +681,6 @@ void MixedSolventElectrolyte::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies",
"Incorrect name for processing this routine: " + xname);
}
string stemp;
size_t nParamsFound;
vector_fp vParams;
string iName = xmLBinarySpecies.attrib("speciesA");
if (iName == "") {
@ -749,11 +716,9 @@ void MixedSolventElectrolyte::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
m_pSpecies_A_ij[iSpot] = iSpecies;
m_pSpecies_B_ij[iSpot] = jSpecies;
size_t num = xmLBinarySpecies.nChildren();
for (size_t iChild = 0; iChild < num; iChild++) {
for (size_t iChild = 0; iChild < xmLBinarySpecies.nChildren(); iChild++) {
XML_Node& xmlChild = xmLBinarySpecies.child(iChild);
stemp = xmlChild.name();
string nodeName = lowercase(stemp);
string nodeName = lowercase(xmlChild.name());
/*
* Process the binary species interaction child elements
*/
@ -762,9 +727,7 @@ void MixedSolventElectrolyte::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
* Get the string containing all of the values
*/
ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessEnthalpy");
nParamsFound = vParams.size();
if (nParamsFound != 2) {
if (vParams.size() != 2) {
throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies::excessEnthalpy for " + ispName
+ "::" + jspName,
"wrong number of params found");
@ -778,9 +741,7 @@ void MixedSolventElectrolyte::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
* Get the string containing all of the values
*/
ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessEntropy");
nParamsFound = vParams.size();
if (nParamsFound != 2) {
if (vParams.size() != 2) {
throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies::excessEntropy for " + ispName
+ "::" + jspName,
"wrong number of params found");
@ -794,9 +755,7 @@ void MixedSolventElectrolyte::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
* Get the string containing all of the values
*/
ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Enthalpy");
nParamsFound = vParams.size();
if (nParamsFound != 2) {
if (vParams.size() != 2) {
throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies::excessVolume_Enthalpy for " + ispName
+ "::" + jspName,
"wrong number of params found");
@ -810,9 +769,7 @@ void MixedSolventElectrolyte::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
* Get the string containing all of the values
*/
ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Entropy");
nParamsFound = vParams.size();
if (nParamsFound != 2) {
if (vParams.size() != 2) {
throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies::excessVolume_Entropy for " + ispName
+ "::" + jspName,
"wrong number of params found");

71
src/thermo/MixtureFugacityTP.cpp
View file

@ -153,8 +153,7 @@ 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());
tmp *= _RT();
double tmp = log(pressure() /m_spthermo->refPressure()) * _RT();
for (size_t k = 0; k < m_kk; k++) {
g[k] += tmp;
}
@ -164,11 +163,8 @@ void MixtureFugacityTP::getIntEnergy_RT(doublereal* urt) const
{
_updateReferenceStateThermo();
copy(m_h0_RT.begin(), m_h0_RT.end(), urt);
doublereal p = pressure();
doublereal tmp = p / _RT();
doublereal v0 = _RT() / p;
for (size_t i = 0; i < m_kk; i++) {
urt[i] -= tmp * v0;
urt[i] -= 1.0;
}
}
@ -220,7 +216,6 @@ void MixtureFugacityTP::getEntropy_R_ref(doublereal* er) const
{
_updateReferenceStateThermo();
copy(m_s0_R.begin(), m_s0_R.end(), er);
return;
}
void MixtureFugacityTP::getCp_R_ref(doublereal* cpr) const
@ -232,8 +227,7 @@ void MixtureFugacityTP::getCp_R_ref(doublereal* cpr) const
void MixtureFugacityTP::getStandardVolumes_ref(doublereal* vol) const
{
_updateReferenceStateThermo();
double pp = refPressure();
doublereal v0 = _RT() / pp;
doublereal v0 = _RT() / refPressure();
for (size_t i = 0; i < m_kk; i++) {
vol[i]= v0;
}
@ -453,31 +447,23 @@ void MixtureFugacityTP::initThermoXML(XML_Node& phaseNode, const std::string& id
doublereal MixtureFugacityTP::z() const
{
doublereal p = pressure();
doublereal rho = density();
doublereal mmw = meanMolecularWeight();
doublereal molarV = mmw / rho;
doublereal rt = _RT();
return p * molarV / rt;
return pressure() * meanMolecularWeight() / (density() * _RT());
}
doublereal MixtureFugacityTP::sresid() const
{
throw CanteraError("MixtureFugacityTP::sresid()", "Base Class: not implemented");
return 0.0;
}
doublereal MixtureFugacityTP::hresid() const
{
throw CanteraError("MixtureFugacityTP::hresid()", "Base Class: not implemented");
return 0.0;
}
doublereal MixtureFugacityTP::psatEst(doublereal TKelvin) const
{
doublereal tcrit = critTemperature();
doublereal pcrit = critPressure();
doublereal tt = tcrit/TKelvin;
doublereal tt = critTemperature() / TKelvin;
if (tt < 1.0) {
return pcrit;
}
@ -488,13 +474,12 @@ doublereal MixtureFugacityTP::psatEst(doublereal TKelvin) const
doublereal MixtureFugacityTP::liquidVolEst(doublereal TKelvin, doublereal& pres) const
{
throw CanteraError("MixtureFugacityTP::liquidVolEst()", "unimplemented");
return 0.0;
}
doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa,
int phase, doublereal rhoguess)
{
double tcrit = critTemperature();
doublereal tcrit = critTemperature();
doublereal mmw = meanMolecularWeight();
if (rhoguess == -1.0) {
if (phase != -1) {
@ -526,8 +511,6 @@ doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa,
* be updated as we go.
*/
double molarVolSpinodal = vc;
doublereal pcheck = 1.0E-30 + 1.0E-8 * presPa;
doublereal presBase, dpdVBase, delMV;
bool conv = false;
/*
* We start on one side of the vc and stick with that side
@ -555,7 +538,8 @@ doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa,
* wrt the molar volume.
* This routine also returns the pressure, presBase
*/
dpdVBase = dpdVCalc(TKelvin, molarVolBase, presBase);
double presBase;
double dpdVBase = dpdVCalc(TKelvin, molarVolBase, presBase);
/*
* If dpdV is positve, then we are in the middle of the
@ -593,7 +577,7 @@ doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa,
/*
* Check for convergence
*/
if (fabs(presBase-presPa) < pcheck) {
if (fabs(presBase-presPa) < 1.0E-30 + 1.0E-8 * presPa) {
conv = true;
break;
}
@ -611,7 +595,7 @@ doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa,
* Newton's method. Then, crop it to a max value
* of 0.1 times the current volume
*/
delMV = - (presBase - presPa) / dpdV;
double delMV = - (presBase - presPa) / dpdV;
if (!gasSide || delMV < 0.0) {
if (fabs(delMV) > 0.2 * molarVolBase) {
delMV = delMV / fabs(delMV) * 0.2 * molarVolBase;
@ -675,7 +659,6 @@ void MixtureFugacityTP::updateMixingExpressions()
}
MixtureFugacityTP::spinodalFunc::spinodalFunc(MixtureFugacityTP* tp) :
ResidEval(),
m_tp(tp)
{
}
@ -683,13 +666,10 @@ MixtureFugacityTP::spinodalFunc::spinodalFunc(MixtureFugacityTP* tp) :
int MixtureFugacityTP::spinodalFunc::evalSS(const doublereal t, const doublereal* const y,
doublereal* const r)
{
int status = 0;
doublereal molarVol = y[0];
doublereal tt = m_tp->temperature();
doublereal pp;
doublereal val = m_tp->dpdVCalc(tt, molarVol, pp);
r[0] = val;
return status;
r[0] = m_tp->dpdVCalc(m_tp->temperature(), molarVol, pp);
return 0;
}
int MixtureFugacityTP::corr0(doublereal TKelvin, doublereal pres, doublereal& densLiqGuess,
@ -767,13 +747,11 @@ int MixtureFugacityTP::phaseState(bool checkState) const
doublereal MixtureFugacityTP::densSpinodalLiquid() const
{
throw CanteraError("", "unimplmented");
return 0.0;
}
doublereal MixtureFugacityTP::densSpinodalGas() const
{
throw CanteraError("", "unimplmented");
return 0.0;
}
doublereal MixtureFugacityTP::satPressure(doublereal TKelvin)
@ -808,24 +786,19 @@ doublereal MixtureFugacityTP::calculatePsat(doublereal TKelvin, doublereal& mola
// we need this because this is a non-const routine that is public
setTemperature(TKelvin);
double tcrit = critTemperature();
double RhoLiquid, RhoGas;
double RhoLiquidGood, RhoGasGood;
double densSave = density();
double tempSave = temperature();
double pres;
doublereal mw = meanMolecularWeight();
if (TKelvin < tcrit) {
if (TKelvin < critTemperature()) {
pres = psatEst(TKelvin);
// trial value = Psat from correlation
doublereal volLiquid = liquidVolEst(TKelvin, pres);
RhoLiquidGood = mw / volLiquid;
RhoGasGood = pres * mw / (GasConstant * TKelvin);
double RhoLiquidGood = mw / volLiquid;
double RhoGasGood = pres * mw / (GasConstant * TKelvin);
doublereal delGRT = 1.0E6;
doublereal liqGRT, gasGRT;
int stab;
doublereal presLast = pres;
/*
* First part of the calculation involves finding a pressure at which the
@ -922,9 +895,9 @@ doublereal MixtureFugacityTP::calculatePsat(doublereal TKelvin, doublereal& mola
}
pres = presGas;
presLast = pres;
RhoGas = RhoGasGood;
RhoLiquid = RhoLiquidGood;
double presLast = pres;
double RhoGas = RhoGasGood;
double RhoLiquid = RhoLiquidGood;
/*
@ -932,8 +905,7 @@ doublereal MixtureFugacityTP::calculatePsat(doublereal TKelvin, doublereal& mola
*/
for (int i = 0; i < 20; i++) {
stab = corr0(TKelvin, pres, RhoLiquid, RhoGas, liqGRT, gasGRT);
int stab = corr0(TKelvin, pres, RhoLiquid, RhoGas, liqGRT, gasGRT);
if (stab == 0) {
presLast = pres;
delGRT = liqGRT - gasGRT;
@ -986,8 +958,7 @@ doublereal MixtureFugacityTP::calculatePsat(doublereal TKelvin, doublereal& mola
} else {
pres = critPressure();
setState_TP(TKelvin, pres);
RhoGas = density();
molarVolGas = mw / RhoGas;
molarVolGas = mw / density();
molarVolLiquid = molarVolGas;
setState_TR(tempSave, densSave);
}
@ -997,13 +968,11 @@ doublereal MixtureFugacityTP::calculatePsat(doublereal TKelvin, doublereal& mola
doublereal MixtureFugacityTP::pressureCalc(doublereal TKelvin, doublereal molarVol) const
{
throw CanteraError("MixtureFugacityTP::pressureCalc", "unimplemented");
return 0.0;
}
doublereal MixtureFugacityTP::dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal& presCalc) const
{
throw CanteraError("MixtureFugacityTP::dpdVCalc", "unimplemented");
return 0.0;
}
void MixtureFugacityTP::_updateReferenceStateThermo() const

41
src/thermo/MolalityVPSSTP.cpp
View file

@ -180,15 +180,13 @@ void MolalityVPSSTP::setMolalitiesByName(const compositionMap& mMap)
* neutrals so that the existing mole fractions are
* preserved.
*/
size_t kk = nSpecies();
/*
* Get a vector of mole fractions
*/
vector_fp mf(kk, 0.0);
vector_fp mf(m_kk, 0.0);
getMoleFractions(DATA_PTR(mf));
double xmolS = mf[m_indexSolvent];
double xmolSmin = std::max(xmolS, m_xmolSolventMIN);
for (size_t k = 0; k < kk; k++) {
double xmolSmin = std::max(mf[m_indexSolvent], m_xmolSolventMIN);
for (size_t k = 0; k < m_kk; k++) {
double mol_k = getValue(mMap, speciesName(k), 0.0);
if (mol_k > 0) {
mf[k] = mol_k * m_Mnaught * xmolSmin;
@ -202,7 +200,7 @@ void MolalityVPSSTP::setMolalitiesByName(const compositionMap& mMap)
size_t largeNeg = npos;
double cNeg = 0.0;
double sum = 0.0;
for (size_t k = 0; k < kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
double ch = charge(k);
if (mf[k] > 0.0) {
if (ch > 0.0) {
@ -239,11 +237,11 @@ void MolalityVPSSTP::setMolalitiesByName(const compositionMap& mMap)
}
sum = 0.0;
for (size_t k = 0; k < kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
sum += mf[k];
}
sum = 1.0/sum;
for (size_t k = 0; k < kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
mf[k] *= sum;
}
setMoleFractions(DATA_PTR(mf));
@ -316,9 +314,8 @@ doublereal MolalityVPSSTP::osmoticCoefficient() const
sum += std::max(m_molalities[k], 0.0);
}
double oc = 1.0;
double lac = log(act[m_indexSolvent]);
if (sum > 1.0E-200) {
oc = - lac / (m_Mnaught * sum);
oc = - log(act[m_indexSolvent]) / (m_Mnaught * sum);
}
return oc;
}
@ -426,9 +423,8 @@ size_t MolalityVPSSTP::findCLMIndex() const
size_t eCl = npos;
size_t eE = npos;
size_t ne = nElements();
string sn;
for (size_t e = 0; e < ne; e++) {
sn = elementName(e);
string sn = elementName(e);
if (sn == "Cl" || sn == "CL") {
eCl = e;
break;
@ -439,7 +435,7 @@ size_t MolalityVPSSTP::findCLMIndex() const
return npos;
}
for (size_t e = 0; e < ne; e++) {
sn = elementName(e);
string sn = elementName(e);
if (sn == "E" || sn == "e") {
eE = e;
break;
@ -466,7 +462,7 @@ size_t MolalityVPSSTP::findCLMIndex() const
}
}
}
sn = speciesName(k);
string sn = speciesName(k);
if (sn != "Cl-" && sn != "CL-") {
continue;
}
@ -522,13 +518,12 @@ std::string MolalityVPSSTP::report(bool show_thermo, doublereal threshold) const
sprintf(p, " potential %12.6g V\n", phi);
s += p;
size_t kk = nSpecies();
vector_fp x(kk);
vector_fp molal(kk);
vector_fp mu(kk);
vector_fp muss(kk);
vector_fp acMolal(kk);
vector_fp actMolal(kk);
vector_fp x(m_kk);
vector_fp molal(m_kk);
vector_fp mu(m_kk);
vector_fp muss(m_kk);
vector_fp acMolal(m_kk);
vector_fp actMolal(m_kk);
getMoleFractions(&x[0]);
getMolalities(&molal[0]);
getChemPotentials(&mu[0]);
@ -590,7 +585,7 @@ std::string MolalityVPSSTP::report(bool show_thermo, doublereal threshold) const
sprintf(p, " ------------- "
" ------------ ------------ ------------ ------------\n");
s += p;
for (size_t k = 0; k < kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
if (x[k] > threshold) {
if (x[k] > SmallNumber) {
sprintf(p, "%18s %12.6g %12.6g %12.6g %12.6g %12.6g\n",
@ -612,7 +607,7 @@ std::string MolalityVPSSTP::report(bool show_thermo, doublereal threshold) const
sprintf(p, " -------------"
" ------------\n");
s += p;
for (size_t k = 0; k < kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
if (x[k] > threshold) {
sprintf(p, "%18s %12.6g %12.6g\n",
speciesName(k).c_str(), x[k], molal[k]);

108
src/thermo/MolarityIonicVPSSTP.cpp
View file

@ -125,7 +125,6 @@ void MolarityIonicVPSSTP::getLnActivityCoefficients(doublereal* lnac) const
void MolarityIonicVPSSTP::getChemPotentials(doublereal* mu) const
{
doublereal xx;
/*
* First get the standard chemical potentials in
* molar form.
@ -142,7 +141,7 @@ void MolarityIonicVPSSTP::getChemPotentials(doublereal* mu) const
*/
doublereal RT = GasConstant * temperature();
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(moleFractions_[k], SmallNumber);
double xx = std::max(moleFractions_[k], SmallNumber);
mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]);
}
}
@ -166,9 +165,8 @@ void MolarityIonicVPSSTP::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] *= GasConstant * T;
}
/*
* Update the activity coefficients, This also update the
@ -176,9 +174,8 @@ void MolarityIonicVPSSTP::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] -= GasConstant * T * T * dlnActCoeffdT_Scaled_[k];
}
}
@ -209,7 +206,6 @@ void MolarityIonicVPSSTP::getPartialMolarCp(doublereal* cpbar) const
void MolarityIonicVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
{
double xx;
/*
* Get the nondimensional standard state entropies
*/
@ -223,7 +219,7 @@ void MolarityIonicVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
s_update_dlnActCoeff_dT();
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(moleFractions_[k], SmallNumber);
double xx = std::max(moleFractions_[k], SmallNumber);
sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
}
/*
@ -247,69 +243,62 @@ void MolarityIonicVPSSTP::getPartialMolarVolumes(doublereal* vbar) const
void MolarityIonicVPSSTP::calcPseudoBinaryMoleFractions() const
{
size_t k;
size_t kCat;
size_t kMax;
doublereal sumCat;
doublereal sumAnion;
doublereal chP, chM;
doublereal sum = 0.0;
doublereal sumMax;
switch (PBType_) {
case PBTYPE_PASSTHROUGH:
for (k = 0; k < m_kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
PBMoleFractions_[k] = moleFractions_[k];
}
break;
case PBTYPE_SINGLEANION:
sumCat = 0.0;
sumAnion = 0.0;
for (k = 0; k < m_kk; k++) {
{
double sumCat = 0.0;
double sumAnion = 0.0;
for (size_t k = 0; k < m_kk; k++) {
moleFractionsTmp_[k] = moleFractions_[k];
}
kMax = npos;
sumMax = 0.0;
for (k = 0; k < cationList_.size(); k++) {
kCat = cationList_[k];
chP = m_speciesCharge[kCat];
size_t kMax = npos;
double sumMax = 0.0;
for (size_t k = 0; k < cationList_.size(); k++) {
size_t kCat = cationList_[k];
double chP = m_speciesCharge[kCat];
if (moleFractions_[kCat] > sumMax) {
kMax = k;
sumMax = moleFractions_[kCat];
}
sumCat += chP * moleFractions_[kCat];
}
k = anionList_[0];
chM = m_speciesCharge[k];
sumAnion = moleFractions_[k] * chM;
sum = sumCat - sumAnion;
size_t ka = anionList_[0];
sumAnion = moleFractions_[ka] * m_speciesCharge[ka];
double sum = sumCat - sumAnion;
if (fabs(sum) > 1.0E-16) {
moleFractionsTmp_[cationList_[kMax]] -= sum / m_speciesCharge[kMax];
sum = 0.0;
for (k = 0; k < numCationSpecies_; k++) {
for (size_t k = 0; k < numCationSpecies_; k++) {
sum += moleFractionsTmp_[k];
}
for (k = 0; k < numCationSpecies_; k++) {
for (size_t k = 0; k < numCationSpecies_; k++) {
moleFractionsTmp_[k]/= sum;
}
}
for (k = 0; k < numCationSpecies_; k++) {
for (size_t k = 0; k < numCationSpecies_; k++) {
PBMoleFractions_[k] = moleFractionsTmp_[cationList_[k]];
}
for (k = 0; k < numPassThroughSpecies_; k++) {
for (size_t k = 0; k < numPassThroughSpecies_; k++) {
PBMoleFractions_[neutralPBindexStart + k] = moleFractions_[passThroughList_[k]];
}
sum = std::max(0.0, PBMoleFractions_[0]);
for (k = 1; k < numPBSpecies_; k++) {
for (size_t k = 1; k < numPBSpecies_; k++) {
sum += PBMoleFractions_[k];
}
for (k = 0; k < numPBSpecies_; k++) {
for (size_t k = 0; k < numPBSpecies_; k++) {
PBMoleFractions_[k] /= sum;
}
break;
}
case PBTYPE_SINGLECATION:
throw CanteraError("eosType", "Unknown type");
@ -348,13 +337,12 @@ void MolarityIonicVPSSTP::initThermo()
/*
* Go find the list of cations and anions
*/
double ch;
numCationSpecies_ = 0;
cationList_.clear();
anionList_.clear();
passThroughList_.clear();
for (size_t k = 0; k < m_kk; k++) {
ch = m_speciesCharge[k];
double ch = m_speciesCharge[k];
if (ch > 0.0) {
cationList_.push_back(k);
numCationSpecies_++;
@ -387,14 +375,9 @@ void MolarityIonicVPSSTP::initLengths()
void MolarityIonicVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
{
std::string subname = "MolarityIonicVPSSTP::initThermoXML";
std::string stemp;
if ((int) id.size() > 0) {
string idp = phaseNode.id();
if (idp != id) {
throw CanteraError(subname, "phasenode and Id are incompatible");
}
if ((int) id.size() > 0 && phaseNode.id() != id) {
throw CanteraError("MolarityIonicVPSSTP::initThermoXML",
"phasenode and Id are incompatible");
}
/*
@ -403,13 +386,14 @@ void MolarityIonicVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string&
* <thermo model="MolarityIonicVPSSTP" />
*/
if (!phaseNode.hasChild("thermo")) {
throw CanteraError(subname, "no thermo XML node");
throw CanteraError("MolarityIonicVPSSTP::initThermoXML",
"no thermo XML node");
}
XML_Node& thermoNode = phaseNode.child("thermo");
std::string mStringa = thermoNode.attrib("model");
std::string mString = lowercase(mStringa);
if (mString != "molarityionicvpss" && mString != "molarityionicvpsstp") {
throw CanteraError(subname.c_str(),
throw CanteraError("MolarityIonicVPSSTP::initThermoXML",
"Unknown thermo model: " + mStringa + " - This object only knows \"MolarityIonicVPSSTP\" ");
}
@ -417,25 +401,14 @@ void MolarityIonicVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string&
* Go get all of the coefficients and factors in the
* activityCoefficients XML block
*/
XML_Node* acNodePtr = 0;
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& acNode = thermoNode.child("activityCoefficients");
acNodePtr = &acNode;
mStringa = acNode.attrib("model");
mString = lowercase(mStringa);
// if (mString != "redlich-kister") {
// throw CanteraError(subname.c_str(),
// "Unknown activity coefficient model: " + mStringa);
//}
size_t n = acNodePtr->nChildren();
for (size_t i = 0; i < n; i++) {
XML_Node& xmlACChild = acNodePtr->child(i);
stemp = xmlACChild.name();
std::string nodeName = lowercase(stemp);
for (size_t i = 0; i < acNode.nChildren(); i++) {
XML_Node& xmlACChild = acNode.child(i);
/*
* Process a binary interaction
*/
if (nodeName == "binaryneutralspeciesparameters") {
if (lowercase(xmlACChild.name()) == "binaryneutralspeciesparameters") {
readXMLBinarySpecies(xmlACChild);
}
}
@ -475,13 +448,12 @@ std::string MolarityIonicVPSSTP::report(bool show_thermo, doublereal threshold)
sprintf(p, " potential %12.6g V\n", phi);
s += p;
size_t kk = nSpecies();
vector_fp x(kk);
vector_fp molal(kk);
vector_fp mu(kk);
vector_fp muss(kk);
vector_fp acMolal(kk);
vector_fp actMolal(kk);
vector_fp x(m_kk);
vector_fp molal(m_kk);
vector_fp mu(m_kk);
vector_fp muss(m_kk);
vector_fp acMolal(m_kk);
vector_fp actMolal(m_kk);
getMoleFractions(&x[0]);
getChemPotentials(&mu[0]);

89
src/thermo/Mu0Poly.cpp
View file

@ -85,10 +85,8 @@ void Mu0Poly::updateProperties(const doublereal* tt, doublereal* cp_R,
}
double T1 = m_t0_int[j];
double cp_Rj = m_cp0_R_int[j];
doublereal rt = 1.0/T;
cp_R[m_index] = cp_Rj;
h_RT[m_index] = rt*(m_h0_R_int[j] + (T - T1) * cp_Rj);
h_RT[m_index] = (m_h0_R_int[j] + (T - T1) * cp_Rj)/T;
s_R[m_index] = m_s0_R_int[j] + cp_Rj * (log(T/T1));
}
@ -128,13 +126,8 @@ void Mu0Poly::modifyParameters(doublereal* coeffs)
Mu0Poly* newMu0ThermoFromXML(const std::string& speciesName,
const XML_Node& Mu0Node)
{
doublereal tmin, tmax;
bool dimensionlessMu0Values = false;
tmin = fpValue(Mu0Node["Tmin"]);
tmax = fpValue(Mu0Node["Tmax"]);
doublereal pref = fpValue(Mu0Node["Pref"]);
doublereal h298 = 0.0;
if (Mu0Node.hasChild("H298")) {
h298 = getFloat(Mu0Node, "H298", "actEnergy");
@ -158,12 +151,10 @@ Mu0Poly* newMu0ThermoFromXML(const std::string& speciesName,
* form. If they were, then the assumed temperature needs to be
* adjusted from the assumed T = 273.15
*/
string uuu = valNode_ptr->attrib("units");
if (uuu == "Dimensionless") {
if (valNode_ptr->attrib("units") == "Dimensionless") {
dimensionlessMu0Values = true;
}
size_t ns = cValues.size();
if (ns != numPoints) {
if (cValues.size() != numPoints) {
throw CanteraError("installMu0ThermoFromXML",
"numPoints inconsistent while processing "
+ speciesName);
@ -177,8 +168,7 @@ Mu0Poly* newMu0ThermoFromXML(const std::string& speciesName,
+ speciesName);
}
getFloatArray(*tempNode_ptr, cTemperatures, false);
ns = cTemperatures.size();
if (ns != numPoints) {
if (cTemperatures.size() != numPoints) {
throw CanteraError("installMu0ThermoFromXML",
"numPoints inconsistent while processing "
+ speciesName);
@ -203,13 +193,12 @@ Mu0Poly* newMu0ThermoFromXML(const std::string& speciesName,
c[2+i*2+1] = cValues[i];
}
return new Mu0Poly(tmin, tmax, pref, &c[0]);
return new Mu0Poly(fpValue(Mu0Node["Tmin"]), fpValue(Mu0Node["Tmax"]),
fpValue(Mu0Node["Pref"]), &c[0]);
}
void Mu0Poly::processCoeffs(const doublereal* coeffs)
{
size_t i, iindex;
double T1, T2;
size_t nPoints = (size_t) coeffs[0];
if (nPoints < 2) {
throw CanteraError("Mu0Poly",
@ -232,8 +221,8 @@ void Mu0Poly::processCoeffs(const doublereal* coeffs)
* Also distribute the data into the internal arrays.
*/
bool ifound = false;
for (i = 0, iindex = 2; i < nPoints; i++) {
T1 = coeffs[iindex];
for (size_t i = 0, iindex = 2; i < nPoints; i++) {
double T1 = coeffs[iindex];
m_t0_int[i] = T1;
m_mu0_R_int[i] = coeffs[iindex+1] / GasConstant;
if (T1 == 298.15) {
@ -241,8 +230,7 @@ void Mu0Poly::processCoeffs(const doublereal* coeffs)
ifound = true;
}
if (i < nPoints - 1) {
T2 = coeffs[iindex+2];
if (T2 <= T1) {
if (coeffs[iindex+2] <= T1) {
throw CanteraError("Mu0Poly",
"Temperatures are not monotonic increasing");
}
@ -257,26 +245,18 @@ void Mu0Poly::processCoeffs(const doublereal* coeffs)
/*
* Starting from the interval with T298, we go up
*/
doublereal mu2, s1, s2, h1, h2, cpi, deltaMu, deltaT;
T1 = m_t0_int[iT298];
doublereal mu1 = m_mu0_R_int[iT298];
m_h0_R_int[iT298] = m_H298;
m_s0_R_int[iT298] = - (mu1 - m_h0_R_int[iT298]) / T1;
for (i = iT298; i < m_numIntervals; i++) {
T1 = m_t0_int[i];
s1 = m_s0_R_int[i];
h1 = m_h0_R_int[i];
mu1 = m_mu0_R_int[i];
T2 = m_t0_int[i+1];
mu2 = m_mu0_R_int[i+1];
deltaMu = mu2 - mu1;
deltaT = T2 - T1;
cpi = (deltaMu - T1 * s1 + T2 * s1) / (deltaT - T2 * log(T2/T1));
h2 = h1 + cpi * deltaT;
s2 = s1 + cpi * log(T2/T1);
m_cp0_R_int[i] = cpi;
m_h0_R_int[i+1] = h2;
m_s0_R_int[i+1] = s2;
m_s0_R_int[iT298] = - (m_mu0_R_int[iT298] - m_h0_R_int[iT298]) / m_t0_int[iT298];
for (size_t i = iT298; i < m_numIntervals; i++) {
double T1 = m_t0_int[i];
double s1 = m_s0_R_int[i];
double T2 = m_t0_int[i+1];
double deltaMu = m_mu0_R_int[i+1] - m_mu0_R_int[i];
double deltaT = T2 - T1;
double cpi = (deltaMu - T1 * s1 + T2 * s1) / (deltaT - T2 * log(T2/T1));
m_cp0_R_int[i] = cpi;
m_h0_R_int[i+1] = m_h0_R_int[i] + cpi * deltaT;
m_s0_R_int[i+1] = s1 + cpi * log(T2/T1);
m_cp0_R_int[i+1] = cpi;
}
@ -284,25 +264,18 @@ void Mu0Poly::processCoeffs(const doublereal* coeffs)
* Starting from the interval with T298, we go down
*/
if (iT298 != 0) {
T2 = m_t0_int[iT298];
mu2 = m_mu0_R_int[iT298];
m_h0_R_int[iT298] = m_H298;
m_s0_R_int[iT298] = - (mu2 - m_h0_R_int[iT298]) / T2;
for (i = iT298 - 1; i != npos; i--) {
T1 = m_t0_int[i];
mu1 = m_mu0_R_int[i];
T2 = m_t0_int[i+1];
mu2 = m_mu0_R_int[i+1];
s2 = m_s0_R_int[i+1];
h2 = m_h0_R_int[i+1];
deltaMu = mu2 - mu1;
deltaT = T2 - T1;
cpi = (deltaMu - T1 * s2 + T2 * s2) / (deltaT - T1 * log(T2/T1));
h1 = h2 - cpi * deltaT;
s1 = s2 - cpi * log(T2/T1);
m_cp0_R_int[i] = cpi;
m_h0_R_int[i] = h1;
m_s0_R_int[i] = s1;
m_s0_R_int[iT298] = - (m_mu0_R_int[iT298] - m_h0_R_int[iT298]) / m_t0_int[iT298];
for (size_t i = iT298 - 1; i != npos; i--) {
double T1 = m_t0_int[i];
double T2 = m_t0_int[i+1];
double s2 = m_s0_R_int[i+1];
double deltaMu = m_mu0_R_int[i+1] - m_mu0_R_int[i];
double deltaT = T2 - T1;
double cpi = (deltaMu - T1 * s2 + T2 * s2) / (deltaT - T1 * log(T2/T1));
m_cp0_R_int[i] = cpi;
m_h0_R_int[i] = m_h0_R_int[i+1] - cpi * deltaT;
m_s0_R_int[i] = s2 - cpi * log(T2/T1);
if (i == (m_numIntervals-1)) {
m_cp0_R_int[i+1] = cpi;
}

10
src/thermo/NasaThermo.cpp
View file

@ -149,8 +149,6 @@ void NasaThermo::update_one(size_t k, doublereal t, doublereal* cp_R,
void NasaThermo::update(doublereal t, doublereal* cp_R,
doublereal* h_RT, doublereal* s_R) const
{
int i;
// load functions of temperature into m_t vector
m_t[0] = t;
m_t[1] = t*t;
@ -161,7 +159,7 @@ void NasaThermo::update(doublereal t, doublereal* cp_R,
// iterate over the groups
std::vector<NasaPoly1>::const_iterator _begin, _end;
for (i = 0; i != m_ngroups; i++) {
for (int i = 0; i != m_ngroups; i++) {
if (t > m_tmid[i]) {
_begin = m_high[i].begin();
_end = m_high[i].end();
@ -224,15 +222,13 @@ doublereal NasaThermo::reportOneHf298(const size_t k) const
const std::vector<NasaPoly1> &mlg = m_low[grp-1];
const NasaPoly1* nlow = &(mlg[pos]);
doublereal tmid = nlow->maxTemp();
double h;
if (298.15 <= tmid) {
h = nlow->reportHf298(0);
return nlow->reportHf298(0);
} else {
const std::vector<NasaPoly1> &mhg = m_high[grp-1];
const NasaPoly1* nhigh = &(mhg[pos]);
h = nhigh->reportHf298(0);
return nhigh->reportHf298(0);
}
return h;
}
void NasaThermo::modifyOneHf298(const size_t k, const doublereal Hf298New)

14
src/thermo/PDSS.cpp
View file

@ -295,29 +295,23 @@ doublereal PDSS::molarVolume_ref() const
doublereal PDSS::enthalpyDelp_mole() const
{
doublereal RT = m_temp * GasConstant;
doublereal tmp = enthalpy_RT_ref();
return enthalpy_mole() - RT * tmp;
return enthalpy_mole() - m_temp * GasConstant * enthalpy_RT_ref();
}
doublereal PDSS::entropyDelp_mole() const
{
doublereal tmp = entropy_R_ref();
return entropy_mole() - GasConstant * tmp;
return entropy_mole() - GasConstant * entropy_R_ref();
}
doublereal PDSS::gibbsDelp_mole() const
{
doublereal RT = m_temp * GasConstant;
doublereal tmp = gibbs_RT_ref();
return gibbs_mole() - RT * tmp;
return gibbs_mole() - m_temp * GasConstant * gibbs_RT_ref();
}
doublereal PDSS::cpDelp_mole() const
{
doublereal tmp = cp_R_ref();
return cp_mole() - GasConstant * tmp;
return cp_mole() - GasConstant * cp_R_ref();
}
doublereal PDSS::pressure() const

26
src/thermo/PDSS_ConstVol.cpp
View file

@ -72,8 +72,7 @@ void PDSS_ConstVol::constructPDSSXML(VPStandardStateTP* tp, size_t spindex,
const XML_Node& phaseNode, bool spInstalled)
{
PDSS::initThermo();
SpeciesThermo& sp = m_tp->speciesThermo();
m_p0 = sp.refPressure(m_spindex);
m_p0 = m_tp->speciesThermo().refPressure(m_spindex);
if (!spInstalled) {
throw CanteraError("PDSS_ConstVol::constructPDSSXML", "spInstalled false not handled");
@ -84,15 +83,12 @@ void PDSS_ConstVol::constructPDSSXML(VPStandardStateTP* tp, size_t spindex,
throw CanteraError("PDSS_ConstVol::constructPDSSXML",
"no standardState Node for species " + speciesNode.name());
}
std::string model = ss->attrib("model");
if (model != "constant_incompressible") {
if (ss->attrib("model") != "constant_incompressible") {
throw CanteraError("PDSS_ConstVol::initThermoXML",
"standardState model for species isn't constant_incompressible: " + speciesNode.name());
}
m_constMolarVolume = ctml::getFloat(*ss, "molarVolume", "toSI");
std::string id = "";
}
void PDSS_ConstVol::constructPDSSFile(VPStandardStateTP* tp, size_t spindex,
@ -126,8 +122,7 @@ void PDSS_ConstVol::constructPDSSFile(VPStandardStateTP* tp, size_t spindex,
XML_Node& speciesList = fxml_phase->child("speciesArray");
XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"],
&(fxml_phase->root()));
const vector<string>&sss = tp->speciesNames();
const XML_Node* s = speciesDB->findByAttr("name", sss[spindex]);
const XML_Node* s = speciesDB->findByAttr("name", tp->speciesName(spindex));
constructPDSSXML(tp, spindex, *s, *fxml_phase, true);
delete fxml;
@ -145,8 +140,7 @@ void PDSS_ConstVol::initThermoXML(const XML_Node& phaseNode, const std::string&
void PDSS_ConstVol::initThermo()
{
PDSS::initThermo();
SpeciesThermo& sp = m_tp->speciesThermo();
m_p0 = sp.refPressure(m_spindex);
m_p0 = m_tp->speciesThermo().refPressure(m_spindex);
m_V0_ptr[m_spindex] = m_constMolarVolume;
m_Vss_ptr[m_spindex] = m_constMolarVolume;
}
@ -160,10 +154,8 @@ PDSS_ConstVol::enthalpy_RT() const
doublereal
PDSS_ConstVol::intEnergy_mole() const
{
doublereal pVRT = (m_pres * m_Vss_ptr[m_spindex]) / (GasConstant * m_temp);
doublereal val = m_h0_RT_ptr[m_spindex] - pVRT;
doublereal RT = GasConstant * m_temp;
return val * RT;
doublereal pV = (m_pres * m_Vss_ptr[m_spindex]);
return m_h0_RT_ptr[m_spindex] * GasConstant * m_temp - pV;
}
doublereal
@ -199,8 +191,7 @@ PDSS_ConstVol::molarVolume() const
doublereal
PDSS_ConstVol::density() const
{
doublereal val = m_Vss_ptr[m_spindex];
return m_mw/val;
return m_mw / m_Vss_ptr[m_spindex];
}
doublereal
@ -221,8 +212,7 @@ doublereal PDSS_ConstVol::entropy_R_ref() const
doublereal PDSS_ConstVol::cp_R_ref() const
{
doublereal val = m_cp0_R_ptr[m_spindex];
return (val);
return m_cp0_R_ptr[m_spindex];
}
doublereal PDSS_ConstVol::molarVolume_ref() const

90
src/thermo/PDSS_HKFT.cpp
View file

@ -212,9 +212,7 @@ PDSS* PDSS_HKFT::duplMyselfAsPDSS() const
doublereal PDSS_HKFT::enthalpy_mole() const
{
// Ok we may change this evaluation method in the future.
doublereal GG = gibbs_mole();
doublereal SS = entropy_mole();
doublereal h = GG + m_temp * SS;
doublereal h = gibbs_mole() + m_temp * entropy_mole();
#ifdef DEBUG_MODE_NOT
doublereal h2 = enthalpy_mole2();
@ -230,29 +228,24 @@ doublereal PDSS_HKFT::enthalpy_mole() const
#ifdef DEBUG_MODE
doublereal PDSS_HKFT::enthalpy_mole2() const
{
doublereal delH = deltaH();
double enthTRPR = m_Mu0_tr_pr + 298.15 * m_Entrop_tr_pr * 1.0E3 * 4.184;
return delH + enthTRPR;
return deltaH() + enthTRPR;
}
#endif
doublereal PDSS_HKFT::intEnergy_mole() const
{
doublereal hh = enthalpy_RT();
doublereal mv = molarVolume();
return hh - mv * m_pres;
return enthalpy_RT() - molarVolume() * m_pres;
}
doublereal PDSS_HKFT::entropy_mole() const
{
doublereal delS = deltaS();
return m_Entrop_tr_pr * 1.0E3 * 4.184 + delS;
return m_Entrop_tr_pr * 1.0E3 * 4.184 + deltaS();
}
doublereal PDSS_HKFT::gibbs_mole() const
{
doublereal delG = deltaG();
return m_Mu0_tr_pr + delG;
return m_Mu0_tr_pr + deltaG();
}
doublereal PDSS_HKFT::cp_mole() const
@ -407,8 +400,7 @@ doublereal PDSS_HKFT::molarVolume() const
doublereal
PDSS_HKFT::density() const
{
doublereal val = molarVolume();
return m_mw/val;
return m_mw / molarVolume();
}
doublereal
@ -570,8 +562,7 @@ void PDSS_HKFT::constructPDSSXML(VPStandardStateTP* tp, size_t spindex,
throw CanteraError("PDSS_HKFT::constructPDSSXML",
"no thermo Node for species " + speciesNode.name());
}
std::string model = lowercase(tn->attrib("model"));
if (model != "hkft") {
if (lowercase(tn->attrib("model")) != "hkft") {
throw CanteraError("PDSS_HKFT::initThermoXML",
"thermo model for species isn't hkft: "
+ speciesNode.name());
@ -622,52 +613,44 @@ void PDSS_HKFT::constructPDSSXML(VPStandardStateTP* tp, size_t spindex,
throw CanteraError("PDSS_HKFT::constructPDSSXML",
"no standardState Node for species " + speciesNode.name());
}
model = lowercase(ss->attrib("model"));
if (model != "hkft") {
if (lowercase(ss->attrib("model")) != "hkft") {
throw CanteraError("PDSS_HKFT::initThermoXML",
"standardState model for species isn't hkft: "
+ speciesNode.name());
}
if (ss->hasChild("a1")) {
doublereal val = getFloat(*ss, "a1");
m_a1 = val;
m_a1 = getFloat(*ss, "a1");
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a1 field");
}
if (ss->hasChild("a2")) {
doublereal val = getFloat(*ss, "a2");
m_a2 = val;
m_a2 = getFloat(*ss, "a2");
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a2 field");
}
if (ss->hasChild("a3")) {
doublereal val = getFloat(*ss, "a3");
m_a3 = val;
m_a3 = getFloat(*ss, "a3");
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a3 field");
}
if (ss->hasChild("a4")) {
doublereal val = getFloat(*ss, "a4");
m_a4 = val;
m_a4 = getFloat(*ss, "a4");
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing a4 field");
}
if (ss->hasChild("c1")) {
doublereal val = getFloat(*ss, "c1");
m_c1 = val;
m_c1 = getFloat(*ss, "c1");
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing c1 field");
}
if (ss->hasChild("c2")) {
doublereal val = getFloat(*ss, "c2");
m_c2 = val;
m_c2 = getFloat(*ss, "c2");
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing c2 field");
}
if (ss->hasChild("omega_Pr_Tr")) {
doublereal val = getFloat(*ss, "omega_Pr_Tr");
m_omega_pr_tr = val;
m_omega_pr_tr = getFloat(*ss, "omega_Pr_Tr");
} else {
throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing omega_Pr_Tr field");
}
@ -738,8 +721,7 @@ void PDSS_HKFT::constructPDSSFile(VPStandardStateTP* tp, size_t spindex,
XML_Node& speciesList = fxml_phase->child("speciesArray");
XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"],
&(fxml_phase->root()));
const vector<string>&sss = tp->speciesNames();
const XML_Node* s = speciesDB->findByAttr("name", sss[spindex]);
const XML_Node* s = speciesDB->findByAttr("name", tp->speciesName(spindex));
constructPDSSXML(tp, spindex, *s, *fxml_phase, true);
delete fxml;
@ -913,8 +895,7 @@ doublereal PDSS_HKFT::ag(const doublereal temp, const int ifunc) const
{
static doublereal ag_coeff[3] = { -2.037662, 5.747000E-3, -6.557892E-6};
if (ifunc == 0) {
doublereal t2 = temp * temp;
return ag_coeff[0] + ag_coeff[1] * temp + ag_coeff[2] * t2;
return ag_coeff[0] + ag_coeff[1] * temp + ag_coeff[2] * temp * temp;
} else if (ifunc == 1) {
return ag_coeff[1] + ag_coeff[2] * 2.0 * temp;
}
@ -928,8 +909,7 @@ doublereal PDSS_HKFT::bg(const doublereal temp, const int ifunc) const
{
static doublereal bg_coeff[3] = { 6.107361, -1.074377E-2, 1.268348E-5};
if (ifunc == 0) {
doublereal t2 = temp * temp;
return bg_coeff[0] + bg_coeff[1] * temp + bg_coeff[2] * t2;
return bg_coeff[0] + bg_coeff[1] * temp + bg_coeff[2] * temp * temp;
} else if (ifunc == 1) {
return bg_coeff[1] + bg_coeff[2] * 2.0 * temp;
}
@ -955,29 +935,23 @@ doublereal PDSS_HKFT::f(const doublereal temp, const doublereal pres, const int
doublereal T1 = (TC-155.0)/300.;
doublereal fac1;
doublereal p2 = (1000. - presBar) * (1000. - presBar);
doublereal p3 = (1000. - presBar) * p2;
doublereal p4 = p2 * p2;
doublereal fac2 = af_coeff[1] * p3 + af_coeff[2] * p4;
if (ifunc == 0) {
fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0);
return fac1 * fac2;
return pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0) * fac2;
} else if (ifunc == 1) {
fac1 = (4.8 * pow(T1,3.8) + 16.0 * af_coeff[0] * pow(T1, 15.0)) / 300.;
return fac1 * fac2;
return (4.8 * pow(T1,3.8) + 16.0 * af_coeff[0] * pow(T1, 15.0)) / 300. * fac2;
} else if (ifunc == 2) {
fac1 = (4.8 * 3.8 * pow(T1,2.8) + 16.0 * 15.0 * af_coeff[0] * pow(T1, 14.0)) / (300. * 300.);
return fac1 * fac2;
return (4.8 * 3.8 * pow(T1,2.8) + 16.0 * 15.0 * af_coeff[0] * pow(T1, 14.0)) / (300. * 300.) * fac2;
} else if (ifunc == 3) {
fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0);
double fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0);
fac2 = - (3.0 * af_coeff[1] * p2 + 4.0 * af_coeff[2] * p3)/ 1.0E5;
return fac1 * fac2;
} else {
throw CanteraError("HKFT_PDSS::gg", "unimplemented");
}
return 0.0;
}
doublereal PDSS_HKFT::g(const doublereal temp, const doublereal pres, const int ifunc) const
@ -1075,8 +1049,7 @@ doublereal PDSS_HKFT::LookupGe(const std::string& elemName)
throw CanteraError("PDSS_HKFT::LookupGe",
"element " + elemName + " does not have a supplied entropy298");
}
geValue *= (-298.15);
return geValue;
return geValue * -298.15;
}
void PDSS_HKFT::convertDGFormation()
@ -1084,25 +1057,16 @@ void PDSS_HKFT::convertDGFormation()
/*
* Ok let's get the element compositions and conversion factors.
*/
size_t ne = m_tp->nElements();
doublereal na;
doublereal ge;
string ename;
doublereal totalSum = 0.0;
for (size_t m = 0; m < ne; m++) {
na = m_tp->nAtoms(m_spindex, m);
for (size_t m = 0; m < m_tp->nElements(); m++) {
double na = m_tp->nAtoms(m_spindex, m);
if (na > 0.0) {
ename = m_tp->elementName(m);
ge = LookupGe(ename);
totalSum += na * ge;
totalSum += na * LookupGe(m_tp->elementName(m));
}
}
// Add in the charge
if (m_charge_j != 0.0) {
ename = "H";
ge = LookupGe(ename);
totalSum -= m_charge_j * ge;
totalSum -= m_charge_j * LookupGe("H");
}
// Ok, now do the calculation. Convert to joules kmol-1
doublereal dg = m_deltaG_formation_tr_pr * 4.184 * 1.0E3;

11
src/thermo/PDSS_IdealGas.cpp
View file

@ -40,8 +40,7 @@ PDSS_IdealGas::PDSS_IdealGas(VPStandardStateTP* tp, size_t spindex, const XML_No
throw CanteraError("PDSS_IdealGas", "sp installing not done yet");
}
m_pdssType = cPDSS_IDEALGAS;
std::string id = "";
constructPDSSXML(tp, spindex, phaseRoot, id);
constructPDSSXML(tp, spindex, phaseRoot, "");
}
PDSS_IdealGas::PDSS_IdealGas(const PDSS_IdealGas& b) :
@ -107,8 +106,7 @@ void PDSS_IdealGas::constructPDSSFile(VPStandardStateTP* tp, size_t spindex,
void PDSS_IdealGas::initThermo()
{
PDSS::initThermo();
SpeciesThermo& sp = m_tp->speciesThermo();
m_p0 = sp.refPressure(m_spindex);
m_p0 = m_tp->speciesThermo().refPressure(m_spindex);
m_minTemp = m_spthermo->minTemp(m_spindex);
m_maxTemp = m_spthermo->maxTemp(m_spindex);
}
@ -122,9 +120,7 @@ PDSS_IdealGas::enthalpy_RT() const
doublereal
PDSS_IdealGas::intEnergy_mole() const
{
doublereal val = m_h0_RT_ptr[m_spindex] - 1.0;
doublereal RT = GasConstant * m_temp;
return val * RT;
return (m_h0_RT_ptr[m_spindex] - 1.0) * GasConstant * m_temp;
}
doublereal
@ -193,7 +189,6 @@ doublereal PDSS_IdealGas::molarVolume_ref() const
doublereal PDSS_IdealGas::pressure() const
{
throw CanteraError("PDSS_IdealGas::pressure()", "unimplemented");
return 0.0;
}
void PDSS_IdealGas::setPressure(doublereal p)

14
src/thermo/PDSS_IonsFromNeutral.cpp
View file

@ -119,8 +119,7 @@ void PDSS_IonsFromNeutral::constructPDSSXML(VPStandardStateTP* tp, size_t spinde
throw CanteraError("PDSS_IonsFromNeutral::constructPDSSXML",
"no thermo Node for species " + speciesNode.name());
}
std::string model = lowercase(tn->attrib("model"));
if (model != "ionfromneutral") {
if (lowercase(tn->attrib("model")) != "ionfromneutral") {
throw CanteraError("PDSS_IonsFromNeutral::constructPDSSXML",
"thermo model for species isn't IonsFromNeutral: "
+ speciesNode.name());
@ -194,9 +193,7 @@ void PDSS_IonsFromNeutral::constructPDSSFile(VPStandardStateTP* tp, size_t spind
XML_Node& speciesList = fxml_phase->child("speciesArray");
XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"],
&(fxml_phase->root()));
const vector<string>&sss = tp->speciesNames();
const XML_Node* s = speciesDB->findByAttr("name", sss[spindex]);
const XML_Node* s = speciesDB->findByAttr("name", tp->speciesName(spindex));
constructPDSSXML(tp, spindex, *s, *fxml_phase, id);
delete fxml;
@ -205,8 +202,7 @@ void PDSS_IonsFromNeutral::constructPDSSFile(VPStandardStateTP* tp, size_t spind
void PDSS_IonsFromNeutral::initThermo()
{
PDSS::initThermo();
SpeciesThermo& sp = m_tp->speciesThermo();
m_p0 = sp.refPressure(m_spindex);
m_p0 = m_tp->speciesThermo().refPressure(m_spindex);
m_minTemp = m_spthermo->minTemp(m_spindex);
m_maxTemp = m_spthermo->maxTemp(m_spindex);
}
@ -226,9 +222,7 @@ PDSS_IonsFromNeutral::enthalpy_RT() const
doublereal
PDSS_IonsFromNeutral::intEnergy_mole() const
{
doublereal val = m_h0_RT_ptr[m_spindex] - 1.0;
doublereal RT = GasConstant * m_temp;
return val * RT;
return (m_h0_RT_ptr[m_spindex] - 1.0) * GasConstant * m_temp;
}
doublereal

19
src/thermo/PDSS_SSVol.cpp
View file

@ -86,8 +86,7 @@ void PDSS_SSVol::constructPDSSXML(VPStandardStateTP* tp, size_t spindex,
const XML_Node& phaseNode, bool spInstalled)
{
PDSS::initThermo();
SpeciesThermo& sp = m_tp->speciesThermo();
m_p0 = sp.refPressure(m_spindex);
m_p0 = m_tp->speciesThermo().refPressure(m_spindex);
if (!spInstalled) {
throw CanteraError("PDSS_SSVol::constructPDSSXML", "spInstalled false not handled");
@ -120,7 +119,6 @@ void PDSS_SSVol::constructPDSSXML(VPStandardStateTP* tp, size_t spindex,
throw CanteraError("PDSS_SSVol::constructPDSSXML",
"standardState model for species isn't constant_incompressible: " + speciesNode.name());
}
std::string id = "";
}
void PDSS_SSVol::constructPDSSFile(VPStandardStateTP* tp, size_t spindex,
@ -153,8 +151,7 @@ void PDSS_SSVol::constructPDSSFile(VPStandardStateTP* tp, size_t spindex,
XML_Node& speciesList = fxml_phase->child("speciesArray");
XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"],
&(fxml_phase->root()));
const vector<string>&sss = tp->speciesNames();
const XML_Node* s = speciesDB->findByAttr("name", sss[spindex]);
const XML_Node* s = speciesDB->findByAttr("name", tp->speciesName(spindex));
constructPDSSXML(tp, spindex, *s, *fxml_phase, true);
delete fxml;
@ -172,8 +169,7 @@ void PDSS_SSVol::initThermoXML(const XML_Node& phaseNode, const std::string& id)
void PDSS_SSVol::initThermo()
{
PDSS::initThermo();
SpeciesThermo& sp = m_tp->speciesThermo();
m_p0 = sp.refPressure(m_spindex);
m_p0 = m_tp->speciesThermo().refPressure(m_spindex);
m_V0_ptr[m_spindex] = m_constMolarVolume;
m_Vss_ptr[m_spindex] = m_constMolarVolume;
}
@ -187,10 +183,8 @@ PDSS_SSVol::enthalpy_RT() const
doublereal
PDSS_SSVol::intEnergy_mole() const
{
doublereal pVRT = (m_pres * m_Vss_ptr[m_spindex]) / (GasConstant * m_temp);
doublereal val = m_h0_RT_ptr[m_spindex] - pVRT;
doublereal RT = GasConstant * m_temp;
return val * RT;
doublereal pV = m_pres * m_Vss_ptr[m_spindex];
return m_h0_RT_ptr[m_spindex] * GasConstant * m_temp - pV;
}
doublereal
@ -226,8 +220,7 @@ PDSS_SSVol::molarVolume() const
doublereal
PDSS_SSVol::density() const
{
doublereal val = m_Vss_ptr[m_spindex];
return m_mw/val;
return m_mw / m_Vss_ptr[m_spindex];
}
doublereal

20
src/thermo/PDSS_Water.cpp
View file

@ -213,26 +213,22 @@ void PDSS_Water::constructSet()
doublereal PDSS_Water::enthalpy_mole() const
{
doublereal h = m_sub.enthalpy();
return h + EW_Offset;
return m_sub.enthalpy() + EW_Offset;
}
doublereal PDSS_Water::intEnergy_mole() const
{
doublereal u = m_sub.intEnergy();
return u + EW_Offset;
return m_sub.intEnergy() + EW_Offset;
}
doublereal PDSS_Water::entropy_mole() const
{
doublereal s = m_sub.entropy();
return s + SW_Offset;
return m_sub.entropy() + SW_Offset;
}
doublereal PDSS_Water::gibbs_mole() const
{
doublereal g = m_sub.Gibbs();
return g + EW_Offset - SW_Offset*m_temp;
return m_sub.Gibbs() + EW_Offset - SW_Offset*m_temp;
}
doublereal PDSS_Water::cp_mole() const
@ -297,9 +293,8 @@ doublereal PDSS_Water::molarVolume_ref() const
doublereal PDSS_Water::pressure() const
{
doublereal p = m_sub.pressure();
m_pres = p;
return p;
m_pres = m_sub.pressure();
return m_pres;
}
void PDSS_Water::setPressure(doublereal p)
@ -388,8 +383,7 @@ doublereal PDSS_Water::density() const
void PDSS_Water::setTemperature(doublereal temp)
{
m_temp = temp;
doublereal dd = m_dens;
m_sub.setState_TR(temp, dd);
m_sub.setState_TR(temp, m_dens);
}
void PDSS_Water::setState_TP(doublereal temp, doublereal pres)

39
src/thermo/Phase.cpp
View file

@ -96,10 +96,9 @@ Phase& Phase::operator=(const Phase& right)
XML_Node *rroot = &(right.m_xml->root());
XML_Node *root_xml = new XML_Node();
(rroot)->copy(root_xml);
string iidd = right.m_xml->id();
m_xml = findXMLPhase(root_xml, iidd);
m_xml = findXMLPhase(root_xml, right.m_xml->id());
if (!m_xml) {
throw CanteraError("Phase::operator=()", "Confused: Couldn't find original phase " + iidd);
throw CanteraError("Phase::operator=()", "Confused: Couldn't find original phase " + right.m_xml->id());
}
if (&(m_xml->root()) != root_xml) {
throw CanteraError("Phase::operator=()", "confused: root changed");
@ -131,13 +130,12 @@ void Phase::setXMLdata(XML_Node& xmlPhase)
XML_Node* xroot = &(xmlPhase.root());
XML_Node *root_xml = new XML_Node();
(xroot)->copy(root_xml);
std::string iidd = xmlPhase.id();
if (m_xml) {
XML_Node *rOld = &(m_xml->root());
delete rOld;
m_xml = 0;
}
m_xml = findXMLPhase(root_xml, iidd);
m_xml = findXMLPhase(root_xml, xmlPhase.id());
if (!m_xml) {
throw CanteraError("Phase::setXMLdata()", "XML 'phase' node not found");
}
@ -300,8 +298,7 @@ void Phase::checkSpeciesArraySize(size_t kk) const
std::string Phase::speciesSPName(int k) const
{
std::string sn = speciesName(k);
return m_name + ":" + sn;
return m_name + ":" + speciesName(k);
}
void Phase::saveState(vector_fp& state) const
@ -372,8 +369,7 @@ void Phase::setMoleFractions(const doublereal* const x)
void Phase::setMoleFractions_NoNorm(const doublereal* const x)
{
m_mmw = dot(x, x + m_kk, m_molwts.begin());
doublereal rmmw = 1.0/m_mmw;
transform(x, x + m_kk, m_ym.begin(), timesConstant<double>(rmmw));
transform(x, x + m_kk, m_ym.begin(), timesConstant<double>(1.0/m_mmw));
transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(),
m_y.begin(), multiplies<double>());
m_stateNum++;
@ -381,9 +377,8 @@ void Phase::setMoleFractions_NoNorm(const doublereal* const x)
void Phase::setMoleFractionsByName(const compositionMap& xMap)
{
size_t kk = nSpecies();
vector_fp mf(kk, 0.0);
for (size_t k = 0; k < kk; k++) {
vector_fp mf(m_kk, 0.0);
for (size_t k = 0; k < m_kk; k++) {
mf[k] = std::max(getValue(xMap, speciesName(k), 0.0), 0.0);
}
setMoleFractions(&mf[0]);
@ -391,8 +386,7 @@ void Phase::setMoleFractionsByName(const compositionMap& xMap)
void Phase::setMoleFractionsByName(const std::string& x)
{
compositionMap c = parseCompString(x, speciesNames());
setMoleFractionsByName(c);
setMoleFractionsByName(parseCompString(x, speciesNames()));
}
void Phase::setMassFractions(const doublereal* const y)
@ -422,9 +416,8 @@ void Phase::setMassFractions_NoNorm(const doublereal* const y)
void Phase::setMassFractionsByName(const compositionMap& yMap)
{
size_t kk = nSpecies();
vector_fp mf(kk, 0.0);
for (size_t k = 0; k < kk; k++) {
vector_fp mf(m_kk, 0.0);
for (size_t k = 0; k < m_kk; k++) {
mf[k] = std::max(getValue(yMap, speciesName(k), 0.0), 0.0);
}
setMassFractions(&mf[0]);
@ -432,8 +425,7 @@ void Phase::setMassFractionsByName(const compositionMap& yMap)
void Phase::setMassFractionsByName(const std::string& y)
{
compositionMap c = parseCompString(y, speciesNames());
setMassFractionsByName(c);
setMassFractionsByName(parseCompString(y, speciesNames()));
}
void Phase::setState_TRX(doublereal t, doublereal dens, const doublereal* x)
@ -534,8 +526,7 @@ void Phase::getMoleFractionsByName(compositionMap& x) const
" 'compositionMap getMoleFractionsByName(double threshold)'"
" instead");
x.clear();
size_t kk = nSpecies();
for (size_t k = 0; k < kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
x[speciesName(k)] = Phase::moleFraction(k);
}
}
@ -659,13 +650,11 @@ doublereal Phase::molarVolume() const
doublereal Phase::chargeDensity() const
{
size_t kk = nSpecies();
doublereal cdens = 0.0;
for (size_t k = 0; k < kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
cdens += charge(k)*moleFraction(k);
}
cdens *= Faraday;
return cdens;
return cdens * Faraday;
}
doublereal Phase::mean_X(const doublereal* const Q) const

220
src/thermo/PhaseCombo_Interaction.cpp
View file

@ -184,7 +184,6 @@ void PhaseCombo_Interaction::getElectrochemPotentials(doublereal* mu) const
void PhaseCombo_Interaction::getChemPotentials(doublereal* mu) const
{
doublereal xx;
/*
* First get the standard chemical potentials in
* molar form.
@ -199,18 +198,17 @@ void PhaseCombo_Interaction::getChemPotentials(doublereal* mu) const
doublereal RT = GasConstant * temperature();
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(moleFractions_[k], SmallNumber);
double xx = std::max(moleFractions_[k], SmallNumber);
mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]);
}
}
doublereal PhaseCombo_Interaction::enthalpy_mole() const
{
size_t kk = nSpecies();
double h = 0;
vector_fp hbar(kk);
vector_fp hbar(m_kk);
getPartialMolarEnthalpies(&hbar[0]);
for (size_t i = 0; i < kk; i++) {
for (size_t i = 0; i < m_kk; i++) {
h += moleFractions_[i]*hbar[i];
}
return h;
@ -218,11 +216,10 @@ doublereal PhaseCombo_Interaction::enthalpy_mole() const
doublereal PhaseCombo_Interaction::entropy_mole() const
{
size_t kk = nSpecies();
double s = 0;
vector_fp sbar(kk);
vector_fp sbar(m_kk);
getPartialMolarEntropies(&sbar[0]);
for (size_t i = 0; i < kk; i++) {
for (size_t i = 0; i < m_kk; i++) {
s += moleFractions_[i]*sbar[i];
}
return s;
@ -230,11 +227,10 @@ doublereal PhaseCombo_Interaction::entropy_mole() const
doublereal PhaseCombo_Interaction::cp_mole() const
{
size_t kk = nSpecies();
double cp = 0;
vector_fp cpbar(kk);
vector_fp cpbar(m_kk);
getPartialMolarCp(&cpbar[0]);
for (size_t i = 0; i < kk; i++) {
for (size_t i = 0; i < m_kk; i++) {
cp += moleFractions_[i]*cpbar[i];
}
return cp;
@ -255,9 +251,8 @@ void PhaseCombo_Interaction::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] *= GasConstant * T;
}
/*
* Update the activity coefficients, This also update the
@ -265,9 +260,8 @@ void PhaseCombo_Interaction::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] -= GasConstant * T * T * dlnActCoeffdT_Scaled_[k];
}
}
@ -298,7 +292,6 @@ void PhaseCombo_Interaction::getPartialMolarCp(doublereal* cpbar) const
void PhaseCombo_Interaction::getPartialMolarEntropies(doublereal* sbar) const
{
double xx;
/*
* Get the nondimensional standard state entropies
*/
@ -312,7 +305,7 @@ void PhaseCombo_Interaction::getPartialMolarEntropies(doublereal* sbar) const
s_update_dlnActCoeff_dT();
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(moleFractions_[k], SmallNumber);
double xx = std::max(moleFractions_[k], SmallNumber);
sbar[k] += - lnActCoeff_Scaled_[k] - log(xx) - T * dlnActCoeffdT_Scaled_[k];
}
/*
@ -325,8 +318,6 @@ void PhaseCombo_Interaction::getPartialMolarEntropies(doublereal* sbar) const
void PhaseCombo_Interaction::getPartialMolarVolumes(doublereal* vbar) const
{
int delAK, delBK;
double XA, XB, g0, g1;
double T = temperature();
/*
@ -335,8 +326,8 @@ void PhaseCombo_Interaction::getPartialMolarVolumes(doublereal* vbar) const
getStandardVolumes(vbar);
for (size_t iK = 0; iK < m_kk; iK++) {
delAK = 0;
delBK = 0;
int delAK = 0;
int delBK = 0;
for (size_t i = 0; i < numBinaryInteractions_; i++) {
size_t iA = m_pSpecies_A_ij[i];
@ -348,11 +339,11 @@ void PhaseCombo_Interaction::getPartialMolarVolumes(doublereal* vbar) const
delBK = 1;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
g0 = (m_VHE_b_ij[i] - T * m_VSE_b_ij[i]);
g1 = (m_VHE_c_ij[i] - T * m_VSE_c_ij[i]);
double g0 = (m_VHE_b_ij[i] - T * m_VSE_b_ij[i]);
double g1 = (m_VHE_c_ij[i] - T * m_VSE_c_ij[i]);
vbar[iK] += XA*XB*(g0+g1*XB)+((delAK-XA)*XB+XA*(delBK-XB))*(g0+g1*XB)+XA*XB*(delBK-XB)*g1;
}
@ -373,14 +364,9 @@ void PhaseCombo_Interaction::initLengths()
void PhaseCombo_Interaction::initThermoXML(XML_Node& phaseNode, const std::string& id)
{
string subname = "PhaseCombo_Interaction::initThermoXML";
string stemp;
if ((int) id.size() > 0) {
string idp = phaseNode.id();
if (idp != id) {
throw CanteraError(subname,
"phasenode and Id are incompatible");
}
if ((int) id.size() > 0 && phaseNode.id() != id) {
throw CanteraError("PhaseCombo_Interaction::initThermoXML",
"phasenode and Id are incompatible");
}
/*
@ -388,14 +374,13 @@ void PhaseCombo_Interaction::initThermoXML(XML_Node& phaseNode, const std::strin
* <thermo model="PhaseCombo_Interaction" />
*/
if (!phaseNode.hasChild("thermo")) {
throw CanteraError(subname,
throw CanteraError("PhaseCombo_Interaction::initThermoXML",
"no thermo XML node");
}
XML_Node& thermoNode = phaseNode.child("thermo");
stemp = thermoNode.attrib("model");
string formString = lowercase(stemp);
string formString = lowercase(thermoNode.attrib("model"));
if (formString != "phasecombo_interaction") {
throw CanteraError(subname,
throw CanteraError("PhaseCombo_Interaction::initThermoXML",
"model name isn't PhaseCombo_Interaction: " + formString);
}
@ -403,27 +388,21 @@ void PhaseCombo_Interaction::initThermoXML(XML_Node& phaseNode, const std::strin
* Go get all of the coefficients and factors in the
* activityCoefficients XML block
*/
XML_Node* acNodePtr = 0;
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& acNode = thermoNode.child("activityCoefficients");
acNodePtr = &acNode;
string mStringa = acNode.attrib("model");
string mString = lowercase(mStringa);
if (mString != "margules") {
throw CanteraError(subname.c_str(),
"Unknown activity coefficient model: " + mStringa);
string mString = acNode.attrib("model");
if (lowercase(mString) != "margules") {
throw CanteraError("PhaseCombo_Interaction::initThermoXML",
"Unknown activity coefficient model: " + mString);
}
size_t n = acNodePtr->nChildren();
for (size_t i = 0; i < n; i++) {
XML_Node& xmlACChild = acNodePtr->child(i);
stemp = xmlACChild.name();
string nodeName = lowercase(stemp);
for (size_t i = 0; i < acNode.nChildren(); i++) {
XML_Node& xmlACChild = acNode.child(i);
/*
* Process a binary salt field, or any of the other XML fields
* that make up the Pitzer Database. Entries will be ignored
* if any of the species in the entry isn't in the solution.
*/
if (nodeName == "binaryneutralspeciesparameters") {
if (lowercase(xmlACChild.name()) == "binaryneutralspeciesparameters") {
readXMLBinarySpecies(xmlACChild);
}
@ -440,18 +419,14 @@ void PhaseCombo_Interaction::initThermoXML(XML_Node& phaseNode, const std::strin
void PhaseCombo_Interaction::s_update_lnActCoeff() const
{
int delAK, delBK;
doublereal XA, XB, g0 , g1;
doublereal xx;
doublereal T = temperature();
doublereal RT = GasConstant*T;
lnActCoeff_Scaled_.assign(m_kk, 0.0);
for (size_t iK = 0; iK < m_kk; iK++) {
/*
* We never sample the end of the mole fraction domains
*/
xx = std::max(moleFractions_[iK], SmallNumber);
double xx = std::max(moleFractions_[iK], SmallNumber);
/*
* First wipe out the ideal solution mixing term
*/
@ -463,17 +438,17 @@ void PhaseCombo_Interaction::s_update_lnActCoeff() const
for (size_t i = 0; i < numBinaryInteractions_; i++) {
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
delAK = 0;
delBK = 0;
int delAK = 0;
int delBK = 0;
if (iA==iK) {
delAK = 1;
} else if (iB==iK) {
delBK = 1;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
double g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / (GasConstant*T);
double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / (GasConstant*T);
lnActCoeff_Scaled_[iK] += (delAK * XB + XA * delBK - XA * XB) * (g0 + g1 * XB) + XA * XB * (delBK - XB) * g1;
}
}
@ -481,27 +456,24 @@ void PhaseCombo_Interaction::s_update_lnActCoeff() const
void PhaseCombo_Interaction::s_update_dlnActCoeff_dT() const
{
int delAK, delBK;
doublereal XA, XB, g0, g1;
doublereal T = temperature();
doublereal RTT = GasConstant*T*T;
dlnActCoeffdT_Scaled_.assign(m_kk, 0.0);
d2lnActCoeffdT2_Scaled_.assign(m_kk, 0.0);
for (size_t iK = 0; iK < m_kk; iK++) {
for (size_t i = 0; i < numBinaryInteractions_; i++) {
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
delAK = 0;
delBK = 0;
int delAK = 0;
int delBK = 0;
if (iA==iK) {
delAK = 1;
} else if (iB==iK) {
delBK = 1;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
g0 = -m_HE_b_ij[i] / RTT;
g1 = -m_HE_c_ij[i] / RTT;
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
double g0 = -m_HE_b_ij[i] / (GasConstant*T*T);
double g1 = -m_HE_c_ij[i] / (GasConstant*T*T);
double temp = (delAK * XB + XA * delBK - XA * XB) * (g0 + g1 * XB) + XA * XB * (delBK - XB) * g1;
dlnActCoeffdT_Scaled_[iK] += temp;
d2lnActCoeffdT2_Scaled_[iK] -= 2.0 * temp / T;
@ -528,18 +500,14 @@ void PhaseCombo_Interaction::getd2lnActCoeffdT2(doublereal* d2lnActCoeffdT2) con
void PhaseCombo_Interaction::getdlnActCoeffds(const doublereal dTds, const doublereal* const dXds,
doublereal* dlnActCoeffds) const
{
int delAK, delBK;
doublereal XA, XB, g0 , g1, dXA, dXB;
doublereal T = temperature();
doublereal RT = GasConstant*T;
doublereal xx;
s_update_dlnActCoeff_dT();
for (size_t iK = 0; iK < m_kk; iK++) {
/*
* We never sample the end of the mole fraction domains
*/
xx = std::max(moleFractions_[iK], SmallNumber);
double xx = std::max(moleFractions_[iK], SmallNumber);
/*
* First wipe out the ideal solution mixing term
*/
@ -551,8 +519,8 @@ void PhaseCombo_Interaction::getdlnActCoeffds(const doublereal dTds, const doub
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
delAK = 0;
delBK = 0;
int delAK = 0;
int delBK = 0;
if (iA==iK) {
delAK = 1;
@ -560,14 +528,14 @@ void PhaseCombo_Interaction::getdlnActCoeffds(const doublereal dTds, const doub
delBK = 1;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
dXA = dXds[iA];
dXB = dXds[iB];
double dXA = dXds[iA];
double dXB = dXds[iB];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
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]) / (GasConstant*T);
double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / (GasConstant*T);
dlnActCoeffds[iK] += ((delBK-XB)*dXA + (delAK-XA)*dXB)*(g0+2*g1*XB) + (delBK-XB)*2*g1*XA*dXB
+ dlnActCoeffdT_Scaled_[iK]*dTds;
@ -577,21 +545,16 @@ void PhaseCombo_Interaction::getdlnActCoeffds(const doublereal dTds, const doub
void PhaseCombo_Interaction::s_update_dlnActCoeff_dlnN_diag() const
{
int delAK, delBK;
doublereal XA, XB, XK, g0 , g1;
doublereal T = temperature();
doublereal RT = GasConstant*T;
doublereal xx;
dlnActCoeffdlnN_diag_.assign(m_kk, 0.0);
for (size_t iK = 0; iK < m_kk; iK++) {
XK = moleFractions_[iK];
double XK = moleFractions_[iK];
/*
* We never sample the end of the mole fraction domains
*/
xx = std::max(moleFractions_[iK], SmallNumber);
double xx = std::max(moleFractions_[iK], SmallNumber);
/*
* First wipe out the ideal solution mixing term
*/
@ -603,8 +566,8 @@ void PhaseCombo_Interaction::s_update_dlnActCoeff_dlnN_diag() const
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
delAK = 0;
delBK = 0;
int delAK = 0;
int delBK = 0;
if (iA==iK) {
delAK = 1;
@ -612,11 +575,11 @@ void PhaseCombo_Interaction::s_update_dlnActCoeff_dlnN_diag() const
delBK = 1;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
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]) / (GasConstant*T);
double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / (GasConstant*T);
dlnActCoeffdlnN_diag_[iK] += 2*(delBK-XB)*(g0*(delAK-XA)+g1*(2*(delAK-XA)*XB+XA*(delBK-XB)));
}
@ -627,14 +590,7 @@ void PhaseCombo_Interaction::s_update_dlnActCoeff_dlnN_diag() const
void PhaseCombo_Interaction::s_update_dlnActCoeff_dlnN() const
{
doublereal delAK, delBK;
double XA, XB, g0, g1, XM;
double xx , delKM;
double T = temperature();
double RT = GasConstant*T;
doublereal delAM, delBM;
dlnActCoeffdlnN_.zero();
/*
@ -644,13 +600,13 @@ void PhaseCombo_Interaction::s_update_dlnActCoeff_dlnN() const
/*
* We never sample the end of the mole fraction domains
*/
xx = std::max(moleFractions_[iK], SmallNumber);
double xx = std::max(moleFractions_[iK], SmallNumber);
for (size_t iM = 0; iM < m_kk; iM++) {
XM = moleFractions_[iM];
double XM = moleFractions_[iM];
if (xx > SmallNumber) {
delKM = 0.0;
double delKM = 0.0;
if (iK == iM) {
delKM = 1.0;
}
@ -662,10 +618,10 @@ void PhaseCombo_Interaction::s_update_dlnActCoeff_dlnN() const
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
delAK = 0.0;
delBK = 0.0;
delAM = 0.0;
delBM = 0.0;
double delAK = 0.0;
double delBK = 0.0;
double delAM = 0.0;
double delBM = 0.0;
if (iA==iK) {
delAK = 1.0;
} else if (iB==iK) {
@ -677,11 +633,11 @@ void PhaseCombo_Interaction::s_update_dlnActCoeff_dlnN() const
delBM = 1.0;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
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]) / (GasConstant*T);
double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / (GasConstant*T);
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);
@ -694,21 +650,17 @@ void PhaseCombo_Interaction::s_update_dlnActCoeff_dlnN() const
void PhaseCombo_Interaction::s_update_dlnActCoeff_dlnX_diag() const
{
doublereal XA, XB, g0 , g1;
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];
size_t iB = m_pSpecies_B_ij[i];
XA = moleFractions_[iA];
XB = moleFractions_[iB];
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
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]) / (GasConstant * T);
double g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / (GasConstant * T);
dlnActCoeffdlnX_diag_[iA] += XA*XB*(2*g1*-2*g0-6*g1*XB);
dlnActCoeffdlnX_diag_[iB] += XA*XB*(2*g1*-2*g0-6*g1*XB);
@ -770,8 +722,6 @@ void PhaseCombo_Interaction::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies",
"Incorrect name for processing this routine: " + xname);
}
string stemp;
size_t nParamsFound;
vector_fp vParams;
string iName = xmLBinarySpecies.attrib("speciesA");
if (iName == "") {
@ -807,11 +757,9 @@ void PhaseCombo_Interaction::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
m_pSpecies_A_ij[iSpot] = iSpecies;
m_pSpecies_B_ij[iSpot] = jSpecies;
size_t num = xmLBinarySpecies.nChildren();
for (size_t iChild = 0; iChild < num; iChild++) {
for (size_t iChild = 0; iChild < xmLBinarySpecies.nChildren(); iChild++) {
XML_Node& xmlChild = xmLBinarySpecies.child(iChild);
stemp = xmlChild.name();
string nodeName = lowercase(stemp);
string nodeName = lowercase(xmlChild.name());
/*
* Process the binary species interaction child elements
*/
@ -820,9 +768,7 @@ void PhaseCombo_Interaction::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
* Get the string containing all of the values
*/
ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessEnthalpy");
nParamsFound = vParams.size();
if (nParamsFound != 2) {
if (vParams.size() != 2) {
throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies::excessEnthalpy for " + ispName
+ "::" + jspName,
"wrong number of params found");
@ -836,9 +782,7 @@ void PhaseCombo_Interaction::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
* Get the string containing all of the values
*/
ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessEntropy");
nParamsFound = vParams.size();
if (nParamsFound != 2) {
if (vParams.size() != 2) {
throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies::excessEntropy for " + ispName
+ "::" + jspName,
"wrong number of params found");
@ -852,9 +796,7 @@ void PhaseCombo_Interaction::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
* Get the string containing all of the values
*/
ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Enthalpy");
nParamsFound = vParams.size();
if (nParamsFound != 2) {
if (vParams.size() != 2) {
throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies::excessVolume_Enthalpy for " + ispName
+ "::" + jspName,
"wrong number of params found");
@ -868,9 +810,7 @@ void PhaseCombo_Interaction::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
* Get the string containing all of the values
*/
ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Entropy");
nParamsFound = vParams.size();
if (nParamsFound != 2) {
if (vParams.size() != 2) {
throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies::excessVolume_Entropy for " + ispName
+ "::" + jspName,
"wrong number of params found");

42
src/thermo/PseudoBinaryVPSSTP.cpp
View file

@ -98,48 +98,45 @@ void PseudoBinaryVPSSTP::getElectrochemPotentials(doublereal* mu) const
void PseudoBinaryVPSSTP::calcPseudoBinaryMoleFractions() const
{
size_t k;
doublereal sumCat;
doublereal sumAnion;
doublereal sum = 0.0;
switch (PBType_) {
case PBTYPE_PASSTHROUGH:
for (k = 0; k < m_kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
PBMoleFractions_[k] = moleFractions_[k];
}
break;
case PBTYPE_SINGLEANION:
sumCat = 0.0;
sumAnion = 0.0;
for (k = 0; k < m_kk; k++) {
{
double sumCat = 0.0;
double sumAnion = 0.0;
for (size_t k = 0; k < m_kk; k++) {
moleFractionsTmp_[k] = moleFractions_[k];
}
for (k = 0; k < cationList_.size(); k++) {
for (size_t k = 0; k < cationList_.size(); k++) {
sumCat += moleFractions_[cationList_[k]];
}
sumAnion = moleFractions_[anionList_[k]];
sumAnion = moleFractions_[anionList_[0]];
PBMoleFractions_[0] = sumCat -sumAnion;
moleFractionsTmp_[indexSpecialSpecies_] -= PBMoleFractions_[0];
for (k = 0; k < numCationSpecies_; k++) {
for (size_t k = 0; k < numCationSpecies_; k++) {
PBMoleFractions_[1+k] = moleFractionsTmp_[cationList_[k]];
}
for (k = 0; k < numPassThroughSpecies_; k++) {
for (size_t k = 0; k < numPassThroughSpecies_; k++) {
PBMoleFractions_[neutralPBindexStart + k] =
moleFractions_[cationList_[k]];
}
sum = std::max(0.0, PBMoleFractions_[0]);
for (k = 1; k < numPBSpecies_; k++) {
double sum = std::max(0.0, PBMoleFractions_[0]);
for (size_t k = 1; k < numPBSpecies_; k++) {
sum += PBMoleFractions_[k];
}
for (k = 0; k < numPBSpecies_; k++) {
for (size_t k = 0; k < numPBSpecies_; k++) {
PBMoleFractions_[k] /= sum;
}
break;
}
case PBTYPE_SINGLECATION:
throw CanteraError("eosType", "Unknown type");
@ -195,13 +192,12 @@ std::string PseudoBinaryVPSSTP::report(bool show_thermo, doublereal threshold) c
sprintf(p, " potential %12.6g V\n", phi);
s += p;
size_t kk = nSpecies();
vector_fp x(kk);
vector_fp molal(kk);
vector_fp mu(kk);
vector_fp muss(kk);
vector_fp acMolal(kk);
vector_fp actMolal(kk);
vector_fp x(m_kk);
vector_fp molal(m_kk);
vector_fp mu(m_kk);
vector_fp muss(m_kk);
vector_fp acMolal(m_kk);
vector_fp actMolal(m_kk);
getMoleFractions(&x[0]);
getChemPotentials(&mu[0]);

18
src/thermo/PureFluidPhase.cpp
View file

@ -70,8 +70,8 @@ void PureFluidPhase::initThermo()
setMolecularWeight(0,m_mw);
double one = 1.0;
setMoleFractions(&one);
double cp0_R, h0_RT, s0_R, T0, p;
T0 = 298.15;
double cp0_R, h0_RT, s0_R, p;
double T0 = 298.15;
if (T0 < m_sub->Tcrit()) {
m_sub->Set(tpx::PropertyPair::TX, T0, 1.0);
p = 0.01*m_sub->P();
@ -218,22 +218,17 @@ void PureFluidPhase::getStandardChemPotentials(doublereal* mu) const
void PureFluidPhase::getEnthalpy_RT(doublereal* hrt) const
{
doublereal rt = _RT();
doublereal h = enthalpy_mole();
hrt[0] = h / rt;
hrt[0] = enthalpy_mole() / _RT();
}
void PureFluidPhase::getEntropy_R(doublereal* sr) const
{
doublereal s = entropy_mole();
sr[0] = s / GasConstant;
sr[0] = entropy_mole() / GasConstant;
}
void PureFluidPhase::getGibbs_RT(doublereal* grt) const
{
doublereal rt = _RT();
doublereal g = gibbs_mole();
grt[0] = g / rt;
grt[0] = gibbs_mole() / _RT();
}
void PureFluidPhase::getEnthalpy_RT_ref(doublereal* hrt) const
@ -327,8 +322,7 @@ void PureFluidPhase::setState_SP(doublereal s, doublereal p,
doublereal PureFluidPhase::satPressure(doublereal t)
{
doublereal vsv = m_sub->v();
Set(tpx::PropertyPair::TV,t,vsv);
Set(tpx::PropertyPair::TV, t, m_sub->v());
return m_sub->Ps();
}

105
src/thermo/RedlichKisterVPSSTP.cpp
View file

@ -154,7 +154,6 @@ void RedlichKisterVPSSTP::getElectrochemPotentials(doublereal* mu) const
void RedlichKisterVPSSTP::getChemPotentials(doublereal* mu) const
{
doublereal xx;
/*
* First get the standard chemical potentials in
* molar form.
@ -169,18 +168,17 @@ void RedlichKisterVPSSTP::getChemPotentials(doublereal* mu) const
doublereal RT = GasConstant * temperature();
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(moleFractions_[k], SmallNumber);
double xx = std::max(moleFractions_[k], SmallNumber);
mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]);
}
}
doublereal RedlichKisterVPSSTP::enthalpy_mole() const
{
size_t kk = nSpecies();
double h = 0;
vector_fp hbar(kk);
vector_fp hbar(m_kk);
getPartialMolarEnthalpies(&hbar[0]);
for (size_t i = 0; i < kk; i++) {
for (size_t i = 0; i < m_kk; i++) {
h += moleFractions_[i]*hbar[i];
}
return h;
@ -188,11 +186,10 @@ doublereal RedlichKisterVPSSTP::enthalpy_mole() const
doublereal RedlichKisterVPSSTP::entropy_mole() const
{
size_t kk = nSpecies();
double s = 0;
vector_fp sbar(kk);
vector_fp sbar(m_kk);
getPartialMolarEntropies(&sbar[0]);
for (size_t i = 0; i < kk; i++) {
for (size_t i = 0; i < m_kk; i++) {
s += moleFractions_[i]*sbar[i];
}
return s;
@ -200,11 +197,10 @@ doublereal RedlichKisterVPSSTP::entropy_mole() const
doublereal RedlichKisterVPSSTP::cp_mole() const
{
size_t kk = nSpecies();
double cp = 0;
vector_fp cpbar(kk);
vector_fp cpbar(m_kk);
getPartialMolarCp(&cpbar[0]);
for (size_t i = 0; i < kk; i++) {
for (size_t i = 0; i < m_kk; i++) {
cp += moleFractions_[i]*cpbar[i];
}
return cp;
@ -225,9 +221,8 @@ void RedlichKisterVPSSTP::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] *= GasConstant * T;
}
/*
* Update the activity coefficients, This also update the
@ -235,9 +230,8 @@ void RedlichKisterVPSSTP::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] -= GasConstant * T * T * dlnActCoeffdT_Scaled_[k];
}
}
@ -268,7 +262,6 @@ void RedlichKisterVPSSTP::getPartialMolarCp(doublereal* cpbar) const
void RedlichKisterVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
{
double xx;
/*
* Get the nondimensional standard state entropies
*/
@ -282,7 +275,7 @@ void RedlichKisterVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
s_update_dlnActCoeff_dT();
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(moleFractions_[k], SmallNumber);
double xx = std::max(moleFractions_[k], SmallNumber);
sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
}
/*
@ -319,14 +312,9 @@ void RedlichKisterVPSSTP::initLengths()
void RedlichKisterVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
std::string subname = "RedlichKisterVPSSTP::initThermoXML";
std::string stemp;
if ((int) id_.size() > 0) {
string idp = phaseNode.id();
if (idp != id_) {
throw CanteraError(subname,
"phasenode and Id are incompatible");
}
if ((int) id_.size() > 0 && phaseNode.id() != id_) {
throw CanteraError("RedlichKisterVPSSTP::initThermoXML",
"phasenode and Id are incompatible");
}
/*
@ -334,14 +322,14 @@ void RedlichKisterVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string&
* <thermo model="Redlich-Kister" />
*/
if (!phaseNode.hasChild("thermo")) {
throw CanteraError(subname, "no thermo XML node");
throw CanteraError("RedlichKisterVPSSTP::initThermoXML",
"no thermo XML node");
}
XML_Node& thermoNode = phaseNode.child("thermo");
std::string mStringa = thermoNode.attrib("model");
std::string mString = lowercase(mStringa);
if (mString != "redlich-kister") {
throw CanteraError(subname.c_str(),
"Unknown thermo model: " + mStringa + " - This object only knows \"Redlich-Kister\" ");
std::string mString = thermoNode.attrib("model");
if (lowercase(mString) != "redlich-kister") {
throw CanteraError("RedlichKisterVPSSTP::initThermoXML",
"Unknown thermo model: " + mString + " - This object only knows \"Redlich-Kister\" ");
}
/*
@ -351,24 +339,19 @@ void RedlichKisterVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string&
XML_Node* acNodePtr = 0;
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& acNode = thermoNode.child("activityCoefficients");
acNodePtr = &acNode;
mStringa = acNode.attrib("model");
mString = lowercase(mStringa);
if (mString != "redlich-kister") {
throw CanteraError(subname.c_str(),
"Unknown activity coefficient model: " + mStringa);
mString = acNode.attrib("model");
if (lowercase(mString) != "redlich-kister") {
throw CanteraError("RedlichKisterVPSSTP::initThermoXML",
"Unknown activity coefficient model: " + mString);
}
size_t n = acNodePtr->nChildren();
for (size_t i = 0; i < n; i++) {
XML_Node& xmlACChild = acNodePtr->child(i);
stemp = xmlACChild.name();
std::string nodeName = lowercase(stemp);
for (size_t i = 0; i < acNode.nChildren(); i++) {
XML_Node& xmlACChild = acNode.child(i);
/*
* Process a binary salt field, or any of the other XML fields
* that make up the Pitzer Database. Entries will be ignored
* if any of the species in the entry isn't in the solution.
*/
if (nodeName == "binaryneutralspeciesparameters") {
if (lowercase(xmlACChild.name()) == "binaryneutralspeciesparameters") {
readXMLBinarySpecies(xmlACChild);
}
}
@ -381,10 +364,7 @@ void RedlichKisterVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string&
void RedlichKisterVPSSTP::s_update_lnActCoeff() const
{
doublereal XA, XB;
doublereal T = temperature();
doublereal RT = GasConstant * T;
lnActCoeff_Scaled_.assign(m_kk, 0.0);
/*
@ -396,8 +376,8 @@ void RedlichKisterVPSSTP::s_update_lnActCoeff() const
for (size_t i = 0; i < numBinaryInteractions_; i++) {
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
XA = moleFractions_[iA];
XB = moleFractions_[iB];
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
doublereal deltaX = XA - XB;
size_t N = m_N_ij[i];
vector_fp& he_vec = m_HE_m_ij[i];
@ -408,7 +388,7 @@ void RedlichKisterVPSSTP::s_update_lnActCoeff() const
doublereal sumMm1 = 0.0;
doublereal sum2 = 0.0;
for (size_t m = 0; m < N; m++) {
doublereal A_ge = (he_vec[m] - T * se_vec[m]) / RT;
doublereal A_ge = (he_vec[m] - T * se_vec[m]) / (GasConstant * T);
sum += A_ge * poly;
sum2 += A_ge * (m + 1) * poly;
poly *= deltaX;
@ -435,7 +415,7 @@ void RedlichKisterVPSSTP::s_update_lnActCoeff() const
double polyk = 1.0;
double fac = 2.0 * XA - 1.0;
for (int m = 0; m < N; m++) {
doublereal A_ge = (he_vec[m] - T * se_vec[m]) / RT;
doublereal A_ge = (he_vec[m] - T * se_vec[m]) / (GasConstant * T);
lnA += A_ge * oneMXA * oneMXA * polyk * (1.0 + 2.0 * XA * m / fac);
lnB += A_ge * XA * XA * polyk * (1.0 - 2.0 * oneMXA * m / fac);
polyk *= fac;
@ -449,15 +429,14 @@ void RedlichKisterVPSSTP::s_update_lnActCoeff() const
void RedlichKisterVPSSTP::s_update_dlnActCoeff_dT() const
{
doublereal XA, XB;
dlnActCoeffdT_Scaled_.assign(m_kk, 0.0);
d2lnActCoeffdT2_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];
XA = moleFractions_[iA];
XB = moleFractions_[iB];
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
doublereal deltaX = XA - XB;
size_t N = m_N_ij[i];
doublereal poly = 1.0;
@ -509,7 +488,6 @@ void RedlichKisterVPSSTP::getd2lnActCoeffdT2(doublereal* d2lnActCoeffdT2) const
void RedlichKisterVPSSTP::s_update_dlnActCoeff_dX_() const
{
doublereal XA, XB;
doublereal T = temperature();
dlnActCoeff_dX_.zero();
@ -517,8 +495,8 @@ void RedlichKisterVPSSTP::s_update_dlnActCoeff_dX_() const
for (size_t i = 0; i < numBinaryInteractions_; i++) {
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
XA = moleFractions_[iA];
XB = moleFractions_[iB];
double XA = moleFractions_[iA];
double XB = moleFractions_[iB];
doublereal deltaX = XA - XB;
size_t N = m_N_ij[i];
doublereal poly = 1.0;
@ -639,9 +617,8 @@ void RedlichKisterVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
throw CanteraError("RedlichKisterVPSSTP::readXMLBinarySpecies",
"Incorrect name for processing this routine: " + xname);
}
std::string stemp;
size_t Npoly = 0;
vector_fp hParams, sParams, vParams;
vector_fp hParams, sParams;
std::string iName = xmLBinarySpecies.attrib("speciesA");
if (iName == "") {
throw CanteraError("RedlichKisterVPSSTP::readXMLBinarySpecies", "no speciesA attrib");
@ -681,11 +658,9 @@ void RedlichKisterVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
m_pSpecies_A_ij[iSpot] = iSpecies;
m_pSpecies_B_ij[iSpot] = jSpecies;
size_t num = xmLBinarySpecies.nChildren();
for (size_t iChild = 0; iChild < num; iChild++) {
for (size_t iChild = 0; iChild < xmLBinarySpecies.nChildren(); iChild++) {
XML_Node& xmlChild = xmLBinarySpecies.child(iChild);
stemp = xmlChild.name();
string nodeName = lowercase(stemp);
string nodeName = lowercase(xmlChild.name());
/*
* Process the binary species interaction child elements
*/
@ -716,16 +691,14 @@ void RedlichKisterVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
#ifdef DEBUG_MODE
void RedlichKisterVPSSTP::Vint(double& VintOut, double& voltsOut)
{
doublereal XA;
doublereal T = temperature();
doublereal RT = GasConstant * T;
double Volts = 0.0;
lnActCoeff_Scaled_.assign(m_kk, 0.0);
for (size_t i = 0; i < numBinaryInteractions_; i++) {
size_t iA = m_pSpecies_A_ij[i];
XA = moleFractions_[iA];
double XA = moleFractions_[iA];
if (XA <= 1.0E-14) {
XA = 1.0E-14;
}
@ -752,7 +725,7 @@ void RedlichKisterVPSSTP::Vint(double& VintOut, double& voltsOut)
}
Volts /= Faraday;
double termp = RT * log((1.0 - XA)/XA) / Faraday;
double termp = GasConstant * T * log((1.0 - XA)/XA) / Faraday;
VintOut = Volts;
voltsOut = Volts + termp;

13
src/thermo/RedlichKwongMFTP.cpp
View file

@ -179,8 +179,7 @@ int RedlichKwongMFTP::eosType() const
doublereal RedlichKwongMFTP::enthalpy_mole() const
{
_updateReferenceStateThermo();
doublereal rt = _RT();
doublereal h_ideal = rt * mean_X(DATA_PTR(m_h0_RT));
doublereal h_ideal = _RT() * mean_X(DATA_PTR(m_h0_RT));
doublereal h_nonideal = hresid();
return h_ideal + h_nonideal;
}
@ -222,10 +221,8 @@ doublereal RedlichKwongMFTP::pressure() const
_updateReferenceStateThermo();
// Get a copy of the private variables stored in the State object
double rho = density();
doublereal T = temperature();
doublereal mmw = meanMolecularWeight();
double molarV = mmw / rho;
double molarV = meanMolecularWeight() / density();
double pp = GasConstant * T/(molarV - m_b_current) - m_a_current/(sqrt(T) * molarV * (molarV + m_b_current));
@ -249,8 +246,7 @@ void RedlichKwongMFTP::calcDensity()
* Set the density in the parent State object directly,
* by calling the Phase::setDensity() function.
*/
double dens = 1.0/invDens;
Phase::setDensity(dens);
Phase::setDensity(1.0/invDens);
}
void RedlichKwongMFTP::setTemperature(const doublereal temp)
@ -379,10 +375,9 @@ void RedlichKwongMFTP::getChemPotentials_RT(doublereal* muRT) const
void RedlichKwongMFTP::getChemPotentials(doublereal* mu) const
{
getGibbs_ref(mu);
doublereal xx;
doublereal rt = temperature() * GasConstant;
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(SmallNumber, moleFraction(k));
double xx = std::max(SmallNumber, moleFraction(k));
mu[k] += rt*(log(xx));
}

9
src/thermo/SemiconductorPhase.cpp
View file

@ -5,13 +5,9 @@ using namespace std;
namespace Cantera
{
const doublereal JD_const1 = 1.0/sqrt(8.0);
const doublereal JD_const2 = 3.0/16.0 - sqrt(3.0)/9.0;
static doublereal JoyceDixon(doublereal r)
{
return log(r) + JD_const1*r - JD_const2*r*r;
return log(r) + 1.0/sqrt(8.0)*r - (3.0/16.0 - sqrt(3.0)/9.0)*r*r;
}
@ -21,8 +17,7 @@ SemiconductorPhase::SemiconductorPhase(std::string infile,
void SemiconductorPhase::getChemPotentials(doublereal* mu) const
{
getActivityConcentrations(DATA_PTR(m_work));
doublereal r = m_work[0]/nc();
mu[0] = ec() + GasConstant*temperature()*(JoyceDixon(r));
mu[0] = ec() + GasConstant*temperature()*(JoyceDixon(m_work[0]/nc()));
mu[1] = ev() + GasConstant*temperature()*(log(m_work[1]/nv()));
}

48
src/thermo/ShomateThermo.h
View file

@ -120,8 +120,7 @@ public:
}
virtual SpeciesThermo* duplMyselfAsSpeciesThermo() const {
ShomateThermo* st = new ShomateThermo(*this);
return (SpeciesThermo*) st;
return new ShomateThermo(*this);
}
//! Install a new species thermodynamic property
@ -233,8 +232,7 @@ public:
const std::vector<ShomatePoly> &mlg = m_low[grp-1];
const ShomatePoly* nlow = &(mlg[pos]);
doublereal tmid = nlow->maxTemp();
if (t < tmid) {
if (t < nlow->maxTemp()) {
nlow->updateProperties(&m_t[0], cp_R, h_RT, s_R);
} else {
const std::vector<ShomatePoly> &mhg = m_high[grp-1];
@ -245,8 +243,6 @@ public:
virtual void update(doublereal t, doublereal* cp_R,
doublereal* h_RT, doublereal* s_R) const {
int i;
doublereal tt = 1.e-3*t;
m_t[0] = tt;
m_t[1] = tt*tt;
@ -257,7 +253,7 @@ public:
m_t[6] = 1.0/(GasConstant * t);
std::vector<ShomatePoly>::const_iterator _begin, _end;
for (i = 0; i != m_ngroups; i++) {
for (int i = 0; i != m_ngroups; i++) {
if (t > m_tmid[i]) {
_begin = m_high[i].begin();
_end = m_high[i].end();
@ -339,47 +335,37 @@ public:
}
virtual doublereal reportOneHf298(const size_t k) const {
doublereal h;
doublereal t = 298.15;
size_t grp = getValue(m_group_map, k);
size_t pos = getValue(m_posInGroup_map, k);
const std::vector<ShomatePoly> &mlg = m_low[grp-1];
const ShomatePoly* nlow = &(mlg[pos]);
const ShomatePoly& nlow = m_low[grp-1][pos];
doublereal tmid = nlow->maxTemp();
if (t <= tmid) {
h = nlow->reportHf298();
if (nlow.maxTemp() > 298.15) {
return nlow.reportHf298();
} else {
const std::vector<ShomatePoly> &mhg = m_high[grp-1];
const ShomatePoly* nhigh = &(mhg[pos]);
h = nhigh->reportHf298();
const ShomatePoly& nhigh = m_high[grp-1][pos];
return nhigh.reportHf298();
}
return h;
}
virtual void modifyOneHf298(const size_t k, const doublereal Hf298New) {
size_t grp = m_group_map[k];
size_t pos = m_posInGroup_map[k];
std::vector<ShomatePoly> &mlg = m_low[grp-1];
ShomatePoly* nlow = &(mlg[pos]);
std::vector<ShomatePoly> &mhg = m_high[grp-1];
ShomatePoly* nhigh = &(mhg[pos]);
doublereal tmid = nlow->maxTemp();
ShomatePoly& nlow = m_low[grp-1][pos];
ShomatePoly& nhigh = m_high[grp-1][pos];
double hnow = reportOneHf298(k);
double delH = Hf298New - hnow;
if (298.15 <= tmid) {
nlow->modifyOneHf298(k, Hf298New);
double h = nhigh->reportHf298(0);
if (nlow.maxTemp() > 298.15) {
nlow.modifyOneHf298(k, Hf298New);
double h = nhigh.reportHf298(0);
double hnew = h + delH;
nhigh->modifyOneHf298(k, hnew);
nhigh.modifyOneHf298(k, hnew);
} else {
nhigh->modifyOneHf298(k, Hf298New);
double h = nlow->reportHf298(0);
nhigh.modifyOneHf298(k, Hf298New);
double h = nlow.reportHf298(0);
double hnew = h + delH;
nlow->modifyOneHf298(k, hnew);
nlow.modifyOneHf298(k, hnew);
}
}

27
src/thermo/SingleSpeciesTP.cpp
View file

@ -118,8 +118,7 @@ doublereal SingleSpeciesTP::cv_mole() const
doublereal cvbar = cp_mole();
doublereal alpha = thermalExpansionCoeff();
doublereal beta = isothermalCompressibility();
doublereal molecW = molecularWeight(0);
doublereal V = molecW/density();
doublereal V = molecularWeight(0)/density();
doublereal T = temperature();
if (beta != 0.0) {
cvbar -= alpha * alpha * V * T / beta;
@ -139,8 +138,7 @@ void SingleSpeciesTP::getChemPotentials(doublereal* mu) const
void SingleSpeciesTP::getChemPotentials_RT(doublereal* murt) const
{
getStandardChemPotentials(murt);
double rt = GasConstant * temperature();
murt[0] /= rt;
murt[0] /= GasConstant * temperature();
}
void SingleSpeciesTP::getElectrochemPotentials(doublereal* mu) const
@ -150,16 +148,14 @@ void SingleSpeciesTP::getElectrochemPotentials(doublereal* mu) const
void SingleSpeciesTP::getPartialMolarEnthalpies(doublereal* hbar) const
{
double _rt = GasConstant * temperature();
getEnthalpy_RT(hbar);
hbar[0] *= _rt;
hbar[0] *= GasConstant * temperature();
}
void SingleSpeciesTP::getPartialMolarIntEnergies(doublereal* ubar) const
{
double _rt = GasConstant * temperature();
getIntEnergy_RT(ubar);
ubar[0] *= _rt;
ubar[0] *= GasConstant * temperature();
}
void SingleSpeciesTP::getPartialMolarEntropies(doublereal* sbar) const
@ -176,9 +172,7 @@ void SingleSpeciesTP::getPartialMolarCp(doublereal* cpbar) const
void SingleSpeciesTP::getPartialMolarVolumes(doublereal* vbar) const
{
double mw = molecularWeight(0);
double dens = density();
vbar[0] = mw / dens;
vbar[0] = molecularWeight(0) / density();
}
/*
@ -193,9 +187,7 @@ void SingleSpeciesTP::getPureGibbs(doublereal* gpure) const
void SingleSpeciesTP::getStandardVolumes(doublereal* vbar) const
{
double mw = molecularWeight(0);
double dens = density();
vbar[0] = mw / dens;
vbar[0] = molecularWeight(0) / density();
}
/*
@ -319,10 +311,9 @@ void SingleSpeciesTP::initThermo()
/*
* Resize temporary arrays.
*/
int leng = 1;
m_h0_RT.resize(leng);
m_cp0_R.resize(leng);
m_s0_R.resize(leng);
m_h0_RT.resize(1);
m_cp0_R.resize(1);
m_s0_R.resize(1);
/*
* Make sure the species mole fraction is equal to 1.0;

50
src/thermo/SpeciesThermoFactory.cpp
View file

@ -55,12 +55,10 @@ static void getSpeciesThermoTypes(std::vector<XML_Node*> & spDataNodeList,
int& has_nasa, int& has_shomate, int& has_simple,
int& has_other)
{
size_t ns = spDataNodeList.size();
for (size_t n = 0; n < ns; n++) {
for (size_t n = 0; n < spDataNodeList.size(); n++) {
XML_Node* spNode = spDataNodeList[n];
if (spNode->hasChild("standardState")) {
const XML_Node& ss = spNode->child("standardState");
string mname = ss["model"];
string mname = spNode->child("standardState")["model"];
if (mname == "water" || mname == "waterIAPWS") {
has_other = 1;
continue;
@ -254,13 +252,11 @@ SpeciesThermoInterpType* newSpeciesThermoInterpType(const std::string& stype,
static SpeciesThermoInterpType* newNasaThermoFromXML(
const std::string& speciesName, vector<XML_Node*> nodes)
{
doublereal tmin0, tmax0, tmin1, tmax1, tmin, tmid, tmax;
const XML_Node& f0 = *nodes[0];
bool dualRange = (nodes.size() > 1);
tmin0 = fpValue(f0["Tmin"]);
tmax0 = fpValue(f0["Tmax"]);
double tmin0 = fpValue(f0["Tmin"]);
double tmax0 = fpValue(f0["Tmax"]);
doublereal p0 = OneAtm;
if (f0.hasAttrib("P0")) {
@ -271,14 +267,15 @@ static SpeciesThermoInterpType* newNasaThermoFromXML(
}
p0 = OneAtm;
tmin1 = tmax0;
tmax1 = tmin1 + 0.0001;
double tmin1 = tmax0;
double tmax1 = tmin1 + 0.0001;
if (dualRange) {
tmin1 = fpValue(nodes[1]->attrib("Tmin"));
tmax1 = fpValue(nodes[1]->attrib("Tmax"));
}
vector_fp c0, c1;
doublereal tmin, tmid, tmax;
if (fabs(tmax0 - tmin1) < 0.01) {
// f0 has the lower T data, and f1 the higher T data
tmin = tmin0;
@ -321,12 +318,9 @@ static SpeciesThermoInterpType* newNasaThermoFromXML(
SpeciesThermoInterpType* newShomateForMineralEQ3(const std::string& name,
const XML_Node& MinEQ3node)
{
std::string astring = MinEQ3node["Tmin"];
doublereal tmin0 = strSItoDbl(astring);
astring = MinEQ3node["Tmax"];
doublereal tmax0 = strSItoDbl(astring);
astring = MinEQ3node["Pref"];
doublereal p0 = strSItoDbl(astring);
doublereal tmin0 = strSItoDbl(MinEQ3node["Tmin"]);
doublereal tmax0 = strSItoDbl(MinEQ3node["Tmax"]);
doublereal p0 = strSItoDbl(MinEQ3node["Pref"]);
doublereal deltaG_formation_pr_tr =
getFloatDefaultUnits(MinEQ3node, "DG0_f_Pr_Tr", "cal/gmol", "actEnergy");
@ -370,13 +364,6 @@ SpeciesThermoInterpType* newShomateForMineralEQ3(const std::string& name,
double Gs = ScalcS - S298smGs;
double c0[7] = {As, Bs, Cs, Ds, Es, Fs, Gs};
c0[0] = As;
c0[1] = Bs;
c0[2] = Cs;
c0[3] = Ds;
c0[4] = Es;
c0[5] = Fs;
c0[6] = Gs;
return newSpeciesThermoInterpType(SHOMATE1, tmin0, tmax0, p0, c0);
}
@ -393,13 +380,12 @@ SpeciesThermoInterpType* newShomateForMineralEQ3(const std::string& name,
static SpeciesThermoInterpType* newShomateThermoFromXML(
const std::string& speciesName, vector<XML_Node*>& nodes)
{
doublereal tmin0, tmax0, tmin1, tmax1, tmin, tmid, tmax;
bool dualRange = false;
if (nodes.size() == 2) {
dualRange = true;
}
tmin0 = fpValue(nodes[0]->attrib("Tmin"));
tmax0 = fpValue(nodes[0]->attrib("Tmax"));
double tmin0 = fpValue(nodes[0]->attrib("Tmin"));
double tmax0 = fpValue(nodes[0]->attrib("Tmax"));
doublereal p0 = OneAtm;
if (nodes[0]->hasAttrib("P0")) {
@ -410,14 +396,15 @@ static SpeciesThermoInterpType* newShomateThermoFromXML(
}
p0 = OneAtm;
tmin1 = tmax0;
tmax1 = tmin1 + 0.0001;
double tmin1 = tmax0;
double tmax1 = tmin1 + 0.0001;
if (dualRange) {
tmin1 = fpValue(nodes[1]->attrib("Tmin"));
tmax1 = fpValue(nodes[1]->attrib("Tmax"));
}
vector_fp c0, c1;
doublereal tmin, tmid, tmax;
if (fabs(tmax0 - tmin1) < 0.01) {
tmin = tmin0;
tmid = tmax0;
@ -581,12 +568,11 @@ static SpeciesThermoInterpType* newAdsorbateThermoFromXML(
if (f.hasChild("floatArray")) {
getFloatArray(f.child("floatArray"), freqs, false);
}
size_t nfreq = freqs.size();
for (size_t n = 0; n < nfreq; n++) {
for (size_t n = 0; n < freqs.size(); n++) {
freqs[n] *= 3.0e10;
}
vector_fp coeffs(nfreq + 2);
coeffs[0] = static_cast<double>(nfreq);
vector_fp coeffs(freqs.size() + 2);
coeffs[0] = static_cast<double>(freqs.size());
coeffs[1] = getFloat(f, "binding_energy", "toSI");
copy(freqs.begin(), freqs.end(), coeffs.begin() + 2);
return new Adsorbate(0, tmin, tmax, pref, &coeffs[0]);

53
src/thermo/StoichSubstanceSSTP.cpp
View file

@ -40,8 +40,7 @@ StoichSubstanceSSTP::StoichSubstanceSSTP(const std::string& infile, std::string
"Couldn't find phase name in file:" + id_);
}
// Check the model name to ensure we have compatibility
const XML_Node& th = xphase->child("thermo");
std::string model = th["model"];
std::string model = xphase->child("thermo")["model"];
if (model != "StoichSubstance" && model != "StoichSubstanceSSTP") {
throw CanteraError("StoichSubstanceSSTP::StoichSubstanceSSTP",
"thermo model attribute must be StoichSubstance");
@ -51,15 +50,11 @@ StoichSubstanceSSTP::StoichSubstanceSSTP(const std::string& infile, std::string
StoichSubstanceSSTP::StoichSubstanceSSTP(XML_Node& xmlphase, const std::string& id_)
{
if (id_ != "") {
std::string idxml = xmlphase["id"];
if (id_ != idxml) {
throw CanteraError("StoichSubstanceSSTP::StoichSubstanceSSTP",
"id's don't match");
}
if (id_ != "" && id_ != xmlphase["id"]) {
throw CanteraError("StoichSubstanceSSTP::StoichSubstanceSSTP",
"id's don't match");
}
const XML_Node& th = xmlphase.child("thermo");
std::string model = th["model"];
std::string model = xmlphase.child("thermo")["model"];
if (model != "StoichSubstance" && model != "StoichSubstanceSSTP") {
throw CanteraError("StoichSubstanceSSTP::StoichSubstanceSSTP",
"thermo model attribute must be StoichSubstance");
@ -158,9 +153,8 @@ void StoichSubstanceSSTP::getStandardChemPotentials(doublereal* mu0) const
void StoichSubstanceSSTP::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 / (GasConstant * temperature());
}
void StoichSubstanceSSTP::getEntropy_R(doublereal* sr) const
@ -183,8 +177,7 @@ void StoichSubstanceSSTP::getCp_R(doublereal* cpr) const
void StoichSubstanceSSTP::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() / (GasConstant * temperature());
}
/*
@ -194,8 +187,7 @@ void StoichSubstanceSSTP::getIntEnergy_RT(doublereal* urt) const
void StoichSubstanceSSTP::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() / (GasConstant * temperature());
}
/*
@ -247,15 +239,13 @@ void StoichSubstanceSSTP::initThermoXML(XML_Node& phaseNode, const std::string&
void StoichSubstanceSSTP::setParameters(int n, doublereal* const c)
{
doublereal rho = c[0];
setDensity(rho);
setDensity(c[0]);
}
void StoichSubstanceSSTP::getParameters(int& n, doublereal* const c) const
{
doublereal rho = density();
n = 1;
c[0] = rho;
c[0] = density();
}
void StoichSubstanceSSTP::setParametersFromXML(const XML_Node& eosdata)
@ -265,8 +255,7 @@ void StoichSubstanceSSTP::setParametersFromXML(const XML_Node& eosdata)
throw CanteraError("StoichSubstanceSSTP::setParametersFromXML",
"thermo model attribute must be StoichSubstance");
}
doublereal rho = ctml::getFloat(eosdata, "density", "toSI");
setDensity(rho);
setDensity(ctml::getFloat(eosdata, "density", "toSI"));
}
// ------ Methods of class electrodeElectron ------
@ -289,9 +278,7 @@ electrodeElectron::electrodeElectron(const std::string& infile, std::string id_)
"Couldn't find phase name in file:" + id_);
}
// Check the model name to ensure we have compatibility
const XML_Node& th = xphase->child("thermo");
std::string model = th["model"];
if (model != "electrodeElectron") {
if (xphase->child("thermo")["model"] != "electrodeElectron") {
throw CanteraError("electrodeElectron::electrodeElectron",
"thermo model attribute must be electrodeElectron");
}
@ -301,16 +288,11 @@ electrodeElectron::electrodeElectron(const std::string& infile, std::string id_)
electrodeElectron::electrodeElectron(XML_Node& xmlphase, const std::string& id_) :
StoichSubstanceSSTP()
{
if (id_ != "") {
std::string idxml = xmlphase["id"];
if (id_ != idxml) {
throw CanteraError("electrodeElectron::electrodeElectron",
"id's don't match");
}
if (id_ != "" && id_ != xmlphase["id"]) {
throw CanteraError("electrodeElectron::electrodeElectron",
"id's don't match");
}
const XML_Node& th = xmlphase.child("thermo");
std::string model = th["model"];
if (model != "electrodeElectron") {
if (xmlphase.child("thermo")["model"] != "electrodeElectron") {
throw CanteraError("electrodeElectron::electrodeElectron",
"thermo model attribute must be electrodeElectron");
}
@ -338,8 +320,7 @@ electrodeElectron::~electrodeElectron()
void electrodeElectron::setParametersFromXML(const XML_Node& eosdata)
{
std::string model = eosdata["model"];
if (model != "electrodeElectron") {
if (eosdata["model"] != "electrodeElectron") {
throw CanteraError("electrodeElectron::setParametersFromXML",
"thermo model attribute must be electrodeElectron");
}

29
src/thermo/SurfPhase.cpp
View file

@ -45,8 +45,7 @@ SurfPhase::SurfPhase(const std::string& infile, std::string id_) :
"Couldn't find phase name in file:" + id_);
}
// Check the model name to ensure we have compatibility
const XML_Node& th = xphase->child("thermo");
string model = th["model"];
string model = xphase->child("thermo")["model"];
if (model != "Surface" && model != "Edge") {
throw CanteraError("SurfPhase::SurfPhase",
"thermo model attribute must be Surface or Edge");
@ -59,8 +58,7 @@ SurfPhase::SurfPhase(XML_Node& xmlphase) :
m_logn0(0.0),
m_press(OneAtm)
{
const XML_Node& th = xmlphase.child("thermo");
string model = th["model"];
string model = xmlphase.child("thermo")["model"];
if (model != "Surface" && model != "Edge") {
throw CanteraError("SurfPhase::SurfPhase",
"thermo model attribute must be Surface or Edge");
@ -216,29 +214,25 @@ void SurfPhase::getPureGibbs(doublereal* g) const
void SurfPhase::getGibbs_RT(doublereal* grt) const
{
_updateThermo();
double rrt = 1.0/(GasConstant*temperature());
scale(m_mu0.begin(), m_mu0.end(), grt, rrt);
scale(m_mu0.begin(), m_mu0.end(), grt, 1.0/(GasConstant*temperature()));
}
void SurfPhase::getEnthalpy_RT(doublereal* hrt) const
{
_updateThermo();
double rrt = 1.0/(GasConstant*temperature());
scale(m_h0.begin(), m_h0.end(), hrt, rrt);
scale(m_h0.begin(), m_h0.end(), hrt, 1.0/(GasConstant*temperature()));
}
void SurfPhase::getEntropy_R(doublereal* sr) const
{
_updateThermo();
double rr = 1.0/GasConstant;
scale(m_s0.begin(), m_s0.end(), sr, rr);
scale(m_s0.begin(), m_s0.end(), sr, 1.0/GasConstant);
}
void SurfPhase::getCp_R(doublereal* cpr) const
{
_updateThermo();
double rr = 1.0/GasConstant;
scale(m_cp0.begin(), m_cp0.end(), cpr, rr);
scale(m_cp0.begin(), m_cp0.end(), cpr, 1.0/GasConstant);
}
void SurfPhase::getStandardVolumes(doublereal* vol) const
@ -341,13 +335,11 @@ void SurfPhase::getCoverages(doublereal* theta) const
void SurfPhase::setCoveragesByName(const std::string& cov)
{
size_t kk = nSpecies();
compositionMap cc = parseCompString(cov, speciesNames());
doublereal c;
vector_fp cv(kk, 0.0);
vector_fp cv(m_kk, 0.0);
bool ifound = false;
for (size_t k = 0; k < kk; k++) {
c = cc[speciesName(k)];
for (size_t k = 0; k < m_kk; k++) {
double c = cc[speciesName(k)];
if (c > 0.0) {
ifound = true;
cv[k] = c;
@ -367,9 +359,8 @@ void SurfPhase::_updateThermo(bool force) const
m_spthermo->update(tnow, DATA_PTR(m_cp0), DATA_PTR(m_h0),
DATA_PTR(m_s0));
m_tlast = tnow;
doublereal rt = GasConstant * tnow;
for (size_t k = 0; k < m_kk; k++) {
m_h0[k] *= rt;
m_h0[k] *= GasConstant * tnow;
m_s0[k] *= GasConstant;
m_cp0[k] *= GasConstant;
m_mu0[k] = m_h0[k] - tnow*m_s0[k];

3
src/thermo/ThermoFactory.cpp
View file

@ -166,13 +166,12 @@ ThermoPhase* ThermoFactory::newThermoPhase(const std::string& model)
std::string eosTypeString(int ieos, int length)
{
std::string ss = "UnknownPhaseType";
for (int n = 0; n < ntypes; n++) {
if (_itypes[n] == ieos) {
return _types[n];
}
}
return ss;
return "UnknownPhaseType";
}
ThermoPhase* newPhase(XML_Node& xmlphase)

26
src/thermo/ThermoPhase.cpp
View file

@ -255,8 +255,6 @@ void ThermoPhase::setState_HPorUV(doublereal Htarget, doublereal p,
double Ttop = Tnew;
double Hbot = Hnew;
double Tbot = Tnew;
double Told = Tnew;
double Hold = Hnew;
bool ignoreBounds = false;
// Unstable phases are those for which
@ -270,8 +268,8 @@ void ThermoPhase::setState_HPorUV(doublereal Htarget, doublereal p,
// Newton iteration
for (int n = 0; n < 500; n++) {
Told = Tnew;
Hold = Hnew;
double Told = Tnew;
double Hold = Hnew;
double cpd = Cpnew;
if (cpd < 0.0) {
unstablePhase = true;
@ -457,8 +455,6 @@ void ThermoPhase::setState_SPorSV(doublereal Starget, doublereal p,
double Ttop = Tnew;
double Sbot = Snew;
double Tbot = Tnew;
double Told = Tnew;
double Sold = Snew;
bool ignoreBounds = false;
// Unstable phases are those for which
@ -470,8 +466,8 @@ void ThermoPhase::setState_SPorSV(doublereal Starget, doublereal p,
// Newton iteration
for (int n = 0; n < 500; n++) {
Told = Tnew;
Sold = Snew;
double Told = Tnew;
double Sold = Snew;
double cpd = Cpnew;
if (cpd < 0.0) {
unstablePhase = true;
@ -653,8 +649,7 @@ void ThermoPhase::initThermoFile(const std::string& inputFile,
void ThermoPhase::initThermoXML(XML_Node& phaseNode, const std::string& id)
{
if (phaseNode.hasChild("state")) {
XML_Node& stateNode = phaseNode.child("state");
setStateFromXML(stateNode);
setStateFromXML(phaseNode.child("state"));
}
setReferenceComposition(0);
}
@ -958,10 +953,9 @@ std::string ThermoPhase::report(bool show_thermo, doublereal threshold) const
}
}
size_t kk = nSpecies();
vector_fp x(kk);
vector_fp y(kk);
vector_fp mu(kk);
vector_fp x(m_kk);
vector_fp y(m_kk);
vector_fp mu(m_kk);
getMoleFractions(&x[0]);
getMassFractions(&y[0]);
getChemPotentials(&mu[0]);
@ -977,7 +971,7 @@ std::string ThermoPhase::report(bool show_thermo, doublereal threshold) const
sprintf(p, " ------------- "
"------------ ------------\n");
s += p;
for (size_t k = 0; k < kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
if (x[k] >= threshold) {
if (x[k] > SmallNumber) {
sprintf(p, "%18s %12.6g %12.6g %12.6g\n",
@ -1000,7 +994,7 @@ std::string ThermoPhase::report(bool show_thermo, doublereal threshold) const
sprintf(p, " -------------"
" ------------\n");
s += p;
for (size_t k = 0; k < kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
if (x[k] >= threshold) {
sprintf(p, "%18s %12.6g %12.6g\n",
speciesName(k).c_str(), x[k], y[k]);

18
src/thermo/VPSSMgr.cpp
View file

@ -150,8 +150,7 @@ VPSSMgr::getStandardChemPotentials(doublereal* mu) const
{
if (m_useTmpStandardStateStorage) {
std::copy(m_gss_RT.begin(), m_gss_RT.end(), mu);
doublereal _rt = GasConstant * m_tlast;
scale(mu, mu+m_kk, mu, _rt);
scale(mu, mu+m_kk, mu, GasConstant * m_tlast);
} else {
throw NotImplementedError("VPSSMgr::getStandardChemPotentials");
}
@ -192,9 +191,8 @@ VPSSMgr::getIntEnergy_RT(doublereal* urt) const
{
if (m_useTmpStandardStateStorage) {
std::copy(m_hss_RT.begin(), m_hss_RT.end(), urt);
doublereal pRT = m_plast / (GasConstant * m_tlast);
for (size_t k = 0; k < m_kk; k++) {
urt[k] -= pRT * m_Vss[k];
urt[k] -= m_plast / (GasConstant * m_tlast) * m_Vss[k];
}
} else {
throw NotImplementedError("VPSSMgr::getEntropy_RT");
@ -255,8 +253,7 @@ VPSSMgr::getGibbs_ref(doublereal* g) const
{
if (m_useTmpRefStateStorage) {
std::copy(m_g0_RT.begin(), m_g0_RT.end(), g);
doublereal _rt = GasConstant * m_tlast;
scale(g, g+m_kk, g, _rt);
scale(g, g+m_kk, g, GasConstant * m_tlast);
} else {
throw NotImplementedError("VPSSMgr::getGibbs_ref");
}
@ -420,8 +417,7 @@ PDSS* VPSSMgr::createInstallPDSS(size_t k, const XML_Node& s,
doublereal VPSSMgr::minTemp(size_t k) const
{
if (k != npos) {
const PDSS* kPDSS = m_vptp_ptr->providePDSS(k);
return kPDSS->minTemp();
return m_vptp_ptr->providePDSS(k)->minTemp();
}
return m_minTemp;
}
@ -429,8 +425,7 @@ doublereal VPSSMgr::minTemp(size_t k) const
doublereal VPSSMgr::maxTemp(size_t k) const
{
if (k != npos) {
const PDSS* kPDSS = m_vptp_ptr->providePDSS(k);
return kPDSS->maxTemp();
return m_vptp_ptr->providePDSS(k)->maxTemp();
}
return m_maxTemp;
}
@ -438,8 +433,7 @@ doublereal VPSSMgr::maxTemp(size_t k) const
doublereal VPSSMgr::refPressure(size_t k) const
{
if (k != npos) {
const PDSS* kPDSS = m_vptp_ptr->providePDSS(k);
return kPDSS->refPressure();
return m_vptp_ptr->providePDSS(k)->refPressure();
}
return m_p0;
}

17
src/thermo/VPSSMgrFactory.cpp
View file

@ -67,15 +67,12 @@ static void getVPSSMgrTypes(std::vector<XML_Node*> & spDataNodeList,
int& has_other)
{
XML_Node* ss_ptr = 0;
string ssModel = "idealGas";
size_t ns = spDataNodeList.size();
for (size_t n = 0; n < ns; n++) {
for (size_t n = 0; n < spDataNodeList.size(); n++) {
bool ifound = false;
XML_Node* spNode = spDataNodeList[n];
if (spNode->hasChild("standardState")) {
const XML_Node& ssN = spNode->child("standardState");
string mm = ssN["model"];
string mm = spNode->child("standardState")["model"];
if (mm == "waterIAPWS" || mm == "waterPDSS") {
has_water++;
ifound = true;
@ -89,8 +86,7 @@ static void getVPSSMgrTypes(std::vector<XML_Node*> & spDataNodeList,
if (spNode->hasChild("thermo")) {
const XML_Node& th = spNode->child("thermo");
if (spNode->hasChild("standardState")) {
ss_ptr = &(spNode->child("standardState"));
ssModel = ss_ptr->attrib("model");
ssModel = spNode->child("standardState")["model"];
}
if (th.hasChild("NASA")) {
if (ssModel == "idealGas") {
@ -211,8 +207,8 @@ VPSSMgrFactory::newVPSSMgr(VPStandardStateTP* vp_ptr,
std::vector<XML_Node*> & spDataNodeList)
{
std::string ssManager="";
std::string vpssManager="";
std::string ssManager;
std::string vpssManager;
// First look for any explicit instructions within the XML Database
// for the standard state manager and the variable pressure
@ -251,8 +247,7 @@ VPSSMgrFactory::newVPSSMgr(VPStandardStateTP* vp_ptr,
int inasaIG = 0, inasaCV = 0, ishomateIG = 0, ishomateCV = 0,
isimpleIG = 0, isimpleCV = 0,
iwater = 0, itpx = 0, iother = 0;
isimpleIG = 0, isimpleCV = 0, iwater = 0, itpx = 0, iother = 0;
int ihptx = 0;
try {

11
src/thermo/VPSSMgr_ConstVol.cpp
View file

@ -104,13 +104,12 @@ VPSSMgr_ConstVol::initThermoXML(XML_Node& phaseNode, const std::string& id)
XML_Node& speciesList = phaseNode.child("speciesArray");
XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"],
&phaseNode.root());
const vector<string>&sss = m_vptp_ptr->speciesNames();
for (size_t k = 0; k < m_kk; k++) {
const XML_Node* s = speciesDB->findByAttr("name", sss[k]);
const XML_Node* s = speciesDB->findByAttr("name", m_vptp_ptr->speciesName(k));
if (!s) {
throw CanteraError("VPSSMgr_ConstVol::initThermoXML",
"no species Node for species " + sss[k]);
"no species Node for species " + m_vptp_ptr->speciesName(k));
}
const XML_Node* ss = s->findByName("standardState");
if (!ss) {
@ -147,11 +146,7 @@ VPSSMgr_ConstVol::createInstallPDSS(size_t k, const XML_Node& speciesNode,
m_Vss[k] = ctml::getFloat(*ss, "molarVolume", "toSI");
installSTSpecies(k, speciesNode, phaseNode_ptr);
PDSS* kPDSS = new PDSS_ConstVol(m_vptp_ptr, k, speciesNode,
*phaseNode_ptr, true);
return kPDSS;
return new PDSS_ConstVol(m_vptp_ptr, k, speciesNode, *phaseNode_ptr, true);
}
PDSS_enumType VPSSMgr_ConstVol::reportPDSSType(int k) const

8
src/thermo/VPSSMgr_General.cpp
View file

@ -120,17 +120,16 @@ void VPSSMgr_General::initThermo()
void VPSSMgr_General::getGibbs_ref(doublereal* g) const
{
doublereal _rt = GasConstant * m_tlast;
if (m_useTmpRefStateStorage) {
std::copy(m_g0_RT.begin(), m_g0_RT.end(), g);
scale(g, g+m_kk, g, _rt);
scale(g, g+m_kk, g, GasConstant * m_tlast);
} else {
for (size_t k = 0; k < m_kk; k++) {
PDSS* kPDSS = m_PDSS_ptrs[k];
kPDSS->setState_TP(m_tlast, m_plast);
double h0_RT = kPDSS->enthalpy_RT_ref();
double s0_R = kPDSS->entropy_R_ref();
g[k] = _rt * (h0_RT - s0_R);
g[k] = GasConstant * m_tlast * (h0_RT - s0_R);
}
}
}
@ -229,8 +228,7 @@ VPSSMgr_General::createInstallPDSS(size_t k, const XML_Node& speciesNode,
PDSS_enumType VPSSMgr_General::reportPDSSType(int k) const
{
PDSS* kPDSS = m_PDSS_ptrs[k];
return kPDSS->reportPDSSType();
return m_PDSS_ptrs[k]->reportPDSSType();
}
VPSSMgr_enumType VPSSMgr_General::reportVPSSMgrType() const

3
src/thermo/VPSSMgr_IdealGas.cpp
View file

@ -93,8 +93,7 @@ VPSSMgr_IdealGas::createInstallPDSS(size_t k, const XML_Node& speciesNode,
{
const XML_Node* ss = speciesNode.findByName("standardState");
if (ss) {
std::string model = ss->attrib("model");
if (model != "ideal_gas") {
if (ss->attrib("model") != "ideal_gas") {
throw CanteraError("VPSSMgr_IdealGas::createInstallPDSS",
"standardState model for species isn't "
"ideal_gas: " + speciesNode["name"]);

38
src/thermo/VPSSMgr_Water_ConstVol.cpp
View file

@ -76,9 +76,8 @@ VPSSMgr_Water_ConstVol::getEnthalpy_RT_ref(doublereal* hrt) const
// Everything should be OK except for the water SS
m_p0 = m_waterSS->pref_safe(m_tlast);
if (m_p0 != m_plast) {
doublereal RT = GasConstant * m_tlast;
m_waterSS->setState_TP(m_tlast, m_p0);
m_h0_RT[0] = (m_waterSS->enthalpy_mole()) / RT;
m_h0_RT[0] = (m_waterSS->enthalpy_mole()) / (GasConstant * m_tlast);
m_waterSS->setState_TP(m_tlast, m_plast);
} else {
m_h0_RT[0] = m_hss_RT[0];
@ -92,9 +91,8 @@ VPSSMgr_Water_ConstVol::getGibbs_RT_ref(doublereal* grt) const
// Everything should be OK except for the water SS
m_p0 = m_waterSS->pref_safe(m_tlast);
if (m_p0 != m_plast) {
doublereal RT = GasConstant * m_tlast;
m_waterSS->setState_TP(m_tlast, m_p0);
m_g0_RT[0] = (m_waterSS->gibbs_mole()) / RT;
m_g0_RT[0] = (m_waterSS->gibbs_mole()) / (GasConstant * m_tlast);
m_waterSS->setState_TP(m_tlast, m_plast);
} else {
m_g0_RT[0] = m_gss_RT[0];
@ -105,10 +103,9 @@ VPSSMgr_Water_ConstVol::getGibbs_RT_ref(doublereal* grt) const
void
VPSSMgr_Water_ConstVol::getGibbs_ref(doublereal* g) const
{
doublereal RT = GasConstant * m_tlast;
getGibbs_RT_ref(g);
for (size_t k = 0; k < m_kk; k++) {
g[k] *= RT;
g[k] *= GasConstant * m_tlast;
}
}
@ -163,12 +160,10 @@ void VPSSMgr_Water_ConstVol::_updateRefStateThermo() const
m_spthermo->update(m_tlast, &m_cp0_R[0], &m_h0_RT[0], &m_s0_R[0]);
for (size_t k = 0; k < m_kk; k++) {
m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
PDSS* kPDSS = m_vptp_ptr->providePDSS(k);
kPDSS->setTemperature(m_tlast);
m_vptp_ptr->providePDSS(k)->setTemperature(m_tlast);
}
doublereal RT = GasConstant * m_tlast;
m_waterSS->setState_TP(m_tlast, m_p0);
m_h0_RT[0] = (m_waterSS->enthalpy_mole())/ RT;
m_h0_RT[0] = (m_waterSS->enthalpy_mole()) / (GasConstant * m_tlast);
m_s0_R[0] = (m_waterSS->entropy_mole()) / GasConstant;
m_cp0_R[0] = (m_waterSS->cp_mole()) / GasConstant;
m_g0_RT[0] = (m_hss_RT[0] - m_sss_R[0]);
@ -178,8 +173,7 @@ void VPSSMgr_Water_ConstVol::_updateRefStateThermo() const
void VPSSMgr_Water_ConstVol::_updateStandardStateThermo()
{
doublereal RT = GasConstant * m_tlast;
doublereal del_pRT = (m_plast - OneAtm) / (RT);
doublereal del_pRT = (m_plast - OneAtm) / (GasConstant * m_tlast);
for (size_t k = 1; k < m_kk; k++) {
m_hss_RT[k] = m_h0_RT[k] + del_pRT * m_Vss[k];
@ -192,7 +186,7 @@ void VPSSMgr_Water_ConstVol::_updateStandardStateThermo()
}
// Do the water
m_waterSS->setState_TP(m_tlast, m_plast);
m_hss_RT[0] = (m_waterSS->enthalpy_mole())/ RT;
m_hss_RT[0] = (m_waterSS->enthalpy_mole()) / (GasConstant * m_tlast);
m_sss_R[0] = (m_waterSS->entropy_mole()) / GasConstant;
m_cpss_R[0] = (m_waterSS->cp_mole()) / GasConstant;
m_gss_RT[0] = (m_hss_RT[0] - m_sss_R[0]);
@ -212,8 +206,6 @@ VPSSMgr_Water_ConstVol::initThermoXML(XML_Node& phaseNode, const std::string& id
XML_Node& speciesList = phaseNode.child("speciesArray");
XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"],
&phaseNode.root());
const vector<string>&sss = m_vptp_ptr->speciesNames();
if (!m_waterSS) {
throw CanteraError("VPSSMgr_Water_ConstVol::initThermoXML",
@ -224,23 +216,20 @@ VPSSMgr_Water_ConstVol::initThermoXML(XML_Node& phaseNode, const std::string& id
m_Vss[0] = (m_waterSS->density()) / m_vptp_ptr->molecularWeight(0);
for (size_t k = 1; k < m_kk; k++) {
const XML_Node* s = speciesDB->findByAttr("name", sss[k]);
const XML_Node* s = speciesDB->findByAttr("name", m_vptp_ptr->speciesName(k));
if (!s) {
throw CanteraError("VPSSMgr_Water_ConstVol::initThermoXML",
"no species Node for species " + sss[k]);
"no species Node for species " + m_vptp_ptr->speciesName(k));
}
const XML_Node* ss = s->findByName("standardState");
if (!ss) {
std::string sName = s->attrib("name");
throw CanteraError("VPSSMgr_Water_ConstVol::initThermoXML",
"no standardState Node for species " + sName);
"no standardState Node for species " + s->attrib("name"));
}
std::string model = ss->attrib("model");
if (model != "constant_incompressible") {
std::string sName = s->attrib("name");
if (ss->attrib("model") != "constant_incompressible") {
throw CanteraError("VPSSMgr_Water_ConstVol::initThermoXML",
"standardState model for species isn't "
"constant_incompressible: " + sName);
"constant_incompressible: " + s->attrib("name"));
}
m_Vss[k] = ctml::getFloat(*ss, "molarVolume", "toSI");
}
@ -284,8 +273,7 @@ VPSSMgr_Water_ConstVol::createInstallPDSS(size_t k, const XML_Node& speciesNode,
throw CanteraError("VPSSMgr_Water_ConstVol::installSpecies",
"no standardState Node for species " + speciesNode.name());
}
std::string model = ss->attrib("model");
if (model != "constant_incompressible") {
if (ss->attrib("model") != "constant_incompressible") {
throw CanteraError("VPSSMgr_Water_ConstVol::initThermoXML",
"standardState model for species isn't "
"constant_incompressible: " + speciesNode.name());

34
src/thermo/VPSSMgr_Water_HKFT.cpp
View file

@ -83,9 +83,8 @@ void
VPSSMgr_Water_HKFT::getGibbs_ref(doublereal* g) const
{
getGibbs_RT_ref(g);
doublereal RT = GasConstant * m_tlast;
for (size_t k = 0; k < m_kk; k++) {
g[k] *= RT;
g[k] *= GasConstant * m_tlast;
}
}
@ -149,9 +148,8 @@ void VPSSMgr_Water_HKFT::updateRefStateThermo() const
void VPSSMgr_Water_HKFT::_updateRefStateThermo() const
{
m_p0 = m_waterSS->pref_safe(m_tlast);
doublereal RT = GasConstant * m_tlast;
m_waterSS->setState_TP(m_tlast, m_p0);
m_h0_RT[0] = (m_waterSS->enthalpy_mole())/ RT;
m_h0_RT[0] = (m_waterSS->enthalpy_mole()) / (GasConstant * m_tlast);
m_s0_R[0] = (m_waterSS->entropy_mole()) / GasConstant;
m_cp0_R[0] = (m_waterSS->cp_mole()) / GasConstant;
m_g0_RT[0] = (m_hss_RT[0] - m_sss_R[0]);
@ -182,10 +180,9 @@ void VPSSMgr_Water_HKFT::_updateRefStateThermo() const
void VPSSMgr_Water_HKFT::_updateStandardStateThermo()
{
doublereal RT = GasConstant * m_tlast;
// Do the water
m_waterSS->setState_TP(m_tlast, m_plast);
m_hss_RT[0] = (m_waterSS->enthalpy_mole())/ RT;
m_hss_RT[0] = (m_waterSS->enthalpy_mole()) / (GasConstant * m_tlast);
m_sss_R[0] = (m_waterSS->entropy_mole()) / GasConstant;
m_cpss_R[0] = (m_waterSS->cp_mole()) / GasConstant;
m_gss_RT[0] = (m_hss_RT[0] - m_sss_R[0]);
@ -216,27 +213,26 @@ VPSSMgr_Water_HKFT::initThermoXML(XML_Node& phaseNode, const std::string& id)
XML_Node& speciesList = phaseNode.child("speciesArray");
XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"],
&phaseNode.root());
const vector<string> &sss = m_vptp_ptr->speciesNames();
m_waterSS->setState_TP(300., OneAtm);
m_Vss[0] = (m_waterSS->density()) / m_vptp_ptr->molecularWeight(0);
for (size_t k = 1; k < m_kk; k++) {
const XML_Node* s = speciesDB->findByAttr("name", sss[k]);
string name = m_vptp_ptr->speciesName(k);
const XML_Node* s = speciesDB->findByAttr("name", name);
if (!s) {
throw CanteraError("VPSSMgr_Water_HKFT::initThermoXML",
"No species Node for species " + sss[k]);
"No species Node for species " + name);
}
const XML_Node* ss = s->findByName("standardState");
if (!ss) {
throw CanteraError("VPSSMgr_Water_HKFT::initThermoXML",
"No standardState Node for species " + sss[k]);
"No standardState Node for species " + name);
}
std::string model = lowercase(ss->attrib("model"));
if (model != "hkft") {
throw CanteraError("VPSSMgr_Water_HKFT::initThermoXML",
"Standard state model for a solute species isn't "
"the HKFT standard state model: " + sss[k]);
"the HKFT standard state model: " + name);
}
}
}
@ -249,18 +245,16 @@ VPSSMgr_Water_HKFT::createInstallPDSS(size_t k, const XML_Node& speciesNode,
const XML_Node* ss = speciesNode.findByName("standardState");
if (!ss) {
std::string sName = speciesNode["name"];
throw CanteraError("VPSSMgr_Water_HKFT::installSpecies",
"No standardState Node for species " + sName);
"No standardState Node for species " + speciesNode["name"]);
}
// Will have to do something for water
// -> make sure it's species 0
// -> make sure it's designated as a real water EOS
if (k == 0) {
string xn = speciesNode["name"];
if (xn != "H2O(L)") {
if (speciesNode["name"] != "H2O(L)") {
throw CanteraError("VPSSMgr_Water_HKFT::installSpecies",
"h2o wrong name: " + xn);
"h2o wrong name: " + speciesNode["name"]);
}
std::string model = ss->attrib("model");
@ -281,12 +275,10 @@ VPSSMgr_Water_HKFT::createInstallPDSS(size_t k, const XML_Node& speciesNode,
kPDSS = m_waterSS;
} else {
std::string model = ss->attrib("model");
if (model != "HKFT") {
std::string sName = speciesNode["name"];
if (ss->attrib("model") != "HKFT") {
throw CanteraError("VPSSMgr_Water_HKFT::initThermoXML",
"standardState model for species isn't "
"HKFT: " + sName);
"HKFT: " + speciesNode["name"]);
}
kPDSS = new PDSS_HKFT(m_vptp_ptr, k, speciesNode, *phaseNode_ptr, true);

6
src/thermo/VPStandardStateTP.cpp
View file

@ -68,8 +68,7 @@ VPStandardStateTP::operator=(const VPStandardStateTP& b)
}
m_PDSS_storage.resize(m_kk);
for (size_t k = 0; k < m_kk; k++) {
PDSS* ptmp = b.m_PDSS_storage[k];
m_PDSS_storage[k] = ptmp->duplMyselfAsPDSS();
m_PDSS_storage[k] = b.m_PDSS_storage[k]->duplMyselfAsPDSS();
}
/*
@ -91,8 +90,7 @@ VPStandardStateTP::operator=(const VPStandardStateTP& b)
* so it occurs after m_VPSS_ptr is set.
*/
for (size_t k = 0; k < m_kk; k++) {
PDSS* ptmp = m_PDSS_storage[k];
ptmp->initAllPtrs(this, m_VPSS_ptr, m_spthermo);
m_PDSS_storage[k]->initAllPtrs(this, m_VPSS_ptr, m_spthermo);
}
/*
* Ok, the VPSSMgr object is ready for business.

26
src/thermo/WaterSSTP.cpp
View file

@ -197,28 +197,23 @@ void WaterSSTP::setParametersFromXML(const XML_Node& eosdata)
void WaterSSTP::getEnthalpy_RT(doublereal* hrt) const
{
double T = temperature();
doublereal h = m_sub->enthalpy();
*hrt = (h + EW_Offset)/(GasConstant*T);
*hrt = (m_sub->enthalpy() + EW_Offset)/(GasConstant*temperature());
}
void WaterSSTP::getIntEnergy_RT(doublereal* ubar) const
{
doublereal u = m_sub->intEnergy();
*ubar = (u + EW_Offset)/GasConstant;
*ubar = (m_sub->intEnergy() + EW_Offset)/GasConstant;
}
void WaterSSTP::getEntropy_R(doublereal* sr) const
{
doublereal s = m_sub->entropy();
sr[0] = (s + SW_Offset) / GasConstant;
sr[0] = (m_sub->entropy() + SW_Offset) / GasConstant;
}
void WaterSSTP::getGibbs_RT(doublereal* grt) const
{
double T = temperature();
doublereal g = m_sub->Gibbs();
*grt = (g + EW_Offset - SW_Offset*T) / (GasConstant * T);
*grt = (m_sub->Gibbs() + EW_Offset - SW_Offset*T) / (GasConstant * T);
if (!m_ready) {
throw CanteraError("waterSSTP::", "Phase not ready");
}
@ -226,9 +221,7 @@ void WaterSSTP::getGibbs_RT(doublereal* grt) const
void WaterSSTP::getStandardChemPotentials(doublereal* gss) const
{
double T = temperature();
doublereal g = m_sub->Gibbs();
*gss = (g + EW_Offset - SW_Offset*T);
*gss = (m_sub->Gibbs() + EW_Offset - SW_Offset*temperature());
if (!m_ready) {
throw CanteraError("waterSSTP::", "Phase not ready");
}
@ -236,8 +229,7 @@ void WaterSSTP::getStandardChemPotentials(doublereal* gss) const
void WaterSSTP::getCp_R(doublereal* cpr) const
{
doublereal cp = m_sub->cp();
cpr[0] = cp / GasConstant;
cpr[0] = m_sub->cp() / GasConstant;
}
doublereal WaterSSTP::cv_mole() const
@ -421,15 +413,13 @@ doublereal WaterSSTP::critDensity() const
void WaterSSTP::setTemperature(const doublereal temp)
{
Phase::setTemperature(temp);
doublereal dd = density();
m_sub->setState_TR(temp, dd);
m_sub->setState_TR(temp, density());
}
void WaterSSTP::setDensity(const doublereal dens)
{
Phase::setDensity(dens);
doublereal temp = temperature();
m_sub->setState_TR(temp, dens);
m_sub->setState_TR(temperature(), dens);
}
doublereal WaterSSTP::satPressure(doublereal t) {