CN110516316A - 一种间断伽辽金法求解欧拉方程的gpu加速方法 - Google Patents

Abstract

本发明属于计算流体力学、高性能计算领域，涉及一种间断伽辽金方法(DGM)的GPU并行加速技术，具体为一种间断伽辽金法求解欧拉方程的GPU加速方法。本发明采用四面体网格对求解区域进行剖分，以基函数、高斯积分、数值通量为基础，GPU为主要计算硬件，CUDA为编程模型建立间断伽辽金方法GPU并行框架。通过CUDA并行框架实现了GPU多线程的管理，通过设计的数据结构和线程访问方式来实现高效的内存访问。在解决面积分需要单元间数据交换而不独立的问题时，采用了按照面网格并行、每个面网格的计算线程处理两个单元的方式，既避开了单元不独立的问题，实现了大规模并行，还减少了计算量。

Description

一种间断伽辽金法求解欧拉方程的GPU加速方法

技术领域

本发明属于计算流体力学、高性能计算领域，涉及一种间断伽辽金方法(DGM)的GPU并行加速技术，具体为一种间断伽辽金法求解欧拉方程的GPU加速方法。

背景技术

欧拉方程是流体力学中描述无粘流体的方程组，其形式如下：

U_t+▽·F＝0 (1)

其中U代表守恒量、U_t代表守恒量对时间t的偏导数，F代表守恒通量，▽·F代表守恒通量的散度，且在三维情况下，有

其中ρ为气体密度，u、v、w为气体的三个速度分量，e为完全气体的单位体积总能量，p为气体压强。

对于上述欧拉方程的数值求解，通常采用以下几种方法：有限差分法、有限体积法、有限元法。其中有限差分法需要采用结构网格，且计算量小，常用于处理结构网格划分的简单几何区域上的求解，而对于复杂几何区域的求解则相对困难。有限体积法可以求解结构或非结构网格，因此可以处理复杂集合区域，应用范围相对较广，但其难以构造高阶格式(需要扩展模板)，且构造的高精度格式要么求解复杂、要么不够紧致。而间断伽辽金方法(DGM)则结合了有限元和有限体积方法，能够处理任意网格和复杂几何区域，而且DGM可以通过简单地增加单元内的解多项式的次数进而增加单元自由度(DOFs)来获得更高的空间精度，是一种高精度的流场求解方法。

间断伽辽金法的数值求解过程如下：将欧拉方程的两端乘以试探函数Φ并在体积Ω上积分，经积分变换可以得到伽辽金弱形式

将积分区域Ω划分为四面体网格，并取其中的一个单元Ω_k来考察上述方程，为了保证单元之间场的连续性，上式的面积分项的被积函数需要使用“数值通量”F*来代替，数值通量由单元两侧的数值计算得到，于是可以将上式重写为

其中u⁺,u^-分别代表积分面两侧的数值。若将守恒量u用基函数φ_i展开，并且试探函数也为φ_i，可得到

其中左端项仅与基函数有关，使用正交的基函数可以得到一个对角的质量矩阵，且仅与本单元相关。右端第一项为体积分项，仅与本单元项相关。以上两项是DG中最直接具有并行性的部分。右端第二项为面积分项，采用数值通量之后，这一项与积分面两侧的场都有关。

虽然DGM有更高的精度，也存在的一个计算量较大的问题，往往需要比其他方法更长的计算时间。但单元DOFs在单元间相对独立的特点使得该方法具有天然的并行性，非常适合于大规模并行计算，利用这一点可以弥补其计算量大的问题。

