From 53faf54ceb1d837029b08b10872f336d8128cd53 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Sat, 2 Mar 2019 14:50:54 -0500
Subject: [PATCH] [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                  |  8 +++----
 include/cantera/thermo/FixedChemPotSSTP.h     |  2 ++
 include/cantera/thermo/GibbsExcessVPSSTP.h    |  1 +
 include/cantera/thermo/IdealMolalSoln.h       |  1 +
 include/cantera/thermo/IdealSolidSolnPhase.h  |  2 ++
 include/cantera/thermo/IdealSolnGasVPSS.h     |  1 +
 include/cantera/thermo/LatticePhase.h         |  1 +
 include/cantera/thermo/LatticeSolidPhase.h    |  2 ++
 .../cantera/thermo/MaskellSolidSolnPhase.h    |  1 +
 include/cantera/thermo/MetalPhase.h           |  4 ++++
 include/cantera/thermo/PureFluidPhase.h       |  1 +
 include/cantera/thermo/StoichSubstance.h      |  2 ++
 include/cantera/thermo/ThermoPhase.h          | 10 ++++++++
 src/kinetics/Reaction.cpp                     | 23 ++++++++-----------
 src/thermo/FixedChemPotSSTP.cpp               |  5 ++++
 src/thermo/GibbsExcessVPSSTP.cpp              |  6 +++++
 src/thermo/IdealMolalSoln.cpp                 | 10 ++++++++
 src/thermo/IdealSolidSolnPhase.cpp            | 10 ++++++++
 src/thermo/IdealSolnGasVPSS.cpp               |  9 ++++++++
 src/thermo/LatticePhase.cpp                   |  5 ++++
 src/thermo/LatticeSolidPhase.cpp              |  5 ++++
 src/thermo/PureFluidPhase.cpp                 |  5 ++++
 src/thermo/StoichSubstance.cpp                |  5 ++++
 src/thermo/ThermoPhase.cpp                    |  6 +++++
 24 files changed, 108 insertions(+), 17 deletions(-)

diff --git a/include/cantera/base/Units.h b/include/cantera/base/Units.h
index 499ebf8d0..71ff11706 100644
--- a/include/cantera/base/Units.h
+++ b/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;
diff --git a/include/cantera/thermo/FixedChemPotSSTP.h b/include/cantera/thermo/FixedChemPotSSTP.h
index bdf8f9e97..08323f270 100644
--- a/include/cantera/thermo/FixedChemPotSSTP.h
+++ b/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
diff --git a/include/cantera/thermo/GibbsExcessVPSSTP.h b/include/cantera/thermo/GibbsExcessVPSSTP.h
index 41f160d8c..272ef1786 100644
--- a/include/cantera/thermo/GibbsExcessVPSSTP.h
+++ b/include/cantera/thermo/GibbsExcessVPSSTP.h
@@ -131,6 +131,7 @@ public:
      * @{
      */
 
+    virtual Units standardConcentrationUnits() const;
     virtual void getActivityConcentrations(doublereal* c) const;
 
     /**
diff --git a/include/cantera/thermo/IdealMolalSoln.h b/include/cantera/thermo/IdealMolalSoln.h
index fe72492e6..5c4bae0c5 100644
--- a/include/cantera/thermo/IdealMolalSoln.h
+++ b/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;
 
diff --git a/include/cantera/thermo/IdealSolidSolnPhase.h b/include/cantera/thermo/IdealSolidSolnPhase.h
index 7de203f6c..c51f42e77 100644
--- a/include/cantera/thermo/IdealSolidSolnPhase.h
+++ b/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
diff --git a/include/cantera/thermo/IdealSolnGasVPSS.h b/include/cantera/thermo/IdealSolnGasVPSS.h
index e9336464e..66b50051b 100644
--- a/include/cantera/thermo/IdealSolnGasVPSS.h
+++ b/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
diff --git a/include/cantera/thermo/LatticePhase.h b/include/cantera/thermo/LatticePhase.h
index 69a0ed604..ae4943d20 100644
--- a/include/cantera/thermo/LatticePhase.h
+++ b/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
diff --git a/include/cantera/thermo/LatticeSolidPhase.h b/include/cantera/thermo/LatticeSolidPhase.h
index 11250940f..1349e9542 100644
--- a/include/cantera/thermo/LatticeSolidPhase.h
+++ b/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;
diff --git a/include/cantera/thermo/MaskellSolidSolnPhase.h b/include/cantera/thermo/MaskellSolidSolnPhase.h
index 4935e1a5a..0fd9f7e09 100644
--- a/include/cantera/thermo/MaskellSolidSolnPhase.h
+++ b/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; }
diff --git a/include/cantera/thermo/MetalPhase.h b/include/cantera/thermo/MetalPhase.h
index 25b618a94..a527927e2 100644
--- a/include/cantera/thermo/MetalPhase.h
+++ b/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;
     }
diff --git a/include/cantera/thermo/PureFluidPhase.h b/include/cantera/thermo/PureFluidPhase.h
index a0de793aa..b5e779c25 100644
--- a/include/cantera/thermo/PureFluidPhase.h
+++ b/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;
 
diff --git a/include/cantera/thermo/StoichSubstance.h b/include/cantera/thermo/StoichSubstance.h
index b2f36ecfd..d9bf2ea54 100644
--- a/include/cantera/thermo/StoichSubstance.h
+++ b/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
diff --git a/include/cantera/thermo/ThermoPhase.h b/include/cantera/thermo/ThermoPhase.h
index 318a9be00..a1b6da7a6 100644
--- a/include/cantera/thermo/ThermoPhase.h
+++ b/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
diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp
index 7ffcdb06b..34455155e 100644
--- a/src/kinetics/Reaction.cpp
+++ b/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,
diff --git a/src/thermo/FixedChemPotSSTP.cpp b/src/thermo/FixedChemPotSSTP.cpp
index 30eac2de6..57ceb3e16 100644
--- a/src/thermo/FixedChemPotSSTP.cpp
+++ b/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;
diff --git a/src/thermo/GibbsExcessVPSSTP.cpp b/src/thermo/GibbsExcessVPSSTP.cpp
index 03b61b412..fc0bf486e 100644
--- a/src/thermo/GibbsExcessVPSSTP.cpp
+++ b/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);
diff --git a/src/thermo/IdealMolalSoln.cpp b/src/thermo/IdealMolalSoln.cpp
index 79c51e188..1167d5e87 100644
--- a/src/thermo/IdealMolalSoln.cpp
+++ b/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) {
diff --git a/src/thermo/IdealSolidSolnPhase.cpp b/src/thermo/IdealSolidSolnPhase.cpp
index 07e10ab96..74a4519a5 100644
--- a/src/thermo/IdealSolidSolnPhase.cpp
+++ b/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();
diff --git a/src/thermo/IdealSolnGasVPSS.cpp b/src/thermo/IdealSolnGasVPSS.cpp
index fead3b757..8997fcb8a 100644
--- a/src/thermo/IdealSolnGasVPSS.cpp
+++ b/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) {
diff --git a/src/thermo/LatticePhase.cpp b/src/thermo/LatticePhase.cpp
index b73d0fe41..acd838f73 100644
--- a/src/thermo/LatticePhase.cpp
+++ b/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);
diff --git a/src/thermo/LatticeSolidPhase.cpp b/src/thermo/LatticeSolidPhase.cpp
index 3720aeb1d..1277e73b3 100644
--- a/src/thermo/LatticeSolidPhase.cpp
+++ b/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();
diff --git a/src/thermo/PureFluidPhase.cpp b/src/thermo/PureFluidPhase.cpp
index 4d7cb31d4..709d99264 100644
--- a/src/thermo/PureFluidPhase.cpp
+++ b/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;
diff --git a/src/thermo/StoichSubstance.cpp b/src/thermo/StoichSubstance.cpp
index 88b8e94f5..a66316dbb 100644
--- a/src/thermo/StoichSubstance.cpp
+++ b/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;
diff --git a/src/thermo/ThermoPhase.cpp b/src/thermo/ThermoPhase.cpp
index 3b6a6cb7a..08d53800d 100644
--- a/src/thermo/ThermoPhase.cpp
+++ b/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));