用于具有不连续系数的椭圆问题不连续Galerkin离散化的p-多网格方法

《Computer Methods in Applied Mechanics and Engineering》：p-multigrid method for the discontinuous Galerkin discretization of elliptic problems with discontinuous coefficients

　　诺磊|张东航|郑伟英 中国科学院数学与系统科学研究院华罗庚数学科学中心，北京，100190，中国

摘要
在本文中，我们扩展了之前关于W循环p多网格方法的研究，利用对称加权内点惩罚离散Galerkin方法进行高阶离散化，从而高效求解具有不连续系数的椭圆问题。我们在解具有分段H2正则性的假设下，针对存在不连续系数和分层Legendre多项式基的情况，给出了p多网格方法的严格收敛性分析。我们的分析表明，收敛速率与网格尺寸、多项式阶数p以及不连续系数的比值无关。理论结果表明，该方法在处理系数有较大跳变的椭圆问题时既高效又稳健。值得注意的是，该方法仅需O(p)次平滑步骤即可达到良好性能。我们通过大量数值实验验证了理论结果，证明了所提方法的有效性与稳健性。

引言
椭圆界面或传输问题在科学与工程应用中十分常见，例如流固耦合、热传导以及复合材料建模等领域。这类问题通常涉及在具有不同材料特性的区域分界面上出现显著不连续性的扩散系数。本文考虑一个矩形或立方域Ω?Rd，重点研究以下二阶椭圆模型问题：
?div(κ?u)=finΩ,u=0on?Ω，
其中f∈L2(Ω)为给定的源项，κ∈L∞(Ω)为扩散系数，其值可能出现较大跳变。这类不连续性给数值方法带来了巨大挑战，使用有限元法或离散Galerkin法离散时往往会导致病态的刚度矩阵。矩阵的条件数会因系数跳变的幅度而受到严重影响，进而降低迭代求解器的性能与稳健性。因此，开发能够应对系数不连续性的预处理器对于实现可扩展且高效的求解器至关重要。

在为问题(1.1)的符合型有限元离散化设计高效且稳健的预处理器方面已经取得了显著进展。诸如Balancing Neumann-Neumann[1]和Bramble-Pasciak-Schatz[2]之类的非重叠域分解方法已展现出对系数跳变和网格尺寸（直至对数级别）的稳健性。然而，这些方法往往依赖于复杂的粗网格求解器，增加了实现难度。相比之下，如Bramble-Pasciak-Xu(BPX)预处理器[3]这样的多层方法则无需使用如此复杂的粗网格求解器。结合共轭梯度法，BPX能为具有不连续系数的问题提供稳健且高效的解决方案。尽管BPX在许多标准问题上表现最优，但对于系数高度不连续的椭圆方程，其性能可能取决于粗网格是否足以处理这些系数跳变[4]。此外，代数多网格法和几何多网格法[5]、[6]、[7]也被开发出来以改善条件数，每种方法都有各自的优缺点。这些发展凸显出根据问题特点定制预处理策略的必要性。

离散Galerkin方法特别适合处理具有不连续系数的问题，为求解(1.1)提供了灵活性与稳健性[8]。通过使用离散多项式空间，离散Galerkin方法可实现高阶精度，非常适合捕捉异质介质中的细尺度解特征。但这一优势也带来了计算复杂度的增加，因为自由度更高且矩阵结构更为密集。不连续系数进一步加剧了这一问题，提高了条件数，给高效迭代求解器带来了挑战。

已有不少有效的预处理器被开发用于离散Galerkin方法的离散化。值得注意的是，Dryja等人[9]、[10]提出了非重叠BDDC预处理器和Neumann-Neumann预处理器，并对其进行了分析，即使在不符合型的子域上，它们也能展现出与网格尺寸和系数跳变无关的准最优条件数。针对具有较大系数对比度的二阶IPDG和H(curl, Ω)-椭圆问题，还有多层预处理器被引入[11]、[12]，进一步推动了该领域的发展。

然而，这些方法大多是为低阶（≤2）的离散Galerkin离散化设计的。对于任意高阶的离散Galerkin方法，虽然也取得了一些进展，但可行的稳健预处理器仍然较少。非重叠Schwarz预处理器[13]和均匀加性Schwarz预处理器[14]已被探索过。在[15]中，一种可扩展且准最优的BDDC预处理器被提出，它对系数跳变具有稳健性，同时还能适应局部网格细化及多项式阶数的变化。尽管有了这些进展，但在存在较大系数不连续性的情况下为高阶离散Galerkin方法开发预处理器仍是一个充满挑战的未解问题。

p-MG方法和hp-MG方法为高阶离散Galerkin方案展现了巨大潜力[16]、[17]、[18]，它们通过利用低阶空间来加速高阶系统的收敛。尽管这些方法在实践中得到广泛应用，但其在离散Galerkin框架下的理论基础最早是在[19]、[20]中建立的。我们之前的工作[21]通过引入分层正交基和多项式平滑器进一步提升了这一框架，从而实现了更好的收敛性。不过，这些改进主要针对的是系数恒定的问题。对于存在不连续系数时的p-MG方法性能的理论理解仍不完善。

在本文中，我们将p-MG技术扩展到具有不连续系数的椭圆界面问题，并在这一更具挑战性的环境中给出了严格的收敛性分析。在我们之前的工作[21]基础上，我们采用了对称加权内点惩罚离散Galerkin方法[8]、[22]，该方法利用与扩散相关的权重来定义一致且对称的双线性形式。针对界面和边界跳变的惩罚项通过扩散系数的调和平均值进行缩放，从而确保在系数出现大幅跳变时仍能保持稳定性与精确性。我们的结果证实，所提出的p-MG方法在系数非恒定情况下依然具备稳健性与高效性。

为进一步展示所提p-MG方法的通用性，我们还在非贴合网格上的椭圆界面问题中对其性能进行了数值研究。这类问题在涉及复杂几何形状的应用中十分常见，因为生成贴合网格往往较为困难。在过去二十年中，离散Galerkin框架下的非贴合有限元方法受到了广泛关注[23]。通过将SWIPDG离散化方法与满足大元素条件的割单元合并算法[24]相结合，我们发现即使在存在较小割单元的情况下，该方法也能保持稳定且精确的数值性能。

需要指出的是，本文中的理论分析是基于在具有张量积多项式空间及分段常数标量扩散系数的均匀笛卡尔网格上进行的SWIPDG离散化。将分析范围扩展到更一般的网格或系数则超出了本文的研究范畴。

本文的主要贡献如下：
1. 我们采用SWIPDG方法，利用高阶精度求解具有较大系数跳变的椭圆问题。
2. 我们建立了分析p-MG方法稳健性与高效性的理论框架。我们的分析表明，收敛速率与网格尺寸h、多项式阶数p以及系数κ均无关。
3. 我们开展了大量的二维和三维数值实验来验证理论结果。特别是，我们证明所需的平滑步骤数量可减少到O(p)，从而提升了计算效率。
4. 我们还通过数值实验证明了所提方法在非贴合网格上、具有C2光滑界面的不连续系数椭圆界面问题中的适用性，其性能在非贴合网格环境下依然稳定且精确。

本文的其余结构如下：第2节介绍带有非继承双线性形式的SWIPDG方案，并给出关键技术结果；第3节基于分层断裂多项式空间定义延拓算子与限制算子，介绍多项式平滑器，并阐述W循环p-MG方法；第4节和第5节建立了对收敛性分析至关重要的逼近性质与平滑性质；第6节给出了两级p-MG方法和W循环p-MG方法的收敛理论；第7节将分析扩展到继承双线性形式的情况；最后，第8节给出大量数值实验结果，第9节总结全文并讨论未来工作方向。

章节节选
SWIPDG离散化与技术结果
在本节中，我们首先在2.1节介绍基本假设并定义计算网格，然后在2.2节构建分层基函数。接着在2.3节介绍SWIPDG方案并阐述关键技术结果。最后，2.4节给出向量化表达式及相应的刚度矩阵。

