Reverse reaction (super-elastic) electron collision reaction

include/cantera/kinetics/BulkKinetics.h

@@ -9,7 +9,6 @@
 #define CT_BULKKINETICS_H
 
 #include "Kinetics.h"
-#include "MultiRate.h"
 #include "ThirdBodyCalc.h"
 
 namespace Cantera
@@ -28,7 +27,6 @@ class BulkKinetics : public Kinetics
         return "bulk";
     }
 
-    bool isReversible(size_t i) override;
     //! @}
 
     //! @name Reaction Mechanism Setup Routines
@@ -148,13 +146,6 @@ class BulkKinetics : public Kinetics
 
     //! @}
 
-    //! Vector of rate handlers
-    vector<unique_ptr<MultiRateBase>> m_bulk_rates;
-    map<string, size_t> m_bulk_types; //!< Mapping of rate handlers
-
-    vector<size_t> m_revindex; //!< Indices of reversible reactions
-    vector<size_t> m_irrev; //!< Indices of irreversible reactions
-
     //! Difference between the global reactants order and the global products
     //! order. Of type "double" to account for the fact that we can have real-
     //! valued stoichiometries.

include/cantera/kinetics/ElectronCollisionPlasmaRate.h

@@ -121,6 +121,9 @@ class ElectronCollisionPlasmaRate : public ReactionRate
      */
     double evalFromStruct(const ElectronCollisionPlasmaData& shared_data);
 
+    void modifyRateConstants(const ElectronCollisionPlasmaData& shared_data,
+                             double& kf, double& kr);
+
     //! Evaluate derivative of reaction rate with respect to temperature
     //! divided by reaction rate
     //! @param shared_data  data shared by all reactions of a given type
@@ -155,6 +158,9 @@ class ElectronCollisionPlasmaRate : public ReactionRate
 
     //! collision cross sections [m2] after interpolation
     vector<double> m_crossSectionsInterpolated;
+
+    //! collision cross section [m2] for super-elastic collision (interpolated)
+    Eigen::ArrayXd m_superElasticCrossSections;
 };
 
 }

include/cantera/kinetics/InterfaceKinetics.h

@@ -9,7 +9,6 @@
 #define CT_IFACEKINETICS_H
 
 #include "Kinetics.h"
