[Thermo] Add ThermoPhase.standardConcentrationUnits method

This method returns the units of the concentration-like terms appearing
in rate expressions, and are needed in order to convert rate constants
from user-specified input units to Cantera's MKS+kmol system.

include/cantera/base/Units.h

@@ -51,14 +51,14 @@ public:
     //! Provide a string representation of these Units
     std::string str() const;
 
-private:
-    //! Scale the unit by the factor `k`
-    void scale(double k) { m_factor *= k; }
-
     //! Raise these Units to a power, changing both the conversion factor and
     //! the dimensions of these Units.
     Units pow(double expoonent) const;
 
+private:
+    //! Scale the unit by the factor `k`
+    void scale(double k) { m_factor *= k; }
+
     double m_factor; //!< conversion factor to Cantera base units
     double m_mass_dim;
     double m_length_dim;

include/cantera/thermo/FixedChemPotSSTP.h

@@ -215,6 +215,8 @@ public:
      * @{
      */
 
+    virtual Units standardConcentrationUnits() const;
+
     //! @copydoc ThermoPhase::getActivityConcentrations
     /*! 
      * For a stoichiometric substance, there is only one species, and the

include/cantera/thermo/GibbsExcessVPSSTP.h

@@ -131,6 +131,7 @@ public:
      * @{
      */
 
+    virtual Units standardConcentrationUnits() const;
     virtual void getActivityConcentrations(doublereal* c) const;
 
     /**

include/cantera/thermo/IdealMolalSoln.h

@@ -270,6 +270,7 @@ public:
      * @{
      */
 
+    virtual Units standardConcentrationUnits() const;
     virtual void getActivityConcentrations(doublereal* c) const;
     virtual doublereal standardConcentration(size_t k=0) const;
 

include/cantera/thermo/IdealSolidSolnPhase.h

@@ -263,6 +263,8 @@ public:
      * @{
      */
 
+    virtual Units standardConcentrationUnits() const;
+
     /**
      * This method returns the array of generalized concentrations. The
      * generalized concentrations are used in the evaluation of the rates of

include/cantera/thermo/IdealSolnGasVPSS.h

@@ -96,6 +96,7 @@ protected:
     //! @}
 
 public:
+    virtual Units standardConcentrationUnits() const;
     virtual void getActivityConcentrations(doublereal* c) const;
 
     //! Returns the standard concentration \f$ C^0_k \f$, which is used to

include/cantera/thermo/LatticePhase.h

@@ -377,6 +377,7 @@ public:
      */
     //@{
 
+    virtual Units standardConcentrationUnits() const;
     virtual void getActivityConcentrations(doublereal* c) const;
 
     //! Return the standard concentration for the kth species

include/cantera/thermo/LatticeSolidPhase.h

@@ -303,6 +303,8 @@ public:
         throw NotImplementedError("LatticeSolidPhase::setConcentrations");
     }
 
+    virtual Units standardConcentrationUnits() const;
+
     virtual void getActivityConcentrations(doublereal* c) const;
 
     virtual void getActivityCoefficients(doublereal* ac) const;

include/cantera/thermo/MaskellSolidSolnPhase.h

@@ -33,6 +33,7 @@ public:
         return "MaskellSolidsoln";
     }
 
+    virtual Units standardConcentrationUnits() const { return Units(1.0); }
     virtual void getActivityConcentrations(doublereal* c) const;
     virtual doublereal standardConcentration(size_t k=0) const { return 1.0; }
     virtual doublereal logStandardConc(size_t k=0) const { return 0.0; }

include/cantera/thermo/MetalPhase.h

@@ -86,6 +86,10 @@ public:
         }
     }
 
+    virtual Units standardConcentrationUnits() const {
+        return Units(1.0);
+    }
+
     virtual doublereal standardConcentration(size_t k=0) const {
         return 1.0;
     }

include/cantera/thermo/PureFluidPhase.h

@@ -81,6 +81,7 @@ public:
     virtual void getPartialMolarCp(doublereal* cpbar) const;
     virtual void getPartialMolarVolumes(doublereal* vbar) const;
 
+    virtual Units standardConcentrationUnits() const;
     virtual void getActivityConcentrations(doublereal* c) const;
     virtual doublereal standardConcentration(size_t k=0) const;
 

