Distinguish between rate units and units in rate parameterization

These differ specifically in the case of sticking reactions, where the
rate coefficient has units like m^3/kmol/s, depending on the reactant
orders while the sticking coefficient itself is dimensionless.

include/cantera/kinetics/Arrhenius.h

@@ -121,24 +121,18 @@ class ArrheniusBase : public ReactionRate
         return m_Ea_R * GasConstant;
     }
 
-    // Return units of the reaction rate expression
-    const Units& rateUnits() const {
-        return m_rate_units;
-    }
-
     //! Return reaction order associated with the reaction rate
     double order() const {
         return m_order;
     }
 
     //! Set units of the reaction rate expression
-    void setRateUnits(const UnitStack& rate_units) {
+    void setRateUnits(const UnitStack& rate_units) override {
+        ReactionRate::setRateUnits(rate_units);
         if (rate_units.size() > 1) {
-            m_rate_units = rate_units.product();
-            m_order = 1 - m_rate_units.dimension("quantity");
+            m_order = 1 - rate_units.product().dimension("quantity");
         } else {
             m_order = NAN;
-            m_rate_units = rate_units.standardUnits();
         }
     }
 
@@ -164,7 +158,6 @@ class ArrheniusBase : public ReactionRate
     std::string m_b_str = "b"; //!< The string for temperature exponent
     std::string m_Ea_str = "Ea"; //!< The string for activation energy
     std::string m_E4_str = ""; //!< The string for an optional 4th parameter
-    Units m_rate_units{0.}; //!< Reaction rate units
 };
 
 //! Arrhenius reaction rate type depends only on temperature

include/cantera/kinetics/ChebyshevRate.h

@@ -232,8 +232,6 @@ class ChebyshevRate final : public ReactionRate
 
     Array2D m_coeffs; //!<< coefficient array
     vector_fp dotProd_; //!< dot product of coeffs with the reduced pressure polynomial
-
-    Units m_rate_units = Units(0.); //!< Reaction rate units
 };
 
 }

include/cantera/kinetics/InterfaceRate.h

@@ -456,12 +456,10 @@ class StickingRate : public RateType, public StickingCoverage
 
     //! Constructor based on AnyMap content
     StickingRate(const AnyMap& node, const UnitStack& rate_units) {
-        // sticking coefficients are dimensionless
-        setParameters(node, Units(1.0));
+        setParameters(node, rate_units);
     }
     explicit StickingRate(const AnyMap& node) {
-        // sticking coefficients are dimensionless
-        setParameters(node, Units(1.0));
+        setParameters(node, {});
     }
 
     unique_ptr<MultiRateBase> newMultiRate() const override {
@@ -473,10 +471,16 @@ class StickingRate : public RateType, public StickingCoverage
         return "sticking-" + RateType::type();
     }
 
+    void setRateUnits(const UnitStack& rate_units) override {
+        // Sticking coefficients are dimensionless
+        RateType::m_conversion_units = Units(1.0);
+    }
+
     virtual void setParameters(
         const AnyMap& node, const UnitStack& rate_units) override
     {
         InterfaceRateBase::setParameters(node);
+        setRateUnits(rate_units);
         RateType::m_negativeA_ok = node.getBool("negative-A", false);
         setStickingParameters(node);
         if (!node.hasKey("sticking-coefficient")) {

include/cantera/kinetics/ReactionRate.h

@@ -53,6 +53,7 @@ class ReactionRate
         : m_input(other.m_input)
         , m_rate_index(other.m_rate_index)
         , m_valid(other.m_valid)
+        , m_conversion_units(other.m_conversion_units)
     {}
 
     ReactionRate& operator=(const ReactionRate& other) {
@@ -62,6 +63,7 @@ class ReactionRate
         m_input = other.m_input;
         m_rate_index = other.m_rate_index;
         m_valid = other.m_valid;
+        m_conversion_units = other.m_conversion_units;
         return *this;
     }
 
@@ -94,6 +96,7 @@ class ReactionRate
     //! @param node  AnyMap object containing reaction rate specification
     //! @param units  unit definitions specific to rate information
     virtual void setParameters(const AnyMap& node, const UnitStack& units) {
+        setRateUnits(units);
         m_input = node;
     }
 
@@ -107,6 +110,30 @@ class ReactionRate
         return out;
     }
 
+    //! Get the units for converting the leading term in the reaction rate expression.
+    //!
+    //! These units are often the same as the units of the rate expression itself, but
+    //! not always; sticking coefficients are a notable exception.
+    //! @since New in Cantera 3.0
+    const Units& conversionUnits() const {
+        return m_conversion_units;
+    }
+
+    //! Set the units of the reaction rate expression
+    //!
+    //! Used to determine the units that should be used for converting terms in the
+    //! reaction rate expression, which often have the same units (for example, the
+    //! Arrhenius pre-exponential) but may also be different (for example, sticking
+    //! coefficients).
+    //! @since New in Cantera 3.0
+    virtual void setRateUnits(const UnitStack& rate_units) {
+        if (rate_units.size() > 1) {
+            m_conversion_units = rate_units.product();
+        } else {
+            m_conversion_units = rate_units.standardUnits();
+        }
+    }
+
     //! Check basic syntax and settings of reaction rate expression
     virtual void check(const std::string& equation) {}
 
@@ -225,6 +252,9 @@ class ReactionRate
     //! Flag indicating composition dependent rate
     bool m_composition_dependent_rate = false;
 
+    //! Units of the leading term in the reaction rate expression
+    Units m_conversion_units{0.};
+
 private:
     //! Return an object that be used to evaluate the rate by converting general input
     //! such as temperature and pressure into the `DataType` struct that is particular

interfaces/cython/cantera/reaction.pxd

@@ -38,6 +38,7 @@ cdef extern from "cantera/kinetics/ReactionRate.h" namespace "Cantera":
         double eval(double, double) except +translate_exception
         double eval(double, vector[double]&) except +translate_exception
         CxxAnyMap parameters() except +translate_exception
+        CxxUnits conversionUnits()
 
 
 cdef extern from "cantera/kinetics/Arrhenius.h" namespace "Cantera":

interfaces/cython/cantera/reaction.pyx

@@ -147,6 +147,14 @@ cdef class ReactionRate:
         def __get__(self):
             return anymap_to_py(self.rate.parameters())
 
+    @property
+    def conversion_units(self) -> Units:
+        """
+        Get the units for converting the leading term in the reaction rate expression
+        to different unit systems.
+        """
+        return Units.copy(self.rate.conversionUnits())
+
 
 cdef class ArrheniusRateBase(ReactionRate):
     """

src/kinetics/Arrhenius.cpp

@@ -23,6 +23,7 @@ ArrheniusBase::ArrheniusBase(double A, double b, double Ea)
 ArrheniusBase::ArrheniusBase(const AnyValue& rate, const UnitSystem& units,
                              const UnitStack& rate_units)
 {
+    setRateUnits(rate_units);
     setRateParameters(rate, units, rate_units);
 }
 
@@ -40,16 +41,14 @@ void ArrheniusBase::setRateParameters(
         m_A = NAN;
         m_b = NAN;
         m_logA = NAN;
-        m_order = NAN;
-        m_rate_units = Units(0.);
+        setRateUnits(Units(0.));
         return;
     }
 
-    setRateUnits(rate_units);
     if (rate.is<AnyMap>()) {
 
         auto& rate_map = rate.as<AnyMap>();
-        m_A = units.convertRateCoeff(rate_map[m_A_str], m_rate_units);
+        m_A = units.convertRateCoeff(rate_map[m_A_str], conversionUnits());
         m_b = rate_map[m_b_str].asDouble();
         if (rate_map.hasKey(m_Ea_str)) {
             m_Ea_R = units.convertActivationEnergy(rate_map[m_Ea_str], "K");
@@ -59,7 +58,7 @@ void ArrheniusBase::setRateParameters(
         }
     } else {
         auto& rate_vec = rate.asVector<AnyValue>(2, 4);
-        m_A = units.convertRateCoeff(rate_vec[0], m_rate_units);
+        m_A = units.convertRateCoeff(rate_vec[0], conversionUnits());
         m_b = rate_vec[1].asDouble();
         if (rate_vec.size() > 2) {
             m_Ea_R = units.convertActivationEnergy(rate_vec[2], "K");
@@ -81,8 +80,8 @@ void ArrheniusBase::getRateParameters(AnyMap& node) const
         return;
     }
 
-    if (m_rate_units.factor() != 0.0) {
-        node[m_A_str].setQuantity(m_A, m_rate_units);
+    if (conversionUnits().factor() != 0.0) {
+        node[m_A_str].setQuantity(m_A, conversionUnits());
     } else {
         node[m_A_str] = m_A;
         // This can't be converted to a different unit system because the dimensions of

src/kinetics/ChebyshevRate.cpp

@@ -63,7 +63,6 @@ ChebyshevRate::ChebyshevRate(const AnyMap& node, const UnitStack& rate_units)
 void ChebyshevRate::setParameters(const AnyMap& node, const UnitStack& rate_units)
 {
     ReactionRate::setParameters(node, rate_units);
-    m_rate_units = rate_units.product();
     const UnitSystem& unit_system = node.units();
     Array2D coeffs(0, 0);
     if (node.hasKey("data")) {
@@ -80,7 +79,7 @@ void ChebyshevRate::setParameters(const AnyMap& node, const UnitStack& rate_unit
                 coeffs(i, j) = vcoeffs[i][j];
             }
         }
-        double offset = unit_system.convertRateCoeff(AnyValue(1.0), m_rate_units);
+        double offset = unit_system.convertRateCoeff(AnyValue(1.0), conversionUnits());
         coeffs(0, 0) += std::log10(offset);
         setLimits(
             unit_system.convert(T_range[0], "K"),
@@ -142,7 +141,7 @@ void ChebyshevRate::getParameters(AnyMap& rateNode) const
     }
     // Unit conversions must take place later, after the destination unit system
     // is known. A lambda function is used here to override the default behavior
-    Units rate_units2 = m_rate_units;
+    Units rate_units2 = conversionUnits();
     auto converter = [rate_units2](AnyValue& coeffs, const UnitSystem& units) {
         if (rate_units2.factor() != 0.0) {
             coeffs.asVector<vector_fp>()[0][0] += \

src/kinetics/Falloff.cpp

@@ -344,7 +344,7 @@ void TroeRate::getParameters(AnyMap& node) const
     AnyMap params;
     if (!valid()) {
         // pass
-    } else if (m_lowRate.rateUnits().factor() != 0.0) {
+    } else if (m_lowRate.conversionUnits().factor() != 0.0) {
         params["A"] = m_a;
         params["T3"].setQuantity(1.0 / m_rt3, "K");
         params["T1"].setQuantity(1.0 / m_rt1, "K");
@@ -472,7 +472,7 @@ void SriRate::getParameters(AnyMap& node) const
     AnyMap params;
     if (!valid()) {
         // pass
-    } else if (m_lowRate.rateUnits().factor() != 0.0) {
+    } else if (m_lowRate.conversionUnits().factor() != 0.0) {
         params["A"] = m_a;
         params["B"].setQuantity(m_b, "K");
         params["C"].setQuantity(m_c, "K");

src/kinetics/Kinetics.cpp

@@ -675,7 +675,7 @@ bool Kinetics::addReaction(shared_ptr<Reaction> r, bool resize)
     // For reactions created outside the context of a Kinetics object, the units
     // of the rate coefficient can't be determined in advance. Do that here.
     if (r->rate_units.factor() == 0) {
-        r->calculateRateCoeffUnits(*this);
+        r->rate()->setRateUnits(r->calculateRateCoeffUnits(*this));
     }
 
     size_t irxn = nReactions(); // index of the new reaction

test/python/test_reaction.py

@@ -1111,6 +1111,10 @@ def test_rate_coeff_units(self):
         rxn = self.from_yaml()
         assert str(rxn.rate_coeff_units) == str(self._rc_units)
 
+    def test_rate_conversion_units(self):
+        rxn = self.from_yaml()
+        assert str(rxn.rate.conversion_units) == str(self._rc_units)
+
     def check_equal(self, one, two):
         # helper function for deprecation tests
         self.assertEqual(type(one), type(two))
@@ -1883,6 +1887,10 @@ def test_site_density(self):
         self.assertEqual(self.soln.site_density,
             self.soln.reaction(self._index).rate.site_density)
 
+    def test_rate_conversion_units(self):
+        rxn = self.from_yaml()
+        assert str(rxn.rate.conversion_units) == str(ct.Units('1'))
+
 
 class TestArrheniusStickReaction(StickReactionTests, utilities.CanteraTest):
     # test interface reaction without coverages