Skip to content

Commit

Permalink
This commit removes corrections made to IdealGasConstPressureMoleReactor
Browse files Browse the repository at this point in the history 
The corrections cause the matrices to become much more dense which
greatly dimishes improvements due to sparsity. It also adds routines
to scale the derivatives by appropriate terms based on the phase.
  • Loading branch information
@anthony-walker @speth
anthony-walker authored and speth committed Apr 11, 2023
1 parent b7ee632 commit 199ee2a
Show file tree
Hide file tree
Showing 2 changed files with 20 additions and 89 deletions.
3 changes: 2 additions & 1 deletion src/kinetics/GasKinetics.cpp
View file
Expand Up @@ -509,7 +509,8 @@ 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, false);
Eigen::SparseMatrix<double> jac = process_derivatives(m_reactantStoich, rop_rates,
false);

// reverse reaction rate coefficients
processEquilibriumConstants(rop_rates.data());
Expand Down
106 changes: 18 additions & 88 deletions src/zeroD/IdealGasConstPressureMoleReactor.cpp
View file
Expand Up @@ -128,20 +128,20 @@ Eigen::SparseMatrix<double> IdealGasConstPressureMoleReactor::jacobian()
// clear former jacobian elements
m_jac_trips.clear();
// determine species derivatives w.r.t concentration
Eigen::SparseMatrix<double> dwdot_dci = m_kin->netProductionRates_ddN();
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
if (!m_surfaces.empty()) {
std::vector<Eigen::Triplet<double>> species_trips(dwdot_dci.nonZeros());
for (int k = 0; k < dwdot_dci.outerSize(); k++) {
for (Eigen::SparseMatrix<double>::InnerIterator it(dwdot_dci, k); it; ++it) {
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());
}
}
mapSurfConcJacobian(species_trips);
dwdot_dci.resize(ssize, ssize);
dwdot_dci.setFromTriplets(species_trips.begin(), species_trips.end());
mapSurfJacobian(species_trips);
dnk_dnj.resize(ssize, ssize);
dnk_dnj.setFromTriplets(species_trips.begin(), species_trips.end());
}
// get net production rates
Eigen::VectorXd netProductionRates = Eigen::VectorXd::Zero(ssize);
Expand All @@ -159,89 +159,15 @@ Eigen::SparseMatrix<double> IdealGasConstPressureMoleReactor::jacobian()
}
}
}
// get current state
vector_fp yCurrent(m_nv);
getState(yCurrent.data());
// get total moles
double N = std::accumulate(yCurrent.data() + m_sidx,
yCurrent.data() + m_sidx + m_nsp, 0.0);
// get d[n_i]/dni
Eigen::SparseMatrix<double> dci_dnj(ssize, ssize);
std::vector<Eigen::Triplet<double>> dci_dnj_trips(ssize);
// gas phase contributions to dci_dni
double dci;
double v_Inv = 1 / m_vol;
double nv_Inv = 1 / (m_vol * N);
double* moles = yCurrent.data() + m_sidx; // kmol per s
for (size_t i = 0; i < m_nsp; i++) {
for (size_t j = 0; j < m_nsp; j++) {
// calculate
if (i == j) {
dci = v_Inv - moles[i] * nv_Inv;
} else {
dci = - moles[i] * nv_Inv;
}
// check greater than zero
if (dci > 0) {
dci_dnj_trips.emplace_back(i, j, dci);
}
}
}
// surface contributions
size_t shift = m_nsp;
for (auto S : m_surfaces) {
for (size_t i = shift; i < S->thermo()->nSpecies() + shift; i++) {
dci_dnj_trips.emplace_back(i, i, 1 / S->area());
}
shift += S->thermo()->nSpecies();
}
// create sparse matrix
dci_dnj.setFromTriplets(dci_dnj_trips.begin(), dci_dnj_trips.end());
// derivatives dwdot / dnj
Eigen::SparseMatrix<double> dnk_dnj = dwdot_dci * dci_dnj;
// get dwdot/dT dT / dnj
Eigen::VectorXd dwdot_dnj_vec(ssize);
// gas phase dwdot / dT
m_kin->getNetProductionRates_ddT(dwdot_dnj_vec.data());
// dT / d[n_i] * d[n_i] / d n_i
double dT_dci = - m_thermo->pressure() /
(GasConstant * m_thermo->molarDensity() * m_thermo->molarDensity());
Eigen::VectorXd dTdnj = dci_dnj.transpose() *
Eigen::VectorXd::Constant(ssize, dT_dci);
// surface values
shift = m_nsp;
for (auto S : m_surfaces) {
for (size_t i = shift; i < S->thermo()->nSpecies() + shift; i++) {
dTdnj[i] = 0;
}
shift += S->thermo()->nSpecies();
}
// get total derivatives
dnk_dnj = dnk_dnj + (dwdot_dnj_vec * dTdnj.transpose());
// convert to dndotk / dnj
double molarVolume = m_thermo->molarVolume();
double molarVol = m_thermo->molarVolume();
// 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) {
// gas phase w.r.t any phase
// gas phase species need the addition of V / N * omega_dot
if (it.row() < m_nsp) {
it.valueRef() = m_vol * it.value() + netProductionRates[it.row()] * molarVolume;
// surf phase w.r.t any phase
} else {
double A = 0.0;
shift = m_nsp;
for (auto &S: m_surfaces) {
size_t nspec = S->thermo()->nSpecies();
if (shift < it.row() && it.row() < nspec + shift) {
A = S->area();
break;
}
shift += nspec;
}
// get appropriate surface area
it.valueRef() = A * it.value();
it.valueRef() = it.value() + netProductionRates[it.row()] * molarVol;
}
m_jac_trips.emplace_back(static_cast<int>(it.row() + m_sidx),
static_cast<int>(it.col() + m_sidx), it.value());
Expand All @@ -252,6 +178,9 @@ Eigen::SparseMatrix<double> IdealGasConstPressureMoleReactor::jacobian()
// getting perturbed state for finite difference
double deltaTemp = m_thermo->temperature()
* std::sqrt(std::numeric_limits<double>::epsilon());
// get current state
vector_fp yCurrent(m_nv);
getState(yCurrent.data());
// finite difference temperature derivatives
vector_fp lhsPerturbed(m_nv, 1.0), lhsCurrent(m_nv, 1.0);
vector_fp rhsPerturbed(m_nv, 0.0), rhsCurrent(m_nv, 0.0);
Expand Down Expand Up @@ -284,7 +213,7 @@ Eigen::SparseMatrix<double> IdealGasConstPressureMoleReactor::jacobian()
netProductionRates[i] *= m_vol;
}
// surface phases
shift = m_nsp;
size_t shift = m_nsp;
for (auto S : m_surfaces) {
S->thermo()->getPartialMolarCp(specificHeat.data() + shift);
S->thermo()->getPartialMolarEnthalpies(enthalpy.data() + shift);
Expand All @@ -297,15 +226,16 @@ Eigen::SparseMatrix<double> IdealGasConstPressureMoleReactor::jacobian()
double qdot = enthalpy.dot(netProductionRates);
// find denominator ahead of time
double NCp = 0.0;
double* moles = yCurrent.data() + m_sidx;
for (size_t i = 0; i < ssize; i++) {
NCp += yCurrent[i + m_sidx] * specificHeat[i];
NCp += moles[i] * specificHeat[i];
}
double denom = 1 / (NCp * NCp);
Eigen::VectorXd hk_dnkdnj_sum = dnk_dnj.transpose() * enthalpy;
Eigen::VectorXd hk_dnkdnj_sums = dnk_dnj.transpose() * enthalpy;
// Add derivatives to jac by spanning columns
for (size_t j = 0; j < ssize; j++) {
m_jac_trips.emplace_back(0, static_cast<int>(j + m_sidx),
(specificHeat[j] * qdot - NCp * hk_dnkdnj_sum[j]) * denom);
(specificHeat[j] * qdot - NCp * hk_dnkdnj_sums[j]) * denom);
}
}
// convert triplets to sparse matrix
Expand Down

0 comments on commit 199ee2a

Please sign in to comment.