CN106126927A - 一种用于反应堆屏蔽计算的并行扫描方法 - Google Patents

Abstract

本发明公开了一种用于核反应堆屏蔽计算的并行扫描方法，该方法是将计算几何用细网格剖分逼近，然后将细网格分成若干个区域，每个区域由一个CPU负责，同时又将一个CPU负责的区域分成若干个并行单元，当一个CPU上的并行单元计算结束后，可以将计算信息传递给相邻的并行单元，从而实现高效的并行屏蔽计算，给核电厂设计的安全性和经济性提供有力的保障。

Description

一种用于反应堆屏蔽计算的并行扫描方法

技术领域

本发明涉及核反应堆屏蔽计算领域，是一种获取反应堆堆芯及其外部三维中子、光子通量密度精细分布的并行计算方法。

背景技术

反应堆中子、光子通量密度分布对整个核电厂设计的安全性和经济性至关重要。中子、光子的辐照会使反应堆厂房内各种设备的材料特性恶化，从而限制整个核电厂的寿命。反应堆屏蔽设计就是为了尽量避免中子、光子从堆芯泄漏到反应厂房中去。因此屏蔽计算也成为了核电厂设计中至关重要的部分，其主要任务就是获取反应堆内的中光子通量分布。

目前广泛应用的方法大体上可以分为：确定论方法和蒙特卡罗方法。其中确定论方法中主要包含了：综合法、点核积分法和三维离散纵标方法，蒙特卡罗方法主要是以MCNP程序为代表的一系列非确定论算法。

由于过去计算条件的限制，综合法和点核积分算法在工程上得到了广泛的应用。所谓“综合法”即通过计算一维和二维的通量分布，从而重构出一个三维通量分布。“点核积分算法”则是将源项分割成许多小源块，计算时每个小源块被当成点源处理，用累积因子的方法计算各点源对剂量点的贡献。由于计算条件的不断改善，“综合法”和“点核积分算法”这类近似较大的算法将逐步退出屏蔽设计的舞台。

毫无疑问，要想得到精确的三维通量分布，就必须进行三维输运计算。目前得到广泛重视的主要算法有三维离散纵标方法和蒙特卡罗方法。由于屏蔽计算通常是深穿透问题，蒙特卡罗方法的收敛性遇到很大的挑战，距离工程应该尚有很长的路要走。三维离散纵标方法的计算量巨大，如果采用串行计算，就无法充分地调用计算资源，难以完成反应堆屏蔽设计的任务。因此，并行三维离散纵标方法是目前离工程应用最近算法。本发明就是一种能够充分调用计算机资源先进并行扫描计算方法。

发明内容

为了克服上述现有技术存在的问题，本发明的目的在于提供一种用于反应堆屏蔽计算的并行扫描方法，充分调用计算机资源以获得堆芯外部精确的三维中子、光子通量分布；该方法是将计算几何用细网格剖分逼近，然后将细网格分成若干个区域，每个区域由一个CPU负责，同时又将一个CPU负责的区域分成若干个并行单元，当一个CPU上的并行单元计算结束后，可以将计算信息传递给相邻的并行单元，从而实现高效的并行屏蔽计算，给核电厂设计的安全性和经济性提供有力的保障。

为了实现上述目的，本发明采取了以下技术方案予以实施：

一种用于核反应堆屏蔽计算的并行扫描方法，该方法包括以下步骤：

步骤一：将反应堆模型沿X方向，Y方向以及Z方向划分成细小网格，用这些细小网格去逼近真实几何，并基于这些细小网格建立三维离散纵标方法的三维细网模型；

步骤二：将三维细网模型沿X方向和Y方向分解成若干二维个子区域，二维子区域的个数就是后续计算中所使用的CPU个数，二维子区域尽量均匀划分，以确保每个CPU的计算量相仿，充分调用计算资源，同时将每个子区域沿Z方向分解成若干并行单元，在每个并行单元内部执行串行计算，并行单元之间执行并行计算；

