diff --git a/.gitignore b/.gitignore index 37905002381e4d47e31fcc450d81bf593ee71a2d..29cf53ef64723c8b6822c87dc308132d623475bf 100644 --- a/.gitignore +++ b/.gitignore @@ -10,3 +10,4 @@ dist tests/resources/test_output +pydoubletcalc_install \ No newline at end of file diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 89e92217fd584f5b4897e288639f6b52b0cce93a..e61c1f5a94a752043177e063d0c84559f862de0f 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -17,6 +17,7 @@ Run tests: - export PATH=$JAVA_HOME/bin:$PATH - export THERMOGIS_JAR=$CI_PROJECT_DIR/resources/thermogis_jar/thermogis-1.7.0-shaded.jar - pixi install + - pixi run pip install pydoubletcalc --index-url https://gitlab-ci-token:$PYDOUBLETCALC_READ@ci.tno.nl/gitlab/api/v4/projects/19719/packages/pypi/simple - pixi run pytest -rx tests/ - rm -rf /var/lib/apt/lists/* only: diff --git a/README.md b/README.md index 9153aea91330a90c9bb7ca77092d0a1a2903874e..32e887b2335775a2e616fd331f718cff1af76da5 100644 --- a/README.md +++ b/README.md @@ -18,6 +18,27 @@ PyThermoGIS has been designed to be used as a python package you import into you It works by creating a python API access to the ThermoGIS techno-economic application, which is written in Java. Because of this dependency you need to Install a Java 17 VM and store the ThermoGIS Jar (a Java executable file) on your computer: +## Installing pydoubletcalc + +[pydoubletcalc](https://ci.tno.nl/gitlab/AGS/Geothermal/thermogis/pydoubletcalc) is the python backend we are writing to replace the java backend for pythermogis. + +It is not yet public, until then you can install the latest version of the code using the following command: + +```bash +pixi run pip install pydoubletcalc --index-url https://token_name:token_password@ci.tno.nl/gitlab/api/v4/projects/19719/packages/pypi/simple +``` +replace token_name and token_password with the credentials from a Personal Access Token, with read access. +(If you don't have a PAT then: Login to gitlab > preferences > Access Tokens) + +To update pydoubletcalc you will first need to uninstall it and then reinstall using the above method. + + +## Installing pydoubletcalc locally + +```bash + pixi add --editable pydoubletcalc@file:C:path\to\pydoubletcalc --pypi +``` + ### 1. Install Java 17 and Download the ThermoGIS JAR This package requires a Java 17 VM (we recommend using [Amazon Corretto 17](https://docs.aws.amazon.com/corretto/latest/corretto-17-ug/downloads-list.html)) and a [ThermoGIS Jar file](https://ci.tno.nl/gitlab/ags_public/pythermogis/-/blob/main/resources/thermogis_jar/thermogis-1.7.0-shaded.jar?ref_type=heads) (Version >=1.7). diff --git a/pixi.lock b/pixi.lock index 0f624af11f635a8bc949d06af2a35a9f94b3aa63..e0f720c2709a21cbfcbcf64c7cbe2b2cc8018d4d 100644 --- a/pixi.lock +++ b/pixi.lock @@ -140,6 +140,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/linux-64/libxcb-1.17.0-h8a09558_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libxml2-2.13.8-h4bc477f_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libzlib-1.3.1-hb9d3cd8_2.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/llvmlite-0.45.1-py313hdd307be_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/locket-1.0.0-pyhd8ed1ab_0.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/linux-64/lz4-4.4.4-py313h8756d67_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/lz4-c-1.10.0-h5888daf_1.conda @@ -162,6 +163,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/linux-64/nh3-0.2.21-py39h77e2912_1.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/nlohmann_json-3.12.0-h3f2d84a_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/nodeenv-1.9.1-pyhd8ed1ab_1.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/numba-0.62.1-py313hd8e3f9f_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/numpy-2.2.5-py313h17eae1a_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/openjpeg-2.5.3-h5fbd93e_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/openpyxl-3.1.5-py313h9c9eb82_1.conda @@ -204,6 +206,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/requests-toolbelt-1.0.0-pyhd8ed1ab_1.conda - conda: https://conda.anaconda.org/conda-forge/noarch/rfc3986-2.0.0-pyhd8ed1ab_1.conda - conda: https://conda.anaconda.org/conda-forge/noarch/rich-14.0.0-pyh29332c3_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/ruff-0.14.8-h813ae00_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/s2n-1.5.21-h7ab7c64_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/secretstorage-3.3.3-py313h78bf25f_3.conda - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-75.8.2-pyhff2d567_0.conda @@ -250,6 +253,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/8f/e9/6a7d025d8da8c4931522922cd706105aa32b3291d1add8c5427cdcd66e63/kiwisolver-1.4.8-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl - pypi: https://files.pythonhosted.org/packages/f5/64/41c4367bcaecbc03ef0d2a3ecee58a7065d0a36ae1aa817fe573a2da66d4/matplotlib-3.10.3-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl - pypi: https://files.pythonhosted.org/packages/d1/80/b9c19f1bb4ac6c5fa6f94a4f278bc68a778473d1814a86a375d7cffa193a/netCDF4-1.7.2-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl + - pypi: https://files.pythonhosted.org/packages/78/f9/690a8600b93c332de3ab4a344a4ac34f00c8f104917061f779db6a918ed6/pathlib-1.0.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/02/65/ad2bc85f7377f5cfba5d4466d5474423a3fb7f6a97fd807c06f92dd3e721/plotly-6.0.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/c9/58/c3bc54c0fad9a82899e9a2703e04ee6e8eaa76caa90c0689fd1b468a4427/pygridsio-1.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/ea/00/d815833441d8c52bf4a6930952e77d3de77d0bf67b3202ccc12dabdae279/pykrige-1.7.2.tar.gz @@ -259,9 +263,11 @@ environments: - pypi: https://files.pythonhosted.org/packages/bd/24/12818598c362d7f300f18e74db45963dbcb85150324092410c8b49405e42/pyproject_hooks-1.2.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/05/19/94d6c66184c7d0f9374330c714f62c147dbb53eda9efdcc8fc6e2ac454c5/rasterio-1.4.3-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl - pypi: https://files.pythonhosted.org/packages/2a/2f/63d2cacc0e525f8e3398bcf32bd3620385f22cd1600834ec49d7f3597a7b/rioxarray-0.19.0-py3-none-any.whl - - pypi: https://files.pythonhosted.org/packages/b5/09/c5b6734a50ad4882432b6bb7c02baf757f5b2f256041da5df242e2d7e6b6/scipy-1.15.3-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl + - pypi: https://files.pythonhosted.org/packages/21/f6/4bfb5695d8941e5c570a04d9fcd0d36bce7511b7d78e6e75c8f9791f82d0/scipy-1.16.3-cp313-cp313-manylinux2014_x86_64.manylinux_2_17_x86_64.whl - pypi: https://files.pythonhosted.org/packages/d1/a7/5c9cb413e4e2ce52c16be717e94abd40ce91b1f8974624d5d56154c5d40b/shapely-2.1.0-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl + - pypi: https://files.pythonhosted.org/packages/d0/30/dc54f88dd4a2b5dc8a0279bdd7270e735851848b762aeb1c1184ed1f6b14/tqdm-4.67.1-py3-none-any.whl - pypi: ./ + - pypi: C:/Users/knappersfy/work/thermogis/pydoubletcalc win-64: - conda: https://conda.anaconda.org/conda-forge/win-64/_openmp_mutex-4.5-2_gnu.conda - conda: https://conda.anaconda.org/conda-forge/noarch/_python_abi3_support-1.0-hd8ed1ab_2.conda @@ -371,6 +377,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/win-64/libxcb-1.17.0-h0e4246c_0.conda - conda: https://conda.anaconda.org/conda-forge/win-64/libxml2-2.13.8-h442d1da_0.conda - conda: https://conda.anaconda.org/conda-forge/win-64/libzlib-1.3.1-h2466b09_2.conda + - conda: https://conda.anaconda.org/conda-forge/win-64/llvmlite-0.45.1-py313h5c49287_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/locket-1.0.0-pyhd8ed1ab_0.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/win-64/lz4-4.4.4-py313h05901a4_0.conda - conda: https://conda.anaconda.org/conda-forge/win-64/lz4-c-1.10.0-h2466b09_1.conda @@ -392,6 +399,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-1.43.1-pyhe01879c_0.conda - conda: https://conda.anaconda.org/conda-forge/win-64/nh3-0.2.21-py39he870945_1.conda - conda: https://conda.anaconda.org/conda-forge/noarch/nodeenv-1.9.1-pyhd8ed1ab_1.conda + - conda: https://conda.anaconda.org/conda-forge/win-64/numba-0.62.1-py313h924e429_0.conda - conda: https://conda.anaconda.org/conda-forge/win-64/numpy-2.2.5-py313hefb8edb_0.conda - conda: https://conda.anaconda.org/conda-forge/win-64/openjpeg-2.5.3-h4d64b90_0.conda - conda: https://conda.anaconda.org/conda-forge/win-64/openpyxl-3.1.5-py313he57e174_1.conda @@ -432,6 +440,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/requests-toolbelt-1.0.0-pyhd8ed1ab_1.conda - conda: https://conda.anaconda.org/conda-forge/noarch/rfc3986-2.0.0-pyhd8ed1ab_1.conda - conda: https://conda.anaconda.org/conda-forge/noarch/rich-14.0.0-pyh29332c3_0.conda + - conda: https://conda.anaconda.org/conda-forge/win-64/ruff-0.14.8-h15e3a1f_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-75.8.2-pyhff2d567_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/shellingham-1.5.4-pyhd8ed1ab_1.conda - conda: https://conda.anaconda.org/conda-forge/noarch/six-1.17.0-pyhd8ed1ab_0.conda @@ -482,6 +491,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/d0/dc/c1abe38c37c071d0fc71c9a474fd0b9ede05d42f5a458d584619cfd2371a/kiwisolver-1.4.8-cp313-cp313-win_amd64.whl - pypi: https://files.pythonhosted.org/packages/b1/0f/eed564407bd4d935ffabf561ed31099ed609e19287409a27b6d336848653/matplotlib-3.10.3-cp313-cp313-win_amd64.whl - pypi: https://files.pythonhosted.org/packages/66/b5/e04550fd53de57001dbd5a87242da7ff784c80790adc48897977b6ccf891/netCDF4-1.7.2-cp313-cp313-win_amd64.whl + - pypi: https://files.pythonhosted.org/packages/78/f9/690a8600b93c332de3ab4a344a4ac34f00c8f104917061f779db6a918ed6/pathlib-1.0.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/02/65/ad2bc85f7377f5cfba5d4466d5474423a3fb7f6a97fd807c06f92dd3e721/plotly-6.0.1-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/c9/58/c3bc54c0fad9a82899e9a2703e04ee6e8eaa76caa90c0689fd1b468a4427/pygridsio-1.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/ea/00/d815833441d8c52bf4a6930952e77d3de77d0bf67b3202ccc12dabdae279/pykrige-1.7.2.tar.gz @@ -491,9 +501,11 @@ environments: - pypi: https://files.pythonhosted.org/packages/bd/24/12818598c362d7f300f18e74db45963dbcb85150324092410c8b49405e42/pyproject_hooks-1.2.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/df/88/9db5f49ebfdd9c12365e4cac76c34ccb1a642b1c8cbab4124b3c681495de/rasterio-1.4.3-cp313-cp313-win_amd64.whl - pypi: https://files.pythonhosted.org/packages/2a/2f/63d2cacc0e525f8e3398bcf32bd3620385f22cd1600834ec49d7f3597a7b/rioxarray-0.19.0-py3-none-any.whl - - pypi: https://files.pythonhosted.org/packages/87/2e/892ad2862ba54f084ffe8cc4a22667eaf9c2bcec6d2bff1d15713c6c0703/scipy-1.15.3-cp313-cp313-win_amd64.whl + - pypi: https://files.pythonhosted.org/packages/cd/01/1204382461fcbfeb05b6161b594f4007e78b6eba9b375382f79153172b4d/scipy-1.16.3-cp313-cp313-win_amd64.whl - pypi: https://files.pythonhosted.org/packages/e3/f0/9f8cdf2258d7aed742459cea51c70d184de92f5d2d6f5f7f1ded90a18c31/shapely-2.1.0-cp313-cp313-win_amd64.whl + - pypi: https://files.pythonhosted.org/packages/d0/30/dc54f88dd4a2b5dc8a0279bdd7270e735851848b762aeb1c1184ed1f6b14/tqdm-4.67.1-py3-none-any.whl - pypi: ./ + - pypi: C:/Users/knappersfy/work/thermogis/pydoubletcalc packages: - conda: https://conda.anaconda.org/conda-forge/linux-64/_libgcc_mutex-0.1-conda_forge.tar.bz2 sha256: fe51de6107f9edc7aa4f786a70f4a883943bc9d39b3bb7307c04c41410990726 @@ -3622,6 +3634,40 @@ packages: purls: [] size: 55476 timestamp: 1727963768015 +- conda: https://conda.anaconda.org/conda-forge/linux-64/llvmlite-0.45.1-py313hdd307be_0.conda + sha256: 5b8e2063e93bea160fd50274d05ce4436f01c383a392f293b769dfb973c4df21 + md5: 5afb15643ef1fcea20798bb6086bb3f9 + depends: + - __glibc >=2.17,<3.0.a0 + - libgcc >=14 + - libstdcxx >=14 + - libzlib >=1.3.1,<2.0a0 + - python >=3.13,<3.14.0a0 + - python_abi 3.13.* *_cp313 + - zstd >=1.5.7,<1.6.0a0 + license: BSD-2-Clause + license_family: BSD + purls: + - pkg:pypi/llvmlite?source=hash-mapping + size: 34153671 + timestamp: 1759394632193 +- conda: https://conda.anaconda.org/conda-forge/win-64/llvmlite-0.45.1-py313h5c49287_0.conda + sha256: 0b63923082e724b2c2939621aef77d9ec65aa468a7b29917a850e47e2083adda + md5: d946ee3e7228e48270589791871a891e + depends: + - libzlib >=1.3.1,<2.0a0 + - python >=3.13,<3.14.0a0 + - python_abi 3.13.* *_cp313 + - ucrt >=10.0.20348.0 + - vc >=14.3,<15 + - vc14_runtime >=14.44.35208 + - zstd >=1.5.7,<1.6.0a0 + license: BSD-2-Clause + license_family: BSD + purls: + - pkg:pypi/llvmlite?source=hash-mapping + size: 22905696 + timestamp: 1759394712687 - conda: https://conda.anaconda.org/conda-forge/noarch/locket-1.0.0-pyhd8ed1ab_0.tar.bz2 sha256: 9afe0b5cfa418e8bdb30d8917c5a6cec10372b037924916f1f85b9f4899a67a6 md5: 91e27ef3d05cc772ce627e51cff111c4 @@ -4081,6 +4127,57 @@ packages: - pkg:pypi/nodeenv?source=hash-mapping size: 34574 timestamp: 1734112236147 +- conda: https://conda.anaconda.org/conda-forge/linux-64/numba-0.62.1-py313hd8e3f9f_0.conda + sha256: 0b0ccf81ecd0cfb934818c3fb88404821904fe84e67815ab9958d37b0dca61e4 + md5: 82ffdc573a667626351f4110605da846 + depends: + - __glibc >=2.17,<3.0.a0 + - _openmp_mutex >=4.5 + - libgcc >=14 + - libstdcxx >=14 + - llvmlite >=0.45.0,<0.46.0a0 + - numpy >=1.23,<3 + - numpy >=1.24,<2.4 + - python >=3.13,<3.14.0a0 + - python_abi 3.13.* *_cp313 + constrains: + - cudatoolkit >=11.2 + - cuda-version >=11.2 + - tbb >=2021.6.0 + - cuda-python >=11.6 + - scipy >=1.0 + - libopenblas !=0.3.6 + license: BSD-2-Clause + license_family: BSD + purls: + - pkg:pypi/numba?source=hash-mapping + size: 5743830 + timestamp: 1759165232580 +- conda: https://conda.anaconda.org/conda-forge/win-64/numba-0.62.1-py313h924e429_0.conda + sha256: 79835953985d64565f76f912517ab5700148e86659b8e79ecd2d0e6d7377ac46 + md5: ae201f33cbcbb2aba93daf4b3263b4a5 + depends: + - llvmlite >=0.45.0,<0.46.0a0 + - numpy >=1.23,<3 + - numpy >=1.24,<2.4 + - python >=3.13,<3.14.0a0 + - python_abi 3.13.* *_cp313 + - ucrt >=10.0.20348.0 + - vc >=14.3,<15 + - vc14_runtime >=14.44.35208 + constrains: + - libopenblas !=0.3.6 + - cudatoolkit >=11.2 + - tbb >=2021.6.0 + - scipy >=1.0 + - cuda-python >=11.6 + - cuda-version >=11.2 + license: BSD-2-Clause + license_family: BSD + purls: + - pkg:pypi/numba?source=hash-mapping + size: 5706597 + timestamp: 1759165298367 - conda: https://conda.anaconda.org/conda-forge/linux-64/numpy-2.2.5-py313h17eae1a_0.conda sha256: c0a200d0e53a1acbfa1d1e2277e3337ea2aa8cb584448790317a98c62dcaebce md5: 6ceeff9ed72e54e4a2f9a1c88f47bdde @@ -4382,6 +4479,10 @@ packages: - pkg:pypi/partd?source=hash-mapping size: 20884 timestamp: 1715026639309 +- pypi: https://files.pythonhosted.org/packages/78/f9/690a8600b93c332de3ab4a344a4ac34f00c8f104917061f779db6a918ed6/pathlib-1.0.1-py3-none-any.whl + name: pathlib + version: 1.0.1 + sha256: f35f95ab8b0f59e6d354090350b44a80a80635d22efdedfa84c7ad1cf0a74147 - conda: https://conda.anaconda.org/conda-forge/noarch/pathspec-0.12.1-pyhd8ed1ab_1.conda sha256: 9f64009cdf5b8e529995f18e03665b03f5d07c0b17445b8badef45bde76249ee md5: 617f15191456cc6a13db418a275435e5 @@ -4663,6 +4764,20 @@ packages: - pkg:pypi/pycparser?source=hash-mapping size: 110100 timestamp: 1733195786147 +- pypi: C:/Users/knappersfy/work/thermogis/pydoubletcalc + name: pydoubletcalc + version: 0.0.3 + sha256: a391e31a4d31cad94ef6af7520fb9418f2c7f1741fa28384734097c9fda1b1fd + requires_dist: + - pandas + - matplotlib + - pathlib + - xarray + - tqdm>=4.67.1 + - scipy>=1.16.2 + - numba>=0.62.1 + requires_python: '>=3.13,<3.14' + editable: true - conda: https://conda.anaconda.org/conda-forge/noarch/pygments-2.19.1-pyhd8ed1ab_0.conda sha256: 28a3e3161390a9d23bc02b4419448f8d27679d9e2c250e29849e37749c8de86b md5: 232fb4577b6687b2d503ef8e254270c9 @@ -4830,7 +4945,7 @@ packages: - pypi: ./ name: pythermogis version: 1.2.0 - sha256: 50c9d857c9314b246b28d7e610647475414c718d839408818aa4f441813b51a7 + sha256: 35de67af8136fb6442b90f5a2956e6ac57f8a6f52f1fbbfad15768e83e08b31f requires_dist: - jpype1>=1.5.2,<2 - xarray==2024.9.0.* @@ -5216,6 +5331,37 @@ packages: - scipy ; extra == 'interp' - scipy ; extra == 'all' requires_python: '>=3.10' +- conda: https://conda.anaconda.org/conda-forge/linux-64/ruff-0.14.8-h813ae00_0.conda + noarch: python + sha256: 4adf379daccb73f03297a6966d1200f6ea65e6a1513d749e7f782e32267fe2bb + md5: 295ce05c06920527a581a5e148a4eec6 + depends: + - python + - __glibc >=2.17,<3.0.a0 + - libgcc >=14 + constrains: + - __glibc >=2.17 + license: MIT + license_family: MIT + purls: + - pkg:pypi/ruff?source=hash-mapping + size: 11340280 + timestamp: 1764866215629 +- conda: https://conda.anaconda.org/conda-forge/win-64/ruff-0.14.8-h15e3a1f_0.conda + noarch: python + sha256: fbcaafffd55c7022464219b95658d38980ee04bb001d35c3d97e2e933d7c6bf7 + md5: 35ec53f16d22dc8b17e17865a98c2120 + depends: + - python + - vc >=14.3,<15 + - vc14_runtime >=14.44.35208 + - ucrt >=10.0.20348.0 + license: MIT + license_family: MIT + purls: + - pkg:pypi/ruff?source=hash-mapping + size: 11874411 + timestamp: 1764866263950 - conda: https://conda.anaconda.org/conda-forge/linux-64/s2n-1.5.21-h7ab7c64_0.conda sha256: c8b252398b502a5cc6ea506fd2fafe7e102e7c9e2ef48b7813566e8a72ce2205 md5: 28b5a7895024a754249b2ad7de372faa @@ -5228,13 +5374,13 @@ packages: purls: [] size: 358164 timestamp: 1749095480268 -- pypi: https://files.pythonhosted.org/packages/87/2e/892ad2862ba54f084ffe8cc4a22667eaf9c2bcec6d2bff1d15713c6c0703/scipy-1.15.3-cp313-cp313-win_amd64.whl +- pypi: https://files.pythonhosted.org/packages/21/f6/4bfb5695d8941e5c570a04d9fcd0d36bce7511b7d78e6e75c8f9791f82d0/scipy-1.16.3-cp313-cp313-manylinux2014_x86_64.manylinux_2_17_x86_64.whl name: scipy - version: 1.15.3 - sha256: b90ab29d0c37ec9bf55424c064312930ca5f4bde15ee8619ee44e69319aab163 + version: 1.16.3 + sha256: 7dc1360c06535ea6116a2220f760ae572db9f661aba2d88074fe30ec2aa1ff88 requires_dist: - - numpy>=1.23.5,<2.5 - - pytest ; extra == 'test' + - numpy>=1.25.2,<2.6 + - pytest>=8.0.0 ; extra == 'test' - pytest-cov ; extra == 'test' - pytest-timeout ; extra == 'test' - pytest-xdist ; extra == 'test' @@ -5245,11 +5391,11 @@ packages: - scikit-umfpack ; extra == 'test' - pooch ; extra == 'test' - hypothesis>=6.30 ; extra == 'test' - - array-api-strict>=2.0,<2.1.1 ; extra == 'test' + - array-api-strict>=2.3.1 ; extra == 'test' - cython ; extra == 'test' - meson ; extra == 'test' - ninja ; sys_platform != 'emscripten' and extra == 'test' - - sphinx>=5.0.0,<8.0.0 ; extra == 'doc' + - sphinx>=5.0.0,<8.2.0 ; extra == 'doc' - intersphinx-registry ; extra == 'doc' - pydata-sphinx-theme>=0.15.2 ; extra == 'doc' - sphinx-copybutton ; extra == 'doc' @@ -5257,10 +5403,11 @@ packages: - matplotlib>=3.5 ; extra == 'doc' - numpydoc ; extra == 'doc' - jupytext ; extra == 'doc' - - myst-nb ; extra == 'doc' + - myst-nb>=1.2.0 ; extra == 'doc' - pooch ; extra == 'doc' - jupyterlite-sphinx>=0.19.1 ; extra == 'doc' - jupyterlite-pyodide-kernel ; extra == 'doc' + - linkify-it-py ; extra == 'doc' - mypy==1.10.0 ; extra == 'dev' - typing-extensions ; extra == 'dev' - types-psutil ; extra == 'dev' @@ -5270,14 +5417,14 @@ packages: - rich-click ; extra == 'dev' - doit>=0.36.0 ; extra == 'dev' - pydevtool ; extra == 'dev' - requires_python: '>=3.10' -- pypi: https://files.pythonhosted.org/packages/b5/09/c5b6734a50ad4882432b6bb7c02baf757f5b2f256041da5df242e2d7e6b6/scipy-1.15.3-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl + requires_python: '>=3.11' +- pypi: https://files.pythonhosted.org/packages/cd/01/1204382461fcbfeb05b6161b594f4007e78b6eba9b375382f79153172b4d/scipy-1.16.3-cp313-cp313-win_amd64.whl name: scipy - version: 1.15.3 - sha256: c9deabd6d547aee2c9a81dee6cc96c6d7e9a9b1953f74850c179f91fdc729cb7 + version: 1.16.3 + sha256: 062246acacbe9f8210de8e751b16fc37458213f124bef161a5a02c7a39284304 requires_dist: - - numpy>=1.23.5,<2.5 - - pytest ; extra == 'test' + - numpy>=1.25.2,<2.6 + - pytest>=8.0.0 ; extra == 'test' - pytest-cov ; extra == 'test' - pytest-timeout ; extra == 'test' - pytest-xdist ; extra == 'test' @@ -5288,11 +5435,11 @@ packages: - scikit-umfpack ; extra == 'test' - pooch ; extra == 'test' - hypothesis>=6.30 ; extra == 'test' - - array-api-strict>=2.0,<2.1.1 ; extra == 'test' + - array-api-strict>=2.3.1 ; extra == 'test' - cython ; extra == 'test' - meson ; extra == 'test' - ninja ; sys_platform != 'emscripten' and extra == 'test' - - sphinx>=5.0.0,<8.0.0 ; extra == 'doc' + - sphinx>=5.0.0,<8.2.0 ; extra == 'doc' - intersphinx-registry ; extra == 'doc' - pydata-sphinx-theme>=0.15.2 ; extra == 'doc' - sphinx-copybutton ; extra == 'doc' @@ -5300,10 +5447,11 @@ packages: - matplotlib>=3.5 ; extra == 'doc' - numpydoc ; extra == 'doc' - jupytext ; extra == 'doc' - - myst-nb ; extra == 'doc' + - myst-nb>=1.2.0 ; extra == 'doc' - pooch ; extra == 'doc' - jupyterlite-sphinx>=0.19.1 ; extra == 'doc' - jupyterlite-pyodide-kernel ; extra == 'doc' + - linkify-it-py ; extra == 'doc' - mypy==1.10.0 ; extra == 'dev' - typing-extensions ; extra == 'dev' - types-psutil ; extra == 'dev' @@ -5313,7 +5461,7 @@ packages: - rich-click ; extra == 'dev' - doit>=0.36.0 ; extra == 'dev' - pydevtool ; extra == 'dev' - requires_python: '>=3.10' + requires_python: '>=3.11' - conda: https://conda.anaconda.org/conda-forge/linux-64/secretstorage-3.3.3-py313h78bf25f_3.conda sha256: 7f548e147e14ce743a796aa5f26aba11f82c14ab53ab25b48f35974ca48f6ac7 md5: 813e01c086f6c4e134e13ef25b02df8c @@ -5525,6 +5673,22 @@ packages: - pkg:pypi/tornado?source=hash-mapping size: 878044 timestamp: 1748003914685 +- pypi: https://files.pythonhosted.org/packages/d0/30/dc54f88dd4a2b5dc8a0279bdd7270e735851848b762aeb1c1184ed1f6b14/tqdm-4.67.1-py3-none-any.whl + name: tqdm + version: 4.67.1 + sha256: 26445eca388f82e72884e0d580d5464cd801a3ea01e63e5601bdff9ba6a48de2 + requires_dist: + - colorama ; sys_platform == 'win32' + - pytest>=6 ; extra == 'dev' + - pytest-cov ; extra == 'dev' + - pytest-timeout ; extra == 'dev' + - pytest-asyncio>=0.24 ; extra == 'dev' + - nbval ; extra == 'dev' + - requests ; extra == 'discord' + - slack-sdk ; extra == 'slack' + - requests ; extra == 'telegram' + - ipywidgets>=6 ; extra == 'notebook' + requires_python: '>=3.7' - conda: https://conda.anaconda.org/conda-forge/noarch/twine-6.1.0-pyh29332c3_0.conda sha256: c5b373f6512b96324c9607d7d91a76bb53c1056cb1012b4f9c86900c6b7f8898 md5: d319066fad04e07a0223bf9936090161 diff --git a/pyproject.toml b/pyproject.toml index 73e4690a10237544fccf468b68b6b676cf48ae6f..b0c22a3131b0acb86fdaf281056289941822f0c8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -36,7 +36,7 @@ markers = [ [tool.setuptools] package-dir = {"" = "src"} -[tool.pixi.project] +[tool.pixi.workspace] channels = ["conda-forge"] platforms = ["win-64", "linux-64"] @@ -69,6 +69,8 @@ dask = ">=2025.5.1,<2026" narwhals = ">=1.43.1,<2" pre-commit = ">=4.3.0,<5" python-dotenv = ">=1.2.1,<2" +numba = ">=0.62.1,<0.63" +ruff = ">=0.14.8,<0.15" [tool.ruff] line-length = 88 diff --git a/src/pythermogis/thermogis_classes/doublet.py b/src/pythermogis/thermogis_classes/doublet.py index 3360a824e7ba0769e48cfff1785798bdf95501d1..50d31633c4abef871eecaac9fed89b9e0e2ed200 100644 --- a/src/pythermogis/thermogis_classes/doublet.py +++ b/src/pythermogis/thermogis_classes/doublet.py @@ -1,6 +1,9 @@ import numpy as np import xarray as xr from jpype import JClass +from pythermogis.workflow.utc.doublet import calculate_doublet_performance, DoubletInput, DoubletOutput +from pythermogis.workflow.utc.utc_properties import UTCConfiguration + def simulate_doublet( output_data: xr.Dataset, @@ -102,27 +105,64 @@ def calculate_performance_of_single_location(mask: float, depth: float, thicknes if not np.isnan(mask) or temperature < utc_properties.minProdTemp(): return (mask_value,) * 14 - # Instantiate ThermoGIS doublet - doublet = instantiate_thermogis_doublet(utc_properties, rng_seed) - - DoubletInput = JClass("thermogis.calc.utc.doublet.records.DoubletInput") - input = DoubletInput( - -9999.0, # unknowninput - thickness, - transmissivity, - transmissivity_with_ntg, - ntg, - depth, - porosity, - temperature, - None, # ates input + use_java_backend = True + if use_java_backend: + doublet = instantiate_thermogis_doublet(utc_properties, rng_seed) + JavaDoubletInput = JClass("thermogis.calc.utc.doublet.records.DoubletInput") + input = JavaDoubletInput( + -9999.0, # unknowninput + thickness, + transmissivity, + transmissivity_with_ntg, + ntg, + depth, + porosity, + temperature, + None, # ates input + ) + # The Java routine which calculates DoubletPerformance, for more detail on the + # simulation inspect the Java source code + results_java = doublet.calculateDoubletPerformance(input) + results = DoubletOutput( + results_java.power(), + results_java.hppower(), + results_java.capex(), + 0.0, + 0.0, + results_java.opex(), + results_java.utc(), + results_java.npv(), + results_java.hprod(), + results_java.cop(), + results_java.cophp(), + results_java.pres(), + results_java.flow(), + results_java.welld(), + results_java.breakthrough(), + results_java.cooling(), + results_java.productionTemp(), + results_java.injectionTemp(), ) - - # The Java routine which calculates DoubletPerformance, for more detail on the simulation inspect the Java source code - results = doublet.calculateDoubletPerformance(input) + else: + # Arrange: instantiate default UTCConfiguration + props = UTCConfiguration(segment_length=1) + # Create a minimal valid DoubletInput + input_data = DoubletInput( + unknown_input_value=-9999.0, + thickness=thickness, + transmissivity=transmissivity, + transmissivity_with_ntg=transmissivity_with_ntg, + ntg=ntg, + depth=depth, + porosity=porosity, + temperature=temperature, + ) + # Act + #perform one simulation to compile all the numba functions before timing + results = calculate_doublet_performance(props, input_data) # If calculation was not successful, return mask value - if results.utc() == -9999.0: + if results.utc == -9999.0: return (mask_value,) * 14 # calculate net-present-value using the utc-cutoffs @@ -131,24 +171,23 @@ def calculate_performance_of_single_location(mask: float, depth: float, thicknes else: utc_cut = utc_properties.utcCutoff() - hprod = results.hprod() - npv = 1e-6 * (utc_cut - results.utc()) * 3.6 * hprod * (1 - utc_properties.taxRate()) + npv = 1e-6 * (utc_cut - results.utc) * 3.6 * results.hprod * (1 - utc_properties.taxRate()) # get values from doublet - output_values = {"power": results.power(), - "heat_pump_power": results.hppower(), - "capex": results.capex(), - "opex": results.opex(), - "utc": results.utc(), + output_values = {"power": results.power, + "heat_pump_power": results.hppower, + "capex": results.capex, + "opex": results.opex, + "utc": results.utc, "npv": npv, - "hprod": hprod, - "cop": results.cop(), - "cophp": results.cophp(), - "pres": results.pres(), - "flow_rate": results.flow(), - "welld": results.welld(), - "inj_temp": results.injectionTemp(), - "prd_temp": results.productionTemp(), + "hprod": results.hprod, + "cop": results.cop, + "cophp": results.cophp, + "pres": results.pres, + "flow_rate": results.flow, + "welld": results.welld, + "inj_temp": results.injection_temp, + "prd_temp": results.production_temp, } # Reset doublet variables for next calculation diff --git a/src/pythermogis/utils/__init__.py b/src/pythermogis/utils/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/src/pythermogis/utils/timer.py b/src/pythermogis/utils/timer.py new file mode 100644 index 0000000000000000000000000000000000000000..818a092f468448ca36da6a0b37deb72d220ab74c --- /dev/null +++ b/src/pythermogis/utils/timer.py @@ -0,0 +1,28 @@ +import sys +import timeit + + +def print_time(start, message: str = "", verbose: bool = True): + if not verbose: + return + + time_taken_seconds = timeit.default_timer() - start + time_taken_message = format_time(time_taken_seconds) + print(f"{message}{time_taken_message}", flush=True, file=sys.stdout) + + return timeit.default_timer() + + +def format_time(seconds: float) -> str: + hours = seconds // 3600 + minutes = (seconds % 3600) // 60 + seconds = seconds % 60 + + if hours > 0: + return f"{hours:.0f}h {minutes:.0f}m {seconds:.0f}s" + elif minutes > 0: + return f"{minutes:.0f}m {seconds:.0f}s" + elif seconds > 10: + return f"{seconds:.1f}s" + else: + return f"{seconds:.3f}s" diff --git a/src/pythermogis/workflow/__init__.py b/src/pythermogis/workflow/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/src/pythermogis/workflow/overview/__init__.py b/src/pythermogis/workflow/overview/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/src/pythermogis/workflow/potential/__init__.py b/src/pythermogis/workflow/potential/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/src/pythermogis/workflow/preprocess/__init__.py b/src/pythermogis/workflow/preprocess/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/src/pythermogis/workflow/utc/__init__.py b/src/pythermogis/workflow/utc/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/src/pythermogis/workflow/utc/constants.py b/src/pythermogis/workflow/utc/constants.py new file mode 100644 index 0000000000000000000000000000000000000000..23f54df698133184f25b126d77885276b43a8c46 --- /dev/null +++ b/src/pythermogis/workflow/utc/constants.py @@ -0,0 +1 @@ +DARCY_SI = 1.0e-12 / 1.01325 diff --git a/src/pythermogis/workflow/utc/cooling_temp.py b/src/pythermogis/workflow/utc/cooling_temp.py new file mode 100644 index 0000000000000000000000000000000000000000..cca4291d4fb5b1c36fb054bbd555817e51725e8b --- /dev/null +++ b/src/pythermogis/workflow/utc/cooling_temp.py @@ -0,0 +1,65 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING + +from pythermogis.workflow.utc.doublet_utils import calc_lifetime +from pythermogis.workflow.utc.flow import calculate_volumetric_flow + +if TYPE_CHECKING: + from pythermogis.workflow.utc.doublet import DoubletInput + from pythermogis.workflow.utc.utc_properties import UTCConfiguration + + +def calculate_cooling_temperature( + props: UTCConfiguration, + input: DoubletInput, + drawdown_pressure: float, + well_distance: float, + injection_temp: float, +) -> float: + results = calculate_volumetric_flow( + props=props, + input_data=input, + original_pressure=drawdown_pressure, + well_distance=well_distance, + injection_temp=injection_temp, + ) + + flowrate = min(results.flowrate, props.max_flow) + + cooling_frac_min = 0.0 + cooling_frac_max = 0.5 + cooling_frac = 0.5 * (cooling_frac_min + cooling_frac_max) + + for _ in range(1000): + if abs(cooling_frac_max - cooling_frac_min) <= 0.001: + break + + cooling_frac = 0.5 * (cooling_frac_min + cooling_frac_max) + + lifetime = calc_lifetime( + well_distance=well_distance, + thickness=input.thickness * input.ntg, + delta_temp_fraction=cooling_frac, + porosity=input.porosity, + flowrate=flowrate, + depth=input.depth, + reservoir_temp=input.temperature, + salinity_surface=props.salinity_surface, + salinity_gradient=props.salinity_gradient, + cp_rock=props.optim_dist_cp_rock, + rho_rock=props.optim_dist_rho_rock, + ) + + if lifetime < props.optim_dist_lifetime: + cooling_frac_min = cooling_frac + else: + cooling_frac_max = cooling_frac + + if (not props.is_ates) and cooling_frac > 0.3: + raise RuntimeError( + f"Large cooling factor ({cooling_frac}) " + f"could result in less reliable calculation" + ) + + return cooling_frac * (input.temperature - injection_temp) diff --git a/src/pythermogis/workflow/utc/doublet.py b/src/pythermogis/workflow/utc/doublet.py new file mode 100644 index 0000000000000000000000000000000000000000..61ecb11951f3a2783aa41e23217ed735ed6f8343 --- /dev/null +++ b/src/pythermogis/workflow/utc/doublet.py @@ -0,0 +1,209 @@ +import timeit +from dataclasses import dataclass +from typing import NamedTuple + +from pythermogis.utils.timer import print_time +from pythermogis.workflow.utc.constants import DARCY_SI +from pythermogis.workflow.utc.cooling_temp import calculate_cooling_temperature +from pythermogis.workflow.utc.doublet_utils import ( + calc_lifetime, + calculate_injection_temp_with_heat_pump, +) +from pythermogis.workflow.utc.pressure import calculate_max_pressure, optimize_pressure +from pythermogis.workflow.utc.utc_properties import UTCConfiguration +from pythermogis.workflow.utc.well_distance import optimize_well_distance + +EUR_PER_CT_PER_KWH = 0.36 +NPV_SCALE = 1e-6 + + +@dataclass(slots=True, frozen=True) +class DoubletInput: + unknown_input_value: float + thickness: float + transmissivity: float + transmissivity_with_ntg: float + ntg: float + depth: float + porosity: float + temperature: float + + @property + def permeability(self) -> float: + return self.transmissivity / self.thickness * 1e-3 * DARCY_SI + + +class DoubletOutput(NamedTuple): + power: float + hppower: float + capex: float + var_opex: float + fixed_opex: float + opex: float + utc: float + npv: float + hprod: float + cop: float + cophp: float + pres: float + flow: float + welld: float + breakthrough: float + cooling: float + production_temp: float + injection_temp: float + + +def calculate_doublet_performance( + props: UTCConfiguration, input: DoubletInput, verbose: bool = False +) -> DoubletOutput | None: + timer = timeit.default_timer() + # determine initial well distance + well_distance = ( + (props.optim_dist_well_dist_min + props.optim_dist_well_dist_max) / 2 + if props.optim_well_dist + else props.default_well_distance + ) + + # determine injection temperature + if props.use_heat_pump: + injection_temperature = calculate_injection_temp_with_heat_pump( + input.temperature, + props.hp_direct_heat_input_temp, + props.dh_return_temp, + input.temperature, + props.max_cooling_temp_range, + props.hp_minimum_injection_temperature, + ) + else: + injection_temperature = max( + input.temperature - props.max_cooling_temp_range, props.dh_return_temp + ) + timer = print_time(timer, "injection temperature: ", verbose=verbose) + + # calculate maximum pressure + drawdown_pressure = calculate_max_pressure( + props, + input, + False, + well_distance, + injection_temperature, + ) + if drawdown_pressure == 0: + return None + timer = print_time(timer, "maximum pressure: ", verbose=verbose) + + # cooling temperature and well distance optimization + if props.optim_well_dist: + end_temperature_p = props.optim_dist_cooling_fraction * ( + input.temperature - injection_temperature + ) + else: + end_temperature_p = calculate_cooling_temperature( + props, + input, + drawdown_pressure, + well_distance, + injection_temperature, + ) + timer = print_time(timer, "cooling temperature: ", verbose=verbose) + + if props.optim_well_dist: + well_distance = optimize_well_distance( + props, + input, + drawdown_pressure, + injection_temperature, + ) + timer = print_time(timer, "well distance optimizer: ", verbose=verbose) + + # stimulation capex + stimulation_capex = ( + 0.0 + if ( + not props.use_stimulation + or input.transmissivity_with_ntg > props.stim_kh_max + ) + else props.stimulation_capex + ) + + # pressure optimizer + pressure_results = optimize_pressure( + props=props, + input=input, + drawdown_pressure=drawdown_pressure, + well_distance=well_distance, + injection_temp=injection_temperature, + stimulation_capex=stimulation_capex, + ) + timer = print_time(timer, "pressure optimizer: ", verbose=verbose) + if pressure_results is None: + return None + + total_variable_opex = sum(pressure_results.variable_opex) + total_fixed_opex = sum(pressure_results.fixed_opex) + + utc_eur_ct_per_kwh = ( + input.unknown_input_value + if pressure_results.utc == input.unknown_input_value + else pressure_results.utc * EUR_PER_CT_PER_KWH + ) + + if abs(input.unknown_input_value) < 10: + if not (abs(utc_eur_ct_per_kwh - input.unknown_input_value) > 1e-4): + return None + else: + test = abs((utc_eur_ct_per_kwh / input.unknown_input_value) - 1) + if not (test > 1e-4): + return None + + breakthrough_years = calc_lifetime( + well_distance, + input.thickness * input.ntg, + 0.001, + input.porosity, + pressure_results.flowrate, + input.depth, + input.temperature, + props.salinity_surface, + props.salinity_gradient, + props.optim_dist_cp_rock, + props.optim_dist_rho_rock, + ) + + utc_cutoff = ( + props.utc_cutoff_deep + if input.depth > props.utc_deep_depth + else props.utc_cutoff + ) + + npv = ( + NPV_SCALE + * (utc_cutoff - utc_eur_ct_per_kwh) + * 3.6 + * pressure_results.discounted_heat_produced + * (1 - props.tax_rate) + ) + opex_first_prod_year = pressure_results.total_opex_ts[props.drilling_time] + print_time(timer, "economics: ", verbose=verbose) + + return DoubletOutput( + power=pressure_results.heat_power_per_doublet, + hppower=pressure_results.hp_added_power, + capex=pressure_results.sum_capex, + var_opex=total_variable_opex, + fixed_opex=total_fixed_opex, + opex=opex_first_prod_year, + utc=utc_eur_ct_per_kwh, + npv=npv, + hprod=pressure_results.discounted_heat_produced, + cop=pressure_results.cop, + cophp=pressure_results.heat_pump_cop, + pres=drawdown_pressure / 1e5, + flow=pressure_results.flowrate, + welld=well_distance, + breakthrough=breakthrough_years, + cooling=end_temperature_p, + production_temp=pressure_results.production_temp, + injection_temp=injection_temperature, + ) diff --git a/src/pythermogis/workflow/utc/doublet_utils.py b/src/pythermogis/workflow/utc/doublet_utils.py new file mode 100644 index 0000000000000000000000000000000000000000..bf3200da7f9ebdccb0447d21d9a3bbf508c476fa --- /dev/null +++ b/src/pythermogis/workflow/utc/doublet_utils.py @@ -0,0 +1,151 @@ +import math + +from numba import njit + +from pythermogis.workflow.utc.water import ( + density, + get_hydrostatic_pressure, + get_salinity, + heat_capacity, +) + + +@njit +def calculate_injection_temp_with_heat_pump( + t_prod: float, + t_goal: float, + dh_return_temp: float, + reservoir_temp: float, + max_cooling_temp_range: float, + hp_minimum_injection_temperature: float, +): + t_kelvin = 273.15 + + eta = 0.5 + delta_t = 3 # deg C + + tc = t_goal + t_kelvin + tr = t_prod + t_kelvin + t_waste = min(dh_return_temp, t_prod) + t_kelvin + + num = ( + -eta * (tc + delta_t) * (t_waste + tr - tc) + - tc * tc + + tr * tc + - 2 * delta_t * (tc - tr) + ) + + denom = tr - tc - eta * (tc + delta_t) + + needed_injection_temp = num / denom - t_kelvin + + return max( + needed_injection_temp, + max(reservoir_temp - max_cooling_temp_range, hp_minimum_injection_temperature), + ) + + +@njit +def get_orc_efficiency(Tx: float, Ts: float, etarel: float) -> float: + return etarel * (Tx - Ts) / (Tx + Ts + 2 * 273.1) + + +@njit +def get_cop_carnot(eta: float, Tout: float, Tin: float) -> float: + if Tout < 0.0 or Tin < 0.0: + return 0.0 + + if eta < 0.0 or eta > 1.0: + return 0.0 + + if Tout < Tin: + return 0.0 + + DHP = 3.0 + TKELVIN = 273.15 + + Tcond = Tout + DHP + Tevap = Tin - DHP + + return eta * (Tcond + TKELVIN) / (Tcond - Tevap) + + +@njit +def calc_lifetime( + well_distance: float, + thickness: float, + delta_temp_fraction: float, + porosity: float, + flowrate: float, + depth: float, + reservoir_temp: float, + salinity_surface: float, + salinity_gradient: float, + cp_rock: float, + rho_rock: float, +) -> float: + # --- Water properties at depth --- + salinity = get_salinity(salinity_surface, salinity_gradient, depth) + pressure = get_hydrostatic_pressure(depth) + rho_water = density(pressure, reservoir_temp, salinity) + cp_water = heat_capacity(pressure, reservoir_temp, salinity) + + dangle = 0.05 * math.pi / 180.0 # radians + halfdist = well_distance * 0.5 + + angleseg = 0.5 * (2 * math.pi * delta_temp_fraction) + 0.5 * dangle + radius = halfdist / math.sin(angleseg) + + Aseg1 = 4 * ( + 0.5 * angleseg * radius**2 + - 0.5 * radius**2 * math.sin(angleseg) * math.cos(angleseg) + ) + + angleseg = 0.5 * (2 * math.pi * delta_temp_fraction) - 0.5 * dangle + radius = halfdist / math.sin(angleseg) + + Aseg2 = 4 * ( + 0.5 * angleseg * radius**2 + - 0.5 * radius**2 * math.sin(angleseg) * math.cos(angleseg) + ) + + Aseg = Aseg1 - Aseg2 + Vseg = Aseg * thickness + + Eseg = ( + 1e-9 + * Vseg + * ((1 - porosity) * cp_rock * rho_rock + porosity * cp_water * rho_water) + ) + + flowseg = flowrate * 365 * 24 * dangle / math.pi + Eflowseg = 1e-9 * flowseg * cp_water * rho_water + + return Eseg / Eflowseg + + +@njit +def get_along_hole_length( + true_vertical_depth: float, + well_distance: float, + curve_scaling_factor: float, + max_true_vertical_depth_stepout_factor: float, +) -> float: + if curve_scaling_factor == 0: + return true_vertical_depth + + stepout = 0.5 * well_distance + + max_stepout = true_vertical_depth * max_true_vertical_depth_stepout_factor + + horizontal_distance = stepout - max_stepout + if horizontal_distance < 0: + horizontal_distance = 0.0 + + stepout -= horizontal_distance + + oblique_distance = ( + math.sqrt(stepout**2 + true_vertical_depth**2) * curve_scaling_factor + ) + + return oblique_distance + horizontal_distance diff --git a/src/pythermogis/workflow/utc/doubletcalc.py b/src/pythermogis/workflow/utc/doubletcalc.py new file mode 100644 index 0000000000000000000000000000000000000000..78fc61f211e1b64e96e389ccde13ae731c51672b --- /dev/null +++ b/src/pythermogis/workflow/utc/doubletcalc.py @@ -0,0 +1,193 @@ +from __future__ import annotations + +import math +from typing import TYPE_CHECKING, NamedTuple + +from numba import njit +from pydoubletcalc import Aquifer, Doublet, Well, WellPipeSegment + +from pythermogis.workflow.utc.rock import get_geothermal_gradient +from pythermogis.workflow.utc.water import get_salinity + +if TYPE_CHECKING: + from pythermogis.workflow.utc.doublet import DoubletInput + from pythermogis.workflow.utc.utc_properties import UTCConfiguration + + +INCH_SI = 0.0254 + + +class Doublet1DResults(NamedTuple): + geothermal_powers: float + cop: float + flowrate: float + pump_power_required: float + production_temp: float + heat_power_produced: list[float] + + +def doubletcalc( + props: UTCConfiguration, + input: DoubletInput, + drawdown_pressure: float, + well_distance: float, + injection_temp: float, +) -> Doublet1DResults | None: + aquifer = Aquifer( + permeability=input.permeability, + porosity=input.porosity, + net_to_gross=input.ntg, + thickness=input.thickness, + salinity=get_salinity( + props.salinity_surface, props.salinity_gradient, input.depth + ), + surface_temperature=props.surface_temperature, + geothermal_gradient=get_geothermal_gradient( + input.depth, input.thickness, input.temperature, props.surface_temperature + ), + rock_heat_capacity=855, + rock_density=2720, + thermal_conductivity=3, + thermal_diffusivity=1.2e-6, + ) + + injector = Well( + aquifer=aquifer, + well_type="injector", + pipe_segments=[ + WellPipeSegment( + ah_depth=get_along_hole_length( + input.depth, + well_distance, + props.well_curv_scaling, + props.max_tvd_stepout_factor, + ), + tv_depth=input.depth, + inner_diameter=props.inner_diameter * INCH_SI, + roughness=props.roughness * 1e-3 * INCH_SI, + ) + ], + aquifer_top_depth=input.depth, + pipe_scaling=0, + target_segment_length=props.segment_length, + outer_diameter=props.outer_diameter * INCH_SI, + skin=get_total_skin_injection(props, input), + pump_present=False, + pump_pressure=0, + pump_depth=0, + pump_efficiency=0, + pressure_balance_tolerance=0.5, + viscosity_mode=props.viscosity_mode, + heat_exchanger_exit_temperature=injection_temp, + ) + + producer = Well( + aquifer=aquifer, + well_type="producer", + pipe_segments=[ + WellPipeSegment( + ah_depth=get_along_hole_length( + input.depth, + well_distance, + props.well_curv_scaling, + props.max_tvd_stepout_factor, + ), + tv_depth=input.depth, + inner_diameter=props.inner_diameter * INCH_SI, + roughness=props.roughness * 1e-3 * INCH_SI, + ) + ], + aquifer_top_depth=input.depth, + pipe_scaling=0, + target_segment_length=props.segment_length, + outer_diameter=props.outer_diameter * INCH_SI, + skin=get_total_skin_production(props, input), + pump_present=True, + pump_pressure=drawdown_pressure, + pump_depth=get_pump_production_depth(props, input.depth), + pump_efficiency=props.pump_efficiency, + pressure_balance_tolerance=0.5, + viscosity_mode=props.viscosity_mode, + heat_exchanger_exit_temperature=injection_temp, + ) + doublet = Doublet( + aquifer=aquifer, + injector=injector, + producer=producer, + well_distance=well_distance, + cooling_fraction=0.1, + yearly_operating_hours=props.load_hours, + pressure_balance_tolerance=1e-1, + viscosity_mode=props.viscosity_mode, + heat_exchanger_exit_temperature=injection_temp, + ) + + # don't calculate lifetime as it is unneded for Doublet1DResults + doublet.calculate_mass_rate() + if not doublet.simulation_success: + return None + doublet.calculate_geothermal_power() + doublet.calculate_power() + + return Doublet1DResults( + geothermal_powers=doublet.geothermal_power, + cop=doublet.cop, + flowrate=doublet.flowrate, + pump_power_required=doublet.pump_power, + production_temp=doublet.production_temp, + heat_power_produced=[ + doublet.geothermal_power, + ] + * props.econ_lifetime_years, + ) + + +def get_total_skin_injection(props: UTCConfiguration, input: DoubletInput) -> float: + if (not props.use_stimulation) or ( + input.transmissivity_with_ntg > props.stim_kh_max + ): + stim_add_skin_inj = 0.0 + else: + stim_add_skin_inj = props.stim_add_skin_inj + + return props.skin_injector + stim_add_skin_inj + + +def get_total_skin_production(props: UTCConfiguration, input: DoubletInput): + if (not props.use_stimulation) or ( + input.transmissivity_with_ntg > props.stim_kh_max + ): + stim_add_skin_prod = 0.0 + else: + stim_add_skin_prod = props.stim_add_skin_prod + + return props.skin_producer + stim_add_skin_prod + + +def get_pump_production_depth(props: UTCConfiguration, depth: float) -> float: + return min(props.pump_depth, depth / 2) + + +@njit +def get_along_hole_length( + true_vertical_depth: float, + well_distance: float, + curve_scaling_factor: float, + max_true_vertical_depth_stepout_factor: float, +) -> float: + if curve_scaling_factor == 0: + # Vertical well: along-hole length equals TVD + return true_vertical_depth + + # Horizontal distance from drilling location + stepout = 0.5 * well_distance + + max_stepout = true_vertical_depth * max_true_vertical_depth_stepout_factor + horizontal_distance = stepout - max_stepout if (stepout - max_stepout) > 0 else 0 + + stepout -= horizontal_distance + + oblique_distance = ( + math.sqrt(stepout**2 + true_vertical_depth**2) * curve_scaling_factor + ) + return oblique_distance + horizontal_distance diff --git a/src/pythermogis/workflow/utc/economics.py b/src/pythermogis/workflow/utc/economics.py new file mode 100644 index 0000000000000000000000000000000000000000..de4329805314623fad11a7aed98ff9de99644f81 --- /dev/null +++ b/src/pythermogis/workflow/utc/economics.py @@ -0,0 +1,498 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING, NamedTuple + +import numpy as np +from numba import njit +from numpy.typing import NDArray + +from pythermogis.workflow.utc.doublet_utils import get_along_hole_length + +if TYPE_CHECKING: + from pythermogis.workflow.utc.doublet import DoubletInput + from pythermogis.workflow.utc.utc_properties import UTCConfiguration + +SECONDS_IN_YEAR = 365 * 24 * 3600 +MJ_TO_GJ = 1e-3 +KWH_TO_MJ = 0.36 * 1e9 +HOURS_IN_YEAR = 8760 + + +class CapexCalculatorResults(NamedTuple): + sum_capex: float + total_capex: float + variable_opex: list[float] + fixed_opex: list[float] + total_opex_ts: list[float] + heat_power_per_year: list[float] + + +class UTCCalculatorResults(NamedTuple): + discounted_heat_produced: float + utc: float + + +class EconomicsResults(NamedTuple): + capex: CapexCalculatorResults + utc: UTCCalculatorResults + + +def calculate_economics( + props: UTCConfiguration, + input: DoubletInput, + well_distance: float, + heat_power_produced: list[float], + stimulation_capex: float, + hp_cop: float, + cop_hp_years: list[float] | None, + season_factor_years: list[float] | None, + hp_added_power_years: list[float] | None, + hp_added_power: float, + pump_power_required: float, +) -> EconomicsResults: + ah_length_array = [ + get_along_hole_length( + true_vertical_depth=input.depth, + well_distance=well_distance, + curve_scaling_factor=props.well_curv_scaling, + max_true_vertical_depth_stepout_factor=props.max_tvd_stepout_factor, + ) + ] + + capex_results = calculate_capex( + props=props, + heat_power_produced=heat_power_produced, + stimulation_capex=stimulation_capex, + hp_cop=hp_cop, + hp_cop_years=cop_hp_years, + season_factor_years=season_factor_years, + hp_added_power_years=hp_added_power_years, + initial_hp_added_power=hp_added_power, + ah_length_array=ah_length_array, + pump_power_required=pump_power_required, + ) + + utc_results = calculate_utc( + props=props, + heat_power_per_year=capex_results.heat_power_per_year, + total_opex_ts=capex_results.total_opex_ts, + total_capex=capex_results.total_capex, + ) + + return EconomicsResults( + capex=capex_results, + utc=utc_results, + ) + + +def calculate_capex( + props: UTCConfiguration, + heat_power_produced: list[float], + stimulation_capex: float, + hp_cop: float, + hp_cop_years: list[float] | None, + season_factor_years: list[float] | None, + hp_added_power_years: list[float] | None, + initial_hp_added_power: float, + ah_length_array: list[float], + pump_power_required: float, +) -> CapexCalculatorResults: + n_years = props.economic_lifetime + + heat_power_per_year = np.zeros(n_years) + variable_opex = np.zeros(n_years) + fixed_opex = np.zeros(n_years) + total_opex_ts = np.zeros(n_years) + total_capex_ts = np.zeros(n_years) + + shifted_heat_power = np.array(heat_power_produced, dtype=np.float64) + shift_time_series(shifted_heat_power, props.drilling_time) + + capex_well = calculate_capex_for_wells( + inj_depth=ah_length_array[0], + prod_depth=ah_length_array[0], + stimulation_capex=stimulation_capex, + well_cost_z2=props.well_cost_z2, + well_cost_z=props.well_cost_z, + well_cost_const=props.well_cost_const, + well_cost_scaling=props.well_cost_scaling, + ) + allocate_capex_for_wells(props.drilling_time, total_capex_ts, capex_well) + + hp_years_arr = ( + np.array(hp_added_power_years, dtype=np.float64) + if hp_added_power_years is not None + else np.zeros(0, dtype=np.float64) + ) + + hp_added_power = get_max_hp_added_power(initial_hp_added_power, hp_years_arr) + hp_after_hp = calculate_hp_added_power_after_heat_pump( + hp_added_power, hp_cop, props.hp_alternative_heating_price + ) + installation_mw = calculate_installation_mw( + float(shifted_heat_power[-1]), hp_added_power + ) + + total_capex_1year = calculate_total_capex_1year( + last_year=max(props.drilling_time - 1, 0), + total_capex_ts=total_capex_ts, + capex_const=props.capex_const, + hp_capex=props.hp_capex, + capex_variable=props.capex_variable, + capex_contingency=props.capex_contingency, + hp_after_hp=hp_after_hp, + installation_mw=installation_mw, + ) + + sum_capex = process_annual_data( + props, + shifted_heat_power, + total_capex_ts, + heat_power_per_year, + variable_opex, + fixed_opex, + total_opex_ts, + total_capex_1year, + installation_mw, + hp_after_hp, + hp_cop, + hp_cop_years, + season_factor_years, + pump_power_required, + ) + + total_capex_sum = float(np.sum(total_capex_ts)) + + return CapexCalculatorResults( + sum_capex=float(sum_capex), + total_capex=float(total_capex_sum), + variable_opex=list(variable_opex), + fixed_opex=list(fixed_opex), + total_opex_ts=list(total_opex_ts), + heat_power_per_year=list(heat_power_per_year), + ) + + +def shift_time_series(series: list[float], shift: int) -> None: + n = len(series) + if shift <= 0 or shift >= n: + for i in range(n): + series[i] = 0 + else: + series[shift:] = series[: n - shift] + for i in range(shift): + series[i] = 0 + + +@njit +def calculate_capex_for_wells( + inj_depth: float, + prod_depth: float, + stimulation_capex: float, + well_cost_z2: float, + well_cost_z: float, + well_cost_const: float, + well_cost_scaling: float, +) -> float: + capex = ( + calculate_well_cost(well_cost_z2, well_cost_z, well_cost_const, inj_depth) + + calculate_well_cost(well_cost_z2, well_cost_z, well_cost_const, prod_depth) + + stimulation_capex + ) + return capex * well_cost_scaling + + +@njit +def calculate_well_cost( + well_cost_z2: float, well_cost_z: float, well_cost_const: float, depth: float +) -> float: + return (well_cost_z2 * depth**2 + well_cost_z * depth) * 1e-6 + well_cost_const + + +@njit +def allocate_capex_for_wells( + drilling_years: int, total_capex_ts: np.ndarray, capex_well: float +) -> None: + yearly_cost = capex_well / max(drilling_years, 1) + for y in range(drilling_years): + total_capex_ts[y] += yearly_cost + + +@njit +def get_max_hp_added_power(initial: float, hp_power_years: np.ndarray) -> float: + if hp_power_years.size == 0: + return initial + max_val = hp_power_years[0] + for i in range(hp_power_years.size): + if hp_power_years[i] > max_val: + max_val = hp_power_years[i] + return max_val + + +@njit +def calculate_hp_added_power_after_heat_pump( + hp_added_power: float, hp_cop: float, hp_alternative_heating_price: float +) -> float: + if hp_alternative_heating_price < 0: + return hp_added_power + return hp_added_power * (1 + 1.0 / (hp_cop - 1.0)) + + +@njit +def calculate_installation_mw( + last_heat_power_value: float, hp_added_power: float +) -> float: + diff = last_heat_power_value - hp_added_power + return diff if diff > 0 else 0.0 + + +@njit +def calculate_total_capex_1year( + last_year: int, + total_capex_ts: np.ndarray, + capex_const: float, + hp_capex: float, + capex_variable: float, + capex_contingency: float, + hp_after_hp: float, + installation_mw: float, +) -> float: + total_capex_ts[last_year] += ( + capex_const + (hp_capex * hp_after_hp + capex_variable * installation_mw) * 1e-3 + ) + + total_capex_ts[last_year] *= 1.0 + capex_contingency / 100.0 + return total_capex_ts[last_year] + + +def process_annual_data( + props: UTCConfiguration, + heat_power_produced: list[float], + total_capex_ts: NDArray[np.float64], + heat_power_per_year: NDArray[np.float64], + variable_opex: NDArray[np.float64], + fixed_opex: NDArray[np.float64], + total_opex_ts: NDArray[np.float64], + total_capex_1year: float, + installation_mw: float, + hp_after_hp: float, + hp_cop: float, + hp_cop_years: float, + season_factor_years: float, + pump_power_required: float, +) -> float: + sum_capex = 0.0 + + for year in range(props.economic_lifetime): + # Season factor + season_factor = ( + season_factor_years[year] + if season_factor_years is not None + else get_heat_exchanger_season_factor(props) + ) + + # COP for that year + year_hp_cop = hp_cop_years[year] if hp_cop_years is not None else hp_cop + + # Annual heat (GJ/yr) + heat_power_per_year[year] = ( + heat_power_produced[year] * season_factor * SECONDS_IN_YEAR * MJ_TO_GJ + ) + + sum_capex += total_capex_ts[year] + total_opex_ts[year] = 0 + + # Skip drilling + negative production years + if year < props.drilling_time or heat_power_produced[year] < 0: + continue + + # Variable OPEX + variable_opex[year] = calculate_variable_opex( + season_factor, + props.elec_purchase_price, + pump_power_required, + heat_power_produced[year], + props.heat_exchanger_parasitic, + heat_power_per_year[year], + year_hp_cop, + hp_after_hp, + props.opex_per_energy, + ) + # Fixed OPEX + fixed_opex[year] = calculate_fixed_opex( + props.opex_base, + props.hp_opex, + props.opex_per_power, + hp_after_hp, + installation_mw, + total_capex_1year, + props.opex_per_capex, + ) + + total_opex_ts[year] = variable_opex[year] + fixed_opex[year] + + return sum_capex + + +@njit +def calculate_variable_opex( + season_factor: float, + elec_price: float, + pump_power: float, + heat_power_produced: float, + heat_exchanger_parasitic: float, + heat_power_gj_per_year: float, + hp_cop: float, + hp_after_hp: float, + opex_per_energy: float, +) -> float: + pump_cost = ( + -(pump_power * 1e3) + * SECONDS_IN_YEAR + * season_factor + * elec_price + / KWH_TO_MJ + * 1e-6 + ) + + parasitic_cost = ( + -(heat_power_produced * heat_exchanger_parasitic * 1e6) + * SECONDS_IN_YEAR + * season_factor + * elec_price + / KWH_TO_MJ + * 1e-6 + ) + + opex_energy_cost = ( + -(1 / 0.36) * opex_per_energy * 1e-2 * heat_power_gj_per_year * 1e4 / 36 * 1e-6 + ) + + heat_pump_cost = 0.0 + if hp_cop != 0: + heat_pump_cost = ( + -(hp_after_hp * 1e6 / hp_cop) + * SECONDS_IN_YEAR + * season_factor + * elec_price + / KWH_TO_MJ + * 1e-6 + ) + + return pump_cost + parasitic_cost + heat_pump_cost + opex_energy_cost + + +@njit +def calculate_fixed_opex( + opex_base: float, + hp_opex: float, + opex_per_power: float, + hp_after_hp: float, + installation_mw: float, + total_capex: float, + opex_per_capex: float, +) -> float: + return ( + -opex_base / 1e6 + - (hp_opex * hp_after_hp + opex_per_power * installation_mw) * 1e-3 + - total_capex * opex_per_capex / 100.0 + ) + + +def get_heat_exchanger_season_factor(props: UTCConfiguration) -> float: + return props.load_hours / HOURS_IN_YEAR + + +def calculate_utc( + props: UTCConfiguration, + heat_power_per_year: float, + total_opex_ts: float, + total_capex: float, +) -> UTCCalculatorResults: + net_cash_worth_eia = calculate_net_cash_worth_eia(props, total_capex) + present_value = total_capex * props.debt_equity - net_cash_worth_eia + yearly_payment = calculate_payment_value(props, present_value) + depreciation_cost = calculate_depreciation_cost(props, total_capex) + + discounted_income = 0.0 + discounted_heat_produced = 0.0 + + for year in range(1, props.economic_lifetime): + inflated_opex = total_opex_ts[year] * ((1 + props.inflation) ** (year - 1)) + + future_value = calculate_future_value( + present_value, yearly_payment, props.interest_loan, year + ) + interest_payment = -(future_value * props.interest_loan) + + taxable_expenses = inflated_opex + depreciation_cost + interest_payment + tax_deduction = -(taxable_expenses * props.tax_rate) + + net_income = inflated_opex + yearly_payment + tax_deduction + + yearly_energy = heat_power_per_year[year] * (1 - props.tax_rate) + + discounted_income += net_income / ((1 + props.equity_return) ** year) + discounted_heat_produced += yearly_energy / ((1 + props.equity_return) ** year) + + utc = calculate_unit_technical_cost( + props, total_capex, discounted_income, discounted_heat_produced + ) + + if discounted_heat_produced < 0 or utc < 0: + raise ValueError( + f"discountedHeatProduced and unitTechnicalCost must be positive. " + f"Values were: {discounted_heat_produced}, {utc}" + ) + + return UTCCalculatorResults( + discounted_heat_produced=discounted_heat_produced, + utc=utc, + ) + + +def calculate_net_cash_worth_eia(props: UTCConfiguration, total_capex: float) -> float: + tax_rate = props.tax_rate + project_interest_rate = calculate_project_interest_rate(props) + return (tax_rate * (props.ecn_eia * total_capex)) / (1 + project_interest_rate) + + +def calculate_payment_value(props: UTCConfiguration, present_value: float) -> float: + if props.interest_loan == 0: + return 0.0 + + factor = (1 + props.interest_loan) ** props.economic_lifetime + return (props.interest_loan / (factor - 1)) * (-(present_value * factor)) + + +def calculate_depreciation_cost(props: UTCConfiguration, total_capex: float) -> float: + return -(total_capex / props.economic_lifetime) + + +def calculate_future_value( + present_value: float, payment: float, interest: float, year: float +) -> float: + future_value = present_value * ((1 + interest) ** (year - 1)) + if payment != 0: + future_value += payment * (((1 + interest) ** (year - 1)) - 1) / interest + return future_value + + +def calculate_unit_technical_cost( + props: UTCConfiguration, + total_capex: float, + discounted_income: float, + discounted_heat_produced: float, +) -> float: + if discounted_heat_produced <= 0: + return 0.0 + return ( + (total_capex * (1 - props.debt_equity) - discounted_income) + / discounted_heat_produced + ) * 1e6 + + +def calculate_project_interest_rate(props: UTCConfiguration) -> float: + return (1 - props.tax_rate) * props.interest_loan * props.debt_equity + ( + 1 - props.debt_equity + ) * props.equity_return diff --git a/src/pythermogis/workflow/utc/flow.py b/src/pythermogis/workflow/utc/flow.py new file mode 100644 index 0000000000000000000000000000000000000000..7fc51d20ea1c6c6cc53a36a3c78889dc9a7293e3 --- /dev/null +++ b/src/pythermogis/workflow/utc/flow.py @@ -0,0 +1,99 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING, NamedTuple + +import numpy as np + +from pythermogis.workflow.utc.doublet_utils import get_orc_efficiency +from pythermogis.workflow.utc.doubletcalc import doubletcalc +from pythermogis.workflow.utc.heatpump import calculate_heat_pump_performance + +if TYPE_CHECKING: + from pythermogis.workflow.utc.doublet import DoubletInput + from pythermogis.workflow.utc.utc_properties import UTCConfiguration + + +class VolumetricFlowResults(NamedTuple): + hp_cop: float + hp_added_power: float + hp_elec_consumption: float + heat_power_per_doublet: float + cop: float + flowrate: float + pump_power_required: float + production_temp: float + heat_power_produced: list[float] + + +def calculate_volumetric_flow( + props: UTCConfiguration, + input_data: DoubletInput, + original_pressure: float, + well_distance: float, + injection_temp: float, +) -> VolumetricFlowResults: + for step in [0, 1e5, -1e5, 2e5, -2e5, 3e5, -3e5]: + d1d_results = doubletcalc( + props, input_data, original_pressure + step, well_distance, injection_temp + ) + if d1d_results is not None: + break + + # if no valid output results return None + if d1d_results is None: + return None + + if props.use_orc: + heat_exchanger_efficiency = get_orc_efficiency( + d1d_results.production_temp, + props.heat_exchanger_basetemp, + props.heat_exchanger_efficiency, + ) + else: + heat_exchanger_efficiency = props.heat_exchanger_efficiency + + heat_power_produced = np.clip( + np.array(d1d_results.heat_power_produced) * heat_exchanger_efficiency, + 1e-6, + None, + ) + + # Calculate cop + heat_power_per_doublet = max( + 1e-6, d1d_results.geothermal_powers * heat_exchanger_efficiency + ) + power_consumption = ( + (heat_power_per_doublet / heat_exchanger_efficiency) / d1d_results.cop + ) + (heat_power_per_doublet * props.heat_exchanger_parasitic) + if props.use_orc: + heat_power_per_doublet -= power_consumption + cop = heat_power_per_doublet / power_consumption + + # Calculate heat pump performance + if props.use_heat_pump: + hp_results = calculate_heat_pump_performance( + props, + input_data, + d1d_results.production_temp, + injection_temp, + d1d_results.flowrate, + ) + hp_cop = hp_results.hp_cop + hp_added_power = hp_results.hp_added_power + hp_elec_consumption = hp_results.hp_elec_consumption + else: + hp_cop = 0.0 + hp_added_power = 0.0 + hp_elec_consumption = 0.0 + + return VolumetricFlowResults( + hp_cop, + hp_added_power, + hp_elec_consumption, + heat_power_per_doublet, + cop, + d1d_results.flowrate, + d1d_results.pump_power_required, + d1d_results.production_temp, + heat_power_produced, + ) diff --git a/src/pythermogis/workflow/utc/heatpump.py b/src/pythermogis/workflow/utc/heatpump.py new file mode 100644 index 0000000000000000000000000000000000000000..b2d7a776fa8bbd7932f6c5a27895842a10633d89 --- /dev/null +++ b/src/pythermogis/workflow/utc/heatpump.py @@ -0,0 +1,88 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING, NamedTuple + +from pythermogis.workflow.utc.doublet_utils import get_cop_carnot +from pythermogis.workflow.utc.water import ( + density, + get_hydrostatic_pressure, + get_salinity, + heat_capacity, +) + +if TYPE_CHECKING: + from pythermogis.workflow.utc.doublet import DoubletInput + from pythermogis.workflow.utc.utc_properties import UTCConfiguration + + +class HeatPumpPerformanceResults(NamedTuple): + hp_cop: float + hp_added_power: float + hp_elec_consumption: float + + +def calculate_heat_pump_performance( + props: UTCConfiguration, + input_data: DoubletInput, + production_temp: float, + injection_temp: float, + flowrate: float, +) -> HeatPumpPerformanceResults: + ETA_CARNOT = 0.6 + + heat_pump_start_temp = calculate_heat_pump_start_temp( + props, production_temp, injection_temp + ) + + hydrostatic_pressure = get_hydrostatic_pressure(0) + + salinity = get_salinity( + props.salinity_surface, + props.salinity_gradient, + input_data.depth, + ) + + rho_water = density(hydrostatic_pressure, props.surface_temperature, salinity) + + cp_water = heat_capacity(hydrostatic_pressure, props.surface_temperature, salinity) + + # calculate the power used by the heat pump + waste_temp = min(heat_pump_start_temp, props.dh_return_temp) + delta_temp = max(0.0, waste_temp - injection_temp) + hp_added_power = (flowrate / 3600.0) * delta_temp * rho_water * cp_water * 1e-6 + if hp_added_power == 0: + return HeatPumpPerformanceResults( + hp_cop=0, + hp_added_power=0, + hp_elec_consumption=0, + ) + + # calculate heat pump coefficient of performance + Tout = props.hp_direct_heat_input_temp + Tin = min(injection_temp, props.dh_return_temp) + hp_cop = get_cop_carnot(ETA_CARNOT, Tout, Tin) + + # calculate hp_elec_consumption + hp_elec_consumption = hp_added_power / (hp_cop - 1) + + return HeatPumpPerformanceResults( + hp_cop=hp_cop, + hp_added_power=hp_added_power, + hp_elec_consumption=hp_elec_consumption, + ) + + +def calculate_heat_pump_start_temp( + props: UTCConfiguration, production_temp: float, injection_temp: float +) -> float: + """ + Python version of calculateHeatPumpStartTemp(...) + """ + delta_temp_geothermal = production_temp - injection_temp + delta_temp_for_hex = 0.0 + + if production_temp > props.dh_return_temp: + delta_temp_for_hex = production_temp - props.dh_return_temp + delta_temp_for_hex = max(0.0, min(delta_temp_for_hex, delta_temp_geothermal)) + + return production_temp - delta_temp_for_hex diff --git a/src/pythermogis/workflow/utc/pressure.py b/src/pythermogis/workflow/utc/pressure.py new file mode 100644 index 0000000000000000000000000000000000000000..244bab4d3ea0b1399ba07be67f5a7152c3d53218 --- /dev/null +++ b/src/pythermogis/workflow/utc/pressure.py @@ -0,0 +1,304 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING, NamedTuple + +from pythermogis.workflow.utc.economics import calculate_economics +from pythermogis.workflow.utc.flow import calculate_volumetric_flow + +if TYPE_CHECKING: + from pythermogis.workflow.utc.doublet import DoubletInput + from pythermogis.workflow.utc.utc_properties import UTCConfiguration + + +def calculate_max_pressure( + props: UTCConfiguration, + input_data: DoubletInput, + use_olsthoorn_max_pressure: bool, + well_distance: float, + injection_temp: float, +) -> float: + if use_olsthoorn_max_pressure: + max_pres = 0.2 * 0.1 * input_data.depth * 100000 + else: + max_pres = input_data.depth * (0.135 - props.hy_gradient) * 100000 + + pressure = min(props.max_pump * 1e5, max_pres) + + results = calculate_volumetric_flow( + props, input_data, pressure, well_distance, injection_temp + ) + + if results is None: + pres_tol = 1e3 + iter_count = 0 + + pres_min = 1e5 + pres_max = pressure + + while iter_count < 1000 and abs(pres_max - pres_min) > pres_tol: + iter_count += 1 + pressure = 0.5 * (pres_min + pres_max) + + results = calculate_volumetric_flow( + props, input_data, pressure, well_distance, injection_temp + ) + + if results is not None: + pres_min = pressure + else: + pres_max = pressure + + if results is None: + pressure -= pres_tol + results = calculate_volumetric_flow( + props, input_data, pressure, well_distance, injection_temp + ) + + if results is None or iter_count >= 1000: + return 0.0 + + return pressure + + if results.heat_power_per_doublet < 0 and not props.is_ates: + pressure /= 2.0 + + results = calculate_volumetric_flow( + props, input_data, pressure, well_distance, injection_temp + ) + + if results.heatPowerPerDoublet() < 0: + return 0.0 + else: + pressure *= 2.0 + + return pressure + + +class PressureOptimizerResults(NamedTuple): + drawdown_pressure: float + flowrate: float + heat_power_per_doublet: float + heat_pump_cop: float + cop: float + utc: float + sum_capex: float + variable_opex: list[float] + fixed_opex: list[float] + total_opex_ts: list[float] + discounted_heat_produced: float + hp_added_power: float + production_temp: float + + +def optimize_pressure( + props: UTCConfiguration, + input: DoubletInput, + drawdown_pressure: float, + well_distance: float, + injection_temp: float, + stimulation_capex: float, +) -> PressureOptimizerResults: + pres_tol = 1e-4 * drawdown_pressure + pres_step = 1e1 * pres_tol + + flow_results = calculate_volumetric_flow( + props=props, + input_data=input, + original_pressure=drawdown_pressure, + well_distance=well_distance, + injection_temp=injection_temp, + ) + if flow_results is None: + return None + + if flow_results.production_temp > injection_temp: + iter_count = 0 + pres_min = 1e5 * props.min_pump + pres_max = drawdown_pressure + + # If slope at initial pressure is positive → begin optimization + if ( + utc_slope( + props, + input, + pres_max, + pres_step, + well_distance, + injection_temp, + stimulation_capex, + ) + > 0 + ): + while iter_count < 100 and abs(pres_max - pres_min) > pres_tol: + iter_count += 1 + pres = 0.5 * (pres_min + pres_max) + drawdown_pressure = pres + + flow_results = calculate_volumetric_flow( + props=props, + input_data=input, + original_pressure=drawdown_pressure, + well_distance=well_distance, + injection_temp=injection_temp, + ) + + if flow_results is not None: + utc_val = find_utc_for_pres( + props, + input, + pres, + well_distance, + injection_temp, + stimulation_capex, + ) + + slope = utc_slope( + props, + input, + pres, + pres_step, + well_distance, + injection_temp, + stimulation_capex, + ) + + if slope < utc_val * 1e-5 * 1e-2 * props.tolerance_utc_increase: + pres_min = pres + else: + pres_max = pres + + else: + # Slightly relax pressure if flow failed + pres_max = pres_max * 1.01 + + if flow_results is not None and flow_results.flowrate > props.max_flow: + iter_count = 0 + pres_min = 0.0 + pres_max = drawdown_pressure + while iter_count < 1000 and abs(pres_max - pres_min) > pres_tol: + iter_count += 1 + pres = 0.5 * (pres_min + pres_max) + drawdown_pressure = pres + + flow_results = calculate_volumetric_flow( + props=props, + input_data=input, + original_pressure=drawdown_pressure, + well_distance=well_distance, + injection_temp=injection_temp, + ) + + if flow_results.flowrate < props.max_flow: + pres_min = pres + else: + pres_max = pres + + # Final adjustment if flow still too high + if flow_results.flowrate > props.max_flow: + drawdown_pressure -= pres_tol + flow_results = calculate_volumetric_flow( + props=props, + input_data=input, + original_pressure=drawdown_pressure, + well_distance=well_distance, + injection_temp=injection_temp, + ) + + if flow_results is None or flow_results.flowrate > props.max_flow: + return None + + cop = (flow_results.heat_power_per_doublet + flow_results.hp_elec_consumption) / ( + flow_results.heat_power_per_doublet / flow_results.cop + + flow_results.hp_elec_consumption + ) + + econ = calculate_economics( + props=props, + input=input, + well_distance=well_distance, + heat_power_produced=flow_results.heat_power_produced, + stimulation_capex=stimulation_capex, + hp_cop=flow_results.hp_cop, + cop_hp_years=None, + season_factor_years=None, + hp_added_power_years=None, + hp_added_power=flow_results.hp_added_power, + pump_power_required=flow_results.pump_power_required, + ) + + return PressureOptimizerResults( + drawdown_pressure=drawdown_pressure, + flowrate=flow_results.flowrate, + heat_power_per_doublet=flow_results.heat_power_per_doublet, + heat_pump_cop=flow_results.hp_cop, + cop=cop, + utc=econ.utc.utc, + sum_capex=econ.capex.sum_capex, + variable_opex=econ.capex.variable_opex, + fixed_opex=econ.capex.fixed_opex, + total_opex_ts=econ.capex.total_opex_ts, + discounted_heat_produced=econ.utc.discounted_heat_produced, + hp_added_power=flow_results.hp_added_power, + production_temp=flow_results.production_temp, + ) + + +def utc_slope( + props: UTCConfiguration, + input: DoubletInput, + pres: float, + pres_step: float, + well_distance: float, + injection_temp: float, + stimulation_capex: float, +) -> float: + utc1 = find_utc_for_pres( + props, + input, + pres - pres_step / 2, + well_distance, + injection_temp, + stimulation_capex, + ) + utc2 = find_utc_for_pres( + props, + input, + pres + pres_step / 2, + well_distance, + injection_temp, + stimulation_capex, + ) + return (utc2 - utc1) / pres_step + + +def find_utc_for_pres( + props: UTCConfiguration, + input: DoubletInput, + pres: float, + well_distance: float, + injection_temp: float, + stimulation_capex: float, +) -> float: + flow_results = calculate_volumetric_flow( + props=props, + input_data=input, + original_pressure=pres, + well_distance=well_distance, + injection_temp=injection_temp, + ) + + econ = calculate_economics( + props=props, + input=input, + well_distance=well_distance, + heat_power_produced=flow_results.heat_power_produced, + stimulation_capex=stimulation_capex, + hp_cop=flow_results.hp_cop, + cop_hp_years=None, + season_factor_years=None, + hp_added_power_years=None, + hp_added_power=flow_results.hp_added_power, + pump_power_required=flow_results.pump_power_required, + ) + + return econ.utc.utc diff --git a/src/pythermogis/workflow/utc/rock.py b/src/pythermogis/workflow/utc/rock.py new file mode 100644 index 0000000000000000000000000000000000000000..6ebd2cf6d585bd30ee770cd60c8173e899e032f8 --- /dev/null +++ b/src/pythermogis/workflow/utc/rock.py @@ -0,0 +1,10 @@ +from numba import njit + + +@njit +def get_geothermal_gradient( + depth: float, thickness: float, temperature: float, surface_temperature: float +) -> float: + d_depth = depth + 0.5 * thickness + dt = temperature - surface_temperature + return dt / d_depth diff --git a/src/pythermogis/workflow/utc/utc_properties.py b/src/pythermogis/workflow/utc/utc_properties.py new file mode 100644 index 0000000000000000000000000000000000000000..b9f939956f915275bdb3d071101dea87f3b6e992 --- /dev/null +++ b/src/pythermogis/workflow/utc/utc_properties.py @@ -0,0 +1,196 @@ +from dataclasses import dataclass, field +from pathlib import Path +from typing import Literal, NamedTuple + + +class PropertyGridInfo(NamedTuple): + name: str + optional: bool + postfix: str + + +class AquiferFile(NamedTuple): + postfix: str + newPostfix: str + + +ViscosityMode = Literal["kestin", "batzlewang"] + + +@dataclass(slots=True, frozen=True) +class UTCConfiguration: + input_data_dir: str = "" + results_dir: str = "" + petrel_output_dir: Path | None = None + n_threads: int = 1 + is_ates: bool = False + check_copied_files: bool = True + validate_input_grids: bool = True + validate_output_grids: bool = True + aquifers: list[str] = field( + default_factory=lambda: [ + "NMVFS", + "NMVFV", + "NMRFT", + "NMRFV", + "NLFFS", + "NLFFD", + "NLLFR", + "NLLFS", + "KNGLG_KNGLS", + "KNNSG", + "KNNSL", + "KNNSY", + "KNNSB", + "KNNSR", + "KNNSF_KNNSP", + "SLDNA", + "SLDND", + "RNROF", + "RNSOB", + "RBMH", + "RBMDU", + "RBMDL", + "RBMVU", + "RBMVL", + "RBSHN", + "ROSL_ROSLU", + "ROSLL", + "DCH", + "DCD", + ] + ) + property_grid_infos: list[PropertyGridInfo] = field( + default_factory=lambda: [ + PropertyGridInfo("Permeability", False, "__k.zmap"), + PropertyGridInfo("PermeabilityLNSD", False, "__k_lnsd.zmap"), + PropertyGridInfo("Porosity", False, "__phi.zmap"), + PropertyGridInfo("Thickness", False, "__thick.zmap"), + PropertyGridInfo("ThicknessSD", False, "__thick_sd.zmap"), + PropertyGridInfo("Depth", False, "__top.zmap"), + PropertyGridInfo("NetToGross", False, "__ntg.zmap"), + PropertyGridInfo("Temperature", True, "__temperature.zmap"), + PropertyGridInfo("HCAccum", True, "__hc_accum.zmap"), + PropertyGridInfo("BoundaryShapefile", True, "__BoundaryShapefile.shp"), + ] + ) + copy_aquifer_files_info: list[AquiferFile] = field( + default_factory=lambda: [ + AquiferFile("__ntg_points.shp", "__ntg_points.shp"), + AquiferFile("__poro_points.shp", "__poro_points.shp"), + AquiferFile("__perm_points.shp", "__perm_points.shp"), + AquiferFile("__points_QC.shp", "__points_QC.shp"), + ] + ) + scenario: str = "basecase" + scen_suffix: str = "" + temp_from_grid: bool = False + exclude_hc_accum: bool = True + use_bounding_shape: bool = False + grid_ext: str = ".nc" + p_values: list[float] = field( + default_factory=lambda: [10.0, 30.0, 50.0, 70.0, 90.0] + ) + # temp_voxet: 'Voxet' = None + surface_temperature: float = 10.0 + temp_gradient: float = 31.0 + undefined_cells_allowed: int = 3 + level_of_detail: int = 4 + remove_padding_from_input_grids: bool = True + use_heat_pump: bool = False + dh_return_temp: float = 30.0 + ates_rosim_template_file: Path | None = None + scale_sd: float = 1.0 + calculate_cop: bool = True + hp_application_mode: bool = False + hp_direct_heat_input_temp: float = 70.0 + utc_cutoff: float = 5.1 + utc_cutoff_deep: float = 6.5 + utc_deep_depth: float = 4000.0 + ates_n_years_to_average: int = 1 + ates_rosim_years: int = 5 + maxdepth: float = 6000.0 + ates_min_depth: float = 0.0 + ates_max_depth: float = 500.0 + min_prod_temp: float = 20.0 + kh_cutoff: float = 1.0 + use_stimulation: bool = False + stim_kh_max: float = 20.0 + max_cooling_temp_range: float = 100.0 + hp_minimum_injection_temperature: float = 15.0 + ates_charge_temp: float = 80.0 + econ_lifetime_years: int = 15 + ates_min_flow: float = 0.0 + hp_include_elect_heat_in_power: bool = False + set_all_values_to_final_rosim_year: bool = False + max_flow: float = 500.0 + charge_temp: float = 80.0 + anisotropy: float = 5.0 + filter_fraction: float = 0.8 + salinity_surface: float = 0.0 + salinity_gradient: float = 46.67 + inj_prod_days: int = 120 + clogging_vel: float = 0.3 + membrane_filter_index: float = 1.0 + depth_mult_factor: float = 0.01 + thermal_radius_factor: float = 2.0 + + # DoubletCalc1DData + max_pump: float = 300.0 + min_pump: float = 1.0 + hy_gradient: float = 0.105 + optim_well_dist: bool = True + optim_dist_well_dist_min: float = 100.0 + optim_dist_well_dist_max: float = 3000.0 + optim_dist_lifetime: int = 50 + max_tvd_stepout_factor: float = 1.0 + optim_dist_cp_rock: float = 2700.0 + optim_dist_rho_rock: float = 1000.0 + optim_dist_cooling_fraction: float = 0.1 + default_well_distance: float = 1500.0 + pump_efficiency: float = 0.6 + pump_depth: float = 300.0 + segment_length: float = 50.0 + outer_diameter: float = 8.5 + inner_diameter: float = 8.5 + roughness: float = 1.38 + skin_injector: float = -1.0 + skin_producer: float = -1.0 + stim_add_skin_inj: float = -3.0 + stim_add_skin_prod: float = -3.0 + stimulation_capex: float = 0.5 + hp_calc_cop: bool = True + imposed_cop: float = 3.0 + hp_capex: float = 600.0 + hp_opex: float = 60.0 + hp_alternative_heating_price: float = 2.8 + viscosity_mode: ViscosityMode = "batzlewang" + + # Economical data + economic_lifetime: int = 15 + drilling_time: int = 2 + tax_rate: float = 0.25 + interest_loan: float = 0.05 + inflation: float = 0.02 + equity_return: float = 0.07 + debt_equity: float = 0.8 + ecn_eia: float = 0.0 + tolerance_utc_increase: float = 0.0 + load_hours: float = 6000.0 + opex_base: float = 0.0 + opex_per_power: float = 100.0 + elec_purchase_price: float = 8.0 + opex_per_energy: float = 0.19 + opex_per_capex: float = 0.0 + well_cost_scaling: float = 1.5 + well_cost_const: float = 0.25 + well_cost_z: float = 700.0 + well_cost_z2: float = 0.2 + capex_const: float = 3.0 + capex_variable: float = 300.0 + capex_contingency: float = 15.0 + well_curv_scaling: float = 1.1 + use_orc: bool = False + heat_exchanger_efficiency: float = 1.0 + heat_exchanger_parasitic: float = 0.0 + heat_exchanger_basetemp: float = 10.0 diff --git a/src/pythermogis/workflow/utc/water.py b/src/pythermogis/workflow/utc/water.py new file mode 100644 index 0000000000000000000000000000000000000000..82b5be2b0b7f56ad76e22a5d6f94cf2a2cf3242e --- /dev/null +++ b/src/pythermogis/workflow/utc/water.py @@ -0,0 +1,92 @@ +from numba import njit + + +@njit +def get_hydrostatic_pressure(depth: float) -> float: + return 9810.0 * depth + 101325.0 + + +@njit +def get_salinity(surface_salinity: float, gradient: float, depth: float) -> float: + return (surface_salinity + gradient * depth) / 1e6 + + +@njit +def density( + P: float, + T: float, + S: float, + pres_input_in_bar: bool = False, + salinity_input_in_ppm: bool = False, +) -> float: + if pres_input_in_bar: + P *= 1e5 # Units.BAR_SI + + if salinity_input_in_ppm: + S *= 1e-6 + + if P < 1.0: + P = 1.0 + if T < 0.0: + T = 0.0 + if S < 0.0: + S = 0.0 + + P_MPa = P * 1e-6 + + density_fresh = 1.0 + 1e-6 * ( + -80.0 * T + - 3.3 * T * T + + 0.00175 * T * T * T + + 489.0 * P_MPa + - 2.0 * T * P_MPa + + 0.016 * T * T * P_MPa + - 1.3e-5 * T * T * T * P_MPa + - 0.333 * P_MPa * P_MPa + - 0.002 * T * P_MPa * P_MPa + ) + + density_val = density_fresh + S * ( + 0.668 + + 0.44 * S + + 1e-6 + * ( + 300.0 * P_MPa + - 2400.0 * P_MPa * S + + T * (80.0 + 3.0 * T - 3300.0 * S - 13.0 * P_MPa + 47.0 * P_MPa * S) + ) + ) + + return density_val * 1000.0 + + +@njit +def heat_capacity( + P: float, T: float, S: float, salinity_input_in_ppm: bool = False +) -> float: + if salinity_input_in_ppm: + S *= 1e-6 + + if P < 1.0: + P = 1.0 + if T < 0.0: + T = 0.0 + if S < 0.0: + S = 0.0 + + T_K = T + 273.15 + S_g_per_kg = S * 1e3 + + heat_capacity_val = ( + (5.328 + -9.76e-2 * S_g_per_kg + 4.04e-4 * S_g_per_kg * S_g_per_kg) + + (-6.913e-3 + 7.351e-4 * S_g_per_kg - 3.15e-6 * S_g_per_kg * S_g_per_kg) * T_K + + (9.6e-6 - 1.927e-6 * S_g_per_kg + 8.23e-9 * S_g_per_kg * S_g_per_kg) + * T_K + * T_K + + (2.5e-9 + 1.666e-9 * S_g_per_kg - 7.125e-12 * S_g_per_kg * S_g_per_kg) + * T_K + * T_K + * T_K + ) + + return heat_capacity_val * 1e3 diff --git a/src/pythermogis/workflow/utc/well_distance.py b/src/pythermogis/workflow/utc/well_distance.py new file mode 100644 index 0000000000000000000000000000000000000000..4a585274297843fa5249167a6b9fc440d8e26528 --- /dev/null +++ b/src/pythermogis/workflow/utc/well_distance.py @@ -0,0 +1,138 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING + +import numpy as np +from scipy.optimize import brentq + +from pythermogis.workflow.utc.doublet_utils import calc_lifetime +from pythermogis.workflow.utc.flow import calculate_volumetric_flow + +if TYPE_CHECKING: + from pythermogis.workflow.utc.doublet import DoubletInput + from pythermogis.workflow.utc.utc_properties import UTCConfiguration + + +def optimize_well_distance_original( + props: UTCConfiguration, + input: DoubletInput, + drawdown_pressure: float, + injection_temp: float, +) -> float: + """ + The original well distance optimization algorithm. This is kept for testing purposes + but has been replaced with the scipy implementation of the brenth algorithm which + does the same thing more accurately. + + Parameters + ---------- + props + input + drawdown_pressure + injection_temp + + Returns + ------- + + """ + + dist_min = props.optim_dist_well_dist_min + dist_max = props.optim_dist_well_dist_max + well_distance = np.mean([dist_min, dist_max]) + + for _iter_count in range(1000): + if abs(dist_max - dist_min) <= 10.0: + return well_distance + well_distance = np.mean([dist_min, dist_max]) + + # --- Compute flow for this distance --- + results = calculate_volumetric_flow( + props=props, + input_data=input, + original_pressure=drawdown_pressure, + well_distance=well_distance, + injection_temp=injection_temp, + ) + + # --- Compute lifetime for this distance --- + lifetime = calc_lifetime( + well_distance=well_distance, + thickness=input.thickness * input.ntg, + delta_temp_fraction=props.optim_dist_cooling_fraction, + porosity=input.porosity, + flowrate=min(results.flowrate, props.max_flow), + depth=input.depth, + reservoir_temp=input.temperature, + salinity_surface=props.salinity_surface, + salinity_gradient=props.salinity_gradient, + cp_rock=props.optim_dist_cp_rock, + rho_rock=props.optim_dist_rho_rock, + ) + + # --- Bisection rule --- + if lifetime < props.optim_dist_lifetime: + dist_min = well_distance + else: + dist_max = well_distance + # If no convergence in 1000 iterations + else: + print( + f"WARNING: Well distance optimization failed to converge." + f" Final dist={well_distance}" + ) + + return well_distance + + +def f1( + well_distance: float, + props: UTCConfiguration, + input: DoubletInput, + drawdown_pressure: float, + injection_temp: float, +) -> float: + # --- Compute flow for this distance --- + results = calculate_volumetric_flow( + props=props, + input_data=input, + original_pressure=drawdown_pressure, + well_distance=well_distance, + injection_temp=injection_temp, + ) + + # --- Compute lifetime for this distance --- + return ( + calc_lifetime( + well_distance=well_distance, + thickness=input.thickness * input.ntg, + delta_temp_fraction=props.optim_dist_cooling_fraction, + porosity=input.porosity, + flowrate=min(results.flowrate, props.max_flow), + depth=input.depth, + reservoir_temp=input.temperature, + salinity_surface=props.salinity_surface, + salinity_gradient=props.salinity_gradient, + cp_rock=props.optim_dist_cp_rock, + rho_rock=props.optim_dist_rho_rock, + ) + - props.optim_dist_lifetime + ) + + +def optimize_well_distance( + props: UTCConfiguration, + input: DoubletInput, + drawdown_pressure: float, + injection_temp: float, +) -> float: + # find the well distance between the min and max which comes closest to the optimal + # doublet lifetime, as defined in props. + well_distance = brentq( + f1, + props.optim_dist_well_dist_min, + props.optim_dist_well_dist_max, + xtol=10, + maxiter=100, + args=(props, input, drawdown_pressure, injection_temp), + ) + return well_distance diff --git a/tests/test_dask_parralelization.py b/tests/test_dask_parralelization.py index 9573cf868cd31705b5b1368b2fe92351353e44a1..138d424cc81b2ea255204e1f79b7fd3294b0c6ab 100644 --- a/tests/test_dask_parralelization.py +++ b/tests/test_dask_parralelization.py @@ -2,15 +2,17 @@ from os import path from pathlib import Path import numpy as np +import pytest import xarray as xr from pythermogis import auto_chunk_dataset, assess_optimal_chunk_size +# @pytest.mark.skip("Useful for making a plot, but not needed for pipeline testing") def test_dask_parralelization(): output_data_path = Path(path.dirname(__file__), "resources") / "test_output" / "parallelization" output_data_path.mkdir(parents=True, exist_ok=True) - assess_optimal_chunk_size(n_simulations = 100, chunk_step_size=50, plot_outfile = output_data_path / "parallelization.png") + assess_optimal_chunk_size(n_simulations = 100, chunk_step_size=10, plot_outfile = output_data_path / "parallelization.png") assert output_data_path.exists() diff --git a/tests/test_doc_examples.py b/tests/test_doc_examples.py index faa28fd92f430322920f52bcfc6dd84e5a46defa..23a2b626576df0e0760ff1fd207c90d27d771a0d 100644 --- a/tests/test_doc_examples.py +++ b/tests/test_doc_examples.py @@ -280,7 +280,7 @@ def test_example8(): output_data_path.mkdir(parents=True, exist_ok=True) # generate simulation samples across desired reservoir properties - Nsamples = 100 + Nsamples = 1000 thickness_samples = np.random.uniform(low=150, high=300, size=Nsamples) porosity_samples = np.random.uniform(low=0.5, high=0.8, size=Nsamples) ntg_samples = np.random.uniform(low=0.25, high=0.5, size=Nsamples) @@ -298,7 +298,7 @@ def test_example8(): ) # For every sample, run a doublet simulation store the output values - simulations = calculate_doublet_performance(reservoir_properties, print_execution_duration=True) + simulations = calculate_doublet_performance(reservoir_properties, chunk_size=100, print_execution_duration=True) # post process the samples to calculate the percentiles of each variable percentiles = np.arange(1,99) diff --git a/tests/test_doublet_speed.py b/tests/test_doublet_speed.py index 36be0e513366e375dc6c38fae7c426740c70226f..c93073de07bd61d9c2028e0a9ba9916b4b07235c 100644 --- a/tests/test_doublet_speed.py +++ b/tests/test_doublet_speed.py @@ -27,15 +27,14 @@ class PyThermoGIS(TestCase): non_nan_data = data[~np.isnan(data)] print(f"Running Performance calculation for ROSL_ROSLU, Pvalues: {p_values}, The number of Non Nan Cells: {len(non_nan_data)}") - # Run calculation across all dimensions of input_grids, and all provided P_values: The Equivalent calculation with 10 cores takes 43 seconds in the Java code + # Run calculation across all dimensions of input_grids, and all provided P_values start = timeit.default_timer() - calculate_doublet_performance_stochastic(input_grids, chunk_size=100, rng_seed=123, p_values=p_values) + calculate_doublet_performance_stochastic(input_grids, chunk_size=10, rng_seed=123, p_values=p_values) time_elapsed = timeit.default_timer() - start print(f"Python calculation took {time_elapsed:.1f} seconds.") - def read_input_grids(self): - new_cellsize=5000 # in m + new_cellsize=10e3 # in m input_grids = resample_grid(read_grid(self.test_files_out_path / "ROSL_ROSLU__thick.zmap"), new_cellsize=new_cellsize).to_dataset(name="thickness_mean") input_grids["thickness_sd"] = resample_grid(read_grid(self.test_files_out_path / "ROSL_ROSLU__thick_sd.zmap"), new_cellsize=new_cellsize) input_grids["ntg"] = resample_grid(read_grid(self.test_files_out_path / "ROSL_ROSLU__ntg.zmap"), new_cellsize=new_cellsize) diff --git a/tests/utc/test_calculate_volumetric_flow.py b/tests/utc/test_calculate_volumetric_flow.py new file mode 100644 index 0000000000000000000000000000000000000000..f3059727e750f4629d42d09c788abb5865c8dd34 --- /dev/null +++ b/tests/utc/test_calculate_volumetric_flow.py @@ -0,0 +1,31 @@ +import numpy as np + +from pythermogis.workflow.utc.doublet import DoubletInput +from pythermogis.workflow.utc.flow import calculate_volumetric_flow +from pythermogis.workflow.utc.utc_properties import UTCConfiguration + + +def test_calculate_volumetric_flow(): + props = UTCConfiguration() + + input_data = DoubletInput( + unknown_input_value=-999.0, + thickness=100.0, + transmissivity=17500.0, + transmissivity_with_ntg=0.0, + ntg=1.0, + depth=2000, + porosity=0.0, + temperature=50.0, + ) + + flow_results = calculate_volumetric_flow(props, input_data, 600000, 1550, 40) + rtol = 0.01 + assert np.isclose(flow_results.hp_cop, 0.0, rtol=rtol) + assert np.isclose(flow_results.hp_added_power, 0.0, rtol=rtol) + assert np.isclose(flow_results.heat_power_per_doublet, 0.027926914290870505, rtol=rtol) + assert np.isclose(flow_results.cop, 5.396068601571214, rtol=rtol) + assert np.isclose(flow_results.flowrate, 18.631507399806665, rtol=rtol) + assert np.isclose(flow_results.pump_power_required, 5.175418689587975, rtol=rtol) + assert np.isclose(flow_results.production_temp, 41.36211427733413, rtol=rtol) + assert np.isclose(flow_results.heat_power_produced[0], 0.027926914290870505, rtol=rtol) \ No newline at end of file diff --git a/tests/utc/test_doublet.py b/tests/utc/test_doublet.py new file mode 100644 index 0000000000000000000000000000000000000000..e42d93acbc41a18effb6747f9451640e18258733 --- /dev/null +++ b/tests/utc/test_doublet.py @@ -0,0 +1,142 @@ +import timeit + +import numpy as np + + +from pythermogis.workflow.utc.doublet import DoubletInput, calculate_doublet_performance +from pythermogis.workflow.utc.utc_properties import UTCConfiguration + + +def test_calculate_doublet_performance_precise(): + # Arrange: instantiate default UTCConfiguration + props = UTCConfiguration( + viscosity_mode="kestin", + dh_return_temp=40, + ) + + # Create a minimal valid DoubletInput + input_data = DoubletInput( + unknown_input_value=-999.0, + thickness=100.0, + transmissivity=17500.0, + transmissivity_with_ntg=0.0, + ntg=1.0, + depth=2000, + porosity=0.0, + temperature=76.0, + ) + + # Act + # perform one simulation to compile all the numba functions before timing + calculate_doublet_performance(props, input_data) + + start = timeit.default_timer() + result = calculate_doublet_performance(props, input_data, verbose=True) + time_elapsed = timeit.default_timer() - start + + print(f"1 simulation took: {time_elapsed:.1f} seconds\n" + f"{1/time_elapsed:.1f} simulations per second") + + + # Assert + rtol = 0.005 # accurate to 0.5% + assert np.isclose(result.flow, 227.2757568359375, rtol=rtol) + assert np.isclose(result.pres, 60, rtol=rtol) + assert np.isclose(result.utc, 6.616096470753937, rtol=rtol) + assert np.isclose(result.welld, 1159.17968, rtol=rtol) + assert np.isclose(result.power, 8.636903762817383, rtol=rtol) + assert np.isclose(result.cop, 13.627557754516602, rtol=rtol) + assert np.isclose(result.var_opex, -7.510325908660889, rtol=rtol) + assert np.isclose(result.fixed_opex, -11.227973937988281, rtol=rtol) + +def test_calculate_doublet_performance_approximate(): + # Arrange: instantiate default UTCConfiguration + props = UTCConfiguration( + viscosity_mode="kestin", + dh_return_temp=40, + segment_length=1.0, + ) + + # Create a minimal valid DoubletInput + input_data = DoubletInput( + unknown_input_value=-999.0, + thickness=100.0, + transmissivity=17500.0, + transmissivity_with_ntg=0.0, + ntg=1.0, + depth=2000, + porosity=0.0, + temperature=76.0, + ) + + # Act + #perform one simulation to compile all the numba functions before timing + result = calculate_doublet_performance(props, input_data) + + start = timeit.default_timer() + n_sims = 1000 + for _ in range(n_sims): + result = calculate_doublet_performance(props, input_data) + time_elapsed = timeit.default_timer() - start + + print(f"{n_sims} simulations took: {time_elapsed:.1f} seconds\n" + f"{n_sims/time_elapsed:.1f} simulations per second") + + # Assert + rtol = 0.01 # accurate to 1% + assert np.isclose(result.flow, 227.2757568359375, rtol=rtol) + assert np.isclose(result.pres, 60, rtol=rtol) + assert np.isclose(result.utc, 6.616096470753937, rtol=rtol) + assert np.isclose(result.welld, 1159.17968, rtol=rtol) + assert np.isclose(result.power, 8.636903762817383, rtol=rtol) + assert np.isclose(result.cop, 13.627557754516602, rtol=rtol) + assert np.isclose(result.var_opex, -7.510325908660889, rtol=rtol) + assert np.isclose(result.fixed_opex, -11.227973937988281, rtol=rtol) + +def test_calculate_doublet_performance(): + # 85ish it/s + props = UTCConfiguration( + viscosity_mode="kestin", + dh_return_temp=40, + segment_length=1.0, + ) + + warmup_input = DoubletInput( + unknown_input_value=-999.0, + thickness=100.0, + transmissivity=17500.0, + transmissivity_with_ntg=0.0, + ntg=1.0, + depth=2000.0, + porosity=0.0, + temperature=76.0, + ) + _ = calculate_doublet_performance(props, warmup_input) + + start = timeit.default_timer() + n_sims = 1000 + + for i in range(n_sims): + transmissivity = 17500.0 + (i % 10) * 5 + temp = 76.0 + (i % 5) * 0.1 + depth = 2000.0 + (i % 4) * 2 + + input_data = DoubletInput( + unknown_input_value=-999.0, + thickness=100.0, + transmissivity=transmissivity, + transmissivity_with_ntg=0.0, + ntg=1.0, + depth=depth, + porosity=0.0, + temperature=temp, + ) + + _ = calculate_doublet_performance(props, input_data) + + time_elapsed = timeit.default_timer() - start + + print( + f"{n_sims} simulations took: {time_elapsed:.1f} seconds\n" + f"{n_sims/time_elapsed:.1f} simulations per second" + ) diff --git a/tests/utc/test_doublet_utils.py b/tests/utc/test_doublet_utils.py new file mode 100644 index 0000000000000000000000000000000000000000..385bac89d561961ee7d50cff0c04c1612b75ea2d --- /dev/null +++ b/tests/utc/test_doublet_utils.py @@ -0,0 +1,65 @@ +import pytest +from pythermogis.workflow.utc.doublet_utils import calculate_injection_temp_with_heat_pump + + +EPSILON = 0.00000001 + + +def test_calculate_injection_temp_with_heat_pump_normal_case(): + t_prod = 50 + t_goal = 70 + dh_return_temp = 40 + reservoir_temp = 60 + max_cooling_temp_range = 15 + hp_minimum_injection_temperature = 30 + + result = calculate_injection_temp_with_heat_pump( + t_prod, + t_goal, + dh_return_temp, + reservoir_temp, + max_cooling_temp_range, + hp_minimum_injection_temperature + ) + + assert result == pytest.approx(45.0, abs=1e-6) + + +def test_calculate_injection_temp_with_heat_pump_enforces_minimum_injection_temp(): + t_prod = 50 + t_goal = 70 + dh_return_temp = 40 + reservoir_temp = 35 + max_cooling_temp_range = 10 + hp_minimum_injection_temperature = 40 + + result = calculate_injection_temp_with_heat_pump( + t_prod, + t_goal, + dh_return_temp, + reservoir_temp, + max_cooling_temp_range, + hp_minimum_injection_temperature + ) + + assert result == pytest.approx(40.0, abs=EPSILON) + + +def test_calculate_injection_temp_with_heat_pump_waste_temp_lower_than_prod_handled(): + t_prod = 60 + t_goal = 80 + dh_return_temp = 50 + reservoir_temp = 70 + max_cooling_temp_range = 20 + hp_minimum_injection_temperature = 35 + + result = calculate_injection_temp_with_heat_pump( + t_prod, + t_goal, + dh_return_temp, + reservoir_temp, + max_cooling_temp_range, + hp_minimum_injection_temperature + ) + + assert result == pytest.approx(50.0, abs=1e-6) diff --git a/tests/utc/test_doubletcalc.py b/tests/utc/test_doubletcalc.py new file mode 100644 index 0000000000000000000000000000000000000000..e49d5d10925f030c41ca77b3e45b39ee32d0e8c7 --- /dev/null +++ b/tests/utc/test_doubletcalc.py @@ -0,0 +1,64 @@ +import timeit + +import numpy as np +import pytest + +from pythermogis.workflow.utc.doublet import DoubletInput +from pythermogis.workflow.utc.doubletcalc import doubletcalc +from pythermogis.workflow.utc.utc_properties import UTCConfiguration + + +def test_doubletcalc(): + props = UTCConfiguration() + + input_data = DoubletInput( + unknown_input_value=-999.0, + thickness=100.0, + transmissivity=17500.0, + transmissivity_with_ntg=17500.0, + ntg=1.0, + depth=2000, + porosity=0.0, + temperature=50.0, + ) + + results = doubletcalc(props, input_data, 600000, 1550, 40) + + assert np.isclose(results.geothermal_powers, 0.027926914290870505, rtol=0.01) + assert np.isclose(results.cop, 5.396068601571214, rtol=0.01) + assert np.isclose(results.flowrate, 18.631507399806665, rtol=0.01) + assert np.isclose(results.pump_power_required, 5.175418689587975, rtol=0.01) + assert np.isclose(results.production_temp, 41.36211427733413, rtol=0.01) + assert np.isclose(results.heat_power_produced[0], 0.027926914290870505, rtol=0.01) + +def test_doubletcalc_performance(): + props = UTCConfiguration(segment_length=1) + + input_data = DoubletInput( + unknown_input_value=-999.0, + thickness=100.0, + transmissivity=17500.0, + transmissivity_with_ntg=17500.0, + ntg=1.0, + depth=2000, + porosity=0.0, + temperature=50.0, + ) + + nsims = 10000 + # one iteration to warm up the JIT compiled functions + doubletcalc(props, input_data, 600000, 1550, 40) + + start = timeit.default_timer() + for _ in range(nsims): + results = doubletcalc(props, input_data, 600000, 1550, 40) + time_elapsed = timeit.default_timer() - start + print(f"{(nsims/time_elapsed):.1f} sims/s") + + rtol=0.25 + assert np.isclose(results.geothermal_powers, 0.027926914290870505, rtol=rtol) + assert np.isclose(results.cop, 5.396068601571214, rtol=rtol) + assert np.isclose(results.flowrate, 18.631507399806665, rtol=rtol) + assert np.isclose(results.pump_power_required, 5.175418689587975, rtol=rtol) + assert np.isclose(results.production_temp, 41.36211427733413, rtol=rtol) + assert np.isclose(results.heat_power_produced[0], 0.027926914290870505, rtol=rtol) \ No newline at end of file diff --git a/tests/utc/test_properties.py b/tests/utc/test_properties.py new file mode 100644 index 0000000000000000000000000000000000000000..ab98eb17c08247970d6bd8f91131fd2721919c52 --- /dev/null +++ b/tests/utc/test_properties.py @@ -0,0 +1,29 @@ +from pythermogis.workflow.utc.utc_properties import UTCConfiguration + + +def test_utc_configuration_instantiation(): + """Test that UTCConfiguration can be instantiated with defaults.""" + cfg = UTCConfiguration() + + assert isinstance(cfg, UTCConfiguration) + assert isinstance(cfg.property_grid_infos, list) + assert isinstance(cfg.copy_aquifer_files_info, list) + assert cfg.n_threads == 1 + assert cfg.is_ates is False + assert cfg.scenario == "basecase" + assert cfg.viscosity_mode == "batzlewang" + assert cfg.surface_temperature == 10.0 + assert cfg.temp_gradient == 31.0 + assert 30.0 in cfg.p_values + assert cfg.p_values == [10.0, 30.0, 50.0, 70.0, 90.0] + assert cfg.property_grid_infos[0].name == "Permeability" + assert cfg.property_grid_infos[0].optional is False + assert cfg.property_grid_infos[0].postfix.endswith(".zmap") + assert len(cfg.aquifers) > 10 + assert "NMRFT" in cfg.aquifers + assert cfg.petrel_output_dir is None + assert cfg.tax_rate == 0.25 + assert cfg.interest_loan == 0.05 + assert cfg.economic_lifetime == 15 + + diff --git a/tests/utc/test_well_distance_optimizer.py b/tests/utc/test_well_distance_optimizer.py new file mode 100644 index 0000000000000000000000000000000000000000..00401228f3d88d7f9051a2bcfe6dba0af8f2f36c --- /dev/null +++ b/tests/utc/test_well_distance_optimizer.py @@ -0,0 +1,46 @@ +import numpy as np + +from pythermogis.workflow.utc.doublet import DoubletInput +from pythermogis.workflow.utc.utc_properties import UTCConfiguration +from pythermogis.workflow.utc.well_distance import ( + optimize_well_distance, + optimize_well_distance_original, +) + + +def test_well_distance_optimizer(): + props = UTCConfiguration() + + input_data = DoubletInput( + unknown_input_value=-999.0, + thickness=100.0, + transmissivity=17500.0, + transmissivity_with_ntg=0.0, + ntg=1.0, + depth=2000, + porosity=0.0, + temperature=50.0, + ) + + well_distance = optimize_well_distance_original(props, input_data, 6000000, 40) + + assert np.isclose(well_distance, 1011.9140625, rtol=0.015) + + +def test_well_distance_optimizer2(): + props = UTCConfiguration() + + input_data = DoubletInput( + unknown_input_value=-999.0, + thickness=100.0, + transmissivity=17500.0, + transmissivity_with_ntg=0.0, + ntg=1.0, + depth=2000, + porosity=0.0, + temperature=50.0, + ) + + well_distance = optimize_well_distance(props, input_data, 6000000, 40) + + assert np.isclose(well_distance, 1011.9140625, rtol=0.015) \ No newline at end of file