Cantera  3.2.0a2
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_rateHandlers) {
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_rateHandlers) {
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 for (auto& rates : m_rateHandlers) {
222 rates->modifyRateConstants(m_ropf.data(), m_ropr.data());
223 }
224
225 // multiply ropf by the activity concentration reaction orders to obtain
226 // the forward rates of progress.
227 m_reactantStoich.multiply(m_actConc.data(), m_ropf.data());
228
229 // For reversible reactions, multiply ropr by the activity concentration
230 // products
231 m_revProductStoich.multiply(m_actConc.data(), m_ropr.data());
232
233 for (size_t j = 0; j != nReactions(); ++j) {
234 m_ropnet[j] = m_ropf[j] - m_ropr[j];
235 }
236
237 // For reactions involving multiple phases, we must check that the phase
238 // being consumed actually exists. This is particularly important for phases
239 // that are stoichiometric phases containing one species with a unity
240 // activity
241 if (m_phaseExistsCheck) {
242 for (size_t j = 0; j != nReactions(); ++j) {
243 if ((m_ropr[j] > m_ropf[j]) && (m_ropr[j] > 0.0)) {
244 for (size_t p = 0; p < nPhases(); p++) {
245 if (m_rxnPhaseIsProduct[j][p] && !m_phaseExists[p]) {
246 m_ropnet[j] = 0.0;
247 m_ropr[j] = m_ropf[j];
248 if (m_ropf[j] > 0.0) {
249 for (size_t rp = 0; rp < nPhases(); rp++) {
250 if (m_rxnPhaseIsReactant[j][rp] && !m_phaseExists[rp]) {
251 m_ropnet[j] = 0.0;
252 m_ropr[j] = m_ropf[j] = 0.0;
253 }
254 }
255 }
256 }
257 if (m_rxnPhaseIsReactant[j][p] && !m_phaseIsStable[p]) {
258 m_ropnet[j] = 0.0;
259 m_ropr[j] = m_ropf[j];
260 }
261 }
262 } else if ((m_ropf[j] > m_ropr[j]) && (m_ropf[j] > 0.0)) {
263 for (size_t p = 0; p < nPhases(); p++) {
264 if (m_rxnPhaseIsReactant[j][p] && !m_phaseExists[p]) {
265 m_ropnet[j] = 0.0;
266 m_ropf[j] = m_ropr[j];
267 if (m_ropf[j] > 0.0) {
268 for (size_t rp = 0; rp < nPhases(); rp++) {
269 if (m_rxnPhaseIsProduct[j][rp] && !m_phaseExists[rp]) {
270 m_ropnet[j] = 0.0;
271 m_ropf[j] = m_ropr[j] = 0.0;
272 }
273 }
274 }
275 }
276 if (m_rxnPhaseIsProduct[j][p] && !m_phaseIsStable[p]) {
277 m_ropnet[j] = 0.0;
278 m_ropf[j] = m_ropr[j];
279 }
280 }
281 }
282 }
283 }
284 m_ROP_ok = true;
285}
286
288{
289 // Get the chemical potentials of the species in the all of the phases used
290 // in the kinetics mechanism
291 for (size_t n = 0; n < nPhases(); n++) {
292 m_thermo[n]->getChemPotentials(m_mu.data() + m_start[n]);
293 }
294
295 // Use the stoichiometric manager to find deltaG for each reaction.
296 getReactionDelta(m_mu.data(), m_rbuf.data());
297 if (deltaG != 0 && (m_rbuf.data() != deltaG)) {
298 for (size_t j = 0; j < nReactions(); ++j) {
299 deltaG[j] = m_rbuf[j];
300 }
301 }
302}
303
305{
306 // Get the chemical potentials of the species
307 for (size_t n = 0; n < nPhases(); n++) {
309 }
310
311 // Use the stoichiometric manager to find deltaG for each reaction.
312 getReactionDelta(m_grt.data(), deltaM);
313}
314
316{
317 // Get the partial molar enthalpy of all species
318 for (size_t n = 0; n < nPhases(); n++) {
320 }
321
322 // Use the stoichiometric manager to find deltaH for each reaction.
323 getReactionDelta(m_grt.data(), deltaH);
324}
325
327{
328 // Get the partial molar entropy of all species in all of the phases
329 for (size_t n = 0; n < nPhases(); n++) {
331 }
332
333 // Use the stoichiometric manager to find deltaS for each reaction.
334 getReactionDelta(m_grt.data(), deltaS);
335}
336
338{
339 // Get the standard state chemical potentials of the species. This is the
340 // array of chemical potentials at unit activity We define these here as the
341 // chemical potentials of the pure species at the temperature and pressure
342 // of the solution.
343 for (size_t n = 0; n < nPhases(); n++) {
345 }
346
347 // Use the stoichiometric manager to find deltaG for each reaction.
348 getReactionDelta(m_mu0.data(), deltaGSS);
349}
350
352{
353 // Get the standard state enthalpies of the species. This is the array of
354 // chemical potentials at unit activity We define these here as the
355 // enthalpies of the pure species at the temperature and pressure of the
356 // solution.
357 for (size_t n = 0; n < nPhases(); n++) {
358 thermo(n).getEnthalpy_RT(m_grt.data() + m_start[n]);
359 }
360 for (size_t k = 0; k < m_kk; k++) {
361 m_grt[k] *= thermo(0).RT();
362 }
363
364 // Use the stoichiometric manager to find deltaH for each reaction.
365 getReactionDelta(m_grt.data(), deltaH);
366}
367
369{
370 // Get the standard state entropy of the species. We define these here as
371 // the entropies of the pure species at the temperature and pressure of the
372 // solution.
373 for (size_t n = 0; n < nPhases(); n++) {
374 thermo(n).getEntropy_R(m_grt.data() + m_start[n]);
375 }
376 for (size_t k = 0; k < m_kk; k++) {
377 m_grt[k] *= GasConstant;
378 }
379
380 // Use the stoichiometric manager to find deltaS for each reaction.
381 getReactionDelta(m_grt.data(), deltaS);
382}
383
384bool InterfaceKinetics::addReaction(shared_ptr<Reaction> r_base, bool resize)
385{
386 if (!m_surf) {
387 init();
388 }
389
390 size_t i = nReactions();
391 shared_ptr<ReactionRate> rate = r_base->rate();
392 if (rate) {
393 rate->setContext(*r_base, *this);
394 }
395
396 bool added = Kinetics::addReaction(r_base, resize);
397 if (!added) {
398 return false;
399 }
400
401 m_rxnPhaseIsReactant.emplace_back(nPhases(), false);
402 m_rxnPhaseIsProduct.emplace_back(nPhases(), false);
403
404 for (const auto& [name, stoich] : r_base->reactants) {
405 size_t k = kineticsSpeciesIndex(name);
406 size_t p = speciesPhaseIndex(k);
407 m_rxnPhaseIsReactant[i][p] = true;
408 }
409 for (const auto& [name, stoich] : r_base->products) {
410 size_t k = kineticsSpeciesIndex(name);
411 size_t p = speciesPhaseIndex(k);
412 m_rxnPhaseIsProduct[i][p] = true;
413 }
414
415 // Set index of rate to number of reaction within kinetics
416 rate->setRateIndex(nReactions() - 1);
417
418 string rtype = rate->subType();
419 if (rtype == "") {
420 rtype = rate->type();
421 }
422
423 // If necessary, add new interface MultiRate evaluator
424 if (m_rateTypes.find(rtype) == m_rateTypes.end()) {
425 m_rateTypes[rtype] = m_rateHandlers.size();
426 m_rateHandlers.push_back(rate->newMultiRate());
427 m_rateHandlers.back()->resize(m_kk, nReactions(), nPhases());
428 }
429
430 // Add reaction rate to evaluator
431 size_t index = m_rateTypes[rtype];
432 m_rateHandlers[index]->add(nReactions() - 1, *rate);
433
434 // Set flag for coverage dependence to true
435 if (rate->compositionDependent()) {
437 }
438
439 // Set flag for electrochemistry to true
440 if (r_base->usesElectrochemistry(*this)) {
442 }
443
444 return true;
445}
446
447void InterfaceKinetics::modifyReaction(size_t i, shared_ptr<Reaction> r_base)
448{
449 Kinetics::modifyReaction(i, r_base);
450
451 shared_ptr<ReactionRate> rate = r_base->rate();
452 rate->setRateIndex(i);
453 rate->setContext(*r_base, *this);
454
455 string rtype = rate->subType();
456 if (rtype == "") {
457 rtype = rate->type();
458 }
459
460 // Ensure that interface MultiRate evaluator is available
461 if (!m_rateTypes.count(rtype)) {
462 throw CanteraError("InterfaceKinetics::modifyReaction",
463 "Interface evaluator not available for type '{}'.", rtype);
464 }
465 // Replace reaction rate evaluator
466 size_t index = m_rateTypes[rate->type()];
467 m_rateHandlers[index]->replace(i, *rate);
468
469 // Invalidate cached data
470 m_redo_rates = true;
471 m_temp += 0.1;
472}
473
474void InterfaceKinetics::setMultiplier(size_t i, double f)
475{
477 m_ROP_ok = false;
478}
479
480void InterfaceKinetics::setIOFlag(int ioFlag)
481{
482 m_ioFlag = ioFlag;
483 if (m_integrator) {
484 m_integrator->setIOFlag(ioFlag);
485 }
486}
487
488void InterfaceKinetics::addThermo(shared_ptr<ThermoPhase> thermo)
489{
491 m_phaseExists.push_back(true);
492 m_phaseIsStable.push_back(true);
493}
494
496{
497 if (thermo(0).nDim() > 2) {
498 throw CanteraError("InterfaceKinetics::init", "no interface phase is present.");
499 }
500}
501
503{
504 size_t kOld = m_kk;
506 if (m_kk != kOld && nReactions()) {
507 throw CanteraError("InterfaceKinetics::resizeSpecies", "Cannot add"
508 " species to InterfaceKinetics after reactions have been added.");
509 }
510 m_actConc.resize(m_kk);
511 m_conc.resize(m_kk);
512 m_mu0.resize(m_kk);
513 m_mu.resize(m_kk);
514 m_mu0_Kc.resize(m_kk);
515 m_grt.resize(m_kk);
516 m_phi.resize(nPhases(), 0.0);
517}
518
519void InterfaceKinetics::advanceCoverages(double tstep, double rtol, double atol,
520 double maxStepSize, size_t maxSteps, size_t maxErrTestFails)
521{
522 if (m_integrator == 0) {
523 vector<InterfaceKinetics*> k{this};
525 }
526 m_integrator->setTolerances(rtol, atol);
527 m_integrator->setMaxStepSize(maxStepSize);
528 m_integrator->setMaxSteps(maxSteps);
529 m_integrator->setMaxErrTestFails(maxErrTestFails);
530 m_integrator->integrate(0.0, tstep);
531 delete m_integrator;
532 m_integrator = 0;
533}
534
536 int ifuncOverride, double timeScaleOverride)
537{
538 // create our own solver object
539 if (m_integrator == 0) {
540 vector<InterfaceKinetics*> k{this};
543 }
544 m_integrator->setIOFlag(m_ioFlag);
545 // New direct method to go here
546 m_integrator->solvePseudoSteadyStateProblem(ifuncOverride, timeScaleOverride);
547}
548
549void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists)
550{
551 checkPhaseIndex(iphase);
552 if (exists) {
553 if (!m_phaseExists[iphase]) {
556 m_phaseExists[iphase] = true;
557 }
558 m_phaseIsStable[iphase] = true;
559 } else {
560 if (m_phaseExists[iphase]) {
562 m_phaseExists[iphase] = false;
563 }
564 m_phaseIsStable[iphase] = false;
565 }
566}
567
568int InterfaceKinetics::phaseExistence(const size_t iphase) const
569{
570 checkPhaseIndex(iphase);
571 return m_phaseExists[iphase];
572}
573
574int InterfaceKinetics::phaseStability(const size_t iphase) const
575{
576 checkPhaseIndex(iphase);
577 return m_phaseIsStable[iphase];
578}
579
580void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStable)
581{
582 checkPhaseIndex(iphase);
583 if (isStable) {
584 m_phaseIsStable[iphase] = true;
585 } else {
586 m_phaseIsStable[iphase] = false;
587 }
588}
589
590double InterfaceKinetics::interfaceCurrent(const size_t iphase)
591{
592 vector<double> charges(m_kk, 0.0);
593 vector<double> netProdRates(m_kk, 0.0);
594 double dotProduct = 0.0;
595
596 thermo(iphase).getCharges(charges.data());
597 getNetProductionRates(netProdRates.data());
598
599 for (size_t k = 0; k < thermo(iphase).nSpecies(); k++)
600 {
601 dotProduct += charges[k] * netProdRates[m_start[iphase] + k];
602 }
603
604 return dotProduct * Faraday;
605}
606
608{
609 // check derivatives are valid
610 assertDerivativesValid("InterfaceKinetics::fwdRatesOfProgress_ddCi");
611 // forward reaction rate coefficients
612 vector<double>& rop_rates = m_rbuf0;
613 getFwdRateConstants(rop_rates.data());
615}
616
618{
619 // check derivatives are valid
620 assertDerivativesValid("InterfaceKinetics::revRatesOfProgress_ddCi");
621 // reverse reaction rate coefficients
622 vector<double>& rop_rates = m_rbuf0;
623 getFwdRateConstants(rop_rates.data());
624 applyEquilibriumConstants(rop_rates.data());
626}
627
629{
630 // check derivatives are valid
631 assertDerivativesValid("InterfaceKinetics::netRatesOfProgress_ddCi");
632 // forward reaction rate coefficients
633 vector<double>& rop_rates = m_rbuf0;
634 getFwdRateConstants(rop_rates.data());
635 Eigen::SparseMatrix<double> jac = calculateCompositionDerivatives(m_reactantStoich,
636 rop_rates);
637
638 // reverse reaction rate coefficients
639 applyEquilibriumConstants(rop_rates.data());
641}
642
644{
645 bool force = settings.empty();
646 if (force || settings.hasKey("skip-coverage-dependence")) {
647 m_jac_skip_coverage_dependence = settings.getBool("skip-coverage-dependence",
648 false);
649 }
650 if (force || settings.hasKey("skip-electrochemistry")) {
651 m_jac_skip_electrochemistry = settings.getBool("skip-electrochemistry",
652 false);
653 }
654 if (force || settings.hasKey("rtol-delta")) {
655 m_jac_rtol_delta = settings.getDouble("rtol-delta", 1e-8);
656 }
657}
658
660{
661 settings["skip-coverage-dependence"] = m_jac_skip_electrochemistry;
662 settings["skip-electrochemistry"] = m_jac_skip_coverage_dependence;
663 settings["rtol-delta"] = m_jac_rtol_delta;
664}
665
667 StoichManagerN& stoich, const vector<double>& in)
668{
669 vector<double>& outV = m_rbuf1;
670 // derivatives handled by StoichManagerN
671 copy(in.begin(), in.end(), outV.begin());
672 return stoich.derivatives(m_actConc.data(), outV.data());
673}
674
676{
678 throw NotImplementedError(name, "Coverage-dependent reactions not supported.");
680 throw NotImplementedError(name, "Electrochemical reactions not supported.");
681 }
682}
683
685{
686 // For reverse rates computed from thermochemistry, multiply the forward
687 // rate coefficients by the reciprocals of the equilibrium constants
688 for (size_t i = 0; i < nReactions(); ++i) {
689 rop[i] *= m_rkcn[i];
690 }
691}
692
693}
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:431
double getDouble(const string &key, double default_) const
If key exists, return it as a double, otherwise return default_.
Definition AnyMap.cpp:1580
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition AnyMap.cpp:1477
bool empty() const
Return boolean indicating whether AnyMap is empty.
Definition AnyMap.cpp:1468
bool getBool(const string &key, bool default_) const
If key exists, return it as a bool, otherwise return default_.
Definition AnyMap.cpp:1575
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.
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.
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< 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.
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:34
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:59
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
Definition Kinetics.h:240
vector< double > m_perturb
Vector of perturbation factors for each reaction's rate of progress vector.
Definition Kinetics.h:1483
vector< double > m_ropf
Forward rate-of-progress for each reaction.
Definition Kinetics.h:1527
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:377
vector< unique_ptr< MultiRateBase > > m_rateHandlers
Vector of rate handlers.
Definition Kinetics.h:1454
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:1504
vector< double > m_rkcn
Reciprocal of the equilibrium constant in concentration units.
Definition Kinetics.h:1524
size_t m_kk
The number of species in all of the phases that participate in this kinetics mechanism.
Definition Kinetics.h:1479
virtual void getReactionDelta(const double *property, double *deltaProperty) const
Change in species properties.
Definition Kinetics.cpp:368
vector< double > m_dH
The enthalpy change for each reaction to calculate Blowers-Masel rates.
Definition Kinetics.h:1539
vector< double > m_ropr
Reverse rate-of-progress for each reaction.
Definition Kinetics.h:1530
size_t nPhases() const
The number of phases participating in the reaction mechanism.
Definition Kinetics.h:185
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:1498
virtual bool addReaction(shared_ptr< Reaction > r, bool resize=true)
Add a single reaction to the mechanism.
Definition Kinetics.cpp:619
vector< size_t > m_revindex
Indices of reversible reactions.
Definition Kinetics.h:1535
virtual void modifyReaction(size_t i, shared_ptr< Reaction > rNew)
Modify the rate expression associated with a reaction.
Definition Kinetics.cpp:712
map< string, size_t > m_rateTypes
Mapping of rate handlers.
Definition Kinetics.h:1455
vector< double > m_ropnet
Net rate-of-progress for each reaction.
Definition Kinetics.h:1533
virtual void addThermo(shared_ptr< ThermoPhase > thermo)
Add a phase to the kinetics manager object.
Definition Kinetics.cpp:568
vector< double > m_rbuf
Buffer used for storage of intermediate reaction-specific results.
Definition Kinetics.h:1542
StoichManagerN m_revProductStoich
Stoichiometry manager for the products of reversible reactions.
Definition Kinetics.h:1471
size_t nReactions() const
Number of reactions in the reaction mechanism.
Definition Kinetics.h:153
vector< size_t > m_irrev
Indices of irreversible reactions.
Definition Kinetics.h:1536
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition Kinetics.h:274
virtual void setMultiplier(size_t i, double f)
Set the multiplier for reaction i to f.
Definition Kinetics.h:1376
vector< double > m_rfn
Forward rate constant for each reaction.
Definition Kinetics.h:1518
StoichManagerN m_reactantStoich
Stoichiometry manager for the reactants for each reaction.
Definition Kinetics.h:1465
virtual void resizeSpecies()
Resize arrays with sizes that depend on the total number of species.
Definition Kinetics.cpp:607
size_t nTotalSpecies() const
The total number of species in all phases participating in the kinetics mechanism.
Definition Kinetics.h:252
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Definition Kinetics.cpp:411
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:329
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:232
double temperature() const
Temperature (K).
Definition Phase.h:563
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:539
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:595
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...