📝 求解不可压缩流体流动的基于有限元的神经网络及相关反问题
摘要
关键词：基于有限元的神经网络、机器学习、预处理、有限元、参数化PDE、Stokes方程、Navier-Stokes方程
引言
理解物理系统在科学和工程的许多应用中至关重要。许多重要的物理系统由偏微分方程（PDEs）建模，然后使用数值方法求解。近年来，一种用于模拟物理系统的较新方法是使用神经网络（NNs）作为替代模型，这种方法与数值方法相比，可以降低复杂度并加快评估速度。这些深度神经网络是通用函数逼近器 [1]，通常通过来自模拟或测量的数据进行训练。然而，生成实验或模拟数据的成本非常高。此外，基于数据驱动的方法不仅不足以充分考虑能量、质量和动量等自然法则，还难以推广到未见的情况。为了解决这些问题，物理信息机器学习（physics-informed machine learning）领域应运而生，结合了（数值）数学、机器学习/人工智能和计算工程学等学科，详见综述 [2, 3, 4, 5]。
在利用神经网络求解微分方程的众多方法中，尤其是所谓的物理信息神经网络（PINNs）[6]近年来得到了广泛关注。PINNs将相应微分方程的残差嵌入到神经网络的损失函数中，并通过自动微分精确地计算所需的导数 [7]。因此，这种方法通常被认为能够避免离散化；特别是，由于物理知识被编码在损失函数中，PINN方法能够解决微分方程而无需进行网格剖分。但是，由于损失函数中更高阶的导数，PINN方法在训练过程中存在高学习复杂度以及附加损失项以结合边界条件，这些都会导致多目标优化问题（至少在标准PINN方法中），因此训练PINNs是具有挑战性的，详见 [8]。
在后续文献调研中，我们专注于基于PINNs的扩展版本或与PINNs有相似之处的方法，如深度Ritz方法、Deep Nitzsche方法和弱对抗神经网络等——所有这些方法的共同点是它们并未显式应用离散化——我们将讨论这些方法。
接下来，我们聚焦于在一定程度上明确使用离散化的物理信息神经网络的文献综述。所谓的变分PINNs（vPINNs），一种类似于Petrov-Galerkin的方法，其中神经网络作为试函数并与有限数量的测试函数进行测试，引入于 [9]。在VarNet方法中，特定的有限元函数作为测试函数，而在-vPINNs [11]中，本质上使用-有限元函数作为测试函数。相比之下，Galerkin神经网络方法 [12] 扩展了经典的Galerkin过程，通常基于从线性空间中选取的试验和测试函数，改为从神经网络参数化的非线性函数集中选取试验和测试函数。PhyGeoNet [13] 将底层问题映射到一个矩形参考域，然后在这个矩形域的网格点上评估PINNs损失，同时结合卷积神经网络。扩展这一思路到较少结构化的网格中，物理信息图网络Galerkin网络 [14] 结合图卷积神经网络与有限元方法（FEM）离散化。在 [15] 中，考虑了有限元离散化，并且对基础网格的细化步骤进行了研究；此外，应用了块Jacobi预处理器。近年来，PINNs还被用于加速某些元素的局部计算 [16]，并且神经网络被用来修正更粗的粗网格的FEM解 [17, 18]。
在本综述的末尾，我们还简要评论了两个与此类似但仅略相关的方法。在 [19] 中，差分方程被使用经典FEM离散化，而神经网络用于参数化差分方程中的未知参数。在 [20] 中，有限元方法被嵌入到神经网络中。
最后，[21, 22] 直接结合了神经网络和有限元：参数化微分方程的离散化被优先由有限元进行，而神经网络用于训练参数（微分方程的解）与相应的离散解系数之间的关系（相对于选定的有限元基）。损失函数由FEM离散化微分方程的残差构成，这使得方法结合了数据驱动和物理信息。在 [21] 中，讨论了一维对流扩散问题的适用性，以及在某些场景中的应用，如压缩问题或Reissner-Mindlin平板问题。基于相同的优势，[21, 22] 建议针对边界或不规则形状的域进行逼近，相较于PINNs和数值参数化，扩展至全局剖分，但更适合后续学习的场景。
综上所述，结合了有限元的神经网络方法，为处理物理系统建模提供了一个新颖的工具。每种方法都有其优缺点，取决于具体问题的需求和约束条件。
首先，我们考虑静止的二维Stokes方程和Navier-Stokes方程，即类型和特征为鞍点结构的方程，对于经典的有限元方法（FEM）来说，这些方程比之前用此方法处理的方程更具挑战性。
其次，我们通过添加预处理器修改了物理信息损失函数。结合L-BFGS，这显著加快了训练过程，并最终提高了所获得解的准确性。
关于现有的物理信息机器学习在流体动力学问题中的应用文献，我们参考了前述综述文章中的非离散化方法，并仅评论了那些明确依赖离散化的技术。已经提到的神经网络多重网格求解器 [17] 已经应用于Navier-Stokes方程。DiscretizationNet [24] 结合了静态Navier-Stokes方程的有限体积离散化以及编码器-解码器卷积神经网络的迭代应用，而 [25] 使用了物理信息损失在网格上进行评估，并结合了神经网络的三维U-Net架构。
最后，我们提到除了这些基于机器学习的方法外，在过去几十年中开发了大量经典的模型降阶技术，特别是在流体动力学应用中；详见 [26]。
论文结构
在第2节中，我们介绍了所提出的方法；首先，我们描述了基础设定和一般思想，然后在第2.3和2.4节中就我们的两个具体模型问题对此方法进行详细讨论。第3节给出了说明我们方法优点的数值示例，第4节提出了对概率逆问题的原型应用。我们在第5节用简短的总结和对未来工作的展望结束论文。最后，附录A包含了一个可以理解为所选方法启发性动机的事实证明。
2. 物理信息、基于FEM的神经网络的预处理损失函数
在本节中，我们介绍我们提出的方法并简要讨论其实现，无论是在一般的高层次上，还是在具体细节上。之后，我们提供了一些更具体的细节，关于其实现和在两个具有代表性的、非平凡模型问题上的应用。
2.1 问题陈述
在以下内容中，我们将考虑可以写成以下抽象变分形式的问题：
其中表示适当的向量空间，在第二个参数上是连续且线性的，但不一定是双线性。事实上，我们甚至会考虑以下参数化版本：
其中表示来自参数区间的一组参数，每个满足（1）中的假设。如记号“”所表示的，我们在接下来的讨论中将假设对于每个，存在唯一的满足（2），这使得映射为从到的良定义映射。
设为有限维子空间，例如，由符合有限元离散化的网格细化参数导致的的子空间；其中表示离散化参数，例如有限元离散化所用网格的网格尺寸。因此，所谓的Galerkin离散化问题（2）如下所示：
同样，我们将在此后进一步假设此离散问题在以下意义上是良定义的：对于每个，存在唯一的满足（2）。为了获得更熟悉的记号，我们可以选择一个的基，其中，并用（参数相关的）系数向量表示。然后，对于每个，，问题（3）变为的个方程的系统：
2.2. 提出的方法
由于（4）是一个可能的非线性方程组，通常包含大量变量，因此对多个求解（4）可能成本较高。为了解决这个问题，我们提出了一种利用神经网络强大逼近能力的方法。更确切地说，我们希望通过神经网络来逼近的映射，其中表示神经网络的网络参数。
训练这个神经网络——换句话说，找到合适的参数，使得成为的一个近似——需要确定一个合适的所谓损失函数。先前的工作 [21] 提出了使用离散方程残差的-范数，例如，确定，其中
这里表示用作训练数据的有限子集。在本文中，我们提出了一种稍微不同的方法，将神经网络的输出表示为，我们通过以下方法确定：
其中分别为左、右预处理器。然后可以通过重新获得所需的。当然，此类预处理器必须仔细考虑，以便考虑到问题的特性。
图1展示了基于有限元的神经网络在实践中实现的示意图。神经网络将参数向量作为输入，并将作为输出。在训练过程中，FEM模块需要直接或分步计算预处理的残差的计算；从这些量中计算得到预处理损失函数。在推理过程中，FEM模块只需提供右侧预处理器，以允许从神经网络输出计算。
为了优化关于神经网络参数的损失函数，通常采用梯度下降法的变体，例如Adam优化器 [27]，或拟牛顿法，例如有限记忆的Broyden-Fletcher-Goldfarb-Shanno (L-BFGS)算法 [28]。这些方法要求计算损失函数关于网络参数的梯度，通过自动微分 [4] 实现；因此FEM模块中的所有计算，例如设置刚度矩阵和/或质量矩阵等，必须自动可微分。在本文中，我们仅选择L-BFGS来实现最小化，因为大多数深度学习中流行的优化器只需要目标函数的一阶导数，即避免了二阶导数的计算，这对于拟牛顿法或牛顿法等二阶方法来说是必要的。作为一种拟牛顿方法，L-BFGS介于一阶和二阶之间，在计算资源和渐近收敛之间提供了折衷。
我们对问题进行预处理的需求很大程度上也是由优化器的选择驱动的。自[29]以来已知，对于二次损失函数，BFGS方法在使用精确线搜索时，产生的迭代结果与共轭梯度（CG）算法相同。由于CG算法的收敛速度高度依赖于Hessian的条件数，因此对优化问题进行预处理是很有意义的。事实上，如果（或者约等于0，对于），则有
(或者分别由这一项支配
因此，我们希望通过选择以某种方式改善在某个极小值点处的Hessian的条件，从而使条件数变得相对较小。即使这些考虑不能构成严格的证明，我们仍然认为它们对我们的方法具有决定性的动机，我们在第3节中的数值观察结果似乎也支持这一点。
2.3 示例1：二维Stokes方程
在以下内容中，我们将展示更多有关所提出的基于FEM的神经网络的实现细节，以解决带有可参数化迎角的翼型周围的Stokes流问题。我们在二维区域中考虑该问题，其中Stokes问题包括找到两个函数，速度和压力，使得
其中函数，以及适当的边界条件将在下文中详细介绍。方程（7）表示动量守恒，而方程（8）表示质量守恒。
使用Galerkin试探空间对系统（7）-（8）的离散化对应形式如下：
其中双线性型和。两个泛函和分别定义为和。相应地定义为方程（9）的残差。为了获得适定问题，离散问题必须满足Ladyzhenskaya-Babuška-Brezzi (LBB) 条件[30, 31, 32]，因此，的选择必须谨慎处理。我们使用Taylor-Hood混合有限元，即对于压力使用连续分片线性试探函数，对于速度使用连续分片二次试探函数。在这种选择下，方程（9）的LBB条件得以满足，从而得到适定的离散问题[33, Lemma 4.23]。
设为的基，其中为对应子空间的维数。当定义且，且时，我们得到残差
其中和。矩阵是对称正定的，系统（10）是对称、非定值的，并表现出典型的鞍点结构。对于这种线性系统，基本的块对角预处理器如下所示：
是所谓的Schur补，详见例如[34]。由于和是对称正定的，我们可以使用它们的Cholesky分解和来定义我们的两个预处理器：
因此，我们得到第2节中引入的预处理残差如下所示：
[35]中的引理2.1表明，简单鞍点问题（10）的条件数被Laplace算子的离散化条件数所控制。因此，经典的FEM理论[33, 定理9.14]表明，对于Lagrange有限元，正如我们选择的速度离散化，条件数与成比例。因此，通常需要对简单鞍点问题进行预处理。事实上，如附录A所示，预处理刚度矩阵的平方条件数（12）与离散化无关地被约束。在附录A中，甚至还展示了该矩阵的不同特征值的数量为三个。
再次考虑BFGS在完美线搜索下等价于用于二次优化问题的CG，[29]，BFGS应用于具有平方预处理矩阵作为其Hessian的此类二次优化问题，最多在三次迭代中收敛。当然，这些性质不直接适用于我们的设置，因为所考虑的问题是高度非线性和非凸的，因为存在神经网络。然而，应用预处理器到问题的线性部分至少会逐步改善该情况下的整体性质，这样的假设似乎是合理的。事实上，我们在第3节中的数值示例支持了这一点。
现在让我们转向显式的数值示例，我们求解具有源项和非齐次Dirichlet边界条件的Stokes方程（7）-(8)。我们定义域，具有长为1的中心翼型NACA 0012[37]。在入口边界上使用与迎角有关的Dirichlet边界条件，在出口边界上使用Neumann边界条件，并在翼型边界上使用无滑移边界条件（见图2）。
在 上，， 在 上，， 在 上，。
请注意，出口边界上的齐次Neumann边界条件会自动将平均出口压力设为零。
图2展示了实验中使用的三角网格。正如我们所预期的，翼型附近的网格比外缘更精细。我们使用了个速度自由度和839个压力自由度，总共对应于7225个自由度。
齐次Neumann边界条件被称为所谓的自然条件，因此通过构造包含在（10）中。然而，在汇总残差之后，仍然需要对速度强制执行非齐次Dirichlet边界条件。为此，我们利用所选基函数是Lagrange型的，即，每个基函数都有一个关联的Lagrange节点，使得
现在通过调整右端向量和（10）中的刚度矩阵来强制所有这些满足Dirichlet边界条件：首先，将刚度矩阵中对应于边界节点的列乘以指定的边界值，然后从中减去结果。接下来，我们将中对应于边界节点的值设置为边界值。在刚度矩阵中，我们修改与边界节点相对应的行和列，使得对角线值为1，非对角项为0。注意，这是在预处理之前完成的，因为这些修改不会影响系统（10）的对称性。
在当前情况下，预处理的刚度矩阵和残差（10）的右端向量与神经网络的预测无关。因此，它们只需计算一次，并且可以在神经网络训练的每次损失计算中重复使用。
2.4 示例2：二维稳态Navier-Stokes方程
为了避免冗余，我们关注从Stokes方程扩展到Navier-Stokes方程所需的变化。因此，未讨论的细节仍然与上一节相同。Navier-Stokes方程如下：
即，相较于Stokes方程（7）-(8)，特别是加入了非线性对流项。此外，我们引入了所谓的运动粘度作为附加参数。离散化得到
其中与之前相同，与之前相同，但乘以，三线性型来源于非线性项。我们再次使用Taylor-Hood单元进行非线性情况下的离散化，因为在线性情况下已经需要它们以确保稳定解。为了形成相应的非线性代数方程，我们利用与（10）相同的Lagrange基函数，并得到离散残差
其中。我们使用与Stokes问题相同的预处理，因此得到预处理系统
同样，我们使用神经网络直接预测和。在我们的数值示例中，我们将求解与上一节相同几何形状、网格和边界条件下的Navier-Stokes方程（14）-(15）。
由于Navier-Stokes问题是非线性的，刚度矩阵依赖于解，即我们在此问题中的神经网络的预测。这意味着与（12）不同，系统矩阵和右端向量可以提前计算一次，而现在在神经网络训练期间的每次L-BFGS迭代中，我们必须重新计算残差（18）的一部分。更确切地说，这适用于对流项和由于包含Dirichlet边界条件而影响的右端向量，如下文所述:
即，可以通过与独立于的3D-张量的收缩来计算，因此可以预先计算并重复使用。由于是稀疏的，可以以一种内存高效的方式完成。
现在我们简要讨论在残差（17）中强制Dirichlet边界条件。除了非线性部分之外，这可以像处理Stokes方程一样处理。唯一的区别是，对于，我们实际上必须计算，其中使用稍微修改过的参数，其中相应的边界条件已经在Dirichlet边界上的所有Lagrange节点上强制执行，即：
其中表示在边界节点处规定的值。对于Stokes方程，边界条件在预处理之前包含在（17）中。通过这些更改，可以按照上述描述进行基于FEM的神经网络的训练。
3. 数值结果
在下文中，我们展示了所提出的基于FEM的神经网络用于求解第2节中描述的带有可参数化迎角的翼型周围Stokes和Navier-Stokes流的结果。特别是，我们将重点关注使用预处理器和不使用预处理器之间的比较。
在PyTorch中实现了各自的神经网络和第2节、2.3节和2.4节中描述的FEM模块[38]。除非另有说明，否则神经网络为普通前馈网络，由五个隐藏层组成，每层有50个神经元，输出层的大小等于离散化的自由度数。除最后的线性层外，我们使用SELU[39]作为激活函数。为了促进高效训练，我们将网络权重设置为Xavier的正态初始化[40]并缩放为。在训练中，我们使用PyTorch的L-BFGS实现。为了使结果更具可比性，我们选择，确保每次迭代只执行一步L-BFGS。此外，我们选择L-BFGS中的线搜索以寻找满足强Wolfe条件的解。
在下文中，“误差”一词不指离散化误差，而是指我们的神经网络解与使用相同网格和Taylor-Hood离散化的纯FEM参考解之间的误差。后者使用FEniCSx计算[41, 42]。
3.1 非参数情况
我们从只有一个输入和一个输出层的线性神经网络进行的具有单个迎角作为训练数据的实验开始。考虑这种非参数化和高度欠定的情况使我们能够评估预处理的好处以及其他效果。
对于和40000次L-BFGS迭代，分别训练和评估预处理和非预处理神经网络，针对单个迎角求解Stokes和Navier-Stokes方程。在后一种情况下，解是通过粘度和计算的。图3表明，在所有这些情况下，预处理都显著加快了训练速度，极大地减少了迭代次数。对于预处理的情况，我们甚至可以观察到在数值精度相关水平上残差的停滞；这可能与第2.3节中（12）下方段落所述的类似原因有关。需要注意的是，训练在运行时间方面也要快得多，因为迭代次数的显著减少远远抵消了由于涉及预处理的几项附加操作而导致的单次迭代运行时间的微小增加。
表1展示了在-范数和-范数下达到的相对误差。同样，预处理的好处在最终误差减少几个数量级方面清晰可见。特别地，注意对于（这对应于增加的Reynolds数），只有使用预处理才能达到令人满意的精度。如图4所示，最终神经网络解和纯FEM解之间几乎没有可见的差异。
3.2 参数情况
在本节中，我们考虑真正的参数化问题，并用多个迎角进行训练，即3, 5, 9, 17个不同的，这些值在范围内等距分布。使用最多4000次L-BFGS迭代来训练各自的网络。
图5显示了在训练过程中针对所考虑的不同配置的损失值。即使在这种情况下，预处理也能够将损失的衰减加速两个数量级以上，无论是对于Stokes方程还是Navier-Stokes方程。对于Stokes问题，我们还在更细的网格上进行了实验，该网格是通过一次统一细化从原始网格获得的，结果是28085个自由度。原始网格和更细网格上的损失行为可以作为网格独立性的指示，这与我们在附录A以及前一节中的讨论是一致的。
接下来，让我们考虑训练误差，即训练数据上的误差，以及泛化误差，即插值输入数据和外推输入数据上的误差。为了评估插值和外推质量，我们测试了在训练范围内未见的迎角以及训练范围外的。图6显示了训练数据（“训练”）上的平均误差，以及-范数和-范数下泛化误差（“插值”和“外推”），对于训练集的不同大小。为了清晰起见，我们只报告速度在方向的误差，因为方向的速度误差和压力误差表现相似。很明显，预处理不仅大大提高了代理模型在训练集本身上的质量，还提高了其插值和外推能力；事实上，对于未经过预处理的问题，随着训练集大小增加而获得的误差远远超过可接受的范围。对于预处理问题，训练集大小增加而导致的误差增加可能是因为使用了相同的网络结构（即相同自由度数）来解决日益复杂的问题。如预期，随着训练数据的增加，泛化误差减小。
我们的研究结果可总结如下：使用预处理的损失函数加速了线性和非线性问题的训练过程，并且提高了所得代理模型在训练和泛化误差方面的质量。
4. 应用：求解与二维稳态Navier-Stokes方程相关的逆问题
在下文中，我们通过将其应用于一个原型逆问题来展示我们方法的相关性，该逆问题与二维稳态Navier-Stokes方程有关。这个问题对应于一个应用场景，即通过翼型上的两个传感器的压力测量值来确定迎角。由于在实际环境中的测量总是存在噪声，能够量化由这些噪声带来的不确定性至关重要。我们的目标是估计与这些观测相关的未知迎角的概率密度。因此，我们采用了一种概率方法，并结合使用马尔可夫链蒙特卡罗（MCMC）方法和完全可微分的基于FEM的神经网络。
设为翼型上两个固定点的压力观测值。我们假设两个压力传感器分别位于翼弦顶部的8.18%和99.94%的位置；见图7。我们还假设测量值是独立的，并根据标准差的正态分布进行随机分布。由此，我们得到一个对数似然函数，除了某些与无关的常数外，给定如下：
我们故意使用来表示仅在这两个测量点处评估的预处理神经网络的结果。将此似然与均匀先验结合起来，我们对结果的后验分布感兴趣。
我们故意使用来表示仅在这两个测量点处评估的预处理神经网络的结果。将此似然与均匀先验结合起来，我们对结果的后验分布感兴趣。
由于的计算对表现出非线性依赖性，我们采用哈密顿蒙特卡罗（HMC）方法[43, 44]，更确切地说是No-U-Turn取样器（NUTS）。后者是一种基于梯度的MCMC方法，能够逼近所需的密度函数。HMC允许通过使用近似哈密顿动力学模拟从此后验分布中进行取样，然后通过Metropolis接受步骤进行校正。为了解决哈密顿微分方程，必须计算目标后验密度函数关于的导数。在这里，基于神经网络的代理建模的一个重要优势显现出来：代理模型不仅是可微的（理论上这也可能适用于其他技术），而且这些导数可以使用AD有效地计算。我们使用了PyTorch库Pyro中的NUTS实现[45]。
接下来，我们希望测试这种方法。实验的详细信息如下：图7（右）显示了针对17个等距分布在范围内的迎角的传感器位置的正向问题。使用此基于FEM的神经网络作为代理模型，我们通过在相应的预测上添加标准差为0.5的高斯噪声生成了10个人造噪声压力“测量值”，这些预测是在参考迎角处进行评估的。“测量值”如图8所示。
图8显示了由MCMC生成的后验分布的1000个样本，采用了1000步的暖启动。
生成的分布的经验均值为5.03，标准差为0.32。因此，MCMC和基于FEM的神经网络的结合能够为此测试用例提供合理的结果。
5. 结论
基于FEM的/受FEM启发的神经网络被认为结合了两种方法的最佳特性：有限元法和神经网络。一方面，所考虑的（参数化）偏微分方程通过一种被很好理解且可靠的方法进行离散化，边界条件和复杂几何的引入相对直接。另一方面，离散化受益于神经网络的强大能力来处理参数依赖性，例如，通过获得廉价评估且完全可微的参数到解的映射，这在逆问题的背景下是有利的。
在本文中，我们首次将这种技术应用于流体动力学中的偏微分方程，这些方程在离散化后产生鞍点问题。特别地，我们介绍了一种修正的方法，其中对残差损失进行预处理有助于提高训练速度和结果精度。我们通过数值示例证明了我们方法的有效性，这些示例包括一个原型但现实的流体动力学模型问题及其相关逆问题。
在本文中为我们的方法提供了概念验证后，一个有趣的未来研究目标是探索包括分布式并行方法或低成本预处理器在内的可能性，这些在大规模有限元实现中是典型的。另一个目标是考虑具有不同特性的其他偏微分方程或偏微分方程系统（包括合适的预处理器）。特别是，扩展到需要稳定化的高雷诺数流动问题或时间相关问题。
💙整理不易，希望各位道友能够多多支持宝库，支持邪云宝库！你的一个点赞、一次转发、 随手分享，都是宝库前进的最大动力~
💛2024，不忘初心，宝库会给大家带来更好的内容，让我们2024，一起暴富！
