CFLを1にすると一万回、100にすると爆発する — Switched Evolution Relaxation

定常の残差はCFLに敏感です。1のままなら一万回回してもなかなか落ちません。100まで上げると、最初の一歩で解が飛びます。その間を残差に選ばせるのがSER(Switched Evolution Relaxation)です。下がれば大きくし、跳ねればまた縮めます。

fixed CFLのスライダーを1から20へ動かすと琥珀色の残差はdivergedへ飛び、SER側は開始CFLと指数だけが動きます。

$\mathrm{CFL}_0=1$、$p=1$ ならシアンの曲線は終盤に折れて下がります。$p$ を1.5に上げるとさらに急になりますが、序盤の残差が揺れた瞬間、シアンも一緒に揺れます。

ここで扱うCFL数(時間刻みを格子幅と波速で無次元化した量)は、定常ソルバーではもはや精度ではなく、収束速度と安定性を同時に握るつまみです。

定常解は偽の時間で追いかける

定常方程式は $R(U) = 0$ の一行です。$R$ は対流・拡散・ソースをまとめた空間残差です。これを直接解くには非線形代数系を解く必要があります。代わりに、ほとんどのCFDコードは偽の時間項を付け足します。

ここで $\tau$ は物理的意味を持たない擬似時間(pseudo-time)です。$\tau \to \infty$ で $\partial U / \partial \tau \to 0$ になれば $R(U)=0$、つまり求めていた定常解です。時間精度は捨て、定常解だけを狙います。これを擬似時間連続法(pseudo-transient continuation)と呼びます。

偽の時間を陰的(backward Euler)に前進させ、残差を線形化すると、各反復はこうなります。

ここで $\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反復で、解の近くで2次収束(誤差が二乗で縮む)します。ただし初期推測が悪いと発散します。

ですから理想の戦略は明確です。序盤は $\Delta\tau$ を小さく保って非線形性に耐え、残差が下がって解に近づいたら $\Delta\tau$ を大きくしてNewtonの2次収束を吸い取る。 問題は「いつ、どれだけ大きくするか」です。SERはその判断を残差に委ねます。

SER — 残差の逆数でCFLを大きくする#

Mulderとvan Leerが1985年に提案した規則は驚くほど単純です。現在のCFLを、初期残差と現在残差の比で大きくします。

$\mathrm{CFL}_0$ は開始CFL、$\lVert R_n \rVert$ は $n$ 番目の残差のL2ノルム、$p$ は通常1前後の指数、$\mathrm{CFL}_{\max}$ は暴走を防ぐ上限です。残差が10倍下がれば($p=1$ のとき)CFLも10倍大きくなります。擬似時間刻みはこのCFLと局所波速から換算します。

分母は対流波速 $|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|$ を使いません。速度の絶対値の最大を、すべてのセルが分け合います。衝撃層の真ん中、速度が 0 に近いセルも、いちばん速いセルの波速を借ります。本文の局所時間前進より一拍だけ保守的です。

Python — 定常Burgersを二つの戦略で解く#

粘性Burgersの定常解 $u u_x = \nu u_{xx}$ を解きます。境界条件は $u(0)=1$、$u(1)=-1$ で、中央に静止衝撃層ができます。擬似時間の陰的前進一歩は、三重対角行列をThomasアルゴリズムで解きます。1D拡散でその後退代入が無条件安定を作る場所と、ここで左辺に時間項を足す場所は同じです。

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とピンクの残差が一緒に動きます。

序盤の残差が平らな区間では、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が逆に下がることがあります。これはバグではなく安全装置です。残差が再び悪化したら、一歩を縮めるのが正しいです。

$p$ はタダの加速ではありません。$p>1$ は残差が少し下がっただけでCFLを大きく上げ、危険を高めます。堅実な出発は $p=1$、$\mathrm{CFL}_0 \approx 1$ です。

Jacobianが不正確(近似Jacobian、または行列フリー)だと、大きなCFLでNewtonの2次収束が1次に落ちます。このときは $\mathrm{CFL}_{\max}$ を下げて安定性を買うほうが良いです。

上限を空けると後半のNewtonが先に死ぬ#

$\mathrm{CFL}_{\max}$ を空けると、残差が一度跳ねた瞬間にCFLが残差比を際限なく追います。擬似時間の陰的反復は $\Delta\tau\to\infty$ でNewtonになります。その2次は解の近くでしか出ません。序盤にその極限へ飛ぶと、行列が解を押し上げます。コードの上限1e5は扉です。残差が十万倍下がったあとでなければ、その領域には入れません。

行列を毎回正確に作れないときも、同じ扉が要ります。大きなCFLの2次が1次に落ちたら、残る選択は上限をさらに下げることです。その行列を組み立てずにスイープで近似する道は、陽解法が飛んだあとのLU-SGSに別にあります。陰的な一歩は残り、行列は保存しなくて済みます。