Make reactions extendable

SConstruct

@@ -1746,6 +1746,7 @@ def postBuildMessage(target, source, env):
     print("*******************************************************")
     print("Compilation completed successfully.\n")
     print("- To run the test suite, type 'scons test'.")
+    print("- To list available tests, type 'scons test-help'.")
     if env['googletest'] == 'none':
         print("  WARNING: You set the 'googletest' to 'none' and all it's tests will be skipped.")
     if os.name == 'nt':

include/cantera/kinetics/BulkKinetics.h

@@ -11,6 +11,7 @@
 
 #include "Kinetics.h"
 #include "RateCoeffMgr.h"
+#include "cantera/kinetics/MultiRate.h"
 
 namespace Cantera
 {
@@ -37,6 +38,8 @@ class BulkKinetics : public Kinetics
                                      bool doIrreversible = false);
 
     virtual bool addReaction(shared_ptr<Reaction> r);
+    virtual void modifyReaction(size_t i, shared_ptr<Reaction> rNew);
+
     virtual void resizeSpecies();
 
     virtual void setMultiplier(size_t i, double f);
@@ -46,6 +49,10 @@ class BulkKinetics : public Kinetics
     virtual void addElementaryReaction(ElementaryReaction& r);
     virtual void modifyElementaryReaction(size_t i, ElementaryReaction& rNew);
 
+    //! Vector of rate handlers
+    std::vector<unique_ptr<MultiRateBase>> m_bulk_rates;
+    std::map<std::string, size_t> m_bulk_types; //!< Mapping of rate handlers
+
     Rate1<Arrhenius> m_rates;
     std::vector<size_t> m_revindex; //!< Indices of reversible reactions
     std::vector<size_t> m_irrev; //!< Indices of irreversible reactions

include/cantera/kinetics/FalloffMgr.h

@@ -9,6 +9,7 @@
 #define CT_FALLOFFMGR_H
 
 #include "Falloff.h"
