基于申威異構(gòu)眾核處理器架構(gòu)的模態(tài)并行算法

喻高遠, 馬志強,3, 李俊杰, 金先龍

(1.上海交通大學(xué) 機械系統(tǒng)與振動國家重點實驗室,上海 200240； 2.上海交通大學(xué) 機械與動力工程學(xué)院,上海 200240;3.中國航發(fā)商用航空發(fā)動機有限責(zé)任公司,上海 200240)

隨著交通運輸業(yè)、能源勘探與開發(fā)業(yè)和航空航天業(yè)等的發(fā)展，對于大型、特大型裝備的需求越來越多，如：高速動車組、3 000 m超深鉆機、大飛機、跨江隧道、跨海大橋等。這些特殊裝備系統(tǒng)的研制往往涉及大規(guī)模復(fù)雜動力學(xué)系統(tǒng)的計算，而模態(tài)分析則是其最耗費時間的計算環(huán)節(jié)，也是其余計算環(huán)節(jié)的基礎(chǔ)，需借助大規(guī)模有限元模型進行高性能計算，故而對傳統(tǒng)串行有限元計算方法和工具形成了挑戰(zhàn)[1-2]。傳統(tǒng)串行計算是以犧牲大型、特大型裝備局部關(guān)鍵細節(jié)進行簡化建模來保證計算效率，因而造成局部關(guān)鍵細節(jié)預(yù)測能力和大量密集模態(tài)的丟失，計算精度較低，無法滿足其系統(tǒng)級高精度高效率數(shù)值分析的需求。隨著并行計算機的快速發(fā)展，利用并行計算機研究和開發(fā)相應(yīng)的并行算法則為大型、特大型裝備系統(tǒng)模態(tài)的求解提供了切實可行的途徑，正逐步成為各國學(xué)者的研究熱點。

在硬件方面，異構(gòu)眾核分布式存儲并行計算機具備計算能力強、性能功耗比高等優(yōu)點，已成為當(dāng)前超級計算機的重要發(fā)展方向,典型的異構(gòu)眾核處理器包括Intel的MIC、Nvidia和AMD的GPU、Godson-T以及申威眾核處理器等[3]。近年來，國內(nèi)外諸多學(xué)者在異構(gòu)眾核分布式存儲并行計算機的基礎(chǔ)上求解各類大規(guī)模、超大規(guī)模有限元系統(tǒng)，來獲取系統(tǒng)的特性，取得了很好的效果。Koric等[4]利用并行SuperLU和PCG(preconditional conjugate gradient)算法基于GPU眾核架構(gòu)完成了某增壓空氣冷卻器的瞬態(tài)動力學(xué)特性分析，其求解自由度超過千萬。Martínez-Frutos等[5]采用CG(conjugate gradient)算法基于GPU眾核架構(gòu)完成了某L型懸臂梁的靜態(tài)特性分析。楊梅芳等[6]使用直接數(shù)值模擬基于MIC眾核架構(gòu)完成了某發(fā)動機燃燒模擬分析。然而，國內(nèi)外學(xué)者關(guān)于模態(tài)并行計算的相關(guān)研究較少，且多是以多核并行計算機和基于GPU眾核架構(gòu)的并行計算機為主。Heng等[7]基于多核并行計算機完成了模態(tài)疊加法并行算法設(shè)計，并將其應(yīng)用于某懸臂梁的模態(tài)并行求解。朱彬等基于GPU眾核架構(gòu)設(shè)計了模態(tài)并行子空間迭代法，基于此完成了某風(fēng)扇結(jié)構(gòu)模態(tài)分分析。目前國內(nèi)外基于申威眾核處理器架構(gòu)的模態(tài)并行求解算法研究相對較少，而基于申威眾核處理器架構(gòu)的并行計算機“神威太湖之光”在峰值性能、持續(xù)性能、性能功耗比3項關(guān)鍵指標(biāo)均居于世界第一[8]。因此，利用基于申威眾核處理器架構(gòu)的并行計算機進行模態(tài)并行計算研究對于提高大型、超大型裝備系統(tǒng)模態(tài)的計算規(guī)模、計算精度和計算效率具有重要意義。考慮到“神威太湖之光”并行計算機核組內(nèi)通信時間遠小于核組間通信，且其訪存能力較弱，故利用“神威太湖之光”并行計算機提高并行效率的關(guān)鍵在于處理好大規(guī)模數(shù)據(jù)的存儲以及各計算核心間的通信和協(xié)作問題。模態(tài)分析的數(shù)學(xué)實質(zhì)可以歸結(jié)為大型稀疏矩陣的廣義特征值問題，該類問題的求解大多基于子空間類投影技術(shù)，主要包括Davidson類子空間方法和Krylov類子空間方法等。Davidson類子空間方法主要用于求解對角占優(yōu)的對稱矩陣特征值問題，其問題適應(yīng)性不如Krylov類子空間方法。Krylov類子空間方法可以追溯到20世紀50年代提出的Lanczos算法和Arnoldi算法[9]。后來國內(nèi)外諸多學(xué)者在Lanczos算法和Arnoldi算法的基礎(chǔ)上進行了一系列重啟動改進，比較著名的是：Sorensen等[10]提出的Arnoldi/Lanczos算法、Stewart等[11]提出的Krylov-Schur算法、Jia等[12-15]提出的加速子空間迭代法等。3種算法在數(shù)學(xué)上具有等價性，是目前Krylov類自子空間算法中的主流算法。與前兩種算法相比，加速子空間迭代法更容易收斂，且代碼實現(xiàn)難度較低，故本文采用加速子空間迭代法進行模態(tài)并行算法設(shè)計。

