基于GPU并行算法的水动力数学模型建立及其效率分析

赵旭东，梁书秀＊，孙昭晨，刘忠波，2，韩松林，任喜峰

（1.大连理工大学 海岸和近海工程国家重点实验室，辽宁 大连 116024；2.大连海事大学 交通管理学院，辽宁 大连 116026）

0 引 言

海洋水动力模型的发展，一方面要求实现不同动力因素的耦合，如目前风－浪－潮耦合模型的需求，这将导致计算区域加大；另一方面需要对所关心区域不同尺度物理现象的精确模拟，这要求局部网格较小.这两方面的要求使采用非结构化网格模拟计算大范围海域的水动力模型成为一种趋势.但是网格密度和计算区域的增大将导致网格数量变得十分庞大，在没有集群机或在集群机计算条件受限情况下，单机计算时间过长，难以在较短时间内获得计算结果，在一定程度上限制了非结构化网格水动力模型在工程领域上的应用.

GPU（graphic processing unit，图形处理器）与CPU（central processing unit，中央处理器）的计算能力的差别主要源于二者的计算架构不同，一般的处理器由控制单元、逻辑单元和存储单元3个部分组成.在CPU 中，控制单元和存储单元占的份额较大，而GPU 中拥有较多的逻辑运算单元，因此相对于CPU，GPU 拥有更强的浮点运算能力.以本文使用的GTX460 显卡为例，其拥有336个流处理器，核心频率为700 MHz，显存容量为1GB，浮点运算能力达到0.961TFLOPS.

早期的GPU 仅仅用来完成图形处理工作，但随着技术的进步，对GPU 的可编程性增强，其能力已经超出了原先设计的功能，尤其是近年来，图形处理器硬件水平得到高速发展［1］，更多的开发者将其强大的并行计算能力应用到通用计算领域.目前GPU 已经应用于n－body模拟、医学CT图像重建、计算流体力学、分子动力学模拟、金融计算和油／气藏模拟等方面.2003年Li等［2］利用API环境实现了格子玻尔兹曼方法在GPU 上的计算.2011年Aaron把GPU 应用在利用大涡模拟解决基于结构化网格的湍流问题上，选用的CPU 是AMD Phenom 2.6GHz四核处理器，在Linux x64操作系统上运行，最多将计算速度提高了15 倍左右.2011年朱小松［3］实现了基于GPU 并行加速的粒子法解决液体晃荡问题，加速能力最多提高到26倍.

总的来说，目前基于GPU 的流体计算，多采用粒子法求解拉格朗日型的N－S方程，而对于求解欧拉型N－S方程，尤其是非结构网格下，利用GPU 解决流体动力学问题的研究很少.针对以上问题，本文采用CUDA 编程模型，基于GPU 架构并行算法建立非结构化网格水动力模型，并对所建模型的计算效率与集群机的性能进行对比分析.

CUDA平台使用带有关键词扩展的标准C开发语言，可以通过类C语言编写在显卡芯片上执行的程序.CUDA 编程模型中，将CPU 作为主机（Host），将GPU 作为协处理器或设备（Device）［4］，二者共同协作完成整个程序.其中，Host端负责执行串行部分的程序，Device端负责并行部分，Device端的程序又称为“kernel”.GPU 执行时的最小单位是thread，CUDA 架构可以拥有海量的thread，非常适合大规模数据并行计算［5］.

GPU 并行计算相对于传统的CPU 有以下几个优势：成本相对较低，一个普通计算机的价格，计算性能达到了TFLOPS，相当于一个高性能计算集群系统［6］；具有可扩展性，可以多个GPU 协同并行计算，支持高效线程，使得应用程序可以提供比现有硬件执行资源更高的并行度［7］；节点之间可以共享存储器数据，显存带宽是内存的10倍左右，可以高速传递数据；有众多的核心，完全可以使用胖节点式的并行算法，而CPU 的大型的并行，都是在集群层面的，而不是在单机上的；CUDA 这一类的GPU 并行编程语言，并行层面在线程级别，相对于CPU 要方便许多.

1 非结构化网格水动力模型

1.1 数学模型

从N－S方程出发，假设海水为不可压缩黏性流体，密度恒定，以连续性方程和动量方程作为二维水动力数值模型的控制方程：

