From 43d224142302f2485d7e246da17c5cadac474064 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 29 Jul 2022 13:26:44 -0400 Subject: [PATCH] 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 | 2 +- include/cantera/thermo/Phase.h | 14 +++++++++---- include/cantera/thermo/SurfPhase.h | 20 +++++++++++++++++-- interfaces/cython/cantera/thermo.pyx | 5 ++++- src/thermo/SurfPhase.cpp | 10 ++++++++-- test/data/consistency-cases.yaml | 2 -- test/python/test_thermo.py | 4 ++++ test/python/utilities.py | 2 ++ test/thermo/thermoToYaml.cpp | 8 ++++++-- test_problems/diamondSurf/runDiamond.cpp | 17 ++-------------- .../diamondSurf/runDiamond_blessed.out | 20 +++++++++---------- test_problems/fortran/f90_demo_blessed.txt | 18 ++++++++--------- 12 files changed, 74 insertions(+), 48 deletions(-) diff --git a/include/cantera/thermo/LatticeSolidPhase.h b/include/cantera/thermo/LatticeSolidPhase.h index 0de31b4279..67a7ac1b2a 100644 --- a/include/cantera/thermo/LatticeSolidPhase.h +++ b/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"); } diff --git a/include/cantera/thermo/Phase.h b/include/cantera/thermo/Phase.h index e6d3e036c4..5d7c64431a 100644 --- a/include/cantera/thermo/Phase.h +++ b/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. diff --git a/include/cantera/thermo/SurfPhase.h b/include/cantera/thermo/SurfPhase.h index a270ef4acc..f191cff2df 100644 --- a/include/cantera/thermo/SurfPhase.h +++ b/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 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; diff --git a/interfaces/cython/cantera/thermo.pyx b/interfaces/cython/cantera/thermo.pyx index 386c1e5988..efb3a890f5 100644 --- a/interfaces/cython/cantera/thermo.pyx +++ b/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): diff --git a/src/thermo/SurfPhase.cpp b/src/thermo/SurfPhase.cpp index bd9fc6ca9d..40ce9be31a 100644 --- a/src/thermo/SurfPhase.cpp +++ b/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(); diff --git a/test/data/consistency-cases.yaml b/test/data/consistency-cases.yaml index 2b6e6f17c3..17024c9038 100644 --- a/test/data/consistency-cases.yaml +++ b/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: diff --git a/test/python/test_thermo.py b/test/python/test_thermo.py index 27fee1f418..0463fb857f 100644 --- a/test/python/test_thermo.py +++ b/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 diff --git a/test/python/utilities.py b/test/python/utilities.py index 6cddb0339e..ec184ff1c9 100644 --- a/test/python/utilities.py +++ b/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) + diff --git a/test/thermo/thermoToYaml.cpp b/test/thermo/thermoToYaml.cpp index 9cbf667e5b..dd158d8779 100644 --- a/test/thermo/thermoToYaml.cpp +++ 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()); diff --git a/test_problems/diamondSurf/runDiamond.cpp b/test_problems/diamondSurf/runDiamond.cpp index c9d75a8e7d..74bcceb258 100644 --- a/test_problems/diamondSurf/runDiamond.cpp +++ b/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; diff --git a/test_problems/diamondSurf/runDiamond_blessed.out b/test_problems/diamondSurf/runDiamond_blessed.out index c3df420d46..90577f52e8 100644 --- a/test_problems/diamondSurf/runDiamond_blessed.out +++ b/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 diff --git a/test_problems/fortran/f90_demo_blessed.txt b/test_problems/fortran/f90_demo_blessed.txt index 1f2428d3f3..945d5cfeee 100644 --- a/test_problems/fortran/f90_demo_blessed.txt +++ b/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