Cantera  2.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
19 #include "cantera/base/ctml.h"
20 
21 #include <cstring>
22 
23 using namespace std;
24 
25 namespace Cantera
26 {
27 
28 ReactionRules::ReactionRules() :
29  skipUndeclaredSpecies(false),
30  skipUndeclaredThirdBodies(false),
31  allowNegativeA(false)
32 {
33 }
34 
36  const ReactionData& rdata, doublereal errorTolerance)
37 {
38  warn_deprecated("checkRxnElementBalance", "Now handled by "
39  "Kinetics::checkReactionBalance. To be removed after Cantera 2.2.");
40  doublereal kstoich;
41 
42  map<string, double> bal, balr, balp;
43  bal.clear();
44  balp.clear();
45  balr.clear();
46  size_t np = rdata.products.size();
47 
48  // iterate over the products
49  for (size_t index = 0; index < np; index++) {
50  size_t kp = rdata.products[index]; // index of the product in 'kin'
51  size_t n = kin.speciesPhaseIndex(kp); // phase this product belongs to
52  size_t klocal = kp - kin.kineticsSpeciesIndex(0,n); // index within this phase
53  kstoich = rdata.pstoich[index]; // product stoichiometric coeff
54  const ThermoPhase& ph = kin.speciesPhase(kp);
55  for (size_t m = 0; m < ph.nElements(); m++) {
56  bal[ph.elementName(m)] += kstoich*ph.nAtoms(klocal,m);
57  balp[ph.elementName(m)] += kstoich*ph.nAtoms(klocal,m);
58  }
59  }
60  for (size_t index = 0; index < rdata.reactants.size(); index++) {
61  size_t kr = rdata.reactants[index];
62  size_t n = kin.speciesPhaseIndex(kr);
63  size_t klocal = kr - kin.kineticsSpeciesIndex(0,n);
64  kstoich = rdata.rstoich[index];
65  const ThermoPhase& ph = kin.speciesPhase(kr);
66  for (size_t m = 0; m < ph.nElements(); m++) {
67  bal[ph.elementName(m)] -= kstoich*ph.nAtoms(klocal,m);
68  balr[ph.elementName(m)] += kstoich*ph.nAtoms(klocal,m);
69  }
70  }
71 
72  map<string, double>::iterator b = bal.begin();
73  string msg = "\n\tElement Reactants Products";
74  bool ok = true;
75  doublereal err, elemsum;
76  for (; b != bal.end(); ++b) {
77  elemsum = fabs(balr[b->first]) + fabs(balp[b->first]);
78  if (elemsum > 0.0) {
79  err = fabs(b->second/elemsum);
80  if (err > errorTolerance) {
81  ok = false;
82  msg += "\n\t"+b->first+" "+ fp2str(balr[b->first])
83  +" "+ fp2str(balp[b->first]);
84  }
85  }
86  }
87  if (!ok) {
88  msg = "The following reaction is unbalanced:\n\t"
89  + rdata.equation + "\n" + msg + "\n";
90  throw CanteraError("checkRxnElementBalance",msg);
91  }
92 }
93 
94 bool getReagents(const XML_Node& rxn, Kinetics& kin, int rp,
95  std::string default_phase, std::vector<size_t>& spnum,
96  vector_fp& stoich, vector_fp& order,
97  const ReactionRules& rules)
98 {
99  warn_deprecated("getReagents", "Now handled through newReaction() and its "
100  "support functions. To be removed after Cantera 2.2.");
101  string rptype;
102 
103  /*
104  * The id of reactants and products are kept in child elements
105  * of reaction, named "reactants" and "products". We search
106  * the XML tree for these children based on the value of rp,
107  * and store the XML element pointer here.
108  */
109  if (rp == 1) {
110  rptype = "reactants";
111  } else {
112  rptype = "products";
113  }
114  const XML_Node& rg = rxn.child(rptype);
115 
116  /*
117  * The species and stoichiometric coefficient for the species
118  * are stored as a colon separated pair. Get all of these
119  * pairs in the reactions/products object.
120  */
121  std::vector<string> key, val;
122  getPairs(rg, key, val);
123 
124  /*
125  * Loop over each of the pairs and process them
126  */
127  doublereal ord, stch;
128  string ph, spName;
129  map<string, size_t> speciesMap;
130  for (size_t n = 0; n < key.size(); n++) {
131  spName = key[n]; // sp is the string name for species
132  ph = "";
133  /*
134  * Search for the species in the kinetics object using the
135  * member function kineticsSpeciesIndex(). We will search
136  * for the species in all phases defined in the kinetics operator.
137  */
138  size_t isp = kin.kineticsSpeciesIndex(spName);
139  if (isp == npos) {
140  if (rules.skipUndeclaredSpecies) {
141  return false;
142  } else {
143  throw CanteraError("getReagents",
144  "Undeclared reactant or product species " + spName);
145  return false;
146  }
147  }
148 
149  /*
150  * For each reagent, we store the the species number, isp
151  * the stoichiometric coefficient, val[n], and the order
152  * species in the reaction rate expression. We assume mass
153  * action kinetics here, but will modify this below for
154  * specified species.
155  */
156  spnum.push_back(isp);
157  stch = fpValue(val[n]);
158  stoich.push_back(stch);
159  ord = doublereal(stch);
160  order.push_back(ord);
161 
162  /*
163  * Needed to process reaction orders below.
164  */
165  speciesMap[spName] = order.size();
166  }
167 
168  /*
169  * Check to see if reaction orders have been specified.
170  */
171 
172  if (rp == 1 && rxn.hasChild("order")) {
173  std::vector<XML_Node*> ord = rxn.getChildren("order");
174  doublereal forder;
175  for (size_t nn = 0; nn < ord.size(); nn++) {
176  const XML_Node& oo = *ord[nn];
177  string sp = oo["species"];
178  size_t loc = speciesMap[sp];
179  if (loc == 0)
180  throw CanteraError("getReagents",
181  "reaction order specified for non-reactant: "
182  +sp);
183  forder = oo.fp_value();
184  if (forder < 0.0) {
185  throw CanteraError("getReagents",
186  "reaction order must be non-negative");
187  }
188  // replace the stoichiometric coefficient
189  // stored above in 'order' with the specified
190  // reaction order
191  order[loc-1] = forder;
192  }
193  }
194  return true;
195 }
196 
197 /**
198  * getArrhenius() parses the XML element called Arrhenius.
199  * The Arrhenius expression is
200  * \f[ k = A T^(b) exp (-E_a / RT). \f]
201  * @deprecated to be removed after Cantera 2.2.
202  */
203 static void getArrhenius(const XML_Node& node, int& labeled,
204  doublereal& A, doublereal& b, doublereal& E)
205 {
206  if (node["name"] == "k0") {
207  labeled = -1;
208  } else if (node["name"] == "kHigh") {
209  labeled = 1;
210  } else {
211  labeled = 0;
212  }
213  /*
214  * We parse the children for the A, b, and E components.
215  */
216  A = getFloat(node, "A", "toSI");
217  b = getFloat(node, "b");
218  E = getFloat(node, "E", "actEnergy");
219  E /= GasConstant;
220 }
221 
222 /**
223  * getStick() processes the XML element called Stick that specifies
224  * the sticking coefficient reaction. This routine will
225  * translate the sticking coefficient value into a "normal"
226  * rate constant for the surface reaction.
227  *
228  * Output
229  * -----------
230  * Output is the normal Arrhenius expressions for a surface
231  * reaction rate constant.
232  *
233  * A - units such that rate of rxn has kmol/m^2/s when
234  * A is multiplied by activity concentrations of
235  * reactants in the normal manner.
236  * n - unitless
237  * E - Units 1/Kelvin
238  * @deprecated to be removed after Cantera 2.2.
239  */
240 static void getStick(const XML_Node& node, Kinetics& kin,
241  ReactionData& r, doublereal& A, doublereal& b, doublereal& E)
242 {
243  size_t nr = r.reactants.size();
244  size_t k, klocal, not_surf = 0;
245  size_t np = 0;
246  doublereal f = 1.0;
247  doublereal order;
248  /*
249  * species is the name of the special reactant whose surface
250  * flux rate will be calculated.
251  * isp = species # in the local phase
252  * ispKinetics = species # in the kinetics object
253  * ispPhaseIndex = phase # of the special species
254  */
255  string spname = node["species"];
256  ThermoPhase& th = kin.speciesPhase(spname);
257  size_t isp = th.speciesIndex(spname);
258  size_t ispKinetics = kin.kineticsSpeciesIndex(spname);
259  size_t ispPhaseIndex = kin.speciesPhaseIndex(ispKinetics);
260 
261  doublereal ispMW = th.molecularWeights()[isp];
262  doublereal sc;
263 
264  // loop over the reactants
265  for (size_t n = 0; n < nr; n++) {
266  k = r.reactants[n];
267  order = r.rorder[n]; // stoich coeff
268 
269  // get the phase species k belongs to
270  np = kin.speciesPhaseIndex(k);
271  const ThermoPhase& p = kin.thermo(np);
272 
273  // get the local index of species k in this phase
274  klocal = p.speciesIndex(kin.kineticsSpeciesName(k));
275 
276  // if it is a surface species, divide f by the standard
277  // concentration for this species, in order to convert
278  // from concentration units used in the law of mass action
279  // to coverages used in the sticking probability expression
280  if (p.eosType() == cSurf || p.eosType() == cEdge) {
281  sc = p.standardConcentration(klocal);
282  f /= pow(sc, order);
283  }
284  // Otherwise:
285  else {
286  // We only allow one species to be in the phase containing the
287  // special sticking coefficient species.
288  if (ispPhaseIndex == np) {
289  not_surf++;
290  }
291  // Other bulk phase species on the other side of ther interface are
292  // treated like surface species.
293  else {
294  sc = p.standardConcentration(klocal);
295  f /= pow(sc, order);
296  }
297  }
298  }
299  if (not_surf != 1) {
300  throw CanteraError("getStick",
301  "reaction probabilities can only be used in "
302  "reactions with exactly 1 gas/liquid species.");
303  }
304 
305  doublereal cbar = sqrt(8.0*GasConstant/(Pi*ispMW));
306  A = 0.25 * getFloat(node, "A", "toSI") * cbar * f;
307  b = getFloat(node, "b") + 0.5;
308  E = getFloat(node, "E", "actEnergy");
309  E /= GasConstant;
310 }
311 
312 //! Read the XML data concerning the coverage dependence of an interfacial reaction
313 /*!
314  * @param node XML node with name reaction containing the reaction information
315  * @param surfphase Surface phase
316  * @param rdata Reaction data for the reaction.
317  *
318  * Example:
319  * @verbatim
320  <coverage species="CH3*">
321  <a> 1.0E-5 </a>
322  <m> 0.0 </m>
323  <actEnergy> 0.0 </actEnergy>
324  </coverage>
325 @endverbatim
326  * @deprecated to be removed after Cantera 2.2.
327  */
328 static void getCoverageDependence(const XML_Node& node,
329  thermo_t& surfphase, ReactionData& rdata)
330 {
331  vector<XML_Node*> cov = node.getChildren("coverage");
332  size_t k, nc = cov.size();
333  doublereal e;
334  string spname;
335  if (nc > 0) {
336  for (size_t n = 0; n < nc; n++) {
337  const XML_Node& cnode = *cov[n];
338  spname = cnode["species"];
339  k = surfphase.speciesIndex(spname);
340  rdata.cov.push_back(doublereal(k));
341  rdata.cov.push_back(getFloat(cnode, "a"));
342  rdata.cov.push_back(getFloat(cnode, "m"));
343  e = getFloat(cnode, "e", "actEnergy");
344  rdata.cov.push_back(e/GasConstant);
345  }
346  }
347 }
348 
349 //! Get falloff parameters for a reaction.
350 /*!
351  * This routine reads the falloff XML node and extracts parameters into a
352  * vector of doubles
353  *
354  * @verbatim
355  <falloff type="Troe"> 0.5 73.2 5000. 9999. </falloff>
356  @endverbatim
357  * @deprecated to be removed after Cantera 2.2.
358 */
359 static void getFalloff(const XML_Node& f, ReactionData& rdata)
360 {
361  string type = f["type"];
362  vector<string> p;
363  getStringArray(f,p);
364  vector_fp c;
365  size_t np = p.size();
366  for (size_t n = 0; n < np; n++) {
367  c.push_back(fpValue(p[n]));
368  }
369  if (type == "Troe") {
370  if (np == 3 || np == 4) {
371  rdata.falloffType = TROE_FALLOFF;
372  } else {
373  throw CanteraError("getFalloff()", "Troe parameterization is specified by number of parameters, "
374  + int2str(np) + ", is not equal to 3 or 4");
375  }
376  } else if (type == "SRI") {
377  if (np == 3 || np == 5) {
378  rdata.falloffType = SRI_FALLOFF;
379  } else {
380  throw CanteraError("getFalloff()", "SRI parameterization is specified by number of parameters, "
381  + int2str(np) + ", is not equal to 3 or 5");
382  }
383  }
384  rdata.falloffParameters = c;
385 }
386 
387 /**
388  * Get the enhanced collision efficiencies. It is assumed that the
389  * reaction mechanism is homogeneous, so that all species belong
390  * to phase(0) of 'kin'.
391  * @deprecated to be removed after Cantera 2.2.
392  */
393 static void getEfficiencies(const XML_Node& eff, Kinetics& kin,
394  ReactionData& rdata, const ReactionRules& rules)
395 {
396  // set the default collision efficiency
397  rdata.default_3b_eff = fpValue(eff["default"]);
398 
399  vector<string> key, val;
400  getPairs(eff, key, val);
401  string nm;
402  string phse = kin.thermo(0).id();
403  for (size_t n = 0; n < key.size(); n++) {
404  nm = key[n];
405  size_t k = kin.kineticsSpeciesIndex(nm, phse);
406  if (k != npos) {
407  rdata.thirdBodyEfficiencies[k] = fpValue(val[n]);
408  } else if (!rules.skipUndeclaredThirdBodies) {
409  throw CanteraError("getEfficiencies", "Encountered third-body "
410  "efficiency for undefined species \"" + nm + "\"\n"
411  "while adding reaction " + int2str(rdata.number+1) + ".");
412  }
413  }
414 }
415 
416 void getRateCoefficient(const XML_Node& kf, Kinetics& kin,
417  ReactionData& rdata, const ReactionRules& rules)
418 {
419  warn_deprecated("getRateCoefficent", "Now handled through newReaction() "
420  "and its support functions. To be removed after Cantera 2.2.");
421  if (rdata.reactionType == PLOG_RXN) {
422  rdata.rateCoeffType = PLOG_REACTION_RATECOEFF_TYPE;
423  for (size_t m = 0; m < kf.nChildren(); m++) {
424  const XML_Node& node = kf.child(m);
425  double p = getFloat(node, "P", "toSI");
426  vector_fp& rate = rdata.plogParameters.insert(
427  std::make_pair(p, vector_fp()))->second;
428  rate.resize(3);
429  rate[0] = getFloat(node, "A", "toSI");
430  rate[1] = getFloat(node, "b");
431  rate[2] = getFloat(node, "E", "actEnergy") / GasConstant;
432  }
433 
434  } else if (rdata.reactionType == CHEBYSHEV_RXN) {
435  rdata.rateCoeffType = CHEBYSHEV_REACTION_RATECOEFF_TYPE;
436  rdata.chebTmin = getFloat(kf, "Tmin", "toSI");
437  rdata.chebTmax = getFloat(kf, "Tmax", "toSI");
438  rdata.chebPmin = getFloat(kf, "Pmin", "toSI");
439  rdata.chebPmax = getFloat(kf, "Pmax", "toSI");
440  const XML_Node& coeffs = kf.child("floatArray");
441  rdata.chebDegreeP = atoi(coeffs["degreeP"].c_str());
442  rdata.chebDegreeT = atoi(coeffs["degreeT"].c_str());
443  getFloatArray(kf, rdata.chebCoeffs, false);
444 
445  } else {
446 
447  string type = kf.attrib("type");
448  if (type == "") {
449  type = "Arrhenius";
450  rdata.rateCoeffType = ARRHENIUS_REACTION_RATECOEFF_TYPE;
451  }
452  if (type == "ExchangeCurrentDensity") {
453  rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE;
454  } else if (type == "Arrhenius") {
455 
456  } else {
457  throw CanteraError("getRateCoefficient", "Unknown type: " + type);
458  }
459 
460  vector_fp c_alt(3,0.0), c_base(3,0.0);
461  for (size_t m = 0; m < kf.nChildren(); m++) {
462  const XML_Node& c = kf.child(m);
463  string nm = c.name();
464  int labeled=0;
465 
466  if (nm == "Arrhenius") {
467  vector_fp coeff(3);
468  if (c["type"] == "stick") {
469  getStick(c, kin, rdata, coeff[0], coeff[1], coeff[2]);
470  c_base = coeff;
471  } else {
472  getArrhenius(c, labeled, coeff[0], coeff[1], coeff[2]);
473  if (labeled == 0 || rdata.reactionType == THREE_BODY_RXN
474  || rdata.reactionType == ELEMENTARY_RXN) {
475  c_base = coeff;
476  } else {
477  c_alt = coeff;
478  }
479  }
480  if (rdata.reactionType == SURFACE_RXN || rdata.reactionType == EDGE_RXN) {
482  kin.thermo(kin.surfacePhaseIndex()), rdata);
483  }
484 
485  if (coeff[0] < 0.0 && !rules.allowNegativeA) {
486  throw CanteraError("getRateCoefficient",
487  "negative A coefficient for reaction "+int2str(rdata.number));
488  }
489  } else if (nm == "Arrhenius_ExchangeCurrentDensity") {
490  vector_fp coeff(3);
491  getArrhenius(c, labeled, coeff[0], coeff[1], coeff[2]);
492  c_base = coeff;
493  rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE;
494  } else if (nm == "falloff") {
495  getFalloff(c, rdata);
496  } else if (nm == "efficiencies") {
497  getEfficiencies(c, kin, rdata, rules);
498  } else if (nm == "electrochem") {
499  rdata.beta = fpValue(c["beta"]);
500  }
501  }
502  /*
503  * Store the coefficients in the ReactionData object for return
504  * from this function.
505  */
506  if (rdata.reactionType == FALLOFF_RXN) {
507  rdata.rateCoeffParameters = c_base;
508  rdata.auxRateCoeffParameters = c_alt;
509  } else if (rdata.reactionType == CHEMACT_RXN) {
510  rdata.rateCoeffParameters = c_alt;
511  rdata.auxRateCoeffParameters = c_base;
512  } else {
513  rdata.rateCoeffParameters = c_base;
514  }
515 
516  }
517 }
518 
519 doublereal isDuplicateReaction(std::map<int, doublereal>& r1,
520  std::map<int, doublereal>& r2)
521 {
522  warn_deprecated("isDuplicateReaction", "Now handled by "
523  "Kinetics::checkDuplicateStoich. To be removed after Cantera 2.2.");
524  map<int, doublereal>::const_iterator b = r1.begin(), e = r1.end();
525  int k1 = b->first;
526  // check for duplicate written in the same direction
527  doublereal ratio = 0.0;
528  if (r1[k1] && r2[k1]) {
529  ratio = r2[k1]/r1[k1];
530  ++b;
531  bool different = false;
532  for (; b != e; ++b) {
533  k1 = b->first;
534  if (!r1[k1] || !r2[k1] || fabs(r2[k1]/r1[k1] - ratio) > 1.e-8) {
535  different = true;
536  break;
537  }
538  }
539  if (!different) {
540  return ratio;
541  }
542  }
543 
544  // check for duplicate written in the reverse direction
545  b = r1.begin();
546  k1 = b->first;
547  if (r1[k1] == 0.0 || r2[-k1] == 0.0) {
548  return 0.0;
549  }
550  ratio = r2[-k1]/r1[k1];
551  ++b;
552  for (; b != e; ++b) {
553  k1 = b->first;
554  if (!r1[k1] || !r2[-k1] || fabs(r2[-k1]/r1[k1] - ratio) > 1.e-8) {
555  return 0.0;
556  }
557  }
558  return ratio;
559 }
560 
562  std::string default_phase, bool check_for_duplicates)
563 {
564  int itot = 0;
565  /*
566  * Search the children of the phase element for the
567  * XML element named reactionArray. If we can't find it,
568  * then return signaling having not found any reactions.
569  * Apparently, we allow multiple reactionArray elements here
570  * Each one will be processed sequentially, with the
571  * end result being purely additive.
572  */
573  vector<XML_Node*> rarrays = p.getChildren("reactionArray");
574  if (rarrays.empty()) {
575  kin.finalize();
576  return false;
577  }
578  for (size_t n = 0; n < rarrays.size(); n++) {
579  /*
580  * Go get a reference to the current XML element,
581  * reactionArray. We will process this element now.
582  */
583  const XML_Node& rxns = *rarrays[n];
584  /*
585  * The reactionArray element has an attribute called,
586  * datasrc. The value of the attribute is the XML
587  * element comprising the top of the
588  * tree of reactions for the phase.
589  * Find this datasrc element starting with the root
590  * of the current XML node.
591  */
592  const XML_Node* rdata = get_XML_Node(rxns["datasrc"], &rxns.root());
593  /*
594  * If the reactionArray element has a child element named "skip", and
595  * if the attribute of skip called "species" has a value of "undeclared",
596  * we will set rxnrule.skipUndeclaredSpecies to 'true'. rxnrule is
597  * passed to the routine that parses each individual reaction so that
598  * the parser will skip all reactions containing an undefined species
599  * without throwing an error.
600  *
601  * Similarly, an attribute named "third_bodies" with the value of
602  * "undeclared" will skip undeclared third body efficiencies (while
603  * retaining the reaction and any other efficiencies).
604  */
605  if (rxns.hasChild("skip")) {
606  const XML_Node& sk = rxns.child("skip");
607  if (sk["species"] == "undeclared") {
608  kin.skipUndeclaredSpecies(true);
609  }
610  if (sk["third_bodies"] == "undeclared") {
611  kin.skipUndeclaredThirdBodies(true);
612  }
613  }
614  /*
615  * Search for child elements called include. We only include
616  * a reaction if it's tagged by one of the include fields.
617  * Or, we include all reactions if there are no include fields.
618  */
619  vector<XML_Node*> incl = rxns.getChildren("include");
620  vector<XML_Node*> allrxns = rdata->getChildren("reaction");
621  // if no 'include' directive, then include all reactions
622  if (incl.empty()) {
623  for (size_t i = 0; i < allrxns.size(); i++) {
624  kin.addReaction(newReaction(*allrxns[i]));
625  ++itot;
626  }
627  } else {
628  for (size_t nii = 0; nii < incl.size(); nii++) {
629  const XML_Node& ii = *incl[nii];
630  string imin = ii["min"];
631  string imax = ii["max"];
632 
633  string::size_type iwild = string::npos;
634  if (imax == imin) {
635  iwild = imin.find("*");
636  if (iwild != string::npos) {
637  imin = imin.substr(0,iwild);
638  imax = imin;
639  }
640  }
641 
642  for (size_t i = 0; i < allrxns.size(); i++) {
643  const XML_Node* r = allrxns[i];
644  string rxid;
645  if (r) {
646  rxid = r->attrib("id");
647  if (iwild != string::npos) {
648  rxid = rxid.substr(0,iwild);
649  }
650  /*
651  * To decide whether the reaction is included or not
652  * we do a lexical min max and operation. This
653  * sometimes has surprising results.
654  */
655  if ((rxid >= imin) && (rxid <= imax)) {
656  kin.addReaction(newReaction(*r));
657  ++itot;
658  }
659  }
660  }
661  }
662  }
663  }
664 
665  if (check_for_duplicates) {
666  kin.checkDuplicates();
667  }
668  /*
669  * Finalize the installation of the kinetics, now that we know
670  * the true number of reactions in the mechanism, itot.
671  */
672  kin.finalize();
673 
674  return true;
675 }
676 
677 bool importKinetics(const XML_Node& phase, std::vector<ThermoPhase*> th,
678  Kinetics* k)
679 {
680  if (k == 0) {
681  return false;
682  }
683 
684  // This phase will be the owning phase for the kinetics operator
685  // For interfaces, it is the surface phase between two volumes.
686  // For homogeneous kinetics, it's the current volumetric phase.
687  string owning_phase = phase["id"];
688 
689  bool check_for_duplicates = false;
690  if (phase.parent()) {
691  if (phase.parent()->hasChild("validate")) {
692  const XML_Node& d = phase.parent()->child("validate");
693  if (d["reactions"] == "yes") {
694  check_for_duplicates = true;
695  }
696  }
697  }
698 
699  // if other phases are involved in the reaction mechanism,
700  // they must be listed in a 'phaseArray' child
701  // element. Homogeneous mechanisms do not need to include a
702  // phaseArray element.
703 
704  vector<string> phase_ids;
705  if (phase.hasChild("phaseArray")) {
706  const XML_Node& pa = phase.child("phaseArray");
707  getStringArray(pa, phase_ids);
708  }
709  phase_ids.push_back(owning_phase);
710 
711  int np = static_cast<int>(phase_ids.size());
712  int nt = static_cast<int>(th.size());
713 
714  // for each referenced phase, attempt to find its id among those
715  // phases specified.
716  bool phase_ok;
717 
718  string phase_id;
719  string msg = "";
720  for (int n = 0; n < np; n++) {
721  phase_id = phase_ids[n];
722  phase_ok = false;
723 
724  // loop over the supplied 'ThermoPhase' objects representing
725  // phases, to find an object with the same id.
726  for (int m = 0; m < nt; m++) {
727  if (th[m]->id() == phase_id) {
728  phase_ok = true;
729 
730  // if no phase with this id has been added to
731  //the kinetics manager yet, then add this one
732  if (k->phaseIndex(phase_id) == npos) {
733  k->addPhase(*th[m]);
734  }
735  }
736  msg += " "+th[m]->id();
737  }
738  if (!phase_ok) {
739  throw CanteraError("importKinetics",
740  "phase "+phase_id+" not found. Supplied phases are:"+msg);
741  }
742  }
743 
744  // allocates arrays, etc. Must be called after the phases have
745  // been added to 'kin', so that the number of species in each
746  // phase is known.
747  k->init();
748 
749  // Install the reactions.
750  return installReactionArrays(phase, *k, owning_phase, check_for_duplicates);
751 }
752 
753 bool buildSolutionFromXML(XML_Node& root, const std::string& id,
754  const std::string& nm, ThermoPhase* th, Kinetics* kin)
755 {
756  XML_Node* x;
757  x = get_XML_NameID(nm, string("#")+id, &root);
758  if (!x) {
759  return false;
760  }
761 
762  /*
763  * Fill in the ThermoPhase object by querying the
764  * const XML_Node tree located at x.
765  */
766  importPhase(*x, th);
767  /*
768  * Create a vector of ThermoPhase pointers of length 1
769  * having the current th ThermoPhase as the entry.
770  */
771  std::vector<ThermoPhase*> phases(1);
772  phases[0] = th;
773  /*
774  * Fill in the kinetics object k, by querying the
775  * const XML_Node tree located by x. The source terms and
776  * eventually the source term vector will be constructed
777  * from the list of ThermoPhases in the vector, phases.
778  */
779  importKinetics(*x, phases, kin);
780  return true;
781 }
782 
783 }
XML_Node * get_XML_Node(const std::string &file_ID, XML_Node *root)
This routine will locate an XML node in either the input XML tree or in another input file specified ...
Definition: global.cpp:214
doublereal fpValue(const std::string &val)
Translate a string into one doublereal value.
doublereal nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
Definition: Phase.cpp:243
const int PLOG_RXN
A pressure-dependent rate expression consisting of several Arrhenius rate expressions evaluated at di...
Definition: reaction_defs.h:50
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:39
static void getCoverageDependence(const XML_Node &node, thermo_t &surfphase, ReactionData &rdata)
Read the XML data concerning the coverage dependence of an interfacial reaction.
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data...
void checkRxnElementBalance(Kinetics &kin, const ReactionData &rdata, doublereal errorTolerance)
This function will check a specific reaction to see if the elements balance.
vector_fp falloffParameters
Values used in the falloff parameterization.
Definition: ReactionData.h:165
std::multimap< double, vector_fp > plogParameters
Arrhenius parameters for P-log reactions.
Definition: ReactionData.h:207
size_t nElements() const
Number of elements.
Definition: Phase.cpp:167
double chebPmin
Minimum pressure for Chebyshev fit.
Definition: ReactionData.h:211
int reactionType
Type of the reaction.
Definition: ReactionData.h:58
std::string attrib(const std::string &attr) const
Function returns the value of an attribute.
Definition: xml.cpp:527
thermo_t & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism...
Definition: Kinetics.h:285
vector_fp auxRateCoeffParameters
Vector of auxiliary rate coefficient parameters.
Definition: ReactionData.h:157
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
int falloffType
Type of falloff parameterization to use.
Definition: ReactionData.h:161
const int EDGE_RXN
A reaction occurring at a one-dimensional interface between two surface phases.
const int CHEBYSHEV_RXN
A general gas-phase pressure-dependent reaction where k(T,P) is defined in terms of a bivariate Cheby...
Definition: reaction_defs.h:56
vector_fp cov
Adjustments to the Arrhenius rate expression dependent on surface species coverages.
Definition: ReactionData.h:186
doublereal beta
Forward value of the apparent Electrochemical transfer coefficient.
Definition: ReactionData.h:201
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
vector_fp pstoich
Product stoichiometric coefficients, in the order given by products.
Definition: ReactionData.h:101
static void getStick(const XML_Node &node, Kinetics &kin, ReactionData &r, doublereal &A, doublereal &b, doublereal &E)
getStick() processes the XML element called Stick that specifies the sticking coefficient reaction...
std::string kineticsSpeciesName(size_t k) const
Return the name of the kth species in the kinetics manager.
Definition: Kinetics.cpp:359
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:100
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:78
static void getArrhenius(const XML_Node &node, int &labeled, doublereal &A, doublereal &b, doublereal &E)
getArrhenius() parses the XML element called Arrhenius.
const doublereal Pi
Pi.
Definition: ct_defs.h:51
const int CHEMACT_RXN
A chemical activation reaction.
Definition: reaction_defs.h:64
vector_fp rateCoeffParameters
Vector of rate coefficient parameters.
Definition: ReactionData.h:152
XML_Node & child(const size_t n) const
Return a changeable reference to the n'th child of the current node.
Definition: xml.cpp:573
std::map< size_t, doublereal > thirdBodyEfficiencies
Map of species index to third body efficiency.
Definition: ReactionData.h:107
virtual doublereal standardConcentration(size_t k=0) const
Return the standard concentration for the kth species.
Definition: ThermoPhase.h:447
bool installReactionArrays(const XML_Node &p, Kinetics &kin, std::string default_phase, bool check_for_duplicates)
Install information about reactions into the kinetics object, kin.
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:97
const int cSurf
A surface phase. Used by class SurfPhase.
Definition: mix_defs.h:40
const int FALLOFF_RXN
The general form for a gas-phase association or dissociation reaction, with a pressure-dependent rate...
Definition: reaction_defs.h:42
virtual void init()
Prepare the class for the addition of reactions.
Definition: Kinetics.h:778
size_t chebDegreeT
Degree of Chebyshev fit in T.
Definition: ReactionData.h:213
std::string equation
The reaction equation. Used only for display purposes.
Definition: ReactionData.h:170
std::vector< size_t > products
Indices of product species.
Definition: ReactionData.h:64
size_t speciesIndex(const std::string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition: Phase.cpp:257
size_t surfacePhaseIndex()
This returns the integer index of the phase which has ThermoPhase type cSurf.
Definition: Kinetics.h:260
bool importPhase(XML_Node &phase, ThermoPhase *th, SpeciesThermoFactory *spfactory)
Import a phase information into an empty ThermoPhase object.
virtual std::pair< size_t, size_t > checkDuplicates(bool throw_err=true) const
Check for duplicate reactions.
Definition: Kinetics.cpp:168
double chebPmax
Maximum pressure for Chebyshev fit.
Definition: ReactionData.h:212
bool getReagents(const XML_Node &rxn, Kinetics &kin, int rp, std::string default_phase, std::vector< size_t > &spnum, vector_fp &stoich, vector_fp &order, const ReactionRules &rules)
Get the reactants or products of a reaction.
bool importKinetics(const XML_Node &phase, std::vector< ThermoPhase * > th, Kinetics *k)
Import a reaction mechanism for a phase or an interface.
doublereal isDuplicateReaction(std::map< int, doublereal > &r1, std::map< int, doublereal > &r2)
This function returns a ratio if two reactions are duplicates of one another, and 0...
static void getFalloff(const XML_Node &f, ReactionData &rdata)
Get falloff parameters for a reaction.
std::string fp2str(const double x, const std::string &fmt)
Convert a double into a c++ string.
Definition: stringUtils.cpp:28
std::string name() const
Returns the name of the XML node.
Definition: xml.h:394
thermo_t & speciesPhase(const std::string &nm)
This function looks up the name of a species and returns a reference to the ThermoPhase object of the...
Definition: Kinetics.cpp:402
Public interface for kinetics managers.
Definition: Kinetics.h:128
int number
Index of this reaction within the mechanism.
Definition: ReactionData.h:61
int getPairs(const XML_Node &node, std::vector< std::string > &key, std::vector< std::string > &val)
This function interprets the value portion of an XML element as a series of "Pairs" separated by whit...
Definition: ctml.cpp:417
Rules for parsing and installing reactions.
doublereal default_3b_eff
The default third body efficiency for species not listed in thirdBodyEfficiencies.
Definition: ReactionData.h:180
size_t chebDegreeP
Degree of Chebyshev fit in P.
Definition: ReactionData.h:214
Intermediate class which stores data about a reaction and its rate parameterization before adding the...
Definition: ReactionData.h:22
std::string id() const
Return the string id for the phase.
Definition: Phase.cpp:147
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:99
shared_ptr< Reaction > newReaction(const XML_Node &rxn_node)
Create a new Reaction object for the reaction defined in rxn_node
Definition: Reaction.cpp:607
const int SURFACE_RXN
A reaction occurring on a surface.
Definition: reaction_defs.h:72
bool hasChild(const std::string &ch) const
Tests whether the current node has a child node with a particular name.
Definition: xml.cpp:563
const int cEdge
An edge between two 2D surfaces.
Definition: mix_defs.h:58
const int THREE_BODY_RXN
A gas-phase reaction that requires a third-body collision partner.
Definition: reaction_defs.h:36
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition: Kinetics.h:323
std::vector< size_t > reactants
Indices of reactant species.
Definition: ReactionData.h:63
vector_fp rstoich
Reactant stoichiometric coefficients, in the order given by reactants.
Definition: ReactionData.h:94
bool buildSolutionFromXML(XML_Node &root, const std::string &id, const std::string &nm, ThermoPhase *th, Kinetics *kin)
Build a single-phase ThermoPhase object with associated kinetics mechanism.
int rateCoeffType
Type of the rate coefficient for the forward rate constant.
Definition: ReactionData.h:147
Definitions of global routines for the importing of data from XML files (see Input File Handling)...
const vector_fp & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition: Phase.cpp:515
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:157
void getStringArray(const XML_Node &node, std::vector< std::string > &v)
This function interprets the value portion of an XML element as a string.
Definition: ctml.cpp:503
virtual int eosType() const
Equation of state type flag.
Definition: ThermoPhase.h:142
size_t speciesPhaseIndex(size_t k)
This function takes as an argument the kineticsSpecies index (i.e., the list index in the list of spe...
Definition: Kinetics.cpp:417
size_t getFloatArray(const XML_Node &node, std::vector< doublereal > &v, const bool convert, const std::string &unitsString, const std::string &nodeName)
This function reads the current node or a child node of the current node with the default name...
Definition: ctml.cpp:323
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:64
const int ELEMENTARY_RXN
A reaction with a rate coefficient that depends only on temperature and voltage that also obeys mass-...
Definition: reaction_defs.h:30
Contains declarations for string manipulation functions within Cantera.
doublereal getFloat(const XML_Node &parent, const std::string &name, const std::string &type)
Get a floating-point value from a child element.
Definition: ctml.cpp:194
std::string elementName(size_t m) const
Name of the element with index m.
Definition: Phase.cpp:186
static void getEfficiencies(const XML_Node &eff, Kinetics &kin, ReactionData &rdata, const ReactionRules &rules)
Get the enhanced collision efficiencies.
void getRateCoefficient(const XML_Node &kf, Kinetics &kin, ReactionData &rdata, const ReactionRules &rules)
Read the rate coefficient data from the XML file.
doublereal fp_value() const
Return the value of an XML node as a single double.
Definition: xml.cpp:481
XML_Node & root() const
Return the root of the current XML_Node tree.
Definition: xml.cpp:1095
virtual void finalize()
Finish adding reactions and prepare for use.
Definition: Kinetics.cpp:548
virtual void addPhase(thermo_t &thermo)
Add a phase to the kinetics manager object.
Definition: Kinetics.cpp:513
XML_Node * parent() const
Returns a pointer to the parent node of the current node.
Definition: xml.cpp:552
double chebTmax
Maximum temperature for Chebyshev fit.
Definition: ReactionData.h:210
vector_fp chebCoeffs
Chebyshev coefficients. length chebDegreeT * chebDegreeP.
Definition: ReactionData.h:217
vector_fp rorder
Reaction order with respect to each reactant species, in the order given by reactants.
Definition: ReactionData.h:72
void skipUndeclaredThirdBodies(bool skip)
Determine behavior when adding a new reaction that contains third-body efficiencies for species not d...
Definition: Kinetics.h:838
XML_Node * get_XML_NameID(const std::string &nameTarget, const std::string &file_ID, XML_Node *root)
This routine will locate an XML node in either the input XML tree or in another input file specified ...
Definition: global.cpp:252
virtual void addReaction(ReactionData &r)
Add a single reaction to the mechanism.
Definition: Kinetics.cpp:557
size_t nChildren(bool discardComments=false) const
Return the number of children.
Definition: xml.cpp:583
void getChildren(const std::string &name, std::vector< XML_Node * > &children) const
Get a vector of pointers to XML_Node containing all of the children of the current node which matches...
Definition: xml.cpp:915
size_t phaseIndex(const std::string &ph)
Return the phase index of a phase in the list of phases defined within the object.
Definition: Kinetics.h:247
void skipUndeclaredSpecies(bool skip)
Determine behavior when adding a new reaction that contains species not defined in any of the phases ...
Definition: Kinetics.h:829
double chebTmin
Minimum temperature for Chebyshev fit.
Definition: ReactionData.h:209