數(shù)值堆熱工流體高精細(xì)并行模擬優(yōu)化技術(shù)研究

董玲玉，周志鋒，戴潮虎，趙 珂，吳宗蕓，劉天才，趙民富，楊 文，胡長軍,*

(1.北京科技大學(xué) 計算機(jī)與通信工程學(xué)院，北京 100083；2.中國原子能科學(xué)研究院 反應(yīng)堆工程技術(shù)研究所，北京 102413)

熱工流體高性能模擬是數(shù)值反應(yīng)堆的重要組成部分。熱工流體計算程序按建模尺度可劃分為系統(tǒng)程序、子通道程序和CFD程序，其中CFD程序能對反應(yīng)堆中復(fù)雜區(qū)域的三維流動過程(包括湍流)進(jìn)行高精細(xì)的模擬，其核心內(nèi)容是求解N-S(Navier-Stokes)方程組。求解N-S方程組的數(shù)值方法包括有限差分法、有限體積法、有限元法、譜方法和譜元法等。其中譜元法是Patera等[1]通過結(jié)合譜方法的高精度和有限元法的區(qū)域適應(yīng)性發(fā)展而來的，符合數(shù)值反應(yīng)堆對高分辨率和高效并行求解的追求。然而堆幾何的復(fù)雜性以及譜元法的高精度特點使得數(shù)值反應(yīng)堆模擬對計算資源、存儲資源的需求巨大。例如在中國實驗快堆中，燃料棒直徑約為繞絲直徑的6.3倍，二者的接觸面需要大量精細(xì)網(wǎng)格來描述。經(jīng)測試，在平均網(wǎng)格質(zhì)量指標(biāo)超過0.6的情況下，單根繞絲棒流道的網(wǎng)格數(shù)量至少需要100萬，快堆全堆的網(wǎng)格數(shù)量則會達(dá)到百億量級。再如Kothe等[2]估算采用大渦模擬的方法進(jìn)行堆芯局部區(qū)域的模擬需要P級的計算能力，進(jìn)行全堆芯的模擬則需要E級的計算能力。高精細(xì)、大規(guī)模數(shù)值模擬對數(shù)據(jù)存儲、區(qū)域分解、并行求解以及計算機(jī)算力等方面提出了巨大挑戰(zhàn)。

以我國神威·太湖之光和天河二號為代表的超級計算機(jī)，計算能力已居世界領(lǐng)先地位。國產(chǎn)超算多為混合架構(gòu)，例如主從核異構(gòu)的神威超算，其98%計算能力分布在從核中，因此想要充分發(fā)揮其計算能力必須提高從核利用率。然而目前的CFD軟件大多為純CPU開發(fā)，亟需開展針對國產(chǎn)異構(gòu)架構(gòu)超算的CFD軟件的優(yōu)化移植工作。

本文以基于譜元法求解N-S方程的數(shù)值方法為研究對象，瞄準(zhǔn)大規(guī)模網(wǎng)格區(qū)域分解時間成本高以及缺少面向國產(chǎn)超算體系結(jié)構(gòu)的優(yōu)化移植方法兩個關(guān)鍵問題，提出一種面向海量精細(xì)網(wǎng)格的混合并行遞歸譜二分法實現(xiàn)的大規(guī)模區(qū)域分解方法，并建立一套以小矩陣乘為核心的申威(SW26010處理器)眾核架構(gòu)并行優(yōu)化技術(shù)。

1 區(qū)域分解與CFD軟件優(yōu)化移植的國內(nèi)外相關(guān)研究

常用區(qū)域分解方法有：Wesseling[3]和Anders[4]提出的遞歸二分法、貪婪算法、Saule等[5]提出的啟發(fā)式算法以及Pothen等[6]根據(jù)遞歸二分法的思想提出的遞歸譜二分法。其中，遞歸譜二分法相對其他算法的劃分效果好，但當(dāng)網(wǎng)格規(guī)模過大時，使用該算法的時間成本不可忽略，針對該問題，Barnard等[7]設(shè)計了多層次遞歸譜二分法(MRSB)、并行多層次遞歸譜二分法(PMRSB)[8-9]，但方案的實現(xiàn)依賴于T3D內(nèi)存共享體系結(jié)構(gòu)且各分區(qū)元素不均衡。王琥等[10]在MRSB的基礎(chǔ)上，提出了多層次遞歸譜二分優(yōu)化算法(MMRSB)，相對MRSB改善了負(fù)載平衡和通信時間兩個方面，但優(yōu)化方案的引入增加了一定的時間成本。目前的研究中缺少在不改變遞歸譜二分法分區(qū)效果的前提下提高分區(qū)效率的算法。

當(dāng)前主流CFD軟件大多是基于MPI并行的純CPU架構(gòu)開發(fā)的，但面對異構(gòu)架構(gòu)的發(fā)展趨勢，涌現(xiàn)出了向不同類型的異構(gòu)架構(gòu)體系進(jìn)行優(yōu)化的研究。Fournier等[11]在IBM Blue Gene/P超算上將Code_Saturne擴(kuò)展到超3.2萬核。Ojha等[12]將基于有限體積法的OpenFOAM移植到采用CPU+Intel Xeon PHI加速卡混合架構(gòu)的Endeavor集群上。面向國內(nèi)超算平臺，孟德龍等[13]使用Athread加速線程庫，將OpenFOAM移植到神威·太湖之光超算上。李芳等[14]將基于有限差分法的OpenCFD移植到申威異構(gòu)眾核處理器，采用直接數(shù)值模擬的方法模擬5.52億網(wǎng)格，在213萬核上并行效率達(dá)到50%。但面向國產(chǎn)異構(gòu)架構(gòu)的優(yōu)化大部分圍繞有限差分、有限體積法為基礎(chǔ)的CFD軟件進(jìn)行開展，并沒有關(guān)于譜元法CFD模擬軟件的優(yōu)化。面向國產(chǎn)異構(gòu)架構(gòu)的移植主要在于發(fā)揮每個計算節(jié)點內(nèi)加速核的算力，常見方法是利用加速核計算程序熱點。傳統(tǒng)CFD中的FVM等方法的程序熱點多為大型稀疏矩陣向量乘(SPMV)，可面向國產(chǎn)異構(gòu)架構(gòu)依據(jù)SPMV經(jīng)典方法進(jìn)行優(yōu)化；而譜元法的計算熱點是小而稠密的矩陣-矩陣乘，由于矩陣規(guī)模較小，數(shù)學(xué)庫和編譯器難以提供最好的內(nèi)核性能，經(jīng)過測試后發(fā)現(xiàn)使用通用優(yōu)化技術(shù)(如OpenACC、OpenMP)或優(yōu)化庫(如BLAS、GEMM)反而會帶來大幅性能下降。

2 數(shù)值反應(yīng)堆熱工流體CFD模擬計算模型

2.1 計算模型

通常，數(shù)值反應(yīng)堆熱工流體CFD模擬的求解對象是不可壓縮N-S方程，如式(1)所示，上下兩個方程分別稱為動量方程和連續(xù)性方程，假設(shè)求解區(qū)域為一維空間Ω：

(1)

其中：u為速度；p為壓力；Re為雷諾數(shù)；f為力源項。N-S方程是一個偏微分方程，時間與空間上均存在偏導(dǎo)數(shù)，而計算機(jī)不能直接計算連續(xù)的偏導(dǎo)數(shù)，因此求解思路是將偏導(dǎo)數(shù)分別在時間與空間上離散為計算機(jī)可計算的格式，譜元法即是一種空間離散的數(shù)值方法。

式(1)應(yīng)用變分原理得到變分公式：

(2)

(3)

使用剛度矩陣、質(zhì)量矩陣、對流算子對式(3)進(jìn)行簡化：

(4)

簡化后得到聚合公式：

(5)

對式(5)中的時間偏導(dǎo)數(shù)使用k階向后差分法(BDFk)進(jìn)行離散得到：

(6)

然而向后差分法為隱式格式，相比顯式格式增加了計算量，因此采用k階外推法(EXTk)對Cun+1項進(jìn)行外推，得到：

(7)

其中，bk、ak分別為隱式時間導(dǎo)數(shù)離散化和顯式外推的系數(shù)，根據(jù)文獻(xiàn)[15]使用的數(shù)值技術(shù)，式(7)最終分裂成兩個部分：

(8)

式(8)即為基于譜元法求解N-S方程的核心計算式，第1項是壓力泊松方程，第2項是亥姆霍茲方程，兩個方程可簡化成如下形式：

(9)

該形式可看作是線性方程組Ax=b的求解。其中：φ為自變量；K為剛度矩陣；M為質(zhì)量矩陣；f為右端項的和；h1和h2為變量系數(shù)，為式(10)的使用條件[16]：

(10)

兩個方程均可使用Krylov子空間迭代法求解。壓力泊松方程的系數(shù)矩陣為非對稱矩陣，可使用廣義最小余量法(GMRES)；速度亥姆霍茲方程的系數(shù)矩陣為對稱矩陣，可使用預(yù)處理共軛梯度法(PCG)。對于求解較復(fù)雜的壓力泊松方程，迭代前可應(yīng)用一些預(yù)處理技術(shù)來加速迭代，例如粗網(wǎng)格技術(shù)。

2.2 基本流程

數(shù)值反應(yīng)堆熱工流體CFD模擬總體上可分為預(yù)處理、主求解與后處理3個步驟。圖1示出基于譜元法的CFD模擬基本流程。在預(yù)處理過程中，需定義流道的幾何形狀(幾何建模)、定義求解域、劃分求解域為譜元(網(wǎng)格生成)、分解區(qū)域、設(shè)置參數(shù)(如流體屬性、邊界條件)等；主求解部分需對N-S方程進(jìn)行時間離散、空間離散與算子分裂，如2.1節(jié)所述將離散后的方程分裂為壓力泊松方程與速度亥姆霍茲方程，隨后在每個譜元內(nèi)對這兩個方程進(jìn)行迭代求解；最后在后處理部分中輸出每個譜元得到的方程解，檢查輸出結(jié)果的正確性，通過適當(dāng)?shù)姆绞斤@示數(shù)據(jù)(如圖形、數(shù)據(jù)對比)。

圖1 熱工流體CFD模擬基本流程Fig.1 Basic flow of thermal fluid CFD simulation

基于譜元法的CFD模擬流程中存在兩個關(guān)鍵點：1) 面對海量精細(xì)網(wǎng)格，區(qū)域分解的時間成本將不可忽略，例如單根繞絲棒流道區(qū)域分解時間為11 h，而并行計算1 s的模擬結(jié)果僅需5 min，為加速程序計算，節(jié)約前處理時間，非常有必要對區(qū)域分解進(jìn)行并行化處理；2) 基于譜元法中張量積的特殊性質(zhì)[17]，小而稠密的矩陣-矩陣乘法是計算核心，以典型基于譜元法的軟件Nek5000為例，矩陣-矩陣乘函數(shù)占整個程序運行時間的60%，因此針對矩陣-矩陣乘的并行優(yōu)化是有意義的。

3 混合并行大規(guī)模區(qū)域分解方法

區(qū)域分解將網(wǎng)格劃分后的元素進(jìn)行聚合，得到粗粒度任務(wù)，是基于譜元法進(jìn)行并行求解的必要過程(圖2)，其最終目的是使各分區(qū)內(nèi)的網(wǎng)格數(shù)量相近以滿足均衡的計算負(fù)載，以及使各分區(qū)之間的共享邊界較少以滿足較低的通信頻率。而找到同時滿足這兩個條件的算法是一個NPC問題[18]，即無法在多項式時間內(nèi)找到最優(yōu)解，遞歸譜二分法是對該問題較好的近似。

圖2 基于譜元法的并行求解算法必要流程Fig.2 Essential process of parallel solution algorithm based on spectral element method

3.1 遞歸譜二分法原理

遞歸譜二分法根據(jù)網(wǎng)格對偶圖的拉普拉斯矩陣對應(yīng)的特征值來進(jìn)行分割，拉普拉斯矩陣為圖的度矩陣與鄰接矩陣的差值。首先使用Lanczos算法計算拉普拉斯矩陣的等效三對角矩陣，根據(jù)三對角矩陣計算費德勒向量，然后對費德勒向量頂點進(jìn)行排序，給每個子域分配近似1/2的頂點，得到兩個子域，在每個子域上遞歸該過程直到分區(qū)數(shù)達(dá)到設(shè)定的值。費德勒向量的含義是圖的拉普拉斯矩陣的次小特征值λ2對應(yīng)的特征向量X2，F(xiàn)iedler[19-20]證明了X2的兩個特殊性質(zhì)。

性質(zhì)1：由費德勒向量劃分的子圖是連通的。子圖連通可一定程度上減小劃分的割數(shù)，如圖3a為子圖連通情況，其割為1；圖3b子圖不連通，其割為2。

a——子圖連通，割為1；b——子圖不連通，割為2圖3 費德勒向量性質(zhì)1Fig.3 The first property of Federer vectors

性質(zhì)2：費德勒向量分量的值近似1/2為正，1/2為負(fù)。該性質(zhì)可基本保持子圖之間的負(fù)載均衡，如圖4a為1對偶圖A，根據(jù)A的拉普拉斯矩陣L(A)可得到圖4c的費德勒向量xA，該向量中每個分量對應(yīng)于圖中1個節(jié)點，分量為正與分量為負(fù)的節(jié)點分別劃分到兩個子圖即得圖4d中的分區(qū)結(jié)果。

a——對偶圖A；b——A的拉普拉斯矩陣L(A)；c——L(A)的費德勒向量xA；d——由xA分量符號得到的分區(qū)結(jié)果圖4 費德勒向量性質(zhì)2Fig.4 The second property of Federer vectors

費德勒向量的兩個特殊性質(zhì)在圖的最小割與負(fù)載均衡之間做了平衡，因此遞歸譜二分法是一種質(zhì)量相對高的區(qū)域分解算法。

3.2 遞歸譜二分法MPI/OpenMP混合并行實現(xiàn)

基于遞歸譜二分法的網(wǎng)格區(qū)域分解適合于數(shù)據(jù)并行，但各分區(qū)操作同一全局?jǐn)?shù)組的過程中產(chǎn)生了數(shù)據(jù)相關(guān)。根據(jù)該特點，設(shè)計并行算法使用MPI通信消除子域之間的數(shù)據(jù)相關(guān)性，算法依功能劃分為3部分：單次圖劃分、進(jìn)程擴(kuò)展、結(jié)果回收。

單次圖劃分的作用是劃分得到兩個子圖，通過多次調(diào)用可完成單個進(jìn)程內(nèi)的圖劃分任務(wù)，劃分過程見3.1節(jié)。圖5示出使用4個進(jìn)程并行地將12個網(wǎng)格劃分到8個子域，其中進(jìn)程0、2、4、6分別進(jìn)行了3、2、1、1次單次圖劃分過程。

圖5 MPI并行遞歸譜二分法流程Fig.5 MPI parallel recursive spectrum bisection algorithm process

進(jìn)程擴(kuò)展的作用是將單次圖劃分得到的子圖分散到其他進(jìn)程，實現(xiàn)并行劃分。程序的開始僅有0號進(jìn)程進(jìn)行單次圖劃分，得到兩個子圖后，保留1個子圖，將相關(guān)的全局?jǐn)?shù)組通過MPI發(fā)送給編號為2depth+Rank(depth為當(dāng)前圖劃分深度，Rank為當(dāng)前進(jìn)程在通信域內(nèi)的編號)的進(jìn)程，隨后啟動進(jìn)程Rank與2depth+Rank開始劃分子圖。圖5第1次擴(kuò)展過程中，0號進(jìn)程擴(kuò)展了2號進(jìn)程。所有擴(kuò)展啟動的進(jìn)程將具有與主進(jìn)程相同的進(jìn)程擴(kuò)展方式，即并行規(guī)模隨著子圖劃分層次深入而增大，活動進(jìn)程的數(shù)量為2的冪次。

結(jié)果回收的作用是將分散到其他進(jìn)程的子圖劃分結(jié)果回收整合，得到整體的劃分結(jié)果。由于同一個進(jìn)程可能多次向編號不同的進(jìn)程發(fā)送數(shù)據(jù)進(jìn)行擴(kuò)展，數(shù)據(jù)的接收也需要根據(jù)發(fā)送的順序反向依次執(zhí)行，因此借助堆棧實現(xiàn)結(jié)果回收：在發(fā)送的過程中逐個將目標(biāo)進(jìn)程的進(jìn)程號壓入棧中，保證回收過程中子圖的劃分層次從下至上產(chǎn)生合并。除0號進(jìn)程與最底層進(jìn)程外，其他進(jìn)程在結(jié)果回收過程中需完成兩個操作：收集下層子圖劃分結(jié)果與向上級進(jìn)程返回結(jié)果。

除粗粒度MPI并行設(shè)計外，對程序計算熱點使用OpenMP進(jìn)行細(xì)粒度并行，可較好地利用兩種模型的優(yōu)點，使程序得到更好的性能。使用性能測試工具gprof對基于譜二分法的Nek5000的區(qū)域分解模塊進(jìn)行執(zhí)行概要分析后，發(fā)現(xiàn)稀疏矩陣向量乘函數(shù)ax占用了程序總體運行時間的46%，且函數(shù)本身執(zhí)行時間較高，因此可并行化函數(shù)ax來加速迭代。函數(shù)ax的調(diào)用主要發(fā)生在Lanczos計算過程中，計算相似三對角矩陣時每次迭代均會出現(xiàn)原對稱矩陣與正交列向量的乘法操作。由于ax處理的是非常稀疏而行內(nèi)元素分布均勻的矩陣，實際非零元素數(shù)量很低，可采用OpenMP共享內(nèi)存多線程并行。

串行稀疏矩陣向量乘函數(shù)ax中矩陣存儲使用CSR(compressed sparse row)格式。為節(jié)省內(nèi)存，CSR通過3個數(shù)組va、ja、ia來存儲矩陣中的非零元素。其中va數(shù)組存儲非零元素的值，ja數(shù)組存儲非零元素列下標(biāo)，ia數(shù)組最后1位存儲非零元素總個數(shù)，其他位是每行第1個非零元素在va數(shù)組中的索引?；贑SR格式的稀疏矩陣向量乘法由兩重循環(huán)構(gòu)成，外層循環(huán)按行遍歷，內(nèi)層循環(huán)讀取ia數(shù)組中的下標(biāo)范圍后按下標(biāo)范圍進(jìn)行遍歷，循環(huán)內(nèi)操作矩陣中的非零元素與對應(yīng)向量中的元素相乘及累加，最終按行編號順序逐步得到結(jié)果向量。

并行的思路在于對循環(huán)進(jìn)行拆解，由于區(qū)域分解問題中處理的稀疏矩陣每一行非零元素數(shù)量對應(yīng)網(wǎng)格相連的其他網(wǎng)格數(shù)量，從網(wǎng)格模型生成角度看這一數(shù)值大致相近，因此拆解外層循環(huán)并給各線程分配相同行數(shù)可保證基本的負(fù)載均衡。圖6為OpenMP任務(wù)劃分示意圖，其中avg=總行數(shù)/線程數(shù)量，即平均每個線程負(fù)責(zé)計算的行數(shù)。另外隨著圖劃分層次的加深，稀疏矩陣向量乘的規(guī)模不斷下降，計算耗時不足以掩蓋線程開銷，并行后的ax函數(shù)反而比串行函數(shù)消耗時間更大，因此設(shè)置了合理的閾值，舍棄小于該閾值的計算加速以在整體上獲得更好的性能。

圖6 OpenMP任務(wù)劃分示意圖Fig.6 OpenMP task division diagram

3.3 并行效果測試

混合并行大規(guī)模區(qū)域分解方法在天河二號超級計算機(jī)上進(jìn)行一系列實驗測試，包括使用MPI并行算法與Nek5000串行遞歸譜二分法的運行時間對比測試、MPI并行算法與MPI+OpenMP混合并行算法的運行時間對比測試。

圖7a示出Nek5000的區(qū)域分解模塊與僅經(jīng)過MPI優(yōu)化的并行譜二分程序運行時間對比。選用從20萬網(wǎng)格到1 000萬網(wǎng)格的8種網(wǎng)格規(guī)模(此處網(wǎng)格指譜元而非配置點)，在8進(jìn)程(8個物理核)運行的條件下，MPI優(yōu)化的并行譜二分程序相對于Nek5000的網(wǎng)格劃分模塊運行速度有明顯提升。劃分1 000萬網(wǎng)格Nek5000需超過11 h，而MPI優(yōu)化后的程序僅需55 min。為測試并行效率，在1 000萬網(wǎng)格下分別使用8進(jìn)程、16進(jìn)程、32進(jìn)程進(jìn)行測試，加速比分別為10.1、16.5、24.9，并行效率分別為1.26、1.03、0.77。并行效率不斷下降的原因是進(jìn)程擴(kuò)展必須在最初的幾次圖劃分結(jié)束后才能完成，而最初的圖劃分矩陣規(guī)模也是最大的，耗時占總時間較大比重，這部分消耗的時間無法隨著進(jìn)程數(shù)的增加而減小。測試過程中出現(xiàn)并行效率超過1的情況是開啟編譯器優(yōu)化帶來的額外性能提升。

在圖7a的基礎(chǔ)上，增加MPI并行算法與MPI+OpenMP混合并行算法的運行時間對比，如圖7b所示。混合并行算法在同樣8進(jìn)程運行條件下，每個進(jìn)程分配12個線程，共96個物理核心，可發(fā)現(xiàn)在原有MPI并行基礎(chǔ)上提升了30%～50%運行效率。劃分1 000萬網(wǎng)格混合并行程序?qū)r間進(jìn)一步縮短到34 min，相較于Nek5000的網(wǎng)格劃分模塊性能提升約95%。但此時并行效率僅0.21，這是由于OpenMP并行優(yōu)化算法存在適用閾值，當(dāng)圖劃分層次的加深到規(guī)模較小的矩陣時停止使用OpenMP，因此處理器占用率實際較低。

a——串行與MPI并行遞歸譜二分程序的運行時間；b——MPI優(yōu)化與混合優(yōu)化的運行時間圖7 不同網(wǎng)格規(guī)模(譜元數(shù))運行時間對比(處理器Xeon E5-2692V2，單節(jié)點內(nèi)存64G)Fig.7 Running time comparison (processor Xeon E5-2692V2, single node memory 64G)

混合并行大規(guī)模區(qū)域分解方法在并行過程中僅利用了主核CPU的計算能力，因此可兼容各超算平臺，但由于上述提到的最初圖劃分矩陣規(guī)模最大以及矩陣規(guī)模隨圖劃分層次加深變小的問題，算法可擴(kuò)展性不高。對于計算熱點稀疏矩陣向量乘的優(yōu)化，即使將OpenMP方法替換為面向異構(gòu)眾核(如神威·太湖之光)或異構(gòu)加速器(如曙光計算云、天河二號)的優(yōu)化，也僅能在一定程度上提高資源利用率，算法可擴(kuò)展性不會得到明顯提升。這是由于主核與其他從核/加速器之間存在通信開銷，同樣需舍棄小于適用閾值的矩陣的計算加速，該問題目前難以避免。

4 面向申威眾核架構(gòu)的并行優(yōu)化技術(shù)

4.1 SW26010介紹

神威·太湖之光由40 960塊自主研發(fā)的SW26010異構(gòu)眾核處理器組成，每個處理器包含4個核組(core group, CG)，260個異構(gòu)核心[21]。如圖8所示，單個核組由1個控制核心(manage processing element, MPE)和運算核心(computing processing element, CPE)陣列(8×8網(wǎng)格結(jié)構(gòu))組成，每個CPE有16 KB的L1指令緩存和64 KB的用戶可編程片上存儲(scratch-pad memory, SPM)，MPE和CPE共同共享8 GB內(nèi)存地址空間。

圖8 SW26010的架構(gòu)Fig.8 Architecture of SW26010

4.2 小矩陣乘面向SW26010的并行優(yōu)化技術(shù)

基于譜元法求解N-S方程的計算核心是底層的矩陣-矩陣乘。以Nek5000為例，使用性能測試工具gprof對整體程序進(jìn)行執(zhí)行概要分析，結(jié)果顯示矩陣-矩陣乘子函數(shù)mxf24、mxf8、mxf10需要60%以上的執(zhí)行時間，而mxfN函數(shù)為Nek5000的計算熱點，其中N(N∈[1,24])代表譜元法的階數(shù)。在計算過程中，僅有以下3種矩陣乘法形式出現(xiàn)(Cn1×n3=An1×n2Bn2×n3)：1)n1=N,n2=N,n3=N；2)n1=N2,n2=N,n3=N；3)n1=N,n2=N,n3=N2。

出于對矩陣規(guī)模的考量，針對情況1在主核進(jìn)行MPE向量化，針對情況2、3考慮使用從核進(jìn)行CPE移植優(yōu)化。

1) MPE向量化

對于情況1，每個矩陣的規(guī)模過小，不適合在CPE上進(jìn)行優(yōu)化，因此在MPE內(nèi)部利用SIMD向量化和指令級并行。MPE和CPE提供256位的SIMD擴(kuò)展部件，所以在MPE或CPE中可同時計算8個整型(32位)或4個雙精度浮點數(shù)(64位)。SW26010處理器提供了直接調(diào)用的對界SIMD接口(simd_load/simd_store)，該接口要求標(biāo)準(zhǔn)類型變量為32字節(jié)對界，而Nek5000經(jīng)sw5編譯后標(biāo)準(zhǔn)類型變量按照16字節(jié)對界，直接使用該接口會導(dǎo)致內(nèi)存訪問異常，處理器產(chǎn)生冗余操作將不對界的數(shù)據(jù)拆分成標(biāo)準(zhǔn)類型數(shù)據(jù)后再進(jìn)行加載/存儲，大幅降低了向量化求解的性能。而使用不對界的SIMD接口(simd_loadu/simd_storeu)可避免引入冗余操作。

Nek5000中矩陣乘代碼由Fortran語言編寫。Fortran中的數(shù)組是按列優(yōu)先存儲的，矩陣B和C的訪問是連續(xù)的，而矩陣A的訪問是不連續(xù)的。因此在計算前將A轉(zhuǎn)置，以便通過提供對內(nèi)存的連續(xù)訪問來提高SIMD性能。另外，矩陣乘法的內(nèi)循環(huán)次數(shù)固定為N，可通過手動展開N個內(nèi)部循環(huán)來激活指令流水線。

2) CPE移植優(yōu)化

情況2中的矩陣A和情況3中的矩陣B規(guī)模足夠大，可以并行處理。通過將矩陣塊的計算分配給M(M∈[1,64])個CPE以利用數(shù)據(jù)并行性，該過程需考慮數(shù)據(jù)的傳遞方式。SW26010在CPE的主存和LDM之間提供了兩個內(nèi)存訪問方式：gload/gstore和DMA。對于第1種方式，gload/gstore接口不考慮數(shù)據(jù)局部性，僅提供8 GB/s的物理帶寬。相比之下，DMA可通過64個CPE連續(xù)訪問大量數(shù)據(jù)來實現(xiàn)30.94 GB/s的有效帶寬。在使用DMA將數(shù)據(jù)移動到LDM之前，需考慮矩陣A如何移動。在文獻(xiàn)[22]的工作中，DMA和CPE的能力未被充分考慮，所以矩陣A在MPE中被轉(zhuǎn)置后才被傳遞到計算中，本文分兩種情況提出新的矩陣傳遞方式。

對于情況2，矩陣A的規(guī)模較大，因此將A劃分為多個數(shù)據(jù)塊并均勻分布到M個CPE上，每個CPE負(fù)責(zé)A的多行以實現(xiàn)負(fù)載均衡。SW26010的DMA有兩種訪問方式：連續(xù)訪問方式與跨步訪問方式。根據(jù)Jiang等[23]的總結(jié)報告，當(dāng)數(shù)據(jù)大于512 B時，跨步訪問方式的帶寬與連續(xù)訪問基本一致。由于Fortran是列優(yōu)先存儲的，而矩陣A劃分后的子數(shù)據(jù)塊基本大于512 B，因此每個CPE使用跨步訪問的DMA來讀取A的子數(shù)據(jù)塊。數(shù)據(jù)傳遞到CPE的LDM后，考慮到計算時能連續(xù)讀取數(shù)據(jù)塊，將傳遞來的A的子數(shù)據(jù)塊進(jìn)行轉(zhuǎn)置以提高性能。矩陣B足夠小，可直接傳遞給每個CPE，并且在計算過程中滿足連續(xù)的內(nèi)存訪問。計算完成后，矩陣C通過跨步訪問DMA方式返回到原地址，傳遞過程如圖9a所示。對于情況3，矩陣B可均分為多個列，這些列通過DMA以連續(xù)內(nèi)存訪問的方式傳遞到每個CPE中。矩陣A加載到LDM中，計算完成后，矩陣C通過連續(xù)訪問DMA方式返回其原始地址，傳遞過程如圖9b所示。

a——情況2的矩陣-矩陣乘優(yōu)化；b——情況3的矩陣-矩陣乘優(yōu)化圖9 數(shù)據(jù)傳遞過程Fig.9 Data transfer process

SW26010的CPE中同樣存在SIMD接口，并且傳遞到LDM的矩陣滿足double類型變量32 B對界，因此從核中可使用對界的SIMD接口來加速計算。對于mxfN(N%4=0)，需在傳輸前對矩陣A進(jìn)行轉(zhuǎn)置，以便連續(xù)加載到矢量寄存器中，提高SIMD性能。

4.3 性能測試

為驗證小矩陣乘優(yōu)化后軟件整體的性能，在神威·太湖之光超算上進(jìn)行單處理器、多處理器性能優(yōu)化測試以及強(qiáng)擴(kuò)展性和弱擴(kuò)展性測試。

1)性能優(yōu)化測試

首先測試了在單處理器中，不同多項式階數(shù)(譜元階數(shù))和譜單元個數(shù)下小矩陣乘優(yōu)化的性能特征。表1列出優(yōu)化前Nek5000在MPE中運行的結(jié)果。從表1可知，多項式階數(shù)(N)和每個CG分配的譜單元個數(shù)(E)均會影響程序執(zhí)行時間。隨著N和E的增加，程序的執(zhí)行時間也會逐漸增加。表2列出經(jīng)過矩陣-矩陣乘優(yōu)化后的執(zhí)行結(jié)果。結(jié)合表1數(shù)據(jù)可得到圖10。圖10a為N固定為12時，性能隨著E增加的提升幅度為6%～8%，提升并不明顯，圖10b為E固定為16時，性能隨著N增加的提升幅度為6%～47%，提升較為明顯。因此矩陣-矩陣乘優(yōu)化方法對多項式階數(shù)更為敏感。另外由圖10b可知，當(dāng)N=10或N=8時，矩陣-矩陣乘優(yōu)化后的性能相比優(yōu)化之前的版本反而有所下降。這是由于該方法依賴矩陣規(guī)模，對于規(guī)模過小的矩陣將其切分并傳輸?shù)紺PE進(jìn)行計算的時間遠(yuǎn)大于直接在MPE上計算的時間，因此該方法適合N≥12的情況。

a——每個CG中不同譜單元個數(shù)對性能的影響；b——多項式階數(shù)對性能的影響圖10 性能提升百分比Fig.10 Percentage performance improvement

表1 優(yōu)化前Nek5000在MPE中執(zhí)行時間Table 1 Execution time of Nek5000 in MPE before optimization

表2 矩陣-矩陣乘優(yōu)化方法的執(zhí)行時間Table 2 Execution time of matrix-matrix multiplication optimization method

圖11a示出單處理器中多種譜元階數(shù)(N)下不同的CPE數(shù)對程序性能的影響，其中譜單元個數(shù)(E)固定為16。圖中橫虛線代表了優(yōu)化前Nek5000在MPE中的執(zhí)行時間，觀察可知當(dāng)開啟的CPE數(shù)量太少時(如4個)，性能比優(yōu)化前反而有所下降；隨著CPE個數(shù)的增加，性能優(yōu)化逐漸達(dá)到瓶頸，因此最佳選擇是使用32個CPE。

