1916 lines
71 KiB
C++
1916 lines
71 KiB
C++
/**
|
|
* @file HMWSoln_input.cpp
|
|
* Definitions for the %HMWSoln ThermoPhase object, which models concentrated
|
|
* electrolyte solutions
|
|
* (see \ref thermoprops and \link Cantera::HMWSoln HMWSoln \endlink) .
|
|
*
|
|
* This file contains definitions for reading in the interaction terms
|
|
* in the formulation.
|
|
*/
|
|
/*
|
|
* Copyright (2006) Sandia Corporation. Under the terms of
|
|
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
|
|
* U.S. Government retains certain rights in this software.
|
|
*/
|
|
#include "cantera/thermo/HMWSoln.h"
|
|
#include "cantera/thermo/ThermoFactory.h"
|
|
#include "cantera/thermo/WaterProps.h"
|
|
#include "cantera/thermo/PDSS_Water.h"
|
|
#include "cantera/base/stringUtils.h"
|
|
|
|
#include <cstring>
|
|
#include <cstdlib>
|
|
#include <cstdio>
|
|
#include <fstream>
|
|
|
|
using namespace std;
|
|
using namespace ctml;
|
|
|
|
namespace Cantera
|
|
{
|
|
|
|
|
|
//! utility function to assign an integer value from a string
|
|
//! for the ElectrolyteSpeciesType field.
|
|
/*!
|
|
* @param estString string name of the electrolyte species type
|
|
*/
|
|
int HMWSoln::interp_est(const std::string& estString)
|
|
{
|
|
const char* cc = estString.c_str();
|
|
string lcs = lowercase(estString);
|
|
const char* ccl = lcs.c_str();
|
|
if (!strcmp(ccl, "solvent")) {
|
|
return cEST_solvent;
|
|
} else if (!strcmp(ccl, "chargedspecies")) {
|
|
return cEST_chargedSpecies;
|
|
} else if (!strcmp(ccl, "weakacidassociated")) {
|
|
return cEST_weakAcidAssociated;
|
|
} else if (!strcmp(ccl, "strongacidassociated")) {
|
|
return cEST_strongAcidAssociated;
|
|
} else if (!strcmp(ccl, "polarneutral")) {
|
|
return cEST_polarNeutral;
|
|
} else if (!strcmp(ccl, "nonpolarneutral")) {
|
|
return cEST_nonpolarNeutral;
|
|
}
|
|
int retn, rval;
|
|
if ((retn = sscanf(cc, "%d", &rval)) != 1) {
|
|
return -1;
|
|
}
|
|
return rval;
|
|
}
|
|
|
|
/*
|
|
* Process an XML node called "SimpleSaltParameters.
|
|
* This node contains all of the parameters necessary to describe
|
|
* the Pitzer model for that particular binary salt.
|
|
* This function reads the XML file and writes the coefficients
|
|
* it finds to an internal data structures.
|
|
*/
|
|
void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
|
|
{
|
|
string xname = BinSalt.name();
|
|
if (xname != "binarySaltParameters") {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt",
|
|
"Incorrect name for processing this routine: " + xname);
|
|
}
|
|
double* charge = DATA_PTR(m_speciesCharge);
|
|
string stemp;
|
|
size_t nParamsFound, i;
|
|
vector_fp vParams;
|
|
string iName = BinSalt.attrib("cation");
|
|
if (iName == "") {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt", "no cation attrib");
|
|
}
|
|
string jName = BinSalt.attrib("anion");
|
|
if (jName == "") {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt", "no anion attrib");
|
|
}
|
|
/*
|
|
* Find the index of the species in the current phase. It's not
|
|
* an error to not find the species
|
|
*/
|
|
size_t iSpecies = speciesIndex(iName);
|
|
if (iSpecies == npos) {
|
|
return;
|
|
}
|
|
string ispName = speciesName(iSpecies);
|
|
if (charge[iSpecies] <= 0) {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt", "cation charge problem");
|
|
}
|
|
size_t jSpecies = speciesIndex(jName);
|
|
if (jSpecies == npos) {
|
|
return;
|
|
}
|
|
string jspName = speciesName(jSpecies);
|
|
if (charge[jSpecies] >= 0) {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt", "anion charge problem");
|
|
}
|
|
|
|
size_t n = iSpecies * m_kk + jSpecies;
|
|
int counter = m_CounterIJ[n];
|
|
for (size_t iChild = 0; iChild < BinSalt.nChildren(); iChild++) {
|
|
XML_Node& xmlChild = BinSalt.child(iChild);
|
|
stemp = xmlChild.name();
|
|
string nodeName = lowercase(stemp);
|
|
/*
|
|
* Process the binary salt child elements
|
|
*/
|
|
if (nodeName == "beta0") {
|
|
/*
|
|
* Get the string containing all of the values
|
|
*/
|
|
getFloatArray(xmlChild, vParams, false, "", "beta0");
|
|
nParamsFound = vParams.size();
|
|
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
|
|
if (nParamsFound != 1) {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt::beta0 for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Beta0MX_ij[counter] = vParams[0];
|
|
m_Beta0MX_ij_coeff(0,counter) = m_Beta0MX_ij[counter];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
|
|
if (nParamsFound != 2) {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt::beta0 for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Beta0MX_ij_coeff(0,counter) = vParams[0];
|
|
m_Beta0MX_ij_coeff(1,counter) = vParams[1];
|
|
m_Beta0MX_ij[counter] = vParams[0];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
|
|
if (nParamsFound != 5) {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt::beta0 for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
for (i = 0; i < nParamsFound; i++) {
|
|
m_Beta0MX_ij_coeff(i, counter) = vParams[i];
|
|
}
|
|
m_Beta0MX_ij[counter] = vParams[0];
|
|
}
|
|
}
|
|
if (nodeName == "beta1") {
|
|
|
|
/*
|
|
* Get the string containing all of the values
|
|
*/
|
|
getFloatArray(xmlChild, vParams, false, "", "beta1");
|
|
nParamsFound = vParams.size();
|
|
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
|
|
if (nParamsFound != 1) {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt::beta1 for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Beta1MX_ij[counter] = vParams[0];
|
|
m_Beta1MX_ij_coeff(0,counter) = m_Beta1MX_ij[counter];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
|
|
if (nParamsFound != 2) {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt::beta1 for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Beta1MX_ij_coeff(0,counter) = vParams[0];
|
|
m_Beta1MX_ij_coeff(1,counter) = vParams[1];
|
|
m_Beta1MX_ij[counter] = vParams[0];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
|
|
if (nParamsFound != 5) {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt::beta1 for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
for (i = 0; i < nParamsFound; i++) {
|
|
m_Beta1MX_ij_coeff(i, counter) = vParams[i];
|
|
}
|
|
m_Beta1MX_ij[counter] = vParams[0];
|
|
}
|
|
}
|
|
if (nodeName == "beta2") {
|
|
getFloatArray(xmlChild, vParams, false, "", "beta2");
|
|
nParamsFound = vParams.size();
|
|
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
|
|
if (nParamsFound != 1) {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt::beta2 for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Beta2MX_ij[counter] = vParams[0];
|
|
m_Beta2MX_ij_coeff(0,counter) = m_Beta2MX_ij[counter];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
|
|
if (nParamsFound != 2) {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt::beta2 for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Beta2MX_ij_coeff(0,counter) = vParams[0];
|
|
m_Beta2MX_ij_coeff(1,counter) = vParams[1];
|
|
m_Beta2MX_ij[counter] = vParams[0];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
|
|
if (nParamsFound != 5) {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt::beta2 for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
for (i = 0; i < nParamsFound; i++) {
|
|
m_Beta2MX_ij_coeff(i, counter) = vParams[i];
|
|
}
|
|
m_Beta2MX_ij[counter] = vParams[0];
|
|
}
|
|
|
|
}
|
|
if (nodeName == "cphi") {
|
|
/*
|
|
* Get the string containing all of the values
|
|
*/
|
|
getFloatArray(xmlChild, vParams, false, "", "Cphi");
|
|
nParamsFound = vParams.size();
|
|
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
|
|
if (nParamsFound != 1) {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt::Cphi for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
m_CphiMX_ij[counter] = vParams[0];
|
|
m_CphiMX_ij_coeff(0,counter) = m_CphiMX_ij[counter];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
|
|
if (nParamsFound != 2) {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt::Cphi for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
m_CphiMX_ij_coeff(0,counter) = vParams[0];
|
|
m_CphiMX_ij_coeff(1,counter) = vParams[1];
|
|
m_CphiMX_ij[counter] = vParams[0];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
|
|
if (nParamsFound != 5) {
|
|
throw CanteraError("HMWSoln::readXMLBinarySalt::Cphi for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
for (i = 0; i < nParamsFound; i++) {
|
|
m_CphiMX_ij_coeff(i, counter) = vParams[i];
|
|
}
|
|
m_CphiMX_ij[counter] = vParams[0];
|
|
}
|
|
}
|
|
|
|
if (nodeName == "alpha1") {
|
|
stemp = xmlChild.value();
|
|
m_Alpha1MX_ij[counter] = atofCheck(stemp.c_str());
|
|
}
|
|
|
|
if (nodeName == "alpha2") {
|
|
stemp = xmlChild.value();
|
|
m_Alpha2MX_ij[counter] = atofCheck(stemp.c_str());
|
|
}
|
|
}
|
|
}
|
|
|
|
/**
|
|
* Process an XML node called "thetaAnion".
|
|
* This node contains all of the parameters necessary to describe
|
|
* the binary interactions between two anions.
|
|
*/
|
|
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);
|
|
}
|
|
double* charge = DATA_PTR(m_speciesCharge);
|
|
string stemp;
|
|
string ispName = BinSalt.attrib("anion1");
|
|
if (ispName == "") {
|
|
throw CanteraError("HMWSoln::readXMLThetaAnion", "no anion1 attrib");
|
|
}
|
|
string jspName = BinSalt.attrib("anion2");
|
|
if (jspName == "") {
|
|
throw CanteraError("HMWSoln::readXMLThetaAnion", "no anion2 attrib");
|
|
}
|
|
/*
|
|
* Find the index of the species in the current phase. It's not
|
|
* an error to not find the species
|
|
*/
|
|
size_t iSpecies = speciesIndex(ispName);
|
|
if (iSpecies == npos) {
|
|
return;
|
|
}
|
|
if (charge[iSpecies] >= 0) {
|
|
throw CanteraError("HMWSoln::readXMLThetaAnion", "anion1 charge problem");
|
|
}
|
|
size_t jSpecies = speciesIndex(jspName);
|
|
if (jSpecies == npos) {
|
|
return;
|
|
}
|
|
if (charge[jSpecies] >= 0) {
|
|
throw CanteraError("HMWSoln::readXMLThetaAnion", "anion2 charge problem");
|
|
}
|
|
|
|
size_t n = iSpecies * m_kk + jSpecies;
|
|
int counter = m_CounterIJ[n];
|
|
for (size_t i = 0; i < BinSalt.nChildren(); i++) {
|
|
XML_Node& xmlChild = BinSalt.child(i);
|
|
stemp = xmlChild.name();
|
|
string nodeName = lowercase(stemp);
|
|
if (nodeName == "theta") {
|
|
getFloatArray(xmlChild, vParams, false, "", stemp);
|
|
nParamsFound = vParams.size();
|
|
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
|
|
if (nParamsFound != 1) {
|
|
throw CanteraError("HMWSoln::readXMLThetaAnion::Theta for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Theta_ij_coeff(0,counter) = vParams[0];
|
|
m_Theta_ij[counter] = vParams[0];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
|
|
if (nParamsFound != 2) {
|
|
throw CanteraError("HMWSoln::readXMLThetaAnion::Theta for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Theta_ij_coeff(0,counter) = vParams[0];
|
|
m_Theta_ij_coeff(1,counter) = vParams[1];
|
|
m_Theta_ij[counter] = vParams[0];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
|
|
if (nParamsFound == 1) {
|
|
vParams.resize(5, 0.0);
|
|
nParamsFound = 5;
|
|
} else if (nParamsFound != 5) {
|
|
throw CanteraError("HMWSoln::readXMLThetaAnion::Theta for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Theta_ij_coeff(j, counter) = vParams[j];
|
|
}
|
|
m_Theta_ij[counter] = vParams[0];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/**
|
|
* Process an XML node called "thetaCation".
|
|
* This node contains all of the parameters necessary to describe
|
|
* the binary interactions between two cation.
|
|
*/
|
|
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);
|
|
}
|
|
double* charge = DATA_PTR(m_speciesCharge);
|
|
string stemp;
|
|
string ispName = BinSalt.attrib("cation1");
|
|
if (ispName == "") {
|
|
throw CanteraError("HMWSoln::readXMLThetaCation", "no cation1 attrib");
|
|
}
|
|
string jspName = BinSalt.attrib("cation2");
|
|
if (jspName == "") {
|
|
throw CanteraError("HMWSoln::readXMLThetaCation", "no cation2 attrib");
|
|
}
|
|
/*
|
|
* Find the index of the species in the current phase. It's not
|
|
* an error to not find the species
|
|
*/
|
|
size_t iSpecies = speciesIndex(ispName);
|
|
if (iSpecies == npos) {
|
|
return;
|
|
}
|
|
if (charge[iSpecies] <= 0) {
|
|
throw CanteraError("HMWSoln::readXMLThetaCation", "cation1 charge problem");
|
|
}
|
|
size_t jSpecies = speciesIndex(jspName);
|
|
if (jSpecies == npos) {
|
|
return;
|
|
}
|
|
if (charge[jSpecies] <= 0) {
|
|
throw CanteraError("HMWSoln::readXMLThetaCation", "cation2 charge problem");
|
|
}
|
|
|
|
size_t n = iSpecies * m_kk + jSpecies;
|
|
int counter = m_CounterIJ[n];
|
|
for (size_t i = 0; i < BinSalt.nChildren(); i++) {
|
|
XML_Node& xmlChild = BinSalt.child(i);
|
|
stemp = xmlChild.name();
|
|
string nodeName = lowercase(stemp);
|
|
if (nodeName == "theta") {
|
|
getFloatArray(xmlChild, vParams, false, "", stemp);
|
|
nParamsFound = vParams.size();
|
|
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
|
|
if (nParamsFound != 1) {
|
|
throw CanteraError("HMWSoln::readXMLThetaCation::Theta for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Theta_ij_coeff(0,counter) = vParams[0];
|
|
m_Theta_ij[counter] = vParams[0];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
|
|
if (nParamsFound != 2) {
|
|
throw CanteraError("HMWSoln::readXMLThetaCation::Theta for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Theta_ij_coeff(0,counter) = vParams[0];
|
|
m_Theta_ij_coeff(1,counter) = vParams[1];
|
|
m_Theta_ij[counter] = vParams[0];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
|
|
if (nParamsFound == 1) {
|
|
vParams.resize(5, 0.0);
|
|
nParamsFound = 5;
|
|
} else if (nParamsFound != 5) {
|
|
throw CanteraError("HMWSoln::readXMLThetaCation::Theta for " + ispName
|
|
+ "::" + jspName,
|
|
"wrong number of params found");
|
|
}
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Theta_ij_coeff(j, counter) = vParams[j];
|
|
}
|
|
m_Theta_ij[counter] = vParams[0];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Process an XML node called "readXMLPsiCommonCation".
|
|
* This node contains all of the parameters necessary to describe
|
|
* the binary interactions between two anions and one common cation.
|
|
*/
|
|
void HMWSoln::readXMLPsiCommonCation(XML_Node& BinSalt)
|
|
{
|
|
string xname = BinSalt.name();
|
|
if (xname != "psiCommonCation") {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonCation",
|
|
"Incorrect name for processing this routine: " + xname);
|
|
}
|
|
double* charge = DATA_PTR(m_speciesCharge);
|
|
string stemp;
|
|
vector_fp vParams;
|
|
size_t nParamsFound = 0;
|
|
string kName = BinSalt.attrib("cation");
|
|
if (kName == "") {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonCation", "no cation attrib");
|
|
}
|
|
string iName = BinSalt.attrib("anion1");
|
|
if (iName == "") {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonCation", "no anion1 attrib");
|
|
}
|
|
string jName = BinSalt.attrib("anion2");
|
|
if (jName == "") {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonCation", "no anion2 attrib");
|
|
}
|
|
/*
|
|
* Find the index of the species in the current phase. It's not
|
|
* an error to not find the species
|
|
*/
|
|
size_t kSpecies = speciesIndex(kName);
|
|
if (kSpecies == npos) {
|
|
return;
|
|
}
|
|
if (charge[kSpecies] <= 0) {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonCation",
|
|
"cation charge problem");
|
|
}
|
|
size_t iSpecies = speciesIndex(iName);
|
|
if (iSpecies == npos) {
|
|
return;
|
|
}
|
|
if (charge[iSpecies] >= 0) {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonCation",
|
|
"anion1 charge problem");
|
|
}
|
|
size_t jSpecies = speciesIndex(jName);
|
|
if (jSpecies == npos) {
|
|
return;
|
|
}
|
|
if (charge[jSpecies] >= 0) {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonCation",
|
|
"anion2 charge problem");
|
|
}
|
|
|
|
size_t n = iSpecies * m_kk + jSpecies;
|
|
int counter = m_CounterIJ[n];
|
|
for (size_t i = 0; i < BinSalt.nChildren(); i++) {
|
|
XML_Node& xmlChild = BinSalt.child(i);
|
|
stemp = xmlChild.name();
|
|
string nodeName = lowercase(stemp);
|
|
if (nodeName == "theta") {
|
|
stemp = xmlChild.value();
|
|
double old = m_Theta_ij[counter];
|
|
m_Theta_ij[counter] = atofCheck(stemp.c_str());
|
|
if (old != 0.0) {
|
|
if (old != m_Theta_ij[counter]) {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonCation",
|
|
"conflicting values");
|
|
}
|
|
}
|
|
}
|
|
if (nodeName == "psi") {
|
|
getFloatArray(xmlChild, vParams, false, "", stemp);
|
|
nParamsFound = vParams.size();
|
|
n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies ;
|
|
|
|
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
|
|
if (nParamsFound != 1) {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonCation::Psi for "
|
|
+ kName + "::" + iName + "::" + jName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Psi_ijk_coeff(0,n) = vParams[0];
|
|
m_Psi_ijk[n] = vParams[0];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
|
|
if (nParamsFound != 2) {
|
|
throw CanteraError("HMWSoln::readXMLPsiCation::Psi for "
|
|
+ kName + "::" + iName + "::" + jName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Psi_ijk_coeff(0,n) = vParams[0];
|
|
m_Psi_ijk_coeff(1,n) = vParams[1];
|
|
m_Psi_ijk[n] = vParams[0];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
|
|
if (nParamsFound == 1) {
|
|
vParams.resize(5, 0.0);
|
|
nParamsFound = 5;
|
|
} else if (nParamsFound != 5) {
|
|
throw CanteraError("HMWSoln::readXMLPsiCation::Psi for "
|
|
+ kName + "::" + iName + "::" + jName,
|
|
"wrong number of params found");
|
|
}
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Psi_ijk_coeff(j, n) = vParams[j];
|
|
}
|
|
m_Psi_ijk[n] = vParams[0];
|
|
}
|
|
|
|
|
|
// fill in the duplicate entries
|
|
n = iSpecies * m_kk *m_kk + kSpecies * m_kk + jSpecies ;
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Psi_ijk_coeff(j, n) = vParams[j];
|
|
}
|
|
m_Psi_ijk[n] = vParams[0];
|
|
|
|
n = jSpecies * m_kk *m_kk + iSpecies * m_kk + kSpecies ;
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Psi_ijk_coeff(j, n) = vParams[j];
|
|
}
|
|
m_Psi_ijk[n] = vParams[0];
|
|
|
|
n = jSpecies * m_kk *m_kk + kSpecies * m_kk + iSpecies ;
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Psi_ijk_coeff(j, n) = vParams[j];
|
|
}
|
|
m_Psi_ijk[n] = vParams[0];
|
|
|
|
n = kSpecies * m_kk *m_kk + jSpecies * m_kk + iSpecies ;
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Psi_ijk_coeff(j, n) = vParams[j];
|
|
}
|
|
m_Psi_ijk[n] = vParams[0];
|
|
|
|
n = kSpecies * m_kk *m_kk + iSpecies * m_kk + jSpecies ;
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Psi_ijk_coeff(j, n) = vParams[j];
|
|
}
|
|
m_Psi_ijk[n] = vParams[0];
|
|
}
|
|
}
|
|
}
|
|
|
|
/**
|
|
* Process an XML node called "PsiCommonAnion".
|
|
* This node contains all of the parameters necessary to describe
|
|
* the binary interactions between two cations and one common anion.
|
|
*/
|
|
void HMWSoln::readXMLPsiCommonAnion(XML_Node& BinSalt)
|
|
{
|
|
string xname = BinSalt.name();
|
|
if (xname != "psiCommonAnion") {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonAnion",
|
|
"Incorrect name for processing this routine: " + xname);
|
|
}
|
|
double* charge = DATA_PTR(m_speciesCharge);
|
|
string stemp;
|
|
vector_fp vParams;
|
|
size_t nParamsFound = 0;
|
|
string kName = BinSalt.attrib("anion");
|
|
if (kName == "") {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "no anion attrib");
|
|
}
|
|
string iName = BinSalt.attrib("cation1");
|
|
if (iName == "") {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "no cation1 attrib");
|
|
}
|
|
string jName = BinSalt.attrib("cation2");
|
|
if (jName == "") {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "no cation2 attrib");
|
|
}
|
|
/*
|
|
* Find the index of the species in the current phase. It's not
|
|
* an error to not find the species
|
|
*/
|
|
size_t kSpecies = speciesIndex(kName);
|
|
if (kSpecies == npos) {
|
|
return;
|
|
}
|
|
if (charge[kSpecies] >= 0) {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "anion charge problem");
|
|
}
|
|
size_t iSpecies = speciesIndex(iName);
|
|
if (iSpecies == npos) {
|
|
return;
|
|
}
|
|
if (charge[iSpecies] <= 0) {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonAnion",
|
|
"cation1 charge problem");
|
|
}
|
|
size_t jSpecies = speciesIndex(jName);
|
|
if (jSpecies == npos) {
|
|
return;
|
|
}
|
|
if (charge[jSpecies] <= 0) {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonAnion",
|
|
"cation2 charge problem");
|
|
}
|
|
|
|
size_t n = iSpecies * m_kk + jSpecies;
|
|
int counter = m_CounterIJ[n];
|
|
for (size_t i = 0; i < BinSalt.nChildren(); i++) {
|
|
XML_Node& xmlChild = BinSalt.child(i);
|
|
stemp = xmlChild.name();
|
|
string nodeName = lowercase(stemp);
|
|
if (nodeName == "theta") {
|
|
stemp = xmlChild.value();
|
|
double old = m_Theta_ij[counter];
|
|
m_Theta_ij[counter] = atofCheck(stemp.c_str());
|
|
if (old != 0.0) {
|
|
if (old != m_Theta_ij[counter]) {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonAnion",
|
|
"conflicting values");
|
|
}
|
|
}
|
|
}
|
|
if (nodeName == "psi") {
|
|
|
|
getFloatArray(xmlChild, vParams, false, "", stemp);
|
|
nParamsFound = vParams.size();
|
|
n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies ;
|
|
|
|
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
|
|
if (nParamsFound != 1) {
|
|
throw CanteraError("HMWSoln::readXMLPsiCommonAnion::Psi for "
|
|
+ kName + "::" + iName + "::" + jName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Psi_ijk_coeff(0,n) = vParams[0];
|
|
m_Psi_ijk[n] = vParams[0];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
|
|
if (nParamsFound != 2) {
|
|
throw CanteraError("HMWSoln::readXMLPsiAnion::Psi for "
|
|
+ kName + "::" + iName + "::" + jName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Psi_ijk_coeff(0,n) = vParams[0];
|
|
m_Psi_ijk_coeff(1,n) = vParams[1];
|
|
m_Psi_ijk[n] = vParams[0];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
|
|
if (nParamsFound == 1) {
|
|
vParams.resize(5, 0.0);
|
|
nParamsFound = 5;
|
|
} else if (nParamsFound != 5) {
|
|
throw CanteraError("HMWSoln::readXMLPsiAnion::Psi for "
|
|
+ kName + "::" + iName + "::" + jName,
|
|
"wrong number of params found");
|
|
}
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Psi_ijk_coeff(j, n) = vParams[j];
|
|
}
|
|
m_Psi_ijk[n] = vParams[0];
|
|
}
|
|
|
|
|
|
// fill in the duplicate entries
|
|
n = iSpecies * m_kk *m_kk + kSpecies * m_kk + jSpecies ;
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Psi_ijk_coeff(j, n) = vParams[j];
|
|
}
|
|
m_Psi_ijk[n] = vParams[0];
|
|
|
|
n = jSpecies * m_kk *m_kk + iSpecies * m_kk + kSpecies ;
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Psi_ijk_coeff(j, n) = vParams[j];
|
|
}
|
|
m_Psi_ijk[n] = vParams[0];
|
|
|
|
n = jSpecies * m_kk *m_kk + kSpecies * m_kk + iSpecies ;
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Psi_ijk_coeff(j, n) = vParams[j];
|
|
}
|
|
m_Psi_ijk[n] = vParams[0];
|
|
|
|
n = kSpecies * m_kk *m_kk + jSpecies * m_kk + iSpecies ;
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Psi_ijk_coeff(j, n) = vParams[j];
|
|
}
|
|
m_Psi_ijk[n] = vParams[0];
|
|
|
|
n = kSpecies * m_kk *m_kk + iSpecies * m_kk + jSpecies ;
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Psi_ijk_coeff(j, n) = vParams[j];
|
|
}
|
|
m_Psi_ijk[n] = vParams[0];
|
|
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
/**
|
|
* Process an XML node called "LambdaNeutral".
|
|
* This node contains all of the parameters necessary to describe
|
|
* the binary interactions between one neutral species and
|
|
* any other species (neutral or otherwise) in the mechanism.
|
|
*/
|
|
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);
|
|
}
|
|
double* charge = DATA_PTR(m_speciesCharge);
|
|
string stemp;
|
|
string iName = BinSalt.attrib("species1");
|
|
if (iName == "") {
|
|
throw CanteraError("HMWSoln::readXMLLambdaNeutral", "no species1 attrib");
|
|
}
|
|
string jName = BinSalt.attrib("species2");
|
|
if (jName == "") {
|
|
throw CanteraError("HMWSoln::readXMLLambdaNeutral", "no species2 attrib");
|
|
}
|
|
/*
|
|
* Find the index of the species in the current phase. It's not
|
|
* an error to not find the species
|
|
*/
|
|
size_t iSpecies = speciesIndex(iName);
|
|
if (iSpecies == npos) {
|
|
return;
|
|
}
|
|
if (charge[iSpecies] != 0) {
|
|
throw CanteraError("HMWSoln::readXMLLambdaNeutral",
|
|
"neutral charge problem");
|
|
}
|
|
size_t jSpecies = speciesIndex(jName);
|
|
if (jSpecies == npos) {
|
|
return;
|
|
}
|
|
|
|
for (size_t i = 0; i < BinSalt.nChildren(); i++) {
|
|
XML_Node& xmlChild = BinSalt.child(i);
|
|
stemp = xmlChild.name();
|
|
string nodeName = lowercase(stemp);
|
|
if (nodeName == "lambda") {
|
|
size_t nCount = iSpecies*m_kk + jSpecies;
|
|
getFloatArray(xmlChild, vParams, false, "", stemp);
|
|
nParamsFound = vParams.size();
|
|
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
|
|
if (nParamsFound != 1) {
|
|
throw CanteraError("HMWSoln::readXMLLambdaNeutral::Lambda for " + iName
|
|
+ "::" + jName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Lambda_nj_coeff(0,nCount) = vParams[0];
|
|
m_Lambda_nj(iSpecies,jSpecies) = vParams[0];
|
|
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
|
|
if (nParamsFound != 2) {
|
|
throw CanteraError("HMWSoln::readXMLLambdaNeutral::Lambda for " + iName
|
|
+ "::" + jName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Lambda_nj_coeff(0,nCount) = vParams[0];
|
|
m_Lambda_nj_coeff(1,nCount) = vParams[1];
|
|
m_Lambda_nj(iSpecies, jSpecies) = vParams[0];
|
|
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
|
|
if (nParamsFound == 1) {
|
|
vParams.resize(5, 0.0);
|
|
nParamsFound = 5;
|
|
} else if (nParamsFound != 5) {
|
|
throw CanteraError("HMWSoln::readXMLLambdaNeutral::Lambda for " + iName
|
|
+ "::" + jName,
|
|
"wrong number of params found");
|
|
}
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Lambda_nj_coeff(j,nCount) = vParams[j];
|
|
}
|
|
m_Lambda_nj(iSpecies, jSpecies) = vParams[0];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/**
|
|
* Process an XML node called "MunnnNeutral".
|
|
* This node contains all of the parameters necessary to describe
|
|
* the self-ternary interactions for one neutral species.
|
|
*/
|
|
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);
|
|
}
|
|
double* charge = DATA_PTR(m_speciesCharge);
|
|
string stemp;
|
|
string iName = BinSalt.attrib("species1");
|
|
if (iName == "") {
|
|
throw CanteraError("HMWSoln::readXMLMunnnNeutral", "no species1 attrib");
|
|
}
|
|
|
|
/*
|
|
* Find the index of the species in the current phase. It's not
|
|
* an error to not find the species
|
|
*/
|
|
size_t iSpecies = speciesIndex(iName);
|
|
if (iSpecies == npos) {
|
|
return;
|
|
}
|
|
if (charge[iSpecies] != 0) {
|
|
throw CanteraError("HMWSoln::readXMLMunnnNeutral",
|
|
"neutral charge problem");
|
|
}
|
|
|
|
for (size_t i = 0; i < BinSalt.nChildren(); i++) {
|
|
XML_Node& xmlChild = BinSalt.child(i);
|
|
stemp = xmlChild.name();
|
|
string nodeName = lowercase(stemp);
|
|
if (nodeName == "munnn") {
|
|
getFloatArray(xmlChild, vParams, false, "", "Munnn");
|
|
nParamsFound = vParams.size();
|
|
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
|
|
if (nParamsFound != 1) {
|
|
throw CanteraError("HMWSoln::readXMLMunnnNeutral::Munnn for " + iName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Mu_nnn_coeff(0,iSpecies) = vParams[0];
|
|
m_Mu_nnn[iSpecies] = vParams[0];
|
|
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
|
|
if (nParamsFound != 2) {
|
|
throw CanteraError("HMWSoln::readXMLMunnnNeutral::Munnn for " + iName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Mu_nnn_coeff(0, iSpecies) = vParams[0];
|
|
m_Mu_nnn_coeff(1, iSpecies) = vParams[1];
|
|
m_Mu_nnn[iSpecies] = vParams[0];
|
|
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
|
|
if (nParamsFound == 1) {
|
|
vParams.resize(5, 0.0);
|
|
nParamsFound = 5;
|
|
} else if (nParamsFound != 5) {
|
|
throw CanteraError("HMWSoln::readXMLMunnnNeutral::Munnn for " + iName,
|
|
"wrong number of params found");
|
|
}
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Mu_nnn_coeff(j, iSpecies) = vParams[j];
|
|
}
|
|
m_Mu_nnn[iSpecies] = vParams[0];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Process an XML node called "readXMLZetaCation".
|
|
* This node contains all of the parameters necessary to describe
|
|
* the ternary interactions between a neutral, a cation and an anion
|
|
*/
|
|
void HMWSoln::readXMLZetaCation(const XML_Node& BinSalt)
|
|
{
|
|
string xname = BinSalt.name();
|
|
if (xname != "zetaCation") {
|
|
throw CanteraError("HMWSoln::readXMLZetaCation",
|
|
"Incorrect name for processing this routine: " + xname);
|
|
}
|
|
double* charge = DATA_PTR(m_speciesCharge);
|
|
string stemp;
|
|
vector_fp vParams;
|
|
size_t nParamsFound = 0;
|
|
|
|
string iName = BinSalt.attrib("neutral");
|
|
if (iName == "") {
|
|
throw CanteraError("HMWSoln::readXMLZetaCation", "no neutral attrib");
|
|
}
|
|
|
|
string jName = BinSalt.attrib("cation1");
|
|
if (jName == "") {
|
|
throw CanteraError("HMWSoln::readXMLZetaCation", "no cation1 attrib");
|
|
}
|
|
|
|
string kName = BinSalt.attrib("anion1");
|
|
if (kName == "") {
|
|
throw CanteraError("HMWSoln::readXMLZetaCation", "no anion1 attrib");
|
|
}
|
|
/*
|
|
* Find the index of the species in the current phase. It's not
|
|
* an error to not find the species
|
|
*/
|
|
size_t iSpecies = speciesIndex(iName);
|
|
if (iSpecies == npos) {
|
|
return;
|
|
}
|
|
if (charge[iSpecies] != 0.0) {
|
|
throw CanteraError("HMWSoln::readXMLZetaCation", "neutral charge problem");
|
|
}
|
|
|
|
size_t jSpecies = speciesIndex(jName);
|
|
if (jSpecies == npos) {
|
|
return;
|
|
}
|
|
if (charge[jSpecies] <= 0.0) {
|
|
throw CanteraError("HMWSoln::readXLZetaCation", "cation1 charge problem");
|
|
}
|
|
|
|
size_t kSpecies = speciesIndex(kName);
|
|
if (kSpecies == npos) {
|
|
return;
|
|
}
|
|
if (charge[kSpecies] >= 0.0) {
|
|
throw CanteraError("HMWSoln::readXMLZetaCation", "anion1 charge problem");
|
|
}
|
|
|
|
for (size_t i = 0; i < BinSalt.nChildren(); i++) {
|
|
XML_Node& xmlChild = BinSalt.child(i);
|
|
stemp = xmlChild.name();
|
|
string nodeName = lowercase(stemp);
|
|
if (nodeName == "zeta") {
|
|
getFloatArray(xmlChild, vParams, false, "", "zeta");
|
|
nParamsFound = vParams.size();
|
|
size_t n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies ;
|
|
|
|
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
|
|
if (nParamsFound != 1) {
|
|
throw CanteraError("HMWSoln::readXMLZetaCation::Zeta for "
|
|
+ iName + "::" + jName + "::" + kName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Psi_ijk_coeff(0,n) = vParams[0];
|
|
m_Psi_ijk[n] = vParams[0];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
|
|
if (nParamsFound != 2) {
|
|
throw CanteraError("HMWSoln::readXMLZetaCation::Zeta for "
|
|
+ iName + "::" + jName + "::" + kName,
|
|
"wrong number of params found");
|
|
}
|
|
m_Psi_ijk_coeff(0,n) = vParams[0];
|
|
m_Psi_ijk_coeff(1,n) = vParams[1];
|
|
m_Psi_ijk[n] = vParams[0];
|
|
} else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
|
|
if (nParamsFound == 1) {
|
|
vParams.resize(5, 0.0);
|
|
nParamsFound = 5;
|
|
} else if (nParamsFound != 5) {
|
|
throw CanteraError("HMWSoln::readXMLZetaCation::Zeta for "
|
|
+ iName + "::" + jName + "::" + kName,
|
|
"wrong number of params found");
|
|
}
|
|
for (size_t j = 0; j < nParamsFound; j++) {
|
|
m_Psi_ijk_coeff(j, n) = vParams[j];
|
|
}
|
|
m_Psi_ijk[n] = vParams[0];
|
|
}
|
|
|
|
// There are no duplicate entries
|
|
}
|
|
}
|
|
}
|
|
|
|
// Process an XML node called "croppingCoefficients"
|
|
// for the cropping coefficients values
|
|
/*
|
|
* @param acNode Activity Coefficient XML Node
|
|
*/
|
|
void HMWSoln::readXMLCroppingCoefficients(const XML_Node& acNode)
|
|
{
|
|
|
|
if (acNode.hasChild("croppingCoefficients")) {
|
|
XML_Node& cropNode = acNode.child("croppingCoefficients");
|
|
if (cropNode.hasChild("ln_gamma_k_min")) {
|
|
XML_Node& gkminNode = cropNode.child("ln_gamma_k_min");
|
|
getOptionalFloat(gkminNode, "pureSolventValue", CROP_ln_gamma_k_min);
|
|
}
|
|
if (cropNode.hasChild("ln_gamma_k_max")) {
|
|
XML_Node& gkmaxNode = cropNode.child("ln_gamma_k_max");
|
|
getOptionalFloat(gkmaxNode, "pureSolventValue", CROP_ln_gamma_k_max);
|
|
}
|
|
|
|
if (cropNode.hasChild("ln_gamma_o_min")) {
|
|
XML_Node& gominNode = cropNode.child("ln_gamma_o_min");
|
|
getOptionalFloat(gominNode, "pureSolventValue", CROP_ln_gamma_o_min);
|
|
}
|
|
|
|
if (cropNode.hasChild("ln_gamma_o_max")) {
|
|
XML_Node& gomaxNode = cropNode.child("ln_gamma_o_max");
|
|
getOptionalFloat(gomaxNode, "pureSolventValue", CROP_ln_gamma_o_max);
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Initialization routine for a HMWSoln phase.
|
|
*
|
|
* This is a virtual routine. This routine will call initThermo()
|
|
* for the parent class as well.
|
|
*/
|
|
void HMWSoln::initThermo()
|
|
{
|
|
MolalityVPSSTP::initThermo();
|
|
initLengths();
|
|
}
|
|
|
|
/*
|
|
* Import, construct, and initialize a HMWSoln phase
|
|
* specification from an XML tree into the current object.
|
|
*
|
|
* This routine is a precursor to constructPhaseXML(XML_Node*)
|
|
* routine, which does most of the work.
|
|
*
|
|
* @param infile XML file containing the description of the
|
|
* phase
|
|
*
|
|
* @param id Optional parameter identifying the name of the
|
|
* phase. If none is given, the first XML
|
|
* phase element will be used.
|
|
*/
|
|
void HMWSoln::constructPhaseFile(std::string inputFile, std::string id) {
|
|
|
|
if (inputFile.size() == 0) {
|
|
throw CanteraError("HMWSoln:constructPhaseFile",
|
|
"input file is null");
|
|
}
|
|
string path = findInputFile(inputFile);
|
|
std::ifstream fin(path.c_str());
|
|
if (!fin) {
|
|
throw CanteraError("HMWSoln:constructPhaseFile","could not open "
|
|
+path+" for reading.");
|
|
}
|
|
/*
|
|
* The phase object automatically constructs an XML object.
|
|
* Use this object to store information.
|
|
*/
|
|
XML_Node &phaseNode_XML = xml();
|
|
XML_Node *fxml = new XML_Node();
|
|
fxml->build(fin);
|
|
XML_Node *fxml_phase = findXMLPhase(fxml, id);
|
|
if (!fxml_phase) {
|
|
throw CanteraError("HMWSoln:constructPhaseFile",
|
|
"ERROR: Can not find phase named " +
|
|
id + " in file named " + inputFile);
|
|
}
|
|
fxml_phase->copy(&phaseNode_XML);
|
|
constructPhaseXML(*fxml_phase, id);
|
|
delete fxml;
|
|
}
|
|
|
|
/*
|
|
* Import, construct, and initialize a HMWSoln phase
|
|
* specification from an XML tree into the current object.
|
|
*
|
|
* Most of the work is carried out by the cantera base
|
|
* routine, importPhase(). That routine imports all of the
|
|
* species and element data, including the standard states
|
|
* of the species.
|
|
*
|
|
* Then, In this routine, we read the information
|
|
* particular to the specification of the activity
|
|
* coefficient model for the Pitzer parameterization.
|
|
*
|
|
* We also read information about the molar volumes of the
|
|
* standard states if present in the XML file.
|
|
*
|
|
* @param phaseNode This object must be the phase node of a
|
|
* complete XML tree
|
|
* description of the phase, including all of the
|
|
* species data. In other words while "phase" must
|
|
* point to an XML phase object, it must have
|
|
* sibling nodes "speciesData" that describe
|
|
* the species in the phase.
|
|
* @param id ID of the phase. If nonnull, a check is done
|
|
* to see if phaseNode is pointing to the phase
|
|
* with the correct id.
|
|
*/
|
|
void HMWSoln::constructPhaseXML(XML_Node& phaseNode, std::string id) {
|
|
string stemp;
|
|
if (id.size() > 0) {
|
|
string idp = phaseNode.id();
|
|
if (idp != id) {
|
|
throw CanteraError("HMWSoln::constructPhaseXML",
|
|
"phasenode and Id are incompatible");
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Find the Thermo XML node
|
|
*/
|
|
if (!phaseNode.hasChild("thermo")) {
|
|
throw CanteraError("HMWSoln::constructPhaseXML",
|
|
"no thermo XML node");
|
|
}
|
|
XML_Node& thermoNode = phaseNode.child("thermo");
|
|
|
|
/*
|
|
* Possibly change the form of the standard concentrations
|
|
*/
|
|
if (thermoNode.hasChild("standardConc")) {
|
|
XML_Node& scNode = thermoNode.child("standardConc");
|
|
m_formGC = 2;
|
|
stemp = scNode.attrib("model");
|
|
string formString = lowercase(stemp);
|
|
if (formString != "") {
|
|
if (formString == "unity") {
|
|
m_formGC = 0;
|
|
printf("exit standardConc = unity not done\n");
|
|
exit(EXIT_FAILURE);
|
|
} else if (formString == "molar_volume") {
|
|
m_formGC = 1;
|
|
printf("exit standardConc = molar_volume not done\n");
|
|
exit(EXIT_FAILURE);
|
|
} else if (formString == "solvent_volume") {
|
|
m_formGC = 2;
|
|
} else {
|
|
throw CanteraError("HMWSoln::constructPhaseXML",
|
|
"Unknown standardConc model: " + formString);
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
* Get the Name of the Solvent:
|
|
* <solvent> solventName </solvent>
|
|
*/
|
|
string solventName = "";
|
|
if (thermoNode.hasChild("solvent")) {
|
|
XML_Node& scNode = thermoNode.child("solvent");
|
|
vector<string> nameSolventa;
|
|
getStringArray(scNode, nameSolventa);
|
|
int nsp = static_cast<int>(nameSolventa.size());
|
|
if (nsp != 1) {
|
|
throw CanteraError("HMWSoln::constructPhaseXML",
|
|
"badly formed solvent XML node");
|
|
}
|
|
solventName = nameSolventa[0];
|
|
}
|
|
|
|
/*
|
|
* Determine the form of the Pitzer model,
|
|
* We will use this information to size arrays below.
|
|
*/
|
|
if (thermoNode.hasChild("activityCoefficients")) {
|
|
XML_Node& scNode = thermoNode.child("activityCoefficients");
|
|
stemp = scNode.attrib("model");
|
|
string formString = lowercase(stemp);
|
|
if (formString != "") {
|
|
if (formString == "pitzer" || formString == "default") {
|
|
m_formPitzer = PITZERFORM_BASE;
|
|
} else if (formString == "base") {
|
|
m_formPitzer = PITZERFORM_BASE;
|
|
} else {
|
|
throw CanteraError("HMWSoln::constructPhaseXML",
|
|
"Unknown Pitzer ActivityCoeff model: "
|
|
+ formString);
|
|
}
|
|
}
|
|
/*
|
|
* Determine the form of the temperature dependence
|
|
* of the Pitzer activity coefficient model.
|
|
*/
|
|
stemp = scNode.attrib("TempModel");
|
|
formString = lowercase(stemp);
|
|
if (formString != "") {
|
|
if (formString == "constant" || formString == "default") {
|
|
m_formPitzerTemp = PITZER_TEMP_CONSTANT;
|
|
} else if (formString == "linear") {
|
|
m_formPitzerTemp = PITZER_TEMP_LINEAR;
|
|
} else if (formString == "complex" || formString == "complex1") {
|
|
m_formPitzerTemp = PITZER_TEMP_COMPLEX1;
|
|
} else {
|
|
throw CanteraError("HMWSoln::constructPhaseXML",
|
|
"Unknown Pitzer ActivityCoeff Temp model: "
|
|
+ formString);
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Determine the reference temperature
|
|
* of the Pitzer activity coefficient model's temperature
|
|
* dependence formulation: defaults to 25C
|
|
*/
|
|
stemp = scNode.attrib("TempReference");
|
|
formString = lowercase(stemp);
|
|
if (formString != "") {
|
|
m_TempPitzerRef = atofCheck(formString.c_str());
|
|
} else {
|
|
m_TempPitzerRef = 273.15 + 25;
|
|
}
|
|
|
|
}
|
|
|
|
/*
|
|
* Call the Cantera importPhase() function. This will import
|
|
* all of the species into the phase. This will also handle
|
|
* all of the solvent and solute standard states
|
|
*/
|
|
bool m_ok = importPhase(phaseNode, this);
|
|
if (!m_ok) {
|
|
throw CanteraError("HMWSoln::constructPhaseXML","importPhase failed ");
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
/**
|
|
* Process the XML file after species are set up.
|
|
*
|
|
* This gets called from importPhase(). It processes the XML file
|
|
* after the species are set up. This is the main routine for
|
|
* reading in activity coefficient parameters.
|
|
*
|
|
* @param phaseNode This object must be the phase node of a
|
|
* complete XML tree
|
|
* description of the phase, including all of the
|
|
* species data. In other words while "phase" must
|
|
* point to an XML phase object, it must have
|
|
* sibling nodes "speciesData" that describe
|
|
* the species in the phase.
|
|
* @param id ID of the phase. If nonnull, a check is done
|
|
* to see if phaseNode is pointing to the phase
|
|
* with the correct id.
|
|
*/
|
|
void HMWSoln::
|
|
initThermoXML(XML_Node& phaseNode, const std::string& id)
|
|
{
|
|
string stemp;
|
|
if (id.size() > 0) {
|
|
string idp = phaseNode.id();
|
|
if (idp != id) {
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"phasenode and Id are incompatible");
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Find the Thermo XML node
|
|
*/
|
|
if (!phaseNode.hasChild("thermo")) {
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"no thermo XML node");
|
|
}
|
|
XML_Node& thermoNode = phaseNode.child("thermo");
|
|
|
|
/*
|
|
* Possibly change the form of the standard concentrations
|
|
*/
|
|
if (thermoNode.hasChild("standardConc")) {
|
|
XML_Node& scNode = thermoNode.child("standardConc");
|
|
m_formGC = 2;
|
|
stemp = scNode.attrib("model");
|
|
string formString = lowercase(stemp);
|
|
if (formString != "") {
|
|
if (formString == "unity") {
|
|
m_formGC = 0;
|
|
printf("exit standardConc = unity not done\n");
|
|
exit(EXIT_FAILURE);
|
|
} else if (formString == "molar_volume") {
|
|
m_formGC = 1;
|
|
printf("exit standardConc = molar_volume not done\n");
|
|
exit(EXIT_FAILURE);
|
|
} else if (formString == "solvent_volume") {
|
|
m_formGC = 2;
|
|
} else {
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"Unknown standardConc model: " + formString);
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Determine the form of the Pitzer model,
|
|
* We will use this information to size arrays below.
|
|
*/
|
|
if (thermoNode.hasChild("activityCoefficients")) {
|
|
XML_Node& scNode = thermoNode.child("activityCoefficients");
|
|
stemp = scNode.attrib("model");
|
|
string formString = lowercase(stemp);
|
|
if (formString != "") {
|
|
if (formString == "pitzer" || formString == "default") {
|
|
m_formPitzer = PITZERFORM_BASE;
|
|
} else if (formString == "base") {
|
|
m_formPitzer = PITZERFORM_BASE;
|
|
} else {
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"Unknown Pitzer ActivityCoeff model: "
|
|
+ formString);
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Determine the form of the temperature dependence
|
|
* of the Pitzer activity coefficient model.
|
|
*/
|
|
stemp = scNode.attrib("TempModel");
|
|
formString = lowercase(stemp);
|
|
if (formString != "") {
|
|
if (formString == "constant" || formString == "default") {
|
|
m_formPitzerTemp = PITZER_TEMP_CONSTANT;
|
|
} else if (formString == "linear") {
|
|
m_formPitzerTemp = PITZER_TEMP_LINEAR;
|
|
} else if (formString == "complex" || formString == "complex1") {
|
|
m_formPitzerTemp = PITZER_TEMP_COMPLEX1;
|
|
} else {
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"Unknown Pitzer ActivityCoeff Temp model: "
|
|
+ formString);
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Determine the reference temperature
|
|
* of the Pitzer activity coefficient model's temperature
|
|
* dependence formulation: defaults to 25C
|
|
*/
|
|
stemp = scNode.attrib("TempReference");
|
|
formString = lowercase(stemp);
|
|
if (formString != "") {
|
|
m_TempPitzerRef = atofCheck(formString.c_str());
|
|
} else {
|
|
m_TempPitzerRef = 273.15 + 25;
|
|
}
|
|
|
|
}
|
|
|
|
/*
|
|
* Get the Name of the Solvent:
|
|
* <solvent> solventName </solvent>
|
|
*/
|
|
string solventName = "";
|
|
if (thermoNode.hasChild("solvent")) {
|
|
XML_Node& scNode = thermoNode.child("solvent");
|
|
vector<string> nameSolventa;
|
|
getStringArray(scNode, nameSolventa);
|
|
int nsp = static_cast<int>(nameSolventa.size());
|
|
if (nsp != 1) {
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"badly formed solvent XML node");
|
|
}
|
|
solventName = nameSolventa[0];
|
|
}
|
|
|
|
/*
|
|
* Initialize all of the lengths of arrays in the object
|
|
* now that we know what species are in the phase.
|
|
*/
|
|
initLengths();
|
|
|
|
/*
|
|
* Reconcile the solvent name and index.
|
|
*/
|
|
for (size_t k = 0; k < m_kk; k++) {
|
|
string sname = speciesName(k);
|
|
if (solventName == sname) {
|
|
setSolvent(k);
|
|
if (k != 0) {
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"Solvent must be species 0 atm");
|
|
}
|
|
m_indexSolvent = k;
|
|
break;
|
|
}
|
|
}
|
|
if (m_indexSolvent == npos) {
|
|
std::cout << "HMWSoln::initThermo: Solvent Name not found"
|
|
<< std::endl;
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"Solvent name not found");
|
|
}
|
|
if (m_indexSolvent != 0) {
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"Solvent " + solventName +
|
|
" should be first species");
|
|
}
|
|
|
|
/*
|
|
* Now go get the specification of the standard states for
|
|
* species in the solution. This includes the molar volumes
|
|
* data blocks for incompressible species.
|
|
*/
|
|
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]);
|
|
if (!s) {
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"Species Data Base " + sss[k] + " not found");
|
|
}
|
|
XML_Node* ss = s->findByName("standardState");
|
|
if (!ss) {
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"Species " + sss[k] +
|
|
" standardState XML block not found");
|
|
}
|
|
string modelStringa = ss->attrib("model");
|
|
if (modelStringa == "") {
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"Species " + sss[k] +
|
|
" standardState XML block model attribute not found");
|
|
}
|
|
string modelString = lowercase(modelStringa);
|
|
if (k == 0) {
|
|
if (modelString == "wateriapws" || modelString == "real_water" ||
|
|
modelString == "waterpdss") {
|
|
/*
|
|
* Store a local pointer to the water standard state model.
|
|
* -> We've hardcoded it to a PDSS_Water model, so this is ok.
|
|
*/
|
|
m_waterSS = dynamic_cast<PDSS_Water*>(providePDSS(0)) ;
|
|
if (!m_waterSS) {
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"Dynamic cast to PDSS_Water failed");
|
|
}
|
|
/*
|
|
* Fill in the molar volume of water (m3/kmol)
|
|
* at standard conditions to fill in the m_speciesSize entry
|
|
* with something reasonable.
|
|
*/
|
|
m_waterSS->setState_TP(300., OneAtm);
|
|
double dens = m_waterSS->density();
|
|
double mw = m_waterSS->molecularWeight();
|
|
m_speciesSize[0] = mw / dens;
|
|
#ifdef DEBUG_HKM_NOT
|
|
cout << "Solvent species " << sss[k] << " has volume " <<
|
|
m_speciesSize[k] << endl;
|
|
#endif
|
|
} else {
|
|
// throw CanteraError("HMWSoln::initThermoXML",
|
|
// "Solvent SS Model \"" + modelStringa +
|
|
// "\" is not allowed, name = " + sss[0]);
|
|
m_waterSS = providePDSS(0);
|
|
m_waterSS->setState_TP(300., OneAtm);
|
|
double dens = m_waterSS->density();
|
|
double mw = m_waterSS->molecularWeight();
|
|
m_speciesSize[0] = mw / dens;
|
|
}
|
|
} else {
|
|
if (modelString != "constant_incompressible" && modelString != "hkft") {
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"Solute SS Model \"" + modelStringa +
|
|
"\" is not known");
|
|
}
|
|
if (modelString == "constant_incompressible") {
|
|
m_speciesSize[k] = getFloat(*ss, "molarVolume", "toSI");
|
|
#ifdef DEBUG_HKM_NOT
|
|
cout << "species " << sss[k] << " has volume " <<
|
|
m_speciesSize[k] << endl;
|
|
#endif
|
|
}
|
|
// HKM Note, have to fill up m_speciesSize[] for HKFT species
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Initialize the water property calculator. It will share
|
|
* the internal eos water calculator.
|
|
*/
|
|
m_waterProps = new WaterProps(dynamic_cast<PDSS_Water*>(m_waterSS));
|
|
|
|
/*
|
|
* Fill in parameters for the calculation of the
|
|
* stoichiometric Ionic Strength
|
|
*
|
|
* The default is that stoich charge is the same as the
|
|
* regular charge.
|
|
*/
|
|
for (size_t k = 0; k < m_kk; k++) {
|
|
m_speciesCharge_Stoich[k] = m_speciesCharge[k];
|
|
}
|
|
|
|
/*
|
|
* 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;
|
|
/*
|
|
* Look for parameters for A_Debye
|
|
*/
|
|
if (acNode.hasChild("A_Debye")) {
|
|
XML_Node& ADebye = acNode.child("A_Debye");
|
|
m_form_A_Debye = A_DEBYE_CONST;
|
|
stemp = "model";
|
|
if (ADebye.hasAttrib(stemp)) {
|
|
string atemp = ADebye.attrib(stemp);
|
|
stemp = lowercase(atemp);
|
|
if (stemp == "water") {
|
|
m_form_A_Debye = A_DEBYE_WATER;
|
|
}
|
|
}
|
|
if (m_form_A_Debye == A_DEBYE_CONST) {
|
|
m_A_Debye = getFloat(acNode, "A_Debye");
|
|
}
|
|
#ifdef DEBUG_HKM_NOT
|
|
cout << "A_Debye = " << m_A_Debye << endl;
|
|
#endif
|
|
}
|
|
|
|
/*
|
|
* Look for Parameters for the Maximum Ionic Strength
|
|
*/
|
|
if (acNode.hasChild("maxIonicStrength")) {
|
|
m_maxIionicStrength = getFloat(acNode, "maxIonicStrength");
|
|
#ifdef DEBUG_HKM_NOT
|
|
cout << "m_maxIionicStrength = "
|
|
<<m_maxIionicStrength << endl;
|
|
#endif
|
|
}
|
|
|
|
|
|
/*
|
|
* Look for parameters for the Ionic radius
|
|
*/
|
|
if (acNode.hasChild("ionicRadius")) {
|
|
XML_Node& irNode = acNode.child("ionicRadius");
|
|
|
|
string Aunits = "";
|
|
double Afactor = 1.0;
|
|
if (irNode.hasAttrib("units")) {
|
|
string Aunits = irNode.attrib("units");
|
|
Afactor = toSI(Aunits);
|
|
}
|
|
|
|
if (irNode.hasAttrib("default")) {
|
|
string ads = irNode.attrib("default");
|
|
double ad = fpValue(ads);
|
|
for (size_t k = 0; k < m_kk; k++) {
|
|
m_Aionic[k] = ad * Afactor;
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
/*
|
|
* First look at the species database.
|
|
* -> Look for the subelement "stoichIsMods"
|
|
* in each of the species SS databases.
|
|
*/
|
|
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++) {
|
|
const XML_Node& sp = *xspecies[j];
|
|
jname = sp["name"];
|
|
if (jname == kname) {
|
|
jmap = j;
|
|
break;
|
|
}
|
|
}
|
|
if (jmap != npos) {
|
|
const XML_Node& sp = *xspecies[jmap];
|
|
getOptionalFloat(sp, "stoichIsMods", m_speciesCharge_Stoich[k]);
|
|
// if (sp.hasChild("stoichIsMods")) {
|
|
// double val = getFloat(sp, "stoichIsMods");
|
|
//m_speciesCharge_Stoich[k] = val;
|
|
//}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Now look at the activity coefficient database
|
|
*/
|
|
if (acNodePtr) {
|
|
if (acNodePtr->hasChild("stoichIsMods")) {
|
|
XML_Node& sIsNode = acNodePtr->child("stoichIsMods");
|
|
|
|
map<string, string> msIs;
|
|
getMap(sIsNode, msIs);
|
|
map<string,string>::const_iterator _b = msIs.begin();
|
|
for (; _b != msIs.end(); ++_b) {
|
|
size_t kk = speciesIndex(_b->first);
|
|
if (kk != npos) {
|
|
double val = fpValue(_b->second);
|
|
m_speciesCharge_Stoich[kk] = val;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
/*
|
|
* Loop through the children getting multiple instances of
|
|
* parameters
|
|
*/
|
|
if (acNodePtr) {
|
|
for (size_t i = 0; i < acNodePtr->nChildren(); i++) {
|
|
XML_Node& xmlACChild = acNodePtr->child(i);
|
|
stemp = xmlACChild.name();
|
|
string nodeName = lowercase(stemp);
|
|
/*
|
|
* 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 == "binarysaltparameters") {
|
|
readXMLBinarySalt(xmlACChild);
|
|
} else if (nodeName == "thetaanion") {
|
|
readXMLThetaAnion(xmlACChild);
|
|
} else if (nodeName == "thetacation") {
|
|
readXMLThetaCation(xmlACChild);
|
|
} else if (nodeName == "psicommonanion") {
|
|
readXMLPsiCommonAnion(xmlACChild);
|
|
} else if (nodeName == "psicommoncation") {
|
|
readXMLPsiCommonCation(xmlACChild);
|
|
} else if (nodeName == "lambdaneutral") {
|
|
readXMLLambdaNeutral(xmlACChild);
|
|
} else if (nodeName == "zetacation") {
|
|
readXMLZetaCation(xmlACChild);
|
|
}
|
|
}
|
|
}
|
|
|
|
// Go look up the optional Cropping parameters
|
|
readXMLCroppingCoefficients(acNode);
|
|
|
|
}
|
|
|
|
/*
|
|
* Fill in the vector specifying the electrolyte species
|
|
* type
|
|
*
|
|
* First fill in default values. Everything is either
|
|
* a charge species, a nonpolar neutral, or the solvent.
|
|
*/
|
|
for (size_t k = 0; k < m_kk; k++) {
|
|
if (fabs(m_speciesCharge[k]) > 0.0001) {
|
|
m_electrolyteSpeciesType[k] = cEST_chargedSpecies;
|
|
if (fabs(m_speciesCharge_Stoich[k] - m_speciesCharge[k])
|
|
> 0.0001) {
|
|
m_electrolyteSpeciesType[k] = cEST_weakAcidAssociated;
|
|
}
|
|
} else if (fabs(m_speciesCharge_Stoich[k]) > 0.0001) {
|
|
m_electrolyteSpeciesType[k] = cEST_weakAcidAssociated;
|
|
} else {
|
|
m_electrolyteSpeciesType[k] = cEST_nonpolarNeutral;
|
|
}
|
|
}
|
|
m_electrolyteSpeciesType[m_indexSolvent] = cEST_solvent;
|
|
/*
|
|
* First look at the species database.
|
|
* -> Look for the subelement "stoichIsMods"
|
|
* 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];
|
|
if (!spPtr) {
|
|
if (spPtr->hasChild("electrolyteSpeciesType")) {
|
|
string est = getChildValue(*spPtr, "electrolyteSpeciesType");
|
|
if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) {
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"Bad electrolyte type: " + est);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
* Then look at the phase thermo specification
|
|
*/
|
|
if (acNodePtr) {
|
|
if (acNodePtr->hasChild("electrolyteSpeciesType")) {
|
|
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) {
|
|
size_t kk = speciesIndex(_b->first);
|
|
if (kk != npos) {
|
|
string est = _b->second;
|
|
if ((m_electrolyteSpeciesType[kk] = interp_est(est)) == -1) {
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"Bad electrolyte type: " + est);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
IMS_typeCutoff_ = 2;
|
|
if (IMS_typeCutoff_ == 2) {
|
|
calcIMSCutoffParams_();
|
|
}
|
|
calcMCCutoffParams_();
|
|
setMoleFSolventMin(1.0E-5);
|
|
|
|
MolalityVPSSTP::initThermoXML(phaseNode, id);
|
|
/*
|
|
* Lastly calculate the charge balance and then add stuff until the charges compensate
|
|
*/
|
|
|
|
vector_fp mf(m_kk, 0.0);
|
|
getMoleFractions(DATA_PTR(mf));
|
|
bool notDone = true;
|
|
|
|
do {
|
|
double sum = 0.0;
|
|
size_t kMaxC = npos;
|
|
double MaxC = 0.0;
|
|
for (size_t k = 0; k < m_kk; k++) {
|
|
sum += mf[k] * m_speciesCharge[k];
|
|
if (fabs(mf[k] * m_speciesCharge[k]) > MaxC) {
|
|
kMaxC = k;
|
|
}
|
|
}
|
|
size_t kHp = speciesIndex("H+");
|
|
size_t kOHm = speciesIndex("OH-");
|
|
|
|
|
|
if (fabs(sum) > 1.0E-30) {
|
|
if (kHp != npos) {
|
|
if (mf[kHp] > sum * 1.1) {
|
|
mf[kHp] -= sum;
|
|
mf[0] += sum;
|
|
notDone = false;
|
|
} else {
|
|
if (sum > 0.0) {
|
|
mf[kHp] *= 0.5;
|
|
mf[0] += mf[kHp];
|
|
sum -= mf[kHp];
|
|
}
|
|
}
|
|
}
|
|
if (notDone) {
|
|
if (kOHm != npos) {
|
|
if (mf[kOHm] > -sum * 1.1) {
|
|
mf[kOHm] += sum;
|
|
mf[0] -= sum;
|
|
notDone = false;
|
|
} else {
|
|
if (sum < 0.0) {
|
|
mf[kOHm] *= 0.5;
|
|
mf[0] += mf[kOHm];
|
|
sum += mf[kOHm];
|
|
}
|
|
}
|
|
}
|
|
if (notDone) {
|
|
if (kMaxC != npos) {
|
|
if (mf[kMaxC] > (1.1 * sum / m_speciesCharge[kMaxC])) {
|
|
mf[kMaxC] -= sum / m_speciesCharge[kMaxC];
|
|
mf[0] += sum / m_speciesCharge[kMaxC];
|
|
} else {
|
|
mf[kMaxC] *= 0.5;
|
|
mf[0] += mf[kMaxC];
|
|
notDone = true;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
setMoleFractions(DATA_PTR(mf));
|
|
} else {
|
|
notDone = false;
|
|
}
|
|
} while (notDone);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
// if (phaseNode.hasChild("state")) {
|
|
// XML_Node& stateNode = phaseNode.child("state");
|
|
// setStateFromXML(stateNode);
|
|
//}
|
|
|
|
}
|
|
//====================================================================================================================
|
|
// Precalculate the IMS Cutoff parameters for typeCutoff = 2
|
|
void HMWSoln::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_;
|
|
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_)
|
|
/
|
|
(IMS_X_o_cutoff_*IMS_X_o_cutoff_/IMS_cCut_ - 2.0 * IMS_X_o_cutoff_));
|
|
double tmp = IMS_afCut_ + IMS_X_o_cutoff_*(IMS_bfCut_ + IMS_dfCut_ *IMS_X_o_cutoff_);
|
|
double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
|
|
IMS_efCut_ = - eterm * (tmp);
|
|
if (fabs(IMS_efCut_ - oldV) < 1.0E-14) {
|
|
converged = true;
|
|
}
|
|
}
|
|
if (!converged) {
|
|
throw CanteraError("HMWSoln::calcIMSCutoffParams_()",
|
|
" failed to converge on the f polynomial");
|
|
}
|
|
converged = false;
|
|
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_;
|
|
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;
|
|
IMS_dgCut_ = ((- IMS_agCut_/IMS_cCut_ + IMS_bgCut_ - IMS_bgCut_*IMS_X_o_cutoff_/IMS_cCut_)
|
|
/
|
|
(IMS_X_o_cutoff_*IMS_X_o_cutoff_/IMS_cCut_ - 2.0 * IMS_X_o_cutoff_));
|
|
double tmp = IMS_agCut_ + IMS_X_o_cutoff_*(IMS_bgCut_ + IMS_dgCut_ *IMS_X_o_cutoff_);
|
|
double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
|
|
IMS_egCut_ = - eterm * (tmp);
|
|
if (fabs(IMS_egCut_ - oldV) < 1.0E-14) {
|
|
converged = true;
|
|
}
|
|
}
|
|
if (!converged) {
|
|
throw CanteraError("HMWSoln::calcIMSCutoffParams_()",
|
|
" failed to converge on the g polynomial");
|
|
}
|
|
}
|
|
|
|
// Precalculate the MC Cutoff parameters
|
|
void HMWSoln::calcMCCutoffParams_()
|
|
{
|
|
MC_X_o_min_ = 0.35;
|
|
MC_X_o_cutoff_ = 0.6;
|
|
MC_slopepCut_ = 0.02;
|
|
MC_cpCut_ = 0.25;
|
|
|
|
// Initial starting values
|
|
MC_apCut_ = MC_X_o_min_;
|
|
MC_epCut_ = 0.0;
|
|
bool converged = false;
|
|
double oldV = 0.0;
|
|
int its;
|
|
double damp = 0.5;
|
|
for (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;
|
|
MC_bpCut_ = damp * MC_bpCutNew + (1-damp) * MC_bpCut_;
|
|
double MC_dpCutNew = ((- MC_apCut_/MC_cpCut_ + MC_bpCut_ - MC_bpCut_ * MC_X_o_cutoff_/MC_cpCut_)
|
|
/
|
|
(MC_X_o_cutoff_ * MC_X_o_cutoff_/MC_cpCut_ - 2.0 * MC_X_o_cutoff_));
|
|
MC_dpCut_ = damp * MC_dpCutNew + (1-damp) * MC_dpCut_;
|
|
double tmp = MC_apCut_ + MC_X_o_cutoff_*(MC_bpCut_ + MC_dpCut_ * MC_X_o_cutoff_);
|
|
double eterm = std::exp(- MC_X_o_cutoff_ / MC_cpCut_);
|
|
MC_epCut_ = - eterm * (tmp);
|
|
double diff = MC_epCut_ - oldV;
|
|
if (fabs(diff) < 1.0E-14) {
|
|
converged = true;
|
|
}
|
|
}
|
|
if (!converged) {
|
|
throw CanteraError("HMWSoln::calcMCCutoffParams_()",
|
|
" failed to converge on the p polynomial");
|
|
}
|
|
}
|
|
|
|
}
|