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

Question

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:

Answer 1

I am unfamiliar with the CTU scheme, but if you are running into problems of numerical dissipation in advection problems, one easy approach might be to do higher order upwinding. That basically means that when you want to calculate the flux at a certain point, you adjust your derivative to be calculated in the "upwind" direction.

If you are on a structured grid like your images suggest, that might be an easy fix. The right/left discretisations for your gradient should keep simple form like this:

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

You might also want to check whether your time-stepping size might be too big. The relevant keyword here is the CFL-condition. Just check that the advection is not so strong, that the transported distance in the time dt exceedes your cell-size. Because if that happens then the discretisation will fail miserably.

hope that helps.

Answer 2

Adding onto @MPIchael good answer, I believe a linear upwind scheme (LUD) will work well. A first-order scheme is generally too diffusive as the discretized version tends to have a second-order term with artificial viscosity while schemes like lax wendroff are second-order though they are quite dispersive. Please refer to "Convection and Diffusion" chapter 5 of "Numerical heat and fluid flow" by Patankar for a detailed discussion of various schemes. For your reference, an oblique advection of a initially non-smooth variable using Lax-Wendroff is shown below. Hope this helps.

Answer 3

Since Wolfgang Bangerths comment below the question might easily be overlooked: The important thing to realize is that you want to solve the very well studied linear transport equation which is usually written as $$ \frac{\partial u(\boldsymbol{x}, t)}{\partial t} + \boldsymbol{a} \cdot \nabla u(\boldsymbol{x}, t) = 0$$ for which a rich literature of schemes exist, mostly the aforementioned finite volume schemes. The fact that most of them are tailored to solving Riemann problems while aiming to surpress the numerical diffusion makes them the recommended choice for this kind of problems.