[論文レビュー] tanhで衝撃波をサブセル単位に — TENO-THINC再構成を実装する

精度を 5 次、6 次、8 次へ上げても、接触不連続は時間とともになまっていきます。WENO でも TENO でも同じです。多項式で跳びを描く限り、そうなります。Takagi らの 2023 年論文は、その道を外します。滑らかな領域は高次 TENO、不連続セルだけが tanh(THINC) です。

再構成の切替を Polynomial、THINC、Hybrid に動かすと、同じセル平均が S 字、垂直な跳び、混合に分かれます。

Polynomial は跳びの上下に少しはみ出します。THINC は 1 セルの中でほぼ垂直です。Hybrid は跳びのセルだけ tanh を使います。

論文を一枚で

• タイトル: High-Order Low-Dissipation Shock-Resolving TENO-THINC Schemes for Hyperbolic Conservation Laws
• 著者: Shinichi Takagi, Hiro Wakimura, Lin Fu, Feng Xiao
• 掲載: Communications in Computational Physics, 2023
• 一言: 偶数点TENO(6点・8点)にTHINC再構成を結合し、滑らかな領域の低散逸性と不連続のサブセル解像度を同時に得ます。

タダで分かるのは、いつ差し替えるかです。TENOはすでに候補ステンシルごとに不連続かどうかを $\delta_k \in \{0,1\}$ で印付けています。この値をそのまま再利用してセルが不連続を含むかを見分け、そのセルだけでTHINCを起動します。新しいパラメータはTHINCの勾配 $\beta$ ただ一つです。

なぜ高次多項式でも接触面はなまるのか

多項式再構成は、本質的に滑らかな曲線を描きます。セル平均という点を通る曲線を、ほぼ垂直な跳びの前で描こうとすると、二択を迫られます。振動(オーバーシュート)を許すか、勾配を殺してなまらせるか。

TENOとWENOは振動を抑える側を選びます。不連続をまたぐステンシルを捨て、ENO性を強制します。安定性は得られますが、毎ステップわずかに数値散逸が積もります。接触不連続(速度・圧力は同じで密度だけが跳ぶ面)は、圧力機構で自ら鋭くなることがありません。そのため長時間移流の後には幅が数十セルに広がります。次数を上げても、このなまりはほとんど減りません。

上で選んだ三つの再構成が、その選択です。セル平均は同じで、曲線だけが違います。TENO が小さな渦をあまり殺さない話は 低散逸 TENO にあります。接触面が数十セルへ広がることは、その低散逸でも止まりません。

THINC — tanhがサブセルの跳びを描く#

THINC(Tangent of Hyperbola for INterface Capturing)は、多項式の代わりに単調増加のtanhを使います。セル $i$ の内部座標 $X \in [-0.5, 0.5]$ に対する再構成関数は

ここで $f_a = (f_{i+1}+f_{i-1})/2$ は近傍平均、$f_d = (f_{i+1}-f_{i-1})/2$ は跳びの半分、$\beta$ はtanhの勾配、$d_i$ は跳び中心の位置です。

$d_i$ は自由変数ではありません。セル平均保存の条件(再構成の積分が $f_i$ に戻ること)が $d_i$ を固定します。その解は解析的に

$T_1 = \tanh(\beta/2)$、$T_2 = \tanh(\alpha\beta/2)$、$\alpha = (f_i - f_a)/f_d$ です。$\alpha$ はセル平均が近傍の間でどこにあるかを表す正規化位置です。セルが単調でない場合($(f_{i+1}-f_i)(f_i-f_{i-1}) \le 0$)、tanhは無意味なので1次再構成($\hat f_{i+1/2} = f_i$)に後退します。

β:鋭さを握るつまみ

$\beta$ 一つがTHINCの性格を決めます。小さければ緩やか、大きければ急峻です。論文は意味のある対応を整理しています。

• $\beta \approx 1.1$ : van Leer制限子並みのTVDスペクトル
• $\beta \approx 1.3$ : Superbee並み(過圧縮傾向)
• $\beta = 1.6 \sim 2.0$ : 鋭い衝撃捕捉の領域

論文は $\beta = 1.8$ を採用します。上げすぎると滑らかな曲線まで階段状に squaring され、人工的な平坦化が生じます。上の図で $\beta$ を 1.0 から 2.4 へ動かすと跳びが立ち、overshoot の表示も一緒に動きます。

$\beta = 1.8$ なら $\beta/2 = 0.9$ です。セル中央の跳びでは $\alpha = 0$ なので $d_i = 0$ です。右端 $X = 0.5$ で $\tanh(\beta(X - d_i)) = \tanh(0.9) \approx 0.716$ です。隣の跳びの約 72% だけが 1 セルの中で戻ります。$\beta = 1.1$ なら $\tanh(0.55) \approx 0.50$ で、van Leer のように緩やかです。$\beta = 2.4$ なら $\tanh(1.2) \approx 0.83$ まで立ちます。1.8 はその間、垂直と squaring の中間です。

$\alpha = (f_i - f_a)/f_d$ です。セル平均が隣の真ん中なら $\alpha = 0$ で、跳びの中心も 0 です。平均が右の隣へ寄ると $\alpha$ は 1 に近づき、対数の中の $(1 - T_2/T_1)$ が縮んで $d_i$ が正の側へずれます。$|\alpha| \ge 1$、または $(f_{i+1} - f_i)(f_i - f_{i-1}) \le 0$ なら、コードは tanh を捨てて $\hat f = f_i$ へ戻ります。山のセルはそのため多項式が引き受けます。

初期 top-hat は $0.35$ と $0.65$ の間です。幅 $0.30$、$N = 100$ なら 30 セルです。$dx = 0.01$、$\Delta t = 0.45 \times dx = 0.0045$ です。五周は $5 / 0.0045 \approx 1111$ ステップで、CFL は 0.45 です。そのあと minmod の幅は 23 セル、ピークは 0.84 です。0.5 を超えるセルが 7 つ減り、高さは 0.16 落ちます。THINC-BVD はピーク 1.00、幅 30 を保ちます。コードコメントの二行が、この計算です。

TENOが見つけ、THINCが埋める — BVDハイブリッド#

THINCを全セルに使うと、滑らかな波まで壊れます。だから不連続セルだけで起動しなければなりません。論文は、TENO重み付けがすでに計算した $\delta_k$ からセルごとの指標 $\zeta_i$ を作ります。

$\zeta_i = 1$ は右に、$-1$ は左に不連続がある合図です。左右の近傍の $\zeta$ の符号がセル $i$ を挟んで向かい合うとき、そのセルが本当に不連続を含むと判定し、THINCに置き換えます。

実用のBVD(Boundary Variation Diminishing)はもっと単純にも使えます。二つの再構成を両方計算し、セル境界での左右の再構成値の差(境界変動)をより減らす側を採用します。下の実装はこのBVD判定に従います。

Python — top-hatを一周させる#

線形移流 $u_t + u_x = 0$ でtop-hatを周期境界で回します。minmod TVDとTHINC-BVDを、同じ格子・同じ時間で比べます。

import numpy as np
 
def minmod(p, q):
    return np.where(p * q <= 0, 0.0, np.where(np.abs(p) < np.abs(q), p, q))
 
def thinc_right_edge(u, beta=1.8):
    # 式2.19~2.20: セルiの右面(i+1/2)の左再構成値
    um, up = np.roll(u, 1), np.roll(u, -1)
    fa, fd = 0.5 * (up + um), 0.5 * (up - um)
    mono = ((up - u) * (u - um) > 0) & (np.abs(fd) > 1e-12)
    alpha = np.where(mono, (u - fa) / np.where(mono, fd, 1.0), 0.0)
    mono = mono & (np.abs(alpha) < 1.0)
    T1, T2 = np.tanh(beta / 2), np.tanh(alpha * beta / 2)
    d = np.where(mono, (1 / (2 * beta)) * np.log((1 - T2 / T1) / (1 + T2 / T1)), 0.0)
    vR = fa + fd * np.tanh(beta * (0.5 - d))
    return np.where(mono, vR, u)  # 非単調なら1次へ後退
 
