基于CUDA的火焰模擬

滿家巨，鄒有，陳傳淼

(湖南師范大學數(shù)學與計算機科學學院高性能計算與隨機信息處理省部共建教育部重點實驗室, 中國 長沙 410081)

CUDA(General Unified Device Architecture)是NVIDIA公司在2007年正式發(fā)布的GPU通用計算開發(fā)環(huán)境和軟件體系，它不但為用戶提供了一系列操作GPU數(shù)據(jù)讀寫和命令流控制的API，讓用戶應用GPU進行通用計算(GPGPU)，還提供了OpenGL和Direct 3D進行互操作的接口，使計算之后的結果可以快速地繪制出來．隨著廉價圖形卡的處理能力日漸強大，許多行業(yè)已經(jīng)開始使用CUDA進行加速計算，如音視頻處理、石油勘探、醫(yī)學成像、科學計算、物理模擬等．

在計算機的計算能力有限時，研究者們只能通過參數(shù)建模的方法來進行流體模擬，如將波浪函數(shù)的線性波組合來模擬浪花[1]、使用拉格朗日粒子模擬波浪[2]等，但這些方法不能模擬比較復雜的、細節(jié)更為豐富的流體效果，于是研究者們開始把注意力轉(zhuǎn)向基于物理的流體模型．流體模擬使用的基于物理的方法主要有：晶格波爾茲曼方法(LBM)[3]、拉格朗日方法[4]以及基于N-S方程的歐拉方法[5]．其中歐拉方法是最常使用的一種計算方法，能夠非常容易地重構光滑的流體表面，借助文獻[6]提出的半拉格朗日法進行求解，可以選取非常大的時間步長并能保持計算的穩(wěn)定，但缺點在于計算量比較大，需要非常高效的方程組求解器才能實現(xiàn)實時模擬．

本文以歐拉方法進行火焰模擬為線索，對模型的計算復雜度進行分析，實現(xiàn)了基于CUDA的線性方程組的求解器；針對火焰模擬中的顏色、漩渦、運動形態(tài)等內(nèi)容，引入漩渦約束因子和隨機風力來制造火焰的搖曳效果，并提出一種簡便的煙霧和焰心效果繪制策略來增加模擬的真實性，最終計算結果通過基于CUDA的光線投射方法實現(xiàn)可視化．

1 CUDA編程模型與圖形API

相較于CPU而言，CUDA主要結構為計算邏輯單元而不是控制邏輯單元，因此其浮點運算能力可達到CPU的幾十倍甚至上百倍．CUDA將對GPU硬件的操作予以屏蔽，讓用戶可以使用簡單的C語言接口實現(xiàn)通用計算．

1.1 CUDA編程模式

一個CUDA程序由主機端程序和設備端程序組成.主機端程序是在CPU上執(zhí)行的串行代碼，而設備端程序是在GPU上執(zhí)行的并行代碼，即kernel程序．使用CUDA進行編程時，首先準備好主機端數(shù)據(jù)，然后將其從主機的內(nèi)存復制到設備內(nèi)存中，調(diào)用kernel程序進行高密度的并行計算.待計算完成再將計算結果從設備端復制回主機內(nèi)存．如圖1，在CUDA架構中，線程被分成許多彼此獨立線程塊(Block)，同一個線程塊內(nèi)的所有線程使用共享內(nèi)存協(xié)調(diào)計算任務，一個GPU上的所有線程塊組織成網(wǎng)格(Grid)．

圖1 線程塊與線程的組織形式[7] 圖2 CUDA存儲模型[7]

1.2 CUDA存儲模式

CUDA為滿足數(shù)據(jù)的不同需求提供了6種存儲器(圖2)：寄存器(Registers)、本地存儲器(Local Memory)、共享存儲器(Shared Memory)、全局存儲器(Global Memory)、常數(shù)存儲器(Constant Memory)和紋理存儲器(Texture Memory)．其各自特點如下：

(1)寄存器的訪問效率最高，每個線程都擁有專屬的寄存器，但它的數(shù)量比較少；

(2)本地存儲器訪問速度較慢，也是線程特有的存儲單元，主要在寄存器空間不足時使用；

(3)共享存儲器訪問速度快，是線程塊所有線程公有的存儲單元，其大小也比較有限.常見的硬件每個線程塊可使用的共享存儲器大小只有16 KB．

(4)全局存儲器的大小與顯卡的顯存大小相當，它可通過主機端程序和設備端程序訪問，充當主機端和設備端的數(shù)據(jù)交互媒介．

(5)常數(shù)存儲器是顯存中的只讀存儲單元，它擁有緩存加速的功能，可以節(jié)約帶寬而加快訪問速度，因此適合存儲程序中需要頻繁訪問的只讀參數(shù)．

(6)紋理存儲器也是只讀存儲器，它由GPU用于紋理渲染的圖形專用單元發(fā)展而來，具有地址映射、數(shù)據(jù)濾波、緩存等功能，適合實現(xiàn)圖像處理、查找表或大量數(shù)據(jù)隨機訪問等操作．

1.3 CUDA與圖形API互操作

CUDA為圖形開發(fā)用戶提供了與OpenGL和Direct 3D這2種主流的圖形API互操作接口．一些OpenGL和Direct 3D的資源可以被映射到CUDA的地址空間，使得計算數(shù)據(jù)可以通過OpenGL或Direct 3D來顯示．本文中使用CUDA與OpenGL的互操作API來實現(xiàn)可視化，將OpenGL緩沖對象映射到CUDA的地址空間，這樣就可以在CUDA中讀取OpenGL寫入的數(shù)據(jù)，也可以用CUDA寫入數(shù)據(jù)待OpenGL使用．

2 基本方程及求解

2.1 基本方程

火焰的燃燒模型通常被當成不可壓縮流來處理[8]，其速度場所對應的動量方程和連續(xù)性方程(Navier-Stoke方程)為(1)和(2)所示，式(3)為火焰的溫度場方程．

(1)

(2)

(3)

根據(jù)Helmholtz-Hodge分解定理對基本方程進行投影操作，可將方程分解為多個步驟分別求解[9]：

(4)

(5)

(6)

(7)

(8)

2.2 計算復雜度分析

根據(jù)半拉格朗日方法，式(4)～(8)的求解過程可分為4步：添加源、對流、擴散和投影．使用邊長為N個劃分的立方體網(wǎng)格進行計算，總的節(jié)點數(shù)為N3．

(1)添加源操作中，需要將速度場和溫度場進行初始化和用戶交互的輸入．給定一個速度場(三維)和溫度場，需對每個節(jié)點的速度和溫度進行賦值操作，共需4N3次乘-加運算．

(2)擴散步中，速度場三維數(shù)據(jù)和溫度場一維數(shù)據(jù)需全部更新，而每節(jié)點數(shù)據(jù)都由其周邊6個位置的數(shù)據(jù)插值所得，假設最大迭代次數(shù)為MAX_ITER，則擴散過程共需要MAX_ITER×4×6N3次浮點運算．

(3)投影由散度計算、線性方程求解和壓力修正3部分組成．首先對每一個網(wǎng)格單元進行速度場的散度計算，散度的計算只針對速度場進行，計算涉及每單元的周邊6個數(shù)值的乘-加操作；線性方程組求解的計算方法和復雜度與擴散步相同；壓力修正即對速度場的3個分量使用對應的2個壓力值進行調(diào)整，各需2次乘加操作．整個投影步需要進行的浮點操作數(shù)為(57+MAX_ITER×6)N3次．

(4)對流步計算也是速度場和溫度場共同參與的計算過程，經(jīng)插值和越界判斷，共需45N3次浮點運算．

綜上所述，在使用穩(wěn)定流方法進行三維的火焰模型計算中，總共需要進行的浮點運算為(106+30MAX_ITER)N3次．

2.3 基于CUDA的線性方程組求解器

由上一小節(jié)可知，模型的求解過程中涉及大量線性方程組的求解，因此加速線性方程組的求解是提高模擬速度的關鍵．分析方程可知，方程組系數(shù)矩陣均為大規(guī)模的稀疏矩陣，其非零元呈對角分布且主對角占優(yōu)，因而非常適合迭代法求解，并可采用一維數(shù)組來存儲以達到高效訪存．

根據(jù)CUDA編程模型，將三維數(shù)據(jù)按高度所在方向分層分配成線程塊，并使用循環(huán)進行遍歷，每一層的網(wǎng)格對應的值可以用線程的索引來獲取：

確定對數(shù)值的訪問模式后，就可以開始設計迭代法求解了．本文實現(xiàn)了3種比較常見的迭代法：Jacobi迭代、Gauss-Seidel和共軛梯度法(CG)，并分別進行求解實驗，表1～3分別展示了3種迭代法在不同問題規(guī)模下相對于純CPU運算的加速情況．本文實驗使用的配置為：CPU：Intel 四核2.34 Hz；內(nèi)存：4.0 GB；顯卡型號：Geforce GTX 285；顯存大小：2 GB；操作系統(tǒng)：Windows 7．

表1 Jacobi迭代(迭代20次)

表2 Gauss-Seidel迭代(迭代20次)

表3 共軛梯度法(迭代20次)

通過以上3個表的數(shù)據(jù)可以看出，在網(wǎng)格規(guī)模較小時，GPU的加速比相對較小，并行加速的優(yōu)勢不能完全發(fā)揮，這主要是因為此時線程數(shù)目較少，GPU的計算和存儲資源不能被充分利用．隨著網(wǎng)格數(shù)目的增大，加速比開始呈上升趨勢．在迭代同樣多的次數(shù)時，Jacobi迭代的計算速度更勝一籌．

3 燃燒效果控制

3.1 漩渦限制因子

基于N-S方程模擬火焰時，通常將其看作是低粘性的流體，其運動會產(chǎn)生一定的旋渦現(xiàn)象．然而，采用半拉格朗日方法進行計算時會引起一定程度的數(shù)值耗散[8]，使得漩渦效果減弱．為了使模擬效果更加逼真，本文使用文獻[10]提出的漩渦限制因子方法來恢復這種細節(jié)，通過捕捉速度場中每一個計算單元上的旋轉(zhuǎn)運動，計算出加強或維持這種旋轉(zhuǎn)所需要的粘滯力，以此來減少由數(shù)值耗散引起的細節(jié)丟失．圖3展示了加入不同大小人工漩渦得到的效果.可以看出，加入人工旋渦后，模擬的結果表現(xiàn)出更多的形態(tài)變化效果，從而更富有表現(xiàn)力．

(a)無人工漩渦 (b)少量人工漩渦 (c)較多人工漩渦圖3 人工漩渦效果

3.2 顏色控制與焰心繪制

真實情況下的火焰是一個非常復雜的系統(tǒng)，火焰的顏色也受到許多因素的影響，如燃料的種類，燃燒的程度，燃燒時的氣壓等，每一個條件的變化都可能引起火焰顏色的變化．本文從火焰顏色的物理模型出發(fā)，使用黑體輻射表來確定火焰的顏色，黑體輻射可以通過文獻[11]中的工具進行計算．

在通常的模擬中，火焰燃燒產(chǎn)生的煙霧場和火焰是作為2個獨立的模型來處理的，這種方法可以很方便地控制煙霧外觀和數(shù)量，然而計算量也大大增加．本文采用一種簡便方法，把煙霧的顏色放到火焰的顏色場之中處理，將煙霧看作是火焰的外延．其基本過程描述如下：

(1)給定溫度范圍[T0,T1]，計算出對應的顏色表；

(2)根據(jù)燃料的性質(zhì)定義一種煙霧的顏色，如木材燃燒一般冒黑煙，可以把煙霧定義成黑色．把定義的煙霧顏色加入顏色表尾部．

(3)對于計算出的溫度場中所有溫度小于T0的點，都使用煙霧的顏色進行繪制．

圖5展示了采用這一方案進行繪制的煙霧效果.

圖4 將煙霧看成是主焰的延伸

(a)無煙霧效果 (b)增加煙霧效果圖5 煙霧效果對比

對于焰心的繪制，本文摒棄了理論和計算都頗為復雜的燃燒模型，直接利用OpenGL進行圖形繪制時的Alpha通道，使處于火焰中心位置的人造焰心與火焰重疊出比較逼真的焰心效果．如圖6所示，對于一個給定的三維火焰效果，從其中心挖去一個平截頭體作為人工焰心，對挖掉的部分填充焰心顏色．在繪制啟用Alpha通道后，視點位置就可以透過外層火焰看到內(nèi)部的焰心，內(nèi)部的焰心顏色會和外層不斷變化的火焰顏色疊加在一起，形成逼真的燃心效果．更進一步，對挖掉的平截頭體的大小和形狀增加一定的隨機性，就可以模擬出隨機跳動的焰心效果了．圖7展示了使用本文方法加入了焰心的繪制效果．

圖6 焰心繪制示意圖

(a)無焰心效果 (b)加入焰心效果圖7 焰心繪制效果

4 模擬效果可視化

4.1 體繪制

體繪制是一個從三維數(shù)據(jù)集生成其在二維平面的投影過程，實現(xiàn)體繪制最常見的方法為光線投射(Ray Marching)．光線投射基礎步驟為：

圖8 繪制流程

(1)投射：視點出發(fā)向所有的像素發(fā)射一條射線；

(2)采樣：沿著射線方向，使用相同的間隔選取采樣點；

(3)插值：對每一個采樣點，使用其周圍的點(通常有8個點)進行顏色值和Alpha值的插值運算，求出該采樣點的顏色值和Alpha值；

