BulkKinetics.cpp

// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.
 
#include "cantera/kinetics/BulkKinetics.h"
#include "cantera/kinetics/Reaction.h"
#include "cantera/thermo/ThermoPhase.h"
#include "cantera/numerics/eigen_dense.h"
 
namespace Cantera
{
 
BulkKinetics::BulkKinetics() {
    setDerivativeSettings(AnyMap()); // use default settings
}
 
bool BulkKinetics::addReaction(shared_ptr<Reaction> r, bool resize)
{
    shared_ptr<ReactionRate> rate = r->rate();
    if (rate) {
        rate->setContext(*r, *this);
    }
 
    bool added = Kinetics::addReaction(r, resize);
    if (!added) {
        // undeclared species, etc.
        return false;
    }
    double dn = 0.0;
    for (const auto& [name, stoich] : r->products) {
        dn += stoich;
    }
    for (const auto& [name, stoich] : r->reactants) {
        dn -= stoich;
    }
 
    m_dn.push_back(dn);
 
    string rtype = rate->subType();
    if (rtype == "") {
        rtype = rate->type();
    }
 
    // If necessary, add new MultiRate evaluator
    if (m_rateTypes.find(rtype) == m_rateTypes.end()) {
        m_rateTypes[rtype] = m_rateHandlers.size();
        m_rateHandlers.push_back(rate->newMultiRate());
        m_rateHandlers.back()->resize(*this);
    }
 
    // Set index of rate to number of reaction within kinetics
    rate->setRateIndex(nReactions() - 1);
 
    // Add reaction rate to evaluator
    size_t index = m_rateTypes[rtype];
    m_rateHandlers[index]->add(nReactions() - 1, *rate);
 
    // Add reaction to third-body evaluator
    if (r->thirdBody() != nullptr) {
        addThirdBody(r);
    }
 
    m_concm.push_back(NAN);
    return true;
}
 
void BulkKinetics::addThirdBody(shared_ptr<Reaction> r)
{
    map<size_t, double> efficiencies;
    for (const auto& [name, efficiency] : r->thirdBody()->efficiencies) {
        size_t k = kineticsSpeciesIndex(name, false);
        if (k != npos) {
            efficiencies[k] = efficiency;
        } else if (!m_skipUndeclaredThirdBodies) {
            throw CanteraError("BulkKinetics::addThirdBody", "Found third-body"
                " efficiency for undefined species '{}' while adding reaction '{}'",
                name, r->equation());
        } else {
            m_hasUndeclaredThirdBodies = true;
        }
    }
    m_multi_concm.install(nReactions() - 1, efficiencies,
                          r->thirdBody()->default_efficiency,
                          r->thirdBody()->mass_action);
}
 
void BulkKinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
{
    // operations common to all reaction types
    Kinetics::modifyReaction(i, rNew);
 
    shared_ptr<ReactionRate> rate = rNew->rate();
    string rtype = rate->subType();
    if (rtype == "") {
        rtype = rate->type();
    }
 
    // Ensure that MultiRate evaluator is available
    if (m_rateTypes.find(rtype) == m_rateTypes.end()) {
        throw CanteraError("BulkKinetics::modifyReaction",
                "Evaluator not available for type '{}'.", rtype);
    }
 
    // Replace reaction rate to evaluator
    size_t index = m_rateTypes[rtype];
    rate->setRateIndex(i);
    rate->setContext(*rNew, *this);
    m_rateHandlers[index]->replace(i, *rate);
    invalidateCache();
}
 
void BulkKinetics::resizeSpecies()
{
    Kinetics::resizeSpecies();
    m_act_conc.resize(m_kk);
    m_phys_conc.resize(m_kk);
    m_grt.resize(m_kk);
    for (auto& rates : m_rateHandlers) {
        rates->resize(*this);
    }
}
 
void BulkKinetics::resizeReactions()
{
    Kinetics::resizeReactions();
    m_rbuf0.resize(nReactions());
    m_rbuf1.resize(nReactions());
    m_rbuf2.resize(nReactions());
    m_rbuf3.assign(nTotalSpecies(), 0.0);
    m_kf0.resize(nReactions());
    m_sbuf0.resize(nTotalSpecies());
    m_state.resize(thermo().stateSize());
    m_multi_concm.resizeCoeffs(nTotalSpecies(), nReactions());
    for (auto& rates : m_rateHandlers) {
        rates->resize(*this);
        // @todo ensure that ReactionData are updated; calling rates->update
        //      blocks correct behavior in update_rates_T
        //      and running updateROP() is premature
    }
}
 
void BulkKinetics::setMultiplier(size_t i, double f)
{
    Kinetics::setMultiplier(i, f);
    m_ROP_ok = false;
}
 
void BulkKinetics::invalidateCache()
{
    Kinetics::invalidateCache();
    m_ROP_ok = false;
}
 
void BulkKinetics::getFwdRateConstants(span<double> kfwd)
{
    checkArraySize("BulkKinetics::getFwdRateConstants", kfwd.size(), nReactions());
    updateROP();
    copy(m_rfn.begin(), m_rfn.end(), kfwd.begin());
    if (legacy_rate_constants_used()) {
        processThirdBodies(kfwd);
    }
}
 
void BulkKinetics::getEquilibriumConstants(span<double> kc)
{
    checkArraySize("BulkKinetics::getEquilibriumConstants", kc.size(), nReactions());
    updateROP();
 
    vector<double>& delta_gibbs0 = m_rbuf0;
    fill(delta_gibbs0.begin(), delta_gibbs0.end(), 0.0);
 
    // compute Delta G^0 for all reactions
    getReactionDelta(m_grt, delta_gibbs0);
 
    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] * logStandConc);
    }
}
 
void BulkKinetics::getRevRateConstants(span<double> krev, bool doIrreversible)
{
    // go get the forward rate constants. -> note, we don't really care about
    // speed or redundancy in these informational routines.
    getFwdRateConstants(krev);
 
    if (doIrreversible) {
        getEquilibriumConstants(m_rbuf0);
        for (size_t i = 0; i < nReactions(); i++) {
            krev[i] /= m_rbuf0[i];
        }
    } else {
        // m_rkcn[] is zero for irreversible reactions
        for (size_t i = 0; i < nReactions(); i++) {
            krev[i] *= m_rkcn[i];
        }
    }
}
 
void BulkKinetics::getDeltaGibbs(span<double> deltaG)
{
    // Get the chemical potentials for each species
    thermo().getChemPotentials(m_sbuf0);
    // Use the stoichiometric manager to find deltaG for each reaction.
    getReactionDelta(m_sbuf0, deltaG);
}
 
void BulkKinetics::getDeltaEnthalpy(span<double> deltaH)
{
    // Get the partial molar enthalpy for each species
    thermo().getPartialMolarEnthalpies(m_sbuf0);
    // Use the stoichiometric manager to find deltaH for each reaction.
    getReactionDelta(m_sbuf0, deltaH);
}
 
void BulkKinetics::getDeltaEntropy(span<double> deltaS)
{
    // Get the partial molar entropy for each species
    thermo().getPartialMolarEntropies(m_sbuf0);
    // Use the stoichiometric manager to find deltaS for each reaction.
    getReactionDelta(m_sbuf0, deltaS);
}
 
void BulkKinetics::getDeltaSSGibbs(span<double> deltaG)
{
    // Get the standard state chemical potentials of the species. This is the
    // array of chemical potentials at unit activity. We define these here as
    // the chemical potentials of the pure species at the temperature and
    // pressure of the solution.
    thermo().getStandardChemPotentials(m_sbuf0);
    // Use the stoichiometric manager to find deltaG for each reaction.
    getReactionDelta(m_sbuf0, deltaG);
}
 
void BulkKinetics::getDeltaSSEnthalpy(span<double> deltaH)
{
    // Get the standard state enthalpies of the species.
    thermo().getEnthalpy_RT(m_sbuf0);
    for (size_t k = 0; k < m_kk; k++) {
        m_sbuf0[k] *= thermo().RT();
    }
    // Use the stoichiometric manager to find deltaH for each reaction.
    getReactionDelta(m_sbuf0, deltaH);
}
 
void BulkKinetics::getDeltaSSEntropy(span<double> deltaS)
{
    // Get the standard state entropy of the species. We define these here as
    // the entropies of the pure species at the temperature and pressure of the
    // solution.
    thermo().getEntropy_R(m_sbuf0);
    for (size_t k = 0; k < m_kk; k++) {
        m_sbuf0[k] *= GasConstant;
    }
    // Use the stoichiometric manager to find deltaS for each reaction.
    getReactionDelta(m_sbuf0, deltaS);
}
 
void BulkKinetics::getDerivativeSettings(AnyMap& settings) const
{
    settings["skip-third-bodies"] = m_jac_skip_third_bodies;
    settings["skip-falloff"] = m_jac_skip_falloff;
    settings["skip-nonideal"] = m_jac_skip_nonideal;
    settings["rtol-delta"] = m_jac_rtol_delta;
}
 
void BulkKinetics::setDerivativeSettings(const AnyMap& settings)
{
    bool force = settings.empty();
    if (force || settings.hasKey("skip-third-bodies")) {
        m_jac_skip_third_bodies = settings.getBool("skip-third-bodies", false);
    }
    if (force || settings.hasKey("skip-falloff")) {
        m_jac_skip_falloff = settings.getBool("skip-falloff", false);
    }
    if (force || settings.hasKey("skip-nonideal")) {
        m_jac_skip_nonideal = settings.getBool("skip-nonideal", false);
    }
    if (force || settings.hasKey("rtol-delta")) {
        m_jac_rtol_delta = settings.getDouble("rtol-delta", 1e-8);
    }
}
 
void BulkKinetics::getFwdRateConstants_ddT(span<double> dkfwd)
{
    assertDerivativesValid("BulkKinetics::getFwdRateConstants_ddT");
    updateROP();
    process_ddT(m_rfn, dkfwd);
}
 
void BulkKinetics::getFwdRatesOfProgress_ddT(span<double> drop)
{
    assertDerivativesValid("BulkKinetics::getFwdRatesOfProgress_ddT");
    updateROP();
    process_ddT(m_ropf, drop);
}
 
void BulkKinetics::getRevRatesOfProgress_ddT(span<double> drop)
{
    assertDerivativesValid("BulkKinetics::getRevRatesOfProgress_ddT");
    updateROP();
    process_ddT(m_ropr, drop);
 
    // reverse rop times scaled inverse equilibrium constant derivatives
    copy(m_ropr.begin(), m_ropr.end(), m_rbuf2.begin());
    applyEquilibriumConstants_ddT(m_rbuf2);
    Eigen::Map<Eigen::VectorXd> dRevRop(drop.data(), nReactions());
    Eigen::Map<Eigen::VectorXd> dRevRop2(m_rbuf2.data(), nReactions());
    dRevRop += dRevRop2;
}
 
void BulkKinetics::getNetRatesOfProgress_ddT(span<double> drop)
{
    assertDerivativesValid("BulkKinetics::getNetRatesOfProgress_ddT");
    updateROP();
    process_ddT(m_ropnet, drop);
 
    // reverse rop times scaled inverse equilibrium constant derivatives
    copy(m_ropr.begin(), m_ropr.end(), m_rbuf2.begin());
    applyEquilibriumConstants_ddT(m_rbuf2);
    Eigen::Map<Eigen::VectorXd> dNetRop(drop.data(), nReactions());
    Eigen::Map<Eigen::VectorXd> dNetRop2(m_rbuf2.data(), nReactions());
    dNetRop -= dNetRop2;
}
 
void BulkKinetics::getFwdRateConstants_ddP(span<double> dkfwd)
{
    assertDerivativesValid("BulkKinetics::getFwdRateConstants_ddP");
    updateROP();
    process_ddP(m_rfn, dkfwd);
}
 
void BulkKinetics::getFwdRatesOfProgress_ddP(span<double> drop)
{
    assertDerivativesValid("BulkKinetics::getFwdRatesOfProgress_ddP");
    updateROP();
    process_ddP(m_ropf, drop);
}
 
void BulkKinetics::getRevRatesOfProgress_ddP(span<double> drop)
{
    assertDerivativesValid("BulkKinetics::getRevRatesOfProgress_ddP");
    updateROP();
    process_ddP(m_ropr, drop);
}
 
void BulkKinetics::getNetRatesOfProgress_ddP(span<double> drop)
{
    assertDerivativesValid("BulkKinetics::getNetRatesOfProgress_ddP");
    updateROP();
    process_ddP(m_ropnet, drop);
}
 
void BulkKinetics::getFwdRateConstants_ddC(span<double> dkfwd)
{
    assertDerivativesValid("BulkKinetics::getFwdRateConstants_ddC");
    updateROP();
    process_ddC(m_reactantStoich, m_rfn, dkfwd, false);
}
 
void BulkKinetics::getFwdRatesOfProgress_ddC(span<double> drop)
{
    assertDerivativesValid("BulkKinetics::getFwdRatesOfProgress_ddC");
    updateROP();
    process_ddC(m_reactantStoich, m_ropf, drop);
}
 
void BulkKinetics::getRevRatesOfProgress_ddC(span<double> drop)
{
    assertDerivativesValid("BulkKinetics::getRevRatesOfProgress_ddC");
    updateROP();
    return process_ddC(m_revProductStoich, m_ropr, drop);
}
 
void BulkKinetics::getNetRatesOfProgress_ddC(span<double> drop)
{
    assertDerivativesValid("BulkKinetics::getNetRatesOfProgress_ddC");
    updateROP();
    process_ddC(m_reactantStoich, m_ropf, drop);
    Eigen::Map<Eigen::VectorXd> dNetRop(drop.data(), nReactions());
 
    process_ddC(m_revProductStoich, m_ropr, m_rbuf2);
    Eigen::Map<Eigen::VectorXd> dNetRop2(m_rbuf2.data(), nReactions());
    dNetRop -= dNetRop2;
}
 
Eigen::SparseMatrix<double> BulkKinetics::fwdRatesOfProgress_ddX()
{
    assertDerivativesValid("BulkKinetics::fwdRatesOfProgress_ddX");
 
    // forward reaction rate coefficients
    vector<double>& rop_rates = m_rbuf0;
    getFwdRateConstants(rop_rates);
    return calculateCompositionDerivatives(m_reactantStoich, rop_rates);
}
 
Eigen::SparseMatrix<double> BulkKinetics::revRatesOfProgress_ddX()
{
    assertDerivativesValid("BulkKinetics::revRatesOfProgress_ddX");
 
    // reverse reaction rate coefficients
    vector<double>& rop_rates = m_rbuf0;
    getFwdRateConstants(rop_rates);
    applyEquilibriumConstants(rop_rates);
    return calculateCompositionDerivatives(m_revProductStoich, rop_rates);
}
 
Eigen::SparseMatrix<double> BulkKinetics::netRatesOfProgress_ddX()
{
    assertDerivativesValid("BulkKinetics::netRatesOfProgress_ddX");
 
    // forward reaction rate coefficients
    vector<double>& rop_rates = m_rbuf0;
    getFwdRateConstants(rop_rates);
    auto jac = calculateCompositionDerivatives(m_reactantStoich, rop_rates);
 
    // reverse reaction rate coefficients
    applyEquilibriumConstants(rop_rates);
    return jac - calculateCompositionDerivatives(m_revProductStoich, rop_rates);
}
 
Eigen::SparseMatrix<double> BulkKinetics::fwdRatesOfProgress_ddCi()
{
    assertDerivativesValid("BulkKinetics::fwdRatesOfProgress_ddCi");
 
    // forward reaction rate coefficients
    vector<double>& rop_rates = m_rbuf0;
    getFwdRateConstants(rop_rates);
    return calculateCompositionDerivatives(m_reactantStoich, rop_rates, false);
}
 
Eigen::SparseMatrix<double> BulkKinetics::revRatesOfProgress_ddCi()
{
    assertDerivativesValid("BulkKinetics::revRatesOfProgress_ddCi");
 
    // reverse reaction rate coefficients
    vector<double>& rop_rates = m_rbuf0;
    getFwdRateConstants(rop_rates);
    applyEquilibriumConstants(rop_rates);
    return calculateCompositionDerivatives(m_revProductStoich, rop_rates, false);
}
 
Eigen::SparseMatrix<double> BulkKinetics::netRatesOfProgress_ddCi()
{
    assertDerivativesValid("BulkKinetics::netRatesOfProgress_ddCi");
 
    // forward reaction rate coefficients
    vector<double>& rop_rates = m_rbuf0;
    getFwdRateConstants(rop_rates);
    auto jac = calculateCompositionDerivatives(m_reactantStoich, rop_rates, false);
 
    // reverse reaction rate coefficients
    applyEquilibriumConstants(rop_rates);
    return jac - calculateCompositionDerivatives(m_revProductStoich, rop_rates, false);
}
 
void BulkKinetics::updateROP()
{
    static const int cacheId = m_cache.getId();
    CachedScalar last = m_cache.getScalar(cacheId);
    double T = thermo().temperature();
    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);
        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, m_delta_gibbs0);
 
        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);
        }
 
        for (size_t i = 0; i != m_irrev.size(); ++i) {
            m_rkcn[ m_irrev[i] ] = 0.0;
        }
    }
 
    if (!last.validate(T, rho, statenum)) {
        // Update terms dependent on species concentrations and temperature
        thermo().getActivityConcentrations(m_act_conc);
        thermo().getConcentrations(m_phys_conc);
        double ctot = thermo().molarDensity();
 
        // Third-body objects interacting with MultiRate evaluator
        m_multi_concm.update(m_phys_conc, ctot, m_concm);
        m_ROP_ok = false;
    }
 
    // loop over MultiRate evaluators for each reaction type
    for (auto& rates : m_rateHandlers) {
        bool changed = rates->update(thermo(), *this);
        if (changed) {
            rates->getRateConstants(m_kf0);
            m_ROP_ok = false;
        }
    }
 
    if (m_ROP_ok) {
        // rates of progress are up-to-date only if both the thermodynamic state
        // and m_perturb are unchanged
        return;
    }
 
    // 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];
    }
 
    copy(m_rfn.begin(), m_rfn.end(), m_ropf.begin());
    processThirdBodies(m_ropf);
    copy(m_ropf.begin(), m_ropf.end(), m_ropr.begin());
 
    // for reversible reactions, multiply ropr by concentration products
    applyEquilibriumConstants(m_ropr);
 
    for (auto& rates : m_rateHandlers) {
        rates->modifyRateConstants(m_ropf, m_ropr);
    }
 
    // multiply ropf and ropr by concentration products
    m_reactantStoich.multiply(m_act_conc, m_ropf);
    m_revProductStoich.multiply(m_act_conc, m_ropr);
 
    for (size_t j = 0; j != nReactions(); ++j) {
        m_ropnet[j] = m_ropf[j] - m_ropr[j];
    }
 
    for (size_t i = 0; i < m_rfn.size(); i++) {
        AssertFinite(m_rfn[i], "BulkKinetics::updateROP",
                     "m_rfn[{}] is not finite.", i);
        AssertFinite(m_ropf[i], "BulkKinetics::updateROP",
                     "m_ropf[{}] is not finite.", i);
        AssertFinite(m_ropr[i], "BulkKinetics::updateROP",
                     "m_ropr[{}] is not finite.", i);
    }
    m_ROP_ok = true;
}
 
