Skip to content

Commit

Permalink
[Kinetics] Consolidate update of rate coefficients and rates of progress
Browse files Browse the repository at this point in the history
  • Loading branch information
@speth @ischoegl
speth authored and ischoegl committed Apr 21, 2023
1 parent b5d4565 commit 80dc64c
Show file tree
Hide file tree
Showing 3 changed files with 69 additions and 96 deletions.
19 changes: 1 addition & 18 deletions include/cantera/kinetics/BulkKinetics.h
View file
Expand Up @@ -88,16 +88,6 @@ class BulkKinetics : public Kinetics
//! @name Rate calculation intermediate methods
//! @{

//! Update temperature-dependent portions of reaction rates and falloff
//! functions.
virtual void update_rates_T();

//! Update properties that depend on concentrations.
//! Currently the enhanced collision partner concentrations are updated
//! here, as well as the pressure-dependent portion of P-log and Chebyshev
//! reactions.
virtual void update_rates_C();

void updateROP() override;

void getThirdBodyConcentrations(double* concm) override;
Expand All @@ -115,15 +105,9 @@ class BulkKinetics : public Kinetics
//! their derivatives.
//! @{

//! Calculate rate coefficients
void processFwdRateCoefficients(double* ropf);

//! Multiply rate with third-body collider concentrations
void processThirdBodies(double* rop);

//! Update the equilibrium constants in molar units.
void updateKc();

//! Multiply rate with inverse equilibrium constant
void applyEquilibriumConstants(double* rop);

Expand Down Expand Up @@ -197,13 +181,12 @@ class BulkKinetics : public Kinetics
double m_jac_rtol_delta;

bool m_ROP_ok = false;
double m_temp = 0.0;
double m_logStandConc = 0.0;

//! Buffers for partial rop results with length nReactions()
vector_fp m_rbuf0;
vector_fp m_rbuf1;
vector_fp m_rbuf2;
vector_fp m_kf0; //!< Forward rate constants without perturbation
vector_fp m_sbuf0;
vector_fp m_state;
vector_fp m_grt;
Expand Down
4 changes: 3 additions & 1 deletion include/cantera/kinetics/Kinetics.h
View file
Expand Up @@ -1347,7 +1347,9 @@ class Kinetics
m_perturb[i] = f;
}

virtual void invalidateCache() {};
virtual void invalidateCache() {
m_cache.clear();
};

//! @}
//! Check for unmarked duplicate reactions and unmatched marked duplicates
Expand Down
142 changes: 65 additions & 77 deletions src/kinetics/BulkKinetics.cpp
View file
Expand Up @@ -137,6 +137,7 @@ void BulkKinetics::resizeReactions()
m_rbuf0.resize(nReactions());
m_rbuf1.resize(nReactions());
m_rbuf2.resize(nReactions());
m_kf0.resize(nReactions());
m_sbuf0.resize(nTotalSpecies());
m_state.resize(thermo().stateSize());
m_multi_concm.resizeCoeffs(nTotalSpecies(), nReactions());
Expand All @@ -158,24 +159,20 @@ void BulkKinetics::invalidateCache()
{
Kinetics::invalidateCache();
m_ROP_ok = false;
m_temp += 0.13579;
}

void BulkKinetics::getFwdRateConstants(double* kfwd)
{
processFwdRateCoefficients(m_ropf.data());

updateROP();
copy(m_rfn.begin(), m_rfn.end(), kfwd);
if (legacy_rate_constants_used()) {
processThirdBodies(m_ropf.data());
processThirdBodies(kfwd);
}

// copy result
copy(m_ropf.begin(), m_ropf.end(), kfwd);
}

void BulkKinetics::getEquilibriumConstants(double* kc)
{
update_rates_T(); // this step ensures that m_grt is updated
updateROP();

vector_fp& delta_gibbs0 = m_rbuf0;
fill(delta_gibbs0.begin(), delta_gibbs0.end(), 0.0);
Expand All @@ -184,8 +181,9 @@ void BulkKinetics::getEquilibriumConstants(double* kc)
getReactionDelta(m_grt.data(), delta_gibbs0.data());

double rrt = 1.0 / thermo().RT();
double logStandConc = log(thermo().standardConcentration());
for (size_t i = 0; i < nReactions(); i++) {
kc[i] = exp(-delta_gibbs0[i] * rrt + m_dn[i] * m_logStandConc);
kc[i] = exp(-delta_gibbs0[i] * rrt + m_dn[i] * logStandConc);
}
}

