基于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