CFL = 1 takes ten thousand steps, CFL = 100 explodes — Switched Evolution Relaxation

Set the CFL number to 1 in a steady-state run and the residual won't drop even after ten thousand iterations. Push it to 100 and the solution blows up on the very first step. Both are everyday sights. This post is about SER (Switched Evolution Relaxation), a trick that walks between those two extremes automatically. It grows the CFL as the residual falls — safe early, then as fast as Newton near the end — and we'll feel it on a small steady Burgers solver.

In a steady solver the CFL number (the time step made dimensionless by grid size and wave speed) is no longer about accuracy. It's a single knob that holds both convergence speed and stability at once.

Steady state is chased with fake time#

A steady equation is one line: $R(U) = 0$. Here $R$ gathers convection, diffusion, and sources into a spatial residual. Solving it directly means cracking a nonlinear algebraic system. Instead, most CFD codes bolt on a fake-time term.

Here $\tau$ is pseudo-time with no physical meaning. As $\tau \to \infty$, $\partial U / \partial \tau \to 0$, so $R(U)=0$ — exactly the steady solution we wanted. We throw away time accuracy and aim only for the steady state. This is pseudo-transient continuation.

March the fake time implicitly (backward Euler) and linearize the residual, and each iteration becomes this:

Here $\partial R/\partial U$ is the flux Jacobian and $\Delta U = U^{n+1}-U^n$. Each iteration solves the left-hand matrix once for $\Delta U$.

Start small, land on Newton#

That single line holds all of SER's intuition. Push $\Delta\tau$ to both extremes.

When $\Delta\tau \to 0$, the term $I/\Delta\tau$ dominates the matrix and $\Delta U \approx -\Delta\tau\, R$. That's plain explicit marching. Rock-solid but slow as a turtle.

When $\Delta\tau \to \infty$, the term $I/\Delta\tau$ vanishes and only $(\partial R/\partial U)\,\Delta U = -R$ remains. That's exactly Newton's method, with quadratic convergence (error shrinks by its square) near the solution. But a bad initial guess makes it diverge.

So the ideal strategy is clear. Keep $\Delta\tau$ small early to survive the nonlinearity, then grow it once the residual drops and you're near the solution, to soak up Newton's quadratic convergence. The catch is "when, and by how much." SER hands that judgment to the residual.

SER — grow CFL by the inverse of the residual#

The rule Mulder and van Leer proposed in 1985 is remarkably simple. Grow the current CFL by the ratio of the initial to the current residual.

$\mathrm{CFL}_0$ is the starting CFL, $\lVert R_n \rVert$ is the L2 norm of the $n$-th residual, $p$ is an exponent usually around 1, and $\mathrm{CFL}_{\max}$ caps runaway. When the residual drops tenfold (with $p=1$), the CFL grows tenfold too. The pseudo-time step is converted from this CFL and the local wave speed.

The denominator sums the convective wave speed $|u_i|$ and the diffusive wave speed $2\nu/\Delta x$. It can differ per cell, so it pairs naturally with local time-stepping.

Play with the parameters in the simulation below. It's the residual history of fixed CFL and SER run side by side on the same steady Burgers problem.

Drop fixed CFL to 1 and the amber curve barely reaches the floor. Push it to 20 and it jumps to "diverged." SER (cyan), by contrast, starts safely at $\mathrm{CFL}_0=1$ yet bends sharply downward late in the run. Raising the exponent $p$ from 1 to 1.5 makes it more aggressive, but push it too far and even SER wobbles early.

Python — solving steady Burgers two ways#

We solve the steady viscous Burgers equation $u u_x = \nu u_{xx}$. The boundary conditions $u(0)=1$, $u(1)=-1$ create a stationary shock layer in the middle. One implicit pseudo-time step solves a tridiagonal matrix with the Thomas algorithm.

import numpy as np
 
def steady_residual(u, nu, dx):
    """R(u) = u u_x - nu u_xx, boundaries u(0)=1, u(1)=-1."""
    um = np.concatenate(([1.0], u[:-1]))   # left neighbor
    up = np.concatenate((u[1:], [-1.0]))   # right neighbor
    conv = u * (up - um) / (2 * dx)
    diff = nu * (up - 2 * u + um) / dx**2
    return conv - diff
 
