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.

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



In [ ]:

    
%pylab inline

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 [ ]:

    
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 [ ]:

    
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 [ ]:

    
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 [ ]:

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

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

Test the upwind method for specific values of the parameters:



In [ ]:

    
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)

Try the following:

Some modifications and extensions will explored during the lab...



In [ ]: