# advection_highresÂ¶

R.J. LeVeque, AMath 574, Winter Quarter 2023

## Advection equation using high resolution methodsÂ¶

This notebook has a function `highres_step`

that currently implements upwind, Lax-Wendroff, and the high-resolution method with the minmod limiter.

InÂ [1]:

```
%matplotlib inline
```

InÂ [2]:

```
from pylab import *
```

InÂ [3]:

```
from IPython.display import HTML
```

InÂ [4]:

```
try:
from clawpack.visclaw import animation_tools
except:
print("Failed to load animation_tools from Clawpack")
```

Hard-wire the advection velocity:

InÂ [5]:

```
u = 2.
```

We now set up a grid with 2 ghost cells on each side, as needed in general for high-resolution methods with limiters.

Note that only one ghost cell is needed for the data that goes into the Riemann problems we need solve for the Godunov update in the interior, but in order to limit the resulting waves based on looking at the neighboring Riemann problem on the upwind side, we need to solve an additional Riemann problem requiring another ghost cell for its data.

InÂ [6]:

```
xlower = 0.
xupper = 1.
num_cells = 20
dx = (xupper - xlower)/num_cells
# cell centers, including two ghost cells on either side:
x = arange(xlower-3*dx/2, xupper+2*dx, dx)
print('Including 4 ghost cells, the grid has %i cells indexed 0 to %i' % (len(x), len(x)-1))
print('cell centers: \n',x)
# interior grid cells are numbered 2,3,...num_cells+1
# It's also useful to define the array of indices of
# interior cells for vectorized methods:
ii = array(range(2,num_cells+2), dtype=int)
print('The indices of interior cells:\n ',ii)
```

InÂ [7]:

```
dt = 1.*dx/u
cfl = u*dt/dx
print('The Courant number is %.3f' % cfl)
```

The Courant number is 1.000

### A simple illustration of a vectorized method:Â¶

InÂ [8]:

```
def upwind_step(x,Qn,dt,method='upwind'):
"""
Illustrates simple vectorized version of the upwind method.
The method parameter is not used, just there for consistency
with highres_step.
"""
num_cells = len(x) - 4
dx = x[1] - x[0] # assuming uniform grid
#indices of interior grid cells (not including ghost cells)
ii = array(range(2,num_cells+2), dtype=int)
cfl = u*dt/dx # assuming u defined globally
Qnp = Qn.copy()
# periodic BCs with two ghost cells on each side
Qn[0] = Qn[-4]
Qn[1] = Qn[-3]
Qn[-2] = Qn[2]
Qn[-1] = Qn[3]
# vectorized version of upwind:
# (equivalent to looping `for i in range(2,num_cells+2)`)
if u>0:
Qnp[ii] = Qn[ii] - cfl*(Qn[ii] - Qn[ii-1])
else:
Qnp[ii] = Qn[ii] - cfl*(Qn[ii+1] - Qn[ii])
return Qnp
```

### The version that allows adding second-order corrections (perhaps with a limiter)Â¶

InÂ [9]:

```
def highres_step(x,Qn,dt,method):
"""
High-resolution method for the advection equation, written
in a more general form similar to the methods described in
Section 12.8 of FVMHP for nonlinear scalar equations.
method can be
'upwind', 'LW', or 'minmod'
in this version.
"""
num_cells = len(Qn) - 4
dx = x[1] - x[0] # assuming uniform grid
#indices of interior grid cells (not including ghost cells)
ii = array(range(2,num_cells+2), dtype=int)
Qnp = Qn.copy()
# periodic BCs with two ghost cells on each side
Qn[0] = Qn[-4]
Qn[1] = Qn[-3]
Qn[-2] = Qn[2]
Qn[-1] = Qn[3]
# waves:
wave = zeros(x.shape)
wave[1:] = Qn[1:] - Qn[0:-1] # for cells that have an interface to left
# speeds:
s = u*ones(x.shape) # all the same for constant-coefficient advection
# fluctuations (12.8) in FVMHP
amdq = where(s<0, s*wave, 0.)
apdq = where(s>0, s*wave, 0.)
# Godunov step: propagate waves to update the proper cell:
Qnp[ii] = Qnp[ii] - dt/dx * apdq[ii] # right-going
Qnp[ii] = Qnp[ii] - dt/dx * amdq[ii+1] # left-going
if method == 'upwind':
# for upwind, zero out all waves in high-resolution correction terms:
wlimiter = zeros(wave.shape)
elif method == 'LW':
# Lax-Wendroff, use full wave:
wlimiter = ones(wave.shape)
else:
# compute things needed for any limiter:
# the wave from the neighboring Riemann problem on upwind side:
wave_upwind = zeros(x.shape)
wave_upwind[2:-1] = where(s[2:-1]>0, wave[1:-2], wave[3:])
# compute ratio theta = wave_upwind/wave, returning 0 where wave==0:
# using the numpy.divide function this way avoids divide-by-zero
theta = divide(wave_upwind, wave, where=wave!=0, out=zeros(wave.shape))
if method == 'minmod':
wlimiter = maximum(0., minimum(1., theta))
else:
print('method %s is not implemented, using upwind' % method)
wlimiter = zeros(wave.shape)
Wtilde = wlimiter * wave
# correction fluxes:
Ftilde = zeros(x.shape)
Ftilde[1:-1] = 0.5*abs(s[1:-1])*(1 - abs(s[1:-1])*dt/dx) * Wtilde[1:-1]
# update cells by flux differencing Ftilde:
Qnp[ii] = Qnp[ii] - dt/dx * (Ftilde[ii+1] - Ftilde[ii])
return Qnp
```

InÂ [10]:

```
def time_stepper(t0, x, Q0, dt, nsteps, one_step=highres_step, method='minmod'):
"""
Take nsteps with time step dt, starting with initial data Q0 at time t0.
To take a single step use the method specified by one_step, which
defaults to Godunov_step since that is the only method currently defined.
But you might want to define a new method for comparison.
"""
Qn = Q0.copy()
for n in range(nsteps):
Qn = one_step(x, Qn, dt, method)
return Qn
```

InÂ [11]:

```
def plotQ(x, Qn, tn):
# only plot the interior points:
plot(x[2:-3], Qn[2:-3], 'bo-', markersize=4)
grid(True)
xlim(xlower, xupper)
title('Time t = %.3f' % tn)
```

InÂ [12]:

```
def make_anim(t0, x, Q0, dt, nsteps, nplot,
one_step=highres_step, method='minmod'):
Qn = Q0.copy()
figsize = (6,2)
figs = [] # to accumulate figures for animation
# plot initial data:
fig = figure(figsize=figsize)
plotQ(x,Q0,t0)
title('t = %.3f, %s' % (t0,method))
figs.append(fig)
close(fig)
for n in range(1,nsteps+1):
# take the next step
Qn = one_step(x,Qn,dt,method)
tn = n*dt
if mod(n,nplot)==0:
fig = figure(figsize=figsize)
plotQ(x,Qn,tn)
title('t = %.3f, %s' % (tn,method))
figs.append(fig)
close(fig)
anim = animation_tools.animate_figs(figs, figsize=figsize)
return anim
```

### Example usageÂ¶

Set up to use the minmod limiter method.

InÂ [13]:

```
num_cells = 100
dx = (xupper - xlower)/num_cells
# cell centers, including two ghost cells on either side:
x = arange(xlower-3*dx/2, xupper+2*dx, dx)
t0 = 0.
Q0 = where(x<0.3, 1., 0.) + exp(-200*(x-0.7)**2)
tn = 0.
figure(figsize=(6,2))
plotQ(x,Q0,t0)
dt = 0.004
cfl = u*dt/dx
nsteps = 20
tn = t0 + nsteps*dt
print('Using dt = %.4f, cfl = %.2f, taking %i steps to time %.3f' \
% (dt,cfl,nsteps,tn))
Qn = time_stepper(t0, x, Q0, dt, nsteps, highres_step, 'minmod')
fig = figure(figsize=(6,2))
plotQ(x,Qn,tn)
```

Using dt = 0.0040, cfl = 0.80, taking 20 steps to time 0.080

### Make an animationÂ¶

InÂ [14]:

```
num_cells = 100
dx = (xupper - xlower)/num_cells
# cell centers, including two ghost cells on either side:
x = arange(xlower-3*dx/2, xupper+2*dx, dx)
t0 = 0.
Q0 = where(x<0.3, 1., 0.) + exp(-200*(x-0.7)**2)
dt = 0.0045
cfl = u*dt/dx
print('Courant number = %.2f' % cfl)
nsteps = 150
nplot = 10
anim = make_anim(t0, x, Q0, dt, nsteps, nplot, highres_step, 'minmod')
HTML(anim.to_jshtml())
```

Courant number = 0.90

Out[14]: