Compare the accuracy and speed of different load aggregation algorithms

This example compares the simulation times and the accuracy of borehole wall temperature predictions of different load aggregation algorithms implemented into the load aggregation module.

The g-function of a single borehole is first calculated. Then, the borehole wall temperature variations are calculated using the load aggregation schemes of Bernier et al. [1], Liu [2], and Claesson and Javed [3],. The time-variation of heat extraction rates is given by the synthetic load profile of Bernier et al. [1].

The following script validates the load aggregation schemes with the exact solution obtained from convolution in the Fourier domain (see ref. [4]).

The script is located in: pygfunction/examples/comparison_load_aggregation.py

  1# -*- coding: utf-8 -*-
  2""" Comparison of the accuracy and computational speed of different load
  3    aggregation algorithms.
  4
  5    The g-function of a single borehole is calculated for boundary condition of
  6    uniform borehole wall temperature along the borehole. Then, the borehole
  7    wall temperature variations resulting from a time-varying load profile
  8    are simulated using the aggregation methods of Bernier et al. (2004),
  9    Liu (2005), and Claesson and Javed (2012). Results are compared to the
 10    exact solution obtained by convolution in the Fourier domain.
 11
 12    Default parameters are used for each of the aggregation schemes.
 13
 14"""
 15from time import perf_counter
 16
 17import matplotlib.pyplot as plt
 18import numpy as np
 19from scipy.interpolate import interp1d
 20from scipy.signal import fftconvolve
 21
 22import pygfunction as gt
 23
 24
 25def main():
 26    # -------------------------------------------------------------------------
 27    # Simulation parameters
 28    # -------------------------------------------------------------------------
 29
 30    # Borehole dimensions
 31    D = 4.0             # Borehole buried depth (m)
 32    H = 150.0           # Borehole length (m)
 33    r_b = 0.075         # Borehole radius (m)
 34
 35    # Ground properties
 36    alpha = 1.0e-6      # Ground thermal diffusivity (m2/s)
 37    k_s = 2.0           # Ground thermal conductivity (W/m.K)
 38    T_g = 10.0          # Undisturbed ground temperature (degC)
 39
 40    # g-Function calculation options
 41    options = {'nSegments': 8,
 42               'disp': True}
 43
 44    # Simulation parameters
 45    dt = 3600.                  # Time step (s)
 46    tmax = 20.*8760. * 3600.    # Maximum time (s)
 47    Nt = int(np.ceil(tmax/dt))  # Number of time steps
 48    time = dt * np.arange(1, Nt+1)
 49
 50    # Evaluate heat extraction rate
 51    Q_b = synthetic_load(time/3600.)
 52
 53    # Load aggregation schemes
 54    ClaessonJaved = gt.load_aggregation.ClaessonJaved(dt, tmax)
 55    MLAA = gt.load_aggregation.MLAA(dt, tmax)
 56    Liu = gt.load_aggregation.Liu(dt, tmax)
 57    LoadAggSchemes = [ClaessonJaved, MLAA, Liu]
 58    loadAgg_labels = ['Claesson and Javed', 'MLAA', 'Liu']
 59    loadAgg_lines = ['b-', 'k--', 'r-.']
 60
 61    # -------------------------------------------------------------------------
 62    # Calculate g-function
 63    # -------------------------------------------------------------------------
 64
 65    # The field contains only one borehole
 66    borehole = gt.boreholes.Borehole(H, D, r_b, x=0., y=0.)
 67    # Evaluate the g-function on a geometrically expanding time grid
 68    time_gFunc = gt.utilities.time_geometric(dt, tmax, 50)
 69    # Calculate g-function
 70    gFunc = gt.gfunction.gFunction(
 71        borehole, alpha, time=time_gFunc, options=options)
 72
 73    # -------------------------------------------------------------------------
 74    # Simulation
 75    # -------------------------------------------------------------------------
 76    nLoadAgg = len(LoadAggSchemes)
 77    T_b = np.zeros((nLoadAgg, Nt))
 78
 79    t_calc = np.zeros(nLoadAgg)
 80    for n, (LoadAgg, label) in enumerate(zip(LoadAggSchemes, loadAgg_labels)):
 81        print(f'Simulation using {label} ...')
 82        # Interpolate g-function at required times
 83        time_req = LoadAgg.get_times_for_simulation()
 84        gFunc_int = interp1d(np.hstack([0., time_gFunc]),
 85                             np.hstack([0., gFunc.gFunc]),
 86                             kind='cubic',
 87                             bounds_error=False,
 88                             fill_value=(0., gFunc.gFunc[-1]))(time_req)
 89        # Initialize load aggregation scheme
 90        LoadAgg.initialize(gFunc_int / (2 * np.pi * k_s))
 91
 92        tic = perf_counter()
 93        for i in range(Nt):
 94            # Increment time step by (1)
 95            LoadAgg.next_time_step(time[i])
 96
 97            # Apply current load
 98            LoadAgg.set_current_load(Q_b[i]/H)
 99
100            # Evaluate borehole wall temperature
101            deltaT_b = LoadAgg.temporal_superposition()
102            T_b[n,i] = T_g - deltaT_b
103        toc = perf_counter()
104        t_calc[n] = toc - tic
105
106    # -------------------------------------------------------------------------
107    # Calculate exact solution from convolution in the Fourier domain
108    # -------------------------------------------------------------------------
109
110    # Heat extraction rate increment
111    dQ = np.zeros(Nt)
112    dQ[0] = Q_b[0]
113    dQ[1:] = Q_b[1:] - Q_b[:-1]
114    # Interpolated g-function
115    g = interp1d(time_gFunc, gFunc.gFunc)(time)
116
117    # Convolution in Fourier domain
118    T_b_exact = T_g - fftconvolve(
119        dQ, g / (2.0 * np.pi * k_s * H), mode='full')[0:Nt]
120
121    # -------------------------------------------------------------------------
122    # plot results
123    # -------------------------------------------------------------------------
124
125    # Configure figure and axes
126    fig = gt.utilities._initialize_figure()
127
128    ax1 = fig.add_subplot(311)
129    # Axis labels
130    ax1.set_xlabel(r'$t$ [hours]')
131    ax1.set_ylabel(r'$Q_b$ [W]')
132    gt.utilities._format_axes(ax1)
133    hours = np.array([(j+1)*dt/3600. for j in range(Nt)])
134    ax1.plot(hours, Q_b)
135
136    ax2 = fig.add_subplot(312)
137    # Axis labels
138    ax2.set_xlabel(r'$t$ [hours]')
139    ax2.set_ylabel(r'$T_b$ [degC]')
140    gt.utilities._format_axes(ax2)
141    for T_b_n, line, label in zip(T_b, loadAgg_lines, loadAgg_labels):
142        ax2.plot(hours, T_b_n, line, label=label)
143    ax2.plot(hours, T_b_exact, 'k.', label='exact')
144    ax2.legend()
145
146    ax3 = fig.add_subplot(313)
147    # Axis labels
148    ax3.set_xlabel(r'$t$ [hours]')
149    ax3.set_ylabel(r'Error [degC]')
150    gt.utilities._format_axes(ax3)
151    for T_b_n, line, label in zip(T_b, loadAgg_lines, loadAgg_labels):
152        ax3.plot(hours, T_b_n - T_b_exact, line, label=label)
153    # Adjust to plot window
154    plt.tight_layout()
155
156    # -------------------------------------------------------------------------
157    # Print performance metrics
158    # -------------------------------------------------------------------------
159
160    # Maximum errors in evaluation of borehole wall temperatures
161    maxError = np.array([np.max(np.abs(T_b_n-T_b_exact)) for T_b_n in T_b])
162    # Print results
163    print('Simulation results')
164    for label, maxError_n, t_calc_n in zip(loadAgg_labels, maxError, t_calc):
165        print()
166        print((f' {label} ').center(60, '-'))
167        print(f'Maximum absolute error : {maxError_n:.3f} degC')
168        print(f'Calculation time : {t_calc_n:.3f} sec')
169
170    return
171
172
173def synthetic_load(x):
174    """
175    Synthetic load profile of Bernier et al. (2004).
176
177    Returns load y (in watts) at time x (in hours).
178    """
179    A = 2000.0
180    B = 2190.0
181    C = 80.0
182    D = 2.0
183    E = 0.01
184    F = 0.0
185    G = 0.95
186
187    func = (168.0-C)/168.0
188    for i in [1, 2, 3]:
189        func += 1.0/(i*np.pi)*(np.cos(C*np.pi*i/84.0) - 1.0) \
190                          *(np.sin(np.pi*i/84.0*(x - B)))
191    func = func*A*np.sin(np.pi/12.0*(x - B)) \
192        *np.sin(np.pi/4380.0*(x - B))
193
194    y = func + (-1.0)**np.floor(D/8760.0*(x - B))*abs(func) \
195        + E*(-1.0)**np.floor(D/8760.0*(x - B)) \
196        /np.sign(np.cos(D*np.pi/4380.0*(x - F)) + G)
197    return -y
198
199
200# Main function
201if __name__ == '__main__':
202    main()

References