劉宏康， 閻超， 林博希， 趙雅甜

北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院， 北京 100083

基于圖剖分的多塊結(jié)構(gòu)網(wǎng)格負(fù)載平衡方法

劉宏康， 閻超*， 林博希， 趙雅甜

北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院， 北京 100083

負(fù)載平衡是影響并行計(jì)算性能的重要因素。針對(duì)多塊結(jié)構(gòu)網(wǎng)格，給出了一種改進(jìn)的多層次圖剖分負(fù)載平衡方法。該方法設(shè)計(jì)了新的網(wǎng)格剖分算法，采用改進(jìn)的子塊分裂方法與圖剖分算法的循環(huán)調(diào)用實(shí)現(xiàn)結(jié)構(gòu)對(duì)接網(wǎng)格剖分，并通過(guò)建立不同物體重疊網(wǎng)格間的連接關(guān)系，實(shí)現(xiàn)了結(jié)構(gòu)重疊網(wǎng)格的負(fù)載平衡。采用2個(gè)典型算例對(duì)方法進(jìn)行了對(duì)比驗(yàn)證，數(shù)值結(jié)果表明，子塊分裂方法對(duì)剖分結(jié)果具有重要影響，采用循環(huán)調(diào)用算法及改進(jìn)的子塊分裂方法能有效地實(shí)現(xiàn)計(jì)算負(fù)載均衡及通信量?jī)?yōu)化，同時(shí)顯著減少了網(wǎng)格塊數(shù)及因虛網(wǎng)格導(dǎo)致的內(nèi)存需求，有利于提高并行效率。該負(fù)載平衡方法與網(wǎng)格拓?fù)錈o(wú)關(guān)，適用于多塊結(jié)構(gòu)對(duì)接網(wǎng)格及重疊網(wǎng)格，且整體型剖分方式對(duì)于多塊結(jié)構(gòu)重疊網(wǎng)格具有更好的剖分效果。

計(jì)算流體力學(xué)； 并行計(jì)算； 結(jié)構(gòu)網(wǎng)格； 負(fù)載平衡； 圖剖分

計(jì)算流體力學(xué)(Computational Fluid Dynamics，CFD)中通常需要將物理流場(chǎng)離散化為計(jì)算網(wǎng)格，其中，多塊結(jié)構(gòu)網(wǎng)格(單塊為結(jié)構(gòu)網(wǎng)格但各塊間的連接關(guān)系可以為一一對(duì)應(yīng)或重疊)技術(shù)扮演著重要角色，在航空航天、氣象預(yù)測(cè)等領(lǐng)域應(yīng)用廣泛。

隨著工程問(wèn)題日益復(fù)雜化，高精度的數(shù)值模擬需求對(duì)計(jì)算格式及網(wǎng)格分辨率提出了更高要求，進(jìn)而導(dǎo)致計(jì)算量急劇增加。受限于串行程序的計(jì)算效率，并行計(jì)算是解決上述問(wèn)題的重要手段，而影響并行性能的關(guān)鍵因素之一是計(jì)算負(fù)載平衡的實(shí)現(xiàn)。

當(dāng)前的CFD并行計(jì)算普遍采用“單程序多數(shù)據(jù)”模式，需要將計(jì)算網(wǎng)格映射到分布式存儲(chǔ)系統(tǒng)的各個(gè)進(jìn)程上。在復(fù)雜實(shí)際問(wèn)題中，計(jì)算區(qū)域一般由多個(gè)大小不等的子塊組成，且存在相互間流場(chǎng)信息傳遞的復(fù)雜問(wèn)題，同時(shí)伴隨著并行規(guī)模的增加，這給結(jié)構(gòu)網(wǎng)格的負(fù)載平衡帶來(lái)了困難。因此，發(fā)展適用于復(fù)雜實(shí)際問(wèn)題的結(jié)構(gòu)網(wǎng)格負(fù)載平衡方法，對(duì)實(shí)現(xiàn)計(jì)算負(fù)載均衡及進(jìn)程間通信成本優(yōu)化具有重要意義。

現(xiàn)階段國(guó)內(nèi)外對(duì)負(fù)載平衡方法的研究較為豐富，提出了一系列啟發(fā)性算法，如貪婪算法[1]、譜方法[2]、幾何剖分方法[3]及圖剖分方法[4]等，并已集成至通用軟件包，如文獻(xiàn)[5-6]等。上述算法對(duì)非結(jié)構(gòu)網(wǎng)格具有天然適用性，但無(wú)法直接適用于結(jié)構(gòu)網(wǎng)格，且在該方面的研究也相對(duì)有限。國(guó)外方面，Ytterstr?m[7]最早提出了遞歸邊對(duì)分(Recursive Edge Bisection，REB)方法進(jìn)行網(wǎng)格細(xì)分以實(shí)現(xiàn)負(fù)載平衡，而Rantakokko[8]進(jìn)一步發(fā)展了多層次剖分的負(fù)載平衡方法，在一定程度上減少了網(wǎng)格塊數(shù)。Ahusborde和Glockner[9]重點(diǎn)研究了非規(guī)則矩形流場(chǎng)的子塊分割算法并進(jìn)行了流場(chǎng)求解，Djomehri和Biswas[10]則采用REB方法進(jìn)行了重疊網(wǎng)格的子塊分割，并基于Overflow求解器開(kāi)展求解，對(duì)比不同算法并行效果。國(guó)內(nèi)方面，司海青和王同光[11]提出了一種無(wú)子塊分割的負(fù)載剖分算法，李桂波和楊國(guó)偉[12]則采用遺傳優(yōu)化算法研究了負(fù)載與通信的并行影響，許正等[13]采用負(fù)載再分配策略來(lái)實(shí)現(xiàn)動(dòng)態(tài)負(fù)載平衡，并提出近似均分子塊分割方法。上述研究往往以對(duì)接結(jié)構(gòu)網(wǎng)格為主，難以適用于多塊結(jié)構(gòu)重疊網(wǎng)格，并且由于對(duì)網(wǎng)格塊數(shù)的關(guān)注不足，存在子塊分割過(guò)多或剖分效果受網(wǎng)格拓?fù)?、進(jìn)程數(shù)影響較大的缺陷。

本文針對(duì)負(fù)載平衡方法進(jìn)行了研究，提出了一種改進(jìn)的結(jié)構(gòu)網(wǎng)格負(fù)載平衡方法。該方法采用網(wǎng)格子塊分割及圖剖分的循環(huán)調(diào)用思路，同時(shí)改進(jìn)了子塊分裂方法，并通過(guò)建立重疊邊界關(guān)聯(lián)，將其拓展至多塊結(jié)構(gòu)重疊網(wǎng)格。最后，對(duì)三段翼30P30N對(duì)接網(wǎng)格及翼身組合體DLR-F6 WB重疊網(wǎng)格算例進(jìn)行數(shù)值對(duì)比分析，驗(yàn)證本文方法的可行性。

1 控制方程及數(shù)值方法

基于格心有限體積法求解控制方程，當(dāng)不考慮外加熱和體力影響時(shí)，直角坐標(biāo)系下三維可壓縮Navier-Stokes方程組的守恒積分形式為

(1)

式中：Ω為控制體；S為控制面；n為微元S的單位法向量；Q為守恒變量；Fc為對(duì)流通量；Fv為黏性通量。

本文采用多塊結(jié)構(gòu)網(wǎng)格，將流場(chǎng)劃分為多個(gè)計(jì)算區(qū)塊，并在每個(gè)區(qū)塊建立局部計(jì)算坐標(biāo)系，而區(qū)塊之間的相鄰關(guān)系采用對(duì)接型或重疊型實(shí)現(xiàn)。

數(shù)值方法包括空間離散方法和時(shí)間推進(jìn)方法，其中，對(duì)流項(xiàng)離散采用Roe格式，通過(guò)MUSCL (Monotone Upstream-centred Schemes for Conservation Laws)插值實(shí)現(xiàn)二階精度，并采用Mixture限制器消除非物理振蕩，黏性項(xiàng)采用二階中心差分格式，時(shí)間推進(jìn)采用LU-SGS (Lower-Upper Symmetric Gauss-Seidel)隱式方法。湍流模型采用兩方程SST (Shear-Stress Transport)模型。

2 負(fù)載平衡問(wèn)題

負(fù)載平衡實(shí)際上是在并行機(jī)或分布式多處理機(jī)的各節(jié)點(diǎn)之間合理分配計(jì)算任務(wù)的過(guò)程，一般分為靜態(tài)負(fù)載平衡問(wèn)題和動(dòng)態(tài)負(fù)載平衡問(wèn)題。本文所研究的結(jié)構(gòu)網(wǎng)格負(fù)載平衡在開(kāi)始數(shù)值計(jì)算前進(jìn)行網(wǎng)格劃分，屬于靜態(tài)負(fù)載平衡問(wèn)題。

為滿足CFD并行效率要求，在把計(jì)算網(wǎng)格映射到各個(gè)進(jìn)程上時(shí)，應(yīng)使得各進(jìn)程的網(wǎng)格單元數(shù)盡可能接近，且滿足進(jìn)程間需要交換的信息量最小，這種映射一般可視為圖剖分問(wèn)題。

所謂圖剖分問(wèn)題[14]，就是把數(shù)據(jù)間的相互依賴關(guān)系抽象成一個(gè)圖，通常用加權(quán)無(wú)向圖來(lái)表示計(jì)算和通信關(guān)系。加權(quán)無(wú)向圖G={E,V,Wv,We}由頂點(diǎn)集V={v1,v2,…,vn}和邊集E={e1,e2,…,em}以及相應(yīng)頂點(diǎn)權(quán)重集Wv={w(v1),w(v2),…,w(vn)}和邊權(quán)重集We={w(e1),w(e2),…,w(em)}構(gòu)成，n為頂點(diǎn)數(shù)，m為邊數(shù)，w(vi)為頂點(diǎn)的權(quán)重，表示頂點(diǎn)的計(jì)算量大小，w(ei)為邊的權(quán)重，表示通信量大小。

圖1 圖的兩種剖分示意 Fig.1 Schematic diagram of two kinds of graph partition

3 負(fù)載平衡方法

現(xiàn)有的圖剖分算法(如多水平算法[15]、多尺度KL/FM算法[16]等)對(duì)于非結(jié)構(gòu)網(wǎng)格具有天然的適用性，但無(wú)法直接用于結(jié)構(gòu)網(wǎng)格的負(fù)載平衡實(shí)現(xiàn)，需要在上述算法基礎(chǔ)上構(gòu)造具有針對(duì)性的負(fù)載平衡策略。

對(duì)于結(jié)構(gòu)網(wǎng)格，當(dāng)原始網(wǎng)格子塊數(shù)目不滿足圖剖分要求或各子塊網(wǎng)格單元數(shù)差異大時(shí)，剖分往往難以獲得理想的負(fù)載平衡效果。解決該問(wèn)題的一種有效途徑是對(duì)原始網(wǎng)格的細(xì)分處理，即通過(guò)分裂網(wǎng)格子塊，增加子塊數(shù)目并使得各子塊的網(wǎng)格規(guī)模減小。

而隨著子塊的不斷細(xì)分，網(wǎng)格塊數(shù)逐漸增加，形成的大量網(wǎng)格子塊將給計(jì)算效率、穩(wěn)定性和精度帶來(lái)影響[17]。在采用虛網(wǎng)格進(jìn)行內(nèi)對(duì)接邊界傳值情況下，大量的內(nèi)對(duì)接邊界將導(dǎo)致內(nèi)存、通信量的快速增加，并且在采用隱式時(shí)間推進(jìn)時(shí)，網(wǎng)格塊的過(guò)于細(xì)化將弱化隱式求解的優(yōu)勢(shì)，不利于收斂。此外，某些高階或緊致格式求解時(shí)，內(nèi)對(duì)接邊界本身也會(huì)帶來(lái)精度損失。因此，結(jié)構(gòu)網(wǎng)格的并行剖分需要權(quán)衡負(fù)載平衡效果與剖分后網(wǎng)格拓?fù)溆绊?，而這與細(xì)分處理方法息息相關(guān)。

針對(duì)上述情況，借鑒Rantakokko[8]的多層次剖分思想，本文提出了一種改進(jìn)的負(fù)載平衡方法，可適用于多塊結(jié)構(gòu)對(duì)接網(wǎng)格，并將其推廣至多塊結(jié)構(gòu)重疊網(wǎng)格。該方法主要包括以下3個(gè)方面：全體塊的圖剖分實(shí)現(xiàn)、子塊的分裂處理以及對(duì)上述2個(gè)過(guò)程的循環(huán)調(diào)用實(shí)現(xiàn)。

3.1 圖剖分實(shí)現(xiàn)

圖剖分是本文方法的重要組成，通過(guò)調(diào)用通用圖剖分軟件包，可直接利用現(xiàn)有成熟算法進(jìn)行網(wǎng)格塊的剖分。因此，本文該部分的工作在于建立適用于圖剖分算法的無(wú)向加權(quán)圖。

針對(duì)多塊結(jié)構(gòu)對(duì)接網(wǎng)格，將每個(gè)計(jì)算網(wǎng)格塊視為圖的頂點(diǎn)vi，頂點(diǎn)的權(quán)重w(vi)取每個(gè)塊的網(wǎng)格單元數(shù)，形成圖的頂點(diǎn)集V及其權(quán)重集Wv；相應(yīng)地，相鄰塊的內(nèi)對(duì)接邊界可建立邊ei，邊的權(quán)重w(ei)取內(nèi)邊界面的單元數(shù)目，形成圖的頂點(diǎn)集E及其權(quán)重集We，從而形成完整連通的無(wú)向加權(quán)圖G。關(guān)于該無(wú)向圖G建立的具體描述可見(jiàn)文獻(xiàn)[18]，此處不贅述。

與對(duì)接網(wǎng)格不同，多物體重疊網(wǎng)格在形成插值邊界前，不同物體網(wǎng)格間不存在連接關(guān)系，此時(shí)所建立的無(wú)向圖是不連通的。

因此，最直觀的方法是將各個(gè)物體計(jì)算網(wǎng)格按照多塊對(duì)接型網(wǎng)格的剖分方式依次處理，此處稱為獨(dú)立型剖分方法。該方法簡(jiǎn)單，算法可通用，但隨進(jìn)程數(shù)增加，存在子塊數(shù)目多且部分權(quán)重過(guò)小的缺陷，對(duì)小網(wǎng)格量物體尤為明顯。

因此，本文提出了另外一種重疊網(wǎng)格剖分方式，即通過(guò)初始網(wǎng)格重疊處理，依據(jù)重疊關(guān)系在兩個(gè)不同物體網(wǎng)格重疊區(qū)域各取一子塊建立連接邊，并按重疊邊界關(guān)系給定邊權(quán)值，如圖2所示，從而保證無(wú)向圖連通并進(jìn)行圖剖分，該方法稱為整體型剖分方法。

圖2 重疊網(wǎng)格邊界單元 Fig.2 Fringe points of overset grid

此外，在整個(gè)負(fù)載平衡實(shí)現(xiàn)過(guò)程中，由于存在網(wǎng)格塊的細(xì)分處理，無(wú)向圖本身是動(dòng)態(tài)的，需要在每次細(xì)分處理后建立新的無(wú)向圖，該建立過(guò)程主要采用遍歷搜索實(shí)現(xiàn)。

3.2 子塊的分裂處理

常用的子塊分裂方法包括結(jié)構(gòu)型剖分方法[19]及REB方法。前者直接均等分裂原始網(wǎng)格塊得到大量子塊，而REB方法則通過(guò)對(duì)最大子塊的最大維數(shù)方向遞歸分裂以獲得均勻子塊。幾何剖分方法簡(jiǎn)單，能快速實(shí)現(xiàn)網(wǎng)格的細(xì)分，但由于無(wú)法預(yù)知分裂后的圖剖分結(jié)果，往往會(huì)為了達(dá)到較好負(fù)載平衡效果而明顯增加分裂子塊數(shù)。

