面向國(guó)產(chǎn)加速器的CFD核心算法并行優(yōu)化

曹義魁，陸忠華，張鑒，劉夏真，袁武，梁姍

1.中國(guó)科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心，北京 100190

2.中國(guó)科學(xué)院大學(xué)，北京 100049

引言

計(jì)算流體動(dòng)力學(xué)（Computational Fluid Dynamics,CFD）是一門新興的交叉型學(xué)科，它結(jié)合數(shù)學(xué)中的離散方法，利用計(jì)算機(jī)強(qiáng)大的算力對(duì)流體力學(xué)中復(fù)雜的微分方程進(jìn)行求解，在科研和工程領(lǐng)域都發(fā)揮著巨大的作用。目前的CFD 商業(yè)軟件、開(kāi)源CFD軟件等，大多數(shù)是基于CPU 開(kāi)發(fā)的。但隨著計(jì)算規(guī)模的不斷擴(kuò)大，CPU 的計(jì)算、訪存及通信性能已經(jīng)不能滿足人們的需要[1]，所以尋求新的方法實(shí)現(xiàn)對(duì)大規(guī)模數(shù)據(jù)進(jìn)行并行處理已經(jīng)成為國(guó)產(chǎn)CFD 軟件發(fā)展的必要條件。

CPU+加速器異構(gòu)架構(gòu)的出現(xiàn)，為人們解決復(fù)雜的大規(guī)模計(jì)算提供了新的方法，它充分集合了CPU 和加速器的各自優(yōu)勢(shì)特點(diǎn)，將具有并行性的計(jì)算密集部分放到加速器上進(jìn)行加速，極大的提高了程序的運(yùn)行速度，這也成為了目前高性能計(jì)算的主流方法[2]。為了能夠充分利用加速器的性能，需要針對(duì)程序中的算法特點(diǎn)設(shè)計(jì)具有高并行度、高訪存帶寬的并行方法。國(guó)內(nèi)外已有較多將CFD 應(yīng)用移植到CPU+加速器異構(gòu)平臺(tái)并進(jìn)行優(yōu)化的相關(guān)研究。

Corrigan 等[3]基于CPU/GPU 異構(gòu)平臺(tái)對(duì)流體力學(xué)中的格子玻爾茲曼方法實(shí)現(xiàn)了并行加速，并采用多種方法實(shí)現(xiàn)了簡(jiǎn)單高效的數(shù)據(jù)訪存模式，使全局訪存帶寬利用率達(dá)到了86%，相應(yīng)程序也得到33.6 倍的加速效果。Jespersen[4]設(shè)計(jì)并實(shí)現(xiàn)了基于對(duì)稱逐次超松弛方法的CFD 軟件的GPU 并行方案，并對(duì)其中的Jacobi 迭代部分進(jìn)行了并行加速，單GPU 卡相對(duì)于單CPU 核心實(shí)現(xiàn)了大約3 倍的加速效果。李大力等[5]針對(duì)高階精度加權(quán)緊致非線性格式的CFD 軟件，實(shí)現(xiàn)并優(yōu)化了Jacobi 解法器的 GPU 并行計(jì)算，與單CPU 核心相比運(yùn)行時(shí)間加速了9.8 倍。董廷星等人[6]利用GPU 加速計(jì)算流體力學(xué)中經(jīng)典的N-S(Navier-Stokes)方程和歐拉(Euler)方程的求解，并采用3 個(gè)測(cè)試算例進(jìn)行實(shí)驗(yàn)，單GPU 卡相對(duì)于單CPU 核心最高能得到33.2 倍的加速效果。鄧亮等人[7]基于ADI 解法器的有限體積CFD 應(yīng)用，設(shè)計(jì)了兩種GPU并行方案，并討論了若干性能優(yōu)化方法，使整個(gè)CFD 應(yīng)用得到了17.3 倍的加速效果。V.Emelyanov[8]在GPU 上開(kāi)發(fā)了基于有限體積法的非定?？蓧嚎s的歐拉方程和N-S 方程求解器，并利用精簡(jiǎn)計(jì)算、訪存優(yōu)化等方法對(duì)GPU 程序進(jìn)行優(yōu)化，最終得到20 ～50 倍的加速效果。此外，他們還在GPU上對(duì)不同聲速下的翼型外流場(chǎng)進(jìn)行了數(shù)值模擬[9]。Lai 等[10]基于CPU/GPU 異構(gòu)平臺(tái)，利用雷諾數(shù)為118 的雙橢球繞流問(wèn)題驗(yàn)證了可壓縮NS 方程求解器的計(jì)算能力，并從數(shù)據(jù)傳輸方面對(duì)程序進(jìn)行了優(yōu)化。隨后，他們又設(shè)計(jì)了高效的多GPU 并行算法對(duì)高超聲速流場(chǎng)進(jìn)行了研究，當(dāng)使用4 張GPU 加速卡時(shí)，可以獲得147 倍的加速[11]。黨冠麟等[12]基于CPU/GPU 異構(gòu)系統(tǒng)，利用自主開(kāi)發(fā)的高精度有限差分CFD 求解器，對(duì)鈍錐邊界層轉(zhuǎn)捩問(wèn)題進(jìn)行了數(shù)值模擬，并對(duì)核函數(shù)進(jìn)行了細(xì)致的優(yōu)化，單GPU 卡相對(duì)單CPU 核心獲得了60 倍的加速。

以上工作都是在CPU+GPU 的異構(gòu)系統(tǒng)上完成。為了打破高性能行業(yè)的技術(shù)封鎖，近年來(lái)國(guó)家高度重視國(guó)產(chǎn)高性能計(jì)算機(jī)的發(fā)展，并取得了突出成就，多次斬獲超級(jí)計(jì)算TOP500 榜單首名，并且我國(guó)在E 級(jí)超級(jí)計(jì)算機(jī)系統(tǒng)的研發(fā)上也處于世界領(lǐng)先地位。使用國(guó)產(chǎn)加速器替代GPU 等加速器已經(jīng)成為一種必然趨勢(shì)[13]。本文利用CPU+國(guó)產(chǎn)加速器的異構(gòu)系統(tǒng)對(duì)國(guó)產(chǎn)自主CFD 軟件進(jìn)行加速，國(guó)產(chǎn)加速器采用類GPU體系架構(gòu)，相關(guān)資料尚未被授權(quán)公開(kāi)，此處對(duì)其架構(gòu)不展開(kāi)介紹。國(guó)產(chǎn)加速器可以完全兼容運(yùn)行在GPU上的CUDA 程序，其優(yōu)化思路也大致和GPU 相同。

1 CFD 軟件介紹

1.1 CFD 軟件計(jì)算流程

本文中用于移植的CFD 軟件的初始版本是CCFD V3.0，此軟件全部用Fortran 語(yǔ)言編寫。CCFD系列軟件是國(guó)家“863 計(jì)劃”“十一五”和“十二五”連續(xù)支持的面向大型飛機(jī)設(shè)計(jì)的大規(guī)模流場(chǎng)數(shù)值模擬軟件[14]，也是“863 計(jì)劃”中我國(guó)集中科研資源重點(diǎn)發(fā)展的高性能計(jì)算軟件之一[15]。課題組在CCFD V1.0 和V2.0 的基礎(chǔ)上，瞄準(zhǔn)國(guó)產(chǎn)高性能計(jì)算平臺(tái)發(fā)展V3.0 版本，并在神威-太湖之光等平臺(tái)上開(kāi)展了相關(guān)研究。CCFD V3.0 軟件的計(jì)算流程主要包括三個(gè)階段，如下圖1所示：

圖1 CCFD V3.0 軟件定常計(jì)算流程圖Fig.1 CCFD V3.0 software steady state calculation flow chart

（1）預(yù)處理階段。主要進(jìn)行網(wǎng)格劃分，設(shè)置計(jì)算參數(shù)，進(jìn)行度量系數(shù)的計(jì)算。

（2）初始化階段。進(jìn)行流場(chǎng)信息的初始化，處理初始邊界條件。

（3）迭代計(jì)算階段。包括邊界更新、時(shí)間推進(jìn)、對(duì)流項(xiàng)的計(jì)算、方程求解等模塊。

1.2 CFD 核心算法

根據(jù)程序特點(diǎn)，本文對(duì)CCFD V3.0 的定常求解器steady 部分做了移植。經(jīng)過(guò)測(cè)試發(fā)現(xiàn)，flux 函數(shù)部分和lu_adi 函數(shù)部分約占整個(gè)程序運(yùn)行時(shí)間的90%，所以本文主要介紹有關(guān)對(duì)流項(xiàng)計(jì)算的flux 函數(shù)部分以及使用LU_ADI 算法求解的lu_adi 函數(shù)部分，這兩部分也是整個(gè)CFD 求解器的核心部分。

對(duì)流項(xiàng)的空間離散方法在CFD 計(jì)算中至關(guān)重要，它不僅會(huì)影響計(jì)算的穩(wěn)定性，對(duì)計(jì)算結(jié)果在精度上也會(huì)有很大影響。CCFD V3.0 軟件集成了多種國(guó)際上評(píng)價(jià)較高的上風(fēng)格式（Roe’s FDS[16]、Van Leer’s FVS[17]、AUSM+[18]）和Jameson 中心格式[19]。其中，本文研究是Roe’s FDS 格式，它是基于黎曼解的通量差分分裂格式，是目前使用最廣泛、評(píng)價(jià)最高的迎風(fēng)格式之一。

時(shí)間推進(jìn)方法上，工程CFD 軟件普遍采用隱式算法，典型的隱式算法有對(duì)角化方法、交替隱式追趕方法（ADI）[20]和LU-SGS[21]方法等。其中，交替方向隱式追趕方法（ADI）對(duì)通過(guò)系數(shù)矩陣在三個(gè)方向上作近似因子（AF）分解，并基于LU 分裂將方程求解分為上下三角的兩個(gè)子過(guò)程，穩(wěn)定性好且計(jì)算效率高，在不要求時(shí)間精度的定常計(jì)算中應(yīng)用廣泛。

flux 函數(shù)的對(duì)流項(xiàng)計(jì)算部分是典型的模版計(jì)算，lu_adi 函數(shù)中的LU_ADI 算法是典型的數(shù)據(jù)依賴方法，這兩部分是CFD 軟件基礎(chǔ)算法的主要部分。由于推導(dǎo)過(guò)程會(huì)占據(jù)大量篇幅，所以此處不再贅述，本文只在下文相應(yīng)章節(jié)處給出具體的相應(yīng)代碼，詳細(xì)的公式推導(dǎo)及具體算法介紹見(jiàn)參考文獻(xiàn)18。

2 核心算法并行方案

國(guó)產(chǎn)加速器架構(gòu)類似于GPU，都是通過(guò)將數(shù)據(jù)分散到成千上萬(wàn)個(gè)小的計(jì)算核心實(shí)現(xiàn)并行加速。本文的移植方案先將軟件移植到CPU+GPU 異構(gòu)系統(tǒng)得到基礎(chǔ)版本，之后再移植到國(guó)產(chǎn)加速器上進(jìn)行優(yōu)化。

GPU 使用CUDA 編程模型[22]，線程按照線程格（Grid）、線程塊（Block）、線程（Thread）的多層次模型進(jìn)行組織。其中，Block 以三維的形式組織在同一個(gè)Grid 內(nèi)，Thread 以三維的形式組織在同一個(gè)Block 內(nèi)，并分別通過(guò)內(nèi)置變量blockIdx(x,y,z)和threadIdx(x,y,z)進(jìn)行標(biāo)識(shí)，從而可以對(duì)每個(gè)線程進(jìn)行索引。

下面將分別介紹flux 函數(shù)和lu_adi 函數(shù)兩部分在GPU 上的并行方案。

2.1 flux 函數(shù)的GPU 并行方案

對(duì)流通量項(xiàng)計(jì)算函數(shù)（flux 函數(shù)）是一個(gè)典型的模版計(jì)算函數(shù)，每個(gè)網(wǎng)格數(shù)據(jù)點(diǎn)的計(jì)算都是獨(dú)立進(jìn)行的，僅需要將相應(yīng)網(wǎng)格點(diǎn)的數(shù)據(jù)映射到GPU 上的單個(gè)線程，可以實(shí)現(xiàn)高效的三維并行，具體的映射關(guān)系如圖2所示。根據(jù)flux 函數(shù)的計(jì)算流程，將每個(gè)方向的計(jì)算都拆分成了4 個(gè)核函數(shù)，分別命名為flux_kernel1-flux_kernel4，拆分的作用是為了實(shí)現(xiàn)對(duì)內(nèi)點(diǎn)和邊界點(diǎn)的分別計(jì)算以及數(shù)據(jù)計(jì)算的最大并行。

圖2 右端項(xiàng)計(jì)算的線程映射Fig.2 Thread mapping of right end item calculation

2.2 lu_adi 函數(shù)的GPU 并行方案

在CPU 版本中，lu_adi 函數(shù)分別在X、Y、Z三個(gè)方向利用LU_ADI 方法求解三維定常可壓縮N-S 方程。該算法最大的特點(diǎn)是在每個(gè)計(jì)算方向上都存在強(qiáng)數(shù)據(jù)依賴性，網(wǎng)格中任一內(nèi)點(diǎn)的更新都需要自身點(diǎn)和同一條網(wǎng)格線上鄰點(diǎn)參與計(jì)算，且需要往返各更新一次，具體的計(jì)算方式如下圖3所示。

圖3 ADI 迭代計(jì)算Fig.3 ADI iterative computation

針對(duì)上述ADI 算法特點(diǎn)，對(duì)lu_adi 函數(shù)中沒(méi)有數(shù)據(jù)依賴的計(jì)算部分仍采用三維并行，網(wǎng)格數(shù)據(jù)與線程的映射和圖3 相同。而對(duì)lu_adi 函數(shù)中具有強(qiáng)數(shù)據(jù)依賴性的計(jì)算部分，只能在沒(méi)有數(shù)據(jù)依賴的另外兩個(gè)方向上實(shí)現(xiàn)兩維并行，有數(shù)據(jù)依賴的方向使用循環(huán)串行執(zhí)行，每個(gè)線程負(fù)責(zé)控制一條網(wǎng)格線上相關(guān)點(diǎn)的計(jì)算，具體的網(wǎng)格數(shù)據(jù)與線程的映射關(guān)系如圖4所示。本文根據(jù)lu_adi 函數(shù)特點(diǎn)，將每個(gè)方向的計(jì)算都拆分成了5 個(gè)核函數(shù)，分別命名為ADI_kernel1-ADI_kernel5，拆分的目的是為了實(shí)現(xiàn)數(shù)據(jù)計(jì)算的最大并行化。

圖4 LU_ADI 迭代計(jì)算的線程映射Fig.4 Thread mapping of LU_ ADI iterative calculation

3 基于國(guó)產(chǎn)加速器的性能優(yōu)化策略

利用上述并行方案實(shí)現(xiàn)了GPU 的加速版本后，再通過(guò)國(guó)產(chǎn)加速器自帶的轉(zhuǎn)碼工具將程序移植到國(guó)產(chǎn)加速器，得到了移植后的基礎(chǔ)版本。使用128 立方的三維網(wǎng)格進(jìn)行測(cè)試，單加速卡相對(duì)于單CPU 核心，移植前后flux 函數(shù)部分的運(yùn)行時(shí)間從2231.21s 降為16.9s，加速了132 倍，lu_adi 函數(shù)部分的運(yùn)行實(shí)踐從1641.93s 降到83.1s，加速了19.7 倍。但此基礎(chǔ)版本仍有很大的優(yōu)化空間，后續(xù)本文會(huì)基于國(guó)產(chǎn)加速器架構(gòu)特點(diǎn)，充分利用國(guó)產(chǎn)加速器上寄存器、共享內(nèi)存等內(nèi)存資源，采用核函數(shù)合并與分解、訪存合并、調(diào)整塊大小等方法對(duì)移植后的程序進(jìn)行優(yōu)化。

3.1 核函數(shù)的分解與合并

在對(duì)核函數(shù)進(jìn)行優(yōu)化時(shí)，應(yīng)充分利用寄存器資源進(jìn)行計(jì)算。但是寄存器資源非常稀少，所以當(dāng)核函數(shù)過(guò)大，中間變量過(guò)多時(shí)，可以采用拆分核函數(shù)的方法，減少同一個(gè)核函數(shù)內(nèi)寄存器的使用數(shù)量，從而增大線程塊并發(fā)的數(shù)量。當(dāng)核函數(shù)規(guī)模較小，中間變量較少，可以將相關(guān)核函數(shù)進(jìn)行合并，在保證線程塊具有較高并行性的情況下，充分利用寄存器資源。本課題中對(duì)lu_adi 函數(shù)進(jìn)行了相關(guān)優(yōu)化，lu_adi 函數(shù)共由五個(gè)核函數(shù)組成，前兩個(gè)核函數(shù)ADI_kernel1 和ADI_kernel2 的計(jì)算中有大量中間數(shù)組的重復(fù)使用。針對(duì)此特點(diǎn)，可以將這些使用到的中間數(shù)組數(shù)據(jù)用寄存器存儲(chǔ)起來(lái)，之后的相關(guān)計(jì)算也都在寄存器上進(jìn)行。此外，ADI_kernel2 核函數(shù)中的一些計(jì)算可以合并到ADI_kernel1 核函數(shù)中，并且不會(huì)對(duì)程序結(jié)果產(chǎn)生影響。通過(guò)測(cè)試發(fā)現(xiàn)，將這兩種優(yōu)化方法結(jié)合起來(lái)使用，這兩個(gè)核函數(shù)的總時(shí)間從16.54 秒降到10.95 秒，速度提升了33.8%。

