TNO Intern

Commit c241f8a5 authored by Florian Knappers's avatar Florian Knappers
Browse files

added potential light (prh and hip) calculation

parent 0a67a8cf
Loading
Loading
Loading
Loading
Loading
+123 −1
Original line number Diff line number Diff line
@@ -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.46.0-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.63.1-py313h5dce7c4_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.10-h4196e79_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
@@ -371,6 +374,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.46.0-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 +396,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.63.1-py313h2da9318_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 +437,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.10-h37e10c4_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
@@ -3622,6 +3628,40 @@ packages:
  purls: []
  size: 55476
  timestamp: 1727963768015
- conda: https://conda.anaconda.org/conda-forge/linux-64/llvmlite-0.46.0-py313hdd307be_0.conda
  sha256: 0e1bc6ee1c7885cc26c37fcd1a2095169a4e4e148860c600d3f685b6a4f32322
  md5: d99ac09b331711fd12e16323ca8d96e4
  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: 34130706
  timestamp: 1765280056189
- conda: https://conda.anaconda.org/conda-forge/win-64/llvmlite-0.46.0-py313h5c49287_0.conda
  sha256: e0da49a4388c6a906a498766b4636a38d8173aee919d12843ae448c44d78384d
  md5: e0908ac3f4ae20a4957fe1898644f42a
  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: 22885319
  timestamp: 1765280139042
- conda: https://conda.anaconda.org/conda-forge/noarch/locket-1.0.0-pyhd8ed1ab_0.tar.bz2
  sha256: 9afe0b5cfa418e8bdb30d8917c5a6cec10372b037924916f1f85b9f4899a67a6
  md5: 91e27ef3d05cc772ce627e51cff111c4
@@ -4081,6 +4121,57 @@ packages:
  - pkg:pypi/nodeenv?source=hash-mapping
  size: 34574
  timestamp: 1734112236147
- conda: https://conda.anaconda.org/conda-forge/linux-64/numba-0.63.1-py313h5dce7c4_0.conda
  sha256: 3ceba93570814df69969edff3156097dc0e86ccefa2ea2bdfe08f84b2023cf04
  md5: dbdae1a85bb346d57fae63269def955a
  depends:
  - __glibc >=2.17,<3.0.a0
  - _openmp_mutex >=4.5
  - libgcc >=14
  - libstdcxx >=14
  - llvmlite >=0.46.0,<0.47.0a0
  - numpy >=1.22.3,<2.4
  - numpy >=1.23,<3
  - python >=3.13,<3.14.0a0
  - python_abi 3.13.* *_cp313
  constrains:
  - cudatoolkit >=11.2
  - cuda-python >=11.6
  - libopenblas !=0.3.6
  - tbb >=2021.6.0
  - cuda-version >=11.2
  - scipy >=1.0
  license: BSD-2-Clause
  license_family: BSD
  purls:
  - pkg:pypi/numba?source=hash-mapping
  size: 5761715
  timestamp: 1765466811957
- conda: https://conda.anaconda.org/conda-forge/win-64/numba-0.63.1-py313h2da9318_0.conda
  sha256: 93fa3a044429298b5437dfca0b30b5834a8ebcbcae9abadba8c7f4fb83b453e5
  md5: 43901612a748ce58c05f2d68a6e0bf46
  depends:
  - llvmlite >=0.46.0,<0.47.0a0
  - numpy >=1.22.3,<2.4
  - numpy >=1.23,<3
  - 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:
  - cuda-python >=11.6
  - cuda-version >=11.2
  - tbb >=2021.6.0
  - scipy >=1.0
  - libopenblas !=0.3.6
  - cudatoolkit >=11.2
  license: BSD-2-Clause
  license_family: BSD
  purls:
  - pkg:pypi/numba?source=hash-mapping
  size: 5720859
  timestamp: 1765808477325
- conda: https://conda.anaconda.org/conda-forge/linux-64/numpy-2.2.5-py313h17eae1a_0.conda
  sha256: c0a200d0e53a1acbfa1d1e2277e3337ea2aa8cb584448790317a98c62dcaebce
  md5: 6ceeff9ed72e54e4a2f9a1c88f47bdde
@@ -4830,7 +4921,7 @@ packages:
- pypi: ./
  name: pythermogis
  version: 1.2.2
  sha256: baf3b6620dc9442685dc46eb20340b3a59ee6defca113b9c467e0338805e9091
  sha256: 2dff903415a0c622cd66d54c9e49503f3dbfc95cfc43bdc51e2d631c652d9136
  requires_dist:
  - jpype1>=1.5.2,<2
  - xarray==2024.9.0.*
@@ -5217,6 +5308,37 @@ packages:
  - scipy ; extra == 'interp'
  - scipy ; extra == 'all'
  requires_python: '>=3.10'
- conda: https://conda.anaconda.org/conda-forge/linux-64/ruff-0.14.10-h4196e79_0.conda
  noarch: python
  sha256: 997b45ce89554f677e4a30cd1d64565949be9d25c806727c3d844fee0d55d7d2
  md5: ddecdee0806589993d96a950ad51b927
  depends:
  - python
  - libgcc >=14
  - __glibc >=2.17,<3.0.a0
  constrains:
  - __glibc >=2.17
  license: MIT
  license_family: MIT
  purls:
  - pkg:pypi/ruff?source=compressed-mapping
  size: 11472103
  timestamp: 1766094973645
- conda: https://conda.anaconda.org/conda-forge/win-64/ruff-0.14.10-h37e10c4_0.conda
  noarch: python
  sha256: 2ca5b7a6ab0a4cddec16f390e1c2d89b5cacd1780963745e463f1bb2e8ab4893
  md5: 987a8290e4bfd9e42e3c58be4481929c
  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=compressed-mapping
  size: 11908812
  timestamp: 1766095035171
- conda: https://conda.anaconda.org/conda-forge/linux-64/s2n-1.5.21-h7ab7c64_0.conda
  sha256: c8b252398b502a5cc6ea506fd2fafe7e102e7c9e2ef48b7813566e8a72ce2205
  md5: 28b5a7895024a754249b2ad7de372faa
+2 −0
Original line number Diff line number Diff line
@@ -70,6 +70,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.63.1,<0.64"
ruff = ">=0.14.10,<0.15"

[tool.ruff]
line-length = 88
+14 −3
Original line number Diff line number Diff line
@@ -5,6 +5,8 @@ import typer
import xarray as xr

from pythermogis import calculate_doublet_performance
from pythermogis.thermogis_classes.potential_properties import instantiate_potential_properties_builder
from pythermogis.potential.potential import calculate_potential

app = typer.Typer(pretty_exceptions_enable=False)

