CN109753682B - 一种基于gpu端的有限元刚度矩阵模拟方法 - Google Patents

Abstract

本发明提供一种基于GPU端的有限元刚度矩阵模拟方法，包括：建立刚度方程Ax＝b并初始化；分裂处理刚度矩阵A；对应分配GPU显存；将CPU端初始化数据复制到GPU上；调用计时函数开始计时；调用设备端多项式预优共轭梯度算法进行计算；将计算结果从GPU上复制回CPU端；释放CPU与GPU显存。本发明基于NVIDIA CUBLAS库以及CUDA语言，结合对角预优的预处理方法的优点，对三对角刚度矩阵先进行了一种较好的分裂，避免了对矩阵的求逆运算，从而大大减少了共轭梯度算法的运算量，从而在GPU上实现了刚度矩阵的多项式预处理。

Description

一种基于GPU端的有限元刚度矩阵模拟方法

技术领域

本发明属于有限元的数值模拟计算与计算机应用领域，涉及预优有限元刚度矩阵的加速求解实现。

背景技术

计算机图形处理单元(简称GPU)运用于通用计算的研究越来越多，特别是在大规模科学和计算领域，因为最初设计用于图形图像处理，所以GPU与生俱来就是一种拥有大量运算单元的并行处理器，并且由GPU提供相同的计算能力，所需要的成本和功耗都要小于基于CPU的系统。

作为一种数值模拟与仿真的基本方法，有限元法以其高度的适应性，成为现代工程设计和结构分析的重要方法之一，并在土木、水利、汽车、机械、航空航天、核工业和大地勘测等众多领域应用广泛。随着科学技术的不断发展，工程问题的规模和复杂程度相应提高，也对有限元计算提出更大规模、更快速度的要求。有限元法的基本思想是“化整为零，积零为整”，与并行计算技术的“分而治之”的基本原则相协调。因此，对于大规模的有限元结构分析，研究基于GPU的加速计算方法具有重要意义。其中，刚度矩阵的求解是加速计算的核心技术。

目前，基于矩阵不完全分解预处理的共轭梯度算法在GPU上的实现加速的研究成果相对较多。例如针对该算法中的三角刚度阵求解的串行性进行优化，采用分层调度的方法提高了其并行性。又如，基于对Krylov子空间方法的研究，在数据存储与迭代计算方面进行创新，实现了Krylov子空间方法的GPU加速实现。目前还有一些研究指出了稀疏刚度方程在GPU上的共轭梯度法求解，并结合空间桁架问题验证了CUDA平台GPU有限元计算的加速性。以及，针对CPU之上的向量内积计算时间问题，提出了一种GPU加速向量乘积的归约策略针对CPU之上的向量内积计算时间问题，提出了一种GPU加速向量乘积的归约策略，并基于OpenMP和MPI编程模型实现了加速计算。

但是，涉及多项式预处理的共轭梯度算法研究相对较少。现有的刚度方程求解采用的方法都是将其离散化为高阶线性方程组，从而将原有问题转化为高阶线性方程组的求解问题。共轭梯度算法作为迭代算法中最为有效的方法，深受研究者的关注。但是在具体实现过程中发现，只有当系数矩阵仅有少数几个互不相同的特征值或者非常良态时，共轭梯度算法才能收敛的非常之快。因此，采用预处理技术的共轭梯度算法在求解刚度方程中，具有很好的适用性。预处理技术一般采用不完全的Cholesky因子预优方法。这种方法虽然是一种非常重要的预优技巧，却也有着明显的缺点：该预优算法需要求解两个三角形方程组，并行效率非常低，无法充分利用GPU加速处理器的并行性能。多项式预优方法由于该算法中只包含矩阵的乘法计算和求逆计算，相对不完全的Cholesky分解预处理方法，具有更好的并行性能。然而求逆运算在具体实现过程中仍然是GPU加速计算的瓶颈。而且，即使是大规模的矩阵、向量的乘法运算，在CPU上实现多项式预优算法也会因为耗时巨大，而无法得到有效应用。

发明内容

本发明所要解决的技术问题是：针对现有对不完全分解法预处理技术的缺陷，以及对多项式预处理方法的研究不足，提供一种基于GPU端的有限元刚度矩阵模拟方法，本发明基于NVIDIA CUBLAS库以及CUDA语言，结合对角预优的预处理方法的优点，对三对角刚度矩阵先进行了一种较好的分裂，避免了对矩阵的求逆运算，从而大大减少了共轭梯度算法的运算量，从而在GPU上实现了刚度矩阵的多项式预处理。

本发明解决上述技术问题所采用的技术方案是：

一种基于GPU端的有限元刚度矩阵模拟方法，包括以下步骤，

(1)建立刚度方程Ax＝b，并进行初始化，其中A表示刚度矩阵；

(2)刚度矩阵A的分裂处理；

(3)对应分配GPU显存；

(4)将CPU端初始化数据复制到GPU上；

(5)调用计时函数开始计时；

(6)调用设备端多项式预优共轭梯度算法进行计算；

(7)将计算结果从GPU上复制回CPU端，用于前端显示；

(8)释放CPU与GPU显存。

进一步地，初始化的实现方法是：

A1.使用CPU端malloc函数对矩阵A，分裂阵M1、N1以及中间矩阵进行动态内存分配以及初始化；

A2.使用CPU端malloc函数对向量动态内存分配以及初始化；

A3.构造单位阵

进一步地，刚度矩阵分裂处理的实现方法是：

B1.选取M1为刚度矩阵A的对角预优矩阵；

B2.对M1的对角线元素取倒数，复制回M1，以节约内存；

B3.构造N1矩阵，使N1＝A-M1。

进一步地，分配GPU显存的实现方法为：使用cudaMalloc函数为上述矩阵和向量分配GPU端显存。

进一步地，将CPU端初始化数据复制到GPU上的实现方法为：使用cublasSetVector函数进行将CPU端初始化数据复制到GPU上。

进一步地，调用计时函数开始计时的实现方法为：

E1.使用CUDA事件来对GPU段运算进行计时，

E2.将待处理数据指针传输到PPcg_device函数的参数里。

进一步地，调用设备端多项式预优共轭梯度算法进行计算的实现方法为：

F1.使用cudaMalloc函数为设备端中间矩阵temp、J、G以及中间向量temp，r动态分配内存；

F2.初始化参数cublas库函数参数alpha，_alpha以及变量浮点型beta,q0,q1，error_norm；

F3.使用cublasScopy函数将b向量复制到r向量；

F4.使用cublasSgemv计算r＝Ax-r；

F5.使用cublasSgemm函数计算G＝M1*N1；

F6.使用cudaMemcpy函数将I矩阵复制到矩阵J；

F7.使用cublasSdot函数计算向量r与自身的内积并将结果复制到error_norm上；

F8.预优共轭梯度算法；

F9.释放GPU显存。

本发明的有益效果是：图形处理器(GPU)作为一种高度并行化的通用计算处理器，可以很好解决大规模科学计算的速度问题。NVIDIA统一计算架构(C UDA)为实现GPU通用计算提供了高效、简便的方法。因此可以用来解决大规模有限元求解的时间过长问题。CUBLAS作为NVIDIA官方提供的线性代数计算加速库，提供了较好的加速策略，降低加速计算门槛，可移植性好。本发明基于GPU的硬件平台与CUBLAS的软件平台，实现了加速有限元刚度矩阵求解这一技术。求解刚度矩阵的共轭梯度算法及其预处理方法研究成果丰富，但是大部分的研究内容都是不完全分解的预处理方法，因此本发明考虑更具有并行性的多项式预处理方法，并选择系数矩阵合适的一种分裂，从而获得该方法在GPU端的加速。

通过本发明实现了对原有刚度矩阵的10倍以上加速求解，解决了多项式预处理共轭梯度算法计算过于缓慢的问题，充分发挥了该算法具有的并行性能，弥补了当前多项式预优算法的研究不足。在三对角对称刚度方程的求解验证中，该算法在系数矩阵阶次较低时(10阶以下)，加速并不明显；但是当系数矩阵阶数上升至几千阶以上时，将会获得几十倍以上的加速说明该方法特别适合求解大规模刚度矩阵计算。

