結(jié)構(gòu)模態(tài)多級分層并行計算方法

喻高遠(yuǎn), 樓云鋒, 李俊杰, 金先龍

(1. 上海交通大學(xué) 機(jī)械與動力工程學(xué)院,上海 200240;2. 上海宇航系統(tǒng)工程研究所,上海 201108)

隨著科學(xué)研究和工程技術(shù)的不斷發(fā)展,出現(xiàn)了諸如高速動車組、3 000 m超深鉆機(jī)、大飛機(jī)、大功率海洋壓縮機(jī)等復(fù)雜裝備系統(tǒng)及跨江隧道、摩天大廈等大型和超大型工程系統(tǒng)[1]。這些系統(tǒng)安全性和穩(wěn)定性等方面的苛刻要求對其動力學(xué)系統(tǒng)性能數(shù)值模擬提出了嚴(yán)峻的挑戰(zhàn),其整體有限元系統(tǒng)建模需要考慮精細(xì)建模以反映局部響應(yīng),致使其有限元系統(tǒng)自由度規(guī)模巨大,求解復(fù)雜[2-4]。在這些復(fù)雜動力學(xué)系統(tǒng)的計算過程中,模態(tài)分析是其最耗費時間的環(huán)節(jié),也是其余計算環(huán)節(jié)的基礎(chǔ),需借助大規(guī)模有限元模型進(jìn)行高精度模態(tài)計算。這些問題在傳統(tǒng)的串行機(jī)上無法得到滿意的解答[5-6]。隨著超級計算機(jī)的快速發(fā)展和大量先進(jìn)算法的涌現(xiàn),利用其研究和開發(fā)相應(yīng)的大規(guī)模模態(tài)并行算法則為這類問題的解決提供了切實可行的方法。

模態(tài)分析的數(shù)學(xué)實質(zhì)可以歸結(jié)為大型稀疏矩陣的廣義特征值問題[7-8],其求解方法主要有Lanczos算法[9]、Arnoldi算法[10]、Krylov-Schur算法[11]等。各國學(xué)者在Lanczos算法、Arnoldi算法、Krylov-Schur算法等的基礎(chǔ)上進(jìn)行了一系列的模態(tài)并行算法研究,主要集中在:從模態(tài)有限元分析中耗時最多的線性方程組求解出發(fā)尋求高效率求解線性方程組的并行計算方法和從有限元問題自身的并行性出發(fā)形成的模態(tài)綜合法兩大類并行計算方法[12]。直接法和迭代法是模態(tài)并行求解的兩種基本算法。直接法通過排序、三角分解和回代計算能夠在預(yù)期內(nèi)得到方程的解。然而,隨著有限元計算規(guī)模的增加,其所需的內(nèi)存空間和計算量也會顯著增加。迭代法通過多次迭代對求解結(jié)果進(jìn)行改善以達(dá)到收斂容差范圍,求解過程中所需內(nèi)存小,容易實現(xiàn)并行化。但迭代法無法保證收斂的合理時間,且對于條件數(shù)很大的病態(tài)問題也可能是不收斂的。在模態(tài)綜合并行計算方面,各國學(xué)者將直接法和迭代法混合并從有限元問題自身的并行出發(fā),將復(fù)雜模態(tài)求解問題拆分為一個個較小的子模態(tài)問題進(jìn)行并行處理,子問題模態(tài)求解過程中方程求解采用直接法,通過2次模態(tài)坐標(biāo)變換得到縮聚后的系統(tǒng)級模態(tài)廣義特征方程以降低其整體規(guī)模,然后縮聚后的系統(tǒng)級模態(tài)廣義特征方程中涉及線性方程組的求解采用迭代法,因而能夠很好地利用直接法和迭代法的優(yōu)點提高并行效率,在結(jié)構(gòu)模態(tài)并行求解領(lǐng)域得到了廣泛的應(yīng)用。Rong等[13]基于模態(tài)綜合法和傳遞矩陣法完成了大規(guī)模特征值的混合并行求解設(shè)計,并將其應(yīng)用于某真空容器的大規(guī)模模態(tài)并行求解。Heng等[14]基于模態(tài)綜合法利用共享式存儲并行計算機(jī)完成了某懸臂梁的大規(guī)模模態(tài)并行求解。Aoyama等[15]通過改進(jìn)模態(tài)綜合法各子系統(tǒng)界面的耦合方式完成了某矩形板的大規(guī)模模態(tài)并行求解。P?rík等[16]采用模態(tài)綜合法基于OpenMP完成了模態(tài)并行計算求解體系設(shè)計,并利用其完成了某軸系的模態(tài)并行求解。然而對于大規(guī)模問題采用模態(tài)綜合法進(jìn)行并行求解時,隨著子區(qū)域數(shù)目的增多,縮聚后所得系統(tǒng)級模態(tài)廣義特征方程的規(guī)模和條件數(shù)也隨之急劇增加,從而給求解帶來了困難。另外,在迭代求解系統(tǒng)級模態(tài)廣義特征方程時,由于所有的進(jìn)程都需要參與全局通信,進(jìn)程間通信和同步開銷的增加也會極大地降低并行效率。

