[论文评述] 当WENO把小涡抹平时 — Fu(2019) TENO低耗散有限体积格式

在斯坦福,Lin Fu 跑湍流直接数值模拟时看见一件怪事。五阶 WENO-JS 把激波抓住了。离激波很远的小涡也一起糊了。阶数写的是五。手里剩下的分辨率更像三阶中心差分。

把 $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$ 是把五阶找回来的最优线性权重。 $\varepsilon$ 在 $10^{-6}$ 附近,用来挡住除以零。这些权重是连续的。模板只轻轻晃一下,本该停在 $w_k=d_k$ 的光滑区里 $w_k$ 也会动,重量漏到更低的阶上。

光滑波走几十步,振幅被削,相位误差攒起来。激波和湍流涡待在同一个 LES 里时,这份耗散是致命的。

TENO的两个思路: 强尺度分离与锐截断#

Fu 改了两处。

模板的组成变了。三个三点模板 $S_0, S_1, S_2$ 旁边,单独放一个大的五点模板 $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$ 的五阶线性式。只有不连续附近才像 ENO 那样挑模板。一个截断能不能扛住所有波,是另一件事。认为每种波都该拆开重构的论点怀疑的就是这次尝试。

从连续权到0/1判定#

$S_3$ 若光滑,也就是 $\delta_3=1$ ,重构就是五阶线性式。

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%。骨架相同,只是权重算法不同。

往权重公式里放一个小晃动。理想线性权重是十分之一、十分之六、十分之三。三个 β 都是 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,大模板活着,五阶线性式原样使用。同样的 10%,WENO 的权重在晃,TENO 一个比特都不翻。

正文公式里的 ε 在 $10^{-6}$ 附近。函数内部是 $10^{-40}$ 。β 在 1 附近时,两者都不出现在结果里。只有 β 贴在 0 上的模板,分母才会分开。跟截断无关,ε 写在哪里,值就不同。五阶系数 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 大多咬住五阶线性式,方波旁边能看见残余振荡。 $C_T$ 是稳定和低耗散之间的秤。

批判性评述 — 调参陷阱

把论文没有大声说死的地方记下来。

$C_T$ 粘在问题上。Fu 推荐 $10^{-5}$ 。强激波和很弱的声波待在一起时,这个值会把声波轻轻抬起,降到 $10^{-6}$ 则在激波附近打转。移植到 OpenFOAM 或 SU2 的小组写道,这个值要按算例重调。

有限体积里, $S_3$ 的光滑指标 $\beta_3$ 很贵。五点多项式的导数积分写满大约六页。写成有限差分就短得多。若要把这篇论文的低耗散 Riemann 通量切换一起用上,有限体积的架子更自然。

和时间积分的配合是空的。SSP-RK3 里,它被空间精度压着,还过得去。和 SSP-RK2 或隐式法混在一起, $\delta_k$ 的间断被事后报告说会碰 Newton 收敛。

拿去用的话,网格收敛用同一个 $C_T$ 来跑。看空间阶是不是真回到 5。回不来, $C_T$ 很可能太大,TENO 一直停在非线性模式。

把C_T固定成一个数时垮掉的验证#

把 $C_T$ 冻在 $10^{-5}$ 再加密网格,光滑单元的 χ 一旦落到截断下面,五阶就没了。上面放进去的 10% 晃动没有越过截断。激波旁边再放一条很弱的波,事情就不一样。

那时阶数检验不是回到 5,而是掉到 1 以下,吃掉阶数的是截断,不是空间算子。同一段代码同时给出 0.79 阶和 2.00 阶的验证问的是,这个算例到底在量什么。 $C_T$ 固定的 TENO,会落到这两个数字中的哪一个。