进行全堆芯三维精细化中子输运计算必须依赖于有效的加速手段，粗网有限差分(coarse mesh finite difference, CMFD) 方法^[1]因其便于实施、加速效果好，被广泛地应用在全堆芯输运计算的加速算法中。然而，三维全堆芯Pin尺度的CMFD是一个大型稀疏非对称线性方程组，其系数矩阵的规模可达到上亿级别，高效求解CMFD线性方程组对实现加速至关重要。广义极小残差算法是求解大型非对称线性方程组的优秀方法，但该方法的优越性取决于良好预处理技术的应用，特别是针对条件数很大或者严重病态的线性方程组，高效的预处理技术尤为重要^[2]。针对于串行计算，不完全LU分解(incomplete LU decomposition, ILU)、对称超松弛(symmetric over-relaxation, SOR)、块预处理、稀疏近似逆等是很好的预处理技术^[3]。但是针对于并行计算，需要采用红黑网格策略，才能充分发挥ILU、SOR的预处理效果^[4]，但同时带来的问题就是引入不必要的计算等待时间，以及使得程序开发变得异常复杂。目前在核反应堆计算领域，Block-Jacobi不完全LU分解(block-jacobi incomplete LU decomposition, BJILU) 因其实施方便而得到了广泛的应用，如美国的MPACT程序^[5]，但其仅对各自CPU中的元素进行了预处理，对不同核心间需要信息传递的元素并未做任何预处理。Xu^[6]提出了一种简化对称超松弛预处理技术(reduced symmetric over-relaxation, RSOR) 与不完全LU分解的混合预处理方法(RSOR and ILU, RSILU)，有效解决了上述问题，但在实际应用中发现该方法有待进一步优化。

本文采用以下2种方法对现有的RSILU预处理进行了优化：1)将不完全LU分解预处理子替换为修正不完全LU分解(modified incomplete LU decomposition, MILU)预处理子，以进一步提高RSILU的预处理效率；2)由严格计算对角块矩阵的逆改为近似计算，以解决RSILU方法在复杂能群结构的多群CMFD问题中，预处理计算耗时过多的问题。利用C5G7-3D基准题和VERA Problem #4基准题搭建CMFD线性方程组对优化后的RSILU方法进行了测试分析。所有测试均基于MPI并行编程模型^[7]，采用空间区域分解的方式进行并行计算。

1 CMFD线性方程组与GMRES算法 1.1 三维多群CMFD线性方程组

建立在Pin尺度上的三维多群CMFD方程组为：

$ \begin{aligned} &\left(\varSigma_{t, g, c}^{k} \phi_{g, c}^{k}-\sum\limits_{g^{\prime}=1}^{G} \varSigma_{s, g^{\prime} \rightarrow g, c}^{k} \phi_{g^{\prime}, c}^{k}\right) V_{c}+\left(\hat{D}_{g, c W}^{k}+\hat{D}_{g, c N}^{k}+\right. \\ &\left.\hat{D}_{g, c E}^{k}+\hat{D}_{g, c S}^{k}+\hat{D}_{g, c T}^{k}+\hat{D}_{g, c B}^{k}\right) \phi_{g, c}^{k}-\left(\hat{D}_{g, W c}^{k} \phi_{g, W}^{k}+\right. \\ &\hat{D}_{g, N c}^{k} \phi_{g, N}^{k}+\hat{D}_{g, E c}^{k} \phi_{g, E}^{k}+\hat{D}_{g, S c}^{k} \phi_{g, S}^{k}+ \\ &\left.\hat{D}_{g, T c}^{k} \phi_{g, T}^{k}+\hat{D}_{g, B c}^{k} \phi_{g, B}^{k}\right)=V_{c} \frac{\chi_{g, c}}{k_{\text {eff }}} \sum\limits_{g=1}^{G} \nu \varSigma_{f, g, c}^{k-1} \phi_{g, c}^{k-1} \end{aligned} $

