劉 俊,高 陽,單桂華,遲學(xué)斌

1(中國科學(xué)院 計算機(jī)網(wǎng)絡(luò)信息中心,北京 100190)

2(中國科學(xué)院大學(xué),北京 100049)

針對CFD 流場模擬結(jié)果的可視化處理是流體動力學(xué)分析的重要一環(huán).常用的流場可視化處理方法包括流線、軌跡線、FTLE 分析等方法.而上述方法的核心算法均為粒子追蹤算法,因此,考慮到粒子追蹤算法之于流場可視化的重要性,為了提高流場可視化分析效率,領(lǐng)域研究人員針對該算法的優(yōu)化和并行化開展了一系列研究工作.現(xiàn)有的并行化算法往往考慮的是比較單一的計算架構(gòu),有的是基于共享存儲的計算環(huán)境,有的是HPC 集群計算環(huán)境,有的是結(jié)合GPU 加速器的異構(gòu)計算環(huán)境.為了兼容各類型計算架構(gòu),充分發(fā)揮異構(gòu)并行計算環(huán)境的計算資源以實現(xiàn)流線快速生成,我們基于自主研發(fā)的GPVis 平臺[1],采用了數(shù)據(jù)并行原語結(jié)合分布式消息通訊的技術(shù)架構(gòu),設(shè)計了一套流線生成系統(tǒng),通過數(shù)據(jù)并行原語實現(xiàn)設(shè)備無關(guān)的并行粒子追蹤算法,并結(jié)合數(shù)據(jù)空間域分解技術(shù)實現(xiàn)整個流場區(qū)域的流場生成任務(wù)的分而治之策略.本文基于這一系列方法在國產(chǎn)超算平臺上完成了該系統(tǒng)的實現(xiàn)和部署,并與已有的流線生成系統(tǒng)在不同參數(shù)條件下進(jìn)行了測試比較.在對測試結(jié)果進(jìn)行分析后,我們發(fā)現(xiàn)新系統(tǒng)可有效改善傳統(tǒng)流線生成方法的負(fù)載不均衡問題,以較高的并行效率實現(xiàn)分布式計算環(huán)境下的可擴(kuò)展計算.

1 研究背景

1.1 并行粒子追蹤

粒子追蹤計算是流場可視化的重要算法之一.而針對粒子追蹤算法并行化的基本方法主要可分為兩類:種子點并行法Parallelize-Over-Seeds (POS)以及數(shù)據(jù)塊并行法Parallelize-Over-Data (POD)[2].其中,POS 算法的基本思想是將初始種子點均勻分布在各處理單元上,并且讓各處理單元進(jìn)行獨立計算.而POD 算法則是通過域分解將整體數(shù)據(jù)分為眾多數(shù)據(jù)塊,并將各數(shù)據(jù)塊分配到各處理單元,然后每個處理單元負(fù)責(zé)完成進(jìn)入該數(shù)據(jù)塊的粒子的追蹤計算,直至粒子達(dá)到終止條件或進(jìn)入其它數(shù)據(jù)塊空間.

POS 方法往往需要較大的計算設(shè)備內(nèi)部存儲空間來支撐,或者需要通過頻繁的I/O 加載操作來消減一部分的內(nèi)部存儲空間需求.這是由于POS 方法并未針對數(shù)據(jù)分塊方面制定域分解策略,它要求每個節(jié)點在執(zhí)行期間可以隨機(jī)訪問全部數(shù)據(jù).這就造成各處理單元要么在啟動時讀入所有數(shù)據(jù),要么將整個數(shù)據(jù)按空間分成小塊再按需讀取.兩種方法中前者的性能容易受到設(shè)備存儲空間(系統(tǒng)內(nèi)存以及GPU 顯存)的瓶頸限制,后者則會帶來較高的I/O 相關(guān)的成本(如大量的I/O 啟動延遲所帶來的時間成本,或者為防止各計算單元的I/O 操作之間的沖突而提供的高帶寬并行I/O 能力支撐而帶來的硬件成本).

POD 方法適用于高度可擴(kuò)展的HPC 環(huán)境,可通過數(shù)據(jù)分塊將單節(jié)點的內(nèi)部存儲空間需求降低到可接受的水平,但同時由于任務(wù)的不可遷移也往往導(dǎo)致嚴(yán)重的負(fù)載不均衡問題.

鑒于兩類基本方法存在的弊端,研究人員在兩類方法基礎(chǔ)上,發(fā)展出任務(wù)竊取法[3]、任務(wù)請求法[4]等混合并行方法.混合并行方法一般會動態(tài)地且冗余地將數(shù)據(jù)塊分配給各處理單元,以實現(xiàn)更優(yōu)的負(fù)載均衡,同時也通過先驗的規(guī)劃或?qū)崟r的策略來降低冗余數(shù)據(jù)讀取所帶來的附加的I/O 操作對算法整體可擴(kuò)展性的影響.

1.2 數(shù)據(jù)并行原語

數(shù)據(jù)并行原語(Data Parallel Primitives,DPP)的設(shè)計思想最早來源于一種用于數(shù)據(jù)并行計算的掃描矢量模型[5].在數(shù)據(jù)結(jié)構(gòu)、計算幾何、圖分析以及數(shù)值分析等眾多算法領(lǐng)域,掃描矢量模型被應(yīng)用于重新設(shè)計各種算法以實現(xiàn)算法的并行化.而數(shù)據(jù)并行原語DPP 作為掃描矢量模型的升級版本,可采用包括Map、Scan、Sort和Reduce 在內(nèi)的一些基礎(chǔ)操作來實現(xiàn)各類并行算法.從目前的研究來看,除了Delaunay、FFT 等少量算法,在DPP 模型基礎(chǔ)上,大部分的可視化算法都可以有相應(yīng)的實現(xiàn)方法[6].雖然并不是所有可視化算法都可以通過DPP 來實現(xiàn),但并不妨礙DPP 模型在海量線程計算設(shè)備中的應(yīng)用.

VTK-m[7]是較典型的且功能相對全面的一套基于DPP的可視化算法庫.為了發(fā)展適用于GPU、MIC等海量線程計算環(huán)境下的可視化算法,Kenneth 等綜合借鑒了PISTON[8]、Dax[9]和EAVL[10]等關(guān)注于海量線程計算的并行可視化架構(gòu)的設(shè)計思想,并在設(shè)備及算法通用性方面進(jìn)行了優(yōu)化,于2016年提出了VTK-m并行可視化計算架構(gòu).VTK-m 采用數(shù)據(jù)并行原語DPP作為算法并行化基本單元,提供了一套與DPP 兼容的計算設(shè)備抽象模型,并針對底層硬件體系架構(gòu)進(jìn)行了兼容適配,以提高VTK-m的兼容性和可移植性.由于用戶無需考慮底層硬件實現(xiàn)細(xì)節(jié)就可以設(shè)計實現(xiàn)各類高效率的算法,因此VTK-m 有助于降低海量線程計算環(huán)境下的并行算法的實現(xiàn)難度.目前,VTK-m的計算設(shè)備抽象模型可兼容包括GPU、MIC 在內(nèi)的多種海量線程計算設(shè)備.

VTK-m 將各類通用算法以基本任務(wù)單元的形式進(jìn)行繼承封裝,而基本任務(wù)單元的運行調(diào)度是通過各類符合DPP 模型的設(shè)備適配器來驅(qū)動的(如圖1所示).這樣,基于VTK-m 設(shè)計的可視化應(yīng)用則可以通過底層設(shè)備適配器的切換而運行于各類海量線程環(huán)境.通過VTK-m 來編寫可視化算法除了可保障算法的設(shè)備通用性之外,在穩(wěn)定性和安全性的保障方面也很便利.例如,如果通過Map 映射操作來跟蹤數(shù)據(jù)的拓?fù)溥B接時,VTK-m 會在計算設(shè)備端的內(nèi)部自動完成拓?fù)浣Y(jié)構(gòu)索引的處理而避免了不安全的操作,這就使得編寫基本任務(wù)單元比直接并行操作整體數(shù)據(jù)結(jié)構(gòu)要更容易且更安全.

圖1 VTK-m 可視化應(yīng)用運行機(jī)制

1.3 基本任務(wù)單元

基本任務(wù)單元是指組成大規(guī)模數(shù)據(jù)并行處理程序的針對小規(guī)模數(shù)據(jù)的無狀態(tài)的細(xì)粒度串行程序片段[11],而采用DPP 設(shè)計并行算法的核心任務(wù)之一就是基本任務(wù)單元的定義.我們采用VTK-m 來設(shè)計海量線程加速器環(huán)境下的并行粒子追蹤算法.由于VTK-m是面向共享內(nèi)存的海量線程運行環(huán)境而設(shè)計的可視化算法庫,沒有分布式存儲及高成本通訊等問題,而且計算資源受到單機(jī)能力頂值的限制,不需要過多考慮算法的大規(guī)模擴(kuò)展問題,因此,較適合采用POS的并行策略來設(shè)計并行粒子追蹤算法.但在實際算法設(shè)計中,VTK-m與傳統(tǒng)的POS 方法也有一定的區(qū)別.比如在定義基本任務(wù)單元的層面,VTK-m 雖然采用了與傳統(tǒng)POS 方法類似的做法,將并行粒子追蹤的基本任務(wù)單元定義為單個粒子的追蹤任務(wù),但在細(xì)節(jié)上仍然與傳統(tǒng)方法有所不同.

很多并行粒子追蹤算法的基本計算單元將固定步數(shù)的單個或一組粒子追蹤作為任務(wù)單元,如文獻(xiàn)[12]中采用的辦法,這樣可以將計算粒度做到足夠小,從而有利于整體的負(fù)載均衡的實現(xiàn).但由于過小的粒度在異構(gòu)計算環(huán)境下容易產(chǎn)生額外的調(diào)度開銷,因此,VTK-m并沒有按這個方法來定義任務(wù)單元,而是采用單個粒子的整條的追蹤計算任務(wù)看作基本任務(wù)單元.這樣做既有利于降低任務(wù)調(diào)度開銷,也有利于降低任務(wù)單元之間的依賴性,從而使得采用Map 操作原語進(jìn)行算法映射更為高效.

