CN109753682B - 一种基于gpu端的有限元刚度矩阵模拟方法 - Google Patents
一种基于gpu端的有限元刚度矩阵模拟方法 Download PDFInfo
- Publication number
- CN109753682B CN109753682B CN201811439272.6A CN201811439272A CN109753682B CN 109753682 B CN109753682 B CN 109753682B CN 201811439272 A CN201811439272 A CN 201811439272A CN 109753682 B CN109753682 B CN 109753682B
- Authority
- CN
- China
- Prior art keywords
- gpu
- matrix
- function
- cpu
- vector
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 68
- 239000011159 matrix material Substances 0.000 title claims abstract description 57
- 238000012545 processing Methods 0.000 title claims abstract description 16
- 238000004088 simulation Methods 0.000 title claims abstract description 8
- 230000006870 function Effects 0.000 claims abstract description 43
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 29
- 238000005457 optimization Methods 0.000 claims abstract description 22
- 238000004364 calculation method Methods 0.000 claims abstract description 20
- HPTJABJPZMULFH-UHFFFAOYSA-N 12-[(Cyclohexylcarbamoyl)amino]dodecanoic acid Chemical compound OC(=O)CCCCCCCCCCCNC(=O)NC1CCCCC1 HPTJABJPZMULFH-UHFFFAOYSA-N 0.000 claims abstract description 5
- 239000013598 vector Substances 0.000 claims description 30
- 238000006073 displacement reaction Methods 0.000 claims description 2
- 230000010076 replication Effects 0.000 claims description 2
- 230000003068 static effect Effects 0.000 claims 2
- 230000005489 elastic deformation Effects 0.000 claims 1
- 238000007781 pre-processing Methods 0.000 abstract description 11
- 230000001133 acceleration Effects 0.000 description 10
- 238000011160 research Methods 0.000 description 9
- 238000000354 decomposition reaction Methods 0.000 description 6
- 238000005516 engineering process Methods 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 2
- 238000002939 conjugate gradient method Methods 0.000 description 2
- 238000012916 structural analysis Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004883 computer application Methods 0.000 description 1
- 238000013500 data storage Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000002203 pretreatment Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
Landscapes
- Complex Calculations (AREA)
Abstract
本发明提供一种基于GPU端的有限元刚度矩阵模拟方法,包括:建立刚度方程Ax=b并初始化;分裂处理刚度矩阵A;对应分配GPU显存;将CPU端初始化数据复制到GPU上;调用计时函数开始计时;调用设备端多项式预优共轭梯度算法进行计算;将计算结果从GPU上复制回CPU端;释放CPU与GPU显存。本发明基于NVIDIA CUBLAS库以及CUDA语言,结合对角预优的预处理方法的优点,对三对角刚度矩阵先进行了一种较好的分裂,避免了对矩阵的求逆运算,从而大大减少了共轭梯度算法的运算量,从而在GPU上实现了刚度矩阵的多项式预处理。
Description
技术领域
本发明属于有限元的数值模拟计算与计算机应用领域,涉及预优有限元刚度矩阵的加速求解实现。
背景技术
计算机图形处理单元(简称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函数对向量动态内存分配以及初始化;
进一步地,刚度矩阵分裂处理的实现方法是:
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阶线性方程组,即
本发明的技术方案主要通过以下方法实现:
A.刚度方程Ax=b的初始化,其中
A1.使用CPU端malloc函数对矩阵A,分裂阵M1、N1以及中间矩阵进行动态内存分配以及初始化,其中分裂阵M1=0,N1=0;
A2.使用CPU端malloc函数对向量动态内存分配以及初始化;
B.刚度矩阵A的分裂处理,实现方法为:
B1.选取M1为系数矩阵A的对角预优矩阵;
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;
F4.使用cublasSgemv计算r=Ax-r;
F5.使用cublasSgemm函数计算G=M1*N1;
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显存;
A1.使用CPU端malloc函数对矩阵A,基于对角预优方法获得分裂阵M1、N1,并对临时矩阵进行动态内存分配以及初始化,作为中间变量;其中分裂阵M1=0,N1=0;
A2.使用CPU端malloc函数对向量动态内存分配以及初始化;
步骤(2)中,刚度矩阵A分裂处理的实现方法是:
B1.选取M1为刚度矩阵A的对角预优矩阵;
步骤(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;
F4.使用cublasSgemv计算r=Ax-r;
F5.使用cublasSgemm函数计算G=M1*N1;
F7.使用cublasSdot函数计算向量r与自身的内积并将结果复制到error_norm上;
F8.预优共轭梯度算法;
F9.释放GPU显存;
步骤(7)的具体实现方法:使用cudaGetVector函数复制会计算结果。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811439272.6A CN109753682B (zh) | 2018-11-29 | 2018-11-29 | 一种基于gpu端的有限元刚度矩阵模拟方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811439272.6A CN109753682B (zh) | 2018-11-29 | 2018-11-29 | 一种基于gpu端的有限元刚度矩阵模拟方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109753682A CN109753682A (zh) | 2019-05-14 |
CN109753682B true CN109753682B (zh) | 2020-12-22 |
Family
ID=66402553
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201811439272.6A Active CN109753682B (zh) | 2018-11-29 | 2018-11-29 | 一种基于gpu端的有限元刚度矩阵模拟方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109753682B (zh) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112016232B (zh) * | 2020-08-31 | 2024-06-14 | 中国原子能科学研究院 | 一种撕裂有限元过程处理方法及系统 |
CN112084650A (zh) * | 2020-09-04 | 2020-12-15 | 杭州百子尖科技股份有限公司 | 基于cuda的化工流程模拟软件提高计算速度的方法 |
CN117473212B (zh) * | 2023-12-27 | 2024-04-16 | 粤港澳大湾区数字经济研究院(福田) | Ntt算法的gpu加速方法、装置、设备及存储介质 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7321365B2 (en) * | 2002-05-31 | 2008-01-22 | Siemens Product Lifecycle Management Software Inc. | Computerized deformation analyzer |
CN105335332A (zh) * | 2015-12-07 | 2016-02-17 | 郑州航空工业管理学院 | 特殊鞍点问题的高效预处理方法 |
CN105808926A (zh) * | 2016-03-02 | 2016-07-27 | 中国地质大学(武汉) | 一种基于gpu并行加速的预条件共轭梯度区域网平差方法 |
CN106126823A (zh) * | 2016-06-23 | 2016-11-16 | 广州中国科学院工业技术研究院 | 一种基于提高迭代法稳定性和收敛性的位移求解方法 |
CN106570204A (zh) * | 2016-09-23 | 2017-04-19 | 西安交通大学 | 一种基于cpu+gpu异构并行计算的透平机械叶片静强度特性分析方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106777769B (zh) * | 2017-01-08 | 2018-04-24 | 浙江大学 | 预测低速冲击下复合材料多层厚板渐进失效的有限元方法 |
-
2018
- 2018-11-29 CN CN201811439272.6A patent/CN109753682B/zh active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7321365B2 (en) * | 2002-05-31 | 2008-01-22 | Siemens Product Lifecycle Management Software Inc. | Computerized deformation analyzer |
CN105335332A (zh) * | 2015-12-07 | 2016-02-17 | 郑州航空工业管理学院 | 特殊鞍点问题的高效预处理方法 |
CN105808926A (zh) * | 2016-03-02 | 2016-07-27 | 中国地质大学(武汉) | 一种基于gpu并行加速的预条件共轭梯度区域网平差方法 |
CN106126823A (zh) * | 2016-06-23 | 2016-11-16 | 广州中国科学院工业技术研究院 | 一种基于提高迭代法稳定性和收敛性的位移求解方法 |
CN106570204A (zh) * | 2016-09-23 | 2017-04-19 | 西安交通大学 | 一种基于cpu+gpu异构并行计算的透平机械叶片静强度特性分析方法 |
Non-Patent Citations (3)
Title |
---|
"基于GPU并行计算及CUDA编程在环境工程中的应用研究";胡兵;《中国优秀硕士学位论文全文数据库 基础科学辑》;20180215;第38-44页 * |
"有限元GPU加速计算的实现方法";张健飞,沈德飞;《计算机辅助工程》;20140430;第23卷(第2期);第41-45页 * |
张健飞,沈德飞."有限元GPU加速计算的实现方法".《计算机辅助工程》.2014,第23卷(第2期), * |
Also Published As
Publication number | Publication date |
---|---|
CN109753682A (zh) | 2019-05-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Fischer | Scaling limits for PDE-based simulation | |
CN109753682B (zh) | 一种基于gpu端的有限元刚度矩阵模拟方法 | |
Chen et al. | An escheduler-based data dependence analysis and task scheduling for parallel circuit simulation | |
Fatemi et al. | AMITIS: A 3D GPU-based hybrid-PIC model for space and plasma physics | |
Ahamed et al. | Iterative methods for sparse linear systems on graphics processing unit | |
Martínez-Frutos et al. | Efficient matrix-free GPU implementation of fixed grid finite element analysis | |
Jespersen | Acceleration of a CFD code with a GPU | |
He et al. | A multiple-GPU based parallel independent coefficient reanalysis method and applications for vehicle design | |
Sooknanan et al. | GPU computing using CUDA in the deployment of smart grids | |
Zhang et al. | Gpu-based implementation of finite element method for elasticity using cuda | |
Xu et al. | Towards a scalable hierarchical high-order CFD solver | |
Wang et al. | A novel parallel algorithm for sparse tensor matrix chain multiplication via tcu-acceleration | |
CN116167304B9 (zh) | 基于神威架构的油藏数值模拟gmres优化方法及系统 | |
Zheng et al. | GPU-based multifrontal optimizing method in sparse Cholesky factorization | |
Oyarzun et al. | Direct numerical simulation of incompressible flows on unstructured meshes using hybrid CPU/GPU supercomputers | |
Lastovetsky | Heterogeneous parallel computing: from clusters of workstations to hierarchical hybrid platforms | |
Zou et al. | Supernodal sparse Cholesky factorization on graphics processing units | |
Jung et al. | Accelerating implicit integration in multi-body dynamics using GPU computing | |
Khan et al. | Analyzing the Implementation of the Newton Raphson Based Power Flow Formulation in CPU+ GPU Computing Environment | |
Zhang et al. | Implementation and efficiency analysis of parallel computation using OpenACC: a case study using flow field simulations | |
Luo et al. | GPU accelerated cell-based adaptive mesh refinement on unstructured quadrilateral grid | |
Malaya et al. | Accelerating matrix processing with GPUs | |
Li et al. | Nonlinear Dynamic Analysis Efficiency by Using a GPU Parallelization. | |
Wang et al. | An efficient architecture for floating-point eigenvalue decomposition | |
Aliaga et al. | Harnessing CUDA dynamic parallelism for the solution of sparse linear systems |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |