#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Thu Apr 13 16:23:46 2017 @author: rg The module which contains static funstions related to nlayer wall calculation """ from nLayerWall_1_0 import EPSILOC12, EPSILOC10, EPSILOC4, EPSILOC6, EPSILOC8, BC_KEYWORDS, NK1, NK2 from nLayerWall_1_0.testcheck.check_fun import check_range from nLayerWall_1_0 import time_list, con_time import numpy as np def temp_interfun( alf_par, wallpl, temp1, temp2, type_unit ): if type_unit == 'steady_temp' : #-------------------------- Check alf_par------------------------------------------------------------------------------- assert isinstance(alf_par, int), 'The input parameter alf_par must be integer; {0:s}'.\ format (wallpl.__name__) if not any ((alf_par == 0, alf_par == 1, alf_par == 2, alf_par == 12, )) : print (alf_par) raise Exception ('The input parameter alf_par could be: 0 or 1 or 2 or or 12; {0:s}'.\ format (temp_interfun.__name__)) #----------------------------------------------------------------------------------------------------------------------- if temp1 != temp2 : tinter = [temp1] t0 = temp1 if alf_par == 0 : uudt = wallpl.uu_noalf * (temp1 - temp2) elif alf_par == 12 : uudt = wallpl.uu_alf12 * (temp1 - temp2) elif alf_par == 1 : uudt = wallpl.uu_alf1 * (temp1 - temp2) elif alf_par == 2 : uudt = wallpl.uu_alf2 * (temp1 - temp2) if alf_par == 12 or alf_par == 1 : tact = t0 - uudt / wallpl.alf1 check_range( 0.0, tact, temp1, temp2 ) tinter.append(tact) t0 = tact # print ('t0, temp1, temp2, wallpl.uu_noalf, uudt, wallpl.nlayer = ', t0, temp1, temp2, wallpl.uu_noalf, # uudt, wallpl.nlayer) # cond_sum = 0.0 for i in range(wallpl.nlayer): if i < wallpl.nlayer -1 : tact = t0 - uudt * ( wallpl.inv_th_res[i] + wallpl.tickradii[i] / wallpl.conductivity[i] ) # cond_sum += wallpl.inv_th_res[i] + wallpl.tickradii[i] / wallpl.conductivity[i] elif i == wallpl.nlayer -1 : tact = t0 - uudt * wallpl.tickradii[i] / wallpl.conductivity[i] # cond_sum += wallpl.tickradii[i] / wallpl.conductivity[i] # uu_act = 1.0/ cond_sum # print ('==============================================================================================') # print ('i, tact, t0, uu_act, cond_sum * uudt =' , i, tact, t0, uu_act, cond_sum * uudt ) # print ('==============================================================================================') check_range( 0.0, tact, temp1, temp2 ) tinter.append(tact) t0 = tact if alf_par == 12 or alf_par == 2 : tact = t0 - uudt / wallpl.alf2 check_range( 0.0, tact, temp1, temp2 ) tinter.append(tact) if alf_par == 12 : assert len(tinter) == wallpl.nlayer + 3 , 'Lenght of output templist in temp_inter for alf_par=12 must'\ 'be nlayer + 3' elif alf_par == 1 or alf_par == 2 : assert len(tinter) == wallpl.nlayer + 2 , 'Lenght of output templist in temp_inter for alf_par=1 or '\ 'alf_par=2 must be nlayer + 2' elif alf_par == 0 : assert len(tinter) == wallpl.nlayer + 1 , 'Lenght of output templist in temp_inter for alf_par= 0 must'\ 'be nlayer + 1' else: tinter = temp1 return tinter elif type_unit == 'uf12' : w_unit = 0.0 c_unit_new = 0.0 d_cond_temp = 0 q_unit_list = [0.0] c_unit_list = [0.0] d_cond_list = [0.0] for dens_i, sheat_i, cond_i, thick_i in zip ( wallpl.density, wallpl.specheat, wallpl.conductivity, wallpl.tickradii) : w_unit += dens_i * sheat_i * thick_i c_unit_new += thick_i * q_unit_list[-1] / cond_i + 0.5 * (thick_i**2) * dens_i * sheat_i / cond_i d_cond_temp += thick_i / cond_i q_unit_list.append(w_unit) c_unit_list.append(c_unit_new) d_cond_list.append(d_cond_temp) w_unit = 1.0 / w_unit n_qunit = len (q_unit_list) n_cunit = len (c_unit_list) q_unit_list = [ w_unit * q_unit_list[i] for i in range(n_qunit) ] c_unit_list = [ w_unit * c_unit_list[i] for i in range(n_cunit) ] assert abs(q_unit_list[-1] - 1.0) <= EPSILOC10 , 'q_unit_list[-1] must be 1.0' return w_unit, q_unit_list, c_unit_list, d_cond_list def tempxfun (x, wallpl, bc_ext, which = None , mode = None, listtemp = None, alf_par = None, cut_fl = False, steady = False ): if wallpl.lsf or wallpl.calc_dict['pure_f12'] : cut_fl = False # if wallpl.pure_f12 : # cut_fl = False if np.allclose (1.0, listtemp[0] , rtol = EPSILOC6, atol = EPSILOC8 ) and np.allclose (0.0, listtemp[-1], rtol = EPSILOC6, atol = EPSILOC8 ) : f12 = 1 elif np.allclose (1.0, listtemp[-1], rtol = EPSILOC6, atol = EPSILOC8 ) and np.allclose (0.0, listtemp[0], rtol = EPSILOC6, atol = EPSILOC8 ) : f12 = 2 #and not wallpl.lsf if steady or bc_ext == BC_KEYWORDS [0] or bc_ext == BC_KEYWORDS [4] : #bc_ext == 'steady_temp' ''' steadu sate temperature distribution inside n-layer wall also this part is used for calculation of U- functions f1 and f2 for 'd' and 'm0' BC ''' #-------------------------- Check alf_par------------------------------------------------------------------------------- assert isinstance(alf_par, int), 'The input parameter alf_par must be integer; {0:s}'.\ format (tempxfun.__name__) if not any ((alf_par == 0, alf_par == 1, alf_par == 2, alf_par == 12, )) : print ('alf_par = ', alf_par) raise Exception ('The input parameter alf_par could be: 0 or 1 or 2 or or 12 ; {0:s}'.\ format (tempxfun.__name__)) #----------------------------------------------------------------------------------------------------------------------- temp_x_out = None assert isinstance(x, int) or isinstance(x, float), 'Coordinate x must be number' x = float(x) assert 0.0 - EPSILOC12 <= x and x <= wallpl.dwall + EPSILOC12, 'Coordinate x is out of range' if alf_par == 12 : iact= 1 uudt = wallpl.uu_alf12 * (listtemp [0] - listtemp[-1]) elif alf_par == 1 : iact= 1 uudt = wallpl.uu_alf1 * (listtemp [0] - listtemp[-1]) elif alf_par == 2 : iact= 0 uudt = wallpl.uu_alf2 * (listtemp [0] - listtemp[-1]) elif alf_par == 0 : iact= 0 uudt = wallpl.uu_noalf * (listtemp [0] - listtemp[-1]) tac1 = listtemp[iact] for i in range(wallpl.nlayer): xact1 = wallpl.xrlist[ i ] xact2 = wallpl.xrlist[ i + 1 ] if xact1 - EPSILOC12 <= x and x <= xact2 + EPSILOC12 : if ( ( i == 0 and f12 == 1 ) or ( i == wallpl.nlayer - 1 and f12 == 2 ) ) and wallpl.lsf and \ not wallpl.calc_dict['pure_f12'] : if f12 == 1 : temp_x_out = np.interp ( x, wallpl.x_lsf1, wallpl.f1_lsf ) elif f12 == 2 : temp_x_out = np.interp ( x, wallpl.x_lsf2, wallpl.f2_lsf ) else : if i == 0 and cut_fl and (bc_ext == BC_KEYWORDS [0] or bc_ext == BC_KEYWORDS [4] or bc_ext == BC_KEYWORDS [2] or bc_ext == BC_KEYWORDS [5] or bc_ext == BC_KEYWORDS [6] or bc_ext == BC_KEYWORDS [7]) : corr = tac1 delta_x = x * NK1 / wallpl.tickradii[0] koeff_x = NK1 / wallpl.tickradii[0] elif i == wallpl.nlayer - 1 and cut_fl and (bc_ext == BC_KEYWORDS [0] or bc_ext == BC_KEYWORDS [4] or bc_ext == BC_KEYWORDS [3] or bc_ext == BC_KEYWORDS [5] or bc_ext == BC_KEYWORDS [6] or bc_ext == BC_KEYWORDS [8]) : corr = tac1 - uudt * wallpl.tickradii[-1] / wallpl.conductivity[-1] delta_x = ( wallpl.xrlist[ -1 ] - x ) * NK2 / wallpl.tickradii[-1] koeff_x = - NK2 / wallpl.tickradii[-1] else : corr = 0.0 delta_x = 1.0 koeff_x = 1.0 if mode == None or mode == 'temp' : temp_x_out = tac1 - uudt * ( x - xact1 ) / wallpl.conductivity[i] - corr * np.exp(- delta_x ) elif mode == 'deriv' : temp_x_out = - uudt / wallpl.conductivity[i] - corr * koeff_x * np.exp(- delta_x ) elif mode == 'flux' : # print ('cut_fl, i = ', cut_fl, i ) temp_x_out = uudt + wallpl.conductivity[i] * corr * koeff_x * np.exp(- delta_x ) break iact += 1 tac1 = listtemp[iact] assert not(temp_x_out == None), 'The calculation is failed' if ( mode == None or mode == 'temp' ) and not cut_fl : check_range( x, temp_x_out, listtemp[0], listtemp[-1] ) # print (bc_ext, cut_fl ) return temp_x_out #!"£$%^&*() #!"£$%^&*() elif wallpl.lsf : # print ('----------------------------------------------------------------------------------------------------' ) # print (' listtemp = ' , listtemp ) # print ('----------------------------------------------------------------------------------------------------' ) # rtol = EPSILOC4, atol = EPSILOC6 #!"£$%^&*() if np.allclose (1.0, listtemp[0] , rtol = EPSILOC6, atol = EPSILOC8 ) and np.allclose (0.0, listtemp[-1], #!"£$%^&*() rtol = EPSILOC6, atol = EPSILOC8 ) : #!"£$%^&*() temp_x_out = np.interp ( x, wallpl.x_lsf1, wallpl.f1_lsf ) #!"£$%^&*() elif np.allclose (1.0, listtemp[-1], rtol = EPSILOC6, atol = EPSILOC8 ) and np.allclose (0.0, listtemp[0], #!"£$%^&*() rtol = EPSILOC6, atol = EPSILOC8 ) : #!"£$%^&*() temp_x_out = np.interp ( x, wallpl.x_lsf2, wallpl.f2_lsf ) #!"£$%^&*() return temp_x_out #!"£$%^&*() elif any ((bc_ext == BC_KEYWORDS [1], bc_ext == BC_KEYWORDS [2], bc_ext == BC_KEYWORDS [3] )) : # and not wallpl.lsf # and and not wallpl.lsf : ''' f1(x): df1/dx (x=0) = -1/lambda1 , df1/dx (x=L) = 0 OR f2(x): df2/dx (x=0) = 0 , df1/dx (x=L) = -1/lambda[-1] this part is used for calculation of U- functions f1 and f2 for 'n', 'tq1' and 'tq2' BC ''' assert isinstance(x, int) or isinstance(x, float), 'Coordinate x must be number' x = float(x) assert 0.0 - EPSILOC12 <= x and x <= wallpl.dwall + EPSILOC12, 'Coordinate x is out of range' f12_out = None for i in range(wallpl.nlayer): xact1 = wallpl.xrlist[i] xact2 = wallpl.xrlist[i+1] if xact1 - EPSILOC12 <= x and x <= xact2 + EPSILOC12 : #----------------------------------------------------------------------------------------------------------------------- if mode == None or mode == 'temp' : f_temp = -0.5 * (x- wallpl.xrlist[i])**2 * wallpl.w_unit * wallpl.density[i] * \ wallpl.specheat[i] / wallpl.conductivity[i] - wallpl.q_unit_list[i] *\ (x- wallpl.xrlist[i]) / wallpl.conductivity[i] - wallpl.c_unit_list[i] elif mode == 'deriv' : f_temp = -(x- wallpl.xrlist[i]) * wallpl.w_unit * wallpl.density[i] * \ wallpl.specheat[i] / wallpl.conductivity[i] - wallpl.q_unit_list[i] / wallpl.conductivity[i] elif mode == 'flux' : f_temp = (x- wallpl.xrlist[i]) * wallpl.w_unit * wallpl.density[i] * \ wallpl.specheat[i] + wallpl.q_unit_list[i] #----------------------------------------------------------------------------------------------------------------------- if mode == None or mode == 'temp' : if which == 1 : f1_temp = - f_temp - wallpl.d_cond_list[i] - (x- wallpl.xrlist[i]) / wallpl.conductivity[i] if bc_ext == BC_KEYWORDS [1] : f12_out = f1_temp elif bc_ext == BC_KEYWORDS [2] : f12_out = f1_temp + 1.0 elif bc_ext == BC_KEYWORDS [3] : f12_out = f1_temp - wallpl.c_unit_list[wallpl.nlayer] + wallpl.d_cond_list[wallpl.nlayer] elif which == 2 : if bc_ext == BC_KEYWORDS [1] or bc_ext == BC_KEYWORDS [2] : f12_out = f_temp elif bc_ext == BC_KEYWORDS [3] : f12_out = f_temp + wallpl.c_unit_list[wallpl.nlayer] + 1.0 elif mode == 'deriv' : if which == 1 : f12_out = - f_temp - 1.0 / wallpl.conductivity[i] elif which == 2 : f12_out = f_temp elif mode == 'flux' : if which == 1 : f12_out = - f_temp + 1.0 elif which == 2 : f12_out = f_temp break assert not(f12_out == None), 'The calculation is failed' # check_range( temp_x_out, listtemp[0], listtemp[-1] ) return f12_out elif bc_ext == BC_KEYWORDS [5] : #and not wallpl.lsf ''' BC - md1 ''' pass elif bc_ext == BC_KEYWORDS [6] : #and not wallpl.lsf ''' BC - md2 ''' pass elif bc_ext == BC_KEYWORDS [7] : #and not wallpl.lsf ''' BC - mn1 ''' pass elif bc_ext == BC_KEYWORDS [8] : #and not wallpl.lsf ''' BC - mn2 ''' pass #======================================================================================================================= def converse (unit_out, unit_read ): if unit_read == '[K]' and unit_out == '[C]' : conv_ax = 1.0 conv_bx = - 273.15 elif unit_read == '[C]' and unit_out == '[K]' : conv_ax = 1.0 conv_bx = 273.15 if unit_out in time_list and unit_read in time_list : for i, time_unit in enumerate(time_list): if unit_read == time_unit: index_read = i if unit_out == time_unit: index_out = i conv_ax = con_time[index_read]/con_time[index_out] conv_bx = 0.0 if unit_read == unit_out : conv_ax = 1.0 conv_bx = 0.0 return conv_ax, conv_bx #======================================================================================================================= def smoothing (dd_series, nn, smoot_type, bc_type ): # -------------- Lanczos smoothing -------------------------------------------------------------------------------------- cor_lan = [] for j in range(nn): if smoot_type == 1 : if bc_type == BC_KEYWORDS [1] : if j != 0 : cor_lan.append( np.sin( j * np.pi / nn ) / ( j * np.pi / nn ) ) else : cor_lan.append( 1.0 ) else : cor_lan.append( np.sin( (j + 1) * np.pi / nn ) / ( (j + 1) * np.pi / nn ) ) elif smoot_type == 2 : if bc_type == BC_KEYWORDS [1] : cor_lan.append( ( 1.0 + np.cos( j * np.pi / nn ) ) / 2.0 ) else : cor_lan.append( ( 1.0 + np.cos( (j + 1) * np.pi / nn ) ) / 2.0 ) cor_lan = np.array( cor_lan ) dd_series = dd_series * cor_lan return dd_series def q_wall (wallpl, flux = None, area_heat = None, temp1 = None) : if flux == None : raise Exception ( 'Missing input parameterflux' ) if area_heat == None and temp1 != None : area_heat = temp1 * wallpl.roce_ekv * wallpl.dwall - flux * wallpl.heat_corr return_quan = area_heat else : temp1 = ( area_heat + flux * wallpl.heat_corr ) / ( wallpl.roce_ekv * wallpl.dwall ) return_quan = temp1 return return_quan # Integration Nodes used in non-uniform 'trapez' integration # Maps x-points (0,dx, 2*dx,..., N*dx ) to integrtated nodes ( 0, t1, t2, .., tM) # time_start&end - statr and end points in time axes def int_nodes (time_step, x_end, time_start, time_end) : x_setp = np.log ( ( np.exp ( x_end ) - 1.0 ) * time_step / ( time_end- time_start ) + 1.0 ) nn_x = int( x_end / x_setp ) + 1 hh = ( time_end - time_start ) / ( np.exp( x_end ) - 1.0 ) bb = time_start - hh x_nodes = np.linspace(0, x_end, nn_x ) return hh * np.exp( x_nodes ) + bb def fun_iter (x, a, nx ) : df = 1.0 / ( 1.0 - np.exp(-x ) * ( a - 1.0 ) / a ) y = ( np.log( ( a * np.exp( x ) - 1.0 ) + 1.0 ) - df * x ) / ( nx - df ) return y def xdelta_fun (nx, a, x_delta0, method, tool ): flag = True x_delta1 = x_delta0 while flag : if method == 1 : x_delta = np.log( ( np.exp(x_delta1) - 1.0 ) * a + 1.0 ) / nx if method == 2 : x_delta = fun_iter (x_delta1, a, nx ) x_delta1 = x_delta if np.abs ((x_delta - x_delta1) / x_delta ) <= tool : flag = False return x_delta def xend_fun (nx, t_start, t_end, t_delta, x_delta0, method, tool ): a = ( t_end- t_start ) / t_delta if a <= nx : print ( 'a = ', a) print ( 'nx = ', nx) raise Exception(' a <= nx, in wall_fun.xend_fun ') x_delta = xdelta_fun (nx, a, x_delta0, method, tool ) # print ( '---------------------------------------------------------------------------------------------------------') # print ( 'x_delta*nx = ', x_delta * nx ) # print ( 'log_check = ', np.log ( ( np.exp(x_delta) - 1 ) * a + 1.0 ) ) # print ( '---------------------------------------------------------------------------------------------------------') return np.log ( ( np.exp(x_delta) - 1 ) * a + 1.0 ) def derivative (x_coord, k ): n_coord = x_coord.shape[0] f_coord = np.array( [ [ -1.0, 1.0] ] * 2**(k + 1) ) de_x = x_coord[1:] - x_coord[0:-1] de_x2 = de_x.reshape( de_x.shape[0], 1 ) de_x2 = np.concatenate ((de_x2, de_x2 ), axis = 1 ) f_coord = f_coord / de_x2 x_copy = x_coord # print ( f_coord.shape ) if n_coord != 2**( k + 1 ) + 1 : raise Exception('Error in input data for derivative n_coord != exp(k+1) + 1 ; {0:s}') # . \format (self.__class__.__name__)) for i in range ( k + 1 ) : add_array = np.array( [0.0]*2**i ) # print (i, x_copy, add_array.shape ) if i == 0 : # print (x_coord.shape) # print (x_coord[:-2:2].shape) # print (x_coord[2::2].shape) de_x = ( x_coord[2::2] - x_coord[:-2:2] )*0.5 # print (de_x.shape) else: de_x = x_copy[1:] - x_copy[:-1] for j , ( f1, f2 ) in enumerate (zip (f_coord[0::2] , f_coord[1::2] ) ) : # print (i, j, de_x, de_x.shape ) f1 = np.concatenate ((f1, add_array ) ) f2 = np.concatenate ((add_array, f2 ) ) # print (i, j, f1.shape, f2.shape) if i == 0 : dd = de_x [j] else : dd = de_x [2*j] f_cor1 = (f2 - f1) / dd f_cor1 = f_cor1.reshape( 1, f_cor1.shape[0]) # print (f_cor1.shape) if j == 0 : f_cor_copy = np.copy ( f_cor1 ) else : f_cor_copy = np.concatenate ( (f_cor_copy, f_cor1 ) ) # print (f_cor_copy.shape) f_coord = f_cor_copy x_copy = x_coord[2**i::2**(i+1)] # print (i, f_coord.shape) # print ( f_coord ) return f_coord [0] def ll_modify (ll_cor, k_start, green_super, ll_p, which_ll = 1, corr = 1.0, call_no = 0, mode = 1 ): # print ('1---------------------------------------------------------------------------------------------') # print ('ll_green.shape = ', ll_green.shape) #----------------------------------------------------------------------------------------------------------------------- if call_no == 0 : ll_green = np.dot( ll_cor [k_start : , : ] , green_super ) kk = corr * ll_green / green_super[ k_start + ll_p : - ll_p ] elif call_no != 0 : ll_green = np.dot( ll_cor , green_super ) kk = corr * ll_green / green_super[ ll_p : - ll_p ] if which_ll == 2 : kk_vec = kk # print ('kk.shape = ', kk.shape) # print ('kk = ', kk) #----------------------------------------------------------------------------------------------------------------------- if call_no == 0 : # kk0 = np.array ( [kk[0]]*k_start) #======================================================================================================================= if mode == 0 : k_mean = np.mean(kk[:k_start]) kk0 = np.array ( [k_mean]*k_start) kk = np.zeros ( ( kk.shape[0] ) ) # kk[:k_start] = # kk1 = np.arrat([0]*kk.shape[0]) else : kk0 = np.array ( [kk[0]]*k_start) #======================================================================================================================= # print ('kk0.shape = ', kk0.shape) kk = np.concatenate ((kk0, kk)) elif call_no != 0 and mode == 0 : kk0 = np.zeros ( ( kk.shape[0] -k_start) ) kk = np.concatenate ((kk[:k_start], kk0)) # print ('2---------------------------------------------------------------------------------------------') # print ('kk.shape = ', kk.shape) kk0 = np.zeros ( ( ll_cor.shape[0], ll_p ) ) # print ('kk0.shape = ', kk0.shape) kk = np.concatenate ( ( kk0 , np.diag ( kk ) , kk0 ), axis = 1 ) # print ('3---------------------------------------------------------------------------------------------') # print ('kk.shape = ', kk.shape) if which_ll == 1 : return kk elif which_ll == 2 : return kk, kk_vec def error_calc (vec_ref = [None], vec_calc = [None], error_type = None, start = 0 ) : l2_norm = None madp = None mae = None ll = min (vec_calc.shape[0] , vec_ref.shape[0] ) vec_1 = vec_ref[start:ll] vec_2 = vec_calc[start:ll] if error_type == 'l2_norm' or error_type == 'all_errors' : l2_norm = 100.0 * np.linalg.norm( vec_1 - vec_2 ) / np.linalg.norm ( vec_1 ) if error_type == 'madp'or error_type == 'all_errors' : madp = 100.0 * np.sum ( abs ( vec_1 - vec_2 ) ) / np.sum ( abs ( vec_1 ) ) if error_type == 'mae' or error_type == 'all_errors' : mae = np.sum ( abs ( vec_1 - vec_2 ) ) / vec_1.shape[0] return l2_norm, madp, mae