11BulkKinetics::BulkKinetics() {
18 "To be removed after Cantera 3.0. Use default constructor instead.");
34 for (
const auto& [name, stoich] : r->products) {
37 for (
const auto& [name, stoich] : r->reactants) {
49 shared_ptr<ReactionRate> rate = r->rate();
50 string rtype = rate->subType();
64 rate->setContext(*r, *
this);
71 if (r->thirdBody() !=
nullptr) {
80void BulkKinetics::addThirdBody(shared_ptr<Reaction> r)
82 map<size_t, double> efficiencies;
83 for (
const auto& [name, efficiency] : r->thirdBody()->efficiencies) {
86 efficiencies[k] = efficiency;
88 throw CanteraError(
"BulkKinetics::addThirdBody",
"Found third-body"
89 " efficiency for undefined species '{}' while adding reaction '{}'",
96 r->thirdBody()->default_efficiency,
97 r->thirdBody()->mass_action);
105 shared_ptr<ReactionRate> rate = rNew->rate();
106 string rtype = rate->subType();
108 rtype = rate->type();
114 "Evaluator not available for type '{}'.", rtype);
119 rate->setRateIndex(i);
120 rate->setContext(*rNew, *
this);
144 m_state.resize(
thermo().stateSize());
160void BulkKinetics::invalidateCache()
162 Kinetics::invalidateCache();
179 vector<double>& delta_gibbs0 =
m_rbuf0;
180 fill(delta_gibbs0.begin(), delta_gibbs0.end(), 0.0);
186 double logStandConc = log(
thermo().standardConcentration());
188 kc[i] = exp(-delta_gibbs0[i] * rrt +
m_dn[i] * logStandConc);
198 if (doIrreversible) {
250 for (
size_t k = 0; k <
m_kk; k++) {
263 for (
size_t k = 0; k <
m_kk; k++) {
273 settings[
"skip-falloff"] = m_jac_skip_falloff;
274 settings[
"rtol-delta"] = m_jac_rtol_delta;
279 bool force = settings.
empty();
280 if (force || settings.
hasKey(
"skip-third-bodies")) {
283 if (force || settings.
hasKey(
"skip-falloff")) {
284 m_jac_skip_falloff = settings.
getBool(
"skip-falloff",
false);
286 if (force || settings.
hasKey(
"rtol-delta")) {
287 m_jac_rtol_delta = settings.
getDouble(
"rtol-delta", 1e-8);
310 Eigen::Map<Eigen::VectorXd> dRevRop(drop,
nReactions());
313 Eigen::Map<Eigen::VectorXd> dRevRop2(m_rbuf2.data(),
nReactions());
324 Eigen::Map<Eigen::VectorXd> dNetRop(drop,
nReactions());
327 Eigen::Map<Eigen::VectorXd> dNetRop2(m_rbuf2.data(),
nReactions());
387 Eigen::Map<Eigen::VectorXd> dNetRop(drop,
nReactions());
390 Eigen::Map<Eigen::VectorXd> dNetRop2(m_rbuf2.data(),
nReactions());
399 vector<double>& rop_rates =
m_rbuf0;
409 vector<double>& rop_rates =
m_rbuf0;
420 vector<double>& rop_rates =
m_rbuf0;
434 vector<double>& rop_rates =
m_rbuf0;
444 vector<double>& rop_rates =
m_rbuf0;
455 vector<double>& rop_rates =
m_rbuf0;
464void BulkKinetics::updateROP()
476 double logStandConc = log(
thermo().standardConcentration());
482 for (
size_t i = 0; i <
m_revindex.size(); i++) {
488 for (
size_t i = 0; i !=
m_irrev.size(); ++i) {
493 if (!last.
validate(T, rho, statenum)) {
504 bool changed = rates->update(
thermo(), *
this);
506 rates->getRateConstants(
m_kf0.data());
537 for (
size_t i = 0; i <
m_rfn.size(); i++) {
539 "m_rfn[{}] is not finite.", i);
541 "m_ropf[{}] is not finite.", i);
543 "m_ropr[{}] is not finite.", i);
577 vector<double>& grt = m_sbuf0;
578 vector<double>& delta_gibbs0 = m_rbuf1;
579 fill(delta_gibbs0.begin(), delta_gibbs0.end(), 0.0);
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++) {
596 drkcn[irxn] *= factor;
599 for (
size_t i = 0; i <
m_irrev.size(); ++i) {
609 copy(in.begin(), in.end(), drop);
611 rates->processRateConstants_ddT(drop,
m_rfn.data(), m_jac_rtol_delta);
618 copy(in.begin(), in.end(), drop);
620 rates->processRateConstants_ddP(drop,
m_rfn.data(), m_jac_rtol_delta);
625 double* drop,
bool mass_action)
627 Eigen::Map<Eigen::VectorXd> out(drop,
nReactions());
633 stoich.
scale(in.data(), out.data(), ctot_inv);
640 Eigen::Map<Eigen::VectorXd> outM(m_rbuf1.data(),
nReactions());
648 if (!m_jac_skip_falloff) {
652 rates->processRateConstants_ddM(
653 outM.data(),
m_rfn.data(), m_jac_rtol_delta);
662 Eigen::SparseMatrix<double> out;
663 vector<double>& scaled = m_rbuf1;
664 vector<double>& outV = m_rbuf2;
667 copy(in.begin(), in.end(), scaled.begin());
676 copy(scaled.begin(), scaled.end(), outV.begin());
684 copy(scaled.begin(), scaled.end(), outV.begin());
685 stoich.multiply(
m_act_conc.data(), outV.data());
688 if (!m_jac_skip_falloff) {
691 rates->processRateConstants_ddM(
692 outV.data(),
m_rfn.data(), m_jac_rtol_delta,
false);
703 if (!
thermo().isIdeal()) {
705 "Not supported for non-ideal ThermoPhase models.");
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.
double getDouble(const string &key, double default_) const
If key exists, return it as a double, otherwise return default_.
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
bool empty() const
Return boolean indicating whether AnyMap is empty.
bool getBool(const string &key, bool default_) const
If key exists, return it as a bool, otherwise return default_.
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.
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
ValueCache m_cache
Cache for saved calculations within each Kinetics object.
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.
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...
vector< double > m_rkcn
Reciprocal of the equilibrium constant in concentration units.
size_t m_kk
The number of species in all of the phases that participate in this kinetics mechanism.
virtual void getReactionDelta(const double *property, double *deltaProperty) const
Change in species properties.
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.
bool m_skipUndeclaredThirdBodies
See skipUndeclaredThirdBodies()
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.
StoichManagerN m_revProductStoich
Stoichiometry manager for the products of reversible reactions.
size_t nReactions() const
Number of reactions in the reaction mechanism.
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.
virtual void setMultiplier(size_t i, double f)
Set the multiplier for reaction i to f.
vector< double > m_rfn
Forward rate constant for each reaction.
StoichManagerN m_reactantStoich
Stoichiometry manager for the reactants for each reaction.
virtual void resizeSpecies()
Resize arrays with sizes that depend on the total number of species.
size_t nTotalSpecies() const
The total number of species in all phases participating in the kinetics mechanism.
vector< double > m_delta_gibbs0
Delta G^0 for all reactions.
An error indicating that an unimplemented function has been called.
virtual void getConcentrations(double *const c) const
Get the species concentrations (kmol/m^3).
virtual double molarDensity() const
Molar density (kmol/m^3).
void restoreState(const vector< double > &state)
Restore a state saved on a previous call to saveState.
double temperature() const
Temperature (K).
void saveState(vector< double > &state) const
Save the current internal state of the phase.
virtual double density() const
Density (kg/m^3).
int stateMFNumber() const
Return the State Mole Fraction Number.
virtual double pressure() const
Return the thermodynamic pressure (Pa).
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.
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.
const double GasConstant
Universal Gas Constant [J/kmol/K].
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.
const double BigNumber
largest number to compare to inf.
A cached property value and the state at which it was evaluated.
double state2
Value of the second state variable for the state at which value was evaluated, for example density or...
bool validate(double state1New)
Check whether the currently cached value is valid based on a single state variable.
double state1
Value of the first state variable for the state at which value was evaluated, for example temperature...