昌繼海,趙銘偉,余 飛,關(guān)振群*
(1.大連理工大學(xué)工程力學(xué)系,遼寧 大連 116024;2.大連理工大學(xué)電子信息與電氣工程學(xué)部,遼寧 大連 116024)
在求解邊界運(yùn)動(dòng)等非定常流動(dòng)問(wèn)題時(shí),針對(duì)網(wǎng)格變形問(wèn)題的并行化研究是求解幾何變形問(wèn)題的重要方向[1]。網(wǎng)格變形方法主要分為插值類方法和彈簧類方法[2]兩大類。
插值類的網(wǎng)格變形方法將流體內(nèi)部網(wǎng)格節(jié)點(diǎn)的位移視為邊界節(jié)點(diǎn)的插值問(wèn)題,該類方法由于在變形過(guò)程中不考慮節(jié)點(diǎn)之間的連接性,直接通過(guò)邊界節(jié)點(diǎn)的位移對(duì)內(nèi)部節(jié)點(diǎn)進(jìn)行插值,不需要迭代,容易實(shí)現(xiàn)并行化,許多研究者都對(duì)插值類方法的并行化進(jìn)行了研究。Boer[3]等人提出的徑向基函數(shù)法(Radial Basis Function, RBF)是一種易于并行的插值方法,該方法通過(guò)選取特定的徑向基函數(shù),利用邊界節(jié)點(diǎn)對(duì)內(nèi)部節(jié)點(diǎn)進(jìn)行插值得到內(nèi)部節(jié)點(diǎn)的位移;Rendall[4]和曾[5]等人在RBF方法的基礎(chǔ)上進(jìn)行了改進(jìn),如通過(guò)貪婪算法減少插值點(diǎn)的數(shù)目等等。Li[6]和Fang[7]等人提出了基于主從架構(gòu)的并行方法來(lái)提高RBF法的并行效率。此外,劉[8]等人提出的背景網(wǎng)格法也是一種易于并行的插值類方法,但是該方法面臨凹邊界的幾何問(wèn)題或者推廣到三維問(wèn)題時(shí)會(huì)遇到困難。
物理類方法主要包括Batina[9]等人提出的彈簧類方法,該方法將網(wǎng)格類比為一個(gè)由彈簧組成的彈性系統(tǒng)。相比于插值類變形方法,對(duì)彈簧類變形方法的并行化研究較少,一方面是由于彈簧類方法變形能力較弱;另一反面,彈簧類方法在計(jì)算過(guò)程中還要始終考慮節(jié)點(diǎn)之間的連接性,需要求解大規(guī)模的線性方程組,導(dǎo)致計(jì)算效率較低。
一些研究者提出了許多改進(jìn)的彈簧類方法,如增加扭轉(zhuǎn)彈簧[10]或垂直彈簧[11]等方式來(lái)提高變形能力?;谶@些改進(jìn)的彈簧類方法,Ma[12]和程[13]等人進(jìn)行了分布式并行化的研究,取得了較好的并行效率。
林[14]等人在Ball-vertex彈簧法的基礎(chǔ)上提出一種局部插值的點(diǎn)球彈簧修勻法(VerBSS)方法,該方法通過(guò)分解彈簧系統(tǒng)的方式對(duì)網(wǎng)格進(jìn)行逐點(diǎn)求解,具有插值類方法的特征,是一種求解效率較高的網(wǎng)格變形方法。
本文在林等人的基礎(chǔ)上,采用分區(qū)和共享單元層的策略,提出了一種分布式的VerBSS并行方案,既保留了該算法較高的變形能力,又實(shí)現(xiàn)了較高加速效率的并行化。最后,通過(guò)算例對(duì)本文提出的VerBSS并行算法和共享式并行方法的變形能力進(jìn)行了對(duì)比,在二維和三維情形下對(duì)該算法的加速性能進(jìn)行了測(cè)試。算例結(jié)果表明這種并行方案變形能力較強(qiáng),具有較高的并行效率,適合推廣到大變形、大規(guī)模的網(wǎng)格變形問(wèn)題。
普通彈簧法的原理,是將網(wǎng)格的邊看作彈簧,其剛度與其邊的長(zhǎng)度成反比。相鄰的節(jié)點(diǎn)i和節(jié)點(diǎn)j之間的力fij可表示為
(1)
ui和uj分別為節(jié)點(diǎn)i和節(jié)點(diǎn)j的位移,kij為ij彈簧的剛度,nij為節(jié)點(diǎn)i到節(jié)點(diǎn)j的單位向量。為避免這種現(xiàn)象的出現(xiàn),通過(guò)在三角形節(jié)點(diǎn)與其對(duì)邊之間插入垂直彈簧來(lái)阻止單元塌陷,如圖1a)所示,p點(diǎn)為i節(jié)點(diǎn)到邊Ejk的垂線的垂足。這種垂直彈簧稱為Ball-vertex彈簧。
圖1 Ball-vertex彈簧法
Ball-vertex彈簧作用在節(jié)點(diǎn)i上的彈簧力可表示為
(2)
其中up和ui分別為節(jié)點(diǎn)p和節(jié)點(diǎn)i的位移,kip為ip彈簧的剛度,nip為節(jié)點(diǎn)i到節(jié)點(diǎn)p的單位向量。p點(diǎn)的位移可以通過(guò)節(jié)點(diǎn)j和節(jié)點(diǎn)k的位移插值獲得,為
up=ξuj+(1-ξ)uk
(3)
其中ξ為節(jié)點(diǎn)在邊Ejk上的插值系數(shù)。因此,ip彈簧作用在p點(diǎn)上的力可分解為作用在j點(diǎn)和k點(diǎn)上的分力。則節(jié)點(diǎn)i上的彈簧系統(tǒng)線性方程組為
(4)
其中K為節(jié)點(diǎn)i的剛度矩陣,u為子彈簧系統(tǒng)各節(jié)點(diǎn)的位移向量,b為i點(diǎn)的載荷向量。
VerBSS方法的求解方式不同于普通的彈簧方法,該方法通過(guò)在每個(gè)內(nèi)部節(jié)點(diǎn)的閉包上形成一個(gè)以該節(jié)點(diǎn)為核心,以所有與該節(jié)點(diǎn)相鄰接的節(jié)點(diǎn)為外圍邊界的多邊形或多面體閉包,如圖1b)所示。
VerBSS算法:
輸入:初始網(wǎng)格M0
輸出:變形網(wǎng)格M1
Step1:導(dǎo)入初始網(wǎng)格,查找所有內(nèi)部節(jié)點(diǎn),建立點(diǎn)及其周圍節(jié)點(diǎn)的拓?fù)潢P(guān)系。
Step2:為每個(gè)內(nèi)部節(jié)點(diǎn)構(gòu)建多面體包腔的子彈簧系統(tǒng),并對(duì)剛度矩陣K進(jìn)行LDLT分解,將分解后的L矩陣和D矩陣,以及轉(zhuǎn)換后的右端項(xiàng)存儲(chǔ)起來(lái)。
Step3:將已知邊界的節(jié)點(diǎn)位移u代入式(4),遍歷網(wǎng)格內(nèi)部節(jié)點(diǎn)進(jìn)行計(jì)算,并統(tǒng)計(jì)內(nèi)部節(jié)點(diǎn)位移量變化,作為收斂條件。
Step4:判斷Step3是否滿足收斂條件,若不收斂,重復(fù)step3,直至收斂。
本文采用METIS分區(qū)工具將初始網(wǎng)格劃分為若干個(gè)分區(qū)[15],進(jìn)行任務(wù)劃分。每個(gè)子分區(qū)相應(yīng)的節(jié)點(diǎn)和單元的數(shù)據(jù)存儲(chǔ)在各自的存儲(chǔ)區(qū)域,通過(guò)拓寬相鄰分區(qū)邊界的方式來(lái)形成共享單元層,使得這些分區(qū)邊界上的節(jié)點(diǎn)被納入到求解序列中,如圖2a)和圖2b)所示為初始網(wǎng)格和分區(qū)后的網(wǎng)格。
圖2 二維機(jī)翼流場(chǎng)網(wǎng)格分區(qū)
選擇一個(gè)分區(qū)將其邊界上所有與邊界線相連接的單元的幾何信息傳遞給其鄰接分區(qū),那么這些共享單元線將拓寬為一層共享單元面,相鄰的兩個(gè)分區(qū)共享這層單元面,共享單元層如圖2c)所示。三維情形中的共享單元面將拓寬為一層共享單元墻,相鄰的兩個(gè)分區(qū)共享這層單元墻。通過(guò)構(gòu)建相鄰分區(qū)之間的共享單元層,每個(gè)計(jì)算分區(qū)的節(jié)點(diǎn)都被分為兩類,一類為普通節(jié)點(diǎn),一類為共享單元層上的節(jié)點(diǎn)。共享單元層上的節(jié)點(diǎn)可分為該分區(qū)的發(fā)送塊和接收塊。通過(guò)與其相鄰的分區(qū)建立對(duì)應(yīng)的發(fā)送和接收關(guān)系,對(duì)通信的發(fā)送地址、接收地址和通信數(shù)據(jù)的類型、數(shù)目以及其它通信的信息進(jìn)行初始化,將相鄰分區(qū)的通信渠道建立起來(lái)。
具體的通信過(guò)程如圖3a)所示,共享單元層上的a節(jié)點(diǎn)和b節(jié)點(diǎn)同時(shí)為左分區(qū)和右分區(qū)所共享。左分區(qū)只計(jì)算a節(jié)點(diǎn),不計(jì)算b節(jié)點(diǎn),而右分區(qū)只計(jì)算b節(jié)點(diǎn),不計(jì)算a節(jié)點(diǎn)。左分區(qū)在求得a節(jié)點(diǎn)的位移之后將a節(jié)點(diǎn)的位移傳遞給右分區(qū),右分區(qū)的計(jì)算同理,如圖3b)所示。
圖3 共享單元層節(jié)點(diǎn)的閉包圖
VerBSS的并行算法步驟如下:
輸入:初始分區(qū)網(wǎng)格
輸出:變形網(wǎng)格
Step1:初始化分區(qū)個(gè)數(shù),每個(gè)網(wǎng)格分區(qū)分配一個(gè)進(jìn)程讀取相應(yīng)的子分區(qū)網(wǎng)格。
Step2:查找相鄰分區(qū)的邊界,并在相鄰分區(qū)中間構(gòu)建共享區(qū)。
Step3:對(duì)每個(gè)進(jìn)程讀取的子分區(qū)網(wǎng)格以及構(gòu)建的共享區(qū)網(wǎng)格進(jìn)行分類,分為發(fā)送塊、接收塊和普通塊,并將每個(gè)進(jìn)程和其它進(jìn)程要通信的發(fā)送地址、接收地址和通信數(shù)據(jù)的類型、數(shù)目等通信信息進(jìn)行初始化。
Step4:每個(gè)進(jìn)程采用VerBSS算法對(duì)該分區(qū)進(jìn)行求解,具體求解過(guò)程參見(jiàn)本文2.2節(jié)中的VerBSS算法流程。
Step5:在各相鄰分區(qū)之間進(jìn)行通信,按照step2建立的通信關(guān)系分別進(jìn)行發(fā)送和接收,每個(gè)進(jìn)程將接收到的位移數(shù)據(jù)整理之后在對(duì)應(yīng)的節(jié)點(diǎn)上更新。
Step6:每個(gè)進(jìn)程對(duì)更新的位移求位移平方差的平均值,判斷其是否滿足給定的收斂條件。若滿足,則停止循環(huán);若不滿足,重復(fù)step4和step5,直到滿足收斂條件為止。
本文算例將在一臺(tái)48核心、128G內(nèi)存的工作站上進(jìn)行測(cè)試(CPU型號(hào)為Inter Xeon(R) E5-2650 v4,主頻2.2G)。
以一個(gè)二維機(jī)翼做俯沖運(yùn)動(dòng)的模型,其節(jié)點(diǎn)數(shù)為603 802,單元數(shù)為1 203 662。機(jī)翼長(zhǎng)度為102個(gè)單位長(zhǎng)度,單元尺寸為1個(gè)單位長(zhǎng)度,單步旋轉(zhuǎn)角度為0.02。如圖4所示為變形前后的對(duì)比圖。
圖4 機(jī)翼俯沖變形前后的對(duì)比圖
采用分區(qū)并行的VerBSS算法單步迭代過(guò)程中各個(gè)步驟的平均耗時(shí)和加速比如表1所示。
表1 二維機(jī)翼旋轉(zhuǎn)分區(qū)并行的單步各步驟耗時(shí)
若采用共享式并行的VerBSS算法,單步迭代過(guò)程中各個(gè)步驟的平均耗時(shí)和加速比如表2所示。
表2 二維機(jī)翼旋轉(zhuǎn)的共享式并行的單步各步驟耗時(shí)
圖5為采用共享式并行方案和分布式并行方案時(shí)的加速效率對(duì)比圖。在單個(gè)變形時(shí)間步的耗時(shí)中,矩陣構(gòu)造時(shí)不需要通信,為完全解耦的過(guò)程,求解步驟和通信步驟交替進(jìn)行,其中求解步驟為并行過(guò)程,通信步驟為串行過(guò)程。在該算例中,當(dāng)分區(qū)數(shù)為16時(shí),加速比為15.58,其并行效率為97.37%。當(dāng)分區(qū)數(shù)為32時(shí),并行的加速比為26.93,并行效率為84.15%。可見(jiàn),由二維模型的算例可知,該并行算法在分區(qū)數(shù)目不超過(guò)16時(shí),并行的加速比和理想加速比相當(dāng)。當(dāng)分區(qū)數(shù)目持續(xù)增大時(shí),達(dá)到收斂條件所需的迭代步數(shù)也相應(yīng)增加,導(dǎo)致并行效率降低,但和共享式的并行方案相比仍然具有較高的加速效率。
圖5 機(jī)翼俯沖模型變形的加速比
對(duì)比普通彈簧法、串行VerBSS算法和本文的并行算法,三種方法變形時(shí)網(wǎng)格的最小夾角隨時(shí)間步的變化如圖6所示,可見(jiàn)在算法并行化的過(guò)程網(wǎng)格單元的質(zhì)量具備良好的魯棒性。
圖6 最小夾角隨時(shí)間步的變化圖
進(jìn)一步驗(yàn)證本文算法在三維情形中的并行效率,一個(gè)三維的機(jī)翼在空間中做翹曲變形的網(wǎng)格模型如圖7所示。圖7a)為8分區(qū)的初始網(wǎng)格內(nèi)部截面圖,圖7b)為變形后的分區(qū)網(wǎng)格內(nèi)部截面圖,其彎曲形狀為拋物線型。
圖7 三維機(jī)翼翹曲變形的流場(chǎng)網(wǎng)格
機(jī)翼周圍的流場(chǎng)長(zhǎng)寬高分別為1 500、500、500個(gè)單位長(zhǎng)度,網(wǎng)格單元尺寸為0.5個(gè)單位長(zhǎng)度,該流場(chǎng)網(wǎng)格包含有132 120個(gè)節(jié)點(diǎn)和656 878個(gè)單元,翼尖變形的單步步長(zhǎng)為0.5個(gè)單位長(zhǎng)度。采用分區(qū)并行的VerBSS算法迭代過(guò)程中各個(gè)步驟的平均耗時(shí)如表3所示。采用共享式并行的VerBSS算法各個(gè)步驟的耗時(shí)如表4所示。
表3 機(jī)翼翹曲分區(qū)并行的單步各步驟耗時(shí)
表4 機(jī)翼翹曲的共享式并行的單步各步驟耗時(shí)
該算例中,當(dāng)分區(qū)數(shù)為16時(shí),加速比為13.23,其并行效率為82.69%。當(dāng)分區(qū)數(shù)為32時(shí),并行的加速比為21.77,并行效率為68.03%。采用共享式和并行方案和分布式并行方案時(shí)的加速效率對(duì)比圖如圖8所示。
圖8 機(jī)翼翹曲并行方法的加速比
三維網(wǎng)格不同分區(qū)之間的拓?fù)溥B接關(guān)系也比二維情況也更加復(fù)雜,通信過(guò)程也更為復(fù)雜,導(dǎo)致額外的通信耗時(shí)。雖然和二維情形相比,VerBSS并行算法推廣到三維情形時(shí),相同的分區(qū)數(shù)目下并行效率低于二維情形的效率,但該并行算法依然能夠保持較高的加速效率,當(dāng)分區(qū)數(shù)目為8時(shí),并行效率仍然超過(guò)80%。
針對(duì)彈簧類網(wǎng)格變形方法求解較慢的特點(diǎn),本文針對(duì)基于彈簧類的VerBSS算法提出了一種分布式的VerBSS算法的并行方案。算例結(jié)果表明: 基于分區(qū)策略的VerBSS并行算法具有良好的并行效率,同時(shí)保持了Ball-vertex彈簧單元的抗壓能力,變形能力較強(qiáng),在二維情形和三維情形下都具備良好的適應(yīng)性。該算法適用于大變形和大規(guī)模的網(wǎng)格變形問(wèn)題。同時(shí),基于分布式的并行方案較共享式的并行架構(gòu)而言,更適合高性能計(jì)算,具有良好的可移植性。