In [1]:
case_name = '3dCavity'
In [2]:
remote_data = True
#data_dir='/gpfs/asrc/home/acimpoeru/cube_JA/CUBE'
data_dir='/gpfs/asrc/home/acimpoeru/3dCavity_Jamil/JA_Re=10000/JA'
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.0032e06 # Reynolds number
reference_area = 1.5 * 1.5 # frontal areas in metres^2
reference_length = 1.5 # metres
U_inf = 0.02144 # 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 [8]:
### Read experimental data of Prasad and Koseff U_function_of_Y
In [9]:
import zutil
import zutil.post as post
reload(zutil)
reload(post)
#from zutil.post import get_csv_data
re_10000_u_of_y = zutil.post.get_csv_data("/gpfs/asrc/home/acimpoeru/3dCavity_Jamil/JA_Re=10000/JA/re_10000_u_of_y.csv",header=True,remote=True,delim=',')
print re_10000_u_of_y.keys()
#print data['y/H'], data['U/Ub']
#y_exp = data['y/H']
#vel_exp = data['U/Ub']
#print data['y']
In [10]:
### Read the experimental data of Prasad and Koseff V_function_of_X
In [11]:
import zutil
import zutil.post as post
reload(zutil)
reload(post)
#from zutil.post import get_csv_data
re_10000_v_of_x = zutil.post.get_csv_data("/gpfs/asrc/home/acimpoeru/3dCavity_Jamil/JA_Re=10000/JA/re_10000_v_of_x.csv",header=True,remote=True,delim=',')
print re_10000_v_of_x.keys()
#print data['y/H'], data['U/Ub']
In [12]:
# Extract the vertical line.......Two cutting planes
In [13]:
from zutil.post import for_each
def plot_velocity_profile(ax,file_root):
sym = PVDReader( FileName=file_root+'.pvd' )
sym.UpdatePipeline()
slice1 = Slice(SliceType = 'Plane',Input=sym)
slice1.SliceType.Normal = [0.0 , 0.0 , 1.0]
slice1.SliceType.Origin = [0.0 , 0.0 , 0.00001]
slice1.UpdatePipeline()
slice2 = Slice(SliceType = "Plane",Input = slice1)
slice2.SliceType.Normal = [1.0 , 0.0 , 0.0]
slice2.SliceType.Origin = [0.0 , 0.0 , 0.0]
slice2.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]))
count+=1
ax.plot(data,pts_array.GetPoints()[:,1]/75, color='k',linestyle='-',linewidth = '5',marker='o',markersize = 20,label='zCFD')
In [38]:
from zutil.post import for_each
import numpy as np
def plot_velocity_profile_1(ax,file_root):
sym = PVDReader( FileName=file_root+'.pvd' )
sym.UpdatePipeline()
slice3 = Slice(SliceType = 'Plane',Input=sym)
slice3.SliceType.Normal = [0.0 , 1.0 , 0.0]
slice3.SliceType.Origin = [0.0 , 0.0 , 0.000001]
slice3.UpdatePipeline()
slice4 = Slice(SliceType = "Plane",Input = slice3)
slice4.SliceType.Normal = [0.0 , 0.0 , 1.0]
slice4.SliceType.Origin = [0.0 , 0.0 , 0.0]
slice4.UpdatePipeline()
point_data_1 = CellDatatoPointData(Input=slice4)
point_data_1.PassCellData = 0
point_data_1.UpdatePipeline()
print point_data_1
point_client_1 = servermanager.Fetch(point_data_1)
for_each(point_client_1,func=plot_array_1,axis=ax)
def plot_array_1(data_array_1,pts_array_1,**kwargs):
ax = kwargs['axis']
data_1 = []
count = 0
v_vel = data_array_1.GetPointData()['V'][:,1]
for p in pts_array_1.GetPoints()[:,0]:
data_1.append(float(v_vel[count]))
count+=1
ax.plot(pts_array_1.GetPoints()[:,0]/75,data_1, color='k',linestyle='-',linewidth = '5',marker = 'o',markersize=20,label='zCFD')
print point_client_1
In [27]:
# calculate the RMS error
import numpy as np
import zutil.post
print re_10000_u_of_y['u/u_b'].shape, pts_
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 Lid Driven Cavity. Re = 10000.Velocity profiles in the symmetry planes (xy,xz)', fontsize=30, 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.2,1.2,-1.1,1.1])
#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_1(ax,get_case_root(case_name,num_procs))
ax.plot(re_10000_u_of_y['u/u_b'],re_10000_u_of_y['y/H'], color='r', label='Exp U_y Prasad and Koseff',linestyle='-',linewidth = '5',marker='o',markersize=20)
ax.plot(re_10000_v_of_x['x/L'],re_10000_v_of_x['v/u_b'], color='r', label='Exp V_x Prasad and Koseff',linestyle='-',linewidth = '5',marker='o',markersize=20)
#ax.plot(les['U'], les['Height'], color='k', label='LES')
ax.legend(loc='upper right', shadow=True)
fig.savefig("Re=10000.pdf")
#from IPython.display import FileLink, display
#display(FileLink('images/ahmed_vel_profile.pdf'))
show()
In [17]:
from zutil.post import residual_plot, get_case_report
residual_plot(get_case_report(case_name))
show()
In [ ]:
if remote_data:
#print 'Disconnecting from remote paraview server connection'
Disconnect()
pass