Introduction

What is a Riemann Problem?

Bernhard Riemann studied many problems in diverse mathematical fields, but for our purposes the term Riemann problem refers to a hyperbolic partial differential equation (PDE) together with special initial data: a piecewise constant function with a single jump discontinuity separating two states. For a hyperbolic system of $m$ PDEs, each of these states would be a vector with $m$ components, say $q_l$ for $x<0$ and $q_r$ for $x \geq 0$.

We will formally define what we mean by a hyperbolic system of PDEs shortly, but for now it suffices to know that these are the PDEs that naturally arise whenever we model wave propagation. If we study small amplitude waves then the equations may be linear (e.g. acoustics or linear elasticity, modeling seismic waves for example). These PDEs arise from linearizing nonlinear hyperbolic systems that must be used when we consider larger amplitude waves (e.g. the Euler equations of compressible gas dynamics, or nonlinear elastic wave equations).

An example: the dam break problem

To get a feel for what a Riemann problem and it's solution looks like, consider a dam break problem in which a dam initially separates water with constant depth $h_l$ to the left of the dam from water with depth $h_r$ to the right of the dam. Suppose the water is initially stationary, with velocity 0 on both sides of the dam. Now suppose we remove the dam at time $t=0$. What happens for $t>0$?

If the water to the left of the dam is deeper than the water to the right (i.e. $h_l > h_r$) then we expect water to flow from left to right in some manner. In reality we might expect the flow to be somewhat turbulent and multidimensional, but suppose that we could idealize this with a one-dimensional model in which the velocity is purely horizontal and is constant throughout the depth of the fluid, so for $t>0$ the velocity varies only with $x$ and $t$, given by some function $u(x,t)$, and the depth is given by $h(x,t)$.

In the notebook Shallow_water.html we will discuss the one-dimensional shallow water equations that can be used to model this idealized situation.

The animation below shows the solution to the shallow water equations in the situation just described. The water is colored red for water that is initially to the left of the dam and blue for water that is initially to the right. The stripes with different shading are to help you visualize the flow. Think of the color as a dye that is simply carried along with the water as it flows. In the notebook Shallow_water_tracer.ipynb we will see that we can in fact model this by adding another equation to the shallow water equations.

In [1]:
%matplotlib inline
from exact_solvers import shallow_water
from ipywidgets import widgets, fixed
from utils.jsanimate_widgets import interact

Will create JSAnimation figures instead of interactive widget
In [2]:
interact(shallow_water.make_demo_plot_function(h_l=3., h_r=1., u_l=0., u_r=0),
         t=widgets.FloatSlider(min=0.001,max=0.6,step=0.1,value=0.), fig=fixed(0))


Once Loop Reflect

Things to note about this animation

  • Initially the fluid velocity $u(x,0) \equiv 0$ at $t=0$. For $t>0$, the water remain stationary far from away from the dam. The region where the water is not stationary grows linearly outward from $x=0$ as $t$ increases.
  • The fluid is accelerated to the right everywhere as water starts to flow.
  • Water upstream from the dam (for $x<0$) starts to accelerate slowly as the depth of the water starts to fall, and both the depth and the velocity remain continuous functions of $x$ in this region for all time. The depth and velocity vary through a rarefaction wave that moves upstream.
  • The water downstream from the dam is stationary until a shock wave passes by, which instantaneously accelerates the fluid from $u_r=0$ to some constant velocity $u_m >0$. At the same time the depth increases from $h_r$ to $h_m$. This middle state depth and velocity is the same as the depth and velocity that the red water upstream from the dam also reaches once it has passed through the rarefaction wave and been fully accelerated.
  • The shock wave moves at some constant velocity $s_r > u_m$.
  • The rarefaction wave spreads out as time evolves and the left and right edges move at constant velocities $\lambda_{l1}$ and $\lambda_{m1}$. We will see later that these speeds are eigenvalues of a certain Jacobian matrices.

A shock wave is a type of wave often seen when solving nonlinear hyperbolic equations. The solution is discontinuous across a shock wave, which means that special mathematical techniques must be used to make sense of functions like this as solutions of a differential equation. It also means that special algorithmic techniques must be used to compute accurate numerical solutions to such equations. Simply replacing derivatives by finite differences typically does not work well.

In the context of water waves, this shock wave is often called a hydraulic jump. In a real dam break problem the solution would not be discontinuous, but there might be turbulent bore that from far enough away looks like a discontinuity. The nonlinear shallow water equations can be used to model such flows in the sense that the jump in depth and velocity and the propagation speed of the shock wave closely approximate the hydraulic jump and speed of the bore.

Why is the Riemann problem important?

Why study the simple case of piecewise constant initial data? Here are some motivations:

  • As we have just seen, even this simple initial data can lead to complicated solutions.

  • However, the solution is simple enough that there is some hope we can determine it exactly, even for nonlinear systems of PDEs.

  • In the animation above we noted that the shock wave and the two edges of the rarefaction wave are all moving at constant speeds. In fact if we plot contours of the depth (or the velocity) in the $x$-$t$ plane, we find that value of each variable is constant along any ray $x = at$ from the origin. In other words $h(x,t) = H(x/t)$ and $u(x,t) = U(x/t)$ for some functions $H$ and $U$ of a single variable. The solution is said to be a similarity solution and this structure greatly facilitates solving the Riemann problem. This is true even for nonlinear hyperbolic problems provided certain conditions are satisfied.

  • Developing an understanding of how the Riemann solution behaves for all possible choices of left and right states for a hyperbolic problem is essential for understanding the behavior of more complicated solutions.

  • The development of mathematical theory of hyperbolic problems has relied extensively on the understanding that comes from studying Riemann problems.

  • Many robust numerical methods for approximating the solution with more general initial data use the Riemann solution as an essential building block. The most basic method of this type is Godunov's method, described below. Numerical methods typically make use of approximate Riemann solvers, which are often computationally cheaper than an exact solver.

The general Riemann problem for shallow water equations

The more general Riemann problem for the shallow water equations has arbitrary depths $h_l, h_r \geq 0$ and also arbitrary velocities $u_l$ and $u_r$ to the left and right of $x=0$. The dam break problem considered above is a special case where the $u_l = u_r = 0$.

In the dam break problem considered above with $h_l > h_r$, the solution consists of a left-going rarefaction wave and a right-going shock wave. If we instead have $h_l < h_r$ with $u_l=u_r=0$, then the water would flow to the left and the Riemann solution would have a left-going shock wave and a right-going rarefaction wave.

In [3]:
interact(shallow_water.make_demo_plot_function(h_l=1., h_r=3., u_l=0., u_r=0),
         t=widgets.FloatSlider(min=0.001,max=0.6,step=0.1,value=0.), fig=fixed(0));


Once Loop Reflect

By choosing non-zero velocities in the initial data, we can also obtain other combinations of shock and rarefaction waves. For example, if the depth is initially the same everywhere ($h_l = h_r$) while the velocities satisfy $u_l>0$ and $u_r<0$ (modeling two opposing streams of water colliding at $x=0$) then the Riemann solution consists of two shock waves. Note that the water is deeper and stationary between the shock waves:

In [4]:
interact(shallow_water.make_demo_plot_function(h_l=1., h_r=1., u_l=0.8, u_r=-0.8),
         t=widgets.FloatSlider(min=0.001,max=0.6,step=0.1,value=0.), fig=fixed(0));


Once Loop Reflect

On the other hand, if $u_l<0$ and $u_r>0$ (so that the water is moving away from $x=0$ on both sides) then the solution consists of two rarefaction waves:

In [5]:
interact(shallow_water.make_demo_plot_function(h_l=1., h_r=1., u_l=-0.8, u_r=0.8),
         t=widgets.FloatSlider(min=0.001,max=0.6,step=0.1,value=0.), fig=fixed(0));


Once Loop Reflect

In the notebook Shallow_water.html we will show how to solve the Riemann problem for arbitrary left and right states and in particular how to determine whether each wave is a shock or rarefaction wave.

Goals of this book

The collection of Jupyter notebooks that make up this book will help you understand how to solve the Riemann problem for general hyperbolic PDE. This will be accomplished by working through examples for several different systems of equations. In the process you should also gain a better understanding of wave propagation phenomena in the context of several important applications. [Expand this?]

When is a PDE hyperbolic?

In this book we consider only first-order hyperbolic equations, meaning they involve only first derivatives in space and time. Many important hyperbolic PDE models can be written as a system of conservation laws:

\begin{align} \label{intro:multiDconslaw} q_t + \nabla \cdot F(q) & = 0. \end{align}

Here $q(x,t)$ is a vector of conserved quantities (e.g. mass, momentum, energy) and each component of $F(q)$ is the corresponding flux. In practice (and in most of this book) we focus on one-dimensional problems, for which (\ref{intro:multiDconslaw}) can be written

