張乾毅， 韋華健， 赫軼男， 李華兵

(桂林電子科技大學(xué) 材料科學(xué)與工程學(xué)院，廣西 桂林 541004)

隨著處理器時鐘頻率極限及“功耗墻”等問題的出現(xiàn)，單核解決方案被摒棄，基于多核架構(gòu)的圖像處理單元(GPU)逐漸出現(xiàn)在人們的視野中。2007年，英偉達公司為GPU增加了一個易用的編程接口——統(tǒng)一計算架構(gòu)(compute unified device architecture，簡稱CUDA)[1]。它以標準C為基礎(chǔ)進行擴展，使程序員可以在C、C++及Fortran等開發(fā)環(huán)境進行并行程序編寫，將利用GPU進行并行計算的門檻大大降低。近年來，人工智能在深度學(xué)習(xí)的不斷推動下高速發(fā)展，這離不開GPU的強大計算能力和優(yōu)秀的GPU編程環(huán)境。CUDA技術(shù)已經(jīng)獲得廣泛關(guān)注，并被應(yīng)用到流體力學(xué)、生物計算、氣象分析、金融分析等眾多領(lǐng)域。

晶格玻爾茲曼方法(LBM)是一種介觀尺度的模擬方法，在流體力學(xué)方面的模擬中被廣泛應(yīng)用，如多孔介質(zhì)流[2]、血液流[3]、多相流[4]淋巴流[5]等。LBM物理背景清晰，易于并行計算，邊界條件處理簡單，因而十分適用于大規(guī)模GPU并行計算。如T?lke等[6]使用CUDA實現(xiàn)了LBM多孔介質(zhì)流動模型的計算加速。鄭彥奎等[7]實現(xiàn)了LBM方腔模型的計算加速。

鑒于此，通過線性尋址和下標尋址方法實現(xiàn)了LBM模型中的泊松流算例。將程序分別用CPU和GPU進行計算，比較2種尋址方法分別在2種處理單元上的計算時間，并算得加速比。

1 GPU與CUDA

基于設(shè)計目標的差異，GPU在設(shè)備架構(gòu)上與CPU存在很大區(qū)別。CPU需要很強的通用性來處理各種不同的數(shù)據(jù)類型，同時邏輯判斷又會引入大量的分支跳轉(zhuǎn)和中斷處理。因此，CPU中分布著多樣的計算、控制單元和占據(jù)大量空間的存儲單元。而GPU的設(shè)計主要是用來解決大量邏輯上相對簡單的任務(wù)，因而GPU采用了數(shù)量眾多的核心和算術(shù)邏輯單元(ALU)，其數(shù)量遠超CPU，通常CPU的核心數(shù)不會超過2位數(shù)，但GPU只配備了簡單的控制邏輯和簡化了的存儲單元。

CUDA編程模型將CPU作為主機端，GPU作為設(shè)備端，二者各自擁有相互獨立的存儲地址空間：主機端內(nèi)存和設(shè)備端顯存。一個典型的CUDA程序包括并行代碼和與其互補的串行代碼。由CPU執(zhí)行復(fù)雜邏輯處理和事務(wù)處理等不適合數(shù)據(jù)并行的串行代碼，并調(diào)用GPU計算密集型的大規(guī)模數(shù)據(jù)并行計算代碼，即內(nèi)核函數(shù)。當設(shè)備端開始執(zhí)行內(nèi)核函數(shù)時，設(shè)備中會產(chǎn)生大量的線程，線程是內(nèi)核函數(shù)的基本單元，負責(zé)執(zhí)行內(nèi)核函數(shù)的指定語句[8]。CUDA程序模型如圖1所示，由一個內(nèi)核啟動所產(chǎn)生的所有線程統(tǒng)稱為一個線程網(wǎng)格，同一網(wǎng)格中的所有線程共享相同的全局內(nèi)存空間。一個網(wǎng)格又由多個線程塊構(gòu)成，同一線程塊內(nèi)的線程協(xié)作可通過同步和共享內(nèi)存的方式實現(xiàn)，不同塊內(nèi)的線程不能協(xié)作。每個GPU設(shè)備包含一組流多處理器(SM)，一個SM又由多個流處理器(SP)和一些其他硬件組成。CUDA程序執(zhí)行過程中會把線程塊中的所有線程以線程束(Warp，含32個線程)為執(zhí)行單位分配給SM中的不同SP來進行計算。

圖1 CUDA程序模型

GPU采用了內(nèi)存分層架構(gòu)，并通過多級緩存機制來提高讀寫效率。CUDA內(nèi)存模型如圖2所示，GPU中主要的內(nèi)存空間包括寄存器、共享內(nèi)存、全局內(nèi)存、常量內(nèi)存、紋理內(nèi)存和本地內(nèi)存。每個SM都擁有獨立的一級緩存，所有SM共享二級緩存。GPU致力于為每個線程分配真實的寄存器，當寄存器使用到達上限時，編譯器會將數(shù)據(jù)放在片外的本地內(nèi)存中。共享內(nèi)存是可受程序員控制的一級緩存，每個SM中的一級緩存與共享內(nèi)存公用同一個內(nèi)存段[9]，它僅對正在執(zhí)行的線程塊中的每個線程可見。

圖2 CUDA內(nèi)存模型

2 晶格玻爾茲曼方法(LBM)

McNamara等[10]用單粒子分布函數(shù)fi取代了格子氣細胞自動機中的布爾變量，其LBM演化方程為

fi(x+eiδt,t+δt)-fi(x,t)=Ω(fi)，

(1)

其中：i為微觀速度的索引，對于D2Q9模型，i取值1～9；fi(x,t)為在x位置t時刻具有ei速度的粒子的分布函數(shù)；δt為時間增量；Ω(fi)為碰撞因子，表示碰撞對于fi的影響。文獻[11-13]采用單弛豫時間來替代碰撞因子項。繼而，LBM方程可寫為

(2)

其中：fi(eq)為局域平衡分布函數(shù)；τ為弛豫時間。求出粒子分布函數(shù)后，則宏觀密度、流體的流速分別為

(3)

(4)

采用圖3所示的D2Q9模型(D指維度，Q指粒子運動方向總數(shù))模型[14]進行計算，其平衡分布函數(shù)具體定義為

圖3 D2Q9模型

(5)

其中：c為基準速度，通常情況下取1；ωi為各方向的權(quán)重系數(shù)。當i=0時，ωi=4/9；當i=1,2,3,4時，ωi=1/9；當i=5,6,7,8時，ωi=1/36。ei的形式為

3 算法實現(xiàn)及優(yōu)化

一個典型的LBM計算案例主要分為3個步驟：1)聲明變量并進行分布函數(shù)的初始化；2)進行格點的碰撞、遷徙流動和邊界處理，并進行宏觀量的計算；3)進行穩(wěn)定性條件判斷，決定是否結(jié)束步驟2)的迭代。由文獻[15]知，步驟2)的迭代計算占整個計算時間的98%，因此該步驟應(yīng)在GPU上進行并行加速。本計算以淋巴管為物理背景，將淋巴管內(nèi)流動的淋巴液視作等截面圓形管道內(nèi)的泊松流，進行模擬計算。嘗試用線性尋址和下標尋址2種不同的尋址方法進行計算，比較這2種尋址方法分別在GPU與CPU上執(zhí)行計算時間的差異。

3.1 線性尋址程序設(shè)計

LBM的基本變量是分布函數(shù)，進而由分布函數(shù)計算求得格子點的密度和x、y方向上的速度等宏觀量。在串行程序的D2Q9模型中，分布函數(shù)通常被設(shè)為三維數(shù)組df[i][j][k]并存儲在主存中，i、j表示整個線程網(wǎng)格中格點的坐標(0≤i≤width，0≤j≤height；width、height均為固定值，分別表示線程網(wǎng)格的寬度和高度)，k為格點運動方向數(shù)，k=0,1,…,8。

在設(shè)備端聲明指針double*df，調(diào)用cudaMallocManaged()函數(shù)在設(shè)備端開辟大小為width*height*9*sizeof(double)的一維線性內(nèi)存空間，用于存放每個格點上9個方向的分布函數(shù)，并將在顯存上獲得的內(nèi)存空間首地址賦值給df。使用cudaMallocManaged()函數(shù)開辟的存儲空間，無論是在串行代碼中還是并行代碼中，都可使用這塊內(nèi)存，因此只需定義一個指針即可在主機和設(shè)備端通用。繼續(xù)開辟多個設(shè)備端內(nèi)存空間，用于存放其他計算變量的數(shù)據(jù)。指針double *dd指向存儲格子點密度的內(nèi)存空間，指針double *dvx、*dvy分別指向存放x和y方向上速度的內(nèi)存空間。

基于上述在全局內(nèi)存中使用三維數(shù)組來儲存分布函數(shù)及其他變量的方式，設(shè)計了3個在CPU端執(zhí)行的collide()、stream()、calculate()函數(shù)，分別代表格點的碰撞、遷徙流動和宏觀量的計算等步驟，并將其命名為方案一。方案二則是將迭代計算過程放在GPU端并行執(zhí)行，對應(yīng)設(shè)計了addKernelCollide()、addKernelCopy()、addKernelStream()、addKernelCalculate()四個內(nèi)核函數(shù)來實現(xiàn)。相較于方案一中執(zhí)行各個步驟都需要嵌套for循環(huán)來遍歷讀寫每個格點上的數(shù)據(jù)，方案二則通過內(nèi)核函數(shù)用GPU映射出大量線程，同時對每個格點數(shù)據(jù)進行操作。

首先，調(diào)用cudaMallocManaged()函數(shù)，在GPU端為分布函數(shù)、密度等變量開辟內(nèi)存空間，再定義一個parameter結(jié)構(gòu)體，為之后格子點遷徙流動做鋪墊。對4個內(nèi)核函數(shù)進行如下討論：

1)碰撞步：啟動addKernelCollide<<>>(df,iSol,…,dev_p)內(nèi)核函數(shù)，其中g(shù)rid為線程網(wǎng)格中線程塊數(shù)，將其設(shè)為width；block為每個線程塊中的線程數(shù)，設(shè)為9*height；iSol為用于判斷格點是否位于邊界上的變量；dev_p為parameter結(jié)構(gòu)體變量。聲明變量i、j、k并賦值，i=blockIdx.x，j=threadIdx.y，k=threadIdx.x。在函數(shù)體內(nèi)分別計算出線程訪問分布函數(shù)一維數(shù)組的下標索引id=adr(i,j,k)=k+(j+i* height)*9和其他變量的下標索引ind=adr(i,j)=j+i* i*height。通過一次if判斷，保留格子邊界以外的所有格子點進行計算，每個線程利用式(2)計算對應(yīng)格子點的分布函數(shù)：

df[id]=df[id]-(df[id]-

feq(k,dd[ind],dvx[ind],dvy[ind]))/Tau。

feq()為平衡分布函數(shù)，Tau為弛豫時間。

2)流動遷徙步：啟動addKernelCopy<<>(df,dfbak)內(nèi)核函數(shù)，聲明變量并賦值，計算出線程訪問的下標索引，方法同上。將df[id]中的數(shù)據(jù)傳遞給dfbak[id]。啟動addKernelStream<<>>(df,dfbak,iSol,…)內(nèi)核函數(shù)，將dfbak[id]的值賦給df[id]，此時dfbak[id]中

id=k+[(j-dev_Prjy[k])+(i-

dev_Prjx[k])*height)]*9。

