When Explicit Blows Up, Solve Implicit Without Inverting the Matrix — LU-SGS

The run died at 3 a.m. Last log line: pressure residual = NaN. CFL went from 0.9 to 0.5, then 0.3. Each step was so short that steady state meant days. Familiar wall in compressible steady work. Implicit time-stepping breaks that wall on the step-size side. The version that never inverts the full matrix is LU-SGS. A Burgers shock still exists at CFL 5.

Toggle the scheme between forward Euler and backward Euler, and drag $\mathrm{Re}(\lambda\Delta t)$ to $-2.5$. The stable set changes from a small disk to the left half-plane.

On forward Euler that point leaves the cyan disk and flips to UNSTABLE. Backward Euler fills the whole left half-plane with cyan, so the same coordinate is stable again.

The stability region opens to infinity#

Time-stepping stability is decided by a single model equation.

Here $\lambda$ is the eigenvalue produced by the spatial discretization (it comes from wave speeds and viscosity), and $u$ is the amplitude of one mode. After one step the amplitude is multiplied by the amplification factor $G$. We need $|G|\le 1$ to avoid blow-up.

Explicit stepping (forward Euler) gives

so $|1+\lambda\Delta t|\le 1$ — the point $\lambda\Delta t$ must sit inside a small disk centered at $-1$ with radius $1$. A viscous grid or a thin cell makes $\lambda$ large, and it leaves that disk fast.

Implicit stepping (backward Euler) is the opposite.

For every physical decaying mode with $\mathrm{Re}(\lambda)<0$, $|G|<1$. It stays stable no matter how large $\Delta t$ grows. The stability constraint dissolves, and now the time step only has to be as small as accuracy demands.

The explicit disk ends at 0 and $-2$ on the real axis. Center $-1$, radius 1. $\mathrm{Re}(\lambda\Delta t)=-2.5$ sits 1.5 from the center, outside the disk, which is why the marker reads UNSTABLE. Backward Euler uses $G=1/(1-\lambda\Delta t)$. When $\mathrm{Re}(\lambda)<0$ the denominator exceeds 1, so $|G|$ stays under 1 as $\Delta t$ grows. The night log cut CFL through 0.9, 0.5, and 0.3 to trap $\lambda\Delta t$ inside that disk. On a viscous grid $\lambda$ grows first. The disk breaks before accuracy asks for a shorter step.

Inverting the full matrix is a cost you can't pay#

Nothing is free. Implicit stepping asks for the residual at the future state.

$R_i$ is the flux divergence (residual) of cell $i$, and $U$ is the conserved variable. Linearize $R(U^{n+1})$ about the current state,

and a linear system for $\Delta U = U^{n+1}-U^n$ remains.

The trouble is the matrix $A$ on the left. With fifty million cells, $A$ has fifty million dimensions. A direct inverse is impossible in both memory and time. A Krylov solver like GMRES works too, but it has to store the matrix explicitly and needs a preconditioner. In the 1990s compressible analysts found a cheaper path: approximate the matrix with no storage, no multiplications, and two sweeps.

LU-SGS — split into D, L, U and sweep twice#

Decompose $A$ into a block lower triangle $L$, a block diagonal $D$, and a block upper triangle $U$.

LU-SGS (Lower-Upper Symmetric Gauss-Seidel) approximates this $A$ by the product

which discards the $L D^{-1} U$ term created by the $LU$ product. When $\Delta U$ is small it's negligible. The solve is now two triangular back-substitutions in sequence.

Forward sweep — lower triangle, front to back:

Backward sweep — upper triangle, back to front:

What fills $D$, $L$, and $U$ decides how good the approximation is. Yoon and Jameson (1988) approximated the nonlinear flux Jacobian with its spectral radius $\rho$. Splitting first-order upwind as $A^{\pm} = \tfrac{1}{2}(a \pm |a|)$ leaves the diagonal $D_i = I/\Delta t + |a_i|/\Delta x$ and off-diagonals that carry only the split face Jacobians. You never build the whole matrix — a few scalars (or a small block) per cell are enough.

The approximate product is $(D+L)\,D^{-1}\,(D+U)$. Expand it and $L D^{-1} U$ sits next to $A$; that term is dropped. Because $D_i = I/\Delta t + |a_i|/\Delta x$, a small $\Delta t$ inflates $I/\Delta t$ and the dropped term shrinks. Raise CFL to 5 and $\Delta t$ grows, so $I/\Delta t$ shrinks and the approximation error grows with it. Explicit stepping still hits NaN in two or three steps on that interval. The printed implicit run crosses it in 49 steps with $u$ closed in $[0.012, 0.988]$, on $n=200$ cells out to $t=0.22$ with $\Delta x=1/200$. A stored-matrix path such as the Krylov convergence comparison keeps the matrix. Here a diagonal division and two sweeps take that seat.

Python — advancing implicit Burgers with one factorization#

The toy problem is shock formation in the 1D Burgers equation. A smooth sine wave steepens into a shock over time. Explicit stepping would pin the CFL below 1. The LU-SGS implicit step survives at CFL 5.

