非結(jié)構(gòu)有限體積CFD計(jì)算的網(wǎng)格重排序優(yōu)化*

張 勇，張 曦，萬(wàn)云博，何先耀，趙 鐘，盧宇彤

(1.中國(guó)空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所，四川 綿陽(yáng) 621000;2.中山大學(xué)計(jì)算機(jī)學(xué)院，廣東 廣州 510006)

1 引言

計(jì)算流體力學(xué)CFD(Computational Fluid Dynamics)通過計(jì)算機(jī)數(shù)值求解流體力學(xué)控制方程以獲得流動(dòng)信息。隨著高性能計(jì)算HPC(High Performance Computing)硬件的快速發(fā)展，CFD在工業(yè)裝備建設(shè)中得到了越來越多的應(yīng)用。

CFD通過計(jì)算網(wǎng)格離散空間，早期主要采用結(jié)構(gòu)網(wǎng)格。隨著幾何外形復(fù)雜性和保真度的增加，結(jié)構(gòu)網(wǎng)格生成在整個(gè)模擬周期中占據(jù)了絕大部分時(shí)間。非結(jié)構(gòu)網(wǎng)格具有自動(dòng)化程度高的優(yōu)點(diǎn)，能夠更好地適應(yīng)諸如飛行器、旋轉(zhuǎn)機(jī)械等復(fù)雜幾何外形，在大氣海洋、航空航天等許多科學(xué)計(jì)算領(lǐng)域運(yùn)用廣泛。由此發(fā)展了眾多數(shù)值模擬軟件，例如，基于有限體積方法的 OpenFOAM(Open Field Operation And Manipulation)[1]、FUN3D(Fully Unstructured Navier-stokes three-Dimensional)[2]、SU2(Stanford University Unstructured code)[3]，基于有限元方法的Fludity[4]和MOOSE(Multiphysics Object Oriented Simulation Environment)[5]，基于通量修正方法的pyFR(Python based Flux Reconstruction)[6]等。然而，由于非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)存儲(chǔ)不規(guī)則，數(shù)據(jù)只能間接訪問[7]，增大了訪存延遲，對(duì)非結(jié)構(gòu)網(wǎng)格計(jì)算效率造成了不利影響。隨著高性能計(jì)算技術(shù)的發(fā)展，科學(xué)計(jì)算進(jìn)入“E級(jí)”時(shí)代，通用圖形加速器GPGPU(General Purpose Graphics Processing Unit)憑借其計(jì)算能力強(qiáng)、功耗低的優(yōu)點(diǎn)，在科學(xué)計(jì)算中得到了越來越多的應(yīng)用。GPU的訪存模式[8]放大了非結(jié)構(gòu)網(wǎng)格不規(guī)則訪存帶來的延遲，如何提高非結(jié)構(gòu)網(wǎng)格的數(shù)據(jù)局部性，充分利用CPU和GPU的計(jì)算性能，提高非結(jié)構(gòu)網(wǎng)格有限體積數(shù)值模擬的計(jì)算性能，已經(jīng)成為高性能計(jì)算領(lǐng)域的重要課題之一。

提高非結(jié)構(gòu)網(wǎng)格的數(shù)據(jù)局部性有許多方法。一些研究人員對(duì)非結(jié)構(gòu)網(wǎng)格的數(shù)據(jù)存儲(chǔ)布局進(jìn)行了研究[9]，但無論是AOS(Array Of Structure)還是SOA(Structure Of Array)的存儲(chǔ)布局都避免不了數(shù)據(jù)的間接訪存。Weng等[10]通過調(diào)整循環(huán)模式改進(jìn)數(shù)據(jù)局部性，但是該方法只優(yōu)化了部分計(jì)算熱點(diǎn)，無法實(shí)現(xiàn)全局性優(yōu)化。

網(wǎng)格重排序是一類能有效提高數(shù)據(jù)局部性的方法，通過對(duì)網(wǎng)格幾何拓?fù)淙珞w、面、點(diǎn)的編號(hào)重排，改變數(shù)據(jù)訪問順序，減少訪存延遲。Lohner[11]針對(duì)有限元方法提出了若干重排序的方法。Andrew等[12]通過一種線排序方法，重排體和面編號(hào)，提升了高階DG(Discontinuous Galerkin)方法計(jì)算Euler方程的性能。Lani等[13]在有限體積方法中使用了RCM(Reverse Cuthill-Mckee)方法，但是性能改善有限。目前，大多數(shù)文獻(xiàn)都是對(duì)網(wǎng)格重排序在CPU或GPU等單一類型硬件上的性能影響展開研究，缺乏對(duì)CPU、GPU等不同硬件類型的統(tǒng)一性研究。

本文根據(jù)非結(jié)構(gòu)網(wǎng)格有限體積CFD計(jì)算的特點(diǎn)，基于風(fēng)雷軟件提出了一種網(wǎng)格重排序策略，以優(yōu)化數(shù)據(jù)局部性。首先,簡(jiǎn)述了非結(jié)構(gòu)有限體積CFD的理論方法和所采用的計(jì)算程序；然后，給出了本文所提出的網(wǎng)格重排序方法,并采用典型算例分析了該方法在CPU和GPU上的數(shù)值計(jì)算結(jié)果。

2 理論基礎(chǔ)與計(jì)算程序

本文研究基于開源CFD軟件“風(fēng)雷(PHengLEI)”[14-18]。風(fēng)雷軟件是中國(guó)空氣動(dòng)力研究與發(fā)展中心研發(fā)的面向流體工程的混合CFD平臺(tái)，面向全國(guó)開源(https://osredm.com/)。風(fēng)雷軟件是全球首款能同時(shí)支持結(jié)構(gòu)網(wǎng)格(算法)和非結(jié)構(gòu)網(wǎng)格(算法)的CFD軟件，也是同時(shí)支持有限體積和有限差分方法的“混合”計(jì)算平臺(tái)。本節(jié)簡(jiǎn)要介紹本文研究涉及到的理論基礎(chǔ)，包括控制方程與數(shù)值方法、網(wǎng)格幾何拓?fù)洹?/p>

2.1 控制方程與數(shù)值方法

NS(Navier-Stokes)流動(dòng)控制方程積分形式可以表示如式(1)所示：

(1)

其中,Fc和Fv分別表示對(duì)流通量和粘性通量(下標(biāo)c代表對(duì)流相關(guān)項(xiàng)，下標(biāo)v代表粘性相關(guān)項(xiàng))；Q表示包含守恒形式的流動(dòng)矢量，可以表示如式(2)所示：

(2)

其中,ρ、u、v、w和E分別表示流體密度、速度在x、y、z這3個(gè)方向的分量及單位質(zhì)量總能量。采用格心式有限體積方法，將這些流場(chǎng)變量定義在控制體(Control Volume)單元Ω中心上。

有限體積方法將物理空間劃分為許多小的控制體(Volume)；體的邊界稱為面(Face)。相應(yīng)地,NS方程可以表示如式(3)所示：

(3)

其中，Ωvol表示控制體的體積；控制體擁有NF個(gè)面，第m個(gè)面的面積為|Sm|；對(duì)流通量Fc采用Roe格式計(jì)算，通過采用限制器抑制流場(chǎng)間斷產(chǎn)生的數(shù)值振蕩；粘性通量Fv采用中心格式計(jì)算；湍流模型采用Splart-Allmaras一方程模型；時(shí)間推進(jìn)采用Runge-Kutta顯式和LU-SGS(Lower Upper-Symmetric Gauss Seidel)隱式方法。程序的具體算法參見文獻(xiàn)[14-18]，在一個(gè)時(shí)間步的流程如圖 1所示。

Figure 1 Flow chart of program in one time step

2.2 網(wǎng)格幾何拓?fù)?/h3>

非結(jié)構(gòu)網(wǎng)格存儲(chǔ)沒有結(jié)構(gòu)網(wǎng)格那樣規(guī)律，其點(diǎn)、面、體(單元)編號(hào)都是無規(guī)律存儲(chǔ)。在非結(jié)構(gòu)網(wǎng)格幾何拓?fù)潢P(guān)系中，面和體之間的連接關(guān)系尤為重要，廣泛應(yīng)用于各類通量計(jì)算中，其具體關(guān)系通過面編號(hào)和體編號(hào)來表示。對(duì)于格心型有限體積方法，體(Volume)即是計(jì)算單元本身，下同。

對(duì)于面編號(hào)，在風(fēng)雷軟件中，每個(gè)面兩側(cè)有且僅有2個(gè)體，分別為左、右單元。左和右的定義是相對(duì)的，軟件中規(guī)定面法向所指向的體為右側(cè)體，另一側(cè)為左側(cè)體，即面法向從左體指向右體。為編程方便，將面分為邊界面和內(nèi)部面2種類型，數(shù)組中先存儲(chǔ)邊界面，再存儲(chǔ)內(nèi)部面。對(duì)于內(nèi)部面，其左、右單元都是內(nèi)部單元；對(duì)于邊界面，左單元為內(nèi)部單元，右單元是左單元以邊界面為對(duì)稱的鏡像虛擬單元。

面法向采用叉乘計(jì)算得到。如圖 2所示，A=P3-P1、B=P2-P1分別是2個(gè)矢量，二者構(gòu)成一個(gè)平面。根據(jù)線性代數(shù)，其單位法向的計(jì)算表達(dá)式如式(4)所示：

(4)

對(duì)于任意凸多邊形的法向量，可以將每條邊與面心連接形成若干子三角形，將各子三角形的法向量以各自面積加權(quán)得到多邊形的法向量。

Figure 2 Diagram of plane normal vector calculation

采用2個(gè)數(shù)組存儲(chǔ)左、右體編號(hào)。如圖 3所示，編號(hào)為O和P的2個(gè)體都包含編號(hào)為f2的面，若規(guī)定面的右單元體編號(hào)存儲(chǔ)在cellRight數(shù)組中，且面的另一側(cè)體編號(hào)存儲(chǔ)在cellLeft數(shù)組中，則可通過數(shù)組直接得到左、右體編號(hào)，即cellRight[f2]=P，cellLeft[f2]=O。

Figure 3 Geometry topology: volume and face

變量ρ、u、v、w、p等存儲(chǔ)在體中心，稱為體數(shù)據(jù)(Volume Data)；對(duì)流通量Fc和粘性通量Fv等變量存儲(chǔ)在面中心，稱為面數(shù)據(jù)(Face Data)。體數(shù)據(jù)和面數(shù)據(jù)由各自的編號(hào)索引，并以SOA方式存儲(chǔ)。體數(shù)據(jù)可以由體編號(hào)直接索引，也可以由面編號(hào)間接索引。例如，表示殘差的體數(shù)據(jù)Res，其體編號(hào)為O,可以通過體索引Res[O]或者面索引Res[cellLeft[f2]]2種方式得到。風(fēng)雷軟件中面、體在數(shù)組中存儲(chǔ)的順序分別如圖4和圖5所示。

Figure 4 Numbering order of face

Figure 5 Numbering order of volume

圖4和圖5中，nCell是內(nèi)部體總數(shù)，nTotalCell表示內(nèi)部體和虛擬體數(shù)量之和，nBoundFace表示邊界面的數(shù)量或者虛擬體的數(shù)量，nTotalFace表示內(nèi)場(chǎng)面和邊界面數(shù)量之和。體編號(hào)在[0,nCell)，虛擬體編號(hào)在[nCell,nCell+nBoundFace)，從內(nèi)部體編號(hào)尾部開始逐一編號(hào)，定義第iFace個(gè)邊界面所對(duì)應(yīng)的虛擬體編號(hào)為iFace+nCell。

3 網(wǎng)格重排序方法

網(wǎng)格重排序分為兩部分：體網(wǎng)格按照RCM(Reverse Cuthill-Mckee)方法進(jìn)行重排，面網(wǎng)格以網(wǎng)格體為單位重新排列。下面分別說明。

3.1 RCM體重排序

RCM方法[19]采用近鄰排序的方法，以保證網(wǎng)格編號(hào)在空間中相對(duì)集中分布，從而使得體數(shù)據(jù)存儲(chǔ)時(shí)在內(nèi)場(chǎng)空間更緊密。通過RCM方法，反映網(wǎng)格鄰接關(guān)系的矩陣帶寬會(huì)減小。如圖 6所示，比較了29萬(wàn)網(wǎng)格規(guī)模下體網(wǎng)格排序前、后的矩陣結(jié)構(gòu),可以發(fā)現(xiàn)體編號(hào)重排序使得鄰接矩陣非零項(xiàng)更加緊湊。

Figure 6 Comparison of matrix structures with and without volume reordering

3.2 面重排序

本文設(shè)計(jì)了一種面重排序算法，針對(duì)風(fēng)雷程序的內(nèi)部面，以體網(wǎng)格擁有的面為基本單位，按照體網(wǎng)格內(nèi)原始的編號(hào)順序，為內(nèi)部面重新排序。這樣可以緩解邊界面編號(hào)造成的內(nèi)部面編號(hào)不連續(xù)問題，提高邊界數(shù)據(jù)訪存效率。具體算法如算法1所述。

算法1面編號(hào)排序算法

輸入：體-面拓?fù)潢P(guān)系offsetCellFace、numFaceOfCell和cellFace。

輸出：排序前、后面編號(hào)映射關(guān)系mapFace,排序后的幾何拓?fù)潢P(guān)系。

第1步升序排列cellFace中每個(gè)體網(wǎng)格對(duì)應(yīng)的面編號(hào)；

第2步初始化數(shù)組mapFace={-1}；

第3步初始化面編號(hào)計(jì)數(shù)器labelFace=nBoundFace-1；

第4步遍歷體網(wǎng)格，重排面編號(hào)：

FORcellID=0 TOnTotalCell-1

獲取體網(wǎng)格cellID在cellFace的偏移量offset=offsetCellFace[cellID];

獲取體網(wǎng)格cellID在cellFace中存儲(chǔ)的面編號(hào)數(shù)量numFaces=numFaceOfCell[cellID];

遍歷組成體網(wǎng)格cellID的全部面：

FORfaceInCell=0 TOnumFaces-1

獲取面編號(hào)faceID=cellFace[offset+faceInCell];

IF(面faceID是內(nèi)部面)&&(面faceID映射關(guān)系仍未改變)THEN

編號(hào)計(jì)數(shù)器labelFace累加1;

更新faceID的映射關(guān)系:

mapFace[faceID]=labelFace;

ENDIF

ENDFOR

ENDFOR

第5步采用mapFace更新網(wǎng)格幾何拓?fù)潢P(guān)系。

該算法首先將cellFace中以體網(wǎng)格為單位存儲(chǔ)的面編號(hào)按照升序排列(第1步)；然后初始化映射數(shù)組mapFace和編號(hào)計(jì)數(shù)器labelFace(第2、3步)，編號(hào)計(jì)數(shù)器從nBoundFace-1開始，只對(duì)內(nèi)部面編號(hào)；以網(wǎng)格體為單位，通過2層嵌套循環(huán)遍歷cellFace中的所有面編號(hào)(第4步)；對(duì)于內(nèi)部面，更新編號(hào)計(jì)數(shù)器，并將該值賦給映射數(shù)組(第2層循環(huán))；網(wǎng)格拓?fù)湎嚓P(guān)的變量，例如cellFace等，根據(jù)映射數(shù)組更新(第5步)。

4 算例測(cè)試與分析

4.1 算例、計(jì)算環(huán)境和測(cè)試方法

為分析網(wǎng)格重排序?qū)Ψ墙Y(jié)構(gòu)網(wǎng)格有限體積CFD計(jì)算的結(jié)果影響和優(yōu)化效果，采用5套不同規(guī)模的ONERA M6機(jī)翼外流場(chǎng)三維非結(jié)構(gòu)網(wǎng)格作為測(cè)試集。如圖 7所示，ONERA M6機(jī)翼三維非結(jié)構(gòu)網(wǎng)格由三棱柱和四面體組成，其中，三棱柱分布在物面附近的邊界層區(qū)域，四面體用于填充剩余計(jì)算域。5套網(wǎng)格的單元數(shù)如表 1所示。

Figure 7 ONERA M6 Mesh

Table 1 5 different meshes

測(cè)試采用風(fēng)雷軟件非結(jié)構(gòu)求解器GPU并行開源版本，該版本程序可在紅山開源平臺(tái)申請(qǐng)下載(https://forge.osredm.com/projects/p68217053/PHengLEI/tree/Branch_TH)。采用的硬件環(huán)境為：基于CUDA 10.0[20]驅(qū)動(dòng)的NVIDIA Tesla V100 GPU，Intel Xeon E5-2692 CPU。

本文的GPU并行優(yōu)化的前提是網(wǎng)格重排序，理論上網(wǎng)格重排序只影響計(jì)算效率，不影響計(jì)算的收斂結(jié)果。對(duì)此，首先對(duì)比分析網(wǎng)格重排序?qū)︼@式計(jì)算和隱式計(jì)算的影響。

4.2 網(wǎng)格重排序?qū)τ?jì)算結(jié)果的影響

4.2.1 對(duì)隱式LU-SGS方法的影響

LU-SGS方法數(shù)據(jù)依賴性強(qiáng)，對(duì)計(jì)算順序較敏感，本節(jié)討論網(wǎng)格重排序?qū)﹄[式計(jì)算LU-SGS方法的影響。

圖 8給出了400萬(wàn)網(wǎng)格規(guī)模(網(wǎng)格3)下平均殘差averageRes(NS方程未知數(shù)在相鄰時(shí)間步二階范數(shù)的平均值)隨迭代步變化的曲線。可以看到，網(wǎng)格編號(hào)排序前后殘差收斂值趨于一致，且重排序的收斂效果優(yōu)于未排序時(shí)的。

Figure 8 Variantion curves of residual error before and after mesh reordering under 4 million grid scale

不同網(wǎng)格規(guī)模下1萬(wàn)步處的平均殘差比較如表 2所示。可見，隨著網(wǎng)格規(guī)模的不斷增大，重排序的計(jì)算平均殘差逐漸減小。與未排序相比，其加速收斂的效果更加明顯。

不同網(wǎng)格規(guī)模下，網(wǎng)格重排序?qū)諝鈩?dòng)力學(xué)系數(shù)的影響如表3所示。由表3可知，網(wǎng)格重排序?qū)ιο禂?shù)的影響在0.24%之內(nèi)，對(duì)阻力系數(shù)的影響在0.5%之內(nèi)，二者均在氣動(dòng)計(jì)算可接受的誤差之內(nèi)，可認(rèn)為重排序不影響計(jì)算結(jié)果。

Table 3 Relative error of aerodynamic coefficients at 10 000 steps for reordering and unordering

綜上，由于LU-SGS方法數(shù)據(jù)依賴性較強(qiáng)，體網(wǎng)格編號(hào)重排改變了計(jì)算順序，提高了數(shù)據(jù)局部性，在一定程度上改善了求解過程的收斂效果，但對(duì)隱式方法的計(jì)算結(jié)果無明顯影響。

4.2.2 對(duì)顯式Runge-Kutta方法的影響

本節(jié)討論網(wǎng)格重排序?qū)︼@式計(jì)算Multi-Stage Runge-Kutta方法的影響。

圖 9給出了40萬(wàn)網(wǎng)格規(guī)模下平均殘差隨迭代步變化的曲線。可以看出，6萬(wàn)步后殘差曲線收斂，未排序和重排序的平均殘差收斂曲線沒有明顯差異。從圖 10中未排序與重排序平均殘差的插值可以看出二者平均殘差是相同的。

Table 2 Comparison of residual values at 10 000 steps for reordering and unordering

Figure 9 Variantion curves of residual error before and after mesh reordering under 400 000 gride scale

Figure 10 Difference of average residual values before and after mesh reordering

可見，網(wǎng)格重排序?qū)︼@式計(jì)算的結(jié)果完全沒有影響，這是因?yàn)椋猴@式方法計(jì)算時(shí)，每個(gè)網(wǎng)格只和局部幾何拓?fù)溆嘘P(guān)，網(wǎng)格重排序雖然改變了體網(wǎng)格、面網(wǎng)格的編號(hào)順序，但是并沒有改變幾何拓?fù)浣Y(jié)構(gòu)，不會(huì)影響計(jì)算結(jié)果。

4.3 CFD程序性能分析

從4.2節(jié)分析可見，網(wǎng)格重排序不會(huì)影響計(jì)算結(jié)果，說明針對(duì)GPU并行算法改造是可行的。本節(jié)分別就CPU和GPU上的計(jì)算來分析影響計(jì)算效率的熱點(diǎn)。

對(duì)于顯式計(jì)算，在CPU計(jì)算時(shí)，發(fā)現(xiàn)耗時(shí)較大的前8個(gè)熱點(diǎn)F0～F7函數(shù)占據(jù)了總計(jì)算時(shí)間的近70%，如圖 11a所示。在GPU計(jì)算時(shí)，發(fā)現(xiàn)耗時(shí)較大的前8個(gè)熱點(diǎn)K0～K7也占據(jù)了總計(jì)算時(shí)間的70%。

對(duì)于隱式計(jì)算，在CPU計(jì)算時(shí)，前8個(gè)熱點(diǎn)函數(shù)FI0～FI7占據(jù)了運(yùn)行總時(shí)間的一半以上。

4.3.1 網(wǎng)格重排序?qū)τ贑PU計(jì)算性能的影響

在CPU上，網(wǎng)格重排序?qū)︼@式計(jì)算和隱式計(jì)算的運(yùn)行時(shí)間都會(huì)產(chǎn)生影響。本節(jié)分別討論網(wǎng)格重排序?qū)︼@式計(jì)算和隱式計(jì)算在不同網(wǎng)格規(guī)模下運(yùn)行時(shí)間的影響。在Intel Xeon E5-2692上記錄的運(yùn)行時(shí)間。

Figure 11 Hot spot functions and kernels of explicit computing and implicit computing on CPU and GPU

表4是網(wǎng)格3(400萬(wàn)網(wǎng)格)規(guī)模下網(wǎng)格排序前、后熱點(diǎn)函數(shù)的運(yùn)行時(shí)間。以未排序的函數(shù)運(yùn)行時(shí)間作為分母，對(duì)運(yùn)行時(shí)間無因次化,結(jié)果如圖 12所示。可見，大部分函數(shù)的運(yùn)行時(shí)間在重排序的網(wǎng)格下縮短了10%～40%。

Table 4 Comparison of hot functions’ explicit executing time before and after mesh reordering

Figure 12 Comparison of hot spot functions’ explicit executing time before and after mesh reordering

針對(duì)不同網(wǎng)格規(guī)模下的顯式計(jì)算，記錄了24進(jìn)程并行計(jì)算的運(yùn)行時(shí)間，結(jié)果如表 5所示。

Table 5 Executing time of explicit computing before and after mesh reordering

對(duì)未排序運(yùn)行時(shí)間作無因次化處理，重排序、未排序顯式計(jì)算的運(yùn)行時(shí)間如圖 13所示。網(wǎng)格重排序縮短了顯式計(jì)算在所有網(wǎng)格規(guī)模下的運(yùn)行時(shí)間,隨著網(wǎng)格規(guī)模擴(kuò)大，重排序能夠節(jié)省15%～20%的運(yùn)行時(shí)間。

重排序改變了數(shù)據(jù)在內(nèi)存中的分布，增強(qiáng)了局部性，提高了Cache命中率，減小了訪存延遲。特別是在網(wǎng)格規(guī)模較大的時(shí)候，內(nèi)存讀取次數(shù)增多，能夠節(jié)省更多的運(yùn)行時(shí)間。

Figure 13 Comparison of explicit parallel computing’s executing time before and after mesh reordering on CPU

表6統(tǒng)計(jì)了CPU隱式計(jì)算在網(wǎng)格3上的熱點(diǎn)函數(shù)網(wǎng)格重排序前、后的運(yùn)行時(shí)間。通過圖 14可知，網(wǎng)格重排序可以使部分熱點(diǎn)函數(shù)的CPU運(yùn)行時(shí)間縮短。FI0梯度計(jì)算的最明顯，運(yùn)行時(shí)間縮短了大約17%；LU-SGS部分相關(guān)的FI3和FI4運(yùn)行時(shí)間縮短了10%以上。

Table 6 Comparison of hot spot functions’ implicit executing time before and after mesh reordering on CPU

Figure 14 Comparison of hot spot functions’ implicit executing time before and after mesh reordering on CPU

本文也對(duì)網(wǎng)格重排序?qū)﹄[式計(jì)算的影響進(jìn)行了評(píng)估。針對(duì)隱式計(jì)算，開展了不同規(guī)模的并行計(jì)算。記錄時(shí)間如表 7所示。

對(duì)未排序運(yùn)行時(shí)間作無因次化處理，重排序、未排序隱式計(jì)算的運(yùn)行時(shí)間如圖 15所示。對(duì)于隱式計(jì)算，網(wǎng)格重排序仍然能夠使縮短運(yùn)行時(shí)間。

Table 7 Comparison of implicit parallel computing’s executing time before and after mesh reordering on CPU

對(duì)于規(guī)模較小的網(wǎng)格1(100萬(wàn)體網(wǎng)格)，運(yùn)行時(shí)間縮短接近10%，比顯式計(jì)算的優(yōu)化效果更明顯。這可能是由于體網(wǎng)格排序?qū)τ诰仃嚽蠼獾膬?yōu)化即使在網(wǎng)格規(guī)模較小的時(shí)候仍然起作用。隨著網(wǎng)格規(guī)模增大，重排序效果也更加明顯，在網(wǎng)格4規(guī)模下，重排序減少了將近20%的運(yùn)行時(shí)間。

Figure 15 Comparison of implicit parallel computing’s executing time before and after mesh reordering on CPU

4.3.2 網(wǎng)格重排序?qū)τ贕PU計(jì)算性能的影響

GPU的SIMT(Single Instruction Multiple Threads)能夠使多個(gè)線程同時(shí)訪問GPU全局內(nèi)存，提高了數(shù)據(jù)吞吐能力。GPU的計(jì)算部件SM訪問片上內(nèi)存的帶寬通常高于CPU的訪存帶寬，因此，有必要研究網(wǎng)格重排序?qū)PU計(jì)算的影響。本節(jié)只討論顯式計(jì)算在GPU上的性能。

顯式計(jì)算在單GPU上400萬(wàn)網(wǎng)格規(guī)模下8個(gè)熱點(diǎn)kernel的運(yùn)行時(shí)間如表 8所示。無因次化(未排序kernel計(jì)算時(shí)間為標(biāo)準(zhǔn))的kernel運(yùn)行時(shí)間比較如圖 16所示。由圖16可知，網(wǎng)格重排序使得主要熱點(diǎn)kernel的運(yùn)行時(shí)間縮短了35%～60%。

Table 8 Executing time of kernels for reordering and unordering

Figure 16 Comparison of hot spot kernels’ executing time before and after mesh reordering on GPU

網(wǎng)格重排序優(yōu)化了數(shù)據(jù)在GPU全局內(nèi)存中的分布，減少了數(shù)據(jù)間接訪問帶來的訪存延遲，能夠縮短大部分熱點(diǎn)kernel的運(yùn)行時(shí)間。

本文還記錄了重排序前、后CFD程序在單GPU上的整體運(yùn)行時(shí)間，如表 9所示。無因次化的運(yùn)行時(shí)間如圖 17所示。網(wǎng)格重排序使得顯式求解CFD程序在GPU上的整體性能得到了提升，運(yùn)行時(shí)間平均縮短了約40%。

Table 9 Comparison of explicit parallel computing’s executing time before and after grid reordering on GPU

Figure 17 Comparison of explicit parallel computing’s executing time before and after mesh reordering on GPU

大部分熱點(diǎn)kernel的運(yùn)行時(shí)間縮短最終使得CFD程序的整體性能得到了提升。因此，網(wǎng)格編號(hào)重排序是一種全局性的優(yōu)化策略。

網(wǎng)格重排序在GPU上的加速效果相比在CPU上的更加明顯。這是GPU的訪存模式造成的：GPU以32個(gè)線程組成1個(gè)warp，作為基本單位連續(xù)訪問全局內(nèi)存，直到訪問完所需數(shù)據(jù)；GPU的L1 Cache只能保證空間局部性，不能保證時(shí)間局部性。因此，數(shù)據(jù)局部性的提高減少了GPU內(nèi)存的訪問頻率，降低了訪存延遲。對(duì)于CPU，數(shù)據(jù)局部性的優(yōu)化只是提升了L1 Cache的空間局部性，對(duì)于時(shí)間局部性改善不大，訪存延遲減少也有限。

5 結(jié)束語(yǔ)

本文采用網(wǎng)格重排序的方法優(yōu)化了數(shù)據(jù)在內(nèi)存中的分布。針對(duì)空氣動(dòng)力學(xué)中典型的三維機(jī)翼網(wǎng)格，分別在CPU和GPU上展開了多種規(guī)模的計(jì)算,發(fā)現(xiàn)網(wǎng)格重排序能提升數(shù)據(jù)局部性，減少訪存延遲，提高計(jì)算效率。具體來說：

(1)網(wǎng)格重排序?qū)τ?jì)算精度的影響較小：排序前、后隱式計(jì)算的差異在0.5%以內(nèi)，可以認(rèn)為對(duì)顯式計(jì)算的結(jié)果幾乎沒有影響。

(2)網(wǎng)格重排序能縮短計(jì)算時(shí)間。在CPU上，當(dāng)體網(wǎng)格規(guī)模大于100萬(wàn)時(shí)，網(wǎng)格重排序能夠縮短15%～20%的計(jì)算時(shí)間；在GPU上，主要計(jì)算熱點(diǎn)在網(wǎng)格重排序后計(jì)算時(shí)間縮短了35%～60%，程序的整體運(yùn)行時(shí)間也相應(yīng)縮短了40%左右。

(3)網(wǎng)格重排序?qū)PU顯式計(jì)算的優(yōu)化效果更明顯。

經(jīng)過重排序后，部分函數(shù)的性能下降，這表明排序算法和具體函數(shù)的數(shù)據(jù)使用密切相關(guān)。未來會(huì)結(jié)合具體函數(shù)的數(shù)據(jù)使用特點(diǎn)開發(fā)、評(píng)估更多的排序算法。

致謝：感謝風(fēng)雷軟件團(tuán)隊(duì)和國(guó)家超級(jí)計(jì)算廣州中心。