附图说明

图1是本发明的整体流程图。

图2是本发明中多项式预调件共轭梯度算法流程图。

图3是本发明的1D拉杆模型图。

图4是图3所示的模型通过本发明的方法进行求解后的加速比图表。

具体实施方式

首先对本发明中的一些技术术语进行解释。

计算机图形处理单元(GPU)：图形处理器(英语：Graphics Processing Unit，缩写：GPU)，又称显示核心、视觉处理器、显示芯片，是一种专门在个人电脑、工作站、游戏机和一些移动设备(如平板电脑、智能手机等)上图像运算工作的微处理器。

中央处理器(CPU)：中央处理器(CPU，Central Processing Unit)是一块超大规模的集成电路，是一台计算机的运算核心(Core)和控制核心(Control Unit)。它的功能主要是解释计算机指令以及处理计算机软件中的数据。

有限元法：有限元法(finite element method)是一种高效能、常用的数值计算方法。科学计算领域，常常需要求解各类微分方程，而许多微分方程的解析解一般很难得到，使用有限元法将微分方程离散化后，可以编制程序，使用计算机辅助求解。

Krylov子空间方法：一种求解大型稀疏线性系统的迭代算法

不完全Cholesky分解预条件共轭梯度法：一种求解稀疏对称正定线性方程组的迭代算法，在求解正定线性方程组时Ax＝b时，先对系数矩阵A进行不完全Cholesky分解，然后再使用共轭梯度算法求解。

以下将结合附图对本发明的具体实施方案做进一步详细说明，应当指出的是，实施例只是对本发明的详细阐述，不应视为对本发明的限定。

在有限元分析领域中，机械结构的应力分析成为一个重要的应用领域，如何加速结构应力的求解速度是研究的重点。其本质是研究刚度矩阵的求解问题。但是随着问题的复杂度上升，求解的刚度矩阵规模也更大，条件数也会更大，导致普通的共轭梯度算法收敛速度缓慢。

考虑到机械结构系统离散化的刚度方程一般均为三对角阵，本发明以多项式预处理的共轭梯度算法为基础，以一维拉杆变形问题为例(如图3)，建立其平衡方程及边界条件，形如

u(x)|_x＝0＝0

其中，u为位移场函数，p为均布载荷，A为横截面面积，E为材料弹性模量，L为拉杆长度。根据有限元理论，对平衡方程离散成1000阶线性方程组，即

为显示本发明的优越性以及可推广性，忽略实际意义，对

随机取值，并充分结合对角预优方法的特点以及GPU并行特点，改善多项式预优算法，优化刚度方程求解。

本发明的技术方案主要通过以下方法实现：

A.刚度方程Ax＝b的初始化，其中

x为待求解向量，初始化为0，

其中取h＝0.002，具体实现方法为：

A1.使用CPU端malloc函数对矩阵A，分裂阵M1、N1以及中间矩阵进行动态内存分配以及初始化，其中分裂阵M1＝0，N1＝0；

A2.使用CPU端malloc函数对向量动态内存分配以及初始化；

A3.初始化单位阵

B.刚度矩阵A的分裂处理，实现方法为：

B1.选取M1为系数矩阵A的对角预优矩阵；

B2.对M1的对角线元素取倒数，复制回M1，

以节约内存；

B3.构造N1矩阵，使N1＝A-M1,

C.对应分配GPU显存，实现方法为：使用cudaMalloc函数为上述矩阵和向量分配GPU端显存。

D.将CPU端初始化数据复制到GPU上，实现方法为：使用cublasSetVector函数进行将CPU端初始化数据复制到GPU上。

调用计时函数开始计时，实现方法为：

E1.使用CUDA事件来对GPU段运算进行计时；

E2.将待处理数据指针传输到PPcg_device函数的参数里。

