基于GPU加速的等幾何拓撲優化高效多重網格求解方法

楊峰 羅世杰 楊江鴻 王英俊

摘要：

針對大規模等幾何拓撲優化（ITO）計算量巨大、傳統求解方法效率低的問題，提出了一種基于樣條h細化的高效多重網格方程求解方法。該方法利用h細化插值得到粗細網格之間的權重信息，然后構造多重網格方法的插值矩陣，獲得更準確的粗細網格映射信息，從而提高求解速度。此外，對多重網格求解過程進行分析，構建其高效GPU并行算法。數值算例表明，所提出的求解方法與線性插值的多重網格共軛梯度法、代數多重網格共軛梯度法和預處理共軛梯度法相比分別取得了最高1.47、11.12和17.02的加速比。GPU并行求解相對于CPU串行求解的加速比高達33.86，顯著提高了大規模線性方程組的求解效率。

關鍵詞：等幾何拓撲優化；方程組求解；h細化；多重網格法；GPU并行計算

中圖分類號：TH122

DOI：10.3969/j.issn.1004132X.2024.04.004

開放科學（資源服務）標識碼（OSID）：

A GPU-accelerated High-efficient Multi-grid Algorithm for ITO

YANG Feng1? LUO Shijie1? YANG Jianghong1? WANG Yingjun1，2

1.National Engineering Research Center of Novel Equipment for Polymer Processing，School of

Mechanical and Automotive Engineering，South China University of Technology，Guangzhou，510641

2.State Key Laboratory of Digital Manufacturing Equipment and Technology，Huazhong University

of Science and Technology，Wuhan，430074

Abstract： An efficient multi-grid equation solving method was proposed based on the h-refinement of splines to address the challenges posed by large-scale ITO computation and low efficiency of traditional solving methods. By the proposed method， the weight information obtained through h-refinement interpolation between coarse and fine grids was used to construct the interpolation matrix of the multi-grid method， thereby enhancing the accuracy of mapping information for both coarse and fine grids and improving computational efficiency. Additionally， a comprehensive analysis of the multi-grid solving process was conducted， culminating in the development of an efficient GPU parallel algorithm. Numerical examples illustrate that the proposed method outperforms existing methods， demonstrating speedup ratios of 1.47， 11.12， and 17.02 in comparison to the linear interpolation multi-grid conjugate gradient method algebraic multi-grid conjugate gradient method， and pre-processing conjugate gradient method respectively. Furthermore， the acceleration rate of GPU parallel solution surpasses that of CPU serial solution by 33.86 times， which significantly enhances the efficiency of solving large-scale linear equations.

Key words： isogeometric topology optimization（ITO）; system of equations; h-refinement; multi-grid method; GPU parallel computing

收稿日期：20230813

基金項目：國家自然科學基金（52075184）；數字制造裝備與技術國家重點實驗室開放基金（DMETKF2021020）

0? 引言

拓撲優化是一種根據給定載荷、約束條件及性能指標，對設計區域內材料按照優化準則自動分布，從而獲得高性能結構的方法，具有較高的設計自由度，是近些年結構設計、增材制造等領域的研究熱點［1］。目前常見的拓撲優化方法有變密度法［2］、均勻化方法［3］、水平集法［4］、漸進結構優化法［5］、移動變形組件法［6］等，這些方法已被廣泛應用于航空航天［7］、船舶汽車［8］、醫療器械［9］、橋梁建筑［10］等領域，具有重要的工程價值。

傳統拓撲優化的結構性能分析采用有限元法，其形函數單元間的C0連續在高階問題上難以滿足求解精度，而且在設計分析優化過程中，CAD數據與CAE數據表達方式不一致將導致模型轉換繁瑣，并導致優化結構模型的邊界曲線輪廓無法精確表達［11］。為了解決求解精度和連續性的要求，HUGHES等［12］采用CAD系統中的B樣條或非均勻有理B樣條（non-uniform rational basis spline， NURBS）基函數替代拉格朗日插值函數作為離散物理場的形函數，實現了CAD模型和CAE模型的統一表達，避免了模型數據轉化，并且提高了單元間連續性。

基于GPU加速的等幾何拓撲優化高效多重網格求解方法——楊? 峰? 羅世杰? 楊江鴻等

中國機械工程 第35卷 第4期 2024年4月

近年來，等幾何拓撲優化（isogeometric topology optimization，ITO）得到了廣泛的研究與發展。SEO等［13］基于等幾何方法提出了一種新型拓撲優化方法（即ITO），為結構設計、分析、優化的一體化研究奠定了基礎。KUMAR等［14］將等幾何分析（isogeometric analysis， IGA）和變密度拓撲優化集成，得到了基于控制點密度的ITO方法。GAO等［15］公開了一套高效、緊湊的ITO的MATLAB代碼，便于更多學者開展研究。關于ITO的介紹詳見WANG等［16］和GAO等［17］的綜述。目前，ITO已廣泛應用于各種工程實際問題的優化，如對梁結構的優化［18］、振動膜的形狀優化［19］、流體的優化［20］等。

ITO中每次迭代都需要對IGA方程進行求解［21］，該部分占據了ITO主要時間成本［22］，并且相比有限元分析，IGA剛度矩陣方程系數更加稠密，計算量更大。因此，研究IGA方程的高效求解方法是提高ITO計算效率的重要途徑。多重網格（multigrid， MG）法是求解偏微分方程數值解的有效方法之一，具有收斂性能好、求解速度快的優勢，通常用作黑箱求解器或者預處理共軛梯度法（preconditioning conjugate gradient， PCG）的預處理器［23］，被廣泛應用于大規模線性方程組的求解。劉石等［24］通過將IGA與多重網格共軛梯度法（multigrid conjugate gradient method， MGCG）結合，用于計算泊松方程，結果表明MGCG比MG法的收斂速度更快。WANG等［25］通過將多級網格、多重網格共軛梯度法和局部更新策略三者相結合，實現了ITO的三重加速，計算時間與傳統ITO相比減少了37%～93%，且優化結果保持一致。

在ITO中，即便使用了高效線性方程組求解方法MGCG，串行程序的計算成本仍很高，故必須尋求更高效率的解決方案。GPU并行計算是一種高效且經濟的加速方法，僅依賴普通顯卡即可實現數十倍甚至上百倍的加速比。因此，眾多學者采用GPU并行計算對線性方程組求解以及拓撲優化流程進行加速。鄧亮等［26］將GPU并行計算應用到流體力學偏微分方程的求解中，相對于串行算法實現了17倍左右的加速比。XIA等［27］將GPU并行和ITO結合，速度相比串行程序提升了兩個數量級。韓琪等［28］對雙向漸進結構拓撲優化算法進行改進，實現了大規模迭代算法的并行計算，驗證了并行計算方法的穩定性和高效性。

在ITO中，采用線性插值的多重網格共軛梯度法（linear interpolation MGCG， L_MGCG）的效果較好，但該方法存在數值不穩定的現象，該現象可能是線性插值關系不能精確表達ITO中樣條粗細網格的映射信息所導致的。基于上述研究，本文提出一種基于h細化插值網格的多重網格共軛梯度法（h-refinement interpolation MGCG， h_MGCG），在ITO中，該方法的插值關系相比L_MGCG更加準確，能更高效地完成方程求解。同時結合GPU并行計算，從硬件和算法兩個方面提高加速比，顯著提高了求解效率。最后，通過二維與三維的數值算例驗證了所提方法在ITO中具有適用性強、迭代收斂速度快、求解效率高以及穩定性能好的優點，并證明了基于GPU并行計算的求解方法與單CPU計算的求解方法具有相同的計算精度。

1? 理論基礎

1.1? NURBS基函數

IGA將NURBS基函數作為單元的形函數，從而實現CAD建模和CAE分析的統一表達。B樣條曲線是在Beziers曲線的基礎上提出的更為靈活的建模方式，它在精確表達模型的同時還能隨意增加控制點個數，從而進行局部修改。NURBS在B樣條的基礎上加入了權重系數，使其能夠準確表示二次規則曲線曲面，在ITO中獲得廣泛應用。NURBS包括節點矢量、控制點向量和權重等要素，其節點矢量可以是非均勻的，即相鄰節點間的距離可以不相等，可記為Ξ=（ξ1，ξ2，…，ξn+P+1），其中，ξi為第i個節點，P是基函數的階次，n是基函數和控制點的個數。B樣條的基函數由Cox-deBoor遞推公式定義為

P=0時：

Ni，0（ξ）=1? ξi≤ξ≤ξi+1

0其他（1）

P≥1時：

Ni，P（ξ）=ξ－ξiξi+P－ξiNi，P－1（ξ）+

ξi+P+1－ξξi+P+1－ξi+1Ni+1，P－1（ξ） （2）

引入權重系數wi，j，得到二維NURBS基函數表達式：

R（P，Q）i，j（ξ，η）=Ni，P（ξ）Ni，P（η）wi，j∑ni=1∑nj=1Ni，P（ξ）Nj，Q（η）wi，j（3）

式中，Ni，P、Nj，P分別為兩個參數方向對應的B樣條基函數；P、Q分別為沿曲面兩個參數方向的基函數階次。

1.2? 基于控制點密度的ITO

在基于密度法的ITO中，采用的材料插值方法為固體各向同性材料懲罰法（solid isotropic material with penalization，SIMP），每個單元將被賦予取值范圍在0～1之間的相對單元密度，而其單元e（e=1，2，…，n）的密度xe和彈性模量Ee之間的關系為

Ee（xe）=Emin+xpe（E0－Emin）? xe∈［0，1］（4）

式中，Emin為彈性模量最小值（Emin>0），其作用為防止剛度矩陣發生奇異；E0為相對單元密度為1的實心單元的彈性模量；p為懲罰系數，通常取3。

經典的以柔度為目標、以體積為約束的拓撲優化問題可表示為

find x=（x1，x2，…，xn）T

min C（x）=UTKU=∑Ne=1Ee（xe）uTekeue

s.t.KU=F

V（x）/V0≤VF

0≤xe≤1（5）

式中，C（x）為目標函數，表示結構柔度值；x為設計變量；N為設計域中設計變量的個數；U為總體位移向量；F為施加載荷向量；K為總體剛度矩陣；ue為單元e的位移向量；ke為單元e的單元剛度矩陣；V0為最初設計域總體積；VF為體積比約束（VF∈［0，1］）。

在IGA中，單元剛度矩陣ke計算公式如下：

ke=∫Ω^eSTDS|J1|d=∫eSTDS|J1||J2|d（6）

式中，S為應變位移矩陣；D為彈性矩陣；、分別為NURBS參數空間Ξ（m）和高斯積分參數域（m）對應的參數空間；J1、J2分別為從NURBS參數空間轉換到物理空間和從高斯積分域轉換到NURBS參數空間的雅可比矩陣。

與傳統的基于單元密度的SIMP法不同，基于控制點密度的ITO可將設計變量作為控制點網格中的高一維坐標。ITO中參數坐標為ξ的點的密度x可由附近控制點的密度求出：

x（ξ）=∑iNi（ξ）xi（7）

式中，Ni為第i個控制點的NURBS基函數；xi為控制點密度，本文以單元中心點的密度作為該單元的密度。

目標函數C對xe的偏導數由下式求得：

Cxe=－p（xe）p-1（E0－Emin）uTekeue（8）

因此，在ITO中計算靈敏度的公式由傳統SIMP法求解靈敏度公式改寫成

Cxi=∑j∈eiCxeijxeijxi=

∑j∈ei－p（xeij）p-1uTeijkeijueijxeijxi（9）

式中，ei為控制點i影響的所有單元；eij為ei中的第j個單元。

1.3? 基于h細化的多層網格策略

在ITO的實際應用中，精確表示物理模型需要使用精細網格。ITO通常采用h細化、p細化和k細化等方法進行細分，這些方法都可以保持原始形狀不變。本文采用h細化，即將多個節點插入NURBS曲線而不改變其階次，該過程只改變向量空間的基底，而曲線的幾何形狀和參數化保持不變。在h細化中，插入節點∈［ξk，ξk+1］到節點向量T=（ξ0，ξ1，…，ξn+k+1）可得到新的節點向量

=（0=ξ0，…，k=ξk，k+1=，k+2=ξk+1，…，n+k+2=ξn+k+1），然后得到新樣條控制點對應的坐標計算公式：

Qnewi=αiQoldi+（1－αi）Qoldi－1（10）

其中，Qnew為新控制點的齊次坐標；Qold為舊控制點的齊次坐標；αi由下式得到：

αi=1????? i≤k－p

－ξiξi+p－ξik－p+1≤i≤k

0i≥k+1（11）

因此，h細化在數學上可表示為矩陣運算：Qnew=P·Qold，P為將舊控制點序列映射到新控制點序列的投影系數矩陣。而細網格到粗網格的限制矩陣求解比較困難，本文直接取投影系數矩陣的轉置R=PT。二維的情況同理，只需在u、v兩個方向分別進行h細化得到映射矩陣Pu、Pv，從舊控制點到新控制點的投影系數矩陣為

Puv=Pu·Pv（12）

以圖1為例，細化前節點向量為（0，0，0，0.5，1，1，1）×（0，0，0，1，1，1），在u、v方向進行h細化，細化后向量為（0，0，0，0.25，0.5，0.75，1，1，1）×（0，0，0，0.5，1，1，1），細化前后各參數方向的NURBS基函數如圖1所示。

2? h_MGCG及其并行算法實現

2.1? 基于h細化的MGCG

多重網格方法主要包括幾何多重網格

（geometric multigrid， GMG）和代數多重網格（algebraic geometric multigrid， AMG）兩類，AMG不依賴實際網格信息與幾何模型，完全根據矩陣的代數信息求得各層級的系數矩陣，GMG依賴幾何信息得到一系列不同規模的矩陣。但是兩者求解的流程基本相同，都可分為前處理階段和求解階段。前處理階段生成多層級矩陣以便于后續的求解；求解階段由共軛梯度（conjugate gradient， CG）和多重網格（MG）兩部分組成，通過反復迭代來獲得滿足條件的最終結果。具體過程包括將網格進行粗化，在細網格上先進行迭代以消除高頻殘差分量，然后將難以消除的低頻分量映射到粗網格，使其變成相對高頻分量，并在粗網格上進行求解以消除相對高頻分量，最后將結果返回到細網格，循環反復前述過程直至達到設定的精度要求。

本文中，位移向量為U，載荷施加向量為F，剛度矩陣為K，對應不同層級的向量和矩陣分別表示為Ul、Fl、Kl，多重網格中的插值矩陣Pl即為前文所述h細化中新舊控制點的映射關系，限制矩陣用Rl=PTl代替，粗化矩陣Kl+1可以通過Galerkin投影方法［29］計算得到：

Kl+1=RlKlPl（13）

MG求解階段可以分為前光滑處理、計算殘差、限制投影、粗矩陣直接求解、插值校正、后光滑處理這6個步驟，MG算法具體求解流程見表1。

前后光滑處理的目的是減小求解中的振蕩誤差，通常采用Jacobi、Gauss-Seidel和SOR平滑處理。本文的光滑處理采用帶權重的Jacobi矩陣：

Ul=Ul+wD-1l（Fl-KlUl）（14）

式中，Dl為矩陣Kl的對角矩陣；w為權重系數，具體細節可參考文獻［30］。

預處理共振梯度（PCG）方法收斂速度快，但需要良好的預條件，而MG本身就是一種有效的預條件。將MG與PCG相結合，通過調用MG_Solve求解一次Kzi=ri方程取得的近似解代替PCG中的zi=M-1ri。MGCG算法流程見表2。

如前文所述，h_MGCG與其他MGCG方法的主要區別在于插值矩陣和限制矩陣的獲取方式不同，其他實現步驟相同。在h細化過程中，盡管提取插值矩陣P需要通過樣條插值完成，但是在整個拓撲優化過程中P保持不變，故只需在拓撲優化前處理階段執行一次提取操作。相比之下，在AMG中每次進行拓撲優化迭代得到新剛度矩陣K后都要重新提取代數網格信息。因此，h_MGCG求解方法能減少拓撲優化過程中MGCG前處理階段的總時間。

2.2? h_MGCG并行計算流程

由2.1所述理論基礎，h_MGCG求解過程中CG迭代和MG部分求解階段，主要進行的是易于并行的矩陣運算，故該部分可采用GPU進行并行計算。而在MG部分的前處理階段，主要工作是通過h細化得到插值矩陣P，該階段只需進行一次且花費時間很少，所以不采用GPU并行計算。本文為了提高粗矩陣部分的求解效率，采用cholesky分解進行直接求解，該求解算法為串行算法，且粗矩陣規模較小，GPU與CPU之間的通信開銷很小，故該部分為CPU串行計算。最后再將結果返回到GPU，進行后續插值修正和后光滑處理。h_MGCG并行求解方法（GPU_h_MGCG）以及整個ITO中的流程如圖2所示。

2.3? 并行計算關鍵核函數

結合MGCG算法流程及圖2所示的流程圖，在h_MGCG求解線性方程組KU=F的過程中，耗時較多的部分有向量的加減運算、標量與向量的乘法運算、矩陣與向量的乘法運算、向量的點積運算，在CPU串行求解方法中這些運算占用大部分時間。針對這些耗時較多的求解過程，利用CUDA架構將其在GPU上實現并行化。向量的加減運算、標量與向量的乘法運算根據向量元素個數n分配線程塊（block）數量num_blocks，將每個Block中線程（thread）數目num_threads取512，則num_blocks可通過下式求得：

num_blocks=（n+num_threads-1）/num_threads

其中，第id個線程將X、Y向量的第id個元素分別和常數a、b相乘并相加，將結果存在Z向量中，偽代碼如下：

if（id

Z［id］=aX［id］+bY［id］;

end

在MG中前后光滑處理的帶權重的Jacobi迭代法的迭代式即為式（14），

其中涉及矩陣求逆，由于該Dl矩陣為對角陣，故對其求逆只需取對角線元素的倒數即可，其線程分配同上，偽代碼如下：

if（id

d［id］［id］=1/D［id］［id］;

end

矩陣與向量的乘法運算實現難度相對較高，也是GPU_h_MGCG并行求解方法中最耗時的部分。為了節省空間和提高運算效率，該求解方法使用了壓縮行存儲（compressed sparse row， CSR）格式，IK、JK、val分別表示矩陣非零元素的行偏移量、列索引和值，例如矩陣K的存儲：

K=1700028050800604

IK=［0? 2? 4? 6? 8］

JK=［0? 1? 1? 2? 0? 2? 1? 3］

val=［1? 7? 2? 8? 5? 8? 6? 4］

矩陣K與向量X的乘法運算的一種簡單實現方式如下：在GPU中分配n個線程，其中第id個線程負責矩陣K中第id行與向量X相乘，并將結果儲存為向量Y的第id的元素。該方法的偽代碼如下：

if（id

for（i=IK［id］; i

Y［id］+=val［i］X［i］;

end

end

但是這種分配方式的缺點是同一線程束（warp）中相鄰線程不能訪問相鄰地址，因而內存訪問效率低下。為了提高內存訪問效率，同時減少線程資源閑置，本文的并行方案采用LightSpMV（light sparse matrix vector multiplication）方法［31］，該方法線程塊和線程調整更合理，從稀疏矩陣平均每行非零元素數量動態分配線程束子塊來求解，通過線程束shuffle指令實現同一線程束的線程交互，減小訪存的延遲（latency）。因為直接訪問寄存器比從共享內存（shared memory）中取值的延遲更低，因此LightSpMV方法實現了內存合并訪問，避免了線程資源的浪費，大幅提高并行求解效率，其線程塊及線程的分配、核函數的并行策略等具體細節可參考文獻［31-32］。

向量的點積運算的速率也是影響求解時間的關鍵因素，其步驟包括將兩個向量的對應元素相乘并將相乘的結果進行累加求和。由于最終結果依賴于每個線程求得的結果，故需要進行歸約求和。其核函數步驟比較復雜，包括申請共享內存空間、向量運算及同步操作。具體流程如下：①每個GPU線程負責兩對元素的相乘并求和，同時對線程塊內進行同步操作，保證向量運算結束時將所得的值存在共享內存之中；②將每個線程塊中的共享內存元素相加，每次分配num_threads/2n個線程，線程編號連續從0到num_threads/（2n-1），避免出現線程分支；③使用一個核函數，分配Block數量為1，負責塊內元素求和，得到最終結果。向量點積并行計算方式如圖3所示，該算法使線程資源得到充分利用，并且合并訪問與存儲，達到最大的帶寬。

結合上述關鍵運算，基于CUDA架構的h_MGCG并行程序整體框架如圖4所示。

3? 數值算例

本節通過1/4圓環和三維懸臂梁這兩種經典的數值算例來驗證ITO高效求解方法的準確性和高效性。1/4圓環用來驗證該求解方法的計算準確性、迭代收斂效果，并與常用求解方法Jacobi共軛梯度法（Jacobi_PCG）、不完全Cholesky分解共軛梯度法（ichol_PCG）、L_MGCG，以及AMG中常用三種不同粗化方式（基于R-S方法的beck_PCG［33］、基于非光滑聚合的HEM_PCG［34］、基于光滑聚合的CW_PCG［35］）的方法進行對比，比較這些求解方法在線性方程組求解時的加速效果。最后，通過三維懸臂梁進一步驗證該高效求解方法的通用性，以及在ITO迭代過程中的有效性和高效性。涉及的MGCG中多級矩陣最大層數均設為10，最粗矩陣規模為2000，即MGCG的網格層數達到10或者粗化的矩陣規模小于2000時，不再粗化網格，所有串行和并行程序都采用雙精度浮點型運算，測試的計算機CPU型號為Intel Xeon Gold 5218，GPU型號為NVIDIA GeForce RTX 3090。

3.1? 1/4圓環

圖5給出了1/4圓環的參數設置，其中圓環算例內圓半徑為1，外圓半徑為2，底邊施加固定約束，在左上角施加一個水平向右的集中載荷F=1。

將1/4圓環離散為64×64規模的二階NURBS單元，在拓撲優化中，體積分數設置為0.5，收斂準則設置為設計變量的最大變化量小于1%或者最大迭代次數為250，不同求解方法的殘差均設置為10-10，帶權重Jacobi的權重固定為w=0.4，在此算例中，MGCG的網格層數為2。圖6所示為用直接法和h_MGCG這兩種不同求解方法的拓撲優化最終迭代結果，可以看出h_MGCG求解方

法得到的最終優化結構和直接法的一致，其中本文的直接法采用Cholesky分解。

拓撲優化過程中兩種求解方法的柔度變化如圖7所示，拓撲優化過程中柔度變化一致，由表3可以看出h_MGCG求得的最終柔度與直接法的最終柔度一致，證明了h_MGCG求解算法的準確性。

將1/4圓環離散為256×256規模的二階NURBS單元，自由度數目為133 128，并以ITO中第12次迭代的線性方程求解為例，觀察求解過程中殘差變化。MGCG中光滑處理的Jacobi權重均選取w=0.4，網格層數為4。

圖8所示為第12次迭代中不同求解方法求解線性方程的殘差變化，可以看出，相比三種基于不同粗化方式的AMG方法，h_MGCG求解方法的收斂速度更快，

by different solving methods

在200次迭代內即可將求解殘差收斂到10-10，特別是與傳統的Jacobi_PCG方法相比，其迭代次數遠遠少于Jacobi_PCG方法的5000余次。此外，h_MGCG方法比L_MGCG方法的迭代次數少29，展現出更好的收斂效果。

各種求解方法花費的時間如圖9所示。結果表明，相較于其他求解方法，h_MGCG方法更加高效，比L_MGCG方法的速度提高了22.8%，與Jacobi_PCG和ichol_PCG方法相比分別具有16.48和2.60的加速比，與HEM_PCG、CW_PCG、beck_PCG方法相比，其加速比分別為6.20、6.27和6.43。

在求解拓撲優化第12次迭代時，不同多重網格求解方法的并行求解時間如圖10所示。結果表明，h_MGCG并行求解方法（GPU_h_MGCG）最為高效，僅需0.328 s即可完成一次求解計算，相比L_MGCG并行求解方法的0.406 s，其求解速度提高了23.8%，也遠超三種粗化方式AMG的并行求解速度。同時，GPU_h_MGCG方法相對于串行h_MGCG方法的加速比達到29.65，相對于傳統的串行Jacobi_PCG方法，加速比更是高達497.53，證明了h_MGCG方法的高效性。

