TNO Intern

Commit 785b4096 authored by Hen Brett's avatar Hen Brett 🐔
Browse files

Merge branch '11-run-speed-and-efficiency-tests' into 'main' 

Resolve "Run speed and efficiency tests"

Closes #11

See merge request AGS/pythermogis!15
parents c6a04dba 5d47f2d8

pythermogis/__init__.py

+1 −1
Original line number Diff line number Diff line 
from thermogis_classes.doublet import calculate_performance_of_single_location, calculate_doublet_performance
 
 No newline at end of file
 
from .thermogis_classes.doublet import calculate_performance_of_single_location, calculate_doublet_performance
 
 No newline at end of file

pythermogis/statistics/calculate_thick_perm_trans.py

+54 −8
Original line number Diff line number Diff line 
import numpy as np

from scipy import stats

import xarray as xr





def generate_thickness_permeability_transmissivity_for_pvalue(thickness_mean: float, thickness_sd: float, ln_permeability_mean: float, ln_permeability_sd: float, Pvalue: float, nSamples: int = 10000) -> float:

def generate_thickness_permeability_transmissivity_for_pvalue(thickness_mean: float, thickness_sd: float, ln_permeability_mean: float, ln_permeability_sd: float, p_values: float, nSamples: int = 10000) -> float:

    """

    Given thickness provided as a normal distribution and ln(permeability) provided as a normal distribution, and a specific p-value, then generate the values for that p-value of the

    -thickness
@@ -20,16 +20,62 @@ def generate_thickness_permeability_transmissivity_for_pvalue(thickness_mean: fl 
    :return:

    thickness, permeability, transmissivity

    """

    p_value_fraction = Pvalue / 100

    if np.isnan(thickness_mean) | np.isnan(thickness_sd) | np.isnan(ln_permeability_mean) | np.isnan(ln_permeability_sd):

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



    p_value_fractions = p_values / 100



    thickness_dist = stats.norm(loc=thickness_mean, scale=thickness_sd)

    thickness_pvalue = thickness_dist.ppf(1 - p_value_fraction)

    thickness_pvalues = thickness_dist.ppf(1 - p_value_fractions)



    ln_permeability_dist = stats.norm(loc=ln_permeability_mean, scale=ln_permeability_sd)

    permeability_pvalue = np.exp(ln_permeability_dist.ppf(1 - p_value_fraction))

    permeability_pvalues = np.exp(ln_permeability_dist.ppf(1 - p_value_fractions))



    # Sampling method for transmissivity

    transmissivity_samples = np.sort(np.exp(ln_permeability_dist.rvs(nSamples) + np.log(thickness_dist.rvs(nSamples))))

    transmissivity_pvalue_sampled = transmissivity_samples[int((1 - p_value_fraction) * nSamples)]

    thickness_samples = thickness_dist.rvs(nSamples)

    thickness_samples = np.where(thickness_samples < 0.0, 0.1, thickness_samples)   # Ensure no negative thicknesses

    transmissivity_samples = np.sort(np.exp(ln_permeability_dist.rvs(nSamples) + np.log(thickness_samples)))



    transmissivity_pvalues_sampled = transmissivity_samples[int((1 - p_value_fractions) * nSamples)]



    return thickness_pvalues, permeability_pvalues, transmissivity_pvalues_sampled



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:

    """

    Given thickness provided as a normal distribution and ln(permeability) provided as a normal distribution, and a specific p-value, then generate the values for that p-value of the

    -thickness

    -permeability

    -transmissivity



    Note: Transmissivity is a function of ln(permeability) * thickness, and so no analytical solution exists to combine these two probability distributions and so the transmissivity distribution has to be sampled.



    :param thickness_mean:

    :param thickness_sd:

    :param ln_permeability_mean:

    :param ln_permeability_sd:

    :param p_values:

    :param nSamples:

    :return:

    thickness, permeability, transmissivity

    """

    if np.isnan(thickness_mean) | np.isnan(thickness_sd) | np.isnan(ln_permeability_mean) | np.isnan(ln_permeability_sd):

        return np.full(len(p_values), np.nan), np.full(len(p_values), np.nan), np.full(len(p_values), np.nan)



    p_value_fractions = 1 - p_values / 100



    thickness_dist = stats.norm(loc=thickness_mean, scale=thickness_sd)

    thickness_pvalues = thickness_dist.ppf(p_value_fractions)



    ln_permeability_dist = stats.norm(loc=ln_permeability_mean, scale=ln_permeability_sd)

    permeability_pvalues = np.exp(ln_permeability_dist.ppf(p_value_fractions))



    # Sampling method for transmissivity

    thickness_samples = thickness_dist.rvs(nSamples)

    thickness_samples = np.where(thickness_samples < 0.0, 0.1, thickness_samples)   # Ensure no negative thicknesses

    transmissivity_samples = np.sort(np.exp(ln_permeability_dist.rvs(nSamples) + np.log(thickness_samples)))



    sample_indexes = p_value_fractions * nSamples

    sample_indexes = sample_indexes.astype(int)

    transmissivity_pvalues_sampled = transmissivity_samples[sample_indexes]



    return thickness_pvalue, permeability_pvalue, transmissivity_pvalue_sampled
 
    return thickness_pvalues, permeability_pvalues, transmissivity_pvalues_sampled
 
 No newline at end of file

