Leapfrog method for advection with outflow boundary

AMath 586, Spring Quarter 2019 at the University of Washington. For other notebooks, see Index.ipynb or the Index of all notebooks on Github.

Sample program to solve the advection equation with the Leapfrog method to explore outflow boundary conditions.

Mathematically we cannot impose the solution at the outflow boudary, which is at $x=1$ since we assume $a>0$ and solve the advection equation on $0\leq x\leq 1$.

But the Leapfrog method cannot be applied at the right boundary and must be augmented with some numerical boundary condition. Often a 1-sided method is used instead at the right-most point instead of a centered method like Leapfrog.

The code below implements several options and results from each are discussed in this notebook:

  • Impose the true solution at the outflow boundary,
  • Set $u=0$ at the boundary,
  • Use the first-order upwind method,
  • Use a modified version of Leapfrog, using 2 time levels but one-sided in space,
  • Use the Beam-Warming method.

Similar boundary condition issues arise whenever the method used in the interior is not a one-sided method. We illustrate the issues with the Leapfrog method because this method is non-dissipative, and so any errors made at the boundary tend to quickly contaminate the interior soluiton in ways that are easy to see.

Note also that because the stencil of Leapfrog involves points to the right as well as to the left, it is possible for information in the computed solution to propagate from right to left even though the exact solution consists only of a advection from left to right. We will see below that errors at the right boundary give rise to waves that propagate into the interior.

In [1]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib
In [2]:
from matplotlib import animation
from IPython.display import HTML

Animation function

This version also plots the error beside the solution, making it easier to see when very small.

In [3]:
def make_animation(adv_input, adv_output, nplot=1):
    
    """
    Plot every `nplot` frames of the solution and turn into
    an animation.    Plot solution and error side-by-side.
    """
    xfine = linspace(adv_input.ax,adv_input.bx,1001)
    fig, ax = plt.subplots(1,2,figsize=(12,6))
    ax1, ax2 = ax

    ax1.set_xlim((adv_input.ax,adv_input.bx))
    ax1.set_ylim((-1.2, 1.2))
    ax2.set_xlim((adv_input.ax,adv_input.bx))
    ax2.set_ylim((-0.1, 0.1))

    line1, = ax1.plot([], [], '+-', color='b', lw=2, label='computed')
    line2, = ax1.plot([], [], color='r', lw=1, label='true')
    ax1.legend(loc='lower left')
    title1 = ax1.set_title('')
    
    line3, = ax2.plot([], [], '+-', color='b')
    title2 = ax2.set_title('')

    def init():
        line1.set_data(adv_output.x_computed, adv_output.u_computed[:,0])
        line2.set_data(xfine, adv_input.utrue(xfine, adv_input.t0))
        title1.set_text('time t = %8.4f' % adv_input.t0)
        #error = adv_output.u_computed[:,0] - adv_input.utrue(adv_output.x_computed, adv_output.t[0])
        error = adv_output.error[:,0]
        line3.set_data(adv_output.x_computed, error)
        title2.set_text('Error at time t = %8.4f' % adv_output.t[0])
        return (line1,line2,title1,line3,title2)

    def animate(n):
        line1.set_data(adv_output.x_computed, adv_output.u_computed[:,n])
        line2.set_data(xfine, adv_input.utrue(xfine, adv_output.t[n]))
        title1.set_text('time t = %8.4f' % adv_output.t[n])
        #error = adv_output.u_computed[:,n] - adv_input.utrue(adv_output.x_computed, adv_output.t[n])
        error = adv_output.error[:,n]
        line3.set_data(adv_output.x_computed, error)
        title2.set_text('Error at time t = %8.4f' % adv_output.t[n])
        return (line1,line2,title1,line3,title2)

    frames = range(0, len(adv_output.t), nplot) # which frames to plot
    anim = animation.FuncAnimation(fig, animate, init_func=init,
                                   frames=frames,
                                   interval=200,
                                   repeat=False,
                                   blit=True)
    close('all')  # so one last frame plot doesn't remain
    return anim

This version of the AdvectionSolutionInput class provides an attribute bc to specify what boundary condition to use at the outflow boundary, for comparison purposes.

In [4]:
class AdvectionSolutionInput(object):
    def __init__(self):
        # inputs:
        self.t0 = 0.
        self.tfinal = 1.
        self.ax = 0.
        self.bx = 1.
        self.mx = 39.
        self.utrue = None
        self.a = 1.
        self.nsteps = 10
        self.bc = 'upwind'
        
class AdvectionSolutionOutput(object):
    def __init__(self):
        # outputs:
        self.h = None
        self.dt = None
        self.t = None
        self.nu = None
        self.x_computed = None
        self.u_computed = None
        self.error = None
In [5]:
def Leapfrog_Advection(advection_solution_input):
    """
    Solve u_t + a*u_x = 0 on [ax,bx] with inflow/outflow boundary conditions,
    and the Leapfrog method.
    
    Input: 
        `advection_solution_input` should be on object of class `AdvectionSolutionInput`
        specifying inputs.
    Output:
        an object of class `AdvectionSolutionOutput` with the solution and other info.

    """

        
    # unpack the inputs for brevity:
    ax = advection_solution_input.ax
    bx = advection_solution_input.bx
    a = advection_solution_input.a
    m = advection_solution_input.mx
    utrue = advection_solution_input.utrue
    t0 = advection_solution_input.t0
    tfinal = advection_solution_input.tfinal
    nsteps = advection_solution_input.nsteps
    
    h = (bx-ax)/float(m+1)    # h = delta x
    x = linspace(ax,bx,m+2)   # note x(1)=0 and x(m+2)=1
                               # u(1)=g0 and u(m+2)=g1 are known from BC's
    dt = tfinal / float(nsteps)
    
    # initial conditions:
    u0 = utrue(x,t0)
    
    # also use for second initial vector needed at time t1:
    u1 = utrue(x, t0+dt)
          
    # initialize u:
    u_nm2 = u0
    u_nm1 = u1
    
    t = empty((nsteps+1,), dtype=float)
    error = zeros((m+2,nsteps+1,), dtype=float)
    u_computed = empty((m+2,nsteps+1), dtype=float)

    t[0] = t0
    t[1] = t0+dt

    u_computed[:,0] = u0
    u_computed[:,1] = u1
    
    u = u1 # initialize
    
    nu = a*dt/h  # Courant number
    
    # main time-stepping loop:
    # now starting at n=2 for this 2-step method!
    
    for n in range(2,nsteps+1):
        t[n] = t[n-1] + dt
        
        # Take a time step:
        J = array(range(1,m+1), dtype=int)
        
        # Leapfrog at interior points:
        u_computed[J,n] = u_computed[J,n-2] - nu*(u_computed[J+1,n-1] - u_computed[J-1,n-1])
        
        # inflow boundary:
        u_computed[0,n] = utrue(ax, t[n])  # exact inflow boundary value
        
        # outflow boundary:
        
        if advection_solution_input.bc == 'utrue':
            # specify exact solution at outflow:
            u_computed[m+1,n] = utrue(bx, t[n])  # exact inflow boundary value
        
        elif advection_solution_input.bc == 'zero':
            # zero at boundary:
            u_computed[m+1,n] = 0.
            
        elif advection_solution_input.bc == 'modLF':
            # modified Leapfrog:
            u_computed[m+1,n] = u_computed[m+1,n-2] - 2*nu*(u_computed[m+1,n-1] - u_computed[m,n-1])
            
        elif advection_solution_input.bc == 'upwind':
            # first order upwind:
            u_computed[m+1,n] = u_computed[m+1,n-1] - nu*(u_computed[m+1,n-1] - u_computed[m,n-1])
            
        elif advection_solution_input.bc == 'BW':
            # Beam-Warming:
            u_computed[m+1,n] = u_computed[m+1,n-1] - 0.5*nu*(3*u_computed[m+1,n-1] - 4*u_computed[m,n-1] \
                                                             + u_computed[m-1,n-1]) \
                                            + 0.5*nu**2 * (u_computed[m+1,n-1] - 2*u_computed[m,n-1] \
                                                             + u_computed[m-1,n-1])           
        else:
            raise IOError('*** unrecognized bc option')
            

        error[:,n] = u_computed[:,n] - utrue(x,t[n])
                
    advection_solution_output = AdvectionSolutionOutput()  # create object for output
    advection_solution_output.dt = dt
    advection_solution_output.h = h
    advection_solution_output.nu = nu
    advection_solution_output.t = t
    advection_solution_output.x_computed = x
    advection_solution_output.u_computed = u_computed
    advection_solution_output.error = error    
        
    return advection_solution_output

Smooth initial conditions:

In [6]:
def eta_gaussian(x):
    """Initial conditions"""
    beta = 600.
    return exp(-beta*(x - 0.5)**2)

ax = 0.
bx = 1.
a = 1. # advection velocity

def utrue_gaussian(x,t):
    """
    True solution for comparison.
    """
    xat = x - a*t
    return eta_gaussian(xat)

Using the true solution at the outflow boundary

Mathematically we cannot impose the solution at the outflow boudary, but in a case where we know the true solution you might think that this would be ok to use for the value at the right boundary, and might give very good results together with Leapfrog elsewhere.

But not so....

In [7]:
advection_solution_input = AdvectionSolutionInput()
advection_solution_input.t0 = 0.
advection_solution_input.tfinal = 2.
advection_solution_input.ax = ax
advection_solution_input.bx = bx
advection_solution_input.utrue = utrue_gaussian
advection_solution_input.a = a

advection_solution_input.mx = 99
advection_solution_input.nsteps = 250
advection_solution_input.bc = 'utrue'

advection_solution_output = Leapfrog_Advection(advection_solution_input)

error_tfinal = abs(advection_solution_output.error[:,-1]).max()
print('Using %i time steps' % advection_solution_input.nsteps)
print('Courant number nu = %.2f' % advection_solution_output.nu)
print('Max-norm Error at t = %6.4f is %12.8f' % (advection_solution_input.tfinal, error_tfinal))

Using 250 time steps
Courant number nu = 0.80
Max-norm Error at t = 2.0000 is   0.12795315
In [8]:
anim = make_animation(advection_solution_input, advection_solution_output, nplot=10)
HTML(anim.to_jshtml())  # or use the line below...
#HTML(anim.to_html5_video())
Out[8]:


Once Loop Reflect

Setting $u = 0$ at the outflow boundary

You shouldn't expect this to work at all well and it doesn't. The method is still stable for $|ak/h| \leq 1$, but no longer consistent with the PDE near the boundary and so even if you refine the grid this would not converge to the right thing. (Note that the reflected waves at the boundary are of magnitude 1.)

In [9]:
advection_solution_input.mx = 99
advection_solution_input.nsteps = 250
advection_solution_input.bc = 'zero'

advection_solution_output = Leapfrog_Advection(advection_solution_input)

error_tfinal = abs(advection_solution_output.error[:,-1]).max()
print('Using %i time steps' % advection_solution_input.nsteps)
print('Courant number nu = %.2f' % advection_solution_output.nu)
print('Max-norm Error at t = %6.4f is %12.8f' % (advection_solution_input.tfinal, error_tfinal))

anim = make_animation(advection_solution_input, advection_solution_output, nplot=10)
HTML(anim.to_jshtml()) 

Using 250 time steps
Courant number nu = 0.80
Max-norm Error at t = 2.0000 is   0.83247932
Out[9]:


Once Loop Reflect

In both of these experiments note that the left-going wave generated at the right boundary is highly oscillatory. When this wave hits the left boundary it interacts with the boundary condition imposed there to generate a right-going error, which is smoother. This relates to the group velocities related to different wave numbers.

Modified Leapfrog

In [10]:
advection_solution_input.mx = 99
advection_solution_input.nsteps = 250
advection_solution_input.bc = 'modLF'

advection_solution_output = Leapfrog_Advection(advection_solution_input)

error_tfinal = abs(advection_solution_output.error[:,-1]).max()
print('Using %i time steps' % advection_solution_input.nsteps)
print('Courant number nu = %.2f' % advection_solution_output.nu)
print('Max-norm Error at t = %6.4f is %12.8f' % (advection_solution_input.tfinal, error_tfinal))

anim = make_animation(advection_solution_input, advection_solution_output, nplot=10)
HTML(anim.to_jshtml()) 

Using 250 time steps
Courant number nu = 0.80
Max-norm Error at t = 2.0000 is   0.02386960
Out[10]:


Once Loop Reflect

First-order upwind

In [11]:
advection_solution_input.mx = 99
advection_solution_input.nsteps = 250
advection_solution_input.bc = 'upwind'

advection_solution_output = Leapfrog_Advection(advection_solution_input)

error_tfinal = abs(advection_solution_output.error[:,-1]).max()
print('Using %i time steps' % advection_solution_input.nsteps)
print('Courant number nu = %.2f' % advection_solution_output.nu)
print('Max-norm Error at t = %6.4f is %12.8f' % (advection_solution_input.tfinal, error_tfinal))

anim = make_animation(advection_solution_input, advection_solution_output, nplot=10)
HTML(anim.to_jshtml()) 

Using 250 time steps
Courant number nu = 0.80
Max-norm Error at t = 2.0000 is   0.00447809
Out[11]:


Once Loop Reflect

Beam-Warming at the outflow boundary

In [12]:
advection_solution_input.mx = 99
advection_solution_input.nsteps = 250
advection_solution_input.bc = 'BW'

advection_solution_output = Leapfrog_Advection(advection_solution_input)

error_tfinal = abs(advection_solution_output.error[:,-1]).max()
print('Using %i time steps' % advection_solution_input.nsteps)
print('Courant number nu = %.2f' % advection_solution_output.nu)
print('Max-norm Error at t = %6.4f is %12.8f' % (advection_solution_input.tfinal, error_tfinal))

anim = make_animation(advection_solution_input, advection_solution_output, nplot=10)
HTML(anim.to_jshtml()) 

Using 250 time steps
Courant number nu = 0.80
Max-norm Error at t = 2.0000 is   0.00198322
Out[12]:


Once Loop Reflect

Order of accuracy

Several of these methods do give convergent solutions and second order accuracy as the grid is refine. Even using the first-order upwind method at the boundary perserves second order accuracy:

In [13]:
advection_solution_input.bc = 'upwind'
advection_solution_input.tfinal = 0.8  # after Gaussian leaves boundary

r_vals = array([1,2,4,8,16,32], dtype=int)
m_vals = 50*r_vals - 1
nsteps_vals = 50*r_vals

E = empty(len(nsteps_vals))
h_vals = empty(len(nsteps_vals))

# print table header:
print("   h         dt      Courant #     error      ratio  estimated order")

for j,nsteps in enumerate(nsteps_vals):
    advection_solution_input.nsteps = nsteps
    advection_solution_input.mx = m_vals[j]  #int(0.5*nsteps) - 1
    advection_solution_output = Leapfrog_Advection(advection_solution_input)
    E[j] = abs(advection_solution_output.error[:,-1]).max()
    h_vals[j] = advection_solution_output.h
    dt = advection_solution_output.dt
    nu = advection_solution_output.nu
    
    if j>0:
        ratio = E[j-1] / E[j]
    else:
        ratio = nan
        
    p = log(ratio)/log(2)
    print("%8.6f  %8.6f  %8.4f  %12.8f    %4.2f        %4.2f" % (h_vals[j], dt, nu, E[j], ratio, p))

loglog(h_vals, E, '-o')
title('Log-log plot of errors')
xlabel('h = Delta x')
ylabel('error')

E_upwind = E  # save to plot again below

   h         dt      Courant #     error      ratio  estimated order
0.020000  0.016000    0.8000    0.09782093     nan         nan
0.010000  0.008000    0.8000    0.00605483    16.16        4.01
0.005000  0.004000    0.8000    0.00152274    3.98        1.99
0.002500  0.002000    0.8000    0.00037997    4.01        2.00
0.001250  0.001000    0.8000    0.00009409    4.04        2.01
0.000625  0.000500    0.8000    0.00002346    4.01        2.00

But even better results are obtained if the second-order Beam-Warming method is used at the boundary. In fact the residual error after the Gaussian hump has passed out of the domain (at tfinal = 0.8 below) decays to zero like $O(h^3)$ with this better boundary condition:

In [14]:
loglog(h_vals, E_upwind, '-o', label='upwind BC')

advection_solution_input.bc = 'BW'
advection_solution_input.tfinal = 0.8

r_vals = array([1,2,4,8,16,32], dtype=int)
m_vals = 50*r_vals - 1
nsteps_vals = 50*r_vals

E = empty(len(nsteps_vals))
h_vals = empty(len(nsteps_vals))

# print table header:
print("   h         dt      Courant #     error      ratio  estimated order")

for j,nsteps in enumerate(nsteps_vals):
    advection_solution_input.nsteps = nsteps
    advection_solution_input.mx = m_vals[j]  #int(0.5*nsteps) - 1
    advection_solution_output = Leapfrog_Advection(advection_solution_input)
    E[j] = abs(advection_solution_output.error[:,-1]).max()
    h_vals[j] = advection_solution_output.h
    dt = advection_solution_output.dt
    nu = advection_solution_output.nu
    
    if j>0:
        ratio = E[j-1] / E[j]
    else:
        ratio = nan
        
    p = log(ratio)/log(2)
    print("%8.6f  %8.6f  %8.4f  %12.8f    %4.2f        %4.2f" % (h_vals[j], dt, nu, E[j], ratio, p))

loglog(h_vals, E, '-o', label='Beam-Warming BC')
title('Log-log plot of errors')
xlabel('h = Delta x')
ylabel('error')
legend()

   h         dt      Courant #     error      ratio  estimated order
0.020000  0.016000    0.8000    0.09486074     nan         nan
0.010000  0.008000    0.8000    0.00264698    35.84        5.16
0.005000  0.004000    0.8000    0.00035276    7.50        2.91
0.002500  0.002000    0.8000    0.00003933    8.97        3.17
0.001250  0.001000    0.8000    0.00000473    8.32        3.06
0.000625  0.000500    0.8000    0.00000059    8.08        3.01
Out[14]:
<matplotlib.legend.Legend at 0x11a06f6d8>