Wave (Wave equations chapter) and diffusion (Diffusion equations chapter)
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), and compact difference notation.

**Remark on terminology.**

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. This is a an excellent resource for further studies on the topic of advection PDEs, with emphasis on gener alizations to real geophysical problems. The book by Fletcher also has a good overview of methods for advection and convection problems.

# One-dimensional time-dependent advection equations¶

We consider the pure advection model

In (1), \(v\) is a given parameter, typically reflecting the transport velocity of a quantity \(u\) with a flow. There is only one boundary condition (3) since the spatial derivative is only first order in the PDE (1). 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 (1) if we remove the boundary condition and consider a process on the infinite domain \((-\infty, \infty)\). The solution is simply

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) is that

if \(x_i=i\Delta x\) and \(t_n=n\Delta t\) are points in a uniform mesh. We see this relation from

provided \(v = \Delta x/\Delta t\). So, whenever we see a scheme that collapses to

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.

## Simplest scheme: forward in time, centered in space¶

### Method¶

A first attempt to solve a PDE like (1) 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:

Written out, we see that this expression implies that

with \(C\) as the Courant number

### Implementation¶

A solver function for our scheme goes as follows.

```
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from devito import Grid, Eq, solve, TimeFunction, Operator
```

```
# %load -s solver_FECS src-advec/advec1D.py
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
# Make sure dx and dt are compatible with x and t
dx = float(x[1] - x[0])
dt = float(t[1] - t[0])
C = v*dt/dx
grid = Grid(shape=(Nx+1,), extent=(L,))
t_s=grid.time_dim
u = TimeFunction(name='u', grid=grid, space_order=2, save=Nt+1)
pde = u.dtr + v*u.dxc
stencil = solve(pde, u.forward)
eq = Eq(u.forward, stencil)
# Set initial condition u(x,0) = I(x)
u.data[1, :] = [I(xi) for xi in x]
# Insert boundary condition
bc = [Eq(u[t_s+1, 0], U0)]
op = Operator([eq] + bc)
op.apply(dt=dt, x_m=1, x_M=Nx-1)
if user_action is not None:
for n in range(0, Nt + 1):
user_action(u.data[n], x, t, n)
```

### 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):

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

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.savefig()
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()
```

### 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.

## Analysis of the scheme¶

It is easy to show that a typical Fourier component

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

with

Note that \(|A_\text{e}| \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 \(A_\text{e}\) (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

Inserting the numerical component in the scheme,

and making use of this equation results in

which implies

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.

## Leapfrog in time, centered differences in space¶

### Method¶

Another explicit scheme is to do a “leapfrog” jump over \(2\Delta t\) in time and combine it with central differences in space:

which results in the updating formula

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`

.

### 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
...
if scheme == 'FE':
u = TimeFunction(name='u', grid=grid, space_order=2, save=Nt+1)
pde = u.dtr + v*u.dxc
...
elif scheme == 'LF':
u = ...
# Use some scheme for the first step
pde0 = ...
...
# Move to the LF scheme after the first step
pde = u.dtc + v*u.dxc
...
```

### 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, which is in fact identical to the exact solution (!).

Exact solution obtained by Leapfrog scheme with \(\Delta t = 0.1\) and \(C=1\).

### Running more test cases¶

We can run two types of initial conditions for \(C=0.8\): one very smooth with a Gaussian function (Figure) and one with a discontinuity in the first derivative (Figure). 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.

Advection of a Gaussian function with a leapfrog scheme and \(C=0.8\), \(\Delta t = 0.001\) (left) and \(\Delta t=0.01\) (right).

Advection of the Gaussian function with a leapfrog scheme, using \(C=0.8\) and \(\Delta t = 0.01\) can be seen in a movie file. Alternatively, with \(\Delta t = 0.001\), we get this movie file.

Advection of half a cosine function with a leapfrog scheme and \(C=0.8\), \(\Delta t = 0.001\) (left) and \(\Delta t=0.01\) (right).

Advection of the cosine hat function with a leapfrog scheme, using \(C=0.8\) and \(\Delta t = 0.01\) can be seen in a movie file. Alternatively, with \(\Delta t = 0.001\), we get this movie file.

### Analysis¶

We can perform a Fourier analysis again. Inserting the numerical Fourier component in the Leapfrog scheme, we get

and

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\).

**Stability.**

For all the working schemes to be presented in this chapter, we get the stability condition \(C\leq 1\):

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

and is to be compared to the exact amplification factor

the section Analysis of dispersion relations compares numerical amplification factors of many schemes with the exact expression.

## 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

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,

written out as

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. This is easy to see since \(C=1\) gives the property (6). However, any \(C<1\) gives a significant reduction in the amplitude of the solution, which is a purely numerical effect, see this Figure and this Figure. Experiments show, however, that reducing \(\Delta t\) or \(\Delta x\), while keeping \(C\) reduces the error.

Advection of a Gaussian function with a forward in time, upwind in space scheme and \(C=0.8\), \(\Delta t = 0.01\) (left) and \(\Delta t=0.001\) (right).

Advection of the Gaussian function with a forward in time, upwind in space scheme, using \(C=0.8\) and \(\Delta t = 0.01\) can be seen in a movie file. Alternatively, with \(\Delta t = 0.005\), we get this movie file.

Advection of half a cosine function with a forward in time, upwind in space scheme and \(C=0.8\), \(\Delta t = 0.001\) (left) and \(\Delta t=0.01\) (right).

Advection of the cosine hat function with a forward in time, upwind in space scheme, using \(C=0.8\) and \(\Delta t = 0.01\) can be seen in a movie file. Alternatively, with \(\Delta t = 0.001\), we get this movie file.

The amplification factor can be computed using this formula.

which means

For \(C<1\) there is, unfortunately, non-physical damping of discrete Fourier components, giving rise to reduced amplitude of \(u^n_i\) as in this Figure and this Figure. The damping seen in these figures is quite severe. Stability requires \(C\leq 1\).

**Interpretation of upwind difference as artificial diffusion.**

One can interpret the upwind difference as extra, artificial diffusion in the equation. Solving

by a forward difference in time and centered differences in space,

actually gives the upwind scheme (10) 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:

## 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

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[t_s, 0]`

is updated as
`u[t_s, 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[t_s, 0]=u[t_s, Nx]`

at a new time level before computing `u[t_s+1, 1]`

. This ensures
that `u[t_s+1, 1]`

becomes right and at the next time level `u[t_s+1, 0]`

at the current
time level is correctly updated.
For the Leapfrog scheme we must update `u[t_s+1, 0]`

and `u[t_s+1, Nx]`

using the scheme:

```
pbc = [Eq(u[t_s+1, 0], u[t_s-1, 0] - C*(u[t_s, 1] - u[t_s, Nx-1]))]
pbc += [Eq(u[t_s+1, Nx], u[t_s+1, 0])]
...
```

## Implementation¶

### Test condition¶

Analytically, we can show that the integral in space under the \(u(x,t)\) curve is constant:

