17 #include <boost/algorithm/string/predicate.hpp>
20 #include <boost/algorithm/string.hpp>
22 namespace 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;
69 Reaction::Reaction(
const string& equation,
70 shared_ptr<ReactionRate> rate_,
71 shared_ptr<ThirdBody> tbody_)
72 : m_third_body(tbody_)
74 setEquation(equation);
76 if (m_third_body && m_third_body->name() !=
"M") {
77 m_third_body->explicit_3rd =
true;
83 string rate_type = node.
getString(
"type",
"Arrhenius");
86 "Cannot instantiate Reaction with empty Kinetics object.");
89 setParameters(node, kin);
97 if (ba::starts_with(rate_type,
"three-body-")) {
99 rateNode[
"type"] = rate_type.substr(11, rate_type.size() - 11);
106 if (rateNode.
hasKey(
"rate-constant")) {
107 if (!ba::starts_with(rate_type,
"interface-")) {
108 rateNode[
"type"] =
"interface-" + rate_type;
110 }
else if (node.
hasKey(
"sticking-coefficient")) {
111 if (!ba::starts_with(rate_type,
"sticking-")) {
112 rateNode[
"type"] =
"sticking-" + rate_type;
116 "Unable to infer interface reaction type.");
123 void Reaction::check()
125 if (!allow_nonreactant_orders) {
126 for (
const auto& [name, order] : orders) {
127 if (reactants.find(name) == reactants.end()) {
129 "Reaction order specified for non-reactant species '{}'", name);
134 if (!allow_negative_orders) {
135 for (
const auto& [name, order] : orders) {
138 "Negative reaction order specified for species '{}'", name);
147 if (reversible && !orders.empty()) {
149 "Reaction orders may only be given for irreversible reactions");
155 m_rate->check(equation());
159 AnyMap Reaction::parameters(
bool withInput)
const
167 static bool reg = AnyMap::addOrderingRules(
"Reaction",
169 {
"head",
"equation"},
170 {
"tail",
"duplicate"},
172 {
"tail",
"negative-orders"},
173 {
"tail",
"nonreactant-orders"}
176 out[
"__type__"] =
"Reaction";
181 void Reaction::getParameters(
AnyMap& reactionNode)
const
185 "Serialization of empty Reaction object is not supported.");
188 reactionNode[
"equation"] = equation();
191 reactionNode[
"duplicate"] =
true;
194 reactionNode[
"orders"] = orders;
196 if (allow_negative_orders) {
197 reactionNode[
"negative-orders"] =
true;
199 if (allow_nonreactant_orders) {
200 reactionNode[
"nonreactant-orders"] =
true;
203 reactionNode.
update(m_rate->parameters());
206 string rtype = reactionNode[
"type"].asString();
207 if (rtype ==
"pressure-dependent-Arrhenius") {
209 }
else if (m_explicit_type && ba::ends_with(rtype,
"Arrhenius")) {
212 reactionNode[
"type"] =
"three-body";
214 reactionNode[
"type"] =
"elementary";
216 }
else if (ba::ends_with(rtype,
"Arrhenius")) {
217 reactionNode.
erase(
"type");
218 }
else if (m_explicit_type) {
219 reactionNode[
"type"] = type();
220 }
else if (ba::ends_with(rtype,
"Blowers-Masel")) {
221 reactionNode[
"type"] =
"Blowers-Masel";
225 m_third_body->getParameters(reactionNode);
233 "Cannot set reaction parameters from empty node.");
238 setEquation(node[
"equation"].asString(), &kin);
240 if (node.
hasKey(
"orders")) {
241 for (
const auto& [name, order] : node[
"orders"].asMap<double>()) {
242 orders[name] = order;
251 duplicate = node.
getBool(
"duplicate",
false);
252 allow_negative_orders = node.
getBool(
"negative-orders",
false);
253 allow_nonreactant_orders = node.
getBool(
"nonreactant-orders",
false);
256 m_third_body->setParameters(node);
257 if (m_third_body->name() ==
"M" && m_third_body->efficiencies.size() == 1) {
258 m_third_body->explicit_3rd =
true;
260 }
else if (node.
hasKey(
"default-efficiency") || node.
hasKey(
"efficiencies")) {
262 "Reaction '{}' specifies efficiency parameters\n"
263 "but does not involve third body colliders.", equation());
267 void Reaction::setRate(shared_ptr<ReactionRate> rate)
271 "Reaction rate for reaction '{}' must not be empty.", equation());
275 string rate_type = m_rate->type();
277 if (rate_type ==
"falloff" || rate_type ==
"chemically-activated") {
278 if (m_third_body->mass_action && !m_from_composition) {
280 "Third-body collider does not use '(+{})' notation.",
281 m_third_body->name());
283 m_third_body->mass_action =
false;
284 }
else if (rate_type ==
"Chebyshev") {
285 if (m_third_body->name() ==
"M") {
286 warn_deprecated(
"Chebyshev reaction equation", input,
"Specifying 'M' "
287 "in the reaction equation for Chebyshev reactions is deprecated.");
288 m_third_body.reset();
290 }
else if (rate_type ==
"pressure-dependent-Arrhenius") {
291 if (m_third_body->name() ==
"M") {
293 "Found superfluous '{}' in pressure-dependent-Arrhenius reaction.",
294 m_third_body->name());
298 if (rate_type ==
"falloff" || rate_type ==
"chemically-activated") {
299 if (!m_from_composition) {
301 "Reaction equation for falloff reaction '{}'\n does not "
302 "contain valid pressure-dependent third body", equation());
304 m_third_body = make_shared<ThirdBody>(
"(+M)");
309 string Reaction::reactantString()
const
311 std::ostringstream result;
312 for (
auto iter = reactants.begin(); iter != reactants.end(); ++iter) {
313 if (iter != reactants.begin()) {
316 if (iter->second != 1.0) {
317 result << iter->second <<
" ";
319 result << iter->first;
322 result << m_third_body->collider();
327 string Reaction::productString()
const
329 std::ostringstream result;
330 for (
auto iter = products.begin(); iter != products.end(); ++iter) {
331 if (iter != products.begin()) {
334 if (iter->second != 1.0) {
335 result << iter->second <<
" ";
337 result << iter->first;
340 result << m_third_body->collider();
345 string Reaction::equation()
const
348 return reactantString() +
" <=> " + productString();
350 return reactantString() +
" => " + productString();
354 void Reaction::setEquation(
const string& equation,
const Kinetics* kin)
358 string rate_type = input.getString(
"type",
"");
359 if (ba::starts_with(rate_type,
"three-body")) {
361 m_explicit_type =
true;
362 }
else if (rate_type ==
"elementary") {
364 m_explicit_type =
true;
366 }
else if (kin && kin->
thermo(0).
nDim() != 3) {
374 for (
const auto& [name, stoich] : reactants) {
376 if (products.count(name)) {
378 size_t generic = third_body ==
"M"
379 || third_body ==
"(+M)" || third_body ==
"(+ M)";
382 if (stoich > 1 && products[third_body] > 1) {
390 if (ba::starts_with(rate_type,
"three-body")) {
392 "Reactants for reaction '{}'\n"
393 "do not contain a valid third body collider", equation);
397 }
else if (countM > 1) {
399 "Multiple generic third body colliders 'M' are not supported", equation);
401 }
else if (count > 1) {
407 }
else if (m_third_body) {
409 auto& effs = m_third_body->efficiencies;
410 if (effs.size() != 1 || !reactants.count(effs.begin()->first)) {
412 "Detected ambiguous third body colliders in reaction '{}'\n"
413 "ThirdBody object needs to specify a single species", equation);
415 third_body = effs.begin()->first;
416 m_third_body->explicit_3rd =
true;
417 }
else if (input.hasKey(
"efficiencies")) {
419 auto effs = input[
"efficiencies"].asMap<
double>();
420 if (effs.size() != 1 || !reactants.count(effs.begin()->first)) {
422 "Detected ambiguous third body colliders in reaction '{}'\n"
423 "Collision efficiencies need to specify single species", equation);
425 third_body = effs.begin()->first;
426 m_third_body = make_shared<ThirdBody>(third_body);
427 m_third_body->explicit_3rd =
true;
428 }
else if (input.hasKey(
"default-efficiency")) {
431 "Detected ambiguous third body colliders in reaction '{}'\n"
432 "Third-body definition requires specification of efficiencies",
434 }
else if (ba::starts_with(rate_type,
"three-body")) {
437 "Detected ambiguous third body colliders in reaction '{}'\n"
438 "A valid ThirdBody or collision efficiency definition is required",
444 }
else if (third_body !=
"M" && !ba::starts_with(rate_type,
"three-body")
445 && !ba::starts_with(third_body,
"(+"))
454 for (
const auto& [name, stoich] : reactants) {
455 if (trunc(stoich) != stoich) {
458 nreac +=
static_cast<size_t>(stoich);
462 for (
const auto& [name, stoich] : products) {
463 if (trunc(stoich) != stoich) {
466 nprod +=
static_cast<size_t>(stoich);
470 if (nreac != 3 && nprod != 3) {
476 string tName = m_third_body->name();
477 if (tName != third_body && third_body !=
"M" && tName !=
"M") {
479 "Detected incompatible third body colliders in reaction '{}'\n"
480 "ThirdBody definition does not match equation", equation);
482 m_third_body->setName(third_body);
484 m_third_body = make_shared<ThirdBody>(third_body);
488 auto reac = reactants.find(third_body);
489 if (trunc(reac->second) != 1) {
492 reactants.erase(reac);
496 auto prod = products.find(third_body);
497 if (trunc(prod->second) != 1) {
500 products.erase(prod);
504 string Reaction::type()
const
507 throw CanteraError(
"Reaction::type",
"Empty Reaction does not have a type");
510 string rate_type = m_rate->type();
511 string sub_type = m_rate->subType();
512 if (sub_type !=
"") {
513 return rate_type +
"-" + sub_type;
517 return "three-body-" + rate_type;
538 for (
const auto& [name, order] : orders) {
541 rate_units.
update(phase.standardConcentrationUnits(), -order);
543 for (
const auto& [name, stoich] : reactants) {
546 if (name ==
"M" || ba::starts_with(name,
"(+")) {
550 }
else if (orders.find(name) == orders.end()) {
553 rate_units.
update(phase.standardConcentrationUnits(), -stoich);
557 if (m_third_body && m_third_body->mass_action) {
562 Reaction::rate_units = rate_units.
product();
566 void updateUndeclared(vector<string>& undeclared,
569 for (
const auto& [name, stoich]: comp) {
571 undeclared.emplace_back(name);
576 void Reaction::checkBalance(
const Kinetics& kin)
const
581 for (
const auto& [name, stoich] : products) {
584 for (
size_t m = 0; m < ph.
nElements(); m++) {
589 for (
const auto& [name, stoich] : reactants) {
592 for (
size_t m = 0; m < ph.
nElements(); m++) {
599 for (
const auto& [elem, balance] : balr) {
600 double elemsum = balr[elem] + balp[elem];
601 double elemdiff = fabs(balp[elem] - balr[elem]);
602 if (elemsum > 0.0 && elemdiff / elemsum > 1e-4) {
604 msg += fmt::format(
" {} {} {}\n",
605 elem, balr[elem], balp[elem]);
610 "The following reaction is unbalanced: {}\n"
611 " Element Reactants Products\n{}",
620 double reac_sites = 0.0;
621 double prod_sites = 0.0;
623 for (
const auto& [name, stoich] : reactants) {
626 reac_sites += stoich * surf.size(k);
629 for (
const auto& [name, stoich] : products) {
630 size_t k = surf.speciesIndex(name);
632 prod_sites += stoich * surf.size(k);
635 if (fabs(reac_sites - prod_sites) > 1e-5 * (reac_sites + prod_sites)) {
637 "Number of surface sites not balanced in reaction {}.\n"
638 "Reactant sites: {}\nProduct sites: {}",
639 equation(), reac_sites, prod_sites);
643 bool Reaction::checkSpecies(
const Kinetics& kin)
const
646 vector<string> undeclared;
647 updateUndeclared(undeclared, reactants, kin);
648 updateUndeclared(undeclared, products, kin);
649 if (!undeclared.empty()) {
653 throw InputFileError(
"Reaction::checkSpecies", input,
"Reaction '{}'\n"
654 "contains undeclared species: '{}'",
655 equation(), ba::join(undeclared,
"', '"));
660 updateUndeclared(undeclared, orders, kin);
661 if (!undeclared.empty()) {
665 if (input.hasKey(
"orders")) {
668 "defines reaction orders for undeclared species: '{}'",
669 equation(), ba::join(undeclared,
"', '"));
672 throw InputFileError(
"Reaction::checkSpecies", input,
"Reaction '{}'\n"
673 "defines reaction orders for undeclared species: '{}'",
674 equation(), ba::join(undeclared,
"', '"));
679 return m_third_body->checkSpecies(*
this, kin);
687 bool Reaction::usesElectrochemistry(
const Kinetics& kin)
const
690 vector<double> e_counter(kin.
nPhases(), 0.0);
693 for (
const auto& [name, stoich] : products) {
701 for (
const auto& [name, stoich] : reactants) {
709 for (
double delta_e : e_counter) {
710 if (std::abs(delta_e) > 1e-4) {
719 ThirdBody::ThirdBody(
const string& third_body)
724 void ThirdBody::setName(
const string& third_body)
726 string name = third_body;
727 if (ba::starts_with(third_body,
"(+ ")) {
729 name = third_body.substr(3, third_body.size() - 4);
730 }
else if (ba::starts_with(third_body,
"(+")) {
732 name = third_body.substr(2, third_body.size() - 3);
735 if (name == m_name) {
738 if (name ==
"M" && efficiencies.size() == 1) {
743 if (efficiencies.size()) {
745 "Conflicting efficiency definition for explicit third body '{}'", name);
748 default_efficiency = 0.;
749 efficiencies[m_name] = 1.;
752 ThirdBody::ThirdBody(
const AnyMap& node)
757 void ThirdBody::setParameters(
const AnyMap& node)
759 if (node.
hasKey(
"default-efficiency")) {
760 double value = node[
"default-efficiency"].asDouble();
761 if (m_name !=
"M" && value != 0.) {
762 throw InputFileError(
"ThirdBody::setParameters", node[
"default-efficiency"],
763 "Invalid default efficiency for explicit collider {};\n"
764 "value is optional and/or needs to be zero", m_name);
766 default_efficiency = value;
768 if (node.
hasKey(
"efficiencies")) {
769 efficiencies = node[
"efficiencies"].asMap<
double>();
772 && (efficiencies.size() != 1 || efficiencies.begin()->first != m_name))
775 "Detected incompatible third body colliders definitions");
779 void ThirdBody::getParameters(
AnyMap& node)
const
781 if (m_name ==
"M" || explicit_3rd) {
782 if (efficiencies.size()) {
783 node[
"efficiencies"] = efficiencies;
786 if (default_efficiency != 1.0 && !explicit_3rd) {
787 node[
"default-efficiency"] = default_efficiency;
792 double ThirdBody::efficiency(
const string& k)
const
794 return getValue(efficiencies, k, default_efficiency);
797 string ThirdBody::collider()
const
800 return " + " + m_name;
802 return " (+" + m_name +
")";
807 vector<string> undeclared;
808 updateUndeclared(undeclared, efficiencies, kin);
810 if (!undeclared.empty()) {
814 rxn.
input[
"efficiencies"],
"Reaction '{}'\n"
815 "defines third-body efficiencies for undeclared species: '{}'",
816 rxn.
equation(), ba::join(undeclared,
"', '"));
820 "is a three-body reaction with undeclared species: '{}'",
821 rxn.
equation(), ba::join(undeclared,
"', '"));
834 return make_unique<Reaction>();
839 return make_unique<Reaction>(rxn_node, kin);
847 vector<string> tokens;
849 tokens.push_back(
"+");
851 size_t last_used =
npos;
852 bool reactants =
true;
853 for (
size_t i = 1; i < tokens.size(); i++) {
854 if (tokens[i] ==
"+" || ba::starts_with(tokens[i],
"(+") ||
855 tokens[i] ==
"<=>" || tokens[i] ==
"=" || tokens[i] ==
"=>") {
856 string species = tokens[i-1];
859 bool mass_action =
true;
860 if (last_used !=
npos && tokens[last_used] ==
"(+"
861 && ba::ends_with(species,
")")) {
864 species =
"(+" + species;
865 }
else if (last_used == i - 1 && ba::starts_with(species,
"(+")
866 && ba::ends_with(species,
")")) {
869 }
else if (last_used == i - 2) {
871 }
else if (last_used == i - 3) {
881 "Error parsing reaction string '{}'.\n"
882 "Current token: '{}'\nlast_used: '{}'",
884 (last_used ==
npos) ?
"n/a" : tokens[last_used]);
887 && mass_action && species !=
"M"))
902 if (tokens[i] ==
"<=>" || tokens[i] ==
"=") {
905 }
else if (tokens[i] ==
"=>") {
914 vector<shared_ptr<Reaction>> all_reactions;
916 auto R = make_shared<Reaction>(node, kinetics);
918 R->validate(kinetics);
919 if (R->valid() && R->checkSpecies(kinetics)) {
920 all_reactions.emplace_back(R);
923 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.
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...
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
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.
unique_ptr< Reaction > newReaction(const AnyMap &rxn_node, const Kinetics &kin)
Create a new Reaction object using the specified parameters.
const U & getValue(const map< T, U > &m, const T &key, const U &default_val)
Const accessor for a value in a map.
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
map< string, double > Composition
Map from string names to doubles.
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...