0
$\begingroup$

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: enter image description here

$\endgroup$
  • $\begingroup$ Have you done a literature search for "numerical methods for advection dominated problems"? $\endgroup$ – Wolfgang Bangerth Mar 20 at 16:10
  • $\begingroup$ 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? $\endgroup$ – Peanutlex Mar 20 at 16:39
  • $\begingroup$ You should look for things such as finite volume methods, or stabilized finite element methods (e.g., the streamline upwind Petrov Galerkin method, SUPG). $\endgroup$ – Wolfgang Bangerth Mar 20 at 18:32
  • $\begingroup$ 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? $\endgroup$ – Peanutlex Mar 20 at 19:19
  • $\begingroup$ No. What you are doing is called "operator splitting", but I don't know what effect that would have in your case. $\endgroup$ – Wolfgang Bangerth Mar 20 at 21:11
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.