Skip to content

Commit

Permalink
Comment changes and typo corrections
Browse files Browse the repository at this point in the history
  • Loading branch information
anthony-walker authored and speth committed Apr 11, 2023
1 parent 199ee2a commit 05055e1
Show file tree
Hide file tree
Showing 13 changed files with 63 additions and 64 deletions.
2 changes: 1 addition & 1 deletion include/cantera/kinetics/GasKinetics.h
Expand Up @@ -139,7 +139,7 @@ class GasKinetics : public BulkKinetics
//! @param ddX true: w.r.t mole fractions false: w.r.t species concentrations
//! @return a sparse matrix of derivative contributions for each reaction of
//! dimensions nTotalReactions by nTotalSpecies
Eigen::SparseMatrix<double> process_derivatives(StoichManagerN& stoich,
Eigen::SparseMatrix<double> calculateCompositionDerivatives(StoichManagerN& stoich,
const vector_fp& in, bool ddX=true);

//! Helper function ensuring that all rate derivatives can be calculated
Expand Down
14 changes: 11 additions & 3 deletions include/cantera/kinetics/InterfaceKinetics.h
Expand Up @@ -307,7 +307,15 @@ class InterfaceKinetics : public Kinetics

/**
* Set/modify derivative settings. For the interface kinetics object this
* includes `"skip-cov-dep"`, `"skip-electrochem"`, and `"rtol-delta"`.
* includes `"skip-coverage-dependence"`, `"skip-electrochem"`, and `"rtol-delta"`.
*
* - `skip-coverage-dependence` (boolean) ... if `true` (default), rate constant
* coverage dependence is not considered when evaluating derivatives.
* - `skip-electrochem` (boolean) ... if `true` (default), electrical charge is
* not considered in evaluating the derivatives and these reactions are treated
* as normal surface reactions.
* - `rtol-delta` (double) ... relative tolerance used to perturb properties
* when calculating numerical derivatives. The default value is 1e-8.
*
* @param settings AnyMap containing settings determining derivative evaluation.
*/
Expand Down Expand Up @@ -335,14 +343,14 @@ class InterfaceKinetics : public Kinetics
//! @param in rate expression used for the derivative calculation
//! @return a sparse matrix of derivative contributions for each reaction of
//! dimensions nTotalReactions by nTotalSpecies
Eigen::SparseMatrix<double> process_derivatives(StoichManagerN& stoich,
Eigen::SparseMatrix<double> calculateCompositionDerivatives(StoichManagerN& stoich,
const vector_fp& in);

//! Helper function ensuring that all rate derivatives can be calculated
//! @param name method name used for error output
//! @throw CanteraError if coverage dependence or electrochemical reactions are
//! included
void assertDerivativesValid(const std::string& name);
void assertDerivativesValid(const string& name);

//! @}

Expand Down
8 changes: 0 additions & 8 deletions include/cantera/kinetics/InterfaceRate.h
Expand Up @@ -44,11 +44,6 @@ struct InterfaceData : public BlowersMaselData

virtual void update(double T) override;

virtual void update(double T, double P) override {
BlowersMaselData::update(T);
pressure = P;
}

virtual void update(double T, const vector_fp& values) override;

using BlowersMaselData::update;
Expand All @@ -72,9 +67,6 @@ struct InterfaceData : public BlowersMaselData
vector_fp electricPotentials; //!< electric potentials of phases
vector_fp standardChemPotentials; //!< standard state chemical potentials
vector_fp standardConcentrations; //!< standard state concentrations

protected:
double m_pressure_buf; //!< buffered pressure
};


Expand Down
2 changes: 1 addition & 1 deletion include/cantera/zeroD/MoleReactor.h
Expand Up @@ -39,7 +39,7 @@ class MoleReactor : public Reactor
protected:
//! Map all surface chemistry derivatives from concentration jacobian to
//! A vector of triplets for the reactor
virtual void mapSurfJacobian(std::vector<Eigen::Triplet<double>> &strips);
virtual void addSurfaceJacobian(vector<Eigen::Triplet<double>> &triplets);

//! Get moles of the system from mass fractions stored by thermo object
//! @param y vector for moles to be put into
Expand Down
2 changes: 1 addition & 1 deletion include/cantera/zeroD/ReactorNet.h
Expand Up @@ -290,7 +290,7 @@ class ReactorNet : public FuncEval
//!
//! @warning This method is an experimental part of the %Cantera
//! API and may be changed or removed without notice.
Eigen::SparseMatrix<double> finiteDifferenceJacobian();
virtual Eigen::SparseMatrix<double> finiteDifferenceJacobian();

//! Method to calculate the system jacobian
//! @warning Depending on the particular implementation, this may return an
Expand Down
3 changes: 3 additions & 0 deletions interfaces/cython/cantera/preconditioners.pyx
Expand Up @@ -96,6 +96,9 @@ cdef class AdaptivePreconditioner(PreconditionerBase):
self.preconditioner.printPreconditioner()

property matrix:
"""
Property to retrieve the latest internal preconditioner matrix.
"""
def __get__(self):
cdef CxxSparseMatrix smat = self.preconditioner.matrix()
return get_from_sparse(smat, smat.rows(), smat.cols())
18 changes: 9 additions & 9 deletions src/kinetics/GasKinetics.cpp
Expand Up @@ -403,7 +403,7 @@ void GasKinetics::getNetRatesOfProgress_ddC(double* drop)
dNetRop -= dNetRop2;
}