\begin{align} \label{intro:conslaw} q_t + f(q)_x & = 0. \end{align}

If the flux in \eqref{intro:conslaw} is a linear function $f(q)= Aq$ for some matrix $A$, then we have

\begin{align} \label{intro:linconslaw} q_t + A q_x & = 0, \end{align}

We say that \eqref{intro:linconslaw} is hyperbolic if $A$ is diagonalizable with real eigenvalues, and the nonlinear system \eqref{intro:conslaw} is hyperbolic at some state $q$ if the Jacobian matrix $f'(q)$ is diagonalizable with real eigenvalues.

A common property of these equations is that information travels at finite speed and they model waves. As we will see in more detail later, the eigenvalues correspond to wave speeds while the eigenvectors tell us about the relation between the components of $q$ in each wave.

Some examples of physical problems that can be modeled with hyperbolic PDEs (and are included as examples in this book) are:

  • Sound waves
  • Surface water waves
  • Traffic flow
  • Gas dynamics

Godunov's method

  • Add a brief description of Godunov's method and how it uses Riemann solutions.
  • Mention the need for approximate Riemann solvers.

A linear example: acoustic waves

We end this introduction with a sample linear system. Acoustic waves in one dimension can be modeled by the system of PDEs

\begin{align} \label{intro:acoustics} p_t + K u_x & = 0 \\ u_t + \frac{1}{\rho} p_x & = 0, \end{align}

where $p$ represents the deviation from ambient pressure, $u$ is velocity, $K$ is the bulk modulus and $\rho$ is the density of the medium. This first-order system is equivalent to the familiar wave equation $u_{tt} = c^2 u_{xx}$ with $c = \sqrt{K/\rho}$; recall that $c$ represents the speed of propagation of waves in the wave equation. Coming back to the first-order form \eqref{intro:acoustics}, we can write it in the canonical form \eqref{intro:linconslaw} by setting

\begin{align} q & = \begin{bmatrix} p \\ u \end{bmatrix} & A & = \begin{bmatrix} 0 & K \\ 1/\rho & 0 \end{bmatrix}. \end{align}

The eigenvalues of $A$ are $\pm \sqrt{K/\rho} = \pm c$; it is true in general that the speed of propagation of waves in \eqref{intro:linconslaw} is given by the eigenvalues of $A$, and this hints at the importance of the hyperbolicity condition (that $A$ be diagonalizable with real eigenvalues).

The Riemann problem for acoustics

The Riemann problem consists of \eqref{intro:acoustics} together with initial data

\begin{align} (p_0,u_0) & = \begin{cases} (p_l, u_l) & \text{ for } x<0 \\ (p_r, u_r) & \text{ for } x>0. \end{cases} \end{align}

The solution of the Riemann problem for this simple linear system consists of two jumps, one proportional to each of the eigenvectors of $A$, and each moving with velocity equal to the corresponding eigenvalue of $A$. Between these two jumps lies a middle state different from the left and right states; we can think of solving the Riemann problem by finding the intersection (in the $p-u$ plane) of lines passing through the left and right states that are each parallel to the appropriate eigenvector of A. This example is explored in more detail in Acoustics.html.

The figure below shows the phase plane (the $p$-$u$ plane for acoustics) and the location of $q_m$ relative to $q_l$ and $q_r$. You can drag the left and right states around to see how this changes the solution. Also shown in the other two plots are the pressure $p(x,t)$ and velocity $u(x,t)$ at one particular time $t=1$ for this Riemann solution.

In [6]:
from IPython.core.display import display, HTML
app = open('phase_plane_acoustics_small.html','r').read()
display(HTML(app))

Acoustics phase plane

Drag the left and right states in the first plot to change them.

Phase plane plot for the shallow water equations

The phase plane plot for the nonlinear shallow water equations is more interesting. Since the eigenvectors of the Jacobian matrix vary with $q$, the curves connecting states separated by a single wave are not straight lines. There is also a distinction between the curves connecting states by a shock wave and those connecting states by a rarefaction wave. This is discussed in much more detail in the notebook Shallow_water.html, but here is an interactive view that allows you to see how the solution changes as you move the left and right states around in the phase plane:

In [7]:
from IPython.core.display import display, HTML
app = open('phase_plane_shallow_water_verysmall.html','r').read()
display(HTML(app))

Shallow water phase plane

Drag the left and right states in the first plot to change them. Drag the circle in the second plot to adjust the time.
In [8]: