An approximate solver for Burgers' equation

As a first example, we return to the inviscid Burgers' equation that we studied in Burgers.ipynb:
\begin{align} \label{burgers} q_t + \left(\frac{1}{2}q^2\right)_x & = 0. \end{align}
Although it is easy to solve the Riemann problem for (\ref{burgers}) exactly, it is nevertheless interesting to consider approximate solvers because a numerical scheme does not make use of the full Riemann solution. Furthermore, Burgers' equation provides a simple setting in which some of the intricacies of more complicated approximate solvers can be introduced. Recall that we are interested in approximate solutions that consist entirely of traveling discontinuities.

Shock wave solutions

Recall that the exact Riemann solution for (\ref{burgers}) consists of a single shock or rarefaction wave. We have a shock if $q_l>q_r$ and we have a rarefaction if $q_l < q_r$. In the case of a shock wave, we can simply use the exact solution as our approximation. We have a single wave of strength $q_r-q_l$ traveling at speed $s=(q_r+q_l)/2$.

In terms of fluxes, the numerical flux is $F=f(q_l)$ if $s>0$ and $F=f(q_r)$ if $s<0$. In the special case $s=0$ we have a stationary shock, and it must be that $f(q_l)=f(q_r) (=F)$.

Rarefaction wave solutions

As discussed in Approximate_solvers.ipynb, for numerical purposes it is convenient to approximate a rarefaction wave by a traveling discontinuity. For Burgers' equation this may seem unnecessary, but for more complicated solvers for systems of equations it will be essential.

We will approximate the effect of the rarefaction wave by a fictitious shock:
$${\cal W} = q_r-q_l$$
whose speed is given by the Rankine-Hugoniot jump condition:
$$s = \frac{f(q_r)-f(q_l)}{q_r-q_l} = \frac{q_r + q_l}{2}.$$
Recall that this is indeed a weak solution of the Riemann problem. This fictitious shock is not entropy-satisfying, but as long as it's traveling entirely to the left or entirely to the right, the effect on the numerical solution will be the same as if we used a (entropy-satisfying) rarefaction wave. The numerical flux is again $F=f(q_l)$ if $s>0$ and $F=f(q_r)$ if $s<0$.

Because this is a scalar equation with convex flux, both the Roe and HLL approaches will simplify to what we have already described. But we briefly discuss them here to further illustrate the main ideas.

A Roe solver

Let us consider a linearized solver, in which we replace our nonlinear hyperbolic system with a linearization about some intermediate state $\hat{q}$. For Burgers' equation, the quasilinear form is $q_t + q q_x = 0$ and the linearization gives
$$q_t + \hat{q}q_x = 0.$$
This is simply the advection equation with velocity $\hat{q}$. The solution of the Riemann problem for this equation consists of a wave ${\cal W} = q_r - q_l$ traveling at speed $\hat{q}$. It remains only to determine the state $\hat{q}$, and thus the wave speed.

The defining feature of a Roe linearization is that it gives the exact solution in case the states $(q_r, q_l)$ lie on a single Hugoniot locus -- i.e., when the solution is a single shock. We can achieve this by choosing
$$\hat{q} = \frac{q_r + q_l}{2}.$$
This is equivalent to using the approximate solver described already in the sections above.

Examples

Below we show solutions for three examples; the first involves a shock, the second a (non-transonic) rarefaction, and the third a transonic rarefaction. In the first row we plot the exact solution (computed using code from exact_solvers/burgers.py) in terms of the waves in the x-t plane; in the second row we plot numerical solutions obtained with Clawpack by using a simple first-order method combined with the Riemann solver. The code for the Riemann solver as used in Clawpack can be examined in the Riemann github repository or by looking at the local copy on your computer if you have installed Clawpack.

The solutions obtained for the shock wave and for the first rarefaction wave are good approximations of the true solution. In the case of the transonic rarefaction, however, we see that part of the rarefaction has been replaced by an entropy-violating shock. At the end of this chapter we will show how to apply an entropy fix so that the solver gives a good approximation also in the transonic case.

Two-wave solvers

For Burgers' equation, the Riemann solution consists only of a single wave. It is thus natural to modify the HLL approach by assuming that one of the waves vanishes, and denote the speed of the other wave simply by $s$. Then the conservation condition discussed in Approximate_solvers.ipynb reduces to
$$f(q_r) - f(q_l) = s (q_r - q_l),$$
which is just the Rankine-Hugoniot condition and again leads to the speed discussed above. Since the solution involves only one wave, that wave must carry the entire jump $q_r - q_l$, so this solver is entirely equivalent to that already described.

It is also possible to follow the full HLL approach, taking \begin{align*} s_1 & = \min f'(q) = \min(q_l, q_r) \\ s_2 & = \max f'(q) = \max(q_l, q_r). \end{align*} Regardless of the values of $q_l$ and $q_r$, this leads to
$$q_m = \frac{q_r + q_l}{2},$$
so that each of the two waves carries half of the jump.

Transonic rarefactions

In the approaches above, the solution was approximated by a single wave traveling either to the left or right. For this scalar problem, this "approximation" is, in fact, an exact weak solution of the Riemann problem. As discussed already, we do not typically need to worry about the fact that it may be entropy-violating, since the effect on the numerical solution (after averaging) is identical to that of the entropy-satisfying solution.

However, if $q_l < 0 < q_r$, then the true solution is a transonic rarefaction, in which part of the wave travels to the left and part travels to the right. In this case, the true Riemann solution would lead to changes to both the left and right adjacent cells, whereas an approximate solution with a single wave will only modify one or the other. This leads to an incorrect numerical solution (on the macroscopic level). It is therefore necessary to modify the Riemann solver, imposing what is known as an entropy fix in the transonic case.

Specifically, we use a solution consisting of two waves, each of which captures the net effect of the corresponding rarefaction wave that appears in the exact solution:

\begin{align} {\cal W}_1 & = q_m - q_l, & s_1 = \frac{q_l + q_m}{2} \\ {\cal W}_2 & = q_r - q_m, & s_2 = \frac{q_m + q_r}{2}. \end{align}

For Burgers' equation, the value $q_s=0$ is the sonic point satisfying $f(q_s)=0$. A transonic rarefaction wave takes the value $q_s$ along $x/t = 0$ and so it makes sense to choose $q_m = 0$. The formulas above then simplify to

\begin{align} {\cal W}_1 & = - q_l, & s_1 = \frac{q_l}{2} \\ {\cal W}_2 & = q_r, & s_2 = \frac{q_r}{2}. \end{align}

Note that this can also be viewed as an HLL solver, although not with the usual choice of wave speeds. Choosing the waves speeds $s^1=q_l/2$ and $s^2=q_r/2$ and then solving for $q_m$ by requiring conservation gives $q_m=0$.

We now repeat the example given above, but with the entropy fix applied.

The entropy fix has no effect on the first two solutions, since it is applied only in the case of a transonic rarefaction. The third solution is greatly improved, and will converge to the correct weak solution as the grid is refined.