1 /**
2 * @file Kinetics.cpp Declarations for the base class for kinetics managers
3 * (see \ref kineticsmgr and class \link Cantera::Kinetics Kinetics \endlink).
4 *
5 * Kinetics managers calculate rates of progress of species due to
6 * homogeneous or heterogeneous kinetics.
7 */
8
9 // This file is part of Cantera. See License.txt in the top-level directory or
10 // at https://cantera.org/license.txt for license and copyright information.
11
12 #include "cantera/kinetics/Kinetics.h"
13 #include "cantera/kinetics/KineticsFactory.h"
14 #include "cantera/kinetics/Reaction.h"
15 #include "cantera/thermo/ThermoPhase.h"
16 #include "cantera/base/stringUtils.h"
17 #include "cantera/base/utilities.h"
18 #include "cantera/base/global.h"
19 #include <unordered_set>
20 #include <boost/algorithm/string/join.hpp>
21
22 using namespace std;
23
24 namespace Cantera
25 {
Kinetics()26 Kinetics::Kinetics() :
27 m_ready(false),
28 m_kk(0),
29 m_thermo(0),
30 m_surfphase(npos),
31 m_rxnphase(npos),
32 m_mindim(4),
33 m_skipUndeclaredSpecies(false),
34 m_skipUndeclaredThirdBodies(false)
35 {
36 }
37
~Kinetics()38 Kinetics::~Kinetics() {}
39
checkReactionIndex(size_t i) const40 void Kinetics::checkReactionIndex(size_t i) const
41 {
42 if (i >= nReactions()) {
43 throw IndexError("Kinetics::checkReactionIndex", "reactions", i,
44 nReactions()-1);
45 }
46 }
47
resizeReactions()48 void Kinetics::resizeReactions()
49 {
50 size_t nRxn = nReactions();
51
52 // Stoichiometry managers
53 m_reactantStoich.resizeCoeffs(m_kk, nRxn);
54 m_productStoich.resizeCoeffs(m_kk, nRxn);
55 m_revProductStoich.resizeCoeffs(m_kk, nRxn);
56
57 m_ready = true;
58 }
59
checkReactionArraySize(size_t ii) const60 void Kinetics::checkReactionArraySize(size_t ii) const
61 {
62 if (nReactions() > ii) {
63 throw ArraySizeError("Kinetics::checkReactionArraySize", ii,
64 nReactions());
65 }
66 }
67
checkPhaseIndex(size_t m) const68 void Kinetics::checkPhaseIndex(size_t m) const
69 {
70 if (m >= nPhases()) {
71 throw IndexError("Kinetics::checkPhaseIndex", "phase", m, nPhases()-1);
72 }
73 }
74
checkPhaseArraySize(size_t mm) const75 void Kinetics::checkPhaseArraySize(size_t mm) const
76 {
77 if (nPhases() > mm) {
78 throw ArraySizeError("Kinetics::checkPhaseArraySize", mm, nPhases());
79 }
80 }
81
checkSpeciesIndex(size_t k) const82 void Kinetics::checkSpeciesIndex(size_t k) const
83 {
84 if (k >= m_kk) {
85 throw IndexError("Kinetics::checkSpeciesIndex", "species", k, m_kk-1);
86 }
87 }
88
checkSpeciesArraySize(size_t kk) const89 void Kinetics::checkSpeciesArraySize(size_t kk) const
90 {
91 if (m_kk > kk) {
92 throw ArraySizeError("Kinetics::checkSpeciesArraySize", kk, m_kk);
93 }
94 }
95
checkDuplicates(bool throw_err) const96 std::pair<size_t, size_t> Kinetics::checkDuplicates(bool throw_err) const
97 {
98 //! Map of (key indicating participating species) to reaction numbers
99 std::map<size_t, std::vector<size_t> > participants;
100 std::vector<std::map<int, double> > net_stoich;
101 std::unordered_set<size_t> unmatched_duplicates;
102 for (size_t i = 0; i < m_reactions.size(); i++) {
103 if (m_reactions[i]->duplicate) {
104 unmatched_duplicates.insert(i);
105 }
106 }
107
108 for (size_t i = 0; i < m_reactions.size(); i++) {
109 // Get data about this reaction
110 unsigned long int key = 0;
111 Reaction& R = *m_reactions[i];
112 net_stoich.emplace_back();
113 std::map<int, double>& net = net_stoich.back();
114 for (const auto& sp : R.reactants) {
115 int k = static_cast<int>(kineticsSpeciesIndex(sp.first));
116 key += k*(k+1);
117 net[-1 -k] -= sp.second;
118 }
119 for (const auto& sp : R.products) {
120 int k = static_cast<int>(kineticsSpeciesIndex(sp.first));
121 key += k*(k+1);
122 net[1+k] += sp.second;
123 }
124
125 // Compare this reaction to others with similar participants
126 vector<size_t>& related = participants[key];
127 for (size_t m : related) {
128 Reaction& other = *m_reactions[m];
129 if (R.duplicate && other.duplicate) {
130 // marked duplicates
131 unmatched_duplicates.erase(i);
132 unmatched_duplicates.erase(m);
133 continue;
134 } else if (R.type() != other.type()) {
135 continue; // different reaction types
136 }
137 doublereal c = checkDuplicateStoich(net_stoich[i], net_stoich[m]);
138 if (c == 0) {
139 continue; // stoichiometries differ (not by a multiple)
140 } else if (c < 0.0 && !R.reversible && !other.reversible) {
141 continue; // irreversible reactions in opposite directions
142 } else if (R.type() == "falloff" ||
143 R.type() == "chemically-activated") {
144 ThirdBody& tb1 = dynamic_cast<FalloffReaction&>(R).third_body;
145 ThirdBody& tb2 = dynamic_cast<FalloffReaction&>(other).third_body;
146 bool thirdBodyOk = true;
147 for (size_t k = 0; k < nTotalSpecies(); k++) {
148 string s = kineticsSpeciesName(k);
149 if (tb1.efficiency(s) * tb2.efficiency(s) != 0.0) {
150 // non-zero third body efficiencies for species `s` in
151 // both reactions
152 thirdBodyOk = false;
153 break;
154 }
155 }
156 if (thirdBodyOk) {
157 continue; // No overlap in third body efficiencies
158 }
159 } else if (R.type() == "three-body") {
160 ThirdBody& tb1 = *(dynamic_cast<ThreeBodyReaction3&>(R).thirdBody());
161 ThirdBody& tb2 = *(dynamic_cast<ThreeBodyReaction3&>(other).thirdBody());
162 bool thirdBodyOk = true;
163 for (size_t k = 0; k < nTotalSpecies(); k++) {
164 string s = kineticsSpeciesName(k);
165 if (tb1.efficiency(s) * tb2.efficiency(s) != 0.0) {
166 // non-zero third body efficiencies for species `s` in
167 // both reactions
168 thirdBodyOk = false;
169 break;
170 }
171 }
172 if (thirdBodyOk) {
173 continue; // No overlap in third body efficiencies
174 }
175 } else if (R.type() == "three-body-legacy") {
176 ThirdBody& tb1 = dynamic_cast<ThreeBodyReaction2&>(R).third_body;
177 ThirdBody& tb2 = dynamic_cast<ThreeBodyReaction2&>(other).third_body;
178 bool thirdBodyOk = true;
179 for (size_t k = 0; k < nTotalSpecies(); k++) {
180 string s = kineticsSpeciesName(k);
181 if (tb1.efficiency(s) * tb2.efficiency(s) != 0.0) {
182 // non-zero third body efficiencies for species `s` in
183 // both reactions
184 thirdBodyOk = false;
185 break;
186 }
187 }
188 if (thirdBodyOk) {
189 continue; // No overlap in third body efficiencies
190 }
191 }
192 if (throw_err) {
193 throw InputFileError("Kinetics::checkDuplicates",
194 R.input, other.input,
195 "Undeclared duplicate reactions detected:\n"
196 "Reaction {}: {}\nReaction {}: {}\n",
197 i+1, other.equation(), m+1, R.equation());
198 } else {
199 return {i,m};
200 }
201 }
202 participants[key].push_back(i);
203 }
204 if (unmatched_duplicates.size()) {
205 size_t i = *unmatched_duplicates.begin();
206 if (throw_err) {
207 throw InputFileError("Kinetics::checkDuplicates",
208 m_reactions[i]->input,
209 "No duplicate found for declared duplicate reaction number {}"
210 " ({})", i, m_reactions[i]->equation());
211 } else {
212 return {i, i};
213 }
214 }
215 return {npos, npos};
216 }
217
checkDuplicateStoich(std::map<int,double> & r1,std::map<int,double> & r2) const218 double Kinetics::checkDuplicateStoich(std::map<int, double>& r1,
219 std::map<int, double>& r2) const
220 {
221 std::unordered_set<int> keys; // species keys (k+1 or -k-1)
222 for (auto& r : r1) {
223 keys.insert(r.first);
224 }
225 for (auto& r : r2) {
226 keys.insert(r.first);
227 }
228 int k1 = r1.begin()->first;
229 // check for duplicate written in the same direction
230 doublereal ratio = 0.0;
231 if (r1[k1] && r2[k1]) {
232 ratio = r2[k1]/r1[k1];
233 bool different = false;
234 for (int k : keys) {
235 if ((r1[k] && !r2[k]) ||
236 (!r1[k] && r2[k]) ||
237 (r1[k] && fabs(r2[k]/r1[k] - ratio) > 1.e-8)) {
238 different = true;
239 break;
240 }
241 }
242 if (!different) {
243 return ratio;
244 }
245 }
246
247 // check for duplicate written in the reverse direction
248 if (r1[k1] == 0.0 || r2[-k1] == 0.0) {
249 return 0.0;
250 }
251 ratio = r2[-k1]/r1[k1];
252 for (int k : keys) {
253 if ((r1[k] && !r2[-k]) ||
254 (!r1[k] && r2[-k]) ||
255 (r1[k] && fabs(r2[-k]/r1[k] - ratio) > 1.e-8)) {
256 return 0.0;
257 }
258 }
259 return ratio;
260 }
261
checkReactionBalance(const Reaction & R)262 void Kinetics::checkReactionBalance(const Reaction& R)
263 {
264 R.checkBalance(*this);
265 warn_deprecated("Kinetics::checkReactionBalance",
266 "To be removed after Cantera 2.6. Replacable by Reaction::checkBalance.");
267 }
268
selectPhase(const double * data,const ThermoPhase * phase,double * phase_data)269 void Kinetics::selectPhase(const double* data, const ThermoPhase* phase,
270 double* phase_data)
271 {
272 for (size_t n = 0; n < nPhases(); n++) {
273 if (phase == m_thermo[n]) {
274 size_t nsp = phase->nSpecies();
275 copy(data + m_start[n],
276 data + m_start[n] + nsp, phase_data);
277 return;
278 }
279 }
280 throw CanteraError("Kinetics::selectPhase", "Phase not found.");
281 }
282
kineticsSpeciesName(size_t k) const283 string Kinetics::kineticsSpeciesName(size_t k) const
284 {
285 for (size_t n = m_start.size()-1; n != npos; n--) {
286 if (k >= m_start[n]) {
287 return thermo(n).speciesName(k - m_start[n]);
288 }
289 }
290 return "<unknown>";
291 }
292
kineticsSpeciesIndex(const std::string & nm) const293 size_t Kinetics::kineticsSpeciesIndex(const std::string& nm) const
294 {
295 for (size_t n = 0; n < m_thermo.size(); n++) {
296 // Check the ThermoPhase object for a match
297 size_t k = thermo(n).speciesIndex(nm);
298 if (k != npos) {
299 return k + m_start[n];
300 }
301 }
302 return npos;
303 }
304
kineticsSpeciesIndex(const std::string & nm,const std::string & ph) const305 size_t Kinetics::kineticsSpeciesIndex(const std::string& nm,
306 const std::string& ph) const
307 {
308 if (ph == "<any>") {
309 return kineticsSpeciesIndex(nm);
310 }
311
312 for (size_t n = 0; n < m_thermo.size(); n++) {
313 string id = thermo(n).name();
314 if (ph == id) {
315 size_t k = thermo(n).speciesIndex(nm);
316 if (k == npos) {
317 return npos;
318 }
319 return k + m_start[n];
320 }
321 }
322 return npos;
323 }
324
speciesPhase(const std::string & nm)325 ThermoPhase& Kinetics::speciesPhase(const std::string& nm)
326 {
327 for (size_t n = 0; n < m_thermo.size(); n++) {
328 size_t k = thermo(n).speciesIndex(nm);
329 if (k != npos) {
330 return thermo(n);
331 }
332 }
333 throw CanteraError("Kinetics::speciesPhase", "unknown species '{}'", nm);
334 }
335
speciesPhase(const std::string & nm) const336 const ThermoPhase& Kinetics::speciesPhase(const std::string& nm) const
337 {
338 for (const auto thermo : m_thermo) {
339 if (thermo->speciesIndex(nm) != npos) {
340 return *thermo;
341 }
342 }
343 throw CanteraError("Kinetics::speciesPhase", "unknown species '{}'", nm);
344 }
345
speciesPhaseIndex(size_t k) const346 size_t Kinetics::speciesPhaseIndex(size_t k) const
347 {
348 for (size_t n = m_start.size()-1; n != npos; n--) {
349 if (k >= m_start[n]) {
350 return n;
351 }
352 }
353 throw CanteraError("Kinetics::speciesPhaseIndex",
354 "illegal species index: {}", k);
355 }
356
reactantStoichCoeff(size_t kSpec,size_t irxn) const357 double Kinetics::reactantStoichCoeff(size_t kSpec, size_t irxn) const
358 {
359 return m_reactantStoich.stoichCoeffs().coeff(kSpec, irxn);
360 }
361
productStoichCoeff(size_t kSpec,size_t irxn) const362 double Kinetics::productStoichCoeff(size_t kSpec, size_t irxn) const
363 {
364 return m_productStoich.stoichCoeffs().coeff(kSpec, irxn);
365 }
366
reactionType(size_t i) const367 int Kinetics::reactionType(size_t i) const {
368 warn_deprecated("Kinetics::reactionType",
369 "To be changed after Cantera 2.6. "
370 "Return string instead of magic number; use "
371 "Kinetics::reactionTypeStr during transition.");
372 return m_reactions[i]->reaction_type;
373 }
374
reactionTypeStr(size_t i) const375 std::string Kinetics::reactionTypeStr(size_t i) const {
376 return m_reactions[i]->type();
377 }
378
reactionString(size_t i) const379 std::string Kinetics::reactionString(size_t i) const
380 {
381 return m_reactions[i]->equation();
382 }
383
384 //! Returns a string containing the reactants side of the reaction equation.
reactantString(size_t i) const385 std::string Kinetics::reactantString(size_t i) const
386 {
387 return m_reactions[i]->reactantString();
388 }
389
390 //! Returns a string containing the products side of the reaction equation.
productString(size_t i) const391 std::string Kinetics::productString(size_t i) const
392 {
393 return m_reactions[i]->productString();
394 }
395
getFwdRatesOfProgress(doublereal * fwdROP)396 void Kinetics::getFwdRatesOfProgress(doublereal* fwdROP)
397 {
398 updateROP();
399 std::copy(m_ropf.begin(), m_ropf.end(), fwdROP);
400 }
401
getRevRatesOfProgress(doublereal * revROP)402 void Kinetics::getRevRatesOfProgress(doublereal* revROP)
403 {
404 updateROP();
405 std::copy(m_ropr.begin(), m_ropr.end(), revROP);
406 }
407
getNetRatesOfProgress(doublereal * netROP)408 void Kinetics::getNetRatesOfProgress(doublereal* netROP)
409 {
410 updateROP();
411 std::copy(m_ropnet.begin(), m_ropnet.end(), netROP);
412 }
413
getReactionDelta(const double * prop,double * deltaProp)414 void Kinetics::getReactionDelta(const double* prop, double* deltaProp)
415 {
416 fill(deltaProp, deltaProp + nReactions(), 0.0);
417 // products add
418 m_productStoich.incrementReactions(prop, deltaProp);
419 // reactants subtract
420 m_reactantStoich.decrementReactions(prop, deltaProp);
421 }
422
getRevReactionDelta(const double * prop,double * deltaProp)423 void Kinetics::getRevReactionDelta(const double* prop, double* deltaProp)
424 {
425 fill(deltaProp, deltaProp + nReactions(), 0.0);
426 // products add
427 m_revProductStoich.incrementReactions(prop, deltaProp);
428 // reactants subtract
429 m_reactantStoich.decrementReactions(prop, deltaProp);
430 }
431
getCreationRates(double * cdot)432 void Kinetics::getCreationRates(double* cdot)
433 {
434 updateROP();
435
436 // zero out the output array
437 fill(cdot, cdot + m_kk, 0.0);
438
439 // the forward direction creates product species
440 m_productStoich.incrementSpecies(m_ropf.data(), cdot);
441
442 // the reverse direction creates reactant species
443 m_reactantStoich.incrementSpecies(m_ropr.data(), cdot);
444 }
445
getDestructionRates(doublereal * ddot)446 void Kinetics::getDestructionRates(doublereal* ddot)
447 {
448 updateROP();
449
450 fill(ddot, ddot + m_kk, 0.0);
451 // the reverse direction destroys products in reversible reactions
452 m_revProductStoich.incrementSpecies(m_ropr.data(), ddot);
453 // the forward direction destroys reactants
454 m_reactantStoich.incrementSpecies(m_ropf.data(), ddot);
455 }
456
getNetProductionRates(doublereal * net)457 void Kinetics::getNetProductionRates(doublereal* net)
458 {
459 updateROP();
460
461 fill(net, net + m_kk, 0.0);
462 // products are created for positive net rate of progress
463 m_productStoich.incrementSpecies(m_ropnet.data(), net);
464 // reactants are destroyed for positive net rate of progress
465 m_reactantStoich.decrementSpecies(m_ropnet.data(), net);
466 }
467
addPhase(ThermoPhase & thermo)468 void Kinetics::addPhase(ThermoPhase& thermo)
469 {
470 // the phase with lowest dimensionality is assumed to be the
471 // phase/interface at which reactions take place
472 if (thermo.nDim() <= m_mindim) {
473 m_mindim = thermo.nDim();
474 m_rxnphase = nPhases();
475 }
476
477 // there should only be one surface phase
478 if (thermo.type() == kineticsType()) {
479 m_surfphase = nPhases();
480 m_rxnphase = nPhases();
481 }
482 m_thermo.push_back(&thermo);
483 m_phaseindex[m_thermo.back()->name()] = nPhases();
484 resizeSpecies();
485 }
486
parameters()487 AnyMap Kinetics::parameters()
488 {
489 AnyMap out;
490 string name = KineticsFactory::factory()->canonicalize(kineticsType());
491 if (name != "none") {
492 out["kinetics"] = name;
493 if (nReactions() == 0) {
494 out["reactions"] = "none";
495 }
496 }
497 return out;
498 }
499
resizeSpecies()500 void Kinetics::resizeSpecies()
501 {
502 m_kk = 0;
503 m_start.resize(nPhases());
504
505 for (size_t i = 0; i < m_thermo.size(); i++) {
506 m_start[i] = m_kk; // global index of first species of phase i
507 m_kk += m_thermo[i]->nSpecies();
508 }
509 invalidateCache();
510 }
511
addReaction(shared_ptr<Reaction> r,bool resize)512 bool Kinetics::addReaction(shared_ptr<Reaction> r, bool resize)
513 {
514 r->validate();
515 if (m_kk == 0) {
516 init();
517 }
518 resizeSpecies();
519
520 // Check validity of reaction within the context of the Kinetics object
521 if (!r->checkSpecies(*this)) {
522 // Do not add reaction
523 return false;
524 }
525
526 // For reactions created outside the context of a Kinetics object, the units
527 // of the rate coefficient can't be determined in advance. Do that here.
528 if (r->rate_units.factor() == 0) {
529 r->calculateRateCoeffUnits(*this);
530 }
531
532 size_t irxn = nReactions(); // index of the new reaction
533
534 // indices of reactant and product species within this Kinetics object
535 std::vector<size_t> rk, pk;
536
537 // Reactant and product stoichiometric coefficients, such that rstoich[i] is
538 // the coefficient for species rk[i]
539 vector_fp rstoich, pstoich;
540
541 for (const auto& sp : r->reactants) {
542 rk.push_back(kineticsSpeciesIndex(sp.first));
543 rstoich.push_back(sp.second);
544 }
545
546 for (const auto& sp : r->products) {
547 pk.push_back(kineticsSpeciesIndex(sp.first));
548 pstoich.push_back(sp.second);
549 }
550
551 // The default order for each reactant is its stoichiometric coefficient,
552 // which can be overridden by entries in the Reaction.orders map. rorder[i]
553 // is the order for species rk[i].
554 vector_fp rorder = rstoich;
555 for (const auto& sp : r->orders) {
556 size_t k = kineticsSpeciesIndex(sp.first);
557 // Find the index of species k within rk
558 auto rloc = std::find(rk.begin(), rk.end(), k);
559 if (rloc != rk.end()) {
560 rorder[rloc - rk.begin()] = sp.second;
561 } else {
562 // If the reaction order involves a non-reactant species, add an
563 // extra term to the reactants with zero stoichiometry so that the
564 // stoichiometry manager can be used to compute the global forward
565 // reaction rate.
566 rk.push_back(k);
567 rstoich.push_back(0.0);
568 rorder.push_back(sp.second);
569 }
570 }
571
572 m_reactantStoich.add(irxn, rk, rorder, rstoich);
573 // product orders = product stoichiometric coefficients
574 m_productStoich.add(irxn, pk, pstoich, pstoich);
575 if (r->reversible) {
576 m_revProductStoich.add(irxn, pk, pstoich, pstoich);
577 }
578
579 m_reactions.push_back(r);
580 m_rfn.push_back(0.0);
581 m_rkcn.push_back(0.0);
582 m_ropf.push_back(0.0);
583 m_ropr.push_back(0.0);
584 m_ropnet.push_back(0.0);
585 m_perturb.push_back(1.0);
586 m_dH.push_back(0.0);
587
588 if (resize) {
589 resizeReactions();
590 } else {
591 m_ready = false;
592 }
593
594 return true;
595 }
596
modifyReaction(size_t i,shared_ptr<Reaction> rNew)597 void Kinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
598 {
599 checkReactionIndex(i);
600 shared_ptr<Reaction>& rOld = m_reactions[i];
601 if (rNew->type() != rOld->type()) {
602 throw CanteraError("Kinetics::modifyReaction",
603 "Reaction types are different: {} != {}.",
604 rOld->type(), rNew->type());
605 }
606
607 if (rNew->reactants != rOld->reactants) {
608 throw CanteraError("Kinetics::modifyReaction",
609 "Reactants are different: '{}' != '{}'.",
610 rOld->reactantString(), rNew->reactantString());
611 }
612
613 if (rNew->products != rOld->products) {
614 throw CanteraError("Kinetics::modifyReaction",
615 "Products are different: '{}' != '{}'.",
616 rOld->productString(), rNew->productString());
617 }
618 m_reactions[i] = rNew;
619 invalidateCache();
620 }
621
reaction(size_t i)622 shared_ptr<Reaction> Kinetics::reaction(size_t i)
623 {
624 checkReactionIndex(i);
625 return m_reactions[i];
626 }
627
reaction(size_t i) const628 shared_ptr<const Reaction> Kinetics::reaction(size_t i) const
629 {
630 checkReactionIndex(i);
631 return m_reactions[i];
632 }
633
634 }
635