diff --git a/Master_process_simu_data.py b/Master_process_simu_data.py new file mode 100644 index 0000000..696f502 --- /dev/null +++ b/Master_process_simu_data.py @@ -0,0 +1,3340 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon Jun 21 13:55:10 2021 + +@author: Rebecca +""" +# %% Import Packages +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import matplotlib.dates as mdates +import matplotlib.ticker as ticker +import math +import re +import os +#import xarray as xr +from datetime import datetime, timedelta +from pathlib import Path +from decimal import Decimal +#import my functions from other Python files in the same folder +import func_readSerrData as func_serr + +# import pdb +# pdb.set_trace() + +plt.close('all') + +#%% DEFINE DIRECTORIES + +#Define Directory Paths +#get current directory +main_dir = Path(os.getcwd()) +#define parent directory +parent = main_dir.parents[0] + +#define directory holding input data +data_dir = parent.joinpath('Exp_Data_for_models') +#define directory where graphs will be saved +graph_dir = main_dir.joinpath('Diff_runs') + +#Read in Headers file +dat_headers = pd.read_csv(data_dir.joinpath('Serrano_dat_headers.csv')) + +#%% FUNCTIONS +def strToArr(comments): + tempstr = comments.split("[",1)[1] #remove everything to the left of '[' + tempstr = tempstr.rsplit("]")[0] #remove everything to the right of ']' + tempstr = tempstr.split() # split the remaining along white spaces + temparr = np.array(tempstr) + return temparr + +def check_mod(rem,divisor): + rem = round(rem,3) + rem = float(rem) + + if rem == divisor: + #if the remainder is equal to the divisor, then the remainder should really + #be zero, and it's a particularity of floating point arithmetic and binary numbers + #and the modulus operator that have made it not equal to zero + rem = 0 + + return rem + +def plotAllSndTemps(arrayname, marker='o'): +#plot all columns on y-axis, with first two columns ignored + # x-axis is number of days starting from start of simulation (column 0 in array) + fig1,ax1 = plt.subplots() + ax1.set_title("Simulated sand temperatures of all TCs") + ax1.set_xlabel("Days since start of simulation") + ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + + for col in range(2,np.size(arrayname,1)): + ax1.plot(arrayname[:,0],arrayname[:,col],marker) + + +def plotTimeSeriesEntireArray(timearray,y_array, fig1, ax1, labels, colors=['g','b','r']): +#plot all columns on y-axis, with first columns as time series + # x-axis is number of days starting from start of simulation (column 0 in array) + + for col in range(np.size(y_array,1)-1,-1,-1): + #ax2.plot(timearray[:,0],y_array[:,col],marker) #Use to get x-axis "days after simulation" + ax1.plot(timearray,y_array[:,col], label=labels[col], color = colors[col], linestyle = '--') + + + return fig1, ax1 + +def DateArrayFromSimDayHoursArray(time_in_days_array,start_date_in_datetime): + #create an empty list for the simulation time array + date_array = [0]*time_in_days_array.shape[0] + #put start date in top row of list + date_array[0] = start_date_in_datetime + for d in range(1,time_in_days_array.shape[0]): + date_array[d] = date_array[d-1] + timedelta(days=(time_in_days_array[d]-time_in_days_array[d-1])) + + return date_array + +def readOutputSimuDatafile(filename, plotall=False): + + run_number = filename.split('_',1) + run_number = run_number[0][1:] + + #read in dimensions of array + with open(filename, 'r') as f: + for line in range(0,2): + comments = f.readline() + if line==0: #read in dimensions of array + dims=comments.split('#')[1] + elif line==1: #read in start and end dates + dates=comments.split('#')[1] + + #Two possible headers for the output files...check the header + if dims[0].isdigit() == False: #If the first character after the # is a digit, then it's the old type of header + c=0 + #find the first character that's a digit + while (dims[c].isdigit()==False): + c=c+1 + #get rid of everything to the left of that character + dims = dims[c:] + + [X,Y,Z] =dims.split(',') + X=int(X) + Y=int(Y) + Z=int(Z) + + if dates[0].isdigit() == False: #If the first character after the # is a digit, then it's the old type of header + c=0 + #find the first character that's a digit + while (dates[c].isdigit()==False): + c=c+1 + #get rid of everything to the left of that character + dates = dates[c:] + + dates = dates.split(',') + start_date = datetime.strptime(dates[0],'%Y-%m-%d') + end_date = dates[1][:dates[1].rindex('\n')] #slices from the beginning to the '\n' + # end_date = datetime.strptime(end_date,'%Y-%m-%d') + end_date = datetime.strptime(end_date,'%Y-%m-%d') + + + #Read data from simulation into holding variable C + C = np.loadtxt(filename, delimiter=',', comments='#') + #Gets number of rows in file + num_rowsC = np.size(C,0) + #extract time data from file + timearrayout = C[:,0:2] + #extract other data from file and "unflatten/reshape" into an array + mtrxOut = np.reshape(C[:,2:],(num_rowsC,X,Y,Z),order='F') + + if plotall == True: + #Plot sand Temps closest to TCs (delete this once the code works) + plotAllSndTemps(C,'-') + + + return mtrxOut, timearrayout, start_date, end_date, run_number + +def readOutputSimuVizDatafile(filename): + + run_number = filename.split('_',1) + run_number = run_number[0][1:] + + #read in dimensions of array + with open(filename, 'r') as f: + df_out = pd.DataFrame() + for line in range(0,11): + comments = f.readline() + if line==0: #read in dimensions of array + dims=comments.split(" ")[1] + elif line==1: #read in start and end dates + dates=comments.split(" ")[1] + elif line==2: #read in number of simulated days + simu_days= int(comments.split(" ")[-1]) + elif line>=3: #read in the rest of the lines into a dataframe + df = pd.DataFrame({ line: strToArr(comments)}) + #put each line from 4-11 into a new dataframe (column of strings) + #make the line number the column heading + df_out = pd.concat([df_out,df], axis=1) + #vertically concatenate each new column to the main DataFrame, df_out + + + + [X,Y,Z] =dims.split(',') + X=int(X) + Y=int(Y) + Z=int(Z) + + dates = dates.split(',') + start_date = datetime.strptime(dates[0],'%Y-%m-%d') + end_date = dates[1][:dates[1].rindex('\n')] #slices from the beginning to the '\n' + end_date = datetime.strptime(end_date,'%Y-%m-%d') + + #Read data from simulation into holding variable C + C = np.loadtxt(filename, delimiter=',', comments='#') + #Gets number of rows in file + num_rowsC = np.size(C,0) + #extract time data from file + timearrayout = C[:,0:2] + #extract other data from file and "unflatten/reshape" into an array + mtrxOut = np.reshape(C[:,2:],(num_rowsC,X,Y,Z),order='F') + + + # xcoords = mtrxOut.shape[1] + # ycoords = mtrxOut.shape[2] + # zcoords = mtrxOut.shape[3] + return mtrxOut, timearrayout, start_date, end_date, df_out[3].dropna(), df_out[4].dropna(), df_out[5].dropna(), df_out[6].dropna(), df_out[7].dropna(), df_out[8].dropna(), df_out[9].dropna(), df_out[10].dropna(), simu_days, run_number + +def readSerranoPerformanceMetricsData(filename1, dat_headers=dat_headers, main_dir=main_dir, data_dir=data_dir): + + + #FILENAMES OF DATA FILES + #(NB! filename1 MUST have a start AND end date, even if they are equal, ie. only one day of data) + + #Parse out dates between which to run the simulation + #find year 1 + f1_name, f1_2, f1_3 = re.split('(_20)',filename1) + + st1 = filename1.find('_20') + date_start = filename1[st1+1:st1+11] + exp_start_date = datetime.date(datetime.strptime(date_start,'%Y-%m-%d')) + + end1 = filename1.find('to-') + date_end = filename1[end1+3:-4] + exp_end_date = datetime.date(datetime.strptime(date_end,'%Y-%m-%d')) + + #Read in sand temps data into a DataFrame + Perf_metrics = func_serr.read_performance_metrics_data(filename1, dat_headers, main_dir, data_dir) + + return Perf_metrics, exp_start_date, exp_end_date + + +def readSerranoSandMultipleDates(filename1, dat_headers=dat_headers, main_dir=main_dir, data_dir=data_dir): + + + #FILENAMES OF DATA FILES + #(NB! filename1 MUST have a start AND end date, even if they are equal, ie. only one day of data) + + #Parse out dates between which to run the simulation + #find year 1 + f1_name, f1_2, f1_3 = re.split('(_20)',filename1) + + st1 = filename1.find('_20') + date_start = filename1[st1+1:st1+11] + exp_start_date = datetime.date(datetime.strptime(date_start,'%Y-%m-%d')) + + end1 = filename1.find('to-') + date_end = filename1[end1+3:-4] + exp_end_date = datetime.date(datetime.strptime(date_end,'%Y-%m-%d')) + + #Read in sand temps data into a DataFrame + Exp_sand_temps = func_serr.read_sand_temps(filename1, dat_headers, main_dir, data_dir) + + return Exp_sand_temps, exp_start_date, exp_end_date + +def readSerranoSandSingleDate(filename1, dat_headers=dat_headers, main_dir=main_dir, data_dir=data_dir): + + #Define Directory Paths + #get current directory + main_dir = Path(os.getcwd()) + #define parent directory + parent = main_dir.parents[0] + + #define directory holding input data + data_dir = parent.joinpath('Exp_Data_for_models') + + #Parse out date from filename + f1_name, f1_2, f1_3 = re.split('(_20)',filename1) + + date_start = '20' + f1_3 + date_start = date_start.rstrip('.dat') + exp_start_date = datetime.strptime(date_start,'%Y-%m-%d') + + + #Read in Headers file + dat_headers = pd.read_csv(data_dir.joinpath('Serrano_dat_headers.csv')) + + #Read in sand temps data into a DataFrame + Exp_sand_temps = func_serr.read_sand_temps(filename1, dat_headers, main_dir, data_dir) + + return Exp_sand_temps, exp_start_date + +def processExpTCDataToArrays(exp_TC_data,exp_start_date): + # Process Experimental Sand Temp data into Numpy arrays + """ + This script processes the experimental data, output by Serrano, from a .dat + or .csv file into matrix form, so that it can be used later + + RECALL!!!! k = 0 is bottom of sand store, while C-layer is bottom-most plane + """ + + #Put time data into a new numpy array (shows the time in hours from the start of the year) + time_h = exp_TC_data['time'].copy().to_numpy() + + #rows and columns of sand temp data + [A, B] = exp_TC_data.shape + + #create 4-D array to hold temp data + Exp_sand_temps_4d = np.zeros(shape=[A,6,6,3]) + + #create loop to map expt'l data to 4-D matrix + for indx1 in range(0,A): + + indx2 = 1 #index for the number of columns in the matrix. It starts at 1 since the first column is the time-stamp (zero-indexed) + + for k in range (2,-1,-1): + for i in range (0,6): + for j in range (0,6): + Exp_sand_temps_4d[indx1,i,j,k] = exp_TC_data.iloc[indx1,indx2] + + #increment the column to be read next time + indx2 = indx2 + 1 + #check that it exists, and that it's not the end of the line + if indx2 >= B: + #break out of the innermost loop + break + + + + #Create a time list for the experimental time array + #create empty list + exp_time_array = [0]*A + #put start date (rounded down to midnight) in top row of list + exp_time_array[0] = datetime(year=exp_start_date.year,month=1,day=1) + timedelta(days=np.floor((time_h[0]/24))) + #get difference between start datetime and midnight of that day + start_offset = timedelta(days=(time_h[0]/24)) - timedelta(days=np.floor((time_h[0]/24))) + + start_year = exp_start_date.year + for d in range(1,A): + #if the next entry is less than the previous entry, and their absolute difference is many months (>5000h) + if (((time_h[d]-time_h[d-1]) < 0) and (abs(time_h[d]-time_h[d-1]) > 8000)): + start_year += 1 + exp_time_array[d] = datetime(year=start_year,month=1,day=1) + timedelta(days=(time_h[d]/24)) - start_offset + else: + exp_time_array[d] = datetime(year=start_year,month=1,day=1) + timedelta(days=(time_h[d]/24)) - start_offset + + return Exp_sand_temps_4d, exp_time_array + +def plotLayersofExpSandTempData_Sbplts(exp_time_array,Exp_sand_temps): + + titles = ['Experimental top layer sand temperatures', + 'Experimental middle layer sand temperatures', + 'Experimental bottom layer sand temperatures'] + layer_name = ['A','B','C'] + + for layer in range(0,3): + fig4,axs = plt.subplots(3,2,sharex=True, sharey=True) #6 SUBPLOTS + fig4.suptitle(titles[layer]) + + sbplt = 0 + #plot top layer of exp sand temps + for a in range(0,2): + for b in range (0,3): + for j in range(0,6): + col = (sbplt*6 + j) + 36*layer + axs[b,a].plot(exp_time_array,Exp_sand_temps.iloc[:,col+1],label=Exp_sand_temps.columns[col+1]) #plot top layer of sand store + axs[b,a].grid(True) + axs[b,a].set_title(layer_name[layer] + '-' + str(sbplt+1) + '-X', fontsize=8) + axs[b,a].set_ylim(20,60) + # axs[b,a].legend() + sbplt = sbplt + 1 + fig4.legend(['X-X-1', + 'X-X-2', + 'X-X-3', + 'X-X-4', + 'X-X-5', + 'X-X-6'],fontsize=8) + fig4.autofmt_xdate() #automatically makes the x-labels rotate + # locator = mdates.AutoDateLocator() + # formatter = mdates.ConciseDateFormatter(locator) + fig4.text(0.06, 0.5, 'Temperature ($\degree$C)', ha='center', va='center', rotation='vertical') + return + +def plotTopLayerofExpSandData_IndPlts(exp_time_array,Exp_sand_temps): +#this function plots 6 graphs with all the indivual TC data, 6 lines per graph, representing each row of TCs from north to south + global exp_start_date,exp_end_date + + #plot top layer of exp sand temps + for a in range(0,6): + fig,ax = plt.subplots() + for j in range(0,6): + col = a*6 + j + ax.plot(exp_time_array,Exp_sand_temps.iloc[:,col+1], label=Exp_sand_temps.columns[col+1]) #plot top layer of sand store + + #locator = mdates.DayLocator(interval=2) + locator = mdates.AutoDateLocator() + formatter = mdates.ConciseDateFormatter(locator) + ax.legend() + ax.grid(True) + ax.set_title("Sand store top layer TCs: "+ str(exp_start_date) + " to " + str(exp_end_date)) + ax.set_ylabel("Temperature ($\degree$C)") + ax.xaxis.set_major_locator(locator) + #ax.xaxis.set_minor_locator(mdates.DayLocator()) + ax.xaxis.set_major_formatter(formatter) + ax.set_ylim(20,60) + fig.autofmt_xdate() #automatically makes the x-labels rotate + return + +def plotMiddleLayerofExpSandData_IndPlts(exp_time_array,Exp_sand_temps): +#this function plots 6 graphs with all the indivual TC data, 6 lines per graph, representing each row of TCs from north to south + + global exp_start_date,exp_end_date + + #plot middle layer of exp sand temps + for a in range(0,6): + fig,ax = plt.subplots() + for j in range(0,6): + col = (a*6 + j) + (36) + ax.plot(exp_time_array,Exp_sand_temps.iloc[:,col+1], label=Exp_sand_temps.columns[col+1]) #plot top layer of sand store + + #locator = mdates.DayLocator(interval=2) + locator = mdates.AutoDateLocator() + formatter = mdates.ConciseDateFormatter(locator) + ax.legend() + ax.grid(True) + ax.set_title("Sand store middle layer TCs: "+ str(exp_start_date) + " to " + str(exp_end_date)) + ax.set_ylabel("Temperature ($\degree$C)") + ax.xaxis.set_major_locator(locator) + #ax.xaxis.set_minor_locator(mdates.DayLocator()) + ax.xaxis.set_major_formatter(formatter) + ax.set_ylim(20,60) + fig.autofmt_xdate() #automatically makes the x-labels rotate + return + +def plotBottomLayerofExpSandData_IndPlts(exp_time_array,Exp_sand_temps): +#this function plots 6 graphs with all the indivual TC data, 6 lines per graph, representing each row of TCs from north to south + + global exp_start_date,exp_end_date + + #plot middle layer of exp sand temps + for a in range(0,6): + fig,ax = plt.subplots() + for j in range(0,6): + col = (a*6 + j) + (36*2) + ax.plot(exp_time_array,Exp_sand_temps.iloc[:,col+1], label=Exp_sand_temps.columns[col+1]) #plot top layer of sand store + + #locator = mdates.DayLocator(interval=2) + locator = mdates.AutoDateLocator() + formatter = mdates.ConciseDateFormatter(locator) + ax.legend() + ax.grid(True) + ax.set_title("Sand store middle layer TCs: "+ str(exp_start_date) + " to " + str(exp_end_date)) + ax.set_ylabel("Temperature ($\degree$C)") + ax.xaxis.set_major_locator(locator) + #ax.xaxis.set_minor_locator(mdates.DayLocator()) + ax.xaxis.set_major_formatter(formatter) + ax.set_ylim(20,60) + fig.autofmt_xdate() #automatically makes the x-labels rotate + return + +def plotExpSandDataAvgofEWRows_IndPlts(exp_time_array,Exp_sand_temps_4d): +#this function plots 1 graph with 6 lines of data, each line is the average TC from each row of TCs (north to south) + + global exp_start_date,exp_end_date + + TC_EW_row_avgs = np.mean(Exp_sand_temps_4d,axis=2) + # layer_name = ["Bottom", + # "Middle", + # "Top"] + layer_name = [", Level C", + ", Level B", + ", Level A"] + + #plot top layer of exp sand temps + for layer in range(2,-1,-1): + fig,ax = plt.subplots() + lgnd = ["Row 1, north edge" + layer_name[layer], + "Row 2"+ layer_name[layer], + "Row 3"+ layer_name[layer], + "Row 4"+ layer_name[layer], + "Row 5"+ layer_name[layer], + "Row 6, south edge"+ layer_name[layer]] + lines = ["--", + "dotted", + "-", + "-", + "-", + "-"] + for j in range(0,6): + + ax.plot(exp_time_array,TC_EW_row_avgs[:,j,layer], label=lgnd[j], linestyle=lines[j]) #plot top layer of sand store + + #locator = mdates.DayLocator(interval=2) + # locator = mdates.AutoDateLocator() + locator = mdates.MonthLocator() + formatter = mdates.ConciseDateFormatter(locator) + ax.legend() + # ax.legend([layer_name[layer] +" Row 1 (north)", + # layer_name[layer] +" Row 2", + # layer_name[layer] +" Row 3", + # layer_name[layer] +" Row 4", + # layer_name[layer] +" Row 5", + # layer_name[layer] +" Row 6 (south)"]) + # ax.legend(["Row 1 (northern edge" + layer_name[layer], + # "Row 2", + # "Row 3", + # "Row 4", + # "Row 5", + # "Row 6 (southern edge"+ layer_name[layer]]) + ax.grid(True) + #ax.set_title("Sand store "+ layer_name[layer] +" layer, average TC row temps: "+ str(exp_start_date) + " to " + str(exp_end_date)) + ax.set_ylabel("Temperature ($\degree$C)") + ax.xaxis.set_major_locator(locator) + # ax.xaxis.set_minor_locator(mdates.MonthLocator()) + ax.xaxis.set_major_formatter(formatter) + ax.set_ylim(20,60) + fig.autofmt_xdate() #automatically makes the x-labels rotate + return + +def plotExpSandDataAvgofNSRows_IndPlts(exp_time_array,Exp_sand_temps_4d): +#this function plots 1 graph with 6 lines of data, each line is the average TC from each row of TCs (north to south) + + global exp_start_date,exp_end_date + + TC_NS_row_avgs = np.mean(Exp_sand_temps_4d,axis=1) + # layer_name = ["Bottom", + # "Middle", + # "Top"] + layer_name = [", Level C", + ", Level B", + ", Level A"] + + #plot top layer of exp sand temps + for layer in range(2,-1,-1): + fig,ax = plt.subplots() + lgnd = ["Col 1, west edge" + layer_name[layer], + "Col 2"+ layer_name[layer], + "Col 3"+ layer_name[layer], + "Col 4"+ layer_name[layer], + "Col 5"+ layer_name[layer], + "Col 6, east edge"+ layer_name[layer]] + lines = ["--", + "dotted", + "-", + "-", + "-", + "-"] + for j in range(0,6): + ax.plot(exp_time_array,TC_NS_row_avgs[:,j,layer], label=lgnd[j], linestyle = lines[j]) #plot top layer of sand store + + #locator = mdates.DayLocator(interval=2) + locator = mdates.AutoDateLocator() + # locator = mdates.MonthLocator() + formatter = mdates.ConciseDateFormatter(locator) + # ax.legend([layer_name[layer] + " Column 1 (west)", + # layer_name[layer] + " Column 2", + # layer_name[layer] + " Column 3", + # layer_name[layer] + " Column 4", + # layer_name[layer] + " Column 5", + # layer_name[layer] + " Column 6 (east)"]) + ax.legend() + ax.grid(True) + #ax.set_title("Sand store "+ layer_name[layer] +" layer, average TC row temps: "+ str(exp_start_date) + " to " + str(exp_end_date)) + ax.set_ylabel("Temperature ($\degree$C)") + ax.xaxis.set_major_locator(locator) + #ax.xaxis.set_minor_locator(mdates.DayLocator()) + ax.xaxis.set_major_formatter(formatter) + ax.set_ylim(20,60) + fig.autofmt_xdate() #automatically makes the x-labels rotate + return + + +def plotExpSandDataXSections(Exp_sand_temps_4d,min_ylim,max_ylim): + + sand_xcoords = [0.5 ,1.5 ,2.5 ,3.5 ,4.5 ,5.5] + sand_ycoords = [0.5 ,1.5 ,2.5 ,3.5 ,4.5 ,5.5] + #plot x-section of all exp sand temp layers + #TOP LAYER + for j in range(0,6): + fig,ax = plt.subplots() + for time in range(0, Exp_sand_temps_4d.shape[0],288): #<- every 288 intervals (5-minutes each) is 1 day + ax.plot(sand_xcoords,Exp_sand_temps_4d[time,:,j,2],label="Day "+str(int(time/288))) #plot k=2 (top of sand store) + + ax.grid(True) + ax.set_title("Expt'l data: Cross-section of top-layer Sand Store TCs, at y=" + str(sand_ycoords[j]) + "m and z=2.2m") + ax.set_ylim(min_ylim,max_ylim) + ax.set_ylabel("Temperature ($\degree$C)") + ax.set_xlabel("Sand store x-coordinate (m), north wall XPS/sand boundary = 0m ") + ax.legend() + + #MIDDLE LAYER + for j in range(0,6): + fig,ax = plt.subplots() + for time in range(0, Exp_sand_temps_4d.shape[0],288): #<- every 288 intervals (5-minutes each) is 1 day + ax.plot(sand_xcoords,Exp_sand_temps_4d[time,:,j,1],label="Day "+str(int(time/288))) #plot k=1 (middle layer of sand store) + + ax.grid(True) + ax.set_title("Expt'l data: Cross-section of middle-layer Sand Store TCs, at y=" + str(sand_ycoords[j]) + "m and z=1.3m") + ax.set_ylim(min_ylim,max_ylim) + ax.set_ylabel("Temperature ($\degree$C)") + ax.set_xlabel("Sand store x-coordinate (m), where north wall XPS/sand boundary = 0m ") + ax.legend() + + #BOTTOM LAYER + for j in range(0,6): + fig,ax = plt.subplots() + for time in range(0, Exp_sand_temps_4d.shape[0],288): #<- every 288 intervals (5-minutes each) is 1 day + ax.plot(sand_xcoords,Exp_sand_temps_4d[time,:,j,0],label="Day "+str(int(time/288))) #plot k=0 (bottom layer of sand store) + + ax.grid(True) + ax.set_title("Expt'l data: Cross-section of bottom-layer Sand Store TCs, at y=" + str(sand_ycoords[j]) + "m and z=0.4m") + ax.set_ylim(min_ylim,max_ylim) + ax.set_ylabel("Temperature ($\degree$C)") + ax.set_xlabel("Sand store x-coordinate (m), where north wall XPS/sand boundary = 0m ") + ax.legend() + + return + +def plotExpSandDataYSections(Exp_sand_temps_4d,min_ylim,max_ylim): + + sand_xcoords = [0.5 ,1.5 ,2.5 ,3.5 ,4.5 ,5.5] + sand_ycoords = [0.5 ,1.5 ,2.5 ,3.5 ,4.5 ,5.5] + #plot y-section of all exp sand temp layers + #TOP LAYER + for i in range(0,6): + fig,ax = plt.subplots() + for time in range(0, Exp_sand_temps_4d.shape[0],288): #<- every 288 intervals (5-minutes each) is 1 day + ax.plot(sand_ycoords,Exp_sand_temps_4d[time,i,:,2],label="Day "+str(int(time/288))) #plot k=2 (top of sand store) + + ax.grid(True) + ax.set_title("Expt'l data: Cross-section of top-layer Sand Store TCs, at x=" + str(sand_xcoords[i]) + "m and z=2.2m") + ax.set_ylim(min_ylim,max_ylim) + ax.set_ylabel("Temperature ($\degree$C)") + ax.set_xlabel("Sand store y-coordinate (m), west wall XPS/sand boundary = 0m ") + ax.legend() + + #MIDDLE LAYER + for i in range(0,6): + fig,ax = plt.subplots() + for time in range(0, Exp_sand_temps_4d.shape[0],288): #<- every 288 intervals (5-minutes each) is 1 day + ax.plot(sand_ycoords,Exp_sand_temps_4d[time,i,:,1],label="Day "+str(int(time/288))) #plot k=1 (middle layer of sand store) + + ax.grid(True) + ax.set_title("Expt'l data: Cross-section of middle-layer Sand Store TCs, at x=" + str(sand_xcoords[i]) + "m and z=1.3m") + ax.set_ylim(min_ylim,max_ylim) + ax.set_ylabel("Temperature ($\degree$C)") + ax.set_xlabel("Sand store y-coordinate (m), where west wall XPS/sand boundary = 0m ") + ax.legend() + + #BOTTOM LAYER + for i in range(0,6): + fig,ax = plt.subplots() + for time in range(0, Exp_sand_temps_4d.shape[0],288): #<- every 288 intervals (5-minutes each) is 1 day + ax.plot(sand_ycoords,Exp_sand_temps_4d[time,i,:,0],label="Day "+str(int(time/288))) #plot k=0 (bottom layer of sand store) + + ax.grid(True) + ax.set_title("Expt'l data: Cross-section of bottom-layer Sand Store TCs, at x=" + str(sand_xcoords[i]) + "m and z=0.4m") + ax.set_ylim(min_ylim,max_ylim) + ax.set_ylabel("Temperature ($\degree$C)") + ax.set_xlabel("Sand store y-coordinate (m), where west wall XPS/sand boundary = 0m ") + ax.legend() + + return + + +def PlotMaxandMean(num_mod_temps,num_mod_time, run_num): + + # Plot max temp over time + + #initialize array + max_temps = np.zeros(shape=(num_mod_time.shape[0])) + + #put max temps at each timestep in the new array + for time in range(0,num_mod_time.shape[0]): + max_temps[time] = np.max(num_mod_temps[time,:,:,:]) + + #Plot the max temps + fig3,ax3 = plt.subplots() + ax3.plot(num_mod_time[:,0],max_temps, label = 'max') + + # Plot mean temp over time + #initialize array + mean_temps = np.zeros(shape=(num_mod_time.shape[0])) + + #put max temps at each timestep in the new array + for time in range(0,num_mod_time.shape[0]): + mean_temps[time] = np.mean(num_mod_temps[time,:,:,:]) + + + ax3.plot(num_mod_time[:,0],mean_temps, label ='mean') + + ax3.set_title("Run "+ run_num + ": Max and Mean temps during cooldown") + ax3.set_xlabel('Time (days)') + ax3.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + ax3.legend() + return max_temps,mean_temps + +def sim_round2hr(dtm, Nmin): + #this function rounds a datetime to the nearest 60 minutes, based on user input. + + rounding = timedelta(minutes=dtm.minute % Nmin, + seconds=dtm.second, + microseconds=dtm.microsecond) + + dtm = dtm - rounding + + #if needed, round up + if rounding >= timedelta(minutes=Nmin/2): + dtm = dtm + timedelta(minutes=Nmin) + + return dtm + + + +def exp_round2hr(): + #this function rounds a datetime to the nearest 60 minutes, based on user input. + #It also puts the corresponding average experimental temps into another array, avg_expmtl_temps4resid + + global exp_time_array, sndTemps_datetimes_rnd, exp_time_array_rn + global avg_expmtl_temps4resid, Exp_sand_temps4resid108, avg_expmtl_temps + + + + #first data point in exp arrays + exp_time_array_rnd[0] = exp_time_array[0] + avg_expmtl_temps4resid[0,:] = avg_expmtl_temps[0,:] + Exp_sand_temps4resid108 [0,:,:,:] = Exp_sand_temps_4d[0,:,:,:] + + idx = 1 + #now, for the rest of the array + for u in range(2,len(sndTemps_datetimes_rnd)): + for v in range(idx,len(exp_time_array)): + #for the very last value in the array + if u == (len(sndTemps_datetimes_rnd) - 1): + if abs(exp_time_array[-1] - sndTemps_datetimes_rnd[-1]) <= timedelta(minutes=15): + exp_time_array_rnd[-1] = exp_time_array[-1] + avg_expmtl_temps4resid[-1,:] = avg_expmtl_temps[-1,:] + Exp_sand_temps4resid108 [-1,:,:,:] = Exp_sand_temps_4d[-1,:,:,:] + + #as the loop is going through all the times in the expmt'l array, is it less than 15min on either side of the hour? + elif abs(exp_time_array[v] - sndTemps_datetimes_rnd[u]) <= timedelta(minutes=15): + #assign difference between experimental timestamp and simultd timestamp + curr_ent = abs(exp_time_array[v] - sndTemps_datetimes_rnd[u]) + next_ent = abs(exp_time_array[v+1] - sndTemps_datetimes_rnd[u]) + if next_ent > curr_ent: + #we want the current v to be stored in the RMSD array + exp_time_array_rnd[u-1] = exp_time_array[v] + avg_expmtl_temps4resid[u-1,:] = avg_expmtl_temps[v,:] + Exp_sand_temps4resid108 [u-1,:,:,:] = Exp_sand_temps_4d[v,:,:,:] + + #and bring up the lower limit that the inner loop loops through + idx = v + 1 + break + + return + + + + + +#%% 0. ASK USER WHICH DATA TO PROCESS +#Excel sheet with different run descriptions +diff_runs_file = graph_dir.joinpath('2021-07-13 Description of different runs.xlsx') +diff_runs = pd.read_excel(diff_runs_file, header=1, index_col=0) + + +#%%% Which run numbers to process? + +#single run +# runs = [399] +# i=0 +# for r in runs: +# runs[i] = str(r) +# i += 1 + +#range of runs +run_alpha = 399 #first run +run_om = 404 #last run +except_these = [] #list of runs to exclude + + +runs = [None]*(run_om - run_alpha + 1 - len(except_these)) +i=0 +for r in range(run_alpha,run_om+1): + + skip = 0 #-> variable to check if run needs to be skipped + + #check if r is any of the exceptions + for unwanted in except_these: + if r == unwanted: + skip = 1 + break + if skip == 1: + pass + #if not, do the usual - writing the run to the list and incrementing the index + else: + runs[i] = str(r) + i=i+1 + + +print('For runs: '); print(runs) +# breakpoint() +#%%% Names of simulation files with import data +df = pd.DataFrame(columns=['runs','fileID','fileID2']) + + +for i in range(0,len(runs)): + df.loc[i,'runs'] = runs[i] + df.loc[i,'fileID'] = 'r'+ runs[i] +'_1Sand_TC_Temps.txt' + df.loc[i,'fileID2'] = 'r'+ runs[i] +'_2Sand_Viz_Temps.txt' + + +#%%% Variables needed for simulation +precon_days = 10 #Number of preconditioning days in each "rewind" + +#%%% Ask user for inputs +#Initialize variables to False +ExpTCimport = 0 +SimTCImport = 0 +VizImport = 0 +SimExpPlot = 0 +VizXSection=0 +ExpTCGraph = 0 +crossrowsNcols = 0 +TdiffExpVSsim = 0 +dE_ExpVSsim = 0 +PerfMetr = 0 +PerfMetrGrTCs = 0 +Save_graphs=0 +plot_individ_SS_graphs = 0 +plot_multiple_on_same = 0 + +#EXPERIMENTAL DATA + +ExpTCGraph = int(input("Only graph experimental TC data? (1/0)")) + + +#VIZ DATA +SimExpPlot = int(input("Import sim Sand TC data and compare to exp data (avg of each level)? (1/0)")) +if SimExpPlot == 1: + SimTCImport = 1 +VizXSection = int(input("Plot temperature cross-sections slices of Viz data? (1/0)")) +if VizXSection == 1: + VizImport = 1 + +#CROSS-ROWS AND COLUMN PLOTS +crossrowsNcols = int(input("Plot cross-rows and cross-columns of sim and exp data for the sand store? (1/0)")) +if crossrowsNcols == 1: + SimTCImport = 1 + +#CALCULATE TEMPERATURE DIFFERENCES B/T EXPERIMENTAL AND SIMULATED DATA +TdiffExpVSsim = int(input("Calc. RMSD b/t sim and expmtl temperatures in sand store and plot entire_SS graph?")) + +if TdiffExpVSsim == 1: + #Plot each entire_SS individually? + plot_individ_SS_graphs = int(input("Plot each entire_SS graph individually?")) + + #Plot multiple graphs of entire_SS on same plot? + plot_multiple_on_same = int(input("Plot multiple graphs of entire_SS on same plot?")) + if plot_multiple_on_same == 1: + fig5, ax5 = plt.subplots() + labels5 = [None]*len(runs) + for run in range(0,len(runs)): + labels5[run] = input("What is the legend label for Run " + runs[run] + "?") + +Save_graphs=int(input("Save output graphs? (1/0)")) + +#CALCULATE DIFFERENCE IN TOTAL ENERGY STORED IN SAND +dE_ExpVSsim = int(input("Calc. KM3 metric?")) + +#PERFORMANCE METRICS +# PerfMetr = int(input("Import Performance Metrics data (ground TCs,etc)? (1/0)")) +# if PerfMetr ==1: +# PerfMetrGrTCs = int(input("Plot ground TCs against simulated ground temp data? (1/0)")) +# if PerfMetrGrTCs ==1: +# VizImport = 1 + +for run in range(0,len(runs)): + #%%% EXPERIMENTAL FILENAMES with import data + + #check the period to get the experimental file names accordingly + if diff_runs.Period[int(runs[run])].rstrip() == 'cooldown June 2021': + + #%%%% Cooldown + mode = 1 + #sand store temperatures file + filename1 = 'Sand_temperatures_2021-05-28-to-2021-07-06.dat' + #loops heat transfer file + filename2 ='Loops_heat_transfer_2021-05-28-to-2021-07-06.dat' + + elif diff_runs.Period[int(runs[run])].rstrip() == 'charging Mar 2021': + + #%%%% Charging + mode = 2 + #sand store temperatures file + filename1 = 'Sand_temperatures_2021-03-21-to-2021-03-23.dat' + #loops heat transfer file + filename2 ='Loops_heat_transfer_2021-03-21-to-2021-03-23.dat' + + elif diff_runs.Period[int(runs[run])].rstrip() == 'discharging Nov 2020': + + #%%%% Discharging + mode = 3 + #sand store temperatures file + filename1 = 'Sand_temperatures_2020-11-25-to-2020-11-27.dat' + #loops heat transfer file + filename2 ='Loops_heat_transfer_2020-11-25-to-2020-11-27.dat' + + elif diff_runs.Period[int(runs[run])].rstrip() == 'charging + discharging Nov 29': + + #%%%% Charging and Discharging (mode 4) + mode = 4 + #sand store temperatures file + filename1 = 'Sand_temperatures_2020-11-28-to-2020-11-30.dat' + #loops heat transfer file + filename2 ='Loops_heat_transfer_2020-11-28-to-2020-11-30.dat' + + elif diff_runs.Period[int(runs[run])].rstrip() == 'charging + discharging Nov 23': + + #%%%% Charging and Discharging 2 (mode 5) + mode = 5 + #sand store temperatures file + filename1 = 'Sand_temperatures_2020-11-23-to-2020-11-23.dat' + #loops heat transfer file + filename2 ='Loops_heat_transfer_2020-11-23-to-2020-11-23.dat' + elif diff_runs.Period[int(runs[run])].rstrip() == 'heating season 2021': + + #%%%% Entire heating season (mode 6) + mode = 6 + #sand store temperatures file + filename1 = 'Sand_temperatures_2020-11-19-to-2021-04-30.dat' + #loops heat transfer file + filename2 ='Loops_heat_transfer_2020-11-19-to-2021-04-30.dat' + + #%% 3. GRAPH EXPERIMENTAL TC DATA + if ExpTCGraph == 1: + + #%%% Import data + Exp_sand_temps, exp_start_date, exp_end_date = readSerranoSandMultipleDates(filename1) + + #%%% Process Experimental Sand Temp data into Numpy arrays + Exp_sand_temps_4d, exp_time_array = processExpTCDataToArrays(Exp_sand_temps,exp_start_date) + + #%%% Plot all layers of TCs over time, in sub plots + plotLayersofExpSandTempData_Sbplts(exp_time_array,Exp_sand_temps) + + #%%%% Plot layers of TCs over time, in individual plots + plotTopLayerofExpSandData_IndPlts(exp_time_array,Exp_sand_temps) + plotMiddleLayerofExpSandData_IndPlts(exp_time_array,Exp_sand_temps) + plotBottomLayerofExpSandData_IndPlts(exp_time_array,Exp_sand_temps) + + #%%%% Plot average of TCs in north/south rows over time, in individual plots, for each layer + plotExpSandDataAvgofNSRows_IndPlts(exp_time_array,Exp_sand_temps_4d) + # Plot average of TCs in west/east rows over time, in individual plots, for each layer + plotExpSandDataAvgofEWRows_IndPlts(exp_time_array,Exp_sand_temps_4d) + + #%%% Plot expmtl data cross-sections + # plotExpSandDataXSections(Exp_sand_temps_4d,0,60) + # plotExpSandDataYSections(Exp_sand_temps_4d,0,60) + + #%%% Plot average temps of sand store data (exp.) only + + #take averages of sand temps (lowest index is bottom-most layer of sand store) + avg_expmtl_temps = np.mean(np.mean(Exp_sand_temps_4d,axis=1),axis=1) + + #plot + fig1, ax1 = plt.subplots() + for col in range(0,np.size(avg_expmtl_temps,1)): + ax1.plot(exp_time_array,avg_expmtl_temps[:,col],marker=',') + + ax1.set_title("Mean temperatures of sand store layers") + ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + + ax1.legend(['Exp-bott','Exp-mid','Exp-top',],loc='best') + + + #Documentation on date locators + #https://matplotlib.org/stable/api/dates_api.html + locator = mdates.AutoDateLocator() + ax1.xaxis.set_major_locator(locator) + ax1.xaxis.set_minor_locator(mdates.DayLocator()) + + loc = ticker.MultipleLocator(base=2.0) + ax1.yaxis.set_major_locator(loc) + ax1.yaxis.set_minor_locator(ticker.MultipleLocator(base=1.0)) + + # #set the major locator on x-axis to every 6 days + # locator.intervald[mdates.DAILY]=[6] #https://matplotlib.org/stable/api/dates_api.html#matplotlib.dates.AutoDateLocator + + #Formatting - make the x-axis labels more consices + #https://matplotlib.org/stable/gallery/ticks_and_spines/date_concise_formatter.html + formatter = mdates.ConciseDateFormatter(locator) + ax1.xaxis.set_major_formatter(formatter) + fig1.autofmt_xdate() #automatically makes the x-labels rotate + + #gridlines on + ax1.grid(True) + #Set max and min values for axes + ax1.set_xlim(exp_time_array[0], exp_time_array[-1]) + if mode == 6: + ax1.set_ylim(20,60) + else: + ax1.set_ylim(30,60) + + + + #%% 4. SIMULATION TC DATA + #%%% Import simulation data into variables + + if SimTCImport == 1: + + #Import sand temp data + sndTemps,sndTemps_time, start_date, end_date, run_number = readOutputSimuDatafile(df.loc[run,'fileID']) + + + + + #%%% Process Simulated Sand Temp data into Numpy arrays + avg_sim_temps_snd_layer = np.mean(np.mean(sndTemps,axis=1),axis=1) + + #take average of all sand temps per timestep + avg_sim_temps_snd = np.mean(np.mean(np.mean(sndTemps,axis=1),axis=1),axis=1) + + + #%%%Print final avg Sand store temps to screen + print("Final avg sim temperatures in sand store layers for Run", run_number, ":\n", avg_sim_temps_snd_layer[-1,:]) + print("Bottom Middle Top") + + #%%% Make list of dates for x-axis of simulation plots + sndTemps_datetimes = DateArrayFromSimDayHoursArray(sndTemps_time[:,0], start_date) + sndTemps_datetimes[1] += timedelta(hours = sndTemps_time[1,1]) #the offset of 0.001h needs to be put back into the array -> it was lost with thw function DateArrayFromSimDayHoursArray + + + # %% 5. EXPERIMENTAL TC DATA + #%%%Import Experimental Data Files + if (ExpTCGraph != 1): + #it's the first run or the experimental data needed isn't the same as the last run + if (run == 0) or (diff_runs.Period[int(runs[run-1])] != diff_runs.Period[int(runs[run])]): + ##IMPORT MULTIPLE DATES OF SERRANO DATA + # import data + Exp_sand_temps, exp_start_date, exp_end_date = readSerranoSandMultipleDates( + filename1) + + # #IMPORT SINGLE DATE OF SERRANO DATA + # filename1 = 'Sand_temperatures_2021-06-04.dat' + # # import data + # Exp_sand_temps, exp_start_date = readSerranoSandSingleDate(filename1) + + # Process Experimental Sand Temp data into Numpy arrays + Exp_sand_temps_4d,exp_time_array = processExpTCDataToArrays(Exp_sand_temps,exp_start_date) + + + #take averages of sand temps (lowest index is bottom-most layer of sand store) + avg_expmtl_temps = np.mean(np.mean(Exp_sand_temps_4d,axis=1),axis=1) + + + if SimExpPlot == 1: + #%% 6. PLOT SIM & EXPT'L TC DATA + #%%% Plot simulated and experimental data + + fig1, ax1 = plt.subplots() + colors=['g','blue','r'] + + # Plot experimental data + labels = ['Exp Level C','Exp Level B','Exp Level A'] + + for col in range(np.size(avg_expmtl_temps,1)-1,-1,-1): + ax1.plot(exp_time_array,avg_expmtl_temps[:,col],color=colors[col], label=labels[col]) + + #Plot simulated data on same figure + labels = ['Sim Level C','Sim Level B','Sim Level A'] + #less than 115 has rewinds, so we don't want to graph them...start at line 964 in TC sand temp data + if int(run_number) < 115: + fig1, ax1 = plotTimeSeriesEntireArray(sndTemps_datetimes[964:], + avg_sim_temps_snd_layer[964:], + fig1, + ax1, + labels, + colors) + else: + fig1, ax1 = plotTimeSeriesEntireArray(sndTemps_datetimes, + avg_sim_temps_snd_layer, + fig1, + ax1, + labels, + colors) + + + #plt.plot(sndTemps_datetimes,avg_sim_temps_snd_layer[:,2],marker='*') + #ax1.set_xlabel("Days since start of simulation") + ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + + #rename figure and change legend to include experimental data + ax1.set_title("Mean temperatures of sand store layers, Run " + run_number) + ax1.legend() + + + #Documentation on date locators + #https://matplotlib.org/stable/api/dates_api.html + locator = mdates.AutoDateLocator() + ax1.xaxis.set_major_locator(locator) + ax1.xaxis.set_minor_locator(mdates.DayLocator()) + + loc = ticker.MultipleLocator(base=2.0) + ax1.yaxis.set_major_locator(loc) + ax1.yaxis.set_minor_locator(ticker.MultipleLocator(base=1.0)) + + # #set the major locator on x-axis to every 6 days + # locator.intervald[mdates.DAILY]=[6] #https://matplotlib.org/stable/api/dates_api.html#matplotlib.dates.AutoDateLocator + + #Formatting - make the x-axis labels more consices + #https://matplotlib.org/stable/gallery/ticks_and_spines/date_concise_formatter.html + formatter = mdates.ConciseDateFormatter(locator) + ax1.xaxis.set_major_formatter(formatter) + fig1.autofmt_xdate() #automatically makes the x-labels rotate + + #gridlines on + ax1.grid(True) + + #Set max and min values for axes + ax1.set_xlim(exp_time_array[0], exp_time_array[-1]) + if mode == 6: + ax1.set_ylim(20,58) + else: + ax1.set_ylim(30,58) + + if Save_graphs == 1: + #save graph to file with proper naming scheme + plt.savefig(graph_dir.joinpath("r" + df.loc[run,'runs'] + ".png"), dpi=150, bbox_inches='tight') + + #%%% Plot expmtl sand TC data (top layer) + + # #SUBPLOTS + # plotLayersofExpSandTempData_Sbplts(exp_time_array,Exp_sand_temps) + + # #INDIVIDUAL PLOTS + # plotTopLayerofExpSandData_IndPlts(exp_time_array,Exp_sand_temps) + + + + + + #%% 7. IMPORT VISUALIZATION DATA + if VizImport == 1: + + #%%% Import sand temp data + sndVizTemps,sndVizTemps_time, start_date_viz, end_date_viz, xcoords, x_viz_indices, ycoords, y_viz_indices, zcoords, z_viz_indices, TCx_coords, TCy_coords, simu_days, run_number_ID2 = readOutputSimuVizDatafile(df.loc[run,'fileID2']) + + + #%%% Convert arrays of strings to integers + x_viz_indices = x_viz_indices.astype(int) + y_viz_indices = y_viz_indices.astype(int) + z_viz_indices = z_viz_indices.astype(int) + + xcoords = xcoords.astype(float) + ycoords = ycoords.astype(float) + zcoords = zcoords.astype(float) + TCx_coords = TCx_coords.astype(float) + TCy_coords = TCy_coords.astype(float) + + #%%% Get Cartesian coords of indices that will be plotted + x_coords2plt = np.zeros(shape=(x_viz_indices.shape[0],1)) #initialize array + for i in range(0,x_coords2plt.shape[0]): + x_coords2plt[i] = xcoords[x_viz_indices[i]] + + y_coords2plt = np.zeros(shape=(y_viz_indices.shape[0],1)) #initialize array + for j in range(0,y_coords2plt.shape[0]): + y_coords2plt[j] = ycoords[y_viz_indices[j]] + + z_coords2plt = np.zeros(shape=(z_viz_indices.shape[0],1)) #initialize array + for k in range(0,z_coords2plt.shape[0]): + z_coords2plt[k] = zcoords[z_viz_indices[k]] + + #%% 8. PLOT DOMAIN VISUALIZATION DATA + + if VizXSection == 1: + + # Plot temperature cross-sections slices + + # %%%% Fixed variables for plotting + # Index values for the 3 axes to plot temps at. + aa=round(len(x_viz_indices)/2) + cc=round(len(z_viz_indices)/2) + #due to simulations with N/S symmetry about y-axis b/t runs 63 and 120 (inclusive) + if int(run_number_ID2) >=63 and int(run_number_ID2) <121: + bb=round(len(y_viz_indices))-1 + else: + bb=round(len(y_viz_indices)/2) + + #data points will be plotted every _(data_pt_freq)__ days + data_pt_freq = 10 + + #x-, y-, and z-coordinates for concrete + x_conc_coord1 = np.mean([x_coords2plt[5,0],x_coords2plt[6,0]]) + x_conc_coord2 = np.mean([x_coords2plt[7,0],x_coords2plt[8,0]]) + x_conc_coord3 = np.mean([x_coords2plt[25,0],x_coords2plt[26,0]]) + x_conc_coord4 = np.mean([x_coords2plt[27,0],x_coords2plt[28,0]]) + + y_conc_coord1 = np.mean([x_coords2plt[5,0],x_coords2plt[6,0]]) + y_conc_coord2 = np.mean([x_coords2plt[7,0],x_coords2plt[8,0]]) + y_conc_coord3 = np.mean([x_coords2plt[25,0],x_coords2plt[26,0]]) + y_conc_coord4 = np.mean([x_coords2plt[27,0],x_coords2plt[28,0]]) + + z_conc_coord1 = np.mean([z_coords2plt[5,0],z_coords2plt[6,0]]) + z_conc_coord2 = np.mean([z_coords2plt[7,0],z_coords2plt[8,0]]) + + #x-, y-, and z-coordinates for XPS + x_XPS_coord1 = np.mean([x_coords2plt[7,0],x_coords2plt[8,0]]) + x_XPS_coord2 = np.mean([x_coords2plt[9,0],x_coords2plt[10,0]]) + x_XPS_coord3 = np.mean([x_coords2plt[23,0],x_coords2plt[24,0]]) + x_XPS_coord4 = np.mean([x_coords2plt[25,0],x_coords2plt[26,0]]) + + y_XPS_coord1 = np.mean([x_coords2plt[7,0],x_coords2plt[8,0]]) + y_XPS_coord2 = np.mean([x_coords2plt[9,0],x_coords2plt[10,0]]) + y_XPS_coord3 = np.mean([x_coords2plt[23,0],x_coords2plt[24,0]]) + y_XPS_coord4 = np.mean([x_coords2plt[25,0],x_coords2plt[26,0]]) + + z_XPS_coord1 = np.mean([z_coords2plt[7,0],z_coords2plt[8,0]]) + z_XPS_coord2 = np.mean([z_coords2plt[9,0],z_coords2plt[10,0]]) + z_XPS_coord3 = np.mean([z_coords2plt[21,0],z_coords2plt[22,0]]) + z_XPS_coord4 = np.mean([z_coords2plt[23,0],z_coords2plt[24,0]]) + + + # %%%% Plot across z, one snapshot in time + # plot one temp at x=0,y=16, vertical cross-section + # fig2,ax2 = plt.subplots() + # ax2.plot(x_coords2plt,sndVizTemps[0,16,:,13]) + + + # %%%% Plot multiple graphs across x-coords over time + # aa=round(len(x_viz_indices)/2) + # cc=round(len(z_viz_indices)/2) + # # due to N/S symmetry about y-axis: + # if int(run_number_ID2) >=63 and int(run_number_ID2) <121: + # bb=round(len(y_viz_indices))-1 + # else: + # bb=round(len(y_viz_indices)/2) + + # #plot temperatures at ~mid-point in z-direction, across all x-coords, and at the ~mid-point of y + # #plot hourly figures at each time step + # for time in range(0,sndVizTemps_time.shape[0]): + # fig2,ax2 = plt.subplots() + # ax2.plot(x_coords2plt[:,0],sndVizTemps[time,:,bb,cc],'+') + # #plot node coordinate spacing on figure + # ax2.plot(x_coords2plt,np.ones(shape=(x_coords2plt.shape),dtype=int)*80,'+') + # ax2.set_title("Run "+run_number_ID2+ ": Temps at t=" + "{:.2f}".format(sndVizTemps_time[time,0]) + "days across x-axis, at y=" + str(y_viz_indices[bb]) + " and z=" + str(z_viz_indices[cc])) + # ax2.set_xlabel('Length of num. domain (m)') + # ax2.set_ylabel('Temperature (\N{DEGREE SIGN}C)') + # ax2.set_ylim(0,70) + # #plot vertical shading for XPS and conc + # ax2.axvspan(x_conc_coord1,x_conc_coord2,color='grey', alpha=0.5, lw=0) + # ax2.axvspan(x_conc_coord3,x_conc_coord4,color='grey', alpha=0.5, lw=0) + # ax2.axvspan(x_XPS_coord1,x_XPS_coord2,color='pink', alpha=0.5, lw=0) + # ax2.axvspan(x_XPS_coord3,x_XPS_coord4,color='pink', alpha=0.5, lw=0) + + # %%%% Plot multiple graphs across y-coords over time + # # plot temperatures at ~mid-point in x-direction, across all y-coords, and at the ~mid-point of z + # #plot hourly figures at each time step + # aa=19 + # cc=22 + # for time in range(0,sndVizTemps_time.shape[0]): + # fig2,ax2 = plt.subplots() + # ax2.plot(y_coords2plt[:,0],sndVizTemps[time,aa,:,cc]) + # ax2.set_title("Run "+run_number_ID2+ ": Temps at t=" + "{:.2f}".format(sndVizTemps_time[time,0]) + "days across y-axis, at x=" + str(x_viz_indices[aa]) + " and z=" + str(z_viz_indices[cc])) + # ax2.set_xlabel('Length of num. domain (m)') + # ax2.set_ylabel('Temperature (\N{DEGREE SIGN}C)') + # ax2.set_ylim(0,70) + # #plot vertical shading for XPS and conc + # ax2.axvspan(x_conc_coord1,x_conc_coord2,color='grey', alpha=0.5, lw=0) + # ax2.axvspan(x_conc_coord3,x_conc_coord4,color='grey', alpha=0.5, lw=0) + # ax2.axvspan(x_XPS_coord1,x_XPS_coord2,color='pink', alpha=0.5, lw=0) + # ax2.axvspan(x_XPS_coord3,x_XPS_coord4,color='pink', alpha=0.5, lw=0) + + + # %%%% Plot multiple graphs across z-coords over time + # #plot temperatures at ~mid-point in x-direction, across all z-coords, and at the ~mid-point of y + # #plot hourly figures at each time step + # for time in range(0,sndVizTemps_time.shape[0]): + # fig2,ax2 = plt.subplots() + # ax2.plot(z_coords2plt[:,0],sndVizTemps[time,aa,bb,:]) + # ax2.set_title("Run "+run_number_ID2+ ": Temps at t=" + "{:.2f}".format(sndVizTemps_time[time,0]) + "days across z-axis, at x=" + str(x_viz_indices[aa]) + " and y=" + str(y_viz_indices[bb])) + # ax2.set_xlabel('Length of num. domain (m)') + # ax2.set_ylabel('Temperature (\N{DEGREE SIGN}C)') + # ax2.set_ylim(0,60) + + + # #plot node coordinate spacing on figure + # ax2.plot(z_coords2plt,np.ones(shape=(z_coords2plt.shape),dtype=int)*50,'+') + + # #plot vertical shading for XPS and conc + # ax2.axvspan(z_conc_coord1,z_conc_coord2,color='grey', alpha=0.5, lw=0) + # ax2.axvspan(z_XPS_coord1,z_XPS_coord2,color='pink', alpha=0.5, lw=0) + # ax2.axvspan(z_XPS_coord3,z_XPS_coord4,color='pink', alpha=0.5, lw=0) + + + # %%%% Plot one temp distribution per day across x-axis + fig1,ax1 = plt.subplots() + + #rewinds no longer being used at run 115 and above + if int(run_number_ID2) >= 115: + for time in range(0,sndVizTemps_time.shape[0],data_pt_freq): + print("Graphed point (day,hour): " + str(sndVizTemps_time[time])) + ax1.plot(x_coords2plt[:,0],sndVizTemps[time,:,bb,cc]) + #plot second temp distribution in array, steady state temp distribution after IC + ax1.plot(x_coords2plt[:,0],sndVizTemps[1,:,bb,cc]) + print("Graphed point (day,hour): " + str(sndVizTemps_time[1])) + else: + for time in range(0,sndVizTemps_time.shape[0],data_pt_freq): + print("Graphed point (day,hour): " + str(sndVizTemps_time[time])) + ax1.plot(x_coords2plt[:,0],sndVizTemps[time,:,bb,cc]) + + #plot last temp distribution in array, just before midnight + ax1.plot(x_coords2plt[:,0],sndVizTemps[-1,:,bb,cc]) + print("Graphed point (day,hour): " + str(sndVizTemps_time[-1])) + # #plot node coordinate spacing on figure + # ax1.plot(x_coords2plt,np.ones(shape=(x_coords2plt.shape),dtype=int)*50,'+') + + #plot vertical shading for XPS and conc + ax1.axvspan(x_conc_coord1,x_conc_coord2,color='grey', alpha=0.5, lw=0) + ax1.axvspan(x_conc_coord3,x_conc_coord4,color='grey', alpha=0.5, lw=0) + ax1.axvspan(x_XPS_coord1,x_XPS_coord2,color='pink', alpha=0.5, lw=0) + ax1.axvspan(x_XPS_coord3,x_XPS_coord4,color='pink', alpha=0.5, lw=0) + + ax1.set_title("Run "+run_number_ID2+ ": Daily temp distribution for " + str(round(sndVizTemps_time[-1,0])) + "days across x-axis, at y=" + str(y_viz_indices[bb]) + " and z=" + str(z_viz_indices[cc])) + ax1.set_ylim(0,60) + ax1.set_xlabel('Length of num. domain (m)') + ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + + if Save_graphs == 1: + #save graph to file with proper naming scheme + plt.savefig(graph_dir.joinpath("r" + df.loc[run,'runs'] + "_viz_x.png"), dpi=150, bbox_inches='tight') + + + # %%%% Plot one temp distribution per day across y-axis + + fig1,ax1 = plt.subplots() + + #rewinds no longer being used at run 115 and above, domain symmetry still used + if int(run_number_ID2) >= 115 and int(run_number_ID2) < 121: + for time in range(0,sndVizTemps_time.shape[0],data_pt_freq): + print("Graphed point (day,hour): " + str(sndVizTemps_time[time])) + ax1.plot(x_coords2plt[:,0], + np.concatenate((sndVizTemps[time,aa,:,cc],np.flip(sndVizTemps[time,aa,:,cc],0)),axis=0)) + #plot second temp distribution in array, steady state temp distribution after IC + ax1.plot(x_coords2plt[:,0], + np.concatenate((sndVizTemps[1,aa,:,cc],np.flip(sndVizTemps[1,aa,:,cc],0)),axis=0)) + print("Graphed point (day,hour): " + str(sndVizTemps_time[1])) + + #rewinds no longer being used at run 115 (and 121) and above, domain symmetry no longer used + elif int(run_number_ID2) >= 121: + for time in range(0,sndVizTemps_time.shape[0],data_pt_freq): #,(precon_days+1)): + print("Graphed point (day,hour): " + str(sndVizTemps_time[time])) + ax1.plot(y_coords2plt[:,0],sndVizTemps[time,aa,:,cc]) + #plot second temp distribution in array, steady state temp distribution after IC + ax1.plot(y_coords2plt[:,0],sndVizTemps[1,aa,:,cc]) + print("Graphed point (day,hour): " + str(sndVizTemps_time[1])) + + + #rewinds being used, domain symmetry used at run 63 and above + elif int(run_number_ID2) >=63: + for time in range(0,sndVizTemps_time.shape[0],data_pt_freq): #,(precon_days+1)): + print("Graphed point (day,hour): " + str(sndVizTemps_time[time])) + ax1.plot(x_coords2plt[:,0], + np.concatenate((sndVizTemps[time,aa,:,cc],np.flip(sndVizTemps[time,aa,:,cc],0)),axis=0)) + #rewinds being used, domain symmetry not used at run 63 and below + else: + for time in range(0,sndVizTemps_time.shape[0],data_pt_freq): #,(precon_days+1)): + print("Graphed point (day,hour): " + str(sndVizTemps_time[time])) + ax1.plot(y_coords2plt[:,0],sndVizTemps[time,aa,:,cc]) + + #plot last temp distribution in array, just before midnight + if int(run_number_ID2) >=63 and int(run_number_ID2) <121: + #due to simulations with N/S symmetry about y-axis starting at run 63 and above + ax1.plot(x_coords2plt[:,0], + np.concatenate((sndVizTemps[-1,aa,:,cc],np.flip(sndVizTemps[-1,aa,:,cc],0)),axis=0)) + print("Graphed point (day,hour): " + str(sndVizTemps_time[-1])) + # #plot node coordinate spacing on figure + # ax1.plot(x_coords2plt,np.ones(shape=(x_coords2plt.shape),dtype=int)*50,'+') + else: + #no N/S symmetry + ax1.plot(y_coords2plt[:,0],sndVizTemps[-1,aa,:,cc]) + print("Graphed point (day,hour): " + str(sndVizTemps_time[-1])) + # #plot node coordinate spacing on figure + # ax1.plot(y_coords2plt,np.ones(shape=(y_coords2plt.shape),dtype=int)*50,'+') + + #plot vertical shading for XPS and conc + ax1.axvspan(x_conc_coord1,x_conc_coord2,color='grey', alpha=0.5, lw=0) + ax1.axvspan(x_conc_coord3,x_conc_coord4,color='grey', alpha=0.5, lw=0) + ax1.axvspan(x_XPS_coord1,x_XPS_coord2,color='pink', alpha=0.5, lw=0) + ax1.axvspan(x_XPS_coord3,x_XPS_coord4,color='pink', alpha=0.5, lw=0) + + ax1.set_title("Run "+run_number_ID2+ ": Daily temp distribution for " + str(round(sndVizTemps_time[-1,0])) + "days across y-axis, at x=" + str(x_viz_indices[aa]) + " and z=" + str(z_viz_indices[cc])) + ax1.set_xlabel('Length of num. domain (m)') + ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + ax1.set_ylim(0,60) + + if Save_graphs == 1: + #save graph to file with proper naming scheme + plt.savefig(graph_dir.joinpath("r" + df.loc[run,'runs'] + "_viz_y.png"), dpi=150, bbox_inches='tight') + + # %%%% Plot one temp distribution per day across z-axis + + fig1,ax1 = plt.subplots() + + #rewinds no longer being used at run 115 and above + if int(run_number_ID2) >= 115: + for time in range(0,sndVizTemps_time.shape[0],data_pt_freq): + print("Graphed point (day,hour): " + str(sndVizTemps_time[time])) + ax1.plot(z_coords2plt[:,0],sndVizTemps[time,aa,bb,:]) + #plot second temp distribution in array, steady state temp distribution after IC + ax1.plot(z_coords2plt[:,0],sndVizTemps[1,aa,bb,:]) + print("Graphed point (day,hour): " + str(sndVizTemps_time[1])) + else: + + for time in range(0,sndVizTemps_time.shape[0],data_pt_freq): + print("Graphed point (day,hour): " + str(sndVizTemps_time[time])) + ax1.plot(z_coords2plt[:,0],sndVizTemps[time,aa,bb,:]) + + + #plot last temp distribution in array, just before midnight + ax1.plot(z_coords2plt[:,0],sndVizTemps[-1,aa,bb,:]) + print("Graphed point (day,hour): " + str(sndVizTemps_time[-1])) + # #plot node coordinate spacing on figure + # ax1.plot(z_coords2plt,np.ones(shape=(z_coords2plt.shape),dtype=int)*50,'+') + + #plot vertical shading for XPS and conc + ax1.axvspan(z_conc_coord1,z_conc_coord2,color='grey', alpha=0.5, lw=0) + ax1.axvspan(z_XPS_coord1,z_XPS_coord2,color='pink', alpha=0.5, lw=0) + ax1.axvspan(z_XPS_coord3,z_XPS_coord4,color='pink', alpha=0.5, lw=0) + + + ax1.set_title("Run "+run_number_ID2+ ": Daily temp distribution for " + str(round(sndVizTemps_time[-1,0])) + "days across z-axis, at x=" + str(x_viz_indices[aa]) + " and y=" + str(y_viz_indices[bb])) + ax1.set_xlabel('Length of num. domain (m)') + ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + ax1.set_ylim(0,60) + + if Save_graphs == 1: + #save graph to file with proper naming scheme + plt.savefig(graph_dir.joinpath("r" + df.loc[run,'runs'] + "_viz_z.png"), dpi=150, bbox_inches='tight') + + #%% 8a. PLOT VIZ DATA AT SPECIFIED TIMESTEP -> A=sndVizTemps[timestep,:,:,:] + #Plot all node temps in y-z cross sectional plane at time = timestep + + # timestep = 0 + # A=sndVizTemps[timestep,:,:,:] + + #%%%% ONE LINE PER GRAPH, MANY GRAPHS + # # To graph N->S temps along x-axis of sand store, at a fixed z-value, going through all y-points + # z= 31 #z_viz_indices index value + # for y in range(0,A.shape[1]): + # fig1,ax1 = plt.subplots() + # ax1.plot(x_coords2plt[:,0],A[:,y,z]) + # ax1.set_title("N->S temps along x-axis, at Viz timestep " +str(timestep)+ ", at y=" + str(y_viz_indices[y]) +", z="+ str(z_viz_indices[z])) + # ax1.set_xlabel('Length of num. domain in x-direction (m)') + # ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + # ax1.set_ylim(0,80) + + # #To graph W->E temps along y-axis of sand store, at a fixed x-value, going up through all z-points + # x=30 + # for z in range(0,A.shape[2]): + # fig1,ax1 = plt.subplots() + # ax1.plot(y_coords2plt[:,0],A[x,:,z]) + # ax1.set_title("W->E temps along y-axis, at Viz timestep " +str(timestep)+ ", at x=" + str(x_coords2plt[x,0]) +"m, z="+ str(z_coords2plt[z,0]) + "m") + # ax1.set_xlabel('Length of num. domain in y-direction(m)') + # ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + # ax1.set_ylim(0,80) + # #To graph N->S temps along x-axis of sand store, at a fixed y-value, going up through all z-points + # y=30 + # for z in range(0,A.shape[2]): + # fig1,ax1 = plt.subplots() + # ax1.plot(x_coords2plt[:,0],A[:,y,z]) + # ax1.set_title("N->S temps along x-axis, at Viz timestep " +str(timestep)+ ", at y=" + str(y_coords2plt[y,0]) +"m, z="+ str(z_coords2plt[z,0]) + "m") + # ax1.set_xlabel('Length of num. domain in x-direction (m)') + # ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + # ax1.set_ylim(0,80) + + #%%%% MANY LINES PER GRAPH, ONE GRAPH + # #To graph N->S temps along x-axis of sand store, for all y-nodes, at a fixed z-value + # z= 24 #z_viz_indices index value + # fig1,ax1 = plt.subplots() + # for y in range(0,A.shape[1]): + # ax1.plot(x_coords2plt[:,0],A[:,y,z]) + # ax1.set_title("N->S temps (x-axis), at Viz timestep " +str(timestep)+ ", at each y output node, at z-node "+ str(z_viz_indices[z])) + # ax1.set_xlabel('Length of num. domain in x-direction (m)') + # ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + # ax1.set_ylim(0,80) + #%%%%MANY LINES PER GRAPH, ONE GRAPH -> used to troubleshoot while running solver + # #To graph N->S temps along x-axis of sand store, for all y-nodes, at a fixed z-value + # A=Temps[:,:,:,0] + # z= 24 #z_viz_indices index value + # fig1,ax1 = plt.subplots() + # for y in range(0,A.shape[1]): + # ax1.plot(node_xcoords,A[:,y,z]) + # ax1.set_title("N->S temps (x-axis), at each y output node, at z-node "+ str(z)) + # ax1.set_xlabel('Length of num. domain in x-direction (m)') + # ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + # ax1.set_ylim(0,80) + + #%%%%MANY LINES PER GRAPH, MANY GRAPHS -> used to troubleshoot while running solver + #To graph N->S temps along x-axis of sand store, for all y-nodes, at each z-value + # A=Temps[:,:,:,0] + # for z in range(0,A.shape[2]): + # fig1,ax1 = plt.subplots() + # for y in range(0,A.shape[1]): + # ax1.plot(node_xcoords,A[:,y,z]) + # ax1.set_title("N->S temps (x-axis), at each y output node, at z-node "+ str(z)) + # ax1.set_xlabel('Length of num. domain in x-direction (m)') + # ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + # ax1.set_ylim(0,80) + + + # for z in range(0,A.shape[2]): + # fig1,ax1 = plt.subplots() + # for x in range(0,A.shape[0]): + # ax1.plot(node_ycoords,A[x,:,z]) + # ax1.set_title("E->W temps (y-axis), at each x output node, at z-node "+ str(z)) + # ax1.set_xlabel('Length of num. domain in y-direction (m)') + # ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + # ax1.set_ylim(0,80) + + #%%%%MANY LINES PER GRAPH, MANY GRAPHS + # #To graph N->S temps along x-axis of sand store, for all y-nodes, at each z-value + # for z in range(0,A.shape[2]): + # fig1,ax1 = plt.subplots() + # for y in range(0,A.shape[1]): + # ax1.plot(x_coords2plt[:,0],A[:,y,z]) + # ax1.set_title("N->S temps (x-axis), at Viz timestep " +str(timestep)+ ", at each y output node, at z-node "+ str(z_viz_indices[z])) + # ax1.set_xlabel('Length of num. domain in x-direction (m)') + # ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + # ax1.set_ylim(0,80) + + + # #MANY LINES PER GRAPH, MANY GRAPHS + # #To graph W->E temps along y-axis of sand store, for all x-nodes, at each z-value + # for z in range(0,A.shape[2]): + # fig1,ax1 = plt.subplots() + # for x in range(0,A.shape[0]): + # ax1.plot(y_coords2plt[:,0],A[x,:,z]) + # ax1.set_title("W->E temps (y-axis), at Viz timestep " +str(timestep)+ ", at each x output node, at z-node "+ str(z_viz_indices[z])) + # ax1.set_xlabel('Length of num. domain in y-direction (m)') + # ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + # ax1.set_ylim(0,80) + + # #MANY LINES PER GRAPH, MANY GRAPHS + # #To graph B->T temps along z-axis of sand store, for all x-nodes, at each y-value + # for y in range(0,A.shape[1]): + # fig1,ax1 = plt.subplots() + # for x in range(0,A.shape[0]): + # ax1.plot(z_coords2plt[:,0],A[x,y,:]) + # ax1.set_title("B->T temps (z-axis), at Viz timestep " +str(timestep)+ ", at each x output node, at y-node "+ str(y_viz_indices[y])) + # ax1.set_xlabel('Height of num. domain in z-direction (m)') + # ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + # ax1.set_ylim(0,80) + + + + + #%% 9. PLOT SAND SIMU AND EXP. DATA N/S AND E/W AVERAGES + + '''this function plots 1 graph with 6 lines of data, each line is the average TC from each row of TCs (north to south) + It will compare the simulated data to the experimental data + ''' + if crossrowsNcols == 1: + TC_exp_EW_row_avgs = np.mean(Exp_sand_temps_4d,axis=2) + TC_sim_EW_row_avgs = np.mean(sndTemps,axis=2) + + TC_exp_NS_row_avgs = np.mean(Exp_sand_temps_4d,axis=1) + TC_sim_NS_row_avgs = np.mean(sndTemps,axis=1) + + layer_name = [", bottom", + ", middle", + ", top"] + + clrs = ['blue', + 'darkorange', + 'g', + 'r', + 'dodgerblue', + 'darkorchid'] + + #%%% Plot all sim and exp rows on one plot + + shift_EW = np.zeros(shape=[3,6]) + #because the averages of the simltd and expmtl rows did not start at the same temperatures, find the difference in starting temperatures b/t these arrays + for layer in range(2,-1,-1): + for i in range(0,6): + shift_EW[layer,i] = TC_exp_EW_row_avgs[0,i,layer] - TC_sim_EW_row_avgs[0,i,layer] + + for layer in range(2,-1,-1): + #E/W rows, from northmost to southmost + fig,ax = plt.subplots() + + + lgnd = ["R1, N",# + layer_name[layer], + "R2",#+ layer_name[layer], + "R3",#+ layer_name[layer], + "R4",#+ layer_name[layer], + "R5",#+ layer_name[layer], + "R6, S"]#+ layer_name[layer]] + + #plot + for i in range(0,6): + ax.plot(exp_time_array,TC_exp_EW_row_avgs[:,i,layer], label="Exp "+ lgnd[i], color=clrs[i]) #plot top layer of sand store + #plot the simulated rows with the shift between starting exp and sim avg temps + ax.plot(sndTemps_datetimes,TC_sim_EW_row_avgs[:,i,layer] + shift_EW[layer,i], label="Sim " + lgnd[i], linestyle='--',color=clrs[i]) #plot top layer of sand store + + ax.set_title("Sand store"+ layer_name[layer] +" layer, avg TC row temps, " + "Run " + str(runs[run]))#: "+ str(exp_start_date) + " to " + str(exp_end_date)) + ax.legend() + ax.set_ylim(20,60) + + #locator = mdates.DayLocator(interval=2) + locator = mdates.AutoDateLocator() + # locator = mdates.MonthLocator() + formatter = mdates.ConciseDateFormatter(locator) + ax.grid(True) + ax.set_ylabel("Temperature ($\degree$C)") + ax.xaxis.set_major_locator(locator) + # ax.xaxis.set_minor_locator(mdates.MonthLocator()) + ax.xaxis.set_major_formatter(formatter) + ax.set_ylim(20,60) + fig.autofmt_xdate() #automatically makes the x-labels rotate + + if Save_graphs == 1: + #save graph to file with proper naming scheme + plt.savefig(graph_dir.joinpath("r" + df.loc[run,'runs'] + "N2S" + layer_name[layer] + ".png"), dpi=150, bbox_inches='tight') + + + + + + shift_NS = np.zeros(shape=[3,6]) + #because the averages of the simltd and expmtl rows did not start at the same temperatures, find the difference in starting temperatures b/t these arrays + for layer in range(2,-1,-1): + for j in range(0,6): + shift_NS[layer,j] = TC_exp_NS_row_avgs[0,j,layer] - TC_sim_NS_row_avgs[0,j,layer] + + + + + #N/S rows, from westmost to eastmost + fig1,ax1 = plt.subplots() + lgnd = ["C1, W",# + layer_name[layer], + "C2",#+ layer_name[layer], + "C3",#+ layer_name[layer], + "C4",#+ layer_name[layer], + "C5",#+ layer_name[layer], + "C6, E"]#+ layer_name[layer]] + + #plot + for j in range(0,6): + ax1.plot(exp_time_array,TC_exp_NS_row_avgs[:,j,layer], label="Exp "+ lgnd[j], color=clrs[j]) #plot top layer of sand store + #plot the simulated rows with the shift between starting exp and sim avg temps + ax1.plot(sndTemps_datetimes,TC_sim_NS_row_avgs[:,j,layer] + shift_NS[layer,j], label="Sim " + lgnd[j], linestyle='--',color=clrs[j]) #plot top layer of sand store + + ax1.set_title("Sand store"+ layer_name[layer] +" layer, avg TC column temps, " + "Run " + str(runs[run]))#: "+ str(exp_start_date) + " to " + str(exp_end_date)) + ax1.legend() + ax1.set_ylim(20,60) + + ax1.grid(True) + ax1.set_ylabel("Temperature ($\degree$C)") + ax1.xaxis.set_major_locator(locator) + # ax.xaxis.set_minor_locator(mdates.MonthLocator()) + ax1.xaxis.set_major_formatter(formatter) + ax1.set_ylim(20,60) + fig1.autofmt_xdate() #automatically makes the x-labels rotate + + if Save_graphs == 1: + #save graph to file with proper naming scheme + plt.savefig(graph_dir.joinpath("r" + df.loc[run,'runs'] + "W2E" + layer_name[layer] + ".png"), dpi=150, bbox_inches='tight') + + #%%% Plot only simulated data, all rows in each layer on one plot + + # for layer in range(2,-1,-1): + # #E/W rows, from northmost to southmost + # fig,ax = plt.subplots() + # lgnd = ["Row 1, northern edge" + layer_name[layer], + # "Row 2"+ layer_name[layer], + # "Row 3"+ layer_name[layer], + # "Row 4"+ layer_name[layer], + # "Row 5"+ layer_name[layer], + # "Row 6, southern edge"+ layer_name[layer]] + # for j in range(0,6): + + # ax.plot(sndTemps_datetimes,TC_sim_EW_row_avgs[:,j,layer], label=lgnd[j], color=clrs[j]) #plot top layer of sand store + + # #N/S rows, from westmost to eastmost + # fig1,ax1 = plt.subplots() + # lgnd = ["Col 1, western edge" + layer_name[layer], + # "Col 2"+ layer_name[layer], + # "Col 3"+ layer_name[layer], + # "Col 4"+ layer_name[layer], + # "Col 5"+ layer_name[layer], + # "Col 6, eastern edge"+ layer_name[layer]] + # for j in range(0,6): + # ax1.plot(sndTemps_datetimes,TC_sim_NS_row_avgs[:,j,layer], label=lgnd[j],color=clrs[j]) #plot top layer of sand store + + # #locator = mdates.DayLocator(interval=2) + # locator = mdates.AutoDateLocator() + # # locator = mdates.MonthLocator() + # formatter = mdates.ConciseDateFormatter(locator) + + + # ax.legend() + # ax.grid(True) + # ax.set_title("Sand store"+ layer_name[layer] +" layer, avg sim TC row temps, " + "Run " + str(runs[run])) #: "+ str(exp_start_date) + " to " + str(exp_end_date)) + # ax.set_ylabel("Temperature ($\degree$C)") + # ax.xaxis.set_major_locator(locator) + # # ax.xaxis.set_minor_locator(mdates.MonthLocator()) + # ax.xaxis.set_major_formatter(formatter) + # ax.set_ylim(20,60) + # fig.autofmt_xdate() #automatically makes the x-labels rotate + + # ax1.legend() + # ax1.grid(True) + # ax1.set_title("Sand store"+ layer_name[layer] +" layer, avg sim TC col temps, " + "Run " + str(runs[run]))#: "+ str(exp_start_date) + " to " + str(exp_end_date)) + # ax1.set_ylabel("Temperature ($\degree$C)") + # ax1.xaxis.set_major_locator(locator) + # #ax1.xaxis.set_minor_locator(mdates.DayLocator()) + # ax1.xaxis.set_major_formatter(formatter) + # ax1.set_ylim(20,60) + # fig1.autofmt_xdate() #automatically makes the x-labels rotate + + + #%%% Plot each row separately (sim and viz) + + # #plot E/W north to south layers of sand temps + # for layer in range(2,-1,-1): + # lgnd = ["Row 1, northern edge" + layer_name[layer], + # "Row 2"+ layer_name[layer], + # "Row 3"+ layer_name[layer], + # "Row 4"+ layer_name[layer], + # "Row 5"+ layer_name[layer], + # "Row 6, southern edge"+ layer_name[layer]] + # for j in range(0,6): + # fig,ax = plt.subplots() + # ax.plot(exp_time_array,TC_exp_EW_row_avgs[:,j,layer], label="Exp",color=clrs[j]) #plot top layer of sand store + # ax.plot(sndTemps_datetimes,TC_sim_EW_row_avgs[:,j,layer], label="Sim", linestyle='--',color=clrs[j]) #plot top layer of sand store + + # #locator = mdates.DayLocator(interval=2) + # locator = mdates.AutoDateLocator() + # # locator = mdates.MonthLocator() + # formatter = mdates.ConciseDateFormatter(locator) + # # ax.legend([layer_name[layer] +" Row 1 (north)", + # # layer_name[layer] +" Row 2", + # # layer_name[layer] +" Row 3", + # # layer_name[layer] +" Row 4", + # # layer_name[layer] +" Row 5", + # # layer_name[layer] +" Row 6 (south)"]) + # ax.legend() + # ax.grid(True) + # ax.set_title("Sand "+ lgnd[j] + " layer, avg TC row temps, " + "Run " + str(runs[run]))#: "+ str(exp_start_date) + " to " + str(exp_end_date)) + # ax.set_ylabel("Temperature ($\degree$C)") + # ax.xaxis.set_major_locator(locator) + # # ax.xaxis.set_minor_locator(mdates.MonthLocator()) + # ax.xaxis.set_major_formatter(formatter) + # ax.set_ylim(20,60) + # fig.autofmt_xdate() #automatically makes the x-labels rotate + + # #N/S rows, from westmost to eastmost + # for layer in range(2,-1,-1): + # lgnd = ["Col 1, western edge" + layer_name[layer], + # "Col 2"+ layer_name[layer], + # "Col 3"+ layer_name[layer], + # "Col 4"+ layer_name[layer], + # "Col 5"+ layer_name[layer], + # "Col 6, eastern edge"+ layer_name[layer]] + # for j in range(0,6): + # fig,ax = plt.subplots() + # ax.plot(exp_time_array,TC_exp_NS_row_avgs[:,j,layer], label="Exp",color=clrs[j]) #plot top layer of sand store + # ax.plot(sndTemps_datetimes,TC_sim_NS_row_avgs[:,j,layer], label="Sim", linestyle='--',color=clrs[j]) #plot top layer of sand store + + # #locator = mdates.DayLocator(interval=2) + # locator = mdates.AutoDateLocator() + # # locator = mdates.MonthLocator() + # formatter = mdates.ConciseDateFormatter(locator) + # # ax.legend([layer_name[layer] +" Row 1 (north)", + # # layer_name[layer] +" Row 2", + # # layer_name[layer] +" Row 3", + # # layer_name[layer] +" Row 4", + # # layer_name[layer] +" Row 5", + # # layer_name[layer] +" Row 6 (south)"]) + # ax.legend() + # ax.grid(True) + # ax.set_title("Sand "+ lgnd[j] + " layer, avg TC col temps, " + "Run " + str(runs[run]))#: "+ str(exp_start_date) + " to " + str(exp_end_date)) + # ax.set_ylabel("Temperature ($\degree$C)") + # ax.xaxis.set_major_locator(locator) + # # ax.xaxis.set_minor_locator(mdates.MonthLocator()) + # ax.xaxis.set_major_formatter(formatter) + # ax.set_ylim(20,60) + # fig.autofmt_xdate() #automatically makes the x-labels rotate + +#%% 10. PLOTS & RMSD OF EXP vs SIM + + '''This will calculate the difference between temperature nodes of experimental vs simulated data, + to give overall metrics for "fit" of the model +''' + + #does this section need to be done? + if TdiffExpVSsim == 1: + ExpTCimport = 0 + + #%%% Check if data has been imported + #were sim TC temps already imported? + if (SimExpPlot == 1 or crossrowsNcols == 1): + pass + else: + + #Check if sim file exists + try: + open(df.loc[run,'fileID'], 'r') + except FileNotFoundError: #if it doesn't exist, print error and move on to next run + print('**************************') + print("File " + df.loc[run,'fileID'] + " does not exist!!!") + print('**************************') + #skip the rest of that run + continue + + #Import sim sand temp data + sndTemps,sndTemps_time, start_date, end_date, run_number = readOutputSimuDatafile(df.loc[run,'fileID']) + + + #Process Simulated Sand Temp data into Numpy arrays + avg_sim_temps_snd_layer = np.mean(np.mean(sndTemps,axis=1),axis=1) + + #take average of all sand temps per timestep + avg_sim_temps_snd = np.mean(np.mean(np.mean(sndTemps,axis=1),axis=1),axis=1) + + + # Make list of dates for x-axis of simulation plots + sndTemps_datetimes = DateArrayFromSimDayHoursArray(sndTemps_time[:,0], start_date) + sndTemps_datetimes[1] += timedelta(hours = sndTemps_time[1,1]) #the offset of 0.001h needs to be put back into the array -> it was lost with thw function DateArrayFromSimDayHoursArray + + #were exp temps already imported? + if SimExpPlot == 1: + pass + else: + #it's the first run or the experimental data needed isn't the same as the last run + if (run == 0) or (diff_runs.Period[int(runs[run-1])] != diff_runs.Period[int(runs[run])]): + ##IMPORT MULTIPLE DATES OF SERRANO DATA + # import data + Exp_sand_temps, exp_start_date, exp_end_date = readSerranoSandMultipleDates(filename1) + + # Process Experimental Sand Temp data into Numpy arrays + Exp_sand_temps_4d,exp_time_array = processExpTCDataToArrays(Exp_sand_temps,exp_start_date) + + #take layer averages of sand temps (lowest index is bottom-most layer of sand store) + avg_expmtl_temps = np.mean(np.mean(Exp_sand_temps_4d,axis=1),axis=1) + + + #%%% Round time arrays to nearest hour + #initialize arrays + sndTemps_datetimes_rnd = [0]*len(sndTemps_datetimes) + exp_time_array_rnd= [0]*(len(sndTemps_datetimes) - 1) + avg_expmtl_SStemp4resid = np.full((len(exp_time_array_rnd)), np.NaN) + avg_expmtl_temps4resid = np.full((len(exp_time_array_rnd),3), np.NaN) + Exp_sand_temps4resid108 = np.full((len(exp_time_array_rnd),6,6,3), np.NaN) + + #round time to nearest hour for simulated data array + for time in range(2,len(sndTemps_datetimes)): # -> it starts at 2 because the value at index 1 is the preconditioning period time value + #sndTemps_datetimes is hourly data, make sure it rounds to the hour + sndTemps_datetimes_rnd[time] = sim_round2hr(sndTemps_datetimes[time], 60) + + #round time to nearest hour for experimental data array + exp_round2hr() + + #get average temps at each of the timesteps for calculating residuals + avg_expmtl_SStemp4resid = np.mean(avg_expmtl_temps4resid, axis=1) + + + #%%% Residuals of 3 avg layer temps and entire SS + + #initialize array + residSimVsExpLayers = pd.DataFrame(exp_time_array_rnd,columns=['datetime']) + col_list = ['LayerC','LayerB','LayerA', 'entireSS'] + residSimVsExpLayers[col_list] = np.nan + #remove first row, bc the 1st row of temperatures in both exp. and sim. arrays should be the same...that's how I initialized the simulation) + residSimVsExpLayers.drop(index=0, axis=0, inplace=True) + + #compare temps residuals and alert if the times aren't exactly the same + for time in range(0,len(sndTemps_datetimes_rnd)-2): + #layers: + residSimVsExpLayers.iloc[time,1] = avg_sim_temps_snd_layer[time+2,0] - avg_expmtl_temps4resid[(time+1),0] + residSimVsExpLayers.iloc[time,2] = avg_sim_temps_snd_layer[time+2,1] - avg_expmtl_temps4resid[(time+1),1] + residSimVsExpLayers.iloc[time,3] = avg_sim_temps_snd_layer[time+2,2] - avg_expmtl_temps4resid[(time+1),2] + #entire sand store: + residSimVsExpLayers.iloc[time,4] = avg_sim_temps_snd[time+2] - avg_expmtl_SStemp4resid[time+1] + + if (sndTemps_datetimes_rnd[time+2] != exp_time_array_rnd[time+1]): + #unless it's the very last timestep + if exp_time_array_rnd[time+1] != exp_time_array_rnd[-1]: + print("ERROR! There is an issue lining up the simulation and experimental temp data when calculating the RMSD in run " +str(runs[run]) + ","+ str(time)) + + #count the number of rows that are NaN (bc of experimental data not existing at the time stamp) + rows_nan = residSimVsExpLayers['LayerC'].isna().sum() + + #%%% RMSD of avg of 3 layers; total RMSD for run + + #using RSMD (forecasted - observed)^2 + residSimVsExpLayers = pd.concat([residSimVsExpLayers, + residSimVsExpLayers[col_list].pow(2, axis=0).add_prefix('sq_')], + axis=1) + + #initialize empty list + sq_col_list =['']*4 + #Make new column names and put into a list + for i in range(0,4): sq_col_list[i] = ('sq_'+col_list[i]) + + #sum the squares and divide by the number of entries, then take the sqrt + RMSD = residSimVsExpLayers[sq_col_list].sum().div(residSimVsExpLayers.shape[0] - rows_nan,axis=0).pow(0.5, axis=0) + + #make a row of the RMSD results -> RMSD of each of Layer C, B, A, RMSD of SSavg temp, and also the average RMSD of all the layers' RMSDs + temp_row = np.array([[int(runs[run]), + RMSD[0], #Layer C + RMSD[1], #Layer B + RMSD[2], #Layer A + RMSD[3], #entire SS RMSD + np.mean(RMSD[0:3])]]) #avg of layer RMSDs + + + + #%%% RMSD for 108 TCs + + # #initial temperature conditions are not the same for the simulated temps as for the exp arrays... + # #so calculate the difference in initial conditions between sim and exp for each TC + # delta_init = np.zeros(shape=(6,6,3)) + # res108 = np.zeros(shape=(len(exp_time_array_rnd)-1,6,6,3)) #make the length one row shorter bc we aren't going to take the resid at t=0s + # RMSD108 = np.zeros(shape=(6,6,3)) + + # for c in range(0,3): + # for b in range (0,6): + # for a in range (0,6): + # delta_init[a,b,c] = sndTemps[1,a,b,c] - Exp_sand_temps4resid108[0,a,b,c] + + # # calculate the residuals for each TC at each time step (except t=0), with the "corrected" sim temp + # # also start calculating the RMSD for each TC by summing the squares of all the residuals as the outer "for loop" runs + # nan_count = 0 + # for time in range(0,len(exp_time_array_rnd)-1): + # for c in range(0,3): + # for b in range (0,6): + # for a in range (0,6): + # res108[time,a,b,c] = (sndTemps[time+2,a,b,c] - delta_init[a,b,c]) - Exp_sand_temps4resid108[time+1,a,b,c] + # #check if residuals for that time step exist or if the values are NaN. Do not add the NaN values to the RMSD108 calcs + # if np.isnan(res108[time,a,b,c]) == True: + # nan_count += 1 + # else: + # RMSD108[a,b,c] = RMSD108[a,b,c] + res108[time,a,b,c]**2 + + # # calculate RMSD for each of the 108 TCs + # RMSD108 = RMSD108/(res108.shape[0] - nan_count/108) + # RMSD108 = np.sqrt(RMSD108) + + # #Calculate average RMSD for 108 TCs + # temp_row = np.append(temp_row ,[[np.mean(RMSD108)]], axis=1) + + + #%%% Plots + + if plot_individ_SS_graphs == 1: #if user wants each graph inidividually + + #%%%% Plots of exp and sim averages of entire SS individually + fig1, ax1 = plt.subplots() + colors='black' + labels = 'Exp data' + + # Plot experimental data + ax1.plot(exp_time_array_rnd,avg_expmtl_SStemp4resid,color=colors, label=labels) + + #Plot simulated data on same figure + colors='crimson' + labels = 'Sim data' + linestyle='--' + #less than r115 has rewinds, so we don't want to graph them...start at line 964 in TC sand temp data + if int(run_number) < 115: + ax1.plot(sndTemps_datetimes_rnd[964:],avg_sim_temps_snd[964:], color=colors, label=labels, linestyle=linestyle) + else: + ax1.plot(sndTemps_datetimes_rnd[2:],avg_sim_temps_snd[2:], color=colors, label=labels, linestyle=linestyle) + + + #plt.plot(sndTemps_datetimes,avg_sim_temps_snd_layer[:,2],marker='*') + #ax1.set_xlabel("Days since start of simulation") + ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + + #rename figure and change legend to include experimental data + # ax1.set_title("Mean temperature of sand store, Run " + run_number) + ax1.legend() + + + #Documentation on date locators + #https://matplotlib.org/stable/api/dates_api.html + locator = mdates.AutoDateLocator() + ax1.xaxis.set_major_locator(locator) + ax1.xaxis.set_minor_locator(mdates.DayLocator()) + + loc = ticker.MultipleLocator(base=2.0) + ax1.yaxis.set_major_locator(loc) + ax1.yaxis.set_minor_locator(ticker.MultipleLocator(base=1.0)) + + # #set the major locator on x-axis to every 6 days + # locator.intervald[mdates.DAILY]=[6] #https://matplotlib.org/stable/api/dates_api.html#matplotlib.dates.AutoDateLocator + + #Formatting - make the x-axis labels more consices + #https://matplotlib.org/stable/gallery/ticks_and_spines/date_concise_formatter.html + formatter = mdates.ConciseDateFormatter(locator) + ax1.xaxis.set_major_formatter(formatter) + fig1.autofmt_xdate() #automatically makes the x-labels rotate + + #gridlines on + ax1.grid(True) + + #Set max and min values for axes + ax1.set_xlim(exp_time_array[0], exp_time_array[-1]) + if mode == 6: + ax1.set_ylim(20,58) + else: + ax1.set_ylim(30,58) + + if Save_graphs == 1: + #save graph to file with proper naming scheme + plt.savefig(graph_dir.joinpath("r" + df.loc[run,'runs'] + "_entireSS.png"), dpi=150, bbox_inches='tight') + + if plot_multiple_on_same == 1: #if user wants each data set on the same figure + #%%%% Plots of exp and sim averages of entire SS + + if run==0: #if it's the first loop through the runs, plot the experimental data + colors='black' + labels = 'Exp data' + # Plot experimental data + # ax5.plot(exp_time_array_rnd,avg_expmtl_SStemp4resid,color=colors, label=labels) + + #Plot simulated data on same figure + # colors='crimson' + labels = 'Sim data' + linestyle=['dashed','dashed','dashed'] + clrs = ['blue', + 'g', + 'r', + 'dodgerblue', + 'darkorchid', + 'blue',] + #less than r115 has rewinds, so we don't want to graph them...start at line 964 in TC sand temp data + if int(run_number) < 115: + ax5.plot(sndTemps_datetimes_rnd[964:],avg_sim_temps_snd[964:], label=labels5[run], linestyle=linestyle[run],color=clrs[run]) + else: + ax5.plot(sndTemps_datetimes_rnd[2:],avg_sim_temps_snd[2:], label=labels5[run], linestyle=linestyle[run],color=clrs[run]) + + + #plt.plot(sndTemps_datetimes,avg_sim_temps_snd_layer[:,2],marker='*') + #ax5.set_xlabel("Days since start of simulation") + ax5.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + + #rename figure and change legend to include experimental data + # ax5.set_title("Comparing mean temperature of SSTES for different runs") + ax5.legend() + + + #Documentation on date locators + #https://matplotlib.org/stable/api/dates_api.html + locator = mdates.AutoDateLocator() + ax5.xaxis.set_major_locator(locator) + ax5.xaxis.set_minor_locator(mdates.DayLocator()) + + loc = ticker.MultipleLocator(base=2.0) + ax5.yaxis.set_major_locator(loc) + ax5.yaxis.set_minor_locator(ticker.MultipleLocator(base=1.0)) + + # #set the major locator on x-axis to every 6 days + # locator.intervald[mdates.DAILY]=[6] #https://matplotlib.org/stable/api/dates_api.html#matplotlib.dates.AutoDateLocator + + #Formatting - make the x-axis labels more consices + #https://matplotlib.org/stable/gallery/ticks_and_spines/date_concise_formatter.html + formatter = mdates.ConciseDateFormatter(locator) + ax5.xaxis.set_major_formatter(formatter) + fig5.autofmt_xdate() #automatically makes the x-labels rotate + + #gridlines on + ax5.grid(True) + + #Set max and min values for axes + ax5.set_xlim(exp_time_array[0], exp_time_array[-1]) + if mode == 6: + ax5.set_ylim(20,58) + else: + ax5.set_ylim(30,58) + + if (runs[run] == runs[-1] and Save_graphs == 1): + #save graph to file with proper naming scheme + plt.savefig(graph_dir.joinpath("r" + df.loc[run,'runs'] + "_entireSS_multiple.png"), dpi=150, bbox_inches='tight') + #%%%% Plots of residuals for averages of 3 layers + + # # Plot of residuals + # sns.lineplot(x=residSimVsExpLayers['datetime'],y=residSimVsExpLayers['LayerC']) + # sns.lineplot(x=residSimVsExpLayers['datetime'],y=residSimVsExpLayers['LayerB']) + # sns.lineplot(x=residSimVsExpLayers['datetime'],y=residSimVsExpLayers['LayerA']) + + + + #%%%% plot 108 residuals over time + # for c in range(0,3): + # fig,ax = plt.subplots() + # for b in range (0,6): + # for a in range (0,6): + # ax.plot(sndTemps_datetimes_rnd[2:-1],res108[:,a,b,c],label = str([a,b])) + # ax.legend() + # ax.set_title(col_list[c] + " residuals") + # ax.set_ylim([-3,3]) + # ax.set_ylabel("Temp, C") + + #%%%%plot 108 RMSD on a heatmap + # fig,axs = plt.subplots(3,1, figsize = (10,10)) + # for c in range(0,3): + # sns.heatmap(RMSD108[:,:,c], ax=axs[c], annot=True, cmap="Blues").invert_yaxis() + # axs[c].set_title(col_list[c] + ' TC RMSDs') + + #%%% Ouput RMSD info + + #Create RMSD array for the runs, or append to it + if run==0: #if the run is the first run processed + RMSD_array = temp_row + else: + RMSD_array = np.vstack((RMSD_array,temp_row)) + + #if last run, convert to Pandas dataframe + if runs[run] == runs[-1]: + RMSD_array_df = pd.DataFrame(RMSD_array,columns=["Run", + "Layer_C", + "Layer_B", + "Layer_A", + "EntireSS_RMSD", + "AvgRMSDofLayers", + ]) + #print dataframe to screen + print(RMSD_array_df) + + + +#%% 11. KM3 calculation + + #does this section need to be done? + if dE_ExpVSsim == 1: + ExpTCimport = 0 + + #%%% Check if temperature data has been imported + #were sim TC temps already imported? + if (SimExpPlot == 1 or crossrowsNcols == 1 or TdiffExpVSsim == 1): + pass + else: + + #Check if sim file exists + try: + open(df.loc[run,'fileID'], 'r') + except FileNotFoundError: #if it doesn't exist, print error and move on to next run + print('**************************') + print("File " + df.loc[run,'fileID'] + " does not exist!!!") + print('**************************') + #skip the rest of that run + continue + + #Import sim sand temp data + sndTemps,sndTemps_time, start_date, end_date, run_number = readOutputSimuDatafile(df.loc[run,'fileID']) + + + #Process Simulated Sand Temp data into Numpy arrays + avg_sim_temps_snd_layer = np.mean(np.mean(sndTemps,axis=1),axis=1) + + #take average of all sand temps per timestep + avg_sim_temps_snd = np.mean(np.mean(np.mean(sndTemps,axis=1),axis=1),axis=1) + + + # Make list of dates for x-axis of simulation plots + sndTemps_datetimes = DateArrayFromSimDayHoursArray(sndTemps_time[:,0], start_date) + sndTemps_datetimes[1] += timedelta(hours = sndTemps_time[1,1]) #the offset of 0.001h needs to be put back into the array -> it was lost with thw function DateArrayFromSimDayHoursArray + + #were exp temps already imported? + if (SimExpPlot == 1 or TdiffExpVSsim == 1): + pass + else: + #it's the first run or the experimental data needed isn't the same as the last run + if (run == 0) or (diff_runs.Period[int(runs[run-1])] != diff_runs.Period[int(runs[run])]): + ##IMPORT MULTIPLE DATES OF SERRANO DATA + # import data + Exp_sand_temps, exp_start_date, exp_end_date = readSerranoSandMultipleDates(filename1) + + # Process Experimental Sand Temp data into Numpy arrays + Exp_sand_temps_4d,exp_time_array = processExpTCDataToArrays(Exp_sand_temps,exp_start_date) + + #take layer averages of sand temps (lowest index is bottom-most layer of sand store) + avg_expmtl_temps = np.mean(np.mean(Exp_sand_temps_4d,axis=1),axis=1) + + #%%% Calc thermal properties of simulated SS + #(code copied from Master_v5_4.py on Jun 15, 2022, svn rev 1133) + + # %%%% 2. Geometry of Domain + # %%%% High-level geometry + + #x-y-z origin lies at north-west corner, bottom of sand store + sand_L = 6 #(m) length of sand store x = 0 m to x = 6 m + sand_W = 6 #(m) width of sand store y = 0 m to y = 6 m + sand_H = 3 #(m) height of sand store, from bottom z = 0 m to top z = 3 m + + XPS_th = 0.35 #(m) spec'd thickness of insulation around sand store (14") + conc_th = 0.10 #(m) spec'd thickness of concrete + + soil_top_th = 1.5 #(m) thickness of soil layer on top of sand store to surface + soil_nrthside_th = 3 #(m) thickness of soil layer on north side of sand store to soil near house + soil_th = 3 #(m) thickness of soil layer on S, E, W, and bottom sides of sand store to far-field temp + + no_TC = 0.5 #(m) buffer of 0.5m between sand store wall and 1st TCs (according to Briana Kemery's Visio document:https://chernode.mae.carleton.ca/sbes_svn/CHEeR_instrumentation/Sensor_connections_to_DAQs/Sand_store_TC_and_moisture_meter_placements.vsdx + node_space = (sand_L - 2*no_TC)/5 #Calculating the distance (m) between TCs in the sand store, assuming that there are 6 x 6 TCs equally spaced on each layer, + + no_pwr = 0.3 #(m) length around outside of sand store where the PEX tubes are not layed down (estimated from photos) + pwr_bott = 0.4 #(m) bottom z-value of x-z plane where PEX tubing is laid down + pwr_mid = 1.3 #(m) middle z-value of x-z plane where PEX tubing is laid down + pwr_upr = 2.2 #(m) top z-value of x-z plane where PEX tubing is laid down + + domain_L = soil_th + soil_nrthside_th + conc_th*2 + XPS_th*2 + sand_L #(m) length of numerical domain, x-axis + domain_W = soil_th*2 + conc_th*2 + XPS_th*2 + sand_W#(m) width of numerical domain, y-axis + domain_H = soil_th + soil_top_th + conc_th + XPS_th*2 + sand_H#(m) height of numerical domain, z-axis + + + #%%%% Node Placement (varying mesh size) + + #In x (or y-direction), from x=0 to the south. + #0m = north far-field temperature + + """ + The following values can be changed to see the effect of number of nodes on the final solution + """ + dx_soil = 0.3 #(m) dx of nodes in coarse grid horiz soil layer: in between the soil far-field boundary and the "increased node spacing layer" + dx_soilconc_boundary = 0.1 #(m) dx of nodes in fine grid horiz soil layer: between the coarse grid and the concrete layer + soilconc_boundary_th = 0.2 #(m) thickness of fine grid horiz soil layer + nodes_conc = 2 #number of nodes to make up concrete layer (horiz and vertical) + nodes_XPS = 2 #number of nodes to make up XPS layer (horiz and vertical) + nodes_no_pwr = 2 # number of horiz. nodes in sand layer, in the "no_pwr" area + dx_sand = 0.2 #(m) dx of horiz. nodes in sand layer, inside the "pwr" area + + dz_soil_top = 0.3 #(m)dz of nodes in top vert. soil layer + dz_pwr_th = 0.03 #(m) the vertical node spacing in the three power slab layers. Each power slab layer has three equi-sized dz layers + #around and on it, with the node of the middle layer directly on the power slab line (i.e. 0.4m, 1.3m, and 2.2m + #above the sand store bottom) + nodes_dz1_sand = 2 #number of nodes in between pwr_bott layer and bottom of sand store (not including the three nodes around the pwr_bott layer) + nodes_dz2_sand = 4 #number of nodes in between pwr_mid layer and pwr_bott layer (not including the three nodes around each pwr layer) + nodes_dz3_sand = 4 #number of nodes in between pwr_upr layer and pwr_mid layer (not including the three nodes around each pwr layer) + nodes_dz4_sand = 4 #number of nodes in between pwr_upr layer and top of sand store(not including the three nodes around pwr_upr layer) + + #%%%% Horizontal Coordinate Vectors + #%%%%% SOIL COORDS + + #(Note:the vector of dx's is the same as the vector for dy's, except for + #the length of the northside soil thickness) + #============================================================ + #============================================================ + #number of nodes in coarser soil grid + QN1 = math.floor(round((soil_nrthside_th - soilconc_boundary_th)/dx_soil,3)) + #remainder of coarse soil length needing to be accounted for + RN1 = (Decimal(soil_nrthside_th) - Decimal(soilconc_boundary_th)) % Decimal(dx_soil) + RN1 = check_mod(RN1,dx_soil) + + + #number of nodes in fine soil grid + Q2 = math.floor(round(soilconc_boundary_th/dx_soilconc_boundary,3)) + #remainder of fine soil length needing to be accounted + R2 = Decimal(soilconc_boundary_th) % Decimal(dx_soilconc_boundary) + R2 = check_mod(R2,dx_soilconc_boundary) + + #SUBVECTOR in which to put X-COORDINATES OF NORTH SOIL NODES + #-------------------------------------------------- + Hsoil_coords_n = np.zeros(shape=[1+QN1+Q2]) + + #far-field north boundary x=0 m + Hsoil_coords_n[0] = 0 + #1st node's x-coordinate, in metres: + Hsoil_coords_n[1] = (dx_soil/2) + RN1 + R2 + #starting at index 3, coarse nodes' x-coordinates in metres + for i in range(2,(1+QN1)): + Hsoil_coords_n[i] = Hsoil_coords_n[i-1] + dx_soil + + #1st node in fine grid area: + i = i+1 + Hsoil_coords_n[i] = Hsoil_coords_n[i-1] + dx_soil/2 + dx_soilconc_boundary/2 + i = i+1 + #fine nodes' x-coordinates in metres + for i in range(i,len(Hsoil_coords_n)): + Hsoil_coords_n[i] = Hsoil_coords_n[i-1] + dx_soilconc_boundary + + + #SUBVECTOR in which to put Y-COORDINATES OF WEST SOIL NODES + #-------------------------------------------------- + #number of nodes in coarser soil grid + Q1 = math.floor((soil_th - soilconc_boundary_th)/dx_soil) + #remainder of coarse soil length needing to be accounted for + R1 = (soil_th - soilconc_boundary_th) % dx_soil + + Hsoil_coords_w = np.zeros(shape=[1+Q1+Q2]) + + #far-field west boundary x=0 m + Hsoil_coords_w[0] = 0 + #1st node's y-coordinate, in metres: + Hsoil_coords_w[1] = (dx_soil/2) + R1 + R2 + #starting at index 3, coarse nodes' x-coordinates in metres + for i in range(2,(1+Q1)): + Hsoil_coords_w[i] = Hsoil_coords_w[i-1] + dx_soil + + #1st node in fine grid area: + i = i+1 + Hsoil_coords_w[i] = Hsoil_coords_w[i-1] + dx_soil/2 + dx_soilconc_boundary/2 + i = i+1 + #fine nodes' x-coordinates in metres + for i in range(i,len(Hsoil_coords_w)): + Hsoil_coords_w[i] = Hsoil_coords_w[i-1] + dx_soilconc_boundary + + + #SUBVECTOR in which to put X-COORDINATES OF SOUTH SOIL NODES + #-------------------------------------------------- + Hsoil_coords_s = np.zeros(shape=[1+Q1+Q2]) + + Hsoil_coords_s[(0)] = soil_nrthside_th + conc_th*2 + XPS_th*2 + sand_L + dx_soilconc_boundary/2 + #fine nodes' x-coordinates in metres + for i in range(1,Q2): + Hsoil_coords_s[i] = Hsoil_coords_s[i-1] + dx_soilconc_boundary + + #1st node in coarse grid area: + i = i+1 + Hsoil_coords_s[i] = Hsoil_coords_s[i-1] + dx_soil/2 + dx_soilconc_boundary/2 + i = i+1 + #coarse nodes' x-coordinates in metres + for i in range(i,(len(Hsoil_coords_s)-1)): + Hsoil_coords_s[i] = Hsoil_coords_s[i-1] + dx_soil + + #far-field south boundary x=domain_L + Hsoil_coords_s[(len(Hsoil_coords_s)-1)] = domain_L + + + #SUBVECTOR in which to put Y-COORDINATES OF EAST SOIL NODES + #-------------------------------------------------- + Hsoil_coords_e = np.zeros(shape=[1+Q1+Q2]) + + Hsoil_coords_e[(0)] = soil_th + conc_th*2 + XPS_th*2 + sand_W + dx_soilconc_boundary/2 + #fine nodes' x-coordinates in metres + for i in range(1,Q2): + Hsoil_coords_e[i] = Hsoil_coords_e[i-1] + dx_soilconc_boundary + + #1st node in coarse grid area: + i = i+1 + Hsoil_coords_e[i] = Hsoil_coords_e[i-1] + dx_soil/2 + dx_soilconc_boundary/2 + i = i+1 + #coarse nodes' x-coordinates in metres + for i in range(i,(len(Hsoil_coords_e)-1)): + Hsoil_coords_e[i] = Hsoil_coords_e[i-1] + dx_soil + + #far-field south boundary y=domain_W + Hsoil_coords_e[(len(Hsoil_coords_e)-1)] = domain_W + + + # %%%%% CONCRETE COORDS + #============================================================ + Hconc_coords_n = np.zeros(shape=[nodes_conc]) + Hconc_coords_w = np.zeros(shape=[nodes_conc]) + Hconc_coords_s = np.zeros(shape=[nodes_conc]) + Hconc_coords_e = np.zeros(shape=[nodes_conc]) + + #dx for concrete: + Q3 = math.floor(conc_th*100/nodes_conc)/100 + #remainder of concrete length needing to be accounted for + R3 = (conc_th*100 % nodes_conc)/100 + + #1st conc node + Hconc_coords_n[0] = soil_nrthside_th + Q3/2 + R3 + Hconc_coords_w[0] = soil_th + Q3/2 + R3 + Hconc_coords_s[0] = soil_nrthside_th + conc_th + XPS_th*2 + sand_L + Q3/2 + Hconc_coords_e[0] = soil_th + conc_th + XPS_th*2 + sand_W + Q3/2 + + #remainder of nodes + for i in range(1,nodes_conc): + Hconc_coords_n[i] = Hconc_coords_n[i-1] + Q3 + Hconc_coords_w[i] = Hconc_coords_w[i-1] + Q3 + Hconc_coords_s[i] = Hconc_coords_s[i-1] + Q3 + Hconc_coords_e[i] = Hconc_coords_e[i-1] + Q3 + + + #%%%%% XPS COORDS + #============================================================ + HXPS_coords_n = np.zeros(shape=[nodes_XPS]) + HXPS_coords_w = np.zeros(shape=[nodes_XPS]) + HXPS_coords_s = np.zeros(shape=[nodes_XPS]) + HXPS_coords_e = np.zeros(shape=[nodes_XPS]) + + #dx for XPS: + Q4 = math.floor(XPS_th*100/nodes_XPS)/100 + #remainder of XPS length needing to be accounted for + R4 = (XPS_th*100 % nodes_XPS)/100 + + #1st XPS node + HXPS_coords_n[0] = soil_nrthside_th + conc_th + Q4/2 + R4 + HXPS_coords_w[0] = soil_th + conc_th + Q4/2 + R4 + HXPS_coords_s[0] = soil_nrthside_th + conc_th + XPS_th + sand_L + Q4/2 + HXPS_coords_e[0] = soil_th + conc_th + XPS_th + sand_W + Q4/2 + + #remainder of nodes + for i in range(1,nodes_XPS): + HXPS_coords_n[i] = HXPS_coords_n[i-1] + Q4 + HXPS_coords_w[i] = HXPS_coords_w[i-1] + Q4 + HXPS_coords_s[i] = HXPS_coords_s[i-1] + Q4 + HXPS_coords_e[i] = HXPS_coords_e[i-1] + Q4 + + + #%%%%% SAND, NO-POWER ZONE COORDS (north and west sides) + #============================================================ + Hnopwr_coords_n = np.zeros(shape=[nodes_no_pwr]) + Hnopwr_coords_w = np.zeros(shape=[nodes_no_pwr]) + Hnopwr_coords_s = np.zeros(shape=[nodes_no_pwr]) + Hnopwr_coords_e = np.zeros(shape=[nodes_no_pwr]) + + #dx for no_pwr: + Q5 = math.floor(no_pwr*100/nodes_no_pwr)/100 + #remainder of no_pwr length needing to be accounted for + R5 = (no_pwr*100 % nodes_no_pwr)/100 + + #1st no_pwr node on north and west sides + Hnopwr_coords_n[0] = soil_nrthside_th + conc_th + XPS_th + Q5/2 + R5 + Hnopwr_coords_w[0] = soil_th + conc_th + XPS_th + Q5/2 + R5 + + #remainder of nodes on north and west sides + for i in range(1,nodes_no_pwr): + Hnopwr_coords_n[i] = Hnopwr_coords_n[i-1] + Q5 + Hnopwr_coords_w[i] = Hnopwr_coords_w[i-1] + Q5 + + + #%%%%% SAND, POWER ZONE COORDS + #============================================================ + #number of nodes in sand power area, x-dir + Q6 = math.floor(round((sand_L - no_pwr*2)/dx_sand,3)) + #remainder of sand length needing to be accounted for + R6 = (Decimal(sand_L) - Decimal(no_pwr)*2) % Decimal(dx_sand) + R6 = check_mod(R6,dx_sand) + + #number of nodes in sand power area, y-dir + Q7 = math.floor(round((sand_W - no_pwr*2)/dx_sand,3)) + #remainder of sand width needing to be accounted for + R7 = (Decimal(sand_W) - Decimal(no_pwr)*2) % Decimal(dx_sand) + R7 = check_mod(R7,dx_sand) + + + Hsand_coords_x = np.zeros(shape=[Q6]) + Hsand_coords_y = np.zeros(shape=[Q7]) + + Hsand_coords_x[0] = soil_nrthside_th + conc_th + XPS_th + no_pwr + dx_sand/2 + R6 + Hsand_coords_y[0] = soil_th + conc_th + XPS_th + no_pwr + dx_sand/2 + R7 + + #remainder of nodes + for i in range(1,Q6): + Hsand_coords_x[i] = Hsand_coords_x[i-1] + dx_sand + + for j in range(1,Q7): + Hsand_coords_y[(j)] = Hsand_coords_y[(j-1)] + dx_sand + + + #%%%%% SAND NO-POWER ZONE COORDS (south and east sides) + #============================================================ + #1st no_pwr node on south and east sides + Hnopwr_coords_s[0] = soil_nrthside_th + conc_th + XPS_th + sand_L - no_pwr + Q5/2 + Hnopwr_coords_e[0] = soil_th + conc_th + XPS_th + sand_W - no_pwr + Q5/2 + + #remainder of nodes on south and east sides + for i in range(1,nodes_no_pwr): + Hnopwr_coords_s[i] = Hnopwr_coords_s[i-1] + Q5 + Hnopwr_coords_e[i] = Hnopwr_coords_e[i-1] + Q5 + + + #%%%%% PUTTING X- AND Y- MESH COORDINATES VECTORS TOGETHER + #============================================================ + node_xcoords = np.concatenate((Hsoil_coords_n,Hconc_coords_n,HXPS_coords_n,Hnopwr_coords_n,Hsand_coords_x,Hnopwr_coords_s,HXPS_coords_s,Hconc_coords_s,Hsoil_coords_s)) + + node_ycoords = np.concatenate((Hsoil_coords_w,Hconc_coords_w,HXPS_coords_w,Hnopwr_coords_w,Hsand_coords_y,Hnopwr_coords_e,HXPS_coords_e,Hconc_coords_e,Hsoil_coords_e)) + + + #%%%% Vertical Coordinate Vectors + #%%%%% SOIL COORDS + + ''' + #Reminder: + #number of nodes in coarser soil grid = Q1 + #remainder of coarse soil length needing to be accounted for = R1 + + #number of nodes in fine soil grid = Q2 + #remainder of fine soil length needing to be accounted for = R2 + ''' + + #SUBVECTOR in which to put Z-COORDINATES OF BOTTOM SOIL NODES + #-------------------------------------------------- + Vsoil_coords_b = np.zeros(shape=[1+Q1+Q2]) + + #far-field bottom boundary z=0 m + Vsoil_coords_b[0] = 0 + #1st node's z-coordinate, in metres: + Vsoil_coords_b[1] = (dx_soil/2) + R1 + R2 + #starting at index 3, coarse nodes' z-coordinates in metres + for i in range(2,(1+Q1)): + Vsoil_coords_b[i] = Vsoil_coords_b[i-1] + dx_soil + + #1st node in fine grid area: + i = i+1 + Vsoil_coords_b[i] = Vsoil_coords_b[i-1] + dx_soil/2 + dx_soilconc_boundary/2 + i = i+1 + #fine nodes' z-coordinates in metres + for i in range(i,len(Vsoil_coords_b)): + Vsoil_coords_b[i] = Vsoil_coords_b[i-1] + dx_soilconc_boundary + + + + #SUBVECTOR in which to put Z-COORDINATES OF TOP SOIL NODES + #-------------------------------------------------- + #number of nodes in coarser soil grid + Q8 = math.floor((soil_top_th - soilconc_boundary_th)/dz_soil_top) + #remainder of coarse soil length needing to be accounted for + R8 = (soil_top_th - soilconc_boundary_th) % dz_soil_top + + Vsoil_coords_t = np.zeros(shape=[1+Q8+Q2]) + + Vsoil_coords_t[0] = soil_th + conc_th + XPS_th*2 + sand_H + dx_soilconc_boundary/2 + #fine nodes' z-coordinates in metres + for i in range(1,Q2): + Vsoil_coords_t[i] = Vsoil_coords_t[i-1] + dx_soilconc_boundary + + #1st node in coarse grid area: + i = i+1 + Vsoil_coords_t[i] = Vsoil_coords_t[i-1] + dz_soil_top/2 + dx_soilconc_boundary/2 + i = i+1 + #coarse nodes' z-coordinates in metres + for i in range(i,len(Vsoil_coords_t)-1): + Vsoil_coords_t[i] = Vsoil_coords_t[i-1] + dz_soil_top + + #far-field south boundary x=domain_L + Vsoil_coords_t[len(Vsoil_coords_t)-1] = domain_H + + + #%%%%% CONCRETE COORDS + #============================================================ + Vconc_coords_b = np.zeros(shape=[nodes_conc]) + + #dx for concrete = Q3: + #remainder of concrete length needing to be accounted for = R3 + + #1st conc node + Vconc_coords_b[0] = soil_th + Q3/2 + R3 + + #remainder of nodes + for i in range(1,nodes_conc): + Vconc_coords_b[i] = Vconc_coords_b[i-1] + Q3 + + + #%%%%% XPS COORDS + #============================================================ + VXPS_coords_b = np.zeros(shape=[nodes_XPS]) + VXPS_coords_t = np.zeros(shape=[nodes_XPS]) + + #dx for XPS = Q4 + #remainder of XPS length needing to be accounted for = R4 + + #1st XPS node + VXPS_coords_b[0] = soil_th + conc_th + Q4/2 + R4 + VXPS_coords_t[0] = soil_th + conc_th + XPS_th + sand_H + Q4/2 + + #remainder of nodes + for i in range(1,nodes_XPS): + VXPS_coords_b[i] = VXPS_coords_b[i-1] + Q4 + VXPS_coords_t[i] = VXPS_coords_t[i-1] + Q4 + + + #%%%%% SAND COORDS + #============================================================ + #dz in space in between pwr_bott layer and bottom of sand store (not including the three nodes around the pwr_bott layer) + dz1 = (pwr_bott - (dz_pwr_th + dz_pwr_th/2))/nodes_dz1_sand + #initialize vector + Vsand1 = np.zeros(shape=[nodes_dz1_sand+3]) + #1st node coordinate + Vsand1[0] = soil_th + conc_th + XPS_th + dz1/2 + #remainder of nodes + for i in range(1,nodes_dz1_sand): + Vsand1[i] = Vsand1[i-1] + dz1 + + #1st node in "3-node pwr layer" + i = i+1 + Vsand1[i] = Vsand1[i-1] + dz1/2 + dz_pwr_th/2 + i = i+1 + #remainder of nodes in "3-node pwr layer" + for i in range(i,len(Vsand1)): + Vsand1[i] = Vsand1[i-1] + dz_pwr_th + + + #********* + + #dz in space #number of nodes in between pwr_mid layer and pwr_bott layer (not including the three nodes around each pwr layer) + dz2 = ((pwr_mid - (dz_pwr_th + dz_pwr_th/2)) - (pwr_bott + (dz_pwr_th + dz_pwr_th/2)))/nodes_dz2_sand + #initialize vector + Vsand2 = np.zeros(shape=[nodes_dz2_sand+3]) + #1st node coordinate + Vsand2[0] = soil_th + conc_th + XPS_th + pwr_bott + (dz_pwr_th + dz_pwr_th/2) + dz2/2 + #remainder of nodes + for i in range(1,nodes_dz2_sand): + Vsand2[i] = Vsand2[i-1] + dz2 + + #1st node in "3-node pwr layer" + i = i+1 + Vsand2[i] = Vsand2[i-1] + dz2/2 + dz_pwr_th/2 + i = i+1 + #remainder of nodes in "3-node pwr layer" + for i in range(i,len(Vsand2)): + Vsand2[i] = Vsand2[i-1] + dz_pwr_th + + + #********* + + #dz in space #number of nodes in between pwr_upr layer and pwr_mid layer (not including the three nodes around each pwr layer) + dz3 = ((pwr_upr - (dz_pwr_th + dz_pwr_th/2)) - (pwr_mid + (dz_pwr_th + dz_pwr_th/2)))/nodes_dz3_sand + #initialize vector + Vsand3 = np.zeros(shape=[nodes_dz3_sand+3]) + #1st node coordinate + Vsand3[0] = soil_th + conc_th + XPS_th + pwr_mid + (dz_pwr_th + dz_pwr_th/2) + dz3/2 + #remainder of nodes + for i in range(1,nodes_dz3_sand): + Vsand3[i] = Vsand3[i-1] + dz3 + + + #1st node in "3-node pwr layer" + i = i+1 + Vsand3[i] = Vsand3[i-1] + dz3/2 + dz_pwr_th/2 + i = i+1 + #remainder of nodes in "3-node pwr layer" + for i in range(i,len(Vsand3)): + Vsand3[i] = Vsand3[i-1] + dz_pwr_th + + + #********* + + #dz in space #number of nodes in between pwr_upr layer and pwr_mid layer (not including the three nodes around each pwr layer) + dz4 = (sand_H - (pwr_upr + (dz_pwr_th + dz_pwr_th/2)))/nodes_dz4_sand + #initialize vector + Vsand4 = np.zeros(shape=[nodes_dz4_sand]) + #1st node coordinate + Vsand4[0] = soil_th + conc_th + XPS_th + pwr_upr + (dz_pwr_th + dz_pwr_th/2) + dz4/2 + #remainder of nodes + for i in range(1,nodes_dz4_sand): + Vsand4[i] = Vsand4[i-1] + dz4 + + + + #%%%%% PUTTING Z- COORDINATE VECTORS TOGETHER + #============================================================ + node_zcoords = np.concatenate((Vsoil_coords_b,Vconc_coords_b,VXPS_coords_b,Vsand1,Vsand2,Vsand3,Vsand4,VXPS_coords_t,Vsoil_coords_t)) + # figure + # plot(ones(length(node_zcoords),1),node_zcoords,'o') + + + + #%%%%VERTICAL index boundaries (6 layers) + indx_bndrz_z = np.array([0, #bottom soil + len(Vsoil_coords_b)-1, #bottom soil + len(Vsoil_coords_b), #bottom conc + len(Vsoil_coords_b)+len(Vconc_coords_b)-1, #bottom conc + len(Vsoil_coords_b)+len(Vconc_coords_b), #bottom XPS + len(Vsoil_coords_b)+len(Vconc_coords_b)+len(VXPS_coords_b)-1, #bottom XPS + len(Vsoil_coords_b)+len(Vconc_coords_b)+len(VXPS_coords_b), #sand + len(Vsoil_coords_b)+len(Vconc_coords_b)+len(VXPS_coords_b)+len(Vsand1)+len(Vsand2)+len(Vsand3)+len(Vsand4)-1, #sand + len(Vsoil_coords_b)+len(Vconc_coords_b)+len(VXPS_coords_b)+len(Vsand1)+len(Vsand2)+len(Vsand3)+len(Vsand4), #top XPS + len(Vsoil_coords_b)+len(Vconc_coords_b)+len(VXPS_coords_b)+len(Vsand1)+len(Vsand2)+len(Vsand3)+len(Vsand4)+len(VXPS_coords_t)-1,#top XPS + len(Vsoil_coords_b)+len(Vconc_coords_b)+len(VXPS_coords_b)+len(Vsand1)+len(Vsand2)+len(Vsand3)+len(Vsand4)+len(VXPS_coords_t),#top soil + len(node_zcoords)-1])#top soil + + #Indices for the top node (inclusive) of each "moisture layer boundary"...appx middle points between moisture sensor z-coordinates + moisture_layer_height = np.zeros(shape=[3], dtype=int) + moisture_layer_height[0] = int(len(Vsoil_coords_b)+len(Vconc_coords_b)+len(VXPS_coords_b)+len(Vsand1)+ math.floor(nodes_dz2_sand/2) - 1) # -1 to acct for zero-based arrays + moisture_layer_height[1] = int(len(Vsoil_coords_b)+len(Vconc_coords_b)+len(VXPS_coords_b)+len(Vsand1)+len(Vsand2)+math.floor(nodes_dz3_sand/2) - 1) # -1 to acct for zero-based arrays + moisture_layer_height[2] = int(indx_bndrz_z[7]) + + vol_slab = (sand_L - 2*no_pwr)*(sand_W - 2*no_pwr)*dz_pwr_th #(m3) volume of each power source slab + + #Define number of nodes in each direction in numerical model domain + + U = len(node_xcoords) #number of nodes in x-direction + V = len(node_ycoords) #number of nodes in y-direction + W = len(node_zcoords) #number of nodes in z-direction + + + pwr_bott_indx = indx_bndrz_z[5] + nodes_dz1_sand + 2 #[m] k-index for bottom layer of PEX tubing + pwr_mid_indx = indx_bndrz_z[5] + len(Vsand1) + nodes_dz2_sand + 2 #[m] k-index for middle layer of PEX tubing + pwr_upr_indx = indx_bndrz_z[5] + len(Vsand1) + len(Vsand2) + nodes_dz3_sand + 2 #[m] k-index for top layer of PEX tubing + + #%%%%HORIZONTAL index boundaries [9 layers] + indx_bndrz_x = np.array([0, #north soil + len(Hsoil_coords_n)-1, #north soil + len(Hsoil_coords_n), #north conc + len(Hsoil_coords_n)+len(Hconc_coords_n)-1, #north conc + len(Hsoil_coords_n)+len(Hconc_coords_n), #north XPS + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)-1, #north XPS + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n), #north sand no power + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+len(Hnopwr_coords_n)-1, #north sand no power + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+len(Hnopwr_coords_n), #sand + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+ + len(Hnopwr_coords_n)+len(Hsand_coords_x)-1,#sand + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+ + len(Hnopwr_coords_n)+len(Hsand_coords_x), #south sand no pwr + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+ + len(Hnopwr_coords_n)+len(Hsand_coords_x)+len(Hnopwr_coords_s)-1, #south sand no pwr + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+ + len(Hnopwr_coords_n)+len(Hsand_coords_x)+len(Hnopwr_coords_s), #south XPS + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+ + len(Hnopwr_coords_n)+len(Hsand_coords_x)+len(Hnopwr_coords_s)+ + len(HXPS_coords_s)-1, #south XPS + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+ + len(Hnopwr_coords_n)+len(Hsand_coords_x)+len(Hnopwr_coords_s)+ + len(HXPS_coords_s), #south conc + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+ + len(Hnopwr_coords_n)+len(Hsand_coords_x)+len(Hnopwr_coords_s)+ + len(HXPS_coords_s)+len(Hconc_coords_s)-1, #south conc + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+ + len(Hnopwr_coords_n)+len(Hsand_coords_x)+len(Hnopwr_coords_s)+ + len(HXPS_coords_s)+len(Hconc_coords_s), #south soil + len(node_xcoords)-1]) #south soil + + + indx_bndrz_y = np.array([0, #west soil + len(Hsoil_coords_w)-1, #west soil + len(Hsoil_coords_w), #west conc + len(Hsoil_coords_w)+len(Hconc_coords_w)-1, #west conc + len(Hsoil_coords_w)+len(Hconc_coords_w), #west XPS + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)-1, #west XPS + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w), #west sand no power + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+len(Hnopwr_coords_w)-1, #west sand no power + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w), #sand + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w)+len(Hsand_coords_y)-1,#sand + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w)+len(Hsand_coords_y), #east sand no pwr + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w)+len(Hsand_coords_y)+len(Hnopwr_coords_e)-1, #east sand no pwr + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w)+len(Hsand_coords_y)+len(Hnopwr_coords_e), #east XPS + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w)+len(Hsand_coords_y)+len(Hnopwr_coords_e)+ + len(HXPS_coords_e)-1, #east XPS + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w)+len(Hsand_coords_y)+len(Hnopwr_coords_e)+ + len(HXPS_coords_e), #east conc + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w)+len(Hsand_coords_y)+len(Hnopwr_coords_e)+ + len(HXPS_coords_e)+len(Hconc_coords_e)-1, #east conc + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w)+len(Hsand_coords_y)+len(Hnopwr_coords_e)+ + len(HXPS_coords_e)+len(Hconc_coords_e), #east soil + len(node_ycoords)-1]) #east soil + + + #%%%%Indices for points within the "power slab" in the x-y plane + x_pwr_low_indx = indx_bndrz_x[8] #inclusive of this index + x_pwr_high_indx = indx_bndrz_x[9]#inclusive of this index + y_pwr_low_indx = indx_bndrz_y[8] #inclusive of this index + y_pwr_high_indx = indx_bndrz_y[9]#inclusive of this index + + #============================================================ + #%%%%Vectors with CARTESIAN COORDS OF LAYER BOUNDARIES + #north-south + N_S_layr_bounds = np.array([0, #north boundary + soil_nrthside_th, #soil/conc + soil_nrthside_th + conc_th, #conc/XPS + soil_nrthside_th + conc_th + XPS_th , #XPS/sand + soil_nrthside_th + conc_th + XPS_th + sand_L, #sand/XPS + soil_nrthside_th + conc_th + XPS_th*2 + sand_L, #XPS/conc + soil_nrthside_th + conc_th*2 + XPS_th*2 + sand_L,#conc/soil + soil_nrthside_th + conc_th*2 + XPS_th*2 + sand_L + soil_th]) #south boundary + + + #west-east + W_E_layr_bounds =np.array([0, #west boundary + soil_th, #soil/conc + soil_th + conc_th, #conc/XPS + soil_th + conc_th + XPS_th , #XPS/sand + soil_th + conc_th + XPS_th + sand_W, #sand/XPS + soil_th + conc_th + XPS_th*2 + sand_W, #XPS/conc + soil_th + conc_th*2 + XPS_th*2 + sand_W,#conc/soil + soil_th*2 + conc_th*2 + XPS_th*2 + sand_W]) #east boundary + + #bottom-top + B_T_layr_bounds =np.array([0, #bottom boundary + soil_th, #soil/conc + soil_th + conc_th, #conc/XPS + soil_th + conc_th + XPS_th , #XPS/sand + soil_th + conc_th + XPS_th + sand_H, #sand/XPS + soil_th + conc_th + XPS_th*2 + sand_H, #XPS/soil + soil_th + conc_th + XPS_th*2 + sand_H + soil_top_th]) #top boundary + + #============================================================ + #%%%%Vectors holding the x and y coordinates of the TCs in the sand store + TC_xcoords = np.zeros(shape=[6]) + TC_ycoords = np.zeros(shape=[6]) + + TC_xcoords[0] = soil_nrthside_th + conc_th + XPS_th + no_TC + TC_ycoords[0] = soil_th + conc_th + XPS_th + no_TC + + for i in range(1,6): + TC_xcoords[i] = TC_xcoords[i-1] + node_space + TC_ycoords[i] = TC_ycoords[i-1] + node_space + + + #Find the indices of the TC locations + sand_TC_node_xindices = np.zeros(shape=[6], dtype=int) + + for TC in range(0,6): + diff = 0 + prev = 100 + for i in range(0,len(node_xcoords)): + diff = abs(TC_xcoords[TC] - node_xcoords[i]) + #this finds the closest nodes to the real locations of the TCs + if round(diff,5) >= round(prev,5): + #then we want the index represented by [i-1] + sand_TC_node_xindices[TC] = i-1 + break + prev = diff + + + sand_TC_node_yindices = np.zeros(shape=[6], dtype=int) + + for TC in range(0,6): + diff = 0 + prev = 100 + for j in range(0,len(node_ycoords)): + diff = abs(TC_ycoords[TC] - node_ycoords[j]) + #this finds the closest nodes to the real locations of the TCs + if round(diff,5) >= round(prev,5): + #then we want the index represented by [j-1] + sand_TC_node_yindices[TC] = j-1 + break + prev = diff + + + # %%%% 2. Material Thermal Properties + ''' + This script defines the thermal properties of the materials used in the simulation and also builds the 3-D thermal property matrices for the simulation. + ''' + #Properties of water + water_Cp = 4.186*1e3 # (J/kgK) specific heat for water + water_rho = 998.2 # (kg/m3) density for water at 293 K https://hypertextbook.com/facts/2007/AllenMa.shtml + + #Properties of sand + drySand_Cp = 830 #(J/kgK) specific heat for dry sand https://www.engineeringtoolbox.com/specific-heat-capacity-d_391.html + drySand_rho = 1850 #(kg/m3) density for dry sand (personal correspondence with Abdulghader Abdulrahman, Research Manager, Dr.Paul Simms' lab) + + #volume percents of water layers in sand store + vol_percent = np.zeros(shape=[3]) + vol_percent[0] = diff_runs['vol% bott'][int(runs[run])] #bottom layer = 33 vol# + vol_percent[1] = diff_runs['vol% mid'][int(runs[run])] #middle layer = 11 vol# + vol_percent[2] = diff_runs['vol% top'][int(runs[run])] #top layer = 7 vol# + + sand_lambda = np.zeros(shape=[3]) + sand_lambda[0] = diff_runs['λ_bott'][int(runs[run])] #(W/mK) bottom layer - thermal conductivity for saturated sand F. Rad (2009) + sand_lambda[1] = diff_runs['λ_mid'][int(runs[run])] #(W/mK) middle layer - moist sand https://www.engineeringtoolbox.com/thermal-conductivity-d_429.html + sand_lambda[2] = diff_runs['λ_top'][int(runs[run])] #(W/mK) top layer - moist sand https://www.engineeringtoolbox.com/thermal-conductivity-d_429.html + + #Properties of concrete at 300 K (Source - Incropera, 6e, p.939) + conc_Cp = 880 #(J/kgK) specific heat + conc_rho = 2300 # (kg/m3) density + conc_lambda = 1.4 #(W/mK) thermal conductivity + + #Properties of XPS at 297 K(Source - Owens Corning data sheets and papers - see ppt presentation MCG5157 final project for sources) + XPS_Cp_orig = 1500 #(J/kgK) specific heat source - article by BASF staff + XPS_rho_orig = 30 # (kg/m3) density + + #Properties of XPS + vol_percent_H2O_SETface_XPS = diff_runs['SET(B) XPS Vol % H2O'][int(runs[run])] + XPS_lambda_SET = diff_runs['SET(B)'][int(runs[run])] #(W/mK) thermal conductivity + XPS_Cp_SET = ((water_Cp * (water_rho*vol_percent_H2O_SETface_XPS)) + (XPS_Cp_orig*XPS_rho_orig)) / (water_rho * vol_percent_H2O_SETface_XPS + XPS_rho_orig) + # (Number of Joules to raise sand in m3 by 1K + Joules to raise water in 1m3 by 1K) / (number of kg in 1m3) + XPS_rho_SET = water_rho * vol_percent_H2O_SETface_XPS + XPS_rho_orig + + + #Properties of soil (soil_rhoCp and soil_lambda) around CHEeR house sand + #store + soil_Cp = 750 #(J/kgK) (A. Wills MASc,Table 6.1) + soil_rho = 2800 #(kg/m3) (A. Wills MASc,Table 6.1) + soil_lambda = 2 #(W/mK) thermal conductivity (A. Wills MASc,Table 6.1) (TRY PARAMETRIZING FOR CLAY VALUES ) + + + #%%%% 2. Build matrices of thermal properties + # (from outside layers inwards) + + U_snd = 31 + V_snd = 31 + W_snd = 23 + z_shift = indx_bndrz_z[6] + x_shift = indx_bndrz_x[6] + y_shift = indx_bndrz_y[6] + #Initialize thermal property matrices + Lambda = np.zeros(shape=[U_snd,V_snd,W_snd]) + Cp = np.zeros(shape=[U_snd,V_snd,W_snd]) + rho = np.zeros(shape=[U_snd,V_snd,W_snd]) + + #----------------------- + + #1st Sand layer [0m to ~halfway between first two PEX layers] + for k in range(indx_bndrz_z[6] - z_shift,moisture_layer_height[0] - z_shift +1): + for i in range(indx_bndrz_x[6] - x_shift,indx_bndrz_x[11] - x_shift+1): + for j in range(indx_bndrz_y[6] - y_shift,indx_bndrz_y[11]+1 - y_shift): + Lambda[i,j,k] = sand_lambda[0] + Cp[i,j,k] = ((water_Cp * (water_rho*vol_percent[0])) + (drySand_Cp*drySand_rho)) / (water_rho * vol_percent[0] + drySand_rho) + # (Number of Joules to raise sand in m3 by 1K + Joules to raise water in 1m3 by 1K) / (number of kg in 1m3) + rho[i,j,k] = water_rho * vol_percent[0] + drySand_rho + + + + + #2nd Sand layer [ of 1st layer to ~halfway between next two PEX layers] + for k in range((moisture_layer_height[0]+1 - z_shift),moisture_layer_height[1]+1 - z_shift): + for i in range(indx_bndrz_x[6] - x_shift,indx_bndrz_x[11]+ 1 - x_shift): + for j in range(indx_bndrz_y[6] - y_shift,indx_bndrz_y[11]+1 - y_shift): + Lambda[i,j,k] = sand_lambda[1] + Cp[i,j,k] = ((water_Cp * (water_rho*vol_percent[1])) + (drySand_Cp*drySand_rho)) / (water_rho * vol_percent[1] + drySand_rho) + rho[i,j,k] = water_rho * vol_percent[1] + drySand_rho + + + + + #3rd Sand layer [ of 2nd layer to top of sand store] + for k in range((moisture_layer_height[1]+1 - z_shift),moisture_layer_height[2]+1 - z_shift): + for i in range(indx_bndrz_x[6] - x_shift,indx_bndrz_x[11]+1 - x_shift): + for j in range(indx_bndrz_y[6] - y_shift,indx_bndrz_y[11]+1 - y_shift): + Lambda[i,j,k] = sand_lambda[2] + Cp[i,j,k] = ((water_Cp * (water_rho*vol_percent[2])) + (drySand_Cp*drySand_rho)) / (water_rho * vol_percent[2] + drySand_rho) + rho[i,j,k] = water_rho * vol_percent[2] + drySand_rho + + + #%%% Calc avg Cp, rho of SSTES + + avg_SSCp = np.mean(Cp) + avg_SSrho = np.mean(rho) + + #%%%Calc KM3 metric + #%%%%Heat input to Sand Store + #Read in Headers file + dat_headers = pd.read_csv(data_dir.joinpath('Serrano_dat_headers.csv')) + + #Read in loops heat transfer data into DataFrame + Exp_heat_transfer_data = func_serr.read_loops_heat_transfer(filename2, dat_headers, main_dir, data_dir) + + #Process heat input data + #Put heat source data into a new numpy array + solar_to_Sand_W = Exp_heat_transfer_data[['time','sand_H2O']].copy().to_numpy() + + #Integrate over the rate of energy going into the sand store to find total energy added + int_energy_added = 0 + for row in range(0,solar_to_Sand_W.shape[0]): + int_energy_added += 5 * 60 * solar_to_Sand_W[row,1] # 5 minute time intervals multiplied by the Watts injected + + int_energy_added = int_energy_added/1e9 #(convert J to GJ) + + # %%%% Heat extracted from Sand Store + # Process discharging data + + #Put discharging data into a new numpy array + Sand_to_SH_W = Exp_heat_transfer_data[['time','sand_2_SH']].copy().to_numpy() + + #Integrate over the rate of energy going out of the sand store to find total energy extracted + int_energy_taken = 0 + for row in range(0,Sand_to_SH_W.shape[0]): + int_energy_taken += 5 * 60 * Sand_to_SH_W[row,1] # 5 minute time intervals multiplied by the Watts injected + + int_energy_taken = int_energy_taken/1e9 #(convert J to GJ) + + + + + #%%% Calc dE_tot + + #let the reference temp equal to 15C. Calculate all dE in reference to it (see p5.53 in notes) + ref_temp = 30 + #pre-calculate mCp, from dE = mCpdT + mCp = (avg_SSrho * 108) * avg_SSCp # mCp = rho*volume*Cp + # 108 m^3 is volume of SS + + #starting temp for both sim and exp temperatures + T1 = avg_sim_temps_snd[0] + + #last temperature in simulated temp array + T2_sim = avg_sim_temps_snd[-1] + + #last temperature in exp temp array + try: #check if Section 10 was run - this pared-down hourly temperature array will exist if it was + T2_exp = avg_expmtl_SStemp4resid[-1] + except NameError: #if not, take the avg of the last timestep in the experimental temp array, which wasn't pared down to hourly data + T2_exp = np.mean(avg_expmtl_temps[-1,:]) + + #Calculate delta energies + dE1 = mCp * (T1 - ref_temp) /1e9 #(convert J to GJ) + dE2_sim = mCp * (T2_sim - ref_temp) /1e9 #(convert J to GJ) + dE2_exp = mCp * (T2_exp - ref_temp) /1e9 #(convert J to GJ) + + dE_tot_sim = dE1 - dE2_sim + dE_tot_exp = dE1 - dE2_exp + + + #Calc dE_heatloss + dE_HL = dE1 + int_energy_added - dE2_exp - int_energy_taken + + #Calc normalized quantities + denom = abs(int_energy_added) + abs(int_energy_taken) + abs(dE_HL) + + dE_norm_exp = dE_tot_exp / denom + dE_norm_sim = dE_tot_sim / denom + + KM3 = (dE2_exp - dE2_sim)/denom *100 + + + #make a row of the results + temp_row1 = np.array([[int(runs[run]), + KM3, + dE1, + dE2_exp, + dE2_sim, + abs(int_energy_added), + abs(int_energy_taken), + abs(dE_HL)]]) + + + # #the simulated change in stored Energy as a percent of the experimental change in stored Energy + # pct_energy_metric = dE_tot_sim/dE_tot_exp * 100 + + # print("\nExperiment: %.2f GJ reduction in STES stored energy" % (dE_tot_exp/1e9)) + # print("Simulation: %.2f GJ reduction in STES stored energy \n" % (dE_tot_sim/1e9)) + + # print("The STES started with %.2f GJ of stored energy, above a zero reference energy at 15\N{degree sign}C." % (dE1/1e9)) + + #%%% Ouput KM3 metric + + #Create RMSD array for the runs, or append to it + if run==0: #if the run is the first run processed + KM3_array = temp_row1 + else: + KM3_array = np.vstack((KM3_array,temp_row1)) + + #if last run, convert to Pandas dataframe + if runs[run] == runs[-1]: + KM3_array_df = pd.DataFrame(KM3_array,columns=["Run", + "KM3", + "Stored_E_init(GJ)", + "Stored_E_fin_exp(GJ)", + "Stored_E_fin_sim(GJ)", + "E_ch (GJ)", + "E_disch (GJ)", + "E_hl (GJ)"]) + #print dataframe to screen + print(KM3_array_df) + + +#%% 12. PLOTS FOR PUBLISHING + +#%%% PCP/initialization + +# #Plot first two data points of viz data, to show SS state after initialization +# fig1,ax1 = plt.subplots() + +# #rewinds no longer being used at run 115 and above +# if int(run_number_ID2) >= 115: +# for time in range(0,1): +# print("Graphed point (day,hour): " + str(sndVizTemps_time[time])) +# ax1.plot(x_coords2plt[:,0],sndVizTemps[time,:,bb,cc]) +# #plot second temp distribution in array, steady state temp distribution after IC +# ax1.plot(x_coords2plt[:,0],sndVizTemps[1,:,bb,cc],linestyle='--') +# print("Graphed point (day,hour): " + str(sndVizTemps_time[1])) + +# ax1.legend(["IC", +# "After PCP"], loc='upper right') + +# #plot vertical shading for XPS and conc +# ax1.axvspan(x_conc_coord1,x_conc_coord2,color='grey', alpha=0.5, lw=0) +# ax1.axvspan(x_conc_coord3,x_conc_coord4,color='grey', alpha=0.5, lw=0) +# ax1.axvspan(x_XPS_coord1,x_XPS_coord2,color='pink', alpha=0.5, lw=0) +# ax1.axvspan(x_XPS_coord3,x_XPS_coord4,color='pink', alpha=0.5, lw=0) + +# # ax1.set_title("Run "+run_number_ID2+ ": Daily temp distribution for " + str(round(sndVizTemps_time[-1,0])) + "days across x-axis, at y=" + str(y_viz_indices[bb]) + " and z=" + str(z_viz_indices[cc])) +# ax1.set_ylim(0,60) +# ax1.set_xlabel('Length of num. domain (m)') +# ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + +#%% 13. PLOT PERFORMANCE METRICS DATA +if PerfMetr == 1: +#%%% Import performance metrics data + Perf_metrics, exp_start_date2, exp_end_date2 = readSerranoPerformanceMetricsData(filename2) + + +#%%% Plot measured ground temps near Sand + + #create empty list + exp_time_array_days2 = [0]*Perf_metrics.shape[0] + #put start date in top row of list + exp_time_array_days2[0] = exp_start_date2 + for d in range(1,Perf_metrics.shape[0]): + exp_time_array_days2[d] = exp_time_array_days2[d-1] + timedelta(days=int(Perf_metrics.iloc[d,0]-Perf_metrics.iloc[d-1,0])) + + fig,ax1 = plt.subplots() + ax1.plot(exp_time_array_days2,Perf_metrics['S_g_2_0'],label='2.0m depth (meas.)') #2.0m temp depth + ax1.plot(exp_time_array_days2,Perf_metrics['S_g_2_7'],label='2.7m depth (meas.)') #2.0m temp depth + ax1.plot(exp_time_array_days2,Perf_metrics['S_g_3_3'],label='3.3m depth (meas.)') #2.0m temp depth + ax1.set_ylabel('Temperature ($\degree$C)') + ax1.set_title("Ground Temperatures ~1m west of sand store, measured") + + locator = mdates.AutoDateLocator() + ax1.xaxis.set_major_locator(locator) + ax1.xaxis.set_minor_locator(mdates.DayLocator()) + + # #Sets major and minor tick spacing + # loc = ticker.MultipleLocator(base=7.0) + # ax1.xaxis.set_major_locator(loc) + # ax1.xaxis.set_minor_locator(ticker.MultipleLocator(base=1.0)) + + #set the major locator on x-axis to every 6 days + locator.intervald[mdates.DAILY]=[7] #https://matplotlib.org/stable/api/dates_api.html#matplotlib.dates.AutoDateLocator + + #Formatting - make the x-axis labels more consices + #https://matplotlib.org/stable/gallery/ticks_and_spines/date_concise_formatter.html + formatter = mdates.ConciseDateFormatter(locator) + ax1.xaxis.set_major_formatter(formatter) + fig.autofmt_xdate() #automatically makes the x-labels rotate + ax1.set_ylim(0,20) + + + ax1.legend() + +#%%% Plot simulated ground temps near ground Sand TCs +if PerfMetrGrTCs == 1: +# Estimate the x-indices closest to where the ground TC hole might be (~4.2m/14ft south of the house wall) + ground_TC_xindices = np.zeros(shape=[2], dtype=int) + + #location on x-axis of approximately 14 ft (4.2m) south of the house's south wall. The concrete wall of + #the SS is appx 3.52m south of the south wall. The ground TC is therefore 0.25m south of the starting portion of the sand store + ground_TC_xcoord = TCx_coords[0] - 0.25 + + #loop to find two closest TCs + diff = 0 + for i in range(0,len(x_coords2plt)): + diff = ground_TC_xcoord - x_coords2plt[i,0] + #this finds the closest nodes to the real locations of the TCs + if diff <=0: + #then we want those two indices, which are sandwiching the ground_TC_xcoord + ground_TC_xindices[0] = i-1 + ground_TC_xindices[1] = i + break + +# Find the two y-indices nearest to (sandwiching) the ground TC location (1.0m west of Sand Store) + ground_TC_yindices = np.zeros(shape=[2], dtype=int) + + #location on y-axis of west-most TC sand node, minus the widths of the no-power-zone, XPS insulation, and concrete + SS_west_wall_y_coord = TCy_coords[0] - 0.5 - 0.35 - 0.1 + + #location of ground TCs are appx 1m west of SS wall + ground_TC_ycoord = SS_west_wall_y_coord - 1 + + #loop to find two closest TCs + diff = 0 + for j in range(0,len(y_coords2plt)): + diff = ground_TC_ycoord - y_coords2plt[j,0] + #this finds the closest nodes to the real locations of the TCs + if diff <=0: + #then we want those two indices, which are sandwiching the ground_TC_ycoord + ground_TC_yindices[0] = j-1 + ground_TC_yindices[1] = j + break + +# Find the z-indices nearest to the ground TC depths (2.0, 2.7, 3.3 below grade) + ground_TC_zindices = np.zeros(shape=[6], dtype=int) + + #location of ground TCs are appx 2.0, 2.7, 3.3m below grade + ground_TC_zcoords = [2.0, 2.7, 3.3] + + #loop to find two closest TCs + diff = 0 + for TC in range(0,3): + for k in range(0,len(z_coords2plt)): + diff = (z_coords2plt[-1,0] - ground_TC_zcoords[TC]) - z_coords2plt[k,0] + #this finds the closest nodes to the real locations of the TCs + if diff <=0: + #then we want those two indices, which are sandwiching the ground_TC_ycoord + ground_TC_zindices[TC*2] = k-1 + ground_TC_zindices[(TC*2)+1] = k + break + +#%% Plot the simulated temps nearest to the ground TCs + + + Perf_metric_dates = DateArrayFromSimDayHoursArray(Perf_metrics['day'].astype('float'), + exp_time_array_days2[0]) + + color_scheme=['blue','red','orange','black'] + #Plot simulated ground temp data with experimental data + + for k in range(0,3): + #plot one figure for each z-depth ground TC + fig1, ax1 = plt.subplots() + c=0 + #plot temps from both indices around the ground temp y-coordinate + for j in range(0,len(ground_TC_yindices)): + #plot temps from both indices around the ground temp x-coordinate + for i in range(0,len(ground_TC_xindices)): + + ax1.plot(Perf_metric_dates[:],sndVizTemps[77:-1, + ground_TC_xindices[i], + ground_TC_yindices[j], + ground_TC_zindices[(k*2)]], + '--', + color=color_scheme[c]) + ax1.plot(Perf_metric_dates[:],sndVizTemps[77:-1, + ground_TC_xindices[i], + ground_TC_yindices[j], + ground_TC_zindices[(k*2)+1]], + color=color_scheme[c]) + + c=c+1 #next colour + + #plt.plot(sndTemps_datetimes,avg_sim_temps_snd_layer[:,2],marker='*') + + ax1.set_title("Ground temperatures at depth of " + str(ground_TC_zcoords[k]) + 'm') + ax1.set_ylabel("Temperature (\N{degree sign}C)") + ax1.set_ylim(20,40) + ax1.legend(['x=3.5m, y=1.75m, z=6.2m','x=3.5m, y=1.75m, z=6.4m', + 'x=3.8m, y=1.75m, z=6.2m','x=3.8m, y=1.75m, z=6.4m', + 'x=3.5m, y=2.35m, z=6.2m','x=3.5m, y=2.35m, z=6.4m', + 'x=3.8m, y=2.35m, z=6.2m','x=3.8m, y=2.35m, z=6.4m'],loc='best') + locator = mdates.AutoDateLocator() + ax1.xaxis.set_major_locator(locator) + ax1.xaxis.set_minor_locator(mdates.DayLocator()) + + # #Sets major and minor tick spacing + # loc = ticker.MultipleLocator(base=7.0) + # ax1.xaxis.set_major_locator(loc) + # ax1.xaxis.set_minor_locator(ticker.MultipleLocator(base=1.0)) + + #set the major locator on x-axis to every 6 days + locator.intervald[mdates.DAILY]=[7] #https://matplotlib.org/stable/api/dates_api.html#matplotlib.dates.AutoDateLocator + + #Formatting - make the x-axis labels more consices + #https://matplotlib.org/stable/gallery/ticks_and_spines/date_concise_formatter.html + formatter = mdates.ConciseDateFormatter(locator) + ax1.xaxis.set_major_formatter(formatter) + fig.autofmt_xdate() #automatically makes the x-labels rotate + + + # # Plot experimental data on same figure + + # for col in range(0,np.size(avg_expmtl_temps,1)): + # ax1.plot(exp_time_array,avg_expmtl_temps[:,col],marker=',') + + # #rename figure and change legend to include experimental data + # ax1.set_title("Mean temperatures of sand store layers, Run " + run_number) + # ax1.legend(['Sim-bott','Sim-mid','Sim-top','Exp-bott','Exp-mid','Exp-top',],loc='best') + + + # #Documentation on date locators + # #https://matplotlib.org/stable/api/dates_api.html + # locator = mdates.AutoDateLocator() + # ax1.xaxis.set_major_locator(locator) + # ax1.xaxis.set_minor_locator(mdates.DayLocator()) + + # loc = ticker.MultipleLocator(base=2.0) + # ax1.yaxis.set_major_locator(loc) + # ax1.yaxis.set_minor_locator(ticker.MultipleLocator(base=1.0)) + + # # #set the major locator on x-axis to every 6 days + # # locator.intervald[mdates.DAILY]=[6] #https://matplotlib.org/stable/api/dates_api.html#matplotlib.dates.AutoDateLocator + + # #Formatting - make the x-axis labels more consices + # #https://matplotlib.org/stable/gallery/ticks_and_spines/date_concise_formatter.html + # formatter = mdates.ConciseDateFormatter(locator) + # ax1.xaxis.set_major_formatter(formatter) + # fig1.autofmt_xdate() #automatically makes the x-labels rotate + + # #gridlines on + # ax1.grid(True) + + # #Set max and min values for axes + # ax1.set_xlim(exp_time_array[0], exp_time_array[-1]) + # #ax1.set_ylim(16,24) + # ax1.set_ylim(30,58) + + + + + +#%%% Plot a sand TC temp with GHI and Ambient temp (Outside Air) + + # df_temp=func_serr.read_performance_metrics_data('Performance_metrics_data_2021-05-28-to-2021-07-06.dat',dat_headers, main_dir, data_dir) + + # #create empty list + # exp_time_array_days2 = [0]*df_temp.shape[0] + # #put start date in top row of list + # exp_time_array_days2[0] = exp_start_date + # for d in range(1,df_temp.shape[0]): + # exp_time_array_days2[d] = exp_time_array_days2[d-1] + timedelta(days=int(df_temp.iloc[d,0]-df_temp.iloc[d-1,0])) + + # fig,ax1 = plt.subplots() + # ax2 = ax1.twinx() #secondary axis + # ax1.plot(exp_time_array_days2,df_temp['amb_temp'],'+',label='Ambient temp', color='b') #Ambient temp + # ax2.plot(exp_time_array_days2,df_temp['GHI'],marker='+',label='GHI', color='g') #GHI + # ax1.plot(exp_time_array,Exp_sand_temps['TC(A-5-1)'],label='TC(A-5-1)', color='orange') #Top layer TC A-5-1 + # ax1.set_ylabel('Temperature ($\degree$C)') + # ax2.set_ylabel('Global Horiz. Irradiance (MJ/m2)') + # fig.legend() + + + + +#%% 14. PLOT ALL NODES IN Temps, WHEN Master_v5_2.py is interrupted during solver-> A=Temps[:,:,:,0] + # ie. The following are to be used with F9, not as part of this current program + +# #Plot all node temps in y-z cross sectional plane +# A=Temps[:,:,:,0] +# import matplotlib.pyplot as plt +# for x in range(0,A.shape[0]): +# fig1,ax1 = plt.subplots() +# ax1.plot(range(0,A.shape[1]),A[x,:,:]) +# ax1.set_title("y-z cross-sections, node "+ str(x)) +# ax1.set_xlabel('Length of num. domain (m)') +# ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") + +# #To graph cross-sectional temps of one x-z-plane at a time, going up through all z-points +# x=37 +# for z in range(0,A.shape[2]): +# fig1,ax1 = plt.subplots() +# ax1.plot(range(0,A.shape[1]),A[x,:,z]) +# ax1.set_title("y cross-section temps, x node "+ str(x) + ", z node " +str(z)) +# ax1.set_xlabel('Nodes across y domain') +# ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") +# ax1.set_ylim(0,60) +# ax1.minorticks_on() +# ax1.grid(which='both',axis='x',color='0.9') +# ax1.axvspan(11.5,13.5, color='lightgray',alpha=0.5) #shading for concrete nodes (12,13) +# ax1.axvspan(48.5,50.5, color='lightgray',alpha=0.5) #shading for concrete nodes (49,50) +# ax1.axvspan(13.5,15.5, color='pink',alpha=0.5) #shading for XPS nodes (14,15) +# ax1.axvspan(46.5,48.5, color='pink',alpha=0.5) #shading for XPS nodes (47,48) + +# #To graph cross-sectional temps of one x-z-plane at a time, going north to south through all x-points +# z=37 +# for x in range(0,A.shape[0]): +# fig1,ax1 = plt.subplots() +# ax1.plot(range(0,A.shape[1]),A[x,:,z]) +# ax1.set_title("y cross-section temps, x node "+ str(x) + ", z node " +str(z)) +# ax1.set_xlabel('Nodes across y domain') +# ax1.set_ylabel("Temperature (\N{DEGREE SIGN}C)") +# ax1.set_ylim(0,60) +# ax1.minorticks_on() +# ax1.grid(which='both',axis='x',color='0.9') +# ax1.axvspan(11.5,13.5, color='lightgray',alpha=0.5) #shading for concrete nodes (12,13) +# ax1.axvspan(48.5,50.5, color='lightgray',alpha=0.5) #shading for concrete nodes (49,50) +# ax1.axvspan(13.5,15.5, color='pink',alpha=0.5) #shading for XPS nodes (14,15) +# ax1.axvspan(46.5,48.5, color='pink',alpha=0.5) #shading for XPS nodes (47,48) + diff --git a/Master_v5_4.py b/Master_v5_4.py new file mode 100644 index 0000000..dab5bf7 --- /dev/null +++ b/Master_v5_4.py @@ -0,0 +1,2564 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue May 19 16:33:00 2020 +@author: Rebecca + + This is a Master_v5_3.py script based on SVN rev 897. + Changes from SVN rev 897: + -I have removed the module matplotlib + -I changed the double backslash \\ that denotes directories in Windows with the single forward slash / used in Linux. + -I removed the math.remainder and replaced it with the modulus operator (%) since it wasn't introduced until Python3.7' + +by: Rebecca PINTO +2020-01-15""" +# %% Import Packages +import pandas as pd +import numpy as np +import re +import os +import math +from datetime import date +from decimal import Decimal +from pathlib import Path + +#import my functions from other Python files in the same folder +#import func_readSerrData as func_serr + +#TO DEBUG, insert "breakpoint()" where you want the simulation to stop + +#%% 0. GET RUN NUMBER +run_num = int(input("Enter the run number: ")) +# run_num = 999 + +#%% 00. USER DECISIONS +# rewind = 0.25 #corresponds to 720 iterations at dt=300s + +# tolerance for steady-state IC while loop, IC_SS function +tolerance = 0.01 #C + +#%%% Is the north wall BC different? +# north_BC_diff = int(input("Is north wall boundary at an elevated temperature (1=Yes/0=No): ")) +north_BC_diff=0 +#%%% Should all nodes in concrete and XPS layers be output for visualization (1=Yes/0=No)? +all_wall_layer_nodes = 1 + +#%% FUNCTIONS from func_readSerrData + +def read_loops_heat_transfer(filename, dat_headers, main_dir, data_dir): + + import pandas as pd + + #import csv file with and assign headers from dat_headers df + col_headers = dat_headers['Loops_heat_transfer'].dropna() + Serr_loops_heat_transfer = pd.read_table(data_dir.joinpath(filename), + header = None, + names = col_headers, + sep = '\s+', + comment="#") + + + return Serr_loops_heat_transfer + +def read_sand_temps(filename, dat_headers, main_dir, data_dir): + + import pandas as pd + + + + #import csv file with TCs and assign headers from dat_headers df + Serr_sand_temps = pd.read_table(data_dir.joinpath(filename), + header = None, + names = dat_headers['Sand_temperatures'].dropna(), + sep = '\s+', + comment="#") + + + return Serr_sand_temps + +#%% FUNCTIONS +def remove_newline_in_string(string): + #remove the \n somehow inserted into the vectors when str() is used + regex = re.compile(r'\n') + string_out = regex.sub("",string) + return string_out + +def check_mod(rem,divisor): + rem = round(rem,3) + rem = float(rem) + + if rem == divisor: + #if the remainder is equal to the divisor, then the remainder should really + #be zero, and it's a particularity of floating point arithmetic and binary numbers + #and the modulus operator that have made it not equal to zero + rem = 0 + + return rem + +def boundaryConditions(startday,day,node_zcoords,domain_H, T_Nbound, T_SEWbound): + global T_avg, ampl, T_avg_N, ampl_N, north_BC_diff + + DOY = startday + day + + #if in higher years + if DOY > 365*4: + #Translate to day in year 0 (initial year) + DOY = DOY-365*4 + elif DOY > 365*3: + DOY = DOY-365*3 + elif DOY > 365*2: + DOY = DOY-365*2 + elif DOY > 365: + DOY = DOY-365 + + #Boundary temperatures of far-field conditions around sand store --> use Kusada eqn: + #T_soil = T_avg_yearly_temp - ampl*exp(-depth_below_surface/damp_dpth)*cos(2*pi*DOY/365 - depth_below_surface/damp_dpth - phase) + + #Fill side wall matrices with boundary conditions + for k in range(0,len(node_zcoords)): + T_SEWbound[k] = T_avg - ampl*math.exp(-(domain_H-node_zcoords[k])/damp_dpth)*math.cos(2*math.pi*DOY/365 - (domain_H-node_zcoords[k])/damp_dpth - phase) + + #different parameters on north domain boundary + if north_BC_diff == 1: + T_Nbound[k] = T_avg_N - ampl_N*math.exp(-(domain_H-node_zcoords[k])/damp_dpth)*math.cos(2*math.pi*DOY/365 - (domain_H-node_zcoords[k])/damp_dpth - phase) + else: + T_Nbound[k] = T_SEWbound[k] + + + return T_Nbound,T_SEWbound + + +def updateBCs(T_Nbound, T_SEWbound,timestep): + global Temps,U,V,W + #APPLY DIRICHLET BOUNDARY CONDITIONS + for k in range(1,W-1): + #NORTH B.C. applies: + Temps[0,:,k,timestep] = T_Nbound[k] + + #SOUTH WALL B.C. applies: + Temps[U-1,:,k,timestep] = T_SEWbound[k] + + #WEST WALL B.C. applies: + Temps[:,0,k,timestep] = T_SEWbound[k] + + #EAST WALL B.C. applies: + Temps[:,V-1,k,timestep] = T_SEWbound[k] + + + #BOTTOM WALL B.C. applies: + Temps[:,:,0,timestep] = T_SEWbound[0] + #TOP WALL B.C. applies: + Temps[:,:,W-1,timestep] = T_SEWbound[W-1] + + +def print2newfile(filename,day,hour,array2print,dim_X, dim_Y, dim_Z, viz = False): + # global start_YYY, start_MM, start_DD, end_YYY, + #note: the brackets around the temparray allow the array to be treated as a list of only one list, instead of as a list of lists....so it prints out the delimiters between the numbers + #https://stackoverflow.com/questions/42068144/numpy-savetxt-is-not-adding-comma-delimiter + temparray = np.concatenate(([day], [hour], array2print.flatten(order='F'))) + with open(filename,'w') as f: + #writes parameters to file - dimensions of array + f.write("#output_array_dims " + str(dim_X) + "," + str(dim_Y) + "," + str(dim_Z)+"\n") + #writes parameters to file - start and end dates + f.write("#start_and_end_dates_of_simulation " + str(start_YYYY) + "-" + str(start_MM) + "-" + str(start_DD)+ "," + str(end_YYYY) + "-" + str(end_MM) + "-" + str(end_DD)+"\n") + f.write("#Number of days of simulation: " + str(numberofdays_simu) + "\n") + if viz == True: + #writes x-, y-, and z-coordinate vectors of visualized nodes to file + f.write("#x-coords " + remove_newline_in_string(str(node_xcoords)) +"\n") + f.write("#x-viz indices " + remove_newline_in_string(str(VIZ_X_index)) +"\n") + f.write("#y-coords " + remove_newline_in_string(str(node_ycoords)) +"\n") + f.write("#y-viz indices " + remove_newline_in_string(str(VIZ_Y_index)) +"\n") + f.write("#z-coords " + remove_newline_in_string(str(node_zcoords)) +"\n") + f.write("#z-viz indices " + remove_newline_in_string(str(VIZ_Z_index)) +"\n") + #writes x- and y-coordinate vectors of location of TC nodes to file + f.write("#TC_xcoords " + remove_newline_in_string(str(TC_xcoords)) +"\n") + f.write("#TC_ycoords " + remove_newline_in_string(str(TC_ycoords)) +"\n") + #write out initial condition in the form of a flattened array + np.savetxt(f,[temparray],newline=' ', delimiter=',',fmt='%.6f') + + +def append2file(filename,day,hour,array2print): + #note: the brackets around the temparray allow the array to be treated as a list of only one list, instead of as a list of lists....so it prints out the delimiters between the numbers + #https://stackoverflow.com/questions/42068144/numpy-savetxt-is-not-adding-comma-delimiter + temparray = np.concatenate(([day], [hour], array2print.flatten(order='F'))) + with open(filename,'a') as f: + f.write("\n") + np.savetxt(f,[temparray],newline=' ', delimiter=',',fmt='%.6f') + +def IC_SS(dt): + global U,V,W, Temps,Lambda, rho, Cp, node_xcoords, node_ycoords, node_zcoords + global pwr_bott_indx,x_pwr_low_indx, x_pwr_high_indx, y_pwr_low_indx, y_pwr_high_indx + global pwr_mid_indx, pwr_upr_indx, nn, tolerance, rewind, diff, get_out + + import matplotlib.pyplot as plt + #initialize + prev=0 + diff=20 + + fig,ax=plt.subplots() + fig1,ax1=plt.subplots() + fig2,ax2=plt.subplots() + + #while (rewind*10*24*3600/300 > nn):#rewinds of 10 days at dt=300s + + while (diff > tolerance): + nn = nn+1 + print(nn,diff) + + # if math.remainder(nn,100) == 0: + # #save the graphs + # ax.set_title("After " + str(nn) + " iterations, max diff=" + str(diff)) + # fig.savefig(main_dir.joinpath("After " + str(nn) + " iterations_x.png"), dpi=150, bbox_inches='tight') + # ax1.set_title("After " + str(nn) + " iterations") + # fig1.savefig(main_dir.joinpath("After " + str(nn) + " iterations_y.png"), dpi=150, bbox_inches='tight') + # ax2.set_title("After " + str(nn) + " iterations") + # fig2.savefig(main_dir.joinpath("After " + str(nn) + " iterations_z.png"), dpi=150, bbox_inches='tight') + + #initialize new figures for the next 1000 lines + # fig,ax=plt.subplots() + # fig1,ax1=plt.subplots() + # fig2,ax2=plt.subplots() + + # ax.plot(range(0,Temps[:,31,24,0].shape[0]),Temps[:,31,24,0]) + # ax1.plot(range(0,Temps[31,:,24,0].shape[0]),Temps[31,:,24,0]) + # ax2.plot(range(0,Temps[31,31,:,0].shape[0]),Temps[31,31,:,0]) + + + #We are not iterating over boundary walls or the sand store + for i in range(1,U-1): + for j in range(1,V-1): + for k in range(1,W-1): + + #cubic volume of the sand store, don't iterate over it + if (k > indx_bndrz_z[6] and k < indx_bndrz_z[7] and + i > indx_bndrz_x[6] and i < indx_bndrz_x[11] and + j > indx_bndrz_y[6] and j < indx_bndrz_y[11]): + + continue #skip the rest of the for loop + else: + #prev time step temp + old_T = Temps[i,j,k,prev] + + #GETTING TEMPS AROUND NODE FROM PREV TIMESTEP AND LAMBDAs, TO USE IN GOV EQN + North = Temps[i-1,j,k,prev] + South = Temps[i+1,j,k,prev] + West = Temps[i,j-1,k,prev] + East = Temps[i,j+1,k,prev] + Bottom = Temps[i,j,k-1,prev] + Top = Temps[i,j,k+1,prev] + + lmbda_S = Lambda[i+1,j,k] + lmbda_E = Lambda[i,j+1,k] + lmbda_T = Lambda[i,j,k+1] + + #APPLY GOVERNING EQUATION + #variable substitution [for easier reading of GE] + #subscripts _ip1/_im1 mean "i plus 1, i minus 1" + g = dt/ (rho[i,j,k] * Cp[i,j,k]) + lmbda = Lambda[i,j,k] + x_i = node_xcoords[i] + x_ip1 = node_xcoords[i+1] #ip1 -> i plus one + x_im1 = node_xcoords[i-1] #im1 -> i minus one + y_j = node_ycoords[j] + y_jp1 = node_ycoords[j+1] + y_jm1 = node_ycoords[j-1] + z_k = node_zcoords[k] + z_kp1 = node_zcoords[k+1] + z_km1 = node_zcoords[k-1] + + #else use GE without Q_source term + Temps[i,j,k,1] = old_T + g *( (lmbda_S-lmbda)*(South-old_T)/((x_ip1-x_i)**2) + + (2*lmbda/(x_ip1-x_im1))*( ((South-old_T)/(x_ip1-x_i))-((North-old_T)/(x_im1-x_i)) ) + + (lmbda_E-lmbda)*(East-old_T)/((y_jp1-y_j)**2) + + (2*lmbda/(y_jp1-y_jm1))*( ((East-old_T)/(y_jp1-y_j))-((West-old_T)/(y_jm1-y_j)) ) + + (lmbda_T-lmbda)*(Top-old_T)/((z_kp1-z_k)**2) + + (2*lmbda/(z_kp1-z_km1))*( ((Top-old_T)/(z_kp1-z_k))-((Bottom-old_T)/(z_km1-z_k)) ) ) + + + #take difference between the maximum temperatures in the new and old timesteps + diff = np.max(abs(Temps[:,:,:,1]-Temps[:,:,:,0])) + + if np.max(Temps)>1000 or np.min(Temps)<-1000 or np.isinf(Temps).sum()>0 or np.isnan(Temps).sum()>0: + get_out=1 + print("Exited for loop - unstable FDM model") + break + + if diff > 20: + print("yes") + if abs(diff) > 20: + print("abs_yes") + #put temps from just-calculated timestep into "previous" timestep + Temps[:,:,:,0] = Temps[:,:,:,1] + + return diff + +def GoverningEqn(dt,n): + global U,V,W, Temps,Lambda, rho, Cp, node_xcoords, node_ycoords, node_zcoords, solar_to_Sand_W, Sand_to_SH_W + global pwr_bott_indx,x_pwr_low_indx, x_pwr_high_indx, y_pwr_low_indx, y_pwr_high_indx + global pwr_mid_indx, pwr_upr_indx + #We are not iterating over boundary nodes, because their temperatures + #are known + prev = 0 + + for i in range(1,U-1): + for j in range(1,V-1): + for k in range(1,W-1): + #prev time step temp + old_T = Temps[i,j,k,prev] + + #GETTING TEMPS AROUND NODE FROM PREV TIMESTEP AND LAMBDAs, TO USE IN GOV EQN + North = Temps[i-1,j,k,prev] + South = Temps[i+1,j,k,prev] + West = Temps[i,j-1,k,prev] + East = Temps[i,j+1,k,prev] + Bottom = Temps[i,j,k-1,prev] + Top = Temps[i,j,k+1,prev] + + lmbda_S = Lambda[i+1,j,k] + lmbda_E = Lambda[i,j+1,k] + lmbda_T = Lambda[i,j,k+1] + + #APPLY GOVERNING EQUATION + #variable substitution [for easier reading of GE] + #subscripts _ip1/_im1 mean "i plus 1, i minus 1" + g = dt/ (rho[i,j,k] * Cp[i,j,k]) + lmbda = Lambda[i,j,k] + x_i = node_xcoords[i] + x_ip1 = node_xcoords[i+1] #ip1 -> i plus one + x_im1 = node_xcoords[i-1] #im1 -> i minus one + y_j = node_ycoords[j] + y_jp1 = node_ycoords[j+1] + y_jm1 = node_ycoords[j-1] + z_k = node_zcoords[k] + z_kp1 = node_zcoords[k+1] + z_km1 = node_zcoords[k-1] + + + #if node being calculated is within the bottom power slab area + if (k==pwr_bott_indx and (i >= x_pwr_low_indx and i <= x_pwr_high_indx) and (j >= y_pwr_low_indx and j <= y_pwr_high_indx)): + #use GE with Q_source term + Temps[i,j,k,1] = old_T + g *( solar_to_Sand_W[n-1,2] - Sand_to_SH_W[n-1,2] + + (lmbda_S-lmbda)*(South-old_T)/((x_ip1-x_i)**2) + + (2*lmbda/(x_ip1-x_im1))*( ((South-old_T)/(x_ip1-x_i))-((North-old_T)/(x_im1-x_i)) ) + + (lmbda_E-lmbda)*(East-old_T)/((y_jp1-y_j)**2) + + (2*lmbda/(y_jp1-y_jm1))*( ((East-old_T)/(y_jp1-y_j))-((West-old_T)/(y_jm1-y_j)) ) + + (lmbda_T-lmbda)*(Top-old_T)/((z_kp1-z_k)**2) + + (2*lmbda/(z_kp1-z_km1))*( ((Top-old_T)/(z_kp1-z_k))-((Bottom-old_T)/(z_km1-z_k)) ) ) + + #if node being calculated is within the middle power slab area + elif (k==pwr_mid_indx and (i >= x_pwr_low_indx and i <= x_pwr_high_indx) + and (j >= y_pwr_low_indx and j <= y_pwr_high_indx)): + #use GE with Q_source term + Temps[i,j,k,1] = old_T + g *( solar_to_Sand_W[n-1,2] - Sand_to_SH_W[n-1,2] + + (lmbda_S-lmbda)*(South-old_T)/((x_ip1-x_i)**2) + + (2*lmbda/(x_ip1-x_im1))*( ((South-old_T)/(x_ip1-x_i))-((North-old_T)/(x_im1-x_i)) ) + + (lmbda_E-lmbda)*(East-old_T)/((y_jp1-y_j)**2) + + (2*lmbda/(y_jp1-y_jm1))*( ((East-old_T)/(y_jp1-y_j))-((West-old_T)/(y_jm1-y_j)) ) + + (lmbda_T-lmbda)*(Top-old_T)/((z_kp1-z_k)**2) + + (2*lmbda/(z_kp1-z_km1))*( ((Top-old_T)/(z_kp1-z_k))-((Bottom-old_T)/(z_km1-z_k)) ) ) + + #if node being calculated is within the top power slab area + elif (k==pwr_upr_indx and (i >= x_pwr_low_indx and i <= x_pwr_high_indx) + and (j >= y_pwr_low_indx and j <= y_pwr_high_indx)): + #use GE with Q_source term + Temps[i,j,k,1] = old_T + g *( solar_to_Sand_W[n-1,2] - Sand_to_SH_W[n-1,2] + + (lmbda_S-lmbda)*(South-old_T)/((x_ip1-x_i)**2) + + (2*lmbda/(x_ip1-x_im1))*( ((South-old_T)/(x_ip1-x_i))-((North-old_T)/(x_im1-x_i)) ) + + (lmbda_E-lmbda)*(East-old_T)/((y_jp1-y_j)**2) + + (2*lmbda/(y_jp1-y_jm1))*( ((East-old_T)/(y_jp1-y_j))-((West-old_T)/(y_jm1-y_j)) ) + + (lmbda_T-lmbda)*(Top-old_T)/((z_kp1-z_k)**2) + + (2*lmbda/(z_kp1-z_km1))*( ((Top-old_T)/(z_kp1-z_k))-((Bottom-old_T)/(z_km1-z_k)) ) ) + else: + #else use GE without Q_source term + Temps[i,j,k,1] = old_T + g *( (lmbda_S-lmbda)*(South-old_T)/((x_ip1-x_i)**2) + + (2*lmbda/(x_ip1-x_im1))*( ((South-old_T)/(x_ip1-x_i))-((North-old_T)/(x_im1-x_i)) ) + + (lmbda_E-lmbda)*(East-old_T)/((y_jp1-y_j)**2) + + (2*lmbda/(y_jp1-y_jm1))*( ((East-old_T)/(y_jp1-y_j))-((West-old_T)/(y_jm1-y_j)) ) + + (lmbda_T-lmbda)*(Top-old_T)/((z_kp1-z_k)**2) + + (2*lmbda/(z_kp1-z_km1))*( ((Top-old_T)/(z_kp1-z_k))-((Bottom-old_T)/(z_km1-z_k)) ) ) + + +# def partialAvg_of_Temps_in_zDir(XYZarray,min_zIndx,max_zIndx): +# #This function calculates the average of a 3-D chunk of an array, from bottom z index, up to and including the top z-index. +# global indx_bndrz_z, moisture_layer_height + +# temp_sum = 0 +# temp_avg = 0 +# for k in range(min_zIndx,max_zIndx + 1): +# temp_sum = temp_sum + np.mean(XYZarray[:,:,k]) +# temp_avg = temp_sum/len(range(min_zIndx,max_zIndx + 1)) + +# return temp_avg + +# %% Define Directory Paths +#get current directory +main_dir = Path(os.getcwd()) +#define parent directory +parent = main_dir.parents[0] + +#define directory holding input data +data_dir = parent.joinpath('Exp_Data_for_models') +#define directory where graphs will be saved +graph_dir = main_dir.joinpath('Diff_runs') + +#Read in Headers file +dat_headers = pd.read_csv(data_dir.joinpath('Serrano_dat_headers.csv')) + + +# %% 1. IMPORT EXPERIMENTAL DATA FILES +#FILENAMES OF DATA FILES +#(NB! filename1 MUST have a start AND end date, even if they are equal, ie. only one day of data) + +#Excel sheet with different run descriptions +diff_runs_file = graph_dir.joinpath('2021-07-13 Description of different runs.xlsx') +diff_runs = pd.read_excel(diff_runs_file, header=1, index_col=0) + +#check the period column to get the experimental file names accordingly +if diff_runs.Period[run_num].rstrip() == 'cooldown June 2021': + + #%%%% Cooldown (mode 1) + #sand store temperatures file + filename1 = 'Sand_temperatures_2021-05-28-to-2021-07-06.dat' + #loops heat transfer file + filename2 ='Loops_heat_transfer_2021-05-28-to-2021-07-06.dat' + +elif diff_runs.Period[run_num].rstrip() == 'charging Mar 2021': + + #%%%% Charging (mode 2) + #sand store temperatures file + filename1 = 'Sand_temperatures_2021-03-21-to-2021-03-23.dat' + #loops heat transfer file + filename2 ='Loops_heat_transfer_2021-03-21-to-2021-03-23.dat' + +elif diff_runs.Period[run_num].rstrip() == 'discharging Nov 2020': + + #%%%% Discharging (mode 3) + #sand store temperatures file + filename1 = 'Sand_temperatures_2020-11-25-to-2020-11-27.dat' + #loops heat transfer file + filename2 ='Loops_heat_transfer_2020-11-25-to-2020-11-27.dat' + +elif diff_runs.Period[run_num].rstrip() == 'charging + discharging Nov 29': + + #%%%% Charging and Discharging (mode 4) + #sand store temperatures file + filename1 = 'Sand_temperatures_2020-11-28-to-2020-11-30.dat' + #loops heat transfer file + filename2 ='Loops_heat_transfer_2020-11-28-to-2020-11-30.dat' + +elif diff_runs.Period[run_num].rstrip() == 'charging + discharging Nov 23': + + #%%%% Charging and Discharging 2 (mode 5) + #sand store temperatures file + filename1 = 'Sand_temperatures_2020-11-23-to-2020-11-23.dat' + #loops heat transfer file + filename2 ='Loops_heat_transfer_2020-11-23-to-2020-11-23.dat' +elif diff_runs.Period[run_num].rstrip() == 'heating season 2021': + + #%%%% Entire heating season (mode 6) + #sand store temperatures file + filename1 = 'Sand_temperatures_2020-11-19-to-2021-04-30.dat' + #loops heat transfer file + filename2 ='Loops_heat_transfer_2020-11-19-to-2021-04-30.dat' +else: + print("Experimental data files not correctly selected. Please check run number and that the Excel input file is correct.") +#%%% Processing of dates +#Parse out dates between which to run the simulation + +#find year 1 +f1_name, f1_2, f1_3 = re.split('(_20)',filename1) + +st1 = filename1.find('_20') +date_start = filename1[st1+1:st1+11] + +end1 = filename1.find('to-') +date_end = filename1[end1+3:-4] + +#convert to datetime format +start_YYYY, start_MM, start_DD = date_start.split('-') +end_YYYY, end_MM, end_DD = date_end.split('-') + +date_start = date(int(start_YYYY),int(start_MM),int(start_DD)) +date_end = date(int(end_YYYY),int(end_MM),int(end_DD)) + +#Read in Headers file +dat_headers = pd.read_csv(data_dir.joinpath('Serrano_dat_headers.csv')) + +#Read in sand temps data into a DataFrame +Exp_sand_temps = read_sand_temps(filename1, dat_headers, main_dir, data_dir) + +#Read in loops heat transfer data into DataFrame +Exp_heat_transfer_data = read_loops_heat_transfer(filename2, dat_headers, main_dir, data_dir) +#Initialize boolean variable +missing_heat_data = False + +# %%1. Timestep Parameters - CAN BE CHANGED + +dt = 300 #time step, in seconds +output_temp_freq_h = 1 #how often the TC temperatures from the 3-D sand store will be output to the data file, in hours +visualization_output_freq_h = 24 #how often the temperatures from the domain nodes (more nodes than TC temps) will be output to the data file, in hours + +output_temp_freq_int = (output_temp_freq_h * 3600)/dt #frequency that temperatures of the 3-D sand store will be output to the data file, in number of timestep intervals +viz_output_freq_int = (visualization_output_freq_h * 3600)/dt #frequency that temperatures of the domain will be output to the data file, in number of timestep intervals + +#How many days in the simulation? +if (date_start == date_end): + numberofdays_simu = 1 +else: + numberofdays_simu = (date_end - date_start).days + 1 + +tot_time = numberofdays_simu*24*3600 #1h = 3600s +timeIntrvls_per_day = int(24*3600/dt) # number of time step intervals per day +totTimeIntrvls = int(tot_time/dt) # total number of time step intervals + +# %%1. Kusuda parameters - CAN BE CHANGED + +#KUSUDA PARAMETERS for Ottawa(ground temp) +T_avg = 8 #(C) +ampl = 18.507 +#depth from surface (m) +damp_dpth = 3 +phase =0.38 + +#different parameters on north domain boundary due to presence of house +if north_BC_diff == 1: + T_avg_N = 8+4 #(C) + ampl_N = 14 +else: + T_avg_N = 8 + ampl_N = 18.507 + +# %% 2. Process Experimental Sand Temp data into Numpy arrays +""" +This script processes the experimental data, output by Serrano, from a .dat +or .csv file into Matlab matrix form, so that it can be used later + +RECALL!!!! k = 0 is bottom of sand store, while C-layer is bottom-most plane +""" + +#Put time data into a new numpy array (shows the time in hours from the start of the year) +time_h = Exp_sand_temps['time'].copy().to_numpy() + +#rows and columns of sand temp data +[A, B] = Exp_sand_temps.shape + +#create 4-D array to hold temp data +Exp_sand_temps_4d = np.zeros(shape=[6,6,3,A]) + +#create loop to map expt'l data to 4-D matrix +for indx1 in range(0,A): + + indx2 = 1 #index for the number of columns in the matrix. It starts at 1 since the first column is the time-stamp (zero-indexed) + + for k in range (2,-1,-1): + for i in range (0,6): + for j in range (0,6): + Exp_sand_temps_4d[i,j,k,indx1] = Exp_sand_temps.iloc[indx1,indx2] + + #increment the column to be read next time + indx2 = indx2 + 1 + #check that it exists, and that it's not the end of the line + if indx2 >= B: + #break out of the innermost loop + break + + +#take averages of sand temps (lowest index is bottom-most layer of sand store) +avg_expmtl_temps = np.mean(np.mean(Exp_sand_temps_4d,axis=0),axis=0).transpose() + +#convert time in hours into minutes from start of year +time_min = time_h*60 + +#Calculate day of year data starts +DOY = math.floor(time_h[0]/24)+1 +# %% 2. Geometry of Domain +# %%% High-level geometry + +#x-y-z origin lies at north-west corner, bottom of sand store +sand_L = 6 #(m) length of sand store x = 0 m to x = 6 m +sand_W = 6 #(m) width of sand store y = 0 m to y = 6 m +sand_H = 3 #(m) height of sand store, from bottom z = 0 m to top z = 3 m + +XPS_th = 0.35 #(m) spec'd thickness of insulation around sand store (14") +conc_th = 0.10 #(m) spec'd thickness of concrete + +soil_top_th = 1.5 #(m) thickness of soil layer on top of sand store to surface +soil_nrthside_th = 3.5 #(m) thickness of soil layer on north side of sand store to soil near house +soil_th = 10 #(m) thickness of soil layer on S, E, W, and bottom sides of sand store to far-field temp + +no_TC = 0.5 #(m) buffer of 0.5m between sand store wall and 1st TCs (according to Briana Kemery's Visio document:https://chernode.mae.carleton.ca/sbes_svn/CHEeR_instrumentation/Sensor_connections_to_DAQs/Sand_store_TC_and_moisture_meter_placements.vsdx +node_space = (sand_L - 2*no_TC)/5 #Calculating the distance (m) between TCs in the sand store, assuming that there are 6 x 6 TCs equally spaced on each layer, + +no_pwr = 0.3 #(m) length around outside of sand store where the PEX tubes are not layed down (estimated from photos) +pwr_bott = 0.4 #(m) bottom z-value of x-z plane where PEX tubing is laid down +pwr_mid = 1.3 #(m) middle z-value of x-z plane where PEX tubing is laid down +pwr_upr = 2.2 #(m) top z-value of x-z plane where PEX tubing is laid down + +domain_L = soil_th + soil_nrthside_th + conc_th*2 + XPS_th*2 + sand_L #(m) length of numerical domain, x-axis +domain_W = soil_th*2 + conc_th*2 + XPS_th*2 + sand_W#(m) width of numerical domain, y-axis +domain_H = soil_th + soil_top_th + conc_th + XPS_th*2 + sand_H#(m) height of numerical domain, z-axis + + +#%%% Node Placement (varying mesh size) + +#In x (or y-direction), from x=0 to the south. +#0m = north far-field temperature + +""" +The following values can be changed to see the effect of number of nodes on the final solution +""" +dx_soil = 0.3 #(m) dx of nodes in coarse grid horiz soil layer: in between the soil far-field boundary and the "increased node spacing layer" +dx_soilconc_boundary = 0.1 #(m) dx of nodes in fine grid horiz soil layer: between the coarse grid and the concrete layer +soilconc_boundary_th = 0.2 #(m) thickness of fine grid horiz soil layer +nodes_conc = 2 #number of nodes to make up concrete layer (horiz and vertical) +nodes_XPS = 2 #number of nodes to make up XPS layer (horiz and vertical) +nodes_no_pwr = 2 # number of horiz. nodes in sand layer, in the "no_pwr" area +dx_sand = 0.2 #(m) dx of horiz. nodes in sand layer, inside the "pwr" area + +dz_soil_top = 0.3 #(m)dz of nodes in top vert. soil layer +dz_pwr_th = 0.03 #(m) the vertical node spacing in the three power slab layers. Each power slab layer has three equi-sized dz layers + #around and on it, with the node of the middle layer directly on the power slab line (i.e. 0.4m, 1.3m, and 2.2m + #above the sand store bottom) +nodes_dz1_sand = 2 #number of nodes in between pwr_bott layer and bottom of sand store (not including the three nodes around the pwr_bott layer) +nodes_dz2_sand = 4 #number of nodes in between pwr_mid layer and pwr_bott layer (not including the three nodes around each pwr layer) +nodes_dz3_sand = 4 #number of nodes in between pwr_upr layer and pwr_mid layer (not including the three nodes around each pwr layer) +nodes_dz4_sand = 4 #number of nodes in between pwr_upr layer and top of sand store(not including the three nodes around pwr_upr layer) + +#%%% Horizontal Coordinate Vectors +#%%%% SOIL COORDS + +#(Note:the vector of dx's is the same as the vector for dy's, except for +#the length of the northside soil thickness) +#============================================================ +#============================================================ +#number of nodes in coarser soil grid +QN1 = math.floor(round((soil_nrthside_th - soilconc_boundary_th)/dx_soil,3)) +#remainder of coarse soil length needing to be accounted for +RN1 = (Decimal(soil_nrthside_th) - Decimal(soilconc_boundary_th)) % Decimal(dx_soil) +RN1 = check_mod(RN1,dx_soil) + + +#number of nodes in fine soil grid +Q2 = math.floor(round(soilconc_boundary_th/dx_soilconc_boundary,3)) +#remainder of fine soil length needing to be accounted +R2 = Decimal(soilconc_boundary_th) % Decimal(dx_soilconc_boundary) +R2 = check_mod(R2,dx_soilconc_boundary) + +#SUBVECTOR in which to put X-COORDINATES OF NORTH SOIL NODES +#-------------------------------------------------- +Hsoil_coords_n = np.zeros(shape=[1+QN1+Q2]) + +#far-field north boundary x=0 m +Hsoil_coords_n[0] = 0 +#1st node's x-coordinate, in metres: +Hsoil_coords_n[1] = (dx_soil/2) + RN1 + R2 +#starting at index 3, coarse nodes' x-coordinates in metres +for i in range(2,(1+QN1)): + Hsoil_coords_n[i] = Hsoil_coords_n[i-1] + dx_soil + +#1st node in fine grid area: +i = i+1 +Hsoil_coords_n[i] = Hsoil_coords_n[i-1] + dx_soil/2 + dx_soilconc_boundary/2 +i = i+1 +#fine nodes' x-coordinates in metres +for i in range(i,len(Hsoil_coords_n)): + Hsoil_coords_n[i] = Hsoil_coords_n[i-1] + dx_soilconc_boundary + + +#SUBVECTOR in which to put Y-COORDINATES OF WEST SOIL NODES +#-------------------------------------------------- +#number of nodes in coarser soil grid +Q1 = math.floor((soil_th - soilconc_boundary_th)/dx_soil) +#remainder of coarse soil length needing to be accounted for +R1 = (soil_th - soilconc_boundary_th) % dx_soil + +Hsoil_coords_w = np.zeros(shape=[1+Q1+Q2]) + +#far-field west boundary x=0 m +Hsoil_coords_w[0] = 0 +#1st node's y-coordinate, in metres: +Hsoil_coords_w[1] = (dx_soil/2) + R1 + R2 +#starting at index 3, coarse nodes' x-coordinates in metres +for i in range(2,(1+Q1)): + Hsoil_coords_w[i] = Hsoil_coords_w[i-1] + dx_soil + +#1st node in fine grid area: +i = i+1 +Hsoil_coords_w[i] = Hsoil_coords_w[i-1] + dx_soil/2 + dx_soilconc_boundary/2 +i = i+1 +#fine nodes' x-coordinates in metres +for i in range(i,len(Hsoil_coords_w)): + Hsoil_coords_w[i] = Hsoil_coords_w[i-1] + dx_soilconc_boundary + + +#SUBVECTOR in which to put X-COORDINATES OF SOUTH SOIL NODES +#-------------------------------------------------- +Hsoil_coords_s = np.zeros(shape=[1+Q1+Q2]) + +Hsoil_coords_s[(0)] = soil_nrthside_th + conc_th*2 + XPS_th*2 + sand_L + dx_soilconc_boundary/2 +#fine nodes' x-coordinates in metres +for i in range(1,Q2): + Hsoil_coords_s[i] = Hsoil_coords_s[i-1] + dx_soilconc_boundary + +#1st node in coarse grid area: +i = i+1 +Hsoil_coords_s[i] = Hsoil_coords_s[i-1] + dx_soil/2 + dx_soilconc_boundary/2 +i = i+1 +#coarse nodes' x-coordinates in metres +for i in range(i,(len(Hsoil_coords_s)-1)): + Hsoil_coords_s[i] = Hsoil_coords_s[i-1] + dx_soil + +#far-field south boundary x=domain_L +Hsoil_coords_s[(len(Hsoil_coords_s)-1)] = domain_L + + +#SUBVECTOR in which to put Y-COORDINATES OF EAST SOIL NODES +#-------------------------------------------------- +Hsoil_coords_e = np.zeros(shape=[1+Q1+Q2]) + +Hsoil_coords_e[(0)] = soil_th + conc_th*2 + XPS_th*2 + sand_W + dx_soilconc_boundary/2 +#fine nodes' x-coordinates in metres +for i in range(1,Q2): + Hsoil_coords_e[i] = Hsoil_coords_e[i-1] + dx_soilconc_boundary + +#1st node in coarse grid area: +i = i+1 +Hsoil_coords_e[i] = Hsoil_coords_e[i-1] + dx_soil/2 + dx_soilconc_boundary/2 +i = i+1 +#coarse nodes' x-coordinates in metres +for i in range(i,(len(Hsoil_coords_e)-1)): + Hsoil_coords_e[i] = Hsoil_coords_e[i-1] + dx_soil + +#far-field south boundary y=domain_W +Hsoil_coords_e[(len(Hsoil_coords_e)-1)] = domain_W + + +# %%%% CONCRETE COORDS +#============================================================ +Hconc_coords_n = np.zeros(shape=[nodes_conc]) +Hconc_coords_w = np.zeros(shape=[nodes_conc]) +Hconc_coords_s = np.zeros(shape=[nodes_conc]) +Hconc_coords_e = np.zeros(shape=[nodes_conc]) + +#dx for concrete: +Q3 = math.floor(conc_th*100/nodes_conc)/100 +#remainder of concrete length needing to be accounted for +R3 = (conc_th*100 % nodes_conc)/100 + +#1st conc node +Hconc_coords_n[0] = soil_nrthside_th + Q3/2 + R3 +Hconc_coords_w[0] = soil_th + Q3/2 + R3 +Hconc_coords_s[0] = soil_nrthside_th + conc_th + XPS_th*2 + sand_L + Q3/2 +Hconc_coords_e[0] = soil_th + conc_th + XPS_th*2 + sand_W + Q3/2 + +#remainder of nodes +for i in range(1,nodes_conc): + Hconc_coords_n[i] = Hconc_coords_n[i-1] + Q3 + Hconc_coords_w[i] = Hconc_coords_w[i-1] + Q3 + Hconc_coords_s[i] = Hconc_coords_s[i-1] + Q3 + Hconc_coords_e[i] = Hconc_coords_e[i-1] + Q3 + + +#%%%% XPS COORDS +#============================================================ +HXPS_coords_n = np.zeros(shape=[nodes_XPS]) +HXPS_coords_w = np.zeros(shape=[nodes_XPS]) +HXPS_coords_s = np.zeros(shape=[nodes_XPS]) +HXPS_coords_e = np.zeros(shape=[nodes_XPS]) + +#dx for XPS: +Q4 = math.floor(XPS_th*100/nodes_XPS)/100 +#remainder of XPS length needing to be accounted for +R4 = (XPS_th*100 % nodes_XPS)/100 + +#1st XPS node +HXPS_coords_n[0] = soil_nrthside_th + conc_th + Q4/2 + R4 +HXPS_coords_w[0] = soil_th + conc_th + Q4/2 + R4 +HXPS_coords_s[0] = soil_nrthside_th + conc_th + XPS_th + sand_L + Q4/2 +HXPS_coords_e[0] = soil_th + conc_th + XPS_th + sand_W + Q4/2 + +#remainder of nodes +for i in range(1,nodes_XPS): + HXPS_coords_n[i] = HXPS_coords_n[i-1] + Q4 + HXPS_coords_w[i] = HXPS_coords_w[i-1] + Q4 + HXPS_coords_s[i] = HXPS_coords_s[i-1] + Q4 + HXPS_coords_e[i] = HXPS_coords_e[i-1] + Q4 + + +#%%%% SAND, NO-POWER ZONE COORDS (north and west sides) +#============================================================ +Hnopwr_coords_n = np.zeros(shape=[nodes_no_pwr]) +Hnopwr_coords_w = np.zeros(shape=[nodes_no_pwr]) +Hnopwr_coords_s = np.zeros(shape=[nodes_no_pwr]) +Hnopwr_coords_e = np.zeros(shape=[nodes_no_pwr]) + +#dx for no_pwr: +Q5 = math.floor(no_pwr*100/nodes_no_pwr)/100 +#remainder of no_pwr length needing to be accounted for +R5 = (no_pwr*100 % nodes_no_pwr)/100 + +#1st no_pwr node on north and west sides +Hnopwr_coords_n[0] = soil_nrthside_th + conc_th + XPS_th + Q5/2 + R5 +Hnopwr_coords_w[0] = soil_th + conc_th + XPS_th + Q5/2 + R5 + +#remainder of nodes on north and west sides +for i in range(1,nodes_no_pwr): + Hnopwr_coords_n[i] = Hnopwr_coords_n[i-1] + Q5 + Hnopwr_coords_w[i] = Hnopwr_coords_w[i-1] + Q5 + + +#%%%% SAND, POWER ZONE COORDS +#============================================================ +#number of nodes in sand power area, x-dir +Q6 = math.floor(round((sand_L - no_pwr*2)/dx_sand,3)) +#remainder of sand length needing to be accounted for +R6 = (Decimal(sand_L) - Decimal(no_pwr)*2) % Decimal(dx_sand) +R6 = check_mod(R6,dx_sand) + +#number of nodes in sand power area, y-dir +Q7 = math.floor(round((sand_W - no_pwr*2)/dx_sand,3)) +#remainder of sand width needing to be accounted for +R7 = (Decimal(sand_W) - Decimal(no_pwr)*2) % Decimal(dx_sand) +R7 = check_mod(R7,dx_sand) + + +Hsand_coords_x = np.zeros(shape=[Q6]) +Hsand_coords_y = np.zeros(shape=[Q7]) + +Hsand_coords_x[0] = soil_nrthside_th + conc_th + XPS_th + no_pwr + dx_sand/2 + R6 +Hsand_coords_y[0] = soil_th + conc_th + XPS_th + no_pwr + dx_sand/2 + R7 + +#remainder of nodes +for i in range(1,Q6): + Hsand_coords_x[i] = Hsand_coords_x[i-1] + dx_sand + +for j in range(1,Q7): + Hsand_coords_y[(j)] = Hsand_coords_y[(j-1)] + dx_sand + + +#%%%% SAND NO-POWER ZONE COORDS (south and east sides) +#============================================================ +#1st no_pwr node on south and east sides +Hnopwr_coords_s[0] = soil_nrthside_th + conc_th + XPS_th + sand_L - no_pwr + Q5/2 +Hnopwr_coords_e[0] = soil_th + conc_th + XPS_th + sand_W - no_pwr + Q5/2 + +#remainder of nodes on south and east sides +for i in range(1,nodes_no_pwr): + Hnopwr_coords_s[i] = Hnopwr_coords_s[i-1] + Q5 + Hnopwr_coords_e[i] = Hnopwr_coords_e[i-1] + Q5 + + +#%%%% PUTTING X- AND Y- MESH COORDINATES VECTORS TOGETHER +#============================================================ +node_xcoords = np.concatenate((Hsoil_coords_n,Hconc_coords_n,HXPS_coords_n,Hnopwr_coords_n,Hsand_coords_x,Hnopwr_coords_s,HXPS_coords_s,Hconc_coords_s,Hsoil_coords_s)) + +node_ycoords = np.concatenate((Hsoil_coords_w,Hconc_coords_w,HXPS_coords_w,Hnopwr_coords_w,Hsand_coords_y,Hnopwr_coords_e,HXPS_coords_e,Hconc_coords_e,Hsoil_coords_e)) + + +#%%% Vertical Coordinate Vectors +#%%%% SOIL COORDS + +''' +#Reminder: +#number of nodes in coarser soil grid = Q1 +#remainder of coarse soil length needing to be accounted for = R1 + +#number of nodes in fine soil grid = Q2 +#remainder of fine soil length needing to be accounted for = R2 +''' + +#SUBVECTOR in which to put Z-COORDINATES OF BOTTOM SOIL NODES +#-------------------------------------------------- +Vsoil_coords_b = np.zeros(shape=[1+Q1+Q2]) + +#far-field bottom boundary z=0 m +Vsoil_coords_b[0] = 0 +#1st node's z-coordinate, in metres: +Vsoil_coords_b[1] = (dx_soil/2) + R1 + R2 +#starting at index 3, coarse nodes' z-coordinates in metres +for i in range(2,(1+Q1)): + Vsoil_coords_b[i] = Vsoil_coords_b[i-1] + dx_soil + +#1st node in fine grid area: +i = i+1 +Vsoil_coords_b[i] = Vsoil_coords_b[i-1] + dx_soil/2 + dx_soilconc_boundary/2 +i = i+1 +#fine nodes' z-coordinates in metres +for i in range(i,len(Vsoil_coords_b)): + Vsoil_coords_b[i] = Vsoil_coords_b[i-1] + dx_soilconc_boundary + + + +#SUBVECTOR in which to put Z-COORDINATES OF TOP SOIL NODES +#-------------------------------------------------- +#number of nodes in coarser soil grid +Q8 = math.floor((soil_top_th - soilconc_boundary_th)/dz_soil_top) +#remainder of coarse soil length needing to be accounted for +R8 = (soil_top_th - soilconc_boundary_th) % dz_soil_top + +Vsoil_coords_t = np.zeros(shape=[1+Q8+Q2]) + +Vsoil_coords_t[0] = soil_th + conc_th + XPS_th*2 + sand_H + dx_soilconc_boundary/2 +#fine nodes' z-coordinates in metres +for i in range(1,Q2): + Vsoil_coords_t[i] = Vsoil_coords_t[i-1] + dx_soilconc_boundary + +#1st node in coarse grid area: +i = i+1 +Vsoil_coords_t[i] = Vsoil_coords_t[i-1] + dz_soil_top/2 + dx_soilconc_boundary/2 +i = i+1 +#coarse nodes' z-coordinates in metres +for i in range(i,len(Vsoil_coords_t)-1): + Vsoil_coords_t[i] = Vsoil_coords_t[i-1] + dz_soil_top + +#far-field south boundary x=domain_L +Vsoil_coords_t[len(Vsoil_coords_t)-1] = domain_H + + +#%%%% CONCRETE COORDS +#============================================================ +Vconc_coords_b = np.zeros(shape=[nodes_conc]) + +#dx for concrete = Q3: +#remainder of concrete length needing to be accounted for = R3 + +#1st conc node +Vconc_coords_b[0] = soil_th + Q3/2 + R3 + +#remainder of nodes +for i in range(1,nodes_conc): + Vconc_coords_b[i] = Vconc_coords_b[i-1] + Q3 + + +#%%%% XPS COORDS +#============================================================ +VXPS_coords_b = np.zeros(shape=[nodes_XPS]) +VXPS_coords_t = np.zeros(shape=[nodes_XPS]) + +#dx for XPS = Q4 +#remainder of XPS length needing to be accounted for = R4 + +#1st XPS node +VXPS_coords_b[0] = soil_th + conc_th + Q4/2 + R4 +VXPS_coords_t[0] = soil_th + conc_th + XPS_th + sand_H + Q4/2 + +#remainder of nodes +for i in range(1,nodes_XPS): + VXPS_coords_b[i] = VXPS_coords_b[i-1] + Q4 + VXPS_coords_t[i] = VXPS_coords_t[i-1] + Q4 + + +#%%%% SAND COORDS +#============================================================ +#dz in space in between pwr_bott layer and bottom of sand store (not including the three nodes around the pwr_bott layer) +dz1 = (pwr_bott - (dz_pwr_th + dz_pwr_th/2))/nodes_dz1_sand +#initialize vector +Vsand1 = np.zeros(shape=[nodes_dz1_sand+3]) +#1st node coordinate +Vsand1[0] = soil_th + conc_th + XPS_th + dz1/2 +#remainder of nodes +for i in range(1,nodes_dz1_sand): + Vsand1[i] = Vsand1[i-1] + dz1 + +#1st node in "3-node pwr layer" +i = i+1 +Vsand1[i] = Vsand1[i-1] + dz1/2 + dz_pwr_th/2 +i = i+1 +#remainder of nodes in "3-node pwr layer" +for i in range(i,len(Vsand1)): + Vsand1[i] = Vsand1[i-1] + dz_pwr_th + + +#********* + +#dz in space #number of nodes in between pwr_mid layer and pwr_bott layer (not including the three nodes around each pwr layer) +dz2 = ((pwr_mid - (dz_pwr_th + dz_pwr_th/2)) - (pwr_bott + (dz_pwr_th + dz_pwr_th/2)))/nodes_dz2_sand +#initialize vector +Vsand2 = np.zeros(shape=[nodes_dz2_sand+3]) +#1st node coordinate +Vsand2[0] = soil_th + conc_th + XPS_th + pwr_bott + (dz_pwr_th + dz_pwr_th/2) + dz2/2 +#remainder of nodes +for i in range(1,nodes_dz2_sand): + Vsand2[i] = Vsand2[i-1] + dz2 + +#1st node in "3-node pwr layer" +i = i+1 +Vsand2[i] = Vsand2[i-1] + dz2/2 + dz_pwr_th/2 +i = i+1 +#remainder of nodes in "3-node pwr layer" +for i in range(i,len(Vsand2)): + Vsand2[i] = Vsand2[i-1] + dz_pwr_th + + +#********* + +#dz in space #number of nodes in between pwr_upr layer and pwr_mid layer (not including the three nodes around each pwr layer) +dz3 = ((pwr_upr - (dz_pwr_th + dz_pwr_th/2)) - (pwr_mid + (dz_pwr_th + dz_pwr_th/2)))/nodes_dz3_sand +#initialize vector +Vsand3 = np.zeros(shape=[nodes_dz3_sand+3]) +#1st node coordinate +Vsand3[0] = soil_th + conc_th + XPS_th + pwr_mid + (dz_pwr_th + dz_pwr_th/2) + dz3/2 +#remainder of nodes +for i in range(1,nodes_dz3_sand): + Vsand3[i] = Vsand3[i-1] + dz3 + + +#1st node in "3-node pwr layer" +i = i+1 +Vsand3[i] = Vsand3[i-1] + dz3/2 + dz_pwr_th/2 +i = i+1 +#remainder of nodes in "3-node pwr layer" +for i in range(i,len(Vsand3)): + Vsand3[i] = Vsand3[i-1] + dz_pwr_th + + +#********* + +#dz in space #number of nodes in between pwr_upr layer and pwr_mid layer (not including the three nodes around each pwr layer) +dz4 = (sand_H - (pwr_upr + (dz_pwr_th + dz_pwr_th/2)))/nodes_dz4_sand +#initialize vector +Vsand4 = np.zeros(shape=[nodes_dz4_sand]) +#1st node coordinate +Vsand4[0] = soil_th + conc_th + XPS_th + pwr_upr + (dz_pwr_th + dz_pwr_th/2) + dz4/2 +#remainder of nodes +for i in range(1,nodes_dz4_sand): + Vsand4[i] = Vsand4[i-1] + dz4 + + + +#%%%% PUTTING Z- COORDINATE VECTORS TOGETHER +#============================================================ +node_zcoords = np.concatenate((Vsoil_coords_b,Vconc_coords_b,VXPS_coords_b,Vsand1,Vsand2,Vsand3,Vsand4,VXPS_coords_t,Vsoil_coords_t)) +# figure +# plot(ones(length(node_zcoords),1),node_zcoords,'o') + + +#%%% Print out mesh size spacing to screen +#------------------------------------------------------------------------ +print('------------') +print('MESH SIZING:') +print('------------') +print('Soil(coarse) dx,dy,dz =',dx_soil, 'm') +print('Soil(fine) dx,dy,dz =',dx_soilconc_boundary, 'm on',soilconc_boundary_th, 'm') +print('Concrete dx,dy,dz =',Q3, 'm') +print('XPS dx,dy,dz =',Q4, 'm') +print('Sand (no power zone) dx =',Q5, 'm') +print('Sand (power zone) dx,dy =',dx_sand, 'm\n') + +print('Sand dz on 0 m < z <',pwr_bott,'m of sand store =',"%.2f" % dz1, 'm') +print('Sand dz on ',pwr_bott,'m < z <',pwr_mid,'m of sand store =',"%.2f" % dz2, 'm') +print('Sand dz on ',pwr_mid,'m < z <',pwr_upr,'m of sand store =',"%.2f" % dz3, 'm') +print('Sand dz on ',pwr_upr,'m < z <',sand_H,'m of sand store =',"%.2f" % dz4, 'm') + + +#%%%VERTICAL index boundaries (6 layers) +indx_bndrz_z = np.array([0, #bottom soil + len(Vsoil_coords_b)-1, #bottom soil + len(Vsoil_coords_b), #bottom conc + len(Vsoil_coords_b)+len(Vconc_coords_b)-1, #bottom conc + len(Vsoil_coords_b)+len(Vconc_coords_b), #bottom XPS + len(Vsoil_coords_b)+len(Vconc_coords_b)+len(VXPS_coords_b)-1, #bottom XPS + len(Vsoil_coords_b)+len(Vconc_coords_b)+len(VXPS_coords_b), #sand + len(Vsoil_coords_b)+len(Vconc_coords_b)+len(VXPS_coords_b)+len(Vsand1)+len(Vsand2)+len(Vsand3)+len(Vsand4)-1, #sand + len(Vsoil_coords_b)+len(Vconc_coords_b)+len(VXPS_coords_b)+len(Vsand1)+len(Vsand2)+len(Vsand3)+len(Vsand4), #top XPS + len(Vsoil_coords_b)+len(Vconc_coords_b)+len(VXPS_coords_b)+len(Vsand1)+len(Vsand2)+len(Vsand3)+len(Vsand4)+len(VXPS_coords_t)-1,#top XPS + len(Vsoil_coords_b)+len(Vconc_coords_b)+len(VXPS_coords_b)+len(Vsand1)+len(Vsand2)+len(Vsand3)+len(Vsand4)+len(VXPS_coords_t),#top soil + len(node_zcoords)-1])#top soil + +#Indices for the top node (inclusive) of each "moisture layer boundary"...appx middle points between moisture sensor z-coordinates +moisture_layer_height = np.zeros(shape=[3], dtype=int) +moisture_layer_height[0] = int(len(Vsoil_coords_b)+len(Vconc_coords_b)+len(VXPS_coords_b)+len(Vsand1)+ math.floor(nodes_dz2_sand/2) - 1) # -1 to acct for zero-based arrays +moisture_layer_height[1] = int(len(Vsoil_coords_b)+len(Vconc_coords_b)+len(VXPS_coords_b)+len(Vsand1)+len(Vsand2)+math.floor(nodes_dz3_sand/2) - 1) # -1 to acct for zero-based arrays +moisture_layer_height[2] = int(indx_bndrz_z[7]) + +vol_slab = (sand_L - 2*no_pwr)*(sand_W - 2*no_pwr)*dz_pwr_th #(m3) volume of each power source slab + +#Define number of nodes in each direction in numerical model domain + +U = len(node_xcoords) #number of nodes in x-direction +V = len(node_ycoords) #number of nodes in y-direction +W = len(node_zcoords) #number of nodes in z-direction + + +pwr_bott_indx = indx_bndrz_z[5] + nodes_dz1_sand + 2 #[m] k-index for bottom layer of PEX tubing +pwr_mid_indx = indx_bndrz_z[5] + len(Vsand1) + nodes_dz2_sand + 2 #[m] k-index for middle layer of PEX tubing +pwr_upr_indx = indx_bndrz_z[5] + len(Vsand1) + len(Vsand2) + nodes_dz3_sand + 2 #[m] k-index for top layer of PEX tubing + +#%%%HORIZONTAL index boundaries [9 layers] +indx_bndrz_x = np.array([0, #north soil + len(Hsoil_coords_n)-1, #north soil + len(Hsoil_coords_n), #north conc + len(Hsoil_coords_n)+len(Hconc_coords_n)-1, #north conc + len(Hsoil_coords_n)+len(Hconc_coords_n), #north XPS + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)-1, #north XPS + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n), #north sand no power + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+len(Hnopwr_coords_n)-1, #north sand no power + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+len(Hnopwr_coords_n), #sand + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+ + len(Hnopwr_coords_n)+len(Hsand_coords_x)-1,#sand + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+ + len(Hnopwr_coords_n)+len(Hsand_coords_x), #south sand no pwr + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+ + len(Hnopwr_coords_n)+len(Hsand_coords_x)+len(Hnopwr_coords_s)-1, #south sand no pwr + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+ + len(Hnopwr_coords_n)+len(Hsand_coords_x)+len(Hnopwr_coords_s), #south XPS + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+ + len(Hnopwr_coords_n)+len(Hsand_coords_x)+len(Hnopwr_coords_s)+ + len(HXPS_coords_s)-1, #south XPS + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+ + len(Hnopwr_coords_n)+len(Hsand_coords_x)+len(Hnopwr_coords_s)+ + len(HXPS_coords_s), #south conc + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+ + len(Hnopwr_coords_n)+len(Hsand_coords_x)+len(Hnopwr_coords_s)+ + len(HXPS_coords_s)+len(Hconc_coords_s)-1, #south conc + len(Hsoil_coords_n)+len(Hconc_coords_n)+len(HXPS_coords_n)+ + len(Hnopwr_coords_n)+len(Hsand_coords_x)+len(Hnopwr_coords_s)+ + len(HXPS_coords_s)+len(Hconc_coords_s), #south soil + len(node_xcoords)-1]) #south soil + + +indx_bndrz_y = np.array([0, #west soil + len(Hsoil_coords_w)-1, #west soil + len(Hsoil_coords_w), #west conc + len(Hsoil_coords_w)+len(Hconc_coords_w)-1, #west conc + len(Hsoil_coords_w)+len(Hconc_coords_w), #west XPS + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)-1, #west XPS + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w), #west sand no power + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+len(Hnopwr_coords_w)-1, #west sand no power + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w), #sand + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w)+len(Hsand_coords_y)-1,#sand + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w)+len(Hsand_coords_y), #east sand no pwr + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w)+len(Hsand_coords_y)+len(Hnopwr_coords_e)-1, #east sand no pwr + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w)+len(Hsand_coords_y)+len(Hnopwr_coords_e), #east XPS + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w)+len(Hsand_coords_y)+len(Hnopwr_coords_e)+ + len(HXPS_coords_e)-1, #east XPS + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w)+len(Hsand_coords_y)+len(Hnopwr_coords_e)+ + len(HXPS_coords_e), #east conc + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w)+len(Hsand_coords_y)+len(Hnopwr_coords_e)+ + len(HXPS_coords_e)+len(Hconc_coords_e)-1, #east conc + len(Hsoil_coords_w)+len(Hconc_coords_w)+len(HXPS_coords_w)+ + len(Hnopwr_coords_w)+len(Hsand_coords_y)+len(Hnopwr_coords_e)+ + len(HXPS_coords_e)+len(Hconc_coords_e), #east soil + len(node_ycoords)-1]) #east soil + + +#%%%Indices for points within the "power slab" in the x-y plane +x_pwr_low_indx = indx_bndrz_x[8] #inclusive of this index +x_pwr_high_indx = indx_bndrz_x[9]#inclusive of this index +y_pwr_low_indx = indx_bndrz_y[8] #inclusive of this index +y_pwr_high_indx = indx_bndrz_y[9]#inclusive of this index + +#============================================================ +#%%%Vectors with CARTESIAN COORDS OF LAYER BOUNDARIES +#north-south +N_S_layr_bounds = np.array([0, #north boundary + soil_nrthside_th, #soil/conc + soil_nrthside_th + conc_th, #conc/XPS + soil_nrthside_th + conc_th + XPS_th , #XPS/sand + soil_nrthside_th + conc_th + XPS_th + sand_L, #sand/XPS + soil_nrthside_th + conc_th + XPS_th*2 + sand_L, #XPS/conc + soil_nrthside_th + conc_th*2 + XPS_th*2 + sand_L,#conc/soil + soil_nrthside_th + conc_th*2 + XPS_th*2 + sand_L + soil_th]) #south boundary + + +#west-east +W_E_layr_bounds =np.array([0, #west boundary + soil_th, #soil/conc + soil_th + conc_th, #conc/XPS + soil_th + conc_th + XPS_th , #XPS/sand + soil_th + conc_th + XPS_th + sand_W, #sand/XPS + soil_th + conc_th + XPS_th*2 + sand_W, #XPS/conc + soil_th + conc_th*2 + XPS_th*2 + sand_W,#conc/soil + soil_th*2 + conc_th*2 + XPS_th*2 + sand_W]) #east boundary + +#bottom-top +B_T_layr_bounds =np.array([0, #bottom boundary + soil_th, #soil/conc + soil_th + conc_th, #conc/XPS + soil_th + conc_th + XPS_th , #XPS/sand + soil_th + conc_th + XPS_th + sand_H, #sand/XPS + soil_th + conc_th + XPS_th*2 + sand_H, #XPS/soil + soil_th + conc_th + XPS_th*2 + sand_H + soil_top_th]) #top boundary + +#============================================================ +#%%%Vectors holding the x and y coordinates of the TCs in the sand store +TC_xcoords = np.zeros(shape=[6]) +TC_ycoords = np.zeros(shape=[6]) + +TC_xcoords[0] = soil_nrthside_th + conc_th + XPS_th + no_TC +TC_ycoords[0] = soil_th + conc_th + XPS_th + no_TC + +for i in range(1,6): + TC_xcoords[i] = TC_xcoords[i-1] + node_space + TC_ycoords[i] = TC_ycoords[i-1] + node_space + + +#Find the indices of the TC locations +sand_TC_node_xindices = np.zeros(shape=[6], dtype=int) + +for TC in range(0,6): + diff = 0 + prev = 100 + for i in range(0,len(node_xcoords)): + diff = abs(TC_xcoords[TC] - node_xcoords[i]) + #this finds the closest nodes to the real locations of the TCs + if round(diff,5) >= round(prev,5): + #then we want the index represented by [i-1] + sand_TC_node_xindices[TC] = i-1 + break + prev = diff + + +sand_TC_node_yindices = np.zeros(shape=[6], dtype=int) + +for TC in range(0,6): + diff = 0 + prev = 100 + for j in range(0,len(node_ycoords)): + diff = abs(TC_ycoords[TC] - node_ycoords[j]) + #this finds the closest nodes to the real locations of the TCs + if round(diff,5) >= round(prev,5): + #then we want the index represented by [j-1] + sand_TC_node_yindices[TC] = j-1 + break + prev = diff + + +# %% 2. Material Thermal Properties +''' +This script defines the thermal properties of the materials used in the simulation and also builds the 3-D thermal property matrices for the simulation. +''' +#Properties of water +water_Cp = 4.186*1e3 # (J/kgK) specific heat for water +water_rho = 998.2 # (kg/m3) density for water at 293 K https://hypertextbook.com/facts/2007/AllenMa.shtml + +#Properties of dry sand -> also used Hamdhan and Clarke (2010) as a reference +drySand_Cp = 830 #(J/kgK) specific heat for dry sand https://www.engineeringtoolbox.com/specific-heat-capacity-d_391.html +drySand_rho = 1850 #(kg/m3) density for dry sand (personal correspondence with Abdulghader Abdulrahman, Research Manager, Dr.Paul Simms' lab) + +#volume percents of water layers in sand store +vol_percent = np.zeros(shape=[3]) +vol_percent[0] = diff_runs['vol% bott'][run_num] #bottom layer = 33 vol# +vol_percent[1] = diff_runs['vol% mid'][run_num] #middle layer = 11 vol# +vol_percent[2] = diff_runs['vol% top'][run_num] #top layer = 7 vol# + +sand_lambda = np.zeros(shape=[3]) +sand_lambda[0] = diff_runs['λ_bott'][run_num] #(W/mK) bottom layer - thermal conductivity for saturated sand F. Rad (2009) +sand_lambda[1] = diff_runs['λ_mid'][run_num] #(W/mK) middle layer - moist sand https://www.engineeringtoolbox.com/thermal-conductivity-d_429.html +sand_lambda[2] = diff_runs['λ_top'][run_num] #(W/mK) top layer - moist sand https://www.engineeringtoolbox.com/thermal-conductivity-d_429.html + +#Properties of concrete at 300 K (Source - Incropera, 6e, p.939) +conc_Cp = 880 #(J/kgK) specific heat +conc_rho = 2300 # (kg/m3) density +conc_lambda = 1.4 #(W/mK) thermal conductivity + +#Properties of XPS at 297 K(Source - Owens Corning data sheets and papers - see ppt presentation MCG5157 final project for sources) +XPS_Cp_orig = 1500 #(J/kgK) specific heat source - article by BASF staff +XPS_rho_orig = 30 # (kg/m3) density + +#Properties of XPS +vol_percent_H2O_SETface_XPS = diff_runs['SET(B) XPS Vol % H2O'][run_num] +XPS_lambda_SET = diff_runs['SET(B)'][run_num] #(W/mK) thermal conductivity +XPS_Cp_SET = ((water_Cp * (water_rho*vol_percent_H2O_SETface_XPS)) + (XPS_Cp_orig*XPS_rho_orig)) / (water_rho * vol_percent_H2O_SETface_XPS + XPS_rho_orig) + # (Number of Joules to raise sand in m3 by 1K + Joules to raise water in 1m3 by 1K) / (number of kg in 1m3) +XPS_rho_SET = water_rho * vol_percent_H2O_SETface_XPS + XPS_rho_orig + + +#Properties of soil (soil_rhoCp and soil_lambda) around CHEeR house sand +#store +soil_Cp = 750 #(J/kgK) (A. Wills MASc,Table 6.1) +soil_rho = 2800 #(kg/m3) (A. Wills MASc,Table 6.1) +soil_lambda = 2 #(W/mK) thermal conductivity (A. Wills MASc,Table 6.1) (TRY PARAMETRIZING FOR CLAY VALUES ) + +#%% 2. Print user inputs to screen: +print("\nCHECK THESE PARAMETERS \n") +print("Start and end dates: " + str(date_start) + " to " + str(date_end) + "\n") +print("Soil width surrounding sand store: " + str(soil_th) + " m and " + str(soil_nrthside_th) + " m \n") +print("Kusuda T_avg: " + str(T_avg) + " oC \n") +print("Kusuda ampl: " + str(ampl) + "\n") + + +#%% 2. Build matrices of thermal properties +# (from outside layers inwards) + +#Initialize thermal property matrices +Lambda = np.zeros(shape=[U,V,W]) +Cp = np.zeros(shape=[U,V,W]) +rho = np.zeros(shape=[U,V,W]) + +#Soil layer - make entire matrix soil +Lambda[:,:,:] = soil_lambda +Cp[:,:,:] = soil_Cp +rho[:,:,:] = soil_rho + +#Concrete layer +for k in range(indx_bndrz_z[2],indx_bndrz_z[9]+1): + for i in range(indx_bndrz_x[2],indx_bndrz_x[15]+1): + for j in range(indx_bndrz_y[2],indx_bndrz_y[15]+1): + Lambda[i,j,k] = conc_lambda + Cp[i,j,k] = conc_Cp + rho[i,j,k] = conc_rho + + + + +#Insulation layer +for k in range(indx_bndrz_z[4],indx_bndrz_z[9]+1): + for i in range(indx_bndrz_x[4],indx_bndrz_x[13]+1): + for j in range(indx_bndrz_y[4],indx_bndrz_y[13]+1): + Lambda[i,j,k] = XPS_lambda_SET + Cp[i,j,k] = XPS_Cp_SET + rho[i,j,k] = XPS_rho_SET + + +#North insulation layer + +#volume percent of water layers in XPS north face +vpH20_NXPS = diff_runs['North XPS Vol % H2O'][run_num] # just an acronym of "vol_percent_H2O_Nface_XPS"; to make it easier to change values when face is homogeneous +vol_percent_H2O_Nface_XPS = [ vpH20_NXPS,vpH20_NXPS ,vpH20_NXPS ] +northXPS_lambda = diff_runs['North'][run_num] +# northXPS_lambda = [2, 2, 2] + +#1st layer [bottom of insulation to ~halfway between first two PEX layers] +for k in range(indx_bndrz_z[4],moisture_layer_height[0]+1): + for i in range(indx_bndrz_x[4],indx_bndrz_x[5]+1): + for j in range(indx_bndrz_y[4],indx_bndrz_y[13]+1): + Lambda[i,j,k] = northXPS_lambda#[0] + Cp[i,j,k] = ((water_Cp * (water_rho*vol_percent_H2O_Nface_XPS[0])) + (XPS_Cp_orig*XPS_rho_orig)) / (water_rho * vol_percent_H2O_Nface_XPS[0] + XPS_rho_orig) + # (Number of Joules to raise sand in m3 by 1K + Joules to raise water in 1m3 by 1K) / (number of kg in 1m3) + rho[i,j,k] = water_rho * vol_percent_H2O_Nface_XPS[0] + XPS_rho_orig + +#2nd layer [ of 1st layer to ~halfway between next two PEX layers] +for k in range((moisture_layer_height[0]+1),moisture_layer_height[1]+1): + for i in range(indx_bndrz_x[4],indx_bndrz_x[5]+1): + for j in range(indx_bndrz_y[4],indx_bndrz_y[13]+1): + Lambda[i,j,k] = northXPS_lambda#[1] + Cp[i,j,k] = ((water_Cp * (water_rho*vol_percent_H2O_Nface_XPS[1])) + (XPS_Cp_orig*XPS_rho_orig)) / (water_rho * vol_percent_H2O_Nface_XPS[1] + XPS_rho_orig) + rho[i,j,k] = water_rho * vol_percent_H2O_Nface_XPS[1] + XPS_rho_orig + +#3rd layer [ of 2nd layer to top of insulation] +for k in range((moisture_layer_height[1]+1),indx_bndrz_z[9]+1): + for i in range(indx_bndrz_x[4],indx_bndrz_x[5]+1): + for j in range(indx_bndrz_y[4],indx_bndrz_y[13]+1): + Lambda[i,j,k] = northXPS_lambda#[2] + Cp[i,j,k] = ((water_Cp * (water_rho*vol_percent_H2O_Nface_XPS[2])) + (XPS_Cp_orig*XPS_rho_orig)) / (water_rho * vol_percent_H2O_Nface_XPS[2] + XPS_rho_orig) + rho[i,j,k] = water_rho * vol_percent_H2O_Nface_XPS[2] + XPS_rho_orig + + +#West insulation layer + +#volume percent of water layers in XPS west face +vpH20_WXPS =diff_runs['West XPS Vol % H2O'][run_num] # just an acronym of "vol_percent_H2O_Nface_XPS"; to make it easier to change values when face is homogeneous +vol_percent_H2O_Wface_XPS = [ vpH20_WXPS, vpH20_WXPS, vpH20_WXPS ] +westXPS_lambda = diff_runs['West'][run_num] +# westXPS_lambda = [0.6, 0.6, 0.6] + +#1st layer [bottom of insulation to ~halfway between first two PEX layers] +for k in range(indx_bndrz_z[4],moisture_layer_height[0]+1): + #this is starting at index 6 instead of 4 to not overwrite the properties already written in the North-face XPS section above + for i in range(indx_bndrz_x[6],indx_bndrz_x[13]+1): + for j in range(indx_bndrz_y[4],indx_bndrz_y[5]+1): + Lambda[i,j,k] = westXPS_lambda#[0] + Cp[i,j,k] = ((water_Cp * (water_rho*vol_percent_H2O_Wface_XPS[0])) + (XPS_Cp_orig*XPS_rho_orig)) / (water_rho * vol_percent_H2O_Wface_XPS[0] + XPS_rho_orig) + # (Number of Joules to raise sand in m3 by 1K + Joules to raise water in 1m3 by 1K) / (number of kg in 1m3) + rho[i,j,k] = water_rho * vol_percent_H2O_Wface_XPS[0] + XPS_rho_orig + +#2nd layer [ of 1st layer to ~halfway between next two PEX layers] +for k in range((moisture_layer_height[0]+1),moisture_layer_height[1]+1): + for i in range(indx_bndrz_x[6],indx_bndrz_x[13]+1): + for j in range(indx_bndrz_y[4],indx_bndrz_y[5]+1): + Lambda[i,j,k] = westXPS_lambda#[1] + Cp[i,j,k] = ((water_Cp * (water_rho*vol_percent_H2O_Wface_XPS[1])) + (XPS_Cp_orig*XPS_rho_orig)) / (water_rho * vol_percent_H2O_Wface_XPS[1] + XPS_rho_orig) + rho[i,j,k] = water_rho * vol_percent_H2O_Wface_XPS[1] + XPS_rho_orig + +#3rd layer [ of 2nd layer to top of insulation] +for k in range((moisture_layer_height[1]+1),indx_bndrz_z[9]+1): + for i in range(indx_bndrz_x[6],indx_bndrz_x[13]+1): + for j in range(indx_bndrz_y[4],indx_bndrz_y[5]+1): + Lambda[i,j,k] = westXPS_lambda#[2] + Cp[i,j,k] = ((water_Cp * (water_rho*vol_percent_H2O_Wface_XPS[2])) + (XPS_Cp_orig*XPS_rho_orig)) / (water_rho * vol_percent_H2O_Wface_XPS[2] + XPS_rho_orig) + rho[i,j,k] = water_rho * vol_percent_H2O_Wface_XPS[2] + XPS_rho_orig + +#Bottom insulation layer + +#volume percent of water layers in XPS west face +vpH20_BXPS = diff_runs['Bott XPS Vol % H2O'][run_num] # just an acronym of "vol_percent_H2O_Nface_XPS"; to make it easier to change values when face is homogeneous +vol_percent_H2O_Bface_XPS = [ vpH20_BXPS, vpH20_BXPS, vpH20_BXPS ] +bottomXPS_lambda = diff_runs['Bott'][run_num] + +for k in range(indx_bndrz_z[4],indx_bndrz_z[5]+1): + #this is starting at index 6 instead of 4 to not overwrite the properties already written in the North-face XPS section above + for i in range(indx_bndrz_x[6],indx_bndrz_x[13]+1): + for j in range(indx_bndrz_y[6],indx_bndrz_y[13]+1): + Lambda[i,j,k] = bottomXPS_lambda#[0] + Cp[i,j,k] = ((water_Cp * (water_rho*vol_percent_H2O_Bface_XPS[0])) + (XPS_Cp_orig*XPS_rho_orig)) / (water_rho * vol_percent_H2O_Bface_XPS[0] + XPS_rho_orig) + # (Number of Joules to raise sand in m3 by 1K + Joules to raise water in 1m3 by 1K) / (number of kg in 1m3) + rho[i,j,k] = water_rho * vol_percent_H2O_Bface_XPS[0] + XPS_rho_orig + + +#----------------------- + +#1st Sand layer [0m to ~halfway between first two PEX layers] +for k in range(indx_bndrz_z[6],moisture_layer_height[0]+1): + for i in range(indx_bndrz_x[6],indx_bndrz_x[11]+1): + for j in range(indx_bndrz_y[6],indx_bndrz_y[11]+1): + Lambda[i,j,k] = sand_lambda[0] + Cp[i,j,k] = ((water_Cp * (water_rho*vol_percent[0])) + (drySand_Cp*drySand_rho)) / (water_rho * vol_percent[0] + drySand_rho) + # (Number of Joules to raise sand in m3 by 1K + Joules to raise water in 1m3 by 1K) / (number of kg in 1m3) + rho[i,j,k] = water_rho * vol_percent[0] + drySand_rho + + + + +#2nd Sand layer [ of 1st layer to ~halfway between next two PEX layers] +for k in range((moisture_layer_height[0]+1),moisture_layer_height[1]+1): + for i in range(indx_bndrz_x[6],indx_bndrz_x[11]+1): + for j in range(indx_bndrz_y[6],indx_bndrz_y[11]+1): + Lambda[i,j,k] = sand_lambda[1] + Cp[i,j,k] = ((water_Cp * (water_rho*vol_percent[1])) + (drySand_Cp*drySand_rho)) / (water_rho * vol_percent[1] + drySand_rho) + rho[i,j,k] = water_rho * vol_percent[1] + drySand_rho + + + + +#3rd Sand layer [ of 2nd layer to top of sand store] +for k in range((moisture_layer_height[1]+1),moisture_layer_height[2]+1): + for i in range(indx_bndrz_x[6],indx_bndrz_x[11]+1): + for j in range(indx_bndrz_y[6],indx_bndrz_y[11]+1): + Lambda[i,j,k] = sand_lambda[2] + Cp[i,j,k] = ((water_Cp * (water_rho*vol_percent[2])) + (drySand_Cp*drySand_rho)) / (water_rho * vol_percent[2] + drySand_rho) + rho[i,j,k] = water_rho * vol_percent[2] + drySand_rho + +# %% 2. Heat input to Sand Store +# %%% Process heat input data + +#Put heat source data into a new numpy array +solar_to_Sand_W = Exp_heat_transfer_data[['time','sand_H2O']].copy().to_numpy() + +#Specific rate of heat transfer (W/m3) to sand store per unit volume of power source slabs +sp_heat_trnsfr=np.array(solar_to_Sand_W[:,1],ndmin=2).transpose()/(3*vol_slab) + +#append specific rate of heat transfer to array +solar_to_Sand_W = np.append(solar_to_Sand_W,sp_heat_trnsfr,axis=1) + +#%%%% Check for missing data +if (solar_to_Sand_W.shape[0] < totTimeIntrvls): + + #set boolean marker to True + missing_heat_data = True + #make new array with orig data + solar_to_Sand_W_orig_data = solar_to_Sand_W + #re-initialize solar_to_Sand_W array + solar_to_Sand_W = np.zeros(shape=(totTimeIntrvls,solar_to_Sand_W.shape[1])) + + new=0 #initialize second counter variable for loop below + solar_to_Sand_W[0,0] = solar_to_Sand_W_orig_data[0,0] #copy first data point + + for orig in range(0,solar_to_Sand_W_orig_data.shape[0]-1): + + #check gaps between data points to not be not larger than 2 time step intervals (dt), try 1.5 intervals + #If gap is larger than 1.5*dt, fill in missing data + if (solar_to_Sand_W_orig_data[orig+1,0] - solar_to_Sand_W_orig_data[orig,0]) > (1.5*dt/3600): + #number of data points to add (fill in) + fill = round((solar_to_Sand_W_orig_data[orig+1,0] - solar_to_Sand_W_orig_data[orig,0])*3600/dt) + + # print(solar_to_Sand_W_orig_data[orig+1,0], + # solar_to_Sand_W_orig_data[orig,0], + # (solar_to_Sand_W_orig_data[orig+1,0]- solar_to_Sand_W_orig_data[orig,0])*60, + # fill, + # new, + # orig, + # new-orig) + + #fill the missing data with average of points on either side + for i in range(0,fill-1): + new = new + 1 #increment second counter variable + #first column data fill + solar_to_Sand_W[new,0] = solar_to_Sand_W[new-1,0] + (dt/3600) + #second column data fill + solar_to_Sand_W[new,1] = (solar_to_Sand_W_orig_data[orig+1,1] + + solar_to_Sand_W_orig_data[orig,1])/2 + #third column data fill + solar_to_Sand_W[new,2] = (solar_to_Sand_W_orig_data[orig+1,2] + + solar_to_Sand_W_orig_data[orig,2])/2 + + #copy the data point on the "far side of the data gap" to the new array + solar_to_Sand_W[new+1,0] = solar_to_Sand_W_orig_data[orig+1,0] + solar_to_Sand_W[new+1,1] = solar_to_Sand_W_orig_data[orig+1,1] + solar_to_Sand_W[new+1,2] = solar_to_Sand_W_orig_data[orig+1,2] + + else: + #copy original data over to new array + solar_to_Sand_W[new+1,0] = solar_to_Sand_W_orig_data[orig+1,0] + solar_to_Sand_W[new+1,1] = solar_to_Sand_W_orig_data[orig+1,1] + solar_to_Sand_W[new+1,2] = solar_to_Sand_W_orig_data[orig+1,2] + + new = new + 1 #increment second counter variable (regardless of if data has been filled or not) + + + + + + + print('*****************************************************************\n') + print('GAPS IN EXPERIMENTAL SAND STORE HEAT DATA - VERIFY solar_to_Sand_W AND solar_to_Sand_W_orig_data ARRAYS MAKE SENSE!!!!\n') + print('*****************************************************************\n\n') + +else: + print('*****************************************************************\n') + print('NO MISSING SAND STORE HEAT DATA\n') + print('*****************************************************************\n\n') + + +# %% 2. Heat extracted from Sand Store +# %%% Process discharging data + +#Put discharging data into a new numpy array +Sand_to_SH_W = Exp_heat_transfer_data[['time','sand_2_SH']].copy().to_numpy() + +#Specific rate of heat extraction (W/m3) from sand store per unit volume of power source slabs +sp_heat_extr=np.array(Sand_to_SH_W[:,1],ndmin=2).transpose()/(3*vol_slab) + +#append specific rate of heat extraction to array +Sand_to_SH_W = np.append(Sand_to_SH_W,sp_heat_extr,axis=1) + + +#%%%% Check for missing data +if (Sand_to_SH_W.shape[0] < totTimeIntrvls): + + #set boolean marker to True + missing_heat_data = True + #make new array with orig data + Sand_to_SH_W_orig_data = Sand_to_SH_W + #re-initialize Sand_to_SH_W array + Sand_to_SH_W = np.zeros(shape=(totTimeIntrvls,Sand_to_SH_W.shape[1])) + + new=0 #initialize second counter variable for loop below + Sand_to_SH_W[0,0] = Sand_to_SH_W_orig_data[0,0] #copy first data point + + for orig in range(0,Sand_to_SH_W_orig_data.shape[0]-1): + + #check gaps between data points to not be not larger than 2 time step intervals (dt), try 1.5 intervals + #If gap is larger than 1.5*dt, fill in missing data + if (Sand_to_SH_W_orig_data[orig+1,0] - Sand_to_SH_W_orig_data[orig,0]) > (1.5*dt/3600): + #number of data points to add (fill in) + + fill = round((Sand_to_SH_W_orig_data[orig+1,0] - Sand_to_SH_W_orig_data[orig,0])*3600/dt) + + #fill the missing data with average of points on either side + for i in range(0,fill-1): + new = new + 1 #increment second counter variable + #first column data fill + Sand_to_SH_W[new,0] = Sand_to_SH_W[new-1,0] + (dt/3600) + #second column data fill + Sand_to_SH_W[new,1] = (Sand_to_SH_W_orig_data[orig+1,1] + + Sand_to_SH_W_orig_data[orig,1])/2 + #third column data fill + Sand_to_SH_W[new,2] = (Sand_to_SH_W_orig_data[orig+1,2] + + Sand_to_SH_W_orig_data[orig,2])/2 + + #copy the data point on the "far side of the data gap" to the new array + Sand_to_SH_W[new+1,0] = Sand_to_SH_W_orig_data[orig+1,0] + Sand_to_SH_W[new+1,1] = Sand_to_SH_W_orig_data[orig+1,1] + Sand_to_SH_W[new+1,2] = Sand_to_SH_W_orig_data[orig+1,2] + + else: + #copy original data over to new array + Sand_to_SH_W[new+1,0] = Sand_to_SH_W_orig_data[orig+1,0] + Sand_to_SH_W[new+1,1] = Sand_to_SH_W_orig_data[orig+1,1] + Sand_to_SH_W[new+1,2] = Sand_to_SH_W_orig_data[orig+1,2] + + new = new + 1 #increment second counter variable (regardless of if data has been filled or not) + + + + + + + print('*****************************************************************\n') + print('GAPS IN EXPERIMENTAL SAND STORE HEAT EXTRACTION DATA - VERIFY solar_to_Sand_W AND solar_to_Sand_W_orig_data ARRAYS MAKE SENSE!!!!\n') + print('*****************************************************************\n\n') + +else: + print('*****************************************************************\n') + print('NO MISSING SAND STORE HEAT EXTRACTION DATA\n') + print('*****************************************************************\n\n') + + +# %%2. Visualization Output Points + +''' +This script is to get the nodes that I will record to visualize the data to check that the FD eqn is performing as expected +''' + +#%%%1st SOIL LAYER - starting at y, or z = 0 (west or bottom) +#every 60 cm, including points +node_dist = 0.6 #(m) +soil_viz_indices = np.zeros(shape=[math.floor(soil_th/node_dist)+1], dtype=int) + +soil_viz_indices[0] = 0 +for indx in range(1,len(soil_viz_indices)): + prev = 100 + for i in range(0,len(Hsoil_coords_w)): + diff = abs((indx)*node_dist - Hsoil_coords_w[i]) #diff = desired position - current index + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index represented by (i-1) + soil_viz_indices[indx] = i-1 + break + #if next iteration would be outside the soil domain, put the last node in the soil domain there + elif i == (len(Hsoil_coords_w)-1): + soil_viz_indices[indx] = i + + prev = diff + + + + +#%%%NORTH SOIL LAYER +soil_nrthside_viz_indices = np.zeros(shape=[math.floor(soil_nrthside_th/node_dist)+1], dtype=int) + +soil_nrthside_viz_indices[0] = 0 +for indx in range(1,len(soil_nrthside_viz_indices)): + prev = 100 + for i in range(0,len(Hsoil_coords_n)): + diff = abs((indx)*node_dist - Hsoil_coords_n[i]) #diff = desired position - current index + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index represented by (i-1) + soil_nrthside_viz_indices[indx] = i-1 + break + #if next iteration would be outside the soil domain, put the last node in the soil domain there + elif i == len(Hsoil_coords_n)-1: + soil_nrthside_viz_indices[indx] = i + + prev = diff + + + + +#%%%EAST SOIL LAYER - ending at y = domain_W +soil_east_viz_indices = np.zeros(shape=[math.floor(soil_th/node_dist)+1], dtype=int) +soil_east_viz_indices[-1] = indx_bndrz_y[-1] #last index in east soil +c = 0 +for indx in range(len(soil_east_viz_indices)-2,-1,-1): + prev = 100 + c = c+1 + d = 0 + for i in range(len(Hsoil_coords_e)-1,-1,-1): + d = d+1 + diff = abs((W_E_layr_bounds[-1] - c*node_dist) - Hsoil_coords_e[i-1]) #diff = desired position - current index + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the previous node + soil_east_viz_indices[indx] = indx_bndrz_y[-1] - d +1 + break + elif i == 1: #If it's at the end of the loop and nothing has been chosen, choose the last node + soil_east_viz_indices[indx] = indx_bndrz_y[-2] + + prev = diff + + + +#%%%SOUTH SOIL LAYER - ending at x = domain_L +soil_south_viz_indices = np.zeros(shape=[math.floor(soil_th/node_dist)+1], dtype=int) +soil_south_viz_indices[-1] = indx_bndrz_x[-1] #last index in south soil +c = 0 +for indx in range(len(soil_south_viz_indices)-2,-1,-1): + prev = 100 + c = c+1 + d = 0 + for i in range(len(Hsoil_coords_s)-1,-1,-1): + d = d+1 + diff = abs((N_S_layr_bounds[-1] - c*node_dist) - Hsoil_coords_s[i-1]) #diff = desired position - current index + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the previous node + soil_south_viz_indices[indx] = indx_bndrz_x[-1] - d +1 + break + elif i == 1: #If it's at the end of the loop and nothing has been chosen, choose the last node + soil_south_viz_indices[indx] = indx_bndrz_x[-2] + + prev = diff + + + +#%%%TOP SOIL LAYER - ending at x = domain_H +#every 50 cm, including points +node_dist = 0.5 #(m) +soil_top_viz_indices = np.zeros(shape=[math.floor(soil_top_th/node_dist)+1], dtype=int) +c = 0 +soil_top_viz_indices[-1] = indx_bndrz_z[-1] #last index in top soil +for indx in range(len(soil_top_viz_indices)-2,-1,-1): + c = c+1 + d = 0 + prev = 100 + for i in range(len(Vsoil_coords_t)-1,-1,-1): + d = d+1 + diff = abs((B_T_layr_bounds[-1] - c*node_dist) - Vsoil_coords_t[i-1]) #diff = desired position - current index + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the previous node + soil_top_viz_indices[indx] = indx_bndrz_z[-1] - d+1 + break + elif i == 1: #If it's at the end of the loop and nothing has been chosen, choose the last node + soil_top_viz_indices[indx] = indx_bndrz_z[-2] + + prev = diff + + + +#%%%1st CONCRETE LAYER (west or bottom) + +if all_wall_layer_nodes==0: + #one node in concrete layer, closest to centre of layer + conc_wb_viz_index = 0 + prev = 100 + for i in range(0,len(Vconc_coords_b)): + diff = abs((W_E_layr_bounds[1]+W_E_layr_bounds[2])/2 - Vconc_coords_b[i]) + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the previous node + conc_wb_viz_index = indx_bndrz_y[2] + (i-1) + break + elif i == len(Vconc_coords_b)-1: #If it's at the end of the loop and nothing has been chosen, choose the last node + conc_wb_viz_index = indx_bndrz_y[3] + + prev = diff + + + #%%%North CONCRETE LAYER + #one node in concrete layer, closest to centre of layer + conc_n_viz_index = 0 + prev = 100 + for i in range(0,len(Hconc_coords_n)): + diff = abs((N_S_layr_bounds[1]+N_S_layr_bounds[2])/2 - Hconc_coords_n[i]) + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the previous node + conc_n_viz_index = indx_bndrz_x[2] + (i-1) + break + elif i == len(Hconc_coords_n)-1: #If it's at the end of the loop and nothing has been chosen, choose the last node + conc_n_viz_index = indx_bndrz_x[3] + + prev = diff + + + #%%%East CONCRETE LAYER + #one node in concrete layer, closest to centre of layer + conc_e_viz_index = 0 + prev = 100 + d=0 + for i in range(len(Hconc_coords_e)-1,-1,-1): + d = d+1 + diff = abs((W_E_layr_bounds[5]+W_E_layr_bounds[6])/2 - Hconc_coords_e[i]) + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the previous node + conc_e_viz_index = indx_bndrz_y[15] - d +2 + break + elif i == 0: #If it's at the end of the loop and nothing has been chosen, choose the last node + conc_e_viz_index = indx_bndrz_y[14] + + prev = diff + + + #%%%South CONCRETE LAYER + #one node in concrete layer, closest to centre of layer + conc_s_viz_index = 0 + prev = 100 + d=0 + for i in range(len(Hconc_coords_s)-1,-1,-1): + d = d+1 + diff = abs((N_S_layr_bounds[5]+N_S_layr_bounds[6])/2 - Hconc_coords_s[i]) + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the previous node + conc_s_viz_index = indx_bndrz_x[15] - d +2 + break + elif i == 0: #If it's at the end of the loop and nothing has been chosen, choose the last node + conc_s_viz_index = indx_bndrz_x[14] + + prev = diff + + + #%%%1st XPS LAYER + #one node in insulation layer, closest to centre of layer + XPS_wb_viz_index = 0 + prev = 100 + for i in range(0,len(VXPS_coords_b)): + diff = abs((W_E_layr_bounds[2]+W_E_layr_bounds[3])/2 - VXPS_coords_b[i]) + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the + XPS_wb_viz_index = indx_bndrz_y[4] + (i-1) + break + elif i == len(VXPS_coords_b)-1: #If it's at the end of the loop and nothing has been chosen, choose the last node + XPS_wb_viz_index = indx_bndrz_y[5] + + prev = diff + + + #%%%North XPS LAYER + #one node in insulation layer, closest to centre of layer + XPS_n_viz_index = 0 + prev = 100 + for i in range(0,len(HXPS_coords_n)): + diff = abs((N_S_layr_bounds[2] + N_S_layr_bounds[3])/2 - HXPS_coords_n[i]) + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the + XPS_n_viz_index = indx_bndrz_x[4] + (i-1) + break + elif i == len(HXPS_coords_n)-1: #If it's at the end of the loop and nothing has been chosen, choose the last node + XPS_n_viz_index = indx_bndrz_x[5] + + prev = diff + + + #%%%Top XPS LAYER + #one node in insulation layer, closest to centre of layer + XPS_t_viz_index = 0 + prev = 100 + d=0 + for i in range(len(VXPS_coords_t)-1,-1,-1): + d=d+1 + diff = abs((B_T_layr_bounds[4] + B_T_layr_bounds[5])/2 - VXPS_coords_t[i]) + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the previous node + XPS_t_viz_index = indx_bndrz_z[9] - d +2 + break + elif i == 0: #If it's at the end of the loop and nothing has been chosen, choose the last node + XPS_t_viz_index = indx_bndrz_z[8] + + prev = diff + + + #%%%East XPS LAYER + #one node in insulation layer, closest to centre of layer + XPS_e_viz_index = 0 + prev = 100 + d=0 + for i in range(len(HXPS_coords_e)-1,-1,-1): + d=d+1 + diff = abs((W_E_layr_bounds[4] + W_E_layr_bounds[5])/2 - HXPS_coords_e[i]) + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the previous node + XPS_e_viz_index = indx_bndrz_y[13] - d +2 + break + elif i == 0: #If it's at the end of the loop and nothing has been chosen, choose the last node + XPS_e_viz_index = indx_bndrz_y[12] + + prev = diff + + + #%%%South XPS LAYER + #one node in insulation layer, closest to centre of layer + XPS_s_viz_index = 0 + prev = 100 + d=0 + for i in range(len(HXPS_coords_s)-1,-1,-1): + d=d+1 + diff = abs((N_S_layr_bounds[4] + N_S_layr_bounds[5])/2 - HXPS_coords_s[i]) + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the previous node + XPS_s_viz_index = indx_bndrz_x[13] - d +2 + break + elif i == 0: #If it's at the end of the loop and nothing has been chosen, choose the last node + XPS_s_viz_index = indx_bndrz_x[12] + + prev = diff + + +#%%%East SAND "no power" node +#Find the indices that correspond to the nodes in the sand store in the no TC zone +sand_noTCe_viz_index = 0 +prev = 100 +d=0 +for i in range(len(Hnopwr_coords_e)-1,-1,-1): + d=d+1 + diff = abs((W_E_layr_bounds[4] + W_E_layr_bounds[4]-no_pwr)/2 - Hnopwr_coords_e[i]) + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the previous node + sand_noTCe_viz_index = indx_bndrz_y[11] - d + 2 + break + elif i == 0: #If it's at the end of the loop and nothing has been chosen, choose the last node + sand_noTCe_viz_index = indx_bndrz_y[10] + + prev = diff + + +#%%%South SAND "no power" node +#Find the indices that correspond to the nodes in the sand store in the no TC zone +sand_noTCs_viz_index = 0 +prev = 100 +d=0 +for i in range(len(Hnopwr_coords_s)-1,-1,-1): + d=d+1 + diff = abs((N_S_layr_bounds[4] + N_S_layr_bounds[4]-no_pwr)/2 - Hnopwr_coords_s[i]) + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the previous node + sand_noTCs_viz_index = indx_bndrz_x[11] - d + 2 + break + elif i == 0: #If it's at the of the loop and nothing has been chosen, choose the last node + sand_noTCs_viz_index = indx_bndrz_x[10] + + prev = diff + + +#%%%North SAND "no power" node +#Find the indices that correspond to the nodes in the sand store in the no TC zone +sand_noTCn_viz_index = 0 +prev = 100 +for i in range(0,len(Hnopwr_coords_n)): + diff = abs((N_S_layr_bounds[3] + N_S_layr_bounds[3]+no_pwr)/2 - Hnopwr_coords_n[i]) + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the + sand_noTCn_viz_index = indx_bndrz_x[6] + (i-1) + break + elif i == len(Hnopwr_coords_n)-1: #If it's at the of the loop and nothing has been chosen, choose the last node + sand_noTCn_viz_index = indx_bndrz_x[6]+ (i) + + prev = diff + + +#%%%WEST SAND "no power" node +#Find the indices that correspond to the nodes in the sand store in the no TC zone +sand_noTCw_viz_index = 0 +prev = 100 +for i in range(0,len(Hnopwr_coords_w)): + diff = abs((W_E_layr_bounds[3] + W_E_layr_bounds[3]+no_pwr)/2 - Hnopwr_coords_w[i]) + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the + sand_noTCw_viz_index = indx_bndrz_y[6] + (i-1) + break + elif i == len(Hnopwr_coords_n)-1: #If it's at the of the loop and nothing has been chosen, choose the last node + sand_noTCw_viz_index = indx_bndrz_y[6]+ i + + prev = diff + + +#%%%N-S SAND "power" nodes +node_dist = 0.5 +sand_x_viz_indices = np.zeros(shape=[math.floor((sand_L-no_pwr*2)/node_dist) + 2], dtype=int) +#If the number of nodes to be visualized is the same as or greater than the number of FDE nodes... +if len(sand_x_viz_indices) >= len(Hsand_coords_x): + sand_x_viz_indices = np.zeros(shape=(len(Hsand_coords_x)), dtype=int) + sand_x_viz_indices[0] = indx_bndrz_x[8] + for i in range(1,len(sand_x_viz_indices)): + sand_x_viz_indices[i] = sand_x_viz_indices[i-1]+1 +else: + prev = 100 + sand_x_viz_indices[0] = indx_bndrz_x[8] + for indx in range(1,len(sand_x_viz_indices)): + prev = 100 + for i in range(0,len(Hsand_coords_x)): + diff = abs(N_S_layr_bounds[3]+no_pwr + indx*node_dist - Hsand_coords_x[i]) #diff = desired position - current index + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index represented by (i-1) + sand_x_viz_indices[indx] = indx_bndrz_x[8]+ i-1 + break + #if next iteration would be outside the soil domain, put the last node in the soil domain there + elif i == len(Hsand_coords_x)-1: + sand_x_viz_indices[indx] = indx_bndrz_x[9] + + prev = diff + + + + +#%%%W-E SAND "power" nodes +sand_y_viz_indices = np.zeros(shape=[math.floor((sand_W-no_pwr*2)/node_dist) + 2], dtype=int) +#If the number of nodes to be visualized is the same as or greater than the number of FDE nodes... +if len(sand_y_viz_indices) >= len(Hsand_coords_y): + sand_y_viz_indices = np.zeros(shape=(len(Hsand_coords_y)), dtype=int) + sand_y_viz_indices[0] = indx_bndrz_y[8] + for i in range(1,len(sand_y_viz_indices)): + sand_y_viz_indices[i] = sand_y_viz_indices[i-1]+1 +else: + prev = 100 + sand_y_viz_indices[0] = indx_bndrz_y[8] + for indx in range(1,len(sand_y_viz_indices)): + prev = 100 + for i in range(0,len(Hsand_coords_y)): + diff = abs(W_E_layr_bounds[3]+no_pwr + indx*node_dist - Hsand_coords_y[i]) #diff = desired position - current index + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index represented by (i-1) + sand_y_viz_indices[indx] = indx_bndrz_y[8]+ i-1 + break + #if next iteration would be outside the soil domain, put the last node in the soil domain there + elif i == len(Hsand_coords_y)-1: + sand_y_viz_indices[indx] = indx_bndrz_y[9] + + prev = diff + + + + +#%%%VERT SAND nodes +#for the nodes between the bottom of sand store and the "first power" +#layer,pick the one closest to centre +sand_z_viz_indices1 = 0 +prev = 100 +for i in range(0,len(Vsand1)-3): + diff = abs((B_T_layr_bounds[3] + B_T_layr_bounds[3]+ pwr_bott)/2 - Vsand1[i]) + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index representing the + sand_z_viz_indices1 = indx_bndrz_z[6] + (i-1) + break + elif i == (len(Vsand1)-4): #If it's at the of the loop and nothing has been chosen, choose the last node not in the power slab + sand_z_viz_indices1 = indx_bndrz_z[6]+ len(Vsand1)- 4 + + prev = diff + + +#for the nodes between the first power layer and the next layer, keep them 0.3m apart +#of the vector +node_dist = 0.3 +sand_z_viz_indices2 = np.zeros(shape=[math.floor((pwr_mid-pwr_bott)/node_dist)-1], dtype=int) +for indx in range(0,len(sand_z_viz_indices2)): + prev = 100 + for i in range (0,(len(Vsand2)-3)): + diff = abs(B_T_layr_bounds[3]+pwr_bott + (indx+1)*node_dist - Vsand2[i]) #diff = desired position - current index + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index represented by (i-1) + sand_z_viz_indices2[indx] = pwr_bott_indx + i+1 + break + #if next iteration would be outside the soil domain, put the last node in the soil domain there + elif i == len(Vsand2)-4: + sand_z_viz_indices2[indx] = pwr_bott_indx + len(Vsand2)- 2 + + prev = diff + + +#keep same distance for next layer (mid and upr pwr layers) +sand_z_viz_indices3 = np.zeros(shape=[math.floor((pwr_upr-pwr_mid)/node_dist)-1], dtype=int) +for indx in range(0,len(sand_z_viz_indices3)): + prev = 100 + for i in range(0,len(Vsand3)-3): + diff = abs(B_T_layr_bounds[3]+pwr_mid+ (indx+1)*node_dist - Vsand3[i]) #diff = desired position - current index + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index represented by (i-1) + sand_z_viz_indices3[indx] = pwr_mid_indx + i+1 + break + #if next iteration would be outside the soil domain, put the last node in the soil domain there + elif i == len(Vsand3)-4: + sand_z_viz_indices3[indx] = pwr_mid_indx + len(Vsand3)- 2 + + prev = diff + + + + +#keep same distance for last layer (upr layer and top of sand store) +sand_z_viz_indices4 = np.zeros(shape=[math.floor((sand_H-pwr_upr)/node_dist)], dtype=int) +for indx in range(0,len(sand_z_viz_indices4)): + prev = 100 + for i in range(0,len(Vsand4)): + diff = abs(B_T_layr_bounds[3] + pwr_upr+(indx+1)*node_dist - Vsand4[i]) #diff = desired position - current index + #this finds the closest node to the desired position + if round(diff,5) >= round(prev,5): + #then we want the index represented by (i-1) + sand_z_viz_indices4[indx] = pwr_upr_indx + i+1 + break + #if next iteration would be outside the soil domain, put the last node in the soil domain there + elif i == len(Vsand4)-1: + sand_z_viz_indices4[indx] = pwr_upr_indx + len(Vsand4)+1 + + prev = diff + + + + +# Put sand indices together +sand_z_viz_indices = np.concatenate(([indx_bndrz_z[6]],[sand_z_viz_indices1],[pwr_bott_indx],sand_z_viz_indices2,[pwr_mid_indx],sand_z_viz_indices3, [pwr_upr_indx],sand_z_viz_indices4,[indx_bndrz_z[7]])) + + +#%%%PUT AXIS INDICES TOGETHER + +#If the nodes from the concrete and XPS layers should all be output: +if all_wall_layer_nodes == 1: + VIZ_X_index = np.concatenate((soil_nrthside_viz_indices,np.arange(soil_nrthside_viz_indices[-1]+1,sand_noTCn_viz_index),[sand_noTCn_viz_index],sand_x_viz_indices,[sand_noTCs_viz_index],np.arange(sand_noTCs_viz_index+1,soil_south_viz_indices[0]),soil_south_viz_indices)) + + VIZ_Y_index = np.concatenate((soil_viz_indices,np.arange(soil_viz_indices[-1]+1,sand_noTCw_viz_index),[sand_noTCw_viz_index],sand_y_viz_indices,[sand_noTCe_viz_index],np.arange(sand_noTCe_viz_index+1,soil_east_viz_indices[0]), soil_east_viz_indices)) + + VIZ_Z_index = np.concatenate((soil_viz_indices,np.arange(soil_viz_indices[-1]+1,sand_z_viz_indices[0]),sand_z_viz_indices,np.arange(sand_z_viz_indices[-1]+1,soil_top_viz_indices[0]),soil_top_viz_indices)) +#If only one node from the concrete and XPS layers should be output: +else: + VIZ_X_index = np.concatenate((soil_nrthside_viz_indices,[conc_n_viz_index],[XPS_n_viz_index],[sand_noTCn_viz_index],sand_x_viz_indices,[sand_noTCs_viz_index],[XPS_s_viz_index],[conc_s_viz_index],soil_south_viz_indices)) + + VIZ_Y_index = np.concatenate((soil_viz_indices,[conc_wb_viz_index],[XPS_wb_viz_index],[sand_noTCw_viz_index],sand_y_viz_indices,[sand_noTCe_viz_index],[XPS_e_viz_index], [conc_e_viz_index], soil_east_viz_indices)) + + VIZ_Z_index = np.concatenate((soil_viz_indices,[conc_wb_viz_index],[XPS_wb_viz_index],sand_z_viz_indices,[XPS_t_viz_index],soil_top_viz_indices)) + + +#%% 3. FDE SOLVER +#%%% Initial conditions +#%%%% Initialize Variables for FDE Solver + +#Initialize temperature matrix +Temps = np.ones(shape=(U,V,W,2)) #Matrix for temperatures (two time points) + +#Side wall boundaries - initialize matrices +T_SEWbound = np.zeros(shape=(len(node_zcoords))) + +#in case the north wall needs to be made different b/c of its proximity to the CHEeR house, an extra variable was created +T_Nbound = np.zeros(shape=(len(node_zcoords))) + +#Initial avg temperatures of 3 sand store layers +IC_SSlayers = np.zeros(shape=(3)) +IC_SSlayers[0]=avg_expmtl_temps[0,0] #Layer C (bottom) +IC_SSlayers[1]=avg_expmtl_temps[0,1] #Layer B (middle) +IC_SSlayers[2]=avg_expmtl_temps[0,2] #Layer A (top) + +#variable that will hold temps of nodes closest to TCs, and will be printed out for comparison with experimental data +toprintTCs = np.zeros(shape=(len(sand_TC_node_xindices), len(sand_TC_node_yindices), 3)) + +# #if graphs of the entire num. domain are wanted, temps for visualization will be printed out +# toprintViz = np.zeros(shape=(len(VIZ_X_index), len(VIZ_Y_index), len(VIZ_Z_index))) + + +#Total time elapsed in days, and hour of the day +day = 0 +hour = 0 + +#%%%% Set Initial Conditions + +#Set initial temperature conditions for the very first timestep (construct layers from outside in) +n = 0 # first timestep + +#Calculate/update boundary condition temps +[T_Nbound,T_SEWbound] = boundaryConditions(DOY,day,node_zcoords, domain_H, T_Nbound, T_SEWbound) + +#Initial temperature of the concrete and XPS layers -> the average of the sand store average temp and the average of the Kusuda array +XPSConc_IC = np.mean([np.mean(IC_SSlayers), np.mean(T_SEWbound)]) + +#Initialize LOWER Soil layers - set temperatures in LOWER soil layer +#gradient in lower soil layer +gradient_z = (XPSConc_IC - T_SEWbound[0]) / (indx_bndrz_z[1]-indx_bndrz_z[0]) + +for k in range(1,indx_bndrz_z[1]+1): + Temps[:,:,k,n] = T_SEWbound[0] + gradient_z*k + + +#Concrete and XPS layers (entire x-y block from bottom of concrete to top of XPS) +Temps[:,:,indx_bndrz_z[2]:indx_bndrz_z[9]+1,0] = XPSConc_IC + +#Initialize TOP Soil layers - set temperatures in TOP soil layer +#gradient in top soil layer +gradient_z = (XPSConc_IC - T_SEWbound[-1]) / (indx_bndrz_z[11]-indx_bndrz_z[10]) + +for k in range(indx_bndrz_z[11],indx_bndrz_z[10]-1,-1): + Temps[:,:,k,n] = T_SEWbound[-1] + gradient_z*(indx_bndrz_z[11]-k) + + +#Initialize MIDDLE Soil layers - set temperatures in MIDDLE soil layer, at elevation of concrete and XPS +for k in range(indx_bndrz_z[2]-2,indx_bndrz_z[9]+3): + + #gradient in x-direction from north boundary wall and from south boundary wall + gradient_x = (XPSConc_IC - T_SEWbound[k]) / (indx_bndrz_x[1]-indx_bndrz_x[0]) + + for i in range(1,indx_bndrz_x[1]+1): + #for temps nearest to north boundary wall + Temps[i,:,k,0] = T_SEWbound[k] + gradient_x*i + #for temps nearest to south boundary wall + Temps[indx_bndrz_x[-1]-i,:,k,0] = T_SEWbound[k] + gradient_x*i + + #gradient in y-direction from west boundary wal and from east boundary walll + gradient_y = (XPSConc_IC - T_SEWbound[k]) / (indx_bndrz_y[1]-indx_bndrz_y[0]) + + for j in range(1,indx_bndrz_y[1]+1): + #for temps nearest to west boundary wall + Temps[indx_bndrz_x[2]-2:indx_bndrz_x[15]+3,j,k,0] = T_SEWbound[k] + gradient_y*j + #for temps nearest to east boundary wall + Temps[indx_bndrz_x[2]-2:indx_bndrz_x[15]+3,indx_bndrz_y[-1]-j,k,0] = T_SEWbound[k] + gradient_y*j + +#Apply BCs again +updateBCs(T_Nbound, T_SEWbound,0) + +#1st Sand layer (0m to 0.9m from bottom of sand store) +for k in range(indx_bndrz_z[6],moisture_layer_height[0]+1): + for i in range(indx_bndrz_x[6],indx_bndrz_x[11]+1): + for j in range(indx_bndrz_y[6],indx_bndrz_y[11]+1): + Temps[i,j,k,n] = IC_SSlayers[0] + + +#2nd Sand layer (0.9m to 1.8m from bottom of sand store) +for k in range((moisture_layer_height[0]+1),moisture_layer_height[1]+1): + for i in range(indx_bndrz_x[6],indx_bndrz_x[11]+1): + for j in range(indx_bndrz_y[6],indx_bndrz_y[11]+1): + Temps[i,j,k,n] = IC_SSlayers[1] + + +#3rd Sand layer (1.8m to 3.0m from bottom of sand store) +for k in range((moisture_layer_height[1]+1),moisture_layer_height[2]+1): + for i in range(indx_bndrz_x[6],indx_bndrz_x[11]+1): + for j in range(indx_bndrz_y[6],indx_bndrz_y[11]+1): + Temps[i,j,k,n] = IC_SSlayers[2] + +# #Copy the initial boundary conditions to next timestep (since now they are "ones" and they are +# #only updated every week in the simulation) +# updateBCs(T_Nbound, T_SEWbound,1) + +#Copy the initial conditions to next timestep (since the IC_SS function will run next, and doesn't iterate +#over the middle of the sand store) +Temps[:,:,:,1] = Temps[:,:,:,0] + + + +#%%%% Print initial conditions to file + +#Copy temperature values that need to be output to compare to experimental data +for i in range(0,np.size(toprintTCs,0)): + for j in range(0,np.size(toprintTCs,1)): + toprintTCs[i,j,0] = Temps[sand_TC_node_xindices[i],sand_TC_node_yindices[j],pwr_bott_indx,n] + toprintTCs[i,j,1] = Temps[sand_TC_node_xindices[i],sand_TC_node_yindices[j],pwr_mid_indx,n] + toprintTCs[i,j,2] = Temps[sand_TC_node_xindices[i],sand_TC_node_yindices[j],pwr_upr_indx,n] + + + +# #Copy intial temperature values to be visualized +# for i in range(0,np.size(toprintViz,0)): +# for j in range(0,np.size(toprintViz,1)): +# for k in range(0,np.size(toprintViz,2)): +# toprintViz[i,j,k] = Temps[VIZ_X_index[i], VIZ_Y_index[j], VIZ_Z_index[k],n] + + +#Output the first timestep of data (initial condition) to output files +#This file will hold temperature values that need to be compared to experimental data +fileID = 'r' + str(run_num) + '_1Sand_TC_Temps.txt' +print2newfile(fileID,day,hour,toprintTCs,len(sand_TC_node_xindices),len(sand_TC_node_yindices),3) + + +# #This file will hold temperature values that will be visualized +# fileID2='r' + str(run_num) + '_2Sand_Viz_Temps.txt' +# temparray = np.concatenate(([day], [hour], toprintViz.flatten(order='F'))) +# print2newfile(fileID2,day,hour,toprintViz,len(VIZ_X_index),len(VIZ_Y_index),len(VIZ_Z_index),True) + + + +#%%% Run function to get initial conditions to steady state (pre-conditioning period) +#create variable to exit the solver in case FDM model becomes unstable +get_out=0 + +#iterations variable, nn +nn=0 +diff_tol = IC_SS(dt) #FUNCTION + +print("IC_SS iterations: nn = " + str(nn) + ", diff = " + str(diff_tol)) + +#%%%% Print steady state initial conditions to file +if get_out==0: #if previous loop had get_out==1, skip this + + n=0 + #Copy temperature values that need to be output to compare to experimental data + for i in range(0,np.size(toprintTCs,0)): + for j in range(0,np.size(toprintTCs,1)): + toprintTCs[i,j,0] = Temps[sand_TC_node_xindices[i],sand_TC_node_yindices[j],pwr_bott_indx,n] + toprintTCs[i,j,1] = Temps[sand_TC_node_xindices[i],sand_TC_node_yindices[j],pwr_mid_indx,n] + toprintTCs[i,j,2] = Temps[sand_TC_node_xindices[i],sand_TC_node_yindices[j],pwr_upr_indx,n] + + + + # #Copy intial temperature values to be visualized + # for i in range(0,np.size(toprintViz,0)): + # for j in range(0,np.size(toprintViz,1)): + # for k in range(0,np.size(toprintViz,2)): + # toprintViz[i,j,k] = Temps[VIZ_X_index[i], VIZ_Y_index[j], VIZ_Z_index[k],n] + + #Output the data with a fake hour value of 0.001h (something close to zero, since we already have values output at t=0s) + append2file(fileID, day, 0.001, toprintTCs) + # append2file(fileID2, day, 0.001, toprintViz) + +#%%% Solver +for n in range(1,totTimeIntrvls): + if get_out==1: #if previous loop had get_out==1 + break + + #Total time elapsed in days, and hour of the day + day = n*dt / 86400 + hour = (day - math.floor(day))*24 + + #%%%% Boundary Conditions Update + if day%7 == 0: #If modulus (remainder) of division of day by 7 is zero + #Calculate/update boundary condition temps every week + [T_Nbound,T_SEWbound] = boundaryConditions(DOY,day,node_zcoords,domain_H, T_Nbound, T_SEWbound) + + #Apply newly calculated B.C.'s to boundary temps + updateBCs(T_Nbound, T_SEWbound,1) + + + + #%%%% Governing Eqn + GoverningEqn(dt,n) + + + #%%%% Print temperatures to output file at requested timestep interval + if n%output_temp_freq_int == 0: #checking modulus (remainder) of operation + #first check if model is stable, if it's not, exit simulation + if np.max(Temps)>1000 or np.min(Temps)<-1000 or np.isinf(Temps).sum()>0 or np.isnan(Temps).sum()>0: + get_out=1 + print("Exited for loop - unstable FDM model") + break + + #Copy temperature values that need to be output to compare to experimental data + for i in range(0,np.size(toprintTCs,0)): + for j in range(0,np.size(toprintTCs,1)): + toprintTCs[i,j,0] = Temps[sand_TC_node_xindices[i],sand_TC_node_yindices[j],pwr_bott_indx,1] + toprintTCs[i,j,1] = Temps[sand_TC_node_xindices[i],sand_TC_node_yindices[j],pwr_mid_indx,1] + toprintTCs[i,j,2] = Temps[sand_TC_node_xindices[i],sand_TC_node_yindices[j],pwr_upr_indx,1] + #write temperatures from that timestep to output file + append2file(fileID, day, hour, toprintTCs) + + #print out last timestep of simulation + elif n == (totTimeIntrvls - 1): + #Copy temperature values that need to be output to compare to experimental data + for i in range(0,np.size(toprintTCs,0)): + for j in range(0,np.size(toprintTCs,1)): + toprintTCs[i,j,0] = Temps[sand_TC_node_xindices[i],sand_TC_node_yindices[j],pwr_bott_indx,1] + toprintTCs[i,j,1] = Temps[sand_TC_node_xindices[i],sand_TC_node_yindices[j],pwr_mid_indx,1] + toprintTCs[i,j,2] = Temps[sand_TC_node_xindices[i],sand_TC_node_yindices[j],pwr_upr_indx,1] + #write temperatures from that timestep to output file + append2file(fileID, day, hour, toprintTCs) + + + # #save variables to output viz file at the requested timestep interval + # if (n%viz_output_freq_int == 0 and n!=0): #modulus + # #Print out temperature values that will be visualized + # for i in range(0,np.size(toprintViz,0)): + # for j in range(0,np.size(toprintViz,1)): + # for k in range(0,np.size(toprintViz,2)): + # toprintViz[i,j,k] = Temps[VIZ_X_index[i], VIZ_Y_index[j], VIZ_Z_index[k],1] + # #print temps to file + # append2file(fileID2, day, hour, toprintViz) + + # elif n == (totTimeIntrvls - 1): #print out last timestep + # #Print out temperature values that will be visualized + # for i in range(0,np.size(toprintViz,0)): + # for j in range(0,np.size(toprintViz,1)): + # for k in range(0,np.size(toprintViz,2)): + # toprintViz[i,j,k] = Temps[VIZ_X_index[i], VIZ_Y_index[j], VIZ_Z_index[k],1] + # #print temps to file + # append2file(fileID2, day, hour, toprintViz) + +#%%%%Prepare for next iteration + #Put temps from just-calculated timestep into "previous" timestep + Temps[:,:,:,0]= Temps[:,:,:,1] + + +#%% 4. PARAMETER OUTPUTS TO FILE + +if get_out ==1: + filename ='r' + str(run_num) +'_SimParameters_MODEL UNSTABLE.txt' + with open(filename,'w') as f: + f.write("Run number "+ str(run_num)+ "\n\n") + f.write("Did not finish simulation, model was unstable. Solver loop exited.") + + f.write("Run number "+ str(run_num)+ "\n\n") + + f.write("Number of pre-conditioning iterations: "+ str(nn)+ "\n") + f.write("Largest temperature difference between final pre-conditioning iterations: "+ str(diff_tol)+ "\n\n") + + f.write(filename1 + "\n") + f.write(filename2 + "\n\n") + + if north_BC_diff == 1: + f.write("North wall BC different? YES \n") + f.write("T_avg_N = " + str(T_avg_N) + "\n") + f.write("ampl_N = " + str(ampl_N) + "\n") + else: + f.write("North wall BC different? NO \n") + f.write("T_avg = " + str(T_avg) + "\n") + f.write("ampl = " + str(ampl) + "\n\n") + + f.write("water_Cp = "+ str(water_Cp)+ "\n") + f.write("water_rho = "+ str(water_rho)+ "\n") + f.write("drySand_Cp = "+ str(drySand_Cp)+ "\n") + f.write("drySand_rho = "+ str(drySand_rho)+ "\n") + f.write("vol_percent = "+ str(vol_percent)+ "\n") + f.write("sand_lambda = "+ str(sand_lambda)+ "\n") + f.write("conc_Cp = "+ str(conc_Cp)+ "\n") + f.write("conc_rho = "+ str(conc_rho)+ "\n") + f.write("conc_lambda = "+ str(conc_lambda)+ "\n") + + f.write("XPS_Cp_orig = "+ str(XPS_Cp_orig)+ "\n") + f.write("XPS_rho_orig = "+ str(XPS_rho_orig)+ "\n") + f.write("XPS_lambda_SET = "+ str(XPS_lambda_SET)+ "\n") + + f.write("vol_percent_H2O_SETface_XPS = "+ str(vol_percent_H2O_SETface_XPS)+ "\n") + f.write("XPS_Cp_SET = "+ str(XPS_Cp_SET)+ "\n") + f.write("XPS_rho_SET = "+ str(XPS_rho_SET)+ "\n") + f.write("vol_percent_H2O_Nface_XPS = "+ str(vol_percent_H2O_Nface_XPS)+ "\n") + f.write("northXPS_lambda = "+ str(northXPS_lambda)+ "\n") + f.write("vol_percent_H2O_Wface_XPS = "+ str(vol_percent_H2O_Wface_XPS)+ "\n") + f.write("westXPS_lambda = "+ str(westXPS_lambda)+ "\n") + f.write("vol_percent_H2O_Bface_XPS = "+ str(vol_percent_H2O_Bface_XPS)+ "\n") + f.write("bottomXPS_lambda = "+ str(bottomXPS_lambda)+ "\n") + f.write("XPS Cp, north wall = "+ str([Cp[indx_bndrz_x[4],indx_bndrz_y[4],indx_bndrz_z[4]], + Cp[indx_bndrz_x[4],indx_bndrz_y[4],(moisture_layer_height[0]+1)], + Cp[indx_bndrz_x[4],indx_bndrz_y[4],(moisture_layer_height[1]+1)]]) + "\n") + f.write("XPS rho, north wall = "+ str([rho[indx_bndrz_x[4],indx_bndrz_y[4],indx_bndrz_z[4]], + rho[indx_bndrz_x[4],indx_bndrz_y[4],(moisture_layer_height[0]+1)], + rho[indx_bndrz_x[4],indx_bndrz_y[4],(moisture_layer_height[1]+1)]]) + "\n") + f.write("XPS Cp, west wall = "+ str([Cp[indx_bndrz_x[6],indx_bndrz_y[4],indx_bndrz_z[4]], + Cp[indx_bndrz_x[6],indx_bndrz_y[4],(moisture_layer_height[0]+1)], + Cp[indx_bndrz_x[6],indx_bndrz_y[4],(moisture_layer_height[1]+1)]]) + "\n") + f.write("XPS rho, west wall = "+ str([rho[indx_bndrz_x[6],indx_bndrz_y[4],indx_bndrz_z[4]], + rho[indx_bndrz_x[6],indx_bndrz_y[4],(moisture_layer_height[0]+1)], + rho[indx_bndrz_x[6],indx_bndrz_y[4],(moisture_layer_height[1]+1)]]) + "\n") + f.write("XPS Cp, bottom face = "+ str(Cp[indx_bndrz_x[6],indx_bndrz_y[6],indx_bndrz_z[4]])+ "\n") + f.write("XPS rho, bottom face = "+ str(rho[indx_bndrz_x[6],indx_bndrz_y[6],indx_bndrz_z[4]]) + "\n") + + f.write("soil_Cp = "+ str(soil_Cp)+ "\n") + f.write("soil_rho = "+ str(soil_rho)+ "\n") + f.write("soil_lambda = "+ str(soil_lambda) + "\n\n") + + f.write("simulation timestep, dt = "+ str(dt)+ "\n") + f.write("output_temp_freq_h = "+ str(output_temp_freq_h)+ "\n") + f.write("visualization_output_freq_h = "+ str(visualization_output_freq_h) +"\n\n") + + f.write("sand_L = "+ str(sand_L)+ "\n") + f.write("sand_W = "+ str(sand_W)+ "\n") + f.write("sand_H = "+ str(sand_H)+ "\n") + f.write("XPS_th = "+ str(XPS_th)+ "\n") + f.write("conc_th = "+ str(conc_th)+ "\n") + f.write("soil_top_th = "+ str(soil_top_th)+ "\n") + f.write("soil_nrthside_th = "+ str(soil_nrthside_th)+ "\n") + f.write("soil_th = "+ str(soil_th)+ "\n") + f.write("no_TC = "+ str(no_TC)+ "\n") + f.write("node_space = "+ str(node_space)+ "\n") + f.write("no_pwr = "+ str(no_pwr)+ "\n") + f.write("pwr_bott = "+ str(pwr_bott)+ "\n") + f.write("pwr_mid = "+ str(pwr_mid)+ "\n") + f.write("pwr_upr = "+ str(pwr_upr)+ "\n\n") + + f.write("------------"+ "\n") + f.write("MESH SIZING:"+ "\n") + f.write("------------"+ "\n") + f.write("Soil(coarse) dx,dy,dz = " + str(dx_soil) + "m" + "\n") + f.write("Soil(fine) dx,dy,dz = " + str(dx_soilconc_boundary)+ "m on " + str(soilconc_boundary_th)+ "m"+ "\n") + f.write("Concrete dx,dy,dz = "+ str(Q3) + "m"+ "\n") + f.write("XPS dx,dy,dz = "+ str(Q4) + "m"+ "\n") + f.write("Sand (no power zone) dx = "+ str(Q5) + "m"+ "\n") + f.write("Sand (power zone) dx,dy = "+ str(dx_sand) + "m\n"+ "\n") + + f.write("Sand dz on 0