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

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

特征線方法(MOC)因其幾何適應(yīng)性和天然的并行能力，被廣泛應(yīng)用于反應(yīng)堆物理數(shù)值計算，眾多知名堆芯中子學(xué)程序均采用了MOC，如CRX[1]、DeCART[2]、nTRACER[3]、MPACT[4]、OpenMOC[5]和NECP-X[6]。但堆芯規(guī)模的MOC輸運計算耗費大量的計算機資源和計算時間，因此有效的加速手段和高效并行算法是全堆芯MOC計算的重點研究方向。

近年來，圖形處理器(GPU)在制造工藝、功耗等方面不斷發(fā)展，借助英偉達(NVIDIA)推出的統(tǒng)一計算設(shè)備構(gòu)架(CUDA)[7]，GPU功能也逐步從圖像處理拓展到了科學(xué)計算領(lǐng)域。相比CPU，GPU兼顧通用性、效能和性能等多方面因素，具有更高的浮點運算能力和顯存帶寬。因此，GPU為加速全堆芯MOC計算提供了新方法和思路。

然而，目前采用GPU加速MOC中子輸運計算的研究較少。國外Boyd等[8]、Tramm等[9]、Chio等[10]先后實現(xiàn)了GPU加速二維和三維MOC中子輸運計算。國內(nèi)韓宇[11]和鄭勇[12]分別實現(xiàn)了GPU加速六角形幾何MOC計算及模塊化MOC計算。異構(gòu)并行方面，張知竹[13]實現(xiàn)了基于多GPU系統(tǒng)的三維MOC計算。以上研究大多側(cè)重于實現(xiàn)GPU對MOC的加速，而進一步的性能分析相對較少。

本文針對二維MOC輸運計算，采用CUDA編程規(guī)范，實現(xiàn)基于NVIDIA GPU加速MOC計算的并行算法。

1 MOC輸運計算

基于平源近似，省略能群標識，特征線形式的穩(wěn)態(tài)多群中子輸運方程為：

(1)

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

對于式(1)，在平源區(qū)內(nèi)沿s進行積分可得：

(2)

引入積分因子，對式(1)進行解析求解可得步特征線法：

(3)

其中，式(3)中指數(shù)項(1-e-Σtls)可通過指數(shù)函數(shù)和線性插值兩種方式進行計算。

采用菱形差分[14]，式(2)中的平均角通量可表示為：

(4)

聯(lián)立式(2)、(4)可得平均角通量：

(5)

實際計算中，有1組平行射線沿角度Ω穿過平源區(qū)n，引入角度離散p，采用數(shù)值積分，可求得平源區(qū)n的平均中子標通量：

(6)

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

2 算法實現(xiàn)

2.1 迭代流程

GPU加速MOC求解流程如圖1所示，其基本思路是將計算密集的輸運掃描部分移植到GPU并行計算，而計算參數(shù)準備和收斂判斷等涉及計算流程控制以及邏輯性強的過程則在CPU進行。計算流程為：1) 在CPU進行計算初始化，生成截面、特征線等信息；2) 將計算所需的信息(截面信息、特征線信息等)從CPU一次性拷貝到GPU；3) 在GPU執(zhí)行平源區(qū)源項計算、輸運掃描和裂變率及通量殘差計算；4) 將裂變率和通量殘差從GPU拷貝到CPU；5) 在CPU計算keff，根據(jù)前后兩次迭代的keff相對變化率和通量殘差進行收斂判斷。在CPU進行keff計算和收斂判斷，相比在GPU進行此過程，多了1個數(shù)據(jù)拷貝過程(過程4)，該數(shù)據(jù)拷貝過程時間占比較小(<計算總時間的0.1%)，因此對程序整體效率影響較小。本文之所以在CPU進行keff計算和收斂判斷，是為了便于后續(xù)研究中多CPU和多GPU異構(gòu)并行的實現(xiàn)。在多GPU異構(gòu)并行中，各GPU只負責(zé)一部分輸運計算任務(wù)，全局keff的求解和收斂判斷，則需CPU將各GPU求得的裂變率和局部通量殘差收集后在CPU進行。

圖1 GPU加速MOC求解流程Fig.1 Flow diagram of MOC on GPU

輸運掃描部分為整個求解流程的重點，本文采用Jacobi MOC掃描算法進行輸運求解，即在輸運掃描前更新所有空間區(qū)域和所有能群的散射源項和裂變源項，而在單次輸運掃描過程中，源項不進行更新?？紤]到MOC計算在特征線和能群層次的并行特性，本文所采用的并行策略為特征線加能群并行方式，即單個GPU線程負責(zé)單條特征線的單個能群的相關(guān)計算。

在CUDA中，將CPU端調(diào)用、GPU端執(zhí)行的函數(shù)稱為核函數(shù)，核函數(shù)在GPU線程上執(zhí)行。GPU線程以線程束形式參與計算，1個線程束中包含若干線程，同一線程束中的線程執(zhí)行相同的指令。線程束的大小由硬件的計算能力版本決定，對于本文所采用的GPU，線程束大小為32，即32條線程組成1個線程束，執(zhí)行相同的指令。本文所涉及的輸運掃描過程在輸運掃描核函數(shù)中進行。由于輸運掃描計算在整個計算過程中最為耗時(時間占比90%以上)，因此輸運掃描核函數(shù)的執(zhí)行性能直接影響程序的整體效率，本文圍繞輸運掃描核函數(shù)的性能進行分析。

2.2 輸運掃描核函數(shù)

輸運掃描核函數(shù)負責(zé)源迭代過程中的輸運掃描計算，其核心思想是將沿某方向貫穿整個求解區(qū)域的1條特征線的1個能群的計算任務(wù)分配給1個GPU線程，算法過程為：1) 為當(dāng)前GPU線程分配計算任務(wù)，即根據(jù)當(dāng)前線程ID(_tID)確定待求解的特征線ID(_trackID)和能群ID(_groupID)；2) 為當(dāng)前GPU線程分配共享內(nèi)存空間；3) 沿當(dāng)前特征線進行輸運掃描。GPU以線程束為單位進行計算，線程束的執(zhí)行由線程束調(diào)度器自動分配，以保證硬件資源的充分利用。因此，開發(fā)人員應(yīng)關(guān)注線程束內(nèi)部的負載均衡。該算法將1條特征線、1個能群的相關(guān)計算分配給1個GPU線程。當(dāng)1個線程束內(nèi)部所有線程都在執(zhí)行同一條特征線的不同能群的計算時，該線程束能保證負載均衡；當(dāng)1個線程束內(nèi)部涉及不同特征線的計算時，由于不同特征線中線段數(shù)量不同，此時就會引入負載不均衡。這種負載不均衡在當(dāng)前算法下不可避免，但由于單個線程計算量本身較小，這種線程束內(nèi)部的負載不均衡對整體計算效率的影響得到了一定程度的緩解。

輸運掃描過程中，當(dāng)前特征線段出/入射點間角通量的變化量采用步特征線法和菱形差分兩種方法求解，而針對步特征線法中指數(shù)項的求解，分別采用指數(shù)函數(shù)計算和線性插值計算。

GPU同時提供了強大的單精度和雙精度浮點運算能力。因此在輸運掃描核函數(shù)的實現(xiàn)過程中，考慮了不同浮點運算精度對程序計算精度和效率的影響。結(jié)合GPU硬件特性，在輸運掃描核函數(shù)實現(xiàn)過程中采用了如下幾項處理方式。

1) 使用GPU共享內(nèi)存存儲角通量信息

輸運掃描過程中，特征線段出射角通量隨計算過程不斷更新，對應(yīng)角通量數(shù)組頻繁的讀/寫操作。從硬件角度來講，相比GPU顯存，作為片上存儲器的共享內(nèi)存具有更高的訪存帶寬和更小的訪存延遲。因此，將角通量數(shù)組存儲在共享內(nèi)存。

2) 采用不同精度進行原子操作

平源區(qū)特定能群的標通量按照數(shù)值積分過程由不同方向上特征線段的角通量加權(quán)求和獲得，該過程涉及不同線程向同一全局內(nèi)存地址寫入數(shù)據(jù)。為保證數(shù)據(jù)寫入的準確性，CUDA提供了讀取-修改-寫入的原子操作。原子操作定義為：不會被線程調(diào)度機制打斷的操作。1個原子操作過程包括3部分：從某個地址讀取數(shù)據(jù)；將該數(shù)據(jù)加某一特定值；將修改后的數(shù)據(jù)寫回原地址。原子操作保證了上述3步操作過程不會被其他線程打斷，即相關(guān)線程順序訪問同一內(nèi)存地址。CUDA在雙精度和單精度下均提供了原子操作功能，且不同浮點運算精度下訪存帶寬不同。因此可在整體采用雙精度的情況下，采用單精度聲明標通量數(shù)組scalarFlux[]，并在原子操作部分采用單精度，即采用混合精度進行計算。

3) 采用不同精度進行指數(shù)函數(shù)計算

對于步特征線-指數(shù)函數(shù)核函數(shù)，其計算過程中涉及指數(shù)函數(shù)計算。CUDA在雙精度和單精度情況下提供了不同的指數(shù)函數(shù)。相對雙精度，單精度指數(shù)函數(shù)計算更快，但精度略有損失。因此，在步特征線-指數(shù)函數(shù)核函數(shù)中，整體采用雙精度計算時，指數(shù)函數(shù)可采用單精度計算，即采用混合精度方式計算。

本文針對菱形差分、步特征線-指數(shù)函數(shù)和步特征線-線性插值3種核函數(shù)，對應(yīng)雙精度、混合精度及單精度分別測試不同組合方式對程序計算精度和計算效率的影響。表1列出不同核函數(shù)和不同精度的組合說明。需要注意的是，為保證特征線信息的準確性，在特征線信息生成過程中采用雙精度進行生成，而對于整體采用單精度的計算過程，特征線信息存儲時采用單精度進行截斷。

3 性能分析

NVIDIA提供了CUDA程序性能分析工具集，便于用戶對GPU程序進行調(diào)試及性能分析。為量化當(dāng)前程序?qū)PU硬件的使用程度并識別當(dāng)前程序的性能瓶頸，本文借助nvprof性能分析工具，對3種輸運掃描核函數(shù)在不同精度下的性能進行分析，性能分析指標列于表2。

分析線程束指令簽發(fā)阻斷原因可輔助識別程序的性能瓶頸。對于每個時鐘周期，線程束調(diào)度器都會嘗試向線程束簽發(fā)1條指令，當(dāng)線程束由于某種原因處于不可用狀態(tài)時，指令簽發(fā)則會被阻斷。造成指令簽發(fā)阻斷的因素主要包括訪存相關(guān)和執(zhí)行相關(guān)。訪存相關(guān)指待簽發(fā)指令正在等待訪存完成，造成簽發(fā)阻斷；執(zhí)行相關(guān)指待簽發(fā)指令正在等待前序指令執(zhí)行完成，以提供待簽發(fā)指令執(zhí)行所需的參數(shù)，造成簽發(fā)阻斷。

4 數(shù)值計算

本文采用二維C5G7基準題[15]進行測試與分析，圖2示出該基準題的1/4堆芯和燃料組件布置。整個計算區(qū)域在x和y方向分別被劃分為51個柵元，每個柵元采用如圖3所示的網(wǎng)格劃分。每個象限選擇14個方位角和3個極角，使用Tabuchi-Yamamoto極角求積組。

表2 性能分析指標Table 2 Metrics for performance analysis

圖2 二維C5G7基準題1/4堆芯布置(a)和燃料組件布置(b)Fig.2 Core configuration (a) and fuel pin composition (b) for 2D C5G7 benchmark problem

keff和通量收斂準則分別設(shè)置為10-6和10-5。計算平臺采用多核工作站，硬件及相關(guān)參數(shù)列于表3。

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

表3 計算平臺參數(shù)Table 3 Parameter of computing platform

4.1 計算精度和效率測試

不同測試情況下，keff、柵元功率及組件功率誤差列于表4，參考解來源于MCNP。表4中，MAX Error為最大棒功率相對誤差，AVG Error為平均棒功率相對誤差，RMS Error為棒功率均方根誤差，MRE Error為平均相對誤差。對于菱形差分和步特征線法，由于方法不同，keff計算結(jié)果相差約8 pcm。由表4可見：本征值和功率分布誤差均在合理范圍內(nèi)；采用混合精度和單精度浮點運算，相對于同等計算條件下的雙精度浮點運算結(jié)果，由于精度損失引入的誤差較小，均低于收斂準則?？梢?，對于二維C5G7基準題，在輸運求解部分，不同的浮點運算精度均能滿足計算精度要求。

表4 二維C5G7基準題特征值及功率分布的結(jié)果Table 4 Result of eigenvalue and power distributions for 2D C5G7 benchmark problem

4.2 輸運掃描核函數(shù)性能分析

本文分析了3種類型輸運掃描核函數(shù)在3種浮點運算精度下的性能特點，識別了GPU硬件使用程度及程序的性能瓶頸。

1) 計算效率分析

圖4 二維C5G7基準題計算時間Fig.4 Runtime of 2D C5G7 benchmark problem

圖5 二維C5G7基準題浮點運算次數(shù)及浮點運算效率Fig.5 Count and efficiency of floating octal point for 2D C5G7 benchmark problem

圖4示出程序整體計算時間，包括輸運掃描核函數(shù)執(zhí)行時間和程序其余部分執(zhí)行時間，其中輸運掃描核函數(shù)占到總計算時間的90%以上。圖5示出程序浮點運算次數(shù)和浮點運算效率，其中柱狀圖表示浮點運算次數(shù)，其頂部數(shù)據(jù)表示浮點運算效率。圖6示出L1和L2緩存命中率。

圖6 二維C5G7基準題L1和L2緩存命中率Fig.6 L1 and L2 cache hit rates for 2D C5G7 benchmark problem

在雙精度情況下，由圖4可知，菱形差分和步特征線-線性插值計算效率相當(dāng)，而步特征線-指數(shù)函數(shù)計算最為耗時。這是由于在雙精度情況下，步特征線-指數(shù)函數(shù)中指數(shù)函數(shù)的計算具有較大的浮點運算量，因此，在如圖5所示浮點運算量中，步特征線-指數(shù)函數(shù)在雙精度情況下的浮點運算量要明顯多于其他核函數(shù)。在雙精度情況下，步特征線-線性插值性能最優(yōu)。

在混合精度情況下，輸運掃描核函數(shù)中標通量數(shù)據(jù)的讀/寫操作(其中寫為原子操作)均在單精度下完成，圖5所示的浮點運算次數(shù)相對于雙精度情況有所下降；同時，由于單精度訪存操作具有更高的訪存效率，因此相較于雙精度情況，L1和L2緩存命中率均有所提升(圖6)。運算量的減少和緩存命中率的提升共同導(dǎo)致了計算效率的提升，因此在圖4中可看到混合精度情況下的計算時間相比雙精度計算有所下降。此外，造成步特征線-指數(shù)函數(shù)性能明顯提升的另一關(guān)鍵因素是單精度形式指數(shù)函數(shù)的使用，這使得浮點運算次數(shù)大幅下降(圖6)，從而導(dǎo)致了計算效率的提升。結(jié)果表明：在混合精度情況下，相對于雙精度計算，引入單精度標通量數(shù)據(jù)讀/寫操作的加速效果約為1.07倍；而在步特征線-指數(shù)函數(shù)中引入單精度指數(shù)函數(shù)計算的加速效果較為明顯，約為2.6倍。在混合精度情況下，步特征線-指數(shù)函數(shù)性能最優(yōu)。

在單精度情況下，3種核函數(shù)的計算效率相對于雙精度和混合精度情況有顯著提升，步特征線-指數(shù)函數(shù)效率優(yōu)于另外兩種核函數(shù)。

2) 性能瓶頸分析

輸運掃描核函數(shù)在不同精度下表現(xiàn)出不同的性能瓶頸。圖7示出硬件浮點運算單元利用率及訪存效率，圖8示出造成線程束指令簽發(fā)阻斷的原因。

圖7 二維C5G7基準題浮點運算單元利用率及訪存效率Fig.7 Arithmetic function unit utilization and DRAM efficiency for 2D C5G7 benchmark problem

圖8 二維C5G7基準題線程束指令簽發(fā)阻斷原因Fig.8 Warp issue stall reason of 2D C5G7 benchmark problem