as long as \(u(0)=u(L)=0\). We can therefore use the property

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.data[n][0] + 0.5*u.data[n][Nx] + np.sum(u.data[n][1:Nx]))
```

if `u`

is a `TimeFunction`

with the `save`

parameter set to \(Nx+1\) and `n`

indicates a current timestep.

### The code¶

An appropriate `solver`

function for multiple schemes may go as shown
below.

```
# %load -s solver src-advec/advec1D.py
def solver(I, U0, v, L, dt, C, T, user_action=None,
scheme='FE', periodic_bc=True):
Nt = int(round(T/np.float64(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
# Make sure dx and dt are compatible with x and t
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))
integral = np.zeros(Nt+1)
grid = Grid(shape=(Nx+1,), extent=(L,), dtype=np.float64)
t_s=grid.time_dim
def u(to=1, so=1):
u = TimeFunction(name='u', grid=grid, time_order=to, space_order=so, save=Nt+1)
return u
if scheme == 'FE':
u = u(so=2)
pde = u.dtr + v*u.dxc
pbc = [Eq(u[t_s+1, 0], u[t_s, 0] - 0.5*C*(u[t_s, 1] - u[t_s, Nx]))]
pbc += [Eq(u[t_s+1, Nx], u[t_s+1, 0])]
elif scheme == 'LF':
# Use UP scheme for the first timestep
u = u(to=2, so=2)
pde0 = u.dtr(fd_order=1) + v*u.dxl(fd_order=1)
stencil0 = solve(pde0, u.forward)
eq0 = Eq(u.forward, stencil0).subs(t_s, 0)
pbc0 = [Eq(u[t_s, 0], u[t_s, Nx]).subs(t_s, 0)]
# Now continue with LF scheme
pde = u.dtc + v*u.dxc
pbc = [Eq(u[t_s+1, 0], u[t_s-1, 0] - C*(u[t_s, 1] - u[t_s, Nx-1]))]
pbc += [Eq(u[t_s+1, Nx], u[t_s+1, 0])]
elif scheme == 'UP':
u = u()
pde = u.dtr + v*u.dxl
pbc = [Eq(u[t_s, 0], u[t_s, Nx])]
elif scheme == 'LW':
u = u(so=2)
pde = u.dtr + v*u.dxc - 0.5*dt*v**2*u.dx2
pbc = [Eq(u[t_s+1, 0], u[t_s, 0] - 0.5*C*(u[t_s, 1] - u[t_s, Nx-1]) + \
0.5*C**2*(u[t_s, 1] - 2*u[t_s, 0] + u[t_s, Nx-1]))]
pbc += [Eq(u[t_s+1, Nx], u[t_s+1, 0])]
else:
raise ValueError('scheme="%s" not implemented' % scheme)
stencil = solve(pde, u.forward)
eq = Eq(u.forward, stencil)
bc_init = [Eq(u[t_s+1, 0], U0).subs(t_s, 0)]
# Set initial condition u(x,0) = I(x)
u.data[0, :] = [I(xi) for xi in x]
# Compute the integral under the curve
integral[0] = dx*(0.5*u.data[0][0] + 0.5*u.data[0][Nx] + np.sum(u.data[0][1:Nx]))
if user_action is not None:
user_action(u.data[0], x, t, 0)
bc = [Eq(u[t_s+1, 0], U0)]
if scheme == 'LF':
op = Operator((pbc0 if periodic_bc else []) + [eq0] + (bc_init if not periodic_bc else []) \
+ (pbc if periodic_bc else []) + [eq] + (bc if not periodic_bc else []))
else:
op = Operator(bc_init + (pbc if periodic_bc else []) + [eq] + (bc if not periodic_bc else []))
op.apply(dt=dt, x_m=1, x_M=Nx if scheme == 'UP' else Nx-1)
for n in range(1, Nt+1):
# Compute the integral under the curve
integral[n] = dx*(0.5*u.data[n][0] + 0.5*u.data[n][Nx] + np.sum(u.data[n][1:Nx]))
if user_action is not None:
user_action(u.data[n], x, t, n)
print('I:', integral[n])
return integral
```

### 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 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.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
# Define video formats and libraries
codecs = dict(flv='flv', mp4='libx264', webm='libvpx',
ogg='libtheora')
# Remove video files
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)
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)
# Finish figure(2)
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()
# Make videos from figure(1) animation files
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`

.

## 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:

It can be nice to include the Backward Euler scheme too, via the \(\theta\)-rule,

When \(\theta\) is different from zero, this gives rise to an *implicit* scheme,

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).

The elements on the diagonal in the matrix become:

On the subdiagonal and superdiagonal we have

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

The dispersion relation follows from inserting \(u^n_q = A^ne^{ikx}\) and using this formula for the spatial differences:

Crank-Nicolson in time, centered in space, Gaussian profile, \(C=0.8\), \(\Delta t = 0.01\) (left) and \(\Delta t=0.005\) (right).

Backward-Euler in time, centered in space, half a cosine profile, \(C=0.8\), \(\Delta t = 0.01\) (left) and \(\Delta t=0.005\) (right).

