[論文レビュー] WENOが小さな渦をなめてしまう時

同じ10%でもχは1/4付近に残り、五点ステンシルは切れません。

スタンフォードで Lin Fu が乱流の直接数値計算を回しているとき、妙なものを見ました。5次の WENO-JS は衝撃をよく捕まえました。衝撃から遠い小さな渦まで一緒にぼけました。次数は5と書いてありました。手元に残る解像度は3次の中心差分側でした。

$C_T$ のスライダーを $10^{-5}$ から $10^{-3}$ や $10^{-7}$ へ動かすと、カットオフは厳しくなったり緩んだりし、同じこぶが橙と青緑で別々に削られます。

N cells: 120

CFL: 0.40

TENO cutoff C_T: 1e-5

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

一周のあと、矩形パルスの角がどちらも一、二セルかを見ます。 $10^{-7}$ で青緑パルスの脇に残振動が出るかどうかも見ます。

論文情報と背景

著者: Lin Fu (Stanford Turbulence Research Center)
誌名: 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)は、小さなステンシル $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$ は $10^{-6}$ 付近の、ゼロ除算を防ぐ定数です。この重みは連続です。ステンシルがほんの少し揺れても、 $w_k=d_k$ であるべき滑らかな区間で $w_k$ が動き、より低い次数へ重みが載ります。

滑らかな波が数十ステップ進むと振幅が削れ、位相誤差が積もります。衝撃と乱流渦が同居する LES では、この散逸が致命的です。

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

Fu は二箇所を変えました。

ステンシルの構成を変えました。三つの3点ステンシル $S_0, S_1, S_2$ に、大きな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次ステンシル再構成#

実際のコードは短いです。五点から平滑度三つと候補多項式三つを作り、TENO は $C_T$ で切ります。WENO-JS は逆二乗で混ぜる側が比較対象です。この関数自体は TENO の枝です。

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% ほど低くなります。骨格は同じで、重みの計算だけが違います。

重みの式に小さな揺れを入れます。理想の線形重みは10分の1、10分の6、10分の3です。βが三つとも1なら、WENO の w はその値に戻ります。三つ目だけ10パーセント上げ、βを1、1、1.1にします。三つ目の未正規化重みは0.3を1.21で割った約0.248です。和で割ると w は約0.105、0.633、0.262です。理想重みとの L1 差は約0.077です。不連続ではありません。平滑度指標が10パーセント動いただけです。

同じ数字を TENO に入れます。タウは0.1、q は6です。コードは大きなステンシル指標を三つのβの最大値にし、四つを一緒に割ります。するとカイはおよそ0.256、0.256、0.244、0.244です。カットオフが $10^{-5}$ なら四桁上です。判定はすべて1のままで、大きなステンシルが生き残り、5次線形式をそのまま使います。同じ10パーセントで WENO の重みは揺れ、TENO は1ビットもひっくり返りません。

本文の式のイプシロンは $10^{-6}$ 付近です。関数の中では $10^{-40}$ です。βが1付近なら、どちらも結果に出ません。βが0に張り付くステンシルだけで分母が分かれます。カットオフとは別に、イプシロンは書かれた場所ごとに違います。5次式の係数2、−13、47、27、−3を足すと60です。60で割れば定数場はそのままです。大きなステンシルが生きているときの床です。

C_Tを動かしたときに残る振幅#

上の移流は同じ初期条件、滑らかなサイン二乗のこぶと矩形パルスを線形に進めます。橙が WENO-JS、青緑が TENO-5 です。 $C_T=10^{-5}$ で一周すると、矩形パルスの角はどちらも一、二セルだけにじます。左のサインのこぶは WENO-JS だけが目に見えて低くなります。 $C_T$ を $10^{-3}$ に上げるとカットが厳しくなり、矩形パルスの近くで小さなステンシルへよく落ち、TENO も散逸が増えます。 $10^{-7}$ に下げると TENO はほとんど5次線形に固執し、矩形パルスの脇に残振動が見えます。 $C_T$ が安定と低散逸の天秤です。

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

論文が自信を持って言わない場所を残しておきます。

$C_T$ は問題に付きます。Fu は $10^{-5}$ を勧めます。強い衝撃ととても弱い音波が一緒にあると、この値では音波が少し浮き、 $10^{-6}$ に下げると衝撃の近くで振動が回ります。OpenFOAM や SU2 へ移したグループは、ケースごとにこの値を合わせ直すと書いています。

有限体積では $S_3$ の平滑度 $\beta_3$ が高いです。5点多項式の微分積分が六ページに及びます。有限差分にするとずっと短いです。この論文の低散逸 Riemann フラックス切替まで一緒に使うなら、有限体積の形が自然です。

時間積分との相性は空いています。SSP-RK3 では空間精度に押されて大丈夫です。SSP-RK2 や陰解法と混ぜると、 $\delta_k$ の不連続が Newton 収束に触れるという後報があります。

実務で使うなら、格子収束を同じ $C_T$ で回します。空間次数が実際に5へ戻るかを見ます。戻らなければ、 $C_T$ が大きすぎて TENO が常時非線形モードである可能性が高いです。

C_Tを一つに固定すると壊れる検証#

$C_T=10^{-5}$ を固定したまま格子を細かくすると、滑らかなセルのカイがカットの下へ落ちた瞬間に5次が消えます。上で入れた10パーセントの揺れはカットを越えませんでした。衝撃ととても弱い波が一緒にあると話が変わります。

そのとき次数検証が5ではなく1を下回って戻るなら、食べられたのは空間演算子ではなくカットオフです。同じコードが0.79次と2.00次を同時に出した検証は、ケースが何を測っているかを問います。 $C_T$ を固定した TENO は、その二つの数字のどちらを出すことになるでしょうか。

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