Ensuring that an independent doublet processes grids the same as the Java-code run from utc (b2cbdec9) · Commits · AGS_public / pythermogis

src/thermogis_classes/doublet.py

+37 −2

Original line number	Diff line number	Diff line
		@@ -3,7 +3,19 @@ import numpy as np

		def calculate_performance(doublet, input_params: dict, hydrocarbons: float, depth: float, thickness: float, porosity: float, ntg: float, temperature: float, transmissivity: float, transmissivity_with_ntg: float) -> dict:
		if hydrocarbons == 0.0:
		output_values = {"flow_rate": 0.0, "cop": 0.0, "power": 0.0}
		output_values = {"power": 0.0,
		"heat_pump_power": 0.0,
		"capex": 0.0,
		"opex": 0.0,
		"utc": 0.0,
		"npv": 0.0,
		"hprod": 0.0,
		"cop": 0.0,
		"cophp": 0.0,
		"pres": 0.0,
		"flow_rate": 0.0,
		"welld": 0.0,
		}
		return output_values

		set_doublet_parameters(doublet, transmissivity_with_ntg, depth, porosity, ntg, temperature, input_params["use_stimulation"], input_params["stimKhMax"], input_params["surface_temperature"],
		@@ -11,14 +23,36 @@ def calculate_performance(doublet, input_params: dict, hydrocarbons: float, dept

		doublet.calculateDoubletPerformance(-99999.0, thickness, transmissivity)

		# 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 = {"flow_rate": doublet.doubletCalc1DData.getFlowrate(), "cop": doublet.doubletCalc1DData.getCop(), "power": doublet.doubletCalc1DData.getHpP()}
		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
		doublet.setProjectVariables(False, 0.0)

		return output_values


		def set_doublet_parameters(doublet, transmissivity_with_ntg, depth, porosity, ntg, temperature, useStimulation, stimKhMax, surface_temperature, return_temperature, use_heat_pump, max_cooling_temperature_range,
		hp_minimum_injection_temperature):
		if not useStimulation or transmissivity_with_ntg > stimKhMax:
		@@ -39,6 +73,7 @@ def set_doublet_parameters(doublet, transmissivity_with_ntg, depth, porosity, nt
		doublet.doubletCalc1DData.setInjectionTemp(injectionTemp)
		doublet.doubletCalc1DData.setDhReturnTemp(return_temperature)


		def get_doublet_output(doublet):
		flowr = doublet.doubletCalc1DData.getFlowrate()
		cop = doublet.doubletCalc1DData.getCop()

tests/tests/test_doublet.py

+43 −27

Original line number	Diff line number	Diff line
		import copy
		import shutil
		import subprocess
		from os import path
		from pathlib import Path
		from unittest import TestCase

		import numpy as np

		from src.thermogis_classes.doublet import set_doublet_parameters, get_doublet_output, calculate_performance
		from src.utils.convert_python_to_java import convert_list_to_ArrayList
		from src.utils.get_repo_root import get_repo_root
		from pygridsio.pygridsio import read_grid
		import xarray as xr

		# Import java classes
		import jpype
		from jpype import JClass
		import jpype.imports
		import numpy as np
		from jpype.types import *
		from pygridsio.pygridsio import read_grid

		from src.thermogis_classes.doublet import set_doublet_parameters, calculate_performance
		from src.utils.get_repo_root import get_repo_root


		class PyThermoGIS(TestCase):
		@@ -100,12 +97,6 @@ class PyThermoGIS(TestCase):

		# Instantiate the UTC properties class
		propsBuilder = UTCPropertiesBuilder()
		propsBuilder.setInputDataDir(str(self.model_path) + "/InputData")
		propsBuilder.setResultsDir(str(self.model_path) + "/Results")
		propsBuilder.setCheckCopiedFiles(False)
		propsBuilder.setAquifers(convert_list_to_ArrayList(ArrayList(), ["simplified", "simplified_stacked"]))
		propsBuilder.setLevelOfDetail(1)
		propsBuilder.setPValues([50.0])
		utc_properties = propsBuilder.build()

		# Create an instance of a ThermoGISDoublet
		@@ -114,7 +105,7 @@ class PyThermoGIS(TestCase):
		# Read grids
		thickness_grid = read_grid(self.model_path / "InputData" / "simplified__thick.zmap")
		ntg_grid = read_grid(self.model_path / "InputData" / "simplified__ntg.zmap")
		porosity_grid = read_grid(self.model_path / "InputData" / "simplified__phi.zmap")
		porosity_grid = read_grid(self.model_path / "InputData" / "simplified__phi.zmap") / 100
		depth_grid = read_grid(self.model_path / "InputData" / "simplified__top.zmap")
		hc_accum_grid = read_grid(self.model_path / "InputData" / "simplified__hc_accum.zmap")

		@@ -129,14 +120,23 @@ class PyThermoGIS(TestCase):

		# Model Parameters from a config
		input_params = {"hp_minimum_injection_temperature": 15, "return_temperature": 30, "surface_temperature": 10, "max_cooling_temperature_range": 100, "stimKhMax": 20, "use_stimulation": False,
		"use_heat_pump": False}
		"use_heat_pump": False, "utc_cutoff_shallow": 5.1, "utc_cutoff_deep": 6.5, "utc_cutoff_depth": 4000.0}

		# Initiate grids for flowr, cop and power
		# Initiate output grids
		dummy_grid = copy.deepcopy(thickness_grid)
		dummy_grid.data[:, :] = np.nan
		flowr_grid = copy.deepcopy(dummy_grid)
		cop_grid = copy.deepcopy(dummy_grid)
		power_grid = copy.deepcopy(dummy_grid)
		heat_pump_power_grid = copy.deepcopy(dummy_grid)
		capex_grid = copy.deepcopy(dummy_grid)
		opex_grid = copy.deepcopy(dummy_grid)
		utc_grid = copy.deepcopy(dummy_grid)
		npv_grid = copy.deepcopy(dummy_grid)
		hprod_grid = copy.deepcopy(dummy_grid)
		cop_grid = copy.deepcopy(dummy_grid)
		cophp_grid = copy.deepcopy(dummy_grid)
		pres_grid = copy.deepcopy(dummy_grid)
		flow_rate_grid = copy.deepcopy(dummy_grid)
		welld_grid = copy.deepcopy(dummy_grid)

		# Set Parameters on doublet
		for i in range(len(thickness_grid.x)):
		@@ -149,16 +149,32 @@ class PyThermoGIS(TestCase):
		thickness = thickness_grid.isel(x=i, y=j).data
		transmissivity = transmissivity_grid.isel(x=i, y=j).data
		hydrocarbons = hc_accum_grid.isel(x=i, y=j).data

		# Run doublet simulation
		output_values = calculate_performance(doublet, input_params, hydrocarbons, depth, thickness, porosity, ntg, temperature, transmissivity, transmissivity_with_ntg)

		flowr_grid.values[j, i] = output_values["flow_rate"]
		cop_grid.values[j, i] = output_values["cop"]
		# Assign output values to grids
		power_grid.values[j, i] = output_values["power"]

		print(flowr_grid)

		benchmark_grid = read_grid(self.test_benchmark_output / "simplified__flowr_P50.nc")
		print(benchmark_grid)
		heat_pump_power_grid.values[j, i] = output_values["heat_pump_power"]
		capex_grid.values[j, i] = output_values["capex"]
		opex_grid.values[j, i] = output_values["opex"]
		utc_grid.values[j, i] = output_values["utc"]
		npv_grid.values[j, i] = output_values["npv"]
		hprod_grid.values[j, i] = output_values["hprod"]
		cop_grid.values[j, i] = output_values["cop"]
		cophp_grid.values[j, i] = output_values["cophp"]
		pres_grid.values[j, i] = output_values["pres"]
		flow_rate_grid.values[j, i] = output_values["flow_rate"]
		welld_grid.values[j, i] = output_values["welld"]

		xr.testing.assert_allclose(power_grid, read_grid(self.test_benchmark_output / "simplified__power_P50.nc"), atol=0.1)
		xr.testing.assert_allclose(utc_grid, read_grid(self.test_benchmark_output / "simplified__utc_P50.nc"), atol=0.1)
		xr.testing.assert_allclose(npv_grid, read_grid(self.test_benchmark_output / "simplified__npv_P50.nc"), atol=0.1)
		xr.testing.assert_allclose(hprod_grid, read_grid(self.test_benchmark_output / "simplified__hprod_P50.nc"), atol=10000)
		xr.testing.assert_allclose(cop_grid, read_grid(self.test_benchmark_output / "simplified__cop_P50.nc"), atol=0.1)
		xr.testing.assert_allclose(pres_grid, read_grid(self.test_benchmark_output / "simplified__pres_P50.nc"), atol=0.1)
		xr.testing.assert_allclose(flow_rate_grid, read_grid(self.test_benchmark_output / "simplified__flowr_P50.nc"), atol=2)
		xr.testing.assert_allclose(welld_grid, read_grid(self.test_benchmark_output / "simplified__welld_P50.nc"), atol=0.1)

		def checkOutputGridValues(self):
		flowRateGrid = read_grid(self.model_path / "Results" / "simplified" / "BaseCase" / "simplified__flowr_P50.nc")