How do these solvers impact the solution accuracy when used within a finite volume discretization?
%matplotlib inline
import numpy as np
from exact_solvers import Euler
from clawpack import riemann
from utils import riemann_tools
import matplotlib.pyplot as plt
from collections import namedtuple
from ipywidgets import interact
from ipywidgets import widgets
Primitive_State = namedtuple('State', Euler.primitive_variables)
gamma = 1.4
problem_data = {}
problem_data['gamma'] = gamma
problem_data['gamma1'] = gamma - 1.0
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)
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(q_l, q_r, N=400, 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)
else:
solver = pyclaw.SharpClawSolver1D(rs)
solver.time_integrator = 'SSP33'
solver.cfl_max = 0.65
solver.cfl_desired = 0.6
solver.lim_type = 1
#solver.char_decomp = 2
solver.kernel_language = 'Python'
solver.bc_lower[0]=pyclaw.BC.extrap
solver.bc_upper[0]=pyclaw.BC.extrap
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]
state.q[density ,:] = 1.
state.q[momentum,:] = 0.
state.q[energy ,:] = ( (xc<0.1)*1.e3 + (0.1<=xc)*(xc<0.9)*1.e-2 + (0.9<=xc)*1.e2 ) / (gamma - 1.)
claw = pyclaw.Controller()
claw.tfinal = 0.038
claw.solution = pyclaw.Solution(state,domain)
claw.solver = solver
claw.num_output_times = 10
claw.keep_copy = True
claw.verbosity=0
return claw
q_l = Euler.conservative_to_primitive(1.,0.,1.)
q_r = Euler.conservative_to_primitive(1./8,0.,1./10)
roe = shocktube(q_l,q_r,N=50,riemann_solver='Roe')
roe.run()
hll = shocktube(q_l,q_r,N=50,riemann_solver='HLL')
hll.run()
fine = shocktube(q_l,q_r,N=1000,riemann_solver='Roe')
fine.run();
xc = roe.solution.state.grid.p_centers[0]
xc_fine = fine.solution.state.grid.p_centers[0]
def plot_frame(i):
fig, ax = plt.subplots(figsize=(12,6))
ax.set_xlim((-1,1)); ax.set_ylim((0,1.1))
ax.plot(xc_fine,fine.frames[i].q[density,:],'-k',lw=1)
ax.plot(xc,hll.frames[i].q[density,:],'-ob',lw=2)
ax.plot(xc,roe.frames[i].q[density,:],'-or',lw=2)
plt.legend(['Fine','HLL','Roe'],loc='best')
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.
q_l = Euler.conservative_to_primitive(1.,0.,1.)
q_r = Euler.conservative_to_primitive(1./8,0.,1./10)
roe = blastwave(q_l,q_r,riemann_solver='Roe')
roe.run()
hll = blastwave(q_l,q_r,riemann_solver='HLL')
hll.run()
fine = blastwave(q_l,q_r,N=2000,riemann_solver='Roe')
fine.run();
xc = roe.solution.state.grid.p_centers[0]
xc_fine = fine.solution.state.grid.p_centers[0]
def plot_frame(i):
fig, ax = plt.subplots(figsize=(12,6))
ax.set_xlim((0.,1)); ax.set_ylim((0,20))
ax.plot(xc_fine,fine.frames[i].q[density,:],'-k',lw=1)
ax.plot(xc,hll.frames[i].q[density,:],'-ob',lw=2)
ax.plot(xc,roe.frames[i].q[density,:],'-or',lw=2)
plt.legend(['Fine','HLL','Roe'],loc='best')
plt.show()
interact(plot_frame, i=widgets.IntSlider(min=0, max=10, description='Frame'));
q_l = Euler.conservative_to_primitive(1.,0.,1.)
q_r = Euler.conservative_to_primitive(1./8,0.,1./10)
roe2 = shocktube(q_l,q_r,N=50,riemann_solver='Roe',solver_type='sharpclaw')
roe2.run()
hll2 = shocktube(q_l,q_r,N=50,riemann_solver='HLL',solver_type='sharpclaw')
hll2.run()
fine2 = shocktube(q_l,q_r,N=1000,riemann_solver='Roe',solver_type='sharpclaw')
fine2.run();
xc = roe2.solution.state.grid.p_centers[0]
xc_fine2 = fine2.solution.state.grid.p_centers[0]
def plot_frame(i):
fig, ax = plt.subplots(figsize=(12,6))
ax.set_xlim((-1,1)); ax.set_ylim((0,1.1))
ax.plot(xc_fine2,fine2.frames[i].q[density,:],'-k',lw=1)
ax.plot(xc,hll2.frames[i].q[density,:],'-ob',lw=2)
ax.plot(xc,roe2.frames[i].q[density,:],'-or',lw=2)
plt.legend(['Fine','HLL','Roe'],loc='best')
plt.show()
interact(plot_frame, i=widgets.IntSlider(min=0, max=10, description='Frame'));
With higher-order discretizations, the difference in the Riemann solvers is less significant.