Where Jacobi Dies on a Boundary-Layer Mesh — the Linelet Line-Implicit Preconditioner

Same equation. Same cell count, same iteration, same convergence criterion. The mesh got squashed toward the wall once, and the iteration count jumped from 18 to 5,000.

None of the physics changed. The only thing that changed is the cell aspect ratio. Yet the linear solver starts behaving as if it were handed a different problem. And you cannot dodge this: hitting $y^+ \approx 1$ in a turbulent run means the first cell off the wall has to be squashed.

That factor of 300 has an address inside the matrix. The linelet preconditioner — a chain of cells threaded along the strongly coupled direction, the one SU2 uses — hands it back. The sweep counter is the check. Short rows are the table. Long rows are the coefficients.

Squash a cell and the matrix becomes a different object#

Discretize the diffusion term over a rectangular cell with finite volumes and the face coefficients come out as

where $\Gamma$ is the diffusion coefficient and $\Delta x$, $\Delta y$ are the cell dimensions. Now take the ratio.

$\mathrm{AR}$ is the cell aspect ratio (streamwise over wall-normal). The coefficient ratio goes as the square of it. A first cell at $\mathrm{AR} = 1000$ has north-south coupling a million times stronger than east-west. Squash the mesh by ten and the matrix skews by a hundred.

Drag the AR slider from 1 to 100 and the vertical stripes die first. Toggle linelet and the leftover horizontal smear loses its lifetime.

When the measured readout under the map sticks near 0.997, the stall is on. At the same AR, linelet clears that smear in a sweep or two.

Add an implicit time derivative and the diagonal becomes

with $V = \Delta x \Delta y$ the cell volume and $\Delta t$ the (pseudo-)time step. That $V/\Delta t$ term is the only thing propping up the diagonal, and it cannot keep pace with $a_N$.

The point-smoother factor has AR in it#

The error decay of a Jacobi sweep is bounded above by the off-diagonal sum over the diagonal.

Set $V/\Delta t = 4a_E$ and this reads 0.5 at $\mathrm{AR} = 1$ and 0.9998 at $\mathrm{AR} = 100$. The trouble is that the bound glues itself to 1 as $\mathrm{AR}$ grows and stops telling you anything. To see what actually survives, go mode by mode. For a Fourier mode $(\theta_x, \theta_y)$ the damping factor is

The slowest mode is the one with the smallest $\theta_y$ — the gentlest variation in the wall-normal direction. Once $a_N$ owns the diagonal, $a_P \approx 2a_N$, and

where $N_y$ is the number of wall-normal layers. The absence of $\mathrm{AR}$ here is bad news, not good: squashing further will not push the factor higher, but the ceiling it saturates at is already useless. With $N_y = 48$ that is 0.99786, or roughly 6,000 sweeps to drive the error down by $10^{-6}$.

What is happening is easy to picture. A point smoother wipes out wall-normal high-frequency error in a sweep or two. Then it is stuck. Killing error that varies along the wall requires information to travel sideways through $a_E$, and since the diagonal is pinned by $a_N$, each step shrinks by $\mathrm{AR}^2$.

Where a diagonal-only iteration stalls, method by method, is in a sweep-count comparison from Jacobi through BiCGSTAB. That piece changes the solver. Here the cell aspect ratio betrays the same solver.

Three remedies on one table#

The linelet's bargain is explicit. It inverts exactly one direction — the strongly coupled one. In exchange, storage is a single three-band array and the cost lands in the same order of magnitude as Jacobi. Because it never touches the full matrix the way ILU(0) does, it is also far less sensitive to MPI partitioning.

How the linelets get built#

1. For each cell, measure the face coefficient ratio $\sigma = a_f / a_{f,\perp}$. For a wall-normal face this is $\sigma = \mathrm{AR}^2$.
2. Mark a face strong when $\sigma \ge \sigma_{\min}$. A common default sits near $\sigma_{\min} = 10$.
3. Start at a cell touching the wall and extend a chain upward along strong faces.
4. Allow at most two strong faces per cell. If a branch appears, keep only the larger coefficient — allowing branches means the result is not tridiagonal.
5. Cut the chain the moment it reaches a cell that already belongs to another linelet.
6. Discard chains of length one and hand those cells back to Jacobi.

