显式炸开时,不求逆矩阵也能解隐式 — LU-SGS

计算在凌晨3点发散了。日志末行总是 pressure residual = NaN。CFL从0.9降到0.5,再降到0.3。一步太短,定常要好几天。可压缩定常计算里常见的墙。隐式时间推进从时间步长这一侧把墙拆开。不求全矩阵逆的那条路是 LU-SGS。CFL等于5时Burgers激波还在。

把格式开关在前向Euler和后向Euler之间拨动,并把 $\mathrm{Re}(\lambda\Delta t)$ 拖到 $-2.5$。稳定区域会从小圆盘变成左半平面。

前向Euler下该点离开青色圆盘,变成UNSTABLE。后向Euler把整个左半平面填成青色,同一坐标重新稳定。

稳定域无限敞开

时间推进的稳定性由一个模型方程决定。

其中 $\lambda$ 是空间离散产生的特征值(来自波速与黏性),$u$ 是某个模态的振幅。一步之后振幅乘上放大因子 $G$。要不发散就需 $|G|\le 1$。

显式推进(前向Euler)给出

即 $|1+\lambda\Delta t|\le 1$,也就是 $\lambda\Delta t$ 必须落在以 $-1$ 为心、半径为 $1$ 的小圆内才稳定。黏性强的网格或薄单元会让 $\lambda$ 变大,很快就跑出这个圆。

隐式推进(后向Euler)正相反。

对所有 $\mathrm{Re}(\lambda)<0$ 的物理衰减模态,$|G|<1$。无论 $\Delta t$ 多大都稳定。稳定性约束消失,时间步长只需小到精度所要求的程度即可。

显式稳定圆盘在实轴上以0和 $-2$ 为两端。圆心 $-1$,半径1。$\mathrm{Re}(\lambda\Delta t)=-2.5$ 离圆心1.5,落在圆外,所以标记变成UNSTABLE。后向Euler的放大因子是 $G=1/(1-\lambda\Delta t)$。$\mathrm{Re}(\lambda)<0$ 时分母大于1,$|G|$ 在 $\Delta t$ 增大后仍不超过1。凌晨日志把CFL从0.9、0.5降到0.3,是想把 $\lambda\Delta t$ 关进这个圆盘。黏性网格上 $\lambda$ 先变大。圆盘在精度要求更小步长之前就会破。

求全矩阵逆的代价付不起

天下没有白吃的午餐。隐式要求未来时刻的残差。

$R_i$ 是单元 $i$ 的通量散度(残差),$U$ 是守恒变量。把 $R(U^{n+1})$ 在当前时刻附近线性化,

便得到关于 $\Delta U = U^{n+1}-U^n$ 的线性系统。

麻烦在左端矩阵 $A$。若有五千万单元,$A$ 就是五千万维。直接求逆在内存和时间上都不可能。像GMRES这样的Krylov求解器也行,但要显式存储矩阵并需要预处理。1990年代的可压缩分析者找到了更便宜的路:不存矩阵、不做乘法,用两次扫掠来近似。

LU-SGS — 拆成 D·L·U,扫两遍#

把 $A$ 分解为块下三角 $L$、块对角 $D$、块上三角 $U$。

LU-SGS(下-上对称Gauss-Seidel)用如下乘积近似这个 $A$。

它丢掉了 $LU$ 乘积中产生的 $L D^{-1} U$ 项。当 $\Delta U$ 较小时可忽略。求解就归结为两个三角系统的顺序回代。

前向扫掠 —— 下三角,从前到后:

后向扫掠 —— 上三角,从后到前:

用什么填 $D$、$L$、$U$,决定近似好不好。Yoon与Jameson(1988)用谱半径 $\rho$ 近似非线性通量雅可比。把一阶迎风拆成 $A^{\pm} = \tfrac{1}{2}(a \pm |a|)$,对角就是 $D_i = I/\Delta t + |a_i|/\Delta x$,非对角只剩拆分后的面雅可比。无需构建整个矩阵,每个单元几个标量(或一个小块)就够了。

近似乘积是 $(D+L)\,D^{-1}\,(D+U)$。展开后 $A$ 旁边多出 $L D^{-1} U$,这一项被丢掉。因为 $D_i = I/\Delta t + |a_i|/\Delta x$,$\Delta t$ 小时 $I/\Delta t$ 把对角撑大,丢掉的项相对变小。CFL升到5,$\Delta t$ 变大,$I/\Delta t$ 变小,近似误差跟着变大。显式在两三步里就出现NaN的那段,输出仍用49步走完,$u$ 落在0.012到0.988之间。格子是 $n=200$,算到 $t=0.22$,$\Delta x=1/200$。与存矩阵的 Krylov求解器收敛比较 不同,这里用一次对角除法和两次扫掠占住那个位置。

Python —— 用一次分解推进隐式Burgers#

玩具问题取1D Burgers方程的激波形成。光滑正弦波随时间卷成激波。显式必须把CFL压在1以下。LU-SGS隐式在CFL 5下推进也不会死。

import numpy as np
 
def upwind_flux_residual(u, dx):
    """用一阶迎风通量计算 d(u^2/2)/dx 的残差 R_i(周期边界)。"""
    f = 0.5 * u * u
    a = u                       # 局部波速 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):
    """一个时间步:用一次LU-SGS应用推进隐式系统。"""
    n = u.size
    a  = u
    ap = 0.5 * (a + np.abs(a))         # A+ (下三角贡献)
    am = 0.5 * (a - np.abs(a))         # A- (上三角贡献)
    D  = 1.0 / dt + np.abs(a) / dx     # 对角块
    R  = upwind_flux_residual(u, dx)
 
    dus = np.zeros(n)                  # 前向扫掠: (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()                   # 后向扫掠: (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)   # 光滑初始条件
    t, steps = 0.0, 0
    while t < tmax:
        dt = cfl * dx / np.max(np.abs(u))    # 允许 CFL>1(隐式)
        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]

同一问题用显式前向Euler在CFL 5下跑,两三步就出NaN。隐式版本以49个稳定步抵达激波。没有逐面的重型Riemann求解器,一次对角除法加两次扫掠就完成了全部工作。

一遍遍扫掠,残差就崩塌

LU-SGS应用一次得到近似解。要解得更准,就在同一系统上重复扫掠。此时收敛速度取决于矩阵的对角占优程度。在黏性主导的情形——黏性雅可比变成1D拉普拉斯算子——对角为 $1+2d$,非对角为 $-d$。其中 $d=\nu\Delta t/\Delta x^2$ 是扩散数(黏性扩散的无量纲时间步长)。

$d=\nu\Delta t/\Delta x^2$ 等于1时,对角是 $1+2d=3$,非对角是 $-d=-1$,左右各一个。对角大于两个非对角之和,几次扫掠残差就掉十个数量级。$d=40$ 时对角是81,非对角是两个 $-40$,行和变成 $81-80=1$,几乎抵消。曲线因此躺平。用大 $\Delta t$ 买来的稳定,和内层松弛需要的扫掠次数,记在同一本账上。

把 $d$ 滑条从1附近拉到40,同样的扫掠次数下残差曲线的斜率会变。

把 $d$ 放在1附近,残差几次扫掠就掉十个数量级。把 $d$ 加到40(巨大的隐式 $\Delta t$),同一条曲线就躺平了。你靠加大步长买到了稳定,但矩阵更硬,内层松弛需要更多扫掠。稳定性与内层收敛不会免费同时到来。

对角里漏掉黏性半径,扫掠就死

对角 $D$ 不加 $\nu/\Delta x^2$,黏性大的单元里 $D$ 会变小。谱半径只剩对流的 $|a|/\Delta x$。扫掠马上发散。黏性雅可比若是一维拉氏算子,对角占优是 $1+2d$ 和 $-d$ 的平衡,抽掉黏性项,平衡的左边就空了。非线性强时不要一次加上整个 $\Delta U$,用 $U^{n+1}=U^n+\omega\,\Delta U$,$\omega\in[0.5,0.8]$。激波附近的振荡会减小。

方向也是一个决定。只重复前向扫掠,信息只往编号增大的方向流。前后交替,$L$ 和 $U$ 才保持对称。MPI分区面上,两次扫掠之间必须插一次通信,邻块的 $\Delta U$ 才会过来。单元编号怎么排,改的就是这次扫的顺序。用RCM打乱编号后同一套LU-SGS变慢的记录是这个决定的另一面。