Merge pull request NREL#193 from softwareengineerprogrammer/main

PR 185 Cleanup

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 3.4.25
+current_version = 3.4.26
 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.25
+    version: 3.4.26
     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/softwareengineerprogrammer/GEOPHIRES-X/v3.4.25.svg
+.. |commits-since| image:: https://img.shields.io/github/commits-since/softwareengineerprogrammer/GEOPHIRES-X/v3.4.26.svg
     :alt: Commits since latest release
-    :target: https://github.com/softwareengineerprogrammer/GEOPHIRES-X/compare/v3.4.25...main
+    :target: https://github.com/softwareengineerprogrammer/GEOPHIRES-X/compare/v3.4.26...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.25'
+version = release = '3.4.26'
 
 pygments_style = 'trac'
 templates_path = ['./templates']

setup.py

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

src/geophires_monte_carlo/MC_GeoPHIRES3.py

@@ -358,6 +358,7 @@ def main(command_line_args=None):
         if 'MC_OUTPUT_FILE' in args and args.MC_OUTPUT_FILE is not None
         else str(Path(Path(args.Input_file).parent, 'MC_Result.txt').absolute())
     )
+    code_file_name = Path(args.Code_File).name
     python_path = 'python'
     html_path = ''
 
@@ -450,6 +451,14 @@ def main(command_line_args=None):
     results_pd = pd.read_csv(output_file)
     df = pd.DataFrame(results_pd)
 
+    if len(results) < 1:
+        # TODO surface actual exceptions instead of giving this generic message
+        raise RuntimeError(
+            'No MC results generated, '
+            f'this is likely caused by {code_file_name} throwing an exception '
+            f'when run with your input file.'
+        )
+
     # Compute the stats along the specified axes.
     mins = np.nanmin(results, 0)
     maxs = np.nanmax(results, 0)

src/geophires_x/AGSWellBores.py

@@ -871,12 +871,12 @@ def CalculateNonverticalPressureDrop(self, model:Model, time_operation: float, t
             # nonvertical wellbore fluid conditions based on current temperature
             rhowater = density_water_kg_per_m3(
                 self.NonverticalProducedTemperature.value[year],
-                pressure=model.reserv.lithostatic_pressure(model.reserv.rhorock.value, model.reserv.InputDepth.quantity().to('m').magnitude)
+                pressure=model.reserv.lithostatic_pressure()
             )
 
             muwater = viscosity_water_Pa_sec(
                 self.NonverticalProducedTemperature.value[year],
-                pressure=model.reserv.lithostatic_pressure(model.reserv.rhorock.value, model.reserv.InputDepth.quantity().to('m').magnitude)
+                pressure=model.reserv.lithostatic_pressure()
             )
             vhoriz = self.q_circulation / rhowater / (math.pi / 4. * self.nonverticalwellborediameter.value ** 2)
 
@@ -957,15 +957,15 @@ def Calculate(self, model: Model) -> None:
                 # MIR figure out how to calculate year and extract Tini from reserv Tresoutput array
                 year = math.trunc(self.time_operation.value / self.al)
                 self.NonverticalProducedTemperature.value[year] = inverselaplace(
-                    self, 16, 0, model.reserv.lithostatic_pressure(model.reserv.rhorock.value, model.reserv.InputDepth.quantity().to('m').magnitude))
+                    self, 16, 0, model.reserv.lithostatic_pressure())
                 # update alpha_fluid value based on next temperature of reservoir
 
                 self.alpha_fluid = self.WaterThermalConductivity.value / density_water_kg_per_m3(
                     self.NonverticalProducedTemperature.value[year],
-                    pressure=model.reserv.lithostatic_pressure(model.reserv.rhorock.value, model.reserv.InputDepth.quantity().to('m').magnitude)
+                    pressure=model.reserv.lithostatic_pressure()
                 ) / heat_capacity_water_J_per_kg_per_K(
                     self.NonverticalProducedTemperature.value[year],
-                    pressure=model.reserv.lithostatic_pressure(model.reserv.rhorock.value, model.reserv.InputDepth.quantity().to('m').magnitude)
+                    pressure=model.reserv.lithostatic_pressure()
                 ) * 24.0 * 3600.0
                 self.time_operation.value += self.al
 
@@ -979,7 +979,7 @@ def Calculate(self, model: Model) -> None:
             self.ProdTempDrop.value = self.tempdropprod.value
             model.reserv.cpwater.value = heat_capacity_water_J_per_kg_per_K(
                 self.NonverticalProducedTemperature.value[0],
-                pressure=model.reserv.lithostatic_pressure(model.reserv.rhorock.value, model.reserv.InputDepth.quantity().to('m').magnitude)
+                pressure=model.reserv.lithostatic_pressure()
             )
             if self.rameyoptionprod.value:
                 self.ProdTempDrop.value = RameyCalc(model.reserv.krock.value,
@@ -1002,13 +1002,13 @@ def Calculate(self, model: Model) -> None:
             if self.productionwellpumping.value:
                 self.rhowaterinj = density_water_kg_per_m3(
                     model.reserv.Tsurf.value,
-                    pressure=model.reserv.lithostatic_pressure(model.reserv.rhorock.value, model.reserv.InputDepth.quantity().to('m').magnitude)
+                    pressure=model.reserv.lithostatic_pressure()
                 ) * np.linspace(1, 1,
                                 len(self.ProducedTemperature.value))
 
                 self.rhowaterprod = density_water_kg_per_m3(
                     model.reserv.Trock.value,
-                    pressure=model.reserv.lithostatic_pressure(model.reserv.rhorock.value, model.reserv.InputDepth.quantity().to('m').magnitude)
+                    pressure=model.reserv.lithostatic_pressure()
                 ) * np.linspace(1, 1, len(self.ProducedTemperature.value))
 
                 self.DPProdWell.value, f3, vprod, self.rhowaterprod = WellPressureDrop(model,
@@ -1113,13 +1113,13 @@ def Calculate(self, model: Model) -> None:
 
             rho_water = density_water_kg_per_m3(
                 self.Tout[0],
-                pressure = model.reserv.lithostatic_pressure(model.reserv.rhorock.value, model.reserv.InputDepth.quantity().to('m').magnitude),
+                pressure = model.reserv.lithostatic_pressure(),
             )
 
 
             model.reserv.cpwater.value = heat_capacity_water_J_per_kg_per_K(
                 self.Tout[0],
-                pressure=model.reserv.lithostatic_pressure(model.reserv.rhorock.value, model.reserv.InputDepth.quantity().to('m').magnitude),
+                pressure=model.reserv.lithostatic_pressure(),
             )  # 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

@@ -98,7 +98,7 @@ def __init__(self, model: Model):
             UnitType=Units.PERCENT,
             PreferredUnits=PercentUnit.TENTH,
             CurrentUnits=PercentUnit.TENTH,
-            ErrMessage="assume default cyclindrical reservoir radius of effect reduction factor (0.1)",
+            ErrMessage="assume default cylindrical reservoir radius of effect reduction factor (0.1)",
             ToolTipText="The radius of effect reduction factor - to account for the fact that we cannot extract 100%"
                         + " of the heat in the cylinder.",
         )
@@ -244,20 +244,17 @@ def Calculate(self, model: Model) -> None:
                                                  ) / 1e15  # 10^15 J
         self.cpwater.value = heat_capacity_water_J_per_kg_per_K(
             model.wellbores.Tinj.value * 0.5 + (self.Trock.value * 0.9 + model.wellbores.Tinj.value * 0.1) * 0.5,
-            pressure=self.lithostatic_pressure(model.reserv.rhorock.value, model.reserv.depth.quantity().to('m').magnitude)
+            pressure=self.lithostatic_pressure()
         )
         self.rhowater.value = density_water_kg_per_m3(
             model.wellbores.Tinj.value * 0.5 + (self.Trock.value * 0.9 + model.wellbores.Tinj.value * 0.1) * 0.5,
-            pressure=self.lithostatic_pressure(model.reserv.rhorock.value, model.reserv.depth.quantity().to('m').magnitude)
+            pressure=self.lithostatic_pressure()
         )
 
         model.logger.info(f'complete {str(__class__)}: {sys._getframe().f_code.co_name}')
 
