11 BulkKinetics::BulkKinetics() {
27 for (
const auto& [name, stoich] : r->products) {
30 for (
const auto& [name, stoich] : r->reactants) {
42 shared_ptr<ReactionRate> rate = r->rate();
43 string rtype = rate->subType();
57 rate->setContext(*r, *
this);
64 if (r->thirdBody() !=
nullptr) {
73 void BulkKinetics::addThirdBody(shared_ptr<Reaction> r)
75 map<size_t, double> efficiencies;
76 for (
const auto& [name, efficiency] : r->thirdBody()->efficiencies) {
79 efficiencies[k] = efficiency;
81 throw CanteraError(
"BulkKinetics::addThirdBody",
"Found third-body"
82 " efficiency for undefined species '{}' while adding reaction '{}'",
89 r->thirdBody()->default_efficiency,
90 r->thirdBody()->mass_action);
98 shared_ptr<ReactionRate> rate = rNew->rate();
99 string rtype = rate->subType();
101 rtype = rate->type();
107 "Evaluator not available for type '{}'.", rtype);
112 rate->setRateIndex(i);
113 rate->setContext(*rNew, *
this);
137 m_state.resize(
thermo().stateSize());
153 void BulkKinetics::invalidateCache()
155 Kinetics::invalidateCache();
172 vector<double>& delta_gibbs0 =
m_rbuf0;
173 fill(delta_gibbs0.begin(), delta_gibbs0.end(), 0.0);
179 double logStandConc = log(
thermo().standardConcentration());
181 kc[i] = exp(-delta_gibbs0[i] * rrt +
m_dn[i] * logStandConc);
191 if (doIrreversible) {
243 for (
size_t k = 0; k <
m_kk; k++) {
256 for (
size_t k = 0; k <
m_kk; k++) {
266 settings[
"skip-falloff"] = m_jac_skip_falloff;
267 settings[
"rtol-delta"] = m_jac_rtol_delta;
272 bool force = settings.
empty();
273 if (force || settings.
hasKey(
"skip-third-bodies")) {
276 if (force || settings.
hasKey(
"skip-falloff")) {
277 m_jac_skip_falloff = settings.
getBool(
"skip-falloff",
false);
279 if (force || settings.
hasKey(
"rtol-delta")) {
280 m_jac_rtol_delta = settings.
getDouble(
"rtol-delta", 1e-8);
303 Eigen::Map<Eigen::VectorXd> dRevRop(drop,
nReactions());
306 Eigen::Map<Eigen::VectorXd> dRevRop2(m_rbuf2.data(),
nReactions());
317 Eigen::Map<Eigen::VectorXd> dNetRop(drop,
nReactions());
320 Eigen::Map<Eigen::VectorXd> dNetRop2(m_rbuf2.data(),
nReactions());
380 Eigen::Map<Eigen::VectorXd> dNetRop(drop,
nReactions());
383 Eigen::Map<Eigen::VectorXd> dNetRop2(m_rbuf2.data(),
nReactions());
392 vector<double>& rop_rates =
m_rbuf0;
402 vector<double>& rop_rates =
m_rbuf0;
413 vector<double>& rop_rates =
m_rbuf0;
427 vector<double>& rop_rates =
m_rbuf0;
437 vector<double>& rop_rates =
m_rbuf0;
448 vector<double>& rop_rates =
m_rbuf0;
457 void BulkKinetics::updateROP()
469 double logStandConc = log(
thermo().standardConcentration());
475 for (
size_t i = 0; i <
m_revindex.size(); i++) {
481 for (
size_t i = 0; i !=
m_irrev.size(); ++i) {
486 if (!last.
validate(T, rho, statenum)) {
497 bool changed = rates->update(
thermo(), *
this);
499 rates->getRateConstants(
m_kf0.data());
530 for (
size_t i = 0; i <
m_rfn.size(); i++) {
532 "m_rfn[{}] is not finite.", i);
534 "m_ropf[{}] is not finite.", i);
536 "m_ropr[{}] is not finite.", i);
570 vector<double>& grt = m_sbuf0;
571 vector<double>& delta_gibbs0 = m_rbuf1;
572 fill(delta_gibbs0.begin(), delta_gibbs0.end(), 0.0);
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++) {
589 drkcn[irxn] *= factor;
592 for (
size_t i = 0; i <
m_irrev.size(); ++i) {
602 copy(in.begin(), in.end(), drop);
604 rates->processRateConstants_ddT(drop,
m_rfn.data(), m_jac_rtol_delta);
611 copy(in.begin(), in.end(), drop);
613 rates->processRateConstants_ddP(drop,
m_rfn.data(), m_jac_rtol_delta);
618 double* drop,
bool mass_action)
620 Eigen::Map<Eigen::VectorXd> out(drop,
nReactions());
626 stoich.
scale(in.data(), out.data(), ctot_inv);
633 Eigen::Map<Eigen::VectorXd> outM(m_rbuf1.data(),
nReactions());
641 if (!m_jac_skip_falloff) {
645 rates->processRateConstants_ddM(
646 outM.data(),
m_rfn.data(), m_jac_rtol_delta);
655 Eigen::SparseMatrix<double> out;
656 vector<double>& scaled = m_rbuf1;
657 vector<double>& outV = m_rbuf2;
660 copy(in.begin(), in.end(), scaled.begin());
669 copy(scaled.begin(), scaled.end(), outV.begin());
677 copy(scaled.begin(), scaled.end(), outV.begin());
678 stoich.multiply(
m_act_conc.data(), outV.data());
681 if (!m_jac_skip_falloff) {
684 rates->processRateConstants_ddM(
685 outV.data(),
m_rfn.data(), m_jac_rtol_delta,
false);
696 if (!
thermo().isIdeal()) {
698 "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_.
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.
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 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.
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
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.
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.
Eigen::SparseMatrix< double > derivatives(const double *product)
Calculate derivatives with respect to species concentrations.
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.
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"
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...