Written Notes

Examples

1D Simple Advection Equation

# forward euler method for time

# second order central method for space

import numpy as np

import matplotlib.pyplot as plt

from random import randrange

# Spatial Parameters

L = 10

n = 301

dx = L/(n-1)

# Initial Conditions

x = np.linspace(0, L, n)

u0 = np.exp(-(x-2)**2/0.1)

# Temporal Parameters

dt = 1e-3

tend = 6

t = 0.0

# Implementing the Formula

c = -1.0 # wave speed [m/s]

cfl = c * dt/2.0/dx

# Solving

sol = [u0]

while t <= tend:

un = sol[-1]

unew = []

for i in range(len(un)-1):

unew.append(un[i] + cfl * (un[i+1] - un[i-1]))

unew.append(0)

sol.append(unew)

t += dt

# Plotting the solution in time

counter = 0

fig = plt.figure('Simple Advection Equation Solution')

ax = fig.add_subplot(111, projection='3d')

colors = ['k', 'b', 'g', 'm', 'r', 'c', 'w']

for i in range(len(sol)-1):

if i % 250 == 0:

z = []

for j in range(len(x)):

z.append(dt * counter)

counter += 1

ax.plot3D(x, z, sol[i], colors[randrange(len(colors)-1)])

ax.set_xlabel('x-direction')

ax.set_ylabel('time Label')

ax.set_zlabel('y-direction')

ax.view_init(elev=20, azim=270)

plt.show()

*As is clear by the plot, the 1D wave is advecting (traveling) in space without diffusion (growing or decaying).

1D Linear Advection-Diffusion Equation

# forward euler method for time

# second order central method for space

import numpy as np

import matplotlib.pyplot as plt

from random import randrange

n = 81

L = 1.0

x = np.linspace(0, L, n)

dx = L/(n-1)

w = 2.0*np.pi

mu = 0.05

c = 1.0

u0 = np.sin(w*x)

dt = 0.001

tend = 1.0

t = 0

cfl = c*dt/2.0/dx

nfl = mu*dt/dx/dx

sol = []

sol.append(u0)

while t < tend:

counter = len(sol)

un = sol[counter-1].copy()

unew = sol[counter-1].copy()

unew[1:-1] = un[1:-1] - cfl * (un[2:] - un[:-2]) + nfl * (un[2:] - 2.0*un[1:-1] + un[:-2])

unew[-1] = un[-1] - cfl * (un[1] - un[-2]) + nfl * (un[-2] - 2.0*un[-1] + un[1])

unew[0] = unew[-1]

sol.append(unew)

t += dt

# Plotting the solution in time

colors = ['k', 'b', 'g', 'm', 'r', 'c', 'w']

for k in range(0, 361):

counter = 0

fig = plt.figure('Simple Advection Equation Solution')

ax = fig.add_subplot(111, projection='3d')

for i in range(len(sol) - 1):

if i % 7 == 0:

z = []

for j in range(len(x)):

z.append(dt * counter)

counter += 1

ax.plot3D(x, z, sol[i], colors[randrange(len(colors) - 1)])

ax.set_xticks([])

ax.set_yticks([])

ax.set_zticks([])

ax.set_ylabel('time')

ax.view_init(elev=20, azim=k)

filename = ['Image/Azim_', str(k), '.png']

plt.savefig("".join(filename), dpi=1000)

plt.close()

*The function is following the expected behavior for an advection-diffusion PDE

2D Linear Advection-Diffusion Equation

# forward euler method for time

# second order central method for space

import numpy as np

import matplotlib.pyplot as plt

nx = 21

ny = 21

Lx = 1.0

Ly = 1.0

x = np.linspace(0, Lx, nx)

y = np.linspace(0, Ly, ny)

dx = Lx/(nx-1)

dy = Ly/(ny-1)

#create a grid of coordinates

xx, yy = np.meshgrid(x, y)

w0 = 2.0*np.pi

u0 = lambda x, y: np.sin(w0*x) * np.sin(w0*y)

# plot the initial condition

# plt.contourf(xx, yy, u0(xx, yy))

# plt.show()

cx = 1.0

cy = 1.0

dt = 0.001

tend = 2.0

t = 0

cflx = cx * dt/dy

cfly = cy * dt/dy

# setup the initial condition

sol = []

u = np.zeros([nx+2, ny+2])

u[1:-1, 1:-1] = u0(xx, yy)

# set periodic boundaries for the initial condition

u[:, 0] = u[:, -3] # x-minus face

u[:, -1] = u[:, 2] # x-plus face

u[0, :] = u[-3, :] # y-minus face

u[-1, :] = u[2, :] # y-plus face

# solving

sol.append(u)

while t < tend:

un = sol[-1]

unew = un.copy()

unew[1:-1, 1:-1] = un[1:-1, 1:-1] - cflx*(un[1:-1, 1:-1] - un[1:-1, :-2]) - cfly * (un[1:-1, 1:-1] - un[:-2, 1:-1])

unew[:, 0] = unew[:, -3] # x-minus face

unew[:, -1] = unew[:, 2] # x-plus face

unew[0, :] = unew[-3, :] # y-minus face

unew[-1, :] = unew[2, :] # y-plus face

sol.append(unew)

t += dt

for i in range(0, 2003):

filename = ['Image/image_', str(i), '.png']

solution = sol[i]

plt.contourf(xx, yy, solution[1:-1, 1:-1], cmap='magma')

cbar = plt.colorbar()

plt.clim(-1, 1)

cbar.set_ticks(np.linspace(-1, 1, 5))

plt.savefig("".join(filename))

plt.close()

*The figure shows that the function is moving and decaying in time, this is the expected behavior for an advection-diffusion PDE.

2D Lid-Driven Cavity flow

# Successive Over Relaxation Method (FTCS)

import numpy as np

import matplotlib.pyplot as plt

# Defining Calculation Domain

nx = 53

ny = 53

lx = 1

ly = 1

dx = lx/(nx-1)

dy = ly/(ny-1)

# apply boundary conditions on omega0

Ut = 10.0 # velocity @ top wall

psi0 = np.zeros([nx, ny])

omega0 = np.zeros([nx, ny])

omega0[:, 0] = 2*(-psi0[:, 1])/dx/dx # left wall

omega0[:, -1] = 2*(-psi0[:, -2])/dx/dx # right wall

omega0[-1, :] = 2*(-psi0[-2, :])/dy/dy - 2*Ut/dy # top wall

omega0[0, :] = 2*(-psi0[1, :])/dy/dy # bottom wall

# Simulation Parameters

beta = 1.5

tol = 1e-3

maxIt = 30

t = 0

nu = 0.05

dt = min(0.25*dx*dx/nu, 4*nu/Ut/Ut) # from Von-Neumann analysis criteria

tend = 1000 * dt

print('dt =', dt, '[sec]')

print('Re =', Ut * lx/nu)

# Initializing solutions list

psisol = []

psisol.append(psi0)

omegasol = []

omegasol.append(omega0)

# Begin solution procedure

time_var = []

while t < tend:

# start by solving the Poisson equation for the stream-function

print(t)

it = 0

err = 1e5

omegan = omegasol[-1]

psin = psisol[-1].copy()

while err > tol and it < maxIt:

psik = np.zeros_like(psin)

psik[1:-1, 1:-1] = psin[1:-1, 1:-1]

# Loop over interior points

for i in range(1, nx-1):

for j in range(1, ny-1):

rhs = (dx*dy)**2*omegan[j, i] + dy**2*(psin[j, i+1]+psin[j, i-1]) + dx**2*(psin[j+1, i]+psin[j-1, i])

rhs *= beta/2/(dx**2 + dy**2)

psin[j, i] = rhs + (1 - beta)*psin[j, i]

error = np.linalg.norm(psin.ravel() - psik.ravel()) # ravel basically unravels a matrix into a vector

it += 1

psisol.append(psin)

# work on omega

omega = np.zeros_like(omegan)

Cx = -(psin[2:, 1:-1] - psin[:-2, 1:-1]) / 2 / dy * (omegan[1:-1, 2:] - omegan[1:-1, :-2]) / 2 / dx # x Convection

Cy = (omegan[2:, 1:-1] - omegan[:-2, 1:-1]) / 2 / dy * (psin[1:-1, 2:] - psin[1:-1, :-2]) / 2 / dx # y Convection

Dx = (omegan[1:-1, 2:] - 2*omegan[1:-1, 1:-1] + omegan[1:-1, :-2])/dx/dx # Diffusion in x

Dy = (omegan[2:, 1:-1] - 2 * omegan[1:-1, 1:-1] + omegan[:-2, 1:-1])/dy/dy # Diffusion in y

rhs = Cx + Cy + nu*(Dx+Dy)

omega[1:-1, 1:-1] = omegan[1:-1, 1:-1] + dt * rhs

# apply boundary conditions on vorticity

omega[:, 0] = 2 * (-psin[:, 1]) / dx / dx # left wall

omega[:, -1] = 2 * (-psin[:, -2]) / dx / dx # right wall

omega[-1, :] = 2 * (-psin[-2, :]) / dy / dy - 2 * Ut / dy # top wall

omega[0, :] = 2 * (-psin[1, :]) / dy / dy # bottom wall

omegasol.append(omega)

time_var.append(t)

t += dt

# Plotting

x = np.linspace(0, 1, nx)

y = np.linspace(0, 1, ny)

xx, yy = np.meshgrid(x, y)

psi = psisol[-1]

u = (psi[2:, 1:-1] - psi[:-2, 1:-1]) / 2 / dy

v = -(psi[1:-1, 2:] - psi[1:-1, :-2]) / 2 / dx

fig = plt.figure()

ax = fig.add_subplot(111)

ax.contourf(xx[1:-1, 1:-1], yy[1:-1, 1:-1], np.sqrt(u * u + v * v), levels=100, cmap='Spectral')

ax.streamplot(xx[1:-1, 1:-1], yy[1:-1, 1:-1], u, v, color=abs(u * u + v * v), cmap='bone', density=3, arrowsize=0.5, linewidth=0.5, arrowstyle='Fancy')

ax.set_xlim([xx[0, 1], xx[0, -2]])

ax.set_aspect(1)

plt.show()

*Be mindful that this code takes a few minutes to run, for purposes of shortening the runtime I would recommend setting nx and ny = 23.

Lid-Driven Cavity flow (2D) #2

# Projection Method (FTCS)

import numpy as np

import matplotlib.pyplot as plt

def ddx(f, dx):

result = np.zeros_like(f)

result[1:-1, 1:-1] = (f[1:-1, 2:] - f[1:-1, :-2]) / 2.0 / dx

return result

def ddy(f, dy):

result = np.zeros_like(f)

result[1:-1, 1:-1] = (f[2:, 1:-1] - f[:-2, 1:-1]) / 2.0 / dy

return result

def laplacian(f, dx, dy):

result = np.zeros_like(f)

result[1:-1, 1:-1] = (f[1:-1, 2:] - 2.0 * f[1:-1, 1:-1] + f[1:-1, :-2]) / dx / dx \

+ (f[2:, 1:-1] - 2.0 * f[1:-1, 1:-1] + f[:-2, 1:-1]) / dy / dy

return result

def div(u, v, dx, dy):

return ddx(u, dx) + ddy(v, dy)

# this is similar to the previous SOR solver

def pressure_poisson(p, dx, dy, b):

pn = np.empty_like(p)

it = 0

err = 1e5

tol = 1e-3

maxit = 50

while it < maxit and err > tol:

pn = p.copy()

p[1:-1, 1:-1] = (((pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy ** 2 +

(pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx ** 2) /

(2 * (dx ** 2 + dy ** 2)) -

dx ** 2 * dy ** 2 / (2 * (dx ** 2 + dy ** 2)) *

b[1:-1, 1:-1])

p[:, -1] = p[:, -2] # dp/dy = 0 at right wall

p[0, :] = p[1, :] # dp/dy = 0 at y = 0

p[:, 0] = p[:, 1] # dp/dx = 0 at x = 0

p[-1, :] = 0 # p = 0 at the top wall

err = np.linalg.norm(p - pn, 2)

it += 1

return p, err

# Defining Calculation Domain

nx = 52

ny = 52

lx = 1

ly = 1

dx = lx/(nx-1)

dy = ly/(ny-1)

nu = 0.025 # kinematic viscosity

Ut = 10.0 # velocity

dt = min(0.25*dx*dx/nu, 4*nu/Ut/Ut)/2 # not exact, but a good approximation

t = 0

tend = 2500*dt

print("Re:", Ut*lx/nu)

u = np.zeros([ny, nx])

v = np.zeros([ny, nx])

uh = np.zeros([ny, nx])

vh = np.zeros([ny, nx])

p = np.zeros([ny, nx])

counter = 0

while t < tend:

# set boundary conditions

# bottom wall

u[0, :] = 0.0

v[0, :] = 0.0

# top wall

u[-1, :] = Ut

v[-1, :] = 0.0

# left wall

u[:, 0] = 0.0

v[:, 0] = 0.0

# right wall

u[:, -1] = 0.0

v[:, -1] = 0.0

# do the x-momentum RHS

# u rhs: - d(uu)/dx - d(vu)/dy + ν d2(u)

uRHS = - ddx(u * u, dx) - ddy(v * u, dy) + nu * laplacian(u, dx, dy)

# v rhs: - d(uv)/dx - d(vv)/dy + ν d2(v)

vRHS = - ddx(u * v, dx) - ddy(v * v, dy) + nu * laplacian(v, dx, dy)

uh = u + dt * uRHS

vh = v + dt * vRHS

# next compute the pressure RHS: prhs = div(un)/dt + div( [urhs, vrhs])

prhs = div(uh, vh, dx, dy) / dt

p, err = pressure_poisson(p, dx, dy, prhs)

# finally compute the true velocities

# u_{n+1} = uh - dt*dpdx

u = uh - dt * ddx(p, dx)

v = vh - dt * ddy(p, dy)

if counter % 5 == 0:

filename = ['Image/image_', str(counter), '.png']

vel = np.sqrt(uh ** 2 + vh ** 2)

plt.contourf(vel, levels=100)

plt.quiver(u, v, headwidth=3, headlength=3, scale=50)

plt.savefig("".join(filename))

plt.close()

print(t / dt)

counter += 1

t += dt

2D Lid Driven Cavity - Finite Volume Method (FVM) with Sparse Solver

import time

import numpy as np

import matplotlib.pyplot as plt

import scipy.sparse

import scipy.sparse.linalg

start = time.time()

# Define Computational Domain and Simulation Parameters

nx = 25

ny = 25

nu = 0.1

lx = 1.0

ly = 1.0

dx = lx / nx

dy = ly / ny

t = 0.0

nsteps = 7500

# cell centered coordinates

xx = np.linspace(dx / 2.0, lx - dx / 2.0, nx, endpoint=True)

yy = np.linspace(dy / 2.0, ly - dy / 2.0, ny, endpoint=True)

xcc, ycc = np.meshgrid(xx, yy)

# x-staggered coordinates

xxs = np.linspace(0, lx, nx + 1, endpoint=True)

xu, yu = np.meshgrid(xxs, yy)

# y-staggered coordinates

yys = np.linspace(0, ly, ny + 1, endpoint=True)

xv, yv = np.meshgrid(xx, yys)

Ut = 75.0

Ub = 0.0

Vl = 0.0

Vr = 0.0

print('Reynolds Number:', Ut * lx / nu)

dt = 0.75 * min(0.25 * dx * dx / nu, 4.0 * nu / Ut / Ut)

print('dt=', dt)

# initialize velocities (stagger everything in the negative direction)

u = np.zeros([ny + 2, nx + 2]) # include ghost cells

# same thing for the y-velocity component

v = np.zeros([ny + 2, nx + 2]) # include ghost cells

ut = np.zeros_like(u)

vt = np.zeros_like(u)

# initialize the pressure

p = np.zeros([nx + 2, ny + 2]) # include ghost cells

# a bunch of lists for animation purposes

usol = [u]

vsol = [v]

psol = [p]

# USE sparse solver

# build pressure coefficient matrix

Ae = 1.0 / dx / dx * np.ones([ny, nx])

As = 1.0 / dy / dy * np.ones([ny, nx])

An = 1.0 / dy / dy * np.ones([ny, nx])

Aw = 1.0 / dx / dx * np.ones([ny, nx])

Aw[:, 0] = 0.0 # set left wall coefficients

Ae[:, -1] = 0.0 # set right wall coefficients

An[-1, :] = 0.0 # set top wall coefficients

As[0, :] = 0.0 # set bottom wall coefficients

Ap = -(Aw + Ae + An + As)

n = nx * ny

d0 = Ap.reshape(n)

de = Ae.reshape(n)[:-1]

dw = Aw.reshape(n)[1:]

ds = As.reshape(n)[nx:]

dn = An.reshape(n)[:-nx]

A1 = scipy.sparse.diags([d0, de, dw, dn, ds], [0, 1, -1, nx, -nx], format='csr')

for n in range(0, nsteps):

u[1:-1, 1] = 0.0 # left wall

u[1:-1, -1] = 0.0 # right wall

u[-1, 1:] = 2.0 * Ut - u[-2, 1:] # top wall

u[0, 1:] = 2.0 * Ub - u[1, 1:] # bottom wall

v[1:, 0] = 2.0 * Vl - v[1:, 1] # left wall

v[1:, -1] = 2.0 * Vr - v[1:, -2] # right wall

v[1, 1:-1] = 0.0 # bottom wall

v[-1, 1:-1] = 0.0 # top wall

# x-momentum (interior points)

for i in range(2, nx + 1):

for j in range(1, ny + 1):

ue = 0.5 * (u[j, i + 1] + u[j, i])

uw = 0.5 * (u[j, i] + u[j, i - 1])

un = 0.5 * (u[j + 1, i] + u[j, i])

us = 0.5 * (u[j, i] + u[j - 1, i])

vn = 0.5 * (v[j + 1, i] + v[j + 1, i - 1])

vs = 0.5 * (v[j, i] + v[j, i - 1])

# convection

convection = - (ue * ue - uw * uw) / dx - (un * vn - us * vs) / dy

# diffusion

diffusion = nu * ((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)

ut[j, i] = u[j, i] + dt * (convection + diffusion)

# y-momentum (interior points)

for i in range(1, nx + 1):

for j in range(2, ny + 1):

ve = 0.5 * (v[j, i + 1] + v[j, i])

vw = 0.5 * (v[j, i] + v[j, i - 1])

ue = 0.5 * (u[j, i + 1] + u[j - 1, i + 1])

uw = 0.5 * (u[j, i] + u[j - 1, i])

vn = 0.5 * (v[j + 1, i] + v[j, i])

vs = 0.5 * (v[j, i] + v[j - 1, i])

# convection = d(uv)/dx + d(vv)/dy

convection = - (ue * ve - uw * vw) / dx - (vn * vn - vs * vs) / dy

# diffusion = d2u/dx2 + d2u/dy2

diffusion = nu * ((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)

vt[j, i] = v[j, i] + dt * (convection + diffusion)

# pressure (interior points)

divut = np.zeros([ny + 2, nx + 2])

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

prhs = 1.0 / dt * divut

# Use the sparse linear solver

pt, info = scipy.sparse.linalg.bicg(A1, prhs[1:-1, 1:-1].ravel(), tol=1e-10)

p = np.zeros([ny + 2, nx + 2])

p[1:-1, 1:-1] = pt.reshape([ny, nx])

# time advance

u[1:-1, 2:-1] = ut[1:-1, 2:-1] - dt * (p[1:-1, 2:-1] - p[1:-1, 1:-2]) / dx

v[2:-1, 1:-1] = vt[2:-1, 1:-1] - dt * (p[2:-1, 1:-1] - p[1:-2, 1:-1]) / dy

# save new solutions

usol.append(u)

vsol.append(v)

t += dt

print("Physical Time:", t, "[sec]")

print("Calculation Time:", time.time()-start, "[sec]")

# Plotting

plt.figure(figsize=[5, 5])

ucc = 0.5 * (u[1:-1, 2:] + u[1:-1, 1:-1])

vcc = 0.5 * (v[2:, 1:-1] + v[1:-1, 1:-1])

x = np.linspace(0, lx, nx)

y = np.linspace(0, ly, ny)

xx, yy = np.meshgrid(x, y)

nn = 1

plt.xlim([xx[0, 0], xx[0, -1]])

plt.ylim([yy[0, 0], yy[-1, 0]])

plt.streamplot(xx, yy, ucc, vcc, color=np.sqrt(ucc * ucc + vcc * vcc), density=5, linewidth=.5, arrowsize=0.01,

cmap='bone')

plt.contourf(xx, yy, np.sqrt(ucc * ucc + vcc * vcc), levels=100, cmap='Spectral')

plt.show()