基于CUDA的有限元矩陣并行裝配算法研究

胡斌星，李新國，孫 鵬

(1.西北工業(yè)大學 航天學院，西安 710072；2.西北工業(yè)大學 空天飛行技術研究所,西安 710072)

1 引 言

隨著航天運載器廣泛采用輕質復合材料和升力體的氣動布局等因素，其剛體和彈性體的模態(tài)耦合問題愈發(fā)突出[1]。傳統(tǒng)的方法是將箭體視為柔性梁分析[2,3]，將助推器也視為一維梁與芯級組成空間梁系，基于特定的小偏差增量模型分析彈性振動對控制系統(tǒng)的影響[4]。而在高超聲速飛行器的分析中也常以均勻分布質量的自由梁[5]或懸臂梁[6]建模來分析其與氣動和推進系統(tǒng)的耦合關系。隨著現(xiàn)代航天理論與相關行業(yè)技術的推進，航天運載器技術從一次性運載火箭向可重復使用運載器過渡，同時對新一代飛行器提出重復使用、健康監(jiān)測、故障診斷和任務重構的要求。將彈性問題從三維空間簡化為平面問題雖克服了計算的困難，在工程上也得到了驗證，但在很多情況下，這種簡化將會導致較大的誤差。如飛行器在飛行過程中產(chǎn)生滾轉時，需考慮結構的扭轉彎曲作用；簡化的非均質梁在載荷不均勻時，按平面柔性梁進行簡化誤差較大[7]；簡單的一維梁模型解析解無法模擬運載器各部位的應力應變以實現(xiàn)健康監(jiān)測，也無法通過傳感器所測的結構損傷程度計算當前極限過載以重構制導指令，無法達成任務重構的目的。

限制采用有限元模型的因素主要是其仿真計算速度較慢，且單元剛度矩陣、組裝和總剛度矩陣的計算均較為耗時，尤其當飛行器本身構建單元數(shù)量較大時，無法在一個或幾個制導周期內(nèi)完成結構強度的校核。目前，GPU計算在有限元分析中的應用尚處于初級階段[8]。許棟等[9]針對淺水問題利用GPU完成了多達112倍的加速比，并對核函數(shù)參數(shù)進行了優(yōu)化分析；曾良等[10]通過共享內(nèi)存邊界的迭代法，實現(xiàn)了流體計算的加速問題；Cecka等[11]針對熱傳導問題的有限元并行組裝進行了深入研究；Dinh等[12]則對GPU的有限元實時仿真問題做了一定研究，并總結了COO(Coordinate)格式下的單元組裝問題；Dziekonski等[13]討論了在多GPU環(huán)境下，大型系數(shù)線性方程的生成問題。

本文結合目前主流飛行器構型，利用GPU浮點運算能力極強的特點，著重從單元剛度矩陣的并行組裝方面進行詳細研究及設計，為實際的飛行仿真及任務提供技術驗證。

2 有限元矩陣并行裝配算法設計

在有限元法求解步驟中，單元剛度矩陣完成計算后下一步就是裝配總剛度矩陣，建立總平衡方程組以求解節(jié)點位移。當節(jié)點數(shù)目較多時，總剛度矩陣的階數(shù)也較大，由于計算機內(nèi)存的容量限制，有可能無法按照傳統(tǒng)方法進行存儲，在ANSYS中采用波前法[13]，即邊分析邊消去以避免總剛度矩陣的生成。不過鑒于波前法需頻繁交換CPU緩存和內(nèi)存的數(shù)據(jù)，限于內(nèi)存帶寬，在求解較大規(guī)模問題時往往較為費時[14]。同時，由于存在大量的內(nèi)存拷貝和索引量操作，裝配步驟在大型有限元模型的整體求解耗時中也占有相當比重。下文在內(nèi)存能夠完全容納總剛度矩陣的假設下，設計總剛度矩陣的并行算法。

2.1 有限元矩陣裝配并行性分析

