基于SW26010處理器的PANDAS眾核并行優化方法及在地質變形分析中的應用

王雪純,邢會林,戴黎明,郭志偉,劉駿標

(1.中國海洋大學 海底科學與探測技術教育部重點實驗室,山東 青島 266100;2.青島海洋科學與技術國家實驗室 深海多學科交叉研究中心,山東 青島 266237;3. 中國海洋大學 海底科學與工程計算國際中心,山東 青島 266100)

地質體變形分析是地球科學問題研究的重要組成部分,包括拉伸、壓縮、剪切等體積和形狀的改變,復雜變形可以表示為多種簡單變形的疊加。有限元方法作為地質體變形分析最重要的方法之一,被廣泛應用于地震預測、油氣資源勘探與開發、土壤與巖石力學行為模擬等諸多領域。宋孝天等[1]利用擴展有限元法計算裂紋擴展路徑,研究了裂隙面變形參數對非貫通裂隙巖體壓縮特性的影響。陳云鍋[2]通過構建青藏高原內部不同震例震黏性有限元模型來分析震后形變動力學過程,定量約束青藏高原殼幔流變屬性以及震后余滑分布和演化。隨著模型規模的增大,單機系統難以滿足算力和存儲的要求,利用并行計算技術和超級計算機系統進行大規模數值模擬計算已成為解決問題的重要手段。為進一步提高計算速度、降低功耗,異構眾核處理器逐漸取代同構多核處理器成為超級計算機系統的首選。但異構眾核處理器上的并行程序設計需要考慮模塊劃分、數據劃分和線程同步等問題,代碼實現相對復雜,已有的同構計算機并行算法無法直接應用于異構眾核架構的計算機系統。

SW26010是我國自主研發的高性能異構眾核處理器,支持消息傳遞接口(message passing interface,MPI)、線程加速庫(Sunway Athread庫)等并行編程模型,能夠提供比傳統高性能平臺更強大的計算能力和訪存帶寬。王澤松[3]基于申威眾核架構特點實現了有限元區域分解方法在超級計算機平臺的基于信息傳遞接口和加速線程庫的高效并行求解策略;傅游等[4]基于“神威·太湖之光”高性能計算平臺的Open ACC編程模型實現了Tend_lin的并行優化;柳安軍等[5]利用共享內存和寄存器通信等技術實現了Palabos軟件在新一代神威超算上兩相流算法的百萬核心規模的并行計算。

并行自適應非線性變形分析軟件(parallel adaptive nonlinear deformation analysis software,PANDAS)是一個基于有限元方法和格子玻爾茲曼方法(lattice boltzmann method, LBM)的數值模擬軟件,已成功應用于含裂隙的非均質多孔材料模擬中,用于解決應力變形、應力破壞、流體流動、熱傳導、化學反應等多物理場高度非線性耦合的問題[6]。目前地質體變形分析領域的眾核并行計算軟件仍存在較大空缺,而傳統的多CPU并行方法不適用于異構眾核處理器,亟需針對異構眾核處理器的體系結構特點研發地質體變形計算的重要算子,即有限元計算的并行方法。本研究基于SW26010異構眾核處理器,利用MPI模型和Sunway Athread庫對PANDAS軟件進行并行優化。通過區域分解、數據壓縮存儲、系數矩陣預處理和并行共軛梯度法求解方程組等方法提升軟件模擬速度,并將眾核并行程序應用于全地球運動模型、多孔介質壓縮變形和接觸體的接觸狀態改變等常見地質變形分析,為復雜地質體的變形模擬分析提供技術支持。

1 SW26010異構眾核處理器及其并行模式

SW26010異構眾核處理器由4個核組組成,每個核組包括1個內存控制器、1個運算控制核心(management processing element,MPE,又稱主核)、64個從核組成的客戶終端設備(customer premise equipment,CPE)、系統管理接口[7],如圖1所示。SW26010單芯片雙精度浮點計算峰值達3.168 TFLOPS。從核通過內存控制器訪問內存并且支持直接內存訪問方式(direct memory access,DMA),每個從核具有一個大小為64 kB的局部數據存儲器(local data memory,LDM),通過DMA方式完成主存向從核LDM空間的數據傳輸,能夠更好地利用帶寬,加快數據傳輸速度。但是從核LDM空間大小有限,無法一次性滿足大規模計算的存儲要求,頻繁啟動DMA將導致程序性能下降,因此需要在每次傳輸時盡可能多的傳輸數據。

圖1 SW26010處理器單核組結構示意圖

SW26010處理器使用MPI+X并行模式完成大規模并行計算,該模式下主核調用MPI接口實現跨核組通信、調用Sunway Athread庫實現從核并行計算以及從核資源控制。主核調用從核資源分為三個過程:從核資源初始化、從核執行數據運算、回收從核資源,分別由主核調用athread_init()、athread_spawn()、athread_halt()三組函數實現。Sunway Athread庫提供了線程級加速方法。

2 SW26010處理器的PANDAS眾核并行優化方法

PANDAS計算主要包括生成單元剛度矩陣、總剛度矩陣組裝以及方程組求解三個過程。求解方程組時采用共軛梯度法經過多次迭代獲得近似值,該算法無多級嵌套循環,易于實現并行,但每次迭代都依賴上次迭代結果,因此迭代過程無法并發。

本研究采取的并行策略為:簡單的邏輯控制和只需少量迭代的步驟仍串行執行,大量迭代的操作并行執行。首先將串行代碼移植到SW26010主核上,然后利用區域分解技術、矩陣壓縮存儲技術、系數矩陣預處理和并行化共軛梯度算法實現PANDAS的眾核并行優化。 通過區域分解將大規模網格數據分解到多個核心上并行計算,調用MPI接口實現數據交換,使用壓縮存儲技術降低內存消耗,利用眾核處理器的主從核異構架構優化共軛梯度算法,加快方程組求解速度,實現更高性能的并行計算。

2.1 區域分解和稀疏矩陣壓縮存儲

區域分解以節點為單位,將大的計算區域分解成多個規模較小的子區域用于并行計算[8],按照選定的軸向和數量劃分區域并記錄各個子區域的通信表,子區域之間節點不重復,劃分時不考慮接觸情況,分割后更新通信表并形成各個子區域的邊界條件和材料屬性等信息。區域劃分后,節點分為三類:邊界節點、內部節點和外部節點。計算過程中,內部節點與外部節點之間需要進行數據交換(圖2),同時邊界節點所在單元被記錄到構成該單元節點的所有分區中。每個主核處理一個子區域,主核間通過MPI_ISEND接口發送數據和MPI_IRECV接口接收數據,實現不同子區域間的數據交互,為分解生成的n個區域分配n+1個(編號0～n)進程共同完成計算。其中,0～(n-1)號進程作為計算器執行運算任務,生成區域剛度矩陣并通過迭代求解器計算推導出的線性方程組;n號進程作為控制器,收集來自其他進程的向量并進行接觸搜索,運算器和控制器二者結合將原問題求解轉為子區域并行求解。

圖2 模型區域分解及數據交換[6]

2.2 系數矩陣選擇性分塊預處理方法

程序中通過定義接觸面組合表示每對接觸面,其中包含一個主接觸面和一個從接觸面。使用罰因子方法處理接觸關系,通過兩個點在接觸面法線方向上的相對位移大小判斷是否接觸。發生接觸時需要滿足法向接觸力F=f·n=E·g,其中E為法向上的罰因子,g為外法線穿透量,f為接觸力,n為接觸面外法線方向[9]。首先通過全局搜索找到可能發生接觸的從節點和主接觸面,然后通過內外算法[10]計算從節點是否存在主接觸面來判斷潛在接觸對間的兩個面有無接觸。接觸問題被歸結為附帶約束條件的線性問題,系數矩陣的條件數通常很大,這將直接影響迭代過程的收斂速度。使用預處理技術會增加部分計算開銷但卻可以加快迭代求解收斂速度,因此通常在求解前對大規模系數矩陣進行預處理。預處理采用選擇性分塊方法,從總剛度矩陣對角線位置的非零節點出發橫向搜索,將被罰因子λ約束的同一接觸組中的節點放在一個以該節點為起點的塊中。若節點之間相鄰則仍保持相鄰位置關系,否則按照搜索順序排列,如圖3所示。重排序后生成的稀疏矩陣的非零項集中在以對角線為中心的帶狀區域并對新矩陣執行LU分解。

