Merge pull request #11 from softwareengineerprogrammer/remove-water-p…

…roperty-fallback-calcs

Remove water property fallback calculations #110

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 3.4.1
+current_version = 3.4.2
 commit = True
 tag = True
 

.cookiecutterrc

@@ -54,7 +54,7 @@ default_context:
     sphinx_doctest: "no"
     sphinx_theme: "sphinx-py3doc-enhanced-theme"
     test_matrix_separate_coverage: "no"
-    version: 3.4.1
+    version: 3.4.2
     version_manager: "bump2version"
     website: "https://github.com/NREL"
     year_from: "2023"

README.rst

@@ -47,9 +47,9 @@ Free software: `MIT license <LICENSE>`__
     :alt: Supported implementations
     :target: https://pypi.org/project/geophires-x
 
-.. |commits-since| image:: https://img.shields.io/github/commits-since/NREL/GEOPHIRES-X/v3.4.1.svg
+.. |commits-since| image:: https://img.shields.io/github/commits-since/NREL/GEOPHIRES-X/v3.4.2.svg
     :alt: Commits since latest release
-    :target: https://github.com/NREL/GEOPHIRES-X/compare/v3.4.1...main
+    :target: https://github.com/NREL/GEOPHIRES-X/compare/v3.4.2...main
 
 .. |docs| image:: https://readthedocs.org/projects/GEOPHIRES-X/badge/?style=flat
     :target: https://nrel.github.io/GEOPHIRES-X

docs/conf.py

@@ -18,7 +18,7 @@
 year = '2023'
 author = 'NREL'
 copyright = f'{year}, {author}'
-version = release = '3.4.1'
+version = release = '3.4.2'
 
 pygments_style = 'trac'
 templates_path = ['./templates']

setup.py