(1)

式中：ϕ、Σ_t、Σ_s、vΣ_f、χ、k_eff分别是中子通量密度、总截面、散射截面、吸收截面、裂变谱和有效增殖因子。式(1)在数学上可表示为线性方程组：

$ \boldsymbol{A x}=\boldsymbol{b} $

(2)

式中：未知量x即为待求的中子通量ϕ，按照“先能群、后空间”的顺序进行排列；系数矩阵A ∈ R^n×n是七对角、稀疏、非对称矩阵，其中主对角线的对角块由能群散射矩阵构成：

$ \boldsymbol{A}=\left[\begin{array}{cccc} \boldsymbol{D}_{1} & \boldsymbol{U}_{A, 1} & & \\ \boldsymbol{L}_{A, 2} & \boldsymbol{D}_{2} & \ddots & \\ & \ddots & \ddots & \boldsymbol{U}_{A, n-1} \\ & & \boldsymbol{L}_{A, n} & \boldsymbol{D}_{n} \end{array}\right]=\boldsymbol{L}_{A}+\boldsymbol{D}+\boldsymbol{U}_{A} $

(3)

式中：D为对角块矩阵；L_A和U_A为严格非对角块矩阵。对于能群、空间耦合求解的CMFD线性方程组；D_i为G×G稠密矩阵(如图 1所示)；G为能群数；L_{A, i}和U_{A, i}为稀疏矩阵。对几何空间进行区域分解，并基于MPI进行分布式内存计算，此时系数矩阵A可进一步表示为：

$ \boldsymbol{A}=\boldsymbol{L}_{I}+\boldsymbol{L}_{P}+\boldsymbol{D}+\boldsymbol{U}_{P}+\boldsymbol{U}_{I} $

(4)

式中：L_P和U_P是存储在当前CPU内的严格非对角块，L_I和U_I是存储在不同CPU内的严格非对角块，如图 2所示。

1.2 预处理GMRES算法

广义极小残差算法(generalized minimum residual method, GMRES)是求解非对称线性方程组的有效方法，该方法是Krylov子空间法的一种，通过在子空间上进行投影以迭代的形式寻找近似解。记r₀ = b - Ax₀为初始残差向量，由r₀生成的Krylov子空间可表示为：

$ K_{m}\left(\boldsymbol{A}, \boldsymbol{r}_{0}\right)=\operatorname{span}\left\{\boldsymbol{r}_{0}, \boldsymbol{A} \boldsymbol{r}_{0}, \boldsymbol{A}^{2} \boldsymbol{r}_{0}, \cdots, \boldsymbol{A}^{m-1} \boldsymbol{r}_{0}\right\} $

(5)

通过K_m(A, r₀) 可得到一组规范正交基V_m+1=[v₁ v₂ … v_m v_m+1]，其中v_m ∈ Rⁿ可通过Gram-Schmidt正交化得到。同时还可得到一个上Hessenberg矩阵${\mathit{\boldsymbol{\hat H}}_{m + 1, m}}$，并满足：

$ \boldsymbol{A} \boldsymbol{V}_{m}=\boldsymbol{V}_{m+1} \hat{\boldsymbol{H}}_{m}=\boldsymbol{V}_{m} \boldsymbol{H}_{m}+\boldsymbol{h}_{m+1, m} \boldsymbol{v}_{m+1} \boldsymbol{e}_{m}^{\mathrm{T}} $

(6)

GMRES方法在Krylov子空间中产生一系列近似解，这些近似解逐步逼近于真解，同时这些近似解的残差向量r_m满足二范数最小的性质。r_m可表示为：

$ \begin{gathered} \left\|\boldsymbol{r}_{m}\right\|_{2}=\left\|\boldsymbol{b}-\boldsymbol{A} \boldsymbol{x}_{m}\right\|_{2}=\left\|\beta \boldsymbol{e}_{1}-\hat{\boldsymbol{H}}_{m} \boldsymbol{y}_{m}\right\|_{2}, \\ \beta=\left\|\boldsymbol{r}_{0}\right\|_{2} \end{gathered} $