3.2 全局訪存的合并與對(duì)齊

通常情況下，計(jì)算數(shù)據(jù)都是存儲(chǔ)在加速器的全局內(nèi)存（Global Memory）上，并且內(nèi)存帶寬往往是影響程序性能的主要因素。因此，設(shè)計(jì)較好的訪存模式實(shí)現(xiàn)較高的訪存帶寬對(duì)程序的優(yōu)化至關(guān)重要。在GPU 上核函數(shù)的內(nèi)存訪問(wèn)是在物理設(shè)備內(nèi)存和片上內(nèi)存間以128 字節(jié)或32 字節(jié)內(nèi)存事物來(lái)實(shí)現(xiàn)的，在國(guó)產(chǎn)加速器上的訪存機(jī)制與GPU 類似。為了最大化帶寬利用率，我們可以對(duì)程序的訪存模式進(jìn)行設(shè)計(jì)，將設(shè)備內(nèi)存事務(wù)的首地址設(shè)置為128 字節(jié)的偶數(shù)倍，并讓同一線程束內(nèi)連續(xù)的線程獲取連續(xù)的內(nèi)存數(shù)據(jù)，從而實(shí)現(xiàn)訪存的合并與對(duì)齊，提高核函數(shù)的運(yùn)行速度。在我們的程序中，以ADI_kernel3_i 核函數(shù)為例,原程序中對(duì)殘差數(shù)組三個(gè)維度的計(jì)算都從第三個(gè)點(diǎn)開(kāi)始，其對(duì)數(shù)據(jù)的訪問(wèn)是非對(duì)齊的。為了實(shí)現(xiàn)合并對(duì)齊訪問(wèn)，我們對(duì)x 方向的計(jì)算仍從第三個(gè)點(diǎn)開(kāi)始，另外兩個(gè)方向的計(jì)算都從第一個(gè)點(diǎn)開(kāi)始，這樣既不會(huì)影響程序的正確性，又可以實(shí)現(xiàn)對(duì)全局內(nèi)存的合并對(duì)齊訪問(wèn)。經(jīng)過(guò)這樣的優(yōu)化后，ADI_kernel3_i 核函數(shù)的運(yùn)行時(shí)間從1.86s 降為1.58s，性能提升了15%。

3.3 共享內(nèi)存的使用

共享內(nèi)存是國(guó)產(chǎn)加速器上的一塊具有固定大小的可編程緩存，與全局內(nèi)存相比，它具有較低的訪存延遲和較高的內(nèi)存帶寬。存儲(chǔ)在共享內(nèi)存上的數(shù)據(jù)可以被同一個(gè)線程塊內(nèi)的線程共享。國(guó)產(chǎn)加速器的一個(gè)計(jì)算卡上有64kB 的共享內(nèi)存，當(dāng)數(shù)據(jù)是雙精度浮點(diǎn)類型時(shí)，一個(gè)核函數(shù)使用的共享內(nèi)存數(shù)組大小不能超過(guò)8192。根據(jù)程序特點(diǎn)，使用共享內(nèi)存存儲(chǔ)全局內(nèi)存數(shù)據(jù)，可以減少對(duì)Global Memory 不必要的頻繁訪問(wèn)，提升程序的整體性能。特別是一些具有數(shù)據(jù)依賴的計(jì)算，使用共享內(nèi)存進(jìn)行優(yōu)化是優(yōu)先選擇。但是使用共享內(nèi)存時(shí)還要注意避免bank 沖突，通常解決bank 沖突的方法是對(duì)數(shù)據(jù)進(jìn)行內(nèi)存填充，改變數(shù)據(jù)到共享內(nèi)存存儲(chǔ)體的映射。此外，還可以利用共享內(nèi)存對(duì)內(nèi)存訪問(wèn)模式進(jìn)行設(shè)計(jì)，以便實(shí)現(xiàn)訪存的合并。

3.3.1 對(duì)lu_adi 函數(shù)使用共享內(nèi)存優(yōu)化

在本課題中，lu_adi 函數(shù)的數(shù)據(jù)依賴部分主要集中在第四個(gè)核函數(shù)ADI_kernel4，并且對(duì)i 方向數(shù)據(jù)的訪問(wèn)是非合并的，具體的計(jì)算形式如下所示。

Specific Code:do k = 2 ,kdim do j = 2 ,jdim do i = 3 ,idim res(i,j,k)=(res(i,j,k)- flmp(i,j,k)* res(i-1,j,k))/diag(i,j,k)end do do i = idim-1 ,2 res(i,j,k)= res(i,j,k)- flmn(i,j,k)* res(i+1,j,k)end do end do end do

針對(duì)i 方向訪存不連續(xù)的問(wèn)題，本文使用共享內(nèi)存設(shè)計(jì)了一種實(shí)現(xiàn)合并的內(nèi)存訪問(wèn)模式。由于一個(gè)核函數(shù)的共享內(nèi)存大小不能超過(guò)64kB，所以沿i正向計(jì)算時(shí)使用三個(gè)大小為4*32*17 的共享內(nèi)存數(shù)組，分別存放res、flmp和diag數(shù)組，并且控制連續(xù)的16 個(gè)線程訪問(wèn)連續(xù)的16 個(gè)全局內(nèi)存數(shù)據(jù)，將其存儲(chǔ)到共享內(nèi)存數(shù)組中，這樣便實(shí)現(xiàn)了合并訪存。i方向的每一條網(wǎng)格線分成n 次來(lái)計(jì)算，n 的大小由網(wǎng)格大小確定。每次循環(huán)沿i 方向取16 個(gè)數(shù)據(jù)，并且在進(jìn)行下次計(jì)算之前將需要用到的數(shù)據(jù)提前存放到共享內(nèi)存數(shù)組中，這也是將共享內(nèi)存數(shù)組的最后一個(gè)維度聲明為17 的原因。n 次循環(huán)之后，一個(gè)線程實(shí)際負(fù)責(zé)了一條網(wǎng)格線上數(shù)據(jù)點(diǎn)的計(jì)算，反向計(jì)算原理相同。此外，n 次循環(huán)中，每次取數(shù)據(jù)之后要進(jìn)行同步，計(jì)算之后也要進(jìn)行同步，否則可能會(huì)出現(xiàn)未知錯(cuò)誤，同步的方法是使用syncthreads()函數(shù)，此函數(shù)的作用是同步同一線程塊內(nèi)的線程操作。設(shè)計(jì)的具體訪存模式如圖5所示。

圖5 ADI_kernel4 核函數(shù)i 方向訪存模式Fig.5 ADI_ Kernel4 kernel function i-direction memory access mode

而對(duì)于j、k 方向數(shù)據(jù)的訪存本就是合并的，使用共享內(nèi)存僅僅只能減少對(duì)全局內(nèi)存的訪存次數(shù)，而使用寄存器可以達(dá)到同樣的效果，并且比共享內(nèi)存更加方便，所以對(duì)j、k 方向使用寄存器進(jìn)行計(jì)算，每次沿相應(yīng)方向連續(xù)更新計(jì)算一條網(wǎng)格線上的所有點(diǎn)。使用此方法進(jìn)行優(yōu)化后，ADI_kernel4 部分實(shí)現(xiàn)了5.4 倍的加速。

3.3.2 對(duì)flux 函數(shù)使用共享內(nèi)存優(yōu)化

flux 函數(shù)中的第一個(gè)核函數(shù)flux_kernel1 的計(jì)算需要用到相同數(shù)組的相鄰兩點(diǎn)，具體計(jì)算形式如下所示。

Specific Code:do k = 2 ,kdim do j = 2 ,jdim do i = 3 ,idim du(i,j,k)= u(i,j,k)- u(i-1,j,k)end do end do end do

并且數(shù)據(jù)的訪問(wèn)已經(jīng)實(shí)現(xiàn)了合并與對(duì)齊，利用共享內(nèi)存可以減少對(duì)全局內(nèi)存的訪問(wèn)次數(shù)，從而減少程序的運(yùn)行時(shí)間。具體方案和對(duì)ADI_kernel4 的優(yōu)化思路大致相同，不同之處在于使用的共享內(nèi)存數(shù)組大小為4*8*33，并且沒(méi)有循環(huán)過(guò)程，也沒(méi)有數(shù)據(jù)依賴。使用共享內(nèi)存進(jìn)行優(yōu)化后，flux_kernel1 部分相比原來(lái)的版本實(shí)現(xiàn)了28%的加速。

