1D Burger's Equation

PDE and Initial Condition:

\begin{aligned} \frac{\partial u}{\partial t} + u\,\frac{\partial u}{\partial x} &= \alpha\frac{\partial^{2} u}{\partial x^{2}} \quad&&\mathrm{;}\quad x\in[0,L],\quad t\in[0,T]\\\\ u(x,0) &= u_{0}(x) \quad&&\mathrm{;}\quad x\in[0,L] \end{aligned}

Upwind Scheme

Under the linearization $u = c + \delta u$, we obtain the Von-Neumann stability condition:

$$ r = \dfrac{c\Delta t}{\Delta x} + \dfrac{\alpha\Delta t}{\Delta x^2} < 1 $$

Defining $r_1 = \dfrac{\Delta t}{\Delta x}$ and $r_2 = \dfrac{\alpha \Delta t}{\Delta x^2}$

The discretized form of the equation is:

$$ u^{n+1}_{i} = u^{n}_{i} - \big(r_2 + r_1 \cdot u^{n}_{i}\big) u^{n}_{i} + \big(r_2 + r_1\cdot u^{n}_{i} \big) u^{n}_{i-1} + \big(r_2 \big) u^{n}_{i+1} \quad;\quad i = 1,2,\ldots, N_{x} $$

• Introduce "ghost" boundary points $u^{n}_{0}$ and $u^{n}_{N_{x} + 1}$ in order to implement boundary conditions.

  • Periodic BC: Update the solution as if the spatial domain wraps around itself, which requires $u^{n}_{1} = u^{n}_{N_{x}}$.

    • The "ghost" point to the left of $u^{n}_{1}$ must be the domain point to the left of $u^{n}_{N_{x}}$:

    $$u^{n}_{0} = u^{n}_{N_{x}-1}$$

    • The "ghost" point to the right of $u^{n}_{N_{x}}$ must be the domain point to the right of $u^{n}_{1}$:

    $$u^{n}_{N_{x}+1} = u^{n}_{2}$$

  • Fixed BC: Update the solution as if the spatial domain does not admit non-zero solutions at particular nodes. For the boundaries, this requires:

    • Setting $u^{n}_{1} = u^{n}_{N_{x}} = 0$, and doing the same for the ghost points.

import time, sys
import numpy as np
import scipy.constants as spc
from scipy import signal as sg
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.animation

from IPython.display import display, Math, HTML
from numpy.linalg import inv
from matplotlib.animation import FuncAnimation
from scipy.ndimage import convolve, generate_binary_structure

plt.rcParams.update(
    {
        'mathtext.fontset':'cm', 
        'font.family':'serif', 
        'font.sans-serif':'Times New Roman', 
        'figure.figsize':[9,4], 
        'figure.titlesize':22,
        "figure.dpi":90, 
        'savefig.dpi':90,
        'axes.titlesize':20, 
        'axes.labelsize':18, 
        'axes.titley': 1.0, 
        'axes.titlepad': 5.0,
        'axes.edgecolor':'black', 
        'axes.grid': False, 
        'grid.alpha': .5,
        'xtick.labelsize':14,'ytick.labelsize':14,
        'xtick.major.size':6,'ytick.major.size':6,
        'xtick.major.width':1.25, 'ytick.major.width':1.25, 
        'xtick.direction':'inout','ytick.direction':'inout',
        'xtick.top':False, 'ytick.right':False,
        'legend.title_fontsize':14, 'legend.fontsize':14,
        'legend.borderaxespad': 1, 'legend.borderpad': 0.5,
        'legend.framealpha': 1,
        'legend.handleheight': 0.5, 'legend.handlelength': 2.0, 'legend.handletextpad':0.5,
        'legend.labelspacing': 0.25,
        'legend.fancybox':False,
        'legend.edgecolor': '0',
        'legend.frameon': True,
        'legend.markerscale': 1.25,
        'animation.embed_limit':2**128,
        'animation.html': 'jshtml'
    }
)
    
