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:
%pylab inline
Populating the interactive namespace from numpy and matplotlib
See also http://faculty.washington.edu/rjl/classes/am574w2015/python.html
A[0], A[1], A[2]
.range(3)
gives the list [0,1,2].range(1,4
) gives the list [1,2,3].A[1:4]
is a new list [A[1], A[2], A[3]]
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.
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:
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.
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.
ubar = 1.
q0 = qtrue(x,0.,ubar)
plot(x,q0, 'bo')
q0fine = qtrue(xfine,0.,ubar)
plot(xfine,q0fine,'r-')
ylim(1,4)
(1, 4)
Test the upwind method for specific values of the parameters:
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
(1, 4)
ubar
is negative and/or if the time interval is longer, so that the periodic boundary conditions play a role.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?)lax_friedrichs
that implements the Lax-Friedrichs method. Try various paramter values and comment on the results.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...
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
(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$.
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
(1, 4)
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...
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
(1, 4)
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.
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)
(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.
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
(1, 4)
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.
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
(1, 4)
The oscillations are more apparent if we switch back to the step function initial data...
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
(1, 4)
Note that the oscillations do not disappear if we reduce the grid size $\Delta x$ and $\Delta t$...
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
(1, 4)