#!/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 def range_brace(x_min, x_max, mid=0.5, beta1=50.0, beta2=100.0, height=1, initial_divisions=11, resolution_factor=1.5): """ http://stackoverflow.com/questions/1289681/drawing-braces-with-pyx x,y = range_brace(0, 100) ax.plot(x, y,'-') ax.plot(y, x,'-') """ # determine x0 adaptively values using second derivitive # could be replaced with less snazzy: # x0 = NP.arange(0, 0.5, .001) x0 = np.array(()) tmpx = np.linspace(0, 0.5, initial_divisions) tmp = beta1**2 * (np.exp(beta1*tmpx)) * (1-np.exp(beta1*tmpx)) / np.power((1+np.exp(beta1*tmpx)),3) tmp += beta2**2 * (np.exp(beta2*(tmpx-0.5))) * (1-np.exp(beta2*(tmpx-0.5))) / np.power((1+np.exp(beta2*(tmpx-0.5))),3) for i in range(0, len(tmpx)-1): t = int(np.ceil(resolution_factor*max(np.abs(tmp[i:i+2]))/float(initial_divisions))) x0 = np.append(x0, np.linspace(tmpx[i],tmpx[i+1],t)) x0 = np.sort(np.unique(x0)) # sort and remove dups # half brace using sum of two logistic functions y0 = mid*2*((1/(1.+np.exp(-1*beta1*x0)))-0.5) y0 += (1-mid)*2*(1/(1.+np.exp(-1*beta2*(x0-0.5)))) # concat and scale x x = np.concatenate((x0, 1-x0[::-1])) * float((x_max-x_min)) + x_min y = np.concatenate((y0, y0[::-1])) * float(height) return (x,y) #try: #from jopy.dataPlotters import BasePlotter #bp = BasePlotter() #except: #bp = None bp = None # The basic settings for he plots xypoints = 1000 loops = 1 repeat = 1 runs = 0 maxruns = 5 plot = True calc = True check = True folder = "dataTTSE" figures = "figuresTTSE" #np.random.seed(1984) fluids = ["CO2","Pentane","R134a","Water","Air","LiBr-0%"] fluids = ["CO2","Pentane","R134a","Water"] fluids = ["Air"] #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"] backends = ["HEOS","REFPROP"] backends = ["HEOS"] #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 T_min = PropsSI('Tmin',"T",0,"D",0,propsfluid)+1 T_max = PropsSI('Tmax',"T",0,"D",0,propsfluid)-1 if backend == "REFPROP": T_min = max(T_min,PropsSI('Ttriple',"T",0,"D",0,propsfluid))+1 p_min = PropsSI( 'P',"T",T_min,"Q",0,propsfluid)+1 p_max = PropsSI( 'pmax',"T",0,"D",0,propsfluid)-1 elif backend == "INCOMP": p_min = 1.5*1e5 p_max = 200.0*1e5 else: T_min = max(T_min,PropsSI('Ttriple',"T",0,"D",0,propsfluid))+1 p_min = PropsSI('ptriple',"T",0,"D",0,propsfluid) p_min = max(p_min,PropsSI('pmin',"T",0,"D",0,propsfluid))+1 p_max = PropsSI( 'pmax',"T",0,"D",0,propsfluid)-1 # One more check to debug things: #p_min = max(p_min,0.01e5) #T_min = max(T_min,200) #p_max = min(p_max,200e5) #T_max = min(T_max,1750) 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: try: state.set_mass_fractions([conc]) except: pass 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.0e5 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 = max(int(xypoints/2),250) T_range_TP = np.linspace(PropsSI('Ttriple',"T",0,"D",0,propsfluid)+1,PropsSI( 'Tcrit',"T",0,"D",0,propsfluid)-0.1,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)