Note

**4/22/2013** A positivity fix has been implemented for the
wave-propagation scheme. The basic idea is to switch to a Lax-flux
based updater when negative density and/or pressure are
detected. Then, the simulation is switched back to using Roe-flux
based updater. All of this is achieved in the Lua script. Results
have been updated accordingly.

Contents

In this entry two algorithms, the wave-propagation algorithm and the MUSCL-Hancock scheme, are benchmarked against a series of 1D shock-tube problems. For most of these problems exact solutions can be computed using an exact Riemann solver, described, for example, in [Kulikovskii2001].

The wave-propagation scheme is implemented in the
`WavePropagationUpdater` in the `slvrs` directory, while the
MUSCL-Hancock scheme is implemented in the `MusclHancock1DUpdater`
in the proto directory. Note that as this point the MUSCL-Hancock
scheme is *not* production quality but is hard-coded to solve just the
1D Euler equations. The code to compute the exact solution is located
in the `sims/code/exactrp` directory.

The problems solved in this entry are quite difficult and really tax the ability of schemes to capture discontinuities accurately. In particular, schemes have a lot of trouble in problems in which vacuum states can form or exist in parts of the domain. Also, ensuring positivity of density and pressure is a difficult problem.

Eight shock-tube problems are solved. Each problem is described below and results of comparing both schemes with the exact solution are shown. Internal energy in the following is computed as \(p/(\gamma-1)\rho\).

Sod-shock with sonic point in rarefaction. Domain is \(x \in [0,1]\), discretized with 100 cells. Gas adiabatic constant of 1.4 was used. Simulation is initialized with a shock at \(x=0.3\), with left and right states

\[\begin{split}\left[
\begin{matrix}
\rho_l \\
u_l \\
p_l
\end{matrix}
\right]
=
\left[
\begin{matrix}
1 \\
0.75 \\
1.0
\end{matrix}
\right],
\qquad
\left[
\begin{matrix}
\rho_r \\
u_r \\
p_r
\end{matrix}
\right]
=
\left[
\begin{matrix}
0.125 \\
0.0 \\
0.1
\end{matrix}
\right].\end{split}\]

and is run to \(t=0.2\).

This problem has a near-vaccum near the location of the discontinuity. Domain is \(x \in [0,1]\), discretized with 100 cells. Gas adiabatic constant of 1.4 is used. Simulation is initialized with a shock at \(x=0.5\), with left and right states

\[\begin{split}\left[
\begin{matrix}
\rho_l \\
u_l \\
p_l
\end{matrix}
\right]
=
\left[
\begin{matrix}
1.0 \\
-2.0 \\
0.4
\end{matrix}
\right],
\qquad
\left[
\begin{matrix}
\rho_r \\
u_r \\
p_r
\end{matrix}
\right]
=
\left[
\begin{matrix}
1.0 \\
2.0 \\
0.4
\end{matrix}
\right].\end{split}\]

and is run to \(t=0.15\).

The second order MUSCL-Hancock **fails** on this problem. The solution
quickly develops negative pressure and density. A positivity fix is
required.

The wave-propagation scheme works with this problem. However, a positivity fix is required. This is implemented by redoing a time-step with Lax fluxes when negative density/pressure is detected and then continuing on with regular Roe fluxes. For this particular problem the negative density/pressure only occurs in the very first time-step and so Roe fluxes can be used for rest of the simulation. The wave-propagation scheme also works with the use of Lax fluxes for the complete simulation. With Lax fluxes used for the complete simulation, the solution is more diffuse, however does not show the strange features around \(x=0.5\).

Results are shown below.

The first-order MUSCL-Hancock also works for this problem. Results are shown below. The wave-propagation scheme seems marginally better than the first-order MUSCL scheme for this problem.

The 1D Noh problem. Domain is \(x \in [0,1]\), discretized with 100 cells. Gas adiabatic constant of \(5/3\) is used. Simulation is initialized with a shock at \(x=0.5\), with left and right states

\[\begin{split}\left[
\begin{matrix}
\rho_l \\
u_l \\
p_l
\end{matrix}
\right]
=
\left[
\begin{matrix}
1.0 \\
1.0 \\
10^{-6}
\end{matrix}
\right],
\qquad
\left[
\begin{matrix}
\rho_r \\
u_r \\
p_r
\end{matrix}
\right]
=
\left[
\begin{matrix}
1.0 \\
-1.0 \\
10^{-6}
\end{matrix}
\right].\end{split}\]

and is run to \(t=1.0\).

The MUSCL-Hancock scheme **fails** on this problem. A positivity fix
needs to be implemented. However, the 1st-order MUSCL-Hancock scheme
works and results are shown below.

1D Euler shock with a stationary contact discontinuity at \(x=0.8\). Domain is \(x \in [0,1]\), discretized with 100 cells. Gas adiabatic constant of \(1.4\) is used. Simulation is initialized with a shock at \(x=0.8\), with left and right states

\[\begin{split}\left[
\begin{matrix}
\rho_l \\
u_l \\
p_l
\end{matrix}
\right]
=
\left[
\begin{matrix}
1.0 \\
-19.59745 \\
1000
\end{matrix}
\right],
\qquad
\left[
\begin{matrix}
\rho_r \\
u_r \\
p_r
\end{matrix}
\right]
=
\left[
\begin{matrix}
1.0 \\
-19.59745 \\
0.01
\end{matrix}
\right].\end{split}\]

and is run to \(t=0.012\).

The MUSCL-Hancock scheme **fails** on this problem. Results with the
1st-order MUSCL-Hancock method is shown below.

1D Euler shock with two strong shocks. Domain is \(x \in [0,1]\), discretized with 100 cells. Gas adiabatic constant of \(1.4\) is used. Simulation is initialized with a shock at \(x=0.4\), with left and right states

\[\begin{split}\left[
\begin{matrix}
\rho_l \\
u_l \\
p_l
\end{matrix}
\right]
=
\left[
\begin{matrix}
5.99924 \\
19.5975 \\
460.894
\end{matrix}
\right],
\qquad
\left[
\begin{matrix}
\rho_r \\
u_r \\
p_r
\end{matrix}
\right]
=
\left[
\begin{matrix}
5.99242 \\
-6.19633 \\
46.0895
\end{matrix}
\right].\end{split}\]

and is run to \(t=0.035\).

1D Euler with a stationary contact discontinuity. Domain is \(x \in [0,1]\), discretized with 100 cells. Gas adiabatic constant of \(1.4\) is used. Simulation is initialized with a shock at \(x=0.5\), with left and right states

\[\begin{split}\left[
\begin{matrix}
\rho_l \\
u_l \\
p_l
\end{matrix}
\right]
=
\left[
\begin{matrix}
1.4 \\
0.0 \\
1.0
\end{matrix}
\right],
\qquad
\left[
\begin{matrix}
\rho_r \\
u_r \\
p_r
\end{matrix}
\right]
=
\left[
\begin{matrix}
1.0 \\
0.0 \\
1.0
\end{matrix}
\right].\end{split}\]

and is run to \(t=2.0\).

1D Euler with a slowly moving contact discontinuity. Domain is \(x \in [0,1]\), discretized with 100 cells. Gas adiabatic constant of \(1.4\) is used. Simulation is initialized with a shock at \(x=0.5\), with left and right states

\[\begin{split}\left[
\begin{matrix}
\rho_l \\
u_l \\
p_l
\end{matrix}
\right]
=
\left[
\begin{matrix}
1.4 \\
0.1 \\
1.0
\end{matrix}
\right],
\qquad
\left[
\begin{matrix}
\rho_r \\
u_r \\
p_r
\end{matrix}
\right]
=
\left[
\begin{matrix}
1.0 \\
0.1 \\
1.0
\end{matrix}
\right].\end{split}\]

and is run to \(t=2.0\).

1D Euler with a sharp peak in density. Domain is \(x \in [0.0,0.5]\), discretized with 800 cells. Gas adiabatic constant of \(1.4\) is used. Simulation is initialized with a shock at \(x=0.4\), with left and right states

\[\begin{split}\left[
\begin{matrix}
\rho_l \\
u_l \\
p_l
\end{matrix}
\right]
=
\left[
\begin{matrix}
0.1261192 \\
8.9047029 \\
782.92899
\end{matrix}
\right],
\qquad
\left[
\begin{matrix}
\rho_r \\
u_r \\
p_r
\end{matrix}
\right]
=
\left[
\begin{matrix}
6.591493 \\
2.2654207 \\
3.1544874
\end{matrix}
\right].\end{split}\]

and is run to \(t=0.0039\).

The MUSCL-Hancock scheme **fails** on this problem. Results with the
1st-order MUSCL-Hancock method is shown below.

The Woodward-Collela blast wave problem consists of two shocks interacting due to reflections off solid walls. The domain is \(x \in [0,1]\), discretized with 400 cells with wall boundary conditions at both ends. Simulation is initialized with two discontinuities, first at \(x_1 = 0.1\) and the other at \(x_2=0.9\). The density and velocity is set everywhere to \(1.0\) and \(0.0\) respectively. The pressure in the three regions, left \(p_l\), middle \(p_m\), and right \(p_r\) are \((p_l,p_m,p_r) = (1000,0.01,100)\). The simulation is run to \(t=0.038\).

In the following, the “exact” solution is computed using wave-propagation method using 2000 cells.

One of the aims of this note was to determine what modifications are needed to the wave-propagation scheme and the MUSCL-Hancock scheme to make them more robust and accurate. Note that the MUSCL-Hancock scheme tested here is only a prototype version and fails on a number of problems. The tests in conducted in this entry will allow a better production quality solver to be developed.

The lessons learned are:

- The wave-propagation scheme needs a positivity fix. For this, a density and pressure floor should be added. More importantly, if the Roe averages lead to a NaN or negative pressure, the Roe fluxes should be replaced (automatically) with a diffusive, but positivity preserving, Rusanov (Lax) flux.
- The MUSCL-Hancock scheme needs a positivity fix also: essentially, if the predicted edge values are negative the slope in the cell should be simply set to zero. This is the main reason why the 2nd order MUSCL-Hancock scheme fails as the predicted edge values do not preserve positivity.
- More accurate (than Rusanov flux) numerical flux needs to be implemented. An HLLC flux will help reduce the diffusion as compared to the wave-propagation scheme.

[Kulikovskii2001] | Andrei G. Kulikoviskii and Nikolai V. Pogorelov
and Andrei Yu. Semenov, Mathematical Aspects of Numerical
Solutions of Hyperbolic Systems, Chapman and Hall/CRC, 2001. |