# How to make a less diffusive code to solve 2D advection equation?

I would like to solve the following differential equation numerically in 2D, $$\frac{\partial z^-}{\partial t}+(\vec{B}\cdot\vec{\nabla})z^-=0,$$ see Wikipedia if you are curious about what the symbols mean, where $$\vec{B}=(B_x,B_y,0).$$ I have attached some code which can solve this to first order accuracy in python by following the corner transport upstream (CTU) algorithm outlined here. The problem is the code is too diffusive for my purposes. Does anyone know of an algorithm I could follow to make the code a higher order of accuracy while still remaining stable? Also, it would nice if the algorithm was an upstream algorithm i.e. each cell update only required information from cells in the upstream direction as it is easier to impose open boundary conditions if this is the case.

Here is the python code:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import os

nx = 64
x_min = -2
x_max = 2
dx = (x_max - x_min) / (nx - 1)
x = np.linspace(x_min, x_max, nx)

ny = 64
y_min = -2
y_max = 2
dy = (y_max - y_min) / (ny - 1)
y = np.linspace(y_min, y_max, ny)

bx = np.zeros((nx,ny)) + 1
by = np.zeros((nx,ny)) + 1

t_max = 1
dt = dx * dy / np.max(abs(bx) * dy + abs(by) * dx) # CFL condition
nt = int(t_max / dt)
t = np.linspace(0, t_max, nt)

zm = np.zeros((nx, ny, nt))
for i in range(nx):
for j in range(ny):
r = np.sqrt(x[i] ** 2 + y[j] ** 2)
if r < 1:
zm[i,j,0] = 2 * np.cos(np.pi * r / 2) ** 2

print(nt)
zm1 = np.zeros((nx, ny))
for n in range(nt - 1):

if n % 100 == 0: print('n =', n)

cx = bx[1:,1:] * dt / dx
cy = by[1:,1:] * dt / dy
zm1[1:,1:]    = (1 - cx) * zm[1:,1:,n] + cx * zm[0:-1,1:,n]
zm[1:,1:,n+1] = (1 - cy) * zm1[1:,1:]  + cy * zm1[1:,0:-1]

zm[0,:,n+1]  = 0
zm[:,0,n+1]  = 0

print('Computation finished, now generating figures...')

os.makedirs('Figures/2d/zm', exist_ok = True)

X = np.zeros((nx,ny))
Y = np.zeros((nx,ny))
for j in range(ny):
X[:,j] = x
for i in range(nx):
Y[i,:] = y
i = -1
for n in range(nt - 1):
# if n % 10 == 0:
i = i + 1
fig = plt.figure()
ax = fig.gca(projection = '3d')
ax.plot_surface(X, Y, zm[:,:,n])
plt.title('t = ' + str(t[n]))
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('zp')
ax.set_zlim(-1, 1)
plt.savefig('Figures/2d/zm/' + "{0:0=4d}".format(i) + '.png')
plt.close(fig)


My attempt at Beam-Warming (see comments):

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import os

nx = 68
x_min = 0
x_max = 2
dx = (x_max - x_min) / (nx - 1)
x = np.linspace(x_min - dx, x_max + dx, nx) # Add ghost cells

ny = 68
y_min = -2
y_max = 0
dy = (y_max - y_min) / (ny - 1)
y = np.linspace(y_min - dy, y_max + dy, ny)

bx = np.zeros((nx,ny)) + 1
by = np.zeros((nx,ny)) + 1

t_max = 10
dt = dx * dy / np.max(abs(bx) * dy + abs(by) * dx) # CFL condition
nt = int(t_max / dt)
t = np.linspace(0, t_max, nt)

zm = np.zeros((nx, ny, nt))
zp = np.zeros((nx, ny, nt))
for i in range(nx):
for j in range(ny):
r = np.sqrt((x[i] - 1) ** 2 + (y[j] + 1) ** 2)
if r < 1:
zm[i,j,0] = 2 * np.cos(np.pi * r / 2) ** 2

print(nt)
for n in range(nt - 1):

if n % 100 == 0: print('n =', n)

zm[2:,2:,n+1] = zm[2:,2:,n] \
- 0.5 * dt * bx[2:,2:] / dx * \
(3 * zm[2:,2:,n] - 4 * zm[1:-1,2:,n] + zm[0:-2,2:,n]) \
+ 0.5 * (dt * bx[2:,2:] / dx) ** 2 * \
(zm[2:,2:,n] - 2 * zm[1:-1,2:,n] + zm[0:-2,2:,n]) \
- 0.5 * dt * by[2:,2:] / dy * \
(3 * zm[2:,2:,n] - 4 * zm[2:,1:-1,n] + zm[2:,0:-2,n]) \
+ 0.5 * (dt * by[2:,2:] / dy) ** 2 * \
(zm[2:,2:,n] - 2 * zm[2:,1:-1,n] + zm[2:,0:-2,n])

zm[0,:,n+1]  = 0
zm[:,0,n+1]  = 0

print('Computation finished, now generating figures...')

os.makedirs('Figures/2d/zm', exist_ok = True)

X = np.zeros((nx,ny))
Y = np.zeros((nx,ny))
for j in range(ny):
X[:,j] = x
for i in range(nx):
Y[i,:] = y
i = -1
for n in range(nt - 1):
if n % 100 == 0:
i = i + 1
fig = plt.figure()
ax = fig.gca(projection = '3d')
ax.plot_surface(X, Y, zm[:,:,n])
plt.title('t = ' + str(t[n]))
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('zp')
ax.set_zlim(-1, 1)
plt.savefig('Figures/2d/zm/' + "{0:0=4d}".format(i) + '.png')
plt.close(fig)


Pictures of the instability:

• Have you done a literature search for "numerical methods for advection dominated problems"? – Wolfgang Bangerth Mar 20 '19 at 16:10
• I have come across a method called the Beam-Warming numerical scheme. However, it appears to be numerically unstable when I extend it to 2D. I have added my attempt to the main post. I was wondering if anyone knew under what conditions the Beam-Warming method is stable in 2D? – Peanutlex Mar 20 '19 at 16:39
• You should look for things such as finite volume methods, or stabilized finite element methods (e.g., the streamline upwind Petrov Galerkin method, SUPG). – Wolfgang Bangerth Mar 20 '19 at 18:32
• Okay, I will look into that. With my attempt at the Beam-Warming method I have noticed that if I update the field like so: $$z_{i,j}^*=z_{i,j}^n+B_x\frac{\partial z}{\partial x},$$ $$z_{i,j}^{n+1}=z_{i,j}^*+B_y\frac{\partial z^*}{\partial y},$$ then the instability does not appear. Do you have any idea why this is? Do you know if there any problems with using this scheme? – Peanutlex Mar 20 '19 at 19:19
• No. What you are doing is called "operator splitting", but I don't know what effect that would have in your case. – Wolfgang Bangerth Mar 20 '19 at 21:11

"left"-leaning discretization: $$\frac{du}{dx} = \frac{3u_i^n - 4u_{i-1}^n + u_{i-2}^n}{2\Delta x}$$ "right"-leaning discretsation: $$\frac{du}{dx} = \frac{-u_{i+2}^n + 4u_{i+1}^n - 3u_i^n}{2\Delta x}$$ Then at every point you check where the wind is blowing from and discretise in the suggested fashion. You will have to be a bit carefull at the border of your domains, as you might not have first and second neighbours in the direction of the border. That should somewhat alleviate your numerical dissipation problems, but some of it will persist.