基于GPU的快速有限元法求解密度场

李果阳，严华，*，张征宇，陈沁梅，祝福顺

(1.四川大学 电子信息学院，成都 610065；2.中国空气动力研究与发展中心 高速空气动力研究所，绵阳 621000）

对于风洞试验，空气流动包含着最重要的特征信息，对处于空气流动形成的流场中物理模型的密度场进行定量分析一直都是流场的重要研究方向［1］。传统的纹影技术由于在定量测量上效果不好而常常作为一种定性分析的方法，干涉技术由于需要严苛的要求才能实现，缺乏实用性。因此，背景纹影技术(Background Oriented Schlieren，BOS）就在这种情况下被Meier［2］提出。运用该方法与风洞中的视频测量技术相结合可以较简单地实现对处于风洞试验中的物理模型进行密度场的视频测量［3-6］，通过测出一束非平行光在流场扰动下的偏折角，进而建立偏折角与扰动流场折射率间的量化关系，得到关于密度场的二阶偏微分方程。利用视频测量法得到关于坐标及坐标偏移量的离散数据，再对于这样的密度场中二阶偏微分方程求解出非平行光经过扰动流场后对应位置的密度投影值。目前有多种常用数值模拟计算方法，经过比较，有限元法因擅长处理各种复杂区域、精度较高及适于并行计算等优点［7-8］，因此本文将采用有限元法求解密度场对应位置的投影值。

从有限元的计算方法可以看出，在用有限元法进行大规模数值计算时，不仅要考虑多种力学行为的全耦合求解问题，计算过程复杂，而且总刚度方程的求解时间占据了整个计算过程很大的比重［9-10］。由于神经网络具有强大的映射、实时和并行计算能力［11］，为了解决有限元中的总刚度矩阵求解耗时，简化有限元法求解的复杂度问题，本文不再采用传统有限元计算方式，而是引入了神经网络进行求解。另外，风洞试验中需要实时分析风洞试验密度场的分布变化，由此需要快速计算密度场。但是，神经网络训练不仅非常耗时且会占据计算机大部分资源，而且有限元计算过程中产生的大规模稀疏线性方程组数据量也很巨大［12］，计算非常耗时，严重影响了密度场的求解效率。对于风洞试验密度场这样的时变系统要求实时分析，其核心就在于能够实时计算，因此，非常有必要探寻求解密度场的快速计算方法。目前主流的高效计算方法主要采用基于GPU的并行计算方式［13-15］，因此本文主要就针对求解过程中耗时的神经网络、刚度矩阵和载荷向量的求解进行并行加速，以实现密度场的快速求解。

为此，本文提出了一种高效求解密度场的思路，给出了神经网络与传统有限元分析融合的有效方法，以及利用GPU快速求解密度场的计算方法和具体步骤，实验结果证明了所提方法的有效性。

1 有限元法求解密度场

1.1 密度场的二阶偏微分方程

利用背景纹影技术的方法测量风洞试验中的流场［16-17］。坐标系的x轴为逆气流方向，y轴与x轴垂直，指向试验段上壁板，z轴与x轴垂直，与光线传输方向相反。为了便于分析，不妨假设流动介质的折射率仅与y轴方向的坐标值有关，且随y轴方向坐标值的增加而增大，沿z轴方向平行射入的光线，在y-z截面内发生偏折，光线偏折角分析如图1所示。

图1 光线的偏折角分析示意图Fig.1 Schematic diagram of ray deflection angle analysis

折射率n与气体密度ρ之间的Gladstone-Dale公式为

由光线传输方程和式(1），可得y轴方向偏折角εy与气体密度ρ之间的关系为

同理可得，x轴方向偏折角为

式中：L为风洞试验段的边长；KGD为Gladstone-Dale常数，在空气中标准状态下取值为0.226 cm3／g。

对式(4）和式(5）中y、x分别微分再相加即可得到密度场应满足的二阶偏微分方程：

式中：ρm为密度场密度。

将式(6）的等式左端记为Δu，等式右端记为

其求解域为空腔模型，记为Ω。

考虑式(6）的变分问题，对于任意的v∈HΩ，同时乘上式(7）两边，使用格林公式，利用边界条件可得

可证明得到，在一定连续性要求下，式(6）可等价于求u∈V，使得a(u，v）＝g(u）对任意的v∈V成立。

考虑上述变分问题的有限维逼近，即在有限维子空间Vh⊂V下，求uh∈Vh，使得a(uh，vh）＝g(uh）对任意的vh∈Vh成立。设在Vh下的一组基函数为，并依次取vh为每个基函数，则可得

分析式(11），其本质是一个线性方程组：

基于上述分析，在初边值的条件约束下，由Lax-M ilgram定理可知，式(6）即可等价于有限元法下求解一个线性方程组［18-19］：

式中：K为单元力和单元位移关系间的系数矩阵，代表了单元的刚度特性，称为单元刚度矩阵；F为结点产生单位位移时，在结点上所需要施加的结点力，称为结点载荷；δ为结点受到扰动流场作用后，结点在其自由度方向偏移的距离，称为结点位移，且K、F、δ的具体表达式如下：

式中：D为弹性矩阵；t为单元厚度；B为单元应变矩阵；J为雅可比矩阵；α1、α4为刚体位移；α2、α3、α5、α6为单元的常应变；x、y为结点坐标。

1.2 密度场的有限元法求解过程

传统的有限元分析是利用变分原理将求解域剖分求解的一种数值计算方法，具体步骤如下：

1）连续体离散化。用有限个离散单元体集合代替原有的连续体，也称网格剖分，即进行单元划分，将全部单元和结点按一定顺序编号，每个单元所受荷载均按静力等效原理作用在结点上，并根据实际情况在位移受约束的结点上设置约束条件。

2）单元分析。建立各个单元的结点位移δe和结点力Fe之间的关系式Fe＝Keδe，Ke表示单元刚度矩阵，即将单元的结点位移作为基本变量，确定一个近似表达式，再应用流体力学理论和虚功原理得到单元中结点力与结点位移的关系式。

3）整体分析。对各个单元组成的整体进行分析，其目的是根据一定法则建立一个线性方程组，得到全部结点载荷F与全部结点位移δ的关系式Kδ＝F，从而求解出结点位移。

本文所实现的快速有限元法采用了神经网络代替单元分析中求解结点力与结点位移关系的复杂计算，从而使得整体分析计算过程简化，可大大缩短整体计算时间。

1.2.1 网格剖分

利用有限元法求解泊松方程首要的就是将求解域进行离散，划分得到有限个互不重叠的单元网格，同时，对这些单元和结点进行编号，这些互不重叠的网格单元通过结点连接相互影响。在划分单元网格时，对计算结果有着重要影响的是网格单元尺寸大小、网格的疏密程度及网格的类型。三角单元网格能够适应复杂的求解域，单元大小划分十分方便，但计算精度相对较低，而四边形单元划分不易描述非正交的直边界和曲边界，但其计算精度相对较高。

本文采用的是四边形单元和三角形单元结合的方式，实现了一种四边形网格单元与三角形网格单元相结合的方式，剖分示意图如图2所示。在求解时为便于网格的剖分实现，采用了先对求解域进行等参四边形单元剖分，再将一个四边形单元一分为二，形成两个等参的三角形单元剖分算法，然后对它们进行单元编号和结点编号。

图2 网格剖分平面示意图Fig.2 Grid subdivision p lan sketch

1.2.2 神经网络拟合

利用神经网络代替传统有限元法中的单元分析，关键在于要清楚单元分析得到的单元刚度矩阵的特性。

对于本文所采用的三角单元剖分，根据弹性力学二维有限元法［18］，可知有限元法中单元刚度矩阵有以下几点性质：

1）单元刚度矩阵是对称方阵。单元刚度矩阵每行元素之和为0，由对称性质可知，每列元素之和也为0。

2）对于平面单元，单元刚度矩阵不会随着几何参数改变。

3）根据单元刚度矩阵的物理意义可以得出其只与单元结点、相对坐标及形函数等相关，与结点绝对坐标无关。

由以上单元刚度矩阵的性质结论可推出，刚度矩阵的计算实质是从单元的尺寸、材料性质到刚度矩阵的映射问题，而神经网络的优势在于：可以方便地建立输入和输出间复杂的映射关系。因此，利用神经网络可以很方便地从单元结点坐标映射到单元刚度矩阵元素。

本文中，由于密度场二阶偏微分方程右端项是关于光线偏折的一阶微分式，且实验数据是使用视频测量方法在没有流场扰动时，拍摄一幅包含高密度标记点的背景板图像作为参考图像，然后加入扰动流场，再拍摄背景板的时序图像，根据视频测量的标定原理测得背景板标记点在坐标系下的(x，y）坐标值及其偏移量Δx和Δy。因此在求解时需要对数据进行预处理，即根据视频测量法原理和光线偏折理论得到坐标偏移量和偏折角的关系，再对其做归一化处理，接着采用神经网络对偏移角拟合得到二阶偏微分方程右端项一阶微分式，即通过对右端项的一阶微分式在点(x，y）处根据泰勒级数展开式可以得到

式中：h为坐标点更新的步长。

神经网络拟合模型如图3所示。将经过预处理的标记点坐标数据(其中80%为训练样本，20%为验证样本）作为网络输入，通过一个含有400个神经元的隐含层，激活函数采用tanh型，优化函数采用目前普遍使用的Adam进行有监督的学习。将网格单元结点通过该训练好的神经网络，最终实现网格单元的结点通过神经网络映射到单元刚度矩阵元素。

图3 神经网络拟合模型Fig.3 Neural network fitting model

1.2.3 整体分析

总刚度矩阵K由单元刚度矩阵Ke组装而成，计算如下：

式中：Ce为单元选择矩阵。根据单元的结点编号和结构中结点编号的对应关系，将单元刚度矩阵进行式(19）迭加组装成总体刚度矩阵，得到结点载荷与结点位移间的关系式，如式(13）所示。

然后进行狄利克雷边界条件约束，采用直接法求解得出结点位移。

2 基于GPU的有限元法并行求解

本文中将有限元法求解密度场的方法采用并行的思想，对求解过程极其耗时且适合并行计算的神经网络、总刚度矩阵和总载荷向量的求解进行并行加速。融合了神经网络的单元分析过程，以时间代价简化了求解复杂度，为进一步加快求解速度采取了数据并行、多线程同时计算的方式。整体分析中，总体刚度矩阵和总载荷向量的求解进行分块多线程并行计算。

2.1 CUDA架构介绍

GPU具有高性能的浮点计算能力，CPU具有低延迟、逻辑计算能力强的特点，CUDA编程正是利用GPU和CPU这种特点来协同处理任务。运行流程则是在CPU端准备数据，再传到GPU端并行计算，最后将计算结果回传到CPU端。CPU和GPU的这种协同工作方式大大提高了计算性能［20-21］。图4为CUDA线程组织结构模型。

为了明白CUDA编程如何进行线程组织，必须了解其结构模型。从图4中可以看到，CUDA中线程层次分为3层，分别是线程网格(grid）、线程块(block）、线程(thread）。1个线程索引(thread Idx）可以是1～3维的，进而多个thread组成的block也可以是1～3维的。如同线程组成线程块，线程块也能组成1～3维的线程网格。执行CUDA应用程序时，整个应用程序的关键是kernel函数，它是以grid的形式组织，以block为单位，各block间并行执行，不同block间的数据共享只能通过全局显存，对于同一block内的不同thread间则可以通过共享内存进行通信，即在kernel函数中存在2个层次的并行，即block之间并行计算及thread之间并行计算。

图4 CUDA线程组织结构模型Fig.4 CUDA thread organization structure model

另外，为了达到更高的效率，在CUDA 编程应格外关注内存的使用。在CUDA编程时要尽量使用寄存器，尽量将数据声明为局部变量。而当存在着数据的重复利用时，可以把数据存放在共享内存里。而对于全局内存，需要注意用一种合理的方式来进行数据的合并访问，以尽量减少设备对内存子系统再次发出访问操作的次数。

2.2 神经网络并行加速方法

本文中将有限元法求解密度场的方法采用并行的思想对求解过程极其耗时的部分进行加速处理，其模型如图5所示。数据拟合的并行化主要是针对神经网络的学习训练过程。首先，将在密度场测量得到的离散数据分为训练样本、验证样本和测试样本，并将其暂放在内存中，然后对神经网络进行初始化参数，并且将训练样本分成多个批量(batch），加载到GPU显存中。在神经网络正向传播过程中，将原本一整个数据分为多片，调用CUDA的kernel函数为每个batch数据分配处理线程。由于在CUDA中一个W rap线程束有32个线程，在进行线程分配时应当尽量选择32的整数倍，但是CUDA中一个线程块里线程数量又不多于1 024，因而可以尽可能提高线程的利用率。从主机端调用kernel函数在GPU设备端对分片好的训练样本在大量线程中被同时并行处理，这就极大地加快了正向传播过程。同理，它的反向传播过程也是各线程中进行一系列的矩阵运算和向量运算。当每个线程所传递的误差信号都到达网络的隐含层，就能得到权值偏导和偏置偏导，而每个batch上得到的是误差信号均值，并且只对其本身参数偏导进行更新，所有线程间互不干扰，所有的参数间也不相关，这就完成了一次并行迭代过程。

图5 神经网络并行模型Fig.5 Neural network parallel model

2.3 总刚度矩阵组装的并行化

2.3.1 总刚度矩阵的并行性分析

对于大规模的有限元分析中，串行的有限元过程中均不是直接在内存中一次建立总的平衡方程组，因为对于大规模的有限元问题，其结点总数也很大，导致总刚度矩阵的阶数太大，从而计算机计算过程中内存不够。因此，串行有限元过程中一般有2种处理方法：一是将总刚度矩阵分块存储在外存中，二是边单元分析边消去，避免总刚度矩阵的生成。但是这2种方法都需要内外存之间频繁地进行数据交换，因此对于结点数较多的大规模问题计算一次是非常耗时的。

总刚度矩阵是将所有单元刚度矩阵按照其结点编号对号入座组装而成，所有单元的单元刚度矩阵不能向总刚度矩阵一次性迭加而成，否则就可能在计算过程中发生计算错误和存取冲突。进一步分析可知，总刚度矩阵组装最大的并行只可能发生在每次同时迭加互相没有共同结点的多个单元的刚度阵所有元素，即实际采用的方案可以是每次同时迭加一个单元刚度矩阵多个位置处的元素，或者每次同时迭加多个单元刚度矩阵在同一位置处的元素。另外，总刚度矩阵的组装与其存储格式密切相关，大型有限元问题产生的总刚度矩阵一般是稀疏对称的。

2.3.2 总刚度矩阵的并行组装

通过上文分析，本文中的问题则是一个稀疏对称的总刚度矩阵，对于这样的总刚度矩阵组装过程很耗时但却非常适合并行计算。总刚度矩阵由每个网格单元的单元刚度矩阵组装而来，每个网格单元的单元刚度矩阵的值只与其形状大小、对应点次序及在整体坐标系的方位有关，因此本次求解划分的矩形网格单元内的三角单元具有相同的单元刚度矩阵。图6表示了划分的任意一个网格位置结点标号。

图6 剖分的任一网格位置结点编号Fig.6 Node number of any grid position of subdivision

本文单元刚度矩阵组装成总刚度矩阵的具体方法是：将单元刚度矩阵Ke的每个子块Kij送到总刚度矩阵的对应位置上去，然后进行迭加即可得到总刚度矩阵K的子块，对每个单元都进行如此操作计算即可得到总刚度矩阵。

在执行并行组装总刚度矩阵时，关键是如何找出Ke中的子块在K中的对应位置，然后为每个Ke的组装分配线程索引，这样在并行化过程中给每个单元分配的线程求解得到的单元刚度矩阵才能正确地组装成总刚度矩阵，解决这个问题需要清楚单元中结点编码(局部码）和整体结构中的结点编码(总码）之间的对应关系。图7中每个单元的3个结点按逆时针方向编为i、j、m，单元刚度矩阵中的子块是按结点的局部码排列的，而总刚度矩阵中的子块是按总码排列的，因此，在单元刚度矩阵中把结点的局部码换成总码，并把其中的子块按照总码次序重新排列，以单元U1为例，局部码i1、j1、m1对应于总码m、m＋1、m＋x＋1，即单元刚度矩阵和总刚度矩阵位置对应关系如图7表示。其中上三角网格单元和下三角网格单元的单元刚度矩阵相同。

图7 单元刚度矩阵与总刚度矩阵位置对应关系Fig.7 Location of element stiffness matrix corresponding to total stiffness matrix

另外，有限元法产生的总刚度矩阵的阶数一般很高，在解算时矩阵的阶数很高和存储容量有限是矛盾的，因此有必要知道总刚度矩阵的特性来解决节约存储容量的问题。经分析，总刚度矩阵的特性如下：

1）对称性。只存储总刚度矩阵的上三角部分，可以节约一半的存储容量。

2）稀疏性。总刚度矩阵的元素绝大部分都是零元素，非零元素只占一小部分。

3）带形分布。总刚度矩阵的非零元素分布在以对角线为中心的带形区域内，在半个带形区域内，包括对角线元素，每行具有的元素个数叫做半带宽，即可以只存储该部分元素值来节约存储空间。

若存储时不采取有效的存储方式，不仅需要大量的存储空间、时间，还会增大计算量。因此，本文在实现过程中根据总刚度矩阵的特性采用了较为常用的CSR矩阵压缩存储格式，用3个一维数组分别存储稀疏矩阵中上三角的半带宽非零元素值、对应的列号及行偏移，使得GPU并行更易实现。

在并行实现时，GPU运行kernel核函数，每个CUDA线程执行一个网格单元刚度矩阵写入总刚度矩阵的任务，主要分为如下步骤：

步骤1根据参数计算出网格单元的刚度矩阵。

步骤2根据网格单元的结点编号，判断网格单元是上三角形单元还是下三角形单元。

步骤3根据判断结果将值写入总刚度矩阵对应的位置。

为能够并行计算刚度矩阵和载荷向量并且不产生线程和存储冲突，本文根据不同单元网格结点编号进行固定空间的索引存取，即每次同时迭加所有单元的单元刚度矩阵同一位置的元素。总刚度矩阵总载荷向量和结点编号密切相关，因此在GPU并行求解总刚度矩阵和总载荷向量时，最为关键的是如何将各结点写入到GPU存储空间。对于本次求解的网格单元，每个结点的载荷在F中最多会迭加6次。为避免同时迭加或者使用CUDA原子操作，CUDA编程时为每个结点分配一个大小为6的空间，即最开始的总载荷向量F大小为N×6。实现中要保证下三角结点的编号顺序始终为(m，m ＋x＋1，m ＋x），上三角顺序(m＋x＋1，m，m＋1），把所有下三角的结点编号按顺序设为索引1、2、3，同理把所有上三角的结点依次设为4、5、6，如图8所示。相同索引在同一结点处只会出现一次或者不出现，因此在求得结点之后先将载荷向量写入对应结点对应的索引处，如对于m点，对应的位置为6×m＋Index(索引），这样即可线程和存储不冲突，从而正确求解。

图8 结点存储位置Fig.8 Node storage location

当所有值全部写入上述数组之后，则每个结点处真正的载荷值为6个值相加，此时每个CUDA线程只需要将6个数据相加便可得到最后的总刚度矩阵和总载荷向量。

3 实验验证

3.1 实验环境

运用上述串并行实现原理方法，进行了密度场有限元法求解。实验平台操作系统采用Ubuntu 16.04，CPU为Intel i7-6800k，GPU为NVIDIA GTX 1080Ti。实验数据由某风洞研究所提供。对比串行和并行实现下不同的网格精度、各主要模块的计算耗时来验证并行快速有限元法的有效性。

3.2 串行速度

表1列出了不同精度大小的网格单元剖分下主要模块串行求解时间。图9展示了有限元法同网格精度下主要模块占总运行时间的比例。根据图9结果进行分析，在不同网格规模下，有限元法求解密度场过程中网络训练和网络预测时间开销相差较小，但是随着网格规模加大，求解总刚度矩阵和总载荷向量的时间开销却成倍增加，从求解总时间来看这整个过程极其耗时，完全无法达到风洞试验分析场密度实时要求。

表1 不同模块下CPU串行求解时间Table 1 CPU serial solution time under different modules

图9 不同网格精度下主要模块占总运行时间的比例Fig.9 Proportion ofmain modules in total running time with different grid precision

3.3 并行速度

为了解决不能实时计算密度场的问题，本文着力于寻求一种基于GPU的快速求解方法，使用现在流行的GPU高效计算硬件平台，极大提高了求解效率。分别进行了如表2所示几组实验。图10展示了有限元法不同网格精度下串行求解和并行求解密度场的加速比；图11展示了有限元法不同网格精度下串行和并行求解密度场主要模块的加速比。本文所指的加速比皆是指串行求解和并行求解运行时间的比值。

图10 不同网格精度下求解密度场的加速比Fig.10 Acceleration ratio of solved density field with different grid precision

图11 不同网格精度下求解密度场主要模块的加速比Fig.11 Acceleration ratio of main modules of solved density field with different grid precision

表2 不同模块下GPU并行求解时间Table 2 GPU parallel solving time under different modules

从不同的网格划分疏密程度的实验结果对比可以看到，剖分的网格较为稀疏时求解密度场总体加速比较小，大概在25.6倍左右，随着剖分网格越密，总体加速比都在逐渐增加，这是由于网格精度增加，有限元法产生的稀疏矩阵规模越大，利用GPU并行求解总刚度矩阵和总载荷向量的优势越明显，而网络训练的加速比改变不大是因为导入的实验数据量并没改变。另外，可以看到求解过程中的求解总刚度矩阵和总载荷向量占总运行时间的比例逐渐上升，导致这一现象的主要原因是随着网格数量增加，总刚度矩阵的组装和总载荷向量的累加操作需要大量的共享存储器，每个流处理器同时运行的block快速减少，从而GPU的利用率下降，运行时间增加。

3.4 精度比较

本文中采用L2范数和L∞范数来计算误差，计算如式(20）和式(21）所示，分别表示本文快速有限元法数值计算结果密度xi与精确结果x＾的均方误差，以及本文快速有限元法数值计算结果密度xi与精确结果绝对值的最大误差。

式中：Δh为剖分的网格单元面积。

从图12可以看出，有限元法串行和并行求解密度场的精度随着网格划分变密其精度增加，但精度增加缓慢，且串行计算精度高于并行计算，这是因为并行计算加入了一些累加计算，所以截断误差影响更显著。另外，该方法精度在绝大多数应用场景中都处于可接受水平，如在网格精度为496×450时，其均方误差和最大绝对值误差均在10－5量级下，在绝大部分场景下，这一误差可以完全忽略。

图12 有限元串并行精度比较Fig.12 Serial and parallel precision comparison of finite element method

4 结论

1）给出了基于GPU实现的快速有限元求解风洞试验密度场的方法。通过对比实验分析对串行方法中耗时较多的神经网络拟合、总刚度矩阵组装及总载荷向量求解实现了GPU并行加速。实验表明，随着网格数量的增加，其加速比可达数百倍以上，且其精度达到了工业要求，提高了有限元法求解密度场的时效性。

2）实现的快速有限元法相对于传统的有限元简化了对刚度矩阵的求解，大大减少了计算时间。

3）随着网格数量增加，加速比相应增加，其精度也逐步提高，但求解更耗时。