基于GPU的三維擴散方程在反應(yīng)堆計算中的應(yīng)用

王飛飛 王海玲 俞海英

基于GPU的三維擴散方程在反應(yīng)堆計算中的應(yīng)用

王飛飛1,2王海玲1,2俞海英3

1（中國科學(xué)院上海應(yīng)用物理研究所 嘉定園區(qū) 上海 201800）2（中國科學(xué)院大學(xué) 北京 100049）3（蘇州工業(yè)職業(yè)技術(shù)學(xué)院 軟件與服務(wù)外包學(xué)院 蘇州 215000）

本文介紹了基于統(tǒng)一計算設(shè)備架構(gòu)(Compute Unified Device Architecture, CUDA)的圖形處理器(Graphic Processing Unit, GPU)計算環(huán)境在釷基熔鹽堆(Thorium Molten Salt Reactor, TMSR)設(shè)計平臺的建立，并將反應(yīng)堆球場計算軟件SRAC (Structure Research and Analysic Corporation)的中子三維擴散計算模塊移植到GPU上進(jìn)行測試及結(jié)果驗證。采用中心點差分方法推導(dǎo)出三維擴散計算的差分方程，并用超松弛迭代法(Successive Over Relaxation Method, SOR)求解φ，研究了SOR迭代算法的并行實現(xiàn)過程。結(jié)果表明，移植的GPU模塊部分計算正確，計算速度得到有效提升，驗證了TMSR設(shè)計平臺在GPU計算環(huán)境下可正常工作。

圖形處理器，統(tǒng)一計算設(shè)備架構(gòu)，釷基熔鹽堆設(shè)計平臺，三維擴散計算，超松弛迭代法

隨著第三代核反應(yīng)堆系統(tǒng)概念設(shè)計的發(fā)展，其堆芯存在結(jié)構(gòu)復(fù)雜、材料繁多、能譜多樣等特點，提高了對反應(yīng)堆物理計算方法的要求。當(dāng)前圖形處理器(Graphic Processing Unit, GPU)計算技術(shù)得到了快速發(fā)展，其在相關(guān)反應(yīng)堆計算如ANSYS[1]、GEANT4[2]、MCNP[3]等軟件中已被廣泛運用，并取得較好結(jié)果。

SRAC (Structure Research and Analysic Corporation)是日本原子能研究結(jié)構(gòu)JAERI開發(fā)的一款Linux下的應(yīng)用于各種類型反應(yīng)堆的中子物理分析程序系統(tǒng)，主要用于堆芯中子擴散計算和輸運計算。目前的多維中子擴散計算，主要依賴有限差分法，它有完整的理論基礎(chǔ)，適用于各類反應(yīng)堆，通過它推導(dǎo)出的差分方程的求解方法是成熟的。對于大型的反應(yīng)堆芯內(nèi)部，由于系數(shù)矩陣階數(shù)相當(dāng)大，利用直接解法比較困難。而迭代法具有程序邏輯簡單、占用存儲空間小等特點，因此通常是利用迭代方法求解差分方程組。當(dāng)前對迭代算法的研究比較成熟，一般用于計算的有雅可比迭代[4]、賽德爾迭代[5]和超松弛法迭代(Successive Over Relaxation Method, SOR)[6]等。當(dāng)差分方程組的階數(shù)比較高時，SOR迭代的收斂速度遠(yuǎn)超過前兩種，因此使用SOR迭代求解三維擴散差分方程會極大縮短計算時間和提高效率。

近年來，由于在GPU中加入了可編程內(nèi)部處理器，使GPU具備很高的運算性能。由于GPU具有較高的浮點計算、內(nèi)存帶寬、多流水線結(jié)構(gòu)及向量處理能力，使其對于計算密集型的科學(xué)領(lǐng)域有非常大的吸引力，越來越成為通用計算的一個有效的并行平臺。

本文討論了在通用超算平臺上的GPU計算環(huán)境的建立，針對釷基熔鹽堆(Thorium Molten Salt Reactor, TMSR)[7]設(shè)計平臺建立了基于統(tǒng)一計算設(shè)備架構(gòu)(Compute Unified Device Architecture, CUDA)的GPU計算環(huán)境，選取了SRAC軟件中的中子多維擴散計算模塊進(jìn)行了GPU計算移植并對結(jié)果進(jìn)行了驗證。

1 GPU計算環(huán)境的建立

1.1 TMSR設(shè)計平臺

TMSR核能系統(tǒng)項目的目標(biāo)是為了研究基于TMSR的先進(jìn)核能系統(tǒng)。熔鹽堆物理設(shè)計與傳統(tǒng)反應(yīng)堆物理在內(nèi)容和方法上都有區(qū)別，具有特殊的結(jié)構(gòu)、動力學(xué)及材料特性。熔鹽堆中燃料組分、燃耗均勻分布、燃料流動和靜態(tài)情況下功率變化、溫度反應(yīng)性系數(shù)等都有很大區(qū)別，堆物理模擬計算與傳統(tǒng)反應(yīng)堆有很大差別。

TMSR設(shè)計平臺作為專項中“反應(yīng)堆物理與工程”項目的課題之一，專門為熔鹽堆的設(shè)計分析提供高效的數(shù)據(jù)存儲服務(wù)和高性能的數(shù)據(jù)處理計算服務(wù)，所有相關(guān)的設(shè)計與計算均在平臺上展開。當(dāng)前在高性能計算領(lǐng)域混合計算成為新的關(guān)鍵詞，位于世界前列的超級計算機有很多都采用了中央處理器(Central Processing Unit, CPU)+GPU的混合計算模式，特別是針對熔鹽堆物理設(shè)計的計算要求，TMSR設(shè)計平臺從開始就設(shè)計為具有CPU+GPU混合結(jié)構(gòu)，具有專門的高密度GPU計算節(jié)點。

1.2 TMSR平臺的硬件環(huán)境

TMSR設(shè)計平臺主要包含計算節(jié)點（刀片GPU計算節(jié)點和專用GPU計算節(jié)點）、登錄節(jié)點、胖節(jié)點、管理服務(wù)器節(jié)點、計算網(wǎng)絡(luò)以及一個并行存儲系統(tǒng)。整個設(shè)計平臺計算集群的節(jié)點類型及數(shù)量如表1所示。

表1 TMSR設(shè)計平臺計算節(jié)點Table 1 Compute nodes of TMSR design platform.

1.3 基于CUDA架構(gòu)的GPU計算環(huán)境

目前主流的GPU計算環(huán)境有OpenCL與NVDIA公司的CUDA開發(fā)環(huán)境，其中CUDA近年來發(fā)展比較成熟，且運行效率更高，而且TMSR設(shè)計平臺使用的GPU顯卡型號是NVDIA公司的M2090，因此我們采用CUDA開發(fā)環(huán)境來構(gòu)建TMSR設(shè)計平臺的GPU計算環(huán)境。

CUDA編程模型對GPU的硬件細(xì)節(jié)進(jìn)行了抽象， CUDA目前已可支持C、C++、Fortran、Python等高級語言，為借助GPU進(jìn)行高性能計算提供了通用開發(fā)環(huán)境，是一種高效的密集計算解決方案。CUDA體系結(jié)構(gòu)包括驅(qū)動、運行環(huán)境和開發(fā)庫(圖1)。驅(qū)動是GPU工作的基礎(chǔ)，不同版本的GPU必須使用相應(yīng)的驅(qū)動程序，為硬件設(shè)備提供抽象訪問接口。運行時環(huán)境提供了應(yīng)用開發(fā)接口和運行時組件，包括基本數(shù)據(jù)類型的定義和各類計算、類型轉(zhuǎn)換、內(nèi)存管理、設(shè)備訪問和執(zhí)行調(diào)度等函數(shù)。開發(fā)庫是基于CUDA技術(shù)所提供的應(yīng)用開發(fā)庫，如CUFFT[8]、CUBLAS[9]等可大幅提高某些應(yīng)用程序的開發(fā)效率。

圖1 CUDA架構(gòu)體系Fig.1 CUDA architecture.

2 SRAC軟件多維計算模塊在GPU中移植

在反應(yīng)堆的物理設(shè)計或臨界計算中，通常需要對多群中子擴散方程[10]進(jìn)行求解。一般情況下，多群擴散方程求解時，解析方法是無法解決的困難問題，而必須借助于數(shù)值方法。本文主要從物理設(shè)計的角度討論利用SOR迭代法求解擴散方程。

2.1 多群中子三維擴散方程理論模型

反應(yīng)堆多群擴散方程可寫成：

式中，φg(r)為第g群的中子注量率；Σt·g為第g群的宏觀總截面；Σg′?g為第g′群到第g群的散射截面；ν為每次裂變產(chǎn)生的中子數(shù);1=′Σg為第g′群的宏觀裂變

截面；Dg為第g群的擴散系數(shù)；Σf為宏觀裂變截面。

對于式(1)的數(shù)值解法，首先把式(1)的微分方程化為便于在計算機上求解的差分方程，然后找求差分方程組的計算方法。本文主要討論三維多群擴散的差分方程推導(dǎo)方法。

2.2 三維差分方程的推導(dǎo)

三維問題的差分解法是把堆芯內(nèi)部分割成許多空間的網(wǎng)點，例如對于xyz幾何形狀情況，則用xi=xi?1+hi、yj=yj?1+lj、zk=zk?1+rk三組直線把堆芯內(nèi)部分成空間網(wǎng)格，可以用和二維類似方法來建立節(jié)點(i, j, k)上的差分方程。三維(x, y, z)幾何情況下擴散式(1)可以寫成：

對于中心點差分情形，討論(i, j, k)節(jié)點，將式(3)在(i, j, k)節(jié)點的網(wǎng)格區(qū)域上積分：

由式(4)得到：

式(5)第一項表示對Vi,j,k的6個側(cè)面求積，式(5)中的導(dǎo)數(shù)應(yīng)用式(6)求近似：

其中：

考慮到D可能不連續(xù)，對于左側(cè)面表面有：

式中，φi?1/2,j,k為交界面上的通量密度；T11為對左側(cè)面求積的值。利用相鄰交界面上中子流連續(xù)條件：

因而：

將式(9)代入式(7)，可得：

其中：

同理，求出T12、T13、T14、T15、T16各項。而

將T11、T12、T13、T14、T15、T16、T2、T3代入式(4)，得到三維差分方程：

其中，各個系數(shù)容易推得為：

由上述討論可以看到，將多維擴散方程化為差分方程組時，最終問題便歸于求解式(15)方程組：式中，A為n階七角矩陣。對于大型的堆芯，n值非常大，因而利用直接解法顯得不是很容易，一般利用迭代方法求解。目前使用最為廣泛的迭代法有雅可比迭代法、賽德爾迭代法和SOR迭代法等。本文采用SOR法，一般如果SOR因子ω選的比較得當(dāng)，方程組的收斂速度會快很多。

2.3 SOR迭代法

在差分方程組中，由于階數(shù)比價高，雅可比和賽德爾迭代方法速度比較緩慢，而SOR通過引用參數(shù)ω，用外推辦法來加速收斂，收斂速度得到極大提高。

SOR迭代格式為：對于三維中子擴散方程，SOR迭代格式為：

將式(16)代入式(18)便可得到：

2.4 并行算法實現(xiàn)方案

設(shè)單塊線程數(shù)為m，塊數(shù)為n，網(wǎng)格數(shù)為l，這樣形成的總線程數(shù)為m×n×l。如圖2所示，在并行計算中，由于求m×n×l個φ的過程中相互獨立，利用GPU的單指令多數(shù)據(jù)的特性，分配m×n×l個處理器，單個處理器計算單個求φ值的迭代過程，處理器之間并行運行。在求誤差時，需要遍歷每個計算結(jié)果，所以不適合做GPU并行執(zhí)行?？紤]到在迭代中，只要誤差值大于ε，所有計算都要重新進(jìn)行，所以在并行程序中需要設(shè)置一個共享變量，當(dāng)誤差值大于ε，則繼續(xù)迭代。迭代結(jié)束后，每個線程需要訪問共享變量，比較誤差值和ε的大小來確定程序繼續(xù)迭代還是輸出結(jié)果。對于SOR迭代法，在計算新值時，使用φi?1,j,k、φi,j?1,k、φi,j,k?1的新值和φi+1,j,k、φi,j+1,k、φi,j,k+1的舊值。由于在計算過程中，新值在不同的處理器中產(chǎn)生，因而考慮將φi,j,k存儲在GPU的全局內(nèi)存中，使每個處理器根據(jù)處理器號順序計算φi,j,k。

圖2 GPU上SOR迭代并行處理流程圖Fig.2 Parallel processing flow chart of SOR on GPU.

