基于局部通信的有限元分析并行算法優(yōu)化研究

吳建平，蔣 濤，彭 軍，銀?？?，楊錦輝

(1.國防科技大學(xué) 氣象海洋學(xué)院，長沙 410073； 2.國防科技大學(xué) 計算機(jī)學(xué)院，長沙 410073)

1 引 言

連續(xù)問題一般難以直接進(jìn)行求解，為此，通常采用有限差分、有限體積及有限元等方法進(jìn)行離散求解。有限元法將連續(xù)的求解區(qū)域離散為一組數(shù)量有限且按一定方式相互聯(lián)結(jié)的單元，典型的單元形式為多邊形，其頂點稱為結(jié)點。各單元之間可以存在面、邊或節(jié)點的重合，同時，單元本身可以采用不同形狀，并在每個單元內(nèi)采用假設(shè)的近似函數(shù)，以分片方式表示整體上的未知函數(shù)，未知函數(shù)及其導(dǎo)數(shù)在各個結(jié)點上的值作為求解時的未知量，即自由度，從而使一個連續(xù)問題轉(zhuǎn)化為自由度有限的離散問題。

由于有限元方法不受計算對象在幾何、物理與邊界條件上的限制，且計算精度高、使用靈活，因而自最初主要針對彈性力學(xué)與結(jié)構(gòu)分析以來，已逐步擴(kuò)展到航空、航天、造船、建筑、機(jī)械、海洋、氣象、水利、核能和地質(zhì)等眾多應(yīng)用領(lǐng)域。在采用有限元離散求解時，如果是采用時間依賴的微分方程建模的問題，且對時間采用隱式時間離散，或者是金屬成形與結(jié)構(gòu)疲勞分析與損傷等非線性與動態(tài)問題，則一般還需要求解線性方程組。因此，有限元離散的一般計算過程可以描述為

uk + 1=uk+(Kk)-1{Nk(uk)+pk}

(1)

式中uk為第k次迭代時的未知向量，pk為第k次迭代時的已知向量，Kk為第k次迭代時的總剛度矩陣，Nk為計算uk影響的線性或非線性算子，且

(2)

(3)

式中Ge為單元e的結(jié)點到全局結(jié)點的映射，如對單元e的第m個結(jié)點，其對應(yīng)的全局結(jié)點編號為n，則Ge第m行第n列上的元素為1，該行中其他元素全為0。

由式(1)可知，有限元分析的一般過程如下。

(3) 計算Nk(uk)，其一般通過先從uk得到每個單元對應(yīng)結(jié)點上的值，再通過單元上的計算，得到每個單元的貢獻(xiàn)，最后對貢獻(xiàn)進(jìn)行累加得到最終的總貢獻(xiàn)。

(4) 求解式(1)右端第二項所對應(yīng)的稀疏線性方程組，并計算得到uk + 1。

前三步為與單元對應(yīng)的計算，最后一步為與結(jié)點對應(yīng)的計算。因此，在具體實現(xiàn)時，采用的數(shù)據(jù)結(jié)構(gòu)可以分為單元與結(jié)點兩類，以適應(yīng)于每類計算各自的需要。在并行計算時，可以在進(jìn)行單元對應(yīng)計算時，按單元進(jìn)行數(shù)據(jù)分布，在進(jìn)行節(jié)點對應(yīng)計算時，按節(jié)點進(jìn)行數(shù)據(jù)分布，并在兩類計算之間進(jìn)行數(shù)據(jù)重分布，來進(jìn)行整體上的并行算法設(shè)計[1,2]。在此整體設(shè)計框架下，在并行計算時，有限元分析中與通信有關(guān)的核心算法如下。

(2) 剛度矩陣的并行裝配。該操作可以通過在本任務(wù)局部，先對其上所有單元完成剛度矩陣的裝配，之后再將每個任務(wù)上的結(jié)果裝配成總剛度矩陣?？倓偠染仃嚤仨毎葱蟹植嫉矫總€任務(wù)上，其分布方式對應(yīng)于結(jié)點在任務(wù)上的分配方式。

(3) 單元貢獻(xiàn)的并行裝配。每個任務(wù)上逐單元計算并通過對其上所有單元進(jìn)行累加，得到部分結(jié)果{Nk(uk)+pk}m，其中m為任務(wù)號。之后，對這些部分結(jié)果進(jìn)行累加，形成總的Nk(uk)+pk。由于{Nk(uk)+pk}m不僅含有分布在本任務(wù)上的分量，還可能含有分布到其他任務(wù)上的分量。因此，同一分量可能在多個任務(wù)上具有局部結(jié)果，這些局部結(jié)果需要進(jìn)行累加。

(4) 稀疏線性方程組的并行求解。由于對三維問題或者規(guī)模很大的問題，采用迭代法求解相對更為有效且容易控制，因此這里考慮迭代法。這種方法的并行計算可以從整體上按矩陣行數(shù)平均分布的方式來進(jìn)行，每個向量的分量也對應(yīng)地進(jìn)行平均分布，其中用到的內(nèi)積與范數(shù)在每個任務(wù)上均進(jìn)行存儲，采用多進(jìn)程規(guī)約從每個任務(wù)上的局部結(jié)果得到最終結(jié)果。一般在迭代中均會采用預(yù)條件技術(shù)，加快收斂速度[3,4]，因此，此操作中除內(nèi)積與范數(shù)計算外，預(yù)條件的應(yīng)用與稀疏矩陣稠密向量乘也可能需要通信。由于預(yù)條件種類繁多，這里不考慮其通信優(yōu)化問題，在實驗中也只采用簡單的分塊型并行預(yù)條件。按數(shù)據(jù)分布規(guī)則，求解所得向量在每個任務(wù)上也平均分布。

有限元分析的并行計算牽涉到諸多問題，吸引了很多學(xué)者進(jìn)行研究。剛度矩陣的裝配是有限元分析并行計算中最復(fù)雜的問題之一，Rezende等[5]對該問題給出了一種共享存儲并行算法，剛度矩陣按有限元網(wǎng)格中每個節(jié)點進(jìn)行分組裝配，每個處理器只裝配與特定節(jié)點組對應(yīng)的行，有效緩解了處理器對同一內(nèi)存單元的同時更新。Unterkircher等[6]提出了一種基于圖論的并行裝配算法與線性方程組求解算法，裝配算法不依賴于維數(shù)、單元類型與模型形狀，且在線性方程組求解器中引入了一種乘性對稱Schwarz預(yù)條件，并通過三維金屬擠壓過程模擬，驗證了所提出方法的有效性。Cecka等[7]針對NVIDIA GPU，提出了采用CUDA，并使用全局、共享和本地內(nèi)存來進(jìn)行并行裝配的多種策略，通過與串行裝配對比進(jìn)行了實驗驗證。Fabrice等[8]對地震波在非線性非彈性地質(zhì)中傳播的有限元數(shù)值模擬，采用基于混合并行的直接求解器Pastix，以及基于網(wǎng)格著色與MPI_Allgather進(jìn)行了并行裝配算法設(shè)計，并以法國里維埃拉地區(qū)尺度模型為對象進(jìn)行了實驗驗證。Jansson[9]對裝配過程采用對總剛度矩陣的每一行采用一個鏈表的數(shù)據(jù)結(jié)構(gòu)，并采用基于PGAS的單邊通信對每個處理器上所得局部裝配結(jié)果，通過遠(yuǎn)程復(fù)制來完成最終裝配。

Johan等[1]對計算流體力學(xué)有限元方法在CM-2與CM-200上的并行計算進(jìn)行了研究，利用預(yù)條件GMRES與數(shù)據(jù)結(jié)構(gòu)進(jìn)行了仔細(xì)考慮，并利用CM上提供的通信原語進(jìn)行了實現(xiàn)。Sonzogni等[10]針對機(jī)群系統(tǒng)，對PETSc-FEM庫進(jìn)行有限元并行計算的代碼實現(xiàn)進(jìn)行了研究，并對laplace方程與NS方程的實現(xiàn)性能進(jìn)行了測試分析。Kennedy等[11]對大規(guī)模基于梯度設(shè)計的結(jié)構(gòu)優(yōu)化問題，描述了一個并行有限元分析框架工具，并著重研究了其中稀疏線性方程組的并行直接求解算法。

Wu等[2]針對混凝土試件的細(xì)觀力學(xué)有限元分析，進(jìn)行了分布存儲并行算法框架整體設(shè)計，并針對其中剛度矩陣裝配、稀疏線性方程組求解及稀疏向量求和等核心問題，采用稀疏數(shù)據(jù)結(jié)構(gòu)及MPI_Alltoall及MPI_Allgather等進(jìn)行了并行算法的具體實現(xiàn)。本文主要以該工作為基礎(chǔ)，對其設(shè)計思想進(jìn)行梳理，并對其中核心算法存在的全局通信等問題，采用局部通信和計算通信重疊等技術(shù)，進(jìn)一步優(yōu)化有限元分析的并行計算效率。

2 剛度矩陣并行裝配算法的優(yōu)化

文獻(xiàn)[2]描述了一種對剛度矩陣進(jìn)行并行裝配的算法，其基本思想可以描述為,先每個任務(wù)按其上單元進(jìn)行局部剛度矩陣裝配，得到局部總剛度矩陣。之后，將各任務(wù)上的局部總剛度矩陣?yán)奂悠饋?，形成總剛度矩陣，并在每個任務(wù)上存儲其若干行，這兩步操作描述如圖1所示。最后，為便于進(jìn)行預(yù)條件構(gòu)造，再對各行中的元素按列號從小到大排序。

圖1 三個任務(wù)時的剛度矩陣并行裝配算法

可以看出，假設(shè)第k個任務(wù)上的局部矩陣記為Ak，Ak用CSR格式[3]存儲在valhbk,idxhbk和iptrhbk中。由于每個任務(wù)上最終只存儲矩陣中給定的若干行，記行數(shù)為nk，所以在任務(wù)k上將Ak(valhbk，idxhbk，iptrhbk)分成p組，記為Ak,j(valhbk,j，idxhbk,j，iptrhbk,j)，j=0～p-1，分別對應(yīng)于Ak的nj行。之后，需要在任務(wù)j上得到所有Ak,j(valhbk,j，idxhbk,j，iptrhbk,j)，k=0～p-1之和。為達(dá)到此目的，可以采用MPI_Reduce_scatter來實現(xiàn)，但需要注意到，這里每個Ak,j都是稀疏矩陣塊，采用CSR數(shù)據(jù)結(jié)構(gòu)valhbk,j,idxhbk,j和iptrhbk,j進(jìn)行存儲，因此，通過采用MPI_Alltoallv將所有valhbk,j,idxhbk,j和 iptrhbk,j收集到任務(wù)j上再進(jìn)行求和來實現(xiàn)。

在將所有valhbk,j,idxhbk,j和iptrhbk,j收集到任務(wù)j上時，假設(shè)用val,idx和iptr三個數(shù)組且采用CSR格式存儲,來一次存儲第k=0,1,…,p-1個任務(wù)上接收到的數(shù)據(jù)，其中val依次存儲非零元素的值，idx依次存儲這些非零元素的列號，iptr記錄每行第一個非零元存儲在val與idx的起始位置。圖2 給出了對3個任務(wù)時，第k個任務(wù)上將接收到的數(shù)據(jù)存儲到val、idx和iptr的示意圖，其中從每個任務(wù)接收到的矩陣塊均為nk行，因此，iptr的維數(shù)為3nk+1，其第一個元素表示矩陣塊A0,k在val與idx的起始位置，第nk+1個元素表示矩陣塊A1,k在val與idx的起始位置，第2nk+1個元素表示矩陣塊A2,k在val與idx的起始位置。本文需要在任務(wù)k上先確定第A0,k,A1,k,A2,k在val與idx的起始位置以及各Aj,k內(nèi)每行的起始位置，后者通過對idxhb進(jìn)行MPI_Alltoallv操作來實現(xiàn)。對前者，先在任務(wù)k上確定各Ak,j的總非零元個數(shù)lengsk,j，再通過MPI_Alltoall得到任務(wù)k上接收到的各Aj,k中非零元個數(shù)lengrk,j。

圖2 三個任務(wù)時任務(wù)k上對所接收矩陣塊進(jìn)行存儲

為對接收到的矩陣塊進(jìn)行求和，引入額外的整型數(shù)組cbpos，cpos和cptr來鏈接val，idx和iptr所記錄元素中所有行號相同者，其中cbpos(m)為第m行第一個非零元素在val與idx的位置，cpos(k)為第m行中當(dāng)前正在操作的非零元素所在位置，cptr(k)為與位置k對應(yīng)元素處于同一行的下一個元素所在位置。采用這種稀疏數(shù)據(jù)結(jié)構(gòu)，便于從val與idx中逐行提取非零元素。在對同一行中具有相同列號的元素進(jìn)行累加后，即可得到總剛度矩陣的該行數(shù)據(jù)。

在經(jīng)過以上操作后，實際上已形成總剛度矩陣，且每個任務(wù)上存儲了應(yīng)該存儲的矩陣塊，但現(xiàn)在每個矩陣行中所存儲的元素是無序的。為使得各行中元素按列號從小到大的順序排列，先將CSR格式的局部矩陣轉(zhuǎn)存為CSC格式，這樣局部矩陣就逐列進(jìn)行存儲，列號從小到大，同時，在同一列內(nèi)，各行上的元素也按行號從小到大的順序進(jìn)行存儲。之后，將CSC格式再次轉(zhuǎn)存為CSR格式，就得到了所需的按行分塊存儲的整體剛度矩陣。

值得注意的是，本文主要考慮離散網(wǎng)格不變且需要反復(fù)進(jìn)行迭代的問題，因此，在前面所述算法中，所有整型信息只需要在起始階段進(jìn)行處理，所得結(jié)果在后續(xù)迭代中直接使用。如無特別說明，這也適用于本文后續(xù)算法的描述。

