摘  要： 提出了一種基于動(dòng)態(tài)網(wǎng)格細(xì)分的煙霧模擬方法，該方法主要采取圖形設(shè)備上的動(dòng)態(tài)網(wǎng)格管理對(duì)煙霧進(jìn)行并行處理以達(dá)到泊松方程的迭代求解。為了實(shí)現(xiàn)高性能,利用高速緩存以提高存取權(quán)限和適應(yīng)硬件的能力。實(shí)驗(yàn)結(jié)果表明，該方法能夠?qū)崿F(xiàn)比較快速的模擬，結(jié)果比較令人滿意。
關(guān)鍵詞： 流體模擬； 煙霧模擬； 網(wǎng)格細(xì)分； 動(dòng)態(tài)網(wǎng)格

    在計(jì)算機(jī)圖形學(xué)動(dòng)畫特效中，具有豐富細(xì)節(jié)流體(如煙霧或流水等)仍然是一個(gè)具有挑戰(zhàn)性的任務(wù)，這是因?yàn)榱黧w運(yùn)動(dòng)具有復(fù)雜性。現(xiàn)實(shí)的動(dòng)畫需要豐富的視覺細(xì)節(jié)(如高分辨率的流體運(yùn)動(dòng)和煙霧現(xiàn)象)以及與環(huán)境的完美融合和高超的照明技術(shù)。然而現(xiàn)有方法的計(jì)算復(fù)雜度都非常高,目前所有GPU的流體模擬都依賴于純粹的基于粒子模塊化的方法(如粒子流體動(dòng)力學(xué)或依賴規(guī)則網(wǎng)格模型)，比較適合快速處理架構(gòu)。規(guī)則網(wǎng)格的主要缺點(diǎn)是在空操作或無用操作上浪費(fèi)時(shí)間。本文方法適用于計(jì)算機(jī)繪圖中的實(shí)際視覺效果，研究重點(diǎn)集中于外觀的表現(xiàn)上，而不考慮嚴(yán)格的物理正確性。
    在計(jì)算流體力學(xué)方面，有自適應(yīng)網(wǎng)格加密(AMR),如伯杰等[1]發(fā)展了一個(gè)基于圖形上下文的自適應(yīng)八叉樹數(shù)據(jù)結(jié)構(gòu)的離散化方法。另一種方式是純粹的以流體模擬粒子為基礎(chǔ)的類似于SPH [2]的方法。對(duì)GPU的流體模擬已有參考文獻(xiàn)[3]和參考文獻(xiàn)[4]的SPH方法。哈里斯等人[5]介紹了3D紋理，以模擬定期笛卡爾網(wǎng)格云動(dòng)態(tài)。
    本文將煙霧看做是特殊的流體，提出了一個(gè)解決方法，充分利用了硬件的具體性能優(yōu)勢(shì)，同時(shí)減少了復(fù)雜的層次遍歷。流體數(shù)量極大地影響著整體性能，因此針對(duì)平流的速率存儲(chǔ)和泊松求解的壓力存儲(chǔ)之間的不同，采用了不同的內(nèi)存訪問模式。

1基于物理的流體模擬
1.1流體動(dòng)力學(xué)方程
    物理模型上的無粘性、不可壓縮的Navier-Stokes方程的質(zhì)量和動(dòng)量守恒方程式為：
  
1.2 對(duì)流（平流）
    動(dòng)量方程已細(xì)分為一個(gè)平流過程和一個(gè)壓力求解過程。使用從塞爾等人的改進(jìn)麥科馬克方法[6]進(jìn)行誤差補(bǔ)償與校正。該方法可以簡單地用在一個(gè)自適應(yīng)網(wǎng)格上，只要確定了粒子軌跡的目的，三線性插值即可用于處理正確的邊界。
1.3離散的壓強(qiáng)方程
    采用洛薩索等與離散的分歧和一個(gè)八叉樹數(shù)據(jù)結(jié)構(gòu)的壓力梯度方法，調(diào)用高斯散度定理的積分形式產(chǎn)生發(fā)散方差：
 
    為了讓并行處理在現(xiàn)代化的圖像處理硬件問題中能夠適當(dāng)?shù)胤纸?，使用了一個(gè)具有43個(gè)細(xì)化的分層網(wǎng)格代替?zhèn)鹘y(tǒng)的具有23個(gè)細(xì)分的八叉樹。通過這種方式更容易求解區(qū)域常數(shù)，另外，層次結(jié)構(gòu)的深度通常比八叉樹小，從而減少昂貴的內(nèi)存遍歷。
　八叉樹方案在實(shí)施時(shí)會(huì)遇到網(wǎng)格占用率不良的問題。例如，當(dāng)被分配的一個(gè)線程塊有8個(gè)節(jié)點(diǎn)或在一個(gè)復(fù)雜任務(wù)的分配方案中進(jìn)行繪制曲線操作時(shí),隨著網(wǎng)格的逐步細(xì)化，64個(gè)線程塊會(huì)被分配到一個(gè)單一亞網(wǎng)格里面。解決的辦法是，對(duì)于每一個(gè)線程計(jì)算網(wǎng)格單元的中間值u*，把其平流速度和壓力值p作為下一次的SIMD方式進(jìn)行推測(cè),結(jié)果，包含一個(gè)子網(wǎng)格的父網(wǎng)格也同時(shí)被模擬,這樣就避免了復(fù)雜的分支操作。這樣處理的時(shí)間花費(fèi)只占總花費(fèi)的1.5%左右。
2.2平流運(yùn)動(dòng)的分層網(wǎng)格結(jié)構(gòu)
    麥科馬克或半拉格朗日平流步法意味著對(duì)有效處理的分層網(wǎng)格系統(tǒng)的數(shù)據(jù)結(jié)構(gòu)有幾個(gè)要求。該主要操作是隨機(jī)的，就像一個(gè)空間內(nèi)讀寫訪問局域網(wǎng)中的每個(gè)單元格周圍目的地的粒子軌跡。二階運(yùn)算三線性插值的速度在任意的網(wǎng)格中都是一樣的，計(jì)算速度比較均勻。這些觀測(cè)值促使紋理存儲(chǔ)讀出速度快，因?yàn)樽x訪問得益于和本地訪問一樣的紋理緩存速度；此外，紋理單元也支持快速三線性插值。
    其基本思想是使用一個(gè)三維紋理容量的金字塔映射網(wǎng)格層次結(jié)構(gòu)的內(nèi)存。為了與一個(gè)三維貼圖過濾形成對(duì)照，只有那些紋理元素在對(duì)應(yīng)的單元格網(wǎng)域包含有效的當(dāng)前模擬信息。圖2所示為一個(gè)二維的分層網(wǎng)格例子，分層網(wǎng)格擁有一個(gè)交錯(cuò)安排和相應(yīng)的紋理布局。除了3個(gè)速度分量,還有一個(gè)類型字段存儲(chǔ)在Alpha通道，用以提供輔助信息的網(wǎng)格單元(如子網(wǎng)格的牽制或流入以及固體邊界條件等)。
    圖2中，灰色紋理元素（實(shí)心箭頭）代表實(shí)際網(wǎng)格，而條紋紋理元素（虛線箭頭）描述父節(jié)點(diǎn)的平均值。在紋理金字塔內(nèi)，活動(dòng)網(wǎng)格塊的安排被列在一個(gè)動(dòng)態(tài)拓?fù)浔碇?。這里有一行為每一個(gè)網(wǎng)格的第一個(gè)單元格分配相應(yīng)的紋理塊坐標(biāo)和限定水平。網(wǎng)格拓?fù)渫耆罁?jù)此表定義，由CPU管理,其中也考慮到GPU分配任務(wù)線程紋理塊。該網(wǎng)格分配給各個(gè)線程,通過閱讀獲得的紋理坐標(biāo)表進(jìn)入,通過增加ID線程作為補(bǔ)償。圖2中白色紋理的數(shù)字元素代表不使用內(nèi)存開銷，但允許在網(wǎng)絡(luò)拓?fù)淇焖僦嘏?，在運(yùn)行時(shí)分配內(nèi)存。另外，它們中的一部分被用來正確處理插值以保持邊界上的良好水平。

    三線性插值需要附加的環(huán)境，一個(gè)網(wǎng)格單元不斷變化的處理水平如圖3所示。在粗糙細(xì)胞的情況下，所使用的平均速度來自黑色箭頭的區(qū)域。通過這種方式，在粗略水平條件下進(jìn)行插值就變得單一化。為了優(yōu)化邊界網(wǎng)格單元，分層更新是必要的。線程塊粗格填充子細(xì)胞插值來自于黑色虛線表示的父細(xì)胞。這個(gè)過程允許使用的快速三線性插值在固有的紋理單元內(nèi)，因?yàn)檎_相鄰關(guān)系在兩個(gè)級(jí)別上都得到了保證。

