19#include <unordered_set>
20#include <boost/algorithm/string.hpp>
30 throw IndexError(
"Kinetics::checkReactionIndex",
"reactions", i,
79 "To be removed after Cantera 3.0. Use reactionPhaseIndex instead.");
88 "Cannot access phases that were not added using smart pointers.");
96 throw IndexError(
"Kinetics::checkSpeciesIndex",
"species", k,
m_kk-1);
110 map<size_t, vector<size_t>> participants;
111 vector<map<int, double>> net_stoich;
112 std::unordered_set<size_t> unmatched_duplicates;
115 unmatched_duplicates.insert(i);
121 unsigned long int key = 0;
123 net_stoich.emplace_back();
124 map<int, double>& net = net_stoich.back();
125 for (
const auto& [name, stoich] : R.
reactants) {
128 net[-1 -k] -= stoich;
130 for (
const auto& [name, stoich] : R.
products) {
137 vector<size_t>& related = participants[key];
138 for (
size_t m : related) {
142 unmatched_duplicates.erase(i);
143 unmatched_duplicates.erase(m);
145 }
else if (R.
type() != other.
type()) {
148 && R.
rate()->type() != other.
rate()->type())
160 bool thirdBodyOk =
true;
177 "Undeclared duplicate reactions detected:\n"
178 "Reaction {}: {}\nReaction {}: {}\n",
184 participants[key].push_back(i);
186 if (unmatched_duplicates.size()) {
187 size_t i = *unmatched_duplicates.begin();
191 "No duplicate found for declared duplicate reaction number {}"
202 std::unordered_set<int> keys;
203 for (
auto& [speciesKey, stoich] : r1) {
204 keys.insert(speciesKey);
206 for (
auto& [speciesKey, stoich] : r2) {
207 keys.insert(speciesKey);
209 int k1 = r1.begin()->first;
212 if (r1[k1] && r2[k1]) {
213 ratio = r2[k1]/r1[k1];
214 bool different =
false;
216 if ((r1[k] && !r2[k]) ||
218 (r1[k] && fabs(r2[k]/r1[k] - ratio) > 1.e-8)) {
229 if (r1[k1] == 0.0 || r2[-k1] == 0.0) {
232 ratio = r2[-k1]/r1[k1];
234 if ((r1[k] && !r2[-k]) ||
235 (!r1[k] && r2[-k]) ||
236 (r1[k] && fabs(r2[-k]/r1[k] - ratio) > 1.e-8)) {
246 warn_deprecated(
"Kinetics::selectPhase",
"To be removed after Cantera 3.0");
247 for (
size_t n = 0; n <
nPhases(); n++) {
251 data +
m_start[n] + nsp, phase_data);
255 throw CanteraError(
"Kinetics::selectPhase",
"Phase not found.");
260 for (
size_t n =
m_start.size()-1; n !=
npos; n--) {
270 for (
size_t n = 0; n <
m_thermo.size(); n++) {
282 warn_deprecated(
"Kinetics::kineticsSpeciesIndex(species_name, phase_name)",
283 "To be removed after Cantera 3.0. Use kineticsSpeciesIndex(species_name).\n"
284 "Species names should be unique across all phases linked to a Kinetics object.");
289 for (
size_t n = 0; n <
m_thermo.size(); n++) {
304 for (
size_t n = 0; n <
m_thermo.size(); n++) {
310 throw CanteraError(
"Kinetics::speciesPhase",
"unknown species '{}'", nm);
320 throw CanteraError(
"Kinetics::speciesPhase",
"unknown species '{}'", nm);
325 for (
size_t n =
m_start.size()-1; n !=
npos; n--) {
331 "illegal species index: {}", k);
345 warn_deprecated(
"Kinetics::reactionType",
"To be removed after Cantera 3.0.");
350 warn_deprecated(
"Kinetics::reactionTypeStr",
"To be removed after Cantera 3.0.");
356 warn_deprecated(
"Kinetics::reactionString",
"To be removed after Cantera 3.0.");
363 warn_deprecated(
"Kinetics::reactantString",
"To be removed after Cantera 3.0.");
370 warn_deprecated(
"Kinetics::productString",
"To be removed after Cantera 3.0.");
394 fill(deltaProp, deltaProp +
nReactions(), 0.0);
403 fill(deltaProp, deltaProp +
nReactions(), 0.0);
415 fill(cdot, cdot +
m_kk, 0.0);
428 fill(ddot, ddot +
m_kk, 0.0);
439 fill(net, net +
m_kk, 0.0);
448 Eigen::Map<Eigen::VectorXd> out(dwdot,
m_kk);
460 Eigen::Map<Eigen::VectorXd> out(dwdot,
m_kk);
472 Eigen::Map<Eigen::VectorXd> out(dwdot,
m_kk);
484 Eigen::SparseMatrix<double> jac;
494 Eigen::SparseMatrix<double> jac;
504 Eigen::Map<Eigen::VectorXd> out(dwdot,
m_kk);
516 Eigen::Map<Eigen::VectorXd> out(dwdot,
m_kk);
528 Eigen::Map<Eigen::VectorXd> out(dwdot,
m_kk);
540 Eigen::SparseMatrix<double> jac;
550 Eigen::SparseMatrix<double> jac;
560 Eigen::Map<Eigen::VectorXd> out(dwdot,
m_kk);
568 Eigen::Map<Eigen::VectorXd> out(dwdot,
m_kk);
576 Eigen::Map<Eigen::VectorXd> out(dwdot,
m_kk);
598 warn_deprecated(
"Kinetics::addThermo",
"The reacting (lowest dimensional) "
599 "phase should be added first. This warning will become an error after "
607 if (boost::algorithm::contains(
thermo->
type(),
"surface")) {
619 "To be removed after Cantera 3.0. Replaced by addThermo.");
640 string name = KineticsFactory::factory()->canonicalize(
kineticsType());
641 if (name !=
"none") {
642 out[
"kinetics"] = name;
644 out[
"reactions"] =
"none";
647 out[
"skip-undeclared-third-bodies"] =
true;
658 for (
size_t i = 0; i <
m_thermo.size(); i++) {
676 if (!r->checkSpecies(*
this)) {
683 if (r->rate_units.factor() == 0) {
684 r->rate()->setRateUnits(r->calculateRateCoeffUnits(*
this));
690 vector<size_t> rk, pk;
694 vector<double> rstoich, pstoich;
696 for (
const auto& [name, stoich] : r->reactants) {
698 rstoich.push_back(stoich);
701 for (
const auto& [name, stoich] : r->products) {
703 pstoich.push_back(stoich);
709 vector<double> rorder = rstoich;
710 for (
const auto& [name, order] : r->orders) {
713 auto rloc = std::find(rk.begin(), rk.end(), k);
714 if (rloc != rk.end()) {
715 rorder[rloc - rk.begin()] = order;
722 rstoich.push_back(0.0);
723 rorder.push_back(order);
735 m_rfn.push_back(0.0);
757 if (rNew->type() != rOld->type()) {
759 "Reaction types are different: {} != {}.",
760 rOld->type(), rNew->type());
763 if (rNew->rate()->type() != rOld->rate()->type()) {
765 "ReactionRate types are different: {} != {}.",
766 rOld->rate()->type(), rNew->rate()->type());
769 if (rNew->reactants != rOld->reactants) {
771 "Reactants are different: '{}' != '{}'.",
772 rOld->reactantString(), rNew->reactantString());
775 if (rNew->products != rOld->products) {
777 "Products are different: '{}' != '{}'.",
778 rOld->productString(), rNew->productString());
798 warn_deprecated(
"Kinetics::reactionEnthalpy",
"To be removed after Cantera 3.0");
800 for (
size_t n = 0; n <
nPhases(); n++) {
803 double rxn_deltaH = 0;
804 for (
const auto& [name, stoich] : products) {
806 rxn_deltaH += hk[k] * stoich;
808 for (
const auto& [name, stoich] : reactants) {
810 rxn_deltaH -= hk[k] * stoich;
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.
virtual void resizeReactions()
Finalize Kinetics object and associated objects.
void checkPhaseIndex(size_t m) const
Check that the specified phase index is in range Throws an exception if m is greater than nPhases()
void checkSpeciesArraySize(size_t mm) const
Check that an array size is at least nSpecies() Throws an exception if kk is less than nSpecies().
virtual void getFwdRatesOfProgress(double *fwdROP)
Return the forward rates of progress of the reactions.
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
vector< shared_ptr< Reaction > > m_reactions
Vector of Reaction objects represented by this Kinetics manager.
void checkSpeciesIndex(size_t k) const
Check that the specified species index is in range Throws an exception if k is greater than nSpecies(...
double checkDuplicateStoich(map< int, double > &r1, map< int, double > &r2) const
Check whether r1 and r2 represent duplicate stoichiometries This function returns a ratio if two reac...
vector< double > m_perturb
Vector of perturbation factors for each reaction's rate of progress vector.
vector< double > m_ropf
Forward rate-of-progress for each reaction.
shared_ptr< Reaction > reaction(size_t i)
Return the Reaction object for reaction i.
bool m_ready
Boolean indicating whether Kinetics object is fully configured.
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...
virtual void getNetRatesOfProgress(double *netROP)
Net rates of progress.
vector< size_t > m_start
m_start is a vector of integers specifying the beginning position for the species vector for the n'th...
virtual string kineticsType() const
Identifies the Kinetics manager type.
string reactionString(size_t i) const
Return a string representing the reaction.
vector< double > m_rkcn
Reciprocal of the equilibrium constant in concentration units.
void selectPhase(const double *data, const ThermoPhase *phase, double *phase_data)
Takes as input an array of properties for all species in the mechanism and copies those values belong...
virtual string reactionTypeStr(size_t i) const
String specifying the type of reaction.
vector< ThermoPhase * > m_thermo
m_thermo is a vector of pointers to ThermoPhase objects that are involved with this kinetics operator
shared_ptr< ThermoPhase > reactionPhase() const
Return pointer to phase where the reactions occur.
size_t m_kk
The number of species in all of the phases that participate in this kinetics mechanism.
size_t m_rxnphase
Phase Index where reactions are assumed to be taking place.
virtual pair< size_t, size_t > checkDuplicates(bool throw_err=true) const
Check for unmarked duplicate reactions and unmatched marked duplicates.
virtual void getReactionDelta(const double *property, double *deltaProperty) const
Change in species properties.
void checkPhaseArraySize(size_t mm) const
Check that an array size is at least nPhases() Throws an exception if mm is less than nPhases().
vector< double > m_dH
The enthalpy change for each reaction to calculate Blowers-Masel rates.
vector< double > m_ropr
Reverse rate-of-progress for each reaction.
size_t nPhases() const
The number of phases participating in the reaction mechanism.
virtual void addPhase(ThermoPhase &thermo)
virtual bool addReaction(shared_ptr< Reaction > r, bool resize=true)
Add a single reaction to the mechanism.
virtual string reactionType(size_t i) const
String specifying the type of reaction.
string kineticsSpeciesName(size_t k) const
Return the name of the kth species in the kinetics manager.
virtual void getDestructionRates(double *ddot)
Species destruction rates [kmol/m^3/s or kmol/m^2/s].
AnyMap parameters()
Return the parameters for a phase definition which are needed to reconstruct an identical object usin...
virtual void init()
Prepare the class for the addition of reactions, after all phases have been added.
string reactantString(size_t i) const
Returns a string containing the reactants side of the reaction equation.
virtual void modifyReaction(size_t i, shared_ptr< Reaction > rNew)
Modify the rate expression associated with a reaction.
vector< double > m_ropnet
Net rate-of-progress for each reaction.
virtual double productStoichCoeff(size_t k, size_t i) const
Stoichiometric coefficient of species k as a product in reaction i.
virtual void addThermo(shared_ptr< ThermoPhase > thermo)
Add a phase to the kinetics manager object.
vector< double > m_rbuf
Buffer used for storage of intermediate reaction-specific results.
Eigen::SparseMatrix< double > m_stoichMatrix
Net stoichiometry (products - reactants)
map< string, size_t > m_phaseindex
Mapping of the phase name to the position of the phase within the kinetics object.
size_t m_mindim
number of spatial dimensions of lowest-dimensional phase.
StoichManagerN m_productStoich
Stoichiometry manager for the products for each reaction.
StoichManagerN m_revProductStoich
Stoichiometry manager for the products of reversible reactions.
virtual double reactionEnthalpy(const Composition &reactants, const Composition &products)
Calculate the reaction enthalpy of a reaction which has not necessarily been added into the Kinetics ...
size_t nReactions() const
Number of reactions in the reaction mechanism.
virtual void getRevRatesOfProgress(double *revROP)
Return the Reverse rates of progress of the reactions.
bool m_hasUndeclaredThirdBodies
Flag indicating whether reactions include undeclared third bodies.
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
size_t m_surfphase
Index in the list of phases of the one surface phase.
vector< double > m_rfn
Forward rate constant for each reaction.
StoichManagerN m_reactantStoich
Stoichiometry manager for the reactants for each reaction.
size_t surfacePhaseIndex() const
This returns the integer index of the phase which has ThermoPhase type cSurf.
virtual void resizeSpecies()
Resize arrays with sizes that depend on the total number of species.
void checkReactionArraySize(size_t ii) const
Check that an array size is at least nReactions() Throws an exception if ii is less than nReactions()...
vector< shared_ptr< ThermoPhase > > m_sharedThermo
vector of shared pointers,
size_t nTotalSpecies() const
The total number of species in all phases participating in the kinetics mechanism.
virtual double reactantStoichCoeff(size_t k, size_t i) const
Stoichiometric coefficient of species k as a reactant in reaction i.
string productString(size_t i) const
Returns a string containing the products side of the reaction equation.
void checkReactionIndex(size_t m) const
Check that the specified reaction index is in range Throws an exception if i is greater than nReactio...
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
ThermoPhase & speciesPhase(const string &nm)
This function looks up the name of a species and returns a reference to the ThermoPhase object of the...
vector< double > m_delta_gibbs0
Delta G^0 for all reactions.
virtual void getCreationRates(double *cdot)
Species creation rates [kmol/m^3/s or kmol/m^2/s].
size_t speciesPhaseIndex(size_t k) const
This function takes as an argument the kineticsSpecies index (that is, the list index in the list of ...
size_t nSpecies() const
Returns the number of species in the phase.
size_t nDim() const
Returns the number of spatial dimensions (1, 2, or 3)
string speciesName(size_t k) const
Name of the species with index k.
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
string name() const
Return the name of the phase.
Abstract base class which stores data about a reaction and its rate parameterization so that it can b...
shared_ptr< ReactionRate > rate()
Get reaction rate pointer.
bool usesThirdBody() const
Check whether reaction involves third body collider.
shared_ptr< ThirdBody > thirdBody()
Get pointer to third-body handler.
bool reversible
True if the current reaction is reversible. False otherwise.
string type() const
The type of reaction, including reaction rate information.
string equation() const
The chemical equation for this reaction.
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.
void resizeCoeffs(size_t nSpc, size_t nRxn)
Resize the sparse coefficient matrix.
void add(size_t rxn, const vector< size_t > &k)
Add a single reaction to the list of reactions that this stoichiometric manager object handles.
const Eigen::SparseMatrix< double > & stoichCoeffs() const
Return matrix containing stoichiometric coefficients.
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.
string type() const override
String indicating the thermodynamic model implemented.
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,...
virtual Eigen::SparseMatrix< double > fwdRatesOfProgress_ddCi()
Calculate derivatives for forward rates-of-progress with respect to species concentration at constant...
Eigen::SparseMatrix< double > creationRates_ddCi()
Calculate derivatives for species creation rates with respect to species concentration at constant te...
virtual Eigen::SparseMatrix< double > netRatesOfProgress_ddCi()
Calculate derivatives for net rates-of-progress with respect to species concentration at constant tem...
void getCreationRates_ddT(double *dwdot)
Calculate derivatives for species creation rates with respect to temperature at constant pressure,...
virtual Eigen::SparseMatrix< double > netRatesOfProgress_ddX()
Calculate derivatives for net rates-of-progress with respect to species mole fractions at constant te...
Eigen::SparseMatrix< double > destructionRates_ddX()
Calculate derivatives for species destruction rates with respect to species mole fractions at constan...
void getCreationRates_ddC(double *dwdot)
Calculate derivatives for species creation rates with respect to molar concentration at constant temp...
void getDestructionRates_ddP(double *dwdot)
Calculate derivatives for species destruction rates with respect to pressure at constant temperature,...
void getNetProductionRates_ddC(double *dwdot)
Calculate derivatives for species net production rates with respect to molar concentration at constan...
void getDestructionRates_ddT(double *dwdot)
Calculate derivatives for species destruction rates with respect to temperature at constant pressure,...
virtual void getFwdRatesOfProgress_ddP(double *drop)
Calculate derivatives for forward rates-of-progress with respect to pressure at constant temperature,...
Eigen::SparseMatrix< double > netProductionRates_ddX()
Calculate derivatives for species net production rates with respect to species mole fractions at cons...
Eigen::SparseMatrix< double > creationRates_ddX()
Calculate derivatives for species creation rates with respect to species mole fractions at constant t...
Eigen::SparseMatrix< double > destructionRates_ddCi()
Calculate derivatives for species destruction rates with respect to species concentration at constant...
void getNetProductionRates_ddT(double *dwdot)
Calculate derivatives for species net production rates with respect to temperature at constant pressu...
virtual Eigen::SparseMatrix< double > fwdRatesOfProgress_ddX()
Calculate derivatives for forward rates-of-progress with respect to species mole fractions at constan...
virtual void getRevRatesOfProgress_ddT(double *drop)
Calculate derivatives for reverse rates-of-progress with respect to temperature at constant pressure,...
void getCreationRates_ddP(double *dwdot)
Calculate derivatives for species creation rates with respect to pressure at constant temperature,...
virtual void getNetRatesOfProgress_ddT(double *drop)
Calculate derivatives for net rates-of-progress with respect to temperature at constant pressure,...
virtual void getRevRatesOfProgress_ddP(double *drop)
Calculate derivatives for reverse rates-of-progress with respect to pressure at constant temperature,...
Eigen::SparseMatrix< double > netProductionRates_ddCi()
Calculate derivatives for species net production rates with respect to species concentration at const...
virtual void getRevRatesOfProgress_ddC(double *drop)
Calculate derivatives for reverse rates-of-progress with respect to molar concentration at constant t...
void getNetProductionRates_ddP(double *dwdot)
Calculate derivatives for species net production rates with respect to pressure at constant temperatu...
virtual void getNetRatesOfProgress_ddP(double *drop)
Calculate derivatives for net rates-of-progress with respect to pressure at constant temperature,...
virtual Eigen::SparseMatrix< double > revRatesOfProgress_ddCi()
Calculate derivatives for forward rates-of-progress with respect to species concentration at constant...
virtual void getFwdRatesOfProgress_ddT(double *drop)
Calculate derivatives for forward rates-of-progress with respect to temperature at constant pressure,...
void getDestructionRates_ddC(double *dwdot)
Calculate derivatives for species destruction rates with respect to molar concentration at constant t...
virtual void getNetRatesOfProgress_ddC(double *drop)
Calculate derivatives for net rates-of-progress with respect to molar concentration at constant tempe...
virtual Eigen::SparseMatrix< double > revRatesOfProgress_ddX()
Calculate derivatives for reverse rates-of-progress with respect to species mole fractions at constan...
virtual void getFwdRatesOfProgress_ddC(double *drop)
Calculate derivatives for forward rates-of-progress with respect to molar concentration at constant t...
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.
map< string, double > Composition
Map from string names to doubles.
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...