一种用于弹性动力学时域边界元素方法的并行算法，采用mCCSR稀疏存储技术

《Advances in Engineering Software》：A parallel algorithm for the time-domain boundary element method in elastodynamics with mCCSR sparse storage

　　周维宇|陈永强
清华大学核能与新能源技术研究所，先进核能技术协同创新中心，教育部先进反应堆工程与安全重点实验室，北京100084，中国

**摘要**
本研究提出了一种基于高性能CPU的并行策略，用于弹性动力学时域边界元方法（TDBEM），该策略采用了混合单元压缩稀疏行（mCCSR）格式。开发了一种定量任务负载模型，以考虑不同单元的积分成本，从而实现动态的超行划分方案，在矩阵生成过程中有效平衡了各个线程的工作负载。在求解阶段，引入了一种统一的、多步的并行方案用于卷积计算，而GMRES迭代求解则采用了基于mCCSR的稀疏矩阵-向量乘法。所提出的框架充分利用了mCCSR格式的单元级稀疏性，最小化了线程调度和同步的开销。通过涉及矩形盒、球体、T梁和复杂机械组件的数值示例，证明了并行效率的显著提升，对于大规模问题，使用48个线程时效率可达到64%。动态划分方法在处理非对称几何体时始终优于传统的平均划分方法，同时保持了稳定的和物理上准确的结果。该方法为大规模弹性动力学模拟提供了一种可扩展且可靠的解决方案，并具有扩展到混合MPI-GPU实现的潜力。

**引言**
工程问题的几何复杂性和计算规模的不断增加，使得传统的串行程序无法满足大规模瞬态弹性动力学分析的需求。硬件架构的进步为并行计算提供了有力支持，使得高效并行算法的开发成为解决大规模结构动力学问题的关键。然而，由于负载不平衡、数据依赖性、并行工作负载组织以及内存和缓存使用的优化等问题，实现高性能并行计算仍然具有挑战性。
在各种数值方法中，边界元方法（BEM）[1,2]因其降维能力、高精度以及对无限和半无限域的自然适应性而受到重视，成为有限元方法（FEM）的宝贵补充。然而，传统BEM生成的系数矩阵通常是密集的、非对称的且满秩的，这给大规模问题带来了巨大的存储和计算负担。为了解决这些问题，已经开发了低秩和分层加速技术，如自适应交叉近似（ACA）[[3], [4], [5], [6]]和快速多极法（FMM）[7,8]。尽管取得了这些进展，但这些方法在单核架构上的性能仍不足以满足复杂工程系统的实时仿真需求。最近，利用密集矩阵子块的低秩特性的并行实现和负载均衡策略已经出现[[9], [10], [11], [12], [13]]，显著提升了计算效率。Wang等人[14]提出了基于CUDA的GPU加速快速多极BEM，用于三维弹性分析。Galvis等人[15]开发了BESLE，这是一种用于三维线性弹性的开源并行边界元软件，支持弹性静态和动态分析，以及各向同性和各向异性材料以及单域或多域分析；Peres等人[16]则提出了一种基于CUDA的弹性固体并行等几何边界元分析框架，该框架基于T样条。
时域边界元方法（TDBEM）已广泛应用于瞬态弹性动力学问题[11],[17],[18],[19]]。虽然标准的显式TDBEM公式受严格的Courant–Friedrichs–Lewy（CFL）条件控制，可能会导致后期数值不稳定，但最近的发展成功引入了适用于三维弹性力学问题的无条件稳定公式[[20],[21],[22]]。这些先进方案允许在长模拟时间内获得稳定结果，而无需执行CFL类型的限制，标志着BEM领域的重要进展。Li等人[5]提出了一种基于GPU的并行方案，包含了矩阵行平均划分策略。Bramas等人[23]通过矩阵切片和重新排序改进了飞机声场的仿真。Saitoh等人[24]利用混合MPI-OpenMP并行方法研究了无限域中的三维多腔粘弹性波传播。
尽管像FM-BEM这样的快速算法显著加速了静态问题，但由于复杂瞬态基本解和多波速的存在，将其扩展到三维弹性动力学TDBEM在数学上较为复杂。另一方面，利用弹性动力学TDBEM的独特物理特性提供了一种更为实用的方法。其时变基本解受波传播因果关系的控制，导致系数矩阵呈稀疏的、块状结构，并具有部分子矩阵对称性。基于这些特性，周维宇和陈永强提出了混合单元压缩稀疏行（mCCSR）格式。该方案避免了复杂的代数近似，在数据结构和硬件层面实现了显著加速，同时严格保持了数值精度。然而，尽管mCCSR在卷积计算和GMRES迭代中显著提升了计算可扩展性，但这些系数矩阵的生成仍然是主要的计算瓶颈，限制了整体效率。
在mCCSR格式中，超行的稀疏性变化很大，使得传统的BEM并行策略[5,23,24]难以直接应用。相比之下，有限元领域中的稀疏矩阵并行技术相对成熟[[26],[27],[28]]。例如，David等人[28]采用了一种动态划分策略，利用前一步的时间步长信息有效平衡了计算负载。高性能计算（HPC）的最新进展成功利用八叉树网格在现代CPU和GPU架构上加速了大规模模拟。例如，Mohammadian等人[29]开发了一个专门用于弹性动力学模拟的多GPU框架。在缩放边界有限元方法（SBFEM）的背景下，Zhan等人[30]提出了一种基于八叉树的网格生成技术，用于大规模并行计算；Liu等人[31]引入了一种高效的隐式-显式时间步进策略，用于三维动态相互作用分析。受这些成功基于域的并行框架的启发，本研究探索了针对TDBEM的先进并行调度策略。
本文研究了使用mCCSR格式的弹性动力学TDBEM，并提出了一种高效的共享内存多核并行方案，针对两个关键阶段：矩阵生成和求解。对于矩阵生成，开发了一种基于三种时域单元的工作负载量化模型；动态负载划分策略利用前一个时间步长的信息，确保了跨时间步长的计算平衡。对于求解阶段，统一调度方案将多时间步长的稀疏矩阵-向量乘法（SpMV）和向量加法整合到单次 passes 中，大幅减少了线程调度和同步的开销。GMRES迭代求解通过SpMV和针对mCCSR结构定制的向量点积操作实现了并行化，从而显著提高了整体速度。此外，为克服数值不稳定性，本文的求解器中内嵌了时间加权稳定策略[32,33]。本文的重点是并行算法设计和性能，而mCCSR格式的详细存储特性已在我们的先前工作中报告。因此，本研究仅从并行工作负载和调度的角度简要讨论了计算方面，主要关注新的并行策略。

