Make volume of SurfPhase consistent with h = u + P*v

Defines molar volume to be zero. Density remains defined as being
per unit area or length, depending on the dimensionality of the phase.

Also fixes setting surface composition using mass/mole fractions so
that the sum of concentrations always equals the site density.

Partially resolves #1312

include/cantera/thermo/LatticeSolidPhase.h

@@ -311,7 +311,7 @@ class LatticeSolidPhase : public ThermoPhase
         throw NotImplementedError("LatticeSolidPhase::getConcentrations");
     }
 
-    virtual doublereal concentration(int k) const {
+    virtual doublereal concentration(size_t k) const {
         throw NotImplementedError("LatticeSolidPhase::concentration");
     }
 

include/cantera/thermo/Phase.h

@@ -65,6 +65,12 @@ class Species;
  * An example of this is the function
  * Phase::setState_TRY(double t, double dens, const double* y).
  *
+ * For bulk (3-dimensional) phases, the mass density has units of kg/m^3, and the molar
+ * density and concentrations have units of kmol/m^3, and the units listed in the
+ * methods of the Phase class assume a bulk phase. However, for surface (2-dimensional)
+ * phases have units of kg/m^2 and kmol/m^2, respectively. And for edge (1-dimensional)
+ * phases, these units kg/m and kmol/m.
+ *
  * Class Phase contains methods for saving and restoring the full internal state
  * of a given phase. These are saveState() and restoreState(). These functions
  * operate on a state vector, which by default uses the first two entries for
@@ -532,7 +538,7 @@ class Phase
      *                  kmol/m^3. The length of the vector must be greater than
      *                  or equal to the number of species within the phase.
      */
-    void getConcentrations(double* const c) const;
+    virtual void getConcentrations(double* const c) const;
 
     //! Concentration of species k.
     //! If k is outside the valid range, an exception will be thrown.
@@ -541,7 +547,7 @@ class Phase
      *
      *    @returns the concentration of species k (kmol m-3).
      */
-    double concentration(const size_t k) const;
+    virtual double concentration(const size_t k) const;
 
     //! Set the concentrations to the specified values within the phase.
     //! We set the concentrations here and therefore we set the overall density
@@ -664,11 +670,11 @@ class Phase
 
     //! Molar density (kmol/m^3).
     //!     @return The molar density of the phase
-    double molarDensity() const;
+    virtual double molarDensity() const;
 
     //! Molar volume (m^3/kmol).
     //!     @return The molar volume of the phase
-    double molarVolume() const;
+    virtual double molarVolume() const;
 
     //! Set the internally stored density (kg/m^3) of the phase.
     //! Note the density of a phase is an independent variable.

include/cantera/thermo/SurfPhase.h

@@ -189,8 +189,8 @@ class SurfPhase : public ThermoPhase
      *
      * @param k Optional parameter indicating the species. The default
      *          is to assume this refers to species 0.
-     * @return
-     *   Returns the standard Concentration in units of m3 kmol-1.
+     * @return the standard concentration in units of kmol/m^2 for surface phases or
+     *     kmol/m for edge phases.
      */
     virtual doublereal standardConcentration(size_t k = 0) const;
     virtual doublereal logStandardConc(size_t k=0) const;
@@ -200,6 +200,20 @@ class SurfPhase : public ThermoPhase
 
     virtual bool addSpecies(shared_ptr<Species> spec);
 
+    //! Since interface phases have no volume, this returns 0.0.
+    virtual double molarVolume() const {
+        return 0.0;
+    }
+
+    //! Since interface phases have no volume, setting this to a value other than 0.0
+    //! raises an exception.
+    virtual void setMolarDensity(const double vm) {
+        if (vm != 0.0) {
+            throw CanteraError("SurfPhase::setMolarDensity",
+                               "The volume of an interface is zero");
+        }
+    }
+
     //! Returns the site density
     /*!
      * Site density kmol m-2
@@ -296,6 +310,8 @@ class SurfPhase : public ThermoPhase
     virtual void setState(const AnyMap& state);
 
 protected:
+    virtual void compositionChanged();
+
     //! Surface site density (kmol m-2)
     doublereal m_n0;
 

interfaces/cython/cantera/thermo.pyx

@@ -691,7 +691,10 @@ cdef class ThermoPhase(_SolutionBase):
                 self._setArray1(thermo_setMoleFractions, X)
 
     property concentrations:
-        """Get/Set the species concentrations [kmol/m^3]."""
+        """
+        Get/Set the species concentrations. Units are kmol/m^3 for bulk phases, kmol/m^2
+        for surface phases, and kmol/m for edge phases.
+        """
         def __get__(self):
             return self._getArray1(thermo_getConcentrations)
         def __set__(self, C):

src/thermo/SurfPhase.cpp

@@ -159,9 +159,8 @@ void SurfPhase::getCp_R(doublereal* cpr) const
 
 void SurfPhase::getStandardVolumes(doublereal* vol) const
 {
-    _updateThermo();
     for (size_t k = 0; k < m_kk; k++) {
-        vol[k] = 1.0/standardConcentration(k);
+        vol[k] = 0.0;
     }
 }
 
@@ -211,6 +210,7 @@ void SurfPhase::setSiteDensity(doublereal n0)
                            "Site density must be positive. Got {}", n0);
     }
     m_n0 = n0;
+    assignDensity(n0 * meanMolecularWeight());
     m_logn0 = log(m_n0);
 }
 
@@ -281,6 +281,12 @@ void SurfPhase::setState(const AnyMap& state) {
     ThermoPhase::setState(state);
 }
 
+void SurfPhase::compositionChanged()
+{
+    ThermoPhase::compositionChanged();
+    assignDensity(m_n0 * meanMolecularWeight());
+}
+
 void SurfPhase::_updateThermo(bool force) const
 {
     doublereal tnow = temperature();

test/data/consistency-cases.yaml

@@ -271,7 +271,6 @@ ideal-surface:
     phase: Pt-surf
     atol_v: 1e-7
     known-failures:
-      h_eq_u_plus_Pv: "Definition of volume is unclear. See GitHub Issue #1312"
       gk_eq_hk_minus_T_sk:
         "Implementation of s_k is incorrect. See GitHub Issue #1313"
       s_eq_sum_sk_Xk:
@@ -292,7 +291,6 @@ ideal-edge:
     phase: TPB
     atol_v: 1e3  # site density of 5e-18 kmol/m = linear molar volume of 2e17 m/kmol
     known-failures:
-      h_eq_u_plus_Pv: "Definition of volume is unclear. See GitHub Issue #1312"
       dSdv_const_T_eq_dPdT_const_V:
         "Compressibility of phase leads to inconsistent results. See GitHub Issue #1312"
   states:

test/python/test_thermo.py

@@ -1113,6 +1113,10 @@ def test_coverages_array(self):
         self.assertNear(C[4], 0.25)
         self.assertNear(sum(C), 1.0)
 
+    def test_mole_fractions(self):
+        self.interface.X = 'c6HM:0.3, c6H*:0.7'
+        self.assertNear(sum(self.interface.concentrations), self.interface.site_density)
+
     def test_coverages_string(self):
         self.interface.coverages = 'c6HM:0.2, c6H*:0.8'
         C = self.interface.coverages

test/python/utilities.py

@@ -94,6 +94,8 @@ def assertIsNaN(self, value):
             self.fail(f"Value '{value}' is a number")
 
     def assertNear(self, a, b, rtol=1e-8, atol=1e-12, msg=None):
+        if a == b:
+            return  # handles case where a == b == inf
         cmp = 2 * abs(a - b)/(abs(a) + abs(b) + 2 * atol / rtol)
         if not cmp < rtol:
             message = ('AssertNear: %.14g - %.14g = %.14g\n' % (a, b, a-b) +

test/thermo/thermoToYaml.cpp

@@ -357,8 +357,12 @@ class ThermoYamlRoundTrip : public testing::Test
             duplicate->setState_TPX(T, P, X);
         }
 
-        EXPECT_NEAR(original->density(), duplicate->density(),
-                    rtol * original->density());
+        double rhoOrig = original->density();
+        double rhoDup = duplicate->density();
+        if (rhoOrig != rhoDup) {
+            EXPECT_NEAR(original->density(), duplicate->density(),
+                        rtol * original->density());
+        }
         if (!skip_cp) {
             EXPECT_NEAR(original->cp_mass(), duplicate->cp_mass(),
                         rtol * original->cp_mass());

test_problems/diamondSurf/runDiamond.cpp

@@ -10,15 +10,6 @@
 using namespace std;
 using namespace Cantera;
 
-void printDbl(double val)
-{
-    if (fabs(val) < 5.0E-17) {
-        cout << " nil";
-    } else {
-        cout << val;
-    }
-}
-
 int main(int argc, char** argv)
 {
 #if defined(_MSC_VER) && _MSC_VER < 1900
@@ -96,15 +87,11 @@ int main(int argc, char** argv)
                 int itp = k - 5;
                 naH = diamond100->nAtoms(itp, 0);
             }
-            cout << k << "  " << naH << "  " ;
-            printDbl(src[k]);
-            cout << endl;
+            writelog("{} {} {}\n", k, naH, src[k]);
             sum += naH * src[k];
         }
 
-        cout << "sum = ";
-        printDbl(sum);
-        cout << endl;
+        writelog("sum = {}\n", sum);
         double mwd = diamond->molecularWeight(0);
         double dens = diamond->density();
         double gr = src[4] * mwd / dens;

test_problems/diamondSurf/runDiamond_blessed.out

@@ -5,17 +5,17 @@ Number of reactions = 20
 0  1  -8.96e-05
 1  2  4.486e-05
 2  3  -3.801e-08
-3  4   nil
+3  4   0.0
 4  0  3.801e-08
-5  2   nil
-6  1   nil
-7  1   nil
-8  0   nil
-9  4   nil
-10  3   nil
-11  3   nil
-12  2   nil
-sum =  nil
+5  2   0.0
+6  1   0.0
+7  1   0.0
+8  0   0.0
+9  4   0.0
+10  3   0.0
+11  3   0.0
+12  2   0.0
+sum =  0.0
 growth rate = 0.4669 microns per hour
 Coverages:
 0   c6HH   0.4622

test_problems/fortran/f90_demo_blessed.txt

@@ -82,27 +82,27 @@ Thermal conductivity:     0.16396      W/m/K
  ********   Interface Kinetics Test   ********
 
  Equation                                     k_fwd        Net rate
- H2 + 2 PT(S) => 2 H(S)                     0.33111E-01   0.33111E-01
+ H2 + 2 PT(S) => 2 H(S)                     0.33954E-01   0.33954E-01
  2 H(S) => H2 + 2 PT(S)                     0.00000E+00   0.00000E+00
  H + PT(S) => H(S)                          0.00000E+00   0.00000E+00
- O2 + 2 PT(S) => 2 O(S)                     0.15290E-01   0.15290E-01
- O2 + 2 PT(S) => 2 O(S)                     0.20585E-01   0.20585E-01
- 2 O(S) => O2 + 2 PT(S)                     0.14336E-07   0.14336E-07
+ O2 + 2 PT(S) => 2 O(S)                     0.16078E-01   0.16078E-01
+ O2 + 2 PT(S) => 2 O(S)                     0.21646E-01   0.21646E-01
+ 2 O(S) => O2 + 2 PT(S)                     0.15097E-07   0.15097E-07
  O + PT(S) => O(S)                          0.00000E+00   0.00000E+00
- H2O + PT(S) => H2O(S)                      0.35283E+00   0.35283E+00
+ H2O + PT(S) => H2O(S)                      0.36181E+00   0.36181E+00
  H2O(S) => H2O + PT(S)                      0.00000E+00   0.00000E+00
  OH + PT(S) => OH(S)                        0.00000E+00   0.00000E+00
  OH(S) => OH + PT(S)                        0.00000E+00   0.00000E+00
  H(S) + O(S) <=> OH(S) + PT(S)              0.00000E+00   0.00000E+00
  H(S) + OH(S) <=> H2O(S) + PT(S)            0.00000E+00   0.00000E+00
  2 OH(S) <=> H2O(S) + O(S)                  0.00000E+00   0.00000E+00
- CO + PT(S) => CO(S)                        0.12687E+00   0.12687E+00
- CO(S) => CO + PT(S)                        0.17276E+00   0.17276E+00
+ CO + PT(S) => CO(S)                        0.13341E+00   0.13341E+00
+ CO(S) => CO + PT(S)                        0.17716E+00   0.17716E+00
  CO2(S) => CO2 + PT(S)                      0.00000E+00   0.00000E+00
- CO(S) + O(S) => CO2(S) + PT(S)             0.13165E-01   0.13165E-01
+ CO(S) + O(S) => CO2(S) + PT(S)             0.13844E-01   0.13844E-01
  CH4 + 2 PT(S) => CH3(S) + H(S)             0.00000E+00   0.00000E+00
  CH3(S) + PT(S) => CH2(S)s + H(S)           0.00000E+00   0.00000E+00
  CH2(S)s + PT(S) => CH(S) + H(S)            0.00000E+00   0.00000E+00
  CH(S) + PT(S) => C(S) + H(S)               0.00000E+00   0.00000E+00
  C(S) + O(S) => CO(S) + PT(S)               0.00000E+00   0.00000E+00
- CO(S) + PT(S) => C(S) + O(S)               0.10366E-06   0.10366E-06
+ CO(S) + PT(S) => C(S) + O(S)               0.10900E-06   0.10900E-06