Riemann Problems

The code solves Riemann Problems for the relativistic Euler equations

\[\begin{split}\partial_t \begin{pmatrix} D \\ S_x \\ S_t \\ \tau \end{pmatrix} + \partial_x \begin{pmatrix} S_x \\ S_x v_x + p \\ S_t v_x \\ (\tau + p) v_x \end{pmatrix} = 0.\end{split}\]

For further details on this system of equations, see the Living Review of Martí and Müller, particularly section 3.1 for the equations and section 8.5 for the solution of the Riemann Problem.

The initial data is piecewise constant: two states \(w_{L, R}\) are specified, each in terms of \(w = (\rho_0, v_x, v_t, \epsilon)\), (the specific rest mass density, normal (\(x\)) and tangential (\(t\)) velocity components, and the specific internal energy). At \(t=0\) the data is set by \(w_L\) for \(x<0\) and \(w_R\) for \(x>0\). Each state has associated with it an equation of state (EOS) to close the set of equations: the EOS does not need to be the same for each state.

Code

To set up a Riemann problem, first set up a left and right state. Each state has its own equation of state. Here we use the first test from the Test Bench section of the Living Review:

from r3d2 import eos_defns, State, RiemannProblem

eos = eos_defns.eos_gamma_law(5.0/3.0)
test_1_U_left = State(10.0, 0.0, 0.0, 2.0, eos, label="L")
test_1_U_right = State(1.0, 0.0, 0.0, 1.5e-6, eos, label="R")
test_1_rp = RiemannProblem(test_1_U_left, test_1_U_right)

The Riemann Problem will produce output directly in the notebook:

test_1_rp

\begin{equation}\begin{cases} \begin{pmatrix} \rho \\ v_x \\ v_t \\ \epsilon \end{pmatrix}_{L} = \begin{pmatrix} 10.0000 \\ 0.0000 \\ 0.0000 \\ 2.0000 \end{pmatrix},\\ \begin{pmatrix} \rho \\ v_x \\ v_t \\ \epsilon \end{pmatrix}_{R} = \begin{pmatrix} 1.0000 \\ 0.0000 \\ 0.0000 \\ 0.0000 \end{pmatrix}, \end{cases} \quad \implies \quad {\cal R}_{\leftarrow}{\cal C}{\cal S}_{\rightarrow}, \quad p_* = 1.4479, \quad\begin{cases} {\cal R}_{\leftarrow}: \lambda^{(0)}\in [-0.7161, 0.1672],\\ {\cal C}: \lambda^{(1)}= 0.7140,\\ {\cal S}_{\rightarrow}: \lambda^{(2)}= 0.8284, \end{cases} \quad \begin{cases} \begin{pmatrix} \rho \\ v_x \\ v_t \\ \epsilon \end{pmatrix}_{\star_L} = \begin{pmatrix} 2.6393 \\ 0.7140 \\ 0.0000 \\ 0.8229 \end{pmatrix},\\ \begin{pmatrix} \rho \\ v_x \\ v_t \\ \epsilon \end{pmatrix}_{\star_R} = \begin{pmatrix} 5.0708 \\ 0.7140 \\ 0.0000 \\ 0.4283 \end{pmatrix}. \end{cases}\end{equation}

This output should be interpreted as follows. First it displays the initial states \({\bf U}_{L, R}\). It then gives the resulting wave pattern: in this case a left going rarefaction, a contact, and a right going shock. It then gives the pressure \(p_*\) in the central states (the pressure and normal velocity do not jump across the contact). Next it gives the speeds of each wave. Finally, it gives the constant states between the left wave and the contact \({\bf U}_{*_L}\) and between the contact and the right state \({\bf U}_{*_R}\).

We can also display the output graphically in the notebook:

from IPython.display import display_png

display_png(test_1_rp)

The first (top left) plot is the wave pattern in the \((x, t)\) plane. Rarefactions are shaded with thin solid lines. Contacts are dashed lines. Discontinuities such as shocks are thick solid lines.

All the other plots show a single quantity as a function of space at a fixed time, or equivalently as a function of the characteristic variable \(\xi = x / t\).

Changing equation of state

There is no need for the two states to have the same equation of state. For example:

eos_air = eos_defns.eos_gamma_law(1.4)
U_vary_eos_L = State(1.0, -0.5, 0.0, 2.0, eos, label="L")
U_vary_eos_R = State(1.0, +0.5, 0.0, 2.0, eos_air, label="R")
test_vary_eos_rp = RiemannProblem(U_vary_eos_L, U_vary_eos_R)
test_vary_eos_rp

\begin{equation}\begin{cases} \begin{pmatrix} \rho \\ v_x \\ v_t \\ \epsilon \end{pmatrix}_{L} = \begin{pmatrix} 1.0000 \\ -0.5000 \\ 0.0000 \\ 2.0000 \end{pmatrix},\\ \begin{pmatrix} \rho \\ v_x \\ v_t \\ \epsilon \end{pmatrix}_{R} = \begin{pmatrix} 1.0000 \\ 0.5000 \\ 0.0000 \\ 2.0000 \end{pmatrix}, \end{cases} \quad \implies \quad {\cal R}_{\leftarrow}{\cal C}{\cal R}_{\rightarrow}, \quad p_* = 0.2598, \quad\begin{cases} {\cal R}_{\leftarrow}: \lambda^{(0)}\in [-0.8955, -0.5734],\\ {\cal C}: \lambda^{(1)}= 0.1224,\\ {\cal R}_{\rightarrow}: \lambda^{(2)}\in [0.6019, 0.8202], \end{cases} \quad \begin{cases} \begin{pmatrix} \rho \\ v_x \\ v_t \\ \epsilon \end{pmatrix}_{\star_L} = \begin{pmatrix} 0.3748 \\ 0.1224 \\ 0.0000 \\ 1.0397 \end{pmatrix},\\ \begin{pmatrix} \rho \\ v_x \\ v_t \\ \epsilon \end{pmatrix}_{\star_R} = \begin{pmatrix} 0.4478 \\ 0.1224 \\ 0.0000 \\ 1.4504 \end{pmatrix}. \end{cases}\end{equation}

display_png(test_vary_eos_rp)

Reactive problems

To include a reaction, make the equation of state of at least one of the states be a reactive EOS:

q_unburnt = 0.1
gamma = 5/3
Cv = 1.0
t_ignition = 2
eos_burnt = eos_defns.eos_gamma_law(gamma)
eos_unburnt = eos_defns.eos_gamma_law_react(gamma, q_unburnt, Cv, t_ignition, eos_burnt)
U_reactive_left = State(1, 0, 0, 2, eos_burnt)
U_reactive_right = State(10, 0, 0, 3, eos_unburnt)
test_reactive_rp = RiemannProblem(U_reactive_left, U_reactive_right)
test_reactive_rp

\begin{equation}\begin{cases} \begin{pmatrix} \rho \\ v_x \\ v_t \\ \epsilon \end{pmatrix}= \begin{pmatrix} 1.0000 \\ 0.0000 \\ 0.0000 \\ 2.0000 \end{pmatrix},\\ \begin{pmatrix} \rho \\ v_x \\ v_t \\ \epsilon \\ q \end{pmatrix}= \begin{pmatrix} 10.0000 \\ 0.0000 \\ 0.0000 \\ 3.0000 \\ 0.1000 \end{pmatrix}, \end{cases} \quad \implies \quad {\cal S}_{\leftarrow}{\cal C}\left({\cal R}_{\rightarrow}{\cal CJDF}_{\rightarrow}\right) , \quad p_* = 5.0558, \quad\begin{cases} {\cal S}_{\leftarrow}: \lambda^{(0)}= -0.8710,\\ {\cal C}: \lambda^{(1)}= -0.5305,\\ \left({\cal R}_{\rightarrow}{\cal CJDF}_{\rightarrow}\right) : \lambda^{(2)}\in [0.2763, 0.6612], \end{cases} \quad \begin{cases} \begin{pmatrix} \rho \\ v_x \\ v_t \\ \epsilon \end{pmatrix}_{\star_L} = \begin{pmatrix} 2.1685 \\ -0.5305 \\ 0.0000 \\ 3.4972 \end{pmatrix},\\ \begin{pmatrix} \rho \\ v_x \\ v_t \\ \epsilon \end{pmatrix}_{\star_R} = \begin{pmatrix} 4.3776 \\ -0.5305 \\ 0.0000 \\ 1.7324 \end{pmatrix}. \end{cases}\end{equation}

display_png(test_reactive_rp)