Skip to content

Commit

Permalink
Merge pull request #240 from MassimoCimmino/issue231_singularMatrix
Browse files Browse the repository at this point in the history 
Issue231 singular matrix
  • Loading branch information
@MassimoCimmino
MassimoCimmino authored Nov 26, 2022
2 parents e942304 + bb9506c commit 8681246
Show file tree
Hide file tree
Showing 3 changed files with 106 additions and 17 deletions.
4 changes: 4 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
96 changes: 79 additions & 17 deletions pygfunction/gfunction.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -174,6 +181,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
----------
Expand Down Expand Up @@ -250,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
Expand Down Expand Up @@ -1411,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.
Expand All @@ -1430,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
Expand Down Expand Up @@ -1524,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, side='right')
if p_long > 0:
time_long = np.concatenate([[time_threshold], self.time[p_long:]])
else:
time_long = self.time
nt_long = len(time_long)
# Initialize g-function
gFunc = np.zeros(nt)
# Initialize segment heat extraction rates
Expand All @@ -1536,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':
Expand All @@ -1548,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
Expand All @@ -1557,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(
Expand All @@ -1598,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
Expand Down Expand Up @@ -1639,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
Expand All @@ -1655,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
Expand Down Expand Up @@ -1900,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.
Expand Down Expand Up @@ -2156,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.
Expand Down Expand Up @@ -3453,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.
Expand Down
23 changes: 23 additions & 0 deletions tests/gfunction_test.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)

0 comments on commit 8681246

Please sign in to comment.