Advecting Without Splitting Directions — The Bell-Colella-Glaz Unsplit Upwind Scheme

In the late 1980s, John Bell and Phillip Colella at Lawrence Berkeley spotted an odd scar in their explosion runs. A flame front set at an angle to the grid axes kept smearing and twisting along the diagonal as time went on. The culprit was not physics but discretization. Dimensional splitting — solving one direction at a time — was distorting diagonal advection. The fix is the Bell-Colella-Glaz (BCG) unsplit second-order upwind scheme. Split the directions and a rotating disk turns into a rhombus. One transverse derivative closes that hole.

Push CFL from 0 to 0.9 and the amber dot slides back toward the cell average. Set the limiter to none and the line beside the jump runs past the neighbor values.

At 0.9 the slope weight nearly vanishes. Put minmod back and that overshoot gets clipped.

The equation is the simplest conservative advection there is. A divergence-free velocity field $\mathbf{u}$ carries a scalar $s$.

Here $s$ is a passive scalar like concentration or temperature, and $\mathbf{u}$ is a prescribed velocity field. Updating cell averages with finite volumes (FV) gives this form.

$F^x = u\,s^{n+1/2}$ is a time-centered flux at the face. All the accuracy comes down to how well you estimate the face value $s^{n+1/2}$.

Splitting one direction at a time twists rotation#

Dimensional splitting breaks a 2D problem into an x-direction 1D update followed by a y-direction 1D update. The implementation is easy — a single 1D solver does it all. The trouble is diagonal flow.

When a scalar moves at $45^\circ$, doing x first then y gives a different answer than the reverse order. Both miss the cross term that lives between the two 1D steps. In a rotating flow this error piles up every step, and a round disk warps into a rhombus. Alternating the order with Strang splitting only cancels the directional bias of the second-order error; the root cause stays.

BCG's prescription is plain: solve x and y at the same time. Within one step, predict every face value at the same time level $t^{n+1/2}$, and bake the transverse change into that prediction.

Turning time into space in the predictor#

The core idea is to push the face value forward by half a step in both time and space with a Taylor expansion. At the right face of cell $i$, the left state reads:

Leaving the time derivative as is would be useless. Substitute the advection equation $\partial_t s = -u\,\partial_x s$ to trade the time derivative for a spatial one (the Cauchy-Kowalewski trick).

$\Delta_x s_{ij}$ is the limited slope in cell $i$, and $u_{i+1/2,j}$ is the advecting face velocity. The $(1 - \mathrm{CFL})$ factor carries the time advance. As CFL approaches 1 the left state's contribution shrinks; at 0 it becomes a plain spatial extrapolation.

The right state seen from the neighboring cell is built symmetrically.

Choosing between the two candidates is Godunov's upwind rule. The sign of the face velocity decides the wind.

In the figure above, CFL 0.9 pulls the amber dot toward the average, and none lets the line overshoot the neighbors. That gap is the time-advance factor in the predictor just above.

The transverse derivative stitches the corners#

BCG injects the coupling that splitting drops straight into the predictor. From the left state at an x-face, subtract half a step of the transverse flux divergence.

$v$ is the y-velocity and $\partial_y(v s)$ is the transverse advective change in cell $ij$. To estimate this term you first need the 1D predicted states $\hat s^y$ on the y-faces. So BCG has two passes. First predict the 1D upwind state in each direction, then reuse that prediction for the transverse correction, then pick the final face value with the same upwind rule. This "corner coupling" carries diagonal advection without the twist.

Putting a leash on the slope — minmod and superbee#

The predictor's accuracy comes from the slope $\Delta_x s$. Using a raw centered slope is second order but overshoots near discontinuities. A slope limiter puts a leash on it. The most conservative one, minmod, takes the smaller of the two one-sided differences, or zero when they disagree in sign.

Superbee keeps the same monotonicity but squeezes the slope as steep as it can, pulling contact fronts sharper. In the 1D demo above, flip between minmod and superbee and you will see the superbee PLM line is steeper.

Putting the same leash on the flux, in one dimension, is the TVD flux limiter on Lax–Wendroff. BCG hangs the limiter on the cell slope, not on the flux.

Spinning a disk in numpy#

Now we run the 2D unsplit BCG. The toy problem is solid-body rotation. Float a scalar disk in a divergence-free field spinning about the center, and compare first-order upwind against BCG on the same grid.

import numpy as np
 
N = 64
dx = 1.0 / N
omega = 2 * np.pi                      # one revolution per unit time
 
def face_velocities():
    j = np.arange(N)
    i = np.arange(N)
    u = -omega * ((j + 0.5) * dx - 0.5)        # x-face velocity (depends on y)
    u = np.broadcast_to(u[None, :], (N, N)).copy()
    v = omega * ((i + 0.5) * dx - 0.5)         # y-face velocity (depends on x)
    v = np.broadcast_to(v[:, None], (N, N)).copy()
    return u, v
 
def minmod_pair(a, b):                          # monotonicity-preserving slope
    return np.where(a * b > 0, np.where(np.abs(a) < np.abs(b), a, b), 0.0)
 
def limited_slopes(s):
    slx = minmod_pair(s - np.roll(s, 1, 0), np.roll(s, -1, 0) - s)
    sly = minmod_pair(s - np.roll(s, 1, 1), np.roll(s, -1, 1) - s)
    return slx, sly
 
def donor_fluxes(s, u, v):                       # first-order upwind (donor cell)
    Fx = np.where(u > 0, u * np.roll(s, 1, 0), u * s)
    Fy = np.where(v > 0, v * np.roll(s, 1, 1), v * s)
    return Fx, Fy
 
