边界层网格上 Jacobi 停摆的地方 — 线元(linelet)线隐式预条件

点光滑器的天花板只剩下层数Ny，线元上界1/3里面没有长宽比。

同一个方程。单元数、迭代法、收敛判据都没有动。只是把网格往壁面压扁了一次，迭代次数就从 18 次跳到了 5,000 次。

物理没有任何改变，变的只有单元的长宽比。可线性求解器却开始表现得像在解另一个问题。而且这一步躲不掉：湍流计算要把 $y^+ \approx 1$ 卡住，壁面第一层单元就必须压扁。

那 300 倍在矩阵的哪一格，用系数去追。SU2 一类代码用的线元(linelet，沿强耦合方向把单元串成的链)能不能把 300 倍还回来，用迭代次数来对。短行是表，长行是系数。

压扁单元，矩阵就变成另一个东西

在矩形单元上用有限体积离散扩散项，面系数是这样的。

a_E = a_W = \frac{\Gamma \Delta y}{\Delta x}, \qquad a_N = a_S = \frac{\Gamma \Delta x}{\Delta y}

$\Gamma$ 是扩散系数， $\Delta x$ 、 $\Delta y$ 是单元尺寸。取两者之比。

\frac{a_N}{a_E} = \left(\frac{\Delta x}{\Delta y}\right)^{2} = \mathrm{AR}^{2}

$\mathrm{AR}$ 是单元长宽比(流向/壁面法向)。系数比按它的平方走。壁面第一层若 $\mathrm{AR} = 1000$ ，南北耦合就是东西耦合的一百万倍。网格压扁十倍，矩阵歪斜一百倍。

把 AR 滑块从 1 推到 100，竖条纹先消失；打开 linelet，横向剩下的色块寿命会缩短。

log₁₀ AR100ω1.00

What to watch: at AR = 1 all three iterations finish in tens of sweeps. Drag AR to 100 and the two point smoothers flatten out — the wall-normal stripes vanish almost at once, but the wall-parallel streaks just sit there and the measured factor climbs to 0.997. Switch to linelet at the same AR and it lands in single digits. The two bounds under the map say why: AR is in the point-smoother one and absent from the linelet one.

图下方 measured 贴到 0.997 附近，就是停滞。同一个 AR 下按下 linelet，色块一两个 sweep 就散掉。

再把时间导数写成隐式，对角项变为

a_P = \frac{V}{\Delta t} + 2a_E + 2a_N

其中 $V = \Delta x \Delta y$ 是单元体积， $\Delta t$ 是(拟)时间步长。这个 $V/\Delta t$ 是唯一支撑对角的机制，但它追不上 $a_N$ 的增长速度。

点光滑器的收敛系数里带着 AR#

Jacobi 迭代的误差衰减，被非对角之和与对角之比从上方压住。

\rho_{\text{point}} \le \frac{2a_E + 2a_N}{V/\Delta t + 2a_E + 2a_N}

取 $V/\Delta t = 4a_E$ ，则 $\mathrm{AR} = 1$ 时为 0.5， $\mathrm{AR} = 100$ 时为 0.9998。问题是这个上界随 $\mathrm{AR}$ 增大会贴到 1，从此不再提供任何信息。要看清真正留下来的是什么，得逐模态分析。对 Fourier 模态 $(\theta_x, \theta_y)$ ，衰减因子为

\rho(\theta_x, \theta_y) = 1 - \frac{V/\Delta t + 2a_E(1 - \cos\theta_x) + 2a_N(1 - \cos\theta_y)}{a_P}

最慢的模态是 $\theta_y$ 最小的那个，也就是壁面法向上变化最平缓的模态。当 $a_N$ 主导对角、 $a_P \approx 2a_N$ 时

\rho_{\max} \approx 1 - \frac{\pi^{2}}{2 N_y^{2}}

$N_y$ 是壁面法向的层数。这里没有 $\mathrm{AR}$ ，这不是好消息而是坏消息：再怎么压扁，系数也不会超过这个值，可这个天花板本身已经高到没用。 $N_y = 48$ 时它是 0.99786，把误差降到 $10^{-6}$ 需要大约 6,000 次。

发生了什么其实很直观。点光滑器一两次就能抹掉壁面法向的高频误差。麻烦在之后：要抹掉沿壁面方向的误差，信息必须顺着 $a_E$ 横向传递，而对角被 $a_N$ 压住，于是每一步的幅度缩小了 $\mathrm{AR}^2$ 倍。

只看对角的迭代，换一个名字会停在哪里，写在从 Jacobi 到 BiCGSTAB 的收敛速度对照里。那篇换的是求解器。这边是单元长宽比背叛了同一个求解器。

把三种处方放在同一张表上

预条件	每次迭代代价	应对各向异性	并行性	卡住的地方
Jacobi(对角)	每单元一次除法	无	完全	$\mathrm{AR} \gtrsim 10$ 就崩
ILU(0)	一次分解 + 前后代入	部分	依赖编号	每换一次区域分解性能就变
线元	线内每单元 Thomas 约 5 flop	精确消去强方向	按线独立	线之外仍然是 Jacobi

线元的交换条件很明确：只把一个强耦合方向精确求逆。代价是一组三对角带的存储，开销与 Jacobi 同一量级。因为不像 ILU(0) 那样动整个矩阵，它对 MPI 分区也不敏感。

线元是怎么搭出来的

在每个单元上量面系数比 $\sigma = a_f / a_{f,\perp}$ 。壁面法向的面上 $\sigma = \mathrm{AR}^2$ 。
把 $\sigma \ge \sigma_{\min}$ 的面标记为强面。工程上常用的默认值在 $\sigma_{\min} = 10$ 附近。
从贴壁单元出发，沿强面向上延伸出一条链。
每个单元最多保留两条强面。出现分叉时只留系数大的一侧——允许分叉，结果就不是三对角。
一旦碰到已属于其他线元的单元，就把链切断。
长度为 1 的链丢弃，把该单元交还给 Jacobi。

只在链内部重新编号，那块子矩阵就成了三对角。三对角就能用 Thomas 算法(先前向消去再回代)在 $O(n)$ 内精确求逆。线元预条件的全部内容就是这一句话。

链在哪里断开，放在迭代表后面看。旋钮是第一层厚度和增长率。

用 Python 数迭代次数#

在 $64 \times 48$ 的网格上只留下误差(右端项 $f = 0$ ，所以精确解是 $u = 0$ )，分别跑点 Jacobi 和线 Gauss-Seidel，数出误差降到 $10^{-6}$ 所需的 sweep 数。

import numpy as np
 
def cell_coefficients(dx, dy, gamma=1.0, diag_factor=4.0):
    """单个矩形单元的 FVM 扩散系数。diag_factor 以 a_E 的倍数指定 V/dt。"""
    a_e = gamma * dy / dx        # 东/西面
    a_n = gamma * dx / dy        # 北/南面
    d0 = diag_factor * a_e       # V/dt 项
    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 边界
    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):
    """三对角精确解法：先前向消去，再回代。"""
    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):
    """把一整条壁面法向的链精确求逆(线 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])   # 线外的耦合放到右端
        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}")

跑出来是这样。

    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

Δy_wall631 µmgrowth1.25log₁₀ σ_min32

What to watch: shrink Δy_wall and the cyan chain grows downward into the layer while the bars on the right cross the red threshold. Raise the growth ratio and the chain gets shorter — the mesh goes isotropic sooner, so there is less for the linelet to own. Push σ_min past a few thousand and the chains disappear entirely: the preconditioner quietly degrades back to plain Jacobi, which is the failure mode nobody notices in the log file.

减小 $\Delta y_{\text{wall}}$ ，链会往边界层深处扎。提高增长率，网格更早回到各向同性，链就变短。链很短，连一半层数都盖不住。线元便宜，是因为这条长度。

压扁 1000 倍后，Jacobi 慢了 300 倍，并如上一节所料在 6,000 次附近撞上天花板。线元反而更快了。相邻链之间的耦合按 $a_E/a_N = \mathrm{AR}^{-2}$ 变小，链与链彼此独立。最坏模态的上界 $2a_E/(V/\Delta t + 2a_E) = 1/3$ 依然成立，只是实测衰减远好于它。

只把 Jacobi 那一列拿来除。18、514、4879、5471。长宽比是 1、10、100、1000。

第一个 10 倍，迭代变成 28.6 倍。下一个 10 倍是 9.5 倍。最后一个 10 倍只有 1.12 倍，也就是 5471 除以 4879。系数比仍按 $\mathrm{AR}^2$ 变大，迭代次数却不再跟着走。 $N_y = 48$ 时模态公式 $1 - \frac{\pi^{2}}{2 N_y^{2}}$ 等于 0.99786。解 $0.99786^{n} = 10^{-6}$ ，n 在 6400 附近。表上的 5471 就在它下面。网格再压扁，这一列也贴在这里。

线元那一列反过来。8、7、4、2。衰减是 0.1839、0.1699、0.0194、0.0002。上界 $2a_E/(V/\Delta t + 2a_E) = 1/3$ 里没有 $\mathrm{AR}$ ，实测全都在它下面。AR 1000 的 0.0002 大约是 1/3 的 0.06%。链与链之间的耦合 $a_E/a_N = \mathrm{AR}^{-2}$ 熄掉以后，一次 Thomas 几乎不留下壁面法向误差。AR 为 1 时是 8 次对 18 次，赚得少。AR 为 10 时是 514 对 7，大约 73 倍。处方在这里分开。

点光滑器的上界用的是同一条对角。 $V/\Delta t = 4a_E$ 时分母是 $6a_E + 2a_N$ 。AR 为 1 时上界 0.5，表里的 rho_J 是 0.4870，几乎相同。AR 为 100 时上界 0.9998，实测 rho_J 是 0.9975。上界贴到 1 之后，模态公式的 0.99786 更诚实。上一节说的大约 6,000 次，和表里的 5471，指的是同一块天花板。

写代码时别掉进去的坑

把线铺反方向。 链必须沿系数大的方向铺，也就是壁面法向。沿壁面铺，收益为零，只多付 Thomas 的开销。网格在流向上看起来"更长"，所以很容易搞反。

把 $\sigma_{\min}$ 留在默认值。 换了网格却不重新审视阈值，链可能整体消失。预条件于是悄悄退化成 Jacobi，日志里不会有任何警告。每次运行都打印线元条数和平均长度，一行就能抓住这个问题。

让 MPI 分区把链切断。 在分区边界被切断的链会变成碎片，核数越多收敛越差。给分区器加上沿线方向的权重，或者直接约束它不要切线。

想用预条件盖住非线性问题。 如果加了线元还是发散，问题通常出在外层循环而非线性求解器。先降 CFL，把亚松弛因子压到 0.7 附近以缩小更新量，再下判断。预条件只会把给定的矩阵解得更快，不会修好一个错的矩阵。

线不沿壁面法向，Thomas 就是白花钱#

链只沿着系数大的一侧铺，也就是壁面法向。顺着流向铺，眼睛只是跟上了长边。Thomas 的开销照付，衰减却和 Jacobi 一样。

日志里留下链条数和平均长度。打出 0，说明 $\sigma_{\min}$ 已经和网格对不上。预条件这时就是对角光滑器。这个状态下只加核，断掉的链会被切得更碎。边加密网格边测，效率掉到 16.3% 的 weak scaling就是这笔账。