Upwinding and stability - Fundamentals of Numerical Computation

Let’s focus on the constant-velocity linear advection equation,

u_t + c u_x = 0, \quad u(0,x)=u_0(x).

(12.2.1)

For now, we suppose there are no boundaries. Keep in mind that $c$ is a velocity, not a speed: if $c>0$ , solutions travel rightward, and if $c<0$ , they travel leftward.

In Traffic flow we argued that $u(x,t)=\psi(x-ct)$ is a solution of (12.2.1). It’s therefore clear that $u(x,t)=u_0(x-ct)$ .

In the advection equation, the domain of dependence at $(x,t)$ is the single point $\{x-ct\}$ , and the upwind direction is to the left or to the right of $x$ if $c$ is positive or negative, respectively.

Any numerical method we choose has an analogous property.

Example 12.2.1

In (12.2.1), suppose we discretize $u_x$ by a centered difference,

u_x(x_i,t_j) \approx \frac{U_{i+1,j}-U_{i-1,j}}{2h}.

(12.2.2)

If we use the Euler time discretization, then

\mathbf{u}_{j+1} = (\mathbf{I} - c \tau \mathbf{D}_x) \mathbf{u}_j,

(12.2.3)

where τ is the time step. Because the matrix in this time step is tridiagonal, the entry $U_{i,j}$ can depend directly only on $U_{i-1,j}$ , $U_{i,j}$ , and $U_{i+1,j}$ . Going back another time step, the dependence extends to space positions $i-2$ and $i+2$ , and so on. When we reach the initial time, the dependence of $U_{i,j}$ reaches from $x_{i-j}$ to $x_{i+j}$ , or between $x_i-jh$ and $x_i+jh$ . If we ignore boundaries, the situation is illustrated in Figure 12.2.1. As $\tau,h\rightarrow 0$ , the numerical domain of dependence fills in the shaded region in the figure, but that region itself does not change.

Numerical domain of dependence for the explicit time stepping scheme in Example {number}. If \tau and h are infinitesimally small, the shaded region is filled in. — Figure 12.2.1:Numerical domain of dependence for the explicit time stepping scheme in Example 12.2.1. If τ and $h$ are infinitesimally small, the shaded region is filled in.

12.2.1The CFL condition¶

We now state an important principle about a required relationship between the domains of dependence.

Although we will not provide the rigor behind this theorem, its conclusion is not difficult to justify. If the CFL condition does not hold, the exact solution at $(x,t)$ could be affected by a change in the initial data while having no effect on the numerical solution. Hence there is no way for the method to get the solution correct for all problems. By contradiction, then, the CFL criterion is necessary for convergence.

Equation (12.2.4) is the implication of the CFL condition for the stated discretization. Notice that $h/\tau$ is the speed at which information moves in the numerical method; thus, it is common to restate the CFL condition in terms of speeds.

We can rearrange (12.2.4) to imply a necessary time step restriction $\tau \le h/|c|$ . This restriction for advection is much less severe than the $\tau = O(h^2)$ restriction we derived for Euler in the heat equation in Stiffness, which is our first indication that advection is less stiff than diffusion.

Example 12.2.3 (The CFL condition in action)

We solve linear advection with velocity $c=2$ and periodic end conditions. The initial condition is numerically, though not mathematically, periodic.

Julia

MATLAB

Python

Example 12.2.3

For time stepping, we use the adaptive explicit method RK4.

using OrdinaryDiffEq
m = 400;
x, Dₓ = FNC.diffper(m, [0, 1])
u_init = x -> exp( -80 * (x - 0.5)^2 )
ode = (u, c, t) -> -c * (Dₓ*u)
IVP = ODEProblem(ode, u_init.(x), (0., 2.), 2.)
u = solve(IVP, RK4());

using Plots
t = range(0, 2, 81);
U = reduce(hcat, u(t) for t in t)
contour(x, t, U'; 
    color=:redsblues,  clims=(-1, 1),
    xaxis=(L"x"),  yaxis=(L"t"),
    title="Linear advection",  right_margin=3Plots.mm)

In the space-time plot above, you can see the initial hump traveling rightward at constant speed. It fully traverses the domain once for each integer multiple of $t=1/2$ .

If we cut $h$ by a factor of 2 (i.e., double $m$ ), then the CFL condition suggests that the time step should be cut by a factor of 2 also.

println("Number of time steps for m = 400: $(length(u.t))")