在硬件方面,科學(xué)和工程計算中使用最廣的并行計算機(jī)為分布式存儲并行計算機(jī),其體系架構(gòu)主要包括純CPU的同構(gòu)型超級計算機(jī)[17]和異構(gòu)型超級計算機(jī)[18]。異構(gòu)型超級計算機(jī)以其高性價比、強(qiáng)計算能力等優(yōu)勢,已成為目前分布式存儲并行計算機(jī)的主流架構(gòu)。對于異構(gòu)型分布式存儲并行計算機(jī)來說,最為重要的部分就是集群環(huán)境下的不同存儲機(jī)制、處理器之間的相互通信與協(xié)作以及各級硬件體系結(jié)構(gòu)的負(fù)載均衡,這也是影響并行效率很重要的因素。因此,利用異構(gòu)型分布式存儲并行計算機(jī)提高并行效率的關(guān)鍵在于處理好計算任務(wù)與異構(gòu)集群硬件拓?fù)潴w系結(jié)構(gòu)適配時的負(fù)載均衡、大規(guī)模數(shù)據(jù)的存儲以及處理器間的相互通信和協(xié)作問題。

為解決采用模態(tài)綜合法利用大量異構(gòu)核組求解大規(guī)模模態(tài)問題并行求解效率低的問題,本文在稀疏存儲和傳統(tǒng)并行模態(tài)綜合法的基礎(chǔ)上提出了一種結(jié)構(gòu)模態(tài)多級分層并行計算方法。采用兩級分區(qū)4次變換,降低系統(tǒng)整體縮減后的廣義特征方程規(guī)模。

1 模態(tài)綜合法

在傳統(tǒng)的并行模態(tài)綜合法中[19-20],有限元網(wǎng)格僅被剖分1次,然后在生成的m個子區(qū)域上按照先內(nèi)部自由度I集后邊界自由度B集的編號原則,組集每個子區(qū)域系統(tǒng)的質(zhì)量矩陣Msub和剛度矩陣Ksub。進(jìn)行第1次坐標(biāo)變換時,各子區(qū)域系統(tǒng)物理坐標(biāo)與模態(tài)坐標(biāo)之間的關(guān)系可以表示為

(1)

式中:UI,UB,{pk},{pB}分別為內(nèi)部節(jié)點和邊界節(jié)點對應(yīng)的位移和模態(tài)坐標(biāo);[Φ]為子區(qū)域的模態(tài)轉(zhuǎn)化矩陣;[ΦIk]為內(nèi)部節(jié)點主模態(tài)矩陣;[ΦIB]為約束模態(tài)。

利用各子區(qū)域系統(tǒng)之間的位移協(xié)調(diào)條件,進(jìn)行第2次坐標(biāo)變換,如式(2)所示

(2)

K*φ*=λ*M*φ*

(3)

縮減剛度矩陣[K*]和縮減質(zhì)量矩陣[M*]為

(4)

采用并行Krylov子空間迭代法求解式(3)即可獲得整體系統(tǒng)所需要的低階模態(tài)頻率和模態(tài)振型,將其代入式(1)即可得各子區(qū)域的模態(tài)振型。

2 分層并行計算方法

分層并行計算方法通過兩級分區(qū)4次變換在分布式數(shù)據(jù)稀疏存儲和傳統(tǒng)并行模態(tài)綜合法的基礎(chǔ)上通過并行任務(wù)映射既改善了不同層級的負(fù)載均衡,又實現(xiàn)了通信分離有效提高了通信效率。此外,它還進(jìn)一步降低了界面方程的規(guī)模,加快了其迭代收斂速度。

2.1 并行任務(wù)映射

通過將計算任務(wù)映射到異構(gòu)眾核超級計算機(jī)的不同硬件層,以實現(xiàn)不同層級的負(fù)載均衡和通信的有效分割。在傳統(tǒng)并行模態(tài)綜合法的基礎(chǔ)上,考慮“神威太湖之光”異構(gòu)眾核超級計算機(jī)的硬件體系架構(gòu),完成結(jié)構(gòu)大規(guī)模模態(tài)分層并行計算任務(wù)映射,如圖1所示。

圖1 結(jié)構(gòu)模態(tài)多級分層并行模態(tài)綜合法計算任務(wù)映射Fig.1 Task mapping of multilevel hierarchical parallel modal synthesis algorithm

在進(jìn)行任務(wù)映射時,將第1級網(wǎng)格子區(qū)域按照節(jié)點順序進(jìn)行映射,第2級網(wǎng)格子區(qū)域按照節(jié)點內(nèi)異構(gòu)群組進(jìn)行映射,各異構(gòu)群組內(nèi)浮點運算按照計算核心進(jìn)行映射。

2.2 兩級分區(qū)

在兩級分區(qū)中,首先采用開源并行有限元分區(qū)程序ParMETIS將有限元網(wǎng)格剖分為p個1級子區(qū)域,然后每個1級子區(qū)域再進(jìn)一步按照異構(gòu)眾核體系架構(gòu)特點剖分為m個2級子區(qū)域。為與“神威太湖之光”異構(gòu)眾核分布式存儲超級計算機(jī)相適配,p應(yīng)為并行計算每次啟動節(jié)點機(jī)總數(shù),m應(yīng)為單個節(jié)點機(jī)內(nèi)異構(gòu)核組的數(shù)量,為4。如圖2所示,若p為2,則有限元網(wǎng)格先被剖分為2個1級子區(qū)域,然后每個1級子區(qū)域又被剖分為4個2級子區(qū)域。

圖2 2級分區(qū)Fig.2 Two-level partitioning

2.3 4次變換

4次變換通過先后在2級子區(qū)域和1級子區(qū)域上應(yīng)用式(2)和式(4)實現(xiàn)變換過程,如圖3所示。首先,形成每個2級子區(qū)域系統(tǒng)的質(zhì)量矩陣和剛度矩陣,并通過第1次坐標(biāo)變換構(gòu)成僅包含內(nèi)部自由度的廣義特征方程,通過子區(qū)域凝聚得到各子區(qū)域的等效剛度矩陣和等效質(zhì)量矩陣;然后,通過組集隸屬于相同1級子區(qū)域內(nèi)的所有2級子區(qū)域等效剛度矩陣和等效質(zhì)量矩陣,并通過第2次坐標(biāo)變換得到僅含有獨立坐標(biāo)的1級子區(qū)域廣義特征方程;最后,通過凝聚、組集和坐標(biāo)變換得到僅含有獨立坐標(biāo)的系統(tǒng)整體的廣義特征方程。

圖3 4次變換Fig.3 Four transformations

為改善計算過程中的內(nèi)存訪問速率,考慮到申威異構(gòu)眾核處理器各核組配置有8 G私有內(nèi)存且均可獨立訪問,在進(jìn)行計算時,各子區(qū)域的數(shù)據(jù)信息均通過多文件流存儲在相應(yīng)的節(jié)點機(jī)內(nèi)各異構(gòu)核組上。此外,相對于傳統(tǒng)并行模態(tài)綜合法,分層并行計算方法通過4次變換2次凝聚進(jìn)一步降低了系統(tǒng)整體廣義特征方程的規(guī)模,加快了其迭代收斂速度。

2.4 三層并行計算機(jī)制

考慮異構(gòu)眾核分布式存儲并行計算機(jī)的通信特點,控制通信協(xié)作以提高通信效率的關(guān)鍵在于實現(xiàn)節(jié)點內(nèi)通信和節(jié)點間通信分離、節(jié)點內(nèi)異構(gòu)核組內(nèi)通信和節(jié)點內(nèi)異構(gòu)核組間通信分離以及減少系統(tǒng)整體求解進(jìn)程間的通信和同步開銷。即要根據(jù)異構(gòu)眾核分布式存儲并行計算機(jī)的通信結(jié)構(gòu)將大量局部通信限制在節(jié)點內(nèi),并最大限度減少不同節(jié)點間的全局通信。而分層并行計算方法剛好滿足這些條件,如圖4所示,它是在兩級分區(qū)4次變換策略的基礎(chǔ)上實現(xiàn)了計算過程的三層并行。

圖4 三層并行機(jī)制Fig.4 Scheme of three-layer parallelization

為節(jié)省內(nèi)存空間和減少計算量,各2級子區(qū)域局部總體剛度矩陣和總體質(zhì)量矩陣均采用列壓縮稀疏存儲格式進(jìn)行存儲。

第2層并行,各節(jié)點機(jī)內(nèi)分別同時進(jìn)行相應(yīng)1級子區(qū)域的組集、并行凝聚和回代求解,通信存在于同一節(jié)點內(nèi)的不同異構(gòu)核組之間及各異構(gòu)核組內(nèi)運算控制核心和計算核心之間。計算過程中的數(shù)據(jù)采用列壓縮稀疏存儲技術(shù)進(jìn)行分布式存儲,在進(jìn)行計算時,由各節(jié)點機(jī)內(nèi)的0號異構(gòu)核組負(fù)責(zé)相應(yīng)1級子區(qū)域的組集、數(shù)據(jù)分發(fā)、并行凝聚結(jié)果的匯總和內(nèi)部回代求解。各節(jié)點機(jī)內(nèi)對應(yīng)1級子區(qū)域的并行計算過程則由同一節(jié)點機(jī)內(nèi)所有的異構(gòu)核組同時參與,1級子區(qū)域?qū)?yīng)的廣義特征方程(KII,MII)采用并行Lanczos方法進(jìn)行求解,求解時涉及的線性方程組采用并行預(yù)處理Cholesky分解。

第3層并行,利用并行Krylov子空間迭代法求解系統(tǒng)整體的縮減廣義特征方程,每個節(jié)點機(jī)僅有1個核組-0號異構(gòu)核組參與求解和通訊,如圖5所示。在進(jìn)行求解時,系統(tǒng)整體的縮減剛度矩陣[K*]和縮減質(zhì)量矩陣[M*]仍分布式存儲在各節(jié)點機(jī)0號異構(gòu)核組對應(yīng)的內(nèi)存空間當(dāng)中,中間計算結(jié)果也以列壓縮格式進(jìn)行分布式存儲。大量的局部通信存在于各異構(gòu)核組內(nèi)運算控制核心和計算核心陣列之間,只有少量的點積操作和整體迭代誤差的計算需要全局通信。因而它能夠有效減少通信量,提高并行計算效率。

圖5 并行Krylov算法及其任務(wù)映射FFig.5 Parallel Krylov algorithm and task mapping

3 算例驗證

為驗證所提算法的正確性和有效性,采用“神威太湖之光”超級計算機(jī)進(jìn)行測試,測試時每個節(jié)點啟動4個異構(gòu)核組。采用所提算法和傳統(tǒng)并行模態(tài)綜合法,進(jìn)行某超深鉆機(jī)盤鼓式制動器轉(zhuǎn)子盤模態(tài)分析,其有限元網(wǎng)格模型如圖6所示,該模型具有2 896 781個實體單元、自由度規(guī)模為6 176 367、剛度矩陣非零元個數(shù)為501 507 352、平均帶寬為260。固定約束其內(nèi)表面8個螺栓孔位置,計算結(jié)構(gòu)的前20階固有頻率,并與經(jīng)典模態(tài)求解算法-Lanczos算法[9,19]的求解結(jié)果進(jìn)行對比。測試所得的前20階固有頻率的最大相對誤差按照式(5)計算

圖6 某轉(zhuǎn)子盤有限元系統(tǒng)Fig.6 Finite element of rotor system

(5)

經(jīng)計算后可知,采用本文算法與傳統(tǒng)并行模態(tài)綜合法所得模態(tài)頻率的計算結(jié)果與經(jīng)典Lanczos算法的相對誤差均小于0.7%,各階振型保持一致,有效驗證了本文計算結(jié)果和傳統(tǒng)并行模態(tài)綜合法計算結(jié)果的正確性。

通過啟動相應(yīng)數(shù)目的節(jié)點機(jī)測試本文分層并行計算方法和傳統(tǒng)并行模態(tài)綜合法的性能。測試時,啟動節(jié)點機(jī)總數(shù)依次為16,32,64,128。由于兩級分區(qū)時相應(yīng)1級子區(qū)域總數(shù)應(yīng)分別等于每次啟動節(jié)點機(jī)總數(shù),故相應(yīng)1級子區(qū)域總數(shù)也依次為16,32,64,128。因“神威太湖之光”超級計算機(jī)每個節(jié)點包含4個異構(gòu)核組,故2級剖分時每個1級子區(qū)域被獨立地剖分為4個2級子區(qū)域。制動盤的并行計算結(jié)果如表1所示。

表1 制動盤并行效率計算結(jié)果Tab.1 Results of parallel computation for rotor of brake system

表1中,并行計算總時間包含各子區(qū)域系統(tǒng)開始計算到各子區(qū)域系統(tǒng)模態(tài)振型求解結(jié)束。1級子區(qū)域求解時間包括組集1級子區(qū)域等效剛度矩陣和等效質(zhì)量矩陣的時間、1級子區(qū)域凝聚縮減的時間、并行求解方程的時間和回代1級子區(qū)域模態(tài)坐標(biāo)的時間。由表1可知,相對于傳統(tǒng)的并行模態(tài)綜合法,本文提出的分層并行計算方法能夠獲得較高的加速比和并行效率。這是由于采用傳統(tǒng)的并行模態(tài)綜合法時,隨著子區(qū)域數(shù)目的增多,形成的整體廣義特征方程的規(guī)模和條件數(shù)也隨之急劇增加,導(dǎo)致了界面方程求解時間的大幅度增加,從而嚴(yán)重影響了系統(tǒng)總體并行效率的提高。從數(shù)學(xué)的角度看,分層并行計算方法1級子區(qū)域的求解實質(zhì)上是在求解和傳統(tǒng)并行模態(tài)綜合法相同規(guī)模的系統(tǒng)廣義特征方程,但其求解時間卻大大縮短了。例如在表1中,512計算核組時1級子區(qū)域的求解時間為732.1 s,而傳統(tǒng)并行模態(tài)綜合法的系統(tǒng)方程求解時間為2 021.2 s,這就節(jié)約了1 289.1 s。主要原因在于:分層并行計算方法通過兩級分區(qū)4次變換不僅將節(jié)點內(nèi)通信與節(jié)點間通信分離以及節(jié)點內(nèi)異構(gòu)核組內(nèi)通信與異構(gòu)核組間通信分離有效改善了通信效率,而且它還進(jìn)一步降低了凝聚后系統(tǒng)方程的規(guī)模,加快了其計算和迭代收斂速度。從負(fù)載均衡與異構(gòu)型超級計算機(jī)適配的角度來看,本文將有限元計算的結(jié)構(gòu)任務(wù)劃分為節(jié)點間并行、節(jié)點內(nèi)核組間并行、核組內(nèi)各計算核心間并行。由于初始1級子區(qū)域計算規(guī)模都相同,且各節(jié)點機(jī)具有相同的計算性能,這就實現(xiàn)了節(jié)點間的負(fù)載均衡和并行計算。2級子區(qū)域按照節(jié)點機(jī)內(nèi)核組的數(shù)量進(jìn)行劃分,與傳統(tǒng)并行模態(tài)綜合法相比,本文直接參與并行計算的子區(qū)域數(shù)目通過兩級分區(qū)降低為其1/4。子區(qū)域數(shù)目的降低一方面有利于降低界面方程的規(guī)模和條件數(shù),從而有效提高系統(tǒng)的數(shù)值收斂性;另一方面有利于大幅度減少參與并行計算的進(jìn)程數(shù)目,從而減少因通信競爭造成的各核組間的計算資源浪費。同時,核組內(nèi)各計算核組間的并行可以充分利用2級子區(qū)域?qū)?yīng)的所有計算核心共享的內(nèi)存空間,有效提高數(shù)據(jù)的訪問速度和計算效率,且能夠避免因分區(qū)過多導(dǎo)致的系統(tǒng)數(shù)值收斂性降低和通信開銷增加帶來的負(fù)載均衡問題。因此,它能夠有效減少1級子區(qū)域求解時間,并獲得良好的加速比和并行效率。

在實際的工程應(yīng)用中,有時復(fù)雜工程結(jié)構(gòu)會包含多種單元類型,為測試多單元混合建模千萬自由度規(guī)模下復(fù)雜工程系統(tǒng)的并行效率,以如圖7所示的某跨江隧道模型為例進(jìn)行分析,該模型具有2 896 781個實體單元、186 121個梁單元、21 685個質(zhì)量單元、自由度規(guī)模為13 167 203、剛度矩陣非零元個數(shù)為1 012 581 369、平均帶寬為412,求解其前20階固有模態(tài)頻率。

圖7 某跨江隧道結(jié)構(gòu)有限元系統(tǒng)Fig.7 FEM of the over-river tunnel structure system

并行計算時啟動節(jié)點機(jī)總數(shù)依次為32,64,128,256。采用傳統(tǒng)并行模態(tài)綜合法和本文多級分層并行計算方法的計算結(jié)果如表2所示。

表2 某跨江隧道的并行計算結(jié)果Tab.2 Results of parallel computation for over-river tunnel

由表2可知,對于傳統(tǒng)的并行模態(tài)綜合法,當(dāng)整體系統(tǒng)方程采用直接法進(jìn)行求解時,隨著子區(qū)域的增多整體系統(tǒng)方程的求解時間大幅度增加,從而嚴(yán)重降低了系統(tǒng)的并行效率。這是由于雖然總體界面方程采用稀疏列壓縮格式存儲,且三角分解時利用并行預(yù)處理Cholesky分解只對下三角部分進(jìn)行分解,但是總體界面方程是高度稠密的,其進(jìn)行三角分解時會增加原有方程的稠密度,致使其組集和三角分解不僅需要申請大量的內(nèi)存,而且需要大量的通信和計算。隨著子區(qū)域的增多,系統(tǒng)整體縮減后的廣義特征方程規(guī)模隨之增大,由此帶來的存儲、通信和計算開銷也就越來越大,系統(tǒng)整體的求解時間也就越長。而采用分層并行計算方法時,由于利用迭代法求解不需要組集形成系統(tǒng)縮減后的廣義特征方程。另外,涉及方程組并行求解時的局部通信僅存在相鄰子區(qū)域之間,只有少量的向量運算和整體迭代誤差的計算需要全局通信。因此,它能夠以較短的求解時間獲得較高的加速比和并行效率。

4 結(jié) 論

(1) 為解決采用模態(tài)綜合法利用異構(gòu)眾核分布式存儲并行計算機(jī)求解大規(guī)模有限元模態(tài)對計算效率造成的損失,在吸收和利用稀疏存儲技術(shù)和傳統(tǒng)并行模態(tài)綜合法優(yōu)點的基礎(chǔ)上提出了一種基于稀疏存儲格式的結(jié)構(gòu)有限元模態(tài)多級分層并行計算方法。

(2) 通過典型數(shù)值算例表明,與傳統(tǒng)的并行模態(tài)綜合法相比,該方法能夠獲得更高的加速度比和并行計算效率,并大幅度節(jié)省內(nèi)存空間。

(3) 本文的研究結(jié)果將為結(jié)構(gòu)有限元軟件在國產(chǎn)異構(gòu)眾核處理器和其他異構(gòu)處理器上的移植及大規(guī)模并行提高參考,對重大裝備系統(tǒng)及復(fù)雜工程系統(tǒng)的研制和使用均具有較強(qiáng)的指導(dǎo)意義和參考價值。