境界層格子でJacobiが止まる場所 — ラインレット線陰的前処理

点スムーザの天井に残るのは層数Nyだけで、ラインレットの上界1/3にアスペクト比はありません。

同じ方程式です。セル数も反復法も収束判定も変えていません。格子を壁側に一度潰しただけで、反復回数が18回から5,000回に跳ね上がりました。

物理は何も変わっていません。変わったのはセルのアスペクト比だけです。それなのに線形ソルバーはまるで別の問題を渡されたかのように振る舞い始めます。しかも避けようがありません。乱流解析で $y^+ \approx 1$ を狙えば、壁面第一セルは潰すしかないからです。

その300倍が行列のどこにいるのか、係数を追います。SU2 などが使うラインレット(linelet・結合の強い方向にセルを連ねた鎖)が300倍を返すかも、反復回数を数えて見ます。短い行は表、長い行は係数です。

セルを潰すと行列は別物になる

矩形セル上で拡散項を有限体積離散化すると、面の係数はこうなります。

a_E = a_W = \frac{\Gamma \Delta y}{\Delta x}, \qquad a_N = a_S = \frac{\Gamma \Delta x}{\Delta y}

$\Gamma$ は拡散係数、 $\Delta x$ ・ $\Delta y$ はセルの寸法です。比を取ってみます。

\frac{a_N}{a_E} = \left(\frac{\Delta x}{\Delta y}\right)^{2} = \mathrm{AR}^{2}

$\mathrm{AR}$ はセルのアスペクト比(流れ方向/壁垂直方向)です。係数比はその二乗で効きます。壁面第一セルが $\mathrm{AR} = 1000$ なら、南北方向の結合は東西方向の百万倍です。格子を10倍潰すと行列は100倍歪みます。

AR スライダーを 1 から 100 へ動かすと縦縞が先に消え、linelet ボタンを入れると横に残った染みの寿命が縮みます。

log₁₀ AR100ω1.00

What to watch: at AR = 1 all three iterations finish in tens of sweeps. Drag AR to 100 and the two point smoothers flatten out — the wall-normal stripes vanish almost at once, but the wall-parallel streaks just sit there and the measured factor climbs to 0.997. Switch to linelet at the same AR and it lands in single digits. The two bounds under the map say why: AR is in the point-smoother one and absent from the linelet one.

図の下の measured が 0.997 付近に張り付いたら停滞です。同じ AR で linelet を押すと、その染みは 1、2 sweep で消えます。

ここに時間微分を陰的に入れると、対角成分はこうなります。

a_P = \frac{V}{\Delta t} + 2a_E + 2a_N

$V = \Delta x \Delta y$ はセル体積、 $\Delta t$ は(擬似)時間刻みです。この $V/\Delta t$ が対角を支える唯一の仕掛けですが、 $a_N$ の増加速度には追いつきません。

点スムーザの収束係数にはARが入っている#

Jacobi 反復の誤差減衰は、対角に対する非対角の総和で上から押さえられます。

\rho_{\text{point}} \le \frac{2a_E + 2a_N}{V/\Delta t + 2a_E + 2a_N}

$V/\Delta t = 4a_E$ と置けば、 $\mathrm{AR} = 1$ で 0.5、 $\mathrm{AR} = 100$ で 0.9998 です。ただしこの上界は $\mathrm{AR}$ が大きくなると 1 に貼り付いてしまい、何も教えてくれなくなります。実際に何が残るのかはモード別に見る必要があります。Fourier モード $(\theta_x, \theta_y)$ に対する減衰因子は

\rho(\theta_x, \theta_y) = 1 - \frac{V/\Delta t + 2a_E(1 - \cos\theta_x) + 2a_N(1 - \cos\theta_y)}{a_P}

最も遅いモードは $\theta_y$ が最小のもの、つまり壁垂直方向に最も緩やかなモードです。 $a_N$ が対角を支配して $a_P \approx 2a_N$ になると

\rho_{\max} \approx 1 - \frac{\pi^{2}}{2 N_y^{2}}

$N_y$ は壁垂直方向の層数です。ここに $\mathrm{AR}$ が現れないのは、良い知らせではなく悪い知らせです。これ以上潰しても係数はこの値を超えませんが、その天井自体がすでに使い物になりません。 $N_y = 48$ なら 0.99786、誤差を $10^{-6}$ まで落とすのに約6,000回です。

何が起きているかは単純です。点スムーザは壁垂直方向の高周波誤差を1、2回で消し去ります。問題はその先です。壁に沿う方向の誤差を消すには情報が $a_E$ を伝って横に進む必要がありますが、対角が $a_N$ に押さえられているため、一歩の大きさが $\mathrm{AR}^2$ 分だけ縮みます。

対角だけを見る反復が手法ごとにどこで止まるかは、Jacobi から BiCGSTAB まで収束速度を比べた記録にあります。あちらはソルバーを替えます。こちらはセルのアスペクト比が同じソルバーを裏切ります。

三つの処方を同じ表に並べる

前処理	反復あたりのコスト	異方性への対応	並列性	引っかかる箇所
Jacobi(対角)	セルあたり除算1回	なし	完全	$\mathrm{AR} \gtrsim 10$ で崩れる
ILU(0)	分解1回 + 前進・後退代入	部分的	番号付けに依存	領域分割ごとに性能が変わる
ラインレット	ライン内セルあたり Thomas 5 flop	強い方向を正確に除去	ライン単位で独立	ラインの外は依然 Jacobi

ラインレットの取引条件は明快です。結合の強い方向を一つだけ正確に反転します。その代わり記憶は3バンド一組、コストは実質的に Jacobi と同じ桁です。ILU(0) のように行列全体に手を入れないので、MPI 分割にも鈍感です。

ラインレットはどう作られるか

各セルで面係数比 $\sigma = a_f / a_{f,\perp}$ を測ります。壁垂直面なら $\sigma = \mathrm{AR}^2$ です。
$\sigma \ge \sigma_{\min}$ の面を強い面として印を付けます。実務での既定値は $\sigma_{\min} = 10$ 付近です。
壁に接するセルから出発し、強い面をたどって上へ鎖を伸ばします。
1セルが持てる強い面は最大2本までです。分岐が生じたら係数の大きい側だけを残します。分岐を許すと結果が三重対角になりません。
既に他のラインレットに属するセルに到達したら鎖を切ります。
長さ1の鎖は捨て、そのセルは Jacobi に戻します。

鎖の中だけで番号を振り直すと、その部分行列は三重対角になります。三重対角なら Thomas 法(前進消去のあと後退代入)で $O(n)$ で正確に反転できます。ラインレット前処理の中身はこの一文がすべてです。

鎖がどこで切れるかは反復の表の次に置きます。手は第一層厚さと成長率です。

Python で反復回数を数える#

$64 \times 48$ の格子に誤差だけを残し(右辺は $f = 0$ なので厳密解は $u = 0$ )、点 Jacobi とライン Gauss-Seidel をそれぞれ回して、誤差が $10^{-6}$ まで落ちるまでの sweep 数を数えます。

import numpy as np
 
def cell_coefficients(dx, dy, gamma=1.0, diag_factor=4.0):
    """矩形セル1個のFVM拡散係数。diag_factor は V/dt を a_E の倍数で指定。"""
    a_e = gamma * dy / dx        # 東/西の面
    a_n = gamma * dx / dy        # 北/南の面
    d0 = diag_factor * a_e       # V/dt 項
    return a_e, a_n, d0
 
def residual_field(u, a_e, a_n, d0):
    r = -d0 * u
    r[1:-1, 1:-1] -= a_e * (2*u[1:-1, 1:-1] - u[:-2, 1:-1] - u[2:, 1:-1])
    r[1:-1, 1:-1] -= a_n * (2*u[1:-1, 1:-1] - u[1:-1, :-2] - u[1:-1, 2:])
    r[0, :] = r[-1, :] = 0.0     # Dirichlet 境界
    r[:, 0] = r[:, -1] = 0.0
    return r
 
def jacobi_sweep(u, a_e, a_n, d0, omega=1.0):
    r = residual_field(u, a_e, a_n, d0)
    u += omega * r / (d0 + 2*a_e + 2*a_n)
 
def thomas(a, b, c, d):
    """三重対角の厳密解法: 前進消去のあと後退代入。"""
    n = len(d)
    cp, dp = np.empty(n), np.empty(n)
    cp[0], dp[0] = c[0]/b[0], d[0]/b[0]
    for k in range(1, n):
        m = b[k] - a[k]*cp[k-1]
        cp[k] = c[k]/m
        dp[k] = (d[k] - a[k]*dp[k-1])/m
    x = np.empty(n)
    x[-1] = dp[-1]
    for k in range(n-2, -1, -1):
        x[k] = dp[k] - cp[k]*x[k+1]
    return x
 
def linelet_sweep(u, a_e, a_n, d0):
    """壁垂直方向の鎖を丸ごと正確に反転する(ライン Gauss-Seidel)。"""
    nx, ny = u.shape
    m = ny - 2
    dg = d0 + 2*a_e + 2*a_n
    a = np.full(m, -a_n)
    b = np.full(m, dg)
    c = np.full(m, -a_n)
    a[0] = c[-1] = 0.0
    for i in range(1, nx-1):
        rhs = a_e * (u[i-1, 1:-1] + u[i+1, 1:-1])   # ライン外の結合は右辺へ
        u[i, 1:-1] = thomas(a, b, c, rhs)
 
def count_sweeps(ar, kind, nx=64, ny=48, tol=1e-6, max_sweeps=20000):
    dx, dy = 1.0/nx, 1.0/nx/ar
    a_e, a_n, d0 = cell_coefficients(dx, dy)
    u = np.random.default_rng(7).standard_normal((nx, ny))
    u[0, :] = u[-1, :] = 0.0
    u[:, 0] = u[:, -1] = 0.0
    e0 = prev = np.linalg.norm(u)
    for k in range(1, max_sweeps + 1):
        jacobi_sweep(u, a_e, a_n, d0) if kind == 'jacobi' else linelet_sweep(u, a_e, a_n, d0)
        e = np.linalg.norm(u)
        rho, prev = e/prev, e
        if e/e0 < tol:
            return k, rho
    return max_sweeps, rho
 
print(f"{'AR':>6} {'Jacobi':>8} {'rho_J':>8} {'linelet':>8} {'rho_L':>8}")
for ar in (1, 10, 100, 1000):
    nj, rj = count_sweeps(ar, 'jacobi')
    nl, rl = count_sweeps(ar, 'linelet')
    print(f"{ar:>6} {nj:>8} {rj:>8.4f} {nl:>8} {rl:>8.4f}")

出力はこうなりました。

    AR   Jacobi    rho_J  linelet    rho_L
     1       18   0.4870        8   0.1839
    10      514   0.9780        7   0.1699
   100     4879   0.9975        4   0.0194
  1000     5471   0.9978        2   0.0002

Δy_wall631 µmgrowth1.25log₁₀ σ_min32

What to watch: shrink Δy_wall and the cyan chain grows downward into the layer while the bars on the right cross the red threshold. Raise the growth ratio and the chain gets shorter — the mesh goes isotropic sooner, so there is less for the linelet to own. Push σ_min past a few thousand and the chains disappear entirely: the preconditioner quietly degrades back to plain Jacobi, which is the failure mode nobody notices in the log file.

$\Delta y_{\text{wall}}$ を小さくすると鎖は境界層の奥へ深くなります。成長率を上げると格子が早く等方に戻り、鎖は短くなります。鎖は短いです。全層の半分も覆いません。ラインレットが安いのはその長さです。