pythermogis/thermogis_classes/doublet.py

+77 −53
Original line number Diff line number Diff line 
from typing import List



import timeit

import numpy as np

import xarray as xr

from jpype import JClass



from pythermogis.physics.temperature_grid_calculation import calculate_temperature_from_gradient

from pythermogis.statistics.calculate_thick_perm_trans import generate_thickness_permeability_transmissivity_for_pvalue

from pythermogis.statistics.calculate_thick_perm_trans import generate_thickness_permeability_transmissivity_for_pvalue, \

    generate_thickness_permeability_transmissivity_for_pvalues

from pythermogis.thermogis_classes.jvm_start import start_jvm
@@ -13,7 +15,27 @@ def calculate_doublet_performance(input_data: xr.Dataset, 
                                  input_params: dict = None,

                                  p_values: List[float] = [50.0]) -> xr.Dataset:

    """

    Perform a ThermoGIS Doublet performance assessment. This will occur across all dimensions of the input_data (ie. input data can have a single value for each required variable, or it can be 1Dimensional or a 2Dimensional grid)

    Perform a ThermoGIS Doublet performance simulation. This will occur across all dimensions of the input_data (ie. input data can have a single value for each required variable, or it can be 1Dimensional or a 2Dimensional grid)



    The default doublet simulation settings are for the ThermoGIS BaseCase; but major input parameters can be set using the input_params variable. This is a dictonary which by default contain the following keys:

    input_params = {"hp_minimum_injection_temperature": 15,

         "return_temperature": 30,

         "surface_temperature": 10,

         "degrees_per_km": 31,

         "max_cooling_temperature_range": 100,

         "stimKhMax": 20,

         "use_stimulation": False,

         "use_heat_pump": False,

         "calculate_cop": True,

         "hp_application_mode": False,

         "hp_direct_heat_input_temp": 70.0,

         "utc_cutoff_shallow": 5.1,

         "utc_cutoff_deep": 6.5,

         "utc_cutoff_depth": 4000.0,

         "rng_seed": np.random.randint(low=0, high=10000)}

    To change only one of these input parameters you can provide a dictionary with a single key to this method, and it will use the default values for the rest of the model parameters.

    e.g. input_params = {"use_stimulation":True}





    :param input_data:

        A xr.Dataset (input_data) with the required variables: "thickness_mean", "thickness_sd", "porosity", "ntg", "depth", "ln_permeability_mean", "ln_permeability_sd", Performance will be calculated across all dimensions.
@@ -29,8 +51,6 @@ def calculate_doublet_performance(input_data: xr.Dataset, 
        "power", "heat_pump_power", "capex", "opex", "utc", "npv", "hprod", "cop", "cophp", "pres", "flow_rate", "welld"



        The output will have the same dimensions as the input_data class, with the additional p_value dimension





    """



    # Check that all essential variables are provided
