Upwinding and stability - Fundamentals of Numerical Computation

Advection phenomena are characterized by the transport of information at finite speed. This implies that the solution at a given point in space and time can be affected only by the initial data in a limited region of space. For now, we continue to ignore the effects of boundaries.

Any numerical method we choose to solve a PDE has analogous property.

Example 12.2.3

In $u_t+cu_x=0$ , 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.1)

If we use the Euler time discretization with step size τ, then

U_{i,j+1} = U_{i,j} - \frac{c\tau}{2h} (U_{i+1,j}-U_{i-1,j}).

(12.2.2)

Starting with $j=0$ , we find that $U_{i,1}$ depends on $U_{i-1,0}$ , $U_{i,0}$ , and $U_{i+1,0}$ . Hence, the numerical domain of dependence at $(x_i,t_1)$ is $\{x_{i-1},x_i,x_{i+1}\}$ .

Now we set $j=1$ . From (12.2.2), we see that $U_{i,2}$ depends on $U_{i-1,1}$ , $U_{i,1}$ , and $U_{i+1,1}$ . In turn, each of these values at time $t_1$ depends on three values at time $t_0$ :

\begin{align*} U_{i-1,1} &\;\text{ uses }\; U_{i-2,0},\, U_{i-1,0},\, U_{i,0}, \\ U_{i,1} &\;\text{ uses }\; U_{i-1,0},\, U_{i,0},\, U_{i+1,0}, \\ U_{i+1,1} &\;\text{ uses }\; U_{i,0},\, U_{i+1,0},\, U_{i+2,0}. \end{align*}

(12.2.3)

Considering only the unique values, the numerical domain of dependence at $(x_i,t_2)$ is $\{x_{i-2},x_{i-1},x_i,x_{i+1},x_{i+2}\}$ . Continuing this kind of reasoning backward from any $U_{i,j}$ , we find that the numerical domain of dependence at $(x_i,t_j)$ is $\{x_{i-j},\ldots,x_{i+j}\}$ .

12.2.1The CFL condition¶

We now state an important principle about a necessary relationship between 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.

Example 12.2.4

In Example 12.2.3 we concluded that the numerical domain of dependence for a centered Euler discretization of the advection equation at $(x_i,t_j)$ is $\{x_{i-j},\ldots,x_{i+j}\}$ . Figure 12.2.1 illustrates what happens as $h$ and τ go to zero in a manner that leaves the ratio $h/\tau$ constant.

$Numerical domain of dependence for the explicit time stepping scheme in Example {number}. As \tau and h approach zero, the shaded region is filled in.$

Figure 12.2.1:Numerical domain of dependence for the explicit time stepping scheme in Example 12.2.3. As τ and $h$ approach zero, the shaded region is filled in.

In order to observe both the exact and numerical domains of dependence at a fixed location $(x,t)$ , we must have $x=ih$ and $t=j\tau$ . The numerical domain of dependence fills in the interval $[x_{i-j},x_{i+j}]$ , which is $[x-jh,x+jh]$ , or

[x - \frac{h}{\tau} t, x + \frac{h}{\tau} t].

(12.2.4)

This interval captures the exact domain of dependence $\{x-ct\}$ only if

|c| \le \frac{h}{\tau}.

(12.2.5)

Equation (12.2.5) 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, which leads to the following restatement.

We can rearrange (12.2.5) 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. This is our first clear indication that advection is less stiff than diffusion.

Example 12.2.5 (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.5

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());

Source

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))")

Number of time steps for m = 400: 565
Number of time steps for m = 800: 1128

Example 12.2.5

[x, Dx] = diffper(400, [0, 1]);
c = 2;
odefun = @(t, u) -c * (Dx*u);
u_init = exp( -80*(x - 0.5).^2 );
sol = ode45(odefun, [0, 2], u_init);
u = @(t) deval(sol, t);

Source

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]);
odefun = @(t, u) -c * (Dx*u);
u_init = exp( -80*(x - 0.5).^2 );
sol = ode45(odefun, [0, 2], u_init);

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

Unrecognized field name "Time".

Example 12.2.5

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

Source

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")

1416 time steps taken for m = 400

2906 time steps taken for m = 800

12.2.2Upwinding¶

The advection equation allows information to propagate in only one direction. This asymmetry can have important repercussions.

Numerical methods can also have a directional preference.

Example 12.2.6

Suppose that in $u_t+cu_x=0$ we use the backward difference

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

(12.2.6)

together with an Euler scheme in time. It should be clear that $U_{i,j}$ depends only on points to the left of $x_i$ , i.e., the upwind direction of the numerical method is to the left. But if $c<0$ , the upwind direction of the PDE is to the right. Hence, it is impossible to satisfy the CFL condition under these choices.

Similarly, the forward difference

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

(12.2.7)

with explicit time stepping leads to the inverse conclusion: its upwind direction is to the right, and it must fail if $c>0$ .

The reasoning of Example 12.2.6 is readily generalized: if the numerical method has an upwind direction, the CFL condition requires that it must agree with the upwind direction of the PDE.

12.2.3Inflow boundary condition¶

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

Upwinding considerations provide the answer. 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. Conversely, at points close to the upwind boundary, the domain of dependence soon wants to move past the boundary, which is impossible, so there has to be information provided there.

For $c>0$ in the advection equation, the inflow is at the left end, and for $c<0$ , it is at the right end.

Example 12.2.7 (Upwind versus downwind)

Julia

MATLAB

Python

Example 12.2.7

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.

Source

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/advection-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);

Source

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.

Source

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/advection-outflow.mp4")

Example 12.2.7

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) chop( -c * Dx * extend(v) );

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

u_init = exp( -80*(x - 0.5).^2 );
sol = ode45(odefun, [0, 1], chop(u_init));
u = @(t) [deval(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.

Source

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)

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];
odefun = @(t, v) chop( -c * Dx * extend(v) );
sol = ode45(odefun, [0, 1], chop(u_init));
u = @(t) [zeros(1, length(t)); deval(sol, t)];    % extend to zero at left

Source

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.

Source

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)

Example 12.2.7

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.

Source

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))

Source

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.

Source

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¶