+#include "cantera/base/global.h"
 
 namespace Cantera
 {
@@ -32,13 +33,35 @@ class FalloffMgr
      *     determine which array entry is modified in method pr_to_falloff.
      * @param reactionType Either `FALLOFF_RXN` or `CHEMACT_RXN`
      * @param f The falloff function.
+     *
+     * @deprecated To be removed after Cantera 2.6.
      */
     void install(size_t rxn, int reactionType, shared_ptr<Falloff> f) {
+        warn_deprecated("FalloffMgr::install()",
+            "To be removed after Cantera 2.6. Specify reaction type using "
+            "string instead.");
+
+        m_rxn.push_back(rxn);
+        m_offset.push_back(m_worksize);
+        m_worksize += f->workSize();
+        m_falloff.push_back(f);
+        m_isfalloff.push_back(reactionType == FALLOFF_RXN);
+        m_indices[rxn] = m_falloff.size()-1;
+    }
+
+    //! Install a new falloff function calculator.
+    /*
+     * @param rxn Index of the falloff reaction. This will be used to
+     *     determine which array entry is modified in method pr_to_falloff.
+     * @param type Reaction type identifier.
+     * @param f The falloff function.
+     */
+     void install(size_t rxn, std::string type, shared_ptr<Falloff> f) {
         m_rxn.push_back(rxn);
         m_offset.push_back(m_worksize);
         m_worksize += f->workSize();
         m_falloff.push_back(f);
-        m_reactionType.push_back(reactionType);
+        m_isfalloff.push_back(type == "falloff");
         m_indices[rxn] = m_falloff.size()-1;
     }
 
@@ -76,14 +99,14 @@ class FalloffMgr
     void pr_to_falloff(doublereal* values, const doublereal* work) {
         for (size_t i = 0; i < m_rxn.size(); i++) {
             double pr = values[m_rxn[i]];
-            if (m_reactionType[i] == FALLOFF_RXN) {
+            if (m_isfalloff[i]) {
                 // Pr / (1 + Pr) * F
                 values[m_rxn[i]] *=
-                    m_falloff[i]->F(pr, work + m_offset[i]) /(1.0 + pr);
+                    m_falloff[i]->F(pr, work + m_offset[i]) / (1.0 + pr);
             } else {
                 // 1 / (1 + Pr) * F
                 values[m_rxn[i]] =
-                    m_falloff[i]->F(pr, work + m_offset[i]) /(1.0 + pr);
+                    m_falloff[i]->F(pr, work + m_offset[i]) / (1.0 + pr);
             }
         }
     }
@@ -96,7 +119,7 @@ class FalloffMgr
     size_t m_worksize;
 
     //! Distinguish between falloff and chemically activated reactions
-    vector_int m_reactionType;
+    std::vector<bool> m_isfalloff;
 
     //! map of external reaction index to local index
     std::map<size_t, size_t> m_indices;

include/cantera/kinetics/Kinetics.h

@@ -13,6 +13,7 @@
 
 #include "StoichManager.h"
 #include "cantera/base/ValueCache.h"
+#include "cantera/kinetics/ReactionFactory.h"
 
 namespace Cantera
 {
@@ -601,9 +602,18 @@ class Kinetics
      * are specific to the particular kinetics manager.
      *
      * @param i   reaction index
+     *
+     * @deprecated To be changed after Cantera 2.6.
      */
     virtual int reactionType(size_t i) const;
 
+    /**
+     * String specifying the type of reaction.
+     *
+     * @param i   reaction index
+     */
+    virtual std::string reactionTypeStr(size_t i) const;
+
     /**
      * True if reaction i has been declared to be reversible. If isReversible(i)
      * is false, then the reverse rate of progress for reaction i is always

include/cantera/kinetics/MultiRate.h

@@ -0,0 +1,127 @@
+/**
+ * @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
+// at https://cantera.org/license.txt for license and copyright information.
+
+#ifndef CT_MULTIRATE_H
+#define CT_MULTIRATE_H
+
+#include "cantera/kinetics/ReactionRate.h"
+
+namespace Cantera
+{
+
+//! An abstract base class for evaluating all reactions of a particular type.
+/**
+ * Because this class has no template parameters, the `Kinetics` object
+ * can store all of these rate coefficient evaluators as a
+ * `vector<shared_ptr<MultiRateBase>>`.
+ *
+ * @todo At the moment, implemented methods are specific to `BulkKinetics`,
+ *     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
+{
+public:
+    virtual ~MultiRateBase() {}
+
+    //! Add reaction rate object to the evaluator
+    //! @param rxn_index  index of reaction
+    //! @param rate  reaction rate object
+    virtual void add(const size_t rxn_index,
+                     ReactionRateBase& rate) = 0;
+
+    //! Replace reaction rate object handled by the evaluator
+    //! @param rxn_index  index of reaction
+    //! @param rate  reaction rate object
+    virtual bool replace(const size_t rxn_index,
+                         ReactionRateBase& rate) = 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;
+
+    //! Update data common to reaction rates of a specific type
+    //! @param bulk  object representing bulk phase
+    virtual void update(const ThermoPhase& bulk) = 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);
+    }
+
+    virtual bool replace(const size_t rxn_index,
+                         ReactionRateBase& rate) override {
+        if (!m_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());
+        }
+        if (m_indices.find(rxn_index) != m_indices.end()) {
+            size_t j = m_indices[rxn_index];
+            m_rates[j] = dynamic_cast<RateType&>(rate);
+            return true;
+        }
+        return false;
+    }
+
+    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);
+        }
+    }
+
+    virtual void update(const ThermoPhase& bulk) override {
+        // update common data once for each reaction type
+        m_shared.update(bulk);
+        if (RateType::uses_update()) {
+            // 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);
+            }
+        }
+    }
+
+protected:
+    std::vector<RateType> m_rates; //! Reaction rate objects
+    std::vector<size_t> m_rxn; //! Index within overall rate vector
+    std::map<size_t, size_t> m_indices; //! Mapping of indices
+    DataType m_shared;
+};
+
+}
+
+#endif