-    #def lithostatic_pressure(self) -> PlainQuantity:
-        """
-        @override
+    def lithostatic_pressure(self) -> PlainQuantity:
+        """@override"""
 
-        Standard reservoir implementation uses depth but CylindricalReservoir sets depth to total drilled length
-        """
-        def lithostatic_pressure(self, rho_rock_kg_per_m3: float, depth_m: float) -> PlainQuantity:
-            return quantity(static_pressure_MPa(rho_rock_kg_per_m3, depth_m), 'MPa')
+        return quantity(static_pressure_MPa(self.rhorock.quantity().to('kg/m**3').magnitude,
+                                            self.InputDepth.quantity().to('m').magnitude), 'MPa')

src/geophires_x/Reservoir.py

@@ -749,14 +749,12 @@ def Calculate(self, model: Model) -> None:
             # calculate reservoir water properties
             self.cpwater.value = heat_capacity_water_J_per_kg_per_K(
                 model.wellbores.Tinj.value * 0.5 + (self.Trock.value * 0.9 + model.wellbores.Tinj.value * 0.1) * 0.5,
-                pressure=self.lithostatic_pressure(self.rhorock.quantity().to('kg/m**3').magnitude,
-                                            self.depth.quantity().to('m').magnitude)
+                pressure=self.lithostatic_pressure()
             )
 
             self.rhowater.value = density_water_kg_per_m3(
                 model.wellbores.Tinj.value * 0.5 + (self.Trock.value * 0.9 + model.wellbores.Tinj.value * 0.1) * 0.5,
-                pressure=self.lithostatic_pressure(self.rhorock.quantity().to('kg/m**3').magnitude,
-                                            self.depth.quantity().to('m').magnitude)
+                pressure=self.lithostatic_pressure()
             )
 
             # temperature gain in injection wells
@@ -768,7 +766,8 @@ def Calculate(self, model: Model) -> None:
 
         model.logger.info(f'complete {str(__class__)}: {sys._getframe().f_code.co_name}')
 
-    def lithostatic_pressure(self, rho_rock_kg_per_m3: float, depth_m: float) -> PlainQuantity:
-        return quantity(static_pressure_MPa(rho_rock_kg_per_m3, depth_m), 'MPa')
+    def lithostatic_pressure(self) -> PlainQuantity:
+        return quantity(static_pressure_MPa(self.rhorock.quantity().to('kg/m**3').magnitude,
+                                            self.depth.quantity().to('m').magnitude), 'MPa')
 
 

src/geophires_x/WellBores.py