void BulkKinetics::getThirdBodyConcentrations(span<double> concm)
{
    checkArraySize("BulkKinetics::getThirdBodyConcentrations", concm.size(), nReactions());
    updateROP();
    std::copy(m_concm.begin(), m_concm.end(), concm.begin());
}
 
void BulkKinetics::processThirdBodies(span<double> rop)
{
    // reactions involving third body
    if (!m_concm.empty()) {
        m_multi_concm.multiply(rop, m_concm);
    }
}
 
void BulkKinetics::applyEquilibriumConstants(span<double> rop)
{
    // For reverse rates computed from thermochemistry, multiply the forward
    // rate coefficients by the reciprocals of the equilibrium constants
    for (size_t i = 0; i < nReactions(); ++i) {
        rop[i] *= m_rkcn[i];
    }
}
 
void BulkKinetics::applyEquilibriumConstants_ddT(span<double> drkcn)
{
    double T = thermo().temperature();
    double P = thermo().pressure();
    double rrt = 1. / thermo().RT();
 
    vector<double>& grt = m_sbuf0;
    vector<double>& delta_gibbs0 = m_rbuf1;
    fill(delta_gibbs0.begin(), delta_gibbs0.end(), 0.0);
 
    // compute perturbed Delta G^0 for all reversible reactions
    thermo().saveState(m_state);
    thermo().setState_TP(T * (1. + m_jac_rtol_delta), P);
    thermo().getStandardChemPotentials(grt);
    getRevReactionDelta(grt, delta_gibbs0);
 
    // apply scaling for derivative of inverse equilibrium constant
    double Tinv = 1. / T;
    double rrt_dTinv = rrt * Tinv / m_jac_rtol_delta;
    double rrtt = rrt * Tinv;
    for (size_t i = 0; i < m_revindex.size(); i++) {
        size_t irxn = m_revindex[i];
        double factor = delta_gibbs0[irxn] - m_delta_gibbs0[irxn];
        factor *= rrt_dTinv;
        factor += m_dn[irxn] * Tinv - m_delta_gibbs0[irxn] * rrtt;
        drkcn[irxn] *= factor;
    }
 
    for (size_t i = 0; i < m_irrev.size(); ++i) {
        drkcn[m_irrev[i]] = 0.0;
    }
 
    thermo().restoreState(m_state);
}
 
void BulkKinetics::process_ddT(span<const double> in, span<double> drop)
{
    // apply temperature derivative. The reference rate constant must be the
    // multiplier-free m_kf0: rate handlers that evaluate the derivative by an internal
    // finite difference compare against a perturbed constant that omits the reaction
    // multiplier
    copy(in.begin(), in.end(), drop.begin());
    for (auto& rates : m_rateHandlers) {
        rates->processRateConstants_ddT(drop, m_kf0, m_jac_rtol_delta);
    }
}
 
void BulkKinetics::process_ddP(span<const double> in, span<double> drop)
{
    // apply pressure derivative
    copy(in.begin(), in.end(), drop.begin());
    for (auto& rates : m_rateHandlers) {
        rates->processRateConstants_ddP(drop, m_kf0, m_jac_rtol_delta);
    }
}
 
void BulkKinetics::process_ddC(StoichManagerN& stoich, span<const double> in,
                               span<double> drop, bool mass_action)
{
    Eigen::Map<Eigen::VectorXd> out(drop.data(), nReactions());
    out.setZero();
    double ctot_inv = 1. / thermo().molarDensity();
 
    // derivatives due to concentrations in law of mass action
    if (mass_action) {
        stoich.scale(in, span<double>(out.data(), nReactions()), ctot_inv);
    }
    if (m_jac_skip_third_bodies || m_multi_concm.empty()) {
        return;
    }
 
    // derivatives due to third-body colliders in law of mass action
    Eigen::Map<Eigen::VectorXd> outM(m_rbuf1.data(), nReactions());
    if (mass_action) {
        outM.fill(0.);
        m_multi_concm.scale(in, asSpan(outM), ctot_inv);
        out += outM;
    }
 
    // derivatives due to reaction rates depending on third-body colliders
    if (!m_jac_skip_falloff) {
        m_multi_concm.scaleM(in, asSpan(outM), m_concm, ctot_inv);
        for (auto& rates : m_rateHandlers) {
            // processing step assigns zeros to entries not dependent on M; use
            // the multiplier-free m_kf0 reference (see process_ddT)
            rates->processRateConstants_ddM(asSpan(outM), m_kf0, m_jac_rtol_delta);
        }
        out += outM;
    }
}
 
Eigen::SparseMatrix<double> BulkKinetics::calculateCompositionDerivatives(
    StoichManagerN& stoich, span<const double> in, bool ddX)
{
    Eigen::SparseMatrix<double> out;
    vector<double>& scaled = m_rbuf1;
    vector<double>& outV = m_rbuf2;
 
    // convert from concentration to mole fraction output
    copy(in.begin(), in.end(), scaled.begin());
    if (ddX) {
        double ctot = thermo().molarDensity();
        for (size_t i = 0; i < nReactions(); ++i) {
            scaled[i] *= ctot;
        }
    }
 
    // derivatives handled by StoichManagerN
    copy(scaled.begin(), scaled.end(), outV.begin());
    processThirdBodies(outV);
    out = stoich.derivatives(m_act_conc, outV);
    if (m_jac_skip_third_bodies || m_multi_concm.empty()) {
        return out;
    }
 
    // derivatives due to law of mass action
    copy(scaled.begin(), scaled.end(), outV.begin());
    stoich.multiply(m_act_conc, outV);
 
    // derivatives due to reaction rates depending on third-body colliders
    if (!m_jac_skip_falloff) {
        for (auto& rates : m_rateHandlers) {
            // processing step does not modify entries not dependent on M
            rates->processRateConstants_ddM(outV, m_kf0, m_jac_rtol_delta, false);
        }
    }
 
    // derivatives handled by ThirdBodyCalc
    out += m_multi_concm.derivatives(outV);
    return out;
}
 