在確定了基本任務(wù)單元后,只要將基本任務(wù)單元映射到一組執(zhí)行線程上,就可以通過海量線程計算設(shè)備的大規(guī)模并行能力來完成總體計算任務(wù).相應(yīng)的,在VTK-m的并行粒子追蹤模塊中,采用設(shè)備無關(guān)的通用方式來描述單個粒子追蹤基本操作的流程,形成Advect函子,然后通過Map 操作來實現(xiàn)以函子形式封裝的基本任務(wù)單元向底層并行計算設(shè)備的映射.Map 操作是VTK-m 中被高度優(yōu)化的DPP 之一,針對底層硬件的特點采用了定制的任務(wù)映射方法以執(zhí)行海量的基礎(chǔ)任務(wù)單元.特別是當(dāng)運行環(huán)境是GPU 這類加速器設(shè)備時,Advect 函子所涉及的輸入輸出數(shù)據(jù)都將自動傳輸?shù)皆O(shè)備內(nèi)存中,或者從設(shè)備內(nèi)存中導(dǎo)出.通過這種方式,基本任務(wù)單元的執(zhí)行位置實際上對開發(fā)人員是不可見的,這也是為什么在編寫基本任務(wù)單元時只能通過數(shù)組句柄來訪問輸入輸出數(shù)據(jù)的主要原因.

2 系統(tǒng)設(shè)計

受到單機(jī)計算存儲能力頂值的限制,大規(guī)模流場數(shù)據(jù)無法直接加載到單臺機(jī)器的主線內(nèi)存或加速器設(shè)備內(nèi)存中.因此,無法僅靠基于GPU、MIC 等海量線程計算硬件的并行粒子追蹤系統(tǒng)來滿足大規(guī)模流場數(shù)據(jù)可視化需求.為了實現(xiàn)大規(guī)模流場的并行可視化,我們需要在其基礎(chǔ)上構(gòu)建適用于異構(gòu)集群環(huán)境的并行粒子追蹤系統(tǒng).

流場可視化方法常用于處理由CFD 數(shù)值模擬產(chǎn)生的結(jié)果數(shù)據(jù)或中間數(shù)據(jù).而CFD 模擬往往會因為計算過程的特點在各類不同的計算環(huán)境下運行,因此,流場可視化軟件也需要適應(yīng)如GPU 加速器、SMP 計算機(jī)、HPC 計算機(jī)集群等各類不同的計算環(huán)境.基于VTK-m的并行粒子追蹤算法可以通過其架構(gòu)中的計算設(shè)備適配器產(chǎn)生可運行于不同的計算硬件上的可執(zhí)行程序,具有很強(qiáng)的可移植性.這是適應(yīng)于上述的科學(xué)可視化運行環(huán)境特點的設(shè)計,但美中不足的是,VTK-m 并沒有解決在大規(guī)模高性能計算機(jī)環(huán)境下開展可視化計算的問題,這是因為在類似于SMP、GPU和MIC 這種采用共享存儲結(jié)合高速總線的體系架構(gòu)中,數(shù)據(jù)訪問速度非常高,無需過多考慮由于數(shù)據(jù)遷移所產(chǎn)生的時間消耗.而HPC 環(huán)境下卻需要重點考慮如I/O、通訊等方面的負(fù)載對計算效率的影響.這也正是本文中系統(tǒng)設(shè)計方面工作所面臨的主要問題.

由于VTK-m 可以在普通計算節(jié)點以及大規(guī)模線程環(huán)境下以粒子為處理單元并行地運行流線生成算法,因此,為了實現(xiàn)異構(gòu)集群計算環(huán)境下的大規(guī)模并行粒子追蹤,我們只需要基于MPI 消息通訊接口,實現(xiàn)相關(guān)算法在HPC 條件下的動態(tài)擴(kuò)展即可.為了達(dá)到上述目標(biāo),我們設(shè)計了的適用于流場可視化的數(shù)據(jù)處理系統(tǒng)框架.如圖2所示,系統(tǒng)框架結(jié)構(gòu)包含5 大部分,分別為:函子計算層,數(shù)據(jù)分塊層,任務(wù)調(diào)度層,進(jìn)程通訊層以及渲染融合層.

圖2 并行粒子追蹤系統(tǒng)結(jié)構(gòu)圖

(1)函子計算層.函子計算層負(fù)責(zé)調(diào)度單節(jié)點中的多線程或大規(guī)模線程資源,以實現(xiàn)分配到單個節(jié)點上的可視化任務(wù)可以得到節(jié)點上足夠全面的計算能力支撐.在本系統(tǒng)中,函子計算層的主要任務(wù)是基于VTK-m的多線程環(huán)境并行處理分配到節(jié)點本地的粒子追蹤基本任務(wù)單元,因此,函子計算層主要是由VTK-m 庫以及之上的算法適配庫組成.通過繼承和調(diào)用按VTK-m設(shè)計思想封裝的基本任務(wù)單元模塊WorkletMapField,系統(tǒng)實現(xiàn)了兼容GPU、MIC 等海量線程計算設(shè)備的并行粒子追蹤模塊.

(2)進(jìn)程通訊層.進(jìn)程通訊層主要負(fù)責(zé)節(jié)點進(jìn)程間的基于MPI的消息接收和發(fā)送,以及通訊相關(guān)的一系列并行算法,主要模塊包括消息通訊模塊、數(shù)據(jù)交換模塊以及進(jìn)程協(xié)調(diào)模塊.其中消息通訊模塊是對MPI接口的封裝,可用于體量小且實時性要求較高的消息的發(fā)送和接收.數(shù)據(jù)交換模塊在MPI 協(xié)議基礎(chǔ)上將節(jié)點間交換的大塊數(shù)據(jù)封裝成MPI 數(shù)據(jù)包,并逐包跟蹤發(fā)送及接收進(jìn)度,以保證數(shù)據(jù)的順序性和完整性.進(jìn)程協(xié)調(diào)模塊包含了一系列分布式并行算法的MPI 實現(xiàn),如掃描、搜索、合并、排序等.在并行流線生成系統(tǒng)中,粒子位置信息的交換以及積分結(jié)果曲線的匯總過程中需要用到較多的通訊資源,此時需調(diào)用進(jìn)程通訊層相關(guān)模塊完成相應(yīng)任務(wù).

(3)任務(wù)分布層.任務(wù)分布層主要負(fù)責(zé)各節(jié)點及進(jìn)程中計算任務(wù)的全生命周期管理.從初期的任務(wù)創(chuàng)建,到中期的任務(wù)執(zhí)行以及任務(wù)遷移,再到最后的任務(wù)回收,任務(wù)分布層各功能模塊以任務(wù)為單元跟蹤其即時狀態(tài),并根據(jù)任務(wù)的起始條件、終止或暫停條件調(diào)整任務(wù)的執(zhí)行狀態(tài),根據(jù)任務(wù)的數(shù)據(jù)需求以及節(jié)點的負(fù)載情況,按總體任務(wù)調(diào)度策略發(fā)送或接收任務(wù),完成任務(wù)負(fù)載在節(jié)點間的流動.由于任務(wù)執(zhí)行過程采用的是基于VTK-m的多線程架構(gòu),因此任務(wù)分布層還需要在線程級實現(xiàn)基本任務(wù)單元的調(diào)度運行,從而實現(xiàn)線程級和進(jìn)程級任務(wù)的有機(jī)融合.

(4)數(shù)據(jù)分塊層.數(shù)據(jù)分塊層主要負(fù)責(zé)實現(xiàn)可視化算法并行化中的域分解策略(如循環(huán)分塊、八叉樹分塊、二分樹分塊等),并在數(shù)據(jù)分布、存儲及I/O 等方面提供可調(diào)用的接口.其中數(shù)據(jù)分布策略在并行算法中的I/O 消耗最小化、負(fù)載均衡化等方面起到非常關(guān)鍵的作用,尤其是對于粒子追蹤這類負(fù)載不可預(yù)測的情況.數(shù)據(jù)分布過程中單元數(shù)據(jù)塊大小、重疊區(qū)域的設(shè)置、數(shù)據(jù)區(qū)塊的節(jié)點分配都對可視化算法的整體運行效率以及負(fù)載均衡有很大影響.在本文中,我們針對并行粒子追蹤算法中區(qū)域劃分策略及冗余擴(kuò)充策略進(jìn)行了討論,并設(shè)計了新算法,具體方法請見第3 節(jié).

(5)渲染融合層.渲染融合層主要負(fù)責(zé)圖像相關(guān)的計算任務(wù)生成、執(zhí)行以及合并操作,主要功能包括圖像渲染、圖層融合、圖元拼合等功能.該層算法主要針對于圖像并行的并行可視化算法,而該類算法在體繪制和光線追蹤法渲染等可視化應(yīng)用中更為常見.在本文所涉及的任務(wù)范圍中,渲染融合層只負(fù)責(zé)流線在并行生成并匯總后的單機(jī)渲染.

3 算法實現(xiàn)及改進(jìn)

正如1.1 節(jié)中分析的那樣,包括POS和POD 在內(nèi)的基礎(chǔ)并行算法均存在較嚴(yán)重的弊端,而現(xiàn)有的混合并行算法有些并不適用于異構(gòu)集群環(huán)境,有些則存在各種其它問題,如基礎(chǔ)任務(wù)單元不兼容、包含整體預(yù)先分析處理步驟等.

為了避免上述問題,我們設(shè)計了一套混合并行粒子追蹤算法,在總體的并行策略上采用了POD 并行策略,而在每個節(jié)點上由于采用了多線程的共享存儲架構(gòu),因此在節(jié)點層面上為了充分利用本地的存儲計算資源采用了基于VTK-m 實現(xiàn)的POS 并行粒子追蹤策略.

以下針對本文提出的混合并行算法中的具體方法進(jìn)行說明如下.

3.1 粒子追蹤算法

設(shè)f為流場,y0為預(yù)設(shè)種子點,則由該流場生成的流線是指如下常微分方程的解:

其中,y(s)是一個關(guān)于沿流線方向的參數(shù)化步長距離s的三維位置函數(shù).求解流線的過程一般采用沿流線方向逐步積分的辦法,這個過程我們稱之為粒子追蹤過程.粒子追蹤過程其實質(zhì)是粒子在流場作用下不斷進(jìn)行積分推進(jìn)的過程,常用的粒子追蹤算法包括2 階、3 階、4 階以及4.5 階的龍格庫塔積分方法.本系統(tǒng)采用的是4 階龍格庫塔積分求解方法.

采用DPP+POD 策略的并行粒子追蹤算法讓不同的進(jìn)程分管不同的數(shù)據(jù)塊區(qū)域,其算法如算法1 所示.因為積分過程會將粒子推向一系列難以預(yù)測的空間位置而形成流線,而組成流線的這些空間位置有可能屬于分布于不同進(jìn)程的不同的數(shù)據(jù)塊區(qū)域,所以各個處理進(jìn)程在計算本地的粒子積分的同時,還需要接收可以落入本地數(shù)據(jù)區(qū)域的粒子,以及將由于積分推進(jìn)跨界進(jìn)入非本地所屬數(shù)據(jù)區(qū)域的粒子準(zhǔn)確發(fā)送到數(shù)據(jù)塊所屬進(jìn)程.由上述流程可看出,各進(jìn)程無法僅根據(jù)自身情況判斷計算過程是否可以終止,而是需要將一個全局的統(tǒng)計信息作為終止條件來進(jìn)行判斷,這就為算法的大規(guī)模并行帶來了難度,是實現(xiàn)算法高可擴(kuò)展性所必須解決的問題.另外,與傳統(tǒng)的POD 策略并行算法不同的是,由于基本任務(wù)單元的差別,積分計算前將按份抽取一定量的粒子進(jìn)行積分,每份粒子的數(shù)量要求與海量線程計算設(shè)備的線程個數(shù)盡量接近.

算法1.基于DPP+POD 策略的并行粒子追蹤算法1)加載本進(jìn)程負(fù)責(zé)的數(shù)據(jù)塊和種子點粒子,已終止粒子累積計數(shù)器置零;2)判斷本地粒子數(shù)是否為0,如為0 則到步驟9);3)從全部本地粒子中取出一份粒子;4)調(diào)用Map 操作原語在多線程設(shè)備中對取出的粒子執(zhí)行粒子追蹤,得到一部分已終止粒子和另一部分跨界粒子;5)將本輪已終止粒子數(shù)廣播到全部進(jìn)程;6)將本輪已終止粒子數(shù)累計入已終止粒子累積計數(shù)器;7)接收由其它進(jìn)程廣播發(fā)出的已終止粒子數(shù),并全部累計入已終止粒子累積計數(shù)器;8)對跨界粒子排序,按粒子當(dāng)前位置所屬進(jìn)程號將全部跨界粒子發(fā)送至所屬進(jìn)程;9)接收由其它進(jìn)程發(fā)出的跨界粒子,存入本地粒子庫;10)判斷當(dāng)前已終止粒子累積計數(shù)器值是否與粒子總數(shù)相等,如相等則退出,否則轉(zhuǎn)到步驟2).

3.2 數(shù)據(jù)分塊

實際上,由于流場本身的不確定性,粒子追蹤的負(fù)荷在空間域的分布是隨機(jī)的,因此,在沒有先驗信息的前提下,采用數(shù)據(jù)無關(guān)的分塊策略肯定是無法達(dá)到最佳的負(fù)載均衡的.這也正是循環(huán)分塊法建議在每個節(jié)點增加數(shù)據(jù)塊數(shù),以逼近更均衡的負(fù)載分布的原因.

粒子追蹤中較常用的數(shù)據(jù)分塊策略包括循環(huán)分塊法[13]、k-D 樹分塊法[14]等.循環(huán)分塊法注重于計算進(jìn)程的負(fù)載均衡,因此該算法在每個進(jìn)程分配散布于空間不同位置的數(shù)據(jù)塊,以緩解由粒子聚集所導(dǎo)致的負(fù)載失衡.為實現(xiàn)進(jìn)程中各個數(shù)據(jù)塊的空間均勻散布,循環(huán)分塊法會要求進(jìn)程個數(shù)不能被任意維度的網(wǎng)格數(shù)整除.但是,循環(huán)分塊法選擇忽略由于粒子流動本身的空間連續(xù)性而造成的算法鄰域空間強(qiáng)相關(guān)性,從而導(dǎo)致較高的節(jié)點間通訊量并極大地影響了整體并行效率.而k-D 樹分塊法想盡量避免粒子的跨節(jié)點漂移以減少節(jié)點間通訊從而提高并行效率,該方法針對粒子分布來進(jìn)行數(shù)據(jù)塊的分配,適用于面向時變流場可視化的動態(tài)分塊策略,但這種策略不適用于將粒子整條追蹤過程設(shè)置為基本任務(wù)單元的算法,因為該方法要求在追蹤過程中檢查粒子位置.另外,循環(huán)分塊法采用固定的塊大小也不利于算法的靈活性和可擴(kuò)展性.

為了改善循環(huán)分塊法的以上問題,我們采用了遞歸二分法進(jìn)行空間分塊.假設(shè)采用N個計算進(jìn)程開展并行計算,每個進(jìn)程要求分配nb個數(shù)據(jù)塊,那么,首先對整個數(shù)據(jù)空間進(jìn)行多次切分,切分原則為在3 個維度中選擇最長軸進(jìn)行切分,在多次切分后,分塊數(shù)總能超過N×nb塊,再將部分多余的分塊與相鄰分塊融合后,分塊數(shù)可達(dá)到N×nb塊,從而滿足數(shù)據(jù)分布化需求.同時,為減少維度分塊數(shù)與進(jìn)程數(shù)之間可能存在的整除關(guān)系所帶來的數(shù)據(jù)塊空間聚集的可能性,各數(shù)據(jù)塊將在隨機(jī)洗牌后再進(jìn)行循環(huán)分配.

3.3 數(shù)據(jù)冗余化

并行算法的可擴(kuò)展性與執(zhí)行過程中的通訊量以及負(fù)載均衡程度緊密相關(guān).針對并行粒子追蹤算法,數(shù)據(jù)通訊主要發(fā)生在粒子跨越數(shù)據(jù)塊邊界進(jìn)入另一個數(shù)據(jù)塊的情況,因此,當(dāng)粒子所屬數(shù)據(jù)塊足夠大時,越界情況發(fā)生的機(jī)率將會減少.另外,并行粒子追蹤過程中負(fù)載失衡的主要原因在于某些區(qū)域存在相對低速區(qū)域、渦旋、臨界點等需要長時間積分才能通過甚至無法通過的區(qū)域(又稱為陷阱區(qū)),盡管可以采用如文獻(xiàn)[13]中所述辦法,在經(jīng)過判斷后發(fā)現(xiàn)這類粒子時立刻中止該粒子積分過程,但這種辦法實際阻斷了后續(xù)結(jié)果的計算,有可能對流場分析過程產(chǎn)生干擾,而且基于VTKm的任務(wù)單元定義方法也決定了不可能在計算中途阻斷粒子積分,而這種情況也為整體系統(tǒng)的并行效率提升增加了難度.

由前面的分析可看出,不管是在降低數(shù)據(jù)通訊量或是提升負(fù)載均衡程度的解決辦法中,在無重疊的數(shù)據(jù)分塊的基礎(chǔ)上,將數(shù)據(jù)區(qū)域進(jìn)行冗余化擴(kuò)充都是有效的辦法之一,我們稱之為數(shù)據(jù)冗余化.在并行算法中,區(qū)域重疊法通過建立重疊區(qū)來擴(kuò)展數(shù)據(jù)區(qū)域是最常用的數(shù)據(jù)冗余化策略.重疊區(qū)又稱為鬼區(qū),在多類并行算法中均有采用,當(dāng)然也包括一些并行粒子追蹤算法.

3.4 進(jìn)程間通訊

并行粒子追蹤算法中主要的進(jìn)程間通訊分為兩類:一類是由于粒子越界所造成的負(fù)載流動,這類通訊是點對點的通訊,但具有隨機(jī)不可控的特性,無法預(yù)測通訊發(fā)生的位置和數(shù)量;另一類是粒子終止?fàn)顟B(tài)匯報所產(chǎn)生的進(jìn)程間通訊,這類通訊是點對全體的通訊,需要將單個進(jìn)程的計算結(jié)果匯報到所有進(jìn)程,讓每個進(jìn)程了解當(dāng)前已完成積分的粒子數(shù),在全部粒子完成積分時終止等待負(fù)載的工作循環(huán),進(jìn)入流線生成的收尾階段.

粒子在流場推動下會在不同的數(shù)據(jù)塊區(qū)域之間移動,由于不同的數(shù)據(jù)塊預(yù)先指定由不同的進(jìn)程負(fù)責(zé),因此積分計算負(fù)載將不可控地且被動地在進(jìn)程間流動,從而造成可能的負(fù)載聚集.通過前面講到的循環(huán)數(shù)據(jù)分塊方法,我們可以緩解這種負(fù)載聚集帶來的通訊量聚集.另外,由于頻繁的發(fā)送越界粒子有可能造成通訊數(shù)量過,多且接收粒子的進(jìn)程有可能在收到不足份的粒子前就開始積分計算從而造成計算時間的浪費,因此,改進(jìn)的算法會在一個時間段里累積一部分越界粒子,當(dāng)粒子數(shù)量達(dá)到足份量時才進(jìn)行發(fā)送,除非在該時間段里無法累積到足量的粒子,此時只能將全部越界粒子發(fā)出.

對于粒子終止?fàn)顟B(tài)匯報機(jī)制,原有的廣播形式的匯報算法實際上是非常昂貴的,尤其是在進(jìn)程數(shù)足夠大時,廣播通訊的規(guī)模以幾何規(guī)模增長,對整體的擴(kuò)展性能產(chǎn)生了非常大的影響.為了改善這一情況,我們設(shè)計了一套樹型匯報機(jī)制,如圖3所示.

圖3 樹型狀態(tài)匯報機(jī)制示意圖

如圖3中,所有進(jìn)程被組織成二叉樹的節(jié)點,每個進(jìn)程最多有一個上游進(jìn)程和兩個下游進(jìn)程,當(dāng)下游進(jìn)程匯報粒子終止情況下,上游進(jìn)程進(jìn)行統(tǒng)計匯總再上報上游進(jìn)程,當(dāng)最上游根進(jìn)程發(fā)現(xiàn)全部粒子積分均終止時,將向下游進(jìn)程發(fā)送退出循環(huán)的信號并退出循環(huán),其它進(jìn)程在收到退出循環(huán)信號后也進(jìn)行同樣操作.整個通訊進(jìn)程是異步且點對點的,因此對整個計算過程的影響被降到最小.經(jīng)改進(jìn)的并行算法如算法2 所示.

算法2.改進(jìn)的基于DPP+POD 策略的并行粒子追蹤算法1)加載本進(jìn)程負(fù)責(zé)的數(shù)據(jù)塊和種子點粒子,已終止粒子累積計數(shù)器置零,打開循環(huán)開關(guān);2)判斷本地粒子數(shù)是否為零,如為零則轉(zhuǎn)到步驟12);3)從全部本地粒子中取出一份粒子;4)調(diào)用Map 操作原語在多線程設(shè)備中對取出的粒子執(zhí)行粒子追蹤,得到一部分已終止粒子和另一部分跨界粒子;5)對跨界粒子排序,按粒子當(dāng)前位置所屬進(jìn)程號將全部跨界粒子發(fā)送至所屬進(jìn)程;6)接收下游進(jìn)程發(fā)來的經(jīng)匯總的已終止粒子數(shù)消息,并與本輪已終止粒子數(shù)相加得到待發(fā)送的經(jīng)匯總的已終止粒子數(shù);7)向上游進(jìn)程發(fā)送經(jīng)匯總的已終止粒子數(shù);8)將經(jīng)匯總的已終止粒子數(shù)累計入已終止粒子累積計數(shù)器;9)判斷當(dāng)前已終止粒子累積計數(shù)器值是否與粒子總數(shù)相等,如相等則轉(zhuǎn)到步驟11),否則轉(zhuǎn)到步驟12);10)接收上游發(fā)送的退出循環(huán)指令,如未收到則轉(zhuǎn)到步驟12);11)關(guān)閉循環(huán)開關(guān),并向下游進(jìn)程發(fā)送退出循環(huán)指令;12)接收由其它進(jìn)程發(fā)出的跨界粒子,存入本地粒子庫;13)判斷循環(huán)開關(guān)是否打開,如關(guān)閉則退出,否則轉(zhuǎn)到步驟2).

4 應(yīng)用測試

4.1 測試環(huán)境

我們在國產(chǎn)高性能計算集群上開展了系統(tǒng)和算法的相關(guān)測試.該集群每個節(jié)點包含32 個海光C86 7185 CPU (32 核,主頻2.0 GHz),128 GB 內(nèi)存.編譯器采用Intel編譯器,啟用OpenMP 與VTK-m 相結(jié)合實現(xiàn)多線程并行,多節(jié)點分布式并行MPI 環(huán)境采用的是OpenMPI,我們根據(jù)實驗需要使用了1 至256 個節(jié)點開展測試.

4.2 測試數(shù)據(jù)及方法

我們采用了由CFD 模擬程序?qū)嶋H產(chǎn)生的數(shù)據(jù)來開展算法測試.該測試數(shù)據(jù)是一套燃燒模擬結(jié)果數(shù)據(jù)[15,16],其網(wǎng)格分辨率為670×459×459,包含億級網(wǎng)格規(guī)模.相關(guān)的CFD 模擬實驗?zāi)M了超音速噴射的可燃?xì)怏w在空氣中混合燃燒產(chǎn)生渦流的過程,以研究燃燒反應(yīng)過程對渦結(jié)構(gòu)的影響.模擬采用的是OpenCFD-Comb 模擬程序.本次測試只采用了一個時間步,以實現(xiàn)基于可視化的穩(wěn)態(tài)流場分析.圖4顯示了該數(shù)據(jù)的流場可視化效果.

我們采用了16 至1024 個進(jìn)程開展并行粒子追蹤計算,單個粒子的追蹤步數(shù)最大值設(shè)定為10 000 步.并行化方法均采用經(jīng)過隨機(jī)排序的遞歸二分法進(jìn)行數(shù)據(jù)分塊,并采用了區(qū)域重疊最大化的數(shù)據(jù)冗余化策略,單進(jìn)程限定的內(nèi)存使用規(guī)模為64 MB.為了測試系統(tǒng)性能,我們在相同的運行環(huán)境下編譯運行了ALPINE 系統(tǒng)[17]中的粒子追蹤模塊,并基于相同的數(shù)據(jù)以及計算目標(biāo)進(jìn)行了對比測試.ALPINE 系統(tǒng)中的粒子追蹤模塊是同樣采用了VTK-m 作為核心計算庫的異構(gòu)并行可視化軟件.

圖4 測試數(shù)據(jù)可視化效果圖

4.3 性能對比分析

我們設(shè)計了一套強(qiáng)擴(kuò)展測試方案,即通過總量固定的任務(wù)量來測試并行算法在不同的節(jié)點進(jìn)程數(shù)的情況下的加速效率.對于粒子追蹤算法來說,在設(shè)置固定的粒子總數(shù)的前提下,將積分任務(wù)均勻分布到各節(jié)點或進(jìn)程上開展強(qiáng)擴(kuò)展測試,我們可以測量并觀察總的執(zhí)行時間是如何隨進(jìn)程數(shù)的變化而發(fā)生變化的.理想狀態(tài)下,加速比將線性增加.而實際上對于進(jìn)程間存在通訊的并行算法來說,在進(jìn)程數(shù)達(dá)到一定規(guī)模的情況下會出現(xiàn)加速比無法提升甚至下降的情況.

本次測試按兩種粒子規(guī)模進(jìn)行了粒子追蹤測試,其粒子規(guī)模n分別為1 兆以及2 兆.測試結(jié)果如圖5所示.其中圖5(a)為跟蹤1 兆粒子時的計算時長與核數(shù)關(guān)系圖,圖5(b)為對應(yīng)的負(fù)載均衡與核數(shù)關(guān)系圖,圖5(c)為跟蹤2 兆粒子時的計算時長與核數(shù)關(guān)系圖,圖5(d)為對應(yīng)的負(fù)載均衡與核數(shù)關(guān)系圖.從測試結(jié)果可看出,ALPINE 系統(tǒng)在核數(shù)達(dá)到512 核之后迅速表現(xiàn)出背離加速方向的發(fā)展趨勢,且在數(shù)據(jù)規(guī)模相對較大時,這種趨勢更為明顯且出現(xiàn)的時間也更早,如圖5(c)所示,針對2 兆粒子進(jìn)行粒子追蹤時,核數(shù)達(dá)到128 核之后,該系統(tǒng)就無法進(jìn)一步實現(xiàn)計算時間的縮減.相比較而言,GPVis 平臺的粒子追蹤系統(tǒng)在核數(shù)達(dá)到1024核的規(guī)模時,仍然表現(xiàn)出加速潛力,且整體加速趨勢趨近于線性加速.

4.4 負(fù)載均衡分析

圖5(b)和圖5(d)通過負(fù)載失衡度這一指標(biāo)說明了計算過程的負(fù)載均衡情況.負(fù)載失衡度是各進(jìn)程中進(jìn)程最大負(fù)載與各進(jìn)程平均負(fù)載的比值.從圖可看出,兩組測試表現(xiàn)出相似的模式,即ALPINE 在核數(shù)達(dá)到256 核之后,負(fù)載失衡的情況迅速惡化,而GPVis 在核數(shù)增加時,負(fù)載失衡度略有提升,但仍在可控范圍.從兩組測試結(jié)果中也可看出不同之處,在數(shù)據(jù)規(guī)模相對較大時,ALPINE的失衡度升高較快,而GPVis的負(fù)載失衡度受到數(shù)據(jù)規(guī)模的影響幾乎可以忽略.

圖5 系統(tǒng)強(qiáng)擴(kuò)展性能測試對比圖

在對運行過程中的具體計算日志進(jìn)行詳細(xì)記錄后,我們發(fā)現(xiàn)了兩者的主要區(qū)別.我們記錄兩組測試中512 核的運行日志,并在512 個進(jìn)程中選取了32 個進(jìn)程進(jìn)行可視化,圖6即為兩組運行測試日志的可視化結(jié)果,圖6(a)為ALPINE 運行日志的抽樣,圖6(b)為GPVis 運行日志的抽樣.圖中每個橫條代表一個進(jìn)程的工作狀態(tài)發(fā)展過程,向右方向為時間方向,兩組橫條的長度被歸一化為相同的長度,以表現(xiàn)整個計算過程.橫條中綠色代表計算時間,紅色代表通訊時間.從圖中可看出,ALPINE 系統(tǒng)由于采用了全局廣播通訊的終止粒子匯報機(jī)制,當(dāng)核數(shù)增加時,相應(yīng)的通訊時長逐漸占據(jù)主導(dǎo)而拖慢了整體性能,而GPVis 粒子追蹤系統(tǒng)采用的優(yōu)化的終止粒子匯報機(jī)制,相應(yīng)的通訊壓力則非常小,從整個計算周期來看,占據(jù)非常小的比例.這也是圖5中兩者負(fù)載失衡度變化趨勢發(fā)生區(qū)別的主要原因.

4.5 大規(guī)模測試分析

當(dāng)采用FTLE 進(jìn)行流場分析處理時,需要開展大規(guī)模的流線生成計算,計算規(guī)模為每個網(wǎng)格點設(shè)立一個被追蹤的粒子.那么對于本案例針對數(shù)據(jù)來說,粒子數(shù)取與網(wǎng)格數(shù)相當(dāng)?shù)囊?guī)模后大小約為140 兆.針對這一規(guī)模,我們開展了相應(yīng)的測試.測試結(jié)果如圖7所示,從中可看出,目標(biāo)系統(tǒng)在較大的計算負(fù)載下且當(dāng)核數(shù)達(dá)到萬核級別時依然表現(xiàn)出可接受的加速趨勢,但并行效率則在1024 核的50%之后出現(xiàn)快速下降,其中在采用2048 核、4096 核和8192 核時的并行效率分別為29%、21.3%和11.1%.

圖6 系統(tǒng)運行甘特圖對比

圖7 大規(guī)模粒子追蹤任務(wù)測試結(jié)果

4.6 測試小結(jié)

從測試結(jié)果可看出,本文提出的系統(tǒng)及算法可適用于異構(gòu)集群計算環(huán)境,且在計算速度、可擴(kuò)展性等方面表現(xiàn)優(yōu)于ALPINE 系統(tǒng)粒子追蹤模塊.經(jīng)過細(xì)致的運行日志分析,我們發(fā)現(xiàn)了文中提出進(jìn)程通訊算法優(yōu)化方法在系統(tǒng)整體的可擴(kuò)展性提升及負(fù)載均衡化方面起到了重要作用.經(jīng)過大規(guī)模測試分析,我們認(rèn)為,本系統(tǒng)可用于大規(guī)模集群環(huán)境下的并行流線生成及流場分析場景.

5 結(jié)束語

本文針對異構(gòu)集群環(huán)境下的并行流線生成系統(tǒng)提出一種改進(jìn)的基于POD+DPP的粒子追蹤算法設(shè)計實現(xiàn)方案,并針對流線生成過程中的粒子追蹤算法的并行化的實現(xiàn)細(xì)節(jié)進(jìn)行了分析和介紹.在針對流線生成系統(tǒng)的并行測試后,我們發(fā)現(xiàn)該系統(tǒng)可實現(xiàn)大規(guī)模CFD 流場的流線并行生成,達(dá)到了理想的可擴(kuò)展能力,并且在并行規(guī)模達(dá)到千核級別時仍表現(xiàn)出較高的并行效率,在核數(shù)提高至萬核級別時保持計算時長持續(xù)減少.

由于這次測試針對的數(shù)據(jù)屬于穩(wěn)態(tài)流場數(shù)據(jù),因此并未涉及非穩(wěn)態(tài)流場的分析問題.穩(wěn)態(tài)流場與非穩(wěn)態(tài)流場的流場或軌跡線生成方法有很大的差別,前者更注重計算的分配,后者更注重數(shù)據(jù)的動態(tài)加載過程,后續(xù)我們將對本系統(tǒng)進(jìn)行進(jìn)一步升級,并應(yīng)用于時變流場的可視化處理.另外,從性能分析結(jié)果中可看出,本文提出的優(yōu)化方法雖然較傳統(tǒng)方法提高了計算時間占用率,但在進(jìn)程數(shù)較多時并行效率仍然處于較低的水平,存在進(jìn)一步優(yōu)化的空間.