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