在并行计算技术方面，当前主要有多核中央处理器(CPU)并行、图形处理器(GPU)与CPU的异构并行两类。其中，多核CPU并行出现最早，当前超算采用最多的架构就是多核CPU架构，通过增加核心数来增加并行任务的数量。这一方法的缺点也很明显，当前单个的CPU线程数只有最多64个，而再需求更多核心则需要增加CPU个数，这会大幅度提高成本，且对于小型工作站而言CPU数目限制比较大，难以实现较大规模的计算。对于GPU并行，以英伟达TITANV为例，其包含80个流处理器簇(SM)，每个流处理器簇又包含64个流处理器(SP)，在英伟达CUDA编程模型下，可以同时处理海量的线程以实现大规模并行计算。当前的英伟达GPU设备的每秒浮点运算次数已达到15.7TFLOPS，是同时期英特尔CPU的10倍，使用一个GPU工作站将能够取代十台CPU工作站，且成本更低。GPU-CPU异构架构则是一种高效的并行计算架构，以当前世界超算天梯榜首Summit为代表的一系列超算就是采用使用这样的架构。异构计算平台的CPU和GPU的数量比较灵活，最简单的应用仅需要一个CPU与一个GPU，即可实现大规模的并行计算。

如果能将GPU并行计算应用于DGM可以很好地弥补DGM计算量偏大的问题，但其中存在若干待解决的问题：在技术上，传统的CPU运行的DGM通常依赖于现成且优化过的线性代数及通信原语，而这些在GPU上都是不可用的；在实施上，虽然半离散格式(5)是基于一个单元展开的，但还是通过右端第二项与其它相邻单元进行了数据交换，也就是这一部分并非是单元独立的，需要通过一定的实施手段来保证其能够并行且效率足够高。

发明内容

针对上述存在问题或不足，为解决现有欧拉方程的数值求解DGM结合多核CPU的方法成本相对较高，并且相对效率低下的问题，本发明提供了一种间断伽辽金法求解欧拉方程的GPU加速方法，基于间断伽辽金法的特点和GPU并行计算的优势，将GPU并行计算应用于间断伽辽金法求解欧拉方程,能最大程度够发挥间断伽辽金法的优势，加速计算流体力学高精度格式的计算。

该间断伽辽金法求解欧拉方程的GPU加速方法的具体技术方案包括以下几个步骤：

步骤1、读取计算网格的信息，并将面网格按照边界类型排序；

步骤2、在CPU端进行几何参数的预计算，并将结果拷贝到GPU显存；预计算并存储几何参数，是因为几何参数在计算全程不发生改变，只需计算一次即可，且几何参数的计算需要访问非连续的内存，无法实现对齐、合并的访问，若在GPU端的DGM计算过程中频繁计算，将大幅降低计算效率。

步骤3、在GPU端完成流场的初始化，并且将时间步初始化为0；初始化时分配与单元数量一致的线程数，每个线程负责初始化单元内的所有场。

步骤4、判断计算时间步是否达到预定的终止时间步，若达到，结束计算，否则继续；

步骤5、在GPU端进入时间推进，具体为：依次启动面上场计算核函数、面积分核函数、体积分核函数、以及时间更新核函数。

步骤6、再次判断当前计算时间步是否达到预定的后处理时间步，若达到，同步GPU线程，并将计算结果拷贝到CPU端，将当前计算得到的流场数据输出；否则不进行任何操作。

步骤7、时间步自增1，转至步骤4。

上述步骤中GPU显存不存储原始的网格信息(如网格结点坐标)而只存储：计算直接需要的且由原始数据计算得到的数据(法向量、体积和面积)，网格拓扑关系(如一个单元的每一个面的全局编号等)，以及流场数据(场、右端项)，并使用数组结构体的方式组织并存储在GPU全局内存中，便于管理与使用，而CPU端仅需在预计算时保留网格原始信息和流场数据，在预计算完成后即可释放，仅需要保留计算的场数据。对于高斯积分点的基函数值、高斯积分权重、常数质量矩阵三类所有线程都会使用到的常数，存储于GPU的常量内存中。

进一步的，所述步骤1中对面网格进行的排序时，以网格编号为关键字，采用桶排序算法进行排序。

进一步的，所述步骤3中使用由文件中导入的流场进行初始化。

进一步的，所述步骤5中，对于所有计算核函数的线程分配，为了保证对齐与合并访问来提高内存访问效率，需要将4个单元或者面划分到一个线程块内，并在同一个线程束内，且对内存中数组的访问正好实现128字节的对齐以及合并访问，线程结构及对内存访问方式如图2所示。

进一步的，所述步骤5中面积分核函数，即式(5)中右端第二项的数值通量，采用计算过程中无需判断的Roe通量格式，以减少线程束分化。

进一步的，所述步骤5中，面积分核函数与体积分核函数的计算相互独立，使用CUDA流，将面上场计算核函数、面积分核函数发布到第一个流，体积分核函数发布到第二个流，再将时间更新核函数发布到第一个流，以实现体积分核函数和面积分核函数的并发，减少计算时间。

本发明采用四面体网格对求解区域进行剖分，以基函数、高斯积分、数值通量为基础，GPU为主要计算硬件，CUDA为编程模型建立间断伽辽金方法GPU并行框架。通过CUDA并行框架实现了GPU多线程的管理，通过设计的数据结构和线程访问方式来实现高效的内存访问。在解决面积分需要单元间数据交换而不独立的问题时，采用了按照面网格并行、每个面网格的计算线程处理两个单元的方式，既避开了单元不独立的问题，实现了大规模并行，还减少了计算量。

综上所述，本发明通过将DGM与GPU结合运用，解决了DGM计算量偏大的问题，相比现有欧拉方程的数值求解采用的DGM结合多核CPU的方法，成本更低且效率更高的问题。

附图说明

图1为本发明的整体流程图；

图2为并行数据的组织方式；

图3为线程结构以及线程访问全局内存方式的示意图；

图4为实施例的GPU计算结果与CPU版本的DGM计算结果的对比图；

图5为实施例实测的时间对比以及加速比。

具体实施方式

下面结合附图和实施例对本发明进行详细说明。

第一部分：采用一阶拉格朗日基作为基函数离散式(5)，并使用显式龙格-库塔(R-K)时间格式

进行时间推进，以数组结构体(SOA)作为并行数据结构，并且对内存布局、线程结构进行设计，具体包括以下几个步骤：

S101、GPU显存不存储原始的网格信息(如网格结点坐标)而只存储计算直接需要且由原始数据计算得到的数据(如场、法向量、体积、面积)以及网格拓扑关系(如一个单元的每一个面的全局编号等)。使用SOA的方式组织存储在GPU显存中的计算数据，并为其分配内存，其中，体网格的数据在结构体ElementArray中，包含：五个守恒量数组U[5]，存储式(5)的五个右端项数组RHS[5]，存储R-K第一步计算值的数组RK[5]，以及体积与单元4个面的面法向量volume，nx[4]，ny[4]，nz[4]；对应于面网格的数据在结构体FaceArray中，包含：面左右两侧的守恒量数组U_L[5]和U_R[5]，存储式(5)的面积分右端项的数组RHS_L[5]和RHS_R[5]，以及面积与指向其左单元的单位法向量area，nx，ny，nz；对于每个非几何参量，在内存中的布局方式为每单元以四个DOFs相邻排列为一组，各单元的组相邻排列，如图1中的U、RHS、RK2等数组，因此每个数组的大小为4*K，其中K为单元或者面网格的总数；对于每个几何参量，一个单元仅保留有一个数据，为了后续的对齐与合并访问，需要在每个数据后进行填充，如图1中的nx、ny、nz等数组，每个数组的大小为4*K，但仅在数组索引为4的倍数的位置存储。这些分别对应于单元和面网格的数据又被分别装在结构体Element和FaceArray中。

S102、计算中使用到的单元高斯积分点坐标，我们直接带入基函数得出一个大小为4*4的数组，代表每个积分点处每个基函数的值，记为tetra_basis_value；对于面高斯积分点坐标，我们同样带图基函数得出一个4*4*4的数组，代表每个面上每个积分点处每个基函数的值，记为face_basis_value。四面体单元和三角面高斯积分点的权重我们分别记作tetra_gauss_weight、face_gauss_weight。式(5)左侧形成的质量矩阵，在拉格朗日基函数下是对角的，因此仅存储其对角元，记为mass_diag。以上涉及到的常数数组将被所有单元共用，因此被布置于GPU的常量内存中。

S103、当前CUDA并行模型中，线程是以线程束为最小单位发布并执行的，且每线程束包含warpsize＝32个线程。本发明采用每线程计算一个DOF的方式，即一个单元或者面需要由4个线程来计算。为了保证对齐与合并访问来提高内存访问效率，需要将4个单元划分到一个线程块内，以保证一个线程块内的所有线程分配到一个线程束内并同时运行，且对内存中数组的访问正好实现128字节的对齐以及合并访问，线程结构及对内存访问如图2所示。

第二部分：包括网格导入及网格数据处理、数据预计算、DGM计算、后处理的整体框架。包括以下几个步骤

S201、将有限元网格的体网格和面网格数据读入第一部分所述的数组中，并以面网格的边界条件类型为依据，使用桶排序算法进行排序，使得具有相同边界类型的面网格相邻，目的是减少不同边界类型的面网格分配到同一个线程束而导致线程束分化的可能性；

S202、如步骤S101中所述GPU计算需要的直接数据，在计算全程中都不再改变，因此在网格信息读入之后，计算出体网格和面网格中的体积、面积和法向量，并拷贝到显存中这里的预计算可以选择由CPU还是GPU完成，若由CPU完成，需在计算完成后将预计算结果拷贝到GPU端，若由GPU完成，则需在预计算前将网格的结点信息和拓扑关系拷贝到GPU端用于预计算，在预计算完成后再将网格结点和拓扑关系数据占用的显存空间释放；之后，使用GPU核函数完成场的初始化；

S203、在主机端判断时间步n是否达到仿真步数上限，若是，结束计算，释放所有已分配的GPU显存空间，否则继续；

S204、DGM的GPU计算划分为体积分计算核函数VolumeKernel、面上场计算核函数FaceFieldKernel、面积分计算核函数FaceKernel、R-K第一步更新核函数RK1Kernel、R-K第二步更新核函数RK2Kernel，VolumeKernel负责计算式(5)中的右端第一项，FaceFieldKernel负责计算式(5)右端第二项中的u⁺，u^-，FaceKernel负责计算式(5)右端第二项，两个R-K更新核函数分别负责组合面积分与体积分，并计算式(6)中的两步时间迭代；

S205、每个时间步计算完成后判断时间步是否达到预先设定好的需要输出后处理文件的时间步，若达到，进行全局的线程同步，然后将GPU端计算的场拷贝到CPU端，并使用CGNS(计算流体力学通用符号系统)格式对场进行输出，并转S204，未达到则不执行同步，直接转S204；

第三部分：包括DGM计算的体积分、面积分和时间更新核函数三个核心计算函数的设计

S301、体积分核函数，采用每线程处理一个DOF、四个线程处理一个单元的方式分配线程。在内存的使用上，为单元的每个DOF、右端项、基函数的梯度、分配共享内存，并将DOF从全局内存读入，将法向量从全局内存读入、节点处基函数的值从常量内存读入，然后计算出梯度值。然后由DOF和常量内存中的高斯积分点计算出积分点处守恒量的值，进而求出体积分值存储在共享内存中，在计算完成之后，同步并转写到全局内存中。

S302、为了将容易产生线程束分化的部分独立出来，本发明中将面上场计算的核函数和面积分计算的核函数分开。面上场计算核函数采用每线程处理一个积分点，四个线程处理一个面的方式分配线程。在内存的使用上，将每个面的每个线程会公共使用的单元场的数据以及法向量存储在共享内存中，使用的高斯积分点坐标从常量内存中读取，法向量从全局内存中读取。计算时，先计算出所有单元计算方法都一致的左侧积分点场值，再根据边界条件计算出右侧积分点场值。这里由于事先为面按边界条件排过序，可以最大程度上避免线程束分化。计算的场值直接存储在全局内存中，供面积分核函数使用。

S303、面积分核函数，采用每线程处理一个DOF、四个线程处理一个面的方式分配线程。在内存的使用上，因为此部分计算有大量数据可以跨线程使用，因此将积分点处的值、数值通量的值、法向量、面积分数值、中间变量全都存放在共享内存中，而面上高斯积分点处基函数的值和高斯积分重权从常量内存中读取。计算时，先计算出所有单元的左侧的面积分值，再判断该面是否是内部面，若是，计算右侧的面积分值。这里的判断是可能导致此核函数线程分化的一个点，也是通过面按边界条件排序来减少线程束分化。在计算完成后，同步并将面积分值写回全局内存。

