陽解法が破綻するとき、行列を解かずに陰解法を解く — LU-SGS

解析は午前3時に発散しました。ログの末尾はいつも pressure residual = NaN です。CFLを0.9から0.5へ、さらに0.3へ下げました。1ステップが短すぎて、定常まで何日もかかります。圧縮性の定常計算でよく当たる壁です。陰的時間進行はその壁を時間刻みの側から壊します。全行列を反転しない計算が LU-SGS です。CFL 5でもBurgersは残ります。

スキームのトグルを前進Eulerと後退Eulerのあいだで切り替え、$\mathrm{Re}(\lambda\Delta t)$ を $-2.5$ まで引いてみてください。安定領域が小さな円盤から左半平面へ変わります。

前進Eulerではその点がシアンの円の外へ出てUNSTABLEに反転します。後退Eulerにすると左半平面全体がシアンになり、同じ座標が再び安定です。

安定領域が無限に開く

時間進行の安定性は、一つのモデル方程式で決まります。

ここで $\lambda$ は空間離散化が生む固有値(波速や粘性から来ます)、$u$ は一つのモードの振幅です。1ステップ後の振幅は増幅率 $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$ です。セルが5千万個なら $A$ は5千万次元です。直接の逆行列はメモリも時間も不可能です。GMRESのようなKrylovソルバーも有効ですが、行列を明示的に格納し前処理が必要です。1990年代の圧縮性解析者はより安い道を見つけました。行列を格納せず、掛け算なしで、二度の掃き出しで近似するのです。

LU-SGS — D・L・Uに分けて二度掃く#

$A$ をブロック下三角 $L$、ブロック対角 $D$、ブロック上三角 $U$ に分解します。

LU-SGS(Lower-Upper Symmetric Gauss-Seidel)は、この $A$ を次の積で近似します。

$LU$ 積から生じる $L D^{-1} U$ 項を捨てたものです。$\Delta U$ が小さければ無視できます。解法は三角システム二つの順次代入で終わります。

前進スイープ — 下三角を前から後ろへ:

後進スイープ — 上三角を後ろから前へ:

$D$、$L$、$U$ を何で埋めるかが近似の質を分けます。YoonとJameson(1988)は非線形フラックスヤコビアンをスペクトル半径 $\rho$ で近似しました。1次風上を $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$ は縮みます。近似誤差が大きくなる場所です。それでも陽解法が2〜3ステップで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):
    """1次風上フラックスで 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):
    """1時間ステップ:陰的システムをLU-SGS 1回適用で進める。"""
    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に回すと、2〜3ステップで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$ が両隣です。対角が二つの非対角の和より大きく、数回のスイープで残差が10桁落ちます。$d=40$ なら対角は81、非対角は $-40$ が二つで、行の和は $81-80=1$ とほぼ打ち消し合います。曲線が寝る理由です。大きな $\Delta t$ で買った安定と、内側緩和が要求するスイープ数は同じ帳簿に並びます。

$d$ のスライダーを1付近から40まで上げると、同じスイープ数でも残差曲線の傾きが変わります。

$d$ を1付近に置くと、残差は数回のスイープで10桁落ちます。$d$ を40まで上げると(巨大な陰的 $\Delta t$)、同じ曲線が寝てしまいます。時間刻みを大きくして安定性を買いましたが、行列が硬くなり、内側の緩和がより多くのスイープを要求します。安定性と内側収束はただで同時には来ません。

対角から粘性半径を抜くとスイープが死ぬ

対角 $D$ に $\nu/\Delta x^2$ を入れないと、粘性の大きいセルで $D$ が縮みます。スペクトル半径を対流の $|a|/\Delta x$ だけで埋めた状態です。スイープはすぐに発散します。粘性ヤコビアンが1Dラプラシアンなら、対角優位は $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が遅くなる記録が、その決定の反対側です。