3.2? 三維懸臂梁

為了驗證基于ITO的高效求解算法在三維問題上的適用性，本文采用三維懸臂梁算例，比較串行求解方法h_MGCG和并行求解方法GPU_h_MGCG的計算效率，并與其他傳統串行求解方法做對比。圖11所示為三維懸臂梁的設計域、邊界條件、外載荷，對其中一個大小為10×4的端面施加完全約束，在另一端施加密度為1、方向垂直向下的均布線載荷，ITO的體積分數設置為0.3，收斂準則設置為設計變量最大變化量小于1%或者最大收斂次數為250。在求解部分，所有MGCG、AMG中Jacobi權重都選擇收斂效果最好的固定權重w=0.3，MGCG層數為3，并將不同求解方法的求解殘差設置為10-8。

將三維懸臂梁離散為規模32×16×4的二階NURBS單元，自由度數目為11 016，不同求解方法在計算線性方程時得到的最終柔度值見表4，可以看出，h_MGCG和GPU_h_MGCG求解方法獲得的最終柔度值與直接法的最終柔度值幾乎一致。圖12所示為三維懸臂梁不同求解方法收斂曲線對比結果，可見h_MGCG方法和直接法求解依然具有幾乎相同的優化收斂曲線和結果，證明在三維結構上h_MGCG方法仍然準確。圖13所示為用直接法和h_MGCG這兩種不同求解方法的三維懸臂梁最終優化結果，可以看到三維算例中h_MGCG求解方法得到的最終優化結果和直接法的最終優化結果仍保持一致。

經過不同算例驗證，在ITO過程中，線性方程組的求解占據了大部分時間，占比甚至達90%以上。特別是就三維算例而言，由于剛度矩陣稠密度更高，故求解時間更長。本文對不同求解方法的效率進行比較。如圖14所示，在每次ITO的線性方程求解中，h_MGCG方法所需迭代次數明顯少于其他傳統求解方法，并優于L_MGCG方法，這再次驗證了由于h_MGCG方法具有更準確的插值關系而帶來的收斂性能的優勢。同時，由圖15可以看出，在每次迭代中線性方程求解所需時間上，h_MGCG方法也表現出比傳統求解方法更高效的特點，特別是相較于Jacobi_PCG方法。

各求解方法在網格規模為32×16×4時ITO的線性方程組總求解時間對比如圖16所示，相較于Jacobi_PCG方法，h_MGCG求解方法速度提高了17.02倍，與三種粗化方式AMG求解方法相比提高了9～11倍，相比L_MGCG求解方法則提高了47%。此外，在GPU并行計算的支持下，h_MGCG方法進一步實現加速，并將串行計算所需244.904 s的時間縮短至52.856 s，但其加速比僅為4.63，加速比不大的原因是規模較小。

為了研究網格規模與并行算法加速比之間的關系，圖17比較了64×32×8、128×32×8和128×64×16三種網格規模下的GPU_h_MGCG相對于h_MGCG的加速比，三種網格規模自由度數目分別為67 320、13 260和463 320，結果表明，相對于各串行算法，總求解時間的加速比分別為19.89、23.75和33.86。這是由于隨著規模的增大，GPU時間開銷占據總時間開銷的占比變大。因此在處理實際大規模問題時，GPU_h_MGCG求解方法能高效、準確地進行求解，其加速比是串行h_MGCG求解方法加速比的數十倍，是串行PCG求解方法加速比的數百倍甚至更高。

4? 結語

本文提出了一種基于h細化的高效多重網格方程求解方法，該方法利用h細化插值得到的粗細網格之間的權重信息構造多重網格方法的插值矩陣，獲得比傳統的線性插值方法更加準確的網格信息。將該方法與ITO相結合，只需在ITO細化過程中提取一次網格信息，可提高求解效率。同時，結合GPU并行計算，大幅提高ITO中線性方程組的求解速度。從硬件和算法兩方面提高求解速度，進而加快ITO求解效率。數值算例表明，h_MGCG求解方法最終的柔度值與直接法求解的結果幾乎一致，表明其求解結果能達到精度要求，同時其拓撲優化最終結果也保持一致。與其他求解方法相比，h_MGCG還具有更快的收斂速度和求解速度，與Jacobi_PCG、AMG和L_MGCG相比，取得了最高17.02、11.12和1.47的加速比。此外，結合GPU并行的GPU_h_MGCG求解方法充分利用了計算機資源，使得h_MGCG并行相對串行獲得了最高33.86的加速比，與傳統求解方法相比提升巨大。在大規模方程求解中，本文方法能夠快速高效、準確求解，對求解大規模ITO問題具有重要意義。

雖然本文方法實現了等幾何方程的高效求解，但仍有一些工作亟待完善。本文方法在進行平滑處理時，Jacobi迭代法權重的選取上采用了固定最佳權重，而對于不同算例，該最佳權重并不總是相同的。今后會針對不同算例自適應地選取Jacobi迭代法的權重，在增強通用性的同時減少迭代次數，以提高求解效率。此外，結合文獻［36］的工作，今后可以較容易地將本文方法擴展到復雜模型的等幾何分析之中，為實際工程應用提供更高效的求解方法。

