李 彪，劉 杰，2

(1.國防科技大學(xué)并行與分布處理國家重點(diǎn)實(shí)驗(yàn)室，湖南 長沙 410073,2.復(fù)雜系統(tǒng)軟件工程湖南省重點(diǎn)實(shí)驗(yàn)室，湖南 長沙 410073)

1 引言

1.1 粒子輸運(yùn)簡介

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 粒子輸運(yùn)MC模擬

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

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

Figure 1 Simulation of particle transport process圖1 模擬粒子輸運(yùn)過程

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

3.1 串行阻塞通信模式

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

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

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

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

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

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

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

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

4.1 Matrix2000加速器

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

算法1二叉樹通信偽代碼

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

1mask=1;

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

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

4dest=myid-1;

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

6myid進(jìn)程SENDmsgtodest;

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

8src=myid+1;

9new_id=myid/2;//重排后新的進(jìn)程號

10myid進(jìn)程RECVmsgfromsrc

11endif

12else

13new_id=myid-res;/*其余的進(jìn)程重排后的進(jìn)程號*/

14endif

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

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

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

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ā)送過消息后，此進(jìn)程就退出通信

34endif

35mask= 2*mask

36endwhile

37endif

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

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

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

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

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

4.3 OpenMP線程級并行

Matrix2000中一個超結(jié)點(diǎn)包含32核，支持64核超線程技術(shù)，為了充分發(fā)揮Matrix2000的性能，通過OpenMP指令實(shí)現(xiàn)線程級細(xì)粒度并行，算法2描述了具體算法偽代碼。加速器端接收CPU端發(fā)送過來的數(shù)據(jù)，接著做初始的工作，針對粒子輸運(yùn)模塊進(jìn)行多線程并行。每個線程跟蹤一個粒子，完成當(dāng)前粒子跟蹤后，進(jìn)入time_to_stop判斷是否所有粒子數(shù)都已完成，如果未完成，當(dāng)前線程臨界更新粒子數(shù)變量nps++，更新后的第nps個粒子分配給當(dāng)前線程進(jìn)行跟蹤模擬。這種調(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è)置鎖，當(dāng)前只有獲得鎖的線程才能更新共享數(shù)據(jù)*/

6 更新共享變量;

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

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

9 更新共享變量;

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

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

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

13endwhile

14 更新共享變量;/*完成M個粒子數(shù)追蹤，進(jìn)行數(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加速器端啟動的是單獨(dú)的一套程序，所有的數(shù)據(jù)傳輸和初始化工作在啟動多線程之前已經(jīng)完成，所以Matrix2000上的OpenMP不需要通過編譯指導(dǎo)語句來傳輸數(shù)據(jù)和設(shè)置各種變量屬性，使得程序結(jié)構(gòu)設(shè)計更簡單、更靈活。

5 實(shí)驗(yàn)及結(jié)果分析

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

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

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

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

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

Table 1 Test results before and after communication optimization

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

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

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

本節(jié)對基于CPU/Matrix2000的異構(gòu)MC程序進(jìn)行弱擴(kuò)展測試，每個結(jié)點(diǎn)啟動10個CPU進(jìn)程，其中8個進(jìn)程分別和結(jié)點(diǎn)內(nèi)的8個Matrix2000超結(jié)點(diǎn)進(jìn)行異構(gòu)計算，剩余2個進(jìn)程協(xié)同計算，每個超結(jié)點(diǎn)使用48核多線程，表2給出了弱擴(kuò)展測試結(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弱擴(kuò)展測試結(jié)果柱狀圖

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

6 結(jié)束語

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