付朝江,王天奇,林悅榮

(1.福建工程學(xué)院 土木工程學(xué)院,福州 350108; 2.福建省土木工程新技術(shù)與信息化重點(diǎn)實(shí)驗(yàn)室,福州 350108)(*通信作者電子郵箱cjfu@163.com)

0 引言

有限元法進(jìn)行結(jié)構(gòu)非線性動(dòng)力反應(yīng)的計(jì)算非常耗時(shí),并行計(jì)算可提供行之有效的方法。顯式時(shí)間離散可避免大量矩陣的存儲(chǔ)和重復(fù)求解線性方程系統(tǒng),但受條件穩(wěn)定性限制,較小的時(shí)間步長導(dǎo)致模擬需要很多時(shí)間步數(shù),采用顯式有限元結(jié)構(gòu)動(dòng)力分析的并行計(jì)算方法可以在較短的時(shí)間內(nèi)求解大規(guī)模問題。

利用有限元法求解結(jié)構(gòu)動(dòng)力問題時(shí),最廣泛采用的顯式時(shí)間積分技術(shù)為中心差分算法,通常為顯式Newmark-β法,可表示為:

(1)

(2)

(3)

(4)

(5)

算法1 中心差分算法。

時(shí)間步循環(huán)

節(jié)點(diǎn)/自由度循環(huán)

1)更新位移和速度向量(式(1)和(2))

2)計(jì)算有效節(jié)點(diǎn)力(式(3)和(4))

3)計(jì)算加速度向量(式(5))

4)計(jì)算穩(wěn)定時(shí)間步

結(jié)束節(jié)點(diǎn)循環(huán)

5)更新時(shí)間t=t+Δt

結(jié)束時(shí)間步循環(huán)

因?yàn)榈?)步涉及材料應(yīng)力更新的單元應(yīng)力計(jì)算，所以第1)和3)步的計(jì)算成本相對(duì)第2)步小。采取區(qū)域分解,各個(gè)計(jì)算步均可進(jìn)行并行化。

算法中的每個(gè)時(shí)間步構(gòu)成一個(gè)自然的同步點(diǎn),為了在t+Δt時(shí)繼續(xù)進(jìn)行,所有節(jié)點(diǎn)和單元在t時(shí)的量需為已知。在下一個(gè)全局同步點(diǎn),負(fù)載平衡使所有的處理器同時(shí)完成各自的工作任務(wù)。在顯式結(jié)構(gòu)動(dòng)力有限元分析中,CPU時(shí)間與單元數(shù)成比例。顯然,顯式非線性動(dòng)力分析的并行實(shí)現(xiàn)主要挑戰(zhàn)在于負(fù)載平衡。

有限元計(jì)算的并行數(shù)值算法大都采取區(qū)域分解技術(shù),將有限元網(wǎng)格劃分為子網(wǎng)格,每個(gè)子網(wǎng)格由單個(gè)處理器獨(dú)立處理。通過采取有效的靜態(tài)網(wǎng)格劃分技術(shù)來實(shí)現(xiàn)處理器間的負(fù)載平衡和通信最小化[1-2]。

動(dòng)態(tài)負(fù)載平衡成為近幾年來并行計(jì)算研究的重要領(lǐng)域,已有的動(dòng)態(tài)負(fù)載平衡算法主要目的是設(shè)計(jì)各種方法來減小計(jì)算過程中出現(xiàn)的負(fù)載不平衡問題,即應(yīng)用自適應(yīng)重網(wǎng)格化技術(shù)[3-4]。這些算法的研究主要是在同構(gòu)并行計(jì)算機(jī)上開展，針對(duì)網(wǎng)絡(luò)機(jī)群環(huán)境的分布式多指令多數(shù)據(jù)流(Multiple Instruction Multiple Data stream, MIMD)并行計(jì)算系統(tǒng)開展的研究甚少。大部分動(dòng)態(tài)負(fù)載平衡算法應(yīng)用于以下兩步[5-7]:第1步,算法確定工作負(fù)載的數(shù)量,在鄰近的處理器中進(jìn)行分配以維持各處理器中負(fù)載均衡,如Cybenko[5]提出的擴(kuò)散算法和Boillat[7]的確定工作流算法。第2步是任務(wù)選定,即算法選擇有限元網(wǎng)格的單元或節(jié)點(diǎn)，該步利用處理器之間單元或節(jié)點(diǎn)的交換來改進(jìn)現(xiàn)有的分解。隨后,將所選定的單元或節(jié)點(diǎn)在鄰近的處理器中進(jìn)行任務(wù)遷移。動(dòng)態(tài)負(fù)載平衡算法確定的分解質(zhì)量很大程度取決于所采用的選擇過程[9]。

已有的動(dòng)態(tài)負(fù)載平衡算法是針對(duì)具體的應(yīng)用來處理計(jì)算過程中出現(xiàn)的計(jì)算不平衡而開展研究,并在專用的同構(gòu)并行計(jì)算機(jī)上進(jìn)行，很少有算法是針對(duì)網(wǎng)絡(luò)機(jī)群環(huán)境的分布式MIMD并行計(jì)算系統(tǒng)開展[10]。本文在網(wǎng)絡(luò)機(jī)群并行計(jì)算環(huán)境下，采取區(qū)域分解、動(dòng)態(tài)任務(wù)分配及計(jì)算與通信重疊的優(yōu)化通信方法，開展非重疊通信區(qū)域分解并行算法、重疊通信區(qū)域分解并行算法、群動(dòng)態(tài)任務(wù)分配算法、動(dòng)態(tài)任務(wù)分配算法及動(dòng)態(tài)負(fù)載平衡算法的研究，探討基于有效并行求解策略的有限元并行算法的有效性，并對(duì)算法的性能進(jìn)行了評(píng)價(jià)分析。

1 并行求解策略

上述中心差分算法可采用重疊區(qū)域分解或非重疊區(qū)域分解技術(shù)實(shí)現(xiàn)。采用重疊區(qū)域分解技術(shù)時(shí),在每個(gè)處理器中位于網(wǎng)格劃分的界面單元被復(fù)制,節(jié)點(diǎn)僅屬于各自區(qū)域。這樣,這些共享單元的內(nèi)力計(jì)算需要在每個(gè)處理器中復(fù)制,該計(jì)算工作非常耗時(shí)。相反,采用非重疊區(qū)域分解技術(shù)時(shí),單元僅屬于各自區(qū)域,界面節(jié)點(diǎn)為各區(qū)域共享,如圖1所示。因而,界面單元力計(jì)算不需復(fù)制,而只是涉及界面節(jié)點(diǎn)位移、速度和加速度向量的計(jì)算需要復(fù)制。在這兩種方法中,處理器間的通信階數(shù)(即重疊區(qū)域分解中的加速度向量和非重疊區(qū)域分解中的內(nèi)力向量)是相同的。因此,針對(duì)顯式非線性動(dòng)力分析的并行實(shí)現(xiàn),采取非重疊區(qū)域分解技術(shù)的方法較為優(yōu)越。該非重疊區(qū)域分解中心差分時(shí)間步并行算法簡稱為非重疊通信區(qū)域分解并行算法(non-overl apped算法)可用算法2描述。

圖1 非重疊區(qū)域分解的網(wǎng)格劃分

算法2 非重疊通信區(qū)域分解并行算法。

時(shí)間步循環(huán)

節(jié)點(diǎn)/自由度循環(huán)

1)更新位移和速度向量(式(1)和(2))

2)計(jì)算有效節(jié)點(diǎn)力(式(3)和(4))

3)通過消息傳遞交換共享節(jié)點(diǎn)的力

4)計(jì)算加速度向量(式(5))

5)計(jì)算穩(wěn)定時(shí)間步

結(jié)束節(jié)點(diǎn)循環(huán)

6)更新時(shí)間t=t+Δt

結(jié)束時(shí)間步循環(huán)

1.1 計(jì)算與通信重疊

在算法2描述的并行算法中, 通過將通信和計(jì)算重疊，實(shí)現(xiàn)通信延遲的最小化，以進(jìn)一步改善算法的性能。該并行算法中,在子區(qū)域中每個(gè)單元的內(nèi)力向量計(jì)算完成之后,通信才開始。然而,如果由共享節(jié)點(diǎn)構(gòu)成的共享單元的內(nèi)力計(jì)算完成,得到所需的信息(共享節(jié)點(diǎn)的內(nèi)力向量)同相鄰節(jié)點(diǎn)進(jìn)行交換,此時(shí)開始通信而不是等待內(nèi)力計(jì)算結(jié)束才開始通信,其余單元的計(jì)算可以獨(dú)立進(jìn)行,從而實(shí)現(xiàn)通信和計(jì)算重疊。通信和計(jì)算重疊可通過將單元力計(jì)算的任務(wù)劃分為兩部分：第一部分計(jì)算共享節(jié)點(diǎn)的內(nèi)力向量,該階段通信發(fā)送和第二部分的計(jì)算同時(shí)進(jìn)行；第二部分計(jì)算內(nèi)部單元的內(nèi)力向量。該優(yōu)化通信的非重疊區(qū)域分解中心差分時(shí)間步并行算法簡稱為重疊通信區(qū)域分解并行算法(overlapped算法)可用算法3描述。

算法3 重疊通信區(qū)域分解并行算法。

時(shí)間步循環(huán)

節(jié)點(diǎn)/自由度循環(huán)

1)更新位移和速度向量(式(1)和(2))

2)計(jì)算有效節(jié)點(diǎn)力(式(3)和(4))

3)設(shè)立非阻塞通信以交換共享節(jié)點(diǎn)的力,同時(shí)進(jìn)行內(nèi)部節(jié)點(diǎn)的有效節(jié)點(diǎn)力的計(jì)算

4)等待處理器間通信完成

5)計(jì)算加速度向量(式(5))

6)計(jì)算穩(wěn)定時(shí)間步

結(jié)束節(jié)點(diǎn)循環(huán)

7)更新時(shí)間t=t+Δt

結(jié)束時(shí)間步循環(huán)

1.2 動(dòng)態(tài)任務(wù)分配

現(xiàn)有的顯式非線性動(dòng)力分析的并行算法是假定整個(gè)計(jì)算過程中計(jì)算狀態(tài)保持不變,據(jù)此假設(shè),計(jì)算負(fù)載按單元或節(jié)點(diǎn)數(shù)相等簡單地分配給各個(gè)處理器，利用圖劃分算法即可實(shí)現(xiàn)[1]；然而,這種劃分不適應(yīng)材料非線性。通常塑性局限于某些單元,而這些單元需要按彈塑性情況處理,可能需要對(duì)單元每個(gè)的積分點(diǎn)按很復(fù)雜的材料準(zhǔn)則進(jìn)行處理；而其他不考慮塑性的單元,涉及的計(jì)算相對(duì)簡單、消耗的計(jì)算時(shí)間較少,因此，并行求解中會(huì)出現(xiàn)負(fù)載不平衡。本節(jié)對(duì)非線性動(dòng)力分析的動(dòng)態(tài)任務(wù)分配方法進(jìn)行探討。

在中心差分法中,動(dòng)態(tài)任務(wù)分配的并行化只在計(jì)算階段考慮[11],可利用主-從(Master-Slave)模式實(shí)現(xiàn)。Master(即主處理器)作為整個(gè)時(shí)間積分過程的協(xié)調(diào)者并分配任務(wù)給Slave(從處理器);Slave從Master接收任務(wù),進(jìn)行必要的計(jì)算并將計(jì)算結(jié)果發(fā)送給Master,等待Master的下一步數(shù)據(jù)。這里無需進(jìn)行有限元網(wǎng)格劃分,位移計(jì)算在Master中進(jìn)行,每個(gè)單元的內(nèi)力計(jì)算分配給Slave。所有的Slave被分配任務(wù)后,Master等待Slave返回的結(jié)果。一旦接收到某個(gè)Slave的結(jié)果數(shù)據(jù)包,一個(gè)新的數(shù)據(jù)包(包含涉及另一個(gè)網(wǎng)格單元的數(shù)據(jù))就會(huì)發(fā)送給Slave。這樣,任務(wù)動(dòng)態(tài)分配給所有處理器,計(jì)算負(fù)載將會(huì)很好地平衡，然而在動(dòng)態(tài)任務(wù)分配中的同步通信會(huì)較高。動(dòng)態(tài)任務(wù)分配(Dynamic Task Allocation, DTA)算法如算法4。

算法4 動(dòng)態(tài)任務(wù)分配(DTA)算法。

主處理器:

時(shí)間步循環(huán)

1)更新位移和速度向量(式(1)和(2))

2)除去所有單元標(biāo)記

3)發(fā)送無標(biāo)記的單元數(shù)據(jù)給空閑的從處理器

4)標(biāo)記信息包已發(fā)送的單元

5)對(duì)單元進(jìn)行循環(huán)

(a)接收來自從處理器的單元內(nèi)力向量

(b)發(fā)送相應(yīng)于無標(biāo)記的單元信息包給從處理器

(c)標(biāo)記相應(yīng)的單元

(d)組集內(nèi)部節(jié)點(diǎn)力向量

結(jié)束單元循環(huán)

6)計(jì)算有效節(jié)點(diǎn)力向量(式(3))

7)計(jì)算加速度向量(式(5))

8)計(jì)算穩(wěn)定時(shí)間步

9)更新時(shí)間t=t+Δt

結(jié)束時(shí)間步循環(huán)

從處理器:

當(dāng)循環(huán)沒有結(jié)束時(shí)

接收來自主處理器的單元信息包

計(jì)算單元內(nèi)力向量

發(fā)送單元內(nèi)力向量給主處理器

結(jié)束循環(huán)

1.3 群動(dòng)態(tài)任務(wù)分配

由DTA算法可知,從處理器需要頻繁地和主處理器進(jìn)行通信,從而增加了算法的通信復(fù)雜度并隨網(wǎng)格中單元數(shù)線性增加。為減小處理器間的通信量,提出群算法。群算法形成多個(gè)群,每個(gè)群由一組單元構(gòu)成,每個(gè)群中的單元數(shù)不必相同。然而,群數(shù)比處理器數(shù)多。DTA中群形成算法(DTA with clustering算法)如算法5。

算法5 DTA中群形成算法。

1)除去所有的單元標(biāo)記

2)計(jì)算網(wǎng)格中每個(gè)單元的度

//度表示連接到該單元的無標(biāo)記單元的數(shù)目

3)對(duì)無標(biāo)記的單元循環(huán)計(jì)算

選擇具有最高度的無標(biāo)記單元

循環(huán)計(jì)算

如果單元的度為零,則

(a)將單元分配給相鄰的群

否則

(a)設(shè)置群數(shù)為:cluster_conter=cluster_conter+1

(b)標(biāo)記單元并選擇與該單元相鄰的所有單元來組成一個(gè)群

(c)標(biāo)記所有選擇的單元

(d)重算所有無標(biāo)記的單元的度

結(jié)束循環(huán)

結(jié)束循環(huán)

DTA算法是對(duì)群操作而不是單元。每次涉及每個(gè)群的信息發(fā)送給空閑的從處理器,從處理器計(jì)算群中所有單元內(nèi)力并組集以形成該群的整體內(nèi)力向量,將其發(fā)送給主處理器并等待下一組群信息;主處理器依次組集由從處理器接收的內(nèi)力向量以形成整個(gè)網(wǎng)格的整體內(nèi)力向量,并進(jìn)一步形成有效的節(jié)點(diǎn)力向量和計(jì)算該時(shí)間步的節(jié)點(diǎn)加速度。群DTA的優(yōu)點(diǎn)就是通信復(fù)雜度減小了NE/NC,其中：NC為形成的群數(shù)、NE為網(wǎng)格中的單元數(shù)。附加的開銷是群的形成,該部分可在預(yù)處理階段完成。

1.4 動(dòng)態(tài)負(fù)載平衡

工作站機(jī)群環(huán)境下有限元計(jì)算的負(fù)載平衡算法可按下述任一方法[12]設(shè)計(jì)：

方法一 將網(wǎng)格劃分成多個(gè)子網(wǎng)格,子網(wǎng)格數(shù)少于處理器數(shù)。首先,每個(gè)子網(wǎng)格分配給各個(gè)處理器。隨著計(jì)算進(jìn)行和負(fù)載不平衡出現(xiàn),調(diào)整處理器數(shù)以減小負(fù)載不平衡。

方法二 將網(wǎng)格劃分成子網(wǎng)格,子網(wǎng)格數(shù)大于處理器數(shù)。將一組子網(wǎng)格分配給每個(gè)處理器。隨著計(jì)算進(jìn)行和負(fù)載不平衡出現(xiàn),在處理器間調(diào)動(dòng)子網(wǎng)格以減小負(fù)載不平衡。

方法三 將網(wǎng)格劃分成子網(wǎng)格,子網(wǎng)格數(shù)等于處理器數(shù)。每個(gè)子網(wǎng)格分配給一個(gè)處理器。隨著計(jì)算進(jìn)行和負(fù)載不平衡出現(xiàn),通過在子網(wǎng)格間轉(zhuǎn)移單元來改進(jìn)初始的劃分以減小負(fù)載不平衡。

方法一為一個(gè)多級(jí)方法,會(huì)導(dǎo)致多個(gè)處理器在同一時(shí)間通信,導(dǎo)致網(wǎng)絡(luò)工作站的性能大大降低。方法二不會(huì)產(chǎn)生這種通信方式；然而,單元塊在處理器間的移動(dòng)不能保證切割邊緣的最小化。方法三不會(huì)產(chǎn)生復(fù)雜的通信方式，因?yàn)閱卧w移代替了單元塊移動(dòng),切割邊緣的最小化可以實(shí)施。由于網(wǎng)絡(luò)工作站并行計(jì)算要求處理器間的通信最小化以充分利用處理器資源,設(shè)計(jì)動(dòng)態(tài)負(fù)載平衡算法時(shí),方法三最優(yōu)。

大多數(shù)動(dòng)態(tài)負(fù)載平衡算法采用擴(kuò)散算法[4]，但這些算法需要較小的迭代次數(shù)收斂,明顯地增加涉及動(dòng)態(tài)負(fù)載平衡的開銷。本文采用方法三實(shí)現(xiàn)負(fù)載動(dòng)態(tài)平衡。動(dòng)態(tài)負(fù)載平衡算法構(gòu)架的設(shè)計(jì)可分為下述5個(gè)階段。

1.4.1 評(píng)價(jià)處理器能力

該階段動(dòng)態(tài)負(fù)載平衡器收集數(shù)據(jù)。將處理器能力CPi定義為每秒內(nèi)處理器能完成的節(jié)點(diǎn)/單元數(shù)的計(jì)算指令。在并行算法執(zhí)行期間,如果一個(gè)處理器同時(shí)正在執(zhí)行各種其他工作,處理器能力隨時(shí)間變化。處理器能力確定是通過將工作負(fù)載除以處理器利用時(shí)間。工作負(fù)載Wi為處理器i的單元/節(jié)點(diǎn)數(shù),處理器利用時(shí)間Ti為處理器i的運(yùn)行真實(shí)時(shí)間，處理器能力CPi可表示為CPi=Wi/Ti。

1.4.2 檢查負(fù)載不平衡

該階段負(fù)載平衡器檢查是否需要重新平衡計(jì)算負(fù)載。為便于判斷,采用負(fù)載不平衡系數(shù)來評(píng)估[13],該系數(shù)定義為:LIF=Tmax/Tideal-1。其中：Tmax為處理器利用時(shí)間的最大值；Tideal為總工作負(fù)載和總處理能力的比值,即表示當(dāng)工作負(fù)載完全平衡和處理器是有效利用時(shí)的理想時(shí)間。

(6)

其中:NP為處理器數(shù),如果總工作負(fù)載是完全平衡,LIF為零;否則將會(huì)大于零。當(dāng)負(fù)載不平衡系數(shù)LIF超過容許值時(shí),將進(jìn)行負(fù)載平衡初始化。

1.4.3 形成負(fù)載轉(zhuǎn)移矩陣

這個(gè)階段基于計(jì)算的負(fù)載不平衡,計(jì)算在處理器間將要轉(zhuǎn)移的工作量,即轉(zhuǎn)移的單元數(shù)，可通過形成負(fù)載轉(zhuǎn)移矩陣來實(shí)現(xiàn)。

依據(jù)LIF,要求負(fù)載平衡,需要形成負(fù)載轉(zhuǎn)移矩陣F,其大小為NP×NP[9]。負(fù)載轉(zhuǎn)移矩陣的每個(gè)元素Fij確定在處理器i和j間進(jìn)行轉(zhuǎn)移的工作量,以實(shí)現(xiàn)計(jì)算負(fù)載平衡。本文負(fù)載轉(zhuǎn)移矩陣F不需顯式形成,而是計(jì)算負(fù)載轉(zhuǎn)移參數(shù),該參數(shù)表示將要發(fā)送或從相鄰處理器接收的單元數(shù)。負(fù)載轉(zhuǎn)移參數(shù)表示為:

ψi=(Ti-Tideal) ×CPi

(7)

ψi為正表明單元數(shù)從子圖i需要發(fā)送給相鄰子網(wǎng)格;類似地,ψi為負(fù)表明從相鄰子網(wǎng)格接收單元數(shù)。根據(jù)ψi值,可容易地構(gòu)成矩陣Fij。為識(shí)別和標(biāo)記每個(gè)子網(wǎng)格i中遷移的單元,本文采用ψ設(shè)計(jì)遷移算法。

1.4.4 選擇遷移單元

根據(jù)計(jì)算的ψi值,選擇相應(yīng)子網(wǎng)格i的單元遷移到子網(wǎng)格j;反之亦然。采用啟發(fā)算法實(shí)現(xiàn)這種遷移[14],如算法6所述。將遷移算法應(yīng)用于相應(yīng)子圖并改變子網(wǎng)格i和子網(wǎng)格j間的頂點(diǎn)。當(dāng)轉(zhuǎn)移子圖的邊界頂點(diǎn)給它相鄰的指定子圖時(shí),遷移算法可最小化邊緣切割。一旦遷移頂點(diǎn)確定,它對(duì)應(yīng)的單元也就標(biāo)定。

算法6 遷移算法。

1)按ψi值,確定主動(dòng)的(ψi為正)、被動(dòng)的(ψi為負(fù))和中性的(ψi為零)子圖。

2)從主動(dòng)子圖收集一組頂點(diǎn),該頂點(diǎn)為界面頂點(diǎn)或子圖不相連部分的頂點(diǎn)或按相同優(yōu)先次序的其他頂點(diǎn)。

3)標(biāo)記所有頂點(diǎn)。

4)從主動(dòng)子圖列表的組中提取一個(gè)標(biāo)記的頂點(diǎn)。

5)通過從主動(dòng)子圖i到相鄰被動(dòng)子圖j交換頂點(diǎn),計(jì)算邊緣切割的減小ECij。

6)計(jì)算所有相鄰被動(dòng)子圖的最大ECij和標(biāo)記相應(yīng)的頂點(diǎn)及被動(dòng)子圖。

7)重復(fù)第4)～7)步，直到一組中所有頂點(diǎn)標(biāo)記完并更新最大ECij值和相應(yīng)的數(shù)據(jù)結(jié)構(gòu)。

8)交換標(biāo)記的頂點(diǎn)給被動(dòng)子圖j和更新界面頂點(diǎn)數(shù)據(jù)結(jié)構(gòu)。

9)更新負(fù)載平衡參數(shù)ψi和ψj。

10)從主動(dòng)頂點(diǎn)組中移動(dòng)標(biāo)記頂點(diǎn)。

11)標(biāo)記組中所有頂點(diǎn)。

12)重復(fù)第4)～11)步,直到平衡參數(shù)ψi為零。

13)重復(fù)第1)～12)步,直到所有子圖的平衡參數(shù)ψ為零。

1.4.5 遷移數(shù)據(jù)

一旦指定為遷移到接收處理器的單元被確定,這些單元連同它們的數(shù)據(jù)一起轉(zhuǎn)移。為了最小化遷移開銷,所有數(shù)據(jù)結(jié)構(gòu)須變換為一個(gè)流。這個(gè)流發(fā)送到接收處理器,它負(fù)責(zé)將接收的數(shù)據(jù)結(jié)構(gòu)恢復(fù)到初始狀態(tài)。

算法7 動(dòng)態(tài)負(fù)載平衡算法。

對(duì)時(shí)間步循環(huán)

對(duì)節(jié)點(diǎn)/自由度循環(huán)

1)更新位移和速度向量(式(1)和式(2))

2)計(jì)算共享節(jié)點(diǎn)的有效節(jié)點(diǎn)力(式(3))

3)建立非阻塞通信調(diào)用來交換共享節(jié)點(diǎn)力和同時(shí)進(jìn)行計(jì)算內(nèi)部點(diǎn)的有效節(jié)點(diǎn)力

4)等待處理器間通信結(jié)束

5)計(jì)算加速度向量(式(4))

6)計(jì)算穩(wěn)定的時(shí)間步

結(jié)束節(jié)點(diǎn)循環(huán)

If(需要檢查負(fù)載不平衡)

檢查負(fù)載不平衡

If(負(fù)載不平衡為真)

計(jì)算數(shù)據(jù)遷移信息

選取將要遷移的單元

遷移選擇的單元

Endif

Endif

7)更新時(shí)間:t=t+Δt

結(jié)束時(shí)間步循環(huán)

本文提出的動(dòng)態(tài)負(fù)載平衡算法應(yīng)用于結(jié)構(gòu)非線性動(dòng)態(tài)分析的并行計(jì)算，在動(dòng)態(tài)分析中，顯式動(dòng)態(tài)分析算法是進(jìn)行時(shí)間積分，當(dāng)負(fù)載平衡在通信后進(jìn)行初始化,負(fù)載評(píng)估發(fā)生在計(jì)算階段。通常,在動(dòng)態(tài)分析中,在每個(gè)時(shí)間步終點(diǎn)之后負(fù)載平衡可初始化[15]。然而,在每個(gè)時(shí)間步終點(diǎn)重新平衡常常導(dǎo)致太多的開銷并可能引起整體計(jì)算過程速度下降。因此,負(fù)載平衡啟動(dòng)通常每隔n個(gè)時(shí)間步進(jìn)行一次,n的最佳值取決于負(fù)載平衡的成本和負(fù)載平衡進(jìn)展的方式。

2 數(shù)值算例

本文提出的非線性動(dòng)力分析并行算法的數(shù)值計(jì)算在DELL工作站機(jī)群上進(jìn)行,該機(jī)群由4臺(tái)雙CPU的DELL工作站通過100.0 Mb/s以太網(wǎng)連接而成,共有8個(gè)CPU(2.4 GHz Xeon chips,512 KB cache),每個(gè)節(jié)點(diǎn)內(nèi)存1.0 GB。每臺(tái)工作站都有真實(shí)的IP地址,使用消息傳遞接口(MPI)[16-17]。操作系統(tǒng)為Red Hat Linux9.0。

圖2(a)為1/4的圓筒形殼,L=1 000 mm,兩直邊固定。受垂直均布荷載p=2.0 kN/m2,施加方式如圖2(b)。殼體厚度2 mm,E=2.06×105MPa,泊松比v=0.3,屈服應(yīng)力fy=235 MPa,質(zhì)量密度為7.8×103kg/m3。進(jìn)行圓筒形殼的非線性分析時(shí)，考慮其幾何非線性。

圖2 計(jì)算模型

將殼劃分為12×24結(jié)構(gòu)網(wǎng)格,采用8節(jié)點(diǎn)的殼單元。通過采用時(shí)間步為0.2 ms來驗(yàn)證算法的穩(wěn)定性和精度。用隱式積分算法(Newmark算法)[18]和本文基于區(qū)域分解的顯式非線性動(dòng)力時(shí)間積分算法(算法4)計(jì)算節(jié)點(diǎn)A位移時(shí)間歷程反應(yīng),計(jì)算結(jié)果如圖3所示?？梢钥闯霰疚牟⑿酗@式時(shí)間積分算法的結(jié)果與Newmark算法的結(jié)果非常吻合,表明本文并行顯式時(shí)間積分算法是有效的,即具有很好的穩(wěn)定性和精度。

圖3 節(jié)點(diǎn)A位移時(shí)程反應(yīng)

為確定動(dòng)態(tài)負(fù)載平衡(DLB)算法中負(fù)載不平衡系數(shù)LIF的容許值,設(shè)定每隔25個(gè)時(shí)間步啟動(dòng)一次負(fù)載平衡,LIF的容許值分別設(shè)為0.05、0.25、0.5時(shí),進(jìn)行算例計(jì)算(自由度為4 710),DLB算法計(jì)算時(shí)間如表1所示?？梢钥闯?容許值為0.05、0.5時(shí)計(jì)算所需的時(shí)間明顯增加,算法的性能下降。由于LIF取為0.05時(shí),該取值小,計(jì)算過程中出現(xiàn)負(fù)載不平衡的次數(shù)多,一旦出現(xiàn)負(fù)載不平衡需要進(jìn)行負(fù)載平衡初始化,數(shù)據(jù)重分配,從而增加了通信開銷;而LIF取為0.5時(shí),由于該取值較大,計(jì)算過程中出現(xiàn)負(fù)載不平衡次數(shù)減少,但由于LIF取值較大,一旦出現(xiàn)負(fù)載不平衡時(shí),負(fù)載不平衡程度高,需要進(jìn)行大量的數(shù)據(jù)重分配,數(shù)據(jù)遷移產(chǎn)生的通信開銷也較大。計(jì)算表明LIF的容許值為0.25時(shí)較優(yōu)。

取LIF的容許值為0.25,負(fù)載平衡啟動(dòng)每隔時(shí)間步n分別設(shè)為15、25、35,進(jìn)行算例計(jì)算(自由度為4 710),DLB算法計(jì)算時(shí)間如表2所示。由表2可見,啟動(dòng)間隔時(shí)間步為15時(shí),因啟動(dòng)負(fù)載不平衡檢查次數(shù)較多而增加時(shí)間；而啟動(dòng)間隔時(shí)間步為35時(shí),因步數(shù)較大,出現(xiàn)負(fù)載不平衡時(shí),負(fù)載不平衡程度較高,進(jìn)行數(shù)據(jù)重分配時(shí)數(shù)據(jù)遷移產(chǎn)生的通信開銷較大,致使計(jì)算時(shí)間比啟動(dòng)間隔時(shí)間步為25時(shí)稍有增加。故本文算例采用啟動(dòng)間隔時(shí)間步為25,LIF的容許值為0.25。

表1 不同 LIF 容許值的動(dòng)態(tài)負(fù)載平衡計(jì)算時(shí)間

表2 不同啟動(dòng)間隔時(shí)間步 n 的動(dòng)態(tài)負(fù)載平衡計(jì)算時(shí)間

為評(píng)估算法的性能加密有限元網(wǎng)格以增加問題大小(自由度)，為保證解的精度采用嚴(yán)格的收斂準(zhǔn)則，分別進(jìn)行非重疊通信區(qū)域分解并行算法、重疊通信區(qū)域分解并行算法、群動(dòng)態(tài)任務(wù)分配算法、動(dòng)態(tài)任務(wù)分配算法及動(dòng)態(tài)負(fù)載平衡算法的并行實(shí)現(xiàn)。為評(píng)價(jià)五種算法的并行性能,分別用1、2、4、6、8個(gè)處理器對(duì)所求問題進(jìn)行求解,用并行加速比對(duì)結(jié)果進(jìn)行性能分析,并行加速比定義為:speedup=用1個(gè)處理器的運(yùn)行時(shí)間/用n個(gè)處理器的運(yùn)行時(shí)間。

五種并行算法的加速比如圖4～5所示?？梢钥闯?隨著問題大小的增加,相應(yīng)算法的加速比提高?；趧?dòng)態(tài)任務(wù)分配的并行算法的性能并不理想，盡管動(dòng)態(tài)任務(wù)分配算法中考慮了動(dòng)態(tài)負(fù)載平衡,其加速比較小，結(jié)果并非滿意。而群動(dòng)態(tài)任務(wù)分配算法性能優(yōu)于動(dòng)態(tài)任務(wù)分配算法,但次于區(qū)域分解算法。這是由于動(dòng)態(tài)任務(wù)分配算法中處理器間的通信開銷很高,因而性能下降。在每個(gè)時(shí)間積分步,動(dòng)態(tài)任務(wù)分配算法中主和從處理器間的通信總數(shù)等于網(wǎng)格中的單元數(shù),因此,通信復(fù)雜度隨網(wǎng)格中的單元數(shù)而增加。再者,所有從處理器和主處理器不斷地進(jìn)行通信,導(dǎo)致較高的通信延遲，因而,這些通信開銷隨著處理器數(shù)的增加而增加。

群動(dòng)態(tài)任務(wù)分配算法減少了動(dòng)態(tài)任務(wù)分配時(shí)處理器間的通信開銷;然而,相對(duì)于優(yōu)化的區(qū)域分解算法,通信開銷仍然較高。群動(dòng)態(tài)任務(wù)分配算法的性能比動(dòng)態(tài)任務(wù)分配算法有所提高,這是由于處理器間的通信減少。由于每個(gè)從處理器處理一組單元,從處理器和主處理器間的通信會(huì)錯(cuò)開,這樣導(dǎo)致較小的通信延遲。區(qū)域分解算法中將部分通信和計(jì)算重疊較非重疊通信和計(jì)算時(shí),其性能得以提高。

動(dòng)態(tài)負(fù)載平衡算法不會(huì)產(chǎn)生復(fù)雜的通信方式,由于單元遷移代替單元塊移動(dòng),切割邊緣的最小化可以實(shí)施,處理器間的通信開銷較少,動(dòng)態(tài)負(fù)載平衡算法最優(yōu)。

本文算法與Newmark算法的CPU時(shí)間(自由度為4 710)如表3所示。可以看出,對(duì)相同規(guī)模的問題本文算法比隱式算法快,相應(yīng)的加速比高;并且算法的計(jì)算性能隨著問題大小的增加而提高,顯示本文并行顯式時(shí)間積分算法具有良好的性能。

圖4 算法的加速比(自由度為4 710)

圖5 算法的加速比(自由度為17 322)

處理器數(shù)CPU耗時(shí)/s本文算法Newmark算法1845.90987.652563.93705.464313.30395.068179.98235.15

3 結(jié)語

對(duì)顯式非線性動(dòng)力有限元分析的并行求解算法進(jìn)行研究。算法采用了非重疊區(qū)域分解技術(shù)、動(dòng)態(tài)任務(wù)分配和動(dòng)態(tài)負(fù)載平衡,通過將通信和計(jì)算重疊,對(duì)區(qū)域分解技術(shù)中的通信進(jìn)行優(yōu)化。數(shù)值算例表明優(yōu)化通信的區(qū)域分解算法的性能優(yōu)于傳統(tǒng)的區(qū)域分解算法;同時(shí)也表明區(qū)域分解技術(shù)遠(yuǎn)優(yōu)于基于動(dòng)態(tài)任務(wù)分配的算法?；趧?dòng)態(tài)任務(wù)分配算法的性能主要受處理器間通信復(fù)雜性的影響,通過采取群方法來降低動(dòng)態(tài)任務(wù)分配的通信復(fù)雜度。動(dòng)態(tài)負(fù)載平衡算法由于單元遷移代替單元塊移動(dòng),可以實(shí)施切割邊緣的最小化,處理器間的通信開銷較少,動(dòng)態(tài)負(fù)載平衡算法最優(yōu)。數(shù)值研究表明群動(dòng)態(tài)任務(wù)分配算法的性能高于動(dòng)態(tài)任務(wù)分配算法,群動(dòng)態(tài)任務(wù)分配算法的性能低于區(qū)域分解算法的性能,動(dòng)態(tài)負(fù)載平衡算法最優(yōu)。

參考文獻(xiàn)(References)

[1] WALSHAW C, CROSS M. Parallel optimisation algorithms for multilevel mesh partitioning[J]. Parallel Computing, 2000, 26(12): 1635-1660.

[2] DEVINE K, FLAHERTY J. Parallel adaptive hp-refinement techniques for conservation laws[J]. Applied Numerical Mathematics, 1996, 20(4): 367-386.

[3] ZHAO B, LIU Y, GOH S H. Parallel finite element analysis of seismic soil structure interaction using a PC cluster[J]. Computers and Geotechnics, 2016, 80: 167-177.

[4] HOU Y R, DU G Z. An expandable local and parallel two-grid finite element scheme[J]. Computers and Mathematics with Applications, 2016, 71(12): 2541-2556.

[5] CYBENKO G. Dynamic load balancing for distributed memory multiprocessors[J]. Journal of Parallel and Distributed Computing, 1989, 7(2): 279-301.

[7] BOILLAT J E. Load balancing and Poisson equation in a graph[J]. Concurrency and Computation: Practice and Experience, 1990, 2(4): 289-313.

[8] DIEKMANN R, FROMMER A, MONIEN B. Efficient schemes for nearest neighbor load balancing[J]. Parallel Computing, 1999, 25(7): 789-812.

[9] SZIVERI J, TOPPING BHV. Transient dynamic nonlinear analysis using MIMD computer architectures[J]. Journal of Computing in Civil Engineering, 2000, 14(2): 79-91.

[10] HU Y F, BLAKE R J. An improved diffusion algorithm for dynamic load balancing[J]. Parallel Computing, 1999, 25(4): 417-444.

[11] TOUHEED N, SELWOOD P, JIMACK P K. A comparison of some dynamic load-balancing algorithms for a parallel adaptive flow solver[J]. Parallel Computing, 2000, 26(12): 1535-1554.

[12] KRYSL P, BITTNAR Z. Parallel explicit finite element solid dynamics with domain decomposition and message passing: dual partitioning scalability[J]. Computers and Structures, 2001, 79(3): 345-360.

[13] RAO A R M. Explicit nonlinear dynamic finite element analysis on homogeneous/heterogeneous parallel computing environment[J]. Advances in Engineering Software, 2006, 37(11): 701-720.

[14] RAO A R M, RAO T, DATTAGURU B. Generating optimised partitions for parallel finite element computations employing float-encoded genetic algorithms[J]. Computer Modeling in Engineering and Sciences, 2004, 5(3): 213-234.

[15] KWAK J Y, CHUN T Y, SHIN S J, et al. Domain decomposition approach to flexible multibody dynamics simulation[J]. Computational Mechanics, 2014, 53(1): 147-158.

[16] 付朝江. 集群MPI環(huán)境下有限元結(jié)構(gòu)分析并行計(jì)算研究[D]. 上海: 上海大學(xué), 2006.(FU C J. The research on parallel computation of finite element structural analysis based on MPI cluster[D]. Shanghai: Shanghai University, 2006.)

[17] CAPODAGLIO G, AULISA E. A particle tracking algorithm for parallel finite element applications[J]. Computers and Fluids, 2017, 159: 338-355.

[18] 付朝江. 隱式非線性動(dòng)力分析有限元并行求解格式[J]. 工程力學(xué), 2010, 27(10): 27-33.(FU C J. Parallel algorithm for implicit nonlinear dynamic finite element analysis[J]. Engineering Mechanics, 2010, 27(10): 27-33.)

This work is partially supported by the National Natural Science Foundation of Chian (51378124).