GRAPES全球模式MPI與Open MP混合并行方案

蔣沁谷金之雁

1）（中國(guó)氣象科學(xué)研究院，北京100081）2）（中國(guó)氣象局?jǐn)?shù)值預(yù)報(bào)中心，北京100081）

GRAPES全球模式MPI與Open MP混合并行方案

蔣沁谷1）金之雁2）*

1）（中國(guó)氣象科學(xué)研究院，北京100081）2）（中國(guó)氣象局?jǐn)?shù)值預(yù)報(bào)中心，北京100081）

隨著多核計(jì)算技術(shù)的發(fā)展，基于多核處理器的集群系統(tǒng)逐漸成為主流架構(gòu)。為適應(yīng)這種既有分布式又有共享內(nèi)存的硬件體系架構(gòu)，使用MPI與Open MP混合編程模型，可以實(shí)現(xiàn)節(jié)點(diǎn)間和節(jié)點(diǎn)內(nèi)兩級(jí)并行，利用消息傳遞與共享并行處理兩種編程方式，MPI用于節(jié)點(diǎn)間通信，Open MP用于節(jié)點(diǎn)內(nèi)并行計(jì)算。該文采用MPI與Open MP混合并行模型，使用區(qū)域分解并行和循環(huán)并行兩種方法，對(duì)GRAPES全球模式進(jìn)行MPI與Open MP混合并行方案設(shè)計(jì)和優(yōu)化。試驗(yàn)結(jié)果表明：MPI與Open MP混合并行方法可以在MPI并行的基礎(chǔ)上提高模式的并行度，在計(jì)算核數(shù)相同的情況下，4個(gè)線程內(nèi)的MPI與Open MP混合并行方案比單一MPI方案效果好，但在線程數(shù)量大于4時(shí)，并行效果顯著下降。

混合并行；數(shù)值天氣預(yù)報(bào)模式；區(qū)域分解；循環(huán)并行

引 言

過(guò)去十幾年里，隨著CPU性能不斷提升及并行計(jì)算編程模型的日趨成熟，氣象數(shù)值模式得到長(zhǎng)足發(fā)展。目前大多氣象數(shù)值模式使用消息傳遞接口（message passing interface，簡(jiǎn)稱(chēng)MPI）編程模型，但隨著氣象模式向高分辨率和模擬更加真實(shí)的物理動(dòng)力過(guò)程方向發(fā)展，數(shù)據(jù)計(jì)算和存儲(chǔ)需求成倍增長(zhǎng)。如氣象模式COSMO（consortium for small-scale modeling）分辨率從2 km×2 km提高到1 km×1 km，計(jì)算量大約增加10倍［1］。數(shù)值天氣預(yù)報(bào)模式為達(dá)到實(shí)時(shí)性要求，使用的處理器數(shù)量不斷增加。但由于受到擴(kuò)展性的限制，MPI應(yīng)用在進(jìn)程數(shù)達(dá)到一定規(guī)模后加速比不再提高，甚至在某些情況下開(kāi)始下降。

隨著多核技術(shù)的發(fā)展，共享內(nèi)存架構(gòu)在高性能計(jì)算機(jī)市場(chǎng)逐漸興起。基于多核處理器的集群系統(tǒng)由于其良好的可擴(kuò)展性及高性?xún)r(jià)比等優(yōu)點(diǎn)逐漸成為主流體系架構(gòu)。如何使應(yīng)用軟件的并行編程模型更加適應(yīng)這種既有分布式內(nèi)存又有共享內(nèi)存特點(diǎn)的硬件體系架構(gòu)，從而發(fā)揮硬件最大性能成為目前研究的一個(gè)重要方向。目前分布式集群下最常用的并行編程模型是MPI。但對(duì)于既有分布式內(nèi)存又有共享內(nèi)存特點(diǎn)的多核處理器集群，采用單一MPI編程模式往往不能取得最理想的性能。針對(duì)集群中共享內(nèi)存特點(diǎn)，可以考慮引入共享內(nèi)存編程的實(shí)際工業(yè)標(biāo)準(zhǔn)Open MP。MPI與Open MP混合編程模式能更好映射多核處理器集群的體系結(jié)構(gòu)，提供節(jié)點(diǎn)間和節(jié)點(diǎn)內(nèi)兩級(jí)并行，充分利用消息傳遞與共享內(nèi)存兩種編程模式的優(yōu)點(diǎn)，MPI用于節(jié)點(diǎn)間粗粒度通信，而Open MP用于節(jié)點(diǎn)內(nèi)計(jì)算，減少節(jié)點(diǎn)內(nèi)通信開(kāi)銷(xiāo)，獲得優(yōu)于單一MPI并行方案應(yīng)用的性能和擴(kuò)展性［2-3］。

氣象數(shù)值模式作為大規(guī)模計(jì)算和存儲(chǔ)密集型的典型應(yīng)用，一直是高性能計(jì)算領(lǐng)域研究的熱點(diǎn)。國(guó)外中尺度數(shù)值模式 WRF（Weather Research and Forecasting model）已實(shí)現(xiàn)MPI與Open MP混合并行，可運(yùn)行于多種計(jì)算架構(gòu)平臺(tái)［4-6］。?ipková等［7］對(duì)WRF模式MPI與Open MP混合并行方案進(jìn)行試驗(yàn)，結(jié)果顯示混合編程模型適合分布式共享內(nèi)存系統(tǒng)。三維海洋模式NEMO（Nucleus for European Modeling of the Ocean）的 MPI并 行是利用對(duì)水平區(qū)域進(jìn)行分解。為了減少NEMO計(jì)算時(shí)間，引入MPI與Open MP混合并行，對(duì)模式垂直層進(jìn)行Open MP線程并行，效果顯著［8］。

張昕等［9］使用Open MP對(duì)中尺度模式MM5進(jìn)行試驗(yàn)，結(jié)果顯示Open MP并行編程實(shí)現(xiàn)比MPI編程簡(jiǎn)單，也能達(dá)到較高的加速比。朱政慧等［10-12］提出了氣象數(shù)值天氣預(yù)報(bào)模式的MPI與Open MP混合并行編程模型，對(duì)高分辨率有限區(qū)同化預(yù)報(bào)系統(tǒng)（HLAFS）進(jìn)行MPI與Open MP混合并行試驗(yàn)，取得較好的并行效果。郭妙等［13］及鄭芳等［14］分別對(duì) GRAPES（Global／Regional Assimilation and Pr Ediction System）輻射過(guò)程運(yùn)用CUDA（Compute Unified Device Architecture）編程在GPGPU（General Purpose Graphic Processing Unit）設(shè)備上進(jìn)行試驗(yàn)，獲得10倍以上的加速比。樊志杰等［3］對(duì)GRAPES四維變分同化系統(tǒng)進(jìn)行MPI與Open MP混合并行，發(fā)現(xiàn)當(dāng)單一MPI實(shí)現(xiàn)程序并行效率下降到90%以下時(shí)，使用混合并行方案實(shí)現(xiàn)的并行效率比單一MPI實(shí)現(xiàn)程序高5%到10%。

中國(guó)氣象局自主研發(fā)數(shù)值天氣預(yù)報(bào)系統(tǒng)GRAPES目前已成功實(shí)現(xiàn)MPI并行方案并投入業(yè)務(wù)使用。為實(shí)現(xiàn)GRAPES在多種高性能計(jì)算平臺(tái)上應(yīng)用，本研究在GRAPES已有粗粒度的MPI并行基礎(chǔ)上，增加細(xì)粒度的Open MP線程級(jí)并行，最終形成GRAPES的MPI與Open MP混合并行方案（以下簡(jiǎn)稱(chēng)混合并行方案），適應(yīng)多核處理器集群架構(gòu)，獲得比單一MPI并行方案更高的擴(kuò)展性。

1 Open MP編程模型

Open MP［15-17］是 由 Open MP Architecture Review Board推出的針對(duì)共享內(nèi)存的并行應(yīng)用程序接口。Open MP由編譯指令、運(yùn)行時(shí)庫(kù)函數(shù)及環(huán)境變量3個(gè)部分組成。Open MP具有可移植性、可擴(kuò)展性、易于使用等優(yōu)點(diǎn)，被廣泛接受，成為共享內(nèi)存體系編程的實(shí)際工業(yè)標(biāo)準(zhǔn)。Open MP使用注釋型指令，對(duì)于不使用Open MP選項(xiàng)編譯或編譯器不支持Open MP，Open MP編譯指令相當(dāng)于程序注釋語(yǔ)句，該特點(diǎn)可以使同一份源代碼保留串、并行兩個(gè)版本。

線程是Open MP中獨(dú)立執(zhí)行任務(wù)的運(yùn)行時(shí)實(shí)體。Open MP使用fork-join并行執(zhí)行模型。一般情況，程序開(kāi)始以主線程串行執(zhí)行，當(dāng)遇到并行域結(jié)構(gòu)指令時(shí)創(chuàng)建一個(gè)線程組，并行域內(nèi)語(yǔ)句以線程并行執(zhí)行。當(dāng)完成并行域內(nèi)語(yǔ)句，退出并行域結(jié)構(gòu)時(shí)，銷(xiāo)毀線程組，只保留主線程繼續(xù)執(zhí)行。線程數(shù)目可由編譯指令子句、運(yùn)行時(shí)庫(kù)函數(shù)或環(huán)境變量指定。

Open MP基于共享內(nèi)存模型。默認(rèn)情況下數(shù)組為共享（shared）屬性，對(duì)同一線程組中所有線程可見(jiàn)。對(duì)于私有（private）變量，每個(gè)線程擁有數(shù)據(jù)對(duì)象各自的副本（copy），且變量可能在不同線程中值不同，如循環(huán)迭代變量一般為私有變量。

2 GRAPES混合并行方案設(shè)計(jì)

2.1 MPI與Open MP混合并行模型

MPI與Open MP混合并行模型提供節(jié)點(diǎn)間以及節(jié)點(diǎn)內(nèi)兩級(jí)并行，MPI位于頂層，Open MP位于底層。從單個(gè)MPI進(jìn)程角度來(lái)看，每個(gè)進(jìn)程內(nèi)部單線程執(zhí)行。在MPI進(jìn)程內(nèi)開(kāi)啟Open MP并行域（！＄OMP PARALLEL指令）產(chǎn)生線程級(jí)并行，而在Open MP并行域外仍為單線程執(zhí)行。理想情況是MPI用于節(jié)點(diǎn)間粗粒度通信而Open MP用于節(jié)點(diǎn)內(nèi)計(jì)算部分。

2.2 GRAPES混合并行方案

GRAPES模式采用水平區(qū)域分解方法劃分水平經(jīng)緯度網(wǎng)格。目前GRAPES模式針對(duì)分布式內(nèi)存計(jì)算機(jī)系統(tǒng)采用消息傳遞接口MPI。將預(yù)報(bào)區(qū)域（domain）按照計(jì)算機(jī)節(jié)點(diǎn)個(gè)數(shù)做劃分，劃分得到的子區(qū)域稱(chēng)為一級(jí)分區(qū)（patch）。每個(gè)MPI任務(wù)完成一個(gè)一級(jí)分區(qū)的計(jì)算，模式計(jì)算過(guò)程中需要相鄰一級(jí)分區(qū)的變量值，每個(gè)一級(jí)分區(qū)變量存儲(chǔ)空間比一級(jí)分區(qū)稍大，存儲(chǔ)相應(yīng)的一級(jí)分區(qū)及重疊區(qū)的變量，當(dāng)計(jì)算涉及到其他一級(jí)分區(qū)區(qū)域格點(diǎn)值時(shí)，必須事先進(jìn)行相應(yīng)重疊區(qū)的通信，獲得重疊區(qū)變量值［18-20］（圖1）。

圖2是GRAPES模式的計(jì)算流程示意圖［19］。GRAPES模式運(yùn)行首先進(jìn)行一些初始化操作及初始資料的讀入。然后設(shè)置預(yù)報(bào)域（domain）、對(duì)MPI進(jìn)程計(jì)算區(qū)間（一級(jí)分區(qū)）進(jìn)行劃分以及分配變量存儲(chǔ)空間。完成初始化操作后，模式通過(guò)積分控制程序啟動(dòng)積分計(jì)算程序，積分計(jì)算程序是GRAPES模式的核心部分，包括動(dòng)力框架和物理過(guò)程兩部分。在1個(gè)積分時(shí)間步內(nèi)，完成模式動(dòng)力框架積分計(jì)算和諸如輻射、陸面、積云降水、微物理計(jì)算等物理過(guò)程的計(jì)算，同時(shí)完成積分時(shí)間遞增，為下一個(gè)積分時(shí)間步計(jì)算做準(zhǔn)備。積分循環(huán)在積分控制程序中完成，該子程序同時(shí)控制程序的歷史數(shù)據(jù)與后處理數(shù)據(jù)輸出。積分計(jì)算是GRAPES模式主要計(jì)算部分，經(jīng)測(cè)試發(fā)現(xiàn)分辨率為1°×1°的GRAPES模式積分48步，使用64個(gè)進(jìn)程進(jìn)行計(jì)算，積分計(jì)算部分占總體計(jì)算時(shí)間的90%以上，這一比例隨積分步數(shù)的增長(zhǎng)會(huì)進(jìn)一步增加。可以看出，積分計(jì)算部分混合并行效果的好壞將決定GRAPES模式整體并行效果。

圖1 GRAPES混合并行水平區(qū)域分解方案［18］Fig.1 The horizontal domain decomposition scheme of GRAPES hybrid parallel（from Reference［18］）

圖2 GRAPES模式計(jì)算流程圖［19］Fig.2 The calculation flow chart of GRAPES model（from Reference［19］）

本文在GRAPES模式已有MPI并行基礎(chǔ)上，進(jìn)一步對(duì)GRAPES模式進(jìn)行Open MP并行。主要工作是對(duì)占主要計(jì)算時(shí)間的積分計(jì)算部分進(jìn)行Open MP并行，有兩種并行思路［21］：①二級(jí)分區(qū)（tile）并行，在每個(gè)一級(jí)分區(qū)內(nèi)根據(jù)可用線程數(shù)進(jìn)行區(qū)域分解成多個(gè)二級(jí)分區(qū)（圖1），對(duì)二級(jí)分區(qū)進(jìn)行高層次較粗粒度并行；②循環(huán)并行，對(duì)主要計(jì)算循環(huán)代碼塊進(jìn)行普通細(xì)粒度并行。

2.2.1 二級(jí)分區(qū)并行

二級(jí)分區(qū)并行具體做法是對(duì)模式預(yù)報(bào)區(qū)域進(jìn)一步進(jìn)行水平區(qū)域分解，將每個(gè)一級(jí)分區(qū)均勻劃分為多個(gè)二級(jí)分區(qū)，二級(jí)分區(qū)數(shù)目由Open MP運(yùn)行時(shí)線程數(shù)決定，每個(gè)二級(jí)分區(qū)塊分配一個(gè)線程進(jìn)行計(jì)算。對(duì)二級(jí)分區(qū)塊進(jìn)行Open MP線程并行使用默認(rèn)的靜態(tài)線程調(diào)度方式。二級(jí)分區(qū)并行代碼如下所示。

二級(jí)分區(qū)并行具有如下優(yōu)點(diǎn)：首先，對(duì)各個(gè)二級(jí)分區(qū)進(jìn)行Open MP并行，并行層次與MPI并行屬于同一層級(jí)，是較粗粒度并行；其次，對(duì)程序二級(jí)分區(qū)并行所需改動(dòng)代碼工作量小，Open MP并行域內(nèi)變量屬性簡(jiǎn)單明了，子程序內(nèi)局部變量默認(rèn)為私有屬性，傳遞參數(shù)變量為共享變量。

雖然二級(jí)分區(qū)并行有這些優(yōu)點(diǎn)，但對(duì)子程序有一定要求，必須是線程安全，即數(shù)據(jù)對(duì)各個(gè)線程是獨(dú)立的，沒(méi)有輸入和輸出，沒(méi)有MPI通信。另外，二級(jí)分區(qū)并行屬于靜態(tài)的區(qū)域分解方法，不能有效處理一些物理過(guò)程的負(fù)載平衡問(wèn)題。所以有些子程序不適合二級(jí)分區(qū)方式并行處理，例如求解Helmholtz方程和有負(fù)載平衡功能的物理過(guò)程。

2.2.2 循環(huán)并行

循環(huán)并行是對(duì)調(diào)用子程序內(nèi)部的循環(huán)塊進(jìn)行Open MP并行。循環(huán)并行需要深入子程序內(nèi)部，對(duì)計(jì)算循環(huán)內(nèi)所有變量的屬性進(jìn)行說(shuō)明，對(duì)沒(méi)有數(shù)據(jù)相關(guān)的可并行的計(jì)算循環(huán)進(jìn)行說(shuō)明等并行化處理，改動(dòng)地方較多。由于是對(duì)子程序內(nèi)部各循環(huán)代碼塊進(jìn)行細(xì)粒度的線程并行，循環(huán)結(jié)構(gòu)越多，需要的Open MP并行指令越多，并行開(kāi)銷(xiāo)也越大。GRAPES模式中循環(huán)結(jié)構(gòu)多為多重循環(huán)嵌套，并行域置于最外層循環(huán)時(shí)并行結(jié)構(gòu)開(kāi)銷(xiāo)最小。由于開(kāi)啟、關(guān)閉并行域開(kāi)銷(xiāo)很大，并行時(shí)盡可能使并行域最大化。在并行域內(nèi)遇到MPI通信時(shí)，使用主線程單獨(dú)執(zhí)行，在主線程執(zhí)行前后需要同步數(shù)據(jù)。線程的調(diào)度策略不同也會(huì)影響循環(huán)并行效果。

2.3 二級(jí)分區(qū)并行效果

二級(jí)分區(qū)并行根據(jù)線程數(shù)目將每個(gè)一級(jí)分區(qū)劃分成多個(gè)二級(jí)分區(qū)，每個(gè)線程分到一個(gè)二級(jí)分區(qū)進(jìn)行計(jì)算，這是區(qū)域分解做法，不同的區(qū)域分解方案下程序并行效果會(huì)有差異。下面以插值程序phy_prep為例，比較不同的二級(jí)分區(qū)劃分方案下混合并行效果。

本文試驗(yàn)環(huán)境為IBM Flex系統(tǒng)。IBM系統(tǒng)節(jié)點(diǎn)配置是4顆8核Power 7處理器。試驗(yàn)采用GRAPES全球水平分辨率1°×1°，計(jì)算網(wǎng)格規(guī)模為360×181×36，開(kāi)啟物理過(guò)程。二級(jí)分區(qū)并行時(shí)二級(jí)分區(qū)數(shù)根據(jù)程序運(yùn)行時(shí)線程數(shù)環(huán)境變量確定。

首先比較垂直層插值程序phy_prep在一維經(jīng)向、一維緯向以及水平二維二級(jí)分區(qū)劃分方案下各自混合并行單個(gè)積分步內(nèi)所用時(shí)間。試驗(yàn)中使用16個(gè)進(jìn)程，8個(gè)線程。表1為以上3種區(qū)域分解方案并行測(cè)試結(jié)果對(duì)比，可以看出，IBM系統(tǒng)表現(xiàn)出經(jīng)向二級(jí)分區(qū)數(shù)少，緯向二級(jí)分區(qū)數(shù)多，計(jì)算時(shí)間短，說(shuō)明一維緯向劃分方案最好。

表1 不同二級(jí)分區(qū)劃分方案并行計(jì)算時(shí)間對(duì)比Table 1 The comparison of parallel computational time with different tile-decomposition schemes

比較在二級(jí)分區(qū)劃分方案固定為一維緯向劃分時(shí)，不同的一級(jí)分區(qū)劃分方法插值程序phy_prep混合并行單個(gè)積分步內(nèi)所用時(shí)間。表2是4個(gè)一維緯向劃分，不同一級(jí)分區(qū)劃分方案下并行測(cè)試結(jié)果對(duì)比，可以看出，一級(jí)分區(qū)不同時(shí)，二級(jí)分區(qū)并行效果會(huì)有差異。經(jīng)向和緯向一級(jí)分區(qū)數(shù)相同時(shí)并行效果最好。對(duì)于IBM系統(tǒng)，由于東西方向?yàn)樽顑?nèi)層計(jì)算循環(huán)，它的長(zhǎng)度直接影響計(jì)算向量的長(zhǎng)度，對(duì)計(jì)算速度有較大影響。

表2 不同一級(jí)分區(qū)劃分方案下計(jì)算時(shí)間對(duì)比Table 2 The comparison of computational time with different patch-decomposition schemes

2.4 循環(huán)并行效果

以插值程序phy_prep為例，對(duì)其進(jìn)行循環(huán)并行，對(duì)比其在靜態(tài)、動(dòng)態(tài)及指導(dǎo)調(diào)度（guided）3種線程調(diào)度方式下隨線程數(shù)增加單個(gè)積分步內(nèi)所用計(jì)算時(shí)間。插值計(jì)算在3種線程調(diào)度策略下的測(cè)試結(jié)果對(duì)比見(jiàn)表3。由表3可以看出，IBM系統(tǒng)環(huán)境下2個(gè)線程內(nèi)動(dòng)態(tài)調(diào)度比靜態(tài)調(diào)度效果好，當(dāng)4個(gè)線程以上時(shí)動(dòng)態(tài)調(diào)度比靜態(tài)調(diào)度差，而指導(dǎo)調(diào)度方式在4個(gè)線程內(nèi)比靜態(tài)調(diào)度用時(shí)少，線程為8時(shí)比靜態(tài)用時(shí)多一些。

以上垂直層插值程序phy_prep中每個(gè)格點(diǎn)計(jì)算量比較均勻，對(duì)于積云對(duì)流參數(shù)化和云微物理過(guò)程等隨天氣情況變化的計(jì)算過(guò)程會(huì)有負(fù)載平衡問(wèn)題［22］。下面測(cè)試格點(diǎn)計(jì)算量不均勻情況，以積云對(duì)流為例，比較靜態(tài)、動(dòng)態(tài)及指導(dǎo)3種線程調(diào)度方式下循環(huán)并行隨線程數(shù)增加單個(gè)積分步內(nèi)所用計(jì)算時(shí)間。積云對(duì)流參數(shù)化計(jì)算時(shí)間在3種線程調(diào)度策略的測(cè)試結(jié)果對(duì)比見(jiàn)表3，可以看出，由于格點(diǎn)的計(jì)算量不太均勻，更適合采用動(dòng)態(tài)和指導(dǎo)調(diào)度方式，指導(dǎo)調(diào)度比動(dòng)態(tài)調(diào)度方式稍好。

表3 3種線程調(diào)度方式下插值與積云對(duì)流參數(shù)化計(jì)算時(shí)間Table 3 The comparison of interpolation and cumulus convection scheme computational time with three thread scheduling policies

從以上試驗(yàn)結(jié)果可以看出，程序循環(huán)并行時(shí)，線程的調(diào)度方式與計(jì)算特性密切相關(guān)，不能一概而論，需要對(duì)于具體計(jì)算過(guò)程具體分析，選擇合適的線程調(diào)度方式。對(duì)于計(jì)算量均勻分布的格點(diǎn)計(jì)算，可選用靜態(tài)線程調(diào)度策略，對(duì)于計(jì)算量非均勻分布的格點(diǎn)計(jì)算，采用動(dòng)態(tài)和指導(dǎo)調(diào)度的效果較好。

2.5 二級(jí)分區(qū)并行和循環(huán)并行對(duì)比

從上面的分析可以看出，二級(jí)分區(qū)并行實(shí)際上是循環(huán)并行的一種特殊情況，針對(duì)一級(jí)分區(qū)進(jìn)行區(qū)域分解，并行層次較高，它要求并行的整個(gè)子程序內(nèi)部數(shù)據(jù)線程安全，不涉及MPI通信。下面對(duì)比垂直層插值程序phy_prep使用二級(jí)分區(qū)和循環(huán)并行單個(gè)積分步內(nèi)所用計(jì)算時(shí)間。由于循環(huán)并行是對(duì)循環(huán)體最外圍的緯向循環(huán)進(jìn)行并行，為了對(duì)比的合理性，二級(jí)分區(qū)并行使用一維緯向二級(jí)分區(qū)劃分方案，循環(huán)并行使用靜態(tài)線程調(diào)度策略。表4是兩種并行方法的測(cè)試結(jié)果對(duì)比?？梢钥闯觯怪睂硬逯党绦騪hy_prep循環(huán)并行效果比二級(jí)分區(qū)并行略好。

綜合考慮，GRAPES模式混合并行原則是對(duì)于計(jì)算量均勻分布，同時(shí)線程安全的格點(diǎn)計(jì)算使用二級(jí)分區(qū)并行，二級(jí)分區(qū)并行使用一維緯向二級(jí)分區(qū)劃分。對(duì)于計(jì)算量不均勻的格點(diǎn)計(jì)算和程序內(nèi)部線程不安全或存在MPI通信，則選擇循環(huán)并行方法。動(dòng)力框架部分線性非線性項(xiàng)計(jì)算、插值計(jì)算可使用二級(jí)分區(qū)并行，也可使用循環(huán)并行，Helmholtz方程求解和拉格朗日上游點(diǎn)的確定及上游點(diǎn)變量插值由于內(nèi)部有MPI通信，使用循環(huán)并行。物理過(guò)程中微物理過(guò)程、長(zhǎng)波輻射、短波輻射、陸面過(guò)程、地形重力波拖曳過(guò)程以及積云對(duì)流過(guò)程為循環(huán)并行。

表4 二級(jí)分區(qū)并行和循環(huán)并行結(jié)果對(duì)比Table 4 The comparison of tile-level and loop-level parallelization results

GRAPES模式混合并行原則適用于積分計(jì)算中大部分的子程序，但仍有一些子程序由于其計(jì)算的特殊性，需要特殊處理。下面分別介紹Helmholtz方程求解和有負(fù)載平衡問(wèn)題的短波輻射過(guò)程的并行處理方法。

3 Helmholtz方程求解混合并行分析

3.1 GCR（generalized conjugate residual）算法

GRAPES模式動(dòng)力框架的主要問(wèn)題之一就是如何有效求解Helmholtz方程［23］。求解Helmholtz方程就是求解形如Ax＝b的方程組，其中系數(shù)矩陣A為一大型稀疏矩陣。求解Helmholtz方程有很多方法，目前GRAPES模式使用GCR算法迭代求解Helmholtz方程。GCR算法計(jì)算流程參見(jiàn)文獻(xiàn)［18］。

下面測(cè)試計(jì)算核為64時(shí)不同線程數(shù)GCR算法混合并行效果，為了對(duì)比合理性，除了線程數(shù)為8時(shí)使用8個(gè)節(jié)點(diǎn)外，其余方案都使用16個(gè)節(jié)點(diǎn)。表5為各線程的GCR算法單個(gè)積分步混合并行效果對(duì)比，可以看出，GCR算法各線程混合并行計(jì)算時(shí)間和單一MPI方案差異不大，線程數(shù)為2時(shí)，GCR算法混合并行效果甚至好于單一MPI方案。

表5 GCR算法混合并行結(jié)果對(duì)比Table 5 The comparison of hybrid parallel results of GCR algorithm

3.2 帶有ILU預(yù)條件子的GCR算法

單一GCR算法收斂速度慢，為了加快收斂速度，GRAPES采用了不完全LU分解（incomplete LU factorization，ILU）預(yù)條件子加速GCR迭代的收斂過(guò)程，稱(chēng)為帶有ILU預(yù)條件子的GCR算法。但I(xiàn)LU方法的計(jì)算相關(guān)性復(fù)雜，并行化程度低，為了減少通信GRAPES模式采用分區(qū)域的局地ILU預(yù)條件子［24］。其分區(qū)域與模式子區(qū)域一級(jí)分區(qū)相同，稱(chēng)為patch-ILU預(yù)條件子，相應(yīng)的全區(qū)域的稱(chēng)為domain-ILU預(yù)條件子。局地預(yù)條件可保持計(jì)算在各子區(qū)域的獨(dú)立性，實(shí)現(xiàn)并行處理，但會(huì)降低預(yù)條件子的性能，迭代次數(shù)略有增加。

實(shí)現(xiàn)混合并行需要實(shí)現(xiàn)patch-ILU的Open MP并行處理。設(shè)計(jì)實(shí)現(xiàn)了兩種方法：一種是根據(jù)ILU計(jì)算方法的數(shù)據(jù)相關(guān)性，用Open MP直接實(shí)現(xiàn)并行處理，稱(chēng)為loop-patch-ILU；另一種方案是對(duì)一級(jí)分區(qū)分區(qū)，采用在更小的二級(jí)分區(qū)區(qū)域內(nèi)的局地預(yù)條件子，即tile-ILU預(yù)條件子，從而實(shí)現(xiàn)混合并行。本試驗(yàn)中tile-ILU預(yù)條件子的GCR算法收斂速度與patch-ILU預(yù)條件子的GCR算法大致相同，其中tile-ILU預(yù)條件子的GCR算法積分36步共迭代1823步，patch-ILU預(yù)條件子的GCR算法積分36步共迭代1826步。圖3是比較以上3種預(yù)條件并行處理方式下單一積分步GCR算法平均使用時(shí)間。由圖3可以看出，使用tile-ILU和loop-patch-ILU預(yù)條件子的GCR方法計(jì)算時(shí)間比patch-ILU的GCR方法少。使用tile-ILU預(yù)條件子的GCR方法在4個(gè)線程數(shù)內(nèi)混合并行效果略微好于loop-patch-ILU預(yù)條件子的GCR方法，線程數(shù)為8時(shí)，混合并行效果明顯好于loop-patch-ILU預(yù)條件子的GCR方法混合并行。

圖3 單一積分步內(nèi)帶有ILU預(yù)條件子的GCR算法平均時(shí)間Fig.3 The average time of GCR algorithm with ILU preconditioner in a single integrate step

4 短波輻射混合并行分析

短波輻射過(guò)程使用RRTMG（Rapid Radiative Transfer Model for GCMs）方案，計(jì)算采用氣柱模型，氣柱之間計(jì)算相互獨(dú)立，所以短波輻射可以使用二級(jí)分區(qū)并行。但由于短波輻射只有在有太陽(yáng)輻射的半球計(jì)算，從而導(dǎo)致負(fù)載不平衡。采用二級(jí)分區(qū)并行不能實(shí)現(xiàn)負(fù)載平衡。針對(duì)短波輻射的負(fù)載不平衡問(wèn)題，GRAPES模式首先通過(guò)根進(jìn)程將收集到的每個(gè)進(jìn)程中天頂角大于零度，即要做短波輻射的氣柱數(shù)由西向東、由南向北排成一維序列，并將其平均分配給每個(gè)進(jìn)程，通過(guò)通信按照新的分配結(jié)果實(shí)現(xiàn)數(shù)據(jù)重新分布，計(jì)算完成后將結(jié)果傳回原來(lái)進(jìn)程。短波輻射混合并行是將平均分配給每個(gè)進(jìn)程的氣柱進(jìn)而分配給線程，這樣短波輻射在每個(gè)線程中的計(jì)算量基本相同。該方法能夠有效解決負(fù)載平衡問(wèn)題。采用負(fù)載平衡的短波輻射過(guò)程64個(gè)進(jìn)程時(shí)通信時(shí)間為0.095 s，通信占短波輻射總時(shí)間的11.43%。

圖4為短波輻射使用無(wú)負(fù)載平衡與有負(fù)載平衡的效果對(duì)比，同時(shí)比較了不同線程調(diào)度方式的并行效果?？梢钥闯?，由于采用了負(fù)載重分配技術(shù)，循環(huán)并行的短波輻射計(jì)算時(shí)間遠(yuǎn)遠(yuǎn)小于無(wú)負(fù)載平衡的二級(jí)分區(qū)并行。比較循環(huán)并行線程不同調(diào)度方式發(fā)現(xiàn)，靜態(tài)調(diào)度方式和指導(dǎo)調(diào)度比動(dòng)態(tài)調(diào)度方式好。

圖4 短波輻射使用多種并行方法時(shí)間對(duì)比Fig.4 The comparsion of different parallel scheme computional time

5 試驗(yàn)結(jié)果和分析

根據(jù)上述方法對(duì)GRAPES模式進(jìn)行混合并行改造，以下是在中國(guó)氣象局的IBM計(jì)算機(jī)上的試驗(yàn)結(jié)果。

5.1 Open MP擴(kuò)展性試驗(yàn)

首先測(cè)試GRAPES模式混合并行方案中時(shí)間積分核心計(jì)算程序在MPI進(jìn)程數(shù)一定的情況下加速比隨Open MP線程數(shù)增加的情況。圖5顯示了GRAPES模式混合并行方案積分計(jì)算程序在進(jìn)程數(shù)為32，64，128時(shí)，線程數(shù)從1增加到8時(shí)的加速比情況。由圖5可以看出，進(jìn)程數(shù)相同時(shí)，隨著線程數(shù)的增加，GRAPES模式混合并行方案加速比增加。線程數(shù)為2時(shí)與理論加速比非常接近，但隨著線程數(shù)的增多，與理論加速比差距加大。線程數(shù)量為8時(shí)僅加速4倍左右。所以GRAPES模式混合并行方案中開(kāi)啟線程數(shù)不宜太多，2個(gè)線程到4個(gè)線程為宜。

圖5 積分計(jì)算程序混合并行加速比情況Fig.5 The speedup of integral computation in hybrid parallelization

5.2 積分計(jì)算主要子程序混合并行試驗(yàn)

下面討論積分計(jì)算程序內(nèi)部主要子程序混合并行效果。測(cè)試在計(jì)算核數(shù)相同情況下，單一MPI方案和混合并行方案并行效果。試驗(yàn)使用64個(gè)計(jì)算核，選取5個(gè)積分計(jì)算程序中計(jì)算時(shí)間最長(zhǎng)的子程序進(jìn)行分析。圖6為5個(gè)子程序在5種試驗(yàn)方案下單步運(yùn)行時(shí)間統(tǒng)計(jì)。由圖5可以看出，陸面過(guò)程、長(zhǎng)波輻射過(guò)程以及微物理過(guò)程混合并行效果好于單一MPI方案。Helmholtz求解子程序和負(fù)載平衡的短波輻射計(jì)算時(shí)間比單一MPI方案多一些。

圖6 積分計(jì)算中主要子程序不同試驗(yàn)方案計(jì)算時(shí)間對(duì)比Fig.6 The comparison of main subroutine integral computation time in each experiment scheme

5.3 積分計(jì)算混合并行整體效果及擴(kuò)展性試驗(yàn)

下面測(cè)試GRAPES模式混合并行方案在大規(guī)模并行下的積分整體效果。為了排除積分第1步數(shù)據(jù)初始化干擾，從第2步積分開(kāi)始計(jì)時(shí)。圖7是IBM系統(tǒng)使用不同計(jì)算核數(shù)GRAPES模式積分35步的計(jì)算時(shí)間。試驗(yàn)測(cè)試到4096核，當(dāng)計(jì)算規(guī)模達(dá)到4096核時(shí)，GRAPES模式單一MPI方案不能運(yùn)行下去，而線程數(shù)為2，4，8時(shí)的混合并行方案仍可運(yùn)行。由圖7可以看出，在計(jì)算核數(shù)一定情況下，線程數(shù)為2時(shí)的混合并行方案計(jì)算時(shí)間比單一MPI方案少。計(jì)算核數(shù)達(dá)到1024及以上時(shí)，4個(gè)線程的混合并行方案計(jì)算時(shí)間少于單一MPI方案。線程數(shù)量大于4時(shí)，混合并行效果顯著下降。

GRAPES模式單一MPI方案在MPI進(jìn)程數(shù)達(dá)到一定規(guī)模后不能運(yùn)行。在IBM系統(tǒng)上，當(dāng)MPI進(jìn)程數(shù)為4096時(shí)，分辨率為1°×1°的GRPAES模式不能運(yùn)行。但GRAPES模式混合并行方案在計(jì)算核數(shù)為4096時(shí)仍能運(yùn)行下去，積分35步計(jì)算時(shí)間比計(jì)算核數(shù)2048時(shí)少。可以看出，GRAPES模式單一MPI方案在達(dá)到擴(kuò)展性限制時(shí)，混合并行方案可以獲得比單一MPI方案更好的擴(kuò)展性。

圖7 不同計(jì)算核數(shù)積分35步計(jì)算時(shí)間情況Fig.7 The integral time of 35 steps with different computing cores

6 結(jié) 論

本試驗(yàn)采用MPI與Open MP混合并行方案，具體采用區(qū)域分解并行和循環(huán)并行兩種方法，對(duì)GRAPES全球模式中計(jì)算時(shí)間最長(zhǎng)的積分部分進(jìn)行混合并行改造，并進(jìn)行初步并行試驗(yàn)，得到以下結(jié)論：

1）GRAPES模式混合并行方案可以適應(yīng)不同計(jì)算架構(gòu)平臺(tái)，可在分布式內(nèi)存系統(tǒng)運(yùn)行，也可運(yùn)行于多核處理器集群系統(tǒng)。

2）在IBM系統(tǒng)上GRAPES模式大規(guī)模運(yùn)行時(shí)，4個(gè)線程內(nèi)的混合并行方案效果好于單一MPI方案，線程數(shù)量大于4時(shí)，混合并行效果顯著下降。當(dāng)GRAPES模式單一MPI方案因擴(kuò)展性限制不能運(yùn)行時(shí)，混合并行方案可以替代它，獲得比單一MPI方案更好的擴(kuò)展性。

3）區(qū)域分解并行僅適合計(jì)算量均勻分布且線程安全的格點(diǎn)計(jì)算，并使用一維緯向二級(jí)分區(qū)劃分，其效果與循環(huán)并行相仿或略差。對(duì)于計(jì)算量不均勻的格點(diǎn)計(jì)算和程序內(nèi)部線程不安全的情況，應(yīng)選擇循環(huán)并行。

4）求解Helmholtz方程的GCR算法本身適合混合并行，但I(xiàn)LU預(yù)條件子不易并行計(jì)算，使用二級(jí)分區(qū)區(qū)域局部預(yù)條件子可有效實(shí)現(xiàn)預(yù)條件子的并行處理。求解Helmholtz算法總體混合并行效果比單一MPI方法略差。

5）負(fù)載平衡對(duì)于短波輻射計(jì)算十分關(guān)鍵。在計(jì)算核數(shù)相同情況下，負(fù)載平衡的短波輻射循環(huán)并行效果接近于單一MPI方案。

6）GRAPES模式中使用循環(huán)并行的長(zhǎng)波輻射方案、微物理過(guò)程、陸面模式過(guò)程并行效果優(yōu)于單一MPI方案。

后續(xù)工作包括選取適合混合并行的求解Helmholtz方程方法，對(duì)混合并行方案進(jìn)一步優(yōu)化，并在計(jì)算資源條件允許下試驗(yàn)高分辨率GRAPES模式混合并行方案在大規(guī)模并行下的運(yùn)行效果。

［1］ Gysi T，F(xiàn)uhrer O，Osuna C，et al.Porting COSMO to Hybrid Architectures.［2013-04-14］.http：∥data1.gfdl.noaa.gov／multi-core／2012／presentations／Session_2_Messmer.pdf.

［2］ 馮云，周淑秋.MPI＋Open MP混合并行編程模型應(yīng)用研究.計(jì)算機(jī)系統(tǒng)應(yīng)用，2006（2）：33-35.

［3］ 樊志杰，趙文濤.GRAPES四維變分同化系統(tǒng)MPI和Open MP混合算法研究.計(jì)算機(jī)光盤(pán)軟件與應(yīng)用，2012（19）：21-23.

［4］ The Weather Research and Forecasting（WRF）Model.［2013-01-09］.http：∥wrf-model.org／.

［5］ The Users Home Page for the Weather Research and Forecasting（WRF）Modeling System.［2013-01-09］.http：∥www.mmm.ucar.edu／wrf／users／.

［6］ Skamarock W C，Klemp J B，Dudhia J，et al.A Description ofthe Advanced Research WRF Version 3.NCAR Tech Note NCAR／TN-475＋STR，2005.

［7］ ?ipkováV，Lúcny A，Gazák M.Experiments with a Hybrid-Parallel Model of Weather Research and Forecasting（WRF）System.GCCP 2010 Book of Abstracts ，2010：37.

［8］ Epicoco I，Mocavero S，Giovanni A.NEMO-Med：Optimization and Improvement of Scalability.CMCC Research Paper，2011.

［9］ 張昕，季仲貞，王斌.Open MP在 MM5中尺度模式中的應(yīng)用試驗(yàn).氣候與環(huán)境研究，2001，6（1）：84-90.

［10］ 朱政慧，施培量，顏宏.用 Open MP并行化氣象預(yù)報(bào)模式試驗(yàn).應(yīng)用氣象學(xué)報(bào)，2002，13（1）：102-108.

［11］ 朱政慧.并行高分辨率有限區(qū)預(yù)報(bào)系統(tǒng)在IBM SP上的建立.應(yīng)用氣象學(xué)報(bào)，2003，14（1）：119-121.

［12］ 朱政慧.一個(gè)數(shù)值天氣預(yù)報(bào)模式的并行混合編程模型及其應(yīng)用.數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用，2005，26（3）：203-204.

［13］ 郭妙，金之雁，周斌.基于通用圖形處理器的GRAPES長(zhǎng)波輻射并行方案.應(yīng)用氣象學(xué)報(bào)，2012，23（3）：348-354.

［14］ 鄭芳，許先斌，向冬冬，等.基于GPU的GRAPES數(shù)值預(yù)報(bào)系統(tǒng)中RRTM模塊的并行化研究.計(jì)算機(jī)科學(xué)，2012，39（6）：370-374.

［15］ Open MP Specications.Open MP Application Programing Interface.V3.0，2008.［2013-01-09］.http：∥www.openmp.org／mp-documents／spec30.pdf.

［16］ Chapman B，Jost G，Van Der Pas R.Using Open MP：Portable Shared Memory Parallel Programming.London：MIT Press，2008.

［17］ Blaise Barney.Open MP.［2013-01-09］.https：∥computing.llnl.gov／tutorials／open MP／.

［18］ 薛紀(jì)善，陳德輝.數(shù)值預(yù)報(bào)系統(tǒng)GRAPES的科學(xué)設(shè)計(jì)與應(yīng)用.北京：科學(xué)出版社，2008.

［19］ 伍湘君.GRAPES高分辨率氣象數(shù)值預(yù)報(bào)模式并行計(jì)算關(guān)鍵技術(shù)研究.北京：國(guó)防科學(xué)技術(shù)大學(xué)，2011.

［20］ 伍湘君，金之雁，黃麗萍，等.GRAPES模式軟件框架與實(shí)現(xiàn).應(yīng)用氣象學(xué)報(bào)，2005，16（4）：539-546.

［21］ Fowler R F，Greenough C.Mixed MPI：Open MP Programming：A Study in Parallelisation of a CFD Multiblock Code.CCLRC Rutherford Appleton Laboratory，2003.

［22］ 金之雁，王鼎興.一種在異構(gòu)系統(tǒng)中實(shí)現(xiàn)負(fù)載平衡的方法.應(yīng)用氣象學(xué)報(bào)，2003，14（4）：410-418.

［23］ 陳德輝，沈?qū)W順.新一代數(shù)值預(yù)報(bào)系統(tǒng)GRAPES研究進(jìn)展.應(yīng)用氣象學(xué)報(bào)，2007，17（6）：773-777.

［24］ 劉宇，曹建文.適用于GRAPES數(shù)值天氣預(yù)報(bào)軟件的ILU預(yù)條件子.計(jì)算機(jī)工程與設(shè)計(jì)，2008，29（3）：731-734.

The Hybrid MPI and Open MP Parallel Scheme of GRAPES_Global Model

Jiang Qingu1）Jin Zhiyan2）

1）（Chinese Academy of Meteorological Sciences，Beijing100081）
2）（Center for Numerical Weather Prediction，CMA，Beijing100081）

Clustered SMP systems are gradually becoming more prominence，as advances in multi-core technology which allows larger numbers of CPUs to have access to a single memory space.To take advantage of benefits of this hardware architecture that combines both distributed and shared memory，utilizing hybrid MPI and Open MP parallel programming model is a good trial.This hierarchical programming model can achieve both inter-node and intra-node parallelization by using a combination of message passing and thread based shared memory parallelization paradigms within the same application.MPI is used to coarse-grained communicate between SMP nodes and Open MP based on threads is used to fine-grained compute within a SMP node.

As a large-scale computing and storage-intensive typical numerical weather forecasting application，GRAPES（Global／Regional Assimilation and Pr Edictions System）has been developed into MPI version and put into operational use.To adapt to SMP cluster systems and achieve higher scalability，a hybrid MPIand Open MP parallel model suitable for GRPAES_Global model is developed with the introduction of horizontal domain decomposition method and loop-level parallelization.In horizontal domain decomposition method，a patch is uniformly divided into several tiles while patches are obtained by dividing the whole forecasting domain.There are two main advantages in performing parallel operations on tiles.Firstly，tilelevel parallelization which applies Open MP at a high level，to some extent，is coarse grained parallelism.Compared to computing work associated with each tile，Open MP thread overhead is negligible.Secondly，implementation of this method is relative simple，and the subroutine thread safety is the only thing to ensure.Loop-level parallelization which can improve load imbalance by adopting different thread scheduling policies is fine grained parallelism.The main computational loops are applied Open MP’s parallel directives in loop-level parallelization method.The preferred method is horizontal domain decomposition for uniform grid computing，while loop-level parallelization method is preferred for non-uniform grid computing and the thread unsafe procedures.Experiments with 1°×1°dataset are performed and timing on main subroutines of integral computation are compared.The hybrid parallel performance is superior to single MPI scheme in terms of long wave radiation process，microphysics and land surface process while Helmholtz equation generalized conjugate residual（GCR）solution has some difficulty in thread parallelism for incomplete LU factorization preconditioner part.ILU part with tile-level parallelization can improve GCR’s hybrid parallelization.Short wave process hybrid parallel performance is close to single MPI scheme under the same computing cores.It requires less elapsed time with increase of the number of threads under certain MPI processes in hybrid parallel scheme.And hybrid parallel scheme within four threads is superior to single MPI scheme under large-scale experiment.Hybrid parallel scheme can also achieve better scalability than single MPI scheme.The experiment shows hybrid MPI and Open MP parallel scheme is suitable for GRAPES_Global model.

hybrid parallel；numerical weather forecasting model；domain decomposition；loop-level parallelization

蔣沁谷，金之雁.GRAPES全球模式MPI與Open MP混合并行方案.應(yīng)用氣象學(xué)報(bào)，2014，25（5）：581-591.

2013-10-25收到，2014-04-30收到再改稿。

國(guó)家自然科學(xué)基金項(xiàng)目（61361120098）

*通信作者，email：jinzy＠cma.gov.cn