因此，借鑒REB方法的遞歸思想，本文對(duì)子塊分裂方法進(jìn)行改進(jìn)，將根據(jù)各子塊權(quán)重與平均權(quán)重關(guān)系來(lái)動(dòng)態(tài)地調(diào)整子塊選取及分裂處理。子塊分裂方法的具體改進(jìn)工作主要圍繞分裂子塊選取原則及分割方式兩個(gè)方面展開(kāi)。

在網(wǎng)格細(xì)化過(guò)程中，為選取更合適的分裂子塊(為便于后文描述，此處稱其為尋優(yōu)子塊)，其選取原則設(shè)計(jì)如下：

1) 求得每個(gè)進(jìn)程的平均網(wǎng)格量Wave，并根據(jù)負(fù)載不均衡率上限ε0，給出子塊理想權(quán)重范圍Wideal=[Wave(2-ε0),Waveε0]，同時(shí)為避免子塊權(quán)重過(guò)小，給定權(quán)重下限值Wmin=Widealσmin。

2) 在每次分裂開(kāi)始前，對(duì)所有子塊編號(hào)及賦權(quán)重，并根據(jù)子塊權(quán)重由大到小進(jìn)行排序{w1(vi),w2(vi),…,wn(vi)}。

3) 根據(jù)步驟2)排序表，由大到小進(jìn)行判斷：當(dāng)w1(vi)∈Wideal時(shí)，該子塊無(wú)需分裂，選擇下一子塊繼續(xù)判斷，依此類推，直到找到某子塊權(quán)重滿足wj(vi)?Wideal且wj(vi)>Wmin，則判定該子塊需要分裂；否則不進(jìn)行子塊分裂。

針對(duì)子塊的分割方式，區(qū)別于REB方法的最大維均等二分處理，此處將依據(jù)子塊拓?fù)浣Y(jié)構(gòu)及權(quán)重關(guān)系進(jìn)行自適應(yīng)分割，其具體過(guò)程為

1) 當(dāng)Wmin

2) 當(dāng)w(vi)>Wave+Wmin時(shí)，按由大到小順序記子塊3個(gè)方向維數(shù)分別為N1、N2、N3，則最大維數(shù)為Nmax=N1。

5) 若步驟4)中條件都不滿足，則按步驟1)的分裂方式進(jìn)行處理。

3.3 圖剖分與子塊分裂的循環(huán)調(diào)用算法

在圖剖分實(shí)現(xiàn)及子塊分裂處理基礎(chǔ)上，建立合理高效的調(diào)用算法是實(shí)現(xiàn)負(fù)載平衡的關(guān)鍵之一。本文方法采用圖剖分算法及子塊分裂處理的循環(huán)調(diào)用，在多層次的細(xì)化過(guò)程中實(shí)現(xiàn)負(fù)載均衡度與通信量的優(yōu)化，且盡可能減少網(wǎng)格子塊數(shù)。

圖3給出了該算法主要流程圖，其中n為網(wǎng)格塊數(shù)，p為計(jì)算進(jìn)程數(shù)，wmax為最大子塊權(quán)重，ε為剖分后負(fù)載不均衡率，定義為單個(gè)進(jìn)程最大負(fù)載量與平均負(fù)載量之比，此處ε0默認(rèn)取為1.05。

圖3 剖分算法流程圖 Fig.3 Flow diagram of partition algorithm

4 算例及結(jié)果分析

為驗(yàn)證本文提出的負(fù)載平衡方法的有效性，采用2個(gè)典型算例進(jìn)行了剖分試驗(yàn)，從負(fù)載不均衡率、網(wǎng)格子塊總數(shù)、邊界面單元總數(shù)(Edge-total)及通信界面單元數(shù)(Edge-cut)4個(gè)方面對(duì)比了不同網(wǎng)格類型下剖分方法的剖分效果。

在對(duì)接網(wǎng)格算例中，所用方法包括圖剖分-最大子塊最大維均等二分法(記為GS1)、圖剖分-尋優(yōu)子塊最大維均等二分法(記為GS2)以及圖剖分-尋優(yōu)子塊自適應(yīng)二分法(記為GS3)，通過(guò)算例分析子塊分裂方法對(duì)負(fù)載平衡的影響，并與采用傳統(tǒng)負(fù)載平衡算法及不采用負(fù)載平衡的結(jié)果進(jìn)行了部分核數(shù)下的并行效率對(duì)比。所采用的傳統(tǒng)負(fù)載平衡算法為裝箱算法[20]，即BPA(Bin-Packing Algorithm)。

在此基礎(chǔ)上，結(jié)合重疊網(wǎng)格剖分特點(diǎn)，比較了獨(dú)立型剖分方法(分別記為GS1-sep、GS3-sep)以及整體型剖分方法(分別記為GS1-int、GS3-int)在多塊結(jié)構(gòu)重疊網(wǎng)格算例中的表現(xiàn)。

4.1 三段翼30P30N翼型

本算例選取高升翼型30P30N[21]進(jìn)行負(fù)載平衡方法驗(yàn)證與對(duì)比，考察不同進(jìn)程數(shù)下的負(fù)載平衡效果，并進(jìn)行了部分狀態(tài)的并行計(jì)算。計(jì)算條件為：馬赫數(shù)Ma=0.2、雷諾數(shù)Re=9×106、迎角α=4°。

