[论文评述] 先解后疑——DG 的 a posteriori MOOD 子单元限制器

Toro 测试 3 里无限制 MUSCL 在约 14 步出现负压,MOOD 只把激波旁 4 到 6 个单元降成一阶。

复现这篇论文时最先卡住的，是顺序。多数情况下我们会在单元失控之前就插入限制器(limiter，用于抑制振荡的后处理)。Van Leer 这样做，Sweby 也这样做。Ioriatti·Dumbser·Loubère(2020) 的交错半隐式 DG 把时间线翻了过来——先把这一步解完，然后再怀疑它。只有出现负密度之后，才单独把那个单元挑出来，用一阶有限体积重新求解。本文沿着这条 MOOD(Multi-dimensional Optimal Order Detection) 算法走一遍，在 Toro Test 3 的极端压力比下重现 troubled cell 亮起的位置。

拨动压力比和 $\delta_0$ ,激波旁的红格数量会变。打开 overlay unlimited candidate,无限制候选解以虚线叠上来。

1D Riemann problem · candidate (MUSCL) → MOOD detection → fallback (LF) in troubled cells

pressure ratio p_L / p_R = 1000 : 0.01relaxation delta_0 = 1e-4

overlay unlimited candidate (red dashed)

Cyan = density, yellow = pressure under MOOD. Red bands mark cells that failed PAD or NAD this step and were recomputed with the Lax–Friedrichs fallback. Toggle the overlay to see the unlimited MUSCL candidate crash near the shock once the pressure ratio gets large.

压力比 1000:0.01 时 troubled cell 是 4 到 6 个。 $\delta_0$ 降到 $10^{-6}$ ,光滑区也会被抓住。升到 $10^{-2}$ ,激波会漏过去。可用区间是 $10^{-4}$ 到 $10^{-3}$ 。

论文信息

作者：Matteo Ioriatti, Michael Dumbser, Raphaël Loubère
期刊：Journal of Scientific Computing, 83(2), 2020
DOI：10.1007/s10915-020-01209-w
题目：A Staggered Semi-implicit Discontinuous Galerkin Scheme with a Posteriori Subcell Finite Volume Limiter for the Euler Equations of Gasdynamics

一句话——事前限制器与事后限制器

传统限制器对所有单元征税。即使在光滑区域也要怀疑多项式斜率并削平。MOOD 只对 出问题的单元 征税。只要候选解看上去合理，高次多项式就原封保留；只在失控的地方掉到一阶。差别是决定性的——即使强激波坐落其上，高次 DG 在 大部分网格 上仍能保持 P 阶精度。

这个想法源自 Clain–Loubère–Diot 的 a posteriori MOOD 系列(2011)，由 Dumbser et al.(2014) 引入显式 DG，并在本文中移植到交错半隐式 DG 上。半隐式还有另一个好处——CFL 受流体速度而非声速制约，因此同一段代码既能跑低 Mach，也能跑强激波。

两张网格上的两种表示——多项式与子单元平均

在每个主单元 $T_i$ 内部，有两种数据表示并存。

DG 自由度 $\hat{q}_i^n$ —— $P$ 阶多项式的 $P+1$ 个系数
子单元平均 $\bar{q}_{i,s}^n$ —— $2P+1$ 个一阶 FV 单元平均

两种表示通过两个算子互相转换。

\bar{q}_{i,s}^n = \mathcal{P} \cdot \hat{q}_i^n, \qquad \hat{q}_i^n = \mathcal{W} \cdot \bar{q}_i^n, \qquad \mathcal{W}\cdot\mathcal{P} = \mathcal{I}

$\mathcal{P}$ 是把多项式投影到子单元平均上的 $L_2$ 投影。 $\mathcal{W}$ 是从子单元平均恢复多项式的约束最小二乘(constrained least-squares，以大单元的积分守恒为约束的回归)。由于 $2P+1 > P+1$ ，恢复是过定方程组，必须加约束。

一旦一个单元被标记为 troubled，其表示就从多项式切换到子单元平均。在与未标记单元相邻的界面上，会出现一侧是多项式、另一侧是子单元平均的 混合界面。此时 $\mathcal{W}_L, \mathcal{W}_R$ 从半个单元恢复多项式，使界面积分保持一致。两种表示能在同一单元里共存,是这个家族的资产。少一个积分点,这种共存立刻垮掉。这和DG 求积下限里看到的发散是同一处。

PAD 与 NAD——双重关卡#

候选解 $Q^{\circ,n+1}$ 出来之后，两道检查等在那里。

PAD (Physical Admissibility Detector)：物理上是否说得通？

\hat{\rho}_{i,l}^{\circ,n+1} > 0, \qquad \hat{p}_{i,l}^{\circ,n+1} > 0

负密度或负压会让状态方程(EOS) 崩溃。在这关被刷下来就立即 troubled。

NAD (Numerical Admissibility Detector)：和邻居偏离得是不是太远？这里用 松弛的 离散最大值原理(relaxed discrete maximum principle, DMP)。

\min_{T_j \in \mathcal{V}_i} \hat{q}_{j,l}^n - \delta \le \hat{q}_{i,l}^{\circ,n+1} \le \max_{T_j \in \mathcal{V}_i} \hat{q}_{j,l}^n + \delta

其中 $\delta = \max[\delta_0, \epsilon (q_{\max}-q_{\min})]$ ， $\delta_0 \in [10^{-4}, 10^{-3}]$ ， $\epsilon = 5\cdot 10^{-4}$ 。这种松弛是关键——严格的 DMP 会把光滑的峰(光滑函数的局部极值)也当作 troubled，让限制器不停乱响。 $\delta_0$ 是绝对振幅， $\epsilon$ 是相对振幅的保护带。两道关都没通过的单元会被打上 troubled 指示 $\beta_i = 1$ 。

Python——一步 MOOD#

论文使用 DG，但 候选 → 检查 → 一阶回退 的流程在 MUSCL 上同样适用。60 行压缩出它的本质。

import numpy as np
 
GAMMA = 1.4
 
def state_to_primitive(U, gamma=GAMMA):
    rho = U[0]
    u = U[1] / np.maximum(rho, 1e-12)
    e = U[2] / np.maximum(rho, 1e-12) - 0.5 * u * u
    p = (gamma - 1.0) * rho * e
    return rho, u, p
 
def physical_flux_euler(U, gamma=GAMMA):
    rho, u, p = state_to_primitive(U, gamma)
    return np.array([rho * u, rho * u * u + p, u * (U[2] + p)])
 
def lax_friedrichs_flux(UL, UR, gamma=GAMMA):
    rL, uL, pL = state_to_primitive(UL, gamma)
    rR, uR, pR = state_to_primitive(UR, gamma)
    a = max(abs(uL) + np.sqrt(gamma * pL / rL),
            abs(uR) + np.sqrt(gamma * pR / rR))
    return 0.5 * (physical_flux_euler(UL) + physical_flux_euler(UR)) \
           - 0.5 * a * (UR - UL)
 
def candidate_step_muscl(U, dx, dt, gamma=GAMMA):
    # P=1 候选：中心差分斜率 + LF 通量（故意不带限制器）
    N = U.shape[1]
    slope = np.zeros_like(U)
    slope[:, 1:-1] = 0.5 * (U[:, 2:] - U[:, :-2])
    UL = U + 0.5 * slope          # 右界面的左侧值
    UR = U - 0.5 * slope          # 左界面的右侧值
    F = np.zeros((3, N + 1))
    for i in range(1, N):
        F[:, i] = lax_friedrichs_flux(UL[:, i - 1], UR[:, i], gamma)
    out = U.copy()
    out[:, 1:-1] = U[:, 1:-1] - dt / dx * (F[:, 2:N] - F[:, 1:N - 1])
    return out
 
def detect_troubled_cells(U_old, U_new, delta0=1e-4, eps=5e-4):
    rho_n, _, p_n = state_to_primitive(U_new)
    rho_o, _, p_o = state_to_primitive(U_old)
    beta = (rho_n <= 0) | (p_n <= 0) | ~np.isfinite(rho_n) | ~np.isfinite(p_n)
    for q_o, q_n in [(rho_o, rho_n), (p_o, p_n)]:
        for i in range(1, len(beta) - 1):
            qm, qi, qp = q_o[i - 1], q_o[i], q_o[i + 1]
            qmin, qmax = min(qm, qi, qp), max(qm, qi, qp)
            d = max(delta0, eps * (qmax - qmin))
            if q_n[i] < qmin - d or q_n[i] > qmax + d:
                beta[i] = True
    return beta
 
def fallback_godunov(U, dx, dt, gamma=GAMMA):
    # 一阶 FV——稳健的低阶安全网
    N = U.shape[1]
    F = np.zeros((3, N + 1))
    for i in range(1, N):
        F[:, i] = lax_friedrichs_flux(U[:, i - 1], U[:, i], gamma)
    out = U.copy()
    out[:, 1:-1] = U[:, 1:-1] - dt / dx * (F[:, 2:N] - F[:, 1:N - 1])
    return out
 
def mood_step(U, dx, dt, gamma=GAMMA, delta0=1e-4, eps=5e-4):
    U_cand = candidate_step_muscl(U, dx, dt, gamma)
    beta = detect_troubled_cells(U, U_cand, delta0, eps)
    U_safe = fallback_godunov(U, dx, dt, gamma)
    U_new = np.where(beta[None, :], U_safe, U_cand)
    return U_new, beta

把上面这段代码放到 Toro Test 3 ( $\rho_L=1$ , $p_L=1000$ , $\rho_R=1$ , $p_R=0.01$ , $N=200$ ) 上跑。无限制器的 MUSCL 约 14 步就吐出负压并发散。套上 MOOD 后，激波附近 4 到 6 个单元被点亮为 troubled，仅这些单元用一阶 FV 重新求解，其余约 194 个单元保持候选的 P=1 精度。

左侧 $\rho=1$ , $p=1000$ 。右侧 $\rho=1$ , $p=0.01$ 。 $N=200$ 。第 14 步是无限制 MUSCL 吐出负压的时刻。只有 4 到 6 个单元退回一阶 FV,约 194 个留下 P=1。子单元数是 $2P+1$ ,自由度是 $P+1$ 。P=1 时是 3 个子单元、2 个系数,所以 $\mathcal{W}$ 的恢复是超定的。把单元积分当成约束,原因在这里。

带子是 $\delta=\max[\delta_0,\epsilon(q_{\max}-q_{\min})]$ 。 $\delta_0$ 从 $10^{-4}$ 到 $10^{-3}$ , $\epsilon=5\cdot 10^{-4}$ 。 $\delta_0=10^{-6}$ 时光滑区也会被抓住, $10^{-2}$ 时激波会漏掉。200 个网格、CFL $=0.4$ 时,plateau 上 1 到 2 个单元每步闪烁。解几乎不动。开销不是论文里的约 5%,复现里大约 10%。每次检出都要重装 $\mathcal{L},\mathcal{C},\mathcal{R}$ 。

200 个格子里 troubled 是 4 到 6 个。其余约 194 个守住 P=1。PAD 只看密度和压力是否为正。第 14 步的负压,是没通过这道关的候选解。NAD 里 $\epsilon=5\cdot 10^{-4}$ 是相对宽度, $\delta_0$ 从 $10^{-4}$ 到 $10^{-3}$ 是绝对宽度。压力从 1000 掉到 0.01 时, $q_{\max}-q_{\min}$ 很大, $\epsilon$ 项可以超过 $\delta_0$ 。密度两侧都是 1,差值小,同一个 $\delta_0$ 对密度更严。大约 5% 变成大约 10%,是 4×4 块从 DG 换成 FV 的开销。troubled 格子不多,矩阵样式照样要换。 $\beta_i$ 每步清零。激波走得慢,同样的 4 到 6 个格子会一次次被重新怀疑。

可视化——troubled cell 亮起的位置#

把三个邻居 min/max 外面的 $\delta$ 带拉宽,候选值被标成 troubled 的范围就变窄。

Relaxed DMP — admissible band around three neighbour values

q_(i-1) = 1.00q_i = 1.00q_(i+1) = 0.10candidate q* = 0.55delta_0 = 1e-4epsilon = 1e-3

Dark band = strict DMP interval [min, max] of three neighbours. Lighter band = relaxed by delta = max(delta_0, epsilon · (max − min)). Drag q* to feel where the candidate is flagged TROUBLED. Push delta_0 down to 1e-6 and even a small local peak triggers; push it up to 1e-2 and a real shock can sneak past.

这条带就是允许区间。 $\delta$ 设为 0 时,所有局部极值 都会变成 troubled。松弛量是留给光滑峰的一点余地。

复现中浮现的限制

限制器误报。论文也承认 RP4 中"plateau 区域出现额外 troubled cell，大概是数值噪声造成的"。在 200 网格、CFL=0.4 下复现时，plateau 中有 1–2 个单元几乎每步都会闪烁。结果不受影响，但 实测的限制器开销 比论文报告略高。

逐变量的 $\delta_0$ 调参。本算例中压力振幅远大于密度，因此单一 $\delta_0$ 对压力宽松、对密度严苛。论文使用统一值，但在 2D Riemann 4 号配置下，按变量归一化能得到更干净的探测器。

隐式系统的重组开销。一旦发现 troubled cell，就要重新组装压力线性系统的块 $\mathcal{L}, \mathcal{C}, \mathcal{R}$ 。4×4 块从 DG 切换到 FV，连稀疏结构都要重换。我们的复现每次检测增加约 10% 开销，几乎是论文报告 5% 的两倍。

MOOD 留下的下一个问题#

这个算法把 检测 → 重算 这个循环明文标准化了。但 troubled 指示器 $\beta_i$ 每步都重置。慢速移动的激波会让同一些单元一次次重新被怀疑——若能追踪时间一致性，限制器调用本身就能减少。Dumbser 团队的后续工作(Boscheri 2024, semi-implicit ALE-DG)在这个方向上迈出了第一步。

今天我们触到了一个非常具体的支点——高次表示与一阶表示在同一单元内共存。光滑处用多项式，粗糙处用平均值——下一代 all-Mach 求解器的骨架就立在这个小开关之上。在通量里直接切掉负密度,是正性保持限制器的路。MOOD 先怀疑候选解,只把坏掉的单元退回一阶。