{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# A Simple Staggered FV Code for the Navier-Stokes Equations\n",
"### Tony Saad
Assistant Professor of Chemical Engineering
University of Utah"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib.animation as animation\n",
"from matplotlib import cm\n",
"plt.rcParams['animation.html'] = 'html5'\n",
"import time"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Reynolds Number: 1000.0\n",
"dt= 0.0004\n"
]
}
],
"source": [
"# define some boiler plate\n",
"nx = 16\n",
"ny = 16\n",
"# dt = 1.0e-3 # timestep size\n",
"ν = 0.01 # dynamic viscosity\n",
"lx = 1.0\n",
"ly = 1.0\n",
"dx = lx/nx\n",
"dy = ly/ny\n",
"t = 0.0\n",
"\n",
"# cell centered coordinates\n",
"xx = np.linspace(dx/2.0,lx - dx/2.0,nx, endpoint=True)\n",
"yy = np.linspace(dy/2.0,ly - dy/2.0,ny, endpoint=True)\n",
"xcc, ycc = np.meshgrid(xx,yy)\n",
"\n",
"# x-staggered coordinates\n",
"xxs = np.linspace(0,lx,nx+1, endpoint=True)\n",
"xu,yu = np.meshgrid(xxs, yy)\n",
"\n",
"# y-staggered coordinates\n",
"yys = np.linspace(0,ly,ny+1, endpoint=True)\n",
"xv,yv = np.meshgrid(xx, yys)\n",
"\n",
"Ut = 10.0\n",
"Ub = 0.0\n",
"Vl = 0.0\n",
"Vr = 0.0\n",
"print('Reynolds Number:', Ut*lx/ν)\n",
"dt = min(0.25*dx*dx/ν, 4.0*ν/Ut/Ut)\n",
"print('dt=', dt)\n",
"\n",
"# initialize velocities - we stagger everything in the negative direction. A scalar cell owns its minus face, only.\n",
"# Then, for example, the u velocity field has a ghost cell at x0 - dx and the plus ghost cell at lx\n",
"u = np.zeros([ny+2, nx+2]) # include ghost cells\n",
"\n",
"# # same thing for the y-velocity component\n",
"v = np.zeros([ny +2, nx+2]) # include ghost cells\n",
"\n",
"ut = np.zeros_like(u)\n",
"vt = np.zeros_like(u) \n",
"\n",
"# initialize the pressure\n",
"p = np.zeros([nx+2,ny+2]); # include ghost cells\n",
"\n",
"# a bunch of lists for animation purposes\n",
"usol=[]\n",
"usol.append(u)\n",
"\n",
"vsol=[]\n",
"vsol.append(v)\n",
"\n",
"psol = []\n",
"psol.append(p)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# build pressure coefficient matrix\n",
"Ap = np.zeros([ny+2,nx+2])\n",
"Ae = 1.0/dx/dx*np.ones([ny+2,nx+2])\n",
"As = 1.0/dy/dy*np.ones([ny+2,nx+2])\n",
"An = 1.0/dy/dy*np.ones([ny+2,nx+2])\n",
"Aw = 1.0/dx/dx*np.ones([ny+2,nx+2])\n",
"# set left wall coefs\n",
"Aw[1:-1,1] = 0.0\n",
"# set right wall coefs\n",
"Ae[1:-1,-2] = 0.0\n",
"# set top wall coefs\n",
"An[-2,1:-1] = 0.0\n",
"# set bottom wall coefs\n",
"As[1,1:-1] = 0.0\n",
"Ap = -(Aw + Ae + An + As)\n",
"\n",
"def pressure_poisson2(p, b, dx, dy):\n",
" pn = np.empty_like(p)\n",
" it = 0\n",
" err = 1e5\n",
" tol = 1e-8\n",
" maxit = 100\n",
" β = 1.1\n",
" while err > tol and it < maxit:\n",
" pn = p.copy() \n",
" for i in range(1,nx+1):\n",
" for j in range(1,ny+1):\n",
" ap = Ap[j,i]\n",
" an = An[j,i]\n",
" aso = As[j,i]\n",
" ae = Ae[j,i]\n",
" aw = Aw[j,i]\n",
" rhs = b[j,i] - 1.0*(ae*p[j,i+1] + aw*p[j,i-1] + an*p[j+1,i] + aso*p[j-1,i])\n",
" p[j,i] = β*rhs/ap + (1-β)*p[j,i]\n",
" err = np.linalg.norm(p - pn)\n",
" it += 1\n",
"# print('Poisson Error:', err) \n",
" return p, err"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# USE sparse solver\n",
"import scipy.linalg\n",
"import scipy.sparse\n",
"import scipy.sparse.linalg\n",
"# build pressure coefficient matrix\n",
"Ap = np.zeros([ny,nx])\n",
"Ae = 1.0/dx/dx*np.ones([ny,nx])\n",
"As = 1.0/dy/dy*np.ones([ny,nx])\n",
"An = 1.0/dy/dy*np.ones([ny,nx])\n",
"Aw = 1.0/dx/dx*np.ones([ny,nx])\n",
"# set left wall coefs\n",
"Aw[:,0] = 0.0\n",
"# set right wall coefs\n",
"Ae[:,-1] = 0.0\n",
"# set top wall coefs\n",
"An[-1,:] = 0.0\n",
"# set bottom wall coefs\n",
"As[0,:] = 0.0\n",
"Ap = -(Aw + Ae + An + As)\n",
"\n",
"n = nx*ny\n",
"d0 = Ap.reshape(n)\n",
"# print(d0)\n",
"de = Ae.reshape(n)[:-1]\n",
"# print(de)\n",
"dw = Aw.reshape(n)[1:]\n",
"# print(dw)\n",
"ds = As.reshape(n)[nx:]\n",
"# print(ds)\n",
"dn = An.reshape(n)[:-nx]\n",
"# print(dn)\n",
"A1 = scipy.sparse.diags([d0, de, dw, dn, ds], [0, 1, -1, nx, -nx], format='csr')\n",
"plt.matshow((A1.toarray()))"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"total time 0.11582800000000049 s\n",
"time per timestep = 0.01158280000000005 s\n",
"mom assembly time 0.04313199999999817 s\n",
"solver time 0.07269600000000231 s\n",
"solver fraction = 63.0 %\n"
]
}
],
"source": [
"# while t < tend:\n",
"t0 = time.clock()\n",
"momtime = 0.0\n",
"solvertime = 0.0\n",
"nsteps = 10\n",
"for n in range(0,nsteps):\n",
" # left wall\n",
" u[1:-1,1] = 0.0\n",
" # right wall\n",
" u[1:-1,-1] = 0.0\n",
" # top wall\n",
" u[-1,1:] = 2.0*Ut - u[-2,1:]\n",
" # bottom wall\n",
" u[0,1:] = 2.0*Ub - u[1,1:]\n",
"\n",
" # left wall\n",
" v[1:,0] = 2.0*Vl - v[1:,1]\n",
" # right wall\n",
" v[1:,-1] = 2.0*Vr - v[1:,-2]\n",
" # bottom wall\n",
" v[1,1:-1] = 0.0\n",
" # top wall\n",
" v[-1,1:-1] = 0.0 \n",
" \n",
"\n",
" # do x-momentum first - u is of size (nx + 2) x (ny + 2) - only need to do the interior points\n",
" tic = time.clock()\n",
" for i in range(2,nx+1):\n",
" for j in range(1,ny+1):\n",
" ue = 0.5*(u[j,i+1] + u[j,i])\n",
" uw = 0.5*(u[j,i] + u[j,i-1]) \n",
" \n",
" un = 0.5*(u[j+1,i] + u[j,i])\n",
" us = 0.5*(u[j,i] + u[j-1,i]) \n",
" \n",
" vn = 0.5*(v[j+1,i] + v[j+1,i-1])\n",
" vs = 0.5*(v[j,i] + v[j,i-1])\n",
" \n",
" # convection = - d(uu)/dx - d(vu)/dy\n",
" convection = - (ue*ue - uw*uw)/dx - (un*vn - us*vs)/dy\n",
" \n",
" # diffusion = d2u/dx2 + d2u/dy2\n",
" diffusion = ν*( (u[j,i+1] - 2.0*u[j,i] + u[j,i-1])/dx/dx + (u[j+1,i] - 2.0*u[j,i] + u[j-1,i])/dy/dy )\n",
" \n",
" ut[j,i] = u[j,i] + dt *(convection + diffusion)\n",
" \n",
" # do y-momentum - only need to do interior points\n",
" for i in range(1,nx+1):\n",
" for j in range(2,ny+1):\n",
" ve = 0.5*(v[j,i+1] + v[j,i])\n",
" vw = 0.5*(v[j,i] + v[j,i-1]) \n",
" \n",
" ue = 0.5*(u[j,i+1] + u[j-1,i+1])\n",
" uw = 0.5*(u[j,i] + u[j-1,i])\n",
" \n",
" vn = 0.5*(v[j+1,i] + v[j,i])\n",
" vs = 0.5*(v[j,i] + v[j-1,i]) \n",
"\n",
" # convection = d(uv)/dx + d(vv)/dy\n",
" convection = - (ue*ve - uw*vw)/dx - (vn*vn - vs*vs)/dy\n",
" \n",
" # diffusion = d2u/dx2 + d2u/dy2\n",
" diffusion = ν*( (v[j,i+1] - 2.0*v[j,i] + v[j,i-1])/dx/dx + (v[j+1,i] - 2.0*v[j,i] + v[j-1,i])/dy/dy )\n",
" \n",
" vt[j,i] = v[j,i] + dt*(convection + diffusion) \n",
" # do pressure - prhs = 1/dt * div(uhat)\n",
" # we will only need to fill the interior points. This size is for convenient indexing\n",
" divut = np.zeros([ny+2,nx+2]) \n",
"# for i in range(1,nx+1):\n",
"# for j in range(1,ny+1):\n",
"# divutilde[j,i] = (ut[j,i+1] - ut[j,i])/dx + (vt[j+1,i] - vt[j,i])/dy\n",
" divut[1:-1,1:-1] = (ut[1:-1,2:] - ut[1:-1,1:-1])/dx + (vt[2:,1:-1] - vt[1:-1,1:-1])/dy\n",
"\n",
" prhs = 1.0/dt * divut\n",
" toc = time.clock()\n",
" momtime += (toc - tic)\n",
" \n",
" ###### Use the direct solver\n",
"# tic=time.clock() \n",
"# pressure_poisson2(p,prhs,dx,dy)\n",
"# toc=time.clock()\n",
"# solvertime += toc - tic\n",
"# psol.append(p)\n",
"\n",
" \n",
" \n",
" ###### Use the sparse linear solver\n",
" tic=time.clock()\n",
"# pt = scipy.sparse.linalg.spsolve(A1,prhs[1:-1,1:-1].ravel()) #theta=sc.linalg.solve_triangular(A,d)\n",
" pt,info = scipy.sparse.linalg.bicg(A1,prhs[1:-1,1:-1].ravel(),tol=1e-10) #theta=sc.linalg.solve_triangular(A,d)\n",
" toc=time.clock()\n",
" solvertime += toc - tic\n",
" p = np.zeros([ny+2,nx+2])\n",
" p[1:-1,1:-1] = pt.reshape([ny,nx])\n",
"\n",
" # time advance\n",
" u[1:-1,2:-1] = ut[1:-1,2:-1] - dt * (p[1:-1,2:-1] - p[1:-1,1:-2])/dx\n",
" v[2:-1,1:-1] = vt[2:-1,1:-1] - dt * (p[2:-1,1:-1] - p[1:-2,1:-1])/dy \n",
" \n",
" \n",
"# # Check mass residual\n",
"# divunp1 = np.zeros([ny+2,nx+2])\n",
"# for i in range(1,nx+1):\n",
"# for j in range(1,ny+1):\n",
"# divunp1[j,i] = (u[j,i+1] - u[j,i])/dx + (v[j+1,i] - v[j,i])/dy\n",
"# residual = np.linalg.norm(divunp1.ravel())\n",
"# if residual > 1e-6:\n",
"# print('Mass residual:',np.linalg.norm(divunp1.ravel()))\n",
" \n",
" # save new solutions\n",
"# usol.append(u)\n",
"# vsol.append(v)\n",
" \n",
" t += dt\n",
"tend = time.clock()\n",
"totaltime = tend - t0\n",
"print('total time', totaltime, 's')\n",
"print('time per timestep =',totaltime/nsteps, 's')\n",
"print('mom assembly time', totaltime - solvertime, 's')\n",
"print('solver time', solvertime, 's')\n",
"print('solver fraction =', np.ceil(solvertime/(tend - t0)*100),'%')"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"divu = (u[1:-1,2:] - u[1:-1,1:-1])/dx + (v[2:,1:-1] - v[1:-1,1:-1])/dy\n",
"plt.imshow(divu,origin='bottom')\n",
"plt.colorbar()"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig = plt.figure(figsize=[6,6],dpi=600)\n",
"ucc = 0.5*(u[1:-1,2:] + u[1:-1,1:-1])\n",
"vcc = 0.5*(v[2:,1:-1] + v[1:-1,1:-1])\n",
"speed = np.sqrt(ucc*ucc + vcc*vcc)\n",
"plt.contourf(speed)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"