無網(wǎng)格法結(jié)構(gòu)拓?fù)鋬?yōu)化模型的GPU 并行加速求解及應(yīng)用

唐 芳，馮應(yīng)朗，盧海山

（1.湖南理工職業(yè)技術(shù)學(xué)院新能源學(xué)院，湖南 湘潭 411105；2.湘潭大學(xué)機(jī)械工程與力學(xué)學(xué)院，湖南 湘潭 411105）

0 引言

結(jié)構(gòu)拓?fù)鋬?yōu)化技術(shù)可在產(chǎn)品的概念設(shè)計(jì)階段提供創(chuàng)新設(shè)計(jì)方案，現(xiàn)已廣泛應(yīng)用于汽車、航空航天以及增材制造等領(lǐng)域[1，2]。目前結(jié)構(gòu)拓?fù)鋬?yōu)化一般利用有限元法進(jìn)行結(jié)構(gòu)分析，然而由于單元的存在，使得基于有限元法的結(jié)構(gòu)拓?fù)鋬?yōu)化模型容易出現(xiàn)棋盤格、單元鉸接等數(shù)值不穩(wěn)定現(xiàn)象[3]，且在處理大變形、含裂紋結(jié)構(gòu)的拓?fù)鋬?yōu)化問題時，面臨網(wǎng)格畸變、不連續(xù)等困難[4]。

無網(wǎng)格法僅需離散節(jié)點(diǎn)信息即可構(gòu)造出高階場函數(shù)，并能有效處理大變形、不連續(xù)等問題。部分學(xué)者將無網(wǎng)格法引入結(jié)構(gòu)拓?fù)鋬?yōu)化，利用無網(wǎng)格不受網(wǎng)格束縛的優(yōu)勢，有效克服了傳統(tǒng)基于有限元法的拓?fù)鋬?yōu)化模型存在的缺點(diǎn)[5]。無網(wǎng)格Galerkin（Element-Free Galerkin，EFG）法是當(dāng)前成熟且應(yīng)用廣泛的無網(wǎng)格法之一，具有收斂快、計(jì)算精度高和穩(wěn)定性好等[6]優(yōu)點(diǎn)，在結(jié)構(gòu)拓?fù)鋬?yōu)化中得到應(yīng)用[7]，有效地抑制了傳統(tǒng)拓?fù)鋬?yōu)化中所出現(xiàn)的棋盤格等數(shù)值不穩(wěn)定現(xiàn)象，同時避免了大變形結(jié)構(gòu)拓?fù)鋬?yōu)化中的網(wǎng)格畸變。

盡管基于無網(wǎng)格法的拓?fù)鋬?yōu)化模型具有上述優(yōu)勢，但由于拓?fù)鋬?yōu)化的求解需要多次迭代，重復(fù)進(jìn)行結(jié)構(gòu)分析計(jì)算，并且無網(wǎng)格法的計(jì)算量大，導(dǎo)致拓?fù)鋬?yōu)化模型的求解極其耗時，難以應(yīng)用于大規(guī)模拓?fù)鋬?yōu)化問題。

近年來，GPU 并行加速技術(shù)因其強(qiáng)大的并行計(jì)算能力在計(jì)算力學(xué)領(lǐng)域得到廣泛應(yīng)用。韓琪等[8]針對大規(guī)模拓?fù)鋬?yōu)化問題計(jì)算量大、計(jì)算效率低的問題，結(jié)合雙向漸進(jìn)結(jié)構(gòu)拓?fù)鋬?yōu)化方法與GPU 并行加速計(jì)算技術(shù)，設(shè)計(jì)了一種并行拓?fù)鋬?yōu)化方法；Xia 等[9]等提出了一種使用GPU 并行策略與等幾何分析的水平集拓?fù)鋬?yōu)化方法，加速比達(dá)到兩個數(shù)量級。而在無網(wǎng)格法的GPU 并行加速研究方面，龔曙光等[10]通過并行化組裝剛度矩陣與求解離散方程，并為了充分發(fā)揮FEM 與EFG 法各自的優(yōu)勢，提出了FE-EFG 耦合法的GPU 并行加速算法。

針對無網(wǎng)格法結(jié)構(gòu)拓?fù)鋬?yōu)化模型的求解計(jì)算存在耗時長的問題，結(jié)合GPU 并行加速技術(shù)，提出了一種求解EFG 法拓?fù)鋬?yōu)化模型的高效計(jì)算方法。該研究對無網(wǎng)格法應(yīng)用于工程結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)具有重要的理論意義與參考價值。

1 EFG 法拓?fù)鋬?yōu)化模型

1.1 EFG 離散方程

位移u（x）的移動最小二乘逼近函數(shù)uh（x）為：

式中，N為節(jié)點(diǎn)數(shù)；U為節(jié)點(diǎn)位移參數(shù)列陣；Φ（x）為形函數(shù)矩陣，且其計(jì)算式為：

式中，A和H分別表示一個矩陣，并有

其中

式中，pm（xN）為基函數(shù)；w（x-xI）為權(quán)函數(shù)。本文采用線性基函數(shù)和三次樣條權(quán)函數(shù)。

利用罰函數(shù)法施加本質(zhì)邊界條件，可得線彈性靜力問題的EFG 離散方程：

式中，K為剛度矩陣；F為節(jié)點(diǎn)力向量；B為應(yīng)變矩陣；α為罰系數(shù)；D為彈性矩陣。

1.2 拓?fù)鋬?yōu)化模型及迭代更新

選擇節(jié)點(diǎn)相對密度參數(shù)ρi作為設(shè)計(jì)變量，則設(shè)計(jì)域內(nèi)任意點(diǎn)的相對密度ρ（x）可由移動最小二乘逼近得到

利用變密度法中的各向同性材料懲罰模型（SIMP）計(jì)算優(yōu)化后的材料彈性模量為

式中，E0為實(shí)體材料的彈性模量，p為懲罰因子。

以結(jié)構(gòu)柔度為目標(biāo)函數(shù)、材料體積為約束條件，建立如下拓?fù)鋬?yōu)化數(shù)學(xué)模型

式中，V0、V分別為優(yōu)化前后的結(jié)構(gòu)體積；ζ為體積保留率。

本文采用OC 準(zhǔn)則法求解上述優(yōu)化模型，其設(shè)計(jì)變量的更新格式為

式中

2 GPU 并行加速算法

2.1 CUDA 架構(gòu)

NVIDA 推出的CUDA 平臺，極大地簡化了GPU并行編程。在CUDA 平臺下，并行程序的執(zhí)行模式為：①在CPU 的內(nèi)存中準(zhǔn)備好數(shù)據(jù)，并復(fù)制到GPU的顯存中；②GPU 執(zhí)行設(shè)備端程序；③將計(jì)算結(jié)果傳回CPU 中的內(nèi)存。GPU 上執(zhí)行核函數(shù)的最小單位為線程，CUDA 將多個線程組成為一個線程塊，多個線程塊則構(gòu)成為一個線程格。線程組織及存儲器的架構(gòu)如圖1 所示。

圖1 CUDA 的線程組織及存儲架構(gòu)

2.2 剛度矩陣計(jì)算

基于交叉節(jié)點(diǎn)對（影響域有交集的兩個節(jié)點(diǎn)）的思想，提出了一種通過交叉節(jié)點(diǎn)對的循環(huán)來計(jì)算剛度矩陣的逐節(jié)點(diǎn)對法[18]。本文借助交叉節(jié)點(diǎn)對思想，結(jié)合上述CUDA 架構(gòu)特點(diǎn)，構(gòu)建了利用GPU 并行加速拓?fù)涞^程中剛度矩陣的并行計(jì)算流程，如圖2 所示。該并行流程有兩個并行層次：第一個層次是交叉節(jié)點(diǎn)對，即每個線程塊處理一個交叉節(jié)點(diǎn)對，第二層次是交叉節(jié)點(diǎn)對的公共積分點(diǎn)，即線程塊中的每個線程處理交叉節(jié)點(diǎn)對的一個公共積分點(diǎn)。與公共積分點(diǎn)相關(guān)的數(shù)據(jù)存儲于對應(yīng)線程的寄存器，與交叉節(jié)點(diǎn)對相關(guān)的數(shù)據(jù)儲存在線程塊的共享內(nèi)存，而最終的剛度矩陣儲存在GPU 的全局內(nèi)存。

圖2 剛度矩陣的GPU 并行計(jì)算

2.3 離散方程求解

結(jié)構(gòu)拓?fù)鋬?yōu)化模型的優(yōu)化迭代過程需要對EFG法形成的離散方程進(jìn)行重復(fù)求解，以獲得結(jié)構(gòu)響應(yīng)。由于直接法所需內(nèi)存大、矩陣分解耗時長，同時由于材料分布的高度非均勻性導(dǎo)致形成的EFG 剛度矩陣條件性態(tài)差，為了提高平衡方程的求解效率，本文采用雅克比預(yù)處理共軛梯度法，并結(jié)合GPU 并行加速技術(shù)求解離散方程，其計(jì)算流程如下：

其中，J為雅克比預(yù)處理矩陣。

2.4 GPU 并行加速流程

無網(wǎng)格法拓?fù)鋬?yōu)化模型的GPU 并行加速求解流程如圖3 所示。在該流程中，耗時量最大的計(jì)算剛度矩陣與求解離散方程兩個部分在GPU 并行計(jì)算，而其余耗時很少的部分則在CPU 串行計(jì)算。此外，由于每一次OC 迭代均需要利用形函數(shù)及其導(dǎo)數(shù)值重新計(jì)算剛度矩陣，為避免反復(fù)計(jì)算形函數(shù)及其導(dǎo)數(shù)值，在前處理計(jì)算部分提前完成形函數(shù)及其導(dǎo)數(shù)值的計(jì)算并存儲至GPU 全局內(nèi)存，從而進(jìn)一步降低計(jì)算耗時。

圖3 GPU 并行加速流程

3 數(shù)值算例

在算例中，材料的彈性模量為E0= 1.0，泊松比為v= 0.3。算例的運(yùn)行平臺配置參數(shù)見表1。

表1 運(yùn)行平臺配置參數(shù)

定義加速比為

式中，tCPU為CPU 串行算法的運(yùn)行耗時，tGPU為GPU 并行加速算法的運(yùn)行耗時。

3.1 二維懸臂梁

圖4 為懸臂梁模型。其中L= 90 mm。懸臂梁左邊界為固定約束，右邊中點(diǎn)承受豎直向下的集中力F=1 N。設(shè)計(jì)域的體積保留率為ζ= 0.1。

圖4 懸臂梁

懸臂梁的拓?fù)鋬?yōu)化結(jié)果如圖5 所示。從圖5 可知，CPU 串行程序的拓?fù)浣Y(jié)果與GPU 并行加速程序的拓?fù)浣Y(jié)果完全吻合，且拓?fù)浣Y(jié)果邊界清晰，無棋盤格等數(shù)值不穩(wěn)定現(xiàn)象。這表明上述無網(wǎng)格拓?fù)鋬?yōu)化模型及其GPU 并行加速求解算法是正確的。

圖5 懸臂梁的拓?fù)浣Y(jié)果

采用數(shù)目分別為1891（31×61）、7381（61×121）、29161（121×241）與115921（241×481）的四組節(jié)點(diǎn)離散該懸臂梁模型。四組節(jié)點(diǎn)規(guī)模下CPU 串行程序的各段耗時及比例見表2。由表2 可知，求解方程耗時占據(jù)了整個CPU 串行程序耗時的大部分，其次是計(jì)算剛度矩陣的耗時，而程序其余部分耗時的占比極低，表明CPU 串行算法的性能瓶頸主要在求解方程與計(jì)算剛度矩陣兩個部分。因此，利用GPU 并行加速剛度矩陣計(jì)算與方程求解，就能夠有效提高拓?fù)鋬?yōu)化模型的求解效率，且能夠避免在內(nèi)存與顯存之間反復(fù)傳輸剛度矩陣等大量的數(shù)據(jù)，從而進(jìn)一步提高拓?fù)鋬?yōu)化模型的求解性能。

表2 CPU 串行程序的各段耗時（s）及比例

四組節(jié)點(diǎn)規(guī)模下GPU 并行算法對CPU 串行算法的加速比如圖6 所示。由圖6 可知，當(dāng)節(jié)點(diǎn)規(guī)模較小時，整體加速比隨節(jié)點(diǎn)規(guī)模增大而增加，但當(dāng)節(jié)點(diǎn)規(guī)模增大到一定數(shù)量后，加速比出現(xiàn)小幅降低。這是由于GPU 的并行線程及寄存器與共享內(nèi)存等資源是有限的，即GPU 能同時并行計(jì)算的任務(wù)量有限。方程求解加速比與整體加速比的變化趨勢較為一致，但當(dāng)節(jié)點(diǎn)規(guī)模較大時，加速比增幅減小。剛度矩陣計(jì)算加速比隨節(jié)點(diǎn)規(guī)模增加呈現(xiàn)先小幅增加而后又小幅降低的趨勢。這表明求解離散方程的加速性能在整個拓?fù)鋬?yōu)化模型的加速求解過程中起主導(dǎo)作用。此外，該算例結(jié)果表明GPU 并行加速算法對于規(guī)模較大的無網(wǎng)格法拓?fù)鋬?yōu)化模型的求解具有更加顯著的加速效果。

圖6 不同節(jié)點(diǎn)數(shù)的加速比

3.2 二維曲形支架

圖7 為曲形支架模型。結(jié)構(gòu)尺寸為R= 100 mm、r= 50 mm。支架右邊頂點(diǎn)承受水平向左的集中力F=1 N，支架底邊施加固定約束。設(shè)計(jì)域的體積保留率為ζ= 0.5。采用8320 個節(jié)點(diǎn)離散該設(shè)計(jì)模型。由圖8 可知，CPU 與GPU 程序的拓?fù)鋬?yōu)化結(jié)果完全一致。

圖7 曲形支架

圖8 曲形支架的拓?fù)浣Y(jié)果

計(jì)算耗時及加速比如圖9 所示。CPU 程序總耗時高達(dá)13371.5 s，而GPU 程序總耗時僅為400.9 s，加速比達(dá)33.4，表明了本文GPU 并行加速算法的強(qiáng)大并行加速能力。

圖9 曲形支架的計(jì)算耗時及加速比

3.3 三維支撐平臺

圖10 為三維支撐平臺模型。結(jié)構(gòu)尺寸如圖11 所示。該支撐平臺由頂部平板與下方的錐形筒體組成。平臺頂面中心承受豎直向下的集中力400 kN，平臺底面均勻分布四處固定約束。該模型具有對稱性，可取整體模型的四分之一，得到平臺的優(yōu)化設(shè)計(jì)分析模型，如圖12 所示。其中，平臺頂部的平板（圖12 中的灰色區(qū)域）劃分為非設(shè)計(jì)域，在優(yōu)化過程中，該部分的材料保持不變，以便于承受載荷。設(shè)計(jì)域的體積保留率為ζ= 0.2。采用13476 個節(jié)點(diǎn)離散該設(shè)計(jì)模型。

圖10 支撐平臺

圖11 支撐平臺模型尺寸（mm）

圖12 支撐平臺設(shè)計(jì)模型

拓?fù)鋬?yōu)化結(jié)果如圖13 所示，該拓?fù)浣Y(jié)果合理、清晰，無棋盤格等現(xiàn)象。這是由于采用移動最小二成逼近構(gòu)造的密度場具有高階連續(xù)性，從而有效地抑制了棋盤格等數(shù)值不穩(wěn)定現(xiàn)象。此外由于劃分了非設(shè)計(jì)域，支撐平臺頂部的平板被保留，因此得到的最優(yōu)結(jié)果便于承受載荷，滿足使用要求。