@@ -164,15 +164,15 @@ def WellPressureDrop(model: Model, Taverage: float, wellflowrate: float, welldia
     rhowater = np.array([
         density_water_kg_per_m3(
             t,
-            pressure=model.reserv.lithostatic_pressure(model.reserv.rhorock.value, model.reserv.depth.quantity().to('m').magnitude),
+            pressure=model.reserv.lithostatic_pressure(),
         )
         for t in Taverage
     ])  # replace with correlation based on Tprodaverage
 
     muwater = np.array([
         viscosity_water_Pa_sec(
             t,
-            pressure=model.reserv.lithostatic_pressure(model.reserv.rhorock.value, model.reserv.depth.quantity().to('m').magnitude),
+            pressure=model.reserv.lithostatic_pressure(),
         )
         for t in Taverage
     ])  # replace with correlation based on Tprodaverage
@@ -228,13 +228,11 @@ def InjectionWellPressureDrop(model: Model, Taverage: float, wellflowrate: float
     """
     # start by calculating wellbore fluid conditions [kPa], noting that most temperature drop happens in
     # upper section (because surrounding rock temperature is lowest in upper section)
-    rhowater = (density_water_kg_per_m3(Taverage, pressure=model.reserv.lithostatic_pressure(model.reserv.rhorock.value,
-                                                                   model.reserv.depth.quantity().to('m').magnitude))
+    rhowater = (density_water_kg_per_m3(Taverage, pressure=model.reserv.lithostatic_pressure())
                 * np.linspace(1, 1, len(model.wellbores.ProducedTemperature.value)))
 
     # replace with correlation based on Tinjaverage
-    muwater = viscosity_water_Pa_sec(Taverage, pressure=model.reserv.lithostatic_pressure(model.reserv.rhorock.value,
-                                                                   model.reserv.depth.quantity().to('m').magnitude)) * np.linspace(1, 1, len(model.wellbores.ProducedTemperature.value))
+    muwater = viscosity_water_Pa_sec(Taverage, pressure=model.reserv.lithostatic_pressure()) * np.linspace(1, 1, len(model.wellbores.ProducedTemperature.value))
     v = nprod / ninj * wellflowrate * (1.0 + waterloss) / rhowater / (math.pi / 4. * welldiam ** 2)
     Rewater = 4. * nprod / ninj * wellflowrate * (1.0 + waterloss) / (
         muwater * math.pi * welldiam)  # laminar or turbulent flow?
@@ -1087,8 +1085,7 @@ def Calculate(self, model: Model) -> None:
         self.production_reservoir_pressure.value = get_hydrostatic_pressure_kPa(model.reserv.Trock.value, model.reserv.Tsurf.value,
                                                                                 model.reserv.depth.quantity().to('m').magnitude,
                                                                                 model.reserv.averagegradient.value,
-                                                                                model.reserv.lithostatic_pressure(model.reserv.rhorock.value,
-                                                                   model.reserv.depth.quantity().to('m').magnitude)) if self.usebuiltinhydrostaticpressurecorrelation else self.Phydrostatic.quantity().to(
+                                                                                model.reserv.lithostatic_pressure()) if self.usebuiltinhydrostaticpressurecorrelation else self.Phydrostatic.quantity().to(
             self.production_reservoir_pressure.CurrentUnits).magnitude
 
         self.production_reservoir_pressure.value = ReservoirPressurePredictor(model.surfaceplant.plant_lifetime.value,
@@ -1114,19 +1111,20 @@ def Calculate(self, model: Model) -> None:
                 if not self.injection_reservoir_temperature.Provided:
                     self.injection_reservoir_temperature.value = (model.reserv.averagegradient.value * self.injection_reservoir_depth.value) + model.reserv.Tsurf.value
 
+                injection_reservoir_static_pressure = quantity(static_pressure_MPa(
+                    model.reserv.rhorock.value, self.injection_reservoir_depth.quantity().to('m').magnitude), 'MPa')
+
                 if self.injection_reservoir_pressure.value < 0: # they didn't provide a pressure so assume hydrostatic.
                     self.injection_reservoir_pressure.value = get_hydrostatic_pressure_kPa(self.injection_reservoir_temperature.value,
                                                                                         model.reserv.Tsurf.value,
                                                                                         self.injection_reservoir_depth.value,
                                                                                         model.reserv.averagegradient.value * 1000.0,
-                                                                                        model.reserv.lithostatic_pressure(model.reserv.rhorock.value,
-                                                                                                               self.injection_reservoir_depth.quantity().to('m').magnitude))
+                                                                                        injection_reservoir_static_pressure)
             self.injection_reservoir_initial_pressure.value = self.injection_reservoir_pressure.value = get_hydrostatic_pressure_kPa(self.injection_reservoir_temperature.value,
                                                                                    model.reserv.Tsurf.value,
                                                                                    self.injection_reservoir_depth.value,
                                                                                    model.reserv.averagegradient.value,
-                                                                                   model.reserv.lithostatic_pressure(model.reserv.rhorock.value,
-                                                                                                                     self.injection_reservoir_depth.quantity().to('m').magnitude))
+                                                                                   injection_reservoir_static_pressure)
 
 #            if not self.injection_reservoir_initial_pressure.Provided:
             self.injection_reservoir_pressure.value = InjectionReservoirPressurePredictor(model.surfaceplant.plant_lifetime.value,

src/geophires_x/__init__.py

@@ -1 +1 @@
-__version__ = '3.4.25'
+__version__ = '3.4.26'