Renumber inside a chain and that submatrix becomes tridiagonal. Tridiagonal means the Thomas algorithm (forward elimination, then back substitution) inverts it exactly in $O(n)$. That sentence is the entire linelet preconditioner.

Where a chain stops sits next to the sweep table. The controls are first-layer thickness and growth ratio.

Counting sweeps in Python#

On a $64 \times 48$ mesh, leave nothing but error behind (the right-hand side is $f = 0$, so the exact answer is $u = 0$), then run point Jacobi and line Gauss-Seidel and count sweeps until the error drops by $10^{-6}$.

import numpy as np
 
def cell_coefficients(dx, dy, gamma=1.0, diag_factor=4.0):
    """FVM diffusion coefficients of one rectangular cell. diag_factor sets V/dt as a multiple of a_E."""
    a_e = gamma * dy / dx        # east/west faces
    a_n = gamma * dx / dy        # north/south faces
    d0 = diag_factor * a_e       # the V/dt term
    return a_e, a_n, d0
 
def residual_field(u, a_e, a_n, d0):
    r = -d0 * u
    r[1:-1, 1:-1] -= a_e * (2*u[1:-1, 1:-1] - u[:-2, 1:-1] - u[2:, 1:-1])
    r[1:-1, 1:-1] -= a_n * (2*u[1:-1, 1:-1] - u[1:-1, :-2] - u[1:-1, 2:])
    r[0, :] = r[-1, :] = 0.0     # Dirichlet boundaries
    r[:, 0] = r[:, -1] = 0.0
    return r
 
def jacobi_sweep(u, a_e, a_n, d0, omega=1.0):
    r = residual_field(u, a_e, a_n, d0)
    u += omega * r / (d0 + 2*a_e + 2*a_n)
 
def thomas(a, b, c, d):
    """Exact tridiagonal solve: forward elimination, then back substitution."""
    n = len(d)
    cp, dp = np.empty(n), np.empty(n)
    cp[0], dp[0] = c[0]/b[0], d[0]/b[0]
    for k in range(1, n):
        m = b[k] - a[k]*cp[k-1]
        cp[k] = c[k]/m
        dp[k] = (d[k] - a[k]*dp[k-1])/m
    x = np.empty(n)
    x[-1] = dp[-1]
    for k in range(n-2, -1, -1):
        x[k] = dp[k] - cp[k]*x[k+1]
    return x
 
def linelet_sweep(u, a_e, a_n, d0):
    """Invert one whole wall-normal chain exactly (line Gauss-Seidel)."""
    nx, ny = u.shape
    m = ny - 2
    dg = d0 + 2*a_e + 2*a_n
    a = np.full(m, -a_n)
    b = np.full(m, dg)
    c = np.full(m, -a_n)
    a[0] = c[-1] = 0.0
    for i in range(1, nx-1):
        rhs = a_e * (u[i-1, 1:-1] + u[i+1, 1:-1])   # off-line coupling goes to the RHS
        u[i, 1:-1] = thomas(a, b, c, rhs)
 
def count_sweeps(ar, kind, nx=64, ny=48, tol=1e-6, max_sweeps=20000):
    dx, dy = 1.0/nx, 1.0/nx/ar
    a_e, a_n, d0 = cell_coefficients(dx, dy)
    u = np.random.default_rng(7).standard_normal((nx, ny))
    u[0, :] = u[-1, :] = 0.0
    u[:, 0] = u[:, -1] = 0.0
    e0 = prev = np.linalg.norm(u)
    for k in range(1, max_sweeps + 1):
        jacobi_sweep(u, a_e, a_n, d0) if kind == 'jacobi' else linelet_sweep(u, a_e, a_n, d0)
        e = np.linalg.norm(u)
        rho, prev = e/prev, e
        if e/e0 < tol:
            return k, rho
    return max_sweeps, rho
 
print(f"{'AR':>6} {'Jacobi':>8} {'rho_J':>8} {'linelet':>8} {'rho_L':>8}")
for ar in (1, 10, 100, 1000):
    nj, rj = count_sweeps(ar, 'jacobi')
    nl, rl = count_sweeps(ar, 'linelet')
    print(f"{ar:>6} {nj:>8} {rj:>8.4f} {nl:>8} {rl:>8.4f}")

