CFD線形ソルバーの選び方：JacobiからBiCGSTABまでの収束速度比較

SPDが保たれるあいだはCGがN回で終わります。対称性が崩れた瞬間、席はBiCGSTABへ移ります。

1952年の学会で、LanczosとHestenes、Stiefelは別々に書いた同じ算法を持って顔を合わせました。対称正定値（SPD）行列を $N$ 回で終える反復です。名前だけが後から共役勾配法（CG）になりました。70年たっても、圧力方程式を解く場所はその部分空間の中にあります。

格子 N を 20 から 80 へ動かすと四本の残差曲線が描き直され、最大反復の上限を上げると床に着く順番が変わります。

N=40 付近では CG が折れて降り、Jacobi は対数軸にほぼ寝ています。N を 80 まで押すと間隔はさらに開きます。BiCGSTAB のギザギザは、残差がなめらかには死なない印です。

なぜCFDでは線形システムを反復法で解くのか#

ナビエ・ストークス方程式を離散化すると、各ステップで $A\mathbf{x}=\mathbf{b}$ が落ちてきます。 $A$ は係数行列、 $\mathbf{x}$ はセル中心の未知数、 $\mathbf{b}$ は右辺です。100万セルなら $A$ は $10^6 \times 10^6$ の疎行列です。LU分解（下三角と上三角の積に割る直接法）は、この大きさでメモリと時間を同時に燃やします。OpenFOAM、SU2、Gerris が反復法に行く理由です。

初期推定 $\mathbf{x}_0$ から、残差 $\mathbf{r} = \mathbf{b} - A\mathbf{x}$ が許容値を下回るまで $\mathbf{x}$ を直します。「少しずつ」の定義がソルバーの名前です。

Krylov部分空間というアイデア#

初期残差 $\mathbf{r}_0$ に $A$ を繰り返し掛けたベクトルが、一つの空間を張ります。

\mathcal{K}_m(A, \mathbf{r}_0) = \text{span}\{\mathbf{r}_0,\, A\mathbf{r}_0,\, A^2 \mathbf{r}_0,\, \dots,\, A^{m-1} \mathbf{r}_0\}

$\mathcal{K}_m$ が Krylov 部分空間、 $m$ は次元、 $\text{span}$ は線形結合で埋めた集合です。更新は $\mathbf{x}_m \in \mathbf{x}_0 + \mathcal{K}_m$ の中で選び、残った残差 $\mathbf{r}_m$ には最適条件を満たさせます。

演算は行列ベクトル積で終わります。疎パターンは壊れません。CG はその中で A ノルム誤差を最小化し、GMRES は 2 ノルム残差を最小化します。後者は Arnoldi がベクトルを積み上げる部分空間の中で起きます。

4つのソルバーの特性比較#

ソルバー	要求される行列	メモリ	長所	短所
Jacobi	対角優位を推奨	極めて低い	並列化が容易、実装が10行程度	遅い（スペクトル半径 $1-O(1/N^2)$ ）
Gauss-Seidel	任意	極めて低い	Jacobiの約2倍速い	本質的に逐次処理、GPUに不向き
CG	SPD必須	ベクトル4本	SPDに対して最適、3項漸化式で軽量	非対称行列には不可
BiCGSTAB	任意	ベクトル7本	非対称でもOK、転置行列が不要	収束が不安定な場合がある、break-downの危険

圧力ポアソン方程式 $-\nabla^2 p = f$ は SPD なので CG が合います。移流が混ざった運動量は非対称なので、BiCGSTAB か GMRES 側です。OpenFOAM の PCG/PBiCGStab、SU2 の BCGSTAB がこの分担です。

Pythonで見る収束速度#

1次元ポアソン $-u'' = 1$ 、境界 $u(0)=u(1)=0$ を有限差分すると三重対角になります。 $N=60$ で Jacobi と CG を同じ許容値 $10^{-10}$ で回しました。

import numpy as np
 
def build_poisson_1d(n):
    A = np.zeros((n, n))
    for i in range(n):
        A[i, i] = 2.0
        if i > 0:     A[i, i - 1] = -1.0
        if i < n - 1: A[i, i + 1] = -1.0
    return A, np.ones(n)
 
def jacobi_sweep(A, b, tol=1e-10, max_iter=2000):
    D  = np.diag(A)
    LU = A - np.diag(D)
    x  = np.zeros_like(b)
    b0 = np.linalg.norm(b)
    history = []
    for _ in range(max_iter):
        x = (b - LU @ x) / D
        res = np.linalg.norm(b - A @ x) / b0
        history.append(res)
        if res < tol:
            break
    return x, history
 
def cg_sweep(A, b, tol=1e-10, max_iter=2000):
    x  = np.zeros_like(b)
    r  = b - A @ x
    p  = r.copy()
    rr = r @ r
    b0 = np.linalg.norm(b)
    history = [np.sqrt(rr) / b0]
    for _ in range(max_iter):
        Ap    = A @ p
        alpha = rr / (p @ Ap)
        x    += alpha * p
        r    -= alpha * Ap
        rr_new = r @ r
        history.append(np.sqrt(rr_new) / b0)
        if history[-1] < tol:
            break
        p  = r + (rr_new / rr) * p
        rr = rr_new
    return x, history
 
