TNO Intern

Commit bbb30ef7 authored by jjflorian's avatar jjflorian
Browse files

change some default settings and add hc accum mask.

parent 6188032e

src/pythermogis/workflow/utc/configuration.py

+19 −19
Original line number Diff line number Diff line
@@ -62,15 +62,15 @@ class UTCConfiguration: 
    )

    property_grid_infos: list[PropertyGridInfo] = field(

        default_factory=lambda: [

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

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

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

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

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

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

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

            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", True, "__hc_accum.nc"),

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

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

        ]

    )
@@ -82,14 +82,14 @@ class UTCConfiguration: 
            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]

        default_factory=lambda: [10.0, 30.0, 50.0, 90.0]

    )

    # temp_voxet: 'Voxet' = None

    surface_temperature: float = 10.0
@@ -119,13 +119,13 @@ class UTCConfiguration: 
    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
@@ -160,19 +160,19 @@ class UTCConfiguration: 
    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
@@ -182,7 +182,7 @@ class UTCConfiguration: 
    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

src/pythermogis/workflow/utc/stats.py

0 → 100644
+10 −0
Original line number Diff line number Diff line 
from scipy.stats import norm

import numpy as np



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
 
 No newline at end of file

src/pythermogis/workflow/utc/utc.py

+46 −11
Original line number Diff line number Diff line
@@ -8,10 +8,12 @@ from pythermogis.workflow.utc.configuration import UTCConfiguration 
from pythermogis.workflow.utc.doublet import DoubletInput, calculate_doublet_performance

from pythermogis.workflow.utc.stats import lognormal_quantile, normal_quantile



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

def run_utc_workflow(config: UTCConfiguration, coarsen_dataset_by: int | None = None ) -> xr.DataTree:

    root = xr.DataTree()

    base = Path(config.input_data_dir)



    temperature_model = xr.open_dataset(base / "NetherlandsTemperatureModel_extended.nc")



    for aquifer in config.aquifers:

        aquifer_ds = xr.Dataset()
@@ -38,7 +40,8 @@ def run_utc_workflow(config: UTCConfiguration) -> xr.DataTree: 


            aquifer_ds[gridinfo.name] = ds[varname]



        aquifer_ds = aquifer_ds.coarsen(x=3, y=3, boundary="trim").mean()

        if coarsen_dataset_by is not None:

            aquifer_ds = aquifer_ds.coarsen(x=coarsen_dataset_by, y=coarsen_dataset_by, boundary="trim").mean()



        root[aquifer] = xr.DataTree(dataset=aquifer_ds, name=aquifer)
@@ -49,7 +52,7 @@ def run_utc_workflow(config: UTCConfiguration) -> xr.DataTree: 
        start = time.perf_counter()



        aquifer_ds = root[aquifer].to_dataset()

        updated_ds = compute_results_for_aquifer(aquifer_ds, aquifer, config)

        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()
@@ -57,14 +60,15 @@ def run_utc_workflow(config: UTCConfiguration) -> xr.DataTree: 


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



    zarr_path = "../ROSL_ROSLU_test.zarr"

    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

        aquifer_ds: xr.Dataset, aquifer_name: str, config: UTCConfiguration, temperature_model: xr.Dataset

):

    output_names = [

        "power",
@@ -133,16 +137,46 @@ def compute_results_for_aquifer( 
    )



    # TODO: loop over pvalues and save grids with pvalue noted

    # Maybe save the results as datasets iso dataarrays, and have p value be a dim?

    # TODO: Maybe save the results as datasets iso dataarrays, and have p value be a dim?



    # calc transmissivity(_with_ntg)

    pvalue = 50

    hydro = compute_hydro_properties(aquifer_ds, pvalue)



    transmissivity = hydro["transmissivity"]

    transmissivity_with_ntg = hydro["transmissivity_with_ntg"]



    aquifer_ds["transmissivity"] = transmissivity

    aquifer_ds["transmissivity_with_ntg"] = transmissivity_with_ntg



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

    )



    # use hc_accum mask

    hc_mask = aquifer_ds["hc_accum"] != 0



    vars_to_mask = [

        "thickness",

        "ntg",

        "depth",

        "porosity",

        "transmissivity",

        "transmissivity_with_ntg",

    ]



    for v in vars_to_mask:

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





    aquifer_ds["temperature"] = temperature



    results = xr.apply_ufunc(

        vectorized_calc,

        -9999,
@@ -152,7 +186,7 @@ def compute_results_for_aquifer( 
        aquifer_ds["ntg"],

        aquifer_ds["depth"],

        aquifer_ds["porosity"],

        75.0,   # TODO: replace when temperature model is available

        temperature,

        output_core_dims=[[] for _ in output_names],

        vectorize=False,

        dask="parallelized",
@@ -196,7 +230,8 @@ def compute_hydro_properties(aquifer_ds: xr.Dataset, pvalue: float) -> xr.Datase 
if __name__ == "__main__":

    config = UTCConfiguration(

        input_data_dir="C:/Users/knappersfy/work/thermogis/pythermogis/workdir/InputData",

        results_dir=".",

        aquifers=["ROSL_ROSLU"],

        segment_length=1,

    )

    run_utc_workflow(config)
 
    run_utc_workflow(config, 2)

tests/utc/test_rosl_roslu.py

0 → 100644
+18 −0
Original line number Diff line number Diff line 
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_ROSLUoutput",

        aquifers=["ROSL_ROSLU"],

        segment_length=1,

    )

    result = run_utc_workflow(config)



    print(result)
 
 No newline at end of file