3D Lid Driven Cavity. Validation Test Case

References

Prasad and Koseff (1989). Reynolds number effect and end wall effects on a lid driven cavity flow. Physics Fluids A (1). February 1989

Define case name

This is the solver case to be analysed


In [1]:
case_name = '3dCavity'

Define Data Location

For remote data the interaction will use ssh to securely interact with the data
This uses the reverse connection capability in paraview so that the paraview server can be submitted to a job scheduler
Note: The default paraview server connection will use port 11111


In [2]:
remote_data = True

#data_dir='/gpfs/asrc/home/acimpoeru/cube_JA/CUBE'

data_dir='/gpfs/asrc/home/acimpoeru/3dCavity_Jamil/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

Initialise Environment


In [3]:
%pylab inline
from paraview.simple import *
paraview.simple._DisableFirstRenderCameraReset()
import pylab as pl
import math
import numpy as np


Populating the interactive namespace from numpy and matplotlib
paraview version 4.1.0

Data Connection

This starts paraview server on remote host and connects


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)


[acimpoeru@login02] Executing task 'pvserver'
[acimpoeru@login02] run: sleep 2;mpiexec /gpfs/thirdparty/zenotech/home/jappa/apps/Paraview/bin/pvserver -rc --client-host=localhost -sp=11115&>/dev/null
[acimpoeru@login02] out: 
[acimpoeru@login02] out: 		   _____ ______ __  __  _____ 
[acimpoeru@login02] out: 		  / ____|  ____|  \/  |/ ____|
[acimpoeru@login02] out: 		 | |    | |__  | \  / | (___  
[acimpoeru@login02] out: 		 | |    |  __| | |\/| |\___ \ 
[acimpoeru@login02] out: 		 | |____| |    | |  | |____) |
[acimpoeru@login02] out: 		  \_____|_|    |_|  |_|_____/ 
[acimpoeru@login02] out:                               
[acimpoeru@login02] out:                               
[acimpoeru@login02] out: 
[acimpoeru@login02] out: ++++++++++++++++++++++++++++: System Data :++++++++++++++++++++++++++++
[acimpoeru@login02] out: + Hostname = login02
[acimpoeru@login02] out: + Kernel = 2.6.32-358.el6.x86_64
[acimpoeru@login02] out: + RHEL Release = Red Hat Enterprise Linux Server release 6.4 (Santiago)
[acimpoeru@login02] out: + Uptime = 14:23:24 up 112 days, 2:07, 51 users,
[acimpoeru@login02] out: + CPU = 2x Intel Xeon X5570 @ 2.93GHz
[acimpoeru@login02] out: + Memory = 49415076 kB
[acimpoeru@login02] out: ++++++++++++++++++++++++++++: User Data :++++++++++++++++++++++++++++++
[acimpoeru@login02] out: + Username = acimpoeru
[acimpoeru@login02] out: +++++++++++++++++++++++: Contact Information :+++++++++++++++++++++++++
[acimpoeru@login02] out: + in case of any problems, contact: support@cfms.org.uk
[acimpoeru@login02] out: + for feedback, contact: feedback@cfms.org.uk 
[acimpoeru@login02] out: +++++++++++++++++++++: Maintenance Information :+++++++++++++++++++++++
[acimpoeru@login02] out: + There is no planned maintenance for the cluster
[acimpoeru@login02] out: +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[acimpoeru@login02] out: 
[acimpoeru@login02] out: 

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'])

Define test conditions


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]:
pressure101325.0
temperature294.5
Reynolds No3200
Ref length150
Speed1.0
Mach No0.0

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

p_k_u_of_y = zutil.post.get_csv_data("/gpfs/asrc/home/acimpoeru/3dCavity_Jamil/JA/p_k_u_of_y.csv",header=True,remote=True,delim=',')
print p_k_u_of_y.keys()


['u/ub', 'y/h']

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

p_k_v_of_x = zutil.post.get_csv_data("/gpfs/asrc/home/acimpoeru/3dCavity_Jamil/JA/p_k_v_of_x.csv",header=True,remote=True,delim=',')
print p_k_v_of_x.keys()
#print data['y/H'], data['U/Ub']


['x/L', 'v/u_b']

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 [14]:
from zutil.post import for_each

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')

Velocity profiles y=0


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 = 3200.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(p_k_u_of_y['u/ub'],p_k_u_of_y['y/h'], color='r', label='Exp U_y Prasad and Koseff',linestyle='-',linewidth = '5',marker='o',markersize=20)
ax.plot(p_k_v_of_x['x/L'],p_k_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("images/ahmed_vel_profile.pdf")
fig.savefig("Re=3200.pdf")
#from IPython.display import FileLink, display 
#display(FileLink('images/ahmed_vel_profile.pdf'))
show()


<paraview.servermanager.CellDatatoPointData object at 0x117ec3f10>
use append poly data filter
<paraview.servermanager.CellDatatoPointData object at 0x117ec3fd0>
use append poly data filter

Calculate the RMS Error

interpolate using numpy.....See reference guide

Convergence


In [17]:
from zutil.post import residual_plot, get_case_report
residual_plot(get_case_report(case_name))
show()


---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-17-4c50d5b3ada1> in <module>()
      1 from zutil.post import residual_plot, get_case_report
----> 2 residual_plot(get_case_report(case_name))
      3 show()

/Users/andrei/Documents/zPost/python/zutil/post.pyc in residual_plot(file)
    207         ax.set_ylabel('l2norm '+my_names[i], multialignment='center')
    208 
--> 209         ax.plot(data[names[0]], data[names[i]], color='r', label=names[i])
    210 
    211 def for_each(surface,func,**kwargs):

//Applications/paraview.app/Contents/Python/matplotlib/axes.pyc in plot(self, *args, **kwargs)
   3896 
   3897 
-> 3898         self.autoscale_view(scalex=scalex, scaley=scaley)
   3899         return lines
   3900 

//Applications/paraview.app/Contents/Python/matplotlib/axes.pyc in autoscale_view(self, tight, scalex, scaley)
   1877                 y1 += delta
   1878             if not _tight:
-> 1879                 y0, y1 = ylocator.view_limits(y0, y1)
   1880             self.set_ybound(y0, y1)
   1881 

//Applications/paraview.app/Contents/Python/matplotlib/ticker.pyc in view_limits(self, vmin, vmax)
   1343         if minpos<=0 or not np.isfinite(minpos):
   1344             raise ValueError(
-> 1345                 "Data has no positive values, and therefore can not be log-scaled.")
   1346 
   1347         if vmin <= minpos:

ValueError: Data has no positive values, and therefore can not be log-scaled.

Cleaning up


In [ ]:
if remote_data:
    #print 'Disconnecting from remote paraview server connection'
    Disconnect()
    pass