12#include "cantera/kinetics/Falloff.h"
19#include <boost/algorithm/string/predicate.hpp>
22#include <boost/algorithm/string.hpp>
24namespace ba = boost::algorithm;
31 shared_ptr<ReactionRate> rate_,
32 shared_ptr<ThirdBody> tbody_)
33 : reactants(reactants_)
35 , m_from_composition(true)
36 , m_third_body(tbody_)
38 if (reactants.count(
"M") || products.count(
"M")) {
39 throw CanteraError(
"Reaction::Reaction",
40 "Third body 'M' must not be included in either reactant or product maps.");
46 for (
const auto& [name, stoich] : reactants) {
47 if (products.count(name)) {
48 third[name] = products.at(name) - stoich;
52 string name = tbody_->name();
53 if (reactants.count(name) && products.count(name)) {
54 throw CanteraError(
"Reaction::Reaction",
55 "'{}' not acting as third body collider must not be included in both "
56 "reactant and product maps.", name);
59 m_third_body->explicit_3rd =
true;
61 }
else if (!tbody_ && third.size() == 1) {
63 string name = third.begin()->first;
64 m_third_body = make_shared<ThirdBody>(name);
66 m_third_body->explicit_3rd =
true;
71Reaction::Reaction(
const string& equation,
72 shared_ptr<ReactionRate> rate_,
73 shared_ptr<ThirdBody> tbody_)
74 : m_third_body(tbody_)
76 setEquation(equation);
78 if (m_third_body && m_third_body->name() !=
"M") {
79 m_third_body->explicit_3rd =
true;
85 string rate_type = node.
getString(
"type",
"Arrhenius");
88 "Cannot instantiate Reaction with empty Kinetics object.");
91 setParameters(node, kin);
99 if (ba::starts_with(rate_type,
"three-body-")) {
101 rateNode[
"type"] = rate_type.substr(11, rate_type.size() - 11);
108 if (rateNode.
hasKey(
"rate-constant")) {
109 if (!ba::starts_with(rate_type,
"interface-")) {
110 rateNode[
"type"] =
"interface-" + rate_type;
112 }
else if (node.
hasKey(
"sticking-coefficient")) {
113 if (!ba::starts_with(rate_type,
"sticking-")) {
114 rateNode[
"type"] =
"sticking-" + rate_type;
118 "Unable to infer interface reaction type.");
125void Reaction::check()
127 if (!allow_nonreactant_orders) {
128 for (
const auto& [name, order] : orders) {
129 if (reactants.find(name) == reactants.end()) {
131 "Reaction order specified for non-reactant species '{}'", name);
136 if (!allow_negative_orders) {
137 for (
const auto& [name, order] : orders) {
140 "Negative reaction order specified for species '{}'", name);
149 if (reversible && !orders.empty()) {
151 "Reaction orders may only be given for irreversible reactions");
157 m_rate->check(equation());
161AnyMap Reaction::parameters(
bool withInput)
const
169 static bool reg = AnyMap::addOrderingRules(
"Reaction",
171 {
"head",
"equation"},
172 {
"tail",
"duplicate"},
174 {
"tail",
"negative-orders"},
175 {
"tail",
"nonreactant-orders"}
178 out[
"__type__"] =
"Reaction";
183void Reaction::getParameters(
AnyMap& reactionNode)
const
187 "Serialization of empty Reaction object is not supported.");
190 reactionNode[
"equation"] = equation();
193 reactionNode[
"duplicate"] =
true;
196 reactionNode[
"orders"] = orders;
198 if (allow_negative_orders) {
199 reactionNode[
"negative-orders"] =
true;
201 if (allow_nonreactant_orders) {
202 reactionNode[
"nonreactant-orders"] =
true;
205 reactionNode.
update(m_rate->parameters());
208 string rtype = reactionNode[
"type"].asString();
209 if (rtype ==
"pressure-dependent-Arrhenius") {
211 }
else if (m_explicit_type && ba::ends_with(rtype,
"Arrhenius")) {
214 reactionNode[
"type"] =
"three-body";
216 reactionNode[
"type"] =
"elementary";
218 }
else if (ba::ends_with(rtype,
"Arrhenius")) {
219 reactionNode.
erase(
"type");
220 }
else if (m_explicit_type) {
221 reactionNode[
"type"] = type();
222 }
else if (ba::ends_with(rtype,
"Blowers-Masel")) {
223 reactionNode[
"type"] =
"Blowers-Masel";
227 m_third_body->getParameters(reactionNode);
235 "Cannot set reaction parameters from empty node.");
240 setEquation(node[
"equation"].asString(), &kin);
242 if (node.
hasKey(
"orders")) {
243 for (
const auto& [name, order] : node[
"orders"].asMap<
double>()) {
244 orders[name] = order;
253 duplicate = node.
getBool(
"duplicate",
false);
254 allow_negative_orders = node.
getBool(
"negative-orders",
false);
255 allow_nonreactant_orders = node.
getBool(
"nonreactant-orders",
false);
258 m_third_body->setParameters(node);
259 if (m_third_body->name() ==
"M" && m_third_body->efficiencies.size() == 1) {
260 m_third_body->explicit_3rd =
true;
262 }
else if (node.
hasKey(
"default-efficiency") || node.
hasKey(
"efficiencies")) {
264 "Reaction '{}' specifies efficiency parameters\n"
265 "but does not involve third body colliders.", equation());
269void Reaction::setRate(shared_ptr<ReactionRate> rate)
273 "Reaction rate for reaction '{}' must not be empty.", equation());
277 string rate_type = m_rate->type();
279 if (rate_type ==
"falloff" || rate_type ==
"chemically-activated") {
280 if (m_third_body->mass_action && !m_from_composition) {
282 "Third-body collider does not use '(+{})' notation.",
283 m_third_body->name());
285 m_third_body->mass_action =
false;
286 }
else if (rate_type ==
"Chebyshev") {
287 if (m_third_body->name() ==
"M") {
288 warn_deprecated(
"Chebyshev reaction equation", input,
"Specifying 'M' "
289 "in the reaction equation for Chebyshev reactions is deprecated.");
290 m_third_body.reset();
292 }
else if (rate_type ==
"pressure-dependent-Arrhenius") {
293 if (m_third_body->name() ==
"M") {
295 "Found superfluous '{}' in pressure-dependent-Arrhenius reaction.",
296 m_third_body->name());
300 if (rate_type ==
"falloff" || rate_type ==
"chemically-activated") {
301 if (!m_from_composition) {
303 "Reaction equation for falloff reaction '{}'\n does not "
304 "contain valid pressure-dependent third body", equation());
306 m_third_body = make_shared<ThirdBody>(
"(+M)");
311string Reaction::reactantString()
const
313 std::ostringstream result;
314 for (
auto iter = reactants.begin(); iter != reactants.end(); ++iter) {
315 if (iter != reactants.begin()) {
318 if (iter->second != 1.0) {
319 result << iter->second <<
" ";
321 result << iter->first;
324 result << m_third_body->collider();
329string Reaction::productString()
const
331 std::ostringstream result;
332 for (
auto iter = products.begin(); iter != products.end(); ++iter) {
333 if (iter != products.begin()) {
336 if (iter->second != 1.0) {
337 result << iter->second <<
" ";
339 result << iter->first;
342 result << m_third_body->collider();
347string Reaction::equation()
const
350 return reactantString() +
" <=> " + productString();
352 return reactantString() +
" => " + productString();
356void Reaction::setEquation(
const string& equation,
const Kinetics* kin)
360 string rate_type = input.getString(
"type",
"");
361 if (ba::starts_with(rate_type,
"three-body")) {
363 m_explicit_type =
true;
364 }
else if (rate_type ==
"elementary") {
366 m_explicit_type =
true;
376 for (
const auto& [name, stoich] : reactants) {
378 if (products.count(name)) {
380 size_t generic = third_body ==
"M"
381 || third_body ==
"(+M)" || third_body ==
"(+ M)";
384 if (stoich > 1 && products[third_body] > 1) {
392 if (ba::starts_with(rate_type,
"three-body")) {
394 "Reactants for reaction '{}'\n"
395 "do not contain a valid third body collider", equation);
399 }
else if (countM > 1) {
401 "Multiple generic third body colliders 'M' are not supported", equation);
403 }
else if (count > 1) {
409 }
else if (m_third_body) {
411 auto& effs = m_third_body->efficiencies;
412 if (effs.size() != 1 || !reactants.count(effs.begin()->first)) {
414 "Detected ambiguous third body colliders in reaction '{}'\n"
415 "ThirdBody object needs to specify a single species", equation);
417 third_body = effs.begin()->first;
418 m_third_body->explicit_3rd =
true;
419 }
else if (input.hasKey(
"efficiencies")) {
421 auto effs = input[
"efficiencies"].asMap<
double>();
422 if (effs.size() != 1 || !reactants.count(effs.begin()->first)) {
424 "Detected ambiguous third body colliders in reaction '{}'\n"
425 "Collision efficiencies need to specify single species", equation);
427 third_body = effs.begin()->first;
428 m_third_body = make_shared<ThirdBody>(third_body);
429 m_third_body->explicit_3rd =
true;
430 }
else if (input.hasKey(
"default-efficiency")) {
433 "Detected ambiguous third body colliders in reaction '{}'\n"
434 "Third-body definition requires specification of efficiencies",
436 }
else if (ba::starts_with(rate_type,
"three-body")) {
439 "Detected ambiguous third body colliders in reaction '{}'\n"
440 "A valid ThirdBody or collision efficiency definition is required",
446 }
else if (third_body !=
"M" && !ba::starts_with(rate_type,
"three-body")
447 && !ba::starts_with(third_body,
"(+"))
456 for (
const auto& [name, stoich] : reactants) {
457 if (trunc(stoich) != stoich) {
460 nreac +=
static_cast<size_t>(stoich);
464 for (
const auto& [name, stoich] : products) {
465 if (trunc(stoich) != stoich) {
468 nprod +=
static_cast<size_t>(stoich);
472 if (nreac != 3 && nprod != 3) {
478 string tName = m_third_body->name();
479 if (tName != third_body && third_body !=
"M" && tName !=
"M") {
481 "Detected incompatible third body colliders in reaction '{}'\n"
482 "ThirdBody definition does not match equation", equation);
484 m_third_body->setName(third_body);
486 m_third_body = make_shared<ThirdBody>(third_body);
490 auto reac = reactants.find(third_body);
491 if (trunc(reac->second) != 1) {
494 reactants.erase(reac);
498 auto prod = products.find(third_body);
499 if (trunc(prod->second) != 1) {
502 products.erase(prod);
506string Reaction::type()
const
509 throw CanteraError(
"Reaction::type",
"Empty Reaction does not have a type");
512 string rate_type = m_rate->type();
513 string sub_type = m_rate->subType();
514 if (sub_type !=
"") {
515 return rate_type +
"-" + sub_type;
519 return "three-body-" + rate_type;
535 UnitStack rate_units(rxn_phase.standardConcentrationUnits());
541 for (
const auto& [name, order] : orders) {
544 rate_units.
update(phase.standardConcentrationUnits(), -order);
546 for (
const auto& [name, stoich] : reactants) {
549 if (name ==
"M" || ba::starts_with(name,
"(+")) {
553 }
else if (orders.find(name) == orders.end()) {
556 rate_units.
update(phase.standardConcentrationUnits(), -stoich);
560 if (m_third_body && m_third_body->mass_action) {
565 Reaction::rate_units = rate_units.
product();
569void updateUndeclared(vector<string>& undeclared,
572 for (
const auto& [name, stoich]: comp) {
574 undeclared.emplace_back(name);
579void Reaction::checkBalance(
const Kinetics& kin)
const
584 for (
const auto& [name, stoich] : products) {
587 for (
size_t m = 0; m < ph.
nElements(); m++) {
592 for (
const auto& [name, stoich] : reactants) {
595 for (
size_t m = 0; m < ph.
nElements(); m++) {
602 for (
const auto& [elem, balance] : balr) {
603 double elemsum = balr[elem] + balp[elem];
604 double elemdiff = fabs(balp[elem] - balr[elem]);
605 if (elemsum > 0.0 && elemdiff / elemsum > 1e-4) {
607 msg += fmt::format(
" {} {} {}\n",
608 elem, balr[elem], balp[elem]);
613 "The following reaction is unbalanced: {}\n"
614 " Element Reactants Products\n{}",
623 double reac_sites = 0.0;
624 double prod_sites = 0.0;
626 for (
const auto& [name, stoich] : reactants) {
629 reac_sites += stoich * surf.size(k);
632 for (
const auto& [name, stoich] : products) {
633 size_t k = surf.speciesIndex(name);
635 prod_sites += stoich * surf.size(k);
638 if (fabs(reac_sites - prod_sites) > 1e-5 * (reac_sites + prod_sites)) {
640 "Number of surface sites not balanced in reaction {}.\n"
641 "Reactant sites: {}\nProduct sites: {}",
642 equation(), reac_sites, prod_sites);
646bool Reaction::checkSpecies(
const Kinetics& kin)
const
649 vector<string> undeclared;
650 updateUndeclared(undeclared, reactants, kin);
651 updateUndeclared(undeclared, products, kin);
652 if (!undeclared.empty()) {
656 throw InputFileError(
"Reaction::checkSpecies", input,
"Reaction '{}'\n"
657 "contains undeclared species: '{}'",
658 equation(), ba::join(undeclared,
"', '"));
663 updateUndeclared(undeclared, orders, kin);
664 if (!undeclared.empty()) {
668 if (input.hasKey(
"orders")) {
671 "defines reaction orders for undeclared species: '{}'",
672 equation(), ba::join(undeclared,
"', '"));
675 throw InputFileError(
"Reaction::checkSpecies", input,
"Reaction '{}'\n"
676 "defines reaction orders for undeclared species: '{}'",
677 equation(), ba::join(undeclared,
"', '"));
682 return m_third_body->checkSpecies(*
this, kin);
690bool Reaction::usesElectrochemistry(
const Kinetics& kin)
const
693 vector<double> e_counter(kin.
nPhases(), 0.0);
696 for (
const auto& [name, stoich] : products) {
704 for (
const auto& [name, stoich] : reactants) {
712 for (
double delta_e : e_counter) {
713 if (std::abs(delta_e) > 1e-4) {
721ThirdBody::ThirdBody(
double default_eff)
722 : default_efficiency(default_eff)
725 "Instantiation with default efficiency is deprecated and will be removed "
726 "after Cantera 3.0. Instantiate with collider name instead.");
729ThirdBody::ThirdBody(
const string& third_body)
736 string name = third_body;
737 if (ba::starts_with(third_body,
"(+ ")) {
739 name = third_body.substr(3, third_body.size() - 4);
740 }
else if (ba::starts_with(third_body,
"(+")) {
742 name = third_body.substr(2, third_body.size() - 3);
755 "Conflicting efficiency definition for explicit third body '{}'",
name);
762ThirdBody::ThirdBody(
const AnyMap& node)
770 "To be removed after Cantera 3.0. Renamed to setParameters");
776 if (node.
hasKey(
"default-efficiency")) {
777 double value = node[
"default-efficiency"].asDouble();
778 if (
m_name !=
"M" && value != 0.) {
779 throw InputFileError(
"ThirdBody::setParameters", node[
"default-efficiency"],
780 "Invalid default efficiency for explicit collider {};\n"
781 "value is optional and/or needs to be zero",
m_name);
785 if (node.
hasKey(
"efficiencies")) {
792 "Detected incompatible third body colliders definitions");
819 return " (+" +
m_name +
")";
824 vector<string> undeclared;
827 if (!undeclared.empty()) {
831 rxn.
input[
"efficiencies"],
"Reaction '{}'\n"
832 "defines third-body efficiencies for undeclared species: '{}'",
833 rxn.
equation(), ba::join(undeclared,
"', '"));
837 "is a three-body reaction with undeclared species: '{}'",
838 rxn.
equation(), ba::join(undeclared,
"', '"));
848ThreeBodyReaction::ThreeBodyReaction()
851 "To be removed after Cantera 3.0. Replaceable with Reaction.");
856ThreeBodyReaction::ThreeBodyReaction(
const Composition& reactants,
858 const ArrheniusRate& rate,
859 const ThirdBody& tbody)
860 : Reaction(reactants,
862 make_shared<ArrheniusRate>(rate),
863 make_shared<ThirdBody>(tbody))
866 "To be removed after Cantera 3.0. Replaceable with Reaction.");
869ThreeBodyReaction::ThreeBodyReaction(
const AnyMap& node,
const Kinetics& kin)
870 : Reaction(node, kin)
873 "To be removed after Cantera 3.0. Replaceable with Reaction.");
876FalloffReaction::FalloffReaction()
879 "To be removed after Cantera 3.0. Replaceable with Reaction.");
884FalloffReaction::FalloffReaction(
const Composition& reactants_,
886 const FalloffRate& rate_,
887 const ThirdBody& tbody_)
890 "To be removed after Cantera 3.0. Replaceable with Reaction.");
895 AnyMap node = rate_.parameters();
897 string rate_type = node[
"type"].asString();
901FalloffReaction::FalloffReaction(
const AnyMap& node,
const Kinetics& kin)
902 : Reaction(node, kin)
905 "To be removed after Cantera 3.0. Replaceable with Reaction.");
907 m_third_body = make_shared<ThirdBody>();
914 return make_unique<Reaction>();
919 return make_unique<Reaction>(rxn_node, kin);
927 vector<string> tokens;
929 tokens.push_back(
"+");
931 size_t last_used =
npos;
932 bool reactants =
true;
933 for (
size_t i = 1; i < tokens.size(); i++) {
934 if (tokens[i] ==
"+" || ba::starts_with(tokens[i],
"(+") ||
935 tokens[i] ==
"<=>" || tokens[i] ==
"=" || tokens[i] ==
"=>") {
936 string species = tokens[i-1];
939 bool mass_action =
true;
940 if (last_used !=
npos && tokens[last_used] ==
"(+"
941 && ba::ends_with(species,
")")) {
944 species =
"(+" + species;
945 }
else if (last_used == i - 1 && ba::starts_with(species,
"(+")
946 && ba::ends_with(species,
")")) {
949 }
else if (last_used == i - 2) {
951 }
else if (last_used == i - 3) {
961 "Error parsing reaction string '{}'.\n"
962 "Current token: '{}'\nlast_used: '{}'",
964 (last_used ==
npos) ?
"n/a" : tokens[last_used]);
967 && mass_action && species !=
"M"))
982 if (tokens[i] ==
"<=>" || tokens[i] ==
"=") {
985 }
else if (tokens[i] ==
"=>") {
994 vector<shared_ptr<Reaction>> all_reactions;
996 auto R = make_shared<Reaction>(node, kinetics);
998 R->validate(kinetics);
999 if (R->valid() && R->checkSpecies(kinetics)) {
1000 all_reactions.emplace_back(R);
1003 return all_reactions;
Header file for class Cantera::Array2D.
Header for reaction rates that involve Arrhenius-type kinetics.
Base class for kinetics managers and also contains the kineticsmgr module documentation (see Kinetics...
Factory class for reaction rate objects.
Header for a simple thermodynamics model of a surface phase derived from ThermoPhase,...
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
Base class defining common data possessed by both AnyMap and AnyValue objects.
A map of string keys to values whose type can vary at runtime.
void copyMetadata(const AnyMap &other)
Copy metadata including input line/column from an existing AnyMap.
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
bool empty() const
Return boolean indicating whether AnyMap is empty.
void setFlowStyle(bool flow=true)
Use "flow" style when outputting this AnyMap to YAML.
bool getBool(const string &key, bool default_) const
If key exists, return it as a bool, otherwise return default_.
const string & getString(const string &key, const string &default_) const
If key exists, return it as a string, otherwise return default_.
void erase(const string &key)
Erase the value held by key.
void update(const AnyMap &other, bool keepExisting=true)
Add items from other to this AnyMap.
A wrapper for a variable whose type is determined at runtime.
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.
Base class for exceptions thrown by Cantera classes.
virtual string getMessage() const
Method overridden by derived classes to format the error message.
Public interface for kinetics managers.
size_t reactionPhaseIndex() const
Phase where the reactions occur.
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
size_t nPhases() const
The number of phases participating in the reaction mechanism.
void skipUndeclaredSpecies(bool skip)
Determine behavior when adding a new reaction that contains species not defined in any of the phases ...
void skipUndeclaredThirdBodies(bool skip)
Determine behavior when adding a new reaction that contains third-body efficiencies for species not d...
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
ThermoPhase & speciesPhase(const string &nm)
This function looks up the name of a species and returns a reference to the ThermoPhase object of the...
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 ...
size_t nDim() const
Returns the number of spatial dimensions (1, 2, or 3)
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
double nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
size_t nElements() const
Number of elements.
string elementName(size_t m) const
Name of the element with index m.
double charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Abstract base class which stores data about a reaction and its rate parameterization so that it can b...
void setRate(shared_ptr< ReactionRate > rate)
Set reaction rate pointer.
shared_ptr< ThirdBody > m_third_body
Relative efficiencies of third-body species in enhancing the reaction rate (if applicable)
bool reversible
True if the current reaction is reversible. False otherwise.
string type() const
The type of reaction, including reaction rate information.
string equation() const
The chemical equation for this reaction.
Composition products
Product species and stoichiometric coefficients.
Composition reactants
Reactant species and stoichiometric coefficients.
AnyMap input
Input data used for specific models.
void setValid(bool valid)
Set validity flag of reaction.
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
Base class for a phase with thermodynamic properties.
bool explicit_3rd
Flag indicating whether third body requires explicit serialization.
void setParameters(const AnyMap &node)
Set third-body efficiencies from AnyMap node
void setEfficiencies(const AnyMap &node)
Set third-body efficiencies from AnyMap node
double efficiency(const string &k) const
Get the third-body efficiency for species k
double default_efficiency
The default third body efficiency for species not listed in efficiencies.
Composition efficiencies
Map of species to third body efficiency.
string collider() const
Suffix representing the third body collider in reaction equation, for example + M or (+M)
bool mass_action
Third body is used by law of mass action (true for three-body reactions, false for falloff reactions)
string m_name
Name of the third body collider.
void getParameters(AnyMap &node) const
Get third-body efficiencies from AnyMap node
void setName(const string &third_body)
Set name of the third body collider.
bool checkSpecies(const Reaction &rxn, const Kinetics &kin) const
Verify that all species involved in collision efficiencies are defined in the Kinetics object.
string name() const
Name of the third body collider.
A representation of the units associated with a dimensional quantity.
double fpValueCheck(const string &val)
Translate a string into one double value, with error checking.
void tokenizeString(const string &in_val, vector< string > &v)
This function separates a string up into tokens according to the location of white space.
shared_ptr< ReactionRate > newReactionRate(const string &type)
Create a new empty ReactionRate object.
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
void parseReactionEquation(Reaction &R, const string &equation, const AnyBase &reactionNode, const Kinetics *kin)
Parse reaction equation.
vector< shared_ptr< Reaction > > getReactions(const AnyValue &items, Kinetics &kinetics)
Create Reaction objects for each item (an AnyMap) in items.
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
const U & getValue(const map< T, U > &m, const T &key, const U &default_val)
Const accessor for a value in a map.
map< string, double > Composition
Map from string names to doubles.
unique_ptr< Reaction > newReaction(const string &type)
Create a new empty Reaction object.
Contains declarations for string manipulation functions within Cantera.
Unit aggregation utility.
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...
Units product() const
Calculate product of units-exponent stack.
void join(double exponent)
Join (update) exponent of standard units, where the updated exponent is the sum of the pre-existing e...
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...