面向天河2A系統(tǒng)的基于蒙特卡羅方法的粒子輸運異構(gòu)協(xié)同計算 *

李 彪，劉 杰，2

(1.國防科技大學并行與分布處理國家重點實驗室，湖南 長沙 410073,2.復雜系統(tǒng)軟件工程湖南省重點實驗室，湖南 長沙 410073)

1 引言

1.1 粒子輸運簡介

輸運理論是研究微觀粒子在介質(zhì)中遷移統(tǒng)計規(guī)律的數(shù)學理論。這里“微觀粒子”指的是中子、光子、分子、電子和離子等。輸運理論研究大量粒子在空間或者某種介質(zhì)中運行時，由于各粒子位置、動量和其他特征量的變化而引起的各種有關(guān)物理量隨時空變化所表現(xiàn)出來的非平衡統(tǒng)計運動規(guī)律[1]。

19世紀中后葉的分子運動論導致粒子輸運理論發(fā)展的開端。1872 年，Boltzmann 推導出了Maxwell分子速度分布函數(shù)隨時空變化的非線性積分-微分方程，該方程也常被稱為粒子輸運方程或者Boltzman方程[1]。粒子輸運方程描述的分布函數(shù)具有7個獨立變量，包括3個空間變量、2個方向角度變量、1個能量變量和1個時間變量。由于實際情況的復雜性，通常不能采用解析方式進行直接求解，只能通過數(shù)值模擬求解。

目前常用的求解粒子輸運問題的數(shù)值方法分成2類[2]：一類是“確定性方法”，這類方法通過對變量進行離散，將粒子輸運方程轉(zhuǎn)化為一組線性代數(shù)方程，然后通過求解此代數(shù)方程組來獲得精確解或近似解；另一類方法是“蒙特卡羅MC(Monte Carlo)方法”，也稱非確定性方法或概率論方法。它是基于統(tǒng)計理論的數(shù)值方法，對所要研究的問題構(gòu)造一個隨機概率模型來模擬大量粒子的運動，以獲得感興趣的物理量分布。

本文研究的是關(guān)于粒子輸運中的非確定性方法。

1.2 蒙特卡羅方法在粒子輸運上的應(yīng)用

蒙特卡羅MC方法亦稱隨機模擬法，是20世紀40年代中期電子計算機誕生后發(fā)展起來的一門計算科學[3]。它是通過計算機來對中子行為進行隨機模擬的數(shù)值方法，是一種以概率統(tǒng)計理論為指導的非常重要的數(shù)值計算方法。

通常蒙特卡羅方法可以粗略地分成2類[3]：一類是所求解的問題本身具有內(nèi)在的隨機性，利用計算機直接模擬這種隨機的過程。例如，在中子輸運理論中，要求建立單個中子在給定幾何系統(tǒng)中的真實運動歷史，通過對大量中子歷史的跟蹤，得到充分的隨機試驗值(或稱抽樣值)，然后用統(tǒng)計的方法得出隨機變量某個數(shù)值特征的估計量，用此估計量作為問題的解。另一類是所求解問題可以轉(zhuǎn)化為某種隨機分布的特征數(shù)，比如隨機事件出現(xiàn)的概率，或者隨機變量的期望值。通過隨機抽樣的方法，以隨機事件出現(xiàn)的頻率估計其概率，或者以抽樣的數(shù)字特征估算隨機變量的數(shù)字特征，并將其作為問題的解。這種方法多用于求解復雜的多維積分問題。

1.3 基于MC方法的異構(gòu)協(xié)同計算

