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

Clean up comments in Thermo class source files

Browse source
This commit is contained in:
Ray Speth 2015-10-20 19:55:49 -04:00
parent 14864e840b
commit 4382ae6c85
52 changed files with 1800 additions and 3438 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

22
include/cantera/numerics/Integrator.h
View file

@ -29,8 +29,8 @@ const int BAND = 32;
* Not all methods are supported by all integrators.
*/
enum MethodType {
BDF_Method, /**< Backward Differentiation */
Adams_Method /**< Adams */
BDF_Method, //!< Backward Differentiation
Adams_Method //! Adams
};
//! Specifies the method used for iteration.
@ -129,51 +129,51 @@ public:
return 0.0;
}
/** The current value of the solution of equation k. */
//! The current value of the solution of equation k.
virtual doublereal& solution(size_t k) {
warn("solution");
return m_dummy;
}
/** The current value of the solution of the system of equations. */
//! The current value of the solution of the system of equations.
virtual doublereal* solution() {
warn("solution");
return 0;
}
/** The number of equations. */
//! The number of equations.
virtual int nEquations() const {
warn("nEquations");
return 0;
}
/** The number of function evaluations. */
//! The number of function evaluations.
virtual int nEvals() const {
warn("nEvals");
return 0;
}
/** Set the maximum integration order that will be used. **/
//! Set the maximum integration order that will be used.
virtual void setMaxOrder(int n) {
warn("setMaxorder");
}
/** Set the solution method */
//! Set the solution method
virtual void setMethod(MethodType t) {
warn("setMethodType");
}
/** Set the linear iterator. */
//! Set the linear iterator.
virtual void setIterator(IterType t) {
warn("setInterator");
}
/** Set the maximum step size */
//! Set the maximum step size
virtual void setMaxStepSize(double hmax) {
warn("setMaxStepSize");
}
/** Set the minimum step size */
//! Set the minimum step size
virtual void setMinStepSize(double hmin) {
warn("setMinStepSize");
}

7
include/cantera/thermo/IdealMolalSoln.h
View file

@ -26,10 +26,6 @@
namespace Cantera
{
/** \addtogroup thermoprops */
/* @{
*/
/**
* This phase is based upon the mixing-rule assumption that
* all molality-based activity coefficients are equal
@ -93,6 +89,8 @@ namespace Cantera
* </idealMolalSolnCutoff>
* </activityCoefficients>
* </thermo>
*
* @ingroup thermoprops
*/
class IdealMolalSoln : public MolalityVPSSTP
{
@ -675,7 +673,6 @@ private:
void calcIMSCutoffParams_();
};
/* @} */
}
#endif

2
include/cantera/zeroD/ReactorBase.h
View file

@ -48,7 +48,7 @@ public:
m_name = name;
}
/** @name Methods to set up a simulation. */
//! @name Methods to set up a simulation.
//@{
/**

2
include/cantera/zeroD/ReactorNet.h
View file

@ -26,7 +26,7 @@ public:
ReactorNet();
virtual ~ReactorNet();
/** @name Methods to set up a simulation. */
//! @name Methods to set up a simulation.
//@{
/**

9
samples/cxx/flamespeed/flamespeed.cpp
View file

@ -161,12 +161,11 @@ int flamespeed(double phi)
int loglevel=1;
bool refine_grid = true;
/* Solve freely propagating flame*/
// Solve freely propagating flame
/* Linearly interpolate to find location where this
temperature would exist. The temperature at this
location will then be fixed for remainder of
calculation.*/
// Linearly interpolate to find location where this temperature would
// exist. The temperature at this location will then be fixed for
// remainder of calculation.
flow.fixTemperature();
refine_grid=false;

36
src/equil/BasisOptimize.cpp
View file