在上述算法流程圖步驟中，將求解tempi,j,k值的過程作為內(nèi)核程序，在GPU上進(jìn)行多線程并行計算，使用C語言編寫的CUDA核來實現(xiàn)。通常情況下，內(nèi)核函數(shù)以網(wǎng)格(grid)的形式執(zhí)行，每個網(wǎng)格由若干個線程塊(block)組成，每一個線程塊又由最多512個線程(thread)組成。由于分配了m×n×l的處理器，因此單個處理器就處理一個線程。

本文所選用的運算實例是SRAC軟件中用于模擬中子三維擴散計算的CITITION程序。該程序先任意假設(shè)初始裂變源分布，并且給定初始的k(0)值，把它們作為初始迭代項代入圖2的第一步迭代求初始值算式中。這樣，由它就能夠求出中子通量密度分布φi,j,k，在判斷誤差和闕值ε大小后，將求得的φi,j,k代入求解tempi,j,k值，求得第二次迭代裂變中子分布，并由它可以求得有效增殖因數(shù)的新的估計值和第二代迭代源，依次類推，逐次地迭代下去。

2.5 實驗結(jié)果

計算機模擬實驗平臺選用TMSR設(shè)計平臺上GPU計算節(jié)點Intel(R)6核CPU處理器和NVIDIA M2090 GPU顯卡。串行計算總共的計算時間為15min，加速比為15/3.4=4.4，并行計算的計算時間為3.4 min。對于SRAC軟件中三維擴散計算模塊，通過對SOR迭代法的并行執(zhí)行，計算速度是串行執(zhí)行的4.4倍。由此可以看出，應(yīng)用GPU實現(xiàn)了更高的計算效率。

表2是功率密度部分參考值與測試值的對比與誤差，可以看出，本文采用的計算方法精確，產(chǎn)生的誤差很小。

表2 功率密度部分參考值與測試值對比與誤差Table 2 Power density reference value compared with the test value and error.

3 結(jié)語

通過對通用GPU計算技術(shù)的調(diào)研及研究，建立了基于CUDA架構(gòu)的TMSR設(shè)計平臺的GPU計算環(huán)境。通過中子物理分析軟件SRAC的三維擴散計算模塊的GPU移植，驗證了整個計算集群的GPU計算環(huán)境可正常工作并可看到移植后的GPU計算部分的計算速度有明顯的提升。目前整個TMSR設(shè)計平臺的GPU計算環(huán)境已經(jīng)投入使用。

另外，國內(nèi)外能用于熔鹽堆物理分析和模擬設(shè)計的相關(guān)軟件很少，需要自主開發(fā)適合熔鹽堆計算的物理設(shè)計軟件。本文從設(shè)計階段就引入GPU計算方法，在很大程度上提高了計算速度，從而提高整個物理設(shè)計及模擬的效率。我們將進(jìn)一步在TMSR設(shè)計平臺上開展GPU計算技術(shù)的應(yīng)用研究，加快GPU計算在反應(yīng)堆設(shè)計軟件中的應(yīng)用步伐。

1 王惠, 郭培卿, 陳小龍. ANSYS和Abaqus軟件加速性能典型算例測試與分析[J]. 計算機工程與科學(xué), 2013, 35(11): 105?110 WANG Hui, GUO Peiqing, CHEN Xiaolong. Testing and analysis of typical examples for ANSYS and Abaqus software GPU-accelerated performance[J]. Computer Engineering and Science, 2013, 35(11): 105?110

2 Wright D H, Koi T, Folge G. Low and high energy modeling in Geant4[J]. AIP Conference Proceedings, 2007, 1: 896

3 胡道未, 張雄杰, 湯彬. 基于MCNP蒙卡軟件的并行運算性能研究[J]. 科技信息, 2013, 1(16): 148?150 HU Daowei, ZHANG Xiongjie, TANG Bin. Parallel computing performance based on MCNP Monte Carlo software[J]. Science and Technology Information, 2013, 1(16): 148?150

4 高靜. 系數(shù)矩陣為L-矩陣的線性方程組的新迭代法[J].蘇州市職業(yè)大學(xué)學(xué)報, 2011, 6(22): 52?55 GAO Jing. A new iteration method for linear equations with the coefficient matrix being L-matrices[J]. Journal of Suzhou Vocational University, 2011, 6(22): 52?55