namespace {
 
//! Accumulate `sign * stoichMat * dRoP` into `out`'s compressed column storage,
//! where stoichMat is K×R (species × reactions) and dRoP is R×K.
//! Each column m of the product is computed via a length-K dense workspace
//! (`scatter`) that is zeroed after each column is folded into `out`.
//! `out` must already have a nonzero slot for every (k, m) that the product
//! touches. The template parameter allows dRoP to be an Eigen::Map (returned
//! by StoichManagerN::derivatives) as well as a concrete SparseMatrix.
template <class Derived>
void accumulateStoichProduct(const Eigen::SparseMatrix<double>& stoichMat,
                              const Eigen::SparseMatrixBase<Derived>& dRoP,
                              Eigen::SparseMatrix<double>& out,
                              double sign,
                              vector<double>& scatter)
{
    const int* outOuter = out.outerIndexPtr();
    const int* outInner = out.innerIndexPtr();
    double* outVal = out.valuePtr();
    for (int m = 0; m < dRoP.outerSize(); m++) {
        // accumulate column m of the product into the dense scatter, touching
        // only species rows k that actually appear
        for (typename Derived::InnerIterator it(dRoP.derived(), m); it; ++it) {
            int i = it.row(); // reaction index
            double v = sign * it.value();
            for (Eigen::SparseMatrix<double>::InnerIterator sit(stoichMat, i);
                 sit; ++sit) {
                scatter[sit.row()] += sit.value() * v;
            }
        }
        // fold the scatter into out's compressed storage for column m and reset
        for (int slot = outOuter[m]; slot < outOuter[m + 1]; slot++) {
            int k = outInner[slot];
            outVal[slot] += scatter[k];
            scatter[k] = 0.0;
        }
#ifndef NDEBUG
        // The fold loop above resets scatter[k] only for rows k that have a slot
        // in out's column m. Any leftover nonzero therefore flags a product entry
        // that out's sparsity pattern cannot store -- it would be silently dropped
        // (and leak into the next column). Catch this in debug builds; reset so
        // the error message is the only consequence.
        for (size_t k = 0; k < scatter.size(); k++) {
            if (scatter[k] != 0.0) {
                int row = static_cast<int>(k);
                scatter[k] = 0.0;
                throw CanteraError("accumulateStoichProduct",
                    "Supplied Jacobian sparsity pattern is missing a required "
                    "nonzero slot at (row={}, col={}); the pattern is incompatible "
                    "with the current mechanism and derivative settings.", row, m);
            }
        }
#endif
    }
}
 
} // anonymous namespace
 
void BulkKinetics::accumulateCompositionDerivatives(
    StoichManagerN& stoich, span<const double> in,
    Eigen::SparseMatrix<double>& out, double sign)
{
    vector<double>& outV = m_rbuf2;
 
    // derivatives due to law of mass action
    copy(in.begin(), in.end(), outV.begin());
    processThirdBodies(outV);
    // stoich.derivatives() returns an Eigen::Map over internal storage that is
    // overwritten on the next call, so pass directly to accumulateStoichProduct
    accumulateStoichProduct(m_stoichMatrix, stoich.derivatives(m_act_conc, outV),
                      out, sign, m_rbuf3);
 
    if (m_jac_skip_third_bodies || m_multi_concm.empty()) {
        return;
    }
 
    // derivatives due to reaction rates depending on third-body colliders
    copy(in.begin(), in.end(), outV.begin());
    stoich.multiply(m_act_conc, outV);
    if (!m_jac_skip_falloff) {
        for (auto& rates : m_rateHandlers) {
            // processing step does not modify entries not dependent on M; use
            // the multiplier-free m_kf0 reference.
            rates->processRateConstants_ddM(outV, m_kf0, m_jac_rtol_delta, false);
        }
    }
    // derivatives handled by ThirdBodyCalc
    accumulateStoichProduct(m_stoichMatrix, m_multi_concm.derivatives(outV),
                      out, sign, m_rbuf3);
}
 
void BulkKinetics::netProductionRates_ddCi(Eigen::SparseMatrix<double>& jac)
{
    assertDerivativesValid("BulkKinetics::netProductionRates_ddCi");
    if (jac.nonZeros() == 0) {
        // build the fixed pattern once; values overwritten below
        jac = Kinetics::netProductionRates_ddCi();
        jac.makeCompressed();
    }
    std::fill(jac.valuePtr(), jac.valuePtr() + jac.nonZeros(), 0.0);
 
    // forward reaction rate coefficients
    getFwdRateConstants(m_rbuf0);
    accumulateCompositionDerivatives(m_reactantStoich, m_rbuf0, jac, 1.0);
 
    // reverse reaction rate coefficients
    applyEquilibriumConstants(m_rbuf0);
    accumulateCompositionDerivatives(m_revProductStoich, m_rbuf0, jac, -1.0);
}
 
void BulkKinetics::assertDerivativesValid(const string& name)
{
    if (!m_jac_skip_nonideal && !thermo().isIdeal()) {
        throw NotImplementedError(name,
            "Not supported for non-ideal ThermoPhase models.");
    }
}
 
}