Cantera  2.0
CKParser.cpp
Go to the documentation of this file.
1 /**
2  * @file CKParser.cpp
3  *
4  */
5 
6 // Copyright 2001 California Institute of Technology
7 
8 #include <numeric>
9 #include <algorithm>
10 #include <fstream>
11 #include <iomanip>
12 #include <math.h>
13 
14 #include "CKParser.h"
15 #include "ckr_utils.h"
16 #include "writelog.h"
17 #include <cstdio>
18 #include <cstdlib>
19 
20 using namespace std;
21 
22 namespace ckr
23 {
24 
25 
26 static string int2s(int n, std::string fmt="%d")
27 {
28  char buf[30];
29  sprintf(buf, fmt.c_str(), n);
30  return string(buf);
31 }
32 
33 /// Exception class for syntax errors.
34 CK_SyntaxError::CK_SyntaxError(std::ostream& f,
35  const std::string& s, int linenum)
36  : m_out(f)
37 {
38  m_msg += "Syntax error: " + s;
39  if (linenum > 0) {
40  m_msg += " (line " + int2s(linenum) + ")\n";
41  }
42 }
43 
44 
45 static int parseGroupString(std::string str, std::vector<std::string>& esyms,
46  vector_int& result);
47 
48 /**
49  * Throw an exception if one of the four lines that must have
50  * 1, 2, 3, or 4 in column 80 do not.
51  */
52 static void illegalThermoLine(std::ostream& f,
53  char n, int linenum = -1)
54 {
55  throw CK_SyntaxError(f, "column 80 must "
56  "contain an integer", linenum);
57 };
58 
59 
60 /**
61  * Throw an exception if number string is bad
62  */
63 static void illegalNumber(std::ostream& f,
64  std::string s, int linenum = -1)
65 {
66  string msg = "illegal number: "+s;
67  throw CK_SyntaxError(f, msg, linenum);
68 };
69 
70 
71 extern void getDefaultAtomicWeights(std::map<std::string,double>& weights);
72 
73 static string d2e(string s)
74 {
75  size_t n;
76  size_t sz = s.size();
77  string r = s;
78  char ch;
79  for (n = 0; n < sz; n++) {
80  ch = s[n];
81  if (ch == 'D') {
82  r[n] = 'E';
83  } else if (ch == 'd') {
84  r[n] = 'e';
85  }
86  }
87  return r;
88 }
89 
90 static double de_atof(std::string s)
91 {
92  string r = d2e(s);
93  //double rval = Cantera::atofCheck(r.c_str());
94  double rval = atof(r.c_str());
95  return rval;
96 }
97 
98 /**
99  * Check validity of the temperatures defining the
100  * temperature ranges for the NASA9 polynomial species thermodynamic
101  * property fits.
102  * @param log log file output stream
103  * @param temp Vector of temperatures
104  */
105 static void checkNASA9Temps(std::ostream& log, vector_fp& temp)
106 {
107  int i;
108  for (i = 1; i < (int) temp.size(); i++) {
109  double tlow = temp[i-1];
110  double thigh = temp[i];
111  if (thigh <= tlow) {
112  string sss = "error reading temperature";
113  throw CK_SyntaxError(log, sss);
114  }
115  }
116 }
117 
118 
119 static double getNumberFromString(std::string s)
120 {
121  bool inexp = false;
122  removeWhiteSpace(s);
123  int sz = static_cast<int>(s.size());
124  char ch;
125  for (int n = 0; n < sz; n++) {
126  ch = s[n];
127  if (!inexp && (ch == 'E' || ch == 'e' || ch == 'D' || ch == 'd')) {
128  inexp = true;
129  } else if (ch == '+' || ch == '-') {
130  if (n > 0 && (s[n-1] != 'E' && s[n-1]
131  != 'e' && s[n-1] != 'd' && s[n-1] != 'D')) {
132  return UNDEF;
133  }
134  } else if (ch != '.' && (ch < '0' || ch > '9')) {
135  return UNDEF;
136  }
137  }
138  return de_atof(s);
139 }
140 
141 /**
142  * Add an element to a species.
143  * @param symbol element symbol
144  * @param atoms number of atoms of this element in the
145  * species (may be non-integral)
146  * @param sp Species object to add element to
147  * @param log log file output stream
148  */
149 static void addElement(std::string symbol, double atoms,
150  Species& sp, std::ostream& log)
151 {
152 
153  if (atoms != 0.0) {
154  Constituent e;
155  e.name = symbol;
156  e.number = atoms;
157  sp.elements.push_back(e);
158  sp.comp[symbol] = atoms;
159  }
160 }
161 
162 
163 /**
164  * Check validity of the three temperatures defining the two
165  * temperature ranges for the NASA polynomial species thermodynamic
166  * property fits.
167  * @param log log file output stream
168  * @param tmin minimum temperature
169  * @param tmid intermediate temperature
170  * @param tmax maximum temperature
171  */
172 static void checkTemps(std::ostream& log, double tmin,
173  double tmid, double tmax)
174 {
175  if (tmin == 0.0 || tmid == 0.0 || tmax == 0.0) {
176  throw CK_SyntaxError(log,
177  "error reading Tmin, Tmid, or Tmax");
178  }
179 }
180 
181 static void getSpecies(std::string s,
182  int n, vector<RxnSpecies>& species, bool debug,
183  std::ostream& log)
184 {
185  removeWhiteSpace(s);
186  // break string into substrings at the '+' characters separating
187  // species symbols
188 
189  bool inplus = true;
190  vector<int> pluses;
191  vector<string> sp;
192  for (int i = n-1; i >= 0; i--) {
193  if (!inplus && s[i] == '+') {
194  pluses.push_back(i);
195  inplus = true;
196  } else if (inplus && s[i] != '+') {
197  inplus = false;
198  }
199  }
200  pluses.push_back(-1);
201  size_t np = pluses.size();
202  size_t loc, nxt;
203  for (size_t nn = 0; nn < np; nn++) {
204  loc = pluses.back();
205  pluses.pop_back();
206  if (nn == np-1) {
207  nxt = s.size();
208  } else {
209  nxt = pluses.back();
210  }
211  sp.push_back(s.substr(loc+1,nxt-loc-1));
212  }
213 
214  string r, num;
215  size_t sz, j, strt=0;
216  RxnSpecies ss;
217  for (size_t nn = 0; nn < sp.size(); nn++) {
218  r = sp[nn];
219  sz = r.size();
220  for (j = 0; j < sz; j++) {
221  if (!((r[j] >= '0' && r[j] <= '9') || r[j] == '.')) {
222  strt = j;
223  break;
224  }
225  }
226  ss.name = r.substr(strt,sz);
227  if (strt == 0) {
228  ss.number = 1.0;
229  } else {
230  ss.number = atof(r.substr(0,strt).c_str());
231  }
232  species.push_back(ss);
233  if (debug) {
234  log << ss.number << " " << ss.name << endl;
235  }
236  }
237 }
238 
239 
240 /**
241  * given a string specifying either the reactant or product side of a
242  * reaction equation, construct a list of Constituent objects
243  * containing the species symbols and stoichiometric coefficients.
244  * @todo allow non-integral stoichiometric coefficients
245  */
246 int getGroups(std::string::const_iterator begin,
247  std::string::const_iterator end, std::vector<std::string>& esyms,
248  std::vector<grouplist_t>& rxngroups)
249 {
250  bool ingroup = false;
251  rxngroups.clear();
252  string g = "";
253  group_t igrp;
254  grouplist_t groups;
255 
256  for (; begin != end; ++begin) {
257  if (*begin == '(') {
258  ingroup = true;
259  g = "";
260  } else if (*begin == ')') {
261  ingroup = false;
262  igrp.clear();
263  if (parseGroupString(g, esyms, igrp) >= 0) {
264  groups.push_back(igrp);
265  } else {
266  return -1;
267  }
268  } else if (*begin == '+') {
269  rxngroups.push_back(groups);
270  groups.clear();
271  } else if (ingroup && *begin != ' ') {
272  g += *begin;
273  }
274  }
275  rxngroups.push_back(groups);
276  return 1;
277 }
278 
279 
280 /**
281  * Constructor. Construct a parser for the specified input file.
282  */
283 CKParser::CKParser(std::istream* infile, const std::string& fname,
284  std::ostream* log) :
285  verbose(true),
286  debug(false),
287  m_line(0),
288  m_nasafmt(false),
289  m_nasa9fmt(false)
290 {
291  m_ckfile = infile;
292  m_ckfilename = fname;
293  m_log = log;
294  m_last_eol = '\n';
295 }
296 
297 
298 
299 // Get a line from the input file, and return it in string s.
300 /*
301  * If the line contains a comment character (!), then return only the
302  * portion preceding it. Non-printing characters are replaced by
303  * spaces.
304  *
305  * The input file is m_ckfile, an istream.
306  *
307  * @param s On return, s contains the line read from the
308  * input file.
309  * @param comment On return, comment contains the text following the
310  * comment character on the line, if any.
311  */
312 void CKParser::getCKLine(std::string& s, std::string& comment)
313 {
314 
315  // Chemkin comment character
316  const char commentChar = '!';
317 
318  // Cantera anti-comment character
319  const char undoCommentChar = '%';
320 
321  // carriage return
322  const char char13 = char(13);
323 
324  // linefeed
325  const char char10 = char(10);
326 
327  istream& f = *m_ckfile;
328 
329  // if putCKLine was called to 'put back' a line, then return this
330  // line, instead of reading a new one
331 
332  if (!m_buf.empty()) {
333  s = m_buf;
334  m_buf = "";
335  comment = m_comment;
336  m_comment = "";
337  return;
338  }
339 
340  // read a line, convert non-printing characters to ' ',
341  // and remove comments
342 
343  comment = "";
344  string line;
345 
346  line = "";
347  char ch = ' ';
348  while (1 > 0) {
349  f.get(ch);
350  if (!f || f.eof()) {
351  break;
352  }
353 
354  // convert tabs to spaces
355  if (ch == '\t') {
356  ch = ' ';
357  }
358 
359  // if an end-of-line character is seen, then break.
360  // Check for all common end-of-line characters.
361  if (ch == char13 || (ch == char10 && (m_last_eol != char13))) {
362  m_last_eol = ch;
363  break;
364  }
365  if (isprint(ch)) {
366  line += ch;
367  }
368  }
369 
370  string::size_type icom = line.find(commentChar);
371 
372  // lines that begin with !% are not comments for Cantera
373  if (icom == 0 && line.size() > 1 && line[1] == undoCommentChar) {
374  line[0] = '%';
375  line[1] = ' ';
376  icom = line.find(commentChar);
377  }
378  int len = static_cast<int>(line.size());
379 
380  for (int i = 0; i < len; i++) if (!isprint(line[i])) {
381  line[i] = ' ';
382  }
383  if (icom != string::npos) {
384  s = line.substr(0, icom);
385  comment = line.substr(icom+1,len-icom-1);
386  } else {
387  s = line;
388  }
389 
390  if (!f || f.eof()) {
391  s = "<EOF>";
392  comment = "";
393  return;
394  }
395  m_line++;
396 }
397 
398 
399 /**
400  *
401  * Put back a line read from the input file. The next call to
402  * getCKLine will return this line.
403  *
404  */
405 
406 void CKParser::putCKLine(std::string& s, std::string& comment)
407 {
408  m_buf = s;
409  m_comment = comment;
410 }
411 
412 
413 bool CKParser::advanceToKeyword(const std::string& kw, const std::string& stop)
414 {
415  string s, c;
416  do {
417  getCKLine(s,c);
418  if (match(s,"<EOF>")) {
419  return false;
420  }
421  if (match(s,kw)) {
422  putCKLine(s,c);
423  return true;
424  }
425  } while (!match(s,stop));
426  putCKLine(s,c);
427  return false;
428 }
429 
430 
431 
432 /**
433  *
434  * Read the element section of the input file, and return
435  * a list of the elements found.
436  *
437  */
438 
439 bool CKParser::readElementSection(elementList& elements)
440 {
441  string s, comment;
442  int firsttok;
443  vector<string> toks;
444 
445  map<string,double> defaultWeights;
446  //ct::ctmap_sd defaultWeights;
447  getDefaultAtomicWeights(defaultWeights);
448 
449  //istream& f = m_ckfile;
450  int ntok;
451 
452  elements.clear();
453  while (1 > 0) {
454 
455 next:
456  if (advanceToKeyword("ELEM", "SPEC")) {
457  firsttok = 1;
458  while (1 > 0) {
459  do {
460  getCKLine(s, comment);
461  getTokens(s, static_cast<int>(s.size()), toks);
462  ntok = static_cast<int>(toks.size());
463  } while (ntok == 0);
464 
465  if (firsttok == 0 && isKeyword(toks[0])) {
466  putCKLine(s,comment);
467  goto next;
468  }
469  for (int i = firsttok; i < ntok; i++) {
470  if (match(toks[i],"END")) {
471  goto next;
472  } else {
473  Element el;
474  string wtstr;
475  el.comment = comment;
476  el.index = static_cast<int>(elements.size());
477  if (extractSlashData(toks[i], el.name, wtstr)) {
478  el.atomicWeight = de_atof(wtstr);
479  el.weightFromDB = false;
480  } else {
481  el.atomicWeight = defaultWeights[capitalize(el.name)];
482  el.weightFromDB = true;
483  }
484  if (el.atomicWeight > 0.0) {
485  el.valid = 1;
486  } else {
487  el.valid = 0;
488  }
489  if (find(elements.begin(),
490  elements.end(), el) < elements.end()) {
491  if (m_log)
492  *m_log << "warning... duplicate element "
493  << el.name << " (ignored)." << endl;
494  } else {
495  elements.push_back(el);
496  }
497  }
498  }
499  firsttok = 0;
500  }
501  } else {
502  if (elements.size() == 0) {
503  *m_log << "no elements found." << endl;
504  return false;
505  } else {
506  return valid(elements);
507  }
508  }
509  }
510  return false;
511 }
512 
513 
514 
515 /**
516  *
517  * Read the SPECIES section of the input file, and return
518  * a list of species names.
519  *
520  */
521 bool CKParser::readSpeciesSection(speciesList& species)
522 {
523  string s, comment;
524  int firsttok;
525  vector<string> toks;
526 
527  int ntok;
528  int nsp = 0;
529 
530  while (1 > 0) {
531 
532 next:
533  if (advanceToKeyword("SPEC", "THER")) {
534  firsttok = 1;
535  while (1 > 0) {
536  do {
537  getCKLine(s, comment);
538  getTokens(s, static_cast<int>(s.size()), toks);
539  ntok = static_cast<int>(toks.size());
540  } while (ntok == 0);
541 
542  if (firsttok == 0 && isKeyword(toks[0])) {
543  putCKLine(s,comment);
544  goto next;
545  }
546  for (int i = firsttok; i < ntok; i++) {
547  if (match(toks[i],"END")) {
548  goto next;
549  } else {
550  Species sp;
551  sp.name = toks[i];
552  if (find(species.begin(), species.end(), sp)
553  < species.end()) {
554  if (m_log)
555  *m_log << "warning... duplicate species "
556  << sp.name << " (ignored)." << endl;
557  } else {
558  nsp++;
559  sp.index = nsp;
560  species.push_back(sp);
561  }
562  }
563  }
564  firsttok = 0;
565  }
566  } else {
567  if (species.size() == 0) {
568  return false;
569  } else {
570  return true;
571  }
572  }
573  }
574  return false;
575 }
576 
577 
578 
579 /**
580  *
581  * Read species data from THERMO section records.
582  *
583  * @param names List of species names (input).
584  * @param species Table of species objects holding data from records
585  * in THERMO section (output).
586  */
587 
588 bool CKParser::readThermoSection(std::vector<std::string>& names,
589  speciesTable& species, vector_fp& temp,
590  int& optionFlag, std::ostream& log)
591 {
592  string s;
593  vector<string> toks;
594 
595  double tmin = -1.0, tmid = -1.0, tmax = -1.0;
596  if (temp.size() == 3) {
597  tmin = temp[0];
598  tmid = temp[1];
599  tmax = temp[2];
600  }
601 
602  int nsp = static_cast<int>(names.size());
603 
604  string comment;
605 
606  // read lines until THERMO section is found. But if EOF reached or
607  // start of REACTIONS section, then there is no THERMO section.
608  do {
609  getCKLine(s,comment);
610  if (match(s,"<EOF>")) {
611  return false;
612  }
613  if (match(s,"REAC")) {
614  putCKLine(s,comment);
615  return false;
616  }
617  } while (!match(s,"THER"));
618 
619  // read the tokens on the THERMO line
620  getTokens(s, static_cast<int>(s.size()), toks);
621  m_nasafmt = false;
622  if (toks.size() >= 2) {
623  unsigned int itt;
624  for (itt = 1; itt < toks.size(); itt++) {
625  if (match(toks[itt],"ALL")) {
626  optionFlag = NoThermoDatabase;
627  } else if (match(toks[itt],"NO_TMID")) {
628  m_nasafmt = true;
629  log << "\nOption 'NO_TMID' specified. Default "
630  "midpoint temperature\n";
631  log << "will be used for all species.\n\n";
632  } else if (match(toks[itt], "NASA9")) {
633  m_nasa9fmt = true;
634  log << "Option NASA9 specified: Use new "
635  "nasa input file format\n\n";
636  } else if (match(toks[itt], "NASA")) {
637  m_nasa9fmt = false;
638  log << "Option NASA specified: Use old "
639  "nasa input file format\n\n";
640  } else throw CK_SyntaxError(log,
641  "unrecognized THERMO option.", m_line);
642  }
643  }
644 
645  // if "THERMO ALL" specified, or if optionFlag is set to HasTempRange,
646  // then the next line must contain the default temperatures
647  // for the database.
648 
649  if (optionFlag == NoThermoDatabase || optionFlag == HasTempRange) {
650  getCKLine(s, comment);
651  getTokens(s, static_cast<int>(s.size()), toks);
652  if (m_nasa9fmt) {
653  //
654  // For NASA9 polynomials, the format is
655  // t1 t2 t3 t4 date
656  // when there are 3 temperature regions
657  //
658  size_t nreg = toks.size() - 2;
659  if (nreg >= 1) {
660  temp.resize(nreg+1);
661  for (size_t i = 0; i <= nreg; i++) {
662  temp[i] = de_atof(toks[i]);
663  }
664  string defaultDate = toks[nreg+1];
665  } else {
666  throw CK_SyntaxError(log, "Default temp region card is bad", m_line);
667  }
668  if (verbose) {
669  log.flags(ios::showpoint | ios::fixed);
670  log.precision(2);
671  log << endl << " Default # of temperature regions: "
672  << nreg << endl;
673  log << " ";
674  for (size_t i = 0; i <= nreg; i++) {
675  log << temp[i] << " ";
676  }
677  log << endl;
678  }
679  checkNASA9Temps(log, temp);
680  } else {
681  //
682  // For NASA polynomials, the format is
683  // tlow tmid thigh
684  // There are always 2 temperature regions
685  //
686  if (toks.size() >= 3) {
687  tmin = de_atof(toks[0]);
688  tmid = de_atof(toks[1]);
689  tmax = de_atof(toks[2]);
690  }
691 
692  if (verbose) {
693  log.flags(ios::showpoint | ios::fixed);
694  log.precision(2);
695  log << endl << " default Tlow, Tmid, Thigh: " << tmin << " "
696  << tmid << " " << tmax << endl;
697  }
698  checkTemps(log, tmin, tmid, tmax);
699  temp.clear();
700  temp.push_back(tmin);
701  temp.push_back(tmid);
702  temp.push_back(tmax);
703  }
704  }
705 
706  /// XXXX BRANCH TO THE DIFFERENT THERMO READERS HERE
707 
708  // Check to see that we expect to be reading a NASA formatted file
709  if (m_nasa9fmt) {
710  bool ok = readNASA9ThermoSection(names, species, temp,
711  optionFlag, log);
712  if (!ok) {
713  throw CK_SyntaxError(log,
714  "In NASA parser. However, we expect a NASA9 file format",
715  -1);
716  }
717  return ok;
718  }
719 
720  // now read in all species records that have names in list 'names'
721 
722  bool getAllSpecies = (nsp > 0 && match(names[0],"<ALL>"));
723  if (getAllSpecies) {
724  names.clear();
725  }
726 
727  map<string, int> dup; // used to check for duplicate THERMO records
728  bool already_read;
729 
730  while (1 > 0) {
731  if (nsp == 0) {
732  break;
733  }
734  already_read = false;
735 
736  Species spec;
737  readThermoRecord(spec);
738 
739  if (spec.name == "<END>") {
740  break;
741  }
742 
743  // check for duplicate thermo data
744  if (dup[spec.name] == 2) {
745  log << "Warning: more than one THERMO record for "
746  << "species " << spec.name << endl;
747  log << "Record at line " << m_line
748  << " of " << m_ckfilename << " ignored." << endl;
749  already_read = true;
750  }
751  dup[spec.name] = 2;
752 
753  if (!already_read && (getAllSpecies
754  || (find(names.begin(), names.end(), spec.name)
755  < names.end()))) {
756 
757  if (spec.tmid == 0.0) {
758  spec.tmid = tmid;
759  log << "Warning: default Tmid used for species "
760  << spec.name << endl;
761  if (spec.tmid < 0.0) {
762  log << "Error: no default Tmid has been entered!"
763  << endl;
764  }
765  }
766  species[spec.name] = spec;
767 
768  if (verbose) {
769  log << endl << "found species " << spec.name;
770  log << " at line " << m_line
771  << " of " << m_ckfilename;
772  writeSpeciesData(log, spec);
773  }
774  checkTemps(log, spec.tlow, spec.tmid, spec.thigh);
775  if (getAllSpecies) {
776  names.push_back(spec.name);
777  nsp = static_cast<int>(names.size());
778  } else {
779  nsp--;
780  }
781  }
782  }
783  return true;
784 }
785 
786 
787 /**
788  *
789  * Read one 4-line species definition record in NASA format.
790  *
791  */
792 
793 void CKParser::readThermoRecord(Species& sp)
794 {
795  string s;
796  string numstr;
797  double cf;
798 
799  // look for line 1, but if a keyword is found first or the end of
800  // the file is reached, return "<END>" as the species name
801  string comment;
802  do {
803  getCKLine(s, comment);
804  if (isKeyword(s) || match(s, "<EOF>")) {
805  sp.name = "<END>";
806  putCKLine(s, comment);
807  return;
808  }
809  } while ((s.size() < 80) || (s[79] != '1'));
810 
811  // next 4 lines must be the NASA-format lines without intervening
812  // comments.
813 
814  //------------- line 1 ---------------------------
815 
816  if (s[79] != '1') {
817  illegalThermoLine(*m_log, s[79], m_line);
818  }
819 
820  // extract the species name and the id string (date)
821  string nameid = s.substr(0,24);
822  vector<string> toks;
823  getTokens(nameid, static_cast<int>(nameid.size()), toks);
824  sp.name = toks[0];
825  sp.id = "";
826  unsigned int j;
827  for (j = 1; j < toks.size(); j++) {
828  if (j > 1) {
829  sp.id += ' ';
830  }
831  sp.id += toks[j];
832  }
833  int iloc;
834  string elementSym;
835  double atoms;
836  int i;
837 
838  // elemental composition (first 4 elements)
839  for (i = 0; i < 4; i++) {
840  elementSym = "";
841  iloc = 24 + 5*i;
842  if (s[iloc] != ' ') {
843  if (s[iloc+1] != ' ') {
844  elementSym = s.substr(iloc,2);
845  } else {
846  elementSym = s.substr(iloc,1);
847  }
848  } else if (s[iloc+1] != ' ') {
849  elementSym = s.substr(iloc+1,1);
850  }
851  atoms = de_atof(s.substr(iloc+2,3));
852  addElement(elementSym, atoms, sp, *m_log);
853  }
854 
855  // single-character phase descriptor
856  sp.phase = s[44];
857 
858  // low, high, and mid temperatures
859  sp.tlow = de_atof(s.substr(45,10));
860  sp.thigh = de_atof(s.substr(55,10));
861 
862  if (!m_nasafmt) {
863  sp.tmid = de_atof(s.substr(65,8));
864 
865  // fifth element, if any
866  elementSym = "";
867  if (s[73] != ' ') {
868  elementSym += s[73];
869  }
870  if (s[74] != ' ') {
871  elementSym += s[74];
872  }
873  atoms = de_atof(s.substr(75,3));
874  addElement(elementSym, atoms, sp, *m_log);
875 
876  // additional elements, if any
877  elementSym = "";
878  int loc = 80;
879  while (loc < (int)(s.size()-9)) {
880  elementSym = "";
881  if (s[loc] != ' ') {
882  elementSym += s[loc];
883  }
884  if (s[loc+1] != ' ') {
885  elementSym += s[loc+1];
886  }
887  atoms = de_atof(s.substr(loc+2,8));
888  addElement(elementSym, atoms, sp, *m_log);
889  loc += 10;
890  }
891  }
892 
893  //-------------- line 2 ----------------------------
894 
895  getCKLine(s, comment);
896  if (s[79] != '2') {
897  illegalThermoLine(*m_log, s[79], m_line);
898  }
899  for (i = 0; i < 5; i++) {
900  numstr = s.substr(i*15, 15);
901  cf = getNumberFromString(numstr);
902  if (cf == UNDEF) {
903  illegalNumber(*m_log, numstr, m_line);
904  }
905  sp.highCoeffs.push_back(cf);
906  }
907 
908  //-------------- line 3 ----------------------------
909 
910  getCKLine(s, comment);
911  if (s[79] != '3') {
912  illegalThermoLine(*m_log, s[79], m_line);
913  }
914  for (i = 0; i < 2; i++) {
915  numstr = s.substr(i*15, 15);
916  cf = getNumberFromString(numstr);
917  if (cf == UNDEF) {
918  illegalNumber(*m_log, numstr, m_line);
919  }
920  sp.highCoeffs.push_back(cf);
921  }
922  for (i = 2; i < 5; i++) {
923  numstr = s.substr(i*15, 15);
924  cf = getNumberFromString(numstr);
925  if (cf == UNDEF) {
926  illegalNumber(*m_log, numstr, m_line);
927  }
928  sp.lowCoeffs.push_back(cf);
929  }
930 
931  //--------------- line 4 ----------------------------
932 
933  getCKLine(s, comment);
934  if (s[79] != '4') {
935  illegalThermoLine(*m_log, s[79], m_line);
936  }
937  for (i = 0; i < 4; i++) {
938  numstr = s.substr(i*15, 15);
939  cf = getNumberFromString(numstr);
940  if (cf == UNDEF) {
941  illegalNumber(*m_log, numstr, m_line);
942  }
943  sp.lowCoeffs.push_back(cf);
944  }
945  sp.valid = 1;
946 }
947 
948 
949 
950 void CKParser::missingAuxData(const std::string& kw)
951 {
952  throw CK_SyntaxError(*m_log, kw +
953  " keyword must be followed by slash-delimited data.", m_line);
954 }
955 
956 
957 /**
958  * Parse the REACTION section of the input file, and return
959  * a list of Reaction objects and the units.
960  */
961 bool CKParser::readReactionSection(const std::vector<std::string>& speciesNames,
962  vector<string>& elementNames, reactionList& reactions,
963  ReactionUnits& units)
964 {
965  string s, comment;
966  vector<string> toks;
967  int nRxns = 0;
968 
969  vector<string> rc, pr;
970  vector_int c;
971 
972  // advance to the beginning of the REACTION section
973  do {
974  getCKLine(s, comment);
975  if (match(s, "<EOF>")) {
976  return false;
977  }
978  } while (!match(s,"REAC"));
979 
980 
981  // look for units specifications
982 
983  getTokens(s, static_cast<int>(s.size()), toks);
984  string tok;
985  units.ActEnergy = Cal_per_Mole;
986  units.Quantity = Moles;
987  unsigned int ir;
988  for (ir = 1; ir < toks.size(); ir++) {
989  tok = toks[ir];
990  if (match(tok,"CAL/MOLE")) {
991  units.ActEnergy = Cal_per_Mole;
992  } else if (match(tok,"KCAL/MOLE")) {
993  units.ActEnergy = Kcal_per_Mole;
994  } else if (match(tok,"JOULES/MOLE")) {
995  units.ActEnergy = Joules_per_Mole;
996  } else if (match(tok,"KJOULES/MOLE")) {
997  units.ActEnergy = Kjoules_per_Mole;
998  } else if (match(tok,"KELVINS")) {
999  units.ActEnergy = Kelvin;
1000  } else if (match(tok,"EVOLTS")) {
1001  units.ActEnergy = Electron_Volts;
1002  } else if (match(tok,"MOLES")) {
1003  units.Quantity = Moles;
1004  } else if (match(tok,"MOLECULES")) {
1005  units.Quantity = Molecules;
1006  }
1007  }
1008 
1009  Reaction rxn;
1010 
1011  vector<string> cm;
1012  bool ok = true;
1013 
1014  if (debug) {
1015  *m_log << "CKParser::readReactions ---> DEBUG MODE" << endl;
1016  }
1017 
1018  while (1 > 0) {
1019 
1020  // skip blank or comment lines
1021  do {
1022  getCKLine(s, comment);
1023  cm.push_back(comment);
1024  } while (s == "" && comment[0] != '%');
1025 
1026 #undef DEBUG_LINE
1027 #ifdef DEBUG_LINE
1028  *m_log << "Line: " << s << endl;
1029 #endif
1030  // end of REACTION section or EOF
1031  /// @todo does this handle case of 1 reaction correctly?
1032  if (isKeyword(s) || s == "<EOF>") {
1033  if (nRxns > 0) {
1034  rxn.number = nRxns;
1035  reactions.push_back(rxn);
1036  //rxn.comment.clear();
1037  }
1038  if (nRxns > 0) {
1039  return ok;
1040  }
1041  return false;
1042  }
1043 
1044  // rxn line
1045  //string::size_type eqloc;
1046  string sleft, sright;
1047  bool auxDataLine, metaDataLine;
1048 
1049 
1050  // if the line contains an '=', it is the start of a new reaction.
1051  // In this case, add the previous reaction to the output list,
1052  // increment the number of reactions, and start processing the
1053  // new reaction.
1054 
1055  size_t eqloc = s.find_first_of("=");
1056  metaDataLine = false;
1057  auxDataLine = false;
1058 
1059  // look for a metadata line
1060  if (s[0] == '%') {
1061  metaDataLine = true;
1062  if (eqloc > 0 && eqloc < s.size()) {
1063  int ierr, ierp;
1064  vector<grouplist_t> rg, pg;
1065  s[eqloc] = ' ';
1066  ierr = getGroups(s.begin(), s.begin() + eqloc,
1067  elementNames, rg);
1068  ierp = getGroups(s.begin() + eqloc, s.end(),
1069  elementNames, pg);
1070  unsigned int nr =
1071  static_cast<unsigned int>(rxn.reactants.size());
1072  unsigned int nratoms = 0;
1073  for (unsigned int ij = 0; ij < nr; ij++) {
1074  nratoms += int(rxn.reactants[ij].number);
1075  }
1076  if (rg.size() != nratoms)
1077  throw CK_SyntaxError(*m_log,
1078  " groups not specified for all reactants", m_line);
1079  else if (ierr < 0)
1080  throw CK_SyntaxError(*m_log,
1081  " error in reactant group specification", m_line);
1082  for (unsigned int ir = 0; ir < nr; ir++) {
1083  rxn.reactants[ir].groups = rg[ir];
1084  }
1085  unsigned int np =
1086  static_cast<unsigned int>(rxn.products.size());
1087  unsigned int npatoms = 0;
1088  for (unsigned int ik = 0; ik < np; ik++) {
1089  npatoms += int(rxn.products[ik].number);
1090  }
1091  if (pg.size() != npatoms)
1092  throw CK_SyntaxError(*m_log,
1093  " groups not specified for all products", m_line);
1094  else if (ierp < 0)
1095  throw CK_SyntaxError(*m_log,
1096  " error in product group specification", m_line);
1097  for (unsigned int ip = 0; ip < np; ip++) {
1098  rxn.products[ip].groups = pg[ip];
1099  }
1100  }
1101  }
1102 
1103  else if (eqloc != string::npos && eqloc < s.size()) {
1104  if (nRxns > 0) {
1105  rxn.number = nRxns;
1106  reactions.push_back(rxn);
1107  }
1108  nRxns++;
1109  rxn = Reaction();
1110  rxn.comment = cm;
1111  cm.clear();
1112  if (debug) {
1113  *m_log << "Parsing reaction " << nRxns << endl;
1114  }
1115  } else {
1116  auxDataLine = true;
1117  }
1118  if (comment != "") {
1119  rxn.lines.push_back(s+'!'+comment);
1120  } else {
1121  rxn.lines.push_back(s);
1122  }
1123 
1124  if (!auxDataLine && !metaDataLine) {
1125 
1126  // depending on the form of the 'equals' symbol,
1127  // determine whether the reaction is reversible or
1128  // irreversible, and separate it into strings for
1129  // each side.
1130 
1131  if (eqloc = s.find("<=>"), eqloc != string::npos) {
1132  rxn.isReversible = true;
1133  sleft = s.substr(0, eqloc);
1134  sright = s.substr(eqloc+3,1000);
1135  } else if (eqloc = s.find("=>"), eqloc != string::npos) {
1136  rxn.isReversible = false;
1137  sleft = s.substr(0, eqloc);
1138  sright = s.substr(eqloc+2,1000);
1139  } else if (eqloc = s.find("="), eqloc != string::npos) {
1140  rxn.isReversible = true;
1141  sleft = s.substr(0, eqloc);
1142  sright = s.substr(eqloc+1,1000);
1143  } else throw CK_SyntaxError(*m_log,
1144  "expected <=>, =>, or =", m_line);
1145 
1146  if (debug) {
1147  *m_log << s << endl;
1148  if (rxn.isReversible) {
1149  *m_log << "Reaction is reversible." << endl;
1150  } else {
1151  *m_log << "Reaction is irreversible." << endl;
1152  }
1153  }
1154 
1155  string::size_type mloc, mloc2;
1156 
1157  // process reactants
1158  if (debug) {
1159  *m_log << "Processing reactants..." << sleft << endl;
1160  }
1161  removeWhiteSpace(sleft);
1162  if (debug) *m_log << "After removing white space: "
1163  << sleft << endl;
1164  rxn.isFalloffRxn = false;
1165 
1166  string sm, mspecies;
1167 
1168  mloc = sleft.find("(+");
1169  if (mloc != string::npos) {
1170  sm = sleft.substr(mloc+2, 1000);
1171  mloc2 = sm.find(")");
1172  if (mloc2 != string::npos) {
1173  mspecies = sm.substr(0,mloc2);
1174  rxn.isFalloffRxn = true;
1175  rxn.type = Falloff;
1176  sleft = sleft.substr(0, mloc);
1177  if (mspecies == "M" || mspecies == "m") {
1178  rxn.thirdBody = "M";
1179  } else {
1180  rxn.thirdBody = mspecies;
1181  }
1182  if (debug) {
1183  *m_log << "Falloff reaction. Third body = "
1184  << rxn.thirdBody << endl;
1185  }
1186  } else throw CK_SyntaxError(*m_log,
1187  "missing )", m_line);
1188  }
1189 
1190  else if ((mloc = sleft.find("+M"), mloc != string::npos) ||
1191  (mloc = sleft.find("+m"), mloc != string::npos)) {
1192 
1193  if (static_cast<int>(mloc) ==
1194  static_cast<int>(sleft.size()) - 2) {
1195  rxn.isThreeBodyRxn = true;
1196  rxn.type = ThreeBody;
1197  sleft = sleft.substr(0, mloc);
1198  rxn.thirdBody = "M";
1199  if (debug) {
1200  *m_log << "Three-body reaction." << endl;
1201  }
1202  } else if (debug) {
1203  *m_log << "Reactant string contains +M or +m, but \n"
1204  << "not last two characters of string: "
1205  << "\"" << sleft << "\"\n"
1206  << "NOT a three-body reaction." << endl;
1207  }
1208  }
1209 
1210  getSpecies(sleft.c_str(),static_cast<int>(sleft.size()),
1211  rxn.reactants, debug, *m_log);
1212  int ir = static_cast<int>(rxn.reactants.size());
1213  for (int iir = 0; iir < ir; iir++) {
1214  if (find(speciesNames.begin(), speciesNames.end(),
1215  rxn.reactants[iir].name) >= speciesNames.end())
1216  throw CK_SyntaxError(*m_log,
1217  "undeclared reactant species "
1218  +rxn.reactants[iir].name, m_line);
1219  }
1220 
1221 
1222  // process Arrhenius coefficients
1223  getTokens(sright, static_cast<int>(sright.size()), toks);
1224  int ntoks = static_cast<int>(toks.size());
1225  if (ntoks < 3) {
1226  throw CK_SyntaxError(*m_log,
1227  "expected 3 Arrhenius parameters", m_line);
1228  }
1229  rxn.kf.A = de_atof(toks[ntoks - 3]);
1230  rxn.kf.n = de_atof(toks[ntoks - 2]);
1231  rxn.kf.E = de_atof(toks[ntoks - 1]);
1232 
1233  // 2/10/03: allow negative prefactor but print a warning
1234  if (rxn.kf.A < 0.0)
1235  *m_log << "Warning: negative prefactor at line "
1236  << m_line << endl;
1237  //throw CK_SyntaxError(*m_log, "negative prefactor", m_line);
1238 
1239  if (debug) {
1240  *m_log << "Processing products..." << sright << endl;
1241  }
1242  sright = sright.substr(0, sright.find(toks[ntoks - 3]) - 1);
1243  if (debug) *m_log << "After removing Arrhenius parameters, "
1244  << "\nproduct string = " << sright << endl;
1245 
1246  removeWhiteSpace(sright);
1247  if (debug) *m_log << "After removing white space: "
1248  << sright << endl;
1249  mloc = sright.find("(+");
1250  if (mloc != string::npos) {
1251  sm = sright.substr(mloc+2, 1000);
1252  mloc2 = sm.find(")");
1253  if (mloc2 != string::npos) {
1254  mspecies = sm.substr(0,mloc2);
1255 
1256  if (rxn.type == ThreeBody)
1257  throw CK_SyntaxError(*m_log,
1258  "mismatched +M or (+M)", m_line);
1259 
1260  rxn.isFalloffRxn = true;
1261  rxn.type = Falloff;
1262  if (debug) {
1263  *m_log << "Falloff reaction. Third body = "
1264  << rxn.thirdBody << endl;
1265  }
1266  } else throw CK_SyntaxError(*m_log,
1267  "missing )", m_line);
1268 
1269  sright = sright.substr(0, mloc);
1270 
1271  if (mspecies == "M" || mspecies == "m") {
1272  rxn.thirdBody = "M";
1273  } else {
1274  if (rxn.thirdBody != mspecies)
1275  throw CK_SyntaxError(*m_log,
1276  "mismatched third body", m_line);
1277  rxn.thirdBody = mspecies;
1278  }
1279  }
1280 
1281  else if ((mloc = sright.find("+M"), mloc != string::npos) ||
1282  (mloc = sright.find("+m"), mloc != string::npos)) {
1283 
1284  if (static_cast<int>(mloc) ==
1285  static_cast<int>(sright.size()) - 2) {
1286 
1287  if (rxn.type == Falloff)
1288  throw CK_SyntaxError(*m_log,
1289  "mismatched +M or (+M)", m_line);
1290  rxn.isThreeBodyRxn = true;
1291  rxn.thirdBody = "M";
1292  sright = sright.substr(0, mloc);
1293  if (debug) {
1294  *m_log << "Three-body reaction." << endl;
1295  }
1296  } else if (debug) {
1297  *m_log << "Product string contains +M or +m, but \n"
1298  << "not last two characters of string: "
1299  << "\"" << sright << "\"\n"
1300  << "NOT a three-body reaction." << endl;
1301  }
1302  }
1303  getSpecies(sright.c_str(),static_cast<int>(sright.size()),
1304  rxn.products, debug, *m_log);
1305  int ip = static_cast<int>(rxn.products.size());
1306  for (int iip = 0; iip < ip; iip++) {
1307  if (find(speciesNames.begin(), speciesNames.end(),
1308  rxn.products[iip].name) >= speciesNames.end())
1309  throw CK_SyntaxError(*m_log,
1310  "undeclared product species "+rxn.products[iip].name, m_line);
1311  }
1312  }
1313 
1314  // auxiliary data line
1315  else if (auxDataLine) {
1316 
1317  bool hasAuxData;
1318  string name, data;
1319  map<string, int> kwindex;
1320  while (1 > 0) {
1321 
1322  hasAuxData = extractSlashData(s, name, data);
1323  if (!hasAuxData && name == "") {
1324  break;
1325  }
1326 
1327  // check for duplicate keyword
1328  if (kwindex[name]) {
1329  throw CK_SyntaxError(*m_log,
1330  "duplicate auxiliary data keyword "
1331  + name, m_line);
1332  } else {
1333  kwindex[name] = 1;
1334  }
1335 
1336 
1337  // low-pressure rate coefficient for falloff rxn
1338 
1339  if (match(name,"LOW")) {
1340  vector<string> klow;
1341  rxn.type = Falloff;
1342  if (hasAuxData) {
1343  getTokens(data, static_cast<int>(data.size()), klow);
1344  if (klow.size() != 3) {
1345  throw CK_SyntaxError(*m_log,
1346  "expected 3 low-pressure Arrhenius parameters", m_line);
1347  }
1348  rxn.kf_aux.A = de_atof(klow[0]);
1349  rxn.kf_aux.n = de_atof(klow[1]);
1350  rxn.kf_aux.E = de_atof(klow[2]);
1351  } else {
1352  missingAuxData("LOW");
1353  }
1354  }
1355 
1356 
1357  // falloff parameters
1358 
1359  else if (match(name,"TROE")) {
1360  vector<string> falloff;
1361  if (kwindex["SRI"] > 0) {
1362  throw CK_SyntaxError(*m_log,
1363  "cannot specify both SRI and TROE", m_line);
1364  }
1365 
1366  if (hasAuxData) {
1367  getTokens(data, static_cast<int>(data.size()), falloff);
1368  int nf = static_cast<int>(falloff.size());
1369  double ff;
1370  rxn.falloffType = Troe;
1371  for (int jf = 0; jf < nf; jf++) {
1372  ff = de_atof(falloff[jf]);
1373  rxn.falloffParameters.push_back(ff);
1374  }
1375  } else {
1376  missingAuxData("TROE");
1377  }
1378  }
1379 
1380  else if (match(name,"SRI")) {
1381  vector<string> falloff;
1382  if (kwindex["TROE"] > 0) {
1383  throw CK_SyntaxError(*m_log,
1384  "cannot specify both SRI and TROE", m_line);
1385  }
1386  if (hasAuxData) {
1387  getTokens(data, static_cast<int>(data.size()), falloff);
1388  int nf = static_cast<int>(falloff.size());
1389  rxn.falloffType = SRI;
1390  double ff;
1391  for (int jf = 0; jf < nf; jf++) {
1392  ff = de_atof(falloff[jf]);
1393  rxn.falloffParameters.push_back(ff);
1394  }
1395  } else {
1396  missingAuxData("SRI");
1397  }
1398  }
1399 
1400 
1401 
1402  // reverse rate coefficient
1403 
1404  else if (match(name,"REV")) {
1405  vector<string> krev;
1406  if (!rxn.isReversible) {
1407  throw CK_SyntaxError(*m_log,
1408  "reverse rate parameters can only be "
1409  "specified for reversible reactions", m_line);
1410  }
1411  if (hasAuxData) {
1412  getTokens(data, static_cast<int>(data.size()), krev);
1413  if (krev.size() != 3) {
1414  throw CK_SyntaxError(*m_log,
1415  "expected 3 Arrhenius parameters", m_line);
1416  }
1417  rxn.krev.A = de_atof(krev[0]);
1418  rxn.krev.n = de_atof(krev[1]);
1419  rxn.krev.E = de_atof(krev[2]);
1420  } else {
1421  missingAuxData("REV");
1422  }
1423  }
1424 
1425 
1426  else if (match(name,"DUP")) {
1427  rxn.isDuplicate = true;
1428  }
1429 
1430  else if (match(name,"END")) {
1431  string c = "";
1432  putCKLine(name,c);
1433  break;
1434  }
1435 
1436 
1437  // Landau-Teller reaction rate parameters
1438 
1439  else if (match(name,"LT")) {
1440  vector<string> bc;
1441  rxn.kf.type = LandauTeller;
1442  if (hasAuxData) {
1443  getTokens(data, static_cast<int>(data.size()), bc);
1444  rxn.kf.B = de_atof(bc[0]);
1445  rxn.kf.C = de_atof(bc[1]);
1446  } else {
1447  missingAuxData("LT");
1448  }
1449  }
1450 
1451  else if (match(name,"RLT")) {
1452  vector<string> bc;
1453  rxn.krev.type = LandauTeller;
1454  if (hasAuxData) {
1455  getTokens(data, static_cast<int>(data.size()), bc);
1456  rxn.krev.B = de_atof(bc[0]);
1457  rxn.krev.C = de_atof(bc[1]);
1458  } else {
1459  missingAuxData("RLT");
1460  }
1461  }
1462 
1463 
1464  // chem activation reactions
1465  else if (match(name,"HIGH")) {
1466  vector<string> khigh;
1467  rxn.type = ChemAct;
1468  if (hasAuxData) {
1469  getTokens(data, static_cast<int>(data.size()), khigh);
1470  rxn.kf_aux.A = de_atof(khigh[0]);
1471  rxn.kf_aux.n = de_atof(khigh[1]);
1472  rxn.kf_aux.E = de_atof(khigh[2]);
1473  } else {
1474  missingAuxData("HIGH");
1475  }
1476  }
1477 
1478  else if (match(name,"FORD")) {
1479  vector<string> nmord;
1480  if (hasAuxData) {
1481  getTokens(data, static_cast<int>(data.size()),
1482  nmord);
1483  rxn.fwdOrder[nmord[0]] = de_atof(nmord[1]);
1484  } else {
1485  missingAuxData("FORD");
1486  }
1487  }
1488 
1489  else if (find(speciesNames.begin(), speciesNames.end(), name)
1490  < speciesNames.end()) {
1491  if (hasAuxData) {
1492  if (rxn.thirdBody == name || rxn.thirdBody == "M") {
1493  rxn.e3b[name] = de_atof(data);
1494  } else if (rxn.thirdBody == "<none>") {
1495  *m_log << "Error in reaction " << nRxns
1496  << ": third-body collision efficiencies cannot be specified"
1497  << " for this reaction type." << endl;
1498  throw CK_SyntaxError(*m_log,
1499  "third-body efficiency error", m_line);
1500  } else {
1501  *m_log << "Reaction " << nRxns << ": illegal species in enhanced "
1502  << "efficiency specification. Species = "
1503  << name << " rxn.thirdBody = "
1504  << rxn.thirdBody << endl;
1505  throw CK_SyntaxError(*m_log,
1506  "third-body efficiency error", m_line);
1507  }
1508  } else {
1509  missingAuxData(name);
1510  }
1511  } else {
1512  Reaction::auxdata vals;
1513  vector<string> toks;
1514  getTokens(data, static_cast<int>(data.size()), toks);
1515  int ntoks = static_cast<int>(toks.size());
1516  for (int itok = 0; itok < ntoks; itok++) {
1517  vals.push_back(de_atof(toks[itok]));
1518  }
1519  rxn.otherAuxData[name] = vals;
1520  }
1521  }
1522  }
1523  }
1524  return false;
1525 }
1526 
1527 
1528 
1529 int parseGroupString(std::string str, std::vector<std::string>& esyms, group_t& result)
1530 {
1531  bool inSymbol=true;
1532  string s = str + '-';
1533  int i;
1534  string num, sym;
1535  int eindx;
1536  string::const_iterator begin = s.begin();
1537  string::const_iterator end = s.end();
1538  vector<string>::iterator e;
1539  result.resize(static_cast<size_t>(esyms.size()),0);
1540  for (; begin != end; ++begin) {
1541 
1542  // new element
1543  if (*begin == '-') {
1544  e = find(esyms.begin(), esyms.end(), sym);
1545  if (e == esyms.end()) {
1546  return -1;
1547  }
1548  eindx = static_cast<int>(e - esyms.begin());
1549  if (num != "") {
1550  i = atoi(num.c_str());
1551  } else {
1552  i = 1;
1553  }
1554  result[eindx] = i;
1555  sym = "";
1556  num = "";
1557  inSymbol = true;
1558  } else if (isdigit(*begin)) {
1559  inSymbol = false;
1560  num += *begin;
1561  } else if (isalpha(*begin) && inSymbol) {
1562  sym += *begin;
1563  }
1564  }
1565  return 1;
1566 }
1567 
1568 } // ckr namespace
1569 
1570 
1571 
1572 
1573 
Generated by   doxygen 1.8.2