Cantera  3.0.0
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,
31 nReactions()-1);
32 }
33}
34
36{
37 size_t nRxn = nReactions();
38
39 // Stoichiometry managers
43
44 m_rbuf.resize(nRxn);
45
46 // products are created for positive net rate of progress
48 // reactants are destroyed for positive net rate of progress
50
51 m_ready = true;
52}
53
55{
56 if (nReactions() > ii) {
57 throw ArraySizeError("Kinetics::checkReactionArraySize", ii,
58 nReactions());
59 }
60}
61
62void Kinetics::checkPhaseIndex(size_t m) const
63{
64 if (m >= nPhases()) {
65 throw IndexError("Kinetics::checkPhaseIndex", "phase", m, nPhases()-1);
66 }
67}
68
69void Kinetics::checkPhaseArraySize(size_t mm) const
70{
71 if (nPhases() > mm) {
72 throw ArraySizeError("Kinetics::checkPhaseArraySize", mm, nPhases());
73 }
74}
75
77{
78 warn_deprecated("Kinetics::surfacePhaseIndex",
79 "To be removed after Cantera 3.0. Use reactionPhaseIndex instead.");
80 return m_surfphase;
81}
82
83shared_ptr<ThermoPhase> Kinetics::reactionPhase() const
84{
85 if (!m_sharedThermo.size()) {
86 // @todo remove after Cantera 3.0
87 throw CanteraError("Kinetics::reactionPhase",
88 "Cannot access phases that were not added using smart pointers.");
89 }
91}
92
93void Kinetics::checkSpeciesIndex(size_t k) const
94{
95 if (k >= m_kk) {
96 throw IndexError("Kinetics::checkSpeciesIndex", "species", k, m_kk-1);
97 }
98}
99
101{
102 if (m_kk > kk) {
103 throw ArraySizeError("Kinetics::checkSpeciesArraySize", kk, m_kk);
104 }
105}
106
107pair<size_t, size_t> Kinetics::checkDuplicates(bool throw_err) const
108{
109 //! Map of (key indicating participating species) to reaction numbers
110 map<size_t, vector<size_t>> participants;
111 vector<map<int, double>> net_stoich;
112 std::unordered_set<size_t> unmatched_duplicates;
113 for (size_t i = 0; i < m_reactions.size(); i++) {
114 if (m_reactions[i]->duplicate) {
115 unmatched_duplicates.insert(i);
116 }
117 }
118
119 for (size_t i = 0; i < m_reactions.size(); i++) {
120 // Get data about this reaction
121 unsigned long int key = 0;
122 Reaction& R = *m_reactions[i];
123 net_stoich.emplace_back();
124 map<int, double>& net = net_stoich.back();
125 for (const auto& [name, stoich] : R.reactants) {
126 int k = static_cast<int>(kineticsSpeciesIndex(name));
127 key += k*(k+1);
128 net[-1 -k] -= stoich;
129 }
130 for (const auto& [name, stoich] : R.products) {
131 int k = static_cast<int>(kineticsSpeciesIndex(name));
132 key += k*(k+1);
133 net[1+k] += stoich;
134 }
135
136 // Compare this reaction to others with similar participants
137 vector<size_t>& related = participants[key];
138 for (size_t m : related) {
139 Reaction& other = *m_reactions[m];
140 if (R.duplicate && other.duplicate) {
141 // marked duplicates
142 unmatched_duplicates.erase(i);
143 unmatched_duplicates.erase(m);
144 continue;
145 } else if (R.type() != other.type()) {
146 continue; // different reaction types
147 } else if (R.rate() && other.rate()
148 && R.rate()->type() != other.rate()->type())
149 {
150 continue; // different rate parameterizations
151 }
152 double c = checkDuplicateStoich(net_stoich[i], net_stoich[m]);
153 if (c == 0) {
154 continue; // stoichiometries differ (not by a multiple)
155 } else if (c < 0.0 && !R.reversible && !other.reversible) {
156 continue; // irreversible reactions in opposite directions
157 } else if (R.usesThirdBody() && other.usesThirdBody()) {
158 ThirdBody& tb1 = *(R.thirdBody());
159 ThirdBody& tb2 = *(other.thirdBody());
160 bool thirdBodyOk = true;
161 for (size_t k = 0; k < nTotalSpecies(); k++) {
162 string s = kineticsSpeciesName(k);
163 if (tb1.efficiency(s) * tb2.efficiency(s) != 0.0) {
164 // non-zero third body efficiencies for species `s` in
165 // both reactions
166 thirdBodyOk = false;
167 break;
168 }
169 }
170 if (thirdBodyOk) {
171 continue; // No overlap in third body efficiencies
172 }
173 }
174 if (throw_err) {
175 throw InputFileError("Kinetics::checkDuplicates",
176 R.input, other.input,
177 "Undeclared duplicate reactions detected:\n"
178 "Reaction {}: {}\nReaction {}: {}\n",
179 i+1, other.equation(), m+1, R.equation());
180 } else {
181 return {i,m};
182 }
183 }
184 participants[key].push_back(i);
185 }
186 if (unmatched_duplicates.size()) {
187 size_t i = *unmatched_duplicates.begin();
188 if (throw_err) {
189 throw InputFileError("Kinetics::checkDuplicates",
190 m_reactions[i]->input,
191 "No duplicate found for declared duplicate reaction number {}"
192 " ({})", i, m_reactions[i]->equation());
193 } else {
194 return {i, i};
195 }
196 }
197 return {npos, npos};
198}
199
200double Kinetics::checkDuplicateStoich(map<int, double>& r1, map<int, double>& r2) const
201{
202 std::unordered_set<int> keys; // species keys (k+1 or -k-1)
203 for (auto& [speciesKey, stoich] : r1) {
204 keys.insert(speciesKey);
205 }
206 for (auto& [speciesKey, stoich] : r2) {
207 keys.insert(speciesKey);
208 }
209 int k1 = r1.begin()->first;
210 // check for duplicate written in the same direction
211 double ratio = 0.0;
212 if (r1[k1] && r2[k1]) {
213 ratio = r2[k1]/r1[k1];
214 bool different = false;
215 for (int k : keys) {
216 if ((r1[k] && !r2[k]) ||
217 (!r1[k] && r2[k]) ||
218 (r1[k] && fabs(r2[k]/r1[k] - ratio) > 1.e-8)) {
219 different = true;
220 break;
221 }
222 }
223 if (!different) {
224 return ratio;
225 }
226 }
227
228 // check for duplicate written in the reverse direction
229 if (r1[k1] == 0.0 || r2[-k1] == 0.0) {
230 return 0.0;
231 }
232 ratio = r2[-k1]/r1[k1];
233 for (int k : keys) {
234 if ((r1[k] && !r2[-k]) ||
235 (!r1[k] && r2[-k]) ||
236 (r1[k] && fabs(r2[-k]/r1[k] - ratio) > 1.e-8)) {
237 return 0.0;
238 }
239 }
240 return ratio;
241}
242
243void Kinetics::selectPhase(const double* data, const ThermoPhase* phase,
244 double* phase_data)
245{
246 warn_deprecated("Kinetics::selectPhase", "To be removed after Cantera 3.0");
247 for (size_t n = 0; n < nPhases(); n++) {
248 if (phase == m_thermo[n]) {
249 size_t nsp = phase->nSpecies();
250 copy(data + m_start[n],
251 data + m_start[n] + nsp, phase_data);
252 return;
253 }
254 }
255 throw CanteraError("Kinetics::selectPhase", "Phase not found.");
256}
257
258string Kinetics::kineticsSpeciesName(size_t k) const
259{
260 for (size_t n = m_start.size()-1; n != npos; n--) {
261 if (k >= m_start[n]) {
262 return thermo(n).speciesName(k - m_start[n]);
263 }
264 }
265 return "<unknown>";
266}
267
268size_t Kinetics::kineticsSpeciesIndex(const string& nm) const
269{
270 for (size_t n = 0; n < m_thermo.size(); n++) {
271 // Check the ThermoPhase object for a match
272 size_t k = thermo(n).speciesIndex(nm);
273 if (k != npos) {
274 return k + m_start[n];
275 }
276 }
277 return npos;
278}
279
280size_t Kinetics::kineticsSpeciesIndex(const string& nm, const string& ph) const
281{
282 warn_deprecated("Kinetics::kineticsSpeciesIndex(species_name, phase_name)",
283 "To be removed after Cantera 3.0. Use kineticsSpeciesIndex(species_name).\n"
284 "Species names should be unique across all phases linked to a Kinetics object.");
285 if (ph == "<any>") {
286 return kineticsSpeciesIndex(nm);
287 }
288
289 for (size_t n = 0; n < m_thermo.size(); n++) {
290 string id = thermo(n).name();
291 if (ph == id) {
292 size_t k = thermo(n).speciesIndex(nm);
293 if (k == npos) {
294 return npos;
295 }
296 return k + m_start[n];
297 }
298 }
299 return npos;
300}
301
303{
304 for (size_t n = 0; n < m_thermo.size(); n++) {
305 size_t k = thermo(n).speciesIndex(nm);
306 if (k != npos) {
307 return thermo(n);
308 }
309 }
310 throw CanteraError("Kinetics::speciesPhase", "unknown species '{}'", nm);
311}
312
313const ThermoPhase& Kinetics::speciesPhase(const string& nm) const
314{
315 for (const auto thermo : m_thermo) {
316 if (thermo->speciesIndex(nm) != npos) {
317 return *thermo;
318 }
319 }
320 throw CanteraError("Kinetics::speciesPhase", "unknown species '{}'", nm);
321}
322
323size_t Kinetics::speciesPhaseIndex(size_t k) const
324{
325 for (size_t n = m_start.size()-1; n != npos; n--) {
326 if (k >= m_start[n]) {
327 return n;
328 }
329 }
330 throw CanteraError("Kinetics::speciesPhaseIndex",
331 "illegal species index: {}", k);
332}
333
334double Kinetics::reactantStoichCoeff(size_t kSpec, size_t irxn) const
335{
336 return m_reactantStoich.stoichCoeffs().coeff(kSpec, irxn);
337}
338
339double Kinetics::productStoichCoeff(size_t kSpec, size_t irxn) const
340{
341 return m_productStoich.stoichCoeffs().coeff(kSpec, irxn);
342}
343
344string Kinetics::reactionType(size_t i) const {
345 warn_deprecated("Kinetics::reactionType", "To be removed after Cantera 3.0.");
346 return m_reactions[i]->type();
347}
348
349string Kinetics::reactionTypeStr(size_t i) const {
350 warn_deprecated("Kinetics::reactionTypeStr", "To be removed after Cantera 3.0.");
351 return reactionType(i);
352}
353
354string Kinetics::reactionString(size_t i) const
355{
356 warn_deprecated("Kinetics::reactionString", "To be removed after Cantera 3.0.");
357 return m_reactions[i]->equation();
358}
359
360//! Returns a string containing the reactants side of the reaction equation.
361string Kinetics::reactantString(size_t i) const
362{
363 warn_deprecated("Kinetics::reactantString", "To be removed after Cantera 3.0.");
364 return m_reactions[i]->reactantString();
365}
366
367//! Returns a string containing the products side of the reaction equation.
368string Kinetics::productString(size_t i) const
369{
370 warn_deprecated("Kinetics::productString", "To be removed after Cantera 3.0.");
371 return m_reactions[i]->productString();
372}
373
375{
376 updateROP();
377 std::copy(m_ropf.begin(), m_ropf.end(), fwdROP);
378}
379
381{
382 updateROP();
383 std::copy(m_ropr.begin(), m_ropr.end(), revROP);
384}
385
387{
388 updateROP();
389 std::copy(m_ropnet.begin(), m_ropnet.end(), netROP);
390}
391
392void Kinetics::getReactionDelta(const double* prop, double* deltaProp) const
393{
394 fill(deltaProp, deltaProp + nReactions(), 0.0);
395 // products add
396 m_productStoich.incrementReactions(prop, deltaProp);
397 // reactants subtract
398 m_reactantStoich.decrementReactions(prop, deltaProp);
399}
400
401void Kinetics::getRevReactionDelta(const double* prop, double* deltaProp) const
402{
403 fill(deltaProp, deltaProp + nReactions(), 0.0);
404 // products add
405 m_revProductStoich.incrementReactions(prop, deltaProp);
406 // reactants subtract
407 m_reactantStoich.decrementReactions(prop, deltaProp);
408}
409
411{
412 updateROP();
413
414 // zero out the output array
415 fill(cdot, cdot + m_kk, 0.0);
416
417 // the forward direction creates product species
418 m_productStoich.incrementSpecies(m_ropf.data(), cdot);
419
420 // the reverse direction creates reactant species
421 m_reactantStoich.incrementSpecies(m_ropr.data(), cdot);
422}
423
425{
426 updateROP();
427
428 fill(ddot, ddot + m_kk, 0.0);
429 // the reverse direction destroys products in reversible reactions
430 m_revProductStoich.incrementSpecies(m_ropr.data(), ddot);
431 // the forward direction destroys reactants
432 m_reactantStoich.incrementSpecies(m_ropf.data(), ddot);
433}
434
436{
437 updateROP();
438
439 fill(net, net + m_kk, 0.0);
440 // products are created for positive net rate of progress
441 m_productStoich.incrementSpecies(m_ropnet.data(), net);
442 // reactants are destroyed for positive net rate of progress
443 m_reactantStoich.decrementSpecies(m_ropnet.data(), net);
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_ddT(buf.data());
452 out = m_productStoich.stoichCoeffs() * buf;
453 // the reverse direction creates reactant species
454 getRevRatesOfProgress_ddT(buf.data());
455 out += m_reactantStoich.stoichCoeffs() * buf;
456}
457
459{
460 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
461 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
462 // the forward direction creates product species
463 getFwdRatesOfProgress_ddP(buf.data());
464 out = m_productStoich.stoichCoeffs() * buf;
465 // the reverse direction creates reactant species
466 getRevRatesOfProgress_ddP(buf.data());
467 out += m_reactantStoich.stoichCoeffs() * buf;
468}
469
471{
472 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
473 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
474 // the forward direction creates product species
475 getFwdRatesOfProgress_ddC(buf.data());
476 out = m_productStoich.stoichCoeffs() * buf;
477 // the reverse direction creates reactant species
478 getRevRatesOfProgress_ddC(buf.data());
479 out += m_reactantStoich.stoichCoeffs() * buf;
480}
481
482Eigen::SparseMatrix<double> Kinetics::creationRates_ddX()
483{
484 Eigen::SparseMatrix<double> jac;
485 // the forward direction creates product species
487 // the reverse direction creates reactant species
489 return jac;
490}
491
492Eigen::SparseMatrix<double> Kinetics::creationRates_ddCi()
493{
494 Eigen::SparseMatrix<double> jac;
495 // the forward direction creates product species
497 // the reverse direction creates reactant species
499 return jac;
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_ddT(buf.data());
508 out = m_revProductStoich.stoichCoeffs() * buf;
509 // the forward direction destroys reactants
510 getFwdRatesOfProgress_ddT(buf.data());
511 out += m_reactantStoich.stoichCoeffs() * buf;
512}
513
515{
516 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
517 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
518 // the reverse direction destroys products in reversible reactions
519 getRevRatesOfProgress_ddP(buf.data());
520 out = m_revProductStoich.stoichCoeffs() * buf;
521 // the forward direction destroys reactants
522 getFwdRatesOfProgress_ddP(buf.data());
523 out += m_reactantStoich.stoichCoeffs() * buf;
524}
525
527{
528 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
529 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
530 // the reverse direction destroys products in reversible reactions
531 getRevRatesOfProgress_ddC(buf.data());
532 out = m_revProductStoich.stoichCoeffs() * buf;
533 // the forward direction destroys reactants
534 getFwdRatesOfProgress_ddC(buf.data());
535 out += m_reactantStoich.stoichCoeffs() * buf;
536}
537
538Eigen::SparseMatrix<double> Kinetics::destructionRates_ddX()
539{
540 Eigen::SparseMatrix<double> jac;
541 // the reverse direction destroys products in reversible reactions
543 // the forward direction destroys reactants
545 return jac;
546}
547
548Eigen::SparseMatrix<double> Kinetics::destructionRates_ddCi()
549{
550 Eigen::SparseMatrix<double> jac;
551 // the reverse direction destroys products in reversible reactions
553 // the forward direction destroys reactants
555 return jac;
556}
557
559{
560 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
561 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
562 getNetRatesOfProgress_ddT(buf.data());
563 out = m_stoichMatrix * buf;
564}
565
567{
568 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
569 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
570 getNetRatesOfProgress_ddP(buf.data());
571 out = m_stoichMatrix * buf;
572}
573
575{
576 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
577 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
578 getNetRatesOfProgress_ddC(buf.data());
579 out = m_stoichMatrix * buf;
580}
581
582Eigen::SparseMatrix<double> Kinetics::netProductionRates_ddX()
583{
585}
586
587Eigen::SparseMatrix<double> Kinetics::netProductionRates_ddCi()
588{
590}
591
592void Kinetics::addThermo(shared_ptr<ThermoPhase> thermo)
593{
594 // the phase with lowest dimensionality is assumed to be the
595 // phase/interface at which reactions take place
596 if (thermo->nDim() <= m_mindim) {
597 if (!m_thermo.empty()) {
598 warn_deprecated("Kinetics::addThermo", "The reacting (lowest dimensional) "
599 "phase should be added first. This warning will become an error after "
600 "Cantera 3.0.");
601 }
602 m_mindim = thermo->nDim();
604 }
605
606 // there should only be one surface phase
607 if (boost::algorithm::contains(thermo->type(), "surface")) {
609 }
610 m_thermo.push_back(thermo.get());
611 m_sharedThermo.push_back(thermo);
612 m_phaseindex[m_thermo.back()->name()] = nPhases();
614}
615
617{
618 warn_deprecated("Kinetics::addPhase",
619 "To be removed after Cantera 3.0. Replaced by addThermo.");
620
621 // the phase with lowest dimensionality is assumed to be the
622 // phase/interface at which reactions take place
623 if (thermo.nDim() <= m_mindim) {
624 m_mindim = thermo.nDim();
626 }
627
628 // there should only be one surface phase
629 if (thermo.type() == kineticsType()) {
631 }
632 m_thermo.push_back(&thermo);
633 m_phaseindex[m_thermo.back()->name()] = nPhases();
635}
636
638{
639 AnyMap out;
640 string name = KineticsFactory::factory()->canonicalize(kineticsType());
641 if (name != "none") {
642 out["kinetics"] = name;
643 if (nReactions() == 0) {
644 out["reactions"] = "none";
645 }
647 out["skip-undeclared-third-bodies"] = true;
648 }
649 }
650 return out;
651}
652
654{
655 m_kk = 0;
656 m_start.resize(nPhases());
657
658 for (size_t i = 0; i < m_thermo.size(); i++) {
659 m_start[i] = m_kk; // global index of first species of phase i
660 m_kk += m_thermo[i]->nSpecies();
661 }
662 invalidateCache();
663}
664
665bool Kinetics::addReaction(shared_ptr<Reaction> r, bool resize)
666{
667 r->check();
668 r->validate(*this);
669
670 if (m_kk == 0) {
671 init();
672 }
674
675 // Check validity of reaction within the context of the Kinetics object
676 if (!r->checkSpecies(*this)) {
677 // Do not add reaction
678 return false;
679 }
680
681 // For reactions created outside the context of a Kinetics object, the units
682 // of the rate coefficient can't be determined in advance. Do that here.
683 if (r->rate_units.factor() == 0) {
684 r->rate()->setRateUnits(r->calculateRateCoeffUnits(*this));
685 }
686
687 size_t irxn = nReactions(); // index of the new reaction
688
689 // indices of reactant and product species within this Kinetics object
690 vector<size_t> rk, pk;
691
692 // Reactant and product stoichiometric coefficients, such that rstoich[i] is
693 // the coefficient for species rk[i]
694 vector<double> rstoich, pstoich;
695
696 for (const auto& [name, stoich] : r->reactants) {
697 rk.push_back(kineticsSpeciesIndex(name));
698 rstoich.push_back(stoich);
699 }
700
701 for (const auto& [name, stoich] : r->products) {
702 pk.push_back(kineticsSpeciesIndex(name));
703 pstoich.push_back(stoich);
704 }
705
706 // The default order for each reactant is its stoichiometric coefficient,
707 // which can be overridden by entries in the Reaction.orders map. rorder[i]
708 // is the order for species rk[i].
709 vector<double> rorder = rstoich;
710 for (const auto& [name, order] : r->orders) {
711 size_t k = kineticsSpeciesIndex(name);
712 // Find the index of species k within rk
713 auto rloc = std::find(rk.begin(), rk.end(), k);
714 if (rloc != rk.end()) {
715 rorder[rloc - rk.begin()] = order;
716 } else {
717 // If the reaction order involves a non-reactant species, add an
718 // extra term to the reactants with zero stoichiometry so that the
719 // stoichiometry manager can be used to compute the global forward
720 // reaction rate.
721 rk.push_back(k);
722 rstoich.push_back(0.0);
723 rorder.push_back(order);
724 }
725 }
726
727 m_reactantStoich.add(irxn, rk, rorder, rstoich);
728 // product orders = product stoichiometric coefficients
729 m_productStoich.add(irxn, pk, pstoich, pstoich);
730 if (r->reversible) {
731 m_revProductStoich.add(irxn, pk, pstoich, pstoich);
732 }
733
734 m_reactions.push_back(r);
735 m_rfn.push_back(0.0);
736 m_delta_gibbs0.push_back(0.0);
737 m_rkcn.push_back(0.0);
738 m_ropf.push_back(0.0);
739 m_ropr.push_back(0.0);
740 m_ropnet.push_back(0.0);
741 m_perturb.push_back(1.0);
742 m_dH.push_back(0.0);
743
744 if (resize) {
746 } else {
747 m_ready = false;
748 }
749
750 return true;
751}
752
753void Kinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
754{
756 shared_ptr<Reaction>& rOld = m_reactions[i];
757 if (rNew->type() != rOld->type()) {
758 throw CanteraError("Kinetics::modifyReaction",
759 "Reaction types are different: {} != {}.",
760 rOld->type(), rNew->type());
761 }
762
763 if (rNew->rate()->type() != rOld->rate()->type()) {
764 throw CanteraError("Kinetics::modifyReaction",
765 "ReactionRate types are different: {} != {}.",
766 rOld->rate()->type(), rNew->rate()->type());
767 }
768
769 if (rNew->reactants != rOld->reactants) {
770 throw CanteraError("Kinetics::modifyReaction",
771 "Reactants are different: '{}' != '{}'.",
772 rOld->reactantString(), rNew->reactantString());
773 }
774
775 if (rNew->products != rOld->products) {
776 throw CanteraError("Kinetics::modifyReaction",
777 "Products are different: '{}' != '{}'.",
778 rOld->productString(), rNew->productString());
779 }
780 m_reactions[i] = rNew;
781 invalidateCache();
782}
783
784shared_ptr<Reaction> Kinetics::reaction(size_t i)
785{
787 return m_reactions[i];
788}
789
790shared_ptr<const Reaction> Kinetics::reaction(size_t i) const
791{
793 return m_reactions[i];
794}
795
796double Kinetics::reactionEnthalpy(const Composition& reactants, const Composition& products)
797{
798 warn_deprecated("Kinetics::reactionEnthalpy", "To be removed after Cantera 3.0");
799 vector<double> hk(nTotalSpecies());
800 for (size_t n = 0; n < nPhases(); n++) {
802 }
803 double rxn_deltaH = 0;
804 for (const auto& [name, stoich] : products) {
805 size_t k = kineticsSpeciesIndex(name);
806 rxn_deltaH += hk[k] * stoich;
807 }
808 for (const auto& [name, stoich] : reactants) {
809 size_t k = kineticsSpeciesIndex(name);
810 rxn_deltaH -= hk[k] * stoich;
811 }
812 return rxn_deltaH;
813}
814
815}
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:427
Array size error.
Base class for exceptions thrown by Cantera classes.
An array index is out of range.
Error thrown for problems processing information contained in an AnyMap or AnyValue.
Definition AnyMap.h:738
virtual void resizeReactions()
Finalize Kinetics object and associated objects.
Definition Kinetics.cpp:35
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:62
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:100
virtual void getFwdRatesOfProgress(double *fwdROP)
Return the forward rates of progress of the reactions.
Definition Kinetics.cpp:374
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
Definition Kinetics.h:254
vector< shared_ptr< Reaction > > m_reactions
Vector of Reaction objects represented by this Kinetics manager.
Definition Kinetics.h:1549
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:93
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:200
vector< double > m_perturb
Vector of perturbation factors for each reaction's rate of progress vector.
Definition Kinetics.h:1546
vector< double > m_ropf
Forward rate-of-progress for each reaction.
Definition Kinetics.h:1608
shared_ptr< Reaction > reaction(size_t i)
Return the Reaction object for reaction i.
Definition Kinetics.cpp:784
bool m_ready
Boolean indicating whether Kinetics object is fully configured.
Definition Kinetics.h:1538
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:401
virtual void getNetRatesOfProgress(double *netROP)
Net rates of progress.
Definition Kinetics.cpp:386
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:1574
virtual string kineticsType() const
Identifies the Kinetics manager type.
Definition Kinetics.h:145
string reactionString(size_t i) const
Return a string representing the reaction.
Definition Kinetics.cpp:354
vector< double > m_rkcn
Reciprocal of the equilibrium constant in concentration units.
Definition Kinetics.h:1605
void selectPhase(const double *data, const ThermoPhase *phase, double *phase_data)
Takes as input an array of properties for all species in the mechanism and copies those values belong...
Definition Kinetics.cpp:243
virtual string reactionTypeStr(size_t i) const
String specifying the type of reaction.
Definition Kinetics.cpp:349
vector< ThermoPhase * > m_thermo
m_thermo is a vector of pointers to ThermoPhase objects that are involved with this kinetics operator
Definition Kinetics.h:1564
shared_ptr< ThermoPhase > reactionPhase() const
Return pointer to phase where the reactions occur.
Definition Kinetics.cpp:83
size_t m_kk
The number of species in all of the phases that participate in this kinetics mechanism.
Definition Kinetics.h:1542
size_t m_rxnphase
Phase Index where reactions are assumed to be taking place.
Definition Kinetics.h:1593
virtual pair< size_t, size_t > checkDuplicates(bool throw_err=true) const
Check for unmarked duplicate reactions and unmatched marked duplicates.
Definition Kinetics.cpp:107
virtual void getReactionDelta(const double *property, double *deltaProperty) const
Change in species properties.
Definition Kinetics.cpp:392
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:69
vector< double > m_dH
The enthalpy change for each reaction to calculate Blowers-Masel rates.
Definition Kinetics.h:1617
vector< double > m_ropr
Reverse rate-of-progress for each reaction.
Definition Kinetics.h:1611
size_t nPhases() const
The number of phases participating in the reaction mechanism.
Definition Kinetics.h:185
virtual void addPhase(ThermoPhase &thermo)
Definition Kinetics.cpp:616
virtual bool addReaction(shared_ptr< Reaction > r, bool resize=true)
Add a single reaction to the mechanism.
Definition Kinetics.cpp:665
virtual string reactionType(size_t i) const
String specifying the type of reaction.
Definition Kinetics.cpp:344
string kineticsSpeciesName(size_t k) const
Return the name of the kth species in the kinetics manager.
Definition Kinetics.cpp:258
virtual void getDestructionRates(double *ddot)
Species destruction rates [kmol/m^3/s or kmol/m^2/s].
Definition Kinetics.cpp:424
AnyMap parameters()
Return the parameters for a phase definition which are needed to reconstruct an identical object usin...
Definition Kinetics.cpp:637
virtual void init()
Prepare the class for the addition of reactions, after all phases have been added.
Definition Kinetics.h:1354
string reactantString(size_t i) const
Returns a string containing the reactants side of the reaction equation.
Definition Kinetics.cpp:361
virtual void modifyReaction(size_t i, shared_ptr< Reaction > rNew)
Modify the rate expression associated with a reaction.
Definition Kinetics.cpp:753
vector< double > m_ropnet
Net rate-of-progress for each reaction.
Definition Kinetics.h:1614
virtual double productStoichCoeff(size_t k, size_t i) const
Stoichiometric coefficient of species k as a product in reaction i.
Definition Kinetics.cpp:339
virtual void addThermo(shared_ptr< ThermoPhase > thermo)
Add a phase to the kinetics manager object.
Definition Kinetics.cpp:592
vector< double > m_rbuf
Buffer used for storage of intermediate reaction-specific results.
Definition Kinetics.h:1620
Eigen::SparseMatrix< double > m_stoichMatrix
Net stoichiometry (products - reactants)
Definition Kinetics.h:1534
map< string, size_t > m_phaseindex
Mapping of the phase name to the position of the phase within the kinetics object.
Definition Kinetics.h:1582
size_t m_mindim
number of spatial dimensions of lowest-dimensional phase.
Definition Kinetics.h:1596
StoichManagerN m_productStoich
Stoichiometry manager for the products for each reaction.
Definition Kinetics.h:1528
StoichManagerN m_revProductStoich
Stoichiometry manager for the products of reversible reactions.
Definition Kinetics.h:1531
virtual double reactionEnthalpy(const Composition &reactants, const Composition &products)
Calculate the reaction enthalpy of a reaction which has not necessarily been added into the Kinetics ...
Definition Kinetics.cpp:796
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:380
bool m_hasUndeclaredThirdBodies
Flag indicating whether reactions include undeclared third bodies.
Definition Kinetics.h:1629
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition Kinetics.h:288
size_t m_surfphase
Index in the list of phases of the one surface phase.
Definition Kinetics.h:1586
vector< double > m_rfn
Forward rate constant for each reaction.
Definition Kinetics.h:1599
StoichManagerN m_reactantStoich
Stoichiometry manager for the reactants for each reaction.
Definition Kinetics.h:1525
size_t surfacePhaseIndex() const
This returns the integer index of the phase which has ThermoPhase type cSurf.
Definition Kinetics.cpp:76
virtual void resizeSpecies()
Resize arrays with sizes that depend on the total number of species.
Definition Kinetics.cpp:653
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:54
vector< shared_ptr< ThermoPhase > > m_sharedThermo
vector of shared pointers,
Definition Kinetics.h:1568
size_t nTotalSpecies() const
The total number of species in all phases participating in the kinetics mechanism.
Definition Kinetics.h:266
virtual double reactantStoichCoeff(size_t k, size_t i) const
Stoichiometric coefficient of species k as a reactant in reaction i.
Definition Kinetics.cpp:334
string productString(size_t i) const
Returns a string containing the products side of the reaction equation.
Definition Kinetics.cpp:368
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:435
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:302
vector< double > m_delta_gibbs0
Delta G^0 for all reactions.
Definition Kinetics.h:1602
virtual void getCreationRates(double *cdot)
Species creation rates [kmol/m^3/s or kmol/m^2/s].
Definition Kinetics.cpp:410
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:323
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:245
size_t nDim() const
Returns the number of spatial dimensions (1, 2, or 3)
Definition Phase.h:646
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:151
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition Phase.cpp:138
string name() const
Return the name of the phase.
Definition Phase.cpp:20
Abstract base class which stores data about a reaction and its rate parameterization so that it can b...
Definition Reaction.h:27
shared_ptr< ReactionRate > rate()
Get reaction rate pointer.
Definition Reaction.h:168
bool usesThirdBody() const
Check whether reaction involves third body collider.
Definition Reaction.h:182
shared_ptr< ThirdBody > thirdBody()
Get pointer to third-body handler.
Definition Reaction.h:176
bool reversible
True if the current reaction is reversible. False otherwise.
Definition Reaction.h:147
string type() const
The type of reaction, including reaction rate information.
Definition Reaction.cpp:506
string equation() const
The chemical equation for this reaction.
Definition Reaction.cpp:347
Composition products
Product species and stoichiometric coefficients.
Definition Reaction.h:135
Composition reactants
Reactant species and stoichiometric coefficients.
Definition Reaction.h:132
AnyMap input
Input data used for specific models.
Definition Reaction.h:160
bool duplicate
True if the current reaction is marked as duplicate.
Definition Reaction.h:150
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.
virtual void getPartialMolarEnthalpies(double *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
string type() const override
String indicating the thermodynamic model implemented.
A class for managing third-body efficiencies, including default values.
Definition Reaction.h:213
double efficiency(const string &k) const
Get the third-body efficiency for species k
Definition Reaction.cpp:809
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:841
Eigen::SparseMatrix< double > creationRates_ddCi()
Calculate derivatives for species creation rates with respect to species concentration at constant te...
Definition Kinetics.cpp:492
virtual Eigen::SparseMatrix< double > netRatesOfProgress_ddCi()
Calculate derivatives for net rates-of-progress with respect to species concentration at constant tem...
Definition Kinetics.h:987
void getCreationRates_ddT(double *dwdot)
Calculate derivatives for species creation rates with respect to temperature at constant pressure,...
Definition Kinetics.cpp:446
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:967
Eigen::SparseMatrix< double > destructionRates_ddX()
Calculate derivatives for species destruction rates with respect to species mole fractions at constan...
Definition Kinetics.cpp:538
void getCreationRates_ddC(double *dwdot)
Calculate derivatives for species creation rates with respect to molar concentration at constant temp...
Definition Kinetics.cpp:470
void getDestructionRates_ddP(double *dwdot)
Calculate derivatives for species destruction rates with respect to pressure at constant temperature,...
Definition Kinetics.cpp:514
void getNetProductionRates_ddC(double *dwdot)
Calculate derivatives for species net production rates with respect to molar concentration at constan...
Definition Kinetics.cpp:574
void getDestructionRates_ddT(double *dwdot)
Calculate derivatives for species destruction rates with respect to temperature at constant pressure,...
Definition Kinetics.cpp:502
virtual void getFwdRatesOfProgress_ddP(double *drop)
Calculate derivatives for forward rates-of-progress with respect to pressure at constant temperature,...
Definition Kinetics.h:790
Eigen::SparseMatrix< double > netProductionRates_ddX()
Calculate derivatives for species net production rates with respect to species mole fractions at cons...
Definition Kinetics.cpp:582
Eigen::SparseMatrix< double > creationRates_ddX()
Calculate derivatives for species creation rates with respect to species mole fractions at constant t...
Definition Kinetics.cpp:482
Eigen::SparseMatrix< double > destructionRates_ddCi()
Calculate derivatives for species destruction rates with respect to species concentration at constant...
Definition Kinetics.cpp:548
void getNetProductionRates_ddT(double *dwdot)
Calculate derivatives for species net production rates with respect to temperature at constant pressu...
Definition Kinetics.cpp:558
virtual Eigen::SparseMatrix< double > fwdRatesOfProgress_ddX()
Calculate derivatives for forward rates-of-progress with respect to species mole fractions at constan...
Definition Kinetics.h:821
virtual void getRevRatesOfProgress_ddT(double *drop)
Calculate derivatives for reverse rates-of-progress with respect to temperature at constant pressure,...
Definition Kinetics.h:852
void getCreationRates_ddP(double *dwdot)
Calculate derivatives for species creation rates with respect to pressure at constant temperature,...
Definition Kinetics.cpp:458
virtual void getNetRatesOfProgress_ddT(double *drop)
Calculate derivatives for net rates-of-progress with respect to temperature at constant pressure,...
Definition Kinetics.h:925
virtual void getRevRatesOfProgress_ddP(double *drop)
Calculate derivatives for reverse rates-of-progress with respect to pressure at constant temperature,...
Definition Kinetics.h:863
Eigen::SparseMatrix< double > netProductionRates_ddCi()
Calculate derivatives for species net production rates with respect to species concentration at const...
Definition Kinetics.cpp:587
virtual void getRevRatesOfProgress_ddC(double *drop)
Calculate derivatives for reverse rates-of-progress with respect to molar concentration at constant t...
Definition Kinetics.h:877
void getNetProductionRates_ddP(double *dwdot)
Calculate derivatives for species net production rates with respect to pressure at constant temperatu...
Definition Kinetics.cpp:566
virtual void getNetRatesOfProgress_ddP(double *drop)
Calculate derivatives for net rates-of-progress with respect to pressure at constant temperature,...
Definition Kinetics.h:936
virtual Eigen::SparseMatrix< double > revRatesOfProgress_ddCi()
Calculate derivatives for forward rates-of-progress with respect to species concentration at constant...
Definition Kinetics.h:914
virtual void getFwdRatesOfProgress_ddT(double *drop)
Calculate derivatives for forward rates-of-progress with respect to temperature at constant pressure,...
Definition Kinetics.h:779
void getDestructionRates_ddC(double *dwdot)
Calculate derivatives for species destruction rates with respect to molar concentration at constant t...
Definition Kinetics.cpp:526
virtual void getNetRatesOfProgress_ddC(double *drop)
Calculate derivatives for net rates-of-progress with respect to molar concentration at constant tempe...
Definition Kinetics.h:950
virtual Eigen::SparseMatrix< double > revRatesOfProgress_ddX()
Calculate derivatives for reverse rates-of-progress with respect to species mole fractions at constan...
Definition Kinetics.h:894
virtual void getFwdRatesOfProgress_ddC(double *drop)
Calculate derivatives for forward rates-of-progress with respect to molar concentration at constant t...
Definition Kinetics.h:804
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:195
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
Definition AnyMap.cpp:1926
map< string, double > Composition
Map from string names to doubles.
Definition ct_defs.h:184
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...