總剛度矩陣根據(jù)單元類型的不同分別計算每個單元局部坐標下的剛度矩陣，經(jīng)局部坐標系下的單元剛度矩陣轉換到整體坐標系下后對號入座完成裝配。

[K]=∑[Ke]

(1)

由于鄰近單元共用某個節(jié)點，且可能在單元剛度矩陣中的位置一樣，則總剛度矩陣的裝配不能單純以單元向量形式并行化，否則可能出現(xiàn)相同地址的訪問沖突，進而造成計算結果錯誤。當且僅當每步并行裝配沒有共同節(jié)點的若干個單元的剛度矩陣，或至少保證每步中所有單元矩陣在總剛度矩陣中的位置不一時，并行裝配才得以實現(xiàn)。從實際出發(fā)，有如下方案。

(1)每次同時裝配一個單元剛度矩陣中多個位置的元素。

(2)每次同時裝配若干個單元剛度矩陣中同一位置的元素。

(3)每次同時裝配若干個單元剛度矩陣中若干位置的元素。

在單元剛度矩陣計算階段，每步均是以長度為r的向量形式對單元剛度矩陣中的某個子塊的某一處元素并行計算，由此可推斷每步裝配能同時計算r個單元的分塊單元剛度矩陣中塊[Ki,j]中(p，q)處的元素。而當r為單元總數(shù)時，并行裝配的每一步即為疊加所有單元的單元剛度矩陣中同一位置的元素。于是可提出，在何種情況下可同時裝配所有單元剛度矩陣同一位置處的元素滿足訪問沖突的要求，即如何保證總剛度矩陣中所有位置的元素在每一步裝配中均只受一個單元影響。結合GPU架構，適用于其SIMD架構的組裝算法應添加以下幾點關鍵思路。

(1)將任務分成若干個相互獨立的子塊，子塊間應不存在對同一地址的訪存沖突。

(2)盡量減少全局內(nèi)存的數(shù)據(jù)傳輸量，且盡量滿足內(nèi)存的合并訪問，提高分級內(nèi)存體系下的數(shù)據(jù)重用和訪問效率。

(3)平衡數(shù)值計算和數(shù)據(jù)流間的效率。

(4)盡量優(yōu)化每個線程塊內(nèi)線程數(shù)和流處理器內(nèi)的可用線程塊數(shù)目。

基于CUDA的并行組裝算法的難點在于，如何保證單元核函數(shù)高效地獲取節(jié)點數(shù)據(jù)并將單元數(shù)據(jù)組裝至系統(tǒng)方程。合并讀取和寫入全局內(nèi)存并且避免共享內(nèi)存訪存沖突在沒有顯著耗時的輸入數(shù)據(jù)重構的條件下幾乎不可避免。此外，通常多個單元共享某一個節(jié)點，則說明其共享自由度，故在整個并行計算中應該盡量避免全局同步機制和原子操作以期避免競態(tài)條件[15]。在此之前，文獻[16，17]介紹了FEM中稀疏矩陣的并行組裝問題,提出每個線程對應計算系統(tǒng)方程Ax=F的每一個非零值，以避免競態(tài)條件；然而，該方法需要每個線程計算單元子程序并從每個貢獻單元的數(shù)據(jù)中提取相關條目，因此，該做法會導致很多的冗余計算。如果ef是每個單元的自由度數(shù)，則單元子程序需要調(diào)用(ef+1)ef次(前一項表示對每個非對稱矩陣 中的非零項，后一項表示在Ae的非零項)。盡管單元特征矩陣的計算過程因單元特性而異可能計算開銷并不大，但多次從全局內(nèi)存中檢索單元子程序所需的節(jié)點數(shù)據(jù)，而非以合并訪問的方式進行，則時間開銷非常大。

考慮組裝算法的通用性，暫且認定單元矩陣的形式不一，僅認為單元子程序作為黑箱在每個線程或線程塊內(nèi)運行，實現(xiàn)組裝過程與單元特征矩陣計算的分割獨立。文獻[18]表明，該方案適宜于低階單元并取得了較好的加速比，同時基于該思路針對高階單元的優(yōu)化在文獻[19]中也已有闡述。則在確定單元特征矩陣與裝配過程解耦后，即確定裝配形式，即GPU內(nèi)的線程或線程塊與矩陣中元素間的關聯(lián)，如某一個某一行(列)、非零元素或非零項子集等。同時，考慮可用的GPU資源，由于其包含了三種分級內(nèi)存系統(tǒng)[20]，可得到如下基于GPU架構的典型并行方案。

表1 數(shù)據(jù)的存儲與硬件資源的關系

數(shù)據(jù)存儲-線程關聯(lián)的組裝策略列入表1，每種方案均有其優(yōu)劣，其中文獻[18]闡述了基于本地寄存器的NZ組裝算法、全局內(nèi)存的行組裝算法以及基于不同顏色區(qū)分的單元組裝算法。故下文主要闡述各并行裝配算法在GPU上的實現(xiàn)及優(yōu)化,其中第一部分詳細描述依照顏色劃分單元的裝配策略，保證兩個線程不會同時寫入某一相同內(nèi)存地址的具體實現(xiàn)；第二部分闡述第二種算法，即將計算分為兩部分，計算所有的單元數(shù)據(jù)并存入全局內(nèi)存，接著組裝所有的單元數(shù)據(jù)至矩陣A；第三種算法較前兩種更復雜，使用了共享內(nèi)存以減少全局內(nèi)存的通信量,節(jié)點和單元皆分割成若干個子集，每個子集調(diào)用若干線程計算單元特征矩陣并計算單元在系統(tǒng)方程內(nèi)相應的位置。方案三能夠減少節(jié)點數(shù)據(jù)的總讀取量，避免了單元數(shù)據(jù)的重復執(zhí)行并減少了全局內(nèi)存的寫入量。

2.2 并行裝配算法

2.2.1 以顏色分組的組裝策略

避免競態(tài)條件的直接方法是在組裝算法前對網(wǎng)格單元進行著色預處理，使得任意兩個相同顏色的單元不共享自由度，如圖1所示。盡管確定一個模型的最小著色數(shù)是一個NP(Non-deterministic Polynomial)完全問題，但已有多種啟發(fā)式算法[21,22]能夠找到近似的最小著色數(shù)。

本文選用網(wǎng)格的對偶圖的頂點著色，其中每個單元代表圖中的節(jié)點，兩個節(jié)點若通過一條邊相連則表示兩個單元共線或者共面。需注意應保證每種顏色的單元數(shù)目盡量相近，因為在CUDA中各流的數(shù)據(jù)量和計算量盡量相同才能發(fā)揮流的最大作用[23]；同時，若某種顏色的單元數(shù)量不足可能啟動較少的線程塊，進而導致GPU資源的低效使用。

圖1 2維單元著色

Fig.1 2-D element coloring

在此，不妨將著色為k的單元集合寫做εk。在集合εk中單元之間不共享節(jié)點，所以不用優(yōu)化從全局顯存中讀取的每個單元的節(jié)點數(shù)據(jù)。相反，每個元素從全局內(nèi)存中檢索自己的節(jié)點數(shù)據(jù)，運行單元核函數(shù)計算單元剛度和質量矩陣，并將該數(shù)據(jù)組裝至系統(tǒng)方程。則預處理每個著色的單元矩陣Ek為

Ek(σk(e),a)=E(e,a)

(2)

σk:εk→{1,…,|εk|}是全局單元號至著色單元號的映射。這些矩陣放置在內(nèi)存中，線程對應單元內(nèi)存的合并訪問。每個著色為k的組裝核函數(shù)流程如下。

(1)Tid←globalThreadID;

(2)For alln∈[1,en]do

(3)Nodes[n]←C(Ek(tid,n))