圖13 支撐平臺的拓?fù)浣Y(jié)果

支撐平臺算例的計(jì)算耗時及加速比，如圖14 所示。CPU 程序總耗時為11235.1 s，而GPU 程序總耗時為518.0 s，加速比達(dá)21.7，表明了GPU 并行加速算法對于三維結(jié)構(gòu)拓?fù)鋬?yōu)化問題同樣具有很好的加速求解性能。

圖14 支撐平臺的計(jì)算耗時及加速比

3.4 多工況固支梁

多載荷工況拓?fù)鋬?yōu)化問題的每步OC 迭代過程均需要多次求解離散方程，以獲得不同載荷工況下的結(jié)構(gòu)響應(yīng)，這將進(jìn)一步增加計(jì)算耗時。本文利用GPU并行加速算法求解多工況固支梁優(yōu)化模型，以測試GPU 算法對多載荷工況拓?fù)鋬?yōu)化問題的加速求解性能。圖15 為二維多工況固支梁模型，結(jié)構(gòu)尺寸如圖15 所示，其中L= 90 mm。固支梁的左右兩端固定，頂邊和底邊中點(diǎn)處分別承受豎直向下與向上的集中力與。設(shè)計(jì)域的體積保留率為ζ= 0.3。采用10011 個節(jié)點(diǎn)離散該設(shè)計(jì)模型。圖16 所示的CPU 與GPU 程序的拓?fù)鋬?yōu)化結(jié)果完全吻合，進(jìn)一步驗(yàn)證了本文所建立的GPU 并行加速求解算法的正確性。

圖15 多工況固支梁

圖16 多工況固支梁的拓?fù)浣Y(jié)果

圖17 給出了多工況固支梁算例的計(jì)算耗時及加速比。CPU 程序總耗時為9383.3 s，GPU 程序總耗時為307.6 s。盡管在兩個載荷工況下每一次OC 迭代均需要求解兩次離散方程，但GPU 程序的耗時仍然遠(yuǎn)遠(yuǎn)小于CPU 程序耗時，且加速比達(dá)到了30.5。

圖17 多工況固支梁的計(jì)算耗時及加速比

4 結(jié)語

針對無網(wǎng)格法結(jié)構(gòu)拓?fù)鋬?yōu)化模型求解計(jì)算耗時長的問題，通過引入GPU 并行加速技術(shù)建立了無網(wǎng)格法拓?fù)鋬?yōu)化模型的并行加速求解算法，以充分發(fā)揮無網(wǎng)格法不受網(wǎng)格束縛、GPU 并行加速計(jì)算效率高的優(yōu)勢。經(jīng)算例分析得到了如下結(jié)論：

（1）無網(wǎng)格拓?fù)鋬?yōu)化模型的GPU 并行加速求解算法的拓?fù)浣Y(jié)果與CPU 串行算法的拓?fù)浣Y(jié)果完全吻合，所得到的最優(yōu)拓?fù)錁?gòu)型清晰，無棋盤格等數(shù)值奇異問題。同時，在結(jié)構(gòu)中劃分了非設(shè)計(jì)域，非設(shè)計(jì)域內(nèi)的材料在優(yōu)化迭代過程中均保持為實(shí)體材料，所得到的拓?fù)浣Y(jié)構(gòu)更為合理，滿足使用要求。

（2）算例的整體加速比最大可達(dá)46.9，表明所提無網(wǎng)格拓?fù)鋬?yōu)化模型的GPU 并行加速求解算法具有優(yōu)良的加速性能。且當(dāng)計(jì)算規(guī)模較大時，加速效果更加顯著，但受限于GPU 自身的計(jì)算資源，加速比存在上限。

（3）盡管多載荷工況拓?fù)鋬?yōu)化問題的每步OC 迭代均需要多次求解離散方程，但相比于CPU 串行算法，GPU 并行加速算法仍能大幅縮短拓?fù)鋬?yōu)化模型的求解耗時，表明所提出的GPU 并行加速算法對于多載荷工況的拓?fù)鋬?yōu)化問題同樣具有很高的求解效率。