import numpy as np
 
def upwind_flux_residual(u, dx):
    """First-order upwind flux for the residual R_i of d(u^2/2)/dx (periodic)."""
    f = 0.5 * u * u
    a = u                       # local wave speed df/du
    f_left  = np.where(a >= 0, np.roll(f, 1), f)
    a_right = np.roll(a, -1)
    f_right = np.where(a_right >= 0, f, np.roll(f, -1))
    return (f_right - f_left) / dx
 
def lu_sgs_burgers_step(u, dt, dx):
    """One time step: advance the implicit system with a single LU-SGS pass."""
    n = u.size
    a  = u
    ap = 0.5 * (a + np.abs(a))         # A+  (lower-triangle contribution)
    am = 0.5 * (a - np.abs(a))         # A-  (upper-triangle contribution)
    D  = 1.0 / dt + np.abs(a) / dx     # diagonal block
    R  = upwind_flux_residual(u, dx)
 
    dus = np.zeros(n)                  # forward sweep: (D+L) dU* = -R
    for i in range(n):
        lower = ap[i - 1] / dx * dus[i - 1]
        dus[i] = (-R[i] + lower) / D[i]
 
    du = dus.copy()                   # backward sweep: (D+U) dU = D dU*
    for i in range(n - 1, -1, -1):
        upper = am[(i + 1) % n] / dx * du[(i + 1) % n]
        du[i] = dus[i] - upper / D[i]
    return u + du
 
def drive_burgers_shock(n=200, cfl=5.0, tmax=0.22):
    dx = 1.0 / n
    x  = (np.arange(n) + 0.5) * dx
    u  = 0.5 + 0.5 * np.sin(2 * np.pi * x)   # smooth initial condition
    t, steps = 0.0, 0
    while t < tmax:
        dt = cfl * dx / np.max(np.abs(u))    # CFL>1 allowed (implicit)
        dt = min(dt, tmax - t)
        u = lu_sgs_burgers_step(u, dt, dx)
        t += dt; steps += 1
    return x, u, steps
 
x, u, steps = drive_burgers_shock(cfl=5.0)
print(f"CFL=5 -> {steps} steps,  u in [{u.min():.3f}, {u.max():.3f}]")
# CFL=5 -> 49 steps,  u in [0.012, 0.988]

Run the same problem with explicit forward Euler at CFL 5 and it produces NaN within two or three steps. The implicit version reaches the shock in 49 stable steps. No heavy Riemann solver per face — one diagonal division and two sweeps did the whole job.

Sweep after sweep, the residual collapses#

A single LU-SGS application gives an approximate solution. To solve more accurately, you repeat the sweep on the same system. The convergence rate then depends on the diagonal dominance of the matrix. In the viscous-dominated case — where the viscous Jacobian becomes a 1D Laplacian — the diagonal is $1+2d$ and the off-diagonals are $-d$. Here $d=\nu\Delta t/\Delta x^2$ is the diffusion number (the nondimensional time step of viscous diffusion).

At $d=\nu\Delta t/\Delta x^2=1$ the diagonal is $1+2d=3$ and each off-diagonal is $-d=-1$. The diagonal beats the sum of the two neighbors, and a few sweeps drop the residual ten orders. At $d=40$ the diagonal is 81 and the off-diagonals are $-40$ twice, so the row sum is $81-80=1$ and nearly cancels. That is why the curve lies down. The stability bought with a large $\Delta t$ and the sweep count the inner relaxation needs are written on the same ledger.

Slide $d$ from near 1 up to 40 and the slope of the residual curve changes at the same sweep count.

Keep $d$ near 1 and the residual drops ten orders in a handful of sweeps. Crank $d$ up to 40 (a huge implicit $\Delta t$) and the same curve flattens. You bought stability by enlarging the step, but the matrix got stiffer, so the inner relaxation needs more sweeps. Stability and inner convergence do not arrive free together.

Drop the viscous radius from the diagonal and the sweep dies#

Leave $\nu/\Delta x^2$ out of the diagonal $D$ and $D$ shrinks in a viscous cell. The spectral radius is then only the convective $|a|/\Delta x$. The sweep diverges at once. If the viscous Jacobian is a 1D Laplacian, diagonal dominance is the balance of $1+2d$ against $-d$, and stripping the viscous term empties the left side of that balance. When the nonlinearity is strong, do not add the full $\Delta U$ at once. Use $U^{n+1}=U^n+\omega\,\Delta U$ with $\omega\in[0.5,0.8]$. Oscillations near shocks shrink.

Direction is a decision too. Repeating only the forward sweep sends information toward increasing indices. Alternating forward and backward keeps $L$ and $U$ symmetric. Across an MPI partition, one communication has to sit between the two sweeps or the neighbor $\Delta U$ never arrives. How the cells are numbered changes that sweep order. The record of LU-SGS slowing after an RCM renumbering is the other side of the same choice.