From 841f7ce657ec3fca11ca8230652d9d3bdac6935b Mon Sep 17 00:00:00 2001 From: MassimoCimmino Date: Sun, 27 Apr 2025 11:37:24 -0400 Subject: [PATCH 1/7] Refactor heat_transfer --- pygfunction/borefield.py | 55 +- pygfunction/boreholes.py | 51 +- pygfunction/heat_transfer.py | 1516 ----------------- pygfunction/heat_transfer/__init__.py | 3 + .../heat_transfer/finite_line_source.py | 244 +++ .../finite_line_source_inclined.py | 696 ++++++++ .../finite_line_source_vertical.py | 637 +++++++ 7 files changed, 1661 insertions(+), 1541 deletions(-) delete mode 100644 pygfunction/heat_transfer.py create mode 100644 pygfunction/heat_transfer/__init__.py create mode 100644 pygfunction/heat_transfer/finite_line_source.py create mode 100644 pygfunction/heat_transfer/finite_line_source_inclined.py create mode 100644 pygfunction/heat_transfer/finite_line_source_vertical.py diff --git a/pygfunction/borefield.py b/pygfunction/borefield.py index 0a6b7dbc..dc5bf5a8 100644 --- a/pygfunction/borefield.py +++ b/pygfunction/borefield.py @@ -1,9 +1,10 @@ # -*- coding: utf-8 -*- +from collections.abc import Callable from typing import Union, List, Dict, Tuple +from typing_extensions import Self # for compatibility with Python <= 3.10 import numpy as np import numpy.typing as npt -from typing_extensions import Self # for compatibility with Python <= 3.10 from .boreholes import Borehole from .utilities import _initialize_figure, _format_axes, _format_axes_3d @@ -107,6 +108,58 @@ def __ne__( check = not self == other_field return check + @property + def is_tilted(self) -> np.ndarray: + """Return an boolean array. True if the corresponding borehole is inclined.""" + return self._is_tilted + + @property + def is_vertical(self) -> np.ndarray: + """Return an boolean array. True if the corresponding borehole is vertical.""" + return np.logical_not(self._is_tilted) + + def distance( + self, other_field: Union[Borehole, Self, List[Borehole]], + outer: bool = True + ) -> np.ndarray: + """ + Evaluate horizontal distances to boreholes of another borefield. + + Parameters + ---------- + other_field : Borehole or Borefield object + Borefield with which to evaluate distances. + outer : , optional + True if the horizontal distance is to be evaluated for all + boreholes of the borefield onto all boreholes of the other + borefield to return a (nBoreholes_other_field, nBoreholes,) + array. If false, the horizontal distance is evaluated pairwise + between boreholes of the borefield and boreholes of the other + borefield. The numbers of boreholes should be the same + (i.e. nBoreholes == nBoreholes_other_field) and a (nBoreholes,) + array is returned. + Default is True. + + Returns + ------- + dis : array, shape (nBoreholes_other_field, nBoreholes,) or (nBoreholes,) + Horizontal distances between boreholes (in meters). + + """ + if outer: + dis = np.maximum( + np.sqrt( + np.subtract.outer(other_field.x, self.x)**2 + + np.subtract.outer(other_field.y, self.y)**2 + ), + self.r_b) + else: + dis = np.maximum( + np.sqrt( + (other_field.x - self.x)**2 + (other_field.y - self.y)**2), + self.r_b) + return dis + def evaluate_g_function( self, alpha: float, diff --git a/pygfunction/boreholes.py b/pygfunction/boreholes.py index 58c10013..80af83ea 100644 --- a/pygfunction/boreholes.py +++ b/pygfunction/boreholes.py @@ -1,4 +1,5 @@ # -*- coding: utf-8 -*- +from typing_extensions import Self # for compatibility with Python <= 3.10 import warnings import numpy as np @@ -54,6 +55,32 @@ def __repr__(self): f' orientation={self.orientation})') return s + @property + def is_tilted(self) -> bool: + """ + Returns true if the borehole is inclined. + + Returns + ------- + bool + True if borehole is inclined. + + """ + return self._is_tilted + + @property + def is_vertical(self) -> bool: + """ + Returns true if the borehole is vertical. + + Returns + ------- + bool + True if borehole is vertical. + + """ + return not self._is_tilted + def distance(self, target): """ Evaluate the distance between the current borehole and a target @@ -86,30 +113,6 @@ def distance(self, target): np.sqrt((self.x - target.x)**2 + (self.y - target.y)**2)) return dis - def is_tilted(self): - """ - Returns true if the borehole is inclined. - - Returns - ------- - bool - True if borehole is inclined. - - """ - return self._is_tilted - - def is_vertical(self): - """ - Returns true if the borehole is vertical. - - Returns - ------- - bool - True if borehole is vertical. - - """ - return not self._is_tilted - def position(self): """ Returns the position of the borehole. diff --git a/pygfunction/heat_transfer.py b/pygfunction/heat_transfer.py deleted file mode 100644 index 55f3cf9c..00000000 --- a/pygfunction/heat_transfer.py +++ /dev/null @@ -1,1516 +0,0 @@ -# -*- coding: utf-8 -*- -import numpy as np -from scipy.integrate import quad, quad_vec -from scipy.special import erfc, erf, roots_legendre - -from .boreholes import Borehole -from .utilities import erfint, exp1, _erf_coeffs - - -def finite_line_source( - time, alpha, borehole1, borehole2, reaSource=True, imgSource=True, - approximation=False, M=11, N=10): - """ - Evaluate the Finite Line Source (FLS) solution. - - This function uses a numerical quadrature to evaluate the one-integral form - of the FLS solution. For vertical boreholes, the FLS solution was proposed - by Claesson and Javed [#FLS-ClaJav2011]_ and extended to boreholes with - different vertical positions by Cimmino and Bernier [#FLS-CimBer2014]_. - The FLS solution is given by: - - .. math:: - h_{1\\rightarrow2}(t) &= \\frac{1}{2H_2} - \\int_{\\frac{1}{\\sqrt{4\\alpha t}}}^{\\infty} - e^{-d_{12}^2s^2}(I_{real}(s)+I_{imag}(s))ds - - - d_{12} &= \\sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2} - - - I_{real}(s) &= erfint((D_2-D_1+H_2)s) - erfint((D_2-D_1)s) - - &+ erfint((D_2-D_1-H_1)s) - erfint((D_2-D_1+H_2-H_1)s) - - I_{imag}(s) &= erfint((D_2+D_1+H_2)s) - erfint((D_2+D_1)s) - - &+ erfint((D_2+D_1+H_1)s) - erfint((D_2+D_1+H_2+H_1)s) - - - erfint(X) &= \\int_{0}^{X} erf(x) dx - - &= Xerf(X) - \\frac{1}{\\sqrt{\\pi}}(1-e^{-X^2}) - - For inclined boreholes, the FLS solution was proposed by Lazzarotto - [#FLS-Lazzar2016]_ and Lazzarotto and Björk [#FLS-LazBjo2016]_. - The FLS solution is given by: - - .. math:: - h_{1\\rightarrow2}(t) &= \\frac{H_1}{2H_2} - \\int_{\\frac{1}{\\sqrt{4\\alpha t}}}^{\\infty} - \\frac{1}{s} - \\int_{0}^{1} (I_{real}(u, s)+I_{imag}(u, s)) du ds - - - I_{real}(u, s) &= - e^{-((x_1 - x_2)^2 + (y_1 - y_2)^2 + (D_1 - D_2)^2) s^2} - - &\\cdot (erf((u H_1 k_{0,real} + k_{2,real}) s) - - erf((u H_1 k_{0,real} + k_{2,real} - H_2) s)) - - &\\cdot e^{(u^2 H_1^2 (k_{0,real}^2 - 1) - + 2 u H_1 (k_{0,real} k_{2,real} - k_{1,real}) + k_{2,real}^2) s^2} - du ds - - - I_{imag}(u, s) &= - -e^{-((x_1 - x_2)^2 + (y_1 - y_2)^2 + (D_1 + D_2)^2) s^2} - - &\\cdot (erf((u H_1 k_{0,imag} + k_{2,imag}) s) - - erf((u H_1 k_{0,imag} + k_{2,imag} - H_2) s)) - - &\\cdot e^{(u^2 H_1^2 (k_{0,imag}^2 - 1) - + 2 u H_1 (k_{0,imag} k_{2,imag} - k_1) + k_{2,imag}^2) s^2} - du ds - - - k_{0,real} &= - sin(\\beta_1) sin(\\beta_2) cos(\\theta_1 - \\theta_2) - + cos(\\beta_1) cos(\\beta_2) - - - k_{0,imag} &= - sin(\\beta_1) sin(\\beta_2) cos(\\theta_1 - \\theta_2) - - cos(\\beta_1) cos(\\beta_2) - - - k_{1,real} &= sin(\\beta_1) - (cos(\\theta_1) (x_1 - x_2) + sin(\\theta_1) (y_1 - y_2)) - + cos(\\beta_1) (D_1 - D_2) - - - k_{1,imag} &= sin(\\beta_1) - (cos(\\theta_1) (x_1 - x_2) + sin(\\theta_1) (y_1 - y_2)) - + cos(\\beta_1) (D_1 + D_2) - - - k_{2,real} &= sin(\\beta_2) - (cos(\\theta_2) (x_1 - x_2) + sin(\\theta_2) (y_1 - y_2)) - + cos(\\beta_2) (D_1 - D_2) - - - k_{2,imag} &= sin(\\beta_2) - (cos(\\theta_2) (x_1 - x_2) + sin(\\theta_2) (y_1 - y_2)) - - cos(\\beta_2) (D_1 + D_2) - - where :math:`\\beta_1` and :math:`\\beta_2` are the tilt angle of the - boreholes (relative to vertical), and :math:`\\theta_1` and - :math:`\\theta_2` are the orientation of the boreholes (relative to the - x-axis). - - .. Note:: - The reciprocal thermal response factor - :math:`h_{2\\rightarrow1}(t)` can be conveniently calculated by: - - .. math:: - h_{2\\rightarrow1}(t) = \\frac{H_2}{H_1} - h_{1\\rightarrow2}(t) - - Parameters - ---------- - time : float or array, shape (K) - Value of time (in seconds) for which the FLS solution is evaluated. - alpha : float - Soil thermal diffusivity (in m2/s). - borehole1 : Borehole object or list of Borehole objects, length (N) - Borehole object of the borehole extracting heat. - borehole2 : Borehole object or list of Borehole objects, length (M) - Borehole object for which the FLS is evaluated. - reaSource : bool - True if the real part of the FLS solution is to be included. - Default is True. - imgSource : bool, optional - True if the image part of the FLS solution is to be included. - Default is True. - approximation : bool, optional - Set to true to use the approximation of the FLS solution of Cimmino - (2021) [#FLS-Cimmin2021]_. This approximation does not require - the numerical evaluation of any integral. - Default is False. - M : int, optional - Number of Gauss-Legendre sample points for the quadrature over - :math:`u`. This is only used for inclined boreholes. - Default is 11. - N : int, optional - Number of terms in the approximation of the FLS solution. This - parameter is unused if `approximation` is set to False. - Default is 10. Maximum is 25. - - Returns - ------- - h : float or array, shape (M, N, K), (M, N) or (K) - Value of the FLS solution. The average (over the length) temperature - drop on the wall of borehole2 due to heat extracted from borehole1 is: - - .. math:: \\Delta T_{b,2} = T_g - \\frac{Q_1}{2\\pi k_s H_2} h - - Notes - ----- - The function returns a float if time is a float and borehole1 and borehole2 - are Borehole objects. If time is a float and any of borehole1 and borehole2 - are lists, the function returns an array, shape (M, N), If time is an array - and borehole1 and borehole2 are Borehole objects, the function returns an - array, shape (K).If time is an array and any of borehole1 and borehole2 are - are lists, the function returns an array, shape (M, N, K). - - Examples - -------- - >>> b1 = gt.boreholes.Borehole(H=150., D=4., r_b=0.075, x=0., y=0.) - >>> b2 = gt.boreholes.Borehole(H=150., D=4., r_b=0.075, x=5., y=0.) - >>> h = gt.heat_transfer.finite_line_source(4*168*3600., 1.0e-6, b1, b2) - h = 0.0110473635393 - >>> h = gt.heat_transfer.finite_line_source( - 4*168*3600., 1.0e-6, b1, b2, approximation=True, N=10) - h = 0.0110474667731 - >>> b3 = gt.boreholes.Borehole( - H=150., D=4., r_b=0.075, x=5., y=0., tilt=3.1415/15, orientation=0.) - >>> h = gt.heat_transfer.finite_line_source( - 4*168*3600., 1.0e-6, b1, b3, M=21) - h = 0.0002017450051 - - References - ---------- - .. [#FLS-ClaJav2011] Claesson, J., & Javed, S. (2011). An analytical - method to calculate borehole fluid temperatures for time-scales from - minutes to decades. ASHRAE Transactions, 117(2), 279-288. - .. [#FLS-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. - .. [#FLS-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. - .. [#FLS-Lazzar2016] Lazzarotto, A. (2016). A methodology for the - calculation of response functions for geothermal fields with - arbitrarily oriented boreholes – Part 1, Renewable Energy, 86, - 1380-1393. - .. [#FLS-LazBjo2016] Lazzarotto, A., & Björk, F. (2016). A methodology for - the calculation of response functions for geothermal fields with - arbitrarily oriented boreholes – Part 2, Renewable Energy, 86, - 1353-1361. - - """ - if isinstance(borehole1, Borehole) and isinstance(borehole2, Borehole): - # Unpack parameters - H1, D1 = borehole1.H, borehole1.D - H2, D2 = borehole2.H, borehole2.D - if borehole1.is_vertical() and borehole2.is_vertical(): - # Boreholes are vertical - dis = borehole1.distance(borehole2) - if time is np.inf: - # Steady-state solution - h = _finite_line_source_steady_state( - dis, H1, D1, H2, D2, reaSource, imgSource) - elif approximation: - # Approximation - h = finite_line_source_approximation( - time, alpha, dis, H1, D1, H2, D2, - reaSource=reaSource, imgSource=imgSource, N=N) - else: - # Integrand of the finite line source solution - f = _finite_line_source_integrand( - dis, H1, D1, H2, D2, reaSource, imgSource) - # Evaluate integral - if isinstance(time, (np.floating, float)): - # Lower bound of integration - a = 1.0 / np.sqrt(4.0*alpha*time) - h = 0.5 / H2 * quad(f, a, np.inf)[0] - else: - # Lower bound of integration - a = 1.0 / np.sqrt(4.0*alpha*time) - # Upper bound of integration - b = np.concatenate(([np.inf], a[:-1])) - h = np.cumsum(np.stack( - [0.5 / H2 * quad(f, a_i, b_i)[0] - for t, a_i, b_i in zip(time, a, b)], - axis=-1), axis=-1) - else: - # At least one borehole is tilted - # Unpack parameters - x1, y1 = borehole1.position() - rb1 = borehole1.r_b - tilt1 = borehole1.tilt - orientation1 = borehole1.orientation - x2, y2 = borehole2.position() - tilt2 = borehole2.tilt - orientation2 = borehole2.orientation - if time is np.inf: - # Steady-state solution - h = _finite_line_source_inclined_steady_state( - rb1, x1, y1, H1, D1, tilt1, orientation1, - x2, y2, H2, D2, tilt2, orientation2, - reaSource, imgSource, M=M) - elif approximation: - # Approximation - h = finite_line_source_inclined_approximation( - time, alpha, rb1, x1, y1, H1, D1, tilt1, orientation1, - x2, y2, H2, D2, tilt2, orientation2, - reaSource=reaSource, imgSource=imgSource, M=M, N=N) - else: - # Integrand of the inclined finite line source solution - f = _finite_line_source_inclined_integrand( - rb1, x1, y1, H1, D1, tilt1, orientation1, - x2, y2, H2, D2, tilt2, orientation2, - reaSource, imgSource, M) - - # Evaluate integral - if isinstance(time, (np.floating, float)): - # Lower bound of integration - a = 1.0 / np.sqrt(4.0*alpha*time) - h = 0.5 / H2 * quad(f, a, np.inf)[0] - else: - # Lower bound of integration - a = 1.0 / np.sqrt(4.0*alpha*time) - # Upper bound of integration - b = np.concatenate(([np.inf], a[:-1])) - h = np.cumsum(np.stack( - [0.5 / H2 * quad(f, a_i, b_i)[0] - for t, a_i, b_i in zip(time, a, b)], - axis=-1), axis=-1) - - else: - # Unpack parameters - if isinstance(borehole1, Borehole): borehole1 = [borehole1] - if isinstance(borehole2, Borehole): borehole2 = [borehole2] - x1 = np.array([b.x for b in borehole1]) - y1 = np.array([b.y for b in borehole1]) - x2 = np.array([b.x for b in borehole2]) - y2 = np.array([b.y for b in borehole2]) - r_b = np.array([b.r_b for b in borehole1]) - dis = np.maximum( - np.sqrt(np.add.outer(x2, -x1)**2 + np.add.outer(y2, -y1)**2), - r_b) - D1 = np.array([b.D for b in borehole1]).reshape(1, -1) - H1 = np.array([b.H for b in borehole1]).reshape(1, -1) - D2 = np.array([b.D for b in borehole2]).reshape(-1, 1) - H2 = np.array([b.H for b in borehole2]).reshape(-1, 1) - - if (np.all([b.is_vertical() for b in borehole1]) - and np.all([b.is_vertical() for b in borehole2])): - # All boreholes are vertical - if time is np.inf: - # Steady-state solution - h = _finite_line_source_steady_state( - dis, H1, D1, H2, D2, reaSource, imgSource) - elif approximation: - # Approximation - h = finite_line_source_approximation( - time, alpha, dis, H1, D1, H2, D2, - reaSource=reaSource, imgSource=imgSource, N=N) - else: - # Evaluate integral - h = finite_line_source_vectorized( - time, alpha, dis, H1, D1, H2, D2, - reaSource=reaSource, imgSource=imgSource) - else: - # At least one borehole is tilted - # Unpack parameters - x1 = x1.reshape(1, -1) - y1 = y1.reshape(1, -1) - tilt1 = np.array([b.tilt for b in borehole1]).reshape(1, -1) - orientation1 = np.array([b.orientation for b in borehole1]).reshape(1, -1) - x2 = x2.reshape(-1, 1) - y2 = y2.reshape(-1, 1) - tilt2 = np.array([b.tilt for b in borehole2]).reshape(-1, 1) - orientation2 = np.array([b.orientation for b in borehole2]).reshape(-1, 1) - r_b = r_b.reshape(1, -1) - if time is np.inf: - # Steady-state solution - h = _finite_line_source_inclined_steady_state( - r_b, x1, y1, H1, D1, tilt1, orientation1, - x2, y2, H2, D2, tilt2, orientation2, - reaSource, imgSource, M=M) - elif approximation: - # Approximation - h = finite_line_source_inclined_approximation( - time, alpha, r_b, x1, y1, H1, D1, tilt1, orientation1, - x2, y2, H2, D2, tilt2, orientation2, - reaSource=reaSource, imgSource=imgSource, M=M, N=N) - else: - # Evaluate integral - h = finite_line_source_inclined_vectorized( - time, alpha, - r_b, x1, y1, H1, D1, tilt1, orientation1, - x2, y2, H2, D2, tilt2, orientation2, - reaSource=reaSource, imgSource=imgSource, M=M) - return h - - -def finite_line_source_approximation( - time, alpha, dis, H1, D1, H2, D2, reaSource=True, imgSource=True, - N=10): - """ - Evaluate the Finite Line Source (FLS) solution using the approximation - of Cimmino (2021) [#FLSApprox-Cimmin2021]_. - - Parameters - ---------- - time : float or array, shape (K) - Value of time (in seconds) for which the FLS solution is evaluated. - alpha : float - Soil thermal diffusivity (in m2/s). - dis : float or array - Radial distances to evaluate the FLS solution. - H1 : float or array - Lengths of the emitting heat sources. - D1 : float or array - Buried depths of the emitting heat sources. - H2 : float or array - Lengths of the receiving heat sources. - D2 : float or array - Buried depths of the receiving heat sources. - reaSource : bool, optional - True if the real part of the FLS solution is to be included. - Default is True. - imgSource : bool, optional - True if the image part of the FLS solution is to be included. - Default is True. - N : int, optional - Number of terms in the approximation of the FLS solution. This - parameter is unused if `approximation` is set to False. - Default is 10. Maximum is 25. - - Returns - ------- - h : float - Value of the FLS solution. The average (over the length) temperature - drop on the wall of borehole2 due to heat extracted from borehole1 is: - - .. math:: \\Delta T_{b,2} = T_g - \\frac{Q_1}{2\\pi k_s H_2} h - - - References - ---------- - .. [#FLSApprox-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. - - """ - - dis = np.divide.outer(dis, np.sqrt(4*alpha*time)) - H1 = np.divide.outer(H1, np.sqrt(4*alpha*time)) - D1 = np.divide.outer(D1, np.sqrt(4*alpha*time)) - H2 = np.divide.outer(H2, np.sqrt(4*alpha*time)) - D2 = np.divide.outer(D2, np.sqrt(4*alpha*time)) - if reaSource and imgSource: - # Full (real + image) FLS solution - p = np.array([1, -1, 1, -1, 1, -1, 1, -1]) - q = np.abs( - np.stack([D2 - D1 + H2, - D2 - D1, - D2 - D1 - H1, - D2 - D1 + H2 - H1, - D2 + D1 + H2, - D2 + D1, - D2 + D1 + H1, - D2 + D1 + H2 + H1], - axis=-1)) - elif reaSource: - # Real FLS solution - p = np.array([1, -1, 1, -1]) - q = np.abs( - np.stack([D2 - D1 + H2, - D2 - D1, - D2 - D1 - H1, - D2 - D1 + H2 - H1], - axis=-1)) - elif imgSource: - # Image FLS solution - p = np.array([1, -1, 1, -1]) - q = np.abs( - np.stack([D2 + D1 + H2, - D2 + D1, - D2 + D1 + H1, - D2 + D1 + H2 + H1], - axis=-1)) - else: - # No heat source - p = np.zeros(1) - q = np.zeros(1) - # Coefficients of the approximation of the error function - a, b = _erf_coeffs(N) - - dd = dis**2 - qq = q**2 - G1 = np.inner( - p, - q * np.inner( - a, - 0.5* exp1(np.expand_dims(dd, axis=(-1, -2)) + np.multiply.outer(qq, b)))) - x3 = np.sqrt(np.expand_dims(dd, axis=-1) + qq) - G3 = np.inner(p, np.exp(-x3**2) / np.sqrt(np.pi) - x3 * erfc(x3)) - - h = 0.5 / H2 * (G1 + G3) - return h - - -def finite_line_source_inclined_approximation( - time, alpha, - rb1, x1, y1, H1, D1, tilt1, orientation1, - x2, y2, H2, D2, tilt2, orientation2, - reaSource=True, imgSource=True, M=11, N=10): - """ - Evaluate the inclined Finite Line Source (FLS) solution using the - approximation method of Cimmino (2021) [#IncFLSApprox-Cimmin2021]_. - - Parameters - ---------- - time : float or array, shape (K) - Value of time (in seconds) for which the FLS solution is evaluated. - alpha : float - Soil thermal diffusivity (in m2/s). - rb1 : array - Radii of the emitting heat sources. - x1 : float or array - x-Positions of the emitting heat sources. - y1 : float or array - y-Positions of the emitting heat sources. - H1 : float or array - Lengths of the emitting heat sources. - D1 : float or array - Buried depths of the emitting heat sources. - tilt1 : float or array - Angles (in radians) from vertical of the emitting heat sources. - orientation1 : float or 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 : float or array - Lengths of the receiving heat sources. - D2 : float or array - Buried depths of the receiving heat sources. - tilt2 : float or array - Angles (in radians) from vertical of the receiving heat sources. - orientation2 : float or array - Directions (in radians) of the tilt the receiving heat sources. - reaSource : bool, optional - True if the real part of the FLS solution is to be included. - Default is True. - imgSource : bool, optional - True if the image part of the FLS solution is to be included. - Default is true. - M : int, optional - Number of points for the Gauss-Legendre quadrature rule along the - receiving heat sources. - Default is 21. - N : int, optional - Number of terms in the approximation of the FLS solution. - Default is 10. Maximum is 25. - - Returns - ------- - h : float - Value of the FLS solution. The average (over the length) temperature - drop on the wall of borehole2 due to heat extracted from borehole1 is: - - .. math:: \\Delta T_{b,2} = T_g - \\frac{Q_1}{2\\pi k_s H_2} h - - References - ---------- - .. [#IncFLSApprox-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. - - """ - # Expected output shape of h, excluding time - output_shape = np.broadcast_shapes( - *[np.shape(arg) for arg in ( - rb1, x1, y1, H1, D1, tilt1, orientation1, - x2, y2, H2, D2, tilt2, orientation2)]) - # Number of dimensions of the output, excluding time - output_ndim = len(output_shape) - # Shape of the time variable - time_shape = np.shape(time) - # Number of dimensions of the time variable - time_ndim = len(time_shape) - # Roots for Gauss-Legendre quadrature - x, w = roots_legendre(M) - u = (0.5 * x + 0.5).reshape((-1, 1) + (1,) * output_ndim) - w = w / 2 - # Coefficients of the approximation of the error function - a, b = _erf_coeffs(N) - b = b.reshape((1, -1) + (1,) * output_ndim) - # Sines and cosines of tilt (b: beta) and orientation (t: theta) - sb1 = np.sin(tilt1) - sb2 = np.sin(tilt2) - cb1 = np.cos(tilt1) - cb2 = np.cos(tilt2) - st1 = np.sin(orientation1) - st2 = np.sin(orientation2) - ct1 = np.cos(orientation1) - ct2 = np.cos(orientation2) - ct12 = np.cos(orientation1 - orientation2) - # Horizontal distances - dx = x1 - x2 - dy = y1 - y2 - rr = dx**2 + dy**2 # Squared radial distance - # Length ratios - H_ratio = H1 / H2 - H_ratio = np.reshape(H_ratio, np.shape(H_ratio) + (1,) * time_ndim) - # Approximation - ss = 1. / (4 * alpha * time) - if reaSource and imgSource: - # Full (real + image) FLS solution - # Axial distances - dzRea = D1 - D2 - dzImg = D1 + D2 - # FLS-inclined coefficients - kRea_0 = sb1 * sb2 * ct12 + cb1 * cb2 - kImg_0 = sb1 * sb2 * ct12 - cb1 * cb2 - kRea_1 = sb1 * (ct1 * dx + st1 * dy) + cb1 * dzRea - kImg_1 = sb1 * (ct1 * dx + st1 * dy) + cb1 * dzImg - kRea_2 = sb2 * (ct2 * dx + st2 * dy) + cb2 * dzRea - kImg_2 = sb2 * (ct2 * dx + st2 * dy) - cb2 * dzImg - dRea_1 = u * H1 * kRea_0 + kRea_2 - dImg_1 = u * H1 * kImg_0 + kImg_2 - dRea_2 = u * H1 * kRea_0 + kRea_2 - H2 - dImg_2 = u * H1 * kImg_0 + kImg_2 - H2 - cRea = np.maximum( - rr + dzRea**2 - dRea_1**2 + (kRea_1 + u * H1)**2 - kRea_1**2, - rb1**2) - cImg = np.maximum( - rr + dzImg**2 - dImg_1**2 + (kImg_1 + u * H1)**2 - kImg_1**2, - rb1**2) - # Signs for summation - pRea_1 = np.sign(dRea_1) - pRea_1 = np.reshape(pRea_1, np.shape(pRea_1) + (1,) * time_ndim) - pRea_2 = np.sign(dRea_2) - pRea_2 = np.reshape(pRea_2, np.shape(pRea_2) + (1,) * time_ndim) - pImg_1 = np.sign(dImg_1) - pImg_1 = np.reshape(pImg_1, np.shape(pImg_1) + (1,) * time_ndim) - pImg_2 = np.sign(dImg_2) - pImg_2 = np.reshape(pImg_2, np.shape(pImg_2) + (1,) * time_ndim) - # FLS-inclined approximation - h = 0.25 * H_ratio * np.einsum('i,j,ij...', w, a, - (pRea_1 * exp1(np.multiply.outer(cRea + b * dRea_1**2, ss)) \ - - pRea_2 * exp1(np.multiply.outer(cRea + b * dRea_2**2, ss)) \ - - pImg_1 * exp1(np.multiply.outer(cImg + b * dImg_1**2, ss)) \ - + pImg_2 * exp1(np.multiply.outer(cImg + b * dImg_2**2, ss))) ) - elif reaSource: - # Real FLS solution - # Axial distance - dzRea = D1 - D2 - # FLS-inclined coefficients - kRea_0 = sb1 * sb2 * ct12 + cb1 * cb2 - kRea_1 = sb1 * (ct1 * dx + st1 * dy) + cb1 * dzRea - kRea_2 = sb2 * (ct2 * dx + st2 * dy) + cb2 * dzRea - dRea_1 = u * H1 * kRea_0 + kRea_2 - dRea_2 = u * H1 * kRea_0 + kRea_2 - H2 - cRea = np.maximum( - rr + dzRea**2 - dRea_1**2 + (kRea_1 + u * H1)**2 - kRea_1**2, - rb1**2) - # Signs for summation - pRea_1 = np.sign(dRea_1) - pRea_1 = np.reshape(pRea_1, np.shape(pRea_1) + (1,) * time_ndim) - pRea_2 = np.sign(dRea_2) - pRea_2 = np.reshape(pRea_2, np.shape(pRea_2) + (1,) * time_ndim) - # FLS-inclined approximation - h = 0.25 * H_ratio * np.einsum('i,j,ij...', w, a, - (pRea_1 * exp1(np.multiply.outer(cRea + b * dRea_1**2, ss)) \ - - pRea_2 * exp1(np.multiply.outer(cRea + b * dRea_2**2, ss))) ) - elif imgSource: - # Image FLS solution - # Axial distance - dzImg = D1 + D2 - # FLS-inclined coefficients - kImg_0 = sb1 * sb2 * ct12 - cb1 * cb2 - kImg_1 = sb1 * (ct1 * dx + st1 * dy) + cb1 * dzImg - kImg_2 = sb2 * (ct2 * dx + st2 * dy) - cb2 * dzImg - dImg_1 = u * H1 * kImg_0 + kImg_2 - dImg_2 = u * H1 * kImg_0 + kImg_2 - H2 - cImg = np.maximum( - rr + dzImg**2 - dImg_1**2 + (kImg_1 + u * H1)**2 - kImg_1**2, - rb1**2) - # Signs for summation - pImg_1 = np.sign(dImg_1) - pImg_1 = np.reshape(pImg_1, np.shape(pImg_1) + (1,) * time_ndim) - pImg_2 = np.sign(dImg_2) - pImg_2 = np.reshape(pImg_2, np.shape(pImg_2) + (1,) * time_ndim) - # FLS-inclined approximation - h = 0.25 * H_ratio * np.einsum('i,j,ij...', w, a, - (-pImg_1 * exp1(np.multiply.outer(cImg + b * dImg_1**2, ss)) \ - + pImg_2 * exp1(np.multiply.outer(cImg + b * dImg_2**2, ss))) ) - else: - # No heat source - h = np.zeros(output_shape + np.shape(time)) - return h - - -def finite_line_source_vectorized( - time, alpha, dis, H1, D1, H2, D2, reaSource=True, imgSource=True, - approximation=False, N=10): - """ - Evaluate the Finite Line Source (FLS) solution. - - This function uses a numerical quadrature to evaluate the one-integral form - of the FLS solution, as proposed by Claesson and Javed - [#FLSVec-ClaJav2011]_ and extended to boreholes with different vertical - positions by Cimmino and Bernier [#FLSVec-CimBer2014]_. The FLS solution - is given by: - - .. math:: - h_{1\\rightarrow2}(t) &= \\frac{1}{2H_2} - \\int_{\\frac{1}{\\sqrt{4\\alpha t}}}^{\\infty} - e^{-d_{12}^2s^2}(I_{real}(s)+I_{imag}(s))ds - - - I_{real}(s) &= erfint((D_2-D_1+H_2)s) - erfint((D_2-D_1)s) - - &+ erfint((D_2-D_1-H_1)s) - erfint((D_2-D_1+H_2-H_1)s) - - I_{imag}(s) &= erfint((D_2+D_1+H_2)s) - erfint((D_2+D_1)s) - - &+ erfint((D_2+D_1+H_1)s) - erfint((D_2+D_1+H_2+H_1)s) - - - erfint(X) &= \\int_{0}^{X} erf(x) dx - - &= Xerf(X) - \\frac{1}{\\sqrt{\\pi}}(1-e^{-X^2}) - - .. Note:: - The reciprocal thermal response factor - :math:`h_{2\\rightarrow1}(t)` can be conveniently calculated by: - - .. math:: - h_{2\\rightarrow1}(t) = \\frac{H_2}{H_1} - h_{1\\rightarrow2}(t) - - Parameters - ---------- - time : float or array, shape (K) - Value of time (in seconds) for which the FLS solution is evaluated. - alpha : float - Soil thermal diffusivity (in m2/s). - dis : float or array - Radial distances to evaluate the FLS solution. - H1 : float or array - Lengths of the emitting heat sources. - D1 : float or array - Buried depths of the emitting heat sources. - H2 : float or array - Lengths of the receiving heat sources. - D2 : float or array - Buried depths of the receiving heat sources. - reaSource : bool - True if the real part of the FLS solution is to be included. - Default is True. - imgSource : bool - True if the image part of the FLS solution is to be included. - Default is True. - approximation : bool, optional - Set to true to use the approximation of the FLS solution of Cimmino - (2021) [#FLSVec-Cimmin2021]_. This approximation does not require - the numerical evaluation of any integral. - Default is False. - N : int, optional - Number of terms in the approximation of the FLS solution. This - parameter is unused if `approximation` is set to False. - Default is 10. Maximum is 25. - - Returns - ------- - h : float - Value of the FLS solution. The average (over the length) temperature - drop on the wall of borehole2 due to heat extracted from borehole1 is: - - .. math:: \\Delta T_{b,2} = T_g - \\frac{Q_1}{2\\pi k_s H_2} h - - Notes - ----- - This is a vectorized version of the :func:`finite_line_source` function - using scipy.integrate.quad_vec to speed up calculations. All arrays - (dis, H1, D1, H2, D2) must follow numpy array broadcasting rules. If time - is an array, the integrals for different time values are stacked on the - last axis. - - - References - ---------- - .. [#FLSVec-ClaJav2011] Claesson, J., & Javed, S. (2011). An analytical - method to calculate borehole fluid temperatures for time-scales from - minutes to decades. ASHRAE Transactions, 117(2), 279-288. - .. [#FLSVec-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. - .. [#FLSVec-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. - - """ - if not approximation: - # Integrand of the finite line source solution - f = _finite_line_source_integrand( - dis, H1, D1, H2, D2, reaSource, imgSource) - - # Evaluate integral - if isinstance(time, (np.floating, float)): - # Lower bound of integration - a = 1.0 / np.sqrt(4.0*alpha*time) - h = 0.5 / H2 * quad_vec(f, a, np.inf)[0] - else: - # Lower bound of integration - a = 1.0 / np.sqrt(4.0*alpha*time) - # Upper bound of integration - b = np.concatenate(([np.inf], a[:-1])) - h = np.cumsum(np.stack( - [0.5 / H2 * quad_vec(f, a_i, b_i)[0] - for t, a_i, b_i in zip(time, a, b)], - axis=-1), axis=-1) - else: - h = finite_line_source_approximation( - time, alpha, dis, H1, D1, H2, D2, reaSource=reaSource, - imgSource=imgSource, N=N) - return h - - -def finite_line_source_equivalent_boreholes_vectorized( - time, alpha, dis, wDis, H1, D1, H2, D2, N2, reaSource=True, imgSource=True): - """ - Evaluate the equivalent Finite Line Source (FLS) solution. - - This function uses a numerical quadrature to evaluate the one-integral form - of the FLS solution, as proposed by Prieto and Cimmino - [#eqFLSVec-PriCim2021]_. The equivalent FLS solution is given by: - - .. math:: - h_{1\\rightarrow2}(t) &= \\frac{1}{2 H_2 N_{b,2}} - \\int_{\\frac{1}{\\sqrt{4\\alpha t}}}^{\\infty} - \\sum_{G_1} \\sum_{G_2} - e^{-d_{12}^2s^2}(I_{real}(s)+I_{imag}(s))ds - - - I_{real}(s) &= erfint((D_2-D_1+H_2)s) - erfint((D_2-D_1)s) - - &+ erfint((D_2-D_1-H_1)s) - erfint((D_2-D_1+H_2-H_1)s) - - I_{imag}(s) &= erfint((D_2+D_1+H_2)s) - erfint((D_2+D_1)s) - - &+ erfint((D_2+D_1+H_1)s) - erfint((D_2+D_1+H_2+H_1)s) - - - erfint(X) &= \\int_{0}^{X} erf(x) dx - - &= Xerf(X) - \\frac{1}{\\sqrt{\\pi}}(1-e^{-X^2}) - - .. Note:: - The reciprocal thermal response factor - :math:`h_{2\\rightarrow1}(t)` can be conveniently calculated by: - - .. math:: - h_{2\\rightarrow1}(t) = \\frac{H_2 N_{b,2}}{H_1 N_{b,1}} - h_{1\\rightarrow2}(t) - - Parameters - ---------- - time : float or array, shape (K) - Value of time (in seconds) for which the FLS solution is evaluated. - alpha : float - Soil thermal diffusivity (in m2/s). - dis : array - Unique radial distances to evaluate the FLS solution. - wDis : array - Number of instances of each unique radial distances. - H1 : float or array - Lengths of the emitting heat sources. - D1 : float or array - Buried depths of the emitting heat sources. - H2 : float or array - Lengths of the receiving heat sources. - D2 : float or array - Buried depths of the receiving heat sources. - N2 : float or array, - Number of segments represented by the receiving heat sources. - reaSource : bool - True if the real part of the FLS solution is to be included. - Default is True. - imgSource : bool - True if the image part of the FLS solution is to be included. - Default is True. - - Returns - ------- - h : float - Value of the FLS solution. The average (over the length) temperature - drop on the wall of borehole2 due to heat extracted from borehole1 is: - - .. math:: \\Delta T_{b,2} = T_g - \\frac{Q_1}{2\\pi k_s H_2} h - - Notes - ----- - This is a vectorized version of the :func:`finite_line_source` function - using scipy.integrate.quad_vec to speed up calculations. All arrays - (dis, H1, D1, H2, D2) must follow numpy array broadcasting rules. If time - is an array, the integrals for different time values are stacked on the - last axis. - - - References - ---------- - .. [#eqFLSVec-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. - - """ - # Integrand of the finite line source solution - f = _finite_line_source_equivalent_boreholes_integrand( - dis, wDis, H1, D1, H2, D2, N2, reaSource, imgSource) - - # Evaluate integral - if isinstance(time, (np.floating, float)): - # Lower bound of integration - a = 1.0 / np.sqrt(4.0*alpha*time) - h = 0.5 / (N2*H2) * quad_vec(f, a, np.inf)[0] - else: - # Lower bound of integration - a = 1.0 / np.sqrt(4.0*alpha*time) - # Upper bound of integration - b = np.concatenate(([np.inf], a[:-1])) - h = np.cumsum(np.stack( - [0.5 / (N2*H2) * quad_vec(f, a_i, b_i)[0] - for t, a_i, b_i in zip(time, a, b)], - axis=-1), axis=-1) - return h - - -def finite_line_source_inclined_vectorized( - time, alpha, - rb1, x1, y1, H1, D1, tilt1, orientation1, - x2, y2, H2, D2, tilt2, orientation2, - reaSource=True, imgSource=True, M=11, approximation=False, N=10): - """ - Evaluate the inclined Finite Line Source (FLS) solution. - - This function uses a numerical quadrature to evaluate the inclined FLS - solution, as proposed by Lazzarotto [#incFLSVec-Lazzar2016]_. - The inclined FLS solution is given by: - - .. math:: - h_{1\\rightarrow2}(t) &= \\frac{H_1}{2H_2} - \\int_{\\frac{1}{\\sqrt{4\\alpha t}}}^{\\infty} - \\frac{1}{s} - \\int_{0}^{1} (I_{real}(u, s)+I_{imag}(u, s)) du ds - - - I_{real}(u, s) &= - e^{-((x_1 - x_2)^2 + (y_1 - y_2)^2 + (D_1 - D_2)^2) s^2} - - &\\cdot (erf((u H_1 k_{0,real} + k_{2,real}) s) - - erf((u H_1 k_{0,real} + k_{2,real} - H_2) s)) - - &\\cdot e^{(u^2 H_1^2 (k_{0,real}^2 - 1) - + 2 u H_1 (k_{0,real} k_{2,real} - k_{1,real}) + k_{2,real}^2) s^2} - du ds - - - I_{imag}(u, s) &= - -e^{-((x_1 - x_2)^2 + (y_1 - y_2)^2 + (D_1 + D_2)^2) s^2} - - &\\cdot (erf((u H_1 k_{0,imag} + k_{2,imag}) s) - - erf((u H_1 k_{0,imag} + k_{2,imag} - H_2) s)) - - &\\cdot e^{(u^2 H_1^2 (k_{0,imag}^2 - 1) - + 2 u H_1 (k_{0,imag} k_{2,imag} - k_1) + k_{2,imag}^2) s^2} - du ds - - - k_{0,real} &= - sin(\\beta_1) sin(\\beta_2) cos(\\theta_1 - \\theta_2) - + cos(\\beta_1) cos(\\beta_2) - - - k_{0,imag} &= - sin(\\beta_1) sin(\\beta_2) cos(\\theta_1 - \\theta_2) - - cos(\\beta_1) cos(\\beta_2) - - - k_{1,real} &= sin(\\beta_1) - (cos(\\theta_1) (x_1 - x_2) + sin(\\theta_1) (y_1 - y_2)) - + cos(\\beta_1) (D_1 - D_2) - - - k_{1,imag} &= sin(\\beta_1) - (cos(\\theta_1) (x_1 - x_2) + sin(\\theta_1) (y_1 - y_2)) - + cos(\\beta_1) (D_1 + D_2) - - - k_{2,real} &= sin(\\beta_2) - (cos(\\theta_2) (x_1 - x_2) + sin(\\theta_2) (y_1 - y_2)) - + cos(\\beta_2) (D_1 - D_2) - - - k_{2,imag} &= sin(\\beta_2) - (cos(\\theta_2) (x_1 - x_2) + sin(\\theta_2) (y_1 - y_2)) - - cos(\\beta_2) (D_1 + D_2) - - where :math:`\\beta_1` and :math:`\\beta_2` are the tilt angle of the - boreholes (relative to vertical), and :math:`\\theta_1` and - :math:`\\theta_2` are the orientation of the boreholes (relative to the - x-axis). - - .. Note:: - The reciprocal thermal response factor - :math:`h_{2\\rightarrow1}(t)` can be conveniently calculated by: - - .. math:: - h_{2\\rightarrow1}(t) = \\frac{H_2}{H_1} - h_{1\\rightarrow2}(t) - - Parameters - ---------- - time : float or array, shape (K) - Value of time (in seconds) for which the FLS solution is evaluated. - alpha : float - Soil thermal diffusivity (in m2/s). - rb1 : array - Radii of the emitting heat sources. - x1 : float or array - x-Positions of the emitting heat sources. - y1 : float or array - y-Positions of the emitting heat sources. - H1 : float or array - Lengths of the emitting heat sources. - D1 : float or array - Buried depths of the emitting heat sources. - tilt1 : float or array - Angles (in radians) from vertical of the emitting heat sources. - orientation1 : float or 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 : float or array - Lengths of the receiving heat sources. - D2 : float or array - Buried depths of the receiving heat sources. - tilt2 : float or array - Angles (in radians) from vertical of the receiving heat sources. - orientation2 : float or array - Directions (in radians) of the tilt the receiving heat sources. - reaSource : bool, optional - True if the real part of the FLS solution is to be included. - Default is True. - imgSource : bool, optional - True if the image part of the FLS solution is to be included. - Default is true. - M : int, optional - Number of points for the Gauss-Legendre quadrature rule along the - receiving heat sources. - Default is 21. - approximation : bool, optional - Set to true to use the approximation of the FLS solution of Cimmino - (2021) [#FLSVec-Cimmin2021]_. This approximation does not require - the numerical evaluation of any integral. - Default is False. - N : int, optional - Number of terms in the approximation of the FLS solution. This - parameter is unused if `approximation` is set to False. - Default is 10. Maximum is 25. - - Returns - ------- - f : callable - Integrand of the finite line source solution. Can be vector-valued. - - Notes - ----- - This is a vectorized version of the :func:`finite_line_source` function - using scipy.integrate.quad_vec to speed up calculations. All arrays - (x1, y1, H1, D1, tilt1, orientation1, x2, y2, H2, D2, tilt2, - orientation2) must follow numpy array broadcasting rules. - - References - ---------- - .. [#incFLSVec-Lazzar2016] Lazzarotto, A. (2016). A methodology for the - calculation of response functions for geothermal fields with - arbitrarily oriented boreholes – Part 1, Renewable Energy, 86, - 1380-1393. - - """ - if not approximation: - # Integrand of the inclined finite line source solution - f = _finite_line_source_inclined_integrand( - rb1, x1, y1, H1, D1, tilt1, orientation1, - x2, y2, H2, D2, tilt2, orientation2, - reaSource, imgSource, M) - - # Evaluate integral - if isinstance(time, (np.floating, float)): - # Lower bound of integration - a = 1.0 / np.sqrt(4.0*alpha*time) - h = 0.5 / H2 * quad_vec(f, a, np.inf, epsabs=1e-4, epsrel=1e-6)[0] - else: - # Lower bound of integration - a = 1.0 / np.sqrt(4.0*alpha*time) - # Upper bound of integration - b = np.concatenate(([np.inf], a[:-1])) - h = np.cumsum( - np.stack( - [0.5 / H2 * quad_vec( - f, a_i, b_i, epsabs=1e-4, epsrel=1e-6)[0] - for i, (a_i, b_i) in enumerate(zip(a, b))], - axis=-1), - axis=-1) - else: - h = finite_line_source_inclined_approximation( - time, alpha, rb1, x1, y1, H1, D1, tilt1, orientation1, - x2, y2, H2, D2, tilt2, orientation2, - reaSource=reaSource, imgSource=imgSource, M=M, N=N) - return h - - -def _finite_line_source_integrand(dis, H1, D1, H2, D2, reaSource, imgSource): - """ - Integrand of the finite line source solution. - - Parameters - ---------- - dis : float or array - Radial distances to evaluate the FLS solution. - H1 : float or array - Lengths of the emitting heat sources. - D1 : float or array - Buried depths of the emitting heat sources. - H2 : float or array - Lengths of the receiving heat sources. - D2 : float or array - Buried depths of the receiving heat sources. - reaSource : bool - True if the real part of the FLS solution is to be included. - imgSource : bool - True if the image part of the FLS solution is to be included. - - Returns - ------- - f : callable - Integrand of the finite line source solution. Can be vector-valued. - - Notes - ----- - All arrays (dis, H1, D1, H2, D2) must follow numpy array broadcasting - rules. - - """ - if reaSource and imgSource: - # Full (real + image) FLS solution - p = np.array([1, -1, 1, -1, 1, -1, 1, -1]) - q = np.stack([D2 - D1 + H2, - D2 - D1, - D2 - D1 - H1, - D2 - D1 + H2 - H1, - D2 + D1 + H2, - D2 + D1, - D2 + D1 + H1, - D2 + D1 + H2 + H1], - axis=-1) - f = lambda s: s**-2 * np.exp(-dis**2*s**2) * np.inner(p, erfint(q*s)) - elif reaSource: - # Real FLS solution - p = np.array([1, -1, 1, -1]) - q = np.stack([D2 - D1 + H2, - D2 - D1, - D2 - D1 - H1, - D2 - D1 + H2 - H1], - axis=-1) - f = lambda s: s**-2 * np.exp(-dis**2*s**2) * np.inner(p, erfint(q*s)) - elif imgSource: - # Image FLS solution - p = np.array([1, -1, 1, -1]) - q = np.stack([D2 + D1 + H2, - D2 + D1, - D2 + D1 + H1, - D2 + D1 + H2 + H1], - axis=-1) - f = lambda s: s**-2 * np.exp(-dis**2*s**2) * np.inner(p, erfint(q*s)) - else: - # No heat source - f = lambda s: np.zeros(np.broadcast_shapes( - *[np.shape(arg) for arg in (dis, H1, D1, H2, D2)])) - return f - - -def _finite_line_source_inclined_integrand( - rb1, x1, y1, H1, D1, tilt1, orientation1, x2, y2, H2, D2, tilt2, orientation2, - reaSource, imgSource, M): - """ - Integrand of the inclined Finite Line Source (FLS) solution. - - Parameters - ---------- - rb1 : array - Radii of the emitting heat sources. - x1 : float or array - x-Positions of the emitting heat sources. - y1 : float or array - y-Positions of the emitting heat sources. - H1 : float or array - Lengths of the emitting heat sources. - D1 : float or array - Buried depths of the emitting heat sources. - tilt1 : float or array - Angles (in radians) from vertical of the emitting heat sources. - orientation1 : float or 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 : float or array - Lengths of the receiving heat sources. - D2 : float or array - Buried depths of the receiving heat sources. - tilt2 : float or array - Angles (in radians) from vertical of the receiving heat sources. - orientation2 : float or array - Directions (in radians) of the tilt the receiving heat sources. - reaSource : bool - True if the real part of the FLS solution is to be included. - Default is True. - imgSource : bool - True if the image part of the FLS solution is to be included. - M : int - Number of points for the Gauss-Legendre quadrature rule along the - receiving heat sources. - - Returns - ------- - f : callable - Integrand of the finite line source solution. Can be vector-valued. - - Notes - ----- - All arrays (x1, y1, H1, D1, tilt1, orientation1, x2, y2, H2, D2, tilt2, - orientation2) must follow numpy array broadcasting rules. - - """ - output_shape = np.broadcast_shapes( - *[np.shape(arg) for arg in ( - rb1, x1, y1, H1, D1, tilt1, orientation1, - x2, y2, H2, D2, tilt2, orientation2)]) - # Roots - x, w = roots_legendre(M) - u = (0.5 * x + 0.5).reshape((-1,) + (1,) * len(output_shape)) - w = w / 2 - # Params - sb1 = np.sin(tilt1) - sb2 = np.sin(tilt2) - cb1 = np.cos(tilt1) - cb2 = np.cos(tilt2) - dx = x1 - x2 - dy = y1 - y2 - if reaSource and imgSource: - # Full (real + image) FLS solution - dzRea = D1 - D2 - dzImg = D1 + D2 - rr = np.maximum(dx**2 + dy**2, rb1**2) - kRea_0 = sb1 * sb2 * np.cos(orientation1 - orientation2) + cb1 * cb2 - kImg_0 = sb1 * sb2 * np.cos(orientation1 - orientation2) - cb1 * cb2 - kRea_1 = sb1 * (np.cos(orientation1) * dx + np.sin(orientation1) * dy) + cb1 * dzRea - kImg_1 = sb1 * (np.cos(orientation1) * dx + np.sin(orientation1) * dy) + cb1 * dzImg - kRea_2 = sb2 * (np.cos(orientation2) * dx + np.sin(orientation2) * dy) + cb2 * dzRea - kImg_2 = sb2 * (np.cos(orientation2) * dx + np.sin(orientation2) * dy) - cb2 * dzImg - f = lambda s: \ - ((H1 / s * ( - np.exp(-(rr + dzRea**2) * s**2 + s**2 * (u**2 * H1**2 * (kRea_0**2 - 1) + 2 * u * H1 * (kRea_0 * kRea_2 - kRea_1) + kRea_2**2)) \ - * (erf((u * H1 * kRea_0 + kRea_2) * s) - erf((u * H1 * kRea_0 + kRea_2 - H2) * s)) - - np.exp(-(rr + dzImg**2) * s**2 + s**2 * (u**2 * H1**2 * (kImg_0**2 - 1) + 2 * u * H1 * (kImg_0 * kImg_2 - kImg_1) + kImg_2**2)) \ - * (erf((u * H1 * kImg_0 + kImg_2) * s) - erf((u * H1 * kImg_0 + kImg_2 - H2) * s)))).T @ w).T - elif reaSource: - # Real FLS solution - dzRea = D1 - D2 - rr = np.maximum(dx**2 + dy**2, rb1**2) - kRea_0 = sb1 * sb2 * np.cos(orientation1 - orientation2) + cb1 * cb2 - kRea_1 = sb1 * (np.cos(orientation1) * dx + np.sin(orientation1) * dy) + cb1 * dzRea - kRea_2 = sb2 * (np.cos(orientation2) * dx + np.sin(orientation2) * dy) + cb2 * dzRea - f = lambda s: \ - ((H1 / s * np.exp(-(rr + dzRea**2) * s**2 + s**2 * (u**2 * H1**2 * (kRea_0**2 - 1) + 2 * u * H1 * (kRea_0 * kRea_2 - kRea_1) + kRea_2**2)) \ - * (erf((u * H1 * kRea_0 + kRea_2) * s) - erf((u * H1 * kRea_0 + kRea_2 - H2) * s))).T @ w).T - elif imgSource: - # Image FLS solution - dzImg = D1 + D2 - kImg_0 = sb1 * sb2 * np.cos(orientation1 - orientation2) - cb1 * cb2 - kImg_1 = sb1 * (np.cos(orientation1) * dx + np.sin(orientation1) * dy) + cb1 * dzImg - kImg_2 = sb2 * (np.cos(orientation2) * dx + np.sin(orientation2) * dy) - cb2 * dzImg - f = lambda s: \ - -((H1 / s * np.exp(-(dx**2 + dy**2 + dzImg**2) * s**2 + s**2 * (u**2 * H1**2 * (kImg_0**2 - 1) + 2 * u * H1 * (kImg_0 * kImg_2 - kImg_1) + kImg_2**2)) \ - * (erf((u * H1 * kImg_0 + kImg_2) * s) - erf((u * H1 * kImg_0 + kImg_2 - H2) * s))).T @ w).T - else: - # No heat source - f = lambda s: np.zeros(output_shape) - return f - - -def _finite_line_source_equivalent_boreholes_integrand(dis, wDis, H1, D1, H2, D2, N2, reaSource, imgSource): - """ - Integrand of the finite line source solution. - - Parameters - ---------- - dis : array - Unique radial distances to evaluate the FLS solution. - wDis : array - Number of instances of each unique radial distances. - H1 : float or array - Lengths of the emitting heat sources. - D1 : float or array - Buried depths of the emitting heat sources. - H2 : float or array - Lengths of the receiving heat sources. - D2 : float or array - Buried depths of the receiving heat sources. - N2 : float or array, - Number of segments represented by the receiving heat sources. - reaSource : bool - True if the real part of the FLS solution is to be included. - imgSource : bool - True if the image part of the FLS solution is to be included. - - Returns - ------- - f : callable - Integrand of the finite line source solution. Can be vector-valued. - - Notes - ----- - All arrays (dis, H1, D1, H2, D2) must follow numpy array broadcasting - rules. - - """ - if reaSource and imgSource: - # Full (real + image) FLS solution - p = np.array([1, -1, 1, -1, 1, -1, 1, -1]) - q = np.stack([D2 - D1 + H2, - D2 - D1, - D2 - D1 - H1, - D2 - D1 + H2 - H1, - D2 + D1 + H2, - D2 + D1, - D2 + D1 + H1, - D2 + D1 + H2 + H1], - axis=-1) - f = lambda s: s**-2 * (np.exp(-dis**2*s**2) @ wDis).T * np.inner(p, erfint(q*s)) - elif reaSource: - # Real FLS solution - p = np.array([1, -1, 1, -1]) - q = np.stack([D2 - D1 + H2, - D2 - D1, - D2 - D1 - H1, - D2 - D1 + H2 - H1], - axis=-1) - f = lambda s: s**-2 * (np.exp(-dis**2*s**2) @ wDis).T * np.inner(p, erfint(q*s)) - elif imgSource: - # Image FLS solution - p = np.array([1, -1, 1, -1]) - q = np.stack([D2 + D1 + H2, - D2 + D1, - D2 + D1 + H1, - D2 + D1 + H2 + H1], - axis=-1) - f = lambda s: s**-2 * (np.exp(-dis**2*s**2) @ wDis).T * np.inner(p, erfint(q*s)) - else: - # No heat source - f = lambda s: np.zeros(np.broadcast_shapes( - *[np.shape(arg) for arg in (H1, D1, H2, D2, N2)])) - return f - - -def _finite_line_source_steady_state(dis, H1, D1, H2, D2, reaSource, imgSource): - """ - Steady-state finite line source solution. - - Parameters - ---------- - dis : float or array - Radial distances to evaluate the FLS solution. - H1 : float or array - Lengths of the emitting heat sources. - D1 : float or array - Buried depths of the emitting heat sources. - H2 : float or array - Lengths of the receiving heat sources. - D2 : float or array - Buried depths of the receiving heat sources. - reaSource : bool - True if the real part of the FLS solution is to be included. - imgSource : bool - True if the image part of the FLS solution is to be included. - - Returns - ------- - h : Steady-state finite line source solution. - - Notes - ----- - All arrays (dis, H1, D1, H2, D2) must follow numpy array broadcasting - rules. - - """ - # Steady-state solution - if reaSource and imgSource: - # Full (real + image) FLS solution - p = np.array([1, -1, 1, -1, 1, -1, 1, -1]) - q = np.stack([D2 - D1 + H2, - D2 - D1, - D2 - D1 - H1, - D2 - D1 + H2 - H1, - D2 + D1 + H2, - D2 + D1, - D2 + D1 + H1, - D2 + D1 + H2 + H1], - axis=-1) - dis = np.expand_dims(dis, axis=-1) - qpd = np.sqrt(q**2 + dis**2) - h = 0.5 / H2 * np.inner(p, q * np.log(q + qpd) - qpd) - elif reaSource: - # Real FLS solution - p = np.array([1, -1, 1, -1]) - q = np.stack([D2 - D1 + H2, - D2 - D1, - D2 - D1 - H1, - D2 - D1 + H2 - H1,], - axis=-1) - dis = np.expand_dims(dis, axis=-1) - qpd = np.sqrt(q**2 + dis**2) - h = 0.5 / H2 * np.inner(p, q * np.log(q + qpd) - qpd) - elif imgSource: - # Image FLS solution - p = np.array([1, -1, 1, -1]) - q = np.stack([D2 + D1 + H2, - D2 + D1, - D2 + D1 + H1, - D2 + D1 + H2 + H1], - axis=-1) - dis = np.expand_dims(dis, axis=-1) - qpd = np.sqrt(q**2 + dis**2) - h = 0.5 / H2 * np.inner(p, q * np.log(q + qpd) - qpd) - else: - # No heat source - h = np.zeros(np.broadcast_shapes( - *[np.shape(arg) for arg in (dis, H1, D1, H2, D2)])) - return h - - -def _finite_line_source_inclined_steady_state( - rb1, x1, y1, H1, D1, tilt1, orientation1, x2, y2, H2, D2, tilt2, - orientation2, reaSource, imgSource, M=11): - """ - Steady-state inclined Finite Line Source (FLS) solution. - - Parameters - ---------- - rb1 : array - Radii of the emitting heat sources. - x1 : float or array - x-Positions of the emitting heat sources. - y1 : float or array - y-Positions of the emitting heat sources. - H1 : float or array - Lengths of the emitting heat sources. - D1 : float or array - Buried depths of the emitting heat sources. - tilt1 : float or array - Angles (in radians) from vertical of the emitting heat sources. - orientation1 : float or 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 : float or array - Lengths of the receiving heat sources. - D2 : float or array - Buried depths of the receiving heat sources. - tilt2 : float or array - Angles (in radians) from vertical of the receiving heat sources. - orientation2 : float or array - Directions (in radians) of the tilt the receiving heat sources. - reaSource : bool - True if the real part of the FLS solution is to be included. - Default is True. - imgSource : bool - True if the image part of the FLS solution is to be included. - M : int - Number of points for the Gauss-Legendre quadrature rule along the - receiving heat sources. - Default is 11. - - Returns - ------- - h : Steady-state inclined finite line source solution. - - Notes - ----- - All arrays (dis, H1, D1, H2, D2) must follow numpy array broadcasting - rules. - - """ - output_shape = np.broadcast_shapes( - *[np.shape(arg) for arg in ( - rb1, x1, y1, H1, D1, tilt1, orientation1, - x2, y2, H2, D2, tilt2, orientation2)]) - # Roots - x, w = roots_legendre(M) - u = (0.5 * x + 0.5).reshape((-1,) + (1,) * len(output_shape)) - w = w / 2 - # Params - sb1 = np.sin(tilt1) - sb2 = np.sin(tilt2) - cb1 = np.cos(tilt1) - cb2 = np.cos(tilt2) - dx = x1 - x2 - dy = y1 - y2 - # Steady-state solution - if reaSource and imgSource: - # Full (real + image) FLS solution - dzRea = D1 - D2 - dzImg = D1 + D2 - rr = np.maximum(dx**2 + dy**2, rb1**2) - kRea_0 = sb1 * sb2 * np.cos(orientation1 - orientation2) + cb1 * cb2 - kImg_0 = sb1 * sb2 * np.cos(orientation1 - orientation2) - cb1 * cb2 - kRea_1 = sb1 * (np.cos(orientation1) * dx + np.sin(orientation1) * dy) + cb1 * dzRea - kImg_1 = sb1 * (np.cos(orientation1) * dx + np.sin(orientation1) * dy) + cb1 * dzImg - kRea_2 = sb2 * (np.cos(orientation2) * dx + np.sin(orientation2) * dy) + cb2 * dzRea - kImg_2 = sb2 * (np.cos(orientation2) * dx + np.sin(orientation2) * dy) - cb2 * dzImg - h = 0.5 * H1 / H2 * (( - np.log((2 * H2 * np.sqrt(rr + dzRea**2 + u**2*H1**2 + H2**2 + 2*u*H1*kRea_1 - 2*H2*kRea_2 - 2*u*H1*H2*kRea_0) + 2*H2**2 - 2*H2*kRea_2 - 2*u*H1*H2*kRea_0) / ((2 * H2 * np.sqrt(rr + dzRea**2 + u**2*H1**2 + 2*u*H1*kRea_1) - 2*H2*kRea_2 - 2*u*H1*H2*kRea_0))) - - np.log((2 * H2 * np.sqrt(dx**2 + dy**2 + dzImg**2 + u**2*H1**2 + H2**2 + 2*u*H1*kImg_1 - 2*H2*kImg_2 - 2*u*H1*H2*kImg_0) + 2*H2**2 - 2*H2*kImg_2 - 2*u*H1*H2*kImg_0) / ((2 * H2 * np.sqrt(dx**2 + dy**2 + dzImg**2 + u**2*H1**2 + 2*u*H1*kImg_1) - 2*H2*kImg_2 - 2*u*H1*H2*kImg_0))) - ).T @ w).T - elif reaSource: - # Real FLS solution - dzRea = D1 - D2 - rr = np.maximum(dx**2 + dy**2, rb1**2) - kRea_0 = sb1 * sb2 * np.cos(orientation1 - orientation2) + cb1 * cb2 - kRea_1 = sb1 * (np.cos(orientation1) * dx + np.sin(orientation1) * dy) + cb1 * dzRea - kRea_2 = sb2 * (np.cos(orientation2) * dx + np.sin(orientation2) * dy) + cb2 * dzRea - h = 0.5 * H1 / H2 * (( - np.log((2 * H2 * np.sqrt(rr + dzRea**2 + u**2*H1**2 + H2**2 + 2*u*H1*kRea_1 - 2*H2*kRea_2 - 2*u*H1*H2*kRea_0) + 2*H2**2 - 2*H2*kRea_2 - 2*u*H1*H2*kRea_0) / ((2 * H2 * np.sqrt(rr + dzRea**2 + u**2*H1**2 + 2*u*H1*kRea_1) - 2*H2*kRea_2 - 2*u*H1*H2*kRea_0))) - ).T @ w).T - elif imgSource: - # Image FLS solution - dzImg = D1 + D2 - kImg_0 = sb1 * sb2 * np.cos(orientation1 - orientation2) - cb1 * cb2 - kImg_1 = sb1 * (np.cos(orientation1) * dx + np.sin(orientation1) * dy) + cb1 * dzImg - kImg_2 = sb2 * (np.cos(orientation2) * dx + np.sin(orientation2) * dy) - cb2 * dzImg - h = 0.5 * H1 / H2 * (( - - np.log((2 * H2 * np.sqrt(dx**2 + dy**2 + dzImg**2 + u**2*H1**2 + H2**2 + 2*u*H1*kImg_1 - 2*H2*kImg_2 - 2*u*H1*H2*kImg_0) + 2*H2**2 - 2*H2*kImg_2 - 2*u*H1*H2*kImg_0) / ((2 * H2 * np.sqrt(dx**2 + dy**2 + dzImg**2 + u**2*H1**2 + 2*u*H1*kImg_1) - 2*H2*kImg_2 - 2*u*H1*H2*kImg_0))) - ).T @ w).T - else: - # No heat source - h = np.zeros(output_shape) - return h diff --git a/pygfunction/heat_transfer/__init__.py b/pygfunction/heat_transfer/__init__.py new file mode 100644 index 00000000..f28e1831 --- /dev/null +++ b/pygfunction/heat_transfer/__init__.py @@ -0,0 +1,3 @@ +from .finite_line_source import finite_line_source +from .finite_line_source_vertical import finite_line_source_vertical +from .finite_line_source_inclined import finite_line_source_inclined diff --git a/pygfunction/heat_transfer/finite_line_source.py b/pygfunction/heat_transfer/finite_line_source.py new file mode 100644 index 00000000..ea508af9 --- /dev/null +++ b/pygfunction/heat_transfer/finite_line_source.py @@ -0,0 +1,244 @@ +# -*- coding: utf-8 -*- +from typing import Union, List + +import numpy as np +import numpy.typing as npt + +from .finite_line_source_inclined import finite_line_source_inclined +from .finite_line_source_vertical import finite_line_source_vertical +from ..borefield import Borefield +from ..boreholes import Borehole + + +def finite_line_source( + time: npt.ArrayLike, + alpha: float, + borefield_j: Union[Borehole, Borefield, List[Borehole]], + borefield_i: Union[Borehole, Borefield, List[Borehole]], + outer: bool = True, + reaSource: bool = True, + imgSource: bool = True, + approximation: bool = False, + M: int = 11, + N: int = 10 + ) -> np.ndarray: + """ + Evaluate the Finite Line Source (FLS) solution. + + This function uses a numerical quadrature to evaluate the one-integral form + of the FLS solution. For vertical boreholes, the FLS solution was proposed + by Claesson and Javed [#FLS-ClaJav2011]_ and extended to boreholes with + different vertical positions by Cimmino and Bernier [#FLS-CimBer2014]_. + The FLS solution is given by: + + .. math:: + h_{ij}(t) &= \\frac{1}{2H_i} + \\int_{\\frac{1}{\\sqrt{4\\alpha t}}}^{\\infty} + e^{-d_{ij}^2s^2}(I_{real}(s)+I_{imag}(s))ds + + + d_{ij} &= \\sqrt{(x_i - x_j)^2 + (y_i - y_j)^2} + + + I_{real}(s) &= erfint((D_i-D_j+H_i)s) - erfint((D_i-D_j)s) + + &+ erfint((D_i-D_j-H_j)s) - erfint((D_i-D_j+H_i-H_j)s) + + I_{imag}(s) &= erfint((D_i+D_j+H_i)s) - erfint((D_i+D_j)s) + + &+ erfint((D_i+D_j+H_j)s) - erfint((D_i+D_j+H_i+H_j)s) + + + erfint(X) &= \\int_{0}^{X} erf(x) dx + + &= Xerf(X) - \\frac{1}{\\sqrt{\\pi}}(1-e^{-X^2}) + + For inclined boreholes, the FLS solution was proposed by Lazzarotto + [#FLS-Lazzar2016]_ and Lazzarotto and Björk [#FLS-LazBjo2016]_. + The FLS solution is given by: + + .. math:: + h_{1\\rightarrow2}(t) &= \\frac{H_j}{2H_i} + \\int_{\\frac{1}{\\sqrt{4\\alpha t}}}^{\\infty} + \\frac{1}{s} + \\int_{0}^{1} (I_{real}(u, s)+I_{imag}(u, s)) du ds + + + I_{real}(u, s) &= + e^{-((x_i - x_j)^2 + (y_i - y_j)^2 + (D_i - D_j)^2) s^2} + + &\\cdot (erf((u H_j k_{0,real} + k_{2,real}) s) + - erf((u H_j k_{0,real} + k_{2,real} - H_i) s)) + + &\\cdot e^{(u^2 H_j^2 (k_{0,real}^2 - 1) + + 2 u H_j (k_{0,real} k_{2,real} - k_{1,real}) + k_{2,real}^2) s^2} + du ds + + + I_{imag}(u, s) &= + -e^{-((x_i - x_j)^2 + (y_i - y_j)^2 + (D_i + D_j)^2) s^2} + + &\\cdot (erf((u H_j k_{0,imag} + k_{2,imag}) s) + - erf((u H_j k_{0,imag} + k_{2,imag} - H_i) s)) + + &\\cdot e^{(u^2 H_j^2 (k_{0,imag}^2 - 1) + + 2 u H_j (k_{0,imag} k_{2,imag} - k_1) + k_{2,imag}^2) s^2} + du ds + + + k_{0,real} &= + sin(\\beta_j) sin(\\beta_i) cos(\\theta_j - \\theta_i) + + cos(\\beta_j) cos(\\beta_i) + + + k_{0,imag} &= + sin(\\beta_j) sin(\\beta_i) cos(\\theta_j - \\theta_i) + - cos(\\beta_j) cos(\\beta_i) + + + k_{1,real} &= sin(\\beta_j) + (cos(\\theta_j) (x_j - x_i) + sin(\\theta_j) (y_j - y_i)) + + cos(\\beta_j) (D_j - D_i) + + + k_{1,imag} &= sin(\\beta_j) + (cos(\\theta_j) (x_j - x_i) + sin(\\theta_j) (y_j - y_i)) + + cos(\\beta_j) (D_i + D_j) + + + k_{2,real} &= sin(\\beta_i) + (cos(\\theta_i) (x_j - x_i) + sin(\\theta_i) (y_j - y_i)) + + cos(\\beta_i) (D_j - D_i) + + + k_{2,imag} &= sin(\\beta_i) + (cos(\\theta_i) (x_j - x_i) + sin(\\theta_i) (y_j - y_i)) + - cos(\\beta_i) (D_i + D_j) + + where :math:`\\beta_j` and :math:`\\beta_i` are the tilt angle of the + boreholes (relative to vertical), and :math:`\\theta_j` and + :math:`\\theta_i` are the orientation of the boreholes (relative to the + x-axis). + + .. Note:: + The reciprocal thermal response factor + :math:`h_{ji}(t)` can be conveniently calculated by: + + .. math:: + h_{ji}(t) = \\frac{H_i}{H_j} + h_{ij}(t) + + Parameters + ---------- + time : float or (nTimes,) array + Value of time (in seconds) for which the FLS solution is evaluated. + alpha : float + Soil thermal diffusivity (in m2/s). + borefield_j : Borehole or Borefield object + Borehole or Borefield object of the boreholes extracting heat. + borefield_i : Borehole or Borefield object + Borehole or Borefield object object for which the FLS is evaluated. + reaSource : bool + True if the real part of the FLS solution is to be included. + Default is True. + imgSource : bool, optional + True if the image part of the FLS solution is to be included. + Default is True. + outer : bool, optional + True if the finite line source is to be evaluated for all boreholes_j + onto all boreholes_i to return a (nBoreholes_i, nBoreholes_j, nTime,) + array. If false, the finite line source is evaluated pairwise between + boreholes_j and boreholes_i. The numbers of boreholes should be the + same (i.e. nBoreholes_j == nBoreholes_i) and a (nBoreholes, nTime,) + array is returned. + Default is True. + approximation : bool, optional + Set to true to use the approximation of the FLS solution of Cimmino + (2021) [#FLS-Cimmin2021]_. This approximation does not require + the numerical evaluation of any integral. + Default is False. + M : int, optional + Number of Gauss-Legendre sample points for the quadrature over + :math:`u`. This is only used for inclined boreholes. + Default is 11. + N : int, optional + Number of terms in the approximation of the FLS solution. This + parameter is unused if `approximation` is set to False. + Default is 10. Maximum is 25. + + Returns + ------- + h : float or array, shape (nBoreholes_i, nBoreholes_j, nTime,), (nBoreholes, nTime,) or (nTime,) + Value of the FLS solution. The average (over the length) temperature + drop on the wall of borehole_i due to heat extracted from borehole_j + is: + + .. math:: \\Delta T_{b,i} = T_g - \\frac{Q_j}{2\\pi k_s H_j} h + + Notes + ----- + The function returns a float if time is a float and both of borehole_j and + borehole_i are Borehole objects. If time is a float and any of borehole_j + and borehole_i are Borefield objects, the function returns an array. If + time is an array and both of borehole_j and borehole_i are Borehole + objects, the function returns an 1d array, shape (nTime,). If time is an + array and any of borehole_j and borehole_i are are Borefield objects, the + function returns an array. + + Examples + -------- + >>> b1 = gt.boreholes.Borehole(H=150., D=4., r_b=0.075, x=0., y=0.) + >>> b2 = gt.boreholes.Borehole(H=150., D=4., r_b=0.075, x=5., y=0.) + >>> h = gt.heat_transfer.finite_line_source(4*168*3600., 1.0e-6, b1, b2) + h = 0.0110473635393 + >>> h = gt.heat_transfer.finite_line_source( + 4*168*3600., 1.0e-6, b1, b2, approximation=True, N=10) + h = 0.0110474667731 + >>> b3 = gt.boreholes.Borehole( + H=150., D=4., r_b=0.075, x=5., y=0., tilt=3.1415/15, orientation=0.) + >>> h = gt.heat_transfer.finite_line_source( + 4*168*3600., 1.0e-6, b1, b3, M=21) + h = 0.0002017450051 + + References + ---------- + .. [#FLS-ClaJav2011] Claesson, J., & Javed, S. (2011). An analytical + method to calculate borehole fluid temperatures for time-scales from + minutes to decades. ASHRAE Transactions, 117(2), 279-288. + .. [#FLS-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. + .. [#FLS-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. + .. [#FLS-Lazzar2016] Lazzarotto, A. (2016). A methodology for the + calculation of response functions for geothermal fields with + arbitrarily oriented boreholes – Part 1, Renewable Energy, 86, + 1380-1393. + .. [#FLS-LazBjo2016] Lazzarotto, A., & Björk, F. (2016). A methodology for + the calculation of response functions for geothermal fields with + arbitrarily oriented boreholes – Part 2, Renewable Energy, 86, + 1353-1361. + + """ + # Convert bore fields to Borefield objects + if isinstance(borefield_j, list): + borefield_j = Borefield.from_boreholes(borefield_j) + if isinstance(borefield_i, list): + borefield_i = Borefield.from_boreholes(borefield_i) + + # Select the analytical solution based on borehole inclination + if np.any(borefield_j.is_tilted) or np.any(borefield_i.is_tilted): + # Inclined boreholes + h = finite_line_source_inclined( + time, alpha, borefield_j, borefield_i, reaSource=reaSource, + imgSource=imgSource, approximation=approximation, outer=outer, M=M, + N=N) + else: + # Vertical boreholes + h = finite_line_source_vertical( + time, alpha, borefield_j, borefield_i, + reaSource=reaSource, imgSource=imgSource, + approximation=approximation, outer=outer, N=N) + return h diff --git a/pygfunction/heat_transfer/finite_line_source_inclined.py b/pygfunction/heat_transfer/finite_line_source_inclined.py new file mode 100644 index 00000000..41085943 --- /dev/null +++ b/pygfunction/heat_transfer/finite_line_source_inclined.py @@ -0,0 +1,696 @@ +# -*- coding: utf-8 -*- +from collections.abc import Callable +from typing import Union, List, Tuple + +import numpy as np +import numpy.typing as npt +from scipy.integrate import quad_vec +from scipy.special import erf, roots_legendre + +from ..borefield import Borefield +from ..boreholes import Borehole +from ..utilities import exp1, _erf_coeffs + + +def finite_line_source_inclined( + time: npt.ArrayLike, + alpha: float, + borefield_j: Union[Borehole, Borefield, List[Borehole]], + borefield_i: Union[Borehole, Borefield, List[Borehole]], + outer: bool = True, + reaSource: bool = True, + imgSource: bool = True, + approximation: bool = False, + M: int = 11, + N: int = 10 + ) -> np.ndarray: + """ + Evaluate the Finite Line Source (FLS) solution for inclined boreholes. + + This function uses a numerical quadrature to evaluate the one-integral form + of the FLS solution. For inclined boreholes, the FLS solution was proposed + by Lazzarotto [#FLSi-Lazzar2016]_ and Lazzarotto and Björk + [#FLSi-LazBjo2016]_. The FLS solution is given by: + + .. math:: + h_{1\\rightarrow2}(t) &= \\frac{H_j}{2H_i} + \\int_{\\frac{1}{\\sqrt{4\\alpha t}}}^{\\infty} + \\frac{1}{s} + \\int_{0}^{1} (I_{real}(u, s)+I_{imag}(u, s)) du ds + + + I_{real}(u, s) &= + e^{-((x_i - x_j)^2 + (y_i - y_j)^2 + (D_i - D_j)^2) s^2} + + &\\cdot (erf((u H_j k_{0,real} + k_{2,real}) s) + - erf((u H_j k_{0,real} + k_{2,real} - H_i) s)) + + &\\cdot e^{(u^2 H_j^2 (k_{0,real}^2 - 1) + + 2 u H_j (k_{0,real} k_{2,real} - k_{1,real}) + k_{2,real}^2) s^2} + du ds + + + I_{imag}(u, s) &= + -e^{-((x_i - x_j)^2 + (y_i - y_j)^2 + (D_i + D_j)^2) s^2} + + &\\cdot (erf((u H_j k_{0,imag} + k_{2,imag}) s) + - erf((u H_j k_{0,imag} + k_{2,imag} - H_i) s)) + + &\\cdot e^{(u^2 H_j^2 (k_{0,imag}^2 - 1) + + 2 u H_j (k_{0,imag} k_{2,imag} - k_1) + k_{2,imag}^2) s^2} + du ds + + + k_{0,real} &= + sin(\\beta_j) sin(\\beta_i) cos(\\theta_j - \\theta_i) + + cos(\\beta_j) cos(\\beta_i) + + + k_{0,imag} &= + sin(\\beta_j) sin(\\beta_i) cos(\\theta_j - \\theta_i) + - cos(\\beta_j) cos(\\beta_i) + + + k_{1,real} &= sin(\\beta_j) + (cos(\\theta_j) (x_j - x_i) + sin(\\theta_j) (y_j - y_i)) + + cos(\\beta_j) (D_j - D_i) + + + k_{1,imag} &= sin(\\beta_j) + (cos(\\theta_j) (x_j - x_i) + sin(\\theta_j) (y_j - y_i)) + + cos(\\beta_j) (D_i + D_j) + + + k_{2,real} &= sin(\\beta_i) + (cos(\\theta_i) (x_j - x_i) + sin(\\theta_i) (y_j - y_i)) + + cos(\\beta_i) (D_j - D_i) + + + k_{2,imag} &= sin(\\beta_i) + (cos(\\theta_i) (x_j - x_i) + sin(\\theta_i) (y_j - y_i)) + - cos(\\beta_i) (D_i + D_j) + + where :math:`\\beta_j` and :math:`\\beta_i` are the tilt angle of the + boreholes (relative to vertical), and :math:`\\theta_j` and + :math:`\\theta_i` are the orientation of the boreholes (relative to the + x-axis). + + .. Note:: + The reciprocal thermal response factor + :math:`h_{ji}(t)` can be conveniently calculated by: + + .. math:: + h_{ji}(t) = \\frac{H_i}{H_j} + h_{ij}(t) + + Parameters + ---------- + time : float or (nTimes,) array + Value of time (in seconds) for which the FLS solution is evaluated. + alpha : float + Soil thermal diffusivity (in m2/s). + borefield_j : Borehole or Borefield object + Borehole or Borefield object of the boreholes extracting heat. + borefield_i : Borehole or Borefield object + Borehole or Borefield object object for which the FLS is evaluated. + reaSource : bool + True if the real part of the FLS solution is to be included. + Default is True. + imgSource : bool, optional + True if the image part of the FLS solution is to be included. + Default is True. + outer : bool, optional + True if the finite line source is to be evaluated for all boreholes_j + onto all boreholes_i to return a (nBoreholes_i, nBoreholes_j, nTime,) + array. If false, the finite line source is evaluated pairwise between + boreholes_j and boreholes_i. The numbers of boreholes should be the + same (i.e. nBoreholes_j == nBoreholes_i) and a (nBoreholes, nTime,) + array is returned. + Default is True. + approximation : bool, optional + Set to true to use the approximation of the FLS solution of Cimmino + (2021) [#FLSi-Cimmin2021]_. This approximation does not require + the numerical evaluation of any integral. + Default is False. + M : int, optional + Number of Gauss-Legendre sample points for the quadrature over + :math:`u`. This is only used for inclined boreholes. + Default is 11. + N : int, optional + Number of terms in the approximation of the FLS solution. This + parameter is unused if `approximation` is set to False. + Default is 10. Maximum is 25. + + Returns + ------- + h : float or array, shape (nBoreholes_i, nBoreholes_j, nTime,), (nBoreholes, nTime,) or (nTime,) + Value of the FLS solution. The average (over the length) temperature + drop on the wall of borehole_i due to heat extracted from borehole_j + is: + + .. math:: \\Delta T_{b,i} = T_g - \\frac{Q_j}{2\\pi k_s H_j} h + + Notes + ----- + The function returns a float if time is a float and both of borehole_j and + borehole_i are Borehole objects. If time is a float and any of borehole_j + and borehole_i are Borefield objects, the function returns an array. If + time is an array and both of borehole_j and borehole_i are Borehole + objects, the function returns an 1d array, shape (nTime,). If time is an + array and any of borehole_j and borehole_i are are Borefield objects, the + function returns an array. + + Examples + -------- + >>> b1 = gt.boreholes.Borehole(H=150., D=4., r_b=0.075, x=0., y=0.) + >>> b2 = gt.boreholes.Borehole(H=150., D=4., r_b=0.075, x=5., y=0.) + >>> h = gt.heat_transfer.finite_line_source(4*168*3600., 1.0e-6, b1, b2) + h = 0.0110473635393 + >>> h = gt.heat_transfer.finite_line_source( + 4*168*3600., 1.0e-6, b1, b2, approximation=True, N=10) + h = 0.0110474667731 + >>> b3 = gt.boreholes.Borehole( + H=150., D=4., r_b=0.075, x=5., y=0., tilt=3.1415/15, orientation=0.) + >>> h = gt.heat_transfer.finite_line_source( + 4*168*3600., 1.0e-6, b1, b3, M=21) + h = 0.0002017450051 + + References + ---------- + .. [#FLSi-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. + .. [#FLSi-Lazzar2016] Lazzarotto, A. (2016). A methodology for the + calculation of response functions for geothermal fields with + arbitrarily oriented boreholes – Part 1, Renewable Energy, 86, + 1380-1393. + .. [#FLSi-LazBjo2016] Lazzarotto, A., & Björk, F. (2016). A methodology for + the calculation of response functions for geothermal fields with + arbitrarily oriented boreholes – Part 2, Renewable Energy, 86, + 1353-1361. + + """ + # Check if both bore fields are Borehole objects + single_pair = ( + isinstance(borefield_j, Borehole) + and isinstance(borefield_i, Borehole) + ) + # Convert bore fields to Borefield objects + if isinstance(borefield_j, (Borehole, list)): + borefield_j = Borefield.from_boreholes(borefield_j) + if isinstance(borefield_i, (Borehole, list)): + borefield_i = Borefield.from_boreholes(borefield_i) + # Convert time to array if it is a list + if isinstance(time, list): + time = np.asarray(time) + + # Evaluate the finite line source solution + if time is np.inf: + # Steady-state finite line source solution + h = _finite_line_source_inclined_steady_state( + borefield_j, borefield_i, reaSource=reaSource, + imgSource=imgSource, outer=outer, M=M) + elif approximation: + # Approximation of the finite line source solution + h = _finite_line_source_inclined_approximation( + time, alpha, borefield_j, borefield_i, reaSource=reaSource, + imgSource=imgSource, outer=outer, M=M, N=N) + else: + # Integrand of the finite line source solution + f = _finite_line_source_inclined_integrand( + borefield_j, borefield_i, reaSource=reaSource, + imgSource=imgSource, outer=outer, M=M) + # Evaluate integral + if isinstance(time, (np.floating, float)): + # Lower bound of integration + a = 1. / np.sqrt(4 * alpha * time) + h = (0.5 / borefield_i.H * quad_vec(f, a, np.inf, epsabs=1e-4, epsrel=1e-6)[0].T).T + else: + # Lower bound of integration + a = 1.0 / np.sqrt(4.0 * alpha * time) + # Upper bound of integration + b = np.concatenate(([np.inf], a[:-1])) + h = np.cumsum(np.stack( + [(0.5 / borefield_i.H * quad_vec(f, a_i, b_i, epsabs=1e-4, epsrel=1e-6)[0].T).T + for t, a_i, b_i in zip(time, a, b)], + axis=-1), axis=-1) + + # Return a 1d array if only Borehole objects were provided + if single_pair: + if outer: + h = h[0, 0, ...] + else: + h = h[0, ...] + # Return a float if time is also a float + if isinstance(time, (np.floating, float)): + h = float(h) + return h + + +def _finite_line_source_inclined_approximation( + time: Union[float, np.ndarray], + alpha: float, + borefield_j: Borefield, + borefield_i: Borefield, + reaSource: bool = True, + imgSource: bool = True, + outer: bool = True, + M: int = 11, + N: int = 10 + ) -> np.ndarray: + """ + Evaluate the inclined Finite Line Source (FLS) solution using the + approximation method of Cimmino (2021) [#IncFLSApprox-Cimmin2021]_. + + Parameters + ---------- + time : float or (nTimes,) array + Value of time (in seconds) for which the FLS solution is evaluated. + alpha : float + Soil thermal diffusivity (in m2/s). + borefield_j : Borehole or Borefield object + Borehole or Borefield object of the boreholes extracting heat. + borefield_i : Borehole or Borefield object + Borehole or Borefield object object for which the FLS is evaluated. + reaSource : bool + True if the real part of the FLS solution is to be included. + Default is True. + imgSource : bool, optional + True if the image part of the FLS solution is to be included. + Default is True. + outer : bool, optional + True if the finite line source is to be evaluated for all boreholes_j + onto all boreholes_i to return a (nBoreholes_i, nBoreholes_j, nTime,) + array. If false, the finite line source is evaluated pairwise between + boreholes_j and boreholes_i. The numbers of boreholes should be the + same (i.e. nBoreholes_j == nBoreholes_i) and a (nBoreholes, nTime,) + array is returned. + Default is True. + M : int, optional + Number of Gauss-Legendre sample points for the quadrature over + :math:`u`. This is only used for inclined boreholes. + Default is 11. + N : int, optional + Number of terms in the approximation of the FLS solution. This + parameter is unused if `approximation` is set to False. + Default is 10. Maximum is 25. + + Returns + ------- + h : float or array, shape (nBoreholes_i, nBoreholes_j, nTime,), (nBoreholes, nTime,) or (nTime,) + Value of the FLS solution. The average (over the length) temperature + drop on the wall of borehole_i due to heat extracted from borehole_j + is: + + .. math:: \\Delta T_{b,i} = T_g - \\frac{Q_j}{2\\pi k_s H_j} h + + References + ---------- + .. [#IncFLSApprox-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. + + """ + # Evaluate coefficients of the FLS solution + p, q, k = _finite_line_source_inclined_coefficients( + borefield_j, borefield_i, reaSource=reaSource, imgSource=imgSource, + outer=outer) + + # Coefficients of the approximation of the error function + a, b = _erf_coeffs(N) + + # Roots for Gauss-Legendre quadrature + x, w = roots_legendre(M) + u = 0.5 * x + 0.5 + w = w / 2 + + # Extract lengths and reshape if outer == True + H_j = borefield_j.H + H_i = borefield_i.H + if outer: + H_i = H_i[..., np.newaxis] + H_ratio = np.divide.outer(borefield_j.H, borefield_i.H).T + else: + H_ratio = borefield_j.H / borefield_i.H + if isinstance(time, np.ndarray): + H_ratio = H_ratio[..., np.newaxis] + + # Additional coefficients for the approximation of the FLS solution + s = 1. / (4 * alpha * time) + d = [ + k[2] + np.multiply.outer(u, H_j * k[0]), + k[2] - H_i + np.multiply.outer(u, H_j * k[0]), + ] + c = np.maximum( + (q - d[0]**2 - k[1]**2) + + (k[1] + np.multiply.outer(u, np.ones_like(k[1]) * H_j))**2, + borefield_j.r_b**2) + + # Approximation of the FLS solution + h = 0.25 * H_ratio * (( + np.sign( + np.multiply.outer(d[0], np.ones_like(s)) + ) * exp1( + np.multiply.outer(c + np.multiply.outer(b, d[0]**2), s)) + - np.sign( + np.multiply.outer(d[1], np.ones_like(s)) + ) * exp1(np.multiply.outer(c + np.multiply.outer(b, d[1]**2), s)) + ).T @ a @ w @ p).T + return h + + +def _finite_line_source_inclined_integrand( + borefield_j: Borefield, + borefield_i: Borefield, + reaSource: bool = True, + imgSource: bool = True, + outer: bool = True, + M: int = 11 + ) -> Callable: + """ + Integrand of the inclined Finite Line Source (FLS) solution. + + Parameters + ---------- + borefield_j : Borehole or Borefield object + Borehole or Borefield object of the boreholes extracting heat. + borefield_i : Borehole or Borefield object + Borehole or Borefield object object for which the FLS is evaluated. + reaSource : bool + True if the real part of the FLS solution is to be included. + Default is True. + imgSource : bool, optional + True if the image part of the FLS solution is to be included. + Default is True. + outer : bool, optional + True if the finite line source is to be evaluated for all boreholes_j + onto all boreholes_i to return a (nBoreholes_i, nBoreholes_j, nTime,) + array. If false, the finite line source is evaluated pairwise between + boreholes_j and boreholes_i. The numbers of boreholes should be the + same (i.e. nBoreholes_j == nBoreholes_i) and a (nBoreholes, nTime,) + array is returned. + Default is True. + M : int, optional + Number of Gauss-Legendre sample points for the quadrature over + :math:`u`. This is only used for inclined boreholes. + Default is 11. + + Returns + ------- + f : callable + Integrand of the finite line source solution. Can be vector-valued. + + """ + # Evaluate coefficients of the FLS solution + p, q, k = _finite_line_source_inclined_coefficients( + borefield_j, borefield_i, reaSource=reaSource, imgSource=imgSource, + outer=outer) + + # Roots for Gauss-Legendre quadrature + x, w = roots_legendre(M) + u = 0.5 * x + 0.5 + w = w / 2 + + # Extract lengths and reshape if outer == True + H_j = borefield_j.H + H_i = borefield_i.H + if outer: + H_i = H_i[..., np.newaxis] + + # Integrand of the inclined finite line source solution + f = lambda s: \ + H_j / s * ((( + np.exp( + s**2 * ( + -q + + np.multiply.outer(u**2, H_j**2 * (k[0]**2 - 1)) + + 2 * np.multiply.outer(u, H_j * (k[0] * k[2] - k[1])) + + k[2]**2 + ) + ) \ + * ( + erf( + s * ( + np.multiply.outer(u, H_j * k[0]) + + k[2] + ) + ) + - erf( + s * ( + np.multiply.outer(u, H_j * k[0]) + + k[2] + - H_i + ) + ) + ) + ).T @ w) @ p).T + return f + + +def _finite_line_source_inclined_steady_state( + borefield_j: Borefield, + borefield_i: Borefield, + reaSource: bool = True, + imgSource: bool = True, + outer: bool = True, + M: int = 11 + ) -> np.ndarray: + """ + Steady-state inclined Finite Line Source (FLS) solution. + + Parameters + ---------- + borefield_j : Borehole or Borefield object + Borehole or Borefield object of the boreholes extracting heat. + borefield_i : Borehole or Borefield object + Borehole or Borefield object object for which the FLS is evaluated. + reaSource : bool + True if the real part of the FLS solution is to be included. + Default is True. + imgSource : bool, optional + True if the image part of the FLS solution is to be included. + Default is True. + outer : bool, optional + True if the finite line source is to be evaluated for all boreholes_j + onto all boreholes_i to return a (nBoreholes_i, nBoreholes_j, nTime,) + array. If false, the finite line source is evaluated pairwise between + boreholes_j and boreholes_i. The numbers of boreholes should be the + same (i.e. nBoreholes_j == nBoreholes_i) and a (nBoreholes, nTime,) + array is returned. + Default is True. + M : int, optional + Number of Gauss-Legendre sample points for the quadrature over + :math:`u`. This is only used for inclined boreholes. + Default is 11. + + Returns + ------- + h : float or array, shape (nBoreholes_i, nBoreholes_j) or (nBoreholes,) + Value of the steady-state FLS solution. The average (over the length) + temperature drop on the wall of borehole_i due to heat extracted from + borehole_j is: + + .. math:: \\Delta T_{b,i} = T_g - \\frac{Q_j}{2\\pi k_s H_j} h + + """ + # Evaluate coefficients of the FLS solution + p, q, k = _finite_line_source_inclined_coefficients( + borefield_j, borefield_i, reaSource=reaSource, imgSource=imgSource, + outer=outer) + + # Roots for Gauss-Legendre quadrature + x, w = roots_legendre(M) + u = 0.5 * x + 0.5 + w = w / 2 + + # Extract lengths and reshape if outer == True + H_j = borefield_j.H + H_i = borefield_i.H + if outer: + H_i = H_i[..., np.newaxis] + H_ratio = np.divide.outer(borefield_j.H, borefield_i.H).T + else: + H_ratio = borefield_j.H / borefield_i.H + + # Steady-state inclined finite line source solution + h = 0.5 * H_ratio * ( + np.log( + ( + 2 * H_i * np.sqrt( + q + + np.multiply.outer(u**2, np.ones_like(k[0]) * H_j**2) + + H_i**2 + + 2 * np.multiply.outer(u, H_j * k[1]) + - 2 * H_i * k[2] + - 2 * H_i * np.multiply.outer(u, H_j * k[0]) + ) + + 2 * H_i**2 + - 2 * H_i * k[2] + - 2 * H_i * np.multiply.outer(u, H_j * k[0]) + ) / ( + 2 * H_i * np.sqrt( + q + + np.multiply.outer(u**2, np.ones_like(k[0]) * H_j**2) + + 2 * np.multiply.outer(u, H_j * k[1]) + ) + - 2 * H_i * k[2] + - 2 * H_i * np.multiply.outer(u, H_j * k[0]) + ) + ).T @ w @ p).T + return h + + +def _finite_line_source_inclined_coefficients( + borefield_j: Borefield, + borefield_i: Borefield, + reaSource: bool = True, + imgSource: bool = True, + outer: bool = True + ) -> Tuple[np.ndarray, np.ndarray]: + """ + Coefficients for the finite line source solutions of inclined boreholes. + + Parameters + ---------- + borefield_j : Borehole or Borefield object + Borehole or Borefield object of the boreholes extracting heat. + borefield_i : Borehole or Borefield object + Borehole or Borefield object object for which the FLS is evaluated. + reaSource : bool + True if the real part of the FLS solution is to be included. + Default is True. + imgSource : bool, optional + True if the image part of the FLS solution is to be included. + Default is True. + outer : bool, optional + True if the finite line source is to be evaluated for all boreholes_j + onto all boreholes_i to return a (nBoreholes_i, nBoreholes_j, nTime,) + array. If false, the finite line source is evaluated pairwise between + boreholes_j and boreholes_i. The numbers of boreholes should be the + same (i.e. nBoreholes_j == nBoreholes_i) and a (nBoreholes, nTime,) + array is returned. + Default is True. + + Returns + ------- + p : array + Weights for the superposition of terms in the integrand of the finite + line source solution. + q : array + Squared distance between top edges of line sources. + k : array + Terms used in arguments of exponential and error functions in the + integrand of the inclined finite line source solution. + + """ + # Sines and cosines of tilt (b: beta) and orientation (t: theta) + sb_j = np.sin(borefield_j.tilt) + sb_i = np.sin(borefield_i.tilt) + cb_j = np.cos(borefield_j.tilt) + cb_i = np.cos(borefield_i.tilt) + st_j = np.sin(borefield_j.orientation) + st_i = np.sin(borefield_i.orientation) + ct_j = np.cos(borefield_j.orientation) + ct_i = np.cos(borefield_i.orientation) + # Horizontal distances between boreholes + dis_ij = borefield_j.distance(borefield_i, outer=outer) + if outer: + dx = np.subtract.outer(borefield_j.x, borefield_i.x).T + dy = np.subtract.outer(borefield_j.y, borefield_i.y).T + ct_ij = np.cos( + np.subtract.outer( + borefield_j.orientation, + borefield_i.orientation) + ).T + else: + dx = borefield_j.x - borefield_i.x + dy = borefield_j.y - borefield_i.y + ct_ij = np.cos(borefield_j.orientation - borefield_i.orientation) + + if reaSource and imgSource: + # Full (real + image) FLS solution + p = np.array([1, -1]) + if outer: + dzRea = np.subtract.outer(borefield_j.D, borefield_i.D).T + dzImg = np.add.outer(borefield_i.D, borefield_j.D) + k = [ + np.stack( + [np.multiply.outer(sb_i, sb_j) * ct_ij + np.multiply.outer(cb_i, cb_j), + np.multiply.outer(sb_i, sb_j) * ct_ij - np.multiply.outer(cb_i, cb_j)], + axis=0), + np.stack( + [sb_j * (ct_j * dx + st_j * dy) + cb_j * dzRea, + sb_j * (ct_j * dx + st_j * dy) + cb_j * dzImg], + axis=0), + np.stack( + [(sb_i * (ct_i * dx.T + st_i * dy.T) + cb_i * dzRea.T).T, + (sb_i * (ct_i * dx.T + st_i * dy.T) - cb_i * dzImg.T).T], + axis=0), + ] + else: + dzRea = borefield_j.D - borefield_i.D + dzImg = borefield_j.D + borefield_i.D + k = [ + np.stack( + [sb_j * sb_i * ct_ij + cb_j * cb_i, + sb_j * sb_i * ct_ij - cb_j * cb_i], + axis=0), + np.stack( + [sb_j * (ct_j * dx + st_j * dy) + cb_j * dzRea, + sb_j * (ct_j * dx + st_j * dy) + cb_j * dzImg], + axis=0), + np.stack( + [sb_i * (ct_i * dx + st_i * dy) + cb_i * dzRea, + sb_i * (ct_i * dx + st_i * dy) - cb_i * dzImg], + axis=0), + ] + q = dis_ij**2 + np.stack( + [dzRea, dzImg], + axis=0)**2 + elif reaSource: + # Real FLS solution + p = np.array([1]) + if outer: + dzRea = np.subtract.outer(borefield_j.D, borefield_i.D).T + k = [ + (np.multiply.outer(sb_i, sb_j) * ct_ij + np.multiply.outer(cb_i, cb_j))[np.newaxis, ...], + (sb_j * (ct_j * dx + st_j * dy) + cb_j * dzRea)[np.newaxis, ...], + (sb_i * (ct_i * dx.T + st_i * dy.T) + cb_i * dzRea.T).T[np.newaxis, ...], + ] + else: + dzRea = borefield_j.D - borefield_i.D + k = [ + (sb_j * sb_i * ct_ij + cb_j * cb_i)[np.newaxis, ...], + (sb_j * (ct_j * dx + st_j * dy) + cb_j * dzRea)[np.newaxis, ...], + (sb_i * (ct_i * dx + st_i * dy) + cb_i * dzRea)[np.newaxis, ...], + ] + q = (dis_ij**2 + dzRea**2)[np.newaxis, ...] + elif imgSource: + # Image FLS solution + p = np.array([-1]) + if outer: + dzImg = np.add.outer(borefield_i.D, borefield_j.D) + k = [ + (np.multiply.outer(sb_i, sb_j) * ct_ij - np.multiply.outer(cb_i, cb_j))[np.newaxis, ...], + (sb_j * (ct_j * dx + st_j * dy) + cb_j * dzImg)[np.newaxis, ...], + (sb_i * (ct_i * dx.T + st_i * dy.T) - cb_i * dzImg.T).T[np.newaxis, ...], + ] + else: + dzImg = borefield_j.D + borefield_i.D + k = [ + (sb_j * sb_i * ct_ij - cb_j * cb_i)[np.newaxis, ...], + (sb_j * (ct_j * dx + st_j * dy) + cb_j * dzImg)[np.newaxis, ...], + (sb_i * (ct_i * dx + st_i * dy) - cb_i * dzImg)[np.newaxis, ...], + ] + q = (dis_ij**2 + dzImg**2)[np.newaxis, ...] + else: + # No heat source + p = np.zeros(0) + if outer: + q = np.zeros((0, len(borefield_i), len(borefield_j))) + else: + q = np.zeros((0, len(borefield_i))) + k = [q, q, q] + return p, q, k diff --git a/pygfunction/heat_transfer/finite_line_source_vertical.py b/pygfunction/heat_transfer/finite_line_source_vertical.py new file mode 100644 index 00000000..006f16ac --- /dev/null +++ b/pygfunction/heat_transfer/finite_line_source_vertical.py @@ -0,0 +1,637 @@ +# -*- coding: utf-8 -*- +from collections.abc import Callable +from typing import Union, List, Tuple + +import numpy as np +import numpy.typing as npt +from scipy.integrate import quad_vec +from scipy.special import erfc + +from ..borefield import Borefield +from ..boreholes import Borehole +from ..utilities import erfint, exp1, _erf_coeffs + + +def finite_line_source_vertical( + time: npt.ArrayLike, + alpha: float, + borefield_j: Union[Borehole, Borefield, List[Borehole]], + borefield_i: Union[Borehole, Borefield, List[Borehole]], + distances: Union[None, npt.ArrayLike] = None, + weights: Union[None, npt.ArrayLike] = None, + outer: bool = True, + reaSource: bool = True, + imgSource: bool = True, + approximation: bool = False, + N: int = 10 + ) -> np.ndarray: + """ + Evaluate the Finite Line Source (FLS) solution for vertical boreholes. + + This function uses a numerical quadrature to evaluate the one-integral form + of the FLS solution. For vertical boreholes, the FLS solution was proposed + by Claesson and Javed [#FLSv-ClaJav2011]_ and extended to boreholes with + different vertical positions by Cimmino and Bernier [#FLSv-CimBer2014]_. + The FLS solution is given by: + + .. math:: + h_{ij}(t) &= \\frac{1}{2H_i} + \\int_{\\frac{1}{\\sqrt{4\\alpha t}}}^{\\infty} + e^{-d_{ij}^2s^2}(I_{real}(s)+I_{imag}(s))ds + + + d_{ij} &= \\sqrt{(x_i - x_j)^2 + (y_i - y_j)^2} + + + I_{real}(s) &= erfint((D_i-D_j+H_i)s) - erfint((D_i-D_j)s) + + &+ erfint((D_i-D_j-H_j)s) - erfint((D_i-D_j+H_i-H_j)s) + + I_{imag}(s) &= erfint((D_i+D_j+H_i)s) - erfint((D_i+D_j)s) + + &+ erfint((D_i+D_j+H_j)s) - erfint((D_i+D_j+H_i+H_j)s) + + + erfint(X) &= \\int_{0}^{X} erf(x) dx + + &= Xerf(X) - \\frac{1}{\\sqrt{\\pi}}(1-e^{-X^2}) + + .. Note:: + The reciprocal thermal response factor + :math:`h_{ji}(t)` can be conveniently calculated by: + + .. math:: + h_{ji}(t) = \\frac{H_i}{H_j} + h_{ij}(t) + + Parameters + ---------- + time : float or (nTime,) array + Value of time (in seconds) for which the FLS solution is evaluated. + alpha : float + Soil thermal diffusivity (in m2/s). + borefield_j : Borehole or Borefield object + Borehole or Borefield object of the boreholes extracting heat. + borefield_i : Borehole or Borefield object + Borehole or Borefield object object for which the FLS is evaluated. + distances : float or (nDistances,) array, optional + If None, distances are evaluated from distances between boreholes in + borefield_j and borefield_i. If not None, distances between boreholes + in boreholes_j and boreholes_i are overwritten and an array of shape + (nBoreholes_i, nBoreholes_j, nDistances, nTime,) (if outer == True) + or (nBoreholes, nDistances, nTime,) (if outer == False) is returned. + Default is None. + weights : (nDistances,) or (nSums, nDistances,) array, optional + If None, the FLS solution is unmodified. If not None, the FLS solution + at each distance is multiplied by the corresponding weight and summed + over the axis. An array of shape + (nBoreholes_i, nBoreholes_j, nSums, nTime,) (if outer == True) + or (nBoreholes, nSums, nTime,) (if outer == False) is returned. Only + applies if distances is not None. + Default is None. + reaSource : bool + True if the real part of the FLS solution is to be included. + Default is True. + imgSource : bool, optional + True if the image part of the FLS solution is to be included. + Default is True. + outer : bool, optional + True if the finite line source is to be evaluated for all boreholes_j + onto all boreholes_i to return a (nBoreholes_i, nBoreholes_j, nTime,) + array. If false, the finite line source is evaluated pairwise between + boreholes_j and boreholes_i. The numbers of boreholes should be the + same (i.e. nBoreholes_j == nBoreholes_i) and a (nBoreholes, nTime,) + array is returned. + Default is True. + approximation : bool, optional + Set to true to use the approximation of the FLS solution of Cimmino + (2021) [#FLSv-Cimmin2021]_. This approximation does not require + the numerical evaluation of any integral. + Default is False. + N : int, optional + Number of terms in the approximation of the FLS solution. This + parameter is unused if `approximation` is set to False. + Default is 10. Maximum is 25. + + Returns + ------- + h : float or array, shape (nBoreholes_i, nBoreholes_j, nTime,), (nBoreholes, nTime,) or (nTime,) + Value of the FLS solution. The average (over the length) temperature + drop on the wall of borehole_i due to heat extracted from borehole_j + is: + + .. math:: \\Delta T_{b,i} = T_g - \\frac{Q_j}{2\\pi k_s H_j} h + + Notes + ----- + The function returns a float if time is a float and both of borehole_j and + borehole_i are Borehole objects. If time is a float and any of borehole_j + and borehole_i are Borefield objects, the function returns an array. If + time is an array and both of borehole_j and borehole_i are Borehole + objects, the function returns an 1d array, shape (nTime,). If time is an + array and any of borehole_j and borehole_i are are Borefield objects, the + function returns an array. + + Examples + -------- + >>> b1 = gt.boreholes.Borehole(H=150., D=4., r_b=0.075, x=0., y=0.) + >>> b2 = gt.boreholes.Borehole(H=150., D=4., r_b=0.075, x=5., y=0.) + >>> h = gt.heat_transfer.finite_line_source(4*168*3600., 1.0e-6, b1, b2) + h = 0.0110473635393 + >>> h = gt.heat_transfer.finite_line_source( + 4*168*3600., 1.0e-6, b1, b2, approximation=True, N=10) + h = 0.0110474667731 + >>> b3 = gt.boreholes.Borehole( + H=150., D=4., r_b=0.075, x=5., y=0., tilt=3.1415/15, orientation=0.) + >>> h = gt.heat_transfer.finite_line_source( + 4*168*3600., 1.0e-6, b1, b3, M=21) + h = 0.0002017450051 + + References + ---------- + .. [#FLSv-ClaJav2011] Claesson, J., & Javed, S. (2011). An analytical + method to calculate borehole fluid temperatures for time-scales from + minutes to decades. ASHRAE Transactions, 117(2), 279-288. + .. [#FLSv-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. + .. [#FLSv-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. + + """ + # Check if both bore fields are Borehole objects + single_pair = ( + isinstance(borefield_j, Borehole) + and isinstance(borefield_i, Borehole) + ) + # Convert bore fields to Borefield objects + if isinstance(borefield_j, (Borehole, list)): + borefield_j = Borefield.from_boreholes(borefield_j) + if isinstance(borefield_i, (Borehole, list)): + borefield_i = Borefield.from_boreholes(borefield_i) + # Convert time to array if it is a list + if isinstance(time, list): + time = np.asarray(time) + # Convert distances to array if it is a list + if isinstance(distances, list): + distances = np.asarray(distances) + # Convert weights to array if it is a list + if isinstance(weights, list): + weights = np.asarray(weights) + + # Evaluate the finite line source solution + if time is np.inf: + # Steady-state finite line source solution + h = _finite_line_source_vertical_steady_state( + borefield_j, borefield_i, distances=distances, weights=weights, + outer=outer, reaSource=reaSource, imgSource=imgSource) + elif approximation: + # Approximation of the finite line source solution + h = _finite_line_source_vertical_approximation( + time, alpha, borefield_j, borefield_i, distances=distances, + weights=weights, outer=outer, reaSource=reaSource, + imgSource=imgSource, N=N) + else: + # Integrand of the finite line source solution + f = _finite_line_source_vertical_integrand( + borefield_j, borefield_i, distances=distances, weights=weights, + outer=outer, reaSource=reaSource, imgSource=imgSource) + # Evaluate integral + if isinstance(time, (np.floating, float)): + # Lower bound of integration + a = 1. / np.sqrt(4 * alpha * time) + h = (0.5 / borefield_i.H * quad_vec(f, a, np.inf, epsabs=1e-4, epsrel=1e-6)[0].T).T + else: + # Lower bound of integration + a = 1.0 / np.sqrt(4.0 * alpha * time) + # Upper bound of integration + b = np.concatenate(([np.inf], a[:-1])) + h = np.cumsum(np.stack( + [(0.5 / borefield_i.H * quad_vec(f, a_i, b_i, epsabs=1e-4, epsrel=1e-6)[0].T).T + for t, a_i, b_i in zip(time, a, b)], + axis=-1), axis=-1) + + # Return a 1d array if only Borehole objects were provided + if single_pair: + if outer: + h = h[0, 0, ...] + else: + h = h[0, ...] + # Return a float if time is also a float + if isinstance(time, (np.floating, float)): + h = float(h) + return h + + +def _finite_line_source_vertical_approximation( + time: Union[float, np.ndarray], + alpha: float, + borefield_j: Borefield, + borefield_i: Borefield, + distances: Union[None, np.ndarray] = None, + weights: Union[None, np.ndarray] = None, + reaSource: bool = True, + imgSource: bool = True, + outer: bool = True, + N: int = 10 + ) -> np.ndarray: + """ + Evaluate the Finite Line Source (FLS) solution using the approximation + of Cimmino (2021) [#FLSApprox-Cimmin2021]_. + + Parameters + ---------- + time : float or (nTimes,) array + Value of time (in seconds) for which the FLS solution is evaluated. + alpha : float + Soil thermal diffusivity (in m2/s). + borefield_j : Borehole or Borefield object + Borehole or Borefield object of the boreholes extracting heat. + borefield_i : Borehole or Borefield object + Borehole or Borefield object object for which the FLS is evaluated. + distances : float or (nDistances,) array, optional + If None, distances are evaluated from distances between boreholes in + borefield_j and borefield_i. If not None, distances between boreholes + in boreholes_j and boreholes_i are overwritten and an array of shape + (nBoreholes_i, nBoreholes_j, nDistances, nTime,) (if outer == True) + or (nBoreholes, nDistances, nTime,) (if outer == False) is returned. + Default is None. + weights : (nDistances,) or (nSums, nDistances,) array, optional + If None, the FLS solution is unmodified. If not None, the FLS solution + at each distance is multiplied by the corresponding weight and summed + over the axis. An array of shape + (nBoreholes_i, nBoreholes_j, nSums, nTime,) (if outer == True) + or (nBoreholes, nSums, nTime,) (if outer == False) is returned. Only + applies if distances is not None. + Default is None. + reaSource : bool + True if the real part of the FLS solution is to be included. + Default is True. + imgSource : bool, optional + True if the image part of the FLS solution is to be included. + Default is True. + outer : bool, optional + True if the finite line source is to be evaluated for all boreholes_j + onto all boreholes_i to return a (nBoreholes_i, nBoreholes_j, nTime,) + array. If false, the finite line source is evaluated pairwise between + boreholes_j and boreholes_i. The numbers of boreholes should be the + same (i.e. nBoreholes_j == nBoreholes_i) and a (nBoreholes, nTime,) + array is returned. + Default is True. + N : int, optional + Number of terms in the approximation of the FLS solution. This + parameter is unused if `approximation` is set to False. + Default is 10. Maximum is 25. + + Returns + ------- + h : float or array, shape (nBoreholes_i, nBoreholes_j, nTime,), (nBoreholes, nTime,) or (nTime,) + Value of the FLS solution. The average (over the length) temperature + drop on the wall of borehole_i due to heat extracted from borehole_j + is: + + .. math:: \\Delta T_{b,i} = T_g - \\frac{Q_j}{2\\pi k_s H_j} h + + References + ---------- + .. [#FLSApprox-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. + + """ + # Evaluate coefficients of the FLS solution + p, q = _finite_line_source_vertical_coefficients( + borefield_j, borefield_i, reaSource=reaSource, imgSource=imgSource, + outer=outer) + # The approximation of the error function is only valid for positive + # arguments and f (= x * erf(x)) is an even function + q = np.abs(q) + q2 = q**2 + + # Coefficients of the approximation of the error function + a, b = _erf_coeffs(N) + + # Distances + if distances is None: + dis = borefield_j.distance(borefield_i, outer=outer) + dis2 = dis**2 + sqrt_dis2_plus_q2 = np.sqrt(dis2 + q2) + dis2_plus_b_q2 = dis2 + np.multiply.outer(b, q2) + else: + dis = distances + dis2 = dis**2 + sqrt_dis2_plus_q2 = np.sqrt(np.add.outer(q2, dis2)) + dis2_plus_b_q2 = np.add.outer( + np.multiply.outer(b, q2), + dis2) + + # Temporal terms + four_alpha_time = 4 * alpha * time + sqrt_four_alpha_time = np.sqrt(four_alpha_time) + + # Term G1 of Cimmino (2021) + G1 = 0.5 * (q.T * ( + exp1( + np.divide.outer( + dis2_plus_b_q2, + four_alpha_time + ) + ).T @ a) @ p).T / sqrt_four_alpha_time + + # Term G3 of Cimmino (2021) + x3 = np.divide.outer(sqrt_dis2_plus_q2, sqrt_four_alpha_time) + G3 = ((np.exp(-x3**2) / np.sqrt(np.pi) - x3 * erfc(x3)).T @ p).T + + # Approximation of the FLS solution + h = 0.5 * ((G1 + G3).T / borefield_i.H).T * sqrt_four_alpha_time + if distances is not None and weights is not None: + if isinstance(time, (np.floating, float)): + axis = -1 + else: + axis = -2 + h = np.tensordot(h, weights, axes=(axis, 0)) + return h + + +def _finite_line_source_vertical_coefficients( + borefield_j: Borefield, + borefield_i: Borefield, + reaSource: bool = True, + imgSource: bool = True, + outer: bool = True + ) -> Tuple[np.ndarray, np.ndarray]: + """ + Coefficients for the finite line source solutions of vertical boreholes. + + Parameters + ---------- + borefield_j : Borehole or Borefield object + Borehole or Borefield object of the boreholes extracting heat. + borefield_i : Borehole or Borefield object + Borehole or Borefield object object for which the FLS is evaluated. + reaSource : bool + True if the real part of the FLS solution is to be included. + Default is True. + imgSource : bool, optional + True if the image part of the FLS solution is to be included. + Default is True. + outer : bool, optional + True if the finite line source is to be evaluated for all boreholes_j + onto all boreholes_i to return a (nBoreholes_i, nBoreholes_j, nTime,) + array. If false, the finite line source is evaluated pairwise between + boreholes_j and boreholes_i. The numbers of boreholes should be the + same (i.e. nBoreholes_j == nBoreholes_i) and a (nBoreholes, nTime,) + array is returned. + Default is True. + + Returns + ------- + p : array + Weights for the superposition of terms in the integrand of the finite + line source solution. + q : array + Arguments of the integral of the error function in the integrand of the + finite line source solution. + + """ + H_j = borefield_j.H + D_j = borefield_j.D + H_i = borefield_i.H + D_i = borefield_i.D + + if reaSource and imgSource: + # Full (real + image) FLS solution + p = np.array([1, -1, 1, -1, 1, -1, 1, -1]) + if outer: + q = np.stack( + [np.subtract.outer(D_i + H_i, D_j), + np.subtract.outer(D_i, D_j), + np.subtract.outer(D_i, D_j + H_j), + np.subtract.outer(D_i + H_i, D_j + H_j), + np.add.outer(D_i + H_i, D_j), + np.add.outer(D_i, D_j), + np.add.outer(D_i, D_j + H_j), + np.add.outer(D_i + H_i, D_j + H_j), + ], + axis=0) + else: + q = np.stack( + [D_i + H_i - D_j, + D_i - D_j, + D_i - (D_j + H_j), + D_i + H_i - (D_j + H_j), + D_i + H_i + D_j, + D_i + D_j, + D_i + D_j + H_j, + D_i + H_i + D_j + H_j, + ], + axis=0) + elif reaSource: + # Real FLS solution + p = np.array([1, -1, 1, -1]) + if outer: + q = np.stack( + [np.subtract.outer(D_i + H_i, D_j), + np.subtract.outer(D_i, D_j), + np.subtract.outer(D_i, D_j + H_j), + np.subtract.outer(D_i + H_i, D_j + H_j), + ], + axis=0) + else: + q = np.stack( + [D_i + H_i - D_j, + D_i - D_j, + D_i - (D_j + H_j), + D_i + H_i - (D_j + H_j), + ], + axis=0) + elif imgSource: + # Image FLS solution + p = np.array([1, -1, 1, -1]) + if outer: + q = np.stack( + [np.add.outer(D_i + H_i, D_j), + np.add.outer(D_i, D_j), + np.add.outer(D_i, D_j + H_j), + np.add.outer(D_i + H_i, D_j + H_j), + ], + axis=0) + else: + q = np.stack( + [D_i + H_i + D_j, + D_i + D_j, + D_i + D_j + H_j, + D_i + H_i + D_j + H_j, + ], + axis=0) + else: + # No heat source + p = np.zeros(0) + if outer: + q = np.zeros((0, len(borefield_i), len(borefield_j))) + else: + q = np.zeros((0, len(borefield_i))) + return p, q + + +def _finite_line_source_vertical_integrand( + borefield_j: Borefield, + borefield_i: Borefield, + distances: Union[None, np.ndarray] = None, + weights: Union[None, np.ndarray] = None, + reaSource: bool = True, + imgSource: bool = True, + outer: bool = True + ) -> Callable: + """ + Integrand of the finite line source solution. + + Parameters + ---------- + borefield_j : Borehole or Borefield object + Borehole or Borefield object of the boreholes extracting heat. + borefield_i : Borehole or Borefield object + Borehole or Borefield object object for which the FLS is evaluated. + distances : float or (nDistances,) array, optional + If None, distances are evaluated from distances between boreholes in + borefield_j and borefield_i. If not None, distances between boreholes + in boreholes_j and boreholes_i are overwritten and an array of shape + (nBoreholes_i, nBoreholes_j, nDistances, nTime,) (if outer == True) + or (nBoreholes, nDistances, nTime,) (if outer == False) is returned. + Default is None. + weights : (nDistances,) or (nSums, nDistances,) array, optional + If None, the FLS solution is unmodified. If not None, the FLS solution + at each distance is multiplied by the corresponding weight and summed + over the axis. An array of shape + (nBoreholes_i, nBoreholes_j, nSums, nTime,) (if outer == True) + or (nBoreholes, nSums, nTime,) (if outer == False) is returned. Only + applies if distances is not None. + Default is None. + reaSource : bool + True if the real part of the FLS solution is to be included. + Default is True. + imgSource : bool, optional + True if the image part of the FLS solution is to be included. + Default is True. + outer : bool, optional + True if the finite line source is to be evaluated for all boreholes_j + onto all boreholes_i to return a (nBoreholes_i, nBoreholes_j, nTime,) + array. If false, the finite line source is evaluated pairwise between + boreholes_j and boreholes_i. The numbers of boreholes should be the + same (i.e. nBoreholes_j == nBoreholes_i) and a (nBoreholes, nTime,) + array is returned. + Default is True. + + Returns + ------- + f : callable + Integrand of the finite line source solution. Can be vector-valued. + + """ + # Evaluate coefficients of the FLS solution + p, q = _finite_line_source_vertical_coefficients( + borefield_j, borefield_i, reaSource=reaSource, imgSource=imgSource, + outer=outer) + # Integrand of the finite line source solution + if distances is None: + dis = borefield_j.distance(borefield_i, outer=outer) + f = lambda s: s**-2 * np.exp(-dis**2 * s**2) * (erfint(q * s).T @ p).T + else: + dis = distances + if weights is None: + f = lambda s: s**-2 * np.multiply.outer( + (erfint(q * s).T @ p).T, np.exp(-dis**2 * s**2) + ) + else: + w = weights + f = lambda s: s**-2 * np.multiply.outer( + (erfint(q * s).T @ p).T, w @ np.exp(-dis**2 * s**2) + ) + return f + + +def _finite_line_source_vertical_steady_state( + borefield_j: Borefield, + borefield_i: Borefield, + distances: Union[None, float, np.ndarray] = None, + weights: Union[None, np.ndarray] = None, + reaSource: bool = True, + imgSource: bool = True, + outer: bool = True + ) -> np.ndarray: + """ + Steady-state finite line source solution. + + Parameters + ---------- + borefield_j : Borehole or Borefield object + Borehole or Borefield object of the boreholes extracting heat. + borefield_i : Borehole or Borefield object + Borehole or Borefield object object for which the FLS is evaluated. + distances : float or (nDistances,) array, optional + If None, distances are evaluated from distances between boreholes in + borefield_j and borefield_i. If not None, distances between boreholes + in boreholes_j and boreholes_i are overwritten and an array of shape + (nBoreholes_i, nBoreholes_j, nDistances, nTime,) (if outer == True) + or (nBoreholes, nDistances, nTime,) (if outer == False) is returned. + Default is None. + weights : (nDistances,) or (nSums, nDistances,) array, optional + If None, the FLS solution is unmodified. If not None, the FLS solution + at each distance is multiplied by the corresponding weight and summed + over the axis. An array of shape + (nBoreholes_i, nBoreholes_j, nSums, nTime,) (if outer == True) + or (nBoreholes, nSums, nTime,) (if outer == False) is returned. Only + applies if distances is not None. + Default is None. + reaSource : bool + True if the real part of the FLS solution is to be included. + Default is True. + imgSource : bool, optional + True if the image part of the FLS solution is to be included. + Default is True. + outer : bool, optional + True if the finite line source is to be evaluated for all boreholes_j + onto all boreholes_i to return a (nBoreholes_i, nBoreholes_j, nTime,) + array. If false, the finite line source is evaluated pairwise between + boreholes_j and boreholes_i. The numbers of boreholes should be the + same (i.e. nBoreholes_j == nBoreholes_i) and a (nBoreholes, nTime,) + array is returned. + Default is True. + + Returns + ------- + h : float or array, shape (nBoreholes_i, nBoreholes_j) or (nBoreholes,) + Value of the steady-state FLS solution. The average (over the length) + temperature drop on the wall of borehole_i due to heat extracted from + borehole_j is: + + .. math:: \\Delta T_{b,i} = T_g - \\frac{Q_j}{2\\pi k_s H_j} h + + """ + # Evaluate coefficients of the FLS solution + p, q = _finite_line_source_vertical_coefficients( + borefield_j, borefield_i, reaSource=reaSource, imgSource=imgSource, + outer=outer) + + # Extract lengths and reshape if outer == True + H_i = borefield_i.H + + # Steady-state finite line source solution + if distances is None: + dis = borefield_j.distance(borefield_i, outer=outer) + q2_plus_dis2 = np.sqrt(q**2 + dis**2) + log_q_q2_plus_dis2 = np.log(q + q2_plus_dis2) + else: + dis = distances + q2_plus_dis2 = np.sqrt(np.add.outer(q**2, dis**2)) + if isinstance(dis, np.ndarray): + q = q[..., np.newaxis] + log_q_q2_plus_dis2 = np.log(q + q2_plus_dis2) + if weights is not None: + q2_plus_dis2 = q2_plus_dis2 @ weights.T + log_q_q2_plus_dis2 = log_q_q2_plus_dis2 @ weights.T + h = 0.5 * (((q * log_q_q2_plus_dis2 - q2_plus_dis2).T @ p) / H_i).T + return h From a9c474be4f174859458f3abbcc4a0350e678b30c Mon Sep 17 00:00:00 2001 From: MassimoCimmino Date: Mon, 28 Apr 2025 15:56:47 -0400 Subject: [PATCH 2/7] Update heat_transfer_test --- tests/heat_transfer_test.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/tests/heat_transfer_test.py b/tests/heat_transfer_test.py index 27cd1954..419acf8d 100644 --- a/tests/heat_transfer_test.py +++ b/tests/heat_transfer_test.py @@ -387,7 +387,7 @@ def test_finite_line_source_Lazzarotto_steady_state( (np.array([0.1, 1., 10.]), True, False, False, 21, 10, np.array( [5.63802126, 6.53283215, 7.03111215])), (np.array([0.1, 1., 10.]), False, True, False, 21, 10, np.array( - [-0.0230257 , -0.14490784, -0.41781906])), + [-0.0230253, -0.14490724, -0.41781839])), (np.array([0.1, 1., 10.]), False, False, False, 21, 10, np.array( [0., 0., 0.])), # Multiple time steps - FLS approximation @@ -396,7 +396,7 @@ def test_finite_line_source_Lazzarotto_steady_state( (np.array([0.1, 1., 10.]), True, False, True, 21, 10, np.array( [5.63803137, 6.53285537, 7.03118607])), (np.array([0.1, 1., 10.]), False, True, True, 21, 10, np.array( - [-0.02302611, -0.14490236, -0.41776766])), + [-0.02302571, -0.14490175, -0.41776699])), (np.array([0.1, 1., 10.]), False, False, True, 21, 10, np.array( [0., 0., 0.])), ]) @@ -479,12 +479,12 @@ def test_finite_line_source_inclined_to_self_steady_state( [[7.03111215, 0.9284482], [0.9284482, 7.03111215]]]).transpose((1, 2, 0))), (np.array([0.1, 1., 10.]), False, True, False, 21, 10, np.array( - [[[-0.02037585, -0.01111298], - [-0.01111298, -0.02037585]], - [[-0.1293867, -0.10593651], - [-0.10593651, -0.1293867]], - [[-0.38493416, -0.34716039], - [-0.34716039, -0.38493416]]]).transpose((1, 2, 0))), + [[[-0.02037549, -0.01111298], + [-0.01111298, -0.02037549]], + [[-0.12938617, -0.10593651], + [-0.10593651, -0.12938617]], + [[-0.38493357, -0.34716039], + [-0.34716039, -0.38493357]]]).transpose((1, 2, 0))), (np.array([0.1, 1., 10.]), False, False, False, 21, 10, np.zeros((2, 2, 3))), # Multiple time steps - FLS approximation (np.array([0.1, 1., 10.]), True, True, True, 21, 10, np.array( @@ -502,12 +502,12 @@ def test_finite_line_source_inclined_to_self_steady_state( [[7.03118607, 0.9283558], [0.9283558, 7.03118607]]]).transpose((1, 2, 0))), (np.array([0.1, 1., 10.]), False, True, True, 21, 10, np.array( - [[[-0.02037611, -0.01111366], - [-0.01111366, -0.02037611]], - [[-0.12938214, -0.10592953], - [-0.10592953, -0.12938214]], - [[-0.38489415, -0.34712584], - [-0.34712584, -0.38489415]]]).transpose((1, 2, 0))), + [[[-0.02037576, -0.01111366], + [-0.01111366, -0.02037576]], + [[-0.1293816, -0.10592953], + [-0.10592953, -0.1293816]], + [[-0.38489356, -0.34712584], + [-0.34712584, -0.38489356]]]).transpose((1, 2, 0))), (np.array([0.1, 1., 10.]), False, False, True, 21, 10, np.zeros((2, 2, 3))), ]) def test_finite_line_source_multiple_inclined_to_multiple_inclined( From f7ff3f5da56522d06f23d15a7c9375afbe2bd8c9 Mon Sep 17 00:00:00 2001 From: MassimoCimmino Date: Tue, 29 Apr 2025 05:34:48 -0400 Subject: [PATCH 3/7] Modify solvers to use the new heat_transfer functions --- pygfunction/borefield.py | 198 +++++++++++- pygfunction/boreholes.py | 178 +++++++---- pygfunction/gfunction.py | 154 ++++++---- pygfunction/networks.py | 6 + pygfunction/solvers/_base_solver.py | 434 ++++++++++++-------------- pygfunction/solvers/detailed.py | 111 +++---- pygfunction/solvers/equivalent.py | 253 ++++++++------- pygfunction/solvers/similarities.py | 457 +++++++++++++++++++--------- pygfunction/utilities.py | 2 +- 9 files changed, 1098 insertions(+), 695 deletions(-) diff --git a/pygfunction/borefield.py b/pygfunction/borefield.py index dc5bf5a8..2019ce33 100644 --- a/pygfunction/borefield.py +++ b/pygfunction/borefield.py @@ -6,7 +6,7 @@ import numpy as np import numpy.typing as npt -from .boreholes import Borehole +from .boreholes import Borehole, _EquivalentBorehole from .utilities import _initialize_figure, _format_axes, _format_axes_3d @@ -44,10 +44,9 @@ def __init__( self, H: npt.ArrayLike, D: npt.ArrayLike, r_b: npt.ArrayLike, x: npt.ArrayLike, y: npt.ArrayLike, tilt: npt.ArrayLike = 0., orientation: npt.ArrayLike = 0.): - # Convert x and y coordinates to arrays - x = np.atleast_1d(x) - y = np.atleast_1d(y) - self.nBoreholes = np.maximum(len(x), len(y)) + self.nBoreholes = np.max( + [len(np.atleast_1d(var)) + for var in (H, D, r_b, x, y, tilt, orientation)]) # Broadcast all variables to arrays of length `nBoreholes` self.H = np.broadcast_to(H, self.nBoreholes) @@ -357,6 +356,77 @@ def evaluate_g_function( return gfunc.gFunc + def segments( + self, nSegments: Union[int, npt.ArrayLike], + segment_ratios: Union[ + None, npt.ArrayLike, List[npt.ArrayLike], Callable[[int], npt.ArrayLike] + ] = None) -> Self: + """ + Split boreholes in the bore field into segments. + + Parameters + ---------- + nSegments : int or (nBoreholes,) array + Number of segments per borehole. + segment_ratios : array or list of arrays, 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,). If all boreholes have the same number of segments + and the same ratios, a single (nSegments,) array can be provided. + Otherwise, a list of (nSegments_i,) arrays with the ratios + associated with each borehole must be provided. If + segment_ratios==None, segments of equal lengths are considered. + Default is None. + + Returns + ------- + segments : Borefield object + Segments in the bore field. + + Examples + -------- + >>> borefield = gt.borefield.Borefield.rectangle_field( + N_1=2, N_2=3, B_1 = 7.5, B_2=7.5, H=150., D=4., r_b=0.075) + >>> borefield.segments(5) + + """ + nSegments = np.broadcast_to(nSegments, self.nBoreholes) + if callable(segment_ratios): + segment_ratios = [segment_ratios(nSeg) for nSeg in nSegments] + if segment_ratios is None: + # Local coordinates of the top-edges of the segments + xi = np.concatenate( + [np.arange(nSeg_i) / nSeg_i for nSeg_i in nSegments]) + # Segment ratios + segment_ratios = np.concatenate( + [np.full(nSeg_i, 1. / nSeg_i) for nSeg_i in nSegments]) + elif isinstance(segment_ratios, list): + # Local coordinates of the top-edges of the segments + xi = np.concatenate( + [np.cumsum(np.concatenate([[0.], seg_rat_i[:-1]])) + for seg_rat_i in segment_ratios]) + # Segment ratios + segment_ratios = np.concatenate(segment_ratios) + else: + # Local coordinates of the top-edges of the segments + xi = np.tile( + np.cumsum(np.concatenate([[0.], segment_ratios[:-1]])), + self.nBoreholes) + # Segment ratios + segment_ratios = np.tile(segment_ratios, self.nBoreholes) + xi = xi * np.repeat(self.H, nSegments) + + # Geometrical parameters of the segments + H = np.repeat(self.H, nSegments) * segment_ratios + r_b = np.repeat(self.r_b, nSegments) + tilt = np.repeat(self.tilt, nSegments) + orientation = np.repeat(self.orientation, nSegments) + D = np.repeat(self.D, nSegments) + xi * np.cos(tilt) + x = np.repeat(self.x, nSegments) + xi * np.sin(tilt) * np.cos(orientation) + y = np.repeat(self.y, nSegments) + xi * np.sin(tilt) * np.sin(orientation) + + return Borefield(H, D, r_b, x, y, tilt=tilt, orientation=orientation) + def visualize_field( self, viewTop: bool = True, view3D: bool = True, labels: bool = True, showTilt: bool = True): @@ -1151,3 +1221,121 @@ def circle_field( # Create the bore field borefield = cls(H, D, r_b, x, y, tilt=tilt, orientation=orientation) return borefield + + +class _EquivalentBorefield(object): + """ + Contains information regarding the dimensions and positions of boreholes within a borefield. + + Attributes + ---------- + H : (nEqBoreholes,) array + Borehole lengths (in meters). + D : (nEqBoreholes,) array + Borehole buried depths (in meters). + r_b : (nEqBoreholes,) array + Borehole radii (in meters). + x : (nEqBoreholes,) list of (nBoreholes_i,) arrays + Position (in meters) of the head of the boreholes along the x-axis. + y : (nEqBoreholes,) list of (nBoreholes_i,) arrays + Position (in meters) of the head of the boreholes along the y-axis. + tilt : (nEqBoreholes,) array, optional + Angle (in radians) from vertical of the axis of the boreholes. + Default is 0. + orientation : (nEqBoreholes,) list of (nBoreholes_i,) arrays, optional + Direction (in radians) of the tilt of the boreholes. Defaults to zero + if the borehole is vertical. + Default is 0. + + """ + + def __init__( + self, H: npt.ArrayLike, D: npt.ArrayLike, r_b: npt.ArrayLike, + x: List[npt.ArrayLike], y: List[npt.ArrayLike], + tilt: npt.ArrayLike, orientation: List[npt.ArrayLike]): + self.H = H + self.D = D + self.r_b = r_b + self.x = x + self.y = y + self.tilt = tilt + self.orientation = orientation + self.nBoreholes = len(x) + + # Identify tilted boreholes + self._is_tilted = np.broadcast_to( + np.greater(np.abs(tilt), 1e-6), + self.nBoreholes) + # Vertical boreholes default to an orientation of zero + if not np.any(self._is_tilted): + self.orientation = [np.broadcast_to(0., len(x_i)) for x_i in x] + else: + self.orientation = [ + np.multiply(orientation_i, _is_tilted_i) + for (orientation_i, _is_tilted_i) + in zip(self.orientation, self._is_tilted)] + + def __getitem__(self, key): + if isinstance(key, (int, np.integer)): + # Returns a borehole object if only one borehole is indexed + output_class = _EquivalentBorehole + else: + # Returns a _EquivalentBorehole object for slices and lists of + # indexes + output_class = _EquivalentBorefield + return output_class( + self.H[key], self.D[key], self.r_b[key], self.x[key], self.y[key], + tilt=self.tilt[key], orientation=self.orientation[key]) + + def __len__(self) -> int: + """Return the number of equivalent boreholes.""" + return self.nBoreholes + + def segments( + self, nSegments: int, + segment_ratios: Union[ + None, npt.ArrayLike, Callable[[int], npt.ArrayLike] + ] = None) -> Self: + """ + Split boreholes in the bore field into segments. + + Parameters + ---------- + nSegments : int or (nBoreholes,) array + Number of segments per borehole. + segment_ratios : array or list of arrays, 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,). If all boreholes have the same number of segments + and the same ratios, a single (nSegments,) array can be provided. + Otherwise, a list of (nSegments_i,) arrays with the ratios + associated with each borehole must be provided. If + segment_ratios==None, segments of equal lengths are considered. + Default is None. + + Returns + ------- + segments : Borefield object + Segments in the bore field. + + Examples + -------- + >>> borefield = gt.borefield.Borefield.rectangle_field( + N_1=2, N_2=3, B_1 = 7.5, B_2=7.5, H=150., D=4., r_b=0.075) + >>> borefield.segments(5) + + """ + return _EquivalentBorefield.from_equivalent_boreholes( + [seg for b in self for seg in b.segments(nSegments, segment_ratios=segment_ratios) + ]) + + @classmethod + def from_equivalent_boreholes(cls, equivalent_boreholes): + H = np.array([b.H for b in equivalent_boreholes]) + D = np.array([b.D for b in equivalent_boreholes]) + r_b = np.array([b.r_b for b in equivalent_boreholes]) + x = [b.x for b in equivalent_boreholes] + y = [b.y for b in equivalent_boreholes] + tilt = np.array([b.tilt for b in equivalent_boreholes]) + orientation = [b.orientation for b in equivalent_boreholes] + return cls(H, D, r_b, x, y, tilt, orientation) diff --git a/pygfunction/boreholes.py b/pygfunction/boreholes.py index 80af83ea..2b6a5d5a 100644 --- a/pygfunction/boreholes.py +++ b/pygfunction/boreholes.py @@ -160,6 +160,8 @@ def segments(self, nSegments, segment_ratios=None): """ if segment_ratios is None: segment_ratios = np.full(nSegments, 1. / nSegments) + elif callable(segment_ratios): + segment_ratios = segment_ratios(nSegments) z = self._segment_edges(nSegments, segment_ratios=segment_ratios)[:-1] boreSegments = [] for z_i, ratios in zip(z, segment_ratios): @@ -208,6 +210,8 @@ def _segment_edges(self, nSegments, segment_ratios=None): """ if segment_ratios is None: segment_ratios = np.full(nSegments, 1. / nSegments) + elif callable(segment_ratios): + segment_ratios = segment_ratios(nSegments) z = np.concatenate(([0.], np.cumsum(segment_ratios))) * self.H return z @@ -241,6 +245,8 @@ def _segment_midpoints(self, nSegments, segment_ratios=None): """ if segment_ratios is None: segment_ratios = np.full(nSegments, 1. / nSegments) + elif callable(segment_ratios): + segment_ratios = segment_ratios(nSegments) z = self._segment_edges(nSegments, segment_ratios=segment_ratios)[:-1] \ + segment_ratios * self.H / 2 return z @@ -290,36 +296,70 @@ class _EquivalentBorehole(object): Simulation, 14 (4), 446-460. """ - def __init__(self, boreholes): - if isinstance(boreholes[0], Borehole): - self.H = boreholes[0].H - self.D = boreholes[0].D - self.r_b = boreholes[0].r_b - self.x = np.array([b.x for b in boreholes]) - self.y = np.array([b.y for b in boreholes]) - self.tilt = boreholes[0].tilt - self.orientation = np.array([b.orientation for b in boreholes]) - elif isinstance(boreholes[0], _EquivalentBorehole): - self.H = boreholes[0].H - self.D = boreholes[0].D - self.r_b = boreholes[0].r_b - self.x = np.concatenate([b.x for b in boreholes]) - self.y = np.concatenate([b.y for b in boreholes]) - self.tilt = boreholes[0].tilt - self.orientation = np.concatenate( - [b.orientation for b in boreholes]) - elif type(boreholes) is tuple: - self.H, self.D, self.r_b, self.x, self.y = boreholes[:5] - self.x = np.atleast_1d(self.x) - self.y = np.atleast_1d(self.y) - if len(boreholes)==7: - self.tilt, self.orientation = boreholes[5:] - - self.nBoreholes = len(self.x) - # Check if borehole is inclined - self._is_tilted = np.abs(self.tilt) > 1.0e-6 + def __init__(self, H, D, r_b, x, y, tilt=0., orientation=0.): + self.nBoreholes = np.maximum(len(x), len(y)) + + # Broadcast all variables to arrays of length `nBoreholes` + self.H = H + self.D = D + self.r_b = r_b + self.x = x + self.y = y + self.tilt = tilt + + # Identify tilted boreholes + self._is_tilted = np.abs(tilt) > 1e-6 + # Vertical boreholes default to an orientation of zero + if not self._is_tilted: + self.orientation = np.broadcast_to(0., self.nBoreholes) + elif np.all(self._is_tilted): + self.orientation = np.broadcast_to(orientation, self.nBoreholes) + else: + self.orientation = np.multiply(orientation, self._is_tilted) - def distance(self, target): + def __getitem__(self, key): + if isinstance(key, (int, np.integer)): + # Returns a borehole object if only one borehole is indexed + output_class = Borehole + else: + # Returns a _EquivalentBorehole object for slices and lists of + # indexes + output_class = _EquivalentBorehole + return output_class( + self.H, self.D, self.r_b, self.x[key], self.y[key], + tilt=self.tilt, orientation=self.orientation[key]) + + def __len__(self) -> int: + """Return the number of boreholes.""" + return self.nBoreholes + + @property + def is_tilted(self): + """ + Returns true if the borehole is inclined. + + Returns + ------- + bool + True if borehole is inclined. + + """ + return self._is_tilted + + @property + def is_vertical(self): + """ + Returns true if the borehole is vertical. + + Returns + ------- + bool + True if borehole is vertical. + + """ + return not self._is_tilted + + def distance(self, other_borehole: Self) -> np.ndarray: """ Evaluate the distance between the current borehole and a target borehole. @@ -349,34 +389,12 @@ def distance(self, target): """ dis = np.maximum( np.sqrt( - np.add.outer(target.x, -self.x)**2 + np.add.outer(target.y, -self.y)**2), + np.add.outer(other_borehole.x, -self.x)**2 + + np.add.outer(other_borehole.y, -self.y)**2 + ), self.r_b) return dis - def is_tilted(self): - """ - Returns true if the borehole is inclined. - - Returns - ------- - bool - True if borehole is inclined. - - """ - return self._is_tilted - - def is_vertical(self): - """ - Returns true if the borehole is vertical. - - Returns - ------- - bool - True if borehole is vertical. - - """ - return not self._is_tilted - def position(self): """ Returns the position of the boreholes represented by the equivalent @@ -423,15 +441,17 @@ def segments(self, nSegments, segment_ratios=None): """ if segment_ratios is None: segment_ratios = np.full(nSegments, 1. / nSegments) + elif callable(segment_ratios): + segment_ratios = segment_ratios(nSegments) z = self._segment_edges(nSegments, segment_ratios=segment_ratios)[:-1] segments = [_EquivalentBorehole( - (ratios * self.H, - self.D + z_i * np.cos(self.tilt), - self.r_b, - self.x + z_i * np.sin(self.tilt) * np.cos(self.orientation), - self.y + z_i * np.sin(self.tilt) * np.sin(self.orientation), - self.tilt, - self.orientation) + ratios * self.H, + self.D + z_i * np.cos(self.tilt), + self.r_b, + self.x + z_i * np.sin(self.tilt) * np.cos(self.orientation), + self.y + z_i * np.sin(self.tilt) * np.sin(self.orientation), + tilt=self.tilt, + orientation=self.orientation ) for z_i, ratios in zip(z, segment_ratios)] return segments @@ -527,6 +547,8 @@ def _segment_edges(self, nSegments, segment_ratios=None): """ if segment_ratios is None: segment_ratios = np.full(nSegments, 1. / nSegments) + elif callable(segment_ratios): + segment_ratios = segment_ratios(nSegments) z = np.concatenate(([0.], np.cumsum(segment_ratios))) * self.H return z @@ -560,10 +582,44 @@ def _segment_midpoints(self, nSegments, segment_ratios=None): """ if segment_ratios is None: segment_ratios = np.full(nSegments, 1. / nSegments) + elif callable(segment_ratios): + segment_ratios = segment_ratios(nSegments) z = self._segment_edges(nSegments, segment_ratios=segment_ratios)[:-1] \ + segment_ratios * self.H / 2 return z + @classmethod + def from_borefield(cls, borefield): + assert np.allclose(borefield.H, borefield.H[0], rtol=1e-6) + assert np.allclose(borefield.D, borefield.D[0], rtol=1e-6) + assert np.allclose(borefield.r_b, borefield.r_b[0], rtol=1e-6) + assert np.allclose(borefield.tilt, borefield.tilt[0], rtol=1e-6) + H = borefield.H[0] + D = borefield.D[0] + r_b = borefield.r_b[0] + x = borefield.x + y = borefield.y + tilt = borefield.tilt[0] + orientation = borefield.orientation + return cls(H, D, r_b, x, y, tilt=tilt, orientation=orientation) + + @classmethod + def from_boreholes(cls, boreholes): + assert np.allclose([b.H for b in boreholes], boreholes[0].H, rtol=1e-6) + assert np.allclose([b.D for b in boreholes], boreholes[0].D, rtol=1e-6) + assert np.allclose( + [b.r_b for b in boreholes], boreholes[0].r_b, rtol=1e-6) + assert np.allclose( + [b.tilt for b in boreholes], boreholes[0].tilt, rtol=1e-6) + H = boreholes[0].H + D = boreholes[0].D + r_b = boreholes[0].r_b + x = np.array([b.x for b in boreholes]) + y = np.array([b.y for b in boreholes]) + tilt = boreholes[0].tilt + orientation = np.array([b.orientation for b in boreholes]) + return cls(H, D, r_b, x, y, tilt=tilt, orientation=orientation) + def find_duplicates(boreField, disp=False): """ diff --git a/pygfunction/gfunction.py b/pygfunction/gfunction.py index 658a2f5c..27fc3fcb 100644 --- a/pygfunction/gfunction.py +++ b/pygfunction/gfunction.py @@ -1,12 +1,15 @@ # -*- coding: utf-8 -*- -import warnings from time import perf_counter +import warnings import numpy as np from scipy.interpolate import interp1d as interp1d from .borefield import Borefield -from .boreholes import Borehole, find_duplicates +from .boreholes import ( + Borehole, + find_duplicates + ) from .networks import Network from .solvers import ( Detailed, @@ -14,10 +17,10 @@ Similarities ) from .utilities import ( - segment_ratios, _initialize_figure, _format_axes ) +from . import utilities class gFunction(object): @@ -33,9 +36,8 @@ class gFunction(object): Attributes ---------- - boreholes_or_network : list of Borehole objects or Network object - List of boreholes included in the bore field, or network of boreholes - and pipes. + borefield_or_network : Borefield object or Network object + The bore field, or network of boreholes and pipes. alpha : float Soil thermal diffusivity (in m2/s). time : float or array, optional @@ -241,7 +243,7 @@ class gFunction(object): rate and reversible flow direction. Renewable Energy, 228, 120661. """ - def __init__(self, boreholes_or_network, alpha, time=None, + def __init__(self, borefield_or_network, alpha, time=None, method='equivalent', boundary_condition=None, m_flow_borehole=None, m_flow_network=None, cp_f=None, options={}): @@ -255,24 +257,24 @@ def __init__(self, boreholes_or_network, alpha, time=None, self.options = options # Format inputs and assign default values where needed - self._format_inputs(boreholes_or_network) + self._format_inputs(borefield_or_network) # Check the validity of inputs self._check_inputs() # Load the chosen solver if self.method.lower()=='similarities': self.solver = Similarities( - self.boreholes, self.network, self.time, + self.borefield, 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.boreholes, self.network, self.time, + self.borefield, 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.boreholes, self.network, self.time, + self.borefield, self.network, self.time, self.boundary_condition, self.m_flow_borehole, self.m_flow_network, self.cp_f, **self.options) else: @@ -347,7 +349,7 @@ def visualize_g_function(self, which=None): _format_axes(ax) # Borefield characteristic time - ts = np.mean([b.H for b in self.boreholes])**2/(9.*self.alpha) + ts = self.borefield.H.mean()**2 / (9 * self.alpha) # Dimensionless time (log) lntts = np.log(self.time/ts) # Draw g-function @@ -408,11 +410,11 @@ def visualize_heat_extraction_rates( # If iBoreholes is None, then plot all boreholes if iBoreholes is None: - iBoreholes = range(len(self.solver.boreholes)) + iBoreholes = range(len(self.solver.borefield)) # Import heat extraction rates Q_t = self._heat_extraction_rates(iBoreholes) # Borefield characteristic time - ts = np.mean([b.H for b in self.solver.boreholes])**2/(9.*self.alpha) + ts = self.borefield.H.mean()**2 / (9 * self.alpha) # Dimensionless time (log) lntts = np.log(self.time/ts) @@ -430,7 +432,7 @@ def visualize_heat_extraction_rates( _format_axes(ax2) # Plot curves for requested boreholes - for i, borehole in enumerate(self.solver.boreholes): + for i, borehole in enumerate(self.solver.borefield): if i in iBoreholes: # Draw heat extraction rate line = ax2.plot(lntts, Q_t[iBoreholes.index(i)]) @@ -485,7 +487,7 @@ def visualize_heat_extraction_rates( _format_axes(ax2) # Plot curves for requested boreholes - for i, borehole in enumerate(self.solver.boreholes): + for i, borehole in enumerate(self.solver.borefield): if i in iBoreholes: # Draw heat extraction rate line = ax2.plot(lntts, Q_t[iBoreholes.index(i)][n]) @@ -559,7 +561,7 @@ def visualize_heat_extraction_rate_profiles( # If iBoreholes is None, then plot all boreholes if iBoreholes is None: - iBoreholes = range(len(self.solver.boreholes)) + iBoreholes = range(len(self.solver.borefield)) # Import heat extraction rate profiles z, Q_b = self._heat_extraction_rate_profiles(time, iBoreholes) @@ -578,7 +580,7 @@ def visualize_heat_extraction_rate_profiles( _format_axes(ax2) # Plot curves for requested boreholes - for i, borehole in enumerate(self.solver.boreholes): + for i, borehole in enumerate(self.solver.borefield): if i in iBoreholes: # Draw heat extraction rate profile line = ax2.plot( @@ -635,7 +637,7 @@ def visualize_heat_extraction_rate_profiles( _format_axes(ax2) # Plot curves for requested boreholes - for i, borehole in enumerate(self.solver.boreholes): + for i, borehole in enumerate(self.solver.borefield): if i in iBoreholes: # Draw heat extraction rate profile line = ax2.plot( @@ -705,11 +707,11 @@ def visualize_temperatures( # If iBoreholes is None, then plot all boreholes if iBoreholes is None: - iBoreholes = range(len(self.solver.boreholes)) + iBoreholes = range(len(self.solver.borefield)) # Import temperatures T_b = self._temperatures(iBoreholes) # Borefield characteristic time - ts = np.mean([b.H for b in self.solver.boreholes])**2/(9.*self.alpha) + ts = self.borefield.H.mean()**2 / (9 * self.alpha) # Dimensionless time (log) lntts = np.log(self.time/ts) @@ -727,7 +729,7 @@ def visualize_temperatures( ax2.set_ylabel(r'$\bar{T}_b$') _format_axes(ax2) # Plot curves for requested boreholes - for i, borehole in enumerate(self.solver.boreholes): + for i, borehole in enumerate(self.solver.borefield): if i in iBoreholes: # Draw borehole wall temperature line = ax2.plot(lntts, T_b[iBoreholes.index(i)]) @@ -781,7 +783,7 @@ def visualize_temperatures( ax2.set_ylabel(r'$\bar{T}_b$') _format_axes(ax2) # Plot curves for requested boreholes - for i, borehole in enumerate(self.solver.boreholes): + for i, borehole in enumerate(self.solver.borefield): if i in iBoreholes: # Draw borehole wall temperature line = ax2.plot(lntts, T_b[iBoreholes.index(i)][n]) @@ -852,7 +854,7 @@ def visualize_temperature_profiles( # If iBoreholes is None, then plot all boreholes if iBoreholes is None: - iBoreholes = range(len(self.boreholes)) + iBoreholes = range(len(self.borefield)) # Import temperature profiles z, T_b = self._temperature_profiles(time, iBoreholes) @@ -871,7 +873,7 @@ def visualize_temperature_profiles( _format_axes(ax2) # Plot curves for requested boreholes - for i, borehole in enumerate(self.solver.boreholes): + for i, borehole in enumerate(self.solver.borefield): if i in iBoreholes: # Draw borehole wall temperature profile line = ax2.plot( @@ -928,7 +930,7 @@ def visualize_temperature_profiles( _format_axes(ax2) # Plot curves for requested boreholes - for i, borehole in enumerate(self.solver.boreholes): + for i, borehole in enumerate(self.solver.borefield): if i in iBoreholes: # Draw borehole wall temperature profile line = ax2.plot( @@ -994,7 +996,13 @@ def _heat_extraction_rates(self, iBoreholes): # heat extraction rate. i0 = self.solver._i0Segments[i] i1 = self.solver._i1Segments[i] - segment_ratios = self.solver.segment_ratios[i] + segment_ratios = self.solver.segment_ratios + if segment_ratios is None: + segment_ratios = np.full(i1 - i0, 1. / (i1 - i0)) + if isinstance(segment_ratios, list): + segment_ratios = segment_ratios[i] + if callable(segment_ratios): + segment_ratios = self.solver.segment_ratios(i1 - i0) if self.solver.nMassFlow == 0: Q_t.append( np.sum( @@ -1033,6 +1041,9 @@ def _heat_extraction_rate_profiles(self, time, iBoreholes): # Initialize lists z = [] Q_b = [] + nBoreSegments = np.broadcast_to( + self.solver.nSegments, + len(self.solver.borefield)) for i in iBoreholes: if self.boundary_condition == 'UHTR': # For the UHTR boundary condition, the solver only returns one @@ -1040,8 +1051,8 @@ def _heat_extraction_rate_profiles(self, time, iBoreholes): # The heat extraction rate is duplicated to draw from # z = D to z = D + H. z.append( - np.array([self.solver.boreholes[i].D, - self.solver.boreholes[i].D + self.solver.boreholes[i].H])) + np.array([self.solver.borefield.D[i], + self.solver.borefield.D[i] + self.solver.borefield.H[i]])) Q_b.append(np.array(2*[self.solver.Q_b])) else: i0 = self.solver._i0Segments[i] @@ -1069,21 +1080,25 @@ def _heat_extraction_rate_profiles(self, time, iBoreholes): kind='linear', copy=False, axis=2)(time) - if self.solver.nBoreSegments[i] > 1: + if nBoreSegments[i] > 1: # Borehole length ratio at the mid-depth of each segment - segment_ratios = self.solver.segment_ratios[i] + segment_ratios = self.solver.segment_ratios + if segment_ratios is None: + segment_ratios = np.full(i1 - i0, 1. / (i1 - i0)) + if isinstance(segment_ratios, list): + segment_ratios = segment_ratios[i] z.append( - self.solver.boreholes[i].D \ - + self.solver.boreholes[i]._segment_midpoints( - self.solver.nBoreSegments[i], + self.solver.borefield.D[i] \ + + self.solver.borefield[i]._segment_midpoints( + nBoreSegments[i], segment_ratios=segment_ratios)) Q_b.append(Q_bi) else: # If there is only one segment, the heat extraction rate is # duplicated to draw from z = D to z = D + H. z.append( - np.array([self.solver.boreholes[i].D, - self.solver.boreholes[i].D + self.solver.boreholes[i].H])) + np.array([self.solver.borefield.D[i], + self.solver.borefield.D[i] + self.solver.borefield.H[i]])) if self.solver.nMassFlow == 0: Q_b.append(np.repeat(Q_bi, 2, axis=0)) else: @@ -1118,7 +1133,11 @@ def _temperatures(self, iBoreholes): # borehole wall temperature. i0 = self.solver._i0Segments[i] i1 = self.solver._i1Segments[i] - segment_ratios = self.solver.segment_ratios[i] + segment_ratios = self.solver.segment_ratios + if segment_ratios is None: + segment_ratios = np.full(i1 - i0, 1. / (i1 - i0)) + if isinstance(segment_ratios, list): + segment_ratios = segment_ratios[i] if self.solver.nMassFlow == 0: T_b.append( np.sum( @@ -1156,6 +1175,9 @@ def _temperature_profiles(self, time, iBoreholes): # Initialize lists z = [] T_b = [] + nBoreSegments = np.broadcast_to( + self.solver.nSegments, + len(self.solver.borefield)) for i in iBoreholes: if self.boundary_condition == 'UBWT': # For the UBWT boundary condition, the solver only returns one @@ -1163,8 +1185,8 @@ def _temperature_profiles(self, time, iBoreholes): # boreholes). The temperature is duplicated to draw from # z = D to z = D + H. z.append( - np.array([self.solver.boreholes[i].D, - self.solver.boreholes[i].D + self.solver.boreholes[i].H])) + np.array([self.solver.borefield.D[i], + self.solver.borefield.D[i] + self.solver.borefield.H[i]])) if time is None: # If time is None, temperatures are extracted at the last # time step. @@ -1208,41 +1230,46 @@ def _temperature_profiles(self, time, iBoreholes): kind='linear', copy=False, axis=2)(time) - if self.solver.nBoreSegments[i] > 1: + if nBoreSegments[i] > 1: # Borehole length ratio at the mid-depth of each segment - segment_ratios = self.solver.segment_ratios[i] + segment_ratios = self.solver.segment_ratios + if segment_ratios is None: + segment_ratios = np.full(i1 - i0, 1. / (i1 - i0)) + if isinstance(segment_ratios, list): + segment_ratios = segment_ratios[i] z.append( - self.solver.boreholes[i].D \ - + self.solver.boreholes[i]._segment_midpoints( - self.solver.nBoreSegments[i], + self.solver.borefield.D[i] \ + + self.solver.borefield[i]._segment_midpoints( + nBoreSegments[i], segment_ratios=segment_ratios)) T_b.append(T_bi) else: # If there is only one segment, the temperature is # duplicated to draw from z = D to z = D + H. z.append( - np.array([self.solver.boreholes[i].D, - self.solver.boreholes[i].D + self.solver.boreholes[i].H])) + np.array([self.solver.borefield.D[i], + self.solver.borefield.D[i] + self.solver.borefield.H[i]])) if self.solver.nMassFlow == 0: T_b.append(np.repeat(T_bi, 2, axis=0)) else: T_b.append(np.repeat(T_bi, 2, axis=1)) return z, T_b - def _format_inputs(self, boreholes_or_network): + def _format_inputs(self, borefield_or_network): """ Process and format the inputs to the gFunction class. """ # Convert borehole to a list if a single borehole is provided - if isinstance(boreholes_or_network, Borehole): - boreholes_or_network = [boreholes_or_network] + if isinstance(borefield_or_network, (Borehole, list)): + borefield_or_network = Borefield.from_boreholes( + borefield_or_network) # Check if a borefield or a network is provided as an input and - # correctly assign the variables self.boreholes and self.network - if isinstance(boreholes_or_network, Network): - self.network = boreholes_or_network - self.boreholes = boreholes_or_network.b + # correctly assign the variables self.borefield and self.network + if isinstance(borefield_or_network, Network): + self.network = borefield_or_network + self.borefield = borefield_or_network.b # If a network is provided and no boundary condition is provided, # use 'MIFT' if self.boundary_condition is None: @@ -1258,7 +1285,7 @@ def _format_inputs(self, boreholes_or_network): self.cp_f = self.network.cp_f else: self.network = None - self.boreholes = boreholes_or_network + self.borefield = borefield_or_network # If a borefield is provided and no boundary condition is provided, # use 'UBWT' if self.boundary_condition is None: @@ -1297,14 +1324,11 @@ def _check_inputs(self): 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 not find_duplicates(self.boreholes), \ + assert isinstance(self.borefield, Borefield), \ + "The borefield is not a valid Borefield object." + assert len(self.borefield) > 0, \ + "The borefield is empty." + assert not find_duplicates(self.borefield), \ "There are duplicate boreholes in the borefield." assert (self.network is None and not self.boundary_condition=='MIFT') or isinstance(self.network, Network), \ "The network is not a valid 'Network' object." @@ -1425,7 +1449,7 @@ def uniform_heat_extraction(boreholes, time, alpha, use_similarities=True, def uniform_temperature(boreholes, time, alpha, nSegments=8, - segment_ratios=segment_ratios, kind='linear', + segment_ratios=utilities.segment_ratios, kind='linear', use_similarities=True, disTol=0.01, tol=1.0e-6, dtype=np.double, disp=False, **kwargs): """ @@ -1535,7 +1559,7 @@ def uniform_temperature(boreholes, time, alpha, nSegments=8, def equal_inlet_temperature( boreholes, UTubes, m_flow_borehole, cp_f, time, alpha, - kind='linear', nSegments=8, segment_ratios=segment_ratios, + kind='linear', nSegments=8, segment_ratios=utilities.segment_ratios, use_similarities=True, disTol=0.01, tol=1.0e-6, dtype=np.double, disp=False, **kwargs): """ @@ -1644,7 +1668,7 @@ def equal_inlet_temperature( def mixed_inlet_temperature( network, m_flow_network, cp_f, time, alpha, kind='linear', - nSegments=8, segment_ratios=segment_ratios, + nSegments=8, segment_ratios=utilities.segment_ratios, use_similarities=True, disTol=0.01, tol=1.0e-6, dtype=np.double, disp=False, **kwargs): """ diff --git a/pygfunction/networks.py b/pygfunction/networks.py index 9fa5884c..8f0322c7 100644 --- a/pygfunction/networks.py +++ b/pygfunction/networks.py @@ -1136,6 +1136,9 @@ def _format_inputs(self, m_flow_network, cp_f, nSegments, segment_ratios): self._segment_ratios = [segment_ratios] * self.nBoreholes elif isinstance(segment_ratios, list): self._segment_ratios = segment_ratios + elif callable(segment_ratios): + self._segment_ratios = [ + segment_ratios(nSegments_i) for nSegments_i in self.nSegments] else: raise ValueError( 'Incorrect format of the segment ratios list.') @@ -1497,6 +1500,9 @@ def _format_inputs(self, m_flow_network, cp_f, nSegments, segment_ratios): self._segment_ratios = [segment_ratios] * self.nBoreholes elif isinstance(segment_ratios, list): self._segment_ratios = segment_ratios + elif callable(segment_ratios): + self._segment_ratios = [ + segment_ratios(nSegments_i) for nSegments_i in self.nSegments] else: raise ValueError( 'Incorrect format of the segment ratios list.') diff --git a/pygfunction/solvers/_base_solver.py b/pygfunction/solvers/_base_solver.py index e73eb2f9..53c6ff04 100644 --- a/pygfunction/solvers/_base_solver.py +++ b/pygfunction/solvers/_base_solver.py @@ -1,19 +1,22 @@ # -*- coding: utf-8 -*- +from abc import ABC, abstractmethod +from collections.abc import Callable +from typing import Union, List from time import perf_counter import numpy as np +import numpy.typing as npt 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 +from .. import utilities -class _BaseSolver: +class _BaseSolver(ABC): """ Template for solver classes. @@ -21,10 +24,10 @@ class _BaseSolver: Attributes ---------- - boreholes : list of Borehole objects - List of boreholes included in the bore field. + borefield : Borefield object + The bore field. network : network object - Model of the network. + The network. time : float or array Values of time (in seconds) for which the g-function is evaluated. boundary_condition : str @@ -111,84 +114,84 @@ class _BaseSolver: 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 + def __init__( + self, + borefield: Borefield, + network: Union[Network, None], + time: npt.ArrayLike, + boundary_condition: str, + m_flow_borehole: Union[npt.ArrayLike, None] = None, + m_flow_network: Union[npt.ArrayLike, None] = None, + cp_f: Union[float, None] = None, + nSegments: int = 8, + segment_ratios: Union[npt.ArrayLike, List[npt.ArrayLike], Callable[[int], npt.ArrayLike]] = utilities.segment_ratios, + approximate_FLS: bool = False, + mQuad: int = 11, + nFLS: int = 10, + linear_threshold: Union[float, None] = None, + disp: bool = False, + profiles: bool = False, + kind: str = 'linear', + dtype: npt.DTypeLike = np.double, + **other_options): + # Input attributes + self.borefield = borefield 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.time = np.asarray(time) 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.nSegments = nSegments + self.segment_ratios = segment_ratios self.approximate_FLS = approximate_FLS self.mQuad = mQuad self.nFLS = nFLS + self.linear_threshold = linear_threshold 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 + self.nMassFlow = 0 + 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 + 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 - def initialize(self, *kwargs): + @property + def segment_lengths(self) -> np.ndarray: + """ + 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). + + """ + return self.segments.H + + @abstractmethod + def initialize(self, *kwargs) -> int: """ 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 @@ -197,14 +200,9 @@ def initialize(self, *kwargs): 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): + ... + + def solve(self, time: npt.ArrayLike, alpha: float) -> np.ndarray: """ Build and solve the system of equations. @@ -222,26 +220,28 @@ def solve(self, time, alpha): """ # Number of time values + scalar_time = isinstance(time, (float, np.floating, int, np.integer)) + if isinstance(time, list) or scalar_time: + time = np.asarray(time) self.time = time - nt = len(self.time) + nTimes = 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) + time_threshold = np.max(self.borefield.r_b)**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') + p0 = np.maximum(0, p_long - 1, dtype=int) if p_long > 0: time_long = np.concatenate([[time_threshold], self.time[p_long:]]) else: time_long = self.time - nt_long = len(time_long) + nTimes_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_b = self.segment_lengths H_tot = np.sum(H_b) if self.disp: print('Building and solving the system of equations ...', end='') @@ -250,56 +250,41 @@ def solve(self, time, alpha): if self.boundary_condition == 'UHTR': # Initialize g-function - gFunc = np.zeros(nt) + gFunc = np.zeros(nTimes) # Initialize segment heat extraction rates - Q_b = 1 + Q_b = np.broadcast_to(1., (self.nSources, nTimes)) # 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 + T_b = np.zeros((self.nSources, nTimes), dtype=self.dtype) + # Evaluate the g-function with uniform heat extraction along + # boreholes + T_b[:, p0:] = np.sum(h_ij.y[:, :, 1:], axis=1) + gFunc[p0:] = (T_b[:, p0:].T @ H_b).T / 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 + 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) + gFunc = np.zeros(nTimes) # Initialize segment heat extraction rates - Q_b = np.zeros((self.nSources, nt), dtype=self.dtype) - T_b = np.zeros(nt, dtype=self.dtype) + Q_b = np.zeros((self.nSources, nTimes), dtype=self.dtype) + T_b = np.zeros(nTimes, 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] + dt = np.concatenate( + [time_long[0:1], time_long[1:] - time_long[:-1]]) + for p in range(nTimes_long): # Thermal response factors evaluated at t=dt - h_dt = h_ij(dt) + h_dt = h_ij(dt[p]) # Reconstructed load history Q_reconstructed = self.load_history_reconstruction( - time_long[0:p+1], Q_b[:,p0:p+p0+1]) + 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) + h_ij.y[:, :, 1:], Q_reconstructed) # Evaluate the g-function with uniform borehole wall # temperature @@ -321,43 +306,45 @@ def solve(self, time, alpha): # Solve the system of equations X = np.linalg.solve(A, B) # Store calculated heat extraction rates - Q_b[:,p+p0] = X[0:self.nSources] + 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] + 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 + 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 + + # Broadcast T_b to expected size + T_b = np.broadcast_to(T_b, (self.nSources, nTimes)) elif self.boundary_condition == 'MIFT': if self.nMassFlow == 0: # Initialize g-function - gFunc = np.zeros((1, 1, nt)) + gFunc = np.zeros((1, 1, nTimes)) # 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) + Q_b = np.zeros((1, self.nSources, nTimes), dtype=self.dtype) + T_b = np.zeros((1, self.nSources, nTimes), dtype=self.dtype) else: # Initialize g-function - gFunc = np.zeros((self.nMassFlow, self.nMassFlow, nt)) + gFunc = np.zeros((self.nMassFlow, self.nMassFlow, nTimes)) # Initialize segment heat extraction rates Q_b = np.zeros( - (self.nMassFlow, self.nSources, nt), dtype=self.dtype) + (self.nMassFlow, self.nSources, nTimes), dtype=self.dtype) T_b = np.zeros( - (self.nMassFlow, self.nSources, nt), dtype=self.dtype) + (self.nMassFlow, self.nSources, nTimes), 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, + self.nSegments, segment_ratios=self.segment_ratios) k_s = self.network.p[0].k_s - for p in range(nt_long): + for p in range(nTimes_long): # Current thermal response factor matrix if p > 0: dt = time_long[p] - time_long[p-1] @@ -393,8 +380,8 @@ def solve(self, time, alpha): -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)], + a_b_j / (2. * np.pi * k_s * np.atleast_2d(self.segments.H).T), + a_in_j / (2. * np.pi * k_s * np.atleast_2d(self.segments.H).T)], [H_b, np.zeros(self.nSources + 1, dtype=self.dtype)]]) B = np.hstack( (-T_b0, @@ -403,8 +390,8 @@ def solve(self, time, alpha): # 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] + 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 @@ -413,13 +400,13 @@ def solve(self, time, alpha): 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) + 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 + 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)): @@ -428,9 +415,9 @@ def solve(self, time, alpha): a_in, a_b = self.network.coefficients_network_heat_extraction_rate( self.m_flow[i], self.cp_f, - self.nBoreSegments, + self.nSegments, segment_ratios=self.segment_ratios) - T_f_in = (-2 * np.pi * k_s * H_tot - a_b @ T_b[j,:,p0:]) / a_in + 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 @@ -440,17 +427,17 @@ def solve(self, time, alpha): 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 + 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 + 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,:,:] + gFunc = gFunc[0, 0, :] + Q_b = Q_b[0, :, :] + T_b = T_b[0, :, :] # Store temperature and heat extraction rate profiles if self.profiles: @@ -460,76 +447,9 @@ def solve(self, time, alpha): 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): + @staticmethod + def load_history_reconstruction( + time: np.ndarray, Q_b: np.ndarray) -> np.ndarray: """ Reconstructs the load history. @@ -552,7 +472,7 @@ def load_history_reconstruction(self, time, Q_b): # Number of heat sources nSources = Q_b.shape[0] # Time step sizes - dt = np.hstack((time[0], time[1:]-time[:-1])) + dt = np.hstack((time[0], time[1:] - time[:-1])) # Time vector t = np.hstack((0., time, time[-1] + time[0])) # Inverted time step sizes @@ -561,30 +481,61 @@ def load_history_reconstruction(self, time, Q_b): t_reconstructed = np.hstack((0., np.cumsum(dt_reconstructed))) # Accumulated heat extracted f = np.hstack( - (np.zeros((nSources, 1), dtype=self.dtype), + (np.zeros((nSources, 1)), np.cumsum(Q_b*dt, axis=1))) - f = np.hstack((f, f[:,-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 + Q_reconstructed = ( + sf(t_reconstructed[1:]) - sf(t_reconstructed[:-1]) + ) / dt_reconstructed return Q_reconstructed + @staticmethod + def temporal_superposition( + h_ij: np.ndarray, Q_reconstructed: np.ndarray) -> np.ndarray: + """ + 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 time steps + nTimes = 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[:,:,:nTimes], dQ) + + return T_b0 + 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 isinstance(self.borefield, Borefield), \ + "The borefield is not a valid 'Borefield' object." + assert len(self.borefield) > 0, \ + "The number of boreholes must be 1 or greater." assert self.network is None or isinstance(self.network, Network), \ "The network is not a valid 'Network' object." if self.boundary_condition == 'MIFT': @@ -597,49 +548,42 @@ def _check_inputs(self): 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), \ + assert not (isinstance(self.m_flow_borehole, 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, \ + assert isinstance(self.time, np.ndarray), \ + "Time should be an array." + assert (isinstance(self.nSegments, (int, np.integer)) + and self.nSegments >= 1) \ + or (isinstance(self.nSegments, (list, np.array)) + and len(self.nSegments) == len(self.borefield) + and np.min(self.nSegments) >=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, \ + assert (isinstance(self.boundary_condition, 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, \ + assert isinstance(self.approximate_FLS, 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, \ + assert isinstance(self.nFLS, int) and 1 <= self.nFLS <= 25, \ + "The option 'nFLS' should be a positive int and lower or equal " \ + "to 25." + assert isinstance(self.disp, bool), \ "The option 'disp' should be set to True or False." - assert type(self.profiles) is bool, \ + assert isinstance(self.profiles, bool), \ "The option 'profiles' should be set to True or False." - assert type(self.kind) is str, \ + assert isinstance(self.kind, 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 index eb66e4b3..66240c97 100644 --- a/pygfunction/solvers/detailed.py +++ b/pygfunction/solvers/detailed.py @@ -2,14 +2,11 @@ from time import perf_counter import numpy as np +import numpy.typing as npt 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 - ) +from ..heat_transfer import finite_line_source class Detailed(_BaseSolver): @@ -23,10 +20,10 @@ class Detailed(_BaseSolver): Parameters ---------- - boreholes : list of Borehole objects - List of boreholes included in the bore field. + borefield : Borefield object + The bore field. network : network object - Model of the network. + The network. time : float or array Values of time (in seconds) for which the g-function is evaluated. boundary_condition : str @@ -139,7 +136,7 @@ class Detailed(_BaseSolver): rate and reversible flow direction. Renewable Energy, 228, 120661. """ - def initialize(self, **kwargs): + def initialize(self, **kwargs) -> int: """ Split boreholes into segments. @@ -152,11 +149,20 @@ def initialize(self, **kwargs): """ # Split boreholes into segments - self.boreSegments = self.borehole_segments() - nSources = len(self.boreSegments) - return nSources + self.nBoreholes = len(self.borefield) + self.segments = self.borefield.segments( + self.nSegments, self.segment_ratios) + self._i1Segments = np.cumsum( + np.broadcast_to(self.nSegments, self.nBoreholes), + dtype=int) + self._i0Segments = np.concatenate( + ([0], self._i1Segments[:-1]), + dtype=int) + return len(self.segments) - def thermal_response_factors(self, time, alpha, kind='linear'): + def thermal_response_factors( + self, time: npt.ArrayLike, alpha: float, kind: str = 'linear' + ) -> interp1d: """ Evaluate the segment-to-segment thermal response factors for all pairs of segments in the borefield at all time steps using the finite line @@ -196,36 +202,38 @@ def thermal_response_factors(self, time, alpha, kind='linear'): 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() + nBoreholes = len(self.borefield) + segment_lengths = self.segment_lengths # --------------------------------------------------------------------- # Segment-to-segment thermal response factors for same-borehole # thermal interactions # --------------------------------------------------------------------- - h, i_segment, j_segment = \ + h, i, j = \ self._thermal_response_factors_borehole_to_self(time, alpha) # Broadcast values to h_ij matrix - h_ij[j_segment, i_segment, 1:] = h + h_ij[i, j, 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)): + i1 = self._i1Segments + i0 = self._i0Segments + for i, (_i0, _i1) in enumerate(zip(i0, i1)): # Segments of the receiving borehole - b2 = self.boreSegments[i0:i1] + b2 = self.segments[_i0:_i1] if i+1 < nBoreholes: # Segments of the emitting borehole - b1 = self.boreSegments[i1:] + b1 = self.segments[_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:] = \ + 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] + segment_lengths[_i0:_i1], + segment_lengths[_i1:]).T[:, :, np.newaxis] # Return 2d array if time is a scalar if np.isscalar(time): @@ -239,7 +247,8 @@ def thermal_response_factors(self, time, alpha, kind='linear'): return h_ij - def _thermal_response_factors_borehole_to_self(self, time, alpha): + def _thermal_response_factors_borehole_to_self( + self, time: npt.ArrayLike, alpha: float) -> np.ndarray: """ Evaluate the segment-to-segment thermal response factors for all pairs of segments between each borehole and itself. @@ -261,40 +270,20 @@ def _thermal_response_factors_borehole_to_self(self, time, alpha): 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) + i1 = self._i1Segments + i0 = self._i0Segments + i, j = zip(*[ + np.meshgrid( + np.arange(_j0, _j1, dtype=int), + np.arange(_j0, _j1, dtype=int), + indexing='ij') + for _j0, _j1 in zip(i0, i1) ]) - 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 + i = np.concatenate([_i.flatten() for _i in i]) + j = np.concatenate([_j.flatten() for _j in j]) + segments_i = self.segments[i] + segments_j = self.segments[j] + h = finite_line_source( + time, alpha, segments_j, segments_i, outer=False, + approximation=self.approximate_FLS, M=self.mQuad, N=self.nFLS) + return h, i, j diff --git a/pygfunction/solvers/equivalent.py b/pygfunction/solvers/equivalent.py index 0a4129a3..01470940 100644 --- a/pygfunction/solvers/equivalent.py +++ b/pygfunction/solvers/equivalent.py @@ -1,7 +1,9 @@ # -*- coding: utf-8 -*- from time import perf_counter +from typing import Tuple, List import numpy as np +import numpy.typing as npt from scipy.cluster.hierarchy import ( cut_tree, dendrogram, @@ -10,11 +12,11 @@ from scipy.interpolate import interp1d as interp1d from ._base_solver import _BaseSolver -from ..boreholes import _EquivalentBorehole +from ..borefield import Borefield, _EquivalentBorefield +from ..boreholes import Borehole, _EquivalentBorehole from ..heat_transfer import ( finite_line_source, - finite_line_source_equivalent_boreholes_vectorized, - finite_line_source_vectorized + finite_line_source_vertical ) from ..networks import _EquivalentNetwork @@ -34,8 +36,8 @@ class Equivalent(_BaseSolver): Parameters ---------- - boreholes : list of Borehole objects - List of boreholes included in the bore field. + borefield : Borefield object + The bore field. network : network object Model of the network. time : float or array @@ -172,7 +174,9 @@ class Equivalent(_BaseSolver): rate and reversible flow direction. Renewable Energy, 228, 120661. """ - def initialize(self, disTol=0.01, tol=1.0e-6, kClusters=1, **kwargs): + def initialize( + self, disTol: float = 0.01, tol: float = 1.0e-6, + kClusters: int = 1, **kwargs) -> int: """ Initialize paramteters. Identify groups for equivalent boreholes. @@ -191,16 +195,22 @@ def initialize(self, disTol=0.01, tol=1.0e-6, kClusters=1, **kwargs): 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)] + # Split boreholes into segments + self.segments = self.borefield.segments( + self.nSegments, self.segment_ratios) + nEqBoreholes = self.nEqBoreholes + nSegments = np.broadcast_to(self.nSegments, nEqBoreholes) + self._i1Segments = np.cumsum( + nSegments, + dtype=int) + self._i0Segments = np.concatenate( + ([0], self._i1Segments[:-1]), + dtype=int) return nSources - def thermal_response_factors(self, time, alpha, kind='linear'): + def thermal_response_factors( + self, time: npt.ArrayLike, alpha: float, kind: str = 'linear' + ) -> interp1d: """ Evaluate the segment-to-segment thermal response factors for all pairs of segments in the borefield at all time steps using the finite line @@ -234,13 +244,17 @@ def thermal_response_factors(self, time, alpha, kind='linear'): if self.disp: print('Calculating segment to segment response factors ...', end='') + nEqBoreholes = self.nEqBoreholes + nSegments = np.broadcast_to(self.nSegments, nEqBoreholes) + i1 = self._i1Segments + i0 = self._i0Segments # 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_lengths = self.segment_lengths # --------------------------------------------------------------------- # Segment-to-segment thermal response factors for borehole-to-borehole @@ -253,23 +267,25 @@ def thermal_response_factors(self, time, alpha, kind='linear'): 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 + [self.borefield[j].nBoreholes for (i, j) in pairs]) # Evaluate FLS at all time steps - h = finite_line_source_equivalent_boreholes_vectorized( - time, alpha, dis, wDis, H1, D1, H2, D2, N2) + borefield_1 = Borefield(H1, D1, self.borefield.r_b[i], 0., 0.) + borefield_2 = Borefield(H2, D2, self.borefield.r_b[j], 0., 0.) + h = finite_line_source_vertical( + time, alpha, borefield_1, borefield_2, distances=dis, + weights=(wDis/N2).T, outer=False) # 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, :] + i_segment = i0[i] + i_pair + j_segment = i0[j] + j_pair + h_ij[j_segment, i_segment, 1:] = h[k_pair, k, :] 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 + h_ij[i_segment, j_segment, 1:] = ( + h[k_pair, k, :] + * segment_lengths[j_segment, np.newaxis] + / segment_lengths[i_segment, np.newaxis] + ) # --------------------------------------------------------------------- # Segment-to-segment thermal response factors for same-borehole thermal @@ -283,20 +299,18 @@ def thermal_response_factors(self, time, alpha, kind='linear'): 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, + dis = self.borefield.r_b[i] + borefield_1 = Borefield(H1, D1, dis, np.zeros_like(H1), 0.) + borefield_2 = Borefield(H2, D2, dis, np.zeros_like(H1), 0.) + h = finite_line_source( + time, alpha, borefield_1, borefield_2, outer=False, 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 + i_segment = i0[i] + i_pair + j_segment = i0[i] + j_pair h_ij[j_segment, i_segment, 1:] = \ - h_ij[j_segment, i_segment, 1:] + h[0, k_pair, :] + h_ij[j_segment, i_segment, 1:] + h[k_pair, :] # Return 2d array if time is a scalar if np.isscalar(time): @@ -310,7 +324,7 @@ def thermal_response_factors(self, time, alpha, kind='linear'): return h_ij - def find_groups(self, tol=1e-6): + def find_groups(self, tol: float = 1e-6) -> int: """ Identify groups of boreholes that can be represented by a single equivalent borehole for the calculation of the g-function. @@ -341,14 +355,18 @@ def find_groups(self, tol=1e-6): tic = perf_counter() # Temperature change of individual boreholes - self.nBoreholes = len(self.boreholes) + self.nBoreholes = len(self.borefield) # Equivalent field formed by all boreholes - eqField = _EquivalentBorehole(self.boreholes) + eqField = _EquivalentBorehole.from_borefield(self.borefield) 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) + data = np.sum( + finite_line_source( + np.inf, 1., self.borefield, self.borefield + ), axis=1 + ).reshape(-1,1) # Split boreholes into groups of same dimensions - unique_boreholes = self._find_unique_boreholes(self.boreholes) + unique_boreholes = self._find_unique_boreholes(self.borefield) # Initialize empty list of clusters self.clusters = [] self.nEqBoreholes = 0 @@ -391,15 +409,16 @@ def find_groups(self, tol=1e-6): 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]) + self.borefield = _EquivalentBorefield.from_equivalent_boreholes( + [_EquivalentBorehole.from_boreholes( + [borehole + for borehole, cluster in zip(self.borefield, self.clusters) + if cluster==i]) + for i in range(self.nEqBoreholes)]) + self.wBoreholes = np.array([b.nBoreholes for b in self.borefield]) # Find similar pairs of boreholes self.borehole_to_self, self.borehole_to_borehole = \ - self._find_axial_borehole_pairs(self.boreholes) + self._find_axial_borehole_pairs(self.borefield) # Store unique distances in the bore field self.dis = eqField.unique_distance(eqField, self.disTol)[0][1:] @@ -407,10 +426,10 @@ def find_groups(self, tol=1e-6): pipes = [self.network.p[self.clusters.index(i)] for i in range(self.nEqBoreholes)] self.network = _EquivalentNetwork( - self.boreholes, + self.borefield, pipes, - nSegments=self.nBoreSegments[0], - segment_ratios=self.segment_ratios[0]) + nSegments=self.nSegments, + segment_ratios=self.segment_ratios) # Stop chrono toc = perf_counter() @@ -419,9 +438,10 @@ def find_groups(self, tol=1e-6): print(f'Calculations will be done using {self.nEqBoreholes} ' f'equivalent boreholes') - return self.nBoreSegments[0]*self.nEqBoreholes + return self.nSegments * self.nEqBoreholes - def segment_lengths(self): + @property + def segment_lengths(self) -> np.ndarray: """ Return the length of all segments in the bore field. @@ -436,17 +456,16 @@ def segment_lengths(self): """ # Borehole lengths - H = np.array([seg.H*seg.nBoreholes - for (borehole, nSegments, ratios) in zip( - self.boreholes, - self.nBoreSegments, - self.segment_ratios) + H = np.array([seg.H * seg.nBoreholes + for borehole in self.borefield for seg in borehole.segments( - nSegments, segment_ratios=ratios)], + self.nSegments, + segment_ratios=self.segment_ratios)], dtype=self.dtype) return H - def _compare_boreholes(self, borehole1, borehole2): + def _compare_boreholes( + self, borehole1: Borehole, borehole2: Borehole) -> bool: """ Compare two boreholes and checks if they have the same dimensions : H, D, and r_b. @@ -473,7 +492,9 @@ def _compare_boreholes(self, borehole1, borehole2): similarity = False return similarity - def _compare_real_pairs(self, pair1, pair2): + def _compare_real_pairs( + self, pair1: Tuple[Borehole, Borehole], + pair2: Tuple[Borehole, Borehole]) -> bool: """ Compare two pairs of boreholes or segments and return True if the two pairs have the same FLS solution for real sources. @@ -510,7 +531,9 @@ def _compare_real_pairs(self, pair1, pair2): similarity = False return similarity - def _compare_image_pairs(self, pair1, pair2): + def _compare_image_pairs( + self, pair1: Tuple[Borehole, Borehole], + pair2: Tuple[Borehole, Borehole]) -> bool: """ Compare two pairs of boreholes or segments and return True if the two pairs have the same FLS solution for mirror sources. @@ -542,7 +565,9 @@ def _compare_image_pairs(self, pair1, pair2): similarity = False return similarity - def _compare_realandimage_pairs(self, pair1, pair2): + def _compare_realandimage_pairs( + self, pair1: Tuple[Borehole, Borehole], + pair2: Tuple[Borehole, Borehole]) -> bool: """ 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. @@ -567,15 +592,18 @@ def _compare_realandimage_pairs(self, pair1, pair2): similarity = False return similarity - def _find_axial_borehole_pairs(self, boreholes): + def _find_axial_borehole_pairs( + self, borefield: _EquivalentBorefield) -> Tuple[ + List[List[int]], List[List[Tuple[int, int]]] + ]: """ 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. + borefield : _EquivalentBorefield object + The equivalent bore field. Returns ------- @@ -590,18 +618,18 @@ def _find_axial_borehole_pairs(self, boreholes): # Compare for the full (real + image) FLS solution compare_pairs = self._compare_realandimage_pairs - nBoreholes = len(boreholes) + nBoreholes = len(borefield) 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): + for i, borehole_i in enumerate(borefield): # 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]): + if self._compare_boreholes(borehole_i, borefield[m]): borehole_set.append(i) break else: @@ -610,13 +638,13 @@ def _find_axial_borehole_pairs(self, boreholes): # 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): + for j, borehole_j in enumerate(borefield[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]) + pair_ref = (borefield[m], borefield[n]) # Add the pair (or the reciprocal pair) to a group # if a similar one is found if compare_pairs(pair0, pair_ref): @@ -634,14 +662,15 @@ def _find_axial_borehole_pairs(self, boreholes): borehole_to_borehole = [[(0, 0)]] return borehole_to_self, borehole_to_borehole - def _find_unique_boreholes(self, boreholes): + def _find_unique_boreholes( + self, borefield: Borefield) -> np.ndarray: """ Find unique sets of dimensions (h, D, r_b) in the bore field. Parameters ---------- - boreholes : list of Borehole objects - Boreholes in the bore field. + borefield : Borefield object + The bore field. Returns ------- @@ -651,9 +680,9 @@ def _find_unique_boreholes(self, boreholes): """ unique_boreholes = [] - for i, borehole_1 in enumerate(boreholes): + for i, borehole_1 in enumerate(borefield): for group in unique_boreholes: - borehole_2 = boreholes[group[0]] + borehole_2 = borefield[group[0]] # Add the borehole to a group if similar dimensions are found if self._compare_boreholes(borehole_1, borehole_2): group.append(i) @@ -662,11 +691,14 @@ def _find_unique_boreholes(self, boreholes): # If no similar boreholes are known, append the groups unique_boreholes.append([i]) - return unique_boreholes + return np.array(unique_boreholes, dtype=int) - def _find_unique_distances(self, dis, indices): + def _find_unique_distances( + self, dis: float, indices: List[Tuple[int, int]]) -> Tuple[ + np.ndarray, np.ndarray + ]: """ - Find the number of occurrences of each unique distances between pairs + Find the number of occurences of each unique distances between pairs of boreholes. Parameters @@ -685,10 +717,10 @@ def _find_unique_distances(self, dis, indices): pair of equivalent boreholes in indices. """ - wDis = np.zeros((len(dis), len(indices)), dtype=np.uint) + wDis = np.zeros((len(dis), len(indices)), dtype=int) for k, pair in enumerate(indices): i, j = pair - b1, b2 = self.boreholes[i], self.boreholes[j] + b1, b2 = self.borefield[i], self.borefield[j] # Generate a flattened array of distances between boreholes i and j if not i == j: dis_ij = b1.distance(b2).flatten() @@ -696,7 +728,7 @@ def _find_unique_distances(self, dis, indices): # 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) + wDis_ij = np.zeros(len(dis), dtype=int) # Get insert positions for the distances iDis = np.searchsorted(dis, dis_ij, side='left') # Find indexes where previous index is closer @@ -705,10 +737,14 @@ def _find_unique_distances(self, dis, indices): np.add.at(wDis_ij, iDis, 1) wDis[:,k] = wDis_ij - return dis.reshape((1, -1)), wDis + return dis, wDis - def _map_axial_segment_pairs(self, iBor, jBor, - reaSource=True, imgSource=True): + def _map_axial_segment_pairs( + self, iBor: int, jBor: int, reaSource: bool = True, + imgSource: bool = True) -> Tuple[ + np.ndarray, np.ndarray, np.ndarray, np.ndarray, + np.ndarray, np.ndarray, np.ndarray + ]: """ Find axial (i.e. disregarding the radial distance) similarities between segment pairs along two boreholes to simplify the evaluation of the @@ -726,7 +762,7 @@ def _map_axial_segment_pairs(self, iBor, jBor, Returns ------- - H1 : float + H1 : array Length of the emitting segments. D1 : array Array of buried depths of the emitting segments. @@ -734,18 +770,18 @@ def _map_axial_segment_pairs(self, iBor, jBor, Length of the receiving segments. D2 : array Array of buried depths of the receiving segments. - i_pair : list + i_pair : array of int Indices of the emitting segments along a borehole. - j_pair : list + j_pair : array of int Indices of the receiving segments along a borehole. - k_pair : list + k_pair : array of int 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] + borehole1 = self.borefield[iBor] + borehole2 = self.borefield[jBor] assert reaSource or imgSource, \ "At least one of reaSource and imgSource must be True." if reaSource and imgSource: @@ -759,9 +795,9 @@ def _map_axial_segment_pairs(self, iBor, jBor, compare_pairs = self._compare_image_pairs # Dive both boreholes into segments segments1 = borehole1.segments( - self.nBoreSegments[iBor], segment_ratios=self.segment_ratios[iBor]) + self.nSegments, segment_ratios=self.segment_ratios) segments2 = borehole2.segments( - self.nBoreSegments[jBor], segment_ratios=self.segment_ratios[jBor]) + self.nSegments, segment_ratios=self.segment_ratios) # Prepare lists of segment lengths H1 = [] H2 = [] @@ -769,13 +805,13 @@ def _map_axial_segment_pairs(self, iBor, jBor, 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]) + i_pair = np.repeat(np.arange(self.nSegments, dtype=int), + self.nSegments) + j_pair = np.tile(np.arange(self.nSegments, dtype=int), + self.nSegments) # Empty list of indices for unique pairs - k_pair = np.empty(self.nBoreSegments[iBor] * self.nBoreSegments[jBor], - dtype=np.uint) + k_pair = np.empty(self.nSegments * self.nSegments, + dtype=int) unique_pairs = [] nPairs = 0 @@ -817,11 +853,18 @@ def _check_solver_specific_inputs(self): "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]), \ + assert (isinstance(self.nSegments, (int, np.integer)) + and self.nSegments >= 1), \ + "The argument for number of segments `nSegments` should be " \ + "of type int." + assert (self.segment_ratios is None + or callable(self.segment_ratios) + or (isinstance(self.segment_ratios, np.ndarray) + and len(self.segment_ratios) == self.nSegments)), \ + "Solver 'equivalent' can only handle identical segment_ratios " \ + "for all boreholes. None or a single array of size " \ + "(nSegments,) must be provided." + assert not np.any([b.is_tilted for b in self.borefield]), \ "Solver 'equivalent' can only handle vertical boreholes." if self.boundary_condition == 'MIFT': assert np.all(np.array(self.network.c, dtype=int) == -1), \ diff --git a/pygfunction/solvers/similarities.py b/pygfunction/solvers/similarities.py index 904de2ca..36a44419 100644 --- a/pygfunction/solvers/similarities.py +++ b/pygfunction/solvers/similarities.py @@ -1,13 +1,17 @@ # -*- coding: utf-8 -*- from time import perf_counter +from typing import Tuple, List import numpy as np +import numpy.typing as npt from scipy.interpolate import interp1d as interp1d from ._base_solver import _BaseSolver +from ..boreholes import Borehole +from ..borefield import Borefield from ..heat_transfer import ( - finite_line_source_inclined_vectorized, - finite_line_source_vectorized + finite_line_source, + finite_line_source_vertical ) @@ -24,10 +28,10 @@ class Similarities(_BaseSolver): Parameters ---------- - boreholes : list of Borehole objects - List of boreholes included in the bore field. + borefield : Borefield object + The bore field. network : network object - Model of the network. + The network. time : float or array Values of time (in seconds) for which the g-function is evaluated. boundary_condition : str @@ -153,7 +157,8 @@ class Similarities(_BaseSolver): rate and reversible flow direction. Renewable Energy, 228, 120661. """ - def initialize(self, disTol=0.01, tol=1.0e-6, **kwargs): + def initialize( + self, disTol: float = 0.01, tol: float = 1e-6, **kwargs) -> int: """ Split boreholes into segments and identify similarities in the borefield. @@ -171,12 +176,55 @@ def initialize(self, disTol=0.01, tol=1.0e-6, **kwargs): # Check the validity of inputs self._check_solver_specific_inputs() # Split boreholes into segments - self.boreSegments = self.borehole_segments() + self.nBoreholes = len(self.borefield) + self.segments = self.borefield.segments( + self.nSegments, self.segment_ratios) + self._i1Segments = np.cumsum( + np.broadcast_to(self.nSegments, self.nBoreholes), + dtype=int) + self._i0Segments = np.concatenate( + ([0], self._i1Segments[:-1]), + dtype=int) + # Shortcut for segment_ratios comparisons + if self.segment_ratios is None: + self._equal_segment_ratios = True + self._uniform_segment_ratios = np.ones( + self.borefield.nBoreholes, dtype=bool) + elif isinstance(self.segment_ratios, np.ndarray): + self._equal_segment_ratios = True + self._uniform_segment_ratios = np.full( + self.nBoreholes, + np.allclose( + self.segment_ratios, + self.segment_ratios[:1], + rtol=1e-6), + dtype=bool) + else: + if callable(self.segment_ratios): + if isinstance(self.nSegments, int): + segment_ratios = [self.segment_ratios(self.nSegments)] + else: + segment_ratios = [self.segment_ratios(nSeg) for nSeg in self.nSegments] + else: + segment_ratios = self.segment_ratios + self._equal_segment_ratios = ( + isinstance(self.nSegments, int) + or np.all(np.equal(self.nSegments, self.nSegments[0])) + and np.all([np.allclose(segment_ratios, segment_ratios[0]) for ratios in segment_ratios]) + ) + # Boreholes with a uniform discretization + self._uniform_segment_ratios = np.array([ + np.allclose(ratios, + ratios[:1], + rtol=1e-6) + for ratios in segment_ratios], dtype=bool) # Initialize similarities self.find_similarities() - return len(self.boreSegments) + return len(self.segments) - def thermal_response_factors(self, time, alpha, kind='linear'): + def thermal_response_factors( + self, time: npt.ArrayLike, alpha: float, kind: str = 'linear' + ) -> interp1d: """ Evaluate the segment-to-segment thermal response factors for all pairs of segments in the borefield at all time steps using the finite line @@ -195,14 +243,14 @@ def thermal_response_factors(self, time, alpha, kind='linear'): Values of time (in seconds) for which the g-function is evaluated. alpha : float Soil thermal diffusivity (in m2/s). - kind : string, optional + kind : str, optional Interpolation method used for segment-to-segment thermal response factors. See documentation for scipy.interpolate.interp1d. Default is 'linear'. Returns ------- - h_ij : interp1d + h_ij : interp1d object interp1d object (scipy.interpolate) of the matrix of segment-to-segment thermal response factors. @@ -210,6 +258,8 @@ def thermal_response_factors(self, time, alpha, kind='linear'): if self.disp: print('Calculating segment to segment response factors ...', end='') + nBoreholes = self.borefield.nBoreholes + nSegments = np.broadcast_to(self.nSegments, nBoreholes) # Number of time values nt = len(np.atleast_1d(time)) # Initialize chrono @@ -222,21 +272,23 @@ def thermal_response_factors(self, time, alpha, kind='linear'): # 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, :] + if np.any(self.borefield.is_vertical): + 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, :] + if np.any(self.borefield.is_tilted): + 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) @@ -255,23 +307,20 @@ def thermal_response_factors(self, time, alpha, kind='linear'): 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, + borefield_1 = Borefield(H1, D1, self.borefield.r_b[i], 0., 0.) + borefield_2 = Borefield(H2, D2, self.borefield.r_b[j], 0., 0.) + h = finite_line_source_vertical( + time, alpha, borefield_1, borefield_2, outer=False, distances=distances, 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 + h_ij[j_segment, i_segment, 1:] = h[k_segment, l_segment, :] + if (self._compare_boreholes(self.borefield[j], self.borefield[i]) and + nSegments[i] == nSegments[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, :] + h_ij[i_segment, j_segment, 1:] = h[k_segment, l_segment, :] else: - h_ij[i_segment, j_segment, 1:] = (h * H2.T / H1.T)[l_segment, k_segment, :] + h_ij[i_segment, j_segment, 1:] = (h.T * H2 / H1).T[k_segment, l_segment, :] # --------------------------------------------------------------------- # Segment-to-segment thermal response factors for borehole-to-borehole # thermal interactions (inclined boreholes) @@ -298,7 +347,8 @@ def thermal_response_factors(self, time, alpha, kind='linear'): return h_ij def _thermal_response_factors_borehole_to_borehole_inclined( - self, time, alpha): + self, time: npt.ArrayLike, alpha: float) -> Tuple[ + np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """ Evaluate the segment-to-segment thermal response factors for all pairs of inclined segments. @@ -316,11 +366,11 @@ def _thermal_response_factors_borehole_to_borehole_inclined( Finite line source solution. hT : array Reciprocal finite line source solution. - i_segment : list + i_segment : array of int Indices of the emitting segments in the bore field. - j_segment : list + j_segment : array of int Indices of the receiving segments in the bore field. - k_segment : list + k_segment : array of int 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. """ @@ -337,9 +387,9 @@ def _thermal_response_factors_borehole_to_borehole_inclined( 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) + i_segment = np.array([], dtype=int) + j_segment = np.array([], dtype=int) + k_segment = np.array([], dtype=int) k0 = 0 for pairs in self.borehole_to_borehole_inclined: # Index of first borehole pair in group @@ -371,14 +421,19 @@ def _thermal_response_factors_borehole_to_borehole_inclined( 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) + borefield_1 = Borefield( + H1, D1, rb1, x1, y1, tilt=tilt1, orientation=orientation1) + borefield_2 = Borefield( + H2, D2, rb1, x2, y2, tilt=tilt2, orientation=orientation2) + h = finite_line_source( + time, alpha, borefield_1, borefield_2, outer=False, + approximation=self.approximate_FLS, M=self.mQuad, 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): + def _thermal_response_factors_borehole_to_self_inclined( + self, time: npt.ArrayLike, alpha: float) -> Tuple[ + np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """ Evaluate the segment-to-segment thermal response factors for all pairs of segments between each inclined borehole and itself. @@ -394,11 +449,11 @@ def _thermal_response_factors_borehole_to_self_inclined(self, time, alpha): ------- h : array Finite line source solution. - i_segment : list + i_segment : array of int Indices of the emitting segments in the bore field. - j_segment : list + j_segment : array of int Indices of the receiving segments in the bore field. - k_segment : list + k_segment : array of int 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. """ @@ -415,9 +470,9 @@ def _thermal_response_factors_borehole_to_self_inclined(self, time, alpha): 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) + i_segment = np.array([], dtype=int) + j_segment = np.array([], dtype=int) + k_segment = np.array([], dtype=int) k0 = 0 for group in self.borehole_to_self_inclined: # Index of first borehole in group @@ -450,13 +505,18 @@ def _thermal_response_factors_borehole_to_self_inclined(self, time, alpha): 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) + borefield_1 = Borefield( + H1, D1, rb1, x1, y1, tilt=tilt1, orientation=orientation1) + borefield_2 = Borefield( + H2, D2, rb1, x2, y2, tilt=tilt2, orientation=orientation2) + h = finite_line_source( + time, alpha, borefield_1, borefield_2, outer=False, + approximation=self.approximate_FLS, M=self.mQuad, N=self.nFLS) return h, i_segment, j_segment, k_segment - def _thermal_response_factors_borehole_to_self_vertical(self, time, alpha): + def _thermal_response_factors_borehole_to_self_vertical( + self, time: npt.ArrayLike, alpha: float) -> Tuple[ + np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """ Evaluate the segment-to-segment thermal response factors for all pairs of segments between each vertical borehole and itself. @@ -472,11 +532,11 @@ def _thermal_response_factors_borehole_to_self_vertical(self, time, alpha): ------- h : array Finite line source solution. - i_segment : list + i_segment : array of int Indices of the emitting segments in the bore field. - j_segment : list + j_segment : array of int Indices of the receiving segments in the bore field. - k_segment : list + k_segment : array of int 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. """ @@ -485,9 +545,9 @@ def _thermal_response_factors_borehole_to_self_vertical(self, time, alpha): 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) + i_segment = np.array([], dtype=int) + j_segment = np.array([], dtype=int) + k_segment = np.array([], dtype=int) k0 = 0 for group in self.borehole_to_self_vertical: # Index of first borehole in group @@ -505,17 +565,21 @@ def _thermal_response_factors_borehole_to_self_vertical(self, time, alpha): 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)) + dis = np.append(dis, np.full(len(H1_i), self.borefield.r_b[i])) else: - dis = self.boreholes[i].r_b + dis = self.borefield.r_b[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 += 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) + borefield_1 = Borefield( + H1, D1, dis, 0., 0.) + borefield_2 = Borefield( + H2, D2, dis, 0., 0.) + h = finite_line_source( + time, alpha, borefield_1, borefield_2, outer=False, + approximation=self.approximate_FLS, M=self.mQuad, N=self.nFLS) return h, i_segment, j_segment, k_segment def find_similarities(self): @@ -534,11 +598,11 @@ def find_similarities(self): # 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) + self._find_axial_borehole_pairs(self.borefield) # 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) + self.borefield, self.borehole_to_borehole_vertical) # Stop chrono toc = perf_counter() @@ -546,7 +610,8 @@ def find_similarities(self): return - def _compare_boreholes(self, borehole1, borehole2): + def _compare_boreholes( + self, borehole1: Borehole, borehole2: Borehole) -> bool: """ Compare two boreholes and checks if they have the same dimensions : H, D, r_b, and tilt. @@ -574,7 +639,9 @@ def _compare_boreholes(self, borehole1, borehole2): similarity = False return similarity - def _compare_real_pairs_vertical(self, pair1, pair2): + def _compare_real_pairs_vertical( + self, pair1: Tuple[Borehole, Borehole], + pair2: Tuple[Borehole, Borehole]) -> bool: """ Compare two pairs of vertical boreholes or segments and return True if the two pairs have the same FLS solution for real sources. @@ -611,7 +678,9 @@ def _compare_real_pairs_vertical(self, pair1, pair2): similarity = False return similarity - def _compare_image_pairs_vertical(self, pair1, pair2): + def _compare_image_pairs_vertical( + self, pair1: Tuple[Borehole, Borehole], + pair2: Tuple[Borehole, Borehole]) -> bool: """ Compare two pairs of vertical boreholes or segments and return True if the two pairs have the same FLS solution for mirror sources. @@ -643,7 +712,9 @@ def _compare_image_pairs_vertical(self, pair1, pair2): similarity = False return similarity - def _compare_realandimage_pairs_vertical(self, pair1, pair2): + def _compare_realandimage_pairs_vertical( + self, pair1: Tuple[Borehole, Borehole], + pair2: Tuple[Borehole, Borehole]) -> bool: """ 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 @@ -669,7 +740,9 @@ def _compare_realandimage_pairs_vertical(self, pair1, pair2): similarity = False return similarity - def _compare_real_pairs_inclined(self, pair1, pair2): + def _compare_real_pairs_inclined( + self, pair1: Tuple[Borehole, Borehole], + pair2: Tuple[Borehole, Borehole]) -> bool: """ Compare two pairs of inclined boreholes or segments and return True if the two pairs have the same FLS solution for real sources. @@ -717,7 +790,9 @@ def _compare_real_pairs_inclined(self, pair1, pair2): similarity = False return similarity - def _compare_image_pairs_inclined(self, pair1, pair2): + def _compare_image_pairs_inclined( + self, pair1: Tuple[Borehole, Borehole], + pair2: Tuple[Borehole, Borehole]) -> bool: """ Compare two pairs of inclined boreholes or segments and return True if the two pairs have the same FLS solution for mirror sources. @@ -765,7 +840,9 @@ def _compare_image_pairs_inclined(self, pair1, pair2): similarity = False return similarity - def _compare_realandimage_pairs_inclined(self, pair1, pair2): + def _compare_realandimage_pairs_inclined( + self, pair1: Tuple[Borehole, Borehole], + pair2: Tuple[Borehole, Borehole]) -> bool: """ 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 @@ -820,7 +897,11 @@ def _compare_realandimage_pairs_inclined(self, pair1, pair2): similarity = False return similarity - def _find_axial_borehole_pairs(self, boreholes): + def _find_axial_borehole_pairs( + self, borefield: Borefield) -> Tuple[ + List[List[int]], List[List[int]], List[List[Tuple[int, int]]], + List[List[Tuple[int, int]]] + ]: """ Find axial (i.e. disregarding the radial distance) similarities between borehole pairs to simplify the evaluation of the FLS solution. @@ -840,7 +921,14 @@ def _find_axial_borehole_pairs(self, boreholes): boreholes that share the same (pairwise) dimensions (H, D). """ - nBoreholes = len(boreholes) + nBoreholes = len(borefield) + nSegments = np.broadcast_to(self.nSegments, nBoreholes) + if self.segment_ratios is None: + segment_ratios = [None] * nBoreholes + elif isinstance(self.segment_ratios, np.ndarray): + segment_ratios = [self.segment_ratios] * nBoreholes + elif callable(self.segment_ratios): + segment_ratios = [self.segment_ratios(nSeg) for nSeg in nSegments] borehole_to_self_vertical = [] borehole_to_self_inclined = [] # Only check for similarities if there is more than one borehole @@ -848,9 +936,9 @@ def _find_axial_borehole_pairs(self, boreholes): 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)): + zip(borefield, nSegments, segment_ratios)): # Compare the borehole to all known unique sets of dimensions - if borehole_i.is_vertical(): + if borehole_i.is_vertical: borehole_to_self = borehole_to_self_vertical compare_pairs = self._compare_realandimage_pairs_vertical else: @@ -860,11 +948,12 @@ def _find_axial_borehole_pairs(self, boreholes): 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 + if (self._compare_boreholes(borehole_i, borefield[m]) and (self._equal_segment_ratios or - (nSegments_i == self.nBoreSegments[m] and + (nSegments_i == nSegments[m] and + segment_ratios[m] is None or np.allclose(ratios_i, - self.segment_ratios[m], + segment_ratios[m], rtol=self.tol)))): borehole_set.append(i) break @@ -873,14 +962,14 @@ def _find_axial_borehole_pairs(self, boreholes): 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:]), + zip(borefield[i+1:], + nSegments[i+1:], + 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(): + 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: @@ -888,30 +977,34 @@ def _find_axial_borehole_pairs(self, boreholes): compare_pairs = self._compare_realandimage_pairs_inclined for pairs in borehole_to_borehole: m, n = pairs[0] - pair_ref = (boreholes[m], boreholes[n]) + pair_ref = (borefield[m], borefield[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 + (nSegments_i == nSegments[m] and + nSegments_j == nSegments[n] and + segment_ratios[m] is None or np.allclose(ratios_i, - self.segment_ratios[m], + segment_ratios[m], rtol=self.tol) and + segment_ratios[n] is None or np.allclose(ratios_j, - self.segment_ratios[n], + 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 + (nSegments_j == nSegments[m] and + nSegments_i == nSegments[n] and + segment_ratios[m] is None or np.allclose(ratios_j, - self.segment_ratios[m], + segment_ratios[m], rtol=self.tol) and + segment_ratios[n] is None or np.allclose(ratios_i, - self.segment_ratios[n], + segment_ratios[n], rtol=self.tol)))): pairs.append((j, i)) break @@ -921,7 +1014,7 @@ def _find_axial_borehole_pairs(self, boreholes): else: # Outputs for a single borehole - if boreholes[0].is_vertical: + if borefield[0].is_vertical: borehole_to_self_vertical = [[0]] borehole_to_self_inclined = [] else: @@ -932,14 +1025,17 @@ def _find_axial_borehole_pairs(self, boreholes): 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): + def _find_distances( + self, borefield: Borefield, + borehole_to_borehole: List[List[Tuple[int, int]]] + ) -> Tuple[List[List[float]], List[np.ndarray]]: """ Find unique distances between pairs of boreholes for each unique pair of boreholes in the bore field. Parameters ---------- - boreholes : list of Borehole objects + borefield : Borefield object Boreholes in the bore field. borehole_to_borehole : list Lists of tuples of borehole indexes for each unique pair of @@ -957,7 +1053,7 @@ def _find_distances(self, boreholes, borehole_to_borehole): 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] + [np.empty(len(group), dtype=int) for group in borehole_to_borehole] # Find unique distances for each group for i, (pairs, distances, distance_indices) in enumerate( zip(borehole_to_borehole, @@ -966,7 +1062,7 @@ def _find_distances(self, boreholes, borehole_to_borehole): nPairs = len(pairs) # Array of all borehole-to-borehole distances within the group all_distances = np.array( - [boreholes[pair[0]].distance(boreholes[pair[1]]) + [borefield[pair[0]].distance(borefield[pair[1]]) for pair in pairs]) # Indices to sort the distance array i_sort = all_distances.argsort() @@ -1000,7 +1096,11 @@ def _find_distances(self, boreholes, borehole_to_borehole): return borehole_to_borehole_distances, borehole_to_borehole_indices def _map_axial_segment_pairs_vertical( - self, i, j, reaSource=True, imgSource=True): + self, i: int, j: int, reaSource: bool = True, + imgSource: bool = True) -> Tuple[ + np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, + np.ndarray, np.ndarray + ]: """ Find axial (i.e. disregarding the radial distance) similarities between segment pairs along two boreholes to simplify the evaluation of the @@ -1026,18 +1126,32 @@ def _map_axial_segment_pairs_vertical( Length of the receiving segments. D2 : array Array of buried depths of the receiving segments. - i_pair : list + i_pair : array of int Indices of the emitting segments along a borehole. - j_pair : list + j_pair : array of int Indices of the receiving segments along a borehole. - k_pair : list + k_pair : array of int Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions corresponding to all pairs in (i_pair, j_pair). """ + nBoreholes = self.borefield.nBoreholes + nSegments = np.broadcast_to(self.nSegments, nBoreholes) + if self.segment_ratios is None: + segment_ratios_i = None + segment_ratios_j = None + elif isinstance(self.segment_ratios, np.ndarray): + segment_ratios_i = self.segment_ratios + segment_ratios_j = self.segment_ratios + elif callable(self.segment_ratios): + segment_ratios_i = self.segment_ratios(nSegments[i]) + segment_ratios_j = self.segment_ratios(nSegments[j]) + else: + segment_ratios_i = self.segment_ratios[i] + segment_ratios_j = self.segment_ratios[j] # Initialize local variables - borehole1 = self.boreholes[i] - borehole2 = self.boreholes[j] + borehole1 = self.borefield[i] + borehole2 = self.borefield[j] assert reaSource or imgSource, \ "At least one of reaSource and imgSource must be True." if reaSource and imgSource: @@ -1051,9 +1165,9 @@ def _map_axial_segment_pairs_vertical( compare_pairs = self._compare_image_pairs_vertical # Dive both boreholes into segments segments1 = borehole1.segments( - self.nBoreSegments[i], segment_ratios=self.segment_ratios[i]) + nSegments[i], segment_ratios=segment_ratios_i) segments2 = borehole2.segments( - self.nBoreSegments[j], segment_ratios=self.segment_ratios[j]) + nSegments[j], segment_ratios=segment_ratios_j) # Prepare lists of segment lengths H1 = [] H2 = [] @@ -1061,13 +1175,13 @@ def _map_axial_segment_pairs_vertical( 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]) + i_pair = np.repeat(np.arange(nSegments[i], dtype=int), + nSegments[j]) + j_pair = np.tile(np.arange(nSegments[j], dtype=int), + nSegments[i]) # Empty list of indices for unique pairs - k_pair = np.empty(self.nBoreSegments[i] * self.nBoreSegments[j], - dtype=np.uint) + k_pair = np.empty(nSegments[i] * nSegments[j], + dtype=int) unique_pairs = [] nPairs = 0 @@ -1098,7 +1212,13 @@ def _map_axial_segment_pairs_vertical( 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): + self, i: int, j: int, reaSource: bool = True, + imgSource: bool = True) -> Tuple[ + np.ndarray, np.ndarray, np.ndarray, np.ndarray, + np.ndarray, np.ndarray, np.ndarray, np.ndarray, + np.ndarray, np.ndarray, np.ndarray, np.ndarray, + np.ndarray, np.ndarray, np.ndarray, np.ndarray + ]: """ Find axial similarities between segment pairs along two boreholes to simplify the evaluation of the FLS solution. @@ -1141,17 +1261,28 @@ def _map_axial_segment_pairs_inclined( 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 + i_pair : array of int Indices of the emitting segments along a borehole. - j_pair : list + j_pair : array of int Indices of the receiving segments along a borehole. - k_pair : list + k_pair : array of int Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions corresponding to all pairs in (i_pair, j_pair). """ + nBoreholes = self.nBoreholes + nSegments = np.broadcast_to(self.nSegments, nBoreholes) + if self.segment_ratios is None: + segment_ratios_i = None + segment_ratios_j = None + elif isinstance(self.segment_ratios, np.ndarray): + segment_ratios_i = self.segment_ratios + segment_ratios_j = self.segment_ratios + else: + segment_ratios_i = self.segment_ratios[i] + segment_ratios_j = self.segment_ratios[j] # Initialize local variables - borehole1 = self.boreholes[i] - borehole2 = self.boreholes[j] + borehole1 = self.borefield[i] + borehole2 = self.borefield[j] assert reaSource or imgSource, \ "At least one of reaSource and imgSource must be True." if reaSource and imgSource: @@ -1165,9 +1296,9 @@ def _map_axial_segment_pairs_inclined( compare_pairs = self._compare_image_pairs_inclined # Dive both boreholes into segments segments1 = borehole1.segments( - self.nBoreSegments[i], segment_ratios=self.segment_ratios[i]) + nSegments[i], segment_ratios=segment_ratios_i) segments2 = borehole2.segments( - self.nBoreSegments[j], segment_ratios=self.segment_ratios[j]) + nSegments[j], segment_ratios=segment_ratios_j) # Prepare lists of FLS-inclined arguments rb1 = [] x1 = [] @@ -1183,13 +1314,13 @@ def _map_axial_segment_pairs_inclined( 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]) + i_pair = np.repeat(np.arange(nSegments[i], dtype=int), + nSegments[j]) + j_pair = np.tile(np.arange(nSegments[j], dtype=int), + nSegments[i]) # Empty list of indices for unique pairs - k_pair = np.empty(self.nBoreSegments[i] * self.nBoreSegments[j], - dtype=np.uint) + k_pair = np.empty(nSegments[i] * nSegments[j], + dtype=int) unique_pairs = [] nPairs = 0 @@ -1232,8 +1363,11 @@ def _map_axial_segment_pairs_inclined( 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): + self, i_pair: npt.ArrayLike, j_pair: npt.ArrayLike, + k_pair: npt.ArrayLike, borehole_to_borehole: List[Tuple[int, int]], + borehole_to_borehole_indices: List[int]) -> Tuple[ + np.ndarray, np.ndarray, np.ndarray, np.ndarray + ]: """ Return the maping of the unique segment-to-segment thermal response factors (h) to the complete h_ij array of the borefield, such that: @@ -1245,11 +1379,11 @@ def _map_segment_pairs_vertical( Parameters ---------- - i_pair : list + i_pair : array of int Indices of the emitting segments. - j_pair : list + j_pair : array of int Indices of the receiving segments. - k_pair : list + k_pair : array of int 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 @@ -1259,29 +1393,40 @@ def _map_segment_pairs_vertical( Returns ------- - i_segment : list + i_segment : array of int Indices of the emitting segments in the bore field. - j_segment : list + j_segment : array of int Indices of the receiving segments in the bore field. - k_segment : list + k_segment : array of int 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 + l_segment : array of int Indices of unique distances for all pairs in (i_pair, j_pair) in the bore field. """ + nBoreholes = self.borefield.nBoreholes + nSegments = np.broadcast_to(self.nSegments, nBoreholes) + i1 = self._i1Segments + i0 = self._i0Segments i_segment = np.concatenate( - [i_pair + self._i0Segments[i] for (i, j) in borehole_to_borehole]) + [i_pair + i0[i] for (i, j) in borehole_to_borehole], + dtype=int) 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)) + [j_pair + i0[j] for (i, j) in borehole_to_borehole], + dtype=int) + 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]) + [np.repeat(i, len(k_pair)) for i in borehole_to_borehole_indices], + dtype=int) return i_segment, j_segment, k_segment, l_segment def _map_segment_pairs_inclined( - self, i_pair, j_pair, k_pair, borehole_to_borehole): + self, i_pair: npt.ArrayLike, j_pair: npt.ArrayLike, + k_pair: npt.ArrayLike, borehole_to_borehole: List[Tuple[int, int]] + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """ Return the maping of the unique segment-to-segment thermal response factors (h) to the complete h_ij array of the borefield, such that: @@ -1293,11 +1438,11 @@ def _map_segment_pairs_inclined( Parameters ---------- - i_pair : list + i_pair : array of int Indices of the emitting segments. - j_pair : list + j_pair : array of int Indices of the receiving segments. - k_pair : list + k_pair : array of int 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 @@ -1305,20 +1450,28 @@ def _map_segment_pairs_inclined( Returns ------- - i_segment : list + i_segment : array of int Indices of the emitting segments in the bore field. - j_segment : list + j_segment : array of int Indices of the receiving segments in the bore field. - k_segment : list + k_segment : array of int 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. """ + nBoreholes = self.borefield.nBoreholes + nSegments = np.broadcast_to(self.nSegments, nBoreholes) + i1 = self._i1Segments + i0 = self._i0Segments i_segment = np.concatenate( - [i_pair + self._i0Segments[i] for (i, j) in borehole_to_borehole]) + [i_pair + i0[i] for (i, j) in borehole_to_borehole], + dtype=int) 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)) + [j_pair + i0[j] for (i, j) in borehole_to_borehole], + dtype=int) + k_segment = np.tile( + k_pair, + len(borehole_to_borehole)) return i_segment, j_segment, k_segment def _check_solver_specific_inputs(self): diff --git a/pygfunction/utilities.py b/pygfunction/utilities.py index 953f8026..34f7cfe1 100644 --- a/pygfunction/utilities.py +++ b/pygfunction/utilities.py @@ -110,7 +110,7 @@ def segment_ratios(nSegments, end_length_ratio=0.02): def is_even(n): "Returns True if n is even." return not(n & 0x1) - assert nSegments >= 1 and isinstance(nSegments, int), \ + assert nSegments >= 1 and isinstance(nSegments, (int, np.integer)), \ "The number of segments `nSegments` should be greater or equal " \ "to 1 and of type int." assert nSegments <= 2 or 0. < end_length_ratio < 0.5 and \ From 2e13a5a71bc2c5118f43b8a26696eca3a6b7c4bc Mon Sep 17 00:00:00 2001 From: MassimoCimmino Date: Tue, 29 Apr 2025 05:38:55 -0400 Subject: [PATCH 4/7] Update doc --- doc/source/modules/heat_transfer.rst | 1 - pygfunction/heat_transfer/__init__.py | 6 ++++++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/doc/source/modules/heat_transfer.rst b/doc/source/modules/heat_transfer.rst index 2ffc7109..94232e96 100644 --- a/doc/source/modules/heat_transfer.rst +++ b/doc/source/modules/heat_transfer.rst @@ -6,5 +6,4 @@ Heat Transfer Module .. automodule:: pygfunction.heat_transfer :members: - :undoc-members: :show-inheritance: diff --git a/pygfunction/heat_transfer/__init__.py b/pygfunction/heat_transfer/__init__.py index f28e1831..8c980898 100644 --- a/pygfunction/heat_transfer/__init__.py +++ b/pygfunction/heat_transfer/__init__.py @@ -1,3 +1,9 @@ from .finite_line_source import finite_line_source from .finite_line_source_vertical import finite_line_source_vertical from .finite_line_source_inclined import finite_line_source_inclined + +__all__ = [ +'finite_line_source', +'finite_line_source_vertical', +'finite_line_source_inclined' +] From 426910744bcaa5843af0b571fd8eb51a5d3cd231 Mon Sep 17 00:00:00 2001 From: MassimoCimmino Date: Tue, 29 Apr 2025 05:40:42 -0400 Subject: [PATCH 5/7] Add borefield to doc --- doc/source/modules.rst | 1 + doc/source/modules/borefield.rst | 10 ++++++++++ 2 files changed, 11 insertions(+) create mode 100644 doc/source/modules/borefield.rst diff --git a/doc/source/modules.rst b/doc/source/modules.rst index bfce80b6..34494e82 100644 --- a/doc/source/modules.rst +++ b/doc/source/modules.rst @@ -7,6 +7,7 @@ Modules .. toctree:: :maxdepth: 2 + modules/borefield modules/boreholes modules/gfunction modules/heat_transfer diff --git a/doc/source/modules/borefield.rst b/doc/source/modules/borefield.rst new file mode 100644 index 00000000..6444214e --- /dev/null +++ b/doc/source/modules/borefield.rst @@ -0,0 +1,10 @@ +.. borefield: + +**************** +Borefield Module +**************** + +.. automodule:: pygfunction.borefield + :members: + :undoc-members: + :show-inheritance: From e9e5f1a45283305d98db7d784c97bc1cda58c824 Mon Sep 17 00:00:00 2001 From: MassimoCimmino Date: Tue, 29 Apr 2025 07:14:32 -0400 Subject: [PATCH 6/7] Broadcast inputs of Borefield class Co-authored-by: j-c-cook --- pygfunction/borefield.py | 30 +++++++++++++++--------------- pygfunction/gfunction.py | 6 +++--- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/pygfunction/borefield.py b/pygfunction/borefield.py index 2019ce33..e8041a07 100644 --- a/pygfunction/borefield.py +++ b/pygfunction/borefield.py @@ -44,29 +44,29 @@ def __init__( self, H: npt.ArrayLike, D: npt.ArrayLike, r_b: npt.ArrayLike, x: npt.ArrayLike, y: npt.ArrayLike, tilt: npt.ArrayLike = 0., orientation: npt.ArrayLike = 0.): - self.nBoreholes = np.max( - [len(np.atleast_1d(var)) - for var in (H, D, r_b, x, y, tilt, orientation)]) + # Extract arguments + arg_names = ("H", "D", "r_b", "x", "y", "tilt", "orientation") + args = (H, D, r_b, x, y, tilt, orientation) + + # Check if arguments are broadcastable to a valid common shape + b = np.broadcast(*args) + assert b.ndim <= 1, "All inputs must be scalars or 1D arrays." + self.nBoreholes = b.size # Broadcast all variables to arrays of length `nBoreholes` - self.H = np.broadcast_to(H, self.nBoreholes) - self.D = np.broadcast_to(D, self.nBoreholes) - self.r_b = np.broadcast_to(r_b, self.nBoreholes) - self.x = np.broadcast_to(x, self.nBoreholes) - self.y = np.broadcast_to(y, self.nBoreholes) - self.tilt = np.broadcast_to(tilt, self.nBoreholes) + for name, value in zip(arg_names, args): + setattr(self, name, np.broadcast_to(value, self.nBoreholes)) # Identify tilted boreholes self._is_tilted = np.broadcast_to( np.greater(np.abs(tilt), 1e-6), self.nBoreholes) # Vertical boreholes default to an orientation of zero - if not np.any(self._is_tilted): - self.orientation = np.broadcast_to(0., self.nBoreholes) - elif np.all(self._is_tilted): - self.orientation = np.broadcast_to(orientation, self.nBoreholes) - else: - self.orientation = np.multiply(orientation, self._is_tilted) + self.orientation = np.where( + np.broadcast_to(self._is_tilted, self.nBoreholes), + orientation, + 0. + ) def __getitem__(self, key): if isinstance(key, (int, np.integer)): diff --git a/pygfunction/gfunction.py b/pygfunction/gfunction.py index 27fc3fcb..d3b43b1d 100644 --- a/pygfunction/gfunction.py +++ b/pygfunction/gfunction.py @@ -262,17 +262,17 @@ def __init__(self, borefield_or_network, alpha, time=None, self._check_inputs() # Load the chosen solver - if self.method.lower()=='similarities': + if self.method.lower() == 'similarities': self.solver = Similarities( self.borefield, 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': + elif self.method.lower() == 'detailed': self.solver = Detailed( self.borefield, 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': + elif self.method.lower() == 'equivalent': self.solver = Equivalent( self.borefield, self.network, self.time, self.boundary_condition, self.m_flow_borehole, From 3cf8b33f4cc1aec6bc8f533b48e612c28303363f Mon Sep 17 00:00:00 2001 From: MassimoCimmino Date: Tue, 29 Apr 2025 07:26:59 -0400 Subject: [PATCH 7/7] Update CHANGELOG.md --- CHANGELOG.md | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 15dc78a0..bd22e7c9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,10 @@ ## Version 2.4 (in development) +### Enhancements + +* [Issue 211](https://github.com/MassimoCimmino/pygfunction/issues/211) - Refactored the `heat_transfer` module. The `finite_line_source`, `finite_line_source_vertical` and `finite_line_source_inclined` functions now use `Borefield` objects as arguments instead of the geometrical parameters of the boreholes. This simplifies calls to these functions. + ### 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. @@ -16,7 +20,7 @@ ### Bug fixes -* [Issue 305](https://github.com/MassimoCimmino/pygfunction/issues/305) - Fixed `ClaessonJaved` to return a float when the *g*-function is a vector (i.e. when there is only one heat source). This is required for compatibility with `numpy` version `2.x`. +* [Issue 307](https://github.com/MassimoCimmino/pygfunction/issues/307) - Fixed `ClaessonJaved` to return a float when the *g*-function is a vector (i.e. when there is only one heat source). This is required for compatibility with `numpy` version `2.x`. ### Other changes