由于功耗和散熱問題日益突出，單塊芯片的性能已經(jīng)不能單純地從提高頻率來獲取提升，摩爾定律將會消亡[4]，于是工業(yè)界逐漸轉(zhuǎn)向多核處理器。但是，隨著技術(shù)的發(fā)展，多核處理器中每個內(nèi)核的設(shè)計越來越復雜，并且頻率較高，功耗較大，嚴重阻礙了多核處理器中核數(shù)的擴展。為了進一步提升性能，同時滿足功耗和散熱的要求，出現(xiàn)了由眾多相對簡單的內(nèi)核構(gòu)建的眾核協(xié)處理器，如 NVIDIA公司推出的 GPU[5]、Intel公司推出的MIC協(xié)處理器[6]和本文所探討的國產(chǎn)加速器Matrix2000等。與多核處理器相比，這些眾核加速器的內(nèi)核相對簡單，且頻率較低，但數(shù)量眾多，能夠以合理的功耗實現(xiàn)更高的性能。

近十年來，異構(gòu)加速器快速發(fā)展，TOP500排行中排名靠前的超算系統(tǒng)的體系結(jié)構(gòu)均采用了異構(gòu)協(xié)同模式，許多研究人員對已有的求解粒子輸運問題的方法在各種加速平臺上做了許多相關(guān)的研究工作。

楊博[7]提出了一種針對CPU/GPU混合異構(gòu)系統(tǒng)的深穿透粒子輸運MC模擬異構(gòu)協(xié)同算法，針對GPU的計算和訪存特點，給出了一種基于粒子數(shù)的任務(wù)劃分方法和高效并行數(shù)據(jù)結(jié)構(gòu)以及線程間的并行規(guī)約方法。相比運行在Intel Xeon X5670上的MC程序，該算法獲得了3.4的加速比，并在TianHe-1A的64個結(jié)點進行了測試，結(jié)果表明該算法具有良好的性能和可擴展性。

崔顯濤[8]基于MC方法，提出了一種面向CPU/MIC混合異構(gòu)系統(tǒng)的粒子輸運并行算法，針對MIC訪存特點提出了適于程序并行的高速數(shù)據(jù)結(jié)構(gòu)和基于靜態(tài)分配的任務(wù)劃分方式。相比運行在CPU端的程序串行程序，改進的MC程序獲得了8.6倍的加速比。

在多能群MC粒子輸運領(lǐng)域，近些年來有許多研究和探索,Bergmann 等人[9]基于GPU開發(fā)了一個MC粒子輸運框架WARP(Weaving All the Random Particles)。WARP的目的是使以前的基于事件的輸運算法適應(yīng)新的基于GPU硬件而開發(fā)的連續(xù)能量Monte Carlo中子輸運代碼。Hamilton等人[10]基于Shift程序包開發(fā)了一個在GPU上運行的連續(xù)能量Monte Carlo中子輸運求解器，并把程序移植到了超算系統(tǒng)Summit上，在耗盡燃料的SMR(Small Modular Reactor)核心配置上，進行了1 024個結(jié)點的弱擴展測試，結(jié)果顯示其并行效率接近93%，這說明GPU對性能提升顯著。

本文的工作與文獻[8,9]的接近，都是只考慮單群的MC程序，不過本文的異構(gòu)算法設(shè)計是基于不同于GPU和MIC體系結(jié)構(gòu)的國產(chǎn)加速器Matrix2000，且研制的程序可以移植到大規(guī)模結(jié)點上，另外，還針對MC程序的數(shù)據(jù)通信進行了優(yōu)化，取得了顯著的性能提升。

2 粒子輸運MC模擬

MC模擬是通過對大量中子歷史的跟蹤來獲取大量的隨機試驗結(jié)果，并對這些結(jié)果進行統(tǒng)計方法處理，最終以得到的統(tǒng)計量作為問題的解。 所謂一個中子的歷史，是指該中子從源出發(fā)，在介質(zhì)中隨機游走，經(jīng)過各種核反應(yīng)作用，直到中子歷史結(jié)束或稱中子“死亡”。所謂“死亡”是指中子被吸收、穿出系統(tǒng)、被熱化，或達到能量權(quán)下限或時間上限。其中時間、能量的載斷是無條件的，而權(quán)截斷是有條件的，由俄羅斯輪盤賭決定[11]。

圖1是一個中子歷史的循環(huán)過程。初始時，從粒子源出發(fā)的一個粒子的位置、方向、能量等參數(shù)是確定的，隨后粒子的運動方向和與原子核的碰撞類型服從特定的概率分布，在經(jīng)歷了一系列碰撞過程之后，粒子歷史趨向于結(jié)束。有幾種不同類型的結(jié)束，比如達到時間與能量界限、逃脫出系統(tǒng)或達到權(quán)下限。在這一個過程中，下一次碰撞僅與當前碰撞后粒子的位置、能量以及方向相關(guān)，完全獨立于先前的碰撞。 因此，這一過程是一個典型的馬爾可夫過程。所以，只要知道粒子與原子核碰撞的規(guī)律，那么粒子的軌跡就可以用蒙特卡羅方法正確地模擬出來，從而得到所關(guān)心的解[12]。

Figure 1 Simulation of particle transport process圖1 模擬粒子輸運過程

3 通信模式優(yōu)化

3.1 串行阻塞通信模式

通過測試MC程序的強擴展性發(fā)現(xiàn)，當進程數(shù)大于256時，效率急速下降，達到1 024進程時基本上已經(jīng)沒有加速效果了，這對移植MC程序到大規(guī)模系統(tǒng)上產(chǎn)生了巨大障礙。經(jīng)過仔細分析程序的通信模式發(fā)現(xiàn)，由于計算結(jié)果數(shù)據(jù)收集的通信模式是主從模式，導致等待時間過長。

首先0號進程初始化任務(wù)，接著把數(shù)據(jù)發(fā)送給其他所有從進程，從進程接收到數(shù)據(jù)后開始計算，并在從進程的計算結(jié)束后設(shè)置同步柵欄，直到所有從進程的計算都結(jié)束后，0號進程才開始收集從進程發(fā)送過來的結(jié)果數(shù)據(jù)。數(shù)據(jù)的收集過程中所有從進程都向主進程阻塞發(fā)送數(shù)據(jù)，主進程進行阻塞接收，先到的信息先被處理。對于主進程來說這是一個串行處理過程，此通信方式會造成通信瓶頸，核數(shù)越多，對程序性能的影響越大。圖2展示了P個進程的通信過程。

Figure 2 Serial data collection mode圖2 串行數(shù)據(jù)收集模式

3.2 優(yōu)化后的二叉樹通信模式

考慮到MC程序0號進程對收集到的數(shù)據(jù)只進行簡單的統(tǒng)計處理，所以可以通過逐層進行兩兩進程通信，實現(xiàn)局部統(tǒng)計處理，最終把統(tǒng)計結(jié)果匯總到0號進程。本文采用二叉樹通信模式來實現(xiàn)這種通信過程，假設(shè)總進程數(shù)為2K+res，其中,K,res都是大于或等于0的整數(shù)，且滿足條件0≤res<2K，為了滿足二叉樹的通信模式，排在前面的2*res個進程先就近奇偶進程號兩兩通信，并把通信過程中已發(fā)送過信息的進程剔除，不參與下次通信，從而下一階段通信的進程數(shù)變?yōu)?K。接著對2K個進程號重排，前2*res個進程號除以2，其余進程依次減去res，得到新的2K個進程的標識號，下一階段就是進行二叉樹通信，算法1描述了具體的通信過程。為了方便描述，圖3展示了7個進程的通信過程，括號里的數(shù)字表示重排后的進程標識號。

Figure 3 Binary tree data collection mode圖3 二叉樹數(shù)據(jù)收集模式

改進后的數(shù)據(jù)收集模式，通信復雜度由2K+res減少為log (2K+res)，極大地減少了通信時間，也避免了大規(guī)模通信時從進程同時向主進程發(fā)送數(shù)據(jù)導致的程序阻塞。

4 基于CPU/Matrix2000的異構(gòu)協(xié)同算法

4.1 Matrix2000加速器

天河2A系統(tǒng)中，每個結(jié)點由2顆Intel Xeon 微處理器和2顆Matrix2000加速器組成，如圖4所示。每個Intel Xeon微處理器包含12核，工作頻率為2.2 GHz，采用英特爾Ivy Bridge微架構(gòu)，峰值性能0.211 2 TFLPOS。每個Matrix2000加速器包含128核，由4個超結(jié)點組成，每個超結(jié)點包含32個計算核，其中超結(jié)點支持64核超線程技術(shù)，峰值性能2.457 6 TFLOPS@1.2 GHz，有8個DDR4內(nèi)存通道，支持×16 PCIE 3.0 EP工作模式。

算法1二叉樹通信偽代碼

//假設(shè)總進程數(shù)為2K+res，0≤res<2K

1mask=1;

2if(myid< 2*res)then/*前2*res個進程先通信，使下一階段通信進程個數(shù)為2K*/

3if(myid%2!=0)then//奇數(shù)號進程

4dest=myid-1;

5new_id=-1/*進程號置為-1，不參與下次通信*/

6myid進程SENDmsgtodest;

7else//偶數(shù)號進程

8src=myid+1;

9new_id=myid/2;//重排后新的進程號

10myid進程RECVmsgfromsrc

11endif

12else

13new_id=myid-res;/*其余的進程重排后的進程號*/

14endif

15if(new_id!=-1)then/*重排后的2K個進程進行通信*/

16while(mask<2K)do//二叉樹通信模式

17if(mask&new_id)then/*獲取重排號之后的偶數(shù)進程號*/

18newsrc=new_id|mask;

19if(newsrc

20src=newsrc*2

21else

22src=newsrc+res

23endif

24 RECVmsgfromsrc

25else

26newdest=new_id& ～mask

27if(newdest

28dest=newdest* 2

29else

30dest=newdest+res

31endif

32 SENDmsgtodest

33 exit//發(fā)送過消息后，此進程就退出通信

34endif

35mask= 2*mask

36endwhile

37endif

Figure 4 Structure of Tianhe-2A system node圖4 天河2A系統(tǒng)結(jié)點結(jié)構(gòu)

4.2 異構(gòu)編程模型

天河2A系統(tǒng)支持3種異構(gòu)通信模式OpenMP4.5、ACL和BCL，其中BCL是一種簡單高效的對稱傳輸接口，CPU和協(xié)處理器之間數(shù)據(jù)通信利用PCIE總線進行傳輸，底層通過SCIF來實現(xiàn)。BCL相比OpenMP4.5更底層，程序移植也更復雜，但是傳輸速率更快，移植后的程序靈活性更好?；贐CL接口的異構(gòu)程序需要編譯2套程序，CPU和加速器分別編譯一套程序，2套程序分別在CPU端和加速器端同時運行，其中運行時加速器端的程序需要通過ACL傳送到加速器端并啟動。

圖5給出了基于CPU/Matrix2000的MC程序的異構(gòu)模式流程。首先，CPU端0號進程啟動MPI進行進程初始化，0號進程負責讀取文件數(shù)據(jù)，并把總計算任務(wù)按照粒子數(shù)均等的方式分配給從進程對應(yīng)的子任務(wù)，再將子任務(wù)的初始化數(shù)據(jù)傳輸給對應(yīng)的從進程。每個從進程控制一個Matrix2000超結(jié)點，通過ACL把加速器端的程序傳輸?shù)組atrix2000端并啟動程序，利用BCL使CPU端與 Matrix2000超結(jié)點建立連接，CPU端把初始化數(shù)據(jù)傳輸?shù)組atrix2000端，加速器端完成CPU端數(shù)據(jù)接收后就利用多線程技術(shù)并行跟蹤每個粒子歷史，直到所有粒子歷史全部完成計算，加速器統(tǒng)計結(jié)果數(shù)據(jù)并把數(shù)據(jù)通過BCL傳輸給對應(yīng)的CPU進程，再由MPI實現(xiàn)進程間的傳輸。整個異構(gòu)控制的過程中，CPU主要負責傳輸數(shù)據(jù)，Matrix2000主要負責計算。

