格子Boltzmann方法在GPU平臺(tái)下對(duì)多孔介質(zhì)流動(dòng)的模擬

摘 要：GPU有著強(qiáng)大的并行計(jì)算能力，現(xiàn)在越來越多的被人應(yīng)用到科學(xué)計(jì)算領(lǐng)域了。格子Boltzmann是一種模擬不可壓縮流動(dòng)的計(jì)算方法，由于其具有天然的并行性，因此本文利用GPU的并行計(jì)算能力，在GPU平臺(tái)上，利用格子Boltzmann模擬周期性的繞方塊多孔介質(zhì)流動(dòng)，在保證程序結(jié)果的計(jì)算精度的情況下，極大的提高了程序計(jì)算效率。同時(shí)將計(jì)算結(jié)果同有限元的結(jié)果進(jìn)行比較，也得到了較好的計(jì)算結(jié)果。

關(guān)鍵詞：GPU；格子Boltzmann；多孔介質(zhì)；并行計(jì)算

在近20年來，格子Boltzmann方法在計(jì)算流體流動(dòng)的模擬中，起到了很大的作用。但是，與此同時(shí)，格子Boltzmann（LBM）方法，在計(jì)算過程中，有著計(jì)算量大，計(jì)算時(shí)間長，內(nèi)存消耗大的特點(diǎn)。但是由于LBM方法是一種顯示計(jì)算格式，在計(jì)算的過程中，只需要利用鄰近網(wǎng)格點(diǎn)的信息，因此該方法有著良好的并行性，尤其是在GPU平臺(tái)上，可以達(dá)到很高的并行效率。

GPU 是一個(gè)大規(guī)模并行計(jì)算架構(gòu)，被廣泛的用于圖像和非圖像計(jì)算領(lǐng)域，機(jī)器學(xué)習(xí)，圖像，語音等領(lǐng)域都有著應(yīng)用。GPU的主要優(yōu)勢(shì)在于單位時(shí)間的浮點(diǎn)計(jì)算能力遠(yuǎn)遠(yuǎn)的高于CPU的浮點(diǎn)計(jì)算能力[ 1 ]。

對(duì)于多孔介質(zhì)流動(dòng)的數(shù)值模擬，大多數(shù)情況下，多孔介質(zhì)流動(dòng)，都屬于緩慢流動(dòng)，而在模擬緩慢流動(dòng)的數(shù)值方法中，LBM方法是其中一種相當(dāng)廣泛的方法。關(guān)于LBM對(duì)于多孔介質(zhì)的流動(dòng)，在很多領(lǐng)域中都有著重要的應(yīng)用：復(fù)合材料、流變學(xué)、地質(zhì)學(xué)、統(tǒng)計(jì)物理、生物科學(xué)等諸多領(lǐng)域都有著廣泛的應(yīng)用。

在本文中，我們實(shí)現(xiàn)了一個(gè)通用的LBM在GPU上運(yùn)行的程序，為了減少數(shù)據(jù)從CPU內(nèi)存上到GPU內(nèi)存上拷貝的消耗，我們將LBM算法中的每一步計(jì)算都放在GPU上面來實(shí)現(xiàn)。

本文的程序是基于nVIDIA公司CUDA C 語言實(shí)現(xiàn)的。一共分為三個(gè)部分，第一部分首先介紹了LBM方法，第二部分介紹了LBM方法在GPU平臺(tái)上的實(shí)現(xiàn)。第三部分，我們利用實(shí)現(xiàn)的程序?qū)Χ嗫捉橘|(zhì)流動(dòng)進(jìn)行模擬，并進(jìn)行相關(guān)的計(jì)算，將計(jì)算得到的結(jié)果與其他文獻(xiàn)結(jié)果進(jìn)行比較，都吻合的很好。

1 格子Boltzmann方法

那么在本文中，該方法由如下幾個(gè)步驟來實(shí)現(xiàn)：

1）對(duì)fi進(jìn)行流動(dòng)演化；

2）通過fi來計(jì)算守恒變量dr和u；

3）利用u*=u+adt/2計(jì)算u*，其中a=F/r0，這里F代表外力；

4）利用計(jì)算出來的u*和dr計(jì)算矩空間的平衡態(tài)分布函數(shù)mieq；

5）利用fi來計(jì)算mi；

6）計(jì)算u**=u*+adt/2=u+adt；

7）將mi映射到分布函數(shù)空間得到分布函數(shù)fi；

8）計(jì)算碰撞后的分布函數(shù)。

圖中的圓點(diǎn)代表每個(gè)粒子，從圖中可以看出，粒子的分布函數(shù)fi只會(huì)遷移到周圍相鄰的幾個(gè)點(diǎn)。在實(shí)際利用格子Boltzmann計(jì)算中，需要計(jì)算的介觀粒子往往很多，并且在模擬多孔介質(zhì)的流動(dòng)過程中，往往需要迭代很多的時(shí)間步，才能得到穩(wěn)態(tài)流動(dòng)的數(shù)值解。因此一個(gè)二維問題的串行程序，往往需要計(jì)算幾天甚至一個(gè)月的時(shí)間才能得到計(jì)算結(jié)果。在GPU上，最基本處理單元是SP（streaming processor），計(jì)算過程中，具體的指令和任務(wù)都是在SP上進(jìn)行處理的。GPU進(jìn)行并行計(jì)算，也就是很多個(gè)SP同時(shí)做處理。

本文程序的編寫，采用的是CUDA平臺(tái)。在CUDA平臺(tái)上，最基本的處理單元是thread，一個(gè)CUDA的并行程序會(huì)有多個(gè)threads來執(zhí)行。一定個(gè)數(shù)的threads會(huì)組成一個(gè)block，同一個(gè)block中的threads可以同步，也可以通過共享內(nèi)存（shared memory）通信。同時(shí)，多個(gè)blocks則會(huì)再構(gòu)成grid。

在GPU中，一個(gè)GPU所擁有的Thread往往有很多個(gè)。在本文所使用的Nvidia-K40的顯卡中，采用的計(jì)算架構(gòu)是Tesla Kepler，該架構(gòu)下面，所能夠調(diào)用的Thread多達(dá)兩千萬個(gè)。

因此，為了提升計(jì)算速度，本文在編寫多孔介質(zhì)流動(dòng)模擬程序的時(shí)候，采用的是一個(gè)Thread處理一個(gè)格子上面的計(jì)算。同時(shí)，因?yàn)橛?jì)算過程中，需要經(jīng)常用到格子的離散速度{ci：i=1，2，…，9}，矩陣M，矩陣M的逆矩陣，所以為了優(yōu)化程序?qū)@些量的訪問速度，本文將這些量放入常量內(nèi)存。

同時(shí)，在每一步計(jì)算過程中，都會(huì)用到矩空間的平衡態(tài)函數(shù)，為了提高函數(shù)的調(diào)用效率，本文對(duì)于矩空間的平衡態(tài)函數(shù)采用宏來進(jìn)行預(yù)定義。

由于在計(jì)算過程中，會(huì)涉及到介觀粒子的遷移，所以對(duì)于介觀粒子的分布函數(shù)，只能存放在GPU的全局內(nèi)存里面（global memory）。但是本文模擬的是多孔介質(zhì)的流動(dòng)，該問題的計(jì)算區(qū)域示意圖如下圖所示：

如上圖中，圖（a）為整個(gè)多孔介質(zhì)流動(dòng)的結(jié)構(gòu)，黑色部分為固體，空白區(qū)域?yàn)橐后w。由于多孔介質(zhì)流動(dòng)是屬于緩慢流動(dòng)，所以在本文中，可以將問題的計(jì)算區(qū)域抽象成如圖（b）所示的計(jì)算區(qū)域。邊界條件，對(duì)于固體邊界采用的是無滑移邊界條件，對(duì)于上下左右的流體邊界，采用的是周期性邊界條件。

3 數(shù)值結(jié)果

3.1 數(shù)值解的收斂性

從文獻(xiàn)[2]可以知道，格子Boltzmann程序的精度是二階精度，也就是說，多孔介質(zhì)的程序計(jì)算結(jié)果是正確的話，那么該程序的收斂階必須是二階收斂。

下面給出不同固體體積分?jǐn)?shù)的情況下的收斂階。在本文中，固體體積分?jǐn)?shù)由如下所示的式子表示：

如下表所示，給出了a=0.9025，0.01這兩種情況下的收斂階。

從上面表可以看出，不管a=0.9025是a=0.01還是的情況下，計(jì)算出來的收斂階都接近2，這也就說明了本文所編寫的程序收斂階為2，所以由此可以看出此多孔介質(zhì)流動(dòng)的模擬程序的結(jié)果與理論相符合。

3.2 滲透率的計(jì)算

同時(shí)，由于多孔介質(zhì)流動(dòng)多為緩慢流動(dòng)，對(duì)于緩慢流動(dòng)，存在達(dá)西定律，達(dá)西定律的公式為如下所示：

其中U為多孔介質(zhì)流動(dòng)的平均速度，K為多孔介質(zhì)的滲透率。

本文將格子Boltzmann方法計(jì)算出來的滲透率同文獻(xiàn)[3]中的結(jié)果進(jìn)行比較，結(jié)果圖如下所示：

從上圖可以看出，本文所采用的數(shù)值方法的GPU程序，與文獻(xiàn)[ 3 ]所采用的有限元方法的計(jì)算結(jié)果吻合的相當(dāng)?shù)暮谩?/p>

4 結(jié)論

從前面的章節(jié)可以看出，本文基于GPU平臺(tái)的多孔介質(zhì)數(shù)值模擬，所得到的結(jié)果與理論吻合的很好，與此同時(shí)，計(jì)算出來的滲透率同有限元計(jì)算得到的結(jié)果也是一樣的，與此同時(shí)，本文的程序是基于GPU平臺(tái)下的大規(guī)模并行程序，在模擬同樣的問題，同樣的網(wǎng)格尺度的情況下，本文所用的時(shí)間是普通串行程序的幾十倍，并且隨著網(wǎng)格的增多，由于GPU里面依然是一個(gè)thread處理一個(gè)格子，但是串行程序卻是一個(gè)線程處理所有的格子，相比之下，在迭代相同的步數(shù)情況下串行程序所消耗的時(shí)間會(huì)越來越多，但是本文的GPU平臺(tái)下的并行程序所消耗的時(shí)間卻不會(huì)有太大的變化。

參考文獻(xiàn)：

[1] Nit C， Itu L M， Suciu C. GPU accelerated blood flow computation using the Lattice Boltzmann Method[C].High Performance Extreme Computing Conference （HPEC），2013 IEEE. IEEE，2013：1-6.

[2] He X，Luo L S.Theory of the lattice Boltzmann method：From the Boltzmann equation to the lattice Boltzmann equation[J].Physical Review E，1997，56（6）：6811.

[3] Yazdchi K， Srivastava S，Luding S. Microstructural effects on the permeability of periodic fibrous porous media[J].International Journal of Multiphase Flow，2011，37（8）：956-966.

作者簡介：

顧超（1991-），男，漢族，湖北仙桃，碩士，研究方向：計(jì)算流體力學(xué)。