<jupyter_start><jupyter_text>Comparison of full discretizations using approximate Riemann solvers In [Euler_approximate_solvers.ipynb](Euler_approximate_solvers.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](http://www.clawpack.org/pyclaw/) to solve two standard test problems for one-dimensional compressible flow:  1. The Sod shocktube problem.  This is a Riemann problem, which we have considered previously.  2. 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, 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](https://peerj.com/articles/cs-68/) and references therein.<jupyter_code>%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 ipywidgets import interact
from ipywidgets import widgets
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
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<jupyter_output><empty_output><jupyter_text>Classic Clawpack algorithm with HLL versus Roe 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.    Sod shock tube problem We first consider the classic shocktube problem proposed by Sod and already discussed in [Euler.ipynb](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.<jupyter_code>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=(12,6))
    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'));<jupyter_output><empty_output><jupyter_text>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.  Woodward-Colella blast wave 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.<jupyter_code>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=(12,6))
    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'));<jupyter_output><empty_output><jupyter_text>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.  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.<jupyter_code>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=(12,6))
    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'));<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11647223518432105, 'fraction_numerical': 0.02053196453569762, 'mean_word_length': 4.164337349397591, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 20, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '13423', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3948, 'n_tokens_neox': 3690, 'n_words': 984, 'path': 'Euler_compare.ipynb', 'repo_name': 'haraldschilly/riemann_book', 'size': '10715'}