基于GPU的嵌套網格裝配方法

摘" 要：采用嵌套網格可以有效地處理大幅運動問題，但隨著網格規模的增大和流動問題復雜度的提高，傳統的基于CPU的嵌套網格裝配方法越來越難以滿足當前的計算需求。針對上述問題，該文基于CUDA平臺，發展一種基于GPU的k-d樹嵌套網格裝配方法，并對k-d樹構建過程和搜索過程進行優化，大大提升貢獻單元搜索效率和物面距計算效率，進而加快嵌套網格裝配速度。

關鍵詞：圖形處理器；嵌套網格；k-d樹；裝配方法；流場計算域

中圖分類號：V211.3" " " 文獻標志碼：A" " " " " 文章編號：2095-2945（2025）01-0177-04

Abstract： Using overset grid can effectively handle problems with large-scale motions. However， as the scale of the grid and the complexity of the flow problem increase， traditional CPU-based overset grid assembly methods are becoming increasingly difficult to meet current computational demands. To address these issues， this paper develops a GPU-based k-d tree nested grid assembly method based on the CUDA platform， and optimizes the k-d tree construction process and search process， which greatly improves the contribution unit search efficiency and object-surface distance calculation efficiency， thereby accelerating the nested grid assembly speed.

Keywords： graphics processor; nested grid; k-d tree; assembly method; flow field computing domain

嵌套網格方法由Steger[1]提出，旨在簡化復雜幾何外形網格的生成。該方法將流場計算域劃分為多個嵌套區域，各區域獨立生成網格并求解，通過網格間插值傳遞流場信息。嵌套網格降低了網格生成難度，提高了靈活性，并保證網格質量，廣泛應用于多體分離數值模擬問題[2-6]。

與傳統單套網格不同，嵌套網格需進行裝配以確定嵌套邊界，對每套網格的單元或節點進行分類（亦即挖洞）。嵌套網格挖洞是嵌套網格裝配過程最主要也是技術難度較高的任務環節。在基于物面距的嵌套網格裝配方法中，物面距計算和宿主單元搜索占據了整個嵌套網格裝配中相當大的一部分時間，成為了嵌套網格裝配流程中極為關鍵的2個步驟，本文將介紹在GPU上實現這2個關鍵步驟的實施細節，關于嵌套網格裝配技術的其他實現細節參見文獻[7]。k-d樹[8]作為一種高查詢效率的二叉樹結構，其優異的平衡性可以顯著降低搜索深度，同時其構建主要基于排序，而GPU在執行排序任務上表現出色，因此本文采用k-d樹來進行物面距計算和宿主單元搜索。

1" 嵌套網格裝配流程

本文中采取的動態嵌套網格方法總體裝配流程為：首先將流場劃分為多個重疊的計算域并各自生成獨立的網格，按照重疊關系進行分層管理，然后在相鄰兩層或同層的各個子網格間，根據壁面距離的大小定義網格間的邊界，隨后對邊界進行拓寬和優化，以保證流場計算的高階精度格式；建立重疊區網格邊界的插值關系，用于流場計算時各區域間的流場信息交換；最后計算得到當前時間步的流場，時間推進一步，如有邊界運動，將貼近該邊界的網格移動到下一時間步的新位置，重新根據壁面距離定義網格邊界。嵌套網格效果如圖1所示，左邊的主翼面為背景網格，前緣襟翼為運動網格，灰色部分表示插值網格單元。

一般而言，網格疏密分布在貼近壁面附近較密，而遠離壁面時逐漸變稀，在本文中用物面距作為判定單元是否激活的依據。

確定嵌套網格邊界的具體算法如下。

循環每個網格簇，計算每個網格點到自身網格所包含的壁面的物面距。對不包含壁面的背景網格，從有壁面網格嵌入的網格層開始，層號從大到小，各網格層的物面距依次設定為Δs，2Δs，3Δs……，Δs的大小由用戶自定義，一般為0.5～2個特征長度。

所有網格點初始化為活動狀態。搜索網格點的宿主單元，比較網格點和宿主單元的物面距，物面距比宿主單元大的網格點，將其改為非激活狀態。

循環各個網格簇的網格單元，根據網格單元中的節點狀態對單元分類。如果單元的所有節點都為激活狀態則該單元為激活單元，如果都為非激活狀態則該單元為非激活單元，既有激活節點又有非激活節點的為邊界插值單元。

該過程不但確定了嵌套網格的邊界，而且還確定了重疊區內網格點的宿主單元，為建立插值關系提供了方便。

2" 物面距計算

物面距是網格單元距離物面邊界的最小距離，是嵌套網格隱式挖洞過程中對網格單元進行分類的主要判據。物面距的計算分為2個階段，物面距k-d樹構建和物面距k-d樹搜索。

2.1" 物面距k-d樹構建

在構建包含n個k維元素的k-d樹時，為了優化性能，避免在每個遞歸細分過程中對當前子樹元素重新排序以查找中值，通常會在構建k-d樹之前，預先對每個維度上的元素進行排序。預排序完成后，在子樹分支上始終保持該排序順序，根據已排序的數組來快速選擇中值點，并遞歸地構建左子樹和右子樹，在最糟糕情況下該算法的復雜度為O（kn log n）[9]。例如對于圖2所示的2維（x，y）元組，按照存放在元組中的位置其元素索引記為indices，依次通過復合鍵xy和yx對元素進行排序并保存排序后的元組索引，對鍵xy排序時依次按鍵x，y進行排序（鍵yx與之類似），此多鍵值排序可以通過按鍵重要性從低到高調用二次基本的穩定排序實現。預排序完成后，以xy為鍵排序的元組中間值為索引5代表的元素（7，2），然后以該元素對元組進行左右子樹劃分（分區），此處為了表述方便劃分方向按照x、y輪流劃分，實際實現中為了優化k-d樹的平衡性，劃分方向判據為使得劃分后元素集合方差最大。第一次以x軸分區后xy索引數組不需要分區，但是需要重新對yx索引數組重新分區。對于yx索引數組，依次取出yx索引數組中對應元素與中值元素（7，2）在鍵xy上進行比較，按照比較結果放入左右子樹分支，此過程需要創建新的yx索引數組。可以看到yx索引數組重新分區后左右分支部分元素間的前后關系與分區前相同，這種分區操作保留了左右分支中初始預排序順序，避免了重復排序。當節點包含元素個數為1、2或3時，分區終止。圖3為圖2中元組對應的k-d樹。

在上述k-d樹構建過程中關鍵的步驟是分區合并過程，為了在GPU上實現比較好的并行效率，需要采取分治的策略并充分利用GPU硬件特性優化合并過程，具體為核函數執行時每個線程束（wrap）讀取n個連續的索引數組元素，并與劃分元素進行比較，根據結果將索引值寫入到2個臨時輸出緩沖區中（分別保存左右分支），該部分調用線程束洗牌函數完成束內元素分區合并工作，然后通過共享內存和線程塊內同步函數__syncthread（）完成線程塊（block）內元素分區合并工作，最后完成所有線程塊內元素分區合并工作。

在上述預排序過程中并沒有對排序方法進行限定，在GPU上不同排序方法的效率存在較大差異，表1為對16 777 216個三維整型元素構建k-d樹的時間開銷，可以看到預排序過程在構建k-d樹的整體流程中占據了相當大的比重，同時基數排序明顯比歸并排序效率更高，在GPU上采用基數排序可大幅提升k-d樹的構建速度。此外如果元組元素在每個維度上的值均無重復則在預排序中無需按復合鍵進行排序，這可以通過對坐標進行一定角度的旋轉方式實現，此種優化策略也可大幅減少預排序時間加快構建速度。

2.2" 物面距k-d樹構建

相比于k-d樹構建過程，搜索過程易于并行，在k-d樹中搜索距離待查點最近鄰點的流程如下。

步驟1：根節點作為當前節點，初始化最小距離dmin（賦極大值），創建棧（保存最近訪問過的節點）。

步驟2：在當前劃分維上對待查點和當前節點進行比較，大于進入右子樹，小于進入左子樹，計算待查點和當前節點的距離dcur，并與最小距離dmin比較，若dcurlt;dmin，更新dmin并將當前節點更新為最近鄰節點Nnearest，將當前節點壓入棧中。

步驟3：對訪問的子樹重復步驟2。

步驟4：當訪問到葉子節點時，進行回溯?；厮輹r，進行出棧操作將彈出節點作為當前節點，判斷以待查點為球心，dmin為半徑的球面是否與當前劃分平面相交，如果相交則進入當前節點的另一個分支重復步驟2。至此得到dmin和Nnearest。

上述搜索方法在GPU中的實施當中主要有兩方面的困難，一是棧的空間開銷很大，二是大量的回溯操作會極大地降低搜索效率。如果在搜索開始時有一個初始的搜索半徑范圍即dmin范圍，則在步驟2中就可以判斷是否要進入左右子樹繼續搜索，可以很大程度減少進棧的節點數量。搜索半徑范圍的估計采用如下方法，物面距搜索采取陣面推進方式從與物面相鄰的網格單元開始一層一層向外搜索，下一層待搜索網格單元C2為上一層搜索網格單元C1的空間拓撲鄰居單元，網格單元Cx在物面上的最近鄰點記為Nearest（Cx），在搜索過程中保存最近鄰點編號，待搜索網格單元C2的初始搜索半徑范圍dmin滿足

dmin=d（C2，Nearest（C2））≤d（C2，Nearest（C1））。

此外，在回溯操作中以球面與劃分平面是否相交作為進一步搜索的依據并不總是合理的，例如在圖4所示的情況下，根據上述判別方法，則需要在圖中灰色區域內繼續搜索，而實際上該圓并未與灰色區域相交，此后的回溯操作均為無效過程。因此，考慮在k-d樹節點上維護一個額外的包圍盒數據結構，保存當前子樹中所有節點在每一維坐標上的最小和最大值，如果查詢點到該包圍盒的最短距離dbox≥dmin，則搜索時不進入該子樹。改進后可以極大地剪枝掉不必要的搜索子樹，節省搜索時間。

3" 宿主單元搜索

與物面距計算類似，在宿主單元搜索中也分為k-d樹的構建和k-d樹的搜索2個階段。

對網格單元構建k-d樹可以分為2種方法，一種是基于網格單元格心坐標直接構建k-d樹，然后采用最近鄰搜索找出距離查詢點最近的網格單元，該查詢點的貢獻單元必然在空間位置上靠近找到的最近鄰網格單元，再結合網格拓撲關系找出貢獻單元即可；另一種是基于網格單元的包圍盒構建k-d樹，搜索到的為可能的貢獻單元集合，再結合精確找重方法得到正確的貢獻單元，此方法相比于上一方法具有更好的可靠性，本文中采用該方法。二維情況下，網格單元i的包圍盒可以描述為超維空間中的超維坐標（ci，di），其中ci=（xmin，ymin），di=（xmax，ymax），分別表示在x、y維度上包圍網格單元的最小坐標值和最大坐標值。此時，可以將任一網格單元描述為位于超維空間區域內，這時樹節點需要8個浮點數據存儲邊界。因為在超維坐標中隱含有條件 ，網格單元對應的包圍盒其超維坐標并不在超維空間區域R內的一些區域內分布，為了構建平衡性較好的k-d樹，在構建階段分區時劃分維的選擇上需要額外判斷，如果當前的劃分維會導致某個子樹對應的超維空間區域根本不會有元素分布，則需要調整劃分維。在GPU上構建k-d樹的其他細節與物面距計算中類似，此處不再贅述。

在宿主單元搜索過程中采取深度優先遍歷的方式，流程如下。

步驟1：根節點作為當前節點，創建數組（保存可能的貢獻單元）。

步驟2：檢查待查點對應包圍盒超維坐標是否與當前節點保存的超維空間區域相交，如果相交，對左右子樹分別遞歸執行上述求相交操作；如果不相交，結束當前節點搜索。

步驟3：對數組中所有可能的貢獻單元進行精確找重（幾何包含）操作，找出真正與待查點元素相交的網格單元。

4" 結束語

本文介紹了一種采用k-d樹進行物面距計算和宿主單元搜索的方法，應用于基于GPU的嵌套網格裝配流程中，以提高嵌套網格裝配的效率和自動化程度。

參考文獻：

[1] STEGER J L， DOUGHERTY F C， BENEK J A. A chimera grid scheme[J]. NASA Technical Reports， 1983.

[2] MEAKIN R. Computations of the unsteady flow about a generic wing/pylon/finned-store configuration[J]. strodynamics Conference， 1992：4568.

[3] AHMAD J， SHANKS S， BUNING P. Aerodynamics of powered missile separation from F/A-18 aircraft[J]. 1st Aerospace Sciences Meeting， 1993：766.

[4] RIZK M， ELLISON S， PREWITT N. Beggar-A store separation predictive tool[C]// 32nd AIAA Fluid Dynamics Conference amp; Exhibit， 2002：3190.

[5] 張培紅，王明，鄧有奇，等.武器分離及艙門開啟過程數值模擬研究[J].空氣動力學學報，2013，31（3）：277-281.

[6] 招啟軍，徐國華.使用高階逆風通量差分裂格式的懸停旋翼流場數值模擬[J].航空動力學報，2005，20（2）：186-191.

[7] 肖天航，支豪林，朱震浩.飛行器非定常氣動計算與優化技術[M].北京：科學出版社，2022：100-130.

[8] BENTLEY J L. Multidimensional binary search trees in database applications[J]. IEEE Transactions on Software Engineering， 1979（4）：333-340.

[9] BROWN R A. Building k-d Tree in O（knlog n） Time[J]. Journal of Computer Graphics Techniques， 2015，4（1）.