Cantera  3.1.0a1
Reaction.cpp
Go to the documentation of this file.
1 /**
2  * @file Reaction.cpp
3  */
4 
5 // This file is part of Cantera. See License.txt in the top-level directory or
6 // at https://cantera.org/license.txt for license and copyright information.
7 
13 #include "cantera/base/Array.h"
14 #include "cantera/base/AnyMap.h"
15 #include "cantera/base/utilities.h"
17 #include <boost/algorithm/string/predicate.hpp>
18 #include <sstream>
19 
20 #include <boost/algorithm/string.hpp>
21 
22 namespace ba = boost::algorithm;
23 
24 namespace Cantera
25 {
26 
27 Reaction::Reaction(const Composition& reactants_,
28  const Composition& products_,
29  shared_ptr<ReactionRate> rate_,
30  shared_ptr<ThirdBody> tbody_)
31  : reactants(reactants_)
32  , products(products_)
33  , m_from_composition(true)
34  , m_third_body(tbody_)
35 {
36  if (reactants.count("M") || products.count("M")) {
37  throw CanteraError("Reaction::Reaction",
38  "Third body 'M' must not be included in either reactant or product maps.");
39  }
40  setRate(rate_);
41 
42  // set flags ensuring correct serialization output
43  Composition third;
44  for (const auto& [name, stoich] : reactants) {
45  if (products.count(name)) {
46  third[name] = products.at(name) - stoich;
47  }
48  }
49  if (tbody_) {
50  string name = tbody_->name();
51  if (reactants.count(name) && products.count(name)) {
52  throw CanteraError("Reaction::Reaction",
53  "'{}' not acting as third body collider must not be included in both "
54  "reactant and product maps.", name);
55  }
56  if (name != "M") {
57  m_third_body->explicit_3rd = true;
58  }
59  } else if (!tbody_ && third.size() == 1) {
60  // implicit third body
61  string name = third.begin()->first;
62  m_third_body = make_shared<ThirdBody>(name);
63  if (name != "M") {
64  m_third_body->explicit_3rd = true;
65  }
66  }
67 }
68 
69 Reaction::Reaction(const string& equation,
70  shared_ptr<ReactionRate> rate_,
71  shared_ptr<ThirdBody> tbody_)
72  : m_third_body(tbody_)
73 {
74  setEquation(equation);
75  setRate(rate_);
76  if (m_third_body && m_third_body->name() != "M") {
77  m_third_body->explicit_3rd = true;
78  }
79 }
80 
81 Reaction::Reaction(const AnyMap& node, const Kinetics& kin)
82 {
83  string rate_type = node.getString("type", "Arrhenius");
84  if (!kin.nPhases()) {
85  throw InputFileError("Reaction", node,
86  "Cannot instantiate Reaction with empty Kinetics object.");
87  }
88 
89  setParameters(node, kin);
90  size_t nDim = kin.thermo(0).nDim();
91  if (!valid()) {
92  // If the reaction isn't valid (for example, contains undefined species),
93  // setting up the rate constant won't work
94  return;
95  }
96  if (nDim == 3) {
97  if (ba::starts_with(rate_type, "three-body-")) {
98  AnyMap rateNode = node;
99  rateNode["type"] = rate_type.substr(11, rate_type.size() - 11);
100  setRate(newReactionRate(rateNode, calculateRateCoeffUnits(kin)));
101  } else {
102  setRate(newReactionRate(node, calculateRateCoeffUnits(kin)));
103  }
104  } else {
105  AnyMap rateNode = node;
106  if (rateNode.hasKey("rate-constant")) {
107  if (!ba::starts_with(rate_type, "interface-")) {
108  rateNode["type"] = "interface-" + rate_type;
109  }
110  } else if (node.hasKey("sticking-coefficient")) {
111  if (!ba::starts_with(rate_type, "sticking-")) {
112  rateNode["type"] = "sticking-" + rate_type;
113  }
114  } else {
115  throw InputFileError("Reaction::Reaction", input,
116  "Unable to infer interface reaction type.");
117  }
118  setRate(newReactionRate(rateNode, calculateRateCoeffUnits(kin)));
119  }
120  check();
121 }
122 
123 void Reaction::check()
124 {
125  if (!allow_nonreactant_orders) {
126  for (const auto& [name, order] : orders) {
127  if (reactants.find(name) == reactants.end()) {
128  throw InputFileError("Reaction::validate", input,
129  "Reaction order specified for non-reactant species '{}'", name);
130  }
131  }
132  }
133 
134  if (!allow_negative_orders) {
135  for (const auto& [name, order] : orders) {
136  if (order < 0.0) {
137  throw InputFileError("Reaction::validate", input,
138  "Negative reaction order specified for species '{}'", name);
139  }
140  }
141  }
142 
143  // If reaction orders are specified, then this reaction does not follow
144  // mass-action kinetics, and is not an elementary reaction. So check that it
145  // is not reversible, since computing the reverse rate from thermochemistry
146  // only works for elementary reactions.
147  if (reversible && !orders.empty()) {
148  throw InputFileError("Reaction::validate", input,
149  "Reaction orders may only be given for irreversible reactions");
150  }
151 
152  // Check reaction rate evaluator to ensure changes introduced after object
153  // instantiation are considered.
154  if (m_rate) {
155  m_rate->check(equation());
156  }
157 }
158 
159 AnyMap Reaction::parameters(bool withInput) const
160 {
161  AnyMap out;
162  getParameters(out);
163  if (withInput) {
164  out.update(input);
165  }
166 
167  static bool reg = AnyMap::addOrderingRules("Reaction",
168  {{"head", "type"},
169  {"head", "equation"},
170  {"tail", "duplicate"},
171  {"tail", "orders"},
172  {"tail", "negative-orders"},
173  {"tail", "nonreactant-orders"}
174  });
175  if (reg) {
176  out["__type__"] = "Reaction";
177  }
178  return out;
179 }
180 
181 void Reaction::getParameters(AnyMap& reactionNode) const
182 {
183  if (!m_rate) {
184  throw CanteraError("Reaction::getParameters",
185  "Serialization of empty Reaction object is not supported.");
186  }
187 
188  reactionNode["equation"] = equation();
189 
190  if (duplicate) {
191  reactionNode["duplicate"] = true;
192  }
193  if (orders.size()) {
194  reactionNode["orders"] = orders;
195  }
196  if (allow_negative_orders) {
197  reactionNode["negative-orders"] = true;
198  }
199  if (allow_nonreactant_orders) {
200  reactionNode["nonreactant-orders"] = true;
201  }
202 
203  reactionNode.update(m_rate->parameters());
204 
205  // strip information not needed for reconstruction
206  string rtype = reactionNode["type"].asString();
207  if (rtype == "pressure-dependent-Arrhenius") {
208  // skip
209  } else if (m_explicit_type && ba::ends_with(rtype, "Arrhenius")) {
210  // retain type information
211  if (m_third_body) {
212  reactionNode["type"] = "three-body";
213  } else {
214  reactionNode["type"] = "elementary";
215  }
216  } else if (ba::ends_with(rtype, "Arrhenius")) {
217  reactionNode.erase("type");
218  } else if (m_explicit_type) {
219  reactionNode["type"] = type();
220  } else if (ba::ends_with(rtype, "Blowers-Masel")) {
221  reactionNode["type"] = "Blowers-Masel";
222  }
223 
224  if (m_third_body) {
225  m_third_body->getParameters(reactionNode);
226  }
227 }
228 
229 void Reaction::setParameters(const AnyMap& node, const Kinetics& kin)
230 {
231  if (node.empty()) {
232  throw InputFileError("Reaction::setParameters", input,
233  "Cannot set reaction parameters from empty node.");
234  }
235 
236  input = node;
237  input.copyMetadata(node);
238  setEquation(node["equation"].asString(), &kin);
239  // Non-stoichiometric reaction orders
240  if (node.hasKey("orders")) {
241  for (const auto& [name, order] : node["orders"].asMap<double>()) {
242  orders[name] = order;
243  if (kin.kineticsSpeciesIndex(name) == npos) {
244  setValid(false);
245  }
246  }
247  }
248 
249  // Flags
250  id = node.getString("id", "");
251  duplicate = node.getBool("duplicate", false);
252  allow_negative_orders = node.getBool("negative-orders", false);
253  allow_nonreactant_orders = node.getBool("nonreactant-orders", false);
254 
255  if (m_third_body) {
256  m_third_body->setParameters(node);
257  if (m_third_body->name() == "M" && m_third_body->efficiencies.size() == 1) {
258  m_third_body->explicit_3rd = true;
259  }
260  } else if (node.hasKey("default-efficiency") || node.hasKey("efficiencies")) {
261  throw InputFileError("Reaction::setParameters", input,
262  "Reaction '{}' specifies efficiency parameters\n"
263  "but does not involve third body colliders.", equation());
264  }
265 }
266 
267 void Reaction::setRate(shared_ptr<ReactionRate> rate)
268 {
269  if (!rate) {
270  throw InputFileError("Reaction::setRate", input,
271  "Reaction rate for reaction '{}' must not be empty.", equation());
272  }
273  m_rate = rate;
274 
275  string rate_type = m_rate->type();
276  if (m_third_body) {
277  if (rate_type == "falloff" || rate_type == "chemically-activated") {
278  if (m_third_body->mass_action && !m_from_composition) {
279  throw InputFileError("Reaction::setRate", input,
280  "Third-body collider does not use '(+{})' notation.",
281  m_third_body->name());
282  }
283  m_third_body->mass_action = false;
284  } else if (rate_type == "Chebyshev") {
285  if (m_third_body->name() == "M") {
286  warn_deprecated("Chebyshev reaction equation", input, "Specifying 'M' "
287  "in the reaction equation for Chebyshev reactions is deprecated.");
288  m_third_body.reset();
289  }
290  } else if (rate_type == "pressure-dependent-Arrhenius") {
291  if (m_third_body->name() == "M") {
292  throw InputFileError("Reaction::setRate", input,
293  "Found superfluous '{}' in pressure-dependent-Arrhenius reaction.",
294  m_third_body->name());
295  }
296  }
297  } else {
298  if (rate_type == "falloff" || rate_type == "chemically-activated") {
299  if (!m_from_composition) {
300  throw InputFileError("Reaction::setRate", input,
301  "Reaction equation for falloff reaction '{}'\n does not "
302  "contain valid pressure-dependent third body", equation());
303  }
304  m_third_body = make_shared<ThirdBody>("(+M)");
305  }
306  }
307 }
308 
309 string Reaction::reactantString() const
310 {
311  std::ostringstream result;
312  for (auto iter = reactants.begin(); iter != reactants.end(); ++iter) {
313  if (iter != reactants.begin()) {
314  result << " + ";
315  }
316  if (iter->second != 1.0) {
317  result << iter->second << " ";
318  }
319  result << iter->first;
320  }
321  if (m_third_body) {
322  result << m_third_body->collider();
323  }
324  return result.str();
325 }
326 
327 string Reaction::productString() const
328 {
329  std::ostringstream result;
330  for (auto iter = products.begin(); iter != products.end(); ++iter) {
331  if (iter != products.begin()) {
332  result << " + ";
333  }
334  if (iter->second != 1.0) {
335  result << iter->second << " ";
336  }
337  result << iter->first;
338  }
339  if (m_third_body) {
340  result << m_third_body->collider();
341  }
342  return result.str();
343 }
344 
345 string Reaction::equation() const
346 {
347  if (reversible) {
348  return reactantString() + " <=> " + productString();
349  } else {
350  return reactantString() + " => " + productString();
351  }
352 }
353 
354 void Reaction::setEquation(const string& equation, const Kinetics* kin)
355 {
356  parseReactionEquation(*this, equation, input, kin);
357 
358  string rate_type = input.getString("type", "");
359  if (ba::starts_with(rate_type, "three-body")) {
360  // state type when serializing
361  m_explicit_type = true;
362  } else if (rate_type == "elementary") {
363  // user override
364  m_explicit_type = true;
365  return;
366  } else if (kin && kin->thermo(0).nDim() != 3) {
367  // interface reactions
368  return;
369  }
370 
371  string third_body;
372  size_t count = 0;
373  size_t countM = 0;
374  for (const auto& [name, stoich] : reactants) {
375  // detect explicitly specified collision partner
376  if (products.count(name)) {
377  third_body = name;
378  size_t generic = third_body == "M"
379  || third_body == "(+M)" || third_body == "(+ M)";
380  count++;
381  countM += generic;
382  if (stoich > 1 && products[third_body] > 1) {
383  count++;
384  countM += generic;
385  }
386  }
387  }
388 
389  if (count == 0) {
390  if (ba::starts_with(rate_type, "three-body")) {
391  throw InputFileError("Reaction::setEquation", input,
392  "Reactants for reaction '{}'\n"
393  "do not contain a valid third body collider", equation);
394  }
395  return;
396 
397  } else if (countM > 1) {
398  throw InputFileError("Reaction::setEquation", input,
399  "Multiple generic third body colliders 'M' are not supported", equation);
400 
401  } else if (count > 1) {
402  // equations with more than one explicit third-body collider are handled as a
403  // regular elementary reaction unless the equation contains a generic third body
404  if (countM) {
405  // generic collider 'M' is selected as third body
406  third_body = "M";
407  } else if (m_third_body) {
408  // third body is defined as explicit object
409  auto& effs = m_third_body->efficiencies;
410  if (effs.size() != 1 || !reactants.count(effs.begin()->first)) {
411  throw InputFileError("Reaction::setEquation", input,
412  "Detected ambiguous third body colliders in reaction '{}'\n"
413  "ThirdBody object needs to specify a single species", equation);
414  }
415  third_body = effs.begin()->first;
416  m_third_body->explicit_3rd = true;
417  } else if (input.hasKey("efficiencies")) {
418  // third body is implicitly defined by efficiency
419  auto effs = input["efficiencies"].asMap<double>();
420  if (effs.size() != 1 || !reactants.count(effs.begin()->first)) {
421  throw InputFileError("Reaction::setEquation", input,
422  "Detected ambiguous third body colliders in reaction '{}'\n"
423  "Collision efficiencies need to specify single species", equation);
424  }
425  third_body = effs.begin()->first;
426  m_third_body = make_shared<ThirdBody>(third_body);
427  m_third_body->explicit_3rd = true;
428  } else if (input.hasKey("default-efficiency")) {
429  // insufficient disambiguation of third bodies
430  throw InputFileError("Reaction::setEquation", input,
431  "Detected ambiguous third body colliders in reaction '{}'\n"
432  "Third-body definition requires specification of efficiencies",
433  equation);
434  } else if (ba::starts_with(rate_type, "three-body")) {
435  // no disambiguation of third bodies
436  throw InputFileError("Reaction::setEquation", input,
437  "Detected ambiguous third body colliders in reaction '{}'\n"
438  "A valid ThirdBody or collision efficiency definition is required",
439  equation);
440  } else {
441  return;
442  }
443 
444  } else if (third_body != "M" && !ba::starts_with(rate_type, "three-body")
445  && !ba::starts_with(third_body, "(+"))
446  {
447  // check for conditions of three-body reactions:
448  // - integer stoichiometric conditions
449  // - either reactant or product side involves exactly three species
450  size_t nreac = 0;
451  size_t nprod = 0;
452 
453  // ensure that all reactants have integer stoichiometric coefficients
454  for (const auto& [name, stoich] : reactants) {
455  if (trunc(stoich) != stoich) {
456  return;
457  }
458  nreac += static_cast<size_t>(stoich);
459  }
460 
461  // ensure that all products have integer stoichiometric coefficients
462  for (const auto& [name, stoich] : products) {
463  if (trunc(stoich) != stoich) {
464  return;
465  }
466  nprod += static_cast<size_t>(stoich);
467  }
468 
469  // either reactant or product side involves exactly three species
470  if (nreac != 3 && nprod != 3) {
471  return;
472  }
473  }
474 
475  if (m_third_body) {
476  string tName = m_third_body->name();
477  if (tName != third_body && third_body != "M" && tName != "M") {
478  throw InputFileError("Reaction::setEquation", input,
479  "Detected incompatible third body colliders in reaction '{}'\n"
480  "ThirdBody definition does not match equation", equation);
481  }
482  m_third_body->setName(third_body);
483  } else {
484  m_third_body = make_shared<ThirdBody>(third_body);
485  }
486 
487  // adjust reactant coefficients
488  auto reac = reactants.find(third_body);
489  if (trunc(reac->second) != 1) {
490  reac->second -= 1.;
491  } else {
492  reactants.erase(reac);
493  }
494 
495  // adjust product coefficients
496  auto prod = products.find(third_body);
497  if (trunc(prod->second) != 1) {
498  prod->second -= 1.;
499  } else {
500  products.erase(prod);
501  }
502 }
503 
504 string Reaction::type() const
505 {
506  if (!m_rate) {
507  throw CanteraError("Reaction::type", "Empty Reaction does not have a type");
508  }
509 
510  string rate_type = m_rate->type();
511  string sub_type = m_rate->subType();
512  if (sub_type != "") {
513  return rate_type + "-" + sub_type;
514  }
515 
516  if (m_third_body) {
517  return "three-body-" + rate_type;
518  }
519 
520  return rate_type;
521 }
522 
523 UnitStack Reaction::calculateRateCoeffUnits(const Kinetics& kin)
524 {
525  if (!valid()) {
526  // If a reaction is invalid because of missing species in the Kinetics
527  // object, determining the units of the rate coefficient is impossible.
528  return UnitStack({});
529  }
530 
531  // Determine the units of the rate coefficient
532  UnitStack rate_units(kin.thermo(0).standardConcentrationUnits());
533 
534  // Set output units to standardConcentrationUnits per second
535  rate_units.join(1.);
536  rate_units.update(Units(1.0, 0, 0, -1), 1.);
537 
538  for (const auto& [name, order] : orders) {
539  const auto& phase = kin.speciesPhase(name);
540  // Account for specified reaction orders
541  rate_units.update(phase.standardConcentrationUnits(), -order);
542  }
543  for (const auto& [name, stoich] : reactants) {
544  // Order for each reactant is the reactant stoichiometric coefficient,
545  // unless already overridden by user-specified orders
546  if (name == "M" || ba::starts_with(name, "(+")) {
547  // calculateRateCoeffUnits may be called before these pseudo-species
548  // have been stripped from the reactants
549  continue;
550  } else if (orders.find(name) == orders.end()) {
551  const auto& phase = kin.speciesPhase(name);
552  // Account for each reactant species
553  rate_units.update(phase.standardConcentrationUnits(), -stoich);
554  }
555  }
556 
557  if (m_third_body && m_third_body->mass_action) {
558  // Account for third-body collision partner as the last entry
559  rate_units.join(-1);
560  }
561 
562  Reaction::rate_units = rate_units.product();
563  return rate_units;
564 }
565 
566 void updateUndeclared(vector<string>& undeclared,
567  const Composition& comp, const Kinetics& kin)
568 {
569  for (const auto& [name, stoich]: comp) {
570  if (kin.kineticsSpeciesIndex(name) == npos) {
571  undeclared.emplace_back(name);
572  }
573  }
574 }
575 
576 void Reaction::checkBalance(const Kinetics& kin) const
577 {
578  Composition balr, balp;
579 
580  // iterate over products and reactants
581  for (const auto& [name, stoich] : products) {
582  const ThermoPhase& ph = kin.speciesPhase(name);
583  size_t k = ph.speciesIndex(name);
584  for (size_t m = 0; m < ph.nElements(); m++) {
585  balr[ph.elementName(m)] = 0.0; // so that balr contains all species
586  balp[ph.elementName(m)] += stoich * ph.nAtoms(k, m);
587  }
588  }
589  for (const auto& [name, stoich] : reactants) {
590  const ThermoPhase& ph = kin.speciesPhase(name);
591  size_t k = ph.speciesIndex(name);
592  for (size_t m = 0; m < ph.nElements(); m++) {
593  balr[ph.elementName(m)] += stoich * ph.nAtoms(k, m);
594  }
595  }
596 
597  string msg;
598  bool ok = true;
599  for (const auto& [elem, balance] : balr) {
600  double elemsum = balr[elem] + balp[elem];
601  double elemdiff = fabs(balp[elem] - balr[elem]);
602  if (elemsum > 0.0 && elemdiff / elemsum > 1e-4) {
603  ok = false;
604  msg += fmt::format(" {} {} {}\n",
605  elem, balr[elem], balp[elem]);
606  }
607  }
608  if (!ok) {
609  throw InputFileError("Reaction::checkBalance", input,
610  "The following reaction is unbalanced: {}\n"
611  " Element Reactants Products\n{}",
612  equation(), msg);
613  }
614 
615  if (kin.thermo(0).nDim() == 3) {
616  return;
617  }
618 
619  // Check that the number of surface sites is balanced
620  double reac_sites = 0.0;
621  double prod_sites = 0.0;
622  auto& surf = dynamic_cast<const SurfPhase&>(kin.thermo(0));
623  for (const auto& [name, stoich] : reactants) {
624  size_t k = surf.speciesIndex(name);
625  if (k != npos) {
626  reac_sites += stoich * surf.size(k);
627  }
628  }
629  for (const auto& [name, stoich] : products) {
630  size_t k = surf.speciesIndex(name);
631  if (k != npos) {
632  prod_sites += stoich * surf.size(k);
633  }
634  }
635  if (fabs(reac_sites - prod_sites) > 1e-5 * (reac_sites + prod_sites)) {
636  throw InputFileError("Reaction::checkBalance", input,
637  "Number of surface sites not balanced in reaction {}.\n"
638  "Reactant sites: {}\nProduct sites: {}",
639  equation(), reac_sites, prod_sites);
640  }
641 }
642 
643 bool Reaction::checkSpecies(const Kinetics& kin) const
644 {
645  // Check for undeclared species
646  vector<string> undeclared;
647  updateUndeclared(undeclared, reactants, kin);
648  updateUndeclared(undeclared, products, kin);
649  if (!undeclared.empty()) {
650  if (kin.skipUndeclaredSpecies()) {
651  return false;
652  } else {
653  throw InputFileError("Reaction::checkSpecies", input, "Reaction '{}'\n"
654  "contains undeclared species: '{}'",
655  equation(), ba::join(undeclared, "', '"));
656  }
657  }
658 
659  undeclared.clear();
660  updateUndeclared(undeclared, orders, kin);
661  if (!undeclared.empty()) {
662  if (kin.skipUndeclaredSpecies()) {
663  return false;
664  } else {
665  if (input.hasKey("orders")) {
666  throw InputFileError("Reaction::checkSpecies", input["orders"],
667  "Reaction '{}'\n"
668  "defines reaction orders for undeclared species: '{}'",
669  equation(), ba::join(undeclared, "', '"));
670  }
671  // Error for empty input AnyMap (that is, XML)
672  throw InputFileError("Reaction::checkSpecies", input, "Reaction '{}'\n"
673  "defines reaction orders for undeclared species: '{}'",
674  equation(), ba::join(undeclared, "', '"));
675  }
676  }
677 
678  if (m_third_body) {
679  return m_third_body->checkSpecies(*this, kin);
680  }
681 
682  checkBalance(kin);
683 
684  return true;
685 }
686 
687 bool Reaction::usesElectrochemistry(const Kinetics& kin) const
688 {
689  // Check electrochemistry
690  vector<double> e_counter(kin.nPhases(), 0.0);
691 
692  // Find the number of electrons in the products for each phase
693  for (const auto& [name, stoich] : products) {
694  size_t kkin = kin.kineticsSpeciesIndex(name);
695  size_t i = kin.speciesPhaseIndex(kkin);
696  size_t kphase = kin.thermo(i).speciesIndex(name);
697  e_counter[i] += stoich * kin.thermo(i).charge(kphase);
698  }
699 
700  // Subtract the number of electrons in the reactants for each phase
701  for (const auto& [name, stoich] : reactants) {
702  size_t kkin = kin.kineticsSpeciesIndex(name);
703  size_t i = kin.speciesPhaseIndex(kkin);
704  size_t kphase = kin.thermo(i).speciesIndex(name);
705  e_counter[i] -= stoich * kin.thermo(i).charge(kphase);
706  }
707 
708  // If the electrons change phases then the reaction is electrochemical
709  for (double delta_e : e_counter) {
710  if (std::abs(delta_e) > 1e-4) {
711  return true;
712  }
713  }
714 
715  return false;
716 }
717 
718 
719 ThirdBody::ThirdBody(const string& third_body)
720 {
721  setName(third_body);
722 }
723 
724 void ThirdBody::setName(const string& third_body)
725 {
726  string name = third_body;
727  if (ba::starts_with(third_body, "(+ ")) {
728  mass_action = false;
729  name = third_body.substr(3, third_body.size() - 4);
730  } else if (ba::starts_with(third_body, "(+")) {
731  mass_action = false;
732  name = third_body.substr(2, third_body.size() - 3);
733  }
734 
735  if (name == m_name) {
736  return;
737  }
738  if (name == "M" && efficiencies.size() == 1) {
739  // revert from explicit name to generic collider
740  m_name = name;
741  return;
742  }
743  if (efficiencies.size()) {
744  throw CanteraError("ThirdBody::setName",
745  "Conflicting efficiency definition for explicit third body '{}'", name);
746  }
747  m_name = name;
748  default_efficiency = 0.;
749  efficiencies[m_name] = 1.;
750 }
751 
752 ThirdBody::ThirdBody(const AnyMap& node)
753 {
754  setParameters(node);
755 }
756 
757 void ThirdBody::setParameters(const AnyMap& node)
758 {
759  if (node.hasKey("default-efficiency")) {
760  double value = node["default-efficiency"].asDouble();
761  if (m_name != "M" && value != 0.) {
762  throw InputFileError("ThirdBody::setParameters", node["default-efficiency"],
763  "Invalid default efficiency for explicit collider {};\n"
764  "value is optional and/or needs to be zero", m_name);
765  }
766  default_efficiency = value;
767  }
768  if (node.hasKey("efficiencies")) {
769  efficiencies = node["efficiencies"].asMap<double>();
770  }
771  if (m_name != "M"
772  && (efficiencies.size() != 1 || efficiencies.begin()->first != m_name))
773  {
774  throw InputFileError("ThirdBody::setParameters", node,
775  "Detected incompatible third body colliders definitions");
776  }
777 }
778 
779 void ThirdBody::getParameters(AnyMap& node) const
780 {
781  if (m_name == "M" || explicit_3rd) {
782  if (efficiencies.size()) {
783  node["efficiencies"] = efficiencies;
784  node["efficiencies"].setFlowStyle();
785  }
786  if (default_efficiency != 1.0 && !explicit_3rd) {
787  node["default-efficiency"] = default_efficiency;
788  }
789  }
790 }
791 
792 double ThirdBody::efficiency(const string& k) const
793 {
794  return getValue(efficiencies, k, default_efficiency);
795 }
796 
797 string ThirdBody::collider() const
798 {
799  if (mass_action) {
800  return " + " + m_name;
801  }
802  return " (+" + m_name + ")";
803 }
804 
805 bool ThirdBody::checkSpecies(const Reaction& rxn, const Kinetics& kin) const
806 {
807  vector<string> undeclared;
808  updateUndeclared(undeclared, efficiencies, kin);
809 
810  if (!undeclared.empty()) {
811  if (!kin.skipUndeclaredThirdBodies()) {
812  if (rxn.input.hasKey("efficiencies")) {
813  throw InputFileError("ThirdBody::checkSpecies",
814  rxn.input["efficiencies"], "Reaction '{}'\n"
815  "defines third-body efficiencies for undeclared species: '{}'",
816  rxn.equation(), ba::join(undeclared, "', '"));
817  }
818  // Error for specified ThirdBody or empty input AnyMap
819  throw InputFileError("ThirdBody::checkSpecies", rxn.input, "Reaction '{}'\n"
820  "is a three-body reaction with undeclared species: '{}'",
821  rxn.equation(), ba::join(undeclared, "', '"));
822  } else if (kin.skipUndeclaredThirdBodies() && m_name != "M") {
823  // Prevent addition of reaction silently as "skip-undeclared-third-bodies"
824  // is set to true
825  return false;
826  }
827  }
828  return true;
829 }
830 
831 
832 unique_ptr<Reaction> newReaction(const string& type)
833 {
834  return make_unique<Reaction>();
835 }
836 
837 unique_ptr<Reaction> newReaction(const AnyMap& rxn_node, const Kinetics& kin)
838 {
839  return make_unique<Reaction>(rxn_node, kin);
840 }
841 
842 void parseReactionEquation(Reaction& R, const string& equation,
843  const AnyBase& reactionNode, const Kinetics* kin)
844 {
845  // Parse the reaction equation to determine participating species and
846  // stoichiometric coefficients
847  vector<string> tokens;
848  tokenizeString(equation, tokens);
849  tokens.push_back("+"); // makes parsing last species not a special case
850 
851  size_t last_used = npos; // index of last-used token
852  bool reactants = true;
853  for (size_t i = 1; i < tokens.size(); i++) {
854  if (tokens[i] == "+" || ba::starts_with(tokens[i], "(+") ||
855  tokens[i] == "<=>" || tokens[i] == "=" || tokens[i] == "=>") {
856  string species = tokens[i-1];
857 
858  double stoich = 1.0;
859  bool mass_action = true;
860  if (last_used != npos && tokens[last_used] == "(+"
861  && ba::ends_with(species, ")")) {
862  // Falloff third body with space, such as "(+ M)"
863  mass_action = false;
864  species = "(+" + species;
865  } else if (last_used == i - 1 && ba::starts_with(species, "(+")
866  && ba::ends_with(species, ")")) {
867  // Falloff 3rd body written without space, such as "(+M)"
868  mass_action = false;
869  } else if (last_used == i - 2) {
870  // Species with no stoich. coefficient
871  } else if (last_used == i - 3) {
872  // Stoich. coefficient and species
873  try {
874  stoich = fpValueCheck(tokens[i-2]);
875  } catch (CanteraError& err) {
876  throw InputFileError("parseReactionEquation", reactionNode,
877  err.getMessage());
878  }
879  } else {
880  throw InputFileError("parseReactionEquation", reactionNode,
881  "Error parsing reaction string '{}'.\n"
882  "Current token: '{}'\nlast_used: '{}'",
883  equation, tokens[i],
884  (last_used == npos) ? "n/a" : tokens[last_used]);
885  }
886  if (!kin || (kin->kineticsSpeciesIndex(species) == npos
887  && mass_action && species != "M"))
888  {
889  R.setValid(false);
890  }
891 
892  if (reactants) {
893  R.reactants[species] += stoich;
894  } else {
895  R.products[species] += stoich;
896  }
897 
898  last_used = i;
899  }
900 
901  // Tokens after this point are part of the products string
902  if (tokens[i] == "<=>" || tokens[i] == "=") {
903  R.reversible = true;
904  reactants = false;
905  } else if (tokens[i] == "=>") {
906  R.reversible = false;
907  reactants = false;
908  }
909  }
910 }
911 
912 vector<shared_ptr<Reaction>> getReactions(const AnyValue& items, Kinetics& kinetics)
913 {
914  vector<shared_ptr<Reaction>> all_reactions;
915  for (const auto& node : items.asVector<AnyMap>()) {
916  auto R = make_shared<Reaction>(node, kinetics);
917  R->check();
918  R->validate(kinetics);
919  if (R->valid() && R->checkSpecies(kinetics)) {
920  all_reactions.emplace_back(R);
921  }
922  }
923  return all_reactions;
924 }
925 
926 }
Header file for class Cantera::Array2D.
Base class for kinetics managers and also contains the kineticsmgr module documentation (see Kinetics...
Factory class for reaction rate objects.
Header for a simple thermodynamics model of a surface phase derived from ThermoPhase,...
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
Base class defining common data possessed by both AnyMap and AnyValue objects.
Definition: AnyMap.h:34
A map of string keys to values whose type can vary at runtime.
Definition: AnyMap.h:427
void copyMetadata(const AnyMap &other)
Copy metadata including input line/column from an existing AnyMap.
Definition: AnyMap.cpp:1493
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition: AnyMap.cpp:1423
bool empty() const
Return boolean indicating whether AnyMap is empty.
Definition: AnyMap.cpp:1418
void setFlowStyle(bool flow=true)
Use "flow" style when outputting this AnyMap to YAML.
Definition: AnyMap.cpp:1726
bool getBool(const string &key, bool default_) const
If key exists, return it as a bool, otherwise return default_.
Definition: AnyMap.cpp:1515
const string & getString(const string &key, const string &default_) const
If key exists, return it as a string, otherwise return default_.
Definition: AnyMap.cpp:1530
void erase(const string &key)
Erase the value held by key.
Definition: AnyMap.cpp:1428
void update(const AnyMap &other, bool keepExisting=true)
Add items from other to this AnyMap.
Definition: AnyMap.cpp:1438
A wrapper for a variable whose type is determined at runtime.
Definition: AnyMap.h:86
const vector< T > & asVector(size_t nMin=npos, size_t nMax=npos) const
Return the held value, if it is a vector of type T.
Definition: AnyMap.inl.h:109
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
virtual string getMessage() const
Method overridden by derived classes to format the error message.
Error thrown for problems processing information contained in an AnyMap or AnyValue.
Definition: AnyMap.h:738
Public interface for kinetics managers.
Definition: Kinetics.h:125
size_t nPhases() const
The number of phases participating in the reaction mechanism.
Definition: Kinetics.h:184
void skipUndeclaredSpecies(bool skip)
Determine behavior when adding a new reaction that contains species not defined in any of the phases ...
Definition: Kinetics.h:1326
void skipUndeclaredThirdBodies(bool skip)
Determine behavior when adding a new reaction that contains third-body efficiencies for species not d...
Definition: Kinetics.h:1338
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition: Kinetics.h:276
ThermoPhase & speciesPhase(const string &nm)
This function looks up the name of a species and returns a reference to the ThermoPhase object of the...
Definition: Kinetics.cpp:260
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
Definition: Kinetics.h:242
size_t speciesPhaseIndex(size_t k) const
This function takes as an argument the kineticsSpecies index (that is, the list index in the list of ...
Definition: Kinetics.cpp:281
size_t nDim() const
Returns the number of spatial dimensions (1, 2, or 3)
Definition: Phase.h:546
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition: Phase.cpp:129
double nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
Definition: Phase.cpp:103
size_t nElements() const
Number of elements.
Definition: Phase.cpp:30
string elementName(size_t m) const
Name of the element with index m.
Definition: Phase.cpp:49
double charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Definition: Phase.h:538
Abstract base class which stores data about a reaction and its rate parameterization so that it can b...
Definition: Reaction.h:25
bool reversible
True if the current reaction is reversible. False otherwise.
Definition: Reaction.h:126
string equation() const
The chemical equation for this reaction.
Definition: Reaction.cpp:345
Composition products
Product species and stoichiometric coefficients.
Definition: Reaction.h:114
Composition reactants
Reactant species and stoichiometric coefficients.
Definition: Reaction.h:111
AnyMap input
Input data used for specific models.
Definition: Reaction.h:139
void setValid(bool valid)
Set validity flag of reaction.
Definition: Reaction.h:89
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
Definition: SurfPhase.h:98
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:390
virtual Units standardConcentrationUnits() const
Returns the units of the "standard concentration" for this phase.
Definition: ThermoPhase.cpp:49
A representation of the units associated with a dimensional quantity.
Definition: Units.h:35
double fpValueCheck(const string &val)
Translate a string into one double value, with error checking.
void tokenizeString(const string &in_val, vector< string > &v)
This function separates a string up into tokens according to the location of white space.
shared_ptr< ReactionRate > newReactionRate(const string &type)
Create a new empty ReactionRate object.
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:180
void parseReactionEquation(Reaction &R, const string &equation, const AnyBase &reactionNode, const Kinetics *kin)
Parse reaction equation.
Definition: Reaction.cpp:842
vector< shared_ptr< Reaction > > getReactions(const AnyValue &items, Kinetics &kinetics)
Create Reaction objects for each item (an AnyMap) in items.
Definition: Reaction.cpp:912
unique_ptr< Reaction > newReaction(const AnyMap &rxn_node, const Kinetics &kin)
Create a new Reaction object using the specified parameters.
Definition: Reaction.cpp:837
const U & getValue(const map< T, U > &m, const T &key, const U &default_val)
Const accessor for a value in a map.
Definition: utilities.h:190
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
Definition: AnyMap.cpp:1926
map< string, double > Composition
Map from string names to doubles.
Definition: ct_defs.h:177
Contains declarations for string manipulation functions within Cantera.
Unit aggregation utility.
Definition: Units.h:105
void update(const Units &units, double exponent)
Update exponent of item with matching units; if it does not exist, add unit-exponent pair at end of s...
Definition: Units.cpp:362
Units product() const
Calculate product of units-exponent stack.
Definition: Units.cpp:377
void join(double exponent)
Join (update) exponent of standard units, where the updated exponent is the sum of the pre-existing e...
Definition: Units.cpp:352
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...