步骤三：根据离散角度的方向确定首个计算的并行单元，对于三维几何，离散角度的方向落在八个象限中：(-，-，-)，(-，-，+)，(-，+，-)，(-，+，+)，(+，-，-)，(+，-，+)，(+，+，-)，(+，+，+)，而空间扫描顺序和离散角度的方向所在的卦限是一一对应的，当离散角度的方向所在的象限确定时，首个计算的并行单元即为该离散方向对应角点处所在的并行单元，当首个并行单元计算结束以后，将其表面上的出射角通量传给相邻的并行单元作为入射角通量；相邻的并行单元接收到入射角通量后，对应的CPU立即执行计算任务，以此类推，所有CPU将很快地被启动起来，对于每一个并行单元采用轴向贯穿流水线式角度顺序，即将八个卦限分成四组：①(-，-，-)和(-，-，+)，②(-，+，-)和(-，+，+)，③(+，-，-)和(+，-，+)，④(+，+，-)和(+，+，+)，先将第一组内两个卦限的所有角度扫描结束再扫描第二组，接着扫描第三组和第四组；

步骤四，当本次迭代过程结束时，每个CPU计算自己对应并行单元每一个网格的迭代误差，并将两两比较每个网格迭代误差，误差绝对值大的保留小的舍去，这样就可以得到每个并行单元的最大迭代误差，记录每一个并行单元的最大迭代误差，并将所有并行单元的最大迭代误差进行比较得到整个三维细网模型的最大迭代误差，当三维细网模型的最大迭代误差大于规定的收敛限时，需要继续迭代计算，当计算模型的最大迭代误差小于规定的收敛限时，迭代结束，并输出反应堆三维中子、光子通量分布，即完成屏蔽计算。

与现有技术相比，本发明有如下突出优点：

1.保持了原有串行计算的高斯赛德尔迭代格式，并行计算与原有串行计算的迭代次数和计算结果完全相同。

2.避免了传统区域分解并行算法由于分区过多而导致的发散。

3.采用了轴向贯穿流水线式角度扫描顺序，在千核量级内保持了较高的并行效率，传统的并行方法在调用100个CPU时并行效率约为60％，而本发明可以达到90％的并行效率。

附图说明

图1是用细网格剖分逼近反应堆模型。

图2是差分细网格微元。

图3是给定角度下的细网格扫描顺序。

图4是子区域和并行单元的划分。

图5是某一角度下直角几何并行单元扫描顺序。

图6是三维圆柱几何下的并行扫描顺序。

图7是P0子区域和PN子区域间的等待时间。

图8是象限扫描顺序。

图9是并行效率随CPU个数的变化图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细说明。

本发明通过细网差分方法将输运方程离散化，同时将细网区域分成若干个子区域，而在每个子区域内部沿轴向又划分成若干并行单元，通过并行单元之间的通信，将所有CPU耦合在一起，且让并行单元的扫描顺序和离散纵标方法中角通量的扫描顺序保持一致，从而保证并行计算的收敛性和高效性，该方法主要包括以下几个方面：

1)通过引入角度空间离散，将一个连续几何上的输运方程离散成每一个细网格上的平衡方程，建立三维细网模型；

2)利用区域分解的思想，将三维的离散网格沿X方向和Y方向划分为若干个子区域，每个子区域贯穿了Z方向的所有网格。同时，为了提高并行效率以及防止高斯赛德尔迭代格式退化成雅格比迭代格式，将每个子区域沿轴向又划分成若干个并行单元，并且让并行单元间的扫描顺序和当前离散角度对应的扫描顺序一致；

3)当首个并行单元的计算完成后，将边界信息传递给相邻的并行单元，相邻单元接收到信息后开始启动计算，同时采用轴向贯穿流水线式角度扫描顺序，进一步提高并行效率；

4)计算每一个并行单元的迭代误差，只要有一个并行单元内的迭代误差未达到收敛准则，就必须重复迭代，直至收敛，并输出三维中、光子通量分布。

如此，对于每一个能群,只有在第一次扫描计算时，各CPU存在等待情况。在第一次扫描结束后，所有CPU都开始计算自己对应的子区域，直至所有能群计算完毕。