def reconstruct(u, beta, use_thinc):
    um, up = np.roll(u, 1), np.roll(u, -1)
    s = minmod(u - um, up - u)
    poly_vR = u + 0.5 * s
    if not use_thinc:
        return poly_vR
    poly_vL = u - 0.5 * s
    t = thinc_right_edge(u, beta)
    # BVD: 右面の跳びをより減らす再構成を採用
    jump_poly = np.abs(poly_vR - np.roll(poly_vL, -1))
    jump_thinc = np.abs(t - np.roll(poly_vL, -1))
    return np.where(jump_thinc < jump_poly, t, poly_vR)
 
def rhs(u, dx, beta, use_thinc):
    vR = reconstruct(u, beta, use_thinc)        # a=1 なので F_{i+1/2} = vR
    return -(vR - np.roll(vR, 1)) / dx
 
def advect_tophat(N=100, revolutions=5.0, beta=1.8, use_thinc=True):
    dx, x = 1.0 / N, (np.arange(N) + 0.5) / N
    dt = 0.45 * dx
    u = ((x > 0.35) & (x < 0.65)).astype(float)   # 初期top-hat
    for _ in range(int(revolutions / dt)):
        u1 = u + dt * rhs(u, dx, beta, use_thinc)            # SSP-RK2
        u = 0.5 * u + 0.5 * (u1 + dt * rhs(u1, dx, beta, use_thinc))
    return x, u
 
x, u_tvd = advect_tophat(use_thinc=False)
_, u_thinc = advect_tophat(use_thinc=True)
print("TVD   peak=%.3f  width(>0.5)=%d" % (u_tvd.max(), (u_tvd > 0.5).sum()))
print("THINC peak=%.3f  width(>0.5)=%d" % (u_thinc.max(), (u_thinc > 0.5).sum()))
# 出力例:
# TVD   peak=0.84  width(>0.5)=23
# THINC peak=1.00  width(>0.5)=30

5周した後、minmod TVDはピークが0.84に沈み、幅が狭まります。THINC-BVDはピーク1.00を保ち、幅も初期値(30セル)をほぼ守ります。より高次のTENOでも、多項式の本性ゆえにこのなまりを完全には止められない、というのが論文の出発点です。

手で触れる再構成と長時間移流

revolutions を動かすと、赤い minmod とシアンの THINC が同じ時間だけ別々に進みます。

revolutions が増えるほど、赤のminmod曲線は両端が崩れて台形に広がります。シアンのTHINC-BVDは灰色の厳密解にほぼ張り付き、長方形を守ります。$\beta$ を1.1まで下げるとTHINCもTVDのように緩やかになり、2.0付近で最も鮮明になります。

再現していて怪しかったこと

約束は明快ですが、実装していて引っかかった点が三つあります。

BVD 判定は検出器の品質に敏感です。滑らかな極値(smooth extremum)を不連続と誤認すると、その峰が不自然に平坦になります。上の単純な BVD は単調条件でこれを防ぎますが、多次元・系のシステムでは特性変数分解が必要です。

$\beta$ は問題とともに動きます。$\beta = 1.8$ は 1D ベンチマークに合わせた値です。強い圧縮性乱流では過圧縮で人工的な階段が生じうるため、適応的な $\beta$ を検討する価値があります。

コストの主張は文脈次第です。「THINC は不連続セルが少ないのでほぼタダ」という記述は、再構成だけを見れば正しいです。しかし分岐の多い BVD 判定は GPU でワープ発散を招くことがあります。OpenFOAM のような非構造 FV コードに移植するには、セルごとの検出ロジックを face flux ループに溶け込ませる必要があり、1D 構造格子ほどきれいにはいきません。

βを1.8に固定するか#

五周のあと THINC-BVD はピーク 1.00、幅 30 セルです。同じ条件の minmod は 0.84、23 セルです。滑らかな極値まで階段に押されるなら、その固定値は過ぎています。つまみ一つを波ごとに疑うのが、THINC と MP5 を分ける選択です。

跳びを tanh で描くか、にじんだあとを掃くかは別の決定です。再構成がすでに接触を守るなら、後処理の先鋭化は要りません。守らないなら、$\beta$ をさらに上げる前に検出器を見ます。