一種分布式的點球彈簧修勻法并行方案

昌繼海，趙銘偉，余 飛，關(guān)振群*

(1.大連理工大學(xué)工程力學(xué)系，遼寧 大連 116024;2.大連理工大學(xué)電子信息與電氣工程學(xué)部，遼寧 大連 116024)

1 引言

在求解邊界運動等非定常流動問題時，針對網(wǎng)格變形問題的并行化研究是求解幾何變形問題的重要方向[1]。網(wǎng)格變形方法主要分為插值類方法和彈簧類方法[2]兩大類。

插值類的網(wǎng)格變形方法將流體內(nèi)部網(wǎng)格節(jié)點的位移視為邊界節(jié)點的插值問題，該類方法由于在變形過程中不考慮節(jié)點之間的連接性，直接通過邊界節(jié)點的位移對內(nèi)部節(jié)點進行插值，不需要迭代，容易實現(xiàn)并行化，許多研究者都對插值類方法的并行化進行了研究。Boer[3]等人提出的徑向基函數(shù)法(Radial Basis Function， RBF)是一種易于并行的插值方法，該方法通過選取特定的徑向基函數(shù)，利用邊界節(jié)點對內(nèi)部節(jié)點進行插值得到內(nèi)部節(jié)點的位移；Rendall[4]和曾[5]等人在RBF方法的基礎(chǔ)上進行了改進，如通過貪婪算法減少插值點的數(shù)目等等。Li[6]和Fang[7]等人提出了基于主從架構(gòu)的并行方法來提高RBF法的并行效率。此外，劉[8]等人提出的背景網(wǎng)格法也是一種易于并行的插值類方法，但是該方法面臨凹邊界的幾何問題或者推廣到三維問題時會遇到困難。

物理類方法主要包括Batina[9]等人提出的彈簧類方法，該方法將網(wǎng)格類比為一個由彈簧組成的彈性系統(tǒng)。相比于插值類變形方法，對彈簧類變形方法的并行化研究較少，一方面是由于彈簧類方法變形能力較弱；另一反面，彈簧類方法在計算過程中還要始終考慮節(jié)點之間的連接性，需要求解大規(guī)模的線性方程組，導(dǎo)致計算效率較低。

一些研究者提出了許多改進的彈簧類方法，如增加扭轉(zhuǎn)彈簧[10]或垂直彈簧[11]等方式來提高變形能力?；谶@些改進的彈簧類方法，Ma[12]和程[13]等人進行了分布式并行化的研究，取得了較好的并行效率。

林[14]等人在Ball-vertex彈簧法的基礎(chǔ)上提出一種局部插值的點球彈簧修勻法(VerBSS)方法，該方法通過分解彈簧系統(tǒng)的方式對網(wǎng)格進行逐點求解，具有插值類方法的特征，是一種求解效率較高的網(wǎng)格變形方法。

本文在林等人的基礎(chǔ)上，采用分區(qū)和共享單元層的策略，提出了一種分布式的VerBSS并行方案，既保留了該算法較高的變形能力，又實現(xiàn)了較高加速效率的并行化。最后，通過算例對本文提出的VerBSS并行算法和共享式并行方法的變形能力進行了對比，在二維和三維情形下對該算法的加速性能進行了測試。算例結(jié)果表明這種并行方案變形能力較強，具有較高的并行效率，適合推廣到大變形、大規(guī)模的網(wǎng)格變形問題。

2 基于彈簧模型的方法

2.1 Ball-vertex彈簧法

普通彈簧法的原理，是將網(wǎng)格的邊看作彈簧，其剛度與其邊的長度成反比。相鄰的節(jié)點i和節(jié)點j之間的力fij可表示為

(1)

ui和uj分別為節(jié)點i和節(jié)點j的位移，kij為ij彈簧的剛度，nij為節(jié)點i到節(jié)點j的單位向量。為避免這種現(xiàn)象的出現(xiàn)，通過在三角形節(jié)點與其對邊之間插入垂直彈簧來阻止單元塌陷，如圖1a)所示，p點為i節(jié)點到邊Ejk的垂線的垂足。這種垂直彈簧稱為Ball-vertex彈簧。

圖1 Ball-vertex彈簧法

Ball-vertex彈簧作用在節(jié)點i上的彈簧力可表示為

(2)

其中up和ui分別為節(jié)點p和節(jié)點i的位移，kip為ip彈簧的剛度，nip為節(jié)點i到節(jié)點p的單位向量。p點的位移可以通過節(jié)點j和節(jié)點k的位移插值獲得，為

up=ξuj+(1-ξ)uk

(3)

其中ξ為節(jié)點在邊Ejk上的插值系數(shù)。因此，ip彈簧作用在p點上的力可分解為作用在j點和k點上的分力。則節(jié)點i上的彈簧系統(tǒng)線性方程組為

(4)

其中K為節(jié)點i的剛度矩陣，u為子彈簧系統(tǒng)各節(jié)點的位移向量，b為i點的載荷向量。

2.2 VerBSS算法的流程

VerBSS方法的求解方式不同于普通的彈簧方法，該方法通過在每個內(nèi)部節(jié)點的閉包上形成一個以該節(jié)點為核心，以所有與該節(jié)點相鄰接的節(jié)點為外圍邊界的多邊形或多面體閉包，如圖1b)所示。

VerBSS算法：

輸入：初始網(wǎng)格M0

輸出：變形網(wǎng)格M1

Step1：導(dǎo)入初始網(wǎng)格，查找所有內(nèi)部節(jié)點，建立點及其周圍節(jié)點的拓撲關(guān)系。

Step2：為每個內(nèi)部節(jié)點構(gòu)建多面體包腔的子彈簧系統(tǒng)，并對剛度矩陣K進行LDLT分解，將分解后的L矩陣和D矩陣，以及轉(zhuǎn)換后的右端項存儲起來。

Step3：將已知邊界的節(jié)點位移u代入式(4)，遍歷網(wǎng)格內(nèi)部節(jié)點進行計算，并統(tǒng)計內(nèi)部節(jié)點位移量變化，作為收斂條件。

Step4：判斷Step3是否滿足收斂條件，若不收斂，重復(fù)step3，直至收斂。

3 VerBSS的分布式并行方案

3.1 VerBSS并行算法的通信方式

本文采用METIS分區(qū)工具將初始網(wǎng)格劃分為若干個分區(qū)[15]，進行任務(wù)劃分。每個子分區(qū)相應(yīng)的節(jié)點和單元的數(shù)據(jù)存儲在各自的存儲區(qū)域，通過拓寬相鄰分區(qū)邊界的方式來形成共享單元層，使得這些分區(qū)邊界上的節(jié)點被納入到求解序列中，如圖2a)和圖2b)所示為初始網(wǎng)格和分區(qū)后的網(wǎng)格。

圖2 二維機翼流場網(wǎng)格分區(qū)

選擇一個分區(qū)將其邊界上所有與邊界線相連接的單元的幾何信息傳遞給其鄰接分區(qū)，那么這些共享單元線將拓寬為一層共享單元面，相鄰的兩個分區(qū)共享這層單元面，共享單元層如圖2c)所示。三維情形中的共享單元面將拓寬為一層共享單元墻，相鄰的兩個分區(qū)共享這層單元墻。通過構(gòu)建相鄰分區(qū)之間的共享單元層，每個計算分區(qū)的節(jié)點都被分為兩類，一類為普通節(jié)點，一類為共享單元層上的節(jié)點。共享單元層上的節(jié)點可分為該分區(qū)的發(fā)送塊和接收塊。通過與其相鄰的分區(qū)建立對應(yīng)的發(fā)送和接收關(guān)系，對通信的發(fā)送地址、接收地址和通信數(shù)據(jù)的類型、數(shù)目以及其它通信的信息進行初始化，將相鄰分區(qū)的通信渠道建立起來。

具體的通信過程如圖3a)所示，共享單元層上的a節(jié)點和b節(jié)點同時為左分區(qū)和右分區(qū)所共享。左分區(qū)只計算a節(jié)點，不計算b節(jié)點，而右分區(qū)只計算b節(jié)點，不計算a節(jié)點。左分區(qū)在求得a節(jié)點的位移之后將a節(jié)點的位移傳遞給右分區(qū)，右分區(qū)的計算同理，如圖3b)所示。

圖3 共享單元層節(jié)點的閉包圖

3.2 VerBSS并行算法的流程

VerBSS的并行算法步驟如下：

輸入：初始分區(qū)網(wǎng)格

輸出：變形網(wǎng)格

Step1：初始化分區(qū)個數(shù)，每個網(wǎng)格分區(qū)分配一個進程讀取相應(yīng)的子分區(qū)網(wǎng)格。

Step2：查找相鄰分區(qū)的邊界，并在相鄰分區(qū)中間構(gòu)建共享區(qū)。

Step3：對每個進程讀取的子分區(qū)網(wǎng)格以及構(gòu)建的共享區(qū)網(wǎng)格進行分類，分為發(fā)送塊、接收塊和普通塊，并將每個進程和其它進程要通信的發(fā)送地址、接收地址和通信數(shù)據(jù)的類型、數(shù)目等通信信息進行初始化。

Step4：每個進程采用VerBSS算法對該分區(qū)進行求解，具體求解過程參見本文2.2節(jié)中的VerBSS算法流程。

Step5：在各相鄰分區(qū)之間進行通信，按照step2建立的通信關(guān)系分別進行發(fā)送和接收，每個進程將接收到的位移數(shù)據(jù)整理之后在對應(yīng)的節(jié)點上更新。

Step6：每個進程對更新的位移求位移平方差的平均值，判斷其是否滿足給定的收斂條件。若滿足，則停止循環(huán)；若不滿足，重復(fù)step4和step5，直到滿足收斂條件為止。

4 VerBSS并行算法的并行效率分析

本文算例將在一臺48核心、128G內(nèi)存的工作站上進行測試(CPU型號為Inter Xeon(R) E5-2650 v4，主頻2.2G)。

4.1 二維模型的算例分析

以一個二維機翼做俯沖運動的模型，其節(jié)點數(shù)為603 802，單元數(shù)為1 203 662。機翼長度為102個單位長度，單元尺寸為1個單位長度，單步旋轉(zhuǎn)角度為0.02。如圖4所示為變形前后的對比圖。

圖4 機翼俯沖變形前后的對比圖

采用分區(qū)并行的VerBSS算法單步迭代過程中各個步驟的平均耗時和加速比如表1所示。

表1 二維機翼旋轉(zhuǎn)分區(qū)并行的單步各步驟耗時

若采用共享式并行的VerBSS算法，單步迭代過程中各個步驟的平均耗時和加速比如表2所示。

表2 二維機翼旋轉(zhuǎn)的共享式并行的單步各步驟耗時

圖5為采用共享式并行方案和分布式并行方案時的加速效率對比圖。在單個變形時間步的耗時中，矩陣構(gòu)造時不需要通信，為完全解耦的過程，求解步驟和通信步驟交替進行，其中求解步驟為并行過程，通信步驟為串行過程。在該算例中，當分區(qū)數(shù)為16時，加速比為15.58，其并行效率為97.37%。當分區(qū)數(shù)為32時，并行的加速比為26.93，并行效率為84.15%。可見，由二維模型的算例可知，該并行算法在分區(qū)數(shù)目不超過16時，并行的加速比和理想加速比相當。當分區(qū)數(shù)目持續(xù)增大時，達到收斂條件所需的迭代步數(shù)也相應(yīng)增加，導(dǎo)致并行效率降低，但和共享式的并行方案相比仍然具有較高的加速效率。

圖5 機翼俯沖模型變形的加速比

對比普通彈簧法、串行VerBSS算法和本文的并行算法，三種方法變形時網(wǎng)格的最小夾角隨時間步的變化如圖6所示，可見在算法并行化的過程網(wǎng)格單元的質(zhì)量具備良好的魯棒性。

圖6 最小夾角隨時間步的變化圖

4.2 三維網(wǎng)格算例分析

進一步驗證本文算法在三維情形中的并行效率，一個三維的機翼在空間中做翹曲變形的網(wǎng)格模型如圖7所示。圖7a)為8分區(qū)的初始網(wǎng)格內(nèi)部截面圖，圖7b)為變形后的分區(qū)網(wǎng)格內(nèi)部截面圖，其彎曲形狀為拋物線型。

圖7 三維機翼翹曲變形的流場網(wǎng)格

機翼周圍的流場長寬高分別為1 500、500、500個單位長度，網(wǎng)格單元尺寸為0.5個單位長度，該流場網(wǎng)格包含有132 120個節(jié)點和656 878個單元，翼尖變形的單步步長為0.5個單位長度。采用分區(qū)并行的VerBSS算法迭代過程中各個步驟的平均耗時如表3所示。采用共享式并行的VerBSS算法各個步驟的耗時如表4所示。

表3 機翼翹曲分區(qū)并行的單步各步驟耗時

表4 機翼翹曲的共享式并行的單步各步驟耗時

該算例中，當分區(qū)數(shù)為16時，加速比為13.23，其并行效率為82.69%。當分區(qū)數(shù)為32時，并行的加速比為21.77，并行效率為68.03%。采用共享式和并行方案和分布式并行方案時的加速效率對比圖如圖8所示。

圖8 機翼翹曲并行方法的加速比

三維網(wǎng)格不同分區(qū)之間的拓撲連接關(guān)系也比二維情況也更加復(fù)雜，通信過程也更為復(fù)雜，導(dǎo)致額外的通信耗時。雖然和二維情形相比，VerBSS并行算法推廣到三維情形時，相同的分區(qū)數(shù)目下并行效率低于二維情形的效率，但該并行算法依然能夠保持較高的加速效率，當分區(qū)數(shù)目為8時，并行效率仍然超過80%。

5 結(jié)語

針對彈簧類網(wǎng)格變形方法求解較慢的特點，本文針對基于彈簧類的VerBSS算法提出了一種分布式的VerBSS算法的并行方案。算例結(jié)果表明： 基于分區(qū)策略的VerBSS并行算法具有良好的并行效率，同時保持了Ball-vertex彈簧單元的抗壓能力，變形能力較強，在二維情形和三維情形下都具備良好的適應(yīng)性。該算法適用于大變形和大規(guī)模的網(wǎng)格變形問題。同時，基于分布式的并行方案較共享式的并行架構(gòu)而言，更適合高性能計算，具有良好的可移植性。