[论文评述] 用 tanh 把激波画到亚网格

若不是只在间断网格换成 tanh，阶数再高，接触面仍会摊到几十个网格。

精度提到五阶、六阶、八阶，接触间断仍会随时间被抹平。WENO 也好，TENO 也罢，一样。用多项式去画跳变，就会这样。Takagi 等人 2023 年的论文岔开了这条路。光滑处留给高阶 TENO，只有间断网格换成 tanh（THINC）。

把重构开关拨到 Polynomial、THINC、Hybrid，同一组网格平均会分成 S 形、竖直跳变或混合。

β = 1.80overshoot: 0.000

Polynomial 在跳变上下略微越界。THINC 在一个网格里几乎竖直。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（六点、八点）与 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]$ 上的重构函数为

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 读数一起动。

$\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% 在一个网格里收回来。 $\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 = \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，因多项式的本性也无法完全阻止这种抹平——这正是论文的出发点。

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

把 revolutions 往前拨，红色 minmod 和青色 THINC 会在同一段时间里分开跑。

β = 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 代码，意味着要把逐网格探测逻辑揉进面通量循环，远不如一维结构网格干净。

要把 β 钉在 1.8 吗#

五圈之后 THINC-BVD 的峰值是 1.00，宽度 30 格。同样条件下 minmod 是 0.84，23 格。光滑极值也被压成台阶时，这个固定值就过了。按波去怀疑那一个旋钮的，是把 THINC 和 MP5 分开的做法。

用 tanh 画跳变，还是事后扫掉涂开的界面，是另一项决定。重构已经守住接触，后处理锐化就闲着。守不住的话，先看探测器，再把 $\beta$ 调高。

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