李果陽(yáng)，嚴(yán)華，*，張征宇，陳沁梅，祝福順

(1.四川大學(xué) 電子信息學(xué)院，成都 610065；2.中國(guó)空氣動(dòng)力研究與發(fā)展中心 高速空氣動(dòng)力研究所，綿陽(yáng) 621000）

對(duì)于風(fēng)洞試驗(yàn)，空氣流動(dòng)包含著最重要的特征信息，對(duì)處于空氣流動(dòng)形成的流場(chǎng)中物理模型的密度場(chǎng)進(jìn)行定量分析一直都是流場(chǎng)的重要研究方向［1］。傳統(tǒng)的紋影技術(shù)由于在定量測(cè)量上效果不好而常常作為一種定性分析的方法，干涉技術(shù)由于需要嚴(yán)苛的要求才能實(shí)現(xiàn)，缺乏實(shí)用性。因此，背景紋影技術(shù)(Background Oriented Schlieren，BOS）就在這種情況下被Meier［2］提出。運(yùn)用該方法與風(fēng)洞中的視頻測(cè)量技術(shù)相結(jié)合可以較簡(jiǎn)單地實(shí)現(xiàn)對(duì)處于風(fēng)洞試驗(yàn)中的物理模型進(jìn)行密度場(chǎng)的視頻測(cè)量［3-6］，通過(guò)測(cè)出一束非平行光在流場(chǎng)擾動(dòng)下的偏折角，進(jìn)而建立偏折角與擾動(dòng)流場(chǎng)折射率間的量化關(guān)系，得到關(guān)于密度場(chǎng)的二階偏微分方程。利用視頻測(cè)量法得到關(guān)于坐標(biāo)及坐標(biāo)偏移量的離散數(shù)據(jù)，再對(duì)于這樣的密度場(chǎng)中二階偏微分方程求解出非平行光經(jīng)過(guò)擾動(dòng)流場(chǎng)后對(duì)應(yīng)位置的密度投影值。目前有多種常用數(shù)值模擬計(jì)算方法，經(jīng)過(guò)比較，有限元法因擅長(zhǎng)處理各種復(fù)雜區(qū)域、精度較高及適于并行計(jì)算等優(yōu)點(diǎn)［7-8］，因此本文將采用有限元法求解密度場(chǎng)對(duì)應(yīng)位置的投影值。

從有限元的計(jì)算方法可以看出，在用有限元法進(jìn)行大規(guī)模數(shù)值計(jì)算時(shí)，不僅要考慮多種力學(xué)行為的全耦合求解問題，計(jì)算過(guò)程復(fù)雜，而且總剛度方程的求解時(shí)間占據(jù)了整個(gè)計(jì)算過(guò)程很大的比重［9-10］。由于神經(jīng)網(wǎng)絡(luò)具有強(qiáng)大的映射、實(shí)時(shí)和并行計(jì)算能力［11］，為了解決有限元中的總剛度矩陣求解耗時(shí)，簡(jiǎn)化有限元法求解的復(fù)雜度問題，本文不再采用傳統(tǒng)有限元計(jì)算方式，而是引入了神經(jīng)網(wǎng)絡(luò)進(jìn)行求解。另外，風(fēng)洞試驗(yàn)中需要實(shí)時(shí)分析風(fēng)洞試驗(yàn)密度場(chǎng)的分布變化，由此需要快速計(jì)算密度場(chǎng)。但是，神經(jīng)網(wǎng)絡(luò)訓(xùn)練不僅非常耗時(shí)且會(huì)占據(jù)計(jì)算機(jī)大部分資源，而且有限元計(jì)算過(guò)程中產(chǎn)生的大規(guī)模稀疏線性方程組數(shù)據(jù)量也很巨大［12］，計(jì)算非常耗時(shí)，嚴(yán)重影響了密度場(chǎng)的求解效率。對(duì)于風(fēng)洞試驗(yàn)密度場(chǎng)這樣的時(shí)變系統(tǒng)要求實(shí)時(shí)分析，其核心就在于能夠?qū)崟r(shí)計(jì)算，因此，非常有必要探尋求解密度場(chǎng)的快速計(jì)算方法。目前主流的高效計(jì)算方法主要采用基于GPU的并行計(jì)算方式［13-15］，因此本文主要就針對(duì)求解過(guò)程中耗時(shí)的神經(jīng)網(wǎng)絡(luò)、剛度矩陣和載荷向量的求解進(jìn)行并行加速，以實(shí)現(xiàn)密度場(chǎng)的快速求解。

為此，本文提出了一種高效求解密度場(chǎng)的思路，給出了神經(jīng)網(wǎng)絡(luò)與傳統(tǒng)有限元分析融合的有效方法，以及利用GPU快速求解密度場(chǎng)的計(jì)算方法和具體步驟，實(shí)驗(yàn)結(jié)果證明了所提方法的有效性。

1 有限元法求解密度場(chǎng)

1.1 密度場(chǎng)的二階偏微分方程

利用背景紋影技術(shù)的方法測(cè)量風(fēng)洞試驗(yàn)中的流場(chǎng)［16-17］。坐標(biāo)系的x軸為逆氣流方向，y軸與x軸垂直，指向試驗(yàn)段上壁板，z軸與x軸垂直，與光線傳輸方向相反。為了便于分析，不妨假設(shè)流動(dòng)介質(zhì)的折射率僅與y軸方向的坐標(biāo)值有關(guān)，且隨y軸方向坐標(biāo)值的增加而增大，沿z軸方向平行射入的光線，在y-z截面內(nèi)發(fā)生偏折，光線偏折角分析如圖1所示。

圖1 光線的偏折角分析示意圖Fig.1 Schematic diagram of ray deflection angle analysis

折射率n與氣體密度ρ之間的Gladstone-Dale公式為

由光線傳輸方程和式(1），可得y軸方向偏折角εy與氣體密度ρ之間的關(guān)系為

同理可得，x軸方向偏折角為

式中：L為風(fēng)洞試驗(yàn)段的邊長(zhǎng)；KGD為Gladstone-Dale常數(shù)，在空氣中標(biāo)準(zhǔn)狀態(tài)下取值為0.226 cm3／g。

對(duì)式(4）和式(5）中y、x分別微分再相加即可得到密度場(chǎng)應(yīng)滿足的二階偏微分方程：

式中：ρm為密度場(chǎng)密度。

將式(6）的等式左端記為Δu，等式右端記為

其求解域?yàn)榭涨荒Ｐ?，記為Ω?/p>

考慮式(6）的變分問題，對(duì)于任意的v∈HΩ，同時(shí)乘上式(7）兩邊，使用格林公式，利用邊界條件可得

可證明得到，在一定連續(xù)性要求下，式(6）可等價(jià)于求u∈V，使得a(u，v）＝g(u）對(duì)任意的v∈V成立。

考慮上述變分問題的有限維逼近，即在有限維子空間Vh?V下，求uh∈Vh，使得a(uh，vh）＝g(uh）對(duì)任意的vh∈Vh成立。設(shè)在Vh下的一組基函數(shù)為，并依次取vh為每個(gè)基函數(shù)，則可得

分析式(11），其本質(zhì)是一個(gè)線性方程組：

基于上述分析，在初邊值的條件約束下，由Lax-M ilgram定理可知，式(6）即可等價(jià)于有限元法下求解一個(gè)線性方程組［18-19］：

式中：K為單元力和單元位移關(guān)系間的系數(shù)矩陣，代表了單元的剛度特性，稱為單元?jiǎng)偠染仃?；F為結(jié)點(diǎn)產(chǎn)生單位位移時(shí)，在結(jié)點(diǎn)上所需要施加的結(jié)點(diǎn)力，稱為結(jié)點(diǎn)載荷；δ為結(jié)點(diǎn)受到擾動(dòng)流場(chǎng)作用后，結(jié)點(diǎn)在其自由度方向偏移的距離，稱為結(jié)點(diǎn)位移，且K、F、δ的具體表達(dá)式如下：

式中：D為彈性矩陣；t為單元厚度；B為單元應(yīng)變矩陣；J為雅可比矩陣；α1、α4為剛體位移；α2、α3、α5、α6為單元的常應(yīng)變；x、y為結(jié)點(diǎn)坐標(biāo)。

1.2 密度場(chǎng)的有限元法求解過(guò)程

傳統(tǒng)的有限元分析是利用變分原理將求解域剖分求解的一種數(shù)值計(jì)算方法，具體步驟如下：

1）連續(xù)體離散化。用有限個(gè)離散單元體集合代替原有的連續(xù)體，也稱網(wǎng)格剖分，即進(jìn)行單元?jiǎng)澐?，將全部單元和結(jié)點(diǎn)按一定順序編號(hào)，每個(gè)單元所受荷載均按靜力等效原理作用在結(jié)點(diǎn)上，并根據(jù)實(shí)際情況在位移受約束的結(jié)點(diǎn)上設(shè)置約束條件。

2）單元分析。建立各個(gè)單元的結(jié)點(diǎn)位移δe和結(jié)點(diǎn)力Fe之間的關(guān)系式Fe＝Keδe，Ke表示單元?jiǎng)偠染仃嚕磳卧慕Y(jié)點(diǎn)位移作為基本變量，確定一個(gè)近似表達(dá)式，再應(yīng)用流體力學(xué)理論和虛功原理得到單元中結(jié)點(diǎn)力與結(jié)點(diǎn)位移的關(guān)系式。

3）整體分析。對(duì)各個(gè)單元組成的整體進(jìn)行分析，其目的是根據(jù)一定法則建立一個(gè)線性方程組，得到全部結(jié)點(diǎn)載荷F與全部結(jié)點(diǎn)位移δ的關(guān)系式Kδ＝F，從而求解出結(jié)點(diǎn)位移。

本文所實(shí)現(xiàn)的快速有限元法采用了神經(jīng)網(wǎng)絡(luò)代替單元分析中求解結(jié)點(diǎn)力與結(jié)點(diǎn)位移關(guān)系的復(fù)雜計(jì)算，從而使得整體分析計(jì)算過(guò)程簡(jiǎn)化，可大大縮短整體計(jì)算時(shí)間。

1.2.1 網(wǎng)格剖分

利用有限元法求解泊松方程首要的就是將求解域進(jìn)行離散，劃分得到有限個(gè)互不重疊的單元網(wǎng)格，同時(shí)，對(duì)這些單元和結(jié)點(diǎn)進(jìn)行編號(hào)，這些互不重疊的網(wǎng)格單元通過(guò)結(jié)點(diǎn)連接相互影響。在劃分單元網(wǎng)格時(shí)，對(duì)計(jì)算結(jié)果有著重要影響的是網(wǎng)格單元尺寸大小、網(wǎng)格的疏密程度及網(wǎng)格的類型。三角單元網(wǎng)格能夠適應(yīng)復(fù)雜的求解域，單元大小劃分十分方便，但計(jì)算精度相對(duì)較低，而四邊形單元?jiǎng)澐植灰酌枋龇钦坏闹边吔绾颓吔?，但其?jì)算精度相對(duì)較高。

本文采用的是四邊形單元和三角形單元結(jié)合的方式，實(shí)現(xiàn)了一種四邊形網(wǎng)格單元與三角形網(wǎng)格單元相結(jié)合的方式，剖分示意圖如圖2所示。在求解時(shí)為便于網(wǎng)格的剖分實(shí)現(xiàn)，采用了先對(duì)求解域進(jìn)行等參四邊形單元剖分，再將一個(gè)四邊形單元一分為二，形成兩個(gè)等參的三角形單元剖分算法，然后對(duì)它們進(jìn)行單元編號(hào)和結(jié)點(diǎn)編號(hào)。

圖2 網(wǎng)格剖分平面示意圖Fig.2 Grid subdivision p lan sketch

1.2.2 神經(jīng)網(wǎng)絡(luò)擬合

利用神經(jīng)網(wǎng)絡(luò)代替?zhèn)鹘y(tǒng)有限元法中的單元分析，關(guān)鍵在于要清楚單元分析得到的單元?jiǎng)偠染仃嚨奶匦浴?/p>

對(duì)于本文所采用的三角單元剖分，根據(jù)彈性力學(xué)二維有限元法［18］，可知有限元法中單元?jiǎng)偠染仃囉幸韵聨c(diǎn)性質(zhì)：

1）單元?jiǎng)偠染仃囀菍?duì)稱方陣。單元?jiǎng)偠染仃嚸啃性刂蜑?，由對(duì)稱性質(zhì)可知，每列元素之和也為0。

2）對(duì)于平面單元，單元?jiǎng)偠染仃嚥粫?huì)隨著幾何參數(shù)改變。

3）根據(jù)單元?jiǎng)偠染仃嚨奈锢硪饬x可以得出其只與單元結(jié)點(diǎn)、相對(duì)坐標(biāo)及形函數(shù)等相關(guān)，與結(jié)點(diǎn)絕對(duì)坐標(biāo)無(wú)關(guān)。

由以上單元?jiǎng)偠染仃嚨男再|(zhì)結(jié)論可推出，剛度矩陣的計(jì)算實(shí)質(zhì)是從單元的尺寸、材料性質(zhì)到剛度矩陣的映射問題，而神經(jīng)網(wǎng)絡(luò)的優(yōu)勢(shì)在于：可以方便地建立輸入和輸出間復(fù)雜的映射關(guān)系。因此，利用神經(jīng)網(wǎng)絡(luò)可以很方便地從單元結(jié)點(diǎn)坐標(biāo)映射到單元?jiǎng)偠染仃囋亍?/p>

本文中，由于密度場(chǎng)二階偏微分方程右端項(xiàng)是關(guān)于光線偏折的一階微分式，且實(shí)驗(yàn)數(shù)據(jù)是使用視頻測(cè)量方法在沒有流場(chǎng)擾動(dòng)時(shí)，拍攝一幅包含高密度標(biāo)記點(diǎn)的背景板圖像作為參考圖像，然后加入擾動(dòng)流場(chǎng)，再拍攝背景板的時(shí)序圖像，根據(jù)視頻測(cè)量的標(biāo)定原理測(cè)得背景板標(biāo)記點(diǎn)在坐標(biāo)系下的(x，y）坐標(biāo)值及其偏移量Δx和Δy。因此在求解時(shí)需要對(duì)數(shù)據(jù)進(jìn)行預(yù)處理，即根據(jù)視頻測(cè)量法原理和光線偏折理論得到坐標(biāo)偏移量和偏折角的關(guān)系，再對(duì)其做歸一化處理，接著采用神經(jīng)網(wǎng)絡(luò)對(duì)偏移角擬合得到二階偏微分方程右端項(xiàng)一階微分式，即通過(guò)對(duì)右端項(xiàng)的一階微分式在點(diǎn)(x，y）處根據(jù)泰勒級(jí)數(shù)展開式可以得到

式中：h為坐標(biāo)點(diǎn)更新的步長(zhǎng)。

神經(jīng)網(wǎng)絡(luò)擬合模型如圖3所示。將經(jīng)過(guò)預(yù)處理的標(biāo)記點(diǎn)坐標(biāo)數(shù)據(jù)(其中80%為訓(xùn)練樣本，20%為驗(yàn)證樣本）作為網(wǎng)絡(luò)輸入，通過(guò)一個(gè)含有400個(gè)神經(jīng)元的隱含層，激活函數(shù)采用tanh型，優(yōu)化函數(shù)采用目前普遍使用的Adam進(jìn)行有監(jiān)督的學(xué)習(xí)。將網(wǎng)格單元結(jié)點(diǎn)通過(guò)該訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)，最終實(shí)現(xiàn)網(wǎng)格單元的結(jié)點(diǎn)通過(guò)神經(jīng)網(wǎng)絡(luò)映射到單元?jiǎng)偠染仃囋亍?/p>

圖3 神經(jīng)網(wǎng)絡(luò)擬合模型Fig.3 Neural network fitting model

1.2.3 整體分析

總剛度矩陣K由單元?jiǎng)偠染仃嘖e組裝而成，計(jì)算如下：

式中：Ce為單元選擇矩陣。根據(jù)單元的結(jié)點(diǎn)編號(hào)和結(jié)構(gòu)中結(jié)點(diǎn)編號(hào)的對(duì)應(yīng)關(guān)系，將單元?jiǎng)偠染仃囘M(jìn)行式(19）迭加組裝成總體剛度矩陣，得到結(jié)點(diǎn)載荷與結(jié)點(diǎn)位移間的關(guān)系式，如式(13）所示。

然后進(jìn)行狄利克雷邊界條件約束，采用直接法求解得出結(jié)點(diǎn)位移。

2 基于GPU的有限元法并行求解

本文中將有限元法求解密度場(chǎng)的方法采用并行的思想，對(duì)求解過(guò)程極其耗時(shí)且適合并行計(jì)算的神經(jīng)網(wǎng)絡(luò)、總剛度矩陣和總載荷向量的求解進(jìn)行并行加速。融合了神經(jīng)網(wǎng)絡(luò)的單元分析過(guò)程，以時(shí)間代價(jià)簡(jiǎn)化了求解復(fù)雜度，為進(jìn)一步加快求解速度采取了數(shù)據(jù)并行、多線程同時(shí)計(jì)算的方式。整體分析中，總體剛度矩陣和總載荷向量的求解進(jìn)行分塊多線程并行計(jì)算。

2.1 CUDA架構(gòu)介紹

GPU具有高性能的浮點(diǎn)計(jì)算能力，CPU具有低延遲、邏輯計(jì)算能力強(qiáng)的特點(diǎn)，CUDA編程正是利用GPU和CPU這種特點(diǎn)來(lái)協(xié)同處理任務(wù)。運(yùn)行流程則是在CPU端準(zhǔn)備數(shù)據(jù)，再傳到GPU端并行計(jì)算，最后將計(jì)算結(jié)果回傳到CPU端。CPU和GPU的這種協(xié)同工作方式大大提高了計(jì)算性能［20-21］。圖4為CUDA線程組織結(jié)構(gòu)模型。

為了明白CUDA編程如何進(jìn)行線程組織，必須了解其結(jié)構(gòu)模型。從圖4中可以看到，CUDA中線程層次分為3層，分別是線程網(wǎng)格(grid）、線程塊(block）、線程(thread）。1個(gè)線程索引(thread Idx）可以是1～3維的，進(jìn)而多個(gè)thread組成的block也可以是1～3維的。如同線程組成線程塊，線程塊也能組成1～3維的線程網(wǎng)格。執(zhí)行CUDA應(yīng)用程序時(shí)，整個(gè)應(yīng)用程序的關(guān)鍵是kernel函數(shù)，它是以grid的形式組織，以block為單位，各block間并行執(zhí)行，不同block間的數(shù)據(jù)共享只能通過(guò)全局顯存，對(duì)于同一block內(nèi)的不同thread間則可以通過(guò)共享內(nèi)存進(jìn)行通信，即在kernel函數(shù)中存在2個(gè)層次的并行，即block之間并行計(jì)算及thread之間并行計(jì)算。

圖4 CUDA線程組織結(jié)構(gòu)模型Fig.4 CUDA thread organization structure model

另外，為了達(dá)到更高的效率，在CUDA 編程應(yīng)格外關(guān)注內(nèi)存的使用。在CUDA編程時(shí)要盡量使用寄存器，盡量將數(shù)據(jù)聲明為局部變量。而當(dāng)存在著數(shù)據(jù)的重復(fù)利用時(shí)，可以把數(shù)據(jù)存放在共享內(nèi)存里。而對(duì)于全局內(nèi)存，需要注意用一種合理的方式來(lái)進(jìn)行數(shù)據(jù)的合并訪問，以盡量減少設(shè)備對(duì)內(nèi)存子系統(tǒng)再次發(fā)出訪問操作的次數(shù)。

2.2 神經(jīng)網(wǎng)絡(luò)并行加速方法

本文中將有限元法求解密度場(chǎng)的方法采用并行的思想對(duì)求解過(guò)程極其耗時(shí)的部分進(jìn)行加速處理，其模型如圖5所示。數(shù)據(jù)擬合的并行化主要是針對(duì)神經(jīng)網(wǎng)絡(luò)的學(xué)習(xí)訓(xùn)練過(guò)程。首先，將在密度場(chǎng)測(cè)量得到的離散數(shù)據(jù)分為訓(xùn)練樣本、驗(yàn)證樣本和測(cè)試樣本，并將其暫放在內(nèi)存中，然后對(duì)神經(jīng)網(wǎng)絡(luò)進(jìn)行初始化參數(shù)，并且將訓(xùn)練樣本分成多個(gè)批量(batch），加載到GPU顯存中。在神經(jīng)網(wǎng)絡(luò)正向傳播過(guò)程中，將原本一整個(gè)數(shù)據(jù)分為多片，調(diào)用CUDA的kernel函數(shù)為每個(gè)batch數(shù)據(jù)分配處理線程。由于在CUDA中一個(gè)W rap線程束有32個(gè)線程，在進(jìn)行線程分配時(shí)應(yīng)當(dāng)盡量選擇32的整數(shù)倍，但是CUDA中一個(gè)線程塊里線程數(shù)量又不多于1 024，因而可以盡可能提高線程的利用率。從主機(jī)端調(diào)用kernel函數(shù)在GPU設(shè)備端對(duì)分片好的訓(xùn)練樣本在大量線程中被同時(shí)并行處理，這就極大地加快了正向傳播過(guò)程。同理，它的反向傳播過(guò)程也是各線程中進(jìn)行一系列的矩陣運(yùn)算和向量運(yùn)算。當(dāng)每個(gè)線程所傳遞的誤差信號(hào)都到達(dá)網(wǎng)絡(luò)的隱含層，就能得到權(quán)值偏導(dǎo)和偏置偏導(dǎo)，而每個(gè)batch上得到的是誤差信號(hào)均值，并且只對(duì)其本身參數(shù)偏導(dǎo)進(jìn)行更新，所有線程間互不干擾，所有的參數(shù)間也不相關(guān)，這就完成了一次并行迭代過(guò)程。

圖5 神經(jīng)網(wǎng)絡(luò)并行模型Fig.5 Neural network parallel model

2.3 總剛度矩陣組裝的并行化

2.3.1 總剛度矩陣的并行性分析

對(duì)于大規(guī)模的有限元分析中，串行的有限元過(guò)程中均不是直接在內(nèi)存中一次建立總的平衡方程組，因?yàn)閷?duì)于大規(guī)模的有限元問題，其結(jié)點(diǎn)總數(shù)也很大，導(dǎo)致總剛度矩陣的階數(shù)太大，從而計(jì)算機(jī)計(jì)算過(guò)程中內(nèi)存不夠。因此，串行有限元過(guò)程中一般有2種處理方法：一是將總剛度矩陣分塊存儲(chǔ)在外存中，二是邊單元分析邊消去，避免總剛度矩陣的生成。但是這2種方法都需要內(nèi)外存之間頻繁地進(jìn)行數(shù)據(jù)交換，因此對(duì)于結(jié)點(diǎn)數(shù)較多的大規(guī)模問題計(jì)算一次是非常耗時(shí)的。

總剛度矩陣是將所有單元?jiǎng)偠染仃嚢凑掌浣Y(jié)點(diǎn)編號(hào)對(duì)號(hào)入座組裝而成，所有單元的單元?jiǎng)偠染仃嚥荒芟蚩倓偠染仃囈淮涡缘佣桑駝t就可能在計(jì)算過(guò)程中發(fā)生計(jì)算錯(cuò)誤和存取沖突。進(jìn)一步分析可知，總剛度矩陣組裝最大的并行只可能發(fā)生在每次同時(shí)迭加互相沒有共同結(jié)點(diǎn)的多個(gè)單元的剛度陣所有元素，即實(shí)際采用的方案可以是每次同時(shí)迭加一個(gè)單元?jiǎng)偠染仃嚩鄠€(gè)位置處的元素，或者每次同時(shí)迭加多個(gè)單元?jiǎng)偠染仃囋谕晃恢锰幍脑?。另外，總剛度矩陣的組裝與其存儲(chǔ)格式密切相關(guān)，大型有限元問題產(chǎn)生的總剛度矩陣一般是稀疏對(duì)稱的。

2.3.2 總剛度矩陣的并行組裝

通過(guò)上文分析，本文中的問題則是一個(gè)稀疏對(duì)稱的總剛度矩陣，對(duì)于這樣的總剛度矩陣組裝過(guò)程很耗時(shí)但卻非常適合并行計(jì)算?？倓偠染仃囉擅總€(gè)網(wǎng)格單元的單元?jiǎng)偠染仃嚱M裝而來(lái)，每個(gè)網(wǎng)格單元的單元?jiǎng)偠染仃嚨闹抵慌c其形狀大小、對(duì)應(yīng)點(diǎn)次序及在整體坐標(biāo)系的方位有關(guān)，因此本次求解劃分的矩形網(wǎng)格單元內(nèi)的三角單元具有相同的單元?jiǎng)偠染仃?。圖6表示了劃分的任意一個(gè)網(wǎng)格位置結(jié)點(diǎn)標(biāo)號(hào)。

圖6 剖分的任一網(wǎng)格位置結(jié)點(diǎn)編號(hào)Fig.6 Node number of any grid position of subdivision

本文單元?jiǎng)偠染仃嚱M裝成總剛度矩陣的具體方法是：將單元?jiǎng)偠染仃嘖e的每個(gè)子塊Kij送到總剛度矩陣的對(duì)應(yīng)位置上去，然后進(jìn)行迭加即可得到總剛度矩陣K的子塊，對(duì)每個(gè)單元都進(jìn)行如此操作計(jì)算即可得到總剛度矩陣。

在執(zhí)行并行組裝總剛度矩陣時(shí)，關(guān)鍵是如何找出Ke中的子塊在K中的對(duì)應(yīng)位置，然后為每個(gè)Ke的組裝分配線程索引，這樣在并行化過(guò)程中給每個(gè)單元分配的線程求解得到的單元?jiǎng)偠染仃嚥拍苷_地組裝成總剛度矩陣，解決這個(gè)問題需要清楚單元中結(jié)點(diǎn)編碼(局部碼）和整體結(jié)構(gòu)中的結(jié)點(diǎn)編碼(總碼）之間的對(duì)應(yīng)關(guān)系。圖7中每個(gè)單元的3個(gè)結(jié)點(diǎn)按逆時(shí)針方向編為i、j、m，單元?jiǎng)偠染仃囍械淖訅K是按結(jié)點(diǎn)的局部碼排列的，而總剛度矩陣中的子塊是按總碼排列的，因此，在單元?jiǎng)偠染仃囍邪呀Y(jié)點(diǎn)的局部碼換成總碼，并把其中的子塊按照總碼次序重新排列，以單元U1為例，局部碼i1、j1、m1對(duì)應(yīng)于總碼m、m＋1、m＋x＋1，即單元?jiǎng)偠染仃嚭涂倓偠染仃囄恢脤?duì)應(yīng)關(guān)系如圖7表示。其中上三角網(wǎng)格單元和下三角網(wǎng)格單元的單元?jiǎng)偠染仃囅嗤?/p>

