無(wú)網(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)用于汽車(chē)、航空航天以及增材制造等領(lǐng)域[1，2]。目前結(jié)構(gòu)拓?fù)鋬?yōu)化一般利用有限元法進(jìn)行結(jié)構(gòu)分析，然而由于單元的存在，使得基于有限元法的結(jié)構(gòu)拓?fù)鋬?yōu)化模型容易出現(xiàn)棋盤(pán)格、單元鉸接等數(shù)值不穩(wěn)定現(xiàn)象[3]，且在處理大變形、含裂紋結(jié)構(gòu)的拓?fù)鋬?yōu)化問(wèn)題時(shí)，面臨網(wǎng)格畸變、不連續(xù)等困難[4]。

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

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

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

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

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

1.1 EFG 離散方程

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

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

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

其中

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

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

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

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

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

利用變密度法中的各向同性材料懲罰模型（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 平臺(tái)，極大地簡(jiǎn)化了GPU并行編程。在CUDA 平臺(tái)下，并行程序的執(zhí)行模式為：①在CPU 的內(nèi)存中準(zhǔn)備好數(shù)據(jù)，并復(fù)制到GPU的顯存中；②GPU 執(zhí)行設(shè)備端程序；③將計(jì)算結(jié)果傳回CPU 中的內(nèi)存。GPU 上執(zhí)行核函數(shù)的最小單位為線程，CUDA 將多個(gè)線程組成為一個(gè)線程塊，多個(gè)線程塊則構(gòu)成為一個(gè)線程格。線程組織及存儲(chǔ)器的架構(gòu)如圖1 所示。

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

2.2 剛度矩陣計(jì)算

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

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

2.3 離散方程求解

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

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

2.4 GPU 并行加速流程

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

圖3 GPU 并行加速流程

3 數(shù)值算例

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

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

定義加速比為

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

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é)果邊界清晰，無(wú)棋盤(pán)格等數(shù)值不穩(wěn)定現(xiàn)象。這表明上述無(wú)網(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 串行程序的各段耗時(shí)及比例見(jiàn)表2。由表2 可知，求解方程耗時(shí)占據(jù)了整個(gè)CPU 串行程序耗時(shí)的大部分，其次是計(jì)算剛度矩陣的耗時(shí)，而程序其余部分耗時(shí)的占比極低，表明CPU 串行算法的性能瓶頸主要在求解方程與計(jì)算剛度矩陣兩個(gè)部分。因此，利用GPU 并行加速剛度矩陣計(jì)算與方程求解，就能夠有效提高拓?fù)鋬?yōu)化模型的求解效率，且能夠避免在內(nèi)存與顯存之間反復(fù)傳輸剛度矩陣等大量的數(shù)據(jù)，從而進(jìn)一步提高拓?fù)鋬?yōu)化模型的求解性能。

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

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

圖7 曲形支架

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

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

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

3.3 三維支撐平臺(tái)

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

圖10 支撐平臺(tái)

圖11 支撐平臺(tái)模型尺寸（mm）

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

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

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

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

圖14 支撐平臺(tái)的計(jì)算耗時(shí)及加速比

3.4 多工況固支梁

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

圖15 多工況固支梁

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

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

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

4 結(jié)語(yǔ)

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

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

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

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