張 倩 張健飛

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

0 引 言

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

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

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

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

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

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

Kδ=F

(1)

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

K(δ)δ=F

(2)

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

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

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

2 光滑聚集代數(shù)多重網(wǎng)格預(yù)條件共軛梯度法

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

2.1 預(yù)條件共軛梯度法

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

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

2) 當(dāng)k=1,2,…時，進(jìn)行如下迭代：

(3)

精確求解Mhk+1=rk+1，接著進(jìn)行如下迭代：

(4)

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

2.2 光滑聚集代數(shù)多重網(wǎng)格法

1) 細(xì)網(wǎng)格前光滑：

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

2) 粗網(wǎng)格校正：

粗網(wǎng)格方程求解：AHuH=fH；

3) 細(xì)網(wǎng)格后光滑：

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

3 程序?qū)崿F(xiàn)

本文程序主要基于Trilinos開發(fā)實(shí)現(xiàn)。Trilinos[16]是由Sandia實(shí)驗(yàn)室開發(fā)的大型數(shù)值軟件包，其致力于更加便利地對數(shù)學(xué)軟件庫進(jìn)行設(shè)計(jì)、開發(fā)、集成和支持， 目標(biāo)是在面向?qū)ο蟮能浖蚣芟麻_發(fā)解決大規(guī)模復(fù)雜物理工程和科學(xué)應(yīng)用的并行算法和數(shù)學(xué)庫。其中：ML庫定義了一類基于多重網(wǎng)格法的預(yù)條件器；AztecOO求解器定義了一系列線性方程組的計(jì)算方法，包括預(yù)條件共軛梯度法；Epetra庫定義了各種矩陣、向量和圖表的構(gòu)造和使用，支持串行、并行計(jì)算和分布式存儲。

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

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

3.1 整體剛度矩陣和整體等效結(jié)點(diǎn)荷載列陣的形成

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

3.2 求解方法的實(shí)現(xiàn)

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

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

4 數(shù)值實(shí)驗(yàn)

為檢驗(yàn)算法和程序的正確性，分析光滑聚集代數(shù)多重網(wǎng)格法的主要參數(shù)對計(jì)算性能的影響，測試程序的并行性和可擴(kuò)展性，在天河二號超級計(jì)算機(jī)上對程序進(jìn)行了算例計(jì)算。

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

圖1 算例示意圖

4.1 算法和程序的驗(yàn)證

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

圖2 Z=50 mm截面

mm

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

4.2 光滑聚集代數(shù)多重網(wǎng)格法的參數(shù)測試與分析

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

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

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

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

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

表2 網(wǎng)格模型一計(jì)算時間 s

表3 網(wǎng)格模型二計(jì)算時間 s

表4 網(wǎng)格模型三計(jì)算時間 s

圖4 網(wǎng)格模型一并行加速比與進(jìn)程數(shù)的關(guān)系

圖5 網(wǎng)格模型二并行加速比與進(jìn)程數(shù)的關(guān)系

圖6 網(wǎng)格模型三并行加速比與進(jìn)程數(shù)的關(guān)系

根據(jù)表2-表4可以看出，數(shù)據(jù)讀入和區(qū)域劃分時間占總時間比例很小。對于相同的網(wǎng)格模型，隨著進(jìn)程數(shù)的增加，其稍有減少，這是因?yàn)檫M(jìn)程間的數(shù)據(jù)通信量減少了。但對于相同的進(jìn)程數(shù)，隨著網(wǎng)格模型的擴(kuò)大，其顯著增多，這是因?yàn)閰^(qū)域分解復(fù)雜度、數(shù)據(jù)讀入量和進(jìn)程間的數(shù)據(jù)通信量都增加了。增量-牛頓法循環(huán)求解時間占總運(yùn)行時間的絕大部分。對于相同的進(jìn)程數(shù)，隨著網(wǎng)格規(guī)模的擴(kuò)大，其顯著增多，這是因?yàn)橛?jì)算量和數(shù)據(jù)通信量都有所增加。但對于相同的網(wǎng)格模型，隨著進(jìn)程數(shù)的增加，其顯著減少，呈明顯下降趨勢，加速效果顯著。

根據(jù)圖4-圖6可以看出，對于相同的網(wǎng)格模型，隨著進(jìn)程數(shù)的增加，其運(yùn)行時間逐漸減少，但并行加速比并不呈線性，并行性能有所下降。這是因?yàn)槊總€進(jìn)程中的計(jì)算量雖然減少了，但其子區(qū)域間的共享信息增加了，故進(jìn)程間的數(shù)據(jù)通信量增加了，從而導(dǎo)致通信開銷變多，致使并行性能逐漸變差。因此，對于特定的問題，計(jì)算時應(yīng)充分考慮多方面因素的影響，從而選擇最優(yōu)的進(jìn)程數(shù)和節(jié)點(diǎn)數(shù)進(jìn)行計(jì)算，以充分發(fā)揮程序的并行性能。對于相同的進(jìn)程數(shù)，隨著網(wǎng)格模型的擴(kuò)大，其運(yùn)行時間雖有所增加，但其增加速度遠(yuǎn)小于網(wǎng)格模型的增大速度，因此，程序具有較好的可擴(kuò)展性。

5 結(jié) 語

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