Cantera  2.0
importKinetics.cpp
Go to the documentation of this file.
1 /**
2  * @file importKinetics.cpp
3  * Declarations of global routines for the importing
4  * of kinetics data from XML files (see \ref inputfiles).
5  *
6  * This file contains routines which are global routines, i.e.,
7  * not part of any object. These routine take as input, ctml
8  * pointers to data, and pointers to %Cantera objects. The purpose
9  * of these routines is to initialize the %Cantera objects with data
10  * from the ctml tree structures.
11  */
12 // Copyright 2002 California Institute of Technology
13 
15 #include "cantera/thermo/mix_defs.h"
16 #include <time.h>
17 #include <memory>
18 
19 // Cantera includes
29 #include "cantera/base/global.h"
31 
32 #include "cantera/base/xml.h"
33 #include "cantera/base/ctml.h"
34 
35 #include <cstdlib>
36 
37 using namespace ctml;
38 using namespace std;
39 
40 namespace Cantera
41 {
42 
43 ReactionRules::ReactionRules() :
44  skipUndeclaredSpecies(false),
45  skipUndeclaredThirdBodies(false),
46  allowNegativeA(false)
47 {
48 }
49 
50 //! these are all used to check for duplicate reactions
51 class rxninfo
52 {
53 public:
54  //! Net stoichiometric coefficients for each reaction
55  std::vector< std::map<int, doublereal> > m_rdata;
56 
57  //! string name (i.e. the reaction equation)
58  std::vector<std::string> m_eqn;
59 
60  //! Indicates whether each reaction is marked "duplicate"
61  std::vector<int> m_dup;
62 
63  //! Number of reactants in each reaction
64  std::vector<size_t> m_nr;
65 
66  //! Indicates "type" of each reaction (see reaction_defs.h)
67  std::vector<int> m_typ;
68 
69  //! Indicates whether each reaction is reversible
70  std::vector<bool> m_rev;
71 
72  //! Map of (vector indicating participating species) to reaction numbers
73  //! Used to speed up duplicate reaction checks.
74  std::map<std::vector<char>, std::vector<size_t> > m_participants;
75 
76  bool installReaction(int i, const XML_Node& r, Kinetics& kin,
77  std::string default_phase, ReactionRules& rule,
78  bool validate_rxn) ;
79 };
80 
81 /*
82  * Check a reaction to see if the elements balance.
83  */
85  const ReactionData& rdata, doublereal errorTolerance)
86 {
87  doublereal kstoich;
88 
89  map<string, double> bal, balr, balp;
90  bal.clear();
91  balp.clear();
92  balr.clear();
93  //cout << "checking " << rdata.equation << endl;
94  size_t np = rdata.products.size();
95 
96  // iterate over the products
97  for (size_t index = 0; index < np; index++) {
98  size_t kp = rdata.products[index]; // index of the product in 'kin'
99  size_t n = kin.speciesPhaseIndex(kp); // phase this product belongs to
100  size_t klocal = kp - kin.kineticsSpeciesIndex(0,n); // index within this phase
101  kstoich = rdata.pstoich[index]; // product stoichiometric coeff
102  const ThermoPhase& ph = kin.speciesPhase(kp);
103  for (size_t m = 0; m < ph.nElements(); m++) {
104  bal[ph.elementName(m)] += kstoich*ph.nAtoms(klocal,m);
105  balp[ph.elementName(m)] += kstoich*ph.nAtoms(klocal,m);
106  //cout << "product species " << ph.speciesName(klocal) << " has " << ph.nAtoms(klocal,m)
107  // << " atoms of " << ph.elementName(m) << " and kstoich = " << kstoich << endl;
108  }
109  }
110  for (size_t index = 0; index < rdata.reactants.size(); index++) {
111  size_t kr = rdata.reactants[index];
112  size_t n = kin.speciesPhaseIndex(kr);
113  //klocal = kr - kin.start(n);
114  size_t klocal = kr - kin.kineticsSpeciesIndex(0,n);
115  kstoich = rdata.rstoich[index];
116  const ThermoPhase& ph = kin.speciesPhase(kr);
117  for (size_t m = 0; m < ph.nElements(); m++) {
118  bal[ph.elementName(m)] -= kstoich*ph.nAtoms(klocal,m);
119  balr[ph.elementName(m)] += kstoich*ph.nAtoms(klocal,m);
120  //cout << "reactant species " << ph.speciesName(klocal) << " has " << ph.nAtoms(klocal,m)
121  // << " atoms of " << ph.elementName(m) << " and kstoich = " << kstoich << endl;
122  }
123  }
124 
125  map<string, double>::iterator b = bal.begin();
126  string msg = "\n\tElement Reactants Products";
127  bool ok = true;
128  doublereal err, elemsum;
129  for (; b != bal.end(); ++b) {
130  elemsum = fabs(balr[b->first]) + fabs(balp[b->first]);
131  if (elemsum > 0.0) {
132  err = fabs(b->second/elemsum);
133  if (err > errorTolerance) {
134  ok = false;
135  msg += "\n\t"+b->first+" "+ fp2str(balr[b->first])
136  +" "+ fp2str(balp[b->first]);
137  }
138  }
139  }
140  if (!ok) {
141  msg = "The following reaction is unbalanced:\n\t"
142  + rdata.equation + "\n" + msg + "\n";
143  throw CanteraError("checkRxnElementBalance",msg);
144  }
145 }
146 
147 
148 /**
149  * Get the reactants or products of a reaction. The information
150  * is returned in the spnum, stoich, and order vectors. The
151  * length of the vectors is the number of different types of
152  * reactants or products found for the reaction.
153  *
154  * Input
155  * --------
156  * rxn -> xml node pointing to the reaction element
157  * in the xml tree.
158  * kin -> Reference to the kinetics object to install
159  * the information into.
160  * rp = 1 -> Go get the reactants for a reaction
161  * -1 -> Go get the products for a reaction
162  * default_phase = String name for the default phase
163  * to loop up species in.
164  * Output
165  * -----------
166  * spnum = vector of species numbers found.
167  * Length is number of reactants or products.
168  * stoich = stoichiometric coefficient of the reactant or product
169  * Length is number of reactants or products.
170  * order = Order of the reactant and product in the reaction
171  * rate expression
172  * rules = If we fail to find a species, we will throw an error
173  * if rule != 1. If rule = 1, we simply return false,
174  * allowing the calling routine to skip this reaction
175  * and continue.
176  */
177 bool getReagents(const XML_Node& rxn, Kinetics& kin, int rp,
178  std::string default_phase, std::vector<size_t>& spnum,
179  vector_fp& stoich, vector_fp& order,
180  const ReactionRules& rules)
181 {
182 
183  string rptype;
184 
185  /*
186  * The id of reactants and products are kept in child elements
187  * of reaction, named "reactants" and "products". We search
188  * the xml tree for these children based on the value of rp,
189  * and store the xml element pointer here.
190  */
191  if (rp == 1) {
192  rptype = "reactants";
193  } else {
194  rptype = "products";
195  }
196  const XML_Node& rg = rxn.child(rptype);
197 
198  /*
199  * The species and stoichiometric coefficient for the species
200  * are stored as a colon separated pair. Get all of these
201  * pairs in the reactions/products object.
202  */
203  vector<string> key, val;
204  getPairs(rg, key, val);
205 
206  /*
207  * Loop over each of the pairs and process them
208  */
209  doublereal ord, stch;
210  string ph, sp;
211  map<string, size_t> speciesMap;
212  for (size_t n = 0; n < key.size(); n++) {
213  sp = key[n]; // sp is the string name for species
214  ph = "";
215  /*
216  * Search for the species in the kinetics object using the
217  * member function kineticsSpeciesIndex(). We will search
218  * for the species in all phases defined in the kinetics operator.
219  */
220  size_t isp = kin.kineticsSpeciesIndex(sp);
221  if (isp == npos) {
222  if (rules.skipUndeclaredSpecies) {
223  return false;
224  } else {
225  throw CanteraError("getReagents",
226  "Undeclared reactant or product species "+sp);
227  return false;
228  }
229  }
230 
231  /*
232  * For each reagent, we store the the species number, isp
233  * the stoichiometric coefficient, val[n], and the order
234  * species in the reaction rate expression. We assume mass
235  * action kinetics here, but will modify this below for
236  * specified species.
237  */
238  spnum.push_back(isp);
239  stch = atof(val[n].c_str());
240  stoich.push_back(stch);
241  ord = doublereal(stch);
242  order.push_back(ord);
243  //cout << key[n] << " " << isp << " " << stch << endl;
244 
245  /*
246  * Needed to process reaction orders below.
247  */
248  speciesMap[sp] = order.size();
249  }
250 
251  /*
252  * Check to see if reaction orders have been specified.
253  */
254  if (rp == 1 && rxn.hasChild("order")) {
255  vector<XML_Node*> ord;
256  rxn.getChildren("order",ord);
257  doublereal forder;
258  for (size_t nn = 0; nn < ord.size(); nn++) {
259  const XML_Node& oo = *ord[nn];
260  string sp = oo["species"];
261  size_t loc = speciesMap[sp];
262  if (loc == 0)
263  throw CanteraError("getReagents",
264  "reaction order specified for non-reactant: "
265  +sp);
266  forder = fpValue(oo());
267  if (forder < 0.0) {
268  throw CanteraError("getReagents",
269  "reaction order must be non-negative");
270  }
271  // replace the stoichiometric coefficient
272  // stored above in 'order' with the specified
273  // reaction order
274  order[loc-1] = forder;
275  }
276  }
277  return true;
278 }
279 
280 
281 /**
282  * getArrhenius() parses the xml element called Arrhenius.
283  * The Arrhenius expression is
284  * \f[ k = A T^(b) exp (-E_a / RT). \f]
285  */
286 static void getArrhenius(const XML_Node& node, int& highlow,
287  doublereal& A, doublereal& b, doublereal& E)
288 {
289 
290  if (node["name"] == "k0") {
291  highlow = 0;
292  } else {
293  highlow = 1;
294  }
295  /*
296  * We parse the children for the A, b, and E components.
297  */
298  A = getFloat(node, "A", "toSI");
299  b = getFloat(node, "b");
300  E = getFloat(node, "E", "actEnergy");
301  E /= GasConstant;
302 }
303 
304 /**
305  * getStick() processes the XML element called Stick that specifies
306  * the sticking coefficient reaction. This routine will
307  * translate the sticking coefficient value into a "normal"
308  * rate constant for the surface reaction.
309  *
310  * Output
311  * -----------
312  * Output is the normal Arrhenius expressions for a surface
313  * reaction rate constant.
314  *
315  * A - units such that rate of rxn has kmol/m^2/s when
316  * A is multiplied by activity concentrations of
317  * reactants in the normal manner.
318  * n - unitless
319  * E - Units 1/Kelvin
320  */
321 static void getStick(const XML_Node& node, Kinetics& kin,
322  ReactionData& r, doublereal& A, doublereal& b, doublereal& E)
323 {
324  size_t nr = r.reactants.size();
325  size_t k, klocal, not_surf = 0;
326  size_t np = 0;
327  doublereal f = 1.0;
328  doublereal order;
329  /*
330  * species is the name of the special reactant whose surface
331  * flux rate will be calculated.
332  * isp = species # in the local phase
333  * ispKinetics = species # in the kinetics object
334  * ispPhaseIndex = phase # of the special species
335  */
336  string spname = node["species"];
337  ThermoPhase& th = kin.speciesPhase(spname);
338  size_t isp = th.speciesIndex(spname);
339  size_t ispKinetics = kin.kineticsSpeciesIndex(spname);
340  size_t ispPhaseIndex = kin.speciesPhaseIndex(ispKinetics);
341 
342  doublereal ispMW = th.molecularWeights()[isp];
343  doublereal sc;
344 
345  // loop over the reactants
346  for (size_t n = 0; n < nr; n++) {
347  k = r.reactants[n];
348  order = r.rorder[n]; // stoich coeff
349 
350  // get the phase species k belongs to
351  np = kin.speciesPhaseIndex(k);
352  const ThermoPhase& p = kin.thermo(np);
353 
354  // get the local index of species k in this phase
355  klocal = p.speciesIndex(kin.kineticsSpeciesName(k));
356 
357  // if it is a surface species, divide f by the standard
358  // concentration for this species, in order to convert
359  // from concentration units used in the law of mass action
360  // to coverages used in the sticking probability
361  // expression
362  if (p.eosType() == cSurf || p.eosType() == cEdge) {
363  sc = p.standardConcentration(klocal);
364  f /= pow(sc, order);
365  }
366  // Otherwise:
367  else {
368  // We only allow one species to be in the phase
369  // containing the special sticking coefficient
370  // species.
371  if (ispPhaseIndex == np) {
372  not_surf++;
373  }
374  // Other bulk phase species on the other side
375  // of ther interface are treated like surface
376  // species.
377  else {
378  sc = p.standardConcentration(klocal);
379  f /= pow(sc, order);
380  }
381  }
382  }
383  if (not_surf != 1) {
384  throw CanteraError("getStick",
385  "reaction probabilities can only be used in "
386  "reactions with exactly 1 gas/liquid species.");
387  }
388 
389  doublereal cbar = sqrt(8.0*GasConstant/(Pi*ispMW));
390  A = 0.25 * getFloat(node, "A", "toSI") * cbar * f;
391  b = getFloat(node, "b") + 0.5;
392  E = getFloat(node, "E", "actEnergy");
393  E /= GasConstant;
394 }
395 
396 static void getCoverageDependence(const XML_Node& node,
397  thermo_t& surfphase, ReactionData& rdata)
398 {
399  vector<XML_Node*> cov;
400  node.getChildren("coverage", cov);
401  size_t k, nc = cov.size();
402  doublereal e;
403  string spname;
404  if (nc > 0) {
405  for (size_t n = 0; n < nc; n++) {
406  const XML_Node& cnode = *cov[n];
407  spname = cnode["species"];
408  k = surfphase.speciesIndex(spname);
409  rdata.cov.push_back(doublereal(k));
410  rdata.cov.push_back(getFloat(cnode, "a"));
411  rdata.cov.push_back(getFloat(cnode, "m"));
412  e = getFloat(cnode, "e", "actEnergy");
413  rdata.cov.push_back(e/GasConstant);
414  }
415  }
416 }
417 
418 
419 //! Get falloff parameters for a reaction.
420 /*!
421  * This routine reads the falloff XML node and extracts parameters into a
422  * vector of doubles
423  *
424  *
425  * @verbatim
426  <falloff type="Troe"> 0.5 73.2 5000. 9999. </falloff>
427  @endverbatim
428 */
429 static void getFalloff(const XML_Node& f, ReactionData& rdata)
430 {
431  string type = f["type"];
432  vector<string> p;
433  getStringArray(f,p);
434  vector_fp c;
435  int np = static_cast<int>(p.size());
436  for (int n = 0; n < np; n++) {
437  c.push_back(fpValue(p[n]));
438  }
439  if (type == "Troe") {
440  if (np == 4) {
441  rdata.falloffType = TROE4_FALLOFF;
442  } else if (np == 3) {
443  rdata.falloffType = TROE3_FALLOFF;
444  } else {
445  throw CanteraError("getFalloff()", "Troe parameterization is specified by number of pararameters, "
446  + int2str(np) + ", is not equal to 3 or 4");
447  }
448  } else if (type == "SRI") {
449  if (np == 5) {
450  rdata.falloffType = SRI5_FALLOFF;
451  if (c[2] < 0.0) {
452  throw CanteraError("getFalloff()", "SRI5 m_c parameter is less than zero: " + fp2str(c[2]));
453  }
454  if (c[3] < 0.0) {
455  throw CanteraError("getFalloff()", "SRI5 m_d parameter is less than zero: " + fp2str(c[3]));
456  }
457  } else if (np == 3) {
458  rdata.falloffType = SRI3_FALLOFF;
459  if (c[2] < 0.0) {
460  throw CanteraError("getFalloff()", "SRI3 m_c parameter is less than zero: " + fp2str(c[2]));
461  }
462  } else {
463  throw CanteraError("getFalloff()", "SRI parameterization is specified by number of pararameters, "
464  + int2str(np) + ", is not equal to 3 or 5");
465  }
466  }
467  rdata.falloffParameters = c;
468 }
469 
470 /**
471  * Get the enhanced collision efficiencies. It is assumed that the
472  * reaction mechanism is homogeneous, so that all species belong
473  * to phase(0) of 'kin'.
474  */
475 static void getEfficiencies(const XML_Node& eff, Kinetics& kin,
476  ReactionData& rdata, const ReactionRules& rules)
477 {
478  // set the default collision efficiency
479  rdata.default_3b_eff = fpValue(eff["default"]);
480 
481  vector<string> key, val;
482  getPairs(eff, key, val);
483  string nm;
484  string phse = kin.thermo(0).id();
485  for (size_t n = 0; n < key.size(); n++) { // ; bb != ee; ++bb) {
486  nm = key[n];// bb->first;
487  size_t k = kin.kineticsSpeciesIndex(nm, phse);
488  if (k != npos) {
489  rdata.thirdBodyEfficiencies[k] = fpValue(val[n]); // bb->second;
490  } else if (!rules.skipUndeclaredThirdBodies) {
491  throw CanteraError("getEfficiencies", "Encountered third-body "
492  "efficiency for undefined species \"" + nm + "\"\n"
493  "while adding reaction " + int2str(rdata.number+1) + ".");
494  }
495  }
496 }
497 
498 /*
499  * Extract the rate coefficient for a reaction from the xml node, kf.
500  * kf should point to a XML element named "rateCoeff".
501  * rdata is the partially filled ReactionData object for the reaction.
502  * This function will fill in more fields in the ReactionData object.
503  *
504  * @param kf Reference to the XML Node named rateCoeff
505  */
506 void getRateCoefficient(const XML_Node& kf, Kinetics& kin,
507  ReactionData& rdata, const ReactionRules& rules)
508 {
509  if (rdata.reactionType == PLOG_RXN) {
510  rdata.rateCoeffType = PLOG_REACTION_RATECOEFF_TYPE;
511  for (size_t m = 0; m < kf.nChildren(); m++) {
512  const XML_Node& node = kf.child(m);
513  double p = getFloat(node, "P", "toSI");
514  vector_fp& rate = rdata.plogParameters.insert(
515  std::make_pair(p, vector_fp()))->second;
516  rate.resize(3);
517  rate[0] = getFloat(node, "A", "toSI");
518  rate[1] = getFloat(node, "b");
519  rate[2] = getFloat(node, "E", "actEnergy") / GasConstant;
520  }
521 
522  } else if (rdata.reactionType == CHEBYSHEV_RXN) {
523  rdata.rateCoeffType = CHEBYSHEV_REACTION_RATECOEFF_TYPE;
524  rdata.chebTmin = getFloat(kf, "Tmin", "toSI");
525  rdata.chebTmax = getFloat(kf, "Tmax", "toSI");
526  rdata.chebPmin = getFloat(kf, "Pmin", "toSI");
527  rdata.chebPmax = getFloat(kf, "Pmax", "toSI");
528  const XML_Node& coeffs = kf.child("floatArray");
529  rdata.chebDegreeP = atoi(coeffs["degreeP"].c_str());
530  rdata.chebDegreeT = atoi(coeffs["degreeT"].c_str());
531  getFloatArray(kf, rdata.chebCoeffs, false);
532 
533  } else {
534 
535  string type = kf.attrib("type");
536  if (type == "") {
537  type = "Arrhenius";
538  rdata.rateCoeffType = ARRHENIUS_REACTION_RATECOEFF_TYPE;
539  }
540  if (type == "ExchangeCurrentDensity") {
541  rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE;
542  } else if (type == "Arrhenius") {
543 
544  } else {
545  throw CanteraError("getRateCoefficient", "Unknown type: " + type);
546  }
547 
548  vector_fp clow(3,0.0), chigh(3,0.0);
549  for (size_t m = 0; m < kf.nChildren(); m++) {
550  const XML_Node& c = kf.child(m);
551  string nm = c.name();
552  int highlow=0;
553 
554  if (nm == "Arrhenius") {
555  vector_fp coeff(3);
556  if (c["type"] == "stick") {
557  getStick(c, kin, rdata, coeff[0], coeff[1], coeff[2]);
558  chigh = coeff;
559  } else {
560  getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]);
561  if (highlow == 1 || rdata.reactionType == THREE_BODY_RXN
562  || rdata.reactionType == ELEMENTARY_RXN) {
563  chigh = coeff;
564  } else {
565  clow = coeff;
566  }
567  }
568  if (rdata.reactionType == SURFACE_RXN) {
569  getCoverageDependence(c,
570  kin.thermo(kin.surfacePhaseIndex()), rdata);
571  }
572 
573  if (coeff[0] <= 0.0 && !rules.allowNegativeA) {
574  throw CanteraError("getRateCoefficient",
575  "negative or zero A coefficient for reaction "+int2str(rdata.number));
576  }
577  } else if (nm == "Arrhenius_ExchangeCurrentDensity") {
578  vector_fp coeff(3);
579  getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]);
580  chigh = coeff;
581  rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE;
582  } else if (nm == "falloff") {
583  getFalloff(c, rdata);
584  } else if (nm == "efficiencies") {
585  getEfficiencies(c, kin, rdata, rules);
586  } else if (nm == "electrochem") {
587  rdata.beta = fpValue(c["beta"]);
588  }
589  }
590  /*
591  * Store the coefficients in the ReactionData object for return
592  * from this function.
593  */
594  if (rdata.reactionType == CHEMACT_RXN) {
595  rdata.rateCoeffParameters = clow;
596  } else {
597  rdata.rateCoeffParameters = chigh;
598  }
599 
600  if (rdata.reactionType == FALLOFF_RXN) {
601  rdata.auxRateCoeffParameters = clow;
602  } else if (rdata.reactionType == CHEMACT_RXN) {
603  rdata.auxRateCoeffParameters = chigh;
604  }
605  }
606 }
607 
608 /*
609  * This function returns true if two reactions are duplicates of
610  * one another, and false otherwise. The input arguments are two
611  * maps from species number to stoichiometric coefficient, one for
612  * each reaction. The reactions are considered duplicates if their
613  * stoichiometric coefficients have the same ratio for all
614  * species.
615  */
616 doublereal isDuplicateReaction(std::map<int, doublereal>& r1,
617  std::map<int, doublereal>& r2)
618 {
619 
620  map<int, doublereal>::const_iterator b = r1.begin(), e = r1.end();
621  int k1 = b->first;
622  doublereal ratio = 0.0;
623  if (r1[k1] == 0.0 || r2[k1] == 0.0) {
624  goto next;
625  }
626  ratio = r2[k1]/r1[k1];
627  ++b;
628  for (; b != e; ++b) {
629  k1 = b->first;
630  if (r1[k1] == 0.0 || r2[k1] == 0.0) {
631  goto next;
632  }
633  if (fabs(r2[k1]/r1[k1] - ratio) > 1.e-8) {
634  goto next;
635  }
636  }
637  return ratio;
638 next:
639  ratio = 0.0;
640  b = r1.begin();
641  k1 = b->first;
642  if (r1[k1] == 0.0 || r2[-k1] == 0.0) {
643  return 0.0;
644  }
645  ratio = r2[-k1]/r1[k1];
646  ++b;
647  for (; b != e; ++b) {
648  k1 = b->first;
649  if (r1[k1] == 0.0 || r2[-k1] == 0.0) {
650  return 0.0;
651  }
652  if (fabs(r2[-k1]/r1[k1] - ratio) > 1.e-8) {
653  return 0.0;
654  }
655  }
656  return ratio;
657 }
658 
659 
660 /**
661  * Install an individual reaction into a kinetics manager. The
662  * data for the reaction is in the xml_node r. In other words, r
663  * points directly to a ctml element named "reaction". i refers
664  * to the number id of the reaction in the kinetics object.
665  *
666  * @param iRxn Reaction number.
667  * @param r XML_Node containing reaction data.
668  * @param kin Kinetics manager to which reaction will be added.
669  * @param default_phase Default phase for locating a species
670  * @param rules Rule for handling reactions with missing species
671  * (skip or flag as error)
672  * @param validate_rxn If true, check that this reaction is not a
673  * duplicate of one already entered, and check that the reaction
674  * balances.
675  *
676  * @ingroup kineticsmgr
677  */
678 bool rxninfo::installReaction(int iRxn, const XML_Node& r, Kinetics& kin,
679  string default_phase, ReactionRules& rules,
680  bool validate_rxn)
681 {
682  // Check to see that we are in fact at a reaction node
683  if (r.name() != "reaction") {
684  throw CanteraError(" rxninfo::installReaction",
685  " expected xml node reaction, got " + r.name());
686  }
687 
688  // We use the ReactionData object to store initial values read in from the
689  // xml data. Then, when we have collected everything we add the reaction to
690  // the kinetics object, kin, at the end of the routine.
691  ReactionData rdata;
692  rdata.validate = validate_rxn;
693 
694  // Check to see if the reaction is specified to be a duplicate of another
695  // reaction. It's an error if the reaction is a duplicate and this is not
696  // set.
697  int dup = (r.hasAttrib("duplicate")) ? 1 : 0;
698 
699  // Check to see if the reaction rate constant can be negative. It's an
700  // error if a negative rate constant is found and this is not set.
701  rules.allowNegativeA = (r.hasAttrib("negative_A")) ? 1 : 0;
702 
703  // Use the contents of the "equation" child element as the reaction's
704  // string representation. Post-process to convert "[" and "]" characters
705  // back into "<" and ">" which cannot easily be stored in an XML file. This
706  // reaction string is used only for display purposes. It is not parsed for
707  // the identities of reactants or products.
708  string eqn = (r.hasChild("equation")) ? r("equation") : "<no equation>";
709  for (size_t nn = 0; nn < eqn.size(); nn++) {
710  if (eqn[nn] == '[') {
711  eqn[nn] = '<';
712  } else if (eqn[nn] == ']') {
713  eqn[nn] = '>';
714  }
715  }
716 
717  // get the reactants
718  bool ok = getReagents(r, kin, 1, default_phase, rdata.reactants,
719  rdata.rstoich, rdata.rorder, rules);
720 
721  // Get the products. We store the id of products in rdata.products
722  ok = ok && getReagents(r, kin, -1, default_phase, rdata.products,
723  rdata.pstoich, rdata.porder, rules);
724 
725  // if there was a problem getting either the reactants or the products,
726  // then abort.
727  if (!ok) {
728  return false;
729  }
730 
731  // check whether the reaction is specified to be
732  // reversible. Default is irreversible.
733  string isrev = r["reversible"];
734  rdata.reversible = (isrev == "yes" || isrev == "true");
735 
736  // If reaction orders are specified, then this reaction does not follow
737  // mass-action kinetics, and is not an elementary reaction. So check that
738  // it is not reversible, since computing the reverse rate from
739  // thermochemistry only works for elementary reactions. Set the type to
740  // global, so that kinetics managers will know to process the reaction
741  // orders.
742  if (r.hasChild("order")) {
743  if (rdata.reversible == true)
744  throw CanteraError("installReaction",
745  "reaction orders may only be given for "
746  "irreversible reactions");
747  rdata.global = true;
748  }
749 
750  // Some reactions can be elementary reactions but have fractional
751  // stoichiometries wrt to some products and reactants. An example of these
752  // are solid reactions involving phase transformations. Species with
753  // fractional stoichiometries must be from single-species phases with
754  // unity activities. For these reactions set the bool isReversibleWithFrac
755  // to true.
756  if (rdata.reversible == true) {
757  for (size_t i = 0; i < rdata.products.size(); i++) {
758  doublereal po = rdata.porder[i];
759  AssertTrace(po == rdata.pstoich[i]);
760  doublereal chk = po - 1.0 * int(po);
761  if (chk != 0.0) {
762  size_t k = rdata.products[i];
763  // Special case when k is a single species phase.
764  if (kin.speciesPhase(k).nSpecies() == 1) {
765  rdata.porder[i] = 0.0;
766  }
767 
768  rdata.isReversibleWithFrac = true;
769  }
770  }
771  for (size_t i = 0; i < rdata.reactants.size(); i++) {
772  doublereal ro = rdata.rorder[i];
773  AssertTrace(ro == rdata.rstoich[i]);
774  doublereal chk = ro - 1.0 * int(ro);
775  if (chk != 0.0) {
776  size_t k = rdata.reactants[i];
777  // Special case when k is a single species phase.
778  if (kin.speciesPhase(k).nSpecies() == 1) {
779  rdata.rorder[i] = 0.0;
780  }
781  rdata.isReversibleWithFrac = true;
782  }
783  }
784  }
785 
786  /*
787  * Search the reaction element for the attribute "type".
788  * If found, then branch on the type, to fill in appropriate
789  * fields in rdata.
790  */
791  rdata.reactionType = ELEMENTARY_RXN;
792  string typ = r["type"];
793  if (typ == "falloff") {
794  rdata.reactionType = FALLOFF_RXN;
795  rdata.falloffType = SIMPLE_FALLOFF;
796  } else if (typ == "chemAct") {
797  rdata.reactionType = CHEMACT_RXN;
798  rdata.falloffType = SIMPLE_FALLOFF;
799  } else if (typ == "threeBody") {
800  rdata.reactionType = THREE_BODY_RXN;
801  } else if (typ == "plog") {
802  rdata.reactionType = PLOG_RXN;
803  } else if (typ == "chebyshev") {
804  rdata.reactionType = CHEBYSHEV_RXN;
805  } else if (typ == "surface") {
806  rdata.reactionType = SURFACE_RXN;
807  } else if (typ == "edge") {
808  rdata.reactionType = EDGE_RXN;
809  } else if (typ != "") {
810  throw CanteraError("installReaction", "Unknown reaction type: " + typ);
811  }
812 
813  // Look for undeclared duplicate reactions.
814  if (validate_rxn) {
815  map<int, doublereal> rxnstoich;
816  vector<char> participants(kin.nTotalSpecies(), 0);
817  for (size_t nn = 0; nn < rdata.reactants.size(); nn++) {
818  rxnstoich[-1 - int(rdata.reactants[nn])] -= rdata.rstoich[nn];
819  participants[rdata.reactants[nn]] += 1;
820  }
821  for (size_t nn = 0; nn < rdata.products.size(); nn++) {
822  rxnstoich[int(rdata.products[nn])+1] += rdata.pstoich[nn];
823  participants[rdata.products[nn]] += 2;
824  }
825 
826  vector<size_t>& related = m_participants[participants];
827  for (size_t mm = 0; mm < related.size(); mm++) {
828  size_t nn = related[mm];
829  if ((rdata.reactants.size() == m_nr[nn])
830  && (rdata.reactionType == m_typ[nn])) {
831  doublereal c = isDuplicateReaction(rxnstoich, m_rdata[nn]);
832  if (c > 0.0
833  || (c < 0.0 && rdata.reversible)
834  || (c < 0.0 && m_rev[nn])) {
835  if ((!dup || !m_dup[nn])) {
836  string msg = string("Undeclared duplicate reactions detected: \n")
837  +"Reaction "+int2str(nn+1)+": "+m_eqn[nn]
838  +"\nReaction "+int2str(iRxn+1)+": "+eqn+"\n";
839  throw CanteraError("installReaction", msg);
840  }
841  }
842  }
843  }
844  m_dup.push_back(dup);
845  m_rev.push_back(rdata.reversible);
846  m_eqn.push_back(eqn);
847  m_nr.push_back(rdata.reactants.size());
848  m_typ.push_back(rdata.reactionType);
849  m_rdata.push_back(rxnstoich);
850  m_participants[participants].push_back(m_rdata.size() - 1);
851  }
852 
853  rdata.equation = eqn;
854  rdata.number = iRxn;
855  rdata.rxn_number = iRxn;
856 
857  // Read the rate coefficient data from the XML file. Trigger an
858  // exception for negative A unless specifically authorized.
859  getRateCoefficient(r.child("rateCoeff"), kin, rdata, rules);
860 
861  // Check to see that the elements balance in the reaction.
862  // Throw an error if they don't
863  if (validate_rxn) {
864  checkRxnElementBalance(kin, rdata);
865  }
866 
867  // Ok we have read everything in about the reaction. Add it to the
868  // kinetics object by calling the Kinetics member function addReaction()
869  kin.addReaction(rdata);
870  return true;
871 }
872 
873 /*
874  * Take information from the XML tree, p, about reactions
875  * and install them into the kinetics object, kin.
876  * default_phase is the default phase to assume when
877  * looking up species.
878  *
879  * At this point, p usually refers to the phase xml element.
880  * One of the children of this element is reactionArray,
881  * the element which determines where in the xml file to
882  * look up the reaction rate data pertaining to the phase.
883  *
884  * On return, if reaction instantiation goes correctly, return true.
885  * If there is a problem, return false.
886  */
888  std::string default_phase, bool check_for_duplicates)
889 {
890  const std::auto_ptr<rxninfo> _rxns(new rxninfo);
891 
892  vector<XML_Node*> rarrays;
893  int itot = 0;
894  /*
895  * Search the children of the phase element for the
896  * xml element named reactionArray. If we can't find it,
897  * then return signaling having not found any reactions.
898  * Apparently, we allow multiple reactionArray elements here
899  * Each one will be processed sequentially, with the
900  * end result being purely additive.
901  */
902  p.getChildren("reactionArray",rarrays);
903  int na = static_cast<int>(rarrays.size());
904  if (na == 0) {
905  kin.finalize();
906  return false;
907  }
908  for (int n = 0; n < na; n++) {
909  /*
910  * Go get a reference to the current xml element,
911  * reactionArray. We will process this element now.
912  */
913  const XML_Node& rxns = *rarrays[n];
914  /*
915  * The reactionArray element has an attribute called,
916  * datasrc. The value of the attribute is the xml
917  * element comprising the top of the
918  * tree of reactions for the phase.
919  * Find this datasrc element starting with the root
920  * of the current xml node.
921  */
922  const XML_Node* rdata = get_XML_Node(rxns["datasrc"], &rxns.root());
923  /*
924  * If the reactionArray element has a child element named "skip", and
925  * if the attribute of skip called "species" has a value of "undeclared",
926  * we will set rxnrule.skipUndeclaredSpecies to 'true'. rxnrule is
927  * passed to the routine that parses each individual reaction so that
928  * the parser will skip all reactions containing an undefined species
929  * without throwing an error.
930  *
931  * Similarly, an attribute named "third_bodies" with the value of
932  * "undeclared" will skip undeclared third body efficiencies (while
933  * retaining the reaction and any other efficiencies).
934  */
935  ReactionRules rxnrule;
936  if (rxns.hasChild("skip")) {
937  const XML_Node& sk = rxns.child("skip");
938  string sskip = sk["species"];
939  if (sskip == "undeclared") {
940  rxnrule.skipUndeclaredSpecies = true;
941  }
942  if (sk["third_bodies"] == "undeclared") {
943  rxnrule.skipUndeclaredThirdBodies = true;
944  }
945  }
946  int i, nrxns = 0;
947  /*
948  * Search for child elements called include. We only include
949  * a reaction if it's tagged by one of the include fields.
950  * Or, we include all reactions if there are no include fields.
951  */
952  vector<XML_Node*> incl;
953  rxns.getChildren("include",incl);
954  int ninc = static_cast<int>(incl.size());
955 
956  vector<XML_Node*> allrxns;
957  rdata->getChildren("reaction",allrxns);
958  nrxns = static_cast<int>(allrxns.size());
959  // if no 'include' directive, then include all reactions
960  if (ninc == 0) {
961  for (i = 0; i < nrxns; i++) {
962  const XML_Node* r = allrxns[i];
963  if (r) {
964  if (_rxns->installReaction(itot, *r, kin,
965  default_phase, rxnrule, check_for_duplicates)) {
966  ++itot;
967  }
968  }
969  }
970  } else {
971  for (int nii = 0; nii < ninc; nii++) {
972  const XML_Node& ii = *incl[nii];
973  string imin = ii["min"];
974  string imax = ii["max"];
975 
976  string::size_type iwild = string::npos;
977  if (imax == imin) {
978  iwild = imin.find("*");
979  if (iwild != string::npos) {
980  imin = imin.substr(0,iwild);
981  imax = imin;
982  }
983  }
984 
985  for (i = 0; i < nrxns; i++) {
986  const XML_Node* r = allrxns[i];
987  string rxid;
988  if (r) {
989  rxid = (*r)["id"];
990  if (iwild != string::npos) {
991  rxid = rxid.substr(0,iwild);
992  }
993  /*
994  * To decide whether the reaction is included or not
995  * we do a lexical min max and operation. This
996  * sometimes has surprising results.
997  */
998  if ((rxid >= imin) && (rxid <= imax)) {
999  if (_rxns->installReaction(itot, *r, kin,
1000  default_phase, rxnrule, check_for_duplicates)) {
1001  ++itot;
1002  }
1003  }
1004  }
1005  }
1006  }
1007  }
1008  }
1009 
1010  /*
1011  * Finalize the installation of the kinetics, now that we know
1012  * the true number of reactions in the mechanism, itot.
1013  */
1014  kin.finalize();
1015 
1016  return true;
1017 }
1018 
1019 
1020 /*
1021  * Import a reaction mechanism for a phase or an interface.
1022  *
1023  * @param phase This is an xml node containing a description
1024  * of a phase. Within the phase is a XML element
1025  * called reactionArray containing the location
1026  * of the description of the reactions that make
1027  * up the kinetics object.
1028  * Also within the phase is an XML element called
1029  * phaseArray containing a listing of other phases
1030  * that participate in the kinetics mechanism.
1031  *
1032  * @param th This is a list of ThermoPhase pointers which must
1033  * include all of
1034  * the phases that participate in the kinetics
1035  * operator. All of the phases must have already
1036  * been initialized and formed within Cantera.
1037  * However, their pointers should not have been
1038  * added to the Kinetics object; this addition
1039  * is carried out here. Additional phases may
1040  * be include; these have no effect.
1041  *
1042  * @param k This is a pointer to the kinetics manager class
1043  * that will be initialized with a kinetics
1044  * mechanism.
1045  */
1046 bool importKinetics(const XML_Node& phase, std::vector<ThermoPhase*> th,
1047  Kinetics* k)
1048 {
1049 
1050  if (k == 0) {
1051  return false;
1052  }
1053 
1054  Kinetics& kin = *k;
1055 
1056  // This phase will be the owning phase for the kinetics operator
1057  // For interfaces, it is the surface phase between two volumes.
1058  // For homogeneous kinetics, it's the current volumetric phase.
1059  string owning_phase = phase["id"];
1060 
1061  bool check_for_duplicates = false;
1062  if (phase.parent()->hasChild("validate")) {
1063  const XML_Node& d = phase.parent()->child("validate");
1064  if (d["reactions"] == "yes") {
1065  check_for_duplicates = true;
1066  }
1067  }
1068 
1069  // if other phases are involved in the reaction mechanism,
1070  // they must be listed in a 'phaseArray' child
1071  // element. Homogeneous mechanisms do not need to include a
1072  // phaseArray element.
1073 
1074  vector<string> phase_ids;
1075  if (phase.hasChild("phaseArray")) {
1076  const XML_Node& pa = phase.child("phaseArray");
1077  getStringArray(pa, phase_ids);
1078  }
1079  phase_ids.push_back(owning_phase);
1080 
1081  int np = static_cast<int>(phase_ids.size());
1082  int nt = static_cast<int>(th.size());
1083 
1084  // for each referenced phase, attempt to find its id among those
1085  // phases specified.
1086  bool phase_ok;
1087 
1088  string phase_id;
1089  string msg = "";
1090  for (int n = 0; n < np; n++) {
1091  phase_id = phase_ids[n];
1092  phase_ok = false;
1093 
1094  // loop over the supplied 'ThermoPhase' objects representing
1095  // phases, to find an object with the same id.
1096  for (int m = 0; m < nt; m++) {
1097  if (th[m]->id() == phase_id) {
1098  phase_ok = true;
1099 
1100  // if no phase with this id has been added to
1101  //the kinetics manager yet, then add this one
1102  if (kin.phaseIndex(phase_id) == npos) {
1103  kin.addPhase(*th[m]);
1104  }
1105  }
1106  msg += " "+th[m]->id();
1107  }
1108  if (!phase_ok) {
1109  throw CanteraError("importKinetics",
1110  "phase "+phase_id+" not found. Supplied phases are:"+msg);
1111  }
1112  }
1113 
1114  // allocates arrays, etc. Must be called after the phases have
1115  // been added to 'kin', so that the number of species in each
1116  // phase is known.
1117  kin.init();
1118 
1119  // Install the reactions.
1120  return installReactionArrays(phase, kin, owning_phase, check_for_duplicates);
1121 }
1122 
1123 /*
1124  * Build a single-phase ThermoPhase object with associated kinetics
1125  * mechanism.
1126  */
1127 bool buildSolutionFromXML(XML_Node& root, std::string id, std::string nm,
1128  ThermoPhase* th, Kinetics* k)
1129 {
1130  XML_Node* x;
1131  x = get_XML_NameID(nm, string("#")+id, &root);
1132  // x = get_XML_Node(string("#")+id, &root);
1133  if (!x) {
1134  return false;
1135  }
1136 
1137  /*
1138  * Fill in the ThermoPhase object by querying the
1139  * const XML_Node tree located at x.
1140  */
1141  importPhase(*x, th);
1142  /*
1143  * Create a vector of ThermoPhase pointers of length 1
1144  * having the current th ThermoPhase as the entry.
1145  */
1146  vector<ThermoPhase*> phases(1);
1147  phases[0] = th;
1148  /*
1149  * Fill in the kinetics object k, by querying the
1150  * const XML_Node tree located by x. The source terms and
1151  * eventually the source term vector will be constructed
1152  * from the list of ThermoPhases in the vector, phases.
1153  */
1154  importKinetics(*x, phases, k);
1155  return true;
1156 }
1157 
1158 }