Cantera  3.1.0a1
Loading...
Searching...
No Matches
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"
17#include <boost/algorithm/string/predicate.hpp>
18#include <sstream>
19
20#include <boost/algorithm/string.hpp>
21
22namespace ba = boost::algorithm;
23
24namespace Cantera
25{
26
27Reaction::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
69Reaction::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
81Reaction::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
123void 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
159AnyMap 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
181void 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
229void 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
267void 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
309string 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
327string 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
345string Reaction::equation() const
346{
347 if (reversible) {
348 return reactantString() + " <=> " + productString();
349 } else {
350 return reactantString() + " => " + productString();
351 }
352}
353
354void 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
504string 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
523UnitStack 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
566void 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
576void 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
643bool 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
687bool 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
719ThirdBody::ThirdBody(const string& third_body)
720{
721 setName(third_body);
722}
723
724void 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
752ThirdBody::ThirdBody(const AnyMap& node)
753{
754 setParameters(node);
755}
756
757void 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
779void 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
792double ThirdBody::efficiency(const string& k) const
793{
794 return getValue(efficiencies, k, default_efficiency);
795}
796
797string ThirdBody::collider() const
798{
799 if (mass_action) {
800 return " + " + m_name;
801 }
802 return " (+" + m_name + ")";
803}
804
805bool 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
832unique_ptr<Reaction> newReaction(const string& type)
833{
834 return make_unique<Reaction>();
835}
836
837unique_ptr<Reaction> newReaction(const AnyMap& rxn_node, const Kinetics& kin)
838{
839 return make_unique<Reaction>(rxn_node, kin);
840}
841
842void 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
912vector<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.
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
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 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
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.
virtual Units standardConcentrationUnits() const
Returns the units of the "standard concentration" for this phase.
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
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
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
map< string, double > Composition
Map from string names to doubles.
Definition ct_defs.h:177
unique_ptr< Reaction > newReaction(const string &type)
Create a new empty Reaction object.
Definition Reaction.cpp:832
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...