# -*- 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)