Advection and Diffusion Equation¶

Jupyter notebook to illustrate numerical integration, using upwind and central finite differences.

Philipp Schlatter, January 2022

Initialisation of the libraries and graphics:

In [14]:

%matplotlib notebook 
# possible options: notebook, inline or widget
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
import matplotlib.animation
from math import pi
params = {'legend.fontsize': 12,
          'legend.loc':'best',
          'figure.figsize': (8,5),
          'lines.markerfacecolor':'none',
          'axes.labelsize': 12,
          'axes.titlesize': 12,
          'xtick.labelsize':12,
          'ytick.labelsize':12,
          'grid.alpha':0.6}
pylab.rcParams.update(params)

Perliminaries¶

Definition of our initial condition, domain size and discretisation points. Here, we take either a wave packet, simply defined as

$$f_1(x) = \sin(2 x) \cdot e^{-x^2/20} \ ,$$

or a square wave (top-hat signal, $f_2(x)$). A third option is a single blop centred around $x=0$,

$$f_3(x) = \cos(x/5) \cdot e^{-x^2/10} \ .$$

In [15]:

L=40
T=40
nx=200
dx=L/nx    
x=np.linspace(-L/2,L/2,num=nx,endpoint=False)

In [16]:

# Wavepacket
f1 = lambda  x: np.sin(2*x)*np.exp(-x**2/20)
# Square wave (top-hat function)
f2 = lambda x: np.array([1 if (t>=-18 and t<=-16) else 0 for t in x])
# one blob
f3 = lambda x: np.cos(x/5)*np.exp(-x**2/10)

# set initial condition (f1, f2 or f3)
u0=f1(x)

In [17]:

# Plot initial condition    
plt.figure()
plt.title('Initial condition')
plt.plot(x,u0,'.-b',lw=1,label='initial condition')
plt.xlabel(r'$x$')
plt.ylabel(r'$f^\prime(x)$')
plt.grid()
plt.legend(loc='upper right')
plt.show()

Advection Equation¶

In a first step, we consider the advection equation, i.e. the differential equation

$$\frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0 \ .$$

The solution can be obtained via the method of characteristics, and yields

$$u(x,t) = u(x-ct,0) = u_0(x-ct)$$

with the initial condition $u_0(x)$. In a periodic domain ($u(-L/2)=u(L/2)$), we do not need to provide addition inflow/outflow conditions (discussion related to type of partical differential equation, see course literature).

Compute the numerical solution, using an upwind scheme in space and time (FTBS - forward time backward space). This scheme is, as discussed later, stable for Courant number $$0<\sigma=c\frac{\Delta t}{\Delta x}<1 \, .$$ Note that $\sigma = 1$ corresponds to an exact integration.

In [18]:

c=1
sigma = 1
dt = sigma*dx/c

nt = int(T/dt)+1
t=np.linspace(0,T,num=nt)

u1=np.zeros((nx,nt))
for i in range(0,nx):
    u1[i,0] = u0[i]
    
# FTBS scheme for the advection equation
for j in range(1,nt):
    for i in range(0,nx):
        u1[i,j] = u1[i,j-1] - c*dt/dx*( u1[i,j-1]-u1[i-1,j-1])

In [20]:

#plt.rcParams["animation.html"] = "jshtml"
#plt.ion()

fig = plt.figure(figsize=(8,5))

def animate(j):
    plt.cla()
    plt.xlabel('$x$')
    plt.ylabel('$f^\prime(x)$')
    plt.grid()
    plt.plot(x,u0,'k--',lw=1,label='initial')
    plt.plot(x,u1[:,4*j],'r',label='advection')
    plt.xlim(-20,20)
    plt.ylim(-1.1,1.1)
    plt.title(f"time={t[4*j]:.2f}  $\sigma$={sigma:.2f}")  
    plt.legend(loc='upper right')
    
ani=matplotlib.animation.FuncAnimation(fig, animate, frames=int((t.size-1)/4+1), repeat=False,interval=2)
#writer = matplotlib.animation.writers['ffmpeg']
#writer = writer(fps=24)
#ani.save('out.mp4', writer=writer)

We can now also show our solution as a so-called space-time diagram, where we can see the complete solution $u(x,t)$.

In [21]:

plt.figure()
plt.pcolor(x,t,u1.T,shading='nearest')
plt.xlabel(r'$x$')
plt.ylabel(r'$t$')
plt.show()

Diffusion equation¶

In a second step, we consider the diffusion equation (heat equation), i.e.

$$ \frac{\partial u}{\partial t} = \nu \frac{\partial^2 u}{\partial x^2} \ .$$

Exponential growth is expected for negative viscosity, $\nu<0$, and the problem is not well-posed anymore.

Here we discretise the equation using explicit Euler in time and central differences in space, which is stable for

$$\beta = \nu \frac{\Delta t}{\Delta x^2} \leq 1/2 \, .$$

In [22]:

nu = 0.1 
beta = 0.45
dt = beta*(dx**2)/nu

T=40
nt = int(T/dt)+1
t=np.linspace(0,T,num=nt)

u2=np.zeros((nx,nt))
for i in range(0,nx):
    u2[i,0] = u0[i]
    
# FTCS scheme for the advection-diffusion equation
for j in range(1,nt):
    for i in range(0,nx):
        u2[i,j] = u2[i,j-1]   \
                 + nu*dt/(dx**2)*( u2[(i+1)%nx,j-1]-2*u2[i,j-1]+u2[i-1,j-1]  )

In [24]:

#plt.rcParams["animation.html"] = "jshtml"
#plt.ioff()

fig = plt.figure(figsize=(8,5))

def animate(j):
    plt.cla()
    plt.xlabel('$x$')
    plt.ylabel('$f^\prime(x)$')
    plt.grid()
    plt.plot(x,u0,'k--',lw=1,label='initial')
    plt.plot(x,u2[:,4*j],'b-',label='diffusion')
    plt.xlim(-20,20)
    plt.ylim(-1.1,1.1)
    plt.title(f"time={t[4*j]:.2f} beta={beta:.2f}")   
    plt.legend(loc='upper right')
    
ani=matplotlib.animation.FuncAnimation(fig, animate, frames=int((t.size-1)/4+1),repeat=False,interval=2)
#writer = matplotlib.animation.writers['ffmpeg']
#writer = writer(fps=24)
#ani.save('out.mp4', writer=writer)

Also here, a space-time diagram gives the full picture:

In [25]:

plt.figure()
plt.pcolor(x,t,u2.T,shading='nearest')
plt.xlabel(r'$x$')
plt.ylabel(r'$t$')
plt.show()

Advection-Diffusion equation¶

In a second step, we consider the advection-diffusion equation, i.e.

$$ \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = \nu \frac{\partial^2 u}{\partial x^2} \ .$$

The solution is obviously pure advection for $\nu=0$, and adds a certain amount of dissipation for $\nu>0$. Exponential growth is expected for $\nu<0$, and the problem is not well-posed anymore.

Here we discretise the equation using explicit Euler in time and central differences in space. For the coupled advection-diffusion equation one gets the following

$$\sigma^2\leq 2\beta \leq 1$$

which can be re-written as two inequalities:

$\beta \leq 1/2$: $\Delta t\leq \frac{1}{2\nu}\Delta x^2$ (as for diffusion equation)
$\sigma^2/\beta \leq 2$: $\Delta t\leq 2\nu/ c^2$ (independent of $\Delta x$)

In [26]:

c=1
nu = 0.1  

beta = 0.45

dt1 = beta*(dx**2)/nu
dt2 = 2*nu/c**2
dt = min(dt1,dt2)

nt = int(T/dt)+1
t=np.linspace(0,T,num=nt)

u3=np.zeros((nx,nt))
for i in range(0,nx):
    u3[i,0] = u0[i]
    
# FTCS scheme for the advection-diffusion equation
for j in range(1,nt):
    for i in range(0,nx):
        u3[i,j] = u3[i,j-1]   \
                 - c*dt/(2*dx)*( u3[(i+1)%nx,j-1]-u3[i-1,j-1] ) \
                 + nu*dt/(dx**2)*( u3[(i+1)%nx,j-1]-2*u3[i,j-1]+u3[i-1,j-1]  )

In [27]:

#plt.rcParams["animation.html"] = "jshtml"
#plt.ioff()

fig = plt.figure(figsize=(8,5))

def animate(j):
    plt.cla()
    plt.xlabel('$x$')
    plt.ylabel('$f^\prime(x)$')
    plt.grid()
    plt.plot(x,u0,'k--',lw=1,label='initial')
    plt.plot(x,u3[:,4*j],'g-',label='advection-diffusion')
    plt.xlim(-20,20)
    plt.ylim(-1.1,1.1)
    plt.title(f"time={t[4*j]:.2f}")   
    plt.legend(loc='upper right')
    
ani=matplotlib.animation.FuncAnimation(fig, animate, frames=int((t.size-1)/4+1), repeat=False,interval=2)
#writer = matplotlib.animation.writers['ffmpeg']
#writer = writer(fps=24)
#ani.save('out.mp4', writer=writer)

Also here, a space-time diagram is instructive:

In [28]:

plt.figure()
plt.pcolor(x,t,u3.T,shading='nearest')
plt.xlabel(r'$x$')
plt.ylabel(r'$t$')
plt.show()

Stability and truncation error¶

The stability for the advection and diffusion equations separately are for a forward-time-central space (FTCS) scheme:

$\sigma = c\frac{\Delta t}{\Delta x} \leq 1$ (note only necessary as based on CFL condition)
$\beta = \nu \frac{\Delta t}{\Delta x^2} \leq 1/2$ (necessary and sufficient)

However, for the coupled advection-diffusion equation one gets the following

$$\sigma^2\leq 2\beta \leq 1$$

which can be re-written as two inequalities:

$\beta \leq 1/2$: $\Delta t\leq \frac{1}{2\nu}\Delta x^2$ (as for diffusion equation)
$\sigma^2/\beta \leq 2$: $\Delta t\leq 2\nu/ c^2$ (independent of $\Delta x$)

One observes the following

for high viscosity the diffusive term is dominant: viscous time-step limit
for low viscosity: advection is stabilised by viscosity
for zero viscosity ($\nu=0$): maximum $\Delta t=0$, i.e. the scheme is unconditionally unstable (central discretisation of advection is always unstable)

The truncation error for the upwind scheme is $$T(x_j,t^n) = \frac{1}{2}c \Delta x u_{xx} - \frac{1}{2} \Delta t u_{tt} $$

and thus the numerical viscosity for the spatial scheme

$$ \nu_{num} = \frac{1}{2}c \Delta x \ , $$

and the numerical viscosity for the temporal scheme (Euler forward) is

$$ \nu_{num,t} = -\frac{1}{2} c^2 \Delta t \ .$$

For the advection-diffusion scheme we have no diffusive errors in space (i.e. $\nu_{num}=0$).

In [ ]: