李華兵， 赫軼男， 張乾毅， 韋華建， 韋國柱

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

近年來，在物理學領域中晶格玻爾茲曼法作為一種流體力學常用方法，越來越受到國內外學者的青睞。傳統(tǒng)基于CPU上的線性運算的運算速度問題也逐漸暴露出來，作為發(fā)明GPU的領頭公司NVIDIA亦認識到這一問題，其在2007年便推出了基于GPU上的并行計算架構(CUDA)，于是便開啟了CUDA渲染時代。

在物理學中，晶格玻爾茲曼方法在計算流體動力學方面的模擬是一個非常好的工具。迄今為止，LBM已成功地應用于多種復雜流體系統(tǒng)，如多孔介質流[1]、多相流[2]、反應擴散流[3]、血液流[4]等。將計算流體動力學與GPU數(shù)值模擬相結合亦取得了眾多成果，如Tolke等[5]利用CUDA的并行優(yōu)勢實現(xiàn)了LBM中的D3Q15模型，張云等[6]實現(xiàn)了多松弛模型的計算加速。這2種方法雖都將CUDA應用在LBM上，但由于LBM規(guī)則格子模型處理難度較大，無法充分展示CUDA架構的加速效果。鑒于此，以D2Q9模型為例，通過改善內存訪問以及數(shù)據(jù)傳輸?shù)姆椒?，進一步提高LBM的計算速度。

1 GPU與CUDA

1.1 GPU與CPU的區(qū)別

GPU與CPU二者之間雖不同，但也存在聯(lián)系。2種器件的相同點是它們都是處理單元；不同點是CPU是“核心的”，而GPU主要用于“圖像”處理。二者具有不同的架構，GPU采用了數(shù)量眾多的計算單元和超長的流水線，但只有非常簡單的控制邏輯，并省去了儲存單元，而CPU不僅被儲存單元占據(jù)了大量空間，且還有復雜的控制邏輯和諸多優(yōu)化電路。雖然在邏輯控制與線性運算方面較好，但相比之下，GPU更擅長大規(guī)模的并行計算。

在最初的設計上，由于CPU僅具備單個DRAM內存塊，其線程之間切換極為麻煩，僅能做到線性切換，而GPU具有多個存儲器控制單元，進而多個線程之間可實現(xiàn)相互的快速轉換。因此，相比于任務并發(fā)型優(yōu)化的CPU，GPU更擅長數(shù)據(jù)的并發(fā)型優(yōu)化。本計算采用Inter(R) Xeon(R) CPU E5-2620 v4，該CPU作為專業(yè)的服務器核心，含有8個內核數(shù)，16個線程，主頻為2.1 GHz。它可通過多種散熱管理功能保護處理器包和系統(tǒng)，以免出現(xiàn)過熱故障。帶有擴展頁表(EPT)的VT-x，也稱為二級地址轉換(SLAT)，可為需要大量內存的虛擬化應用提供加速。至強E5-2620 v4 2.1 GHz 處理器集多種功能于一身，可滿足幾乎所有數(shù)據(jù)處理量大的基礎架構應用程序的需求。而采用的專業(yè)級顯卡型號為Nvidia Quadro GP100，該顯卡具有較高的性價比，擁有16 GiB的顯存內存及3 584個CUDA核心。一個專業(yè)級顯卡所具備的CUDA核心直接反映了計算速度的快慢。由于要應用于LBM計算上，需達到雙精度浮點計算，該型號專業(yè)級顯卡用作LBM運算為較優(yōu)選擇。

GPU的作用已經不僅局限于圖形處理與圖像識別。經實驗證明，在高精度的浮點運算與并行計算領域上，相比于CPU，GPU甚至可以提供數(shù)十倍乃至于上百倍速度比。

1.2 CUDA并行計算架構

在GPU數(shù)值計算領域，通用計算軟件體系主要有OpenMP、CUDA、MPI。本研究采用的是CUDA，其是由NVIDIA公司所推出的并行運算開發(fā)平臺，且隨著時間的推演CUDA從起初僅可使用C語言進行程序編譯，發(fā)展到現(xiàn)在的CUDA 3.0可支持C++以及FORTRAN，展示了其所具有的高兼容性。關于GPU并行運算加速LBM的計算在國內外有很多例子，如Li等[7]將LBM 演化過程映射為對2D紋理單元的光柵化和幀緩存的處理，并采用Nvidia GeForce4 Ti4600顯卡，首次實現(xiàn)了在GPU上加速LBM計算，較傳統(tǒng)的CPU獲得了15.3的加速比。隨后Fan等[8]采用相似編程方法將LBM運行在32個GPU節(jié)點的GPU集群，與CPU集群相比，獲得了21.4的加速比。

并行計算是指同時使用多種計算資源解決計算問題的過程，是提高計算機系統(tǒng)計算速度和處理能力的一種有效手段。并行計算具有分明的并行層次，首先是作為底層的單核指令級并行，即在單個處理器能讓多條指令命令同時并行執(zhí)行；在這之上的是多核并行，由單個芯片構成，其允許多個處理器核心上的多線程之間相互并行；其次是多處理器并行，多個處理器之間通過集成安裝的方式，使線程與進程之間到達并行；最頂為集群或分布式并行，作為一個獨立的并行計算平臺，亦可稱為節(jié)點，各個節(jié)點之間為達到數(shù)據(jù)之間的共享，借助網(wǎng)絡傳輸?shù)姆绞阶罱K達到集群或分布式的并行。

2 晶格玻爾茲曼法

用單粒子分布函數(shù)fi代替布爾變量，則LBM演化方程為

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

(1)

其中：fi(x,t)是在x位置t時刻具有ei速度的粒子的分布函數(shù)；Ω(fi)為碰撞因子，表示碰撞對于fi的影響。

之后引入單弛豫時間近似代替碰撞因子項[9]，碰撞項被簡化為

(2)

其中τ為弛豫時間。因此單馳豫LBM方程可寫為

(3)

平衡分布函數(shù)是將LBM應用在不同問題中的關鍵要素，只要提供一個合適的平衡分布函數(shù)，不同的物理問題都可用LBM來解決。平衡分布函數(shù)的一般形式為

(4)

其中：u為宏觀流動速度矢量；A、B、C、D為常數(shù)，需要通過守恒定理加以確定；φ為標量，如密度(ρ)、溫度或組分濃度等，其等于所有分布函數(shù)的總和。

二維九速(D2Q9) 晶格玻爾茲曼模型如圖 1 所示，其邊界不存在六角形格子中的階梯狀邊界問題，因而得到了廣泛使用。利用Chapman-Enskog多尺度展開[10]，以及質量與動量守恒，得到局域平衡分布函數(shù)：

質量守恒時，

(5)

動量守恒時，

(6)

(7)

圖1 D2Q9晶格玻爾茲曼模型

3 基于CUDA的程序算法優(yōu)化加速

該程序的的執(zhí)行步驟為：

1) GPU開辟內存空間，調用CUDA函數(shù)中的cudaMalloc()函數(shù)。

2) 調用CUDA架構中的cudaMemcpy()函數(shù)，將數(shù)據(jù)傳遞到GPU上，且此函數(shù)中的最后一個參數(shù)應為cudaMemcpyHostToDevice。

3) 在GPU并行執(zhí)行，設計了KernelCollide、KernelStream、KernelCalculate三個核函數(shù)，分別代表每個格子進行碰撞、遷移流動、計算每個格子點的宏觀量。本計算采用CUDA中的共享內存，對KernelStream、KernelCalculate兩個核函數(shù)進行優(yōu)化。

4) 結束迭代循環(huán)，調用CUDA架構中的的cudaMemcpy()函數(shù)，將GPU的數(shù)據(jù)傳回CPU，此函數(shù)中的最后一個參數(shù)應為cudaMemcpyDeviceToHost。

基于CUDA上的LBM其最大的問題在于內存的訪問以及數(shù)據(jù)的傳輸，覃章榮等[11]已證明通過使用CUDA架構中的全局內存、共享內存、紋理內存能夠對LBM進行加速，且共享內存優(yōu)化效果最佳。本研究即針對共享內存的使用再進行優(yōu)化。

3.1 共享內存

CPU與GPU均采用動態(tài)隨機存取儲存器DRAM。與CPU上不可編譯的一級緩存與二級緩存不同的是，GPU上的內存設備有很多種：寄存器、共享內存、全局內存、本地內存、常量內存以及紋理內存，如圖2所示。每個線程都對應一個私有的本地內存；每個塊中都有一個共享內存，且該共享內存對整個塊中的線程可見，同時不同塊的線程無法訪問另一個塊中的共享內存；每個線程均可以訪問網(wǎng)格中的全局內存、常量內存以及紋理內存。

圖2 CUDA內存模型

共享內存具有低延遲、高帶寬的特點。共享內存作為GPU上的一個關鍵部分，其使同一個線程塊中的線程可以相互協(xié)同，最終實現(xiàn)內存利用最大化，進而降低全局內存讀取的延遲。在核函數(shù)中使用_shared_關鍵字修飾變量，定義一維、二維甚至三維數(shù)組作為共享內存緩沖區(qū)。由于共享內存有大小的限制，不可設置太大的數(shù)組。在使用共享內存時還需用到線程同步函數(shù)_syncthreads()，該函數(shù)的調用將確保線程塊中的每個線程都執(zhí)行完_syncthreads()前的語句，才會執(zhí)行下一條語句，最終保證其他線程都能執(zhí)行完。在VS編譯器可能會出現(xiàn)_syncthreads()未定義標識符的警告提醒，這是由于VS編譯器的原因所致，編譯無法識別該函數(shù)，但是nvcc可以識別。

3.2 程序設計

在整個程序中預先定義了一個Lattice_Str結構體，用于存放9個分布函數(shù)、x與y方向上的速度以及格子點的密度等數(shù)據(jù)。再定義一個parameter結構體，為之后格子點遷徙流動做鋪墊。

struct Lattice_Str

{

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

Cuda_Vector m_v; //包含x與y方向上的速度

double m_dDen; // 密度

…

};

struct parameter

{

int iPrjx[9]; //代表 {0,1,-1,0,0,1,-1,1,-1};

int iPrjy[9];//代表 {0,0, 0,1,-1,-1, 1,1,-1};

…

}

針對3個內核函數(shù)分別進行討論。

1) KernelCollide格子碰撞內核函數(shù)。由于在LBM中粒子碰撞過程主要涉及到有關計算，用式 (6) 計算平衡分布函數(shù)，放入_device_關鍵字修飾的成員函數(shù)，該函數(shù)僅可在GPU上運行。

_ device _ void feq(double *dR,int Order,double Den, double Velx,double Vely,parameter *p)

{

double dfeq;

double dDotMet;

dDotMet=p→iPrjx[Order]*Velx+ p→iPrjy[Order] * Vely;

dfeq=Den*p→dCoe[Order]*(1+3*dDotMet + 4.5*dDotMet*dDotMet-1.5*(Velx*Velx+Vely*Vely));

*dR=dfeq;

}

由于GPU具有優(yōu)秀的計算能力，計算速度得到提升。

2) KernelStream格子遷徙流動內核函數(shù)。啟動KernelStream<<>>(…)核函數(shù)，在本計算中先采用500×100的格子，即50×10個Block，且每個Block有10×10個Thread。而設定共享內存緩沖區(qū)為10×10。每個格子點對應一個結構體，該結構體包含9個分布函數(shù)。

struct SharedMemory_data

{

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

};

_ shared _ SharedMemory _ data

La_SharedM[10][10];

計算二維網(wǎng)格點上的索引i=threadIdx.x+blockIdx.x*blockDim.x;j=threadIdx.y+blockIdx.y*blockDim.y。建立一個循環(huán)，將原矩陣格子上8個方向的分布函數(shù)傳給_shared_修飾的共享內存緩沖區(qū)La_SharedM[tthreadIdx.y][threadIdx.y]，按threadIdx.x維劃分，即threadIdx.y固定而threadIdx.x連續(xù)變化，這樣在一個線程塊即可進行一行一行的訪問，從而避免了存儲體沖突。在傳給共享內存緩沖區(qū)時，同時做到傳播到相鄰的格子點，即(0,0)(1,0)(-1,0)(0,1)(0,-1)(1,1)(-1,1)(-1,-1)(1,-1)這幾個方向，需注意的是要添加線程同步函數(shù)_syncthreads()。

_ global _ void KernelStream(…，parameter* p)

{

…

_ shared _ parameter st_p;

st_p=*p;

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

{

La_SharedM[threadIdx.y][threadIdx.x].m_df[k]=Lattice_Str[i- st_p.iPrjx[k]][j-st_p.iPrjy[k]].m_df[k];

}

_ syncthreads();

}

最終將數(shù)據(jù)從共享內存中傳回到原矩陣中。

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

{

Lattice_Str[i][j].m_df[k]=La_SharedM[threadIdx.y][threadIdx.x].m_df[k];

}

3)KernelCalculate計算格子點宏觀量核函數(shù)。啟動KernelCalculate<<>>(…)核函數(shù)，設置線程塊與線程個數(shù)。創(chuàng)建共享內存緩沖區(qū)，將x與y方向上的速度以及格子點的密度進行初始化。

struct SharedMemory_data

{

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

double m_dDen; //格子點密度

double m_vx; //x方向速度

double m_vy; //y方向速度

…

};

_ shared _ SharedMemory_data

Lattice_Str[[10][10];

Lattice_Str[threadIdx.y][threadIdx.x].m_vx=0;

Lattice_Str[threadIdx.y][threadIdx.x].m_vy=0;

Lattice_Str[threadIdx.y][threadIdx.x].m_dDen=0;

類似KernelStream格子遷徙流動內核函數(shù)，計算二維網(wǎng)格點上的索引，通過建立循環(huán)將原矩陣9個方向的分布函數(shù)以及x與y方向上的速度拷入共享內存緩沖區(qū)，再將數(shù)據(jù)從共享內存?zhèn)骰卦仃?，同時也要針對可能產生的錯誤進行if語句判斷，不能一味地為了計算效率而忽略錯誤判斷。

3.3 計算結果及討論

本計算采用的是泊松流模型，在對程序進行優(yōu)化加速的同時，也能保證程序結果的正確性。泊松流水平方向速度分布如圖 3 所示。從圖3可看出，泊松流的水平方向流速為拋物線。本模型進行了10 000步計算(流場還未達到穩(wěn)定)，采用不同共享內存網(wǎng)格大小，與不采用共享內存相比，計算加速比，得到不同的共享內存加速比如圖 4 所示。從圖 4 可看出，最高加速比可達 2.15，隨著共享內存的增加，加速比有所下降，這是因為共享內存的分配是有限的。本計算采用的計時函數(shù)是CPU上的clock()函數(shù)，與傳統(tǒng)上使用cudaEventRecord()函數(shù)單一測核函數(shù)計算時間不同，要測加速比也需要考慮創(chuàng)建內存時間以及釋放內存的耗時，即整個項目的時間，不能只考慮核函數(shù)運行所節(jié)省的時間。

圖3 泊松流水平方向速度分布

圖4 不同共享內存的加速比

4 結束語

針對基于GPU上的CUDA架構，對LBM數(shù)值模擬程序加速進行討論與分析，CUDA架構中的共享內存LBM計算，CUDA的并行計算與LBM天然的并行性能夠很好結合。該程序中由于有針對可能出現(xiàn)的越界現(xiàn)象的判斷，這會影響計算效率，需后續(xù)進行一步優(yōu)化。