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

定常解析でCFL数を1にすると、残差は一万回回しても落ちません。100まで上げると、最初の一歩で解が爆発します。どちらもよく見る光景です。この記事では、その両極端を自動で行き来する手法であるSER(Switched Evolution Relaxation)を扱います。残差が下がるぶんだけCFLを自分で大きくし、序盤は安全に、終盤はNewtonのように速く着地する方法を、小さな定常Burgersソルバーで実際に触ってみます。

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

ですから理想の戦略は明確です。序盤は $\Delta\tau$ を小さく保って非線形性に耐え、残差が下がって解に近づいたら $\Delta\tau$ を大きくしてNewtonの2次収束を吸い取る。 問題は「いつ、どれだけ大きくするか」です。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}$ は暴走を防ぐ上限です。残差が10倍下がれば( $p=1$ のとき)CFLも10倍大きくなります。擬似時間刻みはこの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)と自然に組み合わさります。

下のシミュレーションでパラメータを実際に操作してみましょう。同じ定常Burgers問題に固定CFLとSERを並べて回した残差履歴です。

● 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.

fixed CFLを1に下げると、琥珀色の曲線はなかなか底に届きません。20まで上げると発散(diverged)に跳ねます。一方SER(シアン)は $\mathrm{CFL}_0=1$ で安全に出発しても、終盤に急激に折れて下がります。指数 $p$ を1から1.5に上げるとさらに積極的に速くなりますが、上げすぎるとSERでさえ序盤で揺れます。

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

同じ精度の定常解に到達するのに、固定CFLは数百回、SERは数十回で十分です。差は終盤で開きます。SERがCFLを数千まで引き上げ、事実上Newton反復に変身するからです。

固定CFLとSERを並べると#

SERがどうやってCFLを引き上げるか、スケジュールそのものを見ましょう。下のグラフで操作してみてください。シアンが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加速を失い、反復数が再び増えます。

現場で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}$ を下げて安定性を買うほうが良い。

覚えておくこと#

定常ソルバーでCFLは精度ではなく、収束速度と安定性のつまみです。小さければ亀、大きければ爆発。
SERは $\mathrm{CFL}_n = \min(\mathrm{CFL}_0 (\lVert R_0\rVert/\lVert R_n\rVert)^p,\ \mathrm{CFL}_{\max})$ で、残差が下がれば自動的に一歩を大きくします。
擬似時間の陰的反復は $\Delta\tau\to\infty$ でNewtonになります。SERはその極限へ滑らかに連れて行くスケジュールです。