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

(國防科技大學計算機學院量子信息研究所兼高性能計算國家重點實驗室,湖南 長沙 410073)

1 引言

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

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

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

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

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

2 OPS簡介

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

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

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

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

2.2 自動代碼生成

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

2.3 OPS的優(yōu)點

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

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

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

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

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

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

基于圖1所示的區(qū)域劃分，在每一個網(wǎng)格區(qū)域上，HNSC軟件分別對流場數(shù)據(jù)進行求解，軟件的執(zhí)行過程如圖2所示。首先初始化MPI并行環(huán)境，Grid函數(shù)讀入網(wǎng)格，在進行計算前通過Init函數(shù)對流場的基本未知量進行初始化，初始化之后的基本未知量通過Solvec函數(shù)組裝成方便控制方程求解的守恒變量。隨后進入主循環(huán)迭代，主循環(huán)中Compute_dt函數(shù)基于場變量的當前值計算時間步長，RK4函數(shù)是HNSC軟件計算的核心，包含4階龍格庫塔算法的全部計算，并在一個時間步內(nèi)更新場變量的數(shù)值。4階龍格庫塔計算之后，進行收斂檢驗，進而計算下一個時間步的步長，并將迭代時間值向前推進，如果達到最大迭代步數(shù)，計算結(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計算控制方程中的通量項，F(xiàn)lux_der函數(shù)計算通量的偏導數(shù)。4階龍格庫塔算法的4階計算具有相似的流程，為了簡潔圖2中僅展示了第1階計算的過程，值得注意的是在4階龍格庫塔算法的每一階計算中，都包含2次邊界條件施加函數(shù)BC的執(zhí)行，這對后面將軟件和算例基于OPS實現(xiàn)是十分重要的。

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

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

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

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

(3)定義計算格式；

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

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

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

在添加相應(yīng)的頭文件之后，可以使用OPS提供的編程接口聲明一個名為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的求解最多針對三維空間來求解)，針對HNSC軟件，grid_dim設(shè)為3。字符型參數(shù)“grid_block”是該block的標簽名，用于代碼轉(zhuǎn)換、編譯以及數(shù)據(jù)I/O時的顯示和輸出。

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

HNSC作為CFD軟件所求解的基本未知量包括密度、速度的XYZ3個坐標軸分量、壓力、溫度和內(nèi)能，分別使用符號rho、u、v、w、p、Temp、e來表示。并且軟件還將涉及CFD求解的過程變量，如通量項等。這些變量在軟件代碼中作為雙精度浮點型數(shù)組來計算。基于OPS來實現(xiàn)時，也要定義對應(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”);

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

求解HNSC軟件時，采用守恒形式的控制方程，解向量、3個方向的通量和源項分別用符號Uvec、Fvec、Gvec、Hvec和Jvec來表示，對于每一個網(wǎng)格點，這些變量是一個長度nvar=5的列向量，因此對解向量應(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”);

因為MPI并行是OPS自動實現(xiàn)的，不再需要關(guān)注軟件實現(xiàn)中的任務(wù)劃分，所以在定義數(shù)據(jù)集的時候，數(shù)據(jù)集的規(guī)模nx、ny、nz是圖1中整個計算區(qū)域所有網(wǎng)格點的規(guī)模。OPS實現(xiàn)MPI并行時會根據(jù)數(shù)據(jù)交換區(qū)的信息halo_pos和halo_neg來定義MPI通信緩沖區(qū)。在網(wǎng)格塊的內(nèi)部，OPS將根據(jù)差分計算的循環(huán)自行定義MPI通信行為的觸發(fā)。因此，在定義好數(shù)據(jù)集和數(shù)據(jù)交換區(qū)大小之后，圖1中HNSC軟件原本負責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è)置成同一個網(wǎng)格塊的數(shù)據(jù)集，這樣便可以實現(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);

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

3.2.3 定義計算格式

OPS通過定義計算格式來描述類似于有限差分計算，這樣計算某一個網(wǎng)格點數(shù)據(jù)時，使用周圍特定相鄰網(wǎng)格點的數(shù)據(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");

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

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

OPS對軟件中遍歷網(wǎng)格點的差分計算循環(huán)使用在某個(或多個)數(shù)據(jù)集上、基于相應(yīng)計算格式的訪問來描述，用戶定義一個kernel函數(shù)來描述對計算格式的具體操作，供OPS循環(huán)接口調(diào)用，從而實現(xiàn)原代碼中的計算循環(huán)。HNSC中計算速度分量u對坐標x的偏導數(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改寫該計算循環(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的一個3維(層)循環(huán)。ops_arg_dat傳遞該循環(huán)需要操作的數(shù)據(jù)集，并規(guī)定了針對該數(shù)據(jù)集的計算格式為S3D_000或S3D_25pt，訪問方式為：OPS_WRITE(寫入)、OPS_READ(讀取)或OPS_RW(讀寫)。原循環(huán)體內(nèi)部的具體計算，即對計算格式的具體操作在函數(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提供的一個宏定義，在使用的時候要根據(jù)函數(shù)的形參位置指定序號(如OPS_ACC0)，并給出計算格式中網(wǎng)格點的相對索引(如OPS_ACC0(1,0,0))。需要特別注意的是，ops_par_loop循環(huán)的范圍loop_range為{0,nx,0,ny,0,nz}，遍歷的是整個網(wǎng)格塊，而{istart,iend,jstart,jend,kstart,kend}是HNSC手工進行區(qū)域分解的子區(qū)域的計算循環(huán)邊界。

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

3.2.5 并行代碼編譯生成

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

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

4 性能測試與分析

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

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

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

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

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

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

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

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

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

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

5 結(jié)束語

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