楊峰 羅世杰 楊江鴻 王英俊

摘要：

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

關(guān)鍵詞：等幾何拓?fù)鋬?yōu)化；方程組求解；h細(xì)化；多重網(wǎng)格法；GPU并行計算

中圖分類號：TH122

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

開放科學(xué)（資源服務(wù)）標(biāo)識碼（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

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

0? 引言

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

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

基于GPU加速的等幾何拓?fù)鋬?yōu)化高效多重網(wǎng)格求解方法——楊? 峰? 羅世杰? 楊江鴻等

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

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

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

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

在ITO中，采用線性插值的多重網(wǎng)格共軛梯度法（linear interpolation MGCG， L_MGCG）的效果較好，但該方法存在數(shù)值不穩(wěn)定的現(xiàn)象，該現(xiàn)象可能是線性插值關(guān)系不能精確表達(dá)ITO中樣條粗細(xì)網(wǎng)格的映射信息所導(dǎo)致的?；谏鲜鲅芯?，本文提出一種基于h細(xì)化插值網(wǎng)格的多重網(wǎng)格共軛梯度法（h-refinement interpolation MGCG， h_MGCG），在ITO中，該方法的插值關(guān)系相比L_MGCG更加準(zhǔn)確，能更高效地完成方程求解。同時結(jié)合GPU并行計算，從硬件和算法兩個方面提高加速比，顯著提高了求解效率。最后，通過二維與三維的數(shù)值算例驗證了所提方法在ITO中具有適用性強、迭代收斂速度快、求解效率高以及穩(wěn)定性能好的優(yōu)點，并證明了基于GPU并行計算的求解方法與單CPU計算的求解方法具有相同的計算精度。

1? 理論基礎(chǔ)

1.1? NURBS基函數(shù)

IGA將NURBS基函數(shù)作為單元的形函數(shù)，從而實現(xiàn)CAD建模和CAE分析的統(tǒng)一表達(dá)。B樣條曲線是在Beziers曲線的基礎(chǔ)上提出的更為靈活的建模方式，它在精確表達(dá)模型的同時還能隨意增加控制點個數(shù)，從而進(jìn)行局部修改。NURBS在B樣條的基礎(chǔ)上加入了權(quán)重系數(shù)，使其能夠準(zhǔn)確表示二次規(guī)則曲線曲面，在ITO中獲得廣泛應(yīng)用。NURBS包括節(jié)點矢量、控制點向量和權(quán)重等要素，其節(jié)點矢量可以是非均勻的，即相鄰節(jié)點間的距離可以不相等，可記為Ξ=（ξ1，ξ2，…，ξn+P+1），其中，ξi為第i個節(jié)點，P是基函數(shù)的階次，n是基函數(shù)和控制點的個數(shù)。B樣條的基函數(shù)由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）

引入權(quán)重系數(shù)wi，j，得到二維NURBS基函數(shù)表達(dá)式：

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分別為兩個參數(shù)方向?qū)?yīng)的B樣條基函數(shù)；P、Q分別為沿曲面兩個參數(shù)方向的基函數(shù)階次。

1.2? 基于控制點密度的ITO

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

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

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

經(jīng)典的以柔度為目標(biāo)、以體積為約束的拓?fù)鋬?yōu)化問題可表示為

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）為目標(biāo)函數(shù)，表示結(jié)構(gòu)柔度值；x為設(shè)計變量；N為設(shè)計域中設(shè)計變量的個數(shù)；U為總體位移向量；F為施加載荷向量；K為總體剛度矩陣；ue為單元e的位移向量；ke為單元e的單元剛度矩陣；V0為最初設(shè)計域總體積；VF為體積比約束（VF∈［0，1］）。

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

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

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

與傳統(tǒng)的基于單元密度的SIMP法不同，基于控制點密度的ITO可將設(shè)計變量作為控制點網(wǎng)格中的高一維坐標(biāo)。ITO中參數(shù)坐標(biāo)為ξ的點的密度x可由附近控制點的密度求出：

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

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

目標(biāo)函數(shù)C對xe的偏導(dǎo)數(shù)由下式求得：

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

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

Cxi=∑j∈eiCxeijxeijxi=

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

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

1.3? 基于h細(xì)化的多層網(wǎng)格策略

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

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

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

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

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

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

0i≥k+1（11）

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

Puv=Pu·Pv（12）

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

2? h_MGCG及其并行算法實現(xiàn)

2.1? 基于h細(xì)化的MGCG

多重網(wǎng)格方法主要包括幾何多重網(wǎng)格

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

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

Kl+1=RlKlPl（13）

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

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

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

式中，Dl為矩陣Kl的對角矩陣；w為權(quán)重系數(shù)，具體細(xì)節(jié)可參考文獻(xiàn)［30］。

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

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

