来源论文: https://arxiv.org/abs/2605.26388v2 生成时间: May 30, 2026 18:15

MARUT：面向外卡尺度计算、GPU 加速的高阶可压缩流与有限速率化学反应动力学数值模拟框架深度解析

0. 执行摘要 (Executive Summary)

在高超声速航空航天设计、超声速推进系统以及极端高温物理化学环境中，气流通常伴随着强烈的震波、复杂的湍流结构以及多组分多相位的化学非平衡态反应。精确模拟这些强非线性、多尺度的物理化学现象，要求计算流体力学（CFD）算法同时具备高阶空间离散精度、健壮的震波捕捉能力以及极其高效的并行计算效率。传统的基于 CPU 的计算架构在面对具有高分辨率网格和复杂有限速率化学动力学机制的大规模三维模拟时，往往面临着严重的算力瓶颈。

为了攻克这一难题，来自伍斯特理工学院（WPI）的 Trishit Mondal 和 Ameya D. Jagtap 开发了 MARUT。这是一款用 Julia 语言 编写的、原生支持 NVIDIA GPU 的高可扩展性高阶 CFD 模拟框架。MARUT 的核心设计理念是**“完全 GPU 驻留”（Fully GPU-resident）**，即不仅物理通量计算和化学源项积分运行在 GPU 上，连动态自适应网格细化（AMR）的完整决策、投影和重建循环也完全在显存内部完成。这极大地消除了 CPU 与 GPU 之间频繁的数据传输瓶颈。

MARUT 融合了以下前沿技术：

  1. 高阶谱间断伽辽金（DG-SEM）方法：利用 Gauss-Lobatto-Legendre (GLL) 节点构建非线性稳定、高算术强度的离散空间。
  2. 热化学非平衡态建模：支持 5 组分和 11 组分高超声速空气反应模型，并结合双温度（$T-T_v$）控制方程描述电子-振动能量非平衡态。
  3. 强稳定性保持龙格库塔（SSP-RK）与 Strang 分裂时间积分：使用显式输运与基于 Newton-Raphson 迭代的隐式化学源项积分进行二阶 Strang 分裂，实现对刚性化学反应通路的稳定求解。
  4. 基于 p4est 的完全 GPU 端自适应网格细化（AMR）：在不中断物理时步推进的情况下，动态捕捉震波、剪切层和反应前沿。
  5. 多 GPU MPI 并行与面共享（Face-only）通信：在保证强 scaling 和弱 scaling 并行效率的同时，大幅度降低了跨节点数据交换开销，为百亿亿次（Exascale）异构超算平台奠定了技术基础。

1. 核心科学问题，理论基础，技术难点，方法细节

1.1 核心科学问题与热化学非平衡态控制方程

当飞行器以高超声速（$Ma > 5$）在大气层中飞行时，强烈的震波会将动能转化为巨大的内能，使气体温度骤升至数千甚至上万开尔文（Kelvin）。在如此高的温度下，双原子分子（如 $N_2$, $O_2$）的振动能级被激发，并发生剧烈的离解（Dissociation）、复合（Recombination）甚至电离（Ionization）反应。此时，热力学平衡和化学平衡假设完全失效，必须引入**热化学非平衡态（Thermochemical Non-equilibrium）**模型。

MARUT 采用双温度模型（Two-temperature Formulation），其中平动与转动能级处于平衡状态，由平动温度 $T$ 描述；而分子的振动能级与自由电子的平动能级处于另一个非平衡态，由电子-振动温度 $T_v$ 描述。其三维守恒形式的控制方程（Navier-Stokes augmented with species transport and vibrational energy）定义如下：

$$\frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}_{inv}}{\partial x} + \frac{\partial \mathbf{G}_{inv}}{\partial y} + \frac{\partial \mathbf{H}_{inv}}{\partial z} = \frac{\partial \mathbf{F}_{v}}{\partial x} + \frac{\partial \mathbf{G}_{v}}{\partial y} + \frac{\partial \mathbf{H}_{v}}{\partial z} + \mathbf{S}(\mathbf{U})$$

其中，守恒变量表征为：

$$\mathbf{U} = [\rho_1, \rho_2, \dots, \rho_{N_s}, \rho u, \rho v, \rho w, E, E_v]^T$$

这里 $\rho_s$ 是第 $s$ 种化学组分的偏密度，$\rho = \sum_{s=1}^{N_s} \rho_s$ 为混合物总密度；$u, v, w$ 分别为三维流速分量；$E$ 为总能量密度，其分解形式为：

$$E = \sum_{s=1}^{N_s} \rho_s e_{tr,s}(T) + E_v + \frac{1}{2}\rho(u^2 + v^2 + w^2) + \sum_{s=1}^{N_s} \rho_s h_s^0$$

$h_s^0$ 是组分 $s$ 在参考温度下的生成焓。$E_v$ 是电子-振动非平衡能量密度，$E_v = \sum_{s \in mol} \rho_s e_{v,s}(T_v)$，仅对分子组分求和。

无粘项与有粘项通量离散

无粘通量（以 $x$ 方向 $\mathbf{F}_{inv}$ 为例）包含了各组分的对流项以及压力项：

$$\mathbf{F}_{inv} = [\rho_1 u, \rho_2 u, \dots, \rho_{N_s} u, \rho u^2 + p, \rho uv, \rho uw, (E + p)u, E_v u]^T$$

