[論文レビュー] WENOが小さな渦をなめてしまう時 — Fu(2019) TENO低散逸有限体積法

5次WENO-JSが滑らかな領域でも過剰に散逸する理由と、TENOが鋭いカットオフでそれを解決する方法

Stanford Center for Turbulence ResearchでLin Fuが乱流の直接数値シミュレーションを回していたら、奇妙なことに気づきました。5次WENO-JSは衝撃波をきちんと捉えるものの、衝撃波から離れた小さな渦も時間とともにぼやけていきます。スキームは5次精度を謳っているのに、実際の解像度は3次中心差分のような感触でした。この記事では、その原因を突き止めたFu(2019)論文をレビューし、同じ初期条件でWENO-JSとTENO-5を並列に動かして比較します。結論を先に言えば「非線形重みの連続性が隠れたコストだった」ということです。

論文情報と背景#

著者: Lin Fu(Stanford Center for Turbulence Research)
学術誌: Computer Physics Communications 235 (2019) 25–39
DOI: 10.1016/j.cpc.2018.10.009
キーワード: TENO, High-order schemes, Shock-capturing, Low-dissipation, Finite-volume method

WENO-JSの隠れた損失#

古典的なWENO-JS(Jiang–Shu, 1996)は、3つの小ステンシル $S_0, S_1, S_2$ それぞれで滑らかさ指標 $\beta_k$ を計算し、非線形重みで多項式を混ぜ合わせます。

w_k = \frac{\alpha_k}{\sum_j \alpha_j}, \qquad \alpha_k = \frac{d_k}{(\beta_k + \varepsilon)^2}

ここで $d_k$ は5次精度を回復するための最適線形重み、 $\varepsilon \sim 10^{-6}$ はゼロ除算を防ぐ定数です。問題は、この重みが連続的だという点にあります。ステンシルがほんの少しぶれるだけで、理想的には $w_k = d_k$ を保つべき滑らかな領域でも $w_k$ がわずかに動き、本来より低次のステンシルに重みが乗ります。

結果として、滑らかな波が数十ステップ進むうちに振幅が削られ、位相誤差が蓄積します。LES(大型渦シミュレーション)のように衝撃波と乱流渦が共存するシミュレーションでは、この散逸は致命的です。

TENOの2つのアイデア: 強いスケール分離とシャープカットオフ#

Fuは2つの変更を加えました。

第一に、ステンシル構成を変えました。3つの小3点ステンシル $S_0, S_1, S_2$ に加え、**1つの大きな5点ステンシル $S_3$ **を別に置きます。 $S_3$ がクリーンであればそれで終了 — 他のステンシルは使いません。

第二に、重みを二値化しました。滑らかさの判定はカットオフパラメータ $C_T$ で行います。

\gamma_k = \left(1 + \frac{\tau_K}{\beta_k + \varepsilon}\right)^q, \qquad \chi_k = \frac{\gamma_k}{\sum_j \gamma_j}

\delta_k = \begin{cases} 0, & \chi_k < C_T \\ 1, & \text{otherwise} \end{cases}

$\tau_K = |\beta_0 - \beta_2|$ はグローバル基準指標、 $q=6$ で比率を強調します。要点は $\delta_k$ が0か1しか取らないことです — ステンシルは完全に捨てられるか、最適重みで生かされるかのどちらかになります。連続的な遷移がないため、滑らかな領域では $S_3$ の5次線形スキームをそのまま使い、不連続付近でのみENOライクなステンシル選択に切り替わります。

連続重みから0/1判定へ#

$S_3$ が滑らかと判定されると( $\delta_3 = 1$ )、再構成は単純な5次線形スキームになります。

u_{i+1/2}^L = \frac{1}{60}\bigl(2\bar u_{i-2} - 13\bar u_{i-1} + 47\bar u_i + 27\bar u_{i+1} - 3\bar u_{i+2}\bigr)

$\bar u_j$ はセル平均値です。このスキームは完全に線形なので、 $C_T$ が0に近いと中心差分のように非散逸的に振る舞います。一方 $\delta_3 = 0$ の場合、生き残った小ステンシルだけを最適線形重み $d_k$ で混ぜます。

u_{i+1/2}^L = \sum_{k=0}^{2} w_k\, u_{k,i+1/2}^L, \qquad w_k = \frac{d_k\, \delta_k}{\sum_j d_j\, \delta_j}

これがENO性を保証します。不連続を横切るステンシルがあればその $\delta_k$ が0になり、他のステンシルに滑らかに引き継がれます。

NumPyで5次ステンシル再構成#

実際のコードは短いです。5つのセルから3つの滑らかさ指標と3つの候補多項式を計算し、TENOは $C_T$ で切り、WENO-JSは逆二乗で混ぜます。

import numpy as np
 
def teno_beta_indicators(um2, um1, u0, up1, up2):
    # 3点ステンシル3つのJiang-Shu滑らかさ指標
    b0 = (13/12)*(um2 - 2*um1 + u0)**2 + 0.25*(um2 - 4*um1 + 3*u0)**2
    b1 = (13/12)*(um1 - 2*u0 + up1)**2 + 0.25*(um1 - up1)**2
    b2 = (13/12)*(u0 - 2*up1 + up2)**2 + 0.25*(3*u0 - 4*up1 + up2)**2
    return b0, b1, b2
 
def teno_weighted_reconstruction(u, i, CT=1e-5, q=6):
    N = len(u)
    um2, um1, u0 = u[(i-2) % N], u[(i-1) % N], u[i]
    up1, up2 = u[(i+1) % N], u[(i+2) % N]
    b0, b1, b2 = teno_beta_indicators(um2, um1, u0, up1, up2)
    b3 = max(b0, b1, b2)
    tau = abs(b0 - b2)
    eps = 1e-40
    g = np.array([(1 + tau/(b + eps))**q for b in (b0, b1, b2, b3)])
    chi = g / g.sum()
    delta = (chi >= CT).astype(float)
    if delta[3] == 1:  # S3が滑らか -> 5次線形をそのまま使用
        return (2*um2 - 13*um1 + 47*u0 + 27*up1 - 3*up2) / 60
    d = np.array([1/10, 6/10, 3/10])
    alpha = d * delta[:3]
    if alpha.sum() < 1e-14:
        return u0  # 最後の手段としてのドナーセル
    w = alpha / alpha.sum()
    p0 = (2*um2 - 7*um1 + 11*u0) / 6
    p1 = (-um1 + 5*u0 + 2*up1) / 6
    p2 = (2*u0 + 5*up1 - up2) / 6
    return w[0]*p0 + w[1]*p1 + w[2]*p2

オイラー時間発展と周期境界を組み合わせて回すと、滑らかなガウス突起はTENOでほぼ原形のまま戻ってくるのに対し、WENO-JSでは約30%低くなります。同じコードの骨格で、重み計算だけが違うのに、結果は体感できるほど違います。

下のシミュレーションでC_Tを動かしてみよう#

以下のインタラクティブシミュレーションは、同じ初期条件(滑らかなsin²バンプ+矩形パルス)を線形移流で移動させ、WENO-JS(オレンジ)とTENO-5(シアン)を比較します。

N cells: 120

CFL: 0.40

TENO cutoff C_T: 1e-5

WENO-JSTENO-5L2 err — WENO: 0.0000 | TENO: 0.0000

$C_T = 10^{-5}$ (デフォルト)で一周回してみてください。矩形パルスのエッジは両スキームとも1〜2セル程度しかぼやけませんが、左側のサインバンプはWENO-JSで明らかに低くなっています。 $C_T$ を $10^{-3}$ に上げるとTENOのカットオフが厳しくなり、矩形パルス付近で小ステンシルに落ちる頻度が増えてTENOの散逸も増えます。逆に $10^{-7}$ に下げるとTENOはほぼ常に5次線形スキームを使いますが、この値では矩形パルス周辺にわずかな振動が見えます。 $C_T$ が安定性と低散逸性の天秤であることが直観的に伝わります。

批判的考察 — チューニングの罠#

論文が自信を持って主張していない部分を押さえておきましょう。

第一に、 $C_T$ は問題依存です。Fuは $C_T = 10^{-5}$ を推奨していますが、強い衝撃波と非常に弱い音波が共存する問題では、この値で音波がわずかに浮き、 $10^{-6}$ に下げると衝撃波近傍で振動が出ます。OpenFOAMやSU2にTENOを移植するグループは、ケースごとの再調整が必要だと報告しています。

第二に、有限体積枠組みでの $S_3$ の滑らかさ指標 $\beta_3$ 計算はコストが高いです。5点多項式の導関数の積分が必要で、式は印刷して6ページに及びます。有限差分で書き直せばはるかに簡単ですが、論文の低散逸Riemannフラックス切り替え戦略も一緒に使いたいなら、有限体積構造が自然です。

第三に、時間積分との相互作用が完全には検証されていません。SSP-RK3なら空間精度に押されて問題ありませんが、SSP-RK2や陰解法と組み合わせると $\delta_k$ の不連続性がNewton収束に影響するという事後報告があります。

実務で使うなら: 格子収束研究を同じ $C_T$ で回し、空間次数の再現が実際に5次になるか確認してください。ならない場合、 $C_T$ が大きすぎてTENOが常に非線形モードで動作している可能性が高いです。

覚えておくこと#

WENO-JSの散逸はバグではなく、連続重みが構造的に支払う対価です。小さな $\beta$ の変化でも $w_k$ が微妙に動き、最適線形結合からずれます。
TENOはカットオフ $C_T$ でステンシルを二値化します。滑らかな領域では5次線形スキームをそのまま使い、不連続付近でのみENO選択に切り替わります。
$C_T$ 一つで安定性と低散逸性が天秤にかかるため、格子収束とともに感度分析を必ず行うべきです。