Last active
August 30, 2016 11:38
-
-
Save zhulianhua/3e1c6d551d1fdb25b9de92ff1ff1a432 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| """This file contains code for use with "Think Stats", | |
| Copyright 2016 Lianhua Zhu | |
| License: GNU GPLv3 http://www.gnu.org/licenses/gpl.html | |
| """ | |
| import sys | |
| import gzip | |
| import os | |
| import argparse | |
| import numpy as np | |
| from scipy import interpolate | |
| import matplotlib | |
| import matplotlib.pyplot as plt | |
| import matplotlib.lines as mlines | |
| from scipy.io import FortranFile | |
| import glob | |
| matplotlib.use('Agg') | |
| import os | |
| def splitall(path): | |
| ''' Copied from https://www.safaribooksonline.com/library/view/python-cookbook/0596001673/ch04s16.html''' | |
| allparts = [] | |
| while 1: | |
| parts = os.path.split(path) | |
| if parts[0] == path: # sentinel for absolute paths | |
| allparts.insert(0, parts[0]) | |
| break | |
| elif parts[1] == path: # sentinel for relative paths | |
| allparts.insert(0, parts[1]) | |
| break | |
| else: | |
| path = parts[0] | |
| allparts.insert(0, parts[1]) | |
| return allparts | |
| def get_meshgrid(N): | |
| ''' parepare the mesh grid for output, note that would include the boundary (0,1) | |
| Args: | |
| N grid points used in simulation | |
| Return: | |
| X, Y | |
| ''' | |
| dx = 1.0 / N | |
| x = np.zeros(N+2) | |
| x[0] = 0 | |
| x[N+1] = 1.0 | |
| for i in range(1,N+1): | |
| x[i] = i * dx - 0.5 * dx | |
| return np.meshgrid(x, x) | |
| def read_uvp(datadir, N): | |
| ''' read the u, v from plain files | |
| Returns: | |
| two N2 * N2 arrays, u and v | |
| ''' | |
| u = np.zeros([N + 2, N + 2]) | |
| v = np.zeros([N + 2, N + 2]) | |
| p = np.zeros([N + 2, N + 2]) | |
| u[N+1, :] = 0.1 | |
| u_tmp = np.genfromtxt(os.path.join(datadir, 'u.dat')) | |
| v_tmp = np.genfromtxt(os.path.join(datadir, 'v.dat')) | |
| p_tmp = np.genfromtxt(os.path.join(datadir, 'rho.dat')) | |
| ## some data (bkg code) had been scaled before written, now scaled back | |
| if u_tmp.max() > 0.15 : | |
| u_tmp = 0.1 * u_tmp | |
| v_tmp = 0.1 * v_tmp | |
| assert u_tmp.max() < 0.15 and u_tmp.max() > 0.05 | |
| assert N == u_tmp.shape[0] | |
| u[1:-1, 1:-1] = u_tmp | |
| v[1:-1, 1:-1] = v_tmp | |
| p[1:-1, 1:-1] = (p_tmp - 1.0) / 3.0 / (0.1 * 0.1) # scale by U^2 | |
| # extrapolate pressure to boundaries | |
| p[1:-1,0] = 1.5 * p[1:-1,1] - 0.5 * p[1:-1,2] | |
| p[1:-1,N+1] = 1.5 * p[1:-1,N] - 0.5 * p[1:-1,N-1] | |
| p[0,1:-1] = 1.5 * p[1,1:-1] - 0.5 * p[2,1:-1] | |
| p[N+1,1:-1] = 1.5 * p[N,1:-1] - 0.5 * p[N-1,1:-1] | |
| p[0,0] = 1.5 * p[1,1] - 0.5 * p[2,2] | |
| p[0,N+1] = 1.5 * p[1,N] - 0.5 * p[2,N-1] | |
| p[N+1,0] = 1.5 * p[N,1] - 0.5 * p[N-1,2] | |
| p[N+1,N+1] = 1.5 * p[N,N] - 0.5 * p[N-1, N-1] | |
| return (u,v,p) | |
| def read_uvp_xuao(datadir, N): | |
| ''' read bindary unformated data from Ao Xu's fortran code ''' | |
| u = np.zeros([N + 2, N + 2]) | |
| v = np.zeros([N + 2, N + 2]) | |
| p = np.zeros([N + 2, N + 2]) | |
| u[N+1, :] = 0.1 | |
| files = glob.glob(os.path.join(datadir,"*.bin")) | |
| ff = FortranFile(files[0], 'r') | |
| uf_tmp = ff.read_reals(dtype=float) | |
| vf_tmp = ff.read_reals(dtype=float) | |
| rhof_tmp = ff.read_reals(dtype=float) | |
| u_tmp = np.reshape(uf_tmp, (-1,N), 'C') | |
| v_tmp = np.reshape(vf_tmp, (-1,N), 'C') | |
| rho_tmp = np.reshape(rhof_tmp, (-1,N), 'C') | |
| ff.close() | |
| ## some data (bkg code) had been scaled before written, now scaled back | |
| if u_tmp.max() > 0.15 : | |
| u_tmp = 0.1 * u_tmp | |
| v_tmp = 0.1 * v_tmp | |
| u[1:-1, 1:-1] = u_tmp | |
| v[1:-1, 1:-1] = v_tmp | |
| p[1:-1, 1:-1] = (rho_tmp - 1.0)/ 3.0 / (0.1 * 0.1) # scale by U^2 | |
| p[1:-1,0] = 1.5 * p[1:-1,1] - 0.5 * p[1:-1,2] | |
| p[1:-1,N+1] = 1.5 * p[1:-1,N] - 0.5 * p[1:-1,N-1] | |
| p[0,1:-1] = 1.5 * p[1,1:-1] - 0.5 * p[2,1:-1] | |
| p[N+1,1:-1] = 1.5 * p[N,1:-1] - 0.5 * p[N-1,1:-1] | |
| p[0,0] = 1.5 * p[1,1] - 0.5 * p[2,2] | |
| p[0,N+1] = 1.5 * p[1,N] - 0.5 * p[2,N-1] | |
| p[N+1,0] = 1.5 * p[N,1] - 0.5 * p[N-1,2] | |
| p[N+1,N+1] = 1.5 * p[N,N] - 0.5 * p[N-1, N-1] | |
| return (u,v,p) | |
| def get_streamfunction(u,v): | |
| ''' calculate the streamfunction \phi from u and v | |
| Return: | |
| the streamfunction 2D filed \phi | |
| ''' | |
| phi = np.zeros(u.shape) | |
| N = u.shape[1] - 2 | |
| dx = 1.0 / N | |
| # Near boundary nodes (i = 1), trapezoidal rule (|-|) | |
| for j in range(0,N+2): | |
| phi[j][1] = phi[j][0] - dx / 4.0 * (v[j][1] + v[j][0]) | |
| # Next near boundary nodes (i = 2), three point integration (|-|--|) !!!! Something wrong! | |
| # for j in range(0,N+2): | |
| # a = (v[j][2] - 2.0 * v[j][1] + v[j][0]) / (2.0 * dx * dx) | |
| # b = (4.0*v[j][1] - v[j][2] - 3.0 * v[j][0]) / (2.0 * dx) | |
| # c = v[j][0] | |
| # phi[j][2] = phi[j][0] - (8.0 * a /3.0 * dx * dx * dx + 2.0 * b * dx * dx + c * dx) | |
| # Next near boundary nodes (i = 2), trapezoidal rule (|-|) | |
| for j in range(0,N+2): | |
| #a = (v[j][2] - 2.0 * v[j][1] + v[j][0]) / (2.0 * dx * dx) | |
| #b = (4.0*v[j][1] - v[j][2] - 3.0 * v[j][0]) / (2.0 * dx) | |
| #c = v[j][0] | |
| phi[j][2] = phi[j][1] - dx / 4.0 * (v[j][1] + v[j][2]) | |
| # other nodes (i > 2), Simpson rule (|--|--|) | |
| for j in range(0,N+2): | |
| for i in range(3,N+2): | |
| phi[j][i] = phi[j][i-2] - dx / 6.0 * (v[j][i-2] + 4.0 * v[j][i-1] + v[j][i]) | |
| return phi / 0.1 * 2.0 # scale by U_w * L | |
| def get_vorticity(u,v): | |
| ''' calculate the vorticity \omega from u and v | |
| Return: | |
| the vorticity 2D filed \omega | |
| ''' | |
| omega = np.zeros(u.shape) | |
| N = u.shape[1] - 2 | |
| dx = 1.0 / N | |
| # inner nodes | |
| for j in range(2,N): | |
| for i in range(2, N): | |
| pupy = (u[j+1][i] - u[j-1][i]) / 3.0 + (u[j+1][i-1] - u[j-1][i-1]) / 12.0 + (u[j+1][i+1] - u[j-1][i+1]) / 12.0 | |
| pvpx = (v[j][i+1] - v[j][i-1]) / 3.0 + (v[j+1][i+1] - v[j+1][i-1]) / 12.0 + (v[j-1][i+1] - v[j-1][i-1]) / 12.0 | |
| omega[j][i] = (pvpx - pupy)/dx | |
| # near boundary nodes | |
| for i in range(2, N): | |
| # near bottom | |
| pupy = (3.0 * u[1][i] + u[2][i]) / 3.0 | |
| pvpx = (v[1][i+1] - v[1][i-1]) / 2.0 | |
| omega[1][i] = (pvpx - pupy)/dx | |
| # near top | |
| pupy = (3.0 * u[N][i] + u[N-1][i] - 4.0 * 0.1) / 3.0 | |
| pvpx = (v[N][i+1] - v[N][i-1]) / 2.0 | |
| omega[N][i] =(pvpx - pupy)/dx | |
| for j in range(2,N): | |
| # near left | |
| pvpx = (3.0 * v[j][1] + v[j][2]) / 3.0 | |
| pupy = (u[j+1][1] - u[j-1][1]) / 2.0 | |
| omega[j][1] = (pvpx - pupy)/dx | |
| # near right | |
| pvpx = (3.0 * v[j][N] + v[j][N-1]) / 3.0 | |
| pupy = (u[j+1][N] - u[j-1][N]) / 2.0 | |
| omega[j][N] = (pvpx - pupy)/dx | |
| # near corner nodes | |
| # (1,1) | |
| pupy = (3.0 * u[1][1] + u[2][1]) / 3.0 | |
| pvpx = (3.0 * v[1][1] + v[1][2]) / 3.0 | |
| omega[1][1] = (pvpx - pupy)/dx | |
| # (1,N) | |
| pupy = (3.0 * u[1][N] + u[2][N]) / 3.0 | |
| pvpx = (3.0 * v[1][N] + v[1][N-1]) / 3.0 | |
| omega[1][N] = (pvpx - pupy)/dx | |
| # (N,1) | |
| pupy = (3.0 * u[N][1] + u[N-1][1] - 4.0 * 0.1) / 3.0 | |
| pvpx = (3.0 * v[N][1] + v[N][2]) / 3.0 | |
| omega[N][1] = (pvpx - pupy)/dx | |
| # (N,N) | |
| pupy = (3.0 * u[N][N] + u[N-1][N] - 4.0 * 0.1) / 3.0 | |
| pvpx = (3.0 * v[N][N] + v[N][N-1]) / 3.0 | |
| omega[N][N] = (pvpx - pupy)/dx | |
| # boundary nodes | |
| for i in range(1, N+1): | |
| omega[0][i] = 2.0 * u[1][i] / dx | |
| omega[N+1][i] = 2.0 * (0.1 - u[N][i]) / dx | |
| for j in range(1, N+1): | |
| omega[j][0] = 2.0 * v[j][1] / dx | |
| omega[j][N+1] = -2.0 * v[j][N] / dx | |
| # corner nodes | |
| omega[0][0] = 0 | |
| omega[0][N+1] = 0 | |
| omega[N+1][0] = 0.1 * 2.0 / dx | |
| omega[N+1][N+1] = 0.1 * 2.0 / dx | |
| return omega / 0.1 # scale by U_w / L | |
| def find_vortex_center(xca, yca, vname, X, Y, p, omega, phi): | |
| N = X.shape[1] - 2 | |
| dx = 1.0 / N | |
| ja = int((yca - 0.5 * dx)/dx) + 1 | |
| ia = int((xca - 0.5 * dx)/dx) + 1 | |
| numn = 5 | |
| il = ia-numn | |
| iu = ia+numn+1 | |
| jl = ja-numn | |
| ju = ja+numn+1 | |
| phi_loc = phi[jl:ju,il:iu] | |
| X_loc = X[jl:ju,il:iu] | |
| Y_loc = Y[jl:ju,il:iu] | |
| p_loc = p[jl:ju,il:iu] | |
| omega_loc = omega[jl:ju,il:iu] | |
| fphi = interpolate.interp2d(X_loc, Y_loc, phi_loc, kind='cubic') | |
| fomega = interpolate.interp2d(X_loc, Y_loc, omega_loc, kind='cubic') | |
| fp = interpolate.interp2d(X_loc, Y_loc, p_loc, kind='cubic') | |
| if vname == "BL": | |
| print "=============== Bottom Left Vortex ================" | |
| print "%.6f %.6f %.6f %.6f %.6f"%(xca,yca,fp(xca,yca)*100,fphi(xca,yca)*1.0e4,fomega(xca,yca)*10) | |
| elif vname == "BR": | |
| print "=============== Bottom Right Vortex ================" | |
| print "%.6f %.6f %.6f %.6f %.6f"%(xca,yca,fp(xca,yca)*100,fphi(xca,yca)*1.0e3,fomega(xca,yca)) | |
| else: | |
| print "=============== Primary Vortex ================" | |
| print "%.6f %.6f %.6f %.6f %.6f"%(xca,yca,fp(xca,yca),fphi(xca,yca),fomega(xca,yca)) | |
| def write_tecplot(dir, X, Y, u,v,p,omega,phi): | |
| DATA = np.array([X,Y,u,v,p,omega,phi]) | |
| N2 = X.shape[0] | |
| fp = open(os.path.join(dir,"results.dat"), "w") | |
| fp.write(" VARIABLES = X, Y, U, V, P, omega, phi\n") | |
| fp.write("ZONE I = %d , J = %d DATAPACKING=BLOCK\n"%(N2,N2)) | |
| for dt in DATA: | |
| for j in range(N2): | |
| for i in range(N2): | |
| fp.write("%.18e\t"%dt[j,i]) | |
| fp.write("\n") | |
| def write_plain_data(dir,X,Y,u,v,p,omega,phi): | |
| np.savetxt(os.path.join(dir,"X.dat"), X) | |
| np.savetxt(os.path.join(dir,"Y.dat"), Y) | |
| np.savetxt(os.path.join(dir,"UA.dat"), u) | |
| np.savetxt(os.path.join(dir,"VA.dat"), v) | |
| np.savetxt(os.path.join(dir,"P.dat"), p) | |
| np.savetxt(os.path.join(dir,"omega.dat"), omega) | |
| np.savetxt(os.path.join(dir,"phi.dat"), phi) | |
| def plot_contour(dir, X, Y, omega, phi): | |
| ''' plot streamline and vorticity contours and save to file''' | |
| fig = plt.figure(figsize=(5,5)) | |
| axes = fig.add_subplot(111) | |
| axes.xaxis.set_ticklabels([]) | |
| axes.yaxis.set_ticklabels([]) | |
| phi_level = np.array([-0.1175, -0.115, -0.11, -0.1, -9e-2, -7e-2, -5e-2, -3e-2, -1e-2, \ | |
| -1e-4, -1e-5, -1e-7, -1e-10, 1e-8, 1e-7, 1e-6, 1e-5, 5e-5, 1e-4, 2.5e-4, 5e-4, 1e-3, 1.5e-3, 3e-3]) | |
| omega_level = np.array([-3.0, -2.0, -1.0, -0.5, 0, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0]) | |
| CS = axes.contour(X,Y,-omega, omega_level, colors='k', linewidths=1) | |
| #plt.clabel(CS, inline=1, fontsize=10) | |
| plt.tight_layout() | |
| fig.savefig(os.path.join(dir,'_'.join(splitall(dir))+"_omega.pdf")) | |
| plt.cla() | |
| axes.xaxis.set_ticklabels([]) | |
| axes.yaxis.set_ticklabels([]) | |
| axes.contour(X,Y,phi, phi_level, colors='k', linewidths=1, linestyles='solid') | |
| fig.savefig(os.path.join(dir,'_'.join(splitall(dir))+"_psi.pdf")) | |
| def get_profile(casedir): | |
| #if os.path.exists(os.path.join(casedir, 'uy.dat')): | |
| if False: | |
| uy = np.genfromtxt(os.path.join(casedir, 'uy.dat')) | |
| vx = np.genfromtxt(os.path.join(casedir, 'vx.dat')) | |
| xi = uy[0,:] | |
| ui = uy[1,:] | |
| xi = vx[0,:] | |
| vi = vx[1,:] | |
| else: | |
| U = np.genfromtxt(os.path.join(casedir, 'UA.dat')) | |
| V = np.genfromtxt(os.path.join(casedir, 'VA.dat')) | |
| X = np.genfromtxt(os.path.join(casedir, 'X.dat')) | |
| Y = np.genfromtxt(os.path.join(casedir, 'Y.dat')) | |
| N = U.shape[0] - 2 | |
| Nh = N / 2 | |
| #fU = interpolate.interp2d(X, Y, U, kind='cubic') | |
| #xi = 0.5 | |
| xi = Y[:,0] | |
| #ui = 10.0 * fU(xi,xi) # data .dat file of velocity are non-scaled | |
| ui = 10.0/16 * (-U[:,Nh-1] + 9.0 * U[:,Nh] + 9.0 * U[:,Nh+1] - U[:,Nh+2]) | |
| np.savetxt(os.path.join(casedir, 'uy.dat'), [xi, ui]) | |
| #fV = interpolate.interp2d(X, Y, V, kind='cubic') | |
| xi = X[0,:] | |
| #xi = 0.5 | |
| #vi = 10.0 * fV(xi,xi) # data .dat file of velocity are non-scaled | |
| vi = 10.0/16 * (-V[Nh-1,:] + 9.0 * V[Nh,:] + 9.0 * V[Nh+1,:] - V[Nh+2,:]) | |
| #print vi.shape | |
| #vi = vi[0,:] | |
| np.savetxt(os.path.join(casedir, 'vx.dat'), [xi, vi]) | |
| return ({'name': casedir, 'xi': xi, 'ui':ui}, \ | |
| {'name': casedir, 'xi': xi, 'vi':vi}) | |
| # def get_vprofile(case, X, Y, V): | |
| # N = V.shape[0] - 2 | |
| # fV = interpolate.interp2d(X, Y, V, kind='cubic') | |
| # xi = X[0,:] | |
| # xi = Y[:,0] | |
| # vi = fV(xi,xi) | |
| # return {'name': case, 'xi': xi, 'vi':vi} | |
| def plot_vprofiles(): | |
| # prepare Ghia's data | |
| # prepare Ghia's data | |
| GhiaRe1000yu = np.genfromtxt("ghia_re1000_XC.dat") | |
| GhiaRe1000xv = np.genfromtxt("ghia_re1000_YC.dat") | |
| GhiaRe5000yu = np.genfromtxt("ghia_re5000_XC.dat") | |
| GhiaRe5000xv = np.genfromtxt("ghia_re5000_YC.dat") | |
| # y-u profiles | |
| fig = plt.figure(figsize=(8,3)) | |
| axes = fig.add_subplot(111) | |
| axes.xaxis.set_ticklabels([]) | |
| #axes.yaxis.set_ticklabels([]) | |
| plt.plot([0.6, 0.6], [0, 1.0], 'k--', lw=0.5) | |
| plt.plot([1.2, 1.2], [0, 1.0], 'k--', lw=0.5) | |
| plt.plot([1.8, 1.8], [0, 1.0], 'k--', lw=0.5) | |
| plt.plot([2.4, 2.4], [0, 1.0], 'k--', lw=0.5) | |
| plt.plot([3.0, 3.0], [0, 1.0], 'k--', lw=0.5) | |
| plt.plot([3.6, 3.6], [0, 1.0], 'k--', lw=0.5) | |
| plt.plot([4.2, 4.2], [0, 1.0], 'k--', lw=0.5) | |
| plt.plot([4.8, 4.8], [0, 1.0], 'k--', lw=0.5) | |
| plt.plot(GhiaRe1000xv[:,1]+0.6, GhiaRe1000xv[:,0], 'ko', mfc='none',ms=4) | |
| plt.plot(GhiaRe1000xv[:,1]+1.2, GhiaRe1000xv[:,0], 'ko', mfc='none',ms=4) | |
| plt.plot(GhiaRe1000xv[:,1]+1.8, GhiaRe1000xv[:,0], 'ko', mfc='none',ms=4) | |
| plt.plot(GhiaRe5000xv[:,1]+3.0, GhiaRe5000xv[:,0], 'ko', mfc='none',ms=4) | |
| plt.plot(GhiaRe5000xv[:,1]+3.6, GhiaRe5000xv[:,0], 'ko', mfc='none',ms=4) | |
| plt.plot(GhiaRe5000xv[:,1]+4.2, GhiaRe5000xv[:,0], 'ko', mfc='none',ms=4) | |
| plt.ylabel(r"$x/L$", fontsize=14) | |
| plt.xlabel(r"$v/u_\textrm{\small{w}}$", fontsize=14) | |
| axes.xaxis.set_label_coords(0.95, -0.1) | |
| # Hide the right and top spines | |
| axes.spines['right'].set_visible(False) | |
| axes.spines['top'].set_visible(False) | |
| # Only show ticks on the left and bottom spines | |
| axes.yaxis.set_ticks_position('left') | |
| axes.xaxis.set_ticks_position('bottom') | |
| axes.set_ylim([0,1.05]) | |
| axes.set_xlim([0,4.5]) | |
| plt.yticks([0, 0.2, 0.4, 0.6, 0.8, 1.0], \ | |
| ['$0.0$', '$0.2$', '$0.4$', '$0.6$', '$0.8$', '$1.0$']) | |
| plt.xticks([0.6, 1.2, 1.8, 2.4, 3.0, 3.6, 4.2, 4.8], \ | |
| ['$0$','$0$','$0$','$0$','$0$','$0$','$0$','$0.6$']) | |
| axes.tick_params(axis='both', which='major', labelsize=14) | |
| # add profiles | |
| uypf, vxpf = get_profile("bkg/Re1000/N64") | |
| plt.plot(vxpf['vi']+0.6, vxpf['xi'], 'b-', lw=1) | |
| uypf, vxpf = get_profile("dugks/Re1000/N64") | |
| plt.plot(vxpf['vi']+0.6, vxpf['xi'], 'r--',lw=1) | |
| uypf, vxpf = get_profile("bkg/Re1000/N128") | |
| plt.plot(vxpf['vi']+1.2, vxpf['xi'], 'b-', lw=1) | |
| uypf, vxpf = get_profile("dugks/Re1000/N128") | |
| plt.plot(vxpf['vi']+1.2, vxpf['xi'], 'r--',lw=1) | |
| uypf, vxpf = get_profile("bkg/Re1000/N256") | |
| plt.plot(vxpf['vi']+1.8, vxpf['xi'], 'b-', lw=1) | |
| uypf, vxpf = get_profile("dugks/Re1000/N256") | |
| plt.plot(vxpf['vi']+1.8, vxpf['xi'], 'r--',lw=1) | |
| uypf, vxpf = get_profile("bkg/Re5000/N128") | |
| plt.plot(vxpf['vi']+3.0, vxpf['xi'], 'b-', lw=1) | |
| uypf, vxpf = get_profile("dugks/Re5000/N128") | |
| plt.plot(vxpf['vi']+3.0, vxpf['xi'], 'r--',lw=1) | |
| uypf, vxpf = get_profile("bkg/Re5000/N256") | |
| plt.plot(vxpf['vi']+3.6, vxpf['xi'], 'b-', lw=1) | |
| uypf, vxpf = get_profile("dugks/Re5000/N256") | |
| plt.plot(vxpf['vi']+3.6, vxpf['xi'], 'r--',lw=1) | |
| uypf, vxpf = get_profile("bkg/Re5000/N512") | |
| plt.plot(vxpf['vi']+4.2, vxpf['xi'], 'b-', lw=1) | |
| uypf, vxpf = get_profile("dugks/Re5000/N256") | |
| plt.plot(vxpf['vi']+4.2, vxpf['xi'], 'r--',lw=1) | |
| plt.text(0.9, 1.1, r"$ \textnormal{Re}=1000$",fontsize=14) | |
| plt.text(3.3, 1.1, r"$\textnormal{Re}=5000$",fontsize=14) | |
| props = {'ha':'left', 'va':'top'} | |
| plt.text(0.22, 0.65, r"N=64", rotation=-50,fontsize=14) | |
| plt.text(0.85, 0.65, r"N=128",rotation=-50,fontsize=14) | |
| plt.text(1.45, 0.65, r"N=256",rotation=-50,fontsize=14) | |
| plt.text(2.6, 0.65, r"N=128",rotation=-50,fontsize=14) | |
| plt.text(3.2, 0.65, r"N=256",rotation=-50,fontsize=14) | |
| plt.text(3.8, 0.65, r"N=512",rotation=-50,fontsize=14) | |
| blue_line = mlines.Line2D([], [], color='blue', marker='None', lw=1, label='Bardow') | |
| red_line = mlines.Line2D([], [], color='red', marker='None', lw=1, label='DUGKS', ls='dashed') | |
| black_line = mlines.Line2D([], [], color='black', marker='o', mfc='none', lw=1, ls='None', label='Ghia et al.', ms=4) | |
| legend = plt.legend(loc='upper left', handles=[blue_line, red_line, black_line], ncol=3, fontsize=14, bbox_to_anchor=(0.1, -0.05)) | |
| # legend position | |
| plt.tight_layout() | |
| fig.savefig("cavity_vx.pdf") | |
| def plot_uprofiles(): | |
| # prepare Ghia's data | |
| GhiaRe1000yu = np.genfromtxt("ghia_re1000_XC.dat") | |
| GhiaRe1000xv = np.genfromtxt("ghia_re1000_YC.dat") | |
| GhiaRe5000yu = np.genfromtxt("ghia_re5000_XC.dat") | |
| GhiaRe5000xv = np.genfromtxt("ghia_re5000_YC.dat") | |
| # y-u profiles | |
| fig = plt.figure(figsize=(8,3)) | |
| axes = fig.add_subplot(111) | |
| axes.xaxis.set_ticklabels([]) | |
| #axes.yaxis.set_ticklabels([]) | |
| plt.plot([-0.5, -0.5], [0, 1.0], 'k--', lw=0.5) | |
| plt.plot([0.5, 0.5], [0, 1.0], 'k--', lw=0.5) | |
| plt.plot([1.5, 1.5], [0, 1.0], 'k--', lw=0.5) | |
| plt.plot([2.5, 2.5], [0, 1.0], 'k--', lw=0.5) | |
| plt.plot([3.5, 3.5], [0, 1.0], 'k--', lw=0.5) | |
| plt.plot([4.5, 4.5], [0, 1.0], 'k--', lw=0.5) | |
| plt.plot([5.5, 5.5], [0, 1.0], 'k--', lw=0.5) | |
| plt.plot([6.5, 6.5], [0, 1.0], 'k--', lw=0.5) | |
| plt.plot(GhiaRe1000yu[:,1]-0.5, GhiaRe1000yu[:,0], 'ko', mfc='none',ms=4) | |
| plt.plot(GhiaRe1000yu[:,1]+0.5, GhiaRe1000yu[:,0], 'ko', mfc='none',ms=4) | |
| plt.plot(GhiaRe1000yu[:,1]+1.5, GhiaRe1000yu[:,0], 'ko', mfc='none',ms=4) | |
| plt.plot(GhiaRe5000yu[:,1]+3.5, GhiaRe5000yu[:,0], 'ko', mfc='none',ms=4) | |
| plt.plot(GhiaRe5000yu[:,1]+4.5, GhiaRe5000yu[:,0], 'ko', mfc='none',ms=4) | |
| plt.plot(GhiaRe5000yu[:,1]+5.5, GhiaRe5000yu[:,0], 'ko', mfc='none',ms=4) | |
| plt.ylabel(r"$y/L$", fontsize=14) | |
| plt.xlabel(r"$u/u_\textrm{\small{w}}$", fontsize=14) | |
| axes.xaxis.set_label_coords(0.95, -0.1) | |
| # Hide the right and top spines | |
| axes.spines['right'].set_visible(False) | |
| axes.spines['top'].set_visible(False) | |
| # Only show ticks on the left and bottom spines | |
| axes.yaxis.set_ticks_position('left') | |
| axes.xaxis.set_ticks_position('bottom') | |
| axes.set_ylim([0,1.05]) | |
| axes.set_xlim([-1.2,6.7]) | |
| plt.yticks([0, 0.2, 0.4, 0.6, 0.8, 1.0], \ | |
| ['$0.0$', '$0.2$', '$0.4$', '$0.6$', '$0.8$', '$1.0$']) | |
| plt.xticks([-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5], \ | |
| ['$0$','$0$','$0$','$0$','$0$','$0$','$0$','$1$']) | |
| axes.tick_params(axis='both', which='major', labelsize=14) | |
| # add profiles | |
| uypf, vxpf = get_profile("bkg/Re1000/N64") | |
| plt.plot(uypf['ui']-0.5, uypf['xi'], 'b-', lw=1) | |
| uypf, vxpf = get_profile("dugks/Re1000/N64") | |
| plt.plot(uypf['ui']-0.5, uypf['xi'], 'r--',lw=1) | |
| uypf, vxpf = get_profile("bkg/Re1000/N128") | |
| plt.plot(uypf['ui']+0.5, uypf['xi'], 'b-', lw=1) | |
| uypf, vxpf = get_profile("dugks/Re1000/N128") | |
| plt.plot(uypf['ui']+0.5, uypf['xi'], 'r--',lw=1) | |
| uypf, vxpf = get_profile("bkg/Re1000/N256") | |
| plt.plot(uypf['ui']+1.5, uypf['xi'], 'b-', lw=1) | |
| uypf, vxpf = get_profile("dugks/Re1000/N256") | |
| plt.plot(uypf['ui']+1.5, uypf['xi'], 'r--',lw=1) | |
| uypf, vxpf = get_profile("bkg/Re5000/N128") | |
| plt.plot(uypf['ui']+3.5, uypf['xi'], 'b-', lw=1) | |
| uypf, vxpf = get_profile("dugks/Re5000/N128") | |
| plt.plot(uypf['ui']+3.5, uypf['xi'], 'r--',lw=1) | |
| uypf, vxpf = get_profile("bkg/Re5000/N256") | |
| plt.plot(uypf['ui']+4.5, uypf['xi'], 'b-', lw=1) | |
| uypf, vxpf = get_profile("dugks/Re5000/N256") | |
| plt.plot(uypf['ui']+4.5, uypf['xi'], 'r--',lw=1) | |
| uypf, vxpf = get_profile("bkg/Re5000/N512") | |
| plt.plot(uypf['ui']+5.5, uypf['xi'], 'b-', lw=1) | |
| uypf, vxpf = get_profile("dugks/Re5000/N256") | |
| plt.plot(uypf['ui']+5.5, uypf['xi'], 'r--',lw=1) | |
| plt.text(0.5, 1.1, r"$ \textnormal{Re}=1000$",fontsize=14) | |
| plt.text(4.5, 1.1, r"$\textnormal{Re}=5000$",fontsize=14) | |
| props = {'ha':'left', 'va':'top'} | |
| plt.text(-0.4, 0.75, r"N=64", rotation=60,fontsize=14) | |
| plt.text(0.65, 0.75, r"N=128",rotation=60,fontsize=14) | |
| plt.text(1.65, 0.75, r"N=256",rotation=60,fontsize=14) | |
| plt.text(3.65, 0.75, r"N=128", rotation=60,fontsize=14) | |
| plt.text(4.65, 0.75, r"N=256",rotation=60,fontsize=14) | |
| plt.text(5.65, 0.75, r"N=512",rotation=60,fontsize=14) | |
| blue_line = mlines.Line2D([], [], color='blue', marker='None', lw=1, label='Bardow') | |
| red_line = mlines.Line2D([], [], color='red', marker='None', lw=1, label='DUGKS', ls='dashed') | |
| black_line = mlines.Line2D([], [], color='black', marker='o', mfc='none', lw=1, ls='None', label='Ghia et al.', ms=4) | |
| legend = plt.legend(loc='upper left', handles=[blue_line, red_line, black_line], ncol=3, fontsize=14, bbox_to_anchor=(0.1, -0.05)) | |
| # legend position | |
| plt.tight_layout() | |
| fig.savefig("cavity_uy.pdf") | |
| def main(): | |
| parser = argparse.ArgumentParser() | |
| parser.add_argument("-N", "--ny", type=int, required=True, help="grid number") | |
| parser.add_argument("-C", "--case",required=True, help="case dir") | |
| args = parser.parse_args() | |
| X, Y = get_meshgrid(args.ny) | |
| u, v, p = read_uvp_xuao(args.case, args.ny) | |
| omega = get_vorticity(u,v) | |
| phi = get_streamfunction(u,v) | |
| #find_vortex_center(X, Y, p, omega, phi) | |
| #print "=============== Re=1000 ==============" | |
| #print "=============== Primary Vortex ================" | |
| #find_vortex_center(0.530835, 0.565225, "C", X, Y, p, omega, phi) # DUGKS N 256 | |
| #find_vortex_center(0.531099, 0.565133, "C", X, Y, p, omega, phi) # DUGKS N 128 | |
| #find_vortex_center(0.532409, 0.565339, "C", X, Y, p, omega, phi) # DUGKS N 64 | |
| #find_vortex_center(0.530716, 0.565158, "C", X, Y, p, omega, phi) # Bardow N 256 | |
| #find_vortex_center(0.530339, 0.565303, "C", X, Y, p, omega, phi) # Badow N 128 | |
| #find_vortex_center(0.530152, 0.571904, "C", X, Y, p, omega, phi) # Bardow N 64 | |
| #find_vortex_center(0.531297, 0.565105, "C", X, Y, p, omega, phi) # Ao Xu N 257 | |
| find_vortex_center(0.530854, 0.565227, "C", X, Y, p, omega, phi) # Ao Xu N 2049 | |
| #print "=============== Bottom Left Vortex ================" | |
| #find_vortex_center(0.083303, 0.078047, "BL", X, Y, p, omega, phi) # DUGKS N 256 | |
| #find_vortex_center(0.083462, 0.077816, "BL", X, Y, p, omega, phi) # DUGKS N 128 | |
| #find_vortex_center(0.084498, 0.076337, "BL", X, Y, p, omega, phi) # DUGKS N 64 | |
| #find_vortex_center(0.083272, 0.078094, "BL", X, Y, p, omega, phi) # Bardow N 256 | |
| #find_vortex_center(0.083055, 0.077727, "BL", X, Y, p, omega, phi) # Bardow N 128 | |
| #find_vortex_center(0.078245, 0.068497, "BL", X, Y, p, omega, phi) # Bardow N 64 | |
| #find_vortex_center(0.083021, 0.077554, "BL", X, Y, p, omega, phi) # Ao Xu N 257 | |
| find_vortex_center(0.083248, 0.078037, "BL", X, Y, p, omega, phi) # Ao Xu N 2049 | |
| #print "=============== Bottom Right Vortex ================" | |
| #find_vortex_center(0.863983, 0.111790, "BR", X, Y, p, omega, phi) # DUGKS N 256 | |
| #find_vortex_center(0.865040, 0.111401, "BR", X, Y, p, omega, phi) # DUGKS N 128 | |
| #find_vortex_center(0.871912, 0.110602, "BR", X, Y, p, omega, phi) # DUGKS N 64 | |
| #find_vortex_center(0.863929, 0.111679, "BR", X, Y, p, omega, phi) # Bardow N 256 | |
| #find_vortex_center(0.864855, 0.111114, "BR", X, Y, p, omega, phi) # Bardow N 128 | |
| #find_vortex_center(0.871587, 0.113014, "BR", X, Y, p, omega, phi) # Bardow N 64 | |
| #find_vortex_center(0.864586, 0.112062, "BR", X, Y, p, omega, phi) # Ao Xu N 256 | |
| find_vortex_center(0.864114, 0.111827, "BR", X, Y, p, omega, phi) # Ao Xu N 2049 | |
| #print "=============== Re=5000 ==============" | |
| #print "=============== Primary Vortex ================" | |
| #find_vortex_center(0.515335, 0.535248, "C", X, Y, p, omega, phi) # DUGKS N 256, Re=5000 | |
| #find_vortex_center(0.516536, 0.535408, "C", X, Y, p, omega, phi) # DUGKS N 128 | |
| #find_vortex_center(0.514993, 0.535204, "C", X, Y, p, omega, phi) # Bardow N 512 | |
| #find_vortex_center(0.514735, 0.535563, "C", X, Y, p, omega, phi) # Bardow N 256 | |
| #find_vortex_center(0.514070, 0.538940, "C", X, Y, p, omega, phi) # Bardow N 128 | |
| #write_tecplot(args.case, X,Y,u,v,p,omega,phi) | |
| #write_plain_data(args.case,X,Y,u,v,p,omega,phi) | |
| #print "=============== plot contours ==============" | |
| #plot_contour(args.case, X,Y,omega,phi) | |
| #print "=============== plot velocity profiles =======" | |
| #plot_uprofiles() | |
| #plot_vprofiles() | |
| if __name__ == '__main__': | |
| main() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
To run example: