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

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

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

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

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

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

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

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

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

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

2 SW26010處理器的PANDAS眾核并行優(yōu)化方法

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

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

2.1 區(qū)域分解和稀疏矩陣壓縮存儲

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

圖2 模型區(qū)域分解及數(shù)據(jù)交換[6]

2.2 系數(shù)矩陣選擇性分塊預處理方法

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

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

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

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

圖4 從核任務劃分圖

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

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

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

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

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

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

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

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

3.2 多孔介質壓縮分析

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

表2 測試模型參數(shù)信息

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

圖5 z方向位移場分布

3.3 斷層系統(tǒng)運動變形分析

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

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

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

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

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

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

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

4 結論

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