CFL取1要一万步，取100会爆炸 — Switched Evolution Relaxation

拿掉CFL上限之后，残差稍一下降步长就没有天花板，后段的Newton还没开始，前段就已经吹飞。

定常残差对CFL很娇气。放在1，转上一万步也还悬在上面。拉到100，第一步就把解吹飞。SER(Switched Evolution Relaxation)让残差在这两头之间自己挑步长。降了就放大，跳了就收回。

把固定CFL滑条从1推到20，琥珀色残差会跳成diverged，SER这边只留起始CFL和指数。

● fixed CFL = 3: 265 iters● SER (CFL₀=1, p=1): 194 iters

viscosity ν = 0.020fixed CFL = 3SER start CFL₀ = 1SER exponent p = 1.0SER cap CFLₘₐₓ = ∞

Residual L2 norm (normalized by initial), log scale. Dashed line = convergence tolerance 1e-9.

$\mathrm{CFL}_0=1$ 、 $p=1$ 时青色曲线在后段折下去。把 $p$ 抬到1.5会更陡，可前期残差一晃，青色也跟着晃。

这里讲的CFL数(用网格尺度和波速把时间步长无量纲化后的量),在定常求解器里已经不再关乎精度,而是同时握着收敛速度与稳定性的一个旋钮。

定常解要用假的时间去追

定常方程只有一行: $R(U) = 0$ 。其中 $R$ 把对流、扩散、源项汇成空间残差。直接求解它就要解一个非线性代数系统。于是大多数CFD程序会加上一个假的时间项。

\frac{\partial U}{\partial \tau} + R(U) = 0

这里 $\tau$ 是没有物理意义的伪时间(pseudo-time)。当 $\tau \to \infty$ 时 $\partial U / \partial \tau \to 0$ ,于是 $R(U)=0$ ,正是我们想要的定常解。舍弃时间精度,只盯定常解。这叫伪瞬态延拓(pseudo-transient continuation)。

把假的时间用隐式(backward Euler)推进并对残差线性化,每次迭代就变成这样:

\left( \frac{I}{\Delta\tau} + \frac{\partial R}{\partial U} \right) \Delta U = -R(U^n)

其中 $\partial R/\partial U$ 是通量Jacobian, $\Delta U = U^{n+1}-U^n$ 。每次迭代解一次左端的矩阵,得到 $\Delta U$ 。

小处起步,落在Newton上#

这一行里装着SER的全部直觉。把 $\Delta\tau$ 推向两个极端看看。

当 $\Delta\tau \to 0$ , $I/\Delta\tau$ 压过矩阵,于是 $\Delta U \approx -\Delta\tau\, R$ 。这跟显式推进完全一样。非常稳,但慢如乌龟。

当 $\Delta\tau \to \infty$ , $I/\Delta\tau$ 消失,只剩 $(\partial R/\partial U)\,\Delta U = -R$ 。这正是Newton法,在解附近具有二阶收敛(误差按平方缩小)。但初值太差就会发散。

所以理想策略很清楚。前期把 $\Delta\tau$ 留小以承受非线性,等残差下降、接近解时再把 $\Delta\tau$ 调大,把Newton的二阶收敛吸收过来。 难点在"何时、调多大"。SER把这个判断交给残差。

SER — 用残差的倒数放大CFL#

Mulder与van Leer在1985年提出的规则简单得惊人。用初始残差与当前残差之比放大当前CFL。

\mathrm{CFL}_n = \min\!\left( \mathrm{CFL}_0 \left( \frac{\lVert R_0 \rVert}{\lVert R_n \rVert} \right)^{\!p},\ \mathrm{CFL}_{\max} \right)

$\mathrm{CFL}_0$ 是起始CFL, $\lVert R_n \rVert$ 是第 $n$ 次残差的L2范数, $p$ 是通常约为1的指数, $\mathrm{CFL}_{\max}$ 用来封顶防止失控。残差降十倍(当 $p=1$ 时)CFL也涨十倍。伪时间步长由这个CFL和局部波速换算。

\Delta\tau_i = \mathrm{CFL}_n \cdot \frac{\Delta x}{|u_i| + 2\nu/\Delta x}

分母是对流波速 $|u_i|$ 与扩散波速 $2\nu/\Delta x$ 之和。每个单元可以不同,因此与局部时间推进(local time-stepping)自然结合。

格点数 N=80 时,间距是 1/81。 $\nu=0.02$ 时扩散波速 $2\nu/\Delta x$ 算出来是 3.24。初始猜测的最大值超不过 1,分母里扩散压过对流 $|u_i|$ 。CFL 固定在 3,则 $\Delta\tau$ 是把 1/81 除以略大于 4 的数再乘 3,不到 0.01。同一个分母放进上限 1e5, $\Delta\tau$ 就到了几百。左端的 $I/\Delta\tau$ 几乎消失。代码并不按单元使用 $|u_i|$ 。所有单元共用 $|u|$ 的最大值。激波层正中、速度接近 0 的单元,也借最快单元的波速。比正文里的局部时间推进保守一拍。

Python — 用两种策略解定常Burgers#

求解粘性Burgers的定常解 $u u_x = \nu u_{xx}$ 。边界条件 $u(0)=1$ 、 $u(1)=-1$ 在中间形成一道静止激波层。一步伪时间隐式推进用Thomas算法解三对角矩阵。一维扩散里那段回代保住无条件稳定的位置,和这里在左端加上时间项的位置是同一处。

import numpy as np
 
def steady_residual(u, nu, dx):
    """R(u) = u u_x - nu u_xx, 边界 u(0)=1, u(1)=-1."""
    um = np.concatenate(([1.0], u[:-1]))   # 左邻
    up = np.concatenate((u[1:], [-1.0]))   # 右邻
    conv = u * (up - um) / (2 * dx)
    diff = nu * (up - 2 * u + um) / dx**2
    return conv - diff
 
def thomas_solve(a, b, c, d):
    """解三对角系 (下/主/上对角 a,b,c, 右端 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):
    """一步伪时间隐式推进: (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                       # 下对角
    b = (up - um) / (2 * dx) + 2 * nu / dx**2 + 1.0 / dt  # 主对角
    c = u / (2 * dx) - nu / dx**2                        # 上对角
    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                          # 线性初值
    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")

输出如下。

fixed CFL=3    ->  588 iters
SER  p=1.0     ->   34 iters

把cap滑条往下拉,输出里的34步会在哪里重新变长,青色CFL和粉色残差一起动。

converged in 195 iterations

start CFL₀ = 1exponent p = 1.0cap CFLₘₐₓ = 2000

Cyan = CFL schedule (left log axis), pink = residual (right log axis). Both share the iteration axis.

前期残差平坦的区间里,CFL也几乎停在 $\mathrm{CFL}_0$ 。残差一开始崩塌的瞬间,CFL几乎垂直地窜上去。这就是SER的本质:解越好,允许的步子越大。 把cap降到20,CFL撞到天花板,失去末期的Newton加速,迭代次数又涨回去。

要达到同样精度的定常解,固定CFL需要几百次,SER只要几十次。差距在末期拉开。因为SER把CFL拉到上千,实际上变成了Newton迭代。

把固定CFL和SER并排放#

588除以34大约是17倍。容差是当前残差除以初始残差,小于1e-9的那一点。 $p=1$ 、起始CFL为1、上限为1e5时,残差比达到1e5的瞬间CFL就贴到天花板。那时残差比是1e-5。到容差1e-9还剩四个数量级。那段的CFL钉在1e5,矩阵已经靠在Newton一侧。固定CFL 3没有这段,所以588步全得自己走完。

调度写成一行。 $\mathrm{CFL}_n = \min(\mathrm{CFL}_0 (\lVert R_0\rVert/\lVert R_n\rVert)^p,\ \mathrm{CFL}_{\max})$ 。残差下降,步子变大,到上限停住。

现场打开SER时的陷阱#

$\mathrm{CFL}_{\max}$ 一定要设。留成无穷大,残差一瞬跳动就让CFL失控,把解吹飞。

当残差非单调(non-monotone)地晃动时, $R_0/R_n$ 会小于1,使CFL反而下降。这不是bug,而是安全阀。残差再度变差时,缩小步子才是对的。

$p$ 不是免费的加速。 $p>1$ 会在残差只降一点时就把CFL放大很多,放大风险。稳妥的起步是 $p=1$ 、 $\mathrm{CFL}_0 \approx 1$ 。

如果Jacobian不精确(近似Jacobian,或无矩阵),大CFL下Newton的二阶收敛会掉到一阶。这时降低 $\mathrm{CFL}_{\max}$ 换取稳定性更划算。

上限空着,后段的Newton先死#

$\mathrm{CFL}_{\max}$ 一旦空着,残差跳一下,CFL就顺着残差比没有尽头地涨。伪时间隐式迭代在 $\Delta\tau\to\infty$ 时变成Newton法。那二阶只在解附近出现。前期就跳进这个极限,矩阵会把解往反方向推。代码里的上限1e5是一扇门。残差降到十万分之一以前,不让步子进那一块。

每次都凑不出精确Jacobian时,同一扇门还在。大CFL的二阶已经掉成一阶,剩下的选择就是把上限再压低。用扫描近似那块矩阵、却不把它组装出来,是显式一步炸掉之后的LU-SGS另写的事。隐式步还在,矩阵不必存下来。