底部剪切力

固边界条件

其中（u，v）为水平面流速，（τsx，τsy）为表面风应力，（τbx，τby）为底部剪切应力，f为科氏力系数，η为水面高程，M为曼宁系数，C为谢才系数.根据Smagorinsky涡黏参数法［8］，水平黏性项可以由下式近似表示：

1.2 空间离散方法

本文使用的模型选用非结构化三角形网格，采用有限体积法进行数值求解，在求解过程中将整个数值计算区域划分为相互不重叠的非结构化三角形网格，如图1所示.

图1 数值求解离散化非结构化三角形网格示意图Fig.1 Illustration of numerically solving the discretized unstructured triangular grid

三角形单元τ中的变量为线性分布，满足二阶空间精度，如式（11）～（15）所示：

1.3 时间离散方法

在本文的数学模型中对连续性方程（1）求面积积分，通过修正过的四阶R－K 时间积分进行求解运算，如式（16）～（18）所示：

2 基于GPU并行算法的实现及优化

目前CUDA 的应用程序大多为基于C 语言进行开发的，程序实现相对容易.而现在广泛使用的海洋模型如POM、ECOM、FVCOM 等均为使用FORTRAN 语言进行开发的，虽然可以在原有的模型中调用CUDA 计算核心，但是并行计算部分的设计，与重新建立模型的工作量相差不大，计算效率也略低于C程序，因此本文以C 语言为开发语言，在CUDA 平台下建立基于GPU 并行计算的水动力模型.

在基于CPU 并行计算水动力模型中，由于CPU 核心的数量远小于网格节点个数，传统的方法是利用美国密西根大学Karypis等编写的用于图的分区和稀疏矩阵排序的串行包——METIS库，将计算区域分为N个子区域，分配给N个CPU 进行计算，把子区域的初始流场信息、几何信息分别存储到对应的CPU 内存中，在每一个CPU中启动计算进程，由主进程调度各CPU 的计算［9］.在这一模式下，计算速度受到了多个方面因素的影响，包括CPU 的计算速度、内存带宽、在分区的边界处须保持CPU间数据交换过程等，因此并行计算的效率并不随CPU 核心的个数呈线性增加.

由于GPU 的并行计算架构与CPU 的并行计算架构有较大的区别，在并行化的实现过程中也有很大的不同.GPU 中的计算核心称为流处理器（SP），以本文中使用的GTX460 显卡为例，包含了330个流处理器.GPU 通过同时往每个流处理器上并发执行线程来实现加速计算.若像传统的并行方式一样利用METIS库将计算区域分块计算，一方面如果分块过多会导致计算过程过于复杂，另一方面并行程序会因为每个分块之间进行数据交换而导致并行计算效率的下降，因此在GPU 并行程序设计中尽可能多地分配线程，每个线程执行对一个网格或网格节点的计算.

2.1 GPU 并行算法的实现

在GPU 上实现并行算法需要以下几个步骤：

（1）执行GPU 并行计算前检测设备，根据不同的硬件设备，设置执行环境，并在CPU 上完成程序的前处理.

（2）在显卡上分配存储空间，用cudaMemcpy函数将在文件中读取的初始数据发送到显存中.

（3）将串行计算程序中需要每次对网格循环计算的程序段用对应的GPU 并行子程序实现，在每个并行子程序中设定在GPU 并行计算中的Grid和Block的维度.

（4）输出计算结果，用cudaMemcpy函数将显存中的数据返回到内存中，再输出到结果文件.

2.2 对GPU 并行算法的优化

在GPU 上实现并行计算，核心技术是算法的并行化和程序优化.在本文中按照计算模型的特点进行了适当的优化.

在GPU 并行程序中处理的对象是空间离散后的每个网格以及每个网格节点，这样就避免了由于分区不均衡引起的计算效率降低，同时还节省了分区边界处数据交换的时间.因此在原串行程序中的每次对网格数组或节点数组进行操作中，将其每次循环需要计算的内容对应地分配到GPU 的各个线程中，以海量的线程个数隐藏了数据访问延迟的时间.

