stash - try to compare various methods of calculating vapor pressure

NREL · Feb 13, 2024 · ef9226e · ef9226e
1 parent 2f35fa1
commit ef9226e
Show file tree

Hide file tree

Showing 2 changed files with 21 additions and 13 deletions.
diff --git a/src/geophires_x/GeoPHIRESUtils.py b/src/geophires_x/GeoPHIRESUtils.py
@@ -12,6 +12,7 @@
 
 import pint
 import scipy
+from iapws import IAPWS97
 from pint.facets.plain import PlainQuantity
 from scipy.interpolate import interp1d
 import numpy as np
@@ -298,8 +299,9 @@ def vapor_pressure_water_kPa(
 
     try:
         if pressure is not None:
+            p_iapws = quantity(IAPWS97(T=celsius_to_kelvin(Twater_degC), P=pressure.to('MPa').magnitude).P, 'MPa').to('kPa').magnitude
             return quantity(CP.PropsSI('P', 'T', celsius_to_kelvin(Twater_degC), 'P', pressure.to('Pa').magnitude,
-                                       'Water'), 'Pa').to('kPa').magnitude
+                                       'Water'), 'Pa').to('kPa').magnitude, p_iapws, _antoine_equation_kPa(Twater_degC)
         else:
             _logger.warning(f'vapor_pressure_water: No pressure provided, using vapor quality=0 instead')
             return quantity(CP.PropsSI('P', 'T', celsius_to_kelvin(Twater_degC), 'Q', 0, 'Water'), 'Pa').to(
@@ -309,21 +311,23 @@ def vapor_pressure_water_kPa(
     except (NotImplementedError, ValueError) as e:
         if enable_fallback_calculation:
             _logger.warning(f'vapor_pressure_water: Fallback calculation triggered for {Twater_degC}C')
+            return _antoine_equation_kPa(Twater_degC)
 
-            if Twater_degC < 100:
-                A = 8.07131
-                B = 1730.63
-                C = 233.426
-            else:
-                A = 8.14019
-                B = 1810.94
-                C = 244.485
-            vp = 133.322 * (
-                10 ** (A - B / (C + Twater_degC))) / 1000  # water vapor pressure in kPa using Antione Equation
-            return vp
 
         raise ValueError(f'Input temperature {Twater_degC} is out of range or otherwise not implemented') from e
 
+def _antoine_equation_kPa(Twater_degC: float):
+    if Twater_degC < 100:
+        A = 8.07131
+        B = 1730.63
+        C = 233.426
+    else:
+        A = 8.14019
+        B = 1810.94
+        C = 244.485
+    vp = 133.322 * (
+        10 ** (A - B / (C + Twater_degC))) / 1000  # water vapor pressure in kPa using Antione Equation
+    return vp
 
 @lru_cache
 def entropy_water_kJ_per_kg_per_K(temperature_degC: float) -> float:

diff --git a/tests/test_geophires_utils.py b/tests/test_geophires_utils.py
@@ -434,7 +434,8 @@ def test_below_100_degrees(self):
         result = vapor_pressure_water_kPa(42)
         self.assertAlmostEqual(result, 8.209563332516748, places=3)
 
-        result = vapor_pressure_water_kPa(42, pressure=quantity(100, 'MPa'))
+        # result = vapor_pressure_water_kPa(42, pressure=quantity(100, 'MPa'))
+        result = vapor_pressure_water_kPa(42, pressure=quantity(1, 'atm'))
         self.assertAlmostEqual(result, 8.209563332516748, places=3)
 
     def test_above_100_degrees(self):
@@ -453,6 +454,9 @@ def test_25_degrees(self):
         result = vapor_pressure_water_kPa(25)
         self.assertAlmostEqual(result, 3.1697468549523626, places=3)
 
+        p_cp, p_iapws, p_antoine = vapor_pressure_water_kPa(35, pressure=quantity(1, 'atm'))
+        self.assertAlmostEqual(p_cp, 8.209563332516748, places=3)
+
     def test_value_error(self):
         with self.assertRaises(ValueError):
             vapor_pressure_water_kPa('abc')