王繼州,袁雪霞

（1中原工學(xué)院 信息商務(wù)學(xué)院計科系,河南 鄭州 450007;2.鄭州經(jīng)貿(mào)職業(yè)學(xué)院 計算機系,河南 鄭州 450006）

在計算機圖形學(xué)動畫特效中，具有豐富細(xì)節(jié)流體(如煙霧或流水等)仍然是一個具有挑戰(zhàn)性的任務(wù)，這是因為流體運動具有復(fù)雜性。現(xiàn)實的動畫需要豐富的視覺細(xì)節(jié) (如高分辨率的流體運動和煙霧現(xiàn)象)以及與環(huán)境的完美融合和高超的照明技術(shù)。然而現(xiàn)有方法的計算復(fù)雜度都非常高,目前所有GPU的流體模擬都依賴于純粹的基于粒子模塊化的方法(如粒子流體動力學(xué)或依賴規(guī)則網(wǎng)格模型)，比較適合快速處理架構(gòu)。規(guī)則網(wǎng)格的主要缺點是在空操作或無用操作上浪費時間。本文方法適用于計算機繪圖中的實際視覺效果，研究重點集中于外觀的表現(xiàn)上，而不考慮嚴(yán)格的物理正確性。

在計算流體力學(xué)方面，有自適應(yīng)網(wǎng)格加密(AMR),如伯杰等[1]發(fā)展了一個基于圖形上下文的自適應(yīng)八叉樹數(shù)據(jù)結(jié)構(gòu)的離散化方法。另一種方式是純粹的以流體模擬粒子為基礎(chǔ)的類似于SPH[2]的方法。對GPU的流體模擬已有參考文獻[3]和參考文獻[4]的SPH方法。哈里斯等人[5]介紹了3D紋理，以模擬定期笛卡爾網(wǎng)格云動態(tài)。

本文將煙霧看做是特殊的流體，提出了一個解決方法，充分利用了硬件的具體性能優(yōu)勢，同時減少了復(fù)雜的層次遍歷。流體數(shù)量極大地影響著整體性能，因此針對平流的速率存儲和泊松求解的壓力存儲之間的不同，采用了不同的內(nèi)存訪問模式。

1 基于物理的流體模擬

1.1 流體動力學(xué)方程

物理模型上的無粘性、不可壓縮的Navier-Stokes方程的質(zhì)量和動量守恒方程式為：

該方程描述了液體流動的連續(xù)性，其中u表示該速度場，ρ表示密度，p表示壓力場，f表示外力作用。

1.2 對流（平流）

動量方程已細(xì)分為一個平流過程和一個壓力求解過程。使用從塞爾等人的改進麥科馬克方法[6]進行誤差補償與校正。該方法可以簡單地用在一個自適應(yīng)網(wǎng)格上，只要確定了粒子軌跡的目的，三線性插值即可用于處理正確的邊界。

1.3 離散的壓強方程

采用洛薩索等與離散的分歧和一個八叉樹數(shù)據(jù)結(jié)構(gòu)的壓力梯度方法，調(diào)用高斯散度定理的積分形式產(chǎn)生發(fā)散方差：

這里n是一個向外的單位法線，Aface是表面積單位，Vcell是體積單位。同樣地，該定理是適用于壓力的拉普拉斯算子，它實際上是一個div(grad(p))項:

剩余發(fā)散的壓力梯度(▽p)最終通過一個對稱矩陣實現(xiàn)。結(jié)果表明,梯度計算與應(yīng)用的直接鄰居網(wǎng)格單元作為擾動的壓力位置,只要標(biāo)準(zhǔn)的核心差異在 O(Δx)內(nèi)時，系統(tǒng)仍然能夠得到一個一致的近似值。

圖1顯示了一個二維的例子。壓力梯度大單元為：

圖1 壓力對八叉樹的離散化

2 并行模擬

2.1 分解法

為了讓并行處理在現(xiàn)代化的圖像處理硬件問題中能夠適當(dāng)?shù)胤纸猓褂昧艘粋€具有43個細(xì)化的分層網(wǎng)格代替?zhèn)鹘y(tǒng)的具有23個細(xì)分的八叉樹。通過這種方式更容易求解區(qū)域常數(shù)，另外，層次結(jié)構(gòu)的深度通常比八叉樹小，從而減少昂貴的內(nèi)存遍歷。

八叉樹方案在實施時會遇到網(wǎng)格占用率不良的問題。例如，當(dāng)被分配的一個線程塊有8個節(jié)點或在一個復(fù)雜任務(wù)的分配方案中進行繪制曲線操作時,隨著網(wǎng)格的逐步細(xì)化，64個線程塊會被分配到一個單一亞網(wǎng)格里面。解決的辦法是，對于每一個線程計算網(wǎng)格單元的中間值u*，把其平流速度和壓力值p作為下一次的SIMD方式進行推測,結(jié)果，包含一個子網(wǎng)格的父網(wǎng)格也同時被模擬,這樣就避免了復(fù)雜的分支操作。這樣處理的時間花費只占總花費的1.5%左右。

2.2 平流運動的分層網(wǎng)格結(jié)構(gòu)

麥科馬克或半拉格朗日平流步法意味著對有效處理的分層網(wǎng)格系統(tǒng)的數(shù)據(jù)結(jié)構(gòu)有幾個要求。該主要操作是隨機的，就像一個空間內(nèi)讀寫訪問局域網(wǎng)中的每個單元格周圍目的地的粒子軌跡。二階運算三線性插值的速度在任意的網(wǎng)格中都是一樣的，計算速度比較均勻。這些觀測值促使紋理存儲讀出速度快，因為讀訪問得益于和本地訪問一樣的紋理緩存速度；此外，紋理單元也支持快速三線性插值。

其基本思想是使用一個三維紋理容量的金字塔映射網(wǎng)格層次結(jié)構(gòu)的內(nèi)存。為了與一個三維貼圖過濾形成對照，只有那些紋理元素在對應(yīng)的單元格網(wǎng)域包含有效的當(dāng)前模擬信息。圖2所示為一個二維的分層網(wǎng)格例子，分層網(wǎng)格擁有一個交錯安排和相應(yīng)的紋理布局。除了3個速度分量,還有一個類型字段存儲在Alpha通道，用以提供輔助信息的網(wǎng)格單元(如子網(wǎng)格的牽制或流入以及固體邊界條件等)。

圖2中，灰色紋理元素（實心箭頭）代表實際網(wǎng)格，而條紋紋理元素（虛線箭頭）描述父節(jié)點的平均值。在紋理金字塔內(nèi)，活動網(wǎng)格塊的安排被列在一個動態(tài)拓?fù)浔碇?。這里有一行為每一個網(wǎng)格的第一個單元格分配相應(yīng)的紋理塊坐標(biāo)和限定水平。網(wǎng)格拓?fù)渫耆罁?jù)此表定義，由CPU管理,其中也考慮到GPU分配任務(wù)線程紋理塊。該網(wǎng)格分配給各個線程,通過閱讀獲得的紋理坐標(biāo)表進入,通過增加ID線程作為補償。圖2中白色紋理的數(shù)字元素代表不使用內(nèi)存開銷，但允許在網(wǎng)絡(luò)拓?fù)淇焖僦嘏?，在運行時分配內(nèi)存。另外，它們中的一部分被用來正確處理插值以保持邊界上的良好水平。

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

2.3 壓力分層網(wǎng)格的結(jié)構(gòu)

以紋理為基礎(chǔ)的數(shù)據(jù)存儲的缺點是缺乏直接寫訪問，必須通過設(shè)備到設(shè)備復(fù)制整個網(wǎng)格的操作。在平流步驟，快速讀取操作的優(yōu)勢緩解了這個問題。例如，聯(lián)麥科馬克與第二階龍格-庫塔方法意味著每個網(wǎng)格單元的總的讀操作數(shù)為5，但是寫操作數(shù)只有一個。在迭代求解器的情況下，寫操作具有更多的權(quán)值。一個非自適應(yīng)笛卡爾網(wǎng)格實驗顯示，基于紋理的方法顯然慢于直接進入整體設(shè)備內(nèi)存與聯(lián)合讀寫的模式。

但是，為了實現(xiàn)高性能有幾個限制必須保留，例如數(shù)據(jù)元素必須按順序訪問線程。由于這些限制的存在，壓力值的存儲與數(shù)組分層網(wǎng)格結(jié)構(gòu)之間需要進行映射。數(shù)組分層網(wǎng)格結(jié)構(gòu)的映射關(guān)系如圖4所示。

圖4中，除了在直接鄰節(jié)點位置的細(xì)胞,矩陣中幾乎一行中所有的元素都為零，典型的迭代求解器要求這些非零元素從內(nèi)存讀。對于恒定分辨率的塊，壓力值可以聯(lián)合起來讀到共享內(nèi)存作進一步處理，但常規(guī)模式打破了邊界的分辨率。為了減少過多的非聯(lián)合訪問，計算細(xì)胞的每個面的壓力值（黑色）的方法通過存儲放在一個額外的數(shù)組中進行。使用它們?nèi)〈鷱募?xì)分細(xì)胞中讀取壓力值，此過程隨著壓力梯度的離散化而愈加完美。

2.4 動態(tài)網(wǎng)格細(xì)化

使用三種常見的標(biāo)準(zhǔn)來決定在下一時間步長中一個細(xì)胞是否精確而塊是否粗造。在固體障礙下，當(dāng)流體與物體相互作用時，網(wǎng)格足夠精確地捕捉場景幾何，分解視覺中的人工物體。第二個標(biāo)準(zhǔn)集中在高漩渦區(qū)域。當(dāng)大小達(dá)到臨界值的時候，網(wǎng)格是精確的。第三個規(guī)則細(xì)化網(wǎng)格在一個煙密度帶內(nèi)。下限不包括從高精度鄰域中的非可見密度值，而對于非常密集地區(qū)的上限這是完全不透明的，只貢獻了一些視覺細(xì)節(jié)。圖5顯示了漩渦跟蹤動態(tài)網(wǎng)格細(xì)化的例子。

改進和粗化過程以并行方式運行，但是GPU本身不能管理自己的資源，數(shù)據(jù)必須被轉(zhuǎn)移到CPU以管理網(wǎng)格。基于這個原因，每個線程塊通知CPU帶著與被細(xì)分的細(xì)胞或者被粗化的塊的數(shù)據(jù)。有了這些數(shù)據(jù)，CPU能夠在下一個時間步長中重組動態(tài)拓?fù)浔硪约敖o線程塊分配任務(wù)。對約1283格的尺寸大的測量表明，與不需要這些傳輸?shù)囊?guī)則笛卡爾網(wǎng)格相比，作為結(jié)果的開銷可以忽略不計（小于1%）。

3 模擬煙霧

運用集成可視化自適應(yīng)流體模擬模塊進行模擬應(yīng)用。輻射度網(wǎng)格的射線追蹤K對彼此平行的網(wǎng)格單元格中最明亮的光源，輻射度是陰影值的總和加上其余斑塊的輻射度；密度和顏色值存儲在一個額外3D紋理中供渲染使用；用標(biāo)準(zhǔn)的光線投射在GPU上的程序來呈現(xiàn)密度的三維體積和照明的信息。模擬的效果如圖6所示。

本文采取了動態(tài)網(wǎng)格細(xì)分在平行的SIMD圖形硬件上的流體模擬，減少了不規(guī)則的八叉樹的執(zhí)行情況，以適應(yīng)硬件的性能限制。納入輻射度為基礎(chǔ)的可視化模塊的模擬，取得了與煙一起互動的真實感渲染幀速率。與同等笛卡爾網(wǎng)格比較，開銷是不可避免的，但如果改進的策略是精心挑選的，加速比是合理的。未來可以在多GPU系統(tǒng)上進一步擴展工作。

[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]MüLLER 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.Adaptivesampling 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.