具有非继承双线性形式的p-MG方法
在本节中，我们基于非继承双线性形式（2.3）介绍p-MG方法的W循环机制。该算法建立在延拓算子、限制算子以及多项式平滑器的设计之上。具体而言，3.1节定义网格间传递算子，3.2节介绍多项式平滑器。3.3节给出了p-MG方法的两级算法与W循环算法的误差传播矩阵及主要收敛结果。

逼近性质
本节探讨p-MG方法的逼近性质。对于任意k=2,?,p，断裂多项式空间Uk(Th)的质量矩阵定义为vk?Mkwk=(vk,wk)Ω，其中vk,wk属于Uk(Th)。由于Legendre多项式的L2正交性，质量矩阵Mk为对角矩阵，即Mk=diag{∥ΨTσ∥L2(T)2:∥σ∥∞≤k,T∈Th}，属于RNhNk×NhNk。设Zk为对角矩阵，其元素为扩散系数的取值，即Zk:=diag{κ|T:∥σ∥∞≤k,T∈Th}，同样属于RNhNk×NhNk。受[Lemma 2.9, Theorem 2.9 [30]]的启发，我们得到了相关结论。

平滑性质
本节首先介绍刚度矩阵的稀疏性，以便于后续分析（5.1节）。接着（5.2节）推导平滑算子的谱半径估计。最后在5.3节确立多项式平滑器的平滑性质。

p-MG方法的收敛性分析
在本节中，我们给出主要收敛结果的证明。6.1节聚焦于两级算法，6.2节讨论W循环方法。

具有继承双线性形式的p-MG方法
为简化并优化矩阵组装过程的同时节省存储空间，对于p≥2且k=2,?,p，我们按如下方式定义继承双线性形式：
A?k?1(vk,wk):=A?k(Pk?1kvk?1,Pk?1kwk?1)，其中vk?1,wk?1属于Uk?1(Th)。p阶双线性形式定义为A?p(vp,wp):=Ap(vp,wp)，适用于任意vp,wp属于Up(Th)。继承型SWIPDG方案的目标是找到u?k属于Uk(Th)，使得A?k(u?k,vk)=(f,vk)Ω，对所有vk属于Uk(Th)均成立。刚度矩阵A?k是通过继承双线性形式（7.1）生成的，即vk?A?kwk=A?k(vk,wk)，对所有vk,wk成立。

数值实验
在本节中，我们通过数值实验验证理论结果，并评估所提出的W循环p-MG算法的收敛行为。算法从k=p这一层级开始实现，采用继承双线性形式{A?k}k=1p。在二维测试中，计算在配备2.9 GHz Intel Xeon Gold 6326 CPU（16核）的工作站上进行，采用未并行化的MATLAB实现。在三维测试中，我们使用了……

结论
在我们之前的工作[21]基础上，我们提出了一种高效的多网格方法，用于对具有系数不连续性的椭圆问题进行高阶SWIPDG离散化，该方法采用了四阶Chebyshev多项式平滑器。我们的理论与数值结果表明，该p多网格方法的收敛速率与网格尺寸h无关，且与多项式阶数p也无关。此外，该方法对具有异质扩散特性的问题也表现出良好的稳健性。

附录
按照[8]中的框架，我们在此详细证明了关于SWIPDG双线性形式的连续性与强制性的Lemma 2.5。我们同时分析了相关常数与网格尺寸h、多项式阶数k以及扩散系数κ之间的依赖关系。我们首先估算将在强制性和连续性论证中反复使用的、与一致性相关的面项。

引理10.1
对于任意(v,wk)属于Uk,*(Th)×Uk(Th)，有|∑F属于Fh∫F{{κ?v}}ω·nF[[wk]]|≤(∑T属于Th∑F属于?Thk2∥κ1CRediT作者贡献说明
诺磊：撰写——初稿、验证、软件、方法论、形式分析。张东航：撰写——审阅与编辑、验证、软件、方法论、形式分析。郑伟英：撰写——审阅与编辑、监督、方法论。

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