3)計算宏觀量步：啟動addKernelCalculate<<>>(df,iSol,…)內(nèi)核函數(shù)，聲明變量并放入共享內(nèi)存中。

__shared__ double sdf[9];//9個分布函數(shù)

__shared__ double ddd;//格子點密度

__shared__ double vx;//x方向速度

__shared__ double vy;//y方向速度

……

for(k=0;k<9;k++)

{

sdf[k]=df[adr(i,j,k)];

}

利用sdf[k]計算得到ddd、vx、vy的值，最后將共享內(nèi)存中的數(shù)據(jù)傳回主存中。

3.2 下標尋址程序設(shè)計

基于之前的線性尋址方法對程序做以下修改。

1)聲明指針double *df0、double **df1、double ***df，在GPU上開辟3個內(nèi)存空間，并將首地址分別賦值給df0、df1、df，具體操作如下：

cudaMallocManaged(void(**)&df0,width*height*9*sizeof(double));

cudaMallocManaged(void(**)&df1,width*height*sizeof(double*));

for(i=0;i

{

for(j=0;j

{df1[i*height+j]=&df0[(i*height+j)*9];}

}

cudaMallocManaged(void(**)&df,width*height*sizeof(double**));

for(i=0;i

{df[i]=&df1[i*height];}

2)聲明指針double *dd0、double **dd，在GPU上開辟2次內(nèi)存空間，并將首地址分別賦值給dd0、dd，具體操作如下：

cudaMallocManaged(void(**)&dd0,width*height*sizeof(double));

cudaMallocManaged(void(**)&dd,width*sizeof(double*));

for(i=0;i

{dd[i]=&dd0[i*height];}

3)其他宏觀量修改方法同上。

在各計算步驟中，用多維數(shù)組表示不同格點位置上的分布函數(shù)和其他宏觀量，形如df[i][j][k]、dd[i][j]、dvx[i][j]、dvy[i][j]。將方案一、二按上述內(nèi)容進行修改，并分別命名為方案三、四。

4 計算結(jié)果分析

采用含有8個Nvidia Quadro GP100顯卡集群的服務(wù)器完成計算，GP100顯卡擁有3 584個CUDA并行計算處理核心，處理雙精度浮點數(shù)的能力為5.2 TFlop/s，搭配16 GiB HBM2顯存，理論帶寬高達717 GiB/s。同時服務(wù)器搭載了Intel(R) Xeon(R)CPU E-52620 v4處理器，核心數(shù)為8個，主頻為2.1 GHz，編譯環(huán)境為CUDA10.0，運行環(huán)境為Linux Ubuntu。

表1 4種泊松流模擬方案運行時間

以平面泊松流為算例，在保證計算結(jié)果正確性的前提下驗證2種尋址方法對程序計算時間的影響。圖4為4種方案迭代4 000步時的模擬結(jié)果，4種方案結(jié)果一致，表明了程序修改的正確性。表1為迭代4 000步時，4種方案模擬的計算時間。從表1可看出，基于線性尋址方式的方案一、二，GPU相對CPU的加速比約為71倍，而基于下標尋址方法的加速比只有約25倍；方案一、三都用CPU完成迭代計算步，方案三的計算時間減少了將近三分之二，因而使用下標尋址的方式更適合CPU端的迭代計算；方案二、四計算時間無太大變化，表明下標尋址對GPU端的并行計算無明顯幫助。

圖4 4種泊松流模擬方案結(jié)果

5 結(jié)束語

采用CUDA實現(xiàn)了LBM的泊松流算例。設(shè)計了線性尋址、下標尋址2種尋址方法，并實現(xiàn)了這2種尋址方法分別在主機端和設(shè)備端上的4種LBM計算方案。4種方案計算結(jié)果一致，表明了程序設(shè)計的正確性。用線性尋址方法獲得了71倍的加速比，表明用CUDA對LBM計算效率的提升有很大幫助，也展現(xiàn)了GPU在大規(guī)?？茖W(xué)計算中的巨大潛力。