A, b = build_poisson_1d(n=60)
_, hist_j  = jacobi_sweep(A, b)
_, hist_cg = cg_sweep(A, b)
print(f"Jacobi: {len(hist_j):4d}回, 最終残差 {hist_j[-1]:.2e}")
print(f"CG:     {len(hist_cg):4d}回, 最終残差 {hist_cg[-1]:.2e}")

出力はこう印字されます。

Jacobi: 2000回, 最終残差 2.11e-02
CG:       60回, 最終残差 1.87e-13

対数の条件数スライダーを上げると、同じ反復予算の中で四本の曲線がどこで止まるかが変わります。反復の上限は、その曲線をどこで切るかを決めます。

condition number κ = 10^3.0iters80

조건수 κ를 키우면 (그림에서 오른쪽으로) Jacobi는 거의 멈추고 CG/GMRES는 √κ 한계를 따른다. BiCGSTAB는 비대칭에 강하지만 잔차가 출렁인다.

既定の $\kappa=10^3$ 、反復 80 付近では Jacobi（赤）がほぼ水平で、CG（シアン）はより急です。 $\kappa$ を $10^6$ 側へ押すと、Jacobi は反復を 200 まで与えても床が見えません。上の出力の $N=60$ は、その図の一点です。Jacobi は 2000 回で 2.11e-2、CG は 60 回で 1.87e-13 です。

出力の数字を一歩あたりに割り戻すと、様子が見えます。初期残差を 1 と置くと、2000 回後の 2.11e-2 は平均の減衰 0.998 です。0.0211 の 2000 乗根です。表のスペクトル半径 $1-O(1/N^2)$ と同じ桁です。 $N=60$ では $1/N^2$ は 1/3600、約 0.00028 です。0.998 は、その次数に定数が数倍ついた値として読めます。

この減衰が続くなら、 $10^{-10}$ までは約 1.2 万回です。すでに 2.11e-2 に着いた後でも、約 1 万回がまだ要ります。 $N$ が 2 倍になると $1/N^2$ は 4 分の 1 になり、同じ残差に必要な回数も 4 倍になります。格子を増やす計算で Jacobi を本番の反復のまま残すと、予算はセル数の二乗に沿って増えます。

CG の 60 回は $N$ と正確に一致します。「最大 $N$ 回」という上限を残さず使った回数です。1.87e-13 は 2.11e-2 より約 $1.1\times 10^{11}$ 倍小さい値です。倍精度で 1 の近くに足せる最小の数は $10^{-16}$ のあたりです。60 回の直交化で二、三桁が擦り減りました。回数そのものは $N$ で止まっています。

100 万セルに 8 バイトを掛けると、ベクトル一本は 8 MB です。CG のベクトル 4 本は 32 MB、BiCGSTAB の 7 本は 56 MB です。疎行列の横に置くと、残差の保存は安い側です。高い側は、あの 1 万回の尾です。

前処理を外した残差曲線

先の曲線は、前処理を切った裸の反復です。前処理行列が入ると、横軸の一区切りは別の意味になります。BiCGSTAB が垂直に跳ねる点は、分母の内積が 0 に張り付く break-down です。反復を足すだけでは戻りません。前処理を替えるか、GMRES へ移ります。再スタート周期 $m$ の慣行値 30 は、メモリとその跳ねを一緒に見る妥協です。

実務ソルバー選択ガイド

圧力ポアソンと、拡散だけの熱伝達は CG に ILU(0) です。対称性が生きているのに別のソルバーを選ぶと、ベクトルが増えるだけです。

非対称な運動量は BiCGSTAB が基本です。残差が揺れると GMRES( $m$ ) に乗り換えます。 $m$ が再スタート周期です。

並列では Gauss-Seidel が逐次なので、GPU では力が抜けます。ここ 10 年は Jacobi 平滑化と代数マルチグリッド（AMG）の組合せが標準に近いです。

SU2 の LU-SGS（Lower-Upper Symmetric Gauss-Seidel）は、Krylov の脇にある定常反復です。圧縮性流れの陰的時間積分に付けると、メモリの足跡は小さくなります。Krylov ほど強くはありません。

ソルバーの名前より前処理の方が、回数を大きく分けます。PCG に AMG を付けると、裸の CG より 100 倍速いことが珍しくありません。上の 1 万回の尾を切れるかどうかの差です。

残差2.11e-2を収束と渡した失敗#

2000 回後も 2.11e-2 のまま時間ステップを渡すと、その圧力誤差が次の非線形反復へ入ります。SPD が生きていれば席は CG です。移流が行列を非対称へ傾けると、CG はその場で資格を失い、BiCGSTAB へ渡します。名前を先に選び、条件数を後から測る順番が失敗です。

残る問いは、反復をさらに 1 万回積むかどうかです。粗い格子で誤差を掃く Vサイクルが、この $1-O(1/N^2)$ の尾を切るのか、細かい格子で Jacobi を再び呼ぶのか。そのどちらに倒れるかが、次の選択の基準です。