Eigen::SparseMatrix<double> GasKinetics::process_derivatives(
Eigen::SparseMatrix<double> GasKinetics::calculateCompositionDerivatives(
StoichManagerN& stoich, const vector_fp& in, bool ddX)
{
Eigen::SparseMatrix<double> out;
Expand Down Expand Up @@ -453,7 +453,7 @@ Eigen::SparseMatrix<double> GasKinetics::fwdRatesOfProgress_ddX()
// forward reaction rate coefficients
vector_fp& rop_rates = m_rbuf0;
processFwdRateCoefficients(rop_rates.data());
return process_derivatives(m_reactantStoich, rop_rates);
return calculateCompositionDerivatives(m_reactantStoich, rop_rates);
}

Eigen::SparseMatrix<double> GasKinetics::revRatesOfProgress_ddX()
Expand All @@ -464,7 +464,7 @@ Eigen::SparseMatrix<double> GasKinetics::revRatesOfProgress_ddX()
vector_fp& rop_rates = m_rbuf0;
processFwdRateCoefficients(rop_rates.data());
processEquilibriumConstants(rop_rates.data());
return process_derivatives(m_revProductStoich, rop_rates);
return calculateCompositionDerivatives(m_revProductStoich, rop_rates);
}

Eigen::SparseMatrix<double> GasKinetics::netRatesOfProgress_ddX()
Expand All @@ -474,11 +474,11 @@ Eigen::SparseMatrix<double> GasKinetics::netRatesOfProgress_ddX()
// forward reaction rate coefficients
vector_fp& rop_rates = m_rbuf0;
processFwdRateCoefficients(rop_rates.data());
Eigen::SparseMatrix<double> jac = process_derivatives(m_reactantStoich, rop_rates);
Eigen::SparseMatrix<double> jac = calculateCompositionDerivatives(m_reactantStoich, rop_rates);

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

Eigen::SparseMatrix<double> GasKinetics::fwdRatesOfProgress_ddN()
Expand All @@ -488,7 +488,7 @@ Eigen::SparseMatrix<double> GasKinetics::fwdRatesOfProgress_ddN()
// forward reaction rate coefficients
vector_fp& rop_rates = m_rbuf0;
processFwdRateCoefficients(rop_rates.data());
return process_derivatives(m_reactantStoich, rop_rates, false);
return calculateCompositionDerivatives(m_reactantStoich, rop_rates, false);
}

