適用于異構(gòu)集群的混合并行流線生成系統(tǒng)①

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

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

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

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

1 研究背景

1.1 并行粒子追蹤

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

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

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

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

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

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

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

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

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

1.3 基本任務(wù)單元

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 粒子追蹤算法

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

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

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

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

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

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

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

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

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

并行算法的可擴(kuò)展性與執(zhí)行過程中的通訊量以及負(fù)載均衡程度緊密相關(guān).針對(duì)并行粒子追蹤算法,數(shù)據(jù)通訊主要發(fā)生在粒子跨越數(shù)據(jù)塊邊界進(jìn)入另一個(gè)數(shù)據(jù)塊的情況,因此,當(dāng)粒子所屬數(shù)據(jù)塊足夠大時(shí),越界情況發(fā)生的機(jī)率將會(huì)減少.另外,并行粒子追蹤過程中負(fù)載失衡的主要原因在于某些區(qū)域存在相對(duì)低速區(qū)域、渦旋、臨界點(diǎn)等需要長時(shí)間積分才能通過甚至無法通過的區(qū)域(又稱為陷阱區(qū)),盡管可以采用如文獻(xiàn)[13]中所述辦法,在經(jīng)過判斷后發(fā)現(xiàn)這類粒子時(shí)立刻中止該粒子積分過程,但這種辦法實(shí)際阻斷了后續(xù)結(jié)果的計(jì)算,有可能對(duì)流場分析過程產(chǎn)生干擾,而且基于VTKm的任務(wù)單元定義方法也決定了不可能在計(jì)算中途阻斷粒子積分,而這種情況也為整體系統(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ù)載流動(dòng),這類通訊是點(diǎn)對(duì)點(diǎn)的通訊,但具有隨機(jī)不可控的特性,無法預(yù)測通訊發(fā)生的位置和數(shù)量;另一類是粒子終止?fàn)顟B(tài)匯報(bào)所產(chǎn)生的進(jìn)程間通訊,這類通訊是點(diǎn)對(duì)全體的通訊,需要將單個(gè)進(jìn)程的計(jì)算結(jié)果匯報(bào)到所有進(jìn)程,讓每個(gè)進(jìn)程了解當(dāng)前已完成積分的粒子數(shù),在全部粒子完成積分時(shí)終止等待負(fù)載的工作循環(huán),進(jìn)入流線生成的收尾階段.

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

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

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

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

算法2.改進(jìn)的基于DPP+POD 策略的并行粒子追蹤算法1)加載本進(jìn)程負(fù)責(zé)的數(shù)據(jù)塊和種子點(diǎn)粒子,已終止粒子累積計(jì)數(shù)器置零,打開循環(huán)開關(guān);2)判斷本地粒子數(shù)是否為零,如為零則轉(zhuǎn)到步驟12);3)從全部本地粒子中取出一份粒子;4)調(diào)用Map 操作原語在多線程設(shè)備中對(duì)取出的粒子執(zhí)行粒子追蹤,得到一部分已終止粒子和另一部分跨界粒子;5)對(duì)跨界粒子排序,按粒子當(dāng)前位置所屬進(jìn)程號(hào)將全部跨界粒子發(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ù)累計(jì)入已終止粒子累積計(jì)數(shù)器;9)判斷當(dāng)前已終止粒子累積計(jì)數(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)高性能計(jì)算集群上開展了系統(tǒng)和算法的相關(guān)測試.該集群每個(gè)節(jié)點(diǎn)包含32 個(gè)海光C86 7185 CPU (32 核,主頻2.0 GHz),128 GB 內(nèi)存.編譯器采用Intel編譯器,啟用OpenMP 與VTK-m 相結(jié)合實(shí)現(xiàn)多線程并行,多節(jié)點(diǎn)分布式并行MPI 環(huán)境采用的是OpenMPI,我們根據(jù)實(shí)驗(yàn)需要使用了1 至256 個(gè)節(jié)點(diǎn)開展測試.

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

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

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

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

4.3 性能對(duì)比分析

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

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

4.4 負(fù)載均衡分析

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

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

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

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

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

圖6 系統(tǒng)運(yùn)行甘特圖對(duì)比

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

4.6 測試小結(jié)

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

5 結(jié)束語

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

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