Cantera  3.2.0a2
Loading...
Searching...
No Matches
Kinetics.cpp
Go to the documentation of this file.
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
18#include "cantera/base/global.h"
19#include <unordered_set>
20#include <boost/algorithm/string.hpp>
21
22using namespace std;
23
24namespace Cantera
25{
26
27void Kinetics::checkReactionIndex(size_t i) const
28{
29 if (i >= nReactions()) {
30 throw IndexError("Kinetics::checkReactionIndex", "reactions", i, nReactions());
31 }
32}
33
35{
36 size_t nRxn = nReactions();
37
38 // Stoichiometry managers
42
43 m_rbuf.resize(nRxn);
44
45 // products are created for positive net rate of progress
47 // reactants are destroyed for positive net rate of progress
49}
50
52{
53 if (nReactions() > ii) {
54 throw ArraySizeError("Kinetics::checkReactionArraySize", ii,
55 nReactions());
56 }
57}
58
59void Kinetics::checkPhaseIndex(size_t m) const
60{
61 if (m >= nPhases()) {
62 throw IndexError("Kinetics::checkPhaseIndex", "phase", m, nPhases());
63 }
64}
65
66void Kinetics::checkPhaseArraySize(size_t mm) const
67{
68 if (nPhases() > mm) {
69 throw ArraySizeError("Kinetics::checkPhaseArraySize", mm, nPhases());
70 }
71}
72
73shared_ptr<ThermoPhase> Kinetics::reactionPhase() const
74{
75 return m_thermo[0];
76}
77
78void Kinetics::checkSpeciesIndex(size_t k) const
79{
80 if (k >= m_kk) {
81 throw IndexError("Kinetics::checkSpeciesIndex", "species", k, m_kk);
82 }
83}
84
86{
87 if (m_kk > kk) {
88 throw ArraySizeError("Kinetics::checkSpeciesArraySize", kk, m_kk);
89 }
90}
91
93{
94 if (flag == "warn" || flag == "error" || flag == "mark-duplicate"
95 || flag == "modify-efficiency")
96 {
97 m_explicit_third_body_duplicates = flag;
98 } else {
99 throw CanteraError("Kinetics::setExplicitThirdBodyDuplicateHandling",
100 "Invalid flag '{}'", flag);
101 }
102}
103
104pair<size_t, size_t> Kinetics::checkDuplicates(bool throw_err, bool fix)
105{
106 //! Map of (key indicating participating species) to reaction numbers
107 map<size_t, vector<size_t>> participants;
108 vector<map<int, double>> net_stoich;
109 std::unordered_set<size_t> unmatched_duplicates;
110 for (size_t i = 0; i < m_reactions.size(); i++) {
111 if (m_reactions[i]->duplicate) {
112 unmatched_duplicates.insert(i);
113 }
114 }
115
116 vector<InputFileError> errs;
117 for (size_t i = 0; i < m_reactions.size(); i++) {
118 // Get data about this reaction
119 unsigned long int key = 0;
120 Reaction& R = *m_reactions[i];
121 net_stoich.emplace_back();
122 map<int, double>& net = net_stoich.back();
123 for (const auto& [name, stoich] : R.reactants) {
124 int k = static_cast<int>(kineticsSpeciesIndex(name));
125 key += k*(k+1);
126 net[-1 -k] -= stoich;
127 }
128 for (const auto& [name, stoich] : R.products) {
129 int k = static_cast<int>(kineticsSpeciesIndex(name));
130 key += k*(k+1);
131 net[1+k] += stoich;
132 }
133
134 // Compare this reaction to others with similar participants
135 vector<size_t>& related = participants[key];
136 for (size_t m : related) {
137 Reaction& other = *m_reactions[m];
138 if (R.duplicate && other.duplicate) {
139 // marked duplicates
140 unmatched_duplicates.erase(i);
141 unmatched_duplicates.erase(m);
142 continue;
143 } else if (R.type() != other.type()) {
144 continue; // different reaction types
145 } else if (R.rate() && other.rate()
146 && R.rate()->type() != other.rate()->type())
147 {
148 continue; // different rate parameterizations
149 }
150 double c = checkDuplicateStoich(net_stoich[i], net_stoich[m]);
151 if (c == 0) {
152 continue; // stoichiometries differ (not by a multiple)
153 } else if (c < 0.0 && !R.reversible && !other.reversible) {
154 continue; // irreversible reactions in opposite directions
155 } else if (R.usesThirdBody() && other.usesThirdBody()) {
156 ThirdBody& tb1 = *(R.thirdBody());
157 ThirdBody& tb2 = *(other.thirdBody());
158 bool thirdBodyOk = true;
159 for (size_t k = 0; k < nTotalSpecies(); k++) {
160 string s = kineticsSpeciesName(k);
161 if (tb1.efficiency(s) * tb2.efficiency(s) != 0.0) {
162 // non-zero third body efficiencies for species `s` in
163 // both reactions
164 thirdBodyOk = false;
165 break;
166 }
167 }
168 if (thirdBodyOk) {
169 continue; // No overlap in third body efficiencies
170 } else if ((tb1.name() == "M") != (tb2.name() == "M")) {
171 // Exactly one of the reactions uses M as the third body
172 if (m_explicit_third_body_duplicates == "mark-duplicate") {
173 R.duplicate = true;
174 other.duplicate = true;
175 continue;
176 } else if (m_explicit_third_body_duplicates == "modify-efficiency") {
177 if (tb1.name() == "M") {
178 tb1.efficiencies[tb2.name()] = 0.0;
179 } else {
180 tb2.efficiencies[tb1.name()] = 0.0;
181 }
182 continue;
183 } else if (m_explicit_third_body_duplicates == "warn") {
184 InputFileError msg("Kinetics::checkDuplicates",
185 R.input, other.input,
186 "Undeclared duplicate third body reactions with a common "
187 "third body detected.\nAdd the field "
188 "'explicit-third-body-duplicates: mark-duplicate' or\n"
189 "'explicit-third-body-duplicates: modify-efficiency' to "
190 "the YAML phase entry\nto choose how these reactions "
191 "should be handled and suppress this warning.\n"
192 "Reaction {}: {}\nReaction {}: {}\n",
193 m+1, R.equation(), i+1, other.equation());
194 if (!fix) {
195 warn_user("Kinetics::checkDuplicates", msg.what());
196 }
197 continue;
198 } // else m_explicit_third_body_duplicates == "error"
199 }
200 }
201 if (throw_err) {
202 errs.emplace_back("Kinetics::checkDuplicates",
203 R.input, other.input,
204 "Undeclared duplicate reactions detected:\n"
205 "Reaction {}: {}\nReaction {}: {}\n",
206 m+1, R.equation(), i+1, other.equation());
207 } else if (fix) {
208 R.duplicate = true;
209 other.duplicate = true;
210 } else {
211 return {i,m};
212 }
213 }
214 participants[key].push_back(i);
215 }
216 if (unmatched_duplicates.size()) {
217 for (auto i : unmatched_duplicates) {
218 if (throw_err) {
219 errs.emplace_back("Kinetics::checkDuplicates",
220 m_reactions[i]->input,
221 "No duplicate found for declared duplicate reaction number {}"
222 " ({})", i, m_reactions[i]->equation());
223 } else if (fix) {
224 m_reactions[i]->duplicate = false;
225 } else {
226 return {i, i};
227 }
228 }
229 }
230 if (errs.empty()) {
231 return {npos, npos};
232 } else if (errs.size() == 1) {
233 throw errs[0];
234 } else {
235 fmt::memory_buffer msg;
236 for (const auto& err : errs) {
237 fmt_append(msg, "\n{}\n", err.getMessage());
238 }
239 throw CanteraError("Kinetics::checkDuplicates", to_string(msg));
240 }
241}
242
243double Kinetics::checkDuplicateStoich(map<int, double>& r1, map<int, double>& r2) const
244{
245 std::unordered_set<int> keys; // species keys (k+1 or -k-1)
246 for (auto& [speciesKey, stoich] : r1) {
247 keys.insert(speciesKey);
248 }
249 for (auto& [speciesKey, stoich] : r2) {
250 keys.insert(speciesKey);
251 }
252 int k1 = r1.begin()->first;
253 // check for duplicate written in the same direction
254 double ratio = 0.0;
255 if (r1[k1] && r2[k1]) {
256 ratio = r2[k1]/r1[k1];
257 bool different = false;
258 for (int k : keys) {
259 if ((r1[k] && !r2[k]) ||
260 (!r1[k] && r2[k]) ||
261 (r1[k] && fabs(r2[k]/r1[k] - ratio) > 1.e-8)) {
262 different = true;
263 break;
264 }
265 }
266 if (!different) {
267 return ratio;
268 }
269 }
270
271 // check for duplicate written in the reverse direction
272 if (r1[k1] == 0.0 || r2[-k1] == 0.0) {
273 return 0.0;
274 }
275 ratio = r2[-k1]/r1[k1];
276 for (int k : keys) {
277 if ((r1[k] && !r2[-k]) ||
278 (!r1[k] && r2[-k]) ||
279 (r1[k] && fabs(r2[-k]/r1[k] - ratio) > 1.e-8)) {
280 return 0.0;
281 }
282 }
283 return ratio;
284}
285
286string Kinetics::kineticsSpeciesName(size_t k) const
287{
288 for (size_t n = m_start.size()-1; n != npos; n--) {
289 if (k >= m_start[n]) {
290 return thermo(n).speciesName(k - m_start[n]);
291 }
292 }
293 return "<unknown>";
294}
295
296size_t Kinetics::kineticsSpeciesIndex(const string& nm) const
297{
298 for (size_t n = 0; n < m_thermo.size(); n++) {
299 // Check the ThermoPhase object for a match
300 size_t k = thermo(n).speciesIndex(nm);
301 if (k != npos) {
302 return k + m_start[n];
303 }
304 }
305 return npos;
306}
307
309{
310 for (size_t n = 0; n < m_thermo.size(); n++) {
311 size_t k = thermo(n).speciesIndex(nm);
312 if (k != npos) {
313 return thermo(n);
314 }
315 }
316 throw CanteraError("Kinetics::speciesPhase", "unknown species '{}'", nm);
317}
318
319const ThermoPhase& Kinetics::speciesPhase(const string& nm) const
320{
321 for (const auto& thermo : m_thermo) {
322 if (thermo->speciesIndex(nm) != npos) {
323 return *thermo;
324 }
325 }
326 throw CanteraError("Kinetics::speciesPhase", "unknown species '{}'", nm);
327}
328
329size_t Kinetics::speciesPhaseIndex(size_t k) const
330{
331 for (size_t n = m_start.size()-1; n != npos; n--) {
332 if (k >= m_start[n]) {
333 return n;
334 }
335 }
336 throw CanteraError("Kinetics::speciesPhaseIndex",
337 "illegal species index: {}", k);
338}
339
340double Kinetics::reactantStoichCoeff(size_t kSpec, size_t irxn) const
341{
342 return m_reactantStoich.stoichCoeffs().coeff(kSpec, irxn);
343}
344
345double Kinetics::productStoichCoeff(size_t kSpec, size_t irxn) const
346{
347 return m_productStoich.stoichCoeffs().coeff(kSpec, irxn);
348}
349
351{
352 updateROP();
353 std::copy(m_ropf.begin(), m_ropf.end(), fwdROP);
354}
355
357{
358 updateROP();
359 std::copy(m_ropr.begin(), m_ropr.end(), revROP);
360}
361
363{
364 updateROP();
365 std::copy(m_ropnet.begin(), m_ropnet.end(), netROP);
366}
367
368void Kinetics::getReactionDelta(const double* prop, double* deltaProp) const
369{
370 fill(deltaProp, deltaProp + nReactions(), 0.0);
371 // products add
372 m_productStoich.incrementReactions(prop, deltaProp);
373 // reactants subtract
374 m_reactantStoich.decrementReactions(prop, deltaProp);
375}
376
377void Kinetics::getRevReactionDelta(const double* prop, double* deltaProp) const
378{
379 fill(deltaProp, deltaProp + nReactions(), 0.0);
380 // products add
381 m_revProductStoich.incrementReactions(prop, deltaProp);
382 // reactants subtract
383 m_reactantStoich.decrementReactions(prop, deltaProp);
384}
385
387{
388 updateROP();
389
390 // zero out the output array
391 fill(cdot, cdot + m_kk, 0.0);
392
393 // the forward direction creates product species
394 m_productStoich.incrementSpecies(m_ropf.data(), cdot);
395
396 // the reverse direction creates reactant species
397 m_reactantStoich.incrementSpecies(m_ropr.data(), cdot);
398}
399
401{
402 updateROP();
403
404 fill(ddot, ddot + m_kk, 0.0);
405 // the reverse direction destroys products in reversible reactions
406 m_revProductStoich.incrementSpecies(m_ropr.data(), ddot);
407 // the forward direction destroys reactants
408 m_reactantStoich.incrementSpecies(m_ropf.data(), ddot);
409}
410
412{
413 updateROP();
414
415 fill(net, net + m_kk, 0.0);
416 // products are created for positive net rate of progress
417 m_productStoich.incrementSpecies(m_ropnet.data(), net);
418 // reactants are destroyed for positive net rate of progress
419 m_reactantStoich.decrementSpecies(m_ropnet.data(), net);
420}
421
423{
424 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
425 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
426 // the forward direction creates product species
427 getFwdRatesOfProgress_ddT(buf.data());
428 out = m_productStoich.stoichCoeffs() * buf;
429 // the reverse direction creates reactant species
430 getRevRatesOfProgress_ddT(buf.data());
431 out += m_reactantStoich.stoichCoeffs() * buf;
432}
433
435{
436 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
437 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
438 // the forward direction creates product species
439 getFwdRatesOfProgress_ddP(buf.data());
440 out = m_productStoich.stoichCoeffs() * buf;
441 // the reverse direction creates reactant species
442 getRevRatesOfProgress_ddP(buf.data());
443 out += m_reactantStoich.stoichCoeffs() * buf;
444}
445
447{
448 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
449 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
450 // the forward direction creates product species
451 getFwdRatesOfProgress_ddC(buf.data());
452 out = m_productStoich.stoichCoeffs() * buf;
453 // the reverse direction creates reactant species
454 getRevRatesOfProgress_ddC(buf.data());
455 out += m_reactantStoich.stoichCoeffs() * buf;
456}
457
458Eigen::SparseMatrix<double> Kinetics::creationRates_ddX()
459{
460 Eigen::SparseMatrix<double> jac;
461 // the forward direction creates product species
463 // the reverse direction creates reactant species
465 return jac;
466}
467
468Eigen::SparseMatrix<double> Kinetics::creationRates_ddCi()
469{
470 Eigen::SparseMatrix<double> jac;
471 // the forward direction creates product species
473 // the reverse direction creates reactant species
475 return jac;
476}
477
479{
480 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
481 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
482 // the reverse direction destroys products in reversible reactions
483 getRevRatesOfProgress_ddT(buf.data());
484 out = m_revProductStoich.stoichCoeffs() * buf;
485 // the forward direction destroys reactants
486 getFwdRatesOfProgress_ddT(buf.data());
487 out += m_reactantStoich.stoichCoeffs() * buf;
488}
489
491{
492 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
493 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
494 // the reverse direction destroys products in reversible reactions
495 getRevRatesOfProgress_ddP(buf.data());
496 out = m_revProductStoich.stoichCoeffs() * buf;
497 // the forward direction destroys reactants
498 getFwdRatesOfProgress_ddP(buf.data());
499 out += m_reactantStoich.stoichCoeffs() * buf;
500}
501
503{
504 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
505 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
506 // the reverse direction destroys products in reversible reactions
507 getRevRatesOfProgress_ddC(buf.data());
508 out = m_revProductStoich.stoichCoeffs() * buf;
509 // the forward direction destroys reactants
510 getFwdRatesOfProgress_ddC(buf.data());
511 out += m_reactantStoich.stoichCoeffs() * buf;
512}
513
514Eigen::SparseMatrix<double> Kinetics::destructionRates_ddX()
515{
516 Eigen::SparseMatrix<double> jac;
517 // the reverse direction destroys products in reversible reactions
519 // the forward direction destroys reactants
521 return jac;
522}
523
524Eigen::SparseMatrix<double> Kinetics::destructionRates_ddCi()
525{
526 Eigen::SparseMatrix<double> jac;
527 // the reverse direction destroys products in reversible reactions
529 // the forward direction destroys reactants
531 return jac;
532}
533
535{
536 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
537 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
538 getNetRatesOfProgress_ddT(buf.data());
539 out = m_stoichMatrix * buf;
540}
541
543{
544 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
545 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
546 getNetRatesOfProgress_ddP(buf.data());
547 out = m_stoichMatrix * buf;
548}
549
551{
552 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
553 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
554 getNetRatesOfProgress_ddC(buf.data());
555 out = m_stoichMatrix * buf;
556}
557
558Eigen::SparseMatrix<double> Kinetics::netProductionRates_ddX()
559{
561}
562
563Eigen::SparseMatrix<double> Kinetics::netProductionRates_ddCi()
564{
566}
567
568void Kinetics::addThermo(shared_ptr<ThermoPhase> thermo)
569{
570 // the phase with lowest dimensionality is assumed to be the
571 // phase/interface at which reactions take place
572 if (thermo->nDim() <= m_mindim) {
573 if (!m_thermo.empty()) {
574 throw CanteraError("Kinetics::addThermo",
575 "The reacting (lowest dimensional) phase must be added first.");
576 }
577 m_mindim = thermo->nDim();
578 }
579
580 m_thermo.push_back(thermo);
581 m_phaseindex[m_thermo.back()->name()] = nPhases();
583}
584
586{
587 AnyMap out;
588 string name = KineticsFactory::factory()->canonicalize(kineticsType());
589 if (name != "none") {
590 out["kinetics"] = name;
591 if (nReactions() == 0) {
592 out["reactions"] = "none";
593 }
595 out["skip-undeclared-third-bodies"] = true;
596 }
597 if (m_explicit_third_body_duplicates == "error") {
598 // "warn" is the default, and does not need to be added. "mark-duplicate"
599 // and "modify-efficiency" do not need to be propagated here as their
600 // effects are already applied to the corresponding reactions.
601 out["explicit-third-body-duplicates"] = "error";
602 }
603 }
604 return out;
605}
606
608{
609 m_kk = 0;
610 m_start.resize(nPhases());
611
612 for (size_t i = 0; i < m_thermo.size(); i++) {
613 m_start[i] = m_kk; // global index of first species of phase i
614 m_kk += m_thermo[i]->nSpecies();
615 }
616 invalidateCache();
617}
618
619bool Kinetics::addReaction(shared_ptr<Reaction> r, bool resize)
620{
621 r->check();
622 r->validate(*this);
623
624 if (m_kk == 0) {
625 init();
626 }
628
629 // Check validity of reaction within the context of the Kinetics object
630 if (!r->checkSpecies(*this)) {
631 // Do not add reaction
632 return false;
633 }
634
635 // For reactions created outside the context of a Kinetics object, the units
636 // of the rate coefficient can't be determined in advance. Do that here.
637 if (r->rate_units.factor() == 0) {
638 r->rate()->setRateUnits(r->calculateRateCoeffUnits(*this));
639 }
640
641 size_t irxn = nReactions(); // index of the new reaction
642
643 // indices of reactant and product species within this Kinetics object
644 vector<size_t> rk, pk;
645
646 // Reactant and product stoichiometric coefficients, such that rstoich[i] is
647 // the coefficient for species rk[i]
648 vector<double> rstoich, pstoich;
649
650 for (const auto& [name, stoich] : r->reactants) {
651 rk.push_back(kineticsSpeciesIndex(name));
652 rstoich.push_back(stoich);
653 }
654
655 for (const auto& [name, stoich] : r->products) {
656 pk.push_back(kineticsSpeciesIndex(name));
657 pstoich.push_back(stoich);
658 }
659
660 // The default order for each reactant is its stoichiometric coefficient,
661 // which can be overridden by entries in the Reaction.orders map. rorder[i]
662 // is the order for species rk[i].
663 vector<double> rorder = rstoich;
664 for (const auto& [name, order] : r->orders) {
665 size_t k = kineticsSpeciesIndex(name);
666 // Find the index of species k within rk
667 auto rloc = std::find(rk.begin(), rk.end(), k);
668 if (rloc != rk.end()) {
669 rorder[rloc - rk.begin()] = order;
670 } else {
671 // If the reaction order involves a non-reactant species, add an
672 // extra term to the reactants with zero stoichiometry so that the
673 // stoichiometry manager can be used to compute the global forward
674 // reaction rate.
675 rk.push_back(k);
676 rstoich.push_back(0.0);
677 rorder.push_back(order);
678 }
679 }
680
681 m_reactantStoich.add(irxn, rk, rorder, rstoich);
682 // product orders = product stoichiometric coefficients
683 m_productStoich.add(irxn, pk, pstoich, pstoich);
684 if (r->reversible) {
685 m_revindex.push_back(irxn);
686 m_revProductStoich.add(irxn, pk, pstoich, pstoich);
687 } else {
688 m_irrev.push_back(irxn);
689 }
690
691 m_reactions.push_back(r);
692 m_rfn.push_back(0.0);
693 m_delta_gibbs0.push_back(0.0);
694 m_rkcn.push_back(0.0);
695 m_ropf.push_back(0.0);
696 m_ropr.push_back(0.0);
697 m_ropnet.push_back(0.0);
698 m_perturb.push_back(1.0);
699 m_dH.push_back(0.0);
700
701 if (resize) {
703 }
704
705 for (const auto& [id, callback] : m_reactionAddedCallbacks) {
706 callback();
707 }
708
709 return true;
710}
711
712void Kinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
713{
715 shared_ptr<Reaction>& rOld = m_reactions[i];
716
717 if (rNew->rate()->type() == "electron-collision-plasma") {
718 throw CanteraError("Kinetics::modifyReaction",
719 "Type electron-collision-plasma is not supported. "
720 "Use the rate object of the reaction to modify the data.");
721 }
722
723 if (rNew->type() != rOld->type()) {
724 throw CanteraError("Kinetics::modifyReaction",
725 "Reaction types are different: {} != {}.",
726 rOld->type(), rNew->type());
727 }
728
729 if (rNew->rate()->type() != rOld->rate()->type()) {
730 throw CanteraError("Kinetics::modifyReaction",
731 "ReactionRate types are different: {} != {}.",
732 rOld->rate()->type(), rNew->rate()->type());
733 }
734
735 if (rNew->reactants != rOld->reactants) {
736 throw CanteraError("Kinetics::modifyReaction",
737 "Reactants are different: '{}' != '{}'.",
738 rOld->reactantString(), rNew->reactantString());
739 }
740
741 if (rNew->products != rOld->products) {
742 throw CanteraError("Kinetics::modifyReaction",
743 "Products are different: '{}' != '{}'.",
744 rOld->productString(), rNew->productString());
745 }
746 m_reactions[i] = rNew;
747 invalidateCache();
748}
749
750shared_ptr<Reaction> Kinetics::reaction(size_t i)
751{
753 return m_reactions[i];
754}
755
756shared_ptr<const Reaction> Kinetics::reaction(size_t i) const
757{
759 return m_reactions[i];
760}
761
762}
Base class for kinetics managers and also contains the kineticsmgr module documentation (see Kinetics...
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:431
Array size error.
Base class for exceptions thrown by Cantera classes.
const char * what() const override
Get a description of the error.
An array index is out of range.
Error thrown for problems processing information contained in an AnyMap or AnyValue.
Definition AnyMap.h:749
virtual void resizeReactions()
Finalize Kinetics object and associated objects.
Definition Kinetics.cpp:34
void checkPhaseIndex(size_t m) const
Check that the specified phase index is in range Throws an exception if m is greater than nPhases()
Definition Kinetics.cpp:59
void checkSpeciesArraySize(size_t mm) const
Check that an array size is at least nSpecies() Throws an exception if kk is less than nSpecies().
Definition Kinetics.cpp:85
virtual void getFwdRatesOfProgress(double *fwdROP)
Return the forward rates of progress of the reactions.
Definition Kinetics.cpp:350
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
Definition Kinetics.h:240
vector< shared_ptr< Reaction > > m_reactions
Vector of Reaction objects represented by this Kinetics manager.
Definition Kinetics.h:1486
void checkSpeciesIndex(size_t k) const
Check that the specified species index is in range Throws an exception if k is greater than nSpecies(...
Definition Kinetics.cpp:78
double checkDuplicateStoich(map< int, double > &r1, map< int, double > &r2) const
Check whether r1 and r2 represent duplicate stoichiometries This function returns a ratio if two reac...
Definition Kinetics.cpp:243
vector< double > m_perturb
Vector of perturbation factors for each reaction's rate of progress vector.
Definition Kinetics.h:1483
vector< double > m_ropf
Forward rate-of-progress for each reaction.
Definition Kinetics.h:1527
shared_ptr< Reaction > reaction(size_t i)
Return the Reaction object for reaction i.
Definition Kinetics.cpp:750
virtual void getRevReactionDelta(const double *g, double *dg) const
Given an array of species properties 'g', return in array 'dg' the change in this quantity in the rev...
Definition Kinetics.cpp:377
void setExplicitThirdBodyDuplicateHandling(const string &flag)
Specify how to handle duplicate third body reactions where one reaction has an explicit third body an...
Definition Kinetics.cpp:92
virtual void getNetRatesOfProgress(double *netROP)
Net rates of progress.
Definition Kinetics.cpp:362
vector< size_t > m_start
m_start is a vector of integers specifying the beginning position for the species vector for the n'th...
Definition Kinetics.h:1504
virtual string kineticsType() const
Identifies the Kinetics manager type.
Definition Kinetics.h:145
vector< double > m_rkcn
Reciprocal of the equilibrium constant in concentration units.
Definition Kinetics.h:1524
shared_ptr< ThermoPhase > reactionPhase() const
Return pointer to phase where the reactions occur.
Definition Kinetics.cpp:73
size_t m_kk
The number of species in all of the phases that participate in this kinetics mechanism.
Definition Kinetics.h:1479
virtual pair< size_t, size_t > checkDuplicates(bool throw_err=true, bool fix=false)
Check for unmarked duplicate reactions and unmatched marked duplicates.
Definition Kinetics.cpp:104
virtual void getReactionDelta(const double *property, double *deltaProperty) const
Change in species properties.
Definition Kinetics.cpp:368
void checkPhaseArraySize(size_t mm) const
Check that an array size is at least nPhases() Throws an exception if mm is less than nPhases().
Definition Kinetics.cpp:66
vector< double > m_dH
The enthalpy change for each reaction to calculate Blowers-Masel rates.
Definition Kinetics.h:1539
vector< double > m_ropr
Reverse rate-of-progress for each reaction.
Definition Kinetics.h:1530
size_t nPhases() const
The number of phases participating in the reaction mechanism.
Definition Kinetics.h:185
vector< shared_ptr< ThermoPhase > > m_thermo
m_thermo is a vector of pointers to ThermoPhase objects that are involved with this kinetics operator
Definition Kinetics.h:1498
map< void *, function< void()> > m_reactionAddedCallbacks
Callback functions that are invoked when the reaction is added.
Definition Kinetics.h:1560
virtual bool addReaction(shared_ptr< Reaction > r, bool resize=true)
Add a single reaction to the mechanism.
Definition Kinetics.cpp:619
vector< size_t > m_revindex
Indices of reversible reactions.
Definition Kinetics.h:1535
string kineticsSpeciesName(size_t k) const
Return the name of the kth species in the kinetics manager.
Definition Kinetics.cpp:286
virtual void getDestructionRates(double *ddot)
Species destruction rates [kmol/m^3/s or kmol/m^2/s].
Definition Kinetics.cpp:400
AnyMap parameters()
Return the parameters for a phase definition which are needed to reconstruct an identical object usin...
Definition Kinetics.cpp:585
virtual void init()
Prepare the class for the addition of reactions, after all phases have been added.
Definition Kinetics.h:1278
virtual void modifyReaction(size_t i, shared_ptr< Reaction > rNew)
Modify the rate expression associated with a reaction.
Definition Kinetics.cpp:712
vector< double > m_ropnet
Net rate-of-progress for each reaction.
Definition Kinetics.h:1533
virtual double productStoichCoeff(size_t k, size_t i) const
Stoichiometric coefficient of species k as a product in reaction i.
Definition Kinetics.cpp:345
virtual void addThermo(shared_ptr< ThermoPhase > thermo)
Add a phase to the kinetics manager object.
Definition Kinetics.cpp:568
vector< double > m_rbuf
Buffer used for storage of intermediate reaction-specific results.
Definition Kinetics.h:1542
Eigen::SparseMatrix< double > m_stoichMatrix
Net stoichiometry (products - reactants)
Definition Kinetics.h:1474
map< string, size_t > m_phaseindex
Mapping of the phase name to the position of the phase within the kinetics object.
Definition Kinetics.h:1512
size_t m_mindim
number of spatial dimensions of lowest-dimensional phase.
Definition Kinetics.h:1515
StoichManagerN m_productStoich
Stoichiometry manager for the products for each reaction.
Definition Kinetics.h:1468
StoichManagerN m_revProductStoich
Stoichiometry manager for the products of reversible reactions.
Definition Kinetics.h:1471
size_t nReactions() const
Number of reactions in the reaction mechanism.
Definition Kinetics.h:153
virtual void getRevRatesOfProgress(double *revROP)
Return the Reverse rates of progress of the reactions.
Definition Kinetics.cpp:356
vector< size_t > m_irrev
Indices of irreversible reactions.
Definition Kinetics.h:1536
bool m_hasUndeclaredThirdBodies
Flag indicating whether reactions include undeclared third bodies.
Definition Kinetics.h:1551
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition Kinetics.h:274
vector< double > m_rfn
Forward rate constant for each reaction.
Definition Kinetics.h:1518
StoichManagerN m_reactantStoich
Stoichiometry manager for the reactants for each reaction.
Definition Kinetics.h:1465
virtual void resizeSpecies()
Resize arrays with sizes that depend on the total number of species.
Definition Kinetics.cpp:607
void checkReactionArraySize(size_t ii) const
Check that an array size is at least nReactions() Throws an exception if ii is less than nReactions()...
Definition Kinetics.cpp:51
size_t nTotalSpecies() const
The total number of species in all phases participating in the kinetics mechanism.
Definition Kinetics.h:252
virtual double reactantStoichCoeff(size_t k, size_t i) const
Stoichiometric coefficient of species k as a reactant in reaction i.
Definition Kinetics.cpp:340
void checkReactionIndex(size_t m) const
Check that the specified reaction index is in range Throws an exception if i is greater than nReactio...
Definition Kinetics.cpp:27
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Definition Kinetics.cpp:411
ThermoPhase & speciesPhase(const string &nm)
This function looks up the name of a species and returns a reference to the ThermoPhase object of the...
Definition Kinetics.cpp:308
vector< double > m_delta_gibbs0
Delta G^0 for all reactions.
Definition Kinetics.h:1521
virtual void getCreationRates(double *cdot)
Species creation rates [kmol/m^3/s or kmol/m^2/s].
Definition Kinetics.cpp:386
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 ...
Definition Kinetics.cpp:329
size_t nDim() const
Returns the number of spatial dimensions (1, 2, or 3)
Definition Phase.h:547
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:142
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition Phase.cpp:129
Abstract base class which stores data about a reaction and its rate parameterization so that it can b...
Definition Reaction.h:25
shared_ptr< ReactionRate > rate()
Get reaction rate pointer.
Definition Reaction.h:147
bool usesThirdBody() const
Check whether reaction involves third body collider.
Definition Reaction.h:161
shared_ptr< ThirdBody > thirdBody()
Get pointer to third-body handler.
Definition Reaction.h:155
bool reversible
True if the current reaction is reversible. False otherwise.
Definition Reaction.h:126
string type() const
The type of reaction, including reaction rate information.
Definition Reaction.cpp:541
string equation() const
The chemical equation for this reaction.
Definition Reaction.cpp:380
Composition products
Product species and stoichiometric coefficients.
Definition Reaction.h:114
Composition reactants
Reactant species and stoichiometric coefficients.
Definition Reaction.h:111
AnyMap input
Input data used for specific models.
Definition Reaction.h:139
bool duplicate
True if the current reaction is marked as duplicate.
Definition Reaction.h:129
void resizeCoeffs(size_t nSpc, size_t nRxn)
Resize the sparse coefficient matrix.
void add(size_t rxn, const vector< size_t > &k)
Add a single reaction to the list of reactions that this stoichiometric manager object handles.
const Eigen::SparseMatrix< double > & stoichCoeffs() const
Return matrix containing stoichiometric coefficients.
Base class for a phase with thermodynamic properties.
A class for managing third-body efficiencies, including default values.
Definition Reaction.h:207
double efficiency(const string &k) const
Get the third-body efficiency for species k
Definition Reaction.cpp:829
Composition efficiencies
Map of species to third body efficiency.
Definition Reaction.h:250
string name() const
Name of the third body collider.
Definition Reaction.h:215
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
virtual Eigen::SparseMatrix< double > fwdRatesOfProgress_ddCi()
Calculate derivatives for forward rates-of-progress with respect to species concentration at constant...
Definition Kinetics.h:810
Eigen::SparseMatrix< double > creationRates_ddCi()
Calculate derivatives for species creation rates with respect to species concentration at constant te...
Definition Kinetics.cpp:468
virtual Eigen::SparseMatrix< double > netRatesOfProgress_ddCi()
Calculate derivatives for net rates-of-progress with respect to species concentration at constant tem...
Definition Kinetics.h:956
void getCreationRates_ddT(double *dwdot)
Calculate derivatives for species creation rates with respect to temperature at constant pressure,...
Definition Kinetics.cpp:422
virtual Eigen::SparseMatrix< double > netRatesOfProgress_ddX()
Calculate derivatives for net rates-of-progress with respect to species mole fractions at constant te...
Definition Kinetics.h:936
Eigen::SparseMatrix< double > destructionRates_ddX()
Calculate derivatives for species destruction rates with respect to species mole fractions at constan...
Definition Kinetics.cpp:514
void getCreationRates_ddC(double *dwdot)
Calculate derivatives for species creation rates with respect to molar concentration at constant temp...
Definition Kinetics.cpp:446
void getDestructionRates_ddP(double *dwdot)
Calculate derivatives for species destruction rates with respect to pressure at constant temperature,...
Definition Kinetics.cpp:490
void getNetProductionRates_ddC(double *dwdot)
Calculate derivatives for species net production rates with respect to molar concentration at constan...
Definition Kinetics.cpp:550
void getDestructionRates_ddT(double *dwdot)
Calculate derivatives for species destruction rates with respect to temperature at constant pressure,...
Definition Kinetics.cpp:478
virtual void getFwdRatesOfProgress_ddP(double *drop)
Calculate derivatives for forward rates-of-progress with respect to pressure at constant temperature,...
Definition Kinetics.h:759
Eigen::SparseMatrix< double > netProductionRates_ddX()
Calculate derivatives for species net production rates with respect to species mole fractions at cons...
Definition Kinetics.cpp:558
Eigen::SparseMatrix< double > creationRates_ddX()
Calculate derivatives for species creation rates with respect to species mole fractions at constant t...
Definition Kinetics.cpp:458
Eigen::SparseMatrix< double > destructionRates_ddCi()
Calculate derivatives for species destruction rates with respect to species concentration at constant...
Definition Kinetics.cpp:524
void getNetProductionRates_ddT(double *dwdot)
Calculate derivatives for species net production rates with respect to temperature at constant pressu...
Definition Kinetics.cpp:534
virtual Eigen::SparseMatrix< double > fwdRatesOfProgress_ddX()
Calculate derivatives for forward rates-of-progress with respect to species mole fractions at constan...
Definition Kinetics.h:790
virtual void getRevRatesOfProgress_ddT(double *drop)
Calculate derivatives for reverse rates-of-progress with respect to temperature at constant pressure,...
Definition Kinetics.h:821
void getCreationRates_ddP(double *dwdot)
Calculate derivatives for species creation rates with respect to pressure at constant temperature,...
Definition Kinetics.cpp:434
virtual void getNetRatesOfProgress_ddT(double *drop)
Calculate derivatives for net rates-of-progress with respect to temperature at constant pressure,...
Definition Kinetics.h:894
virtual void getRevRatesOfProgress_ddP(double *drop)
Calculate derivatives for reverse rates-of-progress with respect to pressure at constant temperature,...
Definition Kinetics.h:832
Eigen::SparseMatrix< double > netProductionRates_ddCi()
Calculate derivatives for species net production rates with respect to species concentration at const...
Definition Kinetics.cpp:563
virtual void getRevRatesOfProgress_ddC(double *drop)
Calculate derivatives for reverse rates-of-progress with respect to molar concentration at constant t...
Definition Kinetics.h:846
void getNetProductionRates_ddP(double *dwdot)
Calculate derivatives for species net production rates with respect to pressure at constant temperatu...
Definition Kinetics.cpp:542
virtual void getNetRatesOfProgress_ddP(double *drop)
Calculate derivatives for net rates-of-progress with respect to pressure at constant temperature,...
Definition Kinetics.h:905
virtual Eigen::SparseMatrix< double > revRatesOfProgress_ddCi()
Calculate derivatives for forward rates-of-progress with respect to species concentration at constant...
Definition Kinetics.h:883
virtual void getFwdRatesOfProgress_ddT(double *drop)
Calculate derivatives for forward rates-of-progress with respect to temperature at constant pressure,...
Definition Kinetics.h:748
void getDestructionRates_ddC(double *dwdot)
Calculate derivatives for species destruction rates with respect to molar concentration at constant t...
Definition Kinetics.cpp:502
virtual void getNetRatesOfProgress_ddC(double *drop)
Calculate derivatives for net rates-of-progress with respect to molar concentration at constant tempe...
Definition Kinetics.h:919
virtual Eigen::SparseMatrix< double > revRatesOfProgress_ddX()
Calculate derivatives for reverse rates-of-progress with respect to species mole fractions at constan...
Definition Kinetics.h:863
virtual void getFwdRatesOfProgress_ddC(double *drop)
Calculate derivatives for forward rates-of-progress with respect to molar concentration at constant t...
Definition Kinetics.h:773
void warn_user(const string &method, const string &msg, const Args &... args)
Print a user warning raised from method as CanteraWarning.
Definition global.h:263
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...