基于OPS的計(jì)算流體力學(xué)軟件多平臺(tái)自動(dòng)并行*

王 巍，車永剛，徐傳福，王正華

(國(guó)防科技大學(xué)計(jì)算機(jī)學(xué)院量子信息研究所兼高性能計(jì)算國(guó)家重點(diǎn)實(shí)驗(yàn)室,湖南 長(zhǎng)沙 410073)

1 引言

當(dāng)前，除了通用處理器之外，通用圖形處理單元GPGPU(General-Purpose Graphic Processing Unit)、Intel集成眾核MIC(Many Integrated Core)、通用數(shù)字信號(hào)處理器GPDSP(General- Purpose Digital Signal Processing)等加速部件被廣泛應(yīng)用于高性能計(jì)算機(jī)系統(tǒng)中，高性能計(jì)算機(jī)變得越來越異構(gòu)和多樣化[1 - 3]。因此，高性能計(jì)算應(yīng)用軟件的開發(fā)將面向多樣化的體系結(jié)構(gòu)。而高性能計(jì)算應(yīng)用軟件的生存期通常比硬件長(zhǎng)得多[4]，針對(duì)不同體系結(jié)構(gòu)的并行編程與性能優(yōu)化導(dǎo)致極大的開銷。因此，使所研制的高性能計(jì)算應(yīng)用軟件具有性能可移植性(即代碼無需經(jīng)過大量修改或使用不同的編程語言重寫，就可移植到不同體系結(jié)構(gòu)平臺(tái)上并獲得可接受的性能)，是發(fā)展高性能數(shù)值模擬的迫切需求[5]。

現(xiàn)今主流的并行軟件開發(fā)方法是由開發(fā)人員根據(jù)軟件的算法特點(diǎn)和流程等設(shè)計(jì)程序并行結(jié)構(gòu)，在源代碼中手工添加MPI、OpenMP、CUDA、OpenCL等并行編程語句的代碼，進(jìn)行手工并行編程。手工并行編程要求開發(fā)人員熟練掌握前述等并行編程規(guī)范和語言，并深入了解底層硬件體系結(jié)構(gòu)才能編寫高效率的代碼，且開發(fā)工作量大、周期長(zhǎng)，很難獲得性能可移植性。

基于高性能庫(kù)(如BLAS(Basic Linear Algebra Subroutine)、FFTW(Faster Fourier Transform in the West)等)編程的方法在一些情況下可以方便地獲得性能可移植性，但是可用的高性能可移植庫(kù)非常有限，并且對(duì)實(shí)際應(yīng)用軟件來說庫(kù)調(diào)用通常只能覆蓋其部分計(jì)算任務(wù)?；诳梢浦簿幊炭蚣艿姆椒ㄊ橇硪环N途徑。例如美國(guó)E級(jí)計(jì)算計(jì)劃中Kokkos項(xiàng)目[6]提供了共享存儲(chǔ)和分布存儲(chǔ)編程模型以及相應(yīng)的C++庫(kù)，使用Kokkos編寫的單個(gè)并行代碼經(jīng)過少量面向體系結(jié)構(gòu)的修改，就可在多種體系結(jié)構(gòu)平臺(tái)上編譯執(zhí)行[5]，但可移植編程框架通常只支持C++編程語言，并且其針對(duì)不同體系結(jié)構(gòu)的優(yōu)化能力還比較弱。

基于領(lǐng)域特定語言DSL(Domain Specific Language，或稱領(lǐng)域?qū)ｉT語言)的并行化方法提供了一種較為可行的途徑，得到研究界的關(guān)注。領(lǐng)域特定語言針對(duì)某一領(lǐng)域的特定問題，提供一套簡(jiǎn)潔抽象的高層語義結(jié)構(gòu)，開發(fā)人員使用其來編程描述所要計(jì)算的問題，然后由一個(gè)低層軟件來自動(dòng)將這些描述轉(zhuǎn)換成不同平臺(tái)上的實(shí)現(xiàn)，使軟件具有良好的性能可移植性。目前研究人員提出了許多面向高性能計(jì)算的DSL框架。如牛津大學(xué)開發(fā)的結(jié)構(gòu)網(wǎng)格下并行求解庫(kù)OPS(Oxford Parallel library for Structured mesh solvers)[7,8]和非結(jié)構(gòu)網(wǎng)格問題求解的DSL框架OP2[9]、斯坦福大學(xué)開發(fā)的Liszt[10]、弗吉尼亞理工大學(xué)等開發(fā)的基于網(wǎng)格的語言與自動(dòng)并行化框架GLAF(Grid-based Language and Auto-parallelization Framework)[11]、日本東京工業(yè)大學(xué)等開發(fā)的針對(duì)層次式多體算法的隱式并行編程框架Tapas[12]、路易斯安那州立大學(xué)等開發(fā)的基于高性能計(jì)算機(jī)體系結(jié)構(gòu)的偏微分方程求解DSL框架Chemora[13]等。

本文旨在探索基于領(lǐng)域特定語言O(shè)PS實(shí)現(xiàn)計(jì)算流體動(dòng)力學(xué)CFD(Computational Fluid Dynamics)軟件面向多平臺(tái)自動(dòng)并行編程技術(shù)，評(píng)估OPS針對(duì)CFD應(yīng)用軟件在各種并行體系結(jié)構(gòu)上的并行能力，探索高性能計(jì)算應(yīng)用開發(fā)的新途徑。本文選擇了一個(gè)結(jié)構(gòu)化網(wǎng)格下基于有限差分的高精度CFD應(yīng)用軟件HNSC (High order Navier-Stokes simulator for Compressible flow)，使用OPS API對(duì)其進(jìn)行了轉(zhuǎn)換，基于OPS前端和后端面向純CPU平臺(tái)和GPU平臺(tái)生成了并行代碼，并與其手工實(shí)現(xiàn)的并行版本的性能進(jìn)行了測(cè)試比較。測(cè)試結(jié)果說明，對(duì)所評(píng)價(jià)的應(yīng)用軟件，基于OPS的方法獲得的性能與手工版本性能可比，并且具有跨CPU與GPU可移植的優(yōu)勢(shì)。

2 OPS簡(jiǎn)介

OPS是英國(guó)牛津大學(xué)等開發(fā)的一套嵌入在C/C++及Fortran語言中，可為多區(qū)結(jié)構(gòu)網(wǎng)格應(yīng)用程序生成底層并行代碼的DSL。它包含源到源的轉(zhuǎn)換程序，可將基于OPS框架開發(fā)的源代碼轉(zhuǎn)換成在各種并行編程模式下實(shí)現(xiàn)的并行代碼。用戶提供包含所要進(jìn)行計(jì)算的全部代碼和對(duì)OPS函數(shù)的調(diào)用，而OPS負(fù)責(zé)控制對(duì)給定體系結(jié)構(gòu)上的各個(gè)計(jì)算循環(huán)的并行化。在使用上，OPS類似于傳統(tǒng)的編程框架，包括相應(yīng)的接口以及支持這些接口的庫(kù)。

2.1 結(jié)構(gòu)網(wǎng)格下計(jì)算過程的抽象

OPS將多塊結(jié)構(gòu)網(wǎng)格下的科學(xué)計(jì)算抽象成4個(gè)主要元素，即：網(wǎng)格塊、定義在網(wǎng)格塊上的數(shù)據(jù)集、數(shù)據(jù)交換區(qū)和格式計(jì)算操作。OPS將有序排列的多塊結(jié)構(gòu)網(wǎng)格看作是一系列無序排列的網(wǎng)格塊的組合，通過相應(yīng)的接口定義網(wǎng)格塊之間的數(shù)據(jù)交換行為，使其能正確地描述多塊網(wǎng)格下的計(jì)算，網(wǎng)格塊給出所計(jì)算問題的維度，計(jì)算規(guī)模則是通過數(shù)據(jù)集的定義給出的，數(shù)據(jù)集是計(jì)算的操作對(duì)象，數(shù)據(jù)交換區(qū)則是網(wǎng)格塊之間的數(shù)據(jù)交換行為的具體定義，格式計(jì)算操作是對(duì)數(shù)據(jù)集中網(wǎng)格點(diǎn)數(shù)據(jù)訪問方式和操作的描述。

OPS在結(jié)構(gòu)網(wǎng)格上描述格式計(jì)算操作時(shí)，先定義一個(gè)計(jì)算格式，即計(jì)算一個(gè)網(wǎng)格點(diǎn)數(shù)據(jù)時(shí)，要訪問相鄰的網(wǎng)格點(diǎn)的個(gè)數(shù)以及這些相鄰網(wǎng)格點(diǎn)相對(duì)于要計(jì)算的網(wǎng)格點(diǎn)的索引偏移。這種描述契合有限差分等科學(xué)計(jì)算本身的邏輯特點(diǎn)。

2.2 自動(dòng)代碼生成

OPS提供一個(gè)源對(duì)源的前端代碼轉(zhuǎn)換器，可以將基于OPS接口編制的代碼自動(dòng)轉(zhuǎn)換成多種并行模型和規(guī)范下的代碼，包括針對(duì)CPU平臺(tái)的MPI、OpenMP，針對(duì)GPU平臺(tái)的CUDA、OpenACC和OpenCL。代碼轉(zhuǎn)換器在生成的代碼中自動(dòng)添加上述編程規(guī)范的語句指導(dǎo)程序異構(gòu)和并行，不需要開發(fā)者再手工修改代碼。進(jìn)而基于OPS為不同平臺(tái)提供后端運(yùn)行時(shí)庫(kù)，采用本地編譯器(如GNU編譯器、Intel編譯器等)，將這些代碼編譯生成面向各種硬件平臺(tái)的可執(zhí)行程序。

2.3 OPS的優(yōu)點(diǎn)

OPS對(duì)多塊結(jié)構(gòu)網(wǎng)格下科學(xué)計(jì)算程序的抽象符合科學(xué)計(jì)算本身的邏輯，容易被特定領(lǐng)域開發(fā)者所接受和使用。在當(dāng)前高性能計(jì)算平臺(tái)向多樣化發(fā)展的形勢(shì)下，開發(fā)者采用OPS開發(fā)維護(hù)一套源代碼就可以得到面向多種平臺(tái)的可執(zhí)行代碼，極大地提升了軟件開發(fā)效率。同時(shí)，OPS自動(dòng)實(shí)現(xiàn)多種模式的并行代碼，其前端代碼轉(zhuǎn)換和后端庫(kù)實(shí)現(xiàn)中都包含了多種性能優(yōu)化技術(shù)，使特定領(lǐng)域開發(fā)人員無需掌握并行編程的細(xì)節(jié)就可以開發(fā)出高效的并行代碼，從而使其可以將精力集中在自身領(lǐng)域的問題解決上。

3 HNSC軟件的OPS實(shí)現(xiàn)

HNSC是一個(gè)空氣動(dòng)力學(xué)數(shù)值模擬的In-House軟件，該軟件基于由斯坦福大學(xué)Sv?rd、NASA蘭利研究中心Carpenter等[14]提出的穩(wěn)定有限差分格式，求解全三維可壓縮Navier Stokes方程組，滿足能量估計(jì)和分部求和規(guī)則[16]，在航空航天CFD計(jì)算領(lǐng)域具有良好的應(yīng)用背景。前期工作中，我們已經(jīng)使用手工方式實(shí)現(xiàn)了該軟件面向CPU的MPI/OpenMP混合并行，使用CUDA進(jìn)行了面向GPU平臺(tái)的編程。而本文研究采用OPS 對(duì)HNSC進(jìn)行多平臺(tái)并行編程，以一個(gè)自由剪切流算例來說明HNSC軟件實(shí)現(xiàn)以及如何基于OPS框架來編程。

3.1 HNSC軟件本身的實(shí)現(xiàn)

HNSC軟件基于結(jié)構(gòu)化網(wǎng)格構(gòu)建一種具有高階精度的中心差分格式來計(jì)算控制方程中的各階偏導(dǎo)數(shù)項(xiàng)，其時(shí)間推進(jìn)采用的是4階段龍格庫(kù)塔顯式算法。自由剪切流算例的計(jì)算區(qū)域可以采用一個(gè)二維的矩形區(qū)域來表示，并且軟件已經(jīng)基于區(qū)域劃分手工實(shí)現(xiàn)了MPI并行，手工實(shí)現(xiàn)的MPI并行構(gòu)架如圖1所示，每一個(gè)進(jìn)程計(jì)算一個(gè)子區(qū)域網(wǎng)格上的數(shù)據(jù)，相鄰的網(wǎng)格塊或者進(jìn)程(如P4和P7進(jìn)程)之間通過MPI消息傳遞接口進(jìn)行數(shù)據(jù)交換，為了滿足自由剪切流的邊界條件，邊界上的子區(qū)域基于周期性MPI拓?fù)鋪磉M(jìn)行數(shù)據(jù)交換，如圖1中P2和P8進(jìn)程所計(jì)算的子區(qū)域。

Figure 1 MPI communication involved in HNSC software圖1 HNSC軟件本身實(shí)現(xiàn)的MPI通信行為

基于圖1所示的區(qū)域劃分，在每一個(gè)網(wǎng)格區(qū)域上，HNSC軟件分別對(duì)流場(chǎng)數(shù)據(jù)進(jìn)行求解，軟件的執(zhí)行過程如圖2所示。首先初始化MPI并行環(huán)境，Grid函數(shù)讀入網(wǎng)格，在進(jìn)行計(jì)算前通過Init函數(shù)對(duì)流場(chǎng)的基本未知量進(jìn)行初始化，初始化之后的基本未知量通過Solvec函數(shù)組裝成方便控制方程求解的守恒變量。隨后進(jìn)入主循環(huán)迭代，主循環(huán)中Compute_dt函數(shù)基于場(chǎng)變量的當(dāng)前值計(jì)算時(shí)間步長(zhǎng)，RK4函數(shù)是HNSC軟件計(jì)算的核心，包含4階龍格庫(kù)塔算法的全部計(jì)算，并在一個(gè)時(shí)間步內(nèi)更新場(chǎng)變量的數(shù)值。4階龍格庫(kù)塔計(jì)算之后，進(jìn)行收斂檢驗(yàn)，進(jìn)而計(jì)算下一個(gè)時(shí)間步的步長(zhǎng)，并將迭代時(shí)間值向前推進(jìn)，如果達(dá)到最大迭代步數(shù)，計(jì)算結(jié)束。

Figure 2 Procedure of HNSC software圖2 HNSC軟件的執(zhí)行流程

圖2右側(cè)是RK4函數(shù)內(nèi)部的具體執(zhí)行內(nèi)容，decomp是將基本未知量從守恒變量中還原出來的函數(shù)，F(xiàn)lux_vec計(jì)算控制方程中的通量項(xiàng)，F(xiàn)lux_der函數(shù)計(jì)算通量的偏導(dǎo)數(shù)。4階龍格庫(kù)塔算法的4階計(jì)算具有相似的流程，為了簡(jiǎn)潔圖2中僅展示了第1階計(jì)算的過程，值得注意的是在4階龍格庫(kù)塔算法的每一階計(jì)算中，都包含2次邊界條件施加函數(shù)BC的執(zhí)行，這對(duì)后面將軟件和算例基于OPS實(shí)現(xiàn)是十分重要的。

3.2 HNSC軟件基于OPS的實(shí)現(xiàn)

結(jié)合二維自由剪切流算例下HNSC軟件的行為與流程，基于OPS對(duì)其進(jìn)行實(shí)現(xiàn)包括以下步驟和關(guān)鍵技術(shù)：

(1)網(wǎng)格塊的聲明；

(2)定義數(shù)據(jù)集、數(shù)據(jù)交換區(qū)及通信；

(3)定義計(jì)算格式；

(4)差分計(jì)算循環(huán)的重構(gòu)；

(5)并行代碼編譯生成。

3.2.1 網(wǎng)格塊的聲明

在添加相應(yīng)的頭文件之后，可以使用OPS提供的編程接口聲明一個(gè)名為grid_block的網(wǎng)格塊：

intgrid_dim;

grid_dim= 3;

ops_blockgrid_block=ops_decl_block(grid_dim,“grid_block”);

其類型是ops_block，函數(shù)接口ops_decl_block的參數(shù)grid_dim是網(wǎng)格塊的維度，可以是1,2或3(一般CFD的求解最多針對(duì)三維空間來求解)，針對(duì)HNSC軟件，grid_dim設(shè)為3。字符型參數(shù)“grid_block”是該block的標(biāo)簽名，用于代碼轉(zhuǎn)換、編譯以及數(shù)據(jù)I/O時(shí)的顯示和輸出。

3.2.2 定義數(shù)據(jù)集、數(shù)據(jù)交換區(qū)及通信

HNSC作為CFD軟件所求解的基本未知量包括密度、速度的XYZ3個(gè)坐標(biāo)軸分量、壓力、溫度和內(nèi)能，分別使用符號(hào)rho、u、v、w、p、Temp、e來表示。并且軟件還將涉及CFD求解的過程變量，如通量項(xiàng)等。這些變量在軟件代碼中作為雙精度浮點(diǎn)型數(shù)組來計(jì)算?；贠PS來實(shí)現(xiàn)時(shí)，也要定義對(duì)應(yīng)的數(shù)據(jù)集，例如密度rho的定義：

double *rho;

intbuffer;

intsize[]={nx,ny,nz};

intbase[]={0,0,0};

inthalo_pos[]={buffer,buffer,buffer};

inthalo_neg[]={-buffer,-buffer,-buffer};

ops_datops_rho=ops_decl_dat(grid_block,1,size,base,halo_pos,halo_neg,rho,“double”,“rho”);

這樣，我們?cè)诿麨間rid_block的網(wǎng)格塊上定義了一個(gè)名為ops_rho的數(shù)據(jù)集，元素?cái)?shù)據(jù)類型和指針rho一致。size是這個(gè)數(shù)據(jù)集的大小或各個(gè)坐標(biāo)方向網(wǎng)格點(diǎn)的數(shù)目，注意grid_block本身的維度是3維，所以size應(yīng)定義3個(gè)維度的大小，分別為nx、ny和nz。base是該數(shù)據(jù)集第1個(gè)元素的索引。halo_pos是沿著循環(huán)索引增大的方向數(shù)據(jù)集和數(shù)據(jù)交換區(qū)的大小，halo_neg是沿著循環(huán)索引負(fù)向的數(shù)據(jù)交換區(qū)的大小，其3個(gè)維度的值都是buffer，buffer的取值是由HNSC軟件的差分計(jì)算的格式確定的，其值不小于坐標(biāo)軸方向上的差分點(diǎn)個(gè)數(shù)，一般取等于差分點(diǎn)個(gè)數(shù)即可，對(duì)不參與格式計(jì)算的數(shù)據(jù)集，buffer的值可以是0。值得注意的是，該數(shù)據(jù)集定義接口中的第2個(gè)參數(shù)1，這意味著一個(gè)網(wǎng)格點(diǎn)對(duì)應(yīng)該數(shù)據(jù)集的1個(gè)值。因?yàn)镠NSC軟件的有限差分法計(jì)算中，密度值和網(wǎng)格點(diǎn)是一一對(duì)應(yīng)的。

求解HNSC軟件時(shí)，采用守恒形式的控制方程，解向量、3個(gè)方向的通量和源項(xiàng)分別用符號(hào)Uvec、Fvec、Gvec、Hvec和Jvec來表示，對(duì)于每一個(gè)網(wǎng)格點(diǎn)，這些變量是一個(gè)長(zhǎng)度nvar=5的列向量，因此對(duì)解向量應(yīng)按以下方式定義數(shù)據(jù)集：

double *Uvec,*Fvec,*Gvec,*Hvec,*Jvec;

intnvar=5;

ops_datops_Uvec=ops_decl_dat(grid_block,nvar,size,base,halo_pos,halo_neg,Uvec,“double”,“Uvec”);

因?yàn)镸PI并行是OPS自動(dòng)實(shí)現(xiàn)的，不再需要關(guān)注軟件實(shí)現(xiàn)中的任務(wù)劃分，所以在定義數(shù)據(jù)集的時(shí)候，數(shù)據(jù)集的規(guī)模nx、ny、nz是圖1中整個(gè)計(jì)算區(qū)域所有網(wǎng)格點(diǎn)的規(guī)模。OPS實(shí)現(xiàn)MPI并行時(shí)會(huì)根據(jù)數(shù)據(jù)交換區(qū)的信息halo_pos和halo_neg來定義MPI通信緩沖區(qū)。在網(wǎng)格塊的內(nèi)部，OPS將根據(jù)差分計(jì)算的循環(huán)自行定義MPI通信行為的觸發(fā)。因此，在定義好數(shù)據(jù)集和數(shù)據(jù)交換區(qū)大小之后，圖1中HNSC軟件原本負(fù)責(zé)MPI通信的fill_buffer函數(shù)可以直接去掉。但是，由于自由剪切流算例在區(qū)域邊界存在通信行為，圖1中BC函數(shù)不能去掉，該函數(shù)涉及到OPS中網(wǎng)格塊外部邊界的數(shù)據(jù)交換區(qū)通信，需要顯式地定義和觸發(fā)。OPS提供了函數(shù)ops_halo_transfer來觸發(fā)網(wǎng)格塊外部邊界的數(shù)據(jù)交換行為，該函數(shù)支持將發(fā)送和接收數(shù)據(jù)的網(wǎng)格塊設(shè)置成同一個(gè)網(wǎng)格塊的數(shù)據(jù)集，這樣便可以實(shí)現(xiàn)HNSC軟件中BC函數(shù)所定義的數(shù)據(jù)交換行為：

