Gringarten-Stehfest Precision [v3.9.52]

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 3.9.50
+current_version = 3.9.52
 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.9.50
+    version: 3.9.52
     version_manager: "bump2version"
     website: "https://github.com/NREL"
     year_from: "2023"

.gitignore

@@ -104,6 +104,9 @@ _site/
 /Useful sites for Sphinx docstrings.txt
 /.github/workflows/workflows.7z
 
+# TODO unignore once favicon is correctly configured
+docs/_images/geophires-favicon.png
+
 # Mypy Cache
 .mypy_cache/
 

README.rst

@@ -58,9 +58,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.9.50.svg
+.. |commits-since| image:: https://img.shields.io/github/commits-since/softwareengineerprogrammer/GEOPHIRES-X/v3.9.52.svg
     :alt: Commits since latest release
-    :target: https://github.com/softwareengineerprogrammer/GEOPHIRES-X/compare/v3.9.50...main
+    :target: https://github.com/softwareengineerprogrammer/GEOPHIRES-X/compare/v3.9.52...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 = '2025'
 author = 'NREL'
 copyright = f'{year}, {author}'
-version = release = '3.9.50'
+version = release = '3.9.52'
 
 pygments_style = 'trac'
 templates_path = ['./templates']

setup.py

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

src/geophires_x/MPFReservoir.py

@@ -1,8 +1,12 @@
 import sys
+
+from mpmath import mp, exp, invertlaplace, sqrt, tanh, workdps
 import numpy as np
-from mpmath import *
+
 import geophires_x.Model as Model
+from .Parameter import intParameter
 from .Reservoir import Reservoir
+from .Units import Units
 
 
 class MPFReservoir(Reservoir):
@@ -12,6 +16,7 @@ class MPFReservoir(Reservoir):
     It also has its own methods and attributes that are unique to this class.
     """
 
+    # noinspection PyUnresolvedReferences,PyProtectedMember
     def __init__(self, model: Model):
         """
         The __init__ function is called automatically when a class is instantiated.
@@ -31,15 +36,33 @@ def __init__(self, model: Model):
         :type model: :class:`~geophires_x.Model.Model`
         :return: None
         """
-        model.logger.info("Init " + str(__class__) + ": " + sys._getframe().f_code.co_name)
+        model.logger.info(f'Init {str(__class__)}: {sys._getframe().f_code.co_name}')
+
         super().__init__(model)  # initialize the parent parameters and variables
         sclass = str(__class__).replace("<class \'", "")
         self.MyClass = sclass.replace("\'>", "")
-        model.logger.info("Complete " + str(__class__) + ": " + sys._getframe().f_code.co_name)
+
+        max_allowed_precision = 15
+        self.gringarten_stehfest_precision = self.ParameterDict[self.gringarten_stehfest_precision.Name] = intParameter(
+            'Gringarten-Stehfest Precision',
+            DefaultValue=15,
+            AllowableRange=list(range(8, max_allowed_precision + 1)),
+            UnitType=Units.NONE,
+            Required=False,
+            ToolTipText='Sets the numerical precision (decimal places) for the inverse Laplace transform '
+                        '(Stehfest algorithm) used in the Gringarten calculation for the Multiple Parallel Fractures '
+                        'Reservoir Model. '
+                        'The default value provides maximum result stability; '
+                        'lower values calculate faster but may reduce consistency.'
+        )
+
+
+        model.logger.info(f'Complete {str(__class__)}: {sys._getframe().f_code.co_name}')
 
     def __str__(self):
         return "MPFReservoir"
 
+    # noinspection PyUnresolvedReferences,PyProtectedMember
     def read_parameters(self, model: Model) -> None:
         """
         The read_parameters function reads in the parameters from a dictionary created by reading the user-provided file