pi = spc.pi
%matplotlib inline
%config InlineBackend.figure_format = 'retina'

def animate1d(u, x, t, fps = 24, frames = 144):
    t_steps = int(len(t) - 1)
    i_frames = np.arange(0, t_steps, int(t_steps/frames))
    fig, ax = plt.subplots(1,1, figsize = (6,6), dpi = 300)
    line, = ax.plot([], [], 'k-', lw = 2)
    ax.axis('off')
    ax.set_xlim(x[0], x[-1])
    ax.set_ylim(u[:,1:-1].min()-1,u[:,1:-1].max()+1)
    ax.set_aspect('auto')
    def animate(i):
        line.set_data(x[1:-1],u[i,1:-1])
        return line,
    plt.close()
    anim = FuncAnimation(fig, animate, interval = int(1e3/fps), frames = i_frames, repeat = False)
    return anim

def periodic_boundary_1D(f):
    f.T[0] = f.T[-2]
    f.T[-1] = f.T[2]
    return

def fixed_boundary_1D(f):
    f.T[0] = 0
    f.T[-1] = 0
    f.T[-2] = 0
    f.T[2] = 0
    return

def set_boundary_1D(u, boundary = 'fix'):
    if boundary == 'fix':
        fixed_boundary_1D(u)
    else:
        periodic_boundary_1D(u)
    return

def burgers_1D(u0, x, alpha, r, T, boundary = 'periodic'):
    dx = x[1]-x[0]  
    c = np.abs(np.max(u0))
    dt = r/(c/dx + alpha/dx/dx)/2
    Nt = int(T/dt)
    t = np.arange(Nt)*dt
    u = np.zeros(t.shape + x.shape)
    u[0] = u0

    advect = (dt/dx)
    diffuse = alpha*(dt/dx)/dx
    
    for n in range(Nt-1):
        
        set_boundary_1D(u, boundary)
        
        Un   = np.copy(u[n,1:-1]) # u^{n}_{j}
        Un_w = np.copy(u[n,0:-2]) # u^{n}_{j-1}
        Un_e = np.copy(u[n,2:  ]) # u^{n}_{j+1}
        
        dUn = diffuse*(Un_e + Un_w - 2*Un) - advect*(Un - Un_w)*Un
        
        u[n+1,1:-1] = u[n,1:-1] + dUn
    
    u_fft = np.fft.fft2(u[:,1:-1])
    
    return t, u, u_fft

def smooth_rect(x, x0, y0, s, w, A):
    """
        's' controls how steep the smoothed-discontinuity is.
        'w' controls the width of the pulse.
        'A' is the amplitude of the pulse.
    """
    L = x[-1]
    y = A*(1. - 2*np.exp(-s*w/L) * np.cosh(s*(x-x0)/L))
    return y0 + y * (y >= 0)

L = 10
Nx = 256
dx = L/(Nx-1)
x = np.arange(-1, Nx+1)*dx
u0 = 4  + smooth_rect(x, x0 = 0.1*L, y0 = 0, s = Nx, w = L/20, A = 10) \
        + smooth_rect(x, x0 = 0.3*L, y0 = 0, s = Nx, w = L/20, A = 5)

t, u, u_fft = burgers_1D(u0, x, alpha = 1e-2, r = .9, T = 4, boundary = 'periodic')
print('Creating plots . . .')
time.sleep(2)
fig, ax = plt.subplots(figsize = (8,5))
mp = ax.pcolormesh(x[1:-1], t, u[:,1:-1], cmap = 'tab20c', shading = 'auto')
ax.set_aspect('auto')
fig.colorbar(mp, fraction = 1/20)
plt.show()
time.sleep(1)
print('Creating animation . . .')
time.sleep(2)
ani = animate1d(u, x, t, fps = 24, frames = int(24*5))
# writergif = mpl.animation.PillowWriter(fps = 24)
# ani.save('shock_pbc.gif', writer = writergif)
display(ani)

Creating plots . . .

Creating animation . . .

Once Loop Reflect