Expand All @@ -196,9 +194,9 @@ void BulkKinetics::getRevRateConstants(double* krev, bool doIrreversible)
getFwdRateConstants(krev);

if (doIrreversible) {
getEquilibriumConstants(m_ropnet.data());
getEquilibriumConstants(m_rbuf0.data());
for (size_t i = 0; i < nReactions(); i++) {
krev[i] /= m_ropnet[i];
krev[i] /= m_rbuf0[i];
}
} else {
// m_rkcn[] is zero for irreversible reactions
Expand Down Expand Up @@ -397,7 +395,7 @@ Eigen::SparseMatrix<double> BulkKinetics::fwdRatesOfProgress_ddX()

// forward reaction rate coefficients
vector_fp& rop_rates = m_rbuf0;
processFwdRateCoefficients(rop_rates.data());
getFwdRateConstants(rop_rates.data());
return calculateCompositionDerivatives(m_reactantStoich, rop_rates);
}

Expand All @@ -407,7 +405,7 @@ Eigen::SparseMatrix<double> BulkKinetics::revRatesOfProgress_ddX()

// reverse reaction rate coefficients
vector_fp& rop_rates = m_rbuf0;
processFwdRateCoefficients(rop_rates.data());
getFwdRateConstants(rop_rates.data());
applyEquilibriumConstants(rop_rates.data());
return calculateCompositionDerivatives(m_revProductStoich, rop_rates);
}
Expand All @@ -418,7 +416,7 @@ Eigen::SparseMatrix<double> BulkKinetics::netRatesOfProgress_ddX()

// forward reaction rate coefficients
vector_fp& rop_rates = m_rbuf0;
processFwdRateCoefficients(rop_rates.data());
getFwdRateConstants(rop_rates.data());
auto jac = calculateCompositionDerivatives(m_reactantStoich, rop_rates);

// reverse reaction rate coefficients
Expand All @@ -432,7 +430,7 @@ Eigen::SparseMatrix<double> BulkKinetics::fwdRatesOfProgress_ddCi()

// forward reaction rate coefficients
vector_fp& rop_rates = m_rbuf0;
processFwdRateCoefficients(rop_rates.data());
getFwdRateConstants(rop_rates.data());
return calculateCompositionDerivatives(m_reactantStoich, rop_rates, false);
}

Expand All @@ -442,7 +440,7 @@ Eigen::SparseMatrix<double> BulkKinetics::revRatesOfProgress_ddCi()

// reverse reaction rate coefficients
vector_fp& rop_rates = m_rbuf0;
processFwdRateCoefficients(rop_rates.data());
getFwdRateConstants(rop_rates.data());
applyEquilibriumConstants(rop_rates.data());
return calculateCompositionDerivatives(m_revProductStoich, rop_rates, false);
}
Expand All @@ -453,49 +451,74 @@ Eigen::SparseMatrix<double> BulkKinetics::netRatesOfProgress_ddCi()

// forward reaction rate coefficients
vector_fp& rop_rates = m_rbuf0;
processFwdRateCoefficients(rop_rates.data());
getFwdRateConstants(rop_rates.data());
auto jac = calculateCompositionDerivatives(m_reactantStoich, rop_rates, false);

// reverse reaction rate coefficients
applyEquilibriumConstants(rop_rates.data());
return jac - calculateCompositionDerivatives(m_revProductStoich, rop_rates, false);
}

void BulkKinetics::update_rates_T()
void BulkKinetics::updateROP()
{
static const int cacheId = m_cache.getId();
CachedScalar last = m_cache.getScalar(cacheId);
double T = thermo().temperature();
m_logStandConc = log(thermo().standardConcentration());
double rho = thermo().density();
int statenum = thermo().stateMFNumber();

if (last.state1 != T || last.state2 != rho) {
// Update properties that are independent of the composition
thermo().getStandardChemPotentials(m_grt.data());
fill(m_delta_gibbs0.begin(), m_delta_gibbs0.end(), 0.0);
double logStandConc = log(thermo().standardConcentration());

// compute Delta G^0 for all reversible reactions
getRevReactionDelta(m_grt.data(), m_delta_gibbs0.data());

double rrt = 1.0 / thermo().RT();
for (size_t i = 0; i < m_revindex.size(); i++) {
size_t irxn = m_revindex[i];
m_rkcn[irxn] = std::min(
exp(m_delta_gibbs0[irxn] * rrt - m_dn[irxn] * logStandConc), BigNumber);
}

if (T != m_temp) {
updateKc();
m_ROP_ok = false;
for (size_t i = 0; i != m_irrev.size(); ++i) {
m_rkcn[ m_irrev[i] ] = 0.0;
}
}

// loop over MultiRate evaluators for each reaction type
for (auto& rates : m_bulk_rates) {
bool changed = rates->update(thermo(), *this);
if (changed) {
rates->getRateConstants(m_rfn.data());
m_ROP_ok = false;
if (!last.validate(T, rho, statenum)) {
// Update terms dependent on species concentrations and temperature
thermo().getActivityConcentrations(m_act_conc.data());
thermo().getConcentrations(m_phys_conc.data());
double ctot = thermo().molarDensity();

// Third-body objects interacting with MultiRate evaluator
m_multi_concm.update(m_phys_conc, ctot, m_concm.data());

// loop over MultiRate evaluators for each reaction type
for (auto& rates : m_bulk_rates) {
bool changed = rates->update(thermo(), *this);
if (changed) {
rates->getRateConstants(m_kf0.data());
}
}
m_ROP_ok = false;
}
m_temp = T;
}

void BulkKinetics::update_rates_C()
{
thermo().getActivityConcentrations(m_act_conc.data());
thermo().getConcentrations(m_phys_conc.data());
double ctot = thermo().molarDensity();
if (m_ROP_ok) {
// rates of progress are up-to-date only if both the thermodynamic state
// and m_perturb are unchanged
return;
}

// Third-body objects interacting with MultiRate evaluator
m_multi_concm.update(m_phys_conc, ctot, m_concm.data());
m_ROP_ok = false;
}
// Scale the forward rate coefficient by the perturbation factor
for (size_t i = 0; i < nReactions(); ++i) {
m_rfn[i] = m_kf0[i] * m_perturb[i];
}

void BulkKinetics::updateROP()
{
processFwdRateCoefficients(m_ropf.data());
copy(m_rfn.begin(), m_rfn.end(), m_ropf.data());
processThirdBodies(m_ropf.data());
copy(m_ropf.begin(), m_ropf.end(), m_ropr.begin());

Expand Down Expand Up @@ -526,21 +549,6 @@ void BulkKinetics::getThirdBodyConcentrations(double* concm)
std::copy(m_concm.begin(), m_concm.end(), concm);
}


void BulkKinetics::processFwdRateCoefficients(double* ropf)
{
update_rates_C();
update_rates_T();

// copy rate coefficients into ropf
copy(m_rfn.begin(), m_rfn.end(), ropf);

// Scale the forward rate coefficient by the perturbation factor
for (size_t i = 0; i < nReactions(); ++i) {
ropf[i] *= m_perturb[i];
}
}

void BulkKinetics::processThirdBodies(double* rop)
{
// reactions involving third body
Expand All @@ -549,26 +557,6 @@ void BulkKinetics::processThirdBodies(double* rop)
}
}

void BulkKinetics::updateKc()
{
thermo().getStandardChemPotentials(m_grt.data());
fill(m_delta_gibbs0.begin(), m_delta_gibbs0.end(), 0.0);

// compute Delta G^0 for all reversible reactions
getRevReactionDelta(m_grt.data(), m_delta_gibbs0.data());

double rrt = 1.0 / thermo().RT();
for (size_t i = 0; i < m_revindex.size(); i++) {
size_t irxn = m_revindex[i];
m_rkcn[irxn] = std::min(
exp(m_delta_gibbs0[irxn] * rrt - m_dn[irxn] * m_logStandConc), BigNumber);
}

for (size_t i = 0; i != m_irrev.size(); ++i) {
m_rkcn[ m_irrev[i] ] = 0.0;
}
}

void BulkKinetics::applyEquilibriumConstants(double* rop)
{
// For reverse rates computed from thermochemistry, multiply the forward
Expand Down

0 comments on commit 80dc64c

Please sign in to comment.