diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index fa541b25..828d421e 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -5,10 +5,12 @@ on: branches: - master - maintenance/2.2.x + - maintenance/2.3.x pull_request: branches: - master - maintenance/2.2.x + - maintenance/2.3.x types: [opened, synchronize, reopened] jobs: diff --git a/CHANGELOG.md b/CHANGELOG.md index d7d5b95d..ecf629be 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,16 @@ # History of changes +## Version 2.3.1 (2025-08-04) + +### New features + +* [Pull Request 325](https://github.com/MassimoCimmino/pygfunction/pull/325) - Borefields and boreholes can now be concatenated using the `+` operator, e.g. using `new_field = field_1 + field_2`. +* [Pull Request 326](https://github.com/MassimoCimmino/pygfunction/pull/326) - Introduced `gFunction.from_static_params` and `Network.from_static_params` methods. These methods facilitate the creation of `Network` objects and the evaluation of g-functions by automatically evaluating the required thermal resistances for the creation of `Pipe` objects. + +### Other changes + +* [Issue 319](https://github.com/MassimoCimmino/pygfunction/issues/319) - Created `solvers` module. `Solver` classes are moved out of the `gfunction` module and into the new module. + ## Version 2.3 (2025-04-29) ### New features diff --git a/doc/source/conf.py b/doc/source/conf.py index 95ac20d3..05361428 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -76,7 +76,7 @@ # The short X.Y version. version = u'2.3' # The full version, including alpha/beta/rc tags. -release = u'2.3.1.dev0' +release = u'2.3.1' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/doc/source/modules.rst b/doc/source/modules.rst index ad2f6a80..bfce80b6 100644 --- a/doc/source/modules.rst +++ b/doc/source/modules.rst @@ -14,4 +14,5 @@ Modules modules/media modules/networks modules/pipes + modules/solvers modules/utilities diff --git a/doc/source/modules/solvers.rst b/doc/source/modules/solvers.rst new file mode 100644 index 00000000..7603914a --- /dev/null +++ b/doc/source/modules/solvers.rst @@ -0,0 +1,9 @@ +.. solvers: + +************** +Solvers Module +************** + +.. automodule:: pygfunction.solvers + :members: + :show-inheritance: diff --git a/pygfunction/__init__.py b/pygfunction/__init__.py index 9b25af59..1c44ec66 100644 --- a/pygfunction/__init__.py +++ b/pygfunction/__init__.py @@ -7,5 +7,6 @@ media, networks, pipes, + solvers, utilities, ) diff --git a/pygfunction/borefield.py b/pygfunction/borefield.py index 0a6b7dbc..3a92b888 100644 --- a/pygfunction/borefield.py +++ b/pygfunction/borefield.py @@ -107,6 +107,42 @@ def __ne__( check = not self == other_field return check + def __add__(self, + other_field: Union[Borehole, List[Borehole], Self]) -> Self: + """Add two borefields together""" + if not isinstance(other_field, (Borehole, list, self.__class__)): + raise TypeError( + f'Expected Borefield, list or Borehole input;' + f' got {other_field}' + ) + # List of boreholes + field = self.to_boreholes() + # Convert other_field to a list if it is a Borehole + if isinstance(other_field, Borehole): + other_field = [other_field] + # Convert borefield to a list if it is a Borefield + if isinstance(other_field, self.__class__): + other_field = other_field.to_boreholes() + return Borefield.from_boreholes(field + other_field) + + def __radd__(self, + other_field: Union[Borehole, List[Borehole], Self]) -> Self: + """Add two borefields together""" + if not isinstance(other_field, (Borehole, list, self.__class__)): + raise TypeError( + f'Expected Borefield, list or Borehole input;' + f' got {other_field}' + ) + # List of boreholes + field = self.to_boreholes() + # Convert other_field to a list if it is a Borehole + if isinstance(other_field, Borehole): + other_field = [other_field] + # Convert borefield to a list if it is a Borefield + if isinstance(other_field, self.__class__): + other_field = other_field.to_boreholes() + return Borefield.from_boreholes(other_field + field) + def evaluate_g_function( self, alpha: float, diff --git a/pygfunction/boreholes.py b/pygfunction/boreholes.py index 58c10013..330ed4b7 100644 --- a/pygfunction/boreholes.py +++ b/pygfunction/boreholes.py @@ -1,8 +1,10 @@ # -*- coding: utf-8 -*- +from typing import Union import warnings import numpy as np from scipy.spatial.distance import pdist +from typing_extensions import Self # for compatibility with Python <= 3.10 from .utilities import _initialize_figure, _format_axes, _format_axes_3d @@ -54,6 +56,56 @@ def __repr__(self): f' orientation={self.orientation})') return s + def __add__(self, other: Union[Self, list]): + """ + Adds two boreholes together to form a borefield + """ + if not isinstance(other, (self.__class__, list)): + # Check if other is a borefield and try the operation using + # other.__radd__ + try: + field = other.__radd__(self) + except: + # Invalid input + raise TypeError( + f'Expected Borefield, list or Borehole input;' + f' got {other}' + ) + elif isinstance(other, list): + # Create a borefield from the borehole and a list + from .borefield import Borefield + field = Borefield.from_boreholes([self] + other) + else: + # Create a borefield from the two boreholes + from .borefield import Borefield + field = Borefield.from_boreholes([self, other]) + return field + + def __radd__(self, other: Union[Self, list]): + """ + Adds two boreholes together to form a borefield + """ + if not isinstance(other, (self.__class__, list)): + # Check if other is a borefield and try the operation using + # other.__radd__ + try: + field = other.__add__(self) + except: + # Invalid input + raise TypeError( + f'Expected Borefield, list or Borehole input;' + f' got {other}' + ) + elif isinstance(other, list): + # Create a borefield from the borehole and a list + from .borefield import Borefield + field = Borefield.from_boreholes(other + [self]) + else: + # Create a borefield from the two boreholes + from .borefield import Borefield + field = Borefield.from_boreholes([other, self]) + return field + def distance(self, target): """ Evaluate the distance between the current borehole and a target diff --git a/pygfunction/enums.py b/pygfunction/enums.py new file mode 100644 index 00000000..cf19f2fe --- /dev/null +++ b/pygfunction/enums.py @@ -0,0 +1,21 @@ +from enum import Enum, auto +from typing import Annotated + + +class PipeType(Enum): + """Enumerator for pipe configuration type.""" + COAXIAL_ANNULAR_IN: Annotated[ + int, "Coaxial pipe (annular inlet)" + ] = auto() + COAXIAL_ANNULAR_OUT: Annotated[ + int, "Coaxial pipe (annular outlet)" + ] = auto() + DOUBLE_UTUBE_PARALLEL: Annotated[ + int, "Double U-tube (parallel)" + ] = auto() + DOUBLE_UTUBE_SERIES: Annotated[ + int, "Double U-tube (series)" + ] = auto() + SINGLE_UTUBE: Annotated[ + int, "Single U-tube" + ] = auto() diff --git a/pygfunction/gfunction.py b/pygfunction/gfunction.py index 5cae8bd3..0edad340 100644 --- a/pygfunction/gfunction.py +++ b/pygfunction/gfunction.py @@ -1,19 +1,26 @@ # -*- coding: utf-8 -*- import warnings from time import perf_counter +from typing import Union, List import numpy as np -from scipy.cluster.hierarchy import cut_tree, dendrogram, linkage -from scipy.constants import pi +import numpy.typing as npt from scipy.interpolate import interp1d as interp1d from .borefield import Borefield -from .boreholes import Borehole, _EquivalentBorehole, find_duplicates -from .heat_transfer import finite_line_source, finite_line_source_vectorized, \ - finite_line_source_equivalent_boreholes_vectorized, \ - finite_line_source_inclined_vectorized -from .networks import Network, _EquivalentNetwork, network_thermal_resistance -from .utilities import _initialize_figure, _format_axes, segment_ratios +from .boreholes import Borehole, find_duplicates +from .media import Fluid +from .networks import Network +from .solvers import ( + Detailed, + Equivalent, + Similarities +) +from .utilities import ( + segment_ratios, + _initialize_figure, + _format_axes +) class gFunction(object): @@ -257,17 +264,17 @@ def __init__(self, boreholes_or_network, alpha, time=None, # Load the chosen solver if self.method.lower()=='similarities': - self.solver = _Similarities( + self.solver = Similarities( self.boreholes, self.network, self.time, self.boundary_condition, self.m_flow_borehole, self.m_flow_network, self.cp_f, **self.options) elif self.method.lower()=='detailed': - self.solver = _Detailed( + self.solver = Detailed( self.boreholes, self.network, self.time, self.boundary_condition, self.m_flow_borehole, self.m_flow_network, self.cp_f, **self.options) elif self.method.lower()=='equivalent': - self.solver = _Equivalent( + self.solver = Equivalent( self.boreholes, self.network, self.time, self.boundary_condition, self.m_flow_borehole, self.m_flow_network, self.cp_f, **self.options) @@ -278,6 +285,172 @@ def __init__(self, boreholes_or_network, alpha, time=None, if self.time is not None: self.gFunc = self.evaluate_g_function(self.time) + @classmethod + def from_static_params(cls, + H: npt.ArrayLike, + D: npt.ArrayLike, + r_b: npt.ArrayLike, + x: npt.ArrayLike, + y: npt.ArrayLike, + alpha: float, + time: npt.ArrayLike = None, + method: str = 'equivalent', + m_flow_network: float = None, + options={}, + tilt: npt.ArrayLike = 0., + orientation: npt.ArrayLike = 0., + boundary_condition: str = 'MIFT', + pipe_type_str: str = None, + pos: List[tuple] = None, + r_in: Union[float, tuple, npt.ArrayLike] = None, + r_out: Union[float, tuple, npt.ArrayLike] = None, + k_s: float = None, + k_g: float = None, + k_p: Union[float, tuple, npt.ArrayLike] = None, + fluid_str: str = None, + fluid_concentration_pct: float = None, + fluid_temperature: float = 20, + epsilon: float = None, + reversible_flow: bool = True, + bore_connectivity: list = None, + J: int = 2, + ): + + """ + Constructs the 'gFunction' class from static parameters. + + Parameters + ---------- + H : float or (nBoreholes,) array + Borehole lengths (in meters). + D : float or (nBoreholes,) array + Borehole buried depths (in meters). + r_b : float or (nBoreholes,) array + Borehole radii (in meters). + x : float or (nBoreholes,) array + Position (in meters) of the head of the boreholes along the x-axis. + y : float or (nBoreholes,) array + Position (in meters) of the head of the boreholes along the y-axis. + alpha : float + Soil thermal diffusivity (in m2/s). + time : float or array, optional + Values of time (in seconds) for which the g-function is evaluated. The + g-function is only evaluated at initialization if a value is provided. + Default is None. + method : str, optional + Method for the evaluation of the g-function. Should be one of 'similarities', 'detailed', or 'equivalent'. + Default is 'equivalent'. See 'gFunction' __init__ for more details. + m_flow_network : float, optional + Fluid mass flow rate into the network of boreholes (in kg/s). + Default is None. + options : dict, optional + A dictionary of solver options. See 'gFunction' __init__ for more details. + tilt : float or (nBoreholes,) array, optional + Angle (in radians) from vertical of the axis of the boreholes. + Default is 0. + orientation : float or (nBoreholes,) array, optional + Direction (in radians) of the tilt of the boreholes. Defaults to zero + if the borehole is vertical. + Default is 0. + boundary_condition : str, optional + Boundary condition for the evaluation of the g-function. Should be one of 'UHTR', 'UBWT', or 'MIFT'. + Default is 'MIFT'. + pipe_type_str : str, optional + Pipe type used for 'MIFT' boundary condition. Should be one of 'COAXIAL_ANNULAR_IN', 'COAXIAL_ANNULAR_OUT', + 'DOUBLE_UTUBE_PARALLEL', 'DOUBLE_UTUBE_SERIES', or 'SINGLE_UTUBE'. + pos : list of tuples, optional + Position (x, y) (in meters) of the pipes inside the borehole. + r_in : float, optional + Inner radius (in meters) of the U-Tube pipes. + r_out : float, optional + Outer radius (in meters) of the U-Tube pipes. + k_s : float, optional + Soil thermal conductivity (in W/m-K). + k_g : float, optional + Grout thermal conductivity (in W/m-K). + k_p : float, optional + Pipe thermal conductivity (in W/m-K). + fluid_str: str, optional + The mixer for this application should be one of: + + - 'Water' - Complete water solution + - 'MEG' - Ethylene glycol mixed with water + - 'MPG' - Propylene glycol mixed with water + - 'MEA' - Ethanol mixed with water + - 'MMA' - Methanol mixed with water + + fluid_concentration_pct: float, optional + Mass fraction of the mixing fluid added to water (in %). + Lower bound = 0. Upper bound is dependent on the mixture. + fluid_temperature: float, optional + Temperature used for evaluating fluid properties (in degC). + Default is 20. + epsilon : float, optional + Pipe roughness (in meters). + reversible_flow : bool, optional + True to treat a negative mass flow rate as the reversal of flow + direction within the borehole. If False, the direction of flow is not + reversed when the mass flow rate is negative, and the absolute value is + used for calculations. + Default True. + bore_connectivity : list, optional + Index of fluid inlet into each borehole. -1 corresponds to a borehole + connected to the bore field inlet. If this parameter is not provided, + parallel connections between boreholes is used. + Default is None. + J : int, optional + Number of multipoles per pipe to evaluate the thermal resistances. + J=1 or J=2 usually gives sufficient accuracy. J=0 corresponds to the + line source approximation. + Default is 2. + + Returns + ------- + gFunction : 'gFunction' object + The g-function. + + Notes + ----- + - When using the 'MIFT' boundary condition, the parameters `pipe_type_str`, + `fluid_str`, `fluid_concentration_pct`, `fluid_temperature`, and + `epsilon` are required. + + """ + + if boundary_condition.upper() not in ['UBWT', 'UHTR', 'MIFT']: + raise ValueError(f"'{boundary_condition}' is not a valid boundary condition.") + + # construct all required pieces + borefield = Borefield(H, D, r_b, x, y, tilt, orientation) + + if boundary_condition.upper() == 'MIFT': + boreholes = borefield.to_boreholes() + cp_f = Fluid(fluid_str, fluid_concentration_pct).cp + boreholes_or_network= Network.from_static_params( + boreholes, + pipe_type_str, + pos, + r_in, + r_out, + k_s, + k_g, + k_p, + m_flow_network, + epsilon, + fluid_str, + fluid_concentration_pct, + fluid_temperature, + reversible_flow, + bore_connectivity, + J, + ) + else: + boreholes_or_network = borefield + cp_f = None + + return cls(boreholes_or_network, alpha, time=time, method=method, boundary_condition=boundary_condition, + m_flow_network=m_flow_network, cp_f=cp_f, options=options) + def evaluate_g_function(self, time): """ Evaluate the g-function. @@ -516,7 +689,7 @@ def visualize_heat_extraction_rates( # Adjust figure to window plt.tight_layout() - + return fig def visualize_heat_extraction_rate_profiles( @@ -1762,3061 +1935,3 @@ def mixed_inlet_temperature( boundary_condition=boundary_condition, options=options) return gFunc.gFunc - - -class _BaseSolver(object): - """ - Template for solver classes. - - Solver classes inherit from this class. - - Attributes - ---------- - boreholes : list of Borehole objects - List of boreholes included in the bore field. - network : network object - Model of the network. - time : float or array - Values of time (in seconds) for which the g-function is evaluated. - boundary_condition : str - Boundary condition for the evaluation of the g-function. Should be one - of - - - 'UHTR' : - Uniform heat transfer rate. - - 'UBWT' : - Uniform borehole wall temperature. - - 'MIFT' : - Mixed inlet fluid temperatures. - - nSegments : int or list, optional - Number of line segments used per borehole, or list of number of - line segments used for each borehole. - Default is 8. - segment_ratios : array, list of arrays, or callable, optional - Ratio of the borehole length represented by each segment. The - sum of ratios must be equal to 1. The shape of the array is of - (nSegments,) or list of (nSegments[i],). If segment_ratios==None, - segments of equal lengths are considered. If a callable is provided, it - must return an array of size (nSegments,) when provided with nSegments - (of type int) as an argument, or an array of size (nSegments[i],) when - provided with an element of nSegments (of type list). - Default is :func:`utilities.segment_ratios`. - m_flow_borehole : (nInlets,) array or (nMassFlow, nInlets,) array, optional - Fluid mass flow rate into each circuit of the network. If a - (nMassFlow, nInlets,) array is supplied, the - (nMassFlow, nMassFlow,) variable mass flow rate g-functions - will be evaluated using the method of Cimmino (2024) - [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary - condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be - provided. - Default is None. - m_flow_network : float or (nMassFlow,) array, optional - Fluid mass flow rate into the network of boreholes. If an array - is supplied, the (nMassFlow, nMassFlow,) variable mass flow - rate g-functions will be evaluated using the method of Cimmino - (2024) [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary - condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be - provided. - Default is None. - cp_f : float, optional - Fluid specific isobaric heat capacity (in J/kg.degC). Only required - for the 'MIFT' boundary condition. - Default is None. - approximate_FLS : bool, optional - Set to true to use the approximation of the FLS solution of Cimmino - (2021). This approximation does not require the numerical evaluation of - any integral. When using the 'equivalent' solver, the approximation is - only applied to the thermal response at the borehole radius. Thermal - interaction between boreholes is evaluated using the FLS solution. - Default is False. - nFLS : int, optional - Number of terms in the approximation of the FLS solution. This - parameter is unused if `approximate_FLS` is set to False. - Default is 10. Maximum is 25. - mQuad : int, optional - Number of Gauss-Legendre sample points for the integral over :math:`u` - in the inclined FLS solution. - Default is 11. - linear_threshold : float, optional - Threshold time (in seconds) under which the g-function is - linearized. The g-function value is then interpolated between 0 - and its value at the threshold. If linear_threshold==None, the - g-function is linearized for times - `t < r_b**2 / (25 * self.alpha)`. - Default is None. - disp : bool, optional - Set to true to print progression messages. - Default is False. - profiles : bool, optional - Set to true to keep in memory the temperatures and heat extraction - rates. - Default is False. - kind : string, optional - Interpolation method used for segment-to-segment thermal response - factors. See documentation for scipy.interpolate.interp1d. - Default is 'linear'. - dtype : numpy dtype, optional - numpy data type used for matrices and vectors. Should be one of - numpy.single or numpy.double. - Default is numpy.double. - - """ - def __init__(self, boreholes, network, time, boundary_condition, - m_flow_borehole=None, m_flow_network=None, cp_f=None, - nSegments=8, segment_ratios=segment_ratios, - approximate_FLS=False, mQuad=11, nFLS=10, - linear_threshold=None, disp=False, profiles=False, - kind='linear', dtype=np.double, **other_options): - self.boreholes = boreholes - self.network = network - # Convert time to a 1d array - self.time = np.atleast_1d(time).flatten() - self.linear_threshold = linear_threshold - self.r_b_max = np.max([b.r_b for b in self.boreholes]) - self.boundary_condition = boundary_condition - nBoreholes = len(self.boreholes) - # Format number of segments and segment ratios - if type(nSegments) is int: - self.nBoreSegments = [nSegments] * nBoreholes - else: - self.nBoreSegments = nSegments - if isinstance(segment_ratios, np.ndarray): - segment_ratios = [segment_ratios] * nBoreholes - elif segment_ratios is None: - segment_ratios = [np.full(n, 1./n) for n in self.nBoreSegments] - elif callable(segment_ratios): - segment_ratios = [segment_ratios(n) for n in self.nBoreSegments] - self.segment_ratios = segment_ratios - # Shortcut for segment_ratios comparisons - self._equal_segment_ratios = \ - (np.all(np.array(self.nBoreSegments, dtype=np.uint) == self.nBoreSegments[0]) - and np.all([np.allclose(segment_ratios, self.segment_ratios[0]) for segment_ratios in self.segment_ratios])) - # Boreholes with a uniform discretization - self._uniform_segment_ratios = [ - np.allclose(segment_ratios, - segment_ratios[0:1], - rtol=1e-6) - for segment_ratios in self.segment_ratios] - # Find indices of first and last segments along boreholes - self._i0Segments = [sum(self.nBoreSegments[0:i]) - for i in range(nBoreholes)] - self._i1Segments = [sum(self.nBoreSegments[0:(i + 1)]) - for i in range(nBoreholes)] - self.nMassFlow = 0 - self.m_flow_borehole = m_flow_borehole - if self.m_flow_borehole is not None: - if not self.m_flow_borehole.ndim == 1: - self.nMassFlow = np.size(self.m_flow_borehole, axis=0) - self.m_flow_borehole = np.atleast_2d(self.m_flow_borehole) - self.m_flow = self.m_flow_borehole - self.m_flow_network = m_flow_network - if self.m_flow_network is not None: - if not isinstance(self.m_flow_network, (np.floating, float)): - self.nMassFlow = len(self.m_flow_network) - self.m_flow_network = np.atleast_1d(self.m_flow_network) - self.m_flow = self.m_flow_network - self.cp_f = cp_f - self.approximate_FLS = approximate_FLS - self.mQuad = mQuad - self.nFLS = nFLS - self.disp = disp - self.profiles = profiles - self.kind = kind - self.dtype = dtype - # Check the validity of inputs - self._check_inputs() - # Initialize the solver with solver-specific options - self.nSources = self.initialize(**other_options) - - return - - def initialize(self, *kwargs): - """ - Perform any calculation required at the initialization of the solver - and returns the number of finite line heat sources in the borefield. - - Raises - ------ - NotImplementedError - - Returns - ------- - nSources : int - Number of finite line heat sources in the borefield used to - initialize the matrix of segment-to-segment thermal response - factors (of size: nSources x nSources). - - """ - raise NotImplementedError( - 'initialize class method not implemented, this method should ' - 'return the number of finite line heat sources in the borefield ' - 'used to initialize the matrix of segment-to-segment thermal ' - 'response factors (of size: nSources x nSources)') - return None - - def solve(self, time, alpha): - """ - Build and solve the system of equations. - - Parameters - ---------- - time : float or array - Values of time (in seconds) for which the g-function is evaluated. - alpha : float - Soil thermal diffusivity (in m2/s). - - Returns - ------- - gFunc : float or array - Values of the g-function - - """ - # Number of time values - self.time = time - nt = len(self.time) - # Evaluate threshold time for g-function linearization - if self.linear_threshold is None: - time_threshold = self.r_b_max**2 / (25 * alpha) - else: - time_threshold = self.linear_threshold - # Find the number of g-function values to be linearized - p_long = np.searchsorted(self.time, time_threshold, side='right') - if p_long > 0: - time_long = np.concatenate([[time_threshold], self.time[p_long:]]) - else: - time_long = self.time - nt_long = len(time_long) - # Calculate segment to segment thermal response factors - h_ij = self.thermal_response_factors(time_long, alpha, kind=self.kind) - # Segment lengths - H_b = self.segment_lengths() - if self.boundary_condition == 'MIFT': - Hb_individual = np.array([b.H for b in self.boreSegments], dtype=self.dtype) - H_tot = np.sum(H_b) - if self.disp: print('Building and solving the system of equations ...', - end='') - # Initialize chrono - tic = perf_counter() - - if self.boundary_condition == 'UHTR': - # Initialize g-function - gFunc = np.zeros(nt) - # Initialize segment heat extraction rates - Q_b = 1 - # Initialize borehole wall temperatures - T_b = np.zeros((self.nSources, nt), dtype=self.dtype) - - # Build and solve the system of equations at all times - p0 = max(0, p_long-1) - for p in range(nt_long): - # Evaluate the g-function with uniform heat extraction along - # boreholes - - # Thermal response factors evaluated at time t[p] - h_dt = h_ij.y[:,:,p+1] - # Borehole wall temperatures are calculated by the sum of - # contributions of all segments - T_b[:,p+p0] = np.sum(h_dt, axis=1) - # The g-function is the average of all borehole wall - # temperatures - gFunc[p+p0] = np.sum(T_b[:,p+p0]*H_b)/H_tot - - # Linearize g-function for times under threshold - if p_long > 0: - gFunc[:p_long] = gFunc[p_long-1] * self.time[:p_long] / time_threshold - T_b[:,:p_long] = T_b[:,p_long-1:p_long] * self.time[:p_long] / time_threshold - - elif self.boundary_condition == 'UBWT': - # Initialize g-function - gFunc = np.zeros(nt) - # Initialize segment heat extraction rates - Q_b = np.zeros((self.nSources, nt), dtype=self.dtype) - T_b = np.zeros(nt, dtype=self.dtype) - - # Build and solve the system of equations at all times - p0 = max(0, p_long-1) - for p in range(nt_long): - # Current thermal response factor matrix - if p > 0: - dt = time_long[p] - time_long[p-1] - else: - dt = time_long[p] - # Thermal response factors evaluated at t=dt - h_dt = h_ij(dt) - # Reconstructed load history - Q_reconstructed = self.load_history_reconstruction( - time_long[0:p+1], Q_b[:,p0:p+p0+1]) - # Borehole wall temperature for zero heat extraction at - # current step - T_b0 = self.temporal_superposition( - h_ij.y[:,:,1:], Q_reconstructed) - - # Evaluate the g-function with uniform borehole wall - # temperature - # --------------------------------------------------------- - # Build a system of equation [A]*[X] = [B] for the - # evaluation of the g-function. [A] is a coefficient - # matrix, [X] = [Q_b,T_b] is a state space vector of the - # borehole heat extraction rates and borehole wall - # temperature (equal for all segments), [B] is a - # coefficient vector. - # - # Spatial superposition: [T_b] = [T_b0] + [h_ij_dt]*[Q_b] - # Energy conservation: sum([Q_b*Hb]) = sum([Hb]) - # --------------------------------------------------------- - A = np.block([[h_dt, -np.ones((self.nSources, 1), - dtype=self.dtype)], - [H_b, 0.]]) - B = np.hstack((-T_b0, H_tot)) - # Solve the system of equations - X = np.linalg.solve(A, B) - # Store calculated heat extraction rates - Q_b[:,p+p0] = X[0:self.nSources] - # The borehole wall temperatures are equal for all segments - T_b[p+p0] = X[-1] - gFunc[p+p0] = T_b[p+p0] - - # Linearize g-function for times under threshold - if p_long > 0: - gFunc[:p_long] = gFunc[p_long-1] * self.time[:p_long] / time_threshold - Q_b[:,:p_long] = 1 + (Q_b[:,p_long-1:p_long] - 1) * self.time[:p_long] / time_threshold - T_b[:p_long] = T_b[p_long-1] * self.time[:p_long] / time_threshold - - elif self.boundary_condition == 'MIFT': - if self.nMassFlow == 0: - # Initialize g-function - gFunc = np.zeros((1, 1, nt)) - # Initialize segment heat extraction rates - Q_b = np.zeros((1, self.nSources, nt), dtype=self.dtype) - T_b = np.zeros((1, self.nSources, nt), dtype=self.dtype) - else: - # Initialize g-function - gFunc = np.zeros((self.nMassFlow, self.nMassFlow, nt)) - # Initialize segment heat extraction rates - Q_b = np.zeros( - (self.nMassFlow, self.nSources, nt), dtype=self.dtype) - T_b = np.zeros( - (self.nMassFlow, self.nSources, nt), dtype=self.dtype) - - for j in range(np.maximum(self.nMassFlow, 1)): - # Build and solve the system of equations at all times - p0 = max(0, p_long-1) - a_in_j, a_b_j = self.network.coefficients_borehole_heat_extraction_rate( - self.m_flow[j], - self.cp_f, - self.nBoreSegments, - segment_ratios=self.segment_ratios) - k_s = self.network.p[0].k_s - for p in range(nt_long): - # Current thermal response factor matrix - if p > 0: - dt = time_long[p] - time_long[p-1] - else: - dt = time_long[p] - # Thermal response factors evaluated at t=dt - h_dt = h_ij(dt) - # Reconstructed load history - Q_reconstructed = self.load_history_reconstruction( - time_long[0:p+1], Q_b[j,:,p0:p+p0+1]) - # Borehole wall temperature for zero heat extraction at - # current step - T_b0 = self.temporal_superposition( - h_ij.y[:,:,1:], Q_reconstructed) - - # Evaluate the g-function with mixed inlet fluid - # temperatures - # --------------------------------------------------------- - # Build a system of equation [A]*[X] = [B] for the - # evaluation of the g-function. [A] is a coefficient - # matrix, [X] = [Q_b,T_b,Tf_in] is a state space vector of - # the borehole heat extraction rates, borehole wall - # temperatures and inlet fluid temperature (into the bore - # field), [B] is a coefficient vector. - # - # Spatial superposition: [T_b] = [T_b0] + [h_ij_dt]*[Q_b] - # Heat transfer inside boreholes: - # [Q_{b,i}] = [a_in]*[T_{f,in}] + [a_{b,i}]*[T_{b,i}] - # Energy conservation: sum([Q_b*H_b]) = sum([H_b]) - # --------------------------------------------------------- - A = np.block( - [[h_dt, - -np.eye(self.nSources, dtype=self.dtype), - np.zeros((self.nSources, 1), dtype=self.dtype)], - [np.eye(self.nSources, dtype=self.dtype), - a_b_j/(2.0*pi*k_s*np.atleast_2d(Hb_individual).T), - a_in_j/(2.0*pi*k_s*np.atleast_2d(Hb_individual).T)], - [H_b, np.zeros(self.nSources + 1, dtype=self.dtype)]]) - B = np.hstack( - (-T_b0, - np.zeros(self.nSources, dtype=self.dtype), - H_tot)) - # Solve the system of equations - X = np.linalg.solve(A, B) - # Store calculated heat extraction rates - Q_b[j,:,p+p0] = X[0:self.nSources] - T_b[j,:,p+p0] = X[self.nSources:2*self.nSources] - # Inlet fluid temperature - T_f_in = X[-1] - # The gFunction is equal to the effective borehole wall - # temperature - # Outlet fluid temperature - T_f_out = T_f_in - 2 * pi * k_s * H_tot / ( - np.sum(np.abs(self.m_flow[j]) * self.cp_f)) - # Average fluid temperature - T_f = 0.5*(T_f_in + T_f_out) - # Borefield thermal resistance - R_field = network_thermal_resistance( - self.network, self.m_flow[j], self.cp_f) - # Effective borehole wall temperature - T_b_eff = T_f - 2 * pi * k_s * R_field - gFunc[j,j,p+p0] = T_b_eff - - for i in range(np.maximum(self.nMassFlow, 1)): - for j in range(np.maximum(self.nMassFlow, 1)): - if not i == j: - # Inlet fluid temperature - a_in, a_b = self.network.coefficients_network_heat_extraction_rate( - self.m_flow[i], - self.cp_f, - self.nBoreSegments, - segment_ratios=self.segment_ratios) - T_f_in = (-2 * pi * k_s * H_tot - a_b @ T_b[j,:,p0:]) / a_in - # The gFunction is equal to the effective borehole wall - # temperature - # Outlet fluid temperature - T_f_out = T_f_in - 2 * pi * k_s * H_tot / np.sum(np.abs(self.m_flow[i]) * self.cp_f) - # Borefield thermal resistance - R_field = network_thermal_resistance( - self.network, self.m_flow[i], self.cp_f) - # Effective borehole wall temperature - T_b_eff = 0.5 * (T_f_in + T_f_out) - 2 * pi * k_s * R_field - gFunc[i,j,p0:] = T_b_eff - - # Linearize g-function for times under threshold - if p_long > 0: - gFunc[:,:,:p_long] = gFunc[:,:,p_long-1] * self.time[:p_long] / time_threshold - Q_b[:,:,:p_long] = 1 + (Q_b[:,:,p_long-1:p_long] - 1) * self.time[:p_long] / time_threshold - T_b[:,:,:p_long] = T_b[:,:,p_long-1:p_long] * self.time[:p_long] / time_threshold - if self.nMassFlow == 0: - gFunc = gFunc[0,0,:] - Q_b = Q_b[0,:,:] - T_b = T_b[0,:,:] - - # Store temperature and heat extraction rate profiles - if self.profiles: - self.Q_b = Q_b - self.T_b = T_b - toc = perf_counter() - if self.disp: print(f' {toc - tic:.3f} sec') - return gFunc - - def segment_lengths(self): - """ - Return the length of all segments in the bore field. - - The segments lengths are used for the energy balance in the calculation - of the g-function. - - Returns - ------- - H : array - Array of segment lengths (in m). - - """ - # Borehole lengths - H_b = np.array([b.H for b in self.boreSegments], dtype=self.dtype) - return H_b - - def borehole_segments(self): - """ - Split boreholes into segments. - - This function goes through the list of boreholes and builds a new list, - with each borehole split into nSegments of equal lengths. - - Returns - ------- - boreSegments : list - List of borehole segments. - - """ - boreSegments = [] # list for storage of boreSegments - for b, nSegments, segment_ratios in zip(self.boreholes, self.nBoreSegments, self.segment_ratios): - segments = b.segments(nSegments, segment_ratios=segment_ratios) - boreSegments.extend(segments) - - return boreSegments - - def temporal_superposition(self, h_ij, Q_reconstructed): - """ - Temporal superposition for inequal time steps. - - Parameters - ---------- - h_ij : array - Values of the segment-to-segment thermal response factor increments - at the given time step. - Q_reconstructed : array - Reconstructed heat extraction rates of all segments at all times. - - Returns - ------- - T_b0 : array - Current values of borehole wall temperatures assuming no heat - extraction during current time step. - - """ - # Number of heat sources - nSources = Q_reconstructed.shape[0] - # Number of time steps - nt = Q_reconstructed.shape[1] - # Spatial and temporal superpositions - dQ = np.concatenate( - (Q_reconstructed[:,0:1], - Q_reconstructed[:,1:] - Q_reconstructed[:,0:-1]), axis=1)[:,::-1] - # Borehole wall temperature - T_b0 = np.einsum('ijk,jk', h_ij[:,:,:nt], dQ) - - return T_b0 - - def load_history_reconstruction(self, time, Q_b): - """ - Reconstructs the load history. - - This function calculates an equivalent load history for an inverted - order of time step sizes. - - Parameters - ---------- - time : array - Values of time (in seconds) in the load history. - Q_b : array - Heat extraction rates (in Watts) of all segments at all times. - - Returns - ------- - Q_reconstructed : array - Reconstructed load history. - - """ - # Number of heat sources - nSources = Q_b.shape[0] - # Time step sizes - dt = np.hstack((time[0], time[1:]-time[:-1])) - # Time vector - t = np.hstack((0., time, time[-1] + time[0])) - # Inverted time step sizes - dt_reconstructed = dt[::-1] - # Reconstructed time vector - t_reconstructed = np.hstack((0., np.cumsum(dt_reconstructed))) - # Accumulated heat extracted - f = np.hstack( - (np.zeros((nSources, 1), dtype=self.dtype), - np.cumsum(Q_b*dt, axis=1))) - f = np.hstack((f, f[:,-1:])) - # Create interpolation object for accumulated heat extracted - sf = interp1d(t, f, kind='linear', axis=1) - # Reconstructed load history - Q_reconstructed = (sf(t_reconstructed[1:]) - sf(t_reconstructed[:-1])) \ - / dt_reconstructed - - return Q_reconstructed - - def _check_inputs(self): - """ - This method ensures that the instances filled in the Solver object - are what is expected. - - """ - assert isinstance(self.boreholes, (list, Borefield)), \ - "Boreholes must be provided in a list." - assert len(self.boreholes) > 0, \ - "The list of boreholes is empty." - assert np.all([isinstance(b, Borehole) for b in self.boreholes]), \ - "The list of boreholes contains elements that are not Borehole " \ - "objects." - assert self.network is None or isinstance(self.network, Network), \ - "The network is not a valid 'Network' object." - if self.boundary_condition == 'MIFT': - assert not (self.m_flow_network is None and self.m_flow_borehole is None), \ - "The mass flow rate 'm_flow_borehole' or 'm_flow_network' must " \ - "be provided when using the 'MIFT' boundary condition." - assert not (self.m_flow_network is not None and self.m_flow_borehole is not None), \ - "Only one of 'm_flow_borehole' or 'm_flow_network' can " \ - "be provided when using the 'MIFT' boundary condition." - assert not self.cp_f is None, \ - "The heat capacity 'cp_f' must " \ - "be provided when using the 'MIFT' boundary condition." - assert not (type(self.m_flow_borehole) is np.ndarray and not np.size(self.m_flow_borehole, axis=1)==self.network.nInlets), \ - "The number of mass flow rates in 'm_flow_borehole' must " \ - "correspond to the number of circuits in the network." - assert type(self.time) is np.ndarray or isinstance(self.time, (float, np.floating)) or self.time is None, \ - "Time should be a float or an array." - # self.nSegments can now be an int or list - assert type(self.nBoreSegments) is list and len(self.nBoreSegments) == \ - len(self.nBoreSegments) and min(self.nBoreSegments) >= 1, \ - "The argument for number of segments `nSegments` should be " \ - "of type int or a list of integers. If passed as a list, the " \ - "length of the list should be equal to the number of boreholes" \ - "in the borefield. nSegments >= 1 is/are required." - acceptable_boundary_conditions = ['UHTR', 'UBWT', 'MIFT'] - assert type(self.boundary_condition) is str and self.boundary_condition in acceptable_boundary_conditions, \ - f"Boundary condition '{self.boundary_condition}' is not an " \ - f"acceptable boundary condition. \n" \ - f"Please provide one of the following inputs : " \ - f"{acceptable_boundary_conditions}" - assert type(self.approximate_FLS) is bool, \ - "The option 'approximate_FLS' should be set to True or False." - assert type(self.nFLS) is int and 1 <= self.nFLS <= 25, \ - "The option 'nFLS' should be a positive int and lower or equal to " \ - "25." - assert type(self.disp) is bool, \ - "The option 'disp' should be set to True or False." - assert type(self.profiles) is bool, \ - "The option 'profiles' should be set to True or False." - assert type(self.kind) is str, \ - "The option 'kind' should be set to a valid interpolation kind " \ - "in accordance with scipy.interpolate.interp1d options." - acceptable_dtypes = (np.single, np.double) - assert np.any([self.dtype is dtype for dtype in acceptable_dtypes]), \ - f"Data type '{self.dtype}' is not an acceptable data type. \n" \ - f"Please provide one of the following inputs : {acceptable_dtypes}" - # Check segment ratios - for j, (ratios, nSegments) in enumerate( - zip(self.segment_ratios, self.nBoreSegments)): - assert len(ratios) == nSegments, \ - f"The length of the segment ratios vectors must correspond to " \ - f"the number of segments, check borehole {j}." - error = np.abs(1. - np.sum(ratios)) - assert(error < 1.0e-6), \ - f"Defined segment ratios must add up to 1. " \ - f", check borehole {j}." - - return - - -class _Detailed(_BaseSolver): - """ - Detailed solver for the evaluation of the g-function. - - This solver superimposes the finite line source (FLS) solution to - estimate the g-function of a geothermal bore field. Each borehole is - modeled as a series of finite line source segments, as proposed in - [#Detailed-CimBer2014]_. - - Parameters - ---------- - boreholes : list of Borehole objects - List of boreholes included in the bore field. - network : network object - Model of the network. - time : float or array - Values of time (in seconds) for which the g-function is evaluated. - boundary_condition : str - Boundary condition for the evaluation of the g-function. Should be one - of - - - 'UHTR' : - **Uniform heat transfer rate**. This is corresponds to boundary - condition *BC-I* as defined by Cimmino and Bernier (2014) - [#Detailed-CimBer2014]_. - - 'UBWT' : - **Uniform borehole wall temperature**. This is corresponds to - boundary condition *BC-III* as defined by Cimmino and Bernier - (2014) [#Detailed-CimBer2014]_. - - 'MIFT' : - **Mixed inlet fluid temperatures**. This boundary condition was - introduced by Cimmino (2015) [#gFunction-Cimmin2015]_ for - parallel-connected boreholes and extended to mixed - configurations by Cimmino (2019) [#Detailed-Cimmin2019]_. - - nSegments : int or list, optional - Number of line segments used per borehole, or list of number of - line segments used for each borehole. - Default is 8. - segment_ratios : array, list of arrays, or callable, optional - Ratio of the borehole length represented by each segment. The - sum of ratios must be equal to 1. The shape of the array is of - (nSegments,) or list of (nSegments[i],). If segment_ratios==None, - segments of equal lengths are considered. If a callable is provided, it - must return an array of size (nSegments,) when provided with nSegments - (of type int) as an argument, or an array of size (nSegments[i],) when - provided with an element of nSegments (of type list). - Default is :func:`utilities.segment_ratios`. - m_flow_borehole : (nInlets,) array or (nMassFlow, nInlets,) array, optional - Fluid mass flow rate into each circuit of the network. If a - (nMassFlow, nInlets,) array is supplied, the - (nMassFlow, nMassFlow,) variable mass flow rate g-functions - will be evaluated using the method of Cimmino (2024) - [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary - condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be - provided. - Default is None. - m_flow_network : float or (nMassFlow,) array, optional - Fluid mass flow rate into the network of boreholes. If an array - is supplied, the (nMassFlow, nMassFlow,) variable mass flow - rate g-functions will be evaluated using the method of Cimmino - (2024) [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary - condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be - provided. - Default is None. - cp_f : float, optional - Fluid specific isobaric heat capacity (in J/kg.degC). Only required - for the 'MIFT' boundary condition. - Default is None. - approximate_FLS : bool, optional - Set to true to use the approximation of the FLS solution of Cimmino - (2021) [#Detailed-Cimmin2021]_. This approximation does not require the - numerical evaluation of any integral. - Default is False. - nFLS : int, optional - Number of terms in the approximation of the FLS solution. This - parameter is unused if `approximate_FLS` is set to False. - Default is 10. Maximum is 25. - mQuad : int, optional - Number of Gauss-Legendre sample points for the integral over :math:`u` - in the inclined FLS solution. - Default is 11. - linear_threshold : float, optional - Threshold time (in seconds) under which the g-function is - linearized. The g-function value is then interpolated between 0 - and its value at the threshold. If linear_threshold==None, the - g-function is linearized for times - `t < r_b**2 / (25 * self.alpha)`. - Default is None. - disp : bool, optional - Set to true to print progression messages. - Default is False. - profiles : bool, optional - Set to true to keep in memory the temperatures and heat extraction - rates. - Default is False. - kind : string, optional - Interpolation method used for segment-to-segment thermal response - factors. See documentation for scipy.interpolate.interp1d. - Default is 'linear'. - dtype : numpy dtype, optional - numpy data type used for matrices and vectors. Should be one of - numpy.single or numpy.double. - Default is numpy.double. - - References - ---------- - .. [#Detailed-CimBer2014] Cimmino, M., & Bernier, M. (2014). A - semi-analytical method to generate g-functions for geothermal bore - fields. International Journal of Heat and Mass Transfer, 70, 641-650. - .. [#Detailed-Cimmin2019] Cimmino, M. (2019). Semi-analytical method for - g-function calculation of bore fields with series- and - parallel-connected boreholes. Science and Technology for the Built - Environment, 25 (8), 1007-1022. - .. [#Detailed-Cimmin2021] Cimmino, M. (2021). An approximation of the - finite line source solution to model thermal interactions between - geothermal boreholes. International Communications in Heat and Mass - Transfer, 127, 105496. - - """ - def initialize(self, **kwargs): - """ - Split boreholes into segments. - - Returns - ------- - nSources : int - Number of finite line heat sources in the borefield used to - initialize the matrix of segment-to-segment thermal response - factors (of size: nSources x nSources). - - """ - # Split boreholes into segments - self.boreSegments = self.borehole_segments() - nSources = len(self.boreSegments) - return nSources - - def thermal_response_factors(self, time, alpha, kind='linear'): - """ - Evaluate the segment-to-segment thermal response factors for all pairs - of segments in the borefield at all time steps using the finite line - source solution. - - This method returns a scipy.interpolate.interp1d object of the matrix - of thermal response factors, containing a copy of the matrix accessible - by h_ij.y[:nSources,:nSources,:nt+1]. The first index along the - third axis corresponds to time t=0. The interp1d object can be used to - obtain thermal response factors at any intermediate time by - h_ij(t)[:nSources,:nSources]. - - Attributes - ---------- - time : float or array - Values of time (in seconds) for which the g-function is evaluated. - alpha : float - Soil thermal diffusivity (in m2/s). - kind : string, optional - Interpolation method used for segment-to-segment thermal response - factors. See documentation for scipy.interpolate.interp1d. - Default is 'linear'. - - Returns - ------- - h_ij : interp1d - interp1d object (scipy.interpolate) of the matrix of - segment-to-segment thermal response factors. - - """ - if self.disp: - print('Calculating segment to segment response factors ...', - end='') - # Number of time values - nt = len(np.atleast_1d(time)) - # Initialize chrono - tic = perf_counter() - # Initialize segment-to-segment response factors - h_ij = np.zeros((self.nSources, self.nSources, nt+1), dtype=self.dtype) - nBoreholes = len(self.boreholes) - segment_lengths = self.segment_lengths() - - # --------------------------------------------------------------------- - # Segment-to-segment thermal response factors for same-borehole - # thermal interactions - # --------------------------------------------------------------------- - h, i_segment, j_segment = \ - self._thermal_response_factors_borehole_to_self(time, alpha) - # Broadcast values to h_ij matrix - h_ij[j_segment, i_segment, 1:] = h - # --------------------------------------------------------------------- - # Segment-to-segment thermal response factors for - # borehole-to-borehole thermal interactions - # --------------------------------------------------------------------- - for i, (i0, i1) in enumerate(zip(self._i0Segments, self._i1Segments)): - # Segments of the receiving borehole - b2 = self.boreSegments[i0:i1] - if i+1 < nBoreholes: - # Segments of the emitting borehole - b1 = self.boreSegments[i1:] - h = finite_line_source( - time, alpha, b1, b2, approximation=self.approximate_FLS, - N=self.nFLS, M=self.mQuad) - # Broadcast values to h_ij matrix - h_ij[i0:i1, i1:, 1:] = h - h_ij[i1:, i0:i1, 1:] = \ - np.swapaxes(h, 0, 1) * np.divide.outer( - segment_lengths[i0:i1], - segment_lengths[i1:]).T[:,:,np.newaxis] - - # Return 2d array if time is a scalar - if np.isscalar(time): - h_ij = h_ij[:,:,1] - - # Interp1d object for thermal response factors - h_ij = interp1d(np.hstack((0., time)), h_ij, - kind=kind, copy=True, axis=2) - toc = perf_counter() - if self.disp: print(f' {toc - tic:.3f} sec') - - return h_ij - - def _thermal_response_factors_borehole_to_self(self, time, alpha): - """ - Evaluate the segment-to-segment thermal response factors for all pairs - of segments between each borehole and itself. - - Attributes - ---------- - time : float or array - Values of time (in seconds) for which the g-function is evaluated. - alpha : float - Soil thermal diffusivity (in m2/s). - - Returns - ------- - h : array - Finite line source solution. - i_segment : list - Indices of the emitting segments in the bore field. - j_segment : list - Indices of the receiving segments in the bore field. - """ - # Indices of the thermal response factors into h_ij - i_segment = np.concatenate( - [np.repeat(np.arange(i0, i1), nSegments) - for i0, i1, nSegments in zip( - self._i0Segments, self._i1Segments, self.nBoreSegments) - ]) - j_segment = np.concatenate( - [np.tile(np.arange(i0, i1), nSegments) - for i0, i1, nSegments in zip( - self._i0Segments, self._i1Segments, self.nBoreSegments) - ]) - # Unpack parameters - x = np.array([b.x for b in self.boreSegments]) - y = np.array([b.y for b in self.boreSegments]) - H = np.array([b.H for b in self.boreSegments]) - D = np.array([b.D for b in self.boreSegments]) - r_b = np.array([b.r_b for b in self.boreSegments]) - # Distances between boreholes - dis = np.maximum( - np.sqrt((x[i_segment] - x[j_segment])**2 + (y[i_segment] - y[j_segment])**2), - r_b[i_segment]) - # FLS solution - if np.all([b.is_vertical() for b in self.boreholes]): - h = finite_line_source_vectorized( - time, alpha, - dis, H[i_segment], D[i_segment], H[j_segment], D[j_segment], - approximation=self.approximate_FLS, N=self.nFLS) - else: - tilt = np.array([b.tilt for b in self.boreSegments]) - orientation = np.array([b.orientation for b in self.boreSegments]) - h = finite_line_source_inclined_vectorized( - time, alpha, - r_b[i_segment], x[i_segment], y[i_segment], H[i_segment], - D[i_segment], tilt[i_segment], orientation[i_segment], - x[j_segment], y[j_segment], H[j_segment], D[j_segment], - tilt[j_segment], orientation[j_segment], M=self.mQuad, - approximation=self.approximate_FLS, N=self.nFLS) - return h, i_segment, j_segment - - - -class _Similarities(_BaseSolver): - """ - Similarities solver for the evaluation of the g-function. - - This solver superimposes the finite line source (FLS) solution to - estimate the g-function of a geothermal bore field. Each borehole is - modeled as a series of finite line source segments, as proposed in - [#Similarities-CimBer2014]_. The number of evaluations of the FLS solution - is decreased by identifying similar pairs of boreholes, for which the same - FLS value can be applied [#Similarities-Cimmin2018]_. - - Parameters - ---------- - boreholes : list of Borehole objects - List of boreholes included in the bore field. - network : network object - Model of the network. - time : float or array - Values of time (in seconds) for which the g-function is evaluated. - boundary_condition : str - Boundary condition for the evaluation of the g-function. Should be one - of - - - 'UHTR' : - **Uniform heat transfer rate**. This is corresponds to boundary - condition *BC-I* as defined by Cimmino and Bernier (2014) - [#Similarities-CimBer2014]_. - - 'UBWT' : - **Uniform borehole wall temperature**. This is corresponds to - boundary condition *BC-III* as defined by Cimmino and Bernier - (2014) [#Similarities-CimBer2014]_. - - 'MIFT' : - **Mixed inlet fluid temperatures**. This boundary condition was - introduced by Cimmino (2015) [#Similarities-Cimmin2015]_ for - parallel-connected boreholes and extended to mixed - configurations by Cimmino (2019) [#Similarities-Cimmin2019]_. - - nSegments : int or list, optional - Number of line segments used per borehole, or list of number of - line segments used for each borehole. - Default is 8. - segment_ratios : array, list of arrays, or callable, optional - Ratio of the borehole length represented by each segment. The - sum of ratios must be equal to 1. The shape of the array is of - (nSegments,) or list of (nSegments[i],). If segment_ratios==None, - segments of equal lengths are considered. If a callable is provided, it - must return an array of size (nSegments,) when provided with nSegments - (of type int) as an argument, or an array of size (nSegments[i],) when - provided with an element of nSegments (of type list). - Default is :func:`utilities.segment_ratios`. - m_flow_borehole : (nInlets,) array or (nMassFlow, nInlets,) array, optional - Fluid mass flow rate into each circuit of the network. If a - (nMassFlow, nInlets,) array is supplied, the - (nMassFlow, nMassFlow,) variable mass flow rate g-functions - will be evaluated using the method of Cimmino (2024) - [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary - condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be - provided. - Default is None. - m_flow_network : float or (nMassFlow,) array, optional - Fluid mass flow rate into the network of boreholes. If an array - is supplied, the (nMassFlow, nMassFlow,) variable mass flow - rate g-functions will be evaluated using the method of Cimmino - (2024) [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary - condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be - provided. - Default is None. - cp_f : float, optional - Fluid specific isobaric heat capacity (in J/kg.degC). Only required - for the 'MIFT' boundary condition. - Default is None. - approximate_FLS : bool, optional - Set to true to use the approximation of the FLS solution of Cimmino - (2021) [#Similarities-Cimmin2021]_. This approximation does not require - the numerical evaluation of any integral. - Default is False. - nFLS : int, optional - Number of terms in the approximation of the FLS solution. This - parameter is unused if `approximate_FLS` is set to False. - Default is 10. Maximum is 25. - mQuad : int, optional - Number of Gauss-Legendre sample points for the integral over :math:`u` - in the inclined FLS solution. - Default is 11. - linear_threshold : float, optional - Threshold time (in seconds) under which the g-function is - linearized. The g-function value is then interpolated between 0 - and its value at the threshold. If linear_threshold==None, the - g-function is linearized for times - `t < r_b**2 / (25 * self.alpha)`. - Default is None. - disp : bool, optional - Set to true to print progression messages. - Default is False. - profiles : bool, optional - Set to true to keep in memory the temperatures and heat extraction - rates. - Default is False. - kind : string, optional - Interpolation method used for segment-to-segment thermal response - factors. See documentation for scipy.interpolate.interp1d. - Default is 'linear'. - dtype : numpy dtype, optional - numpy data type used for matrices and vectors. Should be one of - numpy.single or numpy.double. - Default is numpy.double. - disTol : float, optional - Relative tolerance on radial distance. Two distances - (d1, d2) between two pairs of boreholes are considered equal if the - difference between the two distances (abs(d1-d2)) is below tolerance. - Default is 0.01. - tol : float, optional - Relative tolerance on length and depth. Two lengths H1, H2 - (or depths D1, D2) are considered equal if abs(H1 - H2)/H2 < tol. - Default is 1.0e-6. - - References - ---------- - .. [#Similarities-CimBer2014] Cimmino, M., & Bernier, M. (2014). A - semi-analytical method to generate g-functions for geothermal bore - fields. International Journal of Heat and Mass Transfer, 70, 641-650. - .. [#Similarities-Cimmin2015] Cimmino, M. (2015). The effects of borehole - thermal resistances and fluid flow rate on the g-functions of geothermal - bore fields. International Journal of Heat and Mass Transfer, 91, - 1119-1127. - .. [#Similarities-Cimmin2018] Cimmino, M. (2018). Fast calculation of the - g-functions of geothermal borehole fields using similarities in the - evaluation of the finite line source solution. Journal of Building - Performance Simulation, 11 (6), 655-668. - .. [#Similarities-Cimmin2019] Cimmino, M. (2019). Semi-analytical method - for g-function calculation of bore fields with series- and - parallel-connected boreholes. Science and Technology for the Built - Environment, 25 (8), 1007-1022. - .. [#Similarities-Cimmin2021] Cimmino, M. (2021). An approximation of the - finite line source solution to model thermal interactions between - geothermal boreholes. International Communications in Heat and Mass - Transfer, 127, 105496. - - """ - def initialize(self, disTol=0.01, tol=1.0e-6, **kwargs): - """ - Split boreholes into segments and identify similarities in the - borefield. - - Returns - ------- - nSources : int - Number of finite line heat sources in the borefield used to - initialize the matrix of segment-to-segment thermal response - factors (of size: nSources x nSources). - - """ - self.disTol = disTol - self.tol = tol - # Check the validity of inputs - self._check_solver_specific_inputs() - # Split boreholes into segments - self.boreSegments = self.borehole_segments() - # Initialize similarities - self.find_similarities() - return len(self.boreSegments) - - def thermal_response_factors(self, time, alpha, kind='linear'): - """ - Evaluate the segment-to-segment thermal response factors for all pairs - of segments in the borefield at all time steps using the finite line - source solution. - - This method returns a scipy.interpolate.interp1d object of the matrix - of thermal response factors, containing a copy of the matrix accessible - by h_ij.y[:nSources,:nSources,:nt+1]. The first index along the - third axis corresponds to time t=0. The interp1d object can be used to - obtain thermal response factors at any intermediat time by - h_ij(t)[:nSources,:nSources]. - - Attributes - ---------- - time : float or array - Values of time (in seconds) for which the g-function is evaluated. - alpha : float - Soil thermal diffusivity (in m2/s). - kind : string, optional - Interpolation method used for segment-to-segment thermal response - factors. See documentation for scipy.interpolate.interp1d. - Default is 'linear'. - - Returns - ------- - h_ij : interp1d - interp1d object (scipy.interpolate) of the matrix of - segment-to-segment thermal response factors. - - """ - if self.disp: - print('Calculating segment to segment response factors ...', - end='') - # Number of time values - nt = len(np.atleast_1d(time)) - # Initialize chrono - tic = perf_counter() - # Initialize segment-to-segment response factors - h_ij = np.zeros((self.nSources, self.nSources, nt+1), dtype=self.dtype) - - # --------------------------------------------------------------------- - # Segment-to-segment thermal response factors for same-borehole thermal - # interactions (vertical boreholes) - # --------------------------------------------------------------------- - # Evaluate FLS at all time steps - h, i_segment, j_segment, k_segment = \ - self._thermal_response_factors_borehole_to_self_vertical( - time, alpha) - # Broadcast values to h_ij matrix - h_ij[j_segment, i_segment, 1:] = h[k_segment, :] - # --------------------------------------------------------------------- - # Segment-to-segment thermal response factors for same-borehole thermal - # interactions (inclined boreholes) - # --------------------------------------------------------------------- - # Evaluate FLS at all time steps - h, i_segment, j_segment, k_segment = \ - self._thermal_response_factors_borehole_to_self_inclined( - time, alpha) - # Broadcast values to h_ij matrix - h_ij[j_segment, i_segment, 1:] = h[k_segment, :] - # --------------------------------------------------------------------- - # Segment-to-segment thermal response factors for borehole-to-borehole - # thermal interactions (vertical boreholes) - # --------------------------------------------------------------------- - for pairs, distances, distance_indices in zip( - self.borehole_to_borehole_vertical, - self.borehole_to_borehole_distances_vertical, - self.borehole_to_borehole_indices_vertical): - # Index of first borehole pair in group - i, j = pairs[0] - # Find segment-to-segment similarities - H1, D1, H2, D2, i_pair, j_pair, k_pair = \ - self._map_axial_segment_pairs_vertical(i, j) - # Locate thermal response factors in the h_ij matrix - i_segment, j_segment, k_segment, l_segment = \ - self._map_segment_pairs_vertical( - i_pair, j_pair, k_pair, pairs, distance_indices) - # Evaluate FLS at all time steps - dis = np.reshape(distances, (-1, 1)) - H1 = H1.reshape(1, -1) - H2 = H2.reshape(1, -1) - D1 = D1.reshape(1, -1) - D2 = D2.reshape(1, -1) - h = finite_line_source_vectorized( - time, alpha, dis, H1, D1, H2, D2, - approximation=self.approximate_FLS, N=self.nFLS) - # Broadcast values to h_ij matrix - h_ij[j_segment, i_segment, 1:] = h[l_segment, k_segment, :] - if (self._compare_boreholes(self.boreholes[j], self.boreholes[i]) and - self.nBoreSegments[i] == self.nBoreSegments[j] and - self._uniform_segment_ratios[i] and - self._uniform_segment_ratios[j]): - h_ij[i_segment, j_segment, 1:] = h[l_segment, k_segment, :] - else: - h_ij[i_segment, j_segment, 1:] = (h * H2.T / H1.T)[l_segment, k_segment, :] - # --------------------------------------------------------------------- - # Segment-to-segment thermal response factors for borehole-to-borehole - # thermal interactions (inclined boreholes) - # --------------------------------------------------------------------- - # Evaluate FLS at all time steps - h, hT, i_segment, j_segment, k_segment = \ - self._thermal_response_factors_borehole_to_borehole_inclined( - time, alpha) - # Broadcast values to h_ij matrix - h_ij[j_segment, i_segment, 1:] = h[k_segment, :] - h_ij[i_segment, j_segment, 1:] = hT[k_segment, :] - - # Return 2d array if time is a scalar - if np.isscalar(time): - h_ij = h_ij[:,:,1] - - # Interp1d object for thermal response factors - h_ij = interp1d( - np.hstack((0., time)), h_ij, - kind=kind, copy=True, assume_sorted=True, axis=2) - toc = perf_counter() - if self.disp: print(f' {toc - tic:.3f} sec') - - return h_ij - - def _thermal_response_factors_borehole_to_borehole_inclined( - self, time, alpha): - """ - Evaluate the segment-to-segment thermal response factors for all pairs - of inclined segments. - - Attributes - ---------- - time : float or array - Values of time (in seconds) for which the g-function is evaluated. - alpha : float - Soil thermal diffusivity (in m2/s). - - Returns - ------- - h : array - Finite line source solution. - hT : array - Reciprocal finite line source solution. - i_segment : list - Indices of the emitting segments in the bore field. - j_segment : list - Indices of the receiving segments in the bore field. - k_segment : list - Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions - corresponding to all pairs in (i_pair, j_pair) in the bore field. - """ - rb1 = np.array([]) - x1 = np.array([]) - y1 = np.array([]) - H1 = np.array([]) - D1 = np.array([]) - tilt1 = np.array([]) - orientation1 = np.array([]) - x2 = np.array([]) - y2 = np.array([]) - H2 = np.array([]) - D2 = np.array([]) - tilt2 = np.array([]) - orientation2 = np.array([]) - i_segment = np.array([], dtype=np.uint) - j_segment = np.array([], dtype=np.uint) - k_segment = np.array([], dtype=np.uint) - k0 = 0 - for pairs in self.borehole_to_borehole_inclined: - # Index of first borehole pair in group - i, j = pairs[0] - # Find segment-to-segment similarities - rb1_i, x1_i, y1_i, H1_i, D1_i, tilt1_i, orientation1_i, \ - x2_i, y2_i, H2_i, D2_i, tilt2_i, orientation2_i, \ - i_pair, j_pair, k_pair = \ - self._map_axial_segment_pairs_inclined(i, j) - # Locate thermal response factors in the h_ij matrix - i_segment_i, j_segment_i, k_segment_i = \ - self._map_segment_pairs_inclined(i_pair, j_pair, k_pair, pairs) - # Append lists - rb1 = np.append(rb1, rb1_i) - x1 = np.append(x1, x1_i) - y1 = np.append(y1, y1_i) - H1 = np.append(H1, H1_i) - D1 = np.append(D1, D1_i) - tilt1 = np.append(tilt1, tilt1_i) - orientation1 = np.append(orientation1, orientation1_i) - x2 = np.append(x2, x2_i) - y2 = np.append(y2, y2_i) - H2 = np.append(H2, H2_i) - D2 = np.append(D2, D2_i) - tilt2 = np.append(tilt2, tilt2_i) - orientation2 = np.append(orientation2, orientation2_i) - i_segment = np.append(i_segment, i_segment_i) - j_segment = np.append(j_segment, j_segment_i) - k_segment = np.append(k_segment, k_segment_i + k0) - k0 += len(k_pair) - # Evaluate FLS at all time steps - h = finite_line_source_inclined_vectorized( - time, alpha, rb1, x1, y1, H1, D1, tilt1, orientation1, - x2, y2, H2, D2, tilt2, orientation2, M=self.mQuad, - approximation=self.approximate_FLS, N=self.nFLS) - hT = (h.T * H2 / H1).T - return h, hT, i_segment, j_segment, k_segment - - def _thermal_response_factors_borehole_to_self_inclined(self, time, alpha): - """ - Evaluate the segment-to-segment thermal response factors for all pairs - of segments between each inclined borehole and itself. - - Attributes - ---------- - time : float or array - Values of time (in seconds) for which the g-function is evaluated. - alpha : float - Soil thermal diffusivity (in m2/s). - - Returns - ------- - h : array - Finite line source solution. - i_segment : list - Indices of the emitting segments in the bore field. - j_segment : list - Indices of the receiving segments in the bore field. - k_segment : list - Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions - corresponding to all pairs in (i_pair, j_pair) in the bore field. - """ - rb1 = np.array([]) - x1 = np.array([]) - y1 = np.array([]) - H1 = np.array([]) - D1 = np.array([]) - tilt1 = np.array([]) - orientation1 = np.array([]) - x2 = np.array([]) - y2 = np.array([]) - H2 = np.array([]) - D2 = np.array([]) - tilt2 = np.array([]) - orientation2 = np.array([]) - i_segment = np.array([], dtype=np.uint) - j_segment = np.array([], dtype=np.uint) - k_segment = np.array([], dtype=np.uint) - k0 = 0 - for group in self.borehole_to_self_inclined: - # Index of first borehole in group - i = group[0] - # Find segment-to-segment similarities - rb1_i, x1_i, y1_i, H1_i, D1_i, tilt1_i, orientation1_i, \ - x2_i, y2_i, H2_i, D2_i, tilt2_i, orientation2_i, \ - i_pair, j_pair, k_pair = \ - self._map_axial_segment_pairs_inclined(i, i) - # Locate thermal response factors in the h_ij matrix - i_segment_i, j_segment_i, k_segment_i = \ - self._map_segment_pairs_inclined( - i_pair, j_pair, k_pair, [(n, n) for n in group]) - # Append lists - rb1 = np.append(rb1, rb1_i) - x1 = np.append(x1, x1_i) - y1 = np.append(y1, y1_i) - H1 = np.append(H1, H1_i) - D1 = np.append(D1, D1_i) - tilt1 = np.append(tilt1, tilt1_i) - orientation1 = np.append(orientation1, orientation1_i) - x2 = np.append(x2, x2_i) - y2 = np.append(y2, y2_i) - H2 = np.append(H2, H2_i) - D2 = np.append(D2, D2_i) - tilt2 = np.append(tilt2, tilt2_i) - orientation2 = np.append(orientation2, orientation2_i) - i_segment = np.append(i_segment, i_segment_i) - j_segment = np.append(j_segment, j_segment_i) - k_segment = np.append(k_segment, k_segment_i + k0) - k0 += len(k_pair) - # Evaluate FLS at all time steps - h = finite_line_source_inclined_vectorized( - time, alpha, rb1, x1, y1, H1, D1, tilt1, orientation1, - x2, y2, H2, D2, tilt2, orientation2, M=self.mQuad, - approximation=self.approximate_FLS, N=self.nFLS) - return h, i_segment, j_segment, k_segment - - def _thermal_response_factors_borehole_to_self_vertical(self, time, alpha): - """ - Evaluate the segment-to-segment thermal response factors for all pairs - of segments between each vertical borehole and itself. - - Attributes - ---------- - time : float or array - Values of time (in seconds) for which the g-function is evaluated. - alpha : float - Soil thermal diffusivity (in m2/s). - - Returns - ------- - h : array - Finite line source solution. - i_segment : list - Indices of the emitting segments in the bore field. - j_segment : list - Indices of the receiving segments in the bore field. - k_segment : list - Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions - corresponding to all pairs in (i_pair, j_pair) in the bore field. - """ - H1 = np.array([]) - D1 = np.array([]) - H2 = np.array([]) - D2 = np.array([]) - dis = np.array([]) - i_segment = np.array([], dtype=np.uint) - j_segment = np.array([], dtype=np.uint) - k_segment = np.array([], dtype=np.uint) - k0 = 0 - for group in self.borehole_to_self_vertical: - # Index of first borehole in group - i = group[0] - # Find segment-to-segment similarities - H1_i, D1_i, H2_i, D2_i, i_pair, j_pair, k_pair = \ - self._map_axial_segment_pairs_vertical(i, i) - # Locate thermal response factors in the h_ij matrix - i_segment_i, j_segment_i, k_segment_i, l_segment_i = \ - self._map_segment_pairs_vertical( - i_pair, j_pair, k_pair, [(n, n) for n in group], [0]) - # Append lists - H1 = np.append(H1, H1_i) - D1 = np.append(D1, D1_i) - H2 = np.append(H2, H2_i) - D2 = np.append(D2, D2_i) - if len(self.borehole_to_self_vertical) > 1: - dis = np.append(dis, np.full(len(H1_i), self.boreholes[i].r_b)) - else: - dis = self.boreholes[i].r_b - i_segment = np.append(i_segment, i_segment_i) - j_segment = np.append(j_segment, j_segment_i) - k_segment = np.append(k_segment, k_segment_i + k0) - k0 += np.max(k_pair) + 1 - # Evaluate FLS at all time steps - h = finite_line_source_vectorized( - time, alpha, dis, H1, D1, H2, D2, - approximation=self.approximate_FLS, N=self.nFLS) - return h, i_segment, j_segment, k_segment - - def find_similarities(self): - """ - Find similarities in the FLS solution for groups of boreholes. - - This function identifies pairs of boreholes for which the evaluation - of the Finite Line Source (FLS) solution is equivalent. - - """ - if self.disp: print('Identifying similarities ...', end='') - # Initialize chrono - tic = perf_counter() - - # Find similar pairs of boreholes - # Boreholes can only be similar if their segments are similar - self.borehole_to_self_vertical, self.borehole_to_self_inclined, \ - self.borehole_to_borehole_vertical, self.borehole_to_borehole_inclined = \ - self._find_axial_borehole_pairs(self.boreholes) - # Find distances for each similar pairs of vertical boreholes - self.borehole_to_borehole_distances_vertical, self.borehole_to_borehole_indices_vertical = \ - self._find_distances( - self.boreholes, self.borehole_to_borehole_vertical) - - # Stop chrono - toc = perf_counter() - if self.disp: print(f' {toc - tic:.3f} sec') - - return - - def _compare_boreholes(self, borehole1, borehole2): - """ - Compare two boreholes and checks if they have the same dimensions : - H, D, r_b, and tilt. - - Parameters - ---------- - borehole1 : Borehole object - First borehole. - borehole2 : Borehole object - Second borehole. - - Returns - ------- - similarity : bool - True if the two boreholes have the same dimensions. - - """ - # Compare lengths (H), buried depth (D) and radius (r_b) - if (abs((borehole1.H - borehole2.H)/borehole1.H) < self.tol and - abs((borehole1.r_b - borehole2.r_b)/borehole1.r_b) < self.tol and - abs((borehole1.D - borehole2.D)/(borehole1.D + 1e-30)) < self.tol and - abs(abs(borehole1.tilt) - abs(borehole2.tilt))/(abs(borehole1.tilt) + 1e-30) < self.tol): - similarity = True - else: - similarity = False - return similarity - - def _compare_real_pairs_vertical(self, pair1, pair2): - """ - Compare two pairs of vertical boreholes or segments and return True if - the two pairs have the same FLS solution for real sources. - - Parameters - ---------- - pair1 : Tuple of Borehole objects - First pair of boreholes or segments. - pair2 : Tuple of Borehole objects - Second pair of boreholes or segments. - - Returns - ------- - similarity : bool - True if the two pairs have the same FLS solution. - - """ - deltaD1 = pair1[1].D - pair1[0].D - deltaD2 = pair2[1].D - pair2[0].D - - # Equality of lengths between pairs - cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol - and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol) - # Equality of lengths in each pair - equal_H = abs((pair1[0].H - pair1[1].H)/pair1[0].H) < self.tol - # Equality of buried depths differences - cond_deltaD = abs(deltaD1 - deltaD2)/abs(deltaD1 + 1e-30) < self.tol - # Equality of buried depths differences if all boreholes have the same - # length - cond_deltaD_equal_H = abs((abs(deltaD1) - abs(deltaD2))/(abs(deltaD1) + 1e-30)) < self.tol - if cond_H and (cond_deltaD or (equal_H and cond_deltaD_equal_H)): - similarity = True - else: - similarity = False - return similarity - - def _compare_image_pairs_vertical(self, pair1, pair2): - """ - Compare two pairs of vertical boreholes or segments and return True if - the two pairs have the same FLS solution for mirror sources. - - Parameters - ---------- - pair1 : Tuple of Borehole objects - First pair of boreholes or segments. - pair2 : Tuple of Borehole objects - Second pair of boreholes or segments. - - Returns - ------- - similarity : bool - True if the two pairs have the same FLS solution. - - """ - sumD1 = pair1[1].D + pair1[0].D - sumD2 = pair2[1].D + pair2[0].D - - # Equality of lengths between pairs - cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol - and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol) - # Equality of buried depths sums - cond_sumD = abs((sumD1 - sumD2)/(sumD1 + 1e-30)) < self.tol - if cond_H and cond_sumD: - similarity = True - else: - similarity = False - return similarity - - def _compare_realandimage_pairs_vertical(self, pair1, pair2): - """ - Compare two pairs of vertical boreholes or segments and return True if - the two pairs have the same FLS solution for both real and mirror - sources. - - Parameters - ---------- - pair1 : Tuple of Borehole objects - First pair of boreholes or segments. - pair2 : Tuple of Borehole objects - Second pair of boreholes or segments. - - Returns - ------- - similarity : bool - True if the two pairs have the same FLS solution. - - """ - if (self._compare_real_pairs_vertical(pair1, pair2) - and self._compare_image_pairs_vertical(pair1, pair2)): - similarity = True - else: - similarity = False - return similarity - - def _compare_real_pairs_inclined(self, pair1, pair2): - """ - Compare two pairs of inclined boreholes or segments and return True if - the two pairs have the same FLS solution for real sources. - - Parameters - ---------- - pair1 : Tuple of Borehole objects - First pair of boreholes or segments. - pair2 : Tuple of Borehole objects - Second pair of boreholes or segments. - - Returns - ------- - similarity : bool - True if the two pairs have the same FLS solution. - - """ - dx1 = pair1[0].x - pair1[1].x; dx2 = pair2[0].x - pair2[1].x - dy1 = pair1[0].y - pair1[1].y; dy2 = pair2[0].y - pair2[1].y - dis1 = np.sqrt(dx1**2 + dy1**2); dis2 = np.sqrt(dx2**2 + dy2**2) - theta_12_1 = np.arctan2(dy1, dx1); theta_12_2 = np.arctan2(dy2, dx2) - deltaD1 = pair1[0].D - pair1[1].D; deltaD2 = pair2[0].D - pair2[1].D - # Equality of lengths between pairs - cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol - and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol) - # Equality of buried depths differences - cond_deltaD = abs(deltaD1 - deltaD2)/(abs(deltaD1) + 1e-30) < self.tol - # Equality of distances - cond_dis = abs(dis1 - dis2)/(abs(dis1) + 1e-30) < self.disTol - # Equality of tilts - cond_beta = ( - abs(abs(pair1[0].tilt) - abs(pair2[0].tilt))/(abs(pair1[0].tilt) + 1e-30) < self.tol - and abs(abs(pair1[1].tilt) - abs(pair2[1].tilt))/(abs(pair1[1].tilt) + 1e-30) < self.tol) - # Equality of relative orientations - sin_b1_cos_dt1_1 = np.sin(pair1[0].tilt) * np.cos(theta_12_1 - pair1[0].orientation) - sin_b2_cos_dt2_1 = np.sin(pair1[1].tilt) * np.cos(theta_12_1 - pair1[1].orientation) - sin_b1_cos_dt1_2 = np.sin(pair2[0].tilt) * np.cos(theta_12_2 - pair2[0].orientation) - sin_b2_cos_dt2_2 = np.sin(pair2[1].tilt) * np.cos(theta_12_2 - pair2[1].orientation) - cond_theta = ( - abs(sin_b1_cos_dt1_1 - sin_b1_cos_dt1_2) / (abs(sin_b1_cos_dt1_1) + 1e-30) > self.tol - and abs(sin_b2_cos_dt2_1 - sin_b2_cos_dt2_2) / (abs(sin_b2_cos_dt2_1) + 1e-30) > self.tol) - if cond_H and cond_deltaD and cond_dis and cond_beta and cond_theta: - similarity = True - else: - similarity = False - return similarity - - def _compare_image_pairs_inclined(self, pair1, pair2): - """ - Compare two pairs of inclined boreholes or segments and return True if - the two pairs have the same FLS solution for mirror sources. - - Parameters - ---------- - pair1 : Tuple of Borehole objects - First pair of boreholes or segments. - pair2 : Tuple of Borehole objects - Second pair of boreholes or segments. - - Returns - ------- - similarity : bool - True if the two pairs have the same FLS solution. - - """ - dx1 = pair1[0].x - pair1[1].x; dx2 = pair2[0].x - pair2[1].x - dy1 = pair1[0].y - pair1[1].y; dy2 = pair2[0].y - pair2[1].y - dis1 = np.sqrt(dx1**2 + dy1**2); dis2 = np.sqrt(dx2**2 + dy2**2) - theta_12_1 = np.arctan2(dy1, dx1); theta_12_2 = np.arctan2(dy2, dx2) - sumD1 = pair1[0].D + pair1[1].D; sumD2 = pair2[0].D + pair2[1].D - # Equality of lengths between pairs - cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol - and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol) - # Equality of buried depths sums - cond_sumD = abs(sumD1 - sumD2)/(abs(sumD1) + 1e-30) < self.tol - # Equality of distances - cond_dis = abs(dis1 - dis2)/(abs(dis1) + 1e-30) < self.disTol - # Equality of tilts - cond_beta = ( - abs(abs(pair1[0].tilt) - abs(pair2[0].tilt))/(abs(pair1[0].tilt) + 1e-30) < self.tol - and abs(abs(pair1[1].tilt) - abs(pair2[1].tilt))/(abs(pair1[1].tilt) + 1e-30) < self.tol) - # Equality of relative orientations - sin_b1_cos_dt1_1 = np.sin(pair1[0].tilt) * np.cos(theta_12_1 - pair1[0].orientation) - sin_b2_cos_dt2_1 = np.sin(pair1[1].tilt) * np.cos(theta_12_1 - pair1[1].orientation) - sin_b1_cos_dt1_2 = np.sin(pair2[0].tilt) * np.cos(theta_12_2 - pair2[0].orientation) - sin_b2_cos_dt2_2 = np.sin(pair2[1].tilt) * np.cos(theta_12_2 - pair2[1].orientation) - cond_theta = ( - abs(sin_b1_cos_dt1_1 - sin_b1_cos_dt1_2) / (abs(sin_b1_cos_dt1_1) + 1e-30) > self.tol - and abs(sin_b2_cos_dt2_1 - sin_b2_cos_dt2_2) / (abs(sin_b2_cos_dt2_1) + 1e-30) > self.tol) - if cond_H and cond_sumD and cond_dis and cond_beta and cond_theta: - similarity = True - else: - similarity = False - return similarity - - def _compare_realandimage_pairs_inclined(self, pair1, pair2): - """ - Compare two pairs of inclined boreholes or segments and return True if - the two pairs have the same FLS solution for both real and mirror - sources. - - Parameters - ---------- - pair1 : Tuple of Borehole objects - First pair of boreholes or segments. - pair2 : Tuple of Borehole objects - Second pair of boreholes or segments. - - Returns - ------- - similarity : bool - True if the two pairs have the same FLS solution. - - Notes - ----- - For inclined boreholes the similarity condition is the same for real - and image parts of the solution. - - """ - dx1 = pair1[0].x - pair1[1].x; dx2 = pair2[0].x - pair2[1].x - dy1 = pair1[0].y - pair1[1].y; dy2 = pair2[0].y - pair2[1].y - dis1 = np.sqrt(dx1**2 + dy1**2); dis2 = np.sqrt(dx2**2 + dy2**2) - theta_12_1 = np.arctan2(dy1, dx1); theta_12_2 = np.arctan2(dy2, dx2) - # Equality of lengths between pairs - cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol - and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol) - # Equality of buried depths - cond_D = ( - abs(pair1[0].D - pair2[0].D)/(abs(pair1[0].D) + 1e-30) < self.tol - and abs(pair1[1].D - pair2[1].D)/(abs(pair1[1].D) + 1e-30) < self.tol) - # Equality of distances - cond_dis = abs(dis1 - dis2)/(abs(dis1) + 1e-30) < self.disTol - # Equality of tilts - cond_beta = ( - abs(abs(pair1[0].tilt) - abs(pair2[0].tilt))/(abs(pair1[0].tilt) + 1e-30) < self.tol - and abs(abs(pair1[1].tilt) - abs(pair2[1].tilt))/(abs(pair1[1].tilt) + 1e-30) < self.tol) - # Equality of relative orientations - sin_b1_cos_dt1_1 = np.sin(pair1[0].tilt) * np.cos(theta_12_1 - pair1[0].orientation) - sin_b2_cos_dt2_1 = np.sin(pair1[1].tilt) * np.cos(theta_12_1 - pair1[1].orientation) - sin_b1_cos_dt1_2 = np.sin(pair2[0].tilt) * np.cos(theta_12_2 - pair2[0].orientation) - sin_b2_cos_dt2_2 = np.sin(pair2[1].tilt) * np.cos(theta_12_2 - pair2[1].orientation) - cond_theta = ( - abs(sin_b1_cos_dt1_1 - sin_b1_cos_dt1_2) / (abs(sin_b1_cos_dt1_1) + 1e-30) < self.tol - and abs(sin_b2_cos_dt2_1 - sin_b2_cos_dt2_2) / (abs(sin_b2_cos_dt2_1) + 1e-30) < self.tol) - if cond_H and cond_D and cond_dis and cond_beta and cond_theta: - similarity = True - else: - similarity = False - return similarity - - def _find_axial_borehole_pairs(self, boreholes): - """ - Find axial (i.e. disregarding the radial distance) similarities between - borehole pairs to simplify the evaluation of the FLS solution. - - Parameters - ---------- - boreholes : list of Borehole objects - Boreholes in the bore field. - - Returns - ------- - borehole_to_self : list - Lists of borehole indexes for each unique set of borehole - dimensions (H, D, r_b) in the bore field. - borehole_to_borehole : list - Lists of tuples of borehole indexes for each unique pair of - boreholes that share the same (pairwise) dimensions (H, D). - - """ - nBoreholes = len(boreholes) - borehole_to_self_vertical = [] - borehole_to_self_inclined = [] - # Only check for similarities if there is more than one borehole - if nBoreholes > 1: - borehole_to_borehole_vertical = [] - borehole_to_borehole_inclined = [] - for i, (borehole_i, nSegments_i, ratios_i) in enumerate( - zip(boreholes, self.nBoreSegments, self.segment_ratios)): - # Compare the borehole to all known unique sets of dimensions - if borehole_i.is_vertical(): - borehole_to_self = borehole_to_self_vertical - compare_pairs = self._compare_realandimage_pairs_vertical - else: - borehole_to_self = borehole_to_self_inclined - compare_pairs = self._compare_realandimage_pairs_inclined - for k, borehole_set in enumerate(borehole_to_self): - m = borehole_set[0] - # Add the borehole to the group if a similar borehole is - # found - if (self._compare_boreholes(borehole_i, boreholes[m]) and - (self._equal_segment_ratios or - (nSegments_i == self.nBoreSegments[m] and - np.allclose(ratios_i, - self.segment_ratios[m], - rtol=self.tol)))): - borehole_set.append(i) - break - else: - # If no similar boreholes are known, append the groups - borehole_to_self.append([i]) - - for j, (borehole_j, nSegments_j, ratios_j) in enumerate( - zip(boreholes[i+1:], - self.nBoreSegments[i+1:], - self.segment_ratios[i+1:]), - start=i+1): - pair0 = (borehole_i, borehole_j) # pair - pair1 = (borehole_j, borehole_i) # reciprocal pair - # Compare pairs of boreholes to known unique pairs - if borehole_i.is_vertical() and borehole_j.is_vertical(): - borehole_to_borehole = borehole_to_borehole_vertical - compare_pairs = self._compare_realandimage_pairs_vertical - else: - borehole_to_borehole = borehole_to_borehole_inclined - compare_pairs = self._compare_realandimage_pairs_inclined - for pairs in borehole_to_borehole: - m, n = pairs[0] - pair_ref = (boreholes[m], boreholes[n]) - # Add the pair (or the reciprocal pair) to a group - # if a similar one is found - if (compare_pairs(pair0, pair_ref) and - (self._equal_segment_ratios or - (nSegments_i == self.nBoreSegments[m] and - nSegments_j == self.nBoreSegments[n] and - np.allclose(ratios_i, - self.segment_ratios[m], - rtol=self.tol) and - np.allclose(ratios_j, - self.segment_ratios[n], - rtol=self.tol)))): - pairs.append((i, j)) - break - elif (compare_pairs(pair1, pair_ref) and - (self._equal_segment_ratios or - (nSegments_j == self.nBoreSegments[m] and - nSegments_i == self.nBoreSegments[n] and - np.allclose(ratios_j, - self.segment_ratios[m], - rtol=self.tol) and - np.allclose(ratios_i, - self.segment_ratios[n], - rtol=self.tol)))): - pairs.append((j, i)) - break - # If no similar pairs are known, append the groups - else: - borehole_to_borehole.append([(i, j)]) - - else: - # Outputs for a single borehole - if boreholes[0].is_vertical: - borehole_to_self_vertical = [[0]] - borehole_to_self_inclined = [] - else: - borehole_to_self_vertical = [] - borehole_to_self_inclined = [[0]] - borehole_to_borehole_vertical = [] - borehole_to_borehole_inclined = [] - return borehole_to_self_vertical, borehole_to_self_inclined, \ - borehole_to_borehole_vertical, borehole_to_borehole_inclined - - def _find_distances(self, boreholes, borehole_to_borehole): - """ - Find unique distances between pairs of boreholes for each unique pair - of boreholes in the bore field. - - Parameters - ---------- - boreholes : list of Borehole objects - Boreholes in the bore field. - borehole_to_borehole : list - Lists of tuples of borehole indexes for each unique pair of - boreholes that share the same (pairwise) dimensions (H, D). - - Returns - ------- - borehole_to_borehole_distances : list - Sorted lists of borehole-to-borehole radial distances for each - unique pair of boreholes. - borehole_to_borehole_indices : list - Lists of indexes of distances associated with each borehole pair. - - """ - nGroups = len(borehole_to_borehole) - borehole_to_borehole_distances = [[] for i in range(nGroups)] - borehole_to_borehole_indices = \ - [np.empty(len(group), dtype=np.uint) for group in borehole_to_borehole] - # Find unique distances for each group - for i, (pairs, distances, distance_indices) in enumerate( - zip(borehole_to_borehole, - borehole_to_borehole_distances, - borehole_to_borehole_indices)): - nPairs = len(pairs) - # Array of all borehole-to-borehole distances within the group - all_distances = np.array( - [boreholes[pair[0]].distance(boreholes[pair[1]]) - for pair in pairs]) - # Indices to sort the distance array - i_sort = all_distances.argsort() - # Sort the distance array - distances_sorted = all_distances[i_sort] - j0 = 0 - j1 = 1 - nDis = 0 - # For each increasing distance in the sorted array : - # 1 - find all distances that are within tolerance - # 2 - add the average distance in the list of unique distances - # 3 - associate the distance index to all pairs for the identified - # distances - # 4 - re-start at the next distance index not yet accounted for. - while j0 < nPairs and j1 > 0: - # Find the first distance outside tolerance - j1 = np.argmax( - distances_sorted >= (1+self.disTol)*distances_sorted[j0]) - if j1 > j0: - # Average distance between pairs of boreholes - distances.append(np.mean(distances_sorted[j0:j1])) - # Apply distance index to borehole pairs - distance_indices[i_sort[j0:j1]] = nDis - else: - # Average distance between pairs of boreholes - distances.append(np.mean(distances_sorted[j0:])) - # Apply distance index to borehole pairs - distance_indices[i_sort[j0:]] = nDis - j0 = j1 - nDis += 1 - return borehole_to_borehole_distances, borehole_to_borehole_indices - - def _map_axial_segment_pairs_vertical( - self, i, j, reaSource=True, imgSource=True): - """ - Find axial (i.e. disregarding the radial distance) similarities between - segment pairs along two boreholes to simplify the evaluation of the - FLS solution. - - The returned H1, D1, H2, and D2 can be used to evaluate the segment-to- - segment response factors using scipy.integrate.quad_vec. - - Parameters - ---------- - i : int - Index of the first borehole. - j : int - Index of the second borehole. - - Returns - ------- - H1 : array - Length of the emitting segments. - D1 : array - Array of buried depths of the emitting segments. - H2 : array - Length of the receiving segments. - D2 : array - Array of buried depths of the receiving segments. - i_pair : list - Indices of the emitting segments along a borehole. - j_pair : list - Indices of the receiving segments along a borehole. - k_pair : list - Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions - corresponding to all pairs in (i_pair, j_pair). - - """ - # Initialize local variables - borehole1 = self.boreholes[i] - borehole2 = self.boreholes[j] - assert reaSource or imgSource, \ - "At least one of reaSource and imgSource must be True." - if reaSource and imgSource: - # Find segment pairs for the full (real + image) FLS solution - compare_pairs = self._compare_realandimage_pairs_vertical - elif reaSource: - # Find segment pairs for the real FLS solution - compare_pairs = self._compare_real_pairs_vertical - elif imgSource: - # Find segment pairs for the image FLS solution - compare_pairs = self._compare_image_pairs_vertical - # Dive both boreholes into segments - segments1 = borehole1.segments( - self.nBoreSegments[i], segment_ratios=self.segment_ratios[i]) - segments2 = borehole2.segments( - self.nBoreSegments[j], segment_ratios=self.segment_ratios[j]) - # Prepare lists of segment lengths - H1 = [] - H2 = [] - # Prepare lists of segment buried depths - D1 = [] - D2 = [] - # All possible pairs (i, j) of indices between segments - i_pair = np.repeat(np.arange(self.nBoreSegments[i], dtype=np.uint), - self.nBoreSegments[j]) - j_pair = np.tile(np.arange(self.nBoreSegments[j], dtype=np.uint), - self.nBoreSegments[i]) - # Empty list of indices for unique pairs - k_pair = np.empty(self.nBoreSegments[i] * self.nBoreSegments[j], - dtype=np.uint) - unique_pairs = [] - nPairs = 0 - - p = 0 - for ii, segment_i in enumerate(segments1): - for jj, segment_j in enumerate(segments2): - pair = (segment_i, segment_j) - # Compare the segment pairs to all known unique pairs - for k, pair_k in enumerate(unique_pairs): - m, n = pair_k[0], pair_k[1] - pair_ref = (segments1[m], segments2[n]) - # Stop if a similar pair is found and assign the index - if compare_pairs(pair, pair_ref): - k_pair[p] = k - break - # If no similar pair is found : add a new pair, increment the - # number of unique pairs, and extract the associated buried - # depths - else: - k_pair[p] = nPairs - H1.append(segment_i.H) - H2.append(segment_j.H) - D1.append(segment_i.D) - D2.append(segment_j.D) - unique_pairs.append((ii, jj)) - nPairs += 1 - p += 1 - return np.array(H1), np.array(D1), np.array(H2), np.array(D2), i_pair, j_pair, k_pair - - def _map_axial_segment_pairs_inclined( - self, i, j, reaSource=True, imgSource=True): - """ - Find axial similarities between segment pairs along two boreholes to - simplify the evaluation of the FLS solution. - - The returned H1, D1, H2, and D2 can be used to evaluate the segment-to- - segment response factors using scipy.integrate.quad_vec. - - Parameters - ---------- - i : int - Index of the first borehole. - j : int - Index of the second borehole. - - Returns - ------- - rb1 : array - Radii of the emitting heat sources. - x1 : array - x-Positions of the emitting heat sources. - y1 : array - y-Positions of the emitting heat sources. - H1 : array - Lengths of the emitting heat sources. - D1 : array - Buried depths of the emitting heat sources. - tilt1 : array - Angles (in radians) from vertical of the emitting heat sources. - orientation1 : array - Directions (in radians) of the tilt the emitting heat sources. - x2 : array - x-Positions of the receiving heat sources. - y2 : array - y-Positions of the receiving heat sources. - H2 : array - Lengths of the receiving heat sources. - D2 : array - Buried depths of the receiving heat sources. - tilt2 : array - Angles (in radians) from vertical of the receiving heat sources. - orientation2 : array - Directions (in radians) of the tilt the receiving heat sources. - i_pair : list - Indices of the emitting segments along a borehole. - j_pair : list - Indices of the receiving segments along a borehole. - k_pair : list - Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions - corresponding to all pairs in (i_pair, j_pair). - """ - # Initialize local variables - borehole1 = self.boreholes[i] - borehole2 = self.boreholes[j] - assert reaSource or imgSource, \ - "At least one of reaSource and imgSource must be True." - if reaSource and imgSource: - # Find segment pairs for the full (real + image) FLS solution - compare_pairs = self._compare_realandimage_pairs_inclined - elif reaSource: - # Find segment pairs for the real FLS solution - compare_pairs = self._compare_real_pairs_inclined - elif imgSource: - # Find segment pairs for the image FLS solution - compare_pairs = self._compare_image_pairs_inclined - # Dive both boreholes into segments - segments1 = borehole1.segments( - self.nBoreSegments[i], segment_ratios=self.segment_ratios[i]) - segments2 = borehole2.segments( - self.nBoreSegments[j], segment_ratios=self.segment_ratios[j]) - # Prepare lists of FLS-inclined arguments - rb1 = [] - x1 = [] - y1 = [] - H1 = [] - D1 = [] - tilt1 = [] - orientation1 = [] - x2 = [] - y2 = [] - H2 = [] - D2 = [] - tilt2 = [] - orientation2 = [] - # All possible pairs (i, j) of indices between segments - i_pair = np.repeat(np.arange(self.nBoreSegments[i], dtype=np.uint), - self.nBoreSegments[j]) - j_pair = np.tile(np.arange(self.nBoreSegments[j], dtype=np.uint), - self.nBoreSegments[i]) - # Empty list of indices for unique pairs - k_pair = np.empty(self.nBoreSegments[i] * self.nBoreSegments[j], - dtype=np.uint) - unique_pairs = [] - nPairs = 0 - - p = 0 - for ii, segment_i in enumerate(segments1): - for jj, segment_j in enumerate(segments2): - pair = (segment_i, segment_j) - # Compare the segment pairs to all known unique pairs - for k, pair_k in enumerate(unique_pairs): - m, n = pair_k[0], pair_k[1] - pair_ref = (segments1[m], segments2[n]) - # Stop if a similar pair is found and assign the index - if compare_pairs(pair, pair_ref): - k_pair[p] = k - break - # If no similar pair is found : add a new pair, increment the - # number of unique pairs, and extract the associated buried - # depths - else: - k_pair[p] = nPairs - rb1.append(segment_i.r_b) - x1.append(segment_i.x) - y1.append(segment_i.y) - H1.append(segment_i.H) - D1.append(segment_i.D) - tilt1.append(segment_i.tilt) - orientation1.append(segment_i.orientation) - x2.append(segment_j.x) - y2.append(segment_j.y) - H2.append(segment_j.H) - D2.append(segment_j.D) - tilt2.append(segment_j.tilt) - orientation2.append(segment_j.orientation) - unique_pairs.append((ii, jj)) - nPairs += 1 - p += 1 - return np.array(rb1), np.array(x1), np.array(y1), np.array(H1), \ - np.array(D1), np.array(tilt1), np.array(orientation1), \ - np.array(x2), np.array(y2), np.array(H2), np.array(D2), \ - np.array(tilt2), np.array(orientation2), i_pair, j_pair, k_pair - - def _map_segment_pairs_vertical( - self, i_pair, j_pair, k_pair, borehole_to_borehole, - borehole_to_borehole_indices): - """ - Return the maping of the unique segment-to-segment thermal response - factors (h) to the complete h_ij array of the borefield, such that: - - h_ij[j_segment, i_segment, :nt] = h[:nt, l_segment, k_segment].T, - - where h is the array of unique segment-to-segment thermal response - factors for a given unique pair of boreholes at all unique distances. - - Parameters - ---------- - i_pair : list - Indices of the emitting segments. - j_pair : list - Indices of the receiving segments. - k_pair : list - Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions - corresponding to all pairs in (i_pair, j_pair). - borehole_to_borehole : list - Tuples of borehole indexes. - borehole_to_borehole_indices : list - Indexes of distances. - - Returns - ------- - i_segment : list - Indices of the emitting segments in the bore field. - j_segment : list - Indices of the receiving segments in the bore field. - k_segment : list - Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions - corresponding to all pairs in (i_pair, j_pair) in the bore field. - l_segment : list - Indices of unique distances for all pairs in (i_pair, j_pair) - in the bore field. - - """ - i_segment = np.concatenate( - [i_pair + self._i0Segments[i] for (i, j) in borehole_to_borehole]) - j_segment = np.concatenate( - [j_pair + self._i0Segments[j] for (i, j) in borehole_to_borehole]) - k_segment = np.tile(k_pair, len(borehole_to_borehole)) - l_segment = np.concatenate( - [np.repeat(i, len(k_pair)) for i in borehole_to_borehole_indices]) - return i_segment, j_segment, k_segment, l_segment - - def _map_segment_pairs_inclined( - self, i_pair, j_pair, k_pair, borehole_to_borehole): - """ - Return the maping of the unique segment-to-segment thermal response - factors (h) to the complete h_ij array of the borefield, such that: - - h_ij[j_segment, i_segment, :nt] = h[:nt, k_segment].T, - - where h is the array of unique segment-to-segment thermal response - factors for a given unique pair of boreholes at all unique distances. - - Parameters - ---------- - i_pair : list - Indices of the emitting segments. - j_pair : list - Indices of the receiving segments. - k_pair : list - Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions - corresponding to all pairs in (i_pair, j_pair). - borehole_to_borehole : list - Tuples of borehole indexes. - - Returns - ------- - i_segment : list - Indices of the emitting segments in the bore field. - j_segment : list - Indices of the receiving segments in the bore field. - k_segment : list - Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions - corresponding to all pairs in (i_pair, j_pair) in the bore field. - - """ - i_segment = np.concatenate( - [i_pair + self._i0Segments[i] for (i, j) in borehole_to_borehole]) - j_segment = np.concatenate( - [j_pair + self._i0Segments[j] for (i, j) in borehole_to_borehole]) - k_segment = np.tile(k_pair, len(borehole_to_borehole)) - return i_segment, j_segment, k_segment - - def _check_solver_specific_inputs(self): - """ - This method ensures that solver specific inputs to the Solver object - are what is expected. - - """ - assert isinstance(self.disTol, (np.floating, float)) and self.disTol > 0., \ - "The distance tolerance 'disTol' should be a positive float." - assert isinstance(self.tol, (np.floating, float)) and self.tol > 0., \ - "The relative tolerance 'tol' should be a positive float." - return - - -class _Equivalent(_BaseSolver): - """ - Equivalent solver for the evaluation of the g-function. - - This solver uses hierarchical agglomerative clustering to identify groups - of boreholes that are expected to have similar borehole wall temperatures - and heat extraction rates, as proposed by Prieto and Cimmino (2021) - [#Equivalent-PriCim2021]_. Each group of boreholes is represented by a - single equivalent borehole. The FLS solution is adapted to evaluate - thermal interactions between groups of boreholes. This greatly reduces - the number of evaluations of the FLS solution and the size of the system of - equations to evaluate the g-function. - - Parameters - ---------- - boreholes : list of Borehole objects - List of boreholes included in the bore field. - network : network object - Model of the network. - time : float or array - Values of time (in seconds) for which the g-function is evaluated. - boundary_condition : str - Boundary condition for the evaluation of the g-function. Should be one - of - - - 'UHTR' : - **Uniform heat transfer rate**. This is corresponds to boundary - condition *BC-I* as defined by Cimmino and Bernier (2014) - [#Equivalent-CimBer2014]_. - - 'UBWT' : - **Uniform borehole wall temperature**. This is corresponds to - boundary condition *BC-III* as defined by Cimmino and Bernier - (2014) [#Equivalent-CimBer2014]_. - - 'MIFT' : - **Mixed inlet fluid temperatures**. This boundary condition was - introduced by Cimmino (2015) [#Equivalent-Cimmin2015]_ for - parallel-connected boreholes and extended to mixed - configurations by Cimmino (2019) [#Equivalent-Cimmin2019]_. - - nSegments : int or list, optional - Number of line segments used per borehole, or list of number of - line segments used for each borehole. - Default is 8. - segment_ratios : array, list of arrays, or callable, optional - Ratio of the borehole length represented by each segment. The - sum of ratios must be equal to 1. The shape of the array is of - (nSegments,) or list of (nSegments[i],). If segment_ratios==None, - segments of equal lengths are considered. If a callable is provided, it - must return an array of size (nSegments,) when provided with nSegments - (of type int) as an argument, or an array of size (nSegments[i],) when - provided with an element of nSegments (of type list). - Default is :func:`utilities.segment_ratios`. - m_flow_borehole : (nInlets,) array or (nMassFlow, nInlets,) array, optional - Fluid mass flow rate into each circuit of the network. If a - (nMassFlow, nInlets,) array is supplied, the - (nMassFlow, nMassFlow,) variable mass flow rate g-functions - will be evaluated using the method of Cimmino (2024) - [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary - condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be - provided. - Default is None. - m_flow_network : float or (nMassFlow,) array, optional - Fluid mass flow rate into the network of boreholes. If an array - is supplied, the (nMassFlow, nMassFlow,) variable mass flow - rate g-functions will be evaluated using the method of Cimmino - (2024) [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary - condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be - provided. - Default is None. - cp_f : float, optional - Fluid specific isobaric heat capacity (in J/kg.degC). Only required - for the 'MIFT' boundary condition. - Default is None. - approximate_FLS : bool, optional - Set to true to use the approximation of the FLS solution of Cimmino - (2021) [#Equivalent-Cimmin2021]_. This approximation does not require - the numerical evaluation of any integral. When using the 'equivalent' - solver, the approximation is only applied to the thermal response at - the borehole radius. Thermal interaction between boreholes is evaluated - using the FLS solution. - Default is False. - nFLS : int, optional - Number of terms in the approximation of the FLS solution. This - parameter is unused if `approximate_FLS` is set to False. - Default is 10. Maximum is 25. - mQuad : int, optional - Number of Gauss-Legendre sample points for the integral over :math:`u` - in the inclined FLS solution. - Default is 11. - linear_threshold : float, optional - Threshold time (in seconds) under which the g-function is - linearized. The g-function value is then interpolated between 0 - and its value at the threshold. If linear_threshold==None, the - g-function is linearized for times - `t < r_b**2 / (25 * self.alpha)`. - Default is None. - disp : bool, optional - Set to true to print progression messages. - Default is False. - profiles : bool, optional - Set to true to keep in memory the temperatures and heat extraction - rates. - Default is False. - kind : string, optional - Interpolation method used for segment-to-segment thermal response - factors. See documentation for scipy.interpolate.interp1d. - Default is 'linear'. - dtype : numpy dtype, optional - numpy data type used for matrices and vectors. Should be one of - numpy.single or numpy.double. - Default is numpy.double. - disTol : float, optional - Relative tolerance on radial distance. Two distances - (d1, d2) between two pairs of boreholes are considered equal if the - difference between the two distances (abs(d1-d2)) is below tolerance. - Default is 0.01. - tol : float, optional - Relative tolerance on length and depth. Two lengths H1, H2 - (or depths D1, D2) are considered equal if abs(H1 - H2)/H2 < tol. - Default is 1.0e-6. - kClusters : int, optional - Increment on the minimum number of equivalent boreholes determined by - cutting the dendrogram of the bore field given by the hierarchical - agglomerative clustering method. Increasing the value of this parameter - increases the accuracy of the method. - Default is 1. - - References - ---------- - .. [#Equivalent-CimBer2014] Cimmino, M., & Bernier, M. (2014). A - semi-analytical method to generate g-functions for geothermal bore - fields. International Journal of Heat and Mass Transfer, 70, 641-650. - .. [#Equivalent-Cimmin2015] Cimmino, M. (2015). The effects of borehole - thermal resistances and fluid flow rate on the g-functions of geothermal - bore fields. International Journal of Heat and Mass Transfer, 91, - 1119-1127. - .. [#Equivalent-Cimmin2018] Cimmino, M. (2018). Fast calculation of the - g-functions of geothermal borehole fields using similarities in the - evaluation of the finite line source solution. Journal of Building - Performance Simulation, 11 (6), 655-668. - .. [#Equivalent-PriCim2021] Prieto, C., & Cimmino, M. - (2021). Thermal interactions in large irregular fields of geothermal - boreholes: the method of equivalent borehole. Journal of Building - Performance Simulation, 14 (4), 446-460. - .. [#Equivalent-Cimmin2021] Cimmino, M. (2021). An approximation of the - finite line source solution to model thermal interactions between - geothermal boreholes. International Communications in Heat and Mass - Transfer, 127, 105496. - - """ - def initialize(self, disTol=0.01, tol=1.0e-6, kClusters=1, **kwargs): - """ - Initialize paramteters. Identify groups for equivalent boreholes. - - Returns - ------- - nSources : int - Number of finite line heat sources in the borefield used to - initialize the matrix of segment-to-segment thermal response - factors (of size: nSources x nSources). - - """ - self.disTol = disTol - self.tol = tol - self.kClusters = kClusters - # Check the validity of inputs - self._check_solver_specific_inputs() - # Initialize groups for equivalent boreholes - nSources = self.find_groups() - self.nBoreSegments = [self.nBoreSegments[0]] * self.nEqBoreholes - self.segment_ratios = [self.segment_ratios[0]] * self.nEqBoreholes - self.boreSegments = self.borehole_segments() - self._i0Segments = [sum(self.nBoreSegments[0:i]) - for i in range(self.nEqBoreholes)] - self._i1Segments = [sum(self.nBoreSegments[0:(i + 1)]) - for i in range(self.nEqBoreholes)] - return nSources - - def thermal_response_factors(self, time, alpha, kind='linear'): - """ - Evaluate the segment-to-segment thermal response factors for all pairs - of segments in the borefield at all time steps using the finite line - source solution. - - This method returns a scipy.interpolate.interp1d object of the matrix - of thermal response factors, containing a copy of the matrix accessible - by h_ij.y[:nSources,:nSources,:nt+1]. The first index along the - third axis corresponds to time t=0. The interp1d object can be used to - obtain thermal response factors at any intermediate time by - h_ij(t)[:nSources,:nSources]. - - Parameters - ---------- - time : float or array - Values of time (in seconds) for which the g-function is evaluated. - alpha : float - Soil thermal diffusivity (in m2/s). - kind : string, optional - Interpolation method used for segment-to-segment thermal response - factors. See documentation for scipy.interpolate.interp1d. - Default is linear. - - Returns - ------- - h_ij : interp1d - interp1d object (scipy.interpolate) of the matrix of - segment-to-segment thermal response factors. - - """ - if self.disp: - print('Calculating segment to segment response factors ...', - end='') - # Number of time values - nt = len(np.atleast_1d(time)) - # Initialize chrono - tic = perf_counter() - # Initialize segment-to-segment response factors - h_ij = np.zeros((self.nSources, self.nSources, nt+1), dtype=self.dtype) - segment_lengths = self.segment_lengths() - - # --------------------------------------------------------------------- - # Segment-to-segment thermal response factors for borehole-to-borehole - # thermal interactions - # --------------------------------------------------------------------- - # Groups correspond to unique pairs of borehole dimensions - for pairs in self.borehole_to_borehole: - i, j = pairs[0] - # Prepare inputs to the FLS function - dis, wDis = self._find_unique_distances(self.dis, pairs) - H1, D1, H2, D2, i_pair, j_pair, k_pair = \ - self._map_axial_segment_pairs(i, j) - H1 = H1.reshape(1, -1) - H2 = H2.reshape(1, -1) - D1 = D1.reshape(1, -1) - D2 = D2.reshape(1, -1) - N2 = np.array( - [[self.boreholes[j].nBoreholes for (i, j) in pairs]]).T - # Evaluate FLS at all time steps - h = finite_line_source_equivalent_boreholes_vectorized( - time, alpha, dis, wDis, H1, D1, H2, D2, N2) - # Broadcast values to h_ij matrix - for k, (i, j) in enumerate(pairs): - i_segment = self._i0Segments[i] + i_pair - j_segment = self._i0Segments[j] + j_pair - h_ij[j_segment, i_segment, 1:] = h[k, k_pair, :] - if not i == j: - h_ij[i_segment, j_segment, 1:] = (h[k, k_pair, :].T \ - * segment_lengths[j_segment]/segment_lengths[i_segment]).T - - # --------------------------------------------------------------------- - # Segment-to-segment thermal response factors for same-borehole thermal - # interactions - # --------------------------------------------------------------------- - # Groups correspond to unique borehole dimensions - for group in self.borehole_to_self: - # Index of first borehole in group - i = group[0] - # Find segment-to-segment similarities - H1, D1, H2, D2, i_pair, j_pair, k_pair = \ - self._map_axial_segment_pairs(i, i) - # Evaluate FLS at all time steps - dis = self.boreholes[i].r_b - H1 = H1.reshape(1, -1) - H2 = H2.reshape(1, -1) - D1 = D1.reshape(1, -1) - D2 = D2.reshape(1, -1) - h = finite_line_source_vectorized( - time, alpha, dis, H1, D1, H2, D2, - approximation=self.approximate_FLS, N=self.nFLS) - # Broadcast values to h_ij matrix - for i in group: - i_segment = self._i0Segments[i] + i_pair - j_segment = self._i0Segments[i] + j_pair - h_ij[j_segment, i_segment, 1:] = \ - h_ij[j_segment, i_segment, 1:] + h[0, k_pair, :] - - # Return 2d array if time is a scalar - if np.isscalar(time): - h_ij = h_ij[:,:,1] - - # Interp1d object for thermal response factors - h_ij = interp1d(np.hstack((0., time)), h_ij, - kind=kind, copy=True, axis=2) - toc = perf_counter() - if self.disp: print(f' {toc - tic:.3f} sec') - - return h_ij - - def find_groups(self, tol=1e-6): - """ - Identify groups of boreholes that can be represented by a single - equivalent borehole for the calculation of the g-function. - - Hierarchical agglomerative clustering is applied to the superposed - steady-state finite line source solution (i.e. the steady-state - dimensionless borehole wall temperature due to a uniform heat - extraction equal for all boreholes). The number of clusters is - evaluated by cutting the dendrogram at the half-height of the longest - branch and incrementing the number of intercepted branches by the value - of the kClusters parameter. - - Parameters - ---------- - tol : float - Tolerance on the temperature to identify the maxiumum number of - equivalent boreholes. - Default is 1e-6. - - Returns - ------- - nSources : int - Number of heat sources in the bore field. - - """ - if self.disp: print('Identifying equivalent boreholes ...', end='') - # Initialize chrono - tic = perf_counter() - - # Temperature change of individual boreholes - self.nBoreholes = len(self.boreholes) - # Equivalent field formed by all boreholes - eqField = _EquivalentBorehole(self.boreholes) - if self.nBoreholes > 1: - # Spatial superposition of the steady-state FLS solution - data = np.sum(finite_line_source(np.inf, 1., self.boreholes, self.boreholes), axis=1).reshape(-1,1) - # Split boreholes into groups of same dimensions - unique_boreholes = self._find_unique_boreholes(self.boreholes) - # Initialize empty list of clusters - self.clusters = [] - self.nEqBoreholes = 0 - for group in unique_boreholes: - if len(group) > 1: - # Maximum temperature - maxTemp = np.max(data[group]) - # Hierarchical agglomerative clustering based on temperatures - clusterization = linkage(data[group], method='complete') - dcoord = np.array( - dendrogram(clusterization, no_plot=True)['dcoord']) - # Maximum number of clusters - # Height to cut each tree to obtain the minimum number of clusters - disLeft = dcoord[:,1] - dcoord[:,0] - disRight = dcoord[:,2] - dcoord[:,3] - if np.max(disLeft) >= np.max(disRight): - i = disLeft.argmax() - height = 0.5*(dcoord[i,1] + dcoord[i,0]) - else: - i = disRight.argmax() - height = 0.5*(dcoord[i,2] + dcoord[i,3]) - # Find the number of clusters and increment by kClusters - # Maximum number of clusters - nClustersMax = min(np.sum(dcoord[:,1] > tol*maxTemp) + 1, - len(group)) - # Optimal number of cluster - nClusters = np.max( - cut_tree(clusterization, height=height)) + 1 - nClusters = min(nClusters + self.kClusters, nClustersMax) - # Cut the tree to find the borehole groups - clusters = cut_tree( - clusterization, n_clusters=nClusters) - self.clusters = self.clusters + \ - [label + self.nEqBoreholes for label in clusters] - else: - nClusters = 1 - self.clusters.append(self.nEqBoreholes) - self.nEqBoreholes += nClusters - else: - self.nEqBoreholes = self.nBoreholes - self.clusters = range(self.nBoreholes) - # Overwrite boreholes with equivalent boreholes - self.boreholes = [_EquivalentBorehole( - [borehole - for borehole, cluster in zip(self.boreholes, self.clusters) - if cluster==i]) - for i in range(self.nEqBoreholes)] - self.wBoreholes = np.array([b.nBoreholes for b in self.boreholes]) - # Find similar pairs of boreholes - self.borehole_to_self, self.borehole_to_borehole = \ - self._find_axial_borehole_pairs(self.boreholes) - # Store unique distances in the bore field - self.dis = eqField.unique_distance(eqField, self.disTol)[0][1:] - - if self.boundary_condition == 'MIFT': - pipes = [self.network.p[self.clusters.index(i)] - for i in range(self.nEqBoreholes)] - self.network = _EquivalentNetwork( - self.boreholes, - pipes, - nSegments=self.nBoreSegments[0], - segment_ratios=self.segment_ratios[0]) - - # Stop chrono - toc = perf_counter() - if self.disp: - print(f' {toc - tic:.3f} sec') - print(f'Calculations will be done using {self.nEqBoreholes} ' - f'equivalent boreholes') - - return self.nBoreSegments[0]*self.nEqBoreholes - - def segment_lengths(self): - """ - Return the length of all segments in the bore field. - - The segments lengths are used for the energy balance in the calculation - of the g-function. For equivalent boreholes, the length of segments - is multiplied by the number of boreholes in the group. - - Returns - ------- - H : array - Array of segment lengths (in m). - - """ - # Borehole lengths - H = np.array([seg.H*seg.nBoreholes - for (borehole, nSegments, ratios) in zip( - self.boreholes, - self.nBoreSegments, - self.segment_ratios) - for seg in borehole.segments( - nSegments, segment_ratios=ratios)], - dtype=self.dtype) - return H - - def _compare_boreholes(self, borehole1, borehole2): - """ - Compare two boreholes and checks if they have the same dimensions : - H, D, and r_b. - - Parameters - ---------- - borehole1 : Borehole object - First borehole. - borehole2 : Borehole object - Second borehole. - - Returns - ------- - similarity : bool - True if the two boreholes have the same dimensions. - - """ - # Compare lengths (H), buried depth (D) and radius (r_b) - if (abs((borehole1.H - borehole2.H)/borehole1.H) < self.tol and - abs((borehole1.r_b - borehole2.r_b)/borehole1.r_b) < self.tol and - abs((borehole1.D - borehole2.D)/(borehole1.D + 1e-30)) < self.tol): - similarity = True - else: - similarity = False - return similarity - - def _compare_real_pairs(self, pair1, pair2): - """ - Compare two pairs of boreholes or segments and return True if the two - pairs have the same FLS solution for real sources. - - Parameters - ---------- - pair1 : Tuple of Borehole objects - First pair of boreholes or segments. - pair2 : Tuple of Borehole objects - Second pair of boreholes or segments. - - Returns - ------- - similarity : bool - True if the two pairs have the same FLS solution. - - """ - deltaD1 = pair1[1].D - pair1[0].D - deltaD2 = pair2[1].D - pair2[0].D - - # Equality of lengths between pairs - cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol - and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol) - # Equality of lengths in each pair - equal_H = abs((pair1[0].H - pair1[1].H)/pair1[0].H) < self.tol - # Equality of buried depths differences - cond_deltaD = abs(deltaD1 - deltaD2)/abs(deltaD1 + 1e-30) < self.tol - # Equality of buried depths differences if all boreholes have the same - # length - cond_deltaD_equal_H = abs((abs(deltaD1) - abs(deltaD2))/(abs(deltaD1) + 1e-30)) < self.tol - if cond_H and (cond_deltaD or (equal_H and cond_deltaD_equal_H)): - similarity = True - else: - similarity = False - return similarity - - def _compare_image_pairs(self, pair1, pair2): - """ - Compare two pairs of boreholes or segments and return True if the two - pairs have the same FLS solution for mirror sources. - - Parameters - ---------- - pair1 : Tuple of Borehole objects - First pair of boreholes or segments. - pair2 : Tuple of Borehole objects - Second pair of boreholes or segments. - - Returns - ------- - similarity : bool - True if the two pairs have the same FLS solution. - - """ - sumD1 = pair1[1].D + pair1[0].D - sumD2 = pair2[1].D + pair2[0].D - - # Equality of lengths between pairs - cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol - and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol) - # Equality of buried depths sums - cond_sumD = abs((sumD1 - sumD2)/(sumD1 + 1e-30)) < self.tol - if cond_H and cond_sumD: - similarity = True - else: - similarity = False - return similarity - - def _compare_realandimage_pairs(self, pair1, pair2): - """ - Compare two pairs of boreholes or segments and return True if the two - pairs have the same FLS solution for both real and mirror sources. - - Parameters - ---------- - pair1 : Tuple of Borehole objects - First pair of boreholes or segments. - pair2 : Tuple of Borehole objects - Second pair of boreholes or segments. - - Returns - ------- - similarity : bool - True if the two pairs have the same FLS solution. - - """ - if (self._compare_real_pairs(pair1, pair2) - and self._compare_image_pairs(pair1, pair2)): - similarity = True - else: - similarity = False - return similarity - - def _find_axial_borehole_pairs(self, boreholes): - """ - Find axial (i.e. disregarding the radial distance) similarities between - borehole pairs to simplify the evaluation of the FLS solution. - - Parameters - ---------- - boreholes : list of Borehole objects - Boreholes in the bore field. - - Returns - ------- - borehole_to_self : list - Lists of borehole indexes for each unique set of borehole - dimensions (H, D, r_b) in the bore field. - borehole_to_borehole : list - Lists of tuples of borehole indexes for each unique pair of - boreholes that share the same (pairwise) dimensions (H, D). - - """ - # Compare for the full (real + image) FLS solution - compare_pairs = self._compare_realandimage_pairs - - nBoreholes = len(boreholes) - borehole_to_self = [] - # Only check for similarities if there is more than one borehole - if nBoreholes > 1: - borehole_to_borehole = [] - for i, borehole_i in enumerate(boreholes): - # Compare the borehole to all known unique sets of dimensions - for k, borehole_set in enumerate(borehole_to_self): - m = borehole_set[0] - # Add the borehole to the group if a similar borehole is - # found - if self._compare_boreholes(borehole_i, boreholes[m]): - borehole_set.append(i) - break - else: - # If no similar boreholes are known, append the groups - borehole_to_self.append([i]) - # Note : The range is different from similarities since - # an equivalent borehole to itself includes borehole-to- - # borehole thermal interactions - for j, borehole_j in enumerate(boreholes[i:], start=i): - pair0 = (borehole_i, borehole_j) # pair - pair1 = (borehole_j, borehole_i) # reciprocal pair - # Compare pairs of boreholes to known unique pairs - for pairs in borehole_to_borehole: - m, n = pairs[0] - pair_ref = (boreholes[m], boreholes[n]) - # Add the pair (or the reciprocal pair) to a group - # if a similar one is found - if compare_pairs(pair0, pair_ref): - pairs.append((i, j)) - break - elif compare_pairs(pair1, pair_ref): - pairs.append((j, i)) - break - # If no similar pairs are known, append the groups - else: - borehole_to_borehole.append([(i, j)]) - else: - # Outputs for a single borehole - borehole_to_self = [[0]] - borehole_to_borehole = [[(0, 0)]] - return borehole_to_self, borehole_to_borehole - - def _find_unique_boreholes(self, boreholes): - """ - Find unique sets of dimensions (h, D, r_b) in the bore field. - - Parameters - ---------- - boreholes : list of Borehole objects - Boreholes in the bore field. - - Returns - ------- - unique_boreholes : list - List of list of borehole indices that correspond to unique - borehole dimensions (H, D, r_b). - - """ - unique_boreholes = [] - for i, borehole_1 in enumerate(boreholes): - for group in unique_boreholes: - borehole_2 = boreholes[group[0]] - # Add the borehole to a group if similar dimensions are found - if self._compare_boreholes(borehole_1, borehole_2): - group.append(i) - break - else: - # If no similar boreholes are known, append the groups - unique_boreholes.append([i]) - - return unique_boreholes - - def _find_unique_distances(self, dis, indices): - """ - Find the number of occurrences of each unique distances between pairs - of boreholes. - - Parameters - ---------- - dis : array - Array of unique distances (in meters) in the bore field. - indices : list - List of tuples of borehole indices. - - Returns - ------- - dis : array - Array of unique distances (in meters) in the bore field. - wDis : array - Array of number of occurences of each unique distance for each - pair of equivalent boreholes in indices. - - """ - wDis = np.zeros((len(dis), len(indices)), dtype=np.uint) - for k, pair in enumerate(indices): - i, j = pair - b1, b2 = self.boreholes[i], self.boreholes[j] - # Generate a flattened array of distances between boreholes i and j - if not i == j: - dis_ij = b1.distance(b2).flatten() - else: - # Remove the borehole radius from the distances - dis_ij = b1.distance(b2)[ - ~np.eye(b1.nBoreholes, dtype=bool)].flatten() - wDis_ij = np.zeros(len(dis), dtype=np.uint) - # Get insert positions for the distances - iDis = np.searchsorted(dis, dis_ij, side='left') - # Find indexes where previous index is closer - prev_iDis_is_less = ((iDis == len(dis))|(np.fabs(dis_ij - dis[np.maximum(iDis-1, 0)]) < np.fabs(dis_ij - dis[np.minimum(iDis, len(dis)-1)]))) - iDis[prev_iDis_is_less] -= 1 - np.add.at(wDis_ij, iDis, 1) - wDis[:,k] = wDis_ij - - return dis.reshape((1, -1)), wDis - - def _map_axial_segment_pairs(self, iBor, jBor, - reaSource=True, imgSource=True): - """ - Find axial (i.e. disregarding the radial distance) similarities between - segment pairs along two boreholes to simplify the evaluation of the - FLS solution. - - The returned H1, D1, H2, and D2 can be used to evaluate the segment-to- - segment response factors using scipy.integrate.quad_vec. - - Parameters - ---------- - iBor : int - Index of the first borehole. - jBor : int - Index of the second borehole. - - Returns - ------- - H1 : float - Length of the emitting segments. - D1 : array - Array of buried depths of the emitting segments. - H2 : float - Length of the receiving segments. - D2 : array - Array of buried depths of the receiving segments. - i_pair : list - Indices of the emitting segments along a borehole. - j_pair : list - Indices of the receiving segments along a borehole. - k_pair : list - Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions - corresponding to all pairs in (i_pair, j_pair). - - """ - # Initialize local variables - borehole1 = self.boreholes[iBor] - borehole2 = self.boreholes[jBor] - assert reaSource or imgSource, \ - "At least one of reaSource and imgSource must be True." - if reaSource and imgSource: - # Find segment pairs for the full (real + image) FLS solution - compare_pairs = self._compare_realandimage_pairs - elif reaSource: - # Find segment pairs for the real FLS solution - compare_pairs = self._compare_real_pairs - elif imgSource: - # Find segment pairs for the image FLS solution - compare_pairs = self._compare_image_pairs - # Dive both boreholes into segments - segments1 = borehole1.segments( - self.nBoreSegments[iBor], segment_ratios=self.segment_ratios[iBor]) - segments2 = borehole2.segments( - self.nBoreSegments[jBor], segment_ratios=self.segment_ratios[jBor]) - # Prepare lists of segment lengths - H1 = [] - H2 = [] - # Prepare lists of segment buried depths - D1 = [] - D2 = [] - # All possible pairs (i, j) of indices between segments - i_pair = np.repeat(np.arange(self.nBoreSegments[iBor], dtype=np.uint), - self.nBoreSegments[jBor]) - j_pair = np.tile(np.arange(self.nBoreSegments[jBor], dtype=np.uint), - self.nBoreSegments[iBor]) - # Empty list of indices for unique pairs - k_pair = np.empty(self.nBoreSegments[iBor] * self.nBoreSegments[jBor], - dtype=np.uint) - unique_pairs = [] - nPairs = 0 - - p = 0 - for i, segment_i in enumerate(segments1): - for j, segment_j in enumerate(segments2): - pair = (segment_i, segment_j) - # Compare the segment pairs to all known unique pairs - for k, pair_k in enumerate(unique_pairs): - m, n = pair_k[0], pair_k[1] - pair_ref = (segments1[m], segments2[n]) - # Stop if a similar pair is found and assign the index - if compare_pairs(pair, pair_ref): - k_pair[p] = k - break - # If no similar pair is found : add a new pair, increment the - # number of unique pairs, and extract the associated buried - # depths - else: - k_pair[p] = nPairs - H1.append(segment_i.H) - H2.append(segment_j.H) - D1.append(segment_i.D) - D2.append(segment_j.D) - unique_pairs.append((i, j)) - nPairs += 1 - p += 1 - return np.array(H1), np.array(D1), np.array(H2), np.array(D2), i_pair, j_pair, k_pair - - def _check_solver_specific_inputs(self): - """ - This method ensures that solver specific inputs to the Solver object - are what is expected. - - """ - assert type(self.disTol) is float and self.disTol > 0., \ - "The distance tolerance 'disTol' should be a positive float." - assert type(self.tol) is float and self.tol > 0., \ - "The relative tolerance 'tol' should be a positive float." - assert type(self.kClusters) is int and self.kClusters >= 0, \ - "The precision increment 'kClusters' should be a positive int." - assert np.all(np.array(self.nBoreSegments, dtype=np.uint) == self.nBoreSegments[0]), \ - "Solver 'equivalent' can only handle equal numbers of segments." - assert np.all([np.allclose(segment_ratios, self.segment_ratios[0]) for segment_ratios in self.segment_ratios]), \ - "Solver 'equivalent' can only handle identical segment_ratios for all boreholes." - assert not np.any([b.is_tilted() for b in self.boreholes]), \ - "Solver 'equivalent' can only handle vertical boreholes." - if self.boundary_condition == 'MIFT': - assert np.all(np.array(self.network.c, dtype=int) == -1), \ - "Solver 'equivalent' is only valid for parallel-connected " \ - "boreholes." - assert (self.m_flow_borehole is None - or (self.m_flow_borehole.ndim==1 and np.allclose(self.m_flow_borehole, self.m_flow_borehole[0])) - or (self.m_flow_borehole.ndim==2 and np.all([np.allclose(self.m_flow_borehole[:, i], self.m_flow_borehole[0, i]) for i in range(self.nBoreholes)]))), \ - "Mass flow rates into the network must be equal for all " \ - "boreholes." - # Use the total network mass flow rate. - if (type(self.network.m_flow_network) is np.ndarray and \ - len(self.network.m_flow_network)==len(self.network.b)): - self.network.m_flow_network = \ - self.network.m_flow_network[0]*len(self.network.b) - # Verify that all boreholes have the same piping configuration - # This is best done by comparing the matrix of thermal resistances. - assert np.all( - [np.allclose(self.network.p[0]._Rd, pipe._Rd) - for pipe in self.network.p]), \ - "All boreholes must have the same piping configuration." - return diff --git a/pygfunction/networks.py b/pygfunction/networks.py index 9fa5884c..8b5eb0bf 100644 --- a/pygfunction/networks.py +++ b/pygfunction/networks.py @@ -1,7 +1,17 @@ # -*- coding: utf-8 -*- +from typing import Union, List + import numpy as np +import numpy.typing as npt from scipy.linalg import block_diag +from .borefield import Borefield +from .boreholes import Borehole +from .enums import PipeType +from .media import Fluid +from .pipes import SingleUTube, MultipleUTube, Coaxial +from .pipes import fluid_to_pipe_thermal_resistance, fluid_to_fluid_thermal_resistance + class Network(object): """ @@ -9,7 +19,7 @@ class Network(object): connections between the boreholes. Contains information regarding the physical dimensions and thermal - characteristics of the pipes and the grout material in each boreholes, the + characteristics of the pipes and the grout material in each borehole, the topology of the connections between boreholes, as well as methods to evaluate fluid temperatures and heat extraction rates based on the work of Cimmino (2018, 2019, 2024) [#Network-Cimmin2018]_, [#Network-Cimmin2019]_, @@ -93,6 +103,150 @@ def __init__(self, boreholes, pipes, bore_connectivity=None, self._initialize_stored_coefficients( m_flow_network, cp_f, nSegments, segment_ratios) + @classmethod + def from_static_params(cls, + boreholes: Union[List[Borehole], Borefield], + pipe_type_str: str, + pos: List[tuple], + r_in: Union[float, tuple, npt.ArrayLike], + r_out: Union[float, tuple, npt.ArrayLike], + k_s: float, + k_g: float, + k_p: Union[float, tuple, npt.ArrayLike], + m_flow_network: float, + epsilon: float, + fluid_str: str, + fluid_concentraton_percent: float, + fluid_temperature: float, + reversible_flow: bool = True, + bore_connectivity: list = None, + J: int = 2): + """ + Constructs the 'Network' class from static parameters. + + Parameters + ---------- + boreholes : list of Borehole objects + List of boreholes included in the bore field. + pipe_type_str : str + Should be one of 'COAXIAL_ANNULAR_IN', 'COAXIAL_ANNULAR_OUT', + 'DOUBLE_UTUBE_PARALLEL', 'DOUBLE_UTUBE_SERIES', or 'SINGLE_UTUBE'. + pos : list of tuples + Position (x, y) (in meters) of the pipes inside the borehole. + r_in : float + Inner radius (in meters) of the U-Tube pipes. + r_out : float + Outer radius (in meters) of the U-Tube pipes. + k_s : float + Soil thermal conductivity (in W/m-K). + k_g : float + Grout thermal conductivity (in W/m-K). + k_p : float, tuple, or (2,) array + Pipe thermal conductivity (in W/m-K). + m_flow_network : float + Fluid mass flow rate into the network of boreholes (in kg/s). + epsilon : float + Pipe roughness (in meters). + fluid_str: str + The mixer for this application should be one of: + + - 'Water' - Complete water solution + - 'MEG' - Ethylene glycol mixed with water + - 'MPG' - Propylene glycol mixed with water + - 'MEA' - Ethanol mixed with water + - 'MMA' - Methanol mixed with water + + fluid_concentration_pct: float + Mass fraction of the mixing fluid added to water (in %). + Lower bound = 0. Upper bound is dependent on the mixture. + fluid_temperature: float, optional + Temperature used for evaluating fluid properties (in degC). + Default is 20. + reversible_flow : bool, optional + True to treat a negative mass flow rate as the reversal of flow + direction within the borehole. If False, the direction of flow is not + reversed when the mass flow rate is negative, and the absolute value is + used for calculations. + Default is True. + bore_connectivity : list, optional + Index of fluid inlet into each borehole. -1 corresponds to a borehole + connected to the bore field inlet. If this parameter is not provided, + parallel connections between boreholes is used. + Default is None. + J : int, optional + Number of multipoles per pipe to evaluate the thermal resistances. + J=1 or J=2 usually gives sufficient accuracy. J=0 corresponds to the + line source approximation. + Default is 2. + + Returns + ------- + Network : 'Network' object. + The network. + + """ + # Convert borefield to list + if isinstance(boreholes, Borefield): + boreholes = boreholes.to_boreholes() + + # The total fluid mass flow rate is divided equally amongst inlets + if bore_connectivity is None: + m_flow_borehole = abs(m_flow_network / len(boreholes)) + else: + m_flow_borehole = abs(m_flow_network / bore_connectivity.count(-1)) + + # Pipe and fluid types + pipe_type = PipeType[pipe_type_str.upper()] + fluid = Fluid(fluid_str, fluid_concentraton_percent, fluid_temperature) + + if pipe_type == PipeType.SINGLE_UTUBE: + # Single U-tube borehole + R_fp = fluid_to_pipe_thermal_resistance( + pipe_type, m_flow_borehole, r_in, r_out, k_p, epsilon, fluid) + pipes = [ + SingleUTube( + pos, r_in, r_out, borehole, k_s, k_g, R_fp, J, reversible_flow) + for borehole in boreholes + ] + + elif pipe_type == PipeType.DOUBLE_UTUBE_PARALLEL: + # Double U-tube borehole (parallel) + R_fp = fluid_to_pipe_thermal_resistance( + pipe_type, m_flow_borehole, r_in, r_out, k_p, epsilon, fluid) + pipes = [ + MultipleUTube( + pos, r_in, r_out, borehole, k_s, k_g, R_fp, 2, 'parallel', J, reversible_flow) + for borehole in boreholes + ] + + elif pipe_type == PipeType.DOUBLE_UTUBE_SERIES: + # Double U-tube borehole (series) + R_fp = fluid_to_pipe_thermal_resistance( + pipe_type, m_flow_borehole, r_in, r_out, k_p, epsilon, fluid) + pipes = [ + MultipleUTube( + pos, r_in, r_out, borehole, k_s, k_g, R_fp, 2, 'series', J, reversible_flow) + for borehole in boreholes + ] + + elif pipe_type in [PipeType.COAXIAL_ANNULAR_IN, PipeType.COAXIAL_ANNULAR_OUT]: + # Coaxial borehole + R_fp = fluid_to_pipe_thermal_resistance( + pipe_type, m_flow_borehole, r_in, r_out, k_p, epsilon, fluid) + R_ff = fluid_to_fluid_thermal_resistance( + pipe_type, m_flow_borehole, r_in, r_out, k_p, epsilon, fluid) + pipes = [ + Coaxial( + pos, np.array(r_in), np.array(r_out), borehole, k_s, k_g, R_ff, R_fp, J, reversible_flow) + for borehole in boreholes + ] + + else: + raise ValueError(f"Unsupported pipe_type: '{pipe_type_str}'") + + return cls(boreholes=boreholes, pipes=pipes, m_flow_network=m_flow_network, bore_connectivity=bore_connectivity, + cp_f=fluid.cp) + def get_inlet_temperature( self, T_f_in, T_b, m_flow_network, cp_f, nSegments, segment_ratios=None): @@ -577,7 +731,7 @@ def coefficients_outlet_temperature( def coefficients_network_inlet_temperature( self, m_flow_network, cp_f, nSegments, segment_ratios=None): """ - Build coefficient matrices to evaluate intlet fluid temperature of the + Build coefficient matrices to evaluate inlet fluid temperature of the network. Returns coefficients for the relation: @@ -992,11 +1146,11 @@ def _initialize_coefficients_connectivity(self): def _initialize_stored_coefficients( self, m_flow_network, cp_f, nSegments, segment_ratios): nMethods = 7 # Number of class methods - self._stored_coefficients = [() for i in range(nMethods)] + self._stored_coefficients = [() for _ in range(nMethods)] self._stored_m_flow_cp = [np.empty(self.nInlets)*np.nan - for i in range(nMethods)] - self._stored_nSegments = [np.nan for i in range(nMethods)] - self._stored_segment_ratios = [np.nan for i in range(nMethods)] + for _ in range(nMethods)] + self._stored_nSegments = [np.nan for _ in range(nMethods)] + self._stored_segment_ratios = [np.nan for _ in range(nMethods)] self._m_flow_cp_model_variables = np.empty(self.nInlets)*np.nan self._nSegments_model_variables = np.nan self._segment_ratios_model_variables = np.nan diff --git a/pygfunction/pipes.py b/pygfunction/pipes.py index 0d86a526..c5ceb5ef 100644 --- a/pygfunction/pipes.py +++ b/pygfunction/pipes.py @@ -1,13 +1,16 @@ # -*- coding: utf-8 -*- import warnings +from typing import Union import numpy as np +import numpy.typing as npt from scipy.constants import pi from scipy.special import binom +from .enums import PipeType +from .media import Fluid from .utilities import _initialize_figure, _format_axes - class _BasePipe(object): """ Template for pipe classes. @@ -3506,7 +3509,7 @@ def _F_mk(q_p, P, n_p, J, r_b, r_out, z, pikg, sigma): P : array Multipoles. n_p : int - Total numper of pipes. + Total number of pipes. J : int Number of multipoles per pipe to evaluate the thermal resistances. J=1 or J=2 usually gives sufficient accuracy. J=0 corresponds to the @@ -3613,3 +3616,224 @@ def _Nusselt_number_turbulent_flow(Re, Pr, fDarcy): Nu = 0.125 * fDarcy * (Re - 1.0e3) * Pr / \ (1.0 + 12.7 * np.sqrt(0.125*fDarcy) * (Pr**(2.0/3.0) - 1.0)) return Nu + + +def fluid_to_pipe_thermal_resistance( + pipe_type: PipeType, m_flow_borehole: float, + r_in: Union[float, tuple, npt.ArrayLike], r_out: Union[float, tuple, npt.ArrayLike], + k_p: Union[float, tuple, npt.ArrayLike], epsilon: float, + fluid: Fluid) -> float: + """ + Computes the fluid to pipe thermal resistance. + + Parameters + ---------- + pipe_type : PipeType + Should be one of 'PipeType.COAXIAL_ANNULAR_IN', 'PipeType.COAXIAL_ANNULAR_OUT', + 'PipeType.DOUBLE_UTUBE_PARALLEL', 'PipeType.DOUBLE_UTUBE_SERIES', or 'PipeType.SINGLE_UTUBE'. + m_flow_borehole : float + Fluid mass flow rate the borehole (in kg/s). + r_in : float + Inner radius (in meters) of the U-Tube pipes. + r_out : float + Outer radius (in meters) of the U-Tube pipes. + k_p : float + Pipe thermal conductivity (in W/m-K). + epsilon : float + Pipe roughness (in meters). + fluid : Fluid + 'Fluid' class object. Used for evaluating fluid properties + + Returns + ------- + float + fluid to pipe thermal resistance (in m-K/W) + + """ + + if pipe_type in [PipeType.SINGLE_UTUBE, PipeType.DOUBLE_UTUBE_SERIES]: + + # The fluid mass flow rate corresponds to the total flow + m_flow_pipe = m_flow_borehole + + # Pipe thermal resistance + R_p = conduction_thermal_resistance_circular_pipe( + r_in, r_out, k_p) + # Convection heat transfer coefficient [W/m2.K] + h_f = convective_heat_transfer_coefficient_circular_pipe( + m_flow_pipe, r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp, + epsilon) + # Film thermal resistance [m.K/W] + R_f = 1.0 / (h_f * 2 * np.pi * r_in) + + return R_p + R_f + + elif pipe_type == PipeType.DOUBLE_UTUBE_PARALLEL: + + # The fluid mass flow rate is divided into the two parallel pipes + m_flow_pipe = m_flow_borehole / 2 + + # Pipe thermal resistance + R_p = conduction_thermal_resistance_circular_pipe( + r_in, r_out, k_p) + # Convection heat transfer coefficient [W/m2.K] + h_f = convective_heat_transfer_coefficient_circular_pipe( + m_flow_pipe, r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp, + epsilon) + # Film thermal resistance [m.K/W] + R_f = 1.0 / (h_f * 2 * np.pi * r_in) + + return R_p + R_f + + elif pipe_type == PipeType.COAXIAL_ANNULAR_IN: + + # The fluid mass flow rate corresponds to the total flow + m_flow_pipe = m_flow_borehole + + # The annular channel is at index 0 + r_in_out = r_out[1] + r_out_in = r_in[0] + r_out_out = r_out[0] + k_p_out = k_p[0] + + # Outer pipe + R_p_out = conduction_thermal_resistance_circular_pipe( + r_out_in, r_out_out, k_p_out) + + # Outer pipe + h_f_a_in, h_f_a_out = \ + convective_heat_transfer_coefficient_concentric_annulus( + m_flow_pipe, r_in_out, r_out_in, fluid.mu, fluid.rho, fluid.k, + fluid.cp, epsilon) + + # Coaxial GHE in borehole + R_f_out_out = 1.0 / (h_f_a_out * 2 * np.pi * r_out_in) + return R_p_out + R_f_out_out + + elif pipe_type == PipeType.COAXIAL_ANNULAR_OUT: + + # The fluid mass flow rate corresponds to the total flow + m_flow_pipe = m_flow_borehole + + # The annular channel is at index 1 + r_in_out = r_out[0] + r_out_in = r_in[1] + r_out_out = r_out[1] + k_p_out = k_p[1] + + # Outer pipe + R_p_out = conduction_thermal_resistance_circular_pipe( + r_out_in, r_out_out, k_p_out) + # Fluid-to-fluid thermal resistance [m.K/W] + + # Outer pipe + h_f_a_in, h_f_a_out = \ + convective_heat_transfer_coefficient_concentric_annulus( + m_flow_pipe, r_in_out, r_out_in, fluid.mu, fluid.rho, fluid.k, + fluid.cp, epsilon) + + # Coaxial GHE in borehole + R_f_out_out = 1.0 / (h_f_a_out * 2 * np.pi * r_out_in) + + return R_p_out + R_f_out_out + + else: + raise ValueError(f"Unsupported pipe_type: '{pipe_type.name}'") + + +def fluid_to_fluid_thermal_resistance(pipe_type: PipeType, m_flow_borehole: float, + r_in: Union[float, tuple, npt.ArrayLike], + r_out: Union[float, tuple, npt.ArrayLike], + k_p: Union[float, tuple, npt.ArrayLike], epsilon: float, + fluid: Fluid) -> float: + """ + Computes the fluid to fluid thermal resistance. + + Parameters + ---------- + pipe_type : PipeType + Should be one of 'PipeType.COAXIAL_ANNULAR_IN', 'PipeType.COAXIAL_ANNULAR_OUT', + 'PipeType.DOUBLE_UTUBE_PARALLEL', 'PipeType.DOUBLE_UTUBE_SERIES', or 'PipeType.SINGLE_UTUBE'. + m_flow_borehole : float + Fluid mass flow rate the borehole (in kg/s). + r_in : float + Inner radius (in meters) of the U-Tube pipes. + r_out : float + Outer radius (in meters) of the U-Tube pipes. + k_p : float + Pipe thermal conductivity (in W/m-K). + epsilon : float + Pipe roughness (in meters). + fluid : Fluid + 'Fluid' class object. Used for evaluating fluid properties + + Returns + ------- + float + fluid to fluid thermal resistance (in m-K/W) + + """ + + if pipe_type == PipeType.COAXIAL_ANNULAR_IN: + + # The fluid mass flow rate corresponds to the total flow + m_flow_pipe = m_flow_borehole + + # The annular channel is at index 0 + r_in_in = r_in[1] + r_in_out = r_out[1] + r_out_in = r_in[0] + k_p_in = k_p[1] + + # Inner pipe + R_p_in = conduction_thermal_resistance_circular_pipe( + r_in_in, r_in_out, k_p_in) + + # Fluid-to-fluid thermal resistance [m.K/W] + # Inner pipe + h_f_in = convective_heat_transfer_coefficient_circular_pipe( + m_flow_pipe, r_in_in, fluid.mu, fluid.rho, fluid.k, fluid.cp, epsilon) + R_f_in = 1.0 / (h_f_in * 2 * np.pi * r_in_in) + + # Outer pipe + h_f_a_in, h_f_a_out = \ + convective_heat_transfer_coefficient_concentric_annulus( + m_flow_borehole, r_in_out, r_out_in, fluid.mu, fluid.rho, fluid.k, + fluid.cp, epsilon) + R_f_out_in = 1.0 / (h_f_a_in * 2 * np.pi * r_in_out) + + return R_f_in + R_p_in + R_f_out_in + + elif pipe_type == PipeType.COAXIAL_ANNULAR_OUT: + + # The fluid mass flow rate corresponds to the total flow + m_flow_pipe = m_flow_borehole + + # The annular channel is at index 1 + r_in_in = r_in[0] + r_in_out = r_out[0] + r_out_in = r_in[1] + k_p_in = k_p[0] + + # Pipe thermal resistances [m.K/W] + # Inner pipe + R_p_in = conduction_thermal_resistance_circular_pipe( + r_in_in, r_in_out, k_p_in) + + # Fluid-to-fluid thermal resistance [m.K/W] + # Inner pipe + h_f_in = convective_heat_transfer_coefficient_circular_pipe( + m_flow_pipe, r_in_in, fluid.mu, fluid.rho, fluid.k, fluid.cp, epsilon) + R_f_in = 1.0 / (h_f_in * 2 * np.pi * r_in_in) + + # Outer pipe + h_f_a_in, h_f_a_out = \ + convective_heat_transfer_coefficient_concentric_annulus( + m_flow_pipe, r_in_out, r_out_in, fluid.mu, fluid.rho, fluid.k, + fluid.cp, epsilon) + R_f_out_in = 1.0 / (h_f_a_in * 2 * np.pi * r_in_out) + + return R_f_in + R_p_in + R_f_out_in + + else: + raise ValueError(f"Unsupported pipe_type: '{pipe_type.name}'") diff --git a/pygfunction/solvers/__init__.py b/pygfunction/solvers/__init__.py new file mode 100644 index 00000000..4ea6ece6 --- /dev/null +++ b/pygfunction/solvers/__init__.py @@ -0,0 +1,11 @@ +from ._base_solver import _BaseSolver +from .detailed import Detailed +from .equivalent import Equivalent +from .similarities import Similarities + +__all__ = [ +'_BaseSolver', +'Detailed', +'Equivalent', +'Similarities' +] diff --git a/pygfunction/solvers/_base_solver.py b/pygfunction/solvers/_base_solver.py new file mode 100644 index 00000000..e73eb2f9 --- /dev/null +++ b/pygfunction/solvers/_base_solver.py @@ -0,0 +1,645 @@ +# -*- coding: utf-8 -*- +from time import perf_counter + +import numpy as np +from scipy.interpolate import interp1d as interp1d + +from ..borefield import Borefield +from ..boreholes import Borehole +from ..networks import ( + Network, + network_thermal_resistance + ) +from ..utilities import segment_ratios + + +class _BaseSolver: + """ + Template for solver classes. + + Solver classes inherit from this class. + + Attributes + ---------- + boreholes : list of Borehole objects + List of boreholes included in the bore field. + network : network object + Model of the network. + time : float or array + Values of time (in seconds) for which the g-function is evaluated. + boundary_condition : str + Boundary condition for the evaluation of the g-function. Should be one + of + + - 'UHTR' : + Uniform heat transfer rate. + - 'UBWT' : + Uniform borehole wall temperature. + - 'MIFT' : + Mixed inlet fluid temperatures. + + nSegments : int or list, optional + Number of line segments used per borehole, or list of number of + line segments used for each borehole. + Default is 8. + segment_ratios : array, list of arrays, or callable, optional + Ratio of the borehole length represented by each segment. The + sum of ratios must be equal to 1. The shape of the array is of + (nSegments,) or list of (nSegments[i],). If segment_ratios==None, + segments of equal lengths are considered. If a callable is provided, it + must return an array of size (nSegments,) when provided with nSegments + (of type int) as an argument, or an array of size (nSegments[i],) when + provided with an element of nSegments (of type list). + Default is :func:`utilities.segment_ratios`. + m_flow_borehole : (nInlets,) array or (nMassFlow, nInlets,) array, optional + Fluid mass flow rate into each circuit of the network. If a + (nMassFlow, nInlets,) array is supplied, the + (nMassFlow, nMassFlow,) variable mass flow rate g-functions + will be evaluated using the method of Cimmino (2024) + [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary + condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be + provided. + Default is None. + m_flow_network : float or (nMassFlow,) array, optional + Fluid mass flow rate into the network of boreholes. If an array + is supplied, the (nMassFlow, nMassFlow,) variable mass flow + rate g-functions will be evaluated using the method of Cimmino + (2024) [#gFunction-CimBer2024]_. Only required for the 'MIFT' boundary + condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be + provided. + Default is None. + cp_f : float, optional + Fluid specific isobaric heat capacity (in J/kg.degC). Only required + for the 'MIFT' boundary condition. + Default is None. + approximate_FLS : bool, optional + Set to true to use the approximation of the FLS solution of Cimmino + (2021). This approximation does not require the numerical evaluation of + any integral. When using the 'equivalent' solver, the approximation is + only applied to the thermal response at the borehole radius. Thermal + interaction between boreholes is evaluated using the FLS solution. + Default is False. + nFLS : int, optional + Number of terms in the approximation of the FLS solution. This + parameter is unused if `approximate_FLS` is set to False. + Default is 10. Maximum is 25. + mQuad : int, optional + Number of Gauss-Legendre sample points for the integral over :math:`u` + in the inclined FLS solution. + Default is 11. + linear_threshold : float, optional + Threshold time (in seconds) under which the g-function is + linearized. The g-function value is then interpolated between 0 + and its value at the threshold. If linear_threshold==None, the + g-function is linearized for times + `t < r_b**2 / (25 * self.alpha)`. + Default is None. + disp : bool, optional + Set to true to print progression messages. + Default is False. + profiles : bool, optional + Set to true to keep in memory the temperatures and heat extraction + rates. + Default is False. + kind : string, optional + Interpolation method used for segment-to-segment thermal response + factors. See documentation for scipy.interpolate.interp1d. + Default is 'linear'. + dtype : numpy dtype, optional + numpy data type used for matrices and vectors. Should be one of + numpy.single or numpy.double. + Default is numpy.double. + + """ + def __init__(self, boreholes, network, time, boundary_condition, + m_flow_borehole=None, m_flow_network=None, cp_f=None, + nSegments=8, segment_ratios=segment_ratios, + approximate_FLS=False, mQuad=11, nFLS=10, + linear_threshold=None, disp=False, profiles=False, + kind='linear', dtype=np.double, **other_options): + self.boreholes = boreholes + self.network = network + # Convert time to a 1d array + self.time = np.atleast_1d(time).flatten() + self.linear_threshold = linear_threshold + self.r_b_max = np.max([b.r_b for b in self.boreholes]) + self.boundary_condition = boundary_condition + nBoreholes = len(self.boreholes) + # Format number of segments and segment ratios + if type(nSegments) is int: + self.nBoreSegments = [nSegments] * nBoreholes + else: + self.nBoreSegments = nSegments + if isinstance(segment_ratios, np.ndarray): + segment_ratios = [segment_ratios] * nBoreholes + elif segment_ratios is None: + segment_ratios = [np.full(n, 1./n) for n in self.nBoreSegments] + elif callable(segment_ratios): + segment_ratios = [segment_ratios(n) for n in self.nBoreSegments] + self.segment_ratios = segment_ratios + # Shortcut for segment_ratios comparisons + self._equal_segment_ratios = \ + (np.all(np.array(self.nBoreSegments, dtype=np.uint) == self.nBoreSegments[0]) + and np.all([np.allclose(segment_ratios, self.segment_ratios[0]) for segment_ratios in self.segment_ratios])) + # Boreholes with a uniform discretization + self._uniform_segment_ratios = [ + np.allclose(segment_ratios, + segment_ratios[0:1], + rtol=1e-6) + for segment_ratios in self.segment_ratios] + # Find indices of first and last segments along boreholes + self._i0Segments = [sum(self.nBoreSegments[0:i]) + for i in range(nBoreholes)] + self._i1Segments = [sum(self.nBoreSegments[0:(i + 1)]) + for i in range(nBoreholes)] + self.nMassFlow = 0 + self.m_flow_borehole = m_flow_borehole + if self.m_flow_borehole is not None: + if not self.m_flow_borehole.ndim == 1: + self.nMassFlow = np.size(self.m_flow_borehole, axis=0) + self.m_flow_borehole = np.atleast_2d(self.m_flow_borehole) + self.m_flow = self.m_flow_borehole + self.m_flow_network = m_flow_network + if self.m_flow_network is not None: + if not isinstance(self.m_flow_network, (np.floating, float)): + self.nMassFlow = len(self.m_flow_network) + self.m_flow_network = np.atleast_1d(self.m_flow_network) + self.m_flow = self.m_flow_network + self.cp_f = cp_f + self.approximate_FLS = approximate_FLS + self.mQuad = mQuad + self.nFLS = nFLS + self.disp = disp + self.profiles = profiles + self.kind = kind + self.dtype = dtype + # Check the validity of inputs + self._check_inputs() + # Initialize the solver with solver-specific options + self.nSources = self.initialize(**other_options) + + return + + def initialize(self, *kwargs): + """ + Perform any calculation required at the initialization of the solver + and returns the number of finite line heat sources in the borefield. + + Raises + ------ + NotImplementedError + + Returns + ------- + nSources : int + Number of finite line heat sources in the borefield used to + initialize the matrix of segment-to-segment thermal response + factors (of size: nSources x nSources). + + """ + raise NotImplementedError( + 'initialize class method not implemented, this method should ' + 'return the number of finite line heat sources in the borefield ' + 'used to initialize the matrix of segment-to-segment thermal ' + 'response factors (of size: nSources x nSources)') + return None + + def solve(self, time, alpha): + """ + Build and solve the system of equations. + + Parameters + ---------- + time : float or array + Values of time (in seconds) for which the g-function is evaluated. + alpha : float + Soil thermal diffusivity (in m2/s). + + Returns + ------- + gFunc : float or array + Values of the g-function + + """ + # Number of time values + self.time = time + nt = len(self.time) + # Evaluate threshold time for g-function linearization + if self.linear_threshold is None: + time_threshold = self.r_b_max**2 / (25 * alpha) + else: + time_threshold = self.linear_threshold + # Find the number of g-function values to be linearized + p_long = np.searchsorted(self.time, time_threshold, side='right') + if p_long > 0: + time_long = np.concatenate([[time_threshold], self.time[p_long:]]) + else: + time_long = self.time + nt_long = len(time_long) + # Calculate segment to segment thermal response factors + h_ij = self.thermal_response_factors(time_long, alpha, kind=self.kind) + # Segment lengths + H_b = self.segment_lengths() + if self.boundary_condition == 'MIFT': + Hb_individual = np.array([b.H for b in self.boreSegments], dtype=self.dtype) + H_tot = np.sum(H_b) + if self.disp: print('Building and solving the system of equations ...', + end='') + # Initialize chrono + tic = perf_counter() + + if self.boundary_condition == 'UHTR': + # Initialize g-function + gFunc = np.zeros(nt) + # Initialize segment heat extraction rates + Q_b = 1 + # Initialize borehole wall temperatures + T_b = np.zeros((self.nSources, nt), dtype=self.dtype) + + # Build and solve the system of equations at all times + p0 = max(0, p_long-1) + for p in range(nt_long): + # Evaluate the g-function with uniform heat extraction along + # boreholes + + # Thermal response factors evaluated at time t[p] + h_dt = h_ij.y[:,:,p+1] + # Borehole wall temperatures are calculated by the sum of + # contributions of all segments + T_b[:,p+p0] = np.sum(h_dt, axis=1) + # The g-function is the average of all borehole wall + # temperatures + gFunc[p+p0] = np.sum(T_b[:,p+p0]*H_b)/H_tot + + # Linearize g-function for times under threshold + if p_long > 0: + gFunc[:p_long] = gFunc[p_long-1] * self.time[:p_long] / time_threshold + T_b[:,:p_long] = T_b[:,p_long-1:p_long] * self.time[:p_long] / time_threshold + + elif self.boundary_condition == 'UBWT': + # Initialize g-function + gFunc = np.zeros(nt) + # Initialize segment heat extraction rates + Q_b = np.zeros((self.nSources, nt), dtype=self.dtype) + T_b = np.zeros(nt, dtype=self.dtype) + + # Build and solve the system of equations at all times + p0 = max(0, p_long-1) + for p in range(nt_long): + # Current thermal response factor matrix + if p > 0: + dt = time_long[p] - time_long[p-1] + else: + dt = time_long[p] + # Thermal response factors evaluated at t=dt + h_dt = h_ij(dt) + # Reconstructed load history + Q_reconstructed = self.load_history_reconstruction( + time_long[0:p+1], Q_b[:,p0:p+p0+1]) + # Borehole wall temperature for zero heat extraction at + # current step + T_b0 = self.temporal_superposition( + h_ij.y[:,:,1:], Q_reconstructed) + + # Evaluate the g-function with uniform borehole wall + # temperature + # --------------------------------------------------------- + # Build a system of equation [A]*[X] = [B] for the + # evaluation of the g-function. [A] is a coefficient + # matrix, [X] = [Q_b,T_b] is a state space vector of the + # borehole heat extraction rates and borehole wall + # temperature (equal for all segments), [B] is a + # coefficient vector. + # + # Spatial superposition: [T_b] = [T_b0] + [h_ij_dt]*[Q_b] + # Energy conservation: sum([Q_b*Hb]) = sum([Hb]) + # --------------------------------------------------------- + A = np.block([[h_dt, -np.ones((self.nSources, 1), + dtype=self.dtype)], + [H_b, 0.]]) + B = np.hstack((-T_b0, H_tot)) + # Solve the system of equations + X = np.linalg.solve(A, B) + # Store calculated heat extraction rates + Q_b[:,p+p0] = X[0:self.nSources] + # The borehole wall temperatures are equal for all segments + T_b[p+p0] = X[-1] + gFunc[p+p0] = T_b[p+p0] + + # Linearize g-function for times under threshold + if p_long > 0: + gFunc[:p_long] = gFunc[p_long-1] * self.time[:p_long] / time_threshold + Q_b[:,:p_long] = 1 + (Q_b[:,p_long-1:p_long] - 1) * self.time[:p_long] / time_threshold + T_b[:p_long] = T_b[p_long-1] * self.time[:p_long] / time_threshold + + elif self.boundary_condition == 'MIFT': + if self.nMassFlow == 0: + # Initialize g-function + gFunc = np.zeros((1, 1, nt)) + # Initialize segment heat extraction rates + Q_b = np.zeros((1, self.nSources, nt), dtype=self.dtype) + T_b = np.zeros((1, self.nSources, nt), dtype=self.dtype) + else: + # Initialize g-function + gFunc = np.zeros((self.nMassFlow, self.nMassFlow, nt)) + # Initialize segment heat extraction rates + Q_b = np.zeros( + (self.nMassFlow, self.nSources, nt), dtype=self.dtype) + T_b = np.zeros( + (self.nMassFlow, self.nSources, nt), dtype=self.dtype) + + for j in range(np.maximum(self.nMassFlow, 1)): + # Build and solve the system of equations at all times + p0 = max(0, p_long-1) + a_in_j, a_b_j = self.network.coefficients_borehole_heat_extraction_rate( + self.m_flow[j], + self.cp_f, + self.nBoreSegments, + segment_ratios=self.segment_ratios) + k_s = self.network.p[0].k_s + for p in range(nt_long): + # Current thermal response factor matrix + if p > 0: + dt = time_long[p] - time_long[p-1] + else: + dt = time_long[p] + # Thermal response factors evaluated at t=dt + h_dt = h_ij(dt) + # Reconstructed load history + Q_reconstructed = self.load_history_reconstruction( + time_long[0:p+1], Q_b[j,:,p0:p+p0+1]) + # Borehole wall temperature for zero heat extraction at + # current step + T_b0 = self.temporal_superposition( + h_ij.y[:,:,1:], Q_reconstructed) + + # Evaluate the g-function with mixed inlet fluid + # temperatures + # --------------------------------------------------------- + # Build a system of equation [A]*[X] = [B] for the + # evaluation of the g-function. [A] is a coefficient + # matrix, [X] = [Q_b,T_b,Tf_in] is a state space vector of + # the borehole heat extraction rates, borehole wall + # temperatures and inlet fluid temperature (into the bore + # field), [B] is a coefficient vector. + # + # Spatial superposition: [T_b] = [T_b0] + [h_ij_dt]*[Q_b] + # Heat transfer inside boreholes: + # [Q_{b,i}] = [a_in]*[T_{f,in}] + [a_{b,i}]*[T_{b,i}] + # Energy conservation: sum([Q_b*H_b]) = sum([H_b]) + # --------------------------------------------------------- + A = np.block( + [[h_dt, + -np.eye(self.nSources, dtype=self.dtype), + np.zeros((self.nSources, 1), dtype=self.dtype)], + [np.eye(self.nSources, dtype=self.dtype), + a_b_j/(2.0*np.pi*k_s*np.atleast_2d(Hb_individual).T), + a_in_j/(2.0*np.pi*k_s*np.atleast_2d(Hb_individual).T)], + [H_b, np.zeros(self.nSources + 1, dtype=self.dtype)]]) + B = np.hstack( + (-T_b0, + np.zeros(self.nSources, dtype=self.dtype), + H_tot)) + # Solve the system of equations + X = np.linalg.solve(A, B) + # Store calculated heat extraction rates + Q_b[j,:,p+p0] = X[0:self.nSources] + T_b[j,:,p+p0] = X[self.nSources:2*self.nSources] + # Inlet fluid temperature + T_f_in = X[-1] + # The gFunction is equal to the effective borehole wall + # temperature + # Outlet fluid temperature + T_f_out = T_f_in - 2 * np.pi * k_s * H_tot / ( + np.sum(np.abs(self.m_flow[j]) * self.cp_f)) + # Average fluid temperature + T_f = 0.5*(T_f_in + T_f_out) + # Borefield thermal resistance + R_field = network_thermal_resistance( + self.network, self.m_flow[j], self.cp_f) + # Effective borehole wall temperature + T_b_eff = T_f - 2 * np.pi * k_s * R_field + gFunc[j,j,p+p0] = T_b_eff + + for i in range(np.maximum(self.nMassFlow, 1)): + for j in range(np.maximum(self.nMassFlow, 1)): + if not i == j: + # Inlet fluid temperature + a_in, a_b = self.network.coefficients_network_heat_extraction_rate( + self.m_flow[i], + self.cp_f, + self.nBoreSegments, + segment_ratios=self.segment_ratios) + T_f_in = (-2 * np.pi * k_s * H_tot - a_b @ T_b[j,:,p0:]) / a_in + # The gFunction is equal to the effective borehole wall + # temperature + # Outlet fluid temperature + T_f_out = T_f_in - 2 * np.pi * k_s * H_tot / np.sum(np.abs(self.m_flow[i]) * self.cp_f) + # Borefield thermal resistance + R_field = network_thermal_resistance( + self.network, self.m_flow[i], self.cp_f) + # Effective borehole wall temperature + T_b_eff = 0.5 * (T_f_in + T_f_out) - 2 * np.pi * k_s * R_field + gFunc[i,j,p0:] = T_b_eff + + # Linearize g-function for times under threshold + if p_long > 0: + gFunc[:,:,:p_long] = gFunc[:,:,p_long-1] * self.time[:p_long] / time_threshold + Q_b[:,:,:p_long] = 1 + (Q_b[:,:,p_long-1:p_long] - 1) * self.time[:p_long] / time_threshold + T_b[:,:,:p_long] = T_b[:,:,p_long-1:p_long] * self.time[:p_long] / time_threshold + if self.nMassFlow == 0: + gFunc = gFunc[0,0,:] + Q_b = Q_b[0,:,:] + T_b = T_b[0,:,:] + + # Store temperature and heat extraction rate profiles + if self.profiles: + self.Q_b = Q_b + self.T_b = T_b + toc = perf_counter() + if self.disp: print(f' {toc - tic:.3f} sec') + return gFunc + + def segment_lengths(self): + """ + Return the length of all segments in the bore field. + + The segments lengths are used for the energy balance in the calculation + of the g-function. + + Returns + ------- + H : array + Array of segment lengths (in m). + + """ + # Borehole lengths + H_b = np.array([b.H for b in self.boreSegments], dtype=self.dtype) + return H_b + + def borehole_segments(self): + """ + Split boreholes into segments. + + This function goes through the list of boreholes and builds a new list, + with each borehole split into nSegments of equal lengths. + + Returns + ------- + boreSegments : list + List of borehole segments. + + """ + boreSegments = [] # list for storage of boreSegments + for b, nSegments, segment_ratios in zip(self.boreholes, self.nBoreSegments, self.segment_ratios): + segments = b.segments(nSegments, segment_ratios=segment_ratios) + boreSegments.extend(segments) + + return boreSegments + + def temporal_superposition(self, h_ij, Q_reconstructed): + """ + Temporal superposition for inequal time steps. + + Parameters + ---------- + h_ij : array + Values of the segment-to-segment thermal response factor increments + at the given time step. + Q_reconstructed : array + Reconstructed heat extraction rates of all segments at all times. + + Returns + ------- + T_b0 : array + Current values of borehole wall temperatures assuming no heat + extraction during current time step. + + """ + # Number of heat sources + nSources = Q_reconstructed.shape[0] + # Number of time steps + nt = Q_reconstructed.shape[1] + # Spatial and temporal superpositions + dQ = np.concatenate( + (Q_reconstructed[:,0:1], + Q_reconstructed[:,1:] - Q_reconstructed[:,0:-1]), axis=1)[:,::-1] + # Borehole wall temperature + T_b0 = np.einsum('ijk,jk', h_ij[:,:,:nt], dQ) + + return T_b0 + + def load_history_reconstruction(self, time, Q_b): + """ + Reconstructs the load history. + + This function calculates an equivalent load history for an inverted + order of time step sizes. + + Parameters + ---------- + time : array + Values of time (in seconds) in the load history. + Q_b : array + Heat extraction rates (in Watts) of all segments at all times. + + Returns + ------- + Q_reconstructed : array + Reconstructed load history. + + """ + # Number of heat sources + nSources = Q_b.shape[0] + # Time step sizes + dt = np.hstack((time[0], time[1:]-time[:-1])) + # Time vector + t = np.hstack((0., time, time[-1] + time[0])) + # Inverted time step sizes + dt_reconstructed = dt[::-1] + # Reconstructed time vector + t_reconstructed = np.hstack((0., np.cumsum(dt_reconstructed))) + # Accumulated heat extracted + f = np.hstack( + (np.zeros((nSources, 1), dtype=self.dtype), + np.cumsum(Q_b*dt, axis=1))) + f = np.hstack((f, f[:,-1:])) + # Create interpolation object for accumulated heat extracted + sf = interp1d(t, f, kind='linear', axis=1) + # Reconstructed load history + Q_reconstructed = (sf(t_reconstructed[1:]) - sf(t_reconstructed[:-1])) \ + / dt_reconstructed + + return Q_reconstructed + + def _check_inputs(self): + """ + This method ensures that the instances filled in the Solver object + are what is expected. + + """ + assert isinstance(self.boreholes, (list, Borefield)), \ + "Boreholes must be provided in a list." + assert len(self.boreholes) > 0, \ + "The list of boreholes is empty." + assert np.all([isinstance(b, Borehole) for b in self.boreholes]), \ + "The list of boreholes contains elements that are not Borehole " \ + "objects." + assert self.network is None or isinstance(self.network, Network), \ + "The network is not a valid 'Network' object." + if self.boundary_condition == 'MIFT': + assert not (self.m_flow_network is None and self.m_flow_borehole is None), \ + "The mass flow rate 'm_flow_borehole' or 'm_flow_network' must " \ + "be provided when using the 'MIFT' boundary condition." + assert not (self.m_flow_network is not None and self.m_flow_borehole is not None), \ + "Only one of 'm_flow_borehole' or 'm_flow_network' can " \ + "be provided when using the 'MIFT' boundary condition." + assert not self.cp_f is None, \ + "The heat capacity 'cp_f' must " \ + "be provided when using the 'MIFT' boundary condition." + assert not (type(self.m_flow_borehole) is np.ndarray and not np.size(self.m_flow_borehole, axis=1)==self.network.nInlets), \ + "The number of mass flow rates in 'm_flow_borehole' must " \ + "correspond to the number of circuits in the network." + assert type(self.time) is np.ndarray or isinstance(self.time, (float, np.floating)) or self.time is None, \ + "Time should be a float or an array." + # self.nSegments can now be an int or list + assert type(self.nBoreSegments) is list and len(self.nBoreSegments) == \ + len(self.nBoreSegments) and min(self.nBoreSegments) >= 1, \ + "The argument for number of segments `nSegments` should be " \ + "of type int or a list of integers. If passed as a list, the " \ + "length of the list should be equal to the number of boreholes" \ + "in the borefield. nSegments >= 1 is/are required." + acceptable_boundary_conditions = ['UHTR', 'UBWT', 'MIFT'] + assert type(self.boundary_condition) is str and self.boundary_condition in acceptable_boundary_conditions, \ + f"Boundary condition '{self.boundary_condition}' is not an " \ + f"acceptable boundary condition. \n" \ + f"Please provide one of the following inputs : " \ + f"{acceptable_boundary_conditions}" + assert type(self.approximate_FLS) is bool, \ + "The option 'approximate_FLS' should be set to True or False." + assert type(self.nFLS) is int and 1 <= self.nFLS <= 25, \ + "The option 'nFLS' should be a positive int and lower or equal to " \ + "25." + assert type(self.disp) is bool, \ + "The option 'disp' should be set to True or False." + assert type(self.profiles) is bool, \ + "The option 'profiles' should be set to True or False." + assert type(self.kind) is str, \ + "The option 'kind' should be set to a valid interpolation kind " \ + "in accordance with scipy.interpolate.interp1d options." + acceptable_dtypes = (np.single, np.double) + assert np.any([self.dtype is dtype for dtype in acceptable_dtypes]), \ + f"Data type '{self.dtype}' is not an acceptable data type. \n" \ + f"Please provide one of the following inputs : {acceptable_dtypes}" + # Check segment ratios + for j, (ratios, nSegments) in enumerate( + zip(self.segment_ratios, self.nBoreSegments)): + assert len(ratios) == nSegments, \ + f"The length of the segment ratios vectors must correspond to " \ + f"the number of segments, check borehole {j}." + error = np.abs(1. - np.sum(ratios)) + assert(error < 1.0e-6), \ + f"Defined segment ratios must add up to 1. " \ + f", check borehole {j}." + + return diff --git a/pygfunction/solvers/detailed.py b/pygfunction/solvers/detailed.py new file mode 100644 index 00000000..ec9113fe --- /dev/null +++ b/pygfunction/solvers/detailed.py @@ -0,0 +1,300 @@ +# -*- coding: utf-8 -*- +from time import perf_counter + +import numpy as np +from scipy.interpolate import interp1d as interp1d + +from ._base_solver import _BaseSolver +from ..heat_transfer import ( + finite_line_source, + finite_line_source_inclined_vectorized, + finite_line_source_vectorized + ) + + +class Detailed(_BaseSolver): + """ + Detailed solver for the evaluation of the g-function. + + This solver superimposes the finite line source (FLS) solution to + estimate the g-function of a geothermal bore field. Each borehole is + modeled as a series of finite line source segments, as proposed in + [#Detailed-CimBer2014]_. + + Parameters + ---------- + boreholes : list of Borehole objects + List of boreholes included in the bore field. + network : network object + Model of the network. + time : float or array + Values of time (in seconds) for which the g-function is evaluated. + boundary_condition : str + Boundary condition for the evaluation of the g-function. Should be one + of + + - 'UHTR' : + **Uniform heat transfer rate**. This corresponds to boundary + condition *BC-I* as defined by Cimmino and Bernier (2014) + [#Detailed-CimBer2014]_. + - 'UBWT' : + **Uniform borehole wall temperature**. This corresponds to + boundary condition *BC-III* as defined by Cimmino and Bernier + (2014) [#Detailed-CimBer2014]_. + - 'MIFT' : + **Mixed inlet fluid temperatures**. This boundary condition was + introduced by Cimmino (2015) [#Detailed-Cimmin2015]_ for + parallel-connected boreholes and extended to mixed + configurations by Cimmino (2019) [#Detailed-Cimmin2019]_. + + nSegments : int or list, optional + Number of line segments used per borehole, or list of number of + line segments used for each borehole. + Default is 8. + segment_ratios : array, list of arrays, or callable, optional + Ratio of the borehole length represented by each segment. The + sum of ratios must be equal to 1. The shape of the array is of + (nSegments,) or list of (nSegments[i],). If segment_ratios==None, + segments of equal lengths are considered. If a callable is provided, it + must return an array of size (nSegments,) when provided with nSegments + (of type int) as an argument, or an array of size (nSegments[i],) when + provided with an element of nSegments (of type list). + Default is :func:`utilities.segment_ratios`. + m_flow_borehole : (nInlets,) array or (nMassFlow, nInlets,) array, optional + Fluid mass flow rate into each circuit of the network. If a + (nMassFlow, nInlets,) array is supplied, the + (nMassFlow, nMassFlow,) variable mass flow rate g-functions + will be evaluated using the method of Cimmino (2024) + [#Detailed-Cimmin2024]_. Only required for the 'MIFT' boundary + condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be + provided. + Default is None. + m_flow_network : float or (nMassFlow,) array, optional + Fluid mass flow rate into the network of boreholes. If an array + is supplied, the (nMassFlow, nMassFlow,) variable mass flow + rate g-functions will be evaluated using the method of Cimmino + (2024) [#Detailed-Cimmin2024]_. Only required for the 'MIFT' boundary + condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be + provided. + Default is None. + cp_f : float, optional + Fluid specific isobaric heat capacity (in J/kg.degC). Only required + for the 'MIFT' boundary condition. + Default is None. + approximate_FLS : bool, optional + Set to true to use the approximation of the FLS solution of Cimmino + (2021) [#Detailed-Cimmin2021]_. This approximation does not require the + numerical evaluation of any integral. + Default is False. + nFLS : int, optional + Number of terms in the approximation of the FLS solution. This + parameter is unused if `approximate_FLS` is set to False. + Default is 10. Maximum is 25. + mQuad : int, optional + Number of Gauss-Legendre sample points for the integral over :math:`u` + in the inclined FLS solution. + Default is 11. + linear_threshold : float, optional + Threshold time (in seconds) under which the g-function is + linearized. The g-function value is then interpolated between 0 + and its value at the threshold. If linear_threshold==None, the + g-function is linearized for times + `t < r_b**2 / (25 * self.alpha)`. + Default is None. + disp : bool, optional + Set to true to print progression messages. + Default is False. + profiles : bool, optional + Set to true to keep in memory the temperatures and heat extraction + rates. + Default is False. + kind : string, optional + Interpolation method used for segment-to-segment thermal response + factors. See documentation for scipy.interpolate.interp1d. + Default is 'linear'. + dtype : numpy dtype, optional + numpy data type used for matrices and vectors. Should be one of + numpy.single or numpy.double. + Default is numpy.double. + + References + ---------- + .. [#Detailed-CimBer2014] Cimmino, M., & Bernier, M. (2014). A + semi-analytical method to generate g-functions for geothermal bore + fields. International Journal of Heat and Mass Transfer, 70, 641-650. + .. [#Detailed-Cimmin2015] Cimmino, M. (2015). The effects of borehole + thermal resistances and fluid flow rate on the g-functions of geothermal + bore fields. International Journal of Heat and Mass Transfer, 91, + 1119-1127. + .. [#Detailed-Cimmin2019] Cimmino, M. (2019). Semi-analytical method for + g-function calculation of bore fields with series- and + parallel-connected boreholes. Science and Technology for the Built + Environment, 25 (8), 1007-1022. + .. [#Detailed-Cimmin2021] Cimmino, M. (2021). An approximation of the + finite line source solution to model thermal interactions between + geothermal boreholes. International Communications in Heat and Mass + Transfer, 127, 105496. + .. [#Detailed-Cimmin2024] Cimmino, M. (2024). g-Functions for fields of + series- and parallel-connected boreholes with variable fluid mass flow + rate and reversible flow direction. Renewable Energy, 228, 120661. + + """ + def initialize(self, **kwargs): + """ + Split boreholes into segments. + + Returns + ------- + nSources : int + Number of finite line heat sources in the borefield used to + initialize the matrix of segment-to-segment thermal response + factors (of size: nSources x nSources). + + """ + # Split boreholes into segments + self.boreSegments = self.borehole_segments() + nSources = len(self.boreSegments) + return nSources + + def thermal_response_factors(self, time, alpha, kind='linear'): + """ + Evaluate the segment-to-segment thermal response factors for all pairs + of segments in the borefield at all time steps using the finite line + source solution. + + This method returns a scipy.interpolate.interp1d object of the matrix + of thermal response factors, containing a copy of the matrix accessible + by h_ij.y[:nSources,:nSources,:nt+1]. The first index along the + third axis corresponds to time t=0. The interp1d object can be used to + obtain thermal response factors at any intermediate time by + h_ij(t)[:nSources,:nSources]. + + Attributes + ---------- + time : float or array + Values of time (in seconds) for which the g-function is evaluated. + alpha : float + Soil thermal diffusivity (in m2/s). + kind : string, optional + Interpolation method used for segment-to-segment thermal response + factors. See documentation for scipy.interpolate.interp1d. + Default is 'linear'. + + Returns + ------- + h_ij : interp1d + interp1d object (scipy.interpolate) of the matrix of + segment-to-segment thermal response factors. + + """ + if self.disp: + print('Calculating segment to segment response factors ...', + end='') + # Number of time values + nt = len(np.atleast_1d(time)) + # Initialize chrono + tic = perf_counter() + # Initialize segment-to-segment response factors + h_ij = np.zeros((self.nSources, self.nSources, nt+1), dtype=self.dtype) + nBoreholes = len(self.boreholes) + segment_lengths = self.segment_lengths() + + # --------------------------------------------------------------------- + # Segment-to-segment thermal response factors for same-borehole + # thermal interactions + # --------------------------------------------------------------------- + h, i_segment, j_segment = \ + self._thermal_response_factors_borehole_to_self(time, alpha) + # Broadcast values to h_ij matrix + h_ij[j_segment, i_segment, 1:] = h + # --------------------------------------------------------------------- + # Segment-to-segment thermal response factors for + # borehole-to-borehole thermal interactions + # --------------------------------------------------------------------- + for i, (i0, i1) in enumerate(zip(self._i0Segments, self._i1Segments)): + # Segments of the receiving borehole + b2 = self.boreSegments[i0:i1] + if i+1 < nBoreholes: + # Segments of the emitting borehole + b1 = self.boreSegments[i1:] + h = finite_line_source( + time, alpha, b1, b2, approximation=self.approximate_FLS, + N=self.nFLS, M=self.mQuad) + # Broadcast values to h_ij matrix + h_ij[i0:i1, i1:, 1:] = h + h_ij[i1:, i0:i1, 1:] = \ + np.swapaxes(h, 0, 1) * np.divide.outer( + segment_lengths[i0:i1], + segment_lengths[i1:]).T[:,:,np.newaxis] + + # Return 2d array if time is a scalar + if np.isscalar(time): + h_ij = h_ij[:,:,1] + + # Interp1d object for thermal response factors + h_ij = interp1d(np.hstack((0., time)), h_ij, + kind=kind, copy=True, axis=2) + toc = perf_counter() + if self.disp: print(f' {toc - tic:.3f} sec') + + return h_ij + + def _thermal_response_factors_borehole_to_self(self, time, alpha): + """ + Evaluate the segment-to-segment thermal response factors for all pairs + of segments between each borehole and itself. + + Attributes + ---------- + time : float or array + Values of time (in seconds) for which the g-function is evaluated. + alpha : float + Soil thermal diffusivity (in m2/s). + + Returns + ------- + h : array + Finite line source solution. + i_segment : list + Indices of the emitting segments in the bore field. + j_segment : list + Indices of the receiving segments in the bore field. + """ + # Indices of the thermal response factors into h_ij + i_segment = np.concatenate( + [np.repeat(np.arange(i0, i1), nSegments) + for i0, i1, nSegments in zip( + self._i0Segments, self._i1Segments, self.nBoreSegments) + ]) + j_segment = np.concatenate( + [np.tile(np.arange(i0, i1), nSegments) + for i0, i1, nSegments in zip( + self._i0Segments, self._i1Segments, self.nBoreSegments) + ]) + # Unpack parameters + x = np.array([b.x for b in self.boreSegments]) + y = np.array([b.y for b in self.boreSegments]) + H = np.array([b.H for b in self.boreSegments]) + D = np.array([b.D for b in self.boreSegments]) + r_b = np.array([b.r_b for b in self.boreSegments]) + # Distances between boreholes + dis = np.maximum( + np.sqrt((x[i_segment] - x[j_segment])**2 + (y[i_segment] - y[j_segment])**2), + r_b[i_segment]) + # FLS solution + if np.all([b.is_vertical() for b in self.boreholes]): + h = finite_line_source_vectorized( + time, alpha, + dis, H[i_segment], D[i_segment], H[j_segment], D[j_segment], + approximation=self.approximate_FLS, N=self.nFLS) + else: + tilt = np.array([b.tilt for b in self.boreSegments]) + orientation = np.array([b.orientation for b in self.boreSegments]) + h = finite_line_source_inclined_vectorized( + time, alpha, + r_b[i_segment], x[i_segment], y[i_segment], H[i_segment], + D[i_segment], tilt[i_segment], orientation[i_segment], + x[j_segment], y[j_segment], H[j_segment], D[j_segment], + tilt[j_segment], orientation[j_segment], M=self.mQuad, + approximation=self.approximate_FLS, N=self.nFLS) + return h, i_segment, j_segment diff --git a/pygfunction/solvers/equivalent.py b/pygfunction/solvers/equivalent.py new file mode 100644 index 00000000..f23943c9 --- /dev/null +++ b/pygfunction/solvers/equivalent.py @@ -0,0 +1,846 @@ +# -*- coding: utf-8 -*- +from time import perf_counter + +import numpy as np +from scipy.cluster.hierarchy import ( + cut_tree, + dendrogram, + linkage + ) +from scipy.interpolate import interp1d as interp1d + +from ._base_solver import _BaseSolver +from ..boreholes import _EquivalentBorehole +from ..heat_transfer import ( + finite_line_source, + finite_line_source_equivalent_boreholes_vectorized, + finite_line_source_vectorized + ) +from ..networks import _EquivalentNetwork + + +class Equivalent(_BaseSolver): + """ + Equivalent solver for the evaluation of the g-function. + + This solver uses hierarchical agglomerative clustering to identify groups + of boreholes that are expected to have similar borehole wall temperatures + and heat extraction rates, as proposed by Prieto and Cimmino (2021) + [#Equivalent-PriCim2021]_. Each group of boreholes is represented by a + single equivalent borehole. The FLS solution is adapted to evaluate + thermal interactions between groups of boreholes. This greatly reduces + the number of evaluations of the FLS solution and the size of the system of + equations to evaluate the g-function. + + Parameters + ---------- + boreholes : list of Borehole objects + List of boreholes included in the bore field. + network : network object + Model of the network. + time : float or array + Values of time (in seconds) for which the g-function is evaluated. + boundary_condition : str + Boundary condition for the evaluation of the g-function. Should be one + of + + - 'UHTR' : + **Uniform heat transfer rate**. This corresponds to boundary + condition *BC-I* as defined by Cimmino and Bernier (2014) + [#Equivalent-CimBer2014]_. + - 'UBWT' : + **Uniform borehole wall temperature**. This corresponds to + boundary condition *BC-III* as defined by Cimmino and Bernier + (2014) [#Equivalent-CimBer2014]_. + - 'MIFT' : + **Mixed inlet fluid temperatures**. This boundary condition was + introduced by Cimmino (2015) [#Equivalent-Cimmin2015]_ for + parallel-connected boreholes and extended to mixed + configurations by Cimmino (2019) [#Equivalent-Cimmin2019]_. + + nSegments : int or list, optional + Number of line segments used per borehole, or list of number of + line segments used for each borehole. + Default is 8. + segment_ratios : array, list of arrays, or callable, optional + Ratio of the borehole length represented by each segment. The + sum of ratios must be equal to 1. The shape of the array is of + (nSegments,) or list of (nSegments[i],). If segment_ratios==None, + segments of equal lengths are considered. If a callable is provided, it + must return an array of size (nSegments,) when provided with nSegments + (of type int) as an argument, or an array of size (nSegments[i],) when + provided with an element of nSegments (of type list). + Default is :func:`utilities.segment_ratios`. + m_flow_borehole : (nInlets,) array or (nMassFlow, nInlets,) array, optional + Fluid mass flow rate into each circuit of the network. If a + (nMassFlow, nInlets,) array is supplied, the + (nMassFlow, nMassFlow,) variable mass flow rate g-functions + will be evaluated using the method of Cimmino (2024) + [#Equivalent-Cimmin2024]_. Only required for the 'MIFT' boundary + condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be + provided. + Default is None. + m_flow_network : float or (nMassFlow,) array, optional + Fluid mass flow rate into the network of boreholes. If an array + is supplied, the (nMassFlow, nMassFlow,) variable mass flow + rate g-functions will be evaluated using the method of Cimmino + (2024) [#Equivalent-Cimmin2024]_. Only required for the 'MIFT' boundary + condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be + provided. + Default is None. + cp_f : float, optional + Fluid specific isobaric heat capacity (in J/kg.degC). Only required + for the 'MIFT' boundary condition. + Default is None. + approximate_FLS : bool, optional + Set to true to use the approximation of the FLS solution of Cimmino + (2021) [#Equivalent-Cimmin2021]_. This approximation does not require + the numerical evaluation of any integral. When using the 'equivalent' + solver, the approximation is only applied to the thermal response at + the borehole radius. Thermal interaction between boreholes is evaluated + using the FLS solution. + Default is False. + nFLS : int, optional + Number of terms in the approximation of the FLS solution. This + parameter is unused if `approximate_FLS` is set to False. + Default is 10. Maximum is 25. + mQuad : int, optional + Number of Gauss-Legendre sample points for the integral over :math:`u` + in the inclined FLS solution. + Default is 11. + linear_threshold : float, optional + Threshold time (in seconds) under which the g-function is + linearized. The g-function value is then interpolated between 0 + and its value at the threshold. If linear_threshold==None, the + g-function is linearized for times + `t < r_b**2 / (25 * self.alpha)`. + Default is None. + disp : bool, optional + Set to true to print progression messages. + Default is False. + profiles : bool, optional + Set to true to keep in memory the temperatures and heat extraction + rates. + Default is False. + kind : string, optional + Interpolation method used for segment-to-segment thermal response + factors. See documentation for scipy.interpolate.interp1d. + Default is 'linear'. + dtype : numpy dtype, optional + numpy data type used for matrices and vectors. Should be one of + numpy.single or numpy.double. + Default is numpy.double. + disTol : float, optional + Relative tolerance on radial distance. Two distances + (d1, d2) between two pairs of boreholes are considered equal if the + difference between the two distances (abs(d1-d2)) is below tolerance. + Default is 0.01. + tol : float, optional + Relative tolerance on length and depth. Two lengths H1, H2 + (or depths D1, D2) are considered equal if abs(H1 - H2)/H2 < tol. + Default is 1.0e-6. + kClusters : int, optional + Increment on the minimum number of equivalent boreholes determined by + cutting the dendrogram of the bore field given by the hierarchical + agglomerative clustering method. Increasing the value of this parameter + increases the accuracy of the method. + Default is 1. + + References + ---------- + .. [#Equivalent-CimBer2014] Cimmino, M., & Bernier, M. (2014). A + semi-analytical method to generate g-functions for geothermal bore + fields. International Journal of Heat and Mass Transfer, 70, 641-650. + .. [#Equivalent-Cimmin2015] Cimmino, M. (2015). The effects of borehole + thermal resistances and fluid flow rate on the g-functions of geothermal + bore fields. International Journal of Heat and Mass Transfer, 91, + 1119-1127. + .. [#Equivalent-Cimmin2019] Cimmino, M. (2019). Semi-analytical method + for g-function calculation of bore fields with series- and + parallel-connected boreholes. Science and Technology for the Built + Environment, 25 (8), 1007-1022. + .. [#Equivalent-PriCim2021] Prieto, C., & Cimmino, M. + (2021). Thermal interactions in large irregular fields of geothermal + boreholes: the method of equivalent borehole. Journal of Building + Performance Simulation, 14 (4), 446-460. + .. [#Equivalent-Cimmin2021] Cimmino, M. (2021). An approximation of the + finite line source solution to model thermal interactions between + geothermal boreholes. International Communications in Heat and Mass + Transfer, 127, 105496. + .. [#Equivalent-Cimmin2024] Cimmino, M. (2024). g-Functions for fields of + series- and parallel-connected boreholes with variable fluid mass flow + rate and reversible flow direction. Renewable Energy, 228, 120661. + + """ + def initialize(self, disTol=0.01, tol=1.0e-6, kClusters=1, **kwargs): + """ + Initialize paramteters. Identify groups for equivalent boreholes. + + Returns + ------- + nSources : int + Number of finite line heat sources in the borefield used to + initialize the matrix of segment-to-segment thermal response + factors (of size: nSources x nSources). + + """ + self.disTol = disTol + self.tol = tol + self.kClusters = kClusters + # Check the validity of inputs + self._check_solver_specific_inputs() + # Initialize groups for equivalent boreholes + nSources = self.find_groups() + self.nBoreSegments = [self.nBoreSegments[0]] * self.nEqBoreholes + self.segment_ratios = [self.segment_ratios[0]] * self.nEqBoreholes + self.boreSegments = self.borehole_segments() + self._i0Segments = [sum(self.nBoreSegments[0:i]) + for i in range(self.nEqBoreholes)] + self._i1Segments = [sum(self.nBoreSegments[0:(i + 1)]) + for i in range(self.nEqBoreholes)] + return nSources + + def thermal_response_factors(self, time, alpha, kind='linear'): + """ + Evaluate the segment-to-segment thermal response factors for all pairs + of segments in the borefield at all time steps using the finite line + source solution. + + This method returns a scipy.interpolate.interp1d object of the matrix + of thermal response factors, containing a copy of the matrix accessible + by h_ij.y[:nSources,:nSources,:nt+1]. The first index along the + third axis corresponds to time t=0. The interp1d object can be used to + obtain thermal response factors at any intermediate time by + h_ij(t)[:nSources,:nSources]. + + Parameters + ---------- + time : float or array + Values of time (in seconds) for which the g-function is evaluated. + alpha : float + Soil thermal diffusivity (in m2/s). + kind : string, optional + Interpolation method used for segment-to-segment thermal response + factors. See documentation for scipy.interpolate.interp1d. + Default is linear. + + Returns + ------- + h_ij : interp1d + interp1d object (scipy.interpolate) of the matrix of + segment-to-segment thermal response factors. + + """ + if self.disp: + print('Calculating segment to segment response factors ...', + end='') + # Number of time values + nt = len(np.atleast_1d(time)) + # Initialize chrono + tic = perf_counter() + # Initialize segment-to-segment response factors + h_ij = np.zeros((self.nSources, self.nSources, nt+1), dtype=self.dtype) + segment_lengths = self.segment_lengths() + + # --------------------------------------------------------------------- + # Segment-to-segment thermal response factors for borehole-to-borehole + # thermal interactions + # --------------------------------------------------------------------- + # Groups correspond to unique pairs of borehole dimensions + for pairs in self.borehole_to_borehole: + i, j = pairs[0] + # Prepare inputs to the FLS function + dis, wDis = self._find_unique_distances(self.dis, pairs) + H1, D1, H2, D2, i_pair, j_pair, k_pair = \ + self._map_axial_segment_pairs(i, j) + H1 = H1.reshape(1, -1) + H2 = H2.reshape(1, -1) + D1 = D1.reshape(1, -1) + D2 = D2.reshape(1, -1) + N2 = np.array( + [[self.boreholes[j].nBoreholes for (i, j) in pairs]]).T + # Evaluate FLS at all time steps + h = finite_line_source_equivalent_boreholes_vectorized( + time, alpha, dis, wDis, H1, D1, H2, D2, N2) + # Broadcast values to h_ij matrix + for k, (i, j) in enumerate(pairs): + i_segment = self._i0Segments[i] + i_pair + j_segment = self._i0Segments[j] + j_pair + h_ij[j_segment, i_segment, 1:] = h[k, k_pair, :] + if not i == j: + h_ij[i_segment, j_segment, 1:] = (h[k, k_pair, :].T \ + * segment_lengths[j_segment]/segment_lengths[i_segment]).T + + # --------------------------------------------------------------------- + # Segment-to-segment thermal response factors for same-borehole thermal + # interactions + # --------------------------------------------------------------------- + # Groups correspond to unique borehole dimensions + for group in self.borehole_to_self: + # Index of first borehole in group + i = group[0] + # Find segment-to-segment similarities + H1, D1, H2, D2, i_pair, j_pair, k_pair = \ + self._map_axial_segment_pairs(i, i) + # Evaluate FLS at all time steps + dis = self.boreholes[i].r_b + H1 = H1.reshape(1, -1) + H2 = H2.reshape(1, -1) + D1 = D1.reshape(1, -1) + D2 = D2.reshape(1, -1) + h = finite_line_source_vectorized( + time, alpha, dis, H1, D1, H2, D2, + approximation=self.approximate_FLS, N=self.nFLS) + # Broadcast values to h_ij matrix + for i in group: + i_segment = self._i0Segments[i] + i_pair + j_segment = self._i0Segments[i] + j_pair + h_ij[j_segment, i_segment, 1:] = \ + h_ij[j_segment, i_segment, 1:] + h[0, k_pair, :] + + # Return 2d array if time is a scalar + if np.isscalar(time): + h_ij = h_ij[:,:,1] + + # Interp1d object for thermal response factors + h_ij = interp1d(np.hstack((0., time)), h_ij, + kind=kind, copy=True, axis=2) + toc = perf_counter() + if self.disp: print(f' {toc - tic:.3f} sec') + + return h_ij + + def find_groups(self, tol=1e-6): + """ + Identify groups of boreholes that can be represented by a single + equivalent borehole for the calculation of the g-function. + + Hierarchical agglomerative clustering is applied to the superposed + steady-state finite line source solution (i.e. the steady-state + dimensionless borehole wall temperature due to a uniform heat + extraction equal for all boreholes). The number of clusters is + evaluated by cutting the dendrogram at the half-height of the longest + branch and incrementing the number of intercepted branches by the value + of the kClusters parameter. + + Parameters + ---------- + tol : float + Tolerance on the temperature to identify the maxiumum number of + equivalent boreholes. + Default is 1e-6. + + Returns + ------- + nSources : int + Number of heat sources in the bore field. + + """ + if self.disp: print('Identifying equivalent boreholes ...', end='') + # Initialize chrono + tic = perf_counter() + + # Temperature change of individual boreholes + self.nBoreholes = len(self.boreholes) + # Equivalent field formed by all boreholes + eqField = _EquivalentBorehole(self.boreholes) + if self.nBoreholes > 1: + # Spatial superposition of the steady-state FLS solution + data = np.sum(finite_line_source(np.inf, 1., self.boreholes, self.boreholes), axis=1).reshape(-1,1) + # Split boreholes into groups of same dimensions + unique_boreholes = self._find_unique_boreholes(self.boreholes) + # Initialize empty list of clusters + self.clusters = [] + self.nEqBoreholes = 0 + for group in unique_boreholes: + if len(group) > 1: + # Maximum temperature + maxTemp = np.max(data[group]) + # Hierarchical agglomerative clustering based on temperatures + clusterization = linkage(data[group], method='complete') + dcoord = np.array( + dendrogram(clusterization, no_plot=True)['dcoord']) + # Maximum number of clusters + # Height to cut each tree to obtain the minimum number of clusters + disLeft = dcoord[:,1] - dcoord[:,0] + disRight = dcoord[:,2] - dcoord[:,3] + if np.max(disLeft) >= np.max(disRight): + i = disLeft.argmax() + height = 0.5*(dcoord[i,1] + dcoord[i,0]) + else: + i = disRight.argmax() + height = 0.5*(dcoord[i,2] + dcoord[i,3]) + # Find the number of clusters and increment by kClusters + # Maximum number of clusters + nClustersMax = min(np.sum(dcoord[:,1] > tol*maxTemp) + 1, + len(group)) + # Optimal number of cluster + nClusters = np.max( + cut_tree(clusterization, height=height)) + 1 + nClusters = min(nClusters + self.kClusters, nClustersMax) + # Cut the tree to find the borehole groups + clusters = cut_tree( + clusterization, n_clusters=nClusters) + self.clusters = self.clusters + \ + [label + self.nEqBoreholes for label in clusters] + else: + nClusters = 1 + self.clusters.append(self.nEqBoreholes) + self.nEqBoreholes += nClusters + else: + self.nEqBoreholes = self.nBoreholes + self.clusters = range(self.nBoreholes) + # Overwrite boreholes with equivalent boreholes + self.boreholes = [_EquivalentBorehole( + [borehole + for borehole, cluster in zip(self.boreholes, self.clusters) + if cluster==i]) + for i in range(self.nEqBoreholes)] + self.wBoreholes = np.array([b.nBoreholes for b in self.boreholes]) + # Find similar pairs of boreholes + self.borehole_to_self, self.borehole_to_borehole = \ + self._find_axial_borehole_pairs(self.boreholes) + # Store unique distances in the bore field + self.dis = eqField.unique_distance(eqField, self.disTol)[0][1:] + + if self.boundary_condition == 'MIFT': + pipes = [self.network.p[self.clusters.index(i)] + for i in range(self.nEqBoreholes)] + self.network = _EquivalentNetwork( + self.boreholes, + pipes, + nSegments=self.nBoreSegments[0], + segment_ratios=self.segment_ratios[0]) + + # Stop chrono + toc = perf_counter() + if self.disp: + print(f' {toc - tic:.3f} sec') + print(f'Calculations will be done using {self.nEqBoreholes} ' + f'equivalent boreholes') + + return self.nBoreSegments[0]*self.nEqBoreholes + + def segment_lengths(self): + """ + Return the length of all segments in the bore field. + + The segments lengths are used for the energy balance in the calculation + of the g-function. For equivalent boreholes, the length of segments + is multiplied by the number of boreholes in the group. + + Returns + ------- + H : array + Array of segment lengths (in m). + + """ + # Borehole lengths + H = np.array([seg.H*seg.nBoreholes + for (borehole, nSegments, ratios) in zip( + self.boreholes, + self.nBoreSegments, + self.segment_ratios) + for seg in borehole.segments( + nSegments, segment_ratios=ratios)], + dtype=self.dtype) + return H + + def _compare_boreholes(self, borehole1, borehole2): + """ + Compare two boreholes and checks if they have the same dimensions : + H, D, and r_b. + + Parameters + ---------- + borehole1 : Borehole object + First borehole. + borehole2 : Borehole object + Second borehole. + + Returns + ------- + similarity : bool + True if the two boreholes have the same dimensions. + + """ + # Compare lengths (H), buried depth (D) and radius (r_b) + if (abs((borehole1.H - borehole2.H)/borehole1.H) < self.tol and + abs((borehole1.r_b - borehole2.r_b)/borehole1.r_b) < self.tol and + abs((borehole1.D - borehole2.D)/(borehole1.D + 1e-30)) < self.tol): + similarity = True + else: + similarity = False + return similarity + + def _compare_real_pairs(self, pair1, pair2): + """ + Compare two pairs of boreholes or segments and return True if the two + pairs have the same FLS solution for real sources. + + Parameters + ---------- + pair1 : Tuple of Borehole objects + First pair of boreholes or segments. + pair2 : Tuple of Borehole objects + Second pair of boreholes or segments. + + Returns + ------- + similarity : bool + True if the two pairs have the same FLS solution. + + """ + deltaD1 = pair1[1].D - pair1[0].D + deltaD2 = pair2[1].D - pair2[0].D + + # Equality of lengths between pairs + cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol + and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol) + # Equality of lengths in each pair + equal_H = abs((pair1[0].H - pair1[1].H)/pair1[0].H) < self.tol + # Equality of buried depths differences + cond_deltaD = abs(deltaD1 - deltaD2)/abs(deltaD1 + 1e-30) < self.tol + # Equality of buried depths differences if all boreholes have the same + # length + cond_deltaD_equal_H = abs((abs(deltaD1) - abs(deltaD2))/(abs(deltaD1) + 1e-30)) < self.tol + if cond_H and (cond_deltaD or (equal_H and cond_deltaD_equal_H)): + similarity = True + else: + similarity = False + return similarity + + def _compare_image_pairs(self, pair1, pair2): + """ + Compare two pairs of boreholes or segments and return True if the two + pairs have the same FLS solution for mirror sources. + + Parameters + ---------- + pair1 : Tuple of Borehole objects + First pair of boreholes or segments. + pair2 : Tuple of Borehole objects + Second pair of boreholes or segments. + + Returns + ------- + similarity : bool + True if the two pairs have the same FLS solution. + + """ + sumD1 = pair1[1].D + pair1[0].D + sumD2 = pair2[1].D + pair2[0].D + + # Equality of lengths between pairs + cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol + and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol) + # Equality of buried depths sums + cond_sumD = abs((sumD1 - sumD2)/(sumD1 + 1e-30)) < self.tol + if cond_H and cond_sumD: + similarity = True + else: + similarity = False + return similarity + + def _compare_realandimage_pairs(self, pair1, pair2): + """ + Compare two pairs of boreholes or segments and return True if the two + pairs have the same FLS solution for both real and mirror sources. + + Parameters + ---------- + pair1 : Tuple of Borehole objects + First pair of boreholes or segments. + pair2 : Tuple of Borehole objects + Second pair of boreholes or segments. + + Returns + ------- + similarity : bool + True if the two pairs have the same FLS solution. + + """ + if (self._compare_real_pairs(pair1, pair2) + and self._compare_image_pairs(pair1, pair2)): + similarity = True + else: + similarity = False + return similarity + + def _find_axial_borehole_pairs(self, boreholes): + """ + Find axial (i.e. disregarding the radial distance) similarities between + borehole pairs to simplify the evaluation of the FLS solution. + + Parameters + ---------- + boreholes : list of Borehole objects + Boreholes in the bore field. + + Returns + ------- + borehole_to_self : list + Lists of borehole indexes for each unique set of borehole + dimensions (H, D, r_b) in the bore field. + borehole_to_borehole : list + Lists of tuples of borehole indexes for each unique pair of + boreholes that share the same (pairwise) dimensions (H, D). + + """ + # Compare for the full (real + image) FLS solution + compare_pairs = self._compare_realandimage_pairs + + nBoreholes = len(boreholes) + borehole_to_self = [] + # Only check for similarities if there is more than one borehole + if nBoreholes > 1: + borehole_to_borehole = [] + for i, borehole_i in enumerate(boreholes): + # Compare the borehole to all known unique sets of dimensions + for k, borehole_set in enumerate(borehole_to_self): + m = borehole_set[0] + # Add the borehole to the group if a similar borehole is + # found + if self._compare_boreholes(borehole_i, boreholes[m]): + borehole_set.append(i) + break + else: + # If no similar boreholes are known, append the groups + borehole_to_self.append([i]) + # Note : The range is different from similarities since + # an equivalent borehole to itself includes borehole-to- + # borehole thermal interactions + for j, borehole_j in enumerate(boreholes[i:], start=i): + pair0 = (borehole_i, borehole_j) # pair + pair1 = (borehole_j, borehole_i) # reciprocal pair + # Compare pairs of boreholes to known unique pairs + for pairs in borehole_to_borehole: + m, n = pairs[0] + pair_ref = (boreholes[m], boreholes[n]) + # Add the pair (or the reciprocal pair) to a group + # if a similar one is found + if compare_pairs(pair0, pair_ref): + pairs.append((i, j)) + break + elif compare_pairs(pair1, pair_ref): + pairs.append((j, i)) + break + # If no similar pairs are known, append the groups + else: + borehole_to_borehole.append([(i, j)]) + else: + # Outputs for a single borehole + borehole_to_self = [[0]] + borehole_to_borehole = [[(0, 0)]] + return borehole_to_self, borehole_to_borehole + + def _find_unique_boreholes(self, boreholes): + """ + Find unique sets of dimensions (h, D, r_b) in the bore field. + + Parameters + ---------- + boreholes : list of Borehole objects + Boreholes in the bore field. + + Returns + ------- + unique_boreholes : list + List of list of borehole indices that correspond to unique + borehole dimensions (H, D, r_b). + + """ + unique_boreholes = [] + for i, borehole_1 in enumerate(boreholes): + for group in unique_boreholes: + borehole_2 = boreholes[group[0]] + # Add the borehole to a group if similar dimensions are found + if self._compare_boreholes(borehole_1, borehole_2): + group.append(i) + break + else: + # If no similar boreholes are known, append the groups + unique_boreholes.append([i]) + + return unique_boreholes + + def _find_unique_distances(self, dis, indices): + """ + Find the number of occurrences of each unique distances between pairs + of boreholes. + + Parameters + ---------- + dis : array + Array of unique distances (in meters) in the bore field. + indices : list + List of tuples of borehole indices. + + Returns + ------- + dis : array + Array of unique distances (in meters) in the bore field. + wDis : array + Array of number of occurences of each unique distance for each + pair of equivalent boreholes in indices. + + """ + wDis = np.zeros((len(dis), len(indices)), dtype=np.uint) + for k, pair in enumerate(indices): + i, j = pair + b1, b2 = self.boreholes[i], self.boreholes[j] + # Generate a flattened array of distances between boreholes i and j + if not i == j: + dis_ij = b1.distance(b2).flatten() + else: + # Remove the borehole radius from the distances + dis_ij = b1.distance(b2)[ + ~np.eye(b1.nBoreholes, dtype=bool)].flatten() + wDis_ij = np.zeros(len(dis), dtype=np.uint) + # Get insert positions for the distances + iDis = np.searchsorted(dis, dis_ij, side='left') + # Find indexes where previous index is closer + prev_iDis_is_less = ((iDis == len(dis))|(np.fabs(dis_ij - dis[np.maximum(iDis-1, 0)]) < np.fabs(dis_ij - dis[np.minimum(iDis, len(dis)-1)]))) + iDis[prev_iDis_is_less] -= 1 + np.add.at(wDis_ij, iDis, 1) + wDis[:,k] = wDis_ij + + return dis.reshape((1, -1)), wDis + + def _map_axial_segment_pairs(self, iBor, jBor, + reaSource=True, imgSource=True): + """ + Find axial (i.e. disregarding the radial distance) similarities between + segment pairs along two boreholes to simplify the evaluation of the + FLS solution. + + The returned H1, D1, H2, and D2 can be used to evaluate the segment-to- + segment response factors using scipy.integrate.quad_vec. + + Parameters + ---------- + iBor : int + Index of the first borehole. + jBor : int + Index of the second borehole. + + Returns + ------- + H1 : float + Length of the emitting segments. + D1 : array + Array of buried depths of the emitting segments. + H2 : float + Length of the receiving segments. + D2 : array + Array of buried depths of the receiving segments. + i_pair : list + Indices of the emitting segments along a borehole. + j_pair : list + Indices of the receiving segments along a borehole. + k_pair : list + Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions + corresponding to all pairs in (i_pair, j_pair). + + """ + # Initialize local variables + borehole1 = self.boreholes[iBor] + borehole2 = self.boreholes[jBor] + assert reaSource or imgSource, \ + "At least one of reaSource and imgSource must be True." + if reaSource and imgSource: + # Find segment pairs for the full (real + image) FLS solution + compare_pairs = self._compare_realandimage_pairs + elif reaSource: + # Find segment pairs for the real FLS solution + compare_pairs = self._compare_real_pairs + elif imgSource: + # Find segment pairs for the image FLS solution + compare_pairs = self._compare_image_pairs + # Dive both boreholes into segments + segments1 = borehole1.segments( + self.nBoreSegments[iBor], segment_ratios=self.segment_ratios[iBor]) + segments2 = borehole2.segments( + self.nBoreSegments[jBor], segment_ratios=self.segment_ratios[jBor]) + # Prepare lists of segment lengths + H1 = [] + H2 = [] + # Prepare lists of segment buried depths + D1 = [] + D2 = [] + # All possible pairs (i, j) of indices between segments + i_pair = np.repeat(np.arange(self.nBoreSegments[iBor], dtype=np.uint), + self.nBoreSegments[jBor]) + j_pair = np.tile(np.arange(self.nBoreSegments[jBor], dtype=np.uint), + self.nBoreSegments[iBor]) + # Empty list of indices for unique pairs + k_pair = np.empty(self.nBoreSegments[iBor] * self.nBoreSegments[jBor], + dtype=np.uint) + unique_pairs = [] + nPairs = 0 + + p = 0 + for i, segment_i in enumerate(segments1): + for j, segment_j in enumerate(segments2): + pair = (segment_i, segment_j) + # Compare the segment pairs to all known unique pairs + for k, pair_k in enumerate(unique_pairs): + m, n = pair_k[0], pair_k[1] + pair_ref = (segments1[m], segments2[n]) + # Stop if a similar pair is found and assign the index + if compare_pairs(pair, pair_ref): + k_pair[p] = k + break + # If no similar pair is found : add a new pair, increment the + # number of unique pairs, and extract the associated buried + # depths + else: + k_pair[p] = nPairs + H1.append(segment_i.H) + H2.append(segment_j.H) + D1.append(segment_i.D) + D2.append(segment_j.D) + unique_pairs.append((i, j)) + nPairs += 1 + p += 1 + return np.array(H1), np.array(D1), np.array(H2), np.array(D2), i_pair, j_pair, k_pair + + def _check_solver_specific_inputs(self): + """ + This method ensures that solver specific inputs to the Solver object + are what is expected. + + """ + assert type(self.disTol) is float and self.disTol > 0., \ + "The distance tolerance 'disTol' should be a positive float." + assert type(self.tol) is float and self.tol > 0., \ + "The relative tolerance 'tol' should be a positive float." + assert type(self.kClusters) is int and self.kClusters >= 0, \ + "The precision increment 'kClusters' should be a positive int." + assert np.all(np.array(self.nBoreSegments, dtype=np.uint) == self.nBoreSegments[0]), \ + "Solver 'equivalent' can only handle equal numbers of segments." + assert np.all([np.allclose(segment_ratios, self.segment_ratios[0]) for segment_ratios in self.segment_ratios]), \ + "Solver 'equivalent' can only handle identical segment_ratios for all boreholes." + assert not np.any([b.is_tilted() for b in self.boreholes]), \ + "Solver 'equivalent' can only handle vertical boreholes." + if self.boundary_condition == 'MIFT': + assert np.all(np.array(self.network.c, dtype=int) == -1), \ + "Solver 'equivalent' is only valid for parallel-connected " \ + "boreholes." + assert (self.m_flow_borehole is None + or (self.m_flow_borehole.ndim==1 and np.allclose(self.m_flow_borehole, self.m_flow_borehole[0])) + or (self.m_flow_borehole.ndim==2 and np.all([np.allclose(self.m_flow_borehole[:, i], self.m_flow_borehole[0, i]) for i in range(self.nBoreholes)]))), \ + "Mass flow rates into the network must be equal for all " \ + "boreholes." + # Use the total network mass flow rate. + if (type(self.network.m_flow_network) is np.ndarray and \ + len(self.network.m_flow_network)==len(self.network.b)): + self.network.m_flow_network = \ + self.network.m_flow_network[0]*len(self.network.b) + # Verify that all boreholes have the same piping configuration + # This is best done by comparing the matrix of thermal resistances. + assert np.all( + [np.allclose(self.network.p[0]._Rd, pipe._Rd) + for pipe in self.network.p]), \ + "All boreholes must have the same piping configuration." + return diff --git a/pygfunction/solvers/similarities.py b/pygfunction/solvers/similarities.py new file mode 100644 index 00000000..887d1444 --- /dev/null +++ b/pygfunction/solvers/similarities.py @@ -0,0 +1,1334 @@ +# -*- coding: utf-8 -*- +from time import perf_counter + +import numpy as np +from scipy.interpolate import interp1d as interp1d + +from ._base_solver import _BaseSolver +from ..heat_transfer import ( + finite_line_source_inclined_vectorized, + finite_line_source_vectorized + ) + + +class Similarities(_BaseSolver): + """ + Similarities solver for the evaluation of the g-function. + + This solver superimposes the finite line source (FLS) solution to + estimate the g-function of a geothermal bore field. Each borehole is + modeled as a series of finite line source segments, as proposed in + [#Similarities-CimBer2014]_. The number of evaluations of the FLS solution + is decreased by identifying similar pairs of boreholes, for which the same + FLS value can be applied [#Similarities-Cimmin2018]_. + + Parameters + ---------- + boreholes : list of Borehole objects + List of boreholes included in the bore field. + network : network object + Model of the network. + time : float or array + Values of time (in seconds) for which the g-function is evaluated. + boundary_condition : str + Boundary condition for the evaluation of the g-function. Should be one + of + + - 'UHTR' : + **Uniform heat transfer rate**. This corresponds to boundary + condition *BC-I* as defined by Cimmino and Bernier (2014) + [#Similarities-CimBer2014]_. + - 'UBWT' : + **Uniform borehole wall temperature**. This corresponds to + boundary condition *BC-III* as defined by Cimmino and Bernier + (2014) [#Similarities-CimBer2014]_. + - 'MIFT' : + **Mixed inlet fluid temperatures**. This boundary condition was + introduced by Cimmino (2015) [#Similarities-Cimmin2015]_ for + parallel-connected boreholes and extended to mixed + configurations by Cimmino (2019) [#Similarities-Cimmin2019]_. + + nSegments : int or list, optional + Number of line segments used per borehole, or list of number of + line segments used for each borehole. + Default is 8. + segment_ratios : array, list of arrays, or callable, optional + Ratio of the borehole length represented by each segment. The + sum of ratios must be equal to 1. The shape of the array is of + (nSegments,) or list of (nSegments[i],). If segment_ratios==None, + segments of equal lengths are considered. If a callable is provided, it + must return an array of size (nSegments,) when provided with nSegments + (of type int) as an argument, or an array of size (nSegments[i],) when + provided with an element of nSegments (of type list). + Default is :func:`utilities.segment_ratios`. + m_flow_borehole : (nInlets,) array or (nMassFlow, nInlets,) array, optional + Fluid mass flow rate into each circuit of the network. If a + (nMassFlow, nInlets,) array is supplied, the + (nMassFlow, nMassFlow,) variable mass flow rate g-functions + will be evaluated using the method of Cimmino (2024) + [#Similarities-Cimmin2024]_. Only required for the 'MIFT' boundary + condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be + provided. + Default is None. + m_flow_network : float or (nMassFlow,) array, optional + Fluid mass flow rate into the network of boreholes. If an array + is supplied, the (nMassFlow, nMassFlow,) variable mass flow + rate g-functions will be evaluated using the method of Cimmino + (2024) [#Similarities-Cimmin2024]_. Only required for the 'MIFT' boundary + condition. Only one of 'm_flow_borehole' and 'm_flow_network' can be + provided. + Default is None. + cp_f : float, optional + Fluid specific isobaric heat capacity (in J/kg.degC). Only required + for the 'MIFT' boundary condition. + Default is None. + approximate_FLS : bool, optional + Set to true to use the approximation of the FLS solution of Cimmino + (2021) [#Similarities-Cimmin2021]_. This approximation does not require + the numerical evaluation of any integral. + Default is False. + nFLS : int, optional + Number of terms in the approximation of the FLS solution. This + parameter is unused if `approximate_FLS` is set to False. + Default is 10. Maximum is 25. + mQuad : int, optional + Number of Gauss-Legendre sample points for the integral over :math:`u` + in the inclined FLS solution. + Default is 11. + linear_threshold : float, optional + Threshold time (in seconds) under which the g-function is + linearized. The g-function value is then interpolated between 0 + and its value at the threshold. If linear_threshold==None, the + g-function is linearized for times + `t < r_b**2 / (25 * self.alpha)`. + Default is None. + disp : bool, optional + Set to true to print progression messages. + Default is False. + profiles : bool, optional + Set to true to keep in memory the temperatures and heat extraction + rates. + Default is False. + kind : string, optional + Interpolation method used for segment-to-segment thermal response + factors. See documentation for scipy.interpolate.interp1d. + Default is 'linear'. + dtype : numpy dtype, optional + numpy data type used for matrices and vectors. Should be one of + numpy.single or numpy.double. + Default is numpy.double. + disTol : float, optional + Relative tolerance on radial distance. Two distances + (d1, d2) between two pairs of boreholes are considered equal if the + difference between the two distances (abs(d1-d2)) is below tolerance. + Default is 0.01. + tol : float, optional + Relative tolerance on length and depth. Two lengths H1, H2 + (or depths D1, D2) are considered equal if abs(H1 - H2)/H2 < tol. + Default is 1.0e-6. + + References + ---------- + .. [#Similarities-CimBer2014] Cimmino, M., & Bernier, M. (2014). A + semi-analytical method to generate g-functions for geothermal bore + fields. International Journal of Heat and Mass Transfer, 70, 641-650. + .. [#Similarities-Cimmin2015] Cimmino, M. (2015). The effects of borehole + thermal resistances and fluid flow rate on the g-functions of geothermal + bore fields. International Journal of Heat and Mass Transfer, 91, + 1119-1127. + .. [#Similarities-Cimmin2018] Cimmino, M. (2018). Fast calculation of the + g-functions of geothermal borehole fields using similarities in the + evaluation of the finite line source solution. Journal of Building + Performance Simulation, 11 (6), 655-668. + .. [#Similarities-Cimmin2019] Cimmino, M. (2019). Semi-analytical method + for g-function calculation of bore fields with series- and + parallel-connected boreholes. Science and Technology for the Built + Environment, 25 (8), 1007-1022. + .. [#Similarities-Cimmin2021] Cimmino, M. (2021). An approximation of the + finite line source solution to model thermal interactions between + geothermal boreholes. International Communications in Heat and Mass + Transfer, 127, 105496. + .. [#Similarities-Cimmin2024] Cimmino, M. (2024). g-Functions for fields of + series- and parallel-connected boreholes with variable fluid mass flow + rate and reversible flow direction. Renewable Energy, 228, 120661. + + """ + def initialize(self, disTol=0.01, tol=1.0e-6, **kwargs): + """ + Split boreholes into segments and identify similarities in the + borefield. + + Returns + ------- + nSources : int + Number of finite line heat sources in the borefield used to + initialize the matrix of segment-to-segment thermal response + factors (of size: nSources x nSources). + + """ + self.disTol = disTol + self.tol = tol + # Check the validity of inputs + self._check_solver_specific_inputs() + # Split boreholes into segments + self.boreSegments = self.borehole_segments() + # Initialize similarities + self.find_similarities() + return len(self.boreSegments) + + def thermal_response_factors(self, time, alpha, kind='linear'): + """ + Evaluate the segment-to-segment thermal response factors for all pairs + of segments in the borefield at all time steps using the finite line + source solution. + + This method returns a scipy.interpolate.interp1d object of the matrix + of thermal response factors, containing a copy of the matrix accessible + by h_ij.y[:nSources,:nSources,:nt+1]. The first index along the + third axis corresponds to time t=0. The interp1d object can be used to + obtain thermal response factors at any intermediate time by + h_ij(t)[:nSources,:nSources]. + + Attributes + ---------- + time : float or array + Values of time (in seconds) for which the g-function is evaluated. + alpha : float + Soil thermal diffusivity (in m2/s). + kind : string, optional + Interpolation method used for segment-to-segment thermal response + factors. See documentation for scipy.interpolate.interp1d. + Default is 'linear'. + + Returns + ------- + h_ij : interp1d + interp1d object (scipy.interpolate) of the matrix of + segment-to-segment thermal response factors. + + """ + if self.disp: + print('Calculating segment to segment response factors ...', + end='') + # Number of time values + nt = len(np.atleast_1d(time)) + # Initialize chrono + tic = perf_counter() + # Initialize segment-to-segment response factors + h_ij = np.zeros((self.nSources, self.nSources, nt+1), dtype=self.dtype) + + # --------------------------------------------------------------------- + # Segment-to-segment thermal response factors for same-borehole thermal + # interactions (vertical boreholes) + # --------------------------------------------------------------------- + # Evaluate FLS at all time steps + h, i_segment, j_segment, k_segment = \ + self._thermal_response_factors_borehole_to_self_vertical( + time, alpha) + # Broadcast values to h_ij matrix + h_ij[j_segment, i_segment, 1:] = h[k_segment, :] + # --------------------------------------------------------------------- + # Segment-to-segment thermal response factors for same-borehole thermal + # interactions (inclined boreholes) + # --------------------------------------------------------------------- + # Evaluate FLS at all time steps + h, i_segment, j_segment, k_segment = \ + self._thermal_response_factors_borehole_to_self_inclined( + time, alpha) + # Broadcast values to h_ij matrix + h_ij[j_segment, i_segment, 1:] = h[k_segment, :] + # --------------------------------------------------------------------- + # Segment-to-segment thermal response factors for borehole-to-borehole + # thermal interactions (vertical boreholes) + # --------------------------------------------------------------------- + for pairs, distances, distance_indices in zip( + self.borehole_to_borehole_vertical, + self.borehole_to_borehole_distances_vertical, + self.borehole_to_borehole_indices_vertical): + # Index of first borehole pair in group + i, j = pairs[0] + # Find segment-to-segment similarities + H1, D1, H2, D2, i_pair, j_pair, k_pair = \ + self._map_axial_segment_pairs_vertical(i, j) + # Locate thermal response factors in the h_ij matrix + i_segment, j_segment, k_segment, l_segment = \ + self._map_segment_pairs_vertical( + i_pair, j_pair, k_pair, pairs, distance_indices) + # Evaluate FLS at all time steps + dis = np.reshape(distances, (-1, 1)) + H1 = H1.reshape(1, -1) + H2 = H2.reshape(1, -1) + D1 = D1.reshape(1, -1) + D2 = D2.reshape(1, -1) + h = finite_line_source_vectorized( + time, alpha, dis, H1, D1, H2, D2, + approximation=self.approximate_FLS, N=self.nFLS) + # Broadcast values to h_ij matrix + h_ij[j_segment, i_segment, 1:] = h[l_segment, k_segment, :] + if (self._compare_boreholes(self.boreholes[j], self.boreholes[i]) and + self.nBoreSegments[i] == self.nBoreSegments[j] and + self._uniform_segment_ratios[i] and + self._uniform_segment_ratios[j]): + h_ij[i_segment, j_segment, 1:] = h[l_segment, k_segment, :] + else: + h_ij[i_segment, j_segment, 1:] = (h * H2.T / H1.T)[l_segment, k_segment, :] + # --------------------------------------------------------------------- + # Segment-to-segment thermal response factors for borehole-to-borehole + # thermal interactions (inclined boreholes) + # --------------------------------------------------------------------- + # Evaluate FLS at all time steps + h, hT, i_segment, j_segment, k_segment = \ + self._thermal_response_factors_borehole_to_borehole_inclined( + time, alpha) + # Broadcast values to h_ij matrix + h_ij[j_segment, i_segment, 1:] = h[k_segment, :] + h_ij[i_segment, j_segment, 1:] = hT[k_segment, :] + + # Return 2d array if time is a scalar + if np.isscalar(time): + h_ij = h_ij[:,:,1] + + # Interp1d object for thermal response factors + h_ij = interp1d( + np.hstack((0., time)), h_ij, + kind=kind, copy=True, assume_sorted=True, axis=2) + toc = perf_counter() + if self.disp: print(f' {toc - tic:.3f} sec') + + return h_ij + + def _thermal_response_factors_borehole_to_borehole_inclined( + self, time, alpha): + """ + Evaluate the segment-to-segment thermal response factors for all pairs + of inclined segments. + + Attributes + ---------- + time : float or array + Values of time (in seconds) for which the g-function is evaluated. + alpha : float + Soil thermal diffusivity (in m2/s). + + Returns + ------- + h : array + Finite line source solution. + hT : array + Reciprocal finite line source solution. + i_segment : list + Indices of the emitting segments in the bore field. + j_segment : list + Indices of the receiving segments in the bore field. + k_segment : list + Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions + corresponding to all pairs in (i_pair, j_pair) in the bore field. + """ + rb1 = np.array([]) + x1 = np.array([]) + y1 = np.array([]) + H1 = np.array([]) + D1 = np.array([]) + tilt1 = np.array([]) + orientation1 = np.array([]) + x2 = np.array([]) + y2 = np.array([]) + H2 = np.array([]) + D2 = np.array([]) + tilt2 = np.array([]) + orientation2 = np.array([]) + i_segment = np.array([], dtype=np.uint) + j_segment = np.array([], dtype=np.uint) + k_segment = np.array([], dtype=np.uint) + k0 = 0 + for pairs in self.borehole_to_borehole_inclined: + # Index of first borehole pair in group + i, j = pairs[0] + # Find segment-to-segment similarities + rb1_i, x1_i, y1_i, H1_i, D1_i, tilt1_i, orientation1_i, \ + x2_i, y2_i, H2_i, D2_i, tilt2_i, orientation2_i, \ + i_pair, j_pair, k_pair = \ + self._map_axial_segment_pairs_inclined(i, j) + # Locate thermal response factors in the h_ij matrix + i_segment_i, j_segment_i, k_segment_i = \ + self._map_segment_pairs_inclined(i_pair, j_pair, k_pair, pairs) + # Append lists + rb1 = np.append(rb1, rb1_i) + x1 = np.append(x1, x1_i) + y1 = np.append(y1, y1_i) + H1 = np.append(H1, H1_i) + D1 = np.append(D1, D1_i) + tilt1 = np.append(tilt1, tilt1_i) + orientation1 = np.append(orientation1, orientation1_i) + x2 = np.append(x2, x2_i) + y2 = np.append(y2, y2_i) + H2 = np.append(H2, H2_i) + D2 = np.append(D2, D2_i) + tilt2 = np.append(tilt2, tilt2_i) + orientation2 = np.append(orientation2, orientation2_i) + i_segment = np.append(i_segment, i_segment_i) + j_segment = np.append(j_segment, j_segment_i) + k_segment = np.append(k_segment, k_segment_i + k0) + k0 += len(k_pair) + # Evaluate FLS at all time steps + h = finite_line_source_inclined_vectorized( + time, alpha, rb1, x1, y1, H1, D1, tilt1, orientation1, + x2, y2, H2, D2, tilt2, orientation2, M=self.mQuad, + approximation=self.approximate_FLS, N=self.nFLS) + hT = (h.T * H2 / H1).T + return h, hT, i_segment, j_segment, k_segment + + def _thermal_response_factors_borehole_to_self_inclined(self, time, alpha): + """ + Evaluate the segment-to-segment thermal response factors for all pairs + of segments between each inclined borehole and itself. + + Attributes + ---------- + time : float or array + Values of time (in seconds) for which the g-function is evaluated. + alpha : float + Soil thermal diffusivity (in m2/s). + + Returns + ------- + h : array + Finite line source solution. + i_segment : list + Indices of the emitting segments in the bore field. + j_segment : list + Indices of the receiving segments in the bore field. + k_segment : list + Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions + corresponding to all pairs in (i_pair, j_pair) in the bore field. + """ + rb1 = np.array([]) + x1 = np.array([]) + y1 = np.array([]) + H1 = np.array([]) + D1 = np.array([]) + tilt1 = np.array([]) + orientation1 = np.array([]) + x2 = np.array([]) + y2 = np.array([]) + H2 = np.array([]) + D2 = np.array([]) + tilt2 = np.array([]) + orientation2 = np.array([]) + i_segment = np.array([], dtype=np.uint) + j_segment = np.array([], dtype=np.uint) + k_segment = np.array([], dtype=np.uint) + k0 = 0 + for group in self.borehole_to_self_inclined: + # Index of first borehole in group + i = group[0] + # Find segment-to-segment similarities + rb1_i, x1_i, y1_i, H1_i, D1_i, tilt1_i, orientation1_i, \ + x2_i, y2_i, H2_i, D2_i, tilt2_i, orientation2_i, \ + i_pair, j_pair, k_pair = \ + self._map_axial_segment_pairs_inclined(i, i) + # Locate thermal response factors in the h_ij matrix + i_segment_i, j_segment_i, k_segment_i = \ + self._map_segment_pairs_inclined( + i_pair, j_pair, k_pair, [(n, n) for n in group]) + # Append lists + rb1 = np.append(rb1, rb1_i) + x1 = np.append(x1, x1_i) + y1 = np.append(y1, y1_i) + H1 = np.append(H1, H1_i) + D1 = np.append(D1, D1_i) + tilt1 = np.append(tilt1, tilt1_i) + orientation1 = np.append(orientation1, orientation1_i) + x2 = np.append(x2, x2_i) + y2 = np.append(y2, y2_i) + H2 = np.append(H2, H2_i) + D2 = np.append(D2, D2_i) + tilt2 = np.append(tilt2, tilt2_i) + orientation2 = np.append(orientation2, orientation2_i) + i_segment = np.append(i_segment, i_segment_i) + j_segment = np.append(j_segment, j_segment_i) + k_segment = np.append(k_segment, k_segment_i + k0) + k0 += len(k_pair) + # Evaluate FLS at all time steps + h = finite_line_source_inclined_vectorized( + time, alpha, rb1, x1, y1, H1, D1, tilt1, orientation1, + x2, y2, H2, D2, tilt2, orientation2, M=self.mQuad, + approximation=self.approximate_FLS, N=self.nFLS) + return h, i_segment, j_segment, k_segment + + def _thermal_response_factors_borehole_to_self_vertical(self, time, alpha): + """ + Evaluate the segment-to-segment thermal response factors for all pairs + of segments between each vertical borehole and itself. + + Attributes + ---------- + time : float or array + Values of time (in seconds) for which the g-function is evaluated. + alpha : float + Soil thermal diffusivity (in m2/s). + + Returns + ------- + h : array + Finite line source solution. + i_segment : list + Indices of the emitting segments in the bore field. + j_segment : list + Indices of the receiving segments in the bore field. + k_segment : list + Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions + corresponding to all pairs in (i_pair, j_pair) in the bore field. + """ + H1 = np.array([]) + D1 = np.array([]) + H2 = np.array([]) + D2 = np.array([]) + dis = np.array([]) + i_segment = np.array([], dtype=np.uint) + j_segment = np.array([], dtype=np.uint) + k_segment = np.array([], dtype=np.uint) + k0 = 0 + for group in self.borehole_to_self_vertical: + # Index of first borehole in group + i = group[0] + # Find segment-to-segment similarities + H1_i, D1_i, H2_i, D2_i, i_pair, j_pair, k_pair = \ + self._map_axial_segment_pairs_vertical(i, i) + # Locate thermal response factors in the h_ij matrix + i_segment_i, j_segment_i, k_segment_i, l_segment_i = \ + self._map_segment_pairs_vertical( + i_pair, j_pair, k_pair, [(n, n) for n in group], [0]) + # Append lists + H1 = np.append(H1, H1_i) + D1 = np.append(D1, D1_i) + H2 = np.append(H2, H2_i) + D2 = np.append(D2, D2_i) + if len(self.borehole_to_self_vertical) > 1: + dis = np.append(dis, np.full(len(H1_i), self.boreholes[i].r_b)) + else: + dis = self.boreholes[i].r_b + i_segment = np.append(i_segment, i_segment_i) + j_segment = np.append(j_segment, j_segment_i) + k_segment = np.append(k_segment, k_segment_i + k0) + k0 += np.max(k_pair) + 1 + # Evaluate FLS at all time steps + h = finite_line_source_vectorized( + time, alpha, dis, H1, D1, H2, D2, + approximation=self.approximate_FLS, N=self.nFLS) + return h, i_segment, j_segment, k_segment + + def find_similarities(self): + """ + Find similarities in the FLS solution for groups of boreholes. + + This function identifies pairs of boreholes for which the evaluation + of the Finite Line Source (FLS) solution is equivalent. + + """ + if self.disp: print('Identifying similarities ...', end='') + # Initialize chrono + tic = perf_counter() + + # Find similar pairs of boreholes + # Boreholes can only be similar if their segments are similar + self.borehole_to_self_vertical, self.borehole_to_self_inclined, \ + self.borehole_to_borehole_vertical, self.borehole_to_borehole_inclined = \ + self._find_axial_borehole_pairs(self.boreholes) + # Find distances for each similar pairs of vertical boreholes + self.borehole_to_borehole_distances_vertical, self.borehole_to_borehole_indices_vertical = \ + self._find_distances( + self.boreholes, self.borehole_to_borehole_vertical) + + # Stop chrono + toc = perf_counter() + if self.disp: print(f' {toc - tic:.3f} sec') + + return + + def _compare_boreholes(self, borehole1, borehole2): + """ + Compare two boreholes and checks if they have the same dimensions : + H, D, r_b, and tilt. + + Parameters + ---------- + borehole1 : Borehole object + First borehole. + borehole2 : Borehole object + Second borehole. + + Returns + ------- + similarity : bool + True if the two boreholes have the same dimensions. + + """ + # Compare lengths (H), buried depth (D) and radius (r_b) + if (abs((borehole1.H - borehole2.H)/borehole1.H) < self.tol and + abs((borehole1.r_b - borehole2.r_b)/borehole1.r_b) < self.tol and + abs((borehole1.D - borehole2.D)/(borehole1.D + 1e-30)) < self.tol and + abs(abs(borehole1.tilt) - abs(borehole2.tilt))/(abs(borehole1.tilt) + 1e-30) < self.tol): + similarity = True + else: + similarity = False + return similarity + + def _compare_real_pairs_vertical(self, pair1, pair2): + """ + Compare two pairs of vertical boreholes or segments and return True if + the two pairs have the same FLS solution for real sources. + + Parameters + ---------- + pair1 : Tuple of Borehole objects + First pair of boreholes or segments. + pair2 : Tuple of Borehole objects + Second pair of boreholes or segments. + + Returns + ------- + similarity : bool + True if the two pairs have the same FLS solution. + + """ + deltaD1 = pair1[1].D - pair1[0].D + deltaD2 = pair2[1].D - pair2[0].D + + # Equality of lengths between pairs + cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol + and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol) + # Equality of lengths in each pair + equal_H = abs((pair1[0].H - pair1[1].H)/pair1[0].H) < self.tol + # Equality of buried depths differences + cond_deltaD = abs(deltaD1 - deltaD2)/abs(deltaD1 + 1e-30) < self.tol + # Equality of buried depths differences if all boreholes have the same + # length + cond_deltaD_equal_H = abs((abs(deltaD1) - abs(deltaD2))/(abs(deltaD1) + 1e-30)) < self.tol + if cond_H and (cond_deltaD or (equal_H and cond_deltaD_equal_H)): + similarity = True + else: + similarity = False + return similarity + + def _compare_image_pairs_vertical(self, pair1, pair2): + """ + Compare two pairs of vertical boreholes or segments and return True if + the two pairs have the same FLS solution for mirror sources. + + Parameters + ---------- + pair1 : Tuple of Borehole objects + First pair of boreholes or segments. + pair2 : Tuple of Borehole objects + Second pair of boreholes or segments. + + Returns + ------- + similarity : bool + True if the two pairs have the same FLS solution. + + """ + sumD1 = pair1[1].D + pair1[0].D + sumD2 = pair2[1].D + pair2[0].D + + # Equality of lengths between pairs + cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol + and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol) + # Equality of buried depths sums + cond_sumD = abs((sumD1 - sumD2)/(sumD1 + 1e-30)) < self.tol + if cond_H and cond_sumD: + similarity = True + else: + similarity = False + return similarity + + def _compare_realandimage_pairs_vertical(self, pair1, pair2): + """ + Compare two pairs of vertical boreholes or segments and return True if + the two pairs have the same FLS solution for both real and mirror + sources. + + Parameters + ---------- + pair1 : Tuple of Borehole objects + First pair of boreholes or segments. + pair2 : Tuple of Borehole objects + Second pair of boreholes or segments. + + Returns + ------- + similarity : bool + True if the two pairs have the same FLS solution. + + """ + if (self._compare_real_pairs_vertical(pair1, pair2) + and self._compare_image_pairs_vertical(pair1, pair2)): + similarity = True + else: + similarity = False + return similarity + + def _compare_real_pairs_inclined(self, pair1, pair2): + """ + Compare two pairs of inclined boreholes or segments and return True if + the two pairs have the same FLS solution for real sources. + + Parameters + ---------- + pair1 : Tuple of Borehole objects + First pair of boreholes or segments. + pair2 : Tuple of Borehole objects + Second pair of boreholes or segments. + + Returns + ------- + similarity : bool + True if the two pairs have the same FLS solution. + + """ + dx1 = pair1[0].x - pair1[1].x; dx2 = pair2[0].x - pair2[1].x + dy1 = pair1[0].y - pair1[1].y; dy2 = pair2[0].y - pair2[1].y + dis1 = np.sqrt(dx1**2 + dy1**2); dis2 = np.sqrt(dx2**2 + dy2**2) + theta_12_1 = np.arctan2(dy1, dx1); theta_12_2 = np.arctan2(dy2, dx2) + deltaD1 = pair1[0].D - pair1[1].D; deltaD2 = pair2[0].D - pair2[1].D + # Equality of lengths between pairs + cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol + and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol) + # Equality of buried depths differences + cond_deltaD = abs(deltaD1 - deltaD2)/(abs(deltaD1) + 1e-30) < self.tol + # Equality of distances + cond_dis = abs(dis1 - dis2)/(abs(dis1) + 1e-30) < self.disTol + # Equality of tilts + cond_beta = ( + abs(abs(pair1[0].tilt) - abs(pair2[0].tilt))/(abs(pair1[0].tilt) + 1e-30) < self.tol + and abs(abs(pair1[1].tilt) - abs(pair2[1].tilt))/(abs(pair1[1].tilt) + 1e-30) < self.tol) + # Equality of relative orientations + sin_b1_cos_dt1_1 = np.sin(pair1[0].tilt) * np.cos(theta_12_1 - pair1[0].orientation) + sin_b2_cos_dt2_1 = np.sin(pair1[1].tilt) * np.cos(theta_12_1 - pair1[1].orientation) + sin_b1_cos_dt1_2 = np.sin(pair2[0].tilt) * np.cos(theta_12_2 - pair2[0].orientation) + sin_b2_cos_dt2_2 = np.sin(pair2[1].tilt) * np.cos(theta_12_2 - pair2[1].orientation) + cond_theta = ( + abs(sin_b1_cos_dt1_1 - sin_b1_cos_dt1_2) / (abs(sin_b1_cos_dt1_1) + 1e-30) > self.tol + and abs(sin_b2_cos_dt2_1 - sin_b2_cos_dt2_2) / (abs(sin_b2_cos_dt2_1) + 1e-30) > self.tol) + if cond_H and cond_deltaD and cond_dis and cond_beta and cond_theta: + similarity = True + else: + similarity = False + return similarity + + def _compare_image_pairs_inclined(self, pair1, pair2): + """ + Compare two pairs of inclined boreholes or segments and return True if + the two pairs have the same FLS solution for mirror sources. + + Parameters + ---------- + pair1 : Tuple of Borehole objects + First pair of boreholes or segments. + pair2 : Tuple of Borehole objects + Second pair of boreholes or segments. + + Returns + ------- + similarity : bool + True if the two pairs have the same FLS solution. + + """ + dx1 = pair1[0].x - pair1[1].x; dx2 = pair2[0].x - pair2[1].x + dy1 = pair1[0].y - pair1[1].y; dy2 = pair2[0].y - pair2[1].y + dis1 = np.sqrt(dx1**2 + dy1**2); dis2 = np.sqrt(dx2**2 + dy2**2) + theta_12_1 = np.arctan2(dy1, dx1); theta_12_2 = np.arctan2(dy2, dx2) + sumD1 = pair1[0].D + pair1[1].D; sumD2 = pair2[0].D + pair2[1].D + # Equality of lengths between pairs + cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol + and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol) + # Equality of buried depths sums + cond_sumD = abs(sumD1 - sumD2)/(abs(sumD1) + 1e-30) < self.tol + # Equality of distances + cond_dis = abs(dis1 - dis2)/(abs(dis1) + 1e-30) < self.disTol + # Equality of tilts + cond_beta = ( + abs(abs(pair1[0].tilt) - abs(pair2[0].tilt))/(abs(pair1[0].tilt) + 1e-30) < self.tol + and abs(abs(pair1[1].tilt) - abs(pair2[1].tilt))/(abs(pair1[1].tilt) + 1e-30) < self.tol) + # Equality of relative orientations + sin_b1_cos_dt1_1 = np.sin(pair1[0].tilt) * np.cos(theta_12_1 - pair1[0].orientation) + sin_b2_cos_dt2_1 = np.sin(pair1[1].tilt) * np.cos(theta_12_1 - pair1[1].orientation) + sin_b1_cos_dt1_2 = np.sin(pair2[0].tilt) * np.cos(theta_12_2 - pair2[0].orientation) + sin_b2_cos_dt2_2 = np.sin(pair2[1].tilt) * np.cos(theta_12_2 - pair2[1].orientation) + cond_theta = ( + abs(sin_b1_cos_dt1_1 - sin_b1_cos_dt1_2) / (abs(sin_b1_cos_dt1_1) + 1e-30) > self.tol + and abs(sin_b2_cos_dt2_1 - sin_b2_cos_dt2_2) / (abs(sin_b2_cos_dt2_1) + 1e-30) > self.tol) + if cond_H and cond_sumD and cond_dis and cond_beta and cond_theta: + similarity = True + else: + similarity = False + return similarity + + def _compare_realandimage_pairs_inclined(self, pair1, pair2): + """ + Compare two pairs of inclined boreholes or segments and return True if + the two pairs have the same FLS solution for both real and mirror + sources. + + Parameters + ---------- + pair1 : Tuple of Borehole objects + First pair of boreholes or segments. + pair2 : Tuple of Borehole objects + Second pair of boreholes or segments. + + Returns + ------- + similarity : bool + True if the two pairs have the same FLS solution. + + Notes + ----- + For inclined boreholes the similarity condition is the same for real + and image parts of the solution. + + """ + dx1 = pair1[0].x - pair1[1].x; dx2 = pair2[0].x - pair2[1].x + dy1 = pair1[0].y - pair1[1].y; dy2 = pair2[0].y - pair2[1].y + dis1 = np.sqrt(dx1**2 + dy1**2); dis2 = np.sqrt(dx2**2 + dy2**2) + theta_12_1 = np.arctan2(dy1, dx1); theta_12_2 = np.arctan2(dy2, dx2) + # Equality of lengths between pairs + cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol + and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol) + # Equality of buried depths + cond_D = ( + abs(pair1[0].D - pair2[0].D)/(abs(pair1[0].D) + 1e-30) < self.tol + and abs(pair1[1].D - pair2[1].D)/(abs(pair1[1].D) + 1e-30) < self.tol) + # Equality of distances + cond_dis = abs(dis1 - dis2)/(abs(dis1) + 1e-30) < self.disTol + # Equality of tilts + cond_beta = ( + abs(abs(pair1[0].tilt) - abs(pair2[0].tilt))/(abs(pair1[0].tilt) + 1e-30) < self.tol + and abs(abs(pair1[1].tilt) - abs(pair2[1].tilt))/(abs(pair1[1].tilt) + 1e-30) < self.tol) + # Equality of relative orientations + sin_b1_cos_dt1_1 = np.sin(pair1[0].tilt) * np.cos(theta_12_1 - pair1[0].orientation) + sin_b2_cos_dt2_1 = np.sin(pair1[1].tilt) * np.cos(theta_12_1 - pair1[1].orientation) + sin_b1_cos_dt1_2 = np.sin(pair2[0].tilt) * np.cos(theta_12_2 - pair2[0].orientation) + sin_b2_cos_dt2_2 = np.sin(pair2[1].tilt) * np.cos(theta_12_2 - pair2[1].orientation) + cond_theta = ( + abs(sin_b1_cos_dt1_1 - sin_b1_cos_dt1_2) / (abs(sin_b1_cos_dt1_1) + 1e-30) < self.tol + and abs(sin_b2_cos_dt2_1 - sin_b2_cos_dt2_2) / (abs(sin_b2_cos_dt2_1) + 1e-30) < self.tol) + if cond_H and cond_D and cond_dis and cond_beta and cond_theta: + similarity = True + else: + similarity = False + return similarity + + def _find_axial_borehole_pairs(self, boreholes): + """ + Find axial (i.e. disregarding the radial distance) similarities between + borehole pairs to simplify the evaluation of the FLS solution. + + Parameters + ---------- + boreholes : list of Borehole objects + Boreholes in the bore field. + + Returns + ------- + borehole_to_self : list + Lists of borehole indexes for each unique set of borehole + dimensions (H, D, r_b) in the bore field. + borehole_to_borehole : list + Lists of tuples of borehole indexes for each unique pair of + boreholes that share the same (pairwise) dimensions (H, D). + + """ + nBoreholes = len(boreholes) + borehole_to_self_vertical = [] + borehole_to_self_inclined = [] + # Only check for similarities if there is more than one borehole + if nBoreholes > 1: + borehole_to_borehole_vertical = [] + borehole_to_borehole_inclined = [] + for i, (borehole_i, nSegments_i, ratios_i) in enumerate( + zip(boreholes, self.nBoreSegments, self.segment_ratios)): + # Compare the borehole to all known unique sets of dimensions + if borehole_i.is_vertical(): + borehole_to_self = borehole_to_self_vertical + compare_pairs = self._compare_realandimage_pairs_vertical + else: + borehole_to_self = borehole_to_self_inclined + compare_pairs = self._compare_realandimage_pairs_inclined + for k, borehole_set in enumerate(borehole_to_self): + m = borehole_set[0] + # Add the borehole to the group if a similar borehole is + # found + if (self._compare_boreholes(borehole_i, boreholes[m]) and + (self._equal_segment_ratios or + (nSegments_i == self.nBoreSegments[m] and + np.allclose(ratios_i, + self.segment_ratios[m], + rtol=self.tol)))): + borehole_set.append(i) + break + else: + # If no similar boreholes are known, append the groups + borehole_to_self.append([i]) + + for j, (borehole_j, nSegments_j, ratios_j) in enumerate( + zip(boreholes[i+1:], + self.nBoreSegments[i+1:], + self.segment_ratios[i+1:]), + start=i+1): + pair0 = (borehole_i, borehole_j) # pair + pair1 = (borehole_j, borehole_i) # reciprocal pair + # Compare pairs of boreholes to known unique pairs + if borehole_i.is_vertical() and borehole_j.is_vertical(): + borehole_to_borehole = borehole_to_borehole_vertical + compare_pairs = self._compare_realandimage_pairs_vertical + else: + borehole_to_borehole = borehole_to_borehole_inclined + compare_pairs = self._compare_realandimage_pairs_inclined + for pairs in borehole_to_borehole: + m, n = pairs[0] + pair_ref = (boreholes[m], boreholes[n]) + # Add the pair (or the reciprocal pair) to a group + # if a similar one is found + if (compare_pairs(pair0, pair_ref) and + (self._equal_segment_ratios or + (nSegments_i == self.nBoreSegments[m] and + nSegments_j == self.nBoreSegments[n] and + np.allclose(ratios_i, + self.segment_ratios[m], + rtol=self.tol) and + np.allclose(ratios_j, + self.segment_ratios[n], + rtol=self.tol)))): + pairs.append((i, j)) + break + elif (compare_pairs(pair1, pair_ref) and + (self._equal_segment_ratios or + (nSegments_j == self.nBoreSegments[m] and + nSegments_i == self.nBoreSegments[n] and + np.allclose(ratios_j, + self.segment_ratios[m], + rtol=self.tol) and + np.allclose(ratios_i, + self.segment_ratios[n], + rtol=self.tol)))): + pairs.append((j, i)) + break + # If no similar pairs are known, append the groups + else: + borehole_to_borehole.append([(i, j)]) + + else: + # Outputs for a single borehole + if boreholes[0].is_vertical: + borehole_to_self_vertical = [[0]] + borehole_to_self_inclined = [] + else: + borehole_to_self_vertical = [] + borehole_to_self_inclined = [[0]] + borehole_to_borehole_vertical = [] + borehole_to_borehole_inclined = [] + return borehole_to_self_vertical, borehole_to_self_inclined, \ + borehole_to_borehole_vertical, borehole_to_borehole_inclined + + def _find_distances(self, boreholes, borehole_to_borehole): + """ + Find unique distances between pairs of boreholes for each unique pair + of boreholes in the bore field. + + Parameters + ---------- + boreholes : list of Borehole objects + Boreholes in the bore field. + borehole_to_borehole : list + Lists of tuples of borehole indexes for each unique pair of + boreholes that share the same (pairwise) dimensions (H, D). + + Returns + ------- + borehole_to_borehole_distances : list + Sorted lists of borehole-to-borehole radial distances for each + unique pair of boreholes. + borehole_to_borehole_indices : list + Lists of indexes of distances associated with each borehole pair. + + """ + nGroups = len(borehole_to_borehole) + borehole_to_borehole_distances = [[] for i in range(nGroups)] + borehole_to_borehole_indices = \ + [np.empty(len(group), dtype=np.uint) for group in borehole_to_borehole] + # Find unique distances for each group + for i, (pairs, distances, distance_indices) in enumerate( + zip(borehole_to_borehole, + borehole_to_borehole_distances, + borehole_to_borehole_indices)): + nPairs = len(pairs) + # Array of all borehole-to-borehole distances within the group + all_distances = np.array( + [boreholes[pair[0]].distance(boreholes[pair[1]]) + for pair in pairs]) + # Indices to sort the distance array + i_sort = all_distances.argsort() + # Sort the distance array + distances_sorted = all_distances[i_sort] + j0 = 0 + j1 = 1 + nDis = 0 + # For each increasing distance in the sorted array : + # 1 - find all distances that are within tolerance + # 2 - add the average distance in the list of unique distances + # 3 - associate the distance index to all pairs for the identified + # distances + # 4 - re-start at the next distance index not yet accounted for. + while j0 < nPairs and j1 > 0: + # Find the first distance outside tolerance + j1 = np.argmax( + distances_sorted >= (1+self.disTol)*distances_sorted[j0]) + if j1 > j0: + # Average distance between pairs of boreholes + distances.append(np.mean(distances_sorted[j0:j1])) + # Apply distance index to borehole pairs + distance_indices[i_sort[j0:j1]] = nDis + else: + # Average distance between pairs of boreholes + distances.append(np.mean(distances_sorted[j0:])) + # Apply distance index to borehole pairs + distance_indices[i_sort[j0:]] = nDis + j0 = j1 + nDis += 1 + return borehole_to_borehole_distances, borehole_to_borehole_indices + + def _map_axial_segment_pairs_vertical( + self, i, j, reaSource=True, imgSource=True): + """ + Find axial (i.e. disregarding the radial distance) similarities between + segment pairs along two boreholes to simplify the evaluation of the + FLS solution. + + The returned H1, D1, H2, and D2 can be used to evaluate the segment-to- + segment response factors using scipy.integrate.quad_vec. + + Parameters + ---------- + i : int + Index of the first borehole. + j : int + Index of the second borehole. + + Returns + ------- + H1 : array + Length of the emitting segments. + D1 : array + Array of buried depths of the emitting segments. + H2 : array + Length of the receiving segments. + D2 : array + Array of buried depths of the receiving segments. + i_pair : list + Indices of the emitting segments along a borehole. + j_pair : list + Indices of the receiving segments along a borehole. + k_pair : list + Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions + corresponding to all pairs in (i_pair, j_pair). + + """ + # Initialize local variables + borehole1 = self.boreholes[i] + borehole2 = self.boreholes[j] + assert reaSource or imgSource, \ + "At least one of reaSource and imgSource must be True." + if reaSource and imgSource: + # Find segment pairs for the full (real + image) FLS solution + compare_pairs = self._compare_realandimage_pairs_vertical + elif reaSource: + # Find segment pairs for the real FLS solution + compare_pairs = self._compare_real_pairs_vertical + elif imgSource: + # Find segment pairs for the image FLS solution + compare_pairs = self._compare_image_pairs_vertical + # Dive both boreholes into segments + segments1 = borehole1.segments( + self.nBoreSegments[i], segment_ratios=self.segment_ratios[i]) + segments2 = borehole2.segments( + self.nBoreSegments[j], segment_ratios=self.segment_ratios[j]) + # Prepare lists of segment lengths + H1 = [] + H2 = [] + # Prepare lists of segment buried depths + D1 = [] + D2 = [] + # All possible pairs (i, j) of indices between segments + i_pair = np.repeat(np.arange(self.nBoreSegments[i], dtype=np.uint), + self.nBoreSegments[j]) + j_pair = np.tile(np.arange(self.nBoreSegments[j], dtype=np.uint), + self.nBoreSegments[i]) + # Empty list of indices for unique pairs + k_pair = np.empty(self.nBoreSegments[i] * self.nBoreSegments[j], + dtype=np.uint) + unique_pairs = [] + nPairs = 0 + + p = 0 + for ii, segment_i in enumerate(segments1): + for jj, segment_j in enumerate(segments2): + pair = (segment_i, segment_j) + # Compare the segment pairs to all known unique pairs + for k, pair_k in enumerate(unique_pairs): + m, n = pair_k[0], pair_k[1] + pair_ref = (segments1[m], segments2[n]) + # Stop if a similar pair is found and assign the index + if compare_pairs(pair, pair_ref): + k_pair[p] = k + break + # If no similar pair is found : add a new pair, increment the + # number of unique pairs, and extract the associated buried + # depths + else: + k_pair[p] = nPairs + H1.append(segment_i.H) + H2.append(segment_j.H) + D1.append(segment_i.D) + D2.append(segment_j.D) + unique_pairs.append((ii, jj)) + nPairs += 1 + p += 1 + return np.array(H1), np.array(D1), np.array(H2), np.array(D2), i_pair, j_pair, k_pair + + def _map_axial_segment_pairs_inclined( + self, i, j, reaSource=True, imgSource=True): + """ + Find axial similarities between segment pairs along two boreholes to + simplify the evaluation of the FLS solution. + + The returned H1, D1, H2, and D2 can be used to evaluate the segment-to- + segment response factors using scipy.integrate.quad_vec. + + Parameters + ---------- + i : int + Index of the first borehole. + j : int + Index of the second borehole. + + Returns + ------- + rb1 : array + Radii of the emitting heat sources. + x1 : array + x-Positions of the emitting heat sources. + y1 : array + y-Positions of the emitting heat sources. + H1 : array + Lengths of the emitting heat sources. + D1 : array + Buried depths of the emitting heat sources. + tilt1 : array + Angles (in radians) from vertical of the emitting heat sources. + orientation1 : array + Directions (in radians) of the tilt the emitting heat sources. + x2 : array + x-Positions of the receiving heat sources. + y2 : array + y-Positions of the receiving heat sources. + H2 : array + Lengths of the receiving heat sources. + D2 : array + Buried depths of the receiving heat sources. + tilt2 : array + Angles (in radians) from vertical of the receiving heat sources. + orientation2 : array + Directions (in radians) of the tilt the receiving heat sources. + i_pair : list + Indices of the emitting segments along a borehole. + j_pair : list + Indices of the receiving segments along a borehole. + k_pair : list + Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions + corresponding to all pairs in (i_pair, j_pair). + """ + # Initialize local variables + borehole1 = self.boreholes[i] + borehole2 = self.boreholes[j] + assert reaSource or imgSource, \ + "At least one of reaSource and imgSource must be True." + if reaSource and imgSource: + # Find segment pairs for the full (real + image) FLS solution + compare_pairs = self._compare_realandimage_pairs_inclined + elif reaSource: + # Find segment pairs for the real FLS solution + compare_pairs = self._compare_real_pairs_inclined + elif imgSource: + # Find segment pairs for the image FLS solution + compare_pairs = self._compare_image_pairs_inclined + # Dive both boreholes into segments + segments1 = borehole1.segments( + self.nBoreSegments[i], segment_ratios=self.segment_ratios[i]) + segments2 = borehole2.segments( + self.nBoreSegments[j], segment_ratios=self.segment_ratios[j]) + # Prepare lists of FLS-inclined arguments + rb1 = [] + x1 = [] + y1 = [] + H1 = [] + D1 = [] + tilt1 = [] + orientation1 = [] + x2 = [] + y2 = [] + H2 = [] + D2 = [] + tilt2 = [] + orientation2 = [] + # All possible pairs (i, j) of indices between segments + i_pair = np.repeat(np.arange(self.nBoreSegments[i], dtype=np.uint), + self.nBoreSegments[j]) + j_pair = np.tile(np.arange(self.nBoreSegments[j], dtype=np.uint), + self.nBoreSegments[i]) + # Empty list of indices for unique pairs + k_pair = np.empty(self.nBoreSegments[i] * self.nBoreSegments[j], + dtype=np.uint) + unique_pairs = [] + nPairs = 0 + + p = 0 + for ii, segment_i in enumerate(segments1): + for jj, segment_j in enumerate(segments2): + pair = (segment_i, segment_j) + # Compare the segment pairs to all known unique pairs + for k, pair_k in enumerate(unique_pairs): + m, n = pair_k[0], pair_k[1] + pair_ref = (segments1[m], segments2[n]) + # Stop if a similar pair is found and assign the index + if compare_pairs(pair, pair_ref): + k_pair[p] = k + break + # If no similar pair is found : add a new pair, increment the + # number of unique pairs, and extract the associated buried + # depths + else: + k_pair[p] = nPairs + rb1.append(segment_i.r_b) + x1.append(segment_i.x) + y1.append(segment_i.y) + H1.append(segment_i.H) + D1.append(segment_i.D) + tilt1.append(segment_i.tilt) + orientation1.append(segment_i.orientation) + x2.append(segment_j.x) + y2.append(segment_j.y) + H2.append(segment_j.H) + D2.append(segment_j.D) + tilt2.append(segment_j.tilt) + orientation2.append(segment_j.orientation) + unique_pairs.append((ii, jj)) + nPairs += 1 + p += 1 + return np.array(rb1), np.array(x1), np.array(y1), np.array(H1), \ + np.array(D1), np.array(tilt1), np.array(orientation1), \ + np.array(x2), np.array(y2), np.array(H2), np.array(D2), \ + np.array(tilt2), np.array(orientation2), i_pair, j_pair, k_pair + + def _map_segment_pairs_vertical( + self, i_pair, j_pair, k_pair, borehole_to_borehole, + borehole_to_borehole_indices): + """ + Return the maping of the unique segment-to-segment thermal response + factors (h) to the complete h_ij array of the borefield, such that: + + h_ij[j_segment, i_segment, :nt] = h[:nt, l_segment, k_segment].T, + + where h is the array of unique segment-to-segment thermal response + factors for a given unique pair of boreholes at all unique distances. + + Parameters + ---------- + i_pair : list + Indices of the emitting segments. + j_pair : list + Indices of the receiving segments. + k_pair : list + Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions + corresponding to all pairs in (i_pair, j_pair). + borehole_to_borehole : list + Tuples of borehole indexes. + borehole_to_borehole_indices : list + Indexes of distances. + + Returns + ------- + i_segment : list + Indices of the emitting segments in the bore field. + j_segment : list + Indices of the receiving segments in the bore field. + k_segment : list + Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions + corresponding to all pairs in (i_pair, j_pair) in the bore field. + l_segment : list + Indices of unique distances for all pairs in (i_pair, j_pair) + in the bore field. + + """ + i_segment = np.concatenate( + [i_pair + self._i0Segments[i] for (i, j) in borehole_to_borehole]) + j_segment = np.concatenate( + [j_pair + self._i0Segments[j] for (i, j) in borehole_to_borehole]) + k_segment = np.tile(k_pair, len(borehole_to_borehole)) + l_segment = np.concatenate( + [np.repeat(i, len(k_pair)) for i in borehole_to_borehole_indices]) + return i_segment, j_segment, k_segment, l_segment + + def _map_segment_pairs_inclined( + self, i_pair, j_pair, k_pair, borehole_to_borehole): + """ + Return the maping of the unique segment-to-segment thermal response + factors (h) to the complete h_ij array of the borefield, such that: + + h_ij[j_segment, i_segment, :nt] = h[:nt, k_segment].T, + + where h is the array of unique segment-to-segment thermal response + factors for a given unique pair of boreholes at all unique distances. + + Parameters + ---------- + i_pair : list + Indices of the emitting segments. + j_pair : list + Indices of the receiving segments. + k_pair : list + Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions + corresponding to all pairs in (i_pair, j_pair). + borehole_to_borehole : list + Tuples of borehole indexes. + + Returns + ------- + i_segment : list + Indices of the emitting segments in the bore field. + j_segment : list + Indices of the receiving segments in the bore field. + k_segment : list + Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions + corresponding to all pairs in (i_pair, j_pair) in the bore field. + + """ + i_segment = np.concatenate( + [i_pair + self._i0Segments[i] for (i, j) in borehole_to_borehole]) + j_segment = np.concatenate( + [j_pair + self._i0Segments[j] for (i, j) in borehole_to_borehole]) + k_segment = np.tile(k_pair, len(borehole_to_borehole)) + return i_segment, j_segment, k_segment + + def _check_solver_specific_inputs(self): + """ + This method ensures that solver specific inputs to the Solver object + are what is expected. + + """ + assert isinstance(self.disTol, (np.floating, float)) and self.disTol > 0., \ + "The distance tolerance 'disTol' should be a positive float." + assert isinstance(self.tol, (np.floating, float)) and self.tol > 0., \ + "The relative tolerance 'tol' should be a positive float." + return diff --git a/setup.cfg b/setup.cfg index aaec5984..10ab4222 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = pygfunction -version = 2.3.1.dev0 +version = 2.3.1 description = A g-function calculator for Python author = Massimo Cimmino author_email = massimo.cimmino@polymtl.ca @@ -22,7 +22,7 @@ classifiers = Topic :: Utilities [options] -packages = pygfunction +packages = find: install_requires = numpy >= 1.26.4 scipy >= 1.13.1 @@ -43,3 +43,9 @@ testing = pytest >= 7.4.4 pytest-cov >= 6.1.1 tox >= 4.25.0 + +[options.packages.find] +exclude = + doc* + examples* + tests* diff --git a/tests/borefield_test.py b/tests/borefield_test.py index d95e5c86..440db758 100644 --- a/tests/borefield_test.py +++ b/tests/borefield_test.py @@ -32,6 +32,35 @@ def test_borefield_init(field, request): H, D, r_b, x, y, tilt=tilt, orientation=orientation) assert borefield == borefield_from_boreholes + +# Test Borefield.__add__ and Borefield.__radd__ +@pytest.mark.parametrize("field, other_field, field_list, other_field_list, field_borehole, other_field_borehole", [ + # Using Borefield objects + ('ten_boreholes_rectangular', 'two_boreholes_inclined', False, False, False, False), + # Using Borefield objects + ('single_borehole', 'two_boreholes_inclined', False, False, True, False), + ('ten_boreholes_rectangular', 'single_borehole_short', False, False, False, True), + # Using Borefield as lists + ('ten_boreholes_rectangular', 'two_boreholes_inclined', False, True, False, False), + ('ten_boreholes_rectangular', 'two_boreholes_inclined', True, False, False, False), + ]) +def test_borefield_add(field, other_field, field_list, other_field_list, field_borehole, other_field_borehole, request): + field = request.getfixturevalue(field) + other_field = request.getfixturevalue(other_field) + reference_field = gt.borefield.Borefield.from_boreholes( + field.to_boreholes() + other_field.to_boreholes() + ) + if field_list: + field = field.to_boreholes() + if other_field_list: + other_field = other_field.to_boreholes() + if field_borehole: + field = field[0] + if other_field_borehole: + other_field = other_field[0] + assert field + other_field == reference_field + + # Test borefield comparison using __eq__ @pytest.mark.parametrize("field, other_field, expected", [ # Fields that are equal @@ -53,6 +82,7 @@ def test_borefield_eq(field, other_field, expected, request): other_field = request.getfixturevalue(other_field) assert (borefield == other_field) == expected + # Test borefield comparison using __ne__ @pytest.mark.parametrize("field, other_field, expected", [ # Fields that are equal diff --git a/tests/boreholes_test.py b/tests/boreholes_test.py index 4a43a771..64a5ae94 100644 --- a/tests/boreholes_test.py +++ b/tests/boreholes_test.py @@ -33,6 +33,24 @@ def test_borehole_init(): ]) +# Test Borehole.__add__ +@pytest.mark.parametrize("borehole, other_borehole, borehole_list, other_borehole_list", [ + ('single_borehole', 'single_borehole_short', False, False), + ('single_borehole', 'single_borehole_short', True, False), + ('single_borehole', 'single_borehole_short', False, True), + ]) +def test_borehole_add(borehole, other_borehole, borehole_list, other_borehole_list, request): + borehole = request.getfixturevalue(borehole)[0] + other_borehole = request.getfixturevalue(other_borehole)[0] + field = gt.borefield.Borefield.from_boreholes( + [borehole, other_borehole]) + if borehole_list: + borehole = [borehole] + if other_borehole_list: + other_borehole = [other_borehole] + assert field == borehole + other_borehole + + # Test Borehole.distance @pytest.mark.parametrize("borehole1, borehole2", [ # Same borehole diff --git a/tests/gfunction_test.py b/tests/gfunction_test.py index 98255c19..5c64d4fc 100644 --- a/tests/gfunction_test.py +++ b/tests/gfunction_test.py @@ -80,6 +80,30 @@ def test_gfunctions_UBWT(field, method, opts, expected, request): boundary_condition='UBWT') assert np.allclose(gFunc.gFunc, expected) + # test from static parameters + + # convert to lists for testing + H = list(borefield.H) + D = list(borefield.D) + r_b = list(borefield.r_b) + x = list(borefield.x) + y = list(borefield.y) + + # get gFunction object with static parameters + gFunc_from_params = gt.gfunction.gFunction.from_static_params( + H=H, + D=D, + r_b=r_b, + x=x, + y=y, + alpha=alpha, + options=options, + method=method, + boundary_condition='UBWT', + ) + gFunc_from_params.evaluate_g_function(time) + assert np.allclose(gFunc_from_params.gFunc, expected) + # Test 'UHTR' g-functions for different bore fields using all solvers, # unequal/uniform segments, and with/without the FLS approximation @@ -151,6 +175,30 @@ def test_gfunctions_UHTR(field, method, opts, expected, request): boundary_condition='UHTR') assert np.allclose(gFunc.gFunc, expected) + # test from static parameters + + # convert to lists for testing + H = list(borefield.H) + D = list(borefield.D) + r_b = list(borefield.r_b) + x = list(borefield.x) + y = list(borefield.y) + + # get gFunction object with static parameters + gFunc_from_params = gt.gfunction.gFunction.from_static_params( + H=H, + D=D, + r_b=r_b, + x=x, + y=y, + alpha=alpha, + options=options, + method=method, + boundary_condition='UHTR', + ) + gFunc_from_params.evaluate_g_function(time) + assert np.allclose(gFunc_from_params.gFunc, expected) + # Test 'MIFT' g-functions for different bore fields using all solvers, # unequal/uniform segments, and with/without the FLS approximation @@ -234,6 +282,56 @@ def test_gfunctions_MIFT( cp_f=fluid.cp, method=method, options=options, boundary_condition='MIFT') assert np.allclose(gFunc.gFunc, expected) + # test from static parameters + + # convert to lists for testing + H = [bore.H for bore in network.b] + D = [bore.D for bore in network.b] + r_b = [bore.r_b for bore in network.b] + x = [bore.x for bore in network.b] + y = [bore.y for bore in network.b] + + # static params + k_s = 2.0 + k_g = 1.0 + k_p = 0.4 + epsilon = 1e-6 + fluid_name = 'MPG' + fluid_pct = 20. + + # Extract the pipe options from the fixture + pipe = network.p[0] + pos = pipe.pos + r_in = pipe.r_in + r_out = pipe.r_out + + # get gFunction object with static parameters + gFunc_from_params = gt.gfunction.gFunction.from_static_params( + H=H, + D=D, + r_b=r_b, + x=x, + y=y, + alpha=alpha, + options=options, + method=method, + boundary_condition='MIFT', + m_flow_network=m_flow_network, + fluid_str=fluid_name, + fluid_concentration_pct=fluid_pct, + pipe_type_str='single_utube', + pos=pos, + r_in=r_in, + r_out=r_out, + k_s=k_s, + k_g=k_g, + k_p=k_p, + epsilon=epsilon, + bore_connectivity=network.c, + ) + gFunc_from_params.evaluate_g_function(time) + assert np.allclose(gFunc_from_params.gFunc, expected) + # Test 'MIFT' g-functions for different bore fields using all solvers, unequal # segments, and with the FLS approximation for variable mass flow rate @@ -274,6 +372,7 @@ def test_gfunctions_MIFT_variable_mass_flow_rate( network, alpha, time=time, m_flow_network=m_flow_network, cp_f=fluid.cp, method=method, options=options, boundary_condition='MIFT') assert np.allclose(gFunc.gFunc, expected) + # variable mass flow not currently supported from gFunction.from_static_params # ============================================================================= @@ -299,7 +398,7 @@ def test_gfunctions_MIFT_variable_mass_flow_rate( # 'detailed' solver - uniform segments, FLS approximation ('detailed', 'uniform_segments_approx', np.array([5.67916493, 6.7395222 , 7.16339216])), ]) -def test_gfunctions_UBWT(two_boreholes_inclined, method, opts, expected, request): +def test_gfunctions_UBWT_inclined(two_boreholes_inclined, method, opts, expected, request): # Extract the bore field from the fixture borefield = two_boreholes_inclined # Extract the g-function options from the fixture @@ -317,6 +416,34 @@ def test_gfunctions_UBWT(two_boreholes_inclined, method, opts, expected, request boundary_condition='UBWT') assert np.allclose(gFunc.gFunc, expected) + # test from static parameters + + # convert to lists for testing + H = list(borefield.H) + D = list(borefield.D) + r_b = list(borefield.r_b) + x = list(borefield.x) + y = list(borefield.y) + tilt = list(borefield.tilt) + orientation = list(borefield.orientation) + + # get gFunction object with static parameters + gFunc_from_params = gt.gfunction.gFunction.from_static_params( + H=H, + D=D, + r_b=r_b, + x=x, + y=y, + tilt=tilt, + orientation=orientation, + alpha=alpha, + options=options, + method=method, + boundary_condition='UBWT', + ) + gFunc_from_params.evaluate_g_function(time) + assert np.allclose(gFunc_from_params.gFunc, expected) + # ============================================================================= # Test gFunction linearization @@ -339,3 +466,141 @@ def test_gfunctions_UBWT_linearization(field, method, opts, expected, request): borefield, alpha, time=time, method=method, options=options, boundary_condition='UBWT') assert np.allclose(gFunc.gFunc, expected) + + # test from static parameters + + # convert to lists for testing + H = list(borefield.H) + D = list(borefield.D) + r_b = list(borefield.r_b) + x = list(borefield.x) + y = list(borefield.y) + + # get gFunction object with static parameters + gFunc_from_params = gt.gfunction.gFunction.from_static_params( + H=H, + D=D, + r_b=r_b, + x=x, + y=y, + alpha=alpha, + options=options, + method=method, + boundary_condition='UBWT', + ) + gFunc_from_params.evaluate_g_function(time) + assert np.allclose(gFunc_from_params.gFunc, expected) + + +@pytest.mark.parametrize("field, boundary_condition, method, opts, pipe_type, m_flow_network, expected", [ + # 'equivalent' solver - unequal segments - UBWT - single u-tube + ('single_borehole', 'UBWT', 'equivalent', 'unequal_segments', 'single_Utube', 0.05, np.array([5.59717446, 6.36257605, 6.60517223])), + ('single_borehole_short', 'UBWT', 'equivalent', 'unequal_segments', 'single_Utube', 0.05, np.array([4.15784411, 4.98477603, 5.27975732])), + ('ten_boreholes_rectangular', 'UBWT', 'equivalent', 'unequal_segments', 'single_Utube', 0.25, np.array([10.89935004, 17.09864925, 19.0795435])), + # 'equivalent' solver - unequal segments - UHTR - single u-tube + ('single_borehole', 'UHTR', 'equivalent', 'unequal_segments', 'single_Utube', 0.05, np.array([5.61855789, 6.41336758, 6.66933682])), + ('single_borehole_short', 'UHTR', 'equivalent', 'unequal_segments', 'single_Utube', 0.05, np.array([4.18276733, 5.03671562, 5.34369772])), + ('ten_boreholes_rectangular', 'UHTR', 'equivalent', 'unequal_segments', 'single_Utube', 0.25, np.array([11.27831804, 18.48075762, 21.00669237])), + # 'equivalent' solver - unequal segments - MIFT - single u-tube + ('single_borehole', 'MIFT', 'equivalent', 'unequal_segments', 'single_Utube', 0.05, np.array([5.76597302, 6.51058473, 6.73746895])), + ('single_borehole_short', 'MIFT', 'equivalent', 'unequal_segments', 'single_Utube', 0.05, np.array([4.17105954, 5.00930075, 5.30832133])), + ('ten_boreholes_rectangular', 'MIFT', 'equivalent', 'unequal_segments', 'single_Utube', 0.25, np.array([12.66229998, 18.57852681, 20.33535907])), + # 'equivalent' solver - unequal segments - MIFT - double u-tube parallel + ('single_borehole', 'MIFT', 'equivalent', 'unequal_segments', 'double_Utube_parallel', 0.05, np.array([6.47497545, 7.18728277, 7.39167598])), + ('single_borehole_short', 'MIFT', 'equivalent', 'unequal_segments', 'double_Utube_parallel', 0.05, np.array([4.17080765, 5.00341368, 5.2989709])), + ('ten_boreholes_rectangular', 'MIFT', 'equivalent', 'unequal_segments', 'double_Utube_parallel', 0.25, np.array([15.96448954, 21.43320976, 22.90761598])), + # 'equivalent' solver - unequal segments - MIFT - double u-tube series + ('single_borehole', 'MIFT', 'equivalent', 'unequal_segments', 'double_Utube_series', 0.05, np.array([5.69118368, 6.44386342, 6.67721347])), + ('single_borehole_short', 'MIFT','equivalent', 'unequal_segments', 'double_Utube_series', 0.05, np.array([4.16750616, 5.00249502, 5.30038701])), + ('ten_boreholes_rectangular', 'MIFT', 'equivalent', 'unequal_segments', 'double_Utube_series', 0.25, np.array([11.94256058, 17.97858109, 19.83460231])), + # 'equivalent' solver - unequal segments - MIFT - double u-tube series asymmetrical + ('single_borehole', 'MIFT', 'equivalent', 'unequal_segments', 'double_Utube_series_asymmetrical', 0.05, np.array([5.69174709, 6.4441862 , 6.67709693])), + ('single_borehole_short', 'MIFT', 'equivalent', 'unequal_segments', 'double_Utube_series_asymmetrical', 0.05, np.array([4.16851817, 5.00453267, 5.30282913])), + ('ten_boreholes_rectangular', 'MIFT', 'equivalent', 'unequal_segments', 'double_Utube_series_asymmetrical', 0.25, np.array([11.96927941, 18.00481705, 19.856554])), + # 'equivalent' solver - unequal segments - MIFT - double u-tube series asymmetrical + ('single_borehole', 'MIFT', 'equivalent', 'unequal_segments', 'double_Utube_series_asymmetrical', 0.05, np.array([5.69174709, 6.4441862, 6.67709693])), + ('single_borehole_short', 'MIFT', 'equivalent', 'unequal_segments', 'double_Utube_series_asymmetrical', 0.05, np.array([4.16851817, 5.00453267, 5.30282913])), + ('ten_boreholes_rectangular', 'MIFT', 'equivalent', 'unequal_segments', 'double_Utube_series_asymmetrical', 0.25, np.array([11.96927941, 18.00481705, 19.856554])), + # 'equivalent' solver - unequal segments - MIFT - coaxial annular inlet + ('single_borehole', 'MIFT', 'equivalent', 'unequal_segments', 'coaxial_annular_in', 0.05, np.array([6.10236427, 6.77069069, 6.95941276])), + ('single_borehole_short', 'MIFT', 'equivalent', 'unequal_segments', 'coaxial_annular_in', 0.05, np.array([4.06874781, 4.89701125, 5.19157017])), + ('ten_boreholes_rectangular', 'MIFT', 'equivalent', 'unequal_segments', 'coaxial_annular_in', 0.25, np.array([16.03433989, 21.18241954, 22.49479982])), + # 'equivalent' solver - unequal segments - MIFT - coaxial annular outlet + ('single_borehole', 'MIFT', 'equivalent', 'unequal_segments', 'coaxial_annular_out', 0.05, np.array([6.10236427, 6.77069069, 6.95941276])), + ('single_borehole_short', 'MIFT', 'equivalent', 'unequal_segments', 'coaxial_annular_out', 0.05, np.array([4.06874781, 4.89701125, 5.19157017])), + ('ten_boreholes_rectangular', 'MIFT', 'equivalent', 'unequal_segments', 'coaxial_annular_out', 0.25, np.array([16.03433989, 21.18241954, 22.49510883])), + ]) +def test_gfunctions_from_static_params(field, boundary_condition, method, opts, pipe_type, m_flow_network, expected, + request): + # Extract the bore field from the fixture for convenience + borefield = request.getfixturevalue(field) + + # convert to lists for testing + H = list(borefield.H) + D = list(borefield.D) + r_b = list(borefield.r_b) + x = list(borefield.x) + y = list(borefield.y) + + # Extract the g-function options from the fixture + options = request.getfixturevalue(opts) + + # Extract the pipe options from the fixture, if needed + pipe = request.getfixturevalue(pipe_type) + pos = pipe.pos + r_in = pipe.r_in + r_out = pipe.r_out + + # replace pipe_type from fixture + if pipe_type in ['single_Utube', 'double_Utube_parallel', 'double_Utube_series', + 'double_Utube_series_asymmetrical']: + k_p = 0.4 + elif pipe_type in ['coaxial_annular_in', 'coaxial_annular_out']: + k_p = (0.4, 0.4) + else: + raise ValueError(f"test pipe_type not recognized: '{pipe_type}'") + + pipe_type = pipe_type.removesuffix('_asymmetrical') + + # Static params + k_s = 2.0 + k_g = 1.0 + epsilon = 1e-6 + fluid_name = 'MPG' + fluid_pct = 20. + + # Mean borehole length [m] + H_mean = np.mean(H) + alpha = 1e-6 # Ground thermal diffusivity [m2/s] + + gfunc = gt.gfunction.gFunction.from_static_params( + H=H, + D=D, + r_b=r_b, + x=x, + y=y, + alpha=alpha, + options=options, + method=method, + boundary_condition=boundary_condition, + k_p=k_p, + k_s=k_s, + k_g=k_g, + epsilon=epsilon, + fluid_str=fluid_name, + fluid_concentration_pct=fluid_pct, + pos=pos, + r_in=r_in, + r_out=r_out, + pipe_type_str=pipe_type, + m_flow_network=m_flow_network, + ) + + # Bore field characteristic time [s] + ts = H_mean ** 2 / (9 * alpha) + # Times for the g-function [s] + time = np.array([0.1, 1., 10.]) * ts + + # evaluate g-function + gFunc = gfunc.evaluate_g_function(time=time) + assert np.allclose(gFunc, expected)