英语原文共 24 页，剩余内容已隐藏，支付完成后下载完整资料
基于计算机模拟的整体式有限元矩阵构建
Francisco Javier Ramirez-Gila，Marcos de Sales Guerra Tsuzukib，Wilfredo
Montealegre-Rubioa,*
a 哥伦比亚国立大学矿业学院机械工程系，麦德林，050034，哥伦比亚
b 圣保罗大学政治学院机电与机械系系统工程计算几何实验室，梅洛莫拉斯，2231圣保罗，巴西
摘要：有限元法（FEM）有数个计算步骤，以数值求解特定的问题，其中为了加速线性方程组，已经做了许多工作。但是，对非结构化网格来说同样耗时的有限元矩阵构造，却没有得到足够的研究。整体式有限元矩阵的产生分为两个步骤：首先由数值积分计算局部矩阵，然后将它们组装到一个整体式的系统中，这两个步骤按照传统方法是利用串行计算进行的。本论文展示了一种快速构建整体式有限元矩阵的方法，这种方法通过在三维领域解泊松方程产生。由于其内在的并行机会，在图形处理单元（GPU）中计算数值积分，并且由于其内在的串行操作，在中央处理单元（CPU）中计算矩阵组件。在数值积分中，由于其对称性，每个局部强度矩阵中只有下三角部分被计算出来，这节省了GPU内存和计算时间。 由于对称性，整体式稀疏矩阵只在其下三角部分包含非零元素，这减少了组装操作和内存使用。 这种方法允许从GPU上的任何非结构化有限元网格尺寸生成整体式稀疏矩阵，只有很少的内存容量，仅受CPU内存的限制。
关键词：有限元法(FEM) 非结构化网格 矩阵生成 并行计算 图形处理单元（GPU） 异构计算
1.介绍
静止状态下的许多物理现象，如电势和磁势，热传导，流体流动和静态条件下的弹性问题，可以用椭圆偏微分方程来描述(EPDEs)。EPDE不涉及时间变量，并描述了稳态响应。线性EPDE的一般形式如式（1）所示，a，b，c是系数函数，f是源（激励）函数，u是未知变量，所有这些函数都可以在空间上变化（x，y，z）。
(1)
EPDE可以通过数学程序（如傅立叶级数）精确求解[1]。但是，经典解并不经常存在，如果允许使用这些分析方法，问题会有许多简化[2]。因此，已经开发了几种数值方法，如有限元法（FEM）和有限差分法来有效地解决EPDEs。
有限元法与其他方法相比有几个优点。主要优点是对于使用非结构化网格的复杂几何体的问题特别适用[2]。有一种获得解决EPDEs问题的合适框架的方法，那就是通过使用有限元法，将它们表述为变分问题，也称为弱解。
EPDE的变分公式是将强配方转化为弱配方的数学处理，允许元素或子域中的近似值，并且使用所谓的Galerkin方法来解决EPDE变分形式[3]。图1总结了使用FEM数值求解EPDE所需的步骤。有些软件是公开可用的，所有这些步骤都是在计算机中自动求解的[4]，然而，最常见的商业有限元软件主要包含步骤4-8，它们是：
(i) 有限元（FE）的域离散化;
(ii) 所有FE的数值积分。在该步骤中，针对每个有限元素计算局部矩阵ke和局部加载向量fe;
(iii) 局部载荷向量fe来自于局部矩阵ke，全局稀疏矩阵K来自于全局载荷向量F;
(iv) 应用边界条件;
(v) 先前形成的线性方程组的解，KU = F，U表示节点解的向量。
图1：解决EPDE的有限元步骤
通过使用直接[5]和迭代方法[6]已经做出许多有效的加速线性方程求解器；然而，我们很少研究有限元稀疏矩阵结构（图1的步骤5和6），这对于密集的非结构化网格也是耗时的 [7]。因此，在本文中，整体式稀疏矩阵的快速构建被开发出来，这是由有限元法在EPDE数值解中产生的。有效构建全局FE矩阵的一种方法是使用并行计算。并行性是当前计算的趋势，因此，微处理器制造商专注于增加内核，而不是增加单线程性能[8]。虽然有几个并行架构，最普遍的多核架构是GPU（图形处理单元），最常见的多核架构是CPU（中央处理器）[9]。但是，由于GPU与CPU相比具有巨大的并行性，科学计算的趋势是使用GPU加速密集并行计算[10]。此外，GPU比其他架构更受青睐，因为其市场在计算机游戏行业中已得到很好的确立，所以被用作协处理器来加速部分或整个代码。在大多数个人电脑中具有图形处理器已成为事实，这使GPU成为应用程序开发人员的经济和有吸引力的选择[10]。
虽然一些研究主要集中在将有限元子程序移植到GPU中[11,12]，其他作品则侧重于将部分FEM代码移植到GPU上，如数值积分[13,14,15,16]，矩阵集合[17,18]，以及方程的大线性系统的解[19,20,21]。尽管以前的一些工作加速了稀疏FE矩阵的构建[7,22,23]，这些方法需要高GPU存储器，这意味着可以分析的问题的大小是有限的，或者需要昂贵的具有高存储容量的GPU。 但是，在大多数配备图形处理器的个人电脑中，GPU具有低内存容量。
因此，所提出的方法基于这样的假设：单个GPU比CPU具有更少的内存; 从而，开发了CPU-GPU并行。 核心思想是计算GPU中的数值积分，这是一个密集并行子程序，并且在CPU中计算组件，这是一个占用大量内存的主要串行子程序。因此，基于这种方法，我们能够构造有限元中产生的任何大型稀疏矩阵，只受CPU内存限制。本文的其余部分安排如下。 在第2节中，给出了用于检查所提议的实施的示例的泊松方程的有限元公式。 第3节介绍了FE矩阵结构的CPU-GPU实现。 结果将在第4节讨论，最后在第5节中得出结论。
2. 有限元椭圆PDE建模
通过求解三维（3D）域中的泊松方程来显示我们在构建大型稀疏有限元矩阵中所提出实现的计算好处。 泊松方程是一个椭圆型PDE。 这种类型的方程通常用于解决标量问题，如电势和磁势，传热和其他无外部流动的问题[3]。 另外，泊松方程用于建模矢量场问题，如无外力的结构问题。 所有这些问题应该处于静态或静止状态。 在本节中，从有力的制定步骤到构建全局稀疏矩阵步骤（图1的步骤1到6），提出了FEM步骤。
泊松方程的强表述由[3]给出：
(2)
其中表示域上的一般标量变量，即可以是电势或温度，这取决于特定的应用。位置矢量是r并且项c可以表示各向同性材料的电导率或热导率。术语 是边界处的标量变量（Dirichlet条件）的规定值，并且流量（或通量）qn的法向分量为了简单起见仅在边界（Neumann条件）下表示为零，仅用于计算评估的目的。在通过方程式（2）的部分进行整合之后，弱形式被发现为：
（3）
其中是一个任意函数[3]。基于泊松方程的弱形式，可以找到有限元公式。 相应地，域被离散化为FE; 接下来，在每个FE中，标量变量近似为：
（4）
其中n是元素中的节点数量，根据Galerkin原理，是形状函数。用方程 （4）代入方程 （3）并评估所有元素的积分，得到一个方程组如：
（5）
其中K是所谓的全局稀疏矩阵或强度矩阵。变量是节点解向量，全局载荷向量F由边界修改条件，不再是零矢量。K矩阵是作为每个FE的数值积分中获得的矩阵贡献的总和而获得的。 这一步被称为程序集[3]：
（6）
其中nel是FE的数量，ke是元素e（局部刚性矩阵）的局部矩阵，可以得到：
（7）
在我们的方法中，我们认为，作为材料属性数组（MP）的一部分，系数c可以从一个元素到另一个元素不同，但对于一个特定元素它仍然是常量，这意味着它可以在公式（7）积分之外。全局稀疏矩阵结构在方程 （6）和方程 （7）。然而，为了计算效率，公式（7）是通过数值求解而不是精确积分。 有几个数值程序可以求解积分，但在有限元中最常见的是Gauss-Legendre或简单的高斯积分[3]。 对于使用高斯积分，必须应用使用参考元素技术的变换，其中物理坐标（x，y，z）被转换为自然坐标（r，s，t）。在局部矩阵计算中的这些坐标系之间的等价（见方程（7））表示为：
（8）
其中J是雅可比矩阵，而项|J|是其行列式。 之后，将高斯积分应用于归一化域（等式（8）的右侧）如下：
（9）
其中Mi，Mj，Mk，ri，sj，tk是3D域中每个方向上积分点的数量，点和权重。在我们的实现中，我们使用具有8个节点的六面体FE来离散化该域，也被称为8节点砖单元，其是低阶三线性FE。8节点砖元素的选择基于大多数FE代码基于低阶元素的想法，并且很少有作品集中在这类元素上[13]。在[3]中可以参考8节点砖元素的精确数值积分的点和权重。
3. 数值实现
构造全局稀疏FE矩阵的典型方法是使用串行计算。 有几种方法来构造稀疏矩阵; 然而，最快的方法是首先创建triplet1表单，然后使用稀疏函数来压缩三重稀疏格式，从而避免一些内存问题[24]。尽管三元组格式很容易创建，但它很难用于大多数稀疏矩阵算法; 因此需要压缩稀疏格式[5]。常见的压缩稀疏格式是CSC（压缩稀疏列），在MATLAB中默认使用[25]。算法1展示了使用MATLAB语法生成稀疏矩阵的典型实现。该算法需要一个子程序，用于计算局部矩阵ke和本地加载矢量fe的数值积分。但是，这项工作只关注生成矩阵K;负载向量F不考虑用于并行计算（该术语没有考虑，因为它在计算上要求不高[15,16]）。数值积分子程序随PDE，元素类型和基函数而变化;它是节点坐标（C）和单元e [17]的力，边界条件（BC）和材料属性（MP）等任何节点域的函数。在计算存储在数组Kall中的局部矩阵ke之后，需要从局部到全局自由度（DOF）的映射。该映射是考虑连通性矩阵（E），该连通性矩阵给出构成元素的全局节点编号以及每个节点的自由度数（DOFxn）。作为映射函数的结果，获得ke的每个条目的全局稀疏矩阵中的位置，即分别存储在阵列Krow和Kcols中的行（Indrow）和列（Indcol）索引。最后，全局稀疏矩阵是在汇编函数的单个调用（例如MATLAB稀疏函数[5,24]）中获得的。
已经提出了几种方法来并行化全局稀疏矩阵结构。在这项工作中，不是将矩阵生成子程序看作一个单独的块，而是将它分成两个子程序：数值积分子程序和汇编子程序，它允许分析寻找并行机会。数值积分子程序有几种并行机会[16]：
1.并行计算所有高斯积分点;
2.平行化每个ke组件[15,16];
3.最后，也可以并行计算所有元素[13,16]。
集成点上的循环（第一种方法）对于高阶FE是一种很好的策略，但由于固有数据race2更新ke中的条目和来自不同集成点的贡献，所以更难并行化，并且通常并发度 较低[16]。 第二种方法对高阶有限元也有帮助，但对低阶元不方便。低阶元素（第三种方法）的循环并行化是最自然的方式，因为线程可以完全由一个局部矩阵ke [16,17]来充电。因此，在我们的实现中选择了第三种方法，因为使用前两种方法8节点砖FE可以实现较低的并行度。
在ke的计算中，不存在相关性，只需要来自特定8节点砖元素的信息。 因此，不同元素的数值积分也是完全可并行化的[15]。一个线程用于计算每个局部矩阵ke。 然后，使用线程块的一维网格，可以计算网格中的所有FE。
在用于标量问题（每个节点1个DOF）的八节点砖单元中，每个ke是具有64个（8times;8）分量的正方形密集矩阵，但是利用矩阵对称性，只有下或上三角矩阵部分（包括其对角线）应该被考虑。这导致浮点操作和内存需求的显着节省，因为大约需要一半局部矩阵组件，特别是36个组件，而不是64个。另外，对于相同的FE，只需要进行一次计算 类型，例如形状函数及其在物理坐标中的衍生物。 我们使用相同的FE类型来离散整个域。
算法2给出了CUDA内核来计算局部强度矩阵ke。一个线程负责计算ke的36个分量，并将它们作为Kvals中的向量存储在GPU全局内存中。对于每个积分点，除了梯度矩阵（B）和局部僵硬度之外，每个线程还必须连续执行多个运算，例如雅可比矩阵（J）计算，其行列式（|J|）及其逆（J-1）有限矩阵（ke）。在每个积分点评估的形状函数（N）及其导数相对于自然坐标（dN）仅计算一次，并且它们存储在快速GPU存储器中。矩阵J的计算需要三个密集的矩阵 - 向量乘法，即dN（3times;8）乘以笛卡尔坐标（8times;3）中的节点坐标的乘法。|J|的计算是用3x3矩阵的行列式和J-1来执行的，公式用于反转相同大小的矩阵。矩阵B被计算为J-1（3x3）乘以dN（3x8）。最后，元素矩阵ke按照等式（9）以及一些修改只计算下三角矩阵分量。
然而，如果FE网格更密集，数值积分可能不会被执行，因为这是一个耗费大量内存的操作，GPU内存不能完全分配[17]。因此，使用基于GPU存储器可用性的简单域分区，这是由Komatitsch等人提出的修改实现[11]。因此，数值积分是在FE子群（GPUg）中开发的，即数值积分内核被调用来操作FEs集。 元素子组由GPU顺序执行，但是一组内的FE上的数值集成并行执行。采用这种方法，全局FE矩阵的构造可以由任何网格尺寸构成
全文共13476字，剩余内容已隐藏，支付完成后下载完整资料
资料编号：[11188]，资料为PDF文档或Word文档，PDF文档可免费转换为Word
以上是毕业论文外文翻译，课题毕业论文、任务书、文献综述、开题报告、程序设计、图纸设计等资料可联系客服协助查找。
