In the previous few chapters we introduced several approximate Riemann solvers for different nonlinear conservation laws and showed some comparisons of the approximate solution to a single Riemann problem with the true solution, where there often appear to be significant differences. However, their intended use is as a building block for finite volume methods such as Godunov's method or high-resolution extensions.
How do these solvers impact the solution accuracy when used within a finite volume discretization? To investigate, we will use them within PyClaw to solve several test problems.
In particular, for both the shallow water equations and the Euler equations, we will apply the numerical method to solve dam break and shock tube problems, which are themselves Riemann problems and so we can compute the exact solution, but once discretized require the solution of a different set of Riemann problems at each cell interface every time step.
For the Euler equations we will also consider the Woodward-Colella blast wave problem. The initial data consists of two Riemann problems, with resulting shock waves of differing strengths. These shock waves later interact with each other.
In this chapter we include extensive sections of code in the notebook. This is meant to more easily allow the reader to use these as templates for setting up other problems in PyClaw. For the code in this chapter we use approximate Riemann solvers from PyClaw that can be found in these files:
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import numpy as np
from exact_solvers import euler
from clawpack import riemann
import matplotlib.pyplot as plt
from utils.snapshot_widgets import interact
from ipywidgets import widgets
We compare results obtained with using the high-resolution wave propagation method implemented in PyClaw (with limiters to avoid nonphysical oscillations). We compare the results obtained when combined with two different approximate Riemann solvers: the Roe solver and HLLE.
from exact_solvers import shallow_water as sw
from clawpack import riemann
from clawpack.riemann.shallow_roe_with_efix_1D_constants import depth, momentum, num_eqn
def setup(riemann_solver='roe',N=20,IC='dam-break'):
from clawpack import pyclaw
if riemann_solver.lower() == 'roe':
rs = riemann.shallow_roe_with_efix_1D
elif riemann_solver.lower() == 'hlle':
rs = riemann.shallow_hlle_1D
solver = pyclaw.ClawSolver1D(rs)
solver.bc_lower[0] = pyclaw.BC.extrap
solver.bc_upper[0] = pyclaw.BC.extrap
xlower = -5.0
xupper = 5.0
x = pyclaw.Dimension(xlower,xupper,N,name='x')
domain = pyclaw.Domain(x)
state = pyclaw.State(domain,num_eqn)
# Gravitational constant
state.problem_data['grav'] = 1.0
state.problem_data['dry_tolerance'] = 1e-3
state.problem_data['sea_level'] = 0.0
xc = state.grid.x.centers
x0=0.
hl = 10.
ul = 0.
hr = 0.5
ur = 0.
state.q[depth,:] = hl * (xc <= x0) + hr * (xc > x0)
state.q[momentum,:] = hl*ul * (xc <= x0) + hr*ur * (xc > x0)
claw = pyclaw.Controller()
claw.keep_copy = True
claw.output_format = None
claw.tfinal = 1.0
claw.solution = pyclaw.Solution(state,domain)
claw.solver = solver
return claw
N = 40 # number of grid cells to use
q_l = [10,0]
q_r = [0.5,0]
# Roe solution
roe_sw = setup(riemann_solver='roe',N=N)
roe_sw.verbosity = 0
status = roe_sw.run()
xc_sw = roe_sw.grid.x.centers
# HLLE solution
hlle_sw = setup(riemann_solver='hlle',N=N)
hlle_sw.verbosity = 0
status = hlle_sw.run()
# Exact solution
xc_exact_sw = np.linspace(-5,5,2000)
states_sw, speeds_sw, reval_sw, wave_types_sw = \
sw.exact_riemann_solution(q_l, q_r)
def plot_frame(i):
t = roe_sw.frames[i].t+1.e-13
fig, ax = plt.subplots(2,1, sharex=True, figsize=(8,5))
variablenames = ["Depth", "Momentum"]
variables = [depth, momentum]
ylims = [[-1,12], [-5,12]]
plt.subplots_adjust(hspace=0)
ax[0].title.set_text('Solutions at t={:.2f}'.format(t))
ax[0].set_xlim((-5,5))
ax[1].set(xlabel = 'x')
for j, variable in enumerate(variables):
ax[j].set_ylim(ylims[j])
ax[j].plot(xc_exact_sw,reval_sw(xc_exact_sw/(t+1.e-16))[j],'-k',lw=1)
ax[j].plot(xc_sw, hlle_sw.frames[i].q[variable,:],'-ob',lw=2,markersize=4)
ax[j].plot(xc_sw, roe_sw.frames[i].q[variable,:],'-or',lw=0.5,markersize=3)
ax[j].legend(['Exact','HLLE','Roe'],loc='best')
ax[j].set(ylabel=variablenames[j])
plt.show()
interact(plot_frame, i=widgets.IntSlider(min=0,max=10,value=5));
Clearly, both solvers give very similar results for this problem. They converge to the same solution, as you can verify by increasing the value of $N$ above in the live notebook. You can also try varying the left and right states.
We perform a similar test for the Euler equations with shock tube initial data.
We first consider the classic shocktube problem proposed by Sod and already discussed in Euler.html. This is a particular Riemann problem in which the initial velocity is zero on both sides of a discontinuity in pressure and/or density of a gas, and so the exact Riemann solver for the Euler equations would provide the exact solution for all time, consisting of a right-going shock, a left-going rarefaction, and an intermediate contact discontinuity.
In the numerical experiments done in this notebook, we use this initial data for a more general finite volume method that could be used to approximate the solution for any initial data. In the first time step there is a single cell interface with nontrivial Riemann data, but as the solution evolves on the grid the Riemann problems that arise in subsequent time steps are very different from the single problem we started with. Depending on the accuracy of the numerical method, the resolution of the grid, and the choice of approximate Riemann solver to use at each grid cell every time step, the numerical solution may deviate significantly from the exact solution to the original shocktube problem. This makes a good initial test problem for numerical methods because the exact solution can be computed for comparison purposes, and because it clearly shows whether the method introduces oscillations around discontinuities and/or smears them out.
from clawpack.riemann.euler_with_efix_1D_constants \
import density, momentum, energy, num_eqn
def shocktube(q_l, q_r, N=50, riemann_solver='HLL',
solver_type='classic'):
from clawpack import pyclaw
from clawpack import riemann
if riemann_solver == 'Roe':
rs = riemann.euler_1D_py.euler_roe_1D
elif riemann_solver == 'HLL':
rs = riemann.euler_1D_py.euler_hll_1D
if solver_type == 'classic':
solver = pyclaw.ClawSolver1D(rs)
solver.limiters = pyclaw.limiters.tvd.MC
else:
solver = pyclaw.SharpClawSolver1D(rs)
solver.kernel_language = 'Python'
solver.bc_lower[0]=pyclaw.BC.extrap
solver.bc_upper[0]=pyclaw.BC.extrap
x = pyclaw.Dimension(-1.0,1.0,N,name='x')
domain = pyclaw.Domain([x])
state = pyclaw.State(domain,num_eqn)
gamma = 1.4
state.problem_data['gamma']= gamma
state.problem_data['gamma1']= gamma-1.
state.problem_data['efix'] = False
xc = state.grid.p_centers[0]
velocity = (xc<=0)*q_l[1] + (xc>0)*q_r[1]
pressure = (xc<=0)*q_l[2] + (xc>0)*q_r[2]
state.q[density ,:] = (xc<=0)*q_l[0] + (xc>0)*q_r[0]
state.q[momentum,:] = velocity * state.q[density,:]
state.q[energy ,:] = pressure/(gamma - 1.) + \
0.5 * state.q[density,:] * velocity**2
claw = pyclaw.Controller()
claw.tfinal = 0.5
claw.solution = pyclaw.Solution(state,domain)
claw.solver = solver
claw.num_output_times = 10
claw.keep_copy = True
claw.verbosity=0
return claw
N = 80 # number of grid cells to use
prim_l = [1.,0.,1.]
prim_r = [1./8,0.,1./10]
q_l = euler.conservative_to_primitive(*prim_l)
q_r = euler.conservative_to_primitive(*prim_r)
# Roe-based solution
roe_st = shocktube(q_l,q_r,N=N,riemann_solver='Roe')
roe_st.run()
xc_st = roe_st.solution.state.grid.p_centers[0]
# HLL-based solution
hll_st = shocktube(q_l,q_r,N=N,riemann_solver='HLL')
hll_st.run()
# Exact solution
xc_exact_st = np.linspace(-1,1,2000)
states_st, speeds_st, reval_st, wave_types_st = euler.exact_riemann_solution(prim_l, prim_r)
def plot_frame(i):
t = roe_st.frames[i].t
fig, ax = plt.subplots(3,1, sharex=True, figsize=(8,6))
variablenames = ["Density", "Momentum", "Energy"]
variables = [density, momentum, energy]
ylims = [[0,1.1], [-0.05,0.35], [0,1.1]]
plt.subplots_adjust(hspace=0)
ax[0].title.set_text('Solutions at t={:.2f}'.format(t))
ax[0].set_xlim((-1,1))
ax[2].set(xlabel = 'x')
for j, variable in enumerate(variables):
ax[j].set_ylim(ylims[j])
ax[j].plot(xc_exact_st,reval_st(xc_exact_st/(t+1.e-16))[j],'-k',lw=1)
ax[j].plot(xc_st,hll_st.frames[i].q[variable,:],'-ob',lw=2,markersize=4)
ax[j].plot(xc_st,roe_st.frames[i].q[variable,:],'-or',lw=0.5,markersize=3)
ax[j].legend(['Exact','HLL','Roe'],loc='best')
ax[j].set(ylabel=variablenames[j])
plt.show()
interact(plot_frame, i=widgets.IntSlider(min=0, max=10, value=5, description='Frame'));
As you might expect, the HLL solver smears the middle wave (contact discontinuity) significantly more than the Roe solver does. Perhaps surprisingly, it captures the shock just as accurately as the Roe solver does. In the live notebook you can refine the grid and observe the convergence.
Next we look at the difference between the HLL and Roe solution when these solvers are employed within a higher-order method of lines discretization using fifth-order WENO and a 4th-order Runge-Kutta scheme.
N = 80 # number of grid cells to use
prim_l = [1.,0.,1.]
prim_r = [1./8,0.,1./10]
q_l = euler.conservative_to_primitive(*prim_l)
q_r = euler.conservative_to_primitive(*prim_r)
roe_weno = shocktube(q_l,q_r,N=N,riemann_solver='Roe',solver_type='sharpclaw')
roe_weno.run()
hll_weno = shocktube(q_l,q_r,N=N,riemann_solver='HLL',solver_type='sharpclaw')
hll_weno.run()
xc_weno = roe_weno.solution.state.grid.p_centers[0]
# Exact solution
xc_exact_weno = np.linspace(-1,1,2000)
states_weno, speeds_weno, reval_weno, wave_types_weno = euler.exact_riemann_solution(prim_l, prim_r)
def plot_frame(i):
t = roe_weno.frames[i].t
fig, ax = plt.subplots(3,1, sharex=True, figsize=(8,6))
variablenames = ["Density", "Momentum", "Energy"]
variables = [density, momentum, energy]
ylims = [[0,1.1], [-0.05,0.35], [0,1.1]]
plt.subplots_adjust(hspace=0)
ax[0].title.set_text('Solutions at t={:.2f}'.format(t))
ax[0].set_xlim((-1,1))
ax[2].set(xlabel = 'x')
for j, variable in enumerate(variables):
ax[j].set_ylim(ylims[j])
ax[j].plot(xc_exact_weno,reval_weno(xc_exact_weno/(t+1.e-16))[j],'-k',lw=1)
ax[j].plot(xc_weno,hll_weno.frames[i].q[variable,:],'-ob',lw=2,markersize=4)
ax[j].plot(xc_weno,roe_weno.frames[i].q[variable,:],'-or',lw=0.5,markersize=3)
ax[j].legend(['Exact','HLL','Roe'],loc='best')
ax[j].set(ylabel=variablenames[j])
plt.show()
interact(plot_frame, i=widgets.IntSlider(min=0, max=10, value=5, description='Frame'));
With higher-order discretizations, the difference in solutions due to using different Riemann solvers is less significant. This is partly because these high-order schemes use more accurate values as inputs to the Riemann problem, so that in smooth regions the jump between most cells is very small.
Finally, we consider the Woodward-Colella blast wave problem, which is discussed for example in (LeVeque 2002). Here the initial velocity is zero and the density is one everywhere. The pressure is \begin{align} p_0(x) = \begin{cases} 1000 & 0 \le x \le 0.1 \\ 0.01 & 0.1 \le x \le 0.9 \\ 100 & 0.9 \le x \le 1 \end{cases} \end{align} The boundaries at $x=0$ and $x=1$ are solid walls. The solution involves a Riemann problem at $x=0.1$ and another at $x=0.9$. Later, the waves resulting from these Riemann problems interact with each other.
from clawpack.riemann.euler_with_efix_1D_constants \
import density, momentum, energy, num_eqn
def blastwave(N=400, riemann_solver='HLL', solver_type='classic'):
from clawpack import pyclaw
from clawpack import riemann
if riemann_solver == 'Roe':
kernel_language = 'Fortran'
rs = riemann.euler_with_efix_1D
elif riemann_solver == 'HLL':
kernel_language = 'Python'
rs = riemann.euler_1D_py.euler_hll_1D
if solver_type == 'classic':
solver = pyclaw.ClawSolver1D(rs)
solver.limiters = pyclaw.limiters.tvd.MC
else:
solver = pyclaw.SharpClawSolver1D(rs)
solver.kernel_language = kernel_language
solver.bc_lower[0]=pyclaw.BC.wall
solver.bc_upper[0]=pyclaw.BC.wall
x = pyclaw.Dimension(0.0,1.0,N,name='x')
domain = pyclaw.Domain([x])
state = pyclaw.State(domain,num_eqn)
gamma = 1.4
state.problem_data['gamma']= gamma
state.problem_data['gamma1']= gamma-1.
state.problem_data['efix'] = False
xc = state.grid.p_centers[0]
pressure = (xc<0.1)*1.e3 + (0.1<=xc)*(xc<0.9)*1.e-2 + (0.9<=xc)*1.e2
state.q[density ,:] = 1.
state.q[momentum,:] = 0.
state.q[energy ,:] = pressure / (gamma - 1.)
claw = pyclaw.Controller()
claw.tfinal = 0.038
claw.solution = pyclaw.Solution(state,domain)
claw.solver = solver
claw.num_output_times = 30
claw.keep_copy = True
claw.verbosity=0
return claw
N = 400 # number of grid cells to use
roe_bw = blastwave(N=N, riemann_solver='Roe')
roe_bw.run()
hll_bw = blastwave(N=N, riemann_solver='HLL')
hll_bw.run()
fine_bw = blastwave(N=4000,riemann_solver='Roe')
fine_bw.run();
xc_bw = roe_bw.solution.state.grid.p_centers[0]
xc_fine_bw = fine_bw.solution.state.grid.p_centers[0]
def plot_frame(i):
t = roe_bw.frames[i].t
fig, ax = plt.subplots(3,1, sharex=True, figsize=(8,6))
variablenames = ["Density", "Momentum", "Energy"]
variables = [density, momentum, energy]
ylims = [[-1,10], [-50,120], [-300,3000]]
plt.subplots_adjust(hspace=0)
ax[0].title.set_text('Solutions at t={:.3f}'.format(t))
ax[0].set_xlim((0,1))
ax[2].set(xlabel = 'x')
for j, variable in enumerate(variables):
ax[j].set_ylim(ylims[j])
ax[j].plot(xc_fine_bw,fine_bw.frames[i].q[variable,:],'-k',lw=1)
ax[j].plot(xc_bw,hll_bw.frames[i].q[variable,:],'-b',lw=2)
ax[j].plot(xc_bw,roe_bw.frames[i].q[variable,:],'--r',lw=2)
ax[j].legend(['Fine','HLL','Roe'],loc='best')
ax[j].set(ylabel=variablenames[j])
plt.show()
interact(plot_frame, i=widgets.IntSlider(min=0, max=30, value=15, description='Frame'));
Here no exact solution is available, so we compare with a solution computed on a finer grid. Again the solutions are fairly similar, though the HLL solution is a bit more smeared.
One should not conclude from these tests that, for instance, the Roe solver is always better than the HLL solver. Many factors besides accuracy should be considered, including cost and robustness. As we have seen the HLL solver is more robust in the presence of near-vacuum states.