High Resolution Shock Capturing Methods: Lab 2

Directly relevant background material can be found in eg David Ketcheson's HyperPython notebooks.

Systems

In lab 1 the MUSCL scheme to solve scalar conservation laws was introduced. These methods work for discontinuous solutions, capturing the behaviour without introducing spurious oscillations. However, the convergence rate is not high (but still, better slow convergence than no convergence!).

In real life, we mostly want to solve conservation laws in

1. multiple dimensions, and
2. for systems of equations.

Working in multiple dimensions is relatively straightforward - the solution method can simply be applied one dimension at a time (although there exist more specialized methods that may have greater accuracy than this approach).

Moving from a scalar conservation law to a system of conservation laws is also relatively straightforward; the same reconstruction & Riemann problem solution approach will work. The main complication is the solution of the Riemann problem.

To check you can see how this would work for a simple problem, try the system

$$ \begin{equation} \partial_t {\bf q} + \partial_x \left( A {\bf q} \right) = {\bf 0}, \qquad A = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}, \end{equation} $$

which is just a system of (uncoupled) advection equations. Implement a MUSCL scheme to solve this on the domain $x \in [-1, 1]$ with periodic boundaries, using initial data $(\sin(\pi x), \sin(\pi x))^T$, checking the solution at various times and the convergence at $t=2$ (where the exact solution is simply the initial profile). Note that the exact solution to the Riemann problem given left and right states ${\bf q}^{(L,R)} = (q_1^{(L,R)}, q_2^{(L,R)})^T$ is

$$ \begin{equation} {\bf q}^{*} = \begin{pmatrix} q_1^{(L)} \\ q_2^{(R)} \end{pmatrix} \end{equation} $$

as the advection velocity for $q_1$ is positive whilst that for $q_2$ is negative.

Approximate Riemann Solvers

It may not always be possible, and frequently is not practical, to solve the Riemann problem exactly. When this occurs, it is typical to approximate the solution of the Riemann problem, sometimes by directly approximating the intercell flux $f_{i-1/2} = f_{i-1/2}(q^{(L)}_{i-1/2}, q^{(R)}_{i-1/2})$.

The simplest (and least accurate, or most diffusive) approximate solver is given by the Lax-Friedrichs, or HLL flux. Given the left and right states for the Riemann problem, and the time and space steps, this approximates the intercell flux as

$$ \begin{equation} f_{i-1/2} = \frac{1}{2} \left( f \left( q^{(L)}_{i-1/2} \right) + f \left( q^{(R)}_{i-1/2} \right) + \frac{\Delta x}{\Delta t} \left( q^{(L)}_{i-1/2} + - q^{(R)}_{i-1/2} \right) \right). \end{equation} $$

[Note: This is not exactly the Lax-Friedrichs or HLL fluxes, which have parameters controlling the wavespeed. Here I have simplified by setting all parameters to "safe" values.]

Implement this flux, and use it within your MUSCL scheme above. Compare the solution against the scheme using the exact solver.

Shallow water equations

A simple nonlinear system, itself a simplification of the Navier-Stokes equations, is the shallow water equations. In non-dimensional form these are

$$ \begin{equation} \partial_t \begin{pmatrix} \phi \\ \phi u \end{pmatrix} + \partial_x \begin{pmatrix} \phi u \\ \phi u^2 + \tfrac{1}{2} \phi^2 \end{pmatrix} = {\bf 0}. \end{equation} $$

It represents the flow of (typically incompressible) fluids in a shallow channel with a flat bottom topography; multi-dimensional versions that include topography as a source term are a standard way of simulating tsunamis, for example. Here $\phi$ is the geopotential (essentially the depth of the fluid) and $u$ is the velocity.

Use your MUSCL method with the Lax-Friedrichs solver to solve a dam-break problem. That is, solve the shallow water equations on the domain $x \in [-1, 1]$ with fixed boundary conditions (the solutions in the ghost cells can be copied from the previous solution), using the initial data

$$ \begin{equation} {\bf q}(x) = \begin{cases} {\bf q}^{(L)} = \begin{pmatrix} \phi^{(L)} \\ \phi^{(L)} u^{(L)} \end{pmatrix} = \begin{pmatrix} 3 \\ 0 \end{pmatrix} & x < 0, \\ {\bf q}^{(R)} = \begin{pmatrix} \phi^{(R)} \\ \phi^{(R)} u^{(R)} \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \end{pmatrix} & x > 0. \end{cases} \end{equation} $$

That is, $u^{(L,R)} = 0$, $\phi^{(L)} = 3$, $\phi^{(R)} = 1$. Plot the solution at $t = 0.4$, investigating how it depends on resolution.

The solution should consist of a continuous rarefaction wave propagating to the left and a discontinuous shock propagating to the right.

Characteristics

The solution of the Riemann problem is constant along characteristics. Analyzing the result requires calculating the characteristic structure. For this we need the Jacobian matrix $\partial {\bf f} / \partial {\bf q}$, and its eigenvectors and eigenvalues. In order for this information to be practically useful in a real algorithm it needs calculating analytically.

Calculate the Jacobian matrix for the shallow water equations. Also calculate its eigenvalues $\lambda^{\pm}$, and the left (${\bf l}^{\pm}$) and right (${\bf r}^{\pm}$) eigenvectors. The eigenvectors should be normalized so that ${\bf l}^i \cdot {\bf r}^j = \delta^{ij}$; that is, the eigenvectors are orthonormal.

You'll want to look at sympy, and the commands symbols, Matrix, diff, inv and eigenvects.

Calculate the characteristic speeds in the dam-break problem initial data. Use this to confirm that the wave structure is as expected.

Use the equation across the rarefaction curve

$$ \begin{equation} \frac{\partial}{\partial \xi} {\bf q} = \frac{{\bf r}^{\pm}}{{\bf r}^{\pm} \cdot \frac{\partial \lambda^{\pm}}{\partial {\bf q}}} \end{equation} $$

to find the solution across the rarefaction wave.