原始網(wǎng)格總網(wǎng)格量約為855.1萬(wàn)，共包括13個(gè)子塊(拓?fù)浣Y(jié)構(gòu)見(jiàn)圖4)，且各子塊網(wǎng)格量分布如圖5所示，圖中橫軸坐標(biāo)為網(wǎng)格塊編號(hào)，縱軸坐標(biāo)為網(wǎng)格量，其中，最大子塊網(wǎng)格量約占總網(wǎng)格量的18.6%。

圖4 30P30N翼型網(wǎng)格拓?fù)?Fig.4 Topology of 30P30N airfoil grid

圖6給出了GS1、GS2、GS3這3種方法及傳統(tǒng)裝箱算法在不同進(jìn)程數(shù)下的剖分結(jié)果，圖中橫坐標(biāo)為進(jìn)程數(shù)。圖6(a)為負(fù)載不均衡率隨進(jìn)程數(shù)變化曲線，可以看到，本文3種負(fù)載平衡方法均達(dá)到負(fù)載不均衡率要求(ε≤1.05)，其中GS3整體略優(yōu)，而B(niǎo)PA所得負(fù)載不均衡率明顯偏高。圖6(b)則為網(wǎng)格子塊總數(shù)隨進(jìn)程數(shù)變化曲線，可見(jiàn)，隨著進(jìn)程數(shù)增加，GS1所得網(wǎng)格子塊總數(shù)快速增加，GS2方法與BPA表現(xiàn)類似，略有減少，而GS3維持著較少且平緩的增長(zhǎng)趨勢(shì)。圖6(c)為邊界面單元總數(shù)及通信界面單元數(shù)隨進(jìn)程數(shù)變化曲線。由于邊界面單元總數(shù)與網(wǎng)格子塊數(shù)正相關(guān)，因此該量的變化呈現(xiàn)出與子塊總數(shù)類似的變化趨勢(shì)。而在通信界面單元數(shù)方面，如圖6(d)所示，GS1、GS2兩種方法的增長(zhǎng)趨勢(shì)基本一致，GS3方法略有優(yōu)勢(shì)，這預(yù)示著，通過(guò)圖剖分算法的通信優(yōu)化，GS1及GS2方法所產(chǎn)生的大量?jī)?nèi)邊界面單元將分配在同一進(jìn)程，從而避免了進(jìn)程間通信量的急劇增加。相比之下，BPA由于缺少對(duì)進(jìn)程間通信量的優(yōu)化，產(chǎn)生了更多的通信界面單元，即意味著更高的通信成本。

圖5 30P30N翼型網(wǎng)格量分布 Fig.5 Distribution of number of grid points of 30P30N airfoil

圖6 30P30N翼型不同進(jìn)程數(shù)時(shí)網(wǎng)格剖分結(jié)果 Fig.6 Partitioning properties as a function of processors for 30P30N airfoil grids

由此可見(jiàn)，本文3種負(fù)載平衡方法皆能滿足負(fù)載平衡要求，而對(duì)比結(jié)果表明，新的子塊選取原則能一定程度地減少剖分子塊數(shù)，且改進(jìn)后的子塊分割方式體現(xiàn)出了更明顯的作用，能夠極大地避免子塊數(shù)過(guò)多的缺陷，進(jìn)而降低因虛網(wǎng)格而額外增加的內(nèi)存需求。

為比較剖分方法對(duì)計(jì)算效率的影響，且考慮到原始網(wǎng)格塊數(shù)限制，對(duì)小核數(shù)情形進(jìn)行了并行計(jì)算，對(duì)比不同負(fù)載平衡方法及不采用負(fù)載平衡的并行影響。

圖7及圖8分別為小核數(shù)時(shí)不同剖分方法及不采用負(fù)載平衡(記為NO-LB)的負(fù)載不平衡率和并行加速比曲線，其中Np為核數(shù)，Nb為原始網(wǎng)格塊數(shù)。從圖中可見(jiàn)，不采用負(fù)載平衡時(shí)，其負(fù)載不平衡率明顯偏高，并對(duì)并行加速比產(chǎn)生了十分不利的影響。而相比于傳統(tǒng)裝箱算法，本文的3種方法皆能獲得更優(yōu)的負(fù)載不均衡率及并行加速比，其中，GS3方法表現(xiàn)最優(yōu)。由此可見(jiàn)，本文的新方法在并行計(jì)算效率方面具有更佳效果。