@@ -51,15 +74,16 @@ def read_parameters(self, model: Model) -> None:
         :type model: :class:`~geophires_x.Model.Model`
         :return: None
         """
-        model.logger.info("Init " + str(__class__) + ": " + sys._getframe().f_code.co_name)
+        model.logger.info(f'Init {str(__class__)}: {sys._getframe().f_code.co_name}')
         # if we call super, we don't need to deal with setting the parameters here,
         # just deal with the special cases for the variables in this class
         # because the call to the super.readparameters will set all the variables,
         # including the ones that are specific to this class
         super().read_parameters(model)  # read the parameters for the parent.
 
-        model.logger.info("Complete " + str(__class__) + ": " + sys._getframe().f_code.co_name)
+        model.logger.info(f'Complete {str(__class__)}: {sys._getframe().f_code.co_name}')
 
+    # noinspection SpellCheckingInspection,PyUnresolvedReferences,PyProtectedMember
     def Calculate(self, model: Model):
         """
         The Calculate function calculates the values of all the parameters that are calculated by this object.
@@ -87,13 +111,31 @@ def Calculate(self, model: Model):
         # calculate non-dimensional temperature array
         Twnd = []
         try:
-            for t in range(1, len(model.reserv.timevector.value)):
-                Twnd = Twnd + [float(invertlaplace(fp, td[t], method='stehfest'))]
-        except:
-            msg = ('Error: GEOPHIRES could not execute numerical inverse laplace calculation for reservoir model 1. '
-                   'Simulation will abort.')
+            # Determine the required precision for this calculation.
+            precision = (
+                self.gringarten_stehfest_precision.value
+                if self.gringarten_stehfest_precision.Provided
+                else mp.dps
+            )
+
+            # Use the workdps context manager for a robust, thread-safe solution.
+            # It temporarily sets mp.dps for this block and restores it automatically.
+            with workdps(precision):
+                twnd_list = [
+                    float(invertlaplace(fp, td[t], method='stehfest'))
+                    for t in range(1, len(model.reserv.timevector.value))
+                ]
+
+                Twnd = np.asarray(twnd_list)
+
+        except Exception as e_:
+            msg = (
+                f'Error: GEOPHIRES could not execute numerical inverse laplace calculation for reservoir model 1 '
+                f'({self.gringarten_stehfest_precision.Name} = {precision}). '
+                'Simulation will abort.'
+            )
             print(msg)
-            raise RuntimeError(msg)
+            raise RuntimeError(msg) from e_
 
         Twnd = np.asarray(Twnd)
 
@@ -102,3 +144,4 @@ def Calculate(self, model: Model):
         model.reserv.Tresoutput.value = np.append([model.reserv.Trock.value], model.reserv.Tresoutput.value)
 
         model.logger.info(f'Complete {str(__class__)}: {sys._getframe().f_code.co_name}')
+

src/geophires_x/__init__.py

@@ -1 +1 @@
-__version__ = '3.9.50'
+__version__ = '3.9.52'

src/geophires_x_client/geophires_x_result.py

@@ -346,7 +346,10 @@ class GeophiresXResult:
                 'Average Annual Total Heating Production',
                 'Average Annual Electricity Use for Pumping',
             ],
-            'Simulation Metadata': [_StringValueField('GEOPHIRES Version')],
+            'Simulation Metadata': [
+                _StringValueField('GEOPHIRES Version'),
+                'Calculation Time',
+            ],
         }
     )
 

src/geophires_x_schema_generator/geophires-request.json

@@ -412,6 +412,15 @@
       "minimum": 0,
       "maximum": 0.2
     },
+    "Gringarten-Stehfest Precision": {
+      "description": "Sets the numerical precision (decimal places) for the inverse Laplace transform (Stehfest algorithm) used in the Gringarten calculation for the Multiple Parallel Fractures Reservoir Model. The default value provides maximum result stability; lower values calculate faster but may reduce consistency.",
+      "type": "integer",
+      "units": null,
+      "category": "Reservoir",
+      "default": 15,
+      "minimum": 8,
+      "maximum": 15
+    },
     "Cylindrical Reservoir Input Depth": {
       "description": "Depth of the inflow end of a cylindrical reservoir",
       "type": "number",

src/geophires_x_schema_generator/geophires-result.json

@@ -555,7 +555,8 @@
     "Simulation Metadata": {
       "type": "object",
       "properties": {
-        "GEOPHIRES Version": {}
+        "GEOPHIRES Version": {},
+        "Calculation Time": {}
       }
     }
   }

tests/example1_addons.csv

@@ -98,7 +98,8 @@ SURFACE EQUIPMENT SIMULATION RESULTS,Average Annual Net Electricity Generation,,
 SURFACE EQUIPMENT SIMULATION RESULTS,Average Pumping Power,,0.2,MW
 SURFACE EQUIPMENT SIMULATION RESULTS,Initial pumping power/net installed power,,3.82,%
 SURFACE EQUIPMENT SIMULATION RESULTS,Heat to Power Conversion Efficiency,,10.07,%
-Simulation Metadata,GEOPHIRES Version,,3.9.44,
+Simulation Metadata,GEOPHIRES Version,,3.9.49,
+Simulation Metadata,Calculation Time,,0.852,sec
 POWER GENERATION PROFILE,THERMAL DRAWDOWN,1,1.0,
 POWER GENERATION PROFILE,THERMAL DRAWDOWN,2,1.0056,
 POWER GENERATION PROFILE,THERMAL DRAWDOWN,3,1.0073,

tests/examples/example1_addons.out

@@ -4,10 +4,10 @@
 
 Simulation Metadata
 ----------------------
- GEOPHIRES Version: 3.9.44
- Simulation Date: 2025-07-29
- Simulation Time:  06:46
- Calculation Time:      0.898 sec
+ GEOPHIRES Version: 3.9.49
+ Simulation Date: 2025-08-21
+ Simulation Time:  09:19
+ Calculation Time:      0.852 sec
 
                            ***SUMMARY OF RESULTS***
 

tests/geophires-result_example-3.csv

@@ -82,6 +82,7 @@ SURFACE EQUIPMENT SIMULATION RESULTS,Average Annual Heat Production,,87.2,GWh
 SURFACE EQUIPMENT SIMULATION RESULTS,Average Pumping Power,,0.18,MW
 SURFACE EQUIPMENT SIMULATION RESULTS,Initial pumping power/net installed power,,100.0,%
 Simulation Metadata,GEOPHIRES Version,,3.0,
+Simulation Metadata,Calculation Time,,0.057,sec
 POWER GENERATION PROFILE,THERMAL DRAWDOWN,0,1.0,
 POWER GENERATION PROFILE,THERMAL DRAWDOWN,1,1.0081,
 POWER GENERATION PROFILE,THERMAL DRAWDOWN,2,1.0097,

tests/geophires_x_tests/test_reservoir.py

@@ -1,8 +1,10 @@
 from __future__ import annotations
 
+import copy
 import os
 import sys
 from pathlib import Path
+from typing import Any
 
 from pint.facets.plain import PlainQuantity
 
@@ -12,10 +14,12 @@
 from geophires_x_client import GeophiresInputParameters
 from geophires_x_client import GeophiresXClient
 from geophires_x_client import GeophiresXResult
+from geophires_x_client import ImmutableGeophiresInputParameters
 from tests.base_test_case import BaseTestCase
 
 
 class ReservoirTestCase(BaseTestCase):
+
     def test_lithostatic_pressure(self):
         p = static_pressure_MPa(2700, 3000)
         self.assertEqual(79.433865, p)
@@ -32,6 +36,72 @@ def test_reservoir_lithostatic_pressure(self):
         self.assertAlmostEqual(79.433865, p.magnitude, places=3)
         self.assertEqual('megapascal', p.units)
 
+    def test_gringarten_stehfest_precision(self):
+        def _log(msg) -> None:
+            print(f'[DEBUG][test_gringarten_stehfest_precision] {msg}')
+
+        param_name = 'Gringarten-Stehfest Precision'
+
+        def _get_result(gringarten_stehfest_precision: int) -> GeophiresXResult:
+            return GeophiresXClient(enable_caching=False).get_geophires_result(
+                ImmutableGeophiresInputParameters(
+                    from_file_path=self._get_test_file_path('generic-egs-case.txt'),
+                    params={param_name: gringarten_stehfest_precision},
+                )
+            )
+
+        _ = _get_result(15)  # warm up any caching
+        result_15 = _get_result(15)
+        result_8 = _get_result(8)
+
+        def calc_time(r: GeophiresXResult) -> float:
+            return r.result['Simulation Metadata']['Calculation Time']['value']
+
+        calc_time_15_sec = calc_time(result_15)
+        calc_time_8_sec = calc_time(result_8)
+
+        _log(f'calc_time_15_sec={calc_time_15_sec}, calc_time_8_sec={calc_time_8_sec}')
+
+        self.assertLess(calc_time_8_sec, calc_time_15_sec)
+
+        max_calc_time_8_sec = 1.0
+        try:
+            self.assertLessEqual(calc_time_8_sec, 1.0)
+        except AssertionError:
+            _log(
+                f'[WARNING] Calculation time for {param_name}=8 was greater than the expected maximum of '
+                f'{max_calc_time_8_sec} seconds. This may indicate a performance regression, '
+                f'depending on the available compute resources.'
+            )
+
+        speedup_pct = ((calc_time_15_sec - calc_time_8_sec) / calc_time_15_sec) * 100
+        _log(f'Speedup: {speedup_pct:.2f}%')
+
+        min_expected_speedup_pct = 25.0
+        try:
+            self.assertGreaterEqual(min_expected_speedup_pct, min_expected_speedup_pct)
+        except AssertionError:
+            _log(
+                f'[WARNING] Speedup for {param_name}=8 was less than the expected minimum of '
+                f'{min_expected_speedup_pct}%. This may indicate a performance regression.'
+            )
+
+        def no_metadata(r: GeophiresXResult) -> dict[str, Any]:
+            ret = copy.deepcopy(r.result)
+            del ret['Simulation Metadata']
+            del ret['metadata']
+            return ret
+
+        result_12_nm = no_metadata(_get_result(12))
+        result_15_nm = no_metadata(result_15)
+        try:
+            self.assertDictAlmostEqual(result_15_nm, result_12_nm, percent=1)
+        except AssertionError as ae:
+            try:
+                self.assertDictEqual(result_15_nm, result_12_nm)
+            except AssertionError as ae_with_dict_diff:
+                raise ae from ae_with_dict_diff
+
     # noinspection PyMethodMayBeStatic
     def _new_model(self, input_file=None) -> Model:
         stash_cwd = Path.cwd()

tests/test_geophires_x.py

@@ -67,6 +67,9 @@ def test_geophires_x_end_use_direct_use_heat(self):
         if 'metadata' in result_same_input.result:
             del result_same_input.result['metadata']
 
+        del result.result['Simulation Metadata']['Calculation Time']
+        del result_same_input.result['Simulation Metadata']['Calculation Time']
+
         self.assertDictEqual(result.result, result_same_input.result)
 
     # noinspection PyMethodMayBeStatic
@@ -362,6 +365,10 @@ def test_RTES_name(self):
         self.assertEqual(PlantType.RTES.value, 'Reservoir Thermal Energy Storage')
 
     def test_input_unit_conversion(self):
+        def delete_metadata(r: GeophiresXResult) -> None:
+            del r.result['metadata']
+            del r.result['Simulation Metadata']['Calculation Time']
+
         client = GeophiresXClient()
 
         result_meters_input = client.get_geophires_result(
@@ -371,7 +378,7 @@ def test_input_unit_conversion(self):
                 )
             )
         )
-        del result_meters_input.result['metadata']
+        delete_metadata(result_meters_input)
 
         result_kilometers_input = client.get_geophires_result(
             GeophiresInputParameters(
@@ -380,7 +387,7 @@ def test_input_unit_conversion(self):
                 )
             )
         )
-        del result_kilometers_input.result['metadata']
+        delete_metadata(result_kilometers_input)
 
         self.assertDictEqual(result_kilometers_input.result, result_meters_input.result)