CN106951683B - 一种用于核电厂厂房屏蔽计算的高效并行扫描方法 - Google Patents

Abstract

本发明公开了一种用于核电厂厂房屏蔽计算的高效并行扫描方法，通过细网格将连续的厂房屏蔽计算模型离散化，同时将细网区域分配到若干个CPU上，而每个CPU上的网格又沿轴向划分成若干并行单元，通过并行单元之间的通信，将所有CPU耦合在一起，将二维网格子区域进行“米”字形划分，划分后的每个CPU子区域都有自己既定的卦限扫描顺序，将所有卦限分成四组，且让四组卦限中的角度同时从二维网格子区域的四个角点同时向中心扫描，结合内部边界条件充分利用阻塞和非阻塞通信，从而保证并行计算的收敛性和高效性；本发明适用于大规模集群并行计算，以便于快速准确获得核电厂厂房的三维中、光子通量分布。

Description

一种用于核电厂厂房屏蔽计算的高效并行扫描方法

技术领域

本发明属于核电厂厂房屏蔽计算技术领域，具体涉及一种计算核电厂厂房三维中子、光子通量密度精细分布的高效并行计算方法。

背景技术

核电厂厂房的中子、光子通量密度分布对整个核电厂的安全性和经济性至关重要。中子、光子的辐照不但会使反应堆厂房内各种设备的材料特性恶化，而且还影响核电厂工作人员的健康。核电厂厂房屏蔽计算就是为了准确获取厂房的中、光子通量分布。因此屏蔽计算也成为了核电厂设计中至关重要的部分。

由于核电厂厂房体积大，且中、光子通量分布从堆芯到外部厂房变化梯度陡，准确计算厂房的中、光子通量分布难度极大。目前广泛应用的方法大体上可以分为：确定论算法和蒙特卡罗算法。其中确定论算法中主要包含了：综合法、点核积分法和离散纵标方法，蒙特卡罗算法主要是以MCNP程序为代表的一系列非确定论算法。

“综合法”通过计算一维和二维的通量分布，从而重构出一个三维通量分布。“点核积分算法”则是将体源分割成许多小源块，并将每个小源块近似成点源，最后积分各点源对计量点的贡献。这两种方法虽然计算量小，但是都无法准确得到堆芯及其外部厂房三维中、光子通量分布。“离散纵标”可以准确得到三维中、光子通量分布，但是其计算量和内存消耗巨大，串行计算完全无法满足大规模屏蔽计算的要求。“蒙特卡罗”方法几何描述能力强，但是从堆芯到核电厂厂房通量变化梯度大，难以收敛，而且只能得到某些指定计量点的通量。目前，国内外能够一次性准确获得从核电厂厂房三维中、光子通量分布的计算方案非常少，因而需要针对核电厂厂房屏蔽计算中大计算量、深穿透的特点，提出速度快、精度高的计算方法。申请号为201610471985.5的申请提出了一种应用于反应堆堆芯屏蔽计算的并行扫描方法，该方法将三维几何进行二维的区域分解，提高了并行效率，但是每次计算时只含有一个角度，依然无法满足全厂房大规模屏蔽计算的需求。在已申请的专利(申请号为201610471985.5)基础上，本发明克服了单角度扫描的限制，实现了多角度同时扫描的高效并行计算方法，适合于超大规模的并行屏蔽计算。

发明内容

本发明的目的在于针对核电厂厂房屏蔽计算中，大计算量、深穿透的特点，提出一种用于核电厂厂房屏蔽计算的高效并行扫描方法，适用于大规模集群并行计算，以便于快速准确获得核电厂厂房的三维中、光子通量分布。

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

一种用于核电厂厂房屏蔽计算的高效并行扫描方法，该方法包括以下步骤：

步骤一：用细网格去逼近核电厂厂房屏蔽计算模型，得到三维计算网格；

步骤二：将三维计算网格分配到若干个CPU上，进行二维区域分解，每个CPU都贯穿Z方向上的所有计算网格，即得到二维网格子区域；

步骤三：选定二维网格子区域的几何中心点，按“米”字形将二维网格子区域划分成八个CPU区域；

步骤四：将八个卦限：(-,-,-)，(-,-,+)，(-,+,-)，(-,+,+)，(+,-,-)，(+,-,+)，(+,+,-)，(+,+,+)分为四组，即：①(-,-,-)和(-,-,+)，②(-,+,-)和(-,+,+)，③(+,-,-)和(+,-,+)，④(+,+,-)和(+,+,+)，同一组的两个卦限放在一起计算；

步骤五：根据每个CPU区域在二维网格子区域中的相对位置确定卦限的扫描顺序；

步骤六：这四组卦限分别从二维网格子区域的四个对角向中心扫描，当某个CPU的出射角通量在下游CPU中不会被立即使用时，则将该CPU的出射角通量存起来，否则利用阻塞通信将其发送给下游CPU；

步骤七：对于步骤六中未立即传递给下游CPU的出射角通量采取非阻塞通信传输，保证计算的效率和通信安全；

步骤八：记录每一个并行单元的最大迭代误差，并进行收敛性验证，若计算已经收敛则输出三维中、光子通量分布，若未收敛则重复步骤六和步骤七，直至计算收敛。

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

1.将核电厂厂房屏蔽计算需要的大内存分配到各个CPU上，可实现超大规模屏蔽计算。

2.二维网格子区域四个角点同时向中心扫描，提高了并行效率，随着CPU数据的增多，并行效率下降非常平缓，能够充分调用大规模集群计算资源，快速、准确地得到三维中、光子通量分布。

3.结合内部边界条件，利用阻塞与非阻塞混合通信，提高了通信的安全性和效率。

4.在通量变化梯度较大的情况下具有良好的收敛性。

附图说明

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

图2是二维网格子区域。

图3是二维网格子区域米字形划分。

图4是二维网格子区域米字形划分边界调整。

图5是6×5二维网格子区域多角度扫描示意图。

图6是两类CPU间的信息示意图。

图7是扫描过程中的内部边界示意图。

图8是6×5二维网格子区域多角度扫描示意图。

图9是二维网格子区域多方向扫描示意图。

具体实施方式

本发明通过细网格将连续的厂房屏蔽计算模型离散化，同时将细网区域分配到若干个CPU上，而每个CPU上的网格又沿轴向划分成若干并行单元，通过并行单元之间的通信，将所有CPU耦合在一起，将二维网格子区域进行“米”字形划分，划分后的每个CPU子区域都有自己既定的卦限扫描顺序，将所有卦限分成四组，且让四组卦限中的角度同时从二维网格子区域的四个角点同时向中心扫描，结合内部边界条件充分利用阻塞和非阻塞通信，从而保证并行计算的收敛性和高效性，该方法主要包括以下几个方面：

1)沿X方向、Y方向和Z方向划分细网格，要求所划分的细网格足以逼近真实的厂房屏蔽计算模型，如此，得到三维计算网格；

2)将三维计算网格沿X方向和Y方向划分为若干个子区域，每个子区域贯穿了Z方向的所有网格，得到二维网格子区域；同时，为了提高并行效率，将每个子区域沿轴向又划分成若干个并行单元；

3)将二维网格子区域按“米”字形进行分区，当出现分界线将某一网格划破时，需要对分界线进行微调，最后形成8个互不重叠的二维CPU子区域；

4)由于Z方向上的两个卦限在二维CPU子区域上具有完全相同的扫描顺序，因此将Z方向上的两个卦限串在一起扫描，即将八个卦限分成四组分别扫描计算；

5)四组卦限分别从二维网格子区域的四个角点同时向中心扫描，必然会出现冲突现象，因此需要规定每个二维CPU子区域的卦限扫描顺序；

6)设置内部边界来缓存不会被下游CPU立即使用的出射角通量，同时将下游CPU急需使用的出射角通量用阻塞通信发送出去；

7)将内部边界缓存起来的出射角通量用非阻塞通信发送出去，以防止通信的死锁和确保通信的高效性；

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