include/cantera/thermo/StoichSubstance.h

@@ -205,6 +205,8 @@ public:
      * @{
      */
 
+    virtual Units standardConcentrationUnits() const;
+
     //! This method returns an array of generalized concentrations
     /*!
      * \f$ C^a_k\f$ are defined such that \f$ a_k = C^a_k / C^0_k, \f$ where

include/cantera/thermo/ThermoPhase.h

@@ -353,6 +353,16 @@ public:
      */
     virtual int standardStateConvention() const;
 
+    //! Returns the units of the "standard concentration" for this phase
+    /*!
+     * These are the units of the values returned by the functions
+     * getActivityConcentrations() and standardConcentration(), which can
+     * vary between different ThermoPhase-derived classes, or change within
+     * a single class depending on input options. See the documentation for
+     * standardConcentration() for the derived class for specific details.
+     */
+    virtual Units standardConcentrationUnits() const;
+
     //! This method returns an array of generalized concentrations
     /*!
      * \f$ C^a_k\f$ are defined such that \f$ a_k = C^a_k / C^0_k, \f$ where

src/kinetics/Reaction.cpp

@@ -298,31 +298,28 @@ Units rateCoeffUnits(const Reaction& R, const Kinetics& kin,
     }
 
     // Determine the units of the rate coefficient
-    double reaction_phase_ndim = static_cast<double>(
-        kin.thermo(kin.reactionPhaseIndex()).nDim());
-    double len_dim = - reaction_phase_ndim;
-    double quantity_dim = 1.0;
+    Units rxn_phase_units = kin.thermo(kin.reactionPhaseIndex()).standardConcentrationUnits();
+    Units rcUnits = rxn_phase_units;
+    rcUnits *= Units(1.0, 0, 0, -1);
     for (const auto& order : R.orders) {
-        len_dim += order.second * kin.speciesPhase(order.first).nDim();
-        quantity_dim -= order.second;
+        const auto& phase = kin.speciesPhase(order.first);
+        rcUnits *= phase.standardConcentrationUnits().pow(-order.second);
     }
     for (const auto& stoich : R.reactants) {
         // Order for each reactant is the reactant stoichiometric coefficient,
         // unless already overridden by user-specified orders
         if (stoich.first == "M") {
-            len_dim += reaction_phase_ndim;
-            quantity_dim -= 1.0;
+            rcUnits *= rxn_phase_units.pow(-1);
         } else if (R.orders.find(stoich.first) == R.orders.end()) {
-            len_dim += stoich.second * kin.speciesPhase(stoich.first).nDim();
-            quantity_dim -= stoich.second;
+            const auto& phase = kin.speciesPhase(stoich.first);
+            rcUnits *= phase.standardConcentrationUnits().pow(-stoich.second);
         }
     }
 
     // Incorporate pressure dependence for low-pressure falloff and high-
     // pressure chemically-activated reaction limits
-    len_dim += pressure_dependence * reaction_phase_ndim;
-    quantity_dim -= pressure_dependence;
-    return Units(1.0, 0, len_dim, -1, 0, 0, quantity_dim);
+    rcUnits *= rxn_phase_units.pow(-pressure_dependence);
+    return rcUnits;
 }
 
 Arrhenius readArrhenius(const Reaction& R, const AnyValue& rate,

src/thermo/FixedChemPotSSTP.cpp

@@ -96,6 +96,11 @@ doublereal FixedChemPotSSTP::thermalExpansionCoeff() const
 
 // ---- Chemical Potentials and Activities ----
 
+Units FixedChemPotSSTP::standardConcentrationUnits() const
+{
+    return Units(1.0); // dimensionless
+}
+
 void FixedChemPotSSTP::getActivityConcentrations(doublereal* c) const
 {
     c[0] = 1.0;

src/thermo/GibbsExcessVPSSTP.cpp

@@ -41,6 +41,12 @@ void GibbsExcessVPSSTP::calcDensity()
 }
 
 // - Activities, Standard States, Activity Concentrations -----------
+
+Units GibbsExcessVPSSTP::standardConcentrationUnits() const
+{
+    return Units(1.0); // dimensionless
+}
+
 void GibbsExcessVPSSTP::getActivityConcentrations(doublereal* c) const
 {
     getActivities(c);

src/thermo/IdealMolalSoln.cpp

@@ -157,6 +157,16 @@ void IdealMolalSoln::setMolarDensity(const doublereal conc)
 
 // ------- Activities and Activity Concentrations
 
+Units IdealMolalSoln::standardConcentrationUnits() const
+{
+    if (m_formGC == 0) {
+        return Units(1.0); // dimensionless
+    } else {
+        // kmol/m^3 for bulk phases
+        return Units(1.0, 0, -static_cast<double>(nDim()), 0, 0, 0, 1);
+    }
+}
+
 void IdealMolalSoln::getActivityConcentrations(doublereal* c) const
 {
     if (m_formGC != 1) {

src/thermo/IdealSolidSolnPhase.cpp

@@ -124,6 +124,16 @@ void IdealSolidSolnPhase::compositionChanged()
 
 // Chemical Potentials and Activities
 
+Units IdealSolidSolnPhase::standardConcentrationUnits() const
+{
+    if (m_formGC == 0) {
+        return Units(1.0); // dimensionless
+    } else {
+        // kmol/m^3 for bulk phases
+        return Units(1.0, 0, -static_cast<double>(nDim()), 0, 0, 0, 1);
+    }
+}
+
 void IdealSolidSolnPhase::getActivityConcentrations(doublereal* c) const
 {
     const doublereal* const dtmp = moleFractdivMMW();

src/thermo/IdealSolnGasVPSS.cpp

@@ -123,6 +123,15 @@ doublereal IdealSolnGasVPSS::isothermalCompressibility() const
     return 0.0;
 }
 
+Units IdealSolnGasVPSS::standardConcentrationUnits() const
+{
+    if (m_idealGas || m_formGC != 0) {
+        return Units(1.0, 0, -3, 0, 0, 0, 1);
+    } else {
+        return Units(1.0);
+    }
+}
+
 void IdealSolnGasVPSS::getActivityConcentrations(doublereal* c) const
 {
     if (m_idealGas) {

src/thermo/LatticePhase.cpp

@@ -75,6 +75,11 @@ void LatticePhase::compositionChanged()
     calcDensity();
 }
 
+Units LatticePhase::standardConcentrationUnits() const
+{
+    return Units(1.0);
+}
+
 void LatticePhase::getActivityConcentrations(doublereal* c) const
 {
     getMoleFractions(c);

src/thermo/LatticeSolidPhase.cpp

@@ -116,6 +116,11 @@ doublereal LatticeSolidPhase::cp_mole() const
     return sum;
 }
 
+Units LatticeSolidPhase::standardConcentrationUnits() const
+{
+    return Units(1.0);
+}
+
 void LatticeSolidPhase::getActivityConcentrations(doublereal* c) const
 {
     _updateThermo();

src/thermo/PureFluidPhase.cpp

@@ -185,6 +185,11 @@ void PureFluidPhase::getPartialMolarVolumes(doublereal* vbar) const
     vbar[0] = 1.0 / molarDensity();
 }
 
+Units PureFluidPhase::standardConcentrationUnits() const
+{
+    return Units(1.0);
+}
+
 void PureFluidPhase::getActivityConcentrations(doublereal* c) const
 {
     c[0] = 1.0;

src/thermo/StoichSubstance.cpp

@@ -51,6 +51,11 @@ doublereal StoichSubstance::thermalExpansionCoeff() const
 
 // ---- Chemical Potentials and Activities ----
 
+Units StoichSubstance::standardConcentrationUnits() const
+{
+    return Units(1.0);
+}
+
 void StoichSubstance::getActivityConcentrations(doublereal* c) const
 {
     c[0] = 1.0;

src/thermo/ThermoPhase.cpp

@@ -61,6 +61,12 @@ int ThermoPhase::standardStateConvention() const
     return m_ssConvention;
 }
 
+Units ThermoPhase::standardConcentrationUnits() const
+{
+    // kmol/m^3 for bulk phases, kmol/m^2 for surface phases, etc.
+    return Units(1.0, 0, -static_cast<double>(nDim()), 0, 0, 0, 1);
+}
+
 doublereal ThermoPhase::logStandardConc(size_t k) const
 {
     return log(standardConcentration(k));