Cantera  3.1.0a1
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 
8 namespace Cantera
9 {
10 
11 BulkKinetics::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 
19 bool 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 
73 void 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 
93 void BulkKinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
94 {
95  // operations common to all reaction types
96  Kinetics::modifyReaction(i, rNew);
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 
147 void BulkKinetics::setMultiplier(size_t i, double f)
148 {
150  m_ROP_ok = false;
151 }
152 
153 void 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 
185 void 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.
189  getFwdRateConstants(krev);
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 
204 void 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 
212 void BulkKinetics::getDeltaEnthalpy(double* deltaH)
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 
220 void BulkKinetics::getDeltaEntropy(double* deltaS)
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 
228 void BulkKinetics::getDeltaSSGibbs(double* deltaG)
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();
372  return process_ddC(m_revProductStoich, m_ropr, drop);
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 
387 Eigen::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 
397 Eigen::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 
408 Eigen::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());
415  auto jac = calculateCompositionDerivatives(m_reactantStoich, rop_rates);
416 
417  // reverse reaction rate coefficients
418  applyEquilibriumConstants(rop_rates.data());
419  return jac - calculateCompositionDerivatives(m_revProductStoich, rop_rates);
420 }
421 
422 Eigen::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 
432 Eigen::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 
443 Eigen::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 
457 void 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
472  getRevReactionDelta(m_grt.data(), m_delta_gibbs0.data());
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
493  m_multi_concm.update(m_phys_conc, ctot, m_concm.data());
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());
517  processThirdBodies(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 
599 void 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 
608 void 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 
617 void 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 
694 void BulkKinetics::assertDerivativesValid(const string& name)
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.
Definition: BulkKinetics.h:153
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.
Definition: BulkKinetics.h:172
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.
Definition: BulkKinetics.h:161
bool m_jac_skip_third_bodies
Derivative settings.
Definition: BulkKinetics.h:175
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.
Definition: BulkKinetics.h:188
vector< double > m_act_conc
Activity concentrations, as calculated by ThermoPhase::getActivityConcentrations.
Definition: BulkKinetics.h:169
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.
Definition: BulkKinetics.h:155
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.
Definition: BulkKinetics.h:152
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()
Definition: BulkKinetics.h:182
vector< size_t > m_irrev
Indices of irreversible reactions.
Definition: BulkKinetics.h:156
ThirdBodyCalc m_multi_concm
used with MultiRate evaluator
Definition: BulkKinetics.h:163
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.
Definition: BulkKinetics.h:166
vector< double > m_kf0
Forward rate constants without perturbation.
Definition: BulkKinetics.h:185
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.
Definition: ctexceptions.h:66
virtual void resizeReactions()
Finalize Kinetics object and associated objects.
Definition: Kinetics.cpp:35
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
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
An error indicating that an unimplemented function has been called.
Definition: ctexceptions.h:195
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.
Definition: ThermoPhase.h:801
virtual void getEntropy_R(double *sr) const
Get the array of nondimensional Entropy functions for the standard state species at the current T and...
Definition: ThermoPhase.h:880
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.
Definition: ThermoPhase.h:1062
virtual void getActivityConcentrations(double *c) const
This method returns an array of generalized concentrations.
Definition: ThermoPhase.h:697
virtual void getStandardChemPotentials(double *mu) const
Get the array of chemical potentials at unit activity for the species at their standard states at the...
Definition: ThermoPhase.h:860
virtual void getEnthalpy_RT(double *hrt) const
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
Definition: ThermoPhase.h:870
virtual void getChemPotentials(double *mu) const
Get the species chemical potentials. Units: J/kmol.
Definition: ThermoPhase.h:775
virtual void getPartialMolarEntropies(double *sbar) const
Returns an array of partial molar entropies of the species in the solution.
Definition: ThermoPhase.h:811
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.
Definition: ThirdBodyCalc.h:23
Eigen::SparseMatrix< double > derivatives(const double *product)
Calculate derivatives with respect to species concentrations.
Definition: ThirdBodyCalc.h:97
void multiply(double *output, const double *concm)
Multiply output with effective third-body concentration.
Definition: ThirdBodyCalc.h:86
void update(const vector< double > &conc, double ctot, double *concm) const
Update third-body concentrations in full vector.
Definition: ThirdBodyCalc.h:75
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.
Definition: ThirdBodyCalc.h:47
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.
Definition: ValueCache.cpp:21
#define AssertFinite(expr, procedure,...)
Throw an exception if the specified exception is not a finite number.
Definition: ctexceptions.h:286
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