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

A steady residual is touchy about CFL. Leave it at 1 and ten thousand iterations still sit high. Push it to 100 and the first step blows the solution apart. SER (Switched Evolution Relaxation) lets the residual pick the step between those two. Down means a larger step. A spike means a smaller one.

Drag the fixed-CFL slider from 1 to 20 and the amber residual jumps to diverged, while SER only exposes the starting CFL and the exponent.

At $\mathrm{CFL}_0=1$ and $p=1$ the cyan curve bends down late. Raising $p$ to 1.5 steepens that drop, but the moment the early residual wobbles, cyan wobbles with it.

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.

With N=80 the spacing is 1/81. At $\nu=0.02$, $2\nu/\Delta x$ is 3.24, larger than $|u_i|$ on the linear guess. Fixed CFL 3 then keeps $\Delta\tau$ under 0.01. The cap 1e5 pushes $\Delta\tau$ into the hundreds, so $I/\Delta\tau$ nearly vanishes. The code shares the maximum of $|u|$ across every cell, even in the shock where the speed is near 0. That is stricter than the local step in the formula.

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. The same back-substitution that makes 1D diffusion unconditionally stable is where this left-hand side picks up its time term.

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

Drop the cap slider and the schedule shows where the printed 34 iterations start getting long again. Cyan is CFL, 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.

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#

588 divided by 34 is about 17. The tolerance is a residual ratio below 1e-9. With $p=1$, a starting CFL of 1, and a cap of 1e5, the ceiling hits once the residual ratio reaches 1e5, that is 1e-5. Four decades remain until 1e-9, all at CFL 1e5, already on the Newton side. Fixed CFL 3 spends the whole 588 steps without that stretch.

The schedule is one line. $\mathrm{CFL}_n = \min(\mathrm{CFL}_0 (\lVert R_0\rVert/\lVert R_n\rVert)^p,\ \mathrm{CFL}_{\max})$. As the residual falls the step grows, then it stops at the cap.

Pitfalls when you switch SER on#

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

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.

$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$.

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.

An open cap kills the late Newton first#

Leave $\mathrm{CFL}_{\max}$ unset and one residual spike lets the CFL chase the residual ratio without a ceiling. The implicit pseudo-time iteration becomes Newton's method as $\Delta\tau\to\infty$. That quadratic rate shows up only near the solution. Jump there early and the matrix shoves the solution the wrong way. The cap of 1e5 in the code is a door: the run may enter that region only after the residual has fallen by a hundred thousand.

The same door matters when the Jacobian cannot be formed exactly every iteration. Once the large-CFL quadratic rate has already fallen to first order, the remaining choice is a lower cap. Approximating that matrix with sweeps, without assembling it, is the point of LU-SGS after an explicit step blows up. The implicit step stays, and the matrix does not have to be stored.