基于射線穿透法的GPU 并行階梯型有限差分網格生成算法*

李 平，麻鐵昌，許香照，馬天寶

（北京理工大學機電學院，北京 100081）

爆炸力學是武器裝備研發，工業安全防護等問題的理論基礎。爆炸的特點為在極短的時間內，爆炸場中各個物理量發生巨大的變化。爆炸過程的研究十分復雜，通常伴隨有化學變化、物理變化和核變化等。因此在理論研究和實驗研究中都存在諸多困難。數值計算是研究爆炸與沖擊問題的重要手段之一。計算爆炸力學是以計算機為工具，探索爆炸力學的規律，豐富力學數據，解決爆炸力學問題。隨著計算機科學的飛速發展，數值模擬以其精度高，經濟成本低，揭示規律清晰完整等諸多優點越來越成為人們的研究熱點。爆炸力學的數值模擬按照其采用的坐標形式可以分為歐拉方法和拉格朗日方法。基于歐拉坐標的有限差分方法是爆炸與沖擊問題數值模擬的常用方法[1]，其特點為在固定的階梯型網格上進行有限差分計算，網格位置和形狀不會隨著計算而改變，克服了拉格朗日方法中可能出現的網格畸變問題，因此能夠很好的模擬多物質的大變形問題[2-3]。本文涉及的有限差分計算中使用的網格數據是一種階梯型笛卡爾網格，其特點是使用正六面體對整個區域進行離散，不同物質被離散為不同屬性的單元，不同屬性單元構成的不同區域之間的交界面為階梯形[4-5]。相比于傳統貼體網格，階梯型有限差分網格的生成速度快，無需手動劃分區域，生成過程自動化程度高。階梯型有限差分網格是一種結構網格，每個單元之間的連接關系固定并不隨著計算而變化，可以使用像所有結構網格一樣的IJK 下標來訪問網格單元和節點。它適用于所有的有限差分離散方法。

階梯型有限差分網格生成方法主要包含兩種：一種是射線穿透法，這種方法即根據網格步長作出一系列射線直接與實體模型進行求交計算[6]，這種方法的優點是思路清晰，步驟簡單；但是由于空間中的射線與幾何模型的求交計算會消耗大量的時間，因此算法運行效率不高。切片法的思路借鑒了快速成型加工領域的方法[7-8]，用一系列互相平行的平面截取幾何模型，得到一系列二維平面中的封閉輪廓線，進而將三維空間中的幾何問題轉換成二維平面上的問題；此方法省去了部分幾何判斷、查詢占用的時間，網格生成效率相對較高。但是此種算法執行過程中網格線的生成和介質屬性映射兩個主要步驟耦合在一起，算法執行的穩定性和并行能力都有待提高。在三維有限差分計算中，網格尺寸越小，網格構成的計算空間與實際物理空間越接近。因此為了保證數值模擬的精度，進行三維有限差分計算通常需要至少千萬量級的網格[9]，如何高效的生成數量如此巨大的數值網格一直是研究熱點之一。MacGrillivray[10]使用高效的數據存儲方法和高效的空間線面相交判別方法實現了萬億量級的網格生成，網格生成過程內存管理方法合理，使整個網格生成過程可以在小型工作站中完成。Berens 等[11]提出了非均勻笛卡爾網格的生成方法，詳細描述了根據計算需求而對網格尺寸的選取方法。以上都是串行的網格生成方法。并行計算是提高網格生成效率的最常用的方法。在計算集群中運行的程序最常使用的是利用MPI(message passing interface)庫進行并行計算，在小型工作站中常用的是OpenMP 庫實現共享內存的并行計算，Ning 等[12]使用多線程技術實現了并行階梯型有限差分網格生成，但是其采用的是切片法生成網格，單個幾何體分別由單個線程執行計算，并行化程度較低，對大型的網格數據生成效率不高。Ishida 等[13]提出采用OpenMP 并行生成自適應笛卡爾網格。Foteinos 等[14]將并行計算應用在“圖到網格”的網格生成算法中，采用分布式集群計算機對算法進行加速。

