Cantera  2.0
CKReader.cpp
Go to the documentation of this file.
1 /**
2  * @file CKReader.cpp
3  *
4  */
5 
6 // Copyright 2001 California Institute of Technology
7 
8 #include <fstream>
9 #include <string>
10 using namespace std;
11 
12 #include "CKParser.h"
13 #include "CKReader.h"
14 #include "thermoFunctions.h"
15 
16 #include <cstring>
17 #include <cstdlib>
18 #include <ctime>
19 #include <iomanip>
20 
21 #include "writelog.h"
22 #include <cstdio>
23 #include "ckr_defs.h"
24 
25 //#include "cantera/base/global.h"
26 //#define APP Cantera::Application
27 
28 namespace ckr
29 {
30 
31 /**
32  * read and optionally validate an input file in Chemkin format.
33  * @param inputFile path to the input file
34  * @param thermoDatabase path to the species database file
35  * @param logfile path to the file where log messages should be written
36  * @return true if no errors were encountered, false otherwise
37  */
38 bool CKReader::read(const std::string& inputFile, const std::string& thermoDatabase,
39  const std::string& logfile)
40 {
41 
42  clock_t t0, t1;
43 
44  t0 = clock();
45 
46  ifstream ckinfile(inputFile.c_str());
47  ofstream log(logfile.c_str());
48  try {
49  // construct a parser for the input file
50  CKParser parser(&ckinfile, inputFile, &log);
51  parser.verbose = verbose;
52  parser.debug = debug;
53 
54  // write header information to the log file
55  struct tm* newtime;
56  time_t aclock;
57  time(&aclock); /* Get time in seconds */
58  newtime = localtime(&aclock); /* Convert time to struct tm form */
59 
60 
61  log << "CKReader version 1.0" << endl
62  << "http://www.cantera.org" << endl << endl
63  << asctime(newtime) << endl
64  << setw(20) << "input file: "
65  << setw(30) << inputFile << endl;
66 
67  if (thermoDatabase != "")
68  log << setw(20) << "species database: "
69  << setw(30) << thermoDatabase << endl;
70 
71  if (!validate)
72  log << endl << "*************** Warning ***************" << endl
73  << " mechanism validation disabled" << endl
74  << "*****************************************" << endl;
75 
76  if (debug) {
77  log << "*** DEBUG MODE ***" << endl;
78  } else {
79  log << "debugging disabled." << endl;
80  }
81 
82  //----------- process ELEMENT section ----------------------
83 
84  bool elok = parser.readElementSection(elements);
85  int nel = static_cast<int>(elements.size());
86  vector<string> elementSymbols;
87  for (int j = 0; j < nel; j++) {
88  elementSymbols.push_back(elements[j].name);
89  }
90 
91  if (verbose) {
92  log.flags(ios::showpoint);
93  log.precision(6);
94  log.width(0);
95 
96  log << endl << newTask("reading elements") << endl;
97 
98  // write summary to log file
99  for (int i = 0; i < nel; i++) {
100  log << i+1 << ". " << pad(elements[i].name,2) << " ";
101  double wt = elements[i].atomicWeight;
102  if (wt == 0.0) {
103  log << "<error>";
104  } else {
105  log << wt;
106  }
107  if (!elements[i].weightFromDB) {
108  log << " (specified)";
109  }
110  if (elements[i].comment != "") {
111  log << " ! " << elements[i].comment;
112  }
113  log << endl;
114  }
115  }
116  log << "\nread " << nel << " elements." << endl;
117 
118  if (!elok) {
119  log << "\nerrors were encountered." << endl;
120  return false;
121  }
122  if (nel == 0) {
123  return false;
124  }
125 
126 
127  //------------ process SPECIES section ------------------------
128 
129  vector<string> speciesSymbols;
130  bool spok = parser.readSpeciesSection(species);
131  int nsp = static_cast<int>(species.size());
132 
133  if (verbose) {
134  log << newTask("reading species") << endl;
135  }
136 
137  for (int i = 0; i < nsp; i++) {
138  Species& s = species[i];
139  if (verbose) {
140  log << i+1 << ". " << s.name << endl;
141  }
142  speciesSymbols.push_back(s.name);
143  }
144  log << "\nread " << nsp << " species." << endl;
145 
146  if (!spok) {
147  log << "\nerrors were encountered." << endl;
148  return false;
149  }
150  if (nsp == 0) {
151  return false;
152  }
153 
154 
155  //------------- process THERMO section -------------------------
156 
157  if (verbose) {
158  log << newTask("looking up species definitions") << endl;
159  }
160 
161  // if a thermo database is specified, get the default Tmin, Tmid, Tmax
162  vector_fp temp;
163 
164  if (thermoDatabase != "") {
165 
166  if (verbose) log << "reading default temperature ranges from "
167  << thermoDatabase << endl;
168 
169  ifstream thermofile(thermoDatabase.c_str());
170  CKParser thermReader(&thermofile, thermoDatabase, &log);
171  thermReader.verbose = verbose;
172  thermReader.debug = debug;
173  int dbflag = HasTempRange;
174  vector<string> dummy;
175  thermReader.readThermoSection(dummy, speciesData, temp, dbflag, log);
176  }
177 
178 
179  bool hasthermo = parser.advanceToKeyword("THERM","REAC");
180 
181  int k, optionFlag = 0;
182  int undefined = static_cast<int>(species.size());
183  string nm;
184  vector<string> undef;
185  bool allsp = (speciesSymbols[0] == "<ALL>");
186  if (hasthermo &&
187  parser.readThermoSection(speciesSymbols,
188  speciesData, temp, optionFlag, log)) {
189  if (allsp) {
190  nsp = static_cast<int>(speciesData.size());
191  for (k = 0; k < nsp; k++) {
192  Species s;
193  s.name = speciesSymbols[k];
194  species.push_back(s);
195  }
196  }
197  undefined = 0;
198  for (k = 0; k < nsp; k++) {
199  nm = species[k].name;
200  species[k] = speciesData[species[k].name];
201  if (species[k].name == "<empty>") {
202  undefined++;
203  undef.push_back(nm);
204  species[k].name = nm;
205  }
206  }
207  int localdefs = nsp - undefined;
208  if (localdefs > 0 && verbose) log << "found definitions for "
209  << localdefs
210  << " of "
211  << nsp
212  << " species in the input file. "
213  << endl;
214  } else {
215  undef = speciesSymbols;
216  if (verbose) {
217  log << "no THERMO section in input file." << endl;
218  }
219  }
220 
221  if (undefined > 0 && thermoDatabase != ""
222  && optionFlag != NoThermoDatabase) {
223 
224  if (verbose) log << "searching external database "
225  << thermoDatabase << " for species definitions..."
226  << endl;
227 
228  ifstream thermofile(thermoDatabase.c_str());
229  CKParser thermoReader(&thermofile, thermoDatabase, &log);
230  thermoReader.verbose = verbose;
231  thermoReader.debug = debug;
232  int dbflag = HasTempRange;
233  thermoReader.readThermoSection(undef, speciesData, temp, dbflag, log);
234  undefined = 0;
235  if (allsp) {
236  species.clear();
237  nsp = static_cast<int>(speciesData.size());
238  for (k = 0; k < nsp; k++) {
239  Species s;
240  s.name = undef[k];
241  species.push_back(s);
242  }
243  }
244  for (int k = 0; k < nsp; k++) {
245  if (species[k].valid == 0) {
246  nm = species[k].name;
247  species[k] = speciesData[species[k].name];
248  if (species[k].name == "<empty>") {
249  undefined++;
250  species[k].name = nm;
251  }
252  }
253  }
254  }
255 
256  if (validate && !validateSpecies(log)) {
257  //Cantera::setError("read","error in species");
258  return false;
259  }
260 
261 
262  //------------- process REACTIONS section -------------------------
263 
264  if (verbose) {
265  log << newTask("reading reactions") << endl;
266  }
267 
268  ckinfile.close();
269  ifstream ckinfile2(inputFile.c_str());
270 
271  // construct a new parser for the input file
272  CKParser parser2(&ckinfile2, inputFile, &log);
273  parser2.verbose = verbose;
274  parser2.debug = debug;
275 
276  parser2.readReactionSection(speciesSymbols, elementSymbols, reactions, units);
277  log << "\nread " << static_cast<int>(reactions.size())
278  << " reactions." << endl;
279 
280  bool rxnok = true;
281  if (validate) {
282  rxnok = rxnok && validateReactions(log);
283  }
284  bool writeok = true;
285  if (verbose || validate) {
286  writeok = writeReactions(log);
287  }
288  rxnok = rxnok && writeok;
289  if (!rxnok) {
290  return false;
291  }
292 
293  log << "\nSuccess... ";
294  t1 = clock();
295  log << "elapsed CPU time = "
296  << double(t1 - t0)/CLOCKS_PER_SEC
297  << " s" << endl;
298  if (!validate) {
299  log << "*** no validation performed ***" << endl;
300  }
301  }
302 
303  catch (CK_Exception& e) {
304  log << e.errorMessage() << endl;
305  //Cantera::setError("CKReader::read",e.errorMessage());
306  return false;
307  } catch (std::exception& e) {
308  log << "an exception was raised in CKReader:\n";
309  log << e.what() << std::endl;
310  return false;
311  }
312 
313  return true;
314 }
315 
316 
317 /// print a summary of all reactions to the log file
318 bool CKReader::writeReactions(std::ostream& log)
319 {
320 
321  bool ok = true;
322  // int ns = species.size();
323  int nrxns = static_cast<int>(reactions.size());
324  log.flags(ios::unitbuf);
325  log.precision(6);
326 
327  log << endl;
328  for (int n = 0; n < nrxns; n++) {
329  Reaction& r = reactions[n];
330 
331  log << "reaction " << r.number << endl;
332  log << " ";
333  printReactionEquation(log, r);
334  log << endl;
335 
336  // rate coefficient
337  if (r.isFalloffRxn) {
338 
339  log << " high P rate coeff: ";
340  ok = ok && writeRateCoeff(r.kf, log) ;
341 
342  log << " low P rate coeff: ";
343  ok = ok && writeRateCoeff(r.kf_aux, log);
344 
345  ok = ok && writeFalloff(r.falloffType, r.falloffParameters, log);
346 
347  } else {
348  log << " rate coeff: ";
349  ok = ok && writeRateCoeff(r.kf, log);
350  }
351  if (r.isReversible && r.krev.A > 0) {
352  log << " reverse rate coeff: ";
353  ok = ok && writeRateCoeff(r.krev, log);
354  }
355  int ne = static_cast<int>(r.e3b.size());
356 
357  if (ne > 0) {
358  vector<string> enhSpecies;
359  getMapKeys(r.e3b, enhSpecies);
360  log << " enhanced collision efficiencies:" << endl;
361  log << " ";
362  for (int nn = 0; nn < ne; nn++) {
363  log << enhSpecies[nn] << " "
364  << r.e3b[enhSpecies[nn]];
365  if (nn < ne-1) {
366  log << ", ";
367  }
368  }
369  log << endl;
370  }
371  if (r.isDuplicate) log
372  << " declared duplicate reaction. See reaction "
373  << r.duplicate << "." << endl;
374  log << endl;
375  }
376  return ok;
377 }
378 
379 
380 
381 /// validate the species
382 bool CKReader::validateSpecies(std::ostream& log)
383 {
384  size_t nel = elements.size();
385  size_t nsp = species.size();
386  double tol;
387 
388  log << newTask("validating species");
389 
390  // check for undeclared elements
391  vector<string> esyms;
392 
393  log << " checking that all species have been defined... ";
394  for (size_t k = 0; k < nsp; k++) {
395  Species& s = species[k];
396  if (s.valid == 0) {
397  log << endl << " species " << s.name << " undefined ";
398  s.valid = -1;
399  }
400  }
401  if (valid(species)) {
402  log << "OK" << endl;
403  } else {
404  log << endl;
405  return false;
406  }
407 
408  log << " checking that all species elements have been declared... ";
409  for (size_t k = 0; k < nsp; k++) {
410  Species& s = species[k];
411 
412  getMapKeys(s.comp, esyms);
413  size_t nm = esyms.size();
414 
415  for (size_t m = 0; m < nm; m++) {
416  size_t j;
417  for (j = 0; j < nel; j++) {
418  if (esyms[m] == elements[j].name) {
419  break;
420  }
421  }
422  if (j == nel) {
423  log << endl << " species "
424  << s.name << ": undeclared element "
425  << esyms[m];
426  s.valid = -1;
427  }
428  }
429  }
430  if (valid(species)) {
431  log << "OK" << endl;
432  } else {
433  log << endl;
434  return false;
435  }
436 
437  log << " checking consistency of species thermo data... ";
438  tol = 0.01;
439  if (checkThermo(log, species, tol)) {
440  log << "OK" << endl;
441  } else {
442  log << endl;
443  return false;
444  }
445  return true;
446 
447 }
448 
449 
450 /// validate the reactions
451 bool CKReader::validateReactions(std::ostream& log)
452 {
453 
454  bool ok = true;
455  // int ns = species.size();
456  int nrxns = static_cast<int>(reactions.size());
457 
458  vector<int> unbal;
459  log << "checking that all reactions balance...";
460  if (checkBalance(log, speciesData, reactions, unbal)) {
461  log << " OK" << endl;
462  } else {
463  int nu = static_cast<int>(unbal.size());
464  for (int iu = 0; iu < nu; iu++) {
465  log << " error... reaction " << unbal[iu]
466  << " does not balance" << endl;
467  }
468  ok = false;
469  }
470 
471  log << "checking for duplicate reactions...";
472 
473  for (int nn = 0; nn < nrxns; nn++) {
474  Reaction& r1 = reactions[nn];
475  for (int mm = nn + 1; mm < nrxns; mm++) {
476  Reaction& r2 = reactions[mm];
477  if (r1 == r2) {
478  r1.duplicate = mm + 1;
479  r2.duplicate = nn + 1;
480  if (!r1.isDuplicate || !r2.isDuplicate) {
481  log << endl << " error... undeclared duplicate reactions: "
482  << nn + 1 << " and " << mm + 1;
483  ok = false;
484  } else {
485  log << endl << " declared duplicate reactions: "
486  << nn + 1
487  << " and " << mm + 1;
488  }
489  }
490  }
491  }
492  if (ok) {
493  log << "...OK" << endl;
494  }
495 
496  return ok;
497 }
498 
499 
500 /**
501  * Check the thermodynamic property parameterizations for all species.
502  * The following are verified:
503  * - The heat capacity is positive throughout the full temperature range;
504  * - The entropy at Tmin is positive;
505  * - The heat capacity, entropy, and enthalpy evaluated at Tmid using
506  * both the high and low polynomial coefficients are the same to within
507  * relative error tol
508  * - The heat capacity at Tmax is not greater than the equipartition limit
509  * for the number of atoms in the molecule
510  */
511 bool checkThermo(std::ostream& log, speciesList& sp, double tol)
512 {
513  const double dt = 0.0001;
514  double t, cp0, h0, s0, cp1, h1, s1;
515  int nsp = static_cast<int>(sp.size());
516  const int n_points = 20;
517 
518  int k;
519  bool ok = true;
520  for (k = 0; k < nsp; k++) {
521  Species& s = sp[k];
522 
523  if (s.valid <= 0) {
524  ok = false;
525  log << endl << "species " << s.name
526  << " contains an error." << endl;
527  }
528  if (!ok) {
529  return false;
530  }
531  }
532 
533  log << endl << " Checking that cp/R is positive... ";
534 
535  for (k = 0; k < nsp; k++) {
536  Species& s = sp[k];
537 
538  // check that cp is positive
539 
540  cp0 = 0.0;
541  for (int j = 0; j < n_points; j++) {
542  t = j*(s.thigh - s.tlow)/(n_points - 1) + s.tlow;
543 
544  cp0 = cp(t, s);
545  if (cp0 < 0.0) {
546  log << endl << " error... Cp/R < 0 at T = " << t
547  << " for species " << s.name << endl;
548  s.valid = -1;
549  ok = false;
550  }
551  }
552  }
553  if (ok) {
554  log << "ok" << endl;
555  } else {
556  return ok;
557  }
558 
559 
560  // check that S(tlow) > 0
561 
562  log << " Checking that the species entropies are positive... ";
563  for (k = 0; k < nsp; k++) {
564  Species& s = sp[k];
565  if (entropy(s.tlow, s) <= 0.0) {
566  log << endl << " error... negative entropy for species "
567  << s.name << endl;
568  s.valid = -1;
569  ok = false;
570  }
571  }
572  if (ok) {
573  log << "ok" << endl;
574  } else {
575  return ok;
576  }
577 
578  log << " Checking that properties are continuous at the midpoint temperature... ";
579  for (k = 0; k < nsp; k++) {
580  Species& s = sp[k];
581 
582  // check continuity at Tmid
583  t = s.tmid - dt;
584  cp0 = cp(t, s);
585  h0 = enthalpy(t, s) + cp0*dt;
586  s0 = entropy(t, s) + dt*cp0/t;
587 
588  t = s.tmid + dt;
589  cp1 = cp(t, s);
590  h1 = enthalpy(t, s) - cp1*dt;
591  s1 = entropy(t, s) - cp1*dt/t;
592 
593  if (absval((cp0 - cp1)/cp0) > tol) {
594  log << endl << "Warning... species " << s.name
595  << ": discontinuity in Cp at Tmid = "
596  << s.tmid << endl;
597  log << "Cp/R (low, high) = " << cp0
598  << ", " << cp1 << endl;
599  ok = false;
600  }
601 
602  if (absval((h0 - h1)/h0) > tol) {
603  log << endl << "Warning... species " << s.name
604  << ": discontinuity in enthalpy at Tmid = "
605  << s.tmid << endl;
606  log << "h/R (low, high) = "
607  << h0 << ", " << h1 << endl;
608  ok = false;
609  }
610 
611  if (absval((s0 - s1)/s0) > tol) {
612  log << endl << "Warning... species " << s.name
613  << ": discontinuity in entropy at Tmid = "
614  << s.tmid << endl;
615  log << "s/R (low, high) = "
616  << s0 << ", " << s1 << endl;
617  ok = false;
618  }
619  }
620  if (ok) {
621  log << "ok \n\n\n";
622  } else {
623  log << "\n\n\n";
624  }
625 
626  log << " Checking that cp is less that the high-temperature\n"
627  << " limiting value predicted by equipartition of energy.\n";
628  log << " Note that this limit does not account for the electronic\n"
629  << " contribution to cp, and so may be violated in some cases."
630  << endl << endl;
631 
632 
633  for (k = 0; k < nsp; k++) {
634  Species& s = sp[k];
635 
636  // check that cp at Tmax is less than the equipartion value
637  // This does not include any possible electronic contribution
638 
639  cp0 = cp(s.thigh, s);
640  int nel = static_cast<int>(s.elements.size());
641  int i;
642  double na = 0.0;
643  for (i = 0; i < nel; i++)
644  if (s.elements[i].name != "E") {
645  na += s.elements[i].number;
646  }
647  int natoms = int(floor(na));
648  double cpmax;
649  switch (natoms) {
650  case 1:
651  cpmax = 2.5;
652  break;
653  case 2:
654  cpmax = 4.5;
655  break;
656  default:
657  cpmax = 3.0*natoms - 3.0;
658  }
659 
660  if (cp0 > cpmax) {
661  double over = 100.0*(cp0 - cpmax)/cpmax;
662  log << endl << "Warning... species " << s.name
663  << ": cp(Tmax) greater than equipartition value \nby "
664  << over << " percent.";
665  if ((natoms > 2) && (cp0 - cpmax < 0.5))
666  log << endl << " (if molecule is linear, cp is ok)"
667  << endl;
668  }
669  }
670 
671  return valid(sp);
672 }
673 
674 
675 
676 /**
677  * Check that all reactions balance.
678  * @param speciesData species property dataset used to look up
679  * elemental compositions.
680  * @param r list of reactions to check
681  * @param unbalanced list of integers specifying reaction numbers of
682  * unbalanced reactions.
683  * @return true if all reactions balance
684  * @todo use reaction number stored in reaction object
685  */
686 bool checkBalance(std::ostream& f, speciesTable& speciesData,
687  reactionList& r, std::vector<int>& unbalanced, double tolerance)
688 {
689  int nrxn = static_cast<int>(r.size());
690  string rname, pname;
691  vector<string> elementSyms;
692  unsigned int m;
693 
694  unbalanced.clear();
695  map<string, double> atoms;
696 
697  for (int i = 0; i < nrxn; i++) {
698  atoms.clear();
699  int nr = static_cast<int>(r[i].reactants.size());
700  int np = static_cast<int>(r[i].products.size());
701  int j;
702  double stoichCoeff;
703  for (j = 0; j < nr; j++) {
704  rname = r[i].reactants[j].name;
705  stoichCoeff = r[i].reactants[j].number;
706  vector<Constituent>& elements = speciesData[rname].elements;
707  for (m = 0; m < elements.size(); m++) {
708  atoms[elements[m].name] -= stoichCoeff * elements[m].number;
709  }
710  }
711  for (j = 0; j < np; j++) {
712  pname = r[i].products[j].name;
713  stoichCoeff = r[i].products[j].number;
714  vector<Constituent>& elements = speciesData[pname].elements;
715  for (m = 0; m < elements.size(); m++) {
716  atoms[elements[m].name] += stoichCoeff * elements[m].number;
717  }
718  }
719  double atms;
720  getMapKeys(atoms, elementSyms);
721  for (m = 0; m < elementSyms.size(); m++) {
722  atms = atoms[elementSyms[m]];
723  if (fabs(atms) > tolerance) {
724  //if (atoms[elementSyms[m]] != 0.0) {
725  // cout << "Reaction " << i+1 << " has an unbalanced element: "
726  // << elementSyms[m] << " "
727  // << atoms[elementSyms[m]] << endl;
728  unbalanced.push_back(i+1);
729  break;
730  }
731  }
732  }
733  return (unbalanced.empty());
734 }
735 
736 }
737 
738