Refactor reaction rate evaluators

include/cantera/kinetics/BulkKinetics.h

@@ -11,12 +11,13 @@
 
 #include "Kinetics.h"
 #include "RateCoeffMgr.h"
+#include "ThirdBodyCalc.h"
 #include "cantera/kinetics/MultiRate.h"
 
 namespace Cantera
 {
 
-class ElementaryReaction;
+class ElementaryReaction2;
 
 //! Partial specialization of Kinetics for chemistry in a single bulk phase
 class BulkKinetics : public Kinetics
@@ -45,9 +46,11 @@ class BulkKinetics : public Kinetics
     virtual void setMultiplier(size_t i, double f);
     virtual void invalidateCache();
 
+    void addThirdBody(shared_ptr<Reaction> r);
+
 protected:
-    virtual void addElementaryReaction(ElementaryReaction& r);
-    virtual void modifyElementaryReaction(size_t i, ElementaryReaction& rNew);
+    virtual void addElementaryReaction(ElementaryReaction2& r);
+    virtual void modifyElementaryReaction(size_t i, ElementaryReaction2& rNew);
 
     //! Vector of rate handlers
     std::vector<unique_ptr<MultiRateBase>> m_bulk_rates;
@@ -62,12 +65,19 @@ class BulkKinetics : public Kinetics
     //! valued stoichiometries.
     vector_fp m_dn;
 
-    //! Activity concentrations, as calculated by
-    //! ThermoPhase::getActivityConcentrations
+    ThirdBodyCalc m_multi_concm; //!< used with MultiRate evaluator
+    vector_fp concm_multi_values; //!< concentrations of third-body collision partners
+    std::vector<size_t> m_multi_indices; //!< reaction indices
+
+    //! Third body concentrations
+    vector_fp m_concm;
+
+    //! Activity concentrations, as calculated by ThermoPhase::getActivityConcentrations
     vector_fp m_act_conc;
 
     //! Physical concentrations, as calculated by ThermoPhase::getConcentrations
     vector_fp m_phys_conc;
+
     vector_fp m_grt;
 
     bool m_ROP_ok;

include/cantera/kinetics/GasKinetics.h

@@ -10,7 +10,6 @@
 #define CT_GASKINETICS_H
 
 #include "BulkKinetics.h"
-#include "ThirdBodyCalc.h"
 #include "FalloffMgr.h"
 #include "Reaction.h"
 
@@ -87,7 +86,7 @@ class GasKinetics : public BulkKinetics
     ThirdBodyCalc m_falloff_concm;
 
     Rate1<Plog> m_plog_rates;
-    Rate1<ChebyshevRate> m_cheb_rates;
+    Rate1<Chebyshev> m_cheb_rates;
     Rate1<BlowersMasel> m_blowersmasel_rates;
 
     //! @name Reaction rate data
@@ -107,16 +106,25 @@ class GasKinetics : public BulkKinetics
 
     void processFalloffReactions();
 
-    void addThreeBodyReaction(ThreeBodyReaction& r);
+    // functions marked as deprecated below are only used for XML import and
+    // transitional reaction types are marked as '-legacy'
+
+    //! @deprecated To be removed after Cantera 2.6 (replaced by MultiRate approach)
+    void addThreeBodyReaction(ThreeBodyReaction2& r);
     void addFalloffReaction(FalloffReaction& r);
-    void addPlogReaction(PlogReaction& r);
-    void addChebyshevReaction(ChebyshevReaction& r);
+    //! @deprecated To be removed after Cantera 2.6 (replaced by MultiRate approach)
+    void addPlogReaction(PlogReaction2& r);
+    //! @deprecated To be removed after Cantera 2.6 (replaced by MultiRate approach)
+    void addChebyshevReaction(ChebyshevReaction2& r);
     void addBlowersMaselReaction(BlowersMaselReaction& r);
 
-    void modifyThreeBodyReaction(size_t i, ThreeBodyReaction& r);
+    //! @deprecated To be removed after Cantera 2.6 (replaced by MultiRate approach)
+    void modifyThreeBodyReaction(size_t i, ThreeBodyReaction2& r);
     void modifyFalloffReaction(size_t i, FalloffReaction& r);
-    void modifyPlogReaction(size_t i, PlogReaction& r);
-    void modifyChebyshevReaction(size_t i, ChebyshevReaction& r);
+    //! @deprecated To be removed after Cantera 2.6 (replaced by MultiRate approach)
+    void modifyPlogReaction(size_t i, PlogReaction2& r);
+    //! @deprecated To be removed after Cantera 2.6 (replaced by MultiRate approach)
+    void modifyChebyshevReaction(size_t i, ChebyshevReaction2& r);
     void modifyBlowersMaselReaction(size_t i, BlowersMaselReaction& r);
 
     //! Update the equilibrium constants in molar units.

include/cantera/kinetics/MultiRate.h

@@ -1,8 +1,5 @@
 /**
  * @file MultiRate.h
- *
- * @warning This file is an experimental part of the %Cantera API and
- *     may be changed or removed without notice.
  */
 
 // This file is part of Cantera. See License.txt in the top-level directory or
@@ -26,9 +23,6 @@ namespace Cantera
  *     which can be updated using information of a single `ThermoPhase`.
  *     `InterfaceKinetics` will require access to an entire `Kinetics` object
  *     or the underlying `vector<ThermoPhase*>` vector (e.g. `m_thermo`).
- *
- * @warning This class is an experimental part of the %Cantera API and
- *     may be changed or removed without notice.
  */
 class MultiRateBase
 {
@@ -47,77 +41,87 @@ class MultiRateBase
     virtual bool replace(const size_t rxn_index,
                          ReactionRateBase& rate) = 0;
 
+    //! Update number of species
+    //! @param n_species  number of species
+    virtual void resizeSpecies(size_t n_species) = 0;
+
     //! Evaluate all rate constants handled by the evaluator
     //! @param bulk  object representing bulk phase
     //! @param kf  array of rate constants
     virtual void getRateConstants(const ThermoPhase& bulk,
-                                  double* kf) const = 0;
+                                  double* kf, double* concm) const = 0;
 
     //! Update data common to reaction rates of a specific type
     //! @param bulk  object representing bulk phase
-    virtual void update(const ThermoPhase& bulk) = 0;
+    //! @param concm  effective third-body concentrations
+    //! @TODO  enable more generic handling of non-trivial concentration dependencies
+    virtual void update(const ThermoPhase& bulk, double* concm) = 0;
 };
 
 
 //! A class template handling all reaction rates specific to `BulkKinetics`.