@@ -61,14 +81,14 @@ def calculate_doublet_performance(input_data: xr.Dataset, 
    output_data["temperature"] = input_data["temperature"].copy()



    # Calculate Thickness, Permeability and Transmissivity for each P-value

    output_data["thickness"], output_data["permeability"], output_data["transmissivity"] = xr.apply_ufunc(generate_thickness_permeability_transmissivity_for_pvalue,

    output_data["thickness"], output_data["permeability"], output_data["transmissivity"] = xr.apply_ufunc(generate_thickness_permeability_transmissivity_for_pvalues,

                                                                 input_data.thickness_mean,

                                                                 input_data.thickness_sd,

                                                                 input_data.ln_permeability_mean,

                                                                 input_data.ln_permeability_sd,

                                                                 p_values,

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

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

                                                                 input_core_dims=[[], [], [], [], ["p_value"]],

                                                                 output_core_dims=[["p_value"], ["p_value"], ["p_value"]],

                                                                 vectorize=True

                                                             )
@@ -92,7 +112,7 @@ def calculate_doublet_performance(input_data: xr.Dataset, 
                                        kwargs={"doublet": doublet, "input_params": input_params},

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

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

                                        vectorize=True,

                                        vectorize=True

                                        )



    # Assign output DataArrays to the output_data object
@@ -154,38 +174,6 @@ def apply_input_params(input_params: dict) -> dict: 
        return {key: input_params.get(key, input_params_basecase[key]) for key in input_params_basecase}





def instantiate_thermogis_doublet(input_params: dict):

    """

    Instantiate a ThermoGIS Doublet class, with a set random seed if provided

    :param rng_seed:

    :return:

        doublet, a JPype instantiated object of the ThermoGISDoublet class

    """

    # Start the jvm

    start_jvm()

    # Instantiate doublet class

    Logger = JClass("logging.Logger")

    Mockito = JClass("org.mockito.Mockito")

    RNG = JClass("tno.geoenergy.stochastic.RandomNumberGenerator")

    ThermoGISDoublet = JClass("thermogis.calc.doublet.ThermoGisDoublet")

    UTCPropertiesBuilder = JClass("thermogis.properties.builders.UTCPropertiesBuilder")



    # Instantiate the UTC properties class

    propsBuilder = UTCPropertiesBuilder()

    utc_properties = (propsBuilder

                      .setUseStimulation(input_params["use_stimulation"])

                      .setUseHeatPump(input_params["use_heat_pump"])

                      .setCalculateCop(input_params["calculate_cop"])

                      .setHpApplicationMode(input_params["hp_application_mode"])

                      .setHpDirectHeatInputTemp(input_params["hp_direct_heat_input_temp"])

                      .build())



    # Instantiate random number generator:

    rng = RNG(input_params["rng_seed"])



    # Create an instance of a ThermoGISDoublet

    doublet = ThermoGISDoublet(Mockito.mock(Logger), rng, utc_properties)

    return doublet





def calculate_performance_of_single_location(mask: float, depth: float, thickness: float, porosity: float, ntg: float, temperature: float, transmissivity: float, transmissivity_with_ntg: float, doublet=None,
@@ -261,6 +249,46 @@ def calculate_performance_of_single_location(mask: float, depth: float, thicknes 
        "cophp"], output_values["pres"], output_values["flow_rate"], output_values["welld"]





def instantiate_thermogis_doublet(input_params: dict):

    """

    Instantiate a ThermoGIS Doublet class, with a set random seed if provided

    :param rng_seed:

    :return:

        doublet, a JPype instantiated object of the ThermoGISDoublet class

    """

    # Start the jvm

    start_jvm()

    # Instantiate doublet class

    Logger = JClass("logging.Logger")

    Mockito = JClass("org.mockito.Mockito")

    RNG = JClass("tno.geoenergy.stochastic.RandomNumberGenerator")

    ThermoGISDoublet = JClass("thermogis.calc.doublet.ThermoGisDoublet")

    UTCPropertiesBuilder = JClass("thermogis.properties.builders.UTCPropertiesBuilder")



    # Instantiate the UTC properties class

    propsBuilder = UTCPropertiesBuilder()

    utc_properties = (propsBuilder

                      .setUseStimulation(input_params["use_stimulation"])

                      .setUseHeatPump(input_params["use_heat_pump"])

                      .setCalculateCop(input_params["calculate_cop"])

                      .setHpApplicationMode(input_params["hp_application_mode"])

                      .setHpDirectHeatInputTemp(input_params["hp_direct_heat_input_temp"])

                      .build())



    # Instantiate random number generator:

    rng = RNG(input_params["rng_seed"])



    # Create an instance of a ThermoGISDoublet

    doublet = ThermoGISDoublet(Mockito.mock(Logger), rng, utc_properties)



    # Set parameters that do not change across cells

    doublet.doubletCalc1DData.setSurfaceTemperature(input_params["surface_temperature"])

    doublet.doubletCalc1DData.setUseHeatPump(input_params["use_heat_pump"])

    doublet.doubletCalc1DData.setDhReturnTemp(input_params["return_temperature"])





    return doublet



def set_doublet_parameters(doublet, transmissivity_with_ntg: float, depth: float, porosity: float, ntg: float, temperature: float, input_params: dict):

    """

    For a single location sets the necessary data on the doublet class, to then run a doublet simulation
@@ -285,16 +313,12 @@ def set_doublet_parameters(doublet, transmissivity_with_ntg: float, depth: float 
    doublet.doubletCalc1DData.setDepth(depth)

    doublet.doubletCalc1DData.setPorosity(porosity)

    doublet.doubletCalc1DData.setNtg(ntg)

    doublet.doubletCalc1DData.setSurfaceTemperature(input_params["surface_temperature"])

    doublet.doubletCalc1DData.setReservoirTemp(temperature)

    doublet.doubletCalc1DData.setUseHeatPump(input_params["use_heat_pump"])



    if not input_params["use_heat_pump"]:

        injectionTemp = np.max([temperature - input_params["max_cooling_temperature_range"], input_params["return_temperature"]])

        doublet.doubletCalc1DData.setInjectionTemp(np.max([temperature - input_params["max_cooling_temperature_range"], input_params["return_temperature"]]))

    elif input_params["use_heat_pump"] and input_params["calculate_cop"] and not input_params["hp_application_mode"]:

        injectionTemp = doublet.calculateInjectionTempWithHeatPump(temperature, input_params["hp_direct_heat_input_temp"])

        doublet.doubletCalc1DData.setInjectionTemp(doublet.calculateInjectionTempWithHeatPump(temperature, input_params["hp_direct_heat_input_temp"]))

    else:

        injectionTemp = np.max([temperature - input_params["max_cooling_temperature_range"], input_params["hp_minimum_injection_temperature"]])

        doublet.doubletCalc1DData.setInjectionTemp(np.max([temperature - input_params["max_cooling_temperature_range"], input_params["hp_minimum_injection_temperature"]]))

    doublet.doubletCalc1DData.setInjectionTemp(injectionTemp)

    doublet.doubletCalc1DData.setDhReturnTemp(input_params["return_temperature"])

pythermogis/thermogis_classes/jvm_start.py

+1 −0
Original line number Diff line number Diff line 
from importlib.resources import files



import jpype

import logging





def get_jvm_path() -> str:

tests/resources/test_input/SLDNA/SLDNA__hc_accum.nc

0 → 100644
+33.7 KiB

File added.

No diff preview for this file type.

View file