Advection tests

This notebook presents an implementation of the first-order upwind method on the scalar advection equation $q_t + \bar u q_x = 0$ with periodic boundary conditions.

This notebook was developed to accompany AMath 574 Lab 2.

To load numpy and matplotlib and force plots to appear inline:



In [1]:

    
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib

Some Python tips

The extent of for loops, if-then-else blocks, function definitions is determined entirely by indentation.
Array indexes start a 0, so a list or array A with 3 elements would be indexed by A[0], A[1], A[2].
range(3) gives the list [0,1,2].
range(1,4) gives the list [1,2,3].
If A has at least 4 elements then A[1:4] is a new list [A[1], A[2], A[3]]
The linspace command works as in Matlab, and returns a numpy array (not a list)

Define a function implementing the upwind method

This function approximates the solution from time $t_0 = 0$ to some final time tfinal by taking nsteps time steps with the upwind method, applied to the advection equation $q_t + \bar u q_x = 0$.

On input, x is an array of cell centers and q0 should be an array of the same length, containing the initial cell averages at time 0.



In [2]:

    
def upwind(ubar,q0,x,tfinal,nsteps):
    dt = float(tfinal)/nsteps
    dx = x[1] - x[0]  # assume equally spaced
    udtdx = ubar * dt / dx
    cfl = abs(udtdx)
    print "dx = %g,  dt = %g" % (dx,dt)
    print "Courant number is ",cfl
    
    qn = hstack([0, q0, 0])   # add a ghost cell on each end
    mx = len(x)  # number of grid cells
    mx2 = mx+2   # number of cells with ghost cells
    for n in range(nsteps):
        qn[0] = qn[mx]
        qn[mx+1] = qn[1]
        qnp = zeros(qn.shape)   # initialize array
        for i in range(1,mx+1):
            if ubar > 0:
                qnp[i] = qn[i] - udtdx * (qn[i]-qn[i-1])
            else:
                qnp[i] = qn[i] - udtdx * (qn[i+1]-qn[i])
                
        qn = qnp   # for next time step
    qfinal = qnp[1:(mx+1)]   # throw away ghost cells
    return qfinal

Set up the computational grid and also a much finer grid for plotting the "exact" solution:



In [3]:

    
mx = 50
dx = 1./mx
x = linspace(dx/2, 1.-dx/2, mx)  # computational grid
xfine = linspace(0,1,5000)   # fine grid for plotting true solution

Define the true solution as a function of $(x,t)$. It will depend on ubar.

Note how the periodic boundary conditions are imposed.



In [4]:

    
def qtrue(x,t,ubar):
    x0 = x - ubar*t   # trace back characteristic to time 0
    x0 = mod(x0, 1.)  # use periodic boundary conditions to map to [0,1]
    q = where(abs(x0-0.2) < 0.1,  3., 2.)   # piecewise constant with values 2 and 3
    return q

Plot the initial data, with blue dots for the cell averages on the computational grid and a red line for the "exact" data.



In [5]:

    
ubar = 1.
q0 = qtrue(x,0.,ubar)
plot(x,q0, 'bo')

q0fine = qtrue(xfine,0.,ubar)
plot(xfine,q0fine,'r-')
ylim(1,4)









    Out[5]:





(1, 4)

Test the upwind method for specific values of the parameters:



In [6]:

    
tfinal = 0.5
nsteps = 30
q0 = qtrue(x,0.,ubar)
q = upwind(ubar,q0,x,tfinal,nsteps)
qfine = qtrue(xfine,tfinal,ubar)  # "exact" solution
plot(x,q, 'bo')
plot(xfine,qfine,'r-')
ylim(1,4)









    



dx = 0.02,  dt = 0.0166667
Courant number is  0.833333333333






    Out[6]:





(1, 4)

Try the following:

Check that this code works if ubar is negative and/or if the time interval is longer, so that the periodic boundary conditions play a role.
With tfinal = 0.5, try the following and make some observations about the results in each case:
- nsteps = 25 (why is the result so good?)
- nsteps = 24 (why is the result so bad?)
Write a new function lax_friedrichs that implements the Lax-Friedrichs method. Try various paramter values and comment on the results. $$Q_i^{n+1} = \frac 1 2 (Q_{i-1}^n + Q_{i+1}^n) - \frac{\bar u \Delta t}{2\Delta x}(Q_{i+1}^n-Q_{i-1}^n)$$
Try both upwind and Lax-Friedrichs on the initial data $q_0(x) = \exp(-30(x-0.2)^2)$ for $0 \leq x \leq 1$ and comment on what you observe.
Implement the second-order accurate Lax-Wendroff method and try it out on these examples.

$$Q_i^{n+1} = Q_i^n - \frac{\bar u \Delta t}{2\Delta x}(Q_{i+1}^n-Q_{i-1}^n) + \frac 1 2 \left(\frac{\bar u \Delta t}{\Delta x}\right)^2 (Q_{i+1}^n - 2Q_i^n + Q_{i-1}^n)$$

Sample solutions

Try 25 time steps. Note that the solution is exact because the upwind method reduces to $Q_i^{n+1} = Q_{i-1}^n$ in this case, so each cell average simply shifts over one cell...



In [7]:

    
tfinal = 0.5
nsteps = 25
q0 = qtrue(x,0.,ubar)
q = upwind(ubar,q0,x,tfinal,nsteps)
qfine = qtrue(xfine,tfinal,ubar)  # "exact" solution
plot(x,q, 'bo')
plot(xfine,qfine,'r-')
ylim(1,4)









    



dx = 0.02,  dt = 0.02
Courant number is  1.0






    Out[7]:





(1, 4)

Try 24 time steps. Note that in this case the Courant number is greater than 1, so the method is unstable. The square pulse should shift over 1.04167 grid cells each time step, but the numerical domain of dependence allows information to propagate at most one grid cell each time step. Notice that the numerical solution is still $Q_i^n=2$ to the right of $x = 0.3 + 24\Delta x = 0.78$.



In [8]:

    
tfinal = 0.5
nsteps = 24
q0 = qtrue(x,0.,ubar)
q = upwind(ubar,q0,x,tfinal,nsteps)
qfine = qtrue(xfine,tfinal,ubar)  # "exact" solution
plot(x,q, 'bo')
plot(xfine,qfine,'r-')
ylim(1,4)









    



dx = 0.02,  dt = 0.0208333
Courant number is  1.04166666667






    Out[8]:





(1, 4)

Lax-Friedrichs method



In [9]:

    
def lax_friedrichs(ubar,q0,x,tfinal,nsteps):
    dt = float(tfinal)/nsteps
    udtdx = ubar * dt / dx
    cfl = abs(udtdx)
    print "Courant number is ",cfl
    
    qn = hstack([0, q0, 0])   # add a ghost cell on each end
    mx = len(x)  # number of grid cells
    mx2 = mx+2   # number of cells with ghost cells
    for n in range(nsteps):
        qn[0] = qn[mx]
        qn[mx+1] = qn[1]
        qnp = zeros(qn.shape)   # initialize array
        for i in range(1,mx+1):
            qnp[i] = 0.5*(qn[i-1] + qn[i+1]) - 0.5*udtdx*(qn[i+1]-qn[i-1])
        qn = qnp   # for next time step
    qfinal = qnp[1:(mx+1)]   # throw away ghost cells
    return qfinal

Note that since $Q_i^{n+1}$ depends only on $Q_{i-1}^n$ and $Q_{i+1}^n$, the odd numbered grid points depend only on the even numbered points at the previous step, and vice versa. For the discontinous data, this results in the funny looking plots where each $Q$ value is repeated twice in the next test...



In [10]:

    
tfinal = 0.5
nsteps = 30
q0 = qtrue(x,0.,ubar)
q = lax_friedrichs(ubar,q0,x,tfinal,nsteps)
qfine = qtrue(xfine,tfinal,ubar)  # "exact" solution
plot(x,q, 'bo')
plot(xfine,qfine,'r-')
ylim(1,4)









    



Courant number is  0.833333333333






    Out[10]:





(1, 4)

Change the initial data to a Gaussian

