Cantera  3.0.0
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
15BulkKinetics::BulkKinetics(ThermoPhase* thermo) : BulkKinetics()
16{
17 warn_deprecated("BulkKinetics::BulkKinetics(ThermoPhase*)",
18 "To be removed after Cantera 3.0. Use default constructor instead.");
20}
21
23 return std::find(m_revindex.begin(), m_revindex.end(), i) < m_revindex.end();
24}
25
26bool BulkKinetics::addReaction(shared_ptr<Reaction> r, bool resize)
27{
28 bool added = Kinetics::addReaction(r, resize);
29 if (!added) {
30 // undeclared species, etc.
31 return false;
32 }
33 double dn = 0.0;
34 for (const auto& [name, stoich] : r->products) {
35 dn += stoich;
36 }
37 for (const auto& [name, stoich] : r->reactants) {
38 dn -= stoich;
39 }
40
41 m_dn.push_back(dn);
42
43 if (r->reversible) {
44 m_revindex.push_back(nReactions()-1);
45 } else {
46 m_irrev.push_back(nReactions()-1);
47 }
48
49 shared_ptr<ReactionRate> rate = r->rate();
50 string rtype = rate->subType();
51 if (rtype == "") {
52 rtype = rate->type();
53 }
54
55 // If necessary, add new MultiRate evaluator
56 if (m_bulk_types.find(rtype) == m_bulk_types.end()) {
57 m_bulk_types[rtype] = m_bulk_rates.size();
58 m_bulk_rates.push_back(rate->newMultiRate());
59 m_bulk_rates.back()->resize(m_kk, nReactions(), nPhases());
60 }
61
62 // Set index of rate to number of reaction within kinetics
63 rate->setRateIndex(nReactions() - 1);
64 rate->setContext(*r, *this);
65
66 // Add reaction rate to evaluator
67 size_t index = m_bulk_types[rtype];
68 m_bulk_rates[index]->add(nReactions() - 1, *rate);
69
70 // Add reaction to third-body evaluator
71 if (r->thirdBody() != nullptr) {
72 addThirdBody(r);
73 }
74
75 m_concm.push_back(NAN);
76 m_ready = resize;
77 return true;
78}
79
80void BulkKinetics::addThirdBody(shared_ptr<Reaction> r)
81{
82 map<size_t, double> efficiencies;
83 for (const auto& [name, efficiency] : r->thirdBody()->efficiencies) {
84 size_t k = kineticsSpeciesIndex(name);
85 if (k != npos) {
86 efficiencies[k] = efficiency;
87 } else if (!m_skipUndeclaredThirdBodies) {
88 throw CanteraError("BulkKinetics::addThirdBody", "Found third-body"
89 " efficiency for undefined species '{}' while adding reaction '{}'",
90 name, r->equation());
91 } else {
93 }
94 }
95 m_multi_concm.install(nReactions() - 1, efficiencies,
96 r->thirdBody()->default_efficiency,
97 r->thirdBody()->mass_action);
98}
99
100void BulkKinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
101{
102 // operations common to all reaction types
104
105 shared_ptr<ReactionRate> rate = rNew->rate();
106 string rtype = rate->subType();
107 if (rtype == "") {
108 rtype = rate->type();
109 }
110
111 // Ensure that MultiRate evaluator is available
112 if (m_bulk_types.find(rtype) == m_bulk_types.end()) {
113 throw CanteraError("BulkKinetics::modifyReaction",
114 "Evaluator not available for type '{}'.", rtype);
115 }
116
117 // Replace reaction rate to evaluator
118 size_t index = m_bulk_types[rtype];
119 rate->setRateIndex(i);
120 rate->setContext(*rNew, *this);
121 m_bulk_rates[index]->replace(i, *rate);
122 invalidateCache();
123}
124
126{
128 m_act_conc.resize(m_kk);
129 m_phys_conc.resize(m_kk);
130 m_grt.resize(m_kk);
131 for (auto& rates : m_bulk_rates) {
132 rates->resize(m_kk, nReactions(), nPhases());
133 }
134}
135
137{
139 m_rbuf0.resize(nReactions());
140 m_rbuf1.resize(nReactions());
141 m_rbuf2.resize(nReactions());
142 m_kf0.resize(nReactions());
143 m_sbuf0.resize(nTotalSpecies());
144 m_state.resize(thermo().stateSize());
146 for (auto& rates : m_bulk_rates) {
147 rates->resize(nTotalSpecies(), nReactions(), nPhases());
148 // @todo ensure that ReactionData are updated; calling rates->update
149 // blocks correct behavior in update_rates_T
150 // and running updateROP() is premature
151 }
152}
153
154void BulkKinetics::setMultiplier(size_t i, double f)
155{
157 m_ROP_ok = false;
158}
159
160void BulkKinetics::invalidateCache()
161{
162 Kinetics::invalidateCache();
163 m_ROP_ok = false;
164}
165
167{
168 updateROP();
169 copy(m_rfn.begin(), m_rfn.end(), kfwd);
171 processThirdBodies(kfwd);
172 }
173}
174
176{
177 updateROP();
178
179 vector<double>& delta_gibbs0 = m_rbuf0;
180 fill(delta_gibbs0.begin(), delta_gibbs0.end(), 0.0);
181
182 // compute Delta G^0 for all reactions
183 getReactionDelta(m_grt.data(), delta_gibbs0.data());
184
185 double rrt = 1.0 / thermo().RT();
186 double logStandConc = log(thermo().standardConcentration());
187 for (size_t i = 0; i < nReactions(); i++) {
188 kc[i] = exp(-delta_gibbs0[i] * rrt + m_dn[i] * logStandConc);
189 }
190}
191
192void BulkKinetics::getRevRateConstants(double* krev, bool doIrreversible)
193{
194 // go get the forward rate constants. -> note, we don't really care about
195 // speed or redundancy in these informational routines.
197
198 if (doIrreversible) {
200 for (size_t i = 0; i < nReactions(); i++) {
201 krev[i] /= m_rbuf0[i];
202 }
203 } else {
204 // m_rkcn[] is zero for irreversible reactions
205 for (size_t i = 0; i < nReactions(); i++) {
206 krev[i] *= m_rkcn[i];
207 }
208 }
209}
210
211void BulkKinetics::getDeltaGibbs(double* deltaG)
212{
213 // Get the chemical potentials for each species
214 thermo().getChemPotentials(m_sbuf0.data());
215 // Use the stoichiometric manager to find deltaG for each reaction.
216 getReactionDelta(m_sbuf0.data(), deltaG);
217}
218
220{
221 // Get the partial molar enthalpy for each species
222 thermo().getPartialMolarEnthalpies(m_sbuf0.data());
223 // Use the stoichiometric manager to find deltaH for each reaction.
224 getReactionDelta(m_sbuf0.data(), deltaH);
225}
226
228{
229 // Get the partial molar entropy for each species
230 thermo().getPartialMolarEntropies(m_sbuf0.data());
231 // Use the stoichiometric manager to find deltaS for each reaction.
232 getReactionDelta(m_sbuf0.data(), deltaS);
233}
234
236{
237 // Get the standard state chemical potentials of the species. This is the
238 // array of chemical potentials at unit activity. We define these here as
239 // the chemical potentials of the pure species at the temperature and
240 // pressure of the solution.
241 thermo().getStandardChemPotentials(m_sbuf0.data());
242 // Use the stoichiometric manager to find deltaG for each reaction.
243 getReactionDelta(m_sbuf0.data(), deltaG);
244}
245
247{
248 // Get the standard state enthalpies of the species.
249 thermo().getEnthalpy_RT(m_sbuf0.data());
250 for (size_t k = 0; k < m_kk; k++) {
251 m_sbuf0[k] *= thermo().RT();
252 }
253 // Use the stoichiometric manager to find deltaH for each reaction.
254 getReactionDelta(m_sbuf0.data(), deltaH);
255}
256
258{
259 // Get the standard state entropy of the species. We define these here as
260 // the entropies of the pure species at the temperature and pressure of the
261 // solution.
262 thermo().getEntropy_R(m_sbuf0.data());
263 for (size_t k = 0; k < m_kk; k++) {
264 m_sbuf0[k] *= GasConstant;
265 }
266 // Use the stoichiometric manager to find deltaS for each reaction.
267 getReactionDelta(m_sbuf0.data(), deltaS);
268}
269
271{
272 settings["skip-third-bodies"] = m_jac_skip_third_bodies;
273 settings["skip-falloff"] = m_jac_skip_falloff;
274 settings["rtol-delta"] = m_jac_rtol_delta;
275}
276
278{
279 bool force = settings.empty();
280 if (force || settings.hasKey("skip-third-bodies")) {
281 m_jac_skip_third_bodies = settings.getBool("skip-third-bodies", false);
282 }
283 if (force || settings.hasKey("skip-falloff")) {
284 m_jac_skip_falloff = settings.getBool("skip-falloff", false);
285 }
286 if (force || settings.hasKey("rtol-delta")) {
287 m_jac_rtol_delta = settings.getDouble("rtol-delta", 1e-8);
288 }
289}
290
292{
293 assertDerivativesValid("BulkKinetics::getFwdRateConstants_ddT");
294 updateROP();
295 process_ddT(m_rfn, dkfwd);
296}
297
299{
300 assertDerivativesValid("BulkKinetics::getFwdRatesOfProgress_ddT");
301 updateROP();
302 process_ddT(m_ropf, drop);
303}
304
306{
307 assertDerivativesValid("BulkKinetics::getRevRatesOfProgress_ddT");
308 updateROP();
309 process_ddT(m_ropr, drop);
310 Eigen::Map<Eigen::VectorXd> dRevRop(drop, nReactions());
311
312 // reverse rop times scaled inverse equilibrium constant derivatives
313 Eigen::Map<Eigen::VectorXd> dRevRop2(m_rbuf2.data(), nReactions());
314 copy(m_ropr.begin(), m_ropr.end(), m_rbuf2.begin());
315 applyEquilibriumConstants_ddT(dRevRop2.data());
316 dRevRop += dRevRop2;
317}
318
320{
321 assertDerivativesValid("BulkKinetics::getNetRatesOfProgress_ddT");
322 updateROP();
323 process_ddT(m_ropnet, drop);
324 Eigen::Map<Eigen::VectorXd> dNetRop(drop, nReactions());
325
326 // reverse rop times scaled inverse equilibrium constant derivatives
327 Eigen::Map<Eigen::VectorXd> dNetRop2(m_rbuf2.data(), nReactions());
328 copy(m_ropr.begin(), m_ropr.end(), m_rbuf2.begin());
329 applyEquilibriumConstants_ddT(dNetRop2.data());
330 dNetRop -= dNetRop2;
331}
332
334{
335 assertDerivativesValid("BulkKinetics::getFwdRateConstants_ddP");
336 updateROP();
337 process_ddP(m_rfn, dkfwd);
338}
339
341{
342 assertDerivativesValid("BulkKinetics::getFwdRatesOfProgress_ddP");
343 updateROP();
344 process_ddP(m_ropf, drop);
345}
346
348{
349 assertDerivativesValid("BulkKinetics::getRevRatesOfProgress_ddP");
350 updateROP();
351 process_ddP(m_ropr, drop);
352}
353
355{
356 assertDerivativesValid("BulkKinetics::getNetRatesOfProgress_ddP");
357 updateROP();
358 process_ddP(m_ropnet, drop);
359}
360
362{
363 assertDerivativesValid("BulkKinetics::getFwdRateConstants_ddC");
364 updateROP();
365 process_ddC(m_reactantStoich, m_rfn, dkfwd, false);
366}
367
369{
370 assertDerivativesValid("BulkKinetics::getFwdRatesOfProgress_ddC");
371 updateROP();
373}
374
376{
377 assertDerivativesValid("BulkKinetics::getRevRatesOfProgress_ddC");
378 updateROP();
380}
381
383{
384 assertDerivativesValid("BulkKinetics::getNetRatesOfProgress_ddC");
385 updateROP();
387 Eigen::Map<Eigen::VectorXd> dNetRop(drop, nReactions());
388
389 process_ddC(m_revProductStoich, m_ropr, m_rbuf2.data());
390 Eigen::Map<Eigen::VectorXd> dNetRop2(m_rbuf2.data(), nReactions());
391 dNetRop -= dNetRop2;
392}
393
394Eigen::SparseMatrix<double> BulkKinetics::fwdRatesOfProgress_ddX()
395{
396 assertDerivativesValid("BulkKinetics::fwdRatesOfProgress_ddX");
397
398 // forward reaction rate coefficients
399 vector<double>& rop_rates = m_rbuf0;
400 getFwdRateConstants(rop_rates.data());
402}
403
404Eigen::SparseMatrix<double> BulkKinetics::revRatesOfProgress_ddX()
405{
406 assertDerivativesValid("BulkKinetics::revRatesOfProgress_ddX");
407
408 // reverse reaction rate coefficients
409 vector<double>& rop_rates = m_rbuf0;
410 getFwdRateConstants(rop_rates.data());
411 applyEquilibriumConstants(rop_rates.data());
413}
414
415Eigen::SparseMatrix<double> BulkKinetics::netRatesOfProgress_ddX()
416{
417 assertDerivativesValid("BulkKinetics::netRatesOfProgress_ddX");
418
419 // forward reaction rate coefficients
420 vector<double>& rop_rates = m_rbuf0;
421 getFwdRateConstants(rop_rates.data());
423
424 // reverse reaction rate coefficients
425 applyEquilibriumConstants(rop_rates.data());
427}
428
429Eigen::SparseMatrix<double> BulkKinetics::fwdRatesOfProgress_ddCi()
430{
431 assertDerivativesValid("BulkKinetics::fwdRatesOfProgress_ddCi");
432
433 // forward reaction rate coefficients
434 vector<double>& rop_rates = m_rbuf0;
435 getFwdRateConstants(rop_rates.data());
436 return calculateCompositionDerivatives(m_reactantStoich, rop_rates, false);
437}
438
439Eigen::SparseMatrix<double> BulkKinetics::revRatesOfProgress_ddCi()
440{
441 assertDerivativesValid("BulkKinetics::revRatesOfProgress_ddCi");
442
443 // reverse reaction rate coefficients
444 vector<double>& rop_rates = m_rbuf0;
445 getFwdRateConstants(rop_rates.data());
446 applyEquilibriumConstants(rop_rates.data());
447 return calculateCompositionDerivatives(m_revProductStoich, rop_rates, false);
448}
449
450Eigen::SparseMatrix<double> BulkKinetics::netRatesOfProgress_ddCi()
451{
452 assertDerivativesValid("BulkKinetics::netRatesOfProgress_ddCi");
453
454 // forward reaction rate coefficients
455 vector<double>& rop_rates = m_rbuf0;
456 getFwdRateConstants(rop_rates.data());
457 auto jac = calculateCompositionDerivatives(m_reactantStoich, rop_rates, false);
458
459 // reverse reaction rate coefficients
460 applyEquilibriumConstants(rop_rates.data());
461 return jac - calculateCompositionDerivatives(m_revProductStoich, rop_rates, false);
462}
463
464void BulkKinetics::updateROP()
465{
466 static const int cacheId = m_cache.getId();
467 CachedScalar last = m_cache.getScalar(cacheId);
468 double T = thermo().temperature();
469 double rho = thermo().density();
470 int statenum = thermo().stateMFNumber();
471
472 if (last.state1 != T || last.state2 != rho) {
473 // Update properties that are independent of the composition
475 fill(m_delta_gibbs0.begin(), m_delta_gibbs0.end(), 0.0);
476 double logStandConc = log(thermo().standardConcentration());
477
478 // compute Delta G^0 for all reversible reactions
480
481 double rrt = 1.0 / thermo().RT();
482 for (size_t i = 0; i < m_revindex.size(); i++) {
483 size_t irxn = m_revindex[i];
484 m_rkcn[irxn] = std::min(
485 exp(m_delta_gibbs0[irxn] * rrt - m_dn[irxn] * logStandConc), BigNumber);
486 }
487
488 for (size_t i = 0; i != m_irrev.size(); ++i) {
489 m_rkcn[ m_irrev[i] ] = 0.0;
490 }
491 }
492
493 if (!last.validate(T, rho, statenum)) {
494 // Update terms dependent on species concentrations and temperature
497 double ctot = thermo().molarDensity();
498
499 // Third-body objects interacting with MultiRate evaluator
501
502 // loop over MultiRate evaluators for each reaction type
503 for (auto& rates : m_bulk_rates) {
504 bool changed = rates->update(thermo(), *this);
505 if (changed) {
506 rates->getRateConstants(m_kf0.data());
507 }
508 }
509 m_ROP_ok = false;
510 }
511
512 if (m_ROP_ok) {
513 // rates of progress are up-to-date only if both the thermodynamic state
514 // and m_perturb are unchanged
515 return;
516 }
517
518 // Scale the forward rate coefficient by the perturbation factor
519 for (size_t i = 0; i < nReactions(); ++i) {
520 m_rfn[i] = m_kf0[i] * m_perturb[i];
521 }
522
523 copy(m_rfn.begin(), m_rfn.end(), m_ropf.data());
525 copy(m_ropf.begin(), m_ropf.end(), m_ropr.begin());
526
527 // multiply ropf by concentration products
528 m_reactantStoich.multiply(m_act_conc.data(), m_ropf.data());
529
530 // for reversible reactions, multiply ropr by concentration products
532 m_revProductStoich.multiply(m_act_conc.data(), m_ropr.data());
533 for (size_t j = 0; j != nReactions(); ++j) {
534 m_ropnet[j] = m_ropf[j] - m_ropr[j];
535 }
536
537 for (size_t i = 0; i < m_rfn.size(); i++) {
538 AssertFinite(m_rfn[i], "BulkKinetics::updateROP",
539 "m_rfn[{}] is not finite.", i);
540 AssertFinite(m_ropf[i], "BulkKinetics::updateROP",
541 "m_ropf[{}] is not finite.", i);
542 AssertFinite(m_ropr[i], "BulkKinetics::updateROP",
543 "m_ropr[{}] is not finite.", i);
544 }
545 m_ROP_ok = true;
546}
547
549{
550 updateROP();
551 std::copy(m_concm.begin(), m_concm.end(), concm);
552}
553
555{
556 // reactions involving third body
557 if (!m_concm.empty()) {
558 m_multi_concm.multiply(rop, m_concm.data());
559 }
560}
561
563{
564 // For reverse rates computed from thermochemistry, multiply the forward
565 // rate coefficients by the reciprocals of the equilibrium constants
566 for (size_t i = 0; i < nReactions(); ++i) {
567 rop[i] *= m_rkcn[i];
568 }
569}
570
572{
573 double T = thermo().temperature();
574 double P = thermo().pressure();
575 double rrt = 1. / thermo().RT();
576
577 vector<double>& grt = m_sbuf0;
578 vector<double>& delta_gibbs0 = m_rbuf1;
579 fill(delta_gibbs0.begin(), delta_gibbs0.end(), 0.0);
580
581 // compute perturbed Delta G^0 for all reversible reactions
582 thermo().saveState(m_state);
583 thermo().setState_TP(T * (1. + m_jac_rtol_delta), P);
584 thermo().getStandardChemPotentials(grt.data());
585 getRevReactionDelta(grt.data(), delta_gibbs0.data());
586
587 // apply scaling for derivative of inverse equilibrium constant
588 double Tinv = 1. / T;
589 double rrt_dTinv = rrt * Tinv / m_jac_rtol_delta;
590 double rrtt = rrt * Tinv;
591 for (size_t i = 0; i < m_revindex.size(); i++) {
592 size_t irxn = m_revindex[i];
593 double factor = delta_gibbs0[irxn] - m_delta_gibbs0[irxn];
594 factor *= rrt_dTinv;
595 factor += m_dn[irxn] * Tinv - m_delta_gibbs0[irxn] * rrtt;
596 drkcn[irxn] *= factor;
597 }
598
599 for (size_t i = 0; i < m_irrev.size(); ++i) {
600 drkcn[m_irrev[i]] = 0.0;
601 }
602
603 thermo().restoreState(m_state);
604}
605
606void BulkKinetics::process_ddT(const vector<double>& in, double* drop)
607{
608 // apply temperature derivative
609 copy(in.begin(), in.end(), drop);
610 for (auto& rates : m_bulk_rates) {
611 rates->processRateConstants_ddT(drop, m_rfn.data(), m_jac_rtol_delta);
612 }
613}
614
615void BulkKinetics::process_ddP(const vector<double>& in, double* drop)
616{
617 // apply pressure derivative
618 copy(in.begin(), in.end(), drop);
619 for (auto& rates : m_bulk_rates) {
620 rates->processRateConstants_ddP(drop, m_rfn.data(), m_jac_rtol_delta);
621 }
622}
623
624void BulkKinetics::process_ddC(StoichManagerN& stoich, const vector<double>& in,
625 double* drop, bool mass_action)
626{
627 Eigen::Map<Eigen::VectorXd> out(drop, nReactions());
628 out.setZero();
629 double ctot_inv = 1. / thermo().molarDensity();
630
631 // derivatives due to concentrations in law of mass action
632 if (mass_action) {
633 stoich.scale(in.data(), out.data(), ctot_inv);
634 }
636 return;
637 }
638
639 // derivatives due to third-body colliders in law of mass action
640 Eigen::Map<Eigen::VectorXd> outM(m_rbuf1.data(), nReactions());
641 if (mass_action) {
642 outM.fill(0.);
643 m_multi_concm.scale(in.data(), outM.data(), ctot_inv);
644 out += outM;
645 }
646
647 // derivatives due to reaction rates depending on third-body colliders
648 if (!m_jac_skip_falloff) {
649 m_multi_concm.scaleM(in.data(), outM.data(), m_concm.data(), ctot_inv);
650 for (auto& rates : m_bulk_rates) {
651 // processing step assigns zeros to entries not dependent on M
652 rates->processRateConstants_ddM(
653 outM.data(), m_rfn.data(), m_jac_rtol_delta);
654 }
655 out += outM;
656 }
657}
658
660 StoichManagerN& stoich, const vector<double>& in, bool ddX)
661{
662 Eigen::SparseMatrix<double> out;
663 vector<double>& scaled = m_rbuf1;
664 vector<double>& outV = m_rbuf2;
665
666 // convert from concentration to mole fraction output
667 copy(in.begin(), in.end(), scaled.begin());
668 if (ddX) {
669 double ctot = thermo().molarDensity();
670 for (size_t i = 0; i < nReactions(); ++i) {
671 scaled[i] *= ctot;
672 }
673 }
674
675 // derivatives handled by StoichManagerN
676 copy(scaled.begin(), scaled.end(), outV.begin());
677 processThirdBodies(outV.data());
678 out = stoich.derivatives(m_act_conc.data(), outV.data());
680 return out;
681 }
682
683 // derivatives due to law of mass action
684 copy(scaled.begin(), scaled.end(), outV.begin());
685 stoich.multiply(m_act_conc.data(), outV.data());
686
687 // derivatives due to reaction rates depending on third-body colliders
688 if (!m_jac_skip_falloff) {
689 for (auto& rates : m_bulk_rates) {
690 // processing step does not modify entries not dependent on M
691 rates->processRateConstants_ddM(
692 outV.data(), m_rfn.data(), m_jac_rtol_delta, false);
693 }
694 }
695
696 // derivatives handled by ThirdBodyCalc
697 out += m_multi_concm.derivatives(outV.data());
698 return out;
699}
700
702{
703 if (!thermo().isIdeal()) {
704 throw NotImplementedError(name,
705 "Not supported for non-ideal ThermoPhase models.");
706 }
707}
708
709}
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
Specialization of Kinetics for chemistry in a single bulk phase.
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:254
ValueCache m_cache
Cache for saved calculations within each Kinetics object.
Definition Kinetics.h:1497
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
bool m_ready
Boolean indicating whether Kinetics object is fully configured.
Definition Kinetics.h:1538
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< double > m_rkcn
Reciprocal of the equilibrium constant in concentration units.
Definition Kinetics.h:1605
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_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
bool m_skipUndeclaredThirdBodies
See skipUndeclaredThirdBodies()
Definition Kinetics.h:1626
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
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
bool m_hasUndeclaredThirdBodies
Flag indicating whether reactions include undeclared third bodies.
Definition Kinetics.h:1629
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
vector< double > m_delta_gibbs0
Delta G^0 for all reactions.
Definition Kinetics.h:1602
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:595
virtual double molarDensity() const
Molar density (kmol/m^3).
Definition Phase.cpp:689
void restoreState(const vector< double > &state)
Restore a state saved on a previous call to saveState.
Definition Phase.cpp:275
double temperature() const
Temperature (K).
Definition Phase.h:662
void saveState(vector< double > &state) const
Save the current internal state of the phase.
Definition Phase.cpp:251
virtual double density() const
Density (kg/m^3).
Definition Phase.h:687
int stateMFNumber() const
Return the State Mole Fraction Number.
Definition Phase.h:866
virtual double pressure() const
Return the thermodynamic pressure (Pa).
Definition Phase.h:680
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.
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.
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: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
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