5 王珂, 韋璋劍. 基于高斯-賽德爾迭代法的電纜接地電流計算方法[J]. 廣西電力, 2013, 6: 15?17 WANG Ke, WEI Zhangjian. Cable grounding current calculation method based on Gauss-Seidel iterative method[J]. Guangxi Electric Power, 2013, 6: 15?17

6 王詩然. 稀疏線性方程組求解的逐次超松弛迭代法[J].沈陽師范大學(xué)學(xué)報(自然科學(xué)版), 2006, 10(22): 407?410 WANG Shiran. Sparse linear equations successive over relaxation iteration method[J]. Journal of Shenyang Normal University (Natural Science), 2006, 10(22): 407?410

7 江綿恒, 徐洪杰, 戴志敏. 未來先進(jìn)核裂變能——核能系統(tǒng)[J]. 中國科學(xué)院院刊, 2012, 27(3): 366?374 JIANG Mianheng, XU Hongjie, DAI Zhimin. Advanced fission energy program-TMSR nuclear energy system[J]. Bulletin of Chinese Academy of Sciences, 2012, 27(3): 366?374

8 馮煌. GPU圖像處理的FFT和卷積算法及性能分析[J].計算機工程與應(yīng)用, 2008, 44(2): 120?129 FENG Huang. Implementation and performance of FFT and convolution in image filtering on GPU[J]. Computer Engineering and Applications, 2008, 44(2): 120?129

9 Carcenac M. From tile algorithm to stripe algorithm: a CUBLAS-based parallel implementation on GPUs of Gauss method for the resolution of extremely large dense linear systems stored on an array of solid state devices[J]. The Journal of Supercomputing, 2014, 68(1): 365?413

10 屠柱國, 湯裕仁, 黃艾香, 等. 用有限元法求解三維中子擴散方程的研究“三維中子擴散有限元程序3DFEMJS”[J]. 核科學(xué)與工程, 1983, 3(2): 115?127 TU Zhuguo, TANG Yuren, HUANG Aixiang, et al. On solving neutron diffusions with finite element method-“3DFEMJS”: a three-dimension few group neutron diffusion finite element code[J]. Nuclear Science and Engineering, 1983, 3(2): 115?127

CLC TL338

Research on three dimensional diffusion equation in the reactor calculation application based on GPU

WANG Feifei1,2WANG Hailing1YU Haiying3
1(Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Jiading Campus, Shanghai 201800, China) 2(University of Chinese Academy of Sciences, Beijing 100049, China) 3(Suzhou Institute of Industrial Technology School of Software and Service Outsourcing, Suzhou 215000, China)

Background: With the development of the current nuclear reactor system, some characteristics such as complex structure, material variety and diverse spectrum of reactor cores need more and more computation consumption for reactor physics. Graphic Processing Unit (GPU) computing technology is widely used on the reactor calculation due to its rapid development. Purpose: This paper illustrates the GPU computing environment based on Compute Unified Device Architecture (CUDA) architecture, which is established on the design platform of Thorium Molten Salt Reactor (TMSR). And the neutron three-dimensional diffusion calculation module of the reactor would be transplanted into the GPU for testing and verification. Methods: First of all, center difference method was used to deduce the difference equations of three dimensional diffusion calculations, make it suitable for GPU. Then, the over relaxation iteration method was applied to solve these difference equations. Finally, the over relaxation iteration algorithm was studied in the parallel implementation of GPU processes. Results: The results show that the transplantation of GPU computing module section can be used to obtain correct output, and it can effectively improve the computing speed. Conclusion: The GPU computing environment on the TMSR design platform can work normally by carefully transplanting of several key methods for difference equations based on the CUDA of GPU.

Graphic Processing Unit (GPU), Compute Unified Device Architecture (CUDA), Design platform of Thorium Molten Salt Reactor (TMSR), 3D diffusion calculation, Successive Over Relaxation Method (SOR)

TL338

10.11889/j.0253-3219.2015.hjs.38.040502

中國科學(xué)院未來先進(jìn)核裂變能先導(dǎo)專項(No.XD02001004)資助

王飛飛，男，1986年出生，2009年畢業(yè)于哈爾濱工業(yè)大學(xué)，現(xiàn)為碩士研究生，核技術(shù)與應(yīng)用專業(yè)

2014-10-16，

2015-03-19