github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-01-24 13:28:01 -05:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'draw_process'

Browse Source
This commit is contained in:
Jorrit Wronski
2015-10-08 18:22:24 +02:00
parent 43a78b1ec8 6e705f902f
commit 6475d35538
10 changed files with 2548 additions and 821 deletions
Download Patch File Download Diff File Expand all files Collapse all files

0
dev/TTSE/check_TTSE.py → dev/TTSE/check_TTSE_old.py
View File

322
dev/TTSE/check_TTSE_v4.py Normal file
View File

@@ -0,0 +1,322 @@
import CoolProp.CoolProp as CP
import matplotlib
matplotlib.rc('font', family='serif', serif='Times New Roman')
#from matplotlib2tikz import save as tikz_save
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib.cm as cmx
import matplotlib.ticker
from matplotlib.patches import Ellipse
from matplotlib.transforms import ScaledTranslation
import numpy as np
import random
from numpy import linspace, meshgrid
from matplotlib.mlab import griddata
from matplotlib.gridspec import GridSpec
# Create the colourmap
#import numpy as np
#import matplotlib.pyplot as plt
import matplotlib._cm, matplotlib.cm
specs = matplotlib._cm.cubehelix(gamma=1.4,s=0.4,r=-0.8,h=2.0)
specs_r = matplotlib.cm._reverse_cmap_spec(specs)
matplotlib.cm.register_cmap(name="jorrithelix" , data=specs)
matplotlib.cm.register_cmap(name="jorrithelix"+"_r", data=specs_r)
def makeGrid(x, y, z, resX=200, resY=200):
"Convert 3 column data to matplotlib grid"
xi = linspace(min(x), max(x), resX)
yi = linspace(min(y), max(y), resY)
Z = griddata(x, y, z, xi, yi)
X, Y = meshgrid(xi, yi)
return X, Y, Z
def getErrors(p, h, out='D', Ref=''):
"Get the relative errors from table-based interpolation"
errorTTSE = 1e3
errorBICUBIC = 1e3
try:
# Using the EOS
CP.disable_TTSE_LUT(Ref)
EOS = CP.PropsSI(out,'P',p,'H',h,Ref)
# Using the TTSE method
CP.enable_TTSE_LUT(Ref)
CP.set_TTSE_mode(Ref,"TTSE")
TTSE = CP.PropsSI(out,'P',p,'H',h,Ref)
# Using the Bicubic method
CP.enable_TTSE_LUT(Ref)
CP.set_TTSE_mode(Ref,"BICUBIC")
BICUBIC = CP.PropsSI(out,'P',p,'H',h,Ref)
errorTTSE = abs(TTSE /EOS-1.0)*100.0
errorBICUBIC = abs(BICUBIC/EOS-1.0)*100.0
except ValueError as VE:
print VE
pass
return errorTTSE,errorBICUBIC
#['YlOrRd', 'PuBuGn', 'hot', 'cubehelix', 'gnuplot', 'gnuplot2']:
for colourmap in ['jorrithelix']:
for out in ['D']:
## landscape figure
#fig = plt.figure(figsize=(10,5))
#ax1 = fig.add_axes((0.08,0.1,0.32,0.83))
#ax2 = fig.add_axes((0.50,0.1,0.32,0.83))
#cbar_ax = fig.add_axes([0.80, 0.075, 0.05, 0.875])
# portrait figure
#fig = plt.figure(figsize=(5,8))
#ax1 = plt.subplot2grid((2,8), (0,0), colspan=7)
#ax2 = plt.subplot2grid((2,8), (1,0), colspan=7)
#cbar_ax = plt.subplot2grid((2,8), (0,7), colspan=1, rowspan=2)
#fig = plt.figure(figsize=(8,4))
#ax1 = plt.subplot2grid((1,7), (0,0), colspan=3)
#ax2 = plt.subplot2grid((1,7), (0,3), colspan=3)
#cbar_ax = plt.subplot2grid((1,7), (0,6), colspan=1, rowspan=1)
#plt.tight_layout()
fig = plt.figure(figsize=(8,4))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#cbar_ax = plt.subplot2grid((1,7), (0,6), colspan=1, rowspan=1)
#plt.tight_layout()
#Ref = 'R245fa'
#Ref = 'Isopentane'
Ref = 'Air'
T = np.linspace(CP.PropsSI(Ref,'Tmin')+0.1,CP.PropsSI(Ref,'Tcrit')-0.01,300)
pV = CP.PropsSI('P','T',T,'Q',1,Ref)
hL = CP.PropsSI('H','T',T,'Q',0,Ref)
hV = CP.PropsSI('H','T',T,'Q',1,Ref)
hTP= np.append(hL,[hV[::-1]])
pTP= np.append(pV,[pV[::-1]])
HHH1, PPP1, EEE1 = [], [], []
HHH2, PPP2, EEE2 = [], [], []
cNorm = colors.LogNorm(vmin=1e-10, vmax=1e-1)
scalarMap = cmx.ScalarMappable(norm = cNorm, cmap = plt.get_cmap(colourmap))
# Setting the limits for enthalpy and pressure
p_min = CP.PropsSI(Ref,'ptriple')
p_max = 60e5
h_min = CP.PropsSI('H','T',CP.PropsSI(Ref,'Ttriple')+0.5,'Q',0,Ref)
h_max = CP.PropsSI('H','T',500+273.15,'P',p_max,Ref)
# Creating some isotherms for better illustration of the cycle
isoT = np.array([0,100,200,300,400])+273.15
isoP = np.logspace(np.log10(p_min),np.log10(p_max),base=10)
ones = np.ones(isoP.shape)
isoH = [ CP.PropsSI('H','T',T*ones,'P',isoP,Ref) for T in isoT ]
print "Lower left and upper right coordinates: ({0},{1}), ({2},{3})".format(h_min,p_min,h_max,p_max)
CP.set_TTSESinglePhase_LUT_range(Ref,h_min,h_max*1.05,p_min,p_max*1.05)
for a_useless_counter in range(40000):
h = random.uniform(h_min,h_max)
p = 10**random.uniform(np.log10(p_min),np.log10(p_max))
try:
# Using the EOS
CP.disable_TTSE_LUT(Ref)
rhoEOS = CP.PropsSI('D','P',p,'H',h,Ref)
TEOS = CP.PropsSI('T','P',p,'H',h,Ref)
if out =='C': cpEOS = CP.PropsSI('C','P',p,'H',h,Ref)
# Using the TTSE method
CP.enable_TTSE_LUT(Ref)
CP.set_TTSE_mode(Ref,"TTSE")
rhoTTSE = CP.PropsSI('D','P',p,'H',h,Ref)
TTTSE = CP.PropsSI('T','P',p,'H',h,Ref)
if out =='C': cpTTSE = CP.PropsSI('C','P',p,'H',h,Ref)
# Using the Bicubic method
CP.enable_TTSE_LUT(Ref)
CP.set_TTSE_mode(Ref,"BICUBIC")
rhoBICUBIC = CP.PropsSI('D','P',p,'H',h,Ref)
TBICUBIC = CP.PropsSI('T','P',p,'H',h,Ref)
if out =='C': cpBICUBIC = CP.PropsSI('C','P',p,'H',h,Ref)
if out == 'D':
errorTTSE = abs(rhoTTSE/rhoEOS-1)*100
errorBICUBIC = abs(rhoBICUBIC/rhoEOS-1)*100
elif out == 'T':
errorTTSE = abs(TTTSE/TEOS-1)*100
errorBICUBIC = abs(TBICUBIC/TEOS-1)*100
elif out == 'C':
errorTTSE = abs(cpTTSE/cpEOS-1)*100
errorBICUBIC = abs(cpBICUBIC/cpEOS-1)*100
HHH1.append(h)
PPP1.append(p)
EEE1.append(errorTTSE)
HHH2.append(h)
PPP2.append(p)
EEE2.append(errorBICUBIC)
except ValueError as VE:
#print VE
pass
HHH1 = np.array(HHH1)
PPP1 = np.array(PPP1)
SC1 = ax1.scatter(HHH1/1e3, PPP1/1e5, s=8, c=EEE1, edgecolors = 'none', cmap = plt.get_cmap(colourmap), norm = cNorm, rasterized=True)
#X, Y, Z = makeGrid(HHH1, np.log10(PPP1), EEE1)
#SC1 = matplotlib.pyplot.contourf(X, Y, Z,
# alpha=0.75,
# norm=cNorm,
# cmap=matplotlib.pyplot.get_cmap(colourmap))#,
# #rasterized=True)
HHH2 = np.array(HHH2)
PPP2 = np.array(PPP2)
SC2 = ax2.scatter(HHH2/1e3, PPP2/1e5, s=8, c=EEE2, edgecolors = 'none', cmap = plt.get_cmap(colourmap), norm = cNorm, rasterized=True)
if out == 'D':
ax1.set_title('rel. density error, TTSE')
ax2.set_title('rel. density error, bicubic')
elif out == 'T':
ax1.set_title('rel. temperature error, TTSE')
ax2.set_title('rel. temperature error, bicubic')
elif out == 'C':
ax1.set_title('rel. heat capacity error, TTSE')
ax2.set_title('rel. heat capacity error, bicubic')
for ax in [ax1, ax2]:
#h_min = np.ceil(h_min)
delta = 0.1
delta_min = 1.0+delta
delta_max = 1.0-delta
#ax.set_xlim(delta_min*h_min/1e3, delta_max*h_max/1e3)
#ax.set_ylim(delta_min*p_min/1e5, delta_max*p_max/1e5)
ax.set_xlim(-155, 800)
ax.set_ylim(0.025, 58)
ax.set_yscale('log')
#ticks = np.array([0.02,0.05,0.1,0.2,0.5,1,2,5,10,20,50])
ticks = np.array([0.05,0.1,0.2,0.5,1,2,5,10,20,50])
labels = [str(tick) for tick in ticks]
ax.set_yticks(ticks)
ax.set_yticklabels(labels)
ax.get_yaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
#ticks = [150,250,350,450,550]
#labels = [str(tick) for tick in ticks]
#ax.set_xticks(ticks)
#ax.set_xticklabels(labels)
#ax.tick_params(axis='y',which='minor',left='off')
#ax.set_xlabel('Enthalpy [kJ \cdot kg^{-1}]')
ax.set_xlabel('Specific Enthalpy [kJ$\cdot$kg$\mathdefault{^{-1}\!}$]')
ax.set_ylabel('Pressure [bar]')
#ax.plot(hL/1e3,pV/1e5,'k',lw = 4)
#ax.plot(hV/1e3,pV/1e5,'k',lw = 4)
ax.plot(hTP/1e3,pTP/1e5,'k',lw = 3)
for i,T in enumerate(isoT):
ax.plot(isoH[i]/1e3,isoP/1e5,'k',lw = 1)
#CB = fig.colorbar(SC1)
#cbar_ax = fig.add_axes([0.80, 0.075, 0.05, 0.875])
#CB = fig.colorbar(SC1, cax=cbar_ax)
#CB = matplotlib.pyplot.colorbar(SC2)
#CB.solids.set_rasterized(True)
#ax2.yaxis.set_visible(False)
#[x0,y0,width,height]
#cbar_ax = fig.add_axes([0.95, 0.00, 0.05, 1.00])
#CB = fig.colorbar(SC2, ax=[ax1,ax2], cax=cbar_ax)
#CB.solids.set_rasterized(True)
#from mpl_toolkits.axes_grid1 import make_axes_locatable
#divider = make_axes_locatable(ax2)
#cbar_ax = divider.append_axes("right", "5%", pad="0%")
#CB = plt.colorbar(SC2, cax=cbar_ax)
#CB.solids.set_rasterized(True)
#CB = fig.colorbar(SC2)
#CB.solids.set_rasterized(True)
from mpl_toolkits.axes_grid1 import make_axes_locatable
divider = make_axes_locatable(ax2)
ax_cb = divider.new_horizontal(size="5%", pad=0.05)
#fig1 = ax.get_figure()
fig.add_axes(ax_cb)
CB = fig.colorbar(SC2, cax=ax_cb)
#aspect = 5./2.
#ax1.set_aspect(aspect)
#ax2.set_aspect(aspect)
CB.solids.set_rasterized(True)
if out == 'D':
CB.set_label(r'$\|\rho/\rho\mathdefault{_{EOS}-1\|\times 100}$ [%]')
elif out == 'T':
CB.set_label(r'$\|T/T\mathdefault{_{EOS}-1\|\times 100}$ [%]')
elif out == 'C':
CB.set_label(r'$\|c\mathdefault{_p}/c\mathdefault{_{p,EOS}-1\|\times 100}$ [%]')
# The plot is finished, now we add an ellipse
#circle=plt.Circle((5,5),.5,color='b',fill=False)
#A scale-free ellipse.
#xy - center of ellipse
#width - total length (diameter) of horizontal axis
#height - total length (diameter) of vertical axis
#angle - rotation in degrees (anti-clockwise)
p_op_min = 1e5
p_op_max = 3e5
h_op_min = CP.PropsSI('H','T',400+273.15,'P',p_op_max,Ref)
h_op_max = CP.PropsSI('H','T', 25+273.15,'P',p_op_max,Ref)
p_op_cen = (p_op_min + p_op_max) / 2.0
h_op_cen = (h_op_min + h_op_max) / 2.0
p_op_hei = p_op_max - p_op_min
h_op_wid = h_op_max - h_op_min
#for ax in [ax1, ax2]:
##x,y = 10,0
### use the axis scale tform to figure out how far to translate
##circ_offset = ScaledTranslation(x,y,ax.transScale)
### construct the composite tform
##circ_tform = circ_offset + ax.transLimits + ax.transAxes
#ellipse = Ellipse(xy=(h_op_cen,p_op_cen), width=h_op_wid, height=p_op_hei, angle=15, color='black')#, transform=circ_tform)
#ax.add_artist(ellipse)
# font_def = font_manager.FontProperties(family='Helvetica', style='normal',
# size=sizeOfFont, weight='normal', stretch='normal')
#
# for a in fig.axes:
# for label in [a.get_xticklabels(), a.get_yticklabels()]:
# label.set_fontproperties(ticks_font
#plt.savefig(out+'_'+colourmap+'_TTSE_BICUBIC.png', dpi = 300, transparent = True)
#plt.savefig(out+'_'+colourmap+'_TTSE_BICUBIC.eps')
# plt.savefig(out+'_'+colourmap+'_TTSE_BICUBIC.pdf')
plt.tight_layout()
plt.savefig('check_TTSE_'+colourmap+'.pdf' )
#tikz_save( 'check_TTSE.tikz')
#plt.savefig(out+'_'+colourmap+'_TTSE_BICUBIC.jpg', dpi = 1200)
plt.close()

276
dev/TTSE/check_TTSE_v5.py Normal file
View File

@@ -0,0 +1,276 @@
import CoolProp.CoolProp as CP
import matplotlib
matplotlib.rc('font', family='serif', serif='Times New Roman')
#from matplotlib2tikz import save as tikz_save
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib.cm as cmx
import matplotlib.ticker
from matplotlib.patches import Ellipse
from matplotlib.transforms import ScaledTranslation
import numpy as np
import random
from numpy import linspace, meshgrid
from matplotlib.mlab import griddata
from matplotlib.gridspec import GridSpec
def makeGrid(x, y, z, resX=200, resY=200):
"Convert 3 column data to matplotlib grid"
xi = linspace(min(x), max(x), resX)
yi = linspace(min(y), max(y), resY)
Z = griddata(x, y, z, xi, yi)
X, Y = meshgrid(xi, yi)
return X, Y, Z
def getErrors(p, h, out='D', Ref=''):
"Get the relative errors from table-based interpolation"
errorTTSE = 1e3
errorBICUBIC = 1e3
try:
# Using the EOS
CP.disable_TTSE_LUT(Ref)
EOS = CP.PropsSI(out,'P',p,'H',h,Ref)
# Using the TTSE method
CP.enable_TTSE_LUT(Ref)
CP.set_TTSE_mode(Ref,"TTSE")
TTSE = CP.PropsSI(out,'P',p,'H',h,Ref)
# Using the Bicubic method
CP.enable_TTSE_LUT(Ref)
CP.set_TTSE_mode(Ref,"BICUBIC")
BICUBIC = CP.PropsSI(out,'P',p,'H',h,Ref)
errorTTSE = abs(TTSE /EOS-1.0)*100.0
errorBICUBIC = abs(BICUBIC/EOS-1.0)*100.0
except ValueError as VE:
print VE
pass
return errorTTSE,errorBICUBIC
#['YlOrRd', 'PuBuGn', 'hot', 'cubehelix', 'gnuplot', 'gnuplot2']:
for colourmap in ['cubehelix']:
for out in ['D']:
## landscape figure
#fig = plt.figure(figsize=(10,5))
#ax1 = fig.add_axes((0.08,0.1,0.32,0.83))
#ax2 = fig.add_axes((0.50,0.1,0.32,0.83))
#cbar_ax = fig.add_axes([0.80, 0.075, 0.05, 0.875])
# portrait figure
fig = plt.figure(figsize=(5,8))
#ax1 = fig.add_axes((0.175,0.575,0.575,0.375))
#ax2 = fig.add_axes((0.175,0.075,0.575,0.375))
#cbar_ax = fig.add_axes([0.80, 0.075, 0.05, 0.875])
#ax1 = fig.add_subplot(241,colspan=3)
#ax2 = fig.add_subplot(245,colspan=3)
#cbar_ax = fig.add_subplot(244,rowspan=2)
ax1 = plt.subplot2grid((2,8), (0,0), colspan=7)
ax2 = plt.subplot2grid((2,8), (1,0), colspan=7)
cbar_ax = plt.subplot2grid((2,8), (0,7), colspan=1, rowspan=2)
#Ref = 'R245fa'
#Ref = 'Isopentane'
Ref = 'Air'
T = np.linspace(CP.PropsSI(Ref,'Tmin')+0.1,CP.PropsSI(Ref,'Tcrit')-0.01,300)
pV = CP.PropsSI('P','T',T,'Q',1,Ref)
hL = CP.PropsSI('H','T',T,'Q',0,Ref)
hV = CP.PropsSI('H','T',T,'Q',1,Ref)
hTP= np.append(hL,[hV[::-1]])
pTP= np.append(pV,[pV[::-1]])
HHH1, PPP1, EEE1 = [], [], []
HHH2, PPP2, EEE2 = [], [], []
cNorm = colors.LogNorm(vmin=1e-10, vmax=1e-1)
scalarMap = cmx.ScalarMappable(norm = cNorm, cmap = plt.get_cmap(colourmap))
# Setting the limits for enthalpy and pressure
p_min = CP.PropsSI(Ref,'ptriple')
p_max = 60e5
h_min = CP.PropsSI('H','T',CP.PropsSI(Ref,'Ttriple')+0.5,'Q',0,Ref)
h_max = CP.PropsSI('H','T',500+273.15,'P',p_max,Ref)
# Creating some isotherms for better illustration of the cycle
isoT = np.array([0,100,200,300,400])+273.15
isoP = np.logspace(np.log10(p_min),np.log10(p_max),base=10)
ones = np.ones(isoP.shape)
isoH = [ CP.PropsSI('H','T',T*ones,'P',isoP,Ref) for T in isoT ]
print "Lower left and upper right coordinates: ({0},{1}), ({2},{3})".format(h_min,p_min,h_max,p_max)
CP.set_TTSESinglePhase_LUT_range(Ref,h_min,h_max*1.05,p_min,p_max*1.05)
for a_useless_counter in range(40000):
h = random.uniform(h_min,h_max)
p = 10**random.uniform(np.log10(p_min),np.log10(p_max))
try:
# Using the EOS
CP.disable_TTSE_LUT(Ref)
rhoEOS = CP.PropsSI('D','P',p,'H',h,Ref)
TEOS = CP.PropsSI('T','P',p,'H',h,Ref)
if out =='C': cpEOS = CP.PropsSI('C','P',p,'H',h,Ref)
# Using the TTSE method
CP.enable_TTSE_LUT(Ref)
CP.set_TTSE_mode(Ref,"TTSE")
rhoTTSE = CP.PropsSI('D','P',p,'H',h,Ref)
TTTSE = CP.PropsSI('T','P',p,'H',h,Ref)
if out =='C': cpTTSE = CP.PropsSI('C','P',p,'H',h,Ref)
# Using the Bicubic method
CP.enable_TTSE_LUT(Ref)
CP.set_TTSE_mode(Ref,"BICUBIC")
rhoBICUBIC = CP.PropsSI('D','P',p,'H',h,Ref)
TBICUBIC = CP.PropsSI('T','P',p,'H',h,Ref)
if out =='C': cpBICUBIC = CP.PropsSI('C','P',p,'H',h,Ref)
if out == 'D':
errorTTSE = abs(rhoTTSE/rhoEOS-1)*100
errorBICUBIC = abs(rhoBICUBIC/rhoEOS-1)*100
elif out == 'T':
errorTTSE = abs(TTTSE/TEOS-1)*100
errorBICUBIC = abs(TBICUBIC/TEOS-1)*100
elif out == 'C':
errorTTSE = abs(cpTTSE/cpEOS-1)*100
errorBICUBIC = abs(cpBICUBIC/cpEOS-1)*100
HHH1.append(h)
PPP1.append(p)
EEE1.append(errorTTSE)
HHH2.append(h)
PPP2.append(p)
EEE2.append(errorBICUBIC)
except ValueError as VE:
#print VE
pass
HHH1 = np.array(HHH1)
PPP1 = np.array(PPP1)
SC1 = ax1.scatter(HHH1/1e3, PPP1/1e5, s=8, c=EEE1, edgecolors = 'none', cmap = plt.get_cmap(colourmap), norm = cNorm, rasterized=True)
#X, Y, Z = makeGrid(HHH1, np.log10(PPP1), EEE1)
#SC1 = matplotlib.pyplot.contourf(X, Y, Z,
# alpha=0.75,
# norm=cNorm,
# cmap=matplotlib.pyplot.get_cmap(colourmap))#,
# #rasterized=True)
HHH2 = np.array(HHH2)
PPP2 = np.array(PPP2)
SC2 = ax2.scatter(HHH2/1e3, PPP2/1e5, s=8, c=EEE2, edgecolors = 'none', cmap = plt.get_cmap(colourmap), norm = cNorm, rasterized=True)
if out == 'D':
ax1.set_title('rel. density error, TTSE')
ax2.set_title('rel. density error, bicubic')
elif out == 'T':
ax1.set_title('rel. temperature error, TTSE')
ax2.set_title('rel. temperature error, bicubic')
elif out == 'C':
ax1.set_title('rel. heat capacity error, TTSE')
ax2.set_title('rel. heat capacity error, bicubic')
for ax in [ax1, ax2]:
#h_min = np.ceil(h_min)
delta = 0.1
delta_min = 1.0+delta
delta_max = 1.0-delta
#ax.set_xlim(delta_min*h_min/1e3, delta_max*h_max/1e3)
#ax.set_ylim(delta_min*p_min/1e5, delta_max*p_max/1e5)
ax.set_xlim(-155, 800)
ax.set_ylim(0.025, 58)
ax.set_yscale('log')
#ticks = np.array([0.02,0.05,0.1,0.2,0.5,1,2,5,10,20,50])
ticks = np.array([0.05,0.1,0.2,0.5,1,2,5,10,20,50])
labels = [str(tick) for tick in ticks]
ax.set_yticks(ticks)
ax.set_yticklabels(labels)
ax.get_yaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
#ticks = [150,250,350,450,550]
#labels = [str(tick) for tick in ticks]
#ax.set_xticks(ticks)
#ax.set_xticklabels(labels)
#ax.tick_params(axis='y',which='minor',left='off')
#ax.set_xlabel('Enthalpy [kJ \cdot kg^{-1}]')
ax.set_xlabel('Specific Enthalpy [kJ$\cdot$kg$\mathdefault{^{-1}\!}$]')
ax.set_ylabel('Pressure [bar]')
#ax.plot(hL/1e3,pV/1e5,'k',lw = 4)
#ax.plot(hV/1e3,pV/1e5,'k',lw = 4)
ax.plot(hTP/1e3,pTP/1e5,'k',lw = 3)
for i,T in enumerate(isoT):
ax.plot(isoH[i]/1e3,isoP/1e5,'k',lw = 1)
#CB = fig.colorbar(SC1)
#cbar_ax = fig.add_axes([0.80, 0.075, 0.05, 0.875])
CB = fig.colorbar(SC1, cax=cbar_ax)
if out == 'D':
CB.set_label(r'$\|\rho/\rho\mathdefault{_{EOS}-1\|\times 100}$ [%]')
elif out == 'T':
CB.set_label(r'$\|T/T\mathdefault{_{EOS}-1\|\times 100}$ [%]')
elif out == 'C':
CB.set_label(r'$\|c\mathdefault{_p}/c\mathdefault{_{p,EOS}-1\|\times 100}$ [%]')
# The plot is finished, now we add an ellipse
#circle=plt.Circle((5,5),.5,color='b',fill=False)
#A scale-free ellipse.
#xy - center of ellipse
#width - total length (diameter) of horizontal axis
#height - total length (diameter) of vertical axis
#angle - rotation in degrees (anti-clockwise)
p_op_min = 1e5
p_op_max = 3e5
h_op_min = CP.PropsSI('H','T',400+273.15,'P',p_op_max,Ref)
h_op_max = CP.PropsSI('H','T', 25+273.15,'P',p_op_max,Ref)
p_op_cen = (p_op_min + p_op_max) / 2.0
h_op_cen = (h_op_min + h_op_max) / 2.0
p_op_hei = p_op_max - p_op_min
h_op_wid = h_op_max - h_op_min
#for ax in [ax1, ax2]:
##x,y = 10,0
### use the axis scale tform to figure out how far to translate
##circ_offset = ScaledTranslation(x,y,ax.transScale)
### construct the composite tform
##circ_tform = circ_offset + ax.transLimits + ax.transAxes
#ellipse = Ellipse(xy=(h_op_cen,p_op_cen), width=h_op_wid, height=p_op_hei, angle=15, color='black')#, transform=circ_tform)
#ax.add_artist(ellipse)
# font_def = font_manager.FontProperties(family='Helvetica', style='normal',
# size=sizeOfFont, weight='normal', stretch='normal')
#
# for a in fig.axes:
# for label in [a.get_xticklabels(), a.get_yticklabels()]:
# label.set_fontproperties(ticks_font
#plt.savefig(out+'_'+colourmap+'_TTSE_BICUBIC.png', dpi = 300, transparent = True)
#plt.savefig(out+'_'+colourmap+'_TTSE_BICUBIC.eps')
# plt.savefig(out+'_'+colourmap+'_TTSE_BICUBIC.pdf')
plt.tight_layout()
plt.savefig('check_TTSE.pdf' )
#tikz_save( 'check_TTSE.tikz')
#plt.savefig(out+'_'+colourmap+'_TTSE_BICUBIC.jpg', dpi = 1200)
plt.close()