最近，GPU （graphic processing unit）并行技術吸引了各個領域研究者的注意，不僅僅被應用在圖形圖像相關算法，還被廣泛應用于各種通用計算中。相對于CPU （central process unit）并行，GPU 對于大規模密集型數據的并行計算擁有巨大的優勢。許多學者也對將GPU 應用到網格生成領域做了相關研究。Qi 等[15]利用GPU 并行實現了二維Delaunay 非結構網格的生成。Park 等[16]提出了基于GPU 生成3D 自適應笛卡爾網格的方法。Schwarz 等[17]提出了一種方法使用GPU 來生成八叉樹結構，并將其應用于計算機圖形學中的像素化顯示。但是他們使用的基于三角形并行化的方法生成的單元數量較少，不適用于進行有限差分計算。根據文獻調研情況，目前運用GPU 并行技術生成三維笛卡爾網格的研究很少。本文中提出的并行算法能夠應用GPU 并行技術在短時間內生成大規模的階梯型有限差分網格，并且對多個擁有復雜外形的幾何體的大規模網格生成過程有很高的效率和準確性。

本文中，選用STL （stereolithography）文件來存儲幾何實體信息，并用其作為算法的輸入文件。在STL 文件中，幾何實體被離散為三角面片的形式，所有三角面片的三維坐標和法向量都存貯其中[18]。在一些切片算法中，STL 文件需要過濾冗余幾何信息并重構出拓撲信息。對于射線穿透法，STL 文件中的信息不再需要任何的處理，STL 文件中的三角形集合T 可以直接作為算法的輸入。階梯型有限差分網格生成算法的核心部分就是由T 作為輸入生成整個計算域大小的數據場F，F 中每一坐標點的值為介質標志。

1 射線穿透法

1.1 計算域設置

（1）計算域的邊界線和計算域內各個幾何體的AABB 包圍盒的邊線同時也是網格線。三維空間中幾何對象的AABB 包圍盒為包含該對象，且邊平行于坐標軸的最小六面體。

（2）計算域內同一維度方向上的網格尺寸s 相同。

條件1 和條件2 同時限定了整個計算域中網格尺寸大小相同。同時，條件1 消除了在幾何體AABB 包圍盒邊緣存在的狹窄網格單元。為了滿足上述兩個條件，網格單元尺寸需要根據輸入的尺寸進行優化。如下列方程分別為X、Y 和Z 方向上的網格尺寸函數：

式中：sopt為待優化的網格尺寸，nx、ny和nz分別為X、Y 和Z 方向上的網格單元數，ΔXi為計算域被幾何體AABB 包圍盒邊界線分割出的子計算區域。

圖1 三維計算域，內部幾何體和幾何體局部圖Fig. 1 Three-dimensional computational domain, the geometries in the domain and the details of the geometry

優化方程(1)～(3)中的sopt，計算使得3 個方程同時取得最小值同時又不大于初始的網格尺寸s 的sopt，即為優化后的網格單元尺寸。

獲得優化的網格單元尺寸后，可以按照下列公式生成整個計算域中網格線：

1.2 射線與幾何體表面的相交判斷

對于任意一條射線穿過計算區域，可能與一些三角面片相交，也可能不相交。如圖2 中所示，由投影平面出發的射線穿過計算區域與幾何體相交，射線簇與幾何體的相交計算的本質是空間中射線r 與三角面片的相交計算，V1、V2和V3為三角面片3 個頂點。為了計算射線與三角面片的交點，首先要判斷射線與哪些面片相交。為了減少機器誤差，這里采用一種沒有除法計算的判斷方法。對于一條射線和一個三角面片，通過下列方程計算射線與三角面片第i 條邊的相對位置Oi：

圖2 射線穿透法原理（以XY 平面為射線投射平面，Z 軸方向為射線方向）Fig. 2 Principle of ray casting method (set the XY plane as the projection plane, and the Z dimension as the ray direction)

Oi表示射線與三角面片的一條邊的相對位置，如果Oi＞0，則射線r 在邊V1V2的左側，如果Oi＜0，則射線r 在邊V1V2的右側，Oi=0 時，射線與該邊相交。分別計算三角面片三條邊與射線的Oi值，如果Oi全部為正或全部為負值，則射線在3 條邊的同側，表明射線穿過了三角形，否則射線與三角形不相交。