Here is what it prints.

    AR   Jacobi    rho_J  linelet    rho_L
     1       18   0.4870        8   0.1839
    10      514   0.9780        7   0.1699
   100     4879   0.9975        4   0.0194
  1000     5471   0.9978        2   0.0002

Shrink $\Delta y_{\text{wall}}$ and the chains dig deeper into the layer. Raise the growth ratio and the mesh returns to isotropy sooner, so the chains get shorter. The chains are short. They do not cover half the layers. That length is why linelet stays cheap.

Squash by a factor of 1000 and Jacobi slows by 300× before hitting the ceiling near 6,000 that the previous section predicted. Linelet gets faster. The coupling between neighboring chains falls off as $a_E/a_N = \mathrm{AR}^{-2}$, so the chains become independent of one another. The worst-mode bound $2a_E/(V/\Delta t + 2a_E) = 1/3$ still holds; the measured decay is simply much better than it.

Split the Jacobi column. 18, 514, 4879, 5471. The aspect ratios are 1, 10, 100, 1000.

The first factor of ten in aspect ratio multiplies the work by 28.6. The next factor of ten multiplies it by 9.5. The last one multiplies it by 1.12, which is 5471 divided by 4879. The coefficient ratio keeps growing as $\mathrm{AR}^2$, and the sweep count stops following. At $N_y = 48$ the mode formula $1 - \frac{\pi^{2}}{2 N_y^{2}}$ is 0.99786. Solving $0.99786^{n} = 10^{-6}$ puts n near 6400. The printed 5471 sits just under that ceiling. Squash the mesh further and this column stays glued there.

The linelet column runs the other way. 8, 7, 4, 2. The measured factors are 0.1839, 0.1699, 0.0194, 0.0002. The bound $2a_E/(V/\Delta t + 2a_E) = 1/3$ has no $\mathrm{AR}$ in it, and every measured value sits under it. At AR 1000, 0.0002 is about 0.06% of 1/3. Coupling between chains, $a_E/a_N = \mathrm{AR}^{-2}$, shuts off, so one Thomas pass leaves almost none of the wall-normal error. At AR 1 it is 8 sweeps against 18. A small gain. At AR 10 it is 514 against 7, about 73 times. That is where the choice splits.

The point-smoother bound uses the same diagonal. With $V/\Delta t = 4a_E$ the denominator is $6a_E + 2a_N$. At AR 1 the bound is 0.5 and the table's rho_J is 0.4870. Nearly the same number. At AR 100 the bound is 0.9998 and the measured rho_J is 0.9975. Once the bound sticks to 1, the mode formula's 0.99786 is the honest one. The earlier estimate of roughly 6,000 sweeps and the printed 5471 are the same ceiling.

Traps to stay out of when you code this#

Laying the lines the wrong way. Chains must follow the direction with the larger coefficient, which is wall-normal. Lay them along the wall and you gain nothing while paying for Thomas. It is easy to get backwards, because the mesh visually looks "long" in the streamwise direction.

Leaving $\sigma_{\min}$ at its default. Change the mesh without revisiting the threshold and the chains can disappear entirely. The preconditioner then degrades silently into Jacobi, and nothing in the log warns you. Printing the linelet count and mean length on every run catches this in one line.

Letting MPI partitions cut the chains. Chains severed at a partition boundary come out in fragments, and convergence gets worse the more cores you add. Give the partitioner a weight along the line direction, or constrain it not to cut lines at all.

Trying to paper over a nonlinear problem. If it still diverges with linelet in place, the culprit is usually the outer loop rather than the linear solver. Lower the CFL first, pull the under-relaxation factor down toward 0.7 so the updates themselves shrink, and only then judge. A preconditioner solves the matrix you gave it faster; it does not repair a wrong one.

Thomas Is Wasted If the Lines Are Not Wall-Normal#

Chains follow the large coefficient, wall-normal only. Lay them streamwise and the eye has simply tracked the long edge. Thomas still costs something. The damping matches Jacobi.

Log the chain count and the mean length. A printed zero means $\sigma_{\min}$ no longer matches the mesh. The preconditioner is already a diagonal smoother. Adding cores in that state only chops the fragments shorter. Weak scaling that drops to 16.3% efficiency while the mesh is refined is the bill for that cut.