In [1]:
from IPython.core.display import HTML
def css_styling():
styles = open("./styles/custom.css", "r").read()
return HTML(styles)
css_styling()
Out[1]:
As a simple example of solving the Riemann problem for a nonlinear system we look at the shallow water equations. These are a simplification of the Navier-Stokes equations, reduced here to one spatial dimension $x$, which determine the height $h(x, t)$ of the water with respect to some reference location, and its velocity $u(x, t)$. In the simplest case (where the bed of the channel is flat, and the gravitational constant is renormalised to $1$) these can be written in the conservation law form $$ \partial_t \begin{pmatrix} h \\ h u \end{pmatrix} + \partial_x \begin{pmatrix} hu \\ h u^2 + \tfrac{1}{2} h^2 \end{pmatrix} = {\bf 0}. $$ The conserved variables ${\bf q} = (q_1, q_2)^T = (h, h u)^T$ are effectively the total mass and momentum of the fluid.
As seen in the theory lesson, to construct the solution we need the eigenvalues and eigenvectors of the Jacobian matrix. We can construct them directly, by first noting that, in terms of the conserved variables, $$ {\bf f} = \begin{pmatrix} f_1 \\ f_2 \end{pmatrix} = \begin{pmatrix} q_2 \\ \frac{q_2^2}{q_1} + \frac{q_1^2}{2} \end{pmatrix}. $$ Therefore the Jacobian is $$ \frac{\partial {\bf f}}{\partial {\bf q}} = \begin{pmatrix} 0 & 1 \\ -u^2 + h & 2 u \end{pmatrix}. $$ The eigenvalues and eigenvectors follow immediately as $$ \begin{align} \lambda_{1} & = u - \sqrt{h}, & \lambda_{2} & = u + \sqrt{h}, \\ {\bf r}_{1} & = \begin{pmatrix} 1 \\ u - \sqrt{h} \end{pmatrix}, & {\bf r}_{2} & = \begin{pmatrix} 1 \\ u + \sqrt{h} \end{pmatrix} . \end{align} $$ Here we have followed the standard convention $\lambda_1 \le \lambda_2 \le \dots \le \lambda_N$.
An alternative approach that may be considerably easier to apply for more complex problems is to write down a different quasilinear form of the equation, which in this case is in terms of the primitive variables ${\bf w} = (h, u)^T$, $$ \partial_t \begin{pmatrix} h \\ u \end{pmatrix} + \begin{pmatrix} u & h \\ 1 & u \end{pmatrix} \partial_x \begin{pmatrix} h \\ u \end{pmatrix} = {\bf 0}. $$ The general form here would be written $$ \partial_t {\bf w} + B({\bf w}) \partial_x {\bf w} = {\bf 0}. $$ It is straightforward to check that $$ B = \left( \frac{\partial {\bf q}}{\partial {\bf w}} \right)^{-1} \frac{\partial {\bf f}}{\partial {\bf w}} = \left( \frac{\partial {\bf q}}{\partial {\bf w}} \right)^{-1} \frac{\partial {\bf f}}{\partial {\bf q}} \left( \frac{\partial {\bf q}}{\partial {\bf w}} \right). $$ Thus $B$ is similar to the Jacobian, so must have the same eigenvalues, which is straightforward to check. We also have that $$ \begin{align} B \left\{ \left( \frac{\partial {\bf q}}{\partial {\bf w}} \right)^{-1} {\bf r} \right\} = \lambda \left\{ \left( \frac{\partial {\bf q}}{\partial {\bf w}} \right)^{-1} {\bf r} \right\}, \end{align} $$ showing that the eigenvectors of the Jacobian can be straightforwardly found from the eigenvectors of $B$, which for the shallow water case are $$ \begin{align} {\bf \hat{r}}_1 &= \begin{pmatrix} -\sqrt{h} \\ 1 \end{pmatrix} & {\bf \hat{r}}_2 &= \begin{pmatrix} \sqrt{h} \\ 1 \end{pmatrix}. \end{align} $$
The solution across a continuous rarefaction wave is given by the solution of the ordinary differential equation $$ \partial_{\xi} {\bf q} = \frac{{\bf r}}{{\bf r} \cdot \partial_{{\bf q}} \lambda} $$ where $\lambda, {\bf r}$ are the eigenvalues and eigenvectors of the Jacobian matrix. Note that we can change variables to get the (physically equivalent) relation differential equation $$ \partial_{\xi} {\bf w} = \frac{{\bf r}}{{\bf r} \cdot \partial_{{\bf w}} \lambda} $$ where now the eigenvectors are those of the appropriate matrix for the quasilinear form for ${\bf w}$. Where ${\bf w}$ are the primitive variables as above, the matrix is $B$ and the eigenvectors given by ${\bf \hat{r}}$ as above.
For the shallow water equations we will solve this equation for the primitive variables for the first wave only - symmetry gives the other wave straightforwardly. Starting from $$ \lambda_1 = u - \sqrt{h}, \qquad {\bf \hat{r}}_1 = \begin{pmatrix} -\sqrt{h} \\ 1 \end{pmatrix} $$ we have $$ \partial_{{\bf w}} \lambda_1 = \begin{pmatrix} -\frac{1}{2 \sqrt{h}} \\ 1 \end{pmatrix} $$ and hence $$ {\bf \hat{r}}_1 \cdot \partial_{{\bf w}} \lambda_1 = \frac{3}{2} $$ from which we have $$ \partial_{\xi} \begin{pmatrix} h \\ u \end{pmatrix} = \frac{2}{3} \begin{pmatrix} -\sqrt{h} \\ 1 \end{pmatrix}. $$
This is straightforwardly integrated to get $$ \begin{pmatrix} h \\ u \end{pmatrix} = \begin{pmatrix} \left( c_1 - \frac{\xi}{3} \right)^2 \\ \frac{2}{3} \xi + c_2 \end{pmatrix}. $$
To fix the integration constants $c_{1,2}$ we need to say which state the solution is starting from. As we are looking at the left wave, we expect it to start from the left state ${\bf w}_l = (h_l, u_l)^T$. The left state will connect to the rarefaction wave when the characteristic speeds match, i.e. when $\xi = \xi_l = \lambda_1 = u_l - \sqrt{h_l}$. Therefore we have $$ \begin{pmatrix} h_l \\ u_l \end{pmatrix} = \begin{pmatrix} \left( c_1 - \frac{\xi_l}{3} \right)^2 \\ \frac{2}{3} \xi_l + c_2 \end{pmatrix}, $$ from which we determine $$ c_1 = \frac{1}{3} \xi_l + \sqrt{h_l}, \qquad c_2 = u_l - \frac{2}{3} \xi_l. $$
This gives the final solution $$ \begin{pmatrix} h \\ u \end{pmatrix} = \begin{pmatrix} \left( \frac{\xi_l - \xi}{3} + \sqrt{h_l} \right)^2 \\ \frac{2}{3} (\xi - \xi_l) + u_l \end{pmatrix}. $$
Let us look at all points that can be connected to a certain state by a rarefaction. We do this in the phase plane, which is the $(h, u)$ plane. The "known state" will be given by a marker, and all states along the rarefaction curve given by the line, sometimes known as an integral curve.
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
In [2]:
hl = np.linspace(0.1, 10.1)
ul = np.linspace(-1.0, 1.0)
HL, UL = np.meshgrid(hl, ul)
XIL = UL - np.sqrt(HL)
xi_min = np.min(XIL)
xi_max = np.max(XIL)
h_min = np.min(hl)
h_max = np.max(hl)
u_min = np.min(ul)
u_max = np.max(ul)
xi = np.linspace(xi_min, xi_max)
In [3]:
def plot_sw_rarefaction(hl, ul):
"Plot the rarefaction curve through the state (hl, ul)"
xil = ul - np.sqrt(hl)
h = ((xil - xi) / 3.0 + np.sqrt(hl))**2
u = 2.0 * (xi - xil) / 3.0 + ul
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax.plot(hl, ul, 'rx', markersize = 16, markeredgewidth = 3)
ax.plot(h, u, 'k--', linewidth = 2)
ax.set_xlabel(r"$h$")
ax.set_ylabel(r"$u$")
dh = h_max - h_min
du = u_max - u_min
ax.set_xbound(h_min - 0.1 * dh, h_max + 0.1 * dh)
ax.set_ybound(u_min - 0.1 * du, u_max + 0.1 * du)
fig.tight_layout()
In [4]:
from ipywidgets import interactive, FloatSlider
interactive(plot_sw_rarefaction,
hl = FloatSlider(min = 0.1, max = 10.0, value = 1.0),
ul = FloatSlider(min = -1.0, max = 1.0, value = 0.0))
There is a problem with this: we haven't checked if the states on the curve can be physically connected to this point. That is, we haven't checked how the characteristic speed changes along the curve.
Here it is obvious: we know that the characteristics must spread across the rarefaction, so $\lambda$ must increase, and as $\xi = \lambda$ we must have the characteristic coordinate increasing.
In [5]:
def plot_sw_rarefaction_physical(hl, ul):
"Plot the rarefaction curve through the state (hl, ul)"
xil = ul - np.sqrt(hl)
xi_physical = np.linspace(xil, xi_max)
xi_unphysical = np.linspace(xi_min, xil)
h_physical = ((xil - xi_physical) / 3.0 + np.sqrt(hl))**2
u_physical = 2.0 * (xi_physical - xil) / 3.0 + ul
h_unphysical = ((xil - xi_unphysical) / 3.0 + np.sqrt(hl))**2
u_unphysical = 2.0 * (xi_unphysical - xil) / 3.0 + ul
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax.plot(hl, ul, 'rx', markersize = 16, markeredgewidth = 3)
ax.plot(h_physical, u_physical, 'k-', linewidth = 2, label="Physical")
ax.plot(h_unphysical, u_unphysical, 'k--', linewidth = 2, label="Unphysical")
ax.set_xlabel(r"$h$")
ax.set_ylabel(r"$u$")
dh = h_max - h_min
du = u_max - u_min
ax.set_xbound(h_min - 0.1 * dh, h_max + 0.1 * dh)
ax.set_ybound(u_min - 0.1 * du, u_max + 0.1 * du)
ax.legend()
fig.tight_layout()
In [6]:
interactive(plot_sw_rarefaction_physical,
hl = FloatSlider(min = 0.1, max = 10.0, value = 1.0),
ul = FloatSlider(min = -1.0, max = 1.0, value = 0.0))
We see that along the physical part of the rarefaction curve the height $h$ decreases.
Instead of writing the solution in terms of the similarity coordinate $\xi$ we can instead write the solution in terms of any other single parameter. It is useful to write it in terms of the height, which can be done simply by re-arranging the equations giving $u$ and $h$ in terms of $\xi$. So, a state with height $h_m$ to the right of the state $(h_l, u_l)$ can be connected across a rarefaction if $$ u_m = u_l + 2 \left( \sqrt{h_l} - \sqrt{h_m} \right). $$
In this form we will look at the characteristic curves and the behaviour in state space to cross-check.
In [7]:
def plot_sw_rarefaction_physical_characteristics(hl, ul, hm):
"Plot the rarefaction curve through the state (hl, ul) finishing at (hm, um)"
um = ul + 2.0 * (np.sqrt(hl) - np.sqrt(hm))
h_maximum = np.max([h_max, hl, hm])
h_minimum = np.min([h_min, hl, hm])
u_maximum = np.max([u_max, ul, um])
u_minimum = np.min([u_min, ul, um])
dh = h_maximum - h_minimum
du = u_maximum - u_minimum
xi_min = u_minimum - np.sqrt(h_maximum)
xi_max = u_maximum - np.sqrt(h_minimum)
xil = ul - np.sqrt(hl)
xim = um - np.sqrt(hm)
xi_physical = np.linspace(xil, xi_max)
xi_unphysical = np.linspace(xi_min, xil)
h_physical = ((xil - xi_physical) / 3.0 + np.sqrt(hl))**2
u_physical = 2.0 * (xi_physical - xil) / 3.0 + ul
h_unphysical = ((xil - xi_unphysical) / 3.0 + np.sqrt(hl))**2
u_unphysical = 2.0 * (xi_unphysical - xil) / 3.0 + ul
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(121)
ax1.plot(hl, ul, 'rx', markersize = 16, markeredgewidth = 3, label=r"$(h_l, u_l)$")
ax1.plot(hm, um, 'b+', markersize = 16, markeredgewidth = 3, label=r"$(h_m, u_m)$")
ax1.plot(h_physical, u_physical, 'k-', linewidth = 2, label="Physical")
ax1.plot(h_unphysical, u_unphysical, 'k--', linewidth = 2, label="Unphysical")
ax1.set_xlabel(r"$h$")
ax1.set_ylabel(r"$u$")
ax1.set_xbound(h_minimum - 0.1 * dh, h_maximum + 0.1 * dh)
ax1.set_ybound(u_minimum - 0.1 * du, u_maximum + 0.1 * du)
ax1.legend()
ax2 = fig.add_subplot(122)
left_edge = np.min([-1.0, -1.0 - xil])
right_edge = np.max([1.0, 1.0 - xim])
x_start_points_l = np.linspace(left_edge, 0.0, 20)
x_start_points_r = np.linspace(0.0, right_edge, 20)
x_end_points_l = x_start_points_l + xil
x_end_points_r = x_start_points_r + xim
for xs, xe in zip(x_start_points_l, x_end_points_l):
ax2.plot([xs, xe], [0.0, 1.0], 'b-')
for xs, xe in zip(x_start_points_r, x_end_points_r):
ax2.plot([xs, xe], [0.0, 1.0], 'g-')
# Rarefaction wave
if (xim > xil):
xi = np.linspace(xil, xim, 11)
x_end_rarefaction = xi
for xe in x_end_rarefaction:
ax2.plot([0.0, xe], [0.0, 1.0], 'r--')
else:
x_fill = [x_end_points_l[-1], x_start_points_l[-1], x_end_points_r[0]]
t_fill = [1.0, 0.0, 1.0]
ax2.fill_between(x_fill, t_fill, 1.0, facecolor = 'red', alpha = 0.5)
ax2.set_xbound(-1.0, 1.0)
ax2.set_ybound(0.0, 1.0)
ax2.set_xlabel(r"$x$")
ax2.set_ylabel(r"$t$")
fig.tight_layout()
In [8]:
interactive(plot_sw_rarefaction_physical_characteristics,
hl = FloatSlider(min = 0.1, max = 10.0, value = 1.0),
ul = FloatSlider(min = -1.0, max = 1.0, value = 0.0),
hm = FloatSlider(min = 0.1, max = 10.0, value = 0.5))
We clearly see that only if $h_m < h_l$ do the characteristics spread as they should for a rarefaction. This is, in fact, already given by results above: we showed that $\partial_{\xi} h \propto -\sqrt{h}$. As the height $h$ is positive, this means that as $\xi$ increase across the rarefaction, the height must decrease.
The above exercise assumed we knew the left state and found all right states connecting it by a rarefaction. Now we assume we know both left and right states, and assume they connect to a central state, both along rarefactions.
First, we need to find which states will connect to the right state across a rarefaction.
Repeat the above calculations for states connecting to a known right state. That is, show that, given the right state $(h_r, u_r)$, the left state that connects to it across a rarefaction satisfies $$ \begin{pmatrix} h \\ u \end{pmatrix} = \begin{pmatrix} \left( -\frac{\xi_r - \xi}{3} + \sqrt{h_r} \right)^2 \\ \frac{2}{3} (\xi - \xi_r) + u_r \end{pmatrix}. $$ or equivalently, given $h_m$, that $$ u_m = u_r - 2 \left( \sqrt{h_r} - \sqrt{h_m} \right). $$ Also check that $h$ decreases across the rarefaction, so for a physical solution $h_m < h_r$.
Then we can plot the curve of all states that can be connected to $(h_l, u_l)$ across a left rarefaction, and the curve of all states that can be connected to $(h_r, u_r)$ across a right rarefaction. If they intersect along the physical part of the curve, then we have the solution to the Riemann problem. Clearly this only occurs if $h_m < h_l$ and $h_m < h_r$.
In this case (and note that this is a special case!) we can solve it analytically. We note that, using our assumption that both curves are rarefactions, we have that $$ \begin{align} u_m & = u_l + 2 \left( \sqrt{h_l} - \sqrt{h_m} \right) \\ & = u_r - 2 \left( \sqrt{h_r} - \sqrt{h_m} \right) \end{align} $$ Therefore we have $$ h_m = \frac{1}{16} \left( u_l - u_r + 2 \left( \sqrt{h_l} + \sqrt{h_r} \right) \right)^2. $$
In [9]:
def plot_sw_all_rarefaction(hl, ul, hr, ur):
"Plot the all rarefaction solution curve for states (hl, ul) and (hr, ur)"
hm = (ul - ur + 2.0 * (np.sqrt(hl) + np.sqrt(hr)))**2 / 16.0
um = ul + 2.0 * (np.sqrt(hl) - np.sqrt(hm))
h_maximum = np.max([h_max, hl, hr, hm])
h_minimum = np.min([h_min, hl, hr, hm])
u_maximum = np.max([u_max, ul, ur, um])
u_minimum = np.min([u_min, ul, ur, um])
dh = h_maximum - h_minimum
du = u_maximum - u_minimum
xil_min = u_minimum - np.sqrt(h_maximum)
xil_max = u_maximum - np.sqrt(h_minimum)
xir_min = u_minimum + np.sqrt(h_minimum)
xir_max = u_maximum + np.sqrt(h_maximum)
xil = ul - np.sqrt(hl)
xilm = um - np.sqrt(hm)
xil_physical = np.linspace(xil, xil_max)
xil_unphysical = np.linspace(xil_min, xil)
hl_physical = ((xil - xil_physical) / 3.0 + np.sqrt(hl))**2
ul_physical = 2.0 * (xil_physical - xil) / 3.0 + ul
hl_unphysical = ((xil - xil_unphysical) / 3.0 + np.sqrt(hl))**2
ul_unphysical = 2.0 * (xil_unphysical - xil) / 3.0 + ul
xir = ur + np.sqrt(hr)
xirm = um + np.sqrt(hm)
xir_unphysical = np.linspace(xir, xir_max)
xir_physical = np.linspace(xir_min, xir)
hr_physical = (-(xir - xir_physical) / 3.0 + np.sqrt(hr))**2
ur_physical = 2.0 * (xir_physical - xir) / 3.0 + ur
hr_unphysical = (-(xir - xir_unphysical) / 3.0 + np.sqrt(hr))**2
ur_unphysical = 2.0 * (xir_unphysical - xir) / 3.0 + ur
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(111)
if (hm < np.min([hl, hr])):
ax1.plot(hm, um, 'b+', markersize = 16, markeredgewidth = 3,
label=r"$(h_m, u_m)$, physical solution")
else:
ax1.plot(hm, um, 'b+', markersize = 16, markeredgewidth = 3,
label=r"$(h_m, u_m)$, not physical solution")
ax1.plot(hl, ul, 'rx', markersize = 16, markeredgewidth = 3, label=r"$(h_l, u_l)$")
ax1.plot(hr, ur, 'go', markersize = 16, markeredgewidth = 3, label=r"$(h_r, u_r)$")
ax1.plot(hl_physical, ul_physical, 'k-', linewidth = 2, label="Physical (left)")
ax1.plot(hl_unphysical, ul_unphysical, 'k--', linewidth = 2, label="Unphysical (left)")
ax1.plot(hr_physical, ur_physical, 'c-', linewidth = 2, label="Physical (right)")
ax1.plot(hr_unphysical, ur_unphysical, 'c--', linewidth = 2, label="Unphysical (right)")
ax1.set_xlabel(r"$h$")
ax1.set_ylabel(r"$u$")
ax1.set_xbound(h_minimum - 0.1 * dh, h_maximum + 0.1 * dh)
ax1.set_ybound(u_minimum - 0.1 * du, u_maximum + 0.1 * du)
ax1.legend()
fig.tight_layout()
In [10]:
interactive(plot_sw_all_rarefaction,
hl = FloatSlider(min = 0.1, max = 10.0, value = 1.0),
ul = FloatSlider(min = -1.0, max = 1.0, value = -0.5),
hr = FloatSlider(min = 0.1, max = 10.0, value = 1.0),
ur = FloatSlider(min = -1.0, max = 1.0, value = 0.5))
Given the central state and the relation along rarefaction curves, we can then construct the characteristics and the solution in terms of the similarity coordinate (which, given a time $t$, gives the solution as a function of $x$).
In [11]:
def plot_sw_all_rarefaction_solution(hl, ul, hr, ur):
"Plot the all rarefaction solution curve for states (hl, ul) and (hr, ur)"
hm = (ul - ur + 2.0 * (np.sqrt(hl) + np.sqrt(hr)))**2 / 16.0
um = ul + 2.0 * (np.sqrt(hl) - np.sqrt(hm))
xi1l = ul - np.sqrt(hl)
xi1m = um - np.sqrt(hm)
xi1r = ur - np.sqrt(hr)
hl_raref = np.linspace(hl, hm, 20)
ul_raref = ul + 2.0 * (np.sqrt(hl) - np.sqrt(hl_raref))
xil_raref = ul_raref - np.sqrt(hl_raref)
xi2r = ur + np.sqrt(hr)
xi2m = um + np.sqrt(hm)
xi2l = ul + np.sqrt(hl)
hr_raref = np.linspace(hm, hr)
ur_raref = ur - 2.0 * (np.sqrt(hr) - np.sqrt(hr_raref))
xir_raref = ur_raref + np.sqrt(hr_raref)
xi_min = np.min([-1.0, xi1l, xi1m, xi2r, xi2m])
xi_max = np.max([1.0, xi1l, xi1m, xi2r, xi2m])
d_xi = xi_max - xi_min
h_max = np.max([hl, hr, hm])
h_min = np.min([hl, hr, hm])
d_h = h_max - h_min
u_max = np.max([ul, ur, um])
u_min = np.min([ul, ur, um])
d_u = u_max - u_min
xi = np.array([xi_min - 0.1 * d_xi, xi1l])
h = np.array([hl, hl])
u = np.array([ul, ul])
xi = np.append(xi, xil_raref)
h = np.append(h, hl_raref)
u = np.append(u, ul_raref)
xi = np.append(xi, [xi1m, xi2m])
h = np.append(h, [hm, hm])
u = np.append(u, [um, um])
xi = np.append(xi, xir_raref)
h = np.append(h, hr_raref)
u = np.append(u, ur_raref)
xi = np.append(xi, [xi2r, xi_max + 0.1 * d_xi])
h = np.append(h, [hr, hr])
u = np.append(u, [ur, ur])
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(221)
if (hm < np.min([hl, hr])):
ax1.plot(xi, h, 'b-', label = "Physical solution")
else:
ax1.plot(xi, h, 'r--', label = "Unphysical solution")
ax1.set_ybound(h_min - 0.1 * d_h, h_max + 0.1 * d_h)
ax1.set_xlabel(r"$\xi$")
ax1.set_ylabel(r"$h$")
ax1.legend()
ax2 = fig.add_subplot(222)
if (hm < np.min([hl, hr])):
ax2.plot(xi, u, 'b-', label = "Physical solution")
else:
ax2.plot(xi, u, 'r--', label = "Unphysical solution")
ax2.set_ybound(u_min - 0.1 * d_u, u_max + 0.1 * d_u)
ax2.set_xlabel(r"$\xi$")
ax2.set_ylabel(r"$u$")
ax2.legend()
ax3 = fig.add_subplot(223)
left_end = np.min([-1.0, 1.1*xi1l])
right_end = np.max([1.0, 1.1*xi2r])
left_edge = left_end - xi1l
right_edge = right_end - xi1r
x1_start_points_l = np.linspace(np.min([left_edge, left_end]), 0.0, 20)
x1_start_points_r = np.linspace(0.0, np.max([right_edge, right_end]), 20)
x1_end_points_l = x1_start_points_l + xi1l
t1_end_points_r = np.ones_like(x1_start_points_r)
# Look for intersections
t1_end_points_r = np.minimum(t1_end_points_r, x1_start_points_r / (xi2r - xi1r))
x1_end_points_r = x1_start_points_r + xi1r * t1_end_points_r
# Note: here we are cheating, and using the characteristic speed of the middle state,
# ignoring howo it varies across the rarefaction
x1_final_points_r = x1_end_points_r + (1.0 - t1_end_points_r) * xi1m
for xs, xe in zip(x1_start_points_l, x1_end_points_l):
ax3.plot([xs, xe], [0.0, 1.0], 'b-')
for xs, xe, te in zip(x1_start_points_r, x1_end_points_r, t1_end_points_r):
ax3.plot([xs, xe], [0.0, te], 'g-')
for xs, xe, ts in zip(x1_end_points_r, x1_final_points_r, t1_end_points_r):
ax3.plot([xs, xe], [ts, 1.0], 'g-')
# Highlight the edges of both rarefactions
ax3.plot([0.0, xi1l], [0.0, 1.0], 'r-', linewidth=2)
ax3.plot([0.0, xi1m], [0.0, 1.0], 'r-', linewidth=2)
ax3.plot([0.0, xi2m], [0.0, 1.0], 'r-', linewidth=2)
ax3.plot([0.0, xi2r], [0.0, 1.0], 'r-', linewidth=2)
# Rarefaction wave
if (xi1l < xi1m):
xi = np.linspace(xi1l, xi1m, 11)
x_end_rarefaction = xi
for xe in x_end_rarefaction:
ax3.plot([0.0, xe], [0.0, 1.0], 'r--')
else:
x_fill = [xi1l, 0.0, xi1m]
t_fill = [1.0, 0.0, 1.0]
ax3.fill_between(x_fill, t_fill, 1.0, facecolor = 'red', alpha = 0.5)
ax3.set_xlabel(r"$x$")
ax3.set_ylabel(r"$t$")
ax3.set_title("1-characteristics")
ax3.set_xbound(left_end, right_end)
ax4 = fig.add_subplot(224)
left_end = np.min([-1.0, 1.1*xi1l])
right_end = np.max([1.0, 1.1*xi2r])
left_edge = left_end - xi2l
right_edge = right_end - xi2r
x2_start_points_l = np.linspace(np.min([left_edge, left_end]), 0.0, 20)
x2_start_points_r = np.linspace(0.0, np.max([right_edge, right_end]), 20)
x2_end_points_r = x2_start_points_r + xi2r
t2_end_points_l = np.ones_like(x2_start_points_l)
# Look for intersections
t2_end_points_l = np.minimum(t2_end_points_l, x2_start_points_l / (xi1l - xi2r))
x2_end_points_l = x2_start_points_l + xi2r * t2_end_points_l
# Note: here we are cheating, and using the characteristic speed of the middle state,
# ignoring howo it varies across the rarefaction
x2_final_points_l = x2_end_points_l + (1.0 - t2_end_points_l) * xi2m
for xs, xe in zip(x2_start_points_r, x2_end_points_r):
ax4.plot([xs, xe], [0.0, 1.0], 'g-')
for xs, xe, te in zip(x2_start_points_l, x2_end_points_l, t2_end_points_l):
ax4.plot([xs, xe], [0.0, te], 'b-')
for xs, xe, ts in zip(x2_end_points_l, x2_final_points_l, t2_end_points_l):
ax4.plot([xs, xe], [ts, 1.0], 'b-')
# Highlight the edges of both rarefactions
ax4.plot([0.0, xi1l], [0.0, 1.0], 'r-', linewidth=2)
ax4.plot([0.0, xi1m], [0.0, 1.0], 'r-', linewidth=2)
ax4.plot([0.0, xi2m], [0.0, 1.0], 'r-', linewidth=2)
ax4.plot([0.0, xi2r], [0.0, 1.0], 'r-', linewidth=2)
# Rarefaction wave
if (xi2r > xi2m):
xi = np.linspace(xi2m, xi2r, 11)
x_end_rarefaction = xi
for xe in x_end_rarefaction:
ax4.plot([0.0, xe], [0.0, 1.0], 'r--')
else:
x_fill = [xi2m, 0.0, xi2r]
t_fill = [1.0, 0.0, 1.0]
ax4.fill_between(x_fill, t_fill, 1.0, facecolor = 'red', alpha = 0.5)
ax4.set_xlabel(r"$x$")
ax4.set_ylabel(r"$t$")
ax4.set_title("2-characteristics")
ax4.set_xbound(left_end, right_end)
fig.tight_layout()
In [12]:
interactive(plot_sw_all_rarefaction_solution,
hl = FloatSlider(min = 0.1, max = 10.0, value = 1.0),
ul = FloatSlider(min = -1.0, max = 1.0, value = -0.5),
hr = FloatSlider(min = 0.1, max = 10.0, value = 1.0),
ur = FloatSlider(min = -1.0, max = 1.0, value = 0.5))
We note that the general theory tells us that across a shock the Rankine-Hugoniot conditions $$ V_s \left[ {\bf q} \right] = \left[ {\bf f}({\bf q}) \right] $$ must be satisfied.
For the shallow water equations we will start, as with the rarefaction case, by assuming we know the left state ${\bf q}_l = (h_l, u_l)$, and work out which states ${\bf q}_m$ can be connected to it across a shock.
Note here that the procedure is identical for the right state as the direction does not matter. However, there will be multiple solutions, and checking which is physically correct does require checking whether the left or the right state is known
Writing out the conditions in full we see that $$ \begin{align} V_s \left( h_m - h_l \right) & = h_m u_m - h_l u_l \\ V_s \left( h_m u_m - h_l u_l \right) & = h_m u_m^2 + \tfrac{1}{2} h_m^2 - h_l u_l^2 - \tfrac{1}{2} h_l^2 \end{align} $$
Eliminating the shock speed $V_s$ gives, using the second equation, $$ u_m^2 - (2 u_l) u_m + \left[ u_l^2 - \tfrac{1}{2} \left( h_l - h_m \right) \left( \frac{h_l}{h_m} - \frac{h_m}{h_l} \right) \right] = 0. $$ This has the solutions (assuming that $h_m$ is known!) $$ u_m = u_l \pm \sqrt{\tfrac{1}{2} \left( h_l - h_m \right) \left( \frac{h_l}{h_m} - \frac{h_m}{h_l} \right)}. $$
We can again use the Rankine-Hugoniot relations to find the shock speed. $$ V_s = u_l \pm \frac{h_m}{h_m - h_l} \sqrt{\tfrac{1}{2} \left( h_l - h_m \right) \left( \frac{h_l}{h_m} - \frac{h_m}{h_l} \right)}. $$
We should at this point find which sign is appropriate. Comparing the shock speeds against the characteristic speed will show that
However, we can see this by plotting the Hugoniot locus: the curve of all states that can be connected to $(h_l, u_l)$ across a shock.
In [13]:
def plot_sw_shock_physical(hl, ul):
"Plot the shock curve through the state (hl, ul)"
h = np.linspace(h_min, h_max, 500)
u_negative = ul - np.sqrt(0.5 * (hl - h) * (hl / h - h / hl))
u_positive = ul + np.sqrt(0.5 * (hl - h) * (hl / h - h / hl))
vs_negative = ul - h / (h - hl) * np.sqrt(0.5 * (hl - h) * (hl / h - h / hl))
vs_positive = ul + h / (h - hl) * np.sqrt(0.5 * (hl - h) * (hl / h - h / hl))
xi1_negative = u_negative - np.sqrt(h)
xi1_positive = u_positive - np.sqrt(h)
xi2_negative = u_negative + np.sqrt(h)
xi2_positive = u_positive + np.sqrt(h)
xi1_l = ul - np.sqrt(hl)
xi2_l = ul + np.sqrt(hl)
h1_physical = h[np.logical_and(xi1_negative <= vs_negative, xi1_l >= vs_negative)]
u1_physical = u_negative[np.logical_and(xi1_negative <= vs_negative, xi1_l >= vs_negative)]
h2_physical = h[np.logical_and(xi2_positive >= vs_positive, xi2_l <= vs_positive)]
u2_physical = u_positive[np.logical_and(xi2_positive >= vs_positive, xi2_l <= vs_positive)]
h1_unphysical = h[np.logical_or(xi1_negative >= vs_negative, xi1_l <= vs_negative)]
u1_unphysical = u_negative[np.logical_or(xi1_negative >= vs_negative, xi1_l <= vs_negative)]
h2_unphysical = h[np.logical_or(xi2_positive <= vs_positive, xi2_l >= vs_positive)]
u2_unphysical = u_positive[np.logical_or(xi2_positive <= vs_positive, xi2_l >= vs_positive)]
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax.plot(hl, ul, 'rx', markersize = 16, markeredgewidth = 3)
ax.plot(h1_physical, u1_physical, 'b-', linewidth = 2,
label="Physical, 1-shock")
ax.plot(h1_unphysical, u1_unphysical, 'b--', linewidth = 2,
label="Unphysical, 1-shock")
ax.plot(h2_physical, u2_physical, 'g-', linewidth = 2,
label="Physical, 2-shock")
ax.plot(h2_unphysical, u2_unphysical, 'g--', linewidth = 2,
label="Unphysical, 2-shock")
ax.plot(h[::5], u_negative[::5], 'co', markersize = 12, markeredgewidth = 2, alpha = 0.3,
label="Negative branch")
ax.plot(h[::5], u_positive[::5], 'ro', markersize = 12, markeredgewidth = 2, alpha = 0.3,
label="Positive branch")
ax.set_xlabel(r"$h$")
ax.set_ylabel(r"$u$")
dh = h_max - h_min
du = u_max - u_min
ax.set_xbound(h_min, h_max)
ax.set_ybound(u_min, u_max)
ax.legend()
fig.tight_layout()
In [14]:
interactive(plot_sw_shock_physical,
hl = FloatSlider(min = 0.1, max = 10.0, value = 1.0),
ul = FloatSlider(min = -1.0, max = 1.0, value = 0.0))
We see from these results, as claimed above, that
When we assumed the solution contained two rarefactions it was possible to write the full solution in closed form. If we assume the solution contains two shocks then it is not possible to do this. However, it is straightforward to find the solution numerically.
We assume the left state ${\bf w}_l = (h_l, u_l)$ and the right state ${\bf w}_r = (h_r, u_r)$ are known, and that they both connect to the central state ${\bf w}_m = (h_m, u_m)$ through shocks. We know that $$ \begin{align} u_m & = u_l - \sqrt{\tfrac{1}{2} \left( h_l - h_m \right) \left( \frac{h_l}{h_m} - \frac{h_m}{h_l} \right)}, \\ u_m & = u_r + \sqrt{\tfrac{1}{2} \left( h_r - h_m \right) \left( \frac{h_r}{h_m} - \frac{h_m}{h_r} \right)}. \end{align} $$ We schematically write these equations as $$ \begin{align} u_m & = \phi_l \left( h_m; {\bf w}_l \right), \\ u_m & = \phi_r \left( h_m; {\bf w}_r \right), \end{align} $$ to indicate that the velocity in the central state, $u_m$, can be written as a function of the single unknown $h_m$ and known data.
We immediately see that $h_m$ is a root of the nonlinear equation $$ \phi \left( h_m; {\bf w}_l, {\bf w}_r \right) = \phi_l \left( h_m; {\bf w}_l \right) - \phi_r \left( h_m; {\bf w}_r \right) = 0. $$
Finding the roots of scalar nonlinear equations is a standard problem in numerical methods, with methods such as bisection, Newton-Raphson and more being well-known. scipy provides a number of standard algorithms - here we will use the recommended brentq method.
Note that as soon as we have numerically determined $h_m$ then either formula above gives $u_m$, and the shock speeds follow.
In [15]:
def plot_sw_all_shock(hl, ul, hr, ur):
"Plot the all shock solution curve for states (hl, ul) and (hr, ur)"
from scipy.optimize import brentq
def phi(hstar):
"Function defining the root"
phi_l = ul - np.sqrt(0.5 * (hl - hstar) * (hl / hstar - hstar / hl))
phi_r = ur + np.sqrt(0.5 * (hr - hstar) * (hr / hstar - hstar / hr))
return phi_l - phi_r
# There is a solution only in the physical case.
physical_solution = True
try:
hm = brentq(phi, np.max([hl, hr]), 10.0 * h_max)
except ValueError:
physical_solution = False
hm = hl
um = ul - np.sqrt(0.5 * (hl - hm) * (hl / hm - hm / hl))
h = np.linspace(h_min, h_max, 500)
u_negative = ul - np.sqrt(0.5 * (hl - h) * (hl / h - h / hl))
u_positive = ur + np.sqrt(0.5 * (hr - h) * (hr / h - h / hr))
h_maximum = np.max([h_max, hl, hr, hm])
h_minimum = np.min([h_min, hl, hr, hm])
u_maximum = np.max([u_max, ul, ur, um])
u_minimum = np.min([u_min, ul, ur, um])
dh = h_maximum - h_minimum
du = u_maximum - u_minimum
xil_min = u_minimum - np.sqrt(h_maximum)
xil_max = u_maximum - np.sqrt(h_minimum)
xir_min = u_minimum + np.sqrt(h_minimum)
xir_max = u_maximum + np.sqrt(h_maximum)
vs_negative = ul - h / (h - hl) * np.sqrt(0.5 * (hl - h) * (hl / h - h / hl))
vs_positive = ur + h / (h - hr) * np.sqrt(0.5 * (hr - h) * (hr / h - h / hr))
xi1_negative = u_negative - np.sqrt(h)
xi1_positive = u_positive - np.sqrt(h)
xi2_negative = u_negative + np.sqrt(h)
xi2_positive = u_positive + np.sqrt(h)
xi1_l = ul - np.sqrt(hl)
xi2_r = ur + np.sqrt(hr)
h1_physical = h[np.logical_and(xi1_negative <= vs_negative, xi1_l >= vs_negative)]
u1_physical = u_negative[np.logical_and(xi1_negative <= vs_negative, xi1_l >= vs_negative)]
h2_physical = h[np.logical_and(xi2_positive >= vs_positive, xi2_r <= vs_positive)]
u2_physical = u_positive[np.logical_and(xi2_positive >= vs_positive, xi2_r <= vs_positive)]
h1_unphysical = h[np.logical_or(xi1_negative >= vs_negative, xi1_l <= vs_negative)]
u1_unphysical = u_negative[np.logical_or(xi1_negative >= vs_negative, xi1_l <= vs_negative)]
h2_unphysical = h[np.logical_or(xi2_positive <= vs_positive, xi2_r >= vs_positive)]
u2_unphysical = u_positive[np.logical_or(xi2_positive <= vs_positive, xi2_r >= vs_positive)]
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax.plot(hl, ul, 'rx', markersize = 16, markeredgewidth = 3, label = r"${\bf w}_l$")
ax.plot(hr, ur, 'r+', markersize = 16, markeredgewidth = 3, label = r"${\bf w}_r$")
if physical_solution:
ax.plot(hm, um, 'ro', markersize = 16, markeredgewidth = 3, label = r"${\bf w}_m$")
ax.plot(h1_physical, u1_physical, 'b-', linewidth = 2,
label="Physical, 1-shock")
ax.plot(h1_unphysical, u1_unphysical, 'b--', linewidth = 2,
label="Unphysical, 1-shock")
ax.plot(h2_physical, u2_physical, 'g-', linewidth = 2,
label="Physical, 2-shock")
ax.plot(h2_unphysical, u2_unphysical, 'g--', linewidth = 2,
label="Unphysical, 2-shock")
ax.set_xlabel(r"$h$")
ax.set_ylabel(r"$u$")
ax.set_xbound(h_minimum - 0.1 * dh, h_maximum + 0.1 * dh)
ax.set_ybound(u_minimum - 0.1 * du, u_maximum + 0.1 * du)
ax.legend()
fig.tight_layout()
In [16]:
interactive(plot_sw_all_shock,
hl = FloatSlider(min = 0.1, max = 10.0, value = 1.0),
ul = FloatSlider(min = -1.0, max = 1.0, value = 0.2),
hr = FloatSlider(min = 0.1, max = 10.0, value = 1.0),
ur = FloatSlider(min = -1.0, max = 1.0, value = -0.2))
Finally, we can plot the solution in physical space.
In [17]:
def plot_sw_all_shock_solution(hl, ul, hr, ur):
"Plot the all shock solution for states (hl, ul) and (hr, ur)"
from scipy.optimize import brentq
def phi(hstar):
"Function defining the root"
phi_l = ul - np.sqrt(0.5 * (hl - hstar) * (hl / hstar - hstar / hl))
phi_r = ur + np.sqrt(0.5 * (hr - hstar) * (hr / hstar - hstar / hr))
return phi_l - phi_r
# There is a solution only in the physical case.
physical_solution = True
try:
hm = brentq(phi, np.max([hl, hr]), 10.0 * h_max)
except ValueError:
physical_solution = False
hm = hl
um = ul - np.sqrt(0.5 * (hl - hm) * (hl / hm - hm / hl))
xi1l = ul - np.sqrt(hl)
xi1m = um - np.sqrt(hm)
xi1r = ur - np.sqrt(hr)
if physical_solution:
vsl = ul - hm / (hm - hl) * np.sqrt(0.5 * (hl - hm) * (hl / hm - hm / hl))
else:
vsl = xi1l
xi2r = ur + np.sqrt(hr)
xi2m = um + np.sqrt(hm)
xi2l = ul + np.sqrt(hl)
if physical_solution:
vsr = ur + hm / (hm - hr) * np.sqrt(0.5 * (hr - hm) * (hr / hm - hm / hr))
else:
vsr = xi2r
xi_min = np.min([-1.0, xi1l, xi1m, xi2r, xi2m])
xi_max = np.max([1.0, xi1l, xi1m, xi2r, xi2m])
d_xi = xi_max - xi_min
h_maximum = np.max([hl, hr, hm])
h_minimum = np.min([hl, hr, hm])
d_h = h_maximum - h_minimum
u_maximum = np.max([ul, ur, um])
u_minimum = np.min([ul, ur, um])
d_u = u_maximum - u_minimum
xi = np.array([xi_min - 0.1 * d_xi, vsl, vsl, vsr, vsr, xi_max + 0.1 * d_xi])
h = np.array([hl, hl, hm, hm, hr, hr])
u = np.array([ul, ul, um, um, ur, ur])
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(221)
if (hm > np.max([hl, hr])):
ax1.plot(xi, h, 'b-', label = "Physical solution")
else:
ax1.plot(xi, h, 'r--', label = "Unphysical solution")
ax1.set_ybound(h_minimum - 0.1 * d_h, h_maximum + 0.1 * d_h)
ax1.set_xlabel(r"$\xi$")
ax1.set_ylabel(r"$h$")
ax1.legend()
ax2 = fig.add_subplot(222)
if (hm > np.max([hl, hr])):
ax2.plot(xi, u, 'b-', label = "Physical solution")
else:
ax2.plot(xi, u, 'r--', label = "Unphysical solution")
ax2.set_ybound(u_minimum - 0.1 * d_u, u_maximum + 0.1 * d_u)
ax2.set_xlabel(r"$\xi$")
ax2.set_ylabel(r"$u$")
ax2.legend()
ax3 = fig.add_subplot(223)
left_end = np.min([-1.0, 1.1*xi1l])
right_end = np.max([1.0, 1.1*xi2r])
left_edge = left_end - xi1l
right_edge = right_end - xi1r
x1_start_points_l = np.linspace(np.min([left_edge, left_end]), 0.0, 20)
x1_start_points_r = np.linspace(0.0, np.max([right_edge, right_end]), 20)
t1_end_points_l = np.ones_like(x1_start_points_l)
t1_end_points_r = np.ones_like(x1_start_points_r)
# Look for intersections
t1_end_points_l = np.minimum(t1_end_points_l, x1_start_points_l / (vsl - xi1l))
x1_end_points_l = x1_start_points_l + xi1l * t1_end_points_l
t1_end_points_r = np.minimum(t1_end_points_r, x1_start_points_r / (vsr - xi1r))
x1_end_points_r = x1_start_points_r + xi1r * t1_end_points_r
# Note: here we are cheating, and using the characteristic speed of the middle state,
# ignoring how it varies across the rarefaction
t1_final_points_r = np.ones_like(x1_start_points_r)
t1_final_points_r = np.minimum(t1_final_points_r,
(x1_end_points_r - t1_end_points_r * xi1m) / (vsl - xi1m))
x1_final_points_r = x1_end_points_r + (t1_final_points_r - t1_end_points_r) * xi1m
for xs, xe, te in zip(x1_start_points_l, x1_end_points_l, t1_end_points_l):
ax3.plot([xs, xe], [0.0, te], 'b-')
for xs, xe, te in zip(x1_start_points_r, x1_end_points_r, t1_end_points_r):
ax3.plot([xs, xe], [0.0, te], 'g-')
for xs, xe, ts, te in zip(x1_end_points_r, x1_final_points_r, t1_end_points_r,
t1_final_points_r):
ax3.plot([xs, xe], [ts, te], 'g-')
# Highlight the shocks
ax3.plot([0.0, vsl], [0.0, 1.0], 'r-', linewidth=2)
ax3.plot([0.0, vsr], [0.0, 1.0], 'r-', linewidth=2)
# Unphysical shock
if not physical_solution:
x_fill = []
if xi1l < xi1m:
x_fill = [xi1l, 0.0, xi1m]
elif xi1l < vsl:
x_fill = [xi1l, 0.0, vsl]
elif vsl < xi1m:
x_fill = [vsl, 0.0, xi1m]
if len(x_fill) > 0:
t_fill = [1.0, 0.0, 1.0]
ax3.fill_between(x_fill, t_fill, 1.0, facecolor = 'red', alpha = 0.5)
x_fill = []
if xi2r > xi2m:
x_fill = [xi2m, 0.0, xi2r]
elif xi2m < vsr:
x_fill = [xi2m, 0.0, vsr]
elif vsr < xi2r:
x_fill = [vsr, 0.0, xi2r]
if len(x_fill) > 0:
t_fill = [1.0, 0.0, 1.0]
ax3.fill_between(x_fill, t_fill, 1.0, facecolor = 'red', alpha = 0.5)
ax3.set_xlabel(r"$x$")
ax3.set_ylabel(r"$t$")
ax3.set_title("1-characteristics")
ax3.set_xbound(left_end, right_end)
ax4 = fig.add_subplot(224)
left_end = np.min([-1.0, 1.1*xi1l])
right_end = np.max([1.0, 1.1*xi2r])
left_edge = left_end - xi2l
right_edge = right_end - xi2r
x2_start_points_l = np.linspace(np.min([left_edge, left_end]), 0.0, 20)
x2_start_points_r = np.linspace(0.0, np.max([right_edge, right_end]), 20)
x2_end_points_r = x2_start_points_r + xi2r
t2_end_points_l = np.ones_like(x2_start_points_l)
t2_end_points_r = np.ones_like(x2_start_points_r)
# Look for intersections
t2_end_points_r = np.minimum(t2_end_points_r, x2_start_points_r / (vsr - xi2r))
x2_end_points_r = x2_start_points_r + xi2r * t2_end_points_r
t2_end_points_l = np.minimum(t2_end_points_l, x2_start_points_l / (vsl - xi2l))
x2_end_points_l = x2_start_points_l + xi2l * t2_end_points_l
# Note: here we are cheating, and using the characteristic speed of the middle state,
# ignoring how it varies across the rarefaction
t2_final_points_l = np.ones_like(x2_start_points_l)
t2_final_points_l = np.minimum(t2_final_points_l,
(x2_end_points_l - t2_end_points_l * xi2m) / (vsr - xi2m))
x2_final_points_l = x2_end_points_l + (t2_final_points_l - t2_end_points_l) * xi2m
for xs, xe, te in zip(x2_start_points_r, x2_end_points_r, t2_end_points_r):
ax4.plot([xs, xe], [0.0, te], 'b-')
for xs, xe, te in zip(x2_start_points_l, x2_end_points_l, t2_end_points_l):
ax4.plot([xs, xe], [0.0, te], 'g-')
for xs, xe, ts, te in zip(x2_end_points_l, x2_final_points_l, t2_end_points_l,
t2_final_points_l):
ax4.plot([xs, xe], [ts, te], 'g-')
# Highlight the shocks
ax4.plot([0.0, vsl], [0.0, 1.0], 'r-', linewidth=2)
ax4.plot([0.0, vsr], [0.0, 1.0], 'r-', linewidth=2)
# Unphysical shock
if not physical_solution:
x_fill = []
if xi1l < xi1m:
x_fill = [xi1l, 0.0, xi1m]
elif xi1l < vsl:
x_fill = [xi1l, 0.0, vsl]
elif vsl < xi1m:
x_fill = [vsl, 0.0, xi1m]
if len(x_fill) > 0:
t_fill = [1.0, 0.0, 1.0]
ax4.fill_between(x_fill, t_fill, 1.0, facecolor = 'red', alpha = 0.5)
x_fill = []
if xi2r > xi2m:
x_fill = [xi2m, 0.0, xi2r]
elif xi2m < vsr:
x_fill = [xi2m, 0.0, vsr]
elif vsr < xi2r:
x_fill = [vsr, 0.0, xi2r]
if len(x_fill) > 0:
t_fill = [1.0, 0.0, 1.0]
ax4.fill_between(x_fill, t_fill, 1.0, facecolor = 'red', alpha = 0.5)
ax4.set_xlabel(r"$x$")
ax4.set_ylabel(r"$t$")
ax4.set_title("2-characteristics")
ax4.set_xbound(left_end, right_end)
fig.tight_layout()
In [18]:
interactive(plot_sw_all_shock_solution,
hl = FloatSlider(min = 0.1, max = 10.0, value = 1.0),
ul = FloatSlider(min = -1.0, max = 1.0, value = 0.2),
hr = FloatSlider(min = 0.1, max = 10.0, value = 1.0),
ur = FloatSlider(min = -1.0, max = 1.0, value = -0.2))
The all shock solution illustrates how the full solution can be obtained. We know that
So, using the results above, we can find the full solution to the Riemann problem by solving the nonlinear algebraic root-finding problem $$ \Phi \left( h_m ; {\bf w}_l, {\bf w}_r \right) = 0, $$ where $$ \Phi \left( h_m ; {\bf w}_l, {\bf w}_r \right) = \Phi_l \left( h_m ; {\bf w}_l \right) - \Phi_r \left( h_m ; {\bf w}_r \right), $$ and $$ \begin{align} \Phi_l & = u_m \left( h_m ; {\bf w}_l \right) & \Phi_r & = u_m \left( h_m ; {\bf w}_r \right) \\ & = \begin{cases} u_l + 2 \left( \sqrt{h_l} - \sqrt{h_m} \right) & h_l > h_m \\ u_l - \sqrt{\tfrac{1}{2} \left( h_l - h_m \right) \left( \frac{h_l}{h_m} - \frac{h_m}{h_l} \right)} & h_l < h_m \end{cases} & & = \begin{cases} u_r - 2 \left( \sqrt{h_r} - \sqrt{h_m} \right) & h_r > h_m \\ u_r + \sqrt{\tfrac{1}{2} \left( h_r - h_m \right) \left( \frac{h_r}{h_m} - \frac{h_m}{h_r} \right)} & h_r < h_m \end{cases}. \end{align} $$
In [19]:
def plot_sw_Riemann_curves(hl, ul, hr, ur):
"Plot the solution curves for states (hl, ul) and (hr, ur)"
from scipy.optimize import brentq
def phi(hstar):
"Function defining the root"
if hl < hstar:
phi_l = ul - np.sqrt(0.5 * (hl - hstar) * (hl / hstar - hstar / hl))
else:
phi_l = ul + 2.0 * (np.sqrt(hl) - np.sqrt(hstar))
if hr < hstar:
phi_r = ur + np.sqrt(0.5 * (hr - hstar) * (hr / hstar - hstar / hr))
else:
phi_r = ur - 2.0 * (np.sqrt(hr) - np.sqrt(hstar))
return phi_l - phi_r
hm = brentq(phi, 0.1 * h_min, 10.0 * h_max)
if hl < hm:
um = ul - np.sqrt(0.5 * (hl - hm) * (hl / hm - hm / hl))
else:
um = ul + 2.0 * (np.sqrt(hl) - np.sqrt(hm))
h_maximum = np.max([h_max, hl, hr, hm])
h_minimum = np.min([h_min, hl, hr, hm])
u_maximum = np.max([u_max, ul, ur, um])
u_minimum = np.min([u_min, ul, ur, um])
dh = h_maximum - h_minimum
du = u_maximum - u_minimum
# Now plot the rarefaction and shock curves as appropriate
# Here we only plot the physical pieces.
h1_shock = np.linspace(hl, h_max)
u1_shock = ul - np.sqrt(0.5 * (hl - h1_shock) * (hl / h1_shock - h1_shock / hl))
h2_shock = np.linspace(hr, h_max)
u2_shock = ur + np.sqrt(0.5 * (hr - h2_shock) * (hr / h2_shock - h2_shock / hr))
h1_rarefaction = np.linspace(h_min, hl)
u1_rarefaction = ul + 2.0 * (np.sqrt(hl) - np.sqrt(h1_rarefaction))
h2_rarefaction = np.linspace(h_min, hr)
u2_rarefaction = ur - 2.0 * (np.sqrt(hr) - np.sqrt(h2_rarefaction))
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax.plot(hl, ul, 'rx', markersize = 16, markeredgewidth = 3, label = r"${\bf w}_l$")
ax.plot(hr, ur, 'r+', markersize = 16, markeredgewidth = 3, label = r"${\bf w}_r$")
ax.plot(hm, um, 'ro', markersize = 16, markeredgewidth = 3, label = r"${\bf w}_m$")
ax.plot(h1_shock, u1_shock, 'b-', linewidth = 2,
label="1-shock")
ax.plot(h1_rarefaction, u1_rarefaction, 'b-.', linewidth = 2,
label="1-rarefaction")
ax.plot(h2_shock, u2_shock, 'g-', linewidth = 2,
label="2-shock")
ax.plot(h2_rarefaction, u2_rarefaction, 'g-.', linewidth = 2,
label="2-rarefaction")
ax.set_xlabel(r"$h$")
ax.set_ylabel(r"$u$")
ax.set_xbound(h_minimum - 0.1 * dh, h_maximum + 0.1 * dh)
ax.set_ybound(u_minimum - 0.1 * du, u_maximum + 0.1 * du)
ax.legend()
fig.tight_layout()
In [20]:
interactive(plot_sw_Riemann_curves,
hl = FloatSlider(min = 0.1, max = 10.0, value = 1.0),
ul = FloatSlider(min = -1.0, max = 1.0, value = 0.2),
hr = FloatSlider(min = 0.1, max = 10.0, value = 1.0),
ur = FloatSlider(min = -1.0, max = 1.0, value = -0.2))
Finally, we can plot the solution in physical space.
In [21]:
def plot_sw_Riemann_solution(hl, ul, hr, ur):
"Plot the Riemann problem solution for states (hl, ul) and (hr, ur)"
from scipy.optimize import brentq
def phi(hstar):
"Function defining the root"
if hl < hstar:
phi_l = ul - np.sqrt(0.5 * (hl - hstar) * (hl / hstar - hstar / hl))
else:
phi_l = ul + 2.0 * (np.sqrt(hl) - np.sqrt(hstar))
if hr < hstar:
phi_r = ur + np.sqrt(0.5 * (hr - hstar) * (hr / hstar - hstar / hr))
else:
phi_r = ur - 2.0 * (np.sqrt(hr) - np.sqrt(hstar))
return phi_l - phi_r
left_raref = False
left_shock = False
right_raref = False
right_shock = False
hm = brentq(phi, 0.1 * h_min, 10.0 * h_max)
if hl < hm:
um = ul - np.sqrt(0.5 * (hl - hm) * (hl / hm - hm / hl))
else:
um = ul + 2.0 * (np.sqrt(hl) - np.sqrt(hm))
h_maximum = np.max([h_max, hl, hr, hm])
h_minimum = np.min([h_min, hl, hr, hm])
u_maximum = np.max([u_max, ul, ur, um])
u_minimum = np.min([u_min, ul, ur, um])
dh = h_maximum - h_minimum
du = u_maximum - u_minimum
xi1l = ul - np.sqrt(hl)
xi1m = um - np.sqrt(hm)
xi1r = ur - np.sqrt(hr)
if hm > hl:
left_shock = True
vsl = ul - hm / (hm - hl) * np.sqrt(0.5 * (hl - hm) * (hl / hm - hm / hl))
else:
left_raref = True
hl_raref = np.linspace(hl, hm, 20)
ul_raref = ul + 2.0 * (np.sqrt(hl) - np.sqrt(hl_raref))
xil_raref = ul_raref - np.sqrt(hl_raref)
xi2r = ur + np.sqrt(hr)
xi2m = um + np.sqrt(hm)
xi2l = ul + np.sqrt(hl)
if hm > hr:
right_shock = True
vsr = ur + hm / (hm - hr) * np.sqrt(0.5 * (hr - hm) * (hr / hm - hm / hr))
else:
right_raref = True
hr_raref = np.linspace(hm, hr)
ur_raref = ur - 2.0 * (np.sqrt(hr) - np.sqrt(hr_raref))
xir_raref = ur_raref + np.sqrt(hr_raref)
xi_min = np.min([-1.0, xi1l, xi1m, xi2r, xi2m])
xi_max = np.max([1.0, xi1l, xi1m, xi2r, xi2m])
d_xi = xi_max - xi_min
h_maximum = np.max([hl, hr, hm])
h_minimum = np.min([hl, hr, hm])
d_h = h_maximum - h_minimum
u_maximum = np.max([ul, ur, um])
u_minimum = np.min([ul, ur, um])
d_u = u_maximum - u_minimum
xi = np.array([xi_min - 0.1 * d_xi])
h = np.array([hl])
u = np.array([ul])
if left_shock:
xi = np.append(xi, [vsl, vsl])
h = np.append(h, [hl, hm])
u = np.append(u, [ul, um])
else:
xi = np.append(xi, xil_raref)
h = np.append(h, hl_raref)
u = np.append(u, ul_raref)
if right_shock:
xi = np.append(xi, [vsr, vsr])
h = np.append(h, [hm, hr])
u = np.append(u, [um, ur])
else:
xi = np.append(xi, xir_raref)
h = np.append(h, hr_raref)
u = np.append(u, ur_raref)
xi = np.append(xi, [xi_max + 0.1 * d_xi])
h = np.append(h, [hr])
u = np.append(u, [ur])
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(221)
ax1.plot(xi, h, 'b-', label = "True solution")
ax1.set_ybound(h_minimum - 0.1 * d_h, h_maximum + 0.1 * d_h)
ax1.set_xlabel(r"$\xi$")
ax1.set_ylabel(r"$h$")
ax1.legend()
ax2 = fig.add_subplot(222)
ax2.plot(xi, u, 'b-', label = "True solution")
ax2.set_ybound(u_minimum - 0.1 * d_u, u_maximum + 0.1 * d_u)
ax2.set_xlabel(r"$\xi$")
ax2.set_ylabel(r"$u$")
ax2.legend()
ax3 = fig.add_subplot(223)
left_end = np.min([-1.0, 1.1*xi1l])
right_end = np.max([1.0, 1.1*xi2r])
left_edge = left_end - xi1l
right_edge = right_end - xi1r
x1_start_points_l = np.linspace(np.min([left_edge, left_end]), 0.0, 20)
x1_start_points_r = np.linspace(0.0, np.max([right_edge, right_end]), 20)
t1_end_points_l = np.ones_like(x1_start_points_l)
t1_end_points_r = np.ones_like(x1_start_points_r)
# Look for intersections
if left_shock:
t1_end_points_l = np.minimum(t1_end_points_l, x1_start_points_l / (vsl - xi1l))
x1_end_points_l = x1_start_points_l + xi1l * t1_end_points_l
if right_shock:
t1_end_points_r = np.minimum(t1_end_points_r, x1_start_points_r / (vsr - xi1r))
else:
t1_end_points_r = np.minimum(t1_end_points_r, x1_start_points_r / (xi2r - xi1r))
x1_end_points_r = x1_start_points_r + xi1r * t1_end_points_r
# Note: here we are cheating, and using the characteristic speed of the middle state,
# ignoring how it varies across the rarefaction
t1_final_points_r = np.ones_like(x1_start_points_r)
if left_shock:
t1_final_points_r = np.minimum(t1_final_points_r,
(x1_end_points_r - t1_end_points_r * xi1m) /
(vsl - xi1m))
x1_final_points_r = x1_end_points_r + (t1_final_points_r - t1_end_points_r) * xi1m
for xs, xe, te in zip(x1_start_points_l, x1_end_points_l, t1_end_points_l):
ax3.plot([xs, xe], [0.0, te], 'b-')
for xs, xe, te in zip(x1_start_points_r, x1_end_points_r, t1_end_points_r):
ax3.plot([xs, xe], [0.0, te], 'g-')
for xs, xe, ts, te in zip(x1_end_points_r, x1_final_points_r, t1_end_points_r,
t1_final_points_r):
ax3.plot([xs, xe], [ts, te], 'g-')
# Highlight the waves
if left_shock:
ax3.plot([0.0, vsl], [0.0, 1.0], 'r-', linewidth=2)
else:
ax3.plot([0.0, xi1l], [0.0, 1.0], 'r-', linewidth=2)
ax3.plot([0.0, xi1m], [0.0, 1.0], 'r-', linewidth=2)
xi = np.linspace(xi1l, xi1m, 11)
x_end_rarefaction = xi
for xe in x_end_rarefaction:
ax3.plot([0.0, xe], [0.0, 1.0], 'r--')
if right_shock:
ax3.plot([0.0, vsr], [0.0, 1.0], 'r-', linewidth=2)
else:
ax3.plot([0.0, xi2m], [0.0, 1.0], 'r-', linewidth=2)
ax3.plot([0.0, xi2r], [0.0, 1.0], 'r-', linewidth=2)
ax3.set_xlabel(r"$x$")
ax3.set_ylabel(r"$t$")
ax3.set_title("1-characteristics")
ax3.set_xbound(left_end, right_end)
ax4 = fig.add_subplot(224)
left_end = np.min([-1.0, 1.1*xi1l])
right_end = np.max([1.0, 1.1*xi2r])
left_edge = left_end - xi2l
right_edge = right_end - xi2r
x2_start_points_l = np.linspace(np.min([left_edge, left_end]), 0.0, 20)
x2_start_points_r = np.linspace(0.0, np.max([right_edge, right_end]), 20)
x2_end_points_r = x2_start_points_r + xi2r
t2_end_points_l = np.ones_like(x2_start_points_l)
t2_end_points_r = np.ones_like(x2_start_points_r)
# Look for intersections
if right_shock:
t2_end_points_r = np.minimum(t2_end_points_r, x2_start_points_r / (vsr - xi2r))
x2_end_points_r = x2_start_points_r + xi2r * t2_end_points_r
if left_shock:
t2_end_points_l = np.minimum(t2_end_points_l, x2_start_points_l / (vsl - xi2l))
else:
t2_end_points_l = np.minimum(t2_end_points_l, x2_start_points_l / (xi1l - xi2l))
x2_end_points_l = x2_start_points_l + xi2l * t2_end_points_l
# Note: here we are cheating, and using the characteristic speed of the middle state,
# ignoring how it varies across the rarefaction
t2_final_points_l = np.ones_like(x2_start_points_l)
if right_shock:
t2_final_points_l = np.minimum(t2_final_points_l,
(x2_end_points_l - t2_end_points_l * xi2m) /
(vsr - xi2m))
x2_final_points_l = x2_end_points_l + (t2_final_points_l - t2_end_points_l) * xi2m
for xs, xe, te in zip(x2_start_points_r, x2_end_points_r, t2_end_points_r):
ax4.plot([xs, xe], [0.0, te], 'b-')
for xs, xe, te in zip(x2_start_points_l, x2_end_points_l, t2_end_points_l):
ax4.plot([xs, xe], [0.0, te], 'g-')
for xs, xe, ts, te in zip(x2_end_points_l, x2_final_points_l, t2_end_points_l,
t2_final_points_l):
ax4.plot([xs, xe], [ts, te], 'g-')
# Highlight the waves
if left_shock:
ax4.plot([0.0, vsl], [0.0, 1.0], 'r-', linewidth=2)
else:
ax4.plot([0.0, xi1l], [0.0, 1.0], 'r-', linewidth=2)
ax4.plot([0.0, xi1m], [0.0, 1.0], 'r-', linewidth=2)
if right_shock:
ax4.plot([0.0, vsr], [0.0, 1.0], 'r-', linewidth=2)
else:
ax4.plot([0.0, xi2m], [0.0, 1.0], 'r-', linewidth=2)
ax4.plot([0.0, xi2r], [0.0, 1.0], 'r-', linewidth=2)
xi = np.linspace(xi2m, xi2r, 11)
x_end_rarefaction = xi
for xe in x_end_rarefaction:
ax4.plot([0.0, xe], [0.0, 1.0], 'r--')
ax4.set_xlabel(r"$x$")
ax4.set_ylabel(r"$t$")
ax4.set_title("2-characteristics")
ax4.set_xbound(left_end, right_end)
fig.tight_layout()
In [22]:
interactive(plot_sw_Riemann_solution,
hl = FloatSlider(min = 0.1, max = 10.0, value = 1.0),
ul = FloatSlider(min = -1.0, max = 1.0, value = 0.2),
hr = FloatSlider(min = 0.1, max = 10.0, value = 1.0),
ur = FloatSlider(min = -1.0, max = 1.0, value = -0.2))