S304、时间更新核函数，在采用(6)所示的R-K格式时，需要两个时间更新核函数，这些核函数与体积分核函数一样采用每线程处理一个DOF，四个线程处理一个单元的方式分配线程。此部分仅需要很少的中间变量，因此不需要使用到共享内存，每个单元从其每个面上获取左侧和右侧面积分，然后判断此单元是该面的左侧还是右侧单元，再乘以一个权值后累加到本单元的体积分项上。例如，如果当前单元是一个面的右侧单元，那么这个面对此单元贡献的面积分就是：0*左侧面积分+1*右侧面积分，如果是左侧单元，那么这个面对此单元贡献的面积分就是：1*左侧面积分+0*右侧面积分。组合完所有面的积分到体积分上后，从常量内存读取质量矩阵对角元，并计算出更新之后的DOF，直接写回全局内存。

图4展示了采用本说明的方法，计算小球绕流问题的结果对比，计算条件为0.5马赫。采用本发明的GPU版本与CPU版本完全一致。图5展示了计算时间对比和加速比对比,其中的测试环境为：IntelXeonCPUE5-2697v4@2.3GHz四核，NVIDIAGRIDP40-1Q GPU，CPU版本使用openMP四线程并行，计算均使用双精度。

Claims (6)

1.一种间断伽辽金法求解欧拉方程的GPU加速方法，包括以下几个步骤：

步骤1、读取计算网格的信息，并将面网格按照边界类型排序；

步骤2、在CPU端进行几何参数的预计算，并将结果拷贝到GPU显存；

步骤3、在GPU端完成流场的初始化，并且将时间步初始化为0；初始化时分配与单元数量一致的线程数，每个线程负责初始化单元内的所有场；

步骤4、判断计算时间步是否达到预定的终止时间步，若达到，结束计算，否则继续；

步骤5、在GPU端进入时间推进，具体为：依次启动面上场计算核函数、面积分核函数、体积分核函数、以及时间更新核函数。

步骤6、再次判断当前计算时间步是否达到预定的后处理时间步，若达到，同步GPU线程，并将计算结果拷贝到CPU端，将当前计算得到的流场数据输出；否则不进行任何操作。

步骤7、时间步自增1，转至步骤4。

上述步骤中GPU显存不存储原始的网格信息而只存储：计算直接需要的且由原始数据计算得到的数据、网格拓扑关系以及流场数据，使用数组结构体的方式组织并存储在GPU全局内存中，而CPU端仅在预计算时保留网格原始信息和流场数据，在预计算完成后释放，仅保留计算的场数据；对于高斯积分点的基函数值、高斯积分权重、常数质量矩阵三类所有线程都会使用到的常数，存储于GPU的常量内存中。

所述计算直接需要的且由原始数据计算得到的数据为法向量、体积和面积。

2.如权利要求1所述间断伽辽金法求解欧拉方程的GPU加速方法，其特征在于：所述步骤1中对面网格进行的排序时，以网格编号为关键字，采用桶排序算法进行排序。

3.如权利要求1所述间断伽辽金法求解欧拉方程的GPU加速方法，其特征在于：所述步骤3中使用由文件中导入的流场进行初始化。

4.如权利要求1所述间断伽辽金法求解欧拉方程的GPU加速方法，其特征在于：所述步骤5中对于所有计算核函数的线程分配，将4个单元或者面划分到一个线程块内，并在同一个线程束内，且对内存中数组的访问正好实现128字节的对齐以及合并访问。

5.如权利要求1所述间断伽辽金法求解欧拉方程的GPU加速方法，其特征在于：所述步骤5中面积分核函数，即间断伽辽金半离散形式中右端第二项的数值通量，采用Roe通量格式。

6.如权利要求1所述间断伽辽金法求解欧拉方程的GPU加速方法，其特征在于：所述步骤5中，面积分核函数与体积分核函数的计算相互独立，使用CUDA流，将面上场计算核函数、面积分核函数发布到第一个流，体积分核函数发布到第二个流，再将时间更新核函数发布到第一个流。