基于GPU的嵌套網(wǎng)格裝配方法

摘" 要：采用嵌套網(wǎng)格可以有效地處理大幅運(yùn)動(dòng)問(wèn)題，但隨著網(wǎng)格規(guī)模的增大和流動(dòng)問(wèn)題復(fù)雜度的提高，傳統(tǒng)的基于CPU的嵌套網(wǎng)格裝配方法越來(lái)越難以滿足當(dāng)前的計(jì)算需求。針對(duì)上述問(wèn)題，該文基于CUDA平臺(tái)，發(fā)展一種基于GPU的k-d樹嵌套網(wǎng)格裝配方法，并對(duì)k-d樹構(gòu)建過(guò)程和搜索過(guò)程進(jìn)行優(yōu)化，大大提升貢獻(xiàn)單元搜索效率和物面距計(jì)算效率，進(jìn)而加快嵌套網(wǎng)格裝配速度。

關(guān)鍵詞：圖形處理器；嵌套網(wǎng)格；k-d樹；裝配方法；流場(chǎng)計(jì)算域

中圖分類號(hào)：V211.3" " " 文獻(xiàn)標(biāo)志碼：A" " " " " 文章編號(hào)：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

嵌套網(wǎng)格方法由Steger[1]提出，旨在簡(jiǎn)化復(fù)雜幾何外形網(wǎng)格的生成。該方法將流場(chǎng)計(jì)算域劃分為多個(gè)嵌套區(qū)域，各區(qū)域獨(dú)立生成網(wǎng)格并求解，通過(guò)網(wǎng)格間插值傳遞流場(chǎng)信息。嵌套網(wǎng)格降低了網(wǎng)格生成難度，提高了靈活性，并保證網(wǎng)格質(zhì)量，廣泛應(yīng)用于多體分離數(shù)值模擬問(wèn)題[2-6]。

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

1" 嵌套網(wǎng)格裝配流程

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

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

確定嵌套網(wǎng)格邊界的具體算法如下。

循環(huán)每個(gè)網(wǎng)格簇，計(jì)算每個(gè)網(wǎng)格點(diǎn)到自身網(wǎng)格所包含的壁面的物面距。對(duì)不包含壁面的背景網(wǎng)格，從有壁面網(wǎng)格嵌入的網(wǎng)格層開始，層號(hào)從大到小，各網(wǎng)格層的物面距依次設(shè)定為Δs，2Δs，3Δs……，Δs的大小由用戶自定義，一般為0.5～2個(gè)特征長(zhǎng)度。

所有網(wǎng)格點(diǎn)初始化為活動(dòng)狀態(tài)。搜索網(wǎng)格點(diǎn)的宿主單元，比較網(wǎng)格點(diǎn)和宿主單元的物面距，物面距比宿主單元大的網(wǎng)格點(diǎn)，將其改為非激活狀態(tài)。

循環(huán)各個(gè)網(wǎng)格簇的網(wǎng)格單元，根據(jù)網(wǎng)格單元中的節(jié)點(diǎn)狀態(tài)對(duì)單元分類。如果單元的所有節(jié)點(diǎn)都為激活狀態(tài)則該單元為激活單元，如果都為非激活狀態(tài)則該單元為非激活單元，既有激活節(jié)點(diǎn)又有非激活節(jié)點(diǎn)的為邊界插值單元。

該過(guò)程不但確定了嵌套網(wǎng)格的邊界，而且還確定了重疊區(qū)內(nèi)網(wǎng)格點(diǎn)的宿主單元，為建立插值關(guān)系提供了方便。

2" 物面距計(jì)算

物面距是網(wǎng)格單元距離物面邊界的最小距離，是嵌套網(wǎng)格隱式挖洞過(guò)程中對(duì)網(wǎng)格單元進(jìn)行分類的主要判據(jù)。物面距的計(jì)算分為2個(gè)階段，物面距k-d樹構(gòu)建和物面距k-d樹搜索。

2.1" 物面距k-d樹構(gòu)建

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

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

在上述預(yù)排序過(guò)程中并沒(méi)有對(duì)排序方法進(jìn)行限定，在GPU上不同排序方法的效率存在較大差異，表1為對(duì)16 777 216個(gè)三維整型元素構(gòu)建k-d樹的時(shí)間開銷，可以看到預(yù)排序過(guò)程在構(gòu)建k-d樹的整體流程中占據(jù)了相當(dāng)大的比重，同時(shí)基數(shù)排序明顯比歸并排序效率更高，在GPU上采用基數(shù)排序可大幅提升k-d樹的構(gòu)建速度。此外如果元組元素在每個(gè)維度上的值均無(wú)重復(fù)則在預(yù)排序中無(wú)需按復(fù)合鍵進(jìn)行排序，這可以通過(guò)對(duì)坐標(biāo)進(jìn)行一定角度的旋轉(zhuǎn)方式實(shí)現(xiàn)，此種優(yōu)化策略也可大幅減少預(yù)排序時(shí)間加快構(gòu)建速度。

2.2" 物面距k-d樹構(gòu)建

相比于k-d樹構(gòu)建過(guò)程，搜索過(guò)程易于并行，在k-d樹中搜索距離待查點(diǎn)最近鄰點(diǎn)的流程如下。

步驟1：根節(jié)點(diǎn)作為當(dāng)前節(jié)點(diǎn)，初始化最小距離dmin（賦極大值），創(chuàng)建棧（保存最近訪問(wèn)過(guò)的節(jié)點(diǎn)）。

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

步驟3：對(duì)訪問(wèn)的子樹重復(fù)步驟2。

步驟4：當(dāng)訪問(wèn)到葉子節(jié)點(diǎn)時(shí)，進(jìn)行回溯?；厮輹r(shí)，進(jìn)行出棧操作將彈出節(jié)點(diǎn)作為當(dāng)前節(jié)點(diǎn)，判斷以待查點(diǎn)為球心，dmin為半徑的球面是否與當(dāng)前劃分平面相交，如果相交則進(jìn)入當(dāng)前節(jié)點(diǎn)的另一個(gè)分支重復(fù)步驟2。至此得到dmin和Nnearest。

上述搜索方法在GPU中的實(shí)施當(dāng)中主要有兩方面的困難，一是棧的空間開銷很大，二是大量的回溯操作會(huì)極大地降低搜索效率。如果在搜索開始時(shí)有一個(gè)初始的搜索半徑范圍即dmin范圍，則在步驟2中就可以判斷是否要進(jìn)入左右子樹繼續(xù)搜索，可以很大程度減少進(jìn)棧的節(jié)點(diǎn)數(shù)量。搜索半徑范圍的估計(jì)采用如下方法，物面距搜索采取陣面推進(jìn)方式從與物面相鄰的網(wǎng)格單元開始一層一層向外搜索，下一層待搜索網(wǎng)格單元C2為上一層搜索網(wǎng)格單元C1的空間拓?fù)溧従訂卧W(wǎng)格單元Cx在物面上的最近鄰點(diǎn)記為Nearest（Cx），在搜索過(guò)程中保存最近鄰點(diǎn)編號(hào)，待搜索網(wǎng)格單元C2的初始搜索半徑范圍dmin滿足

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

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

3" 宿主單元搜索

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

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

在宿主單元搜索過(guò)程中采取深度優(yōu)先遍歷的方式，流程如下。

步驟1：根節(jié)點(diǎn)作為當(dāng)前節(jié)點(diǎn)，創(chuàng)建數(shù)組（保存可能的貢獻(xiàn)單元）。

步驟2：檢查待查點(diǎn)對(duì)應(yīng)包圍盒超維坐標(biāo)是否與當(dāng)前節(jié)點(diǎn)保存的超維空間區(qū)域相交，如果相交，對(duì)左右子樹分別遞歸執(zhí)行上述求相交操作；如果不相交，結(jié)束當(dāng)前節(jié)點(diǎn)搜索。

步驟3：對(duì)數(shù)組中所有可能的貢獻(xiàn)單元進(jìn)行精確找重（幾何包含）操作，找出真正與待查點(diǎn)元素相交的網(wǎng)格單元。

4" 結(jié)束語(yǔ)

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

參考文獻(xiàn)：

[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] 張培紅，王明，鄧有奇，等.武器分離及艙門開啟過(guò)程數(shù)值模擬研究[J].空氣動(dòng)力學(xué)學(xué)報(bào)，2013，31（3）：277-281.

[6] 招啟軍，徐國(guó)華.使用高階逆風(fēng)通量差分裂格式的懸停旋翼流場(chǎng)數(shù)值模擬[J].航空動(dòng)力學(xué)報(bào)，2005，20（2）：186-191.

[7] 肖天航，支豪林，朱震浩.飛行器非定常氣動(dòng)計(jì)算與優(yōu)化技術(shù)[M].北京：科學(xué)出版社，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）.