def bcg_fluxes(s, u, v, dt):                     # unsplit second-order predictor
    slx, sly = limited_slopes(s)
    # 1) normal-direction 1D predicted states
    sLx = np.roll(s, 1, 0) + 0.5 * (1 - dt/dx * u) * np.roll(slx, 1, 0)
    sRx = s - 0.5 * (1 + dt/dx * u) * slx
    hatX = np.where(u > 0, sLx, np.where(u < 0, sRx, 0.5 * (sLx + sRx)))
    sBy = np.roll(s, 1, 1) + 0.5 * (1 - dt/dx * v) * np.roll(sly, 1, 1)
    sTy = s - 0.5 * (1 + dt/dx * v) * sly
    hatY = np.where(v > 0, sBy, np.where(v < 0, sTy, 0.5 * (sBy + sTy)))
    # 2) transverse (corner) correction — reuse the predicted states
    divVy = (np.roll(v * hatY, -1, 1) - v * hatY) / dx
    divUx = (np.roll(u * hatX, -1, 0) - u * hatX) / dx
    sLx -= 0.5 * dt * np.roll(divVy, 1, 0)
    sRx -= 0.5 * dt * divVy
    sBy -= 0.5 * dt * np.roll(divUx, 1, 1)
    sTy -= 0.5 * dt * divUx
    # 3) pick the final face value with the same upwind rule
    sX = np.where(u > 0, sLx, np.where(u < 0, sRx, 0.5 * (sLx + sRx)))
    sY = np.where(v > 0, sBy, np.where(v < 0, sTy, 0.5 * (sBy + sTy)))
    return u * sX, v * sY
 
def advance(s, u, v, dt, scheme):
    Fx, Fy = bcg_fluxes(s, u, v, dt) if scheme == 'bcg' else donor_fluxes(s, u, v)
    return s - dt * ((np.roll(Fx, -1, 0) - Fx) / dx + (np.roll(Fy, -1, 1) - Fy) / dx)
 
def make_disk():
    X, Y = np.meshgrid((np.arange(N)+0.5)*dx, (np.arange(N)+0.5)*dx, indexing='ij')
    return ((X - 0.5)**2 + (Y - 0.78)**2 < 0.13**2).astype(float)
 
u, v = face_velocities()
Umax = omega * 0.5 * np.sqrt(2)
dt = 0.6 * dx / Umax                    # CFL = 0.6
steps = int(round(1.0 / dt))            # exactly one revolution
 
for scheme in ('upwind1', 'bcg'):
    s = make_disk()
    mass0 = s.sum()
    for _ in range(steps):
        s = advance(s, u, v, dt, scheme)
    print(f"{scheme:8s}  peak={s.max():.3f}  "
          f"mass_err={abs(s.sum()-mass0)/mass0:.2e}")

The output looks like this.

upwind1   peak=0.349  mass_err=3.1e-16
bcg       peak=0.806  mass_err=2.8e-16

Flip the scheme toggle from 1st-order upwind to BCG unsplit and the peak walks from 0.349 to 0.806 on the grid, the same numbers as the printout.

On 1st-order upwind the disk spreads into a ring before it finishes one turn, edge blurred. BCG unsplit keeps the boundary much longer. At CFL 0.9, BCG roughens a little. The gap over first-order upwind stays.

Both schemes conserve mass to machine precision. Conservative flux differencing does that. What differs is the peak. First-order upwind crushes the disk to less than half its height in one revolution, while BCG keeps 80% of the original height.

What is left after two turns#

The code already fixes one turn. N is 64, so the spacing is 1/64. Corner speed is 2π times 0.5 times the square root of 2, about 4.44. At CFL 0.6 the step is about 0.00211, and time 1 takes 474 steps. The disk center (0.5, 0.78) sits 0.28 off (0.5, 0.5) and moves at about 1.76. The step uses the corner speed 4.44, not 1.76.

Diameter 0.26 is 16.6 cells, about 217 cells of scalar. At CFL 0.6 the left slope weight is 0.2 and the right is 0.8. At the cap 0.8 the left weight falls to 0.1. Peaks 0.349 and 0.806 are those weights plus the transverse term, stacked 474 times.

Mass errors 3.1e-16 and 2.8e-16 are roundoff. Lost height is 0.651 against 0.194, about 3.4 times more for upwind than for BCG. The transverse term subtracts about 0.00106 times the transverse divergence, once a step, 474 times a turn. Smearing is how much slope one cell keeps, not the choice of 64 cells.

Where it leaks when you port it#

• The CFL is tighter than a split scheme. Unsplit 2D propagates information along the diagonal, so its stability limit is smaller than the 1D one. In practice, start with $\mathrm{CFL}\lesssim 0.8$.
• Always store velocity on the faces. Interpolating cell-centered velocity to faces breaks the divergence-free property, and the mass you thought was conserved leaks out. A MAC (staggered) layout is the answer.
• Do not drop the transverse term. Implement only the normal prediction and you have plain dimensionally split MUSCL. The code runs, but the rotational twist stays.
• The limiter goes to zero at extrema. At the peak of a hump the slope gets clipped to zero, dropping to locally first order. That is why the peak slowly sags. To preserve it, look at extremum-preserving limiters.

The choice where mass is right and the disk still dies#

A peak of 0.349 is what first-order upwind has left after one turn. BCG sits at 0.806. The mass errors are 3.1e-16 and 2.8e-16, so a conservation log passes both. Look at the disk and only one of them is still a disk.

Whether first-order upwind is acceptable on a diagonal scalar is that gap. If the contact height can be thrown away, the donor cell is enough. Keeping the height means the transverse term has to enter the same step. Move the mesh off Cartesian and minmod does not travel for free. The Barth–Jespersen limiter is the next decision: how a fake extremum gets cut by the neighbor minimum and maximum.