CN109960865A - 薄板无网格Galerkin结构动力响应分析的GPU加速方法 - Google Patents

Abstract

本发明公开了一种薄板无网格Galerkin结构动力响应分析的GPU加速方法。它通过CPU将薄板数据读入到主机内存，根据输入数据计算GPU并行所需相关数据，并将得到的数据复制到GPU全局存储器中；在GPU上设置线程块与线程数量，并行加速计算节点相应的形函数值，然后建立GPU线程块与交叉节点对一一映射模式，并行加速组装薄板的总体刚度矩阵与质量矩阵；在GPU上根据边界条件修正总体刚度矩阵，然后由质量矩阵和修正后的刚度矩阵求得阻尼矩阵、有效刚度矩阵，并将有效刚度矩阵采用三角分解；在GPU中采用Newmark法对薄板进行动力响应分析，最后输出求解的位移、速度和加速度的结果。本发明极大地提高了动力响应分析的求解效率。

Description

薄板无网格Galerkin结构动力响应分析的GPU加速方法

技术领域

本发明属于计算机辅助工程中薄板的结构动力仿真分析技术领域，具体涉及一种薄板无网格Galerkin(伽辽金)结构动力响应分析的GPU(Graphic Processing Unit，图形处理器)加速方法。

背景技术

薄板动力响应分析是研究薄板在动载荷的激励下所产生的位移、速度和加速度等力学参量的动态响应信息，是薄板结构动力设计与灵敏度分析及优化的重要基础之一。在工程中，许多领域都需要采用板壳结构，如车辆、船舶、飞机和压力容器等，然而，这些设备经常处于各种激励的作用，组成部件中的薄板不可避免地产生各种各样的振动。当薄板振动量超出允许范围之后，将会影响设备的工作性能，缩短设备的使用寿命，甚至造成严重的事故。因此，现代工程设计对薄板振动问题的解决提出了更严格的要求。

随着有限元法在工程领域的应用，薄板动力响应可以通过有限元法来模拟分析计算。但有限元法是基于网格的数值方法，网格的质量和尺寸严重影响计算的精度和计算工作量，而且复杂结构的网格生成也极具挑战性。近年来，无网格法得到迅速发展，它不需要划分网格，克服了有限元法对网格的依懒性，在涉及网格畸变时显示出明显优势。目前，无网格Galerkin法是应用最广的无网格计算方法，具有计算精度高、前后处理简便、自适应性好等优点，在结构动力学问题的分析中得到广泛应用。然而薄板的刚度矩阵和质量矩阵的组装耗时长、阻尼矩阵和总体离散系统方程求解速度慢等缺点，严重影响无网格Galerkin法在薄板结构动力分析中的应用。

近年来，随着计算技术和计算机技术的迅猛发展，CPU(Central ProcessingUnit)具有体积小、重量轻结构灵活、价格低廉等优点，但随着晶体管集成工艺进入瓶颈，晶体管密度增长达到物理极限。尽管Intel公司提出了多核CPU，但对CPU提高计算能力帮助有限，与此同时CPU芯片的时钟速率也已经达到了极限，很难再提高。然而，GPU在数据并行处理能力和存储器带宽上逐渐优于CPU，且GPU具有超高的浮点计算能力、性价比高、功耗小等优点，在通用计算领域取得了很大的进展。再者，GPU统一计算设备架构(CUDA)的常用编程模型是直接使用C语言进行开发，大大减少了编程工作。利用CUDA C调用GPU对薄板的刚度矩阵和质量矩阵并行组装，及薄板的阻尼矩阵和总体离散系统方程的求解，极大地提高了结构动力响应分析的计算效率，因此研究基于GPU并行加速的薄板无网格Galerkin法结构动力响应分析具有重要意义。

发明内容

本发明的目的在于针对薄板无网格Galerkin法结构动力响应分析中存在刚度矩阵与质量矩阵组装耗时长、阻尼矩阵和总体离散系统方程求解效率低等问题，提供一种薄板无网格Galerkin结构动力响应分析的GPU加速方法。

本发明的薄板无网格Galerkin结构动力响应分析的GPU加速方法，包括如下顺序的步骤：

(1)通过CPU读取薄板的节点坐标、背景积分网格、材料物理性能、约束条件和动态载荷数据，以及求解设置参数与动力响应的计算时间步长，并存储至主机内存中；在薄板内部和位移边界布置积分点，并求解薄板内部和位移边界上的积分点数据；通过对薄板的节点、积分点循环搜索，确定每个节点影响域内的积分点数据，以及每个积分点定义域内的节点数据；再对薄板的节点循环搜索，确定交叉节点对信息；然后，将薄板内部和位移边界上的积分点数据、交叉节点对信息存储至主机内存中；

(2)在CPU上将薄板的节点坐标、积分点数据、积分点定义域内的节点数据、交叉节点对信息由主机内存复制到GPU全局存储器中，设置GPU中的线程块与线程的数量，并行加速计算节点的形函数值、形函数一阶导数值及形函数二阶导数值；并行加速组装总体刚度矩阵K、总体惩罚刚度矩阵K^α和总体质量矩阵M，且以按行压缩存储CSR格式存储于GPU全局存储器中；

(3)在CPU中根据设置的时间步长与总时间，读取薄板的节点坐标、作用力以及力作用的时间，依次求得Newmark法每一时间步长下的总体外力向量，并对应存储于总体外力向量矩阵F中，最后由主机内存复制至GPU全局存储器中；

(4)读取GPU全局存储器中总体刚度矩阵K和总体惩罚刚度矩阵K^α数据，根据边界条件对总体刚度矩阵进行修正其中，为修正后的总体刚度矩阵；

(5)读取GPU全局存储器中总体质量矩阵M、修正后的总体刚度矩阵设置质量比例系数α、刚度比例系数β，在GPU上求解阻尼矩阵C，即

(6)设置位移向量矩阵a、速度向量矩阵加速度向量矩阵以及参数α₀、δ，并根据参数α₀、δ计算积分常数c₀、c₁、c₂、c₃、c₄、c₅、c₆、c₇；其中，α₀＝1/4、σ＝1/2、 c₆＝Δt(1-δ)、c₇＝δΔt，Δt为时间步长；

(7)读取GPU全局存储器中修正后的总体刚度矩阵总体质量矩阵M、阻尼矩阵C，在GPU上形成有效刚度矩阵即：并对有效刚度矩阵进行三角分解，即：其中，为有效刚度矩阵，L为三角分解后的下三角矩阵，L^T为三角分解后的上三角矩阵；

(8)求解每一时间步长(t＝0,Δt,2Δt,…)的位移、速度、加速度；其中，t为计算时间，Δt为时间步长；

(9)将主机内存中存储的所求节点每一时间步长的位移、速度、加速度的计算结果写入到输出文件中。

具体的，所述步骤(8)包括如下步骤：

(a)在CPU中读取位移向量矩阵a、速度向量矩阵加速度向量矩阵计算质量矩阵系数向量M_x与阻尼矩阵系数向量C_x： 并复制到GPU的全局存储器中；其中，a_t为t时刻位移向量，为t时刻速度向量，为t时刻加速度向量；

(b)读取GPU全局存储器中总体力向量矩阵F、总体质量矩阵M、阻尼矩阵C，计算t+Δt时刻的有效载荷向量即：其中，F_t+Δt为t+Δt时刻的总体力向量；

(c)在GPU上采用回代的方法求解线性方程组，计算t+Δt时刻的位移向量a_t+Δt，即其中，a_t+Δt为t+Δt时刻的位移向量，时刻的有效载荷向量，L为三角分解后的下三角矩阵，L^T为三角分解后的上三角矩阵；

(d)在GPU中计算t+Δt时刻的加速度向量速度向量即

(e)读取t+Δt时刻的位移向量a_t+Δt、速度向量加速度向量并由GPU的全局存储中复制到主机内存中位移向量矩阵a、速度向量矩阵加速度向量矩阵所对应的位置；

(f)判断t+Δt时刻是否达到所设置的总时间，若达到结束步骤(8)，若没达到，转到步骤(a)。

本发明与现有技术相比，具有如下有益效果：

(1)本发明在GPU上同时对总体质量矩阵与总体刚度矩阵进行组装，且阻尼矩阵与有效刚度矩阵的计算也在GPU上实现，因此极大地缩短了薄板动力响应分析过程中矩阵的组装计算时间；同时，本发明采用稀疏格式存储总体质量矩阵、总体刚度矩阵、阻尼矩阵与有效刚度矩阵，减小了存储需求，能够实现较大规模问题的计算。

(2)本发明采用Newmark法对薄板动力响应进行分析。在GPU上对薄板有效刚度矩阵仅进行一次三角分解，位移、速度与加速度响应值均在GPU上进行循环迭代计算，极大地减少了薄板动力响应分析的求解时间。

(3)本发明采用无网格Galerkin法，不需要网格，对薄板动力响应分析结果具有较高精度，而且本发明方法可广泛适用于二维任意几何形状的薄板动力响应分析问题。

附图说明

图1为本发明方法的主流程框图。

图2为本发明实施例的薄板示意图

图3为图2所示实施例的薄板节点示意图。

图4为本发明方法的薄板动力响应分析曲线图。

具体实施方式

下面结合附图和实施例对本发明作进一步详细的描述。

参见图1，本发明薄板无网格Galerkin结构动力响应分析的GPU加速方法，包括如下顺序的步骤：

(1)通过CPU读取薄板的节点坐标、背景积分网格、材料物理性能、约束条件和动态载荷数据，以及求解设置参数与动力响应的计算时间步长，并存储至主机内存中；在薄板内部和位移边界布置积分点，并求解薄板内部和位移边界上的积分点数据；通过对薄板的节点、积分点循环搜索，确定每个节点影响域内的积分点数据，以及每个积分点定义域内的节点数据；再对薄板的节点循环搜索，确定交叉节点对信息；然后，将薄板内部和位移边界上的积分点数据、交叉节点对信息存储至主机内存中。

(2)在CPU上将薄板的节点坐标、积分点数据、积分点定义域内的节点数据、交叉节点对信息由主机内存复制到GPU全局存储器中，设置GPU中的线程块与线程的数量，并行加速计算节点的形函数值、形函数一阶导数值及形函数二阶导数值；并行加速组装总体刚度矩阵K、总体惩罚刚度矩阵K^α和总体质量矩阵M，且以按行压缩存储CSR格式存储于GPU全局存储器中。

(3)在CPU中根据设置的时间步长与总时间，读取薄板的节点坐标、作用力以及力作用的时间，依次求得Newmark法每一时间步长下的总体外力向量，并对应存储于总体外力向量矩阵F中，最后由主机内存复制至GPU全局存储器中。

(4)读取GPU全局存储器中总体刚度矩阵K和总体惩罚刚度矩阵K^α数据，根据边界条件对总体刚度矩阵进行修正其中，为修正后的总体刚度矩阵。

(5)读取GPU全局存储器中总体质量矩阵M、修正后的总体刚度矩阵设置质量比例系数α、刚度比例系数β，在GPU上求解阻尼矩阵C，即

(6)设置位移向量矩阵a、速度向量矩阵加速度向量矩阵以及参数α₀、δ，并根据参数α₀、δ计算积分常数c₀、c₁、c₂、c₃、c₄、c₅、c₆、c₇；其中，α₀＝1/4、σ＝1/2、 c₆＝Δt(1-δ)、c₇＝δΔt，t为时间步长。

(7)读取GPU全局存储器中修正后的总体刚度矩阵总体质量矩阵M、阻尼矩阵C，在GPU上形成有效刚度矩阵即：并对有效刚度矩阵进行三角分解，即：其中，为有效刚度矩阵，L为三角分解后的下三角矩阵，L^T为三角分解后的上三角矩阵。

(8)求解每一时间步长(t＝0,Δt,2Δt,…)的位移、速度、加速度；其中，t为计算时间，Δt为时间步长；

具体包括如下过程：

(a)在CPU中读取位移向量矩阵a、速度向量矩阵加速度向量矩阵计算质量矩阵系数向量M_x与阻尼矩阵系数向量C_x： 并复制到GPU的全局存储器中；其中，a_t为t时刻位移向量，为t时刻速度向量，为t时刻加速度向量；

(b)读取GPU全局存储器中总体力向量矩阵F、总体质量矩阵M、阻尼矩阵C，计算t+Δt时刻的有效载荷向量即：其中，F_t+Δt为t+Δt时刻的总体力向量；

(c)在GPU上采用回代的方法求解线性方程组，计算t+Δt时刻的位移向量a_t+Δt，即其中，a_t+Δt为t+Δt时刻的位移向量，为t+Δt时刻的有效载荷向量，L为三角分解后的下三角矩阵，L^T为三角分解后的上三角矩阵；

(d)在GPU中计算t+Δt时刻的加速度向量速度向量即

(e)读取t+Δt时刻的位移向量a_t+Δt、速度向量加速度向量并由GPU的全局存储中复制到主机内存中位移向量矩阵a、速度向量矩阵加速度向量矩阵所对应的位置；

(f)判断t+Δt时刻是否达到所设置的总时间，若达到结束步骤(8)，若没达到，转到步骤(a)。

(9)将主机内存中存储的所求节点每一时间步长的位移、速度、加速度的计算结果写入到输出文件中。

下面是将本发明方法应用于工程上的一个实例：

参见图2，是一块边长为2m、厚度为0.1m的四边固定约束的方形薄板，材料的弹性模量为2.0×10⁸Pa，泊松比μ＝0.3，密度为ρ＝7800kg/m³，中间承受10kN恒力，力的作用时间是0.5s，时间步长0.001s。无网格Galerkin法动力响应分析的薄板包括9801个节点，其分布如图3所示。本发明方法针对本实施例的具体实施步骤如下：

第一步：将已经准备好的薄板数据(包括薄板的节点坐标、背景积分网格顶点的编号和坐标、边界节点对中两节点的编号、边界信息、惩罚系数)、材料物理参数(包括材料弹性模量、泊松比、材料密度)、节点影响域系数、作用力以及力的作用时间，输入到input_Dynamic.txt文件中，在程序运行时，系统自动将以上数据读入到主机内存中。

一个标准的input_Dynamic.txt文件格式如下表1：

表1

第二步：在CPU中开辟9604(本实施例中背景积分网格个数)*4(每个四边形背景积分网格内布置的积分点数)*4*sizeof(double)字节的存储空间用于存储薄板内部积分点信息的数组，然后根据本实施例四边形背景积分网格顶点的坐标与编号、位移边界节点，在薄板内部布置积分点，并求得积分点的坐标、权系数及雅可比系数。

第三步：在CPU中开辟396(本实施例中位移边界节点对个数)*4(每个单元节点对上布置的积分点数)*4*sizeof(double)字节的存储空间用于存储位移边界积分点数组，然后通过CPU读取主机内存中存储的薄板背景积分网格的顶点坐标与编号，边界节点对中两节点的编号，在边界上布置积分点，求得边界上积分点的坐标、权系数、雅克比系数。

第四步：在CPU中开辟9801(模型节点的个数)*2(x、y两个方向)*sizeof(double)字节的存储空间用于存储节点的动态影响域半径数组，然后对薄板的节点循环，确定其它节点和该节点的最小距离r_min，再利用节点影响域半径＝3.0(节点影响域系数)*r_min，求得该节点的影响域半径，同理求得积分点的定义域半径。

第五步：在CPU中读取薄板的节点坐标、积分点坐标、节点的影响域半径、积分点的定义域半径，根据节点影响域半径，确定组成交叉节点对信息、每个节点影响域内的积分点数据以及每个积分点定义域内的节点数据；并将薄板的节点坐标、积分点坐标、交叉节点对信息、每个节点影响域内的积分点数据以及每个积分点定义域内的节点数据由主机内存复制到GPU全局存储器中，并为它们分配和主机内存中同样大小的存储空间。

第六步：读取GPU全局存储器中存储的薄板的节点坐标、积分点数据、积分点定义域内的节点数据，建立线程与积分点一一映射模式，通过最小二乘逼近，并行求得每个积分点定义域内节点的形函数值、形函数对应于各个坐标轴的二阶导数值，在本实施例中，为形函数对x、y轴的二阶导数值。

第七步：在GPU中建立线程块和交叉节点对一一映射模式，并行加速组装总体质量矩阵M、总体刚度矩阵K和总体惩罚刚度矩阵K^α，并以按行压缩存储CSR格式存储于GPU全局存储器中。

第八步：在CPU中根据设置的时间步长0.001s与总时间0.5s，以及存储在主机中薄板的节点坐标、作用力以及力作用的时间，依次求得Newmark法每一步长下的总体外力向量，并对应存储于总体外力向量矩阵F中，最后复制到GPU全局存储器中。

第九步：在GPU中读取GPU全局存储器中存储的总体刚度矩阵K和总体惩罚刚度矩阵K^α数据，根据边界条件对总体刚度矩阵进行修正

第十步：在GPU中根据总体质量矩阵M、修正后的总体刚度矩阵设置质量比例系数α、刚度比例系数β，在GPU上求解阻尼矩阵C，即

第十一步：在CPU中设置位移向量矩阵a、速度向量矩阵加速度向量矩阵以及参数α₀、δ，并根据参数α₀、δ计算积分常数c₀、c₁、c₂、c₃、c₄、c₅、c₆、c₇；其中，α₀＝1/4、δ＝1/2、 c₆＝Δt(1-δ)、c₇＝δΔt，Δt为时间步长。

第十二步：在GPU中根据修正后的总体刚度矩阵总体质量矩阵M、阻尼矩阵C，在GPU上形成有效的刚度矩阵即：

第十三步：在GPU中对有效刚度矩阵进行三角分解，即：

第十四步：在CPU中读取位移向量矩阵a、速度向量矩阵加速度向量矩阵计算质量矩阵与阻尼矩阵的系数向量并赋值： 并复制到GPU的全局存储器中。

第十五步：在GPU中读取总体力向量矩阵F、质量矩阵M、阻尼矩阵C，计算t+Δt时刻的有效载荷即

第十六步：在GPU中读取有效刚度矩阵三角分解数据L，采用回代的方法计算t+Δt时刻的位移a_t+Δt，即

第十七步：在GPU中计算t+Δt时刻的加速度速度即

第十八步：读取t+Δt时刻的位移向量a_t+Δ_t、速度向量加速度向量并由GPU的全局存储中复制到主机内存中位移向量矩阵a、速度向量矩阵加速度向量矩阵所对应的位置。

第十九步：判断t+Δt时刻是否达到所设置的总时间，若达到转到第二十步，若没达到，转到第十四步。

第二十步：将主机内存中存储的所求节点的每一时间步长的位移、速度、加速度的计算结果写入到输出文件中，本实施例薄板无网格Galerkin法结构动力响应分析的结果如图4表示(注：参考解为有限元在网格极细密情况下所求的解)。

采用上述方法进行薄板动力响应分析，整理CPU与GPU动力响应的计算时间(s)如表2所示：

表2

节点数	CPU	GPU	加速比
				9801	759.713	28.219	26.82

注：

虽然参考优先实施例对本发明进行描述，但以上所述并不限定本发明的保护范围，凡属于本发明思路下的技术方案均属于本发明的保护范围，任何在本发明的精神及原则内的修改、改进等，也应视为本发明的保护范围。

Claims (2)

1.一种薄板无网格Galerkin结构动力响应分析的GPU加速方法，其特征在于包括如下顺序的步骤：

