15 #include <unordered_set> 21 Kinetics::Kinetics() :
27 m_skipUndeclaredSpecies(false),
28 m_skipUndeclaredThirdBodies(false)
32 Kinetics::~Kinetics() {}
79 std::map<size_t, std::vector<size_t> > participants;
80 std::vector<std::map<int, double> > net_stoich;
81 std::unordered_set<size_t> unmatched_duplicates;
84 unmatched_duplicates.insert(i);
90 unsigned long int key = 0;
92 net_stoich.emplace_back();
93 std::map<int, double>& net = net_stoich.back();
97 net[-1 -k] -= sp.second;
102 net[1+k] += sp.second;
106 vector<size_t>& related = participants[key];
107 for (
size_t m : related) {
111 unmatched_duplicates.erase(i);
112 unmatched_duplicates.erase(m);
126 bool thirdBodyOk =
true;
142 bool thirdBodyOk =
true;
158 "Undeclared duplicate reactions detected:\n" 159 "Reaction {}: {}\nReaction {}: {}\n",
165 participants[key].push_back(i);
167 if (unmatched_duplicates.size()) {
168 size_t i = *unmatched_duplicates.begin();
171 "No duplicate found for declared duplicate reaction number {}" 181 std::map<int, double>& r2)
const 183 std::unordered_set<int> keys;
185 keys.insert(r.first);
188 keys.insert(r.first);
190 int k1 = r1.begin()->first;
192 doublereal ratio = 0.0;
193 if (r1[k1] && r2[k1]) {
194 ratio = r2[k1]/r1[k1];
195 bool different =
false;
197 if ((r1[k] && !r2[k]) ||
199 (r1[k] && fabs(r2[k]/r1[k] - ratio) > 1.e-8)) {
210 if (r1[k1] == 0.0 || r2[-k1] == 0.0) {
213 ratio = r2[-k1]/r1[k1];
215 if ((r1[k] && !r2[-k]) ||
216 (!r1[k] && r2[-k]) ||
217 (r1[k] && fabs(r2[-k]/r1[k] - ratio) > 1.e-8)) {
231 double stoich = sp.second;
232 for (
size_t m = 0; m < ph.
nElements(); m++) {
240 double stoich = sp.second;
241 for (
size_t m = 0; m < ph.
nElements(); m++) {
248 for (
const auto& el : balr) {
249 const string& elem = el.first;
250 double elemsum = balr[elem] + balp[elem];
251 double elemdiff = fabs(balp[elem] - balr[elem]);
252 if (elemsum > 0.0 && elemdiff/elemsum > 1e-4) {
254 msg += fmt::format(
" {} {} {}\n",
255 elem, balr[elem], balp[elem]);
259 msg =
"The following reaction is unbalanced: " + R.
equation() +
"\n" +
260 " Element Reactants Products\n" + msg;
266 doublereal* phase_data)
268 for (
size_t n = 0; n <
nPhases(); n++) {
272 data +
m_start[n] + nsp, phase_data);
276 throw CanteraError(
"Kinetics::selectPhase",
"Phase not found.");
281 for (
size_t n =
m_start.size()-1; n !=
npos; n--) {
291 for (
size_t n = 0; n <
m_thermo.size(); n++) {
303 const std::string& ph)
const 309 for (
size_t n = 0; n <
m_thermo.size(); n++) {
324 for (
size_t n = 0; n <
m_thermo.size(); n++) {
330 throw CanteraError(
"speciesPhase",
"unknown species "+nm);
335 for (
size_t n =
m_start.size()-1; n !=
npos; n--) {
340 throw CanteraError(
"speciesPhaseIndex",
"illegal species index: {}", k);
375 fill(deltaProp, deltaProp +
nReactions(), 0.0);
385 fill(deltaProp, deltaProp +
nReactions(), 0.0);
397 fill(cdot, cdot +
m_kk, 0.0);
411 fill(ddot, ddot +
m_kk, 0.0);
422 fill(net, net +
m_kk, 0.0);
454 for (
size_t i = 0; i <
m_thermo.size(); i++) {
473 if (r->reversible && !r->orders.empty()) {
474 throw CanteraError(
"Kinetics::addReaction",
"Reaction orders may only " 475 "be given for irreversible reactions");
479 for (
const auto& sp : r->reactants) {
484 throw CanteraError(
"Kinetics::addReaction",
"Reaction '" +
485 r->equation() +
"' contains the undeclared species '" +
490 for (
const auto& sp : r->products) {
495 throw CanteraError(
"Kinetics::addReaction",
"Reaction '" +
496 r->equation() +
"' contains the undeclared species '" +
501 for (
const auto& sp : r->orders) {
506 throw CanteraError(
"Kinetics::addReaction",
"Reaction '{}' has " 507 "a reaction order specified for the undeclared species '{}'",
508 r->equation(), sp.first);
517 std::vector<size_t> rk, pk;
523 for (
const auto& sp : r->reactants) {
525 rstoich.push_back(sp.second);
528 for (
const auto& sp : r->products) {
530 pstoich.push_back(sp.second);
537 for (
const auto& sp : r->orders) {
540 auto rloc = std::find(rk.begin(), rk.end(), k);
541 if (rloc != rk.end()) {
542 rorder[rloc - rk.begin()] = sp.second;
549 rstoich.push_back(0.0);
550 rorder.push_back(sp.second);
563 m_rfn.push_back(0.0);
576 if (rNew->reaction_type != rOld->reaction_type) {
578 "Reaction types are different: {} != {}.",
579 rOld->reaction_type, rNew->reaction_type);
582 if (rNew->reactants != rOld->reactants) {
584 "Reactants are different: '{}' != '{}'.",
585 rOld->reactantString(), rNew->reactantString());
588 if (rNew->products != rOld->products) {
590 "Products are different: '{}' != '{}'.",
591 rOld->productString(), rNew->productString());
virtual void modifyReaction(size_t i, shared_ptr< Reaction > rNew)
Modify the rate expression associated with a reaction.
std::string kineticsSpeciesName(size_t k) const
Return the name of the kth species in the kinetics manager.
StoichManagerN m_irrevProductStoich
Stoichiometry manager for the products of irreversible reactions.
StoichManagerN m_revProductStoich
Stoichiometry manager for the products of reversible reactions.
vector_fp m_ropr
Reverse rate-of-progress for each reaction.
size_t nElements() const
Number of elements.
std::vector< thermo_t * > m_thermo
m_thermo is a vector of pointers to ThermoPhase objects that are involved with this kinetics operator...
void checkSpeciesIndex(size_t k) const
Check that the specified species index is in range Throws an exception if k is greater than nSpecies(...
int reaction_type
Type of the reaction.
size_t speciesIndex(const std::string &name) const
Returns the index of a species named 'name' within the Phase object.
thermo_t & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism...
size_t m_kk
The number of species in all of the phases that participate in this kinetics mechanism.
std::vector< size_t > m_start
m_start is a vector of integers specifying the beginning position for the species vector for the n'th...
const size_t npos
index returned by functions to indicate "no position"
virtual void getNetRatesOfProgress(doublereal *netROP)
Net rates of progress.
void selectPhase(const doublereal *data, const thermo_t *phase, doublereal *phase_data)
virtual void resizeSpecies()
Resize arrays with sizes that depend on the total number of species.
virtual bool addReaction(shared_ptr< Reaction > r)
Add a single reaction to the mechanism.
A class for managing third-body efficiencies, including default values.
size_t nSpecies() const
Returns the number of species in the phase.
virtual void getRevReactionDelta(const doublereal *g, doublereal *dg)
Given an array of species properties 'g', return in array 'dg' the change in this quantity in the rev...
vector_fp m_ropnet
Net rate-of-progress for each reaction.
const int CHEMACT_RXN
A chemical activation reaction.
size_t nPhases() const
The number of phases participating in the reaction mechanism.
virtual void getDestructionRates(doublereal *ddot)
Species destruction rates [kmol/m^3/s or kmol/m^2/s].
virtual void getCreationRates(doublereal *cdot)
Species creation rates [kmol/m^3/s or kmol/m^2/s].
void checkSpeciesArraySize(size_t mm) const
Check that an array size is at least nSpecies() Throws an exception if kk is less than nSpecies()...
virtual void getNetProductionRates(doublereal *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
vector_fp m_ropf
Forward rate-of-progress for each reaction.
Base class for a phase with thermodynamic properties.
const int FALLOFF_RXN
The general form for a gas-phase association or dissociation reaction, with a pressure-dependent rate...
virtual void init()
Prepare the class for the addition of reactions, after all phases have been added.
virtual double productStoichCoeff(size_t k, size_t i) const
Stoichiometric coefficient of species k as a product in reaction i.
void checkReactionBalance(const Reaction &R)
Check that the specified reaction is balanced (same number of atoms for each element in the reactants...
bool reversible
True if the current reaction is reversible. False otherwise.
virtual double reactantStoichCoeff(size_t k, size_t i) const
Stoichiometric coefficient of species k as a reactant in reaction i.
bool duplicate
True if the current reaction is marked as duplicate.
size_t nTotalSpecies() const
The total number of species in all phases participating in the kinetics mechanism.
std::string speciesName(size_t k) const
Name of the species with index k.
std::map< std::string, doublereal > Composition
Map from string names to doubles.
vector_fp m_rfn
Forward rate constant for each reaction.
StoichManagerN m_reactantStoich
Stoichiometry manager for the reactants for each reaction.
void checkReactionArraySize(size_t ii) const
Check that an array size is at least nReactions() Throws an exception if ii is less than nReactions()...
thermo_t & speciesPhase(const std::string &nm)
This function looks up the name of a species and returns a reference to the ThermoPhase object of the...
virtual void getRevRatesOfProgress(doublereal *revROP)
Return the Reverse rates of progress of the reactions.
vector_fp m_perturb
Vector of perturbation factors for each reaction's rate of progress vector.
A reaction that is first-order in [M] at low pressure, like a third-body reaction, but zeroth-order in [M] as pressure increases.
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
virtual std::string kineticsType() const
Identifies the Kinetics manager type.
const U & getValue(const std::map< T, U > &m, const T &key, const U &default_val)
Const accessor for a value in a std::map.
Base class for exceptions thrown by Cantera classes.
size_t nReactions() const
Number of reactions in the reaction mechanism.
Composition reactants
Reactant species and stoichiometric coefficients.
Base class for kinetics managers and also contains the kineticsmgr module documentation (see Kinetics...
Intermediate class which stores data about a reaction and its rate parameterization so that it can be...
std::map< std::string, size_t > m_phaseindex
Mapping of the phase id, i.e., the id attribute in the XML phase element to the position of the phase...
void checkPhaseArraySize(size_t mm) const
Check that an array size is at least nPhases() Throws an exception if mm is less than nPhases()...
virtual std::pair< size_t, size_t > checkDuplicates(bool throw_err=true) const
Check for unmarked duplicate reactions and unmatched marked duplicates.
double efficiency(const std::string &k) const
Get the third-body efficiency for species k
const int THREE_BODY_RXN
A gas-phase reaction that requires a third-body collision partner.
vector_fp m_rkcn
Reciprocal of the equilibrium constant in concentration units.
virtual std::string type() const
String indicating the thermodynamic model implemented.
double checkDuplicateStoich(std::map< int, double > &r1, std::map< int, double > &r2) const
Check whether r1 and r2 represent duplicate stoichiometries This function returns a ratio if two reac...
std::vector< shared_ptr< Reaction > > m_reactions
Vector of Reaction objects represented by this Kinetics manager.
std::string id() const
Return the string id for the phase.
size_t m_mindim
number of spatial dimensions of lowest-dimensional phase.
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
size_t speciesPhaseIndex(size_t k)
This function takes as an argument the kineticsSpecies index (i.e., the list index in the list of spe...
A reaction with a non-reacting third body "M" that acts to add or remove energy from the reacting spe...
Contains declarations for string manipulation functions within Cantera.
size_t m_surfphase
Index in the list of phases of the one surface phase.
size_t nDim() const
Returns the number of spatial dimensions (1, 2, or 3)
shared_ptr< Reaction > reaction(size_t i)
Return the Reaction object for reaction i.
An array index is out of range.
Composition products
Product species and stoichiometric coefficients.
void checkReactionIndex(size_t m) const
Check that the specified reaction index is in range Throws an exception if i is greater than nReactio...
size_t m_rxnphase
Phase Index where reactions are assumed to be taking place.
Namespace for the Cantera kernel.
virtual void getFwdRatesOfProgress(doublereal *fwdROP)
Return the forward rates of progress of the reactions.
doublereal nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
virtual void addPhase(thermo_t &thermo)
Add a phase to the kinetics manager object.
bool m_skipUndeclaredSpecies
virtual void getReactionDelta(const doublereal *property, doublereal *deltaProperty)
Change in species properties.
std::string elementName(size_t m) const
Name of the element with index m.
void checkPhaseIndex(size_t m) const
Check that the specified phase index is in range Throws an exception if m is greater than nPhases() ...
std::string equation() const
The chemical equation for this reaction.