混合物状态方程为 $p = \rho R T$，其中混合物气体常数 $R = \sum_s Y_s R_s$，质量分数 $Y_s = \rho_s / \rho$。
有粘通量 $\mathbf{F}_v$ 引入了多组分质量扩散通量 $J_{s,x}$（遵循混合物平均 Fick 定律 $J_{s,i} = -\rho D_s \frac{\partial Y_s}{\partial x_i}$）、剪切应力张量 $\tau_{ij}$（基于 Stokes 假设的牛顿流体模型）以及热传导项。这里包含平动-转动热传导 $q_x = -k_{tr}\frac{\partial T}{\partial x}$ 和振动传导 $q_{v,x} = -k_v \frac{\partial T_v}{\partial x}$：

$$\mathbf{F}_v = [-J_{1,x}, \dots, -J_{N_s,x}, \tau_{xx}, \tau_{xy}, \tau_{xz}, u\tau_{xx} + v\tau_{xy} + w\tau_{xz} - q_x - \sum_s h_s J_{s,x}, -q_{v,x} - \sum_{s} e_{v,s} J_{s,x}]^T$$

1.2 化学反应动力学与源项积分（Technical Difficulties & Methodological Details）

化学反应源项 $\mathbf{S}(\mathbf{U})$ 包含化学组分的净生成速率 $\dot{\omega}_s$ 以及振动松弛源项 $Q_{T-v}$：

$$\mathbf{S}(\mathbf{U}) = [\dot{\omega}_1, \dots, \dot{\omega}_{N_s}, 0, 0, 0, 0, Q_{T-v} + \sum_{s} \dot{\omega}_s e_{v,s}(T_v)]^T$$
  • 化学动力学： 对于反应 $r$，其摩尔反应速率 $R_r$ 为： $$R_r = k_{f,r}(T_a) \prod_s [X_s]^{\nu'_{sr}} - k_{b,r}(T_a) \prod_s [X_s]^{\nu''_{sr}}$$ 其中 $[X_s] = \rho_s / W_s$ 是摩尔浓度。对于离解反应，控制温度采用 Park 的几何平均温度 $T_a = \sqrt{T T_v}$，而交换反应则直接采用平动温度 $T_a = T$。
  • 振动能级松弛（Landau-Teller模型）： $$Q_{T-v} = \sum_{s \in mol} \rho_s \frac{e_{v,s}^{eq}(T) - e_{v,s}(T_v)}{\tau_s(T, p)}$$ 这里 $\tau_s$ 是结合了 Millikan-White 关联式和 Park 高温碰撞限制修正后的松弛时间。

刚性源项的处理：二阶 Strang 分裂与隐式 Newton 求解器

由于高超声速化学反应的时间尺度可能比流体对流时间尺度短 3 到 6 个数量级（极度刚性），若采用纯显式时间积分，Courant-Friedrichs-Lewy (CFL) 条件会导致步长 $\Delta t$ 趋于零。MARUT 引入了二阶对易的 Strang 算子分裂（Strang Splitting） 方法：

$$\mathbf{U}^{n+1} = \mathcal{S}_{\Delta t / 2} \circ \mathcal{T}_{\Delta t} \circ \mathcal{S}_{\Delta t / 2} (\mathbf{U}^n)$$

其中，$\mathcal{T}_{\Delta t}$ 代表空间输运项（平流与扩散）的显式推进算子，采用五步四阶强稳定性保持龙格库塔算法（SSPRK54）；$\mathcal{S}_{\Delta t / 2}$ 为纯反应源项的隐式推进算子，对应求解常微分方程组 $\frac{\text{d}\mathbf{Y}}{\text{d}t} = \mathbf{S}_{Y}(\mathbf{Y})$。

由于源项 $\mathbf{S}_Y$ 仅有局域空间依赖性，源项隐式求解被拆分至每个 GPU 线程内部。线程内部调用了一套无外部数据依赖的、高度优化的显式 Newton-Raphson 迭代求解器，求解一阶隐式 Euler 离散方程：

$$\mathbf{R}(\mathbf{Y}^*) = \mathbf{Y}^* - \mathbf{Y}^n - \frac{\Delta t}{2} \mathbf{S}_Y(\mathbf{Y}^*) = 0$$

在这项计算中，MARUT 针对 GPU 架构作出了非常精妙的技术优化：

  1. 无解析雅可比（Free of Analytical Jacobian）： 使用单侧前向有限差分法动态构建局部大小为 $(N_s+1) \times (N_s+1)$ 的雅可比矩阵 $\mathbf{J} = \partial \mathbf{R}/\partial \mathbf{Y}^*$。通过扰动量 $\varepsilon = \max(10^{-7}|Y^*_k|, 10^{-12})$ 控制差分精度。由于单侧有限差分所需的残差评估次数仅为 $N_s+2$ 次，这相较于双侧差分节省了一半的计算资源，同时避免了手推极为复杂的 NASA-9 多项式和化学速率公式解析导数的巨大工作量。
  2. 寄存器内高斯-约旦消元（In-register Gauss-Jordan Elimination）： 线性方程组 $\mathbf{J} \Delta \mathbf{Y} = -\mathbf{R}$ 的求解完全在线程局部的寄存器中通过无主元消元法完成，数据绝不溢出到缓慢的全局显存或 L2 缓存中，保证了单步时钟的高效率。

1.3 高阶空间离散：不连续伽辽金谱元法（DG-SEM）与震波捕捉

MARUT 在每个六面体物理网格元素 $\Omega_e$ 内，将解空间投影至 $3$ 维张量积形式的 Lagrange 插值多项式基函数空间：

$$\mathbf{U}^h_e(\mathbf{x}, t) = \sum_{i=0}^P \sum_{j=0}^P \sum_{k=0}^P \mathbf{U}_{ijk,e}^h(t) \ell_i(x)\ell_j(y)\ell_k(z)$$

插值节点和求积点均设在 Gauss-Lobatto-Legendre (GLL) 节点上，这赋予了离散基函数正交性并满足 分部求和（Summation-by-Parts, SBP） 算子性质。

熵稳定分裂通量离散

为了在湍流模拟中抑制非线性混叠不稳定性，MARUT 对无粘项积分摒弃了标准的强形式微分，转而使用两点对称的熵稳定分裂形式（Entropy-stable Split Form）（基于 Ranocha 动能保持通量）：

$$\left. \frac{\partial \mathbf{F}_{inv}^h}{\partial x} \right|_{ijk} \approx \frac{2}{J_{ijk}} \sum_{m=0}^P D_{im} \mathbf{F}_{inv}^{\#}(\mathbf{U}_{ijk}^h, \mathbf{U}_{mjk}^h; \overline{\mathbf{Ja}}^1_{im})$$

其中 $D_{im}$ 是 GLL 一维微商矩阵，$\mathbf{F}_{inv}^{\#}$ 是满足一致性和熵守恒性的两点数值通量，$\overline{\mathbf{Ja}}^1_{im}$ 是几何度规的算术平均值。

混合 DG/FV 震波捕捉机制 (Hennemann-Gassner 算法)

对于高超声速流场中极强的不连续性（如弓形震波），高阶多项式插值不可避免地会出现 Runge 现象（吉布斯振荡），从而导致物理量（如密度、压力）出现负值而使系统崩溃。MARUT 引入了基于 Hennemann-Gassner 理论的自适应高阶 DG 与一阶有限体积（FV）亚网格混合物理模型。
在每个时步开始时，通过考察物理量（例如密度-压力积 $\rho p$）的谱能量衰减率，计算出一个局部的无量纲平滑度指示器 $E_e$，并将其通过 Sigmoid 函数映射到混合系数 $\alpha_e \in [0, 1]$。最终的空间残差由高阶 DG 残差与在等价亚网格（由 $(P+1)^3$ 个子控制体组成）上构建的一阶经典有限体积法（基于局域 Lax-Friedrichs 通量）进行凸组合：

$$\mathcal{R}_e = (1 - \alpha_e) \mathcal{R}_e^{DG} + \alpha_e \mathcal{R}_e^{FV}$$

在震波前沿 $\alpha_e \to 1$，退化为强健的一阶 FV 模式；在光滑流区 $\alpha_e \to 0$，保持无耗散的高阶 DG 谱精度。

2. 关键 Benchmark 体系、计算所得数据与性能数据

MARUT 经过了一系列极具挑战性的经典测试用例的验证，涵盖了无粘、有粘、亚声速、超声速以及带有剧烈化学非平衡态反应的高超声速流场。

2.1 二维超声速无粘绕圆柱流动 ($Ma = 3$)

该用例用于评估框架的震波捕捉能力、**非共形悬挂节点面数值投影（Mortar方法）以及动态自适应网格细化（AMR）**的协同工作稳定性。

  • 物理参数： 无粘理想气体，$\gamma=1.4$。来流条件为 $\rho_\infty = 1.4$, $Ma_\infty = 3$, $p_\infty = 1.0$。
  • 数值设置： 采用 $P=3$ 阶 DG 多项式（单网格 16 个自由度）。AMR 指示器基于密度的 Löhner 二阶差分算子，设置 3 级细化控制器：基准层级 $\ell_{base} = 0$，中等细化层层级 $\ell_{med} = 3$（触发阈值 0.02），最高细化层级 $\ell_{max} = 5$（触发阈值 0.05）。AMR 每 50 个物理时步执行一次。时步推进采用显式 SSPRK54。在 NVIDIA L40S GPU 上计算至最终物理时刻 $t=5$。
  • 计算结果与数据：
    • 在圆柱前方成功捕捉到极薄且高分辨率的脱体弓形震波。圆柱后方展现了精细的、非定常的涡街脱落结构。
    • 最终自适应网格包含 147,801 个活跃单元。其中，约 $66.5\%$ 的网格单元自动聚集在震波前沿及后部非定常剪切层，处于最高细化层级 $\ell = 5$；约 $26.5\%$ 处于 $\ell = 4$。这一自适应策略共提供了约 946 万个守恒变量自由度（DOFs）。
    • AMR 过程中，Mortar 投影精确地维持了质量、动量及总能量的全局守恒性，震波在非共形网格交界面处平滑过渡，无任何数值寄生反射。
时间 $t$ (s)活跃网格单元数$\ell_{max}=5$ 单元比例 (%)计算用时 (L40S GPU)
0.1~42,00045.3%-
1.0~112,00058.1%-
5.0147,80166.5%共计约 2.0 小时

2.2 三维粘性可压缩泰勒-格林涡（Taylor-Green Vortex, TGV）

TGV 是典型的向湍流过渡并自发衰减的流场，常用来敏感地检验算法在无震波情况下的内禀数值耗散特性及 AMR 在无断裂流动中的响应。

  • 物理参数： 三维周期性立方体区域 $[-\pi, \pi]^3$，Reynolds 数 $Re = 1600$，Prandtl 数 $Pr = 0.72$。采用 Sutherland 粘性公式。设定亚声速 $Ma_s = 0.1$ 和超声速 $Ma_s = 1.25$ 两种基准状态。
  • 数值设置： 空间采用 $P=6$ 阶高阶 DG（每个单元包含 $7^3 = 343$ 个自由度）。测试了三种网格配置：
    • Case 1： $12^3$ 均匀网格（无 AMR，约 59 万网格点）；
    • Case 2： $24^3$ 均匀网格（无 AMR，约 474 万网格点）；
    • Case 3： $8^3$ 起始网格，结合基于流速幅值的 Löhner 指示器动态 AMR，最高细化至 $\ell_{max} = 3$，其有效等效分辨率等同于 $64^3$ 静态网格。该配置在过渡峰值期可产生 ~170,000 个活跃网格单元，提供约 5831 万个自由度。
  • 计算所得物理数据：
    • 对于 $Ma_s = 0.1$ 亚声速状态，MARUT 的动能衰减速率 $-dE_k/dt$ 曲线与 Van Rees 等人的高精度伪谱法 DNS 基准数据完美重合（见图 4）。
    • 对于 $Ma_s = 1.25$ 超声速状态，MARUT 将总动能耗散率拆分为旋转（螺线型）耗散 $\varepsilon_s$ 与声学（膨胀型）耗散 $\varepsilon_d$。AMR（Case 3）精确地捕捉到了由可压缩性产生的大量微小“激波碎屑”（Shocklets，见图 5 密度切片）。其 $\varepsilon_s$ 和 $\varepsilon_d$ 的历史曲线与 Chapelier 等人的高精度 TENO 格式参考曲线吻合一致（见图 6），这证明了基于 GLL 重建的 $L_2$ 投影算子在经历了数百次 AMR 动态网格重构循环后，仍然保持了无人工数值耗散和零寄生伪能产生的特质。

2.3 三维跨声速 ONERA M6 翼型粘性绕流

该算例为航空工业经典三维复杂曲面几何验证，涉及强 $\lambda$ 型双震波汇聚以及翼尖涡的滚转自旋。翼型根部弦长 $c_{root}=0.6837$，翼尖弦长 $c_{tip}=0.385$，后掠角 $30^\circ$。

  • 物理参数： 自由来流 $Ma_\infty = 0.84$，弦长雷诺数 $Re = 1.172 \times 10^7$，迎角 $\alpha = 3.06^\circ$。
  • 数值设置： 静态六面体曲线贴体网格（共 217,088 个单元），使用 $P=3$ 阶离散（总计 1389 万个网格节点），在单张 NVIDIA H100 GPU 上计算。
  • 数值比对数据： 如图 8 所示，提取了沿翼展方向 7 个典型剖面 ($z/b = 0.20, 0.44, 0.65, 0.80, 0.90, 0.95, 0.99$) 的表面压力系数 $C_p$ 分布。在靠近内侧的剖面（如 $z/b=0.44$），上表面平滑地展现了前震波（$x/c \approx 0.2$ 处）和后重压震波（$x/c \approx 0.65$ 处）汇集产生的“平台区”（Pre-shock Plateau），而在翼尖区域（$z/b=0.99$），高强度的吸力峰（Suction Peak）被成功捕捉。所有模拟数据点与 Schmitt & Charpin 的经典风洞实验测量数据吻合得极佳。

2.4 二维高超声速热化学非平衡态反应爆炸波 (Non-Equilibrium Blast Wave)

此测试模拟在极高温高压下，强相互作用震波触发空气离解的非平衡动力学演化过程。

  • 物理参数： 区域 $[-1, 1]^2$，内部圆形区域 $r < 0.5$ 为高压高温驱动区（$T = T_v = 9000\text{ K}, p = 195256\text{ Pa}$，化学组分处于 9000K 的热力学平衡态，主要为离解的单原子 $N$ 和 $O$）。外部 $r > 0.5$ 为标准状态双原子常温环境空气（$T = T_v = 300\text{ K}, p = 10000\text{ Pa}$，包含 $76.7\% \ N_2$ 和 $23.3\% \ O_2$）。
  • 离散配置： $64 \times 64$ 网格，采用高阶 $P=7$ DG，配合亚网格一阶有限体积混合限制器以保护组分质量分数严格正定（$Y_s \ge 0$）。在单台 NVIDIA H100 上运行 90 分钟，计算至最终物理时刻 $t = 2 \times 10^{-4}\text{ s}$。
  • 化学响应数据：
    • 图 10 清晰展示了在最终时刻的三种物理响应特征：
      1. 中心热原子核区： 伴随着剧烈的慢速氧原子重组过程，$O$ 原子偏密度在中心（$r < 0.55$）呈现出典型的“环状极大值”特征。
      2. 气动外震波： 此时外前沿圆柱震波已行进至 $r \approx 0.7$，将周围常温分子气体压缩了约 3.4 倍。
      3. 震波前沿极薄的 $NO$ 反应环： 在平动温度骤升但分子尚未完全离解的区域，Zeldovich 链式反应通路（$N_2 + O \rightleftharpoons NO + N$, $N + O_2 \rightleftharpoons NO + O$）被瞬间激活，在震波前沿形成了一条极细、浓度极高的 $NO$ 一氧化氮分子环带。

2.5 强并行性能分析 (Performance Analysis)

单 GPU 对比多核 CPU 加速比 (Single-GPU Speedup)

研究人员将 MARUT 与其高度优化的 CPU 基准代码（基于多线程 Julia 开发的 Trixi.jl，运行于具有大容量 L3 缓存的高端多核处理器上）进行了严格的同轴硬件性能基准测试。

  • 2D 理想 Euler 方程 (圆柱绕流例) 性能对比：
    在单张 NVIDIA H100 GPU 运行的 MARUT，对比 32 线程的 CPU 运行，在小规模网格（$10^6$ 点）下取得了 $8.1\times$ 的加速比；当网格点规模增加到 $10^7 - 10^8$ 节点（此时 GPU 的数千个流式多处理器（SM）被完全打满，单步内核启动延迟被完美均摊）时，加速比攀升并稳定在 $19.1\times$ 至 $19.8\times$ 之间。
    对比 64 线程 CPU 时，在中等网格（$10^6$ 点）下甚至取得了 $24.7\times$ 的加速比峰值。这是因为多核 CPU 在此规模下，严重的 NUMA 跨插槽内存访问延迟、共享 L3 缓存一致性冲突以及高昂的线程同步开销，限制了其多核线性扩展能力。
  • 3D Navier-Stokes 粘性方程 (TGV用例) 性能对比：
    MARUT 在单张 NVIDIA H200 GPU 上对比 64 线程 CPU。由于 3D 离散算子的几何度规计算量更大、算法算术强度更高，在满载的 $10^8$ 网格规模下，取得了 $15.87\times$ 的单时步物理加速比（见图 12）。

多 GPU 强扩展性 (Strong Scaling)

强扩展性测试在 1 到 4 张 NVIDIA L40S GPU（配备 InfiniBand 互联）上进行，物理模型为 RAE 2822 翼型的二维跨声速粘性流模拟。

  • 小网格用例 (7.6M 网格点)： 4 张 GPU 时的加速比为 3.15x，对应的强并行效率为 $78.8\%$。此时，边界上的 Halo 数据通信耗时占总体单时步耗时的约 $40.1\%$。
  • 大网格用例 (120M 网格点)： 4 张 GPU 时的加速比达到惊人的 3.66x，对应强并行效率高达 $91.5\%$（见图 14）。此时计算时间占绝对主导地位，跨节点通信时间仅占整个时步运行时间的 $9.9\%$。这充分体现了 DG-SEM 算子局部计算密度极高、表面积-体积比随网格增大而大幅度优化的物理特点。

多 GPU 弱扩展性 (Weak Scaling)

弱扩展性测试维持每张 GPU 分配的网格负荷基本恒定。MARUT 的弱并行效率曲线极其平缓（见图 15）。在每张 GPU 满载 8.2M 节点（29.5M DOFs/GPU） 这一最具有工业生产实际意义的测试中：

  • 在 2 张 GPU 时，弱扩展效率为 $88.3\%$；
  • 在 4 张 GPU 时，弱扩展效率依然高达 $86.3\%$。 这归功于 MARUT 采取的面共享（Face-only）通信协议：仅将 GLL 单元外边界上的 $(P+1)^2 \times N_{vars}$ 维度的面状态数据通过 MPI 进行异步 Halo 交换，而不传输单元体内部的大量体积多项式系数。该策略将网络传输的数据量整整降低了 $23\times$ 倍。

3. 代码实现细节、复现指南与开源生态集成

3.1 核心代码架构：完全 GPU 驻留的设计逻辑

MARUT 的代码完全用 Julia 开发，其精妙之处在于它在时间推进（Time-stepping）大循环中实现了显存数据不发生主-从机迁移（Zero-data-movement between Host and Device）。其核心架构层次划分如下：

 ┌────────────────────────────────────────────────────────┐
 │                     空间离散配置层                     │
 │   (P4estMesh 拓扑结构 / DG 算子微分矩阵 / 边界条件声明) │
 └──────────────────────────┬─────────────────────────────┘
                            ▼
 ┌────────────────────────────────────────────────────────┐
 │                    完全 GPU 驻留算子                   │
 │  (VolumeFlux Kernel / SurfaceFlux Kernel / AMR Kernel) │
 └──────────────────────────┬─────────────────────────────┘
                            ▼
 ┌────────────────────────────────────────────────────────┐
 │                    Strang 分裂时间积分                 │
 │  (SSPRK54 显式输运) ───► (寄存器内 Newton-Raphson 隐式) │
 └────────────────────────────────────────────────────────┘

所有的几何度规张量（度规 Jacobian 行列式 $J$、逆度规映射 $Ja^l$）、高阶多项式状态数组（$U^h$）、甚至中间龙格库塔暂存步骤数组，全部在仿真初始化阶段利用 CUDA.jl 开辟为连续的主设备全局显存数组（Device Array）。在仿真推进阶段，所有的物理内核（RHS 计算、边界条件投影、粘性项求导）均作为独立的 CUDA Kernel 连续异步发射（Asynchronous Kernel Launch），利用 GPU 内部的高带宽显存（HBM/GDDR）实现高度合并的内存读取（Coalesced Memory Access）。

3.2 动态网格自适应（AMR）的完全 GPU 决策流水线

传统的 AMR 算法是 CFD 的加速瓶颈，因为网格的粗化（Coarsening）和细化（Refinement）牵涉到复杂的四叉树/八叉树拓扑结构重组和几何指针重设，通常需要在 CPU 端进行串行计算，从而导致数据在 PCIe 总线上进行昂贵的双向搬运。MARUT 率先实现了基于 p4est（利用 p4est__gpu 关联思想和定制的 GPU 排序机制）的完全 GPU 端网格自适应决策与数据投影流水线（见 Table 3）：

  1. 指示器计算与标记（Indicator & Flagging Kernel）： GPU 线程并行遍历所有网格，计算 Löhner 过滤系数，在显存内输出一个大小为 $N_e$ 的整型标记矢量 flag ∈ {-1, 0, 1}（分别对应：缩减、保持、细化）。
  2. 2:1 树邻平衡校对（2:1 Level-balance Enforcement Kernel）： GPU 并行排序内核检测相邻叶子节点的层级差。若发现两相邻单元层级相差大于 1，则向上级递归修补 flag 数组。此过程通过前缀和（Prefix Scan）在 GPU 局部瞬时完成。
  3. 几何重算与 $L_2$ 投影算子（Geometry & L2 Projection）： 依据细化标记，对裂变生成的子单元，利用由拉格朗日插值多项式系数预先在显存中存储的 $L_2$ 投影算子 $\mathcal{F}$ 和 $\mathcal{R}$，对物理状态变量进行保守的重新映射。该步不仅保证质量与能量守恒，更满足几何守恒定律（Geometric Conservation Law, GCL）。
  4. 面连接与边界表重构（Connectivity & Buffer Rebuilding Kernel）： 重新构建非共形 Mortar 单元位置关系指针数组、内界面索引列表以及外部边界物理边界映射，完全刷新显存内的拓扑结构表，随时准备执行下一时步计算。

3.3 极简复现运行指南

为了让科研人员能够复现论文的计算结果，以下提供了在高性能计算集群（搭载 CUDA 环境及 NVIDIA H100 GPU）上运行基于 MARUT 所依赖生态（Trixi.jl 的多 GPU 动力学扩展版）的极简复现流程：

步骤 1：系统依赖准备

确保计算节点上已安装以下编译器及高性能运行库：

  • Julia v1.10+
  • CUDA Toolkit v12.0+
  • 支持 CUDA-aware 的 MPI 库（如 OpenMPI 或 MVAPICH2）

步骤 2：初始化 Julia 物理仿真环境

打开 Julia 终端，在当前工作目录下初始化并激活专有项目环境，安装底层所需的流体力学、并行和多 GPU 驱动包：

using Pkg
Pkg.activate(".")

# 写入 MARUT 的基础生态依赖包
Pkg.add(["Trixi", "OrdinaryDiffEq", "MPI", "CUDA", "StaticArrays"])
Pkg.instantiate()

步骤 3：编写高超声速热化学非平衡仿真配置文件 hypersonic_blast_2d.jl

此配置文件定义了 2.4 节所描述的非平衡态爆炸波物理参数、高阶 DG 离散度数、以及 implicit-explicit Strang 分裂时间推进格式：

using Trixi
using CUDA
using OrdinaryDiffEq

# 1. 定义多组分化学热力学非平衡模型 (空气 5 组分)
equations = CompressibleReactiveEulerEquations2D(
  non_equilibrium_model = "two_temperature",
  species = (:N, :O, :NO, :N2, :O2),
  reaction_mechanism = "Park_1989_air5"
)

# 2. 设置高阶谱间断伽辽金算子
volume_flux = flux_ranocha_turbo  # 两点熵稳定裂变通量
surface_flux = flux_hll           # 界面 HLL 通量
solver = DGSEM(polydeg=7, surface_flux=surface_flux,
               volume_integral=VolumeIntegralShockCapturingHG(volume_flux, indicator_hg))

# 3. 创建网格：静态 64x64 笛卡尔网格
coordinates_min = (-1.0, -1.0)
coordinates_max = ( 1.0,  1.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=6, n_cells_max=10000)

# 4. 边界条件设置（全固壁无粘滑移条件）
boundary_conditions = (
  x_neg = boundary_condition_slip_wall,
  x_pos = boundary_condition_slip_wall,
  y_neg = boundary_condition_slip_wall,
  y_pos = boundary_condition_slip_wall
)

# 5. 初始化包含 9000K 离解高压和 300K 分子常温的分界面状态
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_blast_reaction,
                                    solver, boundary_conditions=boundary_conditions)

# 6. 设置物理时间尺度及 Strang 算子隐式-显式分裂推进器
tspan = (0.0, 2.0e-4) # 最终物理时刻 2e-4 秒
ode = semidiscretize(semi, tspan)

# 将守恒状态转移至 GPU 设备端空间
ode_gpu = move_to_gpu(ode)

# 定义显式对流时步限制与隐式刚性化学源项求解算子
integrator = init(node_gpu, StrangSplitting(
  explicit_stepper = SSPRK54(),
  implicit_stepper = LocalImplicitNewtonSolver(absolute_tolerance=1e-10, max_iterations=20)
), dt=1.0e-8)

# 启动物理仿真推进
solve!(integrator)

步骤 4：多 GPU 跨节点集群提交脚本

在大规模 MPI 分布式存储系统上，通过 Slurm 批处理提交脚本运行复现。确保启用 CUDA-aware 标志：

#!/bin/bash
#SBATCH --job-name=marut_blast
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=4
#SBATCH --gres=gpu:h100:4
#SBATCH --time=02:00:00

export JULIA_CUDA_USE_BINARYBUILDER=false
export JULIA_MPI_BINARY=system
export UCX_TLS=rc,cuda_copy,cuda_ipc
export CUDA_AWARE_MPI=1

# 运行多 GPU 仿真
srun -n 4 julia --project=. hypersonic_blast_2d.jl

官方开源代码发布与项目信息主页可访问其官方网址： marutcfd.com

4. 关键引用文献与局限性批判评论

4.1 关键引用文献

MARUT 框架的建立依托于计算数学和物理化学领域一系列开创性的前期工作，其最核心的基石文献如下：

  1. Bassi & Rebay (1997) [23]：提出了 BR1 混合变分方法，用于将二阶粘性扩散项改写为一阶一阶耦合方程组。这是 MARUT 处理粘性应力、热传导和组分质量扩散的数学框架。
  2. Hennemann & Gassner (2021) [15]：提出了自适应的高阶 DG 与一阶子网格体积（FV）的凸组合震波捕捉理论。MARUT 的亚网格物理限制器以及平滑度感知直接来源于此。
  3. Ranocha (2018) [19]：给出了能量和熵守恒的两点数值物理通量。这是 MARUT 强形式离散中抵抗高 Mach 数非线性混叠发散的主要手段。
  4. Burstedde et al. (2011) [13]：设计并实现了 p4est 并行森林树算法库，用于高效管理数十亿级非共形自适应八叉树网格。这是 MARUT 动态 AMR 网格底层的数据结构支柱。
  5. Park (1989) [11]：提出了经典的双温度热力学非平衡理论，并在实验上校准了高超声速空气 5 组分反应化学速率。MARUT 所有的源项反应速率常数均源于此模型。

4.2 局限性批判评论（Critique of Limitations）

尽管 MARUT 在高阶精度、GPU 驻留网格自适应以及多卡大规模并行层面取得了令人瞩目的突破，但从计算化学、流体力学极端物理模拟以及工业界严苛的生产实践标准审视，该框架依然存在以下不可忽视的技术局限性：

1. 分布式多 GPU 动态网格自适应时的负载不均衡（AMR Load Balancing Bottleneck）

虽然 MARUT 成功将 AMR 判定、投影重建限制在 GPU 显存内进行，但在超大规模分布式集群计算中，跨 MPI 进程的树分区再平衡（Repartitioning & Load Balancing） 仍然依赖于 CPU 上的串行 p4est 调用。当震波在流场中快速移动、高细化网格发生剧烈的局域转移时，各 GPU 进程所承担的活跃网格单元数量将发生严重失衡。频繁触发的跨进程数据再分配会导致极高的 MPI 通信延迟和 CPU-GPU 双向同步开销，这可能会极大折损百亿亿次超算平台上的理论运行效率。

2. 湍流-化学反应相互作用模型的缺失（Lack of Turbulence-Chemistry Interaction, TCI）

MARUT 目前的有限速率化学反应源项计算采用的是**层流有限速率反应（Laminar finite-rate chemistry）**模型，即假定化学反应速率仅由网格求积点上的瞬时大尺度热力学量（平动/振动温度、平均偏密度）直接决定。然而，在实际的高 Reynolds 数高超声速湍流（如超燃冲压发动机燃烧室内部）中，亚网格尺度的脉动（Sub-grid Scale Fluctuation）对化学反应通路具有极强的放大或抑制作用。缺失诸如概率密度函数（PDF）或增厚火焰面（TFM）等湍流-化学反应相互作用（TCI）模型，使得 MARUT 无法用于精确预测诸如燃尽极限、点火延迟以及超声速剧烈湍流燃烧等高复杂度热力学工程问题。

3. 粘性边界层网格处极其苛刻的时间步长制约（Explicit Transport Time-step Constraint）

虽然化学源项采用了隐式 Newton 求解，但空间对流与扩散项仍采用显式 SSPRK54 推进。在进行高 Reynolds 数粘性壁面流模拟（如 ONERA M6 翼型）时，为了解析极薄的粘性底层，壁面法向的 AMR 网格间距 $\Delta x_{min}$ 会被压缩得极小。由于粘性扩散项的稳定性限制条件为 $\Delta t_D \propto \Delta x^2 / (2P+1)^2$，这导致全局物理推进步长 $\Delta t$ 极其微小，需要耗费数百万个物理时步。MARUT 目前缺乏健壮的、支持完全 GPU 并行的隐式算子（如下三角-上三角对称松弛 LU-SGS 或预条件 Krylov 子空间方法 GMRES）来攻克粘性扩散时间尺度瓶颈。

4. 显存容量瓶颈对高阶大模型的制约（GPU Memory Footprint Limit）

高阶间断伽辽金法在提升多项式阶数 $P$ 时，由于其在三维空间下的张量积特性，单元局部的自由度数量呈 $(P+1)^3$ 爆炸式增长。MARUT 的每一步 RK 推进、局部雅可比矩阵构建以及 $L_2$ 投影算子均需要消耗大量的辅存空间。在单个 GPU 显存容量（如 L40S 仅 48GB）有限的物理现实下，这极大地限制了单节点所能承载的物理网格上限，迫使科研人员必须采用大量的并行节点，从而无形中抬高了计算的经济硬件成本。

5. 补充技术洞察：微观量子化学向宏观 CFD 模拟的跨尺度桥梁

对于计算化学和分子物理研究人员而言，MARUT 的真正潜力不仅仅在于其快速求解流体力学方程，更在于它充当了微观量子化学（Microscopic Quantum Chemistry）与宏观高超声速工程（Macroscopic Hypersonic Engineering）之间的跨尺度物理化学计算桥梁。

5.1 微观量子物理常数向宏观热力学多项式的映射

在 MARUT 使用的双温度化学反应源项中，分子的本征平动/转动感应热容、电子-振动能量态分布等宏观性质，本质上是大量微观量子统计状态的系综平均。例如，各化学组分在高温下的平动-转动比热、焓值多项式（NASA-9 拟合曲线）：

$$\frac{C_{p,s}^0(T)}{R} = a_1 T^{-2} + a_2 T^{-1} + a_3 + a_4 T + a_5 T^2 + a_6 T^3 + a_7 T^4$$

这些多项式的系数 $a_1 \sim a_7$ 并不是无源的水，它们直接来自于微观量子化学第一性原理计算（Ab Initio Calculations）。通过求解各分子的电子 Schrödinger 方程，获取电子能级、简谐振动频率 $\omega_e$ 以及转动惯量 $I_{rot}$，进而构建出其量子配分函数（Partition Function）$Q_{tot} = Q_{trans} Q_{rot} Q_{vib} Q_{elec}$。宏观内能和比热即通过配分函数的热力学微商导数直接输出：

$$U_s = R T^2 \left( \frac{\partial \ln Q_{tot}}{\partial T} \right)_{V}$$

MARUT 内部采用的 Park 双温度控制方程，实际上是将配分函数中的振动简谐振子能级单独剥离，使其在特定的慢速电子碰撞松弛时间 $\tau_s$ 下向平动能级靠拢。这体现了量子物理中能量在不同量子自由度之间流转、分配的动态历程。MARUT 的出现，让计算化学家能够在复杂多维流动（例如强震波反射、高速膨胀流喷管）中，直接测试、检验和校准微观碰撞截面（Collision Cross Section）与配分函数截断能级的准确性。

2.2 赋能可微编程与科学人工智能（Differentiable Programming & SciML）

MARUT 基于 Julia 语言 编写的本源特性，赋予了它一个传统基于 C++/Fortran 编写的 CFD 软件所无法企及的生态优势：天然兼容可微编程（Differentiable Programming）和科学机器学习（SciML）生态系统。

 ┌────────────────────────────────────────────────────────┐
 │                 MARUT 物理计算引擎                     │
 │   (高阶间断伽辽金物理输运 + 隐式有限速率化学反应求解)  │
 └──────────────────────────┬─────────────────────────────┘
                            ▲
                            │ (梯度反向传播: ChainRules.jl / Zygote.jl)
                            ▼
 ┌────────────────────────────────────────────────────────┐
 │                 Julia 科学机器学习生态                 │
 │  (SurrogateModel 神经网络 / 反应常数参数辨识)          │
 └────────────────────────────────────────────────────────┘

通过结合 Julia 的自动微分（Automatic Differentiation, AD）库（如 Zygote.jl 和 ChainRules.jl），MARUT 的整个物理仿真流程可以变为一个“可微物理算子”。这为以下两个前沿跨领域研究提供了无限的可能：

  1. 反应速率常数的逆向参数发现（Inverse Parameter Discovery）：
    在高温非平衡态空气动力学中，许多微观电离、非平衡辐射反应的 Arrhenius 常数（$A_r, \beta_r, E_r$）在极高温度下的实验测量值存在巨大的误差和不确定性。科研人员可以利用 MARUT 进行全流场模拟，将高超声速风洞中测得的震波脱体距离（Shock Stand-off Distance）或光谱发射线强度作为目标损失函数 $L$，通过 MARUT 物理计算引擎进行梯度自动反向传播（Backpropagation），在数百维的参数空间中直接计算梯度并逆向更新微观化学反应速率参数，极大地加速了非平衡态物理常数的精准辨识。
  2. 物理神经网络（PINNs）与非平衡热动力学替代模型：
    在 MARUT 内部，最耗时的部分依然是三维极端非平衡源项的隐式 Newton 积分。科研人员可以利用 GPU 上的历史高精度数值轨迹数据，在线训练一个多层深度算子网络（DeepONet）或物理感知神经网络（PINN），将其直接嵌入到 MARUT 的时间积分大循环中，用极其廉价的神经网络前向推理（Inference）直接替代昂贵的多维刚性 Newton 迭代，从而在保证空间高阶精度的同时，实现百倍以上的仿真加速。这种物理机制与深度学习在 GPU 设备端的深度融合，代表了下一代 E 级超算数值仿真不可逆转的技术趋势。