Cantera  3.1.0a1
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 
17 #include "cantera/base/utilities.h"
18 #include "cantera/base/global.h"
19 #include <unordered_set>
20 #include <boost/algorithm/string.hpp>
21 
22 using namespace std;
23 
24 namespace Cantera
25 {
26 
27 void Kinetics::checkReactionIndex(size_t i) const
28 {
29  if (i >= nReactions()) {
30  throw IndexError("Kinetics::checkReactionIndex", "reactions", i,
31  nReactions()-1);
32  }
33 }
34 
35 void Kinetics::resizeReactions()
36 {
37  size_t nRxn = nReactions();
38 
39  // Stoichiometry managers
40  m_reactantStoich.resizeCoeffs(m_kk, nRxn);
41  m_productStoich.resizeCoeffs(m_kk, nRxn);
42  m_revProductStoich.resizeCoeffs(m_kk, nRxn);
43 
44  m_rbuf.resize(nRxn);
45 
46  // products are created for positive net rate of progress
47  m_stoichMatrix = m_productStoich.stoichCoeffs();
48  // reactants are destroyed for positive net rate of progress
49  m_stoichMatrix -= m_reactantStoich.stoichCoeffs();
50 
51  m_ready = true;
52 }
53 
54 void Kinetics::checkReactionArraySize(size_t ii) const
55 {
56  if (nReactions() > ii) {
57  throw ArraySizeError("Kinetics::checkReactionArraySize", ii,
58  nReactions());
59  }
60 }
61 
62 void Kinetics::checkPhaseIndex(size_t m) const
63 {
64  if (m >= nPhases()) {
65  throw IndexError("Kinetics::checkPhaseIndex", "phase", m, nPhases()-1);
66  }
67 }
68 
69 void Kinetics::checkPhaseArraySize(size_t mm) const
70 {
71  if (nPhases() > mm) {
72  throw ArraySizeError("Kinetics::checkPhaseArraySize", mm, nPhases());
73  }
74 }
75 
76 size_t Kinetics::reactionPhaseIndex() const
77 {
78  warn_deprecated("Kinetics::reactionPhaseIndex", "The reacting phase is always "
79  "the first phase. To be removed after Cantera 3.1.");
80  return 0;
81 }
82 
83 shared_ptr<ThermoPhase> Kinetics::reactionPhase() const
84 {
85  return m_thermo[0];
86 }
87 
88 void Kinetics::checkSpeciesIndex(size_t k) const
89 {
90  if (k >= m_kk) {
91  throw IndexError("Kinetics::checkSpeciesIndex", "species", k, m_kk-1);
92  }
93 }
94 
95 void Kinetics::checkSpeciesArraySize(size_t kk) const
96 {
97  if (m_kk > kk) {
98  throw ArraySizeError("Kinetics::checkSpeciesArraySize", kk, m_kk);
99  }
100 }
101 
102 pair<size_t, size_t> Kinetics::checkDuplicates(bool throw_err) const
103 {
104  //! Map of (key indicating participating species) to reaction numbers
105  map<size_t, vector<size_t>> participants;
106  vector<map<int, double>> net_stoich;
107  std::unordered_set<size_t> unmatched_duplicates;
108  for (size_t i = 0; i < m_reactions.size(); i++) {
109  if (m_reactions[i]->duplicate) {
110  unmatched_duplicates.insert(i);
111  }
112  }
113 
114  for (size_t i = 0; i < m_reactions.size(); i++) {
115  // Get data about this reaction
116  unsigned long int key = 0;
117  Reaction& R = *m_reactions[i];
118  net_stoich.emplace_back();
119  map<int, double>& net = net_stoich.back();
120  for (const auto& [name, stoich] : R.reactants) {
121  int k = static_cast<int>(kineticsSpeciesIndex(name));
122  key += k*(k+1);
123  net[-1 -k] -= stoich;
124  }
125  for (const auto& [name, stoich] : R.products) {
126  int k = static_cast<int>(kineticsSpeciesIndex(name));
127  key += k*(k+1);
128  net[1+k] += stoich;
129  }
130 
131  // Compare this reaction to others with similar participants
132  vector<size_t>& related = participants[key];
133  for (size_t m : related) {
134  Reaction& other = *m_reactions[m];
135  if (R.duplicate && other.duplicate) {
136  // marked duplicates
137  unmatched_duplicates.erase(i);
138  unmatched_duplicates.erase(m);
139  continue;
140  } else if (R.type() != other.type()) {
141  continue; // different reaction types
142  } else if (R.rate() && other.rate()
143  && R.rate()->type() != other.rate()->type())
144  {
145  continue; // different rate parameterizations
146  }
147  double c = checkDuplicateStoich(net_stoich[i], net_stoich[m]);
148  if (c == 0) {
149  continue; // stoichiometries differ (not by a multiple)
150  } else if (c < 0.0 && !R.reversible && !other.reversible) {
151  continue; // irreversible reactions in opposite directions
152  } else if (R.usesThirdBody() && other.usesThirdBody()) {
153  ThirdBody& tb1 = *(R.thirdBody());
154  ThirdBody& tb2 = *(other.thirdBody());
155  bool thirdBodyOk = true;
156  for (size_t k = 0; k < nTotalSpecies(); k++) {
157  string s = kineticsSpeciesName(k);
158  if (tb1.efficiency(s) * tb2.efficiency(s) != 0.0) {
159  // non-zero third body efficiencies for species `s` in
160  // both reactions
161  thirdBodyOk = false;
162  break;
163  }
164  }
165  if (thirdBodyOk) {
166  continue; // No overlap in third body efficiencies
167  }
168  }
169  if (throw_err) {
170  throw InputFileError("Kinetics::checkDuplicates",
171  R.input, other.input,
172  "Undeclared duplicate reactions detected:\n"
173  "Reaction {}: {}\nReaction {}: {}\n",
174  i+1, other.equation(), m+1, R.equation());
175  } else {
176  return {i,m};
177  }
178  }
179  participants[key].push_back(i);
180  }
181  if (unmatched_duplicates.size()) {
182  size_t i = *unmatched_duplicates.begin();
183  if (throw_err) {
184  throw InputFileError("Kinetics::checkDuplicates",
185  m_reactions[i]->input,
186  "No duplicate found for declared duplicate reaction number {}"
187  " ({})", i, m_reactions[i]->equation());
188  } else {
189  return {i, i};
190  }
191  }
192  return {npos, npos};
193 }
194 
195 double Kinetics::checkDuplicateStoich(map<int, double>& r1, map<int, double>& r2) const
196 {
197  std::unordered_set<int> keys; // species keys (k+1 or -k-1)
198  for (auto& [speciesKey, stoich] : r1) {
199  keys.insert(speciesKey);
200  }
201  for (auto& [speciesKey, stoich] : r2) {
202  keys.insert(speciesKey);
203  }
204  int k1 = r1.begin()->first;
205  // check for duplicate written in the same direction
206  double ratio = 0.0;
207  if (r1[k1] && r2[k1]) {
208  ratio = r2[k1]/r1[k1];
209  bool different = false;
210  for (int k : keys) {
211  if ((r1[k] && !r2[k]) ||
212  (!r1[k] && r2[k]) ||
213  (r1[k] && fabs(r2[k]/r1[k] - ratio) > 1.e-8)) {
214  different = true;
215  break;
216  }
217  }
218  if (!different) {
219  return ratio;
220  }
221  }
222 
223  // check for duplicate written in the reverse direction
224  if (r1[k1] == 0.0 || r2[-k1] == 0.0) {
225  return 0.0;
226  }
227  ratio = r2[-k1]/r1[k1];
228  for (int k : keys) {
229  if ((r1[k] && !r2[-k]) ||
230  (!r1[k] && r2[-k]) ||
231  (r1[k] && fabs(r2[-k]/r1[k] - ratio) > 1.e-8)) {
232  return 0.0;
233  }
234  }
235  return ratio;
236 }
237 
238 string Kinetics::kineticsSpeciesName(size_t k) const
239 {
240  for (size_t n = m_start.size()-1; n != npos; n--) {
241  if (k >= m_start[n]) {
242  return thermo(n).speciesName(k - m_start[n]);
243  }
244  }
245  return "<unknown>";
246 }
247 
248 size_t Kinetics::kineticsSpeciesIndex(const string& nm) const
249 {
250  for (size_t n = 0; n < m_thermo.size(); n++) {
251  // Check the ThermoPhase object for a match
252  size_t k = thermo(n).speciesIndex(nm);
253  if (k != npos) {
254  return k + m_start[n];
255  }
256  }
257  return npos;
258 }
259 
260 ThermoPhase& Kinetics::speciesPhase(const string& nm)
261 {
262  for (size_t n = 0; n < m_thermo.size(); n++) {
263  size_t k = thermo(n).speciesIndex(nm);
264  if (k != npos) {
265  return thermo(n);
266  }
267  }
268  throw CanteraError("Kinetics::speciesPhase", "unknown species '{}'", nm);
269 }
270 
271 const ThermoPhase& Kinetics::speciesPhase(const string& nm) const
272 {
273  for (const auto& thermo : m_thermo) {
274  if (thermo->speciesIndex(nm) != npos) {
275  return *thermo;
276  }
277  }
278  throw CanteraError("Kinetics::speciesPhase", "unknown species '{}'", nm);
279 }
280 
281 size_t Kinetics::speciesPhaseIndex(size_t k) const
282 {
283  for (size_t n = m_start.size()-1; n != npos; n--) {
284  if (k >= m_start[n]) {
285  return n;
286  }
287  }
288  throw CanteraError("Kinetics::speciesPhaseIndex",
289  "illegal species index: {}", k);
290 }
291 
292 double Kinetics::reactantStoichCoeff(size_t kSpec, size_t irxn) const
293 {
294  return m_reactantStoich.stoichCoeffs().coeff(kSpec, irxn);
295 }
296 
297 double Kinetics::productStoichCoeff(size_t kSpec, size_t irxn) const
298 {
299  return m_productStoich.stoichCoeffs().coeff(kSpec, irxn);
300 }
301 
302 void Kinetics::getFwdRatesOfProgress(double* fwdROP)
303 {
304  updateROP();
305  std::copy(m_ropf.begin(), m_ropf.end(), fwdROP);
306 }
307 
308 void Kinetics::getRevRatesOfProgress(double* revROP)
309 {
310  updateROP();
311  std::copy(m_ropr.begin(), m_ropr.end(), revROP);
312 }
313 
314 void Kinetics::getNetRatesOfProgress(double* netROP)
315 {
316  updateROP();
317  std::copy(m_ropnet.begin(), m_ropnet.end(), netROP);
318 }
319 
320 void Kinetics::getReactionDelta(const double* prop, double* deltaProp) const
321 {
322  fill(deltaProp, deltaProp + nReactions(), 0.0);
323  // products add
324  m_productStoich.incrementReactions(prop, deltaProp);
325  // reactants subtract
326  m_reactantStoich.decrementReactions(prop, deltaProp);
327 }
328 
329 void Kinetics::getRevReactionDelta(const double* prop, double* deltaProp) const
330 {
331  fill(deltaProp, deltaProp + nReactions(), 0.0);
332  // products add
333  m_revProductStoich.incrementReactions(prop, deltaProp);
334  // reactants subtract
335  m_reactantStoich.decrementReactions(prop, deltaProp);
336 }
337 
338 void Kinetics::getCreationRates(double* cdot)
339 {
340  updateROP();
341 
342  // zero out the output array
343  fill(cdot, cdot + m_kk, 0.0);
344 
345  // the forward direction creates product species
346  m_productStoich.incrementSpecies(m_ropf.data(), cdot);
347 
348  // the reverse direction creates reactant species
349  m_reactantStoich.incrementSpecies(m_ropr.data(), cdot);
350 }
351 
352 void Kinetics::getDestructionRates(double* ddot)
353 {
354  updateROP();
355 
356  fill(ddot, ddot + m_kk, 0.0);
357  // the reverse direction destroys products in reversible reactions
358  m_revProductStoich.incrementSpecies(m_ropr.data(), ddot);
359  // the forward direction destroys reactants
360  m_reactantStoich.incrementSpecies(m_ropf.data(), ddot);
361 }
362 
363 void Kinetics::getNetProductionRates(double* net)
364 {
365  updateROP();
366 
367  fill(net, net + m_kk, 0.0);
368  // products are created for positive net rate of progress
369  m_productStoich.incrementSpecies(m_ropnet.data(), net);
370  // reactants are destroyed for positive net rate of progress
371  m_reactantStoich.decrementSpecies(m_ropnet.data(), net);
372 }
373 
374 void Kinetics::getCreationRates_ddT(double* dwdot)
375 {
376  Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
377  Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
378  // the forward direction creates product species
379  getFwdRatesOfProgress_ddT(buf.data());
380  out = m_productStoich.stoichCoeffs() * buf;
381  // the reverse direction creates reactant species
382  getRevRatesOfProgress_ddT(buf.data());
383  out += m_reactantStoich.stoichCoeffs() * buf;
384 }
385 
386 void Kinetics::getCreationRates_ddP(double* dwdot)
387 {
388  Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
389  Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
390  // the forward direction creates product species
391  getFwdRatesOfProgress_ddP(buf.data());
392  out = m_productStoich.stoichCoeffs() * buf;
393  // the reverse direction creates reactant species
394  getRevRatesOfProgress_ddP(buf.data());
395  out += m_reactantStoich.stoichCoeffs() * buf;
396 }
397 
398 void Kinetics::getCreationRates_ddC(double* dwdot)
399 {
400  Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
401  Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
402  // the forward direction creates product species
403  getFwdRatesOfProgress_ddC(buf.data());
404  out = m_productStoich.stoichCoeffs() * buf;
405  // the reverse direction creates reactant species
406  getRevRatesOfProgress_ddC(buf.data());
407  out += m_reactantStoich.stoichCoeffs() * buf;
408 }
409 
410 Eigen::SparseMatrix<double> Kinetics::creationRates_ddX()
411 {
412  Eigen::SparseMatrix<double> jac;
413  // the forward direction creates product species
414  jac = m_productStoich.stoichCoeffs() * fwdRatesOfProgress_ddX();
415  // the reverse direction creates reactant species
416  jac += m_reactantStoich.stoichCoeffs() * revRatesOfProgress_ddX();
417  return jac;
418 }
419 
420 Eigen::SparseMatrix<double> Kinetics::creationRates_ddCi()
421 {
422  Eigen::SparseMatrix<double> jac;
423  // the forward direction creates product species
424  jac = m_productStoich.stoichCoeffs() * fwdRatesOfProgress_ddCi();
425  // the reverse direction creates reactant species
426  jac += m_reactantStoich.stoichCoeffs() * revRatesOfProgress_ddCi();
427  return jac;
428 }
429 
430 void Kinetics::getDestructionRates_ddT(double* dwdot)
431 {
432  Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
433  Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
434  // the reverse direction destroys products in reversible reactions
435  getRevRatesOfProgress_ddT(buf.data());
436  out = m_revProductStoich.stoichCoeffs() * buf;
437  // the forward direction destroys reactants
438  getFwdRatesOfProgress_ddT(buf.data());
439  out += m_reactantStoich.stoichCoeffs() * buf;
440 }
441 
442 void Kinetics::getDestructionRates_ddP(double* dwdot)
443 {
444  Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
445  Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
446  // the reverse direction destroys products in reversible reactions
447  getRevRatesOfProgress_ddP(buf.data());
448  out = m_revProductStoich.stoichCoeffs() * buf;
449  // the forward direction destroys reactants
450  getFwdRatesOfProgress_ddP(buf.data());
451  out += m_reactantStoich.stoichCoeffs() * buf;
452 }
453 
454 void Kinetics::getDestructionRates_ddC(double* dwdot)
455 {
456  Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
457  Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
458  // the reverse direction destroys products in reversible reactions
459  getRevRatesOfProgress_ddC(buf.data());
460  out = m_revProductStoich.stoichCoeffs() * buf;
461  // the forward direction destroys reactants
462  getFwdRatesOfProgress_ddC(buf.data());
463  out += m_reactantStoich.stoichCoeffs() * buf;
464 }
465 
466 Eigen::SparseMatrix<double> Kinetics::destructionRates_ddX()
467 {
468  Eigen::SparseMatrix<double> jac;
469  // the reverse direction destroys products in reversible reactions
470  jac = m_revProductStoich.stoichCoeffs() * revRatesOfProgress_ddX();
471  // the forward direction destroys reactants
472  jac += m_reactantStoich.stoichCoeffs() * fwdRatesOfProgress_ddX();
473  return jac;
474 }
475 
476 Eigen::SparseMatrix<double> Kinetics::destructionRates_ddCi()
477 {
478  Eigen::SparseMatrix<double> jac;
479  // the reverse direction destroys products in reversible reactions
480  jac = m_revProductStoich.stoichCoeffs() * revRatesOfProgress_ddCi();
481  // the forward direction destroys reactants
482  jac += m_reactantStoich.stoichCoeffs() * fwdRatesOfProgress_ddCi();
483  return jac;
484 }
485 
486 void Kinetics::getNetProductionRates_ddT(double* dwdot)
487 {
488  Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
489  Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
490  getNetRatesOfProgress_ddT(buf.data());
491  out = m_stoichMatrix * buf;
492 }
493 
494 void Kinetics::getNetProductionRates_ddP(double* dwdot)
495 {
496  Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
497  Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
498  getNetRatesOfProgress_ddP(buf.data());
499  out = m_stoichMatrix * buf;
500 }
501 
502 void Kinetics::getNetProductionRates_ddC(double* dwdot)
503 {
504  Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
505  Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
506  getNetRatesOfProgress_ddC(buf.data());
507  out = m_stoichMatrix * buf;
508 }
509 
510 Eigen::SparseMatrix<double> Kinetics::netProductionRates_ddX()
511 {
512  return m_stoichMatrix * netRatesOfProgress_ddX();
513 }
514 
515 Eigen::SparseMatrix<double> Kinetics::netProductionRates_ddCi()
516 {
517  return m_stoichMatrix * netRatesOfProgress_ddCi();
518 }
519 
520 void Kinetics::addThermo(shared_ptr<ThermoPhase> thermo)
521 {
522  // the phase with lowest dimensionality is assumed to be the
523  // phase/interface at which reactions take place
524  if (thermo->nDim() <= m_mindim) {
525  if (!m_thermo.empty()) {
526  throw CanteraError("Kinetics::addThermo",
527  "The reacting (lowest dimensional) phase must be added first.");
528  }
529  m_mindim = thermo->nDim();
530  }
531 
532  m_thermo.push_back(thermo);
533  m_phaseindex[m_thermo.back()->name()] = nPhases();
534  resizeSpecies();
535 }
536 
537 AnyMap Kinetics::parameters()
538 {
539  AnyMap out;
540  string name = KineticsFactory::factory()->canonicalize(kineticsType());
541  if (name != "none") {
542  out["kinetics"] = name;
543  if (nReactions() == 0) {
544  out["reactions"] = "none";
545  }
546  if (m_hasUndeclaredThirdBodies) {
547  out["skip-undeclared-third-bodies"] = true;
548  }
549  }
550  return out;
551 }
552 
553 void Kinetics::resizeSpecies()
554 {
555  m_kk = 0;
556  m_start.resize(nPhases());
557 
558  for (size_t i = 0; i < m_thermo.size(); i++) {
559  m_start[i] = m_kk; // global index of first species of phase i
560  m_kk += m_thermo[i]->nSpecies();
561  }
562  invalidateCache();
563 }
564 
565 bool Kinetics::addReaction(shared_ptr<Reaction> r, bool resize)
566 {
567  r->check();
568  r->validate(*this);
569 
570  if (m_kk == 0) {
571  init();
572  }
573  resizeSpecies();
574 
575  // Check validity of reaction within the context of the Kinetics object
576  if (!r->checkSpecies(*this)) {
577  // Do not add reaction
578  return false;
579  }
580 
581  // For reactions created outside the context of a Kinetics object, the units
582  // of the rate coefficient can't be determined in advance. Do that here.
583  if (r->rate_units.factor() == 0) {
584  r->rate()->setRateUnits(r->calculateRateCoeffUnits(*this));
585  }
586 
587  size_t irxn = nReactions(); // index of the new reaction
588 
589  // indices of reactant and product species within this Kinetics object
590  vector<size_t> rk, pk;
591 
592  // Reactant and product stoichiometric coefficients, such that rstoich[i] is
593  // the coefficient for species rk[i]
594  vector<double> rstoich, pstoich;
595 
596  for (const auto& [name, stoich] : r->reactants) {
597  rk.push_back(kineticsSpeciesIndex(name));
598  rstoich.push_back(stoich);
599  }
600 
601  for (const auto& [name, stoich] : r->products) {
602  pk.push_back(kineticsSpeciesIndex(name));
603  pstoich.push_back(stoich);
604  }
605 
606  // The default order for each reactant is its stoichiometric coefficient,
607  // which can be overridden by entries in the Reaction.orders map. rorder[i]
608  // is the order for species rk[i].
609  vector<double> rorder = rstoich;
610  for (const auto& [name, order] : r->orders) {
611  size_t k = kineticsSpeciesIndex(name);
612  // Find the index of species k within rk
613  auto rloc = std::find(rk.begin(), rk.end(), k);
614  if (rloc != rk.end()) {
615  rorder[rloc - rk.begin()] = order;
616  } else {
617  // If the reaction order involves a non-reactant species, add an
618  // extra term to the reactants with zero stoichiometry so that the
619  // stoichiometry manager can be used to compute the global forward
620  // reaction rate.
621  rk.push_back(k);
622  rstoich.push_back(0.0);
623  rorder.push_back(order);
624  }
625  }
626 
627  m_reactantStoich.add(irxn, rk, rorder, rstoich);
628  // product orders = product stoichiometric coefficients
629  m_productStoich.add(irxn, pk, pstoich, pstoich);
630  if (r->reversible) {
631  m_revProductStoich.add(irxn, pk, pstoich, pstoich);
632  }
633 
634  m_reactions.push_back(r);
635  m_rfn.push_back(0.0);
636  m_delta_gibbs0.push_back(0.0);
637  m_rkcn.push_back(0.0);
638  m_ropf.push_back(0.0);
639  m_ropr.push_back(0.0);
640  m_ropnet.push_back(0.0);
641  m_perturb.push_back(1.0);
642  m_dH.push_back(0.0);
643 
644  if (resize) {
645  resizeReactions();
646  } else {
647  m_ready = false;
648  }
649 
650  return true;
651 }
652 
653 void Kinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
654 {
655  checkReactionIndex(i);
656  shared_ptr<Reaction>& rOld = m_reactions[i];
657  if (rNew->type() != rOld->type()) {
658  throw CanteraError("Kinetics::modifyReaction",
659  "Reaction types are different: {} != {}.",
660  rOld->type(), rNew->type());
661  }
662 
663  if (rNew->rate()->type() != rOld->rate()->type()) {
664  throw CanteraError("Kinetics::modifyReaction",
665  "ReactionRate types are different: {} != {}.",
666  rOld->rate()->type(), rNew->rate()->type());
667  }
668 
669  if (rNew->reactants != rOld->reactants) {
670  throw CanteraError("Kinetics::modifyReaction",
671  "Reactants are different: '{}' != '{}'.",
672  rOld->reactantString(), rNew->reactantString());
673  }
674 
675  if (rNew->products != rOld->products) {
676  throw CanteraError("Kinetics::modifyReaction",
677  "Products are different: '{}' != '{}'.",
678  rOld->productString(), rNew->productString());
679  }
680  m_reactions[i] = rNew;
681  invalidateCache();
682 }
683 
684 shared_ptr<Reaction> Kinetics::reaction(size_t i)
685 {
686  checkReactionIndex(i);
687  return m_reactions[i];
688 }
689 
690 shared_ptr<const Reaction> Kinetics::reaction(size_t i) const
691 {
692  checkReactionIndex(i);
693  return m_reactions[i];
694 }
695 
696 }
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.
Definition: ctexceptions.h:135
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
An array index is out of range.
Definition: ctexceptions.h:165
Error thrown for problems processing information contained in an AnyMap or AnyValue.
Definition: AnyMap.h:738
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
bool usesThirdBody() const
Check whether reaction involves third body collider.
Definition: Reaction.h:161
bool reversible
True if the current reaction is reversible. False otherwise.
Definition: Reaction.h:126
shared_ptr< ReactionRate > rate()
Get reaction rate pointer.
Definition: Reaction.h:147
string type() const
The type of reaction, including reaction rate information.
Definition: Reaction.cpp:504
string equation() const
The chemical equation for this reaction.
Definition: Reaction.cpp:345
shared_ptr< ThirdBody > thirdBody()
Get pointer to third-body handler.
Definition: Reaction.h:155
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
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:390
A class for managing third-body efficiencies, including default values.
Definition: Reaction.h:192
double efficiency(const string &k) const
Get the third-body efficiency for species k
Definition: Reaction.cpp:792
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
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:180
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
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...