@@ -13,7 +13,7 @@ def read(*names, **kwargs):
 
 setup(
     name='geophires-x',
-    version='3.4.1',
+    version='3.4.2',
     license='MIT',
     description='GEOPHIRES is a free and open-source geothermal techno-economic simulator.',
     long_description='{}\n{}'.format(

src/geophires_x/AGSWellBores.py

@@ -1044,19 +1044,19 @@ def Calculate(self, model: Model) -> None:
                 else:  # PI is used for both the verticals
                     UpgoingPumpingPower, self.PumpingPowerProd.value, self.DPProdWell.value, self.Pprodwellhead.value = \
                         ProdPressureDropAndPumpingPowerUsingIndexes(
-                            model, self.usebuiltinhydrostaticpressurecorrelation, self.productionwellpumping.value,
+                            model, self.productionwellpumping.value,
                             self.usebuiltinppwellheadcorrelation,
-                            model.reserv.Trock.value, model.reserv.Tsurf.value, model.reserv.depth.value,
-                            model.reserv.averagegradient.value, self.ppwellhead.value, self.PI.value,
+                            model.reserv.Trock.value, model.reserv.depth.value,
+                            self.ppwellhead.value, self.PI.value,
                             self.prodwellflowrate.value, f3, vprod,
                             self.prodwelldiam.value, self.nprod.value, model.surfaceplant.pump_efficiency.value,
                             self.rhowaterprod)
 
                     DowngoingPumpingPower, ppp2, dppw, ppwh = ProdPressureDropAndPumpingPowerUsingIndexes(
-                        model, self.usebuiltinhydrostaticpressurecorrelation, self.productionwellpumping.value,
+                        model, self.productionwellpumping.value,
                         self.usebuiltinppwellheadcorrelation,
-                        model.reserv.Trock.value, model.reserv.Tsurf.value, model.reserv.depth.value,
-                        model.reserv.averagegradient.value, self.ppwellhead.value, self.PI.value,
+                        model.reserv.Trock.value, model.reserv.depth.value,
+                        self.ppwellhead.value, self.PI.value,
                         self.prodwellflowrate.value, f3, vprod,
                         self.injwelldiam.value, self.nprod.value, model.surfaceplant.pump_efficiency.value,
                         self.rhowaterinj)
@@ -1102,7 +1102,10 @@ def Calculate(self, model: Model) -> None:
             self.ProducedTemperature.value = self.InterpolatedTemperatureArray.copy()
 
             tot_length, vert_length, horizontal_lengths = self.calculatedrillinglengths(model)
-            model.reserv.depth.value = model.reserv.InputDepth.value * 1000.0  # in this case, reserv.depth is just the vertical drill depth
+
+            # in this case, reserv.depth is just the vertical drill depth
+            # FIXME earlier calculations use depth before this value is set, meaning they're using the wrong value
+            model.reserv.depth.value = model.reserv.InputDepth.quantity().to(model.reserv.depth.CurrentUnits).magnitude
 
             # getTandP results must be rejiggered to match wellbores expected output. Once done,
             # the surfaceplant and economics models should just work
@@ -1116,19 +1119,12 @@ def Calculate(self, model: Model) -> None:
             rho_water = density_water_kg_per_m3(
                 self.Tout[0],
                 pressure = model.reserv.lithostatic_pressure(),
-
-                # FIXME TODO - get rid of fallback calculations https://github.com/NREL/GEOPHIRES-X/issues/110
-                enable_fallback_calculation=True
             )
 
 
             model.reserv.cpwater.value = heat_capacity_water_J_per_kg_per_K(
                 self.Tout[0],
-
                 pressure=model.reserv.lithostatic_pressure(),
-
-                # FIXME TODO - get rid of fallback calculations https://github.com/NREL/GEOPHIRES-X/issues/110
-                enable_fallback_calculation=True
             )  # Need this for surface plant output calculation
 
             # set pumping power to zero for all times, assuming that the thermosphere wil always

src/geophires_x/CylindricalReservoir.py

@@ -4,8 +4,9 @@
 from functools import lru_cache
 
 import numpy as np
+from pint.facets.plain import PlainQuantity
 
-from geophires_x.GeoPHIRESUtils import density_water_kg_per_m3
+from geophires_x.GeoPHIRESUtils import density_water_kg_per_m3, lithostatic_pressure_MPa, quantity
 
 from geophires_x.GeoPHIRESUtils import heat_capacity_water_J_per_kg_per_K
 import geophires_x.Model as Model
@@ -228,7 +229,9 @@ def Calculate(self, model: Model) -> None:
         # initialize with the Initial reservoir temperature
         self.Tresoutput.value = np.array(len(self.timevector.value) * [self.Trock.value])
         # depth in this case is actually the total length of the drilled assembly
-        self.depth.value = self.InputDepth.value / 1000.0 + self.OutputDepth.value + self.Length.value
+        self.depth.value = ((self.InputDepth.quantity() + self.OutputDepth.quantity() + self.Length.quantity()
+             ).to(self.depth.CurrentUnits).magnitude)
+
         # Total volume of all laterals but hollow cylinder - doesn't include drilled-out area, units = m3
         self.resvolcalc.value = (
             model.wellbores.numnonverticalsections.value
@@ -256,3 +259,12 @@ def Calculate(self, model: Model) -> None:
         )
 
         model.logger.info(f'complete {str(__class__)}: {sys._getframe().f_code.co_name}')
+
+    def lithostatic_pressure(self) -> PlainQuantity:
+        """
+        @override
+
+        Standard reservoir implementation uses depth but CylindricalReservoir sets depth to total drilled length
+        """
+        return quantity(lithostatic_pressure_MPa(self.rhorock.quantity().to('kg/m**3').magnitude,
+                                                 self.InputDepth.quantity().to('m').magnitude), 'MPa')

src/geophires_x/GeoPHIRESUtils.py

@@ -104,15 +104,13 @@ def quantity(value: float, unit: str) -> PlainQuantity:
 @lru_cache
 def density_water_kg_per_m3(
     Twater_degC: float,
-    pressure: Optional[PlainQuantity] = None,
-    enable_fallback_calculation=False) -> float:
+    pressure: Optional[PlainQuantity] = None) -> float:
     """
     Calculate the density of water as a function of temperature.
 
     Args:
         Twater_degC: The temperature of water in degrees C.
         pressure: Pressure - should be provided
-        enable_fallback_calculation: deprecation shim, do not use, see https://github.com/NREL/GEOPHIRES-X/issues/110
     Returns:
         The density of water in kg/m³.
     Raises:
@@ -129,15 +127,8 @@ def density_water_kg_per_m3(
             return CP.PropsSI('D', 'T', celsius_to_kelvin(Twater_degC), 'Q', 0, 'Water')
 
     except (NotImplementedError, ValueError) as e:
-        if enable_fallback_calculation:
-            _logger.warning(f'density_water: Fallback calculation triggered for {Twater_degC}C')
-
-            Twater_K = celsius_to_kelvin(Twater_degC)
-            # water density correlation as used in Geophires v1.2 [kg/m3]
-            rho_water = (.7983223 + (1.50896E-3 - 2.9104E-6 * Twater_K) * Twater_K) * 1E3
-            return rho_water
-
-        raise ValueError(f'Input temperature {Twater_degC} is out of range or otherwise not implemented') from e
+        raise ValueError(f'Input temperature & pressure ({Twater_degC}, {pressure}) '
+                         f'are out of range or otherwise could not be used to calculate') from e
 
 
 def celsius_to_kelvin(celsius: float) -> float:
@@ -161,12 +152,13 @@ def celsius_to_kelvin(celsius: float) -> float:
 @lru_cache
 def viscosity_water_Pa_sec(
         Twater_degC: float,
-        pressure: Optional[PlainQuantity] = None,
-        enable_fallback_calculation=False) -> float:
+    pressure: Optional[PlainQuantity] = None) -> float:
     """
-    The ViscosityWater function is used to calculate the dynamic viscosity of water as a function of temperature.
+    Calculate the dynamic viscosity of water as a function of temperature and pressure.
+
     Args:
         Twater_degC: the temperature of water in degrees C
+        pressure: Pressure - should be provided
     Returns:
         Viscosity of water in Pa·s (Ns/m2)
     Raises:
@@ -181,29 +173,21 @@ def viscosity_water_Pa_sec(
             return CP.PropsSI('V', 'T', celsius_to_kelvin(Twater_degC), 'Q', 0, 'Water')
 
     except (NotImplementedError, ValueError) as e:
-        if enable_fallback_calculation:
-            _logger.warning(f'viscosity_water: Fallback calculation triggered for {Twater_degC}C')
-
-            # accurate to within 2.5% from 0 to 370 degrees C [Ns/m2]
-            muwater = 2.414E-5 * np.power(10, 247.8 / (Twater_degC + 273.15 - 140))
-            return muwater
-
-        raise ValueError(f'Input temperature {Twater_degC} is out of range or otherwise not implemented') from e
+        raise ValueError(f'Input temperature & pressure ({Twater_degC}, {pressure}) '
+                         f'are out of range or otherwise could not be used to calculate') from e
 
 
 @lru_cache
 def heat_capacity_water_J_per_kg_per_K(
     Twater_degC: float,
     pressure: Optional[PlainQuantity] = None,
-    enable_fallback_calculation=False
 ) -> float:
     """
     Calculate the isobaric specific heat capacity (c_p) of water as a function of temperature.
 
     Args:
         Twater_degC: The temperature of water in degrees C.
         pressure: Pressure - should be provided
-        enable_fallback_calculation: deprecation shim, do not use, see https://github.com/NREL/GEOPHIRES-X/issues/110
     Returns:
         The isobaric specific heat capacity of water as a function of temperature in J/(kg·K).
     Raises:
@@ -224,21 +208,8 @@ def heat_capacity_water_J_per_kg_per_K(
             return CP.PropsSI('C', 'T', celsius_to_kelvin(Twater_degC), 'Q', 0, 'Water')
 
     except (NotImplementedError, ValueError) as e:
-        if enable_fallback_calculation:
-            _logger.warning(f'heat_capacity_water: Fallback calculation triggered for {Twater_degC}C')
-
-            Twater = (Twater_degC + 273.15) / 1000
-            A = -203.6060
-            B = 1523.290
-            C = -3196.413
-            D = 2474.455
-            E = 3.855326
-            # water specific heat capacity in J/(kg·K)
-            cpwater = (A + B * Twater + C * Twater ** 2 + D * Twater ** 3 + E / (Twater ** 2)) / 18.02 * 1000
-
-            return cpwater
-
-        raise ValueError(f'Input temperature {Twater_degC} is out of range or otherwise not implemented') from e
+        raise ValueError(f'Input temperature & pressure ({Twater_degC}, {pressure}) '
+                         f'are out of range or otherwise could not be used to calculate') from e
 
 
 @lru_cache
@@ -273,15 +244,15 @@ def RecoverableHeat(Twater_degC: float) -> float:
 
 
 @lru_cache
-def vapor_pressure_water_kPa(Twater_degC: float, enable_fallback_calculation=False) -> float:
+def vapor_pressure_water_kPa(
+    Twater_degC: float,
+    pressure: Optional[PlainQuantity] = None) -> float:
     """
     Calculate the vapor pressure of water as a function of temperature.
 
-    TODO accept pressure as parameter https://github.com/NREL/GEOPHIRES-X/issues/118
-
     Args:
         Twater_degC: the temperature of water in degrees C
-        enable_fallback_calculation: deprecation shim, do not use, see https://github.com/NREL/GEOPHIRES-X/issues/110
+        pressure: Pressure - should be provided
     Returns:
         The vapor pressure of water as a function of temperature in kPa
     Raises:
@@ -295,34 +266,19 @@ def vapor_pressure_water_kPa(Twater_degC: float, enable_fallback_calculation=Fal
         raise ValueError(f'Twater_degC ({Twater_degC}) must be greater than or equal to 0')
 
     try:
-        return (quantity(CP.PropsSI('P', 'T', celsius_to_kelvin(Twater_degC), 'Q', 0, 'Water'), 'Pa')
-                .to('kPa').magnitude)
-    except (NotImplementedError, ValueError) as e:
-        if enable_fallback_calculation:
-            _logger.warning(f'vapor_pressure_water: Fallback calculation triggered for {Twater_degC}C')
-            return _vapor_pressure_antoine_equation_kPa(Twater_degC)
+        if pressure is not None:
+            return (quantity(
+                CP.PropsSI('P', 'T', celsius_to_kelvin(Twater_degC), 'P', pressure.to('Pa').magnitude, 'Water'), 'Pa')
+                    .to('kPa').magnitude)
+        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('kPa').magnitude)
 
 
+    except (NotImplementedError, ValueError) as e:
         raise ValueError(f'Input temperature {Twater_degC} is out of range or otherwise not implemented') from e
 
-def _vapor_pressure_antoine_equation_kPa(Twater_degC: float) -> float:
-    """
-    water vapor pressure in kPa using Antione Equation
-    Do not add additional consumers, use geophires_x.GeoPHIRESUtils.vapor_pressure_water_kPa instead.
-    """
-
-    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
-    return vp
-
 
 @lru_cache
 def entropy_water_kJ_per_kg_per_K(temperature_degC: float) -> float:

src/geophires_x/Parameter.py

@@ -120,7 +120,7 @@ class boolParameter(Parameter):
         DefaultValue (bool, True):  The default value of that parameter
     """
 
-    value: bool = True
+    value: bool = True # FIXME set from DefaultValue
     DefaultValue: bool = value
     json_parameter_type: str = 'boolean'
 
@@ -137,7 +137,7 @@ class intParameter(Parameter):
         AllowableRange (list): A list of the valid values
     """
 
-    value: int = 0
+    value: int = 0 # FIXME set from DefaultValue
     DefaultValue: int = value
     AllowableRange: List[int] = field(default_factory=list)
     json_parameter_type: str = 'integer'
@@ -158,8 +158,13 @@ class floatParameter(Parameter):
                 which means that any value is valid by default
     """
 
-    value: float = 0.0
-    DefaultValue: float = value
+    def __post_init__(self):
+        if self.value is None:
+            self.value = self.DefaultValue
+
+    value: float = None
+
+    DefaultValue: float = 0.0
     Min: float = -1.8e30
     Max: float = 1.8e30
     json_parameter_type: str = 'number'
@@ -176,7 +181,7 @@ class strParameter(Parameter):
         DefaultValue (str, ""):  The default value of that parameter
     """
 
-    value: str = ""
+    value: str = "" # FIXME set from DefaultValue
     DefaultValue: str = value
     json_parameter_type: str = 'string'
 
@@ -196,7 +201,7 @@ class listParameter(Parameter):
             which means that any value is valid by default
     """
 
-    value: List[float] = field(default_factory=list)
+    value: List[float] = field(default_factory=list) # FIXME set from DefaultValue
     DefaultValue: List[float] = field(default_factory=list)
     Min: float = -1.8e308
     Max: float = 1.8e308