2.2? h_MGCG并行計算流程

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

2.3? 并行計算關(guān)鍵核函數(shù)

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

num_blocks=（n+num_threads-1）/num_threads

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

if（id

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

end

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

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

if（id

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

end

矩陣與向量的乘法運算實現(xiàn)難度相對較高，也是GPU_h_MGCG并行求解方法中最耗時的部分。為了節(jié)省空間和提高運算效率，該求解方法使用了壓縮行存儲（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的乘法運算的一種簡單實現(xiàn)方式如下：在GPU中分配n個線程，其中第id個線程負(fù)責(zé)矩陣K中第id行與向量X相乘，并將結(jié)果儲存為向量Y的第id的元素。該方法的偽代碼如下：

if（id

for（i=IK［id］; i

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

end

end

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

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

結(jié)合上述關(guān)鍵運算，基于CUDA架構(gòu)的h_MGCG并行程序整體框架如圖4所示。

3? 數(shù)值算例

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

3.1? 1/4圓環(huán)

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

將1/4圓環(huán)離散為64×64規(guī)模的二階NURBS單元，在拓?fù)鋬?yōu)化中，體積分?jǐn)?shù)設(shè)置為0.5，收斂準(zhǔn)則設(shè)置為設(shè)計變量的最大變化量小于1%或者最大迭代次數(shù)為250，不同求解方法的殘差均設(shè)置為10-10，帶權(quán)重Jacobi的權(quán)重固定為w=0.4，在此算例中，MGCG的網(wǎng)格層數(shù)為2。圖6所示為用直接法和h_MGCG這兩種不同求解方法的拓?fù)鋬?yōu)化最終迭代結(jié)果，可以看出h_MGCG求解方

法得到的最終優(yōu)化結(jié)構(gòu)和直接法的一致，其中本文的直接法采用Cholesky分解。

拓?fù)鋬?yōu)化過程中兩種求解方法的柔度變化如圖7所示，拓?fù)鋬?yōu)化過程中柔度變化一致，由表3可以看出h_MGCG求得的最終柔度與直接法的最終柔度一致，證明了h_MGCG求解算法的準(zhǔn)確性。

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

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

by different solving methods

在200次迭代內(nèi)即可將求解殘差收斂到10-10，特別是與傳統(tǒng)的Jacobi_PCG方法相比，其迭代次數(shù)遠(yuǎn)遠(yuǎn)少于Jacobi_PCG方法的5000余次。此外，h_MGCG方法比L_MGCG方法的迭代次數(shù)少29，展現(xiàn)出更好的收斂效果。

各種求解方法花費的時間如圖9所示。結(jié)果表明，相較于其他求解方法，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。

在求解拓?fù)鋬?yōu)化第12次迭代時，不同多重網(wǎng)格求解方法的并行求解時間如圖10所示。結(jié)果表明，h_MGCG并行求解方法（GPU_h_MGCG）最為高效，僅需0.328 s即可完成一次求解計算，相比L_MGCG并行求解方法的0.406 s，其求解速度提高了23.8%，也遠(yuǎn)超三種粗化方式AMG的并行求解速度。同時，GPU_h_MGCG方法相對于串行h_MGCG方法的加速比達(dá)到29.65，相對于傳統(tǒng)的串行Jacobi_PCG方法，加速比更是高達(dá)497.53，證明了h_MGCG方法的高效性。

3.2? 三維懸臂梁

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

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

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

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

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

4? 結(jié)語

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

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

參考文獻(xiàn)：

［1］? 楊雨豪，鄭偉，王英俊.一種自由度縮減和收斂加速的高效等幾何拓?fù)鋬?yōu)化方法［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單元剛度映射下的高效多重網(wǎng)格等幾何拓?fù)鋬?yōu)化方法［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］? 劉宏亮，祝雪峰，楊迪雄.基于等幾何分析的結(jié)構(gòu)優(yōu)化設(shè)計研究進(jìn)展［J］.固體力學(xué)學(xué)報，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］? 劉石，陳德祥，馮永新，等.等幾何分析的多重網(wǎng)格共軛梯度法［J］.應(yīng)用數(shù)學(xué)和力學(xué)，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并行計算及其優(yōu)化［J］.計算機應(yīng)用，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的大規(guī)模拓?fù)鋬?yōu)化問題并行計算方法［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年生，碩士研究生。研究方向為拓?fù)鋬?yōu)化、GPU并行計算。

王英?。ㄍㄐ抛髡撸?，男，1984年生，教授、博士研究生導(dǎo)師。研究方向為拓?fù)鋬?yōu)化、等幾何分析、CAD/CAE一體化。E-mail：wangyj84@scut.edu.cn。