1输运方程离散化

三维中子输运方程通过角度空间离散化，转化为某个角度某个网格上的平衡方程，具体如下：

稳态多群中子输运方程如公式所示：

Ω·▽φ_g(r,Ω)+Σ_t,g(r)φ_g(r,Ω)＝Q_g(r,Ω),g＝1,...,G 公式(1)

式中：

r——中子或光子的位置向量；

Ω——中子或光子的飞行方向向量；

φ_g(r,Ω)——第g能群中子或光子角通量密度；

g——能群标识；

G——能群总数；

Σ_t,g(r)——第g能群位置r处的总截面；

Q_g(r,Ω)——中子输运方程源项。

在屏蔽计算中，第g能群源项的具体形式为：

$\begin{array}{c}{Q}_g(r,\mathrm{\Ω})=Q_g^e(r,\mathrm{\Ω})+Q_g^s(r,\mathrm{\Ω})\\ =Q_g^e(r,\mathrm{\Ω})+\underset{g^{\′}=1}{\overset{G}{\Σ}}{\∫}_{{\mathrm{\Ω}}^{\′}}{\Σ}_{s,g^{\′}-g}(r,{\mathrm{\Ω}}^{\′}\→\mathrm{\Ω}){\mathrm{\φ}}_{g^{\′}}(r,{\mathrm{\Ω}}^{\′}){\mathrm{d\Ω}}^{\′}\end{array}$

式中：

——第g能群的给定外部源项；

——第g能群的散射源项。

经过角度离散后，在第g群、第m离散方向上的形式如公式所示：

1.1空间离散化

将计算几何用细网格逼近如图1所示，再将公式在如图2所示的一个微小的网格上积分，令V_ijk为x_i-1/2＜x＜x_i+1/2，y_j-1/2＜y＜y_j+1/2，z_k-1/2＜z＜z_k+1/2所定义的区域体积，对于左边的第一项有：

对于左边第二项有：

对于左边第三项有：

于是，式经过积分后，可以写成：

为了进一步简化式，定义如下关系式：

V_ijk＝Δx_iΔy_jΔz_k 公式(12)

式中：

φ_{i±1/2,j,k,g,m}——第g能群，第m个离散角度，±x网格边界面上的平均角通量；

φ_{i,j±1/2,k,g,m}——第g能群，第m个离散角度，±y网格边界面上的平均角通量；

φ_{i,j,k±1/2,g,m}——第g能群，第m个离散角度，±z网格边界面上的平均角通量；

φ_i,j,k,g,m——第g能群，第m个离散角度，网格(i,j,k)的平均角通量；

V_ijk——网格(i,j,k)的体积；

Q_i,j,k,g,m——第g能群，第m个离散角度，网格(i,j,k)的平均源强；于是，式可以简化为：

式具有非常明显物理意义，即：网格(i,j,k)上出射、入射粒子关于网格内部消失项及源项的平衡方程，如此，便建立反应堆屏蔽计算的三维细网模型。

1.2空间扫描顺序

为了保证迭代计算的收敛性，空间扫描顺序的选取务必要满足初始误差能够在扫描过程中不断减小，即扫描顺序对误差具有缩小的作用。本发明中采用了“风向”空间扫描顺序，即网格的扫描方向与离散角度的方向保持一致，如图3所示：

假设当前扫描方向为Ω_m，则扫描起始网格为左上角网格(i1，j1)，其入射中子角通量由边界条件给出。(i1，j1)网格计算完成后，得到i方向的出射中子角通量(图中用虚箭头表示)，(i2，j1)网格开始计算(j方向的入射中子角通量由边界条件给出)。同理，在(i2，j1)网格计算完成后，(i3，j1)网格开始计算；(i3，j1)网格计算完成后，(i4，j1)网格开始计算。j方向的扫描过程与i方向一致：(i1，j1)网格计算完成后，给出j方向的出射中子角通量，之后(i1，j2)开始计算(i方向的入射中子角通量由边界条件给出)；基于同样的过程，(i1，j3)、(i1，j4)依次完成各自网格的差分计算。