通過一條穿過幾何體的射線將會計算出多個交點，將這些交點按照射線方向排序，得到交點集合P。在三維空間中，直線與三角面片的位置關系可分為7 類：穿過、錯過、穿點、穿邊、平行、共面不相交和共面相交。很顯然，穿過和錯過可以由上述方法判斷并計算交點。穿點是直線穿過三角面片任意一頂點，同時由于三角面片組成的是閉合幾何體，此直線一定也穿過其他三角面片的頂點。在這種情況下，在判別式中，O1、O2、O3中有兩個為0，另外一個不為0。交點計算公式仍然成立，但是會計算出多個相同的交點，在排序后所得的集合P 中，相同的交點只保留一個。對于穿邊的情況，依據相同的原理，O1、O2、O3中有一個為0，另外兩個同號。交點計算公式依然成立，在這個位置會有兩個交點，在P 中保留其中的一個。在直線與三角面片平行時，O1、O2、O3一定不會符號相同，因此會判定為錯過，不計算交點。當直線與三角面片共面不相交時，O1、O2、O3全部等于0，會判定為不相交，不計算交點。當直線與三角面片共面并穿過面片時，說明這個面片為幾何物體的邊界，邊界上不會有網格點，因此判別式可以正確將其判斷為不相交。

綜上所述，在穿點與穿邊的情況中，會出現一項或兩項Oi為0，這種情況同樣應視為相交，計算出的交點可能會在其相鄰三角面片中再次出現，使用上述原理僅保留一個交點即可。

2 GPU 并行網格生成策略

三維空間中幾何圖元的搜索、查找和相交計算是射線穿透法的主要耗時部分。同時，這些計算都可以轉化為某一簡單計算流程的疊加。因此，將這些耗時巨大的計算轉移到高度并行化的GPU 中并行執行可以大大提高算法效率。GPU 采用單程序多數據（single program/multiple data，SPMD）模型，其目的就是通過其內部大量的線程、線程束、線程塊和線程網格等并行層級來執行大量的、比較簡單的計算任務。GPU 并行算法的效率主要由2 個部分決定：（1）單個線程中計算時間；（2）主機內存與GPU 內存數據交換時間。

下面，針對上述兩部分分別優化并行算法。

2.1 幾何數據傳輸優化

下面以一個單獨幾何體G 為例說明數據傳輸過程。作為輸入，STL 中的三角面片數據被存儲到內存序列中，以數組T 表示。幾何體G 的AABB 包圍盒被一系列水平平面二次劃分為M 個包圍盒子區域。隨后，將三角形序列T 按照其在空間中所屬的包圍盒子區域分解為M 個子三角形序列，將一個包圍盒子區域的信息SubAABB 和與其對應的子三角形序列Tsub作為一個批次（batch），則輸入數據被劃分為M 個批次，每個批次都包含了處理每個批次中所有射線計算所需的全部數據。隨后，將M 個批次數據逐個傳入GPU 內存中執行并行計算。每一個線程執行計算后的輸出數據為一條射線r 方向上的網格個數序列Cr，Cr中元素的數量表示射線r 方向上的總網格數被幾何體分割出的段數，每段的網格數由Cr中一個元素表示。幾何體G 的AABB 包圍盒內全部射線對應的網格個數序列的集合可以表示為C={Cr1,Cr2,···,Cri,···,Crn}，C 中包含的數據即可表示幾何體G 的網格數據。對于計算域中的每一個幾何體，重復上述過程，將每個幾何體計算所得的網格個數序列集合C 合并在一起，即可得到整個計算域中的階梯型有限差分網格劃分結果。上述數據轉換與計算過程如圖3 所示。圖中左上部分表示待計算區域及計算區域中的幾何體，其中幾何體以三角面片形式被表示。右上圖表示幾何體的AABB 包圍盒構成的子計算區域。圖中下半部分表示幾何數據及邊界條件等被組織為數據批次，并逐批次傳如GPU 進行計算的并行計算流程。

圖3 三維幾何數據轉換與并行計算過程圖Fig. 3 Diagram of three-dimensional geometric data conversion and parallel computing process

2.2 線程內射線生成

GPU 并行計算開始執行之后，大量的線程會同時執行相同的指令，這些指令的設計直接影響并行算法的并行執行效率。對于域中的每個幾何體，網格線和射線的分布可以通過AABB 包圍盒及包圍盒子區域的對角點坐標(x1, y1)和(x2, y2)與網格尺寸sx和sy來確定。當GPU 執行某個計算任務時，通常需要將輸入數據從主機內存復制到GPU 中的設備內存。為了避免主機和設備之間數據傳輸的時間消耗和GPU 內存中的容量消耗，可以在GPU 中的每個線程中自動生成網格線和射線的起點坐標。

假設幾何模型G 的AABB 包圍盒可以表示為兩個對角點坐標(x1, y1)和(x2, y2)。X 與Y 方向上的網格尺寸可以分別表示為sx和sy。在GPU 的每個線程中，當前線程開始計算時，可以通過下列公式生成射線的起點坐標(xr, yr)：

式中：Nindex和nx分別表示GPU 中當前線程的索引序號和當前包圍盒子區域中X 坐標方向的網格數。Nindex和nx分別可以由下式計算：

生成射線起點坐標以后，應用1.2 節中所述方法判斷射線與哪些三角形相交并計算出交點坐標。將各個交點Z 坐標連同計算域Z 方向的最大、最小值按照由小到大的次序排列，可得到Z 坐標序列(0＜i＜2n)。由于計算域中幾何體都是封閉的，交點個數必定為偶數，以2n 表示。則對于射線r 方向上所有的網格單元，第i 個網格單元的屬性可以按照下列公式賦值：

式中：sz為計算域Z 方向上的網格尺寸。這一步驟通常稱為介質屬性映射。本節所述計算全部由單個線程獨立完成，計算流程如圖4 所示。

圖4 單個GPU 線程內數據計算流程Fig. 4 Data computing process in a single GPU thread

3 數值試驗結果和討論

應用文中提出的基于射線穿透法的GPU 并行網格生成方法，使用Visual C++和Nvida CUDA[19]編制基于GPU 的并行階梯型有限差分網格生成程序并對程序進行性能測試。所用測試服務器配置如下：CPU 型號，Intel Xeon E5-2650 v2（2.6 GHz）；GPU 型號，Nvidia Quadro K2200；顯存容量，4 GB。

3.1 網格生成結果與可視化

如圖5(a)所示橋梁模型包含3 種材質，由多個幾何體組合而成。模型總計包含24 184 個三角面片，被劃分為1.15×109個網格單元。總計用時21.45 s。圖5(b)為網格圖，由放大后的局部網格圖可以清晰地看出，三角面片構成的幾何模型被轉變為由不同材質的六面體單元構成的階梯型有限差分網格。

為了驗證并行算法的網格生成效率，將上述模型分別離散為不同網格規模的階梯型有限差分網格，并對比網格生成時間與傳統串行射線穿透法的網格生成時間。如圖6 所示，橋梁模型分別被離散為1.44×108、1.15×109和9.2×109個網格單元。可以看出，對同一計算模型應用本文提出的GPU 并行算法的網格生成效率遠高于傳統串行CPU 算法的執行效率。傳統串行算法在生成網格規模達到1×1010數量級的時候耗費的時間超過2 000 s，這無疑大大影響了建模與計算時的靈活性，在很多需要多次調整計算域參數的數值模擬中是難以接受的。另外，表示這種數量級的網格數據通常可以達到4 GB 以上，本文中提出的分批次的數據處理方法使得算法能夠處理的數據規模不依賴于GPU 內存大小，使程序在常見的4 GB 顯存容量的GPU 中可以高效執行。CPU 的硬件發展已趨于穩定，短時間內難以有巨大的提升。本文中提出的并行算法可以解決這些瓶頸，使得在擁有一顆普通GPU 的PC 機上就足以進行1×1010數量級的網格生成。對于不同的初始幾何模型，網格生成算法的效率應該不受模型的復雜程度的影響。表1 中展示了3 種不同的初始幾何模型分別包含19 202、78 354 和95 062 個三角面片，分別生成相同規模（1×109個網格單元）的三維笛卡爾網格。表中分別給出了CPU 串行算法和GPU 并行算法的生成時間。可以看出，對于復雜程度不同的幾何模型生成相同規模的網格，本文提出的并行算法的效率是傳統串行算法的8～11 倍。

圖5 橋梁模型及其階梯型有限差分網格生成結果圖及細節放大圖Fig. 5 A bridge model, the finite difference mesh generated and the details of the mesh