```
from IPython.display import HTML
_s = """
<div>
<video loop controls width='640' height='365' preload='none'>
<source src='https://raw.githubusercontent.com/hplgit/fdm-book/master/doc/pub/book/html/mov-advec/gaussian/CN/C08_dt0005/movie.mp4' type='video/mp4; codecs="avc1.42E01E, mp4a.40.2"'>
<source src='https://raw.githubusercontent.com/hplgit/fdm-book/master/doc/pub/book/html/mov-advec/gaussian/CN/C08_dt0005/movie.webm' type='video/webm; codecs="vp8, vorbis"'>
<source src='https://raw.githubusercontent.com/hplgit/fdm-book/master/doc/pub/book/html/mov-advec/gaussian/CN/C08_dt0005/movie.ogg' type='video/ogg; codecs="theora, vorbis"'>
</video>
</div>
Crank-Nicolson in time, centered in space, $C=0.8$, $\Delta t = 0.005$. <!-- \label{advec:1D:CN:mov:C08:dt2} -->
<!-- Issue warning if in a Safari browser -->
<script language="javascript">
if (!!(window.safari)) {
document.write("<div style=\"width: 95%%; padding: 10px; border: 1px solid #100; border-radius: 4px;\"><p><font color=\"red\">The above movie will not play in Safari - use Chrome, Firefox, or Opera.</font></p></div>")}
</script>
"""
HTML(_s)
```

```
_s = """
<div>
<video loop controls width='640' height='365' preload='none'>
<source src='https://raw.githubusercontent.com/hplgit/fdm-book/master/doc/pub/book/html/mov-advec/cosinehat/BE/C_08_dt005.mp4' type='video/mp4; codecs="avc1.42E01E, mp4a.40.2"'>
<source src='https://raw.githubusercontent.com/hplgit/fdm-book/master/doc/pub/book/html/mov-advec/cosinehat/BE/C_08_dt005.webm' type='video/webm; codecs="vp8, vorbis"'>
<source src='https://raw.githubusercontent.com/hplgit/fdm-book/master/doc/pub/book/html/mov-advec/cosinehat/BE/C_08_dt005.ogg' type='video/ogg; codecs="theora, vorbis"'>
</video>
</div>
Backward-Euler in time, centered in space, $C=0.8$, $\Delta t = 0.005$. <!-- \label{advec:1D:BE:mov:C08:dt2} -->
<!-- Issue warning if in a Safari browser -->
<script language="javascript">
if (!!(window.safari)) {
document.write("<div style=\"width: 95%%; padding: 10px; border: 1px solid #100; border-radius: 4px;\"><p><font color=\"red\">The above movie will not play in Safari - use Chrome, Firefox, or Opera.</font></p></div>")}
</script>
"""
HTML(_s)
```

This figure 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 this figure 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.

## 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 \(\mathcal{O}{\Delta t^2} + \mathcal{O}{\Delta x^2}\).

Let us follow the recipe. First we have the three-term Taylor polynomial,

From the PDE we have that temporal derivatives can be substituted by spatial derivatives:

and furthermore,

Inserted in the Taylor polynomial formula, we get

To obtain second-order accuracy in space we now use central differences:

or written out,

This is the explicit Lax-Wendroff scheme.

**Lax-Wendroff works because of artificial viscosity.**

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 \(\frac{1}{2}\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

This means that \(|A|=1\) and also that we have an exact solution if \(C=1\)!

## 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

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}\):

for

Now,

so

An appropriate dimensionless quantity to plot is the scaled wave velocity \(c/v\):

Figures from this up to this 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 |

Dispersion relations for \(C=1\).

Dispersion relations for \(C=1\).

Dispersion relations for \(C=0.8\).

Dispersion relations for \(C=0.8\).

Dispersion relations for \(C=0.5\).

Dispersion relations for \(C=0.5\).

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 this figure, 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 this figure, 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.

# One-dimensional stationary advection-diffusion equation¶

Now we pay attention to a physical process where advection (or convection) is in balance with diffusion:

For simplicity, we assume \(v\) and \(\alpha\) 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

Equations of the form (11) or
(12) 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.

## A simple model problem¶

We consider (11) 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

