POM海洋模式的并行算法

李冬，劉璟，韓桂軍，張學(xué)峰，王喜冬

（1. 南開(kāi)大學(xué)信息技術(shù)科學(xué)學(xué)院，計(jì)算機(jī)科學(xué)與技術(shù)系，天津，300071；2. 國(guó)家海洋局國(guó)家海洋信息中心，天津，300171）

POM海洋模式的并行算法

李冬1,2，劉璟1，韓桂軍2，張學(xué)峰2，王喜冬2

（1. 南開(kāi)大學(xué)信息技術(shù)科學(xué)學(xué)院，計(jì)算機(jī)科學(xué)與技術(shù)系，天津，300071；2. 國(guó)家海洋局國(guó)家海洋信息中心，天津，300171）

POM模式目前尚無(wú)正式發(fā)布的并行版本。通過(guò)對(duì)POM串行程序的數(shù)據(jù)流向分析，討論了POM模式并行化所涉及的關(guān)鍵算法和主要技術(shù)問(wèn)題；并基于消息傳遞接口（MPI），研發(fā)了POM模式的并行版本。測(cè)試結(jié)果表明，POM并行軟件效率較高，達(dá)到了業(yè)務(wù)化要求，業(yè)已應(yīng)用于國(guó)家海洋信息中心的再分析業(yè)務(wù)化系統(tǒng)中。

POM模式；并行算法；MPI；再分析

POM (Princeton Ocean Model) 模式[1]是美國(guó)普林斯頓大學(xué)發(fā)展的三維斜壓原始方程海洋模式。國(guó)內(nèi)外許多業(yè)務(wù)化數(shù)值預(yù)報(bào)和再分析系統(tǒng)都是以POM模式為基礎(chǔ)開(kāi)發(fā)的。

近年來(lái)，隨著業(yè)務(wù)化系統(tǒng)的需要及計(jì)算條件的逐步提高，海洋模式的時(shí)空分辨率越來(lái)越高。我們多年的研究實(shí)踐表明，采用POM模式的串行版軟件，能夠較好的滿足潮汐和三維潮流的數(shù)值計(jì)算時(shí)效要求。但若研發(fā)再分析、數(shù)據(jù)同化和數(shù)值預(yù)報(bào)系統(tǒng)，無(wú)論從再分析業(yè)務(wù)化運(yùn)行時(shí)間，數(shù)值預(yù)報(bào)的時(shí)效、還是從模式的深入開(kāi)發(fā)調(diào)試上講，POM模式的現(xiàn)有串行版軟件均難以滿足業(yè)務(wù)化系統(tǒng)的要求。

遺憾的是，POM模式目前尚無(wú)正式發(fā)布的并行版本。雖然Internet上可以下載到并行的Cousins版本[2](該版本采用TOPAZ軟件開(kāi)發(fā)了MP-POM（Massively Parallel-POM，并已移植到Cray T3E和SGI Origin 2000），但用戶的使用情況表明，其計(jì)算結(jié)果與相應(yīng)的串行版本的計(jì)算結(jié)果存在明顯差異，不能滿足實(shí)際應(yīng)用。同時(shí)，三維海洋模式的并行化是極具挑戰(zhàn)性的工作，技術(shù)難度大，復(fù)雜程度高。所有這些，都給使用POM模式的海洋科研工作者帶來(lái)極大的不便，阻礙了其業(yè)務(wù)化系統(tǒng)的建設(shè)和應(yīng)用。

本文將討論P(yáng)OM模式并行化過(guò)程中涉及的關(guān)鍵算法和主要技術(shù)問(wèn)題，在對(duì)POM模式串行程序進(jìn)行數(shù)據(jù)流向分析的基礎(chǔ)上，基于消息傳遞接口（MPI）[3]，依托自主開(kāi)發(fā)的基礎(chǔ)通信模塊，設(shè)計(jì)并完成了POM模式并行軟件的開(kāi)發(fā)。

文章第二部分對(duì)POM海洋模式進(jìn)行簡(jiǎn)單介紹；第三部分詳細(xì)討論P(yáng)OM模式的并行算法、并行程序框架及通信軟件模塊；第四部分介紹并行程序的測(cè)試結(jié)果；第五部分給出結(jié)論。

2 POM海洋模式簡(jiǎn)介

POM模式的主要特點(diǎn)為：

(1) 垂向混合系數(shù)由二階湍流閉合模式確定。

(2) 垂向采用σ坐標(biāo)。

(3) 水平采用正交曲線網(wǎng)格及“Arakawa C”網(wǎng)格。

(4) 水平差分格式采用顯格式，垂向差分格式采用隱格式。

(5) 采用自由表面及時(shí)間分裂格式；外模式(正壓模)為二維，采用較短的時(shí)間步長(zhǎng)（由CFL條件及外重力波波速?zèng)Q定）；內(nèi)模式(斜壓模)為三維，采用較長(zhǎng)的時(shí)間步長(zhǎng)（由CFL條件及內(nèi)波波速?zèng)Q定）。

(6) 模式包含完整的熱力學(xué)過(guò)程。

3 POM模式并行算法及軟件設(shè)計(jì)

3.1 算法及程序框架

由于POM模式的差分格式，水平方向采用顯格式，z方向采用隱格式，因此，在水平方向各格點(diǎn)的計(jì)算僅與相鄰格點(diǎn)的前面時(shí)刻的數(shù)據(jù)相關(guān)。設(shè)計(jì)POM模式并行化的關(guān)鍵是如何有效地進(jìn)行重疊區(qū)數(shù)據(jù)的更新。

重疊區(qū)數(shù)據(jù)的更新需要通過(guò)與周圍相鄰進(jìn)程的通信完成，但若完全依靠通信機(jī)制來(lái)實(shí)現(xiàn)，通信開(kāi)銷太大，嚴(yán)重影響程序的并行效率。我們?cè)趯?duì)POM模式串行程序進(jìn)行數(shù)據(jù)流分析的基礎(chǔ)上，將相鄰進(jìn)程通信和本進(jìn)程的冗余計(jì)算相結(jié)合來(lái)實(shí)現(xiàn)重疊區(qū)數(shù)據(jù)的更新，并采用MPI提供的派生數(shù)據(jù)類型將多變量數(shù)據(jù)打包在一起更新，極大地降低了通信的開(kāi)銷。簡(jiǎn)單起見(jiàn)，以一維情形為例，簡(jiǎn)要說(shuō)明并行程序的設(shè)計(jì)思路。假設(shè)有如下簡(jiǎn)單的串行代碼：

……（此處省略部分代碼）

由于計(jì)算y(i)時(shí)要用到x(i-1)值，計(jì)算z(i)時(shí)用到y(tǒng)(i+1)和y(i-1)，因此改寫(xiě)成并行程序時(shí)，每個(gè)進(jìn)程要進(jìn)行兩次重疊區(qū)的更新操作（在計(jì)算數(shù)組y、z之前各進(jìn)行一次）。

而通過(guò)數(shù)據(jù)流相關(guān)分析可知（圖1），事實(shí)上，每個(gè)進(jìn)程在計(jì)算數(shù)組y之前，只要更新自己重疊區(qū)中的x數(shù)組的3個(gè)元素（在圖中以圓圈標(biāo)記），之后計(jì)算數(shù)組y和z的工作就完全在本進(jìn)程完成，而無(wú)需與其它進(jìn)程通信。與完全通過(guò)通信機(jī)制進(jìn)行重疊區(qū)的更新相比，這雖然略微增加了本地的冗余計(jì)算（計(jì)算y時(shí)有冗余），但卻減少了一次重疊區(qū)更新操作。而通過(guò)通信進(jìn)行的重疊區(qū)更新具有很大的通信開(kāi)銷（算法見(jiàn)3.2.2），在實(shí)際程序中，當(dāng)大量這樣的更新操作被減少時(shí)，并行程序的效率會(huì)得到顯著提高。

圖1 本地冗余計(jì)算與進(jìn)程通信結(jié)合實(shí)現(xiàn)重疊區(qū)更新示意圖Fig. 1 Halo regions updated by the combination of processes communicating and local computing

我們將這樣的一段代碼稱之為一個(gè)“并行區(qū)域”，即通信主要集中于此代碼段的入口，而在代碼段內(nèi)部主要是本地計(jì)算，幾乎不進(jìn)行（或僅有少量）進(jìn)程之間的通信。采用本地冗余計(jì)算與進(jìn)程通信相結(jié)合的方式進(jìn)行重疊區(qū)的更新，其帶來(lái)的另一個(gè)優(yōu)點(diǎn)是：由于各變量均于“并行區(qū)域”的入口處進(jìn)行重疊區(qū)的更新，因此，可以通過(guò)派生數(shù)據(jù)類型將這些變量打包一并更新，再次顯著地減少了并行程序的通信次數(shù)，極大地提高了并行效率。

我們?cè)趯?duì)POM模式串行程序進(jìn)行數(shù)據(jù)流向分析的基礎(chǔ)上，采用本地冗余計(jì)算與通信結(jié)合更新重疊區(qū)的方式，基于MPI開(kāi)發(fā)了POM并行程序，基本流程見(jiàn)圖2。

其中，在內(nèi)（斜壓）、外（正壓）模式之間建立了3個(gè)“并行區(qū)域”，在這些區(qū)域中，僅有少量的進(jìn)程通信，大部分通信集中于它們的入口，通過(guò)采用派生數(shù)據(jù)類型，將需要消息傳輸?shù)淖兞恳黄鸶?，顯著減少了通信次數(shù)，提高了并行效率。

圖2 POM并行程序的基本流程圖Fig. 2 Flow chart of parallel POM

3.2 通信模塊的設(shè)計(jì)和開(kāi)發(fā)

正確、高效地實(shí)現(xiàn)不同進(jìn)程之間的通信，是并行程序設(shè)計(jì)要解決的關(guān)鍵問(wèn)題之一。我們對(duì)MPI標(biāo)準(zhǔn)庫(kù)函數(shù)進(jìn)行了封裝，開(kāi)發(fā)了一系列高層的基礎(chǔ)通信模塊，并依托這些基礎(chǔ)通信模塊進(jìn)行POM模式的并行開(kāi)發(fā)工作。

這些基礎(chǔ)通信模塊包括：區(qū)域劃分及進(jìn)程的虛擬拓?fù)浣Y(jié)構(gòu)模塊、數(shù)據(jù)重疊區(qū)更新模塊、派生數(shù)據(jù)類型模塊、數(shù)據(jù)合并及分解模塊、I/O接口模塊等。

3.2.1 區(qū)域劃分及進(jìn)程的虛擬拓?fù)浣Y(jié)構(gòu)模塊根據(jù)對(duì)POM數(shù)值模式的差分格式（水平方向采用顯格式，z方向采用隱格式）及其串行程序的數(shù)據(jù)流向分析，POM并行軟件將計(jì)算海區(qū)按經(jīng)緯度劃分成多個(gè)子海區(qū)分配給不同的進(jìn)程。

進(jìn)行區(qū)域劃分以后，設(shè)計(jì)進(jìn)程的虛擬拓?fù)浣Y(jié)構(gòu)，對(duì)進(jìn)程進(jìn)行編號(hào)，建立數(shù)據(jù)和任務(wù)與進(jìn)程之間的映射；同時(shí)，為了各進(jìn)程之間通信的方便，建立進(jìn)程之間的聯(lián)系，即每個(gè)進(jìn)程除了記錄自己的信息之外，還要記錄其周圍相鄰進(jìn)程的基本信息。圖3是采用16個(gè)進(jìn)程時(shí)的進(jìn)程虛擬拓?fù)浣Y(jié)構(gòu)示意圖，圖中的數(shù)字表示進(jìn)程的編號(hào)，每個(gè)進(jìn)程只負(fù)責(zé)自己所管轄海區(qū)的計(jì)算任務(wù)。圖中的海區(qū)即為基于POM的中國(guó)海及鄰近海域預(yù)報(bào)模式的計(jì)算區(qū)域。

