面向異構超算的結構分析高效并行計算方法*

丁峻宏,苗新強,李根國

(1.上海超算科技有限公司，上海 201203;2.中國電子科技集團公司第五十一研究所，上海 201802;3.上海超級計算中心，上海 201203)

1 引言

當前，在很多科學和工程研究領域為了提高計算分析的規模和效率，通常采用超級計算機進行并行計算[1 - 3]。其中，有限元和計算機技術的融合和快速發展使其被廣泛應用于航空、汽車、土木和建筑等許多工程領域。在進行大規模并行計算時，往往需要根據超級計算機的硬件體系結構特點設計出合適的并行算法，以充分發揮其高效的計算性能。

傳統的超級計算機在每個結點內配置的計算設備比較單一，一般只配置CPU，被稱為同構型超級計算機[4]。近年來硬件加速技術的飛速發展使得新型研制的超級計算機傾向于在每個結點內部配置不同性能的計算設備，除了配置CPU外，還配置有協處理器，如圖形處理器GPU(Graphic Processing Unit)[5]和眾核處理器MIC(Many Integrated Core)[6]，這也就催生了異構型超級計算機。對于傳統的同構型超級計算機，國內外一些學者提出了不少好的方法。例如，針對共享存儲的同構型超級計算機提出了循環級的細粒度并行計算方法[7]，而針對分布式存儲的同構型超級計算機提出了任務級的粗粒度并行計算方法[8]。在異構超級計算機中最早是采用GPU作為協處理器，很多學者也對其有一定的研究[9 - 11]。異構型超級計算機從采用Intel公司推出的MIC作為協處理器開始，將并行計算從多核時代帶入了眾核時代[12]。對適合異構眾核超算的并行計算方法的研究目前處于不斷探索階段[13 - 15]，相關的研究還不是太多。

本文嘗試探索異構眾核架構下有限元結構分析的并行優化方法，研究有限元結構分析與高性能計算平臺的最優適配方法，為同類應用問題的并行移植和性能優化提供借鑒與參考。

2 結構分析并行計算與異構眾核計算平臺

2.1 結構分析并行計算

結構分析平衡方程可以表達為：

Kx=P

(1)

其中，K為總體剛度矩陣，P為總體外部載荷向量，x為位移向量。

利用區域分解法進行有限元結構分析并行計算的基本原理[16]如下所示：

首先將總體網格模型分為數個子區域；然后按照先內部自由度后邊界自由度的編號原則同時獨立形成每個子區域的系統方程：

(2)

其中，xI,xB,PI,PB分別為內部節點和邊界節點對應的位移和外部載荷向量;K**為剛度矩陣中的分塊矩陣，其中下標I和B分別表示內部自由度和邊界自由度。

通過縮聚同時獨立消去各子區域內部節點自由度后，得到只含邊界節點自由度未知量的界面方程:

(3)

其中有效剛度矩陣:

(4)

有效載荷向量:

(5)

接著將所有子區域的界面方程聯立求解，得到各子區域邊界節點位移。根據求得的邊界節點位移xB，各子區域內部節點位移xI由式(6)回代求解:

xI=(KII)-1(PI-KIBxB)

(6)

計算各子區域的應變值和應力值，并分別輸出結果。

2.2 異構眾核計算平臺

“天河二號”異構眾核超算包含1.6萬計算結點(單結點配置2顆Intel Xeon E5 多核處理器和3塊Intel Xeon Phi眾核處理器)。結點間互連采用自制主干拓撲結構網絡，網絡接口使用NIC控制器。

傳統區域分解法應用在同構型超級計算機上效果很好，但應用在異構眾核超級計算機上時卻會面臨一些問題。這是由于異構眾核超級計算機在大幅度提升系統計算性能的同時也為高效并行計算方法的設計帶來了一定的挑戰：首先MIC的計算性能與CPU相比迥異，因此如何實現CPU與MIC設備之間以及不同結點之間的負載均衡成為迫切需要解決的問題。其次，在異構眾核超級計算機中不同結點、設備以及計算核心間的通信延遲差異顯得更加突出，而通信延遲是制約大規模并行計算效率的主要因素之一，故如何提高系統的通信效率成為并行計算方法設計的一個關鍵問題。最后，每個MIC卡內存空間不大，如何利用有限的內存空間實現大規模并發性計算任務也是需要重點考慮的問題。

3 多層次多粒度協同并行計算方法

從硬件體系結構看，眾核架構下異構超算由多計算單元(多核或眾核)、多設備(CPU或MIC卡)和多結點的層次性元素構成。在結點之間和設備之間可以考慮使用任務級的粗粒度并行計算方法，而多核或眾核的不同計算核心間適合采用循環級的細粒度并行計算方法。為了實現并行計算方法與異構超算平臺的最優適配，本文提出了一種針對結構有限元分析的多層次多粒度協同并行計算方法。

3.1 總體方法

多層次多粒度協同并行計算方法基于對計算任務的層次性和粒度性剖析。如圖1所示，本文將有限元結構分析的計算作業的并行處理劃分為3個不同層面，以與異構眾核超算硬件結構相匹配：結點間、設備間和核間。其中前兩者使用粗粒度并行方法，而后者使用細粒度并行方法。

Figure 1 Task mapping圖1 任務映射

(1)結點間并行。

任務級的粗粒度并行計算方法適合結點間并行，這是由于采用分布式存儲的異構眾核超算的結點之間通信延遲遠遠高于結點內部，因此任務級的粗粒度并行計算方法適合于結點之間的并行。這樣不僅能夠實現大規模并行計算數據的分布式存儲，而且可以有效降低通信開銷。

如圖1所示，結點間并行在區域分解法基礎上開展：首先結構有限元網格被劃分為M個一級子區域，其中M為并行計算所調用的整體結點數量。每個一級子區域被分別分配給一個計算結點，從而實現結點間的并行計算。由于各計算結點的性能都一致，因此為實現結點間的負載均衡，在進行一級剖分時各子區域的規模也應當相同。

(2)設備間并行。

相對于結點間通信來說，同一結點內不同設備間的通信開銷要小一些，但比每個設備內部的核間通信開銷還大很多，所以基于區域分解的粗粒度并行計算方法多用于不同設備之間。為最大程度限制設備間的通信在同一計算結點內，應將相鄰的子區域分配在一起。

如圖1所示，設備間的并行立足于(1)中的一級分區結果：每個一級子區域被進一步分解為(Q+R)個二級子區域，其中Q指的是單個結點內部CPU設備的數目，R則定義為單個結點內部MIC卡設備的數目。然后，每一個設備將單獨處理一個二級子區域，從而實現設備間的并行計算。由于各二級子區域均由同一個一級子區域派生而來，這就將相鄰的子區域分配在了一起，這樣就使得設備之間的通信被限制在同一計算結點內。考慮到MIC和CPU計算性能的表現相差明顯，因此在二級分區時分配給CPU設備和MIC設備的子區域規模應有所不同。為盡量在不同設備間尋求負載均衡，本文將特定的計算任務分別分配給CPU和MIC單獨計算，以確定二者計算性能的比值，然后據此確定分配給CPU設備和MIC設備的子區域規模的大小。

(3)核間并行。

隨著硬件技術的快速發展，異構眾核超級計算機在每個結點內集成的計算核心數也越來越多。例如，“天河二號”上每顆CPU配置12個多核計算核心，每個MIC卡配置61個眾核計算核心(其中57個可供用戶使用)。如此眾多的計算核心一方面具有較高的硬件并發計算特性，另一方面這些計算核心間通過共享系統二級緩存和內存具有很高的數據通信效率。若采用基于區域分解的粗粒度并行計算方法利用這些多核或眾核資源難免會面臨因分區過多而導致系統通信開銷增加的問題。而細粒度并行計算方法本身就具備高度的并發性，能夠有效利用多核或眾核資源較高的硬件并發計算特性。此外，它是基于共享存儲實現的，能充分利用多核或眾核間共享二級緩存和內存的優勢提高系統數據訪問效率[17]。因此，本文采用循環級的細粒度并行計算方法實現核間的并行計算。

核間并行將基于每一個二級子區域的計算任務繼續作分解，這將觸及到處于最底層的計算模塊。基于區域分解法進行有限元結構分析并行計算的主要步驟包括：組集子區域系統平衡方程、縮聚、求解界面方程、回代內部自由度以及計算子區域應變和應力。如圖1所示，核間并行要做的工作就是基于二級子區域的每個計算步驟，搜尋其中的熱點計算程序，即大的循環結構，并將其進一步分解為一批彼此獨立、各自執行的子任務，其后CPU或MIC的一個計算核心將被調配去執行每一個子任務。

3.2 多層次多粒度協同并行計算方法的實現

本文采用offload模式實現有限元結構分析的多層次多粒度協同并行計算方法，在通信層面使用MPI+OpenMP實現。如圖2所示，有限元結構分析的多層次多粒度協同并行計算流程為：

第1步首先實施兩級分區方法，以獲取用于結構分析并行所需的計算數據文件，包括M個子區域的節點、單元、載荷和邊界條件以及相鄰分區信息等。

第2步在CPU端啟動執行M*(Q+R)個MPI進程，其中每個進程負責讀取和處理一個二級子區域數據文件。

第3步在每一個進程的內部再各自進一步派生出K個線程，并加載到CPU端利用多核資源完成相應子區域總體剛度矩陣和外部載荷向量的計算。其中，K等于單個結點內所有CPU設備的總核數除以單個結點內的進程總數并取整。

第4步與MIC設備相關聯的進程首先將子區域的總體剛度矩陣和外部載荷向量加載到MIC卡，在MIC端進一步派生出S個線程并進行縮聚計算，計算結果將最終被傳回至CPU端；而與CPU設備相關聯的進程則直接在CPU端派生T個線程進行縮聚計算。其中，S等于單個MIC卡總核數的4倍；T等于單顆CPU的總核數。

第5步各進程共同利用并行預處理共軛梯度PCG(Preconditioned Conjugate Gradient)算法求解界面方程，得到各子區域邊界節點位移后再回代求解內部位移。其中，求解界面方程時每個進程內部派生出K個線程，以有效利用CPU的多核資源縮短計算量較大的任務的處理時間。

第6步與MIC設備相關聯的進程首先將子區域節點位移加載到MIC卡，然后在MIC端派生S個線程回代求解子區域應變/應力，計算結果將被傳回至CPU端；而與CPU設備相關聯的進程則直接在CPU端派生T個線程計算子區域應變/應力。

第7步如果迭代繼續，則流程將跳轉至第2步再執行一遍;否則該流程結束。

Figure 2 Flowchart of collaborative parallel computing based on multi-layer and multi-grain圖2 多層次多粒度協同并行計算流程

3.3 關鍵技術研究

為保證多層次多粒度并行計算方法的順利實施，本文主要針對MIC卡的內存瓶頸問題、縮聚算法的高度并發性實現和界面方程的并行求解等核心問題進行研究。

3.3.1 MIC卡的內存瓶頸問題

傳統有限元分析程序一般采用變帶寬格式存儲結構剛度矩陣，由于它不能避免對帶寬內大量零元素的存儲，因此隨著有限元規模的擴大會導致系統內存需求的急劇增加。但是，單個MIC卡配置的內存空間不太大，要利用眾核資源進行并行計算首先必須降低系統的內存需求。

近年來，最新發展起來的行壓縮存儲技術通過僅針對剛度矩陣中的非零元素進行存儲，可以大幅度減少系統的內存需求。考慮到CPU設備往往配置有大量的內存空間，并且這些內存空間對同一結點內的所有CPU設備都是共享的，因此本文首先在CPU端采用變帶寬格式完成各二級子區域總體剛度矩陣和外部載荷向量的組集，將其轉化為占用內存空間較小的行壓縮存儲格式后再加載到MIC端利用眾核資源進行并行計算。

3.3.2 縮聚算法的高度并發性實現

本文通過將縮聚計算轉換為一系列相互獨立的線性方程組的求解，并利用OpenMP來開展并行化處理，從而實現了整個計算任務的高度并發性執行。本文設計的縮聚算法的偽代碼如算法1所示。

算法1縮聚算法

1.//calculate the condensed stiffness matrix

2.!$ompparalleldo

3.doi=1,z

4.b=KIB(:,i)

5. solveKIIt1=bwith a direct sparse solver

6. computet2=KBIt1

8.enddo

9.!$ompendparalleldo

10.//calculate the condensed load vector

11.b=PI

12.solveKIIt1=bwith a direct sparse solver

13.computet2=KBIt1

15.//solve interface equations

17.//calculate internal degrees of freedom

18.b=PI-KIBxB

19.solveKIIxI=bwith a direct sparse solver

3.3.3 界面方程的并行求解

在本文中，并行預條件共軛梯度算法被用于求解界面方程。考慮到求解界面方程的過程中往往需要大量的迭代，將計算數據在CPU和MIC設備間頻繁地傳輸會嚴重影響系統通信效率，因此本文在求解界面方程時只利用CPU端的多核處理器完成相關的計算工作。對于計算量較小的任務，如向量與向量的和或點積操作，直接交給每個MPI進程完成；而對于計算量較大的任務，如矩陣和向量的乘積，在每個MPI進程內部將進一步派生出許多線程，這樣可以有效發揮CPU端的多核資源優勢，從而縮短計算時間。

4 數值算例

在廣州超算“天河二號”計算機上開展了一系列測試，以驗證本文所設計的并行計算方法的有效性。該超級計算機的主要配置請參閱2.2節。并行計算每次啟動的結點機總數依次為50,100,150和200，實際參與并行計算的核數依次為9 750,19 500,29 250和39 000。

該異構超算上的每顆CPU有12個計算核心，每個MIC卡可使用57個計算核心，因此配置2顆CPU和3個MIC卡的結點內共有195個計算核心。“天河二號”CPU/MIC單結點內的協同并行加速測試表明，1個MIC卡的計算性能約相當于2路CPU的性能，即1CPU+1MIC可以達到1CPU計算速度的3倍。

4.1 并行計算數據準備

在并行計算開始前先進行兩級分區處理。對于每個計算模型，一級剖分時產生的一級子區域總數應與每次參與并行計算的結點機總數相等，即應依次為50，100，150，200。由于“天河二號”每個結點機內配置5個計算設備(2顆CPU、3個MIC)，因此二級剖分時每個一級子區域被進一步剖分為5個二級子區域。

4.2 結構靜力分析并行計算數值算例

結構靜力分析并行計算的測試模型如圖3所示，主要分析土木領域某盾構隧道模型在土體重力作用下的初始應變和應力。隧道主線長2 792 m，土層長2 000 m,寬300 m，深80 m。網格剖分采用四面體單元，該模型具有120 574 032個單元，23 379 547個節點。

結構靜力分析并行計算的結果數據如表1所示。表1中，并行計算總時間從讀取數據文件開始到計算各子區域應變/應力結束；縮聚時間指的是消去各子區域內部自由度所花費的時間、迭代求解時間包括求解界面方程的時間、回代內部自由度的時間和計算子區域應變和應力的時間。

Figure 3 Tunnel model for parallel computing圖3 隧道計算模型

由表1可見，隨著計算核數的增加，有限元結構靜力分析并行計算的總時間依次減少，加速比呈現逐步增加的趨勢。另據本文相關研究測試數據表明，與傳統迭代計算方法相比，在使用19 500核時，本文方法總計算耗時只有前者的67%，并且總計算時間節省優勢隨著核數的增加不斷擴大。

Table 1 Statistics of time and performance of parallel computing for static analysis表1 靜力學并行計算時間和性能統計

由于本文設計的多層次多粒度協同并行計算方法不但實現了CPU設備和MIC設備間的負載均衡，而且顯著節省了系統所需的通信開銷，故能有效利用異構眾核超級計算機的硬件資源獲取較高的加速比和并行計算效率。

4.3 結構動力分析并行計算數值算例

結構動力分析并行計算的測試模型如圖4所示，主要分析核電領域某防浪堤模型在地震載荷作用下的動力響應行為。防浪堤長2 280 m，土層長1 820 m，寬900 m，深350 m。網格剖分依然采用四面體單元，該模型具有113 492 512個單元，21 669 412個節點。動力分析的時間步長定為0.01 s，求解10 s，總計1 000個時間步。

結構動力分析并行計算的結果如表2所示。表2中，并行計算總時間從程序啟動開始到所有時間步都計算完畢結束；縮聚時間指的是消去各子區域內部自由度所花費的時間、迭代求解是所有時間步的迭代求解時間的總和，每個時間步的迭代求解時間包括求解界面方程的時間、回代內部自由度的時間和計算子區域應變/應力的時間。

Figure 4 Breakwater model for parallel computing圖4 防浪堤計算模型

Table 2 Statistics of time and performance of parallel computing for dynamic analysis表2 動力學并行計算時間和性能統計

由表2可見，利用多層次多粒度協同并行計算方法求解結構動力分析問題時，系統同樣能夠獲取較高的加速比和并行計算效率。另據本文相關研究測試數據表明，與傳統迭代計算方法相比，在使用19 500核時，本文方法總計算耗時只有前者的63%，并且總計算時間節省優勢隨著核數的增加不斷擴大。

相對靜力學分析而言，載荷均衡和最小化系統間的通信開銷對動力學分析顯得更加重要。這是由于動力學分析要對多個時間步進行求解，任何載荷分配的不均衡和通信效率的下降都會在多次求解中被放大，從而極大地影響到系統整體的并行計算效率。

結構動力分析并行計算需要對多個時間步依次求解，每個時間步都類似求解一次結構靜力學問題。因此，結構動力分析時間步越多，所累計節省的絕對計算時間也越多，工程問題計算分析的效率提升優勢也更為明顯。

本文設計的多層次多粒度并行計算方法很好地解決了CPU設備和MIC設備間的負載均衡問題，對于大部分計算任務實現了多核和眾核資源的協同并發執行。計算結果表明顯著降低了在超級計算機上執行計算任務時所產生的通信開銷，因而可以大大提高大規模結構動力分析的并行計算效率。

4.4 單機可擴展性測試

為評估核數一定時不同網格數量的計算規模對程序性能的影響，本文在單個結點機上進行了系統可擴展性測試。測試模型為左端固定右端受豎直載荷的懸臂梁，主要分析懸臂梁在靜力載荷作用下的應變和應力。懸臂梁采用四面體單元進行網格剖分，通過不斷加大計算網格的節點數和單元數測試系統在不同計算規模下的性能，得到的并行計算結果如表3所示。

Table 3 Scalability test on a single machine表3 單機可擴展性測試

由表3可見，隨著懸臂梁模型計算網格節點數和單元數的成倍增加，并行計算總時間也基本呈現出同樣成倍增加的趨勢，這表明本文方法在有限元規模擴大時具有良好的可擴展性。依據此趨勢，在研究其他復雜工程問題時，可以從小規模計算模型建模調試入手，在此基礎上繼續擴大計算模型網格規模、調用更多計算結點，以獲得更為準確的仿真計算結果。

5 結束語

為提高異構眾核超算上有限元結構分析大規模并行計算的效率，設計了一種多層次多粒度協同并行計算方法。通過使每個計算作業被分解映射到異構超算的各硬件層面運行，在有效處理CPU與MIC之間的負載均衡問題的基礎上顯著降低了異構超算的通信成本。因此，它能夠通過充分發揮超算系統的硬件資源規模優勢來表現出最優的計算性能。

采用本文方法在“天河二號”超算系統上開展了幾個有限元結構靜力分析和動力分析所對應的的大規模并行計算數值算例的測試驗證，啟動的CPU+MIC核數總計達39 000個，分析的有限元規模超過1億單元。測試結果表明，本文方法表現出良好的加速比和并行計算效率。未來計劃將在新一代國產異構超級計算機上進一步開展測試和工程應用研究。