基于申威異構眾核處理器架構的模態并行算法

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

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

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

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

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

1 大規模特征值問題求解算法

模態分析的數學的描述為

Kφ=λMφ

(1)

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

(2)

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

(K-σM)x=M

(3)

式(2)可改寫為

Asν=μsν

(4)

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

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

(5)

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

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

(1) 初始化

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

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

whilej

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

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

③計算:AA=Y′×Q；

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

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

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

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

VV為ν構成的向量空間；

⑦對BBB做Cholesky分解并將產生上三角陣LL；

⑧計算:VV=VV/LL′；

⑨計算:Q=Y×VV；

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

end while

(3)檢查

j

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

(4)計算

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

2 模態加速子空間迭代法分析

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

圖1 盤鼓式制動器轉子盤有限元網格模型

(6)

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

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

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

圖3 轉子盤不同測試規模下各步驟的時間比例

3 并行計算實現

3.1 處理器架構

采用用申威眾核處理器進行模態并行加速子空間迭代法設計，其架構如圖4所示。

圖4 申威眾核處理器架構

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

3.2 大規模模態并行求解體系

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

圖5 模態分析并行計算體系

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

(1) 多文件流數據讀取。各核組同步讀取對應的剛度矩陣和質量矩陣數據文件，核組間無數據通信交流。剛度矩陣和質量矩陣是由組集系統模型的各部分結構化網格數據后并行求解獲得[18]。

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

圖6 并行LU求解算法

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

圖7 基于異構眾核加速的矩陣向量乘法

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

圖8 QZ及Cholesky分解算法的并行化實現

4 數值算例

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

4.1 某制動裝備典型應用

對于該超深鉆進盤鼓式制動器轉子盤有限元模型，整體結構測試規模如表2所示。

表2 超深鉆進盤鼓式制動器轉子盤測試規模

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

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

表3 轉子盤不同測試規模下的結果變化

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

圖9 不同自由度規模下模態頻率的相對誤差

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

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

表4 方案6并行計算結果

表5 方案7并行計算結果

表6 方案8并行計算結果

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

4.2 某跨江隧道典型應用

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

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

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

表7 某跨江隧道并行計算結果

5 結 論

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

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