[論文レビュー] GMRESが1000回回っても止まらなかった — 全速度 melt pool のための物理ベース・ブロック前処理

固有値が1に集まるのは厳密なSchur補元のときだけです。AMGでSを粗く解くと外側の反復が再び開きます。

レーザーで金属を溶かす melt pool を陰的に回しました。最初のタイムステップで GMRES が 1000 回回っても、残差は 1e-2 に張り付いたままです。Mach を 1e-3 に下げると音速は流速の 1000 倍になります。ヤコビアンはその比で病みました。Weston・Nourgaliev・Delplanque・Barker（2019）は、速度–圧力 Schur 補元を前処理にしてこの停滞を解きました。今日はその前処理を 2 場の toy に移し、Mach を変えても反復数がどこに張り付くかを記録します。

Mach スライダーを下げると、上のトラックの音響波だけが $c=1/M$ で速くなり、下の対流トレーサは $u=1$ のままです。

Mach number M = 3.2e-2

explicit acoustic CFL forces Δt ∝ M·Δx — the smaller M, the stiffer the coupled system.

$M=10^{-2}$ では音響波がトラックを 100 周する間に、トレーサは一目盛りしか動きません。 $M=10^{-3}$ では比が 1000:1 です。残差が 1e-2 に張り付いた夜と同じ隙間です。

この論文が立つ場所

著者: B. Weston, R. Nourgaliev, J.-P. Delplanque, A. T. Barker (LLNL / UC Davis)
ジャーナル: Journal of Computational Physics, 397, 108847 (2019)
タイトル: Preconditioning a Newton-Krylov solver for all-speed melt pool flow physics
一言: 相変化を含む全速度（all-speed）圧縮性 Navier–Stokes を完全陰的 JFNK で解き、原始変数ブロック Schur 補元前処理で低 Mach の剛性を抑える。

既存の前処理 — element-block SOR、monolithic AMG、block Gauss–Seidel — はいずれも、タイムステップや Mach が極端になると反復数が爆発しました。論文が出したのは vP-vT Schur 補元です。物理結合のブロックの形を、前処理行列へそのまま移したものです。

M がゼロに向かうとヤコビアンが病む理由#

全速度とは、一つのコードで非圧縮（ $M\to 0$ ）から超音速までを解くことです。問題は低 Mach 極限です。音響波速度 $c$ と流速 $u$ の比が $c/u = 1/M$ と広がります。 $M=10^{-3}$ なら、音は流体より 1000 倍速いのです。

完全陰的にすると、この二つのスケールが一つのヤコビアンに混ざります。圧力–速度結合ブロックが $1/M^2$ で大きくなり、行列の条件数が爆発します。Krylov ソルバーは固有値が広く散らばった行列の上で這うように進みます。これが GMRES の止まらなかった理由です。

コードの四点に $M$ を置くと $\beta=1/M^2$ がすぐ出ます。 $M=1$ なら $\beta=1$ 。 $M=10^{-1}$ なら $\beta=10^{2}$ 。 $M=10^{-2}$ なら $\beta=10^{4}$ 。 $M=10^{-3}$ なら $\beta=10^{6}$ です。圧力–速度ブロックの一マスが百万倍になります。格子数は $n=64$ 、間隔は $dx=1/64$ で四つの場合とも同じです。行列の大きさも $2n=128$ のままです。変わるのは $\beta$ を掛けたスペクトルだけです。

音響 CFL で時間刻みを絞ると $\Delta t \propto M\Delta x$ です。 $M$ を 1 から $10^{-3}$ に下げると、許される刻みは 1000 分の 1 になります。完全陰的はその制限を外し、代わりに条件数を $\beta$ へ渡します。同じ剛性を前処理が逆に増やした場所は、音速 CFL 10 のブロック前処理に別に書いてあります。

JFNK（Jacobian-Free Newton–Krylov）はヤコビアンを陽に作らず、行列–ベクトル積だけを Fréchet 差分で近似します。

J\boldsymbol{\kappa} \approx \frac{F(x + h\boldsymbol{\kappa}) - F(x)}{h}

$F$ は非線形残差、 $\boldsymbol{\kappa}$ は Krylov ベクトル、 $h$ は小さな有限数です。行列を保存しないのでメモリは減ります。前処理が無いと、低 Mach ではその差分が散らばった固有値の上をそのまま歩きます。

