Add files via upload

Examples/Calc_FFT_Coeffs.py

@@ -0,0 +1,86 @@
+from __future__ import division
+
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+import sys
+import timeit
+
+from scipy.optimize import curve_fit
+from scipy.interpolate import interp1d as interp1d
+
+sourcepath = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+sys.path.append(sourcepath)
+
+import utilities 
+
+
+def calc_FFT_coeff(nt,gfunc,tlog,lin_gm,nBhe,nt_future,loads,dt):
+
+	nsteps_total = loads.size + nt_future
+	Time = np.linspace(dt,nsteps_total*dt,nsteps_total)
+
+	# convert loads
+	FFT_Loads = []
+	for i in range(0,nBhe):
+		gload = np.concatenate([loads,np.zeros(nt_future)])		
+		gload[1:] -= np.roll(gload,1)[1:]
+
+		FFT_Loads.append(utilities.FourierSci(gload))
+
+	# Temperatures at each borehole
+	T_borehole = np.ones([nBhe,Time.size])*10
+	for i in range(0,nBhe):			
+		for j in range(0,nBhe):
+
+			# FFT Gfunc 
+			fG = interp1d(tlog,gfunc)
+			FFT_Gfunc = utilities.FourierSci(fG(Time))		
+
+			# Calc Tborehole
+			T_borehole[j,:] -= np.real(utilities.invFourierSci(FFT_Loads[i]*FFT_Gfunc)) 
+
+	return True
+
+
+def main():
+
+	dt = 30
+	nmax = 250000		
+	stepsize = 7500	
+
+	steps = np.arange(stepsize,nmax,stepsize)
+	gfunc = np.random.rand(200)*10	
+	times_arr = np.zeros(steps.size)
+
+
+	for i in range(0,steps.size):
+		print(steps[i])
+		nt = steps[i]
+
+		tlog = np.geomspace(dt,3*nt*dt,200)
+		lin_g = interp1d(tlog,gfunc)
+		nBhe = 1	
+		loads = np.random.rand(nt)
+		nt_future = loads.size
+
+		times_arr[i] = np.min(timeit.repeat(lambda :calc_FFT_coeff(nt,gfunc,tlog,lin_g,nBhe,nt_future,loads,dt),repeat = 20, number = 1,globals=globals()))/1
+
+
+	print(steps)
+	popt_l, pcov = curve_fit(utilities.func_lin0,steps, times_arr)	
+
+	print('m: '+ str(popt_l[0]))
+
+	plt.plot(steps,times_arr,label = 'measured',color = 'b')
+	plt.plot(steps,popt_l[0]*steps,label = 'lin fit', color = 'r')
+
+	plt.ylabel('comp time')
+	plt.xlabel('steps')
+	plt.legend()
+	plt.show()
+
+
+# Main function
+if __name__ == '__main__':
+    main()

Examples/Calc_SA_Coeffs.py

@@ -0,0 +1,87 @@
+from __future__ import division
+
+import numpy as np
+import math
+
+import matplotlib.pyplot as plt
+import timeit
+from scipy.optimize import curve_fit
+import os
+import sys
+
+sourcepath = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+sys.path.append(sourcepath)
+
+import utilities 
+
+
+
+def nsec(nt,psi_FFT,psi_sa):
+	return np.sqrt(psi_sa/psi_FFT*nt-1)
+
+
+def calc_sa_coeff(T_borehole,loads,gMatrix,i,nSteps_soil,nBhe):
+	for j in range(0,nBhe):
+		for k in range(0,nBhe):	
+			T_borehole[k][i:nSteps_soil] -= (loads[j][i]-loads[j][i-1]) * gMatrix[j,k,0:nSteps_soil-i]
+
+def main():
+
+	nt = 250000
+	stepsize = 1000
+	nBhe = 1
+	nb = 130000
+
+	# random values for testing
+	T_borehole = np.random.rand(nBhe,nt)
+	gMatrix = np.random.rand(nBhe,nBhe,nt)
+	loads = np.random.rand(nBhe,nt)	
+	times_arr = np.zeros(int(nt/stepsize)-1)
+
+	for j in range(stepsize,nt,stepsize):
+		# measure computational times
+		times_vec = timeit.repeat(lambda: calc_sa_coeff(T_borehole,loads,gMatrix,j,nt,nBhe),repeat =50, number = 10,globals=globals())	
+		times_arr[int(j/stepsize)-1] = np.min(times_vec)/10		
+
+	times_arr = np.flip(times_arr)
+	steps = np.arange(stepsize,nt,stepsize)
+
+	# curve fitting 
+	# fit nr<nb
+	popt_l, pcov = curve_fit(utilities.func_lin,steps[steps<nb], times_arr[steps<nb])
+	m_low = popt_l[0]
+	n_low = popt_l[1]
+	print('m_low:, n_low: %d',m_low,n_low)
+
+	# fit nr>nb
+	popt_l, pcov = curve_fit(utilities.func_lin,steps[steps>nb], times_arr[steps>nb])
+	m_high = popt_l[0]
+	n_high = popt_l[1]
+	print('m_high:, n_high %d',m_high,n_high)
+
+	# - - - - - - - - - - - - - - - - - - - - - - - 
+	#				plot results
+	# - - - - - - - - - - - - - - - - - - - - - - - 
+
+	fig = plt.figure(figsize=utilities.cm2inch(10, 6))
+	ax1 = fig.add_subplot(111)	
+	plt.subplots_adjust(bottom=0.22, top=0.90, left=0.18, right = 0.94) 
+
+	ax1.set_ylabel(r'$\mathrm{Computation \: time \: [} 10^{-3} \: \mathrm{s]}$')
+	ax1.set_xlabel(r'$\mathrm{remaining \: Steps \: [-]}$')
+
+	ax1.plot(steps,1000*times_arr,label = r'$\mathrm{measured}$',color = 'k',linewidth = 0.7)
+
+	ax1.plot(steps[steps<nb],1000*(steps[steps<nb]*m_low + n_low),label = r'$\mathrm{linear \: fit}$',linestyle = '--',color = 'grey',linewidth = 0.7,marker = 'o',markersize = 4,markevery=2,markerfacecolor='none')
+	ax1.plot(steps[steps>nb],1000*(steps[steps>nb]*m_high + n_high),linestyle = '--',color = 'grey',linewidth = 0.7,marker = 'o',markersize = 4,markevery=2,markerfacecolor='none')
+
+	ax1.legend(loc=0,frameon=False)
+	plt.show()
+
+
+
+
+
+# Main function
+if __name__ == '__main__':
+    main()