由于在本文中使用的是非结构化网格，无法按照矩阵方式分配在每个线程中计算的网格，故在GPU 并行程序的线程设计中将Block设置为一维，并且每个Block中包含64个线程，即两个线程束.则整体的Block的个数为N／64＋1（其中N为计算单元个数），这样就将GPU 划分为（N／64＋1）×64个线程，这种线程分配方式最大限度地降低了分支流程和空载线程的个数，实现了GPU 核心的最高利用效率.根据式（1）～（16），每个GPU 的线程计算一个网格或网格节点的数据，将所需要的数据和需要计算的结果作为kernel函数的参数，并在所有节点的数据计算完成后，将计算结果返回到显存的数组中.

为了减少显存和内存之间数据传输速度的限制，在整个计算过程中，完成初始化后将变量传递到显存中，在迭代计算过程中，所有数据均在显存中进行交换，直至需要将中间结果输出，再将显存中的数据返回到内存的数组中，这样减少了由于数据交换引起的时间损耗.

3 模型验证

数值实验的模型水深D＝10 m，长L＝3.7km，宽B＝1.5km，整个流场的右侧为开边界，有作用振幅为0.5m 的M2分潮.为了避免受到边界影响，在该数值实验结果对比中，验证点选择整个流场的中间点，距离左侧边界1.85 km［10］.

根据Ippen 的理论［11］，对于此种水域情况，有如下解析解：

其中A为振幅，ω为角频率，x为该点距离左侧边界 的 距 离，h为 水 深，T为 潮 周 期，η为 潮 波 面 升高，u为纵向流速.

在数值模型中将计算区域剖分成1 620个三角形网格，验证点潮位和纵向流速计算结果历时曲线如图2所示.从图中可以看出，数值计算结果与解析解吻合度较好.

图2 验证点潮位及流速历时曲线计算结果对比Fig.2 The contrast of elevation and current speed duration curve result in the verified points

4 GPU 算法的并行效率分析

为了更好地说明GPU 在水动力模型并行计算过程中起到的作用，数值模型实验采用在上一章模型验证中的计算水域，并设计了4种不同网格密度的工况算例，分别使用CPU 单线程（Intel CoreTMi7 930 2.8GHz）、集群机（使用30核心、24核心、20核心，最大计算节点数为12，CPU 型号为Intel Xeon E5430，主频2.66 GHz），以及GTX460显卡（330 个流处理器，核心频率700 MHz）进行计算，并对计算耗用时间进行对比分析.在集群机计算过程中发现，在工况1～3中，网格数不多时，总的计算核心数保持不变，集群机节点数以及每个节点使用核心数变化对计算时间影响不大，计算所需时间基本不变；在工况4中，网格数超过20万，受集群机节点数以及每个节点使用核心数影响明显，不同节点的分配方式得到的加速比R如图3所示，其中加速比为单线程计算时间与集群机计算使用时间的比值，横坐标为节点数×每个节点使用核心数.

由图3可以看出，在核心总数相同的条件下，使用的节点数越多，计算效率越高，相应的加速比就越高.因此集群机的计算效率受计算核心总数，以及使用节点个数的影响比较明显.

采用集群机的最优化的节点分配方案，并与GPU（Nvidia GTX460）并行计算的加速能力进行对比如表1和图4所示.

图3 集群机在不同节点分配条件下计算工况4得到的加速比Fig.3 The speedup ratio in different distributions of compute nodes when using cluster to compute Case 4

表1 集群机和GPU 不同网格数条件下所用时间Tab.1 The time of cluster and GPU cost in different number of grids

图4 不同网格数条件下并行加速能力对比Fig.4 The contrast of speedup capacity at different number of grids

表1中结果为分别使用单线程CPU 串行算法、3种不同计算核心数的MPI并行算法，以及GPU 并行算法，计算4种不同网格数下算例所使用时间.图4为表1中各并行算法所对应的加速比.在本文中最多只使用了12个计算节点，由图4可以看出，随着网格个数的增加，GPU 并行计算的加速比能够实现相对稳定的提高.计算网格数在工况1～3 时，加速比提高的速度低于集群机，网格数在工况3～4时，由于节点使用个数的限制，使用30 核和20 核的加速比略有下降，GPU 的加速比仍稳定地提高，其加速比提高的速度基本与24核集群机的平行.由加速比随网格数变化曲线的趋势可以看出，若网格数继续增加，集群机的加速能力会受到节点分配方案和等待其他分区计算完成时间的影响而有所降低，而GPU加速能力仍会相对稳定地提高.