Eigen::SparseMatrix<double> GasKinetics::revRatesOfProgress_ddN()
Expand All @@ -499,7 +499,7 @@ Eigen::SparseMatrix<double> GasKinetics::revRatesOfProgress_ddN()
vector_fp& rop_rates = m_rbuf0;
processFwdRateCoefficients(rop_rates.data());
processEquilibriumConstants(rop_rates.data());
return process_derivatives(m_revProductStoich, rop_rates, false);
return calculateCompositionDerivatives(m_revProductStoich, rop_rates, false);
}

Eigen::SparseMatrix<double> GasKinetics::netRatesOfProgress_ddN()
Expand All @@ -509,12 +509,12 @@ Eigen::SparseMatrix<double> GasKinetics::netRatesOfProgress_ddN()
// forward reaction rate coefficients
vector_fp& rop_rates = m_rbuf0;
processFwdRateCoefficients(rop_rates.data());
Eigen::SparseMatrix<double> jac = process_derivatives(m_reactantStoich, rop_rates,
Eigen::SparseMatrix<double> jac = calculateCompositionDerivatives(m_reactantStoich, rop_rates,
false);

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

bool GasKinetics::addReaction(shared_ptr<Reaction> r, bool resize)
Expand Down
26 changes: 13 additions & 13 deletions src/kinetics/InterfaceKinetics.cpp
Expand Up @@ -29,12 +29,12 @@ InterfaceKinetics::~InterfaceKinetics()

void InterfaceKinetics::resizeReactions()
{
Kinetics::resizeReactions();

// resize buffer
m_rbuf0.resize(nReactions());
m_rbuf1.resize(nReactions());

Kinetics::resizeReactions();

for (auto& rates : m_interfaceRates) {
rates->resize(nTotalSpecies(), nReactions(), nPhases());
// @todo ensure that ReactionData are updated; calling rates->update
Expand Down Expand Up @@ -621,7 +621,7 @@ Eigen::SparseMatrix<double> InterfaceKinetics::fwdRatesOfProgress_ddN()
// forward reaction rate coefficients
vector_fp& rop_rates = m_rbuf0;
processFwdRateCoefficients(rop_rates.data());
return process_derivatives(m_reactantStoich, rop_rates);
return calculateCompositionDerivatives(m_reactantStoich, rop_rates);
}

Eigen::SparseMatrix<double> InterfaceKinetics::revRatesOfProgress_ddN()
Expand All @@ -632,7 +632,7 @@ Eigen::SparseMatrix<double> InterfaceKinetics::revRatesOfProgress_ddN()
vector_fp& rop_rates = m_rbuf0;
processFwdRateCoefficients(rop_rates.data());
processEquilibriumConstants(rop_rates.data());
return process_derivatives(m_revProductStoich, rop_rates);
return calculateCompositionDerivatives(m_revProductStoich, rop_rates);
}

Eigen::SparseMatrix<double> InterfaceKinetics::netRatesOfProgress_ddN()
Expand All @@ -642,18 +642,19 @@ Eigen::SparseMatrix<double> InterfaceKinetics::netRatesOfProgress_ddN()
// forward reaction rate coefficients
vector_fp& rop_rates = m_rbuf0;
processFwdRateCoefficients(rop_rates.data());
Eigen::SparseMatrix<double> jac = process_derivatives(m_reactantStoich, rop_rates);
Eigen::SparseMatrix<double> jac = calculateCompositionDerivatives(m_reactantStoich,
rop_rates);

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

void InterfaceKinetics::setDerivativeSettings(const AnyMap& settings)
{
bool force = settings.empty();
if (force || settings.hasKey("skip-cov-dep")) {
m_jac_skip_cov_dependance = settings.getBool("skip-cov-dep",
if (force || settings.hasKey("skip-coverage-dependence")) {
m_jac_skip_cov_dependence = settings.getBool("skip-coverage-dependence",
true);
}
if (force || settings.hasKey("skip-electrochem")) {
Expand All @@ -669,9 +670,8 @@ void InterfaceKinetics::processFwdRateCoefficients(double* ropf)
{
// evaluate rate constants and equilibrium constants at temperature and phi
// (electric potential)
_update_rates_T();
// get updated activities (rates updated below)
_update_rates_C();
updateROP();
// copy rate coefficients into ropf
copy(m_rfn.begin(), m_rfn.end(), ropf);

Expand All @@ -681,7 +681,7 @@ void InterfaceKinetics::processFwdRateCoefficients(double* ropf)
}
}

Eigen::SparseMatrix<double> InterfaceKinetics::process_derivatives(
Eigen::SparseMatrix<double> InterfaceKinetics::calculateCompositionDerivatives(
StoichManagerN& stoich, const vector_fp& in)
{
Eigen::SparseMatrix<double> out;
Expand All @@ -692,9 +692,9 @@ Eigen::SparseMatrix<double> InterfaceKinetics::process_derivatives(
return out;
}

void InterfaceKinetics::assertDerivativesValid(const std::string& name)
void InterfaceKinetics::assertDerivativesValid(const string& name)
{
if (!m_jac_skip_cov_dependance) {
if (!m_jac_skip_cov_dependence) {
throw NotImplementedError(name, "Coverage-dependent reactions not supported.");
}
else if (!m_jac_skip_electrochem) {
Expand Down
5 changes: 0 additions & 5 deletions src/kinetics/InterfaceRate.cpp
Expand Up @@ -48,7 +48,6 @@ bool InterfaceData::update(const ThermoPhase& phase, const Kinetics& kin)
}

double T = phase.temperature();
double P = phase.pressure();
bool changed = false;
const auto& surf = dynamic_cast<const SurfPhase&>(
kin.thermo(kin.reactionPhaseIndex()));
Expand All @@ -57,10 +56,6 @@ bool InterfaceData::update(const ThermoPhase& phase, const Kinetics& kin)
density = surf.siteDensity();
changed = true;
}
if (P != pressure) {
pressure = P;
changed = true;
}
if (T != temperature) {
ReactionData::update(T);
sqrtT = sqrt(T);
Expand Down
2 changes: 1 addition & 1 deletion src/zeroD/IdealGasConstPressureMoleReactor.cpp
Expand Up @@ -139,7 +139,7 @@ Eigen::SparseMatrix<double> IdealGasConstPressureMoleReactor::jacobian()
species_trips.emplace_back(it.row(), it.col(), it.value());
}
}
mapSurfJacobian(species_trips);
addSurfaceJacobian(species_trips);
dnk_dnj.resize(ssize, ssize);
dnk_dnj.setFromTriplets(species_trips.begin(), species_trips.end());
}
Expand Down
25 changes: 13 additions & 12 deletions src/zeroD/IdealGasMoleReactor.cpp
Expand Up @@ -156,7 +156,7 @@ void IdealGasMoleReactor::eval(double time, double* LHS, double* RHS)
Eigen::SparseMatrix<double> IdealGasMoleReactor::jacobian()
{
if (m_nv == 0) {
throw CanteraError("IdealGasConstPressureMoleReactor::jacobian",
throw CanteraError("IdealGasMoleReactor::jacobian",
"Reactor must be initialized first.");
}
// clear former jacobian elements
Expand All @@ -165,24 +165,25 @@ Eigen::SparseMatrix<double> IdealGasMoleReactor::jacobian()
Eigen::SparseMatrix<double> dnk_dnj = m_kin->netProductionRates_ddN();
// species size that accounts for surface species
size_t ssize = m_nv - m_sidx;
// map concentration derivatives from surfaces to same vector
// map concentration derivatives from the surface chemistry jacobian
// to the reactor jacobian
if (!m_surfaces.empty()) {
std::vector<Eigen::Triplet<double>> species_trips(dnk_dnj.nonZeros());
for (int k = 0; k < dnk_dnj.outerSize(); k++) {
for (Eigen::SparseMatrix<double>::InnerIterator it(dnk_dnj, k); it; ++it) {
species_trips.emplace_back(it.row(), it.col(), it.value());
}
}
mapSurfJacobian(species_trips);
addSurfaceJacobian(species_trips);
dnk_dnj.resize(ssize, ssize);
dnk_dnj.setFromTriplets(species_trips.begin(), species_trips.end());
}
// add species to species derivatives elements to the jacobian
// calculate ROP derivatives, excluding the terms where dnk/dnj is zero but
// molarVolume * wdot is not, as it reduces matrix sparsity and diminishes
// performance.
for (int k = 0; k < dnk_dnj.outerSize(); k++) {
for (Eigen::SparseMatrix<double>::InnerIterator it(dnk_dnj, k); it; ++it) {
// scale value by appropriate quantities.
m_jac_trips.emplace_back(static_cast<int>(it.row() + m_sidx),
static_cast<int>(it.col() + m_sidx), it.value());
}
Expand Down Expand Up @@ -223,10 +224,10 @@ Eigen::SparseMatrix<double> IdealGasMoleReactor::jacobian()
m_kin->getNetProductionRates(netProductionRates.data());
m_thermo->getPartialMolarCp(specificHeat.data());
// surface phases
size_t shift = m_nsp;
size_t offset = m_nsp;
for (auto S : m_surfaces) {
S->thermo()->getPartialMolarCp(specificHeat.data() + shift);
S->thermo()->getPartialMolarEnthalpies(internal_energy.data() + shift);
S->thermo()->getPartialMolarCp(specificHeat.data() + offset);
S->thermo()->getPartialMolarEnthalpies(internal_energy.data() + offset);
// get additional net production rates
auto curr_kin = S->kinetics();
vector_fp prod_rates(curr_kin->nTotalSpecies());
Expand All @@ -238,10 +239,10 @@ Eigen::SparseMatrix<double> IdealGasMoleReactor::jacobian()
}
}
// scale production rates by areas
for (size_t i = shift; i < shift + S->thermo()->nSpecies(); i++) {
for (size_t i = offset; i < offset + S->thermo()->nSpecies(); i++) {
netProductionRates[i] *= S->area();
}
shift += S->thermo()->nSpecies();
offset += S->thermo()->nSpecies();
}
// convert Cp to Cv for ideal gas as Cp - Cv = R
for (size_t i = 0; i < m_nsp; i++) {
Expand All @@ -250,16 +251,16 @@ Eigen::SparseMatrix<double> IdealGasMoleReactor::jacobian()
}
// scale net production rates by volume to get molar rate
double qdot = internal_energy.dot(netProductionRates);
// find denominator ahead of time
// find the sum of n_i and cp_i
double NCv = 0.0;
double* moles = yCurrent.data() + m_sidx;
for (size_t i = 0; i < ssize; i++) {
NCv += moles[i] * specificHeat[i];
}
// Make denominator beforehand
// make denominator beforehand
double denom = 1 / (NCv * NCv);
Eigen::VectorXd uk_dnkdnj_sums = dnk_dnj.transpose() * internal_energy;
// Add derivatives to jac by spanning columns
// add derivatives to jacobian
for (size_t j = 0; j < ssize; j++) {
m_jac_trips.emplace_back(0, static_cast<int>(j + m_sidx),
(specificHeat[j] * qdot - NCv * uk_dnkdnj_sums[j]) * denom);
Expand Down
4 changes: 2 additions & 2 deletions src/zeroD/MoleReactor.cpp
Expand Up @@ -84,7 +84,7 @@ void MoleReactor::evalSurfaces(double* LHS, double* RHS, double* sdot)
}
}

void MoleReactor::mapSurfJacobian(std::vector<Eigen::Triplet<double>> &strips)
void MoleReactor::addSurfaceJacobian(vector<Eigen::Triplet<double>> &triplets)
{
for (auto& S : m_surfaces) {
// sync states
Expand All @@ -103,7 +103,7 @@ void MoleReactor::mapSurfJacobian(std::vector<Eigen::Triplet<double>> &strips)
double scalar = (row < m_nsp) ? m_vol : S->area();
scalar /= (col < m_nsp) ? m_vol : S->area();
// add to appropriate row and col
strips.emplace_back(row, col, scalar * it.value());
triplets.emplace_back(row, col, scalar * it.value());
}
}
}
Expand Down

0 comments on commit 05055e1

Please sign in to comment.