基于GPU的无矩阵多重网格稳定广义有限元方法在界面问题中的应用

《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：GPU-based matrix-free multigrid methods of stable generalized finite element methods for interface problems

　　•基于GPU的无矩阵多重网格方法在处理界面问题的半光滑有限元法中的应用。•可通过舒尔补方案将标准有限元多重网格运算应用于半光滑有限元法。•一种用于半光滑有限元法中有限元部分、富集部分及交叉部分的混合无矩阵方案。•整体计算复杂度与传统有限元法相当。

引言
在大多数工程模拟和实际应用中，问题解往往具有非平滑特性，如裂纹、多物理场现象、接触问题以及复杂几何结构[1]。一方面，为非平滑特性生成适配的网格极为困难，尤其是当这些非平滑特性会随时间变化时，比如裂纹扩展或移动界面。另一方面，非平滑问题的解可能呈奇异状态或具有极大梯度，导致传统有限元法无法获得足够精确的近似结果。广义/扩展有限元法是解决非平滑问题的有效工具，在过去二十年里取得了显著进展。这类方法的核心思路是通过引入能够捕捉问题局部特性的附加函数，来提升对非平滑问题的近似精度[1]、[2]。这些局部函数（称为富集函数）是通过单位分解法“粘贴”到有限元网格中的[3]。通常，广义/扩展有限元法使用简单、固定且与非平滑特征无关的未适配网格，例如界面表面、域边界、裂纹面等，因此其网格生成工作相对简单。关于广义/扩展有限元法的更多内容可参考[1]、[2]、[4]。

与大量近似结果及工程应用相比，针对广义/扩展有限元法中方程组高效求解的研究相对较少。这一领域存在若干难点：一方面，有限元函数与富集函数之间可能存在近乎线性的依赖关系，从而导致得到的刚度矩阵条件数不良或呈奇异状态[5]。不良的条件数会给消元法带来严重的舍入误差，也会让迭代法需要更多迭代次数。另一方面，广义/扩展有限元法的刚度矩阵由有限元部分、富集部分及交叉部分构成，并非标准的有限元矩阵，而且粗细网格的试验空间也不具备嵌套结构。因此，要为广义/扩展有限元法开发出如同传统有限元法中的多重网格方法那样的快速求解器颇具挑战性[6]、[7]、[8]。本研究旨在探索适用于广义/扩展有限元法的基于多重网格的快速求解器。现有文献中已有一些针对广义/扩展有限元法的多重网格方法。在[9]中，通过基于伪L2投影的传递算子，从非嵌套网格的层次结构中构建出嵌套的有限元空间层次，从而实现了多重网格方法的应用。有研究针对三维裂纹扩展问题，提出了连接不同尺度的XFEM全局及局部多重网格方法[10]。在[11]中，通过修改延拓算子的稀疏模式以避免跨裂纹插值，提出了一种针对XFEM的准代数多重网格方法。针对断裂问题，[12]提出了XFEM的代数多重网格方法。在[13]中，开发了一种基于多重网格的重分析方法，用于高效预测疲劳裂纹的扩展情况。为了解决界面问题的高阶离散化问题，有研究提出了两级嵌套的几何多重网格方法[14]。在[15]中，设计了基于有限元部分与富集部分的块高斯-赛德尔迭代方案作为外部迭代，再以完整的多重网格方法作为内部迭代来处理有限元块。其他未适配网格有限元法的多重网格方法可参考[16]、[17]、[18]。

上述研究成果均基于CPU架构，而针对广义/扩展有限元法的GPU加速技术则很少被应用。一般来说，GPU加速要求程序和算法具备良好的并行性且内存占用较低。首先，这些方法需要通过组装全局刚度矩阵来实现，而存储矩阵及进行求解会占用大量GPU内存，尤其是在大规模计算中。其次，对于广义/扩展有限元法（以及其他未适配网格方法），构建平滑化算子与传递算子需要复杂的处理方式，这可能会降低算法的并行性，阻碍GPU加速技术的发展。第三，有些针对广义/扩展有限元法的多重网格方法是基于惩罚项或尼奇技术开发的，其中惩罚参数的选择取决于具体问题，且需要通过实验确定。当使用特定配置的NVIDIA AMGX硬件时，复杂的非标准有限元矩阵结构可能会导致收敛行为难以控制。最后，那些适用于传统有限元法GPU加速的无矩阵技术[19]、[20]、[21]，目前尚未在广义/扩展有限元法中得到充分应用。

