19 #include <unordered_set>
20 #include <boost/algorithm/string.hpp>
27 void Kinetics::checkReactionIndex(
size_t i)
const
29 if (i >= nReactions()) {
30 throw IndexError(
"Kinetics::checkReactionIndex",
"reactions", i,
35 void Kinetics::resizeReactions()
37 size_t nRxn = nReactions();
40 m_reactantStoich.resizeCoeffs(m_kk, nRxn);
41 m_productStoich.resizeCoeffs(m_kk, nRxn);
42 m_revProductStoich.resizeCoeffs(m_kk, nRxn);
47 m_stoichMatrix = m_productStoich.stoichCoeffs();
49 m_stoichMatrix -= m_reactantStoich.stoichCoeffs();
54 void Kinetics::checkReactionArraySize(
size_t ii)
const
56 if (nReactions() > ii) {
62 void Kinetics::checkPhaseIndex(
size_t m)
const
65 throw IndexError(
"Kinetics::checkPhaseIndex",
"phase", m, nPhases()-1);
69 void Kinetics::checkPhaseArraySize(
size_t mm)
const
72 throw ArraySizeError(
"Kinetics::checkPhaseArraySize", mm, nPhases());
76 size_t Kinetics::reactionPhaseIndex()
const
78 warn_deprecated(
"Kinetics::reactionPhaseIndex",
"The reacting phase is always "
79 "the first phase. To be removed after Cantera 3.1.");
83 shared_ptr<ThermoPhase> Kinetics::reactionPhase()
const
88 void Kinetics::checkSpeciesIndex(
size_t k)
const
91 throw IndexError(
"Kinetics::checkSpeciesIndex",
"species", k, m_kk-1);
95 void Kinetics::checkSpeciesArraySize(
size_t kk)
const
102 pair<size_t, size_t> Kinetics::checkDuplicates(
bool throw_err)
const
105 map<size_t, vector<size_t>> participants;
106 vector<map<int, double>> net_stoich;
107 std::unordered_set<size_t> unmatched_duplicates;
108 for (
size_t i = 0; i < m_reactions.size(); i++) {
109 if (m_reactions[i]->duplicate) {
110 unmatched_duplicates.insert(i);
114 for (
size_t i = 0; i < m_reactions.size(); i++) {
116 unsigned long int key = 0;
118 net_stoich.emplace_back();
119 map<int, double>& net = net_stoich.back();
120 for (
const auto& [name, stoich] : R.
reactants) {
121 int k =
static_cast<int>(kineticsSpeciesIndex(name));
123 net[-1 -k] -= stoich;
125 for (
const auto& [name, stoich] : R.
products) {
126 int k =
static_cast<int>(kineticsSpeciesIndex(name));
132 vector<size_t>& related = participants[key];
133 for (
size_t m : related) {
137 unmatched_duplicates.erase(i);
138 unmatched_duplicates.erase(m);
140 }
else if (R.
type() != other.
type()) {
143 && R.
rate()->type() != other.
rate()->type())
147 double c = checkDuplicateStoich(net_stoich[i], net_stoich[m]);
155 bool thirdBodyOk =
true;
156 for (
size_t k = 0; k < nTotalSpecies(); k++) {
157 string s = kineticsSpeciesName(k);
172 "Undeclared duplicate reactions detected:\n"
173 "Reaction {}: {}\nReaction {}: {}\n",
179 participants[key].push_back(i);
181 if (unmatched_duplicates.size()) {
182 size_t i = *unmatched_duplicates.begin();
185 m_reactions[i]->input,
186 "No duplicate found for declared duplicate reaction number {}"
187 " ({})", i, m_reactions[i]->equation());
195 double Kinetics::checkDuplicateStoich(map<int, double>& r1, map<int, double>& r2)
const
197 std::unordered_set<int> keys;
198 for (
auto& [speciesKey, stoich] : r1) {
199 keys.insert(speciesKey);
201 for (
auto& [speciesKey, stoich] : r2) {
202 keys.insert(speciesKey);
204 int k1 = r1.begin()->first;
207 if (r1[k1] && r2[k1]) {
208 ratio = r2[k1]/r1[k1];
209 bool different =
false;
211 if ((r1[k] && !r2[k]) ||
213 (r1[k] && fabs(r2[k]/r1[k] - ratio) > 1.e-8)) {
224 if (r1[k1] == 0.0 || r2[-k1] == 0.0) {
227 ratio = r2[-k1]/r1[k1];
229 if ((r1[k] && !r2[-k]) ||
230 (!r1[k] && r2[-k]) ||
231 (r1[k] && fabs(r2[-k]/r1[k] - ratio) > 1.e-8)) {
238 string Kinetics::kineticsSpeciesName(
size_t k)
const
240 for (
size_t n = m_start.size()-1; n !=
npos; n--) {
241 if (k >= m_start[n]) {
242 return thermo(n).speciesName(k - m_start[n]);
248 size_t Kinetics::kineticsSpeciesIndex(
const string& nm)
const
250 for (
size_t n = 0; n < m_thermo.size(); n++) {
252 size_t k = thermo(n).speciesIndex(nm);
254 return k + m_start[n];
262 for (
size_t n = 0; n < m_thermo.size(); n++) {
268 throw CanteraError(
"Kinetics::speciesPhase",
"unknown species '{}'", nm);
271 const ThermoPhase& Kinetics::speciesPhase(
const string& nm)
const
273 for (
const auto& thermo : m_thermo) {
274 if (thermo->speciesIndex(nm) !=
npos) {
278 throw CanteraError(
"Kinetics::speciesPhase",
"unknown species '{}'", nm);
281 size_t Kinetics::speciesPhaseIndex(
size_t k)
const
283 for (
size_t n = m_start.size()-1; n !=
npos; n--) {
284 if (k >= m_start[n]) {
289 "illegal species index: {}", k);
292 double Kinetics::reactantStoichCoeff(
size_t kSpec,
size_t irxn)
const
294 return m_reactantStoich.stoichCoeffs().coeff(kSpec, irxn);
297 double Kinetics::productStoichCoeff(
size_t kSpec,
size_t irxn)
const
299 return m_productStoich.stoichCoeffs().coeff(kSpec, irxn);
302 void Kinetics::getFwdRatesOfProgress(
double* fwdROP)
305 std::copy(m_ropf.begin(), m_ropf.end(), fwdROP);
308 void Kinetics::getRevRatesOfProgress(
double* revROP)
311 std::copy(m_ropr.begin(), m_ropr.end(), revROP);
314 void Kinetics::getNetRatesOfProgress(
double* netROP)
317 std::copy(m_ropnet.begin(), m_ropnet.end(), netROP);
320 void Kinetics::getReactionDelta(
const double* prop,
double* deltaProp)
const
322 fill(deltaProp, deltaProp + nReactions(), 0.0);
324 m_productStoich.incrementReactions(prop, deltaProp);
326 m_reactantStoich.decrementReactions(prop, deltaProp);
329 void Kinetics::getRevReactionDelta(
const double* prop,
double* deltaProp)
const
331 fill(deltaProp, deltaProp + nReactions(), 0.0);
333 m_revProductStoich.incrementReactions(prop, deltaProp);
335 m_reactantStoich.decrementReactions(prop, deltaProp);
338 void Kinetics::getCreationRates(
double* cdot)
343 fill(cdot, cdot + m_kk, 0.0);
346 m_productStoich.incrementSpecies(m_ropf.data(), cdot);
349 m_reactantStoich.incrementSpecies(m_ropr.data(), cdot);
352 void Kinetics::getDestructionRates(
double* ddot)
356 fill(ddot, ddot + m_kk, 0.0);
358 m_revProductStoich.incrementSpecies(m_ropr.data(), ddot);
360 m_reactantStoich.incrementSpecies(m_ropf.data(), ddot);
363 void Kinetics::getNetProductionRates(
double* net)
367 fill(net, net + m_kk, 0.0);
369 m_productStoich.incrementSpecies(m_ropnet.data(), net);
371 m_reactantStoich.decrementSpecies(m_ropnet.data(), net);
374 void Kinetics::getCreationRates_ddT(
double* dwdot)
376 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
377 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
379 getFwdRatesOfProgress_ddT(buf.data());
380 out = m_productStoich.stoichCoeffs() * buf;
382 getRevRatesOfProgress_ddT(buf.data());
383 out += m_reactantStoich.stoichCoeffs() * buf;
386 void Kinetics::getCreationRates_ddP(
double* dwdot)
388 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
389 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
391 getFwdRatesOfProgress_ddP(buf.data());
392 out = m_productStoich.stoichCoeffs() * buf;
394 getRevRatesOfProgress_ddP(buf.data());
395 out += m_reactantStoich.stoichCoeffs() * buf;
398 void Kinetics::getCreationRates_ddC(
double* dwdot)
400 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
401 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
403 getFwdRatesOfProgress_ddC(buf.data());
404 out = m_productStoich.stoichCoeffs() * buf;
406 getRevRatesOfProgress_ddC(buf.data());
407 out += m_reactantStoich.stoichCoeffs() * buf;
410 Eigen::SparseMatrix<double> Kinetics::creationRates_ddX()
412 Eigen::SparseMatrix<double> jac;
414 jac = m_productStoich.stoichCoeffs() * fwdRatesOfProgress_ddX();
416 jac += m_reactantStoich.stoichCoeffs() * revRatesOfProgress_ddX();
420 Eigen::SparseMatrix<double> Kinetics::creationRates_ddCi()
422 Eigen::SparseMatrix<double> jac;
424 jac = m_productStoich.stoichCoeffs() * fwdRatesOfProgress_ddCi();
426 jac += m_reactantStoich.stoichCoeffs() * revRatesOfProgress_ddCi();
430 void Kinetics::getDestructionRates_ddT(
double* dwdot)
432 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
433 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
435 getRevRatesOfProgress_ddT(buf.data());
436 out = m_revProductStoich.stoichCoeffs() * buf;
438 getFwdRatesOfProgress_ddT(buf.data());
439 out += m_reactantStoich.stoichCoeffs() * buf;
442 void Kinetics::getDestructionRates_ddP(
double* dwdot)
444 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
445 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
447 getRevRatesOfProgress_ddP(buf.data());
448 out = m_revProductStoich.stoichCoeffs() * buf;
450 getFwdRatesOfProgress_ddP(buf.data());
451 out += m_reactantStoich.stoichCoeffs() * buf;
454 void Kinetics::getDestructionRates_ddC(
double* dwdot)
456 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
457 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
459 getRevRatesOfProgress_ddC(buf.data());
460 out = m_revProductStoich.stoichCoeffs() * buf;
462 getFwdRatesOfProgress_ddC(buf.data());
463 out += m_reactantStoich.stoichCoeffs() * buf;
466 Eigen::SparseMatrix<double> Kinetics::destructionRates_ddX()
468 Eigen::SparseMatrix<double> jac;
470 jac = m_revProductStoich.stoichCoeffs() * revRatesOfProgress_ddX();
472 jac += m_reactantStoich.stoichCoeffs() * fwdRatesOfProgress_ddX();
476 Eigen::SparseMatrix<double> Kinetics::destructionRates_ddCi()
478 Eigen::SparseMatrix<double> jac;
480 jac = m_revProductStoich.stoichCoeffs() * revRatesOfProgress_ddCi();
482 jac += m_reactantStoich.stoichCoeffs() * fwdRatesOfProgress_ddCi();
486 void Kinetics::getNetProductionRates_ddT(
double* dwdot)
488 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
489 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
490 getNetRatesOfProgress_ddT(buf.data());
491 out = m_stoichMatrix * buf;
494 void Kinetics::getNetProductionRates_ddP(
double* dwdot)
496 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
497 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
498 getNetRatesOfProgress_ddP(buf.data());
499 out = m_stoichMatrix * buf;
502 void Kinetics::getNetProductionRates_ddC(
double* dwdot)
504 Eigen::Map<Eigen::VectorXd> out(dwdot, m_kk);
505 Eigen::Map<Eigen::VectorXd> buf(m_rbuf.data(), nReactions());
506 getNetRatesOfProgress_ddC(buf.data());
507 out = m_stoichMatrix * buf;
510 Eigen::SparseMatrix<double> Kinetics::netProductionRates_ddX()
512 return m_stoichMatrix * netRatesOfProgress_ddX();
515 Eigen::SparseMatrix<double> Kinetics::netProductionRates_ddCi()
517 return m_stoichMatrix * netRatesOfProgress_ddCi();
520 void Kinetics::addThermo(shared_ptr<ThermoPhase> thermo)
524 if (thermo->nDim() <= m_mindim) {
525 if (!m_thermo.empty()) {
527 "The reacting (lowest dimensional) phase must be added first.");
529 m_mindim = thermo->nDim();
532 m_thermo.push_back(thermo);
533 m_phaseindex[m_thermo.back()->name()] = nPhases();
540 string name = KineticsFactory::factory()->canonicalize(kineticsType());
541 if (name !=
"none") {
542 out[
"kinetics"] = name;
543 if (nReactions() == 0) {
544 out[
"reactions"] =
"none";
546 if (m_hasUndeclaredThirdBodies) {
547 out[
"skip-undeclared-third-bodies"] =
true;
553 void Kinetics::resizeSpecies()
556 m_start.resize(nPhases());
558 for (
size_t i = 0; i < m_thermo.size(); i++) {
560 m_kk += m_thermo[i]->nSpecies();
565 bool Kinetics::addReaction(shared_ptr<Reaction> r,
bool resize)
576 if (!r->checkSpecies(*
this)) {
583 if (r->rate_units.factor() == 0) {
584 r->rate()->setRateUnits(r->calculateRateCoeffUnits(*
this));
587 size_t irxn = nReactions();
590 vector<size_t> rk, pk;
594 vector<double> rstoich, pstoich;
596 for (
const auto& [name, stoich] : r->reactants) {
597 rk.push_back(kineticsSpeciesIndex(name));
598 rstoich.push_back(stoich);
601 for (
const auto& [name, stoich] : r->products) {
602 pk.push_back(kineticsSpeciesIndex(name));
603 pstoich.push_back(stoich);
609 vector<double> rorder = rstoich;
610 for (
const auto& [name, order] : r->orders) {
611 size_t k = kineticsSpeciesIndex(name);
613 auto rloc = std::find(rk.begin(), rk.end(), k);
614 if (rloc != rk.end()) {
615 rorder[rloc - rk.begin()] = order;
622 rstoich.push_back(0.0);
623 rorder.push_back(order);
627 m_reactantStoich.add(irxn, rk, rorder, rstoich);
629 m_productStoich.add(irxn, pk, pstoich, pstoich);
631 m_revProductStoich.add(irxn, pk, pstoich, pstoich);
634 m_reactions.push_back(r);
635 m_rfn.push_back(0.0);
636 m_delta_gibbs0.push_back(0.0);
637 m_rkcn.push_back(0.0);
638 m_ropf.push_back(0.0);
639 m_ropr.push_back(0.0);
640 m_ropnet.push_back(0.0);
641 m_perturb.push_back(1.0);
653 void Kinetics::modifyReaction(
size_t i, shared_ptr<Reaction> rNew)
655 checkReactionIndex(i);
656 shared_ptr<Reaction>& rOld = m_reactions[i];
657 if (rNew->type() != rOld->type()) {
659 "Reaction types are different: {} != {}.",
660 rOld->type(), rNew->type());
663 if (rNew->rate()->type() != rOld->rate()->type()) {
665 "ReactionRate types are different: {} != {}.",
666 rOld->rate()->type(), rNew->rate()->type());
669 if (rNew->reactants != rOld->reactants) {
671 "Reactants are different: '{}' != '{}'.",
672 rOld->reactantString(), rNew->reactantString());
675 if (rNew->products != rOld->products) {
677 "Products are different: '{}' != '{}'.",
678 rOld->productString(), rNew->productString());
680 m_reactions[i] = rNew;
684 shared_ptr<Reaction> Kinetics::reaction(
size_t i)
686 checkReactionIndex(i);
687 return m_reactions[i];
690 shared_ptr<const Reaction> Kinetics::reaction(
size_t i)
const
692 checkReactionIndex(i);
693 return m_reactions[i];
Base class for kinetics managers and also contains the kineticsmgr module documentation (see Kinetics...
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.
Base class for exceptions thrown by Cantera classes.
An array index is out of range.
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
Abstract base class which stores data about a reaction and its rate parameterization so that it can b...
bool usesThirdBody() const
Check whether reaction involves third body collider.
bool reversible
True if the current reaction is reversible. False otherwise.
shared_ptr< ReactionRate > rate()
Get reaction rate pointer.
string type() const
The type of reaction, including reaction rate information.
string equation() const
The chemical equation for this reaction.
shared_ptr< ThirdBody > thirdBody()
Get pointer to third-body handler.
Composition products
Product species and stoichiometric coefficients.
Composition reactants
Reactant species and stoichiometric coefficients.
AnyMap input
Input data used for specific models.
bool duplicate
True if the current reaction is marked as duplicate.
Base class for a phase with thermodynamic properties.
A class for managing third-body efficiencies, including default values.
double efficiency(const string &k) const
Get the third-body efficiency for species k
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...