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

This checkin adds the ability to read NASA9 polynomial standard states,

Browse source 
with arbitrary number of temperature intervals,
to Cantera's underlying routines.

Cantera can now read fortran formatted
NASA9 polynomials to create cti files

The cti files can then be translated into xml files.

And, the xml files can now be read into constant pressure standard state
objects and used normally within the guts of Cantera.

Currently, the standard state just links into the GeneralSpeciesThermo
object, which means that it's calculation speed is slow. However, atm
this satisfies the initial use case for this new capability.

In the near future, I'll push this out to the matlab and python interfaces.
This commit is contained in:
Harry Moffat 2007-09-13 15:05:39 +00:00
parent b7e616a2a2
commit f70c394fa6
36 changed files with 2676 additions and 499 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

8
Cantera/src/base/xml.h
View file

@ -124,7 +124,13 @@ namespace Cantera {
XML_Node* parent() const { return m_parent; }
XML_Node* setParent(XML_Node* p) { m_parent = p; return p; }
bool hasChild(std::string ch) const {
//! Tests whether the current node has a child node with a particular name
/*!
* @param ch Name of the child node to test
*
* @return Returns true if the child node exists, false otherwise.
*/
bool hasChild(const std::string ch) const {
return (m_childindex.find(ch) != m_childindex.end());
}
bool hasAttrib(std::string a) const {

321
Cantera/src/converters/CKParser.cpp
View file

@ -5,40 +5,7 @@
// Copyright 2001 California Institute of Technology
//
// $Log$
// Revision 1.20 2005-12-09 17:49:34 dggoodwin
// removed critical and saturation properties from ThermoPhase
//
// Revision 1.19 2005/07/28 23:02:44 hkmoffa
// Got rid of one warning message.
//
// Revision 1.18 2005/07/26 03:56:35 dggoodwin
// cleanup
//
// Revision 1.17 2005/07/25 03:51:21 dggoodwin
// now recognizes the FORD keyword
//
// Revision 1.16 2005/01/07 10:26:43 dggoodwin
// merged changes from branch
//
// Revision 1.15.2.2 2004/12/18 15:16:13 dggoodwin
// minor cleanup
//
// Revision 1.15.2.1 2004/12/18 15:00:02 dggoodwin
// *** empty log message ***
//
//
// Revision 1.6 2004/07/02 16:48:13 hkmoffa
// Moved CK_SyntaxError definition to the .h file. It's used in more
// than one .cpp file.
//
// Revision 1.5 2004/05/13 17:45:05 dggoodwin
// fixed bug in which a species name beginning with M was interpreted as a third body
//
// Revision 1.4 2004/05/13 16:58:33 dggoodwin
// *** empty log message ***
//
//
// $Id$
// turn off warnings about truncating long names under Windows
#ifdef WIN32
@ -62,29 +29,29 @@ using namespace std;
namespace ckr {
static string int2s(int n, string fmt="%d") {
static string int2s(int n, std::string fmt="%d") {
char buf[30];
sprintf(buf, fmt.c_str(), n);
return string(buf);
}
/// Exception class for syntax errors.
CK_SyntaxError::CK_SyntaxError(ostream& f,
const string& s, int linenum)
CK_SyntaxError::CK_SyntaxError(std::ostream& f,
const std::string& s, int linenum)
: m_out(f) {
m_msg += "Syntax error: " + s;
if (linenum > 0) m_msg += " (line " + int2s(linenum) + ")\n";
}
static int parseGroupString(string str, vector<string>& esyms,
static int parseGroupString(std::string str, std::vector<std::string>& esyms,
vector_int& result);
/**
* Throw an exception if one of the four lines that must have
* 1, 2, 3, or 4 in column 80 do not.
*/
static void illegalThermoLine(ostream& f,
static void illegalThermoLine(std::ostream& f,
char n, int linenum = -1)
{
throw CK_SyntaxError(f, "column 80 must "
@ -95,15 +62,15 @@ namespace ckr {
/**
* Throw an exception if number string is bad
*/
static void illegalNumber(ostream& f,
string s, int linenum = -1)
static void illegalNumber(std::ostream& f,
std::string s, int linenum = -1)
{
string msg = "illegal number: "+s;
throw CK_SyntaxError(f, msg, linenum);
};
extern void getDefaultAtomicWeights(map<string,double>& weights);
extern void getDefaultAtomicWeights(std::map<std::string,double>& weights);
static string d2e(string s) {
size_t n;
@ -118,14 +85,34 @@ namespace ckr {
return r;
}
static double de_atof(string s) {
static double de_atof(std::string s) {
string r = d2e(s);
//double rval = Cantera::atofCheck(r.c_str());
double rval = atof(r.c_str());
return rval;
}
static double getNumberFromString(string s) {
/**
* Check validity of the temperatures defining the
* temperature ranges for the NASA9 polynomial species thermodynamic
* property fits.
* @param log log file output stream
* @param temp Vector of temperatures
*/
static void checkNASA9Temps(std::ostream& log, vector_fp &temp) {
int i;
for (i = 1; i < (int) temp.size(); i++) {
double tlow = temp[i-1];
double thigh = temp[i];
if (thigh <= tlow) {
string sss = "error reading temperature";
throw CK_SyntaxError(log, sss);
}
}
}
static double getNumberFromString(std::string s) {
bool inexp = false;
removeWhiteSpace(s);
int sz = static_cast<int>(s.size());
@ -155,8 +142,8 @@ namespace ckr {
* @param sp Species object to add element to
* @param log log file output stream
*/
static void addElement(string symbol, double atoms,
Species& sp, ostream& log) {
static void addElement(std::string symbol, double atoms,
Species& sp, std::ostream& log) {
if (atoms != 0.0) {
Constituent e;
@ -177,8 +164,8 @@ namespace ckr {
* @param tmid intermediate temperature
* @param tmax maximum temperature
*/
static void checkTemps(ostream& log, double tmin,
double tmid, double tmax)
static void checkTemps(std::ostream& log, double tmin,
double tmid, double tmax)
{
if (tmin == 0.0 || tmid == 0.0 || tmax == 0.0) {
throw CK_SyntaxError(log,
@ -186,8 +173,9 @@ namespace ckr {
}
}
static void getSpecies(string s,
int n, vector<RxnSpecies>& species, bool debug, ostream& log) {
static void getSpecies(std::string s,
int n, vector<RxnSpecies>& species, bool debug,
std::ostream& log) {
removeWhiteSpace(s);
// break string into substrings at the '+' characters separating
// species symbols
@ -247,8 +235,9 @@ namespace ckr {
* containing the species symbols and stoichiometric coefficients.
* @todo allow non-integral stoichiometric coefficients
*/
int getGroups(string::const_iterator begin, string::const_iterator end, vector<string>& esyms,
vector<grouplist_t>& rxngroups)
int getGroups(std::string::const_iterator begin,
std::string::const_iterator end, std::vector<std::string>& esyms,
std::vector<grouplist_t>& rxngroups)
{
bool ingroup = false;
rxngroups.clear();
@ -284,29 +273,36 @@ namespace ckr {
/**
* Constructor. Construct a parser for the specified input file.
*/
CKParser::CKParser(istream* infile, const string& fname, ostream* log)
: verbose(true), debug(false), m_line (0) {
m_ckfile = infile;
m_ckfilename = fname;
m_log = log;
m_nasafmt = false;
m_last_eol = '\n';
}
CKParser::CKParser(std::istream* infile, const std::string& fname,
std::ostream* log) :
verbose(true),
debug(false),
m_line (0),
m_nasafmt(false),
m_nasa9fmt(false)
{
m_ckfile = infile;
m_ckfilename = fname;
m_log = log;
m_last_eol = '\n';
}
/**
*
* Get a line from the input file, and return it in string s. If the
* line contains a comment character (!), then return only the
* portion preceding it. Non-printing characters are replaced by
* spaces.
* @param s On return, s contains the line read from the
* input file.
* @param comment On return, comment contains the text following the
* comment character on the line, if any.
*
*/
void CKParser::getCKLine(string& s, string& comment) {
// Get a line from the input file, and return it in string s.
/*
* If the line contains a comment character (!), then return only the
* portion preceding it. Non-printing characters are replaced by
* spaces.
*
* The input file is m_ckfile, an istream.
*
* @param s On return, s contains the line read from the
* input file.
* @param comment On return, comment contains the text following the
* comment character on the line, if any.
*/
void CKParser::getCKLine(std::string& s, std::string& comment) {
// Chemkin comment character
const char commentChar = '!';
@ -398,13 +394,13 @@ namespace ckr {
*
*/
void CKParser::putCKLine(string& s, string& comment) {
void CKParser::putCKLine(std::string& s, std::string& comment) {
m_buf = s;
m_comment = comment;
}
bool CKParser::advanceToKeyword(const string& kw, const string& stop) {
bool CKParser::advanceToKeyword(const std::string& kw, const std::string& stop) {
string s, c;
do {
getCKLine(s,c);
@ -573,9 +569,9 @@ next:
*
*/
bool CKParser::readThermoSection(vector<string>& names,
speciesTable& species, vector_fp& temp,
int& optionFlag, ostream& log) {
bool CKParser::readThermoSection(std::vector<std::string>& names,
speciesTable& species, vector_fp& temp,
int& optionFlag, std::ostream& log) {
string s;
vector<string> toks;
@ -614,39 +610,97 @@ next:
}
else if (match(toks[itt],"NO_TMID")) {
m_nasafmt = true;
log << "\nOption 'NO_TMID' specified. Default midpoint temperature\n";
log << "\nOption 'NO_TMID' specified. Default "
"midpoint temperature\n";
log << "will be used for all species.\n\n";
}
else throw CK_SyntaxError(log,
} else if (match(toks[itt], "NASA9")) {
m_nasa9fmt = true;
log << "Option NASA9 specified: Use new "
"nasa input file format\n\n";
} else if (match(toks[itt], "NASA")) {
m_nasa9fmt = false;
log << "Option NASA specified: Use old "
"nasa input file format\n\n";
} else throw CK_SyntaxError(log,
"unrecognized THERMO option.", m_line);
}
}
// if "THERMO ALL" specified, or if optionFlag is set to HasTempRange,
// then the next line must be the 3 default temperatures for the database.
// then the next line must contain the default temperatures
// for the database.
if (optionFlag == NoThermoDatabase || optionFlag == HasTempRange) {
getCKLine(s, comment);
getTokens(s, static_cast<int>(s.size()), toks);
if (toks.size() >= 3) {
if (m_nasa9fmt) {
//
// For NASA9 polynomials, the format is
// t1 t2 t3 t4 date
// when there are 3 temperature regions
//
int nreg = toks.size() - 2;
if (nreg >= 1) {
temp.resize(nreg+1);
for (int i = 0; i <= nreg; i++) {
temp[i] = de_atof(toks[i]);
}
string defaultDate = toks[nreg+1];
} else {
throw CK_SyntaxError(log, "Default temp region card is bad", m_line);
}
if (verbose) {
log.flags(ios::showpoint | ios::fixed);
log.precision(2);
log << endl << " Default # of temperature regions: "
<< nreg << endl;
log << " ";
for (int i = 0; i <= nreg; i++) {
log << temp[i] << " ";
}
log << endl;
}
checkNASA9Temps(log, temp);
} else {
//
// For NASA polynomials, the format is
// tlow tmid thigh
// There are always 2 temperature regions
//
if (toks.size() >= 3) {
tmin = de_atof(toks[0]);
tmid = de_atof(toks[1]);
tmax = de_atof(toks[2]);
}
}
if (verbose) {
if (verbose) {
log.flags(ios::showpoint | ios::fixed);
log.precision(2);
log << endl << " default Tlow, Tmid, Thigh: " << tmin << " "
<< tmid << " " << tmax << endl;
}
checkTemps(log, tmin, tmid, tmax);
temp.clear();
temp.push_back(tmin);
temp.push_back(tmid);
temp.push_back(tmax);
}
checkTemps(log, tmin, tmid, tmax);
temp.clear();
temp.push_back(tmin);
temp.push_back(tmid);
temp.push_back(tmax);
}
}
/// XXXX BRANCH TO THE DIFFERENT THERMO READERS HERE
// Check to see that we expect to be reading a NASA formatted file
if (m_nasa9fmt) {
bool ok = readNASA9ThermoSection(names, species, temp,
optionFlag, log);
if (!ok) {
throw CK_SyntaxError(log,
"In NASA parser. However, we expect a NASA9 file format",
-1);
}
return ok;
}
// now read in all species records that have names in list 'names'
bool getAllSpecies = (nsp > 0 && match(names[0],"<ALL>"));
@ -844,7 +898,7 @@ next:
void CKParser::missingAuxData(const string& kw) {
void CKParser::missingAuxData(const std::string& kw) {
throw CK_SyntaxError(*m_log, kw +
" keyword must be followed by slash-delimited data.", m_line);
}
@ -854,7 +908,7 @@ next:
* Parse the REACTION section of the input file, and return
* a list of Reaction objects and the units.
*/
bool CKParser::readReactionSection(const vector<string>& speciesNames,
bool CKParser::readReactionSection(const std::vector<std::string>& speciesNames,
vector<string>& elementNames, reactionList& reactions,
ReactionUnits& units)
{
@ -952,44 +1006,44 @@ next:
// look for a metadata line
if (s[0] == '%') {
metaDataLine = true;
if (eqloc > 0 && eqloc < int(s.size())) {
int ierr, ierp;
vector<grouplist_t> rg, pg;
s[eqloc] = ' ';
ierr = getGroups(s.begin(), s.begin() + eqloc,
elementNames, rg);
ierp = getGroups(s.begin() + eqloc, s.end(),
elementNames, pg);
unsigned int nr =
static_cast<unsigned int>(rxn.reactants.size());
unsigned int nratoms = 0;
for (unsigned int ij = 0; ij < nr; ij++)
nratoms += int(rxn.reactants[ij].number);
if (rg.size() != nratoms)
throw CK_SyntaxError(*m_log,
" groups not specified for all reactants", m_line);
else if (ierr < 0)
throw CK_SyntaxError(*m_log,
" error in reactant group specification", m_line);
for (unsigned int ir = 0; ir < nr; ir++) {
rxn.reactants[ir].groups = rg[ir];
}
unsigned int np =
static_cast<unsigned int>(rxn.products.size());
unsigned int npatoms = 0;
for (unsigned int ik = 0; ik < np; ik++)
npatoms += int(rxn.products[ik].number);
if (pg.size() != npatoms)
throw CK_SyntaxError(*m_log,
" groups not specified for all products", m_line);
else if (ierp < 0)
throw CK_SyntaxError(*m_log,
" error in product group specification", m_line);
for (unsigned int ip = 0; ip < np; ip++) {
rxn.products[ip].groups = pg[ip];
}
}
metaDataLine = true;
if (eqloc > 0 && eqloc < int(s.size())) {
int ierr, ierp;
vector<grouplist_t> rg, pg;
s[eqloc] = ' ';
ierr = getGroups(s.begin(), s.begin() + eqloc,
elementNames, rg);
ierp = getGroups(s.begin() + eqloc, s.end(),
elementNames, pg);
unsigned int nr =
static_cast<unsigned int>(rxn.reactants.size());
unsigned int nratoms = 0;
for (unsigned int ij = 0; ij < nr; ij++)
nratoms += int(rxn.reactants[ij].number);
if (rg.size() != nratoms)
throw CK_SyntaxError(*m_log,
" groups not specified for all reactants", m_line);
else if (ierr < 0)
throw CK_SyntaxError(*m_log,
" error in reactant group specification", m_line);
for (unsigned int ir = 0; ir < nr; ir++) {
rxn.reactants[ir].groups = rg[ir];
}
unsigned int np =
static_cast<unsigned int>(rxn.products.size());
unsigned int npatoms = 0;
for (unsigned int ik = 0; ik < np; ik++)
npatoms += int(rxn.products[ik].number);
if (pg.size() != npatoms)
throw CK_SyntaxError(*m_log,
" groups not specified for all products", m_line);
else if (ierp < 0)
throw CK_SyntaxError(*m_log,
" error in product group specification", m_line);
for (unsigned int ip = 0; ip < np; ip++) {
rxn.products[ip].groups = pg[ip];
}
}
}
else if (eqloc >= 0 && eqloc < int(s.size())) {
@ -1217,7 +1271,8 @@ next:
// check for duplicate keyword
if (kwindex[name]) {
throw CK_SyntaxError(*m_log,
"duplicate auxiliary data keyword " + name, m_line);
"duplicate auxiliary data keyword "
+ name, m_line);
}
else
kwindex[name] = 1;
@ -1415,7 +1470,7 @@ next:
int parseGroupString(string str, vector<string>& esyms, group_t& result) {
int parseGroupString(std::string str, std::vector<std::string>& esyms, group_t& result) {
bool inSymbol=true;
string s = str + '-';
int i;

153
Cantera/src/converters/CKParser.h
View file

@ -16,7 +16,7 @@
#include <fstream>
#include <string>
#include <iostream>
using namespace std;
//using namespace std;
#include "ckr_defs.h"
#include "Element.h"
@ -26,70 +26,109 @@ using namespace std;
namespace ckr {
// typedefs
// typedefs
/// readability constants
//@{
const int NoThermoDatabase = 10;
const int HasTempRange = 11;
//@}
/// readability constants
//@{
const int NoThermoDatabase = 10;
const int HasTempRange = 11;
//@}
/// Exception class for syntax errors.
class CK_SyntaxError : public CK_Exception {
public:
CK_SyntaxError(ostream& f, const string& s, int linenum = -1);
ostream& m_out;
};
/// Exception class for syntax errors.
class CK_SyntaxError : public CK_Exception {
public:
CK_SyntaxError(std::ostream& f, const std::string& s, int linenum = -1);
std::ostream& m_out;
};
/**
* Chemkin mechanism file parser. For internal use by class CKReader.
/**
* Chemkin mechanism file parser. For internal use by class CKReader.
*/
class CKParser {
public:
CKParser(std::string ckfile, std::ostream* log);
CKParser(std::istream* infile, const std::string& fname,
std::ostream* log);
/// Destructor.
~CKParser(){}
bool readElementSection(elementList& elements);
bool readSpeciesSection(speciesList& species);
bool readThermoSection(vector<string>& names,
speciesTable& speciesData, vector_fp& temp,
int& optionFlag, ostream& log);
bool readReactionSection(const vector<string>& speciesNames,
vector<string>& elementNames,
reactionList& reactions, ReactionUnits& units);
bool advanceToKeyword(const std::string& kw, const std::string& stop);
bool verbose;
bool debug;
bool readNASA9ThermoSection(std::vector<string>& names,
speciesTable& species, vector_fp& temp,
int& optionFlag, std::ostream& log);
void readNASA9ThermoRecord(Species& sp);
private:
//! Local value of the line number being read
/*!
* This is used for debug IO printout purposes
*/
class CKParser {
public:
CKParser(string ckfile, ostream* log);
CKParser(istream* infile, const string& fname, ostream* log);
/// Destructor.
~CKParser(){}
int m_line;
bool readElementSection(elementList& elements);
bool readSpeciesSection(speciesList& species);
bool readThermoSection(vector<string>& names,
speciesTable& speciesData, vector_fp& temp,
int& optionFlag, ostream& log);
bool readReactionSection(const vector<string>& speciesNames,
vector<string>& elementNames,
reactionList& reactions, ReactionUnits& units);
bool advanceToKeyword(const string& kw, const string& stop);
bool verbose;
bool debug;
private:
std::string m_buf;
std::string m_comment;
//! This is the input file that is read
/*!
* It's an istream
*/
std::istream* m_ckfile;
int m_line;
string m_buf;
string m_comment;
istream* m_ckfile;
string m_ckfilename;
ostream* m_log;
bool m_nasafmt;
char m_last_eol;
void readThermoRecord(Species& sp);
void getCKLine(string& s, string& comment);
void putCKLine(string& s, string& comment);
void missingAuxData(const string& kw);
};
std::string m_ckfilename;
//! Pointer to the ostream for writing debugging output log info
std::ostream* m_log;
bool m_nasafmt;
//! Boolean indicating new NASA input file format
/*!
* If this is true, a completely different input file parser is
* used.
*/
bool m_nasa9fmt;
char m_last_eol;
void readThermoRecord(Species& sp);
//! Get a line from the input file, and return it in string s.
/*!
* If the line contains a comment character (!), then return only the
* portion preceding it. Non-printing characters are replaced by
* spaces.
*
* The input file is m_ckfile, an istream.
*
* @param s On return, s contains the line read from the
* input file.
* @param comment On return, comment contains the text following the
* comment character on the line, if any.
*/
void getCKLine(std::string& s, std::string& comment);
void putCKLine(std::string& s, std::string& comment);
void missingAuxData(const std::string& kw);
void checkSpeciesName(std::string spname);
};
}
#endif

21
Cantera/src/converters/CKReader.cpp
View file

@ -42,8 +42,8 @@ namespace ckr {
* @return true if no errors were encountered, false otherwise
*/
bool CKReader::read(const string& inputFile, const string& thermoDatabase,
const string& logfile) {
bool CKReader::read(const std::string& inputFile, const std::string& thermoDatabase,
const std::string& logfile) {
clock_t t0, t1;
@ -174,8 +174,9 @@ bool CKReader::read(const string& inputFile, const string& thermoDatabase,
string nm;
vector<string> undef;
bool allsp = (speciesSymbols[0] == "<ALL>");
if (hasthermo && parser.readThermoSection(speciesSymbols,
speciesData, temp, optionFlag, log)) {
if (hasthermo &&
parser.readThermoSection(speciesSymbols,
speciesData, temp, optionFlag, log)) {
if (allsp) {
nsp = static_cast<int>(speciesData.size());
for (k = 0; k < nsp; k++) {
@ -294,7 +295,7 @@ bool CKReader::read(const string& inputFile, const string& thermoDatabase,
/// print a summary of all reactions to the log file
bool CKReader::writeReactions(ostream& log) {
bool CKReader::writeReactions(std::ostream& log) {
bool ok = true;
// int ns = species.size();
@ -356,7 +357,7 @@ bool CKReader::writeReactions(ostream& log) {
/// validate the species
bool CKReader::validateSpecies(ostream& log) {
bool CKReader::validateSpecies(std::ostream& log) {
int nel = static_cast<int>(elements.size());
int nsp = static_cast<int>(species.size());
double nm, tol;
@ -420,7 +421,7 @@ bool CKReader::validateSpecies(ostream& log) {
/// validate the reactions
bool CKReader::validateReactions(ostream& log) {
bool CKReader::validateReactions(std::ostream& log) {
bool ok = true;
// int ns = species.size();
@ -479,7 +480,7 @@ bool CKReader::validateReactions(ostream& log) {
* - The heat capacity at Tmax is not greater than the equipartition limit
* for the number of atoms in the molecule
*/
bool checkThermo(ostream& log, speciesList& sp, double tol) {
bool checkThermo(std::ostream& log, speciesList& sp, double tol) {
const double dt = 0.0001;
double t, cp0, h0, s0, cp1, h1, s1;
int nsp = static_cast<int>(sp.size());
@ -635,8 +636,8 @@ bool checkThermo(ostream& log, speciesList& sp, double tol) {
* @return true if all reactions balance
* @todo use reaction number stored in reaction object
*/
bool checkBalance(ostream& f, speciesTable& speciesData,
reactionList& r, vector<int>& unbalanced, double tolerance)
bool checkBalance(std::ostream& f, speciesTable& speciesData,
reactionList& r, std::vector<int>& unbalanced, double tolerance)
{
int nrxn = static_cast<int>(r.size());
string rname, pname;

12
Cantera/src/converters/CKReader.h
View file

@ -33,14 +33,14 @@ namespace ckr {
/// Construct a new empty Group object
Group() : name("<empty>"), index(-1) {}
Group(const string& nm) : name(nm), index(-1) {}
Group(const std::string& nm) : name(nm), index(-1) {}
/// Destructor
~Group() {}
string name; //!< name
std::string name; //!< name
int index; //!< index number
map<string, double> comp; //!< elemental composition
std::map<std::string, double> comp; //!< elemental composition
bool operator==(const Group& g) const {
return (name == g.name);
@ -51,7 +51,7 @@ namespace ckr {
};
/// a list (vector) of Groups
typedef vector<Group> groupList;
typedef std::vector<Group> groupList;
/**
@ -88,8 +88,8 @@ namespace ckr {
* @param logFile file to write logging and error messages to.
* @return true if no errors encountered, false otherwise.
*/
bool read(const string& inputFile,
const string& thermoDatabase, const string& logFile);
bool read(const std::string& inputFile,
const std::string& thermoDatabase, const std::string& logFile);
void write(string outputFile); ///< not implemented.

2
Cantera/src/converters/Makefile.in
View file

@ -20,7 +20,7 @@ PIC_FLAG=@PIC@
CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT) $(PIC_FLAG)
OBJS = atomicWeightDB.o CKParser.o CKReader.o Reaction.o ckr_utils.o \
thermoFunctions.o writelog.o ck2ct.o
thermoFunctions.o writelog.o ck2ct.o Species.o
# ck2ctml.o
CONV_H = CKReader.h thermoFunctions.h \
Element.h Reaction.h CKParser.h \

551
Cantera/src/converters/NASA9Parser.cpp Normal file
View file

@ -0,0 +1,551 @@
/**
* @file CKParser.cpp
*
*/
// Copyright 2001 California Institute of Technology
//
// $Id$
//
#include <numeric>
#include <algorithm>
#include <fstream>
#include <iomanip>
#include <math.h>
#include "CKParser.h"
#include "ckr_utils.h"
#include "writelog.h"
#include <stdio.h>
//#include "../stringUtils.h"
#include <ctype.h>
#include <string>
using namespace std;
namespace ckr {
static std::string int2s(int n, std::string fmt="%d") {
char buf[30];
sprintf(buf, fmt.c_str(), n);
return string(buf);
}
/**
* Add an element to a species.
* @param symbol element symbol
* @param atoms number of atoms of this element in the
* species (may be non-integral)
* @param sp Species object to add element to
* @param log log file output stream
*/
static void addElement(std::string symbol, double atoms,
Species& sp, std::ostream& log) {
if (atoms != 0.0) {
Constituent e;
e.name = symbol;
e.number = atoms;
sp.elements.push_back(e);
sp.comp[symbol] = atoms;
}
}
/**
* Throw an exception if number string is bad
*/
static void illegalNumber(std::ostream& f,
std::string s, int linenum = -1)
{
string msg = "illegal number: "+s;
throw CK_SyntaxError(f, msg, linenum);
};
void CKParser::checkSpeciesName(std::string spname) {
if (spname.size() <= 0) {
string sss = "Empty for string name";
throw CK_SyntaxError(*m_log, sss, m_line);
}
char first = spname[0];
if (isdigit(first)) {
string sss = "First char of string name is number";
throw CK_SyntaxError(*m_log, sss, m_line);
}
if (isspace(first)) {
string sss = "First char of string name is white space";
throw CK_SyntaxError(*m_log, sss, m_line);
}
}
static std::string d2e(std::string s) {
size_t n;
size_t sz = s.size();
string r = s;
char ch;
for (n = 0; n < sz; n++) {
ch = s[n];
if (ch == 'D') r[n] = 'E';
else if (ch == 'd') r[n] = 'e';
}
return r;
}
static double de_atof(std::string s) {
string r = d2e(s);
//double rval = Cantera::atofCheck(r.c_str());
double rval = atof(r.c_str());
return rval;
}
/**
* Check validity of the three temperatures defining the two
* temperature ranges for the NASA polynomial species thermodynamic
* property fits.
* @param log log file output stream
* @param tmin minimum temperature
* @param tmid intermediate temperature
* @param tmax maximum temperature
*/
static void checkTemps(std::ostream& log, double tmin,
double tmid, double tmax)
{
if (tmin == 0.0 || tmid == 0.0 || tmax == 0.0) {
throw CK_SyntaxError(log,
"error reading Tmin, Tmid, or Tmax");
}
}
static double getNumberFromString(std::string s) {
bool inexp = false;
removeWhiteSpace(s);
int sz = static_cast<int>(s.size());
char ch;
for (int n = 0; n < sz; n++) {
ch = s[n];
if (!inexp && (ch == 'E' || ch == 'e' || ch == 'D' || ch == 'd'))
inexp = true;
else if (ch == '+' || ch == '-') {
if (n > 0 && (s[n-1] != 'E' && s[n-1]
!= 'e' && s[n-1] != 'd' && s[n-1] != 'D')) {
return UNDEF;
}
}
else if (ch != '.' && (ch < '0' || ch > '9')) {
return UNDEF;
}
}
return de_atof(s);
}
static int de_atoi(std::ostream &log, std::string s, int line = -1) {
double val = getNumberFromString(s);
int ival = (int) val;
double val2 = (double) ival;
if (fabs(val - val2) >= 0.00001 * (val + val2)) {
string sss = "de_atoi: Conversion of int failed: " + s;
throw CK_SyntaxError(log, sss, line);
}
return ival;
}
/**
* Check validity of the temperatures defining the
* temperature ranges for the NASA9 polynomial species thermodynamic
* property fits.
* @param log log file output stream
* @param temp Vector of temperatures
*/
static void checkNASA9Temps(std::ostream& log, vector_fp &temp) {
int i;
for (i = 1; i <= (int) temp.size(); i++) {
double tlow = temp[i-1];
double thigh = temp[i];
if (thigh <= tlow) {
string sss = "error reading temperature";
throw CK_SyntaxError(log, sss);
}
}
}
/**
*
* Read species data from THERMO section records.
*
* @param names List of species names (input).
* @param species Table of species objects holding data from records
* in THERMO section (output).
* @param temp Devault vector of temperature region boundaries
* There are one more temperatures than there are
* temperature regions.
* @param allowExtThermoData True if 'THERMO' specified, false if
* 'THERMO ALL' specified.
*
* @return True, if the THERMO section exists and the species
* have all been successfully processed. False, if
* the THERMO section doesn't exist or there were
* additional problems.
*/
bool CKParser::readNASA9ThermoSection(std::vector<string>& names,
speciesTable& species, vector_fp& temp,
int& optionFlag, std::ostream& log) {
// String buffer for lines
string s;
vector<string> toks;
string defaultDate="";
int nreg = 2;
int i;
int nsp = static_cast<int>(names.size());
// Comment string
string comment;
// if "THERMO ALL" specified, or if optionFlag is set to HasTempRange,
// then the next line must be the default temperatures for the database.
//
// This line will have nreg+2 tokens on it
// The last token is a date.
if (0) {
if (optionFlag == NoThermoDatabase || optionFlag == HasTempRange) {
getCKLine(s, comment);
getTokens(s, static_cast<int>(s.size()), toks);
nreg = toks.size();
if (nreg >= 1) {
temp.resize(nreg+1);
for (i = 0; i <= nreg; i++) {
temp[i] = de_atof(toks[i]);
}
defaultDate = toks[nreg+1];
}
if (verbose) {
log.flags(ios::showpoint | ios::fixed);
log.precision(2);
log << endl << " Default # of temperature regions: " << nreg << endl;
log << " ";
for (i = 0; i <= nreg; i++) {
log << temp[i] << " ";
}
log << endl;
}
checkNASA9Temps(log, temp);
}
}
// Check to see that we expect to be reading a NASA9 formatted file
if (!m_nasa9fmt) {
throw CK_SyntaxError(log,
"In NASA9 parser. However, we expect a different file format",
-1);
}
// now read in all species records that have names in list 'names'
bool getAllSpecies = (nsp > 0 && match(names[0], "<ALL>"));
if (getAllSpecies) names.clear();
// Map between the number of times a species name appears in the database
map<string, int> dup; // used to check for duplicate THERMO records
bool already_read;
while (1 > 0) {
// If we don't have any more species to read, break
if (nsp == 0) break;
already_read = false;
// Read a new species record from the section
Species spec;
readNASA9ThermoRecord(spec);
// we signal the end of the section by putting <END> as a
// species name. Break if you find the end of the section.
if (spec.name == "<END>") {
break;
}
// check for duplicate thermo data
if (dup[spec.name] == 2) {
log << "Warning: more than one THERMO record for "
<< "species " << spec.name << endl;
log << "Record at line " << m_line
<< " of " << m_ckfilename << " ignored." << endl;
already_read = true;
}
// Set the record in the map to 2 to create a signal for the
// next time.
dup[spec.name] = 2;
// Check to see whether we need this particlar species name
if (!already_read && (getAllSpecies
|| (find(names.begin(), names.end(), spec.name)
< names.end()))) {
// Add the species object to the map. Note we are
// doing a copy constructor here, so we create a
// lasting entry.
species[spec.name] = spec;
if (verbose) {
log << endl << "found species " << spec.name;
log << " at line " << m_line
<< " of " << m_ckfilename;
writeSpeciesData(log, spec);
}
//checkTemps(log, spec.tlow, spec.tmid, spec.thigh);
if (getAllSpecies) {
names.push_back(spec.name);
nsp = static_cast<int>(names.size());
}
else
nsp--;
}
}
return true;
}
/**
*
* Read one species definition in a NASA9 string.
*
*/
void CKParser::readNASA9ThermoRecord(Species& sp) {
string s;
string numstr;
double cf;
// Set to the NASA9 polynomial format
sp.thermoFormatType = 1;
// look for line 1, but if a keyword is found first or the end of
// the file is reached, return "<END>" as the species name
string comment;
// Name of the species
string nameid;
vector<string> toks;
int nToks = 0;
// Loop forward until we get to the next nonempty line.
do {
getCKLine(s, comment);
if (isKeyword(s) || match(s, "<EOF>")) {
sp.name = "<END>";
putCKLine(s, comment);
return;
}
// The first 18 spaces are devoted to the name of the species
string nameid = s.substr(0,18);
getTokens(nameid, static_cast<int>(nameid.size()), toks);
nToks = toks.size();
} while (nToks == 0);
//------------- line 1 ---------------------------
// Everything after the first 18 spaces is a comment.
int nt = s.size();
sp.m_commentsRef = s.substr(18, nt-18);
// Parse the species name
sp.name = toks[0];
sp.id = "";
if (nToks > 1) {
throw CK_SyntaxError(*m_log,
"Illegal number of tokens for string name", m_line);
}
checkSpeciesName(sp.name);
//------------- line 2 ---------------------------
getCKLine(s, comment);
if (s.size() < 79) {
throw CK_SyntaxError(*m_log,
"Size of second line is too small", m_line);
}
// Read the number of temperature regions.
string sN = s.substr(0,2);
sp.nTempRegions = de_atoi(*m_log, sN);
string refDataCode = s.substr(3,6);
// elemental composition (first 5 elements)
for (int i = 0; i < 5; i++) {
string elementSym = "";
int iloc = 10 + 8*i;
if (s[iloc] != ' ') {
if (s[iloc+1] != ' ') {
elementSym = s.substr(iloc,2);
} else {
elementSym = s.substr(iloc,1);
}
} else if (s[iloc+1] != ' ') {
elementSym = s.substr(iloc+1,1);
}
double atoms = de_atof(s.substr(iloc+2,6));
addElement(elementSym, atoms, sp, *m_log);
}
// single-character phase descriptor
sp.phase = s.substr(50,2);
// Molecular weight in gm per gmol
string molecWeight = s.substr(52, 13);
// Heat of formation at 298.15 K in J / gmol
string Hf298_Jgmol = s.substr(65, 15);
vector_fp *coeffs_ptr;
for (int i = 0; i < sp.nTempRegions; i++) {
coeffs_ptr = new vector_fp(9);
vector_fp &coeffs = *coeffs_ptr;
//------------- line 3 ---------------------------
getCKLine(s, comment);
if (s.size() < 79) {
throw CK_SyntaxError(*m_log,
"Size of third line is too small", m_line);
}
string sTlow = s.substr(0, 11);
double tLow = de_atof(sTlow);
string sTHigh = s.substr(11, 11);
double tHigh = de_atof(sTHigh);
string sNCoeff = s.substr(22, 1);
int nCoeff = de_atoi(*m_log, sNCoeff);
if (nCoeff != 7) {
throw CK_SyntaxError(*m_log, "ncoeff ne 7", m_line);
}
string sTCoeff1 = s.substr(24, 5);
double TCoeff1 = de_atof(sTCoeff1);
if (TCoeff1 != -2.0) {
throw CK_SyntaxError(*m_log, "TCoeff1 ne -2.0", m_line);
}
string sTCoeff2 = s.substr(29, 5);
double TCoeff2 = de_atof(sTCoeff2);
if (TCoeff2 != -1.0) {
throw CK_SyntaxError(*m_log, "TCoeff2 ne -1.0", m_line);
}
string sTCoeff3 = s.substr(34, 5);
double TCoeff3 = de_atof(sTCoeff3);
if (TCoeff3 != 0.0) {
throw CK_SyntaxError(*m_log, "TCoeff3 ne 0.0", m_line);
}
string sTCoeff4 = s.substr(39, 5);
double TCoeff4 = de_atof(sTCoeff4);
if (TCoeff4 != 1.0) {
throw CK_SyntaxError(*m_log, "TCoeff4 ne 1.0", m_line);
}
string sTCoeff5 = s.substr(44, 5);
double TCoeff5 = de_atof(sTCoeff5);
if (TCoeff5 != 2.0) {
throw CK_SyntaxError(*m_log, "TCoeff5 ne 2.0", m_line);
}
string sTCoeff6 = s.substr(49, 5);
double TCoeff6 = de_atof(sTCoeff6);
if (TCoeff6 != 3.0) {
throw CK_SyntaxError(*m_log, "TCoeff6 ne 3.0", m_line);
}
string sTCoeff7 = s.substr(54, 5);
double TCoeff7 = de_atof(sTCoeff7);
if (TCoeff7 != 4.0) {
throw CK_SyntaxError(*m_log, "TCoeff7 ne 4.0", m_line);
}
string sHf298mHF0 = s.substr(65, 15);
//------------- line 4 ---------------------------
getCKLine(s, comment);
if (s.size() < 79) {
throw CK_SyntaxError(*m_log,
"Size of third line is too small", m_line);
}
numstr = s.substr(0, 16);
cf = getNumberFromString(numstr);
if (cf == UNDEF) illegalNumber(*m_log, numstr, m_line);
coeffs[0] = cf;
numstr = s.substr(16, 16);
cf = getNumberFromString(numstr);
if (cf == UNDEF) illegalNumber(*m_log, numstr, m_line);
coeffs[1] = cf;
numstr = s.substr(32, 16);
cf = getNumberFromString(numstr);
if (cf == UNDEF) illegalNumber(*m_log, numstr, m_line);
coeffs[2] = cf;
numstr = s.substr(48, 16);
cf = getNumberFromString(numstr);
if (cf == UNDEF) illegalNumber(*m_log, numstr, m_line);
coeffs[3] = cf;
numstr = s.substr(64, 16);
cf = getNumberFromString(numstr);
if (cf == UNDEF) illegalNumber(*m_log, numstr, m_line);
coeffs[4] = cf;
//------------- line 5 ---------------------------
getCKLine(s, comment);
if (s.size() < 79) {
throw CK_SyntaxError(*m_log,
"Size of fourth line is too small", m_line);
}
numstr = s.substr(0, 16);
cf = getNumberFromString(numstr);
if (cf == UNDEF) illegalNumber(*m_log, numstr, m_line);
coeffs[5] = cf;
numstr = s.substr(16, 16);
cf = getNumberFromString(numstr);
if (cf == UNDEF) illegalNumber(*m_log, numstr, m_line);
coeffs[6] = cf;
numstr = s.substr(48, 16);
cf = getNumberFromString(numstr);
if (cf == UNDEF) illegalNumber(*m_log, numstr, m_line);
coeffs[7] = cf;
numstr = s.substr(64, 16);
cf = getNumberFromString(numstr);
if (cf == UNDEF) illegalNumber(*m_log, numstr, m_line);
coeffs[8] = cf;
// Store the coefficients.
sp.minTemps.push_back(tLow);
sp.maxTemps.push_back(tHigh);
sp.region_coeffs.push_back(coeffs_ptr);
}
sp.valid = 1;
}
} // ckr namespace

2
Cantera/src/converters/Reaction.cpp
View file

@ -110,7 +110,7 @@ namespace ckr {
* not participate in the reaction.
*/
double Reaction::stoichCoefficient(const string& s) const {
double Reaction::stoichCoefficient(const std::string& s) const {
int k;
int nr = static_cast<int>(reactants.size());
for (k = 0; k < nr; k++)

101
Cantera/src/converters/Species.cpp Normal file
View file

@ -0,0 +1,101 @@
/**
* @file Species.cpp
*
*/
// Copyright 2001 California Institute of Technology
#include "Species.h"
#include <iostream>
#include <stdio.h>
namespace ckr {
// Construct an empty Species object
Species::Species() :
thermoFormatType(0),
name ("<empty>"),
id ("<none>"),
phase (""),
tlow(0.0),
tmid(0.0),
thigh(0.0),
nTempRegions(2),
valid(0),
index(-1)
{
}
// Destructor
Species::~Species() {
delR();
}
void Species::delR() {
int iReg = region_coeffs.size();
for (int i = 0; i < iReg; i++) {
if (region_coeffs[i]) {
delete region_coeffs[i];
region_coeffs[i] = 0;
}
}
}
//! Copy constructor
Species::Species(const Species& s)
{
/*
* Use the assignment operator to do the brunt
* of the work for the copy construtor.
*/
*this = s;
}
// Assignment operator
Species& Species::operator=(const Species& s) {
if (&s == this) return *this;
thermoFormatType = s.thermoFormatType;
name = s.name;
id = s.id;
phase = s.phase;
tlow = s.tlow;
tmid = s.tmid;
thigh = s.thigh;
nTempRegions = s.nTempRegions;
elements = s.elements;
comp = s.comp;
lowCoeffs = s.lowCoeffs;
highCoeffs = s.highCoeffs;
delR();
int iReg = s.region_coeffs.size();
for (int i = 0; i < iReg; i++) {
region_coeffs.push_back(new vector_fp(*(s.region_coeffs[i])));
}
minTemps = s.minTemps;
maxTemps = s.maxTemps;
m_commentsRef = s.m_commentsRef;
valid = s.valid;
index = s.index;
return *this;
}
// Test for equality based on name only.
bool Species::operator==(const Species& s) const {
return (s.name == name);
}
bool Species::operator!=(const Species& s) const {
return !(*this == s);
}
// Used to sort lists of species by index number.
bool Species::operator<(const Species& s) const {
return (index < s.index);
}
}

96
Cantera/src/converters/Species.h
View file

@ -2,7 +2,9 @@
* @file Species.h
*
*/
/*
* $Id$
*/
// Copyright 2001 California Institute of Technology
@ -16,72 +18,54 @@
//#include "Cantera.h"
using namespace std;
namespace ckr {
/**
* Holds species data read in from entries in the THERMO section of
* the input file.
*/
class Species {
public:
/**
* Holds species data read in from entries in the THERMO section of
* a chemkin or nasa9 fortran formatted input file.
*/
class Species {
public:
/// Construct an empty Species object
Species() :
name ("<empty>"),
id ("<none>"),
phase (""),
tlow(0.0),
tmid(0.0),
thigh(0.0),
valid(0),
index(-1)
{}
Species();
//! Copy constructor
Species(const Species& s);
/// Destructor
~Species() {}
~Species();
/// Assignment operator
Species& operator=(const Species& s) {
if (&s == this) return *this;
name = s.name;
id = s.id;
phase = s.phase;
tlow = s.tlow;
tmid = s.tmid;
thigh = s.thigh;
elements = s.elements;
comp = s.comp;
lowCoeffs = s.lowCoeffs;
highCoeffs = s.highCoeffs;
valid = s.valid;
index = s.index;
return *this;
}
Species& operator=(const Species& s);
/// Test for equality based on name only.
bool operator==(const Species& s) const {
return (s.name == name);
}
bool operator==(const Species& s) const;
bool operator!=(const Species& s) const {
return !(*this == s);
}
bool operator!=(const Species& s) const;
/// Used to sort lists of species by index number.
bool operator<(const Species& s) const {
return (index < s.index);
}
bool operator<(const Species& s) const;
//! Type of thermodynamic representation
/*!
* 0 This is a 2 region NASA polynomial representation
*
* 1 This is a multiple temperature region NASA9 polynomial
* representation.
*/
int thermoFormatType;
string name; //!< Species name
//! Species Name
string name;
string id; //!< ID tag from 'date' field in input
string phase; //!< Phase string. Usually "G", "L", or "S".
double tlow; //!< Min temperature for thermo data fit
double tmid; //!< Mid temperature for thermo data fit
double thigh; //!< Max temperature for thermo data fit
/// list of Constituent objects defining elemental composition
vector<Constituent> elements;
@ -94,19 +78,31 @@ public:
/// polynomial coefficients for the upper temperature range
vector_fp highCoeffs;
//! Number of temperature regions
int nTempRegions;
std::vector<vector_fp *> region_coeffs;
vector_fp minTemps;
vector_fp maxTemps;
/// flag set by the validation routines
int valid;
/// position in the list of species in the input file
int index;
};
/// A list of Species
typedef vector<Species> speciesList;
string m_commentsRef;
/// A map from species names to Species objects
typedef map<string, Species> speciesTable;
private:
//! Delete private data
void delR();
};
//! Shorthand for a list of Species
typedef vector<Species> speciesList;
//! A map from species names to Species objects
typedef map<string, Species> speciesTable;
}
#endif

128
Cantera/src/converters/ck2ct.cpp
View file

@ -1,7 +1,9 @@
/**
* @file ck2ct.cpp
*
* Convert CK-format reaction mechanism files to Cantera input format.
* Convert CK-format reaction mechanism files to Cantera input format.
*/
/*
* $Id$
*
*/
#ifdef WIN32
@ -115,23 +117,67 @@ namespace pip {
}
// add a NASA polynomial parameterization
static void addNASA(FILE* f,
const vector_fp& low, const vector_fp& high,
doublereal minx, doublereal midx,
doublereal maxx) {
// add a NASA polynomial parameterization
static void addNASA(FILE* f,
const vector_fp& low, const vector_fp& high,
doublereal minx, doublereal midx,
doublereal maxx) {
fprintf(f," thermo = (\n");
fprintf(f," NASA( [%8.2f, %8.2f], ", minx, midx);
fprintf(f,"[%17.9E, %17.9E, \n", low[0], low[1]);
fprintf(f," %17.9E, %17.9E, %17.9E,\n", low[2], low[3], low[4]);
fprintf(f," %17.9E, %17.9E] ),\n", low[5], low[6]);
fprintf(f," NASA( [%8.2f, %8.2f], ", midx, maxx);
fprintf(f,"[%17.9E, %17.9E, \n", high[0], high[1]);
fprintf(f," %17.9E, %17.9E, %17.9E,\n", high[2], high[3], high[4]);
fprintf(f," %17.9E, %17.9E] )\n", high[5], high[6]);
fprintf(f," )");
fprintf(f," thermo = (\n");
fprintf(f," NASA( [%8.2f, %8.2f], ", minx, midx);
fprintf(f,"[%17.9E, %17.9E, \n", low[0], low[1]);
fprintf(f," %17.9E, %17.9E, %17.9E,\n",
low[2], low[3], low[4]);
fprintf(f," %17.9E, %17.9E] ),\n", low[5], low[6]);
fprintf(f," NASA( [%8.2f, %8.2f], ", midx, maxx);
fprintf(f,"[%17.9E, %17.9E, \n", high[0], high[1]);
fprintf(f," %17.9E, %17.9E, %17.9E,\n",
high[2], high[3], high[4]);
fprintf(f," %17.9E, %17.9E] )\n", high[5], high[6]);
fprintf(f," )");
}
//! Add a NASA polynomial parameterization to the cti file
/*!
* This little tidbit of code writes out the polynomials to the cti file
*/
static void addNASA9(FILE* f,
const std::vector<vector_fp *> &region_coeffs,
const vector_fp &minTemps, const vector_fp &maxTemps)
{
int nReg = region_coeffs.size();
if ((int) minTemps.size() != nReg) {
throw CanteraError("addNASA9", "incompat");
}
if ((int) maxTemps.size() != nReg) {
throw CanteraError("addNASA9", "incompat");
}
fprintf(f," thermo = (\n");
for (int i = 0; i < nReg; i++) {
double minT = minTemps[i];
double maxT = maxTemps[i];
const vector_fp &coeffs = *(region_coeffs[i]);
if ((int) coeffs.size() != 9) {
throw CanteraError("addNASA9", "incompat");
}
fprintf(f," NASA9( [%8.2f, %8.2f], ", minT, maxT);
fprintf(f,"[%17.9E, %17.9E, %17.9E,\n", coeffs[0],
coeffs[1], coeffs[2]);
fprintf(f," %17.9E, %17.9E, %17.9E,\n",
coeffs[3], coeffs[4], coeffs[5]);
fprintf(f," %17.9E, %17.9E, %17.9E] )",
coeffs[6], coeffs[7], coeffs[8]);
if (i < nReg - 1) {
fprintf(f,",\n");
} else {
fprintf(f,"\n");
}
}
fprintf(f," )");
}
static void addTransportParams(FILE* f, string name) {
@ -230,24 +276,42 @@ namespace pip {
}
fprintf(f," atoms = \"%s\",\n", str.c_str());
if (sp.lowCoeffs.size() == 0) {
throw CanteraError("addSpecies",
"Low Nasa Thermo Polynomial was not found");
}
if (sp.highCoeffs.size() == 0) {
throw CanteraError("addSpecies",
"High Nasa Thermo Polynomial was not found");
}
if (sp.tlow >= sp.thigh) {
// Add the NASA block according to the thermoFormatType value
if (sp.thermoFormatType == 0) {
if (sp.lowCoeffs.size() == 0) {
throw CanteraError("addSpecies",
"Low Nasa Thermo Polynomial was not found");
}
if (sp.highCoeffs.size() == 0) {
throw CanteraError("addSpecies",
"High Nasa Thermo Polynomial was not found");
}
if (sp.tlow >= sp.thigh) {
throw CanteraError("addSpecies",
"Low temp limit is greater or equal to high temp limit");
}
addNASA(f, sp.lowCoeffs, sp.highCoeffs,
sp.tlow, sp.tmid, sp.thigh);
} else if (sp.thermoFormatType == 1) {
// This new typs is a multiregion 9 coefficient formulation
addNASA9(f, sp.region_coeffs, sp.minTemps, sp.maxTemps);
} else {
throw CanteraError("addSpecies", "Unknown thermoFormatType");
}
addNASA(f, sp.lowCoeffs, sp.highCoeffs,
sp.tlow, sp.tmid, sp.thigh);
if (_with_transport)
addTransportParams(f, sp.name);
if (sp.id != "") fprintf(f,",\n note = \"%s\"", sp.id.c_str());
if (_with_transport) {
addTransportParams(f, sp.name);
}
if (sp.id != "" || sp.m_commentsRef != "") {
fprintf(f,",\n note = \"");
if (sp.id != "") {
fprintf(f, "%s", sp.id.c_str());
}
if (sp.m_commentsRef != "") {
fprintf(f, " %s", sp.m_commentsRef.c_str());
}
fprintf(f, "\"");
}
fprintf(f,"\n )\n");
}
@ -419,7 +483,7 @@ namespace pip {
writeline(f);
for (i = 0; i < nsp; i++) {
addSpecies(f, idtag, r.species[i]);
addSpecies(f, idtag, r.species[i]);
}
fprintf(f, "\n\n\n");

10
Cantera/src/converters/ckr_defs.h
View file

@ -14,15 +14,15 @@
#include <string>
#include <iostream>
#include <vector>
using namespace std;
//using namespace std;
//#include "../ctvector.h"
/// the namespace for the CKReader packaage
namespace ckr {
typedef vector<double> vector_fp;
typedef vector<double> vector_int;
typedef std::vector<double> vector_fp;
typedef std::vector<double> vector_int;
//typedef vector<double> vector_fp;
// exceptions
@ -32,11 +32,11 @@ namespace ckr {
m_msg = "";
}
virtual ~CK_Exception() {}
string errorMessage() {
std::string errorMessage() {
return m_msg;
}
protected:
string m_msg;
std::string m_msg;
};
const double UNDEF = -9999.1234;

162
Cantera/src/converters/ckr_utils.cpp
View file

@ -16,131 +16,133 @@
#include <math.h>
#include "ckr_utils.h"
using namespace std;
namespace ckr {
bool match(const string& s1, const string& s2)
{
bool match(const std::string& s1, const std::string& s2)
{
size_t n = s2.size();
if (s1.size() < n) return false;
for (size_t i = 0; i < n; i++)
if (s2[i] != '*' && (toupper(s1[i]) != toupper(s2[i]))) return false;
if (s2[i] != '*' && (toupper(s1[i]) != toupper(s2[i]))) return false;
return true;
}
}
void removeWhiteSpace(string& s) {
void removeWhiteSpace(std::string& s) {
string r;
int ssize = static_cast<int>(s.size());
for (int n = 0; n < ssize; n++) if (s[n] != ' '
&& s[n] != '\t' && s[n] != '\n') r += s[n];
&& s[n] != '\t' && s[n] != '\n') r += s[n];
s = r;
}
}
/**
* Get tokens from char array of length n beginning at 'begin'. Tokens are
* delimited by character 'delim' (space by default), with the exception of
* text contained within forward slashes ('/'), which is treated literally.
*
* code
* vector<string> tokens;
* char line[] = "a b/3.0 txt/ c /4.0 6/ d";
* int n = strlen(line);
* getTokens(line, n, tokens);
* for (int i = 0; i < tokens.size(); i++) cout << tokens[i] << endl;
* endcode
*
*/
/**
* Get tokens from char array of length n beginning at 'begin'. Tokens are
* delimited by character 'delim' (space by default), with the exception of
* text contained within forward slashes ('/'), which is treated literally.
*
* code
* vector<string> tokens;
* char line[] = "a b/3.0 txt/ c /4.0 6/ d";
* int n = strlen(line);
* getTokens(line, n, tokens);
* for (int i = 0; i < tokens.size(); i++) cout << tokens[i] << endl;
* endcode
*
*/
void getTokens(string& s, int n, vector<string>& toks, char delim) {
void getTokens(std::string& s, int n, std::vector<std::string>& toks, char delim) {
string::iterator q, p = s.begin(), end = p + n;
vector<string> tokk;
int inslash = -1;
//p = begin;
while (1 > 0) {
for (; p < end; p++) {
if (*p != delim) break;
}
q = p;
for (; q < end; q++) {
if (*q == '/') inslash *= -1;
if (inslash < 0 && *q == delim) break;
}
if (p != q) {
tokk.push_back(s.substr(p - s.begin(), q - p));
}
p = q;
if (p == end) break;
for (; p < end; p++) {
if (*p != delim) break;
}
q = p;
for (; q < end; q++) {
if (*q == '/') inslash *= -1;
if (inslash < 0 && *q == delim) break;
}
if (p != q) {
tokk.push_back(s.substr(p - s.begin(), q - p));
}
p = q;
if (p == end) break;
}
toks.clear();
int nt = static_cast<int>(tokk.size());
string t = "";
for (int i = 0; i < nt; i++) {
if (tokk[i][0] == '/') t += tokk[i];
else {
if (t != "") toks.push_back(t);
t = tokk[i];
}
if (tokk[i][0] == '/') t += tokk[i];
else {
if (t != "") toks.push_back(t);
t = tokk[i];
}
}
if (t != "") toks.push_back(t);
}
}
/**
* Look for a slash-delimited number in string s,
* and if found set v to the numerical value, and
* set s to the portion of the string before the
* first slash. Return true if slash data found,
* false otherwise.
*/
bool extractSlashData(string& s, string& name, string& data) {
/**
* Look for a slash-delimited number in string s,
* and if found set v to the numerical value, and
* set s to the portion of the string before the
* first slash. Return true if slash data found,
* false otherwise.
*/
bool extractSlashData(std::string& s, std::string& name, std::string& data) {
int slen = static_cast<int>(s.size());
string::size_type n = s.find_first_of("/");
if (n != string::npos && (static_cast<int>(n) < slen)) {
int m;
string::size_type n = s.find_first_of("/");
if (n != string::npos && (static_cast<int>(n) < slen)) {
int m;
for (m = static_cast<int>(n)+1; m < slen; m++) if (s[m] == '/') break;
if (m < slen) {
data = s.substr(n+1,m-n-1);
name = s.substr(0,n);
removeWhiteSpace(name);
s = s.substr(m+1,1000);
return true;
}
else {
name = s;
removeWhiteSpace(name);
data = "";
s = "";
return false;
}
}
else {
if (m < slen) {
data = s.substr(n+1,m-n-1);
name = s.substr(0,n);
removeWhiteSpace(name);
s = s.substr(m+1,1000);
return true;
}
else {
name = s;
removeWhiteSpace(name);
data = "";
s = "";
return false;
}
}
}
else {
name = s;
removeWhiteSpace(name);
data = "";
s = "";
return false;
}
}
/**
* Return a modified version of string word, in which
* the first letter is upper case, and the rest are
* lower case.
*/
string capitalize(const string& word) {
/**
* Return a modified version of string word, in which
* the first letter is upper case, and the rest are
* lower case.
*/
string capitalize(const std::string& word) {
string cap = word;
int n = static_cast<int>(word.size());
if (n > 0) {
cap[0] = toupper(word[0]);
for (int m = 1; m < n; m++) {
cap[m] = tolower(word[m]);
}
cap[0] = toupper(word[0]);
for (int m = 1; m < n; m++) {
cap[m] = tolower(word[m]);
}
}
return cap;
}
}
}

112
Cantera/src/converters/thermoFunctions.cpp
View file

@ -1,6 +1,7 @@
/**
* @file thermoFunctions.cpp
*
* File containing thermo evalulation functions for NASA polynomials,
* which are used in testing the interpolations.
*/
// Copyright 2001 California Institute of Technology
@ -20,66 +21,117 @@ using namespace std;
namespace ckr {
/**
* non-dimensional heat capacity (\f$ C_p/R \f$) at constant P for
* one species @param t temperature @param s species object
*/
double cp(double t, const Species& s) {
/**
* non-dimensional heat capacity (\f$ C_p/R \f$) at constant P for
* one species @param t temperature @param s species object
*/
double cp(double t, const Species& s) {
if (s.thermoFormatType == 1) {
const vector_fp* cpc;
int ireg = -1;
for (int i = 0; i < s.nTempRegions; i++) {
if (t <= s.maxTemps[i]) {
ireg = i;
break;
}
}
cpc = s.region_coeffs[ireg];
const vector_fp& c = *cpc;
double cp0r = c[0]/(t*t) + c[1]/t + c[2] + c[3]*t + c[4]*t*t
+ c[5]*t*t*t + c[6]*t*t*t*t;
return cp0r;
}
const vector_fp* cpc;
if (t > s.tmid) cpc = &s.highCoeffs;
else cpc = &s.lowCoeffs;
const vector_fp& c = *cpc;
double cp0r = c[0] + c[1]*t + c[2]*t*t + c[3]*t*t*t + c[4]*t*t*t*t;
return cp0r;
}
}
/**
* enthalpy in Kelvin (\f$ H/R \f$) for
* one species. @param t temperature @param s species object
*/
double enthalpy(double t, const Species& s) {
/**
* enthalpy in Kelvin (\f$ H/R \f$) for
* one species. @param t temperature @param s species object
*/
double enthalpy(double t, const Species& s) {
if (s.thermoFormatType == 1) {
const vector_fp* cpc;
int ireg = -1;
for (int i = 0; i < s.nTempRegions; i++) {
if (t <= s.maxTemps[i]) {
ireg = i;
break;
}
}
cpc = s.region_coeffs[ireg];
const vector_fp& c = *cpc;
double h0rt = -c[0]/(t*t) + c[1]*log(t)/t
+ c[2] + 0.5*c[3]*t + c[4]*t*t/3.0 + 0.25*c[5]*t*t*t
+ 0.2*c[6]*t*t*t*t + c[7]/t;
return t*h0rt;
}
const vector_fp* cp;
if (t > s.tmid) cp = &s.highCoeffs;
else cp = &s.lowCoeffs;
const vector_fp& c = *cp;
double h0rt = c[0] + 0.5*c[1]*t + c[2]*t*t/3.0 + 0.25*c[3]*t*t*t
+ 0.2*c[4]*t*t*t*t + c[5]/t;
+ 0.2*c[4]*t*t*t*t + c[5]/t;
return t*h0rt;
}
}
/**
* non-dimensional entropy (\f$ S/R \f$) for
* one species @param t temperature @param s species object
*/
double entropy(double t, const Species& s) {
/**
* non-dimensional entropy (\f$ S/R \f$) for
* one species @param t temperature @param s species object
*/
double entropy(double t, const Species& s) {
if (s.thermoFormatType == 1) {
const vector_fp* cpc;
int ireg = -1;
for (int i = 0; i < s.nTempRegions; i++) {
if (t <= s.maxTemps[i]) {
ireg = i;
break;
}
}
cpc = s.region_coeffs[ireg];
const vector_fp& c = *cpc;
double s0r = -0.5*c[0]/(t*t) - c[1]/t
+ c[2]*log(t) + c[3]*t + 0.5*c[4]*t*t + c[5]*t*t*t/3.0
+ 0.25*c[6]*t*t*t*t + c[8];
return t*s0r;
}
const vector_fp* cp;
if (t > s.tmid) cp = &s.highCoeffs;
else cp = &s.lowCoeffs;
const vector_fp& c = *cp;
double s0r = c[0]*log(t) + c[1]*t + 0.5*c[2]*t*t + c[3]*t*t*t/3.0
+ 0.25*c[4]*t*t*t*t + c[6];
+ 0.25*c[4]*t*t*t*t + c[6];
return t*s0r;
}
}
/**
* Gibbs function in Kelvin (\f$ G/R \f$) for
* one species. @param t temperature @param s species object
*/
double gibbs(double t, const Species& s) {
/**
* Gibbs function in Kelvin (\f$ G/R \f$) for
* one species. @param t temperature @param s species object
*/
double gibbs(double t, const Species& s) {
if (s.thermoFormatType == 1) {
double s0r = entropy(t, s);
double h0r = enthalpy(t, s);
return (h0r - s0r * t);
}
const vector_fp* cp;
if (t > s.tmid) cp = &s.highCoeffs;
else cp = &s.lowCoeffs;
const vector_fp& c = *cp;
double h0rt = c[0] + 0.5*c[1]*t + c[2]*t*t/3.0 + 0.25*c[3]*t*t*t
+ 0.2*c[4]*t*t*t*t + c[5]/t;
+ 0.2*c[4]*t*t*t*t + c[5]/t;
double s0r = c[0]*log(t) + c[1]*t + 0.5*c[2]*t*t + c[3]*t*t*t/3.0
+ 0.25*c[4]*t*t*t*t + c[6];
+ 0.25*c[4]*t*t*t*t + c[6];
return t*(h0rt - s0r);
}
}
}

78
Cantera/src/converters/writelog.cpp
View file

@ -14,11 +14,13 @@
#include <stdio.h>
#include "writelog.h"
using namespace std;
namespace ckr {
/// format a string for the log file message printed when starting a new task
string newTask(string msg) {
string newTask(std::string msg) {
string s = "\n";
s += msg + "...\n";
return s;
@ -26,7 +28,7 @@ namespace ckr {
/// print the falloff parameters for a pressure-dependent reaction
bool writeFalloff(int type, const vector_fp& c, ostream& log) {
bool writeFalloff(int type, const vector_fp& c, std::ostream& log) {
log.precision(6);
log.width(0);
@ -86,7 +88,7 @@ namespace ckr {
/// print the rate coefficient parameters
bool writeRateCoeff(const RateCoeff& k, ostream& log) {
bool writeRateCoeff(const RateCoeff& k, std::ostream& log) {
log.precision(10);
log.width(0);
@ -143,7 +145,7 @@ namespace ckr {
/**
* Write onto an output stream the chemical equation for a reaction.
*/
void printReactionEquation(ostream& f, const Reaction& r) {
void printReactionEquation(std::ostream& f, const Reaction& r) {
// r.write(f);
f << reactionEquation(r);
}
@ -200,7 +202,7 @@ namespace ckr {
* @param spec instance of Species class
*/
void writeSpeciesData(ostream& log, const Species& spec) {
void writeSpeciesData(std::ostream& log, const Species& spec) {
if (!spec.id.empty())
log << endl << " id/date: " << spec.id << endl;
@ -212,32 +214,54 @@ namespace ckr {
<< " composition: (";
for (size_t ie = 0; ie < spec.elements.size(); ie++) {
if (!spec.elements[ie].name.empty()) {
log.flags(ios::fixed);
log.precision(0);
if (ie > 0) log << ", ";
log << spec.elements[ie].number << " "
<< spec.elements[ie].name;
}
if (!spec.elements[ie].name.empty()) {
log.flags(ios::fixed);
log.precision(0);
if (ie > 0) log << ", ";
log << spec.elements[ie].number << " "
<< spec.elements[ie].name;
}
}
log << ")";
log.flags(ios::showpoint | ios::fixed);
log.precision(2);
log << endl << " Tlow, Tmid, Thigh: (" << spec.tlow << ", " <<
if (spec.thermoFormatType == 0) {
log.flags(ios::showpoint | ios::fixed);
log.precision(2);
log << endl << " Tlow, Tmid, Thigh: (" << spec.tlow << ", " <<
spec.tmid << ", " << spec.thigh << ")" << endl << endl;
log << " coefficients (low, high):" << endl;
log.flags(ios::scientific | ios::uppercase | ios::internal);
log.precision(8);
for (int j = 0; j < 7; j++) {
log << " a" << j + 1;
log.setf(ios::showpos);
log << " \t" << spec.lowCoeffs[j]
<< " \t" << spec.highCoeffs[j] << endl;
log.unsetf(ios::showpos);
log << " coefficients (low, high):" << endl;
log.flags(ios::scientific | ios::uppercase | ios::internal );
log.precision(8);
for (int j = 0; j < 7; j++) {
log << " a" << j + 1;
log.setf(ios::showpos);
log << " \t" << spec.lowCoeffs[j]
<< " \t" << spec.highCoeffs[j] << endl;
log.unsetf(ios::showpos);
}
log << endl;
} else if (spec.thermoFormatType == 1) {
log.flags(ios::showpoint | ios::fixed);
log.precision(2);
log << endl;
log << "Number of temp regions = " << spec.nTempRegions << endl;
for (int i = 0; i < spec.nTempRegions; i++) {
log << " Tlow, Thigh: (" << spec.minTemps[i] << ", "
<< spec.maxTemps[i] << ")" << endl << endl;
log << " coefficients :" << endl;
log.flags( ios::scientific | ios::uppercase | ios::internal);
log.precision(8);
vector_fp &cc = *spec.region_coeffs[i];
for (int j = 0; j < 9; j++) {
log << " a" << j + 1;
log.setf(ios::showpos);
log << " \t" << cc[j] << endl;
log.unsetf(ios::showpos);
}
log << endl;
}
}
log << endl;
}
}
}

14
Cantera/src/converters/writelog.h
View file

@ -12,7 +12,7 @@
#include <string>
#include <vector>
#include <iostream>
using namespace std;
//using namespace std;
#include "Species.h"
#include "Reaction.h"
@ -20,12 +20,12 @@ using namespace std;
//#include "Cantera.h"
namespace ckr {
string newTask(string msg);
bool writeFalloff(int type, const vector_fp& c, ostream& log);
bool writeRateCoeff(const RateCoeff& k, ostream& log);
void printReactionEquation(ostream& f, const Reaction& r);
void writeSpeciesData(ostream& log, const Species& spec);
string reactionEquation(const Reaction& r);
std::string newTask(string msg);
bool writeFalloff(int type, const vector_fp& c, std::ostream& log);
bool writeRateCoeff(const RateCoeff& k, std::ostream& log);
void printReactionEquation(std::ostream& f, const Reaction& r);
void writeSpeciesData(std::ostream& log, const Species& spec);
std::string reactionEquation(const Reaction& r);
}
#endif

3
Cantera/src/thermo/ConstCpPoly.h
View file

@ -101,6 +101,9 @@ namespace Cantera {
//! Returns an integer representing the type of parameterization
virtual int reportType() const { return CONSTANT_CP; }
//! Returns an integer representing the species index
virtual int speciesIndex() const { return m_index; }
//! Update the properties for this species, given a temperature polynomial
/*!
* This method is called with a pointer to an array containing the functions of

24
Cantera/src/thermo/GeneralSpeciesThermo.cpp
View file

@ -158,6 +158,30 @@ namespace Cantera {
m_thigh_min = min(maxTemp, m_thigh_min);
}
// Install a new species thermodynamic property
// parameterization for one species.
/*
* @param stit_ptr Pointer to the SpeciesThermoInterpType object
* This will set up the thermo for one species
*/
void GeneralSpeciesThermo::install_STIT(SpeciesThermoInterpType *stit_ptr) {
/*
* Resize the arrays if necessary, filling the empty
* slots with the zero pointer.
*/
int index = stit_ptr->speciesIndex();
if (index > m_kk - 1) {
m_sp.resize(index+1, 0);
m_kk = index+1;
}
AssertThrow(m_sp[index] == 0,
"Index position isn't null, duplication of assignment: " + int2str(index));
/*
* Now, simply assign the position
*/
m_sp[index] = stit_ptr;
}
/**
* Update the properties for one species.
*/

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

@ -19,6 +19,7 @@
#include "ct_defs.h"
#include "SpeciesThermoMgr.h"
#include "NasaPoly1.h"
#include "Nasa9Poly1.h"
#include "speciesThermoTypes.h"
//#include "polyfit.h"
@ -88,6 +89,14 @@ namespace Cantera {
doublereal minTemp, doublereal maxTemp,
doublereal refPressure);
//! Install a new species thermodynamic property
//! parameterization for one species.
/*!
* @param stit_ptr Pointer to the SpeciesThermoInterpType object
* This will set up the thermo for one species
*/
virtual void install_STIT(SpeciesThermoInterpType *stit_ptr);
//! Like update(), but only updates the single species k.
/*!
* @param k species index

4
Cantera/src/thermo/Makefile.in
View file

@ -32,14 +32,14 @@ CXX_FLAGS = @CXXFLAGS@ $(LOCAL_DEFS) $(CXX_OPT) $(PIC_FLAG) $(DEBUG_FLAG)
# Basic Cantera Thermodynamics Object Files
THERMO_OBJ = State.o Elements.o Constituents.o Phase.o \
ThermoPhase.o IdealGasPhase.o ConstDensityThermo.o \
SpeciesThermoFactory.o ConstCpPoly.o \
SpeciesThermoFactory.o ConstCpPoly.o Nasa9Poly1.o Nasa9PolyMultiTempRegion.o \
Mu0Poly.o GeneralSpeciesThermo.o SurfPhase.o \
ThermoFactory.o phasereport.o @phase_object_files@
THERMO_H = State.h Elements.h Constituents.h Phase.h mix_defs.h \
ThermoPhase.h IdealGasPhase.h ConstDensityThermo.h \
SpeciesThermoFactory.h ThermoFactory.h \
NasaPoly1.h NasaPoly2.h NasaThermo.h \
NasaPoly1.h NasaPoly2.h NasaThermo.h Nasa9Poly1.h Nasa9PolyMultiTempRegion.h \
ShomateThermo.h ShomatePoly.h ConstCpPoly.h \
SimpleThermo.h SpeciesThermoMgr.h \
SpeciesThermoInterpType.h \

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

@ -139,6 +139,8 @@ namespace Cantera {
//! Returns an integer representing the type of parameterization
virtual int reportType() const { return MU0_INTERP; }
//! Returns an integer representing the species index
virtual int speciesIndex() const { return m_index; }
//! Update the properties for this species, given a temperature polynomial
/*!

298
Cantera/src/thermo/Nasa9Poly1.cpp Normal file
View file

@ -0,0 +1,298 @@
/**
* @file Nasa9Poly1.h
* Header for a single-species standard state object derived
* from
* \link Cantera::SpeciesThermoInterpType SpeciesThermoInterpType\endlink
* based
* on the NASA 9 coefficient temperature polynomial form applied to one temperature region
* (see \ref spthermo and class \link Cantera::Nasa9Poly1 Nasa9Poly1\endlink).
*
* This parameterization has one NASA temperature region.
*/
/* $Author$
* $Revision$
* $Date$
*/
// Copyright 2007 Sandia National Laboratories
//#include "global.h"
#include "Nasa9Poly1.h"
namespace Cantera {
// The NASA 9 polynomial parameterization for one temperature range.
/*
* This parameterization expresses the heat capacity via a
* 7 coefficient polynomial.
*
* Note that this is the form used in the
* 2002 NASA equilibrium program
*
* NASA Glenn Coefficients for Calculating Thermodynamic
* Properties of Individual Species,
* B. J. McBride, M. J. Zehe, S. Gordon
* NASA/TP-2002-211556, Sept. 2002
*
*
* Nine coefficients \f$(a_0,\dots,a_6)\f$ are used to represent
* \f$ C_p^0(T)\f$, \f$ H^0(T)\f$, and \f$ S^0(T) \f$ as
* polynomials in \f$ T \f$ :
* \f[
* \frac{c_p(T)}{R} = a_0 T^{-2} + a_1 T^{-1} + a_2 + a_3 T
* + a_4 T^2 + a_5 T^3 + a_6 T^4
* \f]
*
* \f[
* \frac{H^0(T)}{RT} = - a_0 T^{-2} + a_1 \frac{\ln(T)}{T} + a_2
* + a_3 T + a_4 T^2 + a_5 T^3 + a_6 T^4 + \frac{a_7}{T}
* \f]
*
* \f[
* \frac{s^0(T)}{R} = - \frac{a_0}{2} T^{-2} - a_1 T^{-1} + a_2 \ln(T)
+ + a_3 T \frac{a_4}{2} T^2 + \frac{a_5}{3} T^3 + \frac{a_6}{4} T^4 + a_8
* \f]
*
* The standard state is assumed to be the ideal gas at the
* standard pressure of 1 bar, for gases.
* For condensed species, the standard state is the
* pure cyrstalline or liquid substance at the standard
* pressure of 1 atm.
*
* These NASA representations may have more than 2 temperature regions.
*
* @ingroup spthermo
*/
//! Empty constructor
Nasa9Poly1::Nasa9Poly1()
: m_lowT(0.0), m_highT (0.0),
m_Pref(1.0E5), m_index (0), m_coeff(array_fp(9)) {}
// constructor used in templated instantiations
/*
* @param n Species index
* @param tlow Minimum temperature
* @param thigh Maximum temperature
* @param pref reference pressure (Pa).
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state.
*/
Nasa9Poly1::Nasa9Poly1(int n, doublereal tlow, doublereal thigh,
doublereal pref,
const doublereal* coeffs) :
m_lowT (tlow),
m_highT (thigh),
m_Pref (pref),
m_index (n),
m_coeff (array_fp(9)) {
std::copy(coeffs, coeffs + 9, m_coeff.begin());
}
// copy constructor
/*
* @param b object to be copied
*/
Nasa9Poly1::Nasa9Poly1(const Nasa9Poly1& b) :
m_lowT (b.m_lowT),
m_highT (b.m_highT),
m_Pref (b.m_Pref),
m_index (b.m_index),
m_coeff (array_fp(9)) {
std::copy(b.m_coeff.begin(),
b.m_coeff.begin() + 9,
m_coeff.begin());
}
// assignment operator
/*
* @param b object to be copied
*/
Nasa9Poly1& Nasa9Poly1::operator=(const Nasa9Poly1& b) {
if (&b != this) {
m_lowT = b.m_lowT;
m_highT = b.m_highT;
m_Pref = b.m_Pref;
m_index = b.m_index;
std::copy(b.m_coeff.begin(),
b.m_coeff.begin() + 9,
m_coeff.begin());
}
return *this;
}
// Destructor
Nasa9Poly1::~Nasa9Poly1() {
}
// duplicator
SpeciesThermoInterpType *
Nasa9Poly1::duplMyselfAsSpeciesThermoInterpType() const {
Nasa9Poly1* np = new Nasa9Poly1(*this);
return (SpeciesThermoInterpType *) np;
}
// Returns the minimum temperature that the thermo
// parameterization is valid
doublereal Nasa9Poly1::minTemp() const { return m_lowT;}
// Returns the maximum temperature that the thermo
// parameterization is valid
doublereal Nasa9Poly1::maxTemp() const {
return m_highT;
}
// Returns the reference pressure (Pa)
doublereal Nasa9Poly1::refPressure() const { return m_Pref; }
// Returns an integer representing the type of parameterization
int Nasa9Poly1::reportType() const {
return NASA9;
}
// Returns an integer representing the species index
int Nasa9Poly1::speciesIndex() const {
return m_index;
}
// Update the properties for this species, given a temperature polynomial
/*
* This method is called with a pointer to an array containing the functions of
* temperature needed by this parameterization, and three pointers to arrays where the
* computed property values should be written. This method updates only one value in
* each array.
*
* Temperature Polynomial:
* tt[0] = t;
* tt[1] = t*t;
* tt[2] = t*t*t;
* tt[3] = t*t*t*t;
* tt[4] = 1.0/t;
* tt[5] = 1.0/(t*t);
* tt[6] = std::log(t);
*
* @param tt vector of temperature polynomials
* @param cp_R Vector of Dimensionless heat capacities.
* (length m_kk).
* @param h_RT Vector of Dimensionless enthalpies.
* (length m_kk).
* @param s_R Vector of Dimensionless entropies.
* (length m_kk).
*/
void Nasa9Poly1::updateProperties(const doublereal* tt,
doublereal* cp_R, doublereal* h_RT,
doublereal* s_R) const {
doublereal ct0 = m_coeff[0] * tt[5]; // a0 / (T^2)
doublereal ct1 = m_coeff[1] * tt[4]; // a1 / T
doublereal ct2 = m_coeff[2]; // a2
doublereal ct3 = m_coeff[3] * tt[0]; // a3 * T
doublereal ct4 = m_coeff[4] * tt[1]; // a4 * T^2
doublereal ct5 = m_coeff[5] * tt[2]; // a5 * T^3
doublereal ct6 = m_coeff[6] * tt[3]; // a6 * T^4
doublereal cpdivR = ct0 + ct1 + ct2 + ct3 + ct4 + ct5 + ct6;
doublereal hdivRT = -ct0 + tt[6]*ct1 + ct2 + 0.5*ct3 + OneThird*ct4
+ 0.25*ct5 + 0.2*ct6 + m_coeff[7] * tt[4];
doublereal sdivR = -2.0*ct0 - ct1 + tt[6]*ct2 + ct3 + 0.5*ct4
+ OneThird*ct5 + 0.25*ct6 + m_coeff[8];
// return the computed properties in the location in the output
// arrays for this species
cp_R[m_index] = cpdivR;
h_RT[m_index] = hdivRT;
s_R[m_index] = sdivR;
//writelog("NASA9poly1: for species "+int2str(m_index)+", h_RT = "+
// fp2str(h)+"\n");
}
// Compute the reference-state property of one species
/*
* Given temperature T in K, this method updates the values of
* the non-dimensional heat capacity at constant pressure,
* enthalpy, and entropy, at the reference pressure, Pref
* of one of the species. The species index is used
* to reference into the cp_R, h_RT, and s_R arrays.
*
* Temperature Polynomial:
* tt[0] = t;
* tt[1] = t*t;
* tt[2] = t*t*t;
* tt[3] = t*t*t*t;
* tt[4] = 1.0/t;
* tt[5] = 1.0/(t*t);
* tt[6] = std::log(t);
*
* @param temp Temperature (Kelvin)
* @param cp_R Vector of Dimensionless heat capacities.
* (length m_kk).
* @param h_RT Vector of Dimensionless enthalpies.
* (length m_kk).
* @param s_R Vector of Dimensionless entropies.
* (length m_kk).
*/
void Nasa9Poly1::updatePropertiesTemp(const doublereal temp,
doublereal* cp_R, doublereal* h_RT,
doublereal* s_R) const {
double tPoly[7];
tPoly[0] = temp;
tPoly[1] = temp * temp;
tPoly[2] = tPoly[1] * temp;
tPoly[3] = tPoly[2] * temp;
tPoly[4] = 1.0 / temp;
tPoly[5] = tPoly[4] / temp;
tPoly[6] = std::log(temp);
updateProperties(tPoly, cp_R, h_RT, s_R);
}
//This utility function reports back the type of
// parameterization and all of the parameters for the
// species, index.
/*
* All parameters are output variables
*
* @param n Species index
* @param type Integer type of the standard type
* @param tlow output - Minimum temperature
* @param thigh output - Maximum temperature
* @param pref output - reference pressure (Pa).
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state.
*/
void Nasa9Poly1::reportParameters(int &n, int &type,
doublereal &tlow, doublereal &thigh,
doublereal &pref,
doublereal* const coeffs) const {
n = m_index;
type = NASA9;
tlow = m_lowT;
thigh = m_highT;
pref = m_Pref;
for (int i = 0; i < 9; i++) {
coeffs[i] = m_coeff[i];
}
}
// Modify parameters for the standard state
/*
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state.
*/
void Nasa9Poly1::modifyParameters(doublereal* coeffs) {
for (int i = 0; i < 9; i++) {
m_coeff[i] = coeffs[i];
}
}
}

228
Cantera/src/thermo/Nasa9Poly1.h Normal file
View file

@ -0,0 +1,228 @@
/**
* @file Nasa9Poly1.h
* Header for a single-species standard state object derived
* from
* \link Cantera::SpeciesThermoInterpType SpeciesThermoInterpType\endlink based
* on the NASA 9 coefficient temperature polynomial form applied to
* one temperature region
* (see \ref spthermo and class \link Cantera::Nasa9Poly1 Nasa9Poly1\endlink).
*
* This parameterization has one NASA temperature region.
*/
#ifndef CT_NASA9POLY1_H
#define CT_NASA9POLY1_H
/* $Author$
* $Revision$
* $Date$
*/
// Copyright 2007 Sandia National Laboratories
#include "global.h"
#include "SpeciesThermoInterpType.h"
namespace Cantera {
//! The NASA 9 polynomial parameterization for one temperature range.
/*!
* This parameterization expresses the heat capacity via a
* 7 coefficient polynomial.
*
* Note that this is the form used in the
* 2002 NASA equilibrium program
*
* NASA Glenn Coefficients for Calculating Thermodynamic
* Properties of Individual Species,
* B. J. McBride, M. J. Zehe, S. Gordon
* NASA/TP-2002-211556, Sept. 2002
*
*
* Nine coefficients \f$(a_0,\dots,a_6)\f$ are used to represent
* \f$ C_p^0(T)\f$, \f$ H^0(T)\f$, and \f$ S^0(T) \f$ as
* polynomials in \f$ T \f$ :
* \f[
* \frac{c_p(T)}{R} = a_0 T^{-2} + a_1 T^{-1} + a_2 + a_3 T
* + a_4 T^2 + a_5 T^3 + a_6 T^4
* \f]
*
* \f[
* \frac{H^0(T)}{RT} = - a_0 T^{-2} + a_1 \frac{\ln(T)}{T} + a_2
* + a_3 T + a_4 T^2 + a_5 T^3 + a_6 T^4 + \frac{a_7}{T}
* \f]
*
* \f[
* \frac{s^0(T)}{R} = - \frac{a_0}{2} T^{-2} - a_1 T^{-1} + a_2 \ln(T)
+ + a_3 T \frac{a_4}{2} T^2 + \frac{a_5}{3} T^3 + \frac{a_6}{4} T^4 + a_8
* \f]
*
* The standard state is assumed to be the ideal gas at the
* standard pressure of 1 bar, for gases.
* For condensed species, the standard state is the
* pure cyrstalline or liquid substance at the standard
* pressure of 1 atm.
*
* These NASA representations may have more than 2 temperature regions.
*
* @ingroup spthermo
*/
class Nasa9Poly1 : public SpeciesThermoInterpType {
public:
//! Empty constructor
Nasa9Poly1();
//! constructor used in templated instantiations
/*!
* @param n Species index
* @param tlow Minimum temperature
* @param thigh Maximum temperature
* @param pref reference pressure (Pa).
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state.
*/
Nasa9Poly1(int n, doublereal tlow, doublereal thigh, doublereal pref,
const doublereal* coeffs);
//! copy constructor
/*!
* @param b object to be copied
*/
Nasa9Poly1(const Nasa9Poly1& b);
//! assignment operator
/*!
* @param b object to be copied
*/
Nasa9Poly1& operator=(const Nasa9Poly1& b);
//! Destructor
virtual ~Nasa9Poly1();
//! duplicator
virtual SpeciesThermoInterpType *
duplMyselfAsSpeciesThermoInterpType() const;
//! Returns the minimum temperature that the thermo
//! parameterization is valid
virtual doublereal minTemp() const;
//! Returns the maximum temperature that the thermo
//! parameterization is valid
virtual doublereal maxTemp() const;
//! Returns the reference pressure (Pa)
virtual doublereal refPressure() const;
//! Returns an integer representing the type of parameterization
virtual int reportType() const;
//! Returns an integer representing the species index
virtual int speciesIndex() const;
//! Update the properties for this species, given a temperature polynomial
/*!
* This method is called with a pointer to an array containing the functions of
* temperature needed by this parameterization, and three pointers to arrays where the
* computed property values should be written. This method updates only one value in
* each array.
*
* Temperature Polynomial:
* tt[0] = t;
* tt[1] = t*t;
* tt[2] = t*t*t;
* tt[3] = t*t*t*t;
* tt[4] = 1.0/t;
* tt[5] = 1.0/(t*t);
* tt[6] = std::log(t);
*
* @param tt vector of temperature polynomials
* @param cp_R Vector of Dimensionless heat capacities.
* (length m_kk).
* @param h_RT Vector of Dimensionless enthalpies.
* (length m_kk).
* @param s_R Vector of Dimensionless entropies.
* (length m_kk).
*/
virtual void updateProperties(const doublereal* tt,
doublereal* cp_R, doublereal* h_RT, doublereal* s_R) const;
//! Compute the reference-state property of one species
/*!
* Given temperature T in K, this method updates the values of
* the non-dimensional heat capacity at constant pressure,
* enthalpy, and entropy, at the reference pressure, Pref
* of one of the species. The species index is used
* to reference into the cp_R, h_RT, and s_R arrays.
*
* Temperature Polynomial:
* tt[0] = t;
* tt[1] = t*t;
* tt[2] = t*t*t;
* tt[3] = t*t*t*t;
* tt[4] = 1.0/t;
* tt[5] = 1.0/(t*t);
* tt[6] = std::log(t);
*
* @param temp Temperature (Kelvin)
* @param cp_R Vector of Dimensionless heat capacities.
* (length m_kk).
* @param h_RT Vector of Dimensionless enthalpies.
* (length m_kk).
* @param s_R Vector of Dimensionless entropies.
* (length m_kk).
*/
virtual void updatePropertiesTemp(const doublereal temp,
doublereal* cp_R, doublereal* h_RT,
doublereal* s_R) const;
//!This utility function reports back the type of
//! parameterization and all of the parameters for the
//! species, index.
/*!
* All parameters are output variables
*
* @param n Species index
* @param type Integer type of the standard type
* @param tlow output - Minimum temperature
* @param thigh output - Maximum temperature
* @param pref output - reference pressure (Pa).
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state.
*/
virtual void reportParameters(int &n, int &type,
doublereal &tlow, doublereal &thigh,
doublereal &pref,
doublereal* const coeffs) const;
//! Modify parameters for the standard state
/*!
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state.
*/
virtual void modifyParameters(doublereal* coeffs);
protected:
//! lowest valid temperature
doublereal m_lowT;
//! highest valid temperature
doublereal m_highT;
//! standard-state pressure
doublereal m_Pref;
//! species index
int m_index;
//! array of polynomial coefficients
array_fp m_coeff;
};
}
#endif

379
Cantera/src/thermo/Nasa9PolyMultiTempRegion.cpp Normal file
View file

@ -0,0 +1,379 @@
/**
* @file Nasa9Poly1.h
* Header for a single-species standard state object derived
* from \link Cantera::SpeciesThermoInterpType
* SpeciesThermoInterpType\endlink based
* on the NASA 9 coefficient temperature polynomial form
* applied to one temperature region
* (see \ref spthermo and class
* \link Cantera::Nasa9Poly1 Nasa9Poly1\endlink).
*
* This parameterization has one NASA temperature region.
*/
/* $Author$
* $Revision$
* $Date$
*/
// Copyright 2007 Sandia National Laboratories
#include "global.h"
#include "ctexceptions.h"
#include "Nasa9PolyMultiTempRegion.h"
namespace Cantera {
// The NASA 9 polynomial parameterization for one temperature range.
/*
* This parameterization expresses the heat capacity via a
* 7 coefficient polynomial.
*
* Note that this is the form used in the
* 2002 NASA equilibrium program
*
* NASA Glenn Coefficients for Calculating Thermodynamic
* Properties of Individual Species,
* B. J. McBride, M. J. Zehe, S. Gordon
* NASA/TP-2002-211556, Sept. 2002
*
*
* Nine coefficients \f$(a_0,\dots,a_6)\f$ are used to represent
* \f$ C_p^0(T)\f$, \f$ H^0(T)\f$, and \f$ S^0(T) \f$ as
* polynomials in \f$ T \f$ :
* \f[
* \frac{c_p(T)}{R} = a_0 T^{-2} + a_1 T^{-1} + a_2 + a_3 T
* + a_4 T^2 + a_5 T^3 + a_6 T^4
* \f]
*
* \f[
* \frac{H^0(T)}{RT} = - a_0 T^{-2} + a_1 \frac{\ln(T)}{T} + a_2
* + a_3 T + a_4 T^2 + a_5 T^3 + a_6 T^4 + \frac{a_7}{T}
* \f]
*
* \f[
* \frac{s^0(T)}{R} = - \frac{a_0}{2} T^{-2} - a_1 T^{-1} + a_2 \ln(T)
+ + a_3 T \frac{a_4}{2} T^2 + \frac{a_5}{3} T^3
* + \frac{a_6}{4} T^4 + a_8
* \f]
*
* The standard state is assumed to be the ideal gas at the
* standard pressure of 1 bar, for gases.
* For condensed species, the standard state is the
* pure cyrstalline or liquid substance at the standard
* pressure of 1 atm.
*
* These NASA representations may have more than 2 temperature regions.
*
* @ingroup spthermo
*/
//! Empty constructor
Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion() :
m_lowT(0.0),
m_highT (0.0),
m_Pref(0.0),
m_index (0),
m_numTempRegions(0),
m_currRegion(0)
{
}
// constructor used in templated instantiations
/*
* @param n Species index
* @param tlow Minimum temperature
* @param thigh Maximum temperature
* @param pref reference pressure (Pa).
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state.
*/
Nasa9PolyMultiTempRegion::
Nasa9PolyMultiTempRegion(std::vector<Cantera::Nasa9Poly1 *> &regionPts) :
m_lowT(0.0),
m_highT (0.0),
m_Pref(0.0),
m_index(0),
m_numTempRegions(0),
m_currRegion(0)
{
m_numTempRegions = regionPts.size();
// Do a shallow copy of the pointers. From now on, we will
// own these pointers and be responsible for deleting them.
m_regionPts = regionPts;
m_lowerTempBounds.resize(m_numTempRegions);
m_lowT = m_regionPts[0]->minTemp();
m_highT = m_regionPts[m_numTempRegions-1]->maxTemp();
m_Pref = m_regionPts[0]->refPressure();
m_index = m_regionPts[0]->speciesIndex();
for (int i = 0; i < m_numTempRegions; i++) {
m_lowerTempBounds[i] = m_regionPts[i]->minTemp();
if (m_regionPts[i]->speciesIndex() != m_index) {
throw CanteraError("Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion",
"m_index inconsistency");
}
if (fabs(m_regionPts[i]->refPressure() - m_Pref) > 0.0001) {
throw CanteraError("Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion",
"refPressure inconsistency");
}
if (i > 0) {
if (m_lowerTempBounds[i-1] >= m_lowerTempBounds[i]) {
throw CanteraError("Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion",
"minTemp bounds inconsistency");
}
if (fabs(m_regionPts[i-1]->maxTemp() - m_lowerTempBounds[i]) > 0.0001) {
throw CanteraError("Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion",
"Temp bounds inconsistency");
}
}
}
}
// copy constructor
/*
* @param b object to be copied
*/
Nasa9PolyMultiTempRegion::
Nasa9PolyMultiTempRegion(const Nasa9PolyMultiTempRegion& b) :
m_lowT (b.m_lowT),
m_highT (b.m_highT),
m_Pref (b.m_Pref),
m_index (b.m_index),
m_numTempRegions(b.m_numTempRegions),
m_lowerTempBounds (b.m_lowerTempBounds),
m_currRegion(b.m_currRegion)
{
m_regionPts.resize(m_numTempRegions);
for (int i = 0; i < m_numTempRegions; i++) {
Nasa9Poly1 * dptr = b.m_regionPts[i];
m_regionPts[i] = new Nasa9Poly1(*dptr);
}
}
// assignment operator
/*
* @param b object to be copied
*/
Nasa9PolyMultiTempRegion&
Nasa9PolyMultiTempRegion::operator=(const Nasa9PolyMultiTempRegion& b) {
if (&b != this) {
for (int i = 0; i < m_numTempRegions; i++) {
delete m_regionPts[i];
m_regionPts[i] = 0;
}
m_lowT = b.m_lowT;
m_highT = b.m_highT;
m_Pref = b.m_Pref;
m_index = b.m_index;
m_numTempRegions = b.m_numTempRegions;
m_lowerTempBounds = b.m_lowerTempBounds;
m_currRegion = b.m_currRegion;
m_regionPts.resize(m_numTempRegions);
for (int i = 0; i < m_numTempRegions; i++) {
m_regionPts[i] = new Nasa9Poly1(*(b.m_regionPts[i]));
}
}
return *this;
}
// Destructor
Nasa9PolyMultiTempRegion::~Nasa9PolyMultiTempRegion() {
for (int i = 0; i < m_numTempRegions; i++) {
delete m_regionPts[i];
m_regionPts[i] = 0;
}
}
// duplicator
SpeciesThermoInterpType *
Nasa9PolyMultiTempRegion::duplMyselfAsSpeciesThermoInterpType() const {
Nasa9PolyMultiTempRegion* np = new Nasa9PolyMultiTempRegion(*this);
return (SpeciesThermoInterpType *) np;
}
// Returns the minimum temperature that the thermo
// parameterization is valid
doublereal Nasa9PolyMultiTempRegion::minTemp() const {
return m_lowT;
}
// Returns the maximum temperature that the thermo
// parameterization is valid
doublereal Nasa9PolyMultiTempRegion::maxTemp() const {
return m_highT;
}
// Returns the reference pressure (Pa)
doublereal Nasa9PolyMultiTempRegion::refPressure() const {
return m_Pref;
}
// Returns an integer representing the type of parameterization
int Nasa9PolyMultiTempRegion::reportType() const {
return NASA9MULTITEMP;
}
// Returns an integer representing the species index
int Nasa9PolyMultiTempRegion::speciesIndex() const {
return m_index;
}
// Update the properties for this species, given a temperature polynomial
/*
* This method is called with a pointer to an array containing the functions of
* temperature needed by this parameterization, and three pointers to arrays where the
* computed property values should be written. This method updates only one value in
* each array.
*
* Temperature Polynomial:
* tt[0] = t;
* tt[1] = t*t;
* tt[2] = t*t*t;
* tt[3] = t*t*t*t;
* tt[4] = 1.0/t;
* tt[5] = 1.0/(t*t);
* tt[6] = std::log(t);
*
* @param tt vector of temperature polynomials
* @param cp_R Vector of Dimensionless heat capacities.
* (length m_kk).
* @param h_RT Vector of Dimensionless enthalpies.
* (length m_kk).
* @param s_R Vector of Dimensionless entropies.
* (length m_kk).
*/
void Nasa9PolyMultiTempRegion::updateProperties(const doublereal* tt,
doublereal* cp_R,
doublereal* h_RT,
doublereal* s_R) const {
// Let's put some additional debugging here.
// This is an external routine
#ifdef DEBUG_HKM
double temp = tt[0];
if (temp < m_regionPts[m_currRegion]->minTemp() ) {
if (m_currRegion != 0) {
throw CanteraError("Nasa9PolyMultiTempRegion::updateProperties",
"region problem");
}
}
if (temp > m_regionPts[m_currRegion]->maxTemp() ) {
if (m_currRegion != m_numTempRegions - 1) {
throw CanteraError("Nasa9PolyMultiTempRegion::updateProperties",
"region problem");
}
}
#endif
m_regionPts[m_currRegion]->updateProperties(tt, cp_R, h_RT, s_R);
}
// Compute the reference-state property of one species
/*
* Given temperature T in K, this method updates the values of
* the non-dimensional heat capacity at constant pressure,
* enthalpy, and entropy, at the reference pressure, Pref
* of one of the species. The species index is used
* to reference into the cp_R, h_RT, and s_R arrays.
*
* Temperature Polynomial:
* tt[0] = t;
* tt[1] = t*t;
* tt[2] = t*t*t;
* tt[3] = t*t*t*t;
* tt[4] = 1.0/t;
* tt[5] = 1.0/(t*t);
* tt[6] = std::log(t);
*
* @param temp Temperature (Kelvin)
* @param cp_R Vector of Dimensionless heat capacities.
* (length m_kk).
* @param h_RT Vector of Dimensionless enthalpies.
* (length m_kk).
* @param s_R Vector of Dimensionless entropies.
* (length m_kk).
*/
void Nasa9PolyMultiTempRegion::updatePropertiesTemp(const doublereal temp,
doublereal* cp_R, doublereal* h_RT,
doublereal* s_R) const {
double tPoly[7];
tPoly[0] = temp;
tPoly[1] = temp * temp;
tPoly[2] = tPoly[1] * temp;
tPoly[3] = tPoly[2] * temp;
tPoly[4] = 1.0 / temp;
tPoly[5] = tPoly[4] / temp;
tPoly[6] = std::log(temp);
// Now find the region
m_currRegion = 0;
for (int i = 1; i < m_numTempRegions; i++) {
if (temp < m_lowerTempBounds[i]) {
break;
}
m_currRegion++;
}
updateProperties(tPoly, cp_R, h_RT, s_R);
}
//This utility function reports back the type of
// parameterization and all of the parameters for the
// species, index.
/*
* All parameters are output variables
*
* @param n Species index
* @param type Integer type of the standard type
* @param tlow output - Minimum temperature
* @param thigh output - Maximum temperature
* @param pref output - reference pressure (Pa).
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state.
*/
void Nasa9PolyMultiTempRegion::reportParameters(int &n, int &type,
doublereal &tlow, doublereal &thigh,
doublereal &pref,
doublereal* const coeffs) const {
n = m_index;
type = NASA9;
tlow = m_lowT;
thigh = m_highT;
pref = m_Pref;
coeffs[0] = m_numTempRegions;
int index = 1;
int n_tmp = 0;;
int type_tmp = 0;
double pref_tmp = 0.0;
for (int iReg = 0; iReg < m_numTempRegions; iReg++) {
m_regionPts[iReg]->reportParameters(n_tmp, type_tmp,
coeffs[index], coeffs[index+1],
pref_tmp, coeffs + index + 2);
index += 11;
}
}
// Modify parameters for the standard state
/*
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state.
*/
void Nasa9PolyMultiTempRegion::modifyParameters(doublereal* coeffs) {
int index = 3;
for (int iReg = 0; iReg < m_numTempRegions; iReg++) {
m_regionPts[iReg]->modifyParameters(coeffs + index);
index += 11;
}
}
}

204
Cantera/src/thermo/Nasa9PolyMultiTempRegion.h Normal file
View file

@ -0,0 +1,204 @@
/**
* @file Nasa9Poly1.h
* Header for a single-species standard state object derived
* from \link Cantera::SpeciesThermoInterpType
* SpeciesThermoInterpType\endlink based
* on the NASA 9 coefficient temperature polynomial form
* applied to one temperature region
* (see \ref spthermo and class \link Cantera::Nasa9Poly1 Nasa9Poly1\endlink).
*
* This parameterization has one NASA temperature region.
*/
#ifndef CT_NASA9POLYMULTITEMPREGION_H
#define CT_NASA9POLYMULTITEMPREGION_H
/* $Author$
* $Revision$
* $Date$
*/
// Copyright 2007 Sandia National Laboratories
#include "global.h"
#include "Nasa9Poly1.h"
namespace Cantera {
//! The NASA 9 polynomial parameterization for multiple temperature ranges.
/*!
*
* @ingroup spthermo
*/
class Nasa9PolyMultiTempRegion : public SpeciesThermoInterpType {
public:
//! Empty constructor
Nasa9PolyMultiTempRegion();
//! constructor used in templated instantiations
/*!
* @param n Species index
* @param tlow Minimum temperature
* @param thigh Maximum temperature
* @param pref reference pressure (Pa).
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state.
*/
Nasa9PolyMultiTempRegion(std::vector<Cantera::Nasa9Poly1 *> &regionPts);
//! copy constructor
/*!
* @param b object to be copied
*/
Nasa9PolyMultiTempRegion(const Nasa9PolyMultiTempRegion& b);
//! assignment operator
/*!
* @param b object to be copied
*/
Nasa9PolyMultiTempRegion& operator=(const Nasa9PolyMultiTempRegion& b);
//! Destructor
virtual ~Nasa9PolyMultiTempRegion();
//! duplicator
virtual SpeciesThermoInterpType *
duplMyselfAsSpeciesThermoInterpType() const;
//! Returns the minimum temperature that the thermo
//! parameterization is valid
virtual doublereal minTemp() const;
//! Returns the maximum temperature that the thermo
//! parameterization is valid
virtual doublereal maxTemp() const;
//! Returns the reference pressure (Pa)
virtual doublereal refPressure() const;
//! Returns an integer representing the type of parameterization
virtual int reportType() const;
//! Returns an integer representing the species index
virtual int speciesIndex() const;
//! Update the properties for this species, given a temperature polynomial
/*!
* This method is called with a pointer to an array containing the functions of
* temperature needed by this parameterization, and three pointers to arrays where the
* computed property values should be written. This method updates only one value in
* each array.
*
* Temperature Polynomial:
* tt[0] = t;
* tt[1] = t*t;
* tt[2] = t*t*t;
* tt[3] = t*t*t*t;
* tt[4] = 1.0/t;
* tt[5] = 1.0/(t*t);
* tt[6] = std::log(t);
*
* @param tt vector of temperature polynomials
* @param cp_R Vector of Dimensionless heat capacities.
* (length m_kk).
* @param h_RT Vector of Dimensionless enthalpies.
* (length m_kk).
* @param s_R Vector of Dimensionless entropies.
* (length m_kk).
*/
virtual void updateProperties(const doublereal* tt,
doublereal* cp_R, doublereal* h_RT,
doublereal* s_R) const;
//! Compute the reference-state property of one species
/*!
* Given temperature T in K, this method updates the values of
* the non-dimensional heat capacity at constant pressure,
* enthalpy, and entropy, at the reference pressure, Pref
* of one of the species. The species index is used
* to reference into the cp_R, h_RT, and s_R arrays.
*
* Temperature Polynomial:
* tt[0] = t;
* tt[1] = t*t;
* tt[2] = t*t*t;
* tt[3] = t*t*t*t;
* tt[4] = 1.0/t;
* tt[5] = 1.0/(t*t);
* tt[6] = std::log(t);
*
* @param temp Temperature (Kelvin)
* @param cp_R Vector of Dimensionless heat capacities.
* (length m_kk).
* @param h_RT Vector of Dimensionless enthalpies.
* (length m_kk).
* @param s_R Vector of Dimensionless entropies.
* (length m_kk).
*/
virtual void updatePropertiesTemp(const doublereal temp,
doublereal* cp_R, doublereal* h_RT,
doublereal* s_R) const;
//!This utility function reports back the type of
//! parameterization and all of the parameters for the
//! species, index.
/*!
* All parameters are output variables
*
* @param n Species index
* @param type Integer type of the standard type
* @param tlow output - Minimum temperature
* @param thigh output - Maximum temperature
* @param pref output - reference pressure (Pa).
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state.
*/
virtual void reportParameters(int &n, int &type,
doublereal &tlow, doublereal &thigh,
doublereal &pref,
doublereal* const coeffs) const;
//! Modify parameters for the standard state
/*!
* @param coeffs Vector of coefficients used to set the
* parameters for the standard state.
*/
virtual void modifyParameters(doublereal* coeffs);
protected:
//! lowest valid temperature
doublereal m_lowT;
//! highest valid temperature
doublereal m_highT;
//! standard-state pressure
doublereal m_Pref;
//! species index
int m_index;
//! Number of temperature regions
int m_numTempRegions;
//! Lower boundaries of each temperature regions
vector_fp m_lowerTempBounds;
//! pointers to the objects
/*!
* This object will now own these pointers and delete
* them when the current object is deleted.
*/
std::vector<Nasa9Poly1 *>m_regionPts;
//! current region
mutable int m_currRegion;
};
}
#endif

4
Cantera/src/thermo/NasaPoly1.h
View file

@ -138,7 +138,9 @@ namespace Cantera {
//! Returns an integer representing the type of parameterization
virtual int reportType() const { return NASA1; }
//! Returns an integer representing the species index
virtual int speciesIndex() const { return m_index; }
//! Update the properties for this species, given a temperature polynomial
/*!
* This method is called with a pointer to an array containing the functions of

3
Cantera/src/thermo/NasaPoly2.h
View file

@ -170,6 +170,9 @@ namespace Cantera {
//! Returns an integer representing the type of parameterization
virtual int reportType() const { return NASA2; }
//! Returns an integer representing the species index
virtual int speciesIndex() const { return m_index; }
//! Update the properties for this species, given a temperature polynomial
/*!

11
Cantera/src/thermo/NasaThermo.h
View file

@ -218,6 +218,17 @@ namespace Cantera {
m_p0 = refPressure;
}
//! Install a new species thermodynamic property
//! parameterization for one species.
/*!
* @param stit_ptr Pointer to the SpeciesThermoInterpType object
* This will set up the thermo for one species
*/
virtual void install_STIT(SpeciesThermoInterpType *stit_ptr) {
throw CanteraError("install_STIT", "not implemented");
}
//! Like update(), but only updates the single species k.
/*!
* @param k species index

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

@ -157,7 +157,10 @@ namespace Cantera {
//! Returns an integer representing the type of parameterization
virtual int reportType() const { return SHOMATE; }
//! Returns an integer representing the species index
virtual int speciesIndex() const { return m_index; }
//! Update the properties for this species, given a temperature polynomial
/*!
@ -463,7 +466,9 @@ namespace Cantera {
//! Returns an integer representing the type of parameterization
virtual int reportType() const { return SHOMATE2; }
//! Returns an integer representing the species index
virtual int speciesIndex() const { return m_index; }
//! Update the properties for this species, given a temperature polynomial
/*!
* This method is called with a pointer to an array containing the functions of

12
Cantera/src/thermo/ShomateThermo.h
View file

@ -212,6 +212,16 @@ namespace Cantera {
}
//! Install a new species thermodynamic property
//! parameterization for one species.
/*!
* @param stit_ptr Pointer to the SpeciesThermoInterpType object
* This will set up the thermo for one species
*/
virtual void install_STIT(SpeciesThermoInterpType *stit_ptr) {
throw CanteraError("install_STIT", "not implemented");
}
//! Like update(), but only updates the single species k.
/*!
* @param k species index
@ -352,7 +362,7 @@ namespace Cantera {
* @param index Species index
*/
virtual int reportType(int index) const { return SHOMATE; }
/*!
* This utility function reports back the type of
* parameterization and all of the parameters for the

10
Cantera/src/thermo/SimpleThermo.h
View file

@ -188,6 +188,16 @@ namespace Cantera {
m_p0 = refPressure;
}
//! Install a new species thermodynamic property
//! parameterization for one species.
/*!
* @param stit_ptr Pointer to the SpeciesThermoInterpType object
* This will set up the thermo for one species
*/
virtual void install_STIT(SpeciesThermoInterpType *stit_ptr) {
throw CanteraError("install_STIT", "not implemented");
}
//! Compute the reference-state properties for all species.
/*!
* Given temperature T in K, this method updates the values of

10
Cantera/src/thermo/SpeciesThermo.h
View file

@ -23,6 +23,8 @@
namespace Cantera {
class SpeciesThermoInterpType;
/**
* @defgroup spthermo Species Standard-State Thermodynamic Properties
*
@ -206,6 +208,14 @@ namespace Cantera {
doublereal maxTemp,
doublereal refPressure)=0;
//! Install a new species thermodynamic property
//! parameterization for one species.
/*!
* @param stit_ptr Pointer to the SpeciesThermoInterpType object
* This will set up the thermo for one species
*/
virtual void install_STIT(SpeciesThermoInterpType *stit_ptr) = 0;
//! Compute the reference-state properties for all species.
/*!

173
Cantera/src/thermo/SpeciesThermoFactory.cpp
View file

@ -23,6 +23,8 @@ using namespace std;
#include "SimpleThermo.h"
#include "GeneralSpeciesThermo.h"
#include "Mu0Poly.h"
#include "Nasa9PolyMultiTempRegion.h"
#include "Nasa9Poly1.h"
#include "SpeciesThermoMgr.h"
#include "speciesThermoTypes.h"
@ -64,19 +66,21 @@ namespace Cantera {
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;
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;
if (th.hasChild("NASA9")) has_other = 1;
if (th.hasChild("NASA9MULTITEMP")) has_other = 1;
} else {
throw UnknownSpeciesThermoModel("getSpeciesThermoTypes:",
spNode->attrib("name"), "missing");
throw UnknownSpeciesThermoModel("getSpeciesThermoTypes:",
spNode->attrib("name"), "missing");
}
}
}
@ -429,6 +433,58 @@ namespace Cantera {
sp.install(speciesName, k, SIMPLE, &c[0], tmin, tmax, p0);
}
/**
* Install a NASA9 polynomial thermodynamic property
* parameterization for species k into a SpeciesThermo instance.
* This is called by method installThermoForSpecies if a NASA9
* block is found in the XML input.
*/
static void installNasa9ThermoFromXML(std::string speciesName,
SpeciesThermo& sp, int k,
const std::vector<XML_Node*>& tp)
{
const XML_Node * fptr = tp[0];
int nRegTmp = tp.size();
int nRegions = 0;
vector_fp cPoly;
Nasa9Poly1 *np_ptr = 0;
std::vector<Nasa9Poly1 *> regionPtrs;
doublereal tmin, tmax, pref;
// Loop over all of the possible temperature regions
for (int i = 0; i < nRegTmp; i++) {
fptr = tp[i];
if (fptr) {
if (fptr->name() == "NASA9") {
if (fptr->hasChild("floatArray")) {
tmin = fpValue((*fptr)["Tmin"]);
tmax = fpValue((*fptr)["Tmax"]);
pref = fpValue((*fptr)["P0"]);
getFloatArray(fptr->child("floatArray"), cPoly, false);
if (cPoly.size() != 9) {
throw CanteraError("installNasa9ThermoFromXML",
"Expected 9 coeff polynomial");
}
np_ptr = new Nasa9Poly1(k, tmin, tmax, pref,
DATA_PTR(cPoly));
regionPtrs.push_back(np_ptr);
nRegions++;
}
}
}
}
if (nRegions == 0) {
throw UnknownSpeciesThermoModel("installThermoForSpecies",
speciesName, " " );
} else if (nRegions == 1) {
sp.install_STIT(np_ptr);
} else {
Nasa9PolyMultiTempRegion* npMulti_ptr = new Nasa9PolyMultiTempRegion(regionPtrs);
sp.install_STIT(npMulti_ptr);
}
}
/**
* Install a species thermodynamic property parameterization
* for one species into a species thermo manager.
@ -436,60 +492,71 @@ namespace Cantera {
* @param s XML node specifying species
* @param spthermo species thermo manager
*/
void SpeciesThermoFactory::
installThermoForSpecies(int k, const XML_Node& s,
SpeciesThermo& spthermo) {
void SpeciesThermoFactory::
installThermoForSpecies(int k, const XML_Node& s,
SpeciesThermo& spthermo) {
/*
* Check to see that the species block has a thermo block
* before processing. Throw an error if not there.
*/
if (!(s.hasChild("thermo"))) {
throw UnknownSpeciesThermoModel("installThermoForSpecies",
s["name"], "<nonexistent>");
s["name"], "<nonexistent>");
}
const XML_Node& thermo = s.child("thermo");
const std::vector<XML_Node*>& tp = thermo.children();
int nc = static_cast<int>(tp.size());
if (nc == 1) {
const XML_Node* f = tp[0];
if (f->name() == "Shomate") {
installShomateThermoFromXML(s["name"], spthermo, k, f, 0);
}
else if (f->name() == "const_cp") {
installSimpleThermoFromXML(s["name"], spthermo, k, *f);
}
else if (f->name() == "NASA") {
installNasaThermoFromXML(s["name"], spthermo, k, f, 0);
}
else if (f->name() == "Mu0") {
installMu0ThermoFromXML(s["name"], spthermo, k, f);
}
else {
throw UnknownSpeciesThermoModel("installThermoForSpecies",
s["name"], f->name());
}
const XML_Node* f = tp[0];
if (f->name() == "Shomate") {
installShomateThermoFromXML(s["name"], spthermo, k, f, 0);
}
else if (f->name() == "const_cp") {
installSimpleThermoFromXML(s["name"], spthermo, k, *f);
}
else if (f->name() == "NASA") {
installNasaThermoFromXML(s["name"], spthermo, k, f, 0);
}
else if (f->name() == "Mu0") {
installMu0ThermoFromXML(s["name"], spthermo, k, f);
}
else if (f->name() == "NASA9") {
installNasa9ThermoFromXML(s["name"], spthermo, k, tp);
}
else {
throw UnknownSpeciesThermoModel("installThermoForSpecies",
s["name"], f->name());
}
}
else if (nc == 2) {
const XML_Node* f0 = tp[0];
const XML_Node* f1 = tp[1];
if (f0->name() == "NASA" && f1->name() == "NASA") {
installNasaThermoFromXML(s["name"], spthermo, k, f0, f1);
}
else if (f0->name() == "Shomate" && f1->name() == "Shomate") {
installShomateThermoFromXML(s["name"], spthermo, k, f0, f1);
}
else {
throw UnknownSpeciesThermoModel("installThermoForSpecies", s["name"],
f0->name() + " and "
+ f1->name());
}
}
else {
throw UnknownSpeciesThermoModel("installThermoForSpecies", s["name"],
"multiple");
}
const XML_Node* f0 = tp[0];
const XML_Node* f1 = tp[1];
if (f0->name() == "NASA" && f1->name() == "NASA") {
installNasaThermoFromXML(s["name"], spthermo, k, f0, f1);
}
else if (f0->name() == "Shomate" && f1->name() == "Shomate") {
installShomateThermoFromXML(s["name"], spthermo, k, f0, f1);
}
else if (f0->name() == "NASA9" && f1->name() == "NASA9") {
installNasa9ThermoFromXML(s["name"], spthermo, k, tp);
} else {
throw UnknownSpeciesThermoModel("installThermoForSpecies", s["name"],
f0->name() + " and "
+ f1->name());
}
}
else if (nc >= 2) {
const XML_Node* f0 = tp[0];
if (f0->name() == "NASA9") {
installNasa9ThermoFromXML(s["name"], spthermo, k, tp);
} else {
throw UnknownSpeciesThermoModel("installThermoForSpecies", s["name"],
"multiple");
}
} else {
throw UnknownSpeciesThermoModel("installThermoForSpecies", s["name"],
"multiple");
}
}
}

3
Cantera/src/thermo/SpeciesThermoInterpType.h
View file

@ -72,6 +72,9 @@ namespace Cantera {
//! Returns an integer representing the type of parameterization
virtual int reportType() const = 0;
//! Returns an integer representing the species index
virtual int speciesIndex() const = 0;
//! Update the properties for this species, given a temperature polynomial
/*!

10
Cantera/src/thermo/SpeciesThermoMgr.h
View file

@ -234,6 +234,16 @@ namespace Cantera {
}
}
//! Install a new species thermodynamic property
//! parameterization for one species.
/*!
* @param stit_ptr Pointer to the SpeciesThermoInterpType object
* This will set up the thermo for one species
*/
virtual void install_STIT(SpeciesThermoInterpType *stit_ptr) {
throw CanteraError("install_STIT", "not implemented");
}
//! Compute the reference-state properties for all species.
/*!
* Given temperature T in K, this method updates the values of

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

@ -57,6 +57,14 @@
//! This is implemented in the class NasaPoly1 in NasaPoly1.h
#define NASA1 256
//! 9 coefficient NASA Polynomials
//! This is implemented in the class Nasa9Poly1 in Nasa9Poly1.h
#define NASA9 512
//! 9 coefficient NASA Polynomials in multiple temperature regions
//! This is implemented in the class Nasa9PolyMultiTempRegion in Nasa9Poly1MultiTempRegion
#define NASA9MULTITEMP 513
#include "ct_defs.h"
#include "stringUtils.h"