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

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