Resolve "add potential (light) calculation"

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

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

    )

from jpype import JClass

from pythermogis.thermogis_classes.jvm_start import start_jvm

def instantiate_potential_properties_builder() -> JClass:

    """

    Create and return a ``PotentialPropertiesBuilder`` Java object.

    The builder wraps the ThermoGIS **PotentialPropertiesBuilder** Java class, which can

    * parse ThermoGIS scenario XML files, and

    * programmatically modify individual scenario parameters prior to running a simulation.

    Returns

    -------

    potential_properties_builder : jpype.JClass

        An instantiated ``PotentialPropertiesBuilder`` that can be used to build or modify

        ``PotentialProperties`` objects for ThermoGIS simulations.

    """

    start_jvm()

    return JClass("thermogis.properties.builders.PotentialPropertiesBuilder")()

import json

from pathlib import Path

import timeit

import numpy as np

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)

        else:

            typer.echo(typer.style("Error: Invalid output file type provided, currently supported: .nc, .csv, .json", fg=typer.colors.RED, bold=True))

@app.command()

def potential(

        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."""

    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

import numpy as np

from pythermogis.physics.water import (

    calculate_hydrostatic_pressure,

    calculate_water_salinity,

    calculate_water_heat_capacity,

    calculate_water_density,

)

def test_calculate_hydrostatic_pressure():

    pres = calculate_hydrostatic_pressure(2089.5)

    assert np.isclose(pres, 2.059932e7, 1e-5)

def test_calculate_water_salinity():

    sal = calculate_water_salinity(0, 46.67, 2089.5)

    assert np.isclose(sal, 0.097516965, 1e-5)

def test_calculate_water_heat_capacity():

    whc = calculate_water_heat_capacity(

        74.7745,

        0.097516965,

    )

    assert np.isclose(whc, 3748.4475341351704, 1e-5)

def test_calculate_water_density():

    whc = calculate_water_density(

        2.059932e7,

        74.7745,

        0.097516965,

    )

    assert np.isclose(whc, 1052.8835058924815, 1e-5)