In [1]:
case_name = 'ahmed'
In [2]:
remote_data = True
data_dir='AHMED/Fluent'
data_host='jappa@login02'
remote_server_auto = True
paraview_cmd='mpiexec ~/apps/Paraview/bin/pvserver -rc --client-host=localhost -sp=11113'
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)
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
reference_area = 0.112032 # frontal areas in metres^2
reference_length = 1.044 # metres
pressure = 101325.0 # Pa
temperature = 294.5 # 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*pressure/density) # m/s
mach = 40.0/speed_of_sound
from IPython.display import HTML
HTML(print_html_parameters(parameters))
Out[7]:
In [8]:
velocity_profile_factor = 1.0/20.0
ins = open( "data/Ahmed_25_Cp.dat", "r" )
cp_array = {'x':[],'y':[],'z':[],'cp':[]}
for line in ins:
if line[0] != '#' and line[0] != '%':
# Split the list
var = line.split()
if len(var) == 4 and var[1] == '0.00':
cp_array['x'].append(float(var[0])/1000.0)
cp_array['y'].append(float(var[1])/1000.0)
cp_array['z'].append(float(var[2])/1000.0)
cp_array['cp'].append(float(var[3]))
ins.close()
ins = open( "data/Ahmed_25_y=0_global.dat", "r" )
profile_array = {'x':[],'y':[],'z':[],'u':[],'delta':[]}
for line in ins:
if line[0] != '#' and line[0] != '%':
# Split the list
var = line.split()
profile_array['x'].append(float(var[0])/1000.0)
profile_array['y'].append(float(var[1])/1000.0)
profile_array['z'].append(float(var[2])/1000.0)
profile_array['u'].append(float(var[3]))
profile_array['delta'].append(float(var[0])/1000.0+(float(var[3])-40.0)/40.0*velocity_profile_factor)
ins.close()
profile_location = []
profile_location.append(profile_array['x'][0])
for x in profile_array['x']:
if x != profile_location[-1]:
profile_location.append(x)
In [9]:
from zutil.post import for_each
def plot_velocity_profile(ax,file_root):
sym = PVDReader( FileName=file_root+'_symmetry.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.0]
clip.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=clip2)
point_data.PassCellData = 0
point_data.UpdatePipeline()
for x in profile_location:
wall_slice = Slice(Input=point_data, SliceType="Plane" )
wall_slice.SliceType.Normal = [1.0,0.0,0.0]
wall_slice.SliceType.Origin = [float(x), 0.0, 0.0]
wall_slice.UpdatePipeline()
slice_client = servermanager.Fetch(wall_slice)
for_each(slice_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(p[0][0])+(float(u_vel[count])-40.0)/40.0*velocity_profile_factor)
count+=1
ax.plot(data, pts_array.GetPoints()[:,2], color='g',marker='s',linestyle='none',fillstyle='none')
from zutil.post import cp_profile
def get_cp_profile(file_root):
wall = PVDReader( FileName=file_root+'_wall.pvd' )
wall.UpdatePipeline()
# Remove floor
clip = Clip( ClipType="Plane", Input=wall )
clip.ClipType.Normal = [0.0, 0.0, 1.0]
clip.ClipType.Origin = [0, 0, 0.049]
clip.UpdatePipeline()
merge = MergeBlocks(Input=clip)
merge.UpdatePipeline()
cp_data = {'x':[],'y':[],'z':[],'cp':[]}
cp_profile(merge,[0.0,1.0,0.0],[0.0, -0.001, 0.0],func=save_data,cp_data=cp_data)
return cp_data
def save_data(data_array,pts_array,**kwargs):
cp_data = kwargs['cp_data']
count = 0
#print data_array.GetPointData().keys()
for cp in data_array.GetPointData()['cp']:
cp_data['cp'].append(float(cp))
cp_data['x'].append(float(pts_array.GetPoints().item((count,0))))
cp_data['y'].append(float(pts_array.GetPoints().item((count,1))))
cp_data['z'].append(float(pts_array.GetPoints().item((count,2))))
count+=1
In [10]:
from zutil.post import get_case_root
zcfd_cp_data = get_cp_profile(get_case_root(case_name,num_procs))
fig = pl.figure(figsize=(30, 10),dpi=100, facecolor='w', edgecolor='k')
fig.suptitle('Ahmed Body Velocity Profiles', fontsize=20, fontweight='bold')
ax = fig.add_subplot(1,1,1)
ax.grid(True)
ax.set_xlabel('x [m]')
ax.set_ylabel('z [m]')
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))
ax.legend(loc='upper right', shadow=True)
fig.savefig("images/ahmed_vel_profile.pdf")
from IPython.display import FileLink, display
display(FileLink('images/ahmed_vel_profile.pdf'))
show()
In [11]:
fig = pl.figure(figsize=(30, 10),dpi=100, facecolor='w', edgecolor='k')
fig.suptitle('Ahmed Body $C_p$ y=0', fontsize=20, fontweight='bold')
ax = fig.add_subplot(1,1,1)
ax.grid(True)
ax.set_xlabel('x [m]')
ax.set_ylabel('$C_p$')
ax.axis([-1.25,0.25,1.0,-1.0])
ax.plot(zcfd_cp_data['x'],zcfd_cp_data['cp'], color='b', label='zCFD',linestyle='-',marker='p')
ax.plot(cp_array['x'],cp_array['cp'], color='r', label='Expt',linestyle='none',marker='s')
ax2 = ax.twinx()
ax2.set_ylabel('z [m]')
ax2.axis([-1.25,0.25,0.0,0.5])
ax2.plot(zcfd_cp_data['x'],zcfd_cp_data['z'], color='grey',linestyle='-',marker='None',markersize=4)
ax.legend(loc='upper right', shadow=True)
fig.savefig("images/ahmed_cp_profile.pdf")
from IPython.display import FileLink, display
display(FileLink('images/ahmed_cp_profile.pdf'))
show()
In [12]:
from zutil.post import calc_force
def calc_force_wall(file_root,ignore_zone,half_model=True,filter=None):
wall = PVDReader( FileName=file_root+'_wall.pvd' )
wall.UpdatePipeline()
clip = Clip( ClipType="Plane", Input=wall )
clip.ClipType.Normal = [0.0, 0.0, 1.0]
clip.ClipType.Origin = [0, 0, 0.049]
clip.UpdatePipeline()
return calc_force(clip,ignore_zone,half_model,filter)
In [13]:
force = []
ignore_zones = []
pf,ff = calc_force_wall(get_case_root(case_name,num_procs),set(ignore_zones),False)
for i in range(0,3):
pf[i] /= reference_area
ff[i] /= reference_area
force.append(pf[i] + ff[i])
C_D = force[0]
C_S = force[1]
C_L = force[2]
print 'alpha = ' + ('%.4f' %alpha)
print 'beta = ' + ('%.4f' %beta)
print 'reference area = ' + ('%.4f' %reference_area)
print 'reference length = ' + ('%.4f' %reference_length)
print 'pressure = ' + ('%.4f' %pressure)
print 'temperature = ' + ('%.4f' %temperature)
print 'gas_constant = ' + ('%.4f' %gas_constant)
print 'density = ' + ('%.4f' %density)
print 'gamma = ' + ('%.4f' %gamma)
print 'speed of sound = ' + ('%.4f' %speed_of_sound)
print 'mach = ' + ('%.4f' %mach)
print 'C_S = ' + ('%.4f' %C_S)
print 'C_D = ' + ('%.4f' %C_D)
print 'C_L = ' + ('%.4f' %C_L)
print 'C_D_pressure = ' + ('%.4f' %pf[0])
print 'C_D_friction = ' + ('%.4f' %ff[0])
In [14]:
from zutil.post import residual_plot, get_case_report
residual_plot(get_case_report(case_name))
show()
In [16]:
if remote_data:
#print 'Disconnecting from remote paraview server connection'
Disconnect()
pass
In [15]: