#!/usr/bin/python # -*- coding: ascii -*- # from __future__ import print_function, division import os, sys from os import path import numpy as np import CoolProp import glob from warnings import warn from time import clock import CoolProp.constants from CoolProp.CoolProp import PropsSI,generate_update_pair,get_parameter_index,set_debug_level,get_phase_index from CoolProp import AbstractState as State import matplotlib as mpl import matplotlib.pyplot as plt import matplotlib.patches as mpatches import copy from itertools import cycle from matplotlib import gridspec, ticker from jopy.dataPlotters import roundList, range_brace try: from jopy.dataPlotters import BasePlotter bp = BasePlotter() except: bp = None # The basic settings for he plots xypoints = 1000 loops = 1 repeat = 1 runs = 0 maxruns = 100 plot = True calc = False check = False folder = "dataDESKTOP" figures = "figures" #np.random.seed(1984) fluids = ["CO2","Pentane","R134a","Water","Air","LiBr-0%"] #glskeys = [r"\glsentryshort{co2}",r"\glsentryshort{pentane}",r"\glsentryshort{r134a}",r"\glsentryshort{water}",r"\glsentryshort{air}",r"\glsentryshort{libr} \SI{0}{\percent}"] #glskeys = [r"\ce{CO2}",r"n-Ppentane",r"R134a",r"Water",r"Air",r"\glsentryshort{libr} \SI{0}{\percent}"] repList = [] #for i in range(len(fluids)): # repList.append(fluids[i]) # repList.append(glskeys[i]) backends = ["INCOMP","HEOS","REFPROP"] #repList.append("HEOS") #repList.append(r"\glsentryshort{cp}") #repList.append("REFPROP") #repList.append(r"\glsentryshort{rp}") #CoolProp.CoolProp.set_debug_level(51) pStr = path.dirname(path.abspath(__file__)) fStr = path.splitext(path.basename(__file__))[0] def getFolderName(): folderName = path.join(pStr,folder) if not path.isdir(folderName): print("Creating data directory "+folderName) os.makedirs(folderName) return folderName def getFigureFolder(): folderName = path.join(pStr,figures) if not path.isdir(folderName): print("Creating data directory "+folderName) os.makedirs(folderName) return folderName repList.append("TimeComp-") repList.append("chapters/FluidProperties/"+path.basename(getFigureFolder())+"/TimeComp-") def getFileName(qualifiers=[]): fileName = path.join(getFolderName(),"-".join(qualifiers)) return fileName # Some file handling def loadNpzData(backend,fluid): dicts = {} globber = getFileName([backend,fluid])+'_[0-9][0-9][0-9].npz' for fname in glob.glob(globber): dataDict = dict(np.load(fname)) dicts[str(dataDict["name"])] = dataDict #if len(dicts)<1: # #print("No readable file found for {0}".format(globber)) # dataDict = dict(name=str(0).zfill(3)) # dicts[str(dataDict["name"])] = dataDict return dicts def saveNpzData(backend,fluid,dicts,start=0,stop=-1): keys = dicts.keys() keys.sort() for k in keys[start:stop]: data = dicts[k] fname = getFileName([backend,fluid])+'_{0}.npz'.format(str(data['name']).zfill(3)) np.savez(fname, **data) return True def splitFluid(propsfluid): fld = propsfluid.split("::") if len(fld)==2: backend = fld[0] fld = fld[1] else: backend = None fld = fld[0] fld = fld.split("-") if len(fld)==2: conc = float(fld[1].strip('%'))/100.0 fld = fld[0] else: conc = None fld = fld[0] return backend,fld,conc def getInpList(backend): if backend=="HEOS": return ["DT","HP"] elif backend=="REFPROP": return ["DT","HP"] elif backend=="INCOMP": return ["PT","HP"] else: raise ValueError("Unknown backend.") def getOutList(inp=None): if inp == "HP": return [["Tmax"],["D"],["S"],["T"],["D","S","T"]] elif inp == "DT": return [["Tmax"],["H"],["P"],["S"],["H","P","S"]] elif inp == "PT": return [["Tmax"],["H"],["D"],["S"],["H","D","S"]] else: raise ValueError("Unknown inputs.") #def getPhaseString(iPhase): # for i in range(11): # if getPhaseString(i)==sPhase: # return i # if iPhase==1: return "liquid" # elif iPhase==2: return "supercritical" # elif iPhase==3: return "supercritical_gas" # elif iPhase==4: return "supercritical_liquid" # elif iPhase==5: return "critical_point" # elif iPhase==6: return "gas" # elif iPhase==7: return "twophase" # elif iPhase==8: return "unknown" # elif iPhase==9: return "not_imposed" # else: raise ValueError("Couldn't find phase.") def getPhaseNum(sPhase): return get_phase_index(sPhase) # for i in range(11): # if getPhaseString(i)==sPhase: # return i def getOutKey(out): return "".join(out) def getOutLabel(out): return ",".join(out) def getTimeKey(inp,out): return "_".join([inp,getOutKey(out)]) def getVectorKey(inp,out): return getTimeKey(inp, out)+"_V" def getCriticalProps(propsfluid): backend,_,_ = splitFluid(propsfluid) if backend!="INCOMP": p_crit_m = PropsSI('pcrit',"T",0,"D",0,propsfluid) *0.995 T_crit_m = PropsSI('Tcrit',"T",0,"D",0,propsfluid) *1.005 d_crit_m = PropsSI('rhocrit',"T",0,"D",0,propsfluid)*0.995 h_crit_m = PropsSI('H',"T",T_crit_m,"D",d_crit_m,propsfluid) s_crit_m = PropsSI('H',"T",T_crit_m,"D",d_crit_m,propsfluid) else: p_crit_m = None T_crit_m = None d_crit_m = None h_crit_m = None s_crit_m = None return dict(P=p_crit_m,T=T_crit_m,D=d_crit_m,H=h_crit_m,S=s_crit_m) def getPTRanges(propsfluid): backend,_,_ = splitFluid(propsfluid) # Setting the limits for enthalpy and pressure if backend == "REFPROP": T_min = PropsSI('Ttriple',"T",0,"D",0,propsfluid) T_max = PropsSI( 'Tmax',"T",0,"D",0,propsfluid) p_min = PropsSI( 'P',"T",T_min,"Q",0,propsfluid) p_max = PropsSI( 'pmax',"T",0,"D",0,propsfluid) elif backend == "INCOMP": T_min = PropsSI('Tmin',"T",0,"D",0,propsfluid) T_max = PropsSI('Tmax',"T",0,"D",0,propsfluid) p_min = 1.5*1e5 p_max = 200.0*1e5 else: p_min = PropsSI('ptriple',"T",0,"D",0,propsfluid) p_max = PropsSI( 'pmax',"T",0,"D",0,propsfluid) T_min = PropsSI('Ttriple',"T",0,"D",0,propsfluid) T_max = PropsSI( 'Tmax',"T",0,"D",0,propsfluid) p_range = np.logspace(np.log10(p_min),np.log10(p_max),xypoints) T_range = np.linspace(T_min,T_max,xypoints) return p_range,T_range #p_max = min(PropsSI('pcrit',"T",0,"D",0,fluid)*20, p_max) #T_max = min(PropsSI('Tcrit',"T",0,"D",0,fluid)* 3, T_max) def getLists(propsfluid): backend,_,_ = splitFluid(propsfluid) """Returns randomised lists of all properties within the ranges""" p,T = getPTRanges(propsfluid) p_min = np.min(p) p_max = np.max(p) T_min = np.min(T) T_max = np.max(T) if backend=="INCOMP": h_min = PropsSI('H','T',T_min,'P',p_min,propsfluid) h_max = PropsSI('H','T',T_max,'P',p_max,propsfluid) else: critProps = getCriticalProps(propsfluid) h_min = PropsSI('H','T',T_min,'Q',0,propsfluid) h_max = PropsSI('H','T',T_min,'Q',1,propsfluid) h_max = max(PropsSI('H','T',critProps["T"],'D',critProps["D"],propsfluid),h_max) h_max = (h_max-h_min)*2.0+h_min loop = True count = 0 while loop: count += 1 h_list = np.random.uniform( h_min, h_max,int(xypoints*2.0)) p_list = np.random.uniform(np.log10(p_min),np.log10(p_max),int(xypoints*2.0)) p_list = np.power(10,p_list) out = ["T","D","S"] res = PropsSI(out,"P",p_list,"H",h_list,propsfluid) T_list = res[:,0] d_list = res[:,1] s_list = res[:,2] mask = np.isfinite(T_list) & np.isfinite(d_list) & np.isfinite(s_list) if np.sum(mask) {1:d} - rerunning".format(xypoints,np.sum(mask))) loop = True else: loop = False p_list = p_list[mask][0:xypoints] h_list = h_list[mask][0:xypoints] T_list = T_list[mask][0:xypoints] d_list = d_list[mask][0:xypoints] s_list = s_list[mask][0:xypoints] return dict(P=p_list,T=T_list,D=d_list,H=h_list,S=s_list) maxTries = 4 if count > maxTries: loop = False raise ValueError("{0}: Could not fill the lists in {0} runs, aborting.".format(propsfluid,maxTries)) def getInpValues(inp,dataDict): in1 = inp[0] in2 = dataDict[in1] in3 = inp[1] in4 = dataDict[in3] return in1,in2,in3,in4 def getStateObj(propsfluid): backend,fld,conc = splitFluid(propsfluid) # fluidstr holds the full information and fluid is only the name # Initialise the state object if backend is not None: state = State(backend,fld) else: state = State(fld) #if backend=="INCOMP": # state.set_mass_fractions([0.0]) if conc is not None: state.set_mass_fractions([conc]) return state def getSpeedMeas(out,in1,in2,in3,in4,propsfluid,vector=False): pair, out1, _ = generate_update_pair(get_parameter_index(in1),in2[0],get_parameter_index(in3),in4[0]) if out1==in2[0]: swap=False else: swap=True if swap: input1 = in4 input2 = in2 else: input1 = in2 input2 = in4 state = getStateObj(propsfluid) outList = [get_parameter_index(i) for i in out] resLst = np.empty((repeat,)) timLst = np.empty((repeat,)) if vector: for j in range(repeat): timLst.fill(np.inf) lrange = range(len(input1)) resTmp = np.inf if len(outList)==1 and outList[0]==CoolProp.constants.iT_max: t1=clock() for l in lrange: for o in outList: resTmp = state.keyed_output(o) t2=clock() timLst[j] = (t2-t1)*1e6/float(len(input1)) else: # We have to update before doing other things t1=clock() for l in lrange: state.update(pair,input1[l],input2[l]) for o in outList: resTmp = state.keyed_output(o) t2=clock() timLst[j] = (t2-t1)*1e6/float(len(input1)) res = None tim = np.min(timLst) # Best of (repeat) return res,tim else: res = np.empty_like(input1) res.fill(np.inf) tim = np.empty_like(input1) tim.fill(np.inf) for i in range(len(input1)): resLst.fill(np.inf) timLst.fill(np.inf) for j in range(repeat): lrange = range(loops) resTmp = np.inf if len(outList)==1 and outList[0]==CoolProp.constants.iT_max: t1=clock() for _ in lrange: for o in outList: resTmp = state.keyed_output(o) t2=clock() timLst[j] = (t2-t1)*1e6/float(loops) resLst[j] = resTmp else: # We have to update before doing other things inV1 = input1[i] inV2 = input2[i]#*(1.0+(l/1000.0)*pow(-1,l)) for l in lrange ] t1=clock() for l in lrange: state.update(pair,inV1,inV2) for o in outList: resTmp = state.keyed_output(o) t2=clock() timLst[j] = (t2-t1)*1e6/float(loops) resLst[j] = resTmp if not np.all(resLst==resLst[0]): raise ValueError("Not all results were the same.") res[i] = resLst[0] tim[i] = np.min(timLst) # Best of three (repeat) return res,tim def checkDataSet(propsfluid,dataDict,fill=True,quiet=False): if not check: return backend,_,_ = splitFluid(propsfluid) if not quiet: print("\n\n-- {0:16s} --".format(propsfluid),end="") # Test for required inputs newInputs = False inLists = getLists(propsfluid) for inp in getInpList(backend): if not quiet: print("\n{0:2s}: ".format(inp),end="") for inVal in inp: if inVal not in dataDict: # A problem if not fill: raise ValueError("The input {0:1s} is missing or faulty, cannot continue.".format(inVal)) #dataDict[inVal] = inLists[inVal] dataDict.update(inLists) newInputs = True if not quiet: print("{0:s}*({1:d}),".format(inVal,len(dataDict[inVal])),end="") else: if not quiet: print("{0:s} ({1:d}),".format(inVal,len(dataDict[inVal])),end="") # All inputs are there in1,in2,in3,in4 = getInpValues(inp,dataDict) #in2 = in2[:3] #in4 = in4[:3] if in2.shape != in4.shape: raise ValueError("The stored data for {0:s} and {1:s} do not have the same shape.".format(in1,in3)) if in2.shape != inLists[in1].shape: raise ValueError("The stored data for {0:s} and its list do not have the same shape {1} vs {2}.".format(in1,in2.shape,inLists[in1].shape)) # Check for time data for out in getOutList(inp): key = getTimeKey(inp, out) okey = getOutKey(out) if key not in dataDict or newInputs or not np.all(np.isfinite(dataDict[key])): if not fill: raise ValueError("The time data for {0:s} is missing or faulty, cannot continue.".format(key)) res,tim = getSpeedMeas(out, in1, in2, in3, in4, propsfluid) dataDict[key] = tim dataDict[okey] = res # We calculated in, why not use it here... if not quiet: print("{0:s}*({1:d}),".format(key,len(dataDict[key])),end="") else: if not quiet: print("{0:s} ({1:d}),".format(key,len(dataDict[key])),end="") if dataDict[key].shape != in2.shape or not np.all(np.isfinite(dataDict[key])): raise ValueError("The stored time data for {0:s} does not have the same shape as the inputs.".format(key)) # Check for vectors for out in getOutList(inp): key = getVectorKey(inp, out) if key not in dataDict or not np.all(np.isfinite(dataDict[key])): if not fill: raise ValueError("The fluid data for {0:s} is missing or faulty, cannot continue.".format(key)) res,tim = getSpeedMeas(out, in1, in2, in3, in4, propsfluid, vector=True) dataDict[key] = tim if not quiet: print("{0:s}*({1:d}),".format(key,dataDict[key].size),end="") else: if not quiet: print("{0:s} ({1:d}),".format(key,dataDict[key].size),end="") if dataDict[key].size!=1 or not np.all(np.isfinite(dataDict[key])): raise ValueError("The vector data for {0:s} does not have the correct size {1}..".format(key,dataDict[key].size)) # inp = getInpList(backend)[0] # Explicit calls # # Check for properties # for out in getOutList(inp)[:-1]: # key = getOutKey(out) # if key not in dataDict or not np.all(np.isfinite(dataDict[key])): # if not fill: # raise ValueError("The fluid data for {0:s} is missing or faulty, cannot continue.".format(key)) # res = PropsSI(out,in1,in2,in3,in4,propsfluid) # dataDict[key] = res # if not quiet: print("{0:s}*({1:d}),".format(key,len(dataDict[key])),end="") # else: # if not quiet: print("{0:s} ({1:d}),".format(key,len(dataDict[key])),end="") # if dataDict[key].shape != in2.shape or not np.all(np.isfinite(dataDict[key])): # raise ValueError("The stored data for {0:s} does not have the same shape as the inputs {1} vs {2}..".format(key,dataDict[key].shape,in2.shape)) # # Check for phase # for out in ["Phase"]: # if backend!="HEOS": # dataDict[key] = np.zeros_like(a, dtype, order, subok) # key = getOutKey(out) # if key not in dataDict or newInputs or not np.all(np.isfinite(dataDict[key])): # if not fill: # raise ValueError("The phase data for {0:s} is missing or faulty, cannot continue.".format(key)) # res = np.empty_like(in2) # res.fill(np.inf) # for i in range(len(in2)): # res[i] = PropsSI(out,in1,in2[i],in3,in4[i],propsfluid) # dataDict[key] = res # if not quiet: print("{0:s}*({1:d}),".format(key,len(dataDict[key])),end="") # else: # if not quiet: print("{0:s} ({1:d}),".format(key,len(dataDict[key])),end="") # if dataDict[key].shape != in2.shape or not np.all(np.isfinite(dataDict[key])): # raise ValueError("The stored data for {0:s} does not have the same shape as the inputs {1} vs {2}..".format(key,dataDict[key].shape,in2.shape)) # # # Now we use the vector data # key = getVectorKey(inp, out) # if key not in dataDict or not np.all(np.isfinite(dataDict[key])): # if not fill: # raise ValueError("The vector data for {0:s} is missing or faulty, cannot continue.".format(key)) # dataDict[key] = np.empty_like(in2) # dataDict[key].fill(np.inf) # res = [] # for _ in range(repeat): # t1=clock() # PropsSI(out,in1,in2,in3,in4,propsfluid) # t2=clock() # res.append((t2-t1)/float(len(in2))) # dataDict[key] = np.min(res) # if not quiet: print("{0:s}*({1}),".format(key,dataDict[key]),end="") # else: # if not quiet: print("{0:s} ({1}),".format(key,dataDict[key]),end="") # try: # float(dataDict[key]) # except: # raise ValueError("The stored vector data for {0:s} cannot be casted to float.".format(key)) # if not quiet: print("") # All data is loaded and checked, we can calculate more now def getEntryCount(dicts,backend,fld): return len(fluidData[fld][backend].keys()) def getUKey(fld,bck,inp,out): return "-".join([fld,bck,inp,"".join(out)]) def getData(fld,backend,inp,out,fluidData): inputs1 = [] inputs2 = [] values = [] times = [] i1key = inp[0] i2key = inp[1] vkey = getOutKey(out) tkey = getTimeKey(inp,out) for dkey in fluidData[fld][backend]: cData = fluidData[fld][backend][dkey] inputs1.append(cData[i1key]) inputs2.append(cData[i2key]) values.append( cData[vkey] ) times.append( cData[tkey] ) ret = {} if len(inputs1)>0: ret[i1key] = np.concatenate(inputs1) ret[i2key] = np.concatenate(inputs2) ret[vkey] = np.concatenate(values ) ret[tkey] = np.concatenate(times ) return ret def getSingleData(fld,backend,key,fluidData): #print("Getting: "+fld+", "+backend+", "+key) values = [] for dkey in fluidData[fld][backend]: if key in fluidData[fld][backend][dkey]: if "P" in fluidData[fld][backend][dkey]: # TODO: Fix this, do we need the mask? #mask = fluidData[fld][backend][dkey]["P"]>0.3e5 mask = fluidData[fld][backend][dkey]["P"]>0.3e5 try: values.append( fluidData[fld][backend][dkey][key][mask] ) except Exception as e: values.append( fluidData[fld][backend][dkey][key] ) print(e) pass else: values.append( fluidData[fld][backend][dkey][key] ) if len(values)>0: if np.size(values[0])>1: return np.concatenate(values) else: return np.array(values) return None def fillDict(fld,backend,fluidData,curDict,curKeys): if curDict is None: curDict = {} for key in curKeys: vals = getSingleData(fld, backend, key, fluidData) if vals is not None: curDict[key] = vals return curDict ################################################ fluidData = {} for fluidstr in fluids: _,fld,_ = splitFluid(fluidstr) if fld not in fluidData: fluidData[fld] = {} for backend in backends: if backend not in fluidData[fld]: # Try to add it propsfluid = "::".join([backend,fluidstr]) try: PropsSI('Tmax',"T",0,"D",0,propsfluid) fluidData[fld][backend] = loadNpzData(backend,fld) except: pass if backend in fluidData[fld]: for k in fluidData[fld][backend]: checkDataSet(propsfluid, fluidData[fld][backend][k], fill=False, quiet=True) else: # Already added backend for fluid pass lastSave = 0 while runs999: raise ValueError("Too many dicts: {0}>999".format(newKey)) k = str(newKey).zfill(3) dicts[k] = {} dicts[k]['name']=k try: checkDataSet(propsfluid, dicts[k], fill=True, quiet=False) except Exception as e: print("There was an error, dataset {0} from {1}.might be faulty:\n{2}".format(k,propsfluid,str(e))) pass todel = [] #for k in dicts: try: checkDataSet(propsfluid, dicts[k], fill=False, quiet=True) except Exception as e: print("There was an error, removing dataset {0} from {1}.\n{2}".format(k,propsfluid,str(e))) todel.append(k) for t in todel: del dicts[t] keys = list(dicts.keys()) keys.sort() # Updated all dicts for this backend, saving data if runs>=maxruns or (lastSave+4)", rotation=0, ha='left', va='bottom', fontsize='medium') handles = [] for h in (rects1[0], rects2[1], rects3[2]): handles.append(copy.copy(h)) handles[-1].set_facecolor('white') handles.append(copy.copy(h)) handles[-1].set_hatch('') labels = (r'$p,T$-fit', r'constant', r'CoolProp', r'explicit, $f(p|\rho,T)$', r'REFPROP', r'implicit, $f(h,p)$') if bp is not None: bp.drawLegend(ax=ax1, loc='lower center', bbox_to_anchor=(0.5, 1.05), ncol=3, handles=handles, labels=labels) else: ax1.legend(handles,labels, loc='lower center', bbox_to_anchor=(0.5, 1.05), ncol=3) ax1.set_ylabel(r'Time per explicit call (us)') ax2.set_ylabel(r'Time per implicit call (us)') fg.savefig(path.join(getFigureFolder(),"TimeComp-"+fld.lower()+".pdf")) if bp is not None: ax1.set_ylabel(r'Time per explicit call (\si{\us})') ax2.set_ylabel(r'Time per implicit call (\si{\us})') mpl.rcParams['text.usetex'] = True fg.savefig(path.join(getFigureFolder(),"TimeComp-"+fld.lower()+"-tex.pdf")) mpl.rcParams['text.usetex'] = False # Fix the wrong baseline for tick in ax1.get_xaxis().get_major_ticks(): tick.set_pad(2*tick.get_pad()) tick.label1 = tick._get_text1() for lab in ax1.xaxis.get_ticklabels(): lab.set_verticalalignment("baseline") #lab.set_pad(1.5*lab.get_pad()) #ax1.set_xticklabels(labels) # #for tick in ax1.xaxis.get_major_ticks(): # tick.label1.set_horizontalalignment('center') bp.savepgf(path.join(getFigureFolder(),"TimeComp-"+fld.lower()+".pgf"),fg,repList) plt.close() # for fld in fluids: # try: # if bp is not None: # bp.figure = None # fg = bp.getFigure() # ccycle = bp.getColorCycle(length=3) # else: # fg = plt.figure() # ccycle = cycle(["b","g","r"]) # # #fg = plt.figure() # ax1 = fg.add_subplot(111) # ax2 = ax1.twinx() # # if "INCOMP" in fluidData[fld]: # el = 3 # hatch = ["","//","x"] # else: # el = 2 # hatch = ["//","x"] # # #one standard deviation above and below the mean of the data # width = 0.25 # the width of the bars # step = 1 # offset = -step-0.5*el*width # # lab = [] # rects1 = [] # rects2 = [] # # curList = [] # if el==3: curList.append(getSingleData(fld, "INCOMP" , "PT_Tmax", fluidData)*1e2) # curList.append(getSingleData(fld, "HEOS" , "DT_Tmax", fluidData)*1e2) # curList.append(getSingleData(fld, "REFPROP", "DT_Tmax", fluidData)*1e2) # # lab.extend(["Tmax"]) # # offset += step # curN = len(curList) # curInd = [ offset + i * width for i in range(curN)] # curCol = ccycle.next() # for i in range(len(hatch)): # rects1.extend(ax1.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k')) # autolabel(ax1,rects1[0:],np.mean(curList,axis=1),np.std(curList,axis=1),fac=1e-2) # # curList = [] # if el==3: curList.append(getSingleData(fld, "INCOMP" , "PT_H", fluidData)) # curList.append(getSingleData(fld, "HEOS" , "DT_H", fluidData)) # curList.append(getSingleData(fld, "REFPROP", "DT_H", fluidData)) # lab.extend(["H"]) # offset += step # curN = len(curList) # curInd = [ offset + i * width for i in range(curN)] # curCol = ccycle.next() # for i in range(len(hatch)): # rects1.extend(ax1.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k')) # autolabel(ax1,rects1[el:],np.mean(curList,axis=1),np.std(curList,axis=1)) # # curList = [] # if el==3: curList.append(getSingleData(fld, "INCOMP" , "PT_D", fluidData)) # curList.append(getSingleData(fld, "HEOS" , "DT_P", fluidData)) # curList.append(getSingleData(fld, "REFPROP", "DT_P", fluidData)) # if el==3: lab.extend(["D/P"]) # else: lab.extend(["P"]) # offset += step # curN = len(curList) # curInd = [ offset + i * width for i in range(curN)] # #curCol = "g" # for i in range(len(hatch)): # rects1.extend(ax1.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k')) # autolabel(ax1,rects1[int(2*el):],np.mean(curList,axis=1),np.std(curList,axis=1)) # # curList = [] # if el==3: curList.append(getSingleData(fld, "INCOMP" , "PT_S", fluidData)) # curList.append(getSingleData(fld, "HEOS" , "DT_S", fluidData)) # curList.append(getSingleData(fld, "REFPROP", "DT_S", fluidData)) # lab.extend(["S"]) # offset += step # curN = len(curList) # curInd = [ offset + i * width for i in range(curN)] # #curCol = "g" # for i in range(len(hatch)): # rects1.extend(ax1.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k')) # autolabel(ax1,rects1[int(3*el):],np.mean(curList,axis=1),np.std(curList,axis=1)) # # curList = [] # if el==3: curList.append(getSingleData(fld, "INCOMP" , "PT_HDS", fluidData)) # curList.append(getSingleData(fld, "HEOS" , "DT_HPS", fluidData)) # curList.append(getSingleData(fld, "REFPROP", "DT_HPS", fluidData)) # lab.extend(["all"]) # offset += step # curN = len(curList) # curInd = [ offset + i * width for i in range(curN)] # #curCol = "g" # for i in range(len(hatch)): # rects1.extend(ax1.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k')) # autolabel(ax1,rects1[int(4*el):],np.mean(curList,axis=1),np.std(curList,axis=1)) # # curList = [] # if el==3: curList.append(getSingleData(fld, "INCOMP" , "HP_D", fluidData)) # curList.append(getSingleData(fld, "HEOS" , "HP_D", fluidData)) # curList.append(getSingleData(fld, "REFPROP", "HP_D", fluidData)) # lab.extend(["D"]) # offset += step # curN = len(curList) # curInd = [ offset + i * width for i in range(curN)] # curCol = ccycle.next() # for i in range(len(hatch)): # rects2.extend(ax2.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k')) # autolabel(ax2,rects2[0:],np.mean(curList,axis=1),np.std(curList,axis=1)) # # curList = [] # if el==3: curList.append(getSingleData(fld, "INCOMP" , "HP_T", fluidData)) # curList.append(getSingleData(fld, "HEOS" , "HP_T", fluidData)) # curList.append(getSingleData(fld, "REFPROP", "HP_T", fluidData)) # lab.extend(["T"]) # offset += step # curN = len(curList) # curInd = [ offset + i * width for i in range(curN)] # #curCol = "r" # for i in range(len(hatch)): # rects2.extend(ax2.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k')) # autolabel(ax2,rects2[int(1*el):],np.mean(curList,axis=1),np.std(curList,axis=1)) # # curList = [] # if el==3: curList.append(getSingleData(fld, "INCOMP" , "HP_S", fluidData)) # curList.append(getSingleData(fld, "HEOS" , "HP_S", fluidData)) # curList.append(getSingleData(fld, "REFPROP", "HP_S", fluidData)) # lab.extend(["S"]) # offset += step # curN = len(curList) # curInd = [ offset + i * width for i in range(curN)] # #curCol = "r" # for i in range(len(hatch)): # rects2.extend(ax2.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k')) # autolabel(ax2,rects2[int(2*el):],np.mean(curList,axis=1),np.std(curList,axis=1)) # # curList = [] # if el==3: curList.append(getSingleData(fld, "INCOMP" , "HP_DST", fluidData)) # curList.append(getSingleData(fld, "HEOS" , "HP_DST", fluidData)) # curList.append(getSingleData(fld, "REFPROP", "HP_DST", fluidData)) # lab.extend(["all"]) # offset += step # curN = len(curList) # curInd = [ offset + i * width for i in range(curN)] # #curCol = "r" # for i in range(len(hatch)): # rects2.extend(ax2.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k')) # autolabel(ax2,rects2[int(3*el):],np.mean(curList,axis=1),np.std(curList,axis=1)) # # # ids = np.arange(len(lab)) # # add some text for labels, title and axes ticks # #if backend=="INCOMP": ax1.set_ylabel(r'Time per $f(p,T)$ call (\si{\us})') # if el==3: ax1.set_ylabel(r'Time per \texttt{Tmax} call (\SI{0.01}{\us}) and'+"\n"+r'per $f(p,T)$ and $f(\rho,T)$ call (\si{\us})') # else: ax1.set_ylabel(r'Time per \texttt{Tmax} call (\SI{0.01}{\us})'+"\n"+r'and per $f(\rho,T)$ call (\si{\us})') # ax2.set_ylabel(r'Time per $f(h,p)$ call (\si{\us})') # # ax1.set_xticks(ids) # ax1.set_xticklabels([r"\texttt{"+i+r"}" for i in lab], rotation=0) # ax1.set_xlim([ids.min()-2.5*width,ids.max()+2.5*width]) # # ax1.spines['top'].set_visible(False) # ax2.spines['top'].set_visible(False) # ax1.xaxis.set_ticks_position('bottom') # ax2.xaxis.set_ticks_position('bottom') # # handles = [] # if el==3: # for h in (rects1[0], rects1[4], rects2[2]): # handles.append(copy.copy(h)) # handles[-1].set_facecolor('white') # handles.append(copy.copy(h)) # handles[-1].set_hatch('') # labels = (r'$p,T$-fit', r'\texttt{Tmax}', r'CoolProp', r'explicit, $f(p|\rho,T)$', r'REFPROP', r'implicit, $f(h,p)$') # else: # for h in (rects1[0], rects1[2], rects2[1]): # handles.append(copy.copy(h)) # handles[-1].set_facecolor('white') # handles.append(copy.copy(h)) # handles[-1].set_hatch('') # labels = (r'', r'\texttt{Tmax}', r'CoolProp', r'explicit, $f(\rho,T)$', r'REFPROP', r'implicit, $f(h,p)$') # handles[0] = mpatches.Patch(visible=False) # # if bp is not None: # bp.drawLegend(ax=ax1, # loc='upper center', # bbox_to_anchor=(0.5, 1.4), # ncol=3, # handles=handles, # labels=labels) # else: # ax1.legend(handles,labels, # loc='upper center', # bbox_to_anchor=(0.5, 1.4), # ncol=3) # # fg.savefig(path.join(getFigureFolder(),"TimeComp-"+fld+".pdf")) # if bp is not None: bp.savepgf(path.join(getFigureFolder(),"TimeComp-"+fld+".pgf"),fg,repList) # plt.close() # # except Exception as e: # print(e) # pass ############################################################# # The second figure compares the backend for the full calculation ############################################################# backendExp = [] backendImp = [] fluidLabel = [] hatchLst = ["///" ,"\\\\\\"] backendsLst = ["HEOS","REFPROP"] for fluidstr in fluids[:-1]: _,fld,_ = splitFluid(fluidstr) #if fld=="CO2": fluidLabel.append("\ce{CO2}") #else: fluidLabel.append(fld) for backend in backendsLst: if backend not in fluidData[fld]: backendExp.append([0]) backendImp.append([0]) continue # Backend exists in fluid data try: inp = getInpList(backend) outExp = getTimeKey(inp[0], getOutList(inp[0])[-1]) outImp = getTimeKey(inp[1], getOutList(inp[1])[-1]) backendExp.append(getSingleData(fld, backend, outExp, fluidData)) backendImp.append(getSingleData(fld, backend, outImp, fluidData)) except Exception as e: backendExp.append([0]) backendImp.append([0]) print(e) pass # Data is prepared, we can plot now. if bp is not None: bp.figure = None fg1 = bp.getFigure() bp2 = BasePlotter() fg2 = bp2.getFigure() ccycle = bp.getColorCycle(length=3) else: fg1 = plt.figure() fg2 = plt.figure() ccycle = cycle(["b","g","r"]) fg1.set_size_inches( (fg1.get_size_inches()[0]*1, fg1.get_size_inches()[1]*0.75) ) fg2.set_size_inches( (fg2.get_size_inches()[0]*1, fg2.get_size_inches()[1]*0.75) ) ccycle.next() # No incomp # #ax1 = fg.add_subplot(111) #ax2 = ax1.twinx() ax1 = fg1.add_subplot(111) ax2 = fg2.add_subplot(111) #entries = int(len(backendExp)/len(fluidLabel)) #one standard deviation above and below the mean of the data rects1 = [] rects2 = [] col1 = ccycle.next() col2 = ccycle.next() numFluids = len(fluidLabel) numBackends = len(backendsLst) step = 1 width = step/(numBackends+1) # the width of the bars offset = -0.5*numBackends*width for f in range(numFluids): ids = np.empty((numBackends,)) for b in range(numBackends): i = f*numBackends+b ids[b] = offset + f*step + b*width rects1.extend(ax1.bar(ids[b], np.mean(backendExp[i]), width, color=col1, hatch=hatchLst[b], rasterized=False))#, yerr=np.std(curList[i]), ecolor='k')) rects2.extend(ax2.bar(ids[b], np.mean(backendImp[i]), width, color=col2, hatch=hatchLst[b], rasterized=False))#, yerr=np.std(curList[i]), ecolor='k')) y_max = np.max(np.concatenate((np.ravel(ax1.get_ylim()),np.ravel(ax2.get_ylim())/10.0))) ax1.set_ylim([0,y_max]) ax2.set_ylim([0,y_max*10.0]) for ax in [ax1,ax2]: ax.set_xticks(range(numFluids)) ax.set_xticklabels(fluidLabel, rotation=25) ax.set_xlim([0.0-0.5*step,numFluids-1+0.5*step]) ax.spines['top'].set_visible(False) ax.xaxis.set_ticks_position('bottom') if ax==ax1: rects = rects1 elif ax==ax2: rects = rects2 handles = (rects[0], rects[1]) labels = (r'CoolProp', r'REFPROP') anchor = (0.5, 1.2) if bp is not None: bp.drawLegend(ax=ax, loc='upper center', bbox_to_anchor=anchor, ncol=numBackends, handles=handles, labels=labels) else: ax.legend(handles,labels, loc='upper center', bbox_to_anchor=anchor, ncol=numBackends) ax1.set_ylabel(r'Time per $f(\rho,T)$ call (us)') ax2.set_ylabel(r'Time per $f(h,p)$ call (us)') fg1.savefig(path.join(getFigureFolder(),"TimeComp-backends-exp.pdf")) fg2.savefig(path.join(getFigureFolder(),"TimeComp-backends-imp.pdf")) if bp is not None: ax1.set_ylabel(r'Time per $f(\rho,T)$ call (\si{\us})') ax2.set_ylabel(r'Time per $f(h,p)$ call (\si{\us})') mpl.rcParams['text.usetex'] = True fg1.savefig(path.join(getFigureFolder(),"TimeComp-backends-exp-tex.pdf")) fg2.savefig(path.join(getFigureFolder(),"TimeComp-backends-imp-tex.pdf")) mpl.rcParams['text.usetex'] = False bp.savepgf(path.join(getFigureFolder(),"TimeComp-backends-exp.pgf"),fg1,repList) bp.savepgf(path.join(getFigureFolder(),"TimeComp-backends-imp.pgf"),fg2,repList) plt.close('all') ############################################################# # The third figure is a heat map of the execution times in # log p h diagram ############################################################# for fluidstr in fluids: try: _,fld,_ = splitFluid(fluidstr) for backend in fluidData[fld]: propsfluid = "::".join([backend,fluidstr]) if backend!="INCOMP": TP = {} points = int(xypoints/2) T_range_TP = np.linspace(PropsSI('Ttriple',"T",0,"D",0,propsfluid),PropsSI( 'Tcrit',"T",0,"D",0,propsfluid),points) T_TP = np.append(T_range_TP,T_range_TP[::-1]) Q_TP = np.zeros_like(T_TP) Q_TP[points:] = 1 points *= 2 out = ["D","H","P","S"] res = PropsSI(out,"T",T_TP,"Q",Q_TP,propsfluid) D_TP = res[:,0] H_TP = res[:,1] P_TP = res[:,2] S_TP = res[:,3] mask = np.isfinite(D_TP) if np.sum(mask)