**需要强调的是，本研究的方法论贡献与现有的并行框架不同。** 与为拓扑上静态稀疏矩阵设计的通用HPC或FEM策略不同，弹性动力学TDBEM中的稀疏性随时间传播的波前的变化而动态演变。与传统多域BEM不同，后者主要是通过几何分解和子域耦合引入块结构，而此处考虑的稀疏块模式则是弹性动力学TDBEM中波传播因果关系的自然结果。此外，尽管之前的mCCSR研究[25]主要集中在静态存储格式上，但本研究通过提出一种基于物理的动态调度算法进行了改进。通过基于不同时域单元状态的数学量化工作负载，所提出的方法将物理波因果关系转化为严谨的动态负载平衡模型。相比之下，FM-BEM主要通过快速近似来减少密集BEM的复杂性，而本研究则通过利用因果诱导的稀疏性和基于mCCSR的数据结构来提高弹性动力学TDBEM的并行效率。

**本文结构如下。** 第2节介绍了弹性动力学TDBEM的基本原理。第3节简要概述了mCCSR存储格式。第4节分析了矩阵生成和求解阶段的并行化策略。第5节详细描述了提出的并行实现方案。第6节通过数值示例（包括矩形棱柱、球体、T梁和复杂工程组件）评估了算法的性能。最后，第7节总结了研究内容，并讨论了向MPI/GPU混合并行化的潜在扩展。

**章节片段**
- **弹性动力学时域边界元方法**
本节介绍了在无体力作用和初始位移及速度假设下的弹性动力学时域边界积分公式，随后介绍了其离散化和基于卷积的数值实现。
在没有体力作用的情况下，弹性动力学问题的控制方程为纳维-柯西方程：
$$
\mu ui,jj(x,t) + (\lambda + \mu) u_j,ji(x,t) = \rho u''i(x,t)
$$
其中μ和λ是拉梅参数，ρ是质量密度，ui表示位移。

- **mCCSR方案**
在弹性动力学时域边界元方法（TDBEM）中，基本解包含具有时间依赖性的冲击函数。由于波传播遵循因果关系，源点产生的响应在每个时间步长都会延迟，导致矩阵具有稀疏特性。此外，根据边界积分方程的结构，在源点和场点之间形成了3 × 3的子矩阵，其中G矩阵的子矩阵具有对称性。

- **mCCSR格式在弹性动力学时域边界元方法中的并行性分析**
本节全面分析了使用mCCSR存储格式的弹性动力学时域边界元方法（TDBEM）在系数矩阵生成和求解阶段的并行特性。通过对矩阵生成阶段的单元积分工作负载、稀疏分布模式和步间依赖性的定量研究，识别了并行效率低下和计算不平衡的潜在来源。

- **采用mCCSR格式的弹性动力学时域边界元方法的并行实现方案**
本节提出了一种高性能的并行实现策略，用于弹性动力学时域边界元方法（TDBEM）。该方案涵盖了所有关键计算阶段，包括系数矩阵生成、卷积计算和GMRES迭代求解，并概述了一个全面的并行工作流程。该策略旨在实现计算负载平衡，最小化线程调度和同步开销。

- **数值示例和性能分析**
本节通过几个数值示例验证了所提出的基于mCCSR的弹性动力学时域边界元方法（TDBEM）的并行方案的效率。所有计算均在配备Intel(R) Xeon(R) Platinum 8488C CPU（48核）的服务器上进行。所有数值实验均在单个CPU上完成，禁用了超线程，并使用默认操作系统调度器管理并行任务。

- **结论**
本研究开发了一种基于高性能CPU的并行计算框架，用于弹性动力学时域边界元方法（TDBEM），该框架采用混合单元压缩稀疏行（mCCSR）格式。主要贡献在于建立了一种统一且可扩展的并行化策略，将定量任务建模、动态负载平衡和优化的求解器级调度集成到了mCCSR稀疏块结构中。

**作者贡献声明**
周维宇：撰写——原始草案、可视化、验证、软件开发、资源管理、方法论研究、资金获取、形式化分析、数据整理、概念化。
陈永强：撰写——审阅与编辑、监督、资金获取。