參考文獻：

［1］? 楊雨豪，鄭偉，王英俊.一種自由度縮減和收斂加速的高效等幾何拓撲優化方法［J］.中國機械工程，2022，33（23）：2811-2821.

YANG Yuhao， ZHENG Wei， WANG Yingjun. An Efficient Isogeometric Topology Optimization Methods Using DOF Reduction and Convergence Acceleration［J］. China Mechanical Engineering， 2022， 33（23）：2811-2821.

［2］? BENDSE M P. Optimal Shape Design as a Material Distribution Problem［J］. Structural Optimization， 1989， 1（4）：193-202.

［3］? BENDSE M P， KIKUCHI N. Generating Optimal Topologies in Structural Design Using a Homogenization Method［J］. Computer Methods in Applied Mechanics and Engineering， 1988， 71（2）：197-224.

［4］? ALLAIRE G， JOUVE F， TOADER A-M. ALevel-set Method for Shape Optimization［J］. Comptes Rendus Mathematique， 2002， 334（12）：1125-1130.

［5］? QUERIN O M， STEVEN G P， XIE Y M. Evolutionary Structural Optimisation（ESO） Using a Bidirectional Algorithm［J］. Engineering Computations， 1998， 15（8）：1031-1048.

［6］? GUO X， ZHANG W， ZHANG J， et al. ExplicitStructural Topology Optimization Based on Moving Morphable Components（MMC） with Curved Skeletons［J］. Computer Methods in Applied Mechanics and Engineering， 2016， 310：711-748.

［7］? WANG Chuang， ZHU Jihong， WU Manqiao， et al. Multi-scale Design and Optimization for Solid-lattice Hybrid Structures and Their Application to Aerospace Vehicle Components［J］. Chinese Journal of Aeronautics， 2021， 34（5）：386-398.

［8］? GAO B， MENG D， SHI W， et al. Topology Optimization and the Evolution Trends of Two-speed Transmission of EVs［J］. Renewable and Sustainable Energy Reviews， 2022， 161：112390.

［9］? BRANDO P， SILVA P T， PARENTE M， et al. μSmartScope—Towards a Low-cost Microscopic Medical Device for Cervical Cancer Screening Using Additive Manufacturing and Optimization［J］. Proceedings of the Institution of Mechanical Engineers， Part L：Journal of Materials：Design and Applications， 2022， 236（2）：267-279.

［10］? LI Y， LAI Y， LU G， et al. Innovative Design of Long-span Steel-concrete Composite Bridge Using Multi-material Topology Optimization［J］. Engineering Structures， 2022， 269：114838.

［11］? 丁延冬，羅年猛，楊奧迪，等.Bézier單元剛度映射下的高效多重網格等幾何拓撲優化方法［J］.中國機械工程，2022，33（23）：2801-2810.

DING Yandong， LUO Nianmeng， YANG Aodi， et al. Efficient Multigrid Isogeometric Topology Optimization under Bézier Element Stiffness Mapping［J］. China Mechanical Engineering， 2022， 33（23）：2801-2810.

［12］? HUGHES T J R， COTTRELL J A， BAZILEVS Y. Isogeometric Analysis：CAD， Finite Elements， NURBS， Exact Geometry and Mesh Refinement［J］. Computer Methods in Applied Mechanics and Engineering， 2005， 194（39/41）：4135-4195.

［13］? SEO Y D， KIM H J， YOUN S K. IsogeometricTopology Optimization Using Trimmed Spline Surfaces［J］. Computer Methods in Applied Mechanics and Engineering， 2010， 199（49/52）：3270-3296.

［14］? KUMAR A V， PARTHASARATHY A. Topology Optimization Using B-spline Finite Elements［J］. Structural and Multidisciplinary Optimization， 2011， 44（4）：471-481.

［15］? GAO J， WANG L， LUO Z， et al. IgaTop：an Implementation of Topology Optimization for Structures Using IGA in MATLAB［J］. Structural and Multidisciplinary Optimization， 2021， 64（3）：1669-1700.

［16］? WANG Y， WANG Z， XIA Z， et al. Structural Design Optimization Using Isogeometric Analysis：a Comprehensive Review［J］. Computer Modeling in Engineering & Sciences， 2018， 117（3）：455-507.

［17］? GAO J， XIAO M， ZHANG Y， et al. A Comprehensive Review of Isogeometric Topology Optimization：Methods， Applications and Prospects［J］. Chinese Journal of Mechanical Engineering， 2020， 33（6）：34-47.

［18］? NAGY A P， ABDALLA M M， GURDAL Z. Isogeometric Sizing and Shape Optimization of Beam Structures［J］. Computer Methods in Applied Mechanics and Engineering， 2010，199：1216-1230.

［19］? NGUYEN D M，ANTON E，ALLAN R G，et al. Isogeometric Shape Optimization of Vibrating Membranes［J］. Computer Methods in Applied Mechanics and Engineering， 2011， 200：1343-1353.

