宋佩濤，張志儉，梁 亮，張 乾，趙 強(qiáng)

(哈爾濱工程大學(xué) 核安全與仿真技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，黑龍江 哈爾濱 150001)

特征線方法(MOC)具有強(qiáng)大的幾何處理能力和天然的并行特性，被廣泛應(yīng)用于反應(yīng)堆物理計(jì)算。隨著計(jì)算機(jī)計(jì)算能力的快速發(fā)展，越來(lái)越多的反應(yīng)堆物理計(jì)算程序采用MOC作為其主要的輸運(yùn)計(jì)算工具，如DeCART[1]、nTRACER[2]、MPACT[3]、OpenMOC[4]和NECP-X[5]等程序。但即使采用基于CPU的高性能計(jì)算集群，堆芯規(guī)模的MOC輸運(yùn)計(jì)算仍然耗費(fèi)大量的計(jì)算時(shí)間和計(jì)算資源。

近年來(lái)，基于GPU和CPU的異構(gòu)系統(tǒng)得到了快速發(fā)展[6]。由于異構(gòu)系統(tǒng)強(qiáng)大的浮點(diǎn)運(yùn)算能力，借助英偉達(dá)(NVIDIA)推出的統(tǒng)一計(jì)算設(shè)備構(gòu)架(CUDA)，CPU-GPU異構(gòu)系統(tǒng)被廣泛應(yīng)用于科學(xué)計(jì)算領(lǐng)域。同時(shí)，在MOC輸運(yùn)計(jì)算領(lǐng)域，已開展了采用GPU加速M(fèi)OC計(jì)算[7-11]及小規(guī)模的多GPU并行計(jì)算[12-13]相關(guān)研究。然而，目前基于CPU-GPU異構(gòu)系統(tǒng)，針對(duì)MOC的大規(guī)模CPU-GPU協(xié)同計(jì)算研究較少。

本文采用性能分析模型，分析二維MOC異構(gòu)并行算法的并行效率和性能瓶頸。同時(shí)針對(duì)顯著影響異構(gòu)并行效率的因素，采用輸運(yùn)計(jì)算與數(shù)據(jù)傳遞相互掩蓋的技術(shù)，提升異構(gòu)并行算法的整體性能和并行效率。

1 特征線方法

MOC是求解中子輸運(yùn)方程的一種數(shù)值方法，其求解思路是采用坐標(biāo)變換，將偏微分方程轉(zhuǎn)化為常微分方程。基于平源近似，省略能群標(biāo)識(shí)，特征線形式的穩(wěn)態(tài)多群中子輸運(yùn)方程為：

(1)

其中：φ(s)為中子角通量密度；s為中子穿行方向；Σt為總截面；Q為源項(xiàng)。

(2)

(3)

實(shí)際計(jì)算中，若干平行射線沿角度Ω穿過平源區(qū)n，引入角度離散p，采用數(shù)值積分，可求得平源區(qū)n的平均中子標(biāo)通量密度φ為：

(4)

其中：ωp為求積權(quán)重；dp為Ωp方向的射線間距。

2 性能分析和效率優(yōu)化

2.1 異構(gòu)并行算法

結(jié)合目前CPU-GPU異構(gòu)并行系統(tǒng)的特點(diǎn)，假設(shè)異構(gòu)系統(tǒng)由N個(gè)異構(gòu)節(jié)點(diǎn)組成，每個(gè)節(jié)點(diǎn)包含Nc個(gè)CPU核心及Ng個(gè)GPU，且通常Nc大于Ng。為充分利用異構(gòu)系統(tǒng)的計(jì)算資源，采用MPI+OpenMP/CUDA的混合編程模型來(lái)實(shí)現(xiàn)二維MOC異構(gòu)并行計(jì)算?；谠摦悩?gòu)模型，采用空間區(qū)域分解和MPI規(guī)范，實(shí)現(xiàn)粗粒度并行；而計(jì)算耗時(shí)的輸運(yùn)掃描過程(細(xì)粒度并行)，則通過基于OpenMP/CUDA的特征線并行來(lái)實(shí)現(xiàn)。

圖1示出基于空間區(qū)域分解和特征線并行的MPI+OpenMP/CUDA混合編程模型。假設(shè)單異構(gòu)節(jié)點(diǎn)包含2個(gè)GPU(Ng=2)及14個(gè)CPU核心(Nc=14)。在單異構(gòu)節(jié)點(diǎn)內(nèi)，待求解區(qū)域被分成4個(gè)子區(qū)，對(duì)應(yīng)啟動(dòng)Ng+Nx(Nx為采用MPI+OpenMP模型計(jì)算的進(jìn)程數(shù)，Nx=2)個(gè)MPI進(jìn)程，將每個(gè)子區(qū)分配給各進(jìn)程進(jìn)行計(jì)算。其中Ng個(gè)進(jìn)程由Ng個(gè)CPU核心啟動(dòng)，每個(gè)進(jìn)程通過CUDA環(huán)境驅(qū)動(dòng)1個(gè)GPU，采用MPI+CUDA模型進(jìn)行計(jì)算；另外Nx個(gè)進(jìn)程調(diào)用剩余的CPU核心，啟動(dòng)Nc-Ng(12)個(gè)OpenMP線程，采用MPI+OpenMP模型進(jìn)行計(jì)算，如圖1所示，可為每個(gè)進(jìn)程分配6個(gè)CPU核心參與計(jì)算。

在MPI+CUDA模型中，采用GPU線程并行實(shí)現(xiàn)細(xì)粒度并行。根據(jù)文獻(xiàn)[14]研究結(jié)果，為獲得最優(yōu)并行性能，單個(gè)GPU線程執(zhí)行單條特征線、單個(gè)能群的相關(guān)計(jì)算。如圖1所示，GPU線程t執(zhí)行特征線m第1個(gè)能群的相關(guān)計(jì)算，線程t+1則執(zhí)行特征線m第2個(gè)能群的相關(guān)計(jì)算，依此類推。而對(duì)于MPI+OpenMP模型，為減少線程開辟和銷毀帶來(lái)的額外開銷，單次外迭代對(duì)應(yīng)進(jìn)行1次OpenMP線程的開辟和銷毀。相應(yīng)的計(jì)算任務(wù)顯式地分配給各線程，且單個(gè)線程負(fù)責(zé)特征線的所有能群和極角的輸運(yùn)計(jì)算。如圖1所示，CPU線程s首先執(zhí)行特征線n的相關(guān)計(jì)算，計(jì)算完成后，特征線n+線程數(shù)(線程數(shù)指的是當(dāng)前進(jìn)程中參與計(jì)算的OpenMP線程數(shù))對(duì)應(yīng)的計(jì)算被分配給線程n，依此類推。

由特征線并行的實(shí)現(xiàn)方式可知，本文采用了Jacobi格式的輸運(yùn)掃描算法，減弱了能群間的耦合關(guān)系，相對(duì)于Gauss-Seidel格式的輸運(yùn)掃描算法，會(huì)增加輸運(yùn)計(jì)算外迭代次數(shù)。對(duì)于二維C5G7 1/4堆芯基準(zhǔn)題，外迭代次數(shù)增量約為12.8%。但文獻(xiàn)[14]研究表明，GPU并行計(jì)算需通過設(shè)置大量并行任務(wù)來(lái)掩蓋GPU訪存延遲，而同時(shí)將特征線和能群并行，能夠保證足夠多的并行任務(wù)，獲得最優(yōu)的計(jì)算性能。因此基于文獻(xiàn)[14]的研究結(jié)果，在CPU-GPU異構(gòu)并行計(jì)算中依然采用Jacobi格式的輸運(yùn)掃描算法。

圖1 基于空間區(qū)域分解和特征線并行的混合MPI+OpenMP/CUDA編程模型Fig.1 Hybrid MPI+OpenMP/CUDA programming model with spatial domain decomposition and ray parallelization