在雙精度和混合精度情況下，如圖7所示，硬件雙精度浮點運算單元利用率均達到100%。結(jié)合圖8可知，60%以上的線程束指令簽發(fā)阻斷由執(zhí)行相關(guān)因素引起，即待簽發(fā)指令多數(shù)情況下在等待前序指令執(zhí)行完成。同時分析圖5所示的浮點運算效率可知，在雙精度和混合精度情況下，3種類型核函數(shù)均達到了較高的雙精度浮點運算效率(其中步特征線-指數(shù)函數(shù)在混合精度情況下效率較低，約36%)。因此，核函數(shù)在雙精度和混合精度情況下表現(xiàn)為指令密集型，密集的指令操作是其性能瓶頸。

在單精度情況下，由于本文使用的GPU具有強大的單精度浮點運算能力(單/雙精度峰值浮點運算能力之比為32∶1)，更大的單精度指令吞吐量導(dǎo)致了更高的計算效率。但程序的單精度浮點運算效率不足硬件峰值的10%(圖5)，同時圖7所示的硬件單精度浮點運算單元利用率不足40%，因此，密集的指令不再是制約程序性能的瓶頸。分析圖8可知，單精度情況下造成線程束指令簽發(fā)阻斷的主要因素為訪存相關(guān)(約90%)，而圖7所示的訪存效率尚未達到15%，可見制約程序性能的并不是訪存帶寬，而是訪存延遲。

若以浮點運算效率(實際達到的硬件峰值計算能力的百分比)衡量GPU硬件利用率，則如圖5所示，輸運掃描核函數(shù)在雙精度和混合精度情況下具有較高的硬件利用率(36%～67%)，在單精度情況下硬件利用率較低(4%～9%)。

3) GPU對CPU加速效果分析

為對比GPU相對CPU計算的加速效果，表5列出采用不同類型核函數(shù)及不同浮點運算精度情況下，單CPU線程和單GPU的總計算時間對比，同時給出了相同計算條件下單GPU對單CPU線程的加速效果。從加速效果看，雙精度和單精度情況下，單GPU對單CPU線程加速效果分別達到35倍和100倍以上。

表5 不同類型核函數(shù)及浮點運算精度下的計算時間Table 5 Runtime under different kernel functions and floating-point precisions

本文實現(xiàn)了單CPU綁定單GPU，將計算密集的輸運掃描部分移植到GPU進行并行計算，計算流程控制等邏輯相關(guān)工作仍在CPU進行。在計算過程中，由CPU進程綁定GPU，控制GPU進行相關(guān)計算，GPU以獨立的加速單元形式存在，輸運掃描過程完全在GPU獨立執(zhí)行。當(dāng)涉及多個GPU并行時，這種由CPU進程綁定GPU，控制GPU進行相關(guān)計算的基本模式不會改變。且在多GPU并行時，各GPU的計算理應(yīng)相互獨立，所涉及的耦合部分則通過與GPU綁定的CPU進程來完成。各GPU所執(zhí)行的計算任務(wù)仍為局部問題的輸運掃描計算，因此，本文涉及的針對GPU的算法和相關(guān)性能分析在多GPU并行計算中仍適用。

5 結(jié)論與展望

本文實現(xiàn)了基于GPU的二維MOC并行算法，結(jié)果表明：程序具有良好的計算精度，不同的浮點運算精度均能保證二維C5G7基準題在輸運求解方面的精度要求；采用雙精度運算時，步特征線-線性插值計算性能最優(yōu)；采用混合精度和單精度運算時，步特征線-指數(shù)函數(shù)計算性能最優(yōu)。雙精度及混合精度情況下，輸運掃描核函數(shù)性能瓶頸在于密集的指令運算；而在單精度情況下，輸運掃描核函數(shù)性能主要受訪存延遲的影響。從加速效果來看，雙精度和單精度浮點運算情況下，單GPU對單CPU線程加速效果分別達到了35倍和100倍以上。本文涉及的算法及相關(guān)結(jié)論可推廣到多GPU并行計算，未來可側(cè)重基于大規(guī)模CPUs-GPUs異構(gòu)系統(tǒng)的并行算法研究。