(7)

经过m步GMRES算法形成近似解x_m满足：

$ \boldsymbol{x}_{m}=\boldsymbol{x}_{0}+\boldsymbol{V}_{m} \boldsymbol{y}_{m} $

(8)

其中y_m∈Rⁿ通过极小化式(9)得到，即：

$ \boldsymbol{y}_{m}=\operatorname{arg}\underset{y}{\operatorname{min}}\left\|\beta \boldsymbol{e}_{1}-\hat{\boldsymbol{H}}_{m} \boldsymbol{y}\right\|_{2} $

(9)

GMRES方法的优越性取决于良好的预处理技术的应用，特别是针对条件数很大或者严重病态的线性方程组，高效的预处理技术尤为重要。预处理的本质是对线性方程组(2)作同解变换，以右预处理为例：

$ \boldsymbol{A x}=\boldsymbol{b} \Leftrightarrow \boldsymbol{A M}^{-1} \boldsymbol{M x}=\boldsymbol{b} \Leftrightarrow \tilde{\boldsymbol{A}} \tilde{\boldsymbol{x}}=\boldsymbol{b} $

(10)

式中 M 为预处理子(预处理矩阵)。本文中采用右预处理GMRES算法求解CMFD线性方程组^[8]，右预处理GMRES算法如下：

算法1：右预处理GMRES算法

1)     选取初值x₀ ∈ Rⁿ，计算初始残差r₀ = b - Ax₀，定义β=‖r₀‖₂, v₁ = r₀/β

2)     For j=1, 2, …, m Do

3)     计算w: = AM^-1v_j

4)         For i = 1, 2, …, j Do

5)             计算h_{i, j}∶ = (w, v_i)

6)             计算w∶=w-h_{i, j}v_i

7)         End Do

8)         计算h_{j+1, j}=‖w‖₂, v_j+1= w/h_{j+1, j}

9)         定义V_m=[v₁, v₂, …, v_m], ${\mathit{\boldsymbol{\hat H}}_m} = \left\{ {{h_{i, j}}} \right\}_{1 \le j \le m}^{1 \le i \le j + 1}$

10)     End Do

11)     计算${\mathit{\boldsymbol{y}}_m} = \arg \mathop {\min }\limits_y {\left\| {\beta {\mathit{\boldsymbol{e}}_1} - {{\mathit{\boldsymbol{\hat H}}}_m}\mathit{\boldsymbol{y}}} \right\|_2}$

12)     计算${\mathit{\boldsymbol{x}}_m} = {\mathit{\boldsymbol{x}}_0} + {\mathit{\boldsymbol{M}}^{ - 1}}\left( {{\mathit{\boldsymbol{V}}_m}{\mathit{\boldsymbol{y}}_m}} \right)$

13)     如果满足收敛标准则停止；否则置x₀ : = x_m并转向1。

2 RSILU预处理算法及其优化 2.1 RSILU预处理子

预处理矩阵 M 的选取对GMRES收敛速率影响极大。文献[6]提出了高效的RSILU混合预处理方法，该方法由RSOR和ILU预处理子共同构成，不仅对各自CPU中的元素进行了预处理，还对不同CPU间需要信息传递的元素也进行了预处理，如图 3所示。其中，RSOR高效预处理不同CPU间需要信息传递的元素，RSOR预处理子为：

$ \boldsymbol{M}_{\mathrm{RSOR}}^{-1}=\boldsymbol{D}^{-1}\left[\boldsymbol{I}-\omega\left(\boldsymbol{L}_{I}+\boldsymbol{U}_{I}\right) \boldsymbol{D}^{-1}\right] $

(11)

式中ω为松弛因子。RSOR不需要红黑网格技术，便于实施，并且当核数的增加时GMRES所需的迭代次数保持不变，同时还能保证解的对称性，是一种简单高效的并行预处理算法，更多详细内容可参考文献[6]。

ILU高效预处理各自CPU中的元素，ILU(0)预处理为：

$ \boldsymbol{M}_{\mathrm{ILU}}=\left(\boldsymbol{L}_{P} \tilde{\boldsymbol{D}}^{-1}+\boldsymbol{I}\right)\left(\tilde{\boldsymbol{D}}+\boldsymbol{U}_{P}\right) $

(12)

式中$\mathit{\boldsymbol{\tilde D}}$为修正后的对角块，通过对 D 的对角元素进行修正得到。$\mathit{\boldsymbol{\tilde D}}$的对角元素计算式为：

$ \tilde{d}_{i, \mathrm{ILU}}=a_{i i}-\sum\limits_{j=1}^{i-1}\left[\frac{a_{i j}}{\tilde{d}_{j}} a_{j i}\right] $

(13)

结合式(11)和式(12)，RSILU预处理为：

$ \begin{aligned} \boldsymbol{M}_{\mathrm{RSLU}}^{-1}=&\left(\tilde{\boldsymbol{D}}+\boldsymbol{U}_{P}\right)^{-1}\left(\boldsymbol{L}_{P} \tilde{\boldsymbol{D}}^{-1}+\boldsymbol{I}\right)^{-1}[\boldsymbol{I}-\\ &\left.\omega\left(\boldsymbol{L}_{I}+\boldsymbol{U}_{I}\right) \boldsymbol{D}^{-1}\right] \end{aligned} $

(14)

为减少RSILU预处理实施过程中的存储和计算负担，使用D_E替换上式中的D，D为D_E的对角元素。最终RSILU预处理子可定义为：

$ \begin{aligned} \boldsymbol{M}_{\mathrm{RSLLU}}^{-1}=&\left(\tilde{\boldsymbol{D}}+\boldsymbol{U}_{P}\right)^{-1}\left(\boldsymbol{L}_{P} \tilde{\boldsymbol{D}}^{-1}+\boldsymbol{I}\right)^{-1}[\boldsymbol{I}-\\ &\left.\omega\left(\boldsymbol{L}_{I}+\boldsymbol{U}_{I}\right) \boldsymbol{D}_{E}{}^{-1}\right] \end{aligned} $

(15)

2.2 RSILU预处理算法的实施步骤

从算法1中可以看出，实施预处理的关键一步是如何计算z = M^-1v。为减少RSILU预处理实施过程中的存储和计算负担，RSILU预处理过程中并未显式地构造预处理矩阵M_RSILU，而是只对角块$\mathit{\boldsymbol{\tilde D}}$进行计算、存储，z 算法2过程为：

算法2：z = M^-1v

1) 针对每个CPU的边界元素，计算w₁ = D_E^-1v

2) 将w₁传递给其他CPU，以用于并行计算

3) 计算w₂ = v-ω(L_I+ U_I)w₁

4) 求解(L_P$\mathit{\boldsymbol{\tilde D}}$^-1+I)w₃ = w₂得到w₃

5) 求解($\mathit{\boldsymbol{\tilde D}}$+ U_P)z = w₃得到 z 。

2.3 优化的RSILU预处理算法

通过算法2可以看出，RSILU预处理子具有以下特点：1)串行计算环境下，RSILU便蜕化为ILU(0) 预处理子；2)并行计算环境下，RSILU不需要进行红黑排序，结构化和非结构化网格均适用；3)仅需额外存储对角块$\mathit{\boldsymbol{\tilde D}}$的LU分解矩阵，用于步骤3和步骤4中求解$\mathit{\boldsymbol{\tilde D}}$^-1x；4)仅需进行2次对角块矩阵的求解和1次不同CPU间的信息传递；

2.3.1 MILU替换ILU

针对RSILU的第1个优化是将ILU预处理子替换为改进的不完全LU分解MILU预处理子。在数学上，MILU预处理子具有比ILU更好预处理效果；在方法的实施上，MILU与ILU的唯一区别，仅仅在于$\mathit{\boldsymbol{\tilde D}}$对角元素的计算式略微有所不同，式(13)中的a_ji被替换为了$\sum\limits_{k = j + 1}^N {{a_{jk}}} $。MILU的$\mathit{\boldsymbol{\tilde D}}$对角元素为：

$ \tilde{d}_{i, \mathrm{MLLU}}=a_{i i}-\sum\limits_{j=1}^{i-1}\left[\frac{a_{i j}}{\tilde{d}_{i}} \sum\limits_{k=j+1}^{N} a_{j k}\right] $

(16)

2.3.2 对角块矩阵的近似求逆

RSILU的第2个优化是针对对角块矩阵的求逆运算。对于多群CMFD，在算法2中的步骤3和步骤4需要计算$\mathit{\boldsymbol{\tilde D}}$^-1x，常采用完全LU分解进行求逆运算。当能群结构较为复杂时，$\mathit{\boldsymbol{\tilde D}}$则是由散射矩阵构成的稠密矩阵。图 1中给出了HELIOS截面库生成的燃料栅元的47群散射矩阵，可以看出，若对$\mathit{\boldsymbol{\tilde D}}$进行完全LU分解来计算$\mathit{\boldsymbol{\tilde D}}$^-1x可能导致计算代价过大以及内存消耗过大。

近似计算$\mathit{\boldsymbol{\tilde D}}$^-1x则是一种折中的办法，例如：对$\mathit{\boldsymbol{\tilde D}}$进行不完全LU分解，或者仅保留$\mathit{\boldsymbol{\tilde D}}$的下三角元素，或者仅保留$\mathit{\boldsymbol{\tilde D}}$的对角元素。近似计算$\mathit{\boldsymbol{\tilde D}}$^-1x将减少RSILU单次预处理的计算时间，但也会导致RSILU预处理效果的降低，因此需要找到一个平衡点，这一点将在下一节的数值验证部分作进一步的讨论。

3 预处理效果验证

为验证优化后的RSILU算法的预处理效果，本文基于C语言开发了预处理GMRES求解器，并选用C5G7-3D基准题^[9-11]和VERA Problem #4基准题^[12]搭建pin尺度单群/多群CMFD线性方程组作为测试题(见表 1)，其中单群CMFD通过对多群CMFD进行能群归并得到^[9]，单群CMFD方程组的系数矩阵没有对角块，是一个传统的7对角矩阵。为了后续表述更简便，优化前的RSILU算法用RSILU-old表示，优化后的RSILU算法用RSILU-new表示。

预处理GMRES选用相对残差10^-8作为收敛标准。串行计算环境：Inter(R) Core(TM) i5-7200U CPU@2.71Hz, RAM 8.0GB；并行计算环境：“天河一号”超级计算机，采用商用InfiniteBand网络连接，每个计算节点包含28个计算核心(14个2×Intel Xeon CPU E5-2690 v4 @2.60 GHz)，RAM 128 GB。编译器选用英特尔编译器Intel-16.0.3；MPI编译环境选用mvapich2-2.2。

3.1 C5G7-3D基准题

本节选用C5G7-3D基准题来测试MILU预处理子对方法的优化效果；C5G7-3D基准题能群结构较简单，仅含有7个能群，因此在多群CMFD的预处理中可采用完全LU分解计算 $\mathit{\boldsymbol{\tilde D}}$^-1x；测试环境包括串行环境和并行环境。

3.1.1 串行计算

图 4中给出了串行环境下，RSOR、RSILU-old和RSILU-new预处理GMRES求解一次C5G7-3D单群/多群CMFD线性方程组的迭代次数和计算时间。由图 4可知：

1) 3种的预处理算子中，RSILU-new收敛最快，计算用时也最少，RSILU-old次之，最次是RSOR；

2) RSILU-new针对单群CMFD的预处理效果要优于多群CMFD。单群CMFD下RSILU-new预处理效果显著优于RSILU-old；多群CMFD下RSILU-new的预处理效果略优于RSILU-old。

3.1.2 并行计算

图 5和图 6给出了不同核数下BJILU、RSILU-old和RSILU-new预处理GMRES求解C5G7-3D单群/多群CMFD线性方程组的迭代次数和计算时间。图 7和图 8分别给出了不同核数下RSOR、BJILU、RSILU-old、RSILU-new、和串行ILU(SILU)预处理GMRES求解C5G7-3D单群/多群CMFD线性方程组的残差收敛历史。其中，SILU是指GMRES在串行环境下用ILU(0)作预处理时的收敛结果，用作对照，以显现出各种预处理算法在并行核数逐渐增加时的预处理效果的变化趋势。

不同预处理算法比较如图 5~8所示，可知：

1) 随着并行核数的增加，RSILU-old、RSILU-new和BJILU的预处理GMRES计算时间逐渐减少，GMRES迭代次数逐渐增加，但迭代次数随核数的增长地很缓慢；

2) 随着并行核数的增加，RSILU-old和RSILU-new的预处理效果的差距不断缩小，但RSILU-new效果始终不差于RSILU-old，且始终优于BJILU；

3) RSILU-new针对单群CMFD的预处理效果要优于多群CMFD，且RSILU-new针对单群CMFD的预处理效果甚至可以超过SILU (如图 7中的(a)~(c)所示)。

3.2 VERA problem #4基准题

本文选用VERA problem #4基准题来检验RSILU-new方法针对复杂能群结构多群CMFD的预处理效果，选用HELIOS截面库生成所需的47群截面，测试采用20个计算核心进行计算，下列结果均通过完整求解VERA problem #4获得。首先对比了不同$\mathit{\boldsymbol{\tilde D}}$^-1x计算方法对RSILU-new预处理效果的影响，测试方案包括：1) 对$\mathit{\boldsymbol{\tilde D}}$进行完全LU分解；2) 对$\mathit{\boldsymbol{\tilde D}}$进行不完全LU分解(ILU(0))；3) 仅保留$\mathit{\boldsymbol{\tilde D}}$的下三角元素；4) 仅保留$\mathit{\boldsymbol{\tilde D}}$的对角元素。

测试结果如图 9所示。可以看出方案d效果最佳。一方面，GMRES迭代次数并未显著增加(仅从900次增加到913次)，使得迭代计算的用时相差不大；另一方面，方案4预处理耗时大幅降低，约为方案1的1/2左右。以上2点使得方案4整体的耗时显著减少，因此在后续测试中，RSILU-new均采用仅保留$\mathit{\boldsymbol{\tilde D}}$的对角元素的近似求解方案。

BJILU、RSILU-old和RSILU-new预处理效果的测试结果如表 2和图 10所示。从结果中可以看出：不论是单群CMFD还是多群CMFD、BJILU和RSILU-old相比，RSILU-new预处理效果都是最优的，不仅GMRES的迭代次数是最少的，而且迭代计算和预处理计算的耗时也是最少的。其中，RSILU-new的预处理计算耗时仅是RSILU-old的1/2左右。整体GMRES的计算耗时相比于优化之前减少30%。

4 结论

1) 优化后的RSILU弥补了该方法之前的一些缺陷，进一步提高了RSILU预处理并行GMRES方法求解大规模CMFD线性方程组的计算效率。

2) 本文中的并行GMRES算法仅基于MPI编程模型，未来可开发MPI和OpenMP混合并行编程技术进一步减少处理器间的通信时间，从而提高并行效率。