胡曉東,梁 姍,袁 武,陸忠華

1(中國科學(xué)院 計算機網(wǎng)絡(luò)信息中心,北京 100190)

2(中國科學(xué)院大學(xué),北京 100049)

計算流體力學(xué)方法將連續(xù)的物理空間離散為網(wǎng)格單元構(gòu)成的計算空間,通過對流體控制方程(Navier-Stokes方程)的離散求解達到揭示分析流動特性的目的.計算網(wǎng)格從單元形狀和數(shù)據(jù)結(jié)構(gòu)上可分為結(jié)構(gòu)化網(wǎng)格、非結(jié)構(gòu)化網(wǎng)格和笛卡爾網(wǎng)格等[1].結(jié)構(gòu)網(wǎng)格的特點是拓撲結(jié)構(gòu)有序,邏輯關(guān)系簡單,流場計算精度和效率高,缺點是對復(fù)雜外型的適應(yīng)性差,在確定各種復(fù)雜外型的空間拓撲關(guān)系時顯得非常困難,難以自動化生成,從而需要投入大量人工時間.Steger[2,3]提出了“重疊網(wǎng)格”的概念,將流場模擬區(qū)域分解為若干子區(qū)域,在每個子區(qū)域上分別生成高質(zhì)量網(wǎng)格,各個子區(qū)域之間存在相互重疊關(guān)系.在流場計算過程中,通過重疊網(wǎng)格內(nèi)各套子網(wǎng)格之間的插值實現(xiàn)流場信息傳遞,這樣對于復(fù)雜的外型結(jié)構(gòu),不同的外型區(qū)域可以單獨生成拓撲,不用再進行耗時的復(fù)雜拓撲劃分,彌補了結(jié)構(gòu)網(wǎng)格的缺點,同時保持了結(jié)構(gòu)化網(wǎng)格計算精度高和邊界處理能力強的優(yōu)點,擴展了對復(fù)雜外型和運動問題的適應(yīng)性,減少了生成網(wǎng)格的時間和工作量,因而得到了較廣泛的應(yīng)用.

重疊網(wǎng)格的實現(xiàn)過程分為挖洞、找點和插值3個步驟[4],即首先確定并標(biāo)記出流場計算中無意義的重疊網(wǎng)格單元、然后確定各個網(wǎng)格點在其他網(wǎng)格中的位置對應(yīng)關(guān)系和相應(yīng)的插值貢獻單元,最后在重疊邊界處通過貢獻單元向插值點傳遞流場信息,實現(xiàn)整個流場的同步更新.為了減少重疊網(wǎng)格處理的迭代時間,通常采用陣面推進類方法確定洞邊界和進行洞面優(yōu)化[5],此類方法對網(wǎng)格點數(shù)據(jù)的連續(xù)性要求高,不利于并行擴展和應(yīng)用,尤其對于網(wǎng)格量大的復(fù)雜外型流場模擬.YikLoon Lee[6,7]提出了結(jié)構(gòu)網(wǎng)格的隱式挖洞算法,相比于陣面推進類方法,隱式挖洞算法采用網(wǎng)格密度嚴格對比的找點過程,替換傳統(tǒng)重疊網(wǎng)格方法的網(wǎng)格狀態(tài)標(biāo)識和貢獻單元尋址步驟,保留流場信息插值這一環(huán)節(jié),其最大特點是不需要額外的重疊關(guān)系輸入文件,重疊網(wǎng)格切割和貢獻單元尋址在一個過程中同時完成.Landmann B[8]對隱式挖洞的流程進行了并行化實現(xiàn),分析了網(wǎng)格裝配各個階段的耗時對比以及算法并行擴展性.Kenway GK[9]等對隱式重疊方法的并行特性進行了較為細致的研究,根據(jù)不同網(wǎng)格塊之間的重疊關(guān)系,建立“重疊矩陣”,將網(wǎng)格比較分配到不同進程進行,并根據(jù)重疊網(wǎng)格區(qū)域大小、實際的查找耗時對網(wǎng)格的分配進行優(yōu)化,達到數(shù)千核的并行擴展效果,并指出了隱式挖洞方法負載不均衡的局限性.

隱式挖洞算法中網(wǎng)格之間的比較判斷會占用大量時間開銷,對于工程實際的復(fù)雜構(gòu)型,網(wǎng)格重疊的結(jié)果可能是某些區(qū)域的子網(wǎng)格塊需要和所重疊的其他子網(wǎng)格中大量網(wǎng)格塊進行信息交換.因此在并行算法設(shè)計中,需要綜合考慮網(wǎng)格挖洞插值和計算階段的負載平衡,并盡量減少全局通信.本文在并行隱式挖洞算法實現(xiàn)的基礎(chǔ)上,提出了笛卡爾輔助網(wǎng)格和多塊結(jié)構(gòu)網(wǎng)格的混合網(wǎng)格方法.通過基于搭接邊界的多層背景網(wǎng)格,實現(xiàn)隱式挖洞的快速確定洞邊界和網(wǎng)格間插值關(guān)系的快速建立；通過對重疊區(qū)域網(wǎng)格重疊權(quán)重的定義和應(yīng)用,建立混合網(wǎng)格的并行分配模式和負載均衡模型,有效減少重疊插值信息在各進程間的通信,并實現(xiàn)計算負載和通信負載在各個進程的均勻分配.

1 并行混合重疊網(wǎng)格算法

1.1 并行隱式挖洞算法

并行隱式挖洞算法的實現(xiàn)流程如圖1所示,將讀入的多塊結(jié)構(gòu)網(wǎng)格分為背景網(wǎng)格和子區(qū)域網(wǎng)格2類：背景網(wǎng)格包圍主要部件,覆蓋從部件物面邊界到遠場的整個流場區(qū)域,子區(qū)域網(wǎng)格包圍其他部件,從部件的物面邊界法向向外延伸一定距離,重疊在背景網(wǎng)格的區(qū)域之內(nèi).將背景網(wǎng)格和子區(qū)域網(wǎng)格塊分配到不同的節(jié)點.在各個節(jié)點內(nèi),根據(jù)網(wǎng)格質(zhì)量的定義,分別計算所分配網(wǎng)格點的網(wǎng)格單元體積和到壁面的最短距離,作為網(wǎng)格質(zhì)量的組成因素.然后對每個網(wǎng)格單元進行初始化標(biāo)記,將所有的背景網(wǎng)格標(biāo)記為計算,將所有的子區(qū)域網(wǎng)格標(biāo)記為插值.

圖1 并行隱式挖洞方法流程

之后在每個進程內(nèi)計算各個子區(qū)域的坐標(biāo)范圍,歸約到主進程,得到各個子區(qū)域的包圍盒,再將所有包圍盒信息發(fā)送到各個節(jié)點.后續(xù)僅需要考慮包圍盒內(nèi)的網(wǎng)格單元.對于這些單元需要進行網(wǎng)格標(biāo)記并尋找貢獻單元,在各個節(jié)點進行的網(wǎng)格質(zhì)量比較,如果子區(qū)域網(wǎng)格單元質(zhì)量更高,將其標(biāo)記由改為計算,所對應(yīng)的背景網(wǎng)格單元標(biāo)記為插值,然后將更新后的網(wǎng)格屬性標(biāo)記收集并存儲到主進程,主進程得到所有更新后的網(wǎng)格信息.用于下一階段的流場計算.

網(wǎng)格密度的定義需要結(jié)合網(wǎng)格生成的特點,包括網(wǎng)格單元的體積和網(wǎng)格距離物面邊界的距離,通常在物面附近,網(wǎng)格具有較小的網(wǎng)格單元體積,較好的正交性和合理的過度.本文中采用網(wǎng)格體積和網(wǎng)格到壁面距離的組合作為網(wǎng)格密度的定義,如式(1)所示：

其中,v為網(wǎng)格單元的體積,d為網(wǎng)格距離物面邊界的距離,通過a和b參數(shù)的調(diào)節(jié)改變網(wǎng)格體積和網(wǎng)格距離物面距離兩個因素的影響,取a=1,b=0時僅考慮網(wǎng)格尺度作為網(wǎng)格質(zhì)量,取a=0,b=1時僅考慮網(wǎng)格到物面距離作為網(wǎng)格質(zhì)量.對網(wǎng)格體積的計算可以在各個網(wǎng)格單元所分配的進程內(nèi)分別進行,不需要與其他進程進行數(shù)據(jù)交換.

1.2 笛卡爾背景網(wǎng)格

本文將笛卡爾輔助背景網(wǎng)格引入到重疊網(wǎng)格系統(tǒng),在多塊結(jié)構(gòu)網(wǎng)格的重疊網(wǎng)格中,通常對模擬中主要部件或非定常過程中靜止部件生成覆蓋整個模擬區(qū)域的多塊網(wǎng)格,如多段翼型或機翼的主翼、飛行器掛載分離模擬中的飛行器、子母彈拋撒模擬中的母彈等；對于運動部件或復(fù)雜的構(gòu)件,生成覆蓋流場局部區(qū)域的子網(wǎng)格塊,構(gòu)成重疊網(wǎng)格系統(tǒng).采用笛卡爾輔助網(wǎng)格后,在物體物面到法向向外足夠包含附面層區(qū)域的一定距離內(nèi),各自生成貼體的局部網(wǎng)格作為輸入文件,流場中的其他范圍由在預(yù)處理中生成的笛卡爾網(wǎng)格覆蓋.圖2顯示的是加入笛卡爾輔助網(wǎng)格后30P30N多段翼型模擬需要的網(wǎng)格輸入,即各自圍繞主翼、襟翼和副翼,在法向一定距離內(nèi)生成的貼體曲面網(wǎng)格.在預(yù)處理過程中生成覆蓋從各個物面到遠場模擬區(qū)域的笛卡爾網(wǎng)格作為輔助背景網(wǎng)格.在靠近模擬物體附近的輔助笛卡爾網(wǎng)格中,各個方向采用固定的網(wǎng)格間隔,這樣在背景網(wǎng)格某一區(qū)域內(nèi)可以根據(jù)網(wǎng)格點的編號方便的確定出網(wǎng)格的體積,即隱式挖洞中的網(wǎng)格密度(背景輔助網(wǎng)格中不存在物面邊界,不必考慮網(wǎng)格點到壁面的距離).

在預(yù)處理階段中生成輔助笛卡爾網(wǎng)格需要確定輔助網(wǎng)格單元的尺寸和整個網(wǎng)格區(qū)域的大小,兩者均可以從輸入的貼體網(wǎng)格中得到.為了使重疊邊界產(chǎn)生在遠離固壁邊界的的合適位置,首先計算出子區(qū)域網(wǎng)單元尺寸的最小值和最大值,再根據(jù)具體問題乘以合適的比例系數(shù).得到背景網(wǎng)格的單元尺寸.通過計算各個模擬物體在各方向的長度,在各個方向上擴大一定距離確定遠場邊界的位置.生成的笛卡爾網(wǎng)格按規(guī)則結(jié)構(gòu)網(wǎng)格的數(shù)據(jù)形式組織,與各個部件網(wǎng)格在預(yù)處理階段統(tǒng)一進行網(wǎng)格的分區(qū),并分配至各個進程,圖3展示的是執(zhí)行隱式挖洞后的效果,各段翼型與背景的重疊邊界,各段翼型之間的重疊邊界均清晰的建立起來.

圖2 翼型附近區(qū)域的貼體網(wǎng)格

圖3 挖洞處理后的翼型網(wǎng)格和笛卡爾網(wǎng)格邊界

1.3 基于笛卡爾網(wǎng)格域的找點插值方法

重疊網(wǎng)格找點的過程即是建立網(wǎng)格間插值關(guān)系的過程,這一過程往往是重疊網(wǎng)格關(guān)系中耗時最多的步驟.引入笛卡爾輔助背景網(wǎng)格后,可以利用輔助網(wǎng)格數(shù)據(jù)結(jié)構(gòu)簡單,網(wǎng)格單元標(biāo)記和網(wǎng)格點坐標(biāo)存在直接聯(lián)系的特點,通過背景網(wǎng)格尋找各個網(wǎng)格單元所要進行比較的單元,也同時標(biāo)記為插值的貢獻單元,在不同重疊子網(wǎng)格單元之間建立插值關(guān)系.

對近壁面區(qū)域背景網(wǎng)格采用相同的網(wǎng)格單元,各方向網(wǎng)格單元的間隔是相等的,即所有的網(wǎng)格單元具有相同的Δx,Δy和Δz,因此可以建立網(wǎng)格數(shù)組標(biāo)記和網(wǎng)格坐標(biāo)的直接聯(lián)系,根據(jù)部件網(wǎng)格坐標(biāo)確定背景網(wǎng)格數(shù)組值的步驟為：對某一個部件網(wǎng)格單元,首先根據(jù)輔助背景網(wǎng)格的分塊邊界,確定網(wǎng)格單元位于哪一個背景網(wǎng)格塊之內(nèi),之后根據(jù)該網(wǎng)格單元格心的坐標(biāo)與部件網(wǎng)格邊界坐標(biāo)各方向的差,以及背景各個方向的單元大小,確定出該部件網(wǎng)格單元在背景塊內(nèi)的編號.xo,yo,zo表示部件網(wǎng)格點(l,m,n)的坐標(biāo),xb,yb,zb表示背景網(wǎng)格的坐標(biāo),Δx,Δy和Δz表示在各個方向上背景網(wǎng)格的單元間隔,通過式(2),即可以得到與該部件網(wǎng)格存在插值關(guān)系的背景網(wǎng)格單元的i,j,k編號.

1.4 重疊權(quán)重定義

定義網(wǎng)格塊重疊權(quán)重系數(shù)的概念,將所有的網(wǎng)格被分為2個部分,對于不存在重疊關(guān)系的網(wǎng)格單元,定義為背景區(qū)域,其權(quán)重系數(shù)設(shè)置為1,對于存在重疊關(guān)系的網(wǎng)格單元,定義為混合區(qū)域,如果有n套網(wǎng)格,那么權(quán)重系數(shù)定義為n.在將網(wǎng)格分配到各個進程時,根據(jù)重疊關(guān)系將位于同一坐標(biāo)區(qū)域內(nèi)的網(wǎng)格盡可能分配到同一進程,達到避免部分進程間插值信息通信的目的.

2 通信優(yōu)化算法實現(xiàn)

2.1 網(wǎng)格塊與背景輔助網(wǎng)格塊的綁定

在確立網(wǎng)格插值關(guān)系的流程之后,每個部件塊的網(wǎng)格單元均記錄了相應(yīng)的背景塊中網(wǎng)格單元插值信息,對各個部件網(wǎng)格塊的信息進行遍歷,累加統(tǒng)計出每個部件網(wǎng)格塊在各個背景網(wǎng)格塊內(nèi)插值點的數(shù)目.比較得到與該部件網(wǎng)格塊插值點最多的背景網(wǎng)格塊,即與附件網(wǎng)格塊重疊區(qū)域最大的背景塊,遍歷結(jié)束后,每個部件網(wǎng)格塊綁定唯一的一個背景網(wǎng)格塊編號,在后續(xù)分配到各個計算節(jié)點時,將部件網(wǎng)格塊與其綁定的背景網(wǎng)格塊分配到同一個節(jié)點內(nèi),如圖4所示,藍色網(wǎng)格線表示部件貼體子網(wǎng)格的網(wǎng)格塊,紅色網(wǎng)格線表示背景輔助網(wǎng)格的網(wǎng)格塊,部件網(wǎng)格塊(1234)與背景網(wǎng)格多個網(wǎng)格塊重疊,根據(jù)重疊關(guān)系可分為(156)、(2567D9)、(738D)、(48D)和(94D)5個部分,根據(jù)網(wǎng)格塊(1234)中各個單元記錄的插值關(guān)系,計算出與其重疊的五個背景網(wǎng)格塊內(nèi)重疊的點數(shù),比較得到網(wǎng)格點數(shù)最多的背景網(wǎng)格塊,如圖所示為背景網(wǎng)格塊(ABCD),標(biāo)記為部件網(wǎng)格塊(1234)綁定的背景網(wǎng)格塊,在并行分配過程中分配到同一節(jié)點.

圖4 確定部件網(wǎng)格塊綁定的背景網(wǎng)格塊

在執(zhí)行該流程的同時,對于每一個背景網(wǎng)格塊,統(tǒng)計與其綁定的部件網(wǎng)格塊數(shù)目,作為該背景網(wǎng)格塊分配的重疊權(quán)重系數(shù),并將該背景網(wǎng)格塊標(biāo)記為重疊塊,沒有被標(biāo)記的背景網(wǎng)格塊為非重疊塊.在分配過程中,重疊背景塊與非重疊背景塊分開處理,對于重疊背景塊,其分配的節(jié)點預(yù)留其綁定部件塊的空間,將與其綁定的塊分配在同一個進程.

2.2 網(wǎng)格塊的分配算法

網(wǎng)格預(yù)處理分塊后的塊數(shù)達到數(shù)倍于計算分配的節(jié)點數(shù),實際分配中網(wǎng)格塊會循環(huán)多次.由于分塊分配在預(yù)處理階段進行,所有網(wǎng)格塊最終循環(huán)輪數(shù)是不確定的,每個節(jié)點上最多分配的網(wǎng)格塊數(shù)量也是無法預(yù)知的.所以采取如下處理方法：把重疊的背景塊最先處理,按照重疊權(quán)重降序進行排列,第一次循環(huán)每個網(wǎng)格塊都分配節(jié)點,從第二次循環(huán)開始,直到重疊系數(shù)為n的背景塊循環(huán)到第n+1次之前,跳過該背景網(wǎng)格塊不分配,以預(yù)留給綁定的部件網(wǎng)格塊.當(dāng)所有背景網(wǎng)格塊循環(huán)分配完畢后,把各個部件網(wǎng)格塊分配到所綁定的背景塊所在節(jié)點.

以圖5為例進行說明,重疊網(wǎng)格系統(tǒng)中共有背景網(wǎng)格塊25塊,部件網(wǎng)格塊3塊,分配到8個進程.背景網(wǎng)格中編號4、6、7的塊與部件網(wǎng)格重疊,第4塊與部件網(wǎng)格塊A和B重疊,權(quán)重系數(shù)定義為3,第6塊和第7塊分別與部件網(wǎng)格塊C和D重疊,權(quán)重系數(shù)定義為2.在分配之前把有重疊關(guān)系的背景塊按權(quán)重系數(shù)排在前面,確保重疊的背景塊先分配到進程號,從第2次循環(huán)開始,第4塊輪空兩次,第6塊和第7塊輪空一次,之后的非重疊塊依次分配.待所有背景網(wǎng)格塊分配完后,把部件網(wǎng)格塊分配到相應(yīng)的背景塊所在節(jié)點.

圖5 根據(jù)重疊權(quán)重進行網(wǎng)格分配

3 算法測試與分析

為了驗證本文提出混合重疊網(wǎng)格方法和通信優(yōu)化方法的有效性.選用ANF模型[10,11]對算法進行測試分析.ANF模型由錐頭圓柱彈身和4個尖頭彈翼組成,彈體子網(wǎng)格的網(wǎng)格量為1039萬,每個彈翼網(wǎng)格子網(wǎng)格的網(wǎng)格量為52萬,在預(yù)處理過程中根據(jù)貼體網(wǎng)格信息,生成多層背景網(wǎng)格,組成總計3017萬的重疊網(wǎng)格系統(tǒng).輸入ANF模型的貼體網(wǎng)格拓撲如圖6所示.

圖6 ANF貼體重疊網(wǎng)格拓撲結(jié)構(gòu)

測試在中國科學(xué)院計算機網(wǎng)絡(luò)信息中心的“元”高性能計算系統(tǒng)進行,系統(tǒng)每個刀片節(jié)點配備兩顆Intel Xeon E5-2680V3處理器(2.5 GHz,12核),256 GB內(nèi)存.編譯環(huán)境為Intel Fortran+OpenMPI,使用-O2編譯優(yōu)化選項.分別使用64核、128核、256核、512核及1024核進行測試并分析.

在不同規(guī)模的測試中,統(tǒng)計各個進程內(nèi)需要與其他進程內(nèi)需要交換插值信息的網(wǎng)格單元,作為通信量比較的依據(jù).圖7-圖10展示了不同規(guī)模測試優(yōu)化前后的通信量對比,實線表示采用循環(huán)分配法依此分配各個網(wǎng)格塊時,每個進程內(nèi)的重疊插值信息通信量；虛線表示采用優(yōu)化算法后,每個進程內(nèi)需要通信的插值信息.測試結(jié)果表明本文提出的混合重疊網(wǎng)格方法可擴展應(yīng)用至千核規(guī)模的重疊網(wǎng)格模擬,在不同規(guī)模下,各個進程的通信量基本平衡,通信優(yōu)化算法能夠有效減少進程間插值信息的傳遞.

圖7 128核測試通信量優(yōu)化效果

圖8 256核測試通信量優(yōu)化效果

圖9 512核測試通信量優(yōu)化效果

圖10 1024核測試通信量優(yōu)化效果

在并行測試過程中,將執(zhí)行混合重疊網(wǎng)格方法的預(yù)處理流程分解為7個部分,分別統(tǒng)計這7個部分的運行時間.這7個部分包括：

(1)多層背景網(wǎng)格建立Background time.

(2)計算各個網(wǎng)格單元的體積Volume time.

(3)計算各個單元的到最近物面距離Distance time.

(4)網(wǎng)格塊的剖分Split time.

(5)建立網(wǎng)格間插值關(guān)系和重疊權(quán)重系數(shù)Overset time.

(6)發(fā)送網(wǎng)格塊信息和重疊關(guān)系至各個進程Sendrecv time.

(7)隱式挖洞標(biāo)識IHC time.

統(tǒng)計的時間結(jié)果如表1所示,展示了從128到1024和時不同處理環(huán)節(jié)消耗的時間,從統(tǒng)計結(jié)果看,計算網(wǎng)格標(biāo)記屬性的體積和最短壁面距離所消耗的時間最少,建立插值關(guān)系是整個流程中最耗時的部分,隨著測試規(guī)模的擴大,網(wǎng)格塊規(guī)模的減少而耗時減少；網(wǎng)格信息發(fā)送至各進程的時間耗時較多并隨著測試規(guī)模擴大而減少；隱式挖洞標(biāo)識的時間基本不隨規(guī)模增大而變化.

表1 重疊網(wǎng)格不同環(huán)節(jié)時間對比(單位：s)

4 結(jié)論

本文采用多層嵌套的背景網(wǎng)格替代空間區(qū)域的結(jié)構(gòu)網(wǎng)格,建立了笛卡爾輔助網(wǎng)格+多塊結(jié)構(gòu)網(wǎng)格的混合網(wǎng)格系統(tǒng).通過輔助網(wǎng)格實現(xiàn)洞邊界的快速確定和網(wǎng)格間插值關(guān)系的快速建立,在此基礎(chǔ)上實現(xiàn)了混合網(wǎng)格的并行重疊流程.通過定義重疊區(qū)域網(wǎng)格權(quán)重、部件網(wǎng)格與背景網(wǎng)格綁定的方法,建立了混合網(wǎng)格的并行分配模式,有效減少重疊插值信息在各進程間的通信,并實現(xiàn)計算負載和通信負載在各個進程的均勻分配.測試表明該方法可應(yīng)用于數(shù)千萬量級的重疊網(wǎng)格系統(tǒng),可擴展至千核規(guī)模,高效的實現(xiàn)多個物體構(gòu)成的復(fù)雜網(wǎng)格系統(tǒng)的重疊關(guān)系建立.