技术分享｜并行代数多重网格算法：如何用黑盒求解器攻克复杂工程计算的效率瓶颈？
、
01 多重网格方法介绍
多重网格方法是一种高效求解偏微分方程离散系统的迭代方法，其核心思想是通过不同网格层次的协同作用加速收敛。它分为几何多重网格（Geometric Multigrid Method, GMG）和代数多重网格（Algebraic Multigrid Method, AMG）两类，分别基于几何信息和纯代数结构构建。
传统迭代方法如雅可比(Jacobi)、高斯-赛德尔(Gauss–Seidel)方法虽能在细网格上快速消除高频误差，但对低频误差效果不佳。多重网格方法通过将问题转移到粗网格，使低频误差在粗网格上变为高频，从而被有效消除。它构建一系列从细到粗的网格，通过限制(Restriction)和插值(Interpolation, or Prolongation)在不同网格间传递信息，利用粗网格快速消除低频误差，细网格修正高频误差，以此加速收敛。
下图反映了一个2D泊松方程的迭代求解过程中残差分布的变化（初始随机分布），模型分辨率为100 × 100个网格点，使用的迭代方法为高斯-赛德尔迭代法。可以发现，长波长(低频)残差的衰减速度要比短波长(高频)残差慢得多。
(图片来自文献Introduction to Numerical Geodynamic Modelling)
02 代数多重网格(AMG)方法
代数多重网格方法是一种用于求解稀疏线性方程组的高效数值计算方法，特别适用于工程和科学计算中的复杂问题。它通过将计算区域划分为多个层次化的网格，以提高计算效率和精度。AMG方法的基本思想是利用粗网格和细网格之间的关系，通过在不同层次上进行平滑和残差修正来加速求解过程。它结合了代数方法和多重网格技术，不需要显式的网格生成，而是直接在代数层面上操作，通过层次化拓扑关系的构建得到各层级的稀疏矩阵。这使得AMG方法具有较大的灵活性，可以适应不规则的几何形状和复杂的边界条件。
（一）主要特点
• AMG与GMG的对比
| 特性 | 几何多重网格（GMG） | 代数多重网格（AMG） |
| 依赖信息 | 几何网格结构 | 矩阵的代数结构 |
| 适用场景 | 结构化网格、规则问题 | 非结构化网格、复杂问题 |
| 粗网格构建 | 基于几何层次 | 基于矩阵的强连接性 |
| 插值/限制算子 | 基于几何插值（如线性插值） | 基于矩阵的数值关系构造 |
• 优势：灵活性强，可看作黑盒求解器，无需用户提供网格信息，适用场景广。
• 挑战：粗网格选择和插值算子构造的算法复杂度高，针对不同领域存在收敛性问题，预处理时间更多。
（二）计算流程
AMG的求解流程分为两个部分：
• 启动过程(Setup)：递归构建层次结构（）直至粗网格足够小；
• 求解过程(Solve)：通过V/W/F等循环(Cycle)在层次间迭代，直到收敛。
下图展示了AMG的计算流程，其中绿色箭头表示Setup过程，蓝色箭头表示Solve过程。在Setup阶段，AMG根据用户提供的原始矩阵，自动生成一系列规模逐渐减小的粗矩阵，原始矩阵作为最细层，不同算法在此阶段的表现各异。进入Solve阶段，基于生成的层次矩阵进行迭代计算。以常用的V-cycle为例，如蓝色箭头所示，计算从最细层网格开始，依次向上在各层网格进行，最后反向返回最细层。在每一层网格中，会依次执行smoother计算、残差计算以及层间的限制和插值计算，这些步骤的具体操作已在图中展示。
(图片来自文献FP16 Acceleration in Structured Multigrid Preconditioner for Real-World Applications)
代数多重网格（AMG）无需几何信息，直接基于系数矩阵的代数结构构建层次。其关键组件包括：
1. 延拓矩阵：定义细网格到粗网格的映射；
2. 限制矩阵：通常取；
3. Galerkin积：粗网格矩阵通过构造。
根据Setup阶段的算法差异，可将AMG分为聚合型和经典型两个派别，下面分别简单介绍。
（三）聚合型AMG
聚合(aggregation)型AMG方法采用基于聚类的粗化过程，它利用聚类算法和相似性度量来确定节点之间的关联性，然后将相似节点划分为同一组，即若干个细网格点聚合形成一个粗网格点。聚合型AMG方法的特点为具有经济性、易于控制等，但对于一些复杂问题其求解收敛性可能会不如经典AMG方法。根据粗细网格插值关系的处理方式可将其分为光滑聚合(Smoothed-Aggregation, SA)型AMG方法和非光滑聚合(Unsmoothed-Aggregation, UA)型AMG方法。代表性代数求解库为Trilinos的MueLu。
下图为聚合型AMG生成粗矩阵（用图(Graph)表示）的过程，细层0、1、2节点聚合形成粗层0节点，细层3、4、5节点聚合形成1节点，SA和UA两种AMG形成的延拓矩阵如图所示，延拓矩阵每列反映了每个粗层节点(聚合体)与相应细层节点的包含关系。
（四）经典型AMG
经典(classical)型AMG方法通过粗网格点和细网格点的划分来完成粗化过程。它根据残差向量的分布，通过确定连接矩阵中的强连接来识别网格节点之间的重要关系。通过保留具有强连接的节点和相关信息，可以构建粗网格和粗细网格间的插值关系。经典AMG方法的特点是鲁棒性较好，在各类实际应用中得到了充分发展，算法适应性不断增强，目前已经发展了多种并行的粗化策略，如Falgout算法、HMIS算法等。代表性代数求解库为HYPRE。
下图为经典型AMG生成粗矩阵（用图(Graph)表示）的过程，从细层中挑选出一些节点作为粗层节点(Coarse Points)，其相邻的节点为细层节点(Fine Points)，延拓矩阵每列反映了C-points和F-points间的插值关系。
03 UNAP介绍
UNAP(UNstructured Algebra Package)是无锡超算先进制造部针对国产异构众核平台开发的非结构网格代数求解库。 该代数求解支持亿级阶矩阵的百万核并行求解，包含Krylov子空间方法(CG、BiCGStab、GMRES...)及预条件子(Jacobi、ILU、MG...)、代数多重网格方法、并行直接解法等，软件支持神威、x86、Windows平台。
UNAP使用聚合型的AMG方法，以两点聚合(pairwise-aggregation)为主，相较于经典类型的AMG，聚合类型的AMG具有算法和实现简洁的特征，在Setup、单步solve等方面具有一定速度优势。
（一）使用方法
1. 构建矩阵、向量
向量构建过程比较简单，传入对应数组指针即可，支持多种形式内存管理，接口如下：
// \brief 从已有的数组中复 制指定长度的数据
// \param val 数据指针
// \param length 向量的大小
// \param comm 通信域
// \param reUse 数据是否会被重复使用
// \param dontDel 如果为真，则不会删除指针
template<typename T>
Vector::Vector(const T *val, const label &length, Communicator *comm,
const bool reUse = false, const int dontDel = 0);UNAP中的矩阵数据结构包含CSR与LDU两种(结构示意图如下)，均为分布式，包含多种矩阵构造接口，这里以CSR构造为例：
(图片来自文献An integrated framework for accelerating reactive flow simulation using GPU and machine learning models)
/*!
*通过三个数组直接构造分布式CSR矩阵
*\param rowPtr 反映当前进程矩阵各行元素数量的数组，长localRows_+1
*\param colInd 当前进程矩阵的列标数组，长localnnz_
*\param values 当前进程矩阵的数值数组，长localnnz_
*\param dist 各进程矩阵首行行标数组，长nProcs+1
*\param symm 矩阵结构是否对称（默认结构对称）
* \param comm 通信域
* \param deepCopy 深拷贝时为true,浅拷贝为false
*/
template<typename Cmpt>
CSRMatrix::CSRMatrix(label* rowPtr, label* colInd, Cmpt* values, label* dist, bool symm,
Communicator* comm, bool deepCopy false);2. 初始化求解器
以AMG方法为例，UNAP构造求解器的步骤如下：
// 构造AMG聚合工具类
OFAgglomeration<scalar> aggl(A);
// 设置聚合相关参数
aggl.SET_maxLevels(6);
aggl.SET_rowSizeInCoarsestLevel(50);
// AMG setup开始
aggl.agglomerate();
// AMG参数设置
MGSolver<scalar> MG(A, aggl);
// 收敛残差
MG.SET_relTol(1e-6);
// 最大迭代步
MG.SET_maxIter(100);
// 信息输出
MG.SET_ifPrint(false);
// 前/后光滑次数
MG.SET_nPreSweeps(1);
MG.SET_nPostSweeps(1);
MG.SET_maxPreSweeps(2);
MG.SET_maxPostSweeps(2);
// 光滑器类型，如Cheby、Jacobi、Gauss-Seidel等
MG.SET_smootherType("Jacobi"]);
//-1 为使用默认 V-cycle
MG.SET_KcycleLevel( -1);
// 某些情况加速收敛选项
MG.SET_scaleCorrection(false);
// 初始化光滑器
MG.initSmoothers();
// AMG setup完成
MG.setUp();3. 求解与结果输出
// 迭代数据管理类
SolverPerformance solverPerf = MG.solve(x, a, b);（二）两类AMG测试
对两类AMG方法进行了简单测试，聚合型使用UNAP库，经典型使用HYPRE库。
1. 算例说明
算例来自CFD仿真的实际工程案例所得到的数据，不可压多相湍流绕流计算，涵盖动量3分量（u、v、w）、压力修正值（p）、流体体积分数（s）、湍动能（k）、湍流耗散率（e）共7个变量的线性方程组。
算例1的网格数为589940，变量e矩阵非零元数为4049160、其他变量（k、s、p、u、v、w）的矩阵非零元数为4141924。
算例2的网格数为2849291，变量e矩阵非零元数为19357108、其他变量（k、s、p、u、v、w）的矩阵非零元数为20149067。
2. 测试环境
测试CPU为sw26010处理器（仅使用主核），编译选项只有-O3，UNAP的求解器参数设置如上示例所示，HYPRE的求解器参数设置如下：
HYPRE_BoomerAMGCreate(&solver);
HYPRE_BoomerAMGSetTol(solver, 1.e-6); /* conv. tolerance */
HYPRE_BoomerAMGSetMaxCoarseSize(solver, 50);
HYPRE_BoomerAMGSetMaxIter(solver, 100);
HYPRE_BoomerAMGSetMaxLevels(solver, 10)； /* maximum number of levels */
HYPRE_BoomerAMGSetPrintLevel(solver, 0);
HYPRE_BoomerAMGSetRelaxType(solver, 0); // jacobi
HYPRE_BoomerAMGSetCoarsenType(solver, 8); // PMIS-8 HMIS-10
HYPRE_BoomerAMGSetup(solver, parcsr_A, par_b, par_x);尽管在实际仿真求解过程中，针对 u、v、w 等分量的求解存在速度更快的迭代法可供选择，但为了统一比较代数多重网格（AMG）方法在不同变量求解中的性能表现，这里对所有变量的求解均采用 AMG 方法，同时保持一些共有参数的一致性（如收敛准则、最大层数、最大迭代步数、最粗层粒度、smoother类型等），这里HYPRE的粗化策略选择了PMIS并行粗化。
3. 测试结果
（1）算例1
①对于 u、v、w 等变量的求解，两种方法的收敛速度大致相当。然而，由于聚合型AMG方法具有更低的计算复杂度，因此在求解速度上展现出明显优势。
②在求解 k 变量时，聚合型 AMG 方法的表现显著优于经典型 AMG 方法；而在求解 s 变量时，聚合型 AMG 方法的表现却明显逊于经典型 AMG 方法。
③在给定的设置条件下，经典型 AMG 方法在求解 p 方程时未能实现收敛。
（2）算例2
①对于 u、v、w、e、k 等变量的求解，经典型 AMG 方法的求解耗时异常地长，尽管其迭代步数显示正常，具体原因尚需进一步分析。
②在求解 s 变量时，聚合型 AMG 方法的迭代步数明显多于经典型 AMG 方法；而在其他变量的求解中，两种方法的迭代步数相差不大。
③同样地，在该设置下，经典型 AMG 方法在求解 p 方程时也未能收敛。
（3）总结
通过上述两个算例的简要分析，可以发现聚合型 AMG 方法与经典型 AMG 方法在不同变量的求解过程中各具特点和优势。在某些变量的求解上，聚合型 AMG 方法凭借更快的求解速度和更低的计算复杂度脱颖而出；而在其他变量的求解中，经典型 AMG 方法则可能更胜一筹。
值得注意的是，在本文算例中，经典型 AMG 在求解 p 方程时未能实现收敛。这一现象暗示着针对不同问题的 p 方程求解，可能需要进一步的参数调整、优化乃至方法改进。同时，经典型 AMG 方法在算例 2 中对某些变量求解耗时异常的问题也引起了我们的关注。这一问题值得深入研究和分析，以提高其在不同场景下的稳定性和效率。
此外，我们还发现聚合型 AMG 方法在求解一些 p 方程等问题时，并行效率不尽如人意。这与网格分区和进程间独立聚合的算法有着紧密的联系，在一些问题中也会较大程度地影响收敛效率。因此，如何改进其聚合算法，以提升并行效率和收敛效率，成为了亟待解决的问题，值得我们深入探索和研究。
参考资料： [1] Gerya T. The multigrid method. In: Introduction to Numerical Geodynamic Modelling. Cambridge University Press; 2019:292-318. [2] Zong, Yi, et al. "FP16 Acceleration in Structured Multigrid Preconditioner for Real-World Applications." Proceedings of the 53rd International Conference on Parallel Processing. 2024.