本研究旨在提出基于无矩阵多重网格方案的广义/扩展有限元法GPU加速技术，重点研究界面问题。界面问题广泛应用于流体-结构相互作用、两相流、多孔介质流动以及多材料领域[22]、[23]。为界面问题生成适配的网格极为困难，尤其是对于三维或随时间变化的界面问题。近几十年来，未适配网格有限元法在界面问题研究方面取得了显著进展，包括基于惩罚项或拉格朗日乘数的有限元法[24]、[25]、[26]，割集有限元法[27]、[28]，浸没有限元法[29]、[30]，以及广义/扩展有限元法[15]、[31]、[32]、[33]、[34]、[35]、[36]、[37]。广义/扩展有限元法的优势在于：其一，其形函数无需考虑本构关系和材料参数即可构建；其二，该方法可以采用符合变分原理的形式，无需引入惩罚项和参数。在本文中，我们采用稳定的半光滑有限元法来解决界面问题。半光滑有限元法是广义/扩展有限元法的一种稳定版本，能够改善广义/扩展有限元法的条件数。半光滑有限元法最早由[5]、[38]、[39]提出，其核心思路是通过减去有限元插值项来修改广义/扩展有限元法中的富集函数。已有研究探讨了半光滑有限元法在界面问题、裂纹问题等多种场景中的应用[15]、[36]、[40]、[41]、[42]、[43]、[44]、[45]、[46]。研究表明，半光滑有限元法在解决界面问题时能够以最优速率收敛[36]、[37]、[40]。此外，半光滑有限元法的条件数与普通有限元法处于同一数量级，而且具有良好的鲁棒性，即当界面接近元素边界时，其收敛性和条件数都不会出现恶化。

无矩阵技术是降低GPU内存占用并开发快速多重网格求解器的关键。无矩阵方法通过逐元素汇总局部向量-向量乘法的结果来计算矩阵-向量乘法，从而减少了全局刚度矩阵的存储量。更重要的是，现代GPU架构因其高度的并行性，能够很好地支持无矩阵运算，因此迭代方法从无矩阵方案中受益匪浅。尽管无矩阵技术在基于适配网格的有限元法和离散化伽辽金法中取得了巨大成功[19]、[20]、[21]、[47]、[48]、[49]，但由于上文提到的缺乏嵌套子空间和富集框架的问题，它在半光滑有限元法（以及其他广义/扩展有限元法和未适配网格有限元法）中应用面临诸多困难。本研究通过充分利用半光滑有限元法的结构解决了这一问题。首先，我们利用舒尔补方案将半光滑有限元法的线性系统转化为两个规模较小的系统，对应的矩阵分别为A?11=A11+P和A22。其中A11是基于未适配网格的标准有限元矩阵，P为扰动项；A22则是比A?11少一维的小型稀疏矩阵。因此，大部分计算工作集中在求解包含A?11的系统中。其次，通过分析扰动项P的影响，我们发现A?11的条件数与A11处于同一数量级，这意味着扰动项P的影响很小。这一发现促使我们决定采用与处理A11相同的操作来构建A?11的多重网格结构，因此无需再设计新的延拓算子和限制算子。第三，我们可以为这种多重网格方法开发GPU无矩阵加速方案。新的无矩阵方案由标准的有限元无矩阵运算以及对扰动项P的特殊处理组成。所提出的MFMG方法保留了半光滑有限元法的所有优点，且实现速度较快。此外，它的计算复杂度并未显著增加。对于半光滑有限元法而言，MFMG方法的计算成本与普通有限元法处于同一数量级。我们还通过与传统的GPU求解器进行了比较，证明了所提出的MFMG方法在计算时间和内存占用方面具有优势，包括在NVIDIA AMGX硬件上的表现。