@ -144,8 +144,8 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
* current number of component species found.
*/
while (jr < nComponents) {
/* - Top of another loop point based on finding a linearly */
/* - independent species */
// Top of another loop point based on finding a linearly independent
// species
while (true) {
/*
* Search the remaining part of the mole number vector, molNum
@ -183,9 +183,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
#endif
molNum[kk] = USEDBEFORE;
/* *********************************************************** */
/* **** CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES ****** */
/* *********************************************************** */
// **** CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES ******
/*
* Modified Gram-Schmidt Method, p. 202 Dalquist
* QR factorization of a matrix without row pivoting.
@ -228,16 +226,12 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
tmp = sm[ml + jr*ne];
sa[jr] += tmp * tmp;
}
/* **************************************************** */
/* **** IF NORM OF NEW ROW .LT. 1E-3 REJECT ********** */
/* **************************************************** */
// **** IF NORM OF NEW ROW .LT. 1E-3 REJECT **********
if (sa[jr] > 1.0e-6) {
break;
}
}
/* ****************************************** */
/* **** REARRANGE THE DATA ****************** */
/* ****************************************** */
// **** REARRANGE THE DATA ******************
if (jr != k) {
if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl >= 1) {
kk = orderVectorSpecies[k];
@ -258,9 +252,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
return nComponents;
}
/* ****************************************************** */
/* **** EVALUATE THE STOICHIOMETRY ********************** */
/* ****************************************************** */
// **** EVALUATE THE STOICHIOMETRY **********************
/*
* Formulate the matrix problem for the stoichiometric
* coefficients. CX + B = 0
@ -361,7 +353,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
}
return nComponents;
} /* basopt() ************************************************************/
} // basopt()
static void print_stringTrunc(const char* str, int space, int alignment)
@ -526,10 +518,8 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
*/
eAbund[kk] = test;
/* *********************************************************** */
/* **** CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX */
/* **** LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX ****** */
/* *********************************************************** */
// **** CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX
// **** LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX ******
/*
* Modified Gram-Schmidt Method, p. 202 Dalquist
* QR factorization of a matrix without row pivoting.
@ -580,16 +570,12 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
double tmp = sm[ml + jr*nComponents];
sa[jr] += tmp * tmp;
}
/* **************************************************** */
/* **** IF NORM OF NEW ROW .LT. 1E-6 REJECT ********** */
/* **************************************************** */
// **** IF NORM OF NEW ROW .LT. 1E-6 REJECT **********
if (sa[jr] > 1.0e-6) {
break;
}
}
/* ****************************************** */
/* **** REARRANGE THE DATA ****************** */
/* ****************************************** */
// **** REARRANGE THE DATA ******************
if (jr != k) {
if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl > 0) {
kk = orderVectorElements[k];

440
src/thermo/DebyeHuckel.cpp
View file

@ -88,10 +88,8 @@ DebyeHuckel::DebyeHuckel(const DebyeHuckel& b) :
m_waterSS(0),
m_densWaterSS(1000.)
{
/*
* Use the assignment operator to do the brunt
* of the work for the copy constructor.
*/
// Use the assignment operator to do the brunt of the work for the copy
// constructor.
*this = b;
}
@ -164,9 +162,8 @@ int DebyeHuckel::eosType() const
return res;
}
//
// -------- Molar Thermodynamic Properties of the Solution ---------------
//
doublereal DebyeHuckel::enthalpy_mole() const
{
getPartialMolarEnthalpies(m_tmpV.data());
@ -196,9 +193,7 @@ doublereal DebyeHuckel::cv_mole() const
throw NotImplementedError("DebyeHuckel::cv_mole");
}
//
// ------- Mechanical Equation of State Properties ------------------------
//
doublereal DebyeHuckel::pressure() const
{
@ -213,32 +208,24 @@ void DebyeHuckel::setPressure(doublereal p)
void DebyeHuckel::setState_TP(doublereal t, doublereal p)
{
Phase::setTemperature(t);
/*
* Store the current pressure
*/
// Store the current pressure
m_Pcurrent = p;
/*
* update the standard state thermo
* -> This involves calling the water function and setting the pressure
*/
// update the standard state thermo. This involves calling the water
// function and setting the pressure
_updateStandardStateThermo();
/*
* Calculate all of the other standard volumes
* -> note these are constant for now
*/
// Calculate all of the other standard volumes. Note these are constant for
// now
calcDensity();
}
void DebyeHuckel::calcDensity()
{
if (m_waterSS) {
/*
* Store the internal density of the water SS.
* Note, we would have to do this for all other
* species if they had pressure dependent properties.
*/
// Store the internal density of the water SS. Note, we would have to do
// this for all other species if they had pressure dependent properties.
m_densWaterSS = m_waterSS->density();
}
double* vbar = &m_pp[0];
@ -276,9 +263,7 @@ void DebyeHuckel::setTemperature(const doublereal temp)
setState_TP(temp, m_Pcurrent);
}
//
// ------- Activities and Activity Concentrations
//
void DebyeHuckel::getActivityConcentrations(doublereal* c) const
{
@ -298,10 +283,9 @@ doublereal DebyeHuckel::standardConcentration(size_t k) const
void DebyeHuckel::getActivities(doublereal* ac) const
{
_updateStandardStateThermo();
/*
* Update the molality array, m_molalities()
* This requires an update due to mole fractions
*/
// Update the molality array, m_molalities(). This requires an update due to
// mole fractions
s_update_lnMolalityActCoeff();
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
@ -324,23 +308,18 @@ void DebyeHuckel::getMolalityActivityCoefficients(doublereal* acMolality) const
}
}
//
// ------ Partial Molar Properties of the Solution -----------------
//
void DebyeHuckel::getChemPotentials(doublereal* mu) const
{
double xx;
/*
* First get the standard chemical potentials in
* molar form.
* -> this requires updates of standard state as a function
* of T and P
*/
// First get the standard chemical potentials in molar form. This requires
// updates of standard state as a function of T and P
getStandardChemPotentials(mu);
/*
* Update the activity coefficients
* This also updates the internal molality array.
*/
// Update the activity coefficients. This also updates the internal molality
// array.
s_update_lnMolalityActCoeff();
double xmolSolvent = moleFraction(m_indexSolvent);
for (size_t k = 0; k < m_kk; k++) {
@ -356,27 +335,21 @@ void DebyeHuckel::getChemPotentials(doublereal* mu) const
void DebyeHuckel::getPartialMolarEnthalpies(doublereal* hbar) const
{
/*
* Get the nondimensional standard state enthalpies
*/
// Get the nondimensional standard state enthalpies
getEnthalpy_RT(hbar);
/*
* Dimensionalize it.
*/
// Dimensionalize it.
for (size_t k = 0; k < m_kk; k++) {
hbar[k] *= RT();
}
/*
* Check to see whether activity coefficients are temperature
* dependent. If they are, then calculate the their temperature
* derivatives and add them into the result.
*/
// Check to see whether activity coefficients are temperature
// dependent. If they are, then calculate the their temperature
// derivatives and add them into the result.
double dAdT = dA_DebyedT_TP();
if (dAdT != 0.0) {
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the
// internally stored molalities.
s_update_lnMolalityActCoeff();
s_update_dlnMolalityActCoeff_dT();
for (size_t k = 0; k < m_kk; k++) {
@ -387,26 +360,21 @@ void DebyeHuckel::getPartialMolarEnthalpies(doublereal* hbar) const
void DebyeHuckel::getPartialMolarEntropies(doublereal* sbar) const
{
/*
* Get the standard state entropies at the temperature
* and pressure of the solution.
*/
// Get the standard state entropies at the temperature and pressure of the
// solution.
getEntropy_R(sbar);
/*
* Dimensionalize the entropies
*/
// Dimensionalize the entropies
for (size_t k = 0; k < m_kk; k++) {
sbar[k] *= GasConstant;
}
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnMolalityActCoeff();
/*
* First we will add in the obvious dependence on the T
* term out front of the log activity term
*/
// First we will add in the obvious dependence on the T term out front of
// the log activity term
doublereal mm;
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
@ -417,11 +385,10 @@ void DebyeHuckel::getPartialMolarEntropies(doublereal* sbar) const
double xmolSolvent = moleFraction(m_indexSolvent);
mm = std::max(SmallNumber, xmolSolvent);
sbar[m_indexSolvent] -= GasConstant *(log(mm) + m_lnActCoeffMolal[m_indexSolvent]);
/*
* Check to see whether activity coefficients are temperature
* dependent. If they are, then calculate the their temperature
* derivatives and add them into the result.
*/
// Check to see whether activity coefficients are temperature dependent. If
// they are, then calculate the their temperature derivatives and add them
// into the result.
double dAdT = dA_DebyedT_TP();
if (dAdT != 0.0) {
s_update_dlnMolalityActCoeff_dT();
@ -434,9 +401,8 @@ void DebyeHuckel::getPartialMolarEntropies(doublereal* sbar) const
void DebyeHuckel::getPartialMolarVolumes(doublereal* vbar) const
{
getStandardVolumes(vbar);
/*
* Update the derivatives wrt the activity coefficients.
*/
// Update the derivatives wrt the activity coefficients.
s_update_lnMolalityActCoeff();
s_update_dlnMolalityActCoeff_dP();
for (size_t k = 0; k < m_kk; k++) {
@ -446,27 +412,18 @@ void DebyeHuckel::getPartialMolarVolumes(doublereal* vbar) const
void DebyeHuckel::getPartialMolarCp(doublereal* cpbar) const
{
/*
* Get the nondimensional Gibbs standard state of the
* species at the T and P of the solution.
*/
getCp_R(cpbar);
for (size_t k = 0; k < m_kk; k++) {
cpbar[k] *= GasConstant;
}
/*
* Check to see whether activity coefficients are temperature
* dependent. If they are, then calculate the their temperature
* derivatives and add them into the result.
*/
// Check to see whether activity coefficients are temperature dependent. If
// they are, then calculate the their temperature derivatives and add them
// into the result.
double dAdT = dA_DebyedT_TP();
if (dAdT != 0.0) {
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally
// stored molalities.
s_update_lnMolalityActCoeff();
s_update_dlnMolalityActCoeff_dT();
s_update_d2lnMolalityActCoeff_dT2();
@ -477,9 +434,7 @@ void DebyeHuckel::getPartialMolarCp(doublereal* cpbar) const
}
}
/*
* -------------- Utilities -------------------------------
*/
// -------------- Utilities -------------------------------
void DebyeHuckel::initThermo()
{
@ -491,7 +446,8 @@ void DebyeHuckel::initThermo()
initLengths();
}
//! Utility function to assign an integer value from a string for the ElectrolyteSpeciesType field.
//! Utility function to assign an integer value from a string for the
//! ElectrolyteSpeciesType field.
/*!
* @param estString input string that will be interpreted
*/
@ -530,19 +486,15 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Find the Thermo XML node
*/
// Find the Thermo XML node
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("DebyeHuckel::initThermoXML",
"no thermo XML node");
}
XML_Node& thermoNode = phaseNode.child("thermo");
/*
* Determine the form of the Debye-Huckel model,
* m_formDH. We will use this information to size arrays below.
*/
// Determine the form of the Debye-Huckel model, m_formDH. We will use this
// information to size arrays below.
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& scNode = thermoNode.child("activityCoefficients");
m_formDH = DHFORM_DILUTE_LIMIT;
@ -564,16 +516,12 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
} else {
/*
* If there is no XML node named "activityCoefficients", assume
* that we are doing the extreme dilute limit assumption
*/
// If there is no XML node named "activityCoefficients", assume
// that we are doing the extreme dilute limit assumption
m_formDH = DHFORM_DILUTE_LIMIT;
}
/*
* Possibly change the form of the standard concentrations
*/
// Possibly change the form of the standard concentrations
if (thermoNode.hasChild("standardConc")) {
XML_Node& scNode = thermoNode.child("standardConc");
m_formGC = 2;
@ -596,13 +544,10 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Reconcile the solvent name and index.
*/
/*
* Get the Name of the Solvent:
* <solvent> solventName </solvent>
*/
// Reconcile the solvent name and index.
// Get the Name of the Solvent:
// <solvent> solventName </solvent>
std::string solventName = "";
if (thermoNode.hasChild("solvent")) {
XML_Node& scNode = thermoNode.child("solvent");
@ -633,17 +578,13 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
" should be first species");
}
/*
* Initialize all of the lengths of arrays in the object
* now that we know what species are in the phase.
*/
// Initialize all of the lengths of arrays in the object now that we know
// what species are in the phase.
initThermo();
/*
* Now go get the specification of the standard states for
* species in the solution. This includes the molar volumes
* data blocks for incompressible 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"],
@ -673,19 +614,16 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
if (k == 0) {
if (modelString == "wateriapws" || modelString == "real_water" ||
modelString == "waterpdss") {
/*
* Initialize the water standard state model
*/
// Initialize the water standard state model
m_waterSS = dynamic_cast<PDSS_Water*>(providePDSS(0));
if (!m_waterSS) {
throw CanteraError("HMWSoln::installThermoXML",
"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.
*/
// 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();
@ -707,17 +645,14 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Go get all of the coefficients and factors in the
* activityCoefficients XML block
*/
// 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
*/
// Look for parameters for A_Debye
if (acNode.hasChild("A_Debye")) {
XML_Node* ss = acNode.findByName("A_Debye");
string modelStringa = ss->attrib("model");
@ -735,24 +670,18 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Initialize the water property calculator. It will share
* the internal eos water calculator.
*/
// Initialize the water property calculator. It will share the internal
// eos water calculator.
if (m_form_A_Debye == A_DEBYE_WATER) {
m_waterProps.reset(new WaterProps(m_waterSS));
}
/*
* Look for parameters for B_Debye
*/
// Look for parameters for B_Debye
if (acNode.hasChild("B_Debye")) {
m_B_Debye = getFloat(acNode, "B_Debye");
}
/*
* Look for parameters for B_dot
*/
// Look for parameters for B_dot
if (acNode.hasChild("B_dot")) {
if (m_formDH == DHFORM_BETAIJ ||
m_formDH == DHFORM_DILUTE_LIMIT ||
@ -761,9 +690,7 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
"B_dot entry in the wrong DH form");
}
double bdot_common = getFloat(acNode, "B_dot");
/*
* Set B_dot parameters for charged species
*/
// Set B_dot parameters for charged species
for (size_t k = 0; k < m_kk; k++) {
double z_k = charge(k);
if (fabs(z_k) > 0.0001) {
@ -774,25 +701,19 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Look for Parameters for the Maximum Ionic Strength
*/
// Look for Parameters for the Maximum Ionic Strength
if (acNode.hasChild("maxIonicStrength")) {
m_maxIionicStrength = getFloat(acNode, "maxIonicStrength");
}
/*
* Look for Helgeson Parameters
*/
// Look for Helgeson Parameters
if (acNode.hasChild("UseHelgesonFixedForm")) {
m_useHelgesonFixedForm = true;
} else {
m_useHelgesonFixedForm = false;
}
/*
* Look for parameters for the Ionic radius
*/
// Look for parameters for the Ionic radius
if (acNode.hasChild("ionicRadius")) {
XML_Node& irNode = acNode.child("ionicRadius");
@ -810,45 +731,35 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* If the Debye-Huckel form is BDOT_AK, we can
* have separate values for the denominator's ionic
* size. -> That's how the activity coefficient is
* parameterized. In this case only do we allow the
* code to read in these parameters.
*/
// If the Debye-Huckel form is BDOT_AK, we can have separate values
// for the denominator's ionic size. -> That's how the activity
// coefficient is parameterized. In this case only do we allow the
// code to read in these parameters.
if (m_formDH == DHFORM_BDOT_AK) {
/*
* Define a string-string map, and interpret the
* value of the XML element as binary pairs separated
* by colons, e.g.:
* Na+:3.0
* Cl-:4.0
* H+:9.0
* OH-:3.5
* Read them into the map.
*/
// Define a string-string map, and interpret the value of the
// XML element as binary pairs separated by colons, e.g.:
// Na+:3.0
// Cl-:4.0
// H+:9.0
// OH-:3.5
// Read them into the map.
map<string, string> m;
getMap(irNode, m);
/*
* Iterate over the map pairs, interpreting the
* first string as a species in the current phase.
* If no match is made, silently ignore the
* lack of agreement (HKM -> may be changed in the
* future).
*/
// Iterate over the map pairs, interpreting the first string as
// a species in the current phase. If no match is made, silently
// ignore the lack of agreement (HKM -> may be changed in the
// future).
for (const auto& b : m) {
size_t kk = speciesIndex(b.first);
m_Aionic[kk] = fpValue(b.second) * Afactor;
}
}
}
/*
* Get the matrix of coefficients for the Beta
* binary interaction parameters. We assume here that
* this matrix is symmetric, so that we only have to
* input 1/2 of the values.
*/
// Get the matrix of coefficients for the Beta binary interaction
// parameters. We assume here that this matrix is symmetric, so that we
// only have to input 1/2 of the values.
if (acNode.hasChild("DHBetaMatrix")) {
if (m_formDH == DHFORM_BETAIJ ||
m_formDH == DHFORM_PITZER_BETAIJ) {
@ -861,22 +772,16 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Fill in parameters for the calculation of the
* stoichiometric Ionic Strength
*
* The default is that stoich charge is the same as the
* regular charge.
*/
// Fill in parameters for the calculation of the stoichiometric Ionic
// Strength. The default is that stoich charge is the same as the
// regular charge.
m_speciesCharge_Stoich.resize(m_kk, 0.0);
for (size_t k = 0; k < m_kk; k++) {
m_speciesCharge_Stoich[k] = m_speciesCharge[k];
}
/*
* First look at the species database.
* -> Look for the subelement "stoichIsMods"
* in each of the species SS databases.
*/
// 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();
size_t jj = xspecies.size();
for (size_t k = 0; k < m_kk; k++) {
@ -899,9 +804,7 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Now look at the activity coefficient database
*/
// Now look at the activity coefficient database
if (acNodePtr && acNodePtr->hasChild("stoichIsMods")) {
XML_Node& sIsNode = acNodePtr->child("stoichIsMods");
map<std::string, std::string> msIs;
@ -914,13 +817,9 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* 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.
*/
// 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;
@ -934,11 +833,9 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
m_electrolyteSpeciesType[m_indexSolvent] = cEST_solvent;
/*
* First look at the species database.
* -> Look for the subelement "stoichIsMods"
* in each of the species SS databases.
*/
// 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();
for (size_t k = 0; k < m_kk; k++) {
std::string kname = speciesName(k);
@ -951,9 +848,8 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
}
/*
* Then look at the phase thermo specification
*/
// Then look at the phase thermo specification
if (acNodePtr && acNodePtr->hasChild("electrolyteSpeciesType")) {
XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType");
map<std::string, std::string> msEST;
@ -968,9 +864,7 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Lastly set the state
*/
// Lastly set the state
if (phaseNode.hasChild("state")) {
XML_Node& stateNode = phaseNode.child("state");
setStateFromXML(stateNode);
@ -1075,25 +969,17 @@ double DebyeHuckel::dA_DebyedP_TP(double tempArg, double presArg) const
return dAdP;
}
/*
* ---------- Other Property Functions
*/
// ---------- Other Property Functions
double DebyeHuckel::AionicRadius(int k) const
{
return m_Aionic[k];
}
/*
* ------------ Private and Restricted Functions ------------------
*/
// ------------ Private and Restricted Functions ------------------
void DebyeHuckel::initLengths()
{
/*
* Obtain the limits of the temperature from the species
* thermo handler's limits.
*/
m_electrolyteSpeciesType.resize(m_kk, cEST_polarNeutral);
m_speciesSize.resize(m_kk);
m_Aionic.resize(m_kk, 0.0);
@ -1141,9 +1027,7 @@ double DebyeHuckel::_osmoticCoeffHelgesonFixedForm() const
double DebyeHuckel::_lnactivityWaterHelgesonFixedForm() const
{
/*
* Update the internally stored vector of molalities
*/
// Update the internally stored vector of molalities
calcMolalities();
double oc = _osmoticCoeffHelgesonFixedForm();
double sum = 0.0;
@ -1161,17 +1045,15 @@ double DebyeHuckel::_lnactivityWaterHelgesonFixedForm() const
void DebyeHuckel::s_update_lnMolalityActCoeff() const
{
double z_k, zs_k1, zs_k2;
/*
* Update the internally stored vector of molalities
*/
// Update the internally stored vector of molalities
calcMolalities();
/*
* Calculate the apparent (real) ionic strength.
*
* Note this is not the stoichiometric ionic strengh,
* where reactions of ions forming neutral salts
* are ignorred in calculating the ionic strength.
*/
// Calculate the apparent (real) ionic strength.
//
// Note this is not the stoichiometric ionic strengh, where reactions of
// ions forming neutral salts are ignorred in calculating the ionic
// strength.
m_IionicMolality = 0.0;
for (size_t k = 0; k < m_kk; k++) {
z_k = m_speciesCharge[k];
@ -1180,9 +1062,7 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
m_IionicMolality /= 2.0;
m_IionicMolality = std::min(m_IionicMolality, m_maxIionicStrength);
/*
* Calculate the stoichiometric ionic charge
*/
// Calculate the stoichiometric ionic charge
m_IionicMolalityStoich = 0.0;
for (size_t k = 0; k < m_kk; k++) {
z_k = m_speciesCharge[k];
@ -1197,21 +1077,14 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
m_IionicMolalityStoich /= 2.0;
m_IionicMolalityStoich = std::min(m_IionicMolalityStoich, m_maxIionicStrength);
/*
* Possibly update the stored value of the
* Debye-Huckel parameter A_Debye
* This parameter appears on the top of the activity
* coefficient expression.
* It depends on T (and P), as it depends explicitly
* on the temperature. Also, the dielectric constant
* is usually a fairly strong function of T, also.
*/
// Possibly update the stored value of the Debye-Huckel parameter A_Debye
// This parameter appears on the top of the activity coefficient expression.
// It depends on T (and P), as it depends explicitly on the temperature.
// Also, the dielectric constant is usually a fairly strong function of T,
// also.
m_A_Debye = A_Debye_TP();
/*
* Calculate a safe value for the mole fraction
* of the solvent
*/
// Calculate a safe value for the mole fraction of the solvent
double xmolSolvent = moleFraction(m_indexSolvent);
xmolSolvent = std::max(8.689E-3, xmolSolvent);
@ -1276,10 +1149,8 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
lnActivitySolvent -=
m_Mnaught * log(10.0) * m_IionicMolality * tmp / 2.0;
/*
* Special section to implement the Helgeson fixed form
* for the water brine activity coefficient.
*/
// Special section to implement the Helgeson fixed form for the water
// brine activity coefficient.
if (m_useHelgesonFixedForm) {
lnActivitySolvent = _lnactivityWaterHelgesonFixedForm();
}
@ -1390,12 +1261,10 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const
default:
throw CanteraError("DebyeHuckel::s_update_lnMolalityActCoeff", "ERROR");
}
/*
* Above, we calculated the ln(activitySolvent). Translate that
* into the molar-based activity coefficient by dividing by
* the solvent mole fraction. Solvents are not on the molality
* scale.
*/
// Above, we calculated the ln(activitySolvent). Translate that into the
// molar-based activity coefficient by dividing by the solvent mole
// fraction. Solvents are not on the molality scale.
xmolSolvent = moleFraction(m_indexSolvent);
m_lnActCoeffMolal[m_indexSolvent] =
lnActivitySolvent - log(xmolSolvent);
@ -1412,10 +1281,8 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const
}
return;
}
/*
* Calculate a safe value for the mole fraction
* of the solvent
*/
// Calculate a safe value for the mole fraction of the solvent
double xmolSolvent = moleFraction(m_indexSolvent);
xmolSolvent = std::max(8.689E-3, xmolSolvent);
double sqrtI = sqrt(m_IionicMolality);
@ -1534,10 +1401,7 @@ void DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2() const
return;
}
/*
* Calculate a safe value for the mole fraction
* of the solvent
*/
// Calculate a safe value for the mole fraction of the solvent
double xmolSolvent = moleFraction(m_indexSolvent);
xmolSolvent = std::max(8.689E-3, xmolSolvent);
double sqrtI = sqrt(m_IionicMolality);
@ -1651,10 +1515,8 @@ void DebyeHuckel::s_update_dlnMolalityActCoeff_dP() const
}
return;
}
/*
* Calculate a safe value for the mole fraction
* of the solvent
*/
// Calculate a safe value for the mole fraction of the solvent
double xmolSolvent = moleFraction(m_indexSolvent);
xmolSolvent = std::max(8.689E-3, xmolSolvent);
double sqrtI = sqrt(m_IionicMolality);

37
src/thermo/FixedChemPotSSTP.cpp
View file

@ -21,9 +21,8 @@
namespace Cantera
{
/*
* ---- Constructors -------
*/
// ---- Constructors -------
FixedChemPotSSTP::FixedChemPotSSTP() :
chemPot_(0.0)
@ -98,18 +97,14 @@ ThermoPhase* FixedChemPotSSTP::duplMyselfAsThermoPhase() const
return new FixedChemPotSSTP(*this);
}
/*
* ---- Utilities -----
*/
// ---- Utilities -----
int FixedChemPotSSTP::eosType() const
{
return cFixedChemPot;
}
/*
* ----- Mechanical Equation of State ------
*/
// ----- Mechanical Equation of State ------
doublereal FixedChemPotSSTP::pressure() const
{
@ -131,9 +126,7 @@ doublereal FixedChemPotSSTP::thermalExpansionCoeff() const
return 0.0;
}
/*
* ---- Chemical Potentials and Activities ----
*/
// ---- Chemical Potentials and Activities ----
void FixedChemPotSSTP::getActivityConcentrations(doublereal* c) const
{
@ -150,18 +143,14 @@ doublereal FixedChemPotSSTP::logStandardConc(size_t k) const
return 0.0;
}
/*
* ---- Partial Molar Properties of the Solution ----
*/
// ---- Partial Molar Properties of the Solution ----
void FixedChemPotSSTP::getPartialMolarVolumes(doublereal* vbar) const
{
vbar[0] = 0.0;
}
/*
* Properties of the Standard State of the Species in the Solution
*/
// Properties of the Standard State of the Species in the Solution
void FixedChemPotSSTP::getStandardChemPotentials(doublereal* mu0) const
{
@ -198,9 +187,7 @@ void FixedChemPotSSTP::getStandardVolumes(doublereal* vbar) const
vbar[0] = 0.0;
}
/*
* ---- Thermodynamic Values for the Species Reference States ----
*/
// ---- Thermodynamic Values for the Species Reference States ----
void FixedChemPotSSTP::getIntEnergy_RT_ref(doublereal* urt) const
{
@ -232,15 +219,11 @@ void FixedChemPotSSTP::getCp_R_ref(doublereal* cpr) const
cpr[0] = 0.0;
}
/*
* ---- Initialization and Internal functions
*/
// ---- Initialization and Internal functions
void FixedChemPotSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
/*
* Find the Thermo XML node
*/
// Find the Thermo XML node
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("FixedChemPotSSTP::initThermoXML", "no thermo XML node");
}

43
src/thermo/GibbsExcessVPSSTP.cpp
View file

@ -4,11 +4,10 @@
* employ excess Gibbs free energy formulations
* (see \ref thermoprops and class \link Cantera::GibbsExcessVPSSTP GibbsExcessVPSSTP\endlink).
*
* Header file for a derived class of ThermoPhase that handles
* variable pressure standard state methods for calculating
* thermodynamic properties that are further based upon expressions
* for the excess Gibbs free energy expressed as a function of
* the mole fractions.
* Header file for a derived class of ThermoPhase that handles variable pressure
* standard state methods for calculating thermodynamic properties that are
* further based upon expressions for the excess Gibbs free energy expressed as
* a function of the mole fractions.
*/
/*
* Copyright (2009) Sandia Corporation. Under the terms of
@ -83,9 +82,8 @@ void GibbsExcessVPSSTP::setConcentrations(const doublereal* const c)
getMoleFractions(moleFractions_.data());
}
/*
* ------------ Mechanical Properties ------------------------------
*/
// ------------ Mechanical Properties ------------------------------
void GibbsExcessVPSSTP::setPressure(doublereal p)
{
setState_TP(temperature(), p);
@ -105,25 +103,19 @@ void GibbsExcessVPSSTP::calcDensity()
void GibbsExcessVPSSTP::setState_TP(doublereal t, doublereal p)
{
Phase::setTemperature(t);
/*
* Store the current pressure
*/
// Store the current pressure
m_Pcurrent = p;
/*
* update the standard state thermo
* -> This involves calling the water function and setting the pressure
*/
// update the standard state thermo. This involves calling the water
// function and setting the pressure
updateStandardStateThermo();
/*
* Calculate the partial molar volumes, and then the density of the fluid
*/
// Calculate the partial molar volumes, and then the density of the fluid
calcDensity();
}
/*
* - Activities, Standard States, Activity Concentrations -----------
*/
// - Activities, Standard States, Activity Concentrations -----------
void GibbsExcessVPSSTP::getActivityConcentrations(doublereal* c) const
{
getActivities(c);
@ -171,14 +163,11 @@ void GibbsExcessVPSSTP::getElectrochemPotentials(doublereal* mu) const
}
}
/*
* ------------ Partial Molar Properties of the Solution ------------
*/
// ------------ Partial Molar Properties of the Solution ------------
void GibbsExcessVPSSTP::getPartialMolarVolumes(doublereal* vbar) const
{
/*
* Get the standard state values in m^3 kmol-1
*/
// Get the standard state values in m^3 kmol-1
getStandardVolumes(vbar);
}

1526
src/thermo/HMWSoln.cpp
View file

File diff suppressed because it is too large Load diff

281
src/thermo/HMWSoln_input.cpp
View file

@ -67,10 +67,9 @@ void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
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
*/
// 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;
@ -94,13 +93,10 @@ void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
XML_Node& xmlChild = BinSalt.child(iChild);
string stemp = xmlChild.name();
string nodeName = lowercase(stemp);
/*
* Process the binary salt child elements
*/
// Process the binary salt child elements
if (nodeName == "beta0") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, vParams, false, "", "beta0");
size_t nParamsFound = vParams.size();
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
@ -133,9 +129,7 @@ void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
}
}
if (nodeName == "beta1") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, vParams, false, "", "beta1");
size_t nParamsFound = vParams.size();
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
@ -200,9 +194,7 @@ void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
}
}
if (nodeName == "cphi") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, vParams, false, "", "Cphi");
size_t nParamsFound = vParams.size();
if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
@ -264,10 +256,9 @@ void HMWSoln::readXMLThetaAnion(XML_Node& BinSalt)
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
*/
// 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;
@ -343,10 +334,9 @@ void HMWSoln::readXMLThetaCation(XML_Node& BinSalt)
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
*/
// 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;
@ -426,10 +416,9 @@ void HMWSoln::readXMLPsiCommonCation(XML_Node& BinSalt)
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
*/
// 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;
@ -561,10 +550,9 @@ void HMWSoln::readXMLPsiCommonAnion(XML_Node& BinSalt)
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
*/
// 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;
@ -692,10 +680,9 @@ void HMWSoln::readXMLLambdaNeutral(XML_Node& BinSalt)
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
*/
// 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;
@ -766,10 +753,8 @@ void HMWSoln::readXMLMunnnNeutral(XML_Node& BinSalt)
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
*/
// 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;
@ -841,10 +826,9 @@ void HMWSoln::readXMLZetaCation(const XML_Node& BinSalt)
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
*/
// 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;
@ -961,10 +945,9 @@ void HMWSoln::constructPhaseFile(std::string inputFile, std::string id_)
throw CanteraError("HMWSoln:constructPhaseFile","could not open "
+path+" for reading.");
}
/*
* The phase object automatically constructs an XML object.
* Use this object to store information.
*/
// The phase object automatically constructs an XML object.
// Use this object to store information.
XML_Node fxml;
fxml.build(fin);
XML_Node* fxml_phase = findXMLPhase(&fxml, id_);
@ -987,18 +970,14 @@ void HMWSoln::constructPhaseXML(XML_Node& phaseNode, std::string id_)
}
}
/*
* Find the Thermo XML node
*/
// 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
*/
// Possibly change the form of the standard concentrations
if (thermoNode.hasChild("standardConc")) {
XML_Node& scNode = thermoNode.child("standardConc");
m_formGC = 2;
@ -1021,10 +1000,9 @@ void HMWSoln::constructPhaseXML(XML_Node& phaseNode, std::string id_)
}
}
}
/*
* Get the Name of the Solvent:
* <solvent> solventName </solvent>
*/
// Get the Name of the Solvent:
// <solvent> solventName </solvent>
string solventName = "";
if (thermoNode.hasChild("solvent")) {
XML_Node& scNode = thermoNode.child("solvent");
@ -1037,10 +1015,8 @@ void HMWSoln::constructPhaseXML(XML_Node& phaseNode, std::string id_)
solventName = nameSolventa[0];
}
/*
* Determine the form of the Pitzer model,
* We will use this information to size arrays below.
*/
// 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");
string stemp = scNode.attrib("model");
@ -1056,10 +1032,9 @@ void HMWSoln::constructPhaseXML(XML_Node& phaseNode, std::string id_)
+ formString);
}
}
/*
* Determine the form of the temperature dependence
* of the Pitzer activity coefficient model.
*/
// Determine the form of the temperature dependence of the Pitzer
// activity coefficient model.
stemp = scNode.attrib("TempModel");
formString = lowercase(stemp);
if (formString != "") {
@ -1076,11 +1051,9 @@ void HMWSoln::constructPhaseXML(XML_Node& phaseNode, std::string id_)
}
}
/*
* Determine the reference temperature
* of the Pitzer activity coefficient model's temperature
* dependence formulation: defaults to 25C
*/
// 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 != "") {
@ -1090,11 +1063,9 @@ void HMWSoln::constructPhaseXML(XML_Node& phaseNode, std::string id_)
}
}
/*
* 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
*/
// Call the importPhase() function. This will import all of the species into
// the phase. This will also handle all of the solvent and solute standard
// states
importPhase(phaseNode, this);
}
@ -1108,18 +1079,14 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Find the Thermo XML node
*/
// 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
*/
// Possibly change the form of the standard concentrations
if (thermoNode.hasChild("standardConc")) {
XML_Node& scNode = thermoNode.child("standardConc");
m_formGC = 2;
@ -1143,10 +1110,8 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Determine the form of the Pitzer model,
* We will use this information to size arrays below.
*/
// 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");
string stemp = scNode.attrib("model");
@ -1163,10 +1128,8 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Determine the form of the temperature dependence
* of the Pitzer activity coefficient model.
*/
// Determine the form of the temperature dependence of the Pitzer
// activity coefficient model.
stemp = scNode.attrib("TempModel");
formString = lowercase(stemp);
if (formString != "") {
@ -1183,11 +1146,9 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Determine the reference temperature
* of the Pitzer activity coefficient model's temperature
* dependence formulation: defaults to 25C
*/
// 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 != "") {
@ -1197,10 +1158,8 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Get the Name of the Solvent:
* <solvent> solventName </solvent>
*/
// Get the Name of the Solvent:
// <solvent> solventName </solvent>
string solventName = "";
if (thermoNode.hasChild("solvent")) {
XML_Node& scNode = thermoNode.child("solvent");
@ -1213,15 +1172,11 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
solventName = nameSolventa[0];
}
/*
* Initialize all of the lengths of arrays in the object
* now that we know what species are in the phase.
*/
// 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.
*/
// Reconcile the solvent name and index.
for (size_t k = 0; k < m_kk; k++) {
string sname = speciesName(k);
if (solventName == sname) {
@ -1246,11 +1201,9 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
" 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.
*/
// 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"],
@ -1279,20 +1232,18 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
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.
*/
// 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.
*/
// 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();
@ -1317,34 +1268,25 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Initialize the water property calculator. It will share
* the internal eos water calculator.
*/
// Initialize the water property calculator. It will share the internal eos
// water calculator.
m_waterProps.reset(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.
*/
// 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] = charge(k);
}
/*
* Go get all of the coefficients and factors in the
* activityCoefficients XML block
*/
// 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
*/
// 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;
@ -1361,16 +1303,12 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Look for Parameters for the Maximum Ionic Strength
*/
// Look for Parameters for the Maximum Ionic Strength
if (acNode.hasChild("maxIonicStrength")) {
m_maxIionicStrength = getFloat(acNode, "maxIonicStrength");
}
/*
* Look for parameters for the Ionic radius
*/
// Look for parameters for the Ionic radius
if (acNode.hasChild("ionicRadius")) {
XML_Node& irNode = acNode.child("ionicRadius");
double Afactor = 1.0;
@ -1388,11 +1326,8 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* First look at the species database.
* -> Look for the subelement "stoichIsMods"
* in each of the species SS databases.
*/
// 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();
for (size_t k = 0; k < m_kk; k++) {
size_t jmap = npos;
@ -1411,9 +1346,7 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Now look at the activity coefficient database
*/
// Now look at the activity coefficient database
if (acNodePtr && acNodePtr->hasChild("stoichIsMods")) {
XML_Node& sIsNode = acNodePtr->child("stoichIsMods");
map<string, string> msIs;
@ -1427,20 +1360,16 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Loop through the children getting multiple instances of
* parameters
*/
// 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);
string 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.
*/
// 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") {
@ -1463,13 +1392,10 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
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.
*/
// 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(charge(k)) > 0.0001) {
m_electrolyteSpeciesType[k] = cEST_chargedSpecies;
@ -1483,11 +1409,9 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
m_electrolyteSpeciesType[m_indexSolvent] = cEST_solvent;
/*
* First look at the species database.
* -> Look for the subelement "stoichIsMods"
* in each of the species SS databases.
*/
// 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();
for (size_t k = 0; k < m_kk; k++) {
const XML_Node* spPtr = xspecies[k];
@ -1499,9 +1423,8 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
}
/*
* Then look at the phase thermo specification
*/
// Then look at the phase thermo specification
if (acNodePtr && acNodePtr->hasChild("electrolyteSpeciesType")) {
XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType");
map<string, string> msEST;
@ -1526,9 +1449,9 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
setMoleFSolventMin(1.0E-5);
MolalityVPSSTP::initThermoXML(phaseNode, id_);
/*
* Lastly calculate the charge balance and then add stuff until the charges compensate
*/
// Lastly calculate the charge balance and then add stuff until the charges
// compensate
vector_fp mf(m_kk, 0.0);
getMoleFractions(mf.data());
bool notDone = true;

19
src/thermo/IdealGasPhase.cpp
View file

@ -38,10 +38,8 @@ IdealGasPhase::IdealGasPhase(const IdealGasPhase& right) :
m_p0(right.m_p0),
m_logc0(right.m_logc0)
{
/*
* Use the assignment operator to do the brunt
* of the work for the copy constructor.
*/
// Use the assignment operator to do the brunt of the work for the copy
// constructor.
*this = right;
}
@ -280,14 +278,11 @@ void IdealGasPhase::setToEquilState(const doublereal* mu_RT)
{
const vector_fp& grt = gibbs_RT_ref();
/*
* Within the method, we protect against inf results if the
* exponent is too high.
*
* If it is too low, we set
* the partial pressure to zero. This capability is needed
* by the elemental potential method.
*/
// Within the method, we protect against inf results if the exponent is too
// high.
//
// If it is too low, we set the partial pressure to zero. This capability is
// needed by the elemental potential method.
doublereal pres = 0.0;
for (size_t k = 0; k < m_kk; k++) {
double tmp = -grt[k] + mu_RT[k];

148
src/thermo/IdealMolalSoln.cpp
View file

@ -4,13 +4,12 @@
* state (see \ref thermoprops
* and class \link Cantera::IdealMolalSoln IdealMolalSoln\endlink).
*
* Definition file for a derived class of ThermoPhase that handles
* variable pressure standard state methods for calculating
* thermodynamic properties that are further based upon
* activities on the molality scale. The Ideal molal
* solution assumes that all molality-based activity
* coefficients are equal to one. This turns out, actually, to be
* highly nonlinear when the solvent densities get low.
* Definition file for a derived class of ThermoPhase that handles variable
* pressure standard state methods for calculating thermodynamic properties that
* are further based upon activities on the molality scale. The Ideal molal
* solution assumes that all molality-based activity coefficients are equal to
* one. This turns out, actually, to be highly nonlinear when the solvent
* densities get low.
*/
/*
* Copyright (2006) Sandia Corporation. Under the terms of
@ -48,10 +47,8 @@ IdealMolalSoln::IdealMolalSoln() :
IdealMolalSoln::IdealMolalSoln(const IdealMolalSoln& b) :
MolalityVPSSTP(b)
{
/*
* Use the assignment operator to do the brunt
* of the work for the copy constructor.
*/
// Use the assignment operator to do the brunt of the work for the copy
// constructor.
*this = b;
}
@ -169,9 +166,7 @@ doublereal IdealMolalSoln::cv_mole() const
throw NotImplementedError("IdealMolalSoln::cv_mole");
}
//
// ------- Mechanical Equation of State Properties ------------------------
//
void IdealMolalSoln::setPressure(doublereal p)
{
@ -226,9 +221,7 @@ void IdealMolalSoln::setState_TP(doublereal temp, doublereal pres)
calcDensity();
}
//
// ------- Activities and Activity Concentrations
//
void IdealMolalSoln::getActivityConcentrations(doublereal* c) const
{
@ -266,10 +259,9 @@ doublereal IdealMolalSoln::standardConcentration(size_t k) const
void IdealMolalSoln::getActivities(doublereal* ac) const
{
_updateStandardStateThermo();
/*
* Update the molality array, m_molalities()
* This requires an update due to mole fractions
*/
// Update the molality array, m_molalities(). This requires an update due to
// mole fractions
if (IMS_typeCutoff_ == 0) {
calcMolalities();
for (size_t k = 0; k < m_kk; k++) {
@ -284,9 +276,8 @@ void IdealMolalSoln::getActivities(doublereal* ac) const
} else {
s_updateIMS_lnMolalityActCoeff();
/*
* Now calculate the array of activities.
*/
// Now calculate the array of activities.
for (size_t k = 1; k < m_kk; k++) {
ac[k] = m_molalities[k] * exp(IMS_lnActCoeffMolal_[k]);
}
@ -317,30 +308,23 @@ void IdealMolalSoln::getMolalityActivityCoefficients(doublereal* acMolality) con
}
}
//
// ------ Partial Molar Properties of the Solution -----------------
//
void IdealMolalSoln::getChemPotentials(doublereal* mu) const
{
// Assertion is made for speed
AssertThrow(m_indexSolvent == 0, "solvent not the first species");
/*
* First get the standard chemical potentials
* -> this requires updates of standard state as a function
* of T and P
* These are defined at unit molality.
*/
// First get the standard chemical potentials. This requires updates of
// standard state as a function of T and P These are defined at unit
// molality.
getStandardChemPotentials(mu);
/*
* Update the molality array, m_molalities()
* This requires an update due to mole fractions
*/
// Update the molality array, m_molalities(). This requires an update due to
// mole fractions
calcMolalities();
/*
* get the solvent mole fraction
*/
// get the solvent mole fraction
double xmolSolvent = moleFraction(m_indexSolvent);
if (IMS_typeCutoff_ == 0 || xmolSolvent > 3.* IMS_X_o_cutoff_/2.0) {
@ -348,18 +332,15 @@ void IdealMolalSoln::getChemPotentials(doublereal* mu) const
double xx = std::max(m_molalities[k], SmallNumber);
mu[k] += RT() * log(xx);
}
/*
* Do the solvent
* -> see my notes
*/
// Do the solvent
// -> see my notes
double xx = std::max(xmolSolvent, SmallNumber);
mu[m_indexSolvent] +=
(RT() * (xmolSolvent - 1.0) / xx);
} else {
/*
* Update the activity coefficients
* This also updates the internal molality array.
*/
// Update the activity coefficients. This also updates the internal
// molality array.
s_updateIMS_lnMolalityActCoeff();
for (size_t k = 1; k < m_kk; k++) {
@ -394,15 +375,12 @@ void IdealMolalSoln::getPartialMolarEntropies(doublereal* sbar) const
double xmolSolvent = moleFraction(m_indexSolvent);
sbar[m_indexSolvent] -= (GasConstant * (xmolSolvent - 1.0) / xmolSolvent);
} else {
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally
// stored molalities.
s_updateIMS_lnMolalityActCoeff();
/*
* First we will add in the obvious dependence on the T
* term out front of the log activity term
*/
// First we will add in the obvious dependence on the T term out front
// of the log activity term
doublereal mm;
for (size_t k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
@ -423,19 +401,15 @@ void IdealMolalSoln::getPartialMolarVolumes(doublereal* vbar) const
void IdealMolalSoln::getPartialMolarCp(doublereal* cpbar) const
{
/*
* Get the nondimensional Gibbs standard state of the
* species at the T and P of the solution.
*/
// Get the nondimensional Gibbs standard state of the species at the T and P
// of the solution.
getCp_R(cpbar);
for (size_t k = 0; k < m_kk; k++) {
cpbar[k] *= GasConstant;
}
}
/*
* -------------- Utilities -------------------------------
*/
// -------------- Utilities -------------------------------
void IdealMolalSoln::initThermo()
{
@ -445,17 +419,13 @@ void IdealMolalSoln::initThermo()
void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
/*
* Find the Thermo XML node
*/
// Find the Thermo XML node
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("IdealMolalSoln::initThermoXML",
"no thermo XML node");
}
/*
* Initialize the whole thermo object, using a virtual function.
*/
// Initialize the whole thermo object, using a virtual function.
initThermo();
if (id_.size() > 0 && phaseNode.id() != id_) {
@ -463,18 +433,14 @@ void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
"phasenode and Id are incompatible");
}
/*
* Find the Thermo XML node
*/
// Find the Thermo XML node
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("IdealMolalSoln::initThermo",
"no thermo XML node");
}
XML_Node& thermoNode = phaseNode.child("thermo");
/*
* Possible change the form of the standard concentrations
*/
// Possible change the form of the standard concentrations
if (thermoNode.hasChild("standardConc")) {
XML_Node& scNode = thermoNode.child("standardConc");
m_formGC = 2;
@ -493,10 +459,8 @@ void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Get the Name of the Solvent:
* <solvent> solventName </solvent>
*/
// Get the Name of the Solvent:
// <solvent> solventName </solvent>
std::string solventName = "";
if (thermoNode.hasChild("solvent")) {
std::vector<std::string> nameSolventa;
@ -551,9 +515,7 @@ void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Reconcile the solvent name and index.
*/
// Reconcile the solvent name and index.
for (size_t k = 0; k < m_kk; k++) {
if (solventName == speciesName(k)) {
m_indexSolvent = k;
@ -572,9 +534,7 @@ void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
" should be first species");
}
/*
* Now go get the molar volumes
*/
// Now go get the molar volumes
XML_Node& speciesList = phaseNode.child("speciesArray");
XML_Node* speciesDB =
get_XML_NameID("speciesData", speciesList["datasrc"],
@ -594,26 +554,20 @@ void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
MolalityVPSSTP::initThermoXML(phaseNode, id_);
setMoleFSolventMin(1.0E-5);
/*
* Set the state
*/
// Set the state
if (phaseNode.hasChild("state")) {
XML_Node& stateNode = phaseNode.child("state");
setStateFromXML(stateNode);
}
}
/*
* ------------ Private and Restricted Functions ------------------
*/
// ------------ Private and Restricted Functions ------------------
void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
{
/*
* Calculate the molalities. Currently, the molalities
* may not be current with respect to the contents of the
* State objects' data.
*/
// Calculate the molalities. Currently, the molalities may not be current
// with respect to the contents of the State objects' data.
calcMolalities();
double xmolSolvent = moleFraction(m_indexSolvent);
@ -640,9 +594,7 @@ void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
IMS_lnActCoeffMolal_[m_indexSolvent] = log(IMS_gamma_o_min_);
return;
} else {
/*
* If we are in the middle region, calculate the connecting polynomials
*/
// If we are in the middle region, calculate the connecting polynomials
double xminus = xmolSolvent - IMS_X_o_cutoff_/2.0;
double xminus2 = xminus * xminus;
double xminus3 = xminus2 * xminus;
@ -715,10 +667,6 @@ void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
void IdealMolalSoln::initLengths()
{
/*
* Obtain the limits of the temperature from the species
* thermo handler's limits.
*/
m_pp.resize(m_kk);
m_speciesMolarVolume.resize(m_kk);
m_tmpV.resize(m_kk);

99
src/thermo/IdealSolidSolnPhase.cpp
View file

@ -108,9 +108,7 @@ int IdealSolidSolnPhase::eosType() const
return res;
}
/********************************************************************
* Molar Thermodynamic Properties of the Solution
********************************************************************/
// Molar Thermodynamic Properties of the Solution
doublereal IdealSolidSolnPhase::enthalpy_mole() const
{
@ -133,33 +131,25 @@ doublereal IdealSolidSolnPhase::cp_mole() const
return GasConstant * mean_X(cp_R_ref());
}
/********************************************************************
* Mechanical Equation of State
********************************************************************/
// Mechanical Equation of State
void IdealSolidSolnPhase::calcDensity()
{
/*
* Calculate the molarVolume of the solution (m**3 kmol-1)
*/
// Calculate the molarVolume of the solution (m**3 kmol-1)
const doublereal* const dtmp = moleFractdivMMW();
double invDens = dot(m_speciesMolarVolume.begin(),
m_speciesMolarVolume.end(), dtmp);
/*
* Set the density in the parent State object directly,
* by calling the Phase::setDensity() function.
*/
// Set the density in the parent State object directly, by calling the
// Phase::setDensity() function.
Phase::setDensity(1.0/invDens);
}
void IdealSolidSolnPhase::setDensity(const doublereal rho)
{
/*
* Unless the input density is exactly equal to the density
* calculated and stored in the State object, we throw an
* exception. This is because the density is NOT an
* independent variable.
*/
// Unless the input density is exactly equal to the density calculated and
// stored in the State object, we throw an exception. This is because the
// density is NOT an independent variable.
if (rho != density()) {
throw CanteraError("IdealSolidSolnPhase::setDensity",
"Density is not an independent variable");
@ -208,9 +198,7 @@ void IdealSolidSolnPhase::setConcentrations(const doublereal* const c)
calcDensity();
}
/********************************************************************
* Chemical Potentials and Activities
********************************************************************/
// Chemical Potentials and Activities
void IdealSolidSolnPhase::getActivityConcentrations(doublereal* c) const
{
@ -311,9 +299,7 @@ void IdealSolidSolnPhase::getChemPotentials_RT(doublereal* mu) const
}
}
/********************************************************************
* Partial Molar Properties
********************************************************************/
// Partial Molar Properties
void IdealSolidSolnPhase::getPartialMolarEnthalpies(doublereal* hbar) const
{
@ -343,9 +329,7 @@ void IdealSolidSolnPhase::getPartialMolarVolumes(doublereal* vbar) const
getStandardVolumes(vbar);
}
/*****************************************************************
* Properties of the Standard State of the Species in the Solution
*****************************************************************/
// Properties of the Standard State of the Species in the Solution
void IdealSolidSolnPhase::getPureGibbs(doublereal* gpure) const
{
@ -400,9 +384,7 @@ void IdealSolidSolnPhase::getStandardVolumes(doublereal* vol) const
copy(m_speciesMolarVolume.begin(), m_speciesMolarVolume.end(), vol);
}
/*********************************************************************
* Thermodynamic Values for the Species Reference States
*********************************************************************/
// Thermodynamic Values for the Species Reference States
void IdealSolidSolnPhase::getEnthalpy_RT_ref(doublereal* hrt) const
{
@ -466,9 +448,7 @@ const vector_fp& IdealSolidSolnPhase::entropy_R_ref() const
return m_s0_R;
}
/*********************************************************************
* Utility Functions
*********************************************************************/
// Utility Functions
void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
@ -477,10 +457,8 @@ void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string&
"phasenode and Id are incompatible");
}
/*
* Check on the thermo field. Must have:
* <thermo model="IdealSolidSolution" />
*/
// Check on the thermo field. Must have:
// <thermo model="IdealSolidSolution" />
if (phaseNode.hasChild("thermo")) {
XML_Node& thNode = phaseNode.child("thermo");
string mString = thNode.attrib("model");
@ -493,13 +471,11 @@ void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string&
"Unspecified thermo model");
}
/*
* Form of the standard concentrations. Must have one of:
*
* <standardConc model="unity" />
* <standardConc model="molar_volume" />
* <standardConc model="solvent_volume" />
*/
// Form of the standard concentrations. Must have one of:
//
// <standardConc model="unity" />
// <standardConc model="molar_volume" />
// <standardConc model="solvent_volume" />
if (phaseNode.hasChild("standardConc")) {
XML_Node& scNode = phaseNode.child("standardConc");
string formStringa = scNode.attrib("model");
@ -519,14 +495,11 @@ void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string&
"Unspecified standardConc model");
}
/*
* Initialize all of the lengths now that we know how many species
* there are in the phase.
*/
// Initialize all of the lengths now that we know how many species
// there are in the phase.
initLengths();
/*
* Now go get the molar volumes
*/
// Now go get the molar volumes
XML_Node& speciesList = phaseNode.child("speciesArray");
XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"],
&phaseNode.root());
@ -537,21 +510,15 @@ void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string&
m_speciesMolarVolume[k] = getFloat(*ss, "molarVolume", "toSI");
}
/*
* Call the base initThermo, which handles setting the initial
* state.
*/
// Call the base initThermo, which handles setting the initial state.
ThermoPhase::initThermoXML(phaseNode, id_);
}
void IdealSolidSolnPhase::initLengths()
{
/*
* Obtain the reference pressure by calling the ThermoPhase
* function refPressure, which in turn calls the
* species thermo reference pressure function of the
* same name.
*/
// Obtain the reference pressure by calling the ThermoPhase function
// refPressure, which in turn calls the species thermo reference pressure
// function of the same name.
m_Pref = refPressure();
m_h0_RT.resize(m_kk);
@ -568,8 +535,7 @@ void IdealSolidSolnPhase::setToEquilState(const doublereal* lambda_RT)
{
const vector_fp& grt = gibbs_RT_ref();
// set the pressure and composition to be consistent with
// the temperature,
// set the pressure and composition to be consistent with the temperature
doublereal pres = 0.0;
for (size_t k = 0; k < m_kk; k++) {
m_pp[k] = -grt[k];
@ -596,9 +562,8 @@ void IdealSolidSolnPhase::_updateThermo() const
{
doublereal tnow = temperature();
if (m_tlast != tnow) {
/*
* Update the thermodynamic functions of the reference state.
*/
// Update the thermodynamic functions of the reference state.
m_spthermo->update(tnow, m_cp0_R.data(), m_h0_RT.data(), m_s0_R.data());
m_tlast = tnow;
doublereal rrt = 1.0 / (GasConstant * tnow);

54
src/thermo/IdealSolnGasVPSS.cpp
View file

@ -55,14 +55,11 @@ IdealSolnGasVPSS::IdealSolnGasVPSS(const IdealSolnGasVPSS& b) :
IdealSolnGasVPSS& IdealSolnGasVPSS::operator=(const IdealSolnGasVPSS& b)
{
if (&b != this) {
/*
* Mostly, this is a passthrough to the underlying
* assignment operator for the ThermoPhae parent object.
*/
// Mostly, this is a passthrough to the underlying assignment operator
// for the ThermoPhae parent object.
VPStandardStateTP::operator=(b);
/*
* However, we have to handle data that we own.
*/
// However, we have to handle data that we own.
m_idealGas = b.m_idealGas;
m_formGC = b.m_formGC;
}
@ -82,9 +79,7 @@ int IdealSolnGasVPSS::eosType() const
return cIdealSolnGasVPSS_iscv;
}
/*
* ------------Molar Thermodynamic Properties -------------------------
*/
// ------------Molar Thermodynamic Properties -------------------------
doublereal IdealSolnGasVPSS::enthalpy_mole() const
{
@ -118,9 +113,7 @@ void IdealSolnGasVPSS::setPressure(doublereal p)
void IdealSolnGasVPSS::calcDensity()
{
/*
* Calculate the molarVolume of the solution (m**3 kmol-1)
*/
// Calculate the molarVolume of the solution (m**3 kmol-1)
if (m_idealGas) {
double dens = (m_Pcurrent * meanMolecularWeight()
/(GasConstant * temperature()));
@ -129,10 +122,8 @@ void IdealSolnGasVPSS::calcDensity()
const doublereal* const dtmp = moleFractdivMMW();
const vector_fp& vss = m_VPSS_ptr->getStandardVolumes();
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.
*/
// Set the density in the parent State object directly
Phase::setDensity(dens);
}
}
@ -200,9 +191,7 @@ void IdealSolnGasVPSS::getActivityCoefficients(doublereal* ac) const
}
}
/*
* ---- Partial Molar Properties of the Solution -----------------
*/
// ---- Partial Molar Properties of the Solution -----------------
void IdealSolnGasVPSS::getChemPotentials_RT(doublereal* muRT) const
{
@ -265,14 +254,11 @@ void IdealSolnGasVPSS::setToEquilState(const doublereal* mu_RT)
updateStandardStateThermo();
const vector_fp& grt = m_VPSS_ptr->Gibbs_RT_ref();
/*
* Within the method, we protect against inf results if the
* exponent is too high.
*
* If it is too low, we set
* the partial pressure to zero. This capability is needed
* by the elemental potential method.
*/
// Within the method, we protect against inf results if the exponent is too
// high.
//
// If it is too low, we set the partial pressure to zero. This capability is
// needed by the elemental potential method.
doublereal pres = 0.0;
double m_p0 = m_VPSS_ptr->refPressure();
for (size_t k = 0; k < m_kk; k++) {
@ -314,13 +300,11 @@ void IdealSolnGasVPSS::initThermoXML(XML_Node& phaseNode, const std::string& id_
}
}
/*
* Form of the standard concentrations. Must have one of:
*
* <standardConc model="unity" />
* <standardConc model="molar_volume" />
* <standardConc model="solvent_volume" />
*/
// Form of the standard concentrations. Must have one of:
//
// <standardConc model="unity" />
// <standardConc model="molar_volume" />
// <standardConc model="solvent_volume" />
if (phaseNode.hasChild("standardConc")) {
if (m_idealGas) {
throw CanteraError("IdealSolnGasVPSS::initThermoXML",

216
src/thermo/IonsFromNeutralVPSSTP.cpp
View file

@ -5,11 +5,10 @@
* (see \ref thermoprops
* and class \link Cantera::IonsFromNeutralVPSSTP IonsFromNeutralVPSSTP\endlink).
*
* Header file for a derived class of ThermoPhase that handles
* variable pressure standard state methods for calculating
* thermodynamic properties that are further based upon expressions
* for the excess Gibbs free energy expressed as a function of
* the mole fractions.
* Header file for a derived class of ThermoPhase that handles variable pressure
* standard state methods for calculating thermodynamic properties that are
* further based upon expressions for the excess Gibbs free energy expressed as
* a function of the mole fractions.
*/
/*
* Copyright (2009) Sandia Corporation. Under the terms of
@ -89,13 +88,11 @@ IonsFromNeutralVPSSTP::operator=(const IonsFromNeutralVPSSTP& b)
return *this;
}
/*
* If we own the underlying neutral molecule phase, then we do a deep
* copy. If not, we do a shallow copy. We get a valid pointer for
* neutralMoleculePhase_ first, because we need it to assign the pointers
* within the PDSS_IonsFromNeutral object. which is done in the
* GibbsExcessVPSSTP::operator=(b) step.
*/
// If we own the underlying neutral molecule phase, then we do a deep copy.
// If not, we do a shallow copy. We get a valid pointer for
// neutralMoleculePhase_ first, because we need it to assign the pointers
// within the PDSS_IonsFromNeutral object. which is done in the
// GibbsExcessVPSSTP::operator=(b) step.
if (IOwnNThermoPhase_) {
if (b.neutralMoleculePhase_) {
delete neutralMoleculePhase_;
@ -160,10 +157,9 @@ void IonsFromNeutralVPSSTP::constructPhaseFile(std::string inputFile, std::strin
throw CanteraError("MargulesVPSSTP:constructPhaseFile","could not open "
+path+" for reading.");
}
/*
* The phase object automatically constructs an XML object.
* Use this object to store information.
*/
// The phase object automatically constructs an XML object.
// Use this object to store information.
XML_Node fxml;
fxml.build(fin);
XML_Node* fxml_phase = findXMLPhase(&fxml, id_);
@ -183,27 +179,21 @@ void IonsFromNeutralVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string i
"phasenode and Id are incompatible");
}
/*
* Find the thermo XML node
*/
// Find the thermo XML node
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
"no thermo XML node");
}
XML_Node& thermoNode = phaseNode.child("thermo");
/*
* Make sure that the thermo model is IonsFromNeutralMolecule
*/
// Make sure that the thermo model is IonsFromNeutralMolecule
string formString = lowercase(thermoNode.attrib("model"));
if (formString != "ionsfromneutralmolecule") {
throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
"model name isn't IonsFromNeutralMolecule: " + formString);
}
/*
* Find the Neutral Molecule Phase
*/
// Find the Neutral Molecule Phase
if (!thermoNode.hasChild("neutralMoleculePhase")) {
throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
"no neutralMoleculePhase XML node");
@ -216,33 +206,25 @@ void IonsFromNeutralVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string i
"neut_ptr = 0");
}
/*
* Create the neutralMolecule ThermoPhase if we haven't already
*/
// Create the neutralMolecule ThermoPhase if we haven't already
if (!neutralMoleculePhase_) {
neutralMoleculePhase_ = newPhase(*neut_ptr);
}
/*
* 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
*/
// 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
importPhase(phaseNode, this);
}
/*
* -------------- Utilities -------------------------------
*/
// -------------- Utilities -------------------------------
int IonsFromNeutralVPSSTP::eosType() const
{
return cIonsFromNeutral;
}
/*
* ------------ Molar Thermodynamic Properties ----------------------
*/
// ------------ Molar Thermodynamic Properties ----------------------
doublereal IonsFromNeutralVPSSTP::enthalpy_mole() const
{
@ -275,9 +257,7 @@ doublereal IonsFromNeutralVPSSTP::cv_mole() const
return mean_X(m_pp);
}
/*
* - Activities, Standard States, Activity Concentrations -----------
*/
// -- Activities, Standard States, Activity Concentrations -----------
void IonsFromNeutralVPSSTP::getDissociationCoeffs(vector_fp& coeffs,
vector_fp& charges, std::vector<size_t>& neutMolIndex) const
@ -289,31 +269,23 @@ void IonsFromNeutralVPSSTP::getDissociationCoeffs(vector_fp& coeffs,
void IonsFromNeutralVPSSTP::getActivityCoefficients(doublereal* ac) const
{
/*
* Update the activity coefficients
*/
// Update the activity coefficients
s_update_lnActCoeff();
/*
* take the exp of the internally stored coefficients.
*/
// take the exp of the internally stored coefficients.
for (size_t k = 0; k < m_kk; k++) {
ac[k] = exp(lnActCoeff_Scaled_[k]);
}
}
/*
* --------- Partial Molar Properties of the Solution -------------
*/
// --------- Partial Molar Properties of the Solution -------------
void IonsFromNeutralVPSSTP::getChemPotentials(doublereal* mu) const
{
size_t icat, jNeut;
doublereal xx, fact2;
/*
* Get the standard chemical potentials of netural molecules
*/
// Get the standard chemical potentials of netural molecules
neutralMoleculePhase_->getStandardChemPotentials(muNeutralMolecule_.data());
doublereal RT_ = GasConstant * temperature();
@ -328,7 +300,7 @@ void IonsFromNeutralVPSSTP::getChemPotentials(doublereal* mu) const
// Do the cation list
for (size_t k = 0; k < cationList_.size(); k++) {
//! Get the id for the next cation
// Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
xx = std::max(SmallNumber, moleFractions_[icat]);
@ -364,20 +336,16 @@ void IonsFromNeutralVPSSTP::getChemPotentials(doublereal* mu) const
void IonsFromNeutralVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const
{
/*
* Get the nondimensional standard state enthalpies
*/
// Get the nondimensional standard state enthalpies
getEnthalpy_RT(hbar);
/*
* dimensionalize it.
*/
// dimensionalize it.
for (size_t k = 0; k < m_kk; k++) {
hbar[k] *= RT();
}
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnActCoeff();
s_update_dlnActCoeffdT();
for (size_t k = 0; k < m_kk; k++) {
@ -387,14 +355,11 @@ void IonsFromNeutralVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const
void IonsFromNeutralVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
{
/*
* Get the nondimensional standard state entropies
*/
// Get the nondimensional standard state entropies
getEntropy_R(sbar);
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnActCoeff();
s_update_dlnActCoeffdT();
@ -402,9 +367,8 @@ void IonsFromNeutralVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
double xx = std::max(moleFractions_[k], SmallNumber);
sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - temperature() * dlnActCoeffdT_Scaled_[k];
}
/*
* dimensionalize it.
*/
// dimensionalize it.
for (size_t k = 0; k < m_kk; k++) {
sbar[k] *= GasConstant;
}
@ -454,25 +418,19 @@ void IonsFromNeutralVPSSTP::setPressure(doublereal p)
void IonsFromNeutralVPSSTP::setState_TP(doublereal t, doublereal p)
{
/*
* This is a two phase process. First, we calculate the standard states
* within the neutral molecule phase.
*/
// This is a two phase process. First, we calculate the standard states
// within the neutral molecule phase.
neutralMoleculePhase_->setState_TP(t, p);
VPStandardStateTP::setState_TP(t,p);
/*
* Calculate the partial molar volumes, and then the density of the fluid
*/
// Calculate the partial molar volumes, and then the density of the fluid
Phase::setDensity(neutralMoleculePhase_->density());
}
void IonsFromNeutralVPSSTP::calcIonMoleFractions(doublereal* const mf) const
{
/*
* Download the neutral mole fraction vector into the
* vector, NeutralMolecMoleFractions_[]
*/
// Download the neutral mole fraction vector into the vector,
// NeutralMolecMoleFractions_[]
neutralMoleculePhase_->getMoleFractions(NeutralMolecMoleFractions_.data());
// Zero the mole fractions
@ -480,9 +438,7 @@ void IonsFromNeutralVPSSTP::calcIonMoleFractions(doublereal* const mf) const
mf[k] = 0.0;
}
/*
* Use the formula matrix to calculate the relative mole numbers.
*/
// Use the formula matrix to calculate the relative mole numbers.
for (size_t jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
for (size_t k = 0; k < m_kk; k++) {
double fmij = fm_neutralMolec_ions_[k + jNeut * m_kk];
@ -490,9 +446,7 @@ void IonsFromNeutralVPSSTP::calcIonMoleFractions(doublereal* const mf) const
}
}
/*
* Normalize the new mole fractions
*/
// Normalize the new mole fractions
doublereal sum = 0.0;
for (size_t k = 0; k < m_kk; k++) {
sum += mf[k];
@ -508,7 +462,7 @@ void IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions() const
doublereal fmij;
doublereal sum = 0.0;
//! Zero the vector we are trying to find.
// Zero the vector we are trying to find.
for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
NeutralMolecMoleFractions_[k] = 0.0;
}
@ -536,7 +490,7 @@ void IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions() const
}
for (size_t k = 0; k < cationList_.size(); k++) {
//! Get the id for the next cation
// Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
if (jNeut != npos) {
@ -565,12 +519,12 @@ void IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions() const
}
for (size_t k = 0; k < m_kk; k++) {
if (fabs(moleFractionsTmp_[k]) > 1.0E-13) {
//! Check to see if we have in fact found the inverse.
// Check to see if we have in fact found the inverse.
if (anionList_[0] != k) {
throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions",
"neutral molecule calc error");
} else {
//! For the single anion case, we will allow some slippage
// For the single anion case, we will allow some slippage
if (fabs(moleFractionsTmp_[k]) > 1.0E-5) {
throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions",
"neutral molecule calc error - anion");
@ -606,8 +560,8 @@ void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal* const
{
doublereal sumy, sumdy;
//check sum dx = 0
//! Zero the vector we are trying to find.
// check sum dx = 0
// Zero the vector we are trying to find.
for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
y_[k] = 0.0;
dy[k] = 0.0;
@ -623,7 +577,7 @@ void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal* const
case cIonSolnType_SINGLEANION:
for (size_t k = 0; k < cationList_.size(); k++) {
//! Get the id for the next cation
// Get the id for the next cation
size_t icat = cationList_[k];
size_t jNeut = fm_invert_ionForNeutral[icat];
if (jNeut != npos) {
@ -707,9 +661,7 @@ void IonsFromNeutralVPSSTP::setConcentrations(const doublereal* const c)
neutralMoleculePhase_->setMoleFractions(NeutralMolecMoleFractions_.data());
}
/*
* ------------ Partial Molar Properties of the Solution ------------
*/
// ------------ Partial Molar Properties of the Solution ------------
void IonsFromNeutralVPSSTP::initThermo()
{
@ -779,27 +731,21 @@ void IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string
"phasenode and Id are incompatible");
}
/*
* Find the Thermo XML node
*/
// Find the Thermo XML node
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
"no thermo XML node");
}
XML_Node& thermoNode = phaseNode.child("thermo");
/*
* Make sure that the thermo model is IonsFromNeutralMolecule
*/
// Make sure that the thermo model is IonsFromNeutralMolecule
string formString = lowercase(thermoNode.attrib("model"));
if (formString != "ionsfromneutralmolecule") {
throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
"model name isn't IonsFromNeutralMolecule: " + formString);
}
/*
* Find the Neutral Molecule Phase
*/
// Find the Neutral Molecule Phase
if (!thermoNode.hasChild("neutralMoleculePhase")) {
throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
"no neutralMoleculePhase XML node");
@ -812,9 +758,7 @@ void IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string
"neut_ptr = 0");
}
/*
* Create the neutralMolecule ThermoPhase if we haven't already
*/
// Create the neutralMolecule ThermoPhase if we haven't already
if (!neutralMoleculePhase_) {
neutralMoleculePhase_ = newPhase(*neut_ptr);
}
@ -934,23 +878,14 @@ void IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string
}
}
}
/*
* This includes the setStateFromXML calls
*/
// This includes the setStateFromXML calls
GibbsExcessVPSSTP::initThermoXML(phaseNode, id_);
/*
* There is one extra step here. We assure ourselves that we
* have charge conservation.
*/
}
void IonsFromNeutralVPSSTP::s_update_lnActCoeff() const
{
size_t icat, jNeut;
/*
* Get the activity coefficiens of the neutral molecules
*/
// Get the activity coefficiens of the neutral molecules
neutralMoleculePhase_->getLnActivityCoefficients(lnActCoeff_NeutralMolecule_.data());
switch (ionSolnType_) {
@ -959,7 +894,7 @@ void IonsFromNeutralVPSSTP::s_update_lnActCoeff() const
case cIonSolnType_SINGLEANION:
// Do the cation list
for (size_t k = 0; k < cationList_.size(); k++) {
//! Get the id for the next cation
// Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
double fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
@ -995,9 +930,7 @@ void IonsFromNeutralVPSSTP::getdlnActCoeffds(const doublereal dTds, const double
doublereal* dlnActCoeffds) const
{
size_t icat, jNeut;
/*
* Get the activity coefficients of the neutral molecules
*/
// Get the activity coefficients of the neutral molecules
if (!geThermo) {
for (size_t k = 0; k < m_kk; k++) {
dlnActCoeffds[k] = dXds[k] / moleFractions_[k];
@ -1016,7 +949,7 @@ void IonsFromNeutralVPSSTP::getdlnActCoeffds(const doublereal dTds, const double
case cIonSolnType_SINGLEANION:
// Do the cation list
for (size_t k = 0; k < cationList_.size(); k++) {
//! Get the id for the next cation
// Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
double fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
@ -1051,9 +984,8 @@ void IonsFromNeutralVPSSTP::getdlnActCoeffds(const doublereal dTds, const double
void IonsFromNeutralVPSSTP::s_update_dlnActCoeffdT() const
{
size_t icat, jNeut;
/*
* Get the activity coefficients of the neutral molecules
*/
// Get the activity coefficients of the neutral molecules
if (!geThermo) {
dlnActCoeffdT_Scaled_.assign(m_kk, 0.0);
return;
@ -1102,9 +1034,8 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeffdT() const
void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnX_diag() const
{
size_t icat, jNeut;
/*
* Get the activity coefficients of the neutral molecules
*/
// Get the activity coefficients of the neutral molecules
if (!geThermo) {
dlnActCoeffdlnX_diag_.assign(m_kk, 0.0);
return;
@ -1118,7 +1049,7 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnX_diag() const
case cIonSolnType_SINGLEANION:
// Do the cation list
for (size_t k = 0; k < cationList_.size(); k++) {
//! Get the id for the next cation
// Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
double fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
@ -1153,9 +1084,8 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnX_diag() const
void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN_diag() const
{
size_t icat, jNeut;
/*
* Get the activity coefficients of the neutral molecules
*/
// Get the activity coefficients of the neutral molecules
if (!geThermo) {
dlnActCoeffdlnN_diag_.assign(m_kk, 0.0);
return;
@ -1169,7 +1099,7 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN_diag() const
case cIonSolnType_SINGLEANION:
// Do the cation list
for (size_t k = 0; k < cationList_.size(); k++) {
//! Get the id for the next cation
// Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
double fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
@ -1206,9 +1136,7 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN() const
size_t kcat = 0, kNeut = 0, mcat = 0, mNeut = 0;
doublereal fmij = 0.0;
dlnActCoeffdlnN_.zero();
/*
* Get the activity coefficients of the neutral molecules
*/
// Get the activity coefficients of the neutral molecules
if (!geThermo) {
throw CanteraError("IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN()", "dynamic cast failed");
}

16
src/thermo/LatticePhase.cpp
View file

@ -296,10 +296,8 @@ void LatticePhase::initThermoXML(XML_Node& phaseNode, const std::string& id_)
"ids don't match");
}
/*
* Check on the thermo field. Must have:
* <thermo model="Lattice" />
*/
// Check on the thermo field. Must have:
// <thermo model="Lattice" />
if (phaseNode.hasChild("thermo")) {
XML_Node& thNode = phaseNode.child("thermo");
std::string mString = thNode.attrib("model");
@ -311,9 +309,8 @@ void LatticePhase::initThermoXML(XML_Node& phaseNode, const std::string& id_)
throw CanteraError("LatticePhase::initThermoXML",
"Unspecified thermo model");
}
/*
* Now go get the molar volumes. use the default if not found
*/
// Now go get the molar volumes. use the default if not found
XML_Node& speciesList = phaseNode.child("speciesArray");
XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"], &phaseNode.root());
@ -329,10 +326,7 @@ void LatticePhase::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Call the base initThermo, which handles setting the initial
* state.
*/
// Call the base initThermo, which handles setting the initial state.
ThermoPhase::initThermoXML(phaseNode, id_);
}

11
src/thermo/LatticeSolidPhase.cpp
View file

@ -227,10 +227,9 @@ void LatticeSolidPhase::getMoleFractions(doublereal* const x) const
for (size_t k = 0; k < nsp; k++) {
(x + strt)[k] /= sum;
}
/*
* At this point we can check against the mole fraction vector of the underlying LatticePhase objects and
* get the same answer.
*/
// At this point we can check against the mole fraction vector of the
// underlying LatticePhase objects and get the same answer.
if (DEBUG_MODE_ENABLED) {
m_lattice[n]->getMoleFractions(&m_x[strt]);
for (size_t k = 0; k < nsp; k++) {
@ -347,9 +346,7 @@ void LatticeSolidPhase::installSlavePhases(XML_Node* phaseNode)
addSpecies(lp->species(k));
kk++;
}
/*
* Add in the lattice stoichiometry constraint
*/
// Add in the lattice stoichiometry constraint
if (n > 0) {
string econ = "LC_" + int2str(n) + "_" + id();
size_t m = addElement(econ, 0.0, 0, 0.0, CT_ELEM_TYPE_LATTICERATIO);

159
src/thermo/MargulesVPSSTP.cpp
View file

@ -81,28 +81,20 @@ ThermoPhase* MargulesVPSSTP::duplMyselfAsThermoPhase() const
return new MargulesVPSSTP(*this);
}
/*
* - Activities, Standard States, Activity Concentrations -----------
*/
// -- Activities, Standard States, Activity Concentrations -----------
void MargulesVPSSTP::getLnActivityCoefficients(doublereal* lnac) const
{
/*
* Update the activity coefficients
*/
// Update the activity coefficients
s_update_lnActCoeff();
/*
* take the exp of the internally stored coefficients.
*/
// take the exp of the internally stored coefficients.
for (size_t k = 0; k < m_kk; k++) {
lnac[k] = lnActCoeff_Scaled_[k];
}
}
/*
* ------------ Partial Molar Properties of the Solution ------------
*/
// ------------ Partial Molar Properties of the Solution ------------
void MargulesVPSSTP::getElectrochemPotentials(doublereal* mu) const
{
@ -115,16 +107,11 @@ void MargulesVPSSTP::getElectrochemPotentials(doublereal* mu) const
void MargulesVPSSTP::getChemPotentials(doublereal* mu) const
{
/*
* First get the standard chemical potentials in
* molar form.
* -> this requires updates of standard state as a function
* of T and P
*/
// First get the standard chemical potentials in molar form. This requires
// updates of standard state as a function of T and P
getStandardChemPotentials(mu);
/*
* Update the activity coefficients
*/
// Update the activity coefficients
s_update_lnActCoeff();
for (size_t k = 0; k < m_kk; k++) {
double xx = std::max(moleFractions_[k], SmallNumber);
@ -175,20 +162,16 @@ doublereal MargulesVPSSTP::cv_mole() const
void MargulesVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const
{
/*
* Get the nondimensional standard state enthalpies
*/
// Get the nondimensional standard state enthalpies
getEnthalpy_RT(hbar);
/*
* dimensionalize it.
*/
// dimensionalize it.
for (size_t k = 0; k < m_kk; k++) {
hbar[k] *= RT();
}
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnActCoeff();
s_update_dlnActCoeff_dT();
for (size_t k = 0; k < m_kk; k++) {
@ -198,24 +181,19 @@ void MargulesVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const
void MargulesVPSSTP::getPartialMolarCp(doublereal* cpbar) const
{
/*
* Get the nondimensional standard state entropies
*/
// Get the nondimensional standard state entropies
getCp_R(cpbar);
double T = temperature();
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnActCoeff();
s_update_dlnActCoeff_dT();
for (size_t k = 0; k < m_kk; k++) {
cpbar[k] -= 2 * T * dlnActCoeffdT_Scaled_[k] + T * T * d2lnActCoeffdT2_Scaled_[k];
}
/*
* dimensionalize it.
*/
// dimensionalize it.
for (size_t k = 0; k < m_kk; k++) {
cpbar[k] *= GasConstant;
}
@ -223,15 +201,12 @@ void MargulesVPSSTP::getPartialMolarCp(doublereal* cpbar) const
void MargulesVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
{
/*
* Get the nondimensional standard state entropies
*/
// Get the nondimensional standard state entropies
getEntropy_R(sbar);
double T = temperature();
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnActCoeff();
s_update_dlnActCoeff_dT();
@ -239,9 +214,8 @@ void MargulesVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
double xx = std::max(moleFractions_[k], SmallNumber);
sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
}
/*
* dimensionalize it.
*/
// dimensionalize it.
for (size_t k = 0; k < m_kk; k++) {
sbar[k] *= GasConstant;
}
@ -251,9 +225,7 @@ void MargulesVPSSTP::getPartialMolarVolumes(doublereal* vbar) const
{
double T = temperature();
/*
* Get the standard state values in m^3 kmol-1
*/
// Get the standard state values in m^3 kmol-1
getStandardVolumes(vbar);
for (size_t i = 0; i < numBinaryInteractions_; i++) {
@ -294,28 +266,22 @@ void MargulesVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
}
/*
* Find the Thermo XML node
*/
// Find the Thermo XML node
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("MargulesVPSSTP::initThermoXML",
"no thermo XML node");
}
XML_Node& thermoNode = phaseNode.child("thermo");
/*
* Make sure that the thermo model is Margules
*/
// Make sure that the thermo model is Margules
string formString = lowercase(thermoNode.attrib("model"));
if (formString != "margules") {
throw CanteraError("MargulesVPSSTP::initThermoXML",
"model name isn't Margules: " + formString);
}
/*
* Go get all of the coefficients and factors in the
* activityCoefficients XML block
*/
// Go get all of the coefficients and factors in the activityCoefficients
// XML block
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& acNode = thermoNode.child("activityCoefficients");
string mStringa = acNode.attrib("model");
@ -325,20 +291,17 @@ void MargulesVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
}
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.
*/
// 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 (lowercase(xmlACChild.name()) == "binaryneutralspeciesparameters") {
readXMLBinarySpecies(xmlACChild);
}
}
}
/*
* Go down the chain
*/
// Go down the chain
GibbsExcessVPSSTP::initThermoXML(phaseNode, id_);
}
@ -475,9 +438,7 @@ void MargulesVPSSTP::s_update_dlnActCoeff_dlnN() const
double T = temperature();
dlnActCoeffdlnN_.zero();
/*
* Loop over the activity coefficient gamma_k
*/
// Loop over the activity coefficient gamma_k
for (size_t iK = 0; iK < m_kk; iK++) {
for (size_t iM = 0; iM < m_kk; iM++) {
double XM = moleFractions_[iM];
@ -594,19 +555,18 @@ void MargulesVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
if (bName == "") {
throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", "no speciesB attrib");
}
/*
* Find the index of the species in the current phase. It's not
* an error to not find the species. What this means is that the A-B interaction referred to in this
* block will be ignored.
*/
// Find the index of the species in the current phase. It's not an error to
// not find the species. What this means is that the A-B interaction
// referred to in this block will be ignored.
size_t aSpecies = speciesIndex(aName);
if (aSpecies == npos) {
return;
}
string aspName = speciesName(aSpecies);
// @TODO Figure out what the original reason is for putting an error condition for charged species
// Seems OK to me.
// @TODO Figure out what the original reason is for putting an error
// condition for charged species. Seems OK to me.
if (charge(aSpecies) != 0.0) {
throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies",
"speciesA has a charge: {}", charge(aSpecies));
@ -629,20 +589,17 @@ void MargulesVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
for (size_t iChild = 0; iChild < xmLBinarySpecies.nChildren(); iChild++) {
XML_Node& xmlChild = xmLBinarySpecies.child(iChild);
string nodeName = lowercase(xmlChild.name());
/*
* Process the binary species interaction parameters.
* They are in subblocks labeled:
* excessEnthalpy
* excessEntropy
* excessVolume_Enthalpy
* excessVolume_Entropy
* Other blocks are currently ignored.
* @TODO determine a policy about ignoring blocks that should or shouldn't be there.
*/
// Process the binary species interaction parameters.
// They are in subblocks labeled:
// excessEnthalpy
// excessEntropy
// excessVolume_Enthalpy
// excessVolume_Entropy
// Other blocks are currently ignored.
// @TODO determine a policy about ignoring blocks that should or shouldn't be there.
if (nodeName == "excessenthalpy") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, vParams, true, "toSI", "excessEnthalpy");
if (vParams.size() != 2) {
throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessEnthalpy for " + aspName
@ -654,9 +611,7 @@ void MargulesVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
}
if (nodeName == "excessentropy") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, vParams, true, "toSI", "excessEntropy");
if (vParams.size() != 2) {
throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessEntropy for " + aspName
@ -668,9 +623,7 @@ void MargulesVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
}
if (nodeName == "excessvolume_enthalpy") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Enthalpy");
if (vParams.size() != 2) {
throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessVolume_Enthalpy for " + aspName
@ -682,9 +635,7 @@ void MargulesVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
}
if (nodeName == "excessvolume_entropy") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Entropy");
if (vParams.size() != 2) {
throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessVolume_Entropy for " + aspName

54
src/thermo/MaskellSolidSolnPhase.cpp
View file

@ -66,9 +66,8 @@ void MaskellSolidSolnPhase::getActivityConcentrations(doublereal* c) const
}
}
/********************************************************************
* Molar Thermodynamic Properties of the Solution
********************************************************************/
// Molar Thermodynamic Properties of the Solution
doublereal MaskellSolidSolnPhase::enthalpy_mole() const
{
_updateThermo();
@ -93,18 +92,13 @@ doublereal MaskellSolidSolnPhase::entropy_mole() const
return s0 + GasConstant * (xlogx(1-rfm) - xlogx(rfm) - xlogx(1-r-rfm) - xlogx((1-fmval)*r) - xlogx(1-r) - xlogx(r));
}
/********************************************************************
* Mechanical Equation of State
********************************************************************/
// Mechanical Equation of State
void MaskellSolidSolnPhase::setDensity(const doublereal rho)
{
/*
* Unless the input density is exactly equal to the density
* calculated and stored in the State object, we throw an
* exception. This is because the density is NOT an
* independent variable.
*/
// Unless the input density is exactly equal to the density calculated and
// stored in the State object, we throw an exception. This is because the
// density is NOT an independent variable.
double dens = density();
if (rho != dens) {
throw CanteraError("MaskellSolidSolnPhase::setDensity",
@ -136,9 +130,7 @@ void MaskellSolidSolnPhase::setMolarDensity(const doublereal n)
"Density is not an independent variable");
}
/********************************************************************
* Chemical Potentials and Activities
********************************************************************/
// Chemical Potentials and Activities
void MaskellSolidSolnPhase::getActivityCoefficients(doublereal* ac) const
{
@ -182,9 +174,7 @@ void MaskellSolidSolnPhase::getChemPotentials_RT(doublereal* mu) const
}
}
/********************************************************************
* Partial Molar Properties
********************************************************************/
// Partial Molar Properties
void MaskellSolidSolnPhase::getPartialMolarEnthalpies(doublereal* hbar) const
{
@ -217,13 +207,13 @@ void MaskellSolidSolnPhase::getPureGibbs(doublereal* gpure) const
void MaskellSolidSolnPhase::getStandardChemPotentials(doublereal* mu) const
{
// What is the difference between this and getPureGibbs? IdealSolidSolnPhase gives the same for both
// What is the difference between this and getPureGibbs? IdealSolidSolnPhase
// gives the same for both
getPureGibbs(mu);
}
/*********************************************************************
* Utility Functions
*********************************************************************/
// Utility Functions
void MaskellSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
if (id_.size() > 0 && phaseNode.id() != id_) {
@ -231,10 +221,8 @@ void MaskellSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string
"phasenode and Id are incompatible");
}
/*
* Check on the thermo field. Must have:
* <thermo model="MaskellSolidSolution" />
*/
// Check on the thermo field. Must have:
// <thermo model="MaskellSolidSolution" />
if (phaseNode.hasChild("thermo")) {
XML_Node& thNode = phaseNode.child("thermo");
std::string mString = thNode.attrib("model");
@ -243,9 +231,7 @@ void MaskellSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string
"Unknown thermo model: " + mString);
}
/*
* Parse the enthalpy of mixing constant
*/
// Parse the enthalpy of mixing constant
if (thNode.hasChild("h_mix")) {
set_h_mix(fpValue(thNode.child("h_mix").value()));
} else {
@ -277,10 +263,7 @@ void MaskellSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string
"MaskellSolidSolution model requires exactly 2 species.");
}
/*
* Call the base initThermo, which handles setting the initial
* state.
*/
// Call the base initThermo, which handles setting the initial state.
VPStandardStateTP::initThermoXML(phaseNode, id_);
}
@ -289,9 +272,8 @@ void MaskellSolidSolnPhase::_updateThermo() const
assert(m_kk == 2);
static const int cacheId = m_cache.getId();
CachedScalar cached = m_cache.getScalar(cacheId);
/*
* Update the thermodynamic functions of the reference state.
*/
// Update the thermodynamic functions of the reference state.
doublereal tnow = temperature();
if (!cached.validate(tnow)) {
m_spthermo->update(tnow, m_cp0_R.data(), m_h0_RT.data(), m_s0_R.data());

28
src/thermo/MetalSHEelectrons.cpp
View file

@ -20,9 +20,7 @@
namespace Cantera
{
/*
* ---- Constructors -------
*/
// ---- Constructors -------
MetalSHEelectrons::MetalSHEelectrons()
{
@ -57,18 +55,14 @@ ThermoPhase* MetalSHEelectrons::duplMyselfAsThermoPhase() const
return new MetalSHEelectrons(*this);
}
/*
* ---- Utilities -----
*/
// ---- Utilities -----
int MetalSHEelectrons::eosType() const
{
return cMetalSHEelectrons;
}
/*
* ----- Mechanical Equation of State ------
*/
// ----- Mechanical Equation of State ------
doublereal MetalSHEelectrons::pressure() const
{
@ -90,9 +84,7 @@ doublereal MetalSHEelectrons::thermalExpansionCoeff() const
return 1.0/temperature();
}
/*
* ---- Chemical Potentials and Activities ----
*/
// ---- Chemical Potentials and Activities ----
void MetalSHEelectrons::getActivityConcentrations(doublereal* c) const
{
@ -109,9 +101,7 @@ doublereal MetalSHEelectrons::logStandardConc(size_t k) const
return 0.0;
}
/*
* Properties of the Standard State of the Species in the Solution
*/
// Properties of the Standard State of the Species in the Solution
void MetalSHEelectrons::getStandardChemPotentials(doublereal* mu0) const
{
@ -155,15 +145,11 @@ void MetalSHEelectrons::getIntEnergy_RT_ref(doublereal* urt) const
urt[0] = m_h0_RT[0] - m_p0 / molarDensity() / RT();
}
/*
* ---- Initialization and Internal functions
*/
// ---- Initialization and Internal functions
void MetalSHEelectrons::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
/*
* Find the Thermo XML node
*/
// Find the Thermo XML node
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("MetalSHEelectrons::initThermoXML",
"no thermo XML node");

36
src/thermo/MineralEQ3.cpp
View file

@ -23,9 +23,7 @@ using namespace std;
namespace Cantera
{
/*
* ---- Constructors -------
*/
// ---- Constructors -------
MineralEQ3::MineralEQ3(const std::string& infile, const std::string& id_)
{
@ -66,18 +64,14 @@ ThermoPhase* MineralEQ3::duplMyselfAsThermoPhase() const
return new MineralEQ3(*this);
}
/*
* ---- Utilities -----
*/
// ---- Utilities -----
int MineralEQ3::eosType() const
{
return cStoichSubstance;
}
/*
* ----- Mechanical Equation of State ------
*/
// ----- Mechanical Equation of State ------
doublereal MineralEQ3::pressure() const
{
@ -99,9 +93,7 @@ doublereal MineralEQ3::thermalExpansionCoeff() const
return 0.0;
}
/*
* ---- Chemical Potentials and Activities ----
*/
// ---- Chemical Potentials and Activities ----
void MineralEQ3::getActivityConcentrations(doublereal* c) const
{
@ -118,9 +110,7 @@ doublereal MineralEQ3::logStandardConc(size_t k) const
return 0.0;
}
/*
* Properties of the Standard State of the Species in the Solution
*/
// Properties of the Standard State of the Species in the Solution
void MineralEQ3::getStandardChemPotentials(doublereal* mu0) const
{
@ -158,9 +148,7 @@ void MineralEQ3::getIntEnergy_RT(doublereal* urt) const
urt[0] = m_h0_RT[0] - m_p0 / molarDensity() / RT();
}
/*
* ---- Thermodynamic Values for the Species Reference States ----
*/
// ---- Thermodynamic Values for the Species Reference States ----
void MineralEQ3::getIntEnergy_RT_ref(doublereal* urt) const
{
@ -168,9 +156,7 @@ void MineralEQ3::getIntEnergy_RT_ref(doublereal* urt) const
urt[0] = m_h0_RT[0] - m_p0 / molarDensity() / RT();
}
/*
* ---- Initialization and Internal functions
*/
// ---- Initialization and Internal functions
void MineralEQ3::setParameters(int n, doublereal* const c)
{
@ -185,9 +171,7 @@ void MineralEQ3::getParameters(int& n, doublereal* const c) const
void MineralEQ3::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
/*
* Find the Thermo XML node
*/
// Find the Thermo XML node
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("HMWSoln::initThermoXML",
"no thermo XML node");
@ -262,9 +246,7 @@ doublereal MineralEQ3::LookupGe(const std::string& elemName)
void MineralEQ3::convertDGFormation()
{
/*
* Ok let's get the element compositions and conversion factors.
*/
// Ok let's get the element compositions and conversion factors.
doublereal totalSum = 0.0;
for (size_t m = 0; m < nElements(); m++) {
double na = nAtoms(0, m);

133
src/thermo/MixedSolventElectrolyte.cpp
View file

@ -85,28 +85,20 @@ ThermoPhase* MixedSolventElectrolyte::duplMyselfAsThermoPhase() const
return new MixedSolventElectrolyte(*this);
}
/*
* - Activities, Standard States, Activity Concentrations -----------
*/
// - Activities, Standard States, Activity Concentrations -----------
void MixedSolventElectrolyte::getActivityCoefficients(doublereal* ac) const
{
/*
* Update the activity coefficients
*/
// Update the activity coefficients
s_update_lnActCoeff();
/*
* take the exp of the internally stored coefficients.
*/
// take the exp of the internally stored coefficients.
for (size_t k = 0; k < m_kk; k++) {
ac[k] = exp(lnActCoeff_Scaled_[k]);
}
}
/*
* ------------ Partial Molar Properties of the Solution ------------
*/
// ------------ Partial Molar Properties of the Solution ------------
void MixedSolventElectrolyte::getElectrochemPotentials(doublereal* mu) const
{
@ -119,16 +111,10 @@ void MixedSolventElectrolyte::getElectrochemPotentials(doublereal* mu) const
void MixedSolventElectrolyte::getChemPotentials(doublereal* mu) const
{
/*
* First get the standard chemical potentials in
* molar form.
* -> this requires updates of standard state as a function
* of T and P
*/
// First get the standard chemical potentials in molar form. This requires
// updates of standard state as a function of T and P
getStandardChemPotentials(mu);
/*
* Update the activity coefficients
*/
// Update the activity coefficients
s_update_lnActCoeff();
for (size_t k = 0; k < m_kk; k++) {
double xx = std::max(moleFractions_[k], SmallNumber);
@ -176,20 +162,15 @@ doublereal MixedSolventElectrolyte::cv_mole() const
void MixedSolventElectrolyte::getPartialMolarEnthalpies(doublereal* hbar) const
{
/*
* Get the nondimensional standard state enthalpies
*/
// Get the nondimensional standard state enthalpies
getEnthalpy_RT(hbar);
/*
* dimensionalize it.
*/
// dimensionalize it.
for (size_t k = 0; k < m_kk; k++) {
hbar[k] *= RT();
}
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnActCoeff();
s_update_dlnActCoeff_dT();
for (size_t k = 0; k < m_kk; k++) {
@ -199,24 +180,18 @@ void MixedSolventElectrolyte::getPartialMolarEnthalpies(doublereal* hbar) const
void MixedSolventElectrolyte::getPartialMolarCp(doublereal* cpbar) const
{
/*
* Get the nondimensional standard state entropies
*/
getCp_R(cpbar);
double T = temperature();
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnActCoeff();
s_update_dlnActCoeff_dT();
for (size_t k = 0; k < m_kk; k++) {
cpbar[k] -= 2 * T * dlnActCoeffdT_Scaled_[k] + T * T * d2lnActCoeffdT2_Scaled_[k];
}
/*
* dimensionalize it.
*/
// dimensionalize it.
for (size_t k = 0; k < m_kk; k++) {
cpbar[k] *= GasConstant;
}
@ -224,15 +199,12 @@ void MixedSolventElectrolyte::getPartialMolarCp(doublereal* cpbar) const
void MixedSolventElectrolyte::getPartialMolarEntropies(doublereal* sbar) const
{
/*
* Get the nondimensional standard state entropies
*/
// Get the nondimensional standard state entropies
getEntropy_R(sbar);
double T = temperature();
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the
// internally stored molalities.
s_update_lnActCoeff();
s_update_dlnActCoeff_dT();
@ -240,9 +212,7 @@ void MixedSolventElectrolyte::getPartialMolarEntropies(doublereal* sbar) const
double xx = std::max(moleFractions_[k], SmallNumber);
sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
}
/*
* dimensionalize it.
*/
// dimensionalize it.
for (size_t k = 0; k < m_kk; k++) {
sbar[k] *= GasConstant;
}
@ -252,9 +222,7 @@ void MixedSolventElectrolyte::getPartialMolarVolumes(doublereal* vbar) const
{
double T = temperature();
/*
* Get the standard state values in m^3 kmol-1
*/
// Get the standard state values in m^3 kmol-1
getStandardVolumes(vbar);
for (size_t iK = 0; iK < m_kk; iK++) {
@ -299,10 +267,8 @@ void MixedSolventElectrolyte::initThermoXML(XML_Node& phaseNode, const std::stri
"phasenode and Id are incompatible");
}
/*
* Check on the thermo field. Must have:
* <thermo model="MixedSolventElectrolyte" />
*/
// Check on the thermo field. Must have:
// <thermo model="MixedSolventElectrolyte" />
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("MixedSolventElectrolyte::initThermoXML",
"no thermo XML node");
@ -314,10 +280,8 @@ void MixedSolventElectrolyte::initThermoXML(XML_Node& phaseNode, const std::stri
"Unknown thermo model: " + mString);
}
/*
* Go get all of the coefficients and factors in the
* activityCoefficients XML block
*/
// Go get all of the coefficients and factors in the activityCoefficients
// XML block
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& acNode = thermoNode.child("activityCoefficients");
mString = acNode.attrib("model");
@ -327,20 +291,17 @@ void MixedSolventElectrolyte::initThermoXML(XML_Node& phaseNode, const std::stri
}
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.
*/
// 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 (lowercase(xmlACChild.name()) == "binaryneutralspeciesparameters") {
readXMLBinarySpecies(xmlACChild);
}
}
}
/*
* Go down the chain
*/
// Go down the chain
MolarityIonicVPSSTP::initThermoXML(phaseNode, id_);
}
@ -479,9 +440,7 @@ void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN() const
double T = temperature();
dlnActCoeffdlnN_.zero();
/*
* Loop over the activity coefficient gamma_k
*/
// Loop over the activity coefficient gamma_k
for (size_t iK = 0; iK < m_kk; iK++) {
for (size_t iM = 0; iM < m_kk; iM++) {
double XM = moleFractions_[iM];
@ -594,10 +553,9 @@ void MixedSolventElectrolyte::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
if (jName == "") {
throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies", "no speciesB attrib");
}
/*
* Find the index of the species in the current phase. It's not
* an error to not find the species
*/
// 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;
@ -623,13 +581,10 @@ void MixedSolventElectrolyte::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
for (size_t iChild = 0; iChild < xmLBinarySpecies.nChildren(); iChild++) {
XML_Node& xmlChild = xmLBinarySpecies.child(iChild);
string nodeName = lowercase(xmlChild.name());
/*
* Process the binary species interaction child elements
*/
// Process the binary species interaction child elements
if (nodeName == "excessenthalpy") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, vParams, true, "toSI", "excessEnthalpy");
if (vParams.size() != 2) {
throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies::excessEnthalpy for " + ispName
@ -641,9 +596,7 @@ void MixedSolventElectrolyte::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
}
if (nodeName == "excessentropy") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, vParams, true, "toSI", "excessEntropy");
if (vParams.size() != 2) {
throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies::excessEntropy for " + ispName
@ -655,9 +608,7 @@ void MixedSolventElectrolyte::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
}
if (nodeName == "excessvolume_enthalpy") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Enthalpy");
if (vParams.size() != 2) {
throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies::excessVolume_Enthalpy for " + ispName
@ -669,9 +620,7 @@ void MixedSolventElectrolyte::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
}
if (nodeName == "excessvolume_entropy") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Entropy");
if (vParams.size() != 2) {
throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies::excessVolume_Entropy for " + ispName

164
src/thermo/MixtureFugacityTP.cpp
View file

@ -41,14 +41,10 @@ MixtureFugacityTP::MixtureFugacityTP(const MixtureFugacityTP& b) :
MixtureFugacityTP& MixtureFugacityTP::operator=(const MixtureFugacityTP& b)
{
if (&b != this) {
/*
* Mostly, this is a passthrough to the underlying
* assignment operator for the ThermoPhase parent object.
*/
// Mostly, this is a passthrough to the underlying assignment operator
// for the ThermoPhase parent object.
ThermoPhase::operator=(b);
/*
* However, we have to handle data that we own.
*/
// However, we have to handle data that we own.
m_Pcurrent = b.m_Pcurrent;
moleFractions_ = b.moleFractions_;
iState_ = b.iState_;
@ -88,9 +84,7 @@ int MixtureFugacityTP::reportSolnBranchActual() const
return iState_;
}
/*
* ---- Partial Molar Properties of the Solution -----------------
*/
// ---- Partial Molar Properties of the Solution -----------------
void MixtureFugacityTP::getChemPotentials_RT(doublereal* muRT) const
{
@ -100,9 +94,7 @@ void MixtureFugacityTP::getChemPotentials_RT(doublereal* muRT) const
}
}
/*
* ----- Thermodynamic Values for the Species Standard States States ----
*/
// ----- Thermodynamic Values for the Species Standard States States ----
void MixtureFugacityTP::getStandardChemPotentials(doublereal* g) const
{
@ -178,10 +170,7 @@ void MixtureFugacityTP::getStandardVolumes(doublereal* vol) const
}
}
/*
* ----- Thermodynamic Values for the Species Reference States ----
*/
// ----- Thermodynamic Values for the Species Reference States ----
void MixtureFugacityTP::getEnthalpy_RT_ref(doublereal* hrt) const
{
@ -332,15 +321,11 @@ void MixtureFugacityTP::calcDensity()
void MixtureFugacityTP::setState_TP(doublereal t, doublereal pres)
{
/*
* A pretty tricky algorithm is needed here, due to problems involving
* standard states of real fluids. For those cases you need
* to combine the T and P specification for the standard state, or else
* you may venture into the forbidden zone, especially when nearing the
* triple point.
* Therefore, we need to do the standard state thermo calc with the
* (t, pres) combo.
*/
// A pretty tricky algorithm is needed here, due to problems involving
// standard states of real fluids. For those cases you need to combine the T
// and P specification for the standard state, or else you may venture into
// the forbidden zone, especially when nearing the triple point. Therefore,
// we need to do the standard state thermo calc with the (t, pres) combo.
getMoleFractions(moleFractions_.data());
Phase::setTemperature(t);
@ -479,10 +464,8 @@ doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa,
}
}
} else {
/*
* Assume the Gas phase initial guess, if nothing is
* specified to the routine
*/
// Assume the Gas phase initial guess, if nothing is specified to
// the routine
rhoguess = presPa * mmw / (GasConstant * TKelvin);
}
}
@ -490,15 +473,13 @@ doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa,
double molarVolBase = mmw / rhoguess;
double molarVolLast = molarVolBase;
double vc = mmw / critDensity();
/*
* molar volume of the spinodal at the current temperature and mole fractions. this will
* be updated as we go.
*/
// molar volume of the spinodal at the current temperature and mole
// fractions. this will be updated as we go.
double molarVolSpinodal = vc;
bool conv = false;
/*
* We start on one side of the vc and stick with that side
*/
// We start on one side of the vc and stick with that side
bool gasSide = molarVolBase > vc;
if (gasSide) {
molarVolLast = (GasConstant * TKelvin)/presPa;
@ -506,36 +487,27 @@ doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa,
molarVolLast = liquidVolEst(TKelvin, presPa);
}
/*
* OK, now we do a small solve to calculate the molar volume given the T,P value.
* The algorithm is taken from dfind()
*/
// OK, now we do a small solve to calculate the molar volume given the T,P
// value. The algorithm is taken from dfind()
for (int n = 0; n < 200; n++) {
/*
* Calculate the predicted reduced pressure, pred0, based on the
* current tau and dd.
* Calculate the derivative of the predicted pressure
* wrt the molar volume.
* This routine also returns the pressure, presBase
*/
// Calculate the predicted reduced pressure, pred0, based on the current
// tau and dd. Calculate the derivative of the predicted pressure wrt
// the molar volume. This routine also returns the pressure, presBase
double presBase;
double dpdVBase = dpdVCalc(TKelvin, molarVolBase, presBase);
/*
* If dpdV is positive, then we are in the middle of the
* 2 phase region and beyond the spinodal stability curve. We need to adjust
* the initial guess outwards and start a new iteration.
*/
// If dpdV is positive, then we are in the middle of the 2 phase region
// and beyond the spinodal stability curve. We need to adjust the
// initial guess outwards and start a new iteration.
if (dpdVBase >= 0.0) {
if (TKelvin > tcrit) {
throw CanteraError("MixtureFugacityTP::densityCalc",
"T > tcrit unexpectedly");
}
/*
* TODO Spawn a calculation for the value of the spinodal point that is
* very accurate. Answer the question as to whether a solution is
* possible on the current side of the vapor dome.
*/
// TODO Spawn a calculation for the value of the spinodal point that
// is very accurate. Answer the question as to whether a
// solution is possible on the current side of the vapor dome.
if (gasSide) {
if (molarVolBase >= vc) {
molarVolSpinodal = molarVolBase;
@ -554,34 +526,25 @@ doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa,
continue;
}
/*
* Check for convergence
*/
// Check for convergence
if (fabs(presBase-presPa) < 1.0E-30 + 1.0E-8 * presPa) {
conv = true;
break;
}
/*
* Dampen and crop the update
*/
// Dampen and crop the update
doublereal dpdV = dpdVBase;
if (n < 10) {
dpdV = dpdVBase * 1.5;
}
/*
* Formulate the update to the molar volume by
* Newton's method. Then, crop it to a max value
* of 0.1 times the current volume
*/
// Formulate the update to the molar volume by Newton's method. Then,
// crop it to a max value of 0.1 times the current volume
double delMV = - (presBase - presPa) / dpdV;
if ((!gasSide || delMV < 0.0) && fabs(delMV) > 0.2 * molarVolBase) {
delMV = delMV / fabs(delMV) * 0.2 * molarVolBase;
}
/*
* Only go 1/10 the way towards the spinodal at any one time.
*/
// Only go 1/10 the way towards the spinodal at any one time.
if (TKelvin < tcrit) {
if (gasSide) {
if (delMV < 0.0 && -delMV > 0.5 * (molarVolBase - molarVolSpinodal)) {
@ -593,9 +556,7 @@ doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa,
}
}
}
/*
* updated the molar volume value
*/
// updated the molar volume value
molarVolLast = molarVolBase;
molarVolBase += delMV;
@ -604,17 +565,13 @@ doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa,
break;
}
/*
* Check for negative molar volumes
*/
// Check for negative molar volumes
if (molarVolBase <= 0.0) {
molarVolBase = std::min(1.0E-30, fabs(delMV*1.0E-4));
}
}
/*
* Check for convergence, and return 0.0 if it wasn't achieved.
*/
// Check for convergence, and return 0.0 if it wasn't achieved.
double densBase = 0.0;
if (! conv) {
molarVolBase = 0.0;
@ -733,24 +690,23 @@ doublereal MixtureFugacityTP::satPressure(doublereal TKelvin)
doublereal MixtureFugacityTP::calculatePsat(doublereal TKelvin, doublereal& molarVolGas,
doublereal& molarVolLiquid)
{
/*
* The algorithm for this routine has undergone quite a bit of work. It probably needs more work.
* However, it seems now to be fairly robust.
* The key requirement is to find an initial pressure where both the liquid and the gas exist. This
* is not as easy as it sounds, and it gets exceedingly hard as the critical temperature is approached
* from below.
* Once we have this initial state, then we seek to equilibrate the Gibbs free energies of the
* gas and liquid and use the formula
*
* dp = VdG
*
* to create an update condition for deltaP using
*
* - (Gliq - Ggas) = (Vliq - Vgas) (deltaP)
*
* @TODO Suggestions for the future would be to switch it to an algorithm that uses the gas molar volume
* and the liquid molar volumes as the fundamental unknowns.
*/
// The algorithm for this routine has undergone quite a bit of work. It
// probably needs more work. However, it seems now to be fairly robust. The
// key requirement is to find an initial pressure where both the liquid and
// the gas exist. This is not as easy as it sounds, and it gets exceedingly
// hard as the critical temperature is approached from below. Once we have
// this initial state, then we seek to equilibrate the Gibbs free energies
// of the gas and liquid and use the formula
//
// dp = VdG
//
// to create an update condition for deltaP using
//
// - (Gliq - Ggas) = (Vliq - Vgas) (deltaP)
//
// @TODO Suggestions for the future would be to switch it to an algorithm
// that uses the gas molar volume and the liquid molar volumes as the
// fundamental unknowns.
// we need this because this is a non-const routine that is public
setTemperature(TKelvin);
@ -767,10 +723,8 @@ doublereal MixtureFugacityTP::calculatePsat(doublereal TKelvin, doublereal& mola
doublereal delGRT = 1.0E6;
doublereal liqGRT, gasGRT;
/*
* First part of the calculation involves finding a pressure at which the
* gas and the liquid state coexists.
*/
// First part of the calculation involves finding a pressure at which
// the gas and the liquid state coexists.
doublereal presLiquid = 0.;
doublereal presGas;
doublereal presBase = pres;
@ -864,9 +818,7 @@ doublereal MixtureFugacityTP::calculatePsat(doublereal TKelvin, doublereal& mola
double RhoGas = RhoGasGood;
double RhoLiquid = RhoLiquidGood;
/*
* Now that we have found a good pressure we can proceed with the algorithm.
*/
// Now that we have found a good pressure we can proceed with the algorithm.
for (int i = 0; i < 20; i++) {
int stab = corr0(TKelvin, pres, RhoLiquid, RhoGas, liqGRT, gasGRT);
if (stab == 0) {

82
src/thermo/MolalityVPSSTP.cpp
View file

@ -5,11 +5,10 @@
* (see \ref thermoprops
* and class \link Cantera::MolalityVPSSTP MolalityVPSSTP\endlink).
*
* Header file for a derived class of ThermoPhase that handles
* variable pressure standard state methods for calculating
* thermodynamic properties that are further based upon activities
* based on the molality scale. These include most of the methods for
* calculating liquid electrolyte thermodynamics.
* Header file for a derived class of ThermoPhase that handles variable pressure
* standard state methods for calculating thermodynamic properties that are
* further based upon activities based on the molality scale. These include
* most of the methods for calculating liquid electrolyte thermodynamics.
*/
/*
* Copyright (2005) Sandia Corporation. Under the terms of
@ -36,11 +35,9 @@ MolalityVPSSTP::MolalityVPSSTP() :
m_xmolSolventMIN(0.01),
m_Mnaught(18.01528E-3)
{
/*
* Change the default to be that charge neutrality in the
* phase is necessary condition for the proper specification
* of thermodynamic functions within the phase
*/
// Change the default to be that charge neutrality in the phase is necessary
// condition for the proper specification of thermodynamic functions within
// the phase
m_chargeNeutralityNecessary = true;
}
@ -75,9 +72,7 @@ ThermoPhase* MolalityVPSSTP::duplMyselfAsThermoPhase() const
return new MolalityVPSSTP(*this);
}
/*
* -------------- Utilities -------------------------------
*/
// -------------- Utilities -------------------------------
void MolalityVPSSTP::setpHScale(const int pHscaleType)
{
@ -165,24 +160,18 @@ void MolalityVPSSTP::setMolalities(const doublereal* const molal)
}
}
setMoleFractions(m_molalities.data());
/*
* Essentially we don't trust the input: We calculate
* the molalities from the mole fractions that we
* just obtained.
*/
// Essentially we don't trust the input: We calculate the molalities from
// the mole fractions that we just obtained.
calcMolalities();
}
void MolalityVPSSTP::setMolalitiesByName(const compositionMap& mMap)
{
/*
* HKM -> Might need to be more complicated here, setting
* neutrals so that the existing mole fractions are
* preserved.
*/
/*
* Get a vector of mole fractions
*/
// HKM -> Might need to be more complicated here, setting neutrals so that
// the existing mole fractions are preserved.
// Get a vector of mole fractions
vector_fp mf(m_kk, 0.0);
getMoleFractions(mf.data());
double xmolSmin = std::max(mf[m_indexSolvent], m_xmolSolventMIN);
@ -192,9 +181,8 @@ void MolalityVPSSTP::setMolalitiesByName(const compositionMap& mMap)
mf[k] = mol_k * m_Mnaught * xmolSmin;
}
}
/*
* check charge neutrality
*/
// check charge neutrality
size_t largePos = npos;
double cPos = 0.0;
size_t largeNeg = npos;
@ -240,11 +228,9 @@ void MolalityVPSSTP::setMolalitiesByName(const compositionMap& mMap)
mf[k] *= sum;
}
setMoleFractions(mf.data());
/*
* After we formally set the mole fractions, we
* calculate the molalities again and store it in
* this object.
*/
// After we formally set the mole fractions, we calculate the molalities
// again and store it in this object.
calcMolalities();
}
@ -254,9 +240,7 @@ void MolalityVPSSTP::setMolalitiesByName(const std::string& x)
setMolalitiesByName(xx);
}
/*
* - Activities, Standard States, Activity Concentrations -----------
*/
// - Activities, Standard States, Activity Concentrations -----------
int MolalityVPSSTP::activityConvention() const
{
@ -296,14 +280,11 @@ void MolalityVPSSTP::getMolalityActivityCoefficients(doublereal* acMolality) con
doublereal MolalityVPSSTP::osmoticCoefficient() const
{
/*
* First, we calculate the activities all over again
*/
// First, we calculate the activities all over again
vector_fp act(m_kk);
getActivities(act.data());
/*
* Then, we calculate the sum of the solvent molalities
*/
// Then, we calculate the sum of the solvent molalities
double sum = 0;
for (size_t k = 1; k < m_kk; k++) {
sum += std::max(m_molalities[k], 0.0);
@ -366,13 +347,11 @@ void MolalityVPSSTP::initThermo()
{
initLengths();
VPStandardStateTP::initThermo();
/*
* The solvent defaults to species 0
*/
// The solvent defaults to species 0
setSolvent(0);
/*
* Find the Cl- species
*/
// Find the Cl- species
m_indexCLM = findCLMIndex();
}
@ -452,9 +431,8 @@ void MolalityVPSSTP::initLengths()
void MolalityVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
initLengths();
/*
* The solvent defaults to species 0
*/
// The solvent defaults to species 0
setSolvent(0);
VPStandardStateTP::initThermoXML(phaseNode, id_);
}

113
src/thermo/MolarityIonicVPSSTP.cpp
View file

@ -5,11 +5,10 @@
* (see \ref thermoprops
* and class \link Cantera::MolarityIonicVPSSTP MolarityIonicVPSSTP\endlink).
*
* Header file for a derived class of ThermoPhase that handles
* variable pressure standard state methods for calculating
* thermodynamic properties that are further based upon expressions
* for the excess Gibbs free energy expressed as a function of
* the mole fractions.
* Header file for a derived class of ThermoPhase that handles variable pressure
* standard state methods for calculating thermodynamic properties that are
* further based upon expressions for the excess Gibbs free energy expressed as
* a function of the mole fractions.
*/
/*
* Copyright (2009) Sandia Corporation. Under the terms of
@ -88,20 +87,14 @@ ThermoPhase* MolarityIonicVPSSTP::duplMyselfAsThermoPhase() const
return new MolarityIonicVPSSTP(*this);
}
/*
* - Activities, Standard States, Activity Concentrations -----------
*/
// - Activities, Standard States, Activity Concentrations -----------
void MolarityIonicVPSSTP::getLnActivityCoefficients(doublereal* lnac) const
{
/*
* Update the activity coefficients
*/
// Update the activity coefficients
s_update_lnActCoeff();
/*
* take the exp of the internally stored coefficients.
*/
// take the exp of the internally stored coefficients.
for (size_t k = 0; k < m_kk; k++) {
lnac[k] = lnActCoeff_Scaled_[k];
}
@ -109,16 +102,11 @@ void MolarityIonicVPSSTP::getLnActivityCoefficients(doublereal* lnac) const
void MolarityIonicVPSSTP::getChemPotentials(doublereal* mu) const
{
/*
* First get the standard chemical potentials in
* molar form.
* -> this requires updates of standard state as a function
* of T and P
*/
// First get the standard chemical potentials in molar form. This requires
// updates of standard state as a function of T and P
getStandardChemPotentials(mu);
/*
* Update the activity coefficients
*/
// Update the activity coefficients
s_update_lnActCoeff();
for (size_t k = 0; k < m_kk; k++) {
double xx = std::max(moleFractions_[k], SmallNumber);
@ -137,21 +125,17 @@ void MolarityIonicVPSSTP::getElectrochemPotentials(doublereal* mu) const
void MolarityIonicVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const
{
/*
* Get the nondimensional standard state enthalpies
*/
// Get the nondimensional standard state enthalpies
getEnthalpy_RT(hbar);
/*
* dimensionalize it.
*/
// dimensionalize it.
double T = temperature();
for (size_t k = 0; k < m_kk; k++) {
hbar[k] *= GasConstant * T;
}
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnActCoeff();
s_update_dlnActCoeff_dT();
for (size_t k = 0; k < m_kk; k++) {
@ -161,24 +145,20 @@ void MolarityIonicVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const
void MolarityIonicVPSSTP::getPartialMolarCp(doublereal* cpbar) const
{
/*
* Get the nondimensional standard state entropies
*/
// Get the nondimensional standard state entropies
getCp_R(cpbar);
double T = temperature();
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnActCoeff();
s_update_dlnActCoeff_dT();
for (size_t k = 0; k < m_kk; k++) {
cpbar[k] -= 2 * T * dlnActCoeffdT_Scaled_[k] + T * T * d2lnActCoeffdT2_Scaled_[k];
}
/*
* dimensionalize it.
*/
// dimensionalize it.
for (size_t k = 0; k < m_kk; k++) {
cpbar[k] *= GasConstant;
}
@ -186,15 +166,12 @@ void MolarityIonicVPSSTP::getPartialMolarCp(doublereal* cpbar) const
void MolarityIonicVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
{
/*
* Get the nondimensional standard state entropies
*/
// Get the nondimensional standard state entropies
getEntropy_R(sbar);
double T = temperature();
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnActCoeff();
s_update_dlnActCoeff_dT();
@ -202,9 +179,8 @@ void MolarityIonicVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
double xx = std::max(moleFractions_[k], SmallNumber);
sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
}
/*
* dimensionalize it.
*/
// dimensionalize it.
for (size_t k = 0; k < m_kk; k++) {
sbar[k] *= GasConstant;
}
@ -212,9 +188,7 @@ void MolarityIonicVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
void MolarityIonicVPSSTP::getPartialMolarVolumes(doublereal* vbar) const
{
/*
* Get the standard state values in m^3 kmol-1
*/
// Get the standard state values in m^3 kmol-1
getStandardVolumes(vbar);
for (size_t iK = 0; iK < m_kk; iK++) {
vbar[iK] += 0.0;
@ -305,9 +279,8 @@ void MolarityIonicVPSSTP::initThermo()
{
GibbsExcessVPSSTP::initThermo();
initLengths();
/*
* Go find the list of cations and anions
*/
// Go find the list of cations and anions
cationList_.clear();
anionList_.clear();
passThroughList_.clear();
@ -346,11 +319,9 @@ void MolarityIonicVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string&
"phasenode and Id are incompatible");
}
/*
* Check on the thermo field. Must have one of:
* <thermo model="MolarityIonicVPSS" />
* <thermo model="MolarityIonicVPSSTP" />
*/
// Check on the thermo field. Must have one of:
// <thermo model="MolarityIonicVPSS" />
// <thermo model="MolarityIonicVPSSTP" />
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("MolarityIonicVPSSTP::initThermoXML",
"no thermo XML node");
@ -363,26 +334,20 @@ void MolarityIonicVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string&
"Unknown thermo model: " + mStringa + " - This object only knows \"MolarityIonicVPSSTP\" ");
}
/*
* Go get all of the coefficients and factors in the
* activityCoefficients XML block
*/
// Go get all of the coefficients and factors in the activityCoefficients
// XML block
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& acNode = thermoNode.child("activityCoefficients");
for (size_t i = 0; i < acNode.nChildren(); i++) {
XML_Node& xmlACChild = acNode.child(i);
/*
* Process a binary interaction
*/
// Process a binary interaction
if (lowercase(xmlACChild.name()) == "binaryneutralspeciesparameters") {
readXMLBinarySpecies(xmlACChild);
}
}
}
/*
* Go down the chain
*/
// Go down the chain
GibbsExcessVPSSTP::initThermoXML(phaseNode, id);
}

34
src/thermo/Mu0Poly.cpp
View file

@ -104,11 +104,10 @@ Mu0Poly* newMu0ThermoFromXML(const XML_Node& Mu0Node)
throw CanteraError("installMu0ThermoFromXML", "missing Mu0Values");
}
getFloatArray(*valNode_ptr, cValues, true, "actEnergy");
/*
* Check to see whether the Mu0's were input in a dimensionless
* form. If they were, then the assumed temperature needs to be
* adjusted from the assumed T = 273.15
*/
// Check to see whether the Mu0's were input in a dimensionless form. If
// they were, then the assumed temperature needs to be adjusted from the
// assumed T = 273.15
if (valNode_ptr->attrib("units") == "Dimensionless") {
dimensionlessMu0Values = true;
}
@ -127,9 +126,7 @@ Mu0Poly* newMu0ThermoFromXML(const XML_Node& Mu0Node)
throw CanteraError("installMu0ThermoFromXML", "numPoints inconsistent");
}
/*
* Fix up dimensionless Mu0 values if input
*/
// Fix up dimensionless Mu0 values if input
if (dimensionlessMu0Values) {
for (size_t i = 0; i < numPoints; i++) {
cValues[i] *= cTemperatures[i] / 273.15;
@ -158,19 +155,16 @@ void Mu0Poly::processCoeffs(const doublereal* coeffs)
m_numIntervals = nPoints - 1;
m_H298 = coeffs[1] / GasConstant;
size_t iT298 = 0;
/*
* Resize according to the number of points
*/
// Resize according to the number of points
m_t0_int.resize(nPoints);
m_h0_R_int.resize(nPoints);
m_s0_R_int.resize(nPoints);
m_cp0_R_int.resize(nPoints);
m_mu0_R_int.resize(nPoints);
/*
* Calculate the T298 interval and make sure that
* the temperatures are strictly monotonic.
* Also distribute the data into the internal arrays.
*/
// Calculate the T298 interval and make sure that the temperatures are
// strictly monotonic. Also distribute the data into the internal arrays.
bool ifound = false;
for (size_t i = 0, iindex = 2; i < nPoints; i++) {
double T1 = coeffs[iindex];
@ -191,9 +185,7 @@ void Mu0Poly::processCoeffs(const doublereal* coeffs)
"One temperature has to be 298.15");
}
/*
* Starting from the interval with T298, we go up
*/
// Starting from the interval with T298, we go up
m_h0_R_int[iT298] = m_H298;
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++) {
@ -209,9 +201,7 @@ void Mu0Poly::processCoeffs(const doublereal* coeffs)
m_cp0_R_int[i+1] = cpi;
}
/*
* Starting from the interval with T298, we go down
*/
// Starting from the interval with T298, we go down
if (iT298 != 0) {
m_h0_R_int[iT298] = m_H298;
m_s0_R_int[iT298] = - (m_mu0_R_int[iT298] - m_h0_R_int[iT298]) / m_t0_int[iT298];

12
src/thermo/PDSS.cpp
View file

@ -94,10 +94,8 @@ PDSS::PDSS(const PDSS& b) :
m_gss_RT_ptr(b.m_gss_RT_ptr),
m_Vss_ptr(b.m_Vss_ptr)
{
/*
* Use the assignment operator to do the brunt
* of the work for the copy constructor.
*/
// Use the assignment operator to do the brunt of the work for the copy
// constructor.
*this = b;
}
@ -114,10 +112,8 @@ PDSS& PDSS::operator=(const PDSS& b)
m_minTemp = b.m_minTemp;
m_maxTemp = b.m_maxTemp;
// Pointers which are zero, are properly assigned in the
// function, initAllPtrs(). which must be called after the
// assignment operation.
// Pointers which are zero, are properly assigned in the function,
// initAllPtrs(). which must be called after the assignment operation.
m_tp = 0;
m_vpssmgr_ptr = 0;
m_mw = b.m_mw;

13
src/thermo/PDSS_ConstVol.cpp
View file

@ -47,10 +47,8 @@ PDSS_ConstVol::PDSS_ConstVol(VPStandardStateTP* tp, size_t spindex,
PDSS_ConstVol::PDSS_ConstVol(const PDSS_ConstVol& b) :
PDSS(b)
{
/*
* Use the assignment operator to do the brunt
* of the work for the copy constructor.
*/
// Use the assignment operator to do the brunt of the work for the copy
// constructor.
*this = b;
}
@ -109,10 +107,9 @@ void PDSS_ConstVol::constructPDSSFile(VPStandardStateTP* tp, size_t spindex,
throw CanteraError("PDSS_ConstVol::initThermo","could not open "
+path+" for reading.");
}
/*
* The phase object automatically constructs an XML object.
* Use this object to store information.
*/
// The phase object automatically constructs an XML object. Use this object
// to store information.
XML_Node fxml;
fxml.build(fin);
XML_Node* fxml_phase = findXMLPhase(&fxml, id);

39
src/thermo/PDSS_HKFT.cpp
View file

@ -23,9 +23,7 @@ using namespace std;
namespace Cantera
{
/*
* Set the default to error exit if there is an input file inconsistency
*/
// Set the default to error exit if there is an input file inconsistency
int PDSS_HKFT::s_InputInconsistencyErrorExit = 1;
PDSS_HKFT::PDSS_HKFT(VPStandardStateTP* tp, size_t spindex) :
@ -147,10 +145,9 @@ PDSS_HKFT::PDSS_HKFT(const PDSS_HKFT& b) :
{
m_pdssType = cPDSS_MOLAL_HKFT;
m_presR_bar = OneAtm * 1.0E-5;
/*
* Use the assignment operator to do the brunt
* of the work for the copy constructor.
*/
// Use the assignment operator to do the brunt of the work for the copy
// constructor.
*this = b;
}
@ -159,15 +156,13 @@ PDSS_HKFT& PDSS_HKFT::operator=(const PDSS_HKFT& b)
if (&b == this) {
return *this;
}
/*
* Call the base class operator
*/
// Call the base class operator
PDSS::operator=(b);
//! Need to call initAllPtrs AFTER, to get the correct m_waterSS
// Need to call initAllPtrs AFTER, to get the correct m_waterSS
m_waterSS = 0;
m_densWaterSS = b.m_densWaterSS;
//! Need to call initAllPtrs AFTER, to get the correct m_waterProps
// Need to call initAllPtrs AFTER, to get the correct m_waterProps
m_born_coeff_j = b.m_born_coeff_j;
m_r_e_j = b.m_r_e_j;
m_deltaG_formation_tr_pr = b.m_deltaG_formation_tr_pr;
@ -396,9 +391,8 @@ void PDSS_HKFT::initThermo()
PDSS::initThermo();
m_waterSS = dynamic_cast<PDSS_Water*>(m_tp->providePDSS(0));
/*
* Section to initialize m_Z_pr_tr and m_Y_pr_tr
*/
// Section to initialize m_Z_pr_tr and m_Y_pr_tr
m_temp = 273.15 + 25.;
m_pres = OneAtm;
doublereal relepsilon = m_waterProps->relEpsilon(m_temp, m_pres, 0);
@ -413,7 +407,7 @@ void PDSS_HKFT::initThermo()
m_charge_j = m_tp->charge(m_spindex);
convertDGFormation();
//! Ok, we have mu. Let's check it against the input value
// Ok, we have mu. Let's check it against the input value
// of DH_F to see that we have some internal consistency
doublereal Hcalc = m_Mu0_tr_pr + 298.15 * (m_Entrop_tr_pr * 1.0E3 * 4.184);
doublereal DHjmol = m_deltaH_formation_tr_pr * 1.0E3 * 4.184;
@ -621,10 +615,9 @@ void PDSS_HKFT::constructPDSSFile(VPStandardStateTP* tp, size_t spindex,
throw CanteraError("PDSS_HKFT::initThermo","could not open "
+path+" for reading.");
}
/*
* The phase object automatically constructs an XML object.
* Use this object to store information.
*/
// The phase object automatically constructs an XML object. Use this object
// to store information.
XML_Node fxml;
fxml.build(fin);
XML_Node* fxml_phase = findXMLPhase(&fxml, id);
@ -691,7 +684,7 @@ doublereal PDSS_HKFT::deltaH() const
return deltaH_calgmol * 1.0E3 * 4.184;
}
#endif
//================================================================================================================
doublereal PDSS_HKFT::deltaG() const
{
doublereal pbar = m_pres * 1.0E-5;
@ -906,9 +899,7 @@ doublereal PDSS_HKFT::LookupGe(const std::string& elemName)
void PDSS_HKFT::convertDGFormation()
{
/*
* Ok let's get the element compositions and conversion factors.
*/
// Ok let's get the element compositions and conversion factors.
doublereal totalSum = 0.0;
for (size_t m = 0; m < m_tp->nElements(); m++) {
double na = m_tp->nAtoms(m_spindex, m);

12
src/thermo/PDSS_IdealGas.cpp
View file

@ -48,10 +48,8 @@ PDSS_IdealGas::PDSS_IdealGas(VPStandardStateTP* tp, size_t spindex, const XML_No
PDSS_IdealGas::PDSS_IdealGas(const PDSS_IdealGas& b) :
PDSS(b)
{
/*
* Use the assignment operator to do the brunt
* of the work for the copy constructor.
*/
// Use the assignment operator to do the brunt of the work for the copy
// constructor.
*this = b;
}
@ -90,11 +88,9 @@ void PDSS_IdealGas::constructPDSSFile(VPStandardStateTP* tp, size_t spindex,
throw CanteraError("PDSS_IdealGas::constructPDSSFile","could not open "
+path+" for reading.");
}
/*
* The phase object automatically constructs an XML object.
* Use this object to store information.
*/
// The phase object automatically constructs an XML object. Use this object
// to store information.
XML_Node fxml;
fxml.build(fin);
XML_Node* fxml_phase = findXMLPhase(&fxml, id);

26
src/thermo/PDSS_IonsFromNeutral.cpp
View file

@ -63,10 +63,8 @@ PDSS_IonsFromNeutral::PDSS_IonsFromNeutral(VPStandardStateTP* tp, size_t spindex
PDSS_IonsFromNeutral::PDSS_IonsFromNeutral(const PDSS_IonsFromNeutral& b) :
PDSS(b)
{
/*
* Use the assignment operator to do the brunt
* of the work for the copy constructor.
*/
// Use the assignment operator to do the brunt of the work for the copy
// constructor.
*this = b;
}
@ -78,11 +76,9 @@ PDSS_IonsFromNeutral& PDSS_IonsFromNeutral::operator=(const PDSS_IonsFromNeutral
PDSS::operator=(b);
/*
* The shallow pointer copy in the next step will be insufficient in most cases. However, its
* functionally the best we can do for this assignment operator. We fix up the pointer in the
* initAllPtrs() function.
*/
// The shallow pointer copy in the next step will be insufficient in most
// cases. However, its functionally the best we can do for this assignment
// operator. We fix up the pointer in the initAllPtrs() function.
neutralMoleculePhase_ = b.neutralMoleculePhase_;
numMult_ = b.numMult_;
@ -178,10 +174,9 @@ void PDSS_IonsFromNeutral::constructPDSSFile(VPStandardStateTP* tp, size_t spind
throw CanteraError("PDSS_IonsFromNeutral::constructPDSSFile","could not open "
+path+" for reading.");
}
/*
* The phase object automatically constructs an XML object.
* Use this object to store information.
*/
// The phase object automatically constructs an XML object. Use this object
// to store information.
XML_Node fxml;
fxml.build(fin);
XML_Node* fxml_phase = findXMLPhase(&fxml, id);
@ -340,9 +335,8 @@ doublereal PDSS_IonsFromNeutral::molarVolume_ref() const
doublereal PDSS_IonsFromNeutral::temperature() const
{
/*
* Obtain the temperature from the owning VPStandardStateTP object if you can.
*/
// Obtain the temperature from the owning VPStandardStateTP object if you
// can.
m_temp = m_vpssmgr_ptr->temperature();
return m_temp;
}

13
src/thermo/PDSS_SSVol.cpp
View file

@ -58,10 +58,8 @@ PDSS_SSVol::PDSS_SSVol(const PDSS_SSVol& b) :
volumeModel_(cSSVOLUME_CONSTANT),
m_constMolarVolume(-1.0)
{
/*
* Use the assignment operator to do the brunt
* of the work for the copy constructor.
*/
// Use the assignment operator to do the brunt of the work for the copy
// constructor.
*this = b;
}
@ -137,10 +135,9 @@ void PDSS_SSVol::constructPDSSFile(VPStandardStateTP* tp, size_t spindex,
throw CanteraError("PDSS_SSVol::initThermo","could not open "
+path+" for reading.");
}
/*
* The phase object automatically constructs an XML object.
* Use this object to store information.
*/
// The phase object automatically constructs an XML object. Use this object
// to store information.
XML_Node fxml;
fxml.build(fin);
XML_Node* fxml_phase = findXMLPhase(&fxml, id);

32
src/thermo/PDSS_Water.cpp
View file

@ -99,10 +99,8 @@ PDSS_Water::PDSS_Water(const PDSS_Water& b) :
m_verbose(b.m_verbose),
m_allowGasPhase(b.m_allowGasPhase)
{
/*
* Use the assignment operator to do the brunt
* of the work for the copy constructor.
*/
// Use the assignment operator to do the brunt of the work for the copy
// constructor.
*this = b;
}
@ -111,9 +109,7 @@ PDSS_Water& PDSS_Water::operator=(const PDSS_Water& b)
if (&b == this) {
return *this;
}
/*
* Call the base class operator
*/
// Call the base class operator
PDSS::operator=(b);
m_sub = b.m_sub;
@ -154,10 +150,9 @@ void PDSS_Water::constructPDSSFile(VPStandardStateTP* tp, int spindex,
throw CanteraError("PDSS_Water::initThermo","could not open "
+path+" for reading.");
}
/*
* The phase object automatically constructs an XML object.
* Use this object to store information.
*/
// The phase object automatically constructs an XML object. Use this object
// to store information.
XML_Node fxml;
fxml.build(fin);
XML_Node* fxml_phase = findXMLPhase(&fxml, id);
@ -171,15 +166,11 @@ void PDSS_Water::constructPDSSFile(VPStandardStateTP* tp, int spindex,
void PDSS_Water::constructSet()
{
/*
* Calculate the molecular weight.
* hard coded to Cantera's elements and Water.
*/
// Calculate the molecular weight. hard coded to Cantera's elements and
// Water.
m_mw = 2 * 1.00794 + 15.9994;
/*
* Set the baseline
*/
// Set the baseline
doublereal T = 298.15;
m_p0 = OneAtm;
doublereal presLow = 1.0E-2;
@ -202,10 +193,7 @@ void PDSS_Water::constructSet()
}
h = enthalpy_mole();
/*
* Set the initial state of the system to 298.15 K and
* 1 bar.
*/
// Set the initial state of the system to 298.15 K and 1 bar.
setTemperature(298.15);
m_dens = m_sub.density(298.15, OneAtm, WATER_LIQUID);
m_pres = OneAtm;

34
src/thermo/Phase.cpp
View file

@ -78,12 +78,10 @@ Phase& Phase::operator=(const Phase& right)
m_elementNames = right.m_elementNames;
m_entropy298 = right.m_entropy298;
m_elem_type = right.m_elem_type;
/*
* This is a little complicated. -> Because we delete m_xml
* in the destructor, we own m_xml completely, and we need
* to have our own individual copies of the XML data tree
* in each object
*/
// This is a little complicated. -> Because we delete m_xml in the
// destructor, we own m_xml completely, and we need to have our own
// individual copies of the XML data tree in each object
if (m_xml) {
XML_Node* rroot = &m_xml->root();
delete rroot;
@ -327,9 +325,7 @@ void Phase::setMoleFractions(const doublereal* const x)
{
// Use m_y as a temporary work vector for the non-negative mole fractions
doublereal norm = 0.0;
/*
* sum is calculated below as the unnormalized molecular weight
*/
// sum is calculated below as the unnormalized molecular weight
doublereal sum = 0;
for (size_t k = 0; k < m_kk; k++) {
double xk = std::max(x[k], 0.0); // Ignore negative mole fractions
@ -337,24 +333,22 @@ void Phase::setMoleFractions(const doublereal* const x)
norm += xk;
sum += m_molwts[k] * xk;
}
/*
* Set m_ym_ to the normalized mole fractions divided by the normalized mean molecular weight:
* m_ym_k = X_k / (sum_k X_k M_k)
*/
// Set m_ym_ to the normalized mole fractions divided by the normalized mean
// molecular weight:
// m_ym_k = X_k / (sum_k X_k M_k)
const doublereal invSum = 1.0/sum;
for (size_t k=0; k < m_kk; k++) {
m_ym[k] = m_y[k]*invSum;
}
/*
* Now set m_y to the normalized mass fractions
* m_y = X_k M_k / (sum_k X_k M_k)
*/
// Now set m_y to the normalized mass fractions:
// m_y = X_k M_k / (sum_k X_k M_k)
for (size_t k=0; k < m_kk; k++) {
m_y[k] = m_ym[k] * m_molwts[k];
}
/*
* Calculate the normalized molecular weight
*/
// Calculate the normalized molecular weight
m_mmw = sum/norm;
m_stateNum++;
}

177
src/thermo/PhaseCombo_Interaction.cpp
View file

@ -1,5 +1,5 @@
/**
* @file
* @file PhaseCombo_Interaction.cpp
*/
/*
* Copyright (2009) Sandia Corporation. Under the terms of
@ -80,37 +80,27 @@ ThermoPhase* PhaseCombo_Interaction::duplMyselfAsThermoPhase() const
return new PhaseCombo_Interaction(*this);
}
/*
* -------------- Utilities -------------------------------
*/
// -------------- Utilities -------------------------------
int PhaseCombo_Interaction::eosType() const
{
return cPhaseCombo_Interaction;
}
/*
* - Activities, Standard States, Activity Concentrations -----------
*/
// - Activities, Standard States, Activity Concentrations -----------
void PhaseCombo_Interaction::getActivityCoefficients(doublereal* ac) const
{
/*
* Update the activity coefficients
*/
// Update the activity coefficients
s_update_lnActCoeff();
/*
* take the exp of the internally stored coefficients.
*/
// take the exp of the internally stored coefficients.
for (size_t k = 0; k < m_kk; k++) {
ac[k] = exp(lnActCoeff_Scaled_[k]);
}
}
/*
* ------------ Partial Molar Properties of the Solution ------------
*/
// ------------ Partial Molar Properties of the Solution ------------
void PhaseCombo_Interaction::getElectrochemPotentials(doublereal* mu) const
{
@ -123,16 +113,10 @@ void PhaseCombo_Interaction::getElectrochemPotentials(doublereal* mu) const
void PhaseCombo_Interaction::getChemPotentials(doublereal* mu) const
{
/*
* First get the standard chemical potentials in
* molar form.
* -> this requires updates of standard state as a function
* of T and P
*/
// First get the standard chemical potentials in molar form. This requires
// updates of standard state as a function of T and P
getStandardChemPotentials(mu);
/*
* Update the activity coefficients
*/
// Update the activity coefficients
s_update_lnActCoeff();
for (size_t k = 0; k < m_kk; k++) {
@ -181,21 +165,16 @@ doublereal PhaseCombo_Interaction::cv_mole() const
void PhaseCombo_Interaction::getPartialMolarEnthalpies(doublereal* hbar) const
{
/*
* Get the nondimensional standard state enthalpies
*/
// Get the nondimensional standard state enthalpies
getEnthalpy_RT(hbar);
/*
* dimensionalize it.
*/
// dimensionalize it.
double T = temperature();
for (size_t k = 0; k < m_kk; k++) {
hbar[k] *= GasConstant * T;
}
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnActCoeff();
s_update_dlnActCoeff_dT();
for (size_t k = 0; k < m_kk; k++) {
@ -205,24 +184,20 @@ void PhaseCombo_Interaction::getPartialMolarEnthalpies(doublereal* hbar) const
void PhaseCombo_Interaction::getPartialMolarCp(doublereal* cpbar) const
{
/*
* Get the nondimensional standard state entropies
*/
// Get the nondimensional standard state entropies
getCp_R(cpbar);
double T = temperature();
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnActCoeff();
s_update_dlnActCoeff_dT();
for (size_t k = 0; k < m_kk; k++) {
cpbar[k] -= 2 * T * dlnActCoeffdT_Scaled_[k] + T * T * d2lnActCoeffdT2_Scaled_[k];
}
/*
* dimensionalize it.
*/
// dimensionalize it.
for (size_t k = 0; k < m_kk; k++) {
cpbar[k] *= GasConstant;
}
@ -230,15 +205,12 @@ void PhaseCombo_Interaction::getPartialMolarCp(doublereal* cpbar) const
void PhaseCombo_Interaction::getPartialMolarEntropies(doublereal* sbar) const
{
/*
* Get the nondimensional standard state entropies
*/
// Get the nondimensional standard state entropies
getEntropy_R(sbar);
double T = temperature();
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnActCoeff();
s_update_dlnActCoeff_dT();
@ -246,9 +218,8 @@ void PhaseCombo_Interaction::getPartialMolarEntropies(doublereal* sbar) const
double xx = std::max(moleFractions_[k], SmallNumber);
sbar[k] += - lnActCoeff_Scaled_[k] - log(xx) - T * dlnActCoeffdT_Scaled_[k];
}
/*
* dimensionalize it.
*/
// dimensionalize it.
for (size_t k = 0; k < m_kk; k++) {
sbar[k] *= GasConstant;
}
@ -258,9 +229,7 @@ void PhaseCombo_Interaction::getPartialMolarVolumes(doublereal* vbar) const
{
double T = temperature();
/*
* Get the standard state values in m^3 kmol-1
*/
// Get the standard state values in m^3 kmol-1
getStandardVolumes(vbar);
for (size_t iK = 0; iK < m_kk; iK++) {
@ -303,10 +272,8 @@ void PhaseCombo_Interaction::initThermoXML(XML_Node& phaseNode, const std::strin
"phasenode and Id are incompatible");
}
/*
* Check on the thermo field. Must have:
* <thermo model="PhaseCombo_Interaction" />
*/
// Check on the thermo field. Must have:
// <thermo model="PhaseCombo_Interaction" />
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("PhaseCombo_Interaction::initThermoXML",
"no thermo XML node");
@ -318,10 +285,8 @@ void PhaseCombo_Interaction::initThermoXML(XML_Node& phaseNode, const std::strin
"model name isn't PhaseCombo_Interaction: " + formString);
}
/*
* Go get all of the coefficients and factors in the
* activityCoefficients XML block
*/
// Go get all of the coefficients and factors in the activityCoefficients
// XML block
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& acNode = thermoNode.child("activityCoefficients");
string mString = acNode.attrib("model");
@ -331,20 +296,17 @@ void PhaseCombo_Interaction::initThermoXML(XML_Node& phaseNode, const std::strin
}
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.
*/
// 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 (lowercase(xmlACChild.name()) == "binaryneutralspeciesparameters") {
readXMLBinarySpecies(xmlACChild);
}
}
}
/*
* Go down the chain
*/
// Go down the chain
GibbsExcessVPSSTP::initThermoXML(phaseNode, id);
}
@ -354,18 +316,13 @@ void PhaseCombo_Interaction::s_update_lnActCoeff() const
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
*/
// We never sample the end of the mole fraction domains
double xx = std::max(moleFractions_[iK], SmallNumber);
/*
* First wipe out the ideal solution mixing term
*/
// First wipe out the ideal solution mixing term
lnActCoeff_Scaled_[iK] = - log(xx);
/*
* Then add in the Margules interaction terms. that's it!
*/
// Then add in the Margules interaction terms. that's it!
for (size_t i = 0; i < numBinaryInteractions_; i++) {
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
@ -435,13 +392,10 @@ void PhaseCombo_Interaction::getdlnActCoeffds(const doublereal dTds, const doubl
s_update_dlnActCoeff_dT();
for (size_t iK = 0; iK < m_kk; iK++) {
/*
* We never sample the end of the mole fraction domains
*/
// We never sample the end of the mole fraction domains
double xx = std::max(moleFractions_[iK], SmallNumber);
/*
* First wipe out the ideal solution mixing term
*/
// First wipe out the ideal solution mixing term
if (xx > SmallNumber) {
dlnActCoeffds[iK] += - 1.0 / xx;
}
@ -477,13 +431,10 @@ void PhaseCombo_Interaction::s_update_dlnActCoeff_dlnN_diag() const
for (size_t iK = 0; iK < m_kk; iK++) {
double XK = moleFractions_[iK];
/*
* We never sample the end of the mole fraction domains
*/
// We never sample the end of the mole fraction domains
double xx = std::max(moleFractions_[iK], SmallNumber);
/*
* First wipe out the ideal solution mixing term
*/
// First wipe out the ideal solution mixing term
if (xx > SmallNumber) {
dlnActCoeffdlnN_diag_[iK] = - 1.0 + xx;
}
@ -515,13 +466,9 @@ void PhaseCombo_Interaction::s_update_dlnActCoeff_dlnN() const
double T = temperature();
dlnActCoeffdlnN_.zero();
/*
* Loop over the activity coefficient gamma_k
*/
// Loop over the activity coefficient gamma_k
for (size_t iK = 0; iK < m_kk; iK++) {
/*
* We never sample the end of the mole fraction domains
*/
// We never sample the end of the mole fraction domains
double xx = std::max(moleFractions_[iK], SmallNumber);
for (size_t iM = 0; iM < m_kk; iM++) {
double XM = moleFractions_[iM];
@ -646,10 +593,9 @@ void PhaseCombo_Interaction::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
if (jName == "") {
throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies", "no speciesB attrib");
}
/*
* Find the index of the species in the current phase. It's not
* an error to not find the species
*/
// 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;
@ -675,13 +621,10 @@ void PhaseCombo_Interaction::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
for (size_t iChild = 0; iChild < xmLBinarySpecies.nChildren(); iChild++) {
XML_Node& xmlChild = xmLBinarySpecies.child(iChild);
string nodeName = lowercase(xmlChild.name());
/*
* Process the binary species interaction child elements
*/
// Process the binary species interaction child elements
if (nodeName == "excessenthalpy") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, vParams, true, "toSI", "excessEnthalpy");
if (vParams.size() != 2) {
throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies::excessEnthalpy for " + ispName
@ -693,9 +636,7 @@ void PhaseCombo_Interaction::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
}
if (nodeName == "excessentropy") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, vParams, true, "toSI", "excessEntropy");
if (vParams.size() != 2) {
throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies::excessEntropy for " + ispName
@ -707,9 +648,7 @@ void PhaseCombo_Interaction::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
}
if (nodeName == "excessvolume_enthalpy") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Enthalpy");
if (vParams.size() != 2) {
throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies::excessVolume_Enthalpy for " + ispName
@ -721,9 +660,7 @@ void PhaseCombo_Interaction::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
}
if (nodeName == "excessvolume_entropy") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Entropy");
if (vParams.size() != 2) {
throw CanteraError("PhaseCombo_Interaction::readXMLBinarySpecies::excessVolume_Entropy for " + ispName

136
src/thermo/RedlichKisterVPSSTP.cpp
View file

@ -79,28 +79,19 @@ ThermoPhase* RedlichKisterVPSSTP::duplMyselfAsThermoPhase() const
return new RedlichKisterVPSSTP(*this);
}
/*
* - Activities, Standard States, Activity Concentrations -----------
*/
// - Activities, Standard States, Activity Concentrations -----------
void RedlichKisterVPSSTP::getLnActivityCoefficients(doublereal* lnac) const
{
/*
* Update the activity coefficients
*/
// Update the activity coefficients
s_update_lnActCoeff();
/*
* take the exp of the internally stored coefficients.
*/
for (size_t k = 0; k < m_kk; k++) {
lnac[k] = lnActCoeff_Scaled_[k];
}
}
/*
* ------------ Partial Molar Properties of the Solution ------------
*/
// ------------ Partial Molar Properties of the Solution ------------
void RedlichKisterVPSSTP::getElectrochemPotentials(doublereal* mu) const
{
@ -113,16 +104,10 @@ void RedlichKisterVPSSTP::getElectrochemPotentials(doublereal* mu) const
void RedlichKisterVPSSTP::getChemPotentials(doublereal* mu) const
{
/*
* First get the standard chemical potentials in
* molar form.
* -> this requires updates of standard state as a function
* of T and P
*/
// First get the standard chemical potentials in molar form. This requires
// updates of standard state as a function of T and P
getStandardChemPotentials(mu);
/*
* Update the activity coefficients
*/
// Update the activity coefficients
s_update_lnActCoeff();
for (size_t k = 0; k < m_kk; k++) {
@ -171,21 +156,16 @@ doublereal RedlichKisterVPSSTP::cv_mole() const
void RedlichKisterVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const
{
/*
* Get the nondimensional standard state enthalpies
*/
// Get the nondimensional standard state enthalpies
getEnthalpy_RT(hbar);
/*
* dimensionalize it.
*/
// dimensionalize it.
double T = temperature();
for (size_t k = 0; k < m_kk; k++) {
hbar[k] *= GasConstant * T;
}
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnActCoeff();
s_update_dlnActCoeff_dT();
for (size_t k = 0; k < m_kk; k++) {
@ -195,24 +175,18 @@ void RedlichKisterVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const
void RedlichKisterVPSSTP::getPartialMolarCp(doublereal* cpbar) const
{
/*
* Get the nondimensional standard state entropies
*/
getCp_R(cpbar);
double T = temperature();
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnActCoeff();
s_update_dlnActCoeff_dT();
for (size_t k = 0; k < m_kk; k++) {
cpbar[k] -= 2 * T * dlnActCoeffdT_Scaled_[k] + T * T * d2lnActCoeffdT2_Scaled_[k];
}
/*
* dimensionalize it.
*/
// dimensionalize it.
for (size_t k = 0; k < m_kk; k++) {
cpbar[k] *= GasConstant;
}
@ -220,15 +194,12 @@ void RedlichKisterVPSSTP::getPartialMolarCp(doublereal* cpbar) const
void RedlichKisterVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
{
/*
* Get the nondimensional standard state entropies
*/
// Get the nondimensional standard state entropies
getEntropy_R(sbar);
double T = temperature();
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
// Update the activity coefficients, This also update the internally stored
// molalities.
s_update_lnActCoeff();
s_update_dlnActCoeff_dT();
@ -236,9 +207,7 @@ void RedlichKisterVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
double xx = std::max(moleFractions_[k], SmallNumber);
sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
}
/*
* dimensionalize it.
*/
// dimensionalize it.
for (size_t k = 0; k < m_kk; k++) {
sbar[k] *= GasConstant;
}
@ -246,9 +215,7 @@ void RedlichKisterVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
void RedlichKisterVPSSTP::getPartialMolarVolumes(doublereal* vbar) const
{
/*
* Get the standard state values in m^3 kmol-1
*/
// Get the standard state values in m^3 kmol-1
getStandardVolumes(vbar);
for (size_t iK = 0; iK < m_kk; iK++) {
vbar[iK] += 0.0;
@ -273,10 +240,8 @@ void RedlichKisterVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string&
"phasenode and Id are incompatible");
}
/*
* Check on the thermo field. Must have:
* <thermo model="Redlich-Kister" />
*/
// Check on the thermo field. Must have:
// <thermo model="Redlich-Kister" />
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("RedlichKisterVPSSTP::initThermoXML",
"no thermo XML node");
@ -288,10 +253,8 @@ void RedlichKisterVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string&
"Unknown thermo model: " + mString + " - This object only knows \"Redlich-Kister\" ");
}
/*
* Go get all of the coefficients and factors in the
* activityCoefficients XML block
*/
// Go get all of the coefficients and factors in the activityCoefficients
// XML block
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& acNode = thermoNode.child("activityCoefficients");
mString = acNode.attrib("model");
@ -301,19 +264,16 @@ void RedlichKisterVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string&
}
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.
*/
// 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 (lowercase(xmlACChild.name()) == "binaryneutralspeciesparameters") {
readXMLBinarySpecies(xmlACChild);
}
}
}
/*
* Go down the chain
*/
// Go down the chain
GibbsExcessVPSSTP::initThermoXML(phaseNode, id_);
}
@ -322,11 +282,10 @@ void RedlichKisterVPSSTP::s_update_lnActCoeff() const
doublereal T = temperature();
lnActCoeff_Scaled_.assign(m_kk, 0.0);
/*
* Scaling: I moved the division of RT higher so that we are always dealing with G/RT dimensionless terms
* within the routine. There is a severe problem with roundoff error in these calculations. The
* dimensionless terms help.
*/
// Scaling: I moved the division of RT higher so that we are always dealing
// with G/RT dimensionless terms within the routine. There is a severe
// problem with roundoff error in these calculations. The dimensionless
// terms help.
for (size_t i = 0; i < numBinaryInteractions_; i++) {
size_t iA = m_pSpecies_A_ij[i];
size_t iB = m_pSpecies_B_ij[i];
@ -558,11 +517,10 @@ void RedlichKisterVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
if (jName == "") {
throw CanteraError("RedlichKisterVPSSTP::readXMLBinarySpecies", "no speciesB attrib");
}
/*
* Find the index of the species in the current phase. It's not
* an error to not find the species. This means that the interaction doesn't occur for the current
* implementation of the phase.
*/
// Find the index of the species in the current phase. It's not an error to
// not find the species. This means that the interaction doesn't occur for
// the current implementation of the phase.
size_t iSpecies = speciesIndex(iName);
if (iSpecies == npos) {
return;
@ -579,9 +537,8 @@ void RedlichKisterVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
if (charge(jSpecies) != 0) {
throw CanteraError("RedlichKisterVPSSTP::readXMLBinarySpecies", "speciesB charge problem");
}
/*
* Ok we have found a valid interaction
*/
// Ok we have found a valid interaction
numBinaryInteractions_++;
size_t iSpot = numBinaryInteractions_ - 1;
m_pSpecies_A_ij.resize(numBinaryInteractions_);
@ -592,21 +549,16 @@ void RedlichKisterVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
for (size_t iChild = 0; iChild < xmLBinarySpecies.nChildren(); iChild++) {
XML_Node& xmlChild = xmLBinarySpecies.child(iChild);
string nodeName = lowercase(xmlChild.name());
/*
* Process the binary species interaction child elements
*/
// Process the binary species interaction child elements
if (nodeName == "excessenthalpy") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, hParams, true, "toSI", "excessEnthalpy");
Npoly = std::max(hParams.size(), Npoly);
}
if (nodeName == "excessentropy") {
/*
* Get the string containing all of the values
*/
// Get the string containing all of the values
getFloatArray(xmlChild, sParams, true, "toSI", "excessEntropy");
Npoly = std::max(sParams.size(), Npoly);
}

129
src/thermo/RedlichKwongMFTP.cpp
View file

@ -78,14 +78,11 @@ RedlichKwongMFTP::RedlichKwongMFTP(const RedlichKwongMFTP& b) :
RedlichKwongMFTP& RedlichKwongMFTP::operator=(const RedlichKwongMFTP& b)
{
if (&b != this) {
/*
* Mostly, this is a passthrough to the underlying
* assignment operator for the ThermoPhae parent object.
*/
// Mostly, this is a passthrough to the underlying assignment operator
// for the ThermoPhae parent object.
MixtureFugacityTP::operator=(b);
/*
* However, we have to handle data that we own.
*/
// However, we have to handle data that we own.
m_standardMixingRules = b.m_standardMixingRules;
m_formTempParam = b.m_formTempParam;
m_b_current = b.m_b_current;
@ -121,9 +118,7 @@ int RedlichKwongMFTP::eosType() const
return cRedlichKwongMFTP;
}
/*
* ------------Molar Thermodynamic Properties -------------------------
*/
// ------------Molar Thermodynamic Properties -------------------------
doublereal RedlichKwongMFTP::enthalpy_mole() const
{
@ -190,16 +185,13 @@ doublereal RedlichKwongMFTP::pressure() const
void RedlichKwongMFTP::calcDensity()
{
/*
* Calculate the molarVolume of the solution (m**3 kmol-1)
*/
// Calculate the molarVolume of the solution (m**3 kmol-1)
const doublereal* const dtmp = moleFractdivMMW();
getPartialMolarVolumes(m_tmpV.data());
double invDens = dot(m_tmpV.begin(), m_tmpV.end(), dtmp);
/*
* Set the density in the parent State object directly,
* by calling the Phase::setDensity() function.
*/
// Set the density in the parent State object directly, by calling the
// Phase::setDensity() function.
Phase::setDensity(1.0/invDens);
}
@ -284,9 +276,7 @@ void RedlichKwongMFTP::getActivityCoefficients(doublereal* ac) const
}
}
/*
* ---- Partial Molar Properties of the Solution -----------------
*/
// ---- Partial Molar Properties of the Solution -----------------
void RedlichKwongMFTP::getChemPotentials_RT(doublereal* muRT) const
{
@ -332,15 +322,11 @@ void RedlichKwongMFTP::getChemPotentials(doublereal* mu) const
void RedlichKwongMFTP::getPartialMolarEnthalpies(doublereal* hbar) const
{
/*
* First we get the reference state contributions
*/
// First we get the reference state contributions
getEnthalpy_RT_ref(hbar);
scale(hbar, hbar+m_kk, hbar, RT());
/*
* We calculate dpdni_
*/
// We calculate dpdni_
doublereal TKelvin = temperature();
doublereal mv = molarVolume();
doublereal sqt = sqrt(TKelvin);
@ -572,14 +558,11 @@ void RedlichKwongMFTP::setToEquilState(const doublereal* mu_RT)
_updateReferenceStateThermo();
getGibbs_RT_ref(m_tmpV.data());
/*
* Within the method, we protect against inf results if the
* exponent is too high.
*
* If it is too low, we set
* the partial pressure to zero. This capability is needed
* by the elemental potential method.
*/
// Within the method, we protect against inf results if the exponent is too
// high.
//
// If it is too low, we set the partial pressure to zero. This capability is
// needed by the elemental potential method.
doublereal pres = 0.0;
double m_p0 = refPressure();
for (size_t k = 0; k < m_kk; k++) {
@ -620,12 +603,10 @@ void RedlichKwongMFTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
{
RedlichKwongMFTP::initLengths();
/*
* Check the model parameter for the Redlich-Kwong equation of state
* two are allowed
* RedlichKwong mixture of species, each of which are RK fluids
* RedlichKwongMFTP mixture of species with cross term coefficients
*/
// Check the model parameter for the Redlich-Kwong equation of state
// two are allowed
// RedlichKwong mixture of species, each of which are RK fluids
// RedlichKwongMFTP mixture of species with cross term coefficients
if (phaseNode.hasChild("thermo")) {
XML_Node& thermoNode = phaseNode.child("thermo");
std::string model = thermoNode["model"];
@ -638,29 +619,20 @@ void RedlichKwongMFTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
"Unknown thermo model : " + model);
}
/*
* Go get all of the coefficients and factors in the
* activityCoefficients XML block
*/
// 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;
size_t nC = acNode.nChildren();
/*
* Loop through the children getting multiple instances of
* parameters
*/
// Loop through the children getting multiple instances of
// parameters
for (size_t i = 0; i < nC; i++) {
XML_Node& xmlACChild = acNodePtr->child(i);
string 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 == "purefluidparameters") {
readXMLPureFluid(xmlACChild);
}
@ -668,19 +640,13 @@ void RedlichKwongMFTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
if (m_standardMixingRules == 1) {
applyStandardMixingRules();
}
/*
* Loop through the children getting multiple instances of
* parameters
*/
// Loop through the children getting multiple instances of
// parameters
for (size_t i = 0; i < nC; i++) {
XML_Node& xmlACChild = acNodePtr->child(i);
string 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 == "crossfluidparameters") {
readXMLCrossFluid(xmlACChild);
}
@ -708,10 +674,8 @@ void RedlichKwongMFTP::readXMLPureFluid(XML_Node& pureFluidParam)
"Incorrect name for processing this routine: " + xname);
}
/*
* Read the species
* Find the index of the species in the current phase. It's not an error to not find the species
*/
// Read the species. Find the index of the species in the current phase.
// It's not an error to not find the species
string iName = pureFluidParam.attrib("species");
if (iName == "") {
throw CanteraError("RedlichKwongMFTP::readXMLPureFluid", "no species attribute");
@ -789,10 +753,8 @@ void RedlichKwongMFTP::readXMLCrossFluid(XML_Node& CrossFluidParam)
"Incorrect name for processing this routine: " + xname);
}
/*
* Read the species
* Find the index of the species in the current phase. It's not an error to not find the species
*/
// Read the species. Find the index of the species in the current phase.
// It's not an error to not find the species
string iName = CrossFluidParam.attrib("species1");
if (iName == "") {
throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid", "no species1 attribute");
@ -918,9 +880,7 @@ doublereal RedlichKwongMFTP::liquidVolEst(doublereal TKelvin, doublereal& presGu
doublereal RedlichKwongMFTP::densityCalc(doublereal TKelvin, doublereal presPa, int phaseRequested, doublereal rhoguess)
{
/*
* It's necessary to set the temperature so that m_a_current is set correctly.
*/
// It's necessary to set the temperature so that m_a_current is set correctly.
setTemperature(TKelvin);
double tcrit = critTemperature();
doublereal mmw = meanMolecularWeight();
@ -937,10 +897,8 @@ doublereal RedlichKwongMFTP::densityCalc(doublereal TKelvin, doublereal presPa,
}
}
} else {
/*
* Assume the Gas phase initial guess, if nothing is
* specified to the routine
*/
// Assume the Gas phase initial guess, if nothing is specified to
// the routine
rhoguess = presPa * mmw / (GasConstant * TKelvin);
}
}
@ -1179,9 +1137,8 @@ int RedlichKwongMFTP::NicholsSolve(double TKelvin, double pres, doublereal a, do
if (TKelvin <= 0.0) {
throw CanteraError("RedlichKwongMFTP::NicholsSolve()", "neg temperature");
}
/*
* Derive the coefficients of the cubic polynomial to solve.
*/
// Derive the coefficients of the cubic polynomial to solve.
doublereal an = 1.0;
doublereal bn = - GasConstant * TKelvin / pres;
doublereal sqt = sqrt(TKelvin);
@ -1196,7 +1153,8 @@ int RedlichKwongMFTP::NicholsSolve(double TKelvin, double pres, doublereal a, do
// Derive the center of the cubic, x_N
doublereal xN = - bn /(3 * an);
// Derive the value of delta**2. This is a key quantity that determines the number of turning points
// Derive the value of delta**2. This is a key quantity that determines the
// number of turning points
doublereal delta2 = (bn * bn - 3 * an * cn) / (9 * an * an);
doublereal delta = 0.0;
@ -1248,9 +1206,7 @@ int RedlichKwongMFTP::NicholsSolve(double TKelvin, double pres, doublereal a, do
nSolnValues = 1;
}
/*
* One real root -> have to determine whether gas or liquid is the root
*/
// One real root -> have to determine whether gas or liquid is the root
if (desc > 0.0) {
doublereal tmpD = sqrt(desc);
doublereal tmp1 = (- yN + tmpD) / (2.0 * an);
@ -1325,9 +1281,8 @@ int RedlichKwongMFTP::NicholsSolve(double TKelvin, double pres, doublereal a, do
}
}
/*
* Unfortunately, there is a heavy amount of roundoff error due to bad conditioning in this
*/
// Unfortunately, there is a heavy amount of roundoff error due to bad
// conditioning in this
double res, dresdV = 0.0;
for (int i = 0; i < nSolnValues; i++) {
for (int n = 0; n < 20; n++) {

69
src/thermo/SingleSpeciesTP.cpp
View file

@ -54,9 +54,7 @@ int SingleSpeciesTP::eosType() const
throw NotImplementedError("SingleSpeciesTP::eosType");
}
/*
* ------------ Molar Thermodynamic Properties --------------------
*/
// ------------ Molar Thermodynamic Properties --------------------
doublereal SingleSpeciesTP::enthalpy_mole() const
{
@ -82,11 +80,9 @@ doublereal SingleSpeciesTP::entropy_mole() const
doublereal SingleSpeciesTP::gibbs_mole() const
{
double gbar;
/*
* Get the chemical potential of the first species.
* This is the same as the partial molar Gibbs
* free energy.
*/
// Get the chemical potential of the first species. This is the same as the
// partial molar Gibbs free energy.
getChemPotentials(&gbar);
return gbar;
}
@ -94,11 +90,9 @@ doublereal SingleSpeciesTP::gibbs_mole() const
doublereal SingleSpeciesTP::cp_mole() const
{
double cpbar;
/*
* Really should have a partial molar heat capacity
* function in ThermoPhase. However, the standard
* state heat capacity will do fine here for now.
*/
// Really should have a partial molar heat capacity function in ThermoPhase.
// However, the standard state heat capacity will do fine here for now.
getCp_R(&cpbar);
cpbar *= GasConstant;
return cpbar;
@ -106,15 +100,13 @@ doublereal SingleSpeciesTP::cp_mole() const
doublereal SingleSpeciesTP::cv_mole() const
{
/*
* For single species, we go directory to the general Cp - Cv relation
*
* Cp = Cv + alpha**2 * V * T / beta
*
* where
* alpha = volume thermal expansion coefficient
* beta = isothermal compressibility
*/
// For single species, we go directory to the general Cp - Cv relation
//
// Cp = Cv + alpha**2 * V * T / beta
//
// where
// alpha = volume thermal expansion coefficient
// beta = isothermal compressibility
doublereal cvbar = cp_mole();
doublereal alpha = thermalExpansionCoeff();
doublereal beta = isothermalCompressibility();
@ -126,9 +118,7 @@ doublereal SingleSpeciesTP::cv_mole() const
return cvbar;
}
/*
* ----------- Partial Molar Properties of the Solution -----------------
*/
// ----------- Partial Molar Properties of the Solution -----------------
void SingleSpeciesTP::getChemPotentials(doublereal* mu) const
{
@ -175,9 +165,7 @@ void SingleSpeciesTP::getPartialMolarVolumes(doublereal* vbar) const
vbar[0] = molecularWeight(0) / density();
}
/*
* Properties of the Standard State of the Species in the Solution
*/
// Properties of the Standard State of the Species in the Solution
void SingleSpeciesTP::getPureGibbs(doublereal* gpure) const
{
@ -190,9 +178,7 @@ void SingleSpeciesTP::getStandardVolumes(doublereal* vbar) const
vbar[0] = molecularWeight(0) / density();
}
/*
* ---- Thermodynamic Values for the Species Reference States -------
*/
// ---- Thermodynamic Values for the Species Reference States -------
void SingleSpeciesTP::getEnthalpy_RT_ref(doublereal* hrt) const
{
@ -224,9 +210,7 @@ void SingleSpeciesTP::getCp_R_ref(doublereal* cpr) const
cpr[0] = m_cp0_R[0];
}
/*
* ------------------ Setting the State ------------------------
*/
// ------------------ Setting the State ------------------------
void SingleSpeciesTP::setState_HP(doublereal h, doublereal p,
doublereal tol)
@ -299,29 +283,22 @@ void SingleSpeciesTP::setState_SV(doublereal s, doublereal v,
void SingleSpeciesTP::initThermo()
{
/*
* Make sure there is one and only one species in this phase.
*/
// Make sure there is one and only one species in this phase.
if (nSpecies() != 1) {
throw CanteraError("initThermo",
"stoichiometric substances may only contain one species.");
}
/*
* Resize temporary arrays.
*/
// Resize temporary arrays.
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;
*/
// Make sure the species mole fraction is equal to 1.0;
double x = 1.0;
ThermoPhase::setMoleFractions(&x);
/*
* Call the base class initThermo object.
*/
// Call the base class initThermo object.
ThermoPhase::initThermo();
}

25
src/thermo/SpeciesThermoFactory.cpp
View file

@ -173,16 +173,14 @@ SpeciesThermoInterpType* newShomateForMineralEQ3(const XML_Node& MinEQ3node)
doublereal e = Entrop_pr_tr * 1.0E3 * 4.184;
doublereal Hcalc = Mu0_tr_pr + 298.15 * e;
/*
* Now calculate the shomate polynomials
*
* Cp first
*
* Shomate: (Joules / gmol / K)
* Cp = As + Bs * t + Cs * t*t + Ds * t*t*t + Es / (t*t)
* where
* t = temperature(Kelvin) / 1000
*/
// Now calculate the shomate polynomials
//
// Cp first
//
// Shomate: (Joules / gmol / K)
// Cp = As + Bs * t + Cs * t*t + Ds * t*t*t + Es / (t*t)
// where
// t = temperature(Kelvin) / 1000
double As = a * 4.184;
double Bs = b * 4.184 * 1000.;
double Cs = 0.0;
@ -382,9 +380,10 @@ static SpeciesThermoInterpType* newAdsorbateThermoFromXML(const XML_Node& f)
SpeciesThermoInterpType* newSpeciesThermoInterpType(const XML_Node& thermo)
{
// Get the children of the thermo XML node. In the next bit of code we take out the comments that
// may have been children of the thermo XML node by doing a selective copy.
// These shouldn't interfere with the algorithm at any point.
// Get the children of the thermo XML node. In the next bit of code we take
// out the comments that may have been children of the thermo XML node by
// doing a selective copy. These shouldn't interfere with the algorithm at
// any point.
const std::vector<XML_Node*>& tpWC = thermo.children();
std::vector<XML_Node*> tp;
for (size_t i = 0; i < tpWC.size(); i++) {

50
src/thermo/StoichSubstance.cpp
View file

@ -20,9 +20,7 @@
namespace Cantera
{
/*
* ---- Constructors -------
*/
// ---- Constructors -------
StoichSubstance::StoichSubstance(const std::string& infile, const std::string& id_)
{
@ -53,18 +51,14 @@ ThermoPhase* StoichSubstance::duplMyselfAsThermoPhase() const
return new StoichSubstance(*this);
}
/*
* ---- Utilities -----
*/
// ---- Utilities -----
int StoichSubstance::eosType() const
{
return cStoichSubstance;
}
/*
* ----- Mechanical Equation of State ------
*/
// ----- Mechanical Equation of State ------
doublereal StoichSubstance::pressure() const
{
@ -86,9 +80,7 @@ doublereal StoichSubstance::thermalExpansionCoeff() const
return 0.0;
}
/*
* ---- Chemical Potentials and Activities ----
*/
// ---- Chemical Potentials and Activities ----
void StoichSubstance::getActivityConcentrations(doublereal* c) const
{
@ -105,9 +97,7 @@ doublereal StoichSubstance::logStandardConc(size_t k) const
return 0.0;
}
/*
* Properties of the Standard State of the Species in the Solution
*/
// Properties of the Standard State of the Species in the Solution
void StoichSubstance::getStandardChemPotentials(doublereal* mu0) const
{
@ -145,9 +135,7 @@ void StoichSubstance::getIntEnergy_RT(doublereal* urt) const
urt[0] = m_h0_RT[0] - m_p0 / molarDensity() / (GasConstant * temperature());
}
/*
* ---- Thermodynamic Values for the Species Reference States ----
*/
// ---- Thermodynamic Values for the Species Reference States ----
void StoichSubstance::getIntEnergy_RT_ref(doublereal* urt) const
{
@ -155,42 +143,32 @@ void StoichSubstance::getIntEnergy_RT_ref(doublereal* urt) const
urt[0] = m_h0_RT[0] - m_p0 / molarDensity() / (GasConstant * temperature());
}
/*
* ---- Initialization and Internal functions
*/
// ---- Initialization and Internal functions
void StoichSubstance::initThermo()
{
/*
* Make sure there is one and only one species in this phase.
*/
// Make sure there is one and only one species in this phase.
if (m_kk != 1) {
throw CanteraError("initThermo",
"stoichiometric substances may only contain one species.");
}
/*
* Store the reference pressure in the variables for the class.
*/
// Store the reference pressure in the variables for the class.
m_p0 = refPressure();
/*
* Resize temporary arrays.
*/
// Resize temporary arrays.
int leng = 1;
m_h0_RT.resize(leng);
m_cp0_R.resize(leng);
m_s0_R.resize(leng);
/*
* Call the base class thermo initializer
*/
// Call the base class thermo initializer
SingleSpeciesTP::initThermo();
}
void StoichSubstance::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
/*
* Find the Thermo XML node
*/
// Find the Thermo XML node
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("StoichSubstance::initThermoXML",
"no thermo XML node");

10
src/thermo/SurfPhase.cpp
View file

@ -273,10 +273,7 @@ void SurfPhase::setCoverages(const doublereal* theta)
for (size_t k = 0; k < m_kk; k++) {
m_work[k] = m_n0*theta[k]/(sum*size(k));
}
/*
* Call the Phase:: class function
* setConcentrations.
*/
// Call the Phase:: class function setConcentrations.
setConcentrations(m_work.data());
}
@ -285,10 +282,7 @@ void SurfPhase::setCoveragesNoNorm(const doublereal* theta)
for (size_t k = 0; k < m_kk; k++) {
m_work[k] = m_n0*theta[k]/size(k);
}
/*
* Call the Phase:: class function
* setConcentrations.
*/
// Call the Phase:: class function setConcentrations.
setConcentrations(m_work.data());
}

75
src/thermo/ThermoFactory.cpp
View file

@ -188,14 +188,17 @@ ThermoPhase* newPhase(const std::string& infile, std::string id)
//! Gather a vector of pointers to XML_Nodes for a phase
/*!
* @param spDataNodeList Output vector of pointer to XML_Nodes which contain the species XML_Nodes for the
* species in the current phase.
* @param spNamesList Output Vector of strings, which contain the names of the species in the phase
* @param spRuleList Output Vector of ints, which contain the value of sprule for each species in the phase
* @param spArray_names Vector of pointers to the XML_Nodes which contains the names of the
* species in the phase
* @param spArray_dbases Input vector of pointers to species data bases.
* We search each data base for the required species names
* @param spDataNodeList Output vector of pointer to XML_Nodes which contain
* the species XML_Nodes for the species in the current phase.
* @param spNamesList Output Vector of strings, which contain the names
* of the species in the phase
* @param spRuleList Output Vector of ints, which contain the value of
* sprule for each species in the phase
* @param spArray_names Vector of pointers to the XML_Nodes which contains
* the names of the species in the phase
* @param spArray_dbases Input vector of pointers to species data bases. We
* search each data base for the required species
* names
* @param sprule Input vector of sprule values
*/
static void formSpeciesXMLNodeList(std::vector<XML_Node*> &spDataNodeList,
@ -220,8 +223,8 @@ static void formSpeciesXMLNodeList(std::vector<XML_Node*> &spDataNodeList,
size_t nsp = spnames.size();
// if 'all' is specified as the one and only species in the
// spArray_names field, then add all species
// defined in the corresponding database to the phase
// spArray_names field, then add all species defined in the
// corresponding database to the phase
if (nsp == 1 && spnames[0] == "all") {
std::vector<XML_Node*> allsp = db->getChildren("species");
nsp = allsp.size();
@ -304,13 +307,10 @@ void importPhase(XML_Node& phase, ThermoPhase* th)
", is not a phase element.");
}
/*
* In this section of code, we get the reference to the
* phase XML tree within the ThermoPhase object. Then,
* we clear it and fill it with the current information that
* we are about to use to construct the object. We will then
* be able to resurrect the information later by calling xml().
*/
// In this section of code, we get the reference to the phase XML tree
// within the ThermoPhase object. Then, we clear it and fill it with the
// current information that we are about to use to construct the object. We
// will then be able to resurrect the information later by calling xml().
th->setXMLdata(phase);
// set the id attribute of the phase to the 'id' attribute in the XML tree.
@ -330,9 +330,9 @@ void importPhase(XML_Node& phase, ThermoPhase* th)
th->setNDim(3); // default
}
// Set equation of state parameters. The parameters are
// specific to each subclass of ThermoPhase, so this is done
// by method setParametersFromXML in each subclass.
// Set equation of state parameters. The parameters are specific to each
// subclass of ThermoPhase, so this is done by method setParametersFromXML
// in each subclass.
const XML_Node& eos = phase.child("thermo");
if (phase.hasChild("thermo")) {
th->setParametersFromXML(eos);
@ -352,20 +352,15 @@ void importPhase(XML_Node& phase, ThermoPhase* th)
}
}
/***************************************************************
* Add the elements.
***************************************************************/
// Add the elements.
if (ssConvention != cSS_CONVENTION_SLAVE) {
installElements(*th, phase);
}
/***************************************************************
* Add the species.
*
* Species definitions may be imported from multiple
* sources. For each one, a speciesArray element must be
* present.
***************************************************************/
// Add the species.
//
// Species definitions may be imported from multiple sources. For each one,
// a speciesArray element must be present.
vector<XML_Node*> sparrays = phase.getChildren("speciesArray");
if (ssConvention != cSS_CONVENTION_SLAVE && sparrays.empty()) {
throw CanteraError("importPhase",
@ -404,11 +399,9 @@ void importPhase(XML_Node& phase, ThermoPhase* th)
}
}
// Get a pointer to the node containing the species
// definitions for the species declared in this
// speciesArray element. This may be in the local file
// containing the phase element, or may be in another
// file.
// Get a pointer to the node containing the species definitions for the
// species declared in this speciesArray element. This may be in the
// local file containing the phase element, or may be in another file.
XML_Node* db = get_XML_Node(speciesArray["datasrc"], &phase.root());
if (db == 0) {
throw CanteraError("importPhase()",
@ -420,10 +413,10 @@ void importPhase(XML_Node& phase, ThermoPhase* th)
dbases.push_back(db);
}
// Now, collect all the species names and all the XML_Node * pointers
// for those species in a single vector. This is where we decide what
// species are to be included in the phase.
// The logic is complicated enough that we put it in a separate routine.
// Now, collect all the species names and all the XML_Node * pointers for
// those species in a single vector. This is where we decide what species
// are to be included in the phase. The logic is complicated enough that we
// put it in a separate routine.
std::vector<XML_Node*> spDataNodeList;
std::vector<std::string> spNamesList;
vector_int spRuleList;
@ -463,8 +456,8 @@ void importPhase(XML_Node& phase, ThermoPhase* th)
// initialization.
th->initThermo();
// Perform any required subclass-specific initialization
// that requires the XML phase object
// Perform any required subclass-specific initialization that requires the
// XML phase object
std::string id = "";
th->initThermoXML(phase, id);
}

103
src/thermo/ThermoPhase.cpp
View file

@ -50,43 +50,31 @@ ThermoPhase::ThermoPhase(const ThermoPhase& right) :
m_chargeNeutralityNecessary(false),
m_ssConvention(cSS_CONVENTION_TEMPERATURE)
{
/*
* Call the assignment operator
*/
// Call the assignment operator
*this = right;
}
ThermoPhase& ThermoPhase::operator=(const ThermoPhase& right)
{
/*
* Check for self assignment.
*/
// Check for self assignment.
if (this == &right) {
return *this;
}
/*
* We need to destruct first
*/
// We need to destruct first
for (size_t k = 0; k < m_speciesData.size(); k++) {
delete m_speciesData[k];
}
delete m_spthermo;
/*
* Call the base class assignment operator
*/
// Call the base class assignment operator
Phase::operator=(right);
/*
* Pointer to the species thermodynamic property manager
* We own this, so we need to do a deep copy
*/
// Pointer to the species thermodynamic property manager
// We own this, so we need to do a deep copy
m_spthermo = (right.m_spthermo)->duplMyselfAsSpeciesThermo();
/*
* Do a deep copy of species Data, because we own this
*/
// Do a deep copy of species Data, because we own this
m_speciesData.resize(m_kk);
for (size_t k = 0; k < m_kk; k++) {
m_speciesData[k] = new XML_Node(*(right.m_speciesData[k]));
@ -291,9 +279,8 @@ void ThermoPhase::setState_HPorUV(doublereal Htarget, doublereal p,
double Tbot = Tnew;
bool ignoreBounds = false;
// Unstable phases are those for which
// cp < 0.0. These are possible for cases where
// we have passed the spinodal curve.
// Unstable phases are those for which cp < 0.0. These are possible for
// cases where we have passed the spinodal curve.
bool unstablePhase = false;
// Counter indicating the last temperature point where the
// phase was unstable
@ -315,9 +302,8 @@ void ThermoPhase::setState_HPorUV(doublereal Htarget, doublereal p,
// Calculate the new T
Tnew = Told + dt;
// Limit the step size so that we are convergent
// This is the step that makes it different from a
// Newton's algorithm
// Limit the step size so that we are convergent This is the step that
// makes it different from a Newton's algorithm
if ((dt > 0.0 && unstablePhase) || (dt <= 0.0 && !unstablePhase)) {
if (Hbot < Htarget && Tnew < (0.75 * Tbot + 0.25 * Told)) {
dt = 0.75 * (Tbot - Told);
@ -404,10 +390,9 @@ void ThermoPhase::setState_HPorUV(doublereal Htarget, doublereal p,
}
}
// We are here when there hasn't been convergence
/*
* Formulate a detailed error message, since questions seem to
* arise often about the lack of convergence.
*/
// Formulate a detailed error message, since questions seem to arise often
// about the lack of convergence.
string ErrString = "No convergence in 500 iterations\n";
if (doUV) {
ErrString += fmt::format(
@ -495,9 +480,8 @@ void ThermoPhase::setState_SPorSV(doublereal Starget, doublereal p,
double Tbot = Tnew;
bool ignoreBounds = false;
// Unstable phases are those for which
// Cp < 0.0. These are possible for cases where
// we have passed the spinodal curve.
// Unstable phases are those for which Cp < 0.0. These are possible for
// cases where we have passed the spinodal curve.
bool unstablePhase = false;
double Tunstable = -1.0;
bool unstablePhaseNew = false;
@ -592,10 +576,9 @@ void ThermoPhase::setState_SPorSV(doublereal Starget, doublereal p,
}
}
// We are here when there hasn't been convergence
/*
* Formulate a detailed error message, since questions seem to
* arise often about the lack of convergence.
*/
// Formulate a detailed error message, since questions seem to arise often
// about the lack of convergence.
string ErrString = "No convergence in 500 iterations\n";
if (doSV) {
ErrString += fmt::format(
@ -848,9 +831,7 @@ void ThermoPhase::getdlnActCoeffdlnN_numderiv(const size_t ld, doublereal* const
double deltaMoles_j = 0.0;
double pres = pressure();
/*
* Evaluate the current base activity coefficients if necessary
*/
// Evaluate the current base activity coefficients if necessary
vector_fp ActCoeff_Base(m_kk);
getActivityCoefficients(ActCoeff_Base.data());
vector_fp Xmol_Base(m_kk);
@ -862,56 +843,40 @@ void ThermoPhase::getdlnActCoeffdlnN_numderiv(const size_t ld, doublereal* const
double v_totalMoles = 1.0;
double TMoles_base = v_totalMoles;
/*
* Loop over the columns species to be deltad
*/
// Loop over the columns species to be deltad
for (size_t j = 0; j < m_kk; j++) {
/*
* Calculate a value for the delta moles of species j
* -> NOte Xmol_[] and Tmoles are always positive or zero
* quantities.
* -> experience has shown that you always need to make the deltas greater than needed to
* change the other mole fractions in order to capture some effects.
*/
// Calculate a value for the delta moles of species j
// -> Note Xmol_[] and Tmoles are always positive or zero quantities.
// -> experience has shown that you always need to make the deltas
// greater than needed to change the other mole fractions in order
// to capture some effects.
double moles_j_base = v_totalMoles * Xmol_Base[j];
deltaMoles_j = 1.0E-7 * moles_j_base + v_totalMoles * 1.0E-13 + 1.0E-150;
/*
* Now, update the total moles in the phase and all of the
* mole fractions based on this.
*/
// Now, update the total moles in the phase and all of the mole
// fractions based on this.
v_totalMoles = TMoles_base + deltaMoles_j;
for (size_t k = 0; k < m_kk; k++) {
Xmol[k] = Xmol_Base[k] * TMoles_base / v_totalMoles;
}
Xmol[j] = (moles_j_base + deltaMoles_j) / v_totalMoles;
/*
* Go get new values for the activity coefficients.
* -> Note this calls setState_PX();
*/
// Go get new values for the activity coefficients.
// -> Note this calls setState_PX();
setState_PX(pres, Xmol.data());
getActivityCoefficients(ActCoeff.data());
/*
* Calculate the column of the matrix
*/
// Calculate the column of the matrix
double* const lnActCoeffCol = dlnActCoeffdlnN + ld * j;
for (size_t k = 0; k < m_kk; k++) {
lnActCoeffCol[k] = (2*moles_j_base + deltaMoles_j) *(ActCoeff[k] - ActCoeff_Base[k]) /
((ActCoeff[k] + ActCoeff_Base[k]) * deltaMoles_j);
}
/*
* Revert to the base case Xmol_, v_totalMoles
*/
// Revert to the base case Xmol_, v_totalMoles
v_totalMoles = TMoles_base;
Xmol = Xmol_Base;
}
/*
* Go get base values for the activity coefficients.
* -> Note this calls setState_TPX() again;
* -> Just wanted to make sure that cantera is in sync
* with VolPhase after this call.
*/
setState_PX(pres, Xmol_Base.data());
}

29
src/thermo/VPSSMgr.cpp
View file

@ -63,19 +63,19 @@ VPSSMgr& VPSSMgr::operator=(const VPSSMgr& right)
return *this;
}
m_kk = right.m_kk;
/*
* What we are doing here is to make a shallow copy of the VPStandardStateTP
* pointer in the "new" VPSSMgr object using the value from the "old"
* VPSSMgr object. This is not appropriate if we are making a copy of a ThermoPhase
* object and the VPSSMgr objects are owned by the ThermoPhase object.
*
* The new object will want to have a different value of m_vptp_ptr than the
* value this is being copied here. It will want to refer to the copy of the
* VPStandardStateTP object being made that will own the new VPSSMgr object.
* However, the assignment object is not the place to carry out this fixup.
*
* We will have to "fix" up the shallow copies later.
*/
// What we are doing here is to make a shallow copy of the VPStandardStateTP
// pointer in the "new" VPSSMgr object using the value from the "old"
// VPSSMgr object. This is not appropriate if we are making a copy of a
// ThermoPhase object and the VPSSMgr objects are owned by the ThermoPhase
// object.
//
// The new object will want to have a different value of m_vptp_ptr than the
// value this is being copied here. It will want to refer to the copy of the
// VPStandardStateTP object being made that will own the new VPSSMgr object.
// However, the assignment object is not the place to carry out this fixup.
//
// We will have to "fix" up the shallow copies later.
m_vptp_ptr = right.m_vptp_ptr;
m_spthermo = right.m_spthermo;
m_tlast = -1.0;
@ -212,6 +212,7 @@ const vector_fp& VPSSMgr::getStandardVolumes() const
}
/*****************************************************************/
void VPSSMgr::getEnthalpy_RT_ref(doublereal* hrt) const
{
if (m_useTmpRefStateStorage) {
@ -387,6 +388,7 @@ PDSS* VPSSMgr::createInstallPDSS(size_t k, const XML_Node& s,
}
/*****************************************************************/
doublereal VPSSMgr::minTemp(size_t k) const
{
if (k != npos) {
@ -416,7 +418,6 @@ PDSS_enumType VPSSMgr::reportPDSSType(int index) const
throw NotImplementedError("VPSSMgr::reportPDSSType()");
}
VPSSMgr_enumType VPSSMgr::reportVPSSMgrType() const
{
throw NotImplementedError("VPSSMgr::reportVPSSType()");

63
src/thermo/VPSSMgrFactory.h
View file

@ -41,26 +41,20 @@ public:
//! Factory to build instances of classes that manage the
//! standard-state thermodynamic properties of a set of species.
/*!
* This class is responsible for making the decision concerning
* which derivative of VPSSMgr object to use.
* The VPSSMgr object is used to calculate
* thermodynamic functions for the standard state.
* It queries the database of species to understand what
* the requirements are for the submodels for all of the
* species in the phase. Then, it picks the derived VPSSMgr
* object to use and passes it back to the calling routine.
* It doesn't load any data into the derived
* VPSSMgr object.
* This class is responsible for making the decision concerning which
* derivative of VPSSMgr object to use. The VPSSMgr object is used to calculate
* thermodynamic functions for the standard state. It queries the database of
* species to understand what the requirements are for the submodels for all of
* the species in the phase. Then, it picks the derived VPSSMgr object to use
* and passes it back to the calling routine. It doesn't load any data into the
* derived VPSSMgr object.
*
* Making the choice of VPSSMgr types is the only
* thing this class does.
* Making the choice of VPSSMgr types is the only thing this class does.
*
* This class is implemented as a singleton -- one in which
* only one instance is needed. The recommended way to access
* the factory is to call this static method, which
* instantiates the class if it is the first call, but
* otherwise simply returns the pointer to the existing
* instance.
* This class is implemented as a singleton -- one in which only one instance is
* needed. The recommended way to access the factory is to call this static
* method, which instantiates the class if it is the first call, but otherwise
* simply returns the pointer to the existing instance.
*
* @ingroup mgrpdssthermocalc
*/
@ -69,12 +63,10 @@ class VPSSMgrFactory : public FactoryBase
public:
//! Static method to return an instance of this class
/*!
* This class is implemented as a singleton -- one in which
* only one instance is needed. The recommended way to access
* the factory is to call this static method, which
* instantiates the class if it is the first call, but
* otherwise simply returns the pointer to the existing
* instance.
* This class is implemented as a singleton -- one in which only one
* instance is needed. The recommended way to access the factory is to call
* this static method, which instantiates the class if it is the first call,
* but otherwise simply returns the pointer to the existing instance.
*/
static VPSSMgrFactory* factory() {
std::unique_lock<std::mutex> lock(vpss_species_thermo_mutex);
@ -86,17 +78,16 @@ public:
//! Delete static instance of this class
/*!
* If it is necessary to explicitly delete the factory before
* the process terminates (for example, when checking for
* memory leaks) then this method can be called to delete it.
* If it is necessary to explicitly delete the factory before the process
* terminates (for example, when checking for memory leaks) then this method
* can be called to delete it.
*/
void deleteFactory();
//! String conversion to an enumType
/*!
* This routine is a string conversion. The string is obtained from the
* standardState model attribute and converted to a VPSSMgr_enumType
* type.
* standardState model attribute and converted to a VPSSMgr_enumType type.
*
* @param ssModel String representing the VPSSMGr object
*/
@ -113,18 +104,18 @@ public:
//! Create a new species property manager for a group of species
/*!
* This routine will look through species nodes. It will discover what
* each species needs for its species property managers. Then,
* it will malloc and return the proper species property manager to use.
* This routine will look through species nodes. It will discover what each
* species needs for its species property managers. Then, it will malloc and
* return the proper species property manager to use.
*
* @param vp_ptr Variable pressure standard state ThermoPhase object
* that will be the owner.
* @param phaseNode_ptr Pointer to the ThermoPhase phase XML Node
* @param spDataNodeList Vector of XML_Nodes, each of which is a species XML Node.
* There are m_kk of these.
* @param spDataNodeList Vector of XML_Nodes, each of which is a species XML
* Node. There are m_kk of these.
*
* @return Returns a pointer to a newly malloced species property
* manager object.
* @return Returns a pointer to a newly malloced species
* property manager object.
*/
virtual VPSSMgr* newVPSSMgr(VPStandardStateTP* vp_ptr,
XML_Node* phaseNode_ptr,

19
src/thermo/VPSSMgr_General.cpp
View file

@ -54,12 +54,12 @@ VPSSMgr_General& VPSSMgr_General::operator=(const VPSSMgr_General& b)
return *this;
}
VPSSMgr::operator=(b);
/*
* Must fill in the shallow pointers. These must have already been transfered
* and stored in the owning VPStandardStateTP class. Note we are aware that at this point
* m_vptr_ptr may refer back to the wrong ThermoPhase object. However, the shallow copy
* performed here is consistent with the assignment operator's general functionality.
*/
// Must fill in the shallow pointers. These must have already been
// transfered and stored in the owning VPStandardStateTP class. Note we are
// aware that at this point m_vptr_ptr may refer back to the wrong
// ThermoPhase object. However, the shallow copy performed here is
// consistent with the assignment operator's general functionality.
m_PDSS_ptrs.resize(m_kk);
for (size_t k = 0; k < m_kk; k++) {
m_PDSS_ptrs[k] = m_vptp_ptr->providePDSS(k);
@ -75,10 +75,9 @@ VPSSMgr* VPSSMgr_General::duplMyselfAsVPSSMgr() const
void VPSSMgr_General::initAllPtrs(VPStandardStateTP* vp_ptr, SpeciesThermo* sp_ptr)
{
VPSSMgr::initAllPtrs(vp_ptr, sp_ptr);
/*
* Must fill in the shallow pointers. These must have already been transfered
* and stored in the owning VPStandardStateTP class.
*/
// Must fill in the shallow pointers. These must have already been
// transfered and stored in the owning VPStandardStateTP class.
m_PDSS_ptrs.resize(m_kk);
for (size_t k = 0; k < m_kk; k++) {
m_PDSS_ptrs[k] = m_vptp_ptr->providePDSS(k);

88
src/thermo/VPStandardStateTP.cpp
View file

@ -19,9 +19,6 @@ using namespace std;
namespace Cantera
{
/*
* Default constructor
*/
VPStandardStateTP::VPStandardStateTP() :
m_Pcurrent(OneAtm),
m_Tlast_ss(-1.0),
@ -44,22 +41,17 @@ VPStandardStateTP::VPStandardStateTP(const VPStandardStateTP& b) :
VPStandardStateTP& VPStandardStateTP::operator=(const VPStandardStateTP& b)
{
if (&b != this) {
/*
* Mostly, this is a passthrough to the underlying
* assignment operator for the ThermoPhase parent object.
*/
// Mostly, this is a passthrough to the underlying assignment operator
// for the ThermoPhase parent object.
ThermoPhase::operator=(b);
/*
* However, we have to handle data that we own.
*/
// However, we have to handle data that we own.
m_Pcurrent = b.m_Pcurrent;
m_Tlast_ss = b.m_Tlast_ss;
m_Plast_ss = b.m_Plast_ss;
m_P0 = b.m_P0;
/*
* Duplicate the pdss objects
*/
// Duplicate the pdss objects
if (m_PDSS_storage.size() > 0) {
for (int k = 0; k < (int) m_PDSS_storage.size(); k++) {
delete m_PDSS_storage[k];
@ -70,32 +62,26 @@ VPStandardStateTP& VPStandardStateTP::operator=(const VPStandardStateTP& b)
m_PDSS_storage[k] = b.m_PDSS_storage[k]->duplMyselfAsPDSS();
}
/*
* Duplicate the VPSS Manager object that conducts the calculations
*/
// Duplicate the VPSS Manager object that conducts the calculations
delete m_VPSS_ptr;
m_VPSS_ptr = (b.m_VPSS_ptr)->duplMyselfAsVPSSMgr();
/*
* The VPSSMgr object contains shallow pointers. Whenever you have shallow
* pointers, they have to be fixed up to point to the correct objects referring
* back to this ThermoPhase's properties.
*/
// The VPSSMgr object contains shallow pointers. Whenever you have
// shallow pointers, they have to be fixed up to point to the correct
// objects referring back to this ThermoPhase's properties.
m_VPSS_ptr->initAllPtrs(this, m_spthermo);
/*
* The PDSS objects contains shallow pointers. Whenever you have shallow
* pointers, they have to be fixed up to point to the correct objects referring
* back to this ThermoPhase's properties. This function also sets m_VPSS_ptr
* so it occurs after m_VPSS_ptr is set.
*/
// The PDSS objects contains shallow pointers. Whenever you have shallow
// pointers, they have to be fixed up to point to the correct objects
// referring back to this ThermoPhase's properties. This function also
// sets m_VPSS_ptr so it occurs after m_VPSS_ptr is set.
for (size_t k = 0; k < m_kk; k++) {
m_PDSS_storage[k]->initAllPtrs(this, m_VPSS_ptr, m_spthermo);
}
/*
* Ok, the VPSSMgr object is ready for business.
* We need to resync the temperature and the pressure of the new standard states
* with what is stored in this object.
*/
// Ok, the VPSSMgr object is ready for business. We need to resync the
// temperature and the pressure of the new standard states with what is
// stored in this object.
m_VPSS_ptr->setState_TP(m_Tlast_ss, m_Plast_ss);
}
return *this;
@ -127,9 +113,8 @@ void VPStandardStateTP::getChemPotentials_RT(doublereal* muRT) const
}
}
/*
* ----- Thermodynamic Values for the Species Standard States States ----
*/
// ----- Thermodynamic Values for the Species Standard States States ----
void VPStandardStateTP::getStandardChemPotentials(doublereal* g) const
{
getGibbs_RT(g);
@ -194,9 +179,7 @@ const vector_fp& VPStandardStateTP::getStandardVolumes() const
return m_VPSS_ptr->getStandardVolumes();
}
/*
* ----- Thermodynamic Values for the Species Reference States ----
*/
// ----- Thermodynamic Values for the Species Reference States ----
void VPStandardStateTP::getEnthalpy_RT_ref(doublereal* hrt) const
{
@ -284,26 +267,21 @@ void VPStandardStateTP::calcDensity()
void VPStandardStateTP::setState_TP(doublereal t, doublereal pres)
{
/*
* A pretty tricky algorithm is needed here, due to problems involving
* standard states of real fluids. For those cases you need
* to combine the T and P specification for the standard state, or else
* you may venture into the forbidden zone, especially when nearing the
* triple point.
* Therefore, we need to do the standard state thermo calc with the
* (t, pres) combo.
*/
// A pretty tricky algorithm is needed here, due to problems involving
// standard states of real fluids. For those cases you need to combine the T
// and P specification for the standard state, or else you may venture into
// the forbidden zone, especially when nearing the triple point. Therefore,
// we need to do the standard state thermo calc with the (t, pres) combo.
Phase::setTemperature(t);
m_Pcurrent = pres;
updateStandardStateThermo();
/*
* Now, we still need to do the calculations for general ThermoPhase objects.
* So, we switch back to a virtual function call, setTemperature, and
* setPressure to recalculate stuff for child ThermoPhase objects of
* the VPStandardStateTP object. At this point,
* we haven't touched m_tlast or m_plast, so some calculations may still
* need to be done at the ThermoPhase object level.
*/
// Now, we still need to do the calculations for general ThermoPhase
// objects. So, we switch back to a virtual function call, setTemperature,
// and setPressure to recalculate stuff for child ThermoPhase objects of the
// VPStandardStateTP object. At this point, we haven't touched m_tlast or
// m_plast, so some calculations may still need to be done at the
// ThermoPhase object level.
calcDensity();
}

58
src/thermo/WaterProps.cpp
View file

@ -98,11 +98,8 @@ doublereal WaterProps::density_T(doublereal T, doublereal P, int ifunc)
doublereal tmp3 = Tc + U3;
doublereal rho = 1000. * (1.0 - tmp1*t4t4/(U2 * tmp3));
/*
* Impose an ideal gas lower bound on rho. We need this
* to ensure positivity of rho, even though it is
* grossly unrepresentative.
*/
// Impose an ideal gas lower bound on rho. We need this to ensure positivity
// of rho, even though it is grossly unrepresentative.
doublereal rhomin = P / (GasConstant * T);
if (rho < rhomin) {
rho = rhomin;
@ -205,19 +202,15 @@ doublereal WaterProps::ADebye(doublereal T, doublereal P_input, int ifunc)
doublereal tmp3 = tmp2 * sqrt(tmp2);
doublereal A_Debye = tmp * tmp3 / (8.0 * Pi);
/*
* dAdT = - 3/2 Ad/T + 1/2 Ad/dw d(dw)/dT - 3/2 Ad/eps d(eps)/dT
* dAdT = - 3/2 Ad/T - 1/2 Ad/Vw d(Vw)/dT - 3/2 Ad/eps d(eps)/dT
*/
// dAdT = - 3/2 Ad/T + 1/2 Ad/dw d(dw)/dT - 3/2 Ad/eps d(eps)/dT
// dAdT = - 3/2 Ad/T - 1/2 Ad/Vw d(Vw)/dT - 3/2 Ad/eps d(eps)/dT
if (ifunc == 1 || ifunc == 2) {
doublereal dAdT = - 1.5 * A_Debye / T;
doublereal depsRelWaterdT = relEpsilon(T, P, 1);
dAdT -= A_Debye * (1.5 * depsRelWaterdT / epsRelWater);
/*
* calculate d(lnV)/dT _constantP, i.e., the cte
*/
// calculate d(lnV)/dT _constantP, i.e., the cte
doublereal cte = coeffThermalExp_IAPWS(T, P);
doublereal contrib2 = - A_Debye * (0.5 * cte);
dAdT += contrib2;
@ -227,11 +220,9 @@ doublereal WaterProps::ADebye(doublereal T, doublereal P_input, int ifunc)
}
if (ifunc == 2) {
/*
* Get the second derivative of the dielectric constant wrt T
* -> we will take each of the terms in dAdT and differentiate
* it again.
*/
// Get the second derivative of the dielectric constant wrt T
// -> we will take each of the terms in dAdT and differentiate
// it again.
doublereal d2AdT2 = 1.5 / T * (A_Debye/T - dAdT);
doublereal d2epsRelWaterdT2 = relEpsilon(T, P, 2);
d2AdT2 += 1.5 * (- dAdT * depsRelWaterdT / epsRelWater
@ -246,16 +237,15 @@ doublereal WaterProps::ADebye(doublereal T, doublereal P_input, int ifunc)
return d2AdT2;
}
}
/*
* A_Debye = (1/(8 Pi)) sqrt(2 Na dw / 1000)
* (e e/(epsilon R T))^3/2
*
* dAdP = + 1/2 Ad/dw d(dw)/dP - 3/2 Ad/eps d(eps)/dP
* dAdP = - 1/2 Ad/Vw d(Vw)/dP - 3/2 Ad/eps d(eps)/dP
* dAdP = + 1/2 Ad * kappa - 3/2 Ad/eps d(eps)/dP
*
* where kappa = - 1/Vw d(Vw)/dP_T (isothermal compressibility)
*/
// A_Debye = (1/(8 Pi)) sqrt(2 Na dw / 1000)
// (e e/(epsilon R T))^3/2
//
// dAdP = + 1/2 Ad/dw d(dw)/dP - 3/2 Ad/eps d(eps)/dP
// dAdP = - 1/2 Ad/Vw d(Vw)/dP - 3/2 Ad/eps d(eps)/dP
// dAdP = + 1/2 Ad * kappa - 3/2 Ad/eps d(eps)/dP
//
// where kappa = - 1/Vw d(Vw)/dP_T (isothermal compressibility)
if (ifunc == 3) {
doublereal dAdP = 0.0;
doublereal depsRelWaterdP = relEpsilon(T, P, 3);
@ -437,14 +427,12 @@ doublereal WaterProps::thermalConductivityWater() const
doublereal rho4 = rho2 * rho2;
doublereal temp2 = (tbar - 1.0) * (tbar - 1.0);
/*
* beta = M / (rho * Rgas) (d (pressure) / dT) at constant rho
*
* Note for ideal gases this is equal to one.
*
* beta = delta (phi0_d() + phiR_d())
* - tau delta (phi0_dt() + phiR_dt())
*/
// beta = M / (rho * Rgas) (d (pressure) / dT) at constant rho
//
// Note for ideal gases this is equal to one.
//
// beta = delta (phi0_d() + phiR_d())
// - tau delta (phi0_dt() + phiR_dt())
doublereal beta = m_waterIAPWS->coeffPresExp();
doublereal dpdT_const_rho = beta * GasConstant * dens / 18.015268;
dpdT_const_rho *= Tstar / presstar;

54
src/thermo/WaterPropsIAPWS.cpp
View file

@ -15,9 +15,8 @@
namespace Cantera
{
/*
* Critical Point values of water in mks units
*/
// Critical Point values of water in mks units
//! Critical Temperature value (kelvin)
const doublereal T_c = 647.096;
//! Critical Pressure (Pascals)
@ -68,9 +67,8 @@ void WaterPropsIAPWS::calcDim(doublereal temperature, doublereal rho)
{
tau = T_c / temperature;
delta = rho / Rho_c;
/*
* Determine the internal state
*/
// Determine the internal state
if (temperature > T_c) {
iState = WATER_SUPERCRIT;
} else {
@ -110,10 +108,8 @@ doublereal WaterPropsIAPWS::density(doublereal temperature, doublereal pressure,
if (phase == WATER_GAS || phase == WATER_SUPERCRIT) {
rhoguess = pressure * M_water / (Rgas * temperature);
} else if (phase == WATER_LIQUID) {
/*
* Provide a guess about the liquid density that is
* relatively high -> convergence from above seems robust.
*/
// Provide a guess about the liquid density that is
// relatively high -> convergence from above seems robust.
rhoguess = 1000.;
} else if (phase == WATER_UNSTABLELIQUID || phase == WATER_UNSTABLEGAS) {
throw CanteraError("WaterPropsIAPWS::density",
@ -124,10 +120,8 @@ doublereal WaterPropsIAPWS::density(doublereal temperature, doublereal pressure,
}
}
} else {
/*
* Assume the Gas phase initial guess, if nothing is
* specified to the routine
*/
// Assume the Gas phase initial guess, if nothing is specified to
// the routine
rhoguess = pressure * M_water / (Rgas * temperature);
}
}
@ -139,14 +133,11 @@ doublereal WaterPropsIAPWS::density(doublereal temperature, doublereal pressure,
if (delta_retn >0.0) {
delta = delta_retn;
/*
* Dimensionalize the density before returning
*/
// Dimensionalize the density before returning
density_retn = delta_retn * Rho_c;
/*
* Set the internal state -> this may be
* a duplication. However, let's just be sure.
*/
// Set the internal state -> this may be a duplication. However, let's
// just be sure.
setState_TR(temperature, density_retn);
} else {
density_retn = -1.0;
@ -168,10 +159,8 @@ doublereal WaterPropsIAPWS::density_const(doublereal pressure,
if (phase == WATER_GAS || phase == WATER_SUPERCRIT) {
rhoguess = pressure * M_water / (Rgas * temperature);
} else if (phase == WATER_LIQUID) {
/*
* Provide a guess about the liquid density that is
* relatively high -> convergence from above seems robust.
*/
// Provide a guess about the liquid density that is
// relatively high -> convergence from above seems robust.
rhoguess = 1000.;
} else if (phase == WATER_UNSTABLELIQUID || phase == WATER_UNSTABLEGAS) {
throw CanteraError("WaterPropsIAPWS::density",
@ -182,10 +171,8 @@ doublereal WaterPropsIAPWS::density_const(doublereal pressure,
}
}
} else {
/*
* Assume the Gas phase initial guess, if nothing is
* specified to the routine
*/
// Assume the Gas phase initial guess, if nothing is specified to
// the routine
rhoguess = pressure * M_water / (Rgas * temperature);
}
}
@ -200,9 +187,7 @@ doublereal WaterPropsIAPWS::density_const(doublereal pressure,
if (delta_retn > 0.0) {
delta = delta_retn;
/*
* Dimensionalize the density before returning
*/
// Dimensionalize the density before returning
density_retn = delta_retn * Rho_c;
} else {
@ -252,9 +237,8 @@ doublereal WaterPropsIAPWS::psat_est(doublereal temperature) const
doublereal q = b / v;
ps = 22.093*exp(q);
}
/*
* Original correlation was in cgs. Convert to mks
*/
// Original correlation was in cgs. Convert to mks
ps *= 1.0E6;
return ps;
}

156
src/thermo/WaterPropsIAPWSphi.cpp
View file

@ -458,9 +458,7 @@ doublereal WaterPropsIAPWSphi::phiR() const
doublereal delta = DELTAsave;
int i, j;
/*
* Write out the first seven polynomials in the expression
*/
// Write out the first seven polynomials in the expression
doublereal T375 = pow(tau, 0.375);
doublereal val = (ni[1] * delta / TAUsqrt +
ni[2] * delta * TAUsqrt * T375 +
@ -469,16 +467,12 @@ doublereal WaterPropsIAPWSphi::phiR() const
ni[5] * DELTAp[2] * T375 * T375 +
ni[6] * DELTAp[3] * T375 +
ni[7] * DELTAp[4] * tau);
/*
* Next, do polynomial contributions 8 to 51
*/
// Next, do polynomial contributions 8 to 51
for (i = 8; i <= 51; i++) {
val += (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] * exp(-DELTAp[ciR[i]]));
}
/*
* Next do contributions 52 to 54
*/
// Next do contributions 52 to 54
for (j = 0; j < 3; j++) {
i = 52 + j;
doublereal dtmp = delta - epsi[j];
@ -487,9 +481,7 @@ doublereal WaterPropsIAPWSphi::phiR() const
exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
}
/*
* Next do contributions 55 and 56
*/
// Next do contributions 55 and 56
for (j = 0; j < 2; j++) {
i = 55 + j;
doublereal deltam1 = delta - 1.0;
@ -520,9 +512,7 @@ doublereal WaterPropsIAPWSphi::phiR_d() const
doublereal delta = DELTAsave;
int i, j;
/*
* Write out the first seven polynomials in the expression
*/
// Write out the first seven polynomials in the expression
doublereal T375 = pow(tau, 0.375);
doublereal val = (ni[1] / TAUsqrt +
ni[2] * TAUsqrt * T375 +
@ -531,17 +521,13 @@ doublereal WaterPropsIAPWSphi::phiR_d() const
ni[5] * 2.0 * delta * T375 * T375 +
ni[6] * 3.0 * DELTAp[2] * T375 +
ni[7] * 4.0 * DELTAp[3] * tau);
/*
* Next, do polynomial contributions 8 to 51
*/
// Next, do polynomial contributions 8 to 51
for (i = 8; i <= 51; i++) {
val += ((ni[i] * exp(-DELTAp[ciR[i]]) * DELTAp[diR[i] - 1] *
TAUp[tiR[i]]) * (diR[i] - ciR[i]* DELTAp[ciR[i]]));
}
/*
* Next do contributions 52 to 54
*/
// Next do contributions 52 to 54
for (j = 0; j < 3; j++) {
i = 52 + j;
doublereal dtmp = delta - epsi[j];
@ -551,9 +537,7 @@ doublereal WaterPropsIAPWSphi::phiR_d() const
val += tmp * (diR[i]/delta - 2.0 * alphai[j] * dtmp);
}
/*
* Next do contributions 55 and 56
*/
// Next do contributions 55 and 56
for (j = 0; j < 2; j++) {
i = 55 + j;
doublereal deltam1 = delta - 1.0;
@ -608,17 +592,13 @@ doublereal WaterPropsIAPWSphi::phiR_dd() const
int i, j;
doublereal atmp;
/*
* Write out the first seven polynomials in the expression
*/
// Write out the first seven polynomials in the expression
doublereal T375 = pow(tau, 0.375);
doublereal val = (ni[4] * 2.0 * TAUsqrt +
ni[5] * 2.0 * T375 * T375 +
ni[6] * 6.0 * delta * T375 +
ni[7] * 12.0 * DELTAp[2] * tau);
/*
* Next, do polynomial contributions 8 to 51
*/
// Next, do polynomial contributions 8 to 51
for (i = 8; i <= 51; i++) {
doublereal dtmp = DELTAp[ciR[i]];
doublereal tmp = ni[i] * exp(-dtmp) * TAUp[tiR[i]];
@ -632,9 +612,7 @@ doublereal WaterPropsIAPWSphi::phiR_dd() const
val += tmp;
}
/*
* Next do contributions 52 to 54
*/
// Next do contributions 52 to 54
for (j = 0; j < 3; j++) {
i = 52 + j;
doublereal dtmp = delta - epsi[j];
@ -652,9 +630,7 @@ doublereal WaterPropsIAPWSphi::phiR_dd() const
diR[i] * (diR[i] - 1.0) * deltmpM2);
}
/*
* Next do contributions 55 and 56
*/
// Next do contributions 55 and 56
for (j = 0; j < 2; j++) {
i = 55 + j;
doublereal deltam1 = delta - 1.0;
@ -741,9 +717,7 @@ doublereal WaterPropsIAPWSphi::phiR_t() const
int i, j;
doublereal atmp, tmp;
/*
* Write out the first seven polynomials in the expression
*/
// Write out the first seven polynomials in the expression
doublereal T375 = pow(tau, 0.375);
doublereal val = ((-0.5) *ni[1] * delta / TAUsqrt / tau +
ni[2] * delta * 0.875 / TAUsqrt * T375 +
@ -752,17 +726,13 @@ doublereal WaterPropsIAPWSphi::phiR_t() const
ni[5] * DELTAp[2] * 0.75 * T375 * T375 / tau +
ni[6] * DELTAp[3] * 0.375 * T375 / tau +
ni[7] * DELTAp[4]);
/*
* Next, do polynomial contributions 8 to 51
*/
// Next, do polynomial contributions 8 to 51
for (i = 8; i <= 51; i++) {
tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]-1] * exp(-DELTAp[ciR[i]]));
val += tiR[i] * tmp;
}
/*
* Next do contributions 52 to 54
*/
// Next do contributions 52 to 54
for (j = 0; j < 3; j++) {
i = 52 + j;
doublereal dtmp = delta - epsi[j];
@ -772,9 +742,7 @@ doublereal WaterPropsIAPWSphi::phiR_t() const
val += tmp *(tiR[i]/tau - 2.0 * betai[j]*ttmp);
}
/*
* Next do contributions 55 and 56
*/
// Next do contributions 55 and 56
for (j = 0; j < 2; j++) {
i = 55 + j;
doublereal deltam1 = delta - 1.0;
@ -820,18 +788,14 @@ doublereal WaterPropsIAPWSphi::phiR_tt() const
int i, j;
doublereal atmp, tmp;
/*
* Write out the first seven polynomials in the expression
*/
// Write out the first seven polynomials in the expression
doublereal T375 = pow(tau, 0.375);
doublereal val = ((-0.5) * (-1.5) * ni[1] * delta / (TAUsqrt * tau * tau) +
ni[2] * delta * 0.875 * (-0.125) * T375 / (TAUsqrt * tau) +
ni[4] * DELTAp[2] * 0.5 * (-0.5)/ (TAUsqrt * tau) +
ni[5] * DELTAp[2] * 0.75 *(-0.25) * T375 * T375 / (tau * tau) +
ni[6] * DELTAp[3] * 0.375 *(-0.625) * T375 / (tau * tau));
/*
* Next, do polynomial contributions 8 to 51
*/
// Next, do polynomial contributions 8 to 51
for (i = 8; i <= 51; i++) {
if (tiR[i] > 1) {
tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]-2] * exp(-DELTAp[ciR[i]]));
@ -839,9 +803,7 @@ doublereal WaterPropsIAPWSphi::phiR_tt() const
}
}
/*
* Next do contributions 52 to 54
*/
// Next do contributions 52 to 54
for (j = 0; j < 3; j++) {
i = 52 + j;
doublereal dtmp = delta - epsi[j];
@ -852,9 +814,7 @@ doublereal WaterPropsIAPWSphi::phiR_tt() const
val += tmp *(atmp * atmp - tiR[i]/(tau*tau) - 2.0*betai[j]);
}
/*
* Next do contributions 55 and 56
*/
// Next do contributions 55 and 56
for (j = 0; j < 2; j++) {
i = 55 + j;
doublereal deltam1 = delta - 1.0;
@ -900,9 +860,7 @@ doublereal WaterPropsIAPWSphi::phiR_dt() const
doublereal delta = DELTAsave;
int i, j;
doublereal tmp;
/*
* Write out the first seven polynomials in the expression
*/
// Write out the first seven polynomials in the expression
doublereal T375 = pow(tau, 0.375);
doublereal val = (ni[1] * (-0.5) / (TAUsqrt * tau) +
ni[2] * (0.875) * T375 / TAUsqrt +
@ -911,18 +869,14 @@ doublereal WaterPropsIAPWSphi::phiR_dt() const
ni[5] * 2.0 * delta * (0.75) * T375 * T375 / tau +
ni[6] * 3.0 * DELTAp[2] * 0.375 * T375 / tau +
ni[7] * 4.0 * DELTAp[3]);
/*
* Next, do polynomial contributions 8 to 51
*/
// Next, do polynomial contributions 8 to 51
for (i = 8; i <= 51; i++) {
tmp = (ni[i] * tiR[i] * exp(-DELTAp[ciR[i]]) * DELTAp[diR[i] - 1] *
TAUp[tiR[i] - 1]);
val += tmp * (diR[i] - ciR[i] * DELTAp[ciR[i]]);
}
/*
* Next do contributions 52 to 54
*/
// Next do contributions 52 to 54
for (j = 0; j < 3; j++) {
i = 52 + j;
doublereal dtmp = delta - epsi[j];
@ -933,9 +887,7 @@ doublereal WaterPropsIAPWSphi::phiR_dt() const
(tiR[i]/tau - 2.0 * betai[j] * ttmp));
}
/*
* Next do contributions 55 and 56
*/
// Next do contributions 55 and 56
for (j = 0; j < 2; j++) {
i = 55 + j;
doublereal deltam1 = delta - 1.0;
@ -977,29 +929,24 @@ doublereal WaterPropsIAPWSphi::dfind(doublereal p_red, doublereal tau, doublerea
doublereal deldd = dd;
doublereal pcheck = 1.0E-30 + 1.0E-8 * p_red;
for (int n = 0; n < 200; n++) {
/*
* Calculate the internal polynomials, and then calculate the
* phi deriv functions needed by this routine.
*/
// Calculate the internal polynomials, and then calculate the phi deriv
// functions needed by this routine.
tdpolycalc(tau, dd);
doublereal q1 = phiR_d();
doublereal q2 = phiR_dd();
/*
* Calculate the predicted reduced pressure, pred0, based on the
* current tau and dd.
*/
// Calculate the predicted reduced pressure, pred0, based on the current
// tau and dd.
doublereal pred0 = dd + dd * dd * q1;
/*
* Calculate the derivative of the predicted reduced pressure
* wrt the reduced density, dd, This is dpddelta
*/
// Calculate the derivative of the predicted reduced pressure wrt the
// reduced density, dd, This is dpddelta
doublereal dpddelta = 1.0 + 2.0 * dd * q1 + dd * dd * q2;
/*
* If dpddelta is negative, then we are in the middle of the
* 2 phase region, beyond the stability curve. We need to adjust
* the initial guess outwards and start a new iteration.
*/
// If dpddelta is negative, then we are in the middle of the 2 phase
// region, beyond the stability curve. We need to adjust the initial
// guess outwards and start a new iteration.
if (dpddelta <= 0.0) {
if (deltaGuess > 1.0) {
dd = dd * 1.05;
@ -1009,50 +956,41 @@ doublereal WaterPropsIAPWSphi::dfind(doublereal p_red, doublereal tau, doublerea
}
continue;
}
/*
* Check for convergence
*/
// Check for convergence
if (fabs(pred0-p_red) < pcheck) {
conv = true;
break;
}
/*
* Dampen and crop the update
*/
// Dampen and crop the update
doublereal dpdx = dpddelta;
if (n < 10) {
dpdx = dpddelta * 1.1;
}
dpdx = std::max(dpdx, 0.001);
/*
* Formulate the update to reduced density using
* Newton's method. Then, crop it to a max value
* of 0.02
*/
// Formulate the update to reduced density using Newton's method. Then,
// crop it to a max value of 0.02
deldd = - (pred0 - p_red) / dpdx;
if (fabs(deldd) > 0.05) {
deldd = deldd * 0.05 / fabs(deldd);
}
/*
* updated the reduced density value
*/
// updated the reduced density value
dd += deldd;
if (fabs(deldd/dd) < 1.0E-14) {
conv = true;
break;
}
/*
* Check for negative densities
*/
// Check for negative densities
if (dd <= 0.0) {
dd = 1.0E-24;
}
}
/*
* Check for convergence, and return 0.0 if it wasn't achieved.
*/
// Check for convergence, and return 0.0 if it wasn't achieved.
if (! conv) {
dd = 0.0;
}

45
src/thermo/WaterSSTP.cpp
View file

@ -56,10 +56,8 @@ WaterSSTP::WaterSSTP(const WaterSSTP& b) :
{
m_waterProps.reset(new WaterProps(&m_sub));
/*
* Use the assignment operator to do the brunt
* of the work for the copy constructor.
*/
// Use the assignment operator to do the brunt of the work for the copy
// constructor.
*this = b;
}
@ -84,18 +82,14 @@ ThermoPhase* WaterSSTP::duplMyselfAsThermoPhase() const
void WaterSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
{
/*
* Do initializations that don't depend on knowing the XML file
*/
// Do initializations that don't depend on knowing the XML file
initThermo();
/*
* Calculate the molecular weight. Note while there may
* be a very good calculated weight in the steam table
* class, using this weight may lead to codes exhibiting
* mass loss issues. We need to grab the elemental
* atomic weights used in the Element class and calculate
* a consistent H2O molecular weight based on that.
*/
// Calculate the molecular weight. Note while there may be a very good
// calculated weight in the steam table class, using this weight may lead to
// codes exhibiting mass loss issues. We need to grab the elemental atomic
// weights used in the Element class and calculate a consistent H2O
// molecular weight based on that.
size_t nH = elementIndex("H");
if (nH == npos) {
throw CanteraError("WaterSSTP::initThermo",
@ -113,9 +107,7 @@ void WaterSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
double one = 1.0;
setMoleFractions(&one);
/*
* Set the baseline
*/
// Set the baseline
doublereal T = 298.15;
Phase::setDensity(7.0E-8);
Phase::setTemperature(T);
@ -140,25 +132,18 @@ void WaterSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
}
h = enthalpy_mole();
/*
* Set the initial state of the system to 298.15 K and
* 1 bar.
*/
// Set the initial state of the system to 298.15 K and 1 bar.
setTemperature(298.15);
double rho0 = m_sub.density(298.15, OneAtm, WATER_LIQUID);
setDensity(rho0);
m_waterProps.reset(new WaterProps(&m_sub));
/*
* We have to do something with the thermo function here.
*/
// We have to do something with the thermo function here.
delete m_spthermo;
m_spthermo = 0;
/*
* Set the flag to say we are ready to calculate stuff
*/
// Set the flag to say we are ready to calculate stuff
m_ready = true;
}
@ -408,9 +393,7 @@ doublereal WaterSSTP::vaporFraction() const
}
return 1.0;
}
/*
* If below tcrit we always return 0 from this class
*/
// If below tcrit we always return 0 from this class
return 0.0;
}