高缓钦,陈红全,张加乐,贾雪松

(1.南京航空航天大学 航空学院,南京 210016;2.非定常空气动力学与流动控制工信部重点实验室(南京航空航天大学),南京 210016)

间断Galerkin(Discontinuous Galerkin method,DG)有限元方法源于20世纪70年代求解中子输运方程[1]。后期流行的Runge-Kutta DG(RKDG)有限元方法[2-3],是结合求解非线性守恒律方程(组)提出的。该方法的主要特点是能够构造任意阶精度的格式(基于非结构网格等),但相较于传统的低阶精度格式,同一网格上控制方程的求解,因高阶近似而涉及的未知量更多,消耗的计算量也就更多,影响到实际复杂工程问题的推广应用[4-5]。因此,很有必要提高算法的计算效率。实现算法的并行化是可考虑的可行途径。

现代图形处理器(GPU)通常拥有数以百计或千计的计算核,浮点运算能力强劲,适合并行地处理数据密集且规模大的运算[6-8]。DG有限元方法具有很好的局部紧致性,构造的高阶格式不需要非常宽的模板(stencil),因此,十分适合算法的GPU并行化[9]。早在2009年,Klöckner等[9]就把Nodal DG有限元方法GPU化,求解了三维Maxwell方程。之后,Siebenborn等[10]又发展用于求解三维Euler方程。Karakus等[11]则拓展用于求解二维不可压Navier Stokes方程。与上述基于Nodal DG有限元方法不同,Xia等[12]基于Modal DG有限元方法,实现了三维Euler方程GPU并行加速。Fuhry等[13]则结合求解二维Euler方程,依据单元或单元边界构建所需的线程结构,开展了算法GPU并行化的深入研究,其并行化效率能与Nodal DG GPU算法求解线性问题[9]相当。Gao等[14]进一步结合Navier Stokes方程,发展用于求解二维层流流动。不难看出,上述工作大多关注的是数据存储、线程结构和核函数构建等具体的GPU化问题,尚未涉及网格拓扑存在的单元不规则分布(非结构等)对并行效率的影响问题。……