保存変数を捨て原始変数 $(P,v,T)$ へ#

最初の一手は変数の置換です。保存変数 $U=(\rho, \rho v, E)$ の代わりに原始変数 $W=(P, v, T)$ を解きます。低 Mach では原始変数のほうがはるかに良く条件付けされるからです。

\frac{\partial F(x)}{\partial U}\frac{\partial U}{\partial W}\,\delta W = -F(x)

$\partial U/\partial W$ は変数変換ヤコビアンです。残差 $F$ は依然として保存則を満たすように書かれているので、変数を変えても質量・運動量・エネルギーの保存は壊れません。ここで自由度を場ごとに束ね、 $3\times3$ ブロック行列に整列します。

M = \begin{bmatrix} M_{vv} & M_{vP} & M_{vT} \\ M_{Pv} & M_{PP} & M_{PT} \\ M_{Tv} & M_{TP} & M_{TT} \end{bmatrix}

各ブロックは一つの原始場（速度・圧力・温度）に対応します。この整列が次の分解の足場です。

ブロック LU と速度–圧力 Schur 補元#

圧力–温度結合（ $M_{PT}, M_{TP}$ ）は弱いです。両者を捨てるとブロック LU 分解が綺麗になり、残る核は速度–圧力 Schur 補元です。

S_{vP} = M_{PP} - M_{Pv}\,M_{vv}^{-1}\,M_{vP}

$S_{vP}$ は速度を消去したのちに圧力だけを残した有効方程式です。SIMPLE 系の圧力補正方程式と同じ骨格です。前処理 $P$ をこのブロック上三角の形に取ると、 $P^{-1}A$ の固有値が 1 の近くに集まります。厳密な $S$ を使えば $A=LU$ の $U$ がそのまま $P$ となり、 $P^{-1}A$ は冪零（unipotent）下三角と相似で、すべての固有値が厳密に 1 です。GMRES は二、三回で終わります。

その二、三回は停止条件とも噛み合います。下のコードの許容残差は rtol=1e-8、反復の上限は maxiter=500、restart の長さは 2n=128 です。固有値が一点に集まれば Krylov 部分空間は二歩で閉じます。集まらなければ、先に 500 回の上限へ当たります。

Python — Mach を下げながら GMRES 反復数を数える#

温度を外し、2 場 $(v,P)$ の toy で核だけを再現します。ブロック系は $A=\begin{bmatrix} I & \beta G \\ D_v & I\end{bmatrix}$ 、 $\beta=1/M^2$ 。前処理は厳密な Schur 補元 $S=I-\beta D_v G$ です。

import numpy as np
from scipy.sparse import diags, eye, bmat
from scipy.sparse.linalg import gmres, LinearOperator, splu
 
n = 64
dx = 1.0 / n
 
def fwd_diff():
    # 周期境界の前進差分 (圧力勾配演算子)
    D = diags([-1.0, 1.0], [0, 1], shape=(n, n)).tolil()
    D[-1, 0] = 1.0
    return (D / dx).tocsr()
 
G = fwd_diff()
Dv = (-G.T).tocsr()          # 発散 = -G^T  ->  Dv*G ~ ラプラシアン
 
def all_speed_system(M):
    beta = 1.0 / M**2        # 低 Mach 剛性係数
    A = bmat([[eye(n), beta * G], [Dv, eye(n)]]).tocsr()
    return A, beta
 
def schur_preconditioner(beta):
    # 速度-圧力 Schur 補元 S = M_PP - M_Pv M_vv^{-1} M_vP = I - beta Dv G
    S = (eye(n) - beta * (Dv @ G)).tocsc()
    luS = splu(S)
    def apply(r):
        ru, rp = r[:n], r[n:]
        zp = luS.solve(rp)               # S z_p = r_p
        zu = ru - beta * (G @ zp)        # z_u = r_u - beta G z_p
        return np.concatenate([zu, zp])
    return LinearOperator((2 * n, 2 * n), matvec=apply)
 
def count_iters(M, use_prec):
    A, beta = all_speed_system(M)
    b = np.random.default_rng(0).standard_normal(2 * n)
    P = schur_preconditioner(beta) if use_prec else None
    it = {"k": 0}
    x, info = gmres(A, b, M=P, rtol=1e-8, maxiter=500, restart=2 * n,
                    callback=lambda xk: it.__setitem__("k", it["k"] + 1))
    return it["k"]
 
print(f"{'Mach':>8}  {'no-prec':>8}  {'Schur':>6}")
for M in [1.0, 1e-1, 1e-2, 1e-3]:
    print(f"{M:>8.0e}  {count_iters(M, False):>8d}  {count_iters(M, True):>6d}")

前処理なしでは $M$ が小さくなるほど反復数が数十・数百へと跳ね上がり、 $M=10^{-3}$ では 500 回の上限の中に入りません。Schur 側は $S=I-\beta D_v G$ を LU で一度解いておき、GMRES が二、三回で残差を 1e-11 まで落とします。 $M=1$ （ $\beta=1$ ）と $M=10^{-3}$ （ $\beta=10^{6}$ ）で Schur の反復数が同じなのは、固有値が 1 に集まった数値です。

インタラクティブ — 音響と対流、開いていく二つの速度

トラックはすでに上で回しました。 $M=10^{-2}$ の 100:1 は $c/u=1/M$ そのものです。陽的 CFL が $\Delta t \propto M\Delta x$ で刻みを踏みつぶす場所も、その画面の速度差と同じです。ヤコビアンが受け取るのは、その比を二乗した $\beta$ です。 $M=10^{-2}$ で $\beta=10^{4}$ 、 $M=10^{-3}$ で $\beta=10^{6}$ です。

インタラクティブ — 前処理が残差を崩す場所

ここからは実際の GMRES 残差です。Mach スライダーを動かすと、前処理の有無の二曲線が描き直されます。

Mach number M = 1.0e-2 (β = 1/M² = 1.0e+4)

unpreconditioned: 0 iters to 1.0e+0 · Schur: 0 iters to 1.0e+0

橙色（前処理なし）は $M$ を下げるほど平らになって止まります。シアン（Schur 前処理）は $M$ がどこでも二、三回で 1e-11 まで落ちます。二曲線の差は $M=10^{-3}$ で最も大きくなります。 $\beta=10^{6}$ のマスを、物理結合を写した前処理がまとめて吸収した結果です。

自分で書いて疑わしかった三つ

厳密な $S$ が無料なのは toy だけです。 ここでは $S=I-\beta D_v G$ を LU で厳密に解き、二回で収束しました。実戦では $S_{vP}=M_{PP}-M_{Pv}M_{vv}^{-1}M_{vP}$ の $M_{vv}^{-1}$ は密な逆行列で、作ることができません。論文は $S$ を近似で解き（AMG 一回の V-cycle など）、その近似の質が外側の反復数を左右します。反復数が Mach に張り付かないコストはここに隠れています。

結合をどこまで捨てるか。 論文は $M_{PT}, M_{TP}$ を捨てても頑健性は崩れないと言います。これは melt pool の物性に合わせた仮定です。強い圧縮性熱結合（爆燃・燃焼）では、同じ切り捨てが収束を止めることがあります。前処理を単純にするかは物理が決めます。

JFNK の Fréchet 差分 $h$ 。 $h$ が小さすぎれば丸め誤差が、大きすぎれば打ち切り誤差が行列–ベクトル積を汚します。前処理が良くても、この $h$ が外れると Newton が停滞します。OpenFOAM の分離型（SIMPLE）から完全結合の陰的へ移るとき、最初に出会う場所です。

結び — 前処理はタダではない

GMRES が 1000 回、残差 1e-2 で止まった夜は、ヤコビアンをより強く解いて終わりませんでした。結合構造そのものを写した前処理が要りました。Schur 補元は速度を消した圧力の方程式であり、それを厳密に解けば固有値が 1 に集まります。実戦の $S_{vP}$ の中の $M_{vv}^{-1}$ は密なので、論文は AMG 一回の V-cycle で近似します。その近似が粗いと、外側の反復が二、三回から再び開きます。

ヤコビアンを保存せず行列–ベクトル積だけを使う分かれ道は、演算子オーバーロードの自動微分と接しています。JFNK はその積を Fréchet 差分 $h$ で置き換えます。 $h$ が外れると、前処理が良くても Newton はその場で止まります。同じ toy で $S$ を AMG により粗く解いたとき、外側の反復が二、三回からどれだけ開くかが次の測定です。