Examples/Calc_opt_periods.py

@@ -0,0 +1,23 @@
+from __future__ import division
+import os
+import sys
+
+sourcepath = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+sys.path.append(sourcepath)
+
+import utilities 
+
+
+nb = 130000		# breakpoint calc time sa model
+nt = 2109906	# total number of timesteps
+
+# computational time coefficients derived from measurements
+psi_FFT = 1.58e-6		# coeff for FFT model
+psi_sa_l = 5.8e-10		# coeff for sa model below breakpoint
+psi_sa_h = 3.7e-9		# coeff for sa model above breakpoint
+c_sa_l = 5.4e-6			# coeff for sa model below breakpoint
+c_sa_h = -9.1e-5		# coeff for sa model above breakpoint
+
+
+print(utilities.opt_nper(nt,psi_sa_l,psi_FFT,c_sa_l))
+

Examples/example_error_coupling.py

@@ -0,0 +1,160 @@
+from __future__ import division
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy import interpolate
+import os
+import sys
+import time as time
+from multiprocessing import Pool
+
+sourcepath = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+sys.path.append(sourcepath)
+
+import bhe_models 
+import gfunctions  
+import utilities 
+import hybrid_model
+
+
+
+def main():
+
+
+	BheData = 		{	    'id'    : 'BHE01',
+							'type' : '2-U',
+							'length': 100.,			# [m]		length
+							'diamB':  0.152,		# [m]		borehole diameter
+							'lmG':    2.0,			# [W/m/K]	thermal conductivity grout
+							'capG':   1000000,		# [MJ/m³/K]	volumetric heat capacity grout
+							'lmP':    0.3,			# [W/m/K]	thermal conductivity pipe
+							'odiamP': 0.032,		# [m]		outer diameter pipe
+							'thickP': 0.0029,		# [m]		wall thickness pipe
+							'distP':  0.04,			# [m]		pipe distance 
+							'densF':  1045,			# [kg/m³]	density fluid	
+							'dynviscF': 0.0035,		# [m]		dynamic viscosity fluid
+							'lmF':    0.43,			# [m]		thermal conductivity fluid
+							'capF':   3800000.,		# [m]		volumetric heat capacity fluid
+							'covH':   1.,			# [m]		cover height
+							'linkedHP': 'HP01',		# [m]		linked heat pump
+							'xCoord' : 30000,		# [m]		x coordinate
+							'yCoord' : 40000,		# [m]		y coordinate
+							'Qf' : 1.92/3600.,		# [m³/s]	flow rate		
+							'Tundist': 12.0,		# [°C]		undisturbed ground temperature
+							'lm' : 2.3,				# [W/m/K]	thermal conductivity ground	
+							'Cm' : 2300000.,		# [MJ/m³/K]	volumetric heat capacity ground
+	}
+
+
+	sim_setup = {	'nz': 	5,				# bhe cells vertical
+				'dt_bhe': 	30,				# bhe time step
+				'dt_soil': 	3600,			# soil time step
+				'nt_part':	24,				# soil time steps per period
+				'n_parts':	1,				# number of periods
+				'error':	0.001,
+				'type' : 'forward'			# euler scheme
+	}	
+
+
+
+	#create random Tin and flow rate
+	M_Tin = np.ones(sim_setup['nt_part']*sim_setup['n_parts'])*50
+	M_qf = np.ones(sim_setup['nt_part']*sim_setup['n_parts'])*BheData['Qf']
+
+
+	x_pos = np.array([1])
+	y_pos = np.array([1])
+	nBhe = x_pos.size
+
+	# -------------------------------------------------------------------------
+	# calc gfunction and create time array tlog
+	# -------------------------------------------------------------------------
+
+	tstart = sim_setup['dt_soil']
+	tstop = sim_setup['dt_soil']*sim_setup['nt_part']*sim_setup['n_parts']+1
+	ntgfunc = 200
+
+	tlog = np.geomspace(tstart, tstop, ntgfunc)
+	gfuncs = gfunctions.g_Matrix(x_pos,y_pos,tlog,BheData)
+
+
+	gfunc_data = { 'gfuncs' : gfuncs,
+				   'gtime'  : tlog
+				   }
+
+	# gfunc for semi-analytical model
+	tlin = np.linspace(sim_setup['dt_soil'],sim_setup['dt_soil']*sim_setup['nt_part'],sim_setup['nt_part'])	# time array soil
+	gMatrix = np.zeros([nBhe,nBhe,sim_setup['nt_part']])										# matrix with gfuncs
+	for i in range(0,nBhe):
+		for j in range(0,nBhe):
+			fG = interpolate.interp1d(gfunc_data['gtime'],gfunc_data['gfuncs'][i,j])
+			gMatrix[i,j,:] = fG(tlin)																# gfunc matrix interpolated to linear time grid	
+
+	# -------------------------------------------------------------------------
+	# setup for simulation
+	# -------------------------------------------------------------------------
+
+	res_loads = np.zeros([nBhe,sim_setup['nt_part']*sim_setup['n_parts']])
+	res_Tins = np.zeros([nBhe,sim_setup['nt_part']*sim_setup['n_parts']]) 
+	res_Touts = np.zeros([nBhe,sim_setup['nt_part']*sim_setup['n_parts']])
+
+	M_Tin = np.array([M_Tin])
+	M_qf = np.array([M_qf])
+
+	T_borehole = np.ones([nBhe,sim_setup['nt_part']])*BheData['Tundist']						
+	T_borehole_ini = np.ones(nBhe)*BheData['Tundist']		
+
+	BHEs = hybrid_model.init_BHEs(nBhe ,BheData,sim_setup['dt_bhe'],sim_setup['nz'], sim_setup['type'])
+
+	# num first part
+	start = time.time()
+	(res_Tins[:,0:sim_setup['nt_part']],
+	res_Touts[:,0:sim_setup['nt_part']],
+	res_loads[:,0:sim_setup['nt_part']]) = hybrid_model.calc_sa_sec_U_forw(sim_setup,
+																  gMatrix,
+																  M_Tin[:,0:sim_setup['nt_part']],
+																  M_qf[:,0:sim_setup['nt_part']],
+																  T_borehole,BHEs,T_borehole_ini)
+
+	for p in range(1,sim_setup['n_parts']):
+
+		# analytical firt part
+		T_borehole  = hybrid_model.calc_FFT_sec(gfunc_data['gfuncs'],
+									 gfunc_data['gtime'],
+									 res_loads[:,0:p*sim_setup['nt_part']],
+									 sim_setup['nt_part'],
+									 sim_setup['dt_soil'],
+									 nBhe,BheData['Tundist'])
+
+		T_borehole_ini = [T_borehole[k,p*sim_setup['nt_part']-1]  for k in range(0,nBhe)]		
+
+		# num second part
+		(res_Tins[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']],
+		res_Touts[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']],
+		res_loads[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']]) = hybrid_model.calc_sa_sec_U_forw(sim_setup,
+																								 gMatrix,
+																								 M_Tin[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']],
+																								 M_qf[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']],
+																								 T_borehole[:,p*sim_setup['nt_part']:],BHEs,T_borehole_ini)
+
+	print('time: ',time.time() - start)
+
+
+
+	# -------------------------------------------------------------------------
+	# plot results
+	# -------------------------------------------------------------------------
+
+	plt.plot(M_Tin[0],label = 'Tin')
+	plt.plot(res_Touts[0],label = 'Tout')
+	plt.legend()
+	plt.show()
+
+
+
+
+
+
+# Main function
+if __name__ == '__main__':
+    main()