綜上所述，本文基于國產(chǎn)申威異構(gòu)眾核分布式存儲并行計算機和加速子空間迭代法分析了各計算步驟的計算量，根據(jù)計算結(jié)果構(gòu)建了大規(guī)模模態(tài)分析并行計算體系，并將其應(yīng)用于某超深鉆井制動系統(tǒng)主體機構(gòu)及某跨江隧道模態(tài)并行計算，實現(xiàn)了上千萬自由度的模態(tài)并行求解。同時，該方法不僅通過分層策略實現(xiàn)了計算過程和數(shù)據(jù)通信的分層，有效提高了通信效率；而且通過計算數(shù)據(jù)的分布式存儲，顯著改善了數(shù)據(jù)訪存效率。

1 大規(guī)模特征值問題求解算法

模態(tài)分析的數(shù)學(xué)的描述為

Kφ=λMφ

(1)

式中：K為模態(tài)系統(tǒng)整體剛度矩陣；M為模態(tài)系統(tǒng)整體質(zhì)量矩陣；λ為模態(tài)系統(tǒng)廣義特征值；φ為對應(yīng)振型向量。K和M可以對工程結(jié)構(gòu)進行有限元離散和積分得到，均為大型稀疏、對稱(半)正定矩陣。模態(tài)分析的本質(zhì)即求解式(1)的多個低階特征對。采用子空間迭代法求解式(1)時，由于Krylov類算法大多數(shù)收斂于最大特征值，需采用Shift-Invert變換進行譜變換，其變換形式為

(2)

式中：σ為移位值；(K-σM)-1可通過變換求解線性系統(tǒng)的解獲得

(K-σM)x=M

(3)

式(2)可改寫為

Asν=μsν

(4)

式中：As=(K-σM)-1；μs=1/(λ-σ)。采用加速子空間迭代法求解式(1)的前m個特征值即求解式(4)的前m個特征值時，考慮到As的存儲數(shù)據(jù)量為自由度規(guī)模n×n，為了最大限度降低中間變量As儲存的內(nèi)存占用空間，變量V(V可以取算法過程變量Q或者Y)與As做矩陣運算后的結(jié)果可通過求解式(5)所示的線性系統(tǒng)獲得，其算法具體步驟如下所示。

(K-σM)(xV)=(MV)

(5)

步驟1輸入矩陣K、M，求解特征值個數(shù)m，迭代初始向量Q，外層迭代控制誤差ε，最大循環(huán)次數(shù)Maxcycle。

步驟2輸出m個外部特征值λj和w。

(1) 初始化

隨機生成初始向量Q,j=0,Y=[·],AA=[·],BB=[·],VV=[·],EE=[·],BBB=[·],LL=[·]。

(2) 進入求解m個特征值的循環(huán)

whilej

①計算:求解方程(K-σM)×(As×Q)=(M×Q)并將結(jié)果存儲于Y中；

②計算:求解方程(K-σM)×(As×Y)=(M×Y)并將結(jié)果存儲于Q中；

③計算:AA=Y′×Q；

④計算:BB=Y′×(E×Y)；

⑤QZ法求解子空間上廣義特征值問題：AA×ν=λj×BB×ν,式中,λj為第j次迭代求得的廣義特征值；

⑥檢查λj是否滿足精度要求,若|(λj-λj-1)/λj|≤ε,則轉(zhuǎn)到步驟4;如果不滿足精度要求則作

BBB=((VV′×BB×VV)′+(VV′×BB×VV))/2

VV為ν構(gòu)成的向量空間；

⑦對BBB做Cholesky分解并將產(chǎn)生上三角陣LL；

⑧計算:VV=VV/LL′；

⑨計算:Q=Y×VV；

⑩令j=j+1并返回步驟①；

end while

(3)檢查

j

若滿足，則轉(zhuǎn)到步驟(4);如果不滿足則輸出計算有誤；

(4)計算

w=sqrt(λj)/(2π)并輸出λj和w。

2 模態(tài)加速子空間迭代法分析

采用第1章算法，進行某超深鉆進盤鼓式制動器轉(zhuǎn)子盤模態(tài)分析，其有限元網(wǎng)格模型如圖1所示，彈性模量為210 GPa，密度為7 800 kg/m3，泊松比為0.3。采用不同的自由度規(guī)模，固定約束其內(nèi)表面8個螺栓孔位置，計算結(jié)構(gòu)的前20階固有頻率，并與經(jīng)典模態(tài)求解算法-Lanczos算法[16-17]的求解結(jié)果進行對比，各測試規(guī)模如表1所示，各自由度規(guī)模下20階固有頻率的最大相對誤差按照式(6)計算后如圖2所示。

圖1 盤鼓式制動器轉(zhuǎn)子盤有限元網(wǎng)格模型

(6)

由圖2可知，各自由度規(guī)模下各階固有頻率的誤差均不超過0.11%，表明：基于加速子空間迭代法的模態(tài)分析可以保證計算精度，故可用于模態(tài)并行算法設(shè)計。

圖2 不同自由度規(guī)模下20階固有頻率的最大相對誤差

不同測試規(guī)模下各步驟的時間比例如圖3所示，由圖3可知，隨著自由度規(guī)模的增加，數(shù)據(jù)讀取和特征值計算的時間比例逐步下降，初始化及方程求解的時間逐步增加，故而大規(guī)模模態(tài)并行算法設(shè)計的關(guān)鍵在于方程求解步驟的并行化，而特征值值計算可采用單節(jié)點計算，以減少通信耗時。

圖3 轉(zhuǎn)子盤不同測試規(guī)模下各步驟的時間比例

3 并行計算實現(xiàn)

3.1 處理器架構(gòu)

采用用申威眾核處理器進行模態(tài)并行加速子空間迭代法設(shè)計，其架構(gòu)如圖4所示。

圖4 申威眾核處理器架構(gòu)

每個申威眾核處理器,共計4核組,各核組可共享32 GB內(nèi)存。每個核組包括1個主核(運算控制核心)和64個從核(核心陣列)。核組間通信采用雙向14 Gbits/s通信網(wǎng)絡(luò)帶寬,主核與從核間通信采用DMA(direct memory access)方式批量訪問主存。從核局部存儲空間大小為64 kB,指令存儲空間為16 kB。

3.2 大規(guī)模模態(tài)并行求解體系

基于申威眾核處理器及接口等功能形成的模態(tài)分析求解體系，如圖5所示。

圖5 模態(tài)分析并行計算體系

整個模態(tài)分析并行計算體系分為：多文件流數(shù)據(jù)讀取、變量初始化及并行求解方程、并行求解模態(tài)固有頻率3個部分。具體介紹如下：

(1) 多文件流數(shù)據(jù)讀取。各核組同步讀取對應(yīng)的剛度矩陣和質(zhì)量矩陣數(shù)據(jù)文件，核組間無數(shù)據(jù)通信交流。剛度矩陣和質(zhì)量矩陣是由組集系統(tǒng)模型的各部分結(jié)構(gòu)化網(wǎng)格數(shù)據(jù)后并行求解獲得[18]。

(2) 并行求解式(5)。考慮到大規(guī)模模態(tài)并行求解時需要求解式(5)兩次，為了節(jié)約方程的求解時間，通過在申威眾核處理器上集成并行LU算法來實現(xiàn)線性方程組的求解，在求解過程中組裝的系統(tǒng)的單元剛度矩陣K僅需要進行一次LU分解，因而可以節(jié)約式(5)的總體求解時間。并行LU算法的實現(xiàn)過程如圖6所示，主要包括矩陣并行Cholesky算法和三角線性方程組并行求解算法，其核心運算步驟為矩陣向量運算和數(shù)據(jù)通信。

圖6 并行LU求解算法

數(shù)據(jù)通信包括核組間通信以及核組內(nèi)通信，核組間通信采用MPI庫實現(xiàn)，核組內(nèi)通信采用Athread庫實現(xiàn)。矩陣向量運算主要包含加減乘除，現(xiàn)以向量乘法a=b·c為例(a、b、c為任意矩陣向量運算過程中的存儲數(shù)組)，其實現(xiàn)過程如圖7所示。各核組上64個從核同步從核組內(nèi)存空間中循環(huán)讀取對應(yīng)數(shù)據(jù)，該部分數(shù)據(jù)段內(nèi)存需小于64 kB進行計算后返回計算結(jié)果于指定位置，通信僅存在于各核組主核于從核之間。

圖7 基于異構(gòu)眾核加速的矩陣向量乘法

(3) 并行加速子空間算法求解模態(tài)固有頻率。按照算法操作屬性主要包括：矩陣向量運算、QZ法并行求解廣義特征值問題以及并行Cholesky分解等。進行矩陣向量運算的步驟主要為①～④、⑥、⑧，且QZ法并行求解廣義特征值問題及Cholesky分解中均存在矩陣向量運算，其實現(xiàn)過程同圖7所示。考慮到計算規(guī)模對于加速子空間各步驟時間占比的影響，步驟①～步驟④中的矩陣向量運算需各核組同步并行計算，存在核組間通信和核組內(nèi)通信，QZ法并行求解廣義特征值問題及Cholesky分解中的矩陣向量計算只在指定核組內(nèi)進行運算，僅存在核組內(nèi)通信。QZ法并行求解廣義特征值問題的算法及并行Cholesky分解算法的實現(xiàn)如圖8所示。

圖8 QZ及Cholesky分解算法的并行化實現(xiàn)

4 數(shù)值算例

基于搭建的模態(tài)并行求解體系完成某超深鉆進盤鼓式制動器轉(zhuǎn)子盤及某跨江隧道模態(tài)分析，求解模態(tài)階數(shù)為前10階。

4.1 某制動裝備典型應(yīng)用

對于該超深鉆進盤鼓式制動器轉(zhuǎn)子盤有限元模型，整體結(jié)構(gòu)測試規(guī)模如表2所示。

表2 超深鉆進盤鼓式制動器轉(zhuǎn)子盤測試規(guī)模

為揭示復(fù)雜裝備大規(guī)模模態(tài)計算的必要性及求解精度，對方案6～方案8的模態(tài)頻率結(jié)果進行了比較分析，表3給出了3種規(guī)模下模態(tài)頻率的對比情況。

由表3可知，同類型裝置的各階模態(tài)頻率隨著計算規(guī)模的增加，頻率會逐漸下降，這是由于工程結(jié)構(gòu)的有限元分析存在剛度矩陣的“硬化效應(yīng)”，導(dǎo)致較小自由度規(guī)模計算時得到的模態(tài)頻率偏高。對比方案6～方案8與方案1～方案5的模態(tài)頻率，最大變化率4.39%，而方案6～方案8的模態(tài)頻率變化率相對低的多，這說明對于類似于盤鼓式制動器轉(zhuǎn)子盤這樣的結(jié)構(gòu)，需要提高相應(yīng)的計算規(guī)模以提高其計算精度。

表3 轉(zhuǎn)子盤不同測試規(guī)模下的結(jié)果變化

為了校驗并行計算結(jié)果的正確性，圖9給出了3種規(guī)模下模態(tài)頻率與經(jīng)典Lanczos算法對應(yīng)的相對誤差，由圖9可知，3種規(guī)模下模態(tài)頻率的計算結(jié)果與經(jīng)典Lancozos算法的相對誤差均小于0.687%，各階振型保持一致，這就有效驗證了本文并行計算結(jié)果的正確性。

圖9 不同自由度規(guī)模下模態(tài)頻率的相對誤差

通過啟動相應(yīng)數(shù)目的節(jié)點機測試本文分層并行計算方法的性能，各規(guī)模下的計算結(jié)果如表4～表6所示。

由表4～表6可知，本文基于申威眾核處理器架構(gòu)提出的有限元模態(tài)分層通信并行計算方法能夠獲得較高的加速比和并行效率。這是由于分層通信策略實現(xiàn)了計算過程和數(shù)據(jù)通信的分層，各核組的主核僅負責(zé)數(shù)據(jù)的讀取和全局通信，而各從核負責(zé)計算且僅與對應(yīng)主核之間存在局部通信，因而可獲得良好的加速比和加速效率。同時，為了進一步降低全局通信的次數(shù)和時間，考慮到組集形成的單元整體質(zhì)量矩陣M需進行反復(fù)調(diào)用計算，將質(zhì)量矩陣數(shù)據(jù)直接存儲于各個核組存儲空間中，可使得實際M參與的計算過程中，僅需要中間少量數(shù)據(jù)通信，而組集形成的剛度矩陣K需要進行LU分解，將其分布式存儲于各個核組上，盡管這一過程會使得數(shù)據(jù)讀取的并行可擴展性較差，然而對于大規(guī)模模態(tài)并行求解，數(shù)據(jù)讀取占總時間的比例相對較低，其并行求解的關(guān)鍵在于如何降低全局的通信量，故可提高整體的加速比和并行效率。

表4 方案6并行計算結(jié)果

表5 方案7并行計算結(jié)果

表6 方案8并行計算結(jié)果

在表4～表6中，當(dāng)計算核數(shù)由4 096增加到16 384的過程中，系統(tǒng)的總體并行效率發(fā)生了顯著下降，這是由于模態(tài)計算過程中隨著計算分區(qū)的增加，采用并行LU求解式(5)的過程中占用的內(nèi)存和并行通訊也隨之增加，導(dǎo)致系統(tǒng)的整體并行效率下降較大。

4.2 某跨江隧道典型應(yīng)用

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

圖10 某跨江隧道主體有限元模型

各計算核數(shù)下的測試結(jié)果如表7所示。由表7可知，對于包含多種類型單元千萬自由度規(guī)模下的跨江隧道系統(tǒng)模態(tài)并行求解，本文所提出的模態(tài)并行計算方法仍然具有良好的加速比和并行效率。然而計算核數(shù)從12 544核提升至16 384核時，總體計算時間僅下降了715.2 s，下降幅度較小，這是由于LU分解和下(上)三角部分的并行求解不僅需要申請大量的內(nèi)存，還需要大量的通信和計算。隨著子區(qū)域的增加，盡管單個區(qū)域的計算時間有所降低，然而求解過程中通信開銷會越來越大，致使系統(tǒng)總體計算時間降低較少，并行效率較低。與制動裝備轉(zhuǎn)子結(jié)構(gòu)相比，該跨江隧道系統(tǒng)模型包含了多種類型的單元，可充分構(gòu)建復(fù)雜工程系統(tǒng)的多單元混合模型，且仍然具備良好的并行效率和加速比，因而更加符合復(fù)雜工程系統(tǒng)的實際需求。

表7 某跨江隧道并行計算結(jié)果

5 結(jié) 論

(1) 為解決采用國產(chǎn)申威異構(gòu)眾核處理求解重大裝備系統(tǒng)級模態(tài)并行求解問題，通過對模態(tài)加速子空間迭代法進行研究和集成開發(fā)，構(gòu)建了一套大規(guī)模并行求解體系，該體系在模態(tài)加速子空間迭代法的基礎(chǔ)上，利用分層通信策略不僅實現(xiàn)了計算過程和數(shù)據(jù)通信的分層，有效改善了通信效率，而且實現(xiàn)了計算過程中數(shù)據(jù)的分布式存儲，顯著改善了數(shù)據(jù)訪存效率。

(2) 通過典型數(shù)值算例表明：對于完全實體單元建模和多單元混合建模下的系統(tǒng)級模型，該方法均能夠獲得良好的加速比和并行效率。并且具備上千萬自由度規(guī)模的并行計算能力，基本可以滿足重大裝備和復(fù)雜工程系統(tǒng)利用國產(chǎn)處理器進行模態(tài)分析的需求。本文的研究結(jié)果對于重大裝備和復(fù)雜工程系統(tǒng)的研制和使用均具有較強的指導(dǎo)意義和參考價值。