Figure 5 Flowchart of heterogeneous logic 圖5 異構(gòu)模式邏輯流程圖

4.3 OpenMP線程級并行

Matrix2000中一個超結(jié)點包含32核，支持64核超線程技術(shù)，為了充分發(fā)揮Matrix2000的性能，通過OpenMP指令實現(xiàn)線程級細粒度并行，算法2描述了具體算法偽代碼。加速器端接收CPU端發(fā)送過來的數(shù)據(jù)，接著做初始的工作，針對粒子輸運模塊進行多線程并行。每個線程跟蹤一個粒子，完成當前粒子跟蹤后，進入time_to_stop判斷是否所有粒子數(shù)都已完成，如果未完成，當前線程臨界更新粒子數(shù)變量nps++，更新后的第nps個粒子分配給當前線程進行跟蹤模擬。這種調(diào)度方式類似迭代塊大小等于1的dynamic調(diào)度。

算法2Matrix2000加速器上的OpenMP線程級并行

1 假設(shè)線程數(shù)為T，加速器分配到M個粒子，粒子編號區(qū)間為(nps,nps+M-1);

2 #pragma omp parallel;/*開啟多線程，每個線程跟蹤一個粒子，完成后繼續(xù)跟蹤下一個粒子*/

3forall threadsdo

4 各個線程初始化自己的私有變量;

5omp_set_lock(&lock);/*設(shè)置鎖，當前只有獲得鎖的線程才能更新共享數(shù)據(jù)*/

6 更新共享變量;

7while(!time_to_stop())do/*判斷當前是否完成M個粒子數(shù)追蹤*/

8nps=nps+1;/*未完成，更新nps，追蹤第nps個粒子*/

9 更新共享變量;

10omp_unset_lock(&lock);/*釋放鎖，使得多線程可以調(diào)用hstory并行追蹤粒子*/

11 callhstory();/*線程各自調(diào)用hstory()函數(shù)，互不干涉，完成輸運過程*/

12omp_set_lock(&lock);/*設(shè)置鎖，進入time_to_stop()中訪問共享數(shù)據(jù)*/

13endwhile

14 更新共享變量;/*完成M個粒子數(shù)追蹤，進行數(shù)據(jù)更新*/

15omp_unset_lock(&lock);

16endfor

相比Intel MIC協(xié)處理器Offload模式多線程技術(shù)以及OpenMP4.5提供的異構(gòu)模式，本文的異構(gòu)模式中的CPU與Matrix2000加速器并不是主從關(guān)系，而是對等的關(guān)系。Matrix2000加速器端啟動的是單獨的一套程序，所有的數(shù)據(jù)傳輸和初始化工作在啟動多線程之前已經(jīng)完成，所以Matrix2000上的OpenMP不需要通過編譯指導語句來傳輸數(shù)據(jù)和設(shè)置各種變量屬性，使得程序結(jié)構(gòu)設(shè)計更簡單、更靈活。

5 實驗及結(jié)果分析

5.1 測試環(huán)境及參數(shù)

測試平臺為天河2A系統(tǒng)，由于ACL和BCL指令只提供C/C++的接口，本文需要對MC程序插入C語言的通信控制接口，CPU端的程序采用Intel的編譯器和基于高速網(wǎng)的mpich3.2進行編譯；加速器端采用自定義交叉編譯器進行編譯，支持OpenMP指令。

單結(jié)點硬件配置為2顆Intel E5-2692 v2 @2.20 GHz(12 核/顆)+8個超結(jié)點Matrix2000，每個超結(jié)點通過超線程技術(shù)實現(xiàn)45核多線程并行。 MC程序的輸入文件參數(shù)除問題規(guī)模隨具體實驗設(shè)置外，均保持默認值。

5.2 通信優(yōu)化測試

考慮到數(shù)據(jù)的通信僅在CPU端進程之間進行，而并不涉及到加速器，所以本節(jié)僅在CPU上進行對比測試，表1給出了通信優(yōu)化前后的測試結(jié)果。

Table 1 Test results before and after communication optimization

分析圖6可以發(fā)現(xiàn),原始程序測試結(jié)果中，512核是拐點，小于512核時程序還有加速效果，大于512核時測試時間反而急速上升，已經(jīng)沒有加速效果，表明此時計算能力不是瓶頸，數(shù)據(jù)通信占主要耗時。相反，優(yōu)化了通信模式的程序測試時，運行時間隨進程數(shù)遞增逐漸下降，沒有出現(xiàn)進程數(shù)越多測試時間反而上升的現(xiàn)象，當大于2 000核時耗時趨于穩(wěn)定，這是由于粒子數(shù)計算規(guī)模受限，增加進程數(shù)的優(yōu)勢已經(jīng)不明顯導致的?？梢栽O(shè)想如果計算規(guī)模增加進程數(shù)越多，優(yōu)化通信后的程序會比原始程序效果更好。

Figure 6 Comparison of two communication modes圖6 通信模式對比圖

5.3 異構(gòu)算法測試

本節(jié)對基于CPU/Matrix2000的異構(gòu)MC程序進行弱擴展測試，每個結(jié)點啟動10個CPU進程，其中8個進程分別和結(jié)點內(nèi)的8個Matrix2000超結(jié)點進行異構(gòu)計算，剩余2個進程協(xié)同計算，每個超結(jié)點使用48核多線程，表2給出了弱擴展測試結(jié)果，圖7給出了相應(yīng)測試結(jié)果的柱狀圖。

Table 2 Weakly scalable test

Figure 7 Bar and curve chart of weakly scalabletest result on Tianhe-2A圖7 天河2A弱擴展測試結(jié)果柱狀圖

表2說明弱擴展到45萬核時，相對5萬核并行效率保持在 22.54%。 分析圖7可以得到，雖然研制的程序可以運行到45萬核，但是可以看到，隨著核數(shù)增大并行效率的下降趨勢很明顯。結(jié)合程序與加速器體系結(jié)構(gòu)分析，導致并行效率下降迅速的原因主要有2點：首先，MC程序的粒子輸運部分的計算不是數(shù)據(jù)密集型，在中子的歷史過程中充滿了大量粒子條件狀態(tài)的判斷，程序的邏輯控制語句居多，這與Matrix2000加速器適合程序邏輯控制簡單、數(shù)據(jù)計算密集的條件相沖突；其次，在程序的OpenMP多線程部分，雖然開啟了多線程技術(shù)，但是粒子歷史的模擬追蹤過程中各個線程都涉及更新共享數(shù)據(jù)，需要設(shè)置臨界區(qū)，在臨界區(qū)內(nèi)同時只能有一個線程更新共享數(shù)據(jù)，其他線程只能等待，這會影響多線程的并行效率。

6 結(jié)束語

本文在現(xiàn)有粒子輸運蒙特卡羅模擬程序的基礎(chǔ)上，提出了一種面向CPU/Matrix2000異構(gòu)系統(tǒng)的粒子輸運異構(gòu)協(xié)同算法，基于天河2A系統(tǒng)的異構(gòu)通信模式BCL和ACL，針對CPU和加速器各自設(shè)計了一套不同的代碼，在CPU端通過ACL實現(xiàn)加速器端代碼的啟動，利用BCL進行CPU與Matrix2000之間的通信，進而提出了一種CPU與加速器Matrix2000之間的簡單高效的對稱通信模式。優(yōu)化了原MC程序的串行數(shù)據(jù)收集通信模式，提出了新的二叉樹通信模式，極大地減少了通信時間，加速比可達17.7。通過優(yōu)化通信模式，以及基于MPI-ACL-OpenMP編程框架，本文實現(xiàn)的基于CPU-Matrix2000異構(gòu)協(xié)同計算的并行程序，可以弱擴展到45萬核，相對5萬核并行效率保持在22.54%。作為未來工作，將探索如何繼續(xù)增強該異構(gòu)程序的擴展性，以更好地發(fā)揮出HPC平臺的算力優(yōu)勢。