基于GPU的結(jié)構(gòu)靜力拓?fù)鋬?yōu)化設(shè)計方法

吳超

摘 要：針對連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化存在的計算量大、計算效率低等問題，開展了基于GPU并行計算的大規(guī)模結(jié)構(gòu)靜力拓?fù)鋬?yōu)化方法研究。首先，為了減少有限元分析的迭代次數(shù)，引入了雅可比（Jacobi）對角線預(yù)處理器，研究基于共軛梯度法和預(yù)處理技術(shù)的結(jié)構(gòu)有限元并行計算方法。其次，基于單元免組裝技術(shù)，結(jié)合并行迭代計算方法，研究基于GPU的結(jié)構(gòu)靜力拓?fù)鋬?yōu)化并行計算方法。在完成上述方法的Matlab和C++并行計算核函數(shù)編程后，進(jìn)行了大量的算例考核。通過給出的算例來驗證提出方法的有效性和計算效率，結(jié)果表明，該方法具有重要的理論價值和工程應(yīng)用前景。

關(guān)鍵詞：拓?fù)鋬?yōu)化;GPU并行;免組裝方法;共軛梯度法;預(yù)處理器

中圖分類號：TB34 ? ? 文獻(xiàn)標(biāo)志碼：A ? ? 文章編號：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ìn)求解算法的效率，二是采用并行計算替代傳統(tǒng)串行計算。科學(xué)計算中的單核（串行）計算的執(zhí)行效率低，其對大規(guī)模計算通常需要數(shù)天甚至數(shù)周才能執(zhí)行成功。因此，在計算中更多的是使用多核（并行）操作。由于堆疊更高計算力導(dǎo)致CPU硬件成本增高，且計算效率依舊差強(qiáng)人意，甚至?xí)母嗟哪芰俊Ｒ虼耍贑PU的結(jié)構(gòu)計算與優(yōu)化很難推廣應(yīng)用到大規(guī)模工程結(jié)構(gòu)的設(shè)計中。自從NVIDIA公司成功研發(fā)出圖形處理單元（Graphics Processing Units，GPU）的硬件和統(tǒng)一計算設(shè)備架構(gòu)（Compute Unified Device Architecture，CUDA）的軟件架構(gòu)[1]，利用GPU在計算上的巨大潛力，越來越多的科學(xué)工作者在工程分析和數(shù)值計算中使用GPU技術(shù)[2]。由于GPU擁有數(shù)千個內(nèi)核（多核）用來并行執(zhí)行數(shù)百萬個線程[3]，所以GPU作為一種高性能的計算平臺，以合理的成本擁有巨大的計算能力和更大的內(nèi)存帶寬，從而引起更多的關(guān)注。綜上所述，為了能夠解決大規(guī)模拓?fù)鋬?yōu)化計算效率問題，使用基于CUDA kernel函數(shù)的GPU并行計算方法，并結(jié)合預(yù)處理共軛梯度法的有限元迭代算法，提出一種基于GPU的并行結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計方法。

1 CUDA kernel函數(shù)并行計算方法

CUDA編程模型[4]即各CUDA線程在GPU中各自獨立且互不影響地運行，而CPU與之協(xié)同計算。這表明CUDA kernel內(nèi)核在GPU上運行，其他不能并行操作部分則在CPU上運行，且CUDA編程模型中每個線程都有自己獨有的局部內(nèi)存，每個線程塊都有對該塊所有線程的可見共享內(nèi)存，并且所有線程都可訪問所有線程共有的全局內(nèi)存。基本的GPU模型執(zhí)行步驟如下。

①在設(shè)備（GPU）上分配內(nèi)存。

②將數(shù)據(jù)從主機(jī)（CPU）傳輸?shù)皆O(shè)備。

③內(nèi)核啟動并行執(zhí)行。

④復(fù)制結(jié)果到主機(jī)。

在Matlab軟件中，通過調(diào)用相關(guān)函數(shù)并綁定，可將CUDA C等高級語言編譯成為PTX格式的中間代碼，從而可以在Matlab中調(diào)用kernel函數(shù)，并行計算數(shù)據(jù)。

采用PTX代碼實現(xiàn)GPU并行計算時，可通過以下4個步驟進(jìn)行（具體流程見圖1）。

①用CUDA C語言來編寫具有g(shù)lobal屬性的kernel函數(shù)，并將其保存為.cu格式的文件。其中，同一個.cu格式的文件中可以包含多個kernel函數(shù)。

②調(diào)用nvcc編譯器將其轉(zhuǎn)換成.ptx格式的文件。

③在Matlab中用parallel.gpu.CUDAKernel（*.ptx，*.cu）命令綁定.ptx文件中需要使用的kenrel核函數(shù)對象，并指定執(zhí)行函數(shù)的線程塊和線程參數(shù)。

④在Matlab中通過feval函數(shù)調(diào)用子程序的方式進(jìn)行并行計算。

以有限元計算為例，由于每個單元均在自己的數(shù)據(jù)空間內(nèi)獨立地進(jìn)行相同計算。因此，編寫kernel核函數(shù)時可采用一個CUDA線程對一個單元進(jìn)行計算的策略，實現(xiàn)循環(huán)展開，從而實現(xiàn)對所有單元并行、無序的計算，實現(xiàn)提高計算效率的目的。

2 預(yù)處理共軛梯度法

共軛梯度方法是一種用于求解大型線性方程組的迭代算法，它是一種介于最速下降法和牛頓法之間的方法。對于大規(guī)模拓?fù)鋬?yōu)化中的大型有限元計算問題，由于直接法不適合并行計算，共軛梯度法以增加求解步驟的處理時間為代價，來減少矩陣求逆操作的內(nèi)存需求。由線性彈性問題離散化產(chǎn)生的線性方程組見式（1）[5]。

[Ku=f]? ? ? ? （1）

式中：[K]是總體剛度矩陣;[u]是位移向量;[f]是節(jié)點力向量。

在計算機(jī)中進(jìn)行計算時，由于舍入誤差不存在，使得每次迭代后的殘差[ri]不能保持正交，從而使收斂速度降低。另外，當(dāng)采用懲罰函數(shù)處理邊界條件時，得到的線性方程組的系數(shù)矩陣通常具有很大的條件數(shù)，使其非常病態(tài)，從而導(dǎo)致迭代收斂的速度非常緩慢。尤其是對大型、超大型方程組而言，通常需要進(jìn)行[n]次迭代才能收斂，這是實際計算所不能接受的。目前，有很多學(xué)者采用預(yù)處理共軛梯度法（PCG）以確保能夠以更快的速度收斂。通過引入預(yù)處理矩陣[M]，對系數(shù)矩陣進(jìn)行適當(dāng)變換，以降低其條件數(shù)，從而達(dá)到加速收斂的目的[6]。

預(yù)處理共軛梯度法（PCG）通過引入一個對稱正定矩陣[M=WTW]，將原始線性系統(tǒng)式（1）替換為等效系統(tǒng)式（2）。

[M-1Ku=M-1f]? ? ?（2）

式中：[M-1]為預(yù)處理矩陣的逆。

雅可比對角線預(yù)處理矩陣是迄今為止實現(xiàn)最簡單、計算成本最低的預(yù)處理矩陣。在并行環(huán)境中，與其他預(yù)處理矩陣相比，對角線預(yù)處理矩陣具有較強(qiáng)的競爭力。雅可比預(yù)處理矩陣[M]實際上是由剛度矩陣[K]的對角線元素組成的。由于雅可比預(yù)處理器是對角矩陣，所以其逆也是對角的。雅可比對角線預(yù)處理矩陣的表示見式（3）。

[M=diag[k11，k22，...，knn]]? ? （3）

式中：[k11，k22，...，knn]為剛度矩陣[K]的對角線元素;[n]為剛度矩陣[K]的維度。

在共軛梯度法中，[ui]表示第[i]步優(yōu)化時更新的位移向量，則負(fù)梯度表示見式（4）。

[ri=-（Kui-f）]? ? ? （4）

式中：[ri]表示第[i]步優(yōu)化更新時的負(fù)梯度方向，梯度[ri]也可稱作每一步的殘差向量，誤差定義為[ui]與最優(yōu)解[u*]之間的差值[ei]，即[ei=u*-][ui]。

在預(yù)處理共軛梯度法中，引入一個新的向量[zi]，在每次計算出殘差向量[ri]以后，增加了一個求解步驟，見式（5）。

[zi=M-1ri]? ? ? ?（5）

式中：[zi]表示第[i]步優(yōu)化更新時經(jīng)過預(yù)處理的殘差向量。

在共軛梯度法中，需要解決搜索方向[pi]和優(yōu)化步長[αi]這兩個問題。假定每一步的優(yōu)化方向用[pi]表示（也稱為搜索方向）。假設(shè)初始搜索方向為初始負(fù)梯度方向，則經(jīng)過預(yù)處理的搜索方向[p0]見式（6）。

[p0=z0=M-1（f-Ku0）]? ? （6）

式中：[p0]表示初始的搜索方向;一般取[u0]為零向量。

在共軛梯度法中，有推論1：第[i]步計算的梯度[ri]和前[i-1]步的搜索向量[{pk}i-1k=1]正交。

根據(jù)此推論，在以后的每次迭代中，優(yōu)化步長[αi]更新為式（7）。

[αi=rTizipTiKpi]? ? ? ?（7）

式中：[αi]表示第[i]步優(yōu)化更新時的優(yōu)化步長;[pi]表示第[i]步優(yōu)化更新時的搜索方向。

在共軛梯度法中，有推論2：第[i]步計算的梯度[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è)置初始猜測[u0]，右端項為[f]，最大迭代次數(shù)為[imax]。

②[r0=f-Ku0]，[z0=M-1r0]，[p0=z0]。

③對[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è)置的精度要求）時，結(jié)束迭代。

本研究使用前文所述的預(yù)處理共軛梯度法和GPU并行法對大規(guī)模拓?fù)湓O(shè)計進(jìn)行優(yōu)化，減少有限元求解的迭代次數(shù)，提高計算效率。

3 基于GPU拓?fù)鋬?yōu)化的關(guān)鍵技術(shù)

在共軛梯度法計算過程中，大部分計算時間都花費在對矩陣向量積（SpMV）的求解上。在共軛梯度法中只有一次對稀疏矩陣向量積（SpMV）的運算，其他部分的運算（包括計算內(nèi)積、標(biāo)量向量積、向量加法）都可以快速求解，可以使用Matlab直接進(jìn)行串行計算。

稀疏矩陣向量積的運算，即求解共軛梯度法中優(yōu)化步長[αi]中的[Kpi]。在大規(guī)模系統(tǒng)中，整體剛度矩陣規(guī)模過大，會大幅降低算法的運算速度。本研究采用了一種基于有限元計算的逐單元（EBE）免組裝方案，包括邊界條件和各種載荷的處理。該策略消除了組裝過程在單元級求解線性方程組，從而避免生成總體矩陣及稀疏矩陣格式轉(zhuǎn)換的開銷。在免組裝解決方案中，共軛梯度法中的計算密集型矩陣向量積SpMV被較小的單元級密集MvP所取代，每個線程僅使用單元剛度矩陣[K（e）]或其約束后的變體來進(jìn)行并行計算，因此計算時間將大幅度縮短[7]。計算公式見式（9）。

[b=e?εK（e）p（e）]? ? ? （9）

式中：[b]表示[Kpi]的矩陣向量積運算;[ε]是單元的集合。

矩陣向量積（SpMV）即計算[b=Kp]的CUDA kernel函數(shù)的基本算法如下。

①將線程分配到每個單元上。

②設(shè)id=blockDim.x×blockIdx.x+threadIdx.x，用來固定每個單元的線程索引號。

③For與該單元相關(guān)聯(lián)的單元[e]（各單元的線程開始同時計算）;

讀取單元剛度矩陣[K（e）]和其節(jié)點的位置信息;

計算[be=K（e）p（e）];

累加[be]到[b]。

④End。

⑤同步線程。

此外，在本研究設(shè)計的基于GPU和預(yù)處理共軛梯度法的優(yōu)化拓?fù)浣Y(jié)構(gòu)中，CUDA進(jìn)行原子操作過程中，在并行更新結(jié)果向量時可能會發(fā)生大量內(nèi)存沖突。此外，GPU并行運算僅支持單精度的CUDA原子，不支持雙精度的CUDA原子，所以會影響計算結(jié)果的精確度。為了解決CUDA原子操作的缺點，本研究采用兩階段方法[8]，并對拓?fù)鋬?yōu)化過程中的預(yù)處理矩陣的形成、柔順度計算、靈敏度計算、靈敏度過濾均使用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]個單元的設(shè)計變量;[xmin]為最小設(shè)計變量（非零以避免奇異）;[V（x）]和[V0]分別為材料體積設(shè)計域體積;[vol]是規(guī)定的體積分?jǐn)?shù)。

本研究使用的是SIMP材料懲罰模型，其關(guān)于設(shè)計變量的表達(dá)式見式（11）。

[Ee（xe）=（xe）pE0e]? ? ?（11）

式中：[E0e]為初始材料彈性模量;[p]為懲罰因子;[Ee]表示第[e]個單元的設(shè)計變量所對應(yīng)的單元彈性模量。經(jīng)過有限元離散后的結(jié)構(gòu)，考慮到單元剛度矩陣與彈性模量具有簡單的線性關(guān)系，SIMP材料懲罰模型也可以表示為[ke=（xe）αkk0e]。

5 優(yōu)化算例

為了驗證本研究所提出的優(yōu)化方法，通過上述方法，對圖2所示的MBB梁結(jié)構(gòu)以柔順度最小為目標(biāo)進(jìn)行拓?fù)鋬?yōu)化求解，并使用式（10）的優(yōu)化模型。

MBB梁結(jié)構(gòu)左右兩側(cè)下端兩點支撐在滾輪上，本研究針對MBB梁結(jié)構(gòu)的一半進(jìn)行優(yōu)化設(shè)計。在頂部中點沿豎直向下方向施加集中豎向載荷1 N。MBB梁結(jié)構(gòu)初始設(shè)計域的尺寸為4 m×1 m×0.004 m，其初始材料的彈性模量為E=2×1011 Pa，其泊松比為v=0.3。將圖2所示設(shè)計域離散為800×200=160 000個四節(jié)點矩形平面應(yīng)力單元。設(shè)置密度過濾半徑為[rmin=2D]，其中[D]為最大單元邊長。在此示例中，目標(biāo)體積與設(shè)計域初始體積的比值設(shè)定為0.4。

本算例采用OC算法進(jìn)行求解，其中[Xmin]設(shè)置為0.001，其他參數(shù)根據(jù)本研究中相應(yīng)章節(jié)的建議值選取。

圖3為使用了本研究的預(yù)處理共軛梯度法，并且通過GPU并行對該算例進(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ī)模巨大，有十幾萬個單元的網(wǎng)格劃分，如果使用基于串行計算的大規(guī)模結(jié)構(gòu)拓?fù)鋬?yōu)化方法[10]，會導(dǎo)致其優(yōu)化速度過于緩慢。經(jīng)過試驗驗證，使用傳統(tǒng)的串行計算導(dǎo)致每一次優(yōu)化迭代的時間在30 min以上，從而導(dǎo)致該算例的結(jié)構(gòu)拓?fù)鋬?yōu)化難以進(jìn)行。而本研究使用共軛梯度法進(jìn)行預(yù)處理，并用GPU進(jìn)行計算，使每次優(yōu)化迭代步小于1 min，從而提高了幾十倍的計算速度，所以對大規(guī)模拓?fù)鋬?yōu)化使用此方法是有效的，在很大程度上減少計算量，縮短計算時間。

6 結(jié)語

本研究基于GPU并行技術(shù)，提出了一種使用預(yù)處理共軛梯度迭代法的連續(xù)體靜力拓?fù)鋬?yōu)化方法。并通過一系列具體算例驗證了所提方法的可行性、有效性和其速度優(yōu)勢。可以得出以下結(jié)論：①以Matlab為基礎(chǔ)，通過調(diào)用CUDA C編寫的kernel核函數(shù)來實現(xiàn)拓?fù)鋬?yōu)化的并行計算;②為了便于GPU的并行計算，本研究使用迭代算法，且為了減少迭代法的結(jié)構(gòu)有限元分析迭代次數(shù)，引入了雅可比（Jacobi）對角線預(yù)處理器，研究了基于共軛梯度法和預(yù)處理技術(shù)的結(jié)構(gòu)有限元并行計算方法;③本研究開發(fā)的GPU并行計算方法來處理基于迭代法的有限元計算中的約束條件、矩陣向量積、預(yù)處理矩陣和拓?fù)鋬?yōu)化過程中的柔順度、靈敏度和靈敏度過濾等;④完成上述方法的Matlab和C++并行計算核函數(shù)編程，提出了基于GPU和預(yù)處理共軛梯度法的結(jié)構(gòu)拓?fù)鋬?yōu)化方法，并完成了算例驗證，驗證了本研究提出方法的有效性和計算效率。

綜上所述，本研究所提出的方法有待于推廣到基于GPU并行計算的多工況和多約束的大規(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并行計算技術(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].河南科技，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].長沙：國防科學(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].長沙理工大學(xué)學(xué)報（自然科學(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.