吳超
摘 要:針對(duì)連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化存在的計(jì)算量大、計(jì)算效率低等問題,開展了基于GPU并行計(jì)算的大規(guī)模結(jié)構(gòu)靜力拓?fù)鋬?yōu)化方法研究。首先,為了減少有限元分析的迭代次數(shù),引入了雅可比(Jacobi)對(duì)角線預(yù)處理器,研究基于共軛梯度法和預(yù)處理技術(shù)的結(jié)構(gòu)有限元并行計(jì)算方法。其次,基于單元免組裝技術(shù),結(jié)合并行迭代計(jì)算方法,研究基于GPU的結(jié)構(gòu)靜力拓?fù)鋬?yōu)化并行計(jì)算方法。在完成上述方法的Matlab和C++并行計(jì)算核函數(shù)編程后,進(jìn)行了大量的算例考核。通過給出的算例來驗(yàn)證提出方法的有效性和計(jì)算效率,結(jié)果表明,該方法具有重要的理論價(jià)值和工程應(yīng)用前景。
關(guān)鍵詞:拓?fù)鋬?yōu)化;GPU并行;免組裝方法;共軛梯度法;預(yù)處理器
中圖分類號(hào):TB34 ? ? 文獻(xiàn)標(biāo)志碼:A ? ? 文章編號(hào):1003-5168(2022)10-0011-05
DOI:10.19968/j.cnki.hnkj.1003-5168.2022.10.002
GPU-Based Structural Static Topology Optimization Design Method
WU Chao
(School of Automotive and Mechanical Engineering,Changsha University of Science and Technology,Changsha? 410014, China)
Abstract:Aiming at the problems of large amount of calculation and low efficiency in topology optimization of continuum structures, the static topology optimization method of large-scale structures based on GPU parallel computing is studied.Firstly,in order to reduce the number of iterations of finite element analysis,Jacobi diagonal preprocessor are introduced to study the parallel calculation method of structural finite element based on conjugate gradient method and preprocessing technology.Secondly,based on the element free assembly technology and the parallel iterative calculation method,the parallel calculation method of structural static topology optimization based on GPU is studied.Completed the Matlab and C++parallel computing kernel function programming of the above method,and completed a large number of numerical examples.Numerical examples are given to verify the effectiveness and computational efficiency of the proposed method,and show that the proposed method has important theoretical value and engineering application prospect.
Keywords:topology optimization;GPU parallel;assembly free method;conjugate gradient method;preprocessor
0 引言
提高大規(guī)模工程結(jié)構(gòu)優(yōu)化設(shè)計(jì)的計(jì)算效率的常用方法有兩種。一是改進(jìn)求解算法的效率,二是采用并行計(jì)算替代傳統(tǒng)串行計(jì)算??茖W(xué)計(jì)算中的單核(串行)計(jì)算的執(zhí)行效率低,其對(duì)大規(guī)模計(jì)算通常需要數(shù)天甚至數(shù)周才能執(zhí)行成功。因此,在計(jì)算中更多的是使用多核(并行)操作。由于堆疊更高計(jì)算力導(dǎo)致CPU硬件成本增高,且計(jì)算效率依舊差強(qiáng)人意,甚至?xí)母嗟哪芰?。因此,基于CPU的結(jié)構(gòu)計(jì)算與優(yōu)化很難推廣應(yīng)用到大規(guī)模工程結(jié)構(gòu)的設(shè)計(jì)中。自從NVIDIA公司成功研發(fā)出圖形處理單元(Graphics Processing Units,GPU)的硬件和統(tǒng)一計(jì)算設(shè)備架構(gòu)(Compute Unified Device Architecture,CUDA)的軟件架構(gòu)[1],利用GPU在計(jì)算上的巨大潛力,越來越多的科學(xué)工作者在工程分析和數(shù)值計(jì)算中使用GPU技術(shù)[2]。由于GPU擁有數(shù)千個(gè)內(nèi)核(多核)用來并行執(zhí)行數(shù)百萬個(gè)線程[3],所以GPU作為一種高性能的計(jì)算平臺(tái),以合理的成本擁有巨大的計(jì)算能力和更大的內(nèi)存帶寬,從而引起更多的關(guān)注。綜上所述,為了能夠解決大規(guī)模拓?fù)鋬?yōu)化計(jì)算效率問題,使用基于CUDA kernel函數(shù)的GPU并行計(jì)算方法,并結(jié)合預(yù)處理共軛梯度法的有限元迭代算法,提出一種基于GPU的并行結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)方法。
1 CUDA kernel函數(shù)并行計(jì)算方法
CUDA編程模型[4]即各CUDA線程在GPU中各自獨(dú)立且互不影響地運(yùn)行,而CPU與之協(xié)同計(jì)算。這表明CUDA kernel內(nèi)核在GPU上運(yùn)行,其他不能并行操作部分則在CPU上運(yùn)行,且CUDA編程模型中每個(gè)線程都有自己獨(dú)有的局部?jī)?nèi)存,每個(gè)線程塊都有對(duì)該塊所有線程的可見共享內(nèi)存,并且所有線程都可訪問所有線程共有的全局內(nèi)存?;镜腉PU模型執(zhí)行步驟如下。
①在設(shè)備(GPU)上分配內(nèi)存。
②將數(shù)據(jù)從主機(jī)(CPU)傳輸?shù)皆O(shè)備。
③內(nèi)核啟動(dòng)并行執(zhí)行。
④復(fù)制結(jié)果到主機(jī)。
在Matlab軟件中,通過調(diào)用相關(guān)函數(shù)并綁定,可將CUDA C等高級(jí)語言編譯成為PTX格式的中間代碼,從而可以在Matlab中調(diào)用kernel函數(shù),并行計(jì)算數(shù)據(jù)。
采用PTX代碼實(shí)現(xiàn)GPU并行計(jì)算時(shí),可通過以下4個(gè)步驟進(jìn)行(具體流程見圖1)。
①用CUDA C語言來編寫具有g(shù)lobal屬性的kernel函數(shù),并將其保存為.cu格式的文件。其中,同一個(gè).cu格式的文件中可以包含多個(gè)kernel函數(shù)。
②調(diào)用nvcc編譯器將其轉(zhuǎn)換成.ptx格式的文件。
③在Matlab中用parallel.gpu.CUDAKernel(*.ptx,*.cu)命令綁定.ptx文件中需要使用的kenrel核函數(shù)對(duì)象,并指定執(zhí)行函數(shù)的線程塊和線程參數(shù)。
④在Matlab中通過feval函數(shù)調(diào)用子程序的方式進(jìn)行并行計(jì)算。
以有限元計(jì)算為例,由于每個(gè)單元均在自己的數(shù)據(jù)空間內(nèi)獨(dú)立地進(jìn)行相同計(jì)算。因此,編寫kernel核函數(shù)時(shí)可采用一個(gè)CUDA線程對(duì)一個(gè)單元進(jìn)行計(jì)算的策略,實(shí)現(xiàn)循環(huán)展開,從而實(shí)現(xiàn)對(duì)所有單元并行、無序的計(jì)算,實(shí)現(xiàn)提高計(jì)算效率的目的。
2 預(yù)處理共軛梯度法
共軛梯度方法是一種用于求解大型線性方程組的迭代算法,它是一種介于最速下降法和牛頓法之間的方法。對(duì)于大規(guī)模拓?fù)鋬?yōu)化中的大型有限元計(jì)算問題,由于直接法不適合并行計(jì)算,共軛梯度法以增加求解步驟的處理時(shí)間為代價(jià),來減少矩陣求逆操作的內(nèi)存需求。由線性彈性問題離散化產(chǎn)生的線性方程組見式(1)[5]。
[Ku=f]? ? ? ? (1)
式中:[K]是總體剛度矩陣;[u]是位移向量;[f]是節(jié)點(diǎn)力向量。
在計(jì)算機(jī)中進(jìn)行計(jì)算時(shí),由于舍入誤差不存在,使得每次迭代后的殘差[ri]不能保持正交,從而使收斂速度降低。另外,當(dāng)采用懲罰函數(shù)處理邊界條件時(shí),得到的線性方程組的系數(shù)矩陣通常具有很大的條件數(shù),使其非常病態(tài),從而導(dǎo)致迭代收斂的速度非常緩慢。尤其是對(duì)大型、超大型方程組而言,通常需要進(jìn)行[n]次迭代才能收斂,這是實(shí)際計(jì)算所不能接受的。目前,有很多學(xué)者采用預(yù)處理共軛梯度法(PCG)以確保能夠以更快的速度收斂。通過引入預(yù)處理矩陣[M],對(duì)系數(shù)矩陣進(jìn)行適當(dāng)變換,以降低其條件數(shù),從而達(dá)到加速收斂的目的[6]。
預(yù)處理共軛梯度法(PCG)通過引入一個(gè)對(duì)稱正定矩陣[M=WTW],將原始線性系統(tǒng)式(1)替換為等效系統(tǒng)式(2)。
[M-1Ku=M-1f]? ? ?(2)
式中:[M-1]為預(yù)處理矩陣的逆。
雅可比對(duì)角線預(yù)處理矩陣是迄今為止實(shí)現(xiàn)最簡(jiǎn)單、計(jì)算成本最低的預(yù)處理矩陣。在并行環(huán)境中,與其他預(yù)處理矩陣相比,對(duì)角線預(yù)處理矩陣具有較強(qiáng)的競(jìng)爭(zhēng)力。雅可比預(yù)處理矩陣[M]實(shí)際上是由剛度矩陣[K]的對(duì)角線元素組成的。由于雅可比預(yù)處理器是對(duì)角矩陣,所以其逆也是對(duì)角的。雅可比對(duì)角線預(yù)處理矩陣的表示見式(3)。
[M=diag[k11,k22,...,knn]]? ? (3)
式中:[k11,k22,...,knn]為剛度矩陣[K]的對(duì)角線元素;[n]為剛度矩陣[K]的維度。
在共軛梯度法中,[ui]表示第[i]步優(yōu)化時(shí)更新的位移向量,則負(fù)梯度表示見式(4)。
[ri=-(Kui-f)]? ? ? (4)
式中:[ri]表示第[i]步優(yōu)化更新時(shí)的負(fù)梯度方向,梯度[ri]也可稱作每一步的殘差向量,誤差定義為[ui]與最優(yōu)解[u*]之間的差值[ei],即[ei=u*-][ui]。
在預(yù)處理共軛梯度法中,引入一個(gè)新的向量[zi],在每次計(jì)算出殘差向量[ri]以后,增加了一個(gè)求解步驟,見式(5)。
[zi=M-1ri]? ? ? ?(5)
式中:[zi]表示第[i]步優(yōu)化更新時(shí)經(jīng)過預(yù)處理的殘差向量。
在共軛梯度法中,需要解決搜索方向[pi]和優(yōu)化步長(zhǎng)[αi]這兩個(gè)問題。假定每一步的優(yōu)化方向用[pi]表示(也稱為搜索方向)。假設(shè)初始搜索方向?yàn)槌跏钾?fù)梯度方向,則經(jīng)過預(yù)處理的搜索方向[p0]見式(6)。
[p0=z0=M-1(f-Ku0)]? ? (6)
式中:[p0]表示初始的搜索方向;一般取[u0]為零向量。
在共軛梯度法中,有推論1:第[i]步計(jì)算的梯度[ri]和前[i-1]步的搜索向量[{pk}i-1k=1]正交。
根據(jù)此推論,在以后的每次迭代中,優(yōu)化步長(zhǎng)[αi]更新為式(7)。
[αi=rTizipTiKpi]? ? ? ?(7)
式中:[αi]表示第[i]步優(yōu)化更新時(shí)的優(yōu)化步長(zhǎng);[pi]表示第[i]步優(yōu)化更新時(shí)的搜索方向。
在共軛梯度法中,有推論2:第[i]步計(jì)算的梯度[ri]和前[i-2]步的搜索向量[{pk}i-2k=1]共軛正交。
根據(jù)此推論,每次迭代更新的搜索方向[pi+1]更新為式(8)。
[pi+1=zi+1+βipi]? ? ?(8)
式中:第[i]步優(yōu)化更新用于更新搜索方向[pi+1]的內(nèi)積[βi=rTi+1zi+1rTizi]。
預(yù)處理共軛梯度法的算法流程如下。
①設(shè)置剛度矩陣[K],設(shè)置初始猜測(cè)[u0],右端項(xiàng)為[f],最大迭代次數(shù)為[imax]。
②[r0=f-Ku0],[z0=M-1r0],[p0=z0]。
③對(duì)[i=1,2,...,imax]進(jìn)行如下迭代。
[αi=rTizipTiKpi]
[ui+1=ui+αipi]
[ri+1=ri+αiKpi]
[zi+1=M-1ri+1]
[βi=rTi+1zi+1rTizi]
[pi+1=zi+1+βipi]
④當(dāng)[ri2<ε]([ε]為設(shè)置的精度要求)時(shí),結(jié)束迭代。
本研究使用前文所述的預(yù)處理共軛梯度法和GPU并行法對(duì)大規(guī)模拓?fù)湓O(shè)計(jì)進(jìn)行優(yōu)化,減少有限元求解的迭代次數(shù),提高計(jì)算效率。
3 基于GPU拓?fù)鋬?yōu)化的關(guān)鍵技術(shù)
在共軛梯度法計(jì)算過程中,大部分計(jì)算時(shí)間都花費(fèi)在對(duì)矩陣向量積(SpMV)的求解上。在共軛梯度法中只有一次對(duì)稀疏矩陣向量積(SpMV)的運(yùn)算,其他部分的運(yùn)算(包括計(jì)算內(nèi)積、標(biāo)量向量積、向量加法)都可以快速求解,可以使用Matlab直接進(jìn)行串行計(jì)算。
稀疏矩陣向量積的運(yùn)算,即求解共軛梯度法中優(yōu)化步長(zhǎng)[αi]中的[Kpi]。在大規(guī)模系統(tǒng)中,整體剛度矩陣規(guī)模過大,會(huì)大幅降低算法的運(yùn)算速度。本研究采用了一種基于有限元計(jì)算的逐單元(EBE)免組裝方案,包括邊界條件和各種載荷的處理。該策略消除了組裝過程在單元級(jí)求解線性方程組,從而避免生成總體矩陣及稀疏矩陣格式轉(zhuǎn)換的開銷。在免組裝解決方案中,共軛梯度法中的計(jì)算密集型矩陣向量積SpMV被較小的單元級(jí)密集MvP所取代,每個(gè)線程僅使用單元?jiǎng)偠染仃嘯K(e)]或其約束后的變體來進(jìn)行并行計(jì)算,因此計(jì)算時(shí)間將大幅度縮短[7]。計(jì)算公式見式(9)。
[b=e?εK(e)p(e)]? ? ? (9)
式中:[b]表示[Kpi]的矩陣向量積運(yùn)算;[ε]是單元的集合。
矩陣向量積(SpMV)即計(jì)算[b=Kp]的CUDA kernel函數(shù)的基本算法如下。
①將線程分配到每個(gè)單元上。
②設(shè)id=blockDim.x×blockIdx.x+threadIdx.x,用來固定每個(gè)單元的線程索引號(hào)。
③For與該單元相關(guān)聯(lián)的單元[e](各單元的線程開始同時(shí)計(jì)算);
讀取單元?jiǎng)偠染仃嘯K(e)]和其節(jié)點(diǎn)的位置信息;
計(jì)算[be=K(e)p(e)];
累加[be]到[b]。
④End。
⑤同步線程。
此外,在本研究設(shè)計(jì)的基于GPU和預(yù)處理共軛梯度法的優(yōu)化拓?fù)浣Y(jié)構(gòu)中,CUDA進(jìn)行原子操作過程中,在并行更新結(jié)果向量時(shí)可能會(huì)發(fā)生大量?jī)?nèi)存沖突。此外,GPU并行運(yùn)算僅支持單精度的CUDA原子,不支持雙精度的CUDA原子,所以會(huì)影響計(jì)算結(jié)果的精確度。為了解決CUDA原子操作的缺點(diǎn),本研究采用兩階段方法[8],并對(duì)拓?fù)鋬?yōu)化過程中的預(yù)處理矩陣的形成、柔順度計(jì)算、靈敏度計(jì)算、靈敏度過濾均使用CUDA kernel函數(shù)編寫的并行程序進(jìn)行處理。
4 優(yōu)化模型
基于GPU并行和預(yù)處理共軛梯度法的拓?fù)鋬?yōu)化是以結(jié)構(gòu)柔順度最小為目標(biāo)的函數(shù),采用SIMP材料懲罰模型和OC準(zhǔn)則法,包含體積約束的結(jié)構(gòu)拓?fù)鋬?yōu)化模型,見式(10)[9]。
[xmin:c(x)=uTKu=e=1N(xe)puTekues.t.:V(x)V0=vol? ? ? ? Ku=f? ? ? ? 0<xmin≤xi≤1? ? ?i=1,2,...,N]? (10)
式中:[xi]表示第[i]個(gè)單元的設(shè)計(jì)變量;[xmin]為最小設(shè)計(jì)變量(非零以避免奇異);[V(x)]和[V0]分別為材料體積設(shè)計(jì)域體積;[vol]是規(guī)定的體積分?jǐn)?shù)。
本研究使用的是SIMP材料懲罰模型,其關(guān)于設(shè)計(jì)變量的表達(dá)式見式(11)。
[Ee(xe)=(xe)pE0e]? ? ?(11)
式中:[E0e]為初始材料彈性模量;[p]為懲罰因子;[Ee]表示第[e]個(gè)單元的設(shè)計(jì)變量所對(duì)應(yīng)的單元彈性模量。經(jīng)過有限元離散后的結(jié)構(gòu),考慮到單元?jiǎng)偠染仃嚺c彈性模量具有簡(jiǎn)單的線性關(guān)系,SIMP材料懲罰模型也可以表示為[ke=(xe)αkk0e]。
5 優(yōu)化算例
為了驗(yàn)證本研究所提出的優(yōu)化方法,通過上述方法,對(duì)圖2所示的MBB梁結(jié)構(gòu)以柔順度最小為目標(biāo)進(jìn)行拓?fù)鋬?yōu)化求解,并使用式(10)的優(yōu)化模型。
MBB梁結(jié)構(gòu)左右兩側(cè)下端兩點(diǎn)支撐在滾輪上,本研究針對(duì)MBB梁結(jié)構(gòu)的一半進(jìn)行優(yōu)化設(shè)計(jì)。在頂部中點(diǎn)沿豎直向下方向施加集中豎向載荷1 N。MBB梁結(jié)構(gòu)初始設(shè)計(jì)域的尺寸為4 m×1 m×0.004 m,其初始材料的彈性模量為E=2×1011 Pa,其泊松比為v=0.3。將圖2所示設(shè)計(jì)域離散為800×200=160 000個(gè)四節(jié)點(diǎn)矩形平面應(yīng)力單元。設(shè)置密度過濾半徑為[rmin=2D],其中[D]為最大單元邊長(zhǎng)。在此示例中,目標(biāo)體積與設(shè)計(jì)域初始體積的比值設(shè)定為0.4。
本算例采用OC算法進(jìn)行求解,其中[Xmin]設(shè)置為0.001,其他參數(shù)根據(jù)本研究中相應(yīng)章節(jié)的建議值選取。
圖3為使用了本研究的預(yù)處理共軛梯度法,并且通過GPU并行對(duì)該算例進(jìn)行拓?fù)鋬?yōu)化的拓?fù)錁?gòu)型,其優(yōu)化拓?fù)錁?gòu)型的柔順度為5.8×105 N·m。
圖4為MBB梁在優(yōu)化拓?fù)錁?gòu)型的柔順度優(yōu)化歷程。其中,柔順度在優(yōu)化前期下降速度較快,在優(yōu)化后期趨于穩(wěn)定。
由于該算例的結(jié)構(gòu)規(guī)模巨大,有十幾萬個(gè)單元的網(wǎng)格劃分,如果使用基于串行計(jì)算的大規(guī)模結(jié)構(gòu)拓?fù)鋬?yōu)化方法[10],會(huì)導(dǎo)致其優(yōu)化速度過于緩慢。經(jīng)過試驗(yàn)驗(yàn)證,使用傳統(tǒng)的串行計(jì)算導(dǎo)致每一次優(yōu)化迭代的時(shí)間在30 min以上,從而導(dǎo)致該算例的結(jié)構(gòu)拓?fù)鋬?yōu)化難以進(jìn)行。而本研究使用共軛梯度法進(jìn)行預(yù)處理,并用GPU進(jìn)行計(jì)算,使每次優(yōu)化迭代步小于1 min,從而提高了幾十倍的計(jì)算速度,所以對(duì)大規(guī)模拓?fù)鋬?yōu)化使用此方法是有效的,在很大程度上減少計(jì)算量,縮短計(jì)算時(shí)間。
6 結(jié)語
本研究基于GPU并行技術(shù),提出了一種使用預(yù)處理共軛梯度迭代法的連續(xù)體靜力拓?fù)鋬?yōu)化方法。并通過一系列具體算例驗(yàn)證了所提方法的可行性、有效性和其速度優(yōu)勢(shì)??梢缘贸鲆韵陆Y(jié)論:①以Matlab為基礎(chǔ),通過調(diào)用CUDA C編寫的kernel核函數(shù)來實(shí)現(xiàn)拓?fù)鋬?yōu)化的并行計(jì)算;②為了便于GPU的并行計(jì)算,本研究使用迭代算法,且為了減少迭代法的結(jié)構(gòu)有限元分析迭代次數(shù),引入了雅可比(Jacobi)對(duì)角線預(yù)處理器,研究了基于共軛梯度法和預(yù)處理技術(shù)的結(jié)構(gòu)有限元并行計(jì)算方法;③本研究開發(fā)的GPU并行計(jì)算方法來處理基于迭代法的有限元計(jì)算中的約束條件、矩陣向量積、預(yù)處理矩陣和拓?fù)鋬?yōu)化過程中的柔順度、靈敏度和靈敏度過濾等;④完成上述方法的Matlab和C++并行計(jì)算核函數(shù)編程,提出了基于GPU和預(yù)處理共軛梯度法的結(jié)構(gòu)拓?fù)鋬?yōu)化方法,并完成了算例驗(yàn)證,驗(yàn)證了本研究提出方法的有效性和計(jì)算效率。
綜上所述,本研究所提出的方法有待于推廣到基于GPU并行計(jì)算的多工況和多約束的大規(guī)模三維結(jié)構(gòu)拓?fù)鋬?yōu)化問題。
參考文獻(xiàn):
[1] PIKLE N K,SATHE S R,VYAVAHARE A Y.Low occupancy high performance elemental products in assembly free FEM on GPU[J].Engineering with Computers,2021(5).
[2] 張朝暉,劉俊起,徐勤建.GPU并行計(jì)算技術(shù)分析與應(yīng)用[J].信息技術(shù),2009(11):86-89.
[3] UTPAL K,DEEPAK S,SINGH G S.GPU-Warp based Finite Element Matrices Generation and Assembly using Coloring Method[J].Journal of Computational Design and Engineering,2018(4):4.
[4] KOMATITSCH D,MICHéA D,ERLEBACHER G.Porting a high-order finite-element earthquake modeling application to NVIDIA graphics cards using CUDA[J].Journal of Parallel and Distributed Computing,2009(5):451-460.
[5] 陳成,趙圣佞.基于Heaviside過濾和可行域調(diào)整的SIMP方法拓?fù)鋬?yōu)化設(shè)計(jì)[J].河南科技,2018(34):26-28.
[6] BOLZ J,F(xiàn)ARMER I,GRINSPUN E,et al.Sparse matrix solvers on the GPU: conjugate gradients and multigrid[J].ACM Transactions on Graphics,2003(3):917-924.
[7] 榮見華.漸進(jìn)結(jié)構(gòu)優(yōu)化方法及其應(yīng)用研究[D].長(zhǎng)沙:國(guó)防科學(xué)技術(shù)大學(xué),2006.
[8] PIKLE N K,SATHE S R,VYAVAHARE A Y.High performance iterative elemental product strategy in assembly-free FEM on GPU with improved occupancy[J].Computing,2018(12):1-25.
[9] 榮見華,彭羅,易繼軍,等.一種新的多輸入多輸出柔順機(jī)構(gòu)拓?fù)鋬?yōu)化方法[J].長(zhǎng)沙理工大學(xué)學(xué)報(bào)(自然科學(xué)版),2021(1):66-78.
[10] SIGMUND O.A 99 line topology optimization code written in Matlab[J].Structural & Multidisciplinary Optimization,2001(2):120-127.