基于GPU的结构静力拓扑优化设计方法

吴超

摘 要：针对连续体结构拓扑优化存在的计算量大、计算效率低等问题，开展了基于GPU并行计算的大规模结构静力拓扑优化方法研究。首先，为了减少有限元分析的迭代次数，引入了雅可比（Jacobi）对角线预处理器，研究基于共轭梯度法和预处理技术的结构有限元并行计算方法。其次，基于单元免组装技术，结合并行迭代计算方法，研究基于GPU的结构静力拓扑优化并行计算方法。在完成上述方法的Matlab和C++并行计算核函数编程后，进行了大量的算例考核。通过给出的算例来验证提出方法的有效性和计算效率，结果表明，该方法具有重要的理论价值和工程应用前景。

关键词：拓扑优化;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 引言

提高大規模工程结构优化设计的计算效率的常用方法有两种。一是改进求解算法的效率，二是采用并行计算替代传统串行计算。科学计算中的单核（串行）计算的执行效率低，其对大规模计算通常需要数天甚至数周才能执行成功。因此，在计算中更多的是使用多核（并行）操作。由于堆叠更高计算力导致CPU硬件成本增高，且计算效率依旧差强人意，甚至会消耗更多的能量。因此，基于CPU的结构计算与优化很难推广应用到大规模工程结构的设计中。自从NVIDIA公司成功研发出图形处理单元（Graphics Processing Units，GPU）的硬件和统一计算设备架构（Compute Unified Device Architecture，CUDA）的软件架构[1]，利用GPU在计算上的巨大潜力，越来越多的科学工作者在工程分析和数值计算中使用GPU技术[2]。由于GPU拥有数千个内核（多核）用来并行执行数百万个线程[3]，所以GPU作为一种高性能的计算平台，以合理的成本拥有巨大的计算能力和更大的内存带宽，从而引起更多的关注。综上所述，为了能够解决大规模拓扑优化计算效率问题，使用基于CUDA kernel函数的GPU并行计算方法，并结合预处理共轭梯度法的有限元迭代算法，提出一种基于GPU的并行结构拓扑优化设计方法。

1 CUDA kernel函数并行计算方法

CUDA编程模型[4]即各CUDA线程在GPU中各自独立且互不影响地运行，而CPU与之协同计算。这表明CUDA kernel内核在GPU上运行，其他不能并行操作部分则在CPU上运行，且CUDA编程模型中每个线程都有自己独有的局部内存，每个线程块都有对该块所有线程的可见共享内存，并且所有线程都可访问所有线程共有的全局内存。基本的GPU模型执行步骤如下。

①在设备（GPU）上分配内存。

②将数据从主机（CPU）传输到设备。

③内核启动并行执行。

④复制结果到主机。

在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核函数对象，并指定执行函数的线程块和线程参数。

④在Matlab中通过feval函数调用子程序的方式进行并行计算。

以有限元计算为例，由于每个单元均在自己的数据空间内独立地进行相同计算。因此，编写kernel核函数时可采用一个CUDA线程对一个单元进行计算的策略，实现循环展开，从而实现对所有单元并行、无序的计算，实现提高计算效率的目的。

2 预处理共轭梯度法

共轭梯度方法是一种用于求解大型线性方程组的迭代算法，它是一种介于最速下降法和牛顿法之间的方法。对于大规模拓扑优化中的大型有限元计算问题，由于直接法不适合并行计算，共轭梯度法以增加求解步骤的处理时间为代价，来减少矩阵求逆操作的内存需求。由线性弹性问题离散化产生的线性方程组见式（1）[5]。

[Ku=f]        （1）

式中：[K]是总体刚度矩阵;[u]是位移向量;[f]是节点力向量。

在计算机中进行计算时，由于舍入误差不存在，使得每次迭代后的残差[ri]不能保持正交，从而使收敛速度降低。另外，当采用惩罚函数处理边界条件时，得到的线性方程组的系数矩阵通常具有很大的条件数，使其非常病态，从而导致迭代收敛的速度非常缓慢。尤其是对大型、超大型方程组而言，通常需要进行[n]次迭代才能收敛，这是实际计算所不能接受的。目前，有很多学者采用预处理共轭梯度法（PCG）以确保能够以更快的速度收敛。通过引入预处理矩阵[M]，对系数矩阵进行适当变换，以降低其条件数，从而达到加速收敛的目的[6]。

预处理共轭梯度法（PCG）通过引入一个对称正定矩阵[M=WTW]，将原始线性系统式（1）替换为等效系统式（2）。

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

式中：[M-1]为预处理矩阵的逆。

雅可比对角线预处理矩阵是迄今为止实现最简单、计算成本最低的预处理矩阵。在并行环境中，与其他预处理矩阵相比，对角线预处理矩阵具有较强的竞争力。雅可比预处理矩阵[M]实际上是由刚度矩阵[K]的对角线元素组成的。由于雅可比预处理器是对角矩阵，所以其逆也是对角的。雅可比对角线预处理矩阵的表示见式（3）。

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

式中：[k11，k22，...，knn]为刚度矩阵[K]的对角线元素;[n]为刚度矩阵[K]的维度。

在共轭梯度法中，[ui]表示第[i]步优化时更新的位移向量，则负梯度表示见式（4）。

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

式中：[ri]表示第[i]步优化更新时的负梯度方向，梯度[ri]也可称作每一步的残差向量，误差定义为[ui]与最优解[u*]之间的差值[ei]，即[ei=u*-][ui]。

在預处理共轭梯度法中，引入一个新的向量[zi]，在每次计算出残差向量[ri]以后，增加了一个求解步骤，见式（5）。

[zi=M-1ri]       （5）

式中：[zi]表示第[i]步优化更新时经过预处理的残差向量。

在共轭梯度法中，需要解决搜索方向[pi]和优化步长[αi]这两个问题。假定每一步的优化方向用[pi]表示（也称为搜索方向）。假设初始搜索方向为初始负梯度方向，则经过预处理的搜索方向[p0]见式（6）。

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

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

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

根据此推论，在以后的每次迭代中，优化步长[αi]更新为式（7）。

[αi=rTizipTiKpi]       （7）

式中：[αi]表示第[i]步优化更新时的优化步长;[pi]表示第[i]步优化更新时的搜索方向。

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

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

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

式中：第[i]步优化更新用于更新搜索方向[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并行法对大规模拓扑设计进行优化，减少有限元求解的迭代次数，提高计算效率。

3 基于GPU拓扑优化的关键技术

在共轭梯度法计算过程中，大部分计算时间都花费在对矩阵向量积（SpMV）的求解上。在共轭梯度法中只有一次对稀疏矩阵向量积（SpMV）的运算，其他部分的运算（包括计算内积、标量向量积、向量加法）都可以快速求解，可以使用Matlab直接进行串行计算。

稀疏矩阵向量积的运算，即求解共轭梯度法中优化步长[αi]中的[Kpi]。在大规模系统中，整体刚度矩阵规模过大，会大幅降低算法的运算速度。本研究采用了一种基于有限元计算的逐单元（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）]和其节点的位置信息;

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

累加[be]到[b]。

④End。

⑤同步线程。

此外，在本研究设计的基于GPU和预处理共轭梯度法的优化拓扑结构中，CUDA进行原子操作过程中，在并行更新结果向量时可能会发生大量内存冲突。此外，GPU并行运算仅支持单精度的CUDA原子，不支持双精度的CUDA原子，所以会影响计算结果的精确度。为了解决CUDA原子操作的缺点，本研究采用两阶段方法[8]，并对拓扑优化过程中的预处理矩阵的形成、柔顺度计算、灵敏度计算、灵敏度过滤均使用CUDA kernel函数编写的并行程序进行处理。

4 优化模型

基于GPU并行和预处理共轭梯度法的拓扑优化是以结构柔顺度最小为目标的函数，采用SIMP材料惩罚模型和OC准则法，包含体积约束的结构拓扑优化模型，见式（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]是规定的体积分数。

本研究使用的是SIMP材料惩罚模型，其关于设计变量的表达式见式（11）。

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

式中：[E0e]为初始材料弹性模量;[p]为惩罚因子;[Ee]表示第[e]个单元的设计变量所对应的单元弹性模量。经过有限元离散后的结构，考虑到单元刚度矩阵与弹性模量具有简单的线性关系，SIMP材料惩罚模型也可以表示为[ke=（xe）αkk0e]。

5 优化算例

为了验证本研究所提出的优化方法，通过上述方法，对图2所示的MBB梁结构以柔顺度最小为目标进行拓扑优化求解，并使用式（10）的优化模型。

MBB梁结构左右两侧下端两点支撑在滚轮上，本研究针对MBB梁结构的一半进行优化设计。在顶部中点沿竖直向下方向施加集中竖向载荷1 N。MBB梁结构初始设计域的尺寸为4 m×1 m×0.004 m，其初始材料的弹性模量为E=2×1011 Pa，其泊松比为v=0.3。将图2所示设计域离散为800×200=160 000个四节点矩形平面应力单元。设置密度过滤半径为[rmin=2D]，其中[D]为最大单元边长。在此示例中，目标体积与设计域初始体积的比值设定为0.4。

本算例采用OC算法进行求解，其中[Xmin]设置为0.001，其他参数根据本研究中相应章节的建议值选取。

图3为使用了本研究的预处理共轭梯度法，并且通过GPU并行对该算例进行拓扑优化的拓扑构型，其优化拓扑构型的柔顺度为5.8×105 N·m。

图4为MBB梁在优化拓扑构型的柔顺度优化历程。其中，柔顺度在优化前期下降速度较快，在优化后期趋于稳定。

由于该算例的结构规模巨大，有十几万个单元的网格划分，如果使用基于串行计算的大规模结构拓扑优化方法[10]，会导致其优化速度过于缓慢。经过试验验证，使用传统的串行计算导致每一次优化迭代的时间在30 min以上，从而导致该算例的结构拓扑优化难以进行。而本研究使用共轭梯度法进行预处理，并用GPU进行计算，使每次优化迭代步小于1 min，从而提高了几十倍的计算速度，所以对大规模拓扑优化使用此方法是有效的，在很大程度上减少计算量，缩短计算时间。

6 结语

本研究基于GPU并行技术，提出了一种使用预处理共轭梯度迭代法的连续体静力拓扑优化方法。并通过一系列具体算例验证了所提方法的可行性、有效性和其速度优势。可以得出以下结论：①以Matlab为基础，通过调用CUDA C编写的kernel核函数来实现拓扑优化的并行计算;②为了便于GPU的并行计算，本研究使用迭代算法，且为了减少迭代法的结构有限元分析迭代次数，引入了雅可比（Jacobi）对角线预处理器，研究了基于共轭梯度法和预处理技术的结构有限元并行计算方法;③本研究开发的GPU并行计算方法来处理基于迭代法的有限元计算中的约束条件、矩阵向量积、预处理矩阵和拓扑优化过程中的柔顺度、灵敏度和灵敏度过滤等;④完成上述方法的Matlab和C++并行计算核函数编程，提出了基于GPU和预处理共轭梯度法的结构拓扑优化方法，并完成了算例验证，验证了本研究提出方法的有效性和计算效率。

综上所述，本研究所提出的方法有待于推廣到基于GPU并行计算的多工况和多约束的大规模三维结构拓扑优化问题。

参考文献：

[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方法拓扑优化设计[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] 荣见华.渐进结构优化方法及其应用研究[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] 荣见华，彭罗，易继军，等.一种新的多输入多输出柔顺机构拓扑优化方法[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.