Note that since we are solving the problem with periodic boundary conditions and the values $q(0,0) \neq q(1,0)$, there is effectively a discontinuity at the boundary that will propagate in.



In [11]:

    
def qtrue2(x,t,ubar):
    x0 = x - ubar*t   # trace back characteristic to time 0
    x0 = mod(x0, 1.)  # use periodic boundary conditions to map to [0,1]
    q = 2 + exp(-30*(x0-0.2)**2)
    return q
plot(xfine, qtrue2(xfine,0,ubar), 'r-')
ylim(1,4)









    Out[11]:





(1, 4)

Test both upwind and Lax-Friedrichs on this data. Go out to larger $t$ to better see the fact that Lax-Friedrichs is more diffusive.



In [12]:

    
tfinal = 2.5
nsteps = 150
q0 = qtrue2(x,0.,ubar)
q = lax_friedrichs(ubar,q0,x,tfinal,nsteps)
qfine = qtrue2(xfine,tfinal,ubar)  # "exact" solution
plot(x,q, 'bo', markersize=4, label='Lax-Friedrichs')
plot(xfine,qfine,'r-')
q = upwind(ubar,q0,x,tfinal,nsteps)
plot(x,q,'g+', label='upwind')
legend()
ylim(1,4)









    



Courant number is  0.833333333333
dx = 0.02,  dt = 0.0166667
Courant number is  0.833333333333






    Out[12]:





(1, 4)

Lax-Wendroff



In [13]:

    
def lax_wendroff(ubar,q0,x,tfinal,nsteps):
    dt = float(tfinal)/nsteps
    udtdx = ubar * dt / dx
    cfl = abs(udtdx)
    print "Courant number is ",cfl
    
    qn = hstack([0, q0, 0])   # add a ghost cell on each end
    mx = len(x)  # number of grid cells
    mx2 = mx+2   # number of cells with ghost cells
    for n in range(nsteps):
        qn[0] = qn[mx]
        qn[mx+1] = qn[1]
        qnp = zeros(qn.shape)   # initialize array
        for i in range(1,mx+1):
            qnp[i] = qn[i] - 0.5*udtdx*(qn[i+1]-qn[i-1]) + 0.5*udtdx**2 * (qn[i+1] - 2*qn[i] + qn[i-1])
        qn = qnp   # for next time step
    qfinal = qnp[1:(mx+1)]   # throw away ghost cells
    return qfinal

Test Lax-Wendroff out to time $t=2.5$. Note that it does a better job where the solution is smooth, but introduces oscillations trailing behind the discontinuity.



In [14]:

    
tfinal = 2.5
nsteps = 150
q0 = qtrue2(x,0.,ubar)
q = lax_wendroff(ubar,q0,x,tfinal,nsteps)
qfine = qtrue2(xfine,tfinal,ubar)  # "exact" solution
plot(x,q, 'bo')
plot(xfine,qfine,'r-')
ylim(1,4)









    



Courant number is  0.833333333333






    Out[14]:





(1, 4)

The oscillations are more apparent if we switch back to the step function initial data...



In [15]:

    
tfinal = 2.5
nsteps = 150
q0 = qtrue(x,0.,ubar)
q = lax_wendroff(ubar,q0,x,tfinal,nsteps)
qfine = qtrue(xfine,tfinal,ubar)  # "exact" solution
plot(x,q, 'bo')
plot(xfine,qfine,'r-')
ylim(1,4)









    



Courant number is  0.833333333333






    Out[15]:





(1, 4)

Note that the oscillations do not disappear if we reduce the grid size $\Delta x$ and $\Delta t$...



In [16]:

    
mx = 400
dx = 1./mx
x = linspace(dx/2, 1.-dx/2, mx)  # computational grid
tfinal = 2.5
nsteps = 1200
q0 = qtrue(x,0.,ubar)
q = lax_wendroff(ubar,q0,x,tfinal,nsteps)
qfine = qtrue(xfine,tfinal,ubar)  # "exact" solution
plot(x,q, 'bo', markersize=4)
plot(xfine,qfine,'r-')
ylim(1,4)









    



Courant number is  0.833333333333






    Out[16]:





(1, 4)



In [16]: