近年来，随着计算机性能的不断提升，海洋水动力数值模型得到了快速发展，目前应用较为广泛的海洋数值计算模型主要有POM^[1]，ECOMSED^[2]，FVCOM^[3]，MIKE^[4]等。纵观主要的海洋数值模型可以发现网格非结构化及算法并行化成为近年来模型发展的趋势，例如FVCOM模型采用了非结构网格和VC-CC(vertex-centered and cell-centered)混合格式的有限体积法，通过体积通量的积分方法，求解控制方程组，对于具有复杂地形和岸线的区域，保证了质量、动量、盐度和热量的守恒性，并且该模型基于MPI并行编程，可实现并行计算，运行效率高；丹麦DHI公司开发的MIKE系列软件基于CC(cell-centered)格式有限体积法和ADI(alternative direction iterative)时间积分方法建立了二维和三维的潮流水动力模型，其水动力模型可以实现GPU、OpenMP(Open Multi-Processing)和MPI(Message Passing Interface)3种方式并行计算。相对于结构化网格而言，非结构化网格可以任意加密，对陆地边界拟合更精确，物质量守恒性更好，较易实现对模拟区域不同尺度物理现象的精确计算。并行算法相较于串行算法其计算时间可以大大缩短，适于模拟大范围，高密度网格区域的水动力。

随着图形处理器(GPU)性能的高速发展^[5]，同时伴随CUDA、OPENCL(Open Computing Language)等GPU架构的并行计算语言支持的扩展及日趋成熟，例如NVIDIA公司在CUDA C的编程语言基础上增加了对Fortran，Python，Matlab等的支持，使得并行层面的线程级别、基于CPU运算的程序代码的GPU加速改写，都变得方便易行。目前，GPU的加速能力已被广泛的应用于科学计算领域。如Michalakes^[6]等人利用GPU实现了对气象模型的加速，Solano-Quinde^[7]等人提出了基于GPU的非结构网格并行计算方法，Bleichrodt^[8]等实现了静压海洋模型的GPU加速，DHI MIKE ZERO也实现了GPU并行算法在水动力模型中的应用。由此可见，基于GPU的并行算法可以有效加快海洋模型计算速度，高效率的完成对高分辨率海洋环境水动力的数值模拟，有着十分广阔的应用前景。

本文基于CUDA C语言的GPU加速算法建立了非结构(三角形)网格有限体积法海洋模型(GPU-OM)。GPU-OM在时间积分上采用改进的4步RK积分方法；空间离散采用二阶迎风格式。在此基础上，利用大长山岛海域水文调查实测资料，将GPU-OM用于对长山群岛及附近海域二维海洋环境水动力场数值模拟，验证并分析了该海域的水位，流速的分布特征。同时，与成熟商业软件MIKE21HD模型的计算效率对比。

1   材料与方法

1.1   非结构网格水动力模型控制方程

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

式中:ξ为水位；(u, v)为水平方向流速；D为水深；Ω为控制体面积；Ê和P为蒸发和降雨的通量速率；Q_b为入流的通量速率；(τ_sx, τ_sy)和(τ_bx, τ_by)分别为表面和底部的剪切应力。

当水表面风速|W|∈(11, 25) 时，自由表面剪切力(τ_sx, τ_sy)为^[9]：

底部剪切力(τ_bx, τ_by)为^[4]：

固边界条件：

式中:(u, v)为水平面流速；(τ_sx, τ_sy)为表面风应力；(τ_bx, τ_by)为底部剪切应力；f为科氏力系数；ξ为水面高程；M为曼宁系数；C为谢才系数。根据Smagorinsky涡黏参数法，水平粘性项可以由下式近似表示：

1.2   空间离散方法

本文数值模型网格采用非结构化三角形网格，采用有限体积法进行数值求解，在求解过程中将整个数值计算区域划分为相互不重叠的非结构三角形网格。水位信息存储在网格的节点上，流速信息存储在网格单元的形心(●)上。

在求解连续性方程时，采用VC格式有限体积方法，其控制体由包围控制体的各个网格的形心及其公共边中点连接组成(Control_1阴影部分)；在求解动量方程时，采用CC格式有限体积法，其控制体为网格单元(Control_2阴影部分)如图 2所示^[10]。

图  2  非结构网格空间重构

Fig.  2  Illustration of the spatial reconstruction of the unstructured triangular grid

下载: 全尺寸图片 幻灯片

以连续性方程为例，在控制单元内，对式1进行积分：

再利用格林公式，将对面积积分转换为沿控制单元边界的线积分：

通过同样的方法可以得到空间离散后的动量方程：

在一阶精度的水动力模型，在各个单元内的水平流速分量为常数，容易在计算过程中出现较大的耗散，导致数值格式不稳定或计算结果不符合实际情况。而在二阶精度的数值模型中，水平流速分量变量在网格单元内为线性分布，通过空间重构方法可以有效提高数值格式的稳定性。目前大多数非结构网格水动力模型中采用了MUSCL格式或最小二乘法进行空间重构。其中，最小二乘法实现方法相对简单高效，因此本文选用该方法网格单元进行重构。

在图 1中，Ⅰ、Ⅱ和Ⅲ为与计算单元τ相邻的网格单元，在三角形单元τ中的水平方向流速(u, v)为线性分布，如下式所示：

图  1  非结构化网格

Fig.  1  The sketch map of unstructured grid

下载: 全尺寸图片 幻灯片

式中:(a, b)即为控制单元内水平流速沿(x, y)方向的梯度值。假设上式在控制单元及其相邻单元的范围内处处成立，根据控制单元及其相邻单元的局部坐标和水平流速值，可以得到方程组：

在上述方程组中含有三个方程，但是未知数只有a和b，该方程组为超定常方程组。根据最小二乘法，该单元及其相邻单元水平流速组成样本的均方差为：

为了使利用式(15) 计算的相邻网格单元水平流速分量与其原有的结果误差最小，则式(17) 取最小值，根据极值原理，可以得出：

根据上式，水平流速沿(x, y)方向梯度则可以通过控制单元及其周围相邻单元中心处的流速和局部坐标计算得出：

在相邻单元公共边中点处的梯度值为其两侧梯度的平均值：

1.3   时间离散方法

在目前主流的水动力数值模型中，高阶时间差分格式一般采用蛙跳格式、交替方向隐式法、预测—校正格式、经典的四阶Runge-Kutta法等。其中经典的四阶RK积分方法稳定区域最大，具有较好的稳定性，但是由于在每一步计算过程中都需要对中间通量进行存储，不利于显存容量不可扩展的GPU设备。为此，本文选用了Dick^[11]提出的改进的四步二阶RK积分格式。连续性方程和动量方程的时间离散格式如下：

上式中修正的四步二阶的RK方法与经典的四阶RK方法相比其精度虽然只有二阶，但是两者的稳定区域相同，计算精度也优于交替隐式法、预测-矫正格式的和2步经典RK积分格式；同时对存储空间的要求要小于4步CLRK法^[12]。

1.4   基于GPU并行加速算法实现及优化方法

1.4.1   基于GPU的海洋模型的实现

目前的大部分水动力模型都是由FORTRAN语言编译生成的，如果利用已有模型算法实现GPU并行计算主要有三个方法：(1) 利用PGI FORTRAN编译器实现CUDA并行计算，2012年，陈昞睿等人^[13]已经通过该方法实现了三维非结构网格水动力模型的建立；(2) 利用OpenACC API直接将原始的串行代码转换为CUDA并行程序，2010年，宋振亚等人利用该方法实现了GPU在POP模式中的应用^[14]；(3) 从底层代码着手，将所有程序转换为CUDA-C代码^[15]。其中前两种方法的实现相对简单，但是无法充分利用GPU的性能，对并行程序的优化也不如CUDA C，因此本文GPU-OM模型选择第三种方法。

并行化实现过程中首先需要将串行程序转换为CUDA程序，模型计算流程如下图所示：

在模型建立的初期，采用了较为简单的方法，将串行程序中所有需要循环计算的部分进行转换(如图 4所示)。

图  3  GPU-OM流程

Fig.  3  GPU-OM flowchart

下载: 全尺寸图片 幻灯片

图  4  串行程序与CUDA C代码对比

Fig.  4  Comparison of sequence execute program and CUDA C program

下载: 全尺寸图片 幻灯片

在整个计算流程中，相对于串行算法的水动力数值模型，GPU-OM时间迭代计算开始之前完成对变量数据的初始化工作，并利用CUDA内建存储器管理函数将内存传递到显存中。为提高效率，使所有的计算工作在GPU端运行，因此在时间迭代计算过程中尽可能的减少CPU-GPU数据交互操作，仅在计算结果需要输出的时候，将所需输出的结果再由显存传递回内存并行输出到结果文件中。

由于在建立GPU-OM的初期仅简单的将串行程序移植为CUDA C并行模式，GPU设备的并行能力受到设备存储带宽以及对全局变量访问延迟的限制而没有充分发挥，为了进一步提高并行效率，对GPU-OM从存储器合并访问，只读缓存等多个方面进行了优化。为了减少由于对全局存储器访问导致的访问延迟，在多次使用到全局存储器变量的kernel函数中先将全局存储器变量传递寄存器，待全部计算结束后，再将寄存器变量传递回全局存储器。除此外将串行模型中的二维数组转换为一维数组，其好处是，(1) 线程和一维数组的对应关系简单，易于编程实现；(2) 每个block上线程数目可设置为线程束(每个线程束包含32个线程)的整数倍，而流处理器恰好以恰好是线程束为最小单位实现上并发执行，这样可以更有效的利用GPU的计算资源。以上方法有效的提高了kernel函数对全局存储器变量的读取效率。

1.4.2   GPU并行程序优化方法

CUDA C的并行程序中，大部分是对全局存储器进行读写操作，而存取延时为400-800 clock cycles非常容易成为性能瓶颈，因此对全局存储器的访问优化是提高GPU并行计算性能的有效方法。本文提出以下几个方法对全局存储器存取进行了优化。

利用寄存器优化全局存储器变量。在一个设备端函数中，如果多次调用全局存储变量，则会在每一调用中都会发生存取延时状态。寄存器作为局部存储器，存取延时仅有1 clock cycles。为了能够有效避免由于数据吞吐量造成的性能瓶颈，本文模型在实现过程中，先将重复读取的全局变量赋给寄存器中的局部变量，在之后计算过程中，仅对局部变量进行存取操作。虽然寄存器的存取速度极快，由于硬件的限制，在每个SM中仅为其分配了有64K空间，对于海量线程的GPU并行程序，这明显是不够用的。因此，合理利用其它存储器(如常数存储器和共享存储器)也是十分必要的。

利用只读数据缓存优化全局存储器变量。CUDA C编译器对只读数据缓存的访问延迟也同样几乎可以忽略，相对于存取全局存储器同样具有明显优势。常量内存即为只读存储数据类型，如果能有充分利用常量内存替代对全局存储变量，可以在一定程度上有效减少数据读取延迟。在本文中，在开始时间迭代计算前，将所需的常量如网格数、节点数等变量传递给设备端的常数存储器。利用这一方法，一方面降低了对寄存器的需求，另一方面在每一次调用kernel函数是的数据传递操作。

将在host端动态分配的设备指针写入一个静态的符号。在大多数CUDA程序中，通常是将设备指针以函数参数的方式传递给设备端函数，当函数需要存取的变量过多，则对该函数的调用和参数设置也变得繁琐，如果将函数分割成多个子函数，又将导致指令发送的冗余以及(1) 和(2) 优化措施的效果降低。因此，本文模型在时间迭代计算前，将设备端数据的指针传递给一个静态符号，这样在同一子文件中的kernel均可以直接调用该变量，而不需要以参数的形式将数据指针传递给kernel，有效的提高了提高的程序的简洁性和易读性。

利用原子性操作实现多个线程对同一数据的写入操作。在本文模型计算动量方程的对流项时，平面积分转换成了沿控制单元边的线积分，需要N×3次循环计算工作，而对单元边做循环计算时，一个公共边的计算结果可以给相邻两个控制单元调用，计算量减少了一般。然而在并行计算过程中，容易出现多个线程同时对同一数据进行操作，如果程序设计不合理，将会导致不同线程对数据操作的混乱，返回错误的计算结果或仅返回一个线程计算结果。在MPI和OpenMP并行算法中通常采用的是加锁的方法，而CUDA C中用原子操作替代了加锁的方法，即利用atomicAdd等函数有效的解决该问题。

2   结果与讨论

2.1   计算区域及边界条件

在本文之前的工作文献[16]中，利用上述方法建立了基本的二维水动力数值模型，和典型算例的解析及进行对比，验证了数值模型的准确性。但是在Ippen理论的算例中，地形为平底的矩形水。

为了验证数值模型在自然情况复杂的水域适用性，在这一章将GPU-OM应用于数值模拟岸线和地形复杂水域的水动力问题。岛群由于岸线曲折多变，岛与岛之间多潮汐水道，地形变化剧烈，非常适用于检验数值计算模型的适用性。本文选择长山群岛海域，采用GPU-OM模型数值计算结果与岛屿之间不同时期潮位、流速、流向实测数据及MIKE21软件计算结果进行对比分析，从而验证模型的有效性。同时与MIKE21软件MPI及GPU并行计算效率进行了比较分析，验证模型的加速效率。

模型计算区域为北黄海部分海域，范围包括整个长山群岛，如图 5所示。

图  5  计算海域位置及网格划分

Fig.  5  Position and the mesh of simulate area

下载: 全尺寸图片 幻灯片

网格剖分：由于计算区域内岛屿众多，地形复杂，为了更好的拟合模拟区域岸线，对岛屿及其附近海域的网格进行了加密，网格分辨率在岛屿区域附近为20 m，在开边界位置处为1000 m左右，整个计算域内网格单元数为161288，网格节点数为83526，如图 5所示。

开边界：模型的开边界条件采用潮位控制，计算时给定边界不同节点处的潮位值，潮位值通过中国海域潮汐预报软件Chinatide^[17]获得。

模型采用自适应时间步长的方法计算时间步长(如式28所示，ΔL为单元中心到节点的最小值)。

曼宁系数根据下式计算得出：

其中：d为单元中心水深。

2.2   模型验证

2.2.1   验证资料说明

模型的验证分为潮位验证和流速验证。

由于计算域较大，难以获得整个海域内同步的潮位及流速、流向监测数据，模型验证期间选择不同时期、不同位置实测潮位和流速、流向数据进行模型的验证，验证点位置见图 6。

图  6  观测站位置

Fig.  6  Location of observation station

下载: 全尺寸图片 幻灯片

验证时间段：1#~4#测点测流时间为2010年11月22日~2010年11月23日，即农历十月十七~十八；5#~11#测点测流时间为2013年10月21日~2013年10月22日，即农历八月十七~十八。

每个验证时间段内对计算区域进行10~15 d的流场计算。计算的初始条件为t=0时，ζ=3.5 m，u=v=0 m/s，其中u, v为水平流速的x, y方向的分量，ζ为初始潮位。计算区域由静止开始，大于12 h后全场流态达到稳定。计算时每小时输出一次潮位、流速和流向计算结果。鉴于文章篇幅限制，文中仅给出代表时段、代表验证点潮位、流速和流向的对比结果。

2.2.2   潮位验证

不同时间段的潮位验证过程曲线，如图 7所示。由图可以看出，长山群岛海域的潮汐属于非正规半日潮，一日潮位过程包括两个涨潮、落潮过程，潮位过程的高低潮不等现象明显。且高低潮时间的相位偏差小于0.5 h，最高最低潮位值偏差小于10 cm，计算结果符合《海岸与河口潮流泥沙模拟技术规程》^[18]要求。

图  7  小长山海洋站潮位观测数据与计算结果的对比

Fig.  7  Comparison of observed and simulated tidal elevation at the xiao chang shan tidal station

下载: 全尺寸图片 幻灯片

2.2.3   流速、流向验证

计算时按小时输出验证潮流站的流速、流向变化过程，用于对流速、流向验证。图 8为分别利用本文模型和MIKE 21在各观测点垂线平均流速、流向验证图。由图可见，本文模型在各验证点计算流速和实测资料吻合较好，计算结果与实测憩流时间和最大流速出现的时间偏差小于0.5 h，流速过程线的形态基本一致，涨、落潮段平均流速偏差小于10%。

图  8  不同测站处流速、流向观测数据与计算结果的对比

Fig.  8  Comparison of observed and simulated tidal speed and direction at different stations

下载: 全尺寸图片 幻灯片

从潮位过程及流速、流向验证结果看，数值模型能够客观反映长山群岛及附近海域的潮流运动规律。

2.3   流场分析

图 9给出了长山群岛海域现状条件下大潮涨急、落急时刻流场分布。

图  9  长山群岛海域大潮期涨急、落急时刻流场

Fig.  9  Current field of rose and snap of spring tide

下载: 全尺寸图片 幻灯片

由计算结果可以看出：涨急时刻，大长山岛北侧海域涨潮流自东向西流动，涨潮流绕过大长山岛西北角后流向为SW向，北侧海域的流速分布比较均匀，基本在0.7 m/s左右，大长山岛西北角处的流速较大，最大达1.4 m/s；大长山岛南侧海域涨潮流沿着大长山岛南岸自东向西流动，在大长山岛西南角处与北侧而来的涨潮流相汇，并朝西南方向流去，其中，大长山岛与小长山岛间的水道中流速较大，最大达2 m/s。

落急时刻，落潮流绕过大长山岛西北角后在北侧海域自西向东流动，北侧海域的流速基本在0.6 m/s左右，普遍弱于涨急时刻的流速，大长山岛西北角处的流速较大，最大达1.0 m/s；在大长山岛南侧海域，从小长山岛与哈仙岛间流入的涨潮流分成了两股，分别沿着大长山岛南岸向东和向西流动，南侧海域的落潮流速明显的弱于涨潮流速，大长山岛与小长山岛间的水道中流速较大，最大达1.7 m/s。

2.4   模型的加速分析

本文模型在文献[16]在应用于矩形平底水池潮驱边界数值模拟中，当网格个数增加到20万个以上时，利用GTX 460显卡(拥有336个核心，显存位宽115 GB/s，单精度浮点计算能力1360 GFLOPS)的GPU并行计算能力已经接近使用20核集群机的计算能力，实现了将加速比提高了十倍以上，并行加速效果十分显著。并行加速比的变化随着网格个数的增加成线性增长。

随着硬件水平的发展，本文在NVIDIA的专业计算显卡Tesla K20X下对模型在地形复杂水域的并行计算加速能力进行了测试。每个Tesla K20X拥有14组SMs(Stream Multiprocessors)、2688个流处理器，同时拥有384-bit的显存位宽和1.5MB的二级缓存，搭配6GB GDDR5显存容量，频率方面核心定在732 MHz，单精度浮点性能为3.93 TFLOPS，显存频率则是5.2 GHz，带宽达249.6 GB/s。CPU并行测试平台采用Intel Xeon E5-2690 v3 2.6G Hz，单精度浮点计算性能为0.1856 TFLOPS。

为了验证本文模型的高效性，在已有工作的基础上，将GPU-OM(分别在GTX 460、GTX 780M和Tesla K20X上运行)与FVCOM模型(分别开启了6、10和32线程)在解决矩形水池潮驱边界问题的并行加速能力进行了对比，如下图所示：

图  10  在不同网格数情况下，GPU-OM与FVCOM在不同设备上模拟潮驱算例时的加速能力

Fig.  10  Speedup of GPU-OM and FVCOM simulations for a variety of meshes and devices in the tidal induced case

下载: 全尺寸图片 幻灯片

由图中结果可以看出，随着网格个数的增加，相对于FVCOM模型，桌面显卡GTX 460和笔记本显卡GTX 780M的计算性能高于同期主流CPU的并行计算能力。专业计算显卡Tesla K20X的计算能力明显高于FVCOM在拥有两颗E5-2690的服务器上并行计算的能力。由此可见，GPU-OM相对于在同期设备上利用CPU并行计算的水动力模型，在计算性能方面具有明显优势。

为了进一步验证在地形和岸线复杂水域GPU-OM的高效计算能力，本文对比了GPU-OM与MIKE21基于MPI和GPU两种并行加速模式在模拟长山群岛海域所用时间分别进行了对比，如下图所示：

图  11  GPU-OM与MIKE21的MPI和GPU并行模式计算时间对比

Fig.  11  The comparison of time consuming of GPU-OM, MIKE21 MPI and GPU parallel mode

下载: 全尺寸图片 幻灯片

由上图可以看出，在长山群岛海域水动力情况模拟过程中所用时间不到单线程版本的1/20，即并行加速能力达到20倍之多。GPU-OM的计算效率优于MIKE21软件利用48线程MPI并行计算的结果，略低于MIKE21 GPU并行计算的效率。不考虑本文模型和MIEK21软件在数值算法复杂程度上的不同，单从计算性能上而言本文模型与MIKE21 GPU数值模型并没有过多差距。计算水域被剖分成161288个网格单元，还不不足以充分利用Tesla K20X的加速性能，在更高分辨率的计算水域，GPU并行计算性能将会得到进一步的提升，相对于MPI并行加速的水动力模型有明显优势。

有上述结果可以看出，无论是与开源的FVCOM模型相比，还是与成熟的商业软件MIKE 21相比GPU-OM在计算性能方面毫不逊色。

3   结论

(1) 本文在CUDA框架下实现了基于GPU的非结构网格、有限体积法的二维水动力数值模型GPU-OM，并利用GPU-OM对长山岛附近海域进行了二维水动力数值模拟。通过对实测水文数据的计算值进行比较验证了该模型的精度，并与FVCOM模型和MIKE 21软件计算的效率进行了对比，证明该模型的高效性。从计算效果看，该模型在研究近岸高分辨率水动力数值模型上有着独特的优越性和高效性。在网格个数相对较少的情况下，Tesla K20高性能计算显卡的能力无法得到完全的发挥；随着网格数的继续增加，并行加速效率将得到显著提高，接近单线程模型的27倍。

(2) 本次数值模拟的精度达到了要求，取得了好的模拟效果，并成功的将GPU-OM应用于长山岛附近海域的海洋环境水动力数值模拟中，为该模型的进一步发展和完善提供借鉴，下一步工作将会实现三维的GPU-OM模型应用于地形和条件更加复杂和海域，并对GPU并行算法做进一步的优化。