圖7 負(fù)載不均衡率曲線(Np

圖8 并行加速比曲線(Np

為進(jìn)一步研究網(wǎng)格塊數(shù)對(duì)數(shù)值計(jì)算的影響，對(duì)進(jìn)程數(shù)為128時(shí)的GS1、GS3及BPA 3種方法的剖分網(wǎng)格進(jìn)行了并行計(jì)算，并與原始網(wǎng)格串行結(jié)果進(jìn)行對(duì)比。圖9為30P30N翼面壓力系數(shù)計(jì)算結(jié)果(Orgin代表原始網(wǎng)格)，圖中橫坐標(biāo)為機(jī)翼表面坐標(biāo)x與弦長(zhǎng)c之比，縱坐標(biāo)Cp為壓力系數(shù)，可見(jiàn)，不同網(wǎng)格結(jié)果與實(shí)驗(yàn)值吻合較好，證實(shí)了本文數(shù)值計(jì)算的正確性。

表1為3套網(wǎng)格的負(fù)載不均衡率及子塊總數(shù)，并給出了并行計(jì)算加速比及預(yù)測(cè)阻力系數(shù)CD?？梢钥吹?，GS1和GS3方法所得負(fù)載不均衡率保持一致，但子塊總數(shù)相差較大，而B(niǎo)PA所得負(fù)載不均衡率及子塊總數(shù)均較大。從并行加速比來(lái)看，GS3剖分方法更有優(yōu)勢(shì)，BPA所得并行加速比最低。阻力系數(shù)的預(yù)測(cè)結(jié)果表明網(wǎng)格拓?fù)鋵?duì)計(jì)算產(chǎn)生了影響，過(guò)于細(xì)分的網(wǎng)格使得阻力系數(shù)偏離原始網(wǎng)格串行結(jié)果更明顯。圖10及圖11則分別給出了GS1、GS3 2種方法及原始網(wǎng)格的計(jì)算殘差及阻力系數(shù)收斂曲線，可以看到，網(wǎng)格的細(xì)分使得收斂性有所變差，GS1方法表現(xiàn)得更為明顯。

圖9 30P30N翼面壓力系數(shù)分布 Fig.9 Pressure coefficients distribution on surfaces of 30P30N airfoil

表1 進(jìn)程數(shù)為128時(shí)的剖分結(jié)果Table 1 Partition results for 128 processors

MethodLoadimba?lanceratioNumberofblocksSpeedupCDGS11．04762490．340．0288GS31．04716794．410．0292BPA1．09643786．910．0289Origin131．000．0293

圖10 殘差收斂曲線 Fig.10 Convergence curves of residual

圖11 阻力系數(shù)收斂曲線 Fig.11 Convergence curves of drag coefficient

4.2 翼身組合體DLR-F6構(gòu)型

為進(jìn)一步驗(yàn)證本文剖分算法對(duì)復(fù)雜重疊網(wǎng)格適用性，選取德國(guó)宇航公司的DLR-F6翼身組合體[22](DLR-F6 WB)的重疊網(wǎng)格作為方法驗(yàn)證算例，考察不同剖分方式的負(fù)載平衡效果。

整套重疊網(wǎng)格由3個(gè)部分組成，包括機(jī)身網(wǎng)格、機(jī)翼網(wǎng)格及背景網(wǎng)格，并且每個(gè)部分由多塊結(jié)構(gòu)網(wǎng)格構(gòu)成。圖12為DLR-F6機(jī)翼與機(jī)身網(wǎng)格重疊關(guān)系示意。原始網(wǎng)格總單元數(shù)約為1 700萬(wàn)，包含19個(gè)子塊，且各子塊網(wǎng)格量分布如圖13所示，其中最大子塊網(wǎng)格量約為314萬(wàn)。

圖14為5種負(fù)載平衡方法在不同進(jìn)程數(shù)下的剖分結(jié)果。圖14(a)縱軸坐標(biāo)為負(fù)載不均衡率，可見(jiàn)，除BPA外，其余4種方法均獲得了滿足負(fù)載平衡要求的剖分結(jié)果。而圖14(b)給出了子塊總數(shù)的變化曲線，可明顯看出：子塊分裂方法及整體型剖分方式皆影響顯著，隨進(jìn)程數(shù)增加，GS3-int方法能極大地減少網(wǎng)格塊數(shù)。而B(niǎo)PA雖然具有與GS3-sep方法類似表現(xiàn)，但其負(fù)載不均衡率難以滿足剖分要求。

圖14(c)為邊界面單元數(shù)隨進(jìn)程數(shù)變化曲線，與網(wǎng)格子塊數(shù)趨勢(shì)類似，且GS3-int方法的邊界面單元數(shù)明顯少于其余4種方法。圖14(d)則為進(jìn)程間通信界面單元總數(shù)變化曲線，可以看到，兩種整體型剖分方法表現(xiàn)一致，且明顯優(yōu)于獨(dú)立型剖分方法及傳統(tǒng)裝箱算法。由此可見(jiàn)，剖分方式是影響該量的主導(dǎo)因素，而采用整體型剖分方法能夠更有效地降低進(jìn)程間通信成本及內(nèi)存需求。

圖12 DLR-F6機(jī)翼與機(jī)身構(gòu)型網(wǎng)格重疊 Fig.12 Overset grid for DLR-F6 wing and body configuration

圖13 DLR-F6 WB構(gòu)型網(wǎng)格量分布 Fig.13 Distribution of number of grid points of DLR-F6 WB configuration

圖14 DLR-F6 WB構(gòu)型不同進(jìn)程數(shù)時(shí)網(wǎng)格剖分結(jié)果 Fig.14 Partitioning properties as a function of processors for DLR-F6 WB configuration grids

為直觀地說(shuō)明本文4種新方法的剖分思路，圖15給出了進(jìn)程數(shù)為128時(shí)，不同方法所得網(wǎng)格子塊網(wǎng)格量分布情況，圖中橫軸為單個(gè)子塊網(wǎng)格量與平均負(fù)載網(wǎng)格量之比，按比值范圍歸為4類，縱軸為各類子塊數(shù)占全部子塊數(shù)的比重。從圖中看出，為實(shí)現(xiàn)負(fù)載均衡，GS1類方法剖分出大量小網(wǎng)格量子塊，而GS3-sep方法則保留了更多中等網(wǎng)格量子塊。形成鮮明對(duì)比的是，GS3-int方法通過(guò)對(duì)全部子塊進(jìn)行篩選分裂，保留了多數(shù)滿足平均負(fù)載網(wǎng)格量的大網(wǎng)格量子塊，并分裂未滿足要求的其余子塊，以實(shí)現(xiàn)較少塊數(shù)下的負(fù)載均衡。

圖15 各子塊網(wǎng)格量比重分布 Fig.15 Proportion for subgrids distribution of different sizes

5 結(jié) 論

本文提出了一種改進(jìn)的多層次圖剖分的多塊結(jié)構(gòu)網(wǎng)格負(fù)載平衡方法。該方法設(shè)計(jì)了新的剖分算法，采用網(wǎng)格子塊分割與圖剖分的循環(huán)調(diào)用思想，在逐步細(xì)化過(guò)程中實(shí)現(xiàn)負(fù)載平衡及通信量?jī)?yōu)化，同時(shí)改進(jìn)了子塊分裂方法，避免了子塊數(shù)目過(guò)多的缺陷，并通過(guò)建立不同物體間重疊邊界關(guān)聯(lián)，將其拓展至多塊結(jié)構(gòu)重疊網(wǎng)格。

1) 該方法剖分過(guò)程無(wú)需人工干預(yù)，與網(wǎng)格拓?fù)錈o(wú)關(guān)，適用于復(fù)雜問(wèn)題的多塊結(jié)構(gòu)網(wǎng)格。

2) 基于圖剖分的循環(huán)調(diào)用算法能夠很好地實(shí)現(xiàn)負(fù)載均衡及通信成本優(yōu)化，提高并行效率。

3) 子塊分裂方法對(duì)剖分結(jié)果影響顯著，改進(jìn)的子塊分裂方法能夠極大地減少網(wǎng)格塊數(shù)及由內(nèi)邊界虛網(wǎng)格所帶來(lái)的內(nèi)存需求。

4) 對(duì)于多塊結(jié)構(gòu)重疊網(wǎng)格，采用整體型剖分方法能夠有效降低進(jìn)程間通信成本，顯著減少網(wǎng)格塊數(shù)及內(nèi)邊界面單元數(shù)。

基于現(xiàn)有的工作，下一步考慮將該方法發(fā)展并應(yīng)用于復(fù)雜多體分離等數(shù)值模擬的并行計(jì)算。

[1] CHIEN Y P, CARPENTER F, ECER A, et al. Load-balancing for parallel computation of fluid dynamics problems[J]. Computer Methods in Applied Mechanics & Engineering, 1995, 120(1-2): 119-130.

[2] RANTAKOKKO J. A framework for partitioning structured grids with inhomogeneous workload[J]. Parallel Algorithms and Applications, 1998, 13(2): 135-151.

[3] MILLER G L, TENG S H, VAVASIS S A. A unified geometric approach to graph separators[C]//Proceedings of 32nd Annual Symposium of Foundations of Computer Science. Piscataway, NJ: IEEE Press， 1991: 538-547.

[4] SIMON H. Partitioning unstructured problems for parallel processing[J]. Computing Systems in Engineering, 1991, 11(4): 135-148.

[5] KARYPIS G, KUMAR V. METIS: Unstructured graph partitioning and sparse matrix ordering system, Technical Report[R]. Minneapolis: University of Minnesota, 1995.

[6] BOMAN E G, ?ATALYüREK U V, CHEVALIER C, et al. The Zoltan and Isorropia parallel toolkits for combinatorial scientific computing: Partitioning, ordering and coloring[J]. Scientific Programming, 2012, 20(2): 129-150.

[7] YTTERSTR?M A. A tool for partitioning structured multiblock meshes for parallel computational mechanics[J]. International Journal of High Performance Computing Applications, 1997, 11(4): 336-343.

[8] RANTAKOKKO J. Partitioning strategies for structured multiblock grids[J]. Parallel Computing, 2000, 26(12): 1661-1680.

[9] AHUSBORDE E, GLOCKNER S. A 2D block-structured mesh partitioner for accurate flow simulations on non-rectangular geometries[J]. Computers & Fluids, 2011, 43(1): 2-13.

[10] DJOMEHRI M J, BISWAS R. Performance enhancement strategies for multi-block overset grid CFD applications[J]. Parallel Computing, 2003, 29(11-12): 1791-1810.

[11] 司海青, 王同光. 多塊并行計(jì)算中負(fù)載平衡策略及時(shí)間成本估算方法[J]. 航空學(xué)報(bào), 2007, 28(S1): 57-61.

SI H Q, WANG T G. Load balancing strategy for parallel calculation and time cost estimation[J]. Acta Aeronautica et Astronautica Sinica, 2007, 28(S1): 57-61 (in Chinese).

[12] 李桂波, 楊國(guó)偉. 基于多塊結(jié)構(gòu)網(wǎng)格的并行計(jì)算及負(fù)載平衡研究[J]. 宇航學(xué)報(bào), 2011, 32(6): 1224-1230.

LI G B, YANG G W. Study on parallel computation and load balance strategy based on multiblock structured grid[J]. Journal of Astronautics, 2011, 32(6): 1224-1230 (in Chinese).