前面所述剛度矩陣并行裝配算法在通信時以MPI_Alltoallv為基礎(chǔ)，這一方面將涉及到每個任務(wù)到另外每個任務(wù)之間的通信與同步，因此，所涉及的任務(wù)是全局性的，當(dāng)進(jìn)行大規(guī)模并行計算時，其個數(shù)很多。另一方面，這也不利于隱藏通信開銷。因此，這種方式必將對并行計算的可擴(kuò)展性形成不利影響。這里引入局部通信操作來替代MPI_Alltoallv的功能，基本思想是，在每個任務(wù)k上，依次判斷其上第j=0,1,…,p-1個矩陣塊的非零元個數(shù)，當(dāng)lengsk,j≠0且j≠k時，才啟動到任務(wù)j的發(fā)送操作；當(dāng)lengrk,j≠0且j≠k時，才啟動到從任務(wù)j接收數(shù)據(jù)的操作。這些發(fā)送與接收操作采用非阻塞式函數(shù)MPI_Isend與MPI_Irecv來實現(xiàn)，并在將valhbk,k，idxhbk,k和iptrhbk,k復(fù)制到val，idx和iptr相應(yīng)位置上后，再等待這些通信操作的完成，以實現(xiàn)本地操作與通信的重疊。

3 稀疏矩陣稠密向量并行乘的優(yōu)化

稀疏矩陣與稠密向量相乘在稀疏線性方程組的迭代解法中很重要，其計算時間與稀疏矩陣中的非零元個數(shù)成正比。按前文所述數(shù)據(jù)分布規(guī)則，在進(jìn)行稀疏矩陣A與稠密向量x的乘法，以得到結(jié)果y的過程中，每個任務(wù)上都只具有A的部分行，且存儲x和y的局部分量，文獻(xiàn)[2]給出了與之對應(yīng)的矩陣向量乘并行算法，其基本做法是將A分為p個行塊，x和y也都分成為p塊，其第k塊分別記為Ak，xk和yk，k=0,1,…,p-1，其分配到任務(wù)k上。這樣，在任務(wù)k上計算yk的數(shù)據(jù)相關(guān)性如圖3所示。

圖3 三個任務(wù)時的稀疏矩陣與稠密向量并行乘法

顯然，進(jìn)行yk的計算時，似乎需要用到其他所有任務(wù)上的x分量。但事實上，由于在計算y的某一個分量時，其等于A中對應(yīng)行與x的內(nèi)積，因此，只需要計算該行中每個非零元與x中相應(yīng)分量(標(biāo)號等于非零元列號者)的乘積之和。對任務(wù)k，只需要統(tǒng)計所存儲Ak中所有非零元素的列號，再分析哪些列號對應(yīng)的x分量不在任務(wù)k上，再獲取這些分量，這就是y=Ax的通信結(jié)構(gòu)。

具體地，先分析yk依賴于x的哪些分量，這些分量又處于哪些任務(wù)上。這些數(shù)據(jù)可以采用類似于稀疏矩陣的存儲方法存儲在idxr中，利用iptrr(j)記錄依賴于第j個任務(wù)上的首分量在idxr的位置，并利用lengr(j)記錄所依賴的任務(wù)j上的分量個數(shù)，其值等于iptrr(j+1)-iptrr(j)。如圖4所示，圖上部idxr與lengr的兩個下標(biāo)中，第一個為其所在任務(wù)的標(biāo)號，第二個為需接收分量對應(yīng)的源任務(wù)號。圖下部idxs與lengs的兩個下標(biāo)中，第一個為其所在任務(wù)號，第二個為需發(fā)送分量對應(yīng)的目的任務(wù)號，且括號中所列參數(shù)為其上參數(shù)的對應(yīng)數(shù)據(jù)，即lengsj,k=lengrk,j且idxsj,k=idxrk,j。因此，idxs與lengs信息可以利用mpi_alltoall來獲取。之后，利用lengs即可計算出需要發(fā)送給其他每個任務(wù)的分量在idxs的首地址iptrs。

圖4 三個任務(wù)時稀疏矩陣向量乘中分量對應(yīng)參數(shù)交換

在計算y=Ax之前，需先收集本任務(wù)要用到的x分量，可以先在每個任務(wù)上利用idxs逐一提取需發(fā)送出去的分量并存儲到vals中，之后通過mpi_alltoallv操作，在每個任務(wù)上將接收到的數(shù)據(jù)存儲到valr，再將valr的元素存儲到臨時數(shù)組w中，使得w(idxr(j))=wr(j)。這樣，就可以直接計算yi=Aiw。

與剛度矩陣的裝配類似，在進(jìn)行y=Ax的實際計算時，這里也可以采用局部非阻塞式通信操作替代MPI_Alltoallv。對任務(wù)k，依次判斷其上lengsk,j(j=0,1,…,p-1)是否等于0，當(dāng)其等于0且j≠k時，才啟動到任務(wù)j的發(fā)送操作；之后，再依次判斷其上的lengrk,j(j=0,1,…,p-1)，當(dāng)其不等于0且j≠k時，才啟動到從任務(wù)j接收數(shù)據(jù)的操作。這些操作分別采用非阻塞式函數(shù)MPI_Isend與MPI_Irecv實現(xiàn)。

此外，采用文獻(xiàn)[12]對矩陣向量乘中分量進(jìn)行分類的方法，將每個任務(wù)上yk的分量分為內(nèi)部分量和邊界分量兩類。內(nèi)部分量的計算不依賴于其他任務(wù)上的x分量，僅依賴于本任務(wù)局部分配到的x分量，無需通信即可進(jìn)行計算。邊界分量至少依賴于其他任務(wù)上的一個x分量，需要在接收到所有依賴的x分量之后才能完成計算?？梢詫⑦@種對分量分類的方法與局部非阻塞式通信相結(jié)合，進(jìn)行計算與通信重疊，減少通信開銷的影響?；舅枷胧牵谌蝿?wù)k上，先對所有l(wèi)engsk,j不為0且不等于k的j，啟動到任務(wù)j的非阻塞式發(fā)送操作。對所有l(wèi)engrk,j不為0且不等于k的j，啟動從任務(wù)j的非阻塞式接收操作。其次，進(jìn)行yk內(nèi)部分量的計算。在非阻塞式接收操作完成之后，再進(jìn)行邊界分量的計算。這樣，內(nèi)部分量的計算就能與通信相互重疊，從而一定程度上實現(xiàn)對通信開銷的隱藏。

4 結(jié)點分量的局地化

采用MPI_Allgatherv進(jìn)行編程實現(xiàn)的方法很簡單，其最終使得每個任務(wù)上都得到x的所有分量，但這帶來了一些不必要的額外通信。事實上，每個單元在具體計算時，只需要用到其對應(yīng)結(jié)點的分量，因此，對任務(wù)k上的所有單元進(jìn)行分析，可以獲知該任務(wù)在進(jìn)行其上單元對應(yīng)計算時實際需要依賴的分量，這一般并不會覆蓋到x的全部分量，本文基于該事實，對結(jié)點分量的局地化進(jìn)行優(yōu)化。

將任務(wù)k上各單元在進(jìn)行計算時需要依賴的結(jié)點分量分為p組，即xk,0,xk,1,…,xk,p -1。類似于第3節(jié)，將這些數(shù)據(jù)存儲到valrk與idxrk，利用iptrrk,j記錄xk,j首分量在valrk與idxrk的位置，并利用lengrk,j記錄xk,j的分量個數(shù)，即iptrrk,j + 1- iptrrk,j。圖5給出了3個任務(wù)時的示意圖，顯然，此時任務(wù)0上各單元的計算不會用到任務(wù)2上的結(jié)點分量，即lengr0,2=0，因此，無需從任務(wù)2接收數(shù)據(jù)。

為在各任務(wù)上獲取其所含單元計算時所需數(shù)據(jù)，先對lengrk,j通過MPI_Alltoall操作，得到lengsk,j=lengrj,k，即需要從任務(wù)k發(fā)送給任務(wù)j的分量個數(shù)。之后，在任務(wù)k上，先對所有l(wèi)engrk,j不為0且不等于k的j，將idxrk,j發(fā)送給任務(wù)j。對所有l(wèi)engsk,j不為0且不等于k的j，從任務(wù)j接收數(shù)據(jù)并存儲到idxsk,j。當(dāng)進(jìn)行具體結(jié)點分量的局地化時，以圖5為例，在每個任務(wù)上對應(yīng)的操作可以概述如下。

任務(wù)0需從任務(wù)1接收x由idxr0,1指明的分量；將x由idxs0,1=idxr1,0指明的分量發(fā)給任務(wù)1。

任務(wù)1需從任務(wù)0和2分別接收x由idxr1,0和idxr1,2指明的分量；將x的由idxs1,0=idxr0,1和idxs1,2=idxr2,1指明的分量分別發(fā)給任務(wù)0和2。

任務(wù)2需從任務(wù)1接收x由idxr2,1指明的分量；將x由idxs2,1=idxr1,2指明的分量發(fā)給任務(wù)1。

圖5 三個任務(wù)上結(jié)點分量的局地化

表1給出與圖5對應(yīng)的一個實例，從每行看，對應(yīng)的是每個任務(wù)需要從其他任務(wù)接收的x分量之指標(biāo)，表中對角線位置所示為按數(shù)據(jù)分布規(guī)則分配給每個任務(wù)的分量。此時，局地化操作對應(yīng)的實際上就是要按表中各列所示，從對角線處指明分量中提取各任務(wù)所需的分量。

表1 三個任務(wù)時對應(yīng)于圖5的結(jié)點分量局地化實例

一般地，在任務(wù)k上，可以先提取x按idxsk,j所指明的分量，依次存儲到valsk,j中。其次，對所有l(wèi)engsk,j不為0且不等于k的j，啟動到任務(wù)j的非阻塞式發(fā)送操作，發(fā)送valsk,j。對所有l(wèi)engrk,j不為0且不等于k的j，啟動從任務(wù)j的非阻塞式接收操作，接收valrk,j。之后，逐步等待接收操作的完成，每完成一個接收操作，則將對應(yīng)的valrk,j按照idxrk,j指明的位置復(fù)制到x中。這樣，不僅通過局部通信減少了所涉及的任務(wù)個數(shù)與通信量，而且通過非阻塞式通信，一定程度上實現(xiàn)了通信與局地數(shù)據(jù)復(fù)制操作的重疊，從而也有利于對通信開銷的隱藏。此時，每個任務(wù)進(jìn)行所屬單元的計算時，需要用到的x分量已全部得到。

5 單元貢獻(xiàn)并行裝配的優(yōu)化

在有限元分析中，在計算單元的貢獻(xiàn)時，需要對每個任務(wù)上的貢獻(xiàn)進(jìn)行累加，由于單元幾乎平均地分配到各個任務(wù)，所以在各個任務(wù)上得到的局部向量是稀疏的。文獻(xiàn)[2]給出了一種對有限元分析中單元貢獻(xiàn)按稀疏向量進(jìn)行并行裝配的簡單方法，其基本思想是，假設(shè)第k個任務(wù)上的各單元總貢獻(xiàn)形成的局部向量為xk，用稀疏格式存儲在valsk和idxsk中，之后，通過MPI_Allgatherv可以使得每個任務(wù)上都得到所有valsk和idxsk，再同時在每個任務(wù)上將接收到的valsk和idxsk疊加到本地貢獻(xiàn)上，得到x=x0+x1+…+xp -1。

上述方法雖然思想簡單，但并非每個任務(wù)均需要所有x的分量，實際上只需要得到按數(shù)據(jù)分布規(guī)則應(yīng)該計算的分量即可，其他分量并不需要計算。以表1為例，此時對應(yīng)的操作即是在各列中，對所有分量求和，并將求和結(jié)果存儲在對角線所示任務(wù)上。本文基于該思想給出一個優(yōu)化算法，其基本思想是，在任務(wù)k上，將valsk和idxsk各分成p段，分別記為valsk,j和idxsk,j(j=0～p-1)，使得idxsk,j每個指標(biāo)對應(yīng)分量是分布到第j個任務(wù)。顯然，該操作需要的通信過程可以看成是第4節(jié)通信過程的逆過程。以圖5為例，此時，在每個任務(wù)上對應(yīng)的操作可以概述如下。

任務(wù)0需從任務(wù)1接收x由idxs0,1=idxr1,0指明的分量，并將其與本地x求和；需將x由idxr0,1指明的分量發(fā)給任務(wù)1。

任務(wù)1需從任務(wù)0和2分別接收由idxs1,0=idxr0,1和idxs1,2=idxr2,1所指分量，并將其與本地x求和；將x由idxr1,0和idxr1,2指明的分量分別發(fā)給任務(wù)0和2。

任務(wù)2需從任務(wù)1接收x由idxs2,1=idxr1,2指明的分量，并將其與本地x求和；需將x由idxr2,1指明的分量發(fā)給任務(wù)1。

一般地，在任務(wù)k上，可以先提取x按idxrk,j所指明的分量，依次存儲到valrk,j中。其次，對所有l(wèi)engrk,j不為0且不等于k的j，啟動到任務(wù)j的非阻塞式發(fā)送操作，發(fā)送valrk,j。對所有l(wèi)engsk,j不為0且不等于k的j，啟動從任務(wù)j的非阻塞式接收操作，接收valsk,j。之后，逐步等待接收操作的完成，每完成一個接收操作，則將對應(yīng)的valsk,j按照idxsk,j指明的位置加到x上。這樣，既通過局部通信減少了所涉及的任務(wù)個數(shù)與通信量，又通過非阻塞式通信實現(xiàn)了通信與局地求和操作的重疊，從而也有利于對通信開銷的隱藏。

6 混凝土細(xì)觀數(shù)值模擬的實驗驗證

細(xì)觀力學(xué)分析與研究是當(dāng)前研究混凝土材料特性的熱門方法之一，該方法將混凝土看作由粗骨料、硬化水泥膠體以及兩者之間的界面粘結(jié)帶組成的三相非均質(zhì)復(fù)合材料，數(shù)值模擬是進(jìn)行該型材料研究的一種重要手段。為對混凝土材料進(jìn)行細(xì)觀力學(xué)數(shù)值模擬，對形成的材料可以采用多種模型，這里采用隨機(jī)骨料隨機(jī)參數(shù)模型[15]，并采用有限元法進(jìn)行離散求解，首先在細(xì)觀層次上對試件進(jìn)行有限元剖分，考慮骨料單元、固化水泥砂漿單元及界面單元材料力學(xué)特性的不同，再用簡單的破壞準(zhǔn)則或損傷模型反映單元剛度的退化，之后模擬試件的裂縫擴(kuò)展過程及破壞形態(tài)。

本文針對2個混凝土試件的靜態(tài)加載，在對有限單元進(jìn)行任務(wù)劃分時，將大致相同數(shù)量的連續(xù)元素分配給每個處理器。在對對應(yīng)于結(jié)點的位移變量進(jìn)行任務(wù)劃分時，采用基于連續(xù)編號、基于離散網(wǎng)格和基于METIS等三種任務(wù)劃分方式，實驗中分別記為CONT，MESH和METIS。 對于CONT，大致相同數(shù)量的連續(xù)變量分配給每個處理器。對于MESH，變量根據(jù)相應(yīng)的坐標(biāo)進(jìn)行分區(qū)，以便每個部分的八個部分盡可能是方形的。對METIS，變量通過metis-4.0的接口metis_PartGraphRecursive來進(jìn)行劃分。

所有的實驗都是在64個節(jié)點的集群上進(jìn)行的，每個節(jié)點有2個Intel(R)Xeon(R)CPU E5-2692 v2@2.20 GHz(緩存30720 kb)和64 G內(nèi)存。操作系統(tǒng)為2.6.32-431.17.1-aftms-th，編譯器為Intel Fortran版本15.0.0.090，節(jié)點間采用Infiniband互連，消息傳遞接口為MPICH 3.2.1。在所有的實驗中，所遇到的稀疏線性系統(tǒng)都是對稱正定的。因此，全部使用共軛梯度迭代。初始解為零向量，當(dāng)殘差向量歐幾里德范數(shù)與初始向量歐幾里德范數(shù)之比小于 1e -6 時停止迭代。為了加快收斂速度，使用文獻(xiàn)[13]提供的預(yù)條件(即對稱正定矩陣的ILUT版本)的塊Jacobi并行化。

實驗針對2個混凝土試件進(jìn)行，其分別為三級配大試件和兩級配濕篩試件。三級配試件的尺寸為300 mm×300 mm×1100 mm，總共離散為53200個單元，44117個結(jié)點，記為模型1。濕篩試件的尺寸為150 mm×150 mm×550 mm，總共離散為78800個單元，71013個網(wǎng)格點，記為模型2。關(guān)于兩試件的具體描述可以參見文獻(xiàn)[14,15]，靜態(tài)加載時的有限元分析計算流程可以參見文獻(xiàn)[15]。每步加載時荷載增量為0.25 kn。無損傷單元時，剛度矩陣與前一步相同，不需要重新裝配，先前所得預(yù)條件進(jìn)行重用。對某些單元在某些加載步時出現(xiàn)退化的非線性問題，可以通過求解多個線性方程組來解決。但在每一步的計算過程中，剛度矩陣和預(yù)條件保持不變，并進(jìn)行重用。對模型1，在第439步出現(xiàn)損壞單元，在第567步時試件完全損壞，共求解696個線性方程組。對模型2，第59步出現(xiàn)損壞單元，第94步時試件完全損壞，共求解178個線性方程組。

表2列出了剛度矩陣并行裝配所需要的時間，由表2可知，在任務(wù)個數(shù)不超過128的情況下，優(yōu)化幾乎起不到效果，優(yōu)化后并行裝配所耗費的時間一般反而更長，只有對試件1，在任務(wù)個數(shù)較多且采用CONT時，耗費時間得以減少。這主要是因為兩方面的原因所致,一是在任務(wù)個數(shù)較少時，將全局通信操作分解為局部通信操作帶來的通信所涉任務(wù)個數(shù)減少不明顯;二是局地復(fù)制操作耗時很短，對通信開銷的隱藏不明顯。但是，在任務(wù)總數(shù)不斷增加的情況下，每個任務(wù)進(jìn)行通信實際牽涉到的其他任務(wù)個數(shù)變化不大，因而，在任務(wù)個數(shù)很多時，將全局通信操作分解為局部通信操作，可以有效減少每個任務(wù)在通信操作中所涉及的其他任務(wù)個數(shù)。因此，在任務(wù)個數(shù)不斷增加的情況下，相比優(yōu)化之前的差異程度，優(yōu)化后有減小的趨勢，由此可以預(yù)計，在采用更多任務(wù)時，優(yōu)化算法存在減少執(zhí)行時間的潛力。

表2 混凝土試件靜載實驗中剛度矩陣并行裝配所用總時間(單位：s)

表3列出了稀疏矩陣稠密向量并行乘所用的時間，由表3可知，無論對哪種劃分、哪種試件或哪種任務(wù)規(guī)模，采用優(yōu)化技術(shù)時，稀疏矩陣稠密向量并行乘所用時間均得以顯著減少。同時，隨著總?cè)蝿?wù)個數(shù)的增加，優(yōu)化后并行乘所用時間與優(yōu)化前相比，總體上有不斷減小的趨勢，這說明總體上優(yōu)化效果在任務(wù)個數(shù)較多時更明顯。同時，對相同試件在相同任務(wù)個數(shù)時的模擬，相對其他兩種劃分算法，采用CONT算法所用時間較長。

表3 混凝土試件靜載實驗中稀疏矩陣稠密向量乘所用總時間(單位：s)

表4給出了節(jié)點分量局地化時所用的總時間，由表4可知，相比于優(yōu)化前，優(yōu)化后節(jié)點分量局地化所用時間大幅減少。在采用優(yōu)化算法之前，所用時間隨任務(wù)個數(shù)的增加，呈不斷增加的趨勢，而在采用優(yōu)化算法之后，這種趨勢得到了明顯扭轉(zhuǎn)，隨著任務(wù)個數(shù)的增加，呈現(xiàn)出一定程度的減少趨勢。

表4 混凝土試件靜載實驗中結(jié)點分量局地化所用總時間(單位：s)

表5給出了結(jié)點分量局地化所用的時間，由 表5 可知，相比于優(yōu)化前，優(yōu)化后單元貢獻(xiàn)并行裝配所用時間大幅減少。在采用優(yōu)化算法之前，所用時間隨任務(wù)個數(shù)的增加呈成倍增加的趨勢，而在采用優(yōu)化算法之后，這種趨勢得到了明顯扭轉(zhuǎn)，隨著任務(wù)個數(shù)的增加，呈現(xiàn)優(yōu)于對折減少的趨勢。

表6給出了表2～表5所列剛度矩陣并行裝配、稀疏矩陣稠密向量并行乘、節(jié)點分量局地化及單元貢獻(xiàn)并行裝配所用時間之和，由表6可知，無論是對試件1還是試件2，對三種劃分中的任意一種，當(dāng)采用相同個數(shù)的任務(wù)時，優(yōu)化后這4者的執(zhí)行時間之和都明顯減少。需要注意的是，對剛度矩陣并行裝配，如表2所示，雖然優(yōu)化算法在任務(wù)個數(shù)更多時具有潛在優(yōu)勢，但在所測試的任務(wù)個數(shù)下，優(yōu)化后幾乎沒什么改進(jìn)，甚至有時反而不如優(yōu)化之前，但由于這一過程所用總時間相對很小，因此，對總時間沒有明顯影響。

表5 混凝土試件靜載實驗中單元貢獻(xiàn)并行裝配所用總時間(單位：s)

表6 混凝土試件靜載實驗中表2～表5所列的時間之和(單位：s)

表7給出了對每個試件，從設(shè)置階段到加載全部完成整個過程所用時間的情況。由表7可知，無論對哪種試件、哪種劃分及哪種任務(wù)數(shù)，采用優(yōu)化算法極大地減少了整個數(shù)值模擬的時間。特別是在任務(wù)個數(shù)較多時，改進(jìn)尤為明顯。同時可見，在不采用優(yōu)化算法時，對單元個數(shù)與結(jié)點個數(shù)都較少的試件1，采用64個任務(wù)時，相對32個任務(wù)已經(jīng)幾乎不能減少執(zhí)行時間。對單元個數(shù)與結(jié)點個數(shù)稍多的試件2，雖然執(zhí)行時間稍有減少，但也微乎其微。但采用優(yōu)化算法時，不僅總時間相對優(yōu)化前有明顯減少，而且在采用64個任務(wù)時，相對32個任務(wù)，大部分情況下，執(zhí)行時間得以進(jìn)一步減少。對試件1，在采用CONT時執(zhí)行時間能得到進(jìn)一步減少。對試件2，在采用METIS與MESH時，相對32個任務(wù)，采用64個任務(wù)執(zhí)行時間都得以進(jìn)一步減少，說明采用優(yōu)化算法還改善了整個數(shù)值模擬的可擴(kuò)展性。

表7 混凝土試件靜載實驗所用總時間(單位：s)

7 結(jié) 論

本文針對有限元分析的計算問題，對其所涉核心算法通過局部通信與計算通信重疊的方式進(jìn)行了優(yōu)化設(shè)計，并通過混凝土試件細(xì)觀力學(xué)靜態(tài)加載數(shù)值模擬進(jìn)行了實驗驗證，結(jié)果表明，無論是采用哪種任務(wù)劃分，優(yōu)化后相比現(xiàn)有算法數(shù)值模擬所用時間都得以顯著減小。同時，實驗驗證中發(fā)現(xiàn)，每個任務(wù)在進(jìn)行通信時，所涉及的其他任務(wù)個數(shù)都相對較多，對METIS與MESH尤其如此。這在采用局部通信操作來實現(xiàn)時，必然引起從同一個任務(wù)同時向多個任務(wù)發(fā)送數(shù)據(jù)，以及同一個任務(wù)同時從多個任務(wù)接收數(shù)據(jù)的通信沖突問題，如何解決該問題，對進(jìn)一步采用計算通信重疊減小通信開銷的影響、提高并行執(zhí)行效率具有很大影響，這是下一步進(jìn)行重點研究的問題。