diff --git a/pixi.lock b/pixi.lock index 468a60a20e93323e6d03b4b4eeffcccb494b05d9..fcf6deac9135016b3a06580f1e01cc9b3f62d41d 100644 --- a/pixi.lock +++ b/pixi.lock @@ -206,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 @@ -252,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 @@ -263,7 +265,9 @@ environments: - pypi: https://files.pythonhosted.org/packages/2a/2f/63d2cacc0e525f8e3398bcf32bd3620385f22cd1600834ec49d7f3597a7b/rioxarray-0.19.0-py3-none-any.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 @@ -436,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 @@ -486,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 @@ -497,7 +503,9 @@ environments: - pypi: https://files.pythonhosted.org/packages/2a/2f/63d2cacc0e525f8e3398bcf32bd3620385f22cd1600834ec49d7f3597a7b/rioxarray-0.19.0-py3-none-any.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 @@ -4471,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 @@ -4752,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.2 + sha256: 89fbb20be0f5da6d0e8e543b9160e4ac3aef8369557ca71ae4a96086f521145e + 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 @@ -4919,7 +4945,7 @@ packages: - pypi: ./ name: pythermogis version: 1.2.0 - sha256: 44271718a4a84514ee6d1f6a413b57b060a4868ea6dd8fadf0e1ed09c39dfae1 + sha256: 35de67af8136fb6442b90f5a2956e6ac57f8a6f52f1fbbfad15768e83e08b31f requires_dist: - jpype1>=1.5.2,<2 - xarray==2024.9.0.* @@ -5305,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 @@ -5616,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 281dff81bfc1a613ce9e68882bb7bfac7ce7397b..b0c22a3131b0acb86fdaf281056289941822f0c8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -70,6 +70,7 @@ 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/workflow/utc/cooling_temp.py b/src/pythermogis/workflow/utc/cooling_temp.py index bd62df841e0f1064524435d0ac2239597b7c69bb..7d279ff8dd67c0505c82ca68f7d1894baa3a60d0 100644 --- a/src/pythermogis/workflow/utc/cooling_temp.py +++ b/src/pythermogis/workflow/utc/cooling_temp.py @@ -1,15 +1,10 @@ -from pythermogis.workflow.utc.flow import calculate_volumetric_flow from pythermogis.workflow.utc.doublet_utils import calc_lifetime +from pythermogis.workflow.utc.flow import calculate_volumetric_flow def calculate_cooling_temperature( - props, - input, - drawdown_pressure: float, - well_distance: float, - injection_temp: float + props, input, drawdown_pressure: float, well_distance: float, injection_temp: float ) -> float: - results = calculate_volumetric_flow( props=props, input_data=input, @@ -51,7 +46,8 @@ def calculate_cooling_temperature( if (not props.is_ates) and cooling_frac > 0.3: raise RuntimeError( - f"Large cooling factor ({cooling_frac}) could result in less reliable calculation" + 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 index 71debd9c7b4969ff5d2e3fa9da4f07be49c7494d..7dcc9676734443c2138b600b0d59eb5adbd0b111 100644 --- a/src/pythermogis/workflow/utc/doublet.py +++ b/src/pythermogis/workflow/utc/doublet.py @@ -16,6 +16,7 @@ from pythermogis.workflow.utc.well_distance import optimize_well_distance EUR_PER_CT_PER_KWH = 0.36 NPV_SCALE = 1e-6 + @dataclass class DoubletInput: unknown_input_value: float @@ -32,6 +33,7 @@ class DoubletInput: DARCY_SI = 1.0e-12 / 1.01325 return self.transmissivity / self.thickness * 1e-3 * DARCY_SI + @dataclass class DoubletOutput: power: float @@ -53,7 +55,10 @@ class DoubletOutput: production_temp: float injection_temp: float -def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput, verbose: bool = False) -> DoubletOutput | None: + +def calculate_doublet_performance( + props: UTCConfiguration, input: DoubletInput, verbose: bool = False +) -> DoubletOutput | None: timer = timeit.default_timer() # determine initial well distance well_distance = ( @@ -73,8 +78,9 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput, props.hp_minimum_injection_temperature, ) else: - injection_temperature = max(input.temperature - props.max_cooling_temp_range, - props.dh_return_temp) + 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 @@ -91,9 +97,8 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput, # cooling temperature and well distance optimization if props.optim_well_dist: - end_temperature_p = ( - props.optim_dist_cooling_fraction * - (input.temperature - injection_temperature) + end_temperature_p = props.optim_dist_cooling_fraction * ( + input.temperature - injection_temperature ) else: end_temperature_p = calculate_cooling_temperature( @@ -117,7 +122,10 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput, # stimulation capex stimulation_capex = ( 0.0 - if (not props.use_stimulation or input.transmissivity_with_ntg > props.stim_kh_max) + if ( + not props.use_stimulation + or input.transmissivity_with_ntg > props.stim_kh_max + ) else props.stimulation_capex ) @@ -167,7 +175,7 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput, breakthrough_years = calc_lifetime( well_distance, - input.thickness*input.ntg, + input.thickness * input.ntg, 0.001, input.porosity, flowrate, @@ -176,19 +184,21 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput, props.salinity_surface, props.salinity_gradient, props.optim_dist_cp_rock, - props.optim_dist_rho_rock + props.optim_dist_rho_rock, ) utc_cutoff = ( - props.utc_cutoff_deep if input.depth > props.utc_deep_depth else props.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 - * discounted_heat_produced_p - * (1 - props.tax_rate) + NPV_SCALE + * (utc_cutoff - utc_eur_ct_per_kwh) + * 3.6 + * discounted_heat_produced_p + * (1 - props.tax_rate) ) opex_first_prod_year = total_opex_ts[props.drilling_time] print_time(timer, "economics: ", verbose=verbose) @@ -211,6 +221,5 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput, breakthrough=breakthrough_years, cooling=end_temperature_p, production_temp=production_temp, - injection_temp=injection_temperature + injection_temp=injection_temperature, ) - diff --git a/src/pythermogis/workflow/utc/doublet_utils.py b/src/pythermogis/workflow/utc/doublet_utils.py index 12ac94c9449a8f04cb81e65d1745e06d6cd8634b..bf3200da7f9ebdccb0447d21d9a3bbf508c476fa 100644 --- a/src/pythermogis/workflow/utc/doublet_utils.py +++ b/src/pythermogis/workflow/utc/doublet_utils.py @@ -1,6 +1,14 @@ import math + from numba import njit -from pythermogis.workflow.utc.water import density, heat_capacity, get_salinity, get_hydrostatic_pressure + +from pythermogis.workflow.utc.water import ( + density, + get_hydrostatic_pressure, + get_salinity, + heat_capacity, +) + @njit def calculate_injection_temp_with_heat_pump( @@ -33,13 +41,15 @@ def calculate_injection_temp_with_heat_pump( return max( needed_injection_temp, - max(reservoir_temp - max_cooling_temp_range, hp_minimum_injection_temperature) + 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: @@ -59,6 +69,7 @@ def get_cop_carnot(eta: float, Tout: float, Tin: float) -> float: return eta * (Tcond + TKELVIN) / (Tcond - Tevap) + @njit def calc_lifetime( well_distance: float, @@ -73,7 +84,6 @@ def calc_lifetime( cp_rock: float, rho_rock: float, ) -> float: - # --- Water properties at depth --- salinity = get_salinity(salinity_surface, salinity_gradient, depth) pressure = get_hydrostatic_pressure(depth) @@ -102,9 +112,10 @@ def calc_lifetime( Aseg = Aseg1 - Aseg2 Vseg = Aseg * thickness - Eseg = 1e-9 * Vseg * ( - (1 - porosity) * cp_rock * rho_rock + - porosity * cp_water * rho_water + Eseg = ( + 1e-9 + * Vseg + * ((1 - porosity) * cp_rock * rho_rock + porosity * cp_water * rho_water) ) flowseg = flowrate * 365 * 24 * dangle / math.pi @@ -120,7 +131,6 @@ def get_along_hole_length( curve_scaling_factor: float, max_true_vertical_depth_stepout_factor: float, ) -> float: - if curve_scaling_factor == 0: return true_vertical_depth @@ -134,6 +144,8 @@ def get_along_hole_length( stepout -= horizontal_distance - oblique_distance = math.sqrt(stepout**2 + true_vertical_depth**2) * curve_scaling_factor + oblique_distance = ( + math.sqrt(stepout**2 + true_vertical_depth**2) * curve_scaling_factor + ) - return oblique_distance + horizontal_distance \ No newline at end of file + return oblique_distance + horizontal_distance diff --git a/src/pythermogis/workflow/utc/doubletcalc.py b/src/pythermogis/workflow/utc/doubletcalc.py index e57ad1b0b3c178ab63d5da05a874074006206cd5..6e1194cc927e0c43832a5a9a5697cc37c5ea0816 100644 --- a/src/pythermogis/workflow/utc/doubletcalc.py +++ b/src/pythermogis/workflow/utc/doubletcalc.py @@ -10,6 +10,7 @@ from pythermogis.workflow.utc.water import get_salinity INCH_SI = 0.0254 + @dataclass class Doublet1DResults: geothermal_powers: float @@ -19,21 +20,26 @@ class Doublet1DResults: production_temp: float heat_power_produced: list[float] + def doubletcalc( - props: UTCConfiguration, - input, - drawdown_pressure: float, - well_distance: float, - injection_temp: float, + props: UTCConfiguration, + input, + 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), + 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), + geothermal_gradient=get_geothermal_gradient( + input.depth, input.thickness, input.temperature, props.surface_temperature + ), rock_heat_capacity=855, rock_density=2720, thermal_conductivity=3, @@ -73,12 +79,19 @@ def doubletcalc( 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 - )], + 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, @@ -117,28 +130,39 @@ def doubletcalc( flowrate=doublet.flowrate, pump_power_required=doublet.pump_power, production_temp=doublet.production_temp, - heat_power_produced=[doublet.geothermal_power,] * props.econ_lifetime_years, + heat_power_produced=[ + doublet.geothermal_power, + ] + * props.econ_lifetime_years, ) + def get_total_skin_injection(props, input): - if (not props.use_stimulation) or (input.transmissivity_with_ntg > props.stim_kh_max): + 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, input): - if (not props.use_stimulation) or (input.transmissivity_with_ntg > props.stim_kh_max): + 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, depth: float) -> float: return min(props.pump_depth, depth / 2) + @njit def get_along_hole_length( true_vertical_depth: float, @@ -161,4 +185,4 @@ def get_along_hole_length( oblique_distance = ( math.sqrt(stepout**2 + true_vertical_depth**2) * curve_scaling_factor ) - return oblique_distance + horizontal_distance \ No newline at end of file + return oblique_distance + horizontal_distance diff --git a/src/pythermogis/workflow/utc/economics.py b/src/pythermogis/workflow/utc/economics.py index 958f4af129ed8c0d505be6155b6ce5f55a6a7f4f..82e951486bd66e40d2ffa3a2fb4091fbe0cb74c7 100644 --- a/src/pythermogis/workflow/utc/economics.py +++ b/src/pythermogis/workflow/utc/economics.py @@ -1,6 +1,7 @@ +from dataclasses import dataclass + import numpy as np -from dataclasses import dataclass from pythermogis.workflow.utc.doublet_utils import get_along_hole_length SECONDS_IN_YEAR = 365 * 24 * 3600 @@ -8,6 +9,7 @@ MJ_TO_GJ = 1e-3 KWH_TO_MJ = 0.36 * 1e9 HOURS_IN_YEAR = 8760 + @dataclass class CapexCalculatorResults: sum_capex: float @@ -17,16 +19,19 @@ class CapexCalculatorResults: total_opex_ts: list[float] heat_power_per_year: list[float] + @dataclass class UTCCalculatorResults: discounted_heat_produced: float utc: float + @dataclass class EconomicsResults: capex: CapexCalculatorResults utc: UTCCalculatorResults + def calculate_economics( props, input, @@ -40,8 +45,6 @@ def calculate_economics( hp_added_power: float, pump_power_required: float, ) -> EconomicsResults: - - ah_length_array = [ get_along_hole_length( true_vertical_depth=input.depth, @@ -76,6 +79,7 @@ def calculate_economics( utc=utc_results, ) + def calculate_capex( props, heat_power_produced: list[float], @@ -88,7 +92,6 @@ def calculate_capex( ah_length_array: list[float], pump_power_required: float, ) -> CapexCalculatorResults: - n_years = props.economic_lifetime heat_power_per_year = np.zeros(n_years) @@ -103,8 +106,12 @@ def calculate_capex( capex_well = calculate_capex_for_wells(props, ah_length_array, stimulation_capex) allocate_capex_for_wells(props, total_capex_ts, capex_well) - hp_added_power = get_max_hp_added_power(initial_hp_added_power, hp_added_power_years) - hp_after_hp = calculate_hp_added_power_after_heat_pump(props, hp_added_power, hp_cop) + hp_added_power = get_max_hp_added_power( + initial_hp_added_power, hp_added_power_years + ) + hp_after_hp = calculate_hp_added_power_after_heat_pump( + props, hp_added_power, hp_cop + ) installation_mw = calculate_installation_mw(shifted_heat_power, hp_added_power) total_capex_1year = calculate_total_capex_1year( @@ -150,6 +157,7 @@ def shift_time_series(series: list[float], shift: int): for i in range(shift): series[i] = 0 + def calculate_capex_for_wells(props, ah_length_array, stimulation_capex): inj_depth = ah_length_array[0] prod_depth = ah_length_array[0] @@ -161,12 +169,13 @@ def calculate_capex_for_wells(props, ah_length_array, stimulation_capex): ) return capex * props.well_cost_scaling + def calculate_well_cost(props, depth): return ( - props.well_cost_z2 * depth**2 - + props.well_cost_z * depth + props.well_cost_z2 * depth**2 + props.well_cost_z * depth ) * 1e-6 + props.well_cost_const + def allocate_capex_for_wells(props, total_capex_ts, capex_well): drilling_years = props.drilling_time yearly_cost = capex_well / max(drilling_years, 1) @@ -174,19 +183,23 @@ def allocate_capex_for_wells(props, total_capex_ts, capex_well): for y in range(drilling_years): total_capex_ts[y] += yearly_cost + def get_max_hp_added_power(initial, hp_power_years): if hp_power_years is not None: return max(hp_power_years) return initial + def calculate_hp_added_power_after_heat_pump(props, hp_added_power, hp_cop): alt_price = props.hp_alternative_heating_price return hp_added_power if alt_price < 0 else hp_added_power * (1 + 1 / (hp_cop - 1)) + def calculate_installation_mw(heat_power_produced, hp_added_power): last_val = heat_power_produced[-1] return max(last_val - hp_added_power, 0) + def calculate_total_capex_1year(props, total_capex_ts, hp_after_hp, installation_mw): last_year = max(props.drilling_time - 1, 0) @@ -195,10 +208,11 @@ def calculate_total_capex_1year(props, total_capex_ts, hp_after_hp, installation + (props.hp_capex * hp_after_hp + props.capex_variable * installation_mw) * 1e-3 ) - total_capex_ts[last_year] *= (1 + props.capex_contingency / 100) + total_capex_ts[last_year] *= 1 + props.capex_contingency / 100 return total_capex_ts[last_year] + def process_annual_data( props, heat_power_produced, @@ -215,11 +229,9 @@ def process_annual_data( season_factor_years, pump_power_required, ): - sum_capex = 0.0 for year in range(props.economic_lifetime): - # Season factor season_factor = ( season_factor_years[year] @@ -232,10 +244,7 @@ def process_annual_data( # Annual heat (GJ/yr) heat_power_per_year[year] = ( - heat_power_produced[year] - * season_factor - * SECONDS_IN_YEAR - * MJ_TO_GJ + heat_power_produced[year] * season_factor * SECONDS_IN_YEAR * MJ_TO_GJ ) sum_capex += total_capex_ts[year] @@ -270,6 +279,7 @@ def process_annual_data( return sum_capex + def calculate_variable_opex( props, season_factor, @@ -283,32 +293,46 @@ def calculate_variable_opex( ): pump_cost = ( -(pump_power * 1e3) - * SECONDS_IN_YEAR * season_factor * elec_price - / KWH_TO_MJ * 1e-6 + * 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 + * SECONDS_IN_YEAR + * season_factor + * elec_price + / KWH_TO_MJ + * 1e-6 ) opex_per_energy = ( -(1 / 0.36) - * props.opex_per_energy * 1e-2 - * heat_power_gj_per_year * 1e4 / 36 * 1e-6 + * props.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 + * SECONDS_IN_YEAR + * season_factor + * elec_price + / KWH_TO_MJ + * 1e-6 ) return pump_cost + parasitic_cost + heat_pump_cost + opex_per_energy + def calculate_fixed_opex(props, hp_after_hp, installation_mw, total_capex): return ( -props.opex_base / 1e6 @@ -316,6 +340,7 @@ def calculate_fixed_opex(props, hp_after_hp, installation_mw, total_capex): - total_capex * props.opex_per_capex / 100 ) + def get_heat_exchanger_season_factor(props): return props.load_hours / HOURS_IN_YEAR @@ -326,7 +351,6 @@ def calculate_utc( total_opex_ts, total_capex, ) -> 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) @@ -336,7 +360,6 @@ def calculate_utc( 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( @@ -369,11 +392,13 @@ def calculate_utc( utc=utc, ) + def calculate_net_cash_worth_eia(props, total_capex): 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, present_value): if props.interest_loan == 0: return 0.0 @@ -381,16 +406,21 @@ def calculate_payment_value(props, present_value): factor = (1 + props.interest_loan) ** props.economic_lifetime return (props.interest_loan / (factor - 1)) * (-(present_value * factor)) + def calculate_depreciation_cost(props, total_capex): return -(total_capex / props.economic_lifetime) + def calculate_future_value(present_value, payment, interest, year): 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, total_capex, discounted_income, discounted_heat_produced): + +def calculate_unit_technical_cost( + props, total_capex, discounted_income, discounted_heat_produced +): if discounted_heat_produced <= 0: return 0.0 return ( @@ -398,8 +428,8 @@ def calculate_unit_technical_cost(props, total_capex, discounted_income, discoun / discounted_heat_produced ) * 1e6 + def calculate_project_interest_rate(props): - return ( - (1 - props.tax_rate) * props.interest_loan * props.debt_equity + - (1 - props.debt_equity) * props.equity_return - ) + 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 index 97411a17ea82fc6a67b2188f8ab1233cc6e382f6..fafc33eec677767c6a5635e1b9f5e913f89910d9 100644 --- a/src/pythermogis/workflow/utc/flow.py +++ b/src/pythermogis/workflow/utc/flow.py @@ -1,9 +1,10 @@ from dataclasses import dataclass + import numpy as np -from pythermogis.workflow.utc.heatpump import calculate_heat_pump_performance 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 @dataclass @@ -28,11 +29,7 @@ def calculate_volumetric_flow( ): for step in [0, 1e5, -1e5, 2e5, -2e5, 3e5, -3e5]: d1d_results = doubletcalc( - props, - input_data, - original_pressure + step, - well_distance, - injection_temp + props, input_data, original_pressure + step, well_distance, injection_temp ) if d1d_results is not None: break @@ -94,4 +91,4 @@ def calculate_volumetric_flow( d1d_results.pump_power_required, d1d_results.production_temp, heat_power_produced, - ) \ No newline at end of file + ) diff --git a/src/pythermogis/workflow/utc/heatpump.py b/src/pythermogis/workflow/utc/heatpump.py index 1043160a8663e271d05d4d559e36e5028f65e304..1f31362cc2a30ce0abde2685cb0eb7bcbfc06b4b 100644 --- a/src/pythermogis/workflow/utc/heatpump.py +++ b/src/pythermogis/workflow/utc/heatpump.py @@ -1,7 +1,13 @@ from dataclasses import dataclass from pythermogis.workflow.utc.doublet_utils import get_cop_carnot -from pythermogis.workflow.utc.water import get_hydrostatic_pressure, get_salinity, density, heat_capacity +from pythermogis.workflow.utc.water import ( + density, + get_hydrostatic_pressure, + get_salinity, + heat_capacity, +) + @dataclass class HeatPumpPerformanceResults: @@ -19,7 +25,6 @@ def calculate_heat_pump_performance( ) -> HeatPumpPerformanceResults: ETA_CARNOT = 0.6 - heat_pump_start_temp = calculate_heat_pump_start_temp( props, production_temp, injection_temp ) @@ -32,9 +37,7 @@ def calculate_heat_pump_performance( input_data.depth, ) - rho_water = density( - hydrostatic_pressure, props.surface_temperature, salinity - ) + rho_water = density(hydrostatic_pressure, props.surface_temperature, salinity) cp_water = heat_capacity(hydrostatic_pressure, props.surface_temperature, salinity) @@ -55,16 +58,18 @@ def calculate_heat_pump_performance( hp_cop = get_cop_carnot(ETA_CARNOT, Tout, Tin) # calculate hp_elec_consumption - hp_elec_consumption = hp_added_power / (hp_cop -1) + 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 + hp_elec_consumption=hp_elec_consumption, ) -def calculate_heat_pump_start_temp(props, production_temp: float, injection_temp: float) -> float: +def calculate_heat_pump_start_temp( + props, production_temp: float, injection_temp: float +) -> float: """ Python version of calculateHeatPumpStartTemp(...) """ @@ -75,4 +80,4 @@ def calculate_heat_pump_start_temp(props, production_temp: float, injection_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 \ No newline at end of file + return production_temp - delta_temp_for_hex diff --git a/src/pythermogis/workflow/utc/pressure.py b/src/pythermogis/workflow/utc/pressure.py index 40179fc41694b14df50264fcc128b084941eb74c..3e5ee06fde6c7f989d62c863f5d68978ae8d4cb0 100644 --- a/src/pythermogis/workflow/utc/pressure.py +++ b/src/pythermogis/workflow/utc/pressure.py @@ -1,4 +1,3 @@ - from dataclasses import dataclass from pythermogis.workflow.utc.economics import calculate_economics @@ -10,7 +9,7 @@ def calculate_max_pressure( input_data, use_olsthoorn_max_pressure: bool, well_distance: float, - injection_temp: float + injection_temp: float, ) -> float: if use_olsthoorn_max_pressure: max_pres = 0.2 * 0.1 * input_data.depth * 100000 @@ -68,6 +67,7 @@ def calculate_max_pressure( return pressure + @dataclass class PressureOptimizerResults: drawdown_pressure: float @@ -84,6 +84,7 @@ class PressureOptimizerResults: hp_added_power: float production_temp: float + def optimize_pressure( props, input, @@ -92,7 +93,6 @@ def optimize_pressure( injection_temp: float, stimulation_capex: float, ): - pres_tol = 1e-4 * drawdown_pressure pres_step = 1e1 * pres_tol @@ -114,9 +114,15 @@ def optimize_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 + 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 @@ -133,12 +139,22 @@ def optimize_pressure( if flow_results is not None: utc_val = find_utc_for_pres( - props, input, pres, well_distance, injection_temp, stimulation_capex + props, + input, + pres, + well_distance, + injection_temp, + stimulation_capex, ) slope = utc_slope( - props, input, pres, pres_step, - well_distance, injection_temp, stimulation_capex + props, + input, + pres, + pres_step, + well_distance, + injection_temp, + stimulation_capex, ) if slope < utc_val * 1e-5 * 1e-2 * props.tolerance_utc_increase: @@ -186,9 +202,9 @@ def optimize_pressure( 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) + 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( @@ -232,12 +248,20 @@ def utc_slope( stimulation_capex, ): utc1 = find_utc_for_pres( - props, input, pres - pres_step / 2, - well_distance, injection_temp, stimulation_capex + 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 + props, + input, + pres + pres_step / 2, + well_distance, + injection_temp, + stimulation_capex, ) return (utc2 - utc1) / pres_step @@ -272,4 +296,4 @@ def find_utc_for_pres( pump_power_required=flow_results.pump_power_required, ) - return econ.utc.utc \ No newline at end of file + return econ.utc.utc diff --git a/src/pythermogis/workflow/utc/rock.py b/src/pythermogis/workflow/utc/rock.py index 245ef6d96654fd03dd470d3f0888664cfafb8cb8..6ebd2cf6d585bd30ee770cd60c8173e899e032f8 100644 --- a/src/pythermogis/workflow/utc/rock.py +++ b/src/pythermogis/workflow/utc/rock.py @@ -1,7 +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 +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 index 1b4da8c6c0afb88f1b81bab9c1ec4ff37c66bc78..1a75ce909c0ae2f478595e072d0a82cf4a1d7bda 100644 --- a/src/pythermogis/workflow/utc/utc_properties.py +++ b/src/pythermogis/workflow/utc/utc_properties.py @@ -1,18 +1,22 @@ from dataclasses import dataclass, field from pathlib import Path -from typing import NamedTuple, Literal +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(frozen=True) class UTCConfiguration: input_data_dir: str = "" @@ -23,36 +27,70 @@ class UTCConfiguration: 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") - ]) + 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]) + 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 @@ -155,4 +193,4 @@ class UTCConfiguration: use_orc: bool = False heat_exchanger_efficiency: float = 1.0 heat_exchanger_parasitic: float = 0.0 - heat_exchanger_basetemp: float = 10.0 \ No newline at end of file + heat_exchanger_basetemp: float = 10.0 diff --git a/src/pythermogis/workflow/utc/water.py b/src/pythermogis/workflow/utc/water.py index 87f3abdbb03642fc240919bea8cfc8c9a26aaf93..cf16e4ff5744361a10958c52221cf13235ffd305 100644 --- a/src/pythermogis/workflow/utc/water.py +++ b/src/pythermogis/workflow/utc/water.py @@ -1,18 +1,24 @@ 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: +@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 @@ -43,23 +49,21 @@ def density(P: float, T: float, S: float, density_val = density_fresh + S * ( 0.668 + 0.44 * S - + 1e-6 * ( + + 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 - ) + + 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: +def heat_capacity( + P: float, T: float, S: float, salinity_input_in_ppm: bool = False +) -> float: if salinity_input_in_ppm: S *= 1e-6 @@ -76,12 +80,15 @@ def heat_capacity(P: float, T: float, S: float, salinity_input_in_ppm: bool = Fa 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 + + (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 -@njit -def get_salinity(surface_salinity, gradient, depth): - return (surface_salinity + gradient * depth) / 1e6 \ No newline at end of file + diff --git a/src/pythermogis/workflow/utc/well_distance.py b/src/pythermogis/workflow/utc/well_distance.py index 0b9be4e3de48ce5609b9c1586d4c2377809fb9d4..1f09a68bb1a4c27ec066700b229df6ebc0429d4a 100644 --- a/src/pythermogis/workflow/utc/well_distance.py +++ b/src/pythermogis/workflow/utc/well_distance.py @@ -32,7 +32,7 @@ def optimize_well_distance_original( dist_max = props.optim_dist_well_dist_max well_distance = np.mean([dist_min, dist_max]) - for iter_count in range(1000): + 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]) @@ -69,7 +69,8 @@ def optimize_well_distance_original( # If no convergence in 1000 iterations else: print( - f"WARNING: Well distance optimization failed to converge. Final dist={well_distance}" + f"WARNING: Well distance optimization failed to converge." + f" Final dist={well_distance}" ) return well_distance @@ -88,19 +89,22 @@ def f1( ) # --- 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 + 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( @@ -119,4 +123,4 @@ def optimize_well_distance( maxiter=100, args=(props, input, drawdown_pressure, injection_temp), ) - return well_distance \ No newline at end of file + return well_distance