[13] 許正, 李津, 朱自強(qiáng). 網(wǎng)絡(luò)連接機(jī)群上CFD計(jì)算的一種負(fù)載平衡方法[J]. 航空學(xué)報(bào), 2005, 26(2):129-134.

XU Z, LI J, ZHU Z Q. Load balancing strategy for parallel CFD calculation on cluster[J]. Acta Aeronautica et Astronautica Sinica, 2005, 26(2):129-134 (in Chinese).

[14] BARNARD S T, SIMON H D. Fast multilevel implementation of recursive spectral bisection for partitioning unstructured problems[J]. Concurrency Practice & Experience, 1994, 6(6): 101-117.

[15] KARYPIS G, KUMAR V. A fast and high quality multi-level scheme for partitioning irregular graphs[J]. SIAM Journal on Scientific Computing, 2006, 20(1): 359-392.

[16] KIRICHENKO A V, MASON K, STRAUME M, et al. An efficient K-way graph partitioning algorithm for task allocation in parallel computing systems[C]//Proceedings of the First International Conference on Systems Integration. Piscataway, NJ: IEEE Press, 1990: 748-751.

[17] GOURDAIN N, GICQUEL L, MONTAGNAC M, et al. High performance parallel computing of flows in complex geometries: I. Methods[J]. Computational Science & Discovery, 2009, 339(1): 104-124.

[18] PELLEGRINI F. Graph partitioning based methods and tools for scientific computing[J]. Parallel Computing, 1997, 23(1-2): 153-164.

[19] 劉巍, 張理論, 王勇獻(xiàn), 等. 計(jì)算空氣動(dòng)力學(xué)并行編程基礎(chǔ)[M]. 北京: 國(guó)防工業(yè)出版社, 2013: 250-255.

LIU W, ZHANG L L, WANG Y X, et al. Foundations of computational aerodynamics parallel programming[M]. Beijing: National Defense Industry Press, 2013: 250-255 (in Chinese).

[20] DJOMEHRI M J, BISWAS R, LOPEZ-BENITEZ N. Load balancing strategies for multi-block overset grid applications[C]//Proceedings of the ISCA 18th International Conference Computers and Their Applications, 2003:63-71.

[21] COOK P, MCDONALD M, FIRMIN M. Airfoil REA2822 pressure distributions, and boundary layer and wake measurements, experimental data base for counter program assessment: AGARD Report AR 138[R]. Paris: AGARD, 1979.

[22] 2nd AIAA CFD Drag Prediction Workshop. [2016-06-27]. http://aaac.larc.nasa.gov/ts-ab/cfdlarc-dpw/June 2003.

(責(zé)任編輯： 李明敏)

URL：www.cnki.net/kcms/detail/11.1929.V.20160812.1059.004.html

Loadbalancestrategybasedongraphpartitionformultiblockstructuredgrids

LIUHongkang,YANChao*,LINBoxi,ZHAOYatian

SchoolofAeronauticScienceandEngineering,BeihangUniversity,Beijing100083,China

Loadbalanceisofsignificancetotheperformanceofparallelcomputing.Aimingatagoodloadbalancewithaslessblocksaspossibleinparallelcomputing,anenhancedpartitioningstrategybasedonmultilevelgraphpartitionisproposedformultiblockstructuredgrids,includingarecursivepartitionalgorithmandanimprovedsubgrid-splittingmethod,andthenextendedtotheoverlappingmultiblockgridsbyestablishingtheconnectionbetweensubgridsofdifferentbodieshereafter.Twotypicalapplications,coveringthe1to1coincidentgridandtheoverlappinggrid,areimplementedtocomparethebehaviorsofvariouspartitionstrategies,asregardsloadbalance,edge-cutandblocknumbers.Resultsdemonstratethatthesubgrid-splittingmethodiscriticaltostructuredgrids,andpartitioningover-lappinggridintegrallyisobviouslyabetteralternative.Specifically,thenewpartitioningstrategyshowsagoodperformanceinloadbalanceandcommunicationoverheads,andmeanwhiledecreasestheamountofblocksenormouslyaswellasthememoryrequirementcausedbytheghostcellofedge-cuts,leadingtoabetterparallelefficiency.Theenhancedpartitionstrategyisapplicabletoboththe1to1coincidentgridsandoverlappinggrids.

computationalfluiddynamics;parallelcomputing;structuredgrids;loadbalance;graphpartition

2016-06-27；Revised2016-07-25；Accepted2016-08-08；Publishedonline2016-08-121059

.E-mailyanchao@buaa.edu.cn

2016-06-27；退修日期2016-07-25；錄用日期2016-08-08； < class="emphasis_bold">網(wǎng)絡(luò)出版時(shí)間

時(shí)間：2016-08-121059

www.cnki.net/kcms/detail/11.1929.V.20160812.1059.004.html

.E-mailyanchao@buaa.edu.cn

劉宏康， 閻超， 林博希， 等． 基于圖剖分的多塊結(jié)構(gòu)網(wǎng)格負(fù)載平衡方法J． 航空學(xué)報(bào),2017,38(5):120558.LIUHK,YANC,LINBX,etal.LoadbalancestrategybasedongraphpartitionformultiblockstructuredgridsJ.ActaAeronauticaetAstronauticaSinica,2017,38(5):120558.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0230

V211.3

A

1000-6893(2017)05-120558-10