TNO Intern

Draft: Python rewrite of java code

Florian Knappersrequested to merge
python-rewrite-of-java-code into main

Merge request reports

src/pythermogis/thermogis_classes/doublet.py

+62 −23
Changes for src/pythermogis/thermogis_classes/doublet.py: 62 added lines, 23 removed lines.
Original line number Diff line number Diff line 
import numpy as np

import xarray as xr

from jpype import JClass

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

from pythermogis.workflow.utc.configuration import UTCConfiguration





def simulate_doublet(

    output_data: xr.Dataset,
@@ -102,11 +105,11 @@ def calculate_performance_of_single_location(mask: float, depth: float, thicknes 
    if not np.isnan(mask) or temperature < utc_properties.minProdTemp():

        return (mask_value,) * 14



    # Instantiate ThermoGIS doublet

    use_java_backend = True

    if use_java_backend:

        doublet = instantiate_thermogis_doublet(utc_properties, rng_seed)



    DoubletInput = JClass("thermogis.calc.utc.doublet.records.DoubletInput")

    input = DoubletInput(

        JavaDoubletInput = JClass("thermogis.calc.utc.doublet.records.DoubletInput")

        input = JavaDoubletInput(

                -9999.0, # unknowninput

                thickness,

                transmissivity,
@@ -117,12 +120,49 @@ def calculate_performance_of_single_location(mask: float, depth: float, thicknes 
                temperature,

                None, # ates input

            )



    # The Java routine which calculates DoubletPerformance, for more detail on the simulation inspect the Java source code

    results = doublet.calculateDoubletPerformance(input)

        # The Java routine which calculates DoubletPerformance, for more detail on the

        # simulation inspect the Java source code

        results_java = doublet.calculateDoubletPerformance(input)

        results = DoubletOutput(

            results_java.power(),

            results_java.hppower(),

            results_java.capex(),

            0.0,

            0.0,

            results_java.opex(),

            results_java.utc(),

            results_java.npv(),

            results_java.hprod(),

            results_java.cop(),

            results_java.cophp(),

            results_java.pres(),

            results_java.flow(),

            results_java.welld(),

            results_java.breakthrough(),

            results_java.cooling(),

            results_java.productionTemp(),

            results_java.injectionTemp(),

        )

    else:

        # Arrange: instantiate default UTCConfiguration

        props = UTCConfiguration(segment_length=1)

        # Create a minimal valid DoubletInput

        input_data = DoubletInput(

            unknown_input_value=-9999.0,

            thickness=thickness,

            transmissivity=transmissivity,

            transmissivity_with_ntg=transmissivity_with_ntg,

            ntg=ntg,

            depth=depth,

            porosity=porosity,

            temperature=temperature,

        )

        # Act

        #perform one simulation to compile all the numba functions before timing

        results = calculate_doublet_performance(props, input_data)



    # If calculation was not successful, return mask value

    if results.utc() == -9999.0:

    if results.utc == -9999.0:

        return (mask_value,) * 14



    # calculate net-present-value using the utc-cutoffs
@@ -131,24 +171,23 @@ def calculate_performance_of_single_location(mask: float, depth: float, thicknes 
    else:

        utc_cut = utc_properties.utcCutoff()



    hprod = results.hprod()

    npv = 1e-6 * (utc_cut - results.utc()) * 3.6 * hprod * (1 - utc_properties.taxRate())

    npv = 1e-6 * (utc_cut - results.utc) * 3.6 * results.hprod * (1 - utc_properties.taxRate())



    # get values from doublet

    output_values = {"power": results.power(),

                     "heat_pump_power": results.hppower(),

                     "capex": results.capex(),

                     "opex": results.opex(),

                     "utc": results.utc(),

    output_values = {"power": results.power,

                     "heat_pump_power": results.hppower,

                     "capex": results.capex,

                     "opex": results.opex,

                     "utc": results.utc,

                     "npv": npv,

                     "hprod": hprod,

                     "cop": results.cop(),

                     "cophp": results.cophp(),

                     "pres": results.pres(),

                     "flow_rate": results.flow(),

                     "welld": results.welld(),

                     "inj_temp": results.injectionTemp(),

                     "prd_temp": results.productionTemp(),

                     "hprod": results.hprod,

                     "cop": results.cop,

                     "cophp": results.cophp,

                     "pres": results.pres,

                     "flow_rate": results.flow,

                     "welld": results.welld,

                     "inj_temp": results.injection_temp,

                     "prd_temp": results.production_temp,

                     }



    # Reset doublet variables for next calculation

src/pythermogis/transmissivity/calculate_thick_perm_trans.py

+81 −0
Changes for src/pythermogis/transmissivity/calculate_thick_perm_trans.py: 81 added lines, 0 removed lines.
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:

    """
@@ -53,3 +56,81 @@ def generate_thickness_permeability_transmissivity_for_pvalues(thickness_mean: f 
    transmissivity_pvalues_sampled = transmissivity_samples[sample_indexes.astype(int)]



    return thickness_pvalues, permeability_pvalues, transmissivity_pvalues_sampled



@njit(parallel=True)

def calculate_transmissivity(

        thickness_mean: NDArray[np.float64],

        thickness_sd: NDArray[np.float64],

        ln_perm_mean: NDArray[np.float64],

        ln_perm_sd: NDArray[np.float64],

        quantiles: NDArray[np.float64],

        z_scores: NDArray[np.float64],

        n_samples: int,

):

    shape = thickness_mean.shape

    n_pvalues = len(quantiles)



    # Flatten spatial dimensions

    tm_flat = thickness_mean.ravel()

    ts_flat = thickness_sd.ravel()

    lpm_flat = ln_perm_mean.ravel()

    lps_flat = ln_perm_sd.ravel()



    n_cells = tm_flat.size



    # Calculate indices for each quantile

    indices = np.empty(n_pvalues, dtype=np.int64)

    for p_idx in range(n_pvalues):

        indices[p_idx] = int(quantiles[p_idx] * n_samples)



    # Initialize output arrays with p_value dimension

    thickness_p = np.empty((n_cells, n_pvalues))

    permeability_p = np.empty((n_cells, n_pvalues))

    transmissivity_p = np.empty((n_cells, n_pvalues))



    for i in prange(n_cells):

        tm = tm_flat[i]

        ts = ts_flat[i]

        lpm = lpm_flat[i]

        lps = lps_flat[i]



        if np.isnan(tm) or np.isnan(ts) or np.isnan(lpm) or np.isnan(lps):

            for p_idx in range(n_pvalues):

                thickness_p[i, p_idx] = np.nan

                permeability_p[i, p_idx] = np.nan

                transmissivity_p[i, p_idx] = np.nan

            continue



        # Generate samples for this cell (same samples used for all p_values)

        thickness_samples = np.empty(n_samples)

        ln_perm_samples = np.empty(n_samples)



        for j in range(n_samples):

            if ts == 0:

                thickness_samples[j] = tm

            else:

                thickness_samples[j] = max(0.01, np.random.normal(tm, ts))

            ln_perm_samples[j] = np.random.normal(lpm, lps)



        # Transmissivity samples

        tr_samples = np.exp(ln_perm_samples) * thickness_samples

        tr_samples.sort()



        # Get transmissivity for each quantile

        for p_idx in range(n_pvalues):

            idx = indices[p_idx]

            transmissivity_p[i, p_idx] = tr_samples[idx]



        # Analytical percentiles for all p_values

        for p_idx in range(n_pvalues):

            z_score = z_scores[p_idx]

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

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



    # Reshape to include p_value dimension

    output_shape = (*shape, n_pvalues)

    return (

        thickness_p.reshape(output_shape),

        permeability_p.reshape(output_shape),

        transmissivity_p.reshape(output_shape),

    )
 
 No newline at end of file

src/pythermogis/utils/__init__.py

0 → 100644
+0 −0

Empty file added.

src/pythermogis/utils/timer.py

0 → 100644
+28 −0
Changes for src/pythermogis/utils/timer.py: 28 added lines, 0 removed lines.
Original line number Diff line number Diff line 
import sys

import timeit





def print_time(start, message: str = "", verbose: bool = True):

    if not verbose:

        return



    time_taken_seconds = timeit.default_timer() - start

    time_taken_message = format_time(time_taken_seconds)

    print(f"{message}{time_taken_message}", flush=True, file=sys.stdout)



    return timeit.default_timer()





def format_time(seconds: float) -> str:

    hours = seconds // 3600

    minutes = (seconds % 3600) // 60

    seconds = seconds % 60



    if hours > 0:

        return f"{hours:.0f}h {minutes:.0f}m {seconds:.0f}s"

    elif minutes > 0:

        return f"{minutes:.0f}m {seconds:.0f}s"

    elif seconds > 10:

        return f"{seconds:.1f}s"

    else:

        return f"{seconds:.3f}s"

src/pythermogis/workflow/overview/__init__.py

0 → 100644
+0 −0

Empty file added.