TNO Intern

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

remove p dim from trans calc

parent 9c6c9adc

src/pythermogis/transmissivity/calculate_thick_perm_trans.py

+6 −12
Original line number Diff line number Diff line
@@ -59,8 +59,8 @@ def _tpkt_kernel( 
    thickness_sd,

    ln_perm_mean,

    ln_perm_sd,

    p_values,

    n_samples,

    p_value,

    n_samples=10_000,

):

    if (

        np.isnan(thickness_mean)
@@ -68,19 +68,13 @@ def _tpkt_kernel( 
        or np.isnan(ln_perm_mean)

        or np.isnan(ln_perm_sd)

    ):

        n = len(p_values)

        return (

            np.full(n, np.nan),

            np.full(n, np.nan),

            np.full(n, np.nan),

        )

        return np.nan, np.nan, np.nan



    p = np.asarray(p_values)

    q = 1.0 - p / 100.0

    q = 1.0 - p_value / 100.0



    # Thickness

    if thickness_sd == 0:

        thickness_p = np.full(p.size, thickness_mean)

        thickness_p = thickness_mean

        thickness_samples = np.full(n_samples, thickness_mean)

    else:

        t_dist = stats.norm(thickness_mean, thickness_sd)
@@ -93,7 +87,7 @@ def _tpkt_kernel( 


    # Transmissivity via sampling

    tr_samples = np.sort(np.exp(k_dist.rvs(n_samples) + np.log(thickness_samples)))

    idx = (q * n_samples).astype(int)

    idx = int(q * n_samples)

    transmissivity_p = tr_samples[idx]



    return thickness_p, permeability_p, transmissivity_p
 
 No newline at end of file

src/pythermogis/workflow/utc/utc.py

+7 −8
Original line number Diff line number Diff line
@@ -144,10 +144,9 @@ def compute_results_for_aquifer( 
    # TODO: Maybe save the results as datasets iso dataarrays, and have p value be a dim?



    # calc transmissivity(_with_ntg)

    p_values = xr.DataArray([50], dims="p")

    hydro = compute_hydro_properties_mc(aquifer_ds, p_values)

    transmissivity = hydro["transmissivity"].squeeze("p", drop=True)

    transmissivity_with_ntg = hydro["transmissivity_with_ntg"].squeeze("p", drop=True)

    hydro = compute_hydro_properties_mc(aquifer_ds, 90)

    transmissivity = hydro["transmissivity"]

    transmissivity_with_ntg = hydro["transmissivity_with_ntg"]

    aquifer_ds["transmissivity"] = transmissivity

    aquifer_ds["transmissivity_with_ntg"] = transmissivity_with_ntg
@@ -221,7 +220,7 @@ def compute_results_for_aquifer( 


def compute_hydro_properties_mc(

    aquifer_ds: xr.Dataset,

    p_values: xr.DataArray,

    p_value: float,

    n_samples: int = 10_000,

) -> xr.Dataset:
@@ -236,10 +235,10 @@ def compute_hydro_properties_mc( 
        thickness_sd,

        ln_perm_mean,

        ln_perm_sd,

        p_values,

        p_value,

        n_samples,

        input_core_dims=[[], [], [], [], ["p"], []],

        output_core_dims=[["p"], ["p"], ["p"]],

        input_core_dims=[[], [], [], [], [], []],

        output_core_dims=[[], [], []],

        vectorize=True,

        dask="parallelized",

        output_dtypes=[np.float64, np.float64, np.float64],