圖3 選擇性分塊預處理方式

2.3 SW26010處理器上的并行共軛梯度算法

PANDAS中的大規模迭代求解器基于一般方程Ax=b的共軛梯度法實現。主核調用從核函數對方程組進行求解,具體的迭代求解步驟由從核完成。主核將計算任務平均分配給從核,使用Sunway Athread庫athread_get_tid()接口獲取從核編號,按照編號分配對應區段的數據(圖4),未被分配的數據由主核完成計算。每個主核啟動從核進行主要迭代操作,通過循環將需要計算的數據依次批量讀入從核私有的LDM存儲空間,減少直接存取主存空間的次數,從核對讀入的數據運算后,將處理完的數據批量寫回主存空間。按照串行共軛梯度算法的計算流程編寫并行算法,迭代過程中涉及到的矩陣運算、數組遍歷等大規模的循環操作由每個主核的從核陣列并行完成,主核間通過MPI接口進行數據交換以實現區域之間的數據傳遞。

圖4 從核任務劃分圖

對于每個主核計算的部分方程組,用系數矩陣的行號除以m+1取余數,根據余數分配到對應編號的從核。從核執行向量的乘加運算,所有從核的行為相似,僅處理的數據不同。將共軛梯度法的迭代過程按照數據關聯關系分解為多個從核函數,根據每個函數需要傳遞的參數個數計算每個從核函數每次可以傳輸的最大數據量,從核調用athread_get()和athread_put()函數每次按照最大數據量讀入LDM空間和寫回內存,依次循環直至所有數據計算完畢,釋放從核資源。主核將計算得到的局部解向量集合獲得最終結果。

綜上所述,本研究提出的眾核并行優化方法為:通過區域分解將大模型分解為多個規模較小的模型分配給主核并行計算,降低單個核心計算數據量;利用壓縮存儲技術將原本的大規模稀疏矩陣壓縮為三個一維矩陣,節省存儲空間;通過對系數矩陣的預處理加速迭代求解過程的收斂速度,減少迭代次數、縮短迭代時間;利用SW26010從核陣列實現大規模方程組的迭代求解,將串行共軛梯度算法分解為多個從核函數執行,主核之間通過MPI實現數據在不同區域之間的傳遞。最終0號主核集合來自其他主核的解向量合并輸出,得到方程組的解向量。

3 PANDAS眾核優化方法在地質體變形分析中的應用及性能測評

為了測試上述方法的效果,將PANDAS模擬程序分別應用于SW26010異構眾核處理器(1.54 GHz)的單個核組和Intel Xeon Platinum 8268處理器(2.9 GHz),并對優化并行后PANDAS模擬結果的可靠性和速度提升效果進行測試。

3.1 全地球模型速度場模擬分析

全地球模型從內到外由內核、外核、地幔(下地幔)、轉換帶、軟流圈、巖石圈地幔和地殼組成,其中轉換帶、軟流圈、巖石圈地幔和地殼的深度界面根據LITHO1.0[11]得出。對板塊邊界的網格進行加密處理,將其離散為11 768 300個節點、10 072 140個單元和38個非連續接觸面。板塊邊界每層的屬性取自PREM模型[12],材料參數信息如表1所示。根據NNR-MORVEL56模型[13]中各板塊的歐拉極點角速度對地球表面節點施加速度邊界條件。

表1 地球內部圈層劃分及物理性質

通過區域分解將模型劃分為64個子區域,使用眾核并行PANDAS對模型模擬計算。墨卡托投影得到地球表面的運動速度場。經過對比,計算獲得的板塊速度與實際GPS結果[14]基本一致。實驗結果表明,眾核并行計算的模擬結果與真實情況基本吻合,該方法能夠應用于復雜變形分析。

3.2 多孔介質壓縮分析

巖石中分布著硬度不同的物質,這些物質受到拉伸或者壓縮時產生不同的變形,導致巖石在壓縮過程中位移發生不均勻變化。對取自塔里木盆地的碳酸巖樣本(分辨率為9 μm)的CT掃描結果進行二值化處理[15],設計用于測試的網格模型。模型材料參數見表2,生成節點數為16 974 593、單元數為16 777 216的八節點六面體網格模型。模型底部固定,從頂部施加z軸負方向的壓力。

表2 測試模型參數信息

用眾核并行的PANDAS模擬巖石模型在壓縮過程中的變形情況,將模型按照坐標遞歸二分法在x、y方向四等分,在z方向二等分,共劃分為32個區域。由圖5可以看出,巖石中存在硬度不同的兩種材料,在壓縮時其內部位移發生不均勻變化。經計算,優化后的眾核并行程序速度提升8.1倍。

圖5 z方向位移場分布

3.3 斷層系統運動變形分析

基于巖石夾層摩擦實驗模型[6]對網格進行精細劃分,生成用于測試的千萬網格三維斷層系統模型。采用速率相關摩擦定律[9]描述沿接觸面的非線性摩擦接觸問題。假設材料的性質相同,密度ρ=2 600 kg/m3,楊氏模量E=2 070 MPa,泊松比γ=0.2。模型的節點數為11 054 164、單元數11 033 920,生成兩個接觸體網格G1、G2并設置一組接觸面。邊界條件設置為:第一階段模型左側物體G1沿著x方向加載,G1加載面上所有節點沿x方向固定,右側物體G2沿y方向移動;第二階段G2從相對靜止狀態以1.0 mm/s的速度沿著y軸負方向移動。

使用眾核并行PANDAS程序模擬兩個接觸體的運動狀態變化,在x、z方向二等分、y方向四等分將模型劃分為16個區域,使用選擇性分塊預處理方法對系數矩陣進行變換。y方向的速度變化結果如圖6所示。從圖6(a)可以看出,在第二階段加載剛開始時,y方向的速度場連續分布,兩個物體之間未發生明顯滑動;隨著加載的進行,如圖6(b)所示,G1接觸面上方區域首先由閉鎖狀態進入局部滑動狀態,接觸面出現局部的速度不連續;圖6(c)表明,隨著加載持續進行,兩個接觸體逐漸開始滑動脫離,接觸面兩側發生較大的速度突變。

圖6 不同場景下y方向速度場

設置收斂條件滿足精度高于1.0×10-10,迭代過程中的相對殘差數據如圖7所示?？梢钥闯?使用選擇性分塊預處理方法時,迭代求解過程中的殘差數據基本呈線性分布,收斂效果良好。對模型進行不同數量的區域分解并測試(表3),在最優情況下眾核并行程序收斂一次所需的時長為CPU串行程序的11.6%,速度提升7.6倍,表明眾核程序在分析接觸問題時具有明顯的加速效果。

表3 計算時長對比(眾核并行/CPU單核串行)

圖7 選擇性分塊預處理方法殘差曲線

上述模型的測試結果表明,眾核并行優化后PANDAS程序的計算能力明顯提升。眾核程序能夠正確分析大規模的復雜地質體變形,對大規模單一地質體變形和復雜接觸體變形的模擬均獲得較好加速效果,可大幅縮短模擬時間。在對斷層系統模型的測試中,當使用64個核心并行計算時,加速比開始降低,可能由于當前規模大小的模型,隨著并行核心數量的增加,進程間的通信、等待時間占比增大造成的,說明本程序仍有進一步優化空間。

4 結論

本研究針對PANDAS進行千萬級大規模模型模擬時求解方程組耗時長的問題,從區域分解、矩陣壓縮存儲和優化迭代求解幾個方面,基于SW26010處理器對求解方程組相關部分核心代碼進行并行優化,并采用全地球運動模型、多孔介質壓縮模型和斷層系統模型對優化后的程序進行驗證。相較于傳統單一多線程或多進程的并行方法,本研究的眾核并行方法結合線程、進程兩種并行方式,充分發揮不同處理器的優勢,提高了計算速度。全地球模型速度場模擬結果與實際GPS數據基本吻合,優化方法具有可行性;多孔介質壓縮模型測試數據表明,眾核并行程序在分析單個物體變形問題時計算速度提升了8.1倍;斷層系統模擬測試數據表明,優化后的程序在分析多個物體接觸問題時速度提升了7.6倍,獲得良好的加速效果。以上結果表明,本研究方法可在地質體變形分析的大規模計算中發揮積極作用。然而實際的地質變形問題涉及的網格模型更加復雜,接觸面數量也更多,仍需在現有研究基礎上設計適用不同模型特征方程組的并行求解算法和區域分解方法。