GRAPES (Global/Regional Assimilation and PrEdictions System) 模式是中国气象局自主研发的静力/非静力多尺度通用数值预报系统，该模式作为国家重点科技攻关项目，将成为中国气象与气候研究的基础和核心。模式计算作为模式的核心，其巨大的计算需求，给计算机系统提出了挑战。如何减少模式计算耗时，提高模式的运行效率是一直以来的研究目标。GRAPES模式包括模式动力框架和物理过程两部分，动力框架中计算时间最长的为半拉格朗日计算和求解Helmholtz大型线性代数方程组^[1-2]，物理过程中耗时较大的为辐射和微物理过程，其中辐射过程计算耗时占到这4部分的40%左右，所以提升辐射过程的性能和计算效率，对整个GRAPES模式的性能提升具有重要科研价值和实际意义。

过去十几年中，气象数值预报模型得益于大规模集群系统以及单个CPU性能的提升。但是随着近几年气象和气候模型向着大规模、长时间、高精度发展，数据计算需求也呈现爆发性增长, 但CPU单处理器在最近几年的性能提升受到很多限制。这些限制一方面来自于通用计算过程中低的指令级并行；另一方面来自于集成电路功率消耗的物理限制^[3]，并且大规模的集群系统的通信开销将使得系统的可扩展性受到限制。如Michalakes^[4-5]提到使用15000颗IBM Blue Gene的处理器仅获得了有限的性能扩展。这些局限将迫使气象模型的优化着力于提升处理系统的计算密度，更充分的细粒度并发，更大的吞吐率和系统访存速度。

近几年图形处理器强大的通用计算能力逐渐被人们所重视，它有效利用晶体管资源设计大量的处理单元，提升系统的细粒度并发度。随着NVDIA CUDA技术的发展，可编程的图形处理器已经发展成为一种高度并行化、多线程、多核的通用处理平台，具有杰出的计算功率和极高的存储带宽。目前，通用图形处理器技术已经被应用到数学运算、金融分析、医学检查、气象预测、电子线路设计等科学计算领域，大幅度提高了程序的运行速度。

通过通用图形处理器系统，对气象预报模式进行加速计算是近年来的新兴技术，并取得了一定成果。Govett^[6]及Henderson^[7]采用通用图形处理器技术对NOAA的Non-hydrostatic Icosahedral Model (NIM) 模式进行了并行试验，分别检测了F2C-ACC, HMPP, PGI编译器在Fermi GPGPU和Tesla GPGPU上的加速效果，结果显示总体效果可以达到10倍左右。Ruetsch等^[8]联合NVIDIA公司对Weather Research and Forecast (WRF) 模式中的长波辐射传输方案通过通用图形处理器技术进行了测试，结果显示整个辐射方案计算速度提高了10倍。

本研究通过应用NVIDIA CUDA技术，以GRAPES全球模式长波辐射方案为例，对其在通用图形处理器系统上的大规模细粒度并发设计和优化进行研究。在原有消息传递 (MPI) 编程模型任务粗粒度并行的基础上，通过该技术进行更细粒度的线程级并行，以加快程序的运行速度，提高模式的运行效率。

通用图形处理器 (GPGPU) 系统采用单指令流多线程的 (SIMT) 并行执行模型，它类似于单指令流多数据流 (SIMD) 模型，1条指令处理多个数据单元。但不同于SIMD，SIMT可以针对单个线程控制其执行和分支路径，从而可以同时为线程级并发和数据级并发提供支持，并且硬件支持的多线程执行可以提供零开销的上下文切换。2006年，NVIDIA公司提出了CUDA (Compute Unified Device Architecture) 架构, 作为通用计算的统一平台，提供并行编程模型和指令集，可以利用NVDIA GPGPU的计算能力，解决很多并行计算问题。CUDA平台提供一整套的开发环境，包括C语言风格的CUDA C, 以及相关编译器、库和应用程序编程接口 (API) 支持等。CUDA将程序分为执行于CPU上和GPGPU上两部分，后者通过编译，将并行数据和高密度计算过程装载到设备上，启动大量并发硬件线程，进行高速处理，主机和设备之间提供专用API进行数据传输和访存。

通用图形处理器系统的优势在于大量的硬件计算单元可以提供大量的并发处理能力，比如NVIDIA Tesla C1060可以提供峰值超过933 GFlops的单精度浮点计算。GPGPU的并行执行能力和计算能力与其硬件结构是分不开的。NVIDIA Tesla C1060具有240个CUDA处理核心，每个核心具有完整的整数流水线 (ALU) 和浮点运算单元 (FPU), 可以执行1个线程的运算，每8个核心被组织成1个流多处理器SM (Stream Multiprocessor)。同时，每个SM具有可配置的64 KB共享内存或一级缓存，以及32 KB×32的寄存器堆。整个架构也同时针对双精度浮点计算进行设计，可以提供超过78 GFlops的双精度浮点计算能力。

通用图形处理器系统的特点在于并行执行，通过其浮点运算能力和高度并行化的结构来实现大规模的并行计算。一般来说，只要保持足够多的线程处于运行状态，较传统CPU就会有数十倍加速，如果应用程序的控制流和数据流适合，并且经过充分优化就可以达到数百倍的加速。通用图形处理器系统的计算能力是有针对性的，具有以下特点：大量的平行数据与计算，较少的逻辑分支判断等^[9-11]。这些计算特点与数值天气预报模式数据密集、计算量大、且高度并行的特点相吻合。因此，将GRAPES全球模式用通用图形处理器来实现并行，以加速计算、提高模式效率是可行的。

GRAPES模式的并行设计是在模式软件框架结构基础上进行的，针对采用消息传递的分布式集群系统^[12-14]。GRAPES模式是一个格点模式，它采用区域分解的方法来并行处理，在计算区域的并行分解上采用了通用的水平经纬网格划分。每个节点并行处理各子区域内的计算负载，当最后一个节点处理完成后进行条件判断，进入下一个循环阶段或结束循环^[15]。每个节点的计算程序相同，但数据不同，这些数据格式相同，高度并行，可通过通用图形处理器技术对程序进行更细粒度的线程级并行，即对每个节点挂载1个图形处理器 (图 1)，对每个节点所分配的数据进行更细粒度的并行，这样每个节点的循环次数大大降低，提高了程序的并发度。

图 1

图 1. 通用图形处理器并行系统架构 Fig 1. The architecture of GPGPU parallel system

在通用图形处理器系统中，可以为每个节点挂载1个GPGPU，GPGPU是该节点的协从处理器。一般来说，如果1个GPGPU相对于CPU的1个处理器节点有m倍加速，那么为CPU的n个处理器各配置1个GPGPU，就会有m×n倍的总加速。本文中采用1个GPGPU与CPU的1个处理器节点进行对比试验，以后的工作中可以引进更多的GPGPU配合更多的CPU节点进行试验，将多个GPGPU并行与CPU多个节点的MPI任务并行进行对比。

GRAPES_GLOBAL模式的辐射模型采用欧洲中期天气预报中心 (ECMWF) 的长短波辐射方案^[15-16], 该方案是Morcrette根据法国里尔大学辐射方案发展的更新版本，其中长波方案使用了覆盖6个谱区域 (0~2080 cm^-1) 的宽带通量发射率方法，这6个谱区域分别对应水汽的旋转和振动旋转谱带的中心、二氧化碳的15 μm谱带、大气窗、臭氧的9.6 μm谱带、25 μm窗区和水汽振动旋转普带的翼。该方案在CPU上执行的长波辐射过程计算量庞大，该方案中地球被切分为i×j个网格点，对应i×j个气柱，根据处理器个数可将全球划分为若干子区域，多个处理器并发执行各子区域的辐射方案，(图 2a), 每个处理器计算该区域内各气柱在不同层次的辐射量时串行执行，需循环N=i×j×k次。由于各区域内气柱之间相互独立，各气柱均要执行辐射方案中相同的计算过程，有很大的并行性，利用GPGPU系统进一步并发计算各子区域内每个气柱的辐射量，1个线程计算1个气柱 (图 2b)，此时则只需进行k次循环，即N=k。

图 2

图 2. CPU集群粗粒度 (a) 和GPU细粒度 (b) 并发方案(t0，……，tm为线程数) Fig 2. CPU coarse-grained parallelism (a) and GPU fine-grained parallelism (b) (t0, …, tm denote the number of thread)

对比以上两个方案可以看出GPGPU细粒度并发方案具有明显的优势，表现为以下3个方面：① 并发度增加。并发度表示系统可以同时处理的基本任务数量。在CPU方案中，并发度等于同时执行的CPU进程数据，该数据在一个规模受限的系统中会比较小；而在GPGPU方案中，由于每个GPGPU本身具有很强的多线程并发能力，其并发度将数百倍于CPU系统。当基本任务之间没有任何相关性时，并发度的提升通常可以带来性能的提升。② 充分利用大的显存带宽。由于GPGPU访存带宽远超CPU系统，通过仔细设计的访存模式，可以有效地进行访存合并，将不同线程的数据通过总线同时传输，从而更加提升访存性能。在很多访存受限的任务中，该优势将更加明显。③ 计算和通信密度增加，提升效率。由于GPGPU系统1颗芯片本身集成了数量众多的处理单元，其计算密度和通信密度大大提升，比较传统的分布式集群系统，其通信和核间协作将有一定比例集中于单颗芯片内部，从而为整体效率的提升保留空间。通过仔细设计的通信模式，可以有效地提升计算通信开销比。

图 3为本试验中长波辐射方案执行模型，该模型分为两部分，分别为执行于CPU和GPGPU上的部分。执行于GPGPU上的并行程序为核函数，每个核函数对应1个网格 (grid)，每个网格又可分为多个线程块 (block)，每个线程块可由最多512个线程组成 (Tesla 1060)。由于每个线程块内的线程数量受到寄存器数量及共享内存使用等硬件条件的限制，当block内线程数量较大时有些复杂的核函数就没有执行，故本文中线程块内线程数量设为128个，全球360×180个气柱被划分为507个线程块。线程被组织成网格和线程块的形式，以适应不同的处理维度和硬件配置，并通过线程调度器将线程块分配给不同的SM执行。

图 3

图 3. 长波辐射方案执行模型 Fig 3. The execute model of the long-wave radiation scheme

CPU对5个核函数的调用是串行执行的，启动核函数之前需要在图形处理器上分配显存空间，并进行由CPU至GPGPU的数据传输。程序计算用到的常数需要在计算之前传入GPGPU，存放于固定内存中，固定是只读的，只有64 KB，每个线程都可访问。而大量的输入数据则可通过CUDA API函数传入到全局内存，它是可读可写的，1个核函数执行完成后，结果将自动保存在全局内存中。注意避免分多次传入数据，否则数据传输时间会大大增加。所有气柱计算完成后，将所需数据由显存传输到内存，以便CPU使用，并释放显存空间，具体流程见图 4。

图 4

图 4. 长波辐射方案流程图 Fig 4. The flow chart of the long-wave radiation scheme

本文基于64位的Fedora Linux操作系统，Tesla C1060 GPGPU，CPU型号Intel (R) core (TM) i7-920(主频为2.67 GHz，总线频率为4.8 GT/s，四核八线程，本文采用单核作对比试验)，CUDA 2.3版本，PGI 10.1编译器构成的小型服务器的试验环境。试验采用分辨率为1°×1°的GRAPES全球模式中长波辐射方案，首先采用360个气柱，气柱之间相互独立，1个线程计算1个气柱，对其进行并行测试。通过对程序的深入分析发现，某些子程序在计算时不同气层之间也相互独立，因此可以对气柱与气层同时并行，此时1个线程计算一个空间点。表 1是以上两种方法的测试结果对比。

表 1

表 1 并行结果对比 Table 1 The comparison of parallel results 核函数 并行计算气柱 用时/ms 并行计算空间点 用时/ms inatm_d () 17.41 2.16 cldprmc_d () 16.21 6.63 setcoef_d () 2.65 0.85 taumol_d () 23.06 3.68 rtrnmc_d () 37.14 37.14

表 1 并行结果对比 Table 1 The comparison of parallel results

Tesla C1060具有30个SM, 每个SM理论上可以同时调度8个线程块，或32个线程束 (warp)，线程束是GPGPU调度的基本单元，该数目受到每个线程的寄存器数目、共享内存使用等因素制约。当使用360个气柱时，由于线程数量较少，每个线程块配置64个线程，共计6线程块，从而统计有6×64/32=12个线程束, 而硬件理论最多可以同时调度30×32=960个线程束, 系统的占用率从理论上讲只有1.25%。有些计算气层之间相互独立，对气柱与气层同时并行，配置每个线程块具有64个线程，但是有37×6个线程块 (37为气层数量)，从而原有的线程束数目增加到37×6×64/32=444，于是可以有足够多的线程束使SM处于繁忙状态。当然，系统的实际占用率会随着应用不同而变化，但是提高硬件的占用率在很多情况下可以获益。试验表明，新的方式计算用时减少到50.46 ms。

结果表明，第2种方法较第1种方法有明显提升，比原来的配置获得超过2倍的提升，可以看出保持较多数量的线程执行，通过大量的线程束来隐藏系统的访存延迟以及保持SM的较高占用率，在本应用中会有效提升性能。在采用MPI技术进行大规模数据并行的系统上采用GPGPU加速时，由于气柱数量有限，需要在垂直方向或其他方法进一步增加模式的并发性。但当气柱数目增加到一定程度后，系统的占用率也会随之提高，有足够多的线程束使SM处于繁忙状态，此时对气柱与气层同时并行的加速效果不是很明显。

由于在GPGPU上运行的总线程数量受到显存大小的限制，显卡满负荷运行1次长波辐射方案，最多可执行54×360个气柱，因此本试验将分任务调用4次长波辐射函数，以完成整个长波辐射的运算。图 5显示了气柱数目从0，5×360, ……, 50×360, 54×360时CPU与GPGPU计算耗时的变化情况。可以看出，随着气柱数目的增加CPU计算耗时呈直线迅速增加，大大超过了GPGPU用时的增加程度，线程数量增加的越多，GPGPU的计算优势将越明显。本试验中CPU与GPGPU之间的数据传输量为1.6 GB，它们之间的带宽为2.85 GB/s，由于硬件限制，数据传输用时和GPGPU计算用时相当。

图 5

图 5. 随着气柱数目增加CPU与GPGPU计算耗时比较 Fig 5. The comparison of computational speed between CPU and GPGPU with the increase number of columns

图 6为通用图形处理器系统计算后输出的全球长波辐射通量与模式输出值的对比图，经验证误差在允许范围内，其相对误差的量级在10^-5至10^-4之间。在保证结果正确的前提下，完整运行一次长波辐射方案初步得到11倍加速。原模式中长波辐射方案用时为38817.56 ms；而基于通用图形处理器系统的长波辐射方案用时为3514.72 ms; 数据传输用时为3606.67 ms。

图 6

图 6. CPU计算的全球长波辐射通量 (a) 与GPGPU计算输出值 (b) 对比 Fig 6. The comparison of the long-wave radiation fluxes results between CPU (a) and GPGPU (b)

可以看出，数据传输用时占了GPGPU总用时的50%以上，当然这与所用硬件带宽有很大关系。Govett^[6]提出将模式中的动力框架、物理模型的所有程序都放在GPGPU上运行，这样可以消除每个积分步的数据拷贝，而只需要在输入、输出及进程间的通讯进行数据拷贝，从而提高模式运行效率。异步传输允许在GPGPU进行计算的同时，进行主机和设备间的通信，实现流式处理，可使用该方法隐藏数据传输时间，在下一步的工作中将运用数据异步传输和其他方法对程序进一步优化。

本试验采用NVIDIA CUDA技术，对GRAPES全球模式中长波辐射传输方案在GPGPU上进行了初步并行试验，得到以下结论：

1) 采用GPGPU对数值预报部分计算进行加速是可以实现的，其中辐射计算在GPGPU (Tesla 1060) 上的运算速度可以达到Intel (R) core (TM) i7-920的10倍以上，相对于CPU来说有着明显的计算优势。

2) 基于通用图形处理器的多线程并行在不超过硬件限制的条件下，线程数量越多，速度越快，增加模式的并发度有利于发挥GPGPU的性能。

3) 在核函数的调用中花费大量的时间来进行数据传输，会消耗由通用图形处理器的计算能力带来的优势，可通过异步传输等方法进行优化来隐藏数据传输时间。

后续工作包括采用最新的NVIDIA GPGPU进行移植，采用各种方法继续对程序进行优化以期得到更好的加速效果，并将其推广到模式中的其他物理模块中，以提高整个模式的运行效率，还可以引进更多的GPGPU，将GRAPES模式并行计算部分进行多节点的对比试验。