為保證異構(gòu)并行計(jì)算中CPU與GPU計(jì)算之間的負(fù)載平衡，本文采用文獻(xiàn)[15]所提出的動(dòng)態(tài)任務(wù)分配模型，根據(jù)硬件的實(shí)際計(jì)算性能實(shí)現(xiàn)最優(yōu)任務(wù)分配。在該模型中，總計(jì)算任務(wù)被任意初始分配給MPI進(jìn)程；隨后執(zhí)行預(yù)處理過程，進(jìn)行若干次輸運(yùn)掃描，并獲取各進(jìn)程輸運(yùn)計(jì)算時(shí)間；假設(shè)每個(gè)MPI進(jìn)程計(jì)算能力在整個(gè)計(jì)算過程中保持為常數(shù)，則進(jìn)程的輸運(yùn)計(jì)算時(shí)間與所分配的任務(wù)量呈正比，可采用預(yù)處理過程中初始分配的計(jì)算任務(wù)和獲取的計(jì)算時(shí)間，求解最優(yōu)任務(wù)分配方案，保證各進(jìn)程之間的負(fù)載平衡。

2.2 計(jì)算流程

圖2示出二維MOC異構(gòu)并行計(jì)算流程：1) 并行初始化過程中執(zhí)行空間區(qū)域分解，程序?qū)⒂?jì)算任務(wù)分配到p個(gè)MPI進(jìn)程；2) CPU端執(zhí)行本地初始化，生成截面、特征線等信息；3) 對(duì)于采用MPI+CUDA并行的進(jìn)程(進(jìn)程#1)，將計(jì)算所需的信息(截面、特征線信息等)從CPU一次性拷貝到GPU，隨后在GPU端執(zhí)行源項(xiàng)計(jì)算、輸運(yùn)掃描和反應(yīng)率及通量殘差計(jì)算，再將反應(yīng)率、通量殘差及邊界角通量密度從GPU通過CPU拷貝到系統(tǒng)內(nèi)存；4) 對(duì)于采用MPI+OpenMP并行的進(jìn)程(進(jìn)程#2)，則依次在CPU進(jìn)行源項(xiàng)計(jì)算、輸運(yùn)掃描和反應(yīng)率及殘差計(jì)算，其中輸運(yùn)掃描過程采用OpenMP并行實(shí)現(xiàn)；5) 空間子區(qū)邊界角通量密度交換由MPI通信實(shí)現(xiàn)；6) 在主進(jìn)程計(jì)算keff并進(jìn)行收斂判斷。

圖2 二維MOC異構(gòu)并行計(jì)算流程Fig.2 Flowchart of heterogeneous parallelization of 2D MOC

2.3 性能分析模型

結(jié)合圖2所示的計(jì)算流程可知，對(duì)于異構(gòu)并行計(jì)算，MPI+CUDA模式是影響異構(gòu)并行性能的主要方面。因此，識(shí)別影響MPI+CUDA并行效率的因素并進(jìn)行相應(yīng)的優(yōu)化是提升異構(gòu)并行效率的主要途徑。而影響MPI+CUDA并行效率的主要因素包括：1) 由于引入空間區(qū)域分解，邊界角通量密度更新滯后，導(dǎo)致迭代次數(shù)增加；2) 負(fù)載不均衡導(dǎo)致的性能下降；3) CPU和GPU間的數(shù)據(jù)拷貝過程對(duì)性能的影響；4) MPI通信對(duì)并行效率的影響。本文采用4個(gè)指標(biāo)量化上述因素對(duì)并行效率的影響。程序計(jì)算時(shí)間由式(5)描述。

Ttotal(p)=T(p)+Tdata_copy(p)+TMPI_comm(p)

(5)

其中：p為進(jìn)程編號(hào)；Ttotal為各進(jìn)程整體計(jì)算時(shí)間；Tdata_copy為CPU與GPU間數(shù)據(jù)拷貝時(shí)間；TMPI_comm為MPI通信時(shí)間；T為各進(jìn)程除數(shù)據(jù)傳遞(MPI通信和CPU與GPU間數(shù)據(jù)拷貝)之外的計(jì)算時(shí)間。

則上述4個(gè)影響因素可由式(6)～(9)量化：

(6)

(7)

(8)

(9)

其中：Np為并行進(jìn)程數(shù)；Tserial為串行計(jì)算時(shí)間。

因此，總并行效率則為上述4個(gè)指標(biāo)的乘積，即：

(10)

2.4 輸運(yùn)計(jì)算和數(shù)據(jù)傳遞相互掩蓋

考慮到圖2所示的數(shù)據(jù)傳遞(CPU與GPU間的數(shù)據(jù)拷貝及MPI通信)顯著影響程序并行效率，本文采用CUDA流、輸運(yùn)計(jì)算和MPI通信異步執(zhí)行等方案來(lái)實(shí)現(xiàn)輸運(yùn)計(jì)算掩蓋數(shù)據(jù)傳遞。對(duì)于MPI+CUDA并行，本文采用CUDA流技術(shù)，將CPU與GPU間的數(shù)據(jù)拷貝和輸運(yùn)計(jì)算分成若干部分，各對(duì)應(yīng)部分組成CUDA流，將“分塊”的數(shù)據(jù)拷貝和輸運(yùn)計(jì)算重疊執(zhí)行，如圖3所示。同時(shí)，考慮到CUDA函數(shù)為異步執(zhí)行，即CPU端調(diào)用CUDA函數(shù)后，該函數(shù)在GPU端執(zhí)行，隨后計(jì)算流程控制權(quán)返回CPU。因此在GPU計(jì)算期間，與之綁定的CPU核心可執(zhí)行CPU與GPU間的數(shù)據(jù)拷貝及MPI通信，實(shí)現(xiàn)數(shù)據(jù)傳遞與GPU計(jì)算的相互掩蓋。對(duì)于MPI+OpenMP并行，則從原來(lái)執(zhí)行輸運(yùn)計(jì)算的t個(gè)線程中分出1個(gè)線程，專門執(zhí)行MPI通信，另外的t-1個(gè)線程依然執(zhí)行輸運(yùn)計(jì)算，實(shí)現(xiàn)MPI通信與輸運(yùn)計(jì)算的相互掩蓋。相比原執(zhí)行方式，實(shí)現(xiàn)輸運(yùn)計(jì)算掩蓋數(shù)據(jù)傳遞之后，可提高計(jì)算效率，如圖3所示標(biāo)注了計(jì)算時(shí)間節(jié)省。

圖3 輸運(yùn)計(jì)算掩蓋數(shù)據(jù)傳遞示意圖Fig.3 Scheme of transport sweep and data movement overlapping

3 數(shù)值驗(yàn)證

3.1 精度驗(yàn)證

本文采用二維C5G7基準(zhǔn)題[16]驗(yàn)證異構(gòu)并行程序計(jì)算精度，圖4示出該基準(zhǔn)題的1/4堆芯布置，每個(gè)燃料組件為17×17柵元排列。整個(gè)計(jì)算區(qū)域被劃分為51×51個(gè)柵元，每個(gè)柵元采用如圖5所示的網(wǎng)格劃分。特征線初始間距為0.03 cm，每個(gè)象限選擇14個(gè)方位角和3個(gè)極角，使用Tabuchi-Yamamoto極角求積組。keff和通量收斂準(zhǔn)則分別設(shè)置為10-6和10-5。計(jì)算平臺(tái)采用天津超算中心TianHe-ES異構(gòu)計(jì)算系統(tǒng)，單個(gè)節(jié)點(diǎn)硬件參數(shù)列于表1。

本征值和歸一化功率分布誤差結(jié)果列于表2，參考解來(lái)源于MCNP。表2中：MAX Error為最大棒功率相對(duì)誤差；AVG Error為平均棒功率相對(duì)誤差；RMS Error為棒功率均方根誤差；MRE Error為平均相對(duì)誤差。由表2可見，異構(gòu)并行計(jì)算得到的本征值和功率分布誤差在可接受范圍內(nèi)，除最小柵元功率外，其他結(jié)果相對(duì)誤差均在MCNP統(tǒng)計(jì)不確定度范圍內(nèi)。結(jié)果表明，本文所提出的并行算法具備良好的計(jì)算精度。

圖4 二維C5G7基準(zhǔn)題1/4堆芯布置Fig.4 Quarter core configuration of 2D C5G7 benchmark problem

圖5 柵元組成及網(wǎng)格劃分Fig.5 Fuel pin layout and mesh

項(xiàng)目參數(shù)規(guī)格CPU2×Intel XeonCPU E5-1690 v4,主頻2.60 GHz,14核GPU2×Tesla K80,即4×GK 210,主頻824 MHz,計(jì)算能力3.7操作系統(tǒng)RedHat 7.2編譯環(huán)境Intel MPI5.1.3,CUDA 9.2

3.2 異構(gòu)并行算法效率分析

為分析異構(gòu)并行算法的計(jì)算效率和性能瓶頸，將二維C5G7基準(zhǔn)題進(jìn)行擴(kuò)展，構(gòu)造10×10燃料組件的全堆芯問題，堆芯布置如圖6所示。采用本文提出的性能分析模型測(cè)試了MPI+CUDA并行算法在采用輸運(yùn)計(jì)算掩蓋數(shù)據(jù)傳遞前的強(qiáng)并行效率及算法各部分對(duì)并行效率的影響，如圖7所示。由圖7可知：1) 隨并行進(jìn)程數(shù)的增多，MPI+CUDA并行算法強(qiáng)并行效率逐步降低，當(dāng)采用20進(jìn)程并行時(shí)，算法強(qiáng)并行效率約為78%；2) MPI通信(EMPI_comm)是造成MPI+CUDA并行算法并行效率低下的主要因素，而CPU和GPU間的數(shù)據(jù)拷貝(Edata_copy)是影響并行效率的另一因素；3) 相對(duì)而言，迭代次數(shù)增加(Eincr_iter)對(duì)整體并行效率的影響較小，且由于圖7所示的并行算例基本保持負(fù)載均衡，因此負(fù)載不均衡(Eload_imb)對(duì)并行效率的影響相對(duì)較小。圖8示出MPI通信和CPU與GPU間的數(shù)據(jù)拷貝在MPI+CUDA并行中的時(shí)間占比，在20個(gè)進(jìn)程并行時(shí)，這兩部分時(shí)間占比分別達(dá)到了11%和5%。可見，提升異構(gòu)并行算法并行效率的主要途徑為減小數(shù)據(jù)傳遞過程(MPI通信和CPU與GPU間的數(shù)據(jù)拷貝)對(duì)并行效率的影響。

表2 本征值及功率分布的結(jié)果Table 2 Result of Eigenvalue and power distribution

圖6 虛擬堆芯組件布置Fig.6 Arrangement of mock core assembly

圖7 MPI+CUDA并行算法強(qiáng)并行效率Fig.7 Strong scaling efficiency for MPI+CUDA parallel algorithm

圖8 MPI+CUDA并行中的時(shí)間占比Fig.8 Proportion of execution time in MPI+CUDA parallelization

為提升異構(gòu)并行算法的并行效率，本文實(shí)現(xiàn)了輸運(yùn)計(jì)算與數(shù)據(jù)傳遞相互掩蓋。針對(duì)圖6所示的全堆芯問題，采用1～5個(gè)異構(gòu)節(jié)點(diǎn)，測(cè)試了異構(gòu)并行算法采用輸運(yùn)計(jì)算掩蓋數(shù)據(jù)傳遞前后的效率對(duì)比。表3列出程序計(jì)算收斂所需的迭代次數(shù)。圖9示出程序整體計(jì)算效率。由表3可知，由于MPI通信和輸運(yùn)計(jì)算同時(shí)執(zhí)行，造成迭代格式退化，迭代次數(shù)增多，迭代次數(shù)增量在3%左右。但由于數(shù)據(jù)傳遞時(shí)間在異構(gòu)并行計(jì)算中占有顯著比例，由圖9可知，實(shí)現(xiàn)輸運(yùn)計(jì)算與數(shù)據(jù)傳遞相互掩蓋后，程序整體效率有所提升。采用5節(jié)點(diǎn)計(jì)算時(shí)，程序整體效率提升達(dá)到8%。

表3 輸運(yùn)計(jì)算掩蓋數(shù)據(jù)傳遞前后迭代次數(shù)對(duì)比Table 3 Comparison of iteration number with and without transport sweep and data movement overlapping

圖9 程序整體計(jì)算效率Fig.9 Overall calculating efficiency of program

圖10示出MPI+OpenMP/CUDA異構(gòu)并行模式下的以計(jì)算節(jié)點(diǎn)為參考的強(qiáng)并行效率(單節(jié)點(diǎn)計(jì)算效率為100%)。由圖10可知，未采用輸運(yùn)計(jì)算掩蓋數(shù)據(jù)傳遞時(shí)，強(qiáng)并行效率隨節(jié)點(diǎn)數(shù)的增加下降較快，這是由于在這種情況下，MPI通信和CPU與GPU間的數(shù)據(jù)拷貝會(huì)顯著影響程序的整體并行效率。而采用輸運(yùn)計(jì)算掩蓋數(shù)據(jù)傳遞后，強(qiáng)并行效率有明顯提升，在5節(jié)點(diǎn)并行時(shí)，強(qiáng)并行效率達(dá)到95%。

針對(duì)圖6所示的問題，采用多CPU節(jié)點(diǎn)進(jìn)行并行計(jì)算，其中CPU節(jié)點(diǎn)與異構(gòu)節(jié)點(diǎn)除GPU以外的部分配置相同。圖11示出CPU-GPU異構(gòu)節(jié)點(diǎn)和CPU節(jié)點(diǎn)的性能對(duì)比。其中，異構(gòu)并行計(jì)算基于MPI+OpenMP/CUDA(輸運(yùn)計(jì)算掩蓋數(shù)據(jù)傳遞)并行模式完成；基于MPI的并行計(jì)算則分別采用5、10、15、20個(gè)CPU節(jié)點(diǎn)并行執(zhí)行。由圖11可知：CPU-GPU異構(gòu)節(jié)點(diǎn)計(jì)算性能明顯優(yōu)于CPU節(jié)點(diǎn)，4個(gè)CPU-GPU異構(gòu)節(jié)點(diǎn)計(jì)算能力優(yōu)于20個(gè)CPU節(jié)點(diǎn)計(jì)算能力。

圖10 MPI+OpenMP/CUDA強(qiáng)并行效率對(duì)比Fig.10 Strong scaling efficiency comparison of MPI+OpenMP/CUDA

圖11 CPU-GPU異構(gòu)節(jié)點(diǎn)與CPU節(jié)點(diǎn)計(jì)算效率對(duì)比Fig.11 Calculating efficiency comparison of CPU-GPU heterogeneous node and CPU node

4 結(jié)論

本文基于空間區(qū)域分解和特征線并行，采用MPI+OpenMP/CUDA編程模型，實(shí)現(xiàn)了適用于CPU-GPU異構(gòu)系統(tǒng)的二維MOC異構(gòu)并行算法。同時(shí)，采用性能分析模型，針對(duì)影響異構(gòu)并行算法并行效率的影響因素進(jìn)行了分析，并針對(duì)性地進(jìn)行了優(yōu)化，實(shí)現(xiàn)了輸運(yùn)計(jì)算和數(shù)據(jù)傳遞的相互掩蓋，緩解了MPI通信和CPU與GPU間的數(shù)據(jù)傳遞對(duì)程序整體計(jì)算效率的影響。數(shù)值結(jié)果表明：程序具備良好的計(jì)算精度；MPI+OpenMP/CUDA異構(gòu)并行模式下，數(shù)據(jù)傳遞(MPI通信和CPU與GPU間的數(shù)據(jù)拷貝)是影響并行效率的主要因素；實(shí)現(xiàn)輸運(yùn)計(jì)算與數(shù)據(jù)傳遞相互掩蓋能有效提升程序整體性能；5節(jié)點(diǎn)(包含20個(gè)GPU)并行時(shí)，程序整體效率提升達(dá)到8%，強(qiáng)并行效率達(dá)到95%。同時(shí)，CPU-GPU異構(gòu)節(jié)點(diǎn)性能明顯優(yōu)于CPU節(jié)點(diǎn)，4異構(gòu)節(jié)點(diǎn)計(jì)算能力優(yōu)于20個(gè)CPU節(jié)點(diǎn)。未來(lái)工作將側(cè)重于基于大規(guī)模CPU-GPU異構(gòu)系統(tǒng)的粗網(wǎng)有限差分(CMFD)加速算法研究。