)

    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"),

            PropertyGridInfo("permeability", False, "__k.nc"),

            PropertyGridInfo("permeability_lnsd", False, "__k_ln_sd.nc"),

            PropertyGridInfo("porosity", False, "__phi.nc"),

            PropertyGridInfo("thickness", False, "__thick.nc"),

            PropertyGridInfo("thickness_sd", False, "__thick_sd.nc"),

            PropertyGridInfo("depth", False, "__top.nc"),

            PropertyGridInfo("ntg", False, "__ntg.nc"),

            PropertyGridInfo("temperature", True, "__temperature.nc"),

            PropertyGridInfo("hc_accum", False, "__hc_accum.nc"),

            PropertyGridInfo("boundary_shapefile", True, "__BoundaryShapefile.nc"),

        ]

    )

    copy_aquifer_files_info: list[AquiferFile] = field(

            AquiferFile("__points_QC.shp", "__points_QC.shp"),

        ]

    )

    scenario: str = "basecase"

    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, 90.0])

    # temp_voxet: 'Voxet' = None

    surface_temperature: float = 10.0

    temp_gradient: float = 31.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

    econ_lifetime_years: int = 30

    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

    anisotropy: float = 3.0

    filter_fraction: float = 0.8

    salinity_surface: float = 0.0

    salinity_gradient: float = 46.67

    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

    imposed_cop: float = 4.0

    hp_capex: float = 200.0

    hp_opex: float = 10.0

    hp_alternative_heating_price: float = 2.8

    viscosity_mode: ViscosityMode = "batzlewang"

    # Economical data

    economic_lifetime: int = 15

    economic_lifetime: int = 30

    drilling_time: int = 2

    tax_rate: float = 0.25

    interest_loan: float = 0.05

    inflation: float = 0.02

    equity_return: float = 0.07

    inflation: float = 0.015

    equity_return: float = 0.145

    debt_equity: float = 0.8

    ecn_eia: float = 0.0

    tolerance_utc_increase: float = 0.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_scaling: float = 1.3

    well_cost_const: float = 0.25

    well_cost_z: float = 700.0

    well_cost_z2: float = 0.2

from pythermogis.workflow.utc.flow import calculate_volumetric_flow

if TYPE_CHECKING:

    from pythermogis.workflow.utc.configuration import UTCConfiguration

    from pythermogis.workflow.utc.doublet import DoubletInput

    from pythermogis.workflow.utc.utc_properties import UTCConfiguration

def calculate_cooling_temperature(

from typing import NamedTuple

from pythermogis.utils.timer import print_time

from pythermogis.workflow.utc.configuration import UTCConfiguration

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 (

    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

from pythermogis.workflow.utc.water import get_salinity

if TYPE_CHECKING:

    from pythermogis.workflow.utc.configuration import UTCConfiguration

    from pythermogis.workflow.utc.doublet import DoubletInput

    from pythermogis.workflow.utc.utc_properties import UTCConfiguration

INCH_SI = 0.0254

from pythermogis.workflow.utc.doublet_utils import get_along_hole_length

if TYPE_CHECKING:

    from pythermogis.workflow.utc.configuration import UTCConfiguration

    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

from pythermogis.workflow.utc.heatpump import calculate_heat_pump_performance

if TYPE_CHECKING:

    from pythermogis.workflow.utc.configuration import UTCConfiguration

    from pythermogis.workflow.utc.doublet import DoubletInput

    from pythermogis.workflow.utc.utc_properties import UTCConfiguration

class VolumetricFlowResults(NamedTuple):

)

if TYPE_CHECKING:

    from pythermogis.workflow.utc.configuration import UTCConfiguration

    from pythermogis.workflow.utc.doublet import DoubletInput

    from pythermogis.workflow.utc.utc_properties import UTCConfiguration

class HeatPumpPerformanceResults(NamedTuple):

from pythermogis.workflow.utc.flow import calculate_volumetric_flow

if TYPE_CHECKING:

    from pythermogis.workflow.utc.configuration import UTCConfiguration

    from pythermogis.workflow.utc.doublet import DoubletInput

    from pythermogis.workflow.utc.utc_properties import UTCConfiguration

def calculate_max_pressure(

import numpy as np

from scipy.stats import norm

def lognormal_quantile(mean, sd, p):

    z = norm.ppf(1 - p / 100)

    return mean * np.exp(sd * z)

def normal_quantile(mean, sd, p):

    z = norm.ppf(1 - p / 100)

    return mean + sd * z

import time

from pathlib import Path

import numpy as np

import xarray as xr

from scipy.special import ndtri

from pythermogis.transmissivity.calculate_thick_perm_trans import (

    calculate_transmissivity,

)

from pythermogis.workflow.utc.configuration import UTCConfiguration

from pythermogis.workflow.utc.doublet import (

    DoubletInput,

    DoubletOutput,

    calculate_doublet_performance,

)

OUTPUT_NAMES = list(DoubletOutput.__annotations__.keys())

NAN_OUTPUTS = DoubletOutput(*[np.nan] * len(OUTPUT_NAMES)) # type: ignore[arg-type]

ZERO_OUTPUTS = DoubletOutput(*[0.0] * len(OUTPUT_NAMES)) # type: ignore[arg-type]

def run_utc_workflow(config: UTCConfiguration) -> xr.DataTree:

    root = xr.DataTree()

    base = Path(config.input_data_dir)

    temperature_model = xr.open_dataset(

        base / "NetherlandsTemperatureModel_extended.nc"

    ).transpose("y", "x", "z")

    for aquifer in config.aquifers:

        aquifer_ds = xr.Dataset()

        for gridinfo in config.property_grid_infos:

            if gridinfo.optional:

                continue

            filename = f"{aquifer}{gridinfo.postfix}"

            fullpath = base / filename

            if not fullpath.exists():

                raise FileNotFoundError(f"Required grid not found: {fullpath}")

            ds = xr.open_dataset(fullpath)

            varnames = [v for v in ds.data_vars if v != "oblique_stereographic"]

            if len(varnames) != 1:

                raise ValueError(

                    f"Expected exactly one property variable in"

                    f" {fullpath}, found {varnames}"

                )

            varname = varnames[0]

            aquifer_ds[gridinfo.name] = ds[varname]

        root[aquifer] = xr.DataTree(dataset=aquifer_ds, name=aquifer)

    for aquifer in root.children:

        print(f"\nProcessing aquifer: {aquifer}")

        start = time.perf_counter()

        aquifer_ds = root[aquifer].to_dataset()

        updated_ds = compute_results_for_aquifer(

            aquifer_ds, aquifer, config, temperature_model

        )

        root[aquifer] = xr.DataTree(dataset=updated_ds, name=aquifer)

        end = time.perf_counter()

        elapsed = end - start

        print(f"Time required for aquifer {aquifer}: {elapsed:.2f} seconds.")

    if config.results_dir != "":

        zarr_path = f"{config.results_dir}/ROSL_ROSLU_test.zarr"

        root.to_zarr(zarr_path, mode="w")

    return root

def compute_results_for_aquifer(

    aquifer_ds: xr.Dataset,

    aquifer_name: str,

    config: UTCConfiguration,

    temperature_model: xr.Dataset,

):

    # apply hc mask

    hc_mask = aquifer_ds["hc_accum"] != 0

    vars_to_mask = [

        "thickness",

        "ntg",

        "depth",

        "porosity",

    ]

    for v in vars_to_mask:

        aquifer_ds[v] = aquifer_ds[v].where(hc_mask)

    # calc temperature grid

    mid_depth = -(aquifer_ds["depth"] + 0.5 * aquifer_ds["thickness"])

    temp_xy_aligned = temperature_model["temperature"].interp(

        x=aquifer_ds.x, y=aquifer_ds.y, method="linear"

    )

    temperature = temp_xy_aligned.interp(z=mid_depth, method="linear")

    aquifer_ds["temperature"] = temperature

    for p_value in config.p_values:

        # calc transmissivity(_with_ntg)

        stoch_results = apply_stochastics(aquifer_ds, p_value)

        aquifer_ds["transmissivity"] = stoch_results["transmissivity"]

        aquifer_ds["transmissivity_with_ntg"] = stoch_results["transmissivity_with_ntg"]

        aquifer_ds["thickness"] = stoch_results["thickness"]

        aquifer_ds["permeability"] = stoch_results["permeability"]

        results = xr.apply_ufunc(

            cell_calculation,

            1.0e30,

            aquifer_ds["thickness"],

            aquifer_ds["transmissivity"],

            aquifer_ds["transmissivity_with_ntg"],

            aquifer_ds["ntg"],

            aquifer_ds["depth"],

            aquifer_ds["porosity"],

            aquifer_ds["temperature"],

            kwargs={"config": config},

            input_core_dims=[[]] * 8,

            output_core_dims=[[] for _ in OUTPUT_NAMES],

            vectorize=True,

            dask="parallelized",

            output_dtypes=[np.float64] * len(OUTPUT_NAMES),

        )

        for name, result in zip(OUTPUT_NAMES, results, strict=True):

            aquifer_ds[f"{name}_p{p_value}"] = result

    aquifer_ds.compute()

    print(f"Computed results for {aquifer_name}")

    return aquifer_ds

def apply_stochastics(

    aquifer_ds: xr.Dataset,

    p_value: float,

    n_samples: int = 10_000,

) -> xr.Dataset:

    q = 1.0 - p_value / 100.0

    z = ndtri(q)

    thickness_p, permeability_p, transmissivity_p = calculate_transmissivity(

        aquifer_ds["thickness"].values,

        aquifer_ds["thickness_sd"].values,

        np.log(aquifer_ds["permeability"].values),

        aquifer_ds["permeability_lnsd"].values,

        q,

        z,

        n_samples,

    )

    transmissivity_with_ntg = transmissivity_p * aquifer_ds["ntg"] / 1000.0

    transmissivity_with_ntg = transmissivity_with_ntg.where(

        transmissivity_with_ntg >= 0

    )

    return xr.Dataset(

        {

            "thickness": (aquifer_ds["thickness"].dims, thickness_p),

            "permeability": (aquifer_ds["thickness"].dims, permeability_p),

            "transmissivity": (aquifer_ds["thickness"].dims, transmissivity_p),

            "transmissivity_with_ntg": transmissivity_with_ntg,

        },

        coords=aquifer_ds.coords,

    )

def cell_calculation(

    unknown_input_value: float,

    thickness: float,

    transmissivity: float,

    transmissivity_with_ntg: float,

    ntg: float,

    depth: float,

    porosity: float,

    temperature: float,

    config: UTCConfiguration,

) -> DoubletOutput:

    inputs = [

        unknown_input_value,

        thickness,

        transmissivity,

        transmissivity_with_ntg,

        ntg,

        depth,

        porosity,

        temperature,

    ]

    if any(np.isnan(v) for v in inputs):

        return NAN_OUTPUTS

    if transmissivity_with_ntg < config.kh_cutoff:

        return ZERO_OUTPUTS

    doublet_input = DoubletInput(

        unknown_input_value=unknown_input_value,

        thickness=thickness,

        transmissivity=transmissivity,

        transmissivity_with_ntg=transmissivity_with_ntg,

        ntg=ntg,

        depth=depth,

        porosity=porosity / 100,

        temperature=temperature,

    )

    try:

        result = calculate_doublet_performance(config, doublet_input)

    except ValueError:

        return NAN_OUTPUTS

    if result is None:

        return NAN_OUTPUTS

    if result.utc > 1000:

        result = result._replace(utc=1000)

    return result

from pythermogis.workflow.utc.flow import calculate_volumetric_flow

if TYPE_CHECKING:

    from pythermogis.workflow.utc.configuration import UTCConfiguration

    from pythermogis.workflow.utc.doublet import DoubletInput

    from pythermogis.workflow.utc.utc_properties import UTCConfiguration

def optimize_well_distance_original(

from pythermogis.workflow.utc.doublet import DoubletInput

from pythermogis.workflow.utc.flow import calculate_volumetric_flow

from pythermogis.workflow.utc.utc_properties import UTCConfiguration

from pythermogis.workflow.utc.configuration import UTCConfiguration

def test_calculate_volumetric_flow():

from pythermogis.workflow.utc.doublet import DoubletInput, calculate_doublet_performance

from pythermogis.workflow.utc.utc_properties import UTCConfiguration

from pythermogis.workflow.utc.configuration import UTCConfiguration

def test_calculate_doublet_performance_precise():

    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.utc, 5.885437310055949, 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)

    assert np.isclose(result.var_opex, -16.189716379776353, rtol=rtol)

    assert np.isclose(result.fixed_opex, -24.199008115259968, rtol=rtol)

def test_calculate_doublet_performance_approximate():

    # Arrange: instantiate default UTCConfiguration

    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.utc, 5.885437310055949, 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)

    assert np.isclose(result.var_opex, -16.189716379776353, rtol=rtol)

    assert np.isclose(result.fixed_opex, -24.199008115259968, rtol=rtol)

def test_calculate_doublet_performance():

    # 85ish it/s

        f"{n_sims} simulations took: {time_elapsed:.1f} seconds\n"

        f"{n_sims/time_elapsed:.1f} simulations per second"

    )

def test_why_does_this_fail():

    props = UTCConfiguration(

        segment_length=1,

        dh_return_temp=40

    )

    # Create a minimal valid DoubletInput

    input_data = DoubletInput(

        unknown_input_value=-9999,

        thickness=5.259,

        transmissivity=1239.250,

        transmissivity_with_ntg=1.032,

        ntg=0.8331,

        depth=852.969,

        porosity=0.22168,

        temperature=39.30,

    )

    # Act

    results = calculate_doublet_performance(props, input_data)

    print(results)

from pythermogis.workflow.utc.doublet import DoubletInput

from pythermogis.workflow.utc.doubletcalc import doubletcalc

from pythermogis.workflow.utc.utc_properties import UTCConfiguration

from pythermogis.workflow.utc.configuration import UTCConfiguration

def test_doubletcalc():

from pythermogis.workflow.utc.utc_properties import UTCConfiguration

from pythermogis.workflow.utc.configuration import UTCConfiguration

def test_utc_configuration_instantiation():

    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.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.p_values == [10.0, 30.0, 50.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 cfg.property_grid_infos[0].postfix.endswith(".nc")

    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

    assert cfg.economic_lifetime == 30

# import matplotlib.pyplot as plt

#

# from pythermogis.workflow.utc.configuration import UTCConfiguration

# from pythermogis.workflow.utc.utc import run_utc_workflow

#

# def test_rosl_roslu():

#     """

#     Integration test that runs the whole ROSL_ROSLU aquifer.

#     The input is identical to the java 2.5.3 input.

#     The output is compared to the java output.

#     """

#     config = UTCConfiguration(

#         input_data_dir="C:/Users/knappersfy/work/thermogis/pythermogis/workdir/ROSL_ROSLU/input",

#         results_dir="C:/Users/knappersfy/work/thermogis/pythermogis/workdir/ROSL_ROSLU/output",

#         aquifers=["ROSL_ROSLU"],

#         segment_length=1,

#         p_values=[50] # For now only calculating this value, because it is way faster than 50

#     )

#     result = run_utc_workflow(config)

#

#     result.ROSL_ROSLU.utc_p50.plot() # should match p value in config

#     plt.show()

#

#     print(result)

 No newline at end of file

import numpy as np

from pythermogis.workflow.utc.doublet import DoubletInput

from pythermogis.workflow.utc.utc_properties import UTCConfiguration

from pythermogis.workflow.utc.configuration import UTCConfiguration

from pythermogis.workflow.utc.well_distance import (

    optimize_well_distance,

    optimize_well_distance_original,