-/**
- * @warning This class is an experimental part of the %Cantera API and
- *     may be changed or removed without notice.
- */
 template <class RateType, class DataType>
 class MultiBulkRates final : public MultiRateBase
 {
 public:
     virtual void add(const size_t rxn_index,
-                     ReactionRateBase& rate) override {
-        m_indices[rxn_index] = m_rates.size();
-        m_rates.push_back(dynamic_cast<RateType&>(rate));
-        m_rxn.push_back(rxn_index);
+                     ReactionRateBase& rate) override
+    {
+        m_indices[rxn_index] = m_rxn_rates.size();
+        m_rxn_rates.emplace_back(rxn_index, dynamic_cast<RateType&>(rate));
     }
 
     virtual bool replace(const size_t rxn_index,
-                         ReactionRateBase& rate) override {
-        if (!m_rates.size()) {
+                         ReactionRateBase& rate) override
+    {
+        if (!m_rxn_rates.size()) {
             throw CanteraError("MultiBulkRate::replace",
                  "Invalid operation: cannot replace rate object "
                  "in empty rate handler.");
         } else if (typeid(rate) != typeid(RateType)) {
             throw CanteraError("MultiBulkRate::replace",
                  "Invalid operation: cannot replace rate object of type '{}' "
                  "with a new rate of type '{}'.",
-                 m_rates[0].type(), rate.type());
+                 m_rxn_rates.at(0).second.type(), rate.type());
         }
         if (m_indices.find(rxn_index) != m_indices.end()) {
             size_t j = m_indices[rxn_index];
-            m_rates[j] = dynamic_cast<RateType&>(rate);
+            m_rxn_rates.at(j).second = dynamic_cast<RateType&>(rate);
             return true;
         }
         return false;
     }
 
+    virtual void resizeSpecies(size_t n_species)
+    {
+        m_shared.resizeSpecies(n_species);
+    }
+
     virtual void getRateConstants(const ThermoPhase& bulk,
-                                  double* kf) const override {
-        for (size_t i = 0; i < m_rates.size(); i++) {
-            kf[m_rxn[i]] = m_rates[i].eval(m_shared);
+                                  double* kf, double* concm) const override
+    {
+        for (const auto& rxn : m_rxn_rates) {
+            kf[rxn.first] = rxn.second.eval(m_shared, concm[rxn.first]);
         }
     }
 
-    virtual void update(const ThermoPhase& bulk) override {
+    virtual void update(const ThermoPhase& bulk, double* concm) override
+    {
         // update common data once for each reaction type
         m_shared.update(bulk);
-        if (RateType::uses_update()) {
+        if (RateType::usesUpdate()) {
             // update reaction-specific data for each reaction. This loop
             // is efficient as all function calls are de-virtualized, and
             // all of the rate objects are contiguous in memory
-            for (auto& rate : m_rates) {
-                rate.update(m_shared, bulk);
+            for (auto& rxn : m_rxn_rates) {
+                rxn.second.update(m_shared, concm[rxn.first]);
             }
         }
     }
 
 protected:
-    std::vector<RateType> m_rates; //! Reaction rate objects
-    std::vector<size_t> m_rxn; //! Index within overall rate vector
+    //! Vector of pairs of reaction rates indices and reaction rates
+    std::vector<std::pair<size_t, RateType>> m_rxn_rates;
     std::map<size_t, size_t> m_indices; //! Mapping of indices
     DataType m_shared;
 };