$\mathrm{AR}$ を1000倍潰すと Jacobi は300倍遅くなり、前節の予測どおり6,000回付近で天井に当たります。ラインレットはむしろ速くなります。隣り合う鎖の間の結合が $a_E/a_N = \mathrm{AR}^{-2}$ で小さくなり、鎖同士が独立していくからです。最悪モード基準の上界 $2a_E/(V/\Delta t + 2a_E) = 1/3$ はそのまま成り立ち、実測の減衰はそれよりずっと良いだけです。

Jacobi の列だけ割ってみます。18、514、4879、5471。アスペクト比は 1、10、100、1000 です。

最初の 10 倍では反復が 28.6 倍になります。次の 10 倍は 9.5 倍。最後の 10 倍は 1.12 倍で、5471 を 4879 で割った値です。係数比は $\mathrm{AR}^2$ で伸び続けるのに、反復は追いません。 $N_y = 48$ のモード式 $1 - \frac{\pi^{2}}{2 N_y^{2}}$ は 0.99786 です。 $0.99786^{n} = 10^{-6}$ を解くと n は 6400 付近です。表の 5471 はそのすぐ下です。格子をさらに潰しても、この列はここに張り付きます。

ラインレットの列は逆です。8、7、4、2。減衰は 0.1839、0.1699、0.0194、0.0002。上界 $2a_E/(V/\Delta t + 2a_E) = 1/3$ に $\mathrm{AR}$ はなく、実測はすべてその下です。AR 1000 の 0.0002 は 1/3 の約 0.06% です。鎖の間の結合 $a_E/a_N = \mathrm{AR}^{-2}$ が消えるので、Thomas 一回で壁垂直の誤差はほとんど残りません。AR 1 では 8 回対 18 回。得は小さいです。AR 10 では 514 対 7、約 73 倍。処方が分かれる場所です。

点スムーザの上界も同じ対角です。 $V/\Delta t = 4a_E$ なら分母は $6a_E + 2a_N$ です。AR 1 の上界 0.5 に対し、表の rho_J は 0.4870。ほぼ同じです。AR 100 の上界は 0.9998、実測 rho_J は 0.9975 です。上界が 1 に張り付いたあとは、モード式の 0.99786 のほうが正直です。前節の約 6,000 回と表の 5471 は同じ天井を指しています。

実装で踏みやすい落とし穴

ラインの向きを逆に敷く。 鎖は係数の大きい方向、つまり壁垂直方向に敷きます。壁に沿って敷けば利得はゼロで Thomas の費用だけが増えます。格子は見た目には流れ方向に「長い」ので、取り違えやすい箇所です。

$\sigma_{\min}$ を既定値のまま放置する。 格子を変えたのに閾値を見直さないと、鎖が丸ごと消えることがあります。すると前処理は静かに Jacobi へ退化し、ログには何の警告も出ません。ラインレット本数と平均長さを毎回出力しておけば、これは一行で捕まえられます。

MPI 分割が鎖を切る。 分割境界で切られた鎖は断片になり、コア数を増やすほど収束が悪化します。パーティショナにライン方向の重みを与えるか、ラインを切らない制約を課してください。

前処理で非線形問題を覆い隠そうとする。 ラインレットを入れても発散するなら、原因は線形ソルバーではなく外側ループであることがほとんどです。まず CFL を下げ、緩和係数を 0.7 付近まで落として更新量そのものを縮めてから判断してください。前処理は与えられた行列を速く解くだけで、間違った行列を直してはくれません。

壁垂直でなければ Thomas は無駄金です#

鎖は係数の大きい側、壁垂直にだけ敷きます。流れ方向に敷くと、目が長い辺を追っただけです。Thomas の費用は出て、減衰は Jacobi と同じです。

ログには鎖の本数と平均長さを残します。0 が出たら $\sigma_{\min}$ が格子とずれています。前処理はもう対角スムーザです。その状態でコアだけ増やすと、切れた鎖はさらに細かく切れます。格子を詰めながら測ると効率が 16.3% まで落ちる weak scalingが、その請求書です。