09-burgers-dg.py

"""DG for Burgers with limiting."""
import sys
import numpy as np
from scipy.special import legendre
from scipy.interpolate import lagrange
import matplotlib.pyplot as plt
import modepy

#  An outline of this code is a follows:
#      1. problem defition
#      2. mesh setup
#      3. definition of f, the numerical flux, and the right hand side
#      4. the slop limiter
#      5. some random plotting routines
#      4. the main time stepping loop

# --------------------------------------------------------------------
# 1
# --------------------------------------------------------------------
# Define the problem
# Periodic boundary conditions
xa, xb = (0, 1)        # domain
nel = 100              # number of elements
p = 3                  # poly order
applylimiter = True    # should we slope limit?
T = 0.5                # length of time
ht = 0.001            # time step size


def uinit(z):
    """Set the initial condition.

    z : ndarray
        Assume z is p+1 x nel
    """
    uu = np.zeros_like(z)
    for i in range(uu.shape[1]):
        if np.all(z[:, i] <= 0.5):
            # make sure the entire cell is less than 0.5
            uu[:, i] = 1.0
    J = np.where(z <= 0.25)
    uu[J] = 4 * z[J]
    return uu


# --------------------------------------------------------------------
# 2
# --------------------------------------------------------------------
# Set the mesh points
x = np.linspace(xa, xb, nel + 1)
h = x[1] - x[0]

# Construct the local quadrature points
qx = modepy.quadrature.jacobi_gauss.legendre_gauss_lobatto_nodes(p)

# Full mesh, all nodes
xx = np.zeros((p + 1, nel))
for e in range(nel):
    xx[:, e] = x[e] + (h / 2) * (qx + 1)


# Vandermonde / mass
def nlegendre(qorder):
    """Compute the Legendre with unite L2."""
    return legendre(qorder) / np.sqrt(2/(2*qorder+1))


V = np.zeros((p+1, p+1))      # Vandermonde
dV = np.zeros((p+1, p+1))     # derivative Vandermonde
for q in range(p+1):
    poly = nlegendre(q)
    dpoly = poly.deriv()
    V[:, q] = poly(qx)
    dV[:, q] = dpoly(qx)
Vinv = np.linalg.inv(V)
M = Vinv.T @ Vinv            # Mass
Minv = (2/h) * V @ V.T       # Mass inverwse
S = M @ dV @ Vinv            # S (see text)
D = dV @ Vinv                # Derivative

n = (p+1)*nel                # total number of points

#   v-left-ext      v-right-ext
# ---|-------------|-----------|
#     ^left-int   ^right-int

# set the vectorized indices
I = np.arange(nel)           # index i
Im1 = np.roll(I, 1)          # index i-1
Ip1 = np.roll(I, -1)         # index i+1

# --------------------------------------------------------------------
# 3
# --------------------------------------------------------------------
# f: function in the conservation law u_t + f(u)_x = 0
#
# nflux: numerical flux
#
# L: right hand side


def f(u):
    """Calculate u**2 / 2."""
    return u**2 / 2


def df(u):
    """Calcualte derivative, u."""
    return u


def nflux(u_int, u_ext, n_int, n_ext):
    """Lax-Fredrichs flux.

    u : nel x 1
    """
    avgf = (f(u_int) + f(u_ext)) / 2
    jumpu = n_int * u_int + n_ext * u_ext
    C = np.max([np.abs(df(u_int)), np.abs(df(u_ext))])
    return avgf + (C/2) * jumpu


def L(u):  # noqa: N802
    """Right-hand side.

    dudt = M^-1 S.T @ f(u) - M^-1 f*(xright)[0,0,...,1] + f*(xleft)
    """
    flux = np.zeros_like(u)
    u_left_int = u[0, I]
    u_left_ext = u[p, Im1]
    u_right_int = u[p, I]
    u_right_ext = u[0, Ip1]
    flux[0, :] = nflux(u_left_int, u_left_ext, -1, 1)     # left
    flux[p, :] = -nflux(u_right_int, u_right_ext, 1, -1)  # right
    Lu = Minv @ S.T @ f(u) + Minv @ flux
    return Lu