(1)通过CPU读取薄板的节点坐标、背景积分网格、材料物理性能、约束条件和动态载荷数据，以及求解设置参数与动力响应的计算时间步长，并存储至主机内存中；在薄板内部和位移边界布置积分点，并求解薄板内部和位移边界上的积分点数据；通过对薄板的节点、积分点循环搜索，确定每个节点影响域内的积分点数据，以及每个积分点定义域内的节点数据；再对薄板的节点循环搜索，确定交叉节点对信息；然后，将薄板内部和位移边界上的积分点数据、交叉节点对信息存储至主机内存中；

(2)在CPU上将薄板的节点坐标、积分点数据、积分点定义域内的节点数据、交叉节点对信息由主机内存复制到GPU全局存储器中，设置GPU中的线程块与线程的数量，并行加速计算节点的形函数值、形函数一阶导数值及形函数二阶导数值；并行加速组装总体刚度矩阵K、总体惩罚刚度矩阵K^α和总体质量矩阵M，且以按行压缩存储CSR格式存储于GPU全局存储器中；

(3)在CPU中根据设置的时间步长与总时间，读取薄板的节点坐标、作用力以及力作用的时间，依次求得Newmark法每一时间步长下的总体外力向量，并对应存储于总体外力向量矩阵F中，最后由主机内存复制至GPU全局存储器中；

(4)读取GPU全局存储器中总体刚度矩阵K和总体惩罚刚度矩阵K^α数据，根据边界条件对总体刚度矩阵进行修正其中，为修正后的总体刚度矩阵；

(5)读取GPU全局存储器中总体质量矩阵M、修正后的总体刚度矩阵设置质量比例系数α、刚度比例系数β，在GPU上求解阻尼矩阵C，即

(6)设置位移向量矩阵a、速度向量矩阵加速度向量矩阵以及参数α₀、δ，并根据参数α₀、δ计算积分常数c₀、c₁、c₂、c₃、c₄、c₅、c₆、c₇；其中，α₀＝1/4、σ＝1/2、 c₆＝Δt(1-δ)、c₇＝δΔt，Δt为时间步长；

(7)读取GPU全局存储器中修正后的总体刚度矩阵总体质量矩阵M、阻尼矩阵C，在GPU上形成有效刚度矩阵即：并对有效刚度矩阵进行三角分解，即：其中，为有效刚度矩阵，L为三角分解后的下三角矩阵，L^T为三角分解后的上三角矩阵；

(8)求解每一时间步长(t＝0,Δt,2Δt,…)的位移、速度、加速度；其中，t为计算时间，Δt为时间步长；

(9)将主机内存中存储的所求节点每一时间步长的位移、速度、加速度的计算结果写入到输出文件中。

2.根据权利要求1所述薄板无网格Galerkin结构动力响应分析的GPU加速方法，其特征在于：所述步骤(8)包括如下步骤：

(a)在CPU中读取位移向量矩阵a、速度向量矩阵加速度向量矩阵计算质量矩阵系数向量M_x与阻尼矩阵系数向量C_x： 并复制到GPU的全局存储器中；其中，a_t为t时刻位移向量，为t时刻速度向量，为t时刻加速度向量；

(b)读取GPU全局存储器中总体力向量矩阵F、总体质量矩阵M、阻尼矩阵C，计算t+Δt时刻的有效载荷向量即：其中，F_t+Δt为t+Δt时刻的总体力向量；

(c)在GPU上采用回代的方法求解线性方程组，计算t+Δt时刻的位移向量a_t+Δt，即其中，a_t+Δt为t+Δt时刻的位移向量，为t+Δt时刻的有效载荷向量，L为三角分解后的下三角矩阵，L^T为三角分解后的上三角矩阵；

(d)在GPU中计算t+Δt时刻的加速度向量速度向量即

(e)读取t+Δt时刻的位移向量a_t+Δt、速度向量加速度向量并由GPU的全局存储中复制到主机内存中位移向量矩阵a、速度向量矩阵加速度向量矩阵所对应的位置；

(f)判断t+Δt时刻是否达到所设置的总时间，若达到结束步骤(8)，若没达到，转到步骤(a)。