Merge branch '93-more-workflow-utc-cleanups' into 'python-rewrite-of-java-code' (835de8be) · Commits · AGS_public / pythermogis

src/pythermogis/workflow/utc/constants.py

0 → 100644

+1 −0

Original line number	Diff line number	Diff line
		DARCY_SI = 1.0e-12 / 1.01325
		No newline at end of file

src/pythermogis/workflow/utc/doublet.py

+19 −35

Original line number	Diff line number	Diff line
		import timeit
		from dataclasses import dataclass

		import numpy as np
		from typing import NamedTuple

		from pythermogis.utils.timer import print_time
		from pythermogis.workflow.utc.cooling_temp import calculate_cooling_temperature
		@@ -12,12 +11,13 @@ from pythermogis.workflow.utc.doublet_utils import (
		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
		from pythermogis.workflow.utc.constants import DARCY_SI

		EUR_PER_CT_PER_KWH = 0.36
		NPV_SCALE = 1e-6


		@dataclass
		@dataclass(slots=True, frozen=True)
		class DoubletInput:
		unknown_input_value: float
		thickness: float
		@@ -30,12 +30,10 @@ class DoubletInput:

		@property
		def permeability(self) -> float:
		DARCY_SI = 1.0e-12 / 1.01325
		return self.transmissivity / self.thickness * 1e-3 * DARCY_SI


		@dataclass
		class DoubletOutput:
		class DoubletOutput(NamedTuple):
		power: float
		hppower: float
		capex: float
		@@ -62,7 +60,7 @@ def calculate_doublet_performance(
		timer = timeit.default_timer()
		# determine initial well distance
		well_distance = (
		np.mean([props.optim_dist_well_dist_min + props.optim_dist_well_dist_max])
		(props.optim_dist_well_dist_min + props.optim_dist_well_dist_max) / 2
		if props.optim_well_dist
		else props.default_well_distance
		)
		@@ -142,27 +140,13 @@ def calculate_doublet_performance(
		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
		variable_opex = pressure_results.variable_opex
		fixed_opex = pressure_results.fixed_opex
		total_opex_ts = pressure_results.total_opex_ts
		cop = pressure_results.cop
		drawdown_pressure = pressure_results.drawdown_pressure # overwrite like Java
		sum_capex = pressure_results.sum_capex
		utc = pressure_results.utc
		hp_added_power = pressure_results.hp_added_power
		production_temp = pressure_results.production_temp

		total_variable_opex = sum(variable_opex)
		total_fixed_opex = sum(fixed_opex)
		total_variable_opex = sum(pressure_results.variable_opex)
		total_fixed_opex = sum(pressure_results.fixed_opex)

		utc_eur_ct_per_kwh = (
		input.unknown_input_value
		if utc == input.unknown_input_value
		else utc * EUR_PER_CT_PER_KWH
		if pressure_results.utc == input.unknown_input_value
		else pressure_results.utc * EUR_PER_CT_PER_KWH
		)

		if abs(input.unknown_input_value) < 10:
		@@ -178,7 +162,7 @@ def calculate_doublet_performance(
		input.thickness * input.ntg,
		0.001,
		input.porosity,
		flowrate,
		pressure_results.flowrate,
		input.depth,
		input.temperature,
		props.salinity_surface,
		@@ -197,29 +181,29 @@ def calculate_doublet_performance(
		NPV_SCALE
		* (utc_cutoff - utc_eur_ct_per_kwh)
		* 3.6
		* discounted_heat_produced_p
		* pressure_results.discounted_heat_produced
		* (1 - props.tax_rate)
		)
		opex_first_prod_year = total_opex_ts[props.drilling_time]
		opex_first_prod_year = pressure_results.total_opex_ts[props.drilling_time]
		print_time(timer, "economics: ", verbose=verbose)

		return DoubletOutput(
		power=heat_power_per_doublet,
		hppower=hp_added_power,
		capex=sum_capex,
		power=pressure_results.heat_power_per_doublet,
		hppower=pressure_results.hp_added_power,
		capex=pressure_results.sum_capex,
		var_opex=total_variable_opex,
		fixed_opex=total_fixed_opex,
		opex=opex_first_prod_year,
		utc=utc_eur_ct_per_kwh,
		npv=npv,
		hprod=discounted_heat_produced_p,
		cop=cop,
		hprod=pressure_results.discounted_heat_produced,
		cop=pressure_results.cop,
		cophp=pressure_results.heat_pump_cop,
		pres=drawdown_pressure / 1e5,
		flow=flowrate,
		flow=pressure_results.flowrate,
		welld=well_distance,
		breakthrough=breakthrough_years,
		cooling=end_temperature_p,
		production_temp=production_temp,
		production_temp=pressure_results.production_temp,
		injection_temp=injection_temperature,
		)

src/pythermogis/workflow/utc/doubletcalc.py

+2 −4

Original line number	Diff line number	Diff line
		from __future__ import annotations

		import math
		from dataclasses import dataclass
		from typing import TYPE_CHECKING
		from typing import TYPE_CHECKING, NamedTuple

		from numba import njit
		from pydoubletcalc import Aquifer, Doublet, Well, WellPipeSegment
		@@ -18,8 +17,7 @@ if TYPE_CHECKING:
		INCH_SI = 0.0254


		@dataclass
		class Doublet1DResults:
		class Doublet1DResults(NamedTuple):
		geothermal_powers: float
		cop: float
		flowrate: float

src/pythermogis/workflow/utc/economics.py

+98 −66

Original line number	Diff line number	Diff line
		from __future__ import annotations

		from dataclasses import dataclass
		from typing import TYPE_CHECKING
		from typing import TYPE_CHECKING, NamedTuple

		import numpy as np
		from numba import njit
		from numpy.typing import NDArray

		from pythermogis.workflow.utc.doublet_utils import get_along_hole_length
		@@ -18,8 +18,7 @@ KWH_TO_MJ = 0.36 * 1e9
		HOURS_IN_YEAR = 8760


		@dataclass
		class CapexCalculatorResults:
		class CapexCalculatorResults(NamedTuple):
		sum_capex: float
		total_capex: float
		variable_opex: list[float]
		@@ -28,14 +27,12 @@ class CapexCalculatorResults:
		heat_power_per_year: list[float]


		@dataclass
		class UTCCalculatorResults:
		class UTCCalculatorResults(NamedTuple):
		discounted_heat_produced: float
		utc: float


		@dataclass
		class EconomicsResults:
		class EconomicsResults(NamedTuple):
		capex: CapexCalculatorResults
		utc: UTCCalculatorResults

		@@ -108,22 +105,43 @@ def calculate_capex(
		total_opex_ts = np.zeros(n_years)
		total_capex_ts = np.zeros(n_years)

		shifted_heat_power = heat_power_produced.copy()
		shifted_heat_power = np.array(heat_power_produced, dtype=np.float64)
		shift_time_series(shifted_heat_power, props.drilling_time)

		capex_well = calculate_capex_for_wells(props, ah_length_array, stimulation_capex)
		allocate_capex_for_wells(props, total_capex_ts, capex_well)
		capex_well = calculate_capex_for_wells(
		inj_depth=ah_length_array[0],
		prod_depth=ah_length_array[0],
		stimulation_capex=stimulation_capex,
		well_cost_z2=props.well_cost_z2,
		well_cost_z=props.well_cost_z,
		well_cost_const=props.well_cost_const,
		well_cost_scaling=props.well_cost_scaling,
		)
		allocate_capex_for_wells(props.drilling_time, total_capex_ts, capex_well)

		hp_added_power = get_max_hp_added_power(
		initial_hp_added_power, hp_added_power_years
		hp_years_arr = (
		np.array(hp_added_power_years, dtype=np.float64)
		if hp_added_power_years is not None
		else np.zeros(0, dtype=np.float64)
		)

		hp_added_power = get_max_hp_added_power(initial_hp_added_power, hp_years_arr)
		hp_after_hp = calculate_hp_added_power_after_heat_pump(
		props, hp_added_power, hp_cop
		hp_added_power, hp_cop, props.hp_alternative_heating_price
		)
		installation_mw = calculate_installation_mw(
		float(shifted_heat_power[-1]), hp_added_power
		)
		installation_mw = calculate_installation_mw(shifted_heat_power, hp_added_power)

		total_capex_1year = calculate_total_capex_1year(
		props, total_capex_ts, hp_after_hp, installation_mw
		last_year=max(props.drilling_time - 1, 0),
		total_capex_ts=total_capex_ts,
		capex_const=props.capex_const,
		hp_capex=props.hp_capex,
		capex_variable=props.capex_variable,
		capex_contingency=props.capex_contingency,
		hp_after_hp=hp_after_hp,
		installation_mw=installation_mw,
		)

		sum_capex = process_annual_data(
		@@ -166,71 +184,84 @@ def shift_time_series(series: list[float], shift: int) -> None:
		series[i] = 0


		@njit
		def calculate_capex_for_wells(
		props: UTCConfiguration, ah_length_array: list[float], stimulation_capex: float
		inj_depth: float,
		prod_depth: float,
		stimulation_capex: float,
		well_cost_z2: float,
		well_cost_z: float,
		well_cost_const: float,
		well_cost_scaling: float,
		) -> float:
		inj_depth = ah_length_array[0]
		prod_depth = ah_length_array[0]

		capex = (
		calculate_well_cost(props, inj_depth)
		+ calculate_well_cost(props, prod_depth)
		calculate_well_cost(well_cost_z2, well_cost_z, well_cost_const, inj_depth)
		+ calculate_well_cost(well_cost_z2, well_cost_z, well_cost_const, prod_depth)
		+ stimulation_capex
		)
		return capex * props.well_cost_scaling
		return capex * well_cost_scaling


		def calculate_well_cost(props: UTCConfiguration, depth: float) -> float:
		return (
		props.well_cost_z2 * depth*2 + props.well_cost_z depth
		) * 1e-6 + props.well_cost_const
		@njit
		def calculate_well_cost(
		well_cost_z2: float, well_cost_z: float, well_cost_const: float, depth: float
		) -> float:
		return (well_cost_z2 * depth*2 + well_cost_z depth) * 1e-6 + well_cost_const


		@njit
		def allocate_capex_for_wells(
		props: UTCConfiguration, total_capex_ts: NDArray[np.float64], capex_well: float
		):
		drilling_years = props.drilling_time
		drilling_years: int, total_capex_ts: np.ndarray, capex_well: float
		) -> None:
		yearly_cost = capex_well / max(drilling_years, 1)

		for y in range(drilling_years):
		total_capex_ts[y] += yearly_cost


		def get_max_hp_added_power(initial: float, hp_power_years: list[float] \| None) -> float:
		if hp_power_years is not None:
		return max(hp_power_years)
		@njit
		def get_max_hp_added_power(initial: float, hp_power_years: np.ndarray) -> float:
		if hp_power_years.size == 0:
		return initial
		max_val = hp_power_years[0]
		for i in range(hp_power_years.size):
		if hp_power_years[i] > max_val:
		max_val = hp_power_years[i]
		return max_val


		@njit
		def calculate_hp_added_power_after_heat_pump(
		props: UTCConfiguration, hp_added_power: float, hp_cop: float
		hp_added_power: float, hp_cop: float, hp_alternative_heating_price: float
		) -> float:
		alt_price = props.hp_alternative_heating_price
		return hp_added_power if alt_price < 0 else hp_added_power * (1 + 1 / (hp_cop - 1))
		if hp_alternative_heating_price < 0:
		return hp_added_power
		return hp_added_power * (1 + 1.0 / (hp_cop - 1.0))


		@njit
		def calculate_installation_mw(
		heat_power_produced: list[float], hp_added_power: float
		last_heat_power_value: float, hp_added_power: float
		) -> float:
		last_val = heat_power_produced[-1]
		return max(last_val - hp_added_power, 0)
		diff = last_heat_power_value - hp_added_power
		return diff if diff > 0 else 0.0


		@njit
		def calculate_total_capex_1year(
		props: UTCConfiguration,
		total_capex_ts: NDArray[np.float64],
		last_year: int,
		total_capex_ts: np.ndarray,
		capex_const: float,
		hp_capex: float,
		capex_variable: float,
		capex_contingency: float,
		hp_after_hp: float,
		installation_mw: float,
		) -> float:
		last_year = max(props.drilling_time - 1, 0)

		total_capex_ts[last_year] += (
		props.capex_const
		+ (props.hp_capex * hp_after_hp + props.capex_variable * installation_mw) * 1e-3
		capex_const + (hp_capex * hp_after_hp + capex_variable * installation_mw) * 1e-3
		)

		total_capex_ts[last_year] *= 1 + props.capex_contingency / 100

		total_capex_ts[last_year] *= 1.0 + capex_contingency / 100.0
		return total_capex_ts[last_year]


		@@ -277,7 +308,6 @@ def process_annual_data(

		# Variable OPEX
		variable_opex[year] = calculate_variable_opex(
		props,
		season_factor,
		props.elec_purchase_price,
		pump_power_required,
		@@ -286,14 +316,17 @@ def process_annual_data(
		heat_power_per_year[year],
		year_hp_cop,
		hp_after_hp,
		props.opex_per_energy,
		)

		# Fixed OPEX
		fixed_opex[year] = calculate_fixed_opex(
		props,
		props.opex_base,
		props.hp_opex,
		props.opex_per_power,
		hp_after_hp,
		installation_mw,
		total_capex_1year,
		props.opex_per_capex,
		)

		total_opex_ts[year] = variable_opex[year] + fixed_opex[year]
		@@ -301,8 +334,8 @@ def process_annual_data(
		return sum_capex


		@njit
		def calculate_variable_opex(
		props: UTCConfiguration,
		season_factor: float,
		elec_price: float,
		pump_power: float,
		@@ -311,6 +344,7 @@ def calculate_variable_opex(
		heat_power_gj_per_year: float,
		hp_cop: float,
		hp_after_hp: float,
		opex_per_energy: float,
		) -> float:
		pump_cost = (
		-(pump_power * 1e3)
		@@ -330,14 +364,8 @@ def calculate_variable_opex(
		* 1e-6
		)

		opex_per_energy = (
		-(1 / 0.36)
		* props.opex_per_energy
		* 1e-2
		* heat_power_gj_per_year
		* 1e4
		/ 36
		* 1e-6
		opex_energy_cost = (
		-(1 / 0.36) * opex_per_energy * 1e-2 * heat_power_gj_per_year * 1e4 / 36 * 1e-6
		)

		heat_pump_cost = 0.0
		@@ -351,19 +379,23 @@ def calculate_variable_opex(
		* 1e-6
		)

		return pump_cost + parasitic_cost + heat_pump_cost + opex_per_energy
		return pump_cost + parasitic_cost + heat_pump_cost + opex_energy_cost


		@njit
		def calculate_fixed_opex(
		props: UTCConfiguration,
		opex_base: float,
		hp_opex: float,
		opex_per_power: float,
		hp_after_hp: float,
		installation_mw: float,
		total_capex: float,
		opex_per_capex: float,
		) -> float:
		return (
		-props.opex_base / 1e6
		- (props.hp_opex * hp_after_hp + props.opex_per_power * installation_mw) * 1e-3
		- total_capex * props.opex_per_capex / 100
		-opex_base / 1e6
		- (hp_opex * hp_after_hp + opex_per_power * installation_mw) * 1e-3
		- total_capex * opex_per_capex / 100.0
		)

src/pythermogis/workflow/utc/flow.py

+2 −4

Original line number	Diff line number	Diff line
		from __future__ import annotations

		from dataclasses import dataclass
		from typing import TYPE_CHECKING
		from typing import TYPE_CHECKING, NamedTuple

		import numpy as np

		@@ -14,8 +13,7 @@ if TYPE_CHECKING:
		from pythermogis.workflow.utc.utc_properties import UTCConfiguration


		@dataclass
		class VolumetricFlowResults:
		class VolumetricFlowResults(NamedTuple):
		hp_cop: float
		hp_added_power: float
		hp_elec_consumption: float