def thomas_solve(a, b, c, d):
    """Solve a tridiagonal system (sub/main/super diag a,b,c, rhs d)."""
    n = len(d)
    b, d = b.copy(), d.copy()
    for i in range(1, n):
        m = a[i] / b[i - 1]
        b[i] -= m * c[i - 1]
        d[i] -= m * d[i - 1]
    x = np.empty(n)
    x[-1] = d[-1] / b[-1]
    for i in range(n - 2, -1, -1):
        x[i] = (d[i] - c[i] * x[i + 1]) / b[i]
    return x
 
def ptc_update(u, nu, dx, dt):
    """One implicit pseudo-time step: (I/dt + dR/du) du = -R."""
    n = len(u)
    um = np.concatenate(([1.0], u[:-1]))
    up = np.concatenate((u[1:], [-1.0]))
    a = -u / (2 * dx) - nu / dx**2                       # sub-diagonal
    b = (up - um) / (2 * dx) + 2 * nu / dx**2 + 1.0 / dt  # main diagonal
    c = u / (2 * dx) - nu / dx**2                        # super-diagonal
    du = thomas_solve(a, b, c, -steady_residual(u, nu, dx))
    return u + du
 
def ser_schedule(r0, rn, cfl0, p, cfl_max):
    return min(cfl0 * (r0 / rn) ** p, cfl_max)
 
def converge(nu=0.02, ser=True, cfl0=1.0, p=1.0, cfl_max=1e5,
             N=80, tol=1e-9, max_iter=600):
    dx = 1.0 / (N + 1)
    x = np.linspace(dx, 1 - dx, N)
    u = 1 - 2 * x                          # linear initial guess
    r0 = np.linalg.norm(steady_residual(u, nu, dx)) / np.sqrt(N)
    rn = r0
    for it in range(1, max_iter + 1):
        cfl = ser_schedule(r0, rn, cfl0, p, cfl_max) if ser else cfl0
        dt = cfl * dx / (np.abs(u).max() + 2 * nu / dx)
        u = ptc_update(u, nu, dx, dt)
        rn = np.linalg.norm(steady_residual(u, nu, dx)) / np.sqrt(N)
        if rn / r0 < tol:
            return u, it
    return u, max_iter
 
for tag, kw in [("fixed CFL=3", dict(ser=False, cfl0=3.0)),
                ("SER  p=1.0", dict(ser=True, cfl0=1.0, p=1.0))]:
    _, iters = converge(**kw)
    print(f"{tag:14s} -> {iters:4d} iters")

The output looks like this.

fixed CFL=3    ->  588 iters
SER  p=1.0     ->   34 iters

To reach a steady solution of the same accuracy, fixed CFL needs hundreds of iterations, SER just dozens. The gap opens up in the late phase, because SER drives the CFL into the thousands and effectively turns into Newton's method.

Fixed CFL and SER, side by side#

Let's watch the schedule itself — how SER pulls the CFL up. Play with it below. Cyan is the CFL schedule, pink is the residual.

While the residual is flat early on, the CFL stays near $\mathrm{CFL}_0$. The instant the residual starts collapsing, the CFL shoots up almost vertically. That is the essence of SER: the better the solution gets, the larger the step it allows. Lower the cap to 20 and the CFL hits the ceiling, loses the late Newton acceleration, and the iteration count climbs back up.

Pitfalls when you switch SER on#

First, always set $\mathrm{CFL}_{\max}$. Leave it at infinity and a single residual spike sends the CFL runaway, killing the solution.

Second, when the residual wobbles non-monotonically, $R_0/R_n$ can drop below 1 and shrink the CFL. That's not a bug — it's a safety valve. If the residual gets worse again, shortening the step is the right move.

Third, $p$ is not free acceleration. A $p>1$ grows the CFL hugely for even a small residual drop, raising the risk. A robust start is $p=1$, $\mathrm{CFL}_0 \approx 1$.

Fourth, if the Jacobian is inexact (approximate Jacobian, or matrix-free), Newton's quadratic convergence drops to first order at large CFL. Then it's better to lower $\mathrm{CFL}_{\max}$ and buy stability.

What to remember#

• In a steady solver the CFL is a knob for convergence speed and stability, not accuracy. Small means turtle, large means explosion.
• SER is $\mathrm{CFL}_n = \min(\mathrm{CFL}_0 (\lVert R_0\rVert/\lVert R_n\rVert)^p,\ \mathrm{CFL}_{\max})$, growing the step automatically as the residual drops.
• The implicit pseudo-time iteration becomes Newton's method as $\Delta\tau\to\infty$. SER is the schedule that carries you smoothly to that limit.