2.3 壓力分層網(wǎng)格的結(jié)構(gòu)
    以紋理為基礎(chǔ)的數(shù)據(jù)存儲(chǔ)的缺點(diǎn)是缺乏直接寫訪問，必須通過設(shè)備到設(shè)備復(fù)制整個(gè)網(wǎng)格的操作。在平流步驟，快速讀取操作的優(yōu)勢(shì)緩解了這個(gè)問題。例如，聯(lián)麥科馬克與第二階龍格-庫塔方法意味著每個(gè)網(wǎng)格單元的總的讀操作數(shù)為5，但是寫操作數(shù)只有一個(gè)。在迭代求解器的情況下，寫操作具有更多的權(quán)值。一個(gè)非自適應(yīng)笛卡爾網(wǎng)格實(shí)驗(yàn)顯示，基于紋理的方法顯然慢于直接進(jìn)入整體設(shè)備內(nèi)存與聯(lián)合讀寫的模式。
    但是，為了實(shí)現(xiàn)高性能有幾個(gè)限制必須保留，例如數(shù)據(jù)元素必須按順序訪問線程。由于這些限制的存在，壓力值的存儲(chǔ)與數(shù)組分層網(wǎng)格結(jié)構(gòu)之間需要進(jìn)行映射。數(shù)組分層網(wǎng)格結(jié)構(gòu)的映射關(guān)系如圖4所示。

    圖4中，除了在直接鄰節(jié)點(diǎn)位置的細(xì)胞,矩陣中幾乎一行中所有的元素都為零，典型的迭代求解器要求這些非零元素從內(nèi)存讀。對(duì)于恒定分辨率的塊，壓力值可以聯(lián)合起來讀到共享內(nèi)存作進(jìn)一步處理，但常規(guī)模式打破了邊界的分辨率。為了減少過多的非聯(lián)合訪問，計(jì)算細(xì)胞的每個(gè)面的壓力值（黑色）的方法通過存儲(chǔ)放在一個(gè)額外的數(shù)組中進(jìn)行。使用它們?nèi)〈鷱募?xì)分細(xì)胞中讀取壓力值，此過程隨著壓力梯度的離散化而愈加完美。
