A Fully Matrix-Free Three-Grid Preconditioner for the Time-Harmonic Maxwell Equations at Extreme Scale

技术摘要：一种用于时谐麦克斯韦方程组的全矩阵自由三网格预条件算法

问题陈述
三维时谐麦克斯韦模拟会产生大规模、复杂且具有不定性的线性系统。为了在高频机制下保持相位精度，需要精细的空间分辨率，这导致未知元的数量随传播距离呈立方级增长。由于旋度-旋度算子（curl–curl operator）具有巨大的梯度核，且波项具有不定性，加之完美匹配层（PML）截断引入了复杂的各向异性系数，这些系统在求解上具有本质上的难度。

现有求解器面临着可扩展性的瓶颈。稀疏直接求解器在三维情形下存在超线性内存填充（fill-in）和通信开销问题。虽然领域分解法和扫掠预条件算法（sweeping preconditioners）可以提供较低的迭代次数，但它们的定向依赖性和因子分解存储需求与现代 GPU 的大规模并行吞吐量并不匹配。标准的几何多网格法虽然天然适用于 GPU，但最终会受限于精确粗网格因子分解的内存和通信成本，这在极端规模下是难以承受的。

方法论
作者提出了一种专门针对采用最优混合求积法（optimally blended quadrature）的旋度共形 Nédelec 有限元离散化的全矩阵自由、无因子分解三网格预条件算法。其核心创新在于将物理问题与稳定机制进行策略性解耦：

1. 三网格层级结构： 求解器在三个嵌套的六面体网格（$h$, $2h$, $4h$）上运行。

   • 细网格 ($h$)： 使用外部 FGMRES 求解器求解未移位的物理麦克斯韦方程（$A_h$）。
   • 中间网格 ($2h$)： 通过内部 FGMRES 过程计算不精确的修正。至关重要的是，这里的 Krylov 算子仍保持为未移位的物理算子（$A_{2h}$），但由一个移位的辅助循环进行预处理。
   • 粗网格 ($4h$)： 辅助循环的底层使用移位算子（$B_{4h}$）来抑制欠解析分量。

2. 移位辅助循环： 为了在不改变物理谱的情况下稳定不定波方程，复移位（$\sigma$）仅应用于辅助层级（$B_\ell = C_\ell - (1+i\sigma)k_0^2 M_\ell$）。该移位严格限制在中间网格修正的预条件器中。主要的物理算子（$A_h$）保持不变，从而为外部求解器保留了物理麦克斯韦算子。

3. 矩阵自由实现： 该方法避免了组装全局矩阵或进行精确的粗网格因子分解。

   • 算子： 所有层级均通过最优混合求积法（旨在将色散误差抑制至四阶）进行重新离散化，而非使用 Galerkin 乘积。
   • 传递： 边传递（上升与限制）通过基于坐标的插值进行几何传递，而不形成全局传递矩阵。
   • 平滑器： 各层级均使用固定工作量的 Krylov 平滑器（带有阻尼 Jacobi 的 GMRES），消除了对内部收敛检查的需求，并使内存占用严格与边未知元的数量成正比。

4. 参数选择： 采用局部傅里叶分析（LFA）来推导混合麦克斯韦分支和兼容的边传递。该分析确定了鲁棒的参数，包括复移位（$\sigma = 0.4$）和 Jacobi 阻尼因子，确保了移位辅助逆算子是有界的，且不会出现实频率奇异性，同时将预处理后的谱限制在有利区域内。

核心贡献

• 架构设计： 一种新型三网格预条件器，通过将复移位隔离到辅助内部循环中，允许外部求解器在未移位的物理算子上运行。这解决了稳定不定波方程与保持物理精度之间的矛盾。
• 极高可扩展性： 该算法通过消除因子分解存储和全局组装实现了极高的可扩展性。它完全依赖于几何边传递和固定工作量的 Krylov 平滑器。
• 优化的离散化： 使用最优混合求积法结合最低阶 Nédelec 元，在不增加系统规模或引入辅助未知元的情况下，提供了卓越的相位保真度。
• GPU 高效性： 实现方案针对 GPU 吞吐量进行了优化，利用基于单元的矩阵自由内核，与存储型模板（stored-stencil）方法相比，最大限度地减少了全局内存带宽的使用。

结果
该方法通过解析麦克斯韦格林张量（Green tensors）进行了验证，并在四种不同的介质（均匀介质、收敛透镜、周期介质和随机各向同性介质）上进行了测试。

• 规模： 求解器成功处理了包含约 108.9 亿个复数边未知元 的系统。
• 性能： 使用统一的求解器配置（无需针对特定介质进行调优），最大的系统在仅 64 块 NVIDIA A100 GPU 上运行的情况下，在 42.0 至 72.0 秒 内完成求解。
• 效率：
  • 弱扩展性： 随着问题规模和 GPU 数量增加（最高达 64 块 GPU），效率保持在 69.6% 至 91.0% 之间。
  • 强扩展性： 在固定问题规模下增加一倍 GPU 数量时，并行效率保持稳健（72.0%–93.8%）。
  • 内存： 每块 GPU 的峰值内存使用量始终低于 37.7 GiB。
  • 设置： 矩阵自由设置时间对于所有配置均低于 0.8 秒。
• 精度： 研究发现复单精度足以在保持数值结果与双精度几乎无偏差的同时，提供显著的性能提升。

意义
论文声称，这些计算代表了目前已报道的最大规模 3D 时谐麦克斯韦模拟，且采用了旋度共形有限元。其意义在于证明了通过无需粗网格因子分解或大规模内存占用的方式，可以实现极大规模的电磁模拟。通过将复移位限制在辅助预条件器中，并利用全矩阵自由、无因子分解的架构，该方法在现代 GPU 硬件上实现了前所未有的速度和可扩展性，解决了以往在该离散化保真度水平下难以处理的大规模问题。