3.2.2 數(shù)據(jù)重疊區(qū)更新模塊 基于區(qū)域劃分的海洋模式并行算法，由于數(shù)值計(jì)算格式的需要，每個(gè)進(jìn)程在計(jì)算過(guò)程中不可避免地要用到其周圍相鄰進(jìn)程的數(shù)據(jù)，因此，通常需要在每個(gè)進(jìn)程所負(fù)責(zé)的計(jì)算區(qū)域周圍建立重疊區(qū)（Overlap Areas或Fake Zones、Halo Regions）。本地進(jìn)程重疊區(qū)數(shù)據(jù)的更新由其周圍相鄰進(jìn)程發(fā)送，本機(jī)接收。對(duì)重疊區(qū)數(shù)據(jù)進(jìn)行正確的更新是保證并行程序結(jié)果正確性的關(guān)鍵，重疊區(qū)數(shù)據(jù)更新效率的好壞也直接影響并行程序最終的并行效率。圖4為POM模式并行計(jì)算軟件采用的重疊區(qū)更新算法的流程示意圖：

（1）非阻塞接收相鄰的北部進(jìn)程發(fā)來(lái)的數(shù)據(jù)，以更新本進(jìn)程的重疊區(qū)1；非阻塞接收相鄰的南部進(jìn)程發(fā)來(lái)的數(shù)據(jù)，以更新本進(jìn)程的重疊區(qū)2。

（2）阻塞發(fā)送本進(jìn)程計(jì)算區(qū)域中的北部數(shù)據(jù)至相鄰的北部進(jìn)程的重疊區(qū)2；阻塞發(fā)送本進(jìn)程計(jì)算區(qū)域中的南部數(shù)據(jù)至相鄰的南部進(jìn)程的重疊區(qū)1。

（3）等待本進(jìn)程重疊區(qū)1和重疊區(qū)2數(shù)據(jù)接收的完成（這是必須的，否則，可能會(huì)導(dǎo)致東南、東北、西南、西北四個(gè)角處的重疊區(qū)數(shù)據(jù)的錯(cuò)誤更新）。但這里并不需要進(jìn)程之間的同步操作。

圖3 進(jìn)程虛擬拓?fù)浣Y(jié)構(gòu)示意圖Fig. 3 Virtual topologies of processes

（4）非阻塞接收相鄰的西部進(jìn)程發(fā)來(lái)的數(shù)據(jù)，以更新本進(jìn)程的重疊區(qū)3；非阻塞接收相鄰的東部進(jìn)程發(fā)來(lái)的數(shù)據(jù)，以更新本進(jìn)程的重疊區(qū)4。

（5）阻塞發(fā)送本進(jìn)程計(jì)算區(qū)域中的西部數(shù)據(jù)至相鄰的西部進(jìn)程的重疊區(qū)4；阻塞發(fā)送本進(jìn)程計(jì)算區(qū)域中的東部數(shù)據(jù)至相鄰的東部進(jìn)程的重疊區(qū)3。

（6）等待本進(jìn)程重疊區(qū)3和重疊區(qū)4數(shù)據(jù)接收的完成。

所有的進(jìn)程同時(shí)執(zhí)行這6個(gè)步驟，最后的結(jié)果是每個(gè)進(jìn)程的重疊區(qū)數(shù)據(jù)都得到了正確、有效的更新。

另外，重疊區(qū)更新模塊還包括：沿經(jīng)度方向進(jìn)行重疊區(qū)更新的模塊及沿緯度方向進(jìn)行重疊區(qū)更新的模塊。

圖4 重疊區(qū)更新算法流程示意圖Fig. 4 Updating algorithms for hallo regions

3.2.3 派生數(shù)據(jù)類型模塊 POM并行程序中經(jīng)常要發(fā)送和接收內(nèi)存中非連續(xù)分布的數(shù)據(jù)。如圖5所示，假如進(jìn)程要更新右側(cè)重疊區(qū)切塊中的海溫?cái)?shù)據(jù)(存放于數(shù)組T(i, j, k)中)，由于Fortran語(yǔ)言是列優(yōu)先的，因此，數(shù)組T(i, j, k)重疊區(qū)中同一層中不同行（如圖中標(biāo)記為1和2的區(qū)域）的元素在內(nèi)存中地址是不連續(xù)的，不同層（如圖中標(biāo)記為1、3或2、4的區(qū)域）的數(shù)據(jù)其地址就更不連續(xù)。此外，進(jìn)行重疊區(qū)更新時(shí)，通常需要更新多個(gè)數(shù)組變量，如海溫、鹽度、流速、水位、水深等，不同的變量在內(nèi)存中當(dāng)然也不是連續(xù)分布的，如果對(duì)這些變量分別發(fā)送和接收，則會(huì)使通信次數(shù)顯著增加，降低并行效率。

