This notebook shows the effect of limiters on advection for various initial conditions. It includes the results for first and second order methods from Figures 6.1-6.3 of Finite Volume Methods for Hyperbolic Problems, as well as results for WENO methods.
%matplotlib inline
import numpy as np
from clawpack import pyclaw
from clawpack import riemann
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import animation
from clawpack.visclaw.JSAnimation import IPython_display
The function below sets up an advection simulation with the specified scheme, on the unit interval with periodic boundary conditions. Since output is written at integer times, the solution at each output time has been advected precisely back to its original location. Plotting each frame shows us how the solution is modified by numerical effects as it moves once through the domain.
def setup(scheme='minmod',cfl_max=0.9,IC='gauss_square',mx=100):
if 'weno' in scheme:
solver = pyclaw.SharpClawSolver1D(riemann.advection_1D)
else:
solver = pyclaw.ClawSolver1D(riemann.advection_1D)
solver.bc_lower[0] = pyclaw.BC.periodic
solver.bc_upper[0] = pyclaw.BC.periodic
if scheme in ('minmod','superbee','MC','vanleer'):
solver.limiters = getattr(pyclaw.limiters.tvd,scheme)
#elif scheme == 'CT':
#solver.limiters = pyclaw.limiters.tvd.cada_torrilhon_limiter
elif scheme == 'Lax-Wendroff':
solver.limiters = 0
elif scheme == 'first-order':
solver.order = 1
elif 'weno' in scheme:
solver.weno_order = int(scheme[4:]) #weno5, weno7, ...
else:
raise Exception('Unrecognized limiter')
solver.cfl_max = cfl_max
solver.cfl_desired = cfl_max*0.9
x = pyclaw.Dimension('x',0.0,1.0,mx)
domain = pyclaw.Domain(x)
num_eqn = 1
state = pyclaw.State(domain,num_eqn)
state.problem_data['u']=1.
grid = state.grid
xc = grid.x.centers
if IC=='gauss_square':
beta=200.; x0=0.3
state.q[0,:] = np.exp(-beta * (xc-x0)**2) + (xc>0.6)*(xc<0.8)
elif IC=='wavepacket':
beta=100.; x0=0.5
state.q[0,:] = np.exp(-beta * (xc-x0)**2) * np.sin(80.*xc)
else:
raise Exception('Unrecognized initial condition.')
claw = pyclaw.Controller()
claw.solution = pyclaw.Solution(state,domain)
claw.solver = solver
claw.keep_copy = True
claw.output_format = None
claw.tfinal =10.0
return claw
results = []
schemes = ('first-order','Lax-Wendroff','minmod','superbee','MC','vanleer','weno5','weno7','weno9')
for scheme in schemes:
claw = setup(scheme=scheme)
claw.run()
results.append(claw.frames)
def animate(results,ymin=-0.1):
fig = plt.figure(figsize=(10,6))
N = len(results)
n = int(np.ceil(np.sqrt(N)))
axes = []
gs1 = matplotlib.gridspec.GridSpec(n, n)
gs1.update(wspace=0.,hspace=0.)
for i in range(n):
for j in range(n):
k = n*i + j
if k<N:
axes.append(plt.subplot(gs1[i,j]));
if j>0:
axes[-1].yaxis.set_ticklabels(())
lines = [0]*len(schemes)
for i in range(len(lines)):
lines[i], = axes[i].plot([], [], lw=2)
for i,ax in enumerate(axes):
ax.set_xlim(0,1); ax.set_ylim(ymin,1.3)
#ax.legend( [schemes[i] ] )
ax.set_title(schemes[i], x = 0.5, y=0.85 )
xc = results[0][0].p_centers[0]
def fplot(frame_number):
for i, line in enumerate(lines):
line.set_data(xc,results[i][frame_number].q[0,:])
return lines,
return matplotlib.animation.FuncAnimation(fig, fplot, frames=len(claw.frames), interval=30)
animate(results)
It's easy to see the differences in dissipation and dispersion of the various schemes. Of course, these properties also depend on the CFL number. Try changing the CFL number in the call to setup()
above and see how the results change.
To see how the limiters affect intermediate frequencies, let's advect a wavepacket.
results = []
for scheme in schemes:
claw = setup(scheme=scheme,IC='wavepacket',mx=300)
claw.run()
results.append(claw.frames)
animate(results,ymin=-1.2)