5 结 论

（1）对不同计算核心分配方案的MPI并行计算与单CPU 串行计算性能进行对比分析，得到结论：当网格节点数在工况1～3（即网格数较少）时，不同的节点分配方案对集群机的加速能力影响不大；在网格数量较多时，集群机加速能力受计算节点分配方式影响显著.

（2）通过对不同计算核心数的MPI并行算法以及GPU 并行算法的加速比的对比得出：使用GPU 并行算法，能够在保证较高计算精度的基础上，大幅提高方程求解的速度，尤其是随着网格数的增加，基于GPU 求解方程的计算加速比得到显著的提高.当网格个数增加到20万个以上时，GPU 的计算能力已经接近使用20核集群机的计算能力，加速比提高了10倍以上，并且加速比的变化率已经超过20核和30核的集群机，并行加速效果十分显著.若考虑其他方面的计算成本，如集群机节点的利用率、能耗等方面，使用GPU 进行并行计算的优势更加明显.由此可见，GPU 在大范围海域的水动力模型的数值计算方面具有较好的发展前景.

［1］ NVIDIA.NVIDIA CUDA C programming guide［EB／OL］.［2013－08－19］.https：／／www.nvidia.com.

［2］ LI Wei，WEI Xiao－ming，Kaufman A.Implementing lattice Boltzmann computation on graphics hardware ［J］.Visual Computer，2003，19（7－8）：444－456.

［3］ 朱小松.粒子法的并行加速及在液体晃荡研究中的应用［D］.大连：大连理工大学，2011.ZHU Xiao－song.Parallel acceleration of particle method and its application on the study of liquid sloshing ［D］.Dalian：Dalian University of Technology，2011.（in Chinese）

［4］ 张 舒，褚艳丽.GPU 高性能运算之CUDA［M］.北京：中国水利水电出版社，2009.ZHANG Shu，CHU Yan－li.GPU High Performance Computing：CUDA［M］.Beijing：China Water Power Press，2009.（in Chinese）

［5］ 王英俊，王启富，王 钢，等.CUDA 架构下的三维弹性静力学边界元并行计算［J］.计算机辅助设计与图形学学报，2012，24（1）：112－119.WANG Ying－jun，WANG Qi－fu，WANG Gang，etal.CUDA based parallel computation of BEM for 3Delastostatics problems［J］.Journal of Computer－Aided Design ＆Computer Graphics，2012，24（1）：112－119.（in Chinese）

［6］ 多相复杂系统国家重点实验室.基于GPU 的多尺度离散模拟并行计算［M］.北京：科学出版社，2009.State Key Laboratory of Multiphase Complex Systems.Multi－scale Discrete Simulation Parallel Computing Based on GPU ［M］.Beijing：Science Press，2009.（in Chinese）

［7］ Kirk D B，Hwu Wen－mei W.Programming Massively Parallel Processors［M］.Beijing：Tsinghua University Press，2010.

［8］ 陈长胜.海洋生态系统动力学与模型［M］.北京：高等教育出版社，2004.CHEN Chang－sheng.Ocean Ecosystem Dynamics and the Model ［M］.Beijing：Higher Education Press，2004.（in Chinese）

［9］ 王 敏，王 明，杨 明.小浪底水库三维数学模型并行计算研究［J］.人民黄河，2012，34（5）：25－27.WANG Min，WANG Ming，YANG Ming.Parallel computing research of Xiaolangdi Reservoir threedimensional mathematical model［J］.Yellow River，2012，34（5）：25－27.（in Chinese）

［10］ 张瑞瑾.三维潮流数值模型及其应用［D］.大连：大连理工大学，2002.ZHANG Rui－jin.Establish and the application of a three dimensional numerical tidal model ［D］.Dalian：Dalian University of Technology，2002.（in Chinese）

［11］ Ippen A T.Estuary and Coastline Hydrodynamics［M］.New York：McGraw－Hill Book Co，1966.