a——不同CPE數(shù)量下多種N的執(zhí)行時間；b——不同CG下優(yōu)化前后的執(zhí)行時間圖11 運行時間對比(處理器SW26010，單節(jié)點內(nèi)存32 G)Fig.11 Running time comparison (processor SW26010, single node memory 32 G)

進(jìn)一步使用多處理器對優(yōu)化前與優(yōu)化后的性能進(jìn)行測試。圖11b示出在64個CG(每個CG使用32個CPE)、128個CG、256個CG、512個CG下的性能對比。其中譜單元的個數(shù)是1 872，多項式階數(shù)為24，在每次迭代過程中，壓力步和速度步各迭代50次，共運行150步。相比未優(yōu)化版本，矩陣-矩陣乘優(yōu)化方法在64個CG下性能提升51.9%，隨著CG數(shù)的增加，每個CG中分配到的譜單元減少，擴(kuò)展到512個CG時仍維持49.7%的性能增長。

2) 擴(kuò)展性測試

圖12a示出優(yōu)化后程序的強(qiáng)可擴(kuò)展性。測試規(guī)模為621 610個譜單元，多項式階數(shù)為12，壓力步迭代10次，速度步迭代50次，共迭代次數(shù)為20。擴(kuò)展性的測試范圍為2 560個CG、5 120個CG、10 240個CG、20 480個CG?？蓮膱D中看到，隨著CG的增加，并行效率逐漸下降，擴(kuò)展到20 480個CG(1 331 200核)時并行效率仍有77%。圖12b示出優(yōu)化后程序的弱可擴(kuò)展性，并行效率以160 個CG為基準(zhǔn)，擴(kuò)展至5 120個CG，譜單元的范圍是42 444、84 888、169 776、339 552、679 104、1 358 208。對比強(qiáng)擴(kuò)展性來說，弱可擴(kuò)展性有更高的并行效率，擴(kuò)展到5 120個CG(332 800核)時并行效率可達(dá)92.4%。

a——強(qiáng)擴(kuò)展性測試；b——弱擴(kuò)展性測試圖12 擴(kuò)展性測試(處理器SW26010，單節(jié)點內(nèi)存32 G)Fig.12 Scalability testing (processor SW26010, single node memory 32 G)

5 結(jié)論及下一步工作

熱工流體模擬是數(shù)值反應(yīng)堆的重要組成部分，為實現(xiàn)高精細(xì)、大規(guī)模的熱工流體CFD模擬，本文提出了一種混合并行大規(guī)模區(qū)域分解方法，建立了一套以小矩陣乘為核心的面向申威眾核架構(gòu)的并行優(yōu)化技術(shù)，并分別在天河二號和神威·太湖之光超算上進(jìn)行了測試。實驗結(jié)果顯示，混合并行大規(guī)模區(qū)域分解方法在8進(jìn)程、12線程條件下，劃分同等規(guī)模網(wǎng)格性能較開源CFD軟件Nek5000的串行區(qū)域分解模塊提升約95%；面向申威的小矩陣乘并行優(yōu)化技術(shù)性能在64個CG下較Nek5000提高了51.9%。目前兩種技術(shù)均在中國數(shù)值反應(yīng)堆核心軟件CVR-PACA中得以應(yīng)用。

下一步，面對國際上高端信息技術(shù)創(chuàng)新和競爭的制高點——E級超算，本文中的兩種方法需要做適應(yīng)修改。E級計算系統(tǒng)總體并行度將達(dá)到千萬量級，極大規(guī)模并行所帶來的復(fù)雜性使應(yīng)用程序需能表示所有異構(gòu)層次的內(nèi)在并行性并充分利用多級內(nèi)存層次結(jié)構(gòu)?；旌喜⑿写笠?guī)模區(qū)域分解方法中，MPI+OpenMP的并行模式最大程度利用了主核CPU的計算能力，但未能充分發(fā)揮E級超算中輔助計算核的算力，因此下一步對計算熱點稀疏矩陣向量乘的優(yōu)化可將OpenMP方法替換為面向異構(gòu)眾核(神威E級超算)或異構(gòu)加速器(曙光E級超算、天河E級超算)的優(yōu)化；面向申威眾核架構(gòu)的并行優(yōu)化方法能發(fā)揮神威E級超算的異構(gòu)架構(gòu)優(yōu)勢，使用了矩陣轉(zhuǎn)置、跨步訪問等技巧實現(xiàn)內(nèi)存的連續(xù)訪問，但對于多級內(nèi)存層次的利用不夠充分，可利用L1 Cache(一級緩存)等緩存部件進(jìn)一步加速求解。除此之外，下一步計劃依據(jù)實際的數(shù)值反應(yīng)堆需求設(shè)計自主的基于譜元法的并行CFD軟件系統(tǒng)。