Solving the advection equation¶
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 periodic boundary conditions.
In [1]:
%pylab inline
In [2]:
from matplotlib import animation
from IPython.display import HTML
Function to make animations of solution:¶
In [3]:
def make_animation(adv_input, adv_output, nplot=1):
"""
Plot every `nplot` frames of the solution and turn into
an animation.
"""
xfine = linspace(adv_input.ax,adv_input.bx,1001)
fig, ax = plt.subplots(figsize=(12,6))
ax.set_xlim((adv_input.ax,adv_input.bx))
ax.set_ylim((-1.2, 1.2))
line1, = ax.plot([], [], '+-', color='b', lw=2, label='computed')
line2, = ax.plot([], [], color='r', lw=1, label='true')
ax.legend(loc='lower left')
title1 = ax.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)
return (line1,line2,title1)
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])
return (line1,line2,title1)
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
Define new classes for input and outputs¶
Note that the AdvectionSolutionInput class includes a time_stepper attribute that we will set to different functions below in order to test different methods.
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.time_stepper = None
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.errors = None
General explicit one-step method¶
In [5]:
def ExplicitAdvection(advection_solution_input):
"""
Solve u_t + a*u_x = 0 on [ax,bx] with periodic boundary conditions,
using centered differences in space and the an arbitrary 1-step method for time stepping,
defined by `advection_solution_input.time_stepper`.
with m+2 grid points (m+1 unknowns), taking nsteps time steps.
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.
This routine can be embedded in a loop on m to test the accuracy.
Note: the vector x and rows of `u_computed` are of length m+2 and includes both boundary points.
Only one is computed, and then the solution is extended to the other with the periodic BCs.
"""
# 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)
t = empty((nsteps+1,), dtype=float)
u_computed = empty((m+2,nsteps+1), dtype=float)
t[0] = t0
error = zeros((m+2,nsteps+1,), dtype=float)
u_computed[:,0] = u0
nu = a*dt/h # Courant number
# main time-stepping loop:
for n in range(1,nsteps+1):
t[n] = t[n-1] + dt
# Take a time step:
u_computed[:,n] = advection_solution_input.time_stepper(u_computed[:,n-1], nu, m)
# augment with boundary value:
u_computed[m+1,n] = u_computed[0,n]
# compute error at this time:
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
Define a smooth solution¶
For any given initial data $u(x,0) = \eta(x)$ in our interval $0 \leq x \leq 1$, the true solution of the advection equation $u_t + au_x=0$ with periodic boundary conditions is simply
$$
u(x,t) = \eta(\text{mod}(x-at, 1))
$$
where the mod function takes the fractional part and maps $x-at$ back to the interval $[0,1]$. This is generalized in the function below to an arbitrary interval [ax, bx].
Here we use the Gaussian $\eta(x)=\exp(-\beta(x - 0.5)^2)$, which for $\beta$ is sufficiently large, decays to zero sufficiently fast near $x=0$ and $x=1$ that the periodic extension is very smooth, so we use this as an initial test problem.
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.
For periodic BC's, we need the periodic extension of eta(x).
Map x-a*t-ax back to interval of length bx-ax
and then evaluate initial data at this point.
"""
xat = ax + mod(x - a*t - ax, bx-ax)
return eta_gaussian(xat)
Forward Euler time stepping¶
This function takes the computed solution u at one time and takes a single time step, returning the resulting array. It is assumed that u is an array of length m+2 that also contains both boundary points.
In [7]:
def FE_time_stepper(u, nu, m):
# check input:
assert len(u) == m+2, "Error: u has unexpected length relative to m"
# indices of interior points and one boundary point, values to update,
# as in integer numpy array:
J = array(range(0,m+1), dtype=int)
# Create indices J-1 and J+1 with periodic boundary conditions:
J = array(range(0,m+1), dtype=int)
Jm1 = mod(J-1, m+1)
Jp1 = mod(J+1, m+1)
u_next = empty(u.shape)
u_next[J] = u[J] - 0.5*nu*(u[Jp1] - u[Jm1])
return u_next
Recall that this method is Lax-Richtmyer unstable if we fix the Courant number as $k,h \rightarrow 0$, although it is stable if we set $k = Ch^2$ as $h\rightarrow 0$. This suggests we need to take $k$ small relative to $h$ to get reasonable results, but even then the solution grows slowly with time.
In [8]:
advection_solution_input = AdvectionSolutionInput()
advection_solution_input.t0 = 0.
advection_solution_input.tfinal = 1.
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 = 1000
advection_solution_input.time_stepper = FE_time_stepper
advection_solution_output = ExplicitAdvection(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))
In [9]:
anim = make_animation(advection_solution_input, advection_solution_output, nplot=100)
HTML(anim.to_jshtml()) # or use the line below...
#HTML(anim.to_html5_video())
Out[9]:
In addition to growing with time, the numerical solution also becomes oscillatory because this method is very dispersive.
Lax-Friedrichs¶
The Lax-Friedrichs method is stable provided the Courant number satisfies $|ak/h| \leq 1$. However, it is still only first order accurate.
In [10]:
def LxF_time_stepper(u, nu, m):
# check input:
assert len(u) == m+2, "Error: u has unexpected length relative to m"
# indices of interior points and one boundary point,
# as in integer numpy array:
J = array(range(0,m+1), dtype=int)
# Create indices J-1 and J+1 with periodic boundary conditions:
J = array(range(0,m+1), dtype=int)
Jm1 = mod(J-1, m+1)
Jp1 = mod(J+1, m+1)
u_next = empty(u.shape)
u_next[J] = 0.5*(u[Jm1] + u[Jp1]) - 0.5*nu*(u[Jp1] - u[Jm1])
return u_next
In [11]:
advection_solution_input = AdvectionSolutionInput()
advection_solution_input.t0 = 0.
advection_solution_input.tfinal = 1.
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 = 120
advection_solution_input.time_stepper = LxF_time_stepper
advection_solution_output = ExplicitAdvection(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))
In [12]:
# make an animation of the results, plotting selective frames:
anim = make_animation(advection_solution_input, advection_solution_output, nplot=10)
HTML(anim.to_jshtml())
Out[12]:
Note that this behaves very differently from Forward Euler, and gives a very diffusive or dissipative solution.
Lax-Wendroff¶
This method is second-order accurate in both space and time.
In [13]:
def LW_time_stepper(u, nu, m):
# check input:
assert len(u) == m+2, "Error: u has unexpected length relative to m"
# indices of interior points and one boundary point,
# as in integer numpy array:
J = array(range(0,m+1), dtype=int)
# Create indices J-1 and J+1 with periodic boundary conditions:
J = array(range(0,m+1), dtype=int)
Jm1 = mod(J-1, m+1)
Jp1 = mod(J+1, m+1)
u_next = empty(u.shape)
u_next[J] = u[J] - 0.5*nu*(u[Jp1] - u[Jm1]) \
+ 0.5*nu**2 * (u[Jm1] - 2*u[J] + u[Jp1])
return u_next
In [14]:
advection_solution_input = AdvectionSolutionInput()
advection_solution_input.t0 = 0.
advection_solution_input.tfinal = 1.
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 = 200
advection_solution_input.time_stepper = LW_time_stepper
advection_solution_output = ExplicitAdvection(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))
In [15]:
anim = make_animation(advection_solution_input, advection_solution_output, nplot=20)
HTML(anim.to_jshtml())
Out[15]:
This method is dispersive, but much less so than Forward Euler. It is stable for $|ak/h| \leq 1$ and second-order accurate in both space and time.
Test the order of accuracy¶
We can test for second order accuracy by fixing the courant number and running the code for a sequence of grid resolutions:
In [16]:
r_vals = array([1,2,4,8,16,32], dtype=int)
m_vals = 50*r_vals - 1
nsteps_vals = 75*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]
advection_solution_output = ExplicitAdvection(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')
Out[16]:
Note that second-order accuracy is observed once $h$ is sufficiently small, although for larger $h$ higher order terms in the error expansion dominate.
Discontinuous data¶
In [17]:
def eta_box(x):
"""Initial conditions"""
return where(logical_or(x<0.3, x>0.7), -0.5, 0.5)
ax = 0.
bx = 1.
a = 1. # advection velocity
def utrue_box(x,t):
"""
True solution for comparison.
For periodic BC's, we need the periodic extension of eta(x).
Map x-a*t-ax back to interval of length bx-ax
and then evaluate initial data at this point.
"""
xat = ax + mod(x - a*t - ax, bx-ax)
return eta_box(xat)
In [18]:
advection_solution_input = AdvectionSolutionInput()
advection_solution_input.t0 = 0.
advection_solution_input.tfinal = 1.
advection_solution_input.ax = ax
advection_solution_input.bx = bx
advection_solution_input.utrue = utrue_box
advection_solution_input.a = a
advection_solution_input.mx = 99
advection_solution_input.nsteps = 200
advection_solution_input.time_stepper = LW_time_stepper
advection_solution_output = ExplicitAdvection(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))
In [19]:
anim = make_animation(advection_solution_input, advection_solution_output, nplot=20)
HTML(anim.to_jshtml())
Out[19]:
Oscillations arise near the discontinuities, trailing behind.
Note that if we refine the grid, the oscillations do not disappear. They become narrower and more concentrated near the discontinuities, but the amplitude does not decrease:
In [20]:
advection_solution_input.mx = 499
advection_solution_input.nsteps = 1000
advection_solution_output = ExplicitAdvection(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))
In [21]:
anim = make_animation(advection_solution_input, advection_solution_output, nplot=100)
HTML(anim.to_jshtml())
Out[21]:
Wave packet data¶
Finally we show wave packet initial data, as discussed briefly in Appendix E.3.10 and for Lax-Wendroff in Example 10.10. We take
$$ \eta(x) = e^{-\beta (x-0.5)^2} \cos(\xi_0 x) $$
In [22]:
def eta_wavepacket(x):
"""Initial conditions"""
beta = 300.
xi0 = 2*pi*20
return exp(-beta*(x - 0.5)**2) * cos(xi0*x)
ax = 0.
bx = 1.
a = 1. # advection velocity
def utrue_wavepacket(x,t):
"""
True solution for comparison.
For periodic BC's, we need the periodic extension of eta(x).
Map x-a*t-ax back to interval of length bx-ax
and then evaluate initial data at this point.
"""
xat = ax + mod(x - a*t - ax, bx-ax)
return eta_wavepacket(xat)
We solve over a longer time period to see how the wave packet propagates too slowly:
In [23]:
advection_solution_input = AdvectionSolutionInput()
advection_solution_input.t0 = 0.
advection_solution_input.tfinal = 4.
advection_solution_input.ax = ax
advection_solution_input.bx = bx
advection_solution_input.utrue = utrue_wavepacket
advection_solution_input.a = a
advection_solution_input.mx = 299
advection_solution_input.nsteps = 2000
advection_solution_input.time_stepper = LW_time_stepper
advection_solution_output = ExplicitAdvection(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))
In [24]:
anim = make_animation(advection_solution_input, advection_solution_output, nplot=100)
HTML(anim.to_jshtml())
Out[24]:
We can compute the group velocity predicted for this wave packet with the Lax-Wendroff method (see Example 10.10) and we find that:
In [25]:
a = 1.
nu = 0.6
h = 1/300.
xi0 = 2*pi*20
cg = a - 0.5*a*h**2 * (1-nu**2) * xi0**2
print('The group velocity is %.2f' % cg)
print('At tfinal=4, the wave packet has only traveled distance %.2f' % (4*cg))