圖6 本文提出的并行算法與傳統算法的網格生成時間比較折線圖Fig. 6 Mesh generation time comparison curve between the proposed parallel algorithm and the traditional algorithm

表1 不同模型生成相同數量網格單元(1×109)的執行時間Table 1 Generation times of different models with the cell number of 1×109

綜上所述，本文中提出的并行算法的網格生成效率是傳統串行算法的8～11 倍。隨著網格數量級和模型復雜度的增高，并行算法節省的網格生成時間越來越多。

3.2 數值模擬結果

工廠廠房內的爆炸是工業生產中危害巨大的一種安全事故。廠房合理的結構設計是降低廠房爆炸造成的損失的有效手段。在這一節，通過文中提出的并行階梯型有限差分網格生成方法，對某工廠廠房進行網格劃分，并對廠房內爆炸問題進行數值模擬研究。如圖7 所示，計算域中包括墻壁、立柱、屋頂、炸藥和空氣5 種介質。廠房為一座兩層建筑，上下兩層有樓梯連接，64 kg TNT 炸藥位于廠房一層中間位置。為了應用有限差分法進行數值模擬，將計算域離散為三維階梯型有限差分網格，如圖8 所示。為了清晰顯示網格生成結果，空氣網格被隱藏掉，圖8 中的炸藥，墻壁、立柱和屋頂4 種介質分別以4 種顏色顯示。網格生成單元總數為1.5×108個。并行生成消耗總時間為2.47 s。

圖8 某廠房三維階梯型有限差分網格Fig. 8 Three-dimensional finite difference mesh of a factory model

圖7 某廠房三維幾何模型Fig. 7 A three-dimensional factory model

應用北京理工大學爆炸與科學國家重點實驗室自主開發的三維多介質流體動力學仿真軟件PMMIC-3D[20]對上述網格數據進行數值模擬，得到計算結果如圖9 所示。圖9（a）中左側為爆炸產生后3 個有代表性的時刻的三維可視化結果，右側為與左側對應的時刻的二維剖面可視化結果。從圖9（a）可以清晰地看到爆炸在廠房一層發生后，形成沖擊波（圖中以白色表示）并膨脹擴散到二層的過程。圖9（b）為廠房二層2 維剖面圖，以壓力的變化為可視化屬性顯示了在4 個典型的時刻沖擊波在二層傳播過程。

圖9 廠房爆炸數值模擬結果圖Fig. 9 Numerical simulation results of factory explosion

為了更進一步分析數值模擬的準確性，選取了廠房二層距樓梯口10 m 處的點作為關鍵點M 進行超壓測試。記錄點M 超壓隨時間的變化如圖10 所示。可以看到，點M 處初始超壓為0，隨后沖擊波傳入二樓并傳到點M 處，壓力值開始上升達到第一個峰值。沖擊波傳播到墻壁后反射再次到達點M，壓力值達到第二個峰值。隨后，隨著不斷的反射，M 點超壓值逐漸減小直至衰減為0。可以看出數值模擬與沖擊波傳播反射理論和實際經驗相吻合，說明計算域的網格剖分能夠滿足大規模有限差分計算的需要。

圖10 廠房二層距離樓梯口10 m 處關鍵點的超壓變化曲線Fig. 10 Change of overpressure with time at a key point which is 10 m from stairway entrance on the second floor

4 結 論

本文提出了一種基于傳統射線穿透法的GPU 并行階梯型有限差分網格生成方法。在這種并行方法中，提出了一種分批次的數據傳輸策略，使得算法能夠處理的數據規模不依賴于GPU 內存大小，打破了硬件對網格劃分規模的限制，平衡了數據傳輸效率和網格生成規模之間的關系。為了減少數據傳輸量，本文提出的并行算法可以由GPU 線程相互獨立的生成射線起點坐標，射線相交計算在GPU 的每個線程中獨自計算，進一步提高了并行算法的執行效率，通過數值試驗的對比可以看出，并行算法的執行效率是傳統射線穿透法執行效率的8～11 倍，并且隨著計算規模的提升，并行算法的加速比有上升趨勢。最后，通過有限差分計算實例驗證了應用并行算法生成的階梯型有限差分網格能夠滿足基于有限差分的數值模擬需求，得到了與理論和實驗一致的數值模擬結果。