Burgers with MacCormack

Conside the 1D Burgers equation:

\begin{equation} \frac{\partial u}{\partial t}=-u\frac{\partial u}{\partial x} \end{equation}

We want to represent this in conservative forms so that we can better deal with potential shocks, which gives us:

\begin{equation} \frac{\partial u}{\partial t}=-\frac{\partial }{\partial x}\left(\frac{u^2}{2}\right) \end{equation}

We can also write this as:

\begin{equation} \frac{\partial u}{\partial t}=-\frac{\partial F}{\partial x} \end{equation}

If we take $F=\frac{u^2}{2}$. Start with a Heaviside function (a sep function) with the following values:

\begin{equation} u(x,0)=\left\{\begin{array}{cc} 1 & 0\leq x<2 \\ 0 & 2\leq x\leq 4 \end{array}\right. \end{equation}

import numpy as np
import matplotlib.pyplot as plt
from numba import njit
from matplotlib.animation import FuncAnimation
from jupyterthemes import jtplot
jtplot.style(theme='monokai')

Create the initial conditions for $u(x,0)$.

def IC(Nx,Nt,x):
    u=np.zeros((Nx,Nt)) #Create the matrix
    mask=np.where(x<2) #Create a mask where x is less than 2
    #Now, create a loop to put the initial values.
    for i in range(len(mask[-1])):
        u[i,0]=1
    return u

Now, you can use the MacCormack method to solve the PDE. The MacCormack scheme is:

\begin{equation} \large\begin{array}{cc} u_i^*=u_i^t-\frac{\Delta t}{\Delta x}(F_{i+1}^t+F_i^t) & \text{Predictor} \\
u_i^{t+1}=\frac{1}{2}\left(u_{i}^t+u_i^*-\frac{\Delta t}{\Delta x}(F_i^*-F_{i-1}^*)\right) & \text{Corrector} \end{array} \end{equation}

Remember that

\begin{equation} F=\frac{u^2}{2} \end{equation}

You can create two functions to solve this system

def F(u):
    return u**2/2

def u_star(u,k,i,j):
    return u[i,j]-k*(F(u[i+1,j])-F(u[i,j]))

def solution(Nx,Nt,x,dx,dt):
    u=IC(Nx,Nt,x) #Use the function to create the matrix and put the IC
    k=dt/dx       #For simplicity
    for j in range(Nt-1):
        for i in range(Nx-1):
            u[i,j+1]=0.5*(u[i,j]+u_star(u,k,i,j)-k*(F(u_star(u,k,i,j))-F(u_star(u,k,i-1,j))))
        u[0,j+1]=1
    return u

Use $N_x=81$, $N_t=70$ and $\Delta t=0.5\Delta x$.

Nx=81
Nt=70
x=np.linspace(0,4,Nx)
dx=4/(Nx-1)
dt=0.5*dx

u=solution(Nx,Nt,x,dx,dt)

Create an animate plot.

fig,axs=plt.subplots(1,figsize=(12,5))
line,=axs.plot([],[],color='teal',marker='o',markevery=5)
axs.set_xlim(x.min(),x.max())
axs.set_ylim(u.min()-0.5,u.max()+0.5)
axs.grid(False)

def animate(i):
    global x,u
    #axs.clear_collections()
    line.set_data(x,u[:,i])
    return line,

anim=FuncAnimation(fig,animate,frames=Nt,interval=60,blit=True)
plt.close()
anim.save('burger.mp4',fps=20)

The oscillations at the front of the shock are characteristic of all 2nd-order methods, but there are ways to remove them. One option is numerical damping.

A common damping term used with MacCormack scheme is: \begin{equation} \epsilon(u_{i+1}^t-2u_i^t+u_{i-1}^t) \end{equation}

I’ll add this term to the predictor step so i get this:

\begin{equation} u_i^*=u_i^t-\frac{\Delta t}{\Delta x}(F_{i+1}^t+F_i^t)+\epsilon(u_{i+1}^t-2u_i^t+u_{i-1}^t) \end{equation}

But I have a problem when $i=0$, because i’m gonna get $u_{-1}^{n}$. and that value is not part of my grid. So what I can do? Well, suppose that you have:

\begin{equation} \large\begin{array}{ccc} \frac{\partial u_0^t}{\partial t}=0 & \text{ then $\rightarrow$ } & \frac{\partial u_0^t}{\partial t}=\frac{u_{0}^t-u_{-1}^t}{\Delta t}=0 \end{array} \end{equation}

So, it’s easy to see that I can have $u_{-1}^t=u_0^{t}$. So, for the firs value, I have:

\begin{equation} u_0^*=u_0^t-\frac{\Delta t}{\Delta x}(F_{1}^t+F_0^t)+\epsilon(u_{1}^t-u_0^t) \end{equation}

def u_star(u,k,i,j):
    global epsilon
    if (i==0):
        return u[i,j]-k*(F(u[i+1,j])-F(u[i,j]))+epsilon*(u[i+1,j]-u[i,j])
    else:
        return u[i,j]-k*(F(u[i+1,j])-F(u[i,j]))+epsilon*(u[i+1,j]-2*u[i,j]+u[i-1,j])

epsilon=0.5
u=solution(Nx,Nt,x,dx,dt)

anim=FuncAnimation(fig,animate,frames=Nt,interval=60,blit=True)

anim.save('burger_no.mp4',fps=20)

That’s all for this post! You can check the jupyter notebook from this problem here.