Adding the benchmark data for extra P-values, for HP, STIM and HP_STIM

from pythermogis.thermogis_classes.java_start import start_jvm

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,

                                             input_params: dict = None) -> float:

    """

        Calculate the performance of a doublet at a single location

    :param mask:

        mask value, if not nan set all output to 0.0

    :param depth:

        depth of the aquifer

    :param thickness:

        thickness of the aquifer

    :param porosity:

        porosity of the aquifer

    :param ntg:

        net-to-gross of the aquifer

    :param temperature:

        temperature of the aquifer

    :param transmissivity:

        transmissivity of the aquifer

    :param transmissivity_with_ntg:

        transmissivity * ntg of the aquifer

    :param doublet:

        a ThermoGIS doublet class

    :param input_params:

        a dictionary containing extra parameters for running the doublet simulation

    :return:

        the values "power", "heat_pump_power", "capex", "opex", "utc", "npv", "hprod", "cop", "cophp", "pres", "flow_rate", "welld" from the doublet calculation

    """

    if np.isnan(thickness):

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

    if not np.isnan(mask):

        return 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

    set_doublet_parameters(doublet, transmissivity_with_ntg, depth, porosity, ntg, temperature, input_params)

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

    doublet.calculateDoubletPerformance(-9999.0, thickness, transmissivity)

    if doublet.getUtcPeurctkWh() == -9999.0:  # If calculation was not successful, return all 0.0

        return 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

    # calculate net-present-value using the utc-cutoffs

    if depth > input_params["utc_cutoff_depth"]:

        utc_cut = input_params["utc_cutoff_deep"]

    else:

        utc_cut = input_params["utc_cutoff_shallow"]

    hprod = doublet.economicalData.getDiscountedHeatProducedP()

    npv = 1e-6 * (utc_cut - doublet.getUtcPeurctkWh()) * 3.6 * hprod * (1 - doublet.economicalData.getTaxRate())

    # get values from doublet

    output_values = {"power": doublet.doubletCalc1DData.getHpP(),

                     "heat_pump_power": doublet.doubletCalc1DData.getHpPHeatPump(),

                     "capex": doublet.economicalData.getSumcapex(),

                     "opex": doublet.economicalData.getOpexFirstProdYear(),

                     "utc": doublet.getUtcPeurctkWh(),

                     "npv": npv,

                     "hprod": hprod,

                     "cop": doublet.doubletCalc1DData.getCop(),

                     "cophp": doublet.doubletCalc1DData.getCopHpP(),

                     "pres": doublet.doubletCalc1DData.getPresP() / 1e5,

                     "flow_rate": doublet.doubletCalc1DData.getFlowrate(),

                     "welld": doublet.doubletCalc1DData.getWellDistP(),

                     }

    # Reset doublet variables for next calculation

    doublet.setProjectVariables(False, 0.0)

    return output_values["power"], output_values["heat_pump_power"], output_values["capex"], output_values["opex"], output_values["utc"], output_values["npv"], output_values["hprod"], output_values["cop"], output_values[

        "cophp"], output_values["pres"], output_values["flow_rate"], output_values["welld"]

def calculate_doublet_performance_across_dimensions(input_data: xr.Dataset,

                                                    input_params: dict = None,

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

        p_values = [50.0]

    # Apply provided input_params; or return the base case

    input_params = implement_input_params(input_params)

    input_params = apply_input_params(input_params)

    # Generate temperature values from gradient if no temperature provided

    if "temperature" not in input_data:

    output_data = output_data.expand_dims({"p_value": p_values})

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

    # TODO: Find a less ugly way to handle all these outputs than instantiating different arrays; probably by setting the output_core_dims in the apply_ufunc

    thickness_data = []

    permeability_data = []

    transmissivity_data = []

    pres_data = []

    flow_rate_data = []

    welld_data = []

    # TODO: Find a less ugly way to handle all these outputs than instantiating different arrays; probably by setting the output_core_dims in the apply_ufunc

    for i, p in enumerate(p_values):

        output_data_arrays = xr.apply_ufunc(calculate_performance_of_single_location,

                                            input_data.mask,

    return output_data

def validate_input_data(input_data):

def validate_input_data(input_data: xr.Dataset):

    """

    Ensure that the input_data Dataset contains the minimum required variables

    :param input_data:

    :return:

    """

    missing_variables = []

    for variable in ["thickness_mean", "thickness_sd", "porosity", "ntg", "depth", "ln_permeability_mean", "ln_permeability_sd"]:

        if variable not in input_data:

        raise ValueError(f"provided input Dataset does not contain the following required variables: {missing_variables}")

def implement_input_params(input_params):

    # Input parameters

def apply_input_params(input_params: dict) -> dict:

    """

    if input_params is None, return the basecase input_params, if input_params contains keys, then change these keys in the basecase input params to the values provided by the user

    :param input_params:

    :return:

    """

    input_params_basecase = {"hp_minimum_injection_temperature": 15,

                             "return_temperature": 30,

                             "surface_temperature": 10,

                             "utc_cutoff_deep": 6.5,

                             "utc_cutoff_depth": 4000.0,

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

    if input_params is None:  # If no input_params provided, use the basecase

        input_params = input_params_basecase

    else:

        input_params = {key: input_params.get(key, input_params_basecase[key]) for key in input_params_basecase}

    return input_params

    if input_params is None:  # If no input_params provided, return the basecase

        return input_params_basecase

    else:  # If input_params provided, then go through the keys in the basecase dictionary, and any keys in common can be changed to the user provided input_params

        return {key: input_params.get(key, input_params_basecase[key]) for key in input_params_basecase}

def instantiate_thermogis_doublet(input_params):

    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,

                                             input_params: dict = None) -> float:

    """

        Calculate the performance of a doublet at a single location

    :param mask:

        mask value, if not nan set all output to 0.0

    :param depth:

        depth of the aquifer

    :param thickness:

        thickness of the aquifer

    :param porosity:

        porosity of the aquifer

    :param ntg:

        net-to-gross of the aquifer

    :param temperature:

        temperature of the aquifer

    :param transmissivity:

        transmissivity of the aquifer

    :param transmissivity_with_ntg:

        transmissivity * ntg of the aquifer

    :param doublet:

        a ThermoGIS doublet class

    :param input_params:

        a dictionary containing extra parameters for running the doublet simulation

    :return:

        the values "power", "heat_pump_power", "capex", "opex", "utc", "npv", "hprod", "cop", "cophp", "pres", "flow_rate", "welld" from the doublet calculation

    """

    if np.isnan(thickness):

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

    if not np.isnan(mask):

        return 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

    set_doublet_parameters(doublet, transmissivity_with_ntg, depth, porosity, ntg, temperature, input_params)

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

    doublet.calculateDoubletPerformance(-9999.0, thickness, transmissivity)

    if doublet.getUtcPeurctkWh() == -9999.0:  # If calculation was not successful, return all 0.0

        return 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

    # calculate net-present-value using the utc-cutoffs

    if depth > input_params["utc_cutoff_depth"]:

        utc_cut = input_params["utc_cutoff_deep"]

    else:

        utc_cut = input_params["utc_cutoff_shallow"]

    hprod = doublet.economicalData.getDiscountedHeatProducedP()

    npv = 1e-6 * (utc_cut - doublet.getUtcPeurctkWh()) * 3.6 * hprod * (1 - doublet.economicalData.getTaxRate())

    # get values from doublet

    output_values = {"power": doublet.doubletCalc1DData.getHpP(),

                     "heat_pump_power": doublet.doubletCalc1DData.getHpPHeatPump(),

                     "capex": doublet.economicalData.getSumcapex(),

                     "opex": doublet.economicalData.getOpexFirstProdYear(),

                     "utc": doublet.getUtcPeurctkWh(),

                     "npv": npv,

                     "hprod": hprod,

                     "cop": doublet.doubletCalc1DData.getCop(),

                     "cophp": doublet.doubletCalc1DData.getCopHpP(),

                     "pres": doublet.doubletCalc1DData.getPresP() / 1e5,

                     "flow_rate": doublet.doubletCalc1DData.getFlowrate(),

                     "welld": doublet.doubletCalc1DData.getWellDistP(),

                     }

    # Reset doublet variables for next calculation

    doublet.setProjectVariables(False, 0.0)

    return output_values["power"], output_values["heat_pump_power"], output_values["capex"], output_values["opex"], output_values["utc"], output_values["npv"], output_values["hprod"], output_values["cop"], output_values[

        "cophp"], output_values["pres"], output_values["flow_rate"], output_values["welld"]

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

    """

    For a single location sets the necessary data on the doublet class, to then run a doublet simulation

            self.assert_output_and_benchmark_is_close(self.test_benchmark_output_STIM, output_grids_p_value, p_value, "_STIM")

    def test_doublet_grid_HP_STIM(self):

        # This tests that the python API produces (approximately) the same output as the Java code, when running on the same set of input data for the stimulation scenario

        # This tests that the python API produces (approximately) the same output as the Java code, when running on the same set of input data for the heat pump and stimulation scenario

        # Read Input grids

        input_grids = self.read_input_grids()

        p_values = [10, 50, 90]