TNO Intern

Resolve "Optimize the UTC scropts"

Hen Brettrequested to merge
87-optimize-the-utc-scropts into python-rewrite-of-java-code

Closes #87

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.utc_properties 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/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/utc/cooling_temp.py

+1 −0
Changes for src/pythermogis/workflow/utc/cooling_temp.py: 1 added line, 0 removed lines.
Original line number Diff line number Diff line 
from pythermogis.workflow.utc.flow import calculate_volumetric_flow

from pythermogis.workflow.utc.doublet_utils import calc_lifetime





def calculate_cooling_temperature(

src/pythermogis/workflow/utc/doublet.py

+32 −20
Changes for src/pythermogis/workflow/utc/doublet.py: 32 added lines, 20 removed lines.
Original line number Diff line number Diff line 
import timeit

from dataclasses import dataclass



from pythermogis.workflow.utc.utc_properties import UTCConfiguration

from pythermogis.workflow.utc.doublet_utils import calculate_injection_temp_with_heat_pump, \

    calc_lifetime

from pythermogis.workflow.utc.pressure import calculate_max_pressure, optimize_pressure

import numpy as np



from pythermogis.utils.timer import print_time

from pythermogis.workflow.utc.cooling_temp import calculate_cooling_temperature

from pythermogis.workflow.utc.doublet_utils import (

    calc_lifetime,

    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
@@ -47,31 +53,31 @@ class DoubletOutput: 
    production_temp: float

    injection_temp: float



def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput) -> DoubletOutput | None:

def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput, verbose: bool = False) -> DoubletOutput | None:

    timer = timeit.default_timer()

    # determine initial well distance

    well_distance = (

        (props.optim_dist_well_dist_min + props.optim_dist_well_dist_max) / 2

        np.mean([props.optim_dist_well_dist_min + props.optim_dist_well_dist_max])

        if props.optim_well_dist

        else props.default_well_distance

    )



    injection_temperature = (

        max(input.temperature - props.max_cooling_temp_range,

            props.hp_minimum_injection_temperature)

        if props.use_heat_pump

        else max(input.temperature - props.max_cooling_temp_range,

                 props.dh_return_temp)

    )



    if props.use_heat_pump and props.calculate_cop and not props.hp_application_mode:

    # determine injection temperature

    if props.use_heat_pump:

        injection_temperature = calculate_injection_temp_with_heat_pump(

            input.temperature,

            props.hp_direct_heat_input_temp,

            props.dh_return_temp,

            input.temperature,

            props.max_cooling_temp_range,

            props.hp_minimum_injection_temperature

            props.hp_minimum_injection_temperature,

        )

    else:

        injection_temperature = max(input.temperature - props.max_cooling_temp_range,

            props.dh_return_temp)

    timer = print_time(timer, "injection temperature: ", verbose=verbose)



    # calculate maximum pressure

    drawdown_pressure = calculate_max_pressure(

        props,

        input,
@@ -79,10 +85,11 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput) 
        well_distance,

        injection_temperature,

    )



    if drawdown_pressure == 0:

        return None

    timer = print_time(timer, "maximum pressure: ", verbose=verbose)



    # cooling temperature and well distance optimization

    if props.optim_well_dist:

        end_temperature_p = (

                props.optim_dist_cooling_fraction *
@@ -96,6 +103,7 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput) 
            well_distance,

            injection_temperature,

        )

    timer = print_time(timer, "cooling temperature: ", verbose=verbose)



    if props.optim_well_dist:

        well_distance = optimize_well_distance(
@@ -104,13 +112,16 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput) 
            drawdown_pressure,

            injection_temperature,

        )

    timer = print_time(timer, "well distance optimizer: ", verbose=verbose)



    # stimulation capex

    stimulation_capex = (

        0.0

        if (not props.use_stimulation or input.transmissivity_with_ntg > props.stim_kh_max)

        else props.stimulation_capex

    )



    # pressure optimizer

    pressure_results = optimize_pressure(

        props=props,

        input=input,
@@ -119,10 +130,11 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput) 
        injection_temp=injection_temperature,

        stimulation_capex=stimulation_capex,

    )



    timer = print_time(timer, "pressure optimizer: ", verbose=verbose)

    if pressure_results is None:

        return None



    # everything underneath here is fast, no point in optimizing

    heat_power_per_doublet = pressure_results.heat_power_per_doublet

    flowrate = pressure_results.flowrate

    discounted_heat_produced_p = pressure_results.discounted_heat_produced
@@ -179,7 +191,7 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput) 
            * (1 - props.tax_rate)

    )

    opex_first_prod_year = total_opex_ts[props.drilling_time]

    hp_cop = 3.0

    print_time(timer, "economics: ", verbose=verbose)



    return DoubletOutput(

        power=heat_power_per_doublet,
@@ -192,7 +204,7 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput) 
        npv=npv,

        hprod=discounted_heat_produced_p,

        cop=cop,

        cophp=hp_cop,

        cophp=pressure_results.heat_pump_cop,

        pres=drawdown_pressure / 1e5,

        flow=flowrate,

        welld=well_distance,