01 多重網(wǎng)格方法介紹

多重網(wǎng)格方法是一種高效求解偏微分方程離散系統(tǒng)的迭代方法，其核心思想是通過不同網(wǎng)格層次的協(xié)同作用加速收斂。它分為幾何多重網(wǎng)格（Geometric Multigrid Method, GMG）和代數(shù)多重網(wǎng)格（Algebraic Multigrid Method, AMG）兩類，分別基于幾何信息和純代數(shù)結(jié)構(gòu)構(gòu)建。

傳統(tǒng)迭代方法如雅可比(Jacobi)、高斯-賽德爾(Gauss–Seidel)方法雖能在細(xì)網(wǎng)格上快速消除高頻誤差，但對低頻誤差效果不佳。多重網(wǎng)格方法通過將問題轉(zhuǎn)移到粗網(wǎng)格，使低頻誤差在粗網(wǎng)格上變?yōu)楦哳l，從而被有效消除。它構(gòu)建一系列從細(xì)到粗的網(wǎng)格，通過限制(Restriction)和插值(Interpolation, or Prolongation)在不同網(wǎng)格間傳遞信息，利用粗網(wǎng)格快速消除低頻誤差，細(xì)網(wǎng)格修正高頻誤差，以此加速收斂。

下圖反映了一個2D泊松方程的迭代求解過程中殘差分布的變化（初始隨機(jī)分布），模型分辨率為100 × 100個網(wǎng)格點(diǎn)，使用的迭代方法為高斯-賽德爾迭代法。可以發(fā)現(xiàn)，長波長(低頻)殘差的衰減速度要比短波長(高頻)殘差慢得多。

(圖片來自文獻(xiàn)Introduction to Numerical Geodynamic Modelling)

02 代數(shù)多重網(wǎng)格(AMG)方法

代數(shù)多重網(wǎng)格方法是一種用于求解稀疏線性方程組的高效數(shù)值計算方法，特別適用于工程和科學(xué)計算中的復(fù)雜問題。它通過將計算區(qū)域劃分為多個層次化的網(wǎng)格，以提高計算效率和精度。AMG方法的基本思想是利用粗網(wǎng)格和細(xì)網(wǎng)格之間的關(guān)系，通過在不同層次上進(jìn)行平滑和殘差修正來加速求解過程。它結(jié)合了代數(shù)方法和多重網(wǎng)格技術(shù)，不需要顯式的網(wǎng)格生成，而是直接在代數(shù)層面上操作，通過層次化拓?fù)潢P(guān)系的構(gòu)建得到各層級的稀疏矩陣。這使得AMG方法具有較大的靈活性，可以適應(yīng)不規(guī)則的幾何形狀和復(fù)雜的邊界條件。

（一）主要特點(diǎn)

? AMG與GMG的對比

特性	幾何多重網(wǎng)格（GMG）	代數(shù)多重網(wǎng)格（AMG）
依賴信息	幾何網(wǎng)格結(jié)構(gòu)	矩陣的代數(shù)結(jié)構(gòu)
適用場景	結(jié)構(gòu)化網(wǎng)格、規(guī)則問題	非結(jié)構(gòu)化網(wǎng)格、復(fù)雜問題
粗網(wǎng)格構(gòu)建	基于幾何層次	基于矩陣的強(qiáng)連接性
插值/限制算子	基于幾何插值（如線性插值）	基于矩陣的數(shù)值關(guān)系構(gòu)造

? 優(yōu)勢：靈活性強(qiáng)，可看作黑盒求解器，無需用戶提供網(wǎng)格信息，適用場景廣。

? 挑戰(zhàn)：粗網(wǎng)格選擇和插值算子構(gòu)造的算法復(fù)雜度高，針對不同領(lǐng)域存在收斂性問題，預(yù)處理時間更多。

（二）計算流程

AMG的求解流程分為兩個部分：

? 啟動過程(Setup)：遞歸構(gòu)建層次結(jié)構(gòu)（）直至粗網(wǎng)格足夠小；

? 求解過程(Solve)：通過V/W/F等循環(huán)(Cycle)在層次間迭代，直到收斂。

下圖展示了AMG的計算流程，其中綠色箭頭表示Setup過程，藍(lán)色箭頭表示Solve過程。在Setup階段，AMG根據(jù)用戶提供的原始矩陣，自動生成一系列規(guī)模逐漸減小的粗矩陣，原始矩陣作為最細(xì)層，不同算法在此階段的表現(xiàn)各異。進(jìn)入Solve階段，基于生成的層次矩陣進(jìn)行迭代計算。以常用的V-cycle為例，如藍(lán)色箭頭所示，計算從最細(xì)層網(wǎng)格開始，依次向上在各層網(wǎng)格進(jìn)行，最后反向返回最細(xì)層。在每一層網(wǎng)格中，會依次執(zhí)行smoother計算、殘差計算以及層間的限制和插值計算，這些步驟的具體操作已在圖中展示。

(圖片來自文獻(xiàn)FP16 Acceleration in Structured Multigrid Preconditioner for Real-World Applications)

代數(shù)多重網(wǎng)格（AMG）無需幾何信息，直接基于系數(shù)矩陣的代數(shù)結(jié)構(gòu)構(gòu)建層次。其關(guān)鍵組件包括：

1. 延拓矩陣：定義細(xì)網(wǎng)格到粗網(wǎng)格的映射；

2. 限制矩陣：通常取；

3. Galerkin積：粗網(wǎng)格矩陣通過構(gòu)造。

根據(jù)Setup階段的算法差異，可將AMG分為聚合型和經(jīng)典型兩個派別，下面分別簡單介紹。

（三）聚合型AMG

聚合(aggregation)型AMG方法采用基于聚類的粗化過程，它利用聚類算法和相似性度量來確定節(jié)點(diǎn)之間的關(guān)聯(lián)性，然后將相似節(jié)點(diǎn)劃分為同一組，即若干個細(xì)網(wǎng)格點(diǎn)聚合形成一個粗網(wǎng)格點(diǎn)。聚合型AMG方法的特點(diǎn)為具有經(jīng)濟(jì)性、易于控制等，但對于一些復(fù)雜問題其求解收斂性可能會不如經(jīng)典AMG方法。根據(jù)粗細(xì)網(wǎng)格插值關(guān)系的處理方式可將其分為光滑聚合(Smoothed-Aggregation, SA)型AMG方法和非光滑聚合(Unsmoothed-Aggregation, UA)型AMG方法。代表性代數(shù)求解庫為Trilinos的MueLu。

下圖為聚合型AMG生成粗矩陣（用圖(Graph)表示）的過程，細(xì)層0、1、2節(jié)點(diǎn)聚合形成粗層0節(jié)點(diǎn)，細(xì)層3、4、5節(jié)點(diǎn)聚合形成1節(jié)點(diǎn)，SA和UA兩種AMG形成的延拓矩陣如圖所示，延拓矩陣每列反映了每個粗層節(jié)點(diǎn)(聚合體)與相應(yīng)細(xì)層節(jié)點(diǎn)的包含關(guān)系。

（四）經(jīng)典型AMG

經(jīng)典(classical)型AMG方法通過粗網(wǎng)格點(diǎn)和細(xì)網(wǎng)格點(diǎn)的劃分來完成粗化過程。它根據(jù)殘差向量的分布，通過確定連接矩陣中的強(qiáng)連接來識別網(wǎng)格節(jié)點(diǎn)之間的重要關(guān)系。通過保留具有強(qiáng)連接的節(jié)點(diǎn)和相關(guān)信息，可以構(gòu)建粗網(wǎng)格和粗細(xì)網(wǎng)格間的插值關(guān)系。經(jīng)典AMG方法的特點(diǎn)是魯棒性較好，在各類實(shí)際應(yīng)用中得到了充分發(fā)展，算法適應(yīng)性不斷增強(qiáng)，目前已經(jīng)發(fā)展了多種并行的粗化策略，如Falgout算法、HMIS算法等。代表性代數(shù)求解庫為HYPRE。

下圖為經(jīng)典型AMG生成粗矩陣（用圖(Graph)表示）的過程，從細(xì)層中挑選出一些節(jié)點(diǎn)作為粗層節(jié)點(diǎn)(Coarse Points)，其相鄰的節(jié)點(diǎn)為細(xì)層節(jié)點(diǎn)(Fine Points)，延拓矩陣每列反映了C-points和F-points間的插值關(guān)系。

03 UNAP介紹

UNAP(UNstructured Algebra Package)是無錫超算先進(jìn)制造部針對國產(chǎn)異構(gòu)眾核平臺開發(fā)的非結(jié)構(gòu)網(wǎng)格代數(shù)求解庫。 該代數(shù)求解支持億級階矩陣的百萬核并行求解，包含Krylov子空間方法(CG、BiCGStab、GMRES...)及預(yù)條件子(Jacobi、ILU、MG...)、代數(shù)多重網(wǎng)格方法、并行直接解法等，軟件支持神威、x86、Windows平臺。

UNAP使用聚合型的AMG方法，以兩點(diǎn)聚合(pairwise-aggregation)為主，相較于經(jīng)典類型的AMG，聚合類型的AMG具有算法和實(shí)現(xiàn)簡潔的特征，在Setup、單步solve等方面具有一定速度優(yōu)勢。

（一）使用方法

  1. 構(gòu)建矩陣、向量

向量構(gòu)建過程比較簡單，傳入對應(yīng)數(shù)組指針即可，支持多種形式內(nèi)存管理，接口如下：

// \brief 從已有的數(shù)組中復(fù)制指定長度的數(shù)據(jù)
// \param val 數(shù)據(jù)指針
// \param length 向量的大小
// \param comm 通信域
// \param reUse 數(shù)據(jù)是否會被重復(fù)使用
// \param dontDel 如果為真，則不會刪除指針
template<typename T>
Vector::Vector(const T *val, const label &length, Communicator *comm, 
        const bool reUse = false, const int dontDel = 0);

UNAP中的矩陣數(shù)據(jù)結(jié)構(gòu)包含CSR與LDU兩種(結(jié)構(gòu)示意圖如下)，均為分布式，包含多種矩陣構(gòu)造接口，這里以CSR構(gòu)造為例：

(圖片來自文獻(xiàn)An integrated framework for accelerating reactive flow simulation using GPU and machine learning models)

/*!
 *通過三個數(shù)組直接構(gòu)造分布式CSR矩陣
 *\param rowPtr 反映當(dāng)前進(jìn)程矩陣各行元素數(shù)量的數(shù)組，長localRows_+1
 *\param colInd 當(dāng)前進(jìn)程矩陣的列標(biāo)數(shù)組，長localnnz_
 *\param values 當(dāng)前進(jìn)程矩陣的數(shù)值數(shù)組，長localnnz_
 *\param dist 各進(jìn)程矩陣首行行標(biāo)數(shù)組，長nProcs+1
 *\param symm 矩陣結(jié)構(gòu)是否對稱（默認(rèn)結(jié)構(gòu)對稱）
 * \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構(gòu)造求解器的步驟如下：

// 構(gòu)造AMG聚合工具類
 OFAgglomeration<scalar> aggl(A);
// 設(shè)置聚合相關(guān)參數(shù)
aggl.SET_maxLevels(6);
aggl.SET_rowSizeInCoarsestLevel(50);
// AMG setup開始
aggl.agglomerate();

// AMG參數(shù)設(shè)置
MGSolver<scalar> MG(A,  aggl);
// 收斂殘差
MG.SET_relTol(1e-6);
// 最大迭代步
MG.SET_maxIter(100);
// 信息輸出
MG.SET_ifPrint(false);
// 前/后光滑次數(shù)
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 為使用默認(rèn) V-cycle
MG.SET_KcycleLevel( -1);
// 某些情況加速收斂選項(xiàng)
MG.SET_scaleCorrection(false);
// 初始化光滑器
MG.initSmoothers();
// AMG setup完成
MG.setUp();

3. 求解與結(jié)果輸出

// 迭代數(shù)據(jù)管理類
SolverPerformance solverPerf = MG.solve(x, a, b);

（二）兩類AMG測試

對兩類AMG方法進(jìn)行了簡單測試，聚合型使用UNAP庫，經(jīng)典型使用HYPRE庫。

1. 算例說明

算例來自CFD仿真的實(shí)際工程案例所得到的數(shù)據(jù)，不可壓多相湍流繞流計算，涵蓋動量3分量（u、v、w）、壓力修正值（p）、流體體積分?jǐn)?shù)（s）、湍動能（k）、湍流耗散率（e）共7個變量的線性方程組。

算例1的網(wǎng)格數(shù)為589940，變量e矩陣非零元數(shù)為4049160、其他變量（k、s、p、u、v、w）的矩陣非零元數(shù)為4141924。

算例2的網(wǎng)格數(shù)為2849291，變量e矩陣非零元數(shù)為19357108、其他變量（k、s、p、u、v、w）的矩陣非零元數(shù)為20149067。

2. 測試環(huán)境

測試CPU為sw26010處理器（僅使用主核），編譯選項(xiàng)只有-O3，UNAP的求解器參數(shù)設(shè)置如上示例所示，HYPRE的求解器參數(shù)設(shè)置如下：

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);

盡管在實(shí)際仿真求解過程中，針對 u、v、w 等分量的求解存在速度更快的迭代法可供選擇，但為了統(tǒng)一比較代數(shù)多重網(wǎng)格（AMG）方法在不同變量求解中的性能表現(xiàn)，這里對所有變量的求解均采用 AMG 方法，同時保持一些共有參數(shù)的一致性（如收斂準(zhǔn)則、最大層數(shù)、最大迭代步數(shù)、最粗層粒度、smoother類型等），這里HYPRE的粗化策略選擇了PMIS并行粗化。

3. 測試結(jié)果

（1）算例1

①對于 u、v、w 等變量的求解，兩種方法的收斂速度大致相當(dāng)。然而，由于聚合型AMG方法具有更低的計算復(fù)雜度，因此在求解速度上展現(xiàn)出明顯優(yōu)勢。

②在求解 k 變量時，聚合型 AMG 方法的表現(xiàn)顯著優(yōu)于經(jīng)典型 AMG 方法；而在求解 s 變量時，聚合型 AMG 方法的表現(xiàn)卻明顯遜于經(jīng)典型 AMG 方法。

③在給定的設(shè)置條件下，經(jīng)典型 AMG 方法在求解 p 方程時未能實(shí)現(xiàn)收斂。

（2）算例2

①對于 u、v、w、e、k 等變量的求解，經(jīng)典型 AMG 方法的求解耗時異常地長，盡管其迭代步數(shù)顯示正常，具體原因尚需進(jìn)一步分析。

②在求解 s 變量時，聚合型 AMG 方法的迭代步數(shù)明顯多于經(jīng)典型 AMG 方法；而在其他變量的求解中，兩種方法的迭代步數(shù)相差不大。

③同樣地，在該設(shè)置下，經(jīng)典型 AMG 方法在求解 p 方程時也未能收斂。

（3）總結(jié)

通過上述兩個算例的簡要分析，可以發(fā)現(xiàn)聚合型 AMG 方法與經(jīng)典型 AMG 方法在不同變量的求解過程中各具特點(diǎn)和優(yōu)勢。在某些變量的求解上，聚合型 AMG 方法憑借更快的求解速度和更低的計算復(fù)雜度脫穎而出；而在其他變量的求解中，經(jīng)典型 AMG 方法則可能更勝一籌。

值得注意的是，在本文算例中，經(jīng)典型 AMG 在求解 p 方程時未能實(shí)現(xiàn)收斂。這一現(xiàn)象暗示著針對不同問題的 p 方程求解，可能需要進(jìn)一步的參數(shù)調(diào)整、優(yōu)化乃至方法改進(jìn)。同時，經(jīng)典型 AMG 方法在算例 2 中對某些變量求解耗時異常的問題也引起了我們的關(guān)注。這一問題值得深入研究和分析，以提高其在不同場景下的穩(wěn)定性和效率。

此外，我們還發(fā)現(xiàn)聚合型 AMG 方法在求解一些 p 方程等問題時，并行效率不盡如人意。這與網(wǎng)格分區(qū)和進(jìn)程間獨(dú)立聚合的算法有著緊密的聯(lián)系，在一些問題中也會較大程度地影響收斂效率。因此，如何改進(jìn)其聚合算法，以提升并行效率和收斂效率，成為了亟待解決的問題，值得我們深入探索和研究。

參考資料：

[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.

?