17#include <boost/algorithm/string/predicate.hpp>
20#include <boost/algorithm/string.hpp>
22namespace ba = boost::algorithm;
29 shared_ptr<ReactionRate> rate_,
30 shared_ptr<ThirdBody> tbody_)
31 : reactants(reactants_)
33 , m_from_composition(true)
34 , m_third_body(tbody_)
36 if (reactants.count(
"M") || products.count(
"M")) {
37 throw CanteraError(
"Reaction::Reaction",
38 "Third body 'M' must not be included in either reactant or product maps.");
44 for (
const auto& [name, stoich] : reactants) {
45 if (products.count(name)) {
46 third[name] = products.at(name) - stoich;
50 string name = tbody_->name();
51 if (reactants.count(name) && products.count(name)) {
52 throw CanteraError(
"Reaction::Reaction",
53 "'{}' not acting as third body collider must not be included in both "
54 "reactant and product maps.", name);
57 m_third_body->explicit_3rd =
true;
59 }
else if (!tbody_ && third.size() == 1) {
61 string name = third.begin()->first;
62 m_third_body = make_shared<ThirdBody>(name);
64 m_third_body->explicit_3rd =
true;
70Reaction::Reaction(
const string& equation,
71 shared_ptr<ReactionRate> rate_,
72 shared_ptr<ThirdBody> tbody_)
73 : m_third_body(tbody_)
76 setEquation(equation);
77 if (m_third_body && m_third_body->name() !=
"M") {
78 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;
116 }
else if (rate_type ==
"Arrhenius") {
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");
160 m_rate->check(equation());
162 string rate_type = m_rate->type();
164 if (rate_type ==
"falloff" || rate_type ==
"chemically-activated") {
165 if (m_third_body->mass_action && !m_from_composition) {
167 "Third-body collider does not use '(+{})' notation.",
168 m_third_body->name());
170 m_third_body->mass_action =
false;
171 }
else if (rate_type ==
"Chebyshev") {
172 if (m_third_body->name() ==
"M") {
173 warn_deprecated(
"Chebyshev reaction equation", input,
"Specifying 'M' "
174 "in the reaction equation for Chebyshev reactions is deprecated.");
175 m_third_body.reset();
177 }
else if (rate_type ==
"pressure-dependent-Arrhenius") {
178 if (m_third_body->name() ==
"M") {
180 "Found superfluous '{}' in pressure-dependent-Arrhenius reaction.",
181 m_third_body->name());
185 if (rate_type ==
"falloff" || rate_type ==
"chemically-activated") {
186 if (!m_from_composition) {
188 "Reaction equation for falloff reaction '{}'\n does not "
189 "contain valid pressure-dependent third body", equation());
191 m_third_body = make_shared<ThirdBody>(
"(+M)");
196AnyMap Reaction::parameters(
bool withInput)
const
204 static bool reg = AnyMap::addOrderingRules(
"Reaction",
206 {
"head",
"equation"},
207 {
"tail",
"duplicate"},
209 {
"tail",
"negative-orders"},
210 {
"tail",
"nonreactant-orders"}
213 out[
"__type__"] =
"Reaction";
218void Reaction::getParameters(
AnyMap& reactionNode)
const
222 "Serialization of empty Reaction object is not supported.");
225 reactionNode[
"equation"] = equation();
228 reactionNode[
"duplicate"] =
true;
231 reactionNode[
"orders"] = orders;
233 if (allow_negative_orders) {
234 reactionNode[
"negative-orders"] =
true;
236 if (allow_nonreactant_orders) {
237 reactionNode[
"nonreactant-orders"] =
true;
240 reactionNode.
update(m_rate->parameters());
243 string rtype = reactionNode[
"type"].asString();
244 if (rtype ==
"pressure-dependent-Arrhenius") {
246 }
else if (m_explicit_type && ba::ends_with(rtype,
"Arrhenius")) {
249 reactionNode[
"type"] =
"three-body";
251 reactionNode[
"type"] =
"elementary";
253 }
else if (ba::ends_with(rtype,
"Arrhenius")) {
254 reactionNode.
erase(
"type");
255 }
else if (m_explicit_type) {
256 reactionNode[
"type"] = type();
257 }
else if (ba::ends_with(rtype,
"Blowers-Masel")) {
258 reactionNode[
"type"] =
"Blowers-Masel";
262 m_third_body->getParameters(reactionNode);
270 "Cannot set reaction parameters from empty node.");
275 setEquation(node[
"equation"].asString(), &kin);
277 if (node.
hasKey(
"orders")) {
278 for (
const auto& [name, order] : node[
"orders"].asMap<
double>()) {
279 orders[name] = order;
288 duplicate = node.
getBool(
"duplicate",
false);
289 allow_negative_orders = node.
getBool(
"negative-orders",
false);
290 allow_nonreactant_orders = node.
getBool(
"nonreactant-orders",
false);
293 m_third_body->setParameters(node);
294 if (m_third_body->name() ==
"M" && m_third_body->efficiencies.size() == 1) {
295 m_third_body->explicit_3rd =
true;
297 }
else if (node.
hasKey(
"default-efficiency") || node.
hasKey(
"efficiencies")) {
299 "Reaction '{}' specifies efficiency parameters\n"
300 "but does not involve third body colliders.", equation());
304void Reaction::setRate(shared_ptr<ReactionRate> rate)
308 "Reaction rate for reaction '{}' must not be empty.", equation());
313string Reaction::reactantString()
const
315 std::ostringstream result;
316 for (
auto iter = reactants.begin(); iter != reactants.end(); ++iter) {
317 if (iter != reactants.begin()) {
320 if (iter->second != 1.0) {
321 result << iter->second <<
" ";
323 result << iter->first;
326 result << m_third_body->collider();
331string Reaction::productString()
const
333 std::ostringstream result;
334 for (
auto iter = products.begin(); iter != products.end(); ++iter) {
335 if (iter != products.begin()) {
338 if (iter->second != 1.0) {
339 result << iter->second <<
" ";
341 result << iter->first;
344 result << m_third_body->collider();
349string Reaction::equation()
const
352 return reactantString() +
" <=> " + productString();
354 return reactantString() +
" => " + productString();
358void Reaction::setEquation(
const string& equation,
const Kinetics* kin)
361 string rate_type = (m_rate) ? m_rate->type() : input.getString(
"type",
"");
362 if (ba::starts_with(rate_type,
"three-body")) {
364 m_explicit_type =
true;
365 }
else if (rate_type ==
"elementary") {
367 m_explicit_type =
true;
369 }
else if (kin && kin->
thermo(0).
nDim() != 3) {
372 }
else if (rate_type ==
"electron-collision-plasma") {
380 for (
const auto& [name, stoich] : reactants) {
382 if (products.count(name)) {
384 size_t generic = third_body ==
"M"
385 || third_body ==
"(+M)" || third_body ==
"(+ M)";
388 if (stoich > 1 && products[third_body] > 1) {
396 if (ba::starts_with(rate_type,
"three-body")) {
398 "Reactants for reaction '{}'\n"
399 "do not contain a valid third body collider", equation);
403 }
else if (countM > 1) {
405 "Multiple generic third body colliders 'M' are not supported", equation);
407 }
else if (count > 1) {
413 }
else if (m_third_body) {
415 auto& effs = m_third_body->efficiencies;
416 if (effs.size() != 1 || !reactants.count(effs.begin()->first)) {
418 "Detected ambiguous third body colliders in reaction '{}'\n"
419 "ThirdBody object needs to specify a single species", equation);
421 third_body = effs.begin()->first;
422 m_third_body->explicit_3rd =
true;
423 }
else if (input.hasKey(
"efficiencies")) {
425 auto effs = input[
"efficiencies"].asMap<
double>();
426 if (effs.size() != 1 || !reactants.count(effs.begin()->first)) {
428 "Detected ambiguous third body colliders in reaction '{}'\n"
429 "Collision efficiencies need to specify single species", equation);
431 third_body = effs.begin()->first;
432 m_third_body = make_shared<ThirdBody>(third_body);
433 m_third_body->explicit_3rd =
true;
434 }
else if (input.hasKey(
"default-efficiency")) {
437 "Detected ambiguous third body colliders in reaction '{}'\n"
438 "Third-body definition requires specification of efficiencies",
440 }
else if (ba::starts_with(rate_type,
"three-body")) {
443 "Detected ambiguous third body colliders in reaction '{}'\n"
444 "A valid ThirdBody or collision efficiency definition is required",
450 }
else if (third_body !=
"M" && !ba::starts_with(rate_type,
"three-body")
451 && !ba::starts_with(third_body,
"(+"))
460 for (
const auto& [name, stoich] : reactants) {
461 if (trunc(stoich) != stoich) {
464 nreac +=
static_cast<size_t>(stoich);
468 for (
const auto& [name, stoich] : products) {
469 if (trunc(stoich) != stoich) {
472 nprod +=
static_cast<size_t>(stoich);
476 if (nreac != 3 && nprod != 3) {
482 string tName = m_third_body->name();
483 if (tName != third_body && third_body !=
"M" && tName !=
"M") {
485 "Detected incompatible third body colliders in reaction '{}'\n"
486 "ThirdBody definition does not match equation", equation);
488 m_third_body->setName(third_body);
490 m_third_body = make_shared<ThirdBody>(third_body);
494 auto reac = reactants.find(third_body);
495 if (trunc(reac->second) != 1) {
498 reactants.erase(reac);
502 auto prod = products.find(third_body);
503 if (trunc(prod->second) != 1) {
506 products.erase(prod);
510string Reaction::type()
const
513 throw CanteraError(
"Reaction::type",
"Empty Reaction does not have a type");
516 string rate_type = m_rate->type();
517 string sub_type = m_rate->subType();
518 if (sub_type !=
"") {
519 return rate_type +
"-" + sub_type;
523 return "three-body-" + rate_type;
544 for (
const auto& [name, order] : orders) {
547 rate_units.
update(phase.standardConcentrationUnits(), -order);
549 for (
const auto& [name, stoich] : reactants) {
552 if (name ==
"M" || ba::starts_with(name,
"(+")) {
556 }
else if (orders.find(name) == orders.end()) {
559 rate_units.
update(phase.standardConcentrationUnits(), -stoich);
563 if (m_third_body && m_third_body->mass_action) {
568 Reaction::rate_units = rate_units.
product();
572void updateUndeclared(vector<string>& undeclared,
575 for (
const auto& [name, stoich]: comp) {
577 undeclared.emplace_back(name);
582void Reaction::checkBalance(
const Kinetics& kin)
const
587 for (
const auto& [name, stoich] : products) {
590 for (
size_t m = 0; m < ph.
nElements(); m++) {
595 for (
const auto& [name, stoich] : reactants) {
598 for (
size_t m = 0; m < ph.
nElements(); m++) {
605 for (
const auto& [elem, balance] : balr) {
606 double elemsum = balr[elem] + balp[elem];
607 double elemdiff = fabs(balp[elem] - balr[elem]);
608 if (elemsum > 0.0 && elemdiff / elemsum > 1e-4) {
610 msg += fmt::format(
" {} {} {}\n",
611 elem, balr[elem], balp[elem]);
616 "The following reaction is unbalanced: {}\n"
617 " Element Reactants Products\n{}",
626 double reac_sites = 0.0;
627 double prod_sites = 0.0;
629 for (
const auto& [name, stoich] : reactants) {
632 reac_sites += stoich * surf.size(k);
635 for (
const auto& [name, stoich] : products) {
636 size_t k = surf.speciesIndex(name);
638 prod_sites += stoich * surf.size(k);
641 if (fabs(reac_sites - prod_sites) > 1e-5 * (reac_sites + prod_sites)) {
643 "Number of surface sites not balanced in reaction {}.\n"
644 "Reactant sites: {}\nProduct sites: {}",
645 equation(), reac_sites, prod_sites);
649bool Reaction::checkSpecies(
const Kinetics& kin)
const
652 vector<string> undeclared;
653 updateUndeclared(undeclared, reactants, kin);
654 updateUndeclared(undeclared, products, kin);
655 if (!undeclared.empty()) {
659 throw InputFileError(
"Reaction::checkSpecies", input,
"Reaction '{}'\n"
660 "contains undeclared species: '{}'",
661 equation(), ba::join(undeclared,
"', '"));
666 updateUndeclared(undeclared, orders, kin);
667 if (!undeclared.empty()) {
671 if (input.hasKey(
"orders")) {
674 "defines reaction orders for undeclared species: '{}'",
675 equation(), ba::join(undeclared,
"', '"));
678 throw InputFileError(
"Reaction::checkSpecies", input,
"Reaction '{}'\n"
679 "defines reaction orders for undeclared species: '{}'",
680 equation(), ba::join(undeclared,
"', '"));
685 return m_third_body->checkSpecies(*
this, kin);
693bool Reaction::usesElectrochemistry(
const Kinetics& kin)
const
696 vector<double> e_counter(kin.
nPhases(), 0.0);
699 for (
const auto& [name, stoich] : products) {
707 for (
const auto& [name, stoich] : reactants) {
715 for (
double delta_e : e_counter) {
716 if (std::abs(delta_e) > 1e-4) {
725ThirdBody::ThirdBody(
const string& third_body)
730void ThirdBody::setName(
const string& third_body)
732 string name = third_body;
733 if (ba::starts_with(third_body,
"(+ ")) {
735 name = third_body.substr(3, third_body.size() - 4);
736 }
else if (ba::starts_with(third_body,
"(+")) {
738 name = third_body.substr(2, third_body.size() - 3);
741 if (name == m_name) {
744 if (name ==
"M" && efficiencies.size() == 1) {
749 if (efficiencies.size()) {
751 "Conflicting efficiency definition for explicit third body '{}'", name);
754 default_efficiency = 0.;
755 efficiencies[m_name] = 1.;
758ThirdBody::ThirdBody(
const AnyMap& node)
763void ThirdBody::setParameters(
const AnyMap& node)
765 if (node.
hasKey(
"default-efficiency")) {
766 double value = node[
"default-efficiency"].asDouble();
767 if (m_name !=
"M" && value != 0.) {
768 throw InputFileError(
"ThirdBody::setParameters", node[
"default-efficiency"],
769 "Invalid default efficiency for explicit collider {};\n"
770 "value is optional and/or needs to be zero", m_name);
772 default_efficiency = value;
774 if (node.
hasKey(
"efficiencies")) {
775 efficiencies = node[
"efficiencies"].asMap<
double>();
778 && (efficiencies.size() != 1 || efficiencies.begin()->first != m_name))
781 "Detected incompatible third body colliders definitions");
785void ThirdBody::getParameters(
AnyMap& node)
const
787 if (m_name ==
"M" || explicit_3rd) {
788 if (efficiencies.size()) {
789 node[
"efficiencies"] = efficiencies;
792 if (default_efficiency != 1.0 && !explicit_3rd) {
793 node[
"default-efficiency"] = default_efficiency;
798double ThirdBody::efficiency(
const string& k)
const
800 return getValue(efficiencies, k, default_efficiency);
803string ThirdBody::collider()
const
806 return " + " + m_name;
808 return " (+" + m_name +
")";
813 vector<string> undeclared;
814 updateUndeclared(undeclared, efficiencies, kin);
816 if (!undeclared.empty()) {
820 rxn.
input[
"efficiencies"],
"Reaction '{}'\n"
821 "defines third-body efficiencies for undeclared species: '{}'",
822 rxn.
equation(), ba::join(undeclared,
"', '"));
826 "is a three-body reaction with undeclared species: '{}'",
827 rxn.
equation(), ba::join(undeclared,
"', '"));
840 return make_unique<Reaction>();
845 return make_unique<Reaction>(rxn_node, kin);
853 vector<string> tokens;
855 tokens.push_back(
"+");
857 size_t last_used =
npos;
858 bool reactants =
true;
859 for (
size_t i = 1; i < tokens.size(); i++) {
860 if (tokens[i] ==
"+" || ba::starts_with(tokens[i],
"(+") ||
861 tokens[i] ==
"<=>" || tokens[i] ==
"=" || tokens[i] ==
"=>") {
862 string species = tokens[i-1];
865 bool mass_action =
true;
866 if (last_used !=
npos && tokens[last_used] ==
"(+"
867 && ba::ends_with(species,
")")) {
870 species =
"(+" + species;
871 }
else if (last_used == i - 1 && ba::starts_with(species,
"(+")
872 && ba::ends_with(species,
")")) {
875 }
else if (last_used == i - 2) {
877 }
else if (last_used == i - 3) {
887 "Error parsing reaction string '{}'.\n"
888 "Current token: '{}'\nlast_used: '{}'",
890 (last_used ==
npos) ?
"n/a" : tokens[last_used]);
893 && mass_action && species !=
"M"))
908 if (tokens[i] ==
"<=>" || tokens[i] ==
"=") {
911 }
else if (tokens[i] ==
"=>") {
920 vector<shared_ptr<Reaction>> all_reactions;
922 auto R = make_shared<Reaction>(node, kinetics);
923 R->validate(kinetics);
924 if (R->valid() && R->checkSpecies(kinetics)) {
925 all_reactions.emplace_back(R);
928 return all_reactions;
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.
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.
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...
bool reversible
True if the current reaction is reversible. False otherwise.
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.
virtual Units standardConcentrationUnits() const
Returns the units of the "standard concentration" for this phase.
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...