From bb19de47bce9544afe2859cb661a0e57394531ab Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Wed, 3 May 2006 19:45:39 +0000 Subject: [PATCH] added electric potential --- Cantera/clib/src/ct.cpp | 5 + Cantera/clib/src/ct.h | 1 + Cantera/python/Cantera/Interface.py | 4 +- Cantera/python/Cantera/ThermoPhase.py | 9 +- Cantera/python/ctml_writer.py | 48 ++++++++ Cantera/python/src/ctthermo_methods.cpp | 2 + Cantera/src/GeneralSpeciesThermo.cpp | 14 +-- Cantera/src/GeneralSpeciesThermo.h | 4 +- Cantera/src/Mu0Poly.cpp | 136 ++++++++++----------- Cantera/src/NasaThermo.h | 37 +++--- Cantera/src/SpeciesThermo.h | 34 ++++-- Cantera/src/SpeciesThermoFactory.cpp | 98 ++++++++++----- Cantera/src/SpeciesThermoFactory.h | 35 +++++- Cantera/src/SpeciesThermoInterpType.h | 5 +- Cantera/src/State.cpp | 4 +- Cantera/src/State.h | 4 +- Cantera/src/thermo/MolalityVPSSTP.cpp | 147 ++++++++++++----------- Cantera/src/thermo/VPStandardStateTP.cpp | 12 +- Cantera/src/transport/TransportBase.h | 2 +- 19 files changed, 375 insertions(+), 226 deletions(-) diff --git a/Cantera/clib/src/ct.cpp b/Cantera/clib/src/ct.cpp index be12cf68b..6f8d8ef00 100755 --- a/Cantera/clib/src/ct.cpp +++ b/Cantera/clib/src/ct.cpp @@ -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(); diff --git a/Cantera/clib/src/ct.h b/Cantera/clib/src/ct.h index 38291af0d..730fa04a5 100755 --- a/Cantera/clib/src/ct.h +++ b/Cantera/clib/src/ct.h @@ -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); diff --git a/Cantera/python/Cantera/Interface.py b/Cantera/python/Cantera/Interface.py index 447b19cf0..19073b963 100644 --- a/Cantera/python/Cantera/Interface.py +++ b/Cantera/python/Cantera/Interface.py @@ -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 diff --git a/Cantera/python/Cantera/ThermoPhase.py b/Cantera/python/Cantera/ThermoPhase.py index 7f46334a9..5224ae88b 100644 --- a/Cantera/python/Cantera/ThermoPhase.py +++ b/Cantera/python/Cantera/ThermoPhase.py @@ -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. diff --git a/Cantera/python/ctml_writer.py b/Cantera/python/ctml_writer.py index ef38bffe9..6d7796796 100644 --- a/Cantera/python/ctml_writer.py +++ b/Cantera/python/ctml_writer.py @@ -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.""" diff --git a/Cantera/python/src/ctthermo_methods.cpp b/Cantera/python/src/ctthermo_methods.cpp index 8da921e6a..50dd2fe75 100644 --- a/Cantera/python/src/ctthermo_methods.cpp +++ b/Cantera/python/src/ctthermo_methods.cpp @@ -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: diff --git a/Cantera/src/GeneralSpeciesThermo.cpp b/Cantera/src/GeneralSpeciesThermo.cpp index 0d6bc4958..0b6d24e11 100644 --- a/Cantera/src/GeneralSpeciesThermo.cpp +++ b/Cantera/src/GeneralSpeciesThermo.cpp @@ -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, diff --git a/Cantera/src/GeneralSpeciesThermo.h b/Cantera/src/GeneralSpeciesThermo.h index 8ac2e6d35..a26d43da5 100644 --- a/Cantera/src/GeneralSpeciesThermo.h +++ b/Cantera/src/GeneralSpeciesThermo.h @@ -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 { diff --git a/Cantera/src/Mu0Poly.cpp b/Cantera/src/Mu0Poly.cpp index 57bfeb66a..1cc1f2857 100644 --- a/Cantera/src/Mu0Poly.cpp +++ b/Cantera/src/Mu0Poly.cpp @@ -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]; diff --git a/Cantera/src/NasaThermo.h b/Cantera/src/NasaThermo.h index 22d19b9a0..87ebdccca 100755 --- a/Cantera/src/NasaThermo.h +++ b/Cantera/src/NasaThermo.h @@ -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 &mlg = m_low[grp-1]; - const vector &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 &mlg = m_low[grp-1]; + const vector &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 m_low_map; - // mutable map m_high_map; vector > m_high; vector > m_low; map m_index; diff --git a/Cantera/src/SpeciesThermo.h b/Cantera/src/SpeciesThermo.h index 56c972bb2..c19193cd1 100755 --- a/Cantera/src/SpeciesThermo.h +++ b/Cantera/src/SpeciesThermo.h @@ -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; diff --git a/Cantera/src/SpeciesThermoFactory.cpp b/Cantera/src/SpeciesThermoFactory.cpp index b873cecbd..867a0e1a1 100755 --- a/Cantera/src/SpeciesThermoFactory.cpp +++ b/Cantera/src/SpeciesThermoFactory.cpp @@ -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 sp; + + // get all of the species nodes sparray.getChildren("species",sp); - int ns = static_cast(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 nodes) { int n = static_cast(nodes.size()); @@ -99,6 +107,9 @@ namespace Cantera { } + /** + * @todo is this used? + */ SpeciesThermo* SpeciesThermoFactory:: newSpeciesThermoOpt(vector nodes) { int n = static_cast(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; diff --git a/Cantera/src/SpeciesThermoFactory.h b/Cantera/src/SpeciesThermoFactory.h index d1e0e3fa6..546751fee 100755 --- a/Cantera/src/SpeciesThermoFactory.h +++ b/Cantera/src/SpeciesThermoFactory.h @@ -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 nodes); virtual SpeciesThermo* newSpeciesThermoOpt(vector 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) { diff --git a/Cantera/src/SpeciesThermoInterpType.h b/Cantera/src/SpeciesThermoInterpType.h index eff84b356..9a443e92f 100644 --- a/Cantera/src/SpeciesThermoInterpType.h +++ b/Cantera/src/SpeciesThermoInterpType.h @@ -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: diff --git a/Cantera/src/State.cpp b/Cantera/src/State.cpp index 4d110b0f6..06d952650 100644 --- a/Cantera/src/State.cpp +++ b/Cantera/src/State.cpp @@ -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); diff --git a/Cantera/src/State.h b/Cantera/src/State.h index e17ab737c..d2e914773 100755 --- a/Cantera/src/State.h +++ b/Cantera/src/State.h @@ -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: diff --git a/Cantera/src/thermo/MolalityVPSSTP.cpp b/Cantera/src/thermo/MolalityVPSSTP.cpp index 17636eb51..aef04292a 100644 --- a/Cantera/src/thermo/MolalityVPSSTP.cpp +++ b/Cantera/src/thermo/MolalityVPSSTP.cpp @@ -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); + } /* diff --git a/Cantera/src/thermo/VPStandardStateTP.cpp b/Cantera/src/thermo/VPStandardStateTP.cpp index 0b71ec6c0..870817d2b 100644 --- a/Cantera/src/thermo/VPStandardStateTP.cpp +++ b/Cantera/src/thermo/VPStandardStateTP.cpp @@ -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(); /* diff --git a/Cantera/src/transport/TransportBase.h b/Cantera/src/transport/TransportBase.h index 39fc37305..3d188bd94 100755 --- a/Cantera/src/transport/TransportBase.h +++ b/Cantera/src/transport/TransportBase.h @@ -1,6 +1,6 @@ /** * @file TransportBase.h - * Provides class Transport. + * Provides class Transport. */ /*