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