Lax–Wendroffが段差で震えた — Flux Limiterで2次TVD移流スキームを作る

φ(r)がSweby領域の中なら、段差に新しい極値は出ません。滑らかな波の次数はリミッターごとに分かれます。

あるテストケースで、段差の前後に妙なこぶが立ち上がりました。線形移流方程式(linear advection、定数速度場による純粋な並進)をLax–Wendroffで解いていて、初期条件は単純な矩形パルスでした。2次スキームなので1次風上より鋭いはずでしたが、結果は段差の前でアンダーシュート、後ろでオーバーシュートです。ログのCFLは0.5で、安定条件は満たしています。

問題はアルゴリズムが不安定なことではありません。アルゴリズムが単調でないことでした。

リミッターのボタンは lax-wendroff から始まります。CFLスライダーを0.5に置くとオレンジの曲線が段差の前後にこぶを作り、superbee に変えるとそのこぶが消えます。

CFL0.50steps: 0

オレンジが数値解、シアンの点線が厳密解です。minmod は両側が鈍く、こぶがありません。superbee は段差が元より鋭く見えることさえあります。CFLを0.9付近まで上げると、こぶのないリミッターもより広がります。安定条件の内側でも、数値拡散と噛み合います。

振動の正体 — Godunovが引いた壁#

1次元線形移流方程式は

\partial_t q + u\,\partial_x q = 0

速度 $u$ が一定ならば、解は初期形状がそのまま右へ平行移動する。数値的にはセル境界のフラックス $f_{i-1/2}$ を定義し、

q_i^{n+1} = q_i^n - \frac{\Delta t}{\Delta x}\bigl(f_{i+1/2} - f_{i-1/2}\bigr)

この形式は flux-conserving form と呼ばれ、質量・運動量・エネルギーなどの保存量を機械精度まで保つ。Lax–Wendroffはここに2次テイラー補正を加えた線形スキームだ。

ところが1959年にGodunovが証明した定理は厳しい。線形定係数移流に対して、単調性(monotonicity)を保ちつつ2次以上の精度を持つ線形スキームは存在しない。 どちらかを諦めるしかない。Lax–Wendroffは2次精度を取って単調性を捨てた。だから段差で振動が出る。

TVD条件:Hartenの不等式#

Harten(1983)が示した抜け道は「全変動が減少する」スキームだ。

\mathrm{TV}(q) = \sum_i |q_{i+1} - q_i|, \qquad \mathrm{TV}(q^{n+1}) \le \mathrm{TV}(q^n)

左辺は数値解のうねりの総和。毎ステップこの値が増えなければ新しい極値(オーバーシュート・アンダーシュート)は生じない。この条件を満たすスキームをTVDと呼ぶ。

Godunovの壁を迂回するにはスキームを非線形にする必要がある。解の滑らかさに応じてアルゴリズムが形を変える、というアイデアだ。滑らかな領域では2次、急な勾配の近くでは1次風上。この切り替えを担うのが flux limiter である。

Flux Limiterを一行で理解する#

Lax–Wendroffのフラックスは1次風上に補正項を足した形だ。

f_{i-1/2} = u\,q_{i-1} + \tfrac{1}{2}\,|u|\!\left(1 - \left|\tfrac{u\Delta t}{\Delta x}\right|\right) \varphi(r)\,(q_i - q_{i-1})

ここで $u$ は移流速度、 $\Delta x$ はセル幅、 $\varphi(r)$ がリミッター、 $r$ は前後の勾配比。

r_{i-1/2} = \frac{q_{i-1} - q_{i-2}}{q_i - q_{i-1}}

解が滑らかなら $r \approx 1$ 、極値の近くでは $r < 0$ あるいは $r \gg 1$ 。リミッターはこの $r$ を見て補正の強さを調整する。

特別な選択として、

$\varphi(r) = 0$ :donor-cell(1次風上)
$\varphi(r) = 1$ :Lax–Wendroff
$\varphi(r) = r$ :Beam–Warming

Sweby図で $0 \le \varphi(r) \le \min(2r,\,2)$ を満たせばTVDが保証される。この台形領域に入る代表的なリミッター4つを比較する。 $\varphi(r) = 0$ で切った場所は1次風上です。その風上が振動の代わりに何を鈍らせるかは、FTCSが震えるときの引き換えにあります。