-#include "MultiRate.h"
 
 namespace Cantera
 {
@@ -100,16 +99,6 @@ class InterfaceKinetics : public Kinetics
     //! @{
 
     void getActivityConcentrations(double* const conc) override;
-
-    bool isReversible(size_t i) override {
-        if (std::find(m_revindex.begin(), m_revindex.end(), i)
-                < m_revindex.end()) {
-            return true;
-        } else {
-            return false;
-        }
-    }
-
     void getFwdRateConstants(double* kfwd) override;
     void getRevRateConstants(double* krev, bool doIrreversible=false) override;
 
@@ -328,25 +317,8 @@ class InterfaceKinetics : public Kinetics
     //! Temporary work vector of length m_kk
     vector<double> m_grt;
 
-    //! List of reactions numbers which are reversible reactions
-    /*!
-     *  This is a vector of reaction numbers. Each reaction in the list is
-     *  reversible. Length = number of reversible reactions
-     */
-    vector<size_t> m_revindex;
-
     bool m_redo_rates = false;
 
-    //! Vector of rate handlers for interface reactions
-    vector<unique_ptr<MultiRateBase>> m_interfaceRates;
-    map<string, size_t> m_interfaceTypes; //!< Rate handler mapping
-
-    //! Vector of irreversible reaction numbers
-    /*!
-     * vector containing the reaction numbers of irreversible reactions.
-     */
-    vector<size_t> m_irrev;
-
     //! Array of concentrations for each species in the kinetics mechanism
     /*!
      * An array of generalized concentrations @f$ C_k @f$ that are defined

include/cantera/kinetics/Kinetics.h

@@ -13,6 +13,7 @@
 
 #include "StoichManager.h"
 #include "cantera/base/ValueCache.h"
+#include "MultiRate.h"
 
 namespace Cantera
 {
@@ -1203,7 +1204,7 @@ class Kinetics
      * @param i   reaction index
      */
     virtual bool isReversible(size_t i) {
-        throw NotImplementedError("Kinetics::isReversible");
+        return std::find(m_revindex.begin(), m_revindex.end(), i) < m_revindex.end();
     }
 
     /**
@@ -1446,6 +1447,10 @@ class Kinetics
      */
     double checkDuplicateStoich(map<int, double>& r1, map<int, double>& r2) const;
 
+    //! Vector of rate handlers
+    vector<unique_ptr<MultiRateBase>> m_rateHandlers;
+    map<string, size_t> m_rateTypes; //!< Mapping of rate handlers
+
     //! @name Stoichiometry management
     //!
     //! These objects and functions handle turning reaction extents into species
@@ -1527,6 +1532,9 @@ class Kinetics
     //! Net rate-of-progress for each reaction
     vector<double> m_ropnet;
 
+    vector<size_t> m_revindex; //!< Indices of reversible reactions
+    vector<size_t> m_irrev; //!< Indices of irreversible reactions
+
     //! The enthalpy change for each reaction to calculate Blowers-Masel rates
     vector<double> m_dH;
 

include/cantera/kinetics/MultiRate.h

@@ -21,6 +21,7 @@ template <class RateType, class DataType>
 class MultiRate final : public MultiRateBase
 {
     CT_DEFINE_HAS_MEMBER(has_update, updateFromStruct)
+    CT_DEFINE_HAS_MEMBER(has_modifyRateConstants, modifyRateConstants)
     CT_DEFINE_HAS_MEMBER(has_ddT, ddTScaledFromStruct)
     CT_DEFINE_HAS_MEMBER(has_ddP, perturbPressure)
     CT_DEFINE_HAS_MEMBER(has_ddM, perturbThirdBodies)
@@ -71,6 +72,14 @@ class MultiRate final : public MultiRateBase
         }
     }
 
