Cantera  3.1.0a1
InterfaceKinetics.cpp
Go to the documentation of this file.
1 /**
2  * @file InterfaceKinetics.cpp
3  */
4 
5 // This file is part of Cantera. See License.txt in the top-level directory or
6 // at https://cantera.org/license.txt for license and copyright information.
7 
12 #include "cantera/base/utilities.h"
13 
14 namespace Cantera
15 {
16 
17 InterfaceKinetics::~InterfaceKinetics()
18 {
19  delete m_integrator;
20 }
21 
23 {
25 
26  // resize buffer
27  m_rbuf0.resize(nReactions());
28  m_rbuf1.resize(nReactions());
29 
30  for (auto& rates : m_interfaceRates) {
31  rates->resize(nTotalSpecies(), nReactions(), nPhases());
32  // @todo ensure that ReactionData are updated; calling rates->update
33  // blocks correct behavior in InterfaceKinetics::_update_rates_T
34  // and running updateROP() is premature
35  }
36 }
37 
39 {
41  m_redo_rates = true;
42 }
43 
45 {
46  // First task is update the electrical potentials from the Phases
48 
49  // Go find the temperature from the surface
50  double T = thermo(0).temperature();
51  m_redo_rates = true;
52  if (T != m_temp || m_redo_rates) {
53  // Calculate the forward rate constant by calling m_rates and store it in m_rfn[]
54  for (size_t n = 0; n < nPhases(); n++) {
56  }
57 
58  // Use the stoichiometric manager to find deltaH for each reaction.
59  getReactionDelta(m_grt.data(), m_dH.data());
60 
61  m_temp = T;
62  m_ROP_ok = false;
63  m_redo_rates = false;
64  }
65 
66  // loop over interface MultiRate evaluators for each reaction type
67  for (auto& rates : m_interfaceRates) {
68  bool changed = rates->update(thermo(0), *this);
69  if (changed) {
70  rates->getRateConstants(m_rfn.data());
71  m_ROP_ok = false;
72  m_redo_rates = true;
73  }
74  }
75 
76  if (!m_ROP_ok) {
77  updateKc();
78  }
79 }
80 
82 {
83  // Store electric potentials for each phase in the array m_phi[].
84  for (size_t n = 0; n < nPhases(); n++) {
85  if (thermo(n).electricPotential() != m_phi[n]) {
86  m_phi[n] = thermo(n).electricPotential();
87  m_redo_rates = true;
88  }
89  }
90 }
91 
93 {
94  for (size_t n = 0; n < nPhases(); n++) {
95  const auto& tp = thermo(n);
96  /*
97  * We call the getActivityConcentrations function of each ThermoPhase
98  * class that makes up this kinetics object to obtain the generalized
99  * concentrations for species within that class. This is collected in
100  * the vector m_conc. m_start[] are integer indices for that vector
101  * denoting the start of the species for each phase.
102  */
103  tp.getActivityConcentrations(m_actConc.data() + m_start[n]);
104 
105  // Get regular concentrations too
106  tp.getConcentrations(m_conc.data() + m_start[n]);
107  }
108  m_ROP_ok = false;
109 }
110 
112 {
113  _update_rates_C();
114  copy(m_actConc.begin(), m_actConc.end(), conc);
115 }
116 
118 {
119  fill(m_rkcn.begin(), m_rkcn.end(), 0.0);
120 
121  if (m_revindex.size() > 0) {
122  /*
123  * Get the vector of standard state electrochemical potentials for
124  * species in the Interfacial kinetics object and store it in m_mu0[]
125  * and m_mu0_Kc[]
126  */
127  updateMu0();
128  double rrt = 1.0 / thermo(0).RT();
129 
130  // compute Delta mu^0 for all reversible reactions
131  getRevReactionDelta(m_mu0_Kc.data(), m_rkcn.data());
132 
133  for (size_t i = 0; i < m_revindex.size(); i++) {
134  size_t irxn = m_revindex[i];
135  if (irxn == npos || irxn >= nReactions()) {
136  throw CanteraError("InterfaceKinetics::updateKc",
137  "illegal value: irxn = {}", irxn);
138  }
139  // WARNING this may overflow HKM
140  m_rkcn[irxn] = exp(m_rkcn[irxn]*rrt);
141  }
142  for (size_t i = 0; i != m_irrev.size(); ++i) {
143  m_rkcn[ m_irrev[i] ] = 0.0;
144  }
145  }
146 }
147 
149 {
150  // First task is update the electrical potentials from the Phases
152 
153  // @todo There is significant potential to further simplify calculations
154  // once the old framework is removed
155  size_t ik = 0;
156  for (size_t n = 0; n < nPhases(); n++) {
158  for (size_t k = 0; k < thermo(n).nSpecies(); k++) {
159  m_mu0_Kc[ik] = m_mu0[ik] + Faraday * m_phi[n] * thermo(n).charge(k);
160  m_mu0_Kc[ik] -= thermo(0).RT() * thermo(n).logStandardConc(k);
161  ik++;
162  }
163  }
164 }
165 
167 {
168  updateMu0();
169  double rrt = 1.0 / thermo(0).RT();
170  std::fill(kc, kc + nReactions(), 0.0);
171  getReactionDelta(m_mu0_Kc.data(), kc);
172  for (size_t i = 0; i < nReactions(); i++) {
173  kc[i] = exp(-kc[i]*rrt);
174  }
175 }
176 
178 {
179  updateROP();
180  for (size_t i = 0; i < nReactions(); i++) {
181  // base rate coefficient multiplied by perturbation factor
182  kfwd[i] = m_rfn[i] * m_perturb[i];
183  }
184 }
185 
186 void InterfaceKinetics::getRevRateConstants(double* krev, bool doIrreversible)
187 {
188  getFwdRateConstants(krev);
189  if (doIrreversible) {
191  for (size_t i = 0; i < nReactions(); i++) {
192  krev[i] /= m_ropnet[i];
193  }
194  } else {
195  for (size_t i = 0; i < nReactions(); i++) {
196  krev[i] *= m_rkcn[i];
197  }
198  }
199 }
200 
202 {
203  // evaluate rate constants and equilibrium constants at temperature and phi
204  // (electric potential)
205  _update_rates_T();
206  // get updated activities (rates updated below)
207  _update_rates_C();
208 
209  if (m_ROP_ok) {
210  return;
211  }
212 
213  for (size_t i = 0; i < nReactions(); i++) {
214  // Scale the base forward rate coefficient by the perturbation factor
215  m_ropf[i] = m_rfn[i] * m_perturb[i];
216  // Multiply the scaled forward rate coefficient by the reciprocal of the
217  // equilibrium constant
218  m_ropr[i] = m_ropf[i] * m_rkcn[i];
219  }
220 
221  // multiply ropf by the activity concentration reaction orders to obtain
222  // the forward rates of progress.
223  m_reactantStoich.multiply(m_actConc.data(), m_ropf.data());
224 
225  // For reversible reactions, multiply ropr by the activity concentration
226  // products
227  m_revProductStoich.multiply(m_actConc.data(), m_ropr.data());
228 
229  for (size_t j = 0; j != nReactions(); ++j) {
230  m_ropnet[j] = m_ropf[j] - m_ropr[j];
231  }
232 
233  // For reactions involving multiple phases, we must check that the phase
234  // being consumed actually exists. This is particularly important for phases
235  // that are stoichiometric phases containing one species with a unity
236  // activity
237  if (m_phaseExistsCheck) {
238  for (size_t j = 0; j != nReactions(); ++j) {
239  if ((m_ropr[j] > m_ropf[j]) && (m_ropr[j] > 0.0)) {
240  for (size_t p = 0; p < nPhases(); p++) {
241  if (m_rxnPhaseIsProduct[j][p] && !m_phaseExists[p]) {
242  m_ropnet[j] = 0.0;
243  m_ropr[j] = m_ropf[j];
244  if (m_ropf[j] > 0.0) {
245  for (size_t rp = 0; rp < nPhases(); rp++) {
246  if (m_rxnPhaseIsReactant[j][rp] && !m_phaseExists[rp]) {
247  m_ropnet[j] = 0.0;
248  m_ropr[j] = m_ropf[j] = 0.0;
249  }
250  }
251  }
252  }
253  if (m_rxnPhaseIsReactant[j][p] && !m_phaseIsStable[p]) {
254  m_ropnet[j] = 0.0;
255  m_ropr[j] = m_ropf[j];
256  }
257  }
258  } else if ((m_ropf[j] > m_ropr[j]) && (m_ropf[j] > 0.0)) {
259  for (size_t p = 0; p < nPhases(); p++) {
260  if (m_rxnPhaseIsReactant[j][p] && !m_phaseExists[p]) {
261  m_ropnet[j] = 0.0;
262  m_ropf[j] = m_ropr[j];
263  if (m_ropf[j] > 0.0) {
264  for (size_t rp = 0; rp < nPhases(); rp++) {
265  if (m_rxnPhaseIsProduct[j][rp] && !m_phaseExists[rp]) {
266  m_ropnet[j] = 0.0;
267  m_ropf[j] = m_ropr[j] = 0.0;
268  }
269  }
270  }
271  }
272  if (m_rxnPhaseIsProduct[j][p] && !m_phaseIsStable[p]) {
273  m_ropnet[j] = 0.0;
274  m_ropf[j] = m_ropr[j];
275  }
276  }
277  }
278  }
279  }
280  m_ROP_ok = true;
281 }
282 
284 {
285  // Get the chemical potentials of the species in the all of the phases used
286  // in the kinetics mechanism
287  for (size_t n = 0; n < nPhases(); n++) {
288  m_thermo[n]->getChemPotentials(m_mu.data() + m_start[n]);
289  }
290 
291  // Use the stoichiometric manager to find deltaG for each reaction.
292  getReactionDelta(m_mu.data(), m_rbuf.data());
293  if (deltaG != 0 && (m_rbuf.data() != deltaG)) {
294  for (size_t j = 0; j < nReactions(); ++j) {
295  deltaG[j] = m_rbuf[j];
296  }
297  }
298 }
299 
301 {
302  // Get the chemical potentials of the species
303  for (size_t n = 0; n < nPhases(); n++) {
305  }
306 
307  // Use the stoichiometric manager to find deltaG for each reaction.
308  getReactionDelta(m_grt.data(), deltaM);
309 }
310 
312 {
313  // Get the partial molar enthalpy of all species
314  for (size_t n = 0; n < nPhases(); n++) {
316  }
317 
318  // Use the stoichiometric manager to find deltaH for each reaction.
319  getReactionDelta(m_grt.data(), deltaH);
320 }
321 
323 {
324  // Get the partial molar entropy of all species in all of the phases
325  for (size_t n = 0; n < nPhases(); n++) {
327  }
328 
329  // Use the stoichiometric manager to find deltaS for each reaction.
330  getReactionDelta(m_grt.data(), deltaS);
331 }
332 
334 {
335  // Get the standard state chemical potentials of the species. This is the
336  // array of chemical potentials at unit activity We define these here as the
337  // chemical potentials of the pure species at the temperature and pressure
338  // of the solution.
339  for (size_t n = 0; n < nPhases(); n++) {
341  }
342 
343  // Use the stoichiometric manager to find deltaG for each reaction.
344  getReactionDelta(m_mu0.data(), deltaGSS);
345 }
346 
348 {
349  // Get the standard state enthalpies of the species. This is the array of
350  // chemical potentials at unit activity We define these here as the
351  // enthalpies of the pure species at the temperature and pressure of the
352  // solution.
353  for (size_t n = 0; n < nPhases(); n++) {
354  thermo(n).getEnthalpy_RT(m_grt.data() + m_start[n]);
355  }
356  for (size_t k = 0; k < m_kk; k++) {
357  m_grt[k] *= thermo(0).RT();
358  }
359 
360  // Use the stoichiometric manager to find deltaH for each reaction.
361  getReactionDelta(m_grt.data(), deltaH);
362 }
363 
365 {
366  // Get the standard state entropy of the species. We define these here as
367  // the entropies of the pure species at the temperature and pressure of the
368  // solution.
369  for (size_t n = 0; n < nPhases(); n++) {
370  thermo(n).getEntropy_R(m_grt.data() + m_start[n]);
371  }
372  for (size_t k = 0; k < m_kk; k++) {
373  m_grt[k] *= GasConstant;
374  }
375 
376  // Use the stoichiometric manager to find deltaS for each reaction.
377  getReactionDelta(m_grt.data(), deltaS);
378 }
379 
380 bool InterfaceKinetics::addReaction(shared_ptr<Reaction> r_base, bool resize)
381 {
382  if (!m_surf) {
383  init();
384  }
385 
386  size_t i = nReactions();
387  bool added = Kinetics::addReaction(r_base, resize);
388  if (!added) {
389  return false;
390  }
391 
392  if (r_base->reversible) {
393  m_revindex.push_back(i);
394  } else {
395  m_irrev.push_back(i);
396  }
397 
398  m_rxnPhaseIsReactant.emplace_back(nPhases(), false);
399  m_rxnPhaseIsProduct.emplace_back(nPhases(), false);
400 
401  for (const auto& [name, stoich] : r_base->reactants) {
402  size_t k = kineticsSpeciesIndex(name);
403  size_t p = speciesPhaseIndex(k);
404  m_rxnPhaseIsReactant[i][p] = true;
405  }
406  for (const auto& [name, stoich] : r_base->products) {
407  size_t k = kineticsSpeciesIndex(name);
408  size_t p = speciesPhaseIndex(k);
409  m_rxnPhaseIsProduct[i][p] = true;
410  }
411 
412  // Set index of rate to number of reaction within kinetics
413  shared_ptr<ReactionRate> rate = r_base->rate();
414  rate->setRateIndex(nReactions() - 1);
415  rate->setContext(*r_base, *this);
416 
417  string rtype = rate->subType();
418  if (rtype == "") {
419  rtype = rate->type();
420  }
421 
422  // If necessary, add new interface MultiRate evaluator
423  if (m_interfaceTypes.find(rtype) == m_interfaceTypes.end()) {
424  m_interfaceTypes[rtype] = m_interfaceRates.size();
425  m_interfaceRates.push_back(rate->newMultiRate());
426  m_interfaceRates.back()->resize(m_kk, nReactions(), nPhases());
427  }
428 
429  // Add reaction rate to evaluator
430  size_t index = m_interfaceTypes[rtype];
431  m_interfaceRates[index]->add(nReactions() - 1, *rate);
432 
433  // Set flag for coverage dependence to true
434  if (rate->compositionDependent()) {
436  }
437 
438  // Set flag for electrochemistry to true
439  if (r_base->usesElectrochemistry(*this)) {
440  m_has_electrochemistry = true;
441  }
442 
443  return true;
444 }
445 
446 void InterfaceKinetics::modifyReaction(size_t i, shared_ptr<Reaction> r_base)
447 {
448  Kinetics::modifyReaction(i, r_base);
449 
450  shared_ptr<ReactionRate> rate = r_base->rate();
451  rate->setRateIndex(i);
452  rate->setContext(*r_base, *this);
453 
454  string rtype = rate->subType();
455  if (rtype == "") {
456  rtype = rate->type();
457  }
458 
459  // Ensure that interface MultiRate evaluator is available
460  if (!m_interfaceTypes.count(rtype)) {
461  throw CanteraError("InterfaceKinetics::modifyReaction",
462  "Interface evaluator not available for type '{}'.", rtype);
463  }
464  // Replace reaction rate evaluator
465  size_t index = m_interfaceTypes[rate->type()];
466  m_interfaceRates[index]->replace(i, *rate);
467 
468  // Invalidate cached data
469  m_redo_rates = true;
470  m_temp += 0.1;
471 }
472 
473 void InterfaceKinetics::setMultiplier(size_t i, double f)
474 {
476  m_ROP_ok = false;
477 }
478 
479 void InterfaceKinetics::setIOFlag(int ioFlag)
480 {
481  m_ioFlag = ioFlag;
482  if (m_integrator) {
483  m_integrator->setIOFlag(ioFlag);
484  }
485 }
486 
487 void InterfaceKinetics::addThermo(shared_ptr<ThermoPhase> thermo)
488 {
490  m_phaseExists.push_back(true);
491  m_phaseIsStable.push_back(true);
492 }
493 
495 {
496  if (thermo(0).nDim() > 2) {
497  throw CanteraError("InterfaceKinetics::init", "no interface phase is present.");
498  }
499 }
500 
502 {
503  size_t kOld = m_kk;
505  if (m_kk != kOld && nReactions()) {
506  throw CanteraError("InterfaceKinetics::resizeSpecies", "Cannot add"
507  " species to InterfaceKinetics after reactions have been added.");
508  }
509  m_actConc.resize(m_kk);
510  m_conc.resize(m_kk);
511  m_mu0.resize(m_kk);
512  m_mu.resize(m_kk);
513  m_mu0_Kc.resize(m_kk);
514  m_grt.resize(m_kk);
515  m_phi.resize(nPhases(), 0.0);
516 }
517 
518 void InterfaceKinetics::advanceCoverages(double tstep, double rtol, double atol,
519  double maxStepSize, size_t maxSteps, size_t maxErrTestFails)
520 {
521  if (m_integrator == 0) {
522  vector<InterfaceKinetics*> k{this};
524  }
525  m_integrator->setTolerances(rtol, atol);
526  m_integrator->setMaxStepSize(maxStepSize);
527  m_integrator->setMaxSteps(maxSteps);
528  m_integrator->setMaxErrTestFails(maxErrTestFails);
529  m_integrator->integrate(0.0, tstep);
530  delete m_integrator;
531  m_integrator = 0;
532 }
533 
535  int ifuncOverride, double timeScaleOverride)
536 {
537  // create our own solver object
538  if (m_integrator == 0) {
539  vector<InterfaceKinetics*> k{this};
542  }
543  m_integrator->setIOFlag(m_ioFlag);
544  // New direct method to go here
545  m_integrator->solvePseudoSteadyStateProblem(ifuncOverride, timeScaleOverride);
546 }
547 
548 void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists)
549 {
550  checkPhaseIndex(iphase);
551  if (exists) {
552  if (!m_phaseExists[iphase]) {
554  m_phaseExistsCheck = std::max(m_phaseExistsCheck, 0);
555  m_phaseExists[iphase] = true;
556  }
557  m_phaseIsStable[iphase] = true;
558  } else {
559  if (m_phaseExists[iphase]) {
561  m_phaseExists[iphase] = false;
562  }
563  m_phaseIsStable[iphase] = false;
564  }
565 }
566 
567 int InterfaceKinetics::phaseExistence(const size_t iphase) const
568 {
569  checkPhaseIndex(iphase);
570  return m_phaseExists[iphase];
571 }
572 
573 int InterfaceKinetics::phaseStability(const size_t iphase) const
574 {
575  checkPhaseIndex(iphase);
576  return m_phaseIsStable[iphase];
577 }
578 
579 void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStable)
580 {
581  checkPhaseIndex(iphase);
582  if (isStable) {
583  m_phaseIsStable[iphase] = true;
584  } else {
585  m_phaseIsStable[iphase] = false;
586  }
587 }
588 
589 double InterfaceKinetics::interfaceCurrent(const size_t iphase)
590 {
591  vector<double> charges(m_kk, 0.0);
592  vector<double> netProdRates(m_kk, 0.0);
593  double dotProduct = 0.0;
594 
595  thermo(iphase).getCharges(charges.data());
596  getNetProductionRates(netProdRates.data());
597 
598  for (size_t k = 0; k < thermo(iphase).nSpecies(); k++)
599  {
600  dotProduct += charges[k] * netProdRates[m_start[iphase] + k];
601  }
602 
603  return dotProduct * Faraday;
604 }
605 
606 Eigen::SparseMatrix<double> InterfaceKinetics::fwdRatesOfProgress_ddCi()
607 {
608  // check derivatives are valid
609  assertDerivativesValid("InterfaceKinetics::fwdRatesOfProgress_ddCi");
610  // forward reaction rate coefficients
611  vector<double>& rop_rates = m_rbuf0;
612  getFwdRateConstants(rop_rates.data());
614 }
615 
616 Eigen::SparseMatrix<double> InterfaceKinetics::revRatesOfProgress_ddCi()
617 {
618  // check derivatives are valid
619  assertDerivativesValid("InterfaceKinetics::revRatesOfProgress_ddCi");
620  // reverse reaction rate coefficients
621  vector<double>& rop_rates = m_rbuf0;
622  getFwdRateConstants(rop_rates.data());
623  applyEquilibriumConstants(rop_rates.data());
625 }
626 
627 Eigen::SparseMatrix<double> InterfaceKinetics::netRatesOfProgress_ddCi()
628 {
629  // check derivatives are valid
630  assertDerivativesValid("InterfaceKinetics::netRatesOfProgress_ddCi");
631  // forward reaction rate coefficients
632  vector<double>& rop_rates = m_rbuf0;
633  getFwdRateConstants(rop_rates.data());
634  Eigen::SparseMatrix<double> jac = calculateCompositionDerivatives(m_reactantStoich,
635  rop_rates);
636 
637  // reverse reaction rate coefficients
638  applyEquilibriumConstants(rop_rates.data());
639  return jac - calculateCompositionDerivatives(m_revProductStoich, rop_rates);
640 }
641 
643 {
644  bool force = settings.empty();
645  if (force || settings.hasKey("skip-coverage-dependence")) {
646  m_jac_skip_coverage_dependence = settings.getBool("skip-coverage-dependence",
647  false);
648  }
649  if (force || settings.hasKey("skip-electrochemistry")) {
650  m_jac_skip_electrochemistry = settings.getBool("skip-electrochemistry",
651  false);
652  }
653  if (force || settings.hasKey("rtol-delta")) {
654  m_jac_rtol_delta = settings.getDouble("rtol-delta", 1e-8);
655  }
656 }
657 
659 {
660  settings["skip-coverage-dependence"] = m_jac_skip_electrochemistry;
661  settings["skip-electrochemistry"] = m_jac_skip_coverage_dependence;
662  settings["rtol-delta"] = m_jac_rtol_delta;
663 }
664 
666  StoichManagerN& stoich, const vector<double>& in)
667 {
668  vector<double>& outV = m_rbuf1;
669  // derivatives handled by StoichManagerN
670  copy(in.begin(), in.end(), outV.begin());
671  return stoich.derivatives(m_actConc.data(), outV.data());
672 }
673 
675 {
677  throw NotImplementedError(name, "Coverage-dependent reactions not supported.");
679  throw NotImplementedError(name, "Electrochemical reactions not supported.");
680  }
681 }
682 
684 {
685  // For reverse rates computed from thermochemistry, multiply the forward
686  // rate coefficients by the reciprocals of the equilibrium constants
687  for (size_t i = 0; i < nReactions(); ++i) {
688  rop[i] *= m_rkcn[i];
689  }
690 }
691 
692 }
Declarations for the implicit integration of surface site density equations (see Kinetics Managers an...
Header for a simple thermodynamics model of a surface phase derived from ThermoPhase,...
A map of string keys to values whose type can vary at runtime.
Definition: AnyMap.h:427
double getDouble(const string &key, double default_) const
If key exists, return it as a double, otherwise return default_.
Definition: AnyMap.cpp:1520
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition: AnyMap.cpp:1423
bool empty() const
Return boolean indicating whether AnyMap is empty.
Definition: AnyMap.cpp:1418
bool getBool(const string &key, bool default_) const
If key exists, return it as a bool, otherwise return default_.
Definition: AnyMap.cpp:1515
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
Advances the surface coverages of the associated set of SurfacePhase objects in time.
void initialize(double t0=0.0)
Must be called before calling method 'advance'.
void integrate(double t0, double t1)
Integrate from t0 to t1. The integrator is reinitialized first.
void setTolerances(double rtol=1.e-7, double atol=1.e-14)
Set the relative and absolute integration tolerances.
void setMaxStepSize(double maxstep=0.0)
Set the maximum integration step-size.
void setMaxSteps(size_t maxsteps=20000)
Set the maximum number of CVODES integration steps.
void solvePseudoSteadyStateProblem(int ifuncOverride=-1, double timeScaleOverride=1.0)
Solve for the pseudo steady-state of the surface problem.
void setMaxErrTestFails(size_t maxErrTestFails=7)
Set the maximum number of CVODES error test failures.
bool addReaction(shared_ptr< Reaction > r, bool resize=true) override
Add a single reaction to the mechanism.
int phaseExistence(const size_t iphase) const
Gets the phase existence int for the ith phase.
double interfaceCurrent(const size_t iphase)
Gets the interface current for the ith phase.
SurfPhase * m_surf
Pointer to the single surface phase.
vector< vector< bool > > m_rxnPhaseIsReactant
Vector of vector of booleans indicating whether a phase participates in a reaction as a reactant.
vector< int > m_phaseIsStable
Vector of int indicating whether phases are stable or not.
Eigen::SparseMatrix< double > revRatesOfProgress_ddCi() override
Calculate derivatives for forward rates-of-progress with respect to species concentration at constant...
void getFwdRateConstants(double *kfwd) override
Return the forward rate constants.
vector< bool > m_phaseExists
Vector of booleans indicating whether phases exist or not.
void getDeltaSSGibbs(double *deltaG) override
Return the vector of values for the reaction standard state Gibbs free energy change.
double m_temp
Current temperature of the data.
void setPhaseStability(const size_t iphase, const int isStable)
Set the stability of a phase in the reaction object.
void getDeltaGibbs(double *deltaG) override
Return the vector of values for the reaction Gibbs free energy change.
bool m_jac_skip_coverage_dependence
A flag used to neglect rate coefficient coverage dependence in derivative formation.
map< string, size_t > m_interfaceTypes
Rate handler mapping.
void _update_rates_phi()
Update properties that depend on the electric potential.
vector< double > m_phi
Vector of phase electric potentials.
void assertDerivativesValid(const string &name)
Helper function ensuring that all rate derivatives can be calculated.
double m_jac_rtol_delta
Relative tolerance used in developing numerical portions of specific derivatives.
void setMultiplier(size_t i, double f) override
Set the multiplier for reaction i to f.
void getDeltaSSEntropy(double *deltaS) override
Return the vector of values for the change in the standard state entropies for each reaction.
void advanceCoverages(double tstep, double rtol=1.e-7, double atol=1.e-14, double maxStepSize=0, size_t maxSteps=20000, size_t maxErrTestFails=7)
Advance the surface coverages in time.
void getDeltaSSEnthalpy(double *deltaH) override
Return the vector of values for the change in the standard state enthalpies of reaction.
void resizeSpecies() override
Resize arrays with sizes that depend on the total number of species.
vector< double > m_grt
Temporary work vector of length m_kk.
void getActivityConcentrations(double *const conc) override
Get the vector of activity concentrations used in the kinetics object.
void getRevRateConstants(double *krev, bool doIrreversible=false) override
Return the reverse rate constants.
void getEquilibriumConstants(double *kc) override
Equilibrium constant for all reactions including the voltage term.
vector< size_t > m_revindex
List of reactions numbers which are reversible reactions.
bool m_has_electrochemistry
A flag stating if the object uses electrochemistry.
void updateKc()
Update the equilibrium constants and stored electrochemical potentials in molar units for all reversi...
void setPhaseExistence(const size_t iphase, const int exists)
Set the existence of a phase in the reaction object.
void updateROP() override
Internal routine that updates the Rates of Progress of the reactions.
void applyEquilibriumConstants(double *rop)
Multiply rate with inverse equilibrium constant.
void getDerivativeSettings(AnyMap &settings) const override
Retrieve derivative settings.
void getDeltaEnthalpy(double *deltaH) override
Return the vector of values for the reactions change in enthalpy.
Eigen::SparseMatrix< double > fwdRatesOfProgress_ddCi() override
Calculate derivatives for forward rates-of-progress with respect to species concentration at constant...
void modifyReaction(size_t i, shared_ptr< Reaction > rNew) override
Modify the rate expression associated with a reaction.
bool m_jac_skip_electrochemistry
A flag used to neglect electrochemical contributions in derivative formation.
void _update_rates_C()
Update properties that depend on the species mole fractions and/or concentration,.
vector< double > m_rbuf0
Buffers for partial rop results with length nReactions()
int phaseStability(const size_t iphase) const
Gets the phase stability int for the ith phase.
vector< unique_ptr< MultiRateBase > > m_interfaceRates
Vector of rate handlers for interface reactions.
vector< double > m_mu0_Kc
Vector of standard state electrochemical potentials modified by a standard concentration term.
ImplicitSurfChem * m_integrator
Pointer to the Implicit surface chemistry object.
void _update_rates_T()
Update properties that depend on temperature.
vector< size_t > m_irrev
Vector of irreversible reaction numbers.
void solvePseudoSteadyStateProblem(int ifuncOverride=-1, double timeScaleOverride=1.0)
Solve for the pseudo steady-state of the surface problem.
bool m_has_coverage_dependence
A flag stating if the object has coverage dependent rates.
virtual void updateMu0()
Update the standard state chemical potentials and species equilibrium constant entries.
void addThermo(shared_ptr< ThermoPhase > thermo) override
Add a thermo phase to the kinetics manager object.
vector< double > m_mu
Vector of chemical potentials for all species.
vector< double > m_actConc
Array of activity concentrations for each species in the kinetics object.
vector< double > m_mu0
Vector of standard state chemical potentials for all species.
void getDeltaElectrochemPotentials(double *deltaM) override
Return the vector of values for the reaction electrochemical free energy change.
void setElectricPotential(int n, double V)
Set the electric potential in the nth phase.
void init() override
Prepare the class for the addition of reactions, after all phases have been added.
void resizeReactions() override
Finalize Kinetics object and associated objects.
int m_phaseExistsCheck
Int flag to indicate that some phases in the kinetics mechanism are non-existent.
Eigen::SparseMatrix< double > calculateCompositionDerivatives(StoichManagerN &stoich, const vector< double > &in)
Process mole fraction derivative.
Eigen::SparseMatrix< double > netRatesOfProgress_ddCi() override
Calculate derivatives for net rates-of-progress with respect to species concentration at constant tem...
void getDeltaEntropy(double *deltaS) override
Return the vector of values for the reactions change in entropy.
vector< vector< bool > > m_rxnPhaseIsProduct
Vector of vector of booleans indicating whether a phase participates in a reaction as a product.
void setDerivativeSettings(const AnyMap &settings) override
Set/modify derivative settings.
vector< double > m_conc
Array of concentrations for each species in the kinetics mechanism.
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
vector< double > m_perturb
Vector of perturbation factors for each reaction's rate of progress vector.
Definition: Kinetics.h:1452
vector< double > m_ropf
Forward rate-of-progress for each reaction.
Definition: Kinetics.h:1496
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:329
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:1473
vector< double > m_rkcn
Reciprocal of the equilibrium constant in concentration units.
Definition: Kinetics.h:1493
size_t m_kk
The number of species in all of the phases that participate in this kinetics mechanism.
Definition: Kinetics.h:1448
virtual void getReactionDelta(const double *property, double *deltaProperty) const
Change in species properties.
Definition: Kinetics.cpp:320
vector< double > m_dH
The enthalpy change for each reaction to calculate Blowers-Masel rates.
Definition: Kinetics.h:1505
vector< double > m_ropr
Reverse rate-of-progress for each reaction.
Definition: Kinetics.h:1499
size_t nPhases() const
The number of phases participating in the reaction mechanism.
Definition: Kinetics.h:184
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:1467
virtual bool addReaction(shared_ptr< Reaction > r, bool resize=true)
Add a single reaction to the mechanism.
Definition: Kinetics.cpp:565
virtual void modifyReaction(size_t i, shared_ptr< Reaction > rNew)
Modify the rate expression associated with a reaction.
Definition: Kinetics.cpp:653
vector< double > m_ropnet
Net rate-of-progress for each reaction.
Definition: Kinetics.h:1502
virtual void addThermo(shared_ptr< ThermoPhase > thermo)
Add a phase to the kinetics manager object.
Definition: Kinetics.cpp:520
vector< double > m_rbuf
Buffer used for storage of intermediate reaction-specific results.
Definition: Kinetics.h:1508
StoichManagerN m_revProductStoich
Stoichiometry manager for the products of reversible reactions.
Definition: Kinetics.h:1437
size_t nReactions() const
Number of reactions in the reaction mechanism.
Definition: Kinetics.h:152
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition: Kinetics.h:276
virtual void setMultiplier(size_t i, double f)
Set the multiplier for reaction i to f.
Definition: Kinetics.h:1369
vector< double > m_rfn
Forward rate constant for each reaction.
Definition: Kinetics.h:1487
StoichManagerN m_reactantStoich
Stoichiometry manager for the reactants for each reaction.
Definition: Kinetics.h:1431
virtual void resizeSpecies()
Resize arrays with sizes that depend on the total number of species.
Definition: Kinetics.cpp:553
size_t nTotalSpecies() const
The total number of species in all phases participating in the kinetics mechanism.
Definition: Kinetics.h:254
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Definition: Kinetics.cpp:363
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
Definition: Kinetics.h:242
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:281
An error indicating that an unimplemented function has been called.
Definition: ctexceptions.h:195
void getCharges(double *charges) const
Copy the vector of species charges into array charges.
Definition: Phase.cpp:405
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:231
double temperature() const
Temperature (K).
Definition: Phase.h:562
double charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Definition: Phase.h:538
This class handles operations involving the stoichiometric coefficients on one side of a reaction (re...
Eigen::SparseMatrix< double > derivatives(const double *conc, const double *rates)
Calculate derivatives with respect to species concentrations.
virtual void getPartialMolarEnthalpies(double *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
Definition: ThermoPhase.h:801
double electricPotential() const
Returns the electric potential of this phase (V).
Definition: ThermoPhase.h:621
virtual void getEntropy_R(double *sr) const
Get the array of nondimensional Entropy functions for the standard state species at the current T and...
Definition: ThermoPhase.h:880
virtual double logStandardConc(size_t k=0) const
Natural logarithm of the standard concentration of the kth species.
Definition: ThermoPhase.cpp:55
double RT() const
Return the Gas Constant multiplied by the current temperature.
Definition: ThermoPhase.h:1062
void getElectrochemPotentials(double *mu) const
Get the species electrochemical potentials.
Definition: ThermoPhase.cpp:76
void setElectricPotential(double v)
Set the electric potential of this phase (V).
Definition: ThermoPhase.h:612
virtual void getStandardChemPotentials(double *mu) const
Get the array of chemical potentials at unit activity for the species at their standard states at the...
Definition: ThermoPhase.h:860
virtual void getEnthalpy_RT(double *hrt) const
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
Definition: ThermoPhase.h:870
virtual void getPartialMolarEntropies(double *sbar) const
Returns an array of partial molar entropies of the species in the solution.
Definition: ThermoPhase.h:811
const double Faraday
Faraday constant [C/kmol].
Definition: ct_defs.h:131
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition: ct_defs.h:120
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
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...