龙格库塔间断有限元方法求解二维欧拉方程的多GPU加速实现

Accelerating the Runge-Kutta Discontinuous Galerking Method for Solving Two-Dimensional Flow on Multi GPUs

摘要: 为解决龙格库塔间断有限元方法(RKDG)求解流场耗时的问题，本文应用二维NACA0012翼型作为测试算例，使用多GPU加速求解。将流程网格按照GPU个数进行剖分，每个GPU计算一个网格区域。各计算节点上设置核函数的线程数等于流场网格数，节点间的数据通信使用MPI (Message Passing Interface)。通信过程中采用CUDA流和MPI非阻塞操作以覆盖数据的传输和计算，减少通信代价。结果表明，与CPU串行程序相比，1个、2个、4个GPU上分别获得了33倍、59倍和108倍的加速比。

Abstract: It is time-consuming to use Runge-Kutta discontinuous Galerkin method to solver flow field. In this article, we use multi GPUs to accelerate computing of two-dimension NACA0012 airfoil. The flow mesh is divided into several blocks according to the number of GPUs. We specialize kernels with a one-element-per-thread strategy, and use MPI to communicate data among computing nodes. In the process of communication, we use CUDA stream and MPI non-blocking operation to overlap computation and communication. The result shows when compared with the serial CPU program, the speedup ratio of GPU code running on one, two, four GPUs is around 33, 59 and 108.

1. 引言

计算流体力学(CFD)在航空工业中扮演着越来越重要的角色。随着近些年的发展，CFD能够模拟出许多现实的场景。但是对于复杂系统，需要采用高精度的算法。龙格库塔间断有限元算法是一种高精度的流场算法，具有能够捕捉间断，处理网格悬点等优点 [1] 。然而由于计算量大，往往需要花费很长的时间完成流场求解。

人们提出了许多方法来减少CFD求解时间，并行计算就是其中一种很好的方法。CPU并行加速在工业应用上已经非常成熟，然而建立大规模CPU并行计算集群非常昂贵。因此GPU是一种更好的选择。近年来，GPU在科学计算上表现出了重大的潜力。CUDA是Nvidia公司2007年推出的GPU并行计算平台和编程模型 [2] ，实验表明，能够显著提高计算性能。

Klöckner A等人在单个GPU上构建了高效的高阶间断有限元解算器，并在三维四面体网格上求解了麦斯威尔方程，与同一代的CPU相比达到了近50倍的加速比 [3] ；何晓峰等在单个GPU上，使用RKDG方法对二维欧拉方程进行了加速 [4] ，与CPU相比获得了25倍的加速比；Dawei Mu等人在单个GPU上，使用间断有限元方法对三维弹性地震波方程进行了加速 [5] ，并与CPU并行程序进行了比较，在单精度下相比1个、2个、4个、8个CPU并行程序获得的加速比分别为12.9、6.8和3.6。本文在此基于多GPU对RKDG方法进行加速。论文结构安排如下，我们首先介绍了龙格库塔间断有限元方法和CUDA程序概述，然后介绍了单GPU以及多GPU上基于CUDA的间断有限元方法求解，最后分别在一个、两个和四个GPU上进行了求解，统计了耗时，并与CPU结果进行了比较。

2. 龙格库塔间断有限元求解

本文中，我们考虑二维双曲方程，格式如下：

$\{\begin{array}{l}{U}_t+F{\left(U\right)}_x+G{\left(U\right)}_y=0,\left(x,y\right)\in R^2\\ U\left(x,y,0\right)=U_0\left(x,y\right).\end{array}$

其中

$U=\left[\begin{array}{c}\rho \\ \rho u\\ \rho v\\ \rho E\end{array}\right]$ , $F\left(U\right)=\left[\begin{array}{c}\rho U\\ \rho U^2+P\\ \rho UV\\ \rho UH\end{array}\right]$ , $G\left(U\right)=\left[\begin{array}{c}\rho V\\ \rho UV\\ \rho V^2+P\\ \rho VH\end{array}\right]$

ρ，P，H和E分别为密度，压强，单位体积的总焓和单位体积的总能。

选择解函数和试探函数空间如下：

$U_h=V_h=V_h^N\left\{p\in BV\cap L^1:{p|}_K\in P^N\left(K\right)\right\}$

K表示三角形剖分单元， $P^N\left(K\right)$ 表示单元K上次数不高于N的多项式。在三角形单元K上取基函数 $\left\{v_0,v_1,\cdots ,v_m\right\}$ ,则解函数 $U_h$ 可表示为：

$U_h={\displaystyle \underset{p=0}{\overset{m}{\sum }}{U}_K^{\left(p\right)}\left(t\right)\cdot v_p\left(x,y\right)}$ , $\left(x,y\right)\in K$

本文处理三阶精度的间断有限元方法，选取三角形单元K上的一组正交基函数如下 [6] ：

$v_0\left(x,y\right)=1$

$v_1\left(x,y\right)=\frac{x-x_0}{\sqrt{\left|{\Delta }_0\right|}}$

$v_2\left(x,y\right)=a_{21}\frac{x-x_0}{\sqrt{\left|{\Delta }_0\right|}}+\frac{y-y_0}{\sqrt{\left|{\Delta }_0\right|}}+a_{22}$

$v_4\left(x,y\right)=a_{41}\frac{{\left(x-x_0\right)}^2}{\left|{\Delta }_0\right|}+\frac{\left(x-x_0\right)\left(y-y_0\right)}{\left|{\Delta }_0\right|}+a_{42}\frac{x-x_0}{\sqrt{\left|{\Delta }_0\right|}}+a_{43}\frac{y-y_0}{\sqrt{\left|{\Delta }_0\right|}}+a_{44}$

$v_5\left(x,y\right)=a_{51}\frac{{\left(x-x_0\right)}^2}{\left|{\Delta }_0\right|}+a_{52}\frac{\left(x-x_0\right)\left(y-y_0\right)}{\left|{\Delta }_0\right|}+\frac{{\left(y-y_0\right)}^2}{\left|{\Delta }_0\right|}+a_{53}\frac{x-x_0}{\sqrt{\left|{\Delta }_0\right|}}+a_{54}\frac{y-y_0}{\sqrt{\left|{\Delta }_0\right|}}+a_{55}$

上式中正交基函数的系数 $a_{ij}$ 由待定系数法求得。由于基函数在单元K上是正交的，所以单元质量矩阵 $M_{ij}={\displaystyle {\sum }_{l=1}^6{\displaystyle \iint v_i{v}_j\text{d}x\text{d}y}}$ 为对角阵。计算过程中通量采用间断有限元方法离散，原双曲方程可以写成半离散形式如下：

$\frac{\text{d}}{\text{d}t}{U}^h=M^{-1}L\left(U^h,t\right)$

使用三步Runge-Kutta方法对上式进行迭代求解稳定解。

3. GPU加速

3.1. Cuda概述

CUDA是显卡制造商Nvidia公司2007年推出的并行计算平台和编程模型，开发人员可以使用c/c++语言来为CUDA架构编写程序。CUDA提供host-device的编程模式以及非常多的接口函数和科学计算库，通过执行大量的线程来达到加速效果。CUDA通过线程格来管理线程。线程格包含多个线程块，线程块包含多个线程。线程格和线程块可以是一维、二维或者三维结构。本文使用一维线程块结构，其线程类似于一个二维的像素数组，如图1所示。

线程索引号为：

其中threadIdx.x表示线程块中的线程编号，blockIdx.x表示线程格中线程块的编号，blockDim.x表示每个

Figure 1. Two-dimensional organazition of blocks and thread

图1. 线程块集合与线程集合的二维组织形式

线程块中的线程数。

3.2. 单GPU加速

单GPU加速时，由于CUDA支持大量的线程，因此选取CUDA线程数等于三角形网格个数。每个线程根据自己的全局索引号计算相应的网格单元即可。另外，还需要在GPU与CPU之间进行数据传输，所以程序的流程图如下：

上式中kernelFun即为核函数，将其循环迭代到流场收敛为止。

我们通过核函数依次求解时间步长、计算守恒量、处理边界条件、计算体积分守恒残差、计算通量、计算单元各条边上的积分残差、最后执行Runge-Kutta时间推进。

3.3. 多GPU加速

在多GPU下，首先将网格进行区域划分，划分后一个GPU对应计算一块区域。划分后可以不对网格进行处理直接计算；也可以先处理网格，将本地网格进行重排序，并重新生成邻接信息等。以图2的网格为例进行说明，我们将它沿中间划分成两个区域。

若直接计算，需要在每个计算节点上储存所有网格单元信息，即需要储存全部8个单元。计算时，只处理属于分配给相应节点的网格，第一个计算节点只计算单元编号为1、2、5、6的流场变量值。这样没有打乱网格单元的邻接关系，操作方便。但是，需要储存所有的三角形单元信息并传输到GPU上，浪

Figure 2. Triangle unstructured mesh

图2. 三角形非结构网格示例

费内存。GPU个数的增加并不能增大问题求解的规模，没有充分利用多个GPU计算的优势。同时，划分后每个区域的单元编号不再都是连续的，因此每个GPU上处理的单元在内存上也不是连续的。当warp请求访问位于不同地址的数据时，数据是非连续的，无法合并访问，大大影响访存的效率 [7] 。

若先将网格进行处理，在每个计算节点上，去除掉其他节点对应的网格，比如第一个计算节点上只需要储存编号为1、2、5、6的单元，所需的储存空间小了一倍。将分配到此节点的网格单元重新连续编号，即将5、6号单元重新编号为3、4，这样每个节点上的网格单元都是连续存储的，并重新生成单元其他信息，再进行计算。缺点是处理网格会稍微麻烦，但是节约了内存，还提高了访存效率，总体上，效率会远高于前一种方式。因此本文采用这种方式。

计算过程中，需要在网格边界处进行信息交换，采用MPI进行通信 [8] 。MPI是目前最流行的用于并行编程的消息传递接库标准。该模型底层是一组处理器，每个处理器有独立地址空间而且只能访问本地的指令和数据。各个处理器的不同进程间交互必须通过消息传递操作来实现。MPI定义了点到点通信、集合通信、并行I/O等操作。

在信息交换时，我们先将数据从GPU拷贝到CPU，再通过MPI在CPU之间交换，最后将数据拷贝回GPU。由于本文处理的是定常问题，因此不需要每一个时间步步都进行交换。同时MPI采用非阻塞通信，上一个时间步发送和接收消息，到下一个时间步才使用MPI_WaitAll来同步信息的交换。由于cuda核函数也是异步进行的，因此这样可以较好的覆盖CUDA的计算和MPI消息传递，提高并行效率。

4. 数值实验

求解naca0012翼型在来流马赫数0.4，攻角5˚下流场收敛时的密度。网格规模为37,360个三角形单元，如图3所示。

分别在CPU以及不同个数GPU上进行了计算，CPU使用Intel Xeon X5690 (主频3.47 Hz)，GPU使用Nvidia Tesla k20 GPGPU，含有2496个流处理器。GPU个数等于网格划分后的区域个数，两个GPU

对应的网格划分结果如图4所示，四个GPU对应划分后的四个区域如图5所示。计算耗时统计如表1所示。结果表明，GPU个数并不影响最终的流场解，流场解的密度等值线如图6所示。

结果显示，单个GPU对比单个CPU获得了33倍的加速比。采用不同GPU个数相对CPU加速比如图7所示。发现随着节点个数的增加，加速比逐渐增加，能够很好的处理大规模的计算。但是由于多个

Figure 3. Triangle unstructured mesh

图3. 三角形非结构网格

Figure 4. Mesh divided into two parts

图4. 网格划分为两个区域

Figure 5. Mesh divided into four parts

图5. 网格划分为四个区域

Table 1. System resulting data of standard experiment

表1. CPU及GPU上运行时间统计

Figure 6. Density contour map of solution

图6. 解的密度等值线图

Figure 7. Speedup ratio of different numbers of GPU relative CPU

图7. 不同GPU个数相对CPU加速比

节点计算会需要节点之间的协调和通信等耗时操作，因此加速比并没有和节点个数成正比关系。

5. 结论

本文基于CUDA采用多GPU对龙格库塔间断有限元方法进行了加速。分别测试了一个、两个和四个GPU，并取得了较好的加速效果。充分利用了多个GPU的特性，可以克服规模增大后单个GPU显存不足的情况。同时，本文针对定常问题，采用数据交换的下一个时间步同步并处理交换的信息，减少了通信的代价，增加了并行效率。