面向GPU的非結構網格有限體積計算流體力學的圖染色方法優化*

張 曦，孫 旭，郭曉虎，杜云飛，盧宇彤，劉 楊

(1. 中山大學 計算機學院(軟件學院)， 廣東 廣州 510006；2. 哈璀國家超算中心 達斯伯里實驗室， 英國 沃林頓 WA4 4AD; 3. 中國空氣動力研究與發展中心， 四川 綿陽 621000)

非結構網格有限體積方法在航空航天、海洋環境等許多科學計算領域得到了廣泛的應用，形成了一些優秀的數值模擬軟件，例如OpenFOAM[1]、FUN3D[2]、Fludity[3]、NNW-PHengLEI(風雷)[4]等。隨著高性能計算技術的發展，科學計算進入“E級”時代，數值模擬軟件的計算性能也不斷提高。特別是通用圖形加速器(general purpose graphics processing unit, GPGPU)的快速發展，計算能力強、功耗低的GPU在科學計算中起到了越來越重要的作用。如何充分利用GPU，提高非結構網格有限體積數值模擬的計算性能，已經成為高性能計算領域的重要課題之一。

GPU性能的發揮主要取決于硬件占用率和訪存延遲。非結構網格的數據存儲是不規則的，需要通過幾何拓撲關系訪問數據，例如通過面循環訪問體變量，這將導致頻繁的不規則內存訪問，增大訪存延遲。因此，如何利用GPU多級存儲體系減小訪存延遲是提高應用性能的關鍵。

幾何拓撲依賴讀寫模式在多線程并行計算中，可能造成多個線程改變相同的全局內存或共享內存，導致錯誤的計算結果，即資源競爭問題。圖染色[5]是解決資源競爭問題的一種重要方法。該方法將發生沖突的面通過染色分成若干組，以保證每一組中的任何兩個面都不會更改相同的體變量。多線程并行計算可以在每一組中正確地實現。

圖染色方法在可移植性和計算性能方面具有優勢。圖染色方法通過常規的編程語言即可實現，幾乎可以移植到任何一種硬件平臺。一些學者在有限元、有限體積計算流體力學(computational fluid dynamics, CFD)的實際案例中發現圖染色方法具有較好的性能[6-7]。FUN3D[8]和SU2[9]也采用了圖染色方法解決多線程計算中的資源競爭問題，但都集中于CPU。本文討論了圖染色方法對于GPU計算性能的影響。

圖染色方法在分組之后，每一組的計算仍然存在不規則訪存的問題，需要進一步優化。Giuliani等[10]針對簡單的幾何拓撲關系，提出了一種優化方法，但無法將其應用在未知量定義在體中心的計算。Sulyok等[11]也提出了優化方法，但是只針對一種網格尺寸進行了性能分析。本研究的圖染色方法可應用在多種規模的復雜三維非結構網格。

1 非結構網格有限體積CFD數學模型

本研究采用的非結構網格有限體積CFD程序是NNW-PHengLEI(風雷)。該程序由中國空氣動力研究與發展中心研發。

1.1 控制方程與離散

控制方程為積分形式的NS(Navier-Stokes)方程，表示為

(1)

其中：Fc和Fv分別為對流通量和黏性通量；Q包含5個未知數，可以表示為

(2)

式中，ρ表示流體密度，u、v、w分別表示速度在x、y、z三個方向的分量，e表示內能。 這些物理量定義在控制體Ω中心上。

有限體積方法將物理空間劃分為許多小的控制體(volume)，體的邊界稱為面(face)。相應地， NS方程的離散形式可以表示為

(3)

式中：Ωvol表示控制體的體積；控制體擁有NF個面，每個面的面積為Sm；對流通量Fc采用Roe格式計算，并且采用限制器抑制激波產生的數值振蕩；黏性通量Fv采用中心差分格式計算；Splart Allmaras方程湍流模式用來計算湍流的影響；顯格式龍格庫塔方法用來離散時間項；最終通過對流通量和黏性通量的殘差來更新Q，保證計算的收斂。程序在一個時間步的流程如圖1所示。

