[Kinetics] make RxnRate thread-safe

include/cantera/kinetics/RxnRates.h

@@ -445,59 +445,8 @@ public:
     //! Identifier of reaction type
     virtual std::string type() const = 0;
 
-    //! Get temperature
-    static double temperature() {
-        return m_temperature;
-    }
-
-    //! Update temperature
-    static void updateTemperature(double T) {
-        if (T <= 0.) {
-            throw CanteraError("RxnRate::updateTemperature",
-                               "Temperature has to be positive.");
-        }
-        m_temperature = T;
-        m_logT = std::log(T);
-        m_recipT = 1./T;
-    }
-
-    //! Get pressure
-    static double pressure() {
-        return m_pressure;
-    }
-
-    //! Update pressure
-    static void updatePressure(double P) {
-        if (P <= 0.) {
-            throw CanteraError("RxnRate::updatePressure",
-                               "Pressure has to be positive.");
-        }
-        m_pressure = P;
-        m_logP = log(P);
-        m_log10P = log10(P);
-    }
-
     //! Evaluate reaction rate
-    virtual double eval() const = 0;
-
-protected:
-    //! Temperature used for reaction rate evaluation
-    static double m_temperature;
-
-    //! Logarithm of temperature used for reaction rate evaluation
-    static double m_logT;
-
-    //! Inverse temperature used for reaction rate evaluation
-    static double m_recipT;
-
-    //! Pressure used for reaction rate evaluation
-    static double m_pressure;
-
-    //! Logarithm of pressure used for reaction rate evaluation
-    static double m_logP;
-
-    //! 10-base logarithm of pressure used for reaction rate evaluation
-    static double m_log10P;
+    virtual double eval_T(double T, double logT, double recipT) const = 0;
 };
 
 
@@ -526,7 +475,7 @@ public:
      */
     void setRateFunction(shared_ptr<Func1> f);
 
-    virtual double eval() const;
+    virtual double eval_T(double T, double logT, double recipT) const;
 
 protected:
     shared_ptr<Func1> m_ratefunc;
@@ -552,8 +501,8 @@ public:
         return "ArrheniusRate";
     }
 
-    virtual double eval() const {
-        return updateRC(m_logT, m_recipT);
+    virtual double eval_T(double T, double logT, double recipT) const {
+        return updateRC(logT, recipT);
     }
 };
 

interfaces/cython/cantera/_cantera.pxd

@@ -315,11 +315,7 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera":
     cdef cppclass CxxRxnRate "Cantera::RxnRate":
         CxxRxnRate()
         string type()
-        double temperature()
-        void updateTemperature(double) except +translate_exception
-        double pressure()
-        void updatePressure(double) except +translate_exception
-        double eval() except +translate_exception
+        double eval_T(double, double, double) except +translate_exception
 
     cdef cppclass CxxCustomPyRate "Cantera::CustomPyRate" (CxxRxnRate):
         CxxCustomPyRate()

interfaces/cython/cantera/reaction.pyx

@@ -408,25 +408,10 @@ cdef class _RxnRate:
     def __repr__(self):
         return "<{}>".format(pystr(self.base.type()))
 
-    def __call__(self, temperature=None):
-        if temperature is not None:
-            self.temperature = temperature
-
-        return self.base.eval()
-
-    property temperature:
-        """ Get/Set temperature used for all reaction rate evaluations"""
-        def __get__(self):
-            return self.base.temperature()
-        def __set__(self, float temperature):
-            self.base.updateTemperature(temperature)
-
-    property pressure:
-        """ Get/Set pressure used for all reaction rate evaluations"""
-        def __get__(self):
-            return self.base.pressure()
-        def __set__(self, float pressure):
-            self.base.updatePressure(pressure)
+    def __call__(self, double temperature):
+        cdef double logT = np.log(temperature)
+        cdef double recipT = 1/temperature
+        return self.base.eval_T(temperature, logT, recipT)
 
 
 cdef class CustomRate(_RxnRate):

src/kinetics/GasKinetics.cpp

@@ -27,6 +27,7 @@ void GasKinetics::update_rates_T()
     double P = thermo().pressure();
     m_logStandConc = log(thermo().standardConcentration());
     double logT = log(T);
+    double recipT = 1./T;
 
     if (T != m_temp) {
         if (!m_rfn.empty()) {
@@ -45,14 +46,9 @@ void GasKinetics::update_rates_T()
     }
 
     if (T != m_temp || P != m_pres) {
-        if (!m_rxn_rates.empty()) {
-            RxnRate::updateTemperature(T);
-            RxnRate::updatePressure(P);
-        }
-
         for (auto& rate : m_rxn_rates) {
             // generic reaction rates
-            m_rfn.data()[rate.first] = rate.second->eval();
+            m_rfn.data()[rate.first] = rate.second->eval_T(T, logT, recipT);
         }
 
         if (m_plog_rates.nReactions()) {
@@ -244,8 +240,6 @@ bool GasKinetics::addReaction(shared_ptr<Reaction> r)
     if (r->rxnRate()) {
         // new generic reaction type handler
         m_rxn_rates[nReactions()-1] = r->rxnRate();
-        RxnRate::updateTemperature(thermo().temperature());
-        RxnRate::updatePressure(thermo().pressure());
     } else if (r->type() == "elementary") {
         addElementaryReaction(dynamic_cast<ElementaryReaction&>(*r));
     } else if (r->type() == "three-body") {
@@ -337,8 +331,6 @@ void GasKinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
     if (rNew->rxnRate()) {
         // new generic reaction type handler
         modifyRxnRate(i, rNew->rxnRate());
-        RxnRate::updateTemperature(thermo().temperature());
-        RxnRate::updatePressure(thermo().pressure());
     } else if (rNew->type() == "elementary") {
         modifyElementaryReaction(i, dynamic_cast<ElementaryReaction&>(*rNew));
     } else if (rNew->type() == "three-body") {

src/kinetics/RxnRates.cpp

@@ -159,22 +159,15 @@ ChebyshevRate::ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax,
     }
 }
 
-double RxnRate::m_temperature = 1.;
-double RxnRate::m_logT = 0.;
-double RxnRate::m_recipT = 1.;
-double RxnRate::m_pressure = 1.;
-double RxnRate::m_logP = 0.;
-double RxnRate::m_log10P = 0.;
-
 CustomPyRate::CustomPyRate() : m_ratefunc(0) {}
 
 void CustomPyRate::setRateFunction(shared_ptr<Func1> f) {
     m_ratefunc = f;
 }
 
-double CustomPyRate::eval() const {
+double CustomPyRate::eval_T(double T, double logT, double recipT) const {
     if (m_ratefunc) {
-        return m_ratefunc->eval(m_temperature);
+        return m_ratefunc->eval(T);
     }
     return 0.;
 }