Cantera  3.0.0
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
19{
20 warn_deprecated("InterfaceKinetics::InterfaceKinetics(ThermoPhase*)",
21 "To be removed after Cantera 3.0. Use default constructor instead.");
23}
24
25InterfaceKinetics::~InterfaceKinetics()
26{
27 delete m_integrator;
28}
29
31{
33
34 // resize buffer
35 m_rbuf0.resize(nReactions());
36 m_rbuf1.resize(nReactions());
37
38 for (auto& rates : m_interfaceRates) {
39 rates->resize(nTotalSpecies(), nReactions(), nPhases());
40 // @todo ensure that ReactionData are updated; calling rates->update
41 // blocks correct behavior in InterfaceKinetics::_update_rates_T
42 // and running updateROP() is premature
43 }
44}
45
47{
49 m_redo_rates = true;
50}
51
53{
54 // First task is update the electrical potentials from the Phases
56
57 // Go find the temperature from the surface
59 m_redo_rates = true;
60 if (T != m_temp || m_redo_rates) {
61 // Calculate the forward rate constant by calling m_rates and store it in m_rfn[]
62 for (size_t n = 0; n < nPhases(); n++) {
64 }
65
66 // Use the stoichiometric manager to find deltaH for each reaction.
67 getReactionDelta(m_grt.data(), m_dH.data());
68
69 m_temp = T;
70 m_ROP_ok = false;
71 m_redo_rates = false;
72 }
73
74 // loop over interface MultiRate evaluators for each reaction type
75 for (auto& rates : m_interfaceRates) {
76 bool changed = rates->update(thermo(reactionPhaseIndex()), *this);
77 if (changed) {
78 rates->getRateConstants(m_rfn.data());
79 m_ROP_ok = false;
80 m_redo_rates = true;
81 }
82 }
83
84 if (!m_ROP_ok) {
85 updateKc();
86 }
87}
88
90{
91 // Store electric potentials for each phase in the array m_phi[].
92 for (size_t n = 0; n < nPhases(); n++) {
93 if (thermo(n).electricPotential() != m_phi[n]) {
95 m_redo_rates = true;
96 }
97 }
98}
99
101{
102 for (size_t n = 0; n < nPhases(); n++) {
103 const ThermoPhase* tp = m_thermo[n];
104 /*
105 * We call the getActivityConcentrations function of each ThermoPhase
106 * class that makes up this kinetics object to obtain the generalized
107 * concentrations for species within that class. This is collected in
108 * the vector m_conc. m_start[] are integer indices for that vector
109 * denoting the start of the species for each phase.
110 */
112
113 // Get regular concentrations too
114 tp->getConcentrations(m_conc.data() + m_start[n]);
115 }
116 m_ROP_ok = false;
117}
118
120{
122 copy(m_actConc.begin(), m_actConc.end(), conc);
123}
124
126{
127 fill(m_rkcn.begin(), m_rkcn.end(), 0.0);
128
129 if (m_revindex.size() > 0) {
130 /*
131 * Get the vector of standard state electrochemical potentials for
132 * species in the Interfacial kinetics object and store it in m_mu0[]
133 * and m_mu0_Kc[]
134 */
135 updateMu0();
136 double rrt = 1.0 / thermo(reactionPhaseIndex()).RT();
137
138 // compute Delta mu^0 for all reversible reactions
139 getRevReactionDelta(m_mu0_Kc.data(), m_rkcn.data());
140
141 for (size_t i = 0; i < m_revindex.size(); i++) {
142 size_t irxn = m_revindex[i];
143 if (irxn == npos || irxn >= nReactions()) {
144 throw CanteraError("InterfaceKinetics::updateKc",
145 "illegal value: irxn = {}", irxn);
146 }
147 // WARNING this may overflow HKM
148 m_rkcn[irxn] = exp(m_rkcn[irxn]*rrt);
149 }
150 for (size_t i = 0; i != m_irrev.size(); ++i) {
151 m_rkcn[ m_irrev[i] ] = 0.0;
152 }
153 }
154}
155
157{
158 // First task is update the electrical potentials from the Phases
160
161 // @todo There is significant potential to further simplify calculations
162 // once the old framework is removed
163 size_t ik = 0;
164 for (size_t n = 0; n < nPhases(); n++) {
166 for (size_t k = 0; k < thermo(n).nSpecies(); k++) {
167 m_mu0_Kc[ik] = m_mu0[ik] + Faraday * m_phi[n] * thermo(n).charge(k);
169 * thermo(n).logStandardConc(k);
170 ik++;
171 }
172 }
173}
174
176{
177 updateMu0();
178 double rrt = 1.0 / thermo(reactionPhaseIndex()).RT();
179 std::fill(kc, kc + nReactions(), 0.0);
180 getReactionDelta(m_mu0_Kc.data(), kc);
181 for (size_t i = 0; i < nReactions(); i++) {
182 kc[i] = exp(-kc[i]*rrt);
183 }
184}
185
187{
188 updateROP();
189 for (size_t i = 0; i < nReactions(); i++) {
190 // base rate coefficient multiplied by perturbation factor
191 kfwd[i] = m_rfn[i] * m_perturb[i];
192 }
193}
194
195void InterfaceKinetics::getRevRateConstants(double* krev, bool doIrreversible)
196{
198 if (doIrreversible) {
200 for (size_t i = 0; i < nReactions(); i++) {
201 krev[i] /= m_ropnet[i];
202 }
203 } else {
204 for (size_t i = 0; i < nReactions(); i++) {
205 krev[i] *= m_rkcn[i];
206 }
207 }
208}
209
211{
212 // evaluate rate constants and equilibrium constants at temperature and phi
213 // (electric potential)
215 // get updated activities (rates updated below)
217
218 if (m_ROP_ok) {
219 return;
220 }
221
222 for (size_t i = 0; i < nReactions(); i++) {
223 // Scale the base forward rate coefficient by the perturbation factor
224 m_ropf[i] = m_rfn[i] * m_perturb[i];
225 // Multiply the scaled forward rate coefficient by the reciprocal of the
226 // equilibrium constant
227 m_ropr[i] = m_ropf[i] * m_rkcn[i];
228 }
229
230 // multiply ropf by the activity concentration reaction orders to obtain
231 // the forward rates of progress.
232 m_reactantStoich.multiply(m_actConc.data(), m_ropf.data());
233
234 // For reversible reactions, multiply ropr by the activity concentration
235 // products
236 m_revProductStoich.multiply(m_actConc.data(), m_ropr.data());
237
238 for (size_t j = 0; j != nReactions(); ++j) {
239 m_ropnet[j] = m_ropf[j] - m_ropr[j];
240 }
241
242 // For reactions involving multiple phases, we must check that the phase
243 // being consumed actually exists. This is particularly important for phases
244 // that are stoichiometric phases containing one species with a unity
245 // activity
246 if (m_phaseExistsCheck) {
247 for (size_t j = 0; j != nReactions(); ++j) {
248 if ((m_ropr[j] > m_ropf[j]) && (m_ropr[j] > 0.0)) {
249 for (size_t p = 0; p < nPhases(); p++) {
250 if (m_rxnPhaseIsProduct[j][p] && !m_phaseExists[p]) {
251 m_ropnet[j] = 0.0;
252 m_ropr[j] = m_ropf[j];
253 if (m_ropf[j] > 0.0) {
254 for (size_t rp = 0; rp < nPhases(); rp++) {
255 if (m_rxnPhaseIsReactant[j][rp] && !m_phaseExists[rp]) {
256 m_ropnet[j] = 0.0;
257 m_ropr[j] = m_ropf[j] = 0.0;
258 }
259 }
260 }
261 }
262 if (m_rxnPhaseIsReactant[j][p] && !m_phaseIsStable[p]) {
263 m_ropnet[j] = 0.0;
264 m_ropr[j] = m_ropf[j];
265 }
266 }
267 } else if ((m_ropf[j] > m_ropr[j]) && (m_ropf[j] > 0.0)) {
268 for (size_t p = 0; p < nPhases(); p++) {
269 if (m_rxnPhaseIsReactant[j][p] && !m_phaseExists[p]) {
270 m_ropnet[j] = 0.0;
271 m_ropf[j] = m_ropr[j];
272 if (m_ropf[j] > 0.0) {
273 for (size_t rp = 0; rp < nPhases(); rp++) {
274 if (m_rxnPhaseIsProduct[j][rp] && !m_phaseExists[rp]) {
275 m_ropnet[j] = 0.0;
276 m_ropf[j] = m_ropr[j] = 0.0;
277 }
278 }
279 }
280 }
281 if (m_rxnPhaseIsProduct[j][p] && !m_phaseIsStable[p]) {
282 m_ropnet[j] = 0.0;
283 m_ropf[j] = m_ropr[j];
284 }
285 }
286 }
287 }
288 }
289 m_ROP_ok = true;
290}
291
293{
294 // Get the chemical potentials of the species in the all of the phases used
295 // in the kinetics mechanism
296 for (size_t n = 0; n < nPhases(); n++) {
297 m_thermo[n]->getChemPotentials(m_mu.data() + m_start[n]);
298 }
299
300 // Use the stoichiometric manager to find deltaG for each reaction.
301 getReactionDelta(m_mu.data(), m_rbuf.data());
302 if (deltaG != 0 && (m_rbuf.data() != deltaG)) {
303 for (size_t j = 0; j < nReactions(); ++j) {
304 deltaG[j] = m_rbuf[j];
305 }
306 }
307}
308
310{
311 // Get the chemical potentials of the species
312 for (size_t n = 0; n < nPhases(); n++) {
314 }
315
316 // Use the stoichiometric manager to find deltaG for each reaction.
317 getReactionDelta(m_grt.data(), deltaM);
318}
319
321{
322 // Get the partial molar enthalpy of all species
323 for (size_t n = 0; n < nPhases(); n++) {
325 }
326
327 // Use the stoichiometric manager to find deltaH for each reaction.
328 getReactionDelta(m_grt.data(), deltaH);
329}
330
332{
333 // Get the partial molar entropy of all species in all of the phases
334 for (size_t n = 0; n < nPhases(); n++) {
336 }
337
338 // Use the stoichiometric manager to find deltaS for each reaction.
339 getReactionDelta(m_grt.data(), deltaS);
340}
341
343{
344 // Get the standard state chemical potentials of the species. This is the
345 // array of chemical potentials at unit activity We define these here as the
346 // chemical potentials of the pure species at the temperature and pressure
347 // of the solution.
348 for (size_t n = 0; n < nPhases(); n++) {
350 }
351
352 // Use the stoichiometric manager to find deltaG for each reaction.
353 getReactionDelta(m_mu0.data(), deltaGSS);
354}
355
357{
358 // Get the standard state enthalpies of the species. This is the array of
359 // chemical potentials at unit activity We define these here as the
360 // enthalpies of the pure species at the temperature and pressure of the
361 // solution.
362 for (size_t n = 0; n < nPhases(); n++) {
363 thermo(n).getEnthalpy_RT(m_grt.data() + m_start[n]);
364 }
365 for (size_t k = 0; k < m_kk; k++) {
367 }
368
369 // Use the stoichiometric manager to find deltaH for each reaction.
370 getReactionDelta(m_grt.data(), deltaH);
371}
372
374{
375 // Get the standard state entropy of the species. We define these here as
376 // the entropies of the pure species at the temperature and pressure of the
377 // solution.
378 for (size_t n = 0; n < nPhases(); n++) {
379 thermo(n).getEntropy_R(m_grt.data() + m_start[n]);
380 }
381 for (size_t k = 0; k < m_kk; k++) {
382 m_grt[k] *= GasConstant;
383 }
384
385 // Use the stoichiometric manager to find deltaS for each reaction.
386 getReactionDelta(m_grt.data(), deltaS);
387}
388
389bool InterfaceKinetics::addReaction(shared_ptr<Reaction> r_base, bool resize)
390{
391 if (!m_surf) {
392 init();
393 }
394
395 size_t i = nReactions();
396 bool added = Kinetics::addReaction(r_base, resize);
397 if (!added) {
398 return false;
399 }
400
401 if (r_base->reversible) {
402 m_revindex.push_back(i);
403 } else {
404 m_irrev.push_back(i);
405 }
406
407 m_rxnPhaseIsReactant.emplace_back(nPhases(), false);
408 m_rxnPhaseIsProduct.emplace_back(nPhases(), false);
409
410 for (const auto& [name, stoich] : r_base->reactants) {
411 size_t k = kineticsSpeciesIndex(name);
412 size_t p = speciesPhaseIndex(k);
413 m_rxnPhaseIsReactant[i][p] = true;
414 }
415 for (const auto& [name, stoich] : r_base->products) {
416 size_t k = kineticsSpeciesIndex(name);
417 size_t p = speciesPhaseIndex(k);
418 m_rxnPhaseIsProduct[i][p] = true;
419 }
420
421 // Set index of rate to number of reaction within kinetics
422 shared_ptr<ReactionRate> rate = r_base->rate();
423 rate->setRateIndex(nReactions() - 1);
424 rate->setContext(*r_base, *this);
425
426 string rtype = rate->subType();
427 if (rtype == "") {
428 rtype = rate->type();
429 }
430
431 // If necessary, add new interface MultiRate evaluator
432 if (m_interfaceTypes.find(rtype) == m_interfaceTypes.end()) {
433 m_interfaceTypes[rtype] = m_interfaceRates.size();
434 m_interfaceRates.push_back(rate->newMultiRate());
435 m_interfaceRates.back()->resize(m_kk, nReactions(), nPhases());
436 }
437
438 // Add reaction rate to evaluator
439 size_t index = m_interfaceTypes[rtype];
440 m_interfaceRates[index]->add(nReactions() - 1, *rate);
441
442 // Set flag for coverage dependence to true
443 if (rate->compositionDependent()) {
445 }
446
447 // Set flag for electrochemistry to true
448 if (r_base->usesElectrochemistry(*this)) {
450 }
451
452 return true;
453}
454
455void InterfaceKinetics::modifyReaction(size_t i, shared_ptr<Reaction> r_base)
456{
457 Kinetics::modifyReaction(i, r_base);
458
459 shared_ptr<ReactionRate> rate = r_base->rate();
460 rate->setRateIndex(i);
461 rate->setContext(*r_base, *this);
462
463 string rtype = rate->subType();
464 if (rtype == "") {
465 rtype = rate->type();
466 }
467
468 // Ensure that interface MultiRate evaluator is available
469 if (!m_interfaceTypes.count(rtype)) {
470 throw CanteraError("InterfaceKinetics::modifyReaction",
471 "Interface evaluator not available for type '{}'.", rtype);
472 }
473 // Replace reaction rate evaluator
474 size_t index = m_interfaceTypes[rate->type()];
475 m_interfaceRates[index]->replace(i, *rate);
476
477 // Invalidate cached data
478 m_redo_rates = true;
479 m_temp += 0.1;
480}
481
482void InterfaceKinetics::setMultiplier(size_t i, double f)
483{
485 m_ROP_ok = false;
486}
487
488void InterfaceKinetics::setIOFlag(int ioFlag)
489{
490 m_ioFlag = ioFlag;
491 if (m_integrator) {
492 m_integrator->setIOFlag(ioFlag);
493 }
494}
495
496void InterfaceKinetics::addThermo(shared_ptr<ThermoPhase> thermo)
497{
499 m_phaseExists.push_back(true);
500 m_phaseIsStable.push_back(true);
501}
502
504{
506 m_phaseExists.push_back(true);
507 m_phaseIsStable.push_back(true);
508}
509
511{
512 size_t ks = reactionPhaseIndex();
513 if (ks == npos) {
514 throw CanteraError("InterfaceKinetics::init",
515 "no surface phase is present.");
516 }
517
518 // Check to see that the interface routine has a dimension of 2
519 m_surf = (SurfPhase*)&thermo(ks);
520 if (m_surf->nDim() != m_nDim) {
521 throw CanteraError("InterfaceKinetics::init",
522 "expected interface dimension = 2, but got dimension = {}",
523 m_surf->nDim());
524 }
525}
526
528{
529 size_t kOld = m_kk;
531 if (m_kk != kOld && nReactions()) {
532 throw CanteraError("InterfaceKinetics::resizeSpecies", "Cannot add"
533 " species to InterfaceKinetics after reactions have been added.");
534 }
535 m_actConc.resize(m_kk);
536 m_conc.resize(m_kk);
537 m_mu0.resize(m_kk);
538 m_mu.resize(m_kk);
539 m_mu0_Kc.resize(m_kk);
540 m_grt.resize(m_kk);
541 m_phi.resize(nPhases(), 0.0);
542}
543
544void InterfaceKinetics::advanceCoverages(double tstep, double rtol, double atol,
545 double maxStepSize, size_t maxSteps, size_t maxErrTestFails)
546{
547 if (m_integrator == 0) {
548 vector<InterfaceKinetics*> k{this};
550 }
551 m_integrator->setTolerances(rtol, atol);
552 m_integrator->setMaxStepSize(maxStepSize);
553 m_integrator->setMaxSteps(maxSteps);
554 m_integrator->setMaxErrTestFails(maxErrTestFails);
555 m_integrator->integrate(0.0, tstep);
556 delete m_integrator;
557 m_integrator = 0;
558}
559
561 int ifuncOverride, double timeScaleOverride)
562{
563 // create our own solver object
564 if (m_integrator == 0) {
565 vector<InterfaceKinetics*> k{this};
568 }
569 m_integrator->setIOFlag(m_ioFlag);
570 // New direct method to go here
571 m_integrator->solvePseudoSteadyStateProblem(ifuncOverride, timeScaleOverride);
572}
573
574void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists)
575{
576 checkPhaseIndex(iphase);
577 if (exists) {
578 if (!m_phaseExists[iphase]) {
581 m_phaseExists[iphase] = true;
582 }
583 m_phaseIsStable[iphase] = true;
584 } else {
585 if (m_phaseExists[iphase]) {
587 m_phaseExists[iphase] = false;
588 }
589 m_phaseIsStable[iphase] = false;
590 }
591}
592
593int InterfaceKinetics::phaseExistence(const size_t iphase) const
594{
595 checkPhaseIndex(iphase);
596 return m_phaseExists[iphase];
597}
598
599int InterfaceKinetics::phaseStability(const size_t iphase) const
600{
601 checkPhaseIndex(iphase);
602 return m_phaseIsStable[iphase];
603}
604
605void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStable)
606{
607 checkPhaseIndex(iphase);
608 if (isStable) {
609 m_phaseIsStable[iphase] = true;
610 } else {
611 m_phaseIsStable[iphase] = false;
612 }
613}
614
615double InterfaceKinetics::interfaceCurrent(const size_t iphase)
616{
617 vector<double> charges(m_kk, 0.0);
618 vector<double> netProdRates(m_kk, 0.0);
619 double dotProduct = 0.0;
620
621 thermo(iphase).getCharges(charges.data());
622 getNetProductionRates(netProdRates.data());
623
624 for (size_t k = 0; k < thermo(iphase).nSpecies(); k++)
625 {
626 dotProduct += charges[k] * netProdRates[m_start[iphase] + k];
627 }
628
629 return dotProduct * Faraday;
630}
631
633{
634 // check derivatives are valid
635 assertDerivativesValid("InterfaceKinetics::fwdRatesOfProgress_ddCi");
636 // forward reaction rate coefficients
637 vector<double>& rop_rates = m_rbuf0;
638 getFwdRateConstants(rop_rates.data());
640}
641
643{
644 // check derivatives are valid
645 assertDerivativesValid("InterfaceKinetics::revRatesOfProgress_ddCi");
646 // reverse reaction rate coefficients
647 vector<double>& rop_rates = m_rbuf0;
648 getFwdRateConstants(rop_rates.data());
649 applyEquilibriumConstants(rop_rates.data());
651}
652
654{
655 // check derivatives are valid
656 assertDerivativesValid("InterfaceKinetics::netRatesOfProgress_ddCi");
657 // forward reaction rate coefficients
658 vector<double>& rop_rates = m_rbuf0;
659 getFwdRateConstants(rop_rates.data());
660 Eigen::SparseMatrix<double> jac = calculateCompositionDerivatives(m_reactantStoich,
661 rop_rates);
662
663 // reverse reaction rate coefficients
664 applyEquilibriumConstants(rop_rates.data());
666}
667
669{
670 bool force = settings.empty();
671 if (force || settings.hasKey("skip-coverage-dependence")) {
672 m_jac_skip_coverage_dependence = settings.getBool("skip-coverage-dependence",
673 false);
674 }
675 if (force || settings.hasKey("skip-electrochemistry")) {
676 m_jac_skip_electrochemistry = settings.getBool("skip-electrochemistry",
677 false);
678 }
679 if (force || settings.hasKey("rtol-delta")) {
680 m_jac_rtol_delta = settings.getDouble("rtol-delta", 1e-8);
681 }
682}
683
685{
686 settings["skip-coverage-dependence"] = m_jac_skip_electrochemistry;
687 settings["skip-electrochemistry"] = m_jac_skip_coverage_dependence;
688 settings["rtol-delta"] = m_jac_rtol_delta;
689}
690
692 StoichManagerN& stoich, const vector<double>& in)
693{
694 vector<double>& outV = m_rbuf1;
695 // derivatives handled by StoichManagerN
696 copy(in.begin(), in.end(), outV.begin());
697 return stoich.derivatives(m_actConc.data(), outV.data());
698}
699
701{
703 throw NotImplementedError(name, "Coverage-dependent reactions not supported.");
705 throw NotImplementedError(name, "Electrochemical reactions not supported.");
706 }
707}
708
710{
711 // For reverse rates computed from thermochemistry, multiply the forward
712 // rate coefficients by the reciprocals of the equilibrium constants
713 for (size_t i = 0; i < nReactions(); ++i) {
714 rop[i] *= m_rkcn[i];
715 }
716}
717
718}
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.
A kinetics manager for heterogeneous reaction mechanisms.
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.
void addPhase(ThermoPhase &thermo) override
Eigen::SparseMatrix< double > revRatesOfProgress_ddCi() override
Calculate derivatives for forward rates-of-progress with respect to species concentration at constant...
size_t m_nDim
Number of dimensions of reacting phase (2 for InterfaceKinetics, 1 for EdgeKinetics)
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.
InterfaceKinetics()=default
Constructor.
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.
size_t reactionPhaseIndex() const
Phase where the reactions occur.
Definition Kinetics.h:235
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:254
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
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
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
vector< double > m_rkcn
Reciprocal of the equilibrium constant in concentration units.
Definition Kinetics.h:1605
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
size_t m_kk
The number of species in all of the phases that participate in this kinetics mechanism.
Definition Kinetics.h:1542
virtual void getReactionDelta(const double *property, double *deltaProperty) const
Change in species properties.
Definition Kinetics.cpp:392
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 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 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
StoichManagerN m_revProductStoich
Stoichiometry manager for the products of reversible reactions.
Definition Kinetics.h:1531
size_t nReactions() const
Number of reactions in the reaction mechanism.
Definition Kinetics.h:153
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
virtual void setMultiplier(size_t i, double f)
Set the multiplier for reaction i to f.
Definition Kinetics.h:1443
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
virtual void resizeSpecies()
Resize arrays with sizes that depend on the total number of species.
Definition Kinetics.cpp:653
size_t nTotalSpecies() const
The total number of species in all phases participating in the kinetics mechanism.
Definition Kinetics.h:266
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Definition Kinetics.cpp:435
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
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:511
virtual void getConcentrations(double *const c) const
Get the species concentrations (kmol/m^3).
Definition Phase.cpp:595
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
double temperature() const
Temperature (K).
Definition Phase.h:662
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:638
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.
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
Definition SurfPhase.h:98
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.
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.
virtual void getActivityConcentrations(double *c) const
This method returns an array of generalized concentrations.
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: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
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...