From a965f39fb2f661427e0aad3195886e71df825bb3 Mon Sep 17 00:00:00 2001 From: Massimo Cimmino Date: Sun, 30 Oct 2022 14:17:47 -0400 Subject: [PATCH 1/8] Linearize g-function at short times --- pygfunction/gfunction.py | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/pygfunction/gfunction.py b/pygfunction/gfunction.py index 4a65c0b6..7a51875a 100644 --- a/pygfunction/gfunction.py +++ b/pygfunction/gfunction.py @@ -263,8 +263,24 @@ def evaluate_g_function(self, time): # Initialize chrono tic = perf_counter() - # Evaluate g-function - self.gFunc = self.solver.solve(time, self.alpha) + # Evaluate threshold time for g-function linearization + r_b_max = np.max([b.r_b for b in self.boreholes]) + time_threshold = r_b_max**2 / (25 * self.alpha) + if np.any(time < time_threshold): + if np.all(time < time_threshold): + time_long = np.array([time_threshold]) + else: + time_long = np.concatenate( + [time_threshold, time[time >= time_threshold]]) + # Evaluate g-function + gFunc_long = self.solver.solve(time_long, self.alpha) + # Interpolate for time < time_threshold + time_short = time[time < time_threshold] + gFunc_short = gFunc_long[0] * time_short / time_threshold + self.gFunc = np.concatenate([gFunc_short, gFunc_long[1:]]) + else: + # Evaluate g-function + self.gFunc = self.solver.solve(time, self.alpha) toc = perf_counter() if self.solver.disp: From e60174e84f5d153119b174afb3f4a4143354c7e0 Mon Sep 17 00:00:00 2001 From: Massimo Cimmino Date: Sun, 30 Oct 2022 14:21:33 -0400 Subject: [PATCH 2/8] Update CHANGELOG.md --- CHANGELOG.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index b98c7fee..47abdf94 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,10 @@ * [Issue 204](https://github.com/MassimoCimmino/pygfunction/issues/204) - Added support for Python 3.9 and 3.10. [CoolProp](https://www.coolprop.org/) is removed from the dependencies and replace with [SecondaryCoolantProps](https://github.com/mitchute/SecondaryCoolantProps). +### Bug fixes + +* [Issue 231](https://github.com/MassimoCimmino/pygfunction/issues/231) - Fixed an issue where the evaluation of g-functions at very low times raises an error due a singular matrix. g-Functions below a threshold time value `t=max(r_b)**2/(25*alpha)` are now linearized. + ## Version 2.2.1 (2022-08-12) ### Bug fixes From 4f240a4a75942f1643fc9110d723e1b3a5635629 Mon Sep 17 00:00:00 2001 From: Massimo Cimmino Date: Sun, 30 Oct 2022 14:25:11 -0400 Subject: [PATCH 3/8] The minimum time is zero --- pygfunction/gfunction.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygfunction/gfunction.py b/pygfunction/gfunction.py index 7a51875a..70c99f41 100644 --- a/pygfunction/gfunction.py +++ b/pygfunction/gfunction.py @@ -275,7 +275,7 @@ def evaluate_g_function(self, time): # Evaluate g-function gFunc_long = self.solver.solve(time_long, self.alpha) # Interpolate for time < time_threshold - time_short = time[time < time_threshold] + time_short = np.maximum(time[time < time_threshold], 0.) gFunc_short = gFunc_long[0] * time_short / time_threshold self.gFunc = np.concatenate([gFunc_short, gFunc_long[1:]]) else: From 43979c5459d7bcb316819be1ae3c77f113e6427e Mon Sep 17 00:00:00 2001 From: Massimo Cimmino Date: Wed, 23 Nov 2022 17:09:01 -0500 Subject: [PATCH 4/8] Store time threshold in gFunction object --- pygfunction/gfunction.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/pygfunction/gfunction.py b/pygfunction/gfunction.py index 70c99f41..5eb0d313 100644 --- a/pygfunction/gfunction.py +++ b/pygfunction/gfunction.py @@ -214,6 +214,9 @@ def __init__(self, boreholes_or_network, alpha, time=None, self._format_inputs(boreholes_or_network) # Check the validity of inputs self._check_inputs() + # Evaluate threshold time for g-function linearization + r_b_max = np.max([b.r_b for b in self.boreholes]) + self.time_threshold = r_b_max**2 / (25 * self.alpha) # Load the chosen solver if self.method.lower()=='similarities': @@ -263,20 +266,18 @@ def evaluate_g_function(self, time): # Initialize chrono tic = perf_counter() - # Evaluate threshold time for g-function linearization - r_b_max = np.max([b.r_b for b in self.boreholes]) - time_threshold = r_b_max**2 / (25 * self.alpha) - if np.any(time < time_threshold): - if np.all(time < time_threshold): - time_long = np.array([time_threshold]) + # Linearize g-function for times under threshold + if np.any(time < self.time_threshold): + if np.all(time < self.time_threshold): + time_long = np.array([self.time_threshold]) else: time_long = np.concatenate( - [time_threshold, time[time >= time_threshold]]) + [np.array([self.time_threshold]), time[time >= self.time_threshold]]) # Evaluate g-function gFunc_long = self.solver.solve(time_long, self.alpha) # Interpolate for time < time_threshold - time_short = np.maximum(time[time < time_threshold], 0.) - gFunc_short = gFunc_long[0] * time_short / time_threshold + time_short = np.maximum(time[time < self.time_threshold], 0.) + gFunc_short = gFunc_long[0] * time_short / self.time_threshold self.gFunc = np.concatenate([gFunc_short, gFunc_long[1:]]) else: # Evaluate g-function From 2703b9a235d4b268a822c1ad2a9e72eacb9b75b9 Mon Sep 17 00:00:00 2001 From: Massimo Cimmino Date: Wed, 23 Nov 2022 17:09:12 -0500 Subject: [PATCH 5/8] Document linearized behavior --- pygfunction/gfunction.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pygfunction/gfunction.py b/pygfunction/gfunction.py index 5eb0d313..7e5b1b0d 100644 --- a/pygfunction/gfunction.py +++ b/pygfunction/gfunction.py @@ -174,6 +174,9 @@ class gFunction(object): - The 'equivalent' solver does not support the 'MIFT' boundary condition when boreholes are connected in series. - The 'equivalent' solver does not support inclined boreholes. + - The g-function is linearized for times `t < r_b**2 / (25 * self.alpha)`. + The g-function value is then interpolated between 0 and its value at the + threshold. References ---------- From 2d774b5c8af1f9558abeda8c728298f63f628810 Mon Sep 17 00:00:00 2001 From: Massimo Cimmino Date: Sat, 26 Nov 2022 15:16:49 -0500 Subject: [PATCH 6/8] Move linearization to _BaseSolver.solve --- pygfunction/gfunction.py | 114 ++++++++++++++++++++++++++------------- 1 file changed, 78 insertions(+), 36 deletions(-) diff --git a/pygfunction/gfunction.py b/pygfunction/gfunction.py index 7e5b1b0d..7f87caa0 100644 --- a/pygfunction/gfunction.py +++ b/pygfunction/gfunction.py @@ -115,6 +115,13 @@ class gFunction(object): Number of Gauss-Legendre sample points for the integral over :math:`u` in the inclined FLS solution. Default is 11. + linear_threshold : float, optional + Threshold time (in seconds) under which the g-function is + linearized. The g-function value is then interpolated between 0 + and its value at the threshold. If linear_threshold==None, the + g-function is linearized for times + `t < r_b**2 / (25 * self.alpha)`. + Default is None. disp : bool, optional Set to true to print progression messages. Default is False. @@ -217,9 +224,6 @@ def __init__(self, boreholes_or_network, alpha, time=None, self._format_inputs(boreholes_or_network) # Check the validity of inputs self._check_inputs() - # Evaluate threshold time for g-function linearization - r_b_max = np.max([b.r_b for b in self.boreholes]) - self.time_threshold = r_b_max**2 / (25 * self.alpha) # Load the chosen solver if self.method.lower()=='similarities': @@ -256,7 +260,7 @@ def evaluate_g_function(self, time): Values of the g-function """ - time = np.atleast_1d(time) + time = np.maximum(np.atleast_1d(time), 0.) assert len(time) == 1 or np.all(time[:-1] <= time[1:]), \ "Time values must be provided in increasing order." # Save time values @@ -269,22 +273,8 @@ def evaluate_g_function(self, time): # Initialize chrono tic = perf_counter() - # Linearize g-function for times under threshold - if np.any(time < self.time_threshold): - if np.all(time < self.time_threshold): - time_long = np.array([self.time_threshold]) - else: - time_long = np.concatenate( - [np.array([self.time_threshold]), time[time >= self.time_threshold]]) - # Evaluate g-function - gFunc_long = self.solver.solve(time_long, self.alpha) - # Interpolate for time < time_threshold - time_short = np.maximum(time[time < self.time_threshold], 0.) - gFunc_short = gFunc_long[0] * time_short / self.time_threshold - self.gFunc = np.concatenate([gFunc_short, gFunc_long[1:]]) - else: - # Evaluate g-function - self.gFunc = self.solver.solve(time, self.alpha) + # Evaluate g-function + self.gFunc = self.solver.solve(time, self.alpha) toc = perf_counter() if self.solver.disp: @@ -1431,6 +1421,13 @@ class _BaseSolver(object): Number of Gauss-Legendre sample points for the integral over :math:`u` in the inclined FLS solution. Default is 11. + linear_threshold : float, optional + Threshold time (in seconds) under which the g-function is + linearized. The g-function value is then interpolated between 0 + and its value at the threshold. If linear_threshold==None, the + g-function is linearized for times + `t < r_b**2 / (25 * self.alpha)`. + Default is None. disp : bool, optional Set to true to print progression messages. Default is False. @@ -1450,13 +1447,15 @@ class _BaseSolver(object): """ def __init__(self, boreholes, network, time, boundary_condition, nSegments=8, segment_ratios=utilities.segment_ratios, - approximate_FLS=False, mQuad=11, nFLS=10, disp=False, - profiles=False, kind='linear', dtype=np.double, - **other_options): + approximate_FLS=False, mQuad=11, nFLS=10, + linear_threshold=None, disp=False, profiles=False, + kind='linear', dtype=np.double, **other_options): self.boreholes = boreholes self.network = network # Convert time to a 1d array self.time = np.atleast_1d(time).flatten() + self.linear_threshold = linear_threshold + self.r_b_max = np.max([b.r_b for b in self.boreholes]) self.boundary_condition = boundary_condition nBoreholes = len(self.boreholes) # Format number of segments and segment ratios @@ -1544,6 +1543,18 @@ def solve(self, time, alpha): # Number of time values self.time = time nt = len(self.time) + # Evaluate threshold time for g-function linearization + if self.linear_threshold is None: + time_threshold = self.r_b_max**2 / (25 * alpha) + else: + time_threshold = self.linear_threshold + # Find the number of g-function values to be linearized + p_long = np.searchsorted(self.time, time_threshold) + if p_long > 0: + time_long = np.concatenate([[time_threshold], self.time[p_long:]]) + else: + time_long = self.time + nt_long = len(time_long) # Initialize g-function gFunc = np.zeros(nt) # Initialize segment heat extraction rates @@ -1556,7 +1567,7 @@ def solve(self, time, alpha): else: T_b = np.zeros((self.nSources, nt), dtype=self.dtype) # Calculate segment to segment thermal response factors - h_ij = self.thermal_response_factors(time, alpha, kind=self.kind) + h_ij = self.thermal_response_factors(time_long, alpha, kind=self.kind) # Segment lengths H_b = self.segment_lengths() if self.boundary_condition == 'MIFT': @@ -1568,7 +1579,8 @@ def solve(self, time, alpha): tic = perf_counter() # Build and solve the system of equations at all times - for p in range(nt): + p0 = max(0, p_long-1) + for p in range(nt_long): if self.boundary_condition == 'UHTR': # Evaluate the g-function with uniform heat extraction along # boreholes @@ -1577,21 +1589,21 @@ def solve(self, time, alpha): h_dt = h_ij.y[:,:,p+1] # Borehole wall temperatures are calculated by the sum of # contributions of all segments - T_b[:,p] = np.sum(h_dt, axis=1) + T_b[:,p+p0] = np.sum(h_dt, axis=1) # The g-function is the average of all borehole wall # temperatures - gFunc[p] = np.sum(T_b[:,p]*H_b)/H_tot + gFunc[p+p0] = np.sum(T_b[:,p+p0]*H_b)/H_tot else: # Current thermal response factor matrix if p > 0: - dt = self.time[p] - self.time[p-1] + dt = time_long[p] - time_long[p-1] else: - dt = self.time[p] + dt = time_long[p] # Thermal response factors evaluated at t=dt h_dt = h_ij(dt) # Reconstructed load history Q_reconstructed = self.load_history_reconstruction( - self.time[0:p+1], Q_b[:,0:p+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( @@ -1618,10 +1630,10 @@ def solve(self, time, alpha): # Solve the system of equations X = np.linalg.solve(A, B) # Store calculated heat extraction rates - Q_b[:,p] = X[0:self.nSources] + Q_b[:,p+p0] = X[0:self.nSources] # The borehole wall temperatures are equal for all segments - T_b[p] = X[-1] - gFunc[p] = T_b[p] + T_b[p+p0] = X[-1] + gFunc[p+p0] = T_b[p+p0] elif self.boundary_condition == 'MIFT': # Evaluate the g-function with mixed inlet fluid # temperatures @@ -1659,8 +1671,8 @@ def solve(self, time, alpha): # Solve the system of equations X = np.linalg.solve(A, B) # Store calculated heat extraction rates - Q_b[:,p] = X[0:self.nSources] - T_b[:,p] = X[self.nSources:2*self.nSources] + Q_b[:,p+p0] = X[0:self.nSources] + T_b[:,p+p0] = X[self.nSources:2*self.nSources] T_f_in = X[-1] # The gFunction is equal to the effective borehole wall # temperature @@ -1675,7 +1687,16 @@ def solve(self, time, alpha): self.network.cp_f) # Effective borehole wall temperature T_b_eff = T_f - 2*pi*self.network.p[0].k_s*R_field - gFunc[p] = T_b_eff + gFunc[p+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 + if not self.boundary_condition == 'UHTR': + Q_b[:,:p_long] = 1 + (Q_b[:,p_long-1:p_long] - 1) * self.time[:p_long] / time_threshold + if self.boundary_condition == 'UBWT': + T_b[:p_long] = T_b[p_long-1] * self.time[:p_long] / time_threshold + else: + T_b[:,:p_long] = T_b[:,p_long-1:p_long] * self.time[:p_long] / time_threshold # Store temperature and heat extraction rate profiles if self.profiles: self.Q_b = Q_b @@ -1920,6 +1941,13 @@ class _Detailed(_BaseSolver): Number of Gauss-Legendre sample points for the integral over :math:`u` in the inclined FLS solution. Default is 11. + linear_threshold : float, optional + Threshold time (in seconds) under which the g-function is + linearized. The g-function value is then interpolated between 0 + and its value at the threshold. If linear_threshold==None, the + g-function is linearized for times + `t < r_b**2 / (25 * self.alpha)`. + Default is None. disp : bool, optional Set to true to print progression messages. Default is False. @@ -2176,6 +2204,13 @@ class _Similarities(_BaseSolver): Number of Gauss-Legendre sample points for the integral over :math:`u` in the inclined FLS solution. Default is 11. + linear_threshold : float, optional + Threshold time (in seconds) under which the g-function is + linearized. The g-function value is then interpolated between 0 + and its value at the threshold. If linear_threshold==None, the + g-function is linearized for times + `t < r_b**2 / (25 * self.alpha)`. + Default is None. disp : bool, optional Set to true to print progression messages. Default is False. @@ -3473,6 +3508,13 @@ class _Equivalent(_BaseSolver): Number of Gauss-Legendre sample points for the integral over :math:`u` in the inclined FLS solution. Default is 11. + linear_threshold : float, optional + Threshold time (in seconds) under which the g-function is + linearized. The g-function value is then interpolated between 0 + and its value at the threshold. If linear_threshold==None, the + g-function is linearized for times + `t < r_b**2 / (25 * self.alpha)`. + Default is None. disp : bool, optional Set to true to print progression messages. Default is False. From 05ad820e2e139126d8fd0a21c2c7c3f7c586b6a4 Mon Sep 17 00:00:00 2001 From: Massimo Cimmino Date: Sat, 26 Nov 2022 16:03:50 -0500 Subject: [PATCH 7/8] Use side='right' in numpy.searchsorted --- pygfunction/gfunction.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygfunction/gfunction.py b/pygfunction/gfunction.py index 7f87caa0..2a431544 100644 --- a/pygfunction/gfunction.py +++ b/pygfunction/gfunction.py @@ -1549,7 +1549,7 @@ def solve(self, time, 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) + p_long = np.searchsorted(self.time, time_threshold, side='right') if p_long > 0: time_long = np.concatenate([[time_threshold], self.time[p_long:]]) else: From bb9506cc378f6168a58f6d584fdcb6f1257129f6 Mon Sep 17 00:00:00 2001 From: Massimo Cimmino Date: Sat, 26 Nov 2022 16:04:24 -0500 Subject: [PATCH 8/8] Test g-function linearization at low time values --- tests/gfunction_test.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/tests/gfunction_test.py b/tests/gfunction_test.py index 2cfab47f..151f1d6c 100644 --- a/tests/gfunction_test.py +++ b/tests/gfunction_test.py @@ -272,3 +272,26 @@ def test_gfunctions_UBWT(two_boreholes_inclined, method, opts, expected, request boreholes, alpha, time=time, method=method, options=options, boundary_condition='UBWT') assert np.allclose(gFunc.gFunc, expected) + + +# ============================================================================= +# Test gFunction linearization +# ============================================================================= +# Test 'UBWT' g-functions for a single borehole at low time values +@pytest.mark.parametrize("field, method, opts, expected", [ + # 'equivalent' solver - unequal segments + ('single_borehole', 'equivalent', 'unequal_segments', np.array([1.35223554e-05, 1.35223554e-04, 2.16066268e-01])), + ]) +def test_gfunctions_UBWT_linearization(field, method, opts, expected, request): + # Extract the bore field from the fixture + boreholes = request.getfixturevalue(field) + # Extract the g-function options from the fixture + options = request.getfixturevalue(opts) + alpha = 1e-6 # Ground thermal diffusivity [m2/s] + # Times for the g-function [s] + time = np.array([0.1, 1., 10.]) * boreholes[0].r_b**2 / (25 * alpha) + # g-Function + gFunc = gt.gfunction.gFunction( + boreholes, alpha, time=time, method=method, options=options, + boundary_condition='UBWT') + assert np.allclose(gFunc.gFunc, expected)