Make SurfPhase incompressible

Partially resolves #1312
Cantera · Aug 9, 2022 · eedcbf9 · eedcbf9
1 parent 43d2241
commit eedcbf9
Show file tree

Hide file tree

Showing 6 changed files with 27 additions and 14 deletions.
diff --git a/include/cantera/thermo/SurfPhase.h b/include/cantera/thermo/SurfPhase.h
@@ -110,6 +110,10 @@ class SurfPhase : public ThermoPhase
         return "Surf";
     }
 
+    virtual bool isCompressible() const {
+        return false;
+    }
+
     //! Return the Molar Enthalpy. Units: J/kmol.
     /*!
      * For an ideal solution,

diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py
@@ -563,7 +563,7 @@ def write_hdf(self, filename, *args, group=None, species='X', mode='a',
         `SolutionArray.write_hdf` via `to_solution_array` and requires a working
         installation of *h5py* (``h5py`` can be installed using pip or conda).
         """
-        cols = ('extra', 'T', 'D', species)
+        cols = ('extra', 'T', 'P', species)
         meta = self.settings
         meta['date'] = formatdate(localtime=True)
         meta['cantera_version'] = __version__

diff --git a/src/kinetics/solveSP.cpp b/src/kinetics/solveSP.cpp
@@ -295,9 +295,14 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin,
 
 void solveSP::updateState(const doublereal* CSolnSP)
 {
+    vector_fp X;
     size_t loc = 0;
     for (size_t n = 0; n < m_numSurfPhases; n++) {
-        m_ptrsSurfPhase[n]->setConcentrations(CSolnSP + loc);
+        X.resize(m_nSpeciesSurfPhase[n]);
+        for (size_t k = 0; k < X.size(); k++) {
+            X[k] = CSolnSP[loc + k] / m_ptrsSurfPhase[n]->siteDensity();
+        }
+        m_ptrsSurfPhase[n]->setMoleFractions_NoNorm(X.data());
         loc += m_nSpeciesSurfPhase[n];
     }
 }

diff --git a/src/thermo/SurfPhase.cpp b/src/thermo/SurfPhase.cpp
@@ -218,25 +218,34 @@ void SurfPhase::setCoverages(const doublereal* theta)
 {
     double sum = 0.0;
     for (size_t k = 0; k < m_kk; k++) {
-        sum += theta[k];
+        sum += theta[k] / size(k);
     }
     if (sum <= 0.0) {
         throw CanteraError("SurfPhase::setCoverages",
                            "Sum of Coverage fractions is zero or negative");
     }
     for (size_t k = 0; k < m_kk; k++) {
-        m_work[k] = m_n0*theta[k]/(sum*size(k));
+        m_work[k] = theta[k] / (sum * size(k));
     }
-    // Call the Phase:: class function setConcentrations.
-    setConcentrations(m_work.data());
+    setMoleFractions(m_work.data());
 }
 
 void SurfPhase::setCoveragesNoNorm(const doublereal* theta)
 {
+    double sum = 0.0;
+    double sum2 = 0.0;
+    for (size_t k = 0; k < m_kk; k++) {
+        sum += theta[k] / size(k);
+        sum2 += theta[k];
+    }
+    if (sum <= 0.0) {
+        throw CanteraError("SurfPhase::setCoverages",
+                           "Sum of Coverage fractions is zero or negative");
+    }
     for (size_t k = 0; k < m_kk; k++) {
-        m_work[k] = m_n0*theta[k]/size(k);
+        m_work[k] = theta[k] * sum2 / (sum * size(k));
     }
-    setConcentrationsNoNorm(m_work.data());
+    setMoleFractions_NoNorm(m_work.data());
 }
 
 void SurfPhase::getCoverages(doublereal* theta) const

diff --git a/test/data/consistency-cases.yaml b/test/data/consistency-cases.yaml
@@ -277,8 +277,6 @@ ideal-surface:
         "Implementation of s_k is incorrect. See GitHub Issue #1313"
       g_eq_sum_gk_Xk:
         "chemPotentials does not protect against inf. See GitHub Issue #1314"
-      dSdv_const_T_eq_dPdT_const_V:
-        "Compressibility of phase leads to inconsistent results. See GitHub Issue #1312"
   states:
   - {T: 800, P: 1 atm, coverages: {Pt(s): 0.5, H(s): 0.4, O(s): 0.1}}
   - {T: 800, P: 5 atm, coverages: {H(s): 1.0}}
@@ -290,9 +288,6 @@ ideal-edge:
     file: surface-phases.yaml
     phase: TPB
     atol_v: 1e3  # site density of 5e-18 kmol/m = linear molar volume of 2e17 m/kmol
-    known-failures:
-      dSdv_const_T_eq_dPdT_const_V:
-        "Compressibility of phase leads to inconsistent results. See GitHub Issue #1312"
   states:
   - {T: 300, P: 1 atm}
   - {T: 900, P: 20 atm}

diff --git a/test/python/test_kinetics.py b/test/python/test_kinetics.py
@@ -228,7 +228,7 @@ def test_surface(self):
                              adjacent=[gas])
         surf1.site_density = surf2.site_density = 5e-9
         gas.TP = surf2.TP = surf1.TP = 900, 2*ct.one_atm
-        surf2.concentrations = surf1.concentrations
+        surf2.coverages = surf1.coverages
 
         self.assertEqual(surf1.n_reactions, surf2.n_reactions)