2区域分解及并行单元的划分

所谓区域分解，就是将细网格划分好的三维几何划分成若干个二维子区域(每个子区域内部依然是三维)如图4所示，其中的“body”即为计算几何，“body”上的网格即为实际计算时所采用网格；“processor”即为一个子区域，该区域内的所有细网格由一个CPU负责计算。图4中将计算几何分解成了9个子区域，由于每个子区域都贯穿了轴向(Z)方向的所有网格，因此也称之为二维子区域。为了使各CPU负载平衡，减少CPU间的等待时间，子区域尽量均匀划分。

所谓并行单元，即图4中的“block”，将每个子区域均匀地划分成若干并行单元。实际计算中从哪个并行单元开始扫描与离散角度所在的象限相关，这将在步骤3中详细说明。需要注意的是在每个并行单元内部执行的串行计算，其扫描顺序即为步骤1.2中的空间扫描顺序。

3并行扫描方法

并行扫描方法是本发明的核心。本发明中的扫描策略是以图4中的“block”(并行单元)为基本单元的。

3.1并行单元扫描顺序

如图5所示，在某一给定的离散角度下左上角的第一个并行单元启动计算，即图5中的step 1，此时只有一个CPU在执行计算任务，其他CPU均处于等待状态。当这个并行单元内所有细网格计算完毕之后，其出射面(黑面，共三个，其余两个不可见)上的角通量传递给相邻三个并行单元，作为其入射角通量。当相邻的三个并行单元接收到入射角度的信息时，立即启动各自的计算，即图5中的step 2，此时共有3个CPU在执行计算任务，其他CPU均处于等待状态。当这三个并行单元内所有网格计算完毕时又各自出射面上的角通量传递给相邻并行单元。因此在step 3中就有6个并行单元同时计算，即6个CPU在执行计算任务。以此类推，所有CPU都将会很快的启动起来。图5展示了三维直角几何的并行扫描策略。对于三维圆柱几何，这种并行策略依然适用，三维圆柱的并行单元扫描顺序如图6所示。值得一提的是，只有在第一个角度执行扫描时，CPU存在等待现象，因为当第一个启动的子区域轴向所有并行单元对于某个离散角度扫描计算结束时，同象限内的第二个离散角度立刻执行计算，无需等待。这也是本发明中的一大亮点。

3.2并行效率分析

并行效率的计算表示式如下：

式中：

η——理论并行效率；

T_s——串行计算时间；

T_p——并行计算时间；

P——所用CPU的个数。

并且，定义如下变量：

M——一个象限中离散角度的个数；

I——X方向的细网个数；

J——Y方向的细网个数；

K——Z方向的细网个数；

I_c——每个并行单元(block)在X方向上包含的细网个数；

J_c——每个并行单元(block)在Y方向上包含的细网个数；

K_c——每个并行单元(block)在Z方向上包含的细网个数；

t₀——计算一个细网格所需要的时间；

T_P1——每个CPU的等待时间；

T_P2——每个CPU的计算时间；

因此，对于一个象限而言，串行计算时间为：

T_s＝MIJKt₀ 公式(16)

所用的CPU个数等于所划分的子区域的个数，即：

如图7所示，从P₀号CPU开始启动计算一直到P_N号CPU开始计算，这期间的等待时间可以从X方向和Y方向上需要等待的次数来计算。同时需要注意，由于一个子区域贯穿了整个Z轴，虽然划分成若干并行单元，但在Z方向上并无等待，因为其隶属于同一个CPU。从图7中容易分析出，在X方向和Y方向各自需要等待的次数分别为：和次。每一次等待的时间恰好等于一个并行单元对于一个角度串行计算所需要的时间，即：I_cJ_cK_ct₀。因此，对于一个CPU其等待时间为：

而对于一个CPU，其计算时间为：

T_p2＝MI_cJ_cKt₀ 公式(19)

因此，本发明的并行效率为：