调用设备端多项式预优共轭梯度算法进行计算，实现方法为：

F1.使用cudaMalloc函数为设备端中间矩阵temp、J、G以及中间向量temp，r动态分配内存；

F2.初始化参数cublas库函数参数alpha＝1.0，_alpha＝-1.0以及变量浮点型beta＝0.0,q0＝0.0,q1＝0.0，error_norm＝0；

F3.使用cublasScopy函数将b向量复制到r向量，

F4.使用cublasSgemv计算r＝Ax-r；

F5.使用cublasSgemm函数计算G＝M1*N1；

F6.使用cudaMemcpy函数将I矩阵复制到矩阵

F7.使用cublasSdot函数计算向量r与自身的内积并将结果复制到error_norm上；

F8.预优共轭梯度算法，具体实现方法如图2所示：

F9.释放GPU显存。

G.计时终止。

H.将计算结果从GPU上复制回CPU端，用于前端显示，具体实现方法：使用cudaGetVector函数复制会计算结果。

I.释放CPU与GPU显存，程序终止。

Claims (1)

1.一种基于GPU端的有限元刚度矩阵模拟方法，针对弹性体变形的求解问题，首先依据静力学原理建立其平衡方程，并基于有限元方法对其动力学方程进行离散，从而将微分方程求解转化为线性方程组的求解，该线性方程组即为需要求解的刚度方程Ax＝b，其中A为刚度矩阵，x为位移向量，b为节点力向量，具体包括以下步骤：

(1)建立变形体的静力学方程，基于有限元方法对其离散成刚度方程Ax＝b

(2)刚度矩阵A的分裂处理；

(3)对应分配GPU显存；

(4)将CPU端初始化数据复制到GPU上；

(5)调用计时函数开始计时；

(6)调用设备端多项式预优共轭梯度算法进行计算；

(7)将计算结果从GPU上复制回CPU端，用于前端显示；

(8)释放CPU与GPU显存；

A.刚度方程Ax＝b的初始化，其中

x为待求解向量，初始化为0，

其中取h＝0.002，具体实现方法为：

A1.使用CPU端malloc函数对矩阵A，基于对角预优方法获得分裂阵M1、N1，并对临时矩阵进行动态内存分配以及初始化，作为中间变量；其中分裂阵M1＝0，N1＝0；

A2.使用CPU端malloc函数对向量动态内存分配以及初始化；

A3.构造单位阵

步骤(2)中，刚度矩阵A分裂处理的实现方法是：

B1.选取M1为刚度矩阵A的对角预优矩阵；

B2.对M1的对角线元素取倒数，复制回M1，

以节约内存；

B3.构造N1矩阵，使

步骤(3)中，分配GPU显存的实现方法为：使用cudaMalloc函数为上述矩阵和向量分配GPU端显存；

步骤(4)中，将CPU端初始化数据复制到GPU上的实现方法为：使用cublasSetVector函数进行将CPU端初始化数据复制到GPU上；

步骤(5)中，调用计时函数开始计时的实现方法为：

E1.使用CUDA事件来对GPU段运算进行计时，

E2.将待处理数据指针传输到PPcg_device函数的参数里；

步骤(6)中，调用设备端多项式预优共轭梯度算法进行计算的实现方法为：

F1.使用cudaMalloc函数为设备端中间矩阵temp、J、G以及中间向量temp，r动态分配内存；

F2.初始化参数cublas库函数参数alpha＝1.0，_alpha＝-1.0以及变量浮点型beta＝0.0,q0＝0.0,q1＝0.0，error_norm＝0；

F3.使用cublasScopy函数将b向量复制到r向量，

F4.使用cublasSgemv计算r＝Ax-r；

F5.使用cublasSgemm函数计算G＝M1*N1；

F6.使用cudaMemcpy函数将I矩阵复制到矩阵

F7.使用cublasSdot函数计算向量r与自身的内积并将结果复制到error_norm上；

F8.预优共轭梯度算法；

F9.释放GPU显存；

步骤(7)的具体实现方法：使用cudaGetVector函数复制会计算结果。

Images