63
dev/Tickets/351.py Normal file
View File

@@ -0,0 +1,63 @@
from CoolProp.CoolProp import PropsSI
from CoolProp.Plots import PropertyPlot
import matplotlib.pylab as pl
from numpy import *
import CoolProp
eta_e_s = 0.88
eta_p_s = 0.5
T_max = 550 + 273.15
p_max = 13000.e3
p_cd = 5.e3
def SimpleRankineCycle(T3, p3, p1, epsilon_e, epsilon_p, fluid='water'):
h1 = PropsSI('H', 'P', p1, 'Q', 0., fluid)
s1 = PropsSI('S', 'P', p1, 'Q', 0., fluid)
T1 = PropsSI('T', 'P', p1, 'Q', 0., fluid)
p2 = p3
h2 = h1 + (PropsSI('H', 'P', p2, 'S', s1, fluid) - h1) / epsilon_p
s2 = PropsSI('S', 'H', h2, 'P', p2, fluid)
T2 = PropsSI('T', 'H', h2, 'P', p2, fluid)
h3 = PropsSI('H', 'P', p3, 'T', T3, fluid)
s3 = PropsSI('S', 'H', h3, 'P', p3, fluid)
p4 = p1
h4 = h3 - epsilon_e * (h3 - PropsSI('H', 'P', p4, 'S', s3, fluid))
s4 = PropsSI('S', 'H', h4, 'P', p4, fluid)
T4 = PropsSI('T', 'H', h4, 'P', p4, fluid)
w_net = h3 - h4
q_boiler = h3 - h2
eta_c = w_net / q_boiler
Ts = PropertyPlot(fluid, 'Ts', 'KSI')
Ts.set_axis_limits([0., 12., 200., 900.])
Ts.calc_isolines(CoolProp.iP, [Ts.system.P.from_SI(p1), Ts.system.P.from_SI(p3)], num=10)
Ts.calc_isolines(CoolProp.iQ, [Ts.system.Q.from_SI(0.), Ts.system.Q.from_SI(1.)], num=11)
Ts.draw_isolines()
states = zip(Ts.system.S.from_SI(array([s1,s2,s3,s4,s1])),Ts.system.T.from_SI(array([T1,T2,T3,T4,T1])))
Ts.draw_process(states, iso_types=None, line_opts={'color':'red', 'lw':1.5})
isot = [
None, # non-isentropic pumping from 1 to 2
CoolProp.iP, # p2=p3
None,
CoolProp.iP, # p4=p1
]
Ts.draw_process(states, iso_types=isot, line_opts={'color':'green', 'lw':1.5})
ax = Ts.axis
ax.text(Ts.system.S.from_SI(s1), Ts.system.T.from_SI(T1), ' 1', fontsize=10, rotation=0, color='r')
ax.text(Ts.system.S.from_SI(s2), Ts.system.T.from_SI(T2), ' 2', fontsize=10, rotation=0, color='r')
ax.text(Ts.system.S.from_SI(s3), Ts.system.T.from_SI(T3), ' 3', fontsize=10, rotation=0, color='r')
ax.text(Ts.system.S.from_SI(s4), Ts.system.T.from_SI(T4), ' 4', fontsize=10, rotation=0, color='r')
ax.text(Ts.system.S.from_SI(8e3),Ts.system.T.from_SI(850),"Efficiency: %.1f%%" %(eta_c*100.))
ax.text(Ts.system.S.from_SI(8e3),Ts.system.T.from_SI(800),"Net work: %d kJ/kg" %(w_net/1000))
ax.text(Ts.system.S.from_SI(8e3),Ts.system.T.from_SI(750),"Heat input: %d kJ/kg" %(q_boiler/1000))
return Ts
Ts = SimpleRankineCycle(T_max, p_max, p_cd, eta_e_s, eta_p_s, fluid="water")
Ts.savefig('ticket-351.pdf')

26
wrappers/Python/CoolProp/CoolProp.pyx
View File

