TNO Intern

Commit 7afe7194 authored by Hen Brett's avatar Hen Brett 🐔
Browse files

Trying out tests of the simulations using the python backend

parent 0ad030a6

src/pythermogis/thermogis_classes/doublet.py

+73 −32
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 workflow.utc.utc_properties import UTCConfiguration





def simulate_doublet(

    output_data: xr.Dataset,
@@ -102,11 +106,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 +121,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=-999.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 +172,24 @@ 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())

    hprod = results.hprod

    npv = 1e-6 * (utc_cut - results.utc) * 3.6 * 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(),

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

+1 −1
Original line number Diff line number Diff line
@@ -112,7 +112,7 @@ class UTCConfiguration: 
    default_well_distance: float = 1500.0

    pump_efficiency: float = 0.6

    pump_depth: float = 300.0

    segment_length: float = 1.0

    segment_length: float = 50.0

    outer_diameter: float = 8.5

    inner_diameter: float = 8.5

    roughness: float = 1.38

tests/test_doc_examples.py

+2 −2
Original line number Diff line number Diff line
@@ -280,7 +280,7 @@ def test_example8(): 
    output_data_path.mkdir(parents=True, exist_ok=True)



    # generate simulation samples across desired reservoir properties

    Nsamples = 100

    Nsamples = 1000

    thickness_samples = np.random.uniform(low=150, high=300, size=Nsamples)

    porosity_samples = np.random.uniform(low=0.5, high=0.8, size=Nsamples)

    ntg_samples = np.random.uniform(low=0.25, high=0.5, size=Nsamples)
@@ -298,7 +298,7 @@ def test_example8(): 
    )



    # For every sample, run a doublet simulation store the output values

    simulations = calculate_doublet_performance(reservoir_properties, print_execution_duration=True)

    simulations = calculate_doublet_performance(reservoir_properties, chunk_size=100, print_execution_duration=True)



    # post process the samples to calculate the percentiles of each variable

    percentiles = np.arange(1,99)

tests/test_doublet_speed.py

+3 −4
Original line number Diff line number Diff line
@@ -2,7 +2,7 @@ import shutil 
from os import path

from unittest.case import TestCase



from pygridsio import read_grid, resample_grid

from pygridsio import read_grid, resample_grid, plot_grid



from pythermogis import *
@@ -27,15 +27,14 @@ class PyThermoGIS(TestCase): 
        non_nan_data = data[~np.isnan(data)]

        print(f"Running Performance calculation for ROSL_ROSLU, Pvalues: {p_values}, The number of Non Nan Cells: {len(non_nan_data)}")



        # Run calculation across all dimensions of input_grids, and all provided P_values: The Equivalent calculation with 10 cores takes 43 seconds in the Java code

        # Run calculation across all dimensions of input_grids, and all provided P_values

        start = timeit.default_timer()

        calculate_doublet_performance_stochastic(input_grids, chunk_size=100, rng_seed=123, p_values=p_values)

        time_elapsed = timeit.default_timer() - start

        print(f"Python calculation took {time_elapsed:.1f} seconds.")





    def read_input_grids(self):

        new_cellsize=5000 # in m

        new_cellsize=20e3 # in m

        input_grids = resample_grid(read_grid(self.test_files_out_path / "ROSL_ROSLU__thick.zmap"), new_cellsize=new_cellsize).to_dataset(name="thickness_mean")

        input_grids["thickness_sd"] = resample_grid(read_grid(self.test_files_out_path / "ROSL_ROSLU__thick_sd.zmap"), new_cellsize=new_cellsize)

        input_grids["ntg"] = resample_grid(read_grid(self.test_files_out_path / "ROSL_ROSLU__ntg.zmap"), new_cellsize=new_cellsize)

tests/utc/test_aquifer.py

0 → 100644
+19 −0
Original line number Diff line number Diff line 
import numpy as np

from pygridsio import read_grid,resample_grid

from pathlib import Path



test_files_in_path = Path(__file__).parent.parent / "resources" / "test_input" / "example_data"



def read_input_grids():

    new_cellsize=1000 # in m

    input_grids = resample_grid(read_grid(test_files_in_path / "ROSL_ROSLU__thick.zmap"), new_cellsize=new_cellsize).to_dataset(name="thickness_mean")

    input_grids["thickness_sd"] = resample_grid(read_grid(test_files_in_path / "ROSL_ROSLU__thick_sd.zmap"), new_cellsize=new_cellsize)

    input_grids["ntg"] = resample_grid(read_grid(test_files_in_path / "ROSL_ROSLU__ntg.zmap"), new_cellsize=new_cellsize)

    input_grids["porosity"] = resample_grid(read_grid(test_files_in_path / "ROSL_ROSLU__poro.zmap"), new_cellsize=new_cellsize) / 100

    input_grids["depth"] = resample_grid(read_grid(test_files_in_path / "ROSL_ROSLU__top.zmap"), new_cellsize=new_cellsize)

    input_grids["ln_permeability_mean"] = np.log(resample_grid(read_grid(test_files_in_path / "ROSL_ROSLU__perm.zmap"), new_cellsize=new_cellsize))

    input_grids["ln_permeability_sd"] = resample_grid(read_grid(test_files_in_path / "ROSL_ROSLU__ln_perm_sd.zmap"), new_cellsize=new_cellsize)

    return input_grids



def test_pydoubletcalc_on_aquifer():

    input_grids = read_input_grids()