［20］? YOUN D H， KIM M G， KIM H S. Shape Design Optimization of SPH Fluid-structure Interactions Considering Geometrically Exact Interfaces［J］. Structural and Multidisciplinary Optimization， 2011， 44：319-336.

［21］? 劉宏亮，祝雪峰，楊迪雄.基于等幾何分析的結構優化設計研究進展［J］.固體力學學報，2018，39（3）：248-267.

LIU Hongliang， ZHU Xuefeng， YANG Dixiong. Research Advances in Isogeometric Analysis-based Optimum Design of Structure［J］. Chinese Journal of Solid Mechanics， 2018， 39（3）：248-267.

［22］? AAGE N， ANDREASSEN E， LAZAROV B S， et al.Giga-voxel Computational Morphogenesis for Structural Design［J］. Nature， 2017， 550（7674）：84-86.

［23］? THOMAS S J， ANANTHAN S， YELLAPANTULA S， et al. AComparison of Classical and Aggregation-based Algebraic Multigrid Preconditioners for High-fidelity Simulation of Wind Turbine Incompressible Flows［J］. SIAM Journal on Scientific Computing， 2019， 41（5）：S196-S219.

［24］? 劉石，陳德祥，馮永新，等.等幾何分析的多重網格共軛梯度法［J］.應用數學和力學，2014，35（6）：630-639.

LIU Shi， CHEN Dexiang， FENG Yongxin， et al. A Multigrid Preconditioned Conjugate Gradient Method for Isogeometric Analysis［J］. Applied Mathematics and Mechanics， 2014， 35（6）：630-639.

［25］? WANG Y， LIAO Z， YE M， et al. An Efficient Isogeometric Topology Optimization Using Multilevel Mesh， MGCG and Local-update Strategy［J］. Advances in Engineering Software， 2020， 139：102733.

［26］? 鄧亮，徐傳福，劉巍，等.交替方向隱式CFD解法器的GPU并行計算及其優化［J］.計算機應用，2013，33（10）：2783-2786.

DENG Liang， XU Chuanfu， LIU Wei， et al. Parallelization and Optimization of Alternating Direction Implicit CFD Solver on GPU［J］. Journal of Computer Applications， 2013， 33（10）：2783-2786.

［27］? XIA Z， WANG Y， WANG Q， et al. GPU Parallel Strategy for Parameterized LSM-based Topology Optimization Using Isogeometric Analysis［J］. Structural and Multidisciplinary Optimization， 2017， 56：413-434.

［28］? 韓琪，蔡勇.基于GPU的大規模拓撲優化問題并行計算方法［J］.計算機仿真，2015，32（4）：221-226.

HAN Qi， CAI Yong. An Efficient Parallel Implementation of Large-scale Topology Optimization Problems Using GPU［J］. Computer Simulation， 2015， 32（4）：221-226.

［29］? AVNAT O， YAVNEH I. On the Recursive Structure of Multigrid Cycles［J］. SIAM Journal on Scientific Computing， 2022（V3）：S103-S126.

［30］? MADAMS A， SIFAKIS E， TERAN J. A Parallel Multigrid Poisson Solver for Fluids Simulation on Large Grids［C］∥Proceedings of the 2010 ACM SIGGRAPH/Eurographics Symposium on Computer Animation. Madrid，2010：65-74.

［31］? LIU Y， SCHMIDT B. Lightspmv：Faster Cuda-compatible Sparse Matrix-vector Multiplication Using Compressed Sparse Rows［J］. Journal of Signal Processing Systems， 2018， 90：69-86.

［32］? LIU Y， GONZLEZ-DOMNGUEZ J， SCHMIDT B. Faster Compressed Sparse Row（CSR）-based Sparse Matrix-vector Multiplication Using CUDA［C］∥GPU Technology Conference. Silicon Valley， 2015：17-20.

［33］? BECK R. Algebraic Multigrid by Component Splitting for Edge Elements on Simplicial Triangulations［R］. ZIB-Report（SC-99-40），1999.

［34］? NOTAY Y. Aggregation-based Algebraic Multigrid for Convection-diffusion Equations［J］. SIAM Journal on Scientific Computing， 2012， 34（4）：A2288-A2316.

［35］? BERNASCHI M， DAMBRA P， PASQUINI D. AMG Based on Compatible Weighted Matching for GPUs［J］. Parallel Computing， 2020， 92：102599.

［36］? WANG Y， GAO L， QU J， et al. Isogeometric Analysis Based on Geometric Reconstruction Models［J］. Frontiers of Mechanical Engineering， 2021，16：782-797.

（編輯? 陳? 勇）

作者簡介：

楊? 峰，男，1998年生，碩士研究生。研究方向為拓撲優化、GPU并行計算。

王英俊（通信作者），男，1984年生，教授、博士研究生導師。研究方向為拓撲優化、等幾何分析、CAD/CAE一體化。E-mail：wangyj84@scut.edu.cn。