In [1]:
case_name = 'cube'
In [2]:
remote_data = True
#data_dir='/gpfs/asrc/home/acimpoeru/cube_JA/CUBE'
#data_dir='/gpfs/asrc/home/acimpoeru/Cube/zCFDInterface/zCFDInterface'
data_dir='/gpfs/asrc/home/acimpoeru/cube_fine_mesh'
#data_dir='/gpfs/asrc/home/acimpoeru/cube_kO_sst/Cube_Spalart_OF/cube/VTK'
data_host='acimpoeru@login02'
remote_server_auto = True
paraview_cmd='mpiexec /gpfs/thirdparty/zenotech/home/jappa/apps/Paraview/bin/pvserver -rc --client-host=localhost -sp=11115'
if not remote_server_auto:
paraview_cmd=None
if not remote_data:
data_host='localhost'
paraview_cmd=None
In [3]:
%pylab inline
from paraview.simple import *
paraview.simple._DisableFirstRenderCameraReset()
import pylab as pl
import math
import numpy as np
In [4]:
from zutil.post import pvserver_connect
if remote_data:
pvserver_connect(data_host=data_host,data_dir=data_dir,paraview_cmd=paraview_cmd,paraview_remote_port=11115)
Get control dictionary
In [5]:
from zutil.post import get_case_parameters,print_html_parameters
parameters=get_case_parameters(case_name,data_host=data_host,data_dir=data_dir)
Get status file
In [6]:
from zutil.post import get_status_dict
status=get_status_dict(case_name,data_host=data_host,data_dir=data_dir)
num_procs = str(status['num processor'])
In [7]:
alpha = 0.0 # degrees
beta = 0.0 # degrees
Re = 0.04e06 # Reynolds number
reference_area = 1.0 # frontal areas in metres^2
reference_length = 1.0 # metres
U_inf = 10.0 # reference velocity
pressure = 101325.0 # Pa
temperature = 288.15 # i.e. 15 degrees celcius
gas_constant = 287.0
density = pressure/(gas_constant * temperature) # kg/m^3 - check
gamma = 1.4
speed_of_sound = math.sqrt(gamma * gas_constant * temperature) # m/s
mach = U_inf/speed_of_sound # Mach number of the case
from IPython.display import HTML
HTML(print_html_parameters(parameters))
Out[7]:
In [9]:
import zutil
import zutil.post as post
reload(zutil)
reload(post)
#from zutil.post import get_csv_data
data = zutil.post.get_csv_data("/gpfs/asrc/home/acimpoeru/cube_JA/CUBE/data.csv",header=True,remote=True,delim=',')
#print data.keys()
#print data['y/H'], data['U/Ub']
y_exp = data['y/H']
vel_exp = data['U/Ub']
#print y_exp, vel_exp
In [11]:
import zutil
import zutil.post as post
reload(zutil)
reload(post)
#from zutil.post import get_csv_data
les = zutil.post.get_csv_data("/gpfs/asrc/home/acimpoeru/cube_JA/CUBE/les.csv",header=True,remote=True,delim=',')
print les.keys()
#print data['y/H'], data['U/Ub']
In [13]:
import zutil
import zutil.post as post
reload(zutil)
reload(post)
#from zutil.post import get_csv_data
OF = zutil.post.get_csv_data("/gpfs/asrc/home/acimpoeru/cube_kO_sst/Cube_Spalart_OF/cube/VTK/Open_Foam.csv",header=True,remote=True,delim=',')
print OF.keys()
In [14]:
from zutil.post import for_each
def plot_velocity_profile(ax,file_root,full_path=False):
if full_path:
sym = PVDReader( FileName=file_root )
else:
sym = PVDReader( FileName=file_root+'.pvd' )
sym.UpdatePipeline()
# clip = Clip( ClipType="Plane", Input=sym )
# clip.ClipType.Normal = [0.0, 1.0, 0.0]
# clip.ClipType.Origin = [0, -0.01, 0
# clip.UpdatePipeline()
slice1 = Slice(SliceType = 'Plane',Input=sym)
slice1.SliceType.Normal = [0.0 , 1.0 , 0.0]
slice1.SliceType.Origin = [0.0 , 0.0 , 0.0]
slice1.UpdatePipeline()
slice2 = Slice(SliceType = "Plane",Input = slice1)
slice2.SliceType.Normal = [1.0 , 0.0 , 0.0]
slice2.SliceType.Origin = [0.5 , 0.0 , 0.0]
slice2.UpdatePipeline()
#clip2 = Clip( ClipType="Plane", Input=clip )
#clip2.ClipType.Normal = [0.0, 0.0, -1.0]
#clip2.ClipType.Origin = [0, 0.0, 0.8]
#clip2.UpdatePipeline()
point_data = CellDatatoPointData(Input=slice2)
point_data.PassCellData = 0
point_data.UpdatePipeline()
print point_data
point_client = servermanager.Fetch(point_data)
for_each(point_client,func=plot_array,axis=ax)
def plot_array(data_array,pts_array,**kwargs):
ax = kwargs['axis']
data = []
count = 0
u_vel = data_array.GetPointData()['V'][:,0]
for p in pts_array.GetPoints()[:,0]:
data.append(float(u_vel[count])/10.0)
count+=1
ax.plot(data, pts_array.GetPoints()[:,2], color='g',marker='o',label='zCFD')
In [15]:
from zutil.post import get_case_root
#zcfd_cp_data = get_cp_profile(get_case_root(case_name,num_procs))
fig = pl.figure(figsize=(20, 10),dpi=600, facecolor='w', edgecolor='k')
fig.suptitle('3D wall mounted cube Velocity profiles', fontsize=20, fontweight='bold')
ax = fig.add_subplot(1,1,1)
ax.grid(True)
ax.set_xlabel('U/U_ref')
ax.set_ylabel('z/H')
#ax.axis([-1.6,0.8,0,0.8])
#ax.plot(zcfd_cp_data['x'],zcfd_cp_data['z'], color='grey',linestyle='-',marker='None',markersize=5)
#ax.plot(profile_array['x'],profile_array['z'], color='grey', label='station',linestyle='none',marker='.',markersize=1)
#ax.plot(profile_array['delta'],profile_array['z'], color='b', label='Expt U',marker='o',linestyle='none',fillstyle='none')
plot_velocity_profile(ax,get_case_root(case_name,num_procs))
#plot_velocity_profile(ax,'/gpfs/asrc/home/acimpoeru/cube_kO_sst/Cube_Spalart_OF/cube/VTK/cube_14000.vtk',True)
ax.plot(data['U/Ub'], data['y/H'], color='r', label='exp')
ax.plot(les['U'], les['Height'], color='k', label='LES')
ax.plot(OF['"U:0"']/0.584, OF['"Points:2"'], color='b', label='Open Foam')
ax.legend(loc='upper right', shadow=True)
fig.savefig("cube.pdf")
#from IPython.display import FileLink, display
#display(FileLink('images/ahmed_vel_profile.pdf'))
show()
In [16]:
## Calculate the RMS Error
In [17]:
from zutil.post import residual_plot, get_case_report
residual_plot(get_case_report(case_name))
show()
In [18]:
if remote_data:
#print 'Disconnecting from remote paraview server connection'
Disconnect()
pass