李冬，劉璟，韓桂軍，張學峰，王喜冬

（1. 南開大學信息技術科學學院，計算機科學與技術系，天津，300071；2. 國家海洋局國家海洋信息中心，天津，300171）

POM海洋模式的并行算法

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

（1. 南開大學信息技術科學學院，計算機科學與技術系，天津，300071；2. 國家海洋局國家海洋信息中心，天津，300171）

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

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

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

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

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

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

文章第二部分對POM海洋模式進行簡單介紹；第三部分詳細討論POM模式的并行算法、并行程序框架及通信軟件模塊；第四部分介紹并行程序的測試結果；第五部分給出結論。

2 POM海洋模式簡介

POM模式的主要特點為：

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

(2) 垂向采用σ坐標。

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

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

(5) 采用自由表面及時間分裂格式；外模式(正壓模)為二維，采用較短的時間步長（由CFL條件及外重力波波速決定）；內模式(斜壓模)為三維，采用較長的時間步長（由CFL條件及內波波速決定）。

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

3 POM模式并行算法及軟件設計

3.1 算法及程序框架

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

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

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

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

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

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

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

我們在對POM模式串行程序進行數(shù)據(jù)流向分析的基礎上，采用本地冗余計算與通信結合更新重疊區(qū)的方式，基于MPI開發(fā)了POM并行程序，基本流程見圖2。

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

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

3.2 通信模塊的設計和開發(fā)

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

這些基礎通信模塊包括：區(qū)域劃分及進程的虛擬拓撲結構模塊、數(shù)據(jù)重疊區(qū)更新模塊、派生數(shù)據(jù)類型模塊、數(shù)據(jù)合并及分解模塊、I/O接口模塊等。

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

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

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

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

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

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

圖3 進程虛擬拓撲結構示意圖Fig. 3 Virtual topologies of processes

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

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

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

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

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

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

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

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

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

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

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

(3) 點在整個計算區(qū)域中的全局坐標轉換為其在相應進程中的局地坐標，并返回此進程號及其所管轄區(qū)域的起始點坐標。

(4) 點在進程中的局地位置坐標轉換為其在整個計算區(qū)域中的全局坐標。

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

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

4 測試結果

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

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

在結果準確性方面，并行軟件的計算結果與串行程序計算結果完全一致。

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

5 結 語

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

表1 性能測試結果Tab. 1 Results of performance

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

參考文獻：

[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] 張林波, 遲學斌. 并行計算導論 [M]. 清華大學出版社, 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

國家重點基礎研究發(fā)展計劃課題（2007CB816001）、國家自然科學基金項目（40776016、40906015和40906016）

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