@@ -481,6 +481,32 @@ cpdef int get_debug_level():
# return False
#
cpdef extract_backend(string in_str):
"""
A Python wrapper of C++ function :cpapi:`CoolProp::extract_backend` .
"""
cdef string bck, fld
# Extract the backend and the fluid from the input string
_extract_backend(in_str, bck, fld)
return bck, fld
cpdef extract_fractions(string flds):
"""
A Python wrapper of C++ function :cpapi:`CoolProp::extract_fractions` .
"""
cdef vector[double] frcs
cdef string del_flds
# Extract the fractions
#frcs.clear()
frcs.push_back(1.0)
del_flds = _extract_fractions(flds, frcs)
# Extract the fluids
fluids = del_flds.split('&')
return fluids,frcs
cdef toSI(constants_header.parameters key, double val):
"""
Convert a value in kSI system to SI system (supports a limited subset of variables)

1197
wrappers/Python/CoolProp/Plots/Common.py
View File

File diff suppressed because it is too large Load Diff

879
wrappers/Python/CoolProp/Plots/Plots.py
View File

@@ -7,257 +7,96 @@ from scipy.interpolate import interp1d
import CoolProp.CoolProp as CP
from .Common import BasePlot
from scipy import interpolate
from scipy.spatial.kdtree import KDTree
class IsoLine(object):
def __init__(self):
self.DEBUG = False
# direct geometry
self.X = None #
self.Y = None #
self.type = None #
self.value = None #
self.unit = None #
self.opts = None #
import warnings
from CoolProp.Plots.Common import IsoLine,BasePlot
import CoolProp
import sys
from CoolProp.Plots.SimpleCycles import StatePoint, StateContainer,\
SimpleRankineCycle
def InlineLabel(xv,yv,x = None, y= None, axis = None, fig = None):
"""
This will give the coordinates and rotation required to align a label with
a line on a plot
"""
def ToPixelCoords(xv,yv,axis,fig):
[Axmin,Axmax]=axis.get_xlim()
[Aymin,Aymax]=axis.get_ylim()
DELTAX_axis=Axmax-Axmin
DELTAY_axis=Aymax-Aymin
width=fig.get_figwidth()
height=fig.get_figheight()
pos=axis.get_position().get_points()
[[Fxmin,Fymin],[Fxmax,Fymax]]=pos
DELTAX_fig=width*(Fxmax-Fxmin)
DELTAY_fig=height*(Fymax-Fymin)
#Convert coords to pixels
x=(xv-Axmin)/DELTAX_axis*DELTAX_fig+Fxmin
y=(yv-Aymin)/DELTAY_axis*DELTAY_fig+Fymin
return x,y
def ToDataCoords(xv,yv,axis,fig):
[Axmin,Axmax]=axis.get_xlim()
[Aymin,Aymax]=axis.get_ylim()
DELTAX_axis=Axmax-Axmin
DELTAY_axis=Aymax-Aymin
width=fig.get_figwidth()
height=fig.get_figheight()
pos=axis.get_position().get_points()
[[Fxmin,Fymin],[Fxmax,Fymax]]=pos
DELTAX_fig=(Fxmax-Fxmin)*width
DELTAY_fig=(Fymax-Fymin)*height
#Convert back to measurements
x=(xv-Fxmin)/DELTAX_fig*DELTAX_axis+Axmin
y=(yv-Fymin)/DELTAY_fig*DELTAY_axis+Aymin
return x,y
def get_x_y_dydx(xv,yv,x):
"""Get x and y coordinates and the linear interpolation derivative"""
# Old implementation:
##Get the rotation angle
#f = interp1d(xv, yv)
#y = f(x)
#h = 0.00001*x
#dy_dx = (f(x+h)-f(x-h))/(2*h)
#return x,y,dy_dx
if len(xv)==len(yv)>1: # assure same length
if len(xv)==len(yv)==2: # only two points
if numpy.min(xv)<x<numpy.max(xv):
dx = xv[1] - xv[0]
dy = yv[1] - yv[0]
dydx = dy/dx
y = yv[0] + dydx * (x-xv[0])
return x,y,dydx
else:
raise ValueError("Your coordinate has to be between the input values.")
else:
limit = 1e-10 # avoid hitting a point directly
diff = numpy.array(xv)-x # get differences
index = numpy.argmin(diff*diff) # nearest neighbour
if (xv[index]<x<xv[index+1] # nearest below, positive inclination
or xv[index]>x>xv[index+1]): # nearest above, negative inclination
if diff[index]<limit:
index = [index-1,index+1]
else:
index = [index, index+1]
elif (xv[index-1]<x<xv[index] # nearest above, positive inclination
or xv[index-1]>x>xv[index]): # nearest below, negative inclination
if diff[index]<limit:
index = [index-1,index+1]
else:
index = [index-1,index]
xvnew = xv[index]
yvnew = yv[index]
return get_x_y_dydx(xvnew,yvnew,x) # Allow for a single recursion
else:
raise ValueError("You have to provide the same amount of x- and y-pairs with at least two entries each.")
if axis is None:
axis=matplotlib.pyplot.gca()
if fig is None:
fig=matplotlib.pyplot.gcf()
if y is None and x is not None:
trash=0
(xv,yv)=ToPixelCoords(xv,yv,axis,fig)
#x is provided but y isn't
(x,trash)=ToPixelCoords(x,trash,axis,fig)
#Get the rotation angle and y-value
x,y,dy_dx = get_x_y_dydx(xv,yv,x)
rot = numpy.arctan(dy_dx)/numpy.pi*180.
elif x is None and y is not None:
#y is provided, but x isn't
_xv = xv[::-1]
_yv = yv[::-1]
#Find x by interpolation
x = interp1d(yv, xv)(y)
trash=0
(xv,yv)=ToPixelCoords(xv,yv,axis,fig)
(x,trash)=ToPixelCoords(x,trash,axis,fig)
#Get the rotation angle and y-value
x,y,dy_dx = get_x_y_dydx(xv,yv,x)
rot = numpy.arctan(dy_dx)/numpy.pi*180.
(x,y)=ToDataCoords(x,y,axis,fig)
return (x,y,rot)
def drawLines(Ref,lines,axis,plt_kwargs=None):
"""
Just an internal method to systematically plot values from
the generated 'line' dicts, method is able to cover the whole
saturation curve. Closes the gap at the critical point and
adds a marker between the two last points of bubble and
dew line if they reach up to critical point.
Returns an array of line objects that can be used to change
the colour or style afterwards.
"""
if not plt_kwargs is None:
for line in lines:
line['opts'] = plt_kwargs
plottedLines = []
if len(lines)==2 and (
'q' in str(lines[0]['type']).lower() and 'q' in str(lines[1]['type']).lower()
) and (
( 0 == lines[0]['value'] and 1 == lines[1]['type'] ) or ( 1 == lines[0]['value'] and 0 == lines[1]['type'] ) ):
# We plot the saturation curve
bubble = lines[0]
dew = lines[1]
line, = axis.plot(bubble['x'],bubble['y'],**bubble['opts'])
plottedLines.extend([line])
line, = axis.plot(dew['x'], dew['y'], **dew['opts'])
plottedLines.extend([line])
# Do we need to test if this is T or p?
Tmax = min(bubble['kmax'],dew['kmax'])
if Tmax>CP.PropsSI(Ref,'Tcrit')-2e-5:
axis.plot(numpy.r_[bubble['x'][-1],dew['x'][-1]],numpy.r_[bubble['y'][-1],dew['y'][-1]],**bubble['opts'])
#axis.plot((bubble['x'][-1]+dew['x'][-1])/2.,(bubble['y'][-1]+dew['y'][-1])/2.,'o',color='Tomato')
else:
for line in lines:
line, = axis.plot(line['x'],line['y'],**line['opts'])
plottedLines.extend([line])
return plottedLines
class IsoLines(BasePlot):
def __init__(self, fluid_ref, graph_type, iso_type, unit_system='SI', **kwargs):
BasePlot.__init__(self, fluid_ref, graph_type, unit_system=unit_system,**kwargs)
if not isinstance(iso_type, str):
raise TypeError("Invalid iso_type input, expected a string")
iso_type = iso_type.upper()
if iso_type not in self.COLOR_MAP.keys() and iso_type != 'Q':
raise ValueError('This kind of isoline is not supported for a ' \
+ str(graph_type) + \
' plot. Please choose from '\
+ str(self.COLOR_MAP.keys()) + ' or Q.')
self.iso_type = iso_type
def __set_axis_limits(self, swap_xy):
class PropertyPlot(BasePlot):
def __init__(self, fluid_name, graph_type, **kwargs):
"""
Generates limits for the axes in terms of x,y defined by 'plot'
based on temperature and pressure.
Create graph for the specified fluid properties
Returns a tuple containing ((xmin, xmax), (ymin, ymax))
Parameters
----------
fluid_name : string or AbstractState
The name of the fluid to be plotted or a state instance
graph_type : string
The graph type to be plotted, like \"PH\" or \"TS\"
axis : :func:`matplotlib.pyplot.gca()`, Optional
The current axis system to be plotted to.
Default: create a new axis system
fig : :func:`matplotlib.pyplot.figure()`, Optional
The current figure to be plotted to.
Default: create a new figure
unit_system : string, ['EUR','KSI','SI']
Select the units used for the plotting. 'EUR' is bar, kJ, C; 'KSI' is kPa, kJ, K; 'SI' is Pa, J, K
tp_limits : string, ['NONE','DEF','ACHP','ORC']
Select the limits in T and p.
reciprocal_density : bool
NOT IMPLEMENTED: If True, 1/rho will be plotted instead of rho
Examples
--------
>>> from CoolProp.Plots import PropertyPlot
>>> plot = PropertyPlot('HEOS::Water', 'TS')
>>> plot.calc_isolines()
>>> plot.show()
>>> import CoolProp
>>> from CoolProp.Plots import PropertyPlot
>>> plot = PropertyPlot('HEOS::R134a', 'PH', unit_system='EUR', tp_limits='ACHP')
>>> plot.calc_isolines(CoolProp.iQ, num=11)
>>> plot.calc_isolines(CoolProp.iT, num=25)
>>> plot.calc_isolines(CoolProp.iS, num=15)
>>> plot.show()
>>> import CoolProp
>>> from CoolProp.Plots import PropertyPlot
>>> plot = PropertyPlot('HEOS::R245fa', 'TS', unit_system='EUR', tp_limits='ORC')
>>> plot.calc_isolines(CoolProp.iQ, num=11)
>>> plot.calc_isolines(CoolProp.iP, iso_range=[1,50], num=10, rounding=True)
>>> plot.draw()
>>> plot.isolines.clear()
>>> plot.props[CoolProp.iP]['color'] = 'green'
>>> plot.props[CoolProp.iP]['lw'] = '0.5'
>>> plot.calc_isolines(CoolProp.iP, iso_range=[1,50], num=10, rounding=False)
>>> plot.show()
.. note::
See the online documentation for a list of the available fluids and
graph types
"""
# Get current axis limits, be sure to set those before drawing isolines
# if no limits are set, use triple point and critical conditions
X = [CP.PropsSI(self.graph_type[1],
'T', 1.5*CP.PropsSI(self.fluid_ref, 'Tcrit'),
'P', CP.PropsSI(self.fluid_ref, 'ptriple'),
self.fluid_ref),
CP.PropsSI(self.graph_type[1],
'T', 1.1*CP.PropsSI(self.fluid_ref, 'Tmin'),
'P', 1.5*CP.PropsSI(self.fluid_ref, 'pcrit'),
self.fluid_ref),
CP.PropsSI(self.graph_type[1],
'T', 1.5*CP.PropsSI(self.fluid_ref, 'Tcrit'),
'P', 1.5*CP.PropsSI(self.fluid_ref, 'pcrit'),
self.fluid_ref),
CP.PropsSI(self.graph_type[1],
'T', 1.1*CP.PropsSI(self.fluid_ref, 'Tmin'),
'P', CP.PropsSI(self.fluid_ref, 'ptriple'),
self.fluid_ref)]
Y = [CP.PropsSI(self.graph_type[0],
'T', 1.5*CP.PropsSI(self.fluid_ref, 'Tcrit'),
'P', CP.PropsSI(self.fluid_ref, 'ptriple'),
self.fluid_ref),
CP.PropsSI(self.graph_type[0],
'T', 1.1*CP.PropsSI(self.fluid_ref, 'Tmin') ,
'P', 1.5*CP.PropsSI(self.fluid_ref, 'pcrit'),
self.fluid_ref),
CP.PropsSI(self.graph_type[0],
'T', 1.1*CP.PropsSI(self.fluid_ref, 'Tcrit'),
'P', 1.5*CP.PropsSI(self.fluid_ref, 'pcrit'),
self.fluid_ref),
CP.PropsSI(self.graph_type[0],
'T', 1.5*CP.PropsSI(self.fluid_ref, 'Tmin') ,
'P', CP.PropsSI(self.fluid_ref, 'ptriple'),
self.fluid_ref)]
limits = [[min(X), max(X)], [min(Y), max(Y)]]
if not self.axis.get_autoscalex_on():
limits[0][0] = max([limits[0][0], min(self.axis.get_xlim())])
limits[0][1] = min([limits[0][1], max(self.axis.get_xlim())])
limits[1][0] = max([limits[1][0], min(self.axis.get_ylim())])
limits[1][1] = min([limits[1][1], max(self.axis.get_ylim())])
# Limits correction in case of KSI unit_system
if self.unit_system == 'KSI':
limits[0] = [l*self.KSI_SCALE_FACTOR[self.graph_type[1]] for l in limits[0]]
limits[1] = [l*self.KSI_SCALE_FACTOR[self.graph_type[0]] for l in limits[1]]
self.axis.set_xlim(limits[0])
self.axis.set_ylim(limits[1])
return limits
def __plotRound(self, values):
super(PropertyPlot, self).__init__(fluid_name, graph_type, **kwargs)
self._isolines = {}
#self._plines = {}
#self._ppoints = {}
self.get_axis_limits()
self._plot_default_annotations()
@property
def isolines(self): return self._isolines
#@property
#def plines(self): return self._plines
#@property
#def ppoints(self): return self._ppoints
def show(self):
self.draw()
super(PropertyPlot, self).show()
def savefig(self, *args, **kwargs):
self.draw()
super(PropertyPlot, self).savefig(*args, **kwargs)
def _plotRound(self, values):
"""
A function round an array-like object while maintaining the
amount of entries. This is needed for the isolines since we
@@ -283,343 +122,249 @@ class IsoLines(BasePlot):
output[i] = numpy.around(inVal[i],decimals=int(val[i]))
output = numpy.unique(output)
return output
def get_isolines(self, iso_range=[], num=None, rounding=False):
def calc_isolines(self, iso_type=None, iso_range=None, num=15, rounding=False, points=250):
"""Calculate lines with constant values of type 'iso_type' in terms of x and y as
defined by the plot object. 'iso_range' either is a collection of values or
simply the minimum and maximum value between which 'num' lines get calculated.
The 'rounding' parameter can be used to generate prettier labels if needed.
"""
This is the core method to obtain lines in the dimensions defined
by 'plot' that describe the behaviour of fluid 'Ref'. The constant
value is determined by 'iName' and has the values of 'iValues'.
'iValues' is an array-like object holding at least one element. Lines
are calculated for every entry in 'iValues'. If the input 'num' is
larger than the amount of entries in 'iValues', an internally defined
pattern is used to calculate an appropriate line spacing between the maximum
and minimum values provided in 'iValues'.
Returns lines[num] - an array of dicts containing 'x' and 'y'
coordinates for bubble and dew line. Additionally, the dict holds
the keys 'label' and 'opts', those can be used for plotting as well.
"""
if iso_range is None or (len(iso_range) == 1 and num != 1):
raise ValueError('Automatic interval detection for isoline \
boundaries is not supported yet, use the \
iso_range=[min, max] parameter.')
if len(iso_range) == 2 and num is None:
raise ValueError('Please specify the number of isoline you want \
e.g. num=10')
if iso_type is None or iso_type == 'all':
for i_type in IsoLine.XY_SWITCH:
if IsoLine.XY_SWITCH[i_type].get(self.y_index*10+self.x_index,None) is not None:
self.calc_isolines(i_type, None, num, rounding, points)
return
if iso_range is None:
if iso_type is CoolProp.iQ:
iso_range = [0.0,1.0]
else:
limits = self.get_axis_limits(iso_type, CoolProp.iT)
iso_range = [limits[0],limits[1]]
if len(iso_range) <= 1 and num != 1:
raise ValueError('You have to provide two values for the iso_range, {0} is not valid.'.format(iso_range))
if len(iso_range) == 2 and (num is None or num < 2):
raise ValueError('Please specify the number of isoline you want e.g. num=10.')
iso_range = numpy.sort(numpy.unique(iso_range))
def generate_ranges(xmin, xmax, num):
if self.iso_type in ['P', 'D']:
return numpy.logspace(math.log(xmin, 2.),
math.log(xmax, 2.),
num=num,
base=2.)
return numpy.linspace(xmin, xmax, num=num)
# Generate iso ranges
if len(iso_range) == 2:
iso_range = generate_ranges(iso_range[0], iso_range[1], num)
#iso_range = plotRound(iso_range)
#else:
# TODO: Automatic interval detection
# iVal = [CP.PropsSI(iName,'T',T_c[i],'D',rho_c[i],Ref) for i in range(len(T_c))]
# iVal = patterns[iName]([numpy.min(iVal),numpy.max(iVal),num])
iso_range = self.generate_ranges(iso_type, iso_range[0], iso_range[1], num)
if rounding:
iso_range = self.__plotRound(iso_range)
switch_xy_map = {'D': ['TS', 'PH', 'PS'],
'S': ['PH', 'PD', 'PT'],
'T': ['PH', 'PS'],
'H': ['PD']}
#TS: TD is defined, SD is not
#PH: PD is defined, HD is not
#PS: PD is defined, SD is not
#PH: PS is more stable than HS
#PD: PS is defined, DS is not
#PT: PS is defined, TS is not
#PH: PT is defined, HT is not
#PS: PT is defined, ST is not
#PD: PH is defined, DH is not
iso_error_map = {'TD': ['S', 'H'],
'HS': ['T', 'D'],}
switch_xy = False
if self.iso_type in ['D', 'S', 'T', 'H']:
if self.graph_type in switch_xy_map[self.iso_type]:
switch_xy = True
if self.graph_type in ['TD', 'HS']:
if self.iso_type in iso_error_map[self.graph_type]:
raise ValueError('You should not reach this point!')
axis_limits = self.__set_axis_limits(switch_xy)
req_prop = self.graph_type[0]
prop2_name = self.graph_type[1]
if switch_xy:
axis_limits.reverse()
req_prop = self.graph_type[1]
prop2_name = self.graph_type[0]
# Calculate the points
if self.iso_type == 'Q':
lines = self._get_sat_lines(x=iso_range)
return lines
# TODO: Determine saturation state if two phase region present
x_range = numpy.linspace(axis_limits[0][0], axis_limits[0][1], 1000.)
x_mesh = [x_range for i in iso_range]
plot_data = self._get_fluid_data(req_prop,
self.iso_type, iso_range,
prop2_name, x_mesh)
if switch_xy:
plot_data = plot_data[::-1]
lines = []
for j in range(len(plot_data[0])):
line = {
'x': plot_data[0][j],
'y': plot_data[1][j],
# TODO
'label': "", #_getIsoLineLabel(self.iso_type, iso_range[j]),
'type': self.iso_type,
'opts': {'color': self.COLOR_MAP[self.iso_type], 'lw':0.75, 'alpha':0.5 }
}
lines.append(line)
return lines
def draw_isolines(self, iso_range, num=None, rounding=False):
"""
Draw lines with constant values of type 'which' in terms of x and y as
defined by 'plot'. 'iMin' and 'iMax' are minimum and maximum value between
which 'num' get drawn.
There should also be helpful error messages...
"""
if iso_range is None or (len(iso_range) == 1 and num != 1):
raise ValueError('Automatic interval detection for isoline \
boundaries is not supported yet, use the \
iso_range=[min, max] parameter.')
if len(iso_range) == 2 and num is None:
raise ValueError('Please specify the number of isoline you want \
e.g. num=10')
if self.iso_type == 'all':
raise ValueError('Plotting all lines automatically is not \
supported, yet..')
if self.iso_type != 'all':
lines = self.get_isolines(iso_range, num, rounding)
drawn_lines = drawLines(self.fluid_ref, lines, self.axis)
self._plot_default_annotations()
return drawn_lines
#else:
# # TODO: assign limits to values automatically
# ll = _getIsoLineIds(plot)
# if not len(ll)==len(iValues):
# raise ValueError('Please provide a properly sized array of bounds.')
# for c,l in enumerate(ll):
# drawIsoLines(Ref, plot, l, iValues=iValues[c], num=num, axis=axis, fig=fig)
class PropsPlot(BasePlot):
def __init__(self, fluid_name, graph_type, units = 'KSI', reciprocal_density = False, **kwargs):
"""
Create graph for the specified fluid properties
iso_range = self._plotRound(iso_range)
# Limits are already in SI units
limits = self._get_axis_limits()
ixrange = self.generate_ranges(self._x_index,limits[0],limits[1],points)
iyrange = self.generate_ranges(self._y_index,limits[2],limits[3],points)
dim = self._system[iso_type]
lines = self.isolines.get(iso_type, [])
for i in range(num):
lines.append(IsoLine(iso_type,self._x_index,self._y_index, value=dim.to_SI(iso_range[i]), state=self._state))
lines[-1].calc_range(ixrange,iyrange)
lines[-1].sanitize_data()
self.isolines[iso_type] = lines
return
def draw_isolines(self):
dimx = self._system[self._x_index]
dimy = self._system[self._y_index]
sat_props = self.props[CoolProp.iQ].copy()
if 'lw' in sat_props: sat_props['lw'] *= 2.0
else: sat_props['lw'] = 1.0
if 'alpha' in sat_props: min([sat_props['alpha']*2.0,1.0])
else: sat_props['alpha'] = 1.0
for i in self.isolines:
props = self.props[i]
dew = None; bub = None
xcrit = None; ycrit = None
if i == CoolProp.iQ:
for line in self.isolines[i]:
if line.value == 0.0: bub = line
elif line.value == 1.0: dew = line
if dew is not None and bub is not None:
xmin,xmax,ymin,ymax = self.get_axis_limits()
xmin = dimx.to_SI(xmin)
xmax = dimx.to_SI(xmax)
ymin = dimy.to_SI(ymin)
ymax = dimy.to_SI(ymax)
dx = xmax-xmin
dy = ymax-ymin
dew_filter = numpy.logical_and(numpy.isfinite(dew.x),numpy.isfinite(dew.y))
#dew_filter = numpy.logical_and(dew_filter,dew.x>dew.x[-1])
stp = min([dew_filter.size,10])
dew_filter[0:-stp] = False
bub_filter = numpy.logical_and(numpy.isfinite(bub.x),numpy.isfinite(bub.y))
if self._x_index == CoolProp.iP or self._x_index == CoolProp.iDmass:
filter_x = lambda x: numpy.log10(x)
else:
filter_x = lambda x: x
if self._y_index == CoolProp.iP or self._y_index == CoolProp.iDmass:
filter_y = lambda y: numpy.log10(y)
else:
filter_y = lambda y: y
if (#(filter_x(dew.x[dew_filter][-1])-filter_x(bub.x[bub_filter][-1])) > 0.010*filter_x(dx) and
(filter_x(dew.x[dew_filter][-1])-filter_x(bub.x[bub_filter][-1])) < 0.050*filter_x(dx) or
(filter_y(dew.y[dew_filter][-1])-filter_y(bub.y[bub_filter][-1])) < 0.010*filter_y(dy)):
f = interp1d(numpy.append(bub.x[bub_filter],dew.x[dew_filter][::-1]),numpy.append(bub.y[bub_filter],dew.y[dew_filter][::-1]),kind='cubic')
x = numpy.linspace(bub.x[bub_filter][-1], dew.x[dew_filter][-1], 11)
y = f(x)
self.axis.plot(dimx.from_SI(x),dimy.from_SI(y),**sat_props)
warnings.warn("Detected an incomplete phase envelope, fixing it numerically.")
xcrit = x[5]; ycrit = y[5]
#Tcrit = self.state.trivial_keyed_output(CoolProp.iT_critical)
#Dcrit = self.state.trivial_keyed_output(CoolProp.irhomass_critical)
#try:
# self.state.update(CoolProp.DmassT_INPUTS, Dcrit, Tcrit)
# xcrit = self.state.keyed_output(self._x_index)
# ycrit = self.state.keyed_output(self._y_index)
#except:
# xcrit = x[5]; ycrit = y[5]
# pass
#self.axis.plot(dimx.from_SI(numpy.array([bub.x[bub_filter][-1], dew.x[dew_filter][-1]])),dimy.from_SI(numpy.array([bub.y[bub_filter][-1], dew.y[dew_filter][-1]])),'o')
for line in self.isolines[i]:
if line.i_index == CoolProp.iQ:
if line.value == 0.0 or line.value == 1.0:
self.axis.plot(dimx.from_SI(line.x),dimy.from_SI(line.y),**sat_props)
else:
if xcrit is not None and ycrit is not None:
self.axis.plot(dimx.from_SI(numpy.append(line.x,xcrit)),dimy.from_SI(numpy.append(line.y,ycrit)),**props)
#try:
# x = numpy.append(line.x,[xcrit])
# y = numpy.append(line.y,[ycrit])
# fltr = numpy.logical_and(numpy.isfinite(x),numpy.isfinite(y))
# f = interp1d(x[fltr][-3:],y[fltr][-3:],kind='linear') # could also be quadratic
# x = numpy.linspace(x[fltr][-2], x[fltr][-1], 5)
# y = f(x)
# #f = interp1d(y[fltr][-5:],x[fltr][-5:],kind='cubic')
# #y = numpy.linspace(y[fltr][-2], y[fltr][-1], 5)
# #x = f(y)
# self.axis.plot(dimx.from_SI(numpy.append(line.x,x)),dimy.from_SI(numpy.append(line.y,y)),**props)
#except:
# self.axis.plot(dimx.from_SI(numpy.append(line.x,xcrit)),dimy.from_SI(numpy.append(line.y,ycrit)),**props)
# pass
else:
self.axis.plot(dimx.from_SI(line.x),dimy.from_SI(line.y),**props)
def draw(self):
self.get_axis_limits()
self.draw_isolines()
#def label_isolines(self, dx=0.075, dy=0.100):
# [xmin, xmax, ymin, ymax] = self.get_axis_limits()
# for i in self.isolines:
# for line in self.isolines[i]:
# if self.get_x_y_dydx(xv, yv, x)
def draw_process(self, statecontainer, points=None, line_opts={'color' : 'r', 'lw' : 1.5}):
""" Draw process or cycle from x and y values in axis units
Parameters
----------
fluid_ref : string
The name of the fluid to be plotted
graph_type : string
The graph type to be plotted
axis : :func:`matplotlib.pyplot.gca()`, Optional
The current axis system to be plotted to.
Default: create a new axis system
fig : :func:`matplotlib.pyplot.figure()`, Optional
The current figure to be plotted to.
Default: create a new figure
units : string, ['KSI','SI']
Select the units used for the plotting. 'KSI' is kPa, kJ, K; 'SI' is Pa, J, K
reciprocal_density : bool
If True, 1/rho will be plotted instead of rho
statecontainer : CoolProp.Plots.SimpleCycles.StateContainer()
A state container object that contains all the information required to draw the process.
Note that points that appear several times get added to a special of highlighted points.
line_opts : dict
Line options (please see :func:`matplotlib.pyplot.plot`), optional
Use this parameter to pass a label for the legend.
Examples
--------
>>> from CoolProp.Plots import PropsPlot
>>> plt = PropsPlot('Water', 'Ph')
>>> plt.show()
>>> plt = PropsPlot('n-Pentane', 'Ts')
>>> plt.set_axis_limits([-0.5, 1.5, 300, 530])
>>> plt.draw_isolines('Q', [0.1, 0.9])
>>> plt.draw_isolines('P', [100, 2000])
>>> plt.draw_isolines('D', [2, 600])
>>> plt.show()
.. note::
See the online documentation for a list of the available fluids and
graph types
"""
BasePlot.__init__(self, fluid_name, graph_type, unit_system=units, **kwargs)
self.smin = kwargs.get('smin', None)
self.smax = kwargs.get('smax', None)
>>> import CoolProp
>>> from CoolProp.Plots import PropertyPlot
>>> pp = PropertyPlot('HEOS::Water', 'TS', unit_system='EUR')
>>> pp.calc_isolines(CoolProp.iP )
>>> pp.calc_isolines(CoolProp.iHmass )
>>> pp.calc_isolines(CoolProp.iQ, num=11)
>>> cycle = SimpleRankineCycle('HEOS::Water', 'TS', unit_system='EUR')
>>> T0 = 300
>>> pp.state.update(CoolProp.QT_INPUTS,0.0,T0+15)
>>> p0 = pp.state.keyed_output(CoolProp.iP)
>>> T2 = 700
>>> pp.state.update(CoolProp.QT_INPUTS,1.0,T2-150)
>>> p2 = pp.state.keyed_output(CoolProp.iP)
>>> cycle.simple_solve(T0, p0, T2, p2, 0.7, 0.8, SI=True)
>>> cycle.steps = 50
>>> sc = cycle.get_state_changes()
>>> pp.draw_process(sc)
>>> # The same calculation can be carried out in another unit system:
>>> cycle.simple_solve(T0-273.15-10, p0/1e5, T2-273.15+50, p2/1e5-5, 0.7, 0.8, SI=False)
>>> sc2 = cycle.get_state_changes()
>>> pp.draw_process(sc2, line_opts={'color':'blue', 'lw':1.5})
>>> pp.show()
"""
warnings.warn("You called the function \"draw_process\", which is not tested.",UserWarning)
dimx = self.system[self.x_index]
dimy = self.system[self.y_index]
marker = line_opts.pop('marker','o')
style = line_opts.pop('linestyle','solid')
style = line_opts.pop('ls',style)
if points is None: points = StateContainer()
xdata = []
ydata = []
old = statecontainer[len(statecontainer)-1]
for i in statecontainer:
point = statecontainer[i]
if point == old:
points.append(point)
old = point
continue
xdata.append(point[self.x_index])
ydata.append(point[self.y_index])
old = point
xdata = dimx.from_SI(numpy.asarray(xdata))
ydata = dimy.from_SI(numpy.asarray(ydata))
self.axis.plot(xdata,ydata,marker='None',linestyle=style,**line_opts)
xdata = numpy.empty(len(points))
ydata = numpy.empty(len(points))
for i in points:
point = points[i]
xdata[i] = point[self.x_index]
ydata[i] = point[self.y_index]
xdata = dimx.from_SI(numpy.asarray(xdata))
ydata = dimy.from_SI(numpy.asarray(ydata))
line_opts['label'] = ''
self.axis.plot(xdata,ydata,marker=marker,linestyle='None',**line_opts)
self._draw_graph()
def __draw_region_lines(self):
lines = self._get_sat_lines(kind='T',
smin=self.smin,
smax=self.smax)
drawLines(self.fluid_ref, lines, self.axis)
def InlineLabel(xv,yv,x=None,y=None,axis=None,fig=None):
warnings.warn("You called the deprecated function \"InlineLabel\", use \"BasePlot.inline_label\".",DeprecationWarning)
plot = PropertyPlot("water","TS",figure=fig,axis=axis)
return plot.inline_label(xv,yv,x,y)
def _draw_graph(self):
self.__draw_region_lines()
self._plot_default_annotations()
def draw_isolines(self, iso_type, iso_range, num=10, rounding=False):
iso_lines = IsoLines(self.fluid_ref,
self.graph_type,
iso_type, unit_system = self.unit_system,
axis=self.axis)
iso_lines.draw_isolines(iso_range, num, rounding)
class PropsPlot(PropertyPlot):
def __init__(self, fluid_name, graph_type, units = 'KSI', reciprocal_density = False, **kwargs):
super(PropsPlot, self).__init__(fluid_name, graph_type, units=units, reciprocal_density=reciprocal_density, **kwargs)
warnings.warn("You called the deprecated class \"PropsPlot\", use \"PropertyPlot\".",DeprecationWarning)
def Ts(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
"""
Deprecated. Use :py:func:`CoolProps.Plots.PropsPlot`
"""
plt = PropsPlot(Ref, 'Ts', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
plt._draw_graph()
if show:
plt.show()
else:
plt._draw_graph()
return plt.axis
def Ph(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
"""
Deprecated. Use :py:func:`CoolProps.Plots.PropsPlot`
"""
plt = PropsPlot(Ref, 'Ph', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
plt._draw_graph()
return plt.axis
def Ps(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
"""
Deprecated. Use :py:func:`CoolProps.Plots.PropsPlot`
"""
plt = PropsPlot(Ref, 'Ps', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
plt._draw_graph()
return plt.axis
def PT(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
"""
Deprecated. Use :py:func:`CoolProps.Plots.PropsPlot`
"""
plt = PropsPlot(Ref, 'PT', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
plt._draw_graph()
return plt.axis
def Prho(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
"""
"""
plt = PropsPlot(Ref, 'PD', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
plt._draw_graph()
return plt.axis
def Trho(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
"""
Deprecated. Use :py:func:`CoolProps.Plots.PropsPlot`
"""
plt = PropsPlot(Ref, 'TD', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
plt._draw_graph()
return plt.axis
def hs(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
"""
Deprecated. Use :py:func:`CoolProps.Plots.PropsPlot`
"""
plt = PropsPlot(Ref, 'hs', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
plt._draw_graph()
return plt.axis
def drawIsoLines(Ref, plot, which, iValues=[], num=0, show=False, axis=None):
"""
Draw lines with constant values of type 'which' in terms of x and y as
defined by 'plot'. 'iMin' and 'iMax' are minimum and maximum value
between which 'num' get drawn.
:Note:
:func:`CoolProps.Plots.drawIsoLines` will be depreciated in future
releases and replaced with :func:`CoolProps.Plots.IsoLines`
Parameters
----------
Ref : str
The given reference fluid
plot : str
The plot type used
which : str
The iso line type
iValues : list
The list of constant iso line values
num : int, Optional
The number of iso lines
(Default: 0 - Use iValues list only)
show : bool, Optional
Show the current plot
(Default: False)
axis : :func:`matplotlib.pyplot.gca()`, Optional
The current axis system to be plotted to.
(Default: create a new axis system)
Examples
--------
>>> from matplotlib import pyplot
>>> from CoolProp.Plots import Ts, drawIsoLines
>>>
>>> Ref = 'n-Pentane'
>>> ax = Ts(Ref)
>>> ax.set_xlim([-0.5, 1.5])
>>> ax.set_ylim([300, 530])
>>> quality = drawIsoLines(Ref, 'Ts', 'Q', [0.3, 0.5, 0.7, 0.8], axis=ax)
>>> isobars = drawIsoLines(Ref, 'Ts', 'P', [100, 2000], num=5, axis=ax)
>>> isochores = drawIsoLines(Ref, 'Ts', 'D', [2, 600], num=7, axis=ax)
>>> pyplot.show()
"""
isolines = IsoLines(Ref, plot, which, axis=axis)
lines = isolines.draw_isolines(iValues, num)
if show:
isolines.show()
return lines
if __name__ == "__main__":
plot = PropertyPlot('HEOS::n-Pentane', 'PD', unit_system='EUR')#, reciprocal_density=True)
plot.calc_isolines(CoolProp.iT)
plot.calc_isolines(CoolProp.iQ, num=11)
#plot.calc_isolines(CoolProp.iSmass)
#plot.calc_isolines(CoolProp.iHmass)
plot.show()

598
wrappers/Python/CoolProp/Plots/SimpleCycles.py
View File

@@ -1,7 +1,19 @@
# -*- coding: utf-8 -*-
from __future__ import print_function, division
from six import with_metaclass
import matplotlib,numpy
import numpy as np
import CoolProp
from CoolProp.CoolProp import PropsSI
from scipy.optimize import newton
from .Common import BasePlot, process_fluid_state, PropertyDict, SIunits
import warnings
from abc import ABCMeta
def SimpleCycle(Ref,Te,Tc,DTsh,DTsc,eta_a,Ts_Ph='Ph',skipPlot=False,axis=None):
"""
@@ -23,6 +35,9 @@ def SimpleCycle(Ref,Te,Tc,DTsh,DTsc,eta_a,Ts_Ph='Ph',skipPlot=False,axis=None):
* skipPlot : If True, won't actually plot anything, just print COP
"""
warnings.warn("This function has been deprecated. PLease consider converting it to an object inheriting from \"BaseCycle\".",DeprecationWarning)
T=numpy.zeros((6))
h=numpy.zeros_like(T)
p=numpy.zeros_like(T)
@@ -107,6 +122,8 @@ def TwoStage(Ref,Q,Te,Tc,DTsh,DTsc,eta_oi,f_p,Tsat_ic,DTsh_ic,Ts_Ph='Ph',prints=
* skipPlot : If True, won't actually plot anything, just print COP
"""
warnings.warn("This function has been deprecated. PLease consider converting it to an object inheriting from \"BaseCycle\".",DeprecationWarning)
T=numpy.zeros((8))
h=numpy.zeros_like(T)
@@ -244,6 +261,8 @@ def EconomizedCycle(Ref,Qin,Te,Tc,DTsh,DTsc,eta_oi,f_p,Ti,Ts_Ph='Ts',skipPlot=Fa
* skipPlot : If True, won't actually plot anything, just print COP
"""
warnings.warn("This function has been deprecated. PLease consider converting it to an object inheriting from \"BaseCycle\".",DeprecationWarning)
m=1
@@ -389,34 +408,563 @@ def EconomizedCycle(Ref,Qin,Te,Tc,DTsh,DTsc,eta_oi,f_p,Ti,Ts_Ph='Ts',skipPlot=Fa
print('Vdisp2: ',(mdot+mdot_inj)/(rho[1]*f*eta_v)*1e6,'cm^3')
return COP
if __name__=='__main__':
from CoolProp.Plots import Ph,Ts
Ref='R290'
fig=matplotlib.pyplot.figure(figsize=(4,3))
ax=fig.add_axes((0.15,0.15,0.8,0.8))
Ph(Ref,Tmin=273.15-30,hbounds=[0,600],axis=ax)
COP=TwoStage('Propane',10000,273.15-5,273.15+43.3,5,7,0.7,0.3,15+273.15,3,prints = True)
matplotlib.pyplot.show()
#class SimpleCycle(object):
# """A class that calculates a simple thermodynamic cycle"""
# def __init__(self, *args, **kwargs):
# object.__init__(self, *args, **kwargs)
# (states, steps, fluid):
# Parameters
# ----------
# x_type : int, str
# Either a letter or an integer that specifies the property type for the x-axis
# y_type : int, str
# Either a letter or an integer that specifies the property type for the y-axis
# states : list
# A collection of state points that follows a fixed scheme defined
# in the implementing subclass.
# fluid_ref : str, CoolProp.AbstractState
# The fluid property provider, either a subclass of CoolProp.AbstractState
# or a string that can be used to generate a CoolProp.AbstractState instance
# via :func:`Common.process_fluid_state`.
# steps : int
# The number of steps used for going from one state to another
#
# for more properties, see :class:`CoolProp.Plots.Common.Base2DObject`.
Ref='R290'
fig=matplotlib.pyplot.figure(figsize=(4,3))
ax=fig.add_axes((0.15,0.15,0.8,0.8))
Ph(Ref,Tmin=273.15-30,hbounds=[0,600],axis=ax)
COP=SimpleCycle(Ref,273.15-5,273.15+45,5,7,0.7,Ts_Ph='Ph')
matplotlib.pyplot.show()
# # See http://stackoverflow.com/questions/1061283/lt-instead-of-cmp
# class ComparableMixin:
# """A mixin class that implements all comparing mathods except for __lt__"""
# def __eq__(self, other):
# return not self<other and not other<self
# def __ne__(self, other):
# return self<other or other<self
# def __gt__(self, other):
# return other<self
# def __ge__(self, other):
# return not self<other
# def __le__(self, other):
# return not other<self
Ref='R410A'
fig=matplotlib.pyplot.figure(figsize=(4,3))
ax=fig.add_axes((0.15,0.15,0.8,0.8))
Ts(Ref,Tmin=273.15-100,sbounds=[0,600],axis=ax)
COP=SimpleCycle(Ref,273.15-5,273.15+45,5,7,0.7,Ts_Ph='Ts')
matplotlib.pyplot.show()
class StatePoint(PropertyDict):
"""A simple fixed dimension dict represented by an object with attributes"""
# Significant digits in SI units
ROUND_DECIMALS = {
CoolProp.iDmass : 5,
CoolProp.iHmass : 5,
CoolProp.iP : 2,
CoolProp.iSmass : 5,
CoolProp.iT : 5,
CoolProp.iUmass : 5,
CoolProp.iQ : 5
}
def __iter__(self):
"""Make sure we always iterate in the same order"""
keys = [CoolProp.iDmass,CoolProp.iHmass,CoolProp.iP,CoolProp.iSmass,CoolProp.iT]
for key in sorted(keys):
yield key
def __str__(self):
return str(self.__dict__)
def __prop_compare(self,other,typ):
# TODO
if self[typ] is None and other[typ] is None: return 0
elif self[typ] is None and other[typ] is not None: return -1
elif self[typ] is not None and other[typ] is None: return 1
else:
A = np.round(self[typ] , self.ROUND_DECIMALS[typ])
B = np.round(other[typ], self.ROUND_DECIMALS[typ])
if A>B: return 1
elif A<B: return -1
elif A==B: return 0
else: raise ValueError("Comparison failed.")
def __eq__(self, other):
for i in self:
if not self.__prop_compare(other,i) == 0:
return False
return True
def __hash__(self):
return hash(repr(self))
class StateContainer(object):
"""A collection of values for the main properties, built to mixin with :class:`CoolProp.Plots.Common.PropertyDict`
Examples
--------
This container has overloaded accessor methods. Just pick your own flavour
or mix the styles as you like:
>>> from __future__ import print_function
>>> import CoolProp
>>> from CoolProp.PLots.SimpleCycles import StateContainer
>>> T0 = 300.000; p0 = 200000.000; h0 = 112745.749; s0 = 393.035
>>> cycle_states = StateContainer()
>>> cycle_states[0,'H'] = h0
>>> cycle_states[0]['S'] = s0
>>> cycle_states[0][CoolProp.iP] = p0
>>> cycle_states[0,CoolProp.iT] = T0
>>> cycle_states[1,"T"] = 300.064
>>> print(cycle_states)
Stored State Points:
state T (K) p (Pa) ρ (kg/m³) h (J/kg) s (J/kg/K)
0 300.000 200000.000 996.601 112745.749 393.035
1 300.064 - - - -
"""
def __init__(self,unit_system=SIunits()):
self._points = {}
self._units = unit_system
@property
def points(self): return self._points
@points.setter
def points(self, value): self._points = value
@property
def units(self): return self._units
@units.setter
def units(self, value): self._units = value
def get_point(self, index, SI=True):
if SI:
state = self[index]
else:
state = self[index]
for i in state:
state[i] = self.units[i].from_SI(state[i])
return state
def set_point(self, index, value, SI=True):
if SI:
self._points[index] = value
else:
for i in value:
self._points[index][i] = self.units[i].to_SI(value[i])
def _list_like(self, value):
"""Try to detect a list-like structure excluding strings"""
return (not hasattr(value, "strip") and
(hasattr(value, "__getitem__") or
hasattr(value, "__iter__")))
# return is_sequence(value) # use from pandas.core.common import is_sequence
def __len__(self):
"""Some cheating to get the correct behaviour"""
return len(self._points)
def __iter__(self):
"""Make sure we iterate in the righ order"""
for key in sorted(self._points):
yield key
def __getitem__(self, index):
"""Another tweak that changes the default access path"""
if self._list_like(index):
if len(index)==0: raise IndexError("Received empty index.")
elif len(index)==1: return self._points[index[0]]
elif len(index)==2: return self._points[index[0]][index[1]]
else: raise IndexError("Received too long index.")
return self._points[index]
def __setitem__(self, index, value):
"""Another tweak that changes the default access path"""
if self._list_like(index):
if len(index)==0: raise IndexError("Received empty index.")
elif len(index)==1: self._points[index[0]] = value
elif len(index)==2:
# safeguard against empty entries
if index[0] not in self._points:
self._points[index[0]] = StatePoint()
self._points[index[0]][index[1]] = value
else: raise IndexError("Received too long index.")
else:
self._points[index] = value
def __str__(self):
out = "Stored State Points:\n"
keys = True
for i in self._points:
if keys:
row = ["{0:>5s}".format("state")]
for j in self._points[i]:
label = u"{0:s} ({1:s})".format(self.units[j].symbol,self.units[j].unit)
row.append(u"{0:>11s}".format(label))
out = out + " ".join(row) + "\n"
keys = False
row = ["{0:>5s}".format(str(i))]
for j in self._points[i]:
try:
row.append(u"{0:11.3f}".format(self.units[j].from_SI(self._points[i][j])))
except:
row.append(u"{0:>11s}".format("-"))
out = out + " ".join(row) + "\n"
return out.encode('utf8', 'replace')
def append(self,new):
i = 0 + self.__len__()
for j in new:
self[i,j] = new[j]
return self
def extend(self,new):
i = 0 + self.__len__()
for j in new:
for k in new[j]:
self[i,k] = new[j][k]
i = i +1
return self
@property
def D(self): return np.array([self._points[k].D for k in self])
@property
def H(self): return np.array([self._points[k].H for k in self])
@property
def P(self): return np.array([self._points[k].P for k in self])
@property
def S(self): return np.array([self._points[k].S for k in self])
@property
def T(self): return np.array([self._points[k].T for k in self])
@property
def U(self): return np.array([self._points[k].U for k in self])
@property
def Q(self): return np.array([self._points[k].Q for k in self])
class BaseCycle(BasePlot):
"""A simple thermodynamic cycle, should not be used on its own."""
# Define the iteration keys
PROPERTIES = {
CoolProp.iDmass : 'density',
CoolProp.iHmass : 'specific enthalpy',
CoolProp.iP : 'pressure',
CoolProp.iSmass : 'specific entropy',
CoolProp.iT : 'temperature'
}
STATECOUNT=0
"""A list of accepted numbers of states"""
STATECHANGE=None
"""A list of lists of tuples that defines how the state transitions
behave for the corresponding entry in BaseCycle.STATECOUNT"""
def __init__(self, fluid_ref, graph_type, unit_system='EUR', **kwargs):
"""Initialises a simple cycle calculator
Parameters
----------
fluid_ref : str, CoolProp.AbstractState
The fluid property provider, either a subclass of CoolProp.AbstractState
or a string that can be used to generate a CoolProp.AbstractState instance
via :func:`Common.process_fluid_state`.
graph_type : string
The graph type to be plotted, like \"PH\" or \"TS\"
unit_system : string, ['EUR','KSI','SI']
Select the units used for the plotting. 'EUR' is bar, kJ, C; 'KSI' is kPa, kJ, K; 'SI' is Pa, J, K
for more properties, see :class:`CoolProp.Plots.Common.BasePlot`.
"""
self._cycle_states = StateContainer()
self._steps = 2
BasePlot.__init__(self, fluid_ref, graph_type, unit_system, **kwargs)
@property
def cycle_states(self): return self._cycle_states
@cycle_states.setter
def cycle_states(self, value):
if len(value) != self.STATECOUNT:
raise ValueError("Your number of states ({0:d}) is not in the list of allowed state counts: {1:s}.".format(len(value),str(self.STATECOUNT)))
self._cycle_states = value
@property
def steps(self): return self._steps
@steps.setter
def steps(self, value): self._steps = int(max([value,2]))
@BasePlot.system.setter
def system(self, value):
if value in self.UNIT_SYSTEMS:
self._system = self.UNIT_SYSTEMS[value]
elif isinstance(value, PropertyDict):
self._system = value
else:
raise ValueError("Invalid unit_system input \"{0:s}\", expected a string from {1:s}".format(str(value),str(self.UNIT_SYSTEMS.keys())))
self._cycle_states._units = self._system
def valid_states(self):
"""Check the formats of BaseCycle.STATECOUNT and BaseCycle.STATECHANGE"""
if len(self.STATECHANGE) != self.STATECOUNT:
raise ValueError("Invalid number of states and or state change operations")
return True
def fill_states(self,objs=None):
"""Try to populate all fields in the state objects"""
if objs is None:
objs = self._cycle_states
local = True
else:
local = False
for i in objs:
full = True
for j in objs[i]:
if objs[i][j] is None:
full = False
if full: continue
if (objs[i][CoolProp.iDmass] is not None and
objs[i][CoolProp.iT] is not None):
self._state.update(CoolProp.DmassT_INPUTS, objs[i][CoolProp.iDmass], objs[i][CoolProp.iT])
elif (objs[i][CoolProp.iP] is not None and
objs[i][CoolProp.iHmass] is not None):
self._state.update(CoolProp.HmassP_INPUTS, objs[i][CoolProp.iHmass], objs[i][CoolProp.iP])
elif (objs[i][CoolProp.iP] is not None and
objs[i][CoolProp.iSmass] is not None):
self._state.update(CoolProp.PSmass_INPUTS, objs[i][CoolProp.iP], objs[i][CoolProp.iSmass])
else:
warnings.warn("Please fill the state[{0:s}] manually.".format(str(i)))
continue
for j in objs[i]:
if objs[i][j] is None:
objs[i][j] = self._state.keyed_output(j)
if local: self._cycle_states = objs
return objs
def state_change(self,in1,in2,start,ty1='lin',ty2='lin'):
"""Calculates a state change defined by the properties in1 and in2
Uses self.states[start] and self.states[start+1] (or self.states[0]) to define
the process and interpolates between the values.
Parameters
----------
in1 : int
The index of the first defined property.
in2 : int
The index of the second defined property.
start : int
The index of the start state.
ty1 : str
The key that defines the type of state change for in1, lin or log.
ty2 : str
The key that defines the type of state change for in2, lin or log.
Returns
-------
scalar or array_like
a list of the length of self.steps+1 that describes the process. It includes start and end state.
"""
self.fill_states()
end = start + 1
if end >= len(self.cycle_states): end -= len(self.cycle_states)
start = self.cycle_states[start]
end = self.cycle_states[end]
#
val = []
inv = [in1,in2]
typ = [ty1,ty2]
for i,v in enumerate(inv):
if typ[i] == 'lin':
val.append(np.linspace(start[v], end[v], self.steps))
elif typ[i] == 'log':
val.append(np.logspace(np.log10(start[v]), np.log10(end[v]), self.steps))
else:
raise ValueError("Unknow range generator {0:s}".format(str(typ[i])))
sc = StateContainer(self._system)
for i,_ in enumerate(val[0]):
sc[i,inv[0]] = val[0][i]
sc[i,inv[1]] = val[1][i]
return self.fill_states(sc)
def get_state_change(self, index):
return self.STATECHANGE[index](self)
def get_state_changes(self):
sc = self.get_state_change(0)
for i in range(1,self.STATECOUNT):
sc.extend(self.get_state_change(i))
return sc
class BasePowerCycle(BaseCycle):
"""A thermodynamic cycle for power producing processes.
Defines the basic properties and methods to unify access to
power cycle-related quantities.
"""
def __init__(self, fluid_ref='HEOS::Water', graph_type='TS', **kwargs):
"""see :class:`CoolProp.Plots.SimpleCycles.BaseCycle` for details."""
BaseCycle.__init__(self, fluid_ref, graph_type, **kwargs)
def eta_carnot(self):
"""Carnot efficiency
Calculates the Carnot efficiency for the specified process, :math:`\eta_c = 1 - \frac{T_c}{T_h}`.
Returns
-------
float
"""
Tvector = self._cycle_states.T
return 1. - np.min(Tvector) / np.max(Tvector)
def eta_thermal(self):
"""Thermal efficiency
The thermal efficiency for the specified process(es), :math:`\eta_{th} = \frac{\dot{W}_{exp} - \dot{W}_{pum}}{\dot{Q}_{in}}`.
Returns
-------
float
"""
raise NotImplementedError("Implement it in the subclass.")
class SimpleRankineCycle(BasePowerCycle):
"""A simple Rankine cycle *without* regeneration"""
STATECOUNT=4
STATECHANGE=[
lambda inp: BaseCycle.state_change(inp,'S','P',0,ty1='log',ty2='log'), # Pumping process
lambda inp: BaseCycle.state_change(inp,'H','P',1,ty1='lin',ty2='lin'), # Heat addition
lambda inp: BaseCycle.state_change(inp,'H','P',2,ty1='log',ty2='log'), # Expansion
lambda inp: BaseCycle.state_change(inp,'H','P',3,ty1='lin',ty2='lin') # Heat removal
]
def __init__(self, fluid_ref='HEOS::Water', graph_type='TS', **kwargs):
"""see :class:`CoolProp.Plots.SimpleCycles.BasePowerCycle` for details."""
BasePowerCycle.__init__(self, fluid_ref, graph_type, **kwargs)
def simple_solve(self, T0, p0, T2, p2, eta_exp, eta_pum, fluid=None, SI=True):
""""
A simple Rankine cycle calculation
Parameters
----------
T0 : float
The coldest point, before the pump
p0 : float
The lowest pressure, before the pump
T2 : float
The hottest point, before the expander
p2 : float
The highest pressure, before the expander
eta_exp : float
Isentropic expander efficiency
eta_pum : float
Isentropic pump efficiency
Examples
--------
>>> import CoolProp
>>> from CoolProp.Plots import PropertyPlot
>>> pp = PropertyPlot('HEOS::Water', 'TS', unit_system='EUR')
>>> cycle = SimpleRankineCycle('HEOS::Water', 'TS', unit_system='EUR')
>>> T0 = 300
>>> pp.state.update(CoolProp.QT_INPUTS,0.0,T0+15)
>>> p0 = pp.state.keyed_output(CoolProp.iP)
>>> T2 = 700
>>> pp.state.update(CoolProp.QT_INPUTS,1.0,T2-150)
>>> p2 = pp.state.keyed_output(CoolProp.iP)
>>> cycle.simple_solve(T0, p0, T2, p2, 0.7, 0.8, SI=True)
>>> cycle.steps = 50
>>> sc = cycle.get_state_changes()
>>> pp.draw_process(sc)
"""
if fluid is not None: self.state = process_fluid_state(fluid)
if self._state is None:
raise ValueError("You have specify a fluid before you calculate.")
cycle_states = StateContainer(unit_system=self._system)
if not SI:
Tc = self._system[CoolProp.iT].to_SI
pc = self._system[CoolProp.iP].to_SI
T0 = Tc(T0)
p0 = pc(p0)
T2 = Tc(T2)
p2 = pc(p2)
# Subcooled liquid
self.state.update(CoolProp.PT_INPUTS,p0,T0)
h0 = self.state.hmass()
s0 = self.state.smass()
# Just a showcase for the different accessor methods
cycle_states[0,'H'] = h0
cycle_states[0]['S'] = s0
cycle_states[0][CoolProp.iP] = p0
cycle_states[0,CoolProp.iT] = T0
# Pressurised liquid
p1 = p2
self.state.update(CoolProp.PSmass_INPUTS,p1,s0)
h1 = h0 + (self.state.hmass() - h0) / eta_pum
self.state.update(CoolProp.HmassP_INPUTS,h1,p1)
s1 = self.state.smass()
T1 = self.state.T()
cycle_states[1,'H'] = h1
cycle_states[1,'S'] = s1
cycle_states[1,'P'] = p1
cycle_states[1,'T'] = T1
# Evaporated vapour
self.state.update(CoolProp.PT_INPUTS,p2,T2)
h2 = self.state.hmass()
s2 = self.state.smass()
cycle_states[2,'H'] = h2
cycle_states[2,'S'] = s2
cycle_states[2,'P'] = p2
cycle_states[2,'T'] = T2
# Expanded gas
p3 = p0
self.state.update(CoolProp.PSmass_INPUTS,p3,s2)
h3 = h2 - eta_exp * (h2 - self.state.hmass())
self.state.update(CoolProp.HmassP_INPUTS,h3,p3)
s3 = self.state.smass()
T3 = self.state.T()
cycle_states[3,'H'] = h3
cycle_states[3,'S'] = s3
cycle_states[3,'P'] = p3
cycle_states[3,'T'] = T3
w_net = h2 - h3
q_boiler = h2 - h1
eta_th = w_net / q_boiler
self.cycle_states = cycle_states
self.fill_states()
def eta_thermal(self):
"""Thermal efficiency
The thermal efficiency for the specified process(es), :math:`\eta_{th} = \frac{\dot{W}_{exp} - \dot{W}_{pum}}{\dot{Q}_{in}}`.
Returns
-------
float
"""
w_net = self.cycle_states[2].H - self.cycle_states[3].H - (self.cycle_states[1].H - self.cycle_states[0].H)
q_boiler = self.cycle_states[2].H - self.cycle_states[1].H
return w_net / q_boiler
## for x in numpy.linspace(0,1):
## Ref='REFPROP-MIX:R152A[%g]&R32[%g]' %(x,1-x)
## COP=SimpleCycle(273.15+8,273.15+44,5,7,0.7,skipPlot=True,Ts_Ph='Ph')
## matplotlib.pyplot.show()

6
wrappers/Python/CoolProp/Plots/__init__.py
View File

@@ -1,7 +1,9 @@
#Bring some functions into the Plots namespace for code concision
from __future__ import absolute_import
from .Plots import PropsPlot, IsoLines, drawIsoLines
from .Plots import Ph, Ts, Ps, PT, Prho, Trho, hs
from .Plots import PropertyPlot,PropsPlot
from .Common import IsoLine
#from .Plots import PropsPlot, IsoLines, drawIsoLines
#from .Plots import Ph, Ts, Ps, PT, Prho, Trho, hs
from .SimpleCycles import SimpleCycle, TwoStage, EconomizedCycle

2
wrappers/Python/CoolProp/tests/test_Props.py
View File

@@ -2,7 +2,7 @@ import unittest
from CoolProp.CoolProp import PropsSI
import CoolProp
import numpy as np
def test_input_types():
for Fluid in ['Water']:
for Tvals in [0.5*PropsSI(Fluid,'Tmin')+0.5*PropsSI(Fluid,'Tcrit'),