基于OpenMP的中子輸運方程特征線法并行計算研究

于 銳，趙 強

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

基于OpenMP的中子輸運方程特征線法并行計算研究

于 銳，趙 強*

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

特征線法是目前求解反應堆中子輸運方程的主要計算方法之一。本文開發了基于OpenMP的中子輸運方程特征線法并行計算程序，以提高特征線法的計算效率。OpenMP是共享存儲體系結構上的一個并行編程模型，采用Fork-Join并行執行方式，適合于SMP共享內存多處理系統和多核處理器體系結構。通過相關基準題測試驗證，表明所開發的程序在有效增殖因數以及相對中子通量(歸一化柵元功率)分布等參數上都能取得良好的精度，且使用OpenMP能取得良好的加速效果，使計算時間顯著減少。

中子輸運方程；特征線法；OpenMP；并行計算

特征線法求解中子輸運方程的基本思想是利用中子飛行的軌跡(即特征線)對求解區域進行掃描跟蹤。理論上適用于任意復雜幾何中子輸運問題的求解，可避免柵元和組件均勻化過程引起的精度下降，因而成為目前反應堆物理計算的一個研究熱點和重點，是反應堆物理計算領域目前被廣泛采用的一種求解輸運方程的方法[1]。特征線法為保證足夠的精度需將網格劃分密集并用大量的特征線來進行跟蹤掃描，導致計算速度緩慢[2]，另外特征線法本身具有收斂慢的特性[3]，限制了特征線法在實際工程中的廣泛應用。國內外開展了大量的加速計算方法研究，一般可歸為兩大類：一類是迭代加速技巧，如粗網再平衡(CMR)、粗網有限差分法(CMFD)、廣義粗網有限差分法(GCMFD)等；另一類是利用計算機技術的發展，如利用GPU、CUDA以及MPI等方式進行加速[4]。

OpenMP(open multi-processing)是適用于共享內存多處理器體系結構的可移植并行編程模型，其應用程序接口由SGI公司發起，由一些主要的計算機硬件與軟件廠商指定并得到認可[5]。OpenMP支持多種系統下的C/C++和Fortran，屬于第二類加速計算方法。使用OpenMP進行并行設計，簡單方便，可使單個計算機的多核能被有效利用，提高程序在單個計算機上的計算效率，并為以后使用“多臺計算機之間并行-單個計算機內部并行”的策略做前期準備工作。本文使用Fortran語言建立基于特征線法的中子輸運方程求解程序和基于OpenMP的并行程序。

1 特征線法求解輸運方程

特征線形式的中子輸運方程離散形式可寫成如下方程[6]：

(1)

在假定已知初始入射角通量和平源近似的條件下，對式(1)積分可求得沿某一方向飛行的中子穿出子區的出射角中子通量：

(2)

沿m方向子區i內的第k段特征線的平均角通量可由下式確定：

(3)

對所有同方向穿過子區i內的中子角通量密度進行體積加權，即可獲得該子區內沿m方向角中子通量密度的平均值，即：

(4)

其中，δAm為特征線密度。

計算得到每一子區沿m方向的平均中子通量后，進行方向加權，便可獲得每一子區的中子通量：

(5)

為保證各子區的面積與實際面積一致，一般按下式進行特征線長度的修正：

(6)

其中，Vi為子區實際面積。

2 基于OpenMP的并行計算

圖1 OpenMP并行框架Fig.1 OpenMP parallel framework

OpenMP是共享存儲體系結構上的一個并行編程模型，適合于多核CPU上的并行程序設計。它以線程為基礎，采用Fork-Join并行執行方式。程序開始于一個單獨的主線程，然后主線程一直串行執行，直到遇到需進行并行計算時，創建新線程或喚醒已有線程來執行并行任務，在執行并行任務時，主線程和派生線程共同工作。在執行完并行域后，派生線程退出或掛起，最后回到單獨的主線程中。OpenMP并行框架如圖1所示[7]。

通過直播、微課等推送形式，實現數字化預習和情況反饋，精準掌握來自當前學生的學情情況和分析數據。通過基于大數據的智能評判系統實現預習預設的問題評測和課前作業自動批改，評估學生已掌握的知識多少和理解程度，自動進行數據分析并及時反饋給教師，教師據此實現針對性的備課，安排適合的教學策略，進行合理的教學設計。

特征線求解輸運方程時需構造源迭代算法，內迭代過程使中子通量收斂后，進行外迭代更新有效增殖因數。本文在內迭代過程中采用OpenMP進行并行化設計，以縮短每次迭代計算時間，從而使整體的計算時間減少，以達到加速效果。在求解每一方向的每一段特征線在子區內的平均中子角通量密度時，通過在原Fortran串行程序原代碼中嵌入OpenMP編譯制導語句進行并行化設計。利用OpenMP編譯制導語句“!MYMOMP PARALLEL DO”創建一個并行域，并指定在此語句之后開始執行并行過程，也即在此語句程序之后的每一段特征線在子區內的平均中子角通量密度求解模塊由多線程執行，由編譯制導語句“!MYMOMP END PARALLEL DO”表明并行域的結束，后續程序按照串行方式執行。在上述兩個編譯制導語句之間具體由幾個線程來執行并行域部分，由用戶自行設置，利用環境變量“OMP_NUM_THREADS”設定更改執行線程的數量，設定值若為4，則在并行域中會有4個線程執行任務。基于OpenMP的并行計算程序流程圖示于圖2。

3 數值計算結果及并行性能分析

本文使用Fortran語言編制特征線法求解中子輸運方程串行程序和基于OpenMP的并行計算程序。計算時所使用計算機的相關參數如下：處理器為Intel(R) Core(TM) i5-2400 CPU@ 3.10 GHz，Windows 32位操作系統，安裝內存為4 GB，編譯器為Visual Studio 2010。本文利用典型壓水堆柵元和沸水堆組件基準題驗證了兩個版本程序的正確性，最后進行了基于OpenMP的并行性能分析。

3.1 壓水堆柵元基準題計算

為驗證基于特征線法的中子輸運方程求解程序的正確性，本文對典型壓水堆柵元進行了計算，柵元由燃料區、包殼區以及慢化劑區構成，其尺寸如圖3所示，截面列于表1[8-9]。

進行柵元特征線求解時，將柵元劃分為40個子區，其中燃料區劃分為16個子區，包殼區劃分為8個子區，慢化劑區劃分為16個子區，子區劃分形式示于圖4。在進行計算時每一方向選擇的特征線為43條，每一象限內選擇14個方位角和3個極角進行兩群計算，采用鏡面反射邊界條件。將計算出的增殖因數與參考值進行對比，結果列于表2。

圖2 基于OpenMP的并行計算程序流程圖Fig.2 Flow chart of parallel program based on OpenMP

圖3 PWR柵元基準題尺寸Fig.3 Dimension of PWR cell benchmark

由表2可看出，本文所開發程序的計算結果良好，與參考值的偏差小，約為23 pcm，可達到較好的精度。并行與串行的程序計算結果一致，從而驗證了并行方法的正確性。

表1 PWR柵元基準題宏觀截面Table 1 Macroscopic cross section for PWR cell benchmark

圖4 典型PWR柵元子區劃分Fig.4 Mesh of typical PWR cell

表2 典型柵元基準題增殖因數計算結果Table 2 Computed kinf results of typical PWR cell

3.2 沸水堆組件基準題計算

該基準題是一包含2個釓柵元的4×4沸水堆組件，其幾何尺寸及柵元布置方式示于圖5。其中：編號1～5代表標準燃料柵元，由富集度為3%的UO2組成；編號6為含有毒物的燃料，由富集度為3%的UO2和富集度為3%的Gd2O3組成。包殼材料為鋯-2合金，以水作為慢化劑[10]。

在計算該基準題時，將組件內的每一個柵元劃分為40個子區，具體劃分方式與進行典型壓水堆柵元計算時對柵元子區的劃分一致，每一象限內選擇14個方位角和3個極角，每個方向有43條特征線，進行兩群計算，所采用的邊界條件為反射邊界條件，計算時所用的截面列于表3[10]。

圖5 BWR組件基準題尺寸及柵元布置Fig.5 Dimension of BWR assembly benchmark and layout of cell

在計算該基準題時，有效增殖因數和中子通量的收斂準則均設置為10-6，為便于與DRAGON參考值以及MOCUM程序給出的結果進行對比，利用式(6)對裂變率進行歸一化處理。其中有效增殖因數的計算結果列于表4，歸一化柵元功率分布示于圖6。

(7)

由表4可知，本文所開發的兩個版本的特征線求解程序計算結果一致，驗證了并行方法的正確性。在求解有效增殖因數時，與DRAGON程序的計算結果的偏差為72 pcm，與MOCUM程序的計算值的偏差為50 pcm，可見本文所建立的程序可取得良好的計算精度。

表3 BWR組件基準題宏觀截面Table 3 Macroscopic cross section for BWR assembly benchmark

表4 BWR組件基準題有效增殖因數計算結果Table 4 Computed keff results of BWR assembly benchmark

圖6 BWR組件基準題歸一化柵元功率分布計算值與參考值結果對比Fig.6 Calculated result normalized cell power of BWR assembly benchmark vs reference result

組件計算時并行程序與串行程序計算所得的功率一致，圖6為BWR組件基準題歸一化柵元功率分布本文計算值與DRAGON參考值的對比，可見本文所建立的求解程序可滿足精度要求。其中，最小相對偏差為0.002 8%，最大相對偏差為0.204 4%，最大相對偏差出現在含有毒物的燃料柵元，因為此區域中子通量密度梯度變化大，本文依然采用了與標準UO2柵元相同的16子區的燃料子區劃分方式，為減小偏差可考慮將燃料子區剖分得更密。

3.3 OpenMP的并行性能分析

3.2節沸水堆組件基準題計算表明，利用OpenMP進行并行化設計后，無論是有效增殖因數還是相對中子通量(歸一化柵元功率)分布都與串行時保持一致，從而驗證了并行方法的正確性。

為進行OpenMP并行性能分析，本文以3.2節沸水堆組件基準題為計算算例。在求解每一段特征線在子區內的平均中子角通量密度時，通過在原Fortran串行程序源代碼中嵌入OpenMP編譯制導語句進行并行化設計。利用OpenMP編譯制導語句以及環境變量設置，創建并行執行的并行域，將線程數設置為4。

為保證數據的可靠性，采集數據時每組數測試10次，剔除異常值后取平均值。如表5所列，使用OpenMP對程序進行并行化設計后，計算時間下降。當線程數為4時，總計算時間從原來的503 s降至217 s，其中可并行部分從原串行程序中的384 s降至101 s，加速比達到3.80，大幅節省了計算時間。

表5 基于OpenMP的并行性能Table 5 Parallel performance based on OpenMP

注：1) 并行程序中可并行部分計算內容在串行程序中的運算時間

4 結論

本文建立了基于OpenMP的中子輸運方程特征線并行求解程序，利用典型壓水堆柵元基準題和沸水堆組件基準題進行了程序的驗證。結果表明，無論是有效增殖因數還是相對中子通量(歸一化柵元功率)分布都能取得良好的計算精度，基于OpenMP的并行程序獲得了良好的加速效果，在相同計算精度下，大幅提高了計算速度，節省了計算時間。

[1] 湯春桃. 中子輸運方程特征線解法及嵌入式組件均勻化方法的研究[D]. 上海：上海交通大學，2009.

[2] 湯春桃，張少泓. CMFD 加速在特征線法輸運計算中的應用[J]. 核動力工程，2009，30(5)：8-12.

TANG Chuntao, ZHANG Shaohong. Application of coarse-mesh finite difference acceleration in transportation calculation by method of characteristic[J]. Nuclear Power Engineering, 2009, 30(5): 8-12(in Chinese).

[3] 張知竹，李慶，王侃. 三維特征線的并行方法研究[J]. 原子能科學技術，2013，47(增刊)：38-42.

ZHANG Zhizhu, LI Qing, WANG Kan. Parallelization method for three dimensional MOC calculation[J]. Atomic Energy Science and Technology, 2013, 47(Suppl.): 38-42(in Chinese).

[4] ZHANG Zhizhu, LI Qing, WANG Kan. Accelerating a three-dimensional MOC calculation using GPU with CUDA and two-level GCMFD method[J]. Annals of Nuclear Energy, 2013, 62: 445-451.

[5] 李建江，薛巍，張武生，等. 并行計算機及編程基礎[M]. 北京：清華大學出版社，2011.

[6] KNOTT D, FORSSEN B H, EDENIUS M. CASMO-4: A fuel assembly burnup program methodology[R]. Studsvik: Studsvik Scandpower, Inc., 1995.

[7] 陳國良. 并行計算：結構·算法·編程[M]. 北京：高等教育出版社，2011.

[8] 張知竹，李慶，王侃. GPU加速三維特征線方法的研究[J]. 核動力工程，2013，34(S1)：18-23.

ZHANG Zhizhu, LI Qing, WANG Kan. Study on acceleration of three-dimensional method of characteristics by GPU[J]. Nuclear Power Engineering, 2013, 34(S1): 18-23(in Chinese).

[9] POSTMA T, VUJIC J. The method of characteristics in general geometry with anisotropic scattering[C]∥Proceedings of International Conference on Mathematics and Computation, Reactor Physics and Environmental Analysis in Nuclear Application. Madrid, Spain: [s. n.], 1999.

[10]XUE Yang, STVAT N. MOCUM: A two-dimensional method of characteristics code based on constructive solid geometry and unstructured meshing for general geometries[J]. Annals of Nuclear Energy, 2012, 46: 20-28.

Research on Parallel Computing of Method of Characteristics of Neutron Transport Equation Based on OpenMP

YU Rui, ZHAO Qiang*

(FundamentalScienceonNuclearSafetyandSimulationTechnologyLaboratory,HarbinEngineeringUniversity,Harbin150001,China)

The method of characteristics (MOC) is one of the main methods for solving reactor neutron transport equation currently. A transport theory code based on the method of characteristics and an OpenMP parallel version of the method of characteristics calculation code were developed. OpenMP is a parallel programming model with shared memory architectures, using Fork-Join parallel execution mode, which is suitable for SMP shared memory multi-processor systems and multi-core processors architecture. The code was verified and validated by different benchmarks. The numerical results demonstrate that the code can give excellent accuracy for both the neutron effective multiplication factor and relative neutron flux (normalized cell power) distribution for neutron transport problem. The use of OpenMP can obtain good acceleration effect, making the calculation time significantly reduced.

neutron transport equation; method of characteristics; OpenMP; parallel computing

2014-07-04;

2014-09-26

于 銳(1989—)，男，寧夏涇源人，碩士研究生，核能科學與工程專業

*通信作者：趙 強，E-mail: zhaoqiang@hrbeu.edu.cn

TL329.2

A

1000-6931(2015)10-1833-06

10.7538/yzk.2015.49.10.1833