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

昌繼海，趙銘偉，余 飛，關振群*

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

1 引言

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

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

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

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

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

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

2 基于彈簧模型的方法

2.1 Ball-vertex彈簧法

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

(1)

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

圖1 Ball-vertex彈簧法

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

(2)

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

up=ξuj+(1-ξ)uk

(3)

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

(4)

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

2.2 VerBSS算法的流程

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

VerBSS算法：

輸入：初始網格M0

輸出：變形網格M1

Step1：導入初始網格，查找所有內部節點，建立點及其周圍節點的拓撲關系。

Step2：為每個內部節點構建多面體包腔的子彈簧系統，并對剛度矩陣K進行LDLT分解，將分解后的L矩陣和D矩陣，以及轉換后的右端項存儲起來。

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

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

3 VerBSS的分布式并行方案

3.1 VerBSS并行算法的通信方式

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

圖2 二維機翼流場網格分區

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

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

圖3 共享單元層節點的閉包圖

3.2 VerBSS并行算法的流程

VerBSS的并行算法步驟如下：

輸入：初始分區網格

輸出：變形網格

Step1：初始化分區個數，每個網格分區分配一個進程讀取相應的子分區網格。

Step2：查找相鄰分區的邊界，并在相鄰分區中間構建共享區。

Step3：對每個進程讀取的子分區網格以及構建的共享區網格進行分類，分為發送塊、接收塊和普通塊，并將每個進程和其它進程要通信的發送地址、接收地址和通信數據的類型、數目等通信信息進行初始化。

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

Step5：在各相鄰分區之間進行通信，按照step2建立的通信關系分別進行發送和接收，每個進程將接收到的位移數據整理之后在對應的節點上更新。

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

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

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

4.1 二維模型的算例分析

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

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

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

表1 二維機翼旋轉分區并行的單步各步驟耗時

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

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

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

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

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

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

4.2 三維網格算例分析

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

圖7 三維機翼翹曲變形的流場網格

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

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

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

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

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

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

5 結語

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