Inserted in the governing equation we get

Dropping the bars is common. We can then simplify to

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

The forthcoming plots illustrate this function for various values of \(\epsilon\).

## A centered finite difference scheme¶

The most obvious idea to solve (13) is to apply centered differences:

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

for \(i=1,\ldots,N_x-1\)

The right-hand side of the linear system is zero except \(b_{N_x}=1\).

This figure 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 this figure. 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 1: Analyze 1D stationary convection-diffusion problem develops a precise criterion: \(u\) is oscillation-free if

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!

Comparison of exact and numerical solution for \(\epsilon =0.1\) and \(N_x=20,40\) with centered differences.

Comparison of exact and numerical solution for \(\epsilon =0.01\) and \(N_x=20,40\) with centered differences.

**Solver.**

A suitable solver for doing the experiments is presented below.

```
# TODO: IMPLEMENT THIS IN DEVITO
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
# Make sure dx and dt are compatible with x and t
dx = x[1] - x[0]
u = np.zeros(Nx+1)
# Representation of sparse matrix and right-hand side
diagonal = np.zeros(Nx+1)
lower = np.zeros(Nx)
upper = np.zeros(Nx)
b = np.zeros(Nx+1)
# Precompute sparse matrix (scipy format)
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
# Insert boundary conditions
upper[0] = 0
lower[-1] = 0
diagonal[0] = diagonal[-1] = 1
b[-1] = 1.0
# Set up sparse matrix and solve
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, x
```

## 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

we use the one-sided *upwind* difference

in case \(v>0\). For \(v<0\) we set

On compact operator notation form, our upwind scheme can be expressed as

provided \(v>0\) (and \(\epsilon > 0\)).

We write out the equations and implement them as shown in the program in the section A centered finite difference scheme. The results appear in this figure and this figure: no more oscillations!

Comparison of exact and numerical solution for \(\epsilon =0.1\) and \(N_x=20,40\) with upwind difference.

Comparison of exact and numerical solution for \(\epsilon =0.01\) and \(N_x=20,40\) with upwind difference.

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 1: Analyze 1D stationary convection-diffusion problem. Moreover, the thicker layer (seemingly larger diffusion) can be understood by doing Exercise 2: Interpret upwind difference as artificial diffusion.

**Exact solution for this model problem.**

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.

# Time-dependent convection-diffusion equations¶

Now it is time to combine time-dependency, convection (advection) and diffusion into one equation:

## 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

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.

## 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?

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.

## 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:

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

# 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 the Applications section of the Diffusion equations chapter,
but change the assumption from a solid to fluid medium.
This gives rise to the extra advection or convection term \(\boldsymbol{v}\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

## Transport of a substance¶

The diffusion of a substance in the section Diffusion of a substance 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 \(\boldsymbol{v}\) such that the effective flux now consists of two and not only one component as in this equation:

Inserted in the equation \(\nabla\cdot\boldsymbol{q} = 0\) we get the extra advection term \(\nabla\cdot (\boldsymbol{v}_\xi)\). Very often we deal with incompressible flows, \(\nabla\cdot\boldsymbol{v} = 0\) such that the advective term becomes \(\boldsymbol{v}\cdot\nabla c\). The mass transport equation for a substance then reads

## Transport of heat in fluids¶

The derivation of the heat equation in the section Heat conduction 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 this equation,

and more terms if work by stresses is also included, where

\(\boldsymbol{v}\) being the velocity of the fluid. The convective term \(\boldsymbol{v}\cdot\nabla e\) must therefore be added to the governing equation, resulting typically in

where \(f\) is some external heating inside the medium.

# Exercises¶

## Exercise 1: Analyze 1D stationary convection-diffusion problem¶

Explain the observations in the numerical experiments from the sections A centered finite difference scheme and Remedy: upwind finite difference scheme by finding exact numerical solutions.

**Hint.**
The difference equations allow solutions on the form \(A^i\), where
\(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.

Filename: `twopt_BVP_analysis1`

.

## Exercise 2: 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 \(\mathcal{O}{\Delta x}\).
Filename: `twopt_BVP_analysis2`

.