@@ -68,11 +70,20 @@ def simulate_doublet(

@app.command()
def potential(
        output_file: Path = typer.Option(None, help="A file to output the data to, accepted filetypes are: .nc, .csv, .xlsx, .json"),
        verbose: bool = typer.Option(False, help="print the input data and output data to terminal")
        temperature: float = typer.Option(50, help="The temperature of the aquifer, if not provided this will be calculated using a Temperature gradient of 8°C + 31°C/km with the depth parameter,"),
        thickness: float = typer.Option(500, help="The thickness of the aquifer, units: [m]"),
        ntg: float = typer.Option(0.5, help="The net-to-gross of the aquifer, units: [0-1]"),
        porosity: float = typer.Option(0.5, help="The porosity of the aquifer, between 0-1 (1 = 100%), units: [0-1]"),
        depth: float = typer.Option(500, help="The depth of the aquifer, +ive downwards, units: [m]"),
) -> None:
    """Perform the potential calculation for a given location."""
    print("So much potential..")
    properties = instantiate_potential_properties_builder().build()
    results = calculate_potential(properties, porosity, temperature, depth, thickness, ntg)

    print("\n---potential results---")
    print(results)



if __name__ == "__main__":
    app()
 No newline at end of file
+95 −0
Original line number Diff line number Diff line
def calculate_water_density(
    pressure_pa: float,
    temperature_celsius: float,
    salinity_fraction: float,
) -> float:
    pressure_pa = max(pressure_pa, 1.0)
    temperature_celsius = max(temperature_celsius, 0.0)
    salinity_fraction = max(salinity_fraction, 0.0)

    # Convert pressure to MPa for correlation
    pressure_mpa = pressure_pa * 1e-6

    # Calculate freshwater density (g/cm³)
    temp_squared = temperature_celsius * temperature_celsius
    temp_cubed = temp_squared * temperature_celsius
    pressure_squared = pressure_mpa * pressure_mpa

    density_fresh = 1.0 + 1e-6 * (
        -80.0 * temperature_celsius
        - 3.3 * temp_squared
        + 0.00175 * temp_cubed
        + 489.0 * pressure_mpa
        - 2.0 * temperature_celsius * pressure_mpa
        + 0.016 * temp_squared * pressure_mpa
        - 1.3e-5 * temp_cubed * pressure_mpa
        - 0.333 * pressure_squared
        - 0.002 * temperature_celsius * pressure_squared
    )

    # Calculate seawater density (g/cm³)
    density_g_per_cm3 = density_fresh + salinity_fraction * (
        0.668
        + 0.44 * salinity_fraction  # Changed from salinity_squared
        + 1e-6
        * (
            300.0 * pressure_mpa
            - 2400.0 * pressure_mpa * salinity_fraction
            + temperature_celsius
            * (
                80.0
                + 3.0 * temperature_celsius
                - 3300.0 * salinity_fraction
                - 13.0 * pressure_mpa
                + 47.0 * pressure_mpa * salinity_fraction
            )
        )
    )

    # Convert from g/cm³ to kg/m³
    density_kg_per_m3 = density_g_per_cm3 * 1000.0

    return density_kg_per_m3


def calculate_water_heat_capacity(
    temperature_celsius: float,
    salinity_fraction: float,
) -> float:
    # Apply minimum constraints
    temperature_celsius = max(temperature_celsius, 0.0)
    salinity_fraction = max(salinity_fraction, 0.0)

    # Convert to units used in correlation
    temperature_kelvin = temperature_celsius + 273.15
    salinity_g_per_kg = salinity_fraction * 1e3

    # Calculate heat capacity polynomial
    salinity_squared = salinity_g_per_kg * salinity_g_per_kg
    temp_squared = temperature_kelvin * temperature_kelvin
    temp_cubed = temp_squared * temperature_kelvin

    heat_capacity_kj_per_kg_k = (
        (5.328 - 9.76e-2 * salinity_g_per_kg + 4.04e-4 * salinity_squared)
        + (-6.913e-3 + 7.351e-4 * salinity_g_per_kg - 3.15e-6 * salinity_squared)
        * temperature_kelvin
        + (9.6e-6 - 1.927e-6 * salinity_g_per_kg + 8.23e-9 * salinity_squared)
        * temp_squared
        + (2.5e-9 + 1.666e-9 * salinity_g_per_kg - 7.125e-12 * salinity_squared)
        * temp_cubed
    )

    # Convert from kJ/kg/K to J/kg/K
    heat_capacity_j_per_kg_k = heat_capacity_kj_per_kg_k * 1e3

    return heat_capacity_j_per_kg_k


def calculate_water_salinity(
    surface_salinity: float, gradient: float, depth: float
) -> float:
    return (surface_salinity + gradient * depth) / 1e6


def calculate_hydrostatic_pressure(depth: float) -> float:
    return 9810 * depth + 101325
+94 −0
Original line number Diff line number Diff line
from typing import NamedTuple

from jpype import JClass

from pythermogis.physics.water import (
    calculate_hydrostatic_pressure,
    calculate_water_density,
    calculate_water_heat_capacity,
    calculate_water_salinity,
)


class PotentialOutput(NamedTuple):
    heat_in_place: float
    potential_recoverable_heat: float


def calculate_potential(
    props: JClass,
    porosity: float,
    temperature: float,
    depth: float,
    thickness: float,
    ntg: float,
) -> PotentialOutput:
    """
    Mimics the java potential calculation.
    For now, it just works for a single location.
    In a later stage we might need to make it such that the input is a xr dataarray,
    and perform this calculation for each cell.
    Besides, this potential calculation basically only calculates the hip and prh.
    To calculate the actual potential like in the java code, the input should contain
    multiple p-value values.
    """
    mid_depth = depth + thickness / 2

    hydrostatic_pressure = calculate_hydrostatic_pressure(mid_depth)
    salinity = calculate_water_salinity(
        props.salinitySurface(), props.salinityGradient(), mid_depth
    )
    water_density = calculate_water_density(hydrostatic_pressure, temperature, salinity)
    water_heat_capacity = calculate_water_heat_capacity(temperature, salinity)
    gamma = calculate_gamma(
        porosity, props.cpRock(), props.rhoRock(), water_density, water_heat_capacity
    )

    heat_in_place = calculate_heat_in_place(
        gamma, temperature, props.surfaceTemp(), thickness, ntg
    )
    potential_recoverable_heat = calculate_potential_recoverable_heat(
        gamma, temperature, props.dhReturnTemp(), thickness, ntg
    )

    return PotentialOutput(
        heat_in_place,
        potential_recoverable_heat,
    )


def calculate_heat_in_place(
    gamma: float,
    temperature: float,
    surface_temperature: float,
    thickness: float,
    ntg: float,
) -> float:
    hip = gamma * (temperature - surface_temperature) * thickness * ntg * 1e-9
    hip = max(0, hip)
    return hip


def calculate_potential_recoverable_heat(
    gamma: float,
    temperature: float,
    dh_return_temperature: float,
    thickness: float,
    ntg: float,
) -> float:
    prh = gamma * (temperature - dh_return_temperature) * thickness * ntg * 1e-9 / 3
    prh = max(0, prh)
    return prh


def calculate_gamma(
    porosity: float,
    rock_heat_capacity: float,
    rock_density: float,
    water_density: float,
    water_heat_capacity: float,
) -> float:
    return (
        porosity * water_density * water_heat_capacity
        + (1 - porosity) * rock_heat_capacity * rock_density
    )
Loading