From 7cacc5fcbdabec5943e31b1ceb2711b90519eae2 Mon Sep 17 00:00:00 2001 From: knappersfy Date: Tue, 9 Dec 2025 14:10:44 +0100 Subject: [PATCH 1/3] more efficiency: use named tuple instead of Dataclass --- src/pythermogis/workflow/utc/doublet.py | 6 +++--- src/pythermogis/workflow/utc/doubletcalc.py | 6 ++---- src/pythermogis/workflow/utc/economics.py | 12 ++++-------- src/pythermogis/workflow/utc/flow.py | 6 ++---- src/pythermogis/workflow/utc/heatpump.py | 6 ++---- src/pythermogis/workflow/utc/pressure.py | 6 ++---- src/pythermogis/workflow/utc/utc_properties.py | 2 +- 7 files changed, 16 insertions(+), 28 deletions(-) diff --git a/src/pythermogis/workflow/utc/doublet.py b/src/pythermogis/workflow/utc/doublet.py index 7dcc967..43cb4ff 100644 --- a/src/pythermogis/workflow/utc/doublet.py +++ b/src/pythermogis/workflow/utc/doublet.py @@ -1,5 +1,6 @@ import timeit from dataclasses import dataclass +from typing import NamedTuple import numpy as np @@ -17,7 +18,7 @@ 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 @@ -34,8 +35,7 @@ class DoubletInput: return self.transmissivity / self.thickness * 1e-3 * DARCY_SI -@dataclass -class DoubletOutput: +class DoubletOutput(NamedTuple): power: float hppower: float capex: float diff --git a/src/pythermogis/workflow/utc/doubletcalc.py b/src/pythermogis/workflow/utc/doubletcalc.py index cf9f7c2..78fc61f 100644 --- a/src/pythermogis/workflow/utc/doubletcalc.py +++ b/src/pythermogis/workflow/utc/doubletcalc.py @@ -1,8 +1,7 @@ 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 diff --git a/src/pythermogis/workflow/utc/economics.py b/src/pythermogis/workflow/utc/economics.py index 17d82c9..c959d11 100644 --- a/src/pythermogis/workflow/utc/economics.py +++ b/src/pythermogis/workflow/utc/economics.py @@ -1,7 +1,6 @@ from __future__ import annotations -from dataclasses import dataclass -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, NamedTuple import numpy as np from numpy.typing import NDArray @@ -18,8 +17,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 +26,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 diff --git a/src/pythermogis/workflow/utc/flow.py b/src/pythermogis/workflow/utc/flow.py index 6af7099..7fc51d2 100644 --- a/src/pythermogis/workflow/utc/flow.py +++ b/src/pythermogis/workflow/utc/flow.py @@ -1,7 +1,6 @@ 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 diff --git a/src/pythermogis/workflow/utc/heatpump.py b/src/pythermogis/workflow/utc/heatpump.py index 37ce497..b2d7a77 100644 --- a/src/pythermogis/workflow/utc/heatpump.py +++ b/src/pythermogis/workflow/utc/heatpump.py @@ -1,7 +1,6 @@ from __future__ import annotations -from dataclasses import dataclass -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, NamedTuple from pythermogis.workflow.utc.doublet_utils import get_cop_carnot from pythermogis.workflow.utc.water import ( @@ -16,8 +15,7 @@ if TYPE_CHECKING: from pythermogis.workflow.utc.utc_properties import UTCConfiguration -@dataclass -class HeatPumpPerformanceResults: +class HeatPumpPerformanceResults(NamedTuple): hp_cop: float hp_added_power: float hp_elec_consumption: float diff --git a/src/pythermogis/workflow/utc/pressure.py b/src/pythermogis/workflow/utc/pressure.py index 54a495d..244bab4 100644 --- a/src/pythermogis/workflow/utc/pressure.py +++ b/src/pythermogis/workflow/utc/pressure.py @@ -1,7 +1,6 @@ from __future__ import annotations -from dataclasses import dataclass -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, NamedTuple from pythermogis.workflow.utc.economics import calculate_economics from pythermogis.workflow.utc.flow import calculate_volumetric_flow @@ -75,8 +74,7 @@ def calculate_max_pressure( return pressure -@dataclass -class PressureOptimizerResults: +class PressureOptimizerResults(NamedTuple): drawdown_pressure: float flowrate: float heat_power_per_doublet: float diff --git a/src/pythermogis/workflow/utc/utc_properties.py b/src/pythermogis/workflow/utc/utc_properties.py index 1a75ce9..b9f9399 100644 --- a/src/pythermogis/workflow/utc/utc_properties.py +++ b/src/pythermogis/workflow/utc/utc_properties.py @@ -17,7 +17,7 @@ class AquiferFile(NamedTuple): ViscosityMode = Literal["kestin", "batzlewang"] -@dataclass(frozen=True) +@dataclass(slots=True, frozen=True) class UTCConfiguration: input_data_dir: str = "" results_dir: str = "" -- GitLab From c84a522c2c6058c30f7570fc0b3dd3cfb3bfe2c0 Mon Sep 17 00:00:00 2001 From: knappersfy Date: Tue, 9 Dec 2025 15:01:43 +0100 Subject: [PATCH 2/3] numba in economics --- src/pythermogis/workflow/utc/economics.py | 152 +++++++++++++--------- 1 file changed, 94 insertions(+), 58 deletions(-) diff --git a/src/pythermogis/workflow/utc/economics.py b/src/pythermogis/workflow/utc/economics.py index c959d11..de43298 100644 --- a/src/pythermogis/workflow/utc/economics.py +++ b/src/pythermogis/workflow/utc/economics.py @@ -3,6 +3,7 @@ from __future__ import annotations 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 @@ -104,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( @@ -162,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) - return initial +@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] @@ -273,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, @@ -282,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] @@ -297,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, @@ -307,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) @@ -326,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 @@ -347,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 ) -- GitLab From 66a171363b148616eed3e649ab85dadbd626f5cf Mon Sep 17 00:00:00 2001 From: knappersfy Date: Tue, 9 Dec 2025 15:19:31 +0100 Subject: [PATCH 3/3] start constants, and no redifining results --- src/pythermogis/workflow/utc/constants.py | 1 + src/pythermogis/workflow/utc/doublet.py | 48 ++++++++--------------- 2 files changed, 17 insertions(+), 32 deletions(-) create mode 100644 src/pythermogis/workflow/utc/constants.py diff --git a/src/pythermogis/workflow/utc/constants.py b/src/pythermogis/workflow/utc/constants.py new file mode 100644 index 0000000..5e05f1f --- /dev/null +++ b/src/pythermogis/workflow/utc/constants.py @@ -0,0 +1 @@ +DARCY_SI = 1.0e-12 / 1.01325 \ No newline at end of file diff --git a/src/pythermogis/workflow/utc/doublet.py b/src/pythermogis/workflow/utc/doublet.py index 43cb4ff..33d6441 100644 --- a/src/pythermogis/workflow/utc/doublet.py +++ b/src/pythermogis/workflow/utc/doublet.py @@ -2,8 +2,6 @@ import timeit from dataclasses import dataclass from typing import NamedTuple -import numpy as np - from pythermogis.utils.timer import print_time from pythermogis.workflow.utc.cooling_temp import calculate_cooling_temperature from pythermogis.workflow.utc.doublet_utils import ( @@ -13,6 +11,7 @@ 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 @@ -31,7 +30,6 @@ class DoubletInput: @property def permeability(self) -> float: - DARCY_SI = 1.0e-12 / 1.01325 return self.transmissivity / self.thickness * 1e-3 * DARCY_SI @@ -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, ) -- GitLab