689 lines
20 KiB
C++
Executable file
689 lines
20 KiB
C++
Executable file
/**
|
|
* @file CKReader.cpp
|
|
*
|
|
*/
|
|
|
|
// Copyright 2001 California Institute of Technology
|
|
|
|
#include <fstream>
|
|
#include <string>
|
|
using namespace std;
|
|
|
|
#include "CKParser.h"
|
|
#include "CKReader.h"
|
|
#include "thermoFunctions.h"
|
|
|
|
#include <cstring>
|
|
#include <cstdlib>
|
|
#include <ctime>
|
|
#include <iomanip>
|
|
|
|
#include "writelog.h"
|
|
#include <cstdio>
|
|
#include "ckr_defs.h"
|
|
|
|
//#include "global.h"
|
|
//#define APP Cantera::Application
|
|
|
|
namespace ckr {
|
|
|
|
|
|
/**
|
|
* read and optionally validate an input file in Chemkin format.
|
|
* @param inputFile path to the input file
|
|
* @param thermoDatabase path to the species database file
|
|
* @param log path to the file where log messages should be written
|
|
* @return true if no errors were encountered, false otherwise
|
|
*/
|
|
|
|
bool CKReader::read(const std::string& inputFile, const std::string& thermoDatabase,
|
|
const std::string& logfile) {
|
|
|
|
clock_t t0, t1;
|
|
|
|
t0 = clock();
|
|
|
|
ifstream ckinfile(inputFile.c_str());
|
|
ofstream log(logfile.c_str());
|
|
try {
|
|
// construct a parser for the input file
|
|
CKParser parser(&ckinfile, inputFile, &log);
|
|
parser.verbose = verbose;
|
|
parser.debug = debug;
|
|
|
|
// write header information to the log file
|
|
struct tm *newtime;
|
|
time_t aclock;
|
|
time( &aclock ); /* Get time in seconds */
|
|
newtime = localtime( &aclock ); /* Convert time to struct tm form */
|
|
|
|
|
|
log << "CKReader version 1.0" << endl
|
|
<< "http://www.cantera.org" << endl << endl
|
|
<< asctime(newtime) << endl
|
|
<< setw(20) << "input file: "
|
|
<< setw(30) << inputFile << endl;
|
|
|
|
if (thermoDatabase != "")
|
|
log << setw(20) << "species database: "
|
|
<< setw(30) << thermoDatabase << endl;
|
|
|
|
if (!validate)
|
|
log << endl << "*************** Warning ***************" << endl
|
|
<< " mechanism validation disabled" << endl
|
|
<< "*****************************************" << endl;
|
|
|
|
if (debug) {
|
|
log << "*** DEBUG MODE ***" << endl;
|
|
}
|
|
else {
|
|
log << "debugging disabled." << endl;
|
|
}
|
|
|
|
//----------- process ELEMENT section ----------------------
|
|
|
|
bool elok = parser.readElementSection(elements);
|
|
int nel = static_cast<int>(elements.size());
|
|
vector<string> elementSymbols;
|
|
for (int j = 0; j < nel; j++) elementSymbols.push_back(elements[j].name);
|
|
|
|
if (verbose) {
|
|
log.flags(ios::showpoint);
|
|
log.precision(6);
|
|
log.width(0);
|
|
|
|
log << endl << newTask("reading elements") << endl;
|
|
|
|
// write summary to log file
|
|
for (int i = 0; i < nel; i++) {
|
|
log << i+1 << ". " << pad(elements[i].name,2) << " ";
|
|
double wt = elements[i].atomicWeight;
|
|
if (wt == 0.0) log << "<error>";
|
|
else log << wt;
|
|
if (!elements[i].weightFromDB) log << " (specified)";
|
|
if (elements[i].comment != "")
|
|
log << " ! " << elements[i].comment;
|
|
log << endl;
|
|
}
|
|
}
|
|
log << "\nread " << nel << " elements." << endl;
|
|
|
|
if (!elok) {
|
|
log << "\nerrors were encountered." << endl;
|
|
return false;
|
|
}
|
|
if (nel == 0) {
|
|
return false;
|
|
}
|
|
|
|
|
|
//------------ process SPECIES section ------------------------
|
|
|
|
vector<string> speciesSymbols;
|
|
bool spok = parser.readSpeciesSection(species);
|
|
int nsp = static_cast<int>(species.size());
|
|
|
|
if (verbose) log << newTask("reading species") << endl;
|
|
|
|
for (int i = 0; i < nsp; i++) {
|
|
Species& s = species[i];
|
|
if (verbose) log << i+1 << ". " << s.name << endl;
|
|
speciesSymbols.push_back(s.name);
|
|
}
|
|
log << "\nread " << nsp << " species." << endl;
|
|
|
|
if (!spok) {
|
|
log << "\nerrors were encountered." << endl;
|
|
return false;
|
|
}
|
|
if (nsp == 0) return false;
|
|
|
|
|
|
//------------- process THERMO section -------------------------
|
|
|
|
if (verbose) log << newTask("looking up species definitions") << endl;
|
|
|
|
// if a thermo database is specified, get the default Tmin, Tmid, Tmax
|
|
vector_fp temp;
|
|
|
|
if (thermoDatabase != "") {
|
|
|
|
if (verbose) log << "reading default temperature ranges from "
|
|
<< thermoDatabase << endl;
|
|
|
|
ifstream thermofile(thermoDatabase.c_str());
|
|
CKParser thermReader(&thermofile, thermoDatabase, &log);
|
|
thermReader.verbose = verbose;
|
|
thermReader.debug = debug;
|
|
int dbflag = HasTempRange;
|
|
vector<string> dummy;
|
|
thermReader.readThermoSection(dummy, speciesData, temp, dbflag, log);
|
|
}
|
|
|
|
|
|
bool hasthermo = parser.advanceToKeyword("THERM","REAC");
|
|
|
|
int k, optionFlag = 0;
|
|
int undefined = static_cast<int>(species.size());
|
|
string nm;
|
|
vector<string> undef;
|
|
bool allsp = (speciesSymbols[0] == "<ALL>");
|
|
if (hasthermo &&
|
|
parser.readThermoSection(speciesSymbols,
|
|
speciesData, temp, optionFlag, log)) {
|
|
if (allsp) {
|
|
nsp = static_cast<int>(speciesData.size());
|
|
for (k = 0; k < nsp; k++) {
|
|
Species s;
|
|
s.name = speciesSymbols[k];
|
|
species.push_back(s);
|
|
}
|
|
}
|
|
undefined = 0;
|
|
for (k = 0; k < nsp; k++) {
|
|
nm = species[k].name;
|
|
species[k] = speciesData[species[k].name];
|
|
if (species[k].name == "<empty>") {
|
|
undefined++;
|
|
undef.push_back(nm);
|
|
species[k].name = nm;
|
|
}
|
|
}
|
|
int localdefs = nsp - undefined;
|
|
if (localdefs > 0 && verbose) log << "found definitions for "
|
|
<< localdefs
|
|
<< " of "
|
|
<< nsp
|
|
<< " species in the input file. "
|
|
<< endl;
|
|
}
|
|
else {
|
|
undef = speciesSymbols;
|
|
if (verbose) log << "no THERMO section in input file." << endl;
|
|
}
|
|
|
|
if (undefined > 0 && thermoDatabase != ""
|
|
&& optionFlag != NoThermoDatabase) {
|
|
|
|
if (verbose) log << "searching external database "
|
|
<< thermoDatabase << " for species definitions..."
|
|
<< endl;
|
|
|
|
ifstream thermofile(thermoDatabase.c_str());
|
|
CKParser thermoReader(&thermofile, thermoDatabase, &log);
|
|
thermoReader.verbose = verbose;
|
|
thermoReader.debug = debug;
|
|
int dbflag = HasTempRange;
|
|
thermoReader.readThermoSection(undef, speciesData, temp, dbflag, log);
|
|
undefined = 0;
|
|
if (allsp) {
|
|
species.clear();
|
|
nsp = static_cast<int>(speciesData.size());
|
|
for (k = 0; k < nsp; k++) {
|
|
Species s;
|
|
s.name = undef[k];
|
|
species.push_back(s);
|
|
}
|
|
}
|
|
for (int k = 0; k < nsp; k++) {
|
|
if (species[k].valid == 0) {
|
|
nm = species[k].name;
|
|
species[k] = speciesData[species[k].name];
|
|
if (species[k].name == "<empty>") {
|
|
undefined++;
|
|
species[k].name = nm;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if (validate && !validateSpecies(log)) {
|
|
//Cantera::setError("read","error in species");
|
|
return false;
|
|
}
|
|
|
|
|
|
//------------- process REACTIONS section -------------------------
|
|
|
|
if (verbose) log << newTask("reading reactions") << endl;
|
|
|
|
ckinfile.close();
|
|
ifstream ckinfile2(inputFile.c_str());
|
|
|
|
// construct a new parser for the input file
|
|
CKParser parser2(&ckinfile2, inputFile, &log);
|
|
parser2.verbose = verbose;
|
|
parser2.debug = debug;
|
|
|
|
parser2.readReactionSection(speciesSymbols, elementSymbols, reactions, units);
|
|
log << "\nread " << static_cast<int>(reactions.size())
|
|
<< " reactions." << endl;
|
|
|
|
bool rxnok = true;
|
|
if (validate) rxnok = rxnok && validateReactions(log);
|
|
bool writeok = true;
|
|
if (verbose || validate) writeok = writeReactions(log);
|
|
rxnok = rxnok && writeok;
|
|
if (!rxnok) return false;
|
|
|
|
log << "\nSuccess... ";
|
|
t1 = clock();
|
|
log << "elapsed CPU time = "
|
|
<< double(t1 - t0)/CLOCKS_PER_SEC
|
|
<< " s" << endl;
|
|
if (!validate) log << "*** no validation performed ***" << endl;
|
|
}
|
|
|
|
catch (CK_Exception e) {
|
|
log << e.errorMessage() << endl;
|
|
//Cantera::setError("CKReader::read",e.errorMessage());
|
|
return false;
|
|
}
|
|
catch (...) {
|
|
log << "an exception was raised in CKReader.";
|
|
return false;
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
|
|
/// print a summary of all reactions to the log file
|
|
bool CKReader::writeReactions(std::ostream& log) {
|
|
|
|
bool ok = true;
|
|
// int ns = species.size();
|
|
int nrxns = static_cast<int>(reactions.size());
|
|
log.flags(ios::unitbuf);
|
|
log.precision(6);
|
|
|
|
log << endl;
|
|
for (int n = 0; n < nrxns; n++) {
|
|
Reaction& r = reactions[n];
|
|
|
|
log << "reaction " << r.number << endl;
|
|
log << " ";
|
|
printReactionEquation(log, r);
|
|
log << endl;
|
|
|
|
// rate coefficient
|
|
if (r.isFalloffRxn) {
|
|
|
|
log << " high P rate coeff: ";
|
|
ok = ok && writeRateCoeff(r.kf, log) ;
|
|
|
|
log << " low P rate coeff: ";
|
|
ok = ok && writeRateCoeff(r.kf_aux, log);
|
|
|
|
ok = ok && writeFalloff(r.falloffType, r.falloffParameters, log);
|
|
|
|
}
|
|
else {
|
|
log << " rate coeff: ";
|
|
ok = ok && writeRateCoeff(r.kf, log);
|
|
}
|
|
if (r.isReversible && r.krev.A > 0) {
|
|
log << " reverse rate coeff: ";
|
|
ok = ok && writeRateCoeff(r.krev, log);
|
|
}
|
|
int ne = static_cast<int>(r.e3b.size());
|
|
|
|
if (ne > 0) {
|
|
vector<string> enhSpecies;
|
|
getMapKeys(r.e3b, enhSpecies);
|
|
log << " enhanced collision efficiencies:" << endl;
|
|
log << " ";
|
|
for (int nn = 0; nn < ne; nn++) {
|
|
log << enhSpecies[nn] << " "
|
|
<< r.e3b[enhSpecies[nn]];
|
|
if (nn < ne-1) log << ", ";
|
|
}
|
|
log << endl;
|
|
}
|
|
if (r.isDuplicate) log
|
|
<< " declared duplicate reaction. See reaction "
|
|
<< r.duplicate << "." << endl;
|
|
log << endl;
|
|
}
|
|
return ok;
|
|
}
|
|
|
|
|
|
|
|
/// validate the species
|
|
bool CKReader::validateSpecies(std::ostream& log) {
|
|
int nel = static_cast<int>(elements.size());
|
|
int nsp = static_cast<int>(species.size());
|
|
double nm, tol;
|
|
int j, k, m;
|
|
|
|
log << newTask("validating species");
|
|
|
|
// check for undeclared elements
|
|
vector<string> esyms;
|
|
|
|
log << " checking that all species have been defined... ";
|
|
for (k = 0; k < nsp; k++) {
|
|
Species& s = species[k];
|
|
if (s.valid == 0) {
|
|
log << endl << " species " << s.name << " undefined ";
|
|
s.valid = -1;
|
|
}
|
|
}
|
|
if (valid(species)) log << "OK" << endl;
|
|
else {
|
|
log << endl;
|
|
return false;
|
|
}
|
|
|
|
log << " checking that all species elements have been declared... ";
|
|
for (k = 0; k < nsp; k++) {
|
|
|
|
Species& s = species[k];
|
|
|
|
getMapKeys(s.comp, esyms);
|
|
nm = esyms.size();
|
|
|
|
for (m = 0; m < nm; m++) {
|
|
for (j = 0; j < nel; j++) {
|
|
if (esyms[m] == elements[j].name) break;
|
|
}
|
|
if (j == nel) {
|
|
log << endl << " species "
|
|
<< s.name << ": undeclared element "
|
|
<< esyms[m];
|
|
s.valid = -1;
|
|
}
|
|
}
|
|
}
|
|
if (valid(species)) log << "OK" << endl;
|
|
else {
|
|
log << endl;
|
|
return false;
|
|
}
|
|
|
|
log << " checking consistency of species thermo data... ";
|
|
tol = 0.01;
|
|
if (checkThermo(log, species, tol)) log << "OK" << endl;
|
|
else {
|
|
log << endl;
|
|
return false;
|
|
}
|
|
return true;
|
|
|
|
}
|
|
|
|
|
|
/// validate the reactions
|
|
bool CKReader::validateReactions(std::ostream& log) {
|
|
|
|
bool ok = true;
|
|
// int ns = species.size();
|
|
int nrxns = static_cast<int>(reactions.size());
|
|
|
|
vector<int> unbal;
|
|
log << "checking that all reactions balance...";
|
|
if (checkBalance(log, speciesData, reactions, unbal)) {
|
|
log << " OK" << endl;
|
|
}
|
|
else {
|
|
int nu = static_cast<int>(unbal.size());
|
|
for (int iu = 0; iu < nu; iu++) {
|
|
log << " error... reaction " << unbal[iu]
|
|
<< " does not balance" << endl;
|
|
}
|
|
ok = false;
|
|
}
|
|
|
|
log << "checking for duplicate reactions...";
|
|
|
|
for (int nn = 0; nn < nrxns; nn++) {
|
|
Reaction& r1 = reactions[nn];
|
|
for (int mm = nn + 1; mm < nrxns; mm++) {
|
|
Reaction& r2 = reactions[mm];
|
|
if (r1 == r2) {
|
|
r1.duplicate = mm + 1;
|
|
r2.duplicate = nn + 1;
|
|
if (!r1.isDuplicate || !r2.isDuplicate) {
|
|
log << endl << " error... undeclared duplicate reactions: "
|
|
<< nn + 1 << " and " << mm + 1;
|
|
ok = false;
|
|
}
|
|
else {
|
|
log << endl << " declared duplicate reactions: "
|
|
<< nn + 1
|
|
<< " and " << mm + 1;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
if (ok) log << "...OK" << endl;
|
|
|
|
return ok;
|
|
}
|
|
|
|
|
|
/**
|
|
* Check the thermodynamic property parameterizations for all species.
|
|
* The following are verified:
|
|
* - The heat capacity is positive throughout the full temperature range;
|
|
* - The entropy at Tmin is positive;
|
|
* - The heat capacity, entropy, and enthalpy evaluated at Tmid using
|
|
* both the high and low polynomial coefficients are the same to within
|
|
* relative error tol
|
|
* - The heat capacity at Tmax is not greater than the equipartition limit
|
|
* for the number of atoms in the molecule
|
|
*/
|
|
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());
|
|
const int n_points = 20;
|
|
|
|
int k;
|
|
bool ok = true;
|
|
for (k = 0; k < nsp; k++) {
|
|
Species& s = sp[k];
|
|
|
|
if (s.valid <= 0) {
|
|
ok = false;
|
|
log << endl << "species " << s.name
|
|
<< " contains an error." << endl;
|
|
}
|
|
if (!ok) return false;
|
|
}
|
|
|
|
log << endl << " Checking that cp/R is positive... ";
|
|
|
|
for (k = 0; k < nsp; k++) {
|
|
Species& s = sp[k];
|
|
|
|
// check that cp is positive
|
|
|
|
cp0 = 0.0;
|
|
for (int j = 0; j < n_points; j++) {
|
|
t = j*(s.thigh - s.tlow)/(n_points - 1) + s.tlow;
|
|
|
|
cp0 = cp(t, s);
|
|
if (cp0 < 0.0) {
|
|
log << endl << " error... Cp/R < 0 at T = " << t
|
|
<< " for species " << s.name << endl;
|
|
s.valid = -1;
|
|
ok = false;
|
|
}
|
|
}
|
|
}
|
|
if (ok) log << "ok" << endl;
|
|
else return ok;
|
|
|
|
|
|
// check that S(tlow) > 0
|
|
|
|
log << " Checking that the species entropies are positive... ";
|
|
for (k = 0; k < nsp; k++) {
|
|
Species& s = sp[k];
|
|
if (entropy(s.tlow, s) <= 0.0) {
|
|
log << endl << " error... negative entropy for species "
|
|
<< s.name << endl;
|
|
s.valid = -1;
|
|
ok = false;
|
|
}
|
|
}
|
|
if (ok) log << "ok" << endl;
|
|
else return ok;
|
|
|
|
log << " Checking that properties are continuous at the midpoint temperature... ";
|
|
for (k = 0; k < nsp; k++) {
|
|
Species& s = sp[k];
|
|
|
|
// check continuity at Tmid
|
|
t = s.tmid - dt;
|
|
cp0 = cp(t, s);
|
|
h0 = enthalpy(t, s) + cp0*dt;
|
|
s0 = entropy(t, s) + dt*cp0/t;
|
|
|
|
t = s.tmid + dt;
|
|
cp1 = cp(t, s);
|
|
h1 = enthalpy(t, s) - cp1*dt;
|
|
s1 = entropy(t, s) - cp1*dt/t;
|
|
|
|
if (absval((cp0 - cp1)/cp0) > tol) {
|
|
log << endl << "Warning... species " << s.name
|
|
<< ": discontinuity in Cp at Tmid = "
|
|
<< s.tmid << endl;
|
|
log << "Cp/R (low, high) = " << cp0
|
|
<< ", " << cp1 << endl;
|
|
ok = false;
|
|
}
|
|
|
|
if (absval((h0 - h1)/h0) > tol) {
|
|
log << endl << "Warning... species " << s.name
|
|
<< ": discontinuity in enthalpy at Tmid = "
|
|
<< s.tmid << endl;
|
|
log << "h/R (low, high) = "
|
|
<< h0 << ", " << h1 << endl;
|
|
ok = false;
|
|
}
|
|
|
|
if (absval((s0 - s1)/s0) > tol) {
|
|
log << endl << "Warning... species " << s.name
|
|
<< ": discontinuity in entropy at Tmid = "
|
|
<< s.tmid << endl;
|
|
log << "s/R (low, high) = "
|
|
<< s0 << ", " << s1 << endl;
|
|
ok = false;
|
|
}
|
|
}
|
|
if (ok) log << "ok \n\n\n";
|
|
else log << "\n\n\n";
|
|
|
|
log << " Checking that cp is less that the high-temperature\n"
|
|
<< " limiting value predicted by equipartition of energy.\n";
|
|
log << " Note that this limit does not account for the electronic\n"
|
|
<< " contribution to cp, and so may be violated in some cases."
|
|
<< endl << endl;
|
|
|
|
|
|
for (k = 0; k < nsp; k++) {
|
|
Species& s = sp[k];
|
|
|
|
// check that cp at Tmax is less than the equipartion value
|
|
// This does not include any possible electronic contribution
|
|
|
|
cp0 = cp(s.thigh, s);
|
|
int nel = static_cast<int>(s.elements.size());
|
|
int i;
|
|
double na = 0.0;
|
|
for (i = 0; i < nel; i++)
|
|
if (s.elements[i].name != "E") na += s.elements[i].number;
|
|
int natoms = int(floor(na));
|
|
double cpmax;
|
|
switch (natoms) {
|
|
case 1: cpmax = 2.5; break;
|
|
case 2: cpmax = 4.5; break;
|
|
default: cpmax = 3.0*natoms - 3.0;
|
|
}
|
|
|
|
if (cp0 > cpmax) {
|
|
double over = 100.0*(cp0 - cpmax)/cpmax;
|
|
log << endl << "Warning... species " << s.name
|
|
<< ": cp(Tmax) greater than equipartition value \nby "
|
|
<< over << " percent.";
|
|
if ((natoms > 2) && (cp0 - cpmax < 0.5))
|
|
log << endl << " (if molecule is linear, cp is ok)"
|
|
<< endl;
|
|
}
|
|
}
|
|
|
|
return valid(sp);
|
|
}
|
|
|
|
|
|
|
|
/**
|
|
* Check that all reactions balance.
|
|
* @param speciesData species property dataset used to look up
|
|
* elemental compositions.
|
|
* @param r list of reactions to check
|
|
* @param unbalanced list of integers specifying reaction numbers of
|
|
* unbalanced reactions.
|
|
* @return true if all reactions balance
|
|
* @todo use reaction number stored in reaction object
|
|
*/
|
|
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;
|
|
vector<string> elementSyms;
|
|
unsigned int m;
|
|
|
|
unbalanced.clear();
|
|
map<string, double> atoms;
|
|
|
|
for (int i = 0; i < nrxn; i++)
|
|
{
|
|
atoms.clear();
|
|
int nr = static_cast<int>(r[i].reactants.size());
|
|
int np = static_cast<int>(r[i].products.size());
|
|
int j;
|
|
double stoichCoeff;
|
|
for (j = 0; j < nr; j++)
|
|
{
|
|
rname = r[i].reactants[j].name;
|
|
stoichCoeff = r[i].reactants[j].number;
|
|
vector<Constituent>& elements = speciesData[rname].elements;
|
|
for (m = 0; m < elements.size(); m++)
|
|
{
|
|
atoms[elements[m].name] -= stoichCoeff * elements[m].number;
|
|
}
|
|
}
|
|
for (j = 0; j < np; j++)
|
|
{
|
|
pname = r[i].products[j].name;
|
|
stoichCoeff = r[i].products[j].number;
|
|
vector<Constituent>& elements = speciesData[pname].elements;
|
|
for (m = 0; m < elements.size(); m++)
|
|
{
|
|
atoms[elements[m].name] += stoichCoeff * elements[m].number;
|
|
}
|
|
}
|
|
double atms;
|
|
getMapKeys(atoms, elementSyms);
|
|
for (m = 0; m < elementSyms.size(); m++) {
|
|
atms = atoms[elementSyms[m]];
|
|
if (fabs(atms) > tolerance) {
|
|
//if (atoms[elementSyms[m]] != 0.0) {
|
|
// cout << "Reaction " << i+1 << " has an unbalanced element: "
|
|
// << elementSyms[m] << " "
|
|
// << atoms[elementSyms[m]] << endl;
|
|
unbalanced.push_back(i+1); break;
|
|
}
|
|
}
|
|
}
|
|
return (unbalanced.empty());
|
|
}
|
|
|
|
}
|
|
|
|
|