ops_halogrouprho_halo_xalong;

ops_halopos_rho=ops_decl_halo(ops_rho,ops_rho,xhalo.iter,xhalo.pbase_from,xhalo.pbase_to,xhalo.axes,xhalo.axes);

ops_haloneg_rho=ops_decl_halo(ops_rho,ops_rho,xhalo.iter,xhalo.nbase_from,xhalo.nbase_to,xhalo.axes,xhalo.axes);

ops_halogrp[] ={pos_rho,neg_rho};

rho_halo_xalong=ops_decl_halo_group(2,grp);

ops_halo_transfer(rho_halo_xalong);

上面的代碼定義了對(duì)grid_block上的數(shù)據(jù)集ops_rho自身在網(wǎng)格塊邊界附近X坐標(biāo)軸方向的數(shù)據(jù)交換，在pos_rho=ops_decl_halo接口中定義了一個(gè)從數(shù)據(jù)集ops_rho(第1個(gè)參數(shù))發(fā)送一個(gè)在X軸正向側(cè)邊界附近，大小為xhalo.iter的數(shù)據(jù)塊到數(shù)據(jù)集ops_rho(第2個(gè)參數(shù))，xhalo.pbase_from為數(shù)據(jù)塊發(fā)送的起點(diǎn)，xhalo.pbase_to為數(shù)據(jù)接收的起點(diǎn)，xhalo.axes為發(fā)生數(shù)據(jù)交換行為的數(shù)據(jù)集之間的坐標(biāo)軸對(duì)應(yīng)關(guān)系。同樣，neg_rho=ops_decl_halo定義了數(shù)據(jù)集ops_rho在X軸負(fù)向側(cè)邊界附近的數(shù)據(jù)塊交換。OPS要求將同步執(zhí)行的邊界數(shù)據(jù)交換放在一個(gè)ops_halo_group類型中，再使用ops_halo_transfer函數(shù)來觸發(fā)。Y方向、Z方向以及其它的數(shù)據(jù)集都按照需求定義如上數(shù)據(jù)交換行為后，將這些內(nèi)容寫入BC函數(shù)中，便實(shí)現(xiàn)了自由剪切流施加邊界條件的BC函數(shù)基于OPS的改寫，記為ops_BC。

3.2.3 定義計(jì)算格式

OPS通過定義計(jì)算格式來描述類似于有限差分計(jì)算，這樣計(jì)算某一個(gè)網(wǎng)格點(diǎn)數(shù)據(jù)時(shí)，使用周圍特定相鄰網(wǎng)格點(diǎn)的數(shù)據(jù)的計(jì)算行為。具體如下所示：

ops_stencilS3D_000,S3D_5pt,S3D_25pt;

ints3d_000[] = {0,0,0};

S3D_000 =ops_decl_stencil(3,1,s3d_000,"0,0,0");

ints3d_25pt[] = {0,0,0,1,0,0,-1,0,0,2,0,0,-2,0,0,3,0,0,-3,0,0,4,0,0,-4,0,0,5,0,0,-5,0,0,6,0,0,-6,0,0,0,1,0,0,-1,0,0,2,0,0,-2,0,0,3,0,0,-3,0,0,4,0,0,-4,0,0,5,0,0,-5,0,0,6,0,0,-6,0};

S3D_25pt=ops_decl_stencil(3,25,s3d_25pt,"6,0,0:-6,0,0;0,6,0:0,-6,0");

計(jì)算格式為ops_stencil類型，S3D_000 =ops_decl_stencil(3,1,s3d_000,“0,0,0”)定義了一個(gè)三維網(wǎng)格塊中，包含1個(gè)網(wǎng)格點(diǎn)，索引為s3d_000，標(biāo)簽名為“0,0,0”的計(jì)算格式，{0,0,0}是該計(jì)算格式中心網(wǎng)格點(diǎn)的坐標(biāo)。計(jì)算格式中網(wǎng)格點(diǎn)的坐標(biāo)是以中心網(wǎng)格點(diǎn)為基準(zhǔn)的相對(duì)坐標(biāo)或相對(duì)索引。之所以定義一個(gè)只有1個(gè)網(wǎng)格點(diǎn)的計(jì)算格式S3D_000是為了計(jì)算時(shí)進(jìn)行賦值操作。S3D_25pt是HNSC軟件計(jì)算二維問題時(shí)使用的計(jì)算格式，整型數(shù)組s3d_25pt定義了該計(jì)算格式中各個(gè)網(wǎng)格點(diǎn)相對(duì)中心點(diǎn)的索引偏移值。

3.2.4 差分計(jì)算循環(huán)的重構(gòu)

OPS對(duì)軟件中遍歷網(wǎng)格點(diǎn)的差分計(jì)算循環(huán)使用在某個(gè)(或多個(gè))數(shù)據(jù)集上、基于相應(yīng)計(jì)算格式的訪問來描述，用戶定義一個(gè)kernel函數(shù)來描述對(duì)計(jì)算格式的具體操作，供OPS循環(huán)接口調(diào)用，從而實(shí)現(xiàn)原代碼中的計(jì)算循環(huán)。HNSC中計(jì)算速度分量u對(duì)坐標(biāo)x的偏導(dǎo)數(shù)的原代碼片段如下所示：

for(k=kstart;k<=kend;k++){

for(j=jstart;j<=jend;j++){

for(i=istart;i<=iend;i++){

dudx[i][j][k]=a1_24*(u[i+1][j][k]-u[i-1][j][k])+a2_24*(u[i+2][j][k]-u[i-2][j][k])+a3_24*(u[i+3][j][k]-u[i-3][j][k])+a4_24*(u[i+4][j][k]-u[i-4][j][k])+a5_24*(u[i+5][j][k]-u[i-5][j][k])+a6_24*(u[i+6][j][k]-u[i-6][j][k] )/dx;

}

}

}

使用OPS循環(huán)接口ops_para_loop改寫該計(jì)算循環(huán)如下：

intloop_range[]={0,nx,0,ny,0,nz};

ops_par_loop(kernel_uderivs,"kernels_uderivs",grid_block,3,loop_range,ops_arg_dat(ops_dudx,1,S3D_000,"double",OPS_WRITE),ops_arg_dat(ops_u,1,S3D_25pt,"double",OPS_READ));

OPS提供接口ops_par_loop來改寫需要并行化的循環(huán)，它指定該循環(huán)是定義在grid_block網(wǎng)格塊上執(zhí)行范圍為loop_range的一個(gè)3維(層)循環(huán)。ops_arg_dat傳遞該循環(huán)需要操作的數(shù)據(jù)集，并規(guī)定了針對(duì)該數(shù)據(jù)集的計(jì)算格式為S3D_000或S3D_25pt，訪問方式為：OPS_WRITE(寫入)、OPS_READ(讀取)或OPS_RW(讀寫)。原循環(huán)體內(nèi)部的具體計(jì)算，即對(duì)計(jì)算格式的具體操作在函數(shù)kernel_uderivs中，由用戶依照OPS的編程方式來定義：

voidkernel_uderivs(double *dudx_in,const double *u_in)

{

dudx_in[OPS_ACC0(0,0,0)] = (a1_24*(u_in[OPS_ACC1(1,0,0)]-u_in[OPS_ACC1(-1,0,0)])+a2_24*(u_in[OPS_ACC1(2,0,0)]-u_in[OPS_ACC1(-2,0,0)])+a3_24*(u_in[OPS_ACC1(3,0,0)]-u_in[OPS_ACC1(-3,0,0)])+a4_24*(u_in[OPS_ACC1(4,0,0)]-u_in[OPS_ACC1(-4,0,0)])+a5_24*(u_in[OPS_ACC1(5,0,0)]-u_in[OPS_ACC1(-5,0,0)])+a6_24*(u_in[OPS_ACC1(6,0,0)]-u_in[OPS_ACC1(-6,0,0)]) )/dx;

}

OPS_ACC是OPS提供的一個(gè)宏定義，在使用的時(shí)候要根據(jù)函數(shù)的形參位置指定序號(hào)(如OPS_ACC0)，并給出計(jì)算格式中網(wǎng)格點(diǎn)的相對(duì)索引(如OPS_ACC0(1,0,0))。需要特別注意的是，ops_par_loop循環(huán)的范圍loop_range為{0,nx,0,ny,0,nz}，遍歷的是整個(gè)網(wǎng)格塊，而{istart,iend,jstart,jend,kstart,kend}是HNSC手工進(jìn)行區(qū)域分解的子區(qū)域的計(jì)算循環(huán)邊界。

對(duì)圖2中HNSC軟件中的各個(gè)函數(shù)，都按照這樣的模式來改寫，并對(duì)其進(jìn)行了封裝。本文僅對(duì)軟件主要的實(shí)現(xiàn)技術(shù)進(jìn)行闡述，OPS還提供了相應(yīng)的接口和數(shù)據(jù)類型支持HNSC軟件中的歸約計(jì)算操作等，具體使用規(guī)范可查閱OPS手冊(cè)。在HNSC軟件的所有函數(shù)基于OPS實(shí)現(xiàn)之后，OPS提供ops_init和ops_exit來執(zhí)行相應(yīng)的初始化和結(jié)束工作，MPI并行時(shí)進(jìn)程通信的拓?fù)浣Y(jié)構(gòu)的建立也是OPS自動(dòng)完成的。我們將基于OPS編寫的HNSC軟件稱為HNSC_OPS。

3.2.5 并行代碼編譯生成

