In Euler_approximate.ipynb we introduced several approximate Riemann solvers for the Euler equations of compressible gas dynamics. 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 two standard test problems for one-dimensional compressible flow:
In this chapter, unlike previous chapters, 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. For more information about the software and algorithms used here, see this paper and references therein.
In [ ]:
%matplotlib inline
In [ ]:
%config InlineBackend.figure_format = 'svg'
import numpy as np
from exact_solvers import euler
from clawpack import riemann
import matplotlib.pyplot as plt
from ipywidgets import interact
from ipywidgets import widgets
In [ ]:
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
In [ ]:
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
We compare results obtained with these two solvers within a second-order Lax-Wendroff based scheme (with limiters to avoid nonphysical oscillations). This method is the basis of Clawpack.
We first consider the classic shocktube problem proposed by Sod and already discussed in Euler.ipynb. 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.
In [ ]:
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=50,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=50,riemann_solver='HLL')
hll_st.run()
# Exact solution
xc_exact_st = np.linspace(-1,1,2000)
states, speeds, reval, wave_types = euler.exact_riemann_solution(prim_l, prim_r)
def plot_frame(i):
t = roe_st.frames[i].t
fig, ax = plt.subplots(figsize=(8,4))
ax.set_xlim((-1,1)); ax.set_ylim((0,1.1))
ax.plot(xc_exact_st,reval(xc_exact_st/(t+1.e-16))[0],'-k',lw=1)
ax.plot(xc_st,hll_st.frames[i].q[density,:],'-ob',lw=2)
ax.plot(xc_st,roe_st.frames[i].q[density,:],'-or',lw=2)
plt.legend(['Exact','HLL','Roe'],loc='best')
plt.title('Density at t={:.2f}'.format(t))
plt.show()
interact(plot_frame, i=widgets.IntSlider(min=0, max=10, 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.
Next 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.
In [ ]:
roe_bw = blastwave(riemann_solver='Roe')
roe_bw.run()
hll_bw = blastwave(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(figsize=(8,4))
ax.set_xlim((0.,1)); ax.set_ylim((0,10))
ax.plot(xc_fine_bw,fine_bw.frames[i].q[density,:],'-k',lw=1)
ax.plot(xc_bw,hll_bw.frames[i].q[density,:],'-ob',lw=2)
ax.plot(xc_bw,roe_bw.frames[i].q[density,:],'-or',lw=2)
plt.legend(['Fine','HLL','Roe'],loc='best')
plt.title('Density at t={:.3f}'.format(t))
plt.show()
interact(plot_frame, i=widgets.IntSlider(min=0, max=30, 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 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.
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.
In [ ]:
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=50,riemann_solver='Roe',solver_type='sharpclaw')
roe_weno.run()
hll_weno = shocktube(q_l,q_r,N=50,riemann_solver='HLL',solver_type='sharpclaw')
hll_weno.run()
xc = roe_weno.solution.state.grid.p_centers[0]
# Exact solution
xc_exact = np.linspace(-1,1,2000)
states, speeds, reval, wave_types = euler.exact_riemann_solution(prim_l, prim_r)
def plot_frame(i):
t = roe_weno.frames[i].t
fig, ax = plt.subplots(figsize=(8,4))
ax.set_xlim((-1,1)); ax.set_ylim((0,1.1))
ax.plot(xc_exact,reval(xc_exact/(t+1.e-16))[0],'-k',lw=1)
ax.plot(xc,hll_weno.frames[i].q[density,:],'-ob',lw=2)
ax.plot(xc,roe_weno.frames[i].q[density,:],'-or',lw=2)
plt.legend(['Exact','HLL','Roe'],loc='best')
plt.title('Density at t={:.2f}'.format(t))
plt.show()
interact(plot_frame, i=widgets.IntSlider(min=0, max=10, 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.