如此，对于每一个能群,都可以执行相同扫描迭代计算，直至所有能群计算完毕。

1 厂房屏蔽计算模型离散化

真实的三维核电厂厂房屏蔽计算模型在空间上是连续的，数值模拟计算需要将连续的模型离散化。本发明采用规则网格对厂房屏蔽计算模型离散化，即沿X方向、Y方向、Z方向划分细小的矩形网格，利用这些矩形网格去逼近真实的厂房模型。每一个网内部只允许有一种材料，这样就可以得到离散化后的三维网格计算模型。

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

本发明中的区域分解就是将三维网格划分成若干个二维网格子区域(每个子区域内部依然是三维)如图1所示，其中的“body”即为三维网格，“processor”即为一个二维子区域，该区域内的所有细网格由一个CPU负责计算。图1中将计算几何分解成了9个二维子区域，每个二维子区域都贯穿了轴向(Z方向)的所有网格。本发明中将所有二维子区域构成的二维平面(如图2所示)称之为二维网格子区域。而二维网格子区域中的一个网格定义为：一个CPU网格。为了使各CPU负载平衡，减少CPU间的等待时间，子区域尽量均匀划分。为了提高并行效率，本发明将每个二维子区域(processor)沿Z方向又划分成若干个并行单元，即图1中的“block”。

3 二维网格子区域的米字形划分

为了让四组卦限的角度同时由二维网格子区域的四个角落同时向中心扫描，需要将二维网格子区域进行“米”字形划分，如图3所示，对于一个6×5的二维网格子区域，选定其几何中心，然后分别沿水平方向、竖直方向以及±45度角方向作四条直线，将整个平面分成八个区域。从图6中可以看出，±45角度方向上的两条直线总是会将CPU网格划破，而对于水平和竖直两条直线，则取决于Y方向和X方向的CPU个数的奇偶性，当CPU个数为奇数时，CPU网格会被划破，而当CPU个数为偶数时，直线则恰好落在二维网格子区域CPU网格的内部边界上。

在并行计算时，要求每个CPU网格只属于“米”字形划分中的一个区域。因此需要对被划破的网格做特殊处理。首先，将每个CPU网格的两条对角线相连，两条对角线的交点即为该CPU网格的中心点，中心点落在哪个区域，该CPU网格就属于哪个区域，当某个CPU网格的中心点恰好落在“米”字形分割线上时，则依据如下卦限的正负性法则来判断：将一个卦限中的角度定义为(Ω_x，Ω_y，Ω_z)，当某个CPU网格的中心点恰好落在“米”字形分割线上时，就意味从不同卦限扫描过来的任务会在该CPU网格相遇，当两个不同卦限的任务相遇时，Ω_x>0的具有优先计算的权限，当Ω_x的符号相同时，Ω_y>0的具有优先计算的权限。综合上述的CPU网格中心点位置和卦限正负法则就可以将二维网格子区划分成八个互不重叠的区域，如图4所示，并按照逆时针方向从左下角开始，将“米”字形分区进行编号(1～8)。

4 卦限分组

在本发明中，区域分解本质上是将三维的网格划分若干个二维子区域，因为每个CPU都贯穿了Z方向上的所有网格，因此Z方向上的两个卦限对于二维网格子区域而言是具有完全相同的扫描顺序的，为了便于说明，现将八个卦限进行编号：①(+,+,+)②(+,-,+)③(-,+,+)④(-,-,+)⑤(+,+,-)⑥(+,-,-)⑦(-,+,-)⑧(-,-,-)，其中(+,+,+)和(+,+,-)扫描顺序相同，(+,-,+)和(+,-,-)扫描顺序相同，(-,+,+)和(-,+,-)扫描顺序相同，(-,-,+)和(-,-,-)扫描顺序相同。因此以扫描顺序为依据可以将这八个卦限分为四组，同一组内的所有角度一起扫描计算，直到一组结束时再切换到另一组扫描。

5 二维CPU子区域卦限扫描顺序确定

本发明中按照扫描的方向，下游的网格是否具备计算条件总是需要依赖上游的网格的计算是否完成，同时同一组内的卦限具有相同扫描顺序，按照上述规定的卦限顺序，各“米”字形CPU分区的卦限扫描顺序依次为：

1区：①→⑤→②→⑥→③→⑦→④→⑧；

2区：①→⑤→③→⑦→②→⑥→④→⑧；

3区：③→⑦→①→⑤→④→⑧→②→⑥；

4区：③→⑦→④→⑧→①→⑤→②→⑥；

5区：④→⑧→③→⑦→②→⑥→①→⑤；

6区：④→⑧→②→⑥→③→⑦→①→⑤；

7区：②→⑥→④→⑧→①→⑤→③→⑦；

8区：②→⑥→①→⑤→④→⑧→③→⑦；

图5给出了在每个卦限只有3个角度的情况下，6×5二维CPU子区域的扫描示意图，本发明允许每个卦限中有任意多个角度。依据既定的扫描顺序，可以对并行效率进行定量分析，并行效率的计算表达式如下：

式中：

η——理论并行效率；

T_s——串行计算时间；

T_p——并行计算时间；

P——所用CPU的个数。

并且，定义如下变量：

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

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

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

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

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

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

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

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

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

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

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

T_s＝MIJKt₀ (2)

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

对于每个CPU，都处于两种状态：计算或等待。若将每个CPU的计算步数和等待步数分别定义为N_work和N_idle，那么总的扫描步数为:N_work+N_idle,记为：N_total。在并行计算只有所有CPU都结束计算时，并行计算才算结束，因此每个CPU具有相同的N_total值，又每个CPU的N_work值均为：8M*K/K_c。因此，每个CPU具有相同的N_idle值。本发明所涉及的多角度扫描方法，每个CPU的等待步数可以表示为：

当I/I_c或J/J_c为奇数时，σ等于1，否则等于0，以确保能够被2整除。因此，对于一个CPU而言，其等待时间和计算时间可以分别表示为：

T_p2＝8MI_cJ_cKt₀ (6)

如此，多角度扫描方法并行效率可以表示为：

5 内部边界条件设置

本发明中对各CPU之间需要传递的信息分成了两类：不会被下游CPU立即使用的信息和下游CPU急需使用的信息，如图6所示，空心箭头为下游CPU急需使用的信息，实心箭头为不会被下游CPU立即使用的信息。在本发明中，将不会被下游CPU立即使用的信息存在内部边界中，当一个卦限的所有角度计算完毕时，再与相邻的CPU通信，而对于下游CPU急需使用的信息则在一个角度计算结束时便立即发送出去。

由于本发明中，四组卦限是同时从四个角点向中心扫描的，因此对于每个二维CPU子区域，两类信息都是动态变化的，这必然导致内部边界的设置也是动态变化的，图7展示了内部边界的动态变化过程，其中黑色粗线表示各CPU子区域的边界，点划线则表示内部边界。

6 混合通信模式

在本发明中，保证通信的安全性和高效性，对于两类信息分别采用了不同的通信方法，即对于不会被下游CPU立即使用的信息采用非阻塞通信以实现通信和计算的重叠，而对于下游CPU急需使用的信息采用阻塞通信，以保证迭代格式不会发生退化。阻塞通信方式在每一个角度计算完毕时，即将信息传递给相邻的下游CPU，简单高效，相比于阻塞通信，非阻塞通信则要复杂的多，下面以图8中的2×2的二维网格子区域为例，对本发明中的非阻塞通信方式进行详细描述。

分析图8中的0号CPU，它首先计算(+,+,+)卦限中的所有角度，当(+,+,+)卦限中的所有角度计算完毕时开始计算(+,+,-)卦限中的所有角度，同时接收来自1号CPU和2号CPU发送出来的信息，即来自(-,+,+)卦限和(+,-,+)卦限的入射角通量；接着0号CPU开始计算(+,-,+)卦限中的所有角度，并将(+,+,-)卦限的出射角通量发送给1号CPU和2号CPU，同时接收来自1号CPU和2号CPU发送出来的信息，即(-,+,-)卦限和(+,-,-)卦限的入射角通量；当(+,-,+)卦限中的所有角度计算完毕时开始计算(+,-,-)卦限中的所有角度，同时将(+,-,+)卦限的出射角通量发送给1号CPU；接着0号CPU开始计算(-,+,+)卦限中的所有角度，同时将(+,-,-)卦限的出射角通量发送给1号CPU，并且同时从1号CPU和2号CPU中接收(-,-,+)卦限的入射角通量；当(-,+,+)卦限中的所有角度计算完毕时，0号CPU开始计算(-,+,-)卦限中的所有角度，并且将(-,+,+)卦限的出射角通量发送给2号CPU同时从1号CPU和2号CPU接收(-,-,-)卦限的入射角通量；当(-,+,-)卦限中的所有角度计算完毕时，0号CPU开始计算(-,-,+)卦限中的所有角度，并同时将(-,+,-)卦限的出射角通量发送给2号CPU；最后0号CPU计算完(-,-,-)卦限中的所有角度，如此一次迭代计算完毕。

图9给出了本发明所涉及多角度并行扫描方法和已申请专利(申请号为201610471985.5，后面称之为传统方法)并行扫描方法从4个CPU到400个CPU的并行效率变化图，可以看出本发明所涉及的多角度并行方法随CPU个数的增加，优势越来越明显，在400CPU时，本发明的并行效率比传统方法高出约13％，且随着CPU个数的增多，本发明的并行效率下降的越来越慢，表明了本发明具有更好的并行扩展性，能够胜任超大规模并行计算。

8 收敛性验证

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

Claims (1)

1.一种用于核电厂厂房屏蔽计算的高效并行扫描方法，该方法包括以下步骤：

步骤一：用细网格去逼近核电厂厂房屏蔽计算模型，得到三维计算网格；

步骤二：将三维计算网格分配到若干个CPU上，进行二维区域分解，每个CPU都贯穿Z方向上的所有计算网格，即得到二维网格子区域；

步骤三：选定二维网格子区域的几何中心点，按“米”字形将二维网格子区域划分成八个CPU区域；

步骤四：将八个卦限进行编号：①(+,+,+)②(+,-,+)③(-,+,+)④(-,-,+)⑤(+,+,-)⑥(+,-,-)⑦(-,+,-)⑧(-,-,-)，将这八个卦限分为四组，分别为：(-,-,-)和(-,-,+)，(-,+,-)和(-,+,+)，(+,-,-)和(+,-,+)，(+,+,-)和(+,+,+)，同一组的两个卦限放在一起计算；

步骤五：根据每个CPU区域在二维网格子区域中的相对位置确定卦限的扫描顺序，下游网格是否具备计算条件总是需要依赖上游的网格的计算是否完成，同时同一组内的卦限具有相同扫描顺序，各“米”字形CPU分区的卦限扫描顺序依次为：

1区：①→⑤→②→⑥→③→⑦→④→⑧，

2区：①→⑤→③→⑦→②→⑥→④→⑧，

3区：③→⑦→①→⑤→④→⑧→②→⑥，

4区：③→⑦→④→⑧→①→⑤→②→⑥，

5区：④→⑧→③→⑦→②→⑥→①→⑤，

6区：④→⑧→②→⑥→③→⑦→①→⑤，

7区：②→⑥→④→⑧→①→⑤→③→⑦，

8区：②→⑥→①→⑤→④→⑧→③→⑦；

步骤六：这四组卦限分别从二维网格子区域的四个对角向中心扫描，当某个CPU的出射角通量在下游CPU中不会被立即使用时，则将该CPU的出射角通量存起来，否则利用阻塞通信将其发送给下游CPU；

步骤七：对于步骤六中未立即传递给下游CPU的出射角通量采取非阻塞通信传输，保证计算的效率和通信安全；

步骤八：记录每一个并行单元的最大迭代误差，并进行收敛性验证，若计算已经收敛则输出三维中、光子通量分布，若未收敛则重复步骤六和步骤七，直至计算收敛。

Images