(4)End for

(5)計算(Ae,Fe)=elem(nodes)至系統(tǒng)方程

在上述算法中，步驟(5)取決于稀疏矩陣的存儲格式。為保證程序的通用性，在此提前計算本地單元數(shù)據(jù)至目標系統(tǒng)方程的索引數(shù)組中。為保證讀取的合并訪問，數(shù)據(jù)格式安排如圖2所示。一個列主元矩陣存儲每個單元的節(jié)點號(灰色部分)，其后表示每個單元數(shù)據(jù)與系統(tǒng)方程的映射關系(白色部分)。假定有一個一階線性三角單元，可知矩陣的列數(shù)為三個點的坐標數(shù)據(jù)加上6個單元特征矩陣Ae及3個等式右端的Fe項。

2.2.2 使用全局顯存的非零項組裝策略

該策略采用一個線程(或線程塊)計算一個有限單元數(shù)據(jù)，為避免競態(tài)條件，先將單元數(shù)據(jù)寫入全局顯存以便稍后的系統(tǒng)方程構建。受GPU硬件架構的限制，全局顯存內(nèi)的單元數(shù)據(jù)在構造系統(tǒng)方程時無法執(zhí)行全局的同步操作，故組裝系統(tǒng)方程和計算單元特征矩陣需執(zhí)行兩個不同的核函數(shù)。鑒于篇幅，針對全局顯存的非零項組裝策略僅介紹單元特征矩陣的核函數(shù)并行計算部分。

圖2 單元數(shù)據(jù)及其附加映射數(shù)據(jù)在GPU內(nèi)存存儲格式

Fig.2 Storage format of unit data and its additional mapping data in GPU memory

在一般情況下，因單元內(nèi)節(jié)點編號并不連續(xù)，而節(jié)點結構體數(shù)組通常是按編號排序的，故通常在讀取該部分數(shù)據(jù)時無法實現(xiàn)合并訪問。在此不妨設置一個節(jié)點數(shù)據(jù)矩陣Ce，其中，Ce(e,a)是單元e的本地節(jié)點a的節(jié)點數(shù)據(jù)。如果節(jié)點數(shù)據(jù)不發(fā)生變化，則可實現(xiàn)加速；若節(jié)點數(shù)據(jù)發(fā)生變化，則Ce需在每個時步進行更新。與上述算法對比，本策略中Ce刷新需要額外的全局顯存帶寬，通常會導致計算速度更慢。故可將單元組織為若干個子集，其中子集內(nèi)的單元可以共享節(jié)點數(shù)據(jù)，使得經(jīng)由全局顯存的數(shù)據(jù)傳輸總量通過預取數(shù)據(jù)實現(xiàn)線程塊內(nèi)的數(shù)據(jù)共享。在此，不妨假定εk是第k個線程塊的單元集合，計算節(jié)點集合Nk的公式為

Nk={E(e,a)∶e∈εk(a=1,…,en)}

(3)

線程塊k需檢索其對應的單元子集εk，然后該單元子集所需的節(jié)點數(shù)據(jù)經(jīng)全局顯存存入共享內(nèi)存。由于GPU內(nèi)的共享內(nèi)存有大小限制，為確保核函數(shù)不會調(diào)用失敗，Nk有如下約束。

(4)

SharedMemSize是以字節(jié)數(shù)表示的共享內(nèi)存大小，Dn是單節(jié)點的數(shù)據(jù)大小，無論單雙精度均轉換成字節(jié)即可獲得Nk的上限。對后續(xù)章節(jié)可能存在的非線性及多步插值問題，|Dn|=nd+nf，其中nd表示坐標數(shù)據(jù)nf牛頓-辛普森迭代的中間值或多步插值中出現(xiàn)的半步長中間量。

確定了滿足公式(4)的Nk后，即可將單元分為若干等分子集。選用圖論的劃分方法，保證每個單元是圖中的一個節(jié)點，若相應的單元間共面則在兩節(jié)點間添加一條連線。METIS[24]提供了較優(yōu)的劃分算法，保證得到一個可接受的Nk值。待Nk經(jīng)計算確定，即可將節(jié)點數(shù)據(jù)存儲至共享內(nèi)存并在單元子程序內(nèi)參與計算。對于每個線程塊k，預處理塊單元矩陣Ek定義為

Ek(e,a)=σk(Ek(e,a))(?e∈εk,a=1,…,en)

(5)

式中σk∶Nk→{1,…,|Nk|}表示第k個塊內(nèi)全局節(jié)點號E(e,a)至塊節(jié)點號Ek(e,a)的映射向量，存儲在共享內(nèi)存中并用于索引節(jié)點數(shù)據(jù)。在每個線程塊k內(nèi)，若干個en維元組連續(xù)存儲以便全局顯存讀取的合并訪問；待完成后便得到一個blockSise行的列主元矩陣，如圖3所示為一個en=3的示意圖。一個線程讀取索引表后即可將共享內(nèi)存內(nèi)的節(jié)點數(shù)據(jù)指向所需的單元子程序；隨后，單元子程序計算完單元特征矩陣，并將其以合并訪問的方式回傳至全局顯存。

2.2.3 使用共享內(nèi)存的非零項組裝策略

共享內(nèi)存較寄存器更大，其訪問延遲較全局內(nèi)存低得多，且塊內(nèi)線程間可以共享單元數(shù)據(jù)以減少子程序對數(shù)據(jù)調(diào)用的總時長，故提出塊內(nèi)線程處理某一節(jié)點子集的思路。假設N是所有節(jié)點的集合，有Nk?N，于是有單元集合εk為

εk={e∈ε?a,E(e,a)∈Nk}

(6)

(7)

(8)

式中|Dn|是輸入節(jié)點數(shù)據(jù)的尺寸，|De|是單元子程序生成的單元數(shù)據(jù)Ae和Fe的大小。對于非線性問題,若采用高階項模擬，假定一個單元內(nèi)的任意節(jié)點有nd個自由度和nf個中間量，可推導

|Dn|=en(nd+nf)

(9)

圖3en=3的列主元預處理塊單元矩陣EI

Fig.3 Column primary element matrix ofen=3

圖4 節(jié)點集合Nk與其對應單元集合εk

Fig.4 Node setNk-cell setεk

此外，|De|=ef(ef+1)，該參數(shù)是本地單元矩陣Ae和力向量Fe的尺寸。考慮單元特征矩陣的對稱性，則有

(10)

式中n*因單元維度和形狀而異。

同時，由于單元之間共享節(jié)點，最優(yōu)的也是最重要的優(yōu)化措施是所有的節(jié)點數(shù)據(jù)僅從全局內(nèi)存中讀取一次，否則會出現(xiàn)多個單元多次讀取某節(jié)點數(shù)據(jù)，可能出現(xiàn)訪存沖突，也使得計算的實時性無法保證。因此，與全局內(nèi)存的非零組裝算法不同，在此以節(jié)點數(shù)據(jù)為導向作為源索引，其后存儲共用該節(jié)點的單元，進而形成列主元的歸約矩陣，如圖5所示。其灰色部分表示全局節(jié)點號，白色代表共用該節(jié)點的單元列表。區(qū)分兩種不同索引的方式參考了Cris等[11]的策略,為區(qū)分這兩部分整數(shù)索引，通過反碼再減1的方式表示目標索引，通過加1表示目標源索引；程序內(nèi)每個線程先判斷整數(shù)的正負號，依據(jù)正負號區(qū)分兩種索引類型。隨著線程依次讀取列表內(nèi)數(shù)據(jù)，將所有指向源向量S的數(shù)據(jù)累加并存儲至目標向量T的目標索引，如下，

(1)Fori:NZ_Reduction

(2)Ifi>0

(3)t←t+S[i-1]

(4)else ifi<0

(5)T[-i-1]←t

(6)End if

(7)End for

其組裝方式如圖6所示。對于常用的低階單元,其單元矩陣尺寸較小，輔以節(jié)點為索引作為輸入量，使得劃分的單元集合尺寸較大，理論效率較高。可以看出，在預處理階段如果Nk內(nèi)部分或所有節(jié)點號連續(xù)，有些全局內(nèi)存的讀取也是可以合并訪問的。同時,為避免共享內(nèi)存的存儲體讀取沖突，分配給每個單元的共享內(nèi)存地址不能連續(xù)。

圖5 線程塊內(nèi)blockSize行的列主元形式存儲的歸約矩陣

Fig.5 Reduction matrix stored in column primitive form of block size row

3 并行組裝算法算例對比

3.1 測試平臺與對象

本文測試平臺硬件選用單機多GPU工作站，采用上述算法在自行設計的并行求解器框架下，模擬細長體彈性飛行器飛行過程，并對單GPU及多GPU并行的性能及拓展性進行對比，列入表2。因為CUDA 9.0版本不支持Visual studio新版本的一部分宏，為避免自行配置發(fā)生未知錯誤，在保證支持編譯C++11標準的編譯器中選擇MSVC14.0。為進一步驗證單機多GPU節(jié)點的耗時情況，帶*號的GPU1作為添加的另一部GPU，以期實現(xiàn)單計算節(jié)點的并行任務分割。

同時，鑒于Windows并非真正的實時操作系統(tǒng)，其CPU端計時器精度僅有1 ms，故在進行時間統(tǒng)計時采用CUDA自帶的cudaEvenet_t計時，其精度為0.5 μs。若干組適宜于航天器典型構型或特征的測試算例列入表3。

3.2 計算精度及效率分析

首先，測試METIS對于單元分類的計算耗時情況。在使用Release且不進行優(yōu)化，網(wǎng)格從文件中讀取的情況下,測試顏色分組耗時，各算例情況如圖7所示。可以看出，因為文件讀取內(nèi)容沒有變化，I/O耗時基本沒有太大變化，幾乎是隨著單元數(shù)目的增加而線性增加。但是兩種劃分方式隨單元數(shù)目的變化有著非常大的變化,紅黑樹在前幾個單元數(shù)量較少的算例中表現(xiàn)優(yōu)異，但從編號3開始,K路劃分的耗時增幅相對較少；當單元數(shù)達到80000左右時，兩者在劃分成22種顏色的情況下耗時比已達到2.5倍左右。從劃分顏色數(shù)量的變化趨勢看，劃分的顏色越多，邊切割數(shù)越多，必然導致耗時越多；但是若在只有一部圖形處理器的工作站上，沒有必要劃分很多種顏色進而產(chǎn)生非常多的流，那么在4～10種顏色劃分中其耗時變化不大，尤其在單元數(shù)較多時波動不超過5%，也驗證了Cris[11]的觀點；而在若干臺、每臺均有多于一部圖形處理器的工作環(huán)境下，在編號6中單元數(shù)較多時，其額外耗時不超過6 ms，相較于集群計算的通信延遲是能接受的。可見，在此使用METIS進行顏色分組是非常高效和合適的方法；同時，也能通過曲線擬合的方法，在已知單元數(shù)目條件下，選擇最優(yōu)的劃分算法。

圖6 利用共享內(nèi)存的非零值組裝策略流程

Fig.6 Flow chart of non-zero assembly strategy using shared memory

表2 計算平臺配置

Tab.2 Computing platform configuration

硬件/軟件配置CPUIntelE31230V5內(nèi)存DDR4240016GGPU0GTX10606GGPU1?GTX10603G操作系統(tǒng)WindowsCUDA版本9.0編程環(huán)境MSVC14.0

表3 飛行器典型構型或部件算例

Tab.3 Typical configuration and components list of aircraft

編號單元名稱單元數(shù)量節(jié)點數(shù)算例描述12D4N5066長方形平面薄板2C3D8120220空間梁型結構3S4R-S337463658升力體飛行器機翼42D3N74344038三角面元機翼5C3D82006240736空間組合梁型結構6S4R-S38877588661升力體飛行器完整構型

圖7 各測試用例劃分耗時

Fig.7 Time-consuming for partitioning test cases

為保證裝配算法在常用的二次單元中同樣有較好的通適性，選取表3的編號3模型和編號5模型，分別將各單元替換為二階四邊形-三角形殼元和二階六面體等參單元。并行裝配算法橫向對比結果如圖8所示。以雙精度計算并調(diào)用cudaThreadSynchronize函數(shù)保證線程同步后截止計時，記錄不包括主機-設備端傳輸耗時之外的所有耗時項。由于高階單元需要更多的計算量、共享內(nèi)存及寄存器空間，其耗時相應也較低階單元更多；同時如圖8(b)所示,在單元數(shù)較多時，基于共享內(nèi)存的非零組裝算法在二階單元情況下獲得了約3的加速比，單純從組裝算法耗時看能夠滿足至多 25 ms 步長的實時仿真要求，而顏色重排序算法僅獲得了約2倍的加速效果，無法滿足實時性要求。但是值得注意的是，圖8兩算例的耗時均在單臺工作機上統(tǒng)計，而顏色分割重排序算法可將任務分配至計算集群內(nèi)的多個計算節(jié)點同時組裝，由于GPU組裝相較于基于顏色劃分的組裝算法會產(chǎn)生額外的CPU-GPU內(nèi)存拷貝操作，因此顏色重排序后將任務分配至各計算節(jié)點的策略在此情況下的計算耗時可能比共享內(nèi)存的非零組裝算法更少。

為進一步驗證并行組裝算法在單機多GPU節(jié)點的耗時情況，經(jīng)修改主程序，以單元顏色分組后的單雙號為識別，將單色號單元分配至GPU0內(nèi)組裝，雙色號單元分配至GPU1內(nèi)組裝，獲得如圖9所示的結果。對比在計算平臺1的計算耗時，顏色分割重排序算法沒有發(fā)生明顯的時間變化，其原因是在該算法中，METIS進行顏色分割的耗時占據(jù)了很大一部分，其作為在CPU上運行的串行算法并不能因新GPU的引入而實現(xiàn)性能提升，符合Amdahl定律。但基于共享內(nèi)存和全局顯存的非零組裝算法均有一定的性能提升，尤其編號5的二階單元組裝算例中，因單元數(shù)較多能夠啟動足夠多的核函數(shù)，基于共享內(nèi)存的非零組裝算法幾乎呈線性加速，耗時控制在10 ms以內(nèi)，實現(xiàn)了約6的加速比。

圖8 一階、二階單元不同組裝算法耗時對比

Fig.8 Time comparison of different assembly algorithms for first-order and second-order elements

圖9 雙GPU下不同階數(shù)單元組裝耗時

Fig.9 Time-consuming comparison of different order unit assembly under dual GPU

4 結 論

為了設計適用于彈性飛行器飛行動力學仿真及結構快速故障診斷平臺，作為初步的算法研究，設計了不同種類的有限元組裝算法，以常用的著色算法為基礎，提出了全局顯存和共享內(nèi)存的非零元素組裝算法。以航天器典型部件為算例，提出了適宜于不同階數(shù)單元的組裝算法，在雙精度下采用單機雙GPU達到了約6倍的加速比，用以往25 ms的制導周期作為判據(jù)，二階單元在20000余個時仍能通過共享內(nèi)存的非零組裝模式滿足組裝的實時性要求。但沒有討論單元矩陣存儲模式對組裝算法的影響，也未討論單精度算法能夠在滿足精度的條件下允許更多的單元滿足實時組裝的需求。以上兩點將是未來的研究重點。