added electric potential

This commit is contained in:
Dave Goodwin
2006-05-03 19:45:39 +00:00
parent bd2c97399e
commit bb19de47bc
19 changed files with 375 additions and 226 deletions

View File

@@ -433,6 +433,11 @@ extern "C" {
catch (CanteraError) {return DERR;}
}
double DLL_EXPORT th_electricPotential(int n) {
try {return th(n)->electricPotential();}
catch (CanteraError) {return DERR;}
}
int DLL_EXPORT th_chemPotentials(int n, int lenm, double* murt) {
thermo_t* thrm = th(n);
int nsp = thrm->nSpecies();

View File

@@ -68,6 +68,7 @@ extern "C" {
double DLL_IMPORT th_gibbs_mass(int n);
double DLL_IMPORT th_cp_mass(int n);
double DLL_IMPORT th_cv_mass(int n);
double DLL_IMPORT th_electricPotential(int n);
int DLL_IMPORT th_chemPotentials(int n, int lenm, double* murt);
int DLL_IMPORT th_elementPotentials(int n, int lenm, double* lambda);
int DLL_IMPORT th_getEnthalpies_RT(int n, int lenm, double* h_rt);

View File

@@ -3,7 +3,7 @@ import string
import os
from constants import *
from SurfacePhase import SurfacePhase
from SurfacePhase import SurfacePhase, EdgePhase
from Kinetics import Kinetics
import XML
@@ -48,8 +48,6 @@ class Interface(SurfacePhase, Kinetics):
if len(fn) > 1:
id = fn[1]
fn = fn[0]
#fname = os.path.basename(fn)
#ff = os.path.splitext(fname)
# read in the root element of the tree if not building from
# an already-built XML tree. Enable preprocessing if the film

View File

@@ -68,7 +68,10 @@ class ThermoPhase(Phase):
def name(self):
return self.idtag
def setName(self, name):
self.idtag = name
def refPressure(self):
"""Reference pressure [Pa].
All standard-state thermodynamic properties are for this pressure.
@@ -127,6 +130,10 @@ class ThermoPhase(Phase):
""" The pressure [Pa]."""
return _cantera.thermo_getfp(self._phase_id,7)
def electricPotential(self):
"""Electric potential [V]."""
return _cantera.thermo_getfp(self._phase_id,25)
def chemPotentials(self, species = []):
"""Species chemical potentials.

View File

@@ -516,6 +516,54 @@ class thermo:
def export(self, f, fmt = 'CSV'):
pass
class Mu0_table(thermo):
"""Properties are computed by specifying a table of standard
chemical potentials vs. T."""
def __init__(self, range = (0.0, 0.0),
h298 = 0.0,
mu0 = None,
p0 = -1.0):
self._t = range
self._h298 = h298
self._mu0 = mu0
self._pref = p0
def build(self, t):
n = t.addChild("Mu0")
n['Tmin'] = `self._t[0]`
n['Tmax'] = `self._t[1]`
if self._pref <= 0.0:
n['P0'] = `_pref`
else:
n['P0'] = `self._pref`
energy_units = _uenergy+'/'+_umol
addFloat(n,"H298", self._h298, defunits = energy_units)
n.addChild("numPoints", len(self._mu0))
mustr = ''
tstr = ''
col = 0
for v in self._mu0:
mu0 = v[1]
t = v[0]
tstr += '%17.9E, ' % t
mustr += '%17.9E, ' % mu0
col += 1
if col == 3:
tstr = tstr[:-2]+'\n'
mustr = mustr[:-2]+'\n'
col = 0
u = n.addChild("floatArray", mustr)
u["size"] = "numPoints"
u["name"] = "Mu0Values"
u = n.addChild("floatArray", tstr)
u["size"] = "numPoints"
u["name"] = "Mu0Temperatures"
class NASA(thermo):
"""NASA polynomial parameterization."""

View File

@@ -101,6 +101,8 @@ thermo_getfp(PyObject *self, PyObject *args)
vv = th_cp_mass(th); break;
case 13:
vv = th_cv_mass(th); break;
case 25:
vv = th_electricPotential(th); break;
case 50:
vv = th_critTemperature(th); break;
case 51:

View File

@@ -72,11 +72,9 @@ namespace Cantera {
/**
* Install parameterization for a species.
* @param index Species index
* @param type ignored, since only NASA type is supported
* @param c coefficients. These are
* - c[0] midpoint temperature
* - c[1] - c[7] coefficients for low T range
* - c[8] - c[14] coefficients for high T range
* @param type parameterization type
* @param c coefficients. The meaning of these depends on
* the parameterization.
*/
void GeneralSpeciesThermo::install(string name,
int index,
@@ -117,7 +115,7 @@ namespace Cantera {
break;
case SHOMATE2:
m_sp[index] = new ShomatePoly2(index, minTemp, maxTemp,
refPressure, c);
refPressure, c);
break;
case NASA2:
m_sp[index] = new NasaPoly2(index, minTemp, maxTemp,
@@ -134,7 +132,7 @@ namespace Cantera {
}
/**
* Update the properties for all species;
* Update the properties for one species.
*/
void GeneralSpeciesThermo::
update_one(int k, doublereal t, doublereal* cp_R,
@@ -145,7 +143,7 @@ namespace Cantera {
/**
* Update the properties for all species;
* Update the properties for all species.
*/
void GeneralSpeciesThermo::
update(doublereal t, doublereal* cp_R,

View File

@@ -24,8 +24,8 @@ namespace Cantera {
* A species thermodynamic property manager for a phase.
* This is a general manager that can handle a wide variety
* of species thermodynamic polynomials for individual species.
* It is slow, however.
*
* It is slow, however, because it recomputes the functions of
* temperature needed for each species.
*
*/
class GeneralSpeciesThermo : public SpeciesThermo {

View File

@@ -108,25 +108,25 @@ namespace Cantera {
*/
bool ifound = false;
for (i = 0, iindex = 2; i < nPoints; i++) {
T1 = coeffs[iindex];
m_t0_int[i] = T1;
m_mu0_R_int[i] = coeffs[iindex+1] / GasConstant;
if (T1 == 298.15) {
iT298 = i;
ifound = true;
}
if (i < nPoints - 1) {
T2 = coeffs[iindex+2];
if (T2 <= T1) {
throw CanteraError("Mu0Poly",
"Temperatures are not monotonic increasing");
}
}
iindex += 2;
T1 = coeffs[iindex];
m_t0_int[i] = T1;
m_mu0_R_int[i] = coeffs[iindex+1] / GasConstant;
if (T1 == 298.15) {
iT298 = i;
ifound = true;
}
if (i < nPoints - 1) {
T2 = coeffs[iindex+2];
if (T2 <= T1) {
throw CanteraError("Mu0Poly",
"Temperatures are not monotonic increasing");
}
}
iindex += 2;
}
if (!ifound) {
throw CanteraError("Mu0Poly",
"One temperature has to be 298.15");
throw CanteraError("Mu0Poly",
"One temperature has to be 298.15");
}
/*
@@ -138,61 +138,61 @@ namespace Cantera {
m_h0_R_int[iT298] = m_H298;
m_s0_R_int[iT298] = - (mu1 - m_h0_R_int[iT298]) / T1;
for (i = iT298; i < m_numIntervals; i++) {
T1 = m_t0_int[i];
s1 = m_s0_R_int[i];
h1 = m_h0_R_int[i];
mu1 = m_mu0_R_int[i];
T2 = m_t0_int[i+1];
mu2 = m_mu0_R_int[i+1];
deltaMu = mu2 - mu1;
deltaT = T2 - T1;
cpi = (deltaMu - T1 * s1 + T2 * s1) / (deltaT - T2 * log(T2/T1));
h2 = h1 + cpi * deltaT;
s2 = s1 + cpi * log(T2/T1);
m_cp0_R_int[i] = cpi;
m_h0_R_int[i+1] = h2;
m_s0_R_int[i+1] = s2;
m_cp0_R_int[i+1] = cpi;
T1 = m_t0_int[i];
s1 = m_s0_R_int[i];
h1 = m_h0_R_int[i];
mu1 = m_mu0_R_int[i];
T2 = m_t0_int[i+1];
mu2 = m_mu0_R_int[i+1];
deltaMu = mu2 - mu1;
deltaT = T2 - T1;
cpi = (deltaMu - T1 * s1 + T2 * s1) / (deltaT - T2 * log(T2/T1));
h2 = h1 + cpi * deltaT;
s2 = s1 + cpi * log(T2/T1);
m_cp0_R_int[i] = cpi;
m_h0_R_int[i+1] = h2;
m_s0_R_int[i+1] = s2;
m_cp0_R_int[i+1] = cpi;
}
/*
* Starting from the interval with T298, we go down
*/
if (iT298 > 0) {
T2 = m_t0_int[iT298];
mu2 = m_mu0_R_int[iT298];
m_h0_R_int[iT298] = m_H298;
m_s0_R_int[iT298] = - (mu2 - m_h0_R_int[iT298]) / T2;
for (i = iT298 - 1; i >= 0; i--) {
T1 = m_t0_int[i];
mu1 = m_mu0_R_int[i];
T2 = m_t0_int[i+1];
mu2 = m_mu0_R_int[i+1];
s2 = m_s0_R_int[i+1];
h2 = m_h0_R_int[i+1];
deltaMu = mu2 - mu1;
deltaT = T2 - T1;
cpi = (deltaMu - T1 * s2 + T2 * s2) / (deltaT - T1 * log(T2/T1));
h1 = h2 - cpi * deltaT;
s1 = s2 - cpi * log(T2/T1);
m_cp0_R_int[i] = cpi;
m_h0_R_int[i] = h1;
m_s0_R_int[i] = s1;
if (i == (m_numIntervals-1)) {
m_cp0_R_int[i+1] = cpi;
}
}
T2 = m_t0_int[iT298];
mu2 = m_mu0_R_int[iT298];
m_h0_R_int[iT298] = m_H298;
m_s0_R_int[iT298] = - (mu2 - m_h0_R_int[iT298]) / T2;
for (i = iT298 - 1; i >= 0; i--) {
T1 = m_t0_int[i];
mu1 = m_mu0_R_int[i];
T2 = m_t0_int[i+1];
mu2 = m_mu0_R_int[i+1];
s2 = m_s0_R_int[i+1];
h2 = m_h0_R_int[i+1];
deltaMu = mu2 - mu1;
deltaT = T2 - T1;
cpi = (deltaMu - T1 * s2 + T2 * s2) / (deltaT - T1 * log(T2/T1));
h1 = h2 - cpi * deltaT;
s1 = s2 - cpi * log(T2/T1);
m_cp0_R_int[i] = cpi;
m_h0_R_int[i] = h1;
m_s0_R_int[i] = s1;
if (i == (m_numIntervals-1)) {
m_cp0_R_int[i+1] = cpi;
}
}
}
#ifdef DEBUG_HKM_NOT
printf(" Temp mu0(J/kmol) cp0(J/kmol/K) "
" h0(J/kmol) s0(J/kmol/K) \n");
" h0(J/kmol) s0(J/kmol/K) \n");
for (i = 0; i < nPoints; i++) {
printf("%12.3g %12.5g %12.5g %12.5g %12.5g\n",
m_t0_int[i], m_mu0_R_int[i] * GasConstant,
m_cp0_R_int[i]* GasConstant,
m_h0_R_int[i]* GasConstant,
m_s0_R_int[i]* GasConstant);
fflush(stdout);
printf("%12.3g %12.5g %12.5g %12.5g %12.5g\n",
m_t0_int[i], m_mu0_R_int[i] * GasConstant,
m_cp0_R_int[i]* GasConstant,
m_h0_R_int[i]* GasConstant,
m_s0_R_int[i]* GasConstant);
fflush(stdout);
}
#endif
}
@@ -261,15 +261,15 @@ namespace Cantera {
*/
void Mu0Poly::
updateProperties(const doublereal* tt, doublereal* cp_R,
doublereal* h_RT, doublereal* s_R) const {
doublereal* h_RT, doublereal* s_R) const {
int j = m_numIntervals;
double T = *tt;
for (int i = 0; i < m_numIntervals; i++) {
double T2 = m_t0_int[i+1];
if (T <=T2) {
j = i;
break;
}
double T2 = m_t0_int[i+1];
if (T <=T2) {
j = i;
break;
}
}
double T1 = m_t0_int[j];
double cp_Rj = m_cp0_R_int[j];

View File

@@ -25,6 +25,9 @@ namespace Cantera {
* A species thermodynamic property manager for the NASA
* polynomial parameterization with two temperature ranges.
*
* This class is designed to efficiently evaluate the properties
* of a large number of species with the NASA parameterization.
*
* The original NASA polynomial parameterization expressed the
* heat capacity as a fourth-order polynomial in temperature, with
* separate coefficients for each of two temperature ranges. (The
@@ -71,8 +74,6 @@ namespace Cantera {
doublereal minTemp, doublereal maxTemp,
doublereal refPressure) {
//writelog("installing NASA thermo for species "+name+"\n");
//writelog(" index = "+int2str(index)+"\n");
int imid = int(c[0]); // midpoint temp converted to integer
int igrp = m_index[imid]; // has this value been seen before?
if (igrp == 0) { // if not, prepare new group
@@ -211,42 +212,40 @@ namespace Cantera {
doublereal &refPressure) {
type = reportType(index);
if (type == NASA) {
int grp = m_group_map[index];
int pos = m_posInGroup_map[index];
const vector<NasaPoly1> &mlg = m_low[grp-1];
const vector<NasaPoly1> &mhg = m_high[grp-1];
const NasaPoly1 *lowPoly = &(mlg[pos]);
const NasaPoly1 *highPoly = &(mhg[pos]);
int itype = NASA;
doublereal tmid = lowPoly->maxTemp();
int grp = m_group_map[index];
int pos = m_posInGroup_map[index];
const vector<NasaPoly1> &mlg = m_low[grp-1];
const vector<NasaPoly1> &mhg = m_high[grp-1];
const NasaPoly1 *lowPoly = &(mlg[pos]);
const NasaPoly1 *highPoly = &(mhg[pos]);
int itype = NASA;
doublereal tmid = lowPoly->maxTemp();
c[0] = tmid;
int n;
double ttemp;
lowPoly->reportParameters(n, itype, minTemp, ttemp, refPressure,
c + 1);
c + 1);
if (n != index) {
throw CanteraError(" ", "confused");
throw CanteraError(" ", "confused");
}
if (itype != NASA1) {
throw CanteraError(" ", "confused");
throw CanteraError(" ", "confused");
}
highPoly->reportParameters(n, itype, ttemp, maxTemp, refPressure,
c + 8);
c + 8);
if (n != index) {
throw CanteraError(" ", "confused");
throw CanteraError(" ", "confused");
}
if (itype != NASA1) {
throw CanteraError(" ", "confused");
throw CanteraError(" ", "confused");
}
} else {
throw CanteraError(" ", "confused");
throw CanteraError(" ", "confused");
}
}
protected:
// mutable map<int, NasaPoly1*> m_low_map;
// mutable map<int, NasaPoly1*> m_high_map;
vector<vector<NasaPoly1> > m_high;
vector<vector<NasaPoly1> > m_low;
map<int, int> m_index;

View File

@@ -23,11 +23,30 @@ namespace Cantera {
/**
* @defgroup spthermo Species Standard-State Thermodynamic Properties
*
* Species thermodynamic property managers compute the
* standard-state properties of pure species. They are designed
* for use by thermodynamic property managers (subclasses of
* ThermoPhase) to compute the thermodynamic properties of
* solutions.
* To compute the thermodynamic properties of multicomponent
* solutions, it is necessary to know something about the
* thermodynamic properties of the individual species present in
* the solution. Exactly what sort of species properties are
* required depends on the thermodynamic model for the
* solution. For a gaseous solution (i.e., a gas mixture), the
* species properties required are usually ideal gas properties at
* the mixture temperature and at a reference pressure (often 1
* atm or 1 bar). For other types of solutions, however, it may
* not be possible to isolate the species in a "pure" state. For
* example, the thermodynamic properties of, say, Na+ and Cl- in
* saltwater are not easily determined from data on the properties
* of solid NaCl, or solid Na metal, or chlorine gas. In this
* case, the solvation in water is fundamental to the identity of
* the species, and some other reference state must be used. One
* common convention for liquid solutions is to use thermodynamic
* data for the solutes for the limit of infinite dilution in the
* pure solvent; another convention is to reference all properties
* to unit molality.
*
* Whatever the conventions used by a particular solution model,
* means need to be provided to compute the species properties in
* the reference state. Class SpeciesThermo is the base class
* for a family of classes that compute these properties.
*/
@@ -50,7 +69,8 @@ namespace Cantera {
* install a new species thermodynamic property
* parameterization for one species.
* @param index The 'update' method will update the property
* values at position \i index in the property arrays.
* values for this species
* at position \i index in the property arrays.
* @param type int flag specifying the type of parameterization to be
* installed.
* @param c vector of coefficients for the parameterization.
@@ -125,7 +145,7 @@ namespace Cantera {
/**
* This utility function reports the type of parameterization
* used for the species, index.
* used for the species with index number index.
*/
virtual int reportType(int index) const = 0;

View File

@@ -32,30 +32,37 @@ namespace Cantera {
SpeciesThermoFactory* SpeciesThermoFactory::s_factory = 0;
/**
* Examine the types of species thermo parameterizations,
* and return a SpeciesThermo manager that can handle the
* parameterizations present.
*/
static void getSpeciesThermoTypes(XML_Node* node,
int& has_nasa, int& has_shomate, int& has_simple,
int &has_other) {
int& has_nasa, int& has_shomate, int& has_simple,
int &has_other) {
const XML_Node& sparray = *node;
vector<XML_Node*> sp;
// get all of the species nodes
sparray.getChildren("species",sp);
int ns = static_cast<int>(sp.size());
for (int n = 0; n < ns; n++) {
XML_Node* spNode = sp[n];
if (spNode->hasChild("thermo")) {
const XML_Node& th = sp[n]->child("thermo");
if (th.hasChild("NASA")) has_nasa = 1;
if (th.hasChild("Shomate")) has_shomate = 1;
if (th.hasChild("const_cp")) has_simple = 1;
if (th.hasChild("poly")) {
if (th.child("poly")["order"] == "1") has_simple = 1;
else throw CanteraError("newSpeciesThermo",
"poly with order > 1 not yet supported");
size_t n, ns = sp.size();
for (n = 0; n < ns; n++) {
XML_Node* spNode = sp[n];
if (spNode->hasChild("thermo")) {
const XML_Node& th = sp[n]->child("thermo");
if (th.hasChild("NASA")) has_nasa = 1;
if (th.hasChild("Shomate")) has_shomate = 1;
if (th.hasChild("const_cp")) has_simple = 1;
if (th.hasChild("poly")) {
if (th.child("poly")["order"] == "1") has_simple = 1;
else throw CanteraError("newSpeciesThermo",
"poly with order > 1 not yet supported");
}
if (th.hasChild("Mu0")) has_other = 1;
} else {
throw UnknownSpeciesThermoModel("getSpeciesThermoTypes:",
spNode->attrib("name"), "missing");
}
if (th.hasChild("Mu0")) has_other = 1;
} else {
throw UnknownSpeciesThermoModel("getSpeciesThermoTypes:",
spNode->attrib("name"), "missing");
}
}
}
@@ -67,18 +74,19 @@ namespace Cantera {
SpeciesThermo* SpeciesThermoFactory::newSpeciesThermo(XML_Node* node) {
int inasa = 0, ishomate = 0, isimple = 0, iother = 0;
try {
getSpeciesThermoTypes(node, inasa, ishomate, isimple, iother);
} catch (UnknownSpeciesThermoModel) {
iother = 1;
getSpeciesThermoTypes(node, inasa, ishomate, isimple, iother);
}
catch (UnknownSpeciesThermoModel) {
iother = 1;
}
if (iother) {
writelog("returning new GeneralSpeciesThermo");
return new GeneralSpeciesThermo();
return new GeneralSpeciesThermo();
}
return newSpeciesThermo(NASA*inasa
+ SHOMATE*ishomate + SIMPLE*isimple);
}
SpeciesThermo* SpeciesThermoFactory::
newSpeciesThermo(vector<XML_Node*> nodes) {
int n = static_cast<int>(nodes.size());
@@ -99,6 +107,9 @@ namespace Cantera {
}
/**
* @todo is this used?
*/
SpeciesThermo* SpeciesThermoFactory::
newSpeciesThermoOpt(vector<XML_Node*> nodes) {
int n = static_cast<int>(nodes.size());
@@ -120,6 +131,7 @@ namespace Cantera {
}
SpeciesThermo* SpeciesThermoFactory::newSpeciesThermo(int type) {
switch (type) {
@@ -156,8 +168,10 @@ namespace Cantera {
writelog("\n**** WARNING ****\nFor species "+name+
", discontinuity in cp/R detected at Tmid = "
+fp2str(tmid)+"\n");
writelog("\tValue computed using low-temperature polynomial: "+fp2str(cplow)+".\n");
writelog("\tValue computed using high-temperature polynomial: "+fp2str(cphigh)+".\n");
writelog("\tValue computed using low-temperature polynomial: "
+fp2str(cplow)+".\n");
writelog("\tValue computed using high-temperature polynomial: "
+fp2str(cphigh)+".\n");
}
// enthalpy
@@ -165,10 +179,13 @@ namespace Cantera {
doublereal hrthigh = enthalpy_RT(tmid, chigh);
delta = hrtlow - hrthigh;
if (fabs(delta/hrtlow) > 0.001) {
writelog("\n**** WARNING ****\nFor species "+name+", discontinuity in h/RT detected at Tmid = "
writelog("\n**** WARNING ****\nFor species "+name+
", discontinuity in h/RT detected at Tmid = "
+fp2str(tmid)+"\n");
writelog("\tValue computed using low-temperature polynomial: "+fp2str(hrtlow)+".\n");
writelog("\tValue computed using high-temperature polynomial: "+fp2str(hrthigh)+".\n");
writelog("\tValue computed using low-temperature polynomial: "
+fp2str(hrtlow)+".\n");
writelog("\tValue computed using high-temperature polynomial: "
+fp2str(hrthigh)+".\n");
}
// entropy
@@ -176,10 +193,13 @@ namespace Cantera {
doublereal srhigh = entropy_R(tmid, chigh);
delta = srlow - srhigh;
if (fabs(delta/srlow) > 0.001) {
writelog("\n**** WARNING ****\nFor species "+name+", discontinuity in s/R detected at Tmid = "
writelog("\n**** WARNING ****\nFor species "+name+
", discontinuity in s/R detected at Tmid = "
+fp2str(tmid)+"\n");
writelog("\tValue computed using low-temperature polynomial: "+fp2str(srlow)+".\n");
writelog("\tValue computed using high-temperature polynomial: "+fp2str(srhigh)+".\n");
writelog("\tValue computed using low-temperature polynomial: "
+fp2str(srlow)+".\n");
writelog("\tValue computed using high-temperature polynomial: "
+fp2str(srhigh)+".\n");
}
}
@@ -187,6 +207,8 @@ namespace Cantera {
/**
* Install a NASA polynomial thermodynamic property
* parameterization for species k into a SpeciesThermo instance.
* This is called by method installThermoForSpecies if a NASA
* block is found in the XML input.
*/
static void installNasaThermoFromXML(string speciesName,
SpeciesThermo& sp, int k,
@@ -194,8 +216,14 @@ namespace Cantera {
doublereal tmin0, tmax0, tmin1, tmax1, tmin, tmid, tmax;
const XML_Node& f0 = *f0ptr;
// default to a single temperature range
bool dualRange = false;
// but if f1ptr is suppled, then it is a two-range
// parameterization
if (f1ptr) {dualRange = true;}
tmin0 = fpValue(f0["Tmin"]);
tmax0 = fpValue(f0["Tmax"]);
tmin1 = tmax0;
@@ -207,6 +235,7 @@ namespace Cantera {
vector_fp c0, c1;
if (fabs(tmax0 - tmin1) < 0.01) {
// f0 has the lower T data, and f1 the higher T data
tmin = tmin0;
tmid = tmax0;
tmax = tmax1;
@@ -214,11 +243,13 @@ namespace Cantera {
if (dualRange)
getFloatArray(f1ptr->child("floatArray"), c1, false);
else {
// if there is no higher range data, then copy c0 to c1.
c1.resize(7,0.0);
copy(c0.begin(), c0.end(), c1.begin());
}
}
else if (fabs(tmax1 - tmin0) < 0.01) {
// f1 has the lower T data, and f0 the higher T data
tmin = tmin1;
tmid = tmax1;
tmax = tmax0;
@@ -229,6 +260,9 @@ namespace Cantera {
throw CanteraError("installNasaThermo",
"non-continuous temperature ranges.");
}
// The NasaThermo species property manager expects the
// coefficients in a different order, so rearrange them.
array_fp c(15);
c[0] = tmid;
doublereal p0 = OneAtm;

View File

@@ -44,10 +44,24 @@ namespace Cantera {
public:
/**
* 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 SpeciesThermoFactory* factory() {
if (!s_factory) s_factory = new SpeciesThermoFactory;
return s_factory;
}
/**
* 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.
*/
static void deleteFactory() {
if (s_factory) {
delete s_factory;
@@ -55,8 +69,9 @@ namespace Cantera {
}
}
/**
* Destructor doesn't do anything. We do not delete statically
* Destructor. Doesn't do anything. We do not delete statically
* created single instance of this class here, because it would
* create an infinite loop if destructor is called for that
* single instance.
@@ -69,6 +84,7 @@ namespace Cantera {
* @param type the type to be created.
*/
virtual SpeciesThermo* newSpeciesThermo(int type);
virtual SpeciesThermo* newSpeciesThermo(XML_Node* node);
virtual SpeciesThermo* newSpeciesThermo(vector<XML_Node*> nodes);
virtual SpeciesThermo* newSpeciesThermoOpt(vector<XML_Node*> nodes);
@@ -78,13 +94,28 @@ namespace Cantera {
SpeciesThermo& spthermo);
private:
/// pointer to the sole instance of this class
static SpeciesThermoFactory* s_factory;
/// Constructor. This is made private, so that only the static
/// method factory() can instantiate the class.
SpeciesThermoFactory(){}
};
////////////////////// Convenience functions ////////////////////
//
// These functions allow using a different factory class that
// derives from SpeciesThermoFactory.
//
//////////////////////////////////////////////////////////////////
/**
* Create a new species thermo manager instance.
* Create a new species thermo manager instance, by specifying
* the type and (optionally) a pointer to the factory to use to
* create it.
*/
inline SpeciesThermo* newSpeciesThermoMgr(int type,
SpeciesThermoFactory* f=0) {

View File

@@ -9,12 +9,15 @@
// Copyright 2001 California Institute of Technology
#include "speciesThermoTypes.h"
#ifndef CT_SPECIESTHERMOINTERPTYPE_H
#define CT_SPECIESTHERMOINTERPTYPE_H
namespace Cantera {
/**
* Base class.
*/
class SpeciesThermoInterpType {
public:

View File

@@ -164,9 +164,7 @@ namespace Cantera {
scale(m_ym.begin(), m_ym.end(), c, m_dens);
}
doublereal State::mean_X(const doublereal* Q) const {
return m_mmw*dot(m_ym.begin(), m_ym.end(), Q);
}
doublereal State::mean_Y(const doublereal* Q) const {
return dot(m_y.begin(), m_y.end(), Q);

View File

@@ -193,7 +193,9 @@ namespace Cantera {
* Array Q should contain pure-species molar property
* values.
*/
doublereal mean_X(const doublereal* Q) const;
doublereal mean_X(const doublereal* Q) const {
return m_mmw*inner_product(m_ym.begin(), m_ym.end(), Q, 0.0);
}
/**
* Evaluate the mass-fraction-weighted mean of Q:

View File

@@ -105,6 +105,7 @@ namespace Cantera {
* Utilities for Solvent ID and Molality
* Here we also calculate and store the molecular weight
* of the solvent and the m_Mnaught parameter.
* @param k index of the solvent.
*/
void MolalityVPSSTP::setSolvent(int k) {
if (k < 0 || k >= m_kk) {
@@ -122,8 +123,9 @@ namespace Cantera {
return m_indexSolvent;
}
/*
* Sets the minimum mole fraction in the molality formulation
/**
* Sets the minimum mole fraction in the molality formulation. The
* minimum mole fraction must be in the range 0 to 0.9.
*/
void MolalityVPSSTP::
setMoleFSolventMin(doublereal xmolSolventMIN) {
@@ -135,9 +137,8 @@ namespace Cantera {
m_xmolSolventMIN = xmolSolventMIN;
}
/*
* Returns the minimum mole fraction in the molality
* formulation.
/**
* Returns the minimum mole fraction in the molality formulation.
*/
doublereal MolalityVPSSTP::moleFSolventMin() const {
return m_xmolSolventMIN;
@@ -146,14 +147,16 @@ namespace Cantera {
/**
* getMolalities():
* We calculate the vector of molalities of the species
* in the phase
* m_i = (n_i) / (1000 * M_o * n_o_p)
*
* where M_o is the molecular weight of the solvent
* n_o is the mole fraction of the solvent
* n_i is the mole fraction of the solute.
* n_o_p = max (n_o_min, n_o)
* n_o_min = minimum mole fraction of solvent allowed
* in the phase
* \f[
* m_i = (n_i) / (1000 * M_o * n_{o,p})
* \f]
* where
* - \f$ M_o \f$ is the molecular weight of the solvent
* - \f$ n_o \f$ is the mole fraction of the solvent
* - \f$ n_i \f$ is the mole fraction of the solute.
* - \f$ n_{o,p} = max (n_{o, min}, n_o) \f$
* - \f$ n_{o,min} \f$ = minimum mole fraction of solvent allowed
* in the denominator.
*/
void MolalityVPSSTP::getMolalities(doublereal * const molal) const {
@@ -163,12 +166,12 @@ namespace Cantera {
xmolSolvent = m_xmolSolventMIN;
}
double denomInv = 1.0/
(m_Mnaught * xmolSolvent);
(m_Mnaught * xmolSolvent);
for (int k = 0; k < m_kk; k++) {
molal[k] *= denomInv;
molal[k] *= denomInv;
}
for (int k = 0; k < m_kk; k++) {
m_molalities[k] = molal[k];
m_molalities[k] = molal[k];
}
}
@@ -191,25 +194,25 @@ namespace Cantera {
double Lsum = 1.0 / m_Mnaught;
for (int k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
m_molalities[k] = molal[k];
Lsum += molal[k];
}
if (k != m_indexSolvent) {
m_molalities[k] = molal[k];
Lsum += molal[k];
}
}
double tmp = 1.0 / Lsum;
m_molalities[m_indexSolvent] = tmp / m_Mnaught;
double sum = m_molalities[m_indexSolvent];
for (int k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
m_molalities[k] = tmp * molal[k];
sum += m_molalities[k];
}
if (k != m_indexSolvent) {
m_molalities[k] = tmp * molal[k];
sum += m_molalities[k];
}
}
if (sum != 1.0) {
tmp = 1.0 / sum;
for (int k = 0; k < m_kk; k++) {
m_molalities[k] *= tmp;
}
tmp = 1.0 / sum;
for (int k = 0; k < m_kk; k++) {
m_molalities[k] *= tmp;
}
}
setMoleFractions(DATA_PTR(m_molalities));
/*
@@ -219,7 +222,7 @@ namespace Cantera {
*/
getMolalities(DATA_PTR(m_molalities));
}
/*
* setMolalitiesByName()
*
@@ -257,48 +260,48 @@ namespace Cantera {
double cNeg = 0.0;
double sum = 0.0;
for (int k = 0; k < kk; k++) {
double ch = charge(k);
if (mf[k] > 0.0) {
if (ch > 0.0) {
if (ch * mf[k] > cPos) {
largePos = k;
cPos = ch * mf[k];
}
}
if (ch < 0.0) {
if (fabs(ch) * mf[k] > cNeg) {
largeNeg = k;
cNeg = fabs(ch) * mf[k];
}
}
}
sum += mf[k] * ch;
double ch = charge(k);
if (mf[k] > 0.0) {
if (ch > 0.0) {
if (ch * mf[k] > cPos) {
largePos = k;
cPos = ch * mf[k];
}
}
if (ch < 0.0) {
if (fabs(ch) * mf[k] > cNeg) {
largeNeg = k;
cNeg = fabs(ch) * mf[k];
}
}
}
sum += mf[k] * ch;
}
if (sum != 0.0) {
if (sum > 0.0) {
if (cPos > sum) {
mf[largePos] -= sum / charge(largePos);
} else {
throw CanteraError("MolalityVPSSTP:setMolalitiesbyName",
"unbalanced charges");
}
} else {
if (cNeg > (-sum)) {
mf[largeNeg] -= (-sum) / fabs(charge(largeNeg));
} else {
throw CanteraError("MolalityVPSSTP:setMolalitiesbyName",
"unbalanced charges");
}
}
if (sum > 0.0) {
if (cPos > sum) {
mf[largePos] -= sum / charge(largePos);
} else {
throw CanteraError("MolalityVPSSTP:setMolalitiesbyName",
"unbalanced charges");
}
} else {
if (cNeg > (-sum)) {
mf[largeNeg] -= (-sum) / fabs(charge(largeNeg));
} else {
throw CanteraError("MolalityVPSSTP:setMolalitiesbyName",
"unbalanced charges");
}
}
}
sum = 0.0;
for (int k = 0; k < kk; k++) {
sum += mf[k];
sum += mf[k];
}
sum = 1.0/sum;
for (int k = 0; k < kk; k++) {
mf[k] *= sum;
mf[k] *= sum;
}
setMoleFractions(DATA_PTR(mf));
/*
@@ -314,15 +317,15 @@ namespace Cantera {
*
* Set the molalities of the solutes by name
*/
void MolalityVPSSTP::setMolalitiesByName(const string& x) {
compositionMap xx;
int kk = nSpecies();
for (int k = 0; k < kk; k++) {
xx[speciesName(k)] = -1.0;
}
parseCompString(x, xx);
setMolalitiesByName(xx);
}
void MolalityVPSSTP::setMolalitiesByName(const string& x) {
compositionMap xx;
int kk = nSpecies();
for (int k = 0; k < kk; k++) {
xx[speciesName(k)] = -1.0;
}
parseCompString(x, xx);
setMolalitiesByName(xx);
}
/*

View File

@@ -169,15 +169,15 @@ namespace Cantera {
*/
/**
* Returns the vector of nondimensional
* enthalpies of the reference state at the current temperature
* of the solution and the reference pressure for the species.
* Returns the vector of nondimensional enthalpies of the
* reference state at the current temperature of the solution and
* the reference pressure for the species.
*/
void VPStandardStateTP::getEnthalpy_RT_ref(doublereal *hrt) const {
/*
* Call the function that makes sure the local copy of
* the species reference thermo functions are up to date
* for the current temperature.
* Call the function that makes sure the local copy of the
* species reference thermo functions are up to date for the
* current temperature.
*/
_updateRefStateThermo();
/*

View File

@@ -1,6 +1,6 @@
/**
* @file TransportBase.h
* Provides class Transport.
* Provides class Transport.
*/
/*