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 : m_from_composition(true)
32 , m_third_body(tbody_)
34 for (
auto& [species, stoich] : reactants_) {
36 reactants[species] = stoich;
39 for (
auto& [species, stoich] : products_) {
41 products[species] = stoich;
44 if (reactants.count(
"M") || products.count(
"M")) {
45 throw CanteraError(
"Reaction::Reaction",
46 "Third body 'M' must not be included in either reactant or product maps.");
52 for (
const auto& [name, stoich] : reactants) {
53 if (products.count(name)) {
54 third[name] = products.at(name) - stoich;
58 string name = tbody_->name();
59 if (reactants.count(name) && products.count(name)) {
60 throw CanteraError(
"Reaction::Reaction",
61 "'{}' not acting as third body collider must not be included in both "
62 "reactant and product maps.", name);
65 m_third_body->explicit_3rd =
true;
67 }
else if (!tbody_ && third.size() == 1
68 && m_rate->type() !=
"electron-collision-plasma")
71 string name = third.begin()->first;
72 m_third_body = make_shared<ThirdBody>(name);
74 m_third_body->explicit_3rd =
true;
80Reaction::Reaction(
const string& equation,
81 shared_ptr<ReactionRate> rate_,
82 shared_ptr<ThirdBody> tbody_)
83 : m_third_body(tbody_)
86 setEquation(equation);
87 if (m_third_body && m_third_body->name() !=
"M") {
88 m_third_body->explicit_3rd =
true;
95 string rate_type = node.
getString(
"type",
"Arrhenius");
98 "Cannot instantiate Reaction with empty Kinetics object.");
101 setParameters(node, kin);
109 if (ba::starts_with(rate_type,
"three-body-")) {
111 rateNode[
"type"] = rate_type.substr(11, rate_type.size() - 11);
118 if (rateNode.
hasKey(
"rate-constant")) {
119 if (!ba::starts_with(rate_type,
"interface-")) {
120 rateNode[
"type"] =
"interface-" + rate_type;
122 }
else if (node.
hasKey(
"sticking-coefficient")) {
123 if (!ba::starts_with(rate_type,
"sticking-")) {
124 rateNode[
"type"] =
"sticking-" + rate_type;
126 }
else if (rate_type ==
"Arrhenius") {
128 "Unable to infer interface reaction type.");
135void Reaction::check()
137 if (!allow_nonreactant_orders) {
138 for (
const auto& [name, order] : orders) {
139 if (reactants.find(name) == reactants.end()) {
141 "Reaction order specified for non-reactant species '{}'", name);
146 if (!allow_negative_orders) {
147 for (
const auto& [name, order] : orders) {
150 "Negative reaction order specified for species '{}'", name);
159 if (reversible && !orders.empty()) {
161 "Reaction orders may only be given for irreversible reactions");
170 m_rate->check(equation());
172 string rate_type = m_rate->type();
174 if (rate_type ==
"falloff" || rate_type ==
"chemically-activated") {
175 if (m_third_body->mass_action && !m_from_composition) {
177 "Third-body collider does not use '(+{})' notation.",
178 m_third_body->name());
180 m_third_body->mass_action =
false;
181 }
else if (rate_type ==
"Chebyshev") {
182 if (m_third_body->name() ==
"M") {
183 warn_deprecated(
"Chebyshev reaction equation", input,
"Specifying 'M' "
184 "in the reaction equation for Chebyshev reactions is deprecated.");
185 m_third_body.reset();
187 }
else if (rate_type ==
"pressure-dependent-Arrhenius") {
188 if (m_third_body->name() ==
"M") {
190 "Found superfluous '{}' in pressure-dependent-Arrhenius reaction.",
191 m_third_body->name());
195 if (rate_type ==
"falloff" || rate_type ==
"chemically-activated") {
196 if (!m_from_composition) {
198 "Reaction equation for falloff reaction '{}'\n does not "
199 "contain valid pressure-dependent third body", equation());
201 m_third_body = make_shared<ThirdBody>(
"(+M)");
206AnyMap Reaction::parameters(
bool withInput)
const
214 static bool reg = AnyMap::addOrderingRules(
"Reaction",
215 {{
"head",
"equation"},
217 {
"tail",
"duplicate"},
219 {
"tail",
"negative-orders"},
220 {
"tail",
"nonreactant-orders"}
223 out[
"__type__"] =
"Reaction";
228void Reaction::getParameters(
AnyMap& reactionNode)
const
232 "Serialization of empty Reaction object is not supported.");
235 reactionNode[
"equation"] = equation();
238 reactionNode[
"duplicate"] =
true;
240 reactionNode.
exclude(
"duplicate");
243 reactionNode[
"orders"] = orders;
245 reactionNode.
exclude(
"orders");
247 if (allow_negative_orders) {
248 reactionNode[
"negative-orders"] =
true;
250 reactionNode.
exclude(
"negative-orders");
252 if (allow_nonreactant_orders) {
253 reactionNode[
"nonreactant-orders"] =
true;
255 reactionNode.
exclude(
"nonreactant-orders");
258 reactionNode.
update(m_rate->parameters());
261 string rtype = reactionNode[
"type"].asString();
262 if (rtype ==
"pressure-dependent-Arrhenius") {
264 }
else if (m_explicit_type && ba::ends_with(rtype,
"Arrhenius")) {
267 reactionNode[
"type"] =
"three-body";
269 reactionNode[
"type"] =
"elementary";
271 }
else if (ba::ends_with(rtype,
"Arrhenius")) {
273 }
else if (m_explicit_type) {
274 reactionNode[
"type"] = type();
275 }
else if (ba::ends_with(rtype,
"Blowers-Masel")) {
276 reactionNode[
"type"] =
"Blowers-Masel";
280 m_third_body->getParameters(reactionNode);
288 "Cannot set reaction parameters from empty node.");
293 setEquation(node[
"equation"].asString(), &kin);
295 if (node.
hasKey(
"orders")) {
296 for (
const auto& [name, order] : node[
"orders"].asMap<
double>()) {
297 orders[name] = order;
306 duplicate = node.
getBool(
"duplicate",
false);
307 allow_negative_orders = node.
getBool(
"negative-orders",
false);
308 allow_nonreactant_orders = node.
getBool(
"nonreactant-orders",
false);
311 m_third_body->setParameters(node);
312 if (m_third_body->name() ==
"M" && m_third_body->efficiencies.size() == 1) {
313 m_third_body->explicit_3rd =
true;
315 }
else if (node.
hasKey(
"default-efficiency") || node.
hasKey(
"efficiencies")) {
317 "Reaction '{}' specifies efficiency parameters\n"
318 "but does not involve third body colliders.", equation());
322void Reaction::setRate(shared_ptr<ReactionRate> rate)
326 "Reaction rate for reaction '{}' must not be empty.", equation());
329 for (
const auto& [
id, callback] : m_setRateCallbacks) {
334void Reaction::registerSetRateCallback(
void*
id,
const function<
void()>& callback)
336 m_setRateCallbacks[id] = callback;
339void Reaction::removeSetRateCallback(
void*
id)
341 m_setRateCallbacks.erase(
id);
344string Reaction::reactantString()
const
346 std::ostringstream result;
347 for (
auto iter = reactants.begin(); iter != reactants.end(); ++iter) {
348 if (iter != reactants.begin()) {
351 if (iter->second != 1.0) {
352 result << iter->second <<
" ";
354 result << iter->first;
357 result << m_third_body->collider();
362string Reaction::productString()
const
364 std::ostringstream result;
365 for (
auto iter = products.begin(); iter != products.end(); ++iter) {
366 if (iter != products.begin()) {
369 if (iter->second != 1.0) {
370 result << iter->second <<
" ";
372 result << iter->first;
375 result << m_third_body->collider();
380string Reaction::equation()
const
383 return reactantString() +
" <=> " + productString();
385 return reactantString() +
" => " + productString();
389void Reaction::setEquation(
const string& equation,
const Kinetics* kin)
392 string rate_type = (m_rate) ? m_rate->type() : input.getString(
"type",
"");
393 if (ba::starts_with(rate_type,
"three-body")) {
395 m_explicit_type =
true;
396 }
else if (rate_type ==
"elementary") {
398 m_explicit_type =
true;
400 }
else if (kin && kin->
thermo(0).
nDim() != 3) {
403 }
else if (rate_type ==
"electron-collision-plasma") {
412 for (
const auto& [name, stoich] : reactants) {
414 if (products.count(name)) {
418 size_t generic = third_body ==
"M"
419 || third_body ==
"(+M)" || third_body ==
"(+ M)";
422 if (stoich > 1 && products[third_body] > 1) {
426 if (ba::starts_with(name,
"(+") && ba::ends_with(name,
")")) {
433 if (ba::starts_with(rate_type,
"three-body")) {
435 "Reactants for reaction '{}'\n"
436 "do not contain a valid third body collider", equation);
440 }
else if (countF == 1) {
443 }
else if (countM > 1) {
445 "Multiple generic third body colliders 'M' are not supported", equation);
446 }
else if (count > 1) {
452 }
else if (m_third_body) {
454 auto& effs = m_third_body->efficiencies;
455 if (effs.size() != 1 || !reactants.count(effs.begin()->first)) {
457 "Detected ambiguous third body colliders in reaction '{}'\n"
458 "ThirdBody object needs to specify a single species", equation);
460 third_body = effs.begin()->first;
461 m_third_body->explicit_3rd =
true;
462 }
else if (input.hasKey(
"efficiencies")) {
464 auto effs = input[
"efficiencies"].asMap<
double>();
465 if (effs.size() != 1 || !reactants.count(effs.begin()->first)) {
467 "Detected ambiguous third body colliders in reaction '{}'\n"
468 "Collision efficiencies need to specify single species", equation);
470 third_body = effs.begin()->first;
471 m_third_body = make_shared<ThirdBody>(third_body);
472 m_third_body->explicit_3rd =
true;
473 }
else if (input.hasKey(
"default-efficiency")) {
476 "Detected ambiguous third body colliders in reaction '{}'\n"
477 "Third-body definition requires specification of efficiencies",
479 }
else if (ba::starts_with(rate_type,
"three-body")) {
482 "Detected ambiguous third body colliders in reaction '{}'\n"
483 "A valid ThirdBody or collision efficiency definition is required",
489 }
else if (third_body !=
"M" && !ba::starts_with(rate_type,
"three-body")
490 && !ba::starts_with(third_body,
"(+"))
499 for (
const auto& [name, stoich] : reactants) {
500 if (trunc(stoich) != stoich) {
503 nreac +=
static_cast<size_t>(stoich);
507 for (
const auto& [name, stoich] : products) {
508 if (trunc(stoich) != stoich) {
511 nprod +=
static_cast<size_t>(stoich);
515 if (nreac != 3 && nprod != 3) {
521 string tName = m_third_body->name();
522 if (tName != third_body && third_body !=
"M" && tName !=
"M") {
524 "Detected incompatible third body colliders in reaction '{}'\n"
525 "ThirdBody definition does not match equation", equation);
527 m_third_body->setName(third_body);
529 m_third_body = make_shared<ThirdBody>(third_body);
533 auto reac = reactants.find(third_body);
534 if (reac == reactants.end()) {
536 "Logic error interpreting apparent third-body reaction");
538 if (trunc(reac->second) != 1) {
541 reactants.erase(reac);
545 auto prod = products.find(third_body);
546 if (prod == products.end()) {
548 "Logic error interpreting apparent third-body reaction");
550 if (trunc(prod->second) != 1) {
553 products.erase(prod);
557string Reaction::type()
const
560 throw CanteraError(
"Reaction::type",
"Empty Reaction does not have a type");
563 string rate_type = m_rate->type();
564 string sub_type = m_rate->subType();
565 if (sub_type !=
"") {
566 return rate_type +
"-" + sub_type;
570 return "three-body-" + rate_type;
591 for (
const auto& [name, order] : orders) {
594 rate_units.
update(phase.standardConcentrationUnits(), -order);
596 for (
const auto& [name, stoich] : reactants) {
599 if (name ==
"M" || ba::starts_with(name,
"(+")) {
603 }
else if (orders.find(name) == orders.end()) {
606 rate_units.
update(phase.standardConcentrationUnits(), -stoich);
610 if (m_third_body && m_third_body->mass_action) {
615 Reaction::rate_units = rate_units.
product();
619void updateUndeclared(vector<string>& undeclared,
622 for (
const auto& [name, stoich]: comp) {
624 undeclared.emplace_back(name);
629void Reaction::checkBalance(
const Kinetics& kin)
const
634 for (
const auto& [name, stoich] : products) {
637 for (
size_t m = 0; m < ph.
nElements(); m++) {
642 for (
const auto& [name, stoich] : reactants) {
645 for (
size_t m = 0; m < ph.
nElements(); m++) {
652 for (
const auto& [elem, balance] : balr) {
653 double scale = std::max(std::abs(balr[elem]), std::abs(balp[elem]));
654 double elemdiff = fabs(balp[elem] - balr[elem]);
655 if (elemdiff > 1e-4 *
scale) {
657 msg += fmt::format(
" {} {} {}\n",
658 elem, balr[elem], balp[elem]);
663 "The following reaction is unbalanced: {}\n"
664 " Element Reactants Products\n{}",
673 double reac_sites = 0.0;
674 double prod_sites = 0.0;
676 for (
const auto& [name, stoich] : reactants) {
679 reac_sites += stoich * surf.size(k);
682 for (
const auto& [name, stoich] : products) {
683 size_t k = surf.speciesIndex(name);
685 prod_sites += stoich * surf.size(k);
688 if (fabs(reac_sites - prod_sites) > 1e-5 * (reac_sites + prod_sites)) {
690 "Number of surface sites not balanced in reaction {}.\n"
691 "Reactant sites: {}\nProduct sites: {}",
692 equation(), reac_sites, prod_sites);
696bool Reaction::checkSpecies(
const Kinetics& kin)
const
699 vector<string> undeclared;
700 updateUndeclared(undeclared, reactants, kin);
701 updateUndeclared(undeclared, products, kin);
702 if (!undeclared.empty()) {
706 throw InputFileError(
"Reaction::checkSpecies", input,
"Reaction '{}'\n"
707 "contains undeclared species: '{}'",
708 equation(), ba::join(undeclared,
"', '"));
713 updateUndeclared(undeclared, orders, kin);
714 if (!undeclared.empty()) {
718 if (input.hasKey(
"orders")) {
721 "defines reaction orders for undeclared species: '{}'",
722 equation(), ba::join(undeclared,
"', '"));
725 throw InputFileError(
"Reaction::checkSpecies", input,
"Reaction '{}'\n"
726 "defines reaction orders for undeclared species: '{}'",
727 equation(), ba::join(undeclared,
"', '"));
732 return m_third_body->checkSpecies(*
this, kin);
740bool Reaction::usesElectrochemistry(
const Kinetics& kin)
const
743 vector<double> e_counter(kin.
nPhases(), 0.0);
746 for (
const auto& [name, stoich] : products) {
754 for (
const auto& [name, stoich] : reactants) {
762 for (
double delta_e : e_counter) {
763 if (std::abs(delta_e) > 1e-4) {
772ThirdBody::ThirdBody(
const string& third_body)
777void ThirdBody::setName(
const string& third_body)
779 string name = third_body;
780 if (ba::starts_with(third_body,
"(+ ")) {
782 name = third_body.substr(3, third_body.size() - 4);
783 }
else if (ba::starts_with(third_body,
"(+")) {
785 name = third_body.substr(2, third_body.size() - 3);
788 if (name == m_name) {
791 if (name ==
"M" && efficiencies.size() == 1) {
796 if (efficiencies.size()) {
798 "Conflicting efficiency definition for explicit third body '{}'", name);
801 default_efficiency = 0.;
802 efficiencies[m_name] = 1.;
805ThirdBody::ThirdBody(
const AnyMap& node)
810void ThirdBody::setParameters(
const AnyMap& node)
812 if (node.
hasKey(
"default-efficiency")) {
813 double value = node[
"default-efficiency"].asDouble();
814 if (m_name !=
"M" && value != 0.) {
815 throw InputFileError(
"ThirdBody::setParameters", node[
"default-efficiency"],
816 "Invalid default efficiency for explicit collider {};\n"
817 "value is optional and/or needs to be zero", m_name);
819 default_efficiency = value;
821 if (node.
hasKey(
"efficiencies")) {
822 efficiencies = node[
"efficiencies"].asMap<
double>();
825 && (efficiencies.size() != 1 || efficiencies.begin()->first != m_name))
828 "Detected incompatible third body colliders definitions");
832void ThirdBody::getParameters(
AnyMap& node)
const
834 if (m_name ==
"M" || explicit_3rd) {
835 if (efficiencies.size()) {
836 node[
"efficiencies"] = efficiencies;
839 if (default_efficiency != 1.0 && !explicit_3rd) {
840 node[
"default-efficiency"] = default_efficiency;
845double ThirdBody::efficiency(
const string& k)
const
847 return getValue(efficiencies, k, default_efficiency);
850string ThirdBody::collider()
const
853 return " + " + m_name;
855 return " (+" + m_name +
")";
860 vector<string> undeclared;
861 updateUndeclared(undeclared, efficiencies, kin);
863 if (!undeclared.empty()) {
867 rxn.
input[
"efficiencies"],
"Reaction '{}'\n"
868 "defines third-body efficiencies for undeclared species: '{}'",
869 rxn.
equation(), ba::join(undeclared,
"', '"));
873 "is a three-body reaction with undeclared species: '{}'",
874 rxn.
equation(), ba::join(undeclared,
"', '"));
887 return make_unique<Reaction>();
892 return make_unique<Reaction>(rxn_node, kin);
900 vector<string> tokens;
902 tokens.push_back(
"+");
904 size_t last_used =
npos;
905 bool reactants =
true;
906 for (
size_t i = 1; i < tokens.size(); i++) {
907 if (tokens[i] ==
"+" || ba::starts_with(tokens[i],
"(+") ||
908 tokens[i] ==
"<=>" || tokens[i] ==
"=" || tokens[i] ==
"=>") {
909 string species = tokens[i-1];
912 bool mass_action =
true;
913 if (last_used !=
npos && tokens[last_used] ==
"(+"
914 && ba::ends_with(species,
")")) {
917 species =
"(+" + species;
918 }
else if (last_used == i - 1 && ba::starts_with(species,
"(+")
919 && ba::ends_with(species,
")")) {
922 }
else if (last_used == i - 2) {
924 }
else if (last_used == i - 3) {
934 "Error parsing reaction string '{}'.\n"
935 "Current token: '{}'\nlast_used: '{}'",
937 (last_used ==
npos) ?
"n/a" : tokens[last_used]);
940 && mass_action && species !=
"M"))
955 if (tokens[i] ==
"<=>" || tokens[i] ==
"=") {
958 }
else if (tokens[i] ==
"=>") {
967 vector<shared_ptr<Reaction>> all_reactions;
969 auto R = make_shared<Reaction>(node, kinetics);
970 R->validate(kinetics);
971 if (R->valid() && R->checkSpecies(kinetics)) {
972 all_reactions.emplace_back(R);
975 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 exclude(const string &key)
Mark key as excluded from this map.
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 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.
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
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...