格子Boltzmann方法在GPU平臺下對多孔介質(zhì)流動的模擬

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

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

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

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

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

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

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

1 格子Boltzmann方法

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

1）對fi進行流動演化；

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

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

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

5）利用fi來計算mi；

6）計算u**=u*+adt/2=u+adt；

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 滲透率的計算

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

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

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

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

4 結(jié)論

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

參考文獻：

[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-），男，漢族，湖北仙桃，碩士，研究方向：計算流體力學。