3.4 調(diào)整線程塊的大小

線程塊大小會(huì)影響每個(gè)線程所能分配的內(nèi)核資源，從而影響整個(gè)程序的并發(fā)性。調(diào)整線程塊大小是比較基礎(chǔ)的優(yōu)化方法，當(dāng)已經(jīng)使用了其他優(yōu)化方法或者找不到更好的優(yōu)化方法時(shí)，對(duì)線程塊大小進(jìn)行調(diào)整通?？梢匀〉揭欢ǖ男Ч?。本課題對(duì)lu_adi 函數(shù)中有關(guān)i 方向計(jì)算的兩個(gè)比較簡(jiǎn)單的核函數(shù)ADI_kernel1_i 和ADI_kernel5_i 采用此方法進(jìn)行優(yōu)化，并取得了一定的加速效果。具體的線程塊配置如表1所示。

表1 ADI_kernel1_i 和ADI_kernel5_i 在不同線程塊大小下的運(yùn)行時(shí)間Table 1 The running time of ADI_kernel1_i and ADI_kernel5_i under different thread block sizes

從圖中可以看到，相同的線程塊大小配置對(duì)于不同的核函數(shù)影響是不同的。ADI_kernel1_i 核函數(shù)中計(jì)算操作較多，用到了較多的寄存器，使用32*8*2 的線程網(wǎng)格大小既可以增大每個(gè)線程分配的寄存器的資源，又可以得到最大的并行性，所以其加速效果也最好。而ADI_kernel3_i 核函數(shù)中計(jì)算較少，主要是讀寫內(nèi)存操作，使用32*8*4 的線程網(wǎng)格可以保證其最大并行性。所以，對(duì)每個(gè)核函數(shù)應(yīng)根據(jù)自身特點(diǎn)，合理設(shè)計(jì)線程塊大小，設(shè)計(jì)相應(yīng)核函數(shù)的最佳線程塊配置，而不應(yīng)該采用統(tǒng)一的線程網(wǎng)格大小。

4 實(shí)驗(yàn)結(jié)果與分析

4.1 實(shí)驗(yàn)平臺(tái)

在CPU+國(guó)產(chǎn)加速器異構(gòu)平臺(tái)進(jìn)行測(cè)試，系統(tǒng)環(huán)境配置：國(guó)產(chǎn)x86 處理器(32 核心)，內(nèi)存容量為128GB DDR4，國(guó)產(chǎn)加速器具有16GB HBM2 顯存，帶寬1TB/s。本文采用Intel 編譯器編譯Fortran 程序，使用國(guó)產(chǎn)加速器自帶的編譯器編譯移植后的程序，并且使用“-O2”優(yōu)化，最終將編譯后的文件用ifort 命令進(jìn)行鏈接生成可執(zhí)行文件。在CPU 上所進(jìn)行的運(yùn)算只使用單核，國(guó)產(chǎn)加速器上使用單張加速卡進(jìn)行加速。

4.2 正確性驗(yàn)證

本文的算例是求解三維方腔流問(wèn)題，使用了128*128*128 和256*256*128 兩種規(guī)模的網(wǎng)格模型進(jìn)行測(cè)試，其中網(wǎng)格的最外面兩層為虛網(wǎng)格點(diǎn)，其作用是用來(lái)更新內(nèi)點(diǎn)，網(wǎng)格數(shù)據(jù)全部為雙精度，設(shè)置定常迭代步數(shù)為2000。通過(guò)打印輸出原CPU 程序和移植優(yōu)化后的最終全部流場(chǎng)信息，驗(yàn)證移植與優(yōu)化的正確性。兩種測(cè)試結(jié)果表明，移植前后全部流場(chǎng)信息的絕對(duì)誤差的最大值都控制在小數(shù)點(diǎn)后15 位，完全滿足工程上的精度要求。

4.3 實(shí)驗(yàn)結(jié)果

