TNO Intern

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

rff

parent b91edcdd

src/pythermogis/transmissivity/calculate_thick_perm_trans.py

+13 −11
Original line number Diff line number Diff line 
import numpy as np

from scipy import stats

import xarray as xr

from numba import njit, prange

from numpy.typing import NDArray





def generate_thickness_permeability_transmissivity_for_pvalues(thickness_mean: float, thickness_sd: float, ln_permeability_mean: float, ln_permeability_sd: float, p_values: xr.DataArray, nSamples: int = 10000) -> float:

    """
@@ -54,17 +57,16 @@ def generate_thickness_permeability_transmissivity_for_pvalues(thickness_mean: f 


    return thickness_pvalues, permeability_pvalues, transmissivity_pvalues_sampled



from numba import njit, prange



@njit(parallel=True)

def calculate_transmissivity(

        thickness_mean,

        thickness_sd,

        ln_perm_mean,

        ln_perm_sd,

        q,  # quantile (e.g., 0.05 for p95)

        z,  # precomputed ndtri(q)

        n_samples,

        thickness_mean: NDArray[np.float64],

        thickness_sd: NDArray[np.float64],

        ln_perm_mean: NDArray[np.float64],

        ln_perm_sd: NDArray[np.float64],

        quantile: float,

        z_score: float,

        n_samples: int,

):

    shape = thickness_mean.shape

    tm_flat = thickness_mean.ravel()
@@ -73,7 +75,7 @@ def calculate_transmissivity( 
    lps_flat = ln_perm_sd.ravel()



    n_cells = tm_flat.size

    idx = int(q * n_samples)

    idx = int(quantile * n_samples)



    thickness_p = np.empty(n_cells)

    permeability_p = np.empty(n_cells)
@@ -108,8 +110,8 @@ def calculate_transmissivity( 
        transmissivity_p[i] = tr_samples[idx]



        # Analytical percentiles

        thickness_p[i] = tm + z * ts if ts > 0 else tm

        permeability_p[i] = np.exp(lpm + z * lps)

        thickness_p[i] = tm + z_score * ts if ts > 0 else tm

        permeability_p[i] = np.exp(lpm + z_score * lps)



    return (

        thickness_p.reshape(shape),

src/pythermogis/workflow/utc/utc.py

+6 −5
Original line number Diff line number Diff line
@@ -3,8 +3,11 @@ 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.transmissivity.calculate_thick_perm_trans import (

    calculate_transmissivity,

)

from pythermogis.workflow.utc.configuration import UTCConfiguration

from pythermogis.workflow.utc.doublet import DoubletInput, calculate_doublet_performance
@@ -199,14 +202,12 @@ def compute_results_for_aquifer( 


    return aquifer_ds



from scipy.special import ndtri



def compute_hydro_properties_mc(

    aquifer_ds: xr.Dataset,

    p_value: float,

    n_samples: int = 10_000,

) -> xr.Dataset:



    start = time.perf_counter()



    q = 1.0 - p_value / 100.0
@@ -234,7 +235,7 @@ def compute_hydro_properties_mc( 
            "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,  # this one is already a DataArray

            "transmissivity_with_ntg": transmissivity_with_ntg,

        },

        coords=aquifer_ds.coords,

    )