从并行效率的表达式中，可以看出：当一个象限中离散角度的个数M越大时，其并行效率越高。对于三维问题，离散角度分布在如图8所示的8个象限中(±,±,±)。由于本发明中区域分解是一个三维几何向一个二维子区域的映射，因此(+,+,±)两个卦限对应相同的CPU启动顺序，(+,-,±)两个卦限对应相同的CPU启动顺序，(-,+,±)两个卦限对应相同的CPU启动顺序，(-,-,±)两个卦限对应相同的CPU启动顺序。因此，在本发明中采用图8所示的象限扫描顺序，即：1、5、2、6、3、7、4、8。如此，将Z方向上相同位置处的象限一起扫描，就相当于把8个象限压缩成4个象限，因此每个象限中的离散角度个数翻倍。于是，并行效率公式就变成下列形式：

在本发明中，把这种象限扫描顺序称之为：轴向贯穿流水线式角度扫描顺序。

测试了不同CPU个数并行计算时的并行效率，如图9所示，可以看出，本发明的计算效率很高，当CPU的个数为121时，并行效率依然高达90％，具有很大计算优势。

4收敛性验证

在输运计算中，要求每一个细网的迭代误差都小于规定收敛准则。在本发明中，由于整体几何被分成了若干个子区域，而每个子区域又分为了若干个并行单元。因此，每个CPU在计算完各自负责的一个并行单元时，需要把这个并行单元的最大细网通量误差记录下来。当一个子区域计算完毕时，每个CPU从自己对应的并行单元中找出其中的最大误差，并将该误差发送给其它CPU进行比较，最后得到所有CPU中的最大误差，即为细网通量最大迭代误差。当细网通量最大迭代误差小于规定的收敛准则时，则迭代收敛，否则继续下一次迭代。

Claims (1)

1.一种用于核反应堆屏蔽计算的并行扫描方法，其特征在于：该方法包括以下步骤：

步骤一：将反应堆模型沿X方向，Y方向以及Z方向划分成细小网格，用这些细小网格去逼近真实几何，并基于这些细小网格建立三维离散纵标方法的三维细网模型；

步骤二：将三维细网模型沿X方向和Y方向分解成若干二维个子区域，二维子区域的个数就是后续计算中所使用的CPU个数，二维子区域尽量均匀划分，以确保每个CPU的计算量相仿，充分调用计算资源，同时将每个子区域沿Z方向分解成若干并行单元，在每个并行单元内部执行串行计算，并行单元之间执行并行计算；

步骤三：根据离散角度的方向确定首个计算的并行单元，对于三维几何，离散角度的方向落在八个象限中：(-，-，-)，(-，-，+)，(-，+，-)，(-，+，+)，(+，-，-)，(+，-，+)，(+，+，-)，(+，+，+)，而空间扫描顺序和离散角度的方向所在的卦限是一一对应的，当离散角度的方向所在的象限确定时，首个计算的并行单元即为该离散方向对应角点处所在的并行单元，当首个并行单元计算结束以后，将其表面上的出射角通量传给相邻的并行单元作为入射角通量；相邻的并行单元接收到入射角通量后，对应的CPU立即执行计算任务，以此类推，所有CPU将很快地被启动起来，对于每一个并行单元采用轴向贯穿流水线式角度顺序，即将八个卦限分成四组：①(-，-，-)和(-，-，+)，②(-，+，-)和(-，+，+)，③(+，-，-)和(+，-，+)，④(+，+，-)和(+，+，+)，先将第一组内两个卦限的所有角度扫描结束再扫描第二组，接着扫描第三组和第四组；

步骤四，当本次迭代过程结束时，每个CPU计算自己对应并行单元每一个网格的迭代误差，并将两两比较每个网格迭代误差，误差绝对值大的保留小的舍去，这样就得到每个并行单元的最大迭代误差，记录每一个并行单元的最大迭代误差，并将所有并行单元的最大迭代误差进行比较得到整个三维细网模型的最大迭代误差，当三维细网模型的最大迭代误差大于规定的收敛限时，需要继续迭代计算，当计算模型的最大迭代误差小于规定的收敛限时，迭代结束，并输出反应堆三维中子、光子通量分布，即完成屏蔽计算。