吳超

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

關鍵詞：拓撲優(yōu)化;GPU并行;免組裝方法;共軛梯度法;預處理器

中圖分類號：TB34 ? ? 文獻標志碼：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ī)模工程結構優(yōu)化設計的計算效率的常用方法有兩種。一是改進求解算法的效率，二是采用并行計算替代傳統串行計算?？茖W計算中的單核（串行）計算的執(zhí)行效率低，其對大規(guī)模計算通常需要數天甚至數周才能執(zhí)行成功。因此，在計算中更多的是使用多核（并行）操作。由于堆疊更高計算力導致CPU硬件成本增高，且計算效率依舊差強人意，甚至會消耗更多的能量。因此，基于CPU的結構計算與優(yōu)化很難推廣應用到大規(guī)模工程結構的設計中。自從NVIDIA公司成功研發(fā)出圖形處理單元（Graphics Processing Units，GPU）的硬件和統一計算設備架構（Compute Unified Device Architecture，CUDA）的軟件架構[1]，利用GPU在計算上的巨大潛力，越來越多的科學工作者在工程分析和數值計算中使用GPU技術[2]。由于GPU擁有數千個內核（多核）用來并行執(zhí)行數百萬個線程[3]，所以GPU作為一種高性能的計算平臺，以合理的成本擁有巨大的計算能力和更大的內存帶寬，從而引起更多的關注。綜上所述，為了能夠解決大規(guī)模拓撲優(yōu)化計算效率問題，使用基于CUDA kernel函數的GPU并行計算方法，并結合預處理共軛梯度法的有限元迭代算法，提出一種基于GPU的并行結構拓撲優(yōu)化設計方法。

1 CUDA kernel函數并行計算方法

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

①在設備（GPU）上分配內存。

②將數據從主機（CPU）傳輸到設備。

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

④復制結果到主機。

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

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

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

②調用nvcc編譯器將其轉換成.ptx格式的文件。

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

④在Matlab中通過feval函數調用子程序的方式進行并行計算。

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

2 預處理共軛梯度法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

根據此推論，在以后的每次迭代中，優(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]共軛正交。

根據此推論，每次迭代更新的搜索方向[pi+1]更新為式（8）。

[pi+1=zi+1+βipi]? ? ?（8）

式中：第[i]步優(yōu)化更新用于更新搜索方向[pi+1]的內積[βi=rTi+1zi+1rTizi]。

預處理共軛梯度法的算法流程如下。

①設置剛度矩陣[K]，設置初始猜測[u0]，右端項為[f]，最大迭代次數為[imax]。

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

③對[i=1，2，...，imax]進行如下迭代。

[αi=rTizipTiKpi]

[ui+1=ui+αipi]

[ri+1=ri+αiKpi]

[zi+1=M-1ri+1]

[βi=rTi+1zi+1rTizi]

[pi+1=zi+1+βipi]

④當[ri2<ε]（[ε]為設置的精度要求）時，結束迭代。

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

3 基于GPU拓撲優(yōu)化的關鍵技術

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

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

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

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

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

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

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

③For與該單元相關聯的單元[e]（各單元的線程開始同時計算）;

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

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

累加[be]到[b]。

④End。

⑤同步線程。

此外，在本研究設計的基于GPU和預處理共軛梯度法的優(yōu)化拓撲結構中，CUDA進行原子操作過程中，在并行更新結果向量時可能會發(fā)生大量內存沖突。此外，GPU并行運算僅支持單精度的CUDA原子，不支持雙精度的CUDA原子，所以會影響計算結果的精確度。為了解決CUDA原子操作的缺點，本研究采用兩階段方法[8]，并對拓撲優(yōu)化過程中的預處理矩陣的形成、柔順度計算、靈敏度計算、靈敏度過濾均使用CUDA kernel函數編寫的并行程序進行處理。

4 優(yōu)化模型

基于GPU并行和預處理共軛梯度法的拓撲優(yōu)化是以結構柔順度最小為目標的函數，采用SIMP材料懲罰模型和OC準則法，包含體積約束的結構拓撲優(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]個單元的設計變量;[xmin]為最小設計變量（非零以避免奇異）;[V（x）]和[V0]分別為材料體積設計域體積;[vol]是規(guī)定的體積分數。

本研究使用的是SIMP材料懲罰模型，其關于設計變量的表達式見式（11）。

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

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

5 優(yōu)化算例

為了驗證本研究所提出的優(yōu)化方法，通過上述方法，對圖2所示的MBB梁結構以柔順度最小為目標進行拓撲優(yōu)化求解，并使用式（10）的優(yōu)化模型。

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

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

圖3為使用了本研究的預處理共軛梯度法，并且通過GPU并行對該算例進行拓撲優(yōu)化的拓撲構型，其優(yōu)化拓撲構型的柔順度為5.8×105 N·m。

圖4為MBB梁在優(yōu)化拓撲構型的柔順度優(yōu)化歷程。其中，柔順度在優(yōu)化前期下降速度較快，在優(yōu)化后期趨于穩(wěn)定。

由于該算例的結構規(guī)模巨大，有十幾萬個單元的網格劃分，如果使用基于串行計算的大規(guī)模結構拓撲優(yōu)化方法[10]，會導致其優(yōu)化速度過于緩慢。經過試驗驗證，使用傳統的串行計算導致每一次優(yōu)化迭代的時間在30 min以上，從而導致該算例的結構拓撲優(yōu)化難以進行。而本研究使用共軛梯度法進行預處理，并用GPU進行計算，使每次優(yōu)化迭代步小于1 min，從而提高了幾十倍的計算速度，所以對大規(guī)模拓撲優(yōu)化使用此方法是有效的，在很大程度上減少計算量，縮短計算時間。

6 結語

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

綜上所述，本研究所提出的方法有待于推廣到基于GPU并行計算的多工況和多約束的大規(guī)模三維結構拓撲優(yōu)化問題。

參考文獻：

[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].信息技術，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過濾和可行域調整的SIMP方法拓撲優(yōu)化設計[J].河南科技，2018（34）：26-28.

[6] BOLZ J，FARMER 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] 榮見華.漸進結構優(yōu)化方法及其應用研究[D].長沙：國防科學技術大學，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] 榮見華，彭羅，易繼軍，等.一種新的多輸入多輸出柔順機構拓撲優(yōu)化方法[J].長沙理工大學學報（自然科學版），2021（1）：66-78.

[10] SIGMUND O.A 99 line topology optimization code written in Matlab[J].Structural & Multidisciplinary Optimization，2001（2）：120-127.