Cantera  3.1.0a1
Loading...
Searching...
No Matches
BulkKinetics.cpp
1// This file is part of Cantera. See License.txt in the top-level directory or
2// at https://cantera.org/license.txt for license and copyright information.
3
7
8namespace Cantera
9{
10
11BulkKinetics::BulkKinetics() {
12 setDerivativeSettings(AnyMap()); // use default settings
13}
14
16 return std::find(m_revindex.begin(), m_revindex.end(), i) < m_revindex.end();
17}
18
19bool BulkKinetics::addReaction(shared_ptr<Reaction> r, bool resize)
20{
21 bool added = Kinetics::addReaction(r, resize);
22 if (!added) {
23 // undeclared species, etc.
24 return false;
25 }
26 double dn = 0.0;
27 for (const auto& [name, stoich] : r->products) {
28 dn += stoich;
29 }
30 for (const auto& [name, stoich] : r->reactants) {
31 dn -= stoich;
32 }
33
34 m_dn.push_back(dn);
35
36 if (r->reversible) {
37 m_revindex.push_back(nReactions()-1);
38 } else {
39 m_irrev.push_back(nReactions()-1);
40 }
41
42 shared_ptr<ReactionRate> rate = r->rate();
43 string rtype = rate->subType();
44 if (rtype == "") {
45 rtype = rate->type();
46 }
47
48 // If necessary, add new MultiRate evaluator
49 if (m_bulk_types.find(rtype) == m_bulk_types.end()) {
50 m_bulk_types[rtype] = m_bulk_rates.size();
51 m_bulk_rates.push_back(rate->newMultiRate());
52 m_bulk_rates.back()->resize(m_kk, nReactions(), nPhases());
53 }
54
55 // Set index of rate to number of reaction within kinetics
56 rate->setRateIndex(nReactions() - 1);
57 rate->setContext(*r, *this);
58
59 // Add reaction rate to evaluator
60 size_t index = m_bulk_types[rtype];
61 m_bulk_rates[index]->add(nReactions() - 1, *rate);
62
63 // Add reaction to third-body evaluator
64 if (r->thirdBody() != nullptr) {
65 addThirdBody(r);
66 }
67
68 m_concm.push_back(NAN);
69 m_ready = resize;
70 return true;
71}
72
73void BulkKinetics::addThirdBody(shared_ptr<Reaction> r)
74{
75 map<size_t, double> efficiencies;
76 for (const auto& [name, efficiency] : r->thirdBody()->efficiencies) {
77 size_t k = kineticsSpeciesIndex(name);
78 if (k != npos) {
79 efficiencies[k] = efficiency;
80 } else if (!m_skipUndeclaredThirdBodies) {
81 throw CanteraError("BulkKinetics::addThirdBody", "Found third-body"
82 " efficiency for undefined species '{}' while adding reaction '{}'",
83 name, r->equation());
84 } else {
86 }
87 }
88 m_multi_concm.install(nReactions() - 1, efficiencies,
89 r->thirdBody()->default_efficiency,
90 r->thirdBody()->mass_action);
91}
92
93void BulkKinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
94{
95 // operations common to all reaction types
97
98 shared_ptr<ReactionRate> rate = rNew->rate();
99 string rtype = rate->subType();
100 if (rtype == "") {
101 rtype = rate->type();
102 }
103
104 // Ensure that MultiRate evaluator is available
105 if (m_bulk_types.find(rtype) == m_bulk_types.end()) {
106 throw CanteraError("BulkKinetics::modifyReaction",
107 "Evaluator not available for type '{}'.", rtype);
108 }
109
110 // Replace reaction rate to evaluator
111 size_t index = m_bulk_types[rtype];
112 rate->setRateIndex(i);
113 rate->setContext(*rNew, *this);
114 m_bulk_rates[index]->replace(i, *rate);
115 invalidateCache();
116}
117
119{
121 m_act_conc.resize(m_kk);
122 m_phys_conc.resize(m_kk);
123 m_grt.resize(m_kk);
124 for (auto& rates : m_bulk_rates) {
125 rates->resize(m_kk, nReactions(), nPhases());
126 }
127}
128
130{
132 m_rbuf0.resize(nReactions());
133 m_rbuf1.resize(nReactions());
134 m_rbuf2.resize(nReactions());
135 m_kf0.resize(nReactions());
136 m_sbuf0.resize(nTotalSpecies());
137 m_state.resize(thermo().stateSize());
139 for (auto& rates : m_bulk_rates) {
140 rates->resize(nTotalSpecies(), nReactions(), nPhases());
141 // @todo ensure that ReactionData are updated; calling rates->update
142 // blocks correct behavior in update_rates_T
143 // and running updateROP() is premature
144 }
145}
146
147void BulkKinetics::setMultiplier(size_t i, double f)
148{
150 m_ROP_ok = false;
151}
152
153void BulkKinetics::invalidateCache()
154{
155 Kinetics::invalidateCache();
156 m_ROP_ok = false;
157}
158
160{
161 updateROP();
162 copy(m_rfn.begin(), m_rfn.end(), kfwd);
164 processThirdBodies(kfwd);
165 }
166}
167
169{
170 updateROP();
171
172 vector<double>& delta_gibbs0 = m_rbuf0;
173 fill(delta_gibbs0.begin(), delta_gibbs0.end(), 0.0);
174
175 // compute Delta G^0 for all reactions
176 getReactionDelta(m_grt.data(), delta_gibbs0.data());
177
178 double rrt = 1.0 / thermo().RT();
179 double logStandConc = log(thermo().standardConcentration());
180 for (size_t i = 0; i < nReactions(); i++) {
181 kc[i] = exp(-delta_gibbs0[i] * rrt + m_dn[i] * logStandConc);
182 }
183}
184
185void BulkKinetics::getRevRateConstants(double* krev, bool doIrreversible)
186{
187 // go get the forward rate constants. -> note, we don't really care about
188 // speed or redundancy in these informational routines.
190
191 if (doIrreversible) {
193 for (size_t i = 0; i < nReactions(); i++) {
194 krev[i] /= m_rbuf0[i];
195 }
196 } else {
197 // m_rkcn[] is zero for irreversible reactions
198 for (size_t i = 0; i < nReactions(); i++) {
199 krev[i] *= m_rkcn[i];
200 }
201 }
202}
203
204void BulkKinetics::getDeltaGibbs(double* deltaG)
205{
206 // Get the chemical potentials for each species
207 thermo().getChemPotentials(m_sbuf0.data());
208 // Use the stoichiometric manager to find deltaG for each reaction.
209 getReactionDelta(m_sbuf0.data(), deltaG);
210}
211
213{
214 // Get the partial molar enthalpy for each species
215 thermo().getPartialMolarEnthalpies(m_sbuf0.data());
216 // Use the stoichiometric manager to find deltaH for each reaction.
217 getReactionDelta(m_sbuf0.data(), deltaH);
218}
219
221{
222 // Get the partial molar entropy for each species
223 thermo().getPartialMolarEntropies(m_sbuf0.data());
224 // Use the stoichiometric manager to find deltaS for each reaction.
225 getReactionDelta(m_sbuf0.data(), deltaS);
226}
227
229{
230 // Get the standard state chemical potentials of the species. This is the
231 // array of chemical potentials at unit activity. We define these here as
232 // the chemical potentials of the pure species at the temperature and
233 // pressure of the solution.
234 thermo().getStandardChemPotentials(m_sbuf0.data());
235 // Use the stoichiometric manager to find deltaG for each reaction.
236 getReactionDelta(m_sbuf0.data(), deltaG);
237}
238
240{
241 // Get the standard state enthalpies of the species.
242 thermo().getEnthalpy_RT(m_sbuf0.data());
243 for (size_t k = 0; k < m_kk; k++) {
244 m_sbuf0[k] *= thermo().RT();
245 }
246 // Use the stoichiometric manager to find deltaH for each reaction.
247 getReactionDelta(m_sbuf0.data(), deltaH);
248}
249
251{
252 // Get the standard state entropy of the species. We define these here as
253 // the entropies of the pure species at the temperature and pressure of the
254 // solution.
255 thermo().getEntropy_R(m_sbuf0.data());
256 for (size_t k = 0; k < m_kk; k++) {
257 m_sbuf0[k] *= GasConstant;
258 }
259 // Use the stoichiometric manager to find deltaS for each reaction.
260 getReactionDelta(m_sbuf0.data(), deltaS);
261}
262
264{
265 settings["skip-third-bodies"] = m_jac_skip_third_bodies;
266 settings["skip-falloff"] = m_jac_skip_falloff;
267 settings["rtol-delta"] = m_jac_rtol_delta;
268}
269
271{
272 bool force = settings.empty();
273 if (force || settings.hasKey("skip-third-bodies")) {
274 m_jac_skip_third_bodies = settings.getBool("skip-third-bodies", false);
275 }
276 if (force || settings.hasKey("skip-falloff")) {
277 m_jac_skip_falloff = settings.getBool("skip-falloff", false);
278 }
279 if (force || settings.hasKey("rtol-delta")) {
280 m_jac_rtol_delta = settings.getDouble("rtol-delta", 1e-8);
281 }
282}
283
285{
286 assertDerivativesValid("BulkKinetics::getFwdRateConstants_ddT");
287 updateROP();
288 process_ddT(m_rfn, dkfwd);
289}
290
292{
293 assertDerivativesValid("BulkKinetics::getFwdRatesOfProgress_ddT");
294 updateROP();
295 process_ddT(m_ropf, drop);
296}
297
299{
300 assertDerivativesValid("BulkKinetics::getRevRatesOfProgress_ddT");
301 updateROP();
302 process_ddT(m_ropr, drop);
303 Eigen::Map<Eigen::VectorXd> dRevRop(drop, nReactions());
304
305 // reverse rop times scaled inverse equilibrium constant derivatives
306 Eigen::Map<Eigen::VectorXd> dRevRop2(m_rbuf2.data(), nReactions());
307 copy(m_ropr.begin(), m_ropr.end(), m_rbuf2.begin());
308 applyEquilibriumConstants_ddT(dRevRop2.data());
309 dRevRop += dRevRop2;
310}
311
313{
314 assertDerivativesValid("BulkKinetics::getNetRatesOfProgress_ddT");
315 updateROP();
316 process_ddT(m_ropnet, drop);
317 Eigen::Map<Eigen::VectorXd> dNetRop(drop, nReactions());
318
319 // reverse rop times scaled inverse equilibrium constant derivatives
320 Eigen::Map<Eigen::VectorXd> dNetRop2(m_rbuf2.data(), nReactions());
321 copy(m_ropr.begin(), m_ropr.end(), m_rbuf2.begin());
322 applyEquilibriumConstants_ddT(dNetRop2.data());
323 dNetRop -= dNetRop2;
324}
325
327{
328 assertDerivativesValid("BulkKinetics::getFwdRateConstants_ddP");
329 updateROP();
330 process_ddP(m_rfn, dkfwd);
331}
332
334{
335 assertDerivativesValid("BulkKinetics::getFwdRatesOfProgress_ddP");
336 updateROP();
337 process_ddP(m_ropf, drop);
338}
339
341{
342 assertDerivativesValid("BulkKinetics::getRevRatesOfProgress_ddP");
343 updateROP();
344 process_ddP(m_ropr, drop);
345}
346
348{
349 assertDerivativesValid("BulkKinetics::getNetRatesOfProgress_ddP");
350 updateROP();
351 process_ddP(m_ropnet, drop);
352}
353
355{
356 assertDerivativesValid("BulkKinetics::getFwdRateConstants_ddC");
357 updateROP();
358 process_ddC(m_reactantStoich, m_rfn, dkfwd, false);
359}
360
362{
363 assertDerivativesValid("BulkKinetics::getFwdRatesOfProgress_ddC");
364 updateROP();
366}
367
369{
370 assertDerivativesValid("BulkKinetics::getRevRatesOfProgress_ddC");
371 updateROP();
373}
374
376{
377 assertDerivativesValid("BulkKinetics::getNetRatesOfProgress_ddC");
378 updateROP();
380 Eigen::Map<Eigen::VectorXd> dNetRop(drop, nReactions());
381
382 process_ddC(m_revProductStoich, m_ropr, m_rbuf2.data());
383 Eigen::Map<Eigen::VectorXd> dNetRop2(m_rbuf2.data(), nReactions());
384 dNetRop -= dNetRop2;
385}
386
387Eigen::SparseMatrix<double> BulkKinetics::fwdRatesOfProgress_ddX()
388{
389 assertDerivativesValid("BulkKinetics::fwdRatesOfProgress_ddX");
390
391 // forward reaction rate coefficients
392 vector<double>& rop_rates = m_rbuf0;
393 getFwdRateConstants(rop_rates.data());
395}
396
397Eigen::SparseMatrix<double> BulkKinetics::revRatesOfProgress_ddX()
398{
399 assertDerivativesValid("BulkKinetics::revRatesOfProgress_ddX");
400
401 // reverse reaction rate coefficients
402 vector<double>& rop_rates = m_rbuf0;
403 getFwdRateConstants(rop_rates.data());
404 applyEquilibriumConstants(rop_rates.data());
406}
407
408Eigen::SparseMatrix<double> BulkKinetics::netRatesOfProgress_ddX()
409{
410 assertDerivativesValid("BulkKinetics::netRatesOfProgress_ddX");
411
412 // forward reaction rate coefficients
413 vector<double>& rop_rates = m_rbuf0;
414 getFwdRateConstants(rop_rates.data());
416
417 // reverse reaction rate coefficients
418 applyEquilibriumConstants(rop_rates.data());
420}
421
422Eigen::SparseMatrix<double> BulkKinetics::fwdRatesOfProgress_ddCi()
423{
424 assertDerivativesValid("BulkKinetics::fwdRatesOfProgress_ddCi");
425
426 // forward reaction rate coefficients
427 vector<double>& rop_rates = m_rbuf0;
428 getFwdRateConstants(rop_rates.data());
429 return calculateCompositionDerivatives(m_reactantStoich, rop_rates, false);
430}
431
432Eigen::SparseMatrix<double> BulkKinetics::revRatesOfProgress_ddCi()
433{
434 assertDerivativesValid("BulkKinetics::revRatesOfProgress_ddCi");
435
436 // reverse reaction rate coefficients
437 vector<double>& rop_rates = m_rbuf0;
438 getFwdRateConstants(rop_rates.data());
439 applyEquilibriumConstants(rop_rates.data());
440 return calculateCompositionDerivatives(m_revProductStoich, rop_rates, false);
441}
442
443Eigen::SparseMatrix<double> BulkKinetics::netRatesOfProgress_ddCi()
444{
445 assertDerivativesValid("BulkKinetics::netRatesOfProgress_ddCi");
446
447 // forward reaction rate coefficients
448 vector<double>& rop_rates = m_rbuf0;
449 getFwdRateConstants(rop_rates.data());
450 auto jac = calculateCompositionDerivatives(m_reactantStoich, rop_rates, false);
451
452 // reverse reaction rate coefficients
453 applyEquilibriumConstants(rop_rates.data());
454 return jac - calculateCompositionDerivatives(m_revProductStoich, rop_rates, false);
455}
456
457void BulkKinetics::updateROP()
458{
459 static const int cacheId = m_cache.getId();
460 CachedScalar last = m_cache.getScalar(cacheId);
461 double T = thermo().temperature();
462 double rho = thermo().density();
463 int statenum = thermo().stateMFNumber();
464
465 if (last.state1 != T || last.state2 != rho) {
466 // Update properties that are independent of the composition
468 fill(m_delta_gibbs0.begin(), m_delta_gibbs0.end(), 0.0);
469 double logStandConc = log(thermo().standardConcentration());
470
471 // compute Delta G^0 for all reversible reactions
473
474 double rrt = 1.0 / thermo().RT();
475 for (size_t i = 0; i < m_revindex.size(); i++) {
476 size_t irxn = m_revindex[i];
477 m_rkcn[irxn] = std::min(
478 exp(m_delta_gibbs0[irxn] * rrt - m_dn[irxn] * logStandConc), BigNumber);
479 }
480
481 for (size_t i = 0; i != m_irrev.size(); ++i) {
482 m_rkcn[ m_irrev[i] ] = 0.0;
483 }
484 }
485
486 if (!last.validate(T, rho, statenum)) {
487 // Update terms dependent on species concentrations and temperature
490 double ctot = thermo().molarDensity();
491
492 // Third-body objects interacting with MultiRate evaluator
494
495 // loop over MultiRate evaluators for each reaction type
496 for (auto& rates : m_bulk_rates) {
497 bool changed = rates->update(thermo(), *this);
498 if (changed) {
499 rates->getRateConstants(m_kf0.data());
500 }
501 }
502 m_ROP_ok = false;
503 }
504
505 if (m_ROP_ok) {
506 // rates of progress are up-to-date only if both the thermodynamic state
507 // and m_perturb are unchanged
508 return;
509 }
510
511 // Scale the forward rate coefficient by the perturbation factor
512 for (size_t i = 0; i < nReactions(); ++i) {
513 m_rfn[i] = m_kf0[i] * m_perturb[i];
514 }
515
516 copy(m_rfn.begin(), m_rfn.end(), m_ropf.data());
518 copy(m_ropf.begin(), m_ropf.end(), m_ropr.begin());
519
520 // multiply ropf by concentration products
521 m_reactantStoich.multiply(m_act_conc.data(), m_ropf.data());
522
523 // for reversible reactions, multiply ropr by concentration products
525 m_revProductStoich.multiply(m_act_conc.data(), m_ropr.data());
526 for (size_t j = 0; j != nReactions(); ++j) {
527 m_ropnet[j] = m_ropf[j] - m_ropr[j];
528 }
529
530 for (size_t i = 0; i < m_rfn.size(); i++) {
531 AssertFinite(m_rfn[i], "BulkKinetics::updateROP",
532 "m_rfn[{}] is not finite.", i);
533 AssertFinite(m_ropf[i], "BulkKinetics::updateROP",
534 "m_ropf[{}] is not finite.", i);
535 AssertFinite(m_ropr[i], "BulkKinetics::updateROP",
536 "m_ropr[{}] is not finite.", i);
537 }
538 m_ROP_ok = true;
539}
540
542{
543 updateROP();
544 std::copy(m_concm.begin(), m_concm.end(), concm);
545}
546
548{
549 // reactions involving third body
550 if (!m_concm.empty()) {
551 m_multi_concm.multiply(rop, m_concm.data());
552 }
553}
554
556{
557 // For reverse rates computed from thermochemistry, multiply the forward
558 // rate coefficients by the reciprocals of the equilibrium constants
559 for (size_t i = 0; i < nReactions(); ++i) {
560 rop[i] *= m_rkcn[i];
561 }
562}
563
565{
566 double T = thermo().temperature();
567 double P = thermo().pressure();
568 double rrt = 1. / thermo().RT();
569
570 vector<double>& grt = m_sbuf0;
571 vector<double>& delta_gibbs0 = m_rbuf1;
572 fill(delta_gibbs0.begin(), delta_gibbs0.end(), 0.0);
573
574 // compute perturbed Delta G^0 for all reversible reactions
575 thermo().saveState(m_state);
576 thermo().setState_TP(T * (1. + m_jac_rtol_delta), P);
577 thermo().getStandardChemPotentials(grt.data());
578 getRevReactionDelta(grt.data(), delta_gibbs0.data());
579
580 // apply scaling for derivative of inverse equilibrium constant
581 double Tinv = 1. / T;
582 double rrt_dTinv = rrt * Tinv / m_jac_rtol_delta;
583 double rrtt = rrt * Tinv;
584 for (size_t i = 0; i < m_revindex.size(); i++) {
585 size_t irxn = m_revindex[i];
586 double factor = delta_gibbs0[irxn] - m_delta_gibbs0[irxn];
587 factor *= rrt_dTinv;
588 factor += m_dn[irxn] * Tinv - m_delta_gibbs0[irxn] * rrtt;
589 drkcn[irxn] *= factor;
590 }
591
592 for (size_t i = 0; i < m_irrev.size(); ++i) {
593 drkcn[m_irrev[i]] = 0.0;
594 }
595
596 thermo().restoreState(m_state);
597}
598
599void BulkKinetics::process_ddT(const vector<double>& in, double* drop)
600{
601 // apply temperature derivative
602 copy(in.begin(), in.end(), drop);
603 for (auto& rates : m_bulk_rates) {
604 rates->processRateConstants_ddT(drop, m_rfn.data(), m_jac_rtol_delta);
605 }
606}
607
608void BulkKinetics::process_ddP(const vector<double>& in, double* drop)
609{
610 // apply pressure derivative
611 copy(in.begin(), in.end(), drop);
612 for (auto& rates : m_bulk_rates) {
613 rates->processRateConstants_ddP(drop, m_rfn.data(), m_jac_rtol_delta);
614 }
615}
616
617void BulkKinetics::process_ddC(StoichManagerN& stoich, const vector<double>& in,
618 double* drop, bool mass_action)
619{
620 Eigen::Map<Eigen::VectorXd> out(drop, nReactions());
621 out.setZero();
622 double ctot_inv = 1. / thermo().molarDensity();
623
624 // derivatives due to concentrations in law of mass action
625 if (mass_action) {
626 stoich.scale(in.data(), out.data(), ctot_inv);
627 }
629 return;
630 }
631
632 // derivatives due to third-body colliders in law of mass action
633 Eigen::Map<Eigen::VectorXd> outM(m_rbuf1.data(), nReactions());
634 if (mass_action) {
635 outM.fill(0.);
636 m_multi_concm.scale(in.data(), outM.data(), ctot_inv);
637 out += outM;
638 }
639
640 // derivatives due to reaction rates depending on third-body colliders
641 if (!m_jac_skip_falloff) {
642 m_multi_concm.scaleM(in.data(), outM.data(), m_concm.data(), ctot_inv);
643 for (auto& rates : m_bulk_rates) {
644 // processing step assigns zeros to entries not dependent on M
645 rates->processRateConstants_ddM(
646 outM.data(), m_rfn.data(), m_jac_rtol_delta);
647 }
648 out += outM;
649 }
650}
651
653 StoichManagerN& stoich, const vector<double>& in, bool ddX)
654{
655 Eigen::SparseMatrix<double> out;
656 vector<double>& scaled = m_rbuf1;
657 vector<double>& outV = m_rbuf2;
658
659 // convert from concentration to mole fraction output
660 copy(in.begin(), in.end(), scaled.begin());
661 if (ddX) {
662 double ctot = thermo().molarDensity();
663 for (size_t i = 0; i < nReactions(); ++i) {
664 scaled[i] *= ctot;
665 }
666 }
667
668 // derivatives handled by StoichManagerN
669 copy(scaled.begin(), scaled.end(), outV.begin());
670 processThirdBodies(outV.data());
671 out = stoich.derivatives(m_act_conc.data(), outV.data());
673 return out;
674 }
675
676 // derivatives due to law of mass action
677 copy(scaled.begin(), scaled.end(), outV.begin());
678 stoich.multiply(m_act_conc.data(), outV.data());
679
680 // derivatives due to reaction rates depending on third-body colliders
681 if (!m_jac_skip_falloff) {
682 for (auto& rates : m_bulk_rates) {
683 // processing step does not modify entries not dependent on M
684 rates->processRateConstants_ddM(
685 outV.data(), m_rfn.data(), m_jac_rtol_delta, false);
686 }
687 }
688
689 // derivatives handled by ThirdBodyCalc
690 out += m_multi_concm.derivatives(outV.data());
691 return out;
692}
693
695{
696 if (!thermo().isIdeal()) {
697 throw NotImplementedError(name,
698 "Not supported for non-ideal ThermoPhase models.");
699 }
700}
701
702}
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:427
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
bool addReaction(shared_ptr< Reaction > r, bool resize=true) override
Add a single reaction to the mechanism.
void getNetRatesOfProgress_ddC(double *drop) override
Calculate derivatives for net rates-of-progress with respect to molar concentration at constant tempe...
Eigen::SparseMatrix< double > netRatesOfProgress_ddX() override
Calculate derivatives for net rates-of-progress with respect to species mole fractions at constant te...
void getNetRatesOfProgress_ddP(double *drop) override
Calculate derivatives for net rates-of-progress with respect to pressure at constant temperature,...
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.
void getDeltaSSGibbs(double *deltaG) override
Return the vector of values for the reaction standard state Gibbs free energy change.
map< string, size_t > m_bulk_types
Mapping of rate handlers.
Eigen::SparseMatrix< double > fwdRatesOfProgress_ddX() override
Calculate derivatives for forward rates-of-progress with respect to species mole fractions at constan...
void getDeltaGibbs(double *deltaG) override
Return the vector of values for the reaction Gibbs free energy change.
vector< double > m_phys_conc
Physical concentrations, as calculated by ThermoPhase::getConcentrations.
Eigen::SparseMatrix< double > revRatesOfProgress_ddX() override
Calculate derivatives for reverse rates-of-progress with respect to species mole fractions at constan...
bool isReversible(size_t i) override
True if reaction i has been declared to be reversible.
void getFwdRatesOfProgress_ddT(double *drop) override
Calculate derivatives for forward rates-of-progress with respect to temperature at constant pressure,...
void assertDerivativesValid(const string &name)
Helper function ensuring that all rate derivatives can be calculated.
void getFwdRatesOfProgress_ddP(double *drop) override
Calculate derivatives for forward rates-of-progress with respect to pressure at constant temperature,...
void getFwdRateConstants_ddT(double *dkfwd) override
Calculate derivatives for forward rate constants with respect to temperature at constant pressure,...
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.
vector< double > m_dn
Difference between the global reactants order and the global products order.
bool m_jac_skip_third_bodies
Derivative settings.
void getDeltaSSEnthalpy(double *deltaH) override
Return the vector of values for the change in the standard state enthalpies of reaction.
Eigen::SparseMatrix< double > calculateCompositionDerivatives(StoichManagerN &stoich, const vector< double > &in, bool ddX=true)
Process derivatives.
void resizeSpecies() override
Resize arrays with sizes that depend on the total number of species.
vector< double > m_grt
Standard chemical potentials for each species.
vector< double > m_act_conc
Activity concentrations, as calculated by ThermoPhase::getActivityConcentrations.
void getRevRateConstants(double *krev, bool doIrreversible=false) override
Return the reverse rate constants.
void getEquilibriumConstants(double *kc) override
Return a vector of Equilibrium constants.
void getNetRatesOfProgress_ddT(double *drop) override
Calculate derivatives for net rates-of-progress with respect to temperature at constant pressure,...
vector< size_t > m_revindex
Indices of reversible reactions.
void processThirdBodies(double *rop)
Multiply rate with third-body collider concentrations.
void applyEquilibriumConstants(double *rop)
Multiply rate with inverse equilibrium constant.
void getDerivativeSettings(AnyMap &settings) const override
Retrieve derivative settings.
void process_ddT(const vector< double > &in, double *drop)
Process temperature derivative.
void getDeltaEnthalpy(double *deltaH) override
Return the vector of values for the reactions change in enthalpy.
vector< unique_ptr< MultiRateBase > > m_bulk_rates
Vector of rate handlers.
void getRevRatesOfProgress_ddT(double *drop) override
Calculate derivatives for reverse rates-of-progress with respect to temperature at constant pressure,...
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.
void getFwdRatesOfProgress_ddC(double *drop) override
Calculate derivatives for forward rates-of-progress with respect to molar concentration at constant t...
void process_ddC(StoichManagerN &stoich, const vector< double > &in, double *drop, bool mass_action=true)
Process concentration (molar density) derivative.
void process_ddP(const vector< double > &in, double *drop)
Process pressure derivative.
void getFwdRateConstants_ddP(double *dkfwd) override
Calculate derivatives for forward rate constants with respect to pressure at constant temperature,...
void getRevRatesOfProgress_ddP(double *drop) override
Calculate derivatives for reverse rates-of-progress with respect to pressure at constant temperature,...
void getRevRatesOfProgress_ddC(double *drop) override
Calculate derivatives for reverse rates-of-progress with respect to molar concentration at constant t...
vector< double > m_rbuf0
Buffers for partial rop results with length nReactions()
vector< size_t > m_irrev
Indices of irreversible reactions.
ThirdBodyCalc m_multi_concm
used with MultiRate evaluator
void applyEquilibriumConstants_ddT(double *drkcn)
Multiply rate with scaled temperature derivatives of the inverse equilibrium constant.
void resizeReactions() override
Finalize Kinetics object and associated objects.
Eigen::SparseMatrix< double > netRatesOfProgress_ddCi() override
Calculate derivatives for net rates-of-progress with respect to species concentration at constant tem...
void getThirdBodyConcentrations(double *concm) override
Return a vector of values of effective concentrations of third-body collision partners of any reactio...
void getDeltaEntropy(double *deltaS) override
Return the vector of values for the reactions change in entropy.
vector< double > m_concm
Third body concentrations.
vector< double > m_kf0
Forward rate constants without perturbation.
void setDerivativeSettings(const AnyMap &settings) override
Set/modify derivative settings.
void getFwdRateConstants_ddC(double *dkfwd) override
Calculate derivatives for forward rate constants with respect to molar concentration at constant temp...
Base class for exceptions thrown by Cantera classes.
virtual void resizeReactions()
Finalize Kinetics object and associated objects.
Definition Kinetics.cpp:35
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
ValueCache m_cache
Cache for saved calculations within each Kinetics object.
Definition Kinetics.h:1403
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
bool m_ready
Boolean indicating whether Kinetics object is fully configured.
Definition Kinetics.h:1444
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< 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_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
virtual bool addReaction(shared_ptr< Reaction > r, bool resize=true)
Add a single reaction to the mechanism.
Definition Kinetics.cpp:565
bool m_skipUndeclaredThirdBodies
See skipUndeclaredThirdBodies()
Definition Kinetics.h:1514
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
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
bool m_hasUndeclaredThirdBodies
Flag indicating whether reactions include undeclared third bodies.
Definition Kinetics.h:1517
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
vector< double > m_delta_gibbs0
Delta G^0 for all reactions.
Definition Kinetics.h:1490
An error indicating that an unimplemented function has been called.
virtual void getConcentrations(double *const c) const
Get the species concentrations (kmol/m^3).
Definition Phase.cpp:482
virtual double molarDensity() const
Molar density (kmol/m^3).
Definition Phase.cpp:576
void restoreState(const vector< double > &state)
Restore a state saved on a previous call to saveState.
Definition Phase.cpp:260
double temperature() const
Temperature (K).
Definition Phase.h:562
void saveState(vector< double > &state) const
Save the current internal state of the phase.
Definition Phase.cpp:236
virtual double density() const
Density (kg/m^3).
Definition Phase.h:587
int stateMFNumber() const
Return the State Mole Fraction Number.
Definition Phase.h:761
virtual double pressure() const
Return the thermodynamic pressure (Pa).
Definition Phase.h:580
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.
void scale(const double *in, double *out, double factor) const
Scale input by reaction order and factor.
virtual void getPartialMolarEnthalpies(double *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
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 void setState_TP(double t, double p)
Set the temperature (K) and pressure (Pa)
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.
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 getChemPotentials(double *mu) const
Get the species chemical potentials. Units: J/kmol.
virtual void getPartialMolarEntropies(double *sbar) const
Returns an array of partial molar entropies of the species in the solution.
void scaleM(const double *in, double *out, const double *concm, double factor) const
Scale entries involving third-body collider in rate expression by third-body concentration and factor...
bool empty() const
Return boolean indicating whether ThirdBodyCalc is empty.
void install(size_t rxnNumber, const map< size_t, double > &efficiencies, double default_efficiency, bool mass_action)
Install reaction that uses third-body effects in ThirdBodyCalc manager.
void multiply(double *output, const double *concm)
Multiply output with effective third-body concentration.
void update(const vector< double > &conc, double ctot, double *concm) const
Update third-body concentrations in full vector.
Eigen::SparseMatrix< double > derivatives(const double *product)
Calculate derivatives with respect to species concentrations.
void scale(const double *in, double *out, double factor) const
Scale entries involving third-body collider in law of mass action by factor.
void resizeCoeffs(size_t nSpc, size_t nRxn)
Resize the sparse coefficient matrix.
CachedScalar getScalar(int id)
Get a reference to a CachedValue object representing a scalar (double) with the given id.
Definition ValueCache.h:161
int getId()
Get a unique id for a cached value.
#define AssertFinite(expr, procedure,...)
Throw an exception if the specified exception is not a finite number.
bool legacy_rate_constants_used()
Returns true if legacy rate constant definition is used.
Definition global.cpp:107
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
const double BigNumber
largest number to compare to inf.
Definition ct_defs.h:160
A cached property value and the state at which it was evaluated.
Definition ValueCache.h:33
double state2
Value of the second state variable for the state at which value was evaluated, for example density or...
Definition ValueCache.h:106
bool validate(double state1New)
Check whether the currently cached value is valid based on a single state variable.
Definition ValueCache.h:39
double state1
Value of the first state variable for the state at which value was evaluated, for example temperature...
Definition ValueCache.h:102