圖7 單元?jiǎng)偠染仃嚺c總剛度矩陣位置對(duì)應(yīng)關(guān)系Fig.7 Location of element stiffness matrix corresponding to total stiffness matrix

另外，有限元法產(chǎn)生的總剛度矩陣的階數(shù)一般很高，在解算時(shí)矩陣的階數(shù)很高和存儲(chǔ)容量有限是矛盾的，因此有必要知道總剛度矩陣的特性來(lái)解決節(jié)約存儲(chǔ)容量的問題。經(jīng)分析，總剛度矩陣的特性如下：

1）對(duì)稱性。只存儲(chǔ)總剛度矩陣的上三角部分，可以節(jié)約一半的存儲(chǔ)容量。

2）稀疏性。總剛度矩陣的元素絕大部分都是零元素，非零元素只占一小部分。

3）帶形分布?？倓偠染仃嚨姆橇阍胤植荚谝詫?duì)角線為中心的帶形區(qū)域內(nèi)，在半個(gè)帶形區(qū)域內(nèi)，包括對(duì)角線元素，每行具有的元素個(gè)數(shù)叫做半帶寬，即可以只存儲(chǔ)該部分元素值來(lái)節(jié)約存儲(chǔ)空間。

若存儲(chǔ)時(shí)不采取有效的存儲(chǔ)方式，不僅需要大量的存儲(chǔ)空間、時(shí)間，還會(huì)增大計(jì)算量。因此，本文在實(shí)現(xiàn)過(guò)程中根據(jù)總剛度矩陣的特性采用了較為常用的CSR矩陣壓縮存儲(chǔ)格式，用3個(gè)一維數(shù)組分別存儲(chǔ)稀疏矩陣中上三角的半帶寬非零元素值、對(duì)應(yīng)的列號(hào)及行偏移，使得GPU并行更易實(shí)現(xiàn)。

