基于PETSc的有限元高性能求解方法

高懷玉，張 峰，秦忠國

(1.河海大學力學與材料學院，江蘇 南京 210098;2.南京市水利規(guī)劃設計院有限責任公司工程分院，江蘇南京 210006)

隨著工程技術領域的不斷開拓和發(fā)展，對大型復雜結(jié)構的分析要借助高階數(shù)值分析模型［1］和大規(guī)模數(shù)學計算，以保證數(shù)值解的精度［2］。在普通的計算設備上采用傳統(tǒng)的串行計算面臨著因計算資源不足而無法計算以及計算時間過長、工作效率低的問題。使用高性能計算機進行并行有限元分析是解決上述難題的一個方案［3-4］。自從美國國家航空航天局NASA的Noor等［4］在1975年發(fā)表第1篇有限元并行計算的論文以來，有限元并行算法的研究越來越具體和廣泛，適用于不同類型計算機的有限元并行算法也相繼產(chǎn)生。20世紀90年代以來，單處理機的計算速度不斷提高，并行計算機的體系結(jié)構趨于成熟，數(shù)據(jù)傳輸網(wǎng)絡的標準化和傳輸速率的大幅度提升，使得并行計算機的研制周期從幾年縮短至幾個月，為研制并行計算機系統(tǒng)創(chuàng)造了有利條件［5-7］。高性能計算機的發(fā)展勢頭異常迅猛，擁有并行計算功能的有限元軟件［8］也相繼出現(xiàn)。李根國等［9］將NASTRAN在“神威I”超級計算機系統(tǒng)上進行移植和并行計算功能的二次開發(fā)，獲得了成功。一些優(yōu)秀的并行數(shù)值算法包也日益引起關注。借助并行算法包開發(fā)有限元并行程序無需考慮復雜的數(shù)據(jù)分布和通信，大大降低了有限元并行計算的難度和成本，縮短了開發(fā)周期［10］。

可移植可擴展科學計算工具箱PETSc(portable，extensible toolkit for scientific computation)由一套數(shù)據(jù)結(jié)構和程序組成，由美國Argonne國家實驗室承擔開發(fā)并獲得成功，主要用于分布式存儲環(huán)境下高效求解偏微分方程組及相關問題的求解［3］。PETSc在軟件實現(xiàn)上采用ANSI C和消息傳遞標準接口MPI［11-12］。針對單處理機多核計算機，筆者利用PETSc提供的高層應用程序開發(fā)平臺開發(fā)出一種線性并行求解方法，并應用編寫的相應求解程序?qū)η蠼馄鬟M行了測試。

1 PETSc的數(shù)據(jù)結(jié)構

PETSc的方程組求解器組件定義在各種數(shù)據(jù)結(jié)構組件之上，包括線性方程組求解器KSP、非線性方程組求解器SNES、時間依賴方程求解器TS。PETSc基于基礎線性代數(shù)子程序庫BLAS、線性代數(shù)包LAPACK、MPI等實現(xiàn)計算功能。PETSc為用戶提供了豐富的Krylov子空間迭代方法和預條件子［13-14］，并具有錯誤檢測、性能統(tǒng)計和圖形打印等功能?；赑ETSc提供的大量對象和求解庫，用戶可以靈活地開發(fā)自己的應用程序，還可隨意添加和完善某些功能。

PETSc的數(shù)據(jù)結(jié)構包括向量Vec、矩陣Mat、網(wǎng)格管理對象、索引與排序、結(jié)構網(wǎng)格DA以及無結(jié)構網(wǎng)格的IS對象。用戶在PETSc環(huán)境下編寫串行或并行應用程序時，要按照固定的框架結(jié)構來編寫，即對于一般的應用對象均需要聲明、創(chuàng)建、設置、使用和銷毀。

2 基于PETSc的有限元求解器概況

基于PETSc的有限元并行求解器(以下簡稱求解器)以C++類［15］的形式存在，求解器的類封裝了一般線性有限元分析所必需的數(shù)據(jù)及操作，用戶按照給定的格式提供數(shù)據(jù)，依次調(diào)用成員函數(shù)，即可完成一個有限元的并行分析，得到位移解。使用求解器時，可以通過改變變量nfreedom值來改變有限元分析的類別，如nfreedom值設為1時，對應的是一維有限元問題，nfreedom值為2和3則分別對應平面和空間有限元問題。自定義的單元結(jié)構體StructElement，對外提供一個單元剛度矩陣計算程序接口，用戶可以自由地將自定義單元接入程序。求解器的類能方便地擴展單元類，具有一定的通用性和可擴展性，它繼承了PETSc工具箱所提供的運行時選項(runtime option)、性能分析等功能，即在求解時可以改變求解程序的迭代方法、預條件方法等設置，這些功能為測試程序的效率提供了方便。

2.1 求解器結(jié)構

求解器在耗時較多的分析部分進行并行處理，具體流程如圖1所示。當程序運行到插入節(jié)點力、形成荷載列陣和引入邊界條件、修改方程組模塊式時，進程0要向后層其他各進程廣播這些信息。當程序運行到計算矩陣預分配參數(shù)模塊時，各進程之間要保持相互通信。最后，各進程把各自計算出的結(jié)果返回到進程0，并保存結(jié)果，結(jié)束計算。并行分析部分包括預分配、形成總剛度矩陣及求解方程組，這些部分消耗的時間一般情況下占整個過程的90%以上，而且這一比例會隨著計算規(guī)模的增大而增大，用如上并行方式求解線彈性問題，理論上能取得很好的效果。

2.2 求解器調(diào)用

求解器提供的可選設置選項說明如下:

a.SetElement(intType，intmaxnode)與 SetElement(int Type，proc－ESElementStiffer)是提供接口給用戶接入自己的單元類型。其中參數(shù)Type是單元類型號，求解器默認類型為八節(jié)點六面體空間等參元;參數(shù)maxnode為該類型單元的節(jié)點數(shù);自定義指針函數(shù)類型 proc－ElementStiffer則為指向該單元類型的單剛求解程序。通常這2條設置應該設置在初始化函數(shù)SolverInitialize()之后。

b.SetKSPType(char*ktype)與 SetPCType(char*ptype)語句用于設置子空間迭代法的類型，和PETSc工具箱的命名一樣。子空間迭代法的類型 ktype可以是KSPCG，KSPGMRES，KSPTCQMR，KSPBCGS，KSPTFQMR，KSPCR，KSPLSOR，KSPBICG中的任何一種。預條件法的類型 ktype包括 PCJACOBI，PCBJACOBI，PCSOR，PCEISENSTAT。

圖1 求解器流程Fig.1 Flow chart of solver

c.SetTolerances(double abstol，double rtol，double dtol，int maxits)用于設置迭代收斂參數(shù)，前3個參數(shù)是基于－范數(shù) ( (∑)1/2)的剩余范數(shù)相對下降量rtol，剩余范數(shù)的絕對值atol及剩余中的相對增加量dtol。Maxits則為最大迭代次數(shù)。

d.SetTimeInfo()函數(shù)用于保存求解時間明細，包括各主要部分的時間消耗，以此來分析求解效率。SaveBinary()與SaveInfo()是將剛度矩陣和等效節(jié)點力列陣分別保存為文本形式和二進制形式，但只在單進程下有效。其中二進制保存的格式為PETSc規(guī)定的格式。使用求解器的類，即為依次調(diào)用類中的成員函數(shù)來實現(xiàn)有限元分析的過程。求解器的一個完整的使用示例如下:

SOLVER S;聲明一個求解器對象

S.Solver Initialize(argc，args);求解器初始化

if(S.rank=0)

S.Read Information(Char*path);由根進程讀入數(shù)據(jù)

Set Element(int Type，int maxnode);設置單元節(jié)點數(shù)，Type為單元類型編號

Set Element(int Type，proc－ESElement Stiffer);設置單剛求解程序

S.Preallocation();創(chuàng)建矩陣、向量結(jié)構、矩陣預分配

S.Add KERE();計算單剛、組集總剛和節(jié)點力列陣

Set KSP Type(char*ktype);設置求解迭代方法

Set PC Type(char*ptype);設置預條件方法

S.Set Bound Condition();根據(jù)邊界條件修正支配方程

Set Tolerances(double abstol，double rtol，double dtol，int maxits);設置收斂誤差

Save Binary();保存系數(shù)矩陣和右端頂，僅單進程運行時有效

腎功能指標水平：取上述血清3 ml，用肌氨酸氧化酶檢測血肌酐(Scr)法酶偶聯(lián)速率法檢測血尿素氮(BUN)；采用核素動態(tài)顯像法測定腎小球過濾率(GFR)。

Save ASCII();文本保存系數(shù)矩陣和右端項，僅單進程運行時有效

S.SOLVE();求解方程組

S.Get Answer();保存結(jié)果

Set Time Info();保存時間耗用明細

S.Finalize();結(jié)束

3 求解器的輸入輸出

求解器的輸入目前提供了2種形式:讀入文件和傳入網(wǎng)格信息結(jié)構體指針。用戶在調(diào)用Solver Initialize()函數(shù)進行初始化以后，即可在進程0里接受該結(jié)構體指針或執(zhí)行讀入文件函數(shù)。如以讀入文件的方式輸入，需要讀入節(jié)點信息、單元信息、材料信息、約束信息、節(jié)點力信息。

求解器的輸出包括3個部分:(a)方程構成。求解器能夠在一個進程運行時輸出方程的系數(shù)矩陣和右端項，默認的狀態(tài)是不輸出，用戶需要在求解方程之前將其設置成輸出。輸出文件的格式有文本文檔和二進制格式。文本文檔可供查閱，保存時間較長，存儲較大;二進制格式存儲小，速度快，可由PETSc程序讀入。(b)求解信息。求解網(wǎng)格信息的簡要說明，求解設置，以及可選的時間信息，供性能分析使用。(c)方程的解。包括方程的解和節(jié)點的位移矢量大小，在文件的開始則為位移的最大值及其節(jié)點號。

4 算例分析

筆者以計算某軟土地基上水閘在自重荷載下所產(chǎn)生的位移為算例，對本文所開發(fā)的并行求解器進行測試。水閘結(jié)構由3部分組成:閘身、樁基、軟土地基，材料參數(shù)見表1。本文僅作線彈性分析，故將地基和樁基密度置為零，閘身密度為2500kg/m3，以計算由水閘自身重力產(chǎn)生的位移。因水閘結(jié)構對稱分布，因此取其一半結(jié)構進行建模［16］計算(圖2)。從求解器求解規(guī)模、并行效率、對于病態(tài)方程組的迭代收斂情況來測試本文開發(fā)的并行求解器。求解環(huán)境為一臺四核I5、主頻為2.8 GHz的32位臺式機，可使用內(nèi)存為3 G。

表1 水閘材料參數(shù)Table 1 Material parameters of water lock

4.1 求解規(guī)模

按網(wǎng)格疏密將水閘結(jié)構離散為八節(jié)點六面體單元網(wǎng)格，表2為6組有限元模型的網(wǎng)格信息，圖2為網(wǎng)格單元數(shù)為51112的水閘網(wǎng)格圖。對每一組有限元模型都進行四進程并行求解，計算的時間如表2所示。

第6組模型在計算時由于每個進程的內(nèi)存使用都超過600 M，系統(tǒng)需要使用虛擬內(nèi)存，導致計算時間達到2 409.33 s。在同樣的計算環(huán)境下，單進程能計算節(jié)點數(shù)不超過40萬的網(wǎng)格，Ansys軟件能求解節(jié)點數(shù)不超過20萬的線彈性問題，可見本文求解器能在普通的多核計算機下處理更大規(guī)模的網(wǎng)格。

取節(jié)點數(shù)為58 622的網(wǎng)格模型比較本文求解器和Ansys軟件計算得到的節(jié)點位移，兩者的x，y，z方向最大位移及總位移間的誤差分別為0，0，0.007%，0.001%，最大位移值的誤差不超過1/10000(表3)。將節(jié)點位移導入軟件Tecplot中，繪制位移云圖。本文求解器得到的位移云圖與Ansys結(jié)果分布規(guī)律基本一致。求解器在提高收斂精度下能得到更精確的解，但因迭代次數(shù)增加而消耗更多的時間。

圖2 水閘結(jié)構及整體網(wǎng)格Fig.2 Structure of water locks and overall grid

表2 不同規(guī)模的網(wǎng)格信息Table 2 Grid information with different scales

表3 本文求解器和Ansys軟件計算的部分節(jié)點的結(jié)果對比Table 3 Comparison of node calculation results obtained by solver developed in this study and ANSYS software

不同規(guī)模網(wǎng)格計算得出的最大節(jié)點位移值如圖3所示，可看出隨著單元總數(shù)增加，節(jié)點最大位移值也逐漸增大，此現(xiàn)象符合有限單元法中位移近似解是從下方逼近精確解的理論，因而水閘結(jié)構的最大位移應該在0.067 m附近。

4.2 并行加速效果

為驗證并行加速效果，選取節(jié)點總數(shù)為178 007的網(wǎng)格數(shù)據(jù)，分別進行1～4個進程的并行計算。迭代方法為共軛梯度法(CG)［17］，預條件子為塊雅可比方法(bjacobi)。圖4為程序在單機多進程并行求解時的并行加速比、并行效率與進程數(shù)相關關系曲線。

程序的并行分析部分主要包括矩陣預分配參數(shù)計算、形成總剛度矩陣和迭代求解方程組。表4列出了程序主要并行部分的并行加速比和并行效率。由表4可知，求解器用于計算求解矩陣分配參數(shù)的時間較少，對于一個節(jié)點數(shù)為17萬的網(wǎng)格，計算時間不到1 s。計算單元剛度矩陣、組集總剛度矩陣時，由于并行程度較高，其加速效果比較理想。另外，求解的迭代次數(shù)隨著進程數(shù)目的增加而變大。

圖3 不同網(wǎng)格規(guī)模下最大節(jié)點位移Fig.3 Maximum displacement of node for different scales of grid

求解器由于其總剛度矩陣的計算采取了精確預分配空間的方法，只存非零元素，故程序運行時占用內(nèi)存較小。多個進程運行時，每個進程占用內(nèi)存會明顯減小。

4.3 對于病態(tài)方程組的迭代收斂

在有限元中，若結(jié)構材料參數(shù)相差過大會導致支配方程呈病態(tài)，即出現(xiàn)迭代次數(shù)過多或者迭代不收斂的問題。PETSc提供的迭代法的預條件子，能有效地解決或緩解這一問題。不考慮實際情況，將閘身的彈性模量作表5所示改變，即改變閘身彈性模量使之分別相差10倍、100倍至100萬倍。將這7組模型分別進行4個進程的并行迭代求解，得出的位移與Ansys軟件計算出的位移進行對比(表6)。由表6可得，采用預條件子技術后，求解程序的迭代方法對因材料參數(shù)相差過大(最大相差10萬倍)而產(chǎn)生的病態(tài)方程組具有良好的適應能力。

圖4 1～4個進程計算時的并行加速比與并行效率Fig.4 Parallel speedup ratio and efficiency for one to four processes

表4 并行分析各部分的并行加速比和并行效率Table 4 Parallel speedup ratio and efficiency for various parts

表5 閘身彈性模量變化Table 5 Lock body elastic modulus change Pa

表6 修改材料參數(shù)后求解器的計算結(jié)果Table 6 Calculated results of solver after material parameters were modified

5 結(jié) 語

基于PETSc的數(shù)據(jù)結(jié)構(矩陣Mat和向量Vec)和線性求解器(KSP)組件開發(fā)了有限元求解器，并以水閘在自重荷載下所產(chǎn)生的位移為算例，對該求解器進行了測試。測試結(jié)果表明，該求解器能在普通的內(nèi)存為3G的四核計算機上，計算高達80萬節(jié)點的線彈性有限元問題，并行加速效果也較為明顯(4個進程計算時加速比能達到3.24)。另外由于使用了PETSc提供的較為成熟的預條件子技術，求解器對病態(tài)方程組也有良好的適應能力，為求解由于材料參數(shù)相差10萬倍而產(chǎn)生的病態(tài)方程組提供了新方法。由于現(xiàn)階段各方面的限制，本文開發(fā)的求解器僅僅使用了PETSc工具箱的很小一部分功能，程序也有很多地方需要改進?；赑ETSc的有限元高性能求解方法的研究和開發(fā)還可以進行以下幾方面的深入研究:

a.對程序進一步完善擴展，能夠?qū)€性有限元問題作進一步的求解，如求出應力、應變。

b.將PETSc中的非線性求解器SNES和時間步進器TS對應到非線性有限元、時間依賴的動力學問題的求解上。

c.本文所開發(fā)的求解器是一個沒有界面的程序，使用起來不方便，若能將其與圖像處理軟件結(jié)合起來或自己形成界面將會使該求解器更接近實際應用。

［1］范穎，王磊，章青.多尺度有限元法及其應用研究進展［J］.水利水電科技進展，2012，32(3):90-94.(FAN Ying，WANG Lei，ZHANG Qing.Research progress and application of multiscale finite element method［J］.Advances in Science and Technology of Water Resources，2012，32(3):90-94.(in Chinese))

［2］張汝清.概說并行計算結(jié)構力學［J］.計算結(jié)構力學及其應用，1995，12(4):447-484.(ZHANG Ruqing.Introduction to parallel computational structural mechanics［J］.Computational Structural Mechanics and Applications，1995，12(4):447-484.(in Chinese))

［3］張健飛.采用PETSc的有限元并行計算實現(xiàn)與優(yōu)化［J］.計算機工程與應用，2010，46(10):57-59.(ZHANG Jianfei.Parallel implementation and tuning of finite element computing with PETSc［J］.Computer Engineering and Applications，2010，46(10):57-59.(in Chinese))

［4］姜弘道.水利高性能計算的進展［J］.水利水電科技進展，2006，26(2):70-76.(JIANG Hongdao.Progress in highperformance computing in water science［J］.Advances in Sciences and Technology of Water Resources，2006，26(2):70-76.(in Chinese))

［5］NOOR A K，VOIGT S J.Hypermatrix scheme for finite element system on CDC STAR-100 computer［J］.Computers and Structures，1975，5(5):287-296.

［6］ DIMITRIS Z.Parallel processing techniques for FE analysis:stiffness，loads and stress evaluation［J］.Computers and Structures，1988，28(2):247-260.

［7］ BITZARAKIS S，PAPADRAKAKIS M.Parallel solution techniques in computational structural mechanics［J］.Computer Methods in Applied Mechanics and Engineering，1997，148:75-104.

［8］CHIEN L S，SUNCT.Parallel processing techniques for finite element analysis of nonliear large truss structures［J］.Computers and Structures，1989，31(6):1023-1029.

［9］李麗君，金先龍，李淵印，等.一種新型并行化有限元結(jié)構模態(tài)分析集成系統(tǒng)［J］.計算力學學報，2004，21(5):546-550.(LI Lijun，JIN Xianlong，LI Yuanyin，et al.A new parallel integration system of finite element model analysis［J］.Chinese Journal of Computational Mechanics，2004，21(5):546-550.(in Chinese))

［10］張林波，遲學斌，莫則堯，等.并行計算導論［M］.北京:清華大學出版社，2006:1-487.

［11］Michael JQ.MPI與OpenMP并行程序設計［M］.北京:清華大學出版社，2004:34-40.

［12］都志輝.高性能計算并行編程技術MPI并行程序設計［M］.北京:清華大學出版社，2001:1-336.

［13］KRYLOV A N.On the numerical solution of the equation by which the frequency of small oscillations is determined in technical problems［J］.Izv Akad NaukAzerb Ser Fiz-Tekh Mat Nauk，1931，4:491-539.

［14］李曉梅，吳建平.Krylov子空間方法及其并行計算［J］.計算機科學，2005，32(1):19-20.(LI Xiaomei，WU Jianping.Krylovsubspace methods and parallel computation［J］.Computer Science，2005，32(1):19-20.(in Chinese))

［15］付成華，周洪波.面向?qū)ο蟮蔫F路路基加固設計系統(tǒng)集成與實現(xiàn)［J］.水利水電科技進展，2012，32(4):70-73.(FU Chenghua，ZHOU Hongbo.Development and implementation of object-oriented railway subgrade strengthening design system［J］.Advances in Science and Technology of Water Resources，2012，32(4):70-73.(in Chinese))

［16］林江，母德偉，余葵，等.萬州長江公路大橋防撞裝置結(jié)構穩(wěn)定性及通航凈空尺度［J］.水利水電科技進展，2013，33(2):59-62.(LIN Jiang，MU Dewei，YU Kui，et al.Structural stability and navigable clearance of anti-collision device of Wanzhou Yangtze River Highway Bridge［J］.Advances in Science and Technology of Water Resources，2013，33(2):59-62.(in Chinese))

［17］HESTENSEM R，STIEFEL E.Methods of conjugate gradients for solving linear systerms［J］.Tournal of Reaserch of the Natural Bureau of Standards，1952，49(6):409-436.