github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-01-12 23:48:22 -05:00
Code Issues Packages Projects Releases Wiki Activity
Files
paramstringlen
CoolProp/dev/TTSE/TimeComp.py
Ian Bell 0e7437f95b More pep8 whitespace; see #1612
2019-01-12 20:48:56 -07:00

1665 lines
68 KiB
Python
Raw Permalink Blame History

#!/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 derivative
# 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) < xypoints:
if False:
print(h_list); print(p_list); print(T_list); print(d_list); print(s_list)
print("There were not enough valid entries in your result vector: {0:d} > {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 runs < maxruns and calc:
check = True # force checking for new records
runs += 1
# Now we have the data structure with the precalculated data
for fluidstr in fluids:
_, fld, _ = splitFluid(fluidstr)
for backend in fluidData[fld]:
propsfluid = "::".join([backend, fluidstr])
dicts = fluidData[fld][backend]
keys = list(dicts.keys())
keys.sort()
while len(keys) < runs:
if len(keys) < 1: newKey = 0
else: newKey = int(keys[-1]) + 1
if newKey > 999:
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) < runs:
print("\n\nDone processing fluids, saving data: ")
for fld in fluidData:
for backend in fluidData[fld]:
saveNpzData(backend, fld, fluidData[fld][backend], start=lastSave, stop=runs)
print("{0} ({1})".format(backend + "::" + fld, len(fluidData[fld][backend].keys()[lastSave:runs])), end=", ")
print("")
lastSave = runs
if not plot: sys.exit(0)
# exclLst = [["Tmax"],["H"],["D"],["S"],["H","D","S"],["D"],["S"],["T"],["D","S","T"]]
#
# Start with a temporary dictionary that holds all the data we need
# for fld in fluidData:
# cstData = {} # Data from calls to constants (overhead)
# expData = {} # Data from explicit EOS calls
# impData = {} # Data from EOS calls that require iterations
# for backend in fluidData[fld]:
#
# curKeys = []
# for inp in getInpList(backend):
# for out in getOutList(inp)[:1]:
# curKeys.append(getTimeKey( inp, out))
# curKeys.append(getVectorKey(inp, out))
# curDict = {}
# fillDict(fld,backend,fluidData,curDict,curKeys)
# cstData[backend] = curDict
#
# curKeys = []
# for inp in getInpList(backend)[:1]:
# for out in getOutList(inp)[1:]:
# curKeys.append(getTimeKey( inp, out))
# curKeys.append(getVectorKey(inp, out))
# curDict = {}
# fillDict(fld,backend,fluidData,curDict,curKeys)
# expData[backend] = curDict
#
# curKeys = []
# for inp in getInpList(backend)[1:]:
# for out in getOutList(inp)[1:]:
# curKeys.append(getTimeKey( inp, out))
# curKeys.append(getVectorKey(inp, out))
# curDict = {}
# fillDict(fld,backend,fluidData,curDict,curKeys)
# impData[backend] = curDict
# curDict[backend] = {}
# for key in :
# vals = getSingleData(fld, backend, key, fluidData)
# if vals is not None: curDict[backend][key] = vals
# if curDict
#
# cstData[backend] = {}
# for out in getOutList(inp[0])[0]:
# res = getData(fld,backend,inp[0],out,fluidData)
# cstData[backend].update(res)
# for out in getOutList(inp[1])[0]:
# res = getData(fld,backend,inp[1],out,fluidData)
# cstData[backend].update
#
# expData[backend] = {}
# for out in getOutList(inp[0])[1:]:
# res = getData(fld,backend,inp[0],out,fluidData)
# expData[backend].update(res)
#
# impData[backend] = {}
# for out in getOutList(inp[1])[1:]:
# res = getData(fld,backend,inp[1],out,fluidData)
# impData[backend].update(res)
#############################################################
# All data is available in the dicts now.
#############################################################
# The first thuing to do is to print some statistical
# measures to give you an idea about the data.
# try:
# #dataOHCP = [cstData["HEOS"]["DT_Tmax"] , cstData["HEOS"]["HP_Tmax"] ]
# #dataOHRP = [cstData["REFPROP"]["DT_Tmax"], cstData["REFPROP"]["HP_Tmax"]]
# print("\n{0} - {1} points ".format(fld,np.size(cstData["HEOS"]["DT_Tmax"])))
# print("Overhead CoolProp: {0:5.3f} us".format(np.mean(cstData["HEOS"]["DT_Tmax"])))#, np.mean(cstData["HEOS"]["HP_Tmax"]))
# print("Overhead REFPROP : {0:5.3f} us".format(np.mean(cstData["REFPROP"]["DT_Tmax"])))#,np.mean(cstData["REFPROP"]["HP_Tmax"]))
# print("Mean EOS in CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.mean(expData["HEOS"]["DT_HPS"]),np.mean(impData["HEOS"]["HP_DST"])))
# print("Std. dev. in CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.std(expData["HEOS"]["DT_HPS"]) ,np.std(impData["HEOS"]["HP_DST"])))
# print("Minimum in CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.min(expData["HEOS"]["DT_HPS"]) ,np.min(impData["HEOS"]["HP_DST"])))
# print("Maximum in CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.max(expData["HEOS"]["DT_HPS"]) ,np.max(impData["HEOS"]["HP_DST"])))
# print("Mean EOS in REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.mean(expData["REFPROP"]["DT_HPS"]),np.mean(impData["REFPROP"]["HP_DST"])))
# print("Std. dev. in REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.std(expData["REFPROP"]["DT_HPS"]) ,np.std(impData["REFPROP"]["HP_DST"])))
# print("Minimum in REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.min(expData["REFPROP"]["DT_HPS"]) ,np.min(impData["REFPROP"]["HP_DST"])))
# print("Maximum in REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.max(expData["REFPROP"]["DT_HPS"]) ,np.max(impData["REFPROP"]["HP_DST"])))
# print("")
#
# print("\n{0} - {1} points ".format(fld,np.size(cstData["HEOS"]["DT_Tmax_V"])))
# print("Overhead CoolProp: {0:5.3f} us".format(np.mean(cstData["HEOS"]["DT_Tmax_V"])))#, np.mean(cstData["HEOS"]["HP_Tmax"]))
# print("Overhead REFPROP : {0:5.3f} us".format(np.mean(cstData["REFPROP"]["DT_Tmax_V"])))#,np.mean(cstData["REFPROP"]["HP_Tmax"]))
# print("Mean EOS in CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.mean(expData["HEOS"]["DT_HPS_V"]),np.mean(impData["HEOS"]["HP_DST_V"])))
# print("Std. dev. in CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.std(expData["HEOS"]["DT_HPS_V"]) ,np.std(impData["HEOS"]["HP_DST_V"])))
# print("Minimum in CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.min(expData["HEOS"]["DT_HPS_V"]) ,np.min(impData["HEOS"]["HP_DST_V"])))
# print("Maximum in CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.max(expData["HEOS"]["DT_HPS_V"]) ,np.max(impData["HEOS"]["HP_DST_V"])))
# print("Mean EOS in REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.mean(expData["REFPROP"]["DT_HPS_V"]),np.mean(impData["REFPROP"]["HP_DST_V"])))
# print("Std. dev. in REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.std(expData["REFPROP"]["DT_HPS_V"]) ,np.std(impData["REFPROP"]["HP_DST_V"])))
# print("Minimum in REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.min(expData["REFPROP"]["DT_HPS_V"]) ,np.min(impData["REFPROP"]["HP_DST_V"])))
# print("Maximum in REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.max(expData["REFPROP"]["DT_HPS_V"]) ,np.max(impData["REFPROP"]["HP_DST_V"])))
# print("")
#
# except Exception as e:
# print(str(e.message))
# pass
#
# try:
# fld = 'Water'
# cstData = {} # Data from calls to constants (overhead)
# expData = {} # Data from explicit EOS calls
# impData = {} # Data from EOS calls that require iterations
# for backend in fluidData[fld]:
# curKeys = []
# for inp in getInpList(backend):
# for out in getOutList(inp)[:1]:
# curKeys.append(getTimeKey( inp, out))
# curKeys.append(getVectorKey(inp, out))
# curDict = {}
# fillDict(fld,backend,fluidData,curDict,curKeys)
# cstData[backend] = curDict
#
# curKeys = []
# for inp in getInpList(backend)[:1]:
# for out in getOutList(inp)[1:]:
# curKeys.append(getTimeKey( inp, out))
# curKeys.append(getVectorKey(inp, out))
# curDict = {}
# fillDict(fld,backend,fluidData,curDict,curKeys)
# expData[backend] = curDict
#
# curKeys = []
# for inp in getInpList(backend)[1:]:
# for out in getOutList(inp)[1:]:
# curKeys.append(getTimeKey( inp, out))
# curKeys.append(getVectorKey(inp, out))
# curDict = {}
# fillDict(fld,backend,fluidData,curDict,curKeys)
# impData[backend] = curDict
#
#
# print("Done")
def autolabel(ax, rects, means, stds, lens=0, fac=1):
return
# attach some text labels
yerr = (stds * 100.0) / means
ypos = np.max(means) # + np.max(stds)
for i in range(len(rects)):
xpos = rects[i].get_x() + rects[i].get_width() / 2.
#ax.text(xpos, 1.05*ypos, '{0:s}{1:4.2f}{2:s}{3:3.1f}{4:s}'.format(r'',means[i]*fac,r'us +- ',yerr[i],r'%'), rotation=90, ha='center', va='bottom', fontsize='smaller')
ax.text(xpos, 1.05 * ypos, '{0:s}{1:4.2f}{2:s}{3:3.1f}{4:s}'.format(r'\SI{', means[i] * fac, r'}{\us} (\SI{\pm', yerr[i], r'}{\percent})'), rotation=90, ha='center', va='bottom', fontsize='smaller')
#ax.text(xpos, 0.25*ypos, str(lens[i]), rotation=90, ha='center', va='bottom', fontsize='smaller')
def axislabel(txt, ax, xPos, yPos=-1):
ax.text(xPos, yPos, txt, rotation=45, ha='center', va='top', fontsize='xx-small')
#############################################################
# All data is available in the dicts now.
#############################################################
# The first plot contains the time data, this is averaged and
# plotted as a bar graph with the standard deviation.
hatchLst = ["", "///", "\\\\\\"]
backendsLst = ["INCOMP", "HEOS", "REFPROP"]
for fluidstr in fluids[:-1]:
_, fld, _ = splitFluid(fluidstr)
outCst = []
labCst = []
outExp = []
labExp = []
outImp = []
labImp = []
DEBUG = True
for backend in backendsLst:
# Backend exists in fluid data?
try:
for i, inp in enumerate(getInpList(backend)):
for j, out in enumerate(getOutList(inp)):
if j == 0: # First output is Tmax
if backend not in fluidData[fld]:
outCst.append([0])
labCst.append("Dummy")
if DEBUG: print("Added a dummy for {0} and {1},{2},{3}".format(fld, backend, inp, out))
else:
outCst.append(getSingleData(fld, backend, getTimeKey(inp, out), fluidData))
labCst.append(getTimeKey(inp, out))
continue
elif i == 0: # First input is explicit
if backend not in fluidData[fld]:
outExp.append([0])
labExp.append("Dummy")
if DEBUG: print("Added a dummy for {0} and {1},{2},{3}".format(fld, backend, inp, out))
else:
outExp.append(getSingleData(fld, backend, getTimeKey(inp, out), fluidData))
labExp.append(getTimeKey(inp, out))
continue
elif i == 1:
if backend not in fluidData[fld]:
outImp.append([0])
labImp.append("Dummy")
if DEBUG: print("Added a dummy for {0} and {1},{2},{3}".format(fld, backend, inp, out))
else:
outImp.append(getSingleData(fld, backend, getTimeKey(inp, out), fluidData))
labImp.append(getTimeKey(inp, out))
continue
else:
raise ValueError("Wrong number.")
except Exception as e:
print(e)
sys.exit(1)
# Do the plotting
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()
rects1 = []
labels1 = []
rects2 = []
labels2 = []
rects3 = []
labels3 = []
col1 = ccycle.next()
col2 = ccycle.next()
col3 = ccycle.next()
numBackends = len(backendsLst)
step = 1
width = step / (numBackends + 1) # the width of the bars
offset = -0.5 * numBackends * width
entries = int((len(outCst) + len(outExp) + len(outImp)) / numBackends)
# for o in range(entries):
# ids = np.empty((numBackends,))
# for b in range(numBackends):
# i = o*numBackends+b
# ids[b] = offset + o*step + b*width
# j = i - 0
# if j < len(outCst):
# rects1.extend(ax1.bar(ids[b], np.mean(outCst[j]), width, color=col1, hatch=hatchLst[b]))#, yerr=np.std(curList[i]), ecolor='k'))
#
#
# j = i - 0
# if j < len(outCst):
# rects1.extend(ax1.bar(ids[b], np.mean(outCst[j]), width, color=col1, hatch=hatchLst[b]))#, yerr=np.std(curList[i]), ecolor='k'))
# else:
# j = i-len(outCst)
# if j < len(outExp):
# rects2.extend(ax1.bar(ids[b], np.mean(outExp[j]), width, color=col2, hatch=hatchLst[b]))#, yerr=np.std(curList[i]), ecolor='k'))
# else:
# j = i-len(outCst)-len(outExp)
# if j < len(outImp):
# rects3.extend(ax2.bar(ids[b], np.mean(outImp[j]), width, color=col3, hatch=hatchLst[b]))#, yerr=np.std(curList[i]), ecolor='k'))
# else:
# raise ValueError("Do not go here!")
DEBUG = True
entries = 2
for o in range(entries):
ids = np.empty((numBackends,))
for b in range(numBackends):
i = b * entries + o
try:
ids[b] = offset + o * step + b * width
rects1.extend(ax1.bar(ids[b], np.mean(outCst[i]), width, color=col1, hatch=hatchLst[b], rasterized=False)) # , yerr=np.std(curList[i]), ecolor='k'))
if DEBUG:
print("Plotting {0}: {1:7.1f} - {2} - {3}".format(fld, np.mean(outCst[i]), "cst.", b))
# print(ids[b],labCst[i])
#axislabel(labCst[i], ax1, ids[b]+0.5*width)
except:
pass
offset += entries * step
#entries = int(len(outExp)/numBackends)
entries = 4
for o in range(entries):
ids = np.empty((numBackends,))
for b in range(numBackends):
i = b * entries + o
try:
ids[b] = offset + o * step + b * width
rects2.extend(ax1.bar(ids[b], np.mean(outExp[i]), width, color=col2, hatch=hatchLst[b], rasterized=False)) # , yerr=np.std(curList[i]), ecolor='k'))
if DEBUG:
print("Plotting {0}: {1:7.1f} - {2} - {3}".format(fld, np.mean(outExp[i]), "exp.", b))
# print(ids[b],labExp[i])
#axislabel(labExp[i], ax1, ids[b]+0.5*width)
except:
pass
x_newaxis = np.max(ids) + 1.5 * width
plt.axvline(x_newaxis, color='k', linestyle='dashed')
offset += entries * step
entries = 4
for o in range(entries):
ids = np.empty((numBackends,))
for b in range(numBackends):
i = b * entries + o
try:
ids[b] = offset + o * step + b * width
rects3.extend(ax2.bar(ids[b], np.mean(outImp[i]), width, color=col3, hatch=hatchLst[b], rasterized=False)) # , yerr=np.std(curList[i]), ecolor='k'))
if DEBUG:
print("Plotting {0}: {1:7.1f} - {2} - {3}".format(fld, np.mean(outImp[i]), "imp.", b))
# print(ids[b],labImp[i])
#axislabel(labImp[i], ax1, ids[b]+0.5*width)
except:
pass
# 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')
labels = [r"ex.", r"im.", r"$h$", r"$\rho|p$", r"$s$", r"all", r"$\rho$", r"$s$", r"$T$", r"all"]
ax1.set_xticks(range(len(labels)))
ax1.set_xticklabels(labels)
# ax1.yaxis.get_label().set_verticalalignment("baseline")
x_min = rects1[0].get_x()
dx = rects1[0].get_width()
x_max = rects3[-1].get_x()
x_min = x_min - 1 * dx
x_max = x_max + 2 * dx
ax1.set_xlim([x_min, x_max])
y_min = 0
y_max_c = np.nanmax([a.get_height() for a in rects1])
y_max_e = np.nanmax([a.get_height() for a in rects2])
y_max_i = np.nanmax([a.get_height() for a in rects3])
y_max = np.max([y_max_c, y_max_e, y_max_i / 10.0])
y_max = np.ceil(1.3 * y_max / 10.0) * 10.0
ax1.set_ylim([y_min, y_max])
ax2.set_ylim([y_min, y_max * 10.0])
ratio = 10.0 / 4.0 * y_max / 250.0 # height of 10 for 4 points if y_max==250
x_min = rects1[0].get_x()
x_max = rects1[-1].get_x() + dx
x, y = range_brace(x_min, x_max)
dy = np.ceil(y_max_c / 10.0) * 10.0
y = dy + y * ratio * (x[-1] - x[0])
ax1.plot(x, y, ls='-', color='k')
ax1.text(np.mean(x), np.max(y), "const.", rotation=0, ha='center', va='bottom', fontsize='medium')
x_min = rects2[0].get_x()
x_max = rects2[-1].get_x() + dx
x, y = range_brace(x_min, x_max)
dy = np.ceil(y_max_e / 10.0) * 10.0
y = dy + y * ratio * (x[-1] - x[0])
ax1.plot(x, y, ls='-', color='k')
ax1.text(np.mean(x), np.max(y), "explicit", rotation=0, ha='center', va='bottom', fontsize='medium')
x_min = rects3[0].get_x()
x_max = rects3[-1].get_x() + dx
x, y = range_brace(x_min, x_max)
dy = np.ceil(y_max_i / 100.0) * 10
y = dy + y * ratio * (x[-1] - x[0])
ax1.plot(x, y, ls='-', color='k')
ax1.text(np.mean(x), np.max(y), "implicit", rotation=0, ha='center', va='bottom', fontsize='medium')
#ax1.text(x_newaxis*0.9, y_max*0.9, "<- left axis", rotation=0, ha='right', va='bottom', fontsize='medium')
#ax1.text(x_newaxis*1.1, y_max*0.9, "right axis ->", 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) < points:
warn("There were not enough valid entries in your result vector. Reducing the number of points from {0:d} to {1:d}.".format(points, np.sum(mask)))
points = np.sum(mask)
TP["T"] = T_TP[mask]
TP["D"] = D_TP[mask]
TP["H"] = H_TP[mask]
TP["P"] = P_TP[mask]
TP["S"] = S_TP[mask]
# saveNpzData(TP)
else:
TP = None
state = getStateObj(propsfluid)
if backend == "HEOS" and state.has_melting_line():
p_melt = np.logspace(np.log10(state.melting_line(CoolProp.constants.iP_min, CoolProp.constants.iT, 0)), np.log10(state.melting_line(CoolProp.constants.iP_max, CoolProp.constants.iT, 0)), xypoints)
#p_melt = p_range
ML = dict(T=[], D=[], H=[], S=[], P=p_melt)
for p in p_melt:
try:
ML["T"].append(state.melting_line(CoolProp.constants.iT, CoolProp.constants.iP, p))
except Exception as ve:
ML["T"].append(np.inf)
res = PropsSI(["D", "H", "P", "S", "T"], "T", ML["T"], "P", ML["P"], propsfluid)
ML["D"] = res[:, 0]
ML["H"] = res[:, 1]
ML["P"] = res[:, 2]
ML["S"] = res[:, 3]
ML["T"] = res[:, 4]
mask = np.isfinite(ML["T"])
ML["P"] = ML["P"][mask]
ML["T"] = ML["T"][mask]
ML["D"] = ML["D"][mask]
ML["H"] = ML["H"][mask]
ML["S"] = ML["S"][mask]
else:
ML = None
#ML = {}
IP = {}
p_range, T_range = getPTRanges(propsfluid)
critProps = getCriticalProps(propsfluid)
try:
IP["T"] = T_range
IP["P"] = np.zeros_like(T_range) + critProps["P"]
res = PropsSI(["D", "H"], "T", IP["T"], "P", IP["P"], propsfluid)
IP["D"] = res[:, 0]
IP["H"] = res[:, 1]
except Exception as ve:
IP = None
IT = {}
try:
IT["P"] = p_range
IT["T"] = np.zeros_like(p_range) + critProps["T"]
res = PropsSI(["D", "H"], "T", IT["T"], "P", IT["P"], propsfluid)
IT["D"] = res[:, 0]
IT["H"] = res[:, 1]
except Exception as ve:
IT = None
ID = {}
try:
ID["T"] = T_range
ID["D"] = np.zeros_like(p_range) + critProps["D"]
res = PropsSI(["P", "H"], "T", ID["T"], "D", ID["D"], propsfluid)
ID["P"] = res[:, 0]
ID["H"] = res[:, 1]
except Exception as ve:
ID = None
IH = {}
try:
IH["P"] = p_range
IH["H"] = np.zeros_like(p_range) + critProps["H"]
res = PropsSI(["D", "T"], "P", IH["P"], "H", IH["H"], propsfluid)
IH["D"] = res[:, 0]
IH["T"] = res[:, 1]
except Exception as ve:
IH = None
IS = {}
try:
IS["P"] = p_range
IS["S"] = np.zeros_like(p_range) + critProps["S"]
res = PropsSI(["D", "H", "T"], "P", IS["P"], "S", IS["S"], propsfluid)
IS["D"] = res[:, 0]
IS["H"] = res[:, 1]
IS["T"] = res[:, 2]
except Exception as ve:
IS = None
for I in [IP, IT, ID, IH, IS]:
if I is not None:
mask = np.isfinite(I["D"]) & np.isfinite(I["H"])
if np.sum(mask) < 20: I = None
else:
for k in I:
I[k] = I[k][mask]
for inp in getInpList(backend):
if bp is not None:
bp.figure = None
fg = bp.getFigure()
else:
fg = plt.figure()
kind = getTimeKey(inp, getOutList(inp)[-1])
t_data = getSingleData(fld, backend, kind, fluidData)
x_data = getSingleData(fld, backend, "H", fluidData)
y_data = getSingleData(fld, backend, "P", fluidData)
gs = gridspec.GridSpec(1, 2, wspace=None, hspace=None, width_ratios=[10, 1])
ax1 = fg.add_subplot(gs[0, 0], axisbg='Tan')
ax1.set_yscale('log')
#ax2 = ax1.twinx()
minHP = np.min(t_data)
maxHP = np.max(t_data)
minIT = np.percentile(t_data, 10)
maxIT = np.percentile(t_data, 90)
difIT = np.log10(maxIT / minIT) * 0.25
print(kind, ": {0:7.2f} to {1:7.2f}".format(minHP, maxHP))
if kind == "DT":
if bp is not None:
cx1 = bp.getColourMap(reverse=True)
else:
cx1 = mpl.cm.get_cmap('cubehelix_r')
minHP = minIT
maxHP = np.power(10, np.log10(maxIT) + difIT)
#minHP = np.power(10,np.log10(np.percentile(t_data,10)*1e6))
#maxHP = np.power(10,np.log10(np.percentile(t_data,90)*1e6)*1.10)
#maxHP = np.power(10,1.10*np.log10(maxHP))
#minHP = np.percentile(t_data,10)*1e6
#maxHP = np.percentile(t_data,99)*1e6
#print(kind,": {0:7.2f} to {1:7.2f}".format(minHP,maxHP))
#minHP = 100
#maxHP = 20000
else:
if bp is not None:
cx1 = bp.getColourMap()
else:
cx1 = mpl.cm.get_cmap('cubehelix')
minHP = np.power(10, np.log10(minIT) - difIT)
maxHP = maxIT
#minHP = np.power(10,np.log10(np.percentile(t_data,10)*1e6)*0.90)
#maxHP = np.power(10,np.log10(np.percentile(t_data,90)*1e6))
# minHP = np.percentile(t_data,01)*1e6
#maxHP = np.percentile(t_data,90)*1e6
#print(kind,": {0:7.2f} to {1:7.2f}".format(minHP,maxHP))
#minHP = 100
#maxHP = 20000
#cx1_r = reverse_colourmap(cx1)
cNorm = mpl.colors.LogNorm(vmin=minHP, vmax=maxHP)
#cNorm = mpl.colors.LogNorm(vmin=ceil(minHP/1e1)*1e1, vmax=floor(maxHP/1e2)*1e2)
#cNorm = mpl.colors.Normalize(vmin=round(minHP,-2), vmax=round(maxHP,-2))
colourSettings = dict(c=t_data, edgecolors='none', cmap=cx1, norm=cNorm)
pointSettings = dict(s=6)
scatterSettings = dict(rasterized=True, alpha=0.5)
#scatterSettings = dict(rasterized=False, alpha=0.5)
scatterSettings.update(colourSettings)
scatterSettings.update(pointSettings)
SC = ax1.scatter(x_data / 1e6, y_data / 1e5, **scatterSettings)
for I in [TP, ML]:
if I is not None:
ax1.plot(I["H"] / 1e6, I["P"] / 1e5, lw=1.5, c='k')
for I in [IP, IT, ID, IS, IH]:
if I is not None:
ax1.plot(I["H"] / 1e6, I["P"] / 1e5, lw=1.0, c='k', alpha=1)
# ax1.set_xlim([0e+0,6e1])
# ax1.set_ylim([5e-1,2e4])
ax1.set_xlim([np.percentile(x_data / 1e6, 0.1), np.percentile(x_data / 1e6, 99.9)])
ax1.set_ylim([np.percentile(y_data / 1e5, 0.1), np.percentile(y_data / 1e5, 99.9)])
formatter = ticker.LogFormatter(base=10.0, labelOnlyBase=False)
#formatter = ticker.ScalarFormatter()
#ticks = roundList(np.logspace(np.log10(ax1.get_ylim()[0]), np.log10(ax1.get_ylim()[1]), 5))
#locator = ticker.FixedLocator(ticks)
# ax1.yaxis.set_major_locator(locator)
ax1.yaxis.set_major_formatter(formatter)
cax = fg.add_subplot(gs[0, 1])
formatter = ticker.ScalarFormatter()
CB = fg.colorbar(SC, cax=cax, format=formatter)
CB.set_alpha(1)
CB.locator = ticker.MaxNLocator(nbins=7)
#ticks = roundList(np.logspace(np.log10(minHP), np.log10(maxHP), 5))
#CB.locator = ticker.FixedLocator(ticks)
CB.update_ticks()
CB.draw_all()
# fg.suptitle("f("+inp+")-"+backend.lower()+"-"+fld.lower())
CB.set_label(backend.upper() + "::" + fld + ', execution time per f(' + inp[0] + "," + inp[1] + ') call (us)')
ax1.set_xlabel(r'Specific enthalpy (MJ/kg)')
ax1.set_ylabel(r'Pressure (bar)')
fg.tight_layout()
fg.savefig(path.join(getFigureFolder(), "TimeComp-" + inp + "-" + backend.lower() + "-" + fld.lower() + ".pdf"))
#CB.set_label(r'Execution time per call (\si{\us})')
#ax1.set_xlabel(r'Specific enthalpy (\si{\mega\J\per\kg})')
#ax1.set_ylabel(r'Pressure (\si{\bar})')
# fg.tight_layout()
# bp.savepgf(path.join(getFigureFolder(),"TimeComp-"+inp+"-"+backend.lower()+"-"+fld.lower()+".pgf"),fg,repList)
plt.close()
except Exception as e:
print(e)
pass
plt.close('all')
Reference in New Issue View Git Blame Copy Permalink