4 Advection-Dominated Equations
Wave (Chapter Chapter 2) and diffusion (Chapter Chapter 3) equations are solved reliably by finite difference methods. As soon as we add a first-order derivative in space, representing advective transport (also known as convective transport), the numerics gets more complicated and intuitively attractive methods no longer work well. We shall show how and why such methods fail and provide remedies. The present chapter builds on basic knowledge about finite difference methods for diffusion and wave equations, including the analysis by Fourier components, truncation error analysis (Appendix Chapter 7), and compact difference notation.
It is common to refer to movement of a fluid as convection, while advection is the transport of some material dissolved or suspended in the fluid. We shall mostly choose the word advection here, but both terms are in heavy use, and for mass transport of a substance the PDE has an advection term, while the similar term for the heat equation is a convection term.
Much more comprehensive discussion of dispersion analysis for advection problems can be found in the book by Duran (Duran 2010). This is a an excellent resource for further studies on the topic of advection PDEs, with emphasis on generalizations to real geophysical problems. The book by Fletcher (Fletcher 2013) also has a good overview of methods for advection and convection problems.
4.1 1D linear advection equations with constant velocity
We consider the pure advection model
\[ \frac{\partial u}{\partial t} + v\frac{\partial u}{\partial x} = 0,\quad x\in (0,L),\ t\in (0,T], \tag{4.1}\] \[ u(x,0) = I(x), \quad x\in (0,L), \tag{4.2}\] \[ u(0,t) = U_0, \quad t\in (0,T]. \tag{4.3}\] In (4.3), \(v\) is a given parameter, typically reflecting the transport velocity of a quantity \(u\) with a flow. There is only one boundary condition (4.2) since the spatial derivative is only first order in the PDE (4.3). The information at \(x=0\) and the initial condition get transported in the positive \(x\) direction if \(v>0\) through the domain.
It is easiest to find the solution of (4.3) if we remove the boundary condition and consider a process on the infinite domain \((-\infty, \infty)\). The solution is simply \[ u(x,t) = I(x-vt)\tp \tag{4.4}\] This is also the solution we expect locally in a finite domain before boundary conditions have reflected or modified the wave.
A particular feature of the solution (4.4) is that \[ u(x_i, t_{n+1}) = u(x_{i-1}, t_n), \tag{4.5}\] if \(x_i=i\Delta x\) and \(t_n=n\Delta t\) are points in a uniform mesh. We see this relation from
\[\begin{align} u(i\Delta x, (n+1)\Delta t) &= I(i\Delta x - v(n+1)\Delta t) \nonumber \\ &= I((i-1)\Delta x - vn\Delta t - v\Delta t + \Delta x) \nonumber \\ &= I((i-1)\Delta x - vn\Delta t) \nonumber \\ &= u((i-1)\Delta x, n\Delta t), \nonumber \end{align}\] provided \(v = \Delta x/\Delta t\). So, whenever we see a scheme that collapses to \[ u^{n+1}**i = u**{i-1}^n, \tag{4.6}\] for the PDE in question, we have in fact a scheme that reproduces the analytical solution, and many of the schemes to be presented possess this nice property!
Finally, we add that a discussion of appropriate boundary conditions for the advection PDE in multiple dimensions is a challenging topic beyond the scope of this text.
4.2 Simplest scheme: forward in time, centered in space
4.2.1 Method
A first attempt to solve a PDE like (4.3) will normally be to look for a time-discretization scheme that is explicit so we avoid solving systems of linear equations. In space, we anticipate that centered differences are most accurate and therefore best. These two arguments lead us to a Forward Euler scheme in time and centered differences in space: \[ [D_t^+ u + vD_{2x} u = 0]^n_i \] Written out, we see that this expression implies that \[ u^{n+1} = u^n - \half C (u^n_{i+1}-u_{i-1}^n), \] with \(C\) as the Courant number \[ C = \frac{v\Delta t}{\Delta x}\tp \] ### Implementation
A solver function for our scheme goes as follows.
import numpy as np
def solver_FECS(I, U0, v, L, dt, C, T, user_action=None):
Nt = int(round(T / float(dt)))
t = np.linspace(0, Nt * dt, Nt + 1) # Mesh points in time
dx = v * dt / C
Nx = int(round(L / dx))
x = np.linspace(0, L, Nx + 1) # Mesh points in space
dx = x[1] - x[0]
dt = t[1] - t[0]
C = v * dt / dx
u = np.zeros(Nx + 1)
u_n = np.zeros(Nx + 1)
for i in range(0, Nx + 1):
u_n[i] = I(x[i])
if user_action is not None:
user_action(u_n, x, t, 0)
for n in range(0, Nt):
for i in range(1, Nx):
u[i] = u_n[i] - 0.5 * C * (u_n[i + 1] - u_n[i - 1])
u[0] = U0
if user_action is not None:
user_action(u, x, t, n + 1)
u_n, u = u, u_n
def solver(I, U0, v, L, dt, C, T, user_action=None, scheme="FE", periodic_bc=True):
Nt = int(round(T / float(dt)))
t = np.linspace(0, Nt * dt, Nt + 1) # Mesh points in time
dx = v * dt / C
Nx = int(round(L / dx))
x = np.linspace(0, L, Nx + 1) # Mesh points in space
dx = x[1] - x[0]
dt = t[1] - t[0]
C = v * dt / dx
print("dt=%g, dx=%g, Nx=%d, C=%g" % (dt, dx, Nx, C))
u = np.zeros(Nx + 1)
u_n = np.zeros(Nx + 1)
u_nm1 = np.zeros(Nx + 1)
integral = np.zeros(Nt + 1)
for i in range(0, Nx + 1):
u_n[i] = I(x[i])
u[0] = U0
integral[0] = dx * (0.5 * u_n[0] + 0.5 * u_n[Nx] + np.sum(u_n[1:-1]))
if user_action is not None:
user_action(u_n, x, t, 0)
for n in range(0, Nt):
if scheme == "FE":
if periodic_bc:
i = 0
u[i] = u_n[i] - 0.5 * C * (u_n[i + 1] - u_n[Nx])
u[Nx] = u[0]
for i in range(1, Nx):
u[i] = u_n[i] - 0.5 * C * (u_n[i + 1] - u_n[i - 1])
elif scheme == "LF":
if n == 0:
if periodic_bc:
i = 0
u_n[i] = u_n[Nx]
for i in range(1, Nx + 1):
u[i] = u_n[i] - C * (u_n[i] - u_n[i - 1])
else:
if periodic_bc:
i = 0
u[i] = u_nm1[i] - C * (u_n[i + 1] - u_n[Nx - 1])
for i in range(1, Nx):
u[i] = u_nm1[i] - C * (u_n[i + 1] - u_n[i - 1])
if periodic_bc:
u[Nx] = u[0]
elif scheme == "UP":
if periodic_bc:
u_n[0] = u_n[Nx]
for i in range(1, Nx + 1):
u[i] = u_n[i] - C * (u_n[i] - u_n[i - 1])
elif scheme == "LW":
if periodic_bc:
i = 0
u[i] = (
u_n[i]
- 0.5 * C * (u_n[i + 1] - u_n[Nx - 1])
+ 0.5 * C * (u_n[i + 1] - 2 * u_n[i] + u_n[Nx - 1])
)
for i in range(1, Nx):
u[i] = (
u_n[i]
- 0.5 * C * (u_n[i + 1] - u_n[i - 1])
+ 0.5 * C * (u_n[i + 1] - 2 * u_n[i] + u_n[i - 1])
)
if periodic_bc:
u[Nx] = u[0]
else:
raise ValueError('scheme="%s" not implemented' % scheme)
if not periodic_bc:
u[0] = U0
integral[n + 1] = dx * (0.5 * u[0] + 0.5 * u[Nx] + np.sum(u[1:-1]))
if user_action is not None:
user_action(u, x, t, n + 1)
u_nm1, u_n, u = u_n, u, u_nm1
print("I:", integral[n + 1])
return integral
def run_FECS(case):
"""Special function for the FECS case."""
if case == "gaussian":
def I(x):
return np.exp(-0.5 * ((x - L / 10) / sigma) ** 2)
elif case == "cosinehat":
def I(x):
return np.cos(np.pi * 5 / L * (x - L / 10)) if x < L / 5 else 0
L = 1.0
sigma = 0.02
legends = []
def plot(u, x, t, n):
"""Animate and plot every m steps in the same figure."""
plt.figure(1)
if n == 0:
lines = plot(x, u)
else:
lines[0].set_ydata(u)
plt.draw()
plt.figure(2)
m = 40
if n % m != 0:
return
print(
"t=%g, n=%d, u in [%g, %g] w/%d points" % (t[n], n, u.min(), u.max(), x.size)
)
if np.abs(u).max() > 3: # Instability?
return
plt.plot(x, u)
legends.append("t=%g" % t[n])
plt.ion()
U0 = 0
dt = 0.001
C = 1
T = 1
solver(I=I, U0=U0, v=1.0, L=L, dt=dt, C=C, T=T, user_action=plot)
plt.legend(legends, loc="lower left")
plt.savefig("tmp.png")
plt.savefig("tmp.pdf")
plt.axis([0, L, -0.75, 1.1])
plt.show()
def run(scheme="UP", case="gaussian", C=1, dt=0.01):
"""General admin routine for explicit and implicit solvers."""
if case == "gaussian":
def I(x):
return np.exp(-0.5 * ((x - L / 10) / sigma) ** 2)
elif case == "cosinehat":
def I(x):
return np.cos(np.pi * 5 / L * (x - L / 10)) if 0 < x < L / 5 else 0
L = 1.0
sigma = 0.02
global lines # needs to be saved between calls to plot
def plot(u, x, t, n):
"""Plot t=0 and t=0.6 in the same figure."""
plt.figure(1)
global lines
if n == 0:
lines = plt.plot(x, u)
plt.axis([x[0], x[-1], -0.5, 1.5])
plt.xlabel("x")
plt.ylabel("u")
plt.axes().set_aspect(0.15)
plt.savefig("tmp_%04d.png" % n)
plt.savefig("tmp_%04d.pdf" % n)
else:
lines[0].set_ydata(u)
plt.axis([x[0], x[-1], -0.5, 1.5])
plt.title("C=%g, dt=%g, dx=%g" % (C, t[1] - t[0], x[1] - x[0]))
plt.legend(["t=%.3f" % t[n]])
plt.xlabel("x")
plt.ylabel("u")
plt.draw()
plt.savefig("tmp_%04d.png" % n)
plt.figure(2)
eps = 1e-14
if abs(t[n] - 0.6) > eps and abs(t[n] - 0) > eps:
return
print(
"t=%g, n=%d, u in [%g, %g] w/%d points" % (t[n], n, u.min(), u.max(), x.size)
)
if np.abs(u).max() > 3: # Instability?
return
plt.plot(x, u)
plt.draw()
if n > 0:
y = [I(x_ - v * t[n]) for x_ in x]
plt.plot(x, y, "k--")
if abs(t[n] - 0.6) < eps:
filename = ("tmp_%s_dt%s_C%s" % (scheme, t[1] - t[0], C)).replace(".", "")
np.savez(filename, x=x, u=u, u_e=y)
plt.ion()
U0 = 0
T = 0.7
v = 1
codecs = dict(flv="flv", mp4="libx264", webm="libvpx", ogg="libtheora")
import glob
import os
for name in glob.glob("tmp_*.png"):
os.remove(name)
for ext in codecs:
name = "movie.%s" % ext
if os.path.isfile(name):
os.remove(name)
if scheme == "CN":
integral = solver_theta(I, v, L, dt, C, T, user_action=plot, FE=False)
elif scheme == "BE":
integral = solver_theta(I, v, L, dt, C, T, theta=1, user_action=plot)
else:
integral = solver(
I=I, U0=U0, v=v, L=L, dt=dt, C=C, T=T, scheme=scheme, user_action=plot
)
plt.figure(2)
plt.axis([0, L, -0.5, 1.1])
plt.xlabel("$x$")
plt.ylabel("$u$")
plt.axes().set_aspect(0.5) # no effect
plt.savefig("tmp1.png")
plt.savefig("tmp1.pdf")
plt.show()
for codec in codecs:
cmd = "ffmpeg -i tmp_%%04d.png -r 25 -vcodec %s movie.%s" % (codecs[codec], codec)
os.system(cmd)
print("Integral of u:", integral.max(), integral.min())
def solver_theta(I, v, L, dt, C, T, theta=0.5, user_action=None, FE=False):
"""
Full solver for the model problem using the theta-rule
difference approximation in time (no restriction on F,
i.e., the time step when theta >= 0.5).
Vectorized implementation and sparse (tridiagonal)
coefficient matrix.
"""
import time
t0 = time.perf_counter() # for measuring the CPU time
Nt = int(round(T / float(dt)))
t = np.linspace(0, Nt * dt, Nt + 1) # Mesh points in time
dx = v * dt / C
Nx = int(round(L / dx))
x = np.linspace(0, L, Nx + 1) # Mesh points in space
dx = x[1] - x[0]
dt = t[1] - t[0]
C = v * dt / dx
print("dt=%g, dx=%g, Nx=%d, C=%g" % (dt, dx, Nx, C))
u = np.zeros(Nx + 1)
u_n = np.zeros(Nx + 1)
u_nm1 = np.zeros(Nx + 1)
integral = np.zeros(Nt + 1)
for i in range(0, Nx + 1):
u_n[i] = I(x[i])
integral[0] = dx * (0.5 * u_n[0] + 0.5 * u_n[Nx] + np.sum(u_n[1:-1]))
if user_action is not None:
user_action(u_n, x, t, 0)
diagonal = np.zeros(Nx + 1)
lower = np.zeros(Nx)
upper = np.zeros(Nx)
b = np.zeros(Nx + 1)
diagonal[:] = 1
lower[:] = -0.5 * theta * C
upper[:] = 0.5 * theta * C
if FE:
diagonal[:] += 4.0 / 6
lower[:] += 1.0 / 6
upper[:] += 1.0 / 6
upper[0] = 0
lower[-1] = 0
diags = [0, -1, 1]
import scipy.sparse
import scipy.sparse.linalg
A = scipy.sparse.diags(
diagonals=[diagonal, lower, upper],
offsets=[0, -1, 1],
shape=(Nx + 1, Nx + 1),
format="csr",
)
for n in range(0, Nt):
b[1:-1] = u_n[1:-1] + 0.5 * (1 - theta) * C * (u_n[:-2] - u_n[2:])
if FE:
b[1:-1] += 1.0 / 6 * u_n[:-2] + 1.0 / 6 * u_n[:-2] + 4.0 / 6 * u_n[1:-1]
b[0] = u_n[Nx]
b[-1] = u_n[0] # boundary conditions
b[0] = 0
b[-1] = 0 # boundary conditions
u[:] = scipy.sparse.linalg.spsolve(A, b)
if user_action is not None:
user_action(u, x, t, n + 1)
integral[n + 1] = dx * (0.5 * u[0] + 0.5 * u[Nx] + np.sum(u[1:-1]))
u_n, u = u, u_n
t1 = time.perf_counter()
return integral
if __name__ == "__main__":
run(scheme="LW", case="gaussian", C=1, dt=0.01)4.2.2 Test cases
The typical solution \(u\) has the shape of \(I\) and is transported at velocity \(v\) to the right (if \(v>0\)). Let us consider two different initial conditions, one smooth (Gaussian pulse) and one non-smooth (half-truncated cosine pulse):
\[ u(x,0) = Ae^{-\half\left(\frac{x-L/10}{\sigma}\right)^2}, \] \[ u(x,0) = A\cos\left(\frac{5\pi}{L}\left( x - \frac{L}{10}\right)\right),\quad x < \frac{L}{5} \hbox{ else } 0\tp \tag{4.7}\] The parameter \(A\) is the maximum value of the initial condition.
Before doing numerical simulations, we scale the PDE problem and introduce \(\bar x = x/L\) and \(\bar t= vt/L\), which gives \[ \frac{\partial\bar u}{\partial \bar t} + \frac{\partial\bar u}{\partial\bar x} = 0\tp \] The unknown \(u\) is scaled by the maximum value of the initial condition: \(\bar u = u/\max |I(x)|\) such that \(|\bar u(\bar x, 0)|\in [0,1]\). The scaled problem is solved by setting \(v=1\), \(L=1\), and \(A=1\). From now on we drop the bars.
To run our test cases and plot the solution, we make the function
def run_FECS(case):
"""Special function for the FECS case."""
if case == "gaussian":
def I(x):
return np.exp(-0.5 * ((x - L / 10) / sigma) ** 2)
elif case == "cosinehat":
def I(x):
return np.cos(np.pi * 5 / L * (x - L / 10)) if x < L / 5 else 0
L = 1.0
sigma = 0.02
legends = []
def plot(u, x, t, n):
"""Animate and plot every m steps in the same figure."""
plt.figure(1)
if n == 0:
lines = plot(x, u)
else:
lines[0].set_ydata(u)
plt.draw()
plt.figure(2)
m = 40
if n % m != 0:
return
print(
"t=%g, n=%d, u in [%g, %g] w/%d points" % (t[n], n, u.min(), u.max(), x.size)
)
if np.abs(u).max() > 3: # Instability?
return
plt.plot(x, u)
legends.append("t=%g" % t[n])
plt.ion()
U0 = 0
dt = 0.001
C = 1
T = 1
solver(I=I, U0=U0, v=1.0, L=L, dt=dt, C=C, T=T, user_action=plot)
plt.legend(legends, loc="lower left")
plt.savefig("tmp.png")
plt.savefig("tmp.pdf")
plt.axis([0, L, -0.75, 1.1])
plt.show()
def run(scheme="UP", case="gaussian", C=1, dt=0.01):
"""General admin routine for explicit and implicit solvers."""
if case == "gaussian":
def I(x):
return np.exp(-0.5 * ((x - L / 10) / sigma) ** 2)
elif case == "cosinehat":
def I(x):
return np.cos(np.pi * 5 / L * (x - L / 10)) if 0 < x < L / 5 else 0
L = 1.0
sigma = 0.02
global lines # needs to be saved between calls to plot
def plot(u, x, t, n):
"""Plot t=0 and t=0.6 in the same figure."""
plt.figure(1)
global lines
if n == 0:
lines = plt.plot(x, u)
plt.axis([x[0], x[-1], -0.5, 1.5])
plt.xlabel("x")
plt.ylabel("u")
plt.axes().set_aspect(0.15)
plt.savefig("tmp_%04d.png" % n)
plt.savefig("tmp_%04d.pdf" % n)
else:
lines[0].set_ydata(u)
plt.axis([x[0], x[-1], -0.5, 1.5])
plt.title("C=%g, dt=%g, dx=%g" % (C, t[1] - t[0], x[1] - x[0]))
plt.legend(["t=%.3f" % t[n]])
plt.xlabel("x")
plt.ylabel("u")
plt.draw()
plt.savefig("tmp_%04d.png" % n)
plt.figure(2)
eps = 1e-14
if abs(t[n] - 0.6) > eps and abs(t[n] - 0) > eps:
return
print(
"t=%g, n=%d, u in [%g, %g] w/%d points" % (t[n], n, u.min(), u.max(), x.size)
)
if np.abs(u).max() > 3: # Instability?
return
plt.plot(x, u)
plt.draw()
if n > 0:
y = [I(x_ - v * t[n]) for x_ in x]
plt.plot(x, y, "k--")
if abs(t[n] - 0.6) < eps:
filename = ("tmp_%s_dt%s_C%s" % (scheme, t[1] - t[0], C)).replace(".", "")
np.savez(filename, x=x, u=u, u_e=y)
plt.ion()
U0 = 0
T = 0.7
v = 1
codecs = dict(flv="flv", mp4="libx264", webm="libvpx", ogg="libtheora")
import glob
import os
for name in glob.glob("tmp_*.png"):
os.remove(name)
for ext in codecs:
name = "movie.%s" % ext
if os.path.isfile(name):
os.remove(name)
if scheme == "CN":
integral = solver_theta(I, v, L, dt, C, T, user_action=plot, FE=False)
elif scheme == "BE":
integral = solver_theta(I, v, L, dt, C, T, theta=1, user_action=plot)
else:
integral = solver(
I=I, U0=U0, v=v, L=L, dt=dt, C=C, T=T, scheme=scheme, user_action=plot
)
plt.figure(2)
plt.axis([0, L, -0.5, 1.1])
plt.xlabel("$x$")
plt.ylabel("$u$")
plt.axes().set_aspect(0.5) # no effect
plt.savefig("tmp1.png")
plt.savefig("tmp1.pdf")
plt.show()
for codec in codecs:
cmd = "ffmpeg -i tmp_%%04d.png -r 25 -vcodec %s movie.%s" % (codecs[codec], codec)
os.system(cmd)
print("Integral of u:", integral.max(), integral.min())
def solver_theta(I, v, L, dt, C, T, theta=0.5, user_action=None, FE=False):
"""
Full solver for the model problem using the theta-rule
difference approximation in time (no restriction on F,
i.e., the time step when theta >= 0.5).
Vectorized implementation and sparse (tridiagonal)
coefficient matrix.
"""
import time
t0 = time.perf_counter() # for measuring the CPU time
Nt = int(round(T / float(dt)))
t = np.linspace(0, Nt * dt, Nt + 1) # Mesh points in time
dx = v * dt / C
Nx = int(round(L / dx))
x = np.linspace(0, L, Nx + 1) # Mesh points in space
dx = x[1] - x[0]
dt = t[1] - t[0]
C = v * dt / dx
print("dt=%g, dx=%g, Nx=%d, C=%g" % (dt, dx, Nx, C))
u = np.zeros(Nx + 1)
u_n = np.zeros(Nx + 1)
u_nm1 = np.zeros(Nx + 1)
integral = np.zeros(Nt + 1)
for i in range(0, Nx + 1):
u_n[i] = I(x[i])
integral[0] = dx * (0.5 * u_n[0] + 0.5 * u_n[Nx] + np.sum(u_n[1:-1]))
if user_action is not None:
user_action(u_n, x, t, 0)
diagonal = np.zeros(Nx + 1)
lower = np.zeros(Nx)
upper = np.zeros(Nx)
b = np.zeros(Nx + 1)
diagonal[:] = 1
lower[:] = -0.5 * theta * C
upper[:] = 0.5 * theta * C
if FE:
diagonal[:] += 4.0 / 6
lower[:] += 1.0 / 6
upper[:] += 1.0 / 6
upper[0] = 0
lower[-1] = 0
diags = [0, -1, 1]
import scipy.sparse
import scipy.sparse.linalg
A = scipy.sparse.diags(
diagonals=[diagonal, lower, upper],
offsets=[0, -1, 1],
shape=(Nx + 1, Nx + 1),
format="csr",
)
for n in range(0, Nt):
b[1:-1] = u_n[1:-1] + 0.5 * (1 - theta) * C * (u_n[:-2] - u_n[2:])
if FE:
b[1:-1] += 1.0 / 6 * u_n[:-2] + 1.0 / 6 * u_n[:-2] + 4.0 / 6 * u_n[1:-1]
b[0] = u_n[Nx]
b[-1] = u_n[0] # boundary conditions
b[0] = 0
b[-1] = 0 # boundary conditions
u[:] = scipy.sparse.linalg.spsolve(A, b)
if user_action is not None:
user_action(u, x, t, n + 1)
integral[n + 1] = dx * (0.5 * u[0] + 0.5 * u[Nx] + np.sum(u[1:-1]))
u_n, u = u, u_n
t1 = time.perf_counter()
return integral
if __name__ == "__main__":
run(scheme="LW", case="gaussian", C=1, dt=0.01)4.2.3 Bug?
Running either of the test cases, the plot becomes a mess, and the printout of \(u\) values in the plot function reveals that \(u\) grows very quickly. We may reduce \(\Delta t\) and make it very small, yet the solution just grows. Such behavior points to a bug in the code. However, choosing a coarse mesh and performing one time step by hand calculations produces the same numbers as the code, so the implementation seems to be correct. The hypothesis is therefore that the solution is unstable.
4.3 Analysis of the scheme
It is easy to show that a typical Fourier component \[ u(x,t)= B\sin (k(x-ct)) \] is a solution of our PDE for any spatial wave length \(\lambda = 2\pi /k\) and any amplitude \(B\). (Since the PDE to be investigated by this method is homogeneous and linear, \(B\) will always cancel out, so we tend to skip this amplitude, but keep it here in the beginning for completeness.)
A general solution may be viewed as a collection of long and short waves with different amplitudes. Algebraically, the work simplifies if we introduce the complex Fourier component \[ u(x,t)=\Aex e^{ikx}, \] with \[ \Aex=Be^{-ikv\Delta t} = Be^{-iCk\Delta x}\tp \] Note that \(|\Aex| \leq 1\).
It turns out that many schemes also allow a Fourier wave component as solution, and we can use the numerically computed values of \(\Aex\) (denoted \(A\)) to learn about the quality of the scheme. Hence, to analyze the difference scheme we have just implemented, we look at how it treats the Fourier component \[ u_q^n = A^n e^{ikq\Delta x}\tp \] Inserting the numerical component in the scheme, \[ [D_t^+ A e^{ikq\Delta x} + v D_{2x}A e^{ikq\Delta x} = 0]^n_q, \] and making use of (6.5) results in \[ [e^{ikq\Delta x} (\frac{A-1}{\Delta t} + v\frac{1}{\Delta x}i\sin (k\Delta x)) = 0]^n_q, \] which implies \[ A = 1 - iC\sin(k\Delta x)\tp \] The numerical solution features the formula \(A^n\). To find out whether \(A^n\) means growth in time, we rewrite \(A\) in polar form: \(A=A_re^{i\phi}\), for real numbers \(A_r\) and \(\phi\), since we then have \(A^n = A_r^ne^{i\phi n}\). The magnitude of \(A^n\) is \(A_r^n\). In our case, \(A_r = (1 + C^2\sin^2(kx))^{1/2} > 1\), so \(A_r^n\) will increase in time, whereas the exact solution will not. Regardless of \(\Delta t\), we get unstable numerical solutions.
4.4 Leapfrog in time, centered differences in space
4.4.1 Method
Another explicit scheme is to do a “leapfrog” jump over \(2\Delta t\) in time and combine it with central differences in space: \[
[D_{2t} u + vD_{2x} u = 0]_i^n,
\] which results in the updating formula \[
u^{n+1}_i = u^{n-1}**i - C(u**{i+1}^n-u_{i-1}^n)\tp
\] A special scheme is needed to compute \(u^1\), but we leave that problem for now. Anyway, this special scheme can be found in advec1D.py.
4.4.2 Implementation
We now need to work with three time levels and must modify our solver a bit:
Nt = int(round(T/float(dt)))
t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time
...
u = np.zeros(Nx+1)
u_1 = np.zeros(Nx+1)
u_2 = np.zeros(Nx+1)
...
for n in range(0, Nt):
if scheme == 'FE':
for i in range(1, Nx):
u[i] = u_1[i] - 0.5*C*(u_1[i+1] - u_1[i-1])
elif scheme == 'LF':
if n == 0:
for i in range(1, Nx):
...
else:
for i in range(1, Nx+1):
u[i] = u_2[i] - C*(u_1[i] - u_1[i-1])
u_2, u_1, u = u_1, u, u_24.4.3 Running a test case
Let us try a coarse mesh such that the smooth Gaussian initial condition is represented by 1 at mesh node 1 and 0 at all other nodes. This triangular initial condition should then be advected to the right. Choosing scaled variables as \(\Delta t=0.1\), \(T=1\), and \(C=1\) gives the plot in Figure Figure 4.1, which is in fact identical to the exact solution (!).
4.4.4 Running more test cases
We can run two types of initial conditions for \(C=0.8\): one very smooth with a Gaussian function (Figure Figure 4.2) and one with a discontinuity in the first derivative (Figure Figure 4.3). Unless we have a very fine mesh, as in the left plots in the figures, we get small ripples behind the main wave, and this main wave has the amplitude reduced.
4.4.5 Analysis
We can perform a Fourier analysis again. Inserting the numerical Fourier component in the Leapfrog scheme, we get \[ A^2 - i2C\sin(k\Delta x) A - 1 = 0, \] and \[ A = -iC\sin(k\Delta x) \pm \sqrt{1-C^2\sin^2(k\Delta x)}\tp \] Rewriting to polar form, \(A=A_re^{i\phi}\), we see that \(A_r=1\), so the numerical component is neither increasing nor decreasing in time, which is exactly what we want. However, for \(C>1\), the square root can become complex valued, so stability is obtained only as long as \(C\leq 1\).
For all the working schemes to be presented in this chapter, we get the stability condition \(C\leq 1\): \[ \Delta t \leq \frac{\Delta x}{v}\tp \] This is called the CFL condition and applies almost always to successful schemes for advection problems. Of course, one can use Crank-Nicolson or Backward Euler schemes for increased and even unconditional stability (no \(\Delta t\) restrictions), but these have other less desired damping problems.
We introduce \(p=k\Delta x\). The amplification factor now reads \[ A = -iC\sin p \pm \sqrt{1-C^2\sin^2 p}, \] and is to be compared to the exact amplification factor \[ \Aex = e^{-ikv\Delta t} = e^{-ikC\Delta x} = e^{-iCp}\tp \] Section Section 4.10 compares numerical amplification factors of many schemes with the exact expression.
4.5 Upwind differences in space
Since the PDE reflects transport of information along with a flow in positive \(x\) direction, when \(v>0\), it could be natural to go (what is called) upstream and not downstream in the spatial derivative to collect information about the change of the function. That is, we approximate \[ \frac{\partial u}{\partial x}(x_i,t_n)\approx [D^-_x u]^n_i = \frac{u^n_{i} - u^n_{i-1}}{\Delta x}\tp \] This is called an upwind difference (the corresponding difference in the time direction would be called a backward difference, and we could use that name in space too, but upwind is the common name for a difference against the flow in advection problems). This spatial approximation does magic compared to the scheme we had with Forward Euler in time and centered difference in space. With an upwind difference, \[ [D^+_t u + vD^-_x u = 0]^n_i, \tag{4.8}\] written out as \[ u^{n+1}_i = u^n_i - C(u^{n}**{i}-u^{n}**{i-1}), \] gives a generally popular and robust scheme that is stable if \(C\leq 1\). As with the Leapfrog scheme, it becomes exact if \(C=1\), exactly as shown in Figure Figure 4.1. This is easy to see since \(C=1\) gives the property (4.6). However, any \(C<1\) gives a significant reduction in the amplitude of the solution, which is a purely numerical effect, see Figures Figure 4.4 and Figure 4.5. Experiments show, however, that reducing \(\Delta t\) or \(\Delta x\), while keeping \(C\) reduces the error.
The amplification factor can be computed using the formula (6.4), \[ \frac{A - 1}{\Delta t} + \frac{v}{\Delta x}(1 - e^{-ik\Delta x}) = 0, \] which means \[ A = 1 - C(1 - \cos(p) - i\sin(p))\tp \] For \(C<1\) there is, unfortunately, non-physical damping of discrete Fourier components, giving rise to reduced amplitude of \(u^n_i\) as in Figures Figure 4.4 and Figure 4.5. The damping seen in these figures is quite severe. Stability requires \(C\leq 1\).
One can interpret the upwind difference as extra, artificial diffusion in the equation. Solving \[ \frac{\partial u}{\partial t} + v\frac{\partial u}{\partial x} = \nu\frac{\partial^2 u}{\partial x^2}, \] by a forward difference in time and centered differences in space, \[ D^+**t u + vD**{2x} u = \nu D_xD_x u]^n_i, \] actually gives the upwind scheme (4.8) if \(\nu = v\Delta x/2\). That is, solving the PDE \(u_t + vu_x=0\) by centered differences in space and forward difference in time is unsuccessful, but by adding some artificial diffusion \(\nu u_{xx}\), the method becomes stable: \[ \frac{\partial u}{\partial t} + v \frac{\partial u}{\partial x} = \left(\dfc + \frac{v\Delta x}{2}\right) \frac{\partial^2 u}{\partial x^2}\tp \]
4.6 Periodic boundary conditions
So far, we have given the value on the left boundary, \(u_0^n\), and used the scheme to propagate the solution signal through the domain. Often, we want to follow such signals for long time series, and periodic boundary conditions are then relevant since they enable a signal that leaves the right boundary to immediately enter the left boundary and propagate through the domain again.
The periodic boundary condition is \[
u(0,t) = u(L,t),\quad u_0^n = u_{N_x}^n\tp
\] It means that we in the first equation, involving \(u_0^n\), insert \(u_{N_x}^n\), and that we in the last equation, involving \(u^{n+1}_{N_x}\) insert \(u^{n+1}_0\). Normally, we can do this in the simple way that u_1[0] is updated as u_1[Nx] at the beginning of a new time level.
In some schemes we may need \(u^{n}_{N_x+1}\) and \(u^{n}_{-1}\). Periodicity then means that these values are equal to \(u^n_1\) and \(u^n_{N_x-1}\), respectively. For the upwind scheme, it is sufficient to set u_1[0]=u_1[Nx] at a new time level before computing u[1]. This ensures that u[1] becomes right and at the next time level u[0] at the current time level is correctly updated. For the Leapfrog scheme we must update u[0] and u[Nx] using the scheme:
if periodic_bc:
i = 0
u[i] = u_2[i] - C*(u_1[i+1] - u_1[Nx-1])
for i in range(1, Nx):
u[i] = u_2[i] - C*(u_1[i+1] - u_1[i-1])
if periodic_bc:
u[Nx] = u[0]4.7 Implementation
4.7.1 Test condition
Analytically, we can show that the integral in space under the \(u(x,t)\) curve is constant:
\[\begin{align*} \int_0^L \left(\frac{\partial u}{\partial t} + v\frac{\partial u}{\partial x} \right) dx &= 0\\ \frac{\partial }{\partial t} \int_0^L udx &= - \int_0^L v\frac{\partial u}{\partial x}dx\\ \frac{\partial u}{\partial t} \int_0^L udx &= [v u]_0^L =0 \end{align*}\] as long as \(u(0)=u(L)=0\). We can therefore use the property \[ \int_0^L u(x,t)dx = \hbox{const} \] as a partial verification during the simulation. Now, any numerical method with \(C\neq 1\) will deviate from the constant, expected value, so the integral is a measure of the error in the scheme. The integral can be computed by the Trapezoidal integration rule
dx*(0.5*u[0] + 0.5*u[Nx] + np.sum(u[1:-1]))if u is an array holding the solution.
4.7.2 The code
An appropriate solver function for multiple schemes may go as shown below.
def solver(I, U0, v, L, dt, C, T, user_action=None,
scheme='FE', periodic_bc=True):
Nt = int(round(T/float(dt)))
t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time
dx = v*dt/C
Nx = int(round(L/dx))
x = np.linspace(0, L, Nx+1) # Mesh points in space
dx = x[1] - x[0]
dt = t[1] - t[0]
C = v*dt/dx
print 'dt=%g, dx=%g, Nx=%d, C=%g' % (dt, dx, Nx, C)
u = np.zeros(Nx+1)
u_n = np.zeros(Nx+1)
u_nm1 = np.zeros(Nx+1)
integral = np.zeros(Nt+1)
for i in range(0, Nx+1):
u_n[i] = I(x[i])
u[0] = U0
integral[0] = dx*(0.5*u_n[0] + 0.5*u_n[Nx] + np.sum(u_n[1:-1]))
if user_action is not None:
user_action(u_n, x, t, 0)
for n in range(0, Nt):
if scheme == 'FE':
if periodic_bc:
i = 0
u[i] = u_n[i] - 0.5*C*(u_n[i+1] - u_n[Nx])
u[Nx] = u[0]
for i in range(1, Nx):
u[i] = u_n[i] - 0.5*C*(u_n[i+1] - u_n[i-1])
elif scheme == 'LF':
if n == 0:
if periodic_bc:
i = 0
u_n[i] = u_n[Nx]
for i in range(1, Nx+1):
u[i] = u_n[i] - C*(u_n[i] - u_n[i-1])
else:
if periodic_bc:
i = 0
u[i] = u_nm1[i] - C*(u_n[i+1] - u_n[Nx-1])
for i in range(1, Nx):
u[i] = u_nm1[i] - C*(u_n[i+1] - u_n[i-1])
if periodic_bc:
u[Nx] = u[0]
elif scheme == 'UP':
if periodic_bc:
u_n[0] = u_n[Nx]
for i in range(1, Nx+1):
u[i] = u_n[i] - C*(u_n[i] - u_n[i-1])
else:
raise ValueError('scheme="%s" not implemented' % scheme)
if not periodic_bc:
u[0] = U0
integral[n+1] = dx*(0.5*u[0] + 0.5*u[Nx] + np.sum(u[1:-1]))
if user_action is not None:
user_action(u, x, t, n+1)
u_nm1, u_n, u = u_n, u, u_nm1
return integral4.7.3 Solving a specific problem
We need to call up the solver function in some kind of administering problem solving function that can solve specific problems and make appropriate visualization. The function below makes both static plots, screen animation, and hard copy videos in various formats.
def run(scheme='UP', case='gaussian', C=1, dt=0.01):
"""General admin routine for explicit and implicit solvers."""
if case == 'gaussian':
def I(x):
return np.exp(-0.5*((x-L/10)/sigma)**2)
elif case == 'cosinehat':
def I(x):
return np.cos(np.pi*5/L*(x - L/10)) if x < L/5 else 0
L = 1.0
sigma = 0.02
global lines # needs to be saved between calls to plot
def plot(u, x, t, n):
"""Plot t=0 and t=0.6 in the same figure."""
plt.figure(1)
global lines
if n == 0:
lines = plt.plot(x, u)
plt.axis([x[0], x[-1], -0.5, 1.5])
plt.xlabel('x'); plt.ylabel('u')
plt.axes().set_aspect(0.15)
plt.savefig('tmp_%04d.png' % n)
plt.savefig('tmp_%04d.pdf' % n)
else:
lines[0].set_ydata(u)
plt.axis([x[0], x[-1], -0.5, 1.5])
plt.title('C=%g, dt=%g, dx=%g' %
(C, t[1]-t[0], x[1]-x[0]))
plt.legend(['t=%.3f' % t[n]])
plt.xlabel('x'); plt.ylabel('u')
plt.draw()
plt.savefig('tmp_%04d.png' % n)
plt.figure(2)
eps = 1E-14
if abs(t[n] - 0.6) > eps and abs(t[n] - 0) > eps:
return
print 't=%g, n=%d, u in [%g, %g] w/%d points' % \
(t[n], n, u.min(), u.max(), x.size)
if np.abs(u).max() > 3: # Instability?
return
plt.plot(x, u)
plt.hold('on')
plt.draw()
if n > 0:
y = [I(x_-v*t[n]) for x_ in x]
plt.plot(x, y, 'k--')
if abs(t[n] - 0.6) < eps:
filename = ('tmp_%s_dt%s_C%s' % \
(scheme, t[1]-t[0], C)).replace('.', '')
np.savez(filename, x=x, u=u, u_e=y)
plt.ion()
U0 = 0
T = 0.7
v = 1
codecs = dict(flv='flv', mp4='libx264', webm='libvpx',
ogg='libtheora')
import glob, os
for name in glob.glob('tmp_*.png'):
os.remove(name)
for ext in codecs:
name = 'movie.%s' % ext
if os.path.isfile(name):
os.remove(name)
integral = solver(
I=I, U0=U0, v=v, L=L, dt=dt, C=C, T=T,
scheme=scheme, user_action=plot)
plt.figure(2)
plt.axis([0, L, -0.5, 1.1])
plt.xlabel('$x$'); plt.ylabel('$u$')
plt.savefig('tmp1.png'); plt.savefig('tmp1.pdf')
plt.show()
for codec in codecs:
cmd = 'ffmpeg -i tmp_%%04d.png -r 25 -vcodec %s movie.%s' % \
(codecs[codec], codec)
os.system(cmd)
print 'Integral of u:', integral.max(), integral.min()The complete code is found in the file advec1D.py.
4.8 A Crank-Nicolson discretization in time and centered differences in space
Another obvious candidate for time discretization is the Crank-Nicolson method combined with centered differences in space: \[ [D_t u]^n_i + v\half([D_{2x} u]^{n+1}**i + [D**{2x} u]^{n}_i) = 0\tp \] It can be nice to include the Backward Euler scheme too, via the \(\theta\)-rule, \[ [D_t u]^n_i + v\theta [D_{2x} u]^{n+1}**i + v(1-\theta)[D**{2x} u]^{n}_i = 0\tp \] When \(\theta\) is different from zero, this gives rise to an implicit scheme, \[ u^{n+1}**i + \frac{\theta}{2} C (u^{n+1}**{i+1} - u^{n+1}_{i-1}) = u^n_i - \frac{1-\theta}{2} C (u^{n}**{i+1} - u^{n}**{i-1}) \] for \(i=1,\ldots,N_x-1\). At the boundaries we set \(u=0\) and simulate just to the point of time when the signal hits the boundary (and gets reflected). \[ u^{n+1}**0 = u^{n+1}**{N_x} = 0\tp \] The elements on the diagonal in the matrix become: \[ A_{i,i} = 1,\quad i=0,\ldots,N_x\tp \] On the subdiagonal and superdiagonal we have \[ A_{i-1,i} = -\frac{\theta}{2} C,\quad A_{i+1,i} = \frac{\theta}{2} C,\quad i=1,\ldots,N_x-1, \] with \(A_{0,1}=0\) and \(A_{N_x-1,N_x}=0\) due to the known boundary conditions. And finally, the right-hand side becomes
\[\begin{align*} b_0 &= u^n_{N_x}\\ b_i &= u^n_i - \frac{1-\theta}{2} C (u^{n}**{i+1} - u^{n}**{i-1}),\quad i=1,\ldots,N_x-1\\ b_{N_x} &= u^n_0 \end{align*}\]
The dispersion relation follows from inserting \(u^n_q = A^ne^{ikx}\) and using the formula (6.5) for the spatial differences: \[ A = \frac{1 - (1-\theta) i C\sin p}{1 + \theta i C\sin p}\tp \]
Figure Figure 4.6 depicts a numerical solution for \(C=0.8\) and the Crank-Nicolson with severe oscillations behind the main wave. These oscillations are damped as the mesh is refined. Switching to the Backward Euler scheme removes the oscillations, but the amplitude is significantly reduced. One could expect that the discontinuous derivative in the initial condition of the half a cosine wave would make even stronger demands on producing a smooth profile, but Figure Figure 4.7 shows that also here, Backward-Euler is capable of producing a smooth profile. All in all, there are no major differences between the Gaussian initial condition and the half a cosine condition for any of the schemes.
4.9 The Lax-Wendroff method
The Lax-Wendroff method is based on three ideas:
- Express the new unknown \(u^{n+1}_i\) in terms of known quantities at \(t=t_n\) by means of a Taylor polynomial of second degree.
- Replace time-derivatives at \(t=t_n\) by spatial derivatives, using the PDE.
- Discretize the spatial derivatives by second-order differences so we achieve a scheme of accuracy \(\Oof{\Delta t^2} + \Oof{\Delta x^2}\).
Let us follow the recipe. First we have the three-term Taylor polynomial, \[ u^{n+1}_i = u^n_i + \Delta t\left(\frac{\partial u}{\partial t}\right)^n_i + \frac{1}{2}\Delta t^2\left(\frac{\partial^2 u}{\partial t^2}\right)^n_i\tp \] From the PDE we have that temporal derivatives can be substituted by spatial derivatives: \[ \frac{\partial u}{\partial t} = -v\frac{\partial u}{\partial x}, \] and furthermore, \[ \frac{\partial ^2 u}{\partial t^2} = v^2\frac{\partial^2 u}{\partial x^2}\tp \] Inserted in the Taylor polynomial formula, we get \[ u^{n+1}_i = u^n_i -v \Delta t\left(\frac{\partial u}{\partial x}\right)^n_i + \frac{1}{2}\Delta t^2 v^2 \left(\frac{\partial^2 u}{\partial x^2}\right)^n_i\tp \] To obtain second-order accuracy in space we now use central differences: \[ u^{n+1}_i = u^n_i -v \Delta t [D_{2x} u]^n_i + \frac{1}{2}\Delta t^2 v^2 [D_xD_x u]^n_i, \] or written out, \[ u^{n+1}_i = u^n_i - \frac{1}{2} C (u^{n}**{i+1} - u^{n}**{i-1}) + \frac{1}{2} C^2 (u^{n}_{i+1}-2u^n_i+u^n_{i-1})\tp \] This is the explicit Lax-Wendroff scheme.
From the formulas above, we notice that the Lax-Wendroff method is nothing but a Forward Euler, central difference in space scheme, which we have shown to be useless because of chronic instability, plus an artificial diffusion term of strength \(\half\Delta t v^2\). It means that we can take an unstable scheme and add some diffusion to stabilize it. This is a common trick to deal with advection problems. Sometimes, the real physical diffusion is not sufficiently large to make schemes stable, so then we also add artificial diffusion.
From an analysis similar to the ones carried out above, we get an amplification factor for the Lax-Wendroff method that equals \[ A = 1 - iC\sin p - 2C^2\sin^2 (p/2)\tp \] This means that \(|A|=1\) and also that we have an exact solution if \(C=1\)!
4.10 Analysis of dispersion relations
We have developed expressions for \(A(C,p)\) in the exact solution \(u_q^n=A^ne^{ikq\Delta x}\) of the discrete equations. Note that the Fourier component that solves the original PDE problem has no damping and moves with constant velocity \(v\). There are two basic errors in the numerical Fourier component: there may be damping and the wave velocity may depend on \(C\) and \(p=k\Delta x\).
The shortest wavelength that can be represented is \(\lambda = 2\Delta x\). The corresponding \(k\) is \(k=2\pi/\lambda = \pi/\Delta x\), so \(p=k\Delta x\in (0,\pi]\).
Given a complex \(A\) as a function of \(C\) and \(p\), how can we visualize it? The two key ingredients in \(A\) is the magnitude, reflecting damping or growth of the wave, and the angle, closely related to the velocity of the wave. The Fourier component \[ D^n e^{ik(x-ct)} \] has damping \(D\) and wave velocity \(c\). Let us express our \(A\) in polar form, \(A = A_re^{-i\phi}\), and insert this expression in our discrete component \(u_q^n = A^ne^{ikq\Delta x} = A^ne^{ikx}\): \[ u^n_q = A_r^n e^{-i\phi n} e^{ikx} = A_r^n e^{i(kx - n\phi)} = A_r^ne^{i(k(x - ct))}, \] for \[ c = \frac{\phi}{k\Delta t}\tp \] Now, \[ k\Delta t = \frac{Ck\Delta x}{v}=\frac{Cp}{v}, \] so \[ c = \frac{\phi v}{Cp}\tp \] An appropriate dimensionless quantity to plot is the scaled wave velocity \(c/v\): \[ \frac{c}{v} = \frac{\phi}{Cp}\tp \] Figures Figure 4.8–Figure 4.13 contain dispersion curves, velocity and damping, for various values of \(C\). The horizontal axis shows the dimensionless frequency \(p\) of the wave, while the figures to the left illustrate the error in wave velocity \(c/v\) (should ideally be 1 for all \(p\)), and the figures to the right display the absolute value (magnitude) of the damping factor \(A_r\). The curves are labeled according to the table below.
| Label | Method |
|---|---|
| FE | Forward Euler in time, centered difference in space |
| LF | Leapfrog in time, centered difference in space |
| UP | Forward Euler in time, upwind difference in space |
| CN | Crank-Nicolson in time, centered difference in space |
| LW | Lax-Wendroff’s method |
| BE | Backward Euler in time, centered difference in space |
The total damping after some time \(T=n\Delta t\) is reflected by \(A_r(C,p)^n\). Since normally \(A_r<1\), the damping goes like \(A_r^{1/\Delta t}\) and approaches zero as \(\Delta t\rightarrow 0\). The only way to reduce damping is to increase \(C\) and/or the mesh resolution.
We can learn a lot from the dispersion relation plots. For example, looking at the plots for \(C=1\), the schemes LW, UP, and LF has no amplitude reduction, but LF has wrong phase velocity for the shortest wave in the mesh. This wave does not (normally) have enough amplitude to be seen, so for all practical purposes, there is no damping or wrong velocity of the individual waves, so the total shape of the wave is also correct. For the CN scheme, see Figure Figure 4.6, each individual wave has its amplitude, but they move with different velocities, so after a while, we see some of these waves lagging behind. For the BE scheme, see Figure Figure 4.7, all the shorter waves are so heavily dampened that we cannot see them after a while. We see only the longest waves, which have slightly wrong velocity, but visible amplitudes are sufficiently equal to produce what looks like a smooth profile.
Another feature was that the Leapfrog method produced oscillations, while the upwind scheme did not. Since the Leapfrog method does not dampen the shorter waves, which have wrong wave velocities of order 10 percent, we can see these waves as noise. The upwind scheme, however, dampens these waves. The same effect is also present in the Lax-Wendroff scheme, but the damping of the intermediate waves is hardly present, so there is visible noise in the total signal.
We realize that, compared to pure truncation error analysis, dispersion analysis sheds more light on the behavior of the computational schemes. Truncation analysis just says that Lax-Wendroff is better than upwind, because of the increased order in time, but most people would say upwind is the better one when looking at the plots.
4.11 Stationary 1D advection-diffusion
Now we pay attention to a physical process where advection (or convection) is in balance with diffusion: \[ v\frac{du}{dx} = \dfc\frac{d^2 u}{dx^2}\tp \tag{4.9}\] For simplicity, we assume \(v\) and \(\dfc\) to be constant, but the extension to the variable-coefficient case is trivial. This equation can be viewed as the stationary limit of the corresponding time-dependent problem \[ \frac{\partial u}{\partial t} + v\frac{\partial u}{\partial x} = \dfc\frac{\partial^2 u}{\partial x^2}\tp \tag{4.10}\]
Equations of the form (4.9) or (4.10) arise from transport phenomena, either mass or heat transport. One can also view the equations as a simple model problem for the Navier-Stokes equations. With the chosen boundary conditions, the differential equation problem models the phenomenon of a boundary layer, where the solution changes rapidly very close to the boundary. This is a characteristic of many fluid flow problems, which makes strong demands to numerical methods. The fundamental numerical difficulty is related to non-physical oscillations of the solution (instability) if the first-derivative spatial term dominates over the second-derivative term.
4.12 A simple model problem
We consider (4.9) on \([0,L]\) equipped with the boundary conditions \(u(0)=U_0\), \(u(L)=U_L\). By scaling we can reduce the number of parameters in the problem. We scale \(x\) by \(\bar x = x/L\), and \(u\) by \[ \bar u = \frac{u - U_0}{U_L-U_0}\tp \] Inserted in the governing equation we get \[ \frac{v(U_L-U_0)}{L}\frac{d\bar u}{d\bar x} = \frac{\dfc(U_L-U_0)}{L^2}\frac{d^2\bar u}{d\bar x^2},\quad \bar u(0)=0,\ \bar u(1)=1\tp \] Dropping the bars is common. We can then simplify to \[ \frac{du}{dx} = \epsilon\frac{d^2 u}{d x^2},\quad u(0)=0,\ u(1)=1\tp \tag{4.11}\]
There are two competing effects in this equation: the advection term transports signals to the right, while the diffusion term transports signals to the left and the right. The value \(u(0)=0\) is transported through the domain if \(\epsilon\) is small, and \(u\approx 0\) except in the vicinity of \(x=1\), where \(u(1)=1\) and the diffusion transports some information about \(u(1)=1\) to the left. For large \(\epsilon\), diffusion dominates and the \(u\) takes on the “average” value, i.e., \(u\) gets a linear variation from 0 to 1 throughout the domain.
It turns out that we can find an exact solution to the differential equation problem and also to many of its discretizations. This is one reason why this model problem has been so successful in designing and investigating numerical methods for mixed convection/advection and diffusion. The exact solution reads \[ \uex (x) = \frac{e^{x/\epsilon} - 1}{e^{1/\epsilon} - 1}\tp \] The forthcoming plots illustrate this function for various values of \(\epsilon\).
4.13 A centered finite difference scheme
The most obvious idea to solve (4.11) is to apply centered differences: \[ [D_{2x} u = \epsilon D_xD_x u]_i \] for \(i=1,\ldots,N_x-1\), with \(u_0=0\) and \(u_{N_x}=1\). Note that this is a coupled system of algebraic equations involving \(u_0,\ldots,u_{N_x}\).
Written out, the scheme becomes a tridiagonal system \[ A_{i-1,i}u_{i-1} + A_{i,i}u_i + A_{i+1.i}u_{i+1} = 0, \] for \(i=1,\ldots,N_x-1\)
\[\begin{align*} A_{0,0} &= 1,\\ A_{i-1,i} &= -\frac{1}{\Delta x} -\epsilon\frac{1}{\Delta x^2},\\ A_{i,i} &= 2\epsilon\frac{1}{\Delta x^2},\\ A_{i,i+1} &= \frac{1}{\Delta x} -\epsilon\frac{1}{\Delta x^2},\\ A_{N_x,N_x} &= 1\tp \end{align*}\] The right-hand side of the linear system is zero except \(b_{N_x}=1\).
Figure Figure 4.14 shows reasonably accurate results with $N_x=20 $ and \(N_x=40\) cells in \(x\) direction and a value of \(\epsilon = 0.1\). Decreasing \(\epsilon\) to \(0.01\) leads to oscillatory solutions as depicted in Figure Figure 4.15. This is, unfortunately, a typical phenomenon in this type of problem: non-physical oscillations arise for small \(\epsilon\) unless the resolution \(N_x\) is big enough. Exercise Section 4.18 develops a precise criterion: \(u\) is oscillation-free if \[ \Delta x \leq \frac{2}{\epsilon}\tp \] If we take the present model as a simplified model for a viscous boundary layer in real, industrial fluid flow applications, \(\epsilon\sim 10^{-6}\) and millions of cells are required to resolve the boundary layer. Fortunately, this is not strictly necessary as we have methods in the next section to overcome the problem!
A suitable solver for doing the experiments is presented below.
import numpy as np
def solver(eps, Nx, method="centered"):
"""
Solver for the two point boundary value problem u'=eps*u'',
u(0)=0, u(1)=1.
"""
x = np.linspace(0, 1, Nx + 1) # Mesh points in space
dx = x[1] - x[0]
u = np.zeros(Nx + 1)
diagonal = np.zeros(Nx + 1)
lower = np.zeros(Nx)
upper = np.zeros(Nx)
b = np.zeros(Nx + 1)
if method == "centered":
diagonal[:] = 2 * eps / dx**2
lower[:] = -1 / dx - eps / dx**2
upper[:] = 1 / dx - eps / dx**2
elif method == "upwind":
diagonal[:] = 1 / dx + 2 * eps / dx**2
lower[:] = 1 / dx - eps / dx**2
upper[:] = -eps / dx**2
upper[0] = 0
lower[-1] = 0
diagonal[0] = diagonal[-1] = 1
b[-1] = 1.0
diags = [0, -1, 1]
import scipy.sparse
import scipy.sparse.linalg
A = scipy.sparse.diags(
diagonals=[diagonal, lower, upper],
offsets=[0, -1, 1],
shape=(Nx + 1, Nx + 1),
format="csr",
)
u[:] = scipy.sparse.linalg.spsolve(A, b)
return u, x4.14 Remedy: upwind finite difference scheme
The scheme can be stabilized by letting the advective transport term, which is the dominating term, collect its information in the flow direction, i.e., upstream or upwind of the point in question. So, instead of using a centered difference \[ \frac{du}{dx}**i\approx \frac{u**{i+1}-u_{i-1}}{2\Delta x}, \] we use the one-sided upwind difference \[ \frac{du}{dx}**i\approx \frac{u**{i}-u_{i-1}}{\Delta x}, \] in case \(v>0\). For \(v<0\) we set \[ \frac{du}{dx}**i\approx \frac{u**{i+1}-u_{i}}{\Delta x}, \] On compact operator notation form, our upwind scheme can be expressed as \[ [D^-_x u = \epsilon D_xD_x u]_i \] provided \(v>0\) (and \(\epsilon > 0\)).
We write out the equations and implement them as shown in the program in Section Section 4.13. The results appear in Figures Figure 4.16 and Figure 4.17: no more oscillations!
We see that the upwind scheme is always stable, but it gives a thicker boundary layer when the centered scheme is also stable. Why the upwind scheme is always stable is easy to understand as soon as we undertake the mathematical analysis in Exercise Section 4.18. Moreover, the thicker layer (seemingly larger diffusion) can be understood by doing Exercise Section 4.19.
It turns out that one can introduce a linear combination of the centered and upwind differences for the first-derivative term in this model problem. One can then adjust the weight in the linear combination so that the numerical solution becomes identical to the analytical solution of the differential equation problem at any mesh point.
Now it is time to combine time-dependency, convection (advection) and diffusion into one equation: \[ \frac{\partial u}{\partial t} + v\frac{\partial u}{\partial x} = \dfc\frac{\partial^2 u}{\partial x^2}\tp \tag{4.12}\]
4.14.1 Analytical insight
The diffusion is now dominated by convection, a wave, and diffusion, a loss of amplitude. One possible analytical solution is a traveling Gaussian function \[ u(x,t) = B\exp{\left(-\left(\frac{x - vt}{4at}\right)\right)}\tp \] This function moves with velocity \(v>0\) to the right (\(v<0\) to the left) due to convection, but at the same time we have a damping \(e^{-16a^2t^2}\) from diffusion.
4.15 Forward in time, centered in space scheme
The Forward Euler for the diffusion equation is a successful scheme, but it has a very strict stability condition. The similar Forward in time, centered in space strategy always gives unstable solutions for the advection PDE. What happens when we have both diffusion and advection present at once? \[ [D_t u + vD_{2x} u = \dfc D_xD_x u + f]_i^n\tp \] We expect that diffusion will stabilize the scheme, but that advection will destabilize it.
Another problem is non-physical oscillations, but not growing amplitudes, due to centered differences in the advection term. There will hence be two types of instabilities to consider. Our analysis showed that pure advection with centered differences in space needs some artificial diffusion to become stable (and then it produces upwind differences for the advection term). Adding more physical diffusion should further help the numerics to stabilize the non-physical oscillations.
The scheme is quickly implemented, but suffers from the need for small space and time steps, according to this reasoning. A better approach is to get rid of the non-physical oscillations in space by simply applying an upwind difference on the advection term.
4.16 Forward in time, upwind in space scheme
A good approximation for the pure advection equation is to use upwind discretization of the advection term. We also know that centered differences are good for the diffusion term, so let us combine these two discretizations: \[ [D_t u + vD^-_{x} u = \dfc D_xD_x u + f]_i^n, \] for \(v>0\). Use \(vD^+ u\) if \(v<0\). In this case the physical diffusion and the extra numerical diffusion \(v\Delta x/2\) will stabilize the solution, but give an overall too large reduction in amplitude compared with the exact solution.
We may also interpret the upwind difference as artificial numerical diffusion and centered differences in space everywhere, so the scheme can be expressed as \[ [D_t u + vD^-_{2x} u = \dfc \frac{v\Delta x}{2} D_xD_x u + f]_i^n\tp \]
4.17 Applications of advection equations
There are two major areas where advection and convection applications arise: transport of a substance and heat transport in a fluid. To derive the models, we may look at the similar derivations of diffusion models in Section Section 3.66, but change the assumption from a solid to fluid medium. This gives rise to the extra advection or convection term \(\vv\cdot\nabla u\). We briefly show how this is done.
Normally, transport in a fluid is dominated by the fluid flow and not diffusion, so we can neglect diffusion compared to advection or convection. The end result is anyway an equation of the form \[ \frac{\partial u}{\partial t} + \vv\cdot\nabla u = 0\tp \] ## Transport of a substance {#sec-advec-app-mass}
The diffusion of a substance in Section Section 3.66.1 takes place in a solid medium, but in a fluid we can have two transport mechanisms: one by diffusion and one by advection. The latter arises from the fact that the substance particles are moved with the fluid velocity \(\vv\) such that the effective flux now consists of two and not only one component as in (3.59): \[ \q = -\dfc\nabla c + \vv\c\tp \] Inserted in the equation \(\nabla\cdot\q = 0\) we get the extra advection term \(\nabla\cdot (\vv\c)\). Very often we deal with incompressible flows, \(\nabla\cdot\vv = 0\) such that the advective term becomes \(\vv\cdot\nabla c\). The mass transport equation for a substance then reads \[ \frac{\partial c}{\partial t} + \vv\cdot\nabla c = \dfc\nabla^2 c\tp \] ## Transport of heat in fluids {#sec-advec-app-heat}
The derivation of the heat equation in Section Section 3.66.2 is limited to heat transport in solid bodies. If we turn the attention to heat transport in fluids, we get a material derivative of the internal energy in (3.61), \[ \frac{De}{dt} = - \nabla\cdot\q, \] and more terms if work by stresses is also included, where \[ \frac{De}{dt} = \frac{\partial e}{\partial t} + \vv\cdot\nabla e, \] \(\vv\) being the velocity of the fluid. The convective term \(\vv\cdot\nabla e\) must therefore be added to the governing equation, resulting typically in \[ \varrho c\left(\frac{\partial T}{\partial t} + \vv\cdot\nabla T\right) = \nabla\cdot(k\nabla T) + f, \tag{4.13}\] where \(f\) is some external heating inside the medium.
4.18 Exercise: Analyze 1D stationary convection-diffusion problem
Explain the observations in the numerical experiments from Sections Section 4.13 and Section 4.14 by finding exact numerical solutions.
\(A\) is an unknown constant and \(i\) is a mesh point counter. There are two solutions for \(A\), so the general solution is a linear combination of the two, where the constants in the linear combination are determined from the boundary conditions.
4.19 Exercise: Interpret upwind difference as artificial diffusion
Consider an upwind, one-sided difference approximation to a term \(du/dx\) in a differential equation. Show that this formula can be expressed as a centered difference plus an artificial diffusion term of strength proportional to \(\Delta x\). This means that introducing an upwind difference also means introducing extra diffusion of order \(\Oof{\Delta x}\).
4.20 Advection Schemes with Devito
Having understood the mathematical properties and challenges of advection schemes in the previous sections, we now implement these methods using Devito’s symbolic framework. Devito allows us to write the discrete equations in a form close to the mathematical notation while generating optimized code automatically.
4.20.1 The Advection Equation
The 1D linear advection equation is:
\[ \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0 \tag{4.14}\]
where \(c\) is the advection velocity (assumed constant and positive). The exact solution is:
\[ u(x, t) = I(x - ct) \]
which represents the initial condition \(I(x)\) traveling to the right at velocity \(c\) without change in shape.
4.20.2 Devito Implementation Patterns
Unlike diffusion and wave equations, the advection equation requires careful treatment of the spatial derivative. Centered differences lead to instability (as we saw with the FTCS scheme), so we need alternative approaches:
| Scheme | Spatial Discretization | Order | Key Property |
|---|---|---|---|
| Upwind | Backward difference | 1st | Stable, diffusive |
| Lax-Wendroff | Centered + diffusion | 2nd | Less diffusion, some dispersion |
| Lax-Friedrichs | Averaged neighbors | 1st | Very diffusive but robust |
All schemes require the CFL condition: \(C = c\Delta t/\Delta x \leq 1\).
4.20.3 Comparison with Wave and Diffusion Equations
The advection equation differs fundamentally from the diffusion and wave equations we’ve solved previously:
| Property | Diffusion | Wave | Advection |
|---|---|---|---|
time_order |
1 | 2 | 1 |
| Spatial deriv. | 2nd (.dx2) |
2nd (.laplace) |
1st (.dx) |
| Stability | \(F \leq 0.5\) | \(C \leq 1\) | \(C \leq 1\) |
| Centered space | Stable | Stable | Unstable |
| Information | Spreads both ways | Spreads both ways | One direction |
The key difference is that advection has directional information flow, which requires using upwind differences rather than centered differences.
4.20.4 Upwind Scheme Implementation
The upwind scheme uses a backward difference for the spatial derivative when \(c > 0\):
\[ \frac{u^{n+1}_i - u^n_i}{\Delta t} + c\frac{u^n_i - u^n_{i-1}}{\Delta x} = 0 \]
which gives the update formula:
\[ u^{n+1}_i = u^n_i - C(u^n_i - u^n_{i-1}) \tag{4.15}\]
In Devito, we express this using shifted indexing:
from devito import Grid, TimeFunction, Eq, Operator, Constant
import numpy as np
def solve_advection_upwind(L, c, Nx, T, C, I):
"""Upwind scheme for 1D advection."""
# Grid setup
dx = L / Nx
dt = C * dx / c
grid = Grid(shape=(Nx + 1,), extent=(L,))
x_dim, = grid.dimensions
u = TimeFunction(name='u', grid=grid, time_order=1, space_order=1)
# Set initial condition
x_coords = np.linspace(0, L, Nx + 1)
u.data[0, :] = I(x_coords)
# Courant number as constant
courant = Constant(name='C', value=C)
# Upwind stencil: u^{n+1} = u - C*(u - u[x-dx])
u_minus = u.subs(x_dim, x_dim - x_dim.spacing)
stencil = u - courant * (u - u_minus)
update = Eq(u.forward, stencil)
op = Operator([update])
# ... time stepping loopThe key line is:
u_minus = u.subs(x_dim, x_dim - x_dim.spacing)This creates a reference to \(u^n_{i-1}\) by substituting x_dim - x_dim.spacing for x_dim in the TimeFunction u.
4.20.5 Lax-Wendroff Scheme Implementation
The Lax-Wendroff scheme achieves second-order accuracy by including both a centered advection term and a diffusion-like correction:
\[ u^{n+1}_i = u^n_i - \frac{C}{2}(u^n_{i+1} - u^n_{i-1}) + \frac{C^2}{2}(u^n_{i+1} - 2u^n_i + u^n_{i-1}) \]
This can be written using Devito’s derivative operators:
def solve_advection_lax_wendroff(L, c, Nx, T, C, I):
"""Lax-Wendroff scheme for 1D advection."""
dx = L / Nx
dt = C * dx / c
grid = Grid(shape=(Nx + 1,), extent=(L,))
u = TimeFunction(name='u', grid=grid, time_order=1, space_order=2)
x_coords = np.linspace(0, L, Nx + 1)
u.data[0, :] = I(x_coords)
courant = Constant(name='C', value=C)
# Lax-Wendroff: u - (C/2)*dx*u.dx + (C²/2)*dx²*u.dx2
# u.dx = centered first derivative
# u.dx2 = centered second derivative
stencil = u - 0.5*courant*dx*u.dx + 0.5*courant**2*dx**2*u.dx2
update = Eq(u.forward, stencil)
op = Operator([update])
# ... time stepping loopHere we use Devito’s built-in derivative operators:
u.dxcomputes the centered first derivative \((u_{i+1} - u_{i-1})/(2\Delta x)\)u.dx2computes the centered second derivative \((u_{i+1} - 2u_i + u_{i-1})/\Delta x^2\)
4.20.6 Lax-Friedrichs Scheme Implementation
The Lax-Friedrichs scheme is simpler but more diffusive:
\[ u^{n+1}_i = \frac{1}{2}(u^n_{i+1} + u^n_{i-1}) - \frac{C}{2}(u^n_{i+1} - u^n_{i-1}) \]
def solve_advection_lax_friedrichs(L, c, Nx, T, C, I):
"""Lax-Friedrichs scheme for 1D advection."""
dx = L / Nx
dt = C * dx / c
grid = Grid(shape=(Nx + 1,), extent=(L,))
x_dim, = grid.dimensions
u = TimeFunction(name='u', grid=grid, time_order=1, space_order=1)
x_coords = np.linspace(0, L, Nx + 1)
u.data[0, :] = I(x_coords)
courant = Constant(name='C', value=C)
# Neighbor values
u_plus = u.subs(x_dim, x_dim + x_dim.spacing)
u_minus = u.subs(x_dim, x_dim - x_dim.spacing)
# Lax-Friedrichs stencil
stencil = 0.5*(u_plus + u_minus) - 0.5*courant*(u_plus - u_minus)
update = Eq(u.forward, stencil)
op = Operator([update])
# ... time stepping loop4.20.7 Periodic Boundary Conditions
For advection problems, periodic boundary conditions are often useful to study wave propagation without boundary effects:
t_dim = grid.stepping_dim
# Periodic BC: u[0] wraps to u[Nx], u[Nx] wraps to u[0]
bc_left = Eq(u[t_dim + 1, 0], u[t_dim, Nx])
bc_right = Eq(u[t_dim + 1, Nx], u[t_dim + 1, 0])
op = Operator([update, bc_left, bc_right])4.20.8 Using the Solvers
The complete solver implementation in src/advec/advec1D_devito.py provides convenient interfaces:
from src.advec import (
solve_advection_upwind,
solve_advection_lax_wendroff,
solve_advection_lax_friedrichs,
exact_advection_periodic
)
import numpy as np
# Define initial condition
def I(x):
return np.exp(-0.5*((x - 0.25)/0.05)**2)
# Solve with upwind scheme
result = solve_advection_upwind(
L=1.0, c=1.0, Nx=100, T=0.5, C=0.8, I=I,
periodic_bc=True
)
# Compare with exact solution
u_exact = exact_advection_periodic(result.x, result.t, c=1.0, L=1.0, I=I)
error = np.max(np.abs(result.u - u_exact))
print(f"Max error: {error:.6f}")4.20.9 Scheme Comparison
The three schemes exhibit different numerical behaviors:
import matplotlib.pyplot as plt
from src.advec import (
solve_advection_upwind,
solve_advection_lax_wendroff,
solve_advection_lax_friedrichs,
exact_advection_periodic
)
import numpy as np
def I(x):
return np.exp(-0.5*((x - 0.25)/0.05)**2)
L, c, Nx, T, C = 1.0, 1.0, 50, 0.5, 0.8
# Solve with all three schemes
r_upwind = solve_advection_upwind(L, c, Nx, T, C, I, periodic_bc=True)
r_lw = solve_advection_lax_wendroff(L, c, Nx, T, C, I, periodic_bc=True)
r_lf = solve_advection_lax_friedrichs(L, c, Nx, T, C, I, periodic_bc=True)
# Exact solution
u_exact = exact_advection_periodic(r_upwind.x, r_upwind.t, c, L, I)
plt.figure(figsize=(10, 6))
plt.plot(r_upwind.x, u_exact, 'k-', lw=2, label='Exact')
plt.plot(r_upwind.x, r_upwind.u, 'b--', label='Upwind')
plt.plot(r_lw.x, r_lw.u, 'r-.', label='Lax-Wendroff')
plt.plot(r_lf.x, r_lf.u, 'g:', label='Lax-Friedrichs')
plt.legend()
plt.xlabel('x')
plt.ylabel('u')
plt.title(f'Advection: Nx={Nx}, C={C}, T={T}')
plt.savefig('advec_scheme_comparison.pdf')The Lax-Wendroff scheme typically preserves the wave amplitude better but may show small oscillations. The upwind and Lax-Friedrichs schemes are more diffusive, causing the wave to spread and reduce in amplitude.
4.20.10 Convergence Testing
We can verify the convergence rates of the schemes:
from src.advec import (
solve_advection_upwind,
solve_advection_lax_wendroff,
convergence_test_advection
)
# Test upwind (expect 1st order)
sizes, errors, rate = convergence_test_advection(
solve_advection_upwind,
grid_sizes=[25, 50, 100, 200],
T=0.25, C=0.8
)
print(f"Upwind convergence rate: {rate:.2f}") # ~1.0
# Test Lax-Wendroff (expect 2nd order)
sizes, errors, rate = convergence_test_advection(
solve_advection_lax_wendroff,
grid_sizes=[25, 50, 100, 200],
T=0.25, C=0.8
)
print(f"Lax-Wendroff convergence rate: {rate:.2f}") # ~2.04.20.11 Key Takeaways
Upwind differencing is essential for stable advection schemes—centered differences in space are unconditionally unstable.
The Courant number \(C = c\Delta t/\Delta x\) controls stability; all schemes require \(C \leq 1\).
Trade-offs exist between accuracy and numerical diffusion:
- Upwind: Stable, 1st order, diffusive
- Lax-Wendroff: 2nd order, less diffusion, may have small oscillations
- Lax-Friedrichs: Very stable, very diffusive
Devito’s shifted indexing via
u.subs(x_dim, x_dim - x_dim.spacing)allows expressing upwind differences naturally.Periodic BCs are implemented by explicitly setting boundary equations that copy values from the opposite end of the domain.
4.21 Exercises: Advection with Devito
4.21.1 Exercise 1: Verify CFL Stability Condition
The upwind scheme requires \(C \leq 1\) for stability.
a) Run the upwind solver with \(C = 0.5\), \(C = 0.9\), and \(C = 1.0\) for \(T = 1.0\) with a Gaussian initial condition. Verify that all solutions remain bounded.
b) Try \(C = 1.01\) and observe what happens. How quickly does the instability grow?
c) For \(C = 1.0\) exactly, the upwind scheme should reproduce the exact solution (up to machine precision). Verify this numerically.
from src.advec import solve_advection_upwind, exact_advection_periodic
import numpy as np
def I(x):
return np.exp(-0.5*((x - 0.25)/0.05)**2)
# Part a: Stable Courant numbers
for C in [0.5, 0.9, 1.0]:
result = solve_advection_upwind(
L=1.0, c=1.0, Nx=100, T=1.0, C=C, I=I
)
print(f"C={C}: u in [{result.u.min():.4f}, {result.u.max():.4f}]")
# Part b: Slightly unstable
# This will raise ValueError since C > 1 violates stability
try:
result = solve_advection_upwind(
L=1.0, c=1.0, Nx=100, T=1.0, C=1.01, I=I
)
except ValueError as e:
print(f"Error: {e}")
# Part c: Exact at C=1
result = solve_advection_upwind(
L=1.0, c=1.0, Nx=100, T=0.5, C=1.0, I=I, periodic_bc=True
)
u_exact = exact_advection_periodic(result.x, result.t, 1.0, 1.0, I)
error = np.max(np.abs(result.u - u_exact))
print(f"Error at C=1: {error:.2e}") # Should be ~machine precision4.21.2 Exercise 2: Compare Numerical Diffusion
The upwind scheme introduces numerical diffusion that causes the wave amplitude to decrease over time.
a) Run all three schemes (upwind, Lax-Wendroff, Lax-Friedrichs) with \(C = 0.8\) for \(T = 2.0\) and track the maximum value of \(u\) over time.
b) Plot the amplitude decay for each scheme. Which scheme preserves the amplitude best?
c) For the Gaussian initial condition, measure the “width” of the pulse (e.g., the distance between points where \(u = 0.5 \max(u)\)) at \(T = 0\) and \(T = 2\). How much has each scheme spread the pulse?
from src.advec import (
solve_advection_upwind,
solve_advection_lax_wendroff,
solve_advection_lax_friedrichs
)
import numpy as np
import matplotlib.pyplot as plt
def I(x):
return np.exp(-0.5*((x - 0.25)/0.05)**2)
# Run all schemes with history
L, c, Nx, T, C = 1.0, 1.0, 100, 2.0, 0.8
r_up = solve_advection_upwind(L, c, Nx, T, C, I, save_history=True)
r_lw = solve_advection_lax_wendroff(L, c, Nx, T, C, I, save_history=True)
r_lf = solve_advection_lax_friedrichs(L, c, Nx, T, C, I, save_history=True)
# Part b: Track amplitude decay
max_up = [np.max(u) for u in r_up.u_history]
max_lw = [np.max(u) for u in r_lw.u_history]
max_lf = [np.max(u) for u in r_lf.u_history]
plt.figure()
plt.plot(r_up.t_history, max_up, 'b-', label='Upwind')
plt.plot(r_lw.t_history, max_lw, 'r--', label='Lax-Wendroff')
plt.plot(r_lf.t_history, max_lf, 'g-.', label='Lax-Friedrichs')
plt.axhline(1.0, color='k', linestyle=':', label='Exact')
plt.xlabel('Time')
plt.ylabel('Max amplitude')
plt.legend()
plt.title('Amplitude decay comparison')
plt.savefig('amplitude_decay.pdf')
# Part c: Measure pulse width at half-maximum
def half_width(u, x):
u_max = np.max(u)
half_max = 0.5 * u_max
above = np.where(u >= half_max)[0]
if len(above) > 0:
return x[above[-1]] - x[above[0]]
return 0
print("Initial width:", half_width(I(r_up.x), r_up.x))
print("Upwind width:", half_width(r_up.u, r_up.x))
print("Lax-Wendroff width:", half_width(r_lw.u, r_lw.x))
print("Lax-Friedrichs width:", half_width(r_lf.u, r_lf.x))4.21.3 Exercise 3: Convergence Rate Verification
Verify the theoretical convergence rates: - Upwind: 1st order - Lax-Wendroff: 2nd order - Lax-Friedrichs: 1st order
a) Use the convergence_test_advection function with grid sizes [25, 50, 100, 200, 400] and verify the rates.
b) Create a log-log plot of error vs grid size for all three schemes.
c) What happens to the convergence rate if you use a discontinuous initial condition (step function) instead of the smooth Gaussian?
from src.advec import (
solve_advection_upwind,
solve_advection_lax_wendroff,
solve_advection_lax_friedrichs,
convergence_test_advection
)
import numpy as np
import matplotlib.pyplot as plt
# Part a: Verify rates
grid_sizes = [25, 50, 100, 200, 400]
sizes_up, err_up, rate_up = convergence_test_advection(
solve_advection_upwind, grid_sizes, T=0.25, C=0.8
)
print(f"Upwind rate: {rate_up:.2f}")
sizes_lw, err_lw, rate_lw = convergence_test_advection(
solve_advection_lax_wendroff, grid_sizes, T=0.25, C=0.8
)
print(f"Lax-Wendroff rate: {rate_lw:.2f}")
sizes_lf, err_lf, rate_lf = convergence_test_advection(
solve_advection_lax_friedrichs, grid_sizes, T=0.25, C=0.8
)
print(f"Lax-Friedrichs rate: {rate_lf:.2f}")
# Part b: Log-log plot
plt.figure()
plt.loglog(sizes_up, err_up, 'b-o', label=f'Upwind (rate={rate_up:.2f})')
plt.loglog(sizes_lw, err_lw, 'r-s', label=f'Lax-Wendroff (rate={rate_lw:.2f})')
plt.loglog(sizes_lf, err_lf, 'g-^', label=f'Lax-Friedrichs (rate={rate_lf:.2f})')
# Reference slopes
h = np.array(sizes_up)
plt.loglog(h, err_up[0]*(h[0]/h), 'k--', alpha=0.5, label='O(h)')
plt.loglog(h, err_lw[0]*(h[0]/h)**2, 'k:', alpha=0.5, label='O(h²)')
plt.xlabel('Grid points')
plt.ylabel('L2 Error')
plt.legend()
plt.title('Convergence comparison')
plt.gca().invert_xaxis()
plt.savefig('convergence_advec.pdf')4.21.4 Exercise 4: Step Function Advection
A step (Heaviside) function is a challenging test case for advection schemes because of the discontinuity.
a) Advect a step function from \(x = 0.25\) using all three schemes with \(C = 0.8\) and \(\Delta x = 0.01\). Compare the results at \(T = 0.5\).
b) The Lax-Wendroff scheme may show oscillations near the discontinuity (Gibbs phenomenon). Observe and document this behavior.
c) How does the upwind scheme handle the step? Does it preserve the sharp transition?
from src.advec import (
solve_advection_upwind,
solve_advection_lax_wendroff,
solve_advection_lax_friedrichs,
step_initial_condition,
exact_advection_periodic
)
import numpy as np
import matplotlib.pyplot as plt
def I(x):
return np.where(x < 0.25, 1.0, 0.0)
L, c, Nx, T, C = 1.0, 1.0, 100, 0.5, 0.8
r_up = solve_advection_upwind(L, c, Nx, T, C, I, periodic_bc=True)
r_lw = solve_advection_lax_wendroff(L, c, Nx, T, C, I, periodic_bc=True)
r_lf = solve_advection_lax_friedrichs(L, c, Nx, T, C, I, periodic_bc=True)
u_exact = exact_advection_periodic(r_up.x, r_up.t, c, L, I)
plt.figure(figsize=(10, 6))
plt.plot(r_up.x, u_exact, 'k-', lw=2, label='Exact')
plt.plot(r_up.x, r_up.u, 'b--', label='Upwind')
plt.plot(r_lw.x, r_lw.u, 'r-.', label='Lax-Wendroff')
plt.plot(r_lf.x, r_lf.u, 'g:', label='Lax-Friedrichs')
plt.legend()
plt.xlabel('x')
plt.ylabel('u')
plt.title('Step function advection')
plt.ylim(-0.2, 1.3)
plt.savefig('step_advection.pdf')
# Note Lax-Wendroff oscillations near discontinuity4.21.5 Exercise 5: Long-Time Integration
With periodic boundary conditions, a wave should return to its starting position after traveling one domain length.
a) Advect a Gaussian pulse for \(T = 1.0\) (one complete cycle with \(c = 1\), \(L = 1\)) and compare the final solution to the initial condition.
b) Run for \(T = 10.0\) (10 cycles) and measure how much the amplitude has decayed for each scheme.
c) For each scheme, estimate after how many cycles the peak amplitude drops to 50% of its initial value.
from src.advec import (
solve_advection_upwind,
solve_advection_lax_wendroff,
solve_advection_lax_friedrichs
)
import numpy as np
def I(x):
return np.exp(-0.5*((x - 0.25)/0.05)**2)
L, c, Nx, C = 1.0, 1.0, 100, 0.8
# Part a: One cycle
for T in [1.0, 10.0]:
r_up = solve_advection_upwind(L, c, Nx, T, C, I, periodic_bc=True)
r_lw = solve_advection_lax_wendroff(L, c, Nx, T, C, I, periodic_bc=True)
r_lf = solve_advection_lax_friedrichs(L, c, Nx, T, C, I, periodic_bc=True)
print(f"\nT = {T} ({int(T)} cycles):")
print(f" Upwind: max = {r_up.u.max():.4f}")
print(f" Lax-Wendroff: max = {r_lw.u.max():.4f}")
print(f" Lax-Friedrichs: max = {r_lf.u.max():.4f}")
# Part c: Find half-life
def find_halflife(solver_func, L, c, Nx, C, I, max_cycles=100):
for n in range(1, max_cycles + 1):
T = float(n)
result = solver_func(L, c, Nx, T, C, I, periodic_bc=True)
if result.u.max() < 0.5:
return n
return max_cycles
print("\nCycles to 50% amplitude:")
print(f" Upwind: {find_halflife(solve_advection_upwind, L, c, Nx, C, I)}")
print(f" Lax-Wendroff: {find_halflife(solve_advection_lax_wendroff, L, c, Nx, C, I)}")
print(f" Lax-Friedrichs: {find_halflife(solve_advection_lax_friedrichs, L, c, Nx, C, I)}")4.21.6 Exercise 6: Effect of Courant Number
The Courant number \(C\) affects both stability and accuracy.
a) For the upwind scheme, run with \(C = 0.2\), \(0.5\), \(0.8\), and \(1.0\) for \(T = 1.0\). Plot the final solutions on the same figure.
b) Which value of \(C\) gives the best accuracy? Why?
c) Measure the L2 error for each \(C\) value and create a plot of error vs. \(C\).
from src.advec import solve_advection_upwind, exact_advection_periodic
import numpy as np
import matplotlib.pyplot as plt
def I(x):
return np.exp(-0.5*((x - 0.25)/0.05)**2)
L, c, Nx, T = 1.0, 1.0, 100, 1.0
C_values = [0.2, 0.5, 0.8, 1.0]
plt.figure(figsize=(10, 6))
errors = []
for C in C_values:
result = solve_advection_upwind(L, c, Nx, T, C, I, periodic_bc=True)
plt.plot(result.x, result.u, label=f'C={C}')
u_exact = exact_advection_periodic(result.x, result.t, c, L, I)
dx = L / Nx
error = np.sqrt(dx * np.sum((result.u - u_exact)**2))
errors.append(error)
# Add exact solution
u_exact = exact_advection_periodic(result.x, T, c, L, I)
plt.plot(result.x, u_exact, 'k--', lw=2, label='Exact')
plt.legend()
plt.xlabel('x')
plt.ylabel('u')
plt.title('Effect of Courant number on upwind scheme')
plt.savefig('courant_effect.pdf')
# Error vs C
plt.figure()
plt.plot(C_values, errors, 'bo-')
plt.xlabel('Courant number C')
plt.ylabel('L2 Error')
plt.title('Error vs Courant number (Upwind)')
plt.savefig('error_vs_courant.pdf')
# C=1 gives exact solution for upwind
print("Errors:", dict(zip(C_values, errors)))4.21.7 Exercise 7: Variable Velocity Field
Modify the upwind solver to handle a spatially varying velocity \(c(x)\).
a) Implement an upwind scheme for: \[ \frac{\partial u}{\partial t} + c(x)\frac{\partial u}{\partial x} = 0 \] where the local Courant number varies: \(C_i = c(x_i)\Delta t/\Delta x\).
b) Test with \(c(x) = 1 + 0.5\sin(2\pi x)\) and observe how the wave stretches and compresses as it moves through regions of different velocity.
from devito import Grid, TimeFunction, Function, Eq, Operator, Constant
import numpy as np
import matplotlib.pyplot as plt
def solve_advection_variable_c(L, c_func, Nx, T, dt, I):
"""Upwind scheme with spatially varying velocity."""
grid = Grid(shape=(Nx + 1,), extent=(L,))
x_dim, = grid.dimensions
t_dim = grid.stepping_dim
u = TimeFunction(name='u', grid=grid, time_order=1, space_order=1)
c = Function(name='c', grid=grid)
x_coords = np.linspace(0, L, Nx + 1)
u.data[0, :] = I(x_coords)
c.data[:] = c_func(x_coords)
dx = L / Nx
dt_const = Constant(name='dt', value=dt)
dx_const = Constant(name='dx', value=dx)
# Local Courant number: C_i = c_i * dt / dx
# Upwind: u^{n+1} = u - (c*dt/dx)*(u - u[x-dx])
u_minus = u.subs(x_dim, x_dim - x_dim.spacing)
stencil = u - (c * dt_const / dx_const) * (u - u_minus)
update = Eq(u.forward, stencil)
# Periodic BCs
bc_left = Eq(u[t_dim + 1, 0], u[t_dim, Nx])
bc_right = Eq(u[t_dim + 1, Nx], u[t_dim + 1, 0])
op = Operator([update, bc_left, bc_right])
Nt = int(round(T / dt))
for n in range(Nt):
op.apply(time_m=n, time_M=n, dt=dt)
return u.data[Nt % 2, :].copy(), x_coords
# Test with variable velocity
def I(x):
return np.exp(-0.5*((x - 0.25)/0.05)**2)
def c_var(x):
return 1.0 + 0.5*np.sin(2*np.pi*x)
L, Nx = 1.0, 200
dx = L / Nx
c_max = 1.5 # max of c(x)
dt = 0.5 * dx / c_max # ensure CFL < 1 everywhere
u_final, x = solve_advection_variable_c(L, c_var, Nx, T=1.0, dt=dt, I=I)
plt.figure(figsize=(10, 6))
plt.plot(x, I(x), 'k--', label='Initial')
plt.plot(x, u_final, 'b-', label='Final (T=1)')
plt.xlabel('x')
plt.ylabel('u')
plt.legend()
plt.title('Advection with variable velocity c(x) = 1 + 0.5*sin(2*pi*x)')
plt.savefig('variable_velocity.pdf')4.21.8 Exercise 8: Advection-Diffusion Equation
Combine advection and diffusion: \[ \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = \nu\frac{\partial^2 u}{\partial x^2} \]
a) Implement a solver using upwind for advection and centered differences for diffusion.
b) Compare the behavior for \(\nu = 0\) (pure advection), \(\nu = 0.01\) (advection-dominated), and \(\nu = 0.1\) (diffusion-dominated).
c) What is the stability condition when both advection and diffusion are present?
from devito import Grid, TimeFunction, Eq, Operator, Constant
import numpy as np
import matplotlib.pyplot as plt
def solve_advec_diff(L, c, nu, Nx, T, C, I):
"""Advection-diffusion with upwind advection + centered diffusion."""
dx = L / Nx
# Stability requires both CFL and diffusion conditions
dt_adv = C * dx / c if c > 0 else np.inf
dt_diff = 0.4 * dx**2 / nu if nu > 0 else np.inf
dt = min(dt_adv, dt_diff)
grid = Grid(shape=(Nx + 1,), extent=(L,))
x_dim, = grid.dimensions
t_dim = grid.stepping_dim
u = TimeFunction(name='u', grid=grid, time_order=1, space_order=2)
x_coords = np.linspace(0, L, Nx + 1)
u.data[0, :] = I(x_coords)
C_const = Constant(name='C', value=c * dt / dx)
F_const = Constant(name='F', value=nu * dt / dx**2)
# Upwind advection + centered diffusion
u_minus = u.subs(x_dim, x_dim - x_dim.spacing)
advection = C_const * (u - u_minus)
diffusion = F_const * dx**2 * u.dx2
stencil = u - advection + diffusion
update = Eq(u.forward, stencil)
# Periodic BCs
bc_left = Eq(u[t_dim + 1, 0], u[t_dim, Nx])
bc_right = Eq(u[t_dim + 1, Nx], u[t_dim + 1, 0])
op = Operator([update, bc_left, bc_right])
Nt = int(round(T / dt))
for n in range(Nt):
op.apply(time_m=n, time_M=n, dt=dt)
return u.data[Nt % 2, :].copy(), x_coords
def I(x):
return np.exp(-0.5*((x - 0.25)/0.05)**2)
L, c, Nx, T, C = 1.0, 1.0, 100, 0.5, 0.8
plt.figure(figsize=(10, 6))
for nu, style in [(0.0, 'b-'), (0.01, 'r--'), (0.1, 'g-.')]:
u, x = solve_advec_diff(L, c, nu, Nx, T, C, I)
plt.plot(x, u, style, label=f'nu={nu}')
plt.plot(x, I(x), 'k:', lw=2, label='Initial')
plt.legend()
plt.xlabel('x')
plt.ylabel('u')
plt.title('Advection-diffusion equation')
plt.savefig('advec_diff.pdf')4.21.9 Exercise 9: Cosine Hat Initial Condition
The “cosine hat” is a smoother alternative to the step function: \[ I(x) = \begin{cases} \cos\left(\frac{5\pi}{L}(x - L/10)\right) & \text{if } x < L/5 \\ 0 & \text{otherwise} \end{cases} \]
a) Implement this initial condition and advect it using all three schemes.
b) Compare the behavior at the sharp cutoff (\(x = L/5\)) between schemes.
c) Does the Lax-Wendroff scheme show oscillations for this smoother discontinuity?
from src.advec import (
solve_advection_upwind,
solve_advection_lax_wendroff,
solve_advection_lax_friedrichs
)
import numpy as np
import matplotlib.pyplot as plt
def cosine_hat(x, L=1.0):
"""Cosine hat initial condition."""
result = np.zeros_like(x)
mask = x < L/5
result[mask] = np.cos(5*np.pi/L * (x[mask] - L/10))
return result
def I(x):
return cosine_hat(x, L=1.0)
L, c, Nx, T, C = 1.0, 1.0, 100, 0.5, 0.8
r_up = solve_advection_upwind(L, c, Nx, T, C, I, periodic_bc=True)
r_lw = solve_advection_lax_wendroff(L, c, Nx, T, C, I, periodic_bc=True)
r_lf = solve_advection_lax_friedrichs(L, c, Nx, T, C, I, periodic_bc=True)
plt.figure(figsize=(10, 6))
plt.plot(r_up.x, I(r_up.x - c*T), 'k-', lw=2, label='Exact')
plt.plot(r_up.x, r_up.u, 'b--', label='Upwind')
plt.plot(r_lw.x, r_lw.u, 'r-.', label='Lax-Wendroff')
plt.plot(r_lf.x, r_lf.u, 'g:', label='Lax-Friedrichs')
plt.legend()
plt.xlabel('x')
plt.ylabel('u')
plt.title('Cosine hat advection')
plt.savefig('cosinehat.pdf')4.21.10 Exercise 10: Implement Leapfrog Scheme
The leapfrog scheme uses a two-level time difference: \[ \frac{u^{n+1}_i - u^{n-1}_i}{2\Delta t} + c\frac{u^n_{i+1} - u^n_{i-1}}{2\Delta x} = 0 \]
This is a three-time-level scheme requiring special initialization for \(u^1\).
a) Implement the leapfrog scheme using Devito with time_order=2.
b) Use the upwind scheme to compute \(u^1\) from \(u^0\), then switch to leapfrog.
c) Compare the leapfrog scheme’s dispersion properties with Lax-Wendroff. Does leapfrog preserve amplitude better?
from devito import Grid, TimeFunction, Eq, Operator, Constant
import numpy as np
import matplotlib.pyplot as plt
def solve_advection_leapfrog(L, c, Nx, T, C, I):
"""Leapfrog scheme with upwind initialization."""
dx = L / Nx
dt = C * dx / c
grid = Grid(shape=(Nx + 1,), extent=(L,))
x_dim, = grid.dimensions
t_dim = grid.stepping_dim
# time_order=2 gives access to u, u.forward, u.backward
u = TimeFunction(name='u', grid=grid, time_order=2, space_order=1)
x_coords = np.linspace(0, L, Nx + 1)
# Set u^0
u.data[0, :] = I(x_coords)
# First step: use upwind to get u^1
courant = Constant(name='C', value=C)
u_minus = u.subs(x_dim, x_dim - x_dim.spacing)
upwind_stencil = u - courant * (u - u_minus)
# For leapfrog: u^{n+1} = u^{n-1} - C*(u^n_{i+1} - u^n_{i-1})
u_plus_x = u.subs(x_dim, x_dim + x_dim.spacing)
u_minus_x = u.subs(x_dim, x_dim - x_dim.spacing)
leapfrog_stencil = u.backward - courant * (u_plus_x - u_minus_x)
# Periodic BCs
bc_left = Eq(u[t_dim + 1, 0], u[t_dim, Nx])
bc_right = Eq(u[t_dim + 1, Nx], u[t_dim + 1, 0])
# First step with upwind
update_first = Eq(u.forward, upwind_stencil)
op_first = Operator([update_first, bc_left, bc_right])
op_first.apply(time_m=0, time_M=0, dt=dt)
# Leapfrog for remaining steps
update_lf = Eq(u.forward, leapfrog_stencil)
op_lf = Operator([update_lf, bc_left, bc_right])
Nt = int(round(T / dt))
for n in range(1, Nt):
op_lf.apply(time_m=n, time_M=n, dt=dt)
return u.data[Nt % 3, :].copy(), x_coords
def I(x):
return np.exp(-0.5*((x - 0.25)/0.05)**2)
L, c, Nx, T, C = 1.0, 1.0, 100, 1.0, 0.8
u_lf, x = solve_advection_leapfrog(L, c, Nx, T, C, I)
# Compare with Lax-Wendroff
from src.advec import solve_advection_lax_wendroff, exact_advection_periodic
r_lw = solve_advection_lax_wendroff(L, c, Nx, T, C, I, periodic_bc=True)
u_exact = exact_advection_periodic(x, T, c, L, I)
plt.figure()
plt.plot(x, u_exact, 'k-', lw=2, label='Exact')
plt.plot(x, u_lf, 'b--', label='Leapfrog')
plt.plot(x, r_lw.u, 'r-.', label='Lax-Wendroff')
plt.legend()
plt.xlabel('x')
plt.ylabel('u')
plt.title('Leapfrog vs Lax-Wendroff')
plt.savefig('leapfrog.pdf')
print(f"Leapfrog max: {u_lf.max():.4f}")
print(f"Lax-Wendroff max: {r_lw.u.max():.4f}")