リミッター	定義	特徴
minmod	$\max(0,\min(1,r))$	最も慎重。段差はよく保つが滑らかな領域はやや鈍い
superbee	$\max(0,\min(2r,1),\min(r,2))$	最も鋭い。滑らかな波を階段状に削ることも
van Leer	$(r+	r
MC	$\max(0,\min(\tfrac{1+r}{2},2,2r))$	既定値として最も無難

Pythonで直接比較#

rollベースのベクトル化で40行以内に収まる。200セル、 $\mathrm{CFL} = 0.5$ 、500ステップ(周期境界、波が領域を約2.5周する距離)。

import numpy as np
 
def limiter_phi(name, r):
    r = np.where(np.isfinite(r), r, 0.0)
    if name == 'donor':    return np.zeros_like(r)
    if name == 'lw':       return np.ones_like(r)
    if name == 'minmod':   return np.maximum(0, np.minimum(1, r))
    if name == 'superbee':
        return np.maximum.reduce([np.zeros_like(r),
                                  np.minimum(2*r, 1),
                                  np.minimum(r, 2)])
    if name == 'vanleer':  return (r + np.abs(r)) / (1 + np.abs(r) + 1e-30)
    if name == 'mc':
        return np.maximum(0, np.minimum.reduce([0.5*(1+r),
                                                2*np.ones_like(r),
                                                2*r]))
 
def advect_limited_fvm(q, cfl, name):
    """1D FV 移流、周期境界、u = +1。"""
    dq  = q - np.roll(q, 1)                  # 面 i-1/2: q_i - q_{i-1}
    r   = np.roll(dq, 1) / (dq + 1e-30)      # (q_{i-1}-q_{i-2}) / (q_i-q_{i-1})
    phi = limiter_phi(name, r)
    flx = np.roll(q, 1) + 0.5 * phi * (1 - cfl) * dq
    return q - cfl * (np.roll(flx, -1) - flx)
 
N, cfl, steps = 200, 0.5, 500
x  = np.arange(N) / N
q0 = ((x > 0.2) & (x < 0.4)).astype(float)
for name in ['lw', 'minmod', 'superbee', 'vanleer', 'mc']:
    q = q0.copy()
    for _ in range(steps):
        q = advect_limited_fvm(q, cfl, name)
    l1 = np.abs(q - q0).sum() / N
    print(f'{name:<9s}  L1={l1:.4f}  overshoot={q.max()-1:.3f}  undershoot={q.min():.3f}')

代表的な実行結果:

リミッター	$L_1$ 誤差	オーバーシュート	アンダーシュート
lax-wendroff	0.056	+0.181	−0.133
minmod	0.081	0.000	0.000
superbee	0.042	+0.002	0.000
van Leer	0.061	0.000	0.000
MC	0.051	0.000	0.000

オーバーシュートが20%近く跳ねるのはlax-wendroffだけで、他は全部ゼロ。誤差の順位も面白い。superbeeが $L_1$ で最小だ。段差をにじませないということで、代わりに滑らかなsine波では波頭を階段状に歪める。

滑らかな波のテスト( $q_0 = \sin(2\pi x)$ 、500ステップ)を別に回すと、収束次数はminmodで約1.5、van Leer/MCで約2.0、superbeeは約1.3まで落ちる。「2次精度」のラベルは極値近傍でリミッターが働くたびに部分的に1次へ落ちるため、こうなる。

r=1 では四つのリミッターが同じ φ を返す#

表の定義に $r=1$ を入れると、minmod、superbee、van Leer、MC はどれも 1 になります。 $\varphi=1$ は Lax–Wendroff です。滑らかな真ん中では四つのリミッターが同じ補正を使います。違いは $r$ が 1 を離れたときに出ます。

$r=2$ を入れてみます。minmod は 1 で止まります。superbee は 2 で Sweby の上端に付きます。van Leer は $4/3$ です。MC は 1.5 です。superbee が段差の $L_1$ を 0.042 と最も小さく保つ理由は、この上端です。同じ選択が sine テストの次数を約 1.3 まで下げます。van Leer と MC は約 2.0 に残り、minmod は約 1.5 です。 $r$ が 0 以下の極値では四つとも 0 に落ち、donor-cell になります。

コードの補正係数もすでに決まっています。速度が 1 で $\mathrm{CFL}=0.5$ なら、ジャンプの前に付く係数は 0.25 です。 $\varphi=1$ でも隣のジャンプの 4 分の 1 しか足しません。500 ステップが進む距離は $0.5\times 500=250$ セルです。格子 200 なら周期境界を 1.25 周します。オーバーシュート +0.181 は波高 1 の 18 パーセントです。アンダーシュート −0.133 は床が負まで沈んだ深さです。minmod の $L_1$ 0.081 は superbee の 0.042 の約 1.9 倍です。こぶを消した代償として、段差がその分にじみます。

実務チェックリスト

衝撃捕獲が目的なら minmod または van Leer から。振動ゼロが最優先。
LES・大規模乱流で superbee は危険。滑らかな乱流構造を階段状に削り、エネルギースペクトルを歪める。van Leer または MC が妥協点。
境界セルでステンシルが足りず暫定的に1次風上を使うと、全体の精度が1次まで落ちる。境界でもリミッターのステンシルを確保するか、ゴーストセルを使うこと。
OpenFOAMの limitedLinear 1 はvan Leer系の1パラメータリミッター。係数0で1次風上、1.0でTVD領域全体を使う。
収束次数の検証は必ず滑らかな解で。不連続を含むテストで $L_1$ 傾きを測ると、リミッターではなく不連続そのものが次数を落としてしまう。

段差に superbee を使うと、sine は 1.3 次まで落ちる#

表で superbee の $L_1$ は 0.042 と最小です。オーバーシュートは +0.002 で、段差の上ではほとんど見えません。同じリミッターを sine に載せると次数は約 1.3 まで落ちます。van Leer と MC は約 2.0 に残ります。不連続を守る選択が、滑らかな波の次数を削ります。

スカラーが密度なら、lax-wendroff のアンダーシュート −0.133 は負の密度の入口です。その穴を塞ぐ限定子は正の値が漏れる場所に別にあります。受動スカラーの段差なら minmod か van Leer で足ります。密度なら、その次の扉を開けます。