Cantera  3.1.0a1
Loading...
Searching...
No Matches
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
13
14namespace Cantera
15{
16
17InterfaceKinetics::~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]) {
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{
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
186void InterfaceKinetics::getRevRateConstants(double* krev, bool doIrreversible)
187{
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)
206 // get updated activities (rates updated below)
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
380bool 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)) {
441 }
442
443 return true;
444}
445
446void 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
473void InterfaceKinetics::setMultiplier(size_t i, double f)
474{
476 m_ROP_ok = false;
477}
478
479void InterfaceKinetics::setIOFlag(int ioFlag)
480{
481 m_ioFlag = ioFlag;
482 if (m_integrator) {
483 m_integrator->setIOFlag(ioFlag);
484 }
485}
486
487void 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
518void 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
548void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists)
549{
550 checkPhaseIndex(iphase);
551 if (exists) {
552 if (!m_phaseExists[iphase]) {
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
567int InterfaceKinetics::phaseExistence(const size_t iphase) const
568{
569 checkPhaseIndex(iphase);
570 return m_phaseExists[iphase];
571}
572
573int InterfaceKinetics::phaseStability(const size_t iphase) const
574{
575 checkPhaseIndex(iphase);
576 return m_phaseIsStable[iphase];
577}
578
579void 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
589double 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
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
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
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());
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.
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
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
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
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.
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.
double electricPotential() const
Returns the electric potential of this phase (V).
virtual void getEntropy_R(double *sr) const
Get the array of nondimensional Entropy functions for the standard state species at the current T and...
virtual double logStandardConc(size_t k=0) const
Natural logarithm of the standard concentration of the kth species.
double RT() const
Return the Gas Constant multiplied by the current temperature.
void getElectrochemPotentials(double *mu) const
Get the species electrochemical potentials.
void setElectricPotential(double v)
Set the electric potential of this phase (V).
virtual void getStandardChemPotentials(double *mu) const
Get the array of chemical potentials at unit activity for the species at their standard states at the...
virtual void getEnthalpy_RT(double *hrt) const
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
virtual void getPartialMolarEntropies(double *sbar) const
Returns an array of partial molar entropies of the species in the solution.
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...