Comparison of full finite volume discretizations using approximate Riemann solvers

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:

In [1]:
%matplotlib inline
In [2]:
%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
Will create static figures with single value of parameters

Shallow water equations

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.

In [3]:
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
In [4]:
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.

Euler equations

We perform a similar test for the Euler equations with shock tube initial data.

Sod shock tube problem

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.

In [5]:
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 [6]:
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.

High-order WENO + Runge-Kutta

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 [7]:
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.

Woodward-Colella blast wave

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.

In [8]:
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
In [9]:
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.