格子波爾茲曼方法的GPU并行計(jì)算.doc

摘 要本文基于圖形處理單元（GPU）的格子玻爾茲曼方法（LBM）實(shí)現(xiàn)了一些經(jīng)典流場(chǎng)的計(jì)算。程序由兩個(gè)函數(shù)迭代實(shí)現(xiàn)，碰撞函數(shù)和遷移函數(shù)。碰撞函數(shù)在共享內(nèi)存內(nèi)執(zhí)行，遷移函數(shù)在設(shè)備內(nèi)存內(nèi)執(zhí)行。為了優(yōu)化程序性能，將邊界條件在碰撞函數(shù)內(nèi)實(shí)現(xiàn)，加長(zhǎng)碰撞流程，減少設(shè)備內(nèi)存的存取開(kāi)銷。通過(guò)對(duì)同一程序設(shè)定不同迭代步數(shù)，發(fā)現(xiàn)GPU格子玻爾茲曼（LB）程序隨迭代步數(shù)增大，執(zhí)行時(shí)間在前N步隨迭代步數(shù)先保持線性增長(zhǎng)，隨后會(huì)加速增長(zhǎng)，最后回歸線性增長(zhǎng)，而相應(yīng)的CPU程序則一直保持線性增長(zhǎng)，這導(dǎo)致加速效果會(huì)從某一部開(kāi)始急劇下降。通過(guò)對(duì)對(duì)同一流場(chǎng)設(shè)定不同的格子數(shù)，發(fā)現(xiàn)GPU格子玻爾茲曼（LB）程序相對(duì)同樣格子數(shù)的CPU程序加速比保持不變。通過(guò)對(duì)同一程序設(shè)定不同的Block結(jié)構(gòu)，發(fā)現(xiàn)隨著B(niǎo)lock中Thread數(shù)的增大，迭代時(shí)間線性縮短，在Thread達(dá)到某一值時(shí)，迭代時(shí)間會(huì)劇烈減小。關(guān)鍵詞：格子波爾茲曼 并行計(jì)算 GPUAbstractBased on graphics processing unit (GPU) Lattice Boltzmann method (LBM), the paper implements some classic flow field calculation . Procedure consists of two iterative function , collision function and transfer function. Collision functions performed within the shared memory, while transfer function within the device memory . In order to optimize the application performance, boundary conditions will be performed within the collision function , which lengthen the collision stream process,thus reduce device memory access .By setting up different iteration steps, we found that for the same GPU lattice boltzmann (LB) program ,with the increase of iteration steps, the program execution time will keep linear growth in N steps, then the growth accelerate, and return linear growth at last, while the corresponding CPU program kept always linear growth, resulting in accelerate ratio will fell sharply from some step.By setting different lattice for the same flow ,we found the GPU lattice boltzmann (LB) procedures remains unchanged speed up ratio relative to the CPUs of the same lattice.Through imposing different Block structure on the same program, we found that with Threads in Blocks increasing , iteration time shortened linearly first, when the Threads reaches a certain value, the iteration time will drastically reduce.Key words: lattice Boltzmann parallel computing GPU目錄一 緒論11.引言12.格子Boltzmann方法概述23.GPU概述54.本文的主要工作6二 格子Boltzmann方法的基本原理71. Boltzmann方程72.格子Boltzmann方程93.格子Boltzmann方法的實(shí)現(xiàn)154.基本的邊界條件175.本章小結(jié)19三 CUDA技術(shù)211. CUDA 編程模型212.CUDA硬件223.優(yōu)化準(zhǔn)則234.本章小結(jié)24四 算法實(shí)現(xiàn)251.算法方面252. LBM的GPU實(shí)現(xiàn)263. 本章小結(jié)27五 結(jié)果分析281. 泊肅葉流282.圓柱繞流293.橢圓繞流334.方腔流345.本章小結(jié)35六 總結(jié)與展望3662 一 緒論1.引言近半個(gè)多世紀(jì)以來(lái)，隨著科學(xué)技術(shù)水平的不斷提高，人們需要處理的科學(xué)及工程問(wèn)題也越來(lái)越復(fù)雜。面對(duì)這些問(wèn)題，傳統(tǒng)的實(shí)驗(yàn)和理論方法越來(lái)越難以滿足實(shí)際需要，時(shí)代呼吁新的科學(xué)手段的出現(xiàn)。 20世紀(jì)50年代計(jì)算機(jī)的誕生使這種設(shè)想成為可能。隨著計(jì)算機(jī)和計(jì)算方法的飛速發(fā)展,一種以計(jì)算方法為基礎(chǔ)、運(yùn)用計(jì)算機(jī)對(duì)科學(xué)及工程實(shí)際問(wèn)題進(jìn)行數(shù)值模擬的科學(xué)手段得以廣泛地應(yīng)用。 科學(xué)計(jì)算或曰計(jì)算機(jī)模擬也逐漸成為與傳統(tǒng)的科學(xué)方法理論和實(shí)驗(yàn)相并列的“ 第三種科學(xué)方法”，其發(fā)展水平也越來(lái)越成為某一國(guó)家在現(xiàn)代科學(xué)和高技術(shù)研究與開(kāi)發(fā)中競(jìng)爭(zhēng)能力的重要標(biāo)志。一般來(lái)說(shuō),這種發(fā)展水平的高低主要取決于兩個(gè)主要因素，即高性能計(jì)算工具( 計(jì)算機(jī)) 和高效計(jì)算方法。 前者是科學(xué)計(jì)算的硬件，后者是科學(xué)計(jì)算的軟件，兩者具有同等重要的地位。 因此,科學(xué)計(jì)算的發(fā)展也基本按照這兩條主線不斷發(fā)展。 對(duì)于高性能計(jì)算工具,處理器扮演最為關(guān)鍵的角色。 在過(guò)去的10年中, CPU 設(shè)計(jì)者受功率的限制不再在單一CPU 上添加更多的晶體管，而是把主要精力放在發(fā)展多核CPU 體系結(jié)構(gòu)上。 另一方面，由于圖形處理本身就是一個(gè)并行任務(wù)，采用流處理體系結(jié)構(gòu)的GPU更容易采用多核策略，且更適合于并行計(jì)算。因此，GPU 的通用計(jì)算能力越來(lái)越受到關(guān)注，基于GPU的并行計(jì)算硬件設(shè)計(jì)逐漸成為高性能計(jì)算領(lǐng)域的新寵，代表了高性能計(jì)算硬件發(fā)展的最新趨勢(shì)。 NVIDIA 公司無(wú)疑是這一領(lǐng)域的領(lǐng)導(dǎo)者，其2007年以來(lái)推出的GPU并行計(jì)算開(kāi)發(fā)環(huán)境CUDA(compute unified device architecture)則將這一領(lǐng)域的水平推向了一個(gè)前所未有的高度。 與此同時(shí)，作為高效計(jì)算方法的重要代表，格子Boltzmann 方法(lattice Boltzmann method, LBM)自產(chǎn)生到現(xiàn)在20余年間受到了計(jì)算流體力學(xué)等諸多領(lǐng)域內(nèi)學(xué)者的普遍青睞。 與傳統(tǒng)的建立在連續(xù)介質(zhì)模型上的計(jì)算流體力學(xué)(computational fluid dynamics, CFD)相比，該方法是20世紀(jì)80年代中期發(fā)展起來(lái)的一種流體建模與模擬的新方法， 它不僅具有算法簡(jiǎn)單、易于處理復(fù)雜邊界條件、壓力能夠直接求解、編程容易等優(yōu)勢(shì)，而且還具有天然的并行性,，這些優(yōu)勢(shì)使得LBM非常適合進(jìn)行大規(guī)模的并行計(jì)算。 因此，自GPU開(kāi)始用于通用計(jì)算的那一刻，國(guó)內(nèi)外就不斷有研究者試圖將GPU與LB 方法結(jié)合起來(lái)以便能夠高效處理大規(guī)模流體計(jì)算問(wèn)題。2.格子Boltzmann方法概述我們可以從不同的層次對(duì)流體的運(yùn)動(dòng)進(jìn)行描述。流體是由分子組成的，從微觀入手，可以在分子量級(jí)上以Hamilton方程對(duì)微觀粒子的運(yùn)動(dòng)進(jìn)行描述，可是這種方法很難實(shí)際應(yīng)用，例如，在溫度為0C和壓強(qiáng)為一個(gè)標(biāo)準(zhǔn)大氣壓的情況下，lcm3的空氣含有26910侈個(gè)分子，要求解每一個(gè)粒子的狀態(tài)是不太可能的。從宏觀入手，可以流速、壓力等物理量為基本變量，建立描述流體運(yùn)動(dòng)的宏觀方程，如Navier-Stokes方程，然后對(duì)該方程進(jìn)行求解從而解決具體的問(wèn)題。此外，還可以對(duì)流體在介于宏觀和微觀的一個(gè)尺度上進(jìn)行描述，稱之為介觀尺度(mesoscale)。以粒子分布函數(shù)為基本變量的Boltzmann方程就是在介觀尺度上對(duì)流體的一種描述，但是直接求解完整的Boltzmann方程是相當(dāng)困難的，為此，人們提出了離散速度模型，建立了格子Boltzmann方法。格子Boltzmann方法(Lattice Boltzmann Method，簡(jiǎn)稱LBM)源自LGA(Lattice Gas Automata)7,在1988年首次作為一種獨(dú)立的數(shù)值方法由McNamara 和Zanetti提出15。他們把格子氣自動(dòng)機(jī)中的整數(shù)運(yùn)算變成實(shí)數(shù)運(yùn)算，克服了格子氣自動(dòng)機(jī)的數(shù)值噪聲的缺點(diǎn)，但是，在他們的模型中，仍然采用較為復(fù)雜的LGA碰撞項(xiàng)。后來(lái)Bhatnagar-GrossKrook(BGK)算子被引入LBM，形成了Lattice BGK(LBGK)模型。該模型采用單一時(shí)間松弛方法，碰撞項(xiàng)十分簡(jiǎn)單，這使得LBGK模型開(kāi)始得到廣泛應(yīng)用。同時(shí)，該模型滿足了各項(xiàng)同性，GalIean不變性,并有獨(dú)立于速度的壓力項(xiàng)，在理論分析和數(shù)值模擬方面都具有很大靈活性。格子-Boltzmann模型保留了格子氣自動(dòng)機(jī)的優(yōu)點(diǎn)，克服了格子氣自動(dòng)機(jī)的不足，程序編制簡(jiǎn)單，計(jì)算效率較高。作為一種數(shù)值方法，LBM一般被認(rèn)為有下列優(yōu)點(diǎn)4：(1) 在LBM的控制方程中對(duì)流項(xiàng)為線性的，無(wú)需進(jìn)行特殊處理，算法也比較簡(jiǎn)單，編程相對(duì)容易； (2) 無(wú)需求解壓力方程便可通過(guò)粒子分布函數(shù)得到壓力；(3) LBM適于處理具有復(fù)雜幾何邊界的問(wèn)題，易于引入微觀效應(yīng)；(4) LBM具有天然的并行性，易于實(shí)現(xiàn)并行化。（1）格子Boltzmann方法的思想格子氣自動(dòng)機(jī)的基本思想是，把計(jì)算區(qū)域分成許多均勻的正六邊形（或正方形）的網(wǎng)格，而那些只有質(zhì)量無(wú)體積的粒子只能在網(wǎng)格點(diǎn)上存在，并沿著網(wǎng)格線在網(wǎng)格間運(yùn)動(dòng)。當(dāng)某一個(gè)粒子從某一網(wǎng)格點(diǎn)到鄰近的網(wǎng)格點(diǎn)時(shí)，有可能和從其他網(wǎng)格點(diǎn)到達(dá)該點(diǎn)的粒子相碰撞。然后根據(jù)碰撞規(guī)則再散射出去，演化為新的運(yùn)動(dòng)粒子流向各節(jié)點(diǎn)的鄰居，形成格子氣自動(dòng)機(jī)。源于LGA，格子Boltzmann方法是也一種應(yīng)用非連續(xù)介質(zhì)思想研究宏觀物理現(xiàn)象，求解流體力學(xué)問(wèn)題的方法。該法把流體及其存在的時(shí)間、空間完全離散，把流體看成由許多只有質(zhì)量沒(méi)有體積的微小粒子組成，所有這些粒子同步地隨著離散的時(shí)間步長(zhǎng)，根據(jù)給定碰撞規(guī)則在網(wǎng)格點(diǎn)上相互碰撞，并沿網(wǎng)格線在節(jié)點(diǎn)之間運(yùn)動(dòng)。碰撞規(guī)則遵循質(zhì)量、動(dòng)量和能量守恒定律。流體運(yùn)動(dòng)的宏觀特征則由微觀流體格子在整體上表現(xiàn)出來(lái)的統(tǒng)計(jì)規(guī)律表示。（2）格子Boltzmann方法與其他方法的比較4通過(guò)對(duì)二維剪切流的模擬對(duì)LBM和譜方法(Spectral Method)進(jìn)行比較，發(fā)現(xiàn)在計(jì)算結(jié)果的精度方面兩種方法基本一致，而在計(jì)算的時(shí)間效率方面LBM要優(yōu)于譜方法，并且在程序上LBM易于并行。以柱群繞流問(wèn)題的二維模擬為例對(duì)LBM和有限差分法(Finite Difference Method，簡(jiǎn)稱FDM)進(jìn)行比較，發(fā)現(xiàn)兩種方法的計(jì)算結(jié)果基本一致，但由于無(wú)需求解壓力方程，LBM的計(jì)算效率要優(yōu)于FDM，在對(duì)多孔介質(zhì)中的流動(dòng)進(jìn)行模擬時(shí)也得到了類似的結(jié)論。通過(guò)對(duì)具有復(fù)雜邊界的SMRX混合反應(yīng)器中的流動(dòng)的模擬發(fā)現(xiàn)，在保證同樣的精度的情況下，LBM對(duì)內(nèi)存的要求小于有限單元法(Fillite Element Method，簡(jiǎn)稱FEM)，并且LBM在模擬多相流和程序的可并行性方面更有優(yōu)勢(shì)。利用LBM和有限體積法(Finite Volume Method，簡(jiǎn)稱FVM)對(duì)管道中收縮段的層流流動(dòng)進(jìn)行模擬，發(fā)現(xiàn)以兩種方法得到的流速和應(yīng)力是一致的，但LBM相對(duì)FVM具有算法簡(jiǎn)單、計(jì)算效率高的優(yōu)勢(shì)。以圓柱繞流問(wèn)題的三維模擬為例對(duì)LBM和FVM進(jìn)行了比較，發(fā)現(xiàn)以兩種方法得到的結(jié)果總體上是一致的，但FVM得到的拖曳力系數(shù)稍高于試驗(yàn)值。（3）曲線邊界的處理經(jīng)典的LBM是在均勻的正方形網(wǎng)格上進(jìn)行計(jì)算的，當(dāng)固體邊壁的輪廓線不與網(wǎng)格線重合時(shí)，需要在邊界處進(jìn)行特殊處理以滿足固壁的流速條件，特別是對(duì)于具有曲線形式的固壁邊界。可以將物理邊界定義在節(jié)點(diǎn)的連線上，然后應(yīng)用反彈邊界條件。還可以將物理邊界定義在節(jié)點(diǎn)上，對(duì)于固體節(jié)點(diǎn)上的未知粒子分布函數(shù)，利用流體節(jié)點(diǎn)和邊界節(jié)點(diǎn)以外插的方法得到。采用這兩種方法能滿足一定的精度要求，但曲線邊界將成鋸齒狀。為了在曲線邊界實(shí)際位置處實(shí)現(xiàn)流速邊界條件，對(duì)反彈邊界條件進(jìn)行改進(jìn)，提出了獲得固體節(jié)點(diǎn)上的未知粒子分布函數(shù)的插值方法。該方法通過(guò)將非平衡態(tài)粒子分布函數(shù)外推的方法構(gòu)造固體節(jié)點(diǎn)上的未知粒子分布函數(shù)，并固壁邊界處將粒子分布函數(shù)進(jìn)行線性插值或二次插值，然后應(yīng)用反彈邊界條件。前者具有空間一階精度，后者具有空間二階精度，但是插值會(huì)引起質(zhì)量上的不守恒，為此，提出了一種基于有限體積概念的處理方法，但仍然存在質(zhì)量上的不守恒問(wèn)題。（4）非均勻網(wǎng)格下的LBM對(duì)于流場(chǎng)中速度或壓力變化較大的部分需要更密的網(wǎng)格以提高精度，若采用均勻的正方形網(wǎng)格，整體的計(jì)算量將增加很多，為克服均勻網(wǎng)格的缺點(diǎn)，近年來(lái)發(fā)展非均勻網(wǎng)格下的LBM成，提出的方法主要有三種，即:基于插值的LBM、局部網(wǎng)格加密法和將傳統(tǒng)數(shù)值方法與LBM結(jié)合的方法。通過(guò)插值的方法得到節(jié)點(diǎn)上的粒子分布函數(shù)，可以將LBM在曲線網(wǎng)格上進(jìn)行應(yīng)用。但由于在每一計(jì)算步都需進(jìn)行插值，該方法相對(duì)經(jīng)典的LBM計(jì)算量增大很多，同時(shí)，插值也會(huì)帶來(lái)一些數(shù)值上的粘性。利用泰勒級(jí)數(shù)展開(kāi)得到網(wǎng)格點(diǎn)上的粒子分布函數(shù)，并以最小二乘法進(jìn)行優(yōu)化，這一方法無(wú)需在每一步都進(jìn)行插值，計(jì)算量與經(jīng)典的LBM相當(dāng)。該方法具有無(wú)網(wǎng)格特性，但其更易于在結(jié)構(gòu)化的網(wǎng)格上進(jìn)行應(yīng)用。局部網(wǎng)格加密法是指在采用正方形網(wǎng)格的基礎(chǔ)上，保持網(wǎng)格的結(jié)構(gòu)不變，而在局部對(duì)網(wǎng)格進(jìn)行加密。一般先整個(gè)計(jì)算區(qū)域建立一個(gè)較粗的網(wǎng)格，再在某些局部的地方建立較細(xì)的網(wǎng)格，從而可以采用嵌套計(jì)算的模式，由于粗細(xì)網(wǎng)格的空間和時(shí)間尺度不一致，需要在粗細(xì)網(wǎng)格的交界處對(duì)粒子分布函數(shù)進(jìn)行轉(zhuǎn)換。該方法一般也需要進(jìn)行插值以獲得所需的粒子分布函數(shù)。將其它數(shù)值方法與LBM結(jié)合，即將FDM、FVM、FEM等方法與LBM相結(jié)合，分別形成FDLBM、FVLBM、FELBM。例如通過(guò)坐標(biāo)變換將FDLBM應(yīng)用到極坐標(biāo)，從而能更好地描述曲線邊界。將譜間斷有限元與LBM結(jié)合起來(lái)，精度高且容易實(shí)現(xiàn)程序上的并行。鑒于最小二乘有限元對(duì)數(shù)值振蕩有所改進(jìn)，將其引入LBM，可以加快收斂，而且該方法也容易實(shí)現(xiàn)程序上的并行。值得注意的是，對(duì)某些問(wèn)題以LBM進(jìn)行模擬，采用均勻網(wǎng)格仍是比較方便的，例如在模擬粒子懸浮流時(shí)，由于粒子在流場(chǎng)中的運(yùn)動(dòng)軌跡有時(shí)是很難預(yù)先知道的，所以不便于采用非均勻網(wǎng)格。（5）多松弛時(shí)間模式由于算法簡(jiǎn)單，采用單松弛時(shí)間模式的LBGK模型自提出后得到廣泛的應(yīng)用，為了完善LBM的數(shù)值穩(wěn)定性等問(wèn)題，近幾年，采用多松弛時(shí)間(MultipleRelaxation Time，簡(jiǎn)稱MIT)模式的LBM得到一定的發(fā)展和應(yīng)用。多松弛時(shí)間模式的碰撞過(guò)程不是在速度空間中進(jìn)行的，而是在矩空間中進(jìn)行的，將粒子分布函數(shù)投影到矩空間后，不同的矩分別對(duì)應(yīng)密度、動(dòng)量、能量、熱量流、應(yīng)力張量等不同的物理量，在碰撞過(guò)程中不同的矩可以采用不同的松弛時(shí)間，從而增加了模型的靈活性。多松弛時(shí)間模式的遷移過(guò)程仍是在速度空間中進(jìn)行的。多松弛時(shí)間模式對(duì)單松弛時(shí)間模式的穩(wěn)定性有一定的改進(jìn)，所能模擬的最大雷諾數(shù)一般比相同網(wǎng)格下的單松弛時(shí)間模式大四倍。在計(jì)算量方面，采用多松弛時(shí)間模式的計(jì)算時(shí)間一般要比單松弛時(shí)問(wèn)模式多1020。（6） LBM的應(yīng)用作為一種數(shù)值方法，LBM自提出后得到不斷的發(fā)展和完善，在很多方面都得到了應(yīng)用， 如多相流、粒子懸浮流、湍流、多孔介質(zhì)中的流動(dòng)、微尺度內(nèi)的流動(dòng)、生物體內(nèi)的流動(dòng)、粘彈性流體、室內(nèi)空氣的流動(dòng)等。3.GPU概述（1）GPU 的發(fā)展5 NVIDIA 公司1999 年推出了全球第一顆 GPU GeForce256。此后，NVIDIA 公司相繼推出第一款擁有可編程頂點(diǎn)著色器的GPU GeForce3以及第一款使用32位浮點(diǎn)流水線作為可編程的頂點(diǎn)處理器的GPU GeForce Fx，2006 年的GeForce 8 系列GPU率先采用了統(tǒng)一渲染架構(gòu)，以通用渲染單元替代了原來(lái)分離的著色單元和像素著色單元， 能夠支持DirecrtX 10.0。此后，CUDA 正式發(fā)布，引發(fā)了GPU通用計(jì)算的革命。GPU在處理能力和存儲(chǔ)帶寬上相對(duì)于CPU 有明顯優(yōu)勢(shì)，GPU 可以通過(guò)增加并行處理單元和存儲(chǔ)器控制單元的方式提高處理能力和存儲(chǔ)器帶寬。CPU 和GPU的浮點(diǎn)計(jì)算能力差異大的主要原因是：GPU 是專為密集運(yùn)算、高度并行計(jì)算而設(shè)計(jì)， 與CPU 相比，GPU內(nèi)部更多的晶體管用于數(shù)據(jù)處理而不是數(shù)據(jù)緩存或控制。 （2）CUDA 介紹 CUDA 是由NVIDIA 推出的通用并行計(jì)算架構(gòu)，可以使用GPU來(lái)解決商業(yè)、工業(yè)以及科學(xué)方面的復(fù)雜計(jì)算問(wèn)題。CUDA 采用類C 語(yǔ)言作為編程語(yǔ)言提供了大量的高性能計(jì)算指令開(kāi)發(fā)能力，使開(kāi)發(fā)者能夠在GPU的強(qiáng)大計(jì)算能力的基礎(chǔ)上建立起一種效率更高的密集數(shù)據(jù)計(jì)算解決方案，所編寫出的程序也就可以在支持CUDA 的處理器上以超高性能運(yùn)行。 CUDA 將GPU視為一個(gè)并行數(shù)據(jù)計(jì)算的設(shè)備，對(duì)所進(jìn)行的計(jì)算進(jìn)行分配和管理。 在CUDA 的架構(gòu)中，不需要像過(guò)去的GPGPU計(jì)算那樣將計(jì)算映射到圖形API(OpenGL和Direct 3D)中，從而降低了CUDA的開(kāi)發(fā)門檻。 4.本文的主要工作本文對(duì)一些經(jīng)典流場(chǎng)分別進(jìn)行LBM方法的CPU串行模擬和GPU并行模擬，并對(duì)各自性能進(jìn)行比較，以確定GPU并行計(jì)算的加速效果，并為確定最優(yōu)加速參數(shù)提供參考。為此，安排本文的主要工作及章節(jié)如下：第二章將對(duì)LBM的基本原理進(jìn)行介紹。第三章將對(duì)CUDA技術(shù)進(jìn)行介紹。 第四章將介紹LBM的具體實(shí)現(xiàn)。 第五章將對(duì)一些經(jīng)典流場(chǎng)的模擬結(jié)果進(jìn)行分析比較。第六章則對(duì)全文進(jìn)行總結(jié)。二 格子Boltzmann方法的基本原理LBM作為一種數(shù)值方法，它不同于傳統(tǒng)的有限差分法、有限元法及有限體積法，本章將主要就LBM的基本原理進(jìn)行介紹。1. Boltzmann方程8統(tǒng)計(jì)物理學(xué)中的Boltzmarm方程是在一定的假設(shè)條件下得到的，這些假設(shè)條件包括：(1)只考慮兩體碰撞；(2)兩個(gè)粒子在碰撞前其速度不相關(guān)，即分子混沌假定；(3)碰撞過(guò)程不受外力影響。在這些假設(shè)條件下，可得到一微分積分形式的Boltzmann方程，它描述的是粒子分布函數(shù)的演化過(guò)程，其形式為： （21）其中，為空間位置矢量，為粒子速度，代表粒子質(zhì)量，為一常數(shù)，為系統(tǒng)所受到的外部體積力，為單粒子分布函數(shù)，以下簡(jiǎn)稱粒子分布函數(shù)，的物理意義是在時(shí)刻粒子出現(xiàn)在空間域和速度域內(nèi)的概率，為碰撞項(xiàng)，該項(xiàng)為一積分項(xiàng)，它描述的是由于碰撞而引起的粒子分布函數(shù)的變化，設(shè)兩個(gè)粒子在碰撞前具有速度，碰撞后速度變?yōu)椋瑒t碰撞項(xiàng)的形式為： （22） 其中，為微分碰撞截面。在以后的討論中，無(wú)特別聲明，將忽略方程中的體積力項(xiàng)。概率密度與密度、速度、內(nèi)能等宏觀的物理量可分別采用如下統(tǒng)計(jì)定義： （2.3） 溫度可由內(nèi)能得到 （2.4）其中，為Boltzmalm常數(shù)。對(duì)N個(gè)粒子組成的系統(tǒng)在某一時(shí)刻的狀態(tài)進(jìn)行描述，需要6N個(gè)變量，包括3N個(gè)空間坐標(biāo)和3N個(gè)動(dòng)量，這6N個(gè)變量可以看成某6N維空間的坐標(biāo)，這個(gè)6N維空間在統(tǒng)計(jì)物理學(xué)中被稱為相空間。一個(gè)偏離平衡態(tài)的系統(tǒng)向平衡態(tài)逼近的過(guò)程稱為松弛(relaxation)。Boltzmann方程描述的是一個(gè)偏離平衡態(tài)不遠(yuǎn)的系統(tǒng)的松弛過(guò)程。經(jīng)過(guò)一定的時(shí)問(wèn)，系統(tǒng)將達(dá)到平衡狀態(tài)，此時(shí)進(jìn)入相空間某微元的粒子數(shù)和離開(kāi)此微元的粒子數(shù)將相等。從Boltzmmm方程得到解析解是非常困難的，一個(gè)主要的難點(diǎn)就是碰撞積分項(xiàng)的處理，物理學(xué)家們就此提出了一些簡(jiǎn)化的方法，即碰撞模型(collisionmodel)，其中為人們所熟知的一個(gè)模型是由Bhamagar、Gross和Krook于1954年共同提出的BGK模型，該模型也曾被Welander于同年獨(dú)立地提出，其形式為： （2.5）其中，為碰撞頻率，為碰撞時(shí)間，它是兩次粒子碰撞之間的平均時(shí)間，為Maxwell分布，也稱MaxwellBoltzmann分布，它是系統(tǒng)達(dá)到平衡時(shí)的粒子分布函數(shù)，其表達(dá)形式為： （2.5）其中，D為空間維數(shù)。從BGK模型的表達(dá)式可以看出，碰撞過(guò)程使粒子分布函數(shù)向平衡分布松弛，也被稱為松弛時(shí)間。2.格子Boltzmann方程（1）有限差分格式推導(dǎo)歷史沿革上，LBM是由LGA發(fā)展而來(lái)，但是格子Boltzmann方程還可以從Boltzmann方程得到，下面簡(jiǎn)單介紹文獻(xiàn)【16】給出的方法。取，則可將Boltzmann方程寫成如下形式 （2.6） 其中，為Maxwell分布： （2.7）在后面的討論中，將采用式(213)給出的Boltzmann方程的形式，密度、速度、內(nèi)能的定義將分別采用如下形式 （2.8）將連續(xù)的粒子速度空間離散，使粒子速度為有限個(gè)值，對(duì)應(yīng)的粒子分布函數(shù)為，簡(jiǎn)記為，滿足下面的演化方程 （2.9）假設(shè)密度和速度的計(jì)算采用如下所示的高斯型數(shù)值積分形式 （2.10）其中，為數(shù)值積分權(quán)重系數(shù)，設(shè)則上式可寫為： （2.11）滿足的演化方程為： （2.12）其中是平衡分布函數(shù)，將其變換形式得： （2.13） 當(dāng)流速較低時(shí)，上式乘積第三部分將很小，應(yīng)用泰勒公式展開(kāi)成二階為： （2.14）所以，平衡分布函數(shù)可寫為： （2.15）其中為平衡分布函數(shù)的權(quán)重系數(shù) （2.16）為聲速（2.17）設(shè)，將演化方程離散，可得到有限差分格式 （2.18）其中，為無(wú)量綱松弛時(shí)間，上式即是標(biāo)準(zhǔn)的LBGK方程。以上介紹了從Boltzmann方程出發(fā)如何得到LBGK方程，其中的權(quán)重系數(shù)取決于所采用的LBGK模型，下面本文將介值的推導(dǎo)過(guò)程。（2）權(quán)重系數(shù)的推導(dǎo)在LBKG模型中應(yīng)用較多的是Qian等人提出的DnQb系列模型，D為空間維數(shù)，Q為離散的粒子速度個(gè)數(shù)，較常用的有D2Q7，D2Q9，D3Q15，D3Q19，D3Q27等，其中二維模型中應(yīng)用較多的是D2Q9模型，下面以D2Q9模型為例，介紹權(quán)重系數(shù)推導(dǎo)過(guò)程。D2Q9模型是一種多速度(multispeed)模型，有三種粒子速度，如式(226)所示，模型的格子形狀見(jiàn)圖21。 圖21D2Q9模型的格子 （2.19） 其中，為格子速度。由對(duì)稱性可知，同一類的粒子速度對(duì)應(yīng)的權(quán)重系數(shù)應(yīng)相同，不妨將平衡分布寫成下面的形式 （2.20）其中 （2.21）平衡分布函數(shù)權(quán)重系數(shù)的選取應(yīng)使其滿足質(zhì)量守恒、動(dòng)量守恒、各向同性等約束條件，即 （2.22）其中，為壓力，為Kronecker記號(hào)， （2.23）將式(2.20)代入式(2.22)，得 (2.24) (2.25) （2.26） 考慮到密度應(yīng)和速度不相關(guān)，由（2.24）可得 (2.27) 由（2.25）化簡(jiǎn)得 （2.28）由（2.26）對(duì)應(yīng)項(xiàng)相等得 （2.29） 考慮到壓力應(yīng)和速度不相關(guān)，可得 （2.30）綜上，有 （2.31）解得： （2.32）（3）松弛時(shí)間的推導(dǎo)Navier-Stokes方程是從宏觀角度對(duì)流體的行為進(jìn)行描述的，Boltzmann方程是從微觀角度對(duì)流體的一種描述，其宏觀表現(xiàn)形式與Navier-Stokes方程是一致的，借助ChapmanEnskog展開(kāi)法，可由Boltzmann方程推導(dǎo)出Navier-Stockes方程。ChapmanEnskog展開(kāi)法是由Chapman和Enskog于1910年到1920年間發(fā)展起來(lái)的一種多尺度展開(kāi)法，其基本思想就是將物理量在多級(jí)尺度上進(jìn)行展開(kāi)，同級(jí)尺度上物理量的量級(jí)保持一致。同時(shí)可以得到： （2.33） 其中，為運(yùn)動(dòng)粘度系數(shù)。則解出： （2.34）3.格子Boltzmann方法的實(shí)現(xiàn)上一節(jié)介紹了格子Boltzmann方程及其宏觀表現(xiàn)，格子Boltzmann方程是LBM的演化方程，根據(jù)該方程，可以編寫程序在計(jì)算機(jī)上就具體的問(wèn)題進(jìn)行數(shù)值模擬，本節(jié)將以D2Q9模型為例，就在計(jì)算機(jī)上如何實(shí)現(xiàn)LBM進(jìn)行基本的介紹。在實(shí)際的計(jì)算中人們常常以無(wú)量綱化的方式進(jìn)行計(jì)算，分別格子尺寸、碰撞時(shí)間為無(wú)量綱化參數(shù)單位，即： （2.35）可以得到其它物理量的無(wú)量綱量。 （2.37） （2.38） （2.39）（2.40） （2.41）其它物理量的無(wú)量化過(guò)程類同，這里不再贅述。經(jīng)過(guò)無(wú)量綱化后，LBM模型中物理量的單位為lu(1atticeunit)。當(dāng)我們通過(guò)LBM計(jì)算得到某物理量后，需根據(jù)上面的過(guò)程將其還原成真實(shí)的物理量。另外，演化方程式(225)其實(shí)包含了兩個(gè)基本的過(guò)程，將其寫成如下形式 （2.42）式(291)描述的是碰撞過(guò)程，式(292)描述的是遷移過(guò)程， 為碰撞后的粒子分布函數(shù)，可稱為碰后粒子分布函數(shù)，將演化方程拆分為碰撞和遷移兩個(gè)過(guò)程便于計(jì)算機(jī)程序的實(shí)現(xiàn)。宏觀量的計(jì)算可用以下公式： （2.11）綜上所述，在算法上可以采用下面給出的一個(gè)計(jì)算流程：(1)給定初始的密度和流速，設(shè)定粒子分布函數(shù)的初始值為 。(2)以式(291)進(jìn)行碰撞過(guò)程和遷移過(guò)程。(3)施加邊界條件。(4)以式(220)計(jì)算新時(shí)刻的密度和流速。(5)以式(227)計(jì)算新時(shí)刻平衡分布函數(shù)。(6)重復(fù)第(2)（5）步。通過(guò)上面的討論，可以發(fā)現(xiàn)粒子分布函數(shù)，是計(jì)算中的主要變量，其它的物理量都是以此為基礎(chǔ)得到的，邊界條件的給定也必須最終以粒子分布函數(shù)來(lái)表示，邊界條件的具體施加方法將在下一節(jié)及以后的章節(jié)中進(jìn)行討論。同時(shí)，我們還可以發(fā)現(xiàn)另外一個(gè)特點(diǎn)。對(duì)某一個(gè)節(jié)點(diǎn)來(lái)說(shuō)，其碰撞過(guò)程不需要其它節(jié)點(diǎn)的信息，而遷移過(guò)程只需要和它相鄰的節(jié)點(diǎn)的信息，這使得計(jì)算具有很好的局部性，由此可見(jiàn)LBM具有天然的并行性。4.基本的邊界條件13對(duì)于位于邊界上或邊界附近的節(jié)點(diǎn)，當(dāng)遷移過(guò)程發(fā)生后，一些粒子分布函數(shù)是未知的，此即通過(guò)邊界條件要解決的問(wèn)題之一。由于在LBM的計(jì)算中以粒子分布函數(shù)為基本的變量，所以邊界條件的給出必須以粒子分布函數(shù)來(lái)描述，所有關(guān)于宏觀變量的邊界條件都需以粒子分布函數(shù)來(lái)表示，這不同于傳統(tǒng)的數(shù)值方法。本節(jié)將以D2Q9模型為例介紹幾個(gè)基本的邊界條件，其它邊界條件將在以后的章節(jié)中討論。（1）周期性邊界條件以左右邊界為例，在施加周期性邊界條件時(shí)，使右邊界節(jié)點(diǎn)上的粒子分布函數(shù)Z、石、像內(nèi)部節(jié)點(diǎn)那樣遷移到左邊界相應(yīng)節(jié)點(diǎn)上，如圖22所示，同樣，使左邊界上的相應(yīng)粒子分布函數(shù)遷移到右邊界。圖22周期性邊界條件示意圖（2）反彈邊界條件反彈邊界條件常被用來(lái)實(shí)現(xiàn)無(wú)滑移邊界，粒子分布函數(shù)在遷移中若碰到固邊界便沿與原方向相反的方向返回，據(jù)物理邊界定義的位置，反彈邊界條件又可分為全程反彈(fullwaybounceback)邊界條件和半程反彈(halfwaybounceback)邊界條件，其程序?qū)崿F(xiàn)過(guò)程如圖23所示。圖23半程反彈(a)和全程反彈(b)程序?qū)崿F(xiàn)過(guò)程示意圖。圖中省略了其它粒子分布函數(shù)灰色的方塊表示的是物理邊界的位置。半程反彈邊界條件的物理邊界定義在邊界節(jié)點(diǎn)和與其相鄰的流體節(jié)點(diǎn)連線的中點(diǎn)，而全程反彈邊界條件的物理邊界定義在邊界節(jié)點(diǎn)上，兩種邊界條件都要滿足一點(diǎn)，即粒子分布函數(shù)在內(nèi)實(shí)際遷移的距離在數(shù)值上等于到相應(yīng)鄰節(jié)點(diǎn)的距離。對(duì)于半程反彈邊界條件，可以虛設(shè)一層邊界節(jié)點(diǎn)以在程序上實(shí)現(xiàn)該邊界條件，文獻(xiàn)3提出了一種不需要虛設(shè)邊界節(jié)點(diǎn)的實(shí)現(xiàn)方法。在邊界節(jié)點(diǎn)上，兩種反彈邊界條件都不發(fā)生碰撞過(guò)程，同時(shí)，若在計(jì)算中引入體積力，在邊界節(jié)點(diǎn)上也都不考慮體積力。（3）反射邊界條件反射(specularreflection)邊界條件常被用來(lái)實(shí)現(xiàn)自由滑移邊界條件，即在邊界上只要求法向速度為零。粒子分布函數(shù)在遷移中若碰到固邊界便如光的反射現(xiàn)象那樣沿反射光方向繼續(xù)遷移，如圖24所示。圖24反射邊界條件示意圖（4）非平衡外推邊界條件10 在邊界處，未知的分布函數(shù)可以用邊界處的平衡分布函數(shù)，和其相鄰內(nèi)側(cè)格子的分布函數(shù)表示。算法如下： （2.43）其中 為邊界格子的分布函數(shù)， 為邊界相鄰內(nèi)側(cè)格子的分布函數(shù)。其中 為邊界格子的平衡分布函數(shù)， 為邊界相鄰內(nèi)側(cè)格子的平衡分布函數(shù)。5.本章小結(jié)本章從幾方面對(duì)LBM的基本原理進(jìn)行了較詳細(xì)的介紹。首先，對(duì)格子Boltzmann方程的推導(dǎo)及LBM的算法實(shí)現(xiàn)進(jìn)行了詳細(xì)的介紹。然后，對(duì)在LBM中應(yīng)用的基本的邊界條件進(jìn)行了介紹。三 CUDA技術(shù)隨著顯卡的發(fā)展，GPU越來(lái)越強(qiáng)大，而且GPU為顯示圖像做了優(yōu)化。在計(jì)算上已經(jīng)超越了通用的CPU。如此強(qiáng)大的芯片如果只是作為顯卡就太浪費(fèi)了，因此NVidia推出CUDA，讓顯卡可以用于圖像計(jì)算以外的目的。 1. CUDA 編程模型CUDA1,14編程模型用CUDA C語(yǔ)言實(shí)現(xiàn)，它是C / C + +的一種擴(kuò)展。在CUDA的C程序中函數(shù)屬于以下三個(gè)類別之一：1主機(jī)代碼,即函數(shù)由CPU運(yùn)行。2內(nèi)核,即主機(jī)代碼調(diào)用并由CPU運(yùn)行。3設(shè)備函數(shù),即函數(shù)由GPU運(yùn)行，而且由內(nèi)核函數(shù)或其他設(shè)備函數(shù)調(diào)用。內(nèi)核并行運(yùn)行在GPU上，通過(guò)在運(yùn)行時(shí)指定一個(gè)網(wǎng)格給出執(zhí)行模式。線程被分成相同的數(shù)組，叫做塊（Block），進(jìn)而組裝形成執(zhí)行網(wǎng)格（Grid），如圖3.1。塊可以有三個(gè)尺寸,網(wǎng)格只能是一或二維。每個(gè)線程用預(yù)定義的blockIdx和threadIdx結(jié)構(gòu)的x,y,和z值標(biāo)示。 圖3.1一個(gè)內(nèi)核的局部變量是只屬于該線程的，聲明為共享時(shí)，在一個(gè)塊的所有線程都可以訪問(wèn)。共享變量沒(méi)有互斥機(jī)制，程序員使用塊同步語(yǔ)言來(lái)管理。內(nèi)核也可以訪問(wèn)全局內(nèi)存，它對(duì)執(zhí)行網(wǎng)格的每個(gè)線程都是可見(jiàn)的。這種訪問(wèn)也沒(méi)有保護(hù)機(jī)制。全局內(nèi)存在內(nèi)核執(zhí)行期間一直存在。全局內(nèi)存可被主機(jī)代碼訪問(wèn),它通常也是CPU和GPU之間通信路徑。最后，值得一提的是,線程也可以訪問(wèn)常數(shù)內(nèi)存和紋理內(nèi)存。常數(shù)存儲(chǔ)器是存儲(chǔ)運(yùn)行時(shí)不變參數(shù)的一種方便的方法。2.CUDA硬件一個(gè)能夠支持CUDA 的GPU內(nèi)包含一組流多處理器(SMs)。每個(gè)SM包含幾個(gè)標(biāo)量處理器(SPs)、寄存器、共享內(nèi)存和常量和紋理緩存,見(jiàn)圖3.2。此外, 芯片外GPU還有設(shè)備內(nèi)存。 圖3.2對(duì)于每個(gè)SM,寄存器和共享內(nèi)存都很快,但容量很小。相比之下，設(shè)備內(nèi)存, 容量大，但存取速度慢。表1總結(jié)了計(jì)算本文使用的GeForce GT540M處理器技術(shù)規(guī)格。表3.1.GeForce GT540M處理器技術(shù)規(guī)格SM個(gè)數(shù)96每個(gè)SM中SP數(shù)8每個(gè)SM的寄存器大小32K每個(gè)SM的共享存儲(chǔ)大小48K每個(gè)SM的常數(shù)存儲(chǔ)大小64K設(shè)備內(nèi)存大小1G一個(gè)內(nèi)核的局部變量存儲(chǔ)在寄存器, 數(shù)組除外，因?yàn)榧拇嫫鞑豢蓪ぶ贰１镜財(cái)?shù)組存儲(chǔ)在所謂的本地內(nèi)存,這事實(shí)上是托管在設(shè)備內(nèi)存。3.優(yōu)化準(zhǔn)則執(zhí)行網(wǎng)格的一個(gè)線程塊只能被單個(gè)SM處理。然而，一個(gè)SM可以同時(shí)處理幾個(gè)塊。這就導(dǎo)致了一些有關(guān)網(wǎng)格布局的限制，列在表2。為了利用GPU大規(guī)模并行架構(gòu)，當(dāng)前活動(dòng)線程的數(shù)量應(yīng)盡可能大。它是由程序員定義一個(gè)網(wǎng)格實(shí)現(xiàn)這一目標(biāo)，同時(shí)要避免寄存器溢流。然而，占有率，即活動(dòng)線程與最大可能數(shù)量的比率，通常不是一個(gè)可靠的性能指標(biāo)。數(shù)據(jù)密集型應(yīng)用程序，限制因素更可能是全局內(nèi)存的最大吞吐量。盡管如此, 需要規(guī)定一個(gè)最小的占有率以隱藏全局內(nèi)存延遲。表3.2.計(jì)算能力2.1的網(wǎng)格布局限制每個(gè)SM能同時(shí)處理的最大線程數(shù)1536每個(gè)線程塊最大線程數(shù)1024線程塊維度1024X1024X64網(wǎng)格維度當(dāng)在SM運(yùn)行時(shí),一塊線程被分割成32線程的線程束。為實(shí)現(xiàn)實(shí)際的并行性,在一個(gè)線程束的所有線程經(jīng)必須遵循相同的指令路徑。當(dāng)分支發(fā)散時(shí),執(zhí)行的分支路徑串行化,這會(huì)大大影響性能。所以,分支粒度最好是一個(gè)線程束大小的倍數(shù)。在一個(gè)半束的線程訪問(wèn)同一共享內(nèi)存單元的不同的位置導(dǎo)致單元沖突,這要用事務(wù)串行化的方法解決。然而, 只要不發(fā)生單元沖突，共享內(nèi)存能和寄存器一樣快。共享內(nèi)存的主要目的是用于塊內(nèi)通信。不過(guò),共享內(nèi)存也方便從全局內(nèi)存預(yù)讀取數(shù)據(jù)，存儲(chǔ)小的本地?cái)?shù)組,或避免寄存器溢流。使用共享內(nèi)存有助于減少與設(shè)備內(nèi)存交換數(shù)據(jù),減小數(shù)據(jù)交換延遲，可能對(duì)性能有重大改善。 4.本章小結(jié)在本節(jié)中,我們將首先對(duì)CUDA編程模型給出一個(gè)簡(jiǎn)短的回顧。然后,我們大概地描述了CUDA硬件和用于本文計(jì)算的GeForce GT540M GPU。最后,我們將討論為達(dá)到最佳效率而應(yīng)考慮得因素。四 算法實(shí)現(xiàn)本章用LBM的LBGK 方程按照D2Q9模型計(jì)算流場(chǎng)，完成了圓柱繞流的GPU 程序模擬，以清楚地說(shuō)明LBM的GPU實(shí)現(xiàn)。 1.算法方面在第二章已經(jīng)給出了算法的流程，本節(jié)以圓柱繞流程序?yàn)槔o予詳細(xì)說(shuō)明。 圖4.1 圓柱繞流示意圖如圖4.1所示，圓柱繞流的左邊為速度邊界，右邊自由流出，上下則是墻壁，中間是固體圓。因此決定在四邊采用非平衡外推法，中間使用反射法。為采用均勻的正方形格子劃分流場(chǎng)，流場(chǎng)長(zhǎng)寬格子數(shù)之比應(yīng)等于流場(chǎng)長(zhǎng)寬之比。有運(yùn)動(dòng)粘度計(jì)算出松弛時(shí)間。假設(shè)每個(gè)格子的密度，速度，算出該初始密度，速度下的平衡分布函數(shù)作為初始分布函數(shù)。然后所有節(jié)點(diǎn)用下面式子完成碰撞過(guò)程: （2.42）碰撞之后，非邊界格子由下面式子獲取新的分布函數(shù): （2.42）中間邊界格子由反彈邊界有: （4.1） 其中-表示的反方向。四周邊界格子用非平衡外推法。左邊，上邊，下邊速度已知，密度等于相鄰內(nèi)側(cè)點(diǎn)密度。右邊密度和速度皆等于相鄰內(nèi)側(cè)格子。然后有下面式子求出分布函數(shù)： （2.43）至此，求出全場(chǎng)的新一輪分布函數(shù)，以當(dāng)前值可以類似求下一輪分布函數(shù)，循環(huán)下去直到平衡某次循環(huán)。最后，用下面式子