基于SA-AMG的彈塑性有限元計算的并行實現

張 倩 張健飛

(河海大學力學與材料學院 江蘇 南京 211100)

0 引 言

當前，數值模擬理論已逐步完善并被廣泛應用于科學、工程、生產等領域，而有限元法是其重要組成部分。彈塑性材料在實際中最為普遍，對其進行計算分析主要采用的就是彈塑性有限元法。與彈性有限元法相比，由于彈塑性有限元法的計算通常需要增量加載和迭代求解，因此其所需計算量和存儲量要更多。隨著工程規模的擴大和復雜性的增加，彈塑性有限元法對計算資源的要求越來越高，傳統的串行算法已不再滿足要求。隨著超級計算機的持續開發，我國在硬件技術方面已經達到了世界先進水平，但在軟件技術方面仍存在一定的問題。因此，為了擴大彈塑性有限元法的計算規模，充分利用超級計算機的計算能力，研究并開發了并行性和可擴展性均較高的彈塑性有限元并行程序。

有限元法的并行計算主要分為兩種[1-2]：一是基于系統方程求解的方法，通常需要形成整體系統方程，且只能對求解部分并行化，所以存儲量大、整體并行度不高；二是基于區域分解的方法，可以分為重疊型和非重疊型，無需形成整體系統方程，且各子區域獨立計算，所以存儲量小、整體并行度高。因此，有限元法的并行計算主要采用的是基于區域分解的方法。然而，從計算的角度來講，有限元法的核心部分是求解方程組Ax=b，其并行求解器主要分為兩類[3]：一是并行直接求解器，由于稀疏矩陣分解會導致非零填充，從而引起存儲量和計算量的增加，且并行度有限，所以其只適合于求解中小規模問題；二是并行迭代求解器，由于避免了非零填充，所以其存儲量和計算量小、并行度高。因此，對有限元法方程組的并行計算主要采用的是并行迭代求解器。

目前，彈塑性有限元法的并行計算[4-6]主要采用的求解方法是：先使用增量-牛頓法將非線性求解問題轉換為迭代的線性求解問題，再使用預條件共軛梯度法對線性化后的方程組進行并行計算。其中，PCG的預條件可分為兩種[7-10]：一是根據經典迭代法構造的預條件，如Jacobi預條件、多項式預條件等，其易于并行化，但迭代收斂的效果不理想；二是根據結合圖著色的不完全分解形成的預條件，如ILU預條件、ICC預條件等，其屬于直接法，故需要較大的存儲量、計算量，且并行度較低。SA-AMG作為一種迭代法，其不僅易于并行化，而且具有很好的收斂性和可擴展性。

本文基于Trilinos軟件包，利用增量-牛頓法、PCG和SA-AMG，實現了一種彈塑性問題的有限元并行求解方法，開發了相應的并行程序。通過對算例的計算，驗證了算法和程序的正確性，分析了光滑聚集代數多重網格法不同的聚集策略、光滑算法和粗網格求解方法對計算性能的影響，測試了程序的并行性和可擴展性。

1 有限元法及其并行性分析

當利用有限元法對實際工程問題計算分析時，通常按照離散化、單元分析和整體分析的步驟進行，逐步形成單元層面的剛度矩陣Ke、等效結點荷載列陣Fe、結點位移列陣δe，以及整體層面的剛度矩陣K、等效結點荷載列陣F和結點位移列陣δ等。其中：離散化是將結構體劃分成由有限個單元組成的網格模型，各單元間只通過相交的結點連接，并對其施加材料屬性、荷載和邊界條件等；單元分析是在每個單元中利用虛功原理或靜力等效原則計算出其各自的Ke和Fe；整體分析是在每個結點處利用力的平衡條件將單元按照原始結構重新組合，由Ke和Fe分別累加得到K和F，從而形成有限元方程：

Kδ=F

(1)

對于彈塑性材料模型，由于其具有非線性的應力-應變關系，所以彈塑性有限元法的方程組是非線性的，即：

K(δ)δ=F

(2)

式中：K(δ)表示K中的所有元素均可用δ中的相應元素表示。

對式 (2) 的計算通常使用增量-牛頓法：首先，設置荷載增量，將所施加的總荷載分成幾段，并對其循環；然后，在每個循環中，使用牛頓法迭代；最后，在每次迭代中，對線性方程組并行求解。因此，通過兩層循環迭代，就將求解非線性方程組轉化為了求解線性方程組。求解得到未知量δ后，根據彈塑性力學理論即可求出其應變ε和應力σ。

根據有限元法的基本原理，在單元分析中Ke和Fe的求解僅利用本單元信息，所以只要將同一個單元的信息存儲在同一個進程中，其計算就能夠完全并行；在整體分析中K和F的求解是將單元進行循環并按照結點編號的順序分別累加的，相互之間不需要通信，因此通過合理安排循環順序就可以提高其并行性。

2 光滑聚集代數多重網格預條件共軛梯度法

本文對增量-牛頓法每次循環迭代中線性化后的方程組使用光滑聚集代數多重網格預條件共軛梯度法并行求解。其中，PCG[11]能有效求解線性代數方程組，已被廣泛應用于有限元法的計算中；SA-AMG[12-13]是由多重網格法衍生而來的，但其僅依靠代數信息即可構造多重網格基本組件，核心思想是：在粗細網格層上均求解代數方程組，并分別用于消除低頻和高頻誤差。

2.1 預條件共軛梯度法

對于線性方程組Ax=b使用PCG求解時，按以下步驟進行：

1) 假設解x的初始值為x0，令其殘差為r0=b-Ax0，精確求解Mh0=r0，令p0=h0；

2) 當k=1,2,…時，進行如下迭代：

(3)

精確求解Mhk+1=rk+1，接著進行如下迭代：

(4)

上述算法中，M為預條件矩陣，其通過光滑聚集代數多重網格法近似求解，從而形成了SA-AMG預條件共軛梯度法。

2.2 光滑聚集代數多重網格法

1) 細網格前光滑：

對Ahuh=fh做μ1次松弛迭代：uh←Sμ1uh+

2) 粗網格校正：

粗網格方程求解：AHuH=fH；

3) 細網格后光滑：

對Akuk=fk做μ2次松弛迭代：uh←Sμ2uh+

3 程序實現

本文程序主要基于Trilinos開發實現。Trilinos[16]是由Sandia實驗室開發的大型數值軟件包，其致力于更加便利地對數學軟件庫進行設計、開發、集成和支持， 目標是在面向對象的軟件框架下開發解決大規模復雜物理工程和科學應用的并行算法和數學庫。其中：ML庫定義了一類基于多重網格法的預條件器；AztecOO求解器定義了一系列線性方程組的計算方法，包括預條件共軛梯度法；Epetra庫定義了各種矩陣、向量和圖表的構造和使用，支持串行、并行計算和分布式存儲。

為了充分利用現有的串行彈塑性有限元程序資源，本文程序采用C++與FORTRAN混合編寫。主程序利用C++編寫，用于調用MPI、Metis以及Trilinos相關操作，實現增量-牛頓法對彈塑性有限元問題的循環迭代求解，其中方程組求解部分使用SA-AMG預條件共軛梯度法；子程序利用FORTRAN編寫，用于執行與單元分析相關的計算。

在主程序中，首先讀入有限元網格模型的相關信息，調用Metis將所有單元分解為幾個子區域并分配給各個進程；然后各個進程分別調用彈性有限元FORTRAN子程序，計算彈性狀態下的Ke并形成分布式存儲的K；最后使用增量-牛頓法求解彈塑性有限元問題，進行荷載增量循環、牛頓法迭代：在每次迭代中，先計算Fe并形成分布式存儲的F，調用Trilinos中的ML庫和AztecOO求解器建立SA-AMG預條件共軛梯度法來并行計算；接著更新F并檢驗其是否收斂，如收斂則退出迭代，如不收斂則繼續迭代，各個進程再分別調用彈塑性有限元FORTRAN子程序，計算彈塑性狀態下的Ke并形成分布式存儲的K和F。

3.1 整體剛度矩陣和整體等效結點荷載列陣的形成

在有限元法的單元分析中，其Ke和Fe的求解可以完全并行。利用Metis將所有單元分解成幾個子區域并分配給各個進程后，一個子區域就對應一個進程，各個進程通過對其子區域中的單元進行循環，調用FORTRAN子程序，就可得到Ke、Fe和單元結點自由度列陣G。接著，根據G即可確定Ke和Fe中的元素分別在K和F中的位置，將其累加形成分布式存儲的K和F。整體剛度矩陣K是按照Epetra庫中的矩陣模式Epetra_FECrsMatrix來定義的，其采用了分布式稀疏行存儲格式， 是一種專門針對有限元計算的矩陣存儲格式。 整體等效結點荷載列陣F是按照Epetra庫中的向量模式Epetra_Vector來定義的，其采用了分布式稠密存儲格式。

3.2 求解方法的實現

在有限元法的整體分析中，分布式存儲的K和F均求解完成后，便可形成線性方程組，通過調用Trilinos的ML庫和AztecOO求解器對其進行并行計算。

求解部分程序的主要步驟為：首先，調用Epetra_LinearProblem定義每次迭代中線性化后的方程組Kx=F；然后，調用AztecOO求解器求解；接著，設置由ML庫提供的光滑聚集代數多重網格預條件MLPrec的主要參數，如聚集策略、光滑算法和粗網格求解方法，并調用ML_Epetra::MultiLevelPreconditioner將參數組合；最后，將預條件MLPrec添加到AztecOO求解器中，形成SA-AMG預條件共軛梯度法，對線性化后的方程組進行并行計算，并設置AztecOO求解器的參數，如輸出信息、最大迭代次數和收斂誤差等。

4 數值實驗

為檢驗算法和程序的正確性，分析光滑聚集代數多重網格法的主要參數對計算性能的影響，測試程序的并行性和可擴展性，在天河二號超級計算機上對程序進行了算例計算。

本算例采用如圖1所示的厚壁圓筒三維結構，其內徑10 mm，外徑20 mm，長100 mm，兩端各結點在X、Y、Z三方向均固定，其余結點在Z方向固定。內表面施加130 MPa的壓力荷載，并將其劃分為4個荷載增量。計算時均采用Von.Mises屈服準則，相關材料參數設置如下：屈服強度σY=380 MPa；彈性模量E=200 GPa；泊松比ν=0.3。利用Abaqus進行離散化，將厚壁圓筒結構劃分成3種不同的有限元網格模型，即網格模型一：結點數32 340、單元數28 800、總自由度數62 040；網格模型二：結點數241 920、單元數228 000、總自由度數473 760；網格模型三：結點數1 889 280、單元數1 833 600、總自由度數3 739 202。

圖1 算例示意圖

4.1 算法和程序的驗證

通過比較本文并行彈塑性有限元程序和既有串行彈塑性有限元程序在相同的網格模型下的計算結果，來驗證算法和程序的正確性。采用網格模型一進行計算分析。根據力學基礎知識可知，厚壁圓筒中段各點的位移應最大，所以從Z=50 mm的截面選取了如圖2所示的代表結點，并分別提取了各點在X、Y方向上的位移UX、UY，如表1所示。

圖2 Z=50 mm截面

mm

根據表1可以看出，本文并行彈塑性有限元程序和既有串行彈塑性有限元程序計算結果基本相同，因此算法和程序是正確的。

4.2 光滑聚集代數多重網格法的參數測試與分析

在利用Trilinos的ML庫進行光滑聚集代數多重網格法計算時，其參數主要有3個：(1) 聚集策略，其控制粗化的過程，主要有MIS、Uncoupled和Uncoupled-MIS 3種類型；(2) 光滑算法，其控制光滑計算的過程，主要有Jacobi、Gauss-Seidel、Chebyshev和Aztec 4種類型；(3) 粗網格求解方法，其控制粗網格求解的過程，主要有Jacobi、Gauss-Seidel、Chebyshev和Amesos-KLU 4種類型。為分析以上各個參數對計算性能的影響，采用網格模型一進行測試分析，所需的求解時間變化如圖3所示。由于光滑算法中Jacobi、Gauss-Seidel無論和誰組合，其收斂速度均較慢，在牛頓法迭代過程中，均未收斂，因此未包含在圖3中。

圖3 SA-AMG的不同參數組合對求解時間的影響

根據圖3可以看出，3種聚集策略中所需求解時間最少的是Uncoupled-MIS；4種光滑算法中所需求解時間最少的是Aztec；4種粗網格求解方法中所需求解時間最少的是Amesos-KLU。因此，將3種參數組合計算效果最優的為Uncoupled-MIS+Aztec+Amesos-KLU。

4.3 程序并行性和可擴展性的測試與分析

本程序大致分為3大部分：首先，讀入有限元網格模型的數據并利用Metis對所有單元進行分解，再將數據和分解結果廣播到所有進程上；然后，分別計算各單元彈性狀態下的Ke，并累加形成初始K；最后，在所有進程上進行增量-牛頓法的并行求解。為測試程序的并行性和可擴展性，分別對3種不同的有限元網格模型進行了不同進程數的并行計算，其中SA-AMG的主要參數采用了最優組合Uncoupled-MIS+Aztec+Amesos-KLU。統計了數據讀入和區域劃分時間、增量-牛頓法迭代求解時間和總運行時間，如表2-表4所示。根據表2-表4中的數據，又分別計算并繪制了增量-牛頓法迭代求解和總運行的并行加速比與進程數的關系，如圖4-圖6所示。

表2 網格模型一計算時間 s

表3 網格模型二計算時間 s

表4 網格模型三計算時間 s

圖4 網格模型一并行加速比與進程數的關系

圖5 網格模型二并行加速比與進程數的關系

圖6 網格模型三并行加速比與進程數的關系

根據表2-表4可以看出，數據讀入和區域劃分時間占總時間比例很小。對于相同的網格模型，隨著進程數的增加，其稍有減少，這是因為進程間的數據通信量減少了。但對于相同的進程數，隨著網格模型的擴大，其顯著增多，這是因為區域分解復雜度、數據讀入量和進程間的數據通信量都增加了。增量-牛頓法循環求解時間占總運行時間的絕大部分。對于相同的進程數，隨著網格規模的擴大，其顯著增多，這是因為計算量和數據通信量都有所增加。但對于相同的網格模型，隨著進程數的增加，其顯著減少，呈明顯下降趨勢，加速效果顯著。

根據圖4-圖6可以看出，對于相同的網格模型，隨著進程數的增加，其運行時間逐漸減少，但并行加速比并不呈線性，并行性能有所下降。這是因為每個進程中的計算量雖然減少了，但其子區域間的共享信息增加了，故進程間的數據通信量增加了，從而導致通信開銷變多，致使并行性能逐漸變差。因此，對于特定的問題，計算時應充分考慮多方面因素的影響，從而選擇最優的進程數和節點數進行計算，以充分發揮程序的并行性能。對于相同的進程數，隨著網格模型的擴大，其運行時間雖有所增加，但其增加速度遠小于網格模型的增大速度，因此，程序具有較好的可擴展性。

5 結 語

本文基于Trilinos的ML庫和AztecOO求解器實現了一種彈塑性問題的有限元并行求解方法和程序，使用增量-牛頓法循環迭代求解非線性問題，其中SA-AMG預條件共軛梯度法為線性求解器。通過與既有串行彈塑性有限元程序計算結果的對比，驗證了算法和程序是正確的；通過分析不同的光滑聚集代數多重網格法的主要參數對計算性能的影響，得到了計算效果最佳的組合為Uncoupled-MIS+Aztec+Amesos-KLU；通過測試不同的有限元網格模型的并行計算，證明了程序具有較好的并行性和可擴展性，可應用于大規模實際問題。