基于CPU-GPU異構并行的MOC中子輸運計算并行效率優化研究

宋佩濤，張志儉，梁 亮，張 乾，趙 強

(哈爾濱工程大學 核安全與仿真技術國防重點學科實驗室，黑龍江 哈爾濱 150001)

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

近年來，基于GPU和CPU的異構系統得到了快速發展[6]。由于異構系統強大的浮點運算能力，借助英偉達(NVIDIA)推出的統一計算設備構架(CUDA)，CPU-GPU異構系統被廣泛應用于科學計算領域。同時，在MOC輸運計算領域，已開展了采用GPU加速MOC計算[7-11]及小規模的多GPU并行計算[12-13]相關研究。然而，目前基于CPU-GPU異構系統，針對MOC的大規模CPU-GPU協同計算研究較少。

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

1 特征線方法

MOC是求解中子輸運方程的一種數值方法，其求解思路是采用坐標變換，將偏微分方程轉化為常微分方程?；谄皆唇?，省略能群標識，特征線形式的穩態多群中子輸運方程為：

(1)

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

(2)

(3)

實際計算中，若干平行射線沿角度Ω穿過平源區n，引入角度離散p，采用數值積分，可求得平源區n的平均中子標通量密度φ為：

(4)

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

2 性能分析和效率優化

2.1 異構并行算法

結合目前CPU-GPU異構并行系統的特點，假設異構系統由N個異構節點組成，每個節點包含Nc個CPU核心及Ng個GPU，且通常Nc大于Ng。為充分利用異構系統的計算資源，采用MPI+OpenMP/CUDA的混合編程模型來實現二維MOC異構并行計算?；谠摦悩嬆Ｐ?，采用空間區域分解和MPI規范，實現粗粒度并行；而計算耗時的輸運掃描過程(細粒度并行)，則通過基于OpenMP/CUDA的特征線并行來實現。

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

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

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

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

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

2.2 計算流程

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

圖2 二維MOC異構并行計算流程Fig.2 Flowchart of heterogeneous parallelization of 2D MOC

2.3 性能分析模型

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

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

(5)

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

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

(6)

(7)

(8)

(9)

其中：Np為并行進程數；Tserial為串行計算時間。

因此，總并行效率則為上述4個指標的乘積，即：

(10)

2.4 輸運計算和數據傳遞相互掩蓋

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

圖3 輸運計算掩蓋數據傳遞示意圖Fig.3 Scheme of transport sweep and data movement overlapping

3 數值驗證

3.1 精度驗證

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

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

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

圖5 柵元組成及網格劃分Fig.5 Fuel pin layout and mesh

項目參數規格CPU2×Intel XeonCPU E5-1690 v4,主頻2.60 GHz,14核GPU2×Tesla K80,即4×GK 210,主頻824 MHz,計算能力3.7操作系統RedHat 7.2編譯環境Intel MPI5.1.3,CUDA 9.2

3.2 異構并行算法效率分析

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

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

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

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

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

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

表3 輸運計算掩蓋數據傳遞前后迭代次數對比Table 3 Comparison of iteration number with and without transport sweep and data movement overlapping

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

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

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

圖10 MPI+OpenMP/CUDA強并行效率對比Fig.10 Strong scaling efficiency comparison of MPI+OpenMP/CUDA

圖11 CPU-GPU異構節點與CPU節點計算效率對比Fig.11 Calculating efficiency comparison of CPU-GPU heterogeneous node and CPU node

4 結論

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