[论文评述] 用 tanh 把激波画到亚网格 — 实现 TENO-THINC 重构

把精度阶数提到五阶、六阶、八阶，接触间断仍会随时间被抹平。WENO 也好，TENO 也罢，结果一样。只要用多项式去画跳变，这个命运就注定了。Takagi 等人 2023 年的论文绕开了它。光滑区域交给高阶 TENO；只在间断网格里换成 tanh 曲线（THINC）。本文在单个网格的层面动手实现这一核心思想，并让一个 top-hat 绕网格转一圈，直接比较多项式与 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（六点、八点）与 THINC 重构结合，同时获得光滑区的低耗散性与间断处的亚网格分辨率。

诀窍在于免费判定"何时切换"。TENO 已经为每个候选模板用 $\delta_k \in \{0,1\}$ 标记是否光滑。直接复用这些标志即可分辨网格是否藏着间断，THINC 只在那里启动。唯一的新参数是 THINC 的斜率 $\beta$ 。

为何高阶多项式也抹平接触面#

多项式重构本质上画的是光滑曲线。当一条曲线必须穿过近乎垂直跳变两侧的网格平均点时，它只有两个选择：允许振荡（过冲），或压低斜率去抹平。

TENO 与 WENO 选择抑制振荡。它们丢弃跨越间断的模板，强制 ENO 性质。这换来稳定性，但每一步都积累一点数值耗散。接触间断（速度、压力相同而密度跳变的面）不会通过压力机制自我变陡。于是经过长时间对流后，它的宽度蔓延到数十个网格。提高阶数几乎无法减少这种抹平。

下面亲自看看。保持相同的网格平均，只改变重构方式。

β = 1.80overshoot: 0.000

Polynomial 把跳变化解成光滑的 S 曲线，并在上下略微越界。切到 THINC，同样的数据却在单个网格内几乎竖直。Hybrid (BVD) 只在跳变网格用 tanh，其余保持多项式。

THINC — tanh 画出亚网格的跳变#

THINC（Tangent of Hyperbola for INterface Capturing）用单调递增的 tanh 取代多项式。网格 $i$ 内部坐标 $X \in [-0.5, 0.5]$ 上的重构函数为

h_i(X) = f_a + f_d \tanh\!\big(\beta(X - d_i)\big)

其中 $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 = \frac{1}{2\beta} \ln \frac{1 - T_2/T_1}{1 + T_2/T_1}

其中 $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 失去意义，便退回一阶（ $\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 读数如何变化，那种平衡便触手可及。

TENO 找到它，THINC 填补它 — BVD 混合#

把 THINC 用到每个网格，连光滑波也会被破坏。所以它必须只在间断网格启动。论文用 TENO 加权已经算出的 $\delta_k$ 构造逐网格指标 $\zeta_i$ 。

\zeta_i = \begin{cases} 2, & \delta_{2,i}=0 \text{ 且 } \delta_{1,i}=0 \\ 1, & \delta_{1,i}=0 \text{ 且 } \delta_{2,i}=1 \\ -1, & \delta_{2,i}=0 \text{ 且 } \delta_{1,i}=1 \\ 0, & \text{其他} \end{cases}

$\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)  # 非单调则退回一阶
 
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

转五圈后，minmod TVD 峰值跌到 0.84 并变窄。THINC-BVD 保持峰值 1.00，宽度也几乎守住初始值（30 个网格）。即便更高阶的 TENO，因多项式的本性也无法完全阻止这种抹平——这正是论文的出发点。

亲手把玩的重构与长时间对流#

操作下面的模拟。它实时显示同一个 top-hat 绕网格运行时，两种方式如何分道扬镳。

β = 1.80

随着 revolutions 增大，红色的 minmod 曲线两端塌陷，铺成梯形。青色的 THINC-BVD 紧贴灰色精确解，守住矩形。把 $\beta$ 降到 1.1，THINC 也像 TVD 一样变缓；在 2.0 附近最为清晰。

复现时心存疑虑之处#

承诺很干净，但实现时有三点绊住了我。

其一，BVD 判定对探测器质量敏感。若把光滑极值（smooth extremum）误判为间断，那个峰会被不自然地压平。上面的简单 BVD 用单调性检查来防范，但多维、方程组系统需要特征变量分解。

其二， $\beta$ 与问题相关。 $\beta = 1.8$ 是为一维基准调出来的值。在强可压缩湍流中，过压缩可能造出人为阶梯，因此自适应 $\beta$ 值得考虑。

其三，成本论断要看上下文。"THINC 几乎免费，因为间断网格很少"这句话，只看重构时是对的。但分支密集的 BVD 判定可能在 GPU 上引发 warp 发散。把它移植到 OpenFOAM 这类非结构 FV 代码，意味着要把逐网格探测逻辑揉进面通量循环，远不如一维结构网格干净。

接下来读什么#

若想了解 THINC 的源头，推荐 Xiao 最初的 VOF 版 THINC；BVD 原理本身则看 Deng 的 P $n$ T $m$ -BVD。自适应 $\beta$ 与多相扩展由 Shyue & Xiao 的 THINC with PE 系列接力。今日一句：别用多项式画跳变，用跳变画跳变。