+    void modifyRateConstants(double* kf, double* kr) override {
+        if constexpr (has_modifyRateConstants<RateType>::value) {
+            for (auto& [iRxn, rate] : m_rxn_rates) {
+                rate.modifyRateConstants(m_shared, kf[iRxn], kr[iRxn]);
+            }
+        }
+    }
+
     void processRateConstants_ddT(double* rop, const double* kf, double deltaT) override
     {
         if constexpr (has_ddT<RateType>::value) {

include/cantera/kinetics/MultiRateBase.h

@@ -52,6 +52,19 @@ class MultiRateBase
     //! @param kf  array of rate constants
     virtual void getRateConstants(double* kf) = 0;
 
+    //! For certain reaction types that do not follow mass action kinetics (for example,
+    //! Butler-Volmer), calculate modifications to the forward and reverse rate
+    //! constants.
+    //! @param[in, out] kf  On input, contains the rate constants as computed by
+    //!     getRateConstants(). The output value is updated by the reactant
+    //!     StoichManagerN to determine the forward rates of progress.
+    //! @param[in, out] kr  On input, contains the reverse rate constants computed from
+    //!     the forward rate constants and the equilibrium constants. The output value
+    //!     is updated by the product StoichManagerN to determine the reverse rates of
+    //!     progress.
+    //! @since New in %Cantera 3.1
+    virtual void modifyRateConstants(double* kf, double* kr) = 0;
+
     //! Evaluate all rate constant temperature derivatives handled by the evaluator;
     //! which are multiplied with the array of rate-of-progress variables.
     //! Depending on the implementation of a rate object, either an exact derivative or

include/cantera/kinetics/ReactionRate.h

@@ -26,9 +26,13 @@ class Reaction;
 //! In addition to the pure virtual methods declared in this class, complete derived
 //! classes must implement the method `evalFromStruct(const DataType& shared_data)`,
 //! where `DataType` is a container for parameters needed to evaluate reactions of that
-//! type. In addition, derived classes may also implement the method
-//! `updateFromStruct(const DataType& shared_data)` to update buffered data that
+//! type. In addition, derived classes may also implement the following methods:
+//! - `updateFromStruct(const DataType& shared_data)`, to update buffered data that
 //! is specific to a given reaction rate.
+//! - `modifyRateConstants(const DataType& shared_data, double& kf, double& kr)`, to
+//!   implement modifications to the forward and reverse rate constants for reactions
+//!   that do not follow mass action kinetics or use the equilibrium constant to
+//!   implement reversibility.
 //!
 //! The calculation of derivatives (or Jacobians) relies on the following methods:
 //!  -  Derived classes may implement the method
@@ -44,6 +48,8 @@ class Reaction;
 //!     which allow for the calculation of numerical derivatives.
 //!  -  For additional information, refer to the @ref kinDerivs "Kinetics Derivatives"
 //!     documentation.
+//!
+//! @since The `modifyRateConstants` method is new in Cantera 3.1.
 //! @ingroup reactionGroup
 class ReactionRate
 {

src/kinetics/BulkKinetics.cpp

@@ -12,10 +12,6 @@ BulkKinetics::BulkKinetics() {
     setDerivativeSettings(AnyMap()); // use default settings
 }
 
-bool BulkKinetics::isReversible(size_t i) {
-    return std::find(m_revindex.begin(), m_revindex.end(), i) < m_revindex.end();
-}
-
 bool BulkKinetics::addReaction(shared_ptr<Reaction> r, bool resize)
 {
     bool added = Kinetics::addReaction(r, resize);
@@ -33,32 +29,26 @@ bool BulkKinetics::addReaction(shared_ptr<Reaction> r, bool resize)
 
     m_dn.push_back(dn);
 
-    if (r->reversible) {
-        m_revindex.push_back(nReactions()-1);
-    } else {
-        m_irrev.push_back(nReactions()-1);
-    }
-
     shared_ptr<ReactionRate> rate = r->rate();
     string rtype = rate->subType();
     if (rtype == "") {
         rtype = rate->type();
     }
 
     // If necessary, add new MultiRate evaluator
-    if (m_bulk_types.find(rtype) == m_bulk_types.end()) {
-        m_bulk_types[rtype] = m_bulk_rates.size();
-        m_bulk_rates.push_back(rate->newMultiRate());
-        m_bulk_rates.back()->resize(m_kk, nReactions(), nPhases());
+    if (m_rateTypes.find(rtype) == m_rateTypes.end()) {
+        m_rateTypes[rtype] = m_rateHandlers.size();
+        m_rateHandlers.push_back(rate->newMultiRate());
+        m_rateHandlers.back()->resize(m_kk, nReactions(), nPhases());
     }
 
     // Set index of rate to number of reaction within kinetics
     rate->setRateIndex(nReactions() - 1);
     rate->setContext(*r, *this);
 
     // Add reaction rate to evaluator
-    size_t index = m_bulk_types[rtype];
-    m_bulk_rates[index]->add(nReactions() - 1, *rate);
+    size_t index = m_rateTypes[rtype];
+    m_rateHandlers[index]->add(nReactions() - 1, *rate);
 
     // Add reaction to third-body evaluator
     if (r->thirdBody() != nullptr) {
@@ -102,16 +92,16 @@ void BulkKinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
     }
 
     // Ensure that MultiRate evaluator is available
-    if (m_bulk_types.find(rtype) == m_bulk_types.end()) {
+    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_bulk_types[rtype];
+    size_t index = m_rateTypes[rtype];
     rate->setRateIndex(i);
     rate->setContext(*rNew, *this);
-    m_bulk_rates[index]->replace(i, *rate);
+    m_rateHandlers[index]->replace(i, *rate);
     invalidateCache();
 }
 
@@ -121,7 +111,7 @@ void BulkKinetics::resizeSpecies()
     m_act_conc.resize(m_kk);
     m_phys_conc.resize(m_kk);
     m_grt.resize(m_kk);
-    for (auto& rates : m_bulk_rates) {
+    for (auto& rates : m_rateHandlers) {
         rates->resize(m_kk, nReactions(), nPhases());
     }
 }
@@ -136,7 +126,7 @@ void BulkKinetics::resizeReactions()
     m_sbuf0.resize(nTotalSpecies());
     m_state.resize(thermo().stateSize());
     m_multi_concm.resizeCoeffs(nTotalSpecies(), nReactions());
-    for (auto& rates : m_bulk_rates) {
+    for (auto& rates : m_rateHandlers) {
         rates->resize(nTotalSpecies(), nReactions(), nPhases());
         // @todo ensure that ReactionData are updated; calling rates->update
         //      blocks correct behavior in update_rates_T
@@ -495,7 +485,7 @@ void BulkKinetics::updateROP()
     }
 
     // loop over MultiRate evaluators for each reaction type
-    for (auto& rates : m_bulk_rates) {
+    for (auto& rates : m_rateHandlers) {
         bool changed = rates->update(thermo(), *this);
         if (changed) {
             rates->getRateConstants(m_kf0.data());
@@ -518,12 +508,17 @@ void BulkKinetics::updateROP()
     processThirdBodies(m_ropf.data());
     copy(m_ropf.begin(), m_ropf.end(), m_ropr.begin());
 
-    // multiply ropf by concentration products
-    m_reactantStoich.multiply(m_act_conc.data(), m_ropf.data());
-
     // for reversible reactions, multiply ropr by concentration products
     applyEquilibriumConstants(m_ropr.data());
+
+    for (auto& rates : m_rateHandlers) {
+        rates->modifyRateConstants(m_ropf.data(), m_ropr.data());
+    }
+
+    // multiply ropf and ropr by concentration products
+    m_reactantStoich.multiply(m_act_conc.data(), m_ropf.data());
     m_revProductStoich.multiply(m_act_conc.data(), m_ropr.data());
+
     for (size_t j = 0; j != nReactions(); ++j) {
         m_ropnet[j] = m_ropf[j] - m_ropr[j];
     }
@@ -601,7 +596,7 @@ void BulkKinetics::process_ddT(const vector<double>& in, double* drop)
 {
     // apply temperature derivative
     copy(in.begin(), in.end(), drop);
-    for (auto& rates : m_bulk_rates) {
+    for (auto& rates : m_rateHandlers) {
         rates->processRateConstants_ddT(drop, m_rfn.data(), m_jac_rtol_delta);
     }
 }
@@ -610,7 +605,7 @@ void BulkKinetics::process_ddP(const vector<double>& in, double* drop)
 {
     // apply pressure derivative
     copy(in.begin(), in.end(), drop);
-    for (auto& rates : m_bulk_rates) {
+    for (auto& rates : m_rateHandlers) {
         rates->processRateConstants_ddP(drop, m_rfn.data(), m_jac_rtol_delta);
     }
 }
@@ -641,7 +636,7 @@ void BulkKinetics::process_ddC(StoichManagerN& stoich, const vector<double>& in,
     // derivatives due to reaction rates depending on third-body colliders
     if (!m_jac_skip_falloff) {
         m_multi_concm.scaleM(in.data(), outM.data(), m_concm.data(), ctot_inv);
-        for (auto& rates : m_bulk_rates) {
+        for (auto& rates : m_rateHandlers) {
             // processing step assigns zeros to entries not dependent on M
             rates->processRateConstants_ddM(
                 outM.data(), m_rfn.data(), m_jac_rtol_delta);
@@ -680,7 +675,7 @@ Eigen::SparseMatrix<double> BulkKinetics::calculateCompositionDerivatives(
 
     // derivatives due to reaction rates depending on third-body colliders
     if (!m_jac_skip_falloff) {
-        for (auto& rates : m_bulk_rates) {
+        for (auto& rates : m_rateHandlers) {
             // processing step does not modify entries not dependent on M
             rates->processRateConstants_ddM(
                 outV.data(), m_rfn.data(), m_jac_rtol_delta, false);