OPS約定所有關(guān)于格式計(jì)算具體操作的kernel函數(shù)都以頭文件的形式來編寫，并將項(xiàng)目中的文件分為HEADERS、OPS_FILE和OTHER_FILES。HEADERS即包含用戶kernel函數(shù)的頭文件，OPS_FILES是包含了ops_par_loop接口描述的需要并行化的循環(huán)文件，變量聲明、OPS相關(guān)定義等所在的文件則放在OTHER_FILES中。用戶必須按照這樣的分類來使用OPS框架提供的模板化的Makefile，否則不能正確地進(jìn)行代碼轉(zhuǎn)換與編譯。模板化的Makefile將調(diào)用代碼轉(zhuǎn)換器進(jìn)行代碼轉(zhuǎn)換，代碼轉(zhuǎn)換器以O(shè)PS_FILES作為輸入，并根據(jù)這些文件中對(duì)kernel函數(shù)的調(diào)用，尋找相應(yīng)的頭文件，并對(duì)頭文件中的kernel函數(shù)進(jìn)行相應(yīng)的轉(zhuǎn)換，以生成面向多體系結(jié)構(gòu)平臺(tái)的代碼，再根據(jù)不同的平臺(tái)調(diào)用相應(yīng)的本地編譯器生成面向多平臺(tái)的可執(zhí)行文件。

基于OPS的前后端，本文生成得到了HNSC軟件的純MPI并行版本HNSC_OPS_MPI、純OpenMP并行版本HNSC_OPS_OpenMP、MPI/OpenMP 2級(jí)混合并行版本HNSC_OPS_MPI_OpenMP、CUDA并行版本HNSC_OPS_CUDA。

4 性能測(cè)試與分析

本文在一個(gè)配備有2塊Intel Xeon E5-2660 v3(Haswell)CPU和1塊NVIDIA Tesla K80 GPU的服務(wù)器上進(jìn)行了性能測(cè)試，該服務(wù)器的硬件和軟件配置如表1所示。這里不僅測(cè)試了基于OPS生成的HNSC并行代碼的性能，也測(cè)試了手工實(shí)現(xiàn)的HNSC軟件的多種并行版本的性能，以進(jìn)行比較分析，包括手工實(shí)現(xiàn)的純MPI并行版本HNSC_Manual_MPI、手工實(shí)現(xiàn)的純OpenMP并行版本HNSC_Manual_OpenMP、手工實(shí)現(xiàn)的MPI/OpenMP 2級(jí)混合并行版本HNSC_Manual_MPI_OpenMP，以及手工移植的面向GPU并行的版本HNSC_Manual_CUDA。測(cè)試的網(wǎng)格規(guī)模分別為2M和4M網(wǎng)格點(diǎn)。

Table 1 Configuration of the server表1 服務(wù)器的配置

圖3給出了純MPI并行情況下，OPS生成的MPI并行代碼與手工MPI并行代碼的性能，2個(gè)版本的代碼都采用的是表1中的Intel編譯器，并且采用相同的編譯選項(xiàng)：-no-prec-div-restrict-fno-alias-fma-fp-model fast=2。從圖3可以看到，在1,2和4進(jìn)程情況下，OPS自動(dòng)生成的并行代碼相比手工并行代碼的性能稍低。在8,16和20進(jìn)程情況下，OPS生成代碼的性能優(yōu)于手工并行代碼的。在2M網(wǎng)格規(guī)模下，OPS生成代碼20進(jìn)程相對(duì)單進(jìn)程的加速比為10.5，手工并行代碼20進(jìn)程相對(duì)于單進(jìn)程的并行加速比為4.2。在4M網(wǎng)格規(guī)模下，OPS生成代碼20進(jìn)程相對(duì)于單進(jìn)程的并行加速比為11.3，手工并行代碼20進(jìn)程相對(duì)于單進(jìn)程的并行加速比為4.4。2種網(wǎng)格規(guī)模下，OPS生成的代碼用滿20核時(shí)獲得的并行加速比都比手工并行代碼的高。

Figure 3 MPI parallel performance comparison of OPS code and manual code 圖3 純MPI并行情況下OPS生成并行代碼與手工并行代碼的性能對(duì)比

圖4給出了純OpenMP并行情況下，OPS生成的代碼與手工OpenMP并行代碼的性能，同樣2個(gè)版本的代碼都采用的是表1中的Intel編譯器，編譯選項(xiàng)增加了OpenMP開關(guān)：-no-prec-div- restrict-fno-alias-fma-fp-model fast=2-qopenmp，測(cè)試時(shí)的線程數(shù)都是通過環(huán)境變量來指定的。從圖4中可以看到，在2種網(wǎng)格規(guī)模下，隨著線程數(shù)的增加，OPS生成的并行代碼逐漸表現(xiàn)出比手工OpenMP并行代碼更好的性能。2M網(wǎng)格規(guī)模下OPS生成代碼用滿20核相對(duì)于單核的并行加速比為4.4，高于手工并行代碼的3.8。4M網(wǎng)格規(guī)模下OPS生成代碼用滿20核相對(duì)于單核的并行加速比為5.7，高于手工并行代碼的4.6。

Figure 4 OpenMP parallel performance comparison of OPS code and manual code圖4 純OpenMP并行情況下OPS生成并行代碼與手工并行代碼的性能對(duì)比

