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
8 #include "cantera/kinetics/Reaction.h"
9 #include "cantera/kinetics/ReactionFactory.h"
10 #include "cantera/kinetics/ReactionRateFactory.h"
11 #include "cantera/kinetics/FalloffFactory.h"
12 #include "cantera/kinetics/Kinetics.h"
13 #include "cantera/thermo/ThermoPhase.h"
14 #include "cantera/base/ctml.h"
15 #include "cantera/base/Array.h"
16 #include "cantera/base/AnyMap.h"
17 #include "cantera/base/utilities.h"
18 #include "cantera/base/stringUtils.h"
19 #include <sstream>
20 #include <set>
21
22 #include <boost/algorithm/string.hpp>
23
24 namespace ba = boost::algorithm;
25
26 namespace Cantera
27 {
28
Reaction()29 Reaction::Reaction()
30 : reaction_type(NONE)
31 , reversible(true)
32 , duplicate(false)
33 , allow_nonreactant_orders(false)
34 , allow_negative_orders(false)
35 , rate_units(0.0)
36 , m_valid(true)
37 {
38 }
39
Reaction(const Composition & reactants_,const Composition & products_)40 Reaction::Reaction(const Composition& reactants_,
41 const Composition& products_)
42 : reaction_type(NONE)
43 , reactants(reactants_)
44 , products(products_)
45 , reversible(true)
46 , duplicate(false)
47 , allow_nonreactant_orders(false)
48 , allow_negative_orders(false)
49 , rate_units(0.0)
50 , m_valid(true)
51 {
52 }
53
Reaction(int type)54 Reaction::Reaction(int type)
55 : reaction_type(type)
56 , reversible(true)
57 , duplicate(false)
58 , allow_nonreactant_orders(false)
59 , allow_negative_orders(false)
60 , rate_units(0.0)
61 , m_valid(true)
62 {
63 warn_deprecated("Reaction::Reaction()",
64 "To be removed after Cantera 2.6. Use constructor without parameter "
65 "'type' instead.");
66 }
67
Reaction(int type,const Composition & reactants_,const Composition & products_)68 Reaction::Reaction(int type, const Composition& reactants_,
69 const Composition& products_)
70 : reaction_type(type)
71 , reactants(reactants_)
72 , products(products_)
73 , reversible(true)
74 , duplicate(false)
75 , allow_nonreactant_orders(false)
76 , allow_negative_orders(false)
77 , rate_units(0.0)
78 , m_valid(true)
79 {
80 warn_deprecated("Reaction::Reaction()",
81 "To be removed after Cantera 2.6. Use constructor without parameter "
82 "'type' instead.");
83 }
84
validate()85 void Reaction::validate()
86 {
87 if (!allow_nonreactant_orders) {
88 for (const auto& order : orders) {
89 if (reactants.find(order.first) == reactants.end()) {
90 throw InputFileError("Reaction::validate", input,
91 "Reaction order specified for non-reactant species '{}'",
92 order.first);
93 }
94 }
95 }
96
97 if (!allow_negative_orders) {
98 for (const auto& order : orders) {
99 if (order.second < 0.0) {
100 throw InputFileError("Reaction::validate", input,
101 "Negative reaction order specified for species '{}'",
102 order.first);
103 }
104 }
105 }
106
107 // If reaction orders are specified, then this reaction does not follow
108 // mass-action kinetics, and is not an elementary reaction. So check that it
109 // is not reversible, since computing the reverse rate from thermochemistry
110 // only works for elementary reactions.
111 if (reversible && !orders.empty()) {
112 throw InputFileError("Reaction::validate", input,
113 "Reaction orders may only be given for irreversible reactions");
114 }
115
116 // Call validation of reaction rate evaluator
117 if (!usesLegacy()) {
118 m_rate->validate(equation());
119 }
120 }
121
parameters(bool withInput) const122 AnyMap Reaction::parameters(bool withInput) const
123 {
124 AnyMap out;
125 getParameters(out);
126 if (withInput) {
127 out.update(input);
128 }
129
130 static bool reg = AnyMap::addOrderingRules("Reaction",
131 {{"head", "type"},
132 {"head", "equation"},
133 {"tail", "duplicate"},
134 {"tail", "orders"},
135 {"tail", "negative-orders"},
136 {"tail", "nonreactant-orders"}
137 });
138 if (reg) {
139 out["__type__"] = "Reaction";
140 }
141 return out;
142 }
143
getParameters(AnyMap & reactionNode) const144 void Reaction::getParameters(AnyMap& reactionNode) const
145 {
146 reactionNode["equation"] = equation();
147
148 if (duplicate) {
149 reactionNode["duplicate"] = true;
150 }
151 if (orders.size()) {
152 reactionNode["orders"] = orders;
153 }
154 if (allow_negative_orders) {
155 reactionNode["negative-orders"] = true;
156 }
157 if (allow_nonreactant_orders) {
158 reactionNode["nonreactant-orders"] = true;
159 }
160
161 if (m_rate) {
162 reactionNode.update(m_rate->parameters(rate_units));
163 }
164 }
165
setParameters(const AnyMap & node,const Kinetics & kin)166 void Reaction::setParameters(const AnyMap& node, const Kinetics& kin)
167 {
168 if (node.empty()) {
169 // empty node: used by newReaction() factory loader
170 return;
171 }
172
173 parseReactionEquation(*this, node["equation"], kin);
174 // Non-stoichiometric reaction orders
175 if (node.hasKey("orders")) {
176 for (const auto& order : node["orders"].asMap<double>()) {
177 orders[order.first] = order.second;
178 if (kin.kineticsSpeciesIndex(order.first) == npos) {
179 setValid(false);
180 }
181 }
182 }
183
184 // remove optional third body notation (used by ChebyshevReaction)
185 reactants.erase("(+M)");
186 products.erase("(+M)");
187
188 // Flags
189 id = node.getString("id", "");
190 duplicate = node.getBool("duplicate", false);
191 allow_negative_orders = node.getBool("negative-orders", false);
192 allow_nonreactant_orders = node.getBool("nonreactant-orders", false);
193
194 calculateRateCoeffUnits(kin);
195 input = node;
196 }
197
setRate(shared_ptr<ReactionRateBase> rate)198 void Reaction::setRate(shared_ptr<ReactionRateBase> rate)
199 {
200 if (!rate) {
201 // null pointer
202 m_rate.reset();
203 } else {
204 m_rate = rate;
205 }
206 }
207
reactantString() const208 std::string Reaction::reactantString() const
209 {
210 std::ostringstream result;
211 for (auto iter = reactants.begin(); iter != reactants.end(); ++iter) {
212 if (iter != reactants.begin()) {
213 result << " + ";
214 }
215 if (iter->second != 1.0) {
216 result << iter->second << " ";
217 }
218 result << iter->first;
219 }
220 return result.str();
221 }
222
productString() const223 std::string Reaction::productString() const
224 {
225 std::ostringstream result;
226 for (auto iter = products.begin(); iter != products.end(); ++iter) {
227 if (iter != products.begin()) {
228 result << " + ";
229 }
230 if (iter->second != 1.0) {
231 result << iter->second << " ";
232 }
233 result << iter->first;
234 }
235 return result.str();
236 }
237
equation() const238 std::string Reaction::equation() const
239 {
240 if (reversible) {
241 return reactantString() + " <=> " + productString();
242 } else {
243 return reactantString() + " => " + productString();
244 }
245 }
246
calculateRateCoeffUnits(const Kinetics & kin)247 void Reaction::calculateRateCoeffUnits(const Kinetics& kin)
248 {
249 if (!valid()) {
250 // If a reaction is invalid because of missing species in the Kinetics
251 // object, determining the units of the rate coefficient is impossible.
252 return;
253 }
254
255 // Determine the units of the rate coefficient
256 Units rxn_phase_units = kin.thermo(kin.reactionPhaseIndex()).standardConcentrationUnits();
257 rate_units = rxn_phase_units;
258 rate_units *= Units(1.0, 0, 0, -1);
259 for (const auto& order : orders) {
260 const auto& phase = kin.speciesPhase(order.first);
261 rate_units *= phase.standardConcentrationUnits().pow(-order.second);
262 }
263 for (const auto& stoich : reactants) {
264 // Order for each reactant is the reactant stoichiometric coefficient,
265 // unless already overridden by user-specified orders
266 if (stoich.first == "M" || ba::starts_with(stoich.first, "(+")) {
267 // calculateRateCoeffUnits may be called before these pseudo-species
268 // have been stripped from the reactants
269 continue;
270 } else if (orders.find(stoich.first) == orders.end()) {
271 const auto& phase = kin.speciesPhase(stoich.first);
272 rate_units *= phase.standardConcentrationUnits().pow(-stoich.second);
273 }
274 }
275
276 if (m_rate) {
277 // Ensure that reaction rate object is updated
278 m_rate->setUnits(rate_units);
279 }
280 }
281
updateUndeclared(std::vector<std::string> & undeclared,const Composition & comp,const Kinetics & kin)282 void updateUndeclared(std::vector<std::string>& undeclared,
283 const Composition& comp, const Kinetics& kin)
284 {
285 for (const auto& sp: comp) {
286 if (kin.kineticsSpeciesIndex(sp.first) == npos) {
287 undeclared.emplace_back(sp.first);
288 }
289 }
290 }
291
undeclaredThirdBodies(const Kinetics & kin) const292 std::pair<std::vector<std::string>, bool> Reaction::undeclaredThirdBodies(
293 const Kinetics& kin) const
294 {
295 std::vector<std::string> undeclared;
296 if (m_third_body) {
297 updateUndeclared(undeclared, m_third_body->efficiencies, kin);
298 bool specified_collision_partner = dynamic_cast<const ThreeBodyReaction3*>(
299 this)->specified_collision_partner;
300 return std::make_pair(undeclared, specified_collision_partner);
301 }
302 return std::make_pair(undeclared, false);
303 }
304
checkBalance(const Kinetics & kin) const305 void Reaction::checkBalance(const Kinetics& kin) const
306 {
307 Composition balr, balp;
308
309 // iterate over products and reactants
310 for (const auto& sp : products) {
311 const ThermoPhase& ph = kin.speciesPhase(sp.first);
312 size_t k = ph.speciesIndex(sp.first);
313 double stoich = sp.second;
314 for (size_t m = 0; m < ph.nElements(); m++) {
315 balr[ph.elementName(m)] = 0.0; // so that balr contains all species
316 balp[ph.elementName(m)] += stoich * ph.nAtoms(k, m);
317 }
318 }
319 for (const auto& sp : reactants) {
320 const ThermoPhase& ph = kin.speciesPhase(sp.first);
321 size_t k = ph.speciesIndex(sp.first);
322 double stoich = sp.second;
323 for (size_t m = 0; m < ph.nElements(); m++) {
324 balr[ph.elementName(m)] += stoich * ph.nAtoms(k, m);
325 }
326 }
327
328 std::string msg;
329 bool ok = true;
330 for (const auto& el : balr) {
331 const std::string& elem = el.first;
332 double elemsum = balr[elem] + balp[elem];
333 double elemdiff = fabs(balp[elem] - balr[elem]);
334 if (elemsum > 0.0 && elemdiff / elemsum > 1e-4) {
335 ok = false;
336 msg += fmt::format(" {} {} {}\n",
337 elem, balr[elem], balp[elem]);
338 }
339 }
340 if (!ok) {
341 throw InputFileError("Reaction::checkBalance", input,
342 "The following reaction is unbalanced: {}\n"
343 " Element Reactants Products\n{}",
344 equation(), msg);
345 }
346 }
347
checkSpecies(const Kinetics & kin) const348 bool Reaction::checkSpecies(const Kinetics& kin) const
349 {
350 // Check for undeclared species
351 std::vector<std::string> undeclared;
352 updateUndeclared(undeclared, reactants, kin);
353 updateUndeclared(undeclared, products, kin);
354 if (!undeclared.empty()) {
355 if (kin.skipUndeclaredSpecies()) {
356 return false;
357 } else {
358 throw InputFileError("Reaction::checkSpecies", input, "Reaction '{}'\n"
359 "contains undeclared species: '{}'",
360 equation(), boost::algorithm::join(undeclared, "', '"));
361 }
362 }
363
364 undeclared.clear();
365 updateUndeclared(undeclared, orders, kin);
366 if (!undeclared.empty()) {
367 if (kin.skipUndeclaredSpecies()) {
368 return false;
369 } else {
370 if (input.hasKey("orders")) {
371 throw InputFileError("Reaction::checkSpecies", input["orders"],
372 "Reaction '{}'\n"
373 "defines reaction orders for undeclared species: '{}'",
374 equation(), boost::algorithm::join(undeclared, "', '"));
375 }
376 // Error for empty input AnyMap (e.g. XML)
377 throw InputFileError("Reaction::checkSpecies", input, "Reaction '{}'\n"
378 "defines reaction orders for undeclared species: '{}'",
379 equation(), boost::algorithm::join(undeclared, "', '"));
380 }
381 }
382
383 // Use helper function while there is no uniform handling of third bodies
384 auto third = undeclaredThirdBodies(kin);
385 undeclared = third.first;
386 bool specified_collision_partner_ = third.second;
387 if (!undeclared.empty()) {
388 if (!kin.skipUndeclaredThirdBodies()) {
389 if (input.hasKey("efficiencies")) {
390 throw InputFileError("Reaction::checkSpecies", input["efficiencies"],
391 "Reaction '{}'\n"
392 "defines third-body efficiencies for undeclared species: '{}'",
393 equation(), boost::algorithm::join(undeclared, "', '"));
394 }
395 // Error for specified ThirdBody or empty input AnyMap
396 throw InputFileError("Reaction::checkSpecies", input, "Reaction '{}'\n"
397 "is a three-body reaction with undeclared species: '{}'",
398 equation(), boost::algorithm::join(undeclared, "', '"));
399 } else if (kin.skipUndeclaredSpecies() && specified_collision_partner_) {
400 return false;
401 }
402 }
403
404 checkBalance(kin);
405
406 return true;
407 }
408
ElementaryReaction2(const Composition & reactants_,const Composition products_,const Arrhenius & rate_)409 ElementaryReaction2::ElementaryReaction2(const Composition& reactants_,
410 const Composition products_,
411 const Arrhenius& rate_)
412 : Reaction(reactants_, products_)
413 , rate(rate_)
414 , allow_negative_pre_exponential_factor(false)
415 {
416 reaction_type = ELEMENTARY_RXN;
417 }
418
ElementaryReaction2()419 ElementaryReaction2::ElementaryReaction2()
420 : Reaction()
421 , allow_negative_pre_exponential_factor(false)
422 {
423 reaction_type = ELEMENTARY_RXN;
424 }
425
validate()426 void ElementaryReaction2::validate()
427 {
428 Reaction::validate();
429 if (!allow_negative_pre_exponential_factor &&
430 rate.preExponentialFactor() < 0) {
431 throw InputFileError("ElementaryReaction2::validate", input,
432 "Undeclared negative pre-exponential factor found in reaction '"
433 + equation() + "'");
434 }
435 }
436
getParameters(AnyMap & reactionNode) const437 void ElementaryReaction2::getParameters(AnyMap& reactionNode) const
438 {
439 Reaction::getParameters(reactionNode);
440 if (allow_negative_pre_exponential_factor) {
441 reactionNode["negative-A"] = true;
442 }
443 AnyMap rateNode;
444 rate.getParameters(rateNode, rate_units);
445 reactionNode["rate-constant"] = std::move(rateNode);
446 }
447
ThirdBody(double default_eff)448 ThirdBody::ThirdBody(double default_eff)
449 : default_efficiency(default_eff)
450 {
451 }
452
ThirdBody(const AnyMap & node)453 ThirdBody::ThirdBody(const AnyMap& node)
454 {
455 setEfficiencies(node);
456 }
457
setEfficiencies(const AnyMap & node)458 void ThirdBody::setEfficiencies(const AnyMap& node)
459 {
460 default_efficiency = node.getDouble("default-efficiency", 1.0);
461 if (node.hasKey("efficiencies")) {
462 efficiencies = node["efficiencies"].asMap<double>();
463 }
464 }
465
efficiency(const std::string & k) const466 double ThirdBody::efficiency(const std::string& k) const
467 {
468 return getValue(efficiencies, k, default_efficiency);
469 }
470
ThreeBodyReaction2()471 ThreeBodyReaction2::ThreeBodyReaction2()
472 {
473 reaction_type = THREE_BODY_RXN;
474 }
475
ThreeBodyReaction2(const Composition & reactants_,const Composition & products_,const Arrhenius & rate_,const ThirdBody & tbody)476 ThreeBodyReaction2::ThreeBodyReaction2(const Composition& reactants_,
477 const Composition& products_,
478 const Arrhenius& rate_,
479 const ThirdBody& tbody)
480 : ElementaryReaction2(reactants_, products_, rate_)
481 , third_body(tbody)
482 {
483 reaction_type = THREE_BODY_RXN;
484 }
485
reactantString() const486 std::string ThreeBodyReaction2::reactantString() const
487 {
488 if (specified_collision_partner) {
489 return ElementaryReaction2::reactantString() + " + "
490 + third_body.efficiencies.begin()->first;
491 } else {
492 return ElementaryReaction2::reactantString() + " + M";
493 }
494 }
495
productString() const496 std::string ThreeBodyReaction2::productString() const
497 {
498 if (specified_collision_partner) {
499 return ElementaryReaction2::productString() + " + "
500 + third_body.efficiencies.begin()->first;
501 } else {
502 return ElementaryReaction2::productString() + " + M";
503 }
504 }
505
calculateRateCoeffUnits(const Kinetics & kin)506 void ThreeBodyReaction2::calculateRateCoeffUnits(const Kinetics& kin)
507 {
508 ElementaryReaction2::calculateRateCoeffUnits(kin);
509 bool specified_collision_partner_ = false;
510 for (const auto& reac : reactants) {
511 // While this reaction was already identified as a three-body reaction in a
512 // pre-processing step, this method is often called before a three-body
513 // reaction is fully instantiated. For the determination of the correct units,
514 // it is necessary to check whether the reaction uses a generic 'M' or an
515 // explicitly specified collision partner that may not have been deleted yet.
516 if (reac.first != "M" && products.count(reac.first)) {
517 // detected specified third-body collision partner
518 specified_collision_partner_ = true;
519 }
520 }
521 if (!specified_collision_partner_) {
522 const ThermoPhase& rxn_phase = kin.thermo(kin.reactionPhaseIndex());
523 rate_units *= rxn_phase.standardConcentrationUnits().pow(-1);
524 }
525 }
526
getParameters(AnyMap & reactionNode) const527 void ThreeBodyReaction2::getParameters(AnyMap& reactionNode) const
528 {
529 ElementaryReaction2::getParameters(reactionNode);
530 if (!specified_collision_partner) {
531 reactionNode["type"] = "three-body";
532 reactionNode["efficiencies"] = third_body.efficiencies;
533 reactionNode["efficiencies"].setFlowStyle();
534 if (third_body.default_efficiency != 1.0) {
535 reactionNode["default-efficiency"] = third_body.default_efficiency;
536 }
537 }
538 }
539
undeclaredThirdBodies(const Kinetics & kin) const540 std::pair<std::vector<std::string>, bool> ThreeBodyReaction2::undeclaredThirdBodies(
541 const Kinetics& kin) const
542 {
543 std::vector<std::string> undeclared;
544 updateUndeclared(undeclared, third_body.efficiencies, kin);
545 return std::make_pair(undeclared, specified_collision_partner);
546 }
547
FalloffReaction()548 FalloffReaction::FalloffReaction()
549 : Reaction()
550 , falloff(new Falloff())
551 , allow_negative_pre_exponential_factor(false)
552 , low_rate_units(0.0)
553 {
554 reaction_type = FALLOFF_RXN;
555 }
556
FalloffReaction(const Composition & reactants_,const Composition & products_,const Arrhenius & low_rate_,const Arrhenius & high_rate_,const ThirdBody & tbody)557 FalloffReaction::FalloffReaction(
558 const Composition& reactants_, const Composition& products_,
559 const Arrhenius& low_rate_, const Arrhenius& high_rate_,
560 const ThirdBody& tbody)
561 : Reaction(reactants_, products_)
562 , low_rate(low_rate_)
563 , high_rate(high_rate_)
564 , third_body(tbody)
565 , falloff(new Falloff())
566 , allow_negative_pre_exponential_factor(false)
567 , low_rate_units(0.0)
568 {
569 reaction_type = FALLOFF_RXN;
570 }
571
reactantString() const572 std::string FalloffReaction::reactantString() const
573 {
574 if (third_body.default_efficiency == 0 &&
575 third_body.efficiencies.size() == 1) {
576 return Reaction::reactantString() + " (+" +
577 third_body.efficiencies.begin()->first + ")";
578 } else {
579 return Reaction::reactantString() + " (+M)";
580 }
581 }
582
productString() const583 std::string FalloffReaction::productString() const
584 {
585 if (third_body.default_efficiency == 0 &&
586 third_body.efficiencies.size() == 1) {
587 return Reaction::productString() + " (+" +
588 third_body.efficiencies.begin()->first + ")";
589 } else {
590 return Reaction::productString() + " (+M)";
591 }
592 }
593
validate()594 void FalloffReaction::validate()
595 {
596 Reaction::validate();
597 if (!allow_negative_pre_exponential_factor &&
598 (low_rate.preExponentialFactor() < 0 ||
599 high_rate.preExponentialFactor() < 0)) {
600 throw InputFileError("FalloffReaction::validate", input, "Negative "
601 "pre-exponential factor found for reaction '" + equation() + "'");
602 }
603 if (low_rate.preExponentialFactor() * high_rate.preExponentialFactor() < 0) {
604 throw InputFileError("FalloffReaction::validate", input, "High and "
605 "low rate pre-exponential factors must have the same sign."
606 "Reaction: '{}'", equation());
607 }
608 }
609
calculateRateCoeffUnits(const Kinetics & kin)610 void FalloffReaction::calculateRateCoeffUnits(const Kinetics& kin)
611 {
612 Reaction::calculateRateCoeffUnits(kin);
613 const ThermoPhase& rxn_phase = kin.thermo(kin.reactionPhaseIndex());
614 low_rate_units = rate_units;
615 low_rate_units *= rxn_phase.standardConcentrationUnits().pow(-1);
616 }
617
getParameters(AnyMap & reactionNode) const618 void FalloffReaction::getParameters(AnyMap& reactionNode) const
619 {
620 Reaction::getParameters(reactionNode);
621 reactionNode["type"] = "falloff";
622 AnyMap lowRateNode;
623 low_rate.getParameters(lowRateNode, low_rate_units);
624 reactionNode["low-P-rate-constant"] = std::move(lowRateNode);
625 AnyMap highRateNode;
626 high_rate.getParameters(highRateNode, rate_units);
627 reactionNode["high-P-rate-constant"] = std::move(highRateNode);
628 falloff->getParameters(reactionNode);
629
630 reactionNode["efficiencies"] = third_body.efficiencies;
631 reactionNode["efficiencies"].setFlowStyle();
632 if (third_body.default_efficiency != 1.0) {
633 reactionNode["default-efficiency"] = third_body.default_efficiency;
634 }
635 }
636
undeclaredThirdBodies(const Kinetics & kin) const637 std::pair<std::vector<std::string>, bool> FalloffReaction::undeclaredThirdBodies(
638 const Kinetics& kin) const
639 {
640 std::vector<std::string> undeclared;
641 updateUndeclared(undeclared, third_body.efficiencies, kin);
642 return std::make_pair(undeclared, false);
643 }
644
ChemicallyActivatedReaction()645 ChemicallyActivatedReaction::ChemicallyActivatedReaction()
646 {
647 reaction_type = CHEMACT_RXN;
648 }
649
ChemicallyActivatedReaction(const Composition & reactants_,const Composition & products_,const Arrhenius & low_rate_,const Arrhenius & high_rate_,const ThirdBody & tbody)650 ChemicallyActivatedReaction::ChemicallyActivatedReaction(
651 const Composition& reactants_, const Composition& products_,
652 const Arrhenius& low_rate_, const Arrhenius& high_rate_,
653 const ThirdBody& tbody)
654 : FalloffReaction(reactants_, products_, low_rate_, high_rate_, tbody)
655 {
656 reaction_type = CHEMACT_RXN;
657 }
658
calculateRateCoeffUnits(const Kinetics & kin)659 void ChemicallyActivatedReaction::calculateRateCoeffUnits(const Kinetics& kin)
660 {
661 Reaction::calculateRateCoeffUnits(kin); // Skip FalloffReaction
662 const ThermoPhase& rxn_phase = kin.thermo(kin.reactionPhaseIndex());
663 low_rate_units = rate_units;
664 rate_units *= rxn_phase.standardConcentrationUnits();
665 }
666
getParameters(AnyMap & reactionNode) const667 void ChemicallyActivatedReaction::getParameters(AnyMap& reactionNode) const
668 {
669 FalloffReaction::getParameters(reactionNode);
670 reactionNode["type"] = "chemically-activated";
671 }
672
PlogReaction2()673 PlogReaction2::PlogReaction2()
674 : Reaction()
675 {
676 reaction_type = PLOG_RXN;
677 }
678
PlogReaction2(const Composition & reactants_,const Composition & products_,const Plog & rate_)679 PlogReaction2::PlogReaction2(const Composition& reactants_,
680 const Composition& products_, const Plog& rate_)
681 : Reaction(reactants_, products_)
682 , rate(rate_)
683 {
684 reaction_type = PLOG_RXN;
685 }
686
getParameters(AnyMap & reactionNode) const687 void PlogReaction2::getParameters(AnyMap& reactionNode) const
688 {
689 Reaction::getParameters(reactionNode);
690 reactionNode["type"] = "pressure-dependent-Arrhenius";
691 rate.getParameters(reactionNode, rate_units);
692 }
693
ChebyshevReaction2()694 ChebyshevReaction2::ChebyshevReaction2()
695 : Reaction()
696 {
697 reaction_type = CHEBYSHEV_RXN;
698 }
699
ChebyshevReaction2(const Composition & reactants_,const Composition & products_,const Chebyshev & rate_)700 ChebyshevReaction2::ChebyshevReaction2(const Composition& reactants_,
701 const Composition& products_,
702 const Chebyshev& rate_)
703 : Reaction(reactants_, products_)
704 , rate(rate_)
705 {
706 reaction_type = CHEBYSHEV_RXN;
707 }
708
getParameters(AnyMap & reactionNode) const709 void ChebyshevReaction2::getParameters(AnyMap& reactionNode) const
710 {
711 Reaction::getParameters(reactionNode);
712 reactionNode["type"] = "Chebyshev";
713 rate.getParameters(reactionNode, rate_units);
714 }
715
InterfaceReaction()716 InterfaceReaction::InterfaceReaction()
717 : is_sticking_coefficient(false)
718 , use_motz_wise_correction(false)
719 {
720 reaction_type = INTERFACE_RXN;
721 }
722
InterfaceReaction(const Composition & reactants_,const Composition & products_,const Arrhenius & rate_,bool isStick)723 InterfaceReaction::InterfaceReaction(const Composition& reactants_,
724 const Composition& products_,
725 const Arrhenius& rate_,
726 bool isStick)
727 : ElementaryReaction2(reactants_, products_, rate_)
728 , is_sticking_coefficient(isStick)
729 , use_motz_wise_correction(false)
730 {
731 reaction_type = INTERFACE_RXN;
732 }
733
calculateRateCoeffUnits(const Kinetics & kin)734 void InterfaceReaction::calculateRateCoeffUnits(const Kinetics& kin)
735 {
736 ElementaryReaction2::calculateRateCoeffUnits(kin);
737 if (is_sticking_coefficient || input.hasKey("sticking-coefficient")) {
738 rate_units = Units(1.0); // sticking coefficients are dimensionless
739 }
740 }
741
getParameters(AnyMap & reactionNode) const742 void InterfaceReaction::getParameters(AnyMap& reactionNode) const
743 {
744 ElementaryReaction2::getParameters(reactionNode);
745 if (is_sticking_coefficient) {
746 reactionNode["sticking-coefficient"] = std::move(reactionNode["rate-constant"]);
747 reactionNode.erase("rate-constant");
748 }
749 if (use_motz_wise_correction) {
750 reactionNode["Motz-Wise"] = true;
751 }
752 if (!sticking_species.empty()) {
753 reactionNode["sticking-species"] = sticking_species;
754 }
755 if (!coverage_deps.empty()) {
756 AnyMap deps;
757 for (const auto& d : coverage_deps) {
758 AnyMap dep;
759 dep["a"] = d.second.a;
760 dep["m"] = d.second.m;
761 dep["E"].setQuantity(d.second.E, "K", true);
762 deps[d.first] = std::move(dep);
763 }
764 reactionNode["coverage-dependencies"] = std::move(deps);
765 }
766 }
767
ElectrochemicalReaction()768 ElectrochemicalReaction::ElectrochemicalReaction()
769 : beta(0.5)
770 , exchange_current_density_formulation(false)
771 {
772 }
773
ElectrochemicalReaction(const Composition & reactants_,const Composition & products_,const Arrhenius & rate_)774 ElectrochemicalReaction::ElectrochemicalReaction(const Composition& reactants_,
775 const Composition& products_,
776 const Arrhenius& rate_)
777 : InterfaceReaction(reactants_, products_, rate_)
778 , beta(0.5)
779 , exchange_current_density_formulation(false)
780 {
781 }
782
getParameters(AnyMap & reactionNode) const783 void ElectrochemicalReaction::getParameters(AnyMap& reactionNode) const
784 {
785 InterfaceReaction::getParameters(reactionNode);
786 if (beta != 0.5) {
787 reactionNode["beta"] = beta;
788 }
789 if (exchange_current_density_formulation) {
790 reactionNode["exchange-current-density-formulation"] = true;
791 }
792 }
793
BlowersMaselReaction()794 BlowersMaselReaction::BlowersMaselReaction()
795 : allow_negative_pre_exponential_factor(false)
796 {
797 reaction_type = BLOWERSMASEL_RXN;
798 }
799
validate()800 void BlowersMaselReaction::validate()
801 {
802 Reaction::validate();
803 if (!allow_negative_pre_exponential_factor &&
804 rate.preExponentialFactor() < 0) {
805 throw InputFileError("BlowersMaselReaction::validate", input,
806 "Undeclared negative pre-exponential factor found in reaction '"
807 + equation() + "'");
808 }
809 }
810
BlowersMaselReaction(const Composition & reactants_,const Composition & products_,const BlowersMasel & rate_)811 BlowersMaselReaction::BlowersMaselReaction(const Composition& reactants_,
812 const Composition& products_,
813 const BlowersMasel& rate_)
814 : Reaction(reactants_, products_)
815 , rate(rate_)
816 , allow_negative_pre_exponential_factor(false)
817 {
818 reaction_type = BLOWERSMASEL_RXN;
819 }
820
getParameters(AnyMap & reactionNode) const821 void BlowersMaselReaction::getParameters(AnyMap& reactionNode) const
822 {
823 Reaction::getParameters(reactionNode);
824 reactionNode["type"] = "Blowers-Masel";
825 if (allow_negative_pre_exponential_factor) {
826 reactionNode["negative-A"] = true;
827 }
828 AnyMap rateNode;
829 rate.getParameters(rateNode, rate_units);
830 reactionNode["rate-constant"] = std::move(rateNode);
831 }
832
BlowersMaselInterfaceReaction()833 BlowersMaselInterfaceReaction::BlowersMaselInterfaceReaction()
834 : is_sticking_coefficient(false)
835 , use_motz_wise_correction(false)
836 {
837 reaction_type = BMINTERFACE_RXN;
838 }
839
BlowersMaselInterfaceReaction(const Composition & reactants_,const Composition & products_,const BlowersMasel & rate_,bool isStick)840 BlowersMaselInterfaceReaction::BlowersMaselInterfaceReaction(const Composition& reactants_,
841 const Composition& products_,
842 const BlowersMasel& rate_,
843 bool isStick)
844 : BlowersMaselReaction(reactants_, products_, rate_)
845 , is_sticking_coefficient(isStick)
846 , use_motz_wise_correction(false)
847 {
848 reaction_type = BMINTERFACE_RXN;
849 }
850
calculateRateCoeffUnits(const Kinetics & kin)851 void BlowersMaselInterfaceReaction::calculateRateCoeffUnits(const Kinetics& kin)
852 {
853 BlowersMaselReaction::calculateRateCoeffUnits(kin);
854 if (is_sticking_coefficient || input.hasKey("sticking-coefficient")) {
855 rate_units = Units(1.0); // sticking coefficients are dimensionless
856 }
857 }
858
getParameters(AnyMap & reactionNode) const859 void BlowersMaselInterfaceReaction::getParameters(AnyMap& reactionNode) const
860 {
861 BlowersMaselReaction::getParameters(reactionNode);
862 reactionNode["type"] = "Blowers-Masel";
863 if (is_sticking_coefficient) {
864 reactionNode["sticking-coefficient"] = std::move(reactionNode["rate-constant"]);
865 reactionNode.erase("rate-constant");
866 }
867 if (use_motz_wise_correction) {
868 reactionNode["Motz-Wise"] = true;
869 }
870 if (!sticking_species.empty()) {
871 reactionNode["sticking-species"] = sticking_species;
872 }
873 if (!coverage_deps.empty()) {
874 AnyMap deps;
875 for (const auto& d : coverage_deps) {
876 AnyMap dep;
877 dep["a"] = d.second.a;
878 dep["m"] = d.second.m;
879 dep["E"].setQuantity(d.second.E, "K", true);
880 deps[d.first] = std::move(dep);
881 }
882 reactionNode["coverage-dependencies"] = std::move(deps);
883 }
884 }
885
ElementaryReaction3()886 ElementaryReaction3::ElementaryReaction3()
887 {
888 setRate(newReactionRate(type()));
889 }
890
ElementaryReaction3(const Composition & reactants,const Composition & products,const ArrheniusRate & rate)891 ElementaryReaction3::ElementaryReaction3(const Composition& reactants,
892 const Composition& products,
893 const ArrheniusRate& rate)
894 : Reaction(reactants, products)
895 {
896 m_rate.reset(new ArrheniusRate(rate));
897 }
898
ElementaryReaction3(const AnyMap & node,const Kinetics & kin)899 ElementaryReaction3::ElementaryReaction3(const AnyMap& node, const Kinetics& kin)
900 {
901 if (!node.empty()) {
902 setParameters(node, kin);
903 setRate(newReactionRate(node, rate_units));
904 } else {
905 setRate(newReactionRate(type()));
906 }
907 }
908
ThreeBodyReaction3()909 ThreeBodyReaction3::ThreeBodyReaction3()
910 {
911 m_third_body.reset(new ThirdBody);
912 setRate(newReactionRate(type()));
913 }
914
ThreeBodyReaction3(const Composition & reactants,const Composition & products,const ArrheniusRate & rate,const ThirdBody & tbody)915 ThreeBodyReaction3::ThreeBodyReaction3(const Composition& reactants,
916 const Composition& products,
917 const ArrheniusRate& rate,
918 const ThirdBody& tbody)
919 : ElementaryReaction3(reactants, products, rate)
920 {
921 m_third_body = std::make_shared<ThirdBody>(tbody);
922 }
923
ThreeBodyReaction3(const AnyMap & node,const Kinetics & kin)924 ThreeBodyReaction3::ThreeBodyReaction3(const AnyMap& node, const Kinetics& kin)
925 {
926 m_third_body.reset(new ThirdBody);
927 if (!node.empty()) {
928 setParameters(node, kin);
929 setRate(newReactionRate(node, rate_units));
930 } else {
931 setRate(newReactionRate(type()));
932 }
933 }
934
detectEfficiencies()935 bool ThreeBodyReaction3::detectEfficiencies()
936 {
937 for (const auto& reac : reactants) {
938 // detect explicitly specified collision partner
939 if (products.count(reac.first)) {
940 m_third_body->efficiencies[reac.first] = 1.;
941 }
942 }
943
944 if (m_third_body->efficiencies.size() == 0) {
945 return false;
946 } else if (m_third_body->efficiencies.size() > 1) {
947 throw CanteraError("ThreeBodyReaction3::detectEfficiencies",
948 "Found more than one explicitly specified collision partner\n"
949 "in reaction '{}'.", equation());
950 }
951
952 m_third_body->default_efficiency = 0.;
953 specified_collision_partner = true;
954 auto sp = m_third_body->efficiencies.begin();
955
956 // adjust reactant coefficients
957 auto reac = reactants.find(sp->first);
958 if (trunc(reac->second) != 1) {
959 reac->second -= 1.;
960 } else {
961 reactants.erase(reac);
962 }
963
964 // adjust product coefficients
965 auto prod = products.find(sp->first);
966 if (trunc(prod->second) != 1) {
967 prod->second -= 1.;
968 } else {
969 products.erase(prod);
970 }
971
972 return true;
973 }
974
calculateRateCoeffUnits(const Kinetics & kin)975 void ThreeBodyReaction3::calculateRateCoeffUnits(const Kinetics& kin)
976 {
977 ElementaryReaction3::calculateRateCoeffUnits(kin);
978 bool specified_collision_partner_ = false;
979 for (const auto& reac : reactants) {
980 // While this reaction was already identified as a three-body reaction in a
981 // pre-processing step, this method is often called before a three-body
982 // reaction is fully instantiated. For the determination of the correct units,
983 // it is necessary to check whether the reaction uses a generic 'M' or an
984 // explicitly specified collision partner that may not have been deleted yet.
985 if (reac.first != "M" && products.count(reac.first)) {
986 // detected specified third-body collision partner
987 specified_collision_partner_ = true;
988 }
989 }
990 if (!specified_collision_partner_) {
991 const ThermoPhase& rxn_phase = kin.thermo(kin.reactionPhaseIndex());
992 rate_units *= rxn_phase.standardConcentrationUnits().pow(-1);
993 }
994 }
995
setParameters(const AnyMap & node,const Kinetics & kin)996 void ThreeBodyReaction3::setParameters(const AnyMap& node, const Kinetics& kin)
997 {
998 if (node.empty()) {
999 // empty node: used by newReaction() factory loader
1000 return;
1001 }
1002 Reaction::setParameters(node, kin);
1003 if (reactants.count("M") != 1 || products.count("M") != 1) {
1004 if (!detectEfficiencies()) {
1005 throw InputFileError("ThreeBodyReaction3::setParameters", node["equation"],
1006 "Reaction equation '{}' does not contain third body 'M'",
1007 node["equation"].asString());
1008 }
1009 return;
1010 }
1011
1012 reactants.erase("M");
1013 products.erase("M");
1014 m_third_body->setEfficiencies(node);
1015 }
1016
getParameters(AnyMap & reactionNode) const1017 void ThreeBodyReaction3::getParameters(AnyMap& reactionNode) const
1018 {
1019 Reaction::getParameters(reactionNode);
1020 if (!specified_collision_partner) {
1021 reactionNode["type"] = "three-body";
1022 reactionNode["efficiencies"] = m_third_body->efficiencies;
1023 reactionNode["efficiencies"].setFlowStyle();
1024 if (m_third_body->default_efficiency != 1.0) {
1025 reactionNode["default-efficiency"] = m_third_body->default_efficiency;
1026 }
1027 }
1028 }
1029
reactantString() const1030 std::string ThreeBodyReaction3::reactantString() const
1031 {
1032 if (specified_collision_partner) {
1033 return ElementaryReaction3::reactantString() + " + "
1034 + m_third_body->efficiencies.begin()->first;
1035 } else {
1036 return ElementaryReaction3::reactantString() + " + M";
1037 }
1038 }
1039
productString() const1040 std::string ThreeBodyReaction3::productString() const
1041 {
1042 if (specified_collision_partner) {
1043 return ElementaryReaction3::productString() + " + "
1044 + m_third_body->efficiencies.begin()->first;
1045 } else {
1046 return ElementaryReaction3::productString() + " + M";
1047 }
1048 }
1049
PlogReaction3()1050 PlogReaction3::PlogReaction3()
1051 {
1052 setRate(newReactionRate(type()));
1053 }
1054
PlogReaction3(const Composition & reactants,const Composition & products,const PlogRate & rate)1055 PlogReaction3::PlogReaction3(const Composition& reactants,
1056 const Composition& products, const PlogRate& rate)
1057 : Reaction(reactants, products)
1058 {
1059 m_rate.reset(new PlogRate(rate));
1060 }
1061
PlogReaction3(const AnyMap & node,const Kinetics & kin)1062 PlogReaction3::PlogReaction3(const AnyMap& node, const Kinetics& kin)
1063 {
1064 if (!node.empty()) {
1065 setParameters(node, kin);
1066 setRate(newReactionRate(node, rate_units));
1067 } else {
1068 setRate(newReactionRate(type()));
1069 }
1070 }
1071
ChebyshevReaction3()1072 ChebyshevReaction3::ChebyshevReaction3()
1073 {
1074 setRate(newReactionRate(type()));
1075 }
1076
ChebyshevReaction3(const Composition & reactants,const Composition & products,const ChebyshevRate3 & rate)1077 ChebyshevReaction3::ChebyshevReaction3(const Composition& reactants,
1078 const Composition& products,
1079 const ChebyshevRate3& rate)
1080 : Reaction(reactants, products)
1081 {
1082 m_rate.reset(new ChebyshevRate3(rate));
1083 }
1084
ChebyshevReaction3(const AnyMap & node,const Kinetics & kin)1085 ChebyshevReaction3::ChebyshevReaction3(const AnyMap& node, const Kinetics& kin)
1086 {
1087 if (!node.empty()) {
1088 setParameters(node, kin);
1089 setRate(newReactionRate(node, rate_units));
1090 } else {
1091 setRate(newReactionRate(type()));
1092 }
1093 }
1094
CustomFunc1Reaction()1095 CustomFunc1Reaction::CustomFunc1Reaction()
1096 {
1097 setRate(newReactionRate(type()));
1098 }
1099
CustomFunc1Reaction(const Composition & reactants,const Composition & products,const CustomFunc1Rate & rate)1100 CustomFunc1Reaction::CustomFunc1Reaction(const Composition& reactants,
1101 const Composition& products,
1102 const CustomFunc1Rate& rate)
1103 : Reaction(reactants, products)
1104 {
1105 m_rate.reset(new CustomFunc1Rate(rate));
1106 }
1107
CustomFunc1Reaction(const AnyMap & node,const Kinetics & kin)1108 CustomFunc1Reaction::CustomFunc1Reaction(const AnyMap& node, const Kinetics& kin)
1109 {
1110 if (!node.empty()) {
1111 setParameters(node, kin);
1112 setRate(newReactionRate(node, rate_units));
1113 } else {
1114 setRate(newReactionRate(type()));
1115 }
1116 }
1117
readArrhenius(const XML_Node & arrhenius_node)1118 Arrhenius readArrhenius(const XML_Node& arrhenius_node)
1119 {
1120 return Arrhenius(getFloat(arrhenius_node, "A", "toSI"),
1121 getFloat(arrhenius_node, "b"),
1122 getFloat(arrhenius_node, "E", "actEnergy") / GasConstant);
1123 }
1124
readArrhenius(const Reaction & R,const AnyValue & rate,const Kinetics & kin,const UnitSystem & units,int pressure_dependence=0)1125 Arrhenius readArrhenius(const Reaction& R, const AnyValue& rate,
1126 const Kinetics& kin, const UnitSystem& units,
1127 int pressure_dependence=0)
1128 {
1129 double A, b, Ta;
1130 Units rc_units = R.rate_units;
1131 if (pressure_dependence) {
1132 Units rxn_phase_units = kin.thermo(kin.reactionPhaseIndex()).standardConcentrationUnits();
1133 rc_units *= rxn_phase_units.pow(-pressure_dependence);
1134 }
1135 if (rate.is<AnyMap>()) {
1136 auto& rate_map = rate.as<AnyMap>();
1137 A = units.convert(rate_map["A"], rc_units);
1138 b = rate_map["b"].asDouble();
1139 Ta = units.convertActivationEnergy(rate_map["Ea"], "K");
1140 } else {
1141 auto& rate_vec = rate.asVector<AnyValue>(3);
1142 A = units.convert(rate_vec[0], rc_units);
1143 b = rate_vec[1].asDouble();
1144 Ta = units.convertActivationEnergy(rate_vec[2], "K");
1145 }
1146 return Arrhenius(A, b, Ta);
1147 }
1148
1149 //! Parse falloff parameters, given a rateCoeff node
1150 /*!
1151 * @verbatim
1152 <falloff type="Troe"> 0.5 73.2 5000. 9999. </falloff>
1153 @endverbatim
1154 */
readFalloff(FalloffReaction & R,const XML_Node & rc_node)1155 void readFalloff(FalloffReaction& R, const XML_Node& rc_node)
1156 {
1157 XML_Node& falloff = rc_node.child("falloff");
1158 std::vector<std::string> p;
1159 vector_fp falloff_parameters;
1160 getStringArray(falloff, p);
1161 size_t np = p.size();
1162 for (size_t n = 0; n < np; n++) {
1163 falloff_parameters.push_back(fpValueCheck(p[n]));
1164 }
1165
1166 if (caseInsensitiveEquals(falloff["type"], "lindemann")) {
1167 if (np != 0) {
1168 throw CanteraError("readFalloff", "Lindemann parameterization "
1169 "takes no parameters, but {} were given", np);
1170 }
1171 R.falloff = newFalloff("Lindemann", falloff_parameters);
1172 } else if (caseInsensitiveEquals(falloff["type"], "troe")) {
1173 if (np != 3 && np != 4) {
1174 throw CanteraError("readFalloff", "Troe parameterization takes "
1175 "3 or 4 parameters, but {} were given", np);
1176 }
1177 R.falloff = newFalloff("Troe", falloff_parameters);
1178 } else if (caseInsensitiveEquals(falloff["type"], "sri")) {
1179 if (np != 3 && np != 5) {
1180 throw CanteraError("readFalloff", "SRI parameterization takes "
1181 "3 or 5 parameters, but {} were given", np);
1182 }
1183 R.falloff = newFalloff("SRI", falloff_parameters);
1184 } else if (caseInsensitiveEquals(falloff["type"], "tsang")) {
1185 if (np != 2) {
1186 throw CanteraError("readFalloff", "Tsang parameterization takes "
1187 "2 parameters, but {} were given", np);
1188 }
1189 R.falloff = newFalloff("Tsang", falloff_parameters);
1190 } else {
1191 throw CanteraError("readFalloff", "Unrecognized falloff type: '{}'",
1192 falloff["type"]);
1193 }
1194 }
1195
readFalloff(FalloffReaction & R,const AnyMap & node)1196 void readFalloff(FalloffReaction& R, const AnyMap& node)
1197 {
1198 if (node.hasKey("Troe")) {
1199 auto& f = node["Troe"].as<AnyMap>();
1200 vector_fp params{
1201 f["A"].asDouble(),
1202 f["T3"].asDouble(),
1203 f["T1"].asDouble()
1204 };
1205 if (f.hasKey("T2")) {
1206 params.push_back(f["T2"].asDouble());
1207 }
1208 R.falloff = newFalloff("Troe", params);
1209 } else if (node.hasKey("SRI")) {
1210 auto& f = node["SRI"].as<AnyMap>();
1211 vector_fp params{
1212 f["A"].asDouble(),
1213 f["B"].asDouble(),
1214 f["C"].asDouble()
1215 };
1216 if (f.hasKey("D")) {
1217 params.push_back(f["D"].asDouble());
1218 }
1219 if (f.hasKey("E")) {
1220 params.push_back(f["E"].asDouble());
1221 }
1222 R.falloff = newFalloff("SRI", params);
1223 } else if (node.hasKey("Tsang")) {
1224 auto& f = node["Tsang"].as<AnyMap>();
1225 vector_fp params{
1226 f["A"].asDouble(),
1227 f["B"].asDouble()
1228 };
1229 R.falloff = newFalloff("Tsang", params);
1230 } else {
1231 R.falloff = newFalloff("Lindemann", {});
1232 }
1233 }
1234
readEfficiencies(ThirdBody & tbody,const XML_Node & rc_node)1235 void readEfficiencies(ThirdBody& tbody, const XML_Node& rc_node)
1236 {
1237 if (!rc_node.hasChild("efficiencies")) {
1238 tbody.default_efficiency = 1.0;
1239 return;
1240 }
1241 const XML_Node& eff_node = rc_node.child("efficiencies");
1242 tbody.default_efficiency = fpValue(eff_node["default"]);
1243 tbody.efficiencies = parseCompString(eff_node.value());
1244 }
1245
readEfficiencies(ThirdBody & tbody,const AnyMap & node)1246 void readEfficiencies(ThirdBody& tbody, const AnyMap& node)
1247 {
1248 tbody.setEfficiencies(node);
1249 }
1250
readBlowersMasel(const Reaction & R,const AnyValue & rate,const Kinetics & kin,const UnitSystem & units,int pressure_dependence=0)1251 BlowersMasel readBlowersMasel(const Reaction& R, const AnyValue& rate,
1252 const Kinetics& kin, const UnitSystem& units,
1253 int pressure_dependence=0)
1254 {
1255 double A, b, Ta0, w;
1256 Units rc_units = R.rate_units;
1257 if (pressure_dependence) {
1258 Units rxn_phase_units = kin.thermo(kin.reactionPhaseIndex()).standardConcentrationUnits();
1259 rc_units *= rxn_phase_units.pow(-pressure_dependence);
1260 }
1261 if (rate.is<AnyMap>()) {
1262 auto& rate_map = rate.as<AnyMap>();
1263 A = units.convert(rate_map["A"], rc_units);
1264 b = rate_map["b"].asDouble();
1265 Ta0 = units.convertActivationEnergy(rate_map["Ea0"], "K");
1266 w = units.convertActivationEnergy(rate_map["w"], "K");
1267 } else {
1268 auto& rate_vec = rate.asVector<AnyValue>(4);
1269 A = units.convert(rate_vec[0], rc_units);
1270 b = rate_vec[1].asDouble();
1271 Ta0 = units.convertActivationEnergy(rate_vec[2], "K");
1272 w = units.convertActivationEnergy(rate_vec[3], "K");
1273 }
1274 return BlowersMasel(A, b, Ta0, w);
1275 }
1276
detectEfficiencies(ThreeBodyReaction2 & R)1277 bool detectEfficiencies(ThreeBodyReaction2& R)
1278 {
1279 for (const auto& reac : R.reactants) {
1280 // detect explicitly specified collision partner
1281 if (R.products.count(reac.first)) {
1282 R.third_body.efficiencies[reac.first] = 1.;
1283 }
1284 }
1285
1286 if (R.third_body.efficiencies.size() == 0) {
1287 return false;
1288 } else if (R.third_body.efficiencies.size() > 1) {
1289 throw CanteraError("detectEfficiencies",
1290 "Found more than one explicitly specified collision partner\n"
1291 "in reaction '{}'.", R.equation());
1292 }
1293
1294 R.third_body.default_efficiency = 0.;
1295 R.specified_collision_partner = true;
1296 auto sp = R.third_body.efficiencies.begin();
1297
1298 // adjust reactant coefficients
1299 auto reac = R.reactants.find(sp->first);
1300 if (trunc(reac->second) != 1) {
1301 reac->second -= 1.;
1302 } else {
1303 R.reactants.erase(reac);
1304 }
1305
1306 // adjust product coefficients
1307 auto prod = R.products.find(sp->first);
1308 if (trunc(prod->second) != 1) {
1309 prod->second -= 1.;
1310 } else {
1311 R.products.erase(prod);
1312 }
1313
1314 return true;
1315 }
1316
setupReaction(Reaction & R,const XML_Node & rxn_node)1317 void setupReaction(Reaction& R, const XML_Node& rxn_node)
1318 {
1319 // Reactant and product stoichiometries
1320 R.reactants = parseCompString(rxn_node.child("reactants").value());
1321 R.products = parseCompString(rxn_node.child("products").value());
1322
1323 // Non-stoichiometric reaction orders
1324 std::vector<XML_Node*> orders = rxn_node.getChildren("order");
1325 for (size_t i = 0; i < orders.size(); i++) {
1326 R.orders[orders[i]->attrib("species")] = orders[i]->fp_value();
1327 }
1328
1329 // Flags
1330 R.id = rxn_node.attrib("id");
1331 R.duplicate = rxn_node.hasAttrib("duplicate");
1332 const std::string& rev = rxn_node["reversible"];
1333 R.reversible = (rev == "true" || rev == "yes");
1334 }
1335
parseReactionEquation(Reaction & R,const AnyValue & equation,const Kinetics & kin)1336 void parseReactionEquation(Reaction& R, const AnyValue& equation,
1337 const Kinetics& kin) {
1338 // Parse the reaction equation to determine participating species and
1339 // stoichiometric coefficients
1340 std::vector<std::string> tokens;
1341 tokenizeString(equation.asString(), tokens);
1342 tokens.push_back("+"); // makes parsing last species not a special case
1343
1344 size_t last_used = npos; // index of last-used token
1345 bool reactants = true;
1346 for (size_t i = 1; i < tokens.size(); i++) {
1347 if (tokens[i] == "+" || ba::starts_with(tokens[i], "(+") ||
1348 tokens[i] == "<=>" || tokens[i] == "=" || tokens[i] == "=>") {
1349 std::string species = tokens[i-1];
1350
1351 double stoich;
1352 if (last_used != npos && tokens[last_used] == "(+") {
1353 // Falloff third body with space, e.g. "(+ M)"
1354 species = "(+" + species;
1355 stoich = -1;
1356 } else if (last_used == i-1 && ba::starts_with(species, "(+")
1357 && ba::ends_with(species, ")")) {
1358 // Falloff 3rd body written without space, e.g. "(+M)"
1359 stoich = -1;
1360 } else if (last_used == i-2) { // Species with no stoich. coefficient
1361 stoich = 1.0;
1362 } else if (last_used == i-3) { // Stoich. coefficient and species
1363 try {
1364 stoich = fpValueCheck(tokens[i-2]);
1365 } catch (CanteraError& err) {
1366 throw InputFileError("parseReactionEquation", equation,
1367 err.getMessage());
1368 }
1369 } else {
1370 throw InputFileError("parseReactionEquation", equation,
1371 "Error parsing reaction string '{}'.\n"
1372 "Current token: '{}'\nlast_used: '{}'",
1373 equation.asString(), tokens[i],
1374 (last_used == npos) ? "n/a" : tokens[last_used]);
1375 }
1376 if (kin.kineticsSpeciesIndex(species) == npos
1377 && stoich != -1 && species != "M") {
1378 R.setValid(false);
1379 }
1380
1381 if (reactants) {
1382 R.reactants[species] += stoich;
1383 } else {
1384 R.products[species] += stoich;
1385 }
1386
1387 last_used = i;
1388 }
1389
1390 // Tokens after this point are part of the products string
1391 if (tokens[i] == "<=>" || tokens[i] == "=") {
1392 R.reversible = true;
1393 reactants = false;
1394 } else if (tokens[i] == "=>") {
1395 R.reversible = false;
1396 reactants = false;
1397 }
1398 }
1399 }
1400
setupReaction(Reaction & R,const AnyMap & node,const Kinetics & kin)1401 void setupReaction(Reaction& R, const AnyMap& node, const Kinetics& kin)
1402 {
1403 parseReactionEquation(R, node["equation"], kin);
1404 // Non-stoichiometric reaction orders
1405 if (node.hasKey("orders")) {
1406 for (const auto& order : node["orders"].asMap<double>()) {
1407 R.orders[order.first] = order.second;
1408 if (kin.kineticsSpeciesIndex(order.first) == npos) {
1409 R.setValid(false);
1410 }
1411 }
1412 }
1413
1414 //Flags
1415 R.id = node.getString("id", "");
1416 R.duplicate = node.getBool("duplicate", false);
1417 R.allow_negative_orders = node.getBool("negative-orders", false);
1418 R.allow_nonreactant_orders = node.getBool("nonreactant-orders", false);
1419
1420 R.input = node;
1421 R.calculateRateCoeffUnits(kin);
1422 }
1423
setupElementaryReaction(ElementaryReaction2 & R,const XML_Node & rxn_node)1424 void setupElementaryReaction(ElementaryReaction2& R, const XML_Node& rxn_node)
1425 {
1426 const XML_Node& rc_node = rxn_node.child("rateCoeff");
1427 if (rc_node.hasChild("Arrhenius")) {
1428 R.rate = readArrhenius(rc_node.child("Arrhenius"));
1429 } else if (rc_node.hasChild("Arrhenius_ExchangeCurrentDensity")) {
1430 R.rate = readArrhenius(rc_node.child("Arrhenius_ExchangeCurrentDensity"));
1431 } else {
1432 throw CanteraError("setupElementaryReaction", "Couldn't find Arrhenius node");
1433 }
1434 if (rxn_node["negative_A"] == "yes") {
1435 R.allow_negative_pre_exponential_factor = true;
1436 }
1437 if (rxn_node["negative_orders"] == "yes") {
1438 R.allow_negative_orders = true;
1439 }
1440 if (rxn_node["nonreactant_orders"] == "yes") {
1441 R.allow_nonreactant_orders = true;
1442 }
1443 setupReaction(R, rxn_node);
1444 }
1445
setupElementaryReaction(ElementaryReaction2 & R,const AnyMap & node,const Kinetics & kin)1446 void setupElementaryReaction(ElementaryReaction2& R, const AnyMap& node,
1447 const Kinetics& kin)
1448 {
1449 setupReaction(R, node, kin);
1450 R.allow_negative_pre_exponential_factor = node.getBool("negative-A", false);
1451 R.rate = readArrhenius(R, node["rate-constant"], kin, node.units());
1452 }
1453
setupThreeBodyReaction(ThreeBodyReaction2 & R,const XML_Node & rxn_node)1454 void setupThreeBodyReaction(ThreeBodyReaction2& R, const XML_Node& rxn_node)
1455 {
1456 readEfficiencies(R.third_body, rxn_node.child("rateCoeff"));
1457 setupElementaryReaction(R, rxn_node);
1458 if (R.third_body.efficiencies.size() == 0) {
1459 detectEfficiencies(R);
1460 }
1461 }
1462
setupThreeBodyReaction(ThreeBodyReaction2 & R,const AnyMap & node,const Kinetics & kin)1463 void setupThreeBodyReaction(ThreeBodyReaction2& R, const AnyMap& node,
1464 const Kinetics& kin)
1465 {
1466 setupElementaryReaction(R, node, kin);
1467 if (R.reactants.count("M") != 1 || R.products.count("M") != 1) {
1468 if (!detectEfficiencies(R)) {
1469 throw InputFileError("setupThreeBodyReaction", node["equation"],
1470 "Reaction equation '{}' does not contain third body 'M'",
1471 node["equation"].asString());
1472 }
1473 } else {
1474 R.reactants.erase("M");
1475 R.products.erase("M");
1476 readEfficiencies(R.third_body, node);
1477 }
1478 }
1479
setupFalloffReaction(FalloffReaction & R,const XML_Node & rxn_node)1480 void setupFalloffReaction(FalloffReaction& R, const XML_Node& rxn_node)
1481 {
1482 XML_Node& rc_node = rxn_node.child("rateCoeff");
1483 std::vector<XML_Node*> rates = rc_node.getChildren("Arrhenius");
1484 int nLow = 0;
1485 int nHigh = 0;
1486 for (size_t i = 0; i < rates.size(); i++) {
1487 XML_Node& node = *rates[i];
1488 if (node["name"] == "") {
1489 R.high_rate = readArrhenius(node);
1490 nHigh++;
1491 } else if (node["name"] == "k0") {
1492 R.low_rate = readArrhenius(node);
1493 nLow++;
1494 } else {
1495 throw CanteraError("setupFalloffReaction", "Found an Arrhenius XML "
1496 "node with an unexpected type '" + node["name"] + "'");
1497 }
1498 }
1499 if (nLow != 1 || nHigh != 1) {
1500 throw CanteraError("setupFalloffReaction", "Did not find the correct "
1501 "number of Arrhenius rate expressions");
1502 }
1503 if (rxn_node["negative_A"] == "yes") {
1504 R.allow_negative_pre_exponential_factor = true;
1505 }
1506 readFalloff(R, rc_node);
1507 readEfficiencies(R.third_body, rc_node);
1508 setupReaction(R, rxn_node);
1509 }
1510
setupFalloffReaction(FalloffReaction & R,const AnyMap & node,const Kinetics & kin)1511 void setupFalloffReaction(FalloffReaction& R, const AnyMap& node,
1512 const Kinetics& kin)
1513 {
1514 setupReaction(R, node, kin);
1515 // setupReaction sets the stoichiometric coefficient for the falloff third
1516 // body to -1.
1517 std::string third_body;
1518 for (auto& reactant : R.reactants) {
1519 if (reactant.second == -1 && ba::starts_with(reactant.first, "(+")) {
1520 third_body = reactant.first;
1521 break;
1522 }
1523 }
1524
1525 // Equation must contain a third body, and it must appear on both sides
1526 if (third_body == "") {
1527 throw InputFileError("setupFalloffReaction", node["equation"],
1528 "Reactants for reaction '{}' do not contain a pressure-dependent "
1529 "third body", node["equation"].asString());
1530 } else if (R.products.count(third_body) == 0) {
1531 throw InputFileError("setupFalloffReaction", node["equation"],
1532 "Unable to match third body '{}' in reactants and products of "
1533 "reaction '{}'", third_body, node["equation"].asString());
1534 }
1535
1536 // Remove the dummy species
1537 R.reactants.erase(third_body);
1538 R.products.erase(third_body);
1539
1540 R.allow_negative_pre_exponential_factor = node.getBool("negative-A", false);
1541 if (third_body == "(+M)") {
1542 readEfficiencies(R.third_body, node);
1543 } else {
1544 // Specific species is listed as the third body
1545 R.third_body.default_efficiency = 0;
1546 R.third_body.efficiencies[third_body.substr(2, third_body.size() - 3)] = 1.0;
1547 }
1548
1549 R.low_rate = readArrhenius(R, node["low-P-rate-constant"], kin,
1550 node.units(), 1);
1551 R.high_rate = readArrhenius(R, node["high-P-rate-constant"], kin,
1552 node.units());
1553
1554 readFalloff(R, node);
1555 }
1556
setupChemicallyActivatedReaction(ChemicallyActivatedReaction & R,const XML_Node & rxn_node)1557 void setupChemicallyActivatedReaction(ChemicallyActivatedReaction& R,
1558 const XML_Node& rxn_node)
1559 {
1560 XML_Node& rc_node = rxn_node.child("rateCoeff");
1561 std::vector<XML_Node*> rates = rc_node.getChildren("Arrhenius");
1562 int nLow = 0;
1563 int nHigh = 0;
1564 for (size_t i = 0; i < rates.size(); i++) {
1565 XML_Node& node = *rates[i];
1566 if (node["name"] == "kHigh") {
1567 R.high_rate = readArrhenius(node);
1568 nHigh++;
1569 } else if (node["name"] == "") {
1570 R.low_rate = readArrhenius(node);
1571 nLow++;
1572 } else {
1573 throw CanteraError("setupChemicallyActivatedReaction", "Found an "
1574 "Arrhenius XML node with an unexpected type '" + node["name"] + "'");
1575 }
1576 }
1577 if (nLow != 1 || nHigh != 1) {
1578 throw CanteraError("setupChemicallyActivatedReaction", "Did not find "
1579 "the correct number of Arrhenius rate expressions");
1580 }
1581 readFalloff(R, rc_node);
1582 readEfficiencies(R.third_body, rc_node);
1583 setupReaction(R, rxn_node);
1584 }
1585
setupPlogReaction(PlogReaction2 & R,const XML_Node & rxn_node)1586 void setupPlogReaction(PlogReaction2& R, const XML_Node& rxn_node)
1587 {
1588 XML_Node& rc = rxn_node.child("rateCoeff");
1589 std::multimap<double, Arrhenius> rates;
1590 for (size_t m = 0; m < rc.nChildren(); m++) {
1591 const XML_Node& node = rc.child(m);
1592 rates.insert({getFloat(node, "P", "toSI"), readArrhenius(node)});
1593 }
1594 R.rate = Plog(rates);
1595 setupReaction(R, rxn_node);
1596 }
1597
setupPlogReaction(PlogReaction2 & R,const AnyMap & node,const Kinetics & kin)1598 void setupPlogReaction(PlogReaction2& R, const AnyMap& node, const Kinetics& kin)
1599 {
1600 setupReaction(R, node, kin);
1601 std::multimap<double, Arrhenius> rates;
1602 for (const auto& rate : node.at("rate-constants").asVector<AnyMap>()) {
1603 rates.insert({rate.convert("P", "Pa"),
1604 readArrhenius(R, AnyValue(rate), kin, node.units())});
1605 }
1606 R.rate = Plog(rates);
1607 }
1608
validate()1609 void PlogReaction2::validate()
1610 {
1611 Reaction::validate();
1612 rate.validate(equation());
1613 }
1614
setupChebyshevReaction(ChebyshevReaction2 & R,const XML_Node & rxn_node)1615 void setupChebyshevReaction(ChebyshevReaction2& R, const XML_Node& rxn_node)
1616 {
1617 XML_Node& rc = rxn_node.child("rateCoeff");
1618 const XML_Node& coeff_node = rc.child("floatArray");
1619 size_t nP = atoi(coeff_node["degreeP"].c_str());
1620 size_t nT = atoi(coeff_node["degreeT"].c_str());
1621
1622 vector_fp coeffs_flat;
1623 getFloatArray(rc, coeffs_flat, false);
1624 Array2D coeffs(nT, nP);
1625 for (size_t t = 0; t < nT; t++) {
1626 for (size_t p = 0; p < nP; p++) {
1627 coeffs(t,p) = coeffs_flat[nP*t + p];
1628 }
1629 }
1630 R.rate = Chebyshev(getFloat(rc, "Tmin", "toSI"),
1631 getFloat(rc, "Tmax", "toSI"),
1632 getFloat(rc, "Pmin", "toSI"),
1633 getFloat(rc, "Pmax", "toSI"),
1634 coeffs);
1635 setupReaction(R, rxn_node);
1636 }
1637
setupChebyshevReaction(ChebyshevReaction2 & R,const AnyMap & node,const Kinetics & kin)1638 void setupChebyshevReaction(ChebyshevReaction2&R, const AnyMap& node,
1639 const Kinetics& kin)
1640 {
1641 setupReaction(R, node, kin);
1642 R.reactants.erase("(+M)"); // remove optional third body notation
1643 R.products.erase("(+M)");
1644 const auto& T_range = node["temperature-range"].asVector<AnyValue>(2);
1645 const auto& P_range = node["pressure-range"].asVector<AnyValue>(2);
1646 auto& vcoeffs = node["data"].asVector<vector_fp>();
1647 Array2D coeffs(vcoeffs.size(), vcoeffs[0].size());
1648 for (size_t i = 0; i < coeffs.nRows(); i++) {
1649 if (vcoeffs[i].size() != vcoeffs[0].size()) {
1650 throw InputFileError("setupChebyshevReaction", node["data"],
1651 "Inconsistent number of coefficients in row {} of matrix", i + 1);
1652 }
1653 for (size_t j = 0; j < coeffs.nColumns(); j++) {
1654 coeffs(i, j) = vcoeffs[i][j];
1655 }
1656 }
1657 const UnitSystem& units = node.units();
1658 coeffs(0, 0) += std::log10(units.convertTo(1.0, R.rate_units));
1659 R.rate = Chebyshev(units.convert(T_range[0], "K"),
1660 units.convert(T_range[1], "K"),
1661 units.convert(P_range[0], "Pa"),
1662 units.convert(P_range[1], "Pa"),
1663 coeffs);
1664 }
1665
setupInterfaceReaction(InterfaceReaction & R,const XML_Node & rxn_node)1666 void setupInterfaceReaction(InterfaceReaction& R, const XML_Node& rxn_node)
1667 {
1668 XML_Node& arr = rxn_node.child("rateCoeff").child("Arrhenius");
1669 if (caseInsensitiveEquals(arr["type"], "stick")) {
1670 R.is_sticking_coefficient = true;
1671 R.sticking_species = arr["species"];
1672
1673 if (caseInsensitiveEquals(arr["motz_wise"], "true")) {
1674 R.use_motz_wise_correction = true;
1675 } else if (caseInsensitiveEquals(arr["motz_wise"], "false")) {
1676 R.use_motz_wise_correction = false;
1677 } else {
1678 // Default value for all reactions
1679 XML_Node* parent = rxn_node.parent();
1680 if (parent && parent->name() == "reactionData"
1681 && caseInsensitiveEquals((*parent)["motz_wise"], "true")) {
1682 R.use_motz_wise_correction = true;
1683 }
1684 }
1685 }
1686 std::vector<XML_Node*> cov = arr.getChildren("coverage");
1687 for (const auto& node : cov) {
1688 CoverageDependency& cdep = R.coverage_deps[node->attrib("species")];
1689 cdep.a = getFloat(*node, "a", "toSI");
1690 cdep.m = getFloat(*node, "m");
1691 cdep.E = getFloat(*node, "e", "actEnergy") / GasConstant;
1692 }
1693 setupElementaryReaction(R, rxn_node);
1694 }
1695
setupInterfaceReaction(InterfaceReaction & R,const AnyMap & node,const Kinetics & kin)1696 void setupInterfaceReaction(InterfaceReaction& R, const AnyMap& node,
1697 const Kinetics& kin)
1698 {
1699 setupReaction(R, node, kin);
1700 R.allow_negative_pre_exponential_factor = node.getBool("negative-A", false);
1701
1702 if (node.hasKey("rate-constant")) {
1703 R.rate = readArrhenius(R, node["rate-constant"], kin, node.units());
1704 } else if (node.hasKey("sticking-coefficient")) {
1705 R.is_sticking_coefficient = true;
1706 R.rate = readArrhenius(R, node["sticking-coefficient"], kin, node.units());
1707 R.use_motz_wise_correction = node.getBool("Motz-Wise",
1708 kin.thermo().input().getBool("Motz-Wise", false));
1709 R.sticking_species = node.getString("sticking-species", "");
1710 } else {
1711 throw InputFileError("setupInterfaceReaction", node,
1712 "Reaction must include either a 'rate-constant' or"
1713 " 'sticking-coefficient' node.");
1714 }
1715
1716 if (node.hasKey("coverage-dependencies")) {
1717 for (const auto& item : node["coverage-dependencies"].as<AnyMap>()) {
1718 double a, E, m;
1719 if (item.second.is<AnyMap>()) {
1720 auto& cov_map = item.second.as<AnyMap>();
1721 a = cov_map["a"].asDouble();
1722 m = cov_map["m"].asDouble();
1723 E = node.units().convertActivationEnergy(cov_map["E"], "K");
1724 } else {
1725 auto& cov_vec = item.second.asVector<AnyValue>(3);
1726 a = cov_vec[0].asDouble();
1727 m = cov_vec[1].asDouble();
1728 E = node.units().convertActivationEnergy(cov_vec[2], "K");
1729 }
1730 R.coverage_deps[item.first] = CoverageDependency(a, E, m);
1731 }
1732 }
1733 }
1734
setupElectrochemicalReaction(ElectrochemicalReaction & R,const XML_Node & rxn_node)1735 void setupElectrochemicalReaction(ElectrochemicalReaction& R,
1736 const XML_Node& rxn_node)
1737 {
1738 XML_Node& rc = rxn_node.child("rateCoeff");
1739 std::string rc_type = toLowerCopy(rc["type"]);
1740 if (rc_type == "exchangecurrentdensity") {
1741 R.exchange_current_density_formulation = true;
1742 } else if (rc_type != "" && rc_type != "arrhenius") {
1743 throw CanteraError("setupElectrochemicalReaction",
1744 "Unknown rate coefficient type: '" + rc_type + "'");
1745 }
1746 if (rc.hasChild("Arrhenius_ExchangeCurrentDensity")) {
1747 R.exchange_current_density_formulation = true;
1748 }
1749
1750 if (rc.hasChild("electrochem") && rc.child("electrochem").hasAttrib("beta")) {
1751 R.beta = fpValueCheck(rc.child("electrochem")["beta"]);
1752 }
1753
1754 setupInterfaceReaction(R, rxn_node);
1755
1756 // Override orders based on the <orders> node
1757 if (rxn_node.hasChild("orders")) {
1758 Composition orders = parseCompString(rxn_node.child("orders").value());
1759 for (const auto& order : orders) {
1760 R.orders[order.first] = order.second;
1761 }
1762 }
1763 }
1764
setupElectrochemicalReaction(ElectrochemicalReaction & R,const AnyMap & node,const Kinetics & kin)1765 void setupElectrochemicalReaction(ElectrochemicalReaction& R,
1766 const AnyMap& node, const Kinetics& kin)
1767 {
1768 setupInterfaceReaction(R, node, kin);
1769 R.beta = node.getDouble("beta", 0.5);
1770 R.exchange_current_density_formulation = node.getBool(
1771 "exchange-current-density-formulation", false);
1772 }
1773
setupBlowersMaselReaction(BlowersMaselReaction & R,const AnyMap & node,const Kinetics & kin)1774 void setupBlowersMaselReaction(BlowersMaselReaction& R, const AnyMap& node,
1775 const Kinetics& kin)
1776 {
1777 setupReaction(R, node, kin);
1778 R.allow_negative_pre_exponential_factor = node.getBool("negative-A", false);
1779 R.rate = readBlowersMasel(R, node["rate-constant"], kin, node.units());
1780 }
1781
setupBlowersMaselInterfaceReaction(BlowersMaselInterfaceReaction & R,const AnyMap & node,const Kinetics & kin)1782 void setupBlowersMaselInterfaceReaction(BlowersMaselInterfaceReaction& R, const AnyMap& node,
1783 const Kinetics& kin)
1784 {
1785 setupReaction(R, node, kin);
1786 R.allow_negative_pre_exponential_factor = node.getBool("negative-A", false);
1787
1788 if (node.hasKey("rate-constant")) {
1789 R.rate = readBlowersMasel(R, node["rate-constant"], kin, node.units());
1790 } else if (node.hasKey("sticking-coefficient")) {
1791 R.is_sticking_coefficient = true;
1792 R.rate = readBlowersMasel(R, node["sticking-coefficient"], kin, node.units());
1793 R.use_motz_wise_correction = node.getBool("Motz-Wise",
1794 kin.thermo().input().getBool("Motz-Wise", false));
1795 R.sticking_species = node.getString("sticking-species", "");
1796 } else {
1797 throw InputFileError("setupBlowersMaselInterfaceReaction", node,
1798 "Reaction must include either a 'rate-constant' or"
1799 " 'sticking-coefficient' node.");
1800 }
1801
1802 if (node.hasKey("coverage-dependencies")) {
1803 for (const auto& item : node["coverage-dependencies"].as<AnyMap>()) {
1804 double a, E, m;
1805 if (item.second.is<AnyMap>()) {
1806 auto& cov_map = item.second.as<AnyMap>();
1807 a = cov_map["a"].asDouble();
1808 m = cov_map["m"].asDouble();
1809 E = node.units().convertActivationEnergy(cov_map["E"], "K");
1810 } else {
1811 auto& cov_vec = item.second.asVector<AnyValue>(3);
1812 a = cov_vec[0].asDouble();
1813 m = cov_vec[1].asDouble();
1814 E = node.units().convertActivationEnergy(cov_vec[2], "K");
1815 }
1816 R.coverage_deps[item.first] = CoverageDependency(a, E, m);
1817 }
1818 }
1819 }
1820
getReactions(const XML_Node & node)1821 std::vector<shared_ptr<Reaction> > getReactions(const XML_Node& node)
1822 {
1823 std::vector<shared_ptr<Reaction> > all_reactions;
1824 for (const auto& rxnnode : node.child("reactionData").getChildren("reaction")) {
1825 all_reactions.push_back(newReaction(*rxnnode));
1826 }
1827 return all_reactions;
1828 }
1829
getReactions(const AnyValue & items,Kinetics & kinetics)1830 std::vector<shared_ptr<Reaction>> getReactions(const AnyValue& items,
1831 Kinetics& kinetics)
1832 {
1833 std::vector<shared_ptr<Reaction>> all_reactions;
1834 for (const auto& node : items.asVector<AnyMap>()) {
1835 shared_ptr<Reaction> R(newReaction(node, kinetics));
1836 R->validate();
1837 if (R->valid() && R->checkSpecies(kinetics)) {
1838 all_reactions.emplace_back(R);
1839 }
1840 }
1841 return all_reactions;
1842 }
1843
1844 }
1845