為此，POM并行程序定義了一系列MPI派生數(shù)據(jù)類型，用以描述程序中經(jīng)常進(jìn)行通信的非連續(xù)數(shù)據(jù)在內(nèi)存中的分布。通過(guò)使用這些派生數(shù)據(jù)類型，可以方便地對(duì)數(shù)據(jù)進(jìn)行抽取，合并和通信，并且可以將多個(gè)不同數(shù)據(jù)類型的變量一起發(fā)送和接收，提高程序的并行效率。此外，在并行I/O中，也經(jīng)常要利用派生數(shù)據(jù)類型進(jìn)行文件的讀寫(xiě)。

3.2.4 數(shù)據(jù)合并及分解模塊 設(shè)計(jì)POM模式并行程序時(shí)，需將計(jì)算區(qū)域進(jìn)行劃分，把數(shù)據(jù)和計(jì)算任務(wù)分配給不同的進(jìn)程。然而，在模式開(kāi)發(fā)及程序的調(diào)試過(guò)程中，開(kāi)發(fā)者經(jīng)常需要了解區(qū)域數(shù)據(jù)與全局?jǐn)?shù)據(jù)之間的關(guān)系；同時(shí)，模式的輸出也要求對(duì)各進(jìn)程的數(shù)據(jù)進(jìn)行合并。因此，設(shè)計(jì)了數(shù)據(jù)的合并和分解模塊，主要包括如下功能：

(1) 將各進(jìn)程上的局部數(shù)組數(shù)據(jù)合并為整個(gè)計(jì)算區(qū)域的整體數(shù)組數(shù)據(jù)。這通過(guò)對(duì)局部數(shù)組定義派生數(shù)據(jù)類型，并利用集合通信中的收集調(diào)用來(lái)實(shí)現(xiàn)。

(2) 將整個(gè)計(jì)算區(qū)域的整體數(shù)組數(shù)據(jù)分解為各個(gè)進(jìn)程上的局部數(shù)組數(shù)據(jù)。這通過(guò)對(duì)局部數(shù)組定義派生數(shù)據(jù)類型，并利用集合通信中的散發(fā)調(diào)用來(lái)實(shí)現(xiàn)。

(3) 點(diǎn)在整個(gè)計(jì)算區(qū)域中的全局坐標(biāo)轉(zhuǎn)換為其在相應(yīng)進(jìn)程中的局地坐標(biāo)，并返回此進(jìn)程號(hào)及其所管轄區(qū)域的起始點(diǎn)坐標(biāo)。

(4) 點(diǎn)在進(jìn)程中的局地位置坐標(biāo)轉(zhuǎn)換為其在整個(gè)計(jì)算區(qū)域中的全局坐標(biāo)。

3.2.5 并行I/O模塊 通過(guò)使用自定義派生數(shù)據(jù)類型、定義文件視圖并使用聚合I/O函數(shù)實(shí)現(xiàn)了POM模式的并行I/O接口模塊。

圖5 派生數(shù)據(jù)類型示意圖Fig. 5 Derived data type

4 測(cè)試結(jié)果

我們對(duì)中國(guó)海及鄰近海域的POMgcs模式（POM模式的變種，與POM的區(qū)別在于其垂向采用 及z坐標(biāo)的混合坐標(biāo)，但并行方案與POM模式完全相同）并行軟件在國(guó)家海洋信息中心的高性能機(jī)群系統(tǒng)上進(jìn)行了測(cè)試。模式的水平網(wǎng)格分辨率為285×307，垂向35層，模式積分時(shí)間為3天。表1為性能測(cè)試結(jié)果；圖6為加速比曲線，其中，虛線代表線性（理想）加速比，實(shí)線代表實(shí)際加速比。

