馮春生,　舒　適,　梁文濤

(湘潭大學(xué) 數(shù)學(xué)與計算科學(xué)學(xué)院，科學(xué)工程計算與數(shù)值仿真湖南省重點實驗室，湖南 湘潭 411105)

彈性力學(xué)問題是固體力學(xué)中的一種基本模型，它廣泛應(yīng)用于建筑、機械、化工、航天等領(lǐng)域. 有限元方法是求解彈性力學(xué)問題的最常用的離散化方法之一.由于有限元方程系數(shù)矩陣的條件數(shù)會隨著網(wǎng)格尺寸的減小而迅速增大，用經(jīng)典迭代法求解時收斂較慢.因此，研究求解大規(guī)模彈性問題有限元離散系統(tǒng)的快速算法非常必要.

非重疊DDM 是為大規(guī)模偏微分方程離散系統(tǒng)構(gòu)造有效并行預(yù)條件子的一類非常重要的方法[1-3].標(biāo)量橢圓方程方面研究工作已比較成熟，[1]針對H(grad)型橢圓問題設(shè)計了非重疊DDM 預(yù)條件子.非重疊DDM預(yù)條件子經(jīng)過近幾十年的發(fā)展，先后形成了如BPS型[1]、Neumann型[2]和FETI型[3]等構(gòu)造預(yù)條件子的模式.最近，[4]為標(biāo)量橢圓問題線性元方程設(shè)計了一種基于簡單粗空間的非重疊DDM 高效預(yù)條件子.多層網(wǎng)格法[5-7]是高效求解規(guī)模偏微分方程離散系統(tǒng)的另一類重要方法.代數(shù)多層網(wǎng)格(AMG)法誕生于20世紀(jì)80年代初[8]，它是求解大規(guī)模偏微分方程離散化系統(tǒng)的最為有效的迭代方法之一.20世紀(jì)90年代以來，受實際應(yīng)用問題的驅(qū)動，AMG 法得以迅速發(fā)展[9-10].目前常用的AMG 法包括經(jīng)典AMG法(簡稱C-AMG法)[11]，基于能量極小基的AMG法[12]，自適應(yīng)AMG法[13]和GAMG 法等，其中GAMG法是指針對不同特性的問題，通過利用其中部分容易獲得的幾何或分析信息而設(shè)計的依賴于問題的AMG 法，如AMGe法[14]和基于光滑聚集的AMG法[15]等.GAMG法由于具有更高的運算效率，已成為當(dāng)前國際AMG 法研究領(lǐng)域中的熱點.

1　線彈性模型問題與線性有限元方程

考慮如下三維線彈性問題：

(1)

(2)

其中

(3)

(4)

變分問題式(2)的線性有限元近似變分問題為：求uh∈Vh(Ω)，使得

a(uh,vh)=(f,vh),?vh∈Vh(Ω),

(5)

則式(5)的代數(shù)形式為

Au=f.

(6)

下面我們利用代數(shù)多層網(wǎng)格法與非重疊DDM預(yù)條件子[4]， 為求解式(6)的PCG法構(gòu)造一種并行預(yù)條件子.

2　一種并行預(yù)條件子

為了給出方程式(6)中系數(shù)矩陣A的一種基于非重疊DDM的預(yù)條件行為B的構(gòu)造算法，只需給出求向量w=Bg的算法， 其中g(shù)為任給的已知向量，其計算公式為

(7)

其中，關(guān)于求解式(7)右端中的各個向量的算法可描述如算法1 (這里設(shè)wd,wW,wl,g分別為以wd,wW,wl和g的分量為自由度的線性有限元函數(shù))：

算法1非重疊DDM的預(yù)條件子

(1) 求wd∈Vd(Ω)，滿足如下變分子問題：

a(wd,vd)=(g,vd),?vd∈Vd(Ω).

(8)

(2) 求wW∈VW(Ω)，滿足如下變分子問題：

a(wW,vW)=(g,vW),?vW∈VW(Ω).

(9)

(3) 對內(nèi)交界面循環(huán)，設(shè)當(dāng)前交界面號為l(它對應(yīng)內(nèi)交界面Γij)，求wl∈Vl(Ωij)，滿足如下變分子問題：

a(wl,vl)=(g,vl),?vl∈Vl(Ωij).

(10)

第一類子系統(tǒng)式(8)：由于其自由度較少，使用經(jīng)典迭代法CG 或C-AMG都可求解；第二類子系統(tǒng)式(9)：由于其系數(shù)矩陣性態(tài)較好，使用CG求解即可；第三類子系統(tǒng)式(10)：每個內(nèi)界面對應(yīng)一個，此類子問題求解性態(tài)的好壞，將直接影響整個預(yù)條件子系統(tǒng)求解的效率.針對第三類子系統(tǒng)式(10)的求解，研究發(fā)現(xiàn)隨著子系統(tǒng)規(guī)模的增大，使用經(jīng)典迭代法如CG求解時，迭代次數(shù)隨子問題規(guī)模的擴大增長很快，而使用C-AMG求解該類子系統(tǒng)時，迭代次數(shù)不穩(wěn)定，且求解效率較低.因此，需要設(shè)計新的高效AMG 法.特別地，通過改變C-AMG 中的粗化策略和提升算子的構(gòu)造方法，我們?yōu)榈谌愖酉到y(tǒng)設(shè)計了一種新的AMG (簡記為AMG-T)法.

3　并行實現(xiàn)

記Ad:Vd(Ω)→Vd(Ω)，AW:VW(Ω)→VW(Ω) 和Al：Vl(Ωij)→Vl(Ωij)為三個對稱正定算子；記Qd，QW和Ql分別表示Vh(Ω) 到相應(yīng)子空間Vd(Ω)，VW(Ω)和Vl(Ωij)的L2投影算子.記變分子問題式(8)、(9)、(10)對應(yīng)的離散系統(tǒng)分別為:

Adud=Qdg，

(11)

AWuW=QWg，

(12)

Alul=Qlg.

(13)

這里，在不與算子Ad,AW,Al產(chǎn)生混淆的情況下，記其對應(yīng)的矩陣為Ad,AW,Al(l=1,…,M).

接下來，針對AW規(guī)模較大，且子系統(tǒng)式(13) 的求解占主要工作量的情形簡要給出MPI + OpenMP二級并行實現(xiàn)策略.關(guān)于第一類子系統(tǒng)式(11)：由于其自由度較少, 通?？梢栽诘?#進程上采用BoomerAMG[16]進行OpenMP 并行求解.如果直接調(diào)用BoomerAMG,求解效率不夠理想，可以采用我們新設(shè)計的一種解法器AMG-T(與前者不同之處在于粗化策略考慮方向性).由表1可知, 當(dāng)md較大時，AMG-T的迭代次數(shù)比較穩(wěn)定.

關(guān)于第二類子系統(tǒng)式(12)：由于其系數(shù)矩陣性態(tài)較好，調(diào)用OpenMP 并行CG在第1#進程上求解即可.由表2可知，并行CG的迭代次數(shù)與網(wǎng)格規(guī)模和進程數(shù)均無關(guān).

關(guān)于第三類子系統(tǒng)式(13)，注意以下兩個事實：

表1　求解第一類子系統(tǒng)(11)的平均迭代次數(shù)

表2　并行 CG 求解第二類子系統(tǒng)(12)的迭代次數(shù)

事實2該類子系統(tǒng)的求解可同步進行，且可自由組合聚集到不同進程并行計算.

因此，我們將第三類子系統(tǒng)平衡分配到余下進程并行計算.進一步，在各進程上將分配到該進程上的第三類子系統(tǒng)再平均分配到各線程上，采用OpenMP并行求解.

(14)

(15)

|Ni-Nj|≤1, |Ni, l-Ni, m|≤1, ?i≠j,?l≠m.

(16)

m=1(1)Nmyid, l,分別表示第myid號進程上的第l個線程上的第m個轉(zhuǎn)換算子對應(yīng)的矩陣和相應(yīng)子系統(tǒng)的系數(shù)矩陣.

算法2MPI+OpenMP并行非重疊DDM的預(yù)條件子

Step 1初始化并行向量z.

Step 2將并行向量r轉(zhuǎn)為串行向量， 不妨仍記為r.

Step 3若myid=0，則zmyid:= (PdTBdPd)r；

若myid=1， 則zmyid:= (PWTBWPW)r；

若myid≠0 且myid≠1并行做以下循環(huán)

form=1,…,nts

Step 3.1求Ft(l,m):=Pmyid(l,m)r.

Step 3.2求yt(l,m)=Bmyid(l,m)Ft(l,m)

Step 3.3求zmyid:=zmyid+(Pmyid(l,m))Tyt(l,m).

end for

注1若子問題規(guī)模 (存儲規(guī)模) 能由一個計算結(jié)點容納時，子問題求解調(diào)用HYPRE[16]中串行AMG 或AMG-T；否則，當(dāng)子問題規(guī)模超過一個計算結(jié)點的容納能力時，子問題求解調(diào)用HYPRE中并行AMG或并行AMG-T求解.本文數(shù)值實驗部分的子問題求解調(diào)用HYPRE中的串行AMG 或AMG-T.

4　數(shù)值實驗

在線彈性模型問題式(1)中， 取單位立方體Ω:=(0,1)×(0,1)×(0,1)；E=3 000,v=

0.3；f=(f1,f2,f3)T， 其中

f1(x,y,z)=(λ+4μ)π2sin(πx)sin(πy)sin(πz)-(λ+μ)(3x2-4x3)πcos(πy)sin(πz)-

(λ+μ)π2y3(1-y)cos(πx)cos(πz),

f2(x,y,z)=-(λ+μ)π2cos(πx)cos(πy)sin(πz)-(λ+μ)(3y2-4y3)πsin(πx)cos(πz)+

[(12x2-6x)μ+π2x3(1-x)(λ+3μ)]sin(πy)sin(πz),

f3(x,y,z)=(λ+μ)π2cos(πx)sin(πy)cos(πz)-(λ+μ)π2x3(1-x)cos(πx)cos(πz)+

[(12y2-6y)μ+π2y3(1-y)(λ+3μ)]sin(πx)sin(πz).

數(shù)值實驗硬件環(huán)境為某2路4核工作站(CPU： Intel Xeon W5590*2， Mem：24 GB)； 軟件環(huán)境為Linux Redhat 5.3，gcc v4.1.2，gfortran v4.1.2，mpich2-1.3.2p1，hypre-2.7.0b.由于該工作站總共只有8 個核,因此我們只針對進程數(shù)和線程數(shù)不大于8 的情形進行數(shù)值實驗，這樣可以保證進程或線程與核一一對應(yīng).實驗中取PCG迭代控制精度為10-6.

算法可擴展性表4給出了不同網(wǎng)格剖分Τmd,k和不同進程數(shù)np下PCG算法的外迭代次數(shù).

由表4可知，當(dāng)d/h不變時，條件數(shù)κ(BA) 為常量，表明該預(yù)條件子的算法可擴展性達到最優(yōu)；當(dāng)d不變，h變小時，條件數(shù)κ(BA) 僅弱依賴于d/h=k.

第三類子系統(tǒng)的并行可擴展性表5給出了相應(yīng)的CPU時間與并行效率.

由表5可知，當(dāng)進程數(shù)np=8 時加速比達到95%以上而未達到100%.可能的原因為：負載不夠平衡.可以預(yù)測，在更多進程參與計算的情況下，該類子問題的求解加速性能不差.若進程數(shù)足夠多，該類子問題的求解時間將大幅下降.

表4　不同Τmd,k與np下PCG的外迭代次數(shù)

表5　不同Τmd,k與np下求解第三類子系統(tǒng)的CPU時間與并行效率

表6　在不同并行模式下求解的時間比值

實驗結(jié)果表明在曙光5000上求解該DDM預(yù)條件系統(tǒng)時，采用MPI+OpenMP 模式比MPI模式計算效率更高.

感謝湘潭大學(xué)王俊仙博士在并行DDM 算法設(shè)計方面提供的良好建議.

[1]DRYJA M, SMITH F,WIDLUND O.Schwarz analysis of iterative substructuring algorithms for elliptic problems in three dimensions[J].SIAM J Numer Anal, 1994,31(6):1662-1694.

[2]DRYJA M,WIDLUND O. Schwarz methods of neumann type for three-dimensional elliptic finite element problems[J].Comm Pure Appl Math, 1995, 48:121-155.

[3]HU Q,SHI Z,YU D.Efficient solvers for saddle point problems arising from domain decompositions with lagrange multipliers[J].SIAM J Numer Anal, 2004,42(3):905-933.

[4]HU Q,SHU S,WANG J.Nonoverlapping domain decomposition methods with simplified coarse spaces for solving three-dimensional elliptic problems[J].Math Comp, 2010, 79(272): 2059-2078.

[5]FENG C,SHU S, XU J,et al.Numerical study of geometric multigrid methods on CPU-GPU heterogeneous computers[J].AAMM, 2014, 6(1):1-23.

[6]FENG C,ZHANG S.Optimal solver for morley element discretization of biharmonic equation on shape-regular grids[J].JCM, 2016, 34(2):159-173.

[7]譚林, 江軍, 舒適, 一種三角形網(wǎng)格下保對稱有限體元方法的預(yù)條件技術(shù)[J].湘潭大學(xué)自然科學(xué)學(xué)報, 2006,28(1):1-6.

[8]BRANDT A, MCCORMICK S,RUGE J.Algebraic multigrid (AMG) for automatic multigrid solution with application to geodetic computations[P].Institute for Computational Studies, POB 1852, Fort Collins, Colorado, 1982.

[9]STERCK H D,YANG U M,HEYS J J.Reducing complextity in parallel algebraic multigrid preconditioners[J].SIAM J Mat Anal Appl,2006,27(4):1019-1039.

[10]BIENZ A,FALGOUT R D,GROPP W,et al.Reducing parallel communication in algebraic multigrid through sparsification[J].SIAM J Sci Comp, 2016, 38(5): S332-S357.

[11]RUGE J W, STüBEN K. Algebraic multigrid[M]//Multigrid Methods.Philadelphia:SIAM,1987:73-130.

[12]XU J,ZIKATANOV L.On the energy minimizing base for algebraic multigrid methods[J].Computing and Visualization in Science, 2004, 7:121-127.

[13]BREZINA M,FALGOUT R D,MACLACHLAN S,et al.Addaptive smoothed aggregation(αSA)[J].SIAM Review, 2005, 47(2):317-346.

[14]CHARTIER T,FALGOUT R D,HENSON V E,et al.Spectral AMGe(ρAMGe)[J].SIAM J Sci Comp,2003,20(1):1-20.

[15]VANEK P,MANDEL J,BREZINA M.Algebraic multigrid by smoothedaggregation for second and fourth order elliptic problems[J].Computing, 1996, 56:179-196.

[16]趙永華, 遲學(xué)斌.基于SMP集群的MPI+OpenMP混合編程模型及有效實現(xiàn)[J].微電子學(xué)與計算機,2005,22(10):7-11.

[17]FALGOUT R D,BAKER A,CHOW E,et al.HYPRE: High performance preconditioner[OL].https://computation.llnl.gov/casc/hypre/software.html.