基于OpenMP的中子輸運(yùn)方程特征線法并行計(jì)算研究

于 銳，趙 強(qiáng)

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

基于OpenMP的中子輸運(yùn)方程特征線法并行計(jì)算研究

于 銳，趙 強(qiáng)*

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

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

中子輸運(yùn)方程；特征線法；OpenMP；并行計(jì)算

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

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

1 特征線法求解輸運(yùn)方程

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

(1)

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

(2)

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

(3)

對所有同方向穿過子區(qū)i內(nèi)的中子角通量密度進(jìn)行體積加權(quán)，即可獲得該子區(qū)內(nèi)沿m方向角中子通量密度的平均值，即：

(4)

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

計(jì)算得到每一子區(qū)沿m方向的平均中子通量后，進(jìn)行方向加權(quán)，便可獲得每一子區(qū)的中子通量：

(5)

為保證各子區(qū)的面積與實(shí)際面積一致，一般按下式進(jìn)行特征線長度的修正：

(6)

其中，Vi為子區(qū)實(shí)際面積。

2 基于OpenMP的并行計(jì)算

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

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

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

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

3 數(shù)值計(jì)算結(jié)果及并行性能分析

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

3.1 壓水堆柵元基準(zhǔn)題計(jì)算

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

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

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

圖3 PWR柵元基準(zhǔn)題尺寸Fig.3 Dimension of PWR cell benchmark

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

表1 PWR柵元基準(zhǔn)題宏觀截面Table 1 Macroscopic cross section for PWR cell benchmark

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

表2 典型柵元基準(zhǔn)題增殖因數(shù)計(jì)算結(jié)果Table 2 Computed kinf results of typical PWR cell

3.2 沸水堆組件基準(zhǔn)題計(jì)算

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

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

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

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

(7)

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

表3 BWR組件基準(zhǔn)題宏觀截面Table 3 Macroscopic cross section for BWR assembly benchmark

表4 BWR組件基準(zhǔn)題有效增殖因數(shù)計(jì)算結(jié)果Table 4 Computed keff results of BWR assembly benchmark

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

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

3.3 OpenMP的并行性能分析

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

為進(jìn)行OpenMP并行性能分析，本文以3.2節(jié)沸水堆組件基準(zhǔn)題為計(jì)算算例。在求解每一段特征線在子區(qū)內(nèi)的平均中子角通量密度時，通過在原Fortran串行程序源代碼中嵌入OpenMP編譯制導(dǎo)語句進(jìn)行并行化設(shè)計(jì)。利用OpenMP編譯制導(dǎo)語句以及環(huán)境變量設(shè)置，創(chuàng)建并行執(zhí)行的并行域，將線程數(shù)設(shè)置為4。

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

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

注：1) 并行程序中可并行部分計(jì)算內(nèi)容在串行程序中的運(yùn)算時間

4 結(jié)論

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

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

[2] 湯春桃，張少泓. CMFD 加速在特征線法輸運(yùn)計(jì)算中的應(yīng)用[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]. 原子能科學(xué)技術(shù)，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] 李建江，薛巍，張武生，等. 并行計(jì)算機(jī)及編程基礎(chǔ)[M]. 北京：清華大學(xué)出版社，2011.

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

[7] 陳國良. 并行計(jì)算：結(jié)構(gòu)·算法·編程[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—)，男，寧夏涇源人，碩士研究生，核能科學(xué)與工程專業(yè)

*通信作者：趙 強(qiáng)，E-mail: zhaoqiang@hrbeu.edu.cn

TL329.2

A

1000-6931(2015)10-1833-06

10.7538/yzk.2015.49.10.1833