蔣沁谷金之雁

1）（中國氣象科學研究院，北京100081）2）（中國氣象局數(shù)值預報中心，北京100081）

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

蔣沁谷1）金之雁2）*

1）（中國氣象科學研究院，北京100081）2）（中國氣象局數(shù)值預報中心，北京100081）

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

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

引 言

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

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

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

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

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

1 Open MP編程模型

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

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

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

2 GRAPES混合并行方案設計

2.1 MPI與Open MP混合并行模型

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

2.2 GRAPES混合并行方案

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

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

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

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

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

2.2.1 二級分區(qū)并行

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

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

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

2.2.2 循環(huán)并行

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 GCR（generalized conjugate residual）算法

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

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

表5 GCR算法混合并行結果對比Table 5 The comparison of hybrid parallel results of GCR algorithm

3.2 帶有ILU預條件子的GCR算法

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

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

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

4 短波輻射混合并行分析

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

圖4為短波輻射使用無負載平衡與有負載平衡的效果對比，同時比較了不同線程調度方式的并行效果。可以看出，由于采用了負載重分配技術，循環(huán)并行的短波輻射計算時間遠遠小于無負載平衡的二級分區(qū)并行。比較循環(huán)并行線程不同調度方式發(fā)現(xiàn)，靜態(tài)調度方式和指導調度比動態(tài)調度方式好。

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

5 試驗結果和分析

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

5.1 Open MP擴展性試驗

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

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

5.2 積分計算主要子程序混合并行試驗

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

圖6 積分計算中主要子程序不同試驗方案計算時間對比Fig.6 The comparison of main subroutine integral computation time in each experiment scheme

5.3 積分計算混合并行整體效果及擴展性試驗

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

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

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

6 結 論

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

1）GRAPES模式混合并行方案可以適應不同計算架構平臺，可在分布式內存系統(tǒng)運行，也可運行于多核處理器集群系統(tǒng)。

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

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

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

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

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

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

［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混合并行編程模型應用研究.計算機系統(tǒng)應用，2006（2）：33-35.

［3］ 樊志杰，趙文濤.GRAPES四維變分同化系統(tǒng)MPI和Open MP混合算法研究.計算機光盤軟件與應用，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中尺度模式中的應用試驗.氣候與環(huán)境研究，2001，6（1）：84-90.

［10］ 朱政慧，施培量，顏宏.用 Open MP并行化氣象預報模式試驗.應用氣象學報，2002，13（1）：102-108.

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

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

［13］ 郭妙，金之雁，周斌.基于通用圖形處理器的GRAPES長波輻射并行方案.應用氣象學報，2012，23（3）：348-354.

［14］ 鄭芳，許先斌，向冬冬，等.基于GPU的GRAPES數(shù)值預報系統(tǒng)中RRTM模塊的并行化研究.計算機科學，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］ 薛紀善，陳德輝.數(shù)值預報系統(tǒng)GRAPES的科學設計與應用.北京：科學出版社，2008.

［19］ 伍湘君.GRAPES高分辨率氣象數(shù)值預報模式并行計算關鍵技術研究.北京：國防科學技術大學，2011.

［20］ 伍湘君，金之雁，黃麗萍，等.GRAPES模式軟件框架與實現(xiàn).應用氣象學報，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］ 金之雁，王鼎興.一種在異構系統(tǒng)中實現(xiàn)負載平衡的方法.應用氣象學報，2003，14（4）：410-418.

［23］ 陳德輝，沈學順.新一代數(shù)值預報系統(tǒng)GRAPES研究進展.應用氣象學報，2007，17（6）：773-777.

［24］ 劉宇，曹建文.適用于GRAPES數(shù)值天氣預報軟件的ILU預條件子.計算機工程與設計，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混合并行方案.應用氣象學報，2014，25（5）：581-591.

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

國家自然科學基金項目（61361120098）

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