基于動態網格細分的煙霧模擬

王繼州,袁雪霞

（1中原工學院 信息商務學院計科系,河南 鄭州 450007;2.鄭州經貿職業學院 計算機系,河南 鄭州 450006）

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

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

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

1 基于物理的流體模擬

1.1 流體動力學方程

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

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

1.2 對流（平流）

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

1.3 離散的壓強方程

采用洛薩索等與離散的分歧和一個八叉樹數據結構的壓力梯度方法，調用高斯散度定理的積分形式產生發散方差：

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

剩余發散的壓力梯度(▽p)最終通過一個對稱矩陣實現。結果表明,梯度計算與應用的直接鄰居網格單元作為擾動的壓力位置,只要標準的核心差異在 O(Δx)內時，系統仍然能夠得到一個一致的近似值。

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

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

2 并行模擬

2.1 分解法

為了讓并行處理在現代化的圖像處理硬件問題中能夠適當地分解，使用了一個具有43個細化的分層網格代替傳統的具有23個細分的八叉樹。通過這種方式更容易求解區域常數，另外，層次結構的深度通常比八叉樹小，從而減少昂貴的內存遍歷。

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

2.2 平流運動的分層網格結構

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

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

圖2中，灰色紋理元素（實心箭頭）代表實際網格，而條紋紋理元素（虛線箭頭）描述父節點的平均值。在紋理金字塔內，活動網格塊的安排被列在一個動態拓撲表中。這里有一行為每一個網格的第一個單元格分配相應的紋理塊坐標和限定水平。網格拓撲完全依據此表定義，由CPU管理,其中也考慮到GPU分配任務線程紋理塊。該網格分配給各個線程,通過閱讀獲得的紋理坐標表進入,通過增加ID線程作為補償。圖2中白色紋理的數字元素代表不使用內存開銷，但允許在網絡拓撲快速重排，在運行時分配內存。另外，它們中的一部分被用來正確處理插值以保持邊界上的良好水平。

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

2.3 壓力分層網格的結構

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

但是，為了實現高性能有幾個限制必須保留，例如數據元素必須按順序訪問線程。由于這些限制的存在，壓力值的存儲與數組分層網格結構之間需要進行映射。數組分層網格結構的映射關系如圖4所示。

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

2.4 動態網格細化

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

改進和粗化過程以并行方式運行，但是GPU本身不能管理自己的資源，數據必須被轉移到CPU以管理網格。基于這個原因，每個線程塊通知CPU帶著與被細分的細胞或者被粗化的塊的數據。有了這些數據，CPU能夠在下一個時間步長中重組動態拓撲表以及給線程塊分配任務。對約1283格的尺寸大的測量表明，與不需要這些傳輸的規則笛卡爾網格相比，作為結果的開銷可以忽略不計（小于1%）。

3 模擬煙霧

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

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

[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.