先进中子学栅格程序KYLIN-Ⅱ输运模块并行优化开发

涂晓兰，柴晓明，刘 东，芦 韡，王 鑫，李勋昭，付元光，郭凤晨

(1.中国核动力研究设计院 核反应堆系统设计技术重点实验室，四川 成都 610213；2.北京应用物理与计算数学研究所，北京 100094)

中国核动力研究设计院不仅承担了压水堆的设计研发，同时承担大量新型反应堆的研发，不同反应堆的燃料组件的结构、燃料成分复杂，为能满足设计研发的需求，中国核动力研究设计院开发形成了先进中子学栅格程序KYLIN-Ⅱ[1-3]，KYLIN-Ⅱ采用子群方法[4]求解共振核素的共振截面，采用特征线方法[5-11]进行复杂几何中子输运计算，并采用广义粗网格有限差分加速方法(GCMFD)[12]来加速中子输运求解流程，采用基于改进预估修正临界-燃耗迭代方法(PPC)的切比雪夫方法[13]求解复杂燃耗链，同时，为方便用户使用，开发形成了支持复杂组件的图形化建模工具[14]和后处理显示工具。针对1个典型的二维AFA3G组件，其网格规模为4 000，能群为190群，特征线模块计算时间可达7 000 s，计算时间较长，需对输运模块进行并行优化，提升组件计算的效率。本文首先介绍输运计算模型，然后对KYLIN-Ⅱ程序的特征线输运模块进行效率分析和并行特征分析，针对热点计算模块，提出相应的并行策略及优化策略，最后通过基准题和典型的压水堆问题验证本文方法的正确性和并行优化的效率。

1 二维几何下特征线法求解输运方程

真实情况下中子输运问题是三维问题，对于轴向均匀的问题，可将计算问题退化为二维问题，然而在求解二维问题时，必须保证计算结果与三维计算结果等效。三维情况下特征线方法及求解公式见文献[5]。

二维情况下沿射线的通量ψ求解公式为：

(1)

其中：g为能群编号；i为网格编号；m为方位角编号；n为极角编号；k为特征线编号；s′n为射线在二维情况下的长度；θn为极角；Q为源项；Σt为总截面。

网格i的标通量计算公式为：

ψi,m,n,k(si,m,n,k))Δdi,m,k

(2)

其中：ωn为极角权重；ωm为方位角权重；Σt,i为总截面；Ai为网格面积；Δdi,m,k为射线段的宽度。

在求解中子输运问题时，各向异性散射会对计算结果有较大的影响[15]，KYLIN-Ⅱ程序截面库包含了高阶Pn散射截面，在特征线中子输运模块中，考虑了各向异性散射，故式(2)中包含了各向异性散射源项，KYLIN-Ⅱ程序对散射源项中的散射截面中的角度应用勒让德多项式进行展开，展开后的散射源项公式为：

Qscat,g,i,k=

(3)

其中：Σs,n,i,g′→g为散射截面；Pl为勒让德函数；Ψi,g′为角通量；Ω为角度。

射线布置及追踪方式影响特征线方法处理复杂几何问题的能力，KYLIN-Ⅱ程序综合应用了循环射线布置及追踪技术以及反向延长追踪技术，循环射线布置及追踪技术用于处理特定的问题，计算效率高；反向延长追踪技术用于处理任意几何问题，几何适应性强。

2 特征线输运模块并行优化方案设计

针对图1所示的典型压水堆二维AFA3G组件例题，对KYLIN-Ⅱ程序的输运模块进行效率分析，KYLIN-Ⅱ程序特征线输运模块中特征线扫描、高阶源项更新较耗时，因此本文主要针对特征线扫描、高阶源项更新进行优化。

图1 AFA3G组件几何网格划分Fig.1 Mesh division of AFA3G assembly

2.1 高阶散射源计算优化策略

KYLIN-Ⅱ程序应用式(3)计算高阶散射源，从公式可知能群g、方向k的散射源来自其他所有能群g′和方向k′角通量的贡献，将是一个关于g、k、g′、k′的4重循环，循环次数较多，使得计算效率较低，同时必须存储角通量数据以备计算，但会使用大量内存，如网格数10 000、能群数200、方向数100，用双精度存储角通量会耗去约1.5 GB内存；且在能群并行中，需通信角通量，而由于数据量大，会显著增加通信时间，降低整体并行效率，因此本文考虑将角通量应用球谐函数展开。

角通量应用球谐函数展开的公式为：

(4)

(5)

式(3)中的勒让德函数应用球谐函数表示：

(6)

将式(4)和(6)代入式(3)，利用球谐函数的正交性，可得：

(7)

利用正交性可得通量矩的计算公式：

(8)

因此在特征线扫描的过程中不需存储角通量，仅需存储通量矩即可，与角通量相比，通量矩的散射阶远小于角度离散数，对于L阶散射，计算共需通量矩的个数为(L+1)2，一般情况下，(L+1)2较K小1～2个量级，可节省内存，同时减少散射源计算的循环次数。

cosj(φk-φk′)

(9)

可见结果一定是实数，没有虚部。因此构造一个函数Zlj，满足：

(10)

则有：

(11)

通过构造Zlj函数，在计算过程中不需计算和存储复数。

2.2 基于MPI的能群并行策略

由于输运模块完全通过外迭代驱动收敛，即取消了散射源内迭代，所有能群1次特征线扫描结束后同时更新散射源和裂变源，因此能群之间没有数据依赖，且各能群之间的计算量严格相等，可采用能群并行，并建立能群并行通信域，同时在高阶散射源更新过程中也涉及能群的循环，因此在特征线扫描和高阶散射源更新部分均应用MPI进行能群的并行，并行流程如图2所示，首先将宏观截面等数据从主CPU广播到其他CPU；然后进行任务分配，当核数等于能群数或核数被能群数整除时，将能群均分到核上，达到理想的负载均衡；特征线扫描完成后对网格标通量、通量矩以及粗网净中子流进行数据通信；GCMFD求解过程中调用PETSC函数库进行并行求解，求解结束后需对粗网格通量和粗网格源项进行数据通信。

图2 特征线模块计算流程图Fig.2 Calculation flow of method of characteristics module

对特征线扫描部分进行性能分析，E指数(式(1))计算时间占比较高，约占特征线扫描时间的40%，因此本文对E指数进行优化，即在第1次计算时，将E指数进行存储，由于在特征线扫描过程中应用MPI并行模式对能群并行，E指数也将分能群存储，即每个CPU只存储当前分配能群的E指数，因此并行计算时不会有太大的内存消耗。

3 数值验证

3.1 C5G7 UO2组件基准题

C5G7 UO2组件基准题是从C5G7 2D基准题中提取的UO2组件，几何布置如图3所示，它的群常数参见OECD发布的基准题文档。计算时的网格规模为6 936，能群数为7，临界计算keff收敛准则是10-5，通量收敛准则是10-4，特征线方法计算条件为射线间距为0.01，0～2π区间方位角数为40，0～π/2区间方位角满足式(12)条件进行布置，0～π区间6个高斯勒让德极角。

m=1,2,…,N

φm∈(0,π/2)

(12)

其中：a为问题长度；b为问题宽度；N为0～π/2区间方位角数目。

图3 C5G7基准题的UO2组件Fig.3 Geometry of C5G7 UO2 assembly

表1所列为并行优化后的计算结果，由于其截面不含高阶散射，因此此基准题角通量球谐函数展开不起作用，其kinf计算结果小数点后前12位数与串行计算结果保持一致，说明本文的E指数优化和能群并行计算精度很高。应用E指数优化具有较高的计算效率，将特征线扫描速度提升75%。并行核数为7，加速比为6.37，并行的效率为90%，加E指数优化后的并行效率提升更好，较大减少了特征线扫描的时间，对闭循环特征线的扫描过程效率提升明显。

表1 C5G7 UO2组件计算结果Table 1 Calculation result of C5G7 UO2 assembly

表1列出的并行前、后迭代次数保持一致，由于输运模块计算过程中取消了散射源内迭代，所有能群一次特征线扫描结束后同时更新散射源和裂变源，能群之间没有数据依赖，因此能群的并行并不会导致迭代的退化。

3.2 AFA3G燃料组件基准题

针对AFA3G燃料组件进行正确性测试和性能测试，AFA3G燃料组件含有呈17×17方形排列的289个栅元，包括264根燃料棒、24根锆合金导向管和1根锆合金测量管。燃料棒栅距为1.26 cm，燃料棒芯块半径为0.409 6 cm，燃料包壳外半径为0.475 cm；导向管内半径为0.560 5 cm，外半径为0.622 5 cm；中心导向管内半径为0.572 5 cm，外半径为0.622 5 cm。本算例仅考虑热态、不含可燃毒物的情况下，燃料富集度为3.0%时的计算，导向管和测量管内部全部填充轻水慢化剂。

计算时的网格规模为4 105，能群数为190，临界计算keff收敛准则是10-5，通量收敛准则是10-4，并行计算的核数为190，特征线方法计算条件为射线间距为0.01，0～2π区间方位角数目为40，0～π/2区间方位角满足式(12)条件进行布置，0～π区间6个高斯勒让德极角。

表2列出了串行优化计算结果，kinf计算结果基本一致，角通量球谐函数展开对高阶散射源更新效率提升非常好，将高阶散射源速度提升87%，E指数优化提升了特征线扫描的效率，将特征线扫描速度提升77%，E指数优化和角通量球谐函数展开同时应用，整体计算时间效率得到了较好的提升，特征线模块的整体速度提升83%。

表3列出了并行计算结果，kinf计算结果基本一致，并行前、后迭代次数保持一致，若只是对能群并行计算效率并不高，加E指数后特征线扫描并行效率有了更好的提升，加角通量球谐函数展开对散射源更新的并行效率有了更大的提升，综合应用E指数优化和球谐函数展开后的并行效果较好，使原计算时间7 148.76 s缩短为38.85 s。

表2 AFA3G燃料组件串行优化计算结果Table 2 Serial optimization result of AFA3G fuel assembly

表3 AFA3G燃料组件并行优化计算结果Table 3 Parallel optimization result of AFA3G fuel assembly

4 结论

本文针对KYLIN-Ⅱ程序输运模块中的特征线扫描和高阶散射源计算应用MPI进行了能群并行，应用角通量球谐函数展开方法对高阶散射源计算进行了优化，应用通量矩取代了角通量，提高了并行通信效率以及高阶散射源的计算效率，并针对特征线扫描计算的E指数进行了优化存储，提高了特征线计算效率，由于每个CPU只存储当前分配能群的E指数，因此并行计算时不会有太大的内存消耗。并对多个基准题进行了验证，结果表明，E指数优化对特征线扫描加速效果显著，角通量球谐函数展开对高阶散射源计算加速效果显著，通过能群并行再结合E指数优化和角通量球谐函数展开，输运计算模块加速效果显著，在相同计算精度下，大幅提高了计算速度。