2.4 動(dòng)態(tài)網(wǎng)格細(xì)化
    使用三種常見的標(biāo)準(zhǔn)來決定在下一時(shí)間步長中一個(gè)細(xì)胞是否精確而塊是否粗造。在固體障礙下，當(dāng)流體與物體相互作用時(shí)，網(wǎng)格足夠精確地捕捉場(chǎng)景幾何，分解視覺中的人工物體。第二個(gè)標(biāo)準(zhǔn)集中在高漩渦區(qū)域。當(dāng)大小達(dá)到臨界值的時(shí)候，網(wǎng)格是精確的。第三個(gè)規(guī)則細(xì)化網(wǎng)格在一個(gè)煙密度帶內(nèi)。下限不包括從高精度鄰域中的非可見密度值，而對(duì)于非常密集地區(qū)的上限這是完全不透明的，只貢獻(xiàn)了一些視覺細(xì)節(jié)。圖5顯示了漩渦跟蹤動(dòng)態(tài)網(wǎng)格細(xì)化的例子。

    改進(jìn)和粗化過程以并行方式運(yùn)行，但是GPU本身不能管理自己的資源，數(shù)據(jù)必須被轉(zhuǎn)移到CPU以管理網(wǎng)格。基于這個(gè)原因，每個(gè)線程塊通知CPU帶著與被細(xì)分的細(xì)胞或者被粗化的塊的數(shù)據(jù)。有了這些數(shù)據(jù)，CPU能夠在下一個(gè)時(shí)間步長中重組動(dòng)態(tài)拓?fù)浔硪约敖o線程塊分配任務(wù)。對(duì)約1283格的尺寸大的測(cè)量表明，與不需要這些傳輸?shù)囊?guī)則笛卡爾網(wǎng)格相比，作為結(jié)果的開銷可以忽略不計(jì)（小于1％）。
3 模擬煙霧
    運(yùn)用集成可視化自適應(yīng)流體模擬模塊進(jìn)行模擬應(yīng)用。輻射度網(wǎng)格的射線追蹤K對(duì)彼此平行的網(wǎng)格單元格中最明亮的光源，輻射度是陰影值的總和加上其余斑塊的輻射度；密度和顏色值存儲(chǔ)在一個(gè)額外3D紋理中供渲染使用；用標(biāo)準(zhǔn)的光線投射在GPU上的程序來呈現(xiàn)密度的三維體積和照明的信息。模擬的效果如圖6所示。

    本文采取了動(dòng)態(tài)網(wǎng)格細(xì)分在平行的SIMD圖形硬件上的流體模擬，減少了不規(guī)則的八叉樹的執(zhí)行情況，以適應(yīng)硬件的性能限制。納入輻射度為基礎(chǔ)的可視化模塊的模擬，取得了與煙一起互動(dòng)的真實(shí)感渲染幀速率。與同等笛卡爾網(wǎng)格比較，開銷是不可避免的，但如果改進(jìn)的策略是精心挑選的，加速比是合理的。未來可以在多GPU系統(tǒng)上進(jìn)一步擴(kuò)展工作。
參考文獻(xiàn)
[1] LOSASSO F, GIBOU F, FEDKIW R. Simulating waterand smoke with an octree data structure[C].In SIGGRAPH ′04: ACM SIGGRAPH 2004 Papers,New York, NY, USA,2004:457-462.
[2] MLLER M, CHARYPAR D, GROSS M. Particlebasedfluid simulation for interactive applications. In SCA′03: Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation,Aire-la-Ville, Switzerland,Switzerland, Eurographics Association,2003:154-159.
[3] KOLB A, CUNTZ N. Dynamic particle coupling forgpubased fluid simulation[C].Proc.18th Symposium on Simulation Technique, 2005:722-727.
[4] ZHANG Y, SOLENTHALER B, PAJAROLA R. Adaptive sampling and rendering of fluids on the gpu. In Proceedings Symposium on Point-Based Graphics,2008:137-146.
[5] HARRIS M J, BAXTER W V, SCHEUERMANNT, et al. Simulation of cloud dynamics on graphicshardware. In  HWWS ′03: Proceedings of the ACM SIGGRAPH/EUROGRAPHICS conference on Graphics hardware, Airela Ville, Switzerland, Switzerland, Eurographics Association,2003:92-101.
[6] SELLE A, FEDKIW R, KIM B, LIU Y, et al. An unconditionally stable maccormack method[J]. J. Sci. Comput, 2008,35(2-3):350-371.