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.

%matplotlib inline

from pylab import *

from IPython.display import HTML

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

Hard-wire the advection velocity:

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.

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)

Including 4 ghost cells, the grid has 24 cells indexed 0 to 23
cell centers: 
 [-0.075 -0.025  0.025  0.075  0.125  0.175  0.225  0.275  0.325  0.375
  0.425  0.475  0.525  0.575  0.625  0.675  0.725  0.775  0.825  0.875
  0.925  0.975  1.025  1.075]
The indices of interior cells:
   [ 2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21]

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:

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)

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

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

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)

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.

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

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

Once Loop Reflect