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

A reworking of the LatticePhase and LatticeSolidPhase classes.

Browse source
This commit is contained in:
Harry Moffat 2011-03-11 00:10:50 +00:00
parent c3847347e8
commit f3bec02808
24 changed files with 1436 additions and 291 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
Cantera/src/base/xml.cpp
View file

@ -52,9 +52,8 @@ namespace Cantera {
*/
XML_Error(int line=0) :
m_line(line),
m_msg(0)
m_msg("Error in XML file")
{
m_msg = "Error in XML file";
if (line > 0) {
m_msg += " at line " + int2str(line+1);
}

1
Cantera/src/equil/vcs_VolPhase.cpp
View file

@ -1571,6 +1571,7 @@ namespace VCSnonideal {
for (eT = 0; eT < nebase; eT++) {
ename = tPhase->elementName(eT);
m_elementNames[e] = ename;
m_elementType[e] = tPhase->elementType(eT);
e++;
}

33
Cantera/src/equil/vcs_defs.h
View file

@ -314,6 +314,11 @@ namespace VCSnonideal {
* constraint to one category.
* @{
*/
//! An element constraint that is current turned off
#define VCS_ELEM_TYPE_TURNEDOFF -1
//! Normal element constraint consisting of positive coefficients for the
//! formula matrix.
/*!
@ -335,11 +340,37 @@ namespace VCSnonideal {
*/
#define VCS_ELEM_TYPE_CHARGENEUTRALITY 2
//! Constraint associated with maintaing a fixed lattice stoichiometry int eh
//! solids
/*!
* The constraint may have positive or negative values. The lattice 0 species will
* have negative values while higher lattices will have positive values
*/
#define VCS_ELEM_TYPE_LATTICERATIO 3
//! Constraint associated with maintaining frozen kinetic equilibria in
//! some functional groups within molecules
/*!
* We seek here to say that some functional groups or ionic states should be
* treated as if they are separate elements given the time scale of the problem.
* This will be abs positive constraint. We have not implemented any examples yet.
* A requirement will be that we must be able to add and subtract these contraints.
*/
#define VCS_ELEM_TYPE_KINETICFROZEN 4
//! Constraint associated with the maintenance of a surface phase
/*!
* We don't have any examples of this yet either. However, surfaces only exist
* because they are interfaces between bulk layers. If we want to treat surfaces
* within thermodynamic systems we must come up with a way to constrain their total
* number.
*/
#define VCS_ELEM_TYPE_SURFACECONSTRAINT 5
//! Other constraint equations
/*!
* currently there are none
*/
#define VCS_ELEM_TYPE_OTHERCONSTRAINT 3
#define VCS_ELEM_TYPE_OTHERCONSTRAINT 6
//@}
/*!

25
Cantera/src/equil/vcs_solve.cpp
View file

@ -518,9 +518,18 @@ namespace VCSnonideal {
* FormulaMatrix[] -> Copy the formula matrix over
*/
for (i = 0; i < nspecies; i++) {
bool nonzero = false;
for (j = 0; j < nelements; j++) {
if (pub->FormulaMatrix[j][i] != 0.0) {
nonzero = true;
}
m_formulaMatrix[j][i] = pub->FormulaMatrix[j][i];
}
if (!nonzero) {
plogf("vcs_prob_specifyFully:: species %d %s has a zero formula matrix!\n", i,
pub->SpName[i].c_str());
return VCS_PUB_BAD;
}
}
/*
@ -573,17 +582,31 @@ namespace VCSnonideal {
/*
* Formulate the Goal Element Abundance Vector
*/
double sum;
if (pub->gai.size() != 0) {
for (i = 0; i < nelements; i++) m_elemAbundancesGoal[i] = pub->gai[i];
for (i = 0; i < nelements; i++) {
m_elemAbundancesGoal[i] = pub->gai[i];
if (pub->m_elType[i] == VCS_ELEM_TYPE_LATTICERATIO) {
if (m_elemAbundancesGoal[i] < 1.0E-10) {
m_elemAbundancesGoal[i] = 0.0;
}
}
}
} else {
if (m_doEstimateEquil == 0) {
for (j = 0; j < nelements; j++) {
m_elemAbundancesGoal[j] = 0.0;
for (kspec = 0; kspec < nspecies; kspec++) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
sum += m_molNumSpecies_old[kspec];
m_elemAbundancesGoal[j] += m_formulaMatrix[j][kspec] * m_molNumSpecies_old[kspec];
}
}
if (pub->m_elType[j] == VCS_ELEM_TYPE_LATTICERATIO) {
if (m_elemAbundancesGoal[j] < 1.0E-10 * sum) {
m_elemAbundancesGoal[j] = 0.0;
}
}
}
} else {
plogf("%sElement Abundances, m_elemAbundancesGoal[], not specified\n", ser);

29
Cantera/src/numerics/ResidEval.h
View file

@ -31,7 +31,7 @@ namespace Cantera {
* Virtual base class for DAE residual function evaluators.
* Classes derived from ResidEval evaluate the residual function
* \f[
\vec{F}(t,\vec{y}, \vec{y^\prime})
* \vec{F}(t,\vec{y}, \vec{y^\prime})
* \f]
* The DAE solver attempts to find a solution y(t) such that F = 0.
* @ingroup DAE_Group
@ -59,13 +59,28 @@ namespace Cantera {
return c_NONE;
}
//! Initialization function
virtual void initSizes()
{
int neq = nEquations();
m_alg.resize(neq, 0);
}
/**
* Specify that solution component k is purely algebraic -
* that is, the derivative of this component does not appear
* in the residual function.
*/
virtual void setAlgebraic(const int k) { m_alg[k] = 1; }
virtual bool isAlgebraic(const int k) {return (m_alg[k] == 1); }
virtual void setAlgebraic(const int k) {
if ((int) m_alg.size() < (k+1)) {
initSizes();
}
m_alg[k] = 1;
}
virtual bool isAlgebraic(const int k) {
return (m_alg[k] == 1);
}
/**
@ -106,6 +121,7 @@ namespace Cantera {
*/
virtual int getInitialConditions(const doublereal t0, doublereal * const y,
doublereal * const ydot) {
initSizes();
throw CanteraError("ResidEval::GetInitialConditions()", "base class called");
return 1;
}
@ -143,7 +159,12 @@ namespace Cantera {
protected:
std::map<int, int> m_alg;
//! Mapping vector that stores whether a degree of freedom is a DAE or not
/*!
* The first index is the equation number. The second index is 1 if it is a DAE,
* and zero if it is not.
*/
std::vector<int> m_alg;
std::map<int, int> m_constrain;
private:

6
Cantera/src/numerics/ResidJacEval.h
View file

@ -60,7 +60,6 @@ namespace Cantera {
public:
//!Default constructor
/*!
* @param atol Initial value of the global tolerance (defaults to 1.0E-13)
@ -208,7 +207,6 @@ namespace Cantera {
const doublereal * const y,
const doublereal * const ydot);
//! Return a vector of delta y's for calculation of the numerical Jacobian
/*!
* There is a default algorithm provided.
@ -233,7 +231,6 @@ namespace Cantera {
doublereal * const delta_y,
const doublereal * const solnWeights = 0);
//! Returns a vector of column scale factors that can be used to column scale Jacobians.
/*!
* Default to yScales[] = 1.0
@ -322,7 +319,6 @@ namespace Cantera {
SquareMatrix &J,
doublereal * const resid);
protected:
//! constant value of atol
@ -335,5 +331,3 @@ namespace Cantera {
#endif

71
Cantera/src/thermo/Constituents.cpp
View file

@ -118,6 +118,9 @@ namespace Cantera {
return m_Elements->atomicNumber(m);
}
int Constituents::elementType(int m) const{
return m_Elements->elementType(m);
}
/*
* Add an element to the set.
@ -155,9 +158,9 @@ namespace Cantera {
*/
void Constituents::
addUniqueElement(const std::string& symbol, doublereal weight,
int atomicNumber, doublereal entropy298)
int atomicNumber, doublereal entropy298, int elem_type)
{
m_Elements->addUniqueElement(symbol, weight, atomicNumber, entropy298);
m_Elements->addUniqueElement(symbol, weight, atomicNumber, entropy298, elem_type);
}
void Constituents::
@ -291,10 +294,10 @@ namespace Cantera {
m_speciesNames.push_back(name);
m_speciesCharge.push_back(charge);
m_speciesSize.push_back(size);
int m_mm = m_Elements->nElements();
int ne = m_Elements->nElements();
// Create a changeable copy of the element composition. We now change the charge potentially
vector_fp compNew(m_mm);
for (int m = 0; m < m_mm; m++) {
vector_fp compNew(ne);
for (int m = 0; m < ne; m++) {
compNew[m] = comp[m];
}
double wt = 0.0;
@ -313,30 +316,17 @@ namespace Cantera {
}
}
} else {
m_Elements->m_elementsFrozen = false;
addUniqueElement("E", 0.000545, 0, 0.0);
m_Elements->m_elementsFrozen = true;
m_mm = m_Elements->nElements();
if (m_kk > 0) {
vector_fp old(m_speciesComp);
m_speciesComp.resize(m_kk*m_mm, 0.0);
for (int k = 0; k < m_kk; k++) {
int m_old = m_mm - 1;
for (int m = 0; m < m_old; m++) {
m_speciesComp[k * m_mm + m] = old[k * (m_old) + m];
}
m_speciesComp[k * (m_mm) + (m_mm-1)] = 0.0;
}
}
addUniqueElementAfterFreeze("E", 0.000545, 0, 0.0, CT_ELEM_TYPE_ELECTRONCHARGE);
ne = m_Elements->nElements();
eindex = m_Elements->elementIndex("E");
compNew.resize(m_mm);
compNew[m_mm-1] = - charge;
compNew.resize(ne);
compNew[ne - 1] = - charge;
//comp[eindex] = -charge;
// throw CanteraError("Constituents::addSpecies",
// "Element List doesn't include E, yet this species has charge:" + name);
}
}
for (int m = 0; m < m_mm; m++) {
for (int m = 0; m < ne; m++) {
m_speciesComp.push_back(compNew[m]);
wt += compNew[m] * aw[m];
}
@ -470,6 +460,8 @@ namespace Cantera {
return m_speciesComp[m_mm * k + m];
}
//====================================================================================================================
/*
*
* getAtoms()
@ -485,6 +477,39 @@ namespace Cantera {
}
}
//====================================================================================================================
int Constituents::addUniqueElementAfterFreeze(const std::string& symbol, doublereal weight, int atomicNumber,
doublereal entropy298, int elem_type)
{
int ii = elementIndex(symbol);
if (ii != -1) {
return ii;
}
// Check to see that the element isn't really in the list
m_Elements->m_elementsFrozen = false;
addUniqueElement(symbol, weight, atomicNumber, entropy298, elem_type);
m_Elements->m_elementsFrozen = true;
int m_mm = m_Elements->nElements();
ii = elementIndex(symbol);
if (ii != m_mm-1) {
throw CanteraError("Constituents::addElementAfterFreeze()", "confused");
}
if (m_kk > 0) {
vector_fp old(m_speciesComp);
m_speciesComp.resize(m_kk*m_mm, 0.0);
for (int k = 0; k < m_kk; k++) {
int m_old = m_mm - 1;
for (int m = 0; m < m_old; m++) {
m_speciesComp[k * m_mm + m] = old[k * (m_old) + m];
}
m_speciesComp[k * (m_mm) + (m_mm-1)] = 0.0;
}
}
return ii;
}
//====================================================================================================================
/*
* This copy constructor just calls the assignment operator
* for this class.

22
Cantera/src/thermo/Constituents.h
View file

@ -110,7 +110,6 @@ namespace Cantera {
/// exception, ElementRangeError, is thrown.
std::string elementName(int m) const;
/// Index of element named 'name'.
/// The index is an integer
/// assigned to each element in the order it was added,
@ -141,6 +140,8 @@ namespace Cantera {
*/
int atomicNumber(int m) const;
int elementType(int m) const;
/// Return a read-only reference to the vector of element names.
const std::vector<std::string>& elementNames() const;
@ -193,7 +194,7 @@ namespace Cantera {
*/
void addUniqueElement(const std::string& symbol, doublereal weight,
int atomicNumber = 0,
doublereal entropy298 = ENTROPY298_UNKNOWN);
doublereal entropy298 = ENTROPY298_UNKNOWN, int elem_type = CT_ELEM_TYPE_ABSPOS);
//! Adde an element, checking for uniqueness
/*!
@ -216,6 +217,23 @@ namespace Cantera {
/// True if freezeElements has been called.
bool elementsFrozen();
//! Add an element after the elements have been frozen, checking for uniqueness
/*!
* The uniqueness is checked by comparing the string symbol. If
* not unique, nothing is done.
*
* @param symbol String symbol of the element
* @param weight Atomic weight of the element (kg kmol-1).
* @param atomicNumber Atomic number of the element (unitless)
* @param entropy298 Entropy of the element at 298 K and 1 bar
* in its most stable form. The default is
* the value ENTROPY298_UNKNOWN, which is
* interpreted as an unknown, and if used
* will cause Cantera to throw an error.
*/
int addUniqueElementAfterFreeze(const std::string& symbol, doublereal weight, int atomicNumber,
doublereal entropy298 = ENTROPY298_UNKNOWN, int elem_type = CT_ELEM_TYPE_ABSPOS);
//@}
/// Returns the number of species in the phase

7
Cantera/src/thermo/DebyeHuckel.h
View file

@ -1255,13 +1255,6 @@ namespace Cantera {
* -------------- Utilities -------------------------------
*/
/**
* Return a reference to the species thermodynamic property
* manager. @todo This method will fail if no species thermo
* manager has been installed.
*/
SpeciesThermo& speciesThermo() { return *m_spthermo; }
//! Initialize the object's internal lengths after species are set
/**

74
Cantera/src/thermo/Elements.cpp
View file

@ -222,6 +222,7 @@ namespace Cantera {
Elements::Elements() :
m_mm(0),
m_elementsFrozen(false),
m_elem_type(0),
numSubscribers(0)
{
}
@ -256,7 +257,7 @@ namespace Cantera {
m_atomicNumbers = right.m_atomicNumbers;
m_elementNames = right.m_elementNames;
m_entropy298 = right.m_entropy298;
m_elem_type = right.m_elem_type;
numSubscribers = 0;
return *this;
@ -335,7 +336,43 @@ namespace Cantera {
AssertTrace(m >= 0 && m < m_mm);
return (m_entropy298[m]);
}
//====================================================================================================================
//! Return the element constraint type
/*!
* Possible types include:
*
* CT_ELEM_TYPE_TURNEDOFF -1
* CT_ELEM_TYPE_ABSPOS 0
* CT_ELEM_TYPE_ELECTRONCHARGE 1
* CT_ELEM_TYPE_CHARGENEUTRALITY 2
* CT_ELEM_TYPE_LATTICERATIO 3
* CT_ELEM_TYPE_KINETICFROZEN 4
* CT_ELEM_TYPE_SURFACECONSTRAINT 5
* CT_ELEM_TYPE_OTHERCONSTRAINT 6
*
* The default is CT_ELEM_TYPE_ABSPOS
*/
int Elements::elementType(int m) const
{
return m_elem_type[m];
}
//====================================================================================================================
// Change the element type of the mth constraint
/*
* Reassigns an element type
*
* @param m Element index
* @param elem_type New elem type to be assigned
*
* @return Returns the old element type
*/
int Elements::changeElementType(int m, int elem_type)
{
int old = m_elem_type[m];
m_elem_type[m] = elem_type;
return old;
}
//====================================================================================================================
/*
*
* Add an element to the current set of elements in the current object.
@ -367,16 +404,22 @@ namespace Cantera {
#ifdef USE_DGG_CODE
m_definedElements[symbol] = nElements() + 1;
#endif
if (symbol == "E") {
m_elem_type.push_back(CT_ELEM_TYPE_ELECTRONCHARGE);
} else {
m_elem_type.push_back(CT_ELEM_TYPE_ABSPOS);
}
m_mm++;
}
//===========================================================================================================
void Elements::
addElement(const XML_Node& e) {
doublereal weight = atof(e["atomicWt"].c_str());
string symbol = e["name"];
addElement(symbol, weight);
}
//===========================================================================================================
/*
* addUniqueElement():
*
@ -393,7 +436,7 @@ namespace Cantera {
#ifdef USE_DGG_CODE
void Elements::
addUniqueElement(const std::string& symbol, doublereal weight, int atomicNumber,
doublereal entropy298)
doublereal entropy298, int elem_type)
{
if (m_elementsFrozen)
throw ElementsFrozen("addElement");
@ -413,13 +456,17 @@ namespace Cantera {
m_elementNames.push_back(symbol);
m_atomicNumbers.push_back(atomicNumber);
m_entropy298.push_back(entropy298);
if (symbol == "E") {
m_elem_type.push_back(CT_ELEM_TYPE_ELECTRONCHARGE);
} else {
m_elem_type.push_back(elem_type);
}
m_mm++;
}
else {
if (m_atomicWeights[i] != weight) {
throw CanteraError("AddUniqueElement",
"Duplicate Elements (" + symbol +
") have different weights");
"Duplicate Elements (" + symbol + ") have different weights");
}
}
}
@ -427,7 +474,8 @@ namespace Cantera {
#else
void Elements::
addUniqueElement(const std::string& symbol,
doublereal weight, int atomicNumber, doublereal entropy298)
doublereal weight, int atomicNumber, doublereal entropy298,
int elem_type)
{
if (weight == -12345.0) {
weight = LookupWtElements(symbol);
@ -458,12 +506,16 @@ namespace Cantera {
m_elementNames.push_back(symbol);
m_atomicNumbers.push_back(atomicNumber);
m_entropy298.push_back(entropy298);
if (symbol == "E") {
m_elem_type.push_back(CT_ELEM_TYPE_ELECTRONCHARGE);
} else {
m_elem_type.push_back(elem_type);
}
m_mm++;
} else {
if (m_atomicWeights[i] != weight) {
throw CanteraError("AddUniqueElement",
"Duplicate Elements (" + symbol +
") have different weights");
"Duplicate Elements (" + symbol + ") have different weights");
}
}
}
@ -506,6 +558,8 @@ namespace Cantera {
m_mm = 0;
m_atomicWeights.resize(0);
m_elementNames.resize(0);
m_entropy298.resize(0);
m_elem_type.resize(0);
m_elementsFrozen = false;
}

106
Cantera/src/thermo/Elements.h
View file

@ -26,6 +26,74 @@ namespace Cantera {
class XML_Node;
class ElementRangeError;
/*!
* @name Types of Element Constraint Equations
*
* There may be several different types of element constraints handled
* by the equilibrium program and by Cantera in other contexts.
* These defines are used to assign each constraint to one category.
* @{
*/
//! An element constraint that is current turned off
#define CT_ELEM_TYPE_TURNEDOFF -1
//! Normal element constraint consisting of positive coefficients for the
//! formula matrix.
/*!
* All species have positive coefficients within the formula matrix.
* With this constraint, we may employ various strategies to handle
* small values of the element number successfully.
*/
#define CT_ELEM_TYPE_ABSPOS 0
//! This refers to conservation of electrons
/*!
* Electrons may have positive or negative values in the Formula matrix.
*/
#define CT_ELEM_TYPE_ELECTRONCHARGE 1
//! This refers to a charge neutrality of a single phase
/*!
* Charge neutrality may have positive or negative values in the Formula matrix.
*/
#define CT_ELEM_TYPE_CHARGENEUTRALITY 2
//! Constraint associated with maintaing a fixed lattice stoichiometry int eh
//! solids
/*!
* The constraint may have positive or negative values. The lattice 0 species will
* have negative values while higher lattices will have positive values
*/
#define CT_ELEM_TYPE_LATTICERATIO 3
//! Constraint associated with maintaining frozen kinetic equilibria in
//! some functional groups within molecules
/*!
* We seek here to say that some functional groups or ionic states should be
* treated as if they are separate elements given the time scale of the problem.
* This will be abs positive constraint. We have not implemented any examples yet.
* A requirement will be that we must be able to add and subtract these contraints.
*/
#define CT_ELEM_TYPE_KINETICFROZEN 4
//! Constraint associated with the maintenance of a surface phase
/*!
* We don't have any examples of this yet either. However, surfaces only exist
* because they are interfaces between bulk layers. If we want to treat surfaces
* within thermodynamic systems we must come up with a way to constrain their total
* number.
*/
#define CT_ELEM_TYPE_SURFACECONSTRAINT 5
//! Other constraint equations
/*!
* currently there are none
*/
#define CT_ELEM_TYPE_OTHERCONSTRAINT 6
//@}
//! Positive number indicating we don't know the gibbs free energy
//! of the element in its most stable state at 298.15 K and 1 bar.
//#define GIBSSFE298_UNKNOWN 123456789.
@ -120,6 +188,37 @@ namespace Cantera {
*/
doublereal entropyElement298(int m) const;
//! Return the element constraint type
/*!
* Possible types include:
*
* CT_ELEM_TYPE_ABSPOS 0
* CT_ELEM_TYPE_ELECTRONCHARGE 1
* CT_ELEM_TYPE_CHARGENEUTRALITY 2
* CT_ELEM_TYPE_LATTICERATIO 3
* CT_ELEM_TYPE_KINETICFROZEN 4
* CT_ELEM_TYPE_SURFACECONSTRAINT 5
* CT_ELEM_TYPE_OTHERCONSTRAINT 6
*
* The default is CT_ELEM_TYPE_ABSPOS
*
* @param m Element index
*
* @return Returns the element type
*/
int elementType(int m) const;
//! Change the element type of the mth constraint
/*!
* Reassigns an element type
*
* @param m Element index
* @param elem_type New elem type to be assigned
*
* @return Returns the old element type
*/
int changeElementType(int m, int elem_type);
/// vector of element atomic weights
const vector_fp& atomicWeights() const { return m_atomicWeights; }
@ -199,7 +298,7 @@ namespace Cantera {
*/
void addUniqueElement(const std::string& symbol,
doublereal weight = -12345.0, int atomicNumber = 0,
doublereal entropy298 = ENTROPY298_UNKNOWN);
doublereal entropy298 = ENTROPY298_UNKNOWN, int elem_type = CT_ELEM_TYPE_ABSPOS);
//! Add an element to the current set of elements in the current object.
/*!
@ -258,7 +357,7 @@ namespace Cantera {
* If this is true, then no elements may be added to the
* object.
*/
bool m_elementsFrozen;
bool m_elementsFrozen;
/**
* Vector of element atomic weights:
@ -285,6 +384,9 @@ namespace Cantera {
*/
vector_fp m_entropy298;
//! Vector of element types
vector_int m_elem_type;
/**
* Number of Constituents Objects that use this object
*

2
Cantera/src/thermo/FixedChemPotSSTP.cpp
View file

@ -144,6 +144,8 @@ namespace Cantera {
ss.addChild("h", sval);
ss.addChild("s", "0.0");
saveSpeciesData(0, s);
delete s;
s = 0;
}
//====================================================================================================================

8
Cantera/src/thermo/HMWSoln.h
View file

@ -2053,14 +2053,6 @@ namespace Cantera {
* -------------- Utilities -------------------------------
*/
/**
* Return a reference to the species thermodynamic property
* manager.
*
* @todo This method will fail if no species thermo
* manager has been installed.
*/
SpeciesThermo& speciesThermo() { return *m_spthermo; }
//! Initialization of a HMWSoln phase using an xml file
/*!

8
Cantera/src/thermo/IdealMolalSoln.h
View file

@ -777,14 +777,6 @@ namespace Cantera {
/*
* -------------- Utilities -------------------------------
*/
/*!
* Return a reference to the species thermodynamic property
* manager. @todo This method will fail if no species thermo
* manager has been installed.
*/
SpeciesThermo& speciesThermo() { return *m_spthermo; }
//! Initialization routine for an IdealMolalSoln phase.
/*!

11
Cantera/src/thermo/IdealSolidSolnPhase.cpp
View file

@ -117,7 +117,7 @@ namespace Cantera {
IdealSolidSolnPhase *ii = new IdealSolidSolnPhase(*this);
return (ThermoPhase*) ii;
}
//====================================================================================================================
//====================================================================================================================
/**
* Equation of state flag. Returns the value cIdealGas, defined
* in mix_defs.h.
@ -679,8 +679,7 @@ namespace Cantera {
* property manager. They are polynomial functions of temperature.
* @see SpeciesThermo
*/
void IdealSolidSolnPhase::
getPartialMolarEnthalpies(doublereal* hbar) const {
void IdealSolidSolnPhase::getPartialMolarEnthalpies(doublereal* hbar) const {
const array_fp& _h = enthalpy_RT_ref();
doublereal rt = GasConstant * temperature();
scale(_h.begin(), _h.end(), hbar, rt);
@ -892,8 +891,7 @@ namespace Cantera {
* units = m^3 / kmol
*/
void IdealSolidSolnPhase::getStandardVolumes(doublereal *vol) const {
copy(m_speciesMolarVolume.begin(),
m_speciesMolarVolume.end(), vol);
copy(m_speciesMolarVolume.begin(), m_speciesMolarVolume.end(), vol);
}
@ -1227,8 +1225,7 @@ namespace Cantera {
"Unknown standardConc model: " + formStringa);
}
} else {
throw CanteraError(subname.c_str(),
"Unspecified standardConc model");
throw CanteraError(subname.c_str(), "Unspecified standardConc model");
}
/*

9
Cantera/src/thermo/IdealSolidSolnPhase.h
View file

@ -579,9 +579,9 @@ namespace Cantera {
/// @name Partial Molar Properties of the Solution -----------------------------
//@{
/**
* Returns an array of partial molar enthalpies for the species
* in the mixture.
//! Returns an array of partial molar enthalpies for the species in the mixture.
/*!
* Units (J/kmol)
* For this phase, the partial molar enthalpies are equal to the
* pure species enthalpies
@ -1057,8 +1057,9 @@ namespace Cantera {
*/
doublereal m_Pcurrent;
//! Vector of molar volumes for each species in the solution
/**
* Species molar volume \f$ m^3 kmol^-1 \f$
* Species molar volumes \f$ m^3 kmol^-1 \f$
*/
array_fp m_speciesMolarVolume;

263
Cantera/src/thermo/LatticePhase.cpp
View file

@ -28,11 +28,20 @@
#include <cmath>
#include <string>
using namespace std;
namespace Cantera {
// Base Empty constructor
LatticePhase::LatticePhase() :
m_tlast(0.0)
m_mm(0),
m_tmin(0.0),
m_tmax(0.0),
m_Pref(OneAtm),
m_Pcurrent(OneAtm),
m_tlast(0.0),
m_speciesMolarVolume(0),
m_site_density(0.0)
{
}
@ -41,7 +50,14 @@ namespace Cantera {
* @param right Object to be copied
*/
LatticePhase::LatticePhase(const LatticePhase &right) :
m_tlast(0.0)
m_mm(0),
m_tmin(0.0),
m_tmax(0.0),
m_Pref(OneAtm),
m_Pcurrent(OneAtm),
m_tlast(0.0),
m_speciesMolarVolume(0),
m_site_density(0.0)
{
*this = operator=(right);
}
@ -56,15 +72,16 @@ namespace Cantera {
m_mm = right.m_mm;
m_tmin = right.m_tmin;
m_tmax = right.m_tmax;
m_p0 = right.m_p0;
m_Pref = right.m_Pref;
m_Pcurrent = right.m_Pcurrent;
m_tlast = right.m_tlast;
m_h0_RT = right.m_h0_RT;
m_cp0_R = right.m_cp0_R;
m_g0_RT = right.m_g0_RT;
m_s0_R = right.m_s0_R;
m_press = right.m_press;
m_vacancy = right.m_vacancy;
m_molar_density = right.m_molar_density;
m_speciesMolarVolume = right.m_speciesMolarVolume;
m_site_density = right.m_site_density;
}
return *this;
}
@ -195,99 +212,176 @@ namespace Cantera {
return GasConstant * (mean_X(&entropy_R_ref()[0]) -
sum_xlogx());
}
//====================================================================================================================
doublereal LatticePhase::gibbs_mole() const {
return enthalpy_mole() - temperature() * entropy_mole();
}
//====================================================================================================================
doublereal LatticePhase::cp_mole() const {
return GasConstant * mean_X(&cp_R_ref()[0]);
}
//====================================================================================================================
doublereal LatticePhase::cv_mole() const {
return cp_mole();
}
void LatticePhase::setPressure(doublereal p) {
m_press = p;
setMolarDensity(m_molar_density);
//====================================================================================================================
doublereal LatticePhase::calcDensity() {
setMolarDensity(m_site_density);
doublereal mw = meanMolecularWeight();
doublereal dens = mw * m_site_density;
/*
* Calculate the molarVolume of the solution (m**3 kmol-1)
*/
// const doublereal * const dtmp = moleFractdivMMW();
// doublereal invDens = dot(m_speciesMolarVolume.begin(), m_speciesMolarVolume.end(), dtmp);
/*
* Set the density in the parent State object directly,
* by calling the State::setDensity() function.
*/
// doublereal dens = 1.0/invDens;
// State::setDensity(dens);
return dens;
}
//====================================================================================================================
void LatticePhase::setPressure(doublereal p) {
m_Pcurrent = p;
calcDensity();
}
//====================================================================================================================
void LatticePhase::setMoleFractions(const doublereal * const x) {
State::setMoleFractions(x);
calcDensity();
}
//====================================================================================================================
void LatticePhase::setMoleFractions_NoNorm(const doublereal * const x) {
State::setMoleFractions(x);
calcDensity();
}
//====================================================================================================================
void LatticePhase::setMassFractions(const doublereal * const y) {
State::setMassFractions(y);
calcDensity();
}
//====================================================================================================================
void LatticePhase::setMassFractions_NoNorm(const doublereal * const y) {
State::setMassFractions_NoNorm(y);
calcDensity();
}
//====================================================================================================================
void LatticePhase::setConcentrations(const doublereal * const c) {
State::setConcentrations(c);
calcDensity();
}
//====================================================================================================================
void LatticePhase::getActivityConcentrations(doublereal* c) const {
getMoleFractions(c);
}
//====================================================================================================================
void LatticePhase::getActivityCoefficients(doublereal* ac) const {
for (int k = 0; k < m_kk; k++) {
ac[k] = 1.0;
}
}
//====================================================================================================================
doublereal LatticePhase::standardConcentration(int k) const {
return 1.0;
}
//====================================================================================================================
doublereal LatticePhase::logStandardConc(int k) const {
return 0.0;
}
//====================================================================================================================
void LatticePhase::getChemPotentials(doublereal* mu) const {
doublereal vdp = ((pressure() - m_spthermo->refPressure())/
molarDensity());
doublereal delta_p = m_Pcurrent - m_Pref;
doublereal xx;
doublereal rt = temperature() * GasConstant;
doublereal RT = temperature() * GasConstant;
const array_fp& g_RT = gibbs_RT_ref();
for (int k = 0; k < m_kk; k++) {
xx = fmaxx(SmallNumber, moleFraction(k));
mu[k] = rt*(g_RT[k] + log(xx)) + vdp;
mu[k] = RT * (g_RT[k] + log(xx))
+ delta_p * m_speciesMolarVolume[k];
}
}
//====================================================================================================================
void LatticePhase::getPartialMolarEnthalpies(doublereal* hbar) const {
const array_fp& _h = enthalpy_RT_ref();
doublereal rt = GasConstant * temperature();
scale(_h.begin(), _h.end(), hbar, rt);
}
//====================================================================================================================
void LatticePhase::getPartialMolarEntropies(doublereal* sbar) const {
const array_fp& _s = entropy_R_ref();
doublereal r = GasConstant;
doublereal xx;
for (int k = 0; k < m_kk; k++) {
xx = fmaxx(SmallNumber, moleFraction(k));
sbar[k] = r * (_s[k] - log(xx));
}
}
//====================================================================================================================
void LatticePhase::getPartialMolarCp(doublereal* cpbar) const {
getCp_R(cpbar);
for (int k = 0; k < m_kk; k++) {
cpbar[k] *= GasConstant;
}
}
//====================================================================================================================
void LatticePhase::getPartialMolarVolumes(doublereal* vbar) const {
getStandardVolumes(vbar);
}
//====================================================================================================================
void LatticePhase::getStandardChemPotentials(doublereal* mu0) const {
const array_fp& gibbsrt = gibbs_RT_ref();
scale(gibbsrt.begin(), gibbsrt.end(), mu0, _RT());
}
//====================================================================================================================
void LatticePhase::getPureGibbs(doublereal* gpure) const {
const array_fp& gibbsrt = gibbs_RT_ref();
scale(gibbsrt.begin(), gibbsrt.end(), gpure, _RT());
}
void LatticePhase::getEnthalpy_RT(doublereal* hrt) const {
const array_fp& _h = enthalpy_RT_ref();
std::copy(_h.begin(), _h.end(), hrt);
doublereal tmp = (pressure() - m_p0) / (molarDensity() * GasConstant * temperature());
doublereal delta_p = (m_Pcurrent - m_Pref);
double RT = GasConstant * temperature();
for (int k = 0; k < m_kk; k++) {
hrt[k] += tmp;
gpure[k] = RT * gibbsrt[k] + delta_p * m_speciesMolarVolume[k];
}
}
//====================================================================================================================
void LatticePhase::getEnthalpy_RT(doublereal* hrt) const {
const array_fp& _h = enthalpy_RT_ref();
doublereal delta_prt = ((m_Pcurrent - m_Pref) / (GasConstant * temperature()));
for (int k = 0; k < m_kk; k++) {
hrt[k] = _h[k] + delta_prt * m_speciesMolarVolume[k];
}
}
//====================================================================================================================
void LatticePhase::getEntropy_R(doublereal* sr) const {
const array_fp& _s = entropy_R_ref();
std::copy(_s.begin(), _s.end(), sr);
}
//====================================================================================================================
void LatticePhase::getGibbs_RT(doublereal* grt) const {
const array_fp& gibbsrt = gibbs_RT_ref();
std::copy(gibbsrt.begin(), gibbsrt.end(), grt);
doublereal RT = _RT();
doublereal delta_prt = (m_Pcurrent - m_Pref)/ RT;
for (int k = 0; k < m_kk; k++) {
grt[k] = gibbsrt[k] + delta_prt * m_speciesMolarVolume[k];
}
}
//====================================================================================================================
void LatticePhase::getGibbs_ref(doublereal *g) const {
getGibbs_RT_ref(g);
for (int k = 0; k < m_kk; k++) {
g[k] *= GasConstant * temperature();
}
}
//===================================================================================================================
void LatticePhase::getCp_R(doublereal* cpr) const {
const array_fp& _cpr = cp_R_ref();
std::copy(_cpr.begin(), _cpr.end(), cpr);
}
//===================================================================================================================
void LatticePhase::getStandardVolumes(doublereal* vbar) const {
doublereal vv = 1.0/m_molar_density;
for (int k = 0; k < m_kk; k++) {
vbar[k] = vv;
}
copy(m_speciesMolarVolume.begin(), m_speciesMolarVolume.end(), vbar);
}
//=======================================================================================================
// Returns the vector of nondimensional Enthalpies of the reference state at the current temperature
@ -310,6 +404,13 @@ namespace Cantera {
_updateThermo();
return m_g0_RT;
}
//====================================================================================================================
void LatticePhase::getGibbs_RT_ref(doublereal *grt) const {
_updateThermo();
for (int k = 0; k < m_kk; k++) {
grt[k] = m_g0_RT[k];
}
}
//=======================================================================================================
// Returns a reference to the dimensionless reference state Entropy vector.
/*
@ -330,7 +431,21 @@ namespace Cantera {
_updateThermo();
return m_cp0_R;
}
//=======================================================================================================
//====================================================================================================================
// Initialize the ThermoPhase object after all species have been set up
/*
* @internal Initialize.
*
* This method performs any initialization required after all
* species have been added. For example, it is used to
* resize internal work arrays that must have an entry for
* each species.
* This method is called from ThermoPhase::initThermoXML(),
* which is called from importPhase(),
* just prior to returning from the function, importPhase().
*
* @see importCTML.cpp
*/
void LatticePhase::initThermo() {
m_kk = nSpecies();
m_mm = nElements();
@ -338,14 +453,62 @@ namespace Cantera {
doublereal tmax = m_spthermo->maxTemp();
if (tmin > 0.0) m_tmin = tmin;
if (tmax > 0.0) m_tmax = tmax;
m_p0 = refPressure();
m_Pref = refPressure();
int leng = m_kk;
m_h0_RT.resize(leng);
m_g0_RT.resize(leng);
m_cp0_R.resize(leng);
m_s0_R.resize(leng);
setMolarDensity(m_molar_density);
m_speciesMolarVolume.resize(leng, 0.0);
ThermoPhase::initThermo();
}
//====================================================================================================================
void LatticePhase::initThermoXML(XML_Node& phaseNode, std::string id) {
std::string subname = "LatticePhase::initThermoXML";
/*
* Check on the thermo field. Must have:
* <thermo model="Lattice" />
*/
if (phaseNode.hasChild("thermo")) {
XML_Node& thNode = phaseNode.child("thermo");
std::string mStringa = thNode.attrib("model");
std::string mString = lowercase(mStringa);
if (mString != "lattice") {
throw CanteraError(subname.c_str(),
"Unknown thermo model: " + mStringa);
}
} else {
throw CanteraError(subname.c_str(),
"Unspecified thermo model");
}
/*
* 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());
const std::vector<std::string> &sss = speciesNames();
for (int k = 0; k < m_kk; k++) {
m_speciesMolarVolume[k] = m_site_density;
XML_Node* s = speciesDB->findByAttr("name", sss[k]);
if (!s) {
throw CanteraError(" LatticePhase::initThermoXML", "database problems");
}
XML_Node *ss = s->findByName("standardState");
if (ss) {
if (ss->findByName("molarVolume")) {
m_speciesMolarVolume[k] = getFloat(*ss, "molarVolume", "toSI");
}
}
}
/*
* Call the base initThermo, which handles setting the initial
* state.
*/
ThermoPhase::initThermoXML(phaseNode, id);
}
//=====================================================================================================
// Update the species reference state thermodynamic functions
@ -355,10 +518,6 @@ namespace Cantera {
*/
void LatticePhase::_updateThermo() const {
doublereal tnow = temperature();
if (fabs(molarDensity() - m_molar_density)/m_molar_density > 0.0001) {
throw CanteraError("_updateThermo","molar density changed from "
+fp2str(m_molar_density)+" to "+fp2str(molarDensity()));
}
if (m_tlast != tnow) {
m_spthermo->update(tnow, &m_cp0_R[0], &m_h0_RT[0], &m_s0_R[0]);
m_tlast = tnow;
@ -370,8 +529,8 @@ namespace Cantera {
}
//=====================================================================================================
void LatticePhase::setParameters(int n, doublereal* const c) {
m_molar_density = c[0];
setMolarDensity(m_molar_density);
m_site_density = c[0];
setMolarDensity(m_site_density);
}
//=====================================================================================================
void LatticePhase::getParameters(int &n, doublereal * const c) const {
@ -382,7 +541,7 @@ namespace Cantera {
//=====================================================================================================
void LatticePhase::setParametersFromXML(const XML_Node& eosdata) {
eosdata._require("model", "Lattice");
m_molar_density = getFloat(eosdata, "site_density", "toSI");
m_site_density = getFloat(eosdata, "site_density", "toSI");
m_vacancy = getChildValue(eosdata, "vacancy_species");
}
//=====================================================================================================

213
Cantera/src/thermo/LatticePhase.h
View file

@ -453,7 +453,7 @@ namespace Cantera {
* independent value of the pressure.
*/
virtual doublereal pressure() const {
return m_press;
return m_Pcurrent;
}
//! Set the internally storred pressure (Pa) at constant
@ -465,7 +465,67 @@ namespace Cantera {
* @param p Input Pressure (Pa)
*/
virtual void setPressure(doublereal p);
//! Calculate the density of the mixture using the partial
//! molar volumes and mole fractions as input
/*!
* The formula for this is
*
* \f[
* \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}}
* \f]
*
* where \f$X_k\f$ are the mole fractions, \f$W_k\f$ are
* the molecular weights, and \f$V_k\f$ are the pure species
* molar volumes.
*
* Note, the basis behind this formula is that in an ideal
* solution the partial molar volumes are equal to the pure
* species molar volumes. We have additionally specified
* in this class that the pure species molar volumes are
* independent of temperature and pressure.
*
* NOTE: This is a non-virtual function, which is not a
* member of the ThermoPhase base class.
*/
doublereal calcDensity();
//! Set the mole fractions
/*!
* @param x Input vector of mole fractions.
* Length: m_kk.
*/
virtual void setMoleFractions(const doublereal * const x);
//! Set the mole fractions, but don't normalize them to one.
/*!
* @param x Input vector of mole fractions.
* Length: m_kk.
*/
virtual void setMoleFractions_NoNorm(const doublereal * const x);
//! Set the mass fractions, and normalize them to one.
/*!
* @param y Input vector of mass fractions.
* Length: m_kk.
*/
virtual void setMassFractions(const doublereal * const y);
//! Set the mass fractions, but don't normalize them to one
/*!
* @param y Input vector of mass fractions.
* Length: m_kk.
*/
virtual void setMassFractions_NoNorm(const doublereal * const y);
//! Set the concentration,
/*!
* @param c Input vector of concentrations.
* Length: m_kk.
*/
virtual void setConcentrations(const doublereal * const c);
//@}
/// @name Activities, Standard States, and Activity Concentrations
/**
@ -547,6 +607,71 @@ namespace Cantera {
*/
virtual void getChemPotentials(doublereal* mu) const;
//@}
/// @name Partial Molar Properties of the Solution -----------------------------
//@{
/**
* Returns an array of partial molar enthalpies for the species
* in the mixture.
* Units (J/kmol)
* For this phase, the partial molar enthalpies are equal to the
* pure species enthalpies
* \f[
* \bar h_k(T,P) = \hat h^{ref}_k(T) + (P - P_{ref}) \hat V^0_k
* \f]
* The reference-state pure-species enthalpies, \f$ \hat h^{ref}_k(T) \f$,
* at the reference pressure,\f$ P_{ref} \f$,
* are computed by the species thermodynamic
* property manager. They are polynomial functions of temperature.
* @see SpeciesThermo
*
* @param hbar Output vector containing partial molar enthalpies.
* Length: m_kk.
*/
virtual void getPartialMolarEnthalpies(doublereal* hbar) const;
/**
* Returns an array of partial molar entropies of the species in the
* solution. Units: J/kmol/K.
* For this phase, the partial molar entropies are equal to the
* pure species entropies plus the ideal solution contribution.
* \f[
* \bar s_k(T,P) = \hat s^0_k(T) - R log(X_k)
* \f]
* The reference-state pure-species entropies,\f$ \hat s^{ref}_k(T) \f$,
* at the reference pressure, \f$ P_{ref} \f$, are computed by the
* species thermodynamic
* property manager. They are polynomial functions of temperature.
* @see SpeciesThermo
*
* @param sbar Output vector containing partial molar entropies.
* Length: m_kk.
*/
virtual void getPartialMolarEntropies(doublereal* sbar) const;
/**
* Returns an array of partial molar Heat Capacities at constant
* pressure of the species in the
* solution. Units: J/kmol/K.
* For this phase, the partial molar heat capacities are equal
* to the standard state heat capacities.
*
* @param cpbar Output vector of partial heat capacities. Length: m_kk.
*/
virtual void getPartialMolarCp(doublereal* cpbar) const;
//! Return an array of partial molar volumes for the
//! species in the mixture. Units: m^3/kmol.
/*!
* @param vbar Output vector of speciar partial molar volumes.
* Length = m_kk. units are m^3/kmol.
*/
virtual void getPartialMolarVolumes(doublereal* vbar) const;
//! Get the array of chemical potentials at unit activity for the
//! species standard states at the current <I>T</I> and <I>P</I> of the solution.
/*!
@ -568,14 +693,7 @@ namespace Cantera {
*/
virtual void getPureGibbs(doublereal* gpure) const;
//! Return an array of partial molar volumes for the
//! species in the mixture. Units: m^3/kmol.
/*!
* @param vbar Output vector of speciar partial molar volumes.
* Length = m_kk. units are m^3/kmol.
*/
virtual void getPartialMolarVolumes(doublereal* vbar) const;
//@}
/// @name Properties of the Standard State of the Species in the Solution
//@{
@ -702,6 +820,25 @@ namespace Cantera {
*/
const array_fp& gibbs_RT_ref() const;
//! Returns the vector of nondimensional
//! Gibbs Free Energies of the reference state at the current temperature
//! of the solution and the reference pressure for the species.
/*!
* @param grt Output vector containing the nondimensional reference state
* Gibbs Free energies. Length: m_kk.
*/
virtual void getGibbs_RT_ref(doublereal *grt) const;
//! Returns the vector of the gibbs function of the reference state at the current temperature
//! of the solution and the reference pressure for the species.
/*!
* units = J/kmol
*
* @param g Output vector containing the reference state
* Gibbs Free energies. Length: m_kk. Units: J/kmol.
*/
virtual void getGibbs_ref(doublereal *g) const;
//! Returns a reference to the dimensionless reference state Entropy vector.
/*!
* This function is part of the layer that checks/recalculates the reference
@ -736,6 +873,33 @@ namespace Cantera {
*/
virtual void initThermo();
//! Import and initialize a ThermoPhase object using an XML tree.
/*!
* Here we read extra information about the XML description
* of a phase. Regular information about elements and species
* and their reference state thermodynamic information
* have already been read at this point.
* For example, we do not need to call this function for
* ideal gas equations of state.
* This function is called from importPhase()
* after the elements and the
* species are initialized with default ideal solution
* level data.
*
* @param phaseNode This object must be the phase node of a
* complete XML tree
* description of the phase, including all of the
* species data. In other words while "phase" must
* point to an XML phase object, it must have
* sibling nodes "speciesData" that describe
* the species in the phase.
* @param id ID of the phase. If nonnull, a check is done
* to see if phaseNode is pointing to the phase
* with the correct id.
*/
virtual void initThermoXML(XML_Node& phaseNode, std::string id);
//! Set the equation of state parameters from the argument list
/*!
* @internal
@ -792,6 +956,7 @@ namespace Cantera {
protected:
//! Number of elements
int m_mm;
@ -810,7 +975,17 @@ namespace Cantera {
doublereal m_tmax;
//! Reference state pressure
doublereal m_p0;
doublereal m_Pref;
//! The current pressure
/*!
* Since the density isn't a function of pressure, but only of the
* mole fractions, we need to independently specify the pressure.
* The density variable which is inherited as part of the State class,
* m_dens, is always kept current whenever T, P, or X[] change.
*/
doublereal m_Pcurrent;
//! Current value of the temperature (Kelvin)
mutable doublereal m_tlast;
@ -827,8 +1002,6 @@ namespace Cantera {
//! Temporary storage for the reference state entropies at the current temperature
mutable array_fp m_s0_R;
//! Current value of the pressure (Pa)
doublereal m_press;
//! String name for the species which represents a vacency
//! in the lattice
@ -837,13 +1010,21 @@ namespace Cantera {
*/
std::string m_vacancy;
//! Molar density of the lattice solid
//! Vector of molar volumes for each species in the solution
/**
* Species molar volumes \f$ m^3 kmol^-1 \f$
*/
array_fp m_speciesMolarVolume;
//! Site Density of the lattice solid
/*!
* Currently, this does not change as a function of T, P or composition
* Currently, this is imposed as a function of T, P or composition
*
* units are kmol m-3
*/
doublereal m_molar_density;
doublereal m_site_density;
// doublereal m_molar_lattice_volume;
private:

467
Cantera/src/thermo/LatticeSolidPhase.cpp
View file

@ -23,8 +23,16 @@
#include "LatticePhase.h"
#include "SpeciesThermo.h"
#include "ThermoFactory.h"
#include "SpeciesThermoFactory.h"
#include "GeneralSpeciesThermo.h"
#include <string>
#ifndef MIN
# define MIN(x,y) (( (x) < (y) ) ? (x) : (y))
#endif
#ifndef MAX
# define MAX(x,y) (( (x) > (y) ) ? (x) : (y))
#endif
using namespace std;
//======================================================================================================================
@ -39,7 +47,9 @@ namespace Cantera {
m_molar_density(0.0),
m_nlattice(0),
m_lattice(0),
m_x(0)
m_x(0),
theta_(0),
tmpV_(0)
{
}
//====================================================================================================================
@ -54,7 +64,9 @@ namespace Cantera {
m_molar_density(0.0),
m_nlattice(0),
m_lattice(0),
m_x(0)
m_x(0),
theta_(0),
tmpV_(0)
{
*this = operator=(right);
}
@ -74,12 +86,15 @@ namespace Cantera {
m_nlattice = right.m_nlattice;
deepStdVectorPointerCopy<LatticePhase>(right.m_lattice, m_lattice);
m_x = right.m_x;
theta_ = right.theta_;
tmpV_ = right.tmpV_;
}
return *this;
}
//====================================================================================================================
// Destructor
LatticeSolidPhase::~LatticeSolidPhase() {
// We own the sublattices. So we have to delete the sublattices
for (int n = 0; n < m_nlattice; n++) {
delete m_lattice[n];
m_lattice[n] = 0;
@ -98,59 +113,121 @@ namespace Cantera {
LatticeSolidPhase *igp = new LatticeSolidPhase(*this);
return (ThermoPhase *) igp;
}
//====================================================================================================================
// Minimum temperature for which the thermodynamic data for the species
// or phase are valid.
/*
* If no argument is supplied, the
* value returned will be the lowest temperature at which the
* data for \e all species are valid. Otherwise, the value
* will be only for species \a k. This function is a wrapper
* that calls the species thermo minTemp function.
*
* @param k index of the species. Default is -1, which will return the max of the min value
* over all species.
*/
doublereal LatticeSolidPhase::minTemp(int k) const {
if (k >= 0) {
for (int n = 0; n < m_nlattice; n++) {
if (lkstart_[n+1] < k) {
double ml = (m_lattice[n])->minTemp(k-lkstart_[n]);
return ml;
}
}
}
doublereal mm = 1.0E300;
for (int n = 0; n < m_nlattice; n++) {
double ml = (m_lattice[n])->minTemp(-1);
mm = MIN(mm, ml);
}
return mm;
}
//====================================================================================================================
// Maximum temperature for which the thermodynamic data for the species
// or phase are valid.
/*
* If no argument is supplied, the
* value returned will be the lowest temperature at which the
* data for \e all species are valid. Otherwise, the value
* will be only for species \a k. This function is a wrapper
* that calls the species thermo minTemp function.
*
* @param k index of the species. Default is -1, which will return the max of the min value
* over all species.
*/
doublereal LatticeSolidPhase::maxTemp(int k) const {
if (k >= 0) {
for (int n = 0; n < m_nlattice; n++) {
if (lkstart_[n+1] < k) {
double ml = (m_lattice[n])->maxTemp(k - lkstart_[n]);
return ml;
}
}
}
doublereal mm = -1.0E300;
for (int n = 0; n < m_nlattice; n++) {
double ml = (m_lattice[n])->maxTemp(-1);
mm = MAX(mm, ml);
}
return mm;
}
//====================================================================================================================
/*
* Returns the reference pressure in Pa. This function is a wrapper
* that calls the species thermo refPressure function.
*/
doublereal LatticeSolidPhase::refPressure() const {
return m_lattice[0]->refPressure();
}
//====================================================================================================================
doublereal LatticeSolidPhase::
enthalpy_mole() const {
_updateThermo();
doublereal ndens, sum = 0.0;
doublereal sum = 0.0;
int n;
for (n = 0; n < m_nlattice; n++) {
ndens = m_lattice[n]->molarDensity();
sum += ndens * m_lattice[n]->enthalpy_mole();
sum += theta_[n] * m_lattice[n]->enthalpy_mole();
}
return sum/molarDensity();
return sum;
}
//====================================================================================================================
doublereal LatticeSolidPhase::intEnergy_mole() const {
_updateThermo();
doublereal ndens, sum = 0.0;
doublereal sum = 0.0;
int n;
for (n = 0; n < m_nlattice; n++) {
ndens = m_lattice[n]->molarDensity();
sum += ndens * m_lattice[n]->intEnergy_mole();
sum += theta_[n] * m_lattice[n]->intEnergy_mole();
}
return sum/molarDensity();
return sum;
}
//====================================================================================================================
doublereal LatticeSolidPhase::entropy_mole() const {
_updateThermo();
doublereal ndens, sum = 0.0;
doublereal sum = 0.0;
int n;
for (n = 0; n < m_nlattice; n++) {
ndens = m_lattice[n]->molarDensity();
sum += ndens * m_lattice[n]->entropy_mole();
sum += theta_[n] * m_lattice[n]->entropy_mole();
}
return sum/molarDensity();
return sum;
}
//====================================================================================================================
doublereal LatticeSolidPhase::gibbs_mole() const {
_updateThermo();
doublereal ndens, sum = 0.0;
doublereal sum = 0.0;
for (int n = 0; n < m_nlattice; n++) {
ndens = m_lattice[n]->molarDensity();
sum += ndens * m_lattice[n]->gibbs_mole();
sum += theta_[n] * m_lattice[n]->gibbs_mole();
}
return sum/molarDensity();
return sum;
}
//====================================================================================================================
doublereal LatticeSolidPhase::cp_mole() const {
_updateThermo();
doublereal sum = 0.0;
for (int n = 0; n < m_nlattice; n++) {
doublereal ndens = m_lattice[n]->molarDensity();
sum += ndens * m_lattice[n]->cp_mole();
sum += theta_[n] * m_lattice[n]->cp_mole();
}
return sum/molarDensity();
return sum;
}
//====================================================================================================================
void LatticeSolidPhase::getActivityConcentrations(doublereal* c) const {
@ -176,12 +253,49 @@ namespace Cantera {
return 0.0;
}
//====================================================================================================================
// Set the pressure at constant temperature. Units: Pa.
/*
*
* @param p Pressure (units - Pa)
*/
void LatticeSolidPhase::setPressure(doublereal p) {
m_press = p;
for (int n = 0; n < m_nlattice; n++) {
m_lattice[n]->setPressure(m_press);
}
calcDensity();
}
//====================================================================================================================
// Calculate the density of the solid mixture
/*
* The formula for this is
*
* \f[
* \rho = \sum_n{ \rho_n \theta_n }
* \f]
*
* where \f$ \rho_n \f$ is the density of the nth sublattice
*
* Note this is a nonvirtual function.
*/
doublereal LatticeSolidPhase::calcDensity() {
double sum = 0.0;
for (int n = 0; n < m_nlattice; n++) {
sum += theta_[n] * m_lattice[n]->density();
}
State::setDensity(sum);
return sum;
}
//====================================================================================================================
// Set the mole fractions to the specified values, and then
// normalize them so that they sum to 1.0 for each of the subphases
/*
* On input, the mole fraction vector is assumed to sum to one for each of the sublattices. The sublattices
* are updated with this mole fraction vector. The mole fractions are also storred within this object, after
* they are normalized to one by dividing by the number of sublattices.
*
* @param x Input vector of mole fractions. There is no restriction
* @param x Input vector of mole fractions. There is no restriction
* on the sum of the mole fraction vector. Internally,
* this object will pass portions of this vector to the sublattices which assume that the portions
* individually sum to one.
@ -189,19 +303,16 @@ namespace Cantera {
*/
void LatticeSolidPhase::setMoleFractions(const doublereal* const x) {
int nsp, strt = 0;
doublereal sum = 0.0;
for (int n = 0; n < m_nlattice; n++) {
nsp = m_lattice[n]->nSpecies();
m_lattice[n]->setMoleFractions(x+strt);
for (int k = 0; k < nsp; k++) {
sum += x[strt + k];
}
m_lattice[n]->setMoleFractions(x + strt);
strt += nsp;
}
for (int k = 0; k < strt; k++) {
m_x[k] = x[k] / sum;
m_x[k] = x[k] / m_nlattice;
}
State::setMoleFractions(DATA_PTR(m_x));
calcDensity();
}
//====================================================================================================================
// Get the species mole fraction vector.
@ -213,6 +324,7 @@ namespace Cantera {
*/
void LatticeSolidPhase::getMoleFractions(doublereal* const x) const {
int nsp, strt = 0;
// the ifdef block should be the way we calculate this.!!!!!
State::getMoleFractions(x);
doublereal sum;
for (int n = 0; n < m_nlattice; n++) {
@ -241,28 +353,103 @@ namespace Cantera {
}
}
//====================================================================================================================
// Get the species chemical potentials. Units: J/kmol.
/*
* This function returns a vector of chemical potentials of the
* species in solution at the current temperature, pressure
* and mole fraction of the solution.
*
* This returns the underlying lattice chemical potentials
*
* @param mu Output vector of species chemical
* potentials. Length: m_kk. Units: J/kmol
*/
void LatticeSolidPhase::getChemPotentials(doublereal* mu) const {
_updateThermo();
int strt = 0;
for (int n = 0; n < m_nlattice; n++) {
doublereal dratio = m_lattice[n]->molarDensity()/molarDensity();
int nlsp = m_lattice[n]->nSpecies();
m_lattice[n]->getChemPotentials(mu+strt);
scale(mu + strt, mu + strt + m_lattice[n]->nSpecies(), mu + strt, dratio);
strt += m_lattice[n]->nSpecies();
strt += nlsp;
}
}
//====================================================================================================================
void LatticeSolidPhase::getPartialMolarEnthalpies(doublereal* hbar) const {
_updateThermo();
int strt = 0;
for (int n = 0; n < m_nlattice; n++) {
int nlsp = m_lattice[n]->nSpecies();
m_lattice[n]->getPartialMolarEnthalpies(hbar + strt);
strt += nlsp;
}
}
//====================================================================================================================
void LatticeSolidPhase::getPartialMolarEntropies(doublereal* sbar) const {
_updateThermo();
int strt = 0;
for (int n = 0; n < m_nlattice; n++) {
int nlsp = m_lattice[n]->nSpecies();
m_lattice[n]->getPartialMolarEntropies(sbar + strt);
strt += nlsp;
}
}
//====================================================================================================================
void LatticeSolidPhase::getPartialMolarCp(doublereal* cpbar) const {
_updateThermo();
int strt = 0;
for (int n = 0; n < m_nlattice; n++) {
int nlsp = m_lattice[n]->nSpecies();
m_lattice[n]->getPartialMolarCp(cpbar + strt);
strt += nlsp;
}
}
//====================================================================================================================
void LatticeSolidPhase::getPartialMolarVolumes(doublereal* vbar) const {
_updateThermo();
int strt = 0;
for (int n = 0; n < m_nlattice; n++) {
int nlsp = m_lattice[n]->nSpecies();
m_lattice[n]->getPartialMolarVolumes(vbar + strt);
strt += nlsp;
}
}
//====================================================================================================================
// Get the array of standard state chemical potentials at unit activity for the species
// at their standard states at the current <I>T</I> and <I>P</I> of the solution.
/*
* These are the standard state chemical potentials \f$ \mu^0_k(T,P)
* \f$. The values are evaluated at the current
* temperature and pressure of the solution.
*
* This returns the underlying lattice standard chemical potentials, as the units are kmol-1 of
* the sublattice species.
*
* @param mu0 Output vector of chemical potentials.
* Length: m_kk. Units: J/kmol
*/
void LatticeSolidPhase::getStandardChemPotentials(doublereal* mu0) const {
_updateThermo();
int strt = 0;
for (int n = 0; n < m_nlattice; n++) {
doublereal dratio = m_lattice[n]->molarDensity()/molarDensity();
m_lattice[n]->getStandardChemPotentials(mu0+strt);
scale(mu0 + strt, mu0 + strt + m_lattice[n]->nSpecies(), mu0 + strt, dratio);
strt += m_lattice[n]->nSpecies();
}
}
//====================================================================================================================
void LatticeSolidPhase::getGibbs_RT_ref(doublereal *grt) const {
_updateThermo();
for (int n = 0; n < m_nlattice; n++) {
m_lattice[n]->getGibbs_RT_ref(grt + lkstart_[n]);
}
}
//====================================================================================================================
void LatticeSolidPhase::getGibbs_ref(doublereal *g) const {
getGibbs_RT_ref(g);
for (int k = 0; k < m_kk; k++) {
g[k] *= GasConstant * temperature();
}
}
//====================================================================================================================
// Add in species from Slave phases
/*
* This hook is used for cSS_CONVENTION_SLAVE phases
@ -272,84 +459,127 @@ namespace Cantera {
void LatticeSolidPhase::installSlavePhases(Cantera::XML_Node* phaseNode)
{
int m, k;
int kk = 0;
int kstart = 0;
SpeciesThermoFactory* spFactory = SpeciesThermoFactory::factory();
SpeciesThermo * spthermo_ptr = new GeneralSpeciesThermo();
setSpeciesThermo(spthermo_ptr);
m_speciesData.clear();
XML_Node& eosdata = phaseNode->child("thermo");
XML_Node& la = eosdata.child("LatticeArray");
std::vector<XML_Node*> lattices;
la.getChildren("phase",lattices);
for (int n = 0; n < m_nlattice; n++) {
LatticePhase *lp = m_lattice[n];
XML_Node* phaseNode_ptr = lattices[n];
int nsp = lp->nSpecies();
vector<doublereal> constArr(lp->nElements());
const vector_fp& aws = lp->atomicWeights();
for (int es = 0; es < lp->nElements(); es++) {
string esName = lp->elementName(es);
double wt = aws[es];
int an = lp->atomicNumber(es);
int e298 = lp->entropyElement298(es);
int et = lp->elementType(es);
addUniqueElementAfterFreeze(esName, wt, an, e298, et);
}
const std::vector<const XML_Node *> & spNode = lp->speciesData();
kstart = kk;
for (k = 0; k < nsp; k++) {
std::string sname = lp->speciesName(k);
std::map<std::string, double> comp;
lp->getAtoms(k, DATA_PTR(constArr));
int nel = nElements();
vector_fp ecomp(nel, 0.0);
for (m = 0; m < lp->nElements(); m++) {
if (constArr[m] != 0.0) {
std::string ename = lp->elementName(m);
comp[ename] = constArr[m];
}
}
int nel = nElements();
vector_fp ecomp(nel, 0.0);
for (m = 0; m < nel; m++) {
double anum = comp[elementName(m)];
if (anum != 0.0) {
ecomp[m] = anum;
std::string oldEname = lp->elementName(m);
int newIndex = elementIndex(oldEname);
if (newIndex < 0) {
throw CanteraError("LatticeSolidPhase::installSlavePhases", "confused");
}
ecomp[newIndex] = constArr[m];
}
}
double chrg = lp->charge(k);
double sz = lp->size(k);
addUniqueSpecies(sname, &ecomp[0], chrg, sz);
spFactory->installThermoForSpecies(kk, *(spNode[k]), this, *m_spthermo, phaseNode_ptr);
m_speciesData.push_back(new XML_Node(*(spNode[k])));
kk++;
}
/*
* Add in the lattice stoichiometry constraint
*/
if (n > 0) {
string econ = "LC_";
econ += int2str(n);
econ += "_" + id();
int m = addUniqueElementAfterFreeze(econ, 0.0, 0, 0.0, CT_ELEM_TYPE_LATTICERATIO);
m_mm = nElements();
LatticePhase *lp0 = m_lattice[0];
int nsp0 = lp0->nSpecies();
for (k = 0; k < nsp0; k++) {
m_speciesComp[k * m_mm + m] = -theta_[0];
}
for (k = 0; k < nsp; k++) {
int ks = kstart + k;
m_speciesComp[ks * m_mm + m] = theta_[n];
}
}
}
}
//====================================================================================================================
// Initialize the ThermoPhase object after all species have been set up
/*
* @internal Initialize.
*
* This method is provided to allow subclasses to perform any initialization required after all
* species have been added. For example, it might be used to
* resize internal work arrays that must have an entry for
* each species. The base class implementation does nothing,
* and subclasses that do not require initialization do not
* need to overload this method. When importing a CTML phase
* description, this method is called from ThermoPhase::initThermoXML(),
* which is called from importPhase(), just prior to returning from function importPhase().
*
* @see importCTML.cpp
*/
void LatticeSolidPhase::initThermo() {
m_kk = nSpecies();
m_mm = nElements();
m_x.resize(m_kk);
initLengths();
int nsp, k, loc = 0;
doublereal ndens;
m_molar_density = 0.0;
for (int n = 0; n < m_nlattice; n++) {
nsp = m_lattice[n]->nSpecies();
ndens = m_lattice[n]->molarDensity();
lkstart_[n] = loc;
nspLattice_[n] = nsp;
for (k = 0; k < nsp; k++) {
m_x[loc] = ndens * m_lattice[n]->moleFraction(k);
m_x[loc] =m_lattice[n]->moleFraction(k) / (double) m_nlattice;
loc++;
}
m_molar_density += ndens;
lkstart_[n+1] = loc;
}
setMoleFractions(DATA_PTR(m_x));
// const vector<string>& spnames = speciesNames();
// int n, k, kl, namesize;
// int nl = m_sitedens.size();
// string s;
// m_lattice.resize(m_kk,-1);
// vector_fp conc(m_kk, 0.0);
// compositionMap xx;
// for (n = 0; n < nl; n++) {
// for (k = 0; k < m_kk; k++) {
// xx[speciesName(k)] = -1.0;
// }
// parseCompString(m_sp[n], xx);
// for (k = 0; k < m_kk; k++) {
// if (xx[speciesName(k)] != -1.0) {
// conc[k] = m_sitedens[n]*xx[speciesName(k)];
// m_lattice[k] = n;
// }
// }
// }
// for (k = 0; k < m_kk; k++) {
// if (m_lattice[k] == -1) {
// throw CanteraError("LatticeSolidPhase::"
// "setParametersFromXML","Species "+speciesName(k)
// +" not a member of any lattice.");
// }
// }
// setMoleFractions(DATA_PTR(conc));
ThermoPhase::initThermo();
}
//====================================================================================================================
// Initialize vectors that depend on the number of species and sublattices
/*
*
*/
void LatticeSolidPhase::initLengths() {
theta_.resize(m_nlattice,0);
nspLattice_.resize(m_nlattice);
lkstart_.resize(m_nlattice+1);
m_x.resize(m_kk, 0.0);
tmpV_.resize(m_kk, 0.0);
}
//====================================================================================================================
void LatticeSolidPhase::_updateThermo() const {
doublereal tnow = temperature();
@ -386,6 +616,14 @@ namespace Cantera {
setMoleFractions(DATA_PTR(m_x));
}
//====================================================================================================================
//====================================================================================================================
// Set the parameters from the XML file
/*!
* Currently, this is the spot that we read in all of the sublattice phases.
* The SetParametersFromXML() call is carried out at
*/
void LatticeSolidPhase::setParametersFromXML(const XML_Node& eosdata) {
eosdata._require("model","LatticeSolid");
XML_Node& la = eosdata.child("LatticeArray");
@ -398,9 +636,78 @@ namespace Cantera {
XML_Node& i = *lattices[n];
m_lattice.push_back((LatticePhase*)newPhase(i));
}
std::vector<string> pnam;
std::vector<string> pval;
XML_Node& ls = eosdata.child("LatticeStoichiometry");
int np = getPairs(ls, pnam, pval);
theta_.resize(nl);
for (int i = 0; i < np; i++) {
double val = fpValueCheck(pval[i]);
bool found = false;
for (int j = 0; j < nl; j++) {
ThermoPhase &tp = *(m_lattice[j]);
string idj = tp.id();
if (idj == pnam[i]) {
theta_[j] = val;
found = true;
break;
}
}
if (!found) {
throw CanteraError("", "not found");
}
}
}
//====================================================================================================================
// Return a changeable reference to the calculation manager
// for species reference-state thermodynamic properties
/*
*
* @param k Speices id. The default is -1, meaning return the default
*
* @internal
*/
SpeciesThermo& LatticeSolidPhase::speciesThermo(int k) {
return *m_spthermo;
/*
int kk;
if (k >= 0) {
for (int n = 0; n < m_nlattice; n++) {
if (lkstart_[n+1] < k) {
kk = k - lkstart_[n];
return m_lattice[n]->speciesThermo(kk);
}
}
}
return m_lattice[0]->speciesThermo(-1);
*/
}
//====================================================================================================================
#ifdef H298MODIFY_CAPABILITY
//! Modify the value of the 298 K Heat of Formation of one species in the phase (J kmol-1)
/*!
* The 298K heat of formation is defined as the enthalpy change to create the standard state
* of the species from its constituent elements in their standard states at 298 K and 1 bar.
*
* @param k Species k
* @param Hf298New Specify the new value of the Heat of Formation at 298K and 1 bar
*/
void LatticeSolidPhase::modifyOneHf298SS(const int k, const doublereal Hf298New) {
for (int n = 0; n < m_nlattice; n++) {
if (lkstart_[n+1] < k) {
int kk = k-lkstart_[n];
SpeciesThermo& l_spthermo = m_lattice[n]->speciesThermo();
l_spthermo.modifyOneHf298(kk, Hf298New);
}
}
m_tlast += 0.0001234;
_updateThermo();
}
#endif
//====================================================================================================================
doublereal LatticeSolidPhase::err(std::string msg) const {
throw CanteraError("LatticeSolidPhase","Unimplemented " + msg);

299
Cantera/src/thermo/LatticeSolidPhase.h
View file

@ -34,17 +34,49 @@
namespace Cantera {
//! A phase that is comprised of an additive combination of other lattice phases
//! A phase that is comprised of a fixed additive combination of other lattice phases
/*!
* This is the main way Cantera describes semiconductors and other solid phases.
* This %ThermoPhase object calculates its properties as a sum over other LatticePhase objects. Each of the %LatticePhase
* objects is a ThermoPhase object by itself.
* This is the main way %Cantera describes semiconductors and other solid phases.
* This %ThermoPhase object calculates its properties as a sum over other %LatticePhase objects. Each of the %LatticePhase
* objects is a %ThermoPhase object by itself.
*
* The results from this LatticeSolidPhase model reduces to the LatticePhase model when there is one
* lattice phase and the molar densities of the sublattice and the molar density within the LatticeSolidPhase
* have the same values.
*
*
* The mole fraction vector is redefined witin the the LatticeSolidPhase object. Each of the mole
* fractions sum to one on each of the sublattices. The routine getMoleFraction() and setMoleFraction()
* have been redefined to use this convention.
*
* <HR>
* <H2> Specification of Species Standard %State Properties </H2>
* <HR>
*
* The standard state properties are calculated in the normal way for each of the sublattices. The normal way
* here means that a thermodynamic polynomial in temperature is developed. Also, a constant volume approximation
* for the pressure dependence is assumed. All of these properties are on a Joules per kmol of sublattice
* constituent basis.
*
* <HR>
* <H2> Specification of Solution Thermodynamic Properties </H2>
* <HR>
* The sum over the %LatticePhase objects is carried out by weighting each %LatticePhase object
* value with the molarDensity of the LatticePhase. Then the resulting quantity is divided by
* value with the molar density (kmol m-3) of its %LatticePhase. Then the resulting quantity is divided by
* the molar density of the total compound. The LatticeSolidPhase object therefore only contains a
* listing of the number of Lattice Phases
* that comprises the solid and it contains a value for the molar density of the entire mixture.
* listing of the number of %LatticePhase object
* that comprises the solid, and it contains a value for the molar density of the entire mixture.
* This is the same thing as saying that
*
* \f[
* L_i = L^{solid} \theta_i
* \f]
*
* \f$ L_i \f$ is the molar volume of the ith lattice. \f$ L^{solid} \f$ is the molar volume of the entire
* solid. \f$ \theta_i \f$ is a fixed weighting factor for the ith lattice representing the lattice
* stoichiometric coefficient. For this object the \f$ \theta_i \f$ values are fixed.
*
*
* Let's take FeS2 as an example, which may be thought of as a combination of two lattices: Fe and S lattice.
* The Fe sublattice has a molar density of 1 gmol cm-3. The S sublattice has a molar density of 2 gmol cm-3.
@ -53,14 +85,34 @@ namespace Cantera {
* associated with the sublattices. The Fe sublattice will have a weight of 1.0 associated with it. The
* S sublattice will have a weight of 2.0 associated with it.
*
* Currently, the molar density is set to a constant.
*
* <HR>
* <H3> Specification of Solution Density Properties </H3>
* <HR>
*
* The results from this LatticeSolidPhase model reduces to the LatticePhase model when there is one
* lattice phase and the molar densities of the sublattice and the molar density within the LatticeSolidPhase
* have the same values.
* Currently, molar density is not a constant within the object, even though the species molar volumes are a
* constant. The basic idea is that a swelling of one of the sublattices will result in a swelling of
* of all of the lattices. Therefore, the molar volumes of the individual lattices are not independent of
* one another.
*
* The mole fraction vector has been redefined within the LatticeSolidPhase object. The mole fractions sum
* to one within each of the individual lattice phases. The routine getMoleFraction() and setMoleFraction()
* The molar volume of the Lattice solid is calculated from the following formula
*
* \f[
* V = \sum_i{ \theta_i V_i^{lattice}}
* \f]
*
* where \f$ V_i^{lattice} \f$ is the molar volume of the ith sublattice. This is calculated from the
* following standard formula.
*
*
* \f[
* V_i = \sum_k{ \X_k V_k}
* \f]
*
* where k is a species in the ith sublattice.
*
* The mole fraction vector is redefined witin the the LatticeSolidPhase object. Each of the mole
* fractions sum to one on each of the sublattices. The routine getMoleFraction() and setMoleFraction()
* have been redefined to use this convention.
*
* (This object is still under construction)
@ -104,6 +156,39 @@ namespace Cantera {
*/
virtual int eosType() const { return cLatticeSolid; }
//! Minimum temperature for which the thermodynamic data for the species
//! or phase are valid.
/*!
* If no argument is supplied, the
* value returned will be the lowest temperature at which the
* data for \e all species are valid. Otherwise, the value
* will be only for species \a k. This function is a wrapper
* that calls the species thermo minTemp function.
*
* @param k index of the species. Default is -1, which will return the max of the min value
* over all species.
*/
virtual doublereal minTemp(int k = -1) const;
//! Maximum temperature for which the thermodynamic data for the species
//! are valid.
/*!
* If no argument is supplied, the
* value returned will be the highest temperature at which the
* data for \e all species are valid. Otherwise, the value
* will be only for species \a k. This function is a wrapper
* that calls the species thermo maxTemp function.
*
* @param k index of the species. Default is -1, which will return the min of the max value
* over all species.
*/
virtual doublereal maxTemp(int k = -1) const;
//! Returns the reference pressure in Pa. This function is a wrapper
//! that calls the species thermo refPressure function.
virtual doublereal refPressure() const ;
//! This method returns the convention used in specification
//! of the standard state, of which there are currently two,
//! temperature based, and variable pressure based.
@ -116,12 +201,11 @@ namespace Cantera {
//! Return the Molar Enthalpy. Units: J/kmol.
/*!
* The molar enthalpy is determined by the following formula, where \f$ C_n \f$ is the
* lattice molar density of the nth lattice, and \f$ C_T \f$ is the molar density
* of the solid compound.
* The molar enthalpy is determined by the following formula, where \f$ \theta_n \f$ is the
* lattice stoichiometric coefficient of the nth lattice
*
* \f[
* \tilde h(T,P) = \frac{\sum_n C_n \tilde h_n(T,P) }{C_T},
* \tilde h(T,P) = {\sum_n \theta_n \tilde h_n(T,P) }
* \f]
*
* \f$ \tilde h_n(T,P) \f$ is the enthalpy of the n<SUP>th</SUP> lattice.
@ -133,12 +217,11 @@ namespace Cantera {
//! Return the Molar Internal Energy. Units: J/kmol.
/*!
* The molar internal energy is determined by the following formula, where \f$ C_n \f$ is the
* lattice molar density of the nth lattice, and \f$ C_T \f$ is the molar density
* of the solid compound.
* The molar enthalpy is determined by the following formula, where \f$ \theta_n \f$ is the
* lattice stoichiometric coefficient of the nth lattice
*
* \f[
* \tilde u(T,P) = \frac{\sum_n C_n \tilde u_n(T,P) }{C_T},
* \tilde u(T,P) = {\sum_n \theta_n \tilde u_n(T,P) }
* \f]
*
* \f$ \tilde u_n(T,P) \f$ is the internal energy of the n<SUP>th</SUP> lattice.
@ -148,13 +231,12 @@ namespace Cantera {
virtual doublereal intEnergy_mole() const;
//! Return the Molar Entropy. Units: J/kmol/K.
/*!
* The molar entropy is determined by the following formula, where \f$ C_n \f$ is the
* lattice molar density of the nth lattice, and \f$ C_T \f$ is the molar density
* of the solid compound.
/*!
* The molar enthalpy is determined by the following formula, where \f$ \theta_n \f$ is the
* lattice stoichiometric coefficient of the nth lattice
*
* \f[
* \tilde s(T,P) = \frac{\sum_n C_n \tilde s_n(T,P) }{C_T},
* \tilde s(T,P) = \sum_n \theta_n \tilde s_n(T,P)
* \f]
*
* \f$ \tilde s_n(T,P) \f$ is the molar entropy of the n<SUP>th</SUP> lattice.
@ -163,14 +245,13 @@ namespace Cantera {
*/
virtual doublereal entropy_mole() const;
//! Return the Molar Enthalpy. Units: J/kmol.
//! Return the Molar Gibbs energy. Units: J/kmol.
/*!
* The molar enthalpy is determined by the following formula, where \f$ C_n \f$ is the
* lattice molar density of the nth lattice, and \f$ C_T \f$ is the molar density
* of the solid compound.
* The molar gibbs free energy is determined by the following formula, where \f$ \theta_n \f$ is the
* lattice stoichiometric coefficient of the nth lattice
*
* \f[
* \tilde h(T,P) = \frac{\sum_n C_n \tilde h_n(T,P) }{C_T},
* \tilde h(T,P) = {\sum_n \theta_n \tilde h_n(T,P) }
* \f]
*
* \f$ \tilde h_n(T,P) \f$ is the enthalpy of the n<SUP>th</SUP> lattice.
@ -226,16 +307,30 @@ namespace Cantera {
*
* @param p Pressure (units - Pa)
*/
virtual void setPressure(doublereal p) {
m_press = p;
setMolarDensity(m_molar_density);
}
virtual void setPressure(doublereal p);
//! Calculate the density of the solid mixture
/*!
* The formula for this is
*
* \f[
* \rho = \sum_n{ \rho_n \theta_n }
* \f]
*
* where \f$ \rho_n \f$ is the density of the nth sublattice
*
* Note this is a nonvirtual function.
*/
doublereal calcDensity();
//! Set the mole fractions to the specified values, and then
//! normalize them so that they sum to 1.0 for each of the subphases
/*!
/*
* On input, the mole fraction vector is assumed to sum to one for each of the sublattices. The sublattices
* are updated with this mole fraction vector. The mole fractions are also storred within this object, after
* they are normalized to one by dividing by the number of sublattices.
*
* @param x Input vector of mole fractions. There is no restriction
* @param x Input vector of mole fractions. There is no restriction
* on the sum of the mole fraction vector. Internally,
* this object will pass portions of this vector to the sublattices which assume that the portions
* individually sum to one.
@ -358,20 +453,88 @@ namespace Cantera {
* species in solution at the current temperature, pressure
* and mole fraction of the solution.
*
* This returns the underlying lattice chemical potentials, as the units are kmol-1 of
* the sublattice species.
*
* @param mu Output vector of species chemical
* potentials. Length: m_kk. Units: J/kmol
*/
virtual void getChemPotentials(doublereal* mu) const;
//! Returns an array of partial molar enthalpies for the species in the mixture.
/*!
* Units (J/kmol)
* For this phase, the partial molar enthalpies are equal to the
* pure species enthalpies
* \f[
* \bar h_k(T,P) = \hat h^{ref}_k(T) + (P - P_{ref}) \hat V^0_k
* \f]
* The reference-state pure-species enthalpies, \f$ \hat h^{ref}_k(T) \f$,
* at the reference pressure,\f$ P_{ref} \f$,
* are computed by the species thermodynamic
* property manager. They are polynomial functions of temperature.
* @see SpeciesThermo
*
* @param hbar Output vector containing partial molar enthalpies.
* Length: m_kk.
*/
virtual void getPartialMolarEnthalpies(doublereal* hbar) const;
/**
* Returns an array of partial molar entropies of the species in the
* solution. Units: J/kmol/K.
* For this phase, the partial molar entropies are equal to the
* pure species entropies plus the ideal solution contribution.
* \f[
* \bar s_k(T,P) = \hat s^0_k(T) - R log(X_k)
* \f]
* The reference-state pure-species entropies,\f$ \hat s^{ref}_k(T) \f$,
* at the reference pressure, \f$ P_{ref} \f$, are computed by the
* species thermodynamic
* property manager. They are polynomial functions of temperature.
* @see SpeciesThermo
*
* @param sbar Output vector containing partial molar entropies.
* Length: m_kk.
*/
virtual void getPartialMolarEntropies(doublereal* sbar) const;
/**
* Returns an array of partial molar Heat Capacities at constant
* pressure of the species in the
* solution. Units: J/kmol/K.
* For this phase, the partial molar heat capacities are equal
* to the standard state heat capacities.
*
* @param cpbar Output vector of partial heat capacities. Length: m_kk.
*/
virtual void getPartialMolarCp(doublereal* cpbar) const;
/**
* returns an array of partial molar volumes of the species
* in the solution. Units: m^3 kmol-1.
*
* For this solution, thepartial molar volumes are equal to the
* constant species molar volumes.
*
* @param vbar Output vector of partial molar volumes. Length: m_kk.
*/
virtual void getPartialMolarVolumes(doublereal* vbar) const;
//! Get the array of standard state chemical potentials at unit activity for the species
//! at their standard states at the current <I>T</I> and <I>P</I> of the solution.
/*!
* These are the standard state chemical potentials \f$ \mu^0_k(T,P)
* \f$. The values are evaluated at the current
* temperature and pressure of the solution
* temperature and pressure of the solution.
*
*
* This returns the underlying lattice standard chemical potentials, as the units are kmol-1 of
* the sublattice species.
*
* @param mu0 Output vector of chemical potentials.
* Length: m_kk.
* Length: m_kk. Units: J/kmol
*/
virtual void getStandardChemPotentials(doublereal* mu0) const;
@ -400,6 +563,36 @@ namespace Cantera {
* @param k index of the species (defaults to zero)
*/
virtual doublereal logStandardConc(int k=0) const;
//@}
/// @name Thermodynamic Values for the Species Reference States --------------------
//@{
/**
* Returns the vector of nondimensional
* enthalpies of the reference state at the current temperature
* of the solution and the reference pressure for the species.
*
* This function fills in its one entry in hrt[] by calling
* the underlying species thermo function for the
* dimensionless gibbs free energy, calculated from the
* dimensionless enthalpy and entropy.
*/
virtual void getGibbs_RT_ref(doublereal *grt) const;
/**
* Returns the vector of the
* gibbs function of the reference state at the current temperature
* of the solution and the reference pressure for the species.
* units = J/kmol
*
* This function fills in its one entry in g[] by calling
* the underlying species thermo functions for the
* gibbs free energy, calculated from enthalpy and the
* entropy, and the multiplying by RT.
*/
virtual void getGibbs_ref(doublereal *g) const;
//! Initialize the ThermoPhase object after all species have been set up
/*!
@ -420,6 +613,9 @@ namespace Cantera {
*/
virtual void initThermo();
//! Initialize vectors that depend on the number of species and sublattices
void initLengths();
//! Add in species from Slave phases
/*!
* This hook is used for cSS_CONVENTION_SLAVE phases
@ -454,6 +650,16 @@ namespace Cantera {
*/
void setLatticeMoleFractionsByName(int n, std::string x);
//! Return a changeable reference to the calculation manager
//! for species reference-state thermodynamic properties
/*!
* This routine returns the calculation manager for the sublattice
*
* @param k Speices id. The default is -1, meaning return the default
*
* @internal
*/
virtual SpeciesThermo& speciesThermo(int k = -1);
#ifdef H298MODIFY_CAPABILITY
@ -465,10 +671,7 @@ namespace Cantera {
* @param k Species k
* @param Hf298New Specify the new value of the Heat of Formation at 298K and 1 bar
*/
virtual void modifyOneHf298SS(const int k, const doublereal Hf298New) {
m_spthermo->modifyOneHf298(k, Hf298New);
m_tlast += 0.0001234;
}
virtual void modifyOneHf298SS(const int k, const doublereal Hf298New);
#endif
private:
@ -507,6 +710,16 @@ namespace Cantera {
*/
mutable vector_fp m_x;
//! Lattice stoichiometric coefficients
std::vector<doublereal> theta_;
//! Temporary vector
mutable vector_fp tmpV_;
std::vector<doublereal> nspLattice_;
std::vector<int> lkstart_;
private:
//! Update the reference thermodynamic functions

14
Cantera/src/thermo/ThermoFactory.cpp
View file

@ -525,7 +525,9 @@ namespace Cantera {
/***************************************************************
* Add the elements.
***************************************************************/
th->addElementsFromXML(phase);
if (ssConvention != cSS_CONVENTION_SLAVE) {
th->addElementsFromXML(phase);
}
/***************************************************************
* Add the species.
@ -538,11 +540,13 @@ namespace Cantera {
vector<XML_Node*> sparrays;
phase.getChildren("speciesArray", sparrays);
int jsp, nspa = static_cast<int>(sparrays.size());
if (nspa == 0) {
throw CanteraError("importPhase",
"phase, " + th->id() + ", has zero \"speciesArray\" XML nodes.\n"
if (ssConvention != cSS_CONVENTION_SLAVE) {
if (nspa == 0) {
throw CanteraError("importPhase",
"phase, " + th->id() + ", has zero \"speciesArray\" XML nodes.\n"
+ " There must be at least one speciesArray nodes "
"with one or more species");
}
}
vector<XML_Node*> dbases;
vector_int sprule(nspa,0);
@ -608,7 +612,7 @@ namespace Cantera {
// If the phase has a species thermo manager already installed,
// delete it since we are adding new species.
delete &th->speciesThermo();
//delete &th->speciesThermo();
// Decide whether the the phase has a variable pressure ss or not
SpeciesThermo* spth = 0;

40
Cantera/src/thermo/ThermoPhase.cpp
View file

@ -822,9 +822,47 @@ namespace Cantera {
if (i == 4) uA[4] = 0.0;
if (i == 5) uA[5] = 0.0;
}
}
//=================================================================================================================
// Install a species thermodynamic property manager.
/*
* The species thermodynamic property manager
* computes properties of the pure species for use in
* constructing solution properties. It is meant for internal
* use, and some classes derived from ThermoPhase may not use
* any species thermodynamic property manager. This method is
* called by function importPhase() in importCTML.cpp.
*
* @param spthermo input pointer to the species thermodynamic property
* manager.
*
* @internal
*/
void ThermoPhase::setSpeciesThermo(SpeciesThermo* spthermo) {
if (m_spthermo) {
if (m_spthermo != spthermo) {
delete m_spthermo;
}
}
m_spthermo = spthermo;
}
//=================================================================================================================
// Return a changeable reference to the calculation manager
// for species reference-state thermodynamic properties
/*
*
* @param k Speices id. The default is -1, meaning return the default
*
* @internal
*/
SpeciesThermo& ThermoPhase::speciesThermo(int k) {
if (!m_spthermo) {
throw CanteraError("ThermoPhase::speciesThermo()",
"species reference state thermo manager was not set");
}
return *m_spthermo;
}
//=================================================================================================================
/*
* initThermoFile():
*

14
Cantera/src/thermo/ThermoPhase.h
View file

@ -758,7 +758,7 @@ namespace Cantera {
* Returns the reference pressure in Pa. This function is a wrapper
* that calls the species thermo refPressure function.
*/
doublereal refPressure() const {
virtual doublereal refPressure() const {
return m_spthermo->refPressure();
}
@ -775,7 +775,7 @@ namespace Cantera {
* @param k index of the species. Default is -1, which will return the max of the min value
* over all species.
*/
doublereal minTemp(int k = -1) const {
virtual doublereal minTemp(int k = -1) const {
return m_spthermo->minTemp(k);
}
@ -844,7 +844,7 @@ namespace Cantera {
* @param k index of the species. Default is -1, which will return the min of the max value
* over all species.
*/
doublereal maxTemp(int k = -1) const {
virtual doublereal maxTemp(int k = -1) const {
return m_spthermo->maxTemp(k);
}
@ -1883,19 +1883,17 @@ namespace Cantera {
*
* @internal
*/
void setSpeciesThermo(SpeciesThermo* spthermo)
{ m_spthermo = spthermo; }
void setSpeciesThermo(SpeciesThermo* spthermo);
//! Return a changeable reference to the calculation manager
//! for species reference-state thermodynamic properties
/*!
*
* @todo This method will fail if no species thermo
* manager has been installed.
* @param k Speices id. The default is -1, meaning return the default
*
* @internal
*/
SpeciesThermo& speciesThermo() { return *m_spthermo; }
virtual SpeciesThermo& speciesThermo(int k = -1);
/**
* @internal

2
Cantera/src/thermo/mix_defs.h
View file

@ -47,7 +47,7 @@ namespace Cantera {
const int cMetalSHEelectrons = 9; // SHE electrode electrons
const int cLatticeSolid = 20; // LatticeSolidPhase.h
const int cLattice = 21;
const int cLattice = 21; //LatticePhase.h
// pure fluids with liquid/vapor eqs of state
const int cPureFluid = 10;