m = 800;
x, Dₓ = FNC.diffper(m, [0, 1])
IVP = ODEProblem(ode, u_init.(x), (0., 2.), 2.)
u = solve(IVP, RK4())
println("Number of time steps for m = 800: $(length(u.t))")

Example 12.2.3

[x, Dx] = diffper(400, [0, 1]);
c = 2;  
ivp = ode(ODEFcn = @(t, u) -c * (Dx*u));
ivp.RelativeTolerance = 1e-5;
ivp.InitialTime = 0;
u_init = exp( -80*(x - 0.5).^2 );
ivp.InitialValue = u_init;

[u, sol] = solutionFcn(ivp, 0, 2);

clf
t = linspace(0, 2, 81);
contour(x, t, u(t)', 12)
xlabel('x'),  ylabel('t')
title('Linear advection')

In the space-time plot above, you can see the initial hump traveling rightward at constant speed. It fully traverses the domain once for each integer multiple of $t=1/2$ .

If we cut $h$ by a factor of 2 (i.e., double $m$ ), then the CFL condition suggests that the time step should be cut by a factor of 2 also.

num_steps_400 = length(sol.Time) - 1

[x, Dx] = diffper(800, [0, 1]);
ivp.ODEFcn = @(t, u) -c * (Dx*u);
ivp.InitialValue = exp( -80*(x - 0.5).^2 );
[u, sol] = solutionFcn(ivp, 0, 2);

num_steps_800 = length(sol.Time) - 1
ratio = num_steps_800 / num_steps_400

Example 12.2.3

For time stepping, we use the adaptive explicit method RK45.

x, Dx, Dxx = FNC.diffper(400, [0, 1])
u_init = exp(-80 * (x - 0.5) ** 2)
c = 2
ode = lambda t, u: -c * (Dx @ u)
sol = solve_ivp(ode, (0, 2), u_init, method="RK45", dense_output=True)
u = sol.sol

t = linspace(0, 2, 81)
U = vstack([u(tj) for tj in t])
contour(x, t, U, levels=arange(0.15, 1.0, 0.2))
xlabel("$x$"),  ylabel("$t$")
title("Linear advection");

In the space-time plot above, you can see the initial hump traveling rightward at constant speed. It fully traverses the domain once for each integer multiple of $t=1/2$ .

If we cut $h$ by a factor of 2 (i.e., double $m$ ), then the CFL condition suggests that the time step should be cut by a factor of 2 also.

print(f"{len(sol.t) - 1} time steps taken for m = 400")

x, Dx, Dxx = FNC.diffper(800, [0, 1])
u_init = exp(-80 * (x - 0.5) ** 2)
sol = solve_ivp(ode, (0, 2), u_init, method="RK45", dense_output=True)
print(f"{len(sol.t) - 1} time steps taken for m = 800")

Example 12.2.4 is a special case of a more general conclusion.

12.2.2Upwinding¶

There are other ways to discretize the $u_x$ term in the advection equation (12.2.1). The implications of the CFL criterion may differ greatly depending on which is chosen.

It’s clear that when the domain of dependence and the numerical method both have a directional preference, they must agree.

It might seem like one should always use a centered difference scheme so that upwinding is not an issue. However, at a shock front, this requires differencing across a jump in the solution, which causes its own difficulties.

12.2.3Inflow boundary condition¶

Now suppose that (12.2.1) is posed on the finite domain $x \in [a,b]$ . Since the PDE has only a first-order derivative in $x$ , we should have only one boundary condition. Should it be specified at the left end, or the right end?

If we impose a condition at the downwind side of the domain, there is no way for that boundary information to propagate into the interior of the domain as time advances. On the other hand, for points close to the upwind boundary, the domain of dependence eventually wants to move past the left boundary. This is impossible, so instead the domain of dependence has to stay there.

In summary, we require an inflow condition on the PDE. For $c>0$ this is at the left end, and for $c<0$ it is at the right end. This requirement is true of the exact PDE as well as any discretization of it.

Example 12.2.6 (Upwind versus downwind)

Julia

MATLAB

Python

Example 12.2.6

If we solve advection over $[0,1]$ with velocity $c=-1$ , the right boundary is in the upwind/inflow direction. Thus a well-posed boundary condition is $u(1,t)=0$ .

We’ll pattern a solution after Function 11.5.2. Since $u(x_m,t)=0$ , we define the ODE interior problem (11.5.4) for $\mathbf{v}$ without $u_m$ . For each evaluation of $\mathbf{v}'$ , we must extend the data back to $x_m$ first.

m = 100
x, Dₓ = FNC.diffmat2(m, [0, 1])

interior = 1:m
extend = v -> [v; 0]

function ode!(f, v, c, t)
    u = extend(v)
    uₓ = Dₓ * u
    @. f = -c * uₓ[interior]
end;

Now we solve for an initial condition that has a single hump.

init = @. exp( -80*(x[interior] - 0.5)^2 )
ivp = ODEProblem(ode!, init, (0., 1), -1)
u = solve(ivp);

t = range(0, 0.75, 80)
U = reduce(hcat, extend(u(t)) for t in t)
contour(x, t, U';
    color=:blues,  clims=(0, 1), 
    xaxis=(L"x"),  yaxis=(L"t"),
    title="Advection with inflow BC")

We find that the hump gracefully exits out the downwind end.

anim = @animate for t in range(0, 1, 161) 
    plot(x, extend(u(t));
        label=@sprintf("t = %.4f", t), 
        xaxis=(L"x"),  yaxis=(L"u(x, t)", (0, 1)), 
        title="Advection equation with inflow BC",  dpi=150)
end
mp4(anim,"figures/upwind-inflow.mp4")

If instead of $u(1,t)=0$ we were to try to impose the downwind condition $u(0,t)=0$ , we only need to change the index of the interior nodes and where to append the zero value.

interior = 2:m+1
extend = v -> [0; v]

init = @. exp( -80*(x[interior] - 0.5)^2 )
ivp = ODEProblem(ode!, init, (0., 0.25), -1)
u = solve(ivp);

t = range(0, 0.2, 61)
U = reduce(hcat, extend(u(t)) for t in t)
contour(x, t, U'; 
    color=:redsblues,  clims=(-1, 1),
    xaxis=(L"x"),  yaxis=(L"t"), 
    title="Advection with outflow BC",  right_margin=3Plots.mm)

This time, the solution blows up as soon as the hump runs into the boundary because there are conflicting demands there.

anim = @animate for t in range(0, 0.2, 41) 
    plot(x, extend(u(t));
        label=@sprintf("t = %.4f", t),
        xaxis=(L"x"),  yaxis=(L"u(x,t)", (0, 1)), 
        title="Advection equation with outflow BC",  dpi=150)
end
mp4(anim,"figures/upwind-outflow.mp4")

Example 12.2.6

If we solve advection over $[0,1]$ with velocity $c=-1$ , the right boundary is in the upwind/inflow direction. Thus a well-posed boundary condition is $u(1,t)=0$ .

m = 100;  c = -1;
[x, Dx] = diffmat2(m, [0, 1]);
chop = @(u) u(1:m);  
extend = @(v) [v; 0];
odefun = @(t, v) -c * chop( Dx * extend(v) );
ivp = ode(ODEFcn = odefun);
ivp.RelativeTolerance = 1e-5;
ivp.InitialTime = 0;

Now we solve for an initial condition that has a single hump.

u_init = exp( -80*(x - 0.5).^2 );
ivp.InitialValue = chop(u_init);
sol = solutionFcn(ivp, 0, 1);
u = @(t) [sol(t); zeros(1, length(t))];    % extend to zero at right

t = linspace(0, 1, 81);
clf,  contour(x, t, u(t)', 0.15:0.2:1)
xlabel x,  ylabel t
title('Advection with inflow BC')

We find that the hump gracefully exits out the downwind end.

clf
plot(x, u(0))
hold on
axis([0, 1, -0.05, 1.05])
title('Advection with inflow BC') 
xlabel('x'),  ylabel('u(x,t)')
vid = VideoWriter("figures/advection-inflow.mp4", "MPEG-4");
vid.Quality = 85;
open(vid);
for frame = 1:length(t)
    cla, plot(x, u(t(frame)))
    str = sprintf("t = %.2f", t(frame));
    text(0.08, 0.85, str);
    writeVideo(vid, frame2im(getframe(gcf)));
end
close(vid)
close(gcf)

If instead of $u(1,t)=0$ we were to try to impose the downwind condition $u(0,t)=0$ , we only need to change the index of the interior nodes and where to append the zero value.

chop = @(u) u(2:m+1);  
extend = @(v) [0; v];
ivp.ODEFcn = @(t, v) -c * chop( Dx * extend(v) );
ivp.InitialValue = chop(u_init);
sol = solutionFcn(ivp, 0, 1);
u = @(t) [zeros(1, length(t)); sol(t)];

FNC_init
contour(x, t, u(t)', 0.15:0.2:1)
xlabel x,  ylabel t
title('Advection with outflow BC')

This time, the solution blows up as soon as the hump runs into the boundary because there are conflicting demands there.

clf
plot(x, u(0))
hold on
axis([0, 1, -0.05, 1.05])
title('Advection with outflow BC') 
xlabel('x'),  ylabel('u(x,t)')
vid = VideoWriter("figures/advection-outflow.mp4", "MPEG-4");
vid.Quality = 85;
open(vid);
for frame = 1:45
    cla, plot(x, u(t(frame)))
    str = sprintf("t = %.2f", t(frame));
    text(0.08, 0.85, str);
    writeVideo(vid, frame2im(getframe(gcf)));
end
close(vid)
close(gcf)

Example 12.2.6

If we solve advection over $[0,1]$ with velocity $c=-1$ , the right boundary is in the upwind/inflow direction. Thus a well-posed boundary condition is $u(1,t)=0$ .

m = 80
x, Dx, Dxx = FNC.diffmat2(m, [0, 1])

chop = lambda u : u[:-1]
extend = lambda v: hstack([v, 0])

ode = lambda t, v: -c * chop( Dx @ extend(v) )
c = -1

Now we solve for an initial condition that has a single hump.

u_init = exp(-80 * (x - 0.5) ** 2)
sol = solve_ivp(ode, (0, 1), chop(u_init), method="RK45", dense_output=True)
u = lambda t: extend(sol.sol(t))

t = linspace(0, 1, 80)
U = [u(tj) for tj in t]
contour(x, t, U, levels=arange(0.15, 1.0, 0.2))
xlabel("$x$"),  ylabel("$t$")
title("Advection with inflow BC");

We find that the hump gracefully exits out the downwind end.

from matplotlib import animation
fig, ax = subplots()
curve = ax.plot(x, u_init)[0]
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
ax.set_xlabel("$x$")
ax.set_ylabel("$u(x,t)$")
ax.set_ylim(-0.1, 1.1)
ax.set_title("Advection with inflow BC")

def snapshot(t):
    curve.set_ydata(u(t))
    time_text.set_text(f"t = {t:.2f}")

anim = animation.FuncAnimation(
    fig, snapshot, frames=linspace(0, 1, 101)
    )
anim.save("figures/advection-inflow.mp4", fps=30)
close()

If, instead of $u(1,t)=0$ , we were to try to impose the downwind condition $u(0,t)=0$ , we only need to change the index of the interior nodes and where to append the zero value.

chop = lambda u : u[1:]
extend = lambda v: hstack([0, v])

sol = solve_ivp(ode, (0, 1), chop(u_init), method="RK45", dense_output=True)
u = lambda t: extend(sol.sol(t))

U = [u(tj) for tj in t]
clf
contour(x, t, U, levels=arange(0.15, 1.0, 0.2))
xlabel("$x$"),  ylabel("$t$")
title("Outflow boundary condition");

This time, the solution blows up as soon as the hump runs into the boundary because there are conflicting demands there.

fig, ax = subplots()
curve = ax.plot(x, u_init)[0]
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
ax.set_xlabel("$x$")
ax.set_ylabel("$u(x,t)$")
ax.set_ylim(-0.1, 1.1)
ax.set_title("Advection with outflow BC")

anim = animation.FuncAnimation(
    fig, snapshot, frames=linspace(0, 0.5, 51)
    )
anim.save("figures/advection-outflow.mp4", fps=30)
close()

12.2.4Exercises¶

✍ Suppose you want to model the weather, including winds up to speed 200 km/hr, using an explicit method with a second-order centered spatial discretization. If the shortest time step you can take is 4 hr, what is the CFL limit on the spatial resolution of the model? Is this a lower bound or an upper bound?
✍ Suppose you want to model the traffic on a high-speed freeway using an explicit method with a second-order centered spatial discretization. Derive a CFL condition on the allowable time step, stating your assumptions carefully.
✍ For the heat equation, the domain of dependence at any $(x,t)$ with $t>0$ is all of $x \in (-\infty,\infty)$ . Show that the CFL condition implies that $\tau/h\to 0$ is required for convergence as $h\to 0$ .
✍ Suppose you wish to solve $u_t = u u_x$ for $x\in[-1,1]$ .
(a) If $u(x,0) = -2+\sin(\pi x)$ , which end of the domain is the inflow?
(b) Does the answer to part (a) change if $u(x,0) = 1 + e^{-16x^2}$ ?