本文采用兩套不同網(wǎng)格規(guī)模進(jìn)行測(cè)試，圖6 和圖7 給出了兩種網(wǎng)格移植前后CCFD V3.0 程序各部分的運(yùn)行時(shí)間對(duì)比。通過(guò)對(duì)測(cè)試結(jié)果的初步分析，發(fā)現(xiàn)128 立方大小的網(wǎng)格和256*256*128 大小的網(wǎng)格的優(yōu)化效果基本一致。限于篇幅，本文只給出128 立方大小網(wǎng)格詳細(xì)的優(yōu)化結(jié)果。表2 是對(duì)ADI_kernel1 和ADI_kernel2 采用核函數(shù)的合并與分解方法的具體優(yōu)化效果，表3、表4 分別是對(duì)ADI_kernel4、flux_kernel1 核函數(shù)使用共享內(nèi)存方法優(yōu)化的三個(gè)方向優(yōu)化效果。圖8 和圖9 分別給出了移植的初始版本和綜合使用上述幾種優(yōu)化方法后，flux 函數(shù)和lu_adi 函數(shù)各個(gè)方向計(jì)算的最終時(shí)間對(duì)比。

圖8 flux 函數(shù)優(yōu)化前后時(shí)間對(duì)比Fig.8 Comparison of time before and after optimization of flux function

圖9 lu_adi 函數(shù)優(yōu)化前后時(shí)間對(duì)比Fig.9 Comparison of time before and after optimization of lu_adi function

表2 ADI_kernel1 和ADI_kernel2 優(yōu)化效果Table 2 The optimization effect of ADI_ Kernel1 and ADI_kernel2

表3 ADI_kernel4 核函數(shù)三個(gè)方向優(yōu)化效果Table 3 The optimization effect of ADI_kernel4 kernel function in three directions

表4 flux_kernel1 核函數(shù)三個(gè)方向優(yōu)化效果Table 4 The optimization effect of flux_kernel1 kernel function in three directions

圖6 CCFD V3.0 三個(gè)版本程序時(shí)間對(duì)比，網(wǎng)格規(guī)模：128＊128＊128Fig.6 Time comparison of three versions of CCFD V3.0,grid size:128＊128＊128

圖7 CCFD V3.0 三個(gè)版本程序時(shí)間對(duì)比，網(wǎng)格規(guī)模：256＊256＊128Fig.7 Time comparison of three versions of CCFD V3.0,grid size:256＊256＊128

從上述圖表可以看出將移植到國(guó)產(chǎn)異構(gòu)系統(tǒng)可以對(duì)程序的性能有明顯的提升，并且采用上述優(yōu)化方法進(jìn)行優(yōu)化后，程序性能會(huì)進(jìn)一步得到提升。其中,對(duì)lu_adi 函數(shù)i 方向的相關(guān)計(jì)算在優(yōu)化后的性能提升最明顯，這是因?yàn)閕 方向的數(shù)據(jù)在初步移植后的訪存是非合并的，優(yōu)化后實(shí)現(xiàn)了訪存的合并，性能得到了很好的提升。

5 實(shí)驗(yàn)結(jié)果與分析

本文將CCFD V3.0 軟件的定常求解器部分移植到國(guó)產(chǎn)異構(gòu)平臺(tái)，并針對(duì)CFD 軟件的主要核心算法，對(duì)flux 函數(shù)部分和lu_adi 函數(shù)部分實(shí)現(xiàn)了在國(guó)產(chǎn)異構(gòu)平臺(tái)上的并行優(yōu)化。根據(jù)異構(gòu)平臺(tái)特點(diǎn)，使用了寄存器、共享內(nèi)存、對(duì)核函數(shù)進(jìn)行合并與分解、實(shí)現(xiàn)訪存的合并與對(duì)齊、調(diào)整線程塊大小等優(yōu)化方法，對(duì)移植后的程序進(jìn)行優(yōu)化，最終取得了很好的加速效果，flux 函數(shù)部分實(shí)現(xiàn)了166 倍的加速，lu_adi 函數(shù)部分實(shí)現(xiàn)了59 倍的加速，整個(gè)CCFD V3.0 程序?qū)崿F(xiàn)了90 倍的加速（國(guó)產(chǎn)加速器上的單個(gè)加速卡相對(duì)于單個(gè)CPU 核心）。

下一步將整個(gè)CCFD V3.0 程序移植到CPU+國(guó)產(chǎn)加速器的異構(gòu)平臺(tái)，并將上述優(yōu)化方法應(yīng)用到軟件的其它模塊中，從而對(duì)整個(gè)CCFD V3.0 應(yīng)用程序?qū)崿F(xiàn)更好的加速效果。

利益沖突聲明

所有作者聲明不存在利益沖突關(guān)系。