測(cè)試結(jié)果表明，在并行效果方面，當(dāng)進(jìn)程數(shù)目不超過(guò)48個(gè)，POMgcs并行程序的并行效率基本在70%以上；當(dāng)進(jìn)程數(shù)目不超過(guò)20個(gè)時(shí)，POMgcs并行程序可以保持線性加速比，特別當(dāng)進(jìn)程數(shù)目在10個(gè)以內(nèi)時(shí)出現(xiàn)了超線性加速比現(xiàn)象，這可能是由于高速緩存(cache)的影響，操作系統(tǒng)開(kāi)銷的均攤等原因造成的。當(dāng)使用的進(jìn)程數(shù)目超過(guò)50個(gè)后，雖然并行程序仍然保持加速，但不再保持線性加速比，這與測(cè)試問(wèn)題的規(guī)模有關(guān)。

在結(jié)果準(zhǔn)確性方面，并行軟件的計(jì)算結(jié)果與串行程序計(jì)算結(jié)果完全一致。

圖6 并行程序加速比Fig. 6 Speedup of parallel POM

5 結(jié) 語(yǔ)

對(duì)POM模式串行程序進(jìn)行數(shù)據(jù)流向分析的基礎(chǔ)上，采用本地冗余計(jì)算與通信結(jié)合更新重疊區(qū)的方式，基于MPI開(kāi)發(fā)了POM并行程序。經(jīng)測(cè)試，并行軟件效率較高，達(dá)到了業(yè)務(wù)化應(yīng)用要求。目前，國(guó)家海洋信息中心已經(jīng)發(fā)布了中國(guó)海及鄰近海域的23年海洋再分析產(chǎn)品，這在中國(guó)海洋界尚數(shù)首次，其中POM海洋模式的并行軟件的研制成功為保證模式的高效運(yùn)行起了重要的作用。

表1 性能測(cè)試結(jié)果Tab. 1 Results of performance

目前POM模式串行程序，對(duì)陸地點(diǎn)和水點(diǎn)做同樣的計(jì)算，只是將陸地點(diǎn)的計(jì)算結(jié)果扣除。這其實(shí)做了許多冗余計(jì)算，使得由此改寫(xiě)的并行程序中許多進(jìn)程的計(jì)算亦有冗余，甚至有的進(jìn)程的工作完全是冗余的。我們進(jìn)一步的工作，將只考慮水點(diǎn)的計(jì)算，但要設(shè)計(jì)更好的區(qū)域劃分算法，以保證進(jìn)程的負(fù)載平衡。

參考文獻(xiàn)：

[1] Blumberg A F, Mellor G L. A description of a three-dimensional coastal ocean model [G]. Three dimensional coastal ocean models,N. S. Heaps, Editor, American Geophysical Union, Washington D C,1987: 1-16.

[2] Cousins S, Xue H. Running the Princeton Ocean Model on a Beowulf Cluster [R]. Terrain-Following Coordinates User’s Workshop, Boulder, Colorado, 2001: August 20-22.

[3] MPI: A Message-Passing Interface Standard [S]. Message Passing Interface Forum, 2003.

[4] 張林波, 遲學(xué)斌. 并行計(jì)算導(dǎo)論 [M]. 清華大學(xué)出版社, 2006:268-279.

Parallel algorithms for Princeton Ocean Model

LI Dong1,2, LIU Jing1, HAN Gui-jun2, ZHANG Xue-feng2, WANG Xi-dong2
(1. Department of Computer Science and Technology, College of Information Technical Science, Nankai University, Tianjin 300071, China;2. National Marine Data and Information Service, SOA, Tianjin 300171, China)

There is no officially published edition of parallel POM (Princeton Ocean Model) at the present time. This paper intends to investigate the key algorithms and major techniques in parallelizing the serial code of POM based on data-flow analysis. The parallel edition of POM has been developed using MPI (message passing interface) and its high performance is verified by our experiments. And now, it has been successfully applied to the reanalysis system in NMDIS (National Marine Data and Information Service)

POM; parallel algorithm; MPI; reanalysis

P717；P731.2

A

1001-6932(2010)03-0329-06

2009-11-19；

2009-12-21

國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃課題（2007CB816001）、國(guó)家自然科學(xué)基金項(xiàng)目（40776016、40906015和40906016）

李冬(1974－)，男，博士研究生，副研究員，主要從事海洋數(shù)據(jù)同化、并行計(jì)算及科學(xué)數(shù)據(jù)可視化研究。電子郵箱：lidong2003@gmail.com