結構模態多級分層并行計算方法

喻高遠, 樓云鋒, 李俊杰, 金先龍

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

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

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

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

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

1 模態綜合法

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

(1)

式中:UI,UB,{pk},{pB}分別為內部節點和邊界節點對應的位移和模態坐標;[Φ]為子區域的模態轉化矩陣;[ΦIk]為內部節點主模態矩陣;[ΦIB]為約束模態。

利用各子區域系統之間的位移協調條件,進行第2次坐標變換,如式(2)所示

(2)

K*φ*=λ*M*φ*

(3)

縮減剛度矩陣[K*]和縮減質量矩陣[M*]為

(4)

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

2 分層并行計算方法

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

2.1 并行任務映射

通過將計算任務映射到異構眾核超級計算機的不同硬件層,以實現不同層級的負載均衡和通信的有效分割。在傳統并行模態綜合法的基礎上,考慮“神威太湖之光”異構眾核超級計算機的硬件體系架構,完成結構大規模模態分層并行計算任務映射,如圖1所示。

圖1 結構模態多級分層并行模態綜合法計算任務映射Fig.1 Task mapping of multilevel hierarchical parallel modal synthesis algorithm

在進行任務映射時,將第1級網格子區域按照節點順序進行映射,第2級網格子區域按照節點內異構群組進行映射,各異構群組內浮點運算按照計算核心進行映射。

2.2 兩級分區

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

圖2 2級分區Fig.2 Two-level partitioning

2.3 4次變換

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

圖3 4次變換Fig.3 Four transformations

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

2.4 三層并行計算機制

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

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

為節省內存空間和減少計算量,各2級子區域局部總體剛度矩陣和總體質量矩陣均采用列壓縮稀疏存儲格式進行存儲。

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

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

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

3 算例驗證

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

圖6 某轉子盤有限元系統Fig.6 Finite element of rotor system

(5)

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

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

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

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

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

圖7 某跨江隧道結構有限元系統Fig.7 FEM of the over-river tunnel structure system

并行計算時啟動節點機總數依次為32,64,128,256。采用傳統并行模態綜合法和本文多級分層并行計算方法的計算結果如表2所示。

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

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

4 結 論

(1) 為解決采用模態綜合法利用異構眾核分布式存儲并行計算機求解大規模有限元模態對計算效率造成的損失,在吸收和利用稀疏存儲技術和傳統并行模態綜合法優點的基礎上提出了一種基于稀疏存儲格式的結構有限元模態多級分層并行計算方法。

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

(3) 本文的研究結果將為結構有限元軟件在國產異構眾核處理器和其他異構處理器上的移植及大規模并行提高參考,對重大裝備系統及復雜工程系統的研制和使用均具有較強的指導意義和參考價值。