From 29ba6cf86fad95ace32a20853dfbdf1cbc3f57b8 Mon Sep 17 00:00:00 2001 From: jowr Date: Thu, 29 May 2014 20:06:07 +0200 Subject: [PATCH 1/4] Added more plotting routines --- dev/TTSE/TTSE_ranges.py | 562 +++++++++++++++++++++++++++++++--------- 1 file changed, 439 insertions(+), 123 deletions(-) diff --git a/dev/TTSE/TTSE_ranges.py b/dev/TTSE/TTSE_ranges.py index 5836852e..867b7a51 100644 --- a/dev/TTSE/TTSE_ranges.py +++ b/dev/TTSE/TTSE_ranges.py @@ -3,10 +3,13 @@ import CoolProp.CoolProp as CP import matplotlib as mpl import matplotlib.pyplot as plt import matplotlib.colors as colors -import matplotlib.cm as cmx import matplotlib.ticker import numpy as np import random +import scipy.ndimage +import scipy.interpolate +from mpl_toolkits.mplot3d import Axes3D +import matplotlib._pylab_helpers #fig = plt.figure(figsize=(10,5)) #ax1 = fig.add_axes((0.08,0.1,0.32,0.83)) @@ -21,168 +24,481 @@ import random #plt.savefig('TTSE_BICUBIC.pdf') #plt.close() -fluid = 'water' +def fill_nan(A): + ''' + interpolate to fill nan values + ''' + inds = np.arange(A.shape[0]) + good = np.where(np.isfinite(A)) + f = scipy.interpolate.interp1d(inds[good], A[good],bounds_error=False) + B = np.where(np.isfinite(A),A,f(inds)) + return B -PRINT = True + +def get_Z(X_in,Y_in,fluid): + ''' + Just a wrapper to call CoolProp + ''' + return False + + +def fill_Z(X,Y): + ''' + Use 2D arrays X and Y to calculate a Z value + ''' + Z = np.empty_like(X) + for i in range(len(X)): + for j in range(len(X[i])): + Z[i,j] = get_Z(X[i,j],Y[i,j],fluid) + Z[i] = fill_nan(Z[i]) + tr = np.transpose(Z) + for i in range(len(tr)): + tr[i] = fill_nan(tr[i]) + return np.transpose(tr) + + +def plotLineAndProjection(ax,X,Y,Z,draw='CXYZ',color='black'): + alpha = 0.5 + + if 'C' in draw: + ax.plot(X,Y,zs=Z,color=color) + #else: + # ax.plot(X,Y,Z,color=color,alpha=0) + + xlim = ax.get_xlim() + ylim = ax.get_ylim() + zlim = ax.get_zlim() + + if 'X' in draw: + constArr = np.ones_like(X) * xlim[0] + ax.plot(constArr,Y,zs=Z,color=color,alpha=alpha) + if 'Y' in draw: + constArr = np.ones_like(Y) * ylim[1] + ax.plot(X,constArr,zs=Z,color=color,alpha=alpha) + if 'Z' in draw: + constArr = np.ones_like(Z) * zlim[0] + ax.plot(X,Y,zs=constArr,color=color,alpha=alpha) + + +def make3Dlpot(X,Y,Z=None,ax=None,invert='',draw='CXYZ',color='blue'): + ''' + A simple wrapper around the plotting routines + + invert: which axes to invert could be 'X' or 'YZ' or ... + draw: which lines to draw could be 'C' or 'YZ' or ... + C for curve and the rest are projections in the different dimensions + ''' + + if Z == None: Z = fill_Z(X, Y) + + ## Now we have all data and need to resample it for a smoother plot + resFactor = 10 + Xr = scipy.ndimage.zoom(X, resFactor, order=1) + Yr = scipy.ndimage.zoom(Y, resFactor, order=1) + Zr = scipy.ndimage.zoom(Z, resFactor, order=1) + + + ## Make the plot + if ax==None: + fig = plt.figure() + ax = fig.gca(projection='3d') + else: + fig = '' + + if 'X' in invert: ax.invert_xaxis() + if 'Y' in invert: ax.invert_yaxis() + if 'Z' in invert: ax.invert_zaxis() + + ## Reduce data again and call the plotting wrapper + stride=np.round(len(Xr)/3.0) + + for i in range(len(Xr)): + if np.mod(i,stride)==0: + plotLineAndProjection(ax,Xr[i],Yr[i],Zr[i],draw=draw,color=color) + plotLineAndProjection(ax,Xr[-1],Yr[-1],Zr[-1],draw=draw,color=color) + + for i in range(len(Xr[0])): + if np.mod(i,stride)==0: + Xi = [row[i] for row in Xr] + Yi = [row[i] for row in Yr] + Zi = [row[i] for row in Zr] + plotLineAndProjection(ax,Xi,Yi,Zi,draw=draw,color=color) + Xi = [row[-1] for row in Xr] + Yi = [row[-1] for row in Yr] + Zi = [row[-1] for row in Zr] + plotLineAndProjection(ax,Xi,Yi,Zi,draw=draw,color=color) + + #cset = ax.contour(X, Y, Z, cmap=cm.coolwarm) + #ax.clabel(cset, fontsize=9, inline=1) + +# cmap = plt.get_cmap('jet') +# ax.plot_surface( +# X, Y, Z, +# rstride=np.round(resFactor*points/dpi), +# cstride=np.round(resFactor*points/dpi), +# alpha=0.5, +# cmap=cmap, +# linewidth=0, +# #antialiased=False +# ) + +# +# stride=np.round(len(X)/7) +# +# if 'C' in draw: +# ax.plot_wireframe( +# X, Y, Z, +# rstride=stride, +# cstride=stride, +# alpha=0.5, +# color = color, +# #cmap=cmap, +# #linewidth=0, +# #antialiased=False +# ) +# else: +# ax.plot_wireframe( +# X, Y, Z, +# #rstride=stride, +# #cstride=stride, +# #alpha=0.5, +# #cmap=cmap, +# linewidth=0, +# #antialiased=False +# ) + ## In case we need a surface plot + #from matplotlib import cm + #from matplotlib.ticker import LinearLocator, FormatStrFormatter + # + #surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm, + # linewidth=0, antialiased=False) + #ax.set_zlim(-1.01, 1.01) + # + #ax.zaxis.set_major_locator(LinearLocator(10)) + #ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f')) + # + #fig.colorbar(surf, shrink=0.5, aspect=5) + + + ## Alternative plotting solution + #cmap = plt.get_cmap('jet') + # + #if PRINT: + # ax.plot_surface( + # X, Y, Z, + # rstride=100, + # cstride=100, + # alpha=0.5, + # cmap=cmap, + # linewidth=0, + # antialiased=True + # ) + #else: + # #ax.plot_surface( + # ax.plot_wireframe( + # X, Y, Z, + # rstride=100, + # cstride=100, + # #alpha=0.5, + # #cmap=cmap, + # #linewidth=0, + # #antialiased=False + # ) + + + ### Make the plot + #fig = plt.figure() + #ax = fig.gca(projection='3d') + ##cset = ax.contour(X, Y, Z, cmap=cm.coolwarm) + ##ax.clabel(cset, fontsize=9, inline=1) + # + #ax.plot_wireframe( + # X, Y, Z, + # rstride=np.round(resFactor*10), + # cstride=np.round(resFactor*10), + # #alpha=0.5, + # #cmap=cmap, + # #linewidth=0, + # #antialiased=False + # ) + + + ## Plot the two-phase dome and its projections + #ax.plot(X_TP,Y_TP,Z_TP,color='black') + + #xlim = ax.get_xlim() + #ylim = ax.get_ylim() + #zlim = ax.get_zlim() + # + #constArr = np.ones_like(X_TP) * xlim[0] + #ax.plot(constArr,Y_TP,Z_TP,color='black') + #constArr = np.ones_like(Y_TP) * ylim[1] + #ax.plot(X_TP,constArr,Z_TP,color='black') + #constArr = np.ones_like(Z_TP) * zlim[0] + #ax.plot(X_TP,Y_TP,constArr,color='black') + + ## Or do we want contour lines? + ##cset = ax.contour(X, Y, Z, zdir='z', offset=zlim[0], cmap=cmap) + ##cset = ax.contour(X, Y, Z, zdir='x', offset=xlim[0], cmap=cmap) + ##cset = ax.contour(X, Y, Z, zdir='y', offset=ylim[1], cmap=cmap) + # + #ax.set_xlim(xlim) + #ax.set_ylim(ylim) + #ax.set_zlim(zlim) + + #majorFormatter = matplotlib.ticker.EngFormatter("J/kg") + #majorFormatter = matplotlib.ticker. + majorFormatter = matplotlib.ticker.ScalarFormatter() + majorFormatter.set_scientific(True) + majorFormatter.set_powerlimits((2,2)) + ax.xaxis.set_major_formatter(majorFormatter) + ax.yaxis.set_major_formatter(majorFormatter) + ax.zaxis.set_major_formatter(majorFormatter) + + return fig, ax, Z + + + + + + +fluid = 'water' +CP.enable_TTSE_LUT(fluid) + +PRINT = False if PRINT: - points = 500 - rstride = 1 - cstride = 1 - antiAli = True + points = 200 + dpi = 300 toFile = True else: - points = 200 - rstride = 5 - cstride = 5 - antiAli = False + points = 100 + dpi = 75 toFile = False -Tmin = CP.PropsSI('Tmin','T',0,'P',0,fluid) +Tmin = CP.PropsSI('Tmin' ,'T',0,'P',0,fluid) + 5.0 Tcri = CP.PropsSI('Tcrit','T',0,'P',0,fluid) Tmax = CP.PropsSI('Tcrit','T',0,'P',0,fluid) * 2.0 ## Get phase boundary -Y_TP = np.linspace(Tmin, Tcri-0.5, 0.5*points) +T_TP = np.linspace(Tmin, Tcri-0.5, 0.5*points) +D_TP = CP.PropsSI('D','T',T_TP,'Q',0,fluid) +P_TP = CP.PropsSI('P','T',T_TP,'Q',0,fluid) +H_TP = CP.PropsSI('H','T',T_TP,'Q',0,fluid) +S_TP = CP.PropsSI('S','T',T_TP,'Q',0,fluid) -X_TP = CP.PropsSI('D','T',Y_TP,'Q',0,fluid) -Z_TP = CP.PropsSI('P','T',Y_TP,'Q',0,fluid) +D_TP = np.append(D_TP, CP.PropsSI('D','T',T_TP[::-1],'Q',1,fluid)) +P_TP = np.append(P_TP, CP.PropsSI('P','T',T_TP[::-1],'Q',1,fluid)) +H_TP = np.append(H_TP, CP.PropsSI('H','T',T_TP[::-1],'Q',1,fluid)) +S_TP = np.append(S_TP, CP.PropsSI('S','T',T_TP[::-1],'Q',1,fluid)) -X_TP = np.append(X_TP, CP.PropsSI('D','T',Y_TP[::-1],'Q',1,fluid)) -Z_TP = np.append(Z_TP, CP.PropsSI('P','T',Y_TP[::-1],'Q',1,fluid)) - -X_TP = np.log10(1.0/X_TP) -Z_TP = np.log10(1.0*Z_TP) - -Y_TP = np.append(Y_TP, Y_TP[::-1]) +logV_TP = np.log10(1.0/D_TP) +logP_TP = np.log10(1.0*P_TP) +T_TP = np.append(T_TP, T_TP[::-1]) -pmin = CP.PropsSI('ptriple','T',0,'P',0,fluid) + 1 +## These are the default TTSE ranges +D_L = CP.PropsSI('D','T',Tmin,'Q',0,fluid) +D_V = CP.PropsSI('D','T',Tmin,'Q',1,fluid) +H_L = CP.PropsSI('H','T',Tmin,'Q',0,fluid) +H_V = CP.PropsSI('H','T',Tmin,'Q',1,fluid) +P_L = CP.PropsSI('P','T',Tmin,'Q',0,fluid) +P_V = CP.PropsSI('P','T',Tmin,'Q',1,fluid) -vmin = np.log10(1.0/CP.PropsSI('D','T',Tmin,'P',pmin,fluid)) -vmax = np.max(X_TP) +Hmin = H_L +Hmax = H_L + (H_V-H_L)*2.0 + +#Pmin = CP.PropsSI('ptriple','T',0,'P',0,fluid) + 1 +Pmin = np.min([P_L,P_V]) +Pmax = CP.PropsSI('pcrit','T',0,'P',0,fluid) * 2.0 # should be p_reduce + +Dmin = CP.PropsSI('D','H',Hmin,'P',Pmax,fluid) +Dmax = CP.PropsSI('D','H',Hmax,'P',Pmin,fluid) +## Set the ranges for the plot +X_TP = H_TP +Y_TP = logP_TP +Z_TP = S_TP -X = np.linspace(vmin,vmax,points) # Volume -Y = np.linspace(Tmin, Tmax,points) # Temperature +Xmin = Hmin +Xmax = Hmax +Ymin = np.log10(Pmin) +Ymax = np.log10(Pmax) + +X = np.linspace(Xmin, Xmax, points) # Enthalpy +Y = np.linspace(Ymin, Ymax, points) # Pressure X, Y = np.meshgrid(X, Y) -#R = np.sqrt(X**2 + Y**2) -#Z = np.sin(R) - - - - -def fillZ(X,Y): - Z = np.empty_like(X) - for i in range(len(X)): - Z[i] = np.log10(CP.PropsSI('P','T',Y[i],'D',1.0/np.power(10.0,X[i]),fluid)) +def get_Z(X_in,Y_in,fluid,out='S'): + ''' + Just a wrapper to call CoolProp + ''' + X = X_in + Y = np.power(10.0,Y_in) + Z = np.NAN + try: + Z = CP.PropsSI(out,'H',X,'P',Y,fluid) + except(ValueError): + print "CoolProp failed, returning NAN" + Z = np.NAN return Z +figHPS, axHPS, Z = make3Dlpot(X,Y,invert='Z',draw='XZ') +## Plot the two-phase dome and its projections +plotLineAndProjection(axHPS,X_TP,Y_TP,Z_TP) + +make3Dlpot(X,Y,Z=Z,ax=axHPS,draw='Y') + +axHPS.set_xlabel(r'$h$') +axHPS.set_ylabel(r'log $p$') +axHPS.set_zlabel(r'$s$') + +## Now we also need thet other variables +get_Z_old = get_Z + +def get_Z(X_in,Y_in,fluid): + return get_Z_old(X_in,Y_in,fluid,out='D') +logVdata = np.log10(1.0/fill_Z(X,Y)) + +def get_Z(X_in,Y_in,fluid): + return get_Z_old(X_in,Y_in,fluid,out='T') +Tdata = fill_Z(X,Y) + +HPSdict = {'H': X, 'P': Y, 'S': Z, 'V': logVdata, 'T': Tdata} -Z = fillZ(X, Y) - -from mpl_toolkits.mplot3d import Axes3D -#fig = plt.figure() -#ax = fig.add_subplot(111, projection='3d') - - -#mpl.rcParams['legend.fontsize'] = 10 +######################################################################### +# Start with the next diagram # -#fig = plt.figure() -##ax = fig.gca(projection='3d') -#ax = Axes3D(fig) -#theta = np.linspace(-4 * np.pi, 4 * np.pi, 100) -#z = np.linspace(-2, 2, 100) -#r = z**2 + 1 -#x = r * np.sin(theta) -#y = r * np.cos(theta) -#ax.plot(x, y, z, label='parametric curve') -#ax.legend() - - -from matplotlib import cm -from matplotlib.ticker import LinearLocator, FormatStrFormatter # -#fig = plt.figure() -#ax = fig.gca(projection='3d') -#X = np.arange(-5, 5, 0.25) -#Y = np.arange(-5, 5, 0.25) -#X, Y = np.meshgrid(X, Y) -#R = np.sqrt(X**2 + Y**2) -#Z = np.sin(R) -#surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm, -# linewidth=0, antialiased=False) -#ax.set_zlim(-1.01, 1.01) # -#ax.zaxis.set_major_locator(LinearLocator(10)) -#ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f')) +######################################################################### + +## Set the ranges for the plot +X_TP = logV_TP +Y_TP = T_TP +Z_TP = logP_TP + +#Xmin = np.log10(1.0/Dmax) +#Xmax = np.log10(1.0/Dmin) +Xmin = np.min(X_TP) +Xmax = np.max(X_TP) +Ymin = np.min(Y_TP) +Ymax = Tmax + +X = np.linspace(Xmin, Xmax, points) # Volume +Y = np.linspace(Ymin, Ymax, points) # Temperature +X, Y = np.meshgrid(X, Y) + +def get_Z(X_in,Y_in,fluid,out='P'): + ''' + Just a wrapper to call CoolProp + ''' + X = 1.0/np.power(10.0,X_in) + Y = Y_in + Z = np.NAN + try: + Z = np.log10(CP.PropsSI(out,'D',X,'T',Y,fluid)) + except(ValueError): + print "CoolProp failed, returning NAN" + Z = np.NAN + return Z + +figVTP, axVTP, Z = make3Dlpot(X,Y,draw='XZ') +## Plot the two-phase dome and its projections +plotLineAndProjection(axVTP,X_TP,Y_TP,Z_TP) + +make3Dlpot(X,Y,Z=Z,ax=axVTP,draw='Y') + +axVTP.set_xlabel(r'log $v$') +axVTP.set_ylabel(r'$T$') +axVTP.set_zlabel(r'log $p$') + +## Now we also need thet other variables +get_Z_old = get_Z + +def get_Z(X_in,Y_in,fluid): + return get_Z_old(X_in,Y_in,fluid,out='H') + +Hdata = fill_Z(X,Y) + +def get_Z(X_in,Y_in,fluid): + return get_Z_old(X_in,Y_in,fluid,out='S') + +Sdata = fill_Z(X,Y) + +VTPdict = {'V': X, 'T': Y, 'P': Z, 'H': Hdata, 'S': Sdata} + + + +######################################################################### +# Now we have all the data and can start mixing the +# different definitions # -#fig.colorbar(surf, shrink=0.5, aspect=5) - - - - -fig = plt.figure() -ax = fig.gca(projection='3d') -#cset = ax.contour(X, Y, Z, cmap=cm.coolwarm) -#ax.clabel(cset, fontsize=9, inline=1) - -ax.plot(X_TP,Y_TP,Z_TP,color='black') - -cmap = plt.get_cmap('jet') - -ax.plot_surface( - X, Y, Z, - rstride=rstride, - cstride=cstride, - alpha=0.5, - cmap=cmap, - linewidth=0, - antialiased=antiAli - ) - -#xlim = ax.get_xlim() -#ylim = ax.get_ylim() -#zlim = ax.get_zlim() # -##cset = ax.contour(X, Y, Z, zdir='z', offset=zlim[0], cmap=cmap) -##cset = ax.contour(X, Y, Z, zdir='x', offset=xlim[0], cmap=cmap) -##cset = ax.contour(X, Y, Z, zdir='y', offset=ylim[1], cmap=cmap) +######################################################################### + +#make3Dlpot(fig=figHPS,HPSdict,Y,invert='Z',draw='XZ') + +make3Dlpot(HPSdict['V'],HPSdict['T'],Z=HPSdict['P'],ax=axVTP,draw='XYZ',color='red') + +make3Dlpot(VTPdict['H'],VTPdict['P'],Z=VTPdict['S'],ax=axHPS,draw='C',color='red') + +#make3Dlpot(X,Y,Z=Z,ax=axVTP,draw='C') + +#plotLineAndProjection(axVTP,HPSdict['V'][0],HPSdict['T'][0],HPSdict['P'][0]) + # -#ax.set_xlim(xlim) -#ax.set_ylim(ylim) -#ax.set_zlim(zlim) +#H_data = HPSgrid[0] +#logP_data = HPSgrid[1] +#logV_data = np.empty_like(H_data) +#T_data = np.empty_like(H_data) +# +# +#for i in range(len(H_data)): +# for j in range(len(H_data[0])): +# logV_data[i,j] = np.log10(1.0/get_Z_old(H_data[i,j],logP_data[i,j],fluid,out='D')) +# T_data[i,j] = get_Z_old(H_data[i,j],logP_data[i,j],fluid,out='T') +# +# +#rstride=np.round(len(H_data)/3.0) +#for i in range(len(H_data)): +# if np.mod(i,rstride)==0: +# #ax.plot(logV_data[i],T_data[i],logP_data[i],color='red') +# plotLineAndProjection(ax,logV_data[i],T_data[i],logP_data[i],color='red') +#plotLineAndProjection(ax,logV_data[-1],T_data[-1],logP_data[-1],color='red') +# +#for i in range(len(H_data[0])): +# if np.mod(i,rstride)==0: +# X = [row[i] for row in logV_data] +# Y = [row[i] for row in T_data] +# Z = [row[i] for row in logP_data] +# plotLineAndProjection(ax,X,Y,Z,color='red') +#X = [row[-1] for row in logV_data] +#Y = [row[-1] for row in T_data] +#Z = [row[-1] for row in logP_data] +#plotLineAndProjection(ax,X,Y,Z,color='red') -ax.set_xlabel(r'log $v$') -ax.set_ylabel(r'$T$') -ax.set_zlabel(r'log $p$') +#figures=[manager.canvas.figure +# for manager in matplotlib._pylab_helpers.Gcf.get_all_fig_managers()] +# +#for i, figure in enumerate(figures): +# figure.savefig('figure%d.png' % i) -#fig = plt.figure() -#ax = fig.gca(projection='3d') -#surf = ax.plot_surface( -# X, Y, Z, -# #rstride=1, -# #cstride=1, -# #cmap=plt.get_cmap('jet'), -# linewidth=0, -# antialiased=False -# ) -#ax.set_zlim(-1.01, 1.01) - -#ax.zaxis.set_major_locator(LinearLocator(10)) -#ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f')) #fig.colorbar(surf, shrink=0.5, aspect=5) if toFile: - plt.savefig('TTSE_RANGES.png', dpi = 300, transparent = True) + figHPS.savefig('phs.png', dpi = dpi, transparent = True) + #figHPS.close() + figVTP.savefig('pvT.png', dpi = dpi, transparent = True) + #figVTP.close() #plt.savefig('TTSE_RANGES.pdf') - plt.close() else: plt.show() + #figHPS.show() + #figVTP.show() From b5f0227392f0511355468fd3204ada35fdd9d29b Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Fri, 30 May 2014 11:35:16 +0200 Subject: [PATCH 2/4] Update README.rst Instructions for how to get VxWorks PPC cross-compiler --- wrappers/Labview/vxWorks/README.rst | 102 +++++++++++++++++++++++++++- 1 file changed, 101 insertions(+), 1 deletion(-) diff --git a/wrappers/Labview/vxWorks/README.rst b/wrappers/Labview/vxWorks/README.rst index 06844bbd..61a848e6 100644 --- a/wrappers/Labview/vxWorks/README.rst +++ b/wrappers/Labview/vxWorks/README.rst @@ -26,4 +26,104 @@ To be determined: What do do about _Dtest, _Nan, _Inf ?? See Also -------- -https://decibel.ni.com/content/docs/DOC-13537 \ No newline at end of file +https://decibel.ni.com/content/docs/DOC-13537 + +Use pre-built binaries from FirstForge +-------------------------------------- + +Instructions from http://firstforge.wpi.edu/sf/wiki/do/viewPage/projects.c--11_toochain/wiki/BinaryInstall + +Debian GNU/Linux (Wheezy/Testing) + +Note: Because of minimal dependencies this may work on other Debian based distributions (e.g. Ubuntu). This is, however, experimental. + +This currently uses GCC 4.8.0 and Binutils 2.22-8 (default in debian testing). + +Add this line to /etc/apt/sources.list to add the repository: + +deb http://debian.repo.frc.s3.amazonaws.com jessie main + +Add the maintainer key for the repository + + sudo wget -O - http://debian.repo.frc.s3.amazonaws.com/rbmj.gpg.key | apt-key add - + +Run the following commands: + + sudo apt-get update + sudo apt-get install gcc-powerpc-wrs-vxwork + + +Building your own version of GCC for VxWorks Target +--------------------------------------------------- + +Instructions from http://firstforge.wpi.edu/sf/wiki/do/viewPage/projects.c--11_toochain/wiki/ManualInstall + +Generic Linux/UNIX + +These are the instructions to build from source manually. Note that some lines begin with $ and some begin with #. Lines that start with a $ can be run as a regular user. Lines that start with # must be run as root (on most distributions, just prefix the command with 'sudo'). + +1: Download all the components. + + $ wget http://ftp.gnu.org/gnu/gcc/gcc-4.8.2/gcc-4.8.2.tar.bz2 + $ wget http://ftp.gnu.org/gnu/binutils/binutils-2.23.1.tar.bz2 + $ wget ftp://ftp.ni.com/pub/devzone/tut/updated_vxworks63gccdist.zip + +2: Set WIND_BASE + + # echo 'export WIND_BASE=/usr/powerpc-wrs-vxworks/wind_base' >> /etc/profile + $ source /etc/profile + +3: Install the WindRiver headers and development resources. + +$ unzip updated_vxworks63gccdist.zip + # mkdir -p /usr/powerpc-wrs-vxworks/wind_base/target + # mkdir -p /usr/powerpc-wrs-vxworks/share/ldscripts + # cp -R gccdist/WindRiver/vxworks-6.3/host /usr/powerpc-wrs-vxworks/wind_base + # cp -R gccdist/WindRiver/vxworks-6.3/target/h/. /usr/powerpc-wrs-vxworks/sys-include + # ln -fsT /usr/powerpc-wrs-vxworks/sys-include/wrn/coreip /usr/powerpc-wrs-vxworks/wind_base/target/h + # sed '/ENTRY(_start)/d' < /usr/powerpc-wrs-vxworks/wind_base/target/h/tool/gnu/ldscripts/link.OUT > /usr/powerpc-wrs-vxworks/share/ldscripts/dkm.ld + +4: extract binutils and gcc, and the dependency libraries + + $ tar -jxf gcc-4.8.2.tar.bz2 + $ tar -jxf binutils-2.23.1.tar.bz2 + $ cd gcc-4.8.0 + $ ./contrib/download_prerequisites + $ cd .. + +5: Build & install binutils + + $ mkdir binutils-build + $ cd binutils-build + $ ../binutils-2.23.1/configure --prefix=/usr --target=powerpc-wrs-vxworks --disable-nls + $ make -j4 + # make install + $ cd .. + +6: Build & install gcc + + $ mkdir gcc-build + $ cd gcc-build + $ ../gcc-4.8.2/configure \ + --prefix=/usr \ + --target=powerpc-wrs-vxworks \ + --with-gnu-as \ + --with-gnu-ld \ + --with-headers \ + --disable-shared \ + --disable-libssp \ + --disable-multilib \ + --with-float=hard \ + --enable-languages=c,c++ \ + --enable-libstdcxx \ + --enable-threads=vxworks \ + --without-gconv \ + --disable-libgomp \ + --disable-nls \ + --disable-libmudflap \ + --with-cpu-PPC603 \ + --disable-symvers \ + CFLAGS_FOR_TARGET='-mstrict-align -mlongcall -g -O2' \ + CXXFLAGS_FOR_TARGET='-mstrict-align -mlongcall -g -O2' + $ make -j4 + # make install From 9ac42beb4c0e4e590d86f751e74eaef634c732ef Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Fri, 30 May 2014 11:36:17 +0200 Subject: [PATCH 3/4] Update README.rst --- wrappers/Labview/vxWorks/README.rst | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/wrappers/Labview/vxWorks/README.rst b/wrappers/Labview/vxWorks/README.rst index 61a848e6..daa27e79 100644 --- a/wrappers/Labview/vxWorks/README.rst +++ b/wrappers/Labview/vxWorks/README.rst @@ -65,45 +65,64 @@ These are the instructions to build from source manually. Note that some lines b 1: Download all the components. $ wget http://ftp.gnu.org/gnu/gcc/gcc-4.8.2/gcc-4.8.2.tar.bz2 + $ wget http://ftp.gnu.org/gnu/binutils/binutils-2.23.1.tar.bz2 + $ wget ftp://ftp.ni.com/pub/devzone/tut/updated_vxworks63gccdist.zip 2: Set WIND_BASE # echo 'export WIND_BASE=/usr/powerpc-wrs-vxworks/wind_base' >> /etc/profile + $ source /etc/profile 3: Install the WindRiver headers and development resources. $ unzip updated_vxworks63gccdist.zip # mkdir -p /usr/powerpc-wrs-vxworks/wind_base/target + # mkdir -p /usr/powerpc-wrs-vxworks/share/ldscripts + # cp -R gccdist/WindRiver/vxworks-6.3/host /usr/powerpc-wrs-vxworks/wind_base + # cp -R gccdist/WindRiver/vxworks-6.3/target/h/. /usr/powerpc-wrs-vxworks/sys-include + # ln -fsT /usr/powerpc-wrs-vxworks/sys-include/wrn/coreip /usr/powerpc-wrs-vxworks/wind_base/target/h + # sed '/ENTRY(_start)/d' < /usr/powerpc-wrs-vxworks/wind_base/target/h/tool/gnu/ldscripts/link.OUT > /usr/powerpc-wrs-vxworks/share/ldscripts/dkm.ld 4: extract binutils and gcc, and the dependency libraries $ tar -jxf gcc-4.8.2.tar.bz2 + $ tar -jxf binutils-2.23.1.tar.bz2 + $ cd gcc-4.8.0 + $ ./contrib/download_prerequisites + $ cd .. 5: Build & install binutils $ mkdir binutils-build + $ cd binutils-build + $ ../binutils-2.23.1/configure --prefix=/usr --target=powerpc-wrs-vxworks --disable-nls + $ make -j4 + # make install + $ cd .. 6: Build & install gcc $ mkdir gcc-build + $ cd gcc-build + $ ../gcc-4.8.2/configure \ --prefix=/usr \ --target=powerpc-wrs-vxworks \ @@ -125,5 +144,7 @@ $ unzip updated_vxworks63gccdist.zip --disable-symvers \ CFLAGS_FOR_TARGET='-mstrict-align -mlongcall -g -O2' \ CXXFLAGS_FOR_TARGET='-mstrict-align -mlongcall -g -O2' + $ make -j4 + # make install From 10aaae464961bb7036856d8ebdc232fdb9c2768b Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Fri, 30 May 2014 12:02:32 +0200 Subject: [PATCH 4/4] Added tr1/memory include for linux compilation --- src/Backends/Helmholtz/ExcessHEFunction.cpp | 28 ++-- src/Backends/Helmholtz/ExcessHEFunction.h | 11 +- .../Helmholtz/HelmholtzEOSMixtureBackend.cpp | 128 +++++++++--------- .../Helmholtz/HelmholtzEOSMixtureBackend.h | 2 +- src/CoolProp.cpp | 3 + src/HumidAirProp.cpp | 4 + src/SpeedTest.cpp | 7 +- src/main.cxx | 49 +++---- 8 files changed, 122 insertions(+), 110 deletions(-) diff --git a/src/Backends/Helmholtz/ExcessHEFunction.cpp b/src/Backends/Helmholtz/ExcessHEFunction.cpp index e442f298..164e4c91 100644 --- a/src/Backends/Helmholtz/ExcessHEFunction.cpp +++ b/src/Backends/Helmholtz/ExcessHEFunction.cpp @@ -1,3 +1,5 @@ +#include + #include "ExcessHEFunction.h" #include "mixture_excess_term_JSON.h" @@ -44,12 +46,12 @@ MixtureExcessHELibrary::MixtureExcessHELibrary() for (unsigned int i = 0; i < CAS1V.size(); ++i) { - // Get the vector of CAS numbers + // Get the vector of CAS numbers std::vector CAS; CAS.resize(2); CAS[0] = CAS1V[i]; CAS[1] = CAS2V[i]; - + // Sort the CAS number vector std::sort(CAS.begin(), CAS.end()); @@ -78,7 +80,7 @@ MixtureExcessHELibrary::MixtureExcessHELibrary() { throw ValueError(); } - + // If not in map, add new entry to map with dictionary if (excess_map.find(CAS) == excess_map.end()) { @@ -93,10 +95,10 @@ MixtureExcessHELibrary::MixtureExcessHELibrary() { // Append dictionary to listing excess_map[CAS].push_back(d); - } + } } } - } + } } void ExcessTerm::construct(const std::vector &components) @@ -105,7 +107,7 @@ void ExcessTerm::construct(const std::vector &components) N = components.size(); F.resize(N, std::vector(N, 0)); - DepartureFunctionMatrix.resize(N, std::vector(N, NULL)); + DepartureFunctionMatrix.resize(N); for (unsigned int i = 0; i < N; ++i) { @@ -162,7 +164,7 @@ double ExcessTerm::alphar(double tau, double delta, const std::vectoralphar(tau,delta); } } @@ -174,7 +176,7 @@ double ExcessTerm::dalphar_dTau(double tau, double delta, const std::vectordalphar_dTau(tau,delta); } } @@ -186,7 +188,7 @@ double ExcessTerm::dalphar_dDelta(double tau, double delta, const std::vectordalphar_dDelta(tau,delta); } } @@ -198,7 +200,7 @@ double ExcessTerm::d2alphar_dDelta2(double tau, double delta, const std::vector< for (unsigned int i = 0; i < N-1; i++) { for (unsigned int j = i + 1; j < N; j++) - { + { summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dDelta2(tau,delta); } } @@ -210,7 +212,7 @@ double ExcessTerm::d2alphar_dTau2(double tau, double delta, const std::vectord2alphar_dTau2(tau,delta); } } @@ -222,7 +224,7 @@ double ExcessTerm::d2alphar_dDelta_dTau(double tau, double delta, const std::vec for (unsigned int i = 0; i < N-1; i++) { for (unsigned int j = i + 1; j < N; j++) - { + { summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dDelta_dTau(tau,delta); } } @@ -280,7 +282,7 @@ GERG2008DepartureFunction::GERG2008DepartureFunction(const std::vector & const std::vector &eta,const std::vector &epsilon,const std::vector &beta, const std::vector &gamma, int Npower) { - + /// Break up into power and gaussian terms { std::vector _n(n.begin(), n.begin()+Npower); diff --git a/src/Backends/Helmholtz/ExcessHEFunction.h b/src/Backends/Helmholtz/ExcessHEFunction.h index df1ceec7..99e68f93 100644 --- a/src/Backends/Helmholtz/ExcessHEFunction.h +++ b/src/Backends/Helmholtz/ExcessHEFunction.h @@ -1,6 +1,7 @@ #ifndef EXCESSHE_FUNCTIONS_H #define EXCESSHE_FUNCTIONS_H +#include #include #include "CoolPropFluid.h" @@ -35,7 +36,7 @@ public: d.add_number("F", F[i]); } d.add_number("Npower", cpjson::get_double(val,"Npower")); - + // Terms for the power d.add_double_vector("n", cpjson::get_double_array(val["n"])); d.add_double_vector("d", cpjson::get_double_array(val["d"])); @@ -59,7 +60,7 @@ public: }; }; -/*! +/*! The abstract base class for departure functions for the excess part of the Helmholtz energy */ class DepartureFunction @@ -67,7 +68,7 @@ class DepartureFunction public: DepartureFunction(){}; virtual ~DepartureFunction(){}; - + /// The excess Helmholtz energy of the binary pair /// Pure-virtual function (must be implemented in derived class virtual double alphar(double tau, double delta) = 0; @@ -86,7 +87,7 @@ public: unsigned int N; std::vector > DepartureFunctionMatrix; std::vector > F; - + ExcessTerm(){}; void construct(const std::vector &components); ~ExcessTerm(); @@ -124,4 +125,4 @@ public: }; } /* namespace CoolProp */ -#endif \ No newline at end of file +#endif diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 007ce30c..89fc7dde 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -5,6 +5,8 @@ * Author: jowr */ +#include + #if defined(_MSC_VER) #define _CRTDBG_MAP_ALLOC #define _CRT_SECURE_NO_WARNINGS @@ -67,7 +69,7 @@ void HelmholtzEOSMixtureBackend::set_components(std::vector comp imposed_phase_index = -1; - // Top-level class can hold copies of the base saturation classes, + // Top-level class can hold copies of the base saturation classes, // saturation classes cannot hold copies of the saturation classes if (generate_SatL_and_SatV) { @@ -76,10 +78,6 @@ void HelmholtzEOSMixtureBackend::set_components(std::vector comp SatV.reset(new HelmholtzEOSMixtureBackend(components, false)); SatV->specify_phase(iphase_gas); } - else - { - SatL = NULL; SatV = NULL; - } } void HelmholtzEOSMixtureBackend::set_mole_fractions(const std::vector &mole_fractions) { @@ -154,7 +152,7 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity_dilute(void) { throw NotImplementedError(format("dilute viscosity not implemented for mixtures")); } - + } long double HelmholtzEOSMixtureBackend::calc_viscosity_background() { @@ -203,8 +201,9 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity(void) { // Get reference fluid name std::string fluid_name = component.transport.viscosity_ecs.reference_fluid; + std::vector names(1, fluid_name); // Get a managed pointer to the reference fluid for ECS - std::tr1::shared_ptr ref_fluid(new HelmholtzEOSMixtureBackend(std::vector(1, fluid_name))); + std::tr1::shared_ptr ref_fluid(new HelmholtzEOSMixtureBackend(names)); // Get the viscosity using ECS return TransportRoutines::viscosity_ECS(*this, *(ref_fluid.get())); } @@ -266,8 +265,9 @@ long double HelmholtzEOSMixtureBackend::calc_conductivity(void) { // Get reference fluid name std::string fluid_name = component.transport.conductivity_ecs.reference_fluid; + std::vector name(1, fluid_name); // Get a managed pointer to the reference fluid for ECS - std::tr1::shared_ptr ref_fluid(new HelmholtzEOSMixtureBackend(std::vector(1, fluid_name))); + std::tr1::shared_ptr ref_fluid(new HelmholtzEOSMixtureBackend(name,false)); // Get the viscosity using ECS return TransportRoutines::conductivity_ECS(*this, *(ref_fluid.get())); } @@ -304,7 +304,7 @@ long double HelmholtzEOSMixtureBackend::calc_conductivity(void) } long double lambda_residual = calc_conductivity_background(); - + // Critical part long double lambda_critical = _HUGE; switch(component.transport.conductivity_critical.type) @@ -412,7 +412,7 @@ void HelmholtzEOSMixtureBackend::mass_to_molar_inputs(long &input_pair, double & case PSmass_INPUTS: ///< Pressure in Pa, Entropy in J/kg/K case PUmass_INPUTS: ///< Pressure in Pa, Internal energy in J/kg case HmassSmass_INPUTS: ///< Enthalpy in J/kg, Entropy in J/kg/K - case SmassUmass_INPUTS: ///< Entropy in J/kg/K, Internal energy in J/kg + case SmassUmass_INPUTS: ///< Entropy in J/kg/K, Internal energy in J/kg case DmassHmass_INPUTS: ///< Mass density in kg/m^3, Enthalpy in J/kg case DmassSmass_INPUTS: ///< Mass density in kg/m^3, Entropy in J/kg/K case DmassUmass_INPUTS: ///< Mass density in kg/m^3, Internal energy in J/kg @@ -447,15 +447,15 @@ void HelmholtzEOSMixtureBackend::mass_to_molar_inputs(long &input_pair, double & } void HelmholtzEOSMixtureBackend::update(long input_pair, double value1, double value2 ) { - clear(); + clear(); - if (is_pure_or_pseudopure == false && mole_fractions.size() == 0) { - throw ValueError("Mole fractions must be set"); + if (is_pure_or_pseudopure == false && mole_fractions.size() == 0) { + throw ValueError("Mole fractions must be set"); } mass_to_molar_inputs(input_pair, value1, value2); - // Set the mole-fraction weighted gas constant for the mixture + // Set the mole-fraction weighted gas constant for the mixture // (or the pure/pseudo-pure fluid) if it hasn't been set yet gas_constant(); @@ -598,7 +598,7 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot { long double p_vap = 0.98*static_cast(_pVanc); long double p_liq = 1.02*static_cast(_pLanc); - + if (value < p_vap){ this->_phase = iphase_gas; _Q = -1000; return; } @@ -636,7 +636,7 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot switch (other) { case iSmolar: - { + { if (value > SatV->calc_smolar()){ this->_phase = iphase_gas; return; } @@ -671,7 +671,7 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot } } } - + // Determine Q based on the input provided if (!is_pure_or_pseudopure){throw ValueError("possibly two-phase inputs not supported for pseudo-pure for now");} @@ -683,7 +683,7 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot SaturationSolvers::saturation_T_pure(&HEOS, _T, options); long double Q; - + if (other == iP) { if (value > 100*DBL_EPSILON + HEOS.SatL->p()){ @@ -721,7 +721,7 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot this->_phase = iphase_twophase; } _Q = Q; - // Load the outputs + // Load the outputs _p = _Q*HEOS.SatV->p() + (1-_Q)*HEOS.SatL->p(); _rhomolar = 1/(_Q/HEOS.SatV->rhomolar() + (1-_Q)/HEOS.SatL->rhomolar()); return; @@ -747,7 +747,7 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot _pVanc = components[0]->ancillaries.pV.evaluate(_T); long double p_vap = 0.98*static_cast(_pVanc); long double p_liq = 1.02*static_cast(_pLanc); - + if (value < p_vap){ this->_phase = iphase_gas; _Q = -1000; return; } @@ -785,7 +785,7 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot switch (other) { case iSmolar: - { + { if (value > SatV->calc_smolar()){ this->_phase = iphase_gas; return; } @@ -820,7 +820,7 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot } } } - + // Determine Q based on the input provided if (!is_pure_or_pseudopure){throw ValueError("possibly two-phase inputs not supported for pseudo-pure for now");} @@ -832,7 +832,7 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot SaturationSolvers::saturation_T_pure(&HEOS, _T, options); long double Q; - + if (other == iP) { if (value > 100*DBL_EPSILON + HEOS.SatL->p()){ @@ -870,7 +870,7 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot this->_phase = iphase_twophase; } _Q = Q; - // Load the outputs + // Load the outputs _p = _Q*HEOS.SatV->p() + (1-_Q)*HEOS.SatL->p(); _rhomolar = 1/(_Q/HEOS.SatV->rhomolar() + (1-_Q)/HEOS.SatL->rhomolar()); return; @@ -962,10 +962,10 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot // this->_phase = iphase_liquid; // } // else if (Q > 1+100*DBL_EPSILON){ -// this->_phase = iphase_gas; +// this->_phase = iphase_gas; // } // else{ -// this->_phase = iphase_twophase; +// this->_phase = iphase_twophase; // } // _Q = Q; // // Load the outputs @@ -1070,12 +1070,12 @@ void get_dtau_ddelta(HelmholtzEOSMixtureBackend *HEOS, long double T, long doubl // dp/ddelta|tau ddelta = rhor*R*T*(1+2*delta*dalphar_dDelta+pow(delta, 2)*d2alphar_dDelta2); // dp/dtau|delta - dtau = dT_dtau*rho*R*(1+delta*dalphar_dDelta-tau*delta*d2alphar_dDelta_dTau); + dtau = dT_dtau*rho*R*(1+delta*dalphar_dDelta-tau*delta*d2alphar_dDelta_dTau); break; } case iHmolar: // dh/dtau|delta - dtau = dT_dtau*R*(-pow(tau,2)*(HEOS->d2alpha0_dTau2()+HEOS->d2alphar_dTau2()) + (1+delta*HEOS->dalphar_dDelta()-tau*delta*HEOS->d2alphar_dDelta_dTau())); + dtau = dT_dtau*R*(-pow(tau,2)*(HEOS->d2alpha0_dTau2()+HEOS->d2alphar_dTau2()) + (1+delta*HEOS->dalphar_dDelta()-tau*delta*HEOS->d2alphar_dDelta_dTau())); // dh/ddelta|tau ddelta = rhor*T*R/rho*(tau*delta*HEOS->d2alphar_dDelta_dTau()+delta*HEOS->dalphar_dDelta()+pow(delta,2)*HEOS->d2alphar_dDelta2()); break; @@ -1102,7 +1102,7 @@ void get_dtau_ddelta(HelmholtzEOSMixtureBackend *HEOS, long double T, long doubl long double HelmholtzEOSMixtureBackend::calc_first_partial_deriv_nocache(long double T, long double rhomolar, int Of, int Wrt, int Constant) { long double dOf_dtau, dOf_ddelta, dWrt_dtau, dWrt_ddelta, dConstant_dtau, dConstant_ddelta; - + get_dtau_ddelta(this, T, rhomolar, Of, dOf_dtau, dOf_ddelta); get_dtau_ddelta(this, T, rhomolar, Wrt, dWrt_dtau, dWrt_ddelta); get_dtau_ddelta(this, T, rhomolar, Constant, dConstant_dtau, dConstant_ddelta); @@ -1119,7 +1119,7 @@ long double HelmholtzEOSMixtureBackend::calc_pressure_nocache(long double T, lon SimpleState reducing = calc_reducing_state_nocache(mole_fractions); long double delta = rhomolar/reducing.rhomolar; long double tau = reducing.T/T; - + // Calculate derivative if needed int nTau = 0, nDelta = 1; long double dalphar_dDelta = calc_alphar_deriv_nocache(nTau, nDelta, mole_fractions, tau, delta); @@ -1143,10 +1143,10 @@ long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long do long double T, value, r, eos, rhomolar; HelmholtzEOSMixtureBackend *HEOS; - solver_resid(HelmholtzEOSMixtureBackend *HEOS, long double T, long double value, int other){ + solver_resid(HelmholtzEOSMixtureBackend *HEOS, long double T, long double value, int other){ this->HEOS = HEOS; this->T = T; this->value = value; this->other = other; }; - double call(double rhomolar){ + double call(double rhomolar){ this->rhomolar = rhomolar; switch(other) { @@ -1159,7 +1159,7 @@ long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long do default: throw ValueError(format("Input not supported")); } - + r = eos-value; return r; }; @@ -1176,7 +1176,7 @@ long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long do switch(other) { - + case iSmolar: { ymelt = calc_smolar_nocache(_T, rhomelt); @@ -1204,7 +1204,7 @@ long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long do default: throw ValueError(); } - + if (is_in_closed_range(ymelt, yc, y)) { long double rhomolar = Brent(resid, rhomelt, rhoc, LDBL_EPSILON, 1e-12, 100, errstring); @@ -1216,7 +1216,7 @@ long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long do return rhomolar; } else - { + { throw ValueError(); } } @@ -1270,7 +1270,7 @@ long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long do long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double p, long double rhomolar_guess) { int phase; - + // Define the residual to be driven to zero class solver_TP_resid : public FuncWrapper1D { @@ -1278,11 +1278,11 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double long double T, p, r, peos, rhomolar, rhor, tau, R_u, delta, dalphar_dDelta; HelmholtzEOSMixtureBackend *HEOS; - solver_TP_resid(HelmholtzEOSMixtureBackend *HEOS, long double T, long double p){ - this->HEOS = HEOS; this->T = T; this->p = p; this->rhor = HEOS->get_reducing().rhomolar; + solver_TP_resid(HelmholtzEOSMixtureBackend *HEOS, long double T, long double p){ + this->HEOS = HEOS; this->T = T; this->p = p; this->rhor = HEOS->get_reducing().rhomolar; this->tau = HEOS->get_reducing().T/T; this->R_u = HEOS->gas_constant(); }; - double call(double rhomolar){ + double call(double rhomolar){ this->rhomolar = rhomolar; delta = rhomolar/rhor; dalphar_dDelta = HEOS->calc_alphar_deriv_nocache(0, 1, HEOS->get_mole_fractions(), tau, delta); @@ -1306,7 +1306,7 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double if (rhomolar_guess < 0) // Not provided { rhomolar_guess = solver_rho_Tp_SRK(T, p, phase); - + if (phase == iphase_gas && rhomolar_guess < 0)// If the guess is bad, probably high temperature, use ideal gas { rhomolar_guess = p/(gas_constant()*T); @@ -1363,7 +1363,7 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp_SRK(long double T, long do long double m_j = 0.480+1.574*accentric_j-0.176*pow(accentric_j, 2); long double a_j = 0.42747*pow(R_u*Tcj,2)/pcj*pow(1+m_j*(1-sqrt(T/Tcj)),2); - + if (i == j){ k_ij = 0; } @@ -1390,12 +1390,12 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp_SRK(long double T, long do long double rhomolar0 = p/(Z0*R_u*T); long double rhomolar1 = p/(Z1*R_u*T); long double rhomolar2 = p/(Z2*R_u*T); - + // Check if only one solution is positive, return the solution if that is the case if (rhomolar0 > 0 && rhomolar1 <= 0 && rhomolar2 <= 0){ return rhomolar0; } if (rhomolar0 <= 0 && rhomolar1 > 0 && rhomolar2 <= 0){ return rhomolar1; } - if (rhomolar0 <= 0 && rhomolar1 <= 0 && rhomolar2 > 0){ return rhomolar2; } - + if (rhomolar0 <= 0 && rhomolar1 <= 0 && rhomolar2 > 0){ return rhomolar2; } + switch(phase) { case iphase_liquid: @@ -1410,7 +1410,7 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp_SRK(long double T, long do } long double HelmholtzEOSMixtureBackend::calc_pressure(void) -{ +{ // Calculate the reducing parameters _delta = _rhomolar/_reducing.rhomolar; _tau = _reducing.T/_T; @@ -1606,7 +1606,7 @@ SimpleState HelmholtzEOSMixtureBackend::calc_reducing_state_nocache(const std::v SimpleState reducing; if (is_pure_or_pseudopure){ reducing = components[0]->pEOS->reduce; - + } else{ reducing.T = Reducing.p->Tr(mole_fractions); @@ -1653,13 +1653,13 @@ long double HelmholtzEOSMixtureBackend::calc_alphar_deriv_nocache(const int nTau else if (nTau == 3 && nDelta == 0){ return components[0]->pEOS->d3alphar_dTau3(tau, delta); } - else + else { throw ValueError(); } } else{ - + std::size_t N = mole_fractions.size(); long double summer = 0; if (nTau == 0 && nDelta == 0){ @@ -1702,13 +1702,13 @@ long double HelmholtzEOSMixtureBackend::calc_alphar_deriv_nocache(const int nTau for (unsigned int i = 0; i < N; ++i){ summer += mole_fractions[i]*components[i]->pEOS->d3alphar_dTau3(tau, delta); } return summer + pExcess.d3alphar_dTau3(tau, delta); }*/ - else + else { throw ValueError(); } } } -long double HelmholtzEOSMixtureBackend::calc_alpha0_deriv_nocache(const int nTau, const int nDelta, const std::vector &mole_fractions, +long double HelmholtzEOSMixtureBackend::calc_alpha0_deriv_nocache(const int nTau, const int nDelta, const std::vector &mole_fractions, const long double &tau, const long double &delta, const long double &Tr, const long double &rhor) { if (is_pure_or_pseudopure) @@ -1743,7 +1743,7 @@ long double HelmholtzEOSMixtureBackend::calc_alpha0_deriv_nocache(const int nTau else if (nTau == 3 && nDelta == 0){ return components[0]->pEOS->d3alpha0_dTau3(tau, delta); } - else + else { throw ValueError(); } @@ -1753,31 +1753,31 @@ long double HelmholtzEOSMixtureBackend::calc_alpha0_deriv_nocache(const int nTau std::size_t N = mole_fractions.size(); long double summer = 0; long double tau_i, delta_i, rho_ci, T_ci; - for (unsigned int i = 0; i < N; ++i){ - rho_ci = components[i]->pEOS->reduce.rhomolar; + for (unsigned int i = 0; i < N; ++i){ + rho_ci = components[i]->pEOS->reduce.rhomolar; T_ci = components[i]->pEOS->reduce.T; tau_i = T_ci*tau/Tr; delta_i = delta*rhor/rho_ci; - if (nTau == 0 && nDelta == 0){ - summer += mole_fractions[i]*(components[i]->pEOS->base0(tau_i, delta_i)+log(mole_fractions[i])); + if (nTau == 0 && nDelta == 0){ + summer += mole_fractions[i]*(components[i]->pEOS->base0(tau_i, delta_i)+log(mole_fractions[i])); } else if (nTau == 0 && nDelta == 1){ - summer += mole_fractions[i]*rhor/rho_ci*components[i]->pEOS->dalpha0_dDelta(tau_i, delta_i); + summer += mole_fractions[i]*rhor/rho_ci*components[i]->pEOS->dalpha0_dDelta(tau_i, delta_i); } else if (nTau == 1 && nDelta == 0){ - summer += mole_fractions[i]*T_ci/Tr*components[i]->pEOS->dalpha0_dTau(tau_i, delta_i); + summer += mole_fractions[i]*T_ci/Tr*components[i]->pEOS->dalpha0_dTau(tau_i, delta_i); } else if (nTau == 0 && nDelta == 2){ - summer += mole_fractions[i]*pow(rhor/rho_ci,2)*components[i]->pEOS->d2alpha0_dDelta2(tau_i, delta_i); + summer += mole_fractions[i]*pow(rhor/rho_ci,2)*components[i]->pEOS->d2alpha0_dDelta2(tau_i, delta_i); } else if (nTau == 1 && nDelta == 1){ - summer += mole_fractions[i]*rhor/rho_ci*T_ci/Tr*components[i]->pEOS->d2alpha0_dDelta_dTau(tau_i, delta_i); + summer += mole_fractions[i]*rhor/rho_ci*T_ci/Tr*components[i]->pEOS->d2alpha0_dDelta_dTau(tau_i, delta_i); } else if (nTau == 2 && nDelta == 0){ - summer += mole_fractions[i]*pow(T_ci/Tr,2)*components[i]->pEOS->d2alpha0_dTau2(tau_i, delta_i); + summer += mole_fractions[i]*pow(T_ci/Tr,2)*components[i]->pEOS->d2alpha0_dTau2(tau_i, delta_i); } - else + else { throw ValueError(); } @@ -1940,7 +1940,7 @@ long double HelmholtzEOSMixtureBackend::mixderiv_d_ndalphardni_dxj__constT_V_xi( { // Gernert 3.118 return mixderiv_d_ndalphardni_dxj__constdelta_tau_xi(i,j) - + mixderiv_ddelta_dxj__constT_V_xi(j)*mixderiv_d_ndalphardni_dDelta(i) + + mixderiv_ddelta_dxj__constT_V_xi(j)*mixderiv_d_ndalphardni_dDelta(i) + mixderiv_dtau_dxj__constT_V_xi(j)*mixderiv_d_ndalphardni_dTau(i); } long double HelmholtzEOSMixtureBackend::mixderiv_ddelta_dxj__constT_V_xi(int j) @@ -2025,7 +2025,7 @@ long double HelmholtzEOSMixtureBackend::mixderiv_d_ndalphardni_dxj__constdelta_t return line1+line2+line3+line4+line5; } long double HelmholtzEOSMixtureBackend::mixderiv_nd2nalphardnidnj__constT_V(int i, int j) -{ +{ double line0 = mixderiv_ndalphar_dni__constT_V_nj(j); // First term from 7.46 double line1 = mixderiv_d_ndalphardni_dDelta(i)*mixderiv_nddeltadni__constT_V_nj(j); double line2 = mixderiv_d_ndalphardni_dTau(i)*mixderiv_ndtaudni__constT_V_nj(j); diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index eba921a1..5a47040e 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -28,7 +28,7 @@ protected: SimpleState _crit; int imposed_phase_index; public: - HelmholtzEOSMixtureBackend(){SatL = NULL; SatV = NULL; imposed_phase_index = -1;}; + HelmholtzEOSMixtureBackend(){imposed_phase_index = -1;}; HelmholtzEOSMixtureBackend(std::vector components, bool generate_SatL_and_SatV = true); HelmholtzEOSMixtureBackend(std::vector &component_names, bool generate_SatL_and_SatV = true); virtual ~HelmholtzEOSMixtureBackend(){}; diff --git a/src/CoolProp.cpp b/src/CoolProp.cpp index ba09aedc..989656dd 100644 --- a/src/CoolProp.cpp +++ b/src/CoolProp.cpp @@ -16,6 +16,8 @@ #endif #endif +#include + #include #include #include @@ -27,6 +29,7 @@ #include "Backends/Helmholtz/Fluids/FluidLibrary.h" #include "Backends/Helmholtz/HelmholtzEOSBackend.h" + namespace CoolProp { diff --git a/src/HumidAirProp.cpp b/src/HumidAirProp.cpp index 16242419..c2f9fad7 100644 --- a/src/HumidAirProp.cpp +++ b/src/HumidAirProp.cpp @@ -2,6 +2,8 @@ #define _CRT_SECURE_NO_WARNINGS #endif +#include + #include "HumidAirProp.h" #include "AbstractState.h" #include "Solvers.h" @@ -9,6 +11,8 @@ #include "Ice.h" #include "CoolProp.h" + + #include #include "math.h" #include "time.h" diff --git a/src/SpeedTest.cpp b/src/SpeedTest.cpp index c37499e1..b6a192c3 100644 --- a/src/SpeedTest.cpp +++ b/src/SpeedTest.cpp @@ -1,3 +1,4 @@ +#include #include "SpeedTest.h" #include "AbstractState.h" #include "DataStructures.h" @@ -10,7 +11,7 @@ void compare_REFPROP_and_CoolProp(std::string fluid, int inputs, double val1, do { time_t t1,t2; double dx = 1/((double)N); - + std::tr1::shared_ptr State(AbstractState::factory("HEOS", fluid)); t1 = clock(); for (std::size_t ii = 0; ii < N; ++ii) @@ -18,7 +19,7 @@ void compare_REFPROP_and_CoolProp(std::string fluid, int inputs, double val1, do State->update(inputs, val1 + ii*d1, val2 + ii*d2); } t2 = clock(); - + double elap = ((double)(t2-t1))/CLOCKS_PER_SEC/((double)N)*1e6; printf("Elapsed time for CoolProp is %g us/call\n",elap); @@ -33,4 +34,4 @@ void compare_REFPROP_and_CoolProp(std::string fluid, int inputs, double val1, do printf("Elapsed time for REFPROP is %g us/call\n",elap); } -} /* namespace CoolProp */ \ No newline at end of file +} /* namespace CoolProp */ diff --git a/src/main.cxx b/src/main.cxx index 389fd7d8..046efb98 100644 --- a/src/main.cxx +++ b/src/main.cxx @@ -17,12 +17,12 @@ using namespace CoolProp; #endif #include "SpeedTest.h" -#include +//#include void generate_melting_curve_data(const char* file_name, const char *fluid_name, double Tmin, double Tmax) { - + FILE *fp; fp = fopen(file_name,"w"); AbstractState *State = AbstractState::factory(std::string("REFPROP"),std::string(fluid_name)); @@ -43,7 +43,7 @@ void generate_melting_curve_data(const char* file_name, const char *fluid_name, } catch(std::exception &e) { - + std::cout << fluid_name << " " << e.what() << std::endl; break; } @@ -51,8 +51,13 @@ void generate_melting_curve_data(const char* file_name, const char *fluid_name, fclose(fp); delete State; } +struct element + { + double d,t,ld; + int l; + }; int main() -{ +{ if (0) { generate_melting_curve_data("Ethylene-I.mlt","ethylene",103.989,110.369); @@ -61,7 +66,7 @@ int main() generate_melting_curve_data("Propylene-II.mlt","propylen",109.6,575); generate_melting_curve_data("ParaHydrogen-I.mlt","parahyd",13.8033,22); generate_melting_curve_data("ParaHydrogen-II.mlt","parahyd",22,2000); - + generate_melting_curve_data("n-Propane.mlt","propane",85.53,2000); generate_melting_curve_data("n-Butane.mlt","butane",134.9,2000); generate_melting_curve_data("n-Pentane.mlt","pentane",143.5,2000); @@ -103,12 +108,12 @@ int main() if (0) { std::vector ss = strsplit(get_global_param_string("FluidsList"),','); - + for (std::vector::iterator it = ss.begin(); it != ss.end(); ++it) { AbstractState *S = AbstractState::factory("HEOS", (*it)); S->update(QT_INPUTS, 0, S->Ttriple()); - std::cout << format("%s %17.15g\n", S->name(), S->p()); + std::cout << format("%s %17.15g\n", S->name().c_str(), S->p()); } } if (1) @@ -137,7 +142,7 @@ int main() double p2 = AS->umolar(); double d2 = AS->rhomolar(); double T2 = AS->delta(); - + double dpdT_constrho2 = (p2-p1)/(T2-T1); AS->update(PT_INPUTS, 101000, 300); @@ -179,13 +184,9 @@ int main() } if (0) { - struct element - { - double d,t,ld; - int l; - }; - double n[] = {0.0125335479355233, 7.8957634722828, -8.7803203303561, 0.31802509345418, -0.26145533859358, -0.0078199751687981, 0.0088089493102134, -0.66856572307965, 0.20433810950965, -6.621260503968699e-005, -0.19232721156002, -0.25709043003438, 0.16074868486251, -0.040092828925807, 3.9343422603254e-007, -7.5941377088144e-006, 0.00056250979351888, -1.5608652257135e-005, 1.1537996422951e-009, 3.6582165144204e-007, -1.3251180074668e-012, -6.2639586912454e-010, -0.10793600908932, 0.017611491008752, 0.22132295167546, -0.40247669763528, 0.58083399985759, 0.0049969146990806, -0.031358700712549, -0.74315929710341, 0.4780732991548, 0.020527940895948, -0.13636435110343, 0.014180634400617, 0.008332650488071301, -0.029052336009585, 0.038615085574206, -0.020393486513704, -0.0016554050063734, 0.0019955571979541, 0.00015870308324157, -1.638856834253e-005, 0.043613615723811, 0.034994005463765, -0.076788197844621, 0.022446277332006, -6.2689710414685e-005, -5.5711118565645e-010, -0.19905718354408, 0.31777497330738, -0.11841182425981 }; - double d[] = {1, 1, 1, 2, 2, 3, 4, 1, 1, 1, 2, 2, 3, 4, 4, 5, 7, 9, 10, 11, 13, 15, 1, 2, 2, 2, 3, 4, 4, 4, 5, 6, 6, 7, 9, 9, 9, 9, 9, 10, 10, 12, 3, 4, 4, 5, 14, 3, 6, 6, 6 }; + + double n[] = {0.0125335479355233, 7.8957634722828, -8.7803203303561, 0.31802509345418, -0.26145533859358, -0.0078199751687981, 0.0088089493102134, -0.66856572307965, 0.20433810950965, -6.621260503968699e-005, -0.19232721156002, -0.25709043003438, 0.16074868486251, -0.040092828925807, 3.9343422603254e-007, -7.5941377088144e-006, 0.00056250979351888, -1.5608652257135e-005, 1.1537996422951e-009, 3.6582165144204e-007, -1.3251180074668e-012, -6.2639586912454e-010, -0.10793600908932, 0.017611491008752, 0.22132295167546, -0.40247669763528, 0.58083399985759, 0.0049969146990806, -0.031358700712549, -0.74315929710341, 0.4780732991548, 0.020527940895948, -0.13636435110343, 0.014180634400617, 0.008332650488071301, -0.029052336009585, 0.038615085574206, -0.020393486513704, -0.0016554050063734, 0.0019955571979541, 0.00015870308324157, -1.638856834253e-005, 0.043613615723811, 0.034994005463765, -0.076788197844621, 0.022446277332006, -6.2689710414685e-005, -5.5711118565645e-010, -0.19905718354408, 0.31777497330738, -0.11841182425981 }; + double d[] = {1, 1, 1, 2, 2, 3, 4, 1, 1, 1, 2, 2, 3, 4, 4, 5, 7, 9, 10, 11, 13, 15, 1, 2, 2, 2, 3, 4, 4, 4, 5, 6, 6, 7, 9, 9, 9, 9, 9, 10, 10, 12, 3, 4, 4, 5, 14, 3, 6, 6, 6 }; double t[] = {-0.5, 0.875, 1, 0.5, 0.75, 0.375, 1, 4, 6, 12, 1, 5, 4, 2, 13, 9, 3, 4, 11, 4, 13, 1, 7, 1, 9, 10, 10, 3, 7, 10, 10, 6, 10, 10, 1, 2, 3, 4, 8, 6, 9, 8, 16, 22, 23, 23, 10, 50, 44, 46, 50 }; double l[] = {0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 6, 6, 6, 6 }; double summer = 0; @@ -216,14 +217,14 @@ int main() summer += (di-lid*pow_delta_li)*exp(ti*log_tau+(di-1)*log_delta-pow_delta_li); } else{ - summer += di*exp(ti*log_tau+(di-1)*log_delta); + summer += di*exp(ti*log_tau+(di-1)*log_delta); } } } double t2 = clock(); double elap = (t2-t1)/CLOCKS_PER_SEC/((double)N)*1e6; printf("%g %g\n",elap, summer); - + } if (0) { @@ -240,7 +241,7 @@ int main() { double T = 300; - + AbstractState *MixRP = AbstractState::factory(std::string("REFPROP"), std::string("propane")); { long N = 100000; @@ -324,9 +325,9 @@ int main() AbstractState *Mix = AbstractState::factory(std::string("CORE"),std::string("Ethane,n-Propane")); Mix->set_mole_fractions(z); - + for (double T = 210; ;T += 0.1) - { + { Mix->update(QT_INPUTS, Q, T); std::cout << format(" %g %g\n",Mix->p(),Mix->T()); } @@ -335,7 +336,7 @@ int main() if(0) { time_t t1,t2; - + std::size_t N = 1000000; AbstractState *State = AbstractState::factory(std::string("CORE"), std::string("Water")); double p = State->p(); @@ -363,7 +364,7 @@ int main() } - + if (0) { @@ -414,5 +415,5 @@ int main() //double sigma = State->surface_tension(); delete State; - } -} \ No newline at end of file + } +}