(4)累加：把在射線上的所有采樣點的顏色值和Alpha值進行由后向前的合成，即可獲得每條光線對應的像素點處的顏色值和Alpha值．

顯然，要通過光線投射獲得最終的畫面，需要對窗口上的每一個像素點都進行一次計算和顏色累加，對于規(guī)模較大的模擬，其計算量相當大．

4.2 基于CUDA的光線投射算法

光線投射的計算量雖然龐大，但任意2條光線的投射過程沒有任何關聯(lián)，非常適合使用GPU來并行計算．使用CUDA實現(xiàn)光線投射大致包括綁定像素緩存、光線投射計算、解決緩存綁定、繪制結果4個步驟[12]，其中最重要的部分為光線投射計算．其基本實現(xiàn)思路如下：

(1)設定繪制窗口大小，盡量使用16的整數(shù)倍；

(2)將繪制窗口中的像素映射到16×16的線程塊中，每個線程塊處理256個像素點；

(3)計算每個像素點在世界坐標中的位置，計算從視點指向像素點的射線方向；

(4)求解每條射線與計算空間的6個面的交點；

(5)在最遠的交點和最近的交點之間以某固定距離采樣，計算采樣點對應的顏色值；

(6)將結果寫入像素緩存中用于最終的顯示．

4.3 可視化效果展示

到此為止，本文已經(jīng)完成了一整套火焰模擬的計算和繪制方案．通過圖9，可以明顯地看出原始模擬在依次加入人工漩渦、煙霧和焰心繪制后得到的不同效果．圖10展示了加入場景繪制的效果圖，其可視化效果基本可以滿足一般游戲或普通動畫要求．

(a)原始效果 (b)加入漩渦 (c)繪制煙霧 (d)繪制焰心圖9 本文實現(xiàn)的效果

圖10 燃燒的火把和蠟燭

5 結語

本文基于CUDA實現(xiàn)了火焰模擬的計算、形狀顏色控制以及可視化全過程，工作包括線性方程組求解器的實現(xiàn)、燃燒效果控制和可視化3部分，針對方程組求解和可視化，實現(xiàn)了基于CUDA的3種線性求解器和光線投射算法，針對燃燒效果的控制提出了簡化的煙霧和焰心繪制方法，總體達到了比較可觀的模擬結果．由于本文的立足點并不是模型改進或算法改良方面的工作，因此后續(xù)的工作還有很多，比如網(wǎng)格改進、模擬對象改變等都是比較好的研究方向．

參考文獻：

[1] DARWYN R, PEACHEY. Modeling waves and surf [J]. Computer Graphics, 1986,20(4): 65-74.

[2] ALAIN F, WILLIAM T. A simple model of ocean waves [J]. Computer Graphics,1986,20(4):75-84.

[3] THüREY N, RüDE U. Free surface lattice-Boltzmann fluid simulations with and without level sets[C]//Proceedings of Workshop on Vision, Modeling, and Visualization (Stanford, CA). Amsterdam:IOS Press, 2004,199-208.

[4] REEVES W T. Particle systems: a technique for modeling a class of fuzzy objects[C]//Proceeding of the 10th annual conference on computer graphics and interactive techniques. New York, NY: ACM, USA, 1983.

[5] SCHLATTER B. A pedagogical tool using smoothed particle hydrodynamics to model fluid flow past a system of cylinders[R]. Oregon:Oregon State University, 1999.

[6] STAM J. Stable fluids [C]//Proceedings of SIGGRAPH. Los Angeles, CA:ACM Press, 1999,121-128.

[7] NVIDIA C. Programming Guid: version 2.1[C]//Proceeding of Symposium on Interactive Ray Tracing. Los Angeles, USA: IEEE Press, 2008:185-190.

[8] 李建明，吳云龍，遲忠先，等.基于流體模型和GPU 加速的火焰實時仿真[J].系統(tǒng)仿真學報, 2007,19(19):4382-4385.

[9] BRIDSON R. Fluid simulation for computer graphics [M].Wellesley, MA：A K Peters, Ltd, 2008.

[10] FEDKIW R, STAM J. Visual simulation of smoke [C]//Proceedings of SIGGRAPH 2001, New York: ACM, 2011:15-22.

[11] Black body spectrum[EB/OL].http://www.mi.infm.it/manini/dida/BlackBody.html.

[12] MARSALEK L, HAUBER A. High-speed volume ray casting with CUDA[C]//Interactive Ray Tracing, 9-10 Aug., 2008, Los Angeles, CA. IEEE Symposium on, IEEE: IEEE Press, 2008, 185.