Implement the Algorithm in Python¶
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
def convection (nt, nx, tmax, xmax, c):
"""
returns the velocity field and distance
"""
# Increments in space and time
dt = tmax/(nt-1)
dx = xmax/(nx-1)
# Initialization of data sructures
u = np.zeros((nx, nt))
x = np.zeros(nx)
# BCs
u[0, :] = u[nx-1, :] = 1
#ICs - 1st way
for i in range(0, nx-1):
if (i > (nx-1)/4 and i < (nx-1)/2):
u[i, 0] = 2
else:
u[i, 0] = 1
# Loop
for n in range(0, nt-1):
for i in range(1, nx-1):
u[i, n+1] = u[i, n]-c*(dt/dx)*(u[i, n]-u[i-1, n])
# X loop
for i in range(0, nx):
x[i] = i*dx
return u, x
def plot(u, x, nt, title):
"""
Plots the Velocity Field Results in 1D
"""
plt.figure()
plt.show()
for i in range(0, nt, 10):
plt.plot(x, u[:,i], 'r')
plt.xlabel('x [m]')
plt.ylabel('u [m/s]')
plt.ylim([0, 2.2])
plt.title(title)
plt.grid(True, linestyle='-.', linewidth='0.5', color='black')
Run¶
u,x = convection(20, 41, 2, 2.0, 0.5)
plot(u,x, 20,'Figure 1: c=0.5m/s, nt=20, nx=41, tmax=0.5s')
u,x = convection(70, 41, 3, 2.0, 0.5)
plot(u,x, 20,'Figure 2: c=0.5m/s, nt=70, nx=41, tmax=0.5s')
Instabilities?¶
To answer that question, we have to think a little bit about what we're actually implementing in code.
In each iteration of our time loop, we use the existing data about our wave to estimate the speed of the wave in the subsequent time step. Initially, the increase in the number of grid points returned more accurate answers. There was less numerical diffusion and the square wave looked much more like a square wave than it did in the initial runs.
Each iteration of our time loop covers a time-step of length Δt, which we have been defining as 0.025
During this iteration, we evaluate the speed of the wave at each of the x points we've created. In the last plot, something has clearly gone wrong.
What has happened is that over the time period Δt, the wave is travelling a distance which is greater than dx. The length dx of each grid box is related to the number of total points $n_x$, so stability can be enforced if the $\delta{t}$ step size is calculated with respect to the size of dx.
$$ \sigma = \frac{u\Delta{t}}{\Delta{x}}\le\sigma_{max} $$where $u$ is the speed of the wave; $\sigma$ is called the Courant number and the value of $\sigma_{max}$ that will ensure stability depends on the discretization used.
In the following version of the code, I will use the CFL number to calculate the appropriate time-step dt depending on the size of dx.
def convection (nt, nx, tmax, xmax, c):
"""
returns the velocity field and distance
"""
# Increments in space and time
#dt = tmax/(nt-1)
dx = xmax/(nx-1)
sigma = 0.5
dt = sigma * dx
# Initialization of data sructures
u = np.zeros((nx, nt))
x = np.zeros(nx)
# BCs
u[0, :] = u[nx-1, :] = 1
#ICs - 1st way
for i in range(0, nx-1):
if (i > (nx-1)/4 and i < (nx-1)/2):
u[i, 0] = 2
else:
u[i, 0] = 1
# Loop
for n in range(0, nt-1):
for i in range(1, nx-1):
u[i, n+1] = u[i, n]-c*(dt/dx)*(u[i, n]-u[i-1, n])
# X loop
for i in range(0, nx):
x[i] = i*dx
return u, x
def plot(u, x, nt, title):
"""
Plots the Velocity Field Results in 1D
"""
plt.figure()
#plt.show()
for i in range(0, nt, 10):
plt.plot(x, u[:,i], 'g')
plt.xlabel('x [m]')
plt.ylabel('u [m/s]')
plt.ylim([0, 2.2])
plt.title(title)
plt.grid(True, linestyle='-.', linewidth='0.5', color='black')
plt.savefig('plot.png')
u,x = convection(20, 41, 2, 2.0, 0.5)
plot(u,x, 20,'Figure 1: c=0.5m/s, nt=20, nx=41, tmax=0.5s')
u,x = convection(20, 81, 2, 2.0, 0.5)
plot(u,x, 20,'Figure 1: c=0.5m/s, nt=20, nx=61, tmax=0.5s')
u,x = convection(20, 200, 2, 2.0, 0.5)
plot(u,x, 20,'Figure 1: c=0.5m/s, nt=20, nx=61, tmax=0.5s')
Notice that as the number of points nx increases, the wave convects a shorter and shorter distance. The number of time iterations we have advanced the solution at is held constant at $nt = 20$, but depending on the value of $nx$ and the corresponding values of $dx$ and $dt$, a shorter time window is being examined overall.