圖1 程序在一個時間步的流程Fig.1 Flow chart in one time step

1.2 網格拓撲和數據結構

非結構網格幾何拓撲關系由體編號和面編號來表示。非結構網格通常采用面的編號順序來存儲體編號。例如將每個面對應的較小的體編號存儲在owner中，將較大的體編號存儲在neighbor中。如圖2所示，編號為O和P的2個體(O

圖2 幾何拓撲關系：體和面Fig.2 Geometry topology: volume cell and face

有限體積CFD程序中，變量ρ、u、v、w、e等存儲在體中心，稱為體數據；對流通量Fc和黏性通量Fv等變量則存儲在面中心，稱為面數據。體數據和面數據由各自的編號索引，并以結構體數組(structure of array，SOA)方式存儲。

非結構網格的幾何拓撲關系和數據存儲方式使得體數據可以由體編號直接索引，也可以由面編號間接索引。例如表示殘差的體變量res，其編號O的數據可以表示為res[O]或者res[owner[f2]]。

1.3 資源競爭

在NNW-PHengLEI(風雷)程序中，絕大多數的體數據以1.2節中提到的間接索引方式，通過面循環實現對體數據的修改，從而導致多個不同的面對同一個體數據進行修改。在GPU計算中，多個線程對應的面很可能同時對相同的體數據進行修改操作，導致錯誤的計算結果。這類計算需要首先解決資源競爭問題，才能在GPU上得到正確的計算結果。本文統計了需要解決資源競爭問題的計算占比，如表1所示，資源競爭問題的計算時間占總計算時間的50%以上，極大地影響了計算性能。圖染色方法可以有效地解決資源競爭問題，需研究其性能及其優化。本文從NNW-PHengLEI(風雷)程序中抽取了兩種典型的算法討論。

表1 資源競爭問題計算時間占總時間的比例

算法1描述了由通量flux(面數據)累加計算殘差res(體數據)的過程，屬于非定常項計算中的通量累加更新殘差計算部分(如圖1、表1中的C6所示)。該算法通過遍歷所有的面(0～numFaces-1)，將面數據flux累加到owner、neighbor間接索引的體變量res。結合圖2可見，GPU的多個線程在處理不同的面(f1、f2、f3、f4)時，有可能同時修改res[O]，造成資源競爭。

算法1 通量累加計算

算法2描述了局部最大壓力plmax(體數據)的計算過程，屬于對流通量計算中的限制器計算部分(如圖1、表1中C5所示)。該算法通過遍歷所有的面，借助owner、neighbor間接索引plmax、ps，完成比較。由圖2可知，plmax[O]有可能被GPU處理f1、f2、f3、f4面的多個線程同時修改，造成資源競爭。

算法2 局部最大壓力計算

2 圖染色算法及其應用

算法1和算法2在GPU計算中的資源競爭問題可以通過圖染色方法解決。首先，將存在沖突的面進行染色分組，保證每一組中的任何兩個面對應的計算不會修改相同的體數據；其次，將每一組計算分別放在GPU上，并通過全局內存讀寫的優化，降低訪存延遲。

2.1 基于貪婪算法的面染色

面染色是一個NP完全問題[12]。在計算流體力學應用中，最優的染色結果只存在于簡單的二維幾何拓撲中[10]。對于復雜的三維幾何拓撲，僅能夠得到一般的染色結果[13]。本文采用貪婪算法對應用在空氣動力學中的復雜三維幾何拓撲進行面染色分組。

首先，計算面的沖突關系。位于同一個體中的面相互沖突，計算面的沖突關系就是收集與該面同體的其他面。計算過程中除了利用owner、neighbor，還要利用表征體-面拓撲的變量volFaces、volFacesPosi。volFaces存儲了每個體包含的面數量以及面編號；volFacesPosi存儲了每個體的拓撲信息在volFaces中的位置。如圖2所示，編號為O的體包含編號為f1、f2、f3、f4的4個面。O對應的拓撲信息存儲在volFaces中，起始位置為volFacesPosi[O]。O包含的面數量可以表示為volFaces[volFacesPosi[O]]=4，對應的面編號可以表示為volFaces[volFacesPosi[O]+m]=fm(m=1,2,3,4)。因此，在包含f2的編號為owner[f2]=O的體中，存在與f2沖突的面f1、f3、f4；類似地，編號neighbor[f2]=P的體包含的面分別為f5、f2、f6、f7，其中f5、f6、f7與f2沖突。最終，收集f2的沖突關系，包括f1、f3、f4、f5、f6、f7。匯集所有面的沖突關系，存儲在變量nfcf、faceConflict、faceConflictPosi中，分別對應沖突面個數、沖突面ID、沖突面的起始位置。f2的沖突信息表示為nfcf[f2]=6，faceConflict[faceConflictPosi[f2]+m] (m=0,1,2,3,4,5)分別對應了f1、f3、f4、f5、f6、f7。具體描述見算法3。

算法3 面的沖突關系計算

接著，利用沖突關系，通過基于貪婪算法的面染色方法將所有面進行染色，避免有沖突關系的面顏色相同，具體描述見算法4。對于編號faceID的面，通過faceConflict[faceConflictPosi[faceID]+m] (m=0, …,nfcf[faceID]-1)遍歷所有與faceID沖突的面，根據這些面已經染的顏色判斷出面faceID可以染色的最小值。染色的結果存儲在變量numFaceGroups和colorFaces中。numFaceGroups表示顏色的種類；colorFaces存儲了每個面的顏色。

算法4 基于貪婪算法的面染色

最后，將不同顏色的面分組。利用面染色結果numFaceGroups和colorFaces，將相同顏色的面連續排列，存儲在變量faceGroup中。同時記錄每種顏色對應的面數量groupNum和每種顏色的面在變量faceGroup中的起始位置colorPosi。

2.2 圖染色分組計算的GPU算法

在面染色分組之后，每個組可以分別在GPU上計算。為了提高數據的局部性，減少數據讀寫的內存訪問延遲，分別采用GPU共享內存、體和面編號重排、面數據重排的策略優化圖染色方法。

算法5實現了通量累加(算法1)的GPU計算。通過染色分組數量numFaceGroups確定循環次數，將每個分組依次放入GPU進行計算，避免了資源競爭。借助面染色分組信息groupNum、colorPosi，可以得到每個分組對應的面數量numFacesInGroup和分組在faceGroup的起始位置colorStart。接著，GPU線程通過faceGroup得到面編號faceID，再借助owner和neighbor得到相應的體編號ownVolID和ngbVolID，將面上的通量flux[faceID]累加到相應的體數據res[ownVolID]和res[ngbVolID]中。

算法5 通量累加的圖染色算法

算法6實現了局部最大壓力(算法2)的GPU計算。該算法借助染色分組信息，將分組依次放入GPU計算，避免了資源競爭。GPU線程編號faceInGroupID通過faceGroup和colorStart，實現了與面編號faceID的映射，并通過owner和neighbor，得到相應的體編號ownVolID和ngbVolID。最終，實現了體數據ps[ownVolID]與plmax[ngbVolID]的比較以及ps[ngbVolID]與plmax[ownVolID]的比較。

算法6 局部最大壓力計算的圖染色算法

在算法5和算法6中，GPU線程根據面編號faceID，通過owner和neighbor間接索引到編號分別為ownVolID和ngbVolID的體數據。ownVolID和ngbVolID很難同時保證體數據讀寫的對齊。此外，染色分組造成了分組中面數據的非對齊訪問。因此，加劇了GPU全局內存訪問的延遲。

采用共享內存優化圖染色法：針對間接索引帶來的體數據無法對齊讀取的問題，利用共享內存減小體數據讀取延遲，相應的描述在算法7、算法8中。

算法7采用GPU共享內存對通量累加的圖染色算法(算法5)進行優化，以減小數據讀取延遲。該算法首先在GPU中定義共享內存變量fluxShare。接著，GPU線程先將面數據flux讀取到fluxShare中，再將fluxShare累加到體變量res中。

算法7 共享內存優化圖染色算法(通量累加)

算法8利用GPU共享內存優化局部最值的圖染色算法(算法6)，以減小數據讀取延遲。該算法先將體變量ps讀取到共享內存psShare中，再利用psShare和體數據plmax進行比較。

算法8 共享內存優化圖染色算法(局部最值)

針對體數據的非對齊讀寫，可以通過重排體編號、面編號的策略[14]進行優化。重排后的編號使得面數量最大(即numFaceGroups最大)的分組實現體數據的對齊讀寫。通過提高數據的局部性，以減少GPU全局內存的訪問延遲，相應算法描述在算法9、算法10中。

算法9 體、面編號重排優化圖染色算法(通量累加)

算法10 體、面編號重排優化圖染色算法(局部最值)

算法9對通量累加的圖染色法GPU計算(算法5)進行了優化。與算法5相比，對owner、neighbor、faceGroup進行了重排。重排后的ownerRe、neighborRe、faceGroupRe可以保證體編號ownVolID和ngbVolID在計算量最大分組中的連續性，使得體數據res的讀寫對齊。

算法10對局部最值的圖染色法GPU計算(算法6)進行了優化。通過體、面編號排，實現了面數量最多分組中體數據plmax、ps讀寫的對齊。

最后，針對面數據的非對齊讀取，重排面數據，使得面數據按照染色分組后的面編號重新排列，實現了面數據的對齊讀取。相應的描述見算法11。

算法11對通量累加的染色法GPU計算(算法5)進行了優化。面數據flux是按照面標號遞增的順序排列的，但是經過染色分組后，分組中的面編號排序(faceGroup)并不是按順序的，造成了面數據的非對齊讀寫。按照faceGroup中的面編號順序重排flux，得到面數據fluxRe，將fluxRe累加到相應的體數據res中。

算法11 面數據重排優化圖染色算法(通量累加)

3 實驗結果與分析

3.1 算例和計算環境

為了分析圖染色及其優化算法的性能，采用ONERA M6外流場CFD計算的三維非結構網格展開性能測試。如圖3所示，以正四面體和棱柱體網格填充三維計算域，棱柱體網格存在于M6機翼附近。并采用如表2所示的5種不同網格規模測試算法的性能。

圖3 ONERA M6網格Fig.3 ONERA M6 mesh

表2 5種網格規模

算法5～11通過CUDA C實現。原CFD程序可以進行單、雙兩種精度的計算，因此本文分別針對雙精度操作數和單精度操作數進行測試。每個kernel都重復運行1 000次并記錄運行時間。

測試平臺包括CUDA 10.0驅動的Nvidia Tesla V100 GPU以及CUDA 8.0驅動的Nvidia Tesla K80 GPU。

3.2 面染色結果

基于貪婪算法的面染色結果如表3所示，可見所有網格都被分為了8組。混合網格中含有棱柱體，棱柱體網格具有5個面，因此，面染色結果至少存在5個分組。對于復雜的三維幾何拓撲，比理想化的分組數量多2～3個分組是可以接受的[15]。

表3 面染色分組

3.3 算例運行時間

算法5～11對應的GPU kernel在V100上的運行時間如表4和表5所示。

表4 運行時間(V100，雙精度)

表5 運行時間(V100，單精度)

算法5～11對應的GPU kernel在K80上的運行時間如表6和表7所示。

表6 運行時間(K80，雙精度)

表7 運行時間(K80，單精度)

3.4 圖染色方法在通量累加計算的性能比較

通過面染色分組，圖染色方法可以消除GPU計算中的資源競爭。在面染色分組后，為了提高數據在GPU內存中的局部性，降低內存讀寫的延遲。針對累加操作，本研究采用共享內存、體編號和面編號重排、面數據重排對圖染色方法進行了優化。下面進一步探討這些優化方法的效果。

為了比較各種圖染色方法在通量累加操作中的性能，以未優化的算法5運行時間為分母，得到算法5(未優化)、算法7(共享內存)、算法9(體、面編號重排)、算法11(面數據重排)的運行時間比值。

圖4所示為V100的比較結果：共享內存的優化效果不明顯，性能提高在5%以內；體編號和面編號重排使得執行時間增大；只有面數據flux重排的優化方法(算法11)有效地減少了運行時間。雙精度和單精度的計算結果近似，相比未加優化的圖染色方法，面數據重排可以帶來20%左右的性能提高。

(a) 雙精度 (a) Double precision

圖5所示為K80的比較結果，和V100類似，只有面數據重排起到了優化效果：雙精度計算的性能提升在15%左右，單精度計算的性能提升在10%左右。

(a) 雙精度(a) Double precision

針對體編號、面編號重排造成計算性能下降的結果作了進一步分析。以網格4(392萬個體、807萬個面)對應的8個染色組為例，分別在未優化(算法5)和體、面編號重排(算法9)2種算法下計時，結果如表8所示。

表8 網格4分組運行時間(V100，單精度)

可見，經過體編號、面編號重排后，第1組的運行時間由0.640 s降低到0.241 s，這與第一組的面數據完全實現對齊相吻合。但是，第2～7組的運行時間全部增加，致使體、面編號重排后的總運行時間超過了未優化算法的運行時間。由此可見，體、面編號重排使得第1組數據的局部性得到提高的同時，惡化了其他幾組的數據局部性。

3.5 圖染色方法在局部最大值計算的性能比較

NNW-PHengLEI(風雷)程序中存在一些局部最值操作，同樣面臨資源競爭問題。針對局部最大值計算，本研究采用共享內存、體編號和面編號重排這兩種方法對圖染色方法進行優化。

為了比較各種圖染色方法在局部壓力最大值計算中的性能，以未優化的算法6的運行時間為分母，得到算法6(未優化)、算法8(共享內存)、算法10(體、面編號重排)的運行時間比值。

圖6所示為V100的比較結果，共享內存優化方法的效果不明顯；體編號、面編號重排使得運行時間增加；未優化的圖染色法(算法6)性能最佳。圖7顯示的K80比較結果與V100結果一致。

(a) 雙精度(a) Double precision

(a) 雙精度(a) Double precision

造成這一結果的原因主要在于局部最值運算僅操作體數據。面循環使得體數據訪問的局部性降低。體、面編號重排或者共享內存都未能提高體數據的局部性。

3.6 風雷程序GPU計算整體性能

本研究已經完成了NNW-PHengLEI(風雷)程序非結構網格程序的異構開發和優化。圖8展示了風雷程序GPU版本在1塊GPU(Nvidia Tesla V100)和CPU版本在28個CPU核(Intel Gold Xeon 6132)并行計算的加速比(雙精度)。可見，在不同網格規模下，平均加速比達到19。

圖8 風雷程序1塊GPU計算和28 CPU核計算的加速比Fig.8 Speedup of 1 GPU compared with 28 CPU cores

4 結論

本文采用圖染色方法解決非結構網格有限體積計算流體力學在GPU計算中的資源競爭問題，并采用了共享內存、體編號和面編號重排、面數據重排的優化策略對圖染色方法進行優化。通過圖染色方法及其優化算法在V100、K80兩種GPU上，針對應用在空氣動力學中不同規模的復雜三維幾何拓撲網格，測試通量累加計算和局部最大壓力計算，得到如下結論：

1)對于通量累加計算，面數據重排可以將圖染色方法在GPU上的計算性能提高10%～20%；共享內存的優化作用小于5%；體編號和面編號重排在提高1個分組數據局部性的同時惡化了其他分組的數據局部性，造成計算性能下降。

2)對于局部最大壓力計算，共享內存和體、面編號重排都無法提高體數據的局部性。

今后，將針對不同的硬件架構，例如AMD GPU、ARM、RISC-V等展開圖染色方法進行研究。