在并行實(shí)現(xiàn)時(shí)，GPU運(yùn)行kernel核函數(shù)，每個(gè)CUDA線程執(zhí)行一個(gè)網(wǎng)格單元?jiǎng)偠染仃噷懭肟倓偠染仃嚨娜蝿?wù)，主要分為如下步驟：

步驟1根據(jù)參數(shù)計(jì)算出網(wǎng)格單元的剛度矩陣。

步驟2根據(jù)網(wǎng)格單元的結(jié)點(diǎn)編號(hào)，判斷網(wǎng)格單元是上三角形單元還是下三角形單元。

步驟3根據(jù)判斷結(jié)果將值寫入總剛度矩陣對(duì)應(yīng)的位置。

為能夠并行計(jì)算剛度矩陣和載荷向量并且不產(chǎn)生線程和存儲(chǔ)沖突，本文根據(jù)不同單元網(wǎng)格結(jié)點(diǎn)編號(hào)進(jìn)行固定空間的索引存取，即每次同時(shí)迭加所有單元的單元?jiǎng)偠染仃囃晃恢玫脑亍？倓偠染仃嚳傒d荷向量和結(jié)點(diǎn)編號(hào)密切相關(guān)，因此在GPU并行求解總剛度矩陣和總載荷向量時(shí)，最為關(guān)鍵的是如何將各結(jié)點(diǎn)寫入到GPU存儲(chǔ)空間。對(duì)于本次求解的網(wǎng)格單元，每個(gè)結(jié)點(diǎn)的載荷在F中最多會(huì)迭加6次。為避免同時(shí)迭加或者使用CUDA原子操作，CUDA編程時(shí)為每個(gè)結(jié)點(diǎn)分配一個(gè)大小為6的空間，即最開始的總載荷向量F大小為N×6。實(shí)現(xiàn)中要保證下三角結(jié)點(diǎn)的編號(hào)順序始終為(m，m ＋x＋1，m ＋x），上三角順序(m＋x＋1，m，m＋1），把所有下三角的結(jié)點(diǎn)編號(hào)按順序設(shè)為索引1、2、3，同理把所有上三角的結(jié)點(diǎn)依次設(shè)為4、5、6，如圖8所示。相同索引在同一結(jié)點(diǎn)處只會(huì)出現(xiàn)一次或者不出現(xiàn)，因此在求得結(jié)點(diǎn)之后先將載荷向量寫入對(duì)應(yīng)結(jié)點(diǎn)對(duì)應(yīng)的索引處，如對(duì)于m點(diǎn)，對(duì)應(yīng)的位置為6×m＋I(xiàn)ndex(索引），這樣即可線程和存儲(chǔ)不沖突，從而正確求解。

圖8 結(jié)點(diǎn)存儲(chǔ)位置Fig.8 Node storage location

當(dāng)所有值全部寫入上述數(shù)組之后，則每個(gè)結(jié)點(diǎn)處真正的載荷值為6個(gè)值相加，此時(shí)每個(gè)CUDA線程只需要將6個(gè)數(shù)據(jù)相加便可得到最后的總剛度矩陣和總載荷向量。

3 實(shí)驗(yàn)驗(yàn)證

3.1 實(shí)驗(yàn)環(huán)境

運(yùn)用上述串并行實(shí)現(xiàn)原理方法，進(jìn)行了密度場(chǎng)有限元法求解。實(shí)驗(yàn)平臺(tái)操作系統(tǒng)采用Ubuntu 16.04，CPU為Intel i7-6800k，GPU為NVIDIA GTX 1080Ti。實(shí)驗(yàn)數(shù)據(jù)由某風(fēng)洞研究所提供。對(duì)比串行和并行實(shí)現(xiàn)下不同的網(wǎng)格精度、各主要模塊的計(jì)算耗時(shí)來(lái)驗(yàn)證并行快速有限元法的有效性。

3.2 串行速度

表1列出了不同精度大小的網(wǎng)格單元剖分下主要模塊串行求解時(shí)間。圖9展示了有限元法同網(wǎng)格精度下主要模塊占總運(yùn)行時(shí)間的比例。根據(jù)圖9結(jié)果進(jìn)行分析，在不同網(wǎng)格規(guī)模下，有限元法求解密度場(chǎng)過(guò)程中網(wǎng)絡(luò)訓(xùn)練和網(wǎng)絡(luò)預(yù)測(cè)時(shí)間開銷相差較小，但是隨著網(wǎng)格規(guī)模加大，求解總剛度矩陣和總載荷向量的時(shí)間開銷卻成倍增加，從求解總時(shí)間來(lái)看這整個(gè)過(guò)程極其耗時(shí)，完全無(wú)法達(dá)到風(fēng)洞試驗(yàn)分析場(chǎng)密度實(shí)時(shí)要求。

表1 不同模塊下CPU串行求解時(shí)間Table 1 CPU serial solution time under different modules

圖9 不同網(wǎng)格精度下主要模塊占總運(yùn)行時(shí)間的比例Fig.9 Proportion ofmain modules in total running time with different grid precision

3.3 并行速度

為了解決不能實(shí)時(shí)計(jì)算密度場(chǎng)的問題，本文著力于尋求一種基于GPU的快速求解方法，使用現(xiàn)在流行的GPU高效計(jì)算硬件平臺(tái)，極大提高了求解效率。分別進(jìn)行了如表2所示幾組實(shí)驗(yàn)。圖10展示了有限元法不同網(wǎng)格精度下串行求解和并行求解密度場(chǎng)的加速比；圖11展示了有限元法不同網(wǎng)格精度下串行和并行求解密度場(chǎng)主要模塊的加速比。本文所指的加速比皆是指串行求解和并行求解運(yùn)行時(shí)間的比值。

圖10 不同網(wǎng)格精度下求解密度場(chǎng)的加速比Fig.10 Acceleration ratio of solved density field with different grid precision

圖11 不同網(wǎng)格精度下求解密度場(chǎng)主要模塊的加速比Fig.11 Acceleration ratio of main modules of solved density field with different grid precision

表2 不同模塊下GPU并行求解時(shí)間Table 2 GPU parallel solving time under different modules

從不同的網(wǎng)格劃分疏密程度的實(shí)驗(yàn)結(jié)果對(duì)比可以看到，剖分的網(wǎng)格較為稀疏時(shí)求解密度場(chǎng)總體加速比較小，大概在25.6倍左右，隨著剖分網(wǎng)格越密，總體加速比都在逐漸增加，這是由于網(wǎng)格精度增加，有限元法產(chǎn)生的稀疏矩陣規(guī)模越大，利用GPU并行求解總剛度矩陣和總載荷向量的優(yōu)勢(shì)越明顯，而網(wǎng)絡(luò)訓(xùn)練的加速比改變不大是因?yàn)閷?dǎo)入的實(shí)驗(yàn)數(shù)據(jù)量并沒改變。另外，可以看到求解過(guò)程中的求解總剛度矩陣和總載荷向量占總運(yùn)行時(shí)間的比例逐漸上升，導(dǎo)致這一現(xiàn)象的主要原因是隨著網(wǎng)格數(shù)量增加，總剛度矩陣的組裝和總載荷向量的累加操作需要大量的共享存儲(chǔ)器，每個(gè)流處理器同時(shí)運(yùn)行的block快速減少，從而GPU的利用率下降，運(yùn)行時(shí)間增加。

3.4 精度比較

本文中采用L2范數(shù)和L∞范數(shù)來(lái)計(jì)算誤差，計(jì)算如式(20）和式(21）所示，分別表示本文快速有限元法數(shù)值計(jì)算結(jié)果密度xi與精確結(jié)果x＾的均方誤差，以及本文快速有限元法數(shù)值計(jì)算結(jié)果密度xi與精確結(jié)果絕對(duì)值的最大誤差。

式中：Δh為剖分的網(wǎng)格單元面積。

從圖12可以看出，有限元法串行和并行求解密度場(chǎng)的精度隨著網(wǎng)格劃分變密其精度增加，但精度增加緩慢，且串行計(jì)算精度高于并行計(jì)算，這是因?yàn)椴⑿杏?jì)算加入了一些累加計(jì)算，所以截?cái)嗾`差影響更顯著。另外，該方法精度在絕大多數(shù)應(yīng)用場(chǎng)景中都處于可接受水平，如在網(wǎng)格精度為496×450時(shí)，其均方誤差和最大絕對(duì)值誤差均在10－5量級(jí)下，在絕大部分場(chǎng)景下，這一誤差可以完全忽略。

圖12 有限元串并行精度比較Fig.12 Serial and parallel precision comparison of finite element method

4 結(jié)論

1）給出了基于GPU實(shí)現(xiàn)的快速有限元求解風(fēng)洞試驗(yàn)密度場(chǎng)的方法。通過(guò)對(duì)比實(shí)驗(yàn)分析對(duì)串行方法中耗時(shí)較多的神經(jīng)網(wǎng)絡(luò)擬合、總剛度矩陣組裝及總載荷向量求解實(shí)現(xiàn)了GPU并行加速。實(shí)驗(yàn)表明，隨著網(wǎng)格數(shù)量的增加，其加速比可達(dá)數(shù)百倍以上，且其精度達(dá)到了工業(yè)要求，提高了有限元法求解密度場(chǎng)的時(shí)效性。

2）實(shí)現(xiàn)的快速有限元法相對(duì)于傳統(tǒng)的有限元簡(jiǎn)化了對(duì)剛度矩陣的求解，大大減少了計(jì)算時(shí)間。

3）隨著網(wǎng)格數(shù)量增加，加速比相應(yīng)增加，其精度也逐步提高，但求解更耗時(shí)。