圖5給出了MPI/OpenMP 2級(jí)混合并行情況下，OPS生成的并行代碼和手工并行代碼的性能對(duì)比，也是使用表1中的Intel編譯器進(jìn)行編譯，編譯選項(xiàng)相同：-no-prec-div-restrict-fno-alias-fma-fp-model fast=2-qopenmp。在2M網(wǎng)格規(guī)模下，OPS生成并行代碼與手工并行代碼的性能在某些進(jìn)程/線程組合時(shí)性能不如手工并行代碼的，但整體幾乎相當(dāng)。而在4M網(wǎng)格規(guī)模下，OPS生成的并行代碼在各種進(jìn)程/線程組合下的性能都優(yōu)于手工并行代碼的。

Figure 5 MPI/OpenMP parallel performance comparison of OPS code and manual code圖5 MPI/OpenMP混合并行情況下OPS生成并行代碼和手工并行代碼的性能對(duì)比

圖6給出了OPS自動(dòng)生成的CUDA代碼與手工編寫的CUDA代碼在NVIDIA K80 GPU上的性能對(duì)比，編譯時(shí)，對(duì)于GPU端執(zhí)行的核函數(shù)OPS版本和手工版本都采用表1中CUDA版本中的nvcc編譯器，測(cè)試的GPU為Kepler架構(gòu)，所以關(guān)于架構(gòu)信息的編譯選項(xiàng)為：-gencode arch=compute_37,code=sm_37，另外添加了性能選項(xiàng)：--use_fast_math，主編譯器及選項(xiàng)和MPI版本代碼相同：-no-prec-div-restrict-fno-alias-fma-fp-model fast=2。作為參考，圖6中也給出了OPS生成MPI并行代碼在用滿20核時(shí)的性能?？梢钥吹?，OPS生成的GPU并行代碼相比手工實(shí)現(xiàn)的GPU并行代碼性能更高，在2M和4M網(wǎng)格規(guī)模下分別達(dá)到手工GPU并行代碼的1.6倍和1.5倍。此外，在2M和4M網(wǎng)格規(guī)模下，OPS生成GPU并行代碼相對(duì)于其MPI并行代碼分別獲得了2.2倍和2.0倍的性能加速。

Figure 6 Performance comparison of OPS generated GPU code,manual GPU code and OPS generated CPU code圖6 OPS生成的GPU代碼和手工GPU代碼以及OPS生成的CPU代碼的性能對(duì)比

對(duì)于HNSC軟件，OPS代碼的性能在純MPI并行、純OpenMP并行、MPI/OpenMP 2級(jí)混合并行以及GPU并行4種并行模式下相較于手工并行代碼能夠分別體現(xiàn)出一定或明顯的優(yōu)勢(shì)，原因包括2個(gè)方面：一方面，OPS自動(dòng)并行和手工并行化是在同一基準(zhǔn)代碼上進(jìn)行的，基準(zhǔn)代碼已經(jīng)包含了如訪存優(yōu)化、數(shù)學(xué)計(jì)算優(yōu)化等措施；另一方面，OPS中的并行代碼生成過程實(shí)現(xiàn)了一系列高級(jí)的性能優(yōu)化。例如，HNSC采用高精度差分格式，差分計(jì)算所需的通信緩沖區(qū)較寬，通信開銷對(duì)性能影響較大，而OPS在MPI并行時(shí)除了對(duì)block進(jìn)行任務(wù)劃分并平衡負(fù)載以外，還從跨循環(huán)的全局層面分析數(shù)據(jù)依賴關(guān)系，以使MPI通信的數(shù)量和大小降到最低。OPS中還采取了一種遲滯執(zhí)行技術(shù)(Lazy Execution Technique)來優(yōu)化全局消息通信[9]。同時(shí)，OPS在代碼生成時(shí)會(huì)從循環(huán)展開、線程調(diào)度和數(shù)據(jù)訪問方面進(jìn)行OpenMP多線程性能優(yōu)化。此外，OPS在進(jìn)行GPU并行時(shí)，提供了對(duì)Cache數(shù)據(jù)載入、紋理內(nèi)存使用等的描述，OPS可以基于這些描述信息來優(yōu)化GPU上的存儲(chǔ)層次訪問，而手工版的GPU代碼并沒有針對(duì)存儲(chǔ)層次訪問進(jìn)行深入的優(yōu)化。

5 結(jié)束語

本文采用OPS領(lǐng)域特定語言，對(duì)高階精度計(jì)算流體力學(xué)軟件HNSC進(jìn)行了重構(gòu)，并面向CPU和GPU平臺(tái)自動(dòng)生成了多個(gè)版本的并行代碼。在一個(gè)有雙Intel Haswell CPU和1塊NVIDIA Tesla K80 GPU的服務(wù)器上的性能測(cè)試表明，基于OPS自動(dòng)生成的并行代碼性能與手工并行代碼的性能可比，有些情況下其性能甚至優(yōu)于手工并行代碼的。并且OPS自動(dòng)生成的GPU并行代碼相對(duì)于其CPU并行代碼有明顯的性能加速，表明了基于OPS方法自動(dòng)并行代碼生成的跨平臺(tái)性能可移植性。這些結(jié)果說明，使用OPS等領(lǐng)域特定語言進(jìn)行計(jì)算流體力學(xué)并行軟件開發(fā)的效率高，并且面向多平臺(tái)的性能可期，是面向E級(jí)計(jì)算時(shí)代一個(gè)很有前景的方向。