本文的其余结构如下：第2节介绍模型问题以及传统的广义/扩展有限元法和半光滑有限元法；第3节阐述传统有限元法的多重网格与无矩阵方案，同时介绍现有的半光滑有限元法多重网格方法；第4节提出针对半光滑有限元法的MFMG方法，并分析该方法的计算复杂度；第5节进行数值实验，并与NVIDIA AMGX硬件进行对比；第6节给出结论与未来研究方向。

模型问题、传统有限元法、半光滑有限元法
对于R2空间中的区域Λ，整数m以及1≤q≤+∞，我们用Wm,q(Λ)表示常用的索伯列夫空间，其范数为‖?‖Wm,q(Λ)，半范数为|?|Wm,q(Λ)。当q=2时，Wm,q(Λ)对应于Hm(Λ)空间；当m=0时，则对应于Lq(Λ)空间。我们考虑一个有界且单连通的开区域Ω?R2，其边界?Ω为分段光滑的曲线。该区域被界面Γ划分为两个子区域（开集），即Ω0和Ω1，满足Ωˉ=Ωˉ0∪Ωˉ1，Ω0∩Ω1=0?，且Γ=Ωˉ0∩Ωˉ1，其中D?表示开集D的闭包。

有限元法的多重网格与无矩阵方案以及现有的半光滑有限元法多重网格方案
为更好地说明半光滑有限元法MFMG方法的思路，我们先简要介绍传统有限元法的多重网格与无矩阵方案。多重网格方法[6]、[7]、[8]、[54]是解决由有限元离散化产生的大型稀疏线性方程组的最有效求解工具之一。这种高效性源于离散化误差的尺度分离特性：简单的松弛算法能够快速抑制高频误差分量，而剩余的平滑低频误差则只能以较慢的速度被消除。

针对界面问题的半光滑有限元法MFMG方法
我们同样受到式(2.15)中矩阵A的块结构启发，但不会像式(3.8)那样引入外部迭代。为此，我们将线性方程式(2.13)用舒尔补技术改写为A?11u1=f?1，A22u2=f2?A12Tu1，其中A?11=A11?A12A22?1A12T，f?1=f1?A12A22?1f2。我们的算法采用直接求解的方式：首先求解方程式(4.1)得到u1，再将u1代入方程式(4.2)求解u2，这样就避免了外部迭代。方程式(4.2)可以通过简单的迭代法或消元法求解。

数值实验
在接下来的实验中，我们首先介绍所使用的硬件和软件环境。实验平台包括一台NVIDIA GeForce RTX 4090 GPU和一台Intel(R) Xeon(R) Platinum 8360Y CPU，操作系统为Ubuntu 22.04。有限元运算是在deal.II FEM库中实现的，该库具备自由度管理、线性代数运算功能，以及构建有限元多重网格层次结构的能力。按照第4节中描述的异构计算策略，批量运算在GPU上以无矩阵方式执行，而稀疏运算则……

结论与评论
在本文中，我们提出了一种面向GPU的无矩阵多重网格方法，用于在未适配笛卡尔网格上的界面问题半光滑有限元法计算。该方法的核心思路是通过舒尔补方案重新构造半光滑有限元法的线性方程组，从而将主要的计算任务集中在与标准有限元自由度相关的子系统上。我们设计了一种部分组装式的多重网格策略，将该策略结合了体区域内的无矩阵运算与界面附近的局部代数处理方法。所提出的……

CRediT作者贡献说明
郭志厚：撰写——初稿、验证、软件实现、方法设计。崔楚：验证、软件实现、方法设计、概念构思。张庆辉：撰写——审阅与编辑、撰写——初稿、验证、监督、资源协调、项目管理、方法设计、研究实施、资金筹集、概念构思。

利益冲突声明
作者声明不存在任何可能影响本文研究结果的已知财务利益或个人关系。

郭志厚|崔楚|张庆辉