# --------------------------------------------------------------------
# 4
# --------------------------------------------------------------------
def minmodv(a, b, c):
    """Vector-based minmod.

    minmod(a,b,c) = { sign(a) min(|a|, |b|, |c|) if sign(a)=sign(b)=sign(c)
                    {                            else 0
    """
    a = a.ravel()
    b = b.ravel()
    c = c.ravel()
    I = np.where(np.abs(np.sign(a) + np.sign(b) + np.sign(c))
                 == 3                          # all three values are the same
                 )[0]
    mm = np.zeros_like(a)
    mm[I] = np.sign(a[I]) * np.vstack((np.abs(a[I]),
                                       np.abs(b[I]),
                                       np.abs(c[I]))).min(axis=0)
    return mm


def slopelimiter(u):
    """Slopelimiters."""
    uhat = Vinv @ u                      # make it modal

    u1 = uhat.copy()
    u1[2:, :] = 0                        # strip to linear
    u1 = V @ u1                          # make it nodal
    ux = (u1[-1, :] - u1[0, :]) / h      # get the slope for each element
    ux = ux.ravel()

    uhat[1:, :] = 0                      # strip to constant (mean)
    ubar = (V @ uhat)[0, :].ravel()      # convert to nodal, and take the left value
    ubarm1 = ubar[Im1]                   # -1
    ubarp1 = ubar[Ip1]                   # +1

    newslope = minmodv(ux, (ubarp1 - ubar)/h, (ubar - ubarm1)/h)

    xbar = (xx[-1, :] + xx[0, :]) / 2    # midpoint of each element

    # make everything p+1 x nel
    bshape = (p+1, nel)                  # broadcast shape
    xcenter = xx - xbar                  # (broadcast) center x to the element
    ubar = np.broadcast_to(ubar, bshape)
    newslope = np.broadcast_to(newslope, bshape)

    # construct the new linear
    newu = ubar + xcenter * newslope
    return newu


def blanklimiter(u):
    """Apply no limiter."""
    return 1.0 * u


# --------------------------------------------------------------------
# 5
# --------------------------------------------------------------------
def plotnodal(xx, u, ax=None):
    """Plot all of the nodal values."""
    if not ax:
        _, ax = plt.subplots()
    ax.plot(xx.ravel(order='F'), u.ravel(order='F'), '-o', ms=3)


def plotmodal(xx, u, u2=None):
    """High fidelity modal."""
    c = Vinv @ u
    m = 40
    v = np.zeros((m, nel))
    zz = np.zeros((m, nel))
    z = np.linspace(-1, 1, m)
    for k in range(nel):
        zz[:, k] = xx[0, k] + (xx[p, k] - xx[0, k])*(z+1)/2
        for q in range(p+1):
            v[:, k] += c[q, k] * nlegendre(q)(z)

    zz = zz.ravel(order='F')
    plt.plot(zz, v.ravel(order='F'), '-o', ms=3)
    if u2 is not None:
        plt.plot(zz, u2(zz), color='tab:red', lw=1)
    plt.show()


# --------------------------------------------------------------------
# 6
# --------------------------------------------------------------------
# initial condition
u = uinit(xx)
u0 = u.copy()
nstep = int(np.ceil(T/ht))
step = 0

if applylimiter:
    limiter = slopelimiter
else:
    limiter = blanklimiter

plt.ion()
fig, ax = plt.subplots()
uline, = ax.plot(xx.ravel(order='F'), u0.ravel(order='F'), '-', color='tab:red', lw=1)
txt = ax.text(1.0 / 3.0, 0.8, 't=%g, i=%g' % (0.0, 0), fontsize=16)

while step * ht < T:

    # SSP 3 time stepping
    u1 = u + ht * L(u)
    u1 = limiter(u1)
    u2 = (1/4) * (3 * u + u1 + ht * L(u1))
    u2 = limiter(u2)
    u[:] = (1/3) * (u + 2 * u2 + 2 * ht * L(u2))
    u = limiter(u)

    step += 1

    #color = adjust_lightness('tab:blue', amount=amount[i])
    uline.set_ydata(u.ravel(order='F'))
    #uline.set_color('tab:blue')
    txt.set_text('t=%g, i=%g' % ((n + 1) * ht, n))

    fig.canvas.draw()
    fig.canvas.flush_events()
plt.show(block=True)