基于GPU-Warp的着色法有限元矩阵生成与装配

可在ScienceDirect上获得目录列表计算设计与工程杂志首页：www.elsevier.com/locate/jcde计算设计与工程学报6（2019）705基于GPU-Warp的着色法有限元矩阵生成与装配Utpal Kiran，Deepak SharmaChang，Sachin Singh Gautam印度理工学院机械工程系，Guwahati，Assam 781039，印度阿提奇莱因福奥文章历史记录：在线发售2018年保留字：有限元法数值积分装配GPUCUDA着色方法A B S T R A C T有限元法已成功地实现了图形处理单元，以实现显着减少模拟时间。本文提出了一种新的有限元矩阵生成策略，包括数值积分和装配。这些策略是使用众所周知的着色方法。所提出的策略使用一个专门的算法来实现细粒度的并行和有效地利用片上存储器资源。计算统一设备架构（CUDA）的warp shuffle功能用于加速数值积分。采用部分并行的数值积分方法，进一步优化了单元刚度矩阵的求解。用每个节点具有三个自由度的8节点六面体单元对三维弹性问题进行了性能评估在NVIDIA Tesla K40 GPU上，我们获得了高达8.2倍的加速比，超过了基于着色的元素组装策略（每个元素使用一个线程此外，所提出的策略实现更好的算法吞吐量和带宽。©2018计算设计与工程学会Elsevier的出版服务这是一个开放在CC BY-NC-ND许可证（http://creativecommons.org/licenses/by-nc-nd/4.0/）下访问文章1. 介绍在通用图形处理单元（GPGPU）上的计算已被证明在加速计算密集型以及存储器密集型科学代码方面非常有效GPU与CPU（中央处理器）协同工作，以加速世界各地的科学、分析、工程和企业应用GPU和CPU的架构设计存在根本差异GPU基本上是一种多核处理器，具有数千个设计用于提供高计算吞吐量的内核，而CPU牺牲处理器的计算吞吐量来提高单核的性能这种不同的设计理念需要开发专 门 的 算 法 ， 以 利 用 GPU 硬 件 的 潜 在 性 能 （ Garland Kirk ，2010）。GPU最适合于面向吞吐量的数值应用，这些应用对大量独立的数据集进行操作。有限元法是求解偏微分方程组的最常用的数值方法由计算设计与工程学会负责进行同行评审。*通讯作者。电 子 邮 件地 址 ：ukiran@iitg.ernet.in（ 美 国 ）Kiran ） ，dsharma@iitg.ernet.in（D.Sharma），ssg@iitg.ernet.in（S.S. Gautam）。选项。对于具有数千或数百万个网格单元的复杂问题，FEM模拟是高度计算和内存密集型的（Filipaly，Cardellini，Barbieri，&Fanfarillo，2017; Mossaiby，Rossi，Dadvand，&Idelsohn，2012）。对于3D（三维）中的问题，计算负荷进一步增加。FEM主要涉及计算网格的所有元素的元素矩阵和力矢量，以及将它们组装成线性方程组对于具有复杂弱形式或高阶元素的3D问题，数值积分和装配时间可能是相当大的。在需要迭代组装的问题中，如非线性弹性，性能可能会降低更多GPU在加速基于FEM的仿真的几乎每一步方面都非常有效（Georgescu，Chow，&Okuda，2013）。几个数量级的加速在用于不同应用的元素矩阵的组装中被报道（Cecka，Lew，&Darve，2011; Fu ， Lewis ， Kirby ， Whitaker ， 2014; Komatitsch ，Michéa，Erlebacher，2009）。虽然数值积分不那么容易并行化，但 Macioclave ， Pjaszewski 和 Banas ' （ 2010 ） ， Dziekonski ，Sypek ， Lamecki 和 Mrozowski （ 2012 ） 以 及 Zhang 和 Shen（2013）中观察到了显着的加速。本文旨在探索一种更有效的GPU实现三维弹性力学问题数值积分和装配的策略。元素矩阵的计算https://doi.org/10.1016/j.jcde.2018.11.0012288-4300/©2018计算设计与工程学会Elsevier的出版服务这是一个在CC BY-NC-ND许可证下的开放获取文章（http://creativecommons.org/licenses/by-nc-nd/4.0/）。××××706U. Kiran等人/Journal of Computational Design and Engineering 6（2019）705-718在GPU上的是有问题的，因为它由许多连续的步骤组成。我们提出的汇编策略仔细使用了一个特殊的算法，合适的数据结构，并在编程语言中最近添加的功能，以获得更好的性能。基于着色法，给出了八节点六面体单元的三种不同装配策略。所有策略都为网格的每个元素分配一个扭曲，以进行完全或部分的数值积分和组装CUDA的经纱洗牌功能已被用于所有的战略，以计算雅可比矩阵在一个非常有效的方式。通过预先计算索引，以CSR（压缩稀疏行）格式完成对全局矩阵的组装完成所提出的策略的剖析，并且将性能与所建立的按元素策略（每个元素单个线程）的组装进行比较（Cecka等人，2011;Fu等人，2014; Reguly& Giles，2013）。本文的组织结构如下。第2节处理的细节，以前的作品，解决了问题，ING有限元GPU上。在第3节中，介绍了CUDA编程模型的关键组件和每个线程可用的内存层次结构。讨论了三维弹性力学问题的有限元列式，以及有限元汇编的串行和并行实现并行实现是基于着色方法，并使用组件的策略（每个元素的单线程）进行计算。第4节提出了新的战略。在第5节中，详细分析了所提出的策略，并将性能相互比较，本文在第6节中结束。2. 以前的作品在GPU上实现有限元模拟的早期研究主要集中在稀疏线性方程组的解决方案这是因为求解器是FEM中最耗时的步骤，并且可以并行完成 有大量的文献讨论了GPU上稀疏线性解算器的许多有效实现（Bolz，Farmer，Grinspun，&Schröoder，2003; Fu等人，2014;Li &Saad，2013; Reguly &Giles，2013）。 几种有效的文库，如CUSP（Dalton，Bell，Olson，&Garland，2014），MAGMA（Bosma，Cannon，&Playoust，1997），AmgX（Naumov等人， 2015）等。也存在。硬件的新功能和更灵活的编程环境帮助研究人员专注于有限元的其他步骤，如元素矩阵的计算和装配。FEM的这 些 步 骤 对 于 涉 及 更 精 细 网 格 （ München &Sirouspour ，2014）、更高阶有限元（Komatitsch等人， 2009），非线性弹性（ Cai ， Wang ， Li ， &Wang ， 2015 ） 和 拓 扑 优 化 （ Ram&Sharma，2017）仅举几例。在Maciovirk et al. （2010），讨论了电磁场数值积分的GPU实现。这里，单个有限元的计算被分配给单个线程块。线程块的各个线程使用积分点上的一个外部循环和形状函数上的两个内部循环来计算元素刚度矩阵的一组条目如果问题很大，则为整个网格或子网格启动线程块网格。实验结果表明，该算法的加速比为3对于NVIDIA GeForce 8800 GTX上的高阶棱柱单元，为19。同一作者的另一项研究调查了使用OpenCL对具有弯曲几何形状的四边形单元的数值积分，Pamaszewski，Banas'和Maciovich（2010）。在Dziekonski等人（2012）中实现了高阶四面体单元的数值积分，其中并行计算了32个单元的单元矩阵。每个元素使用81个线程块来执行81个高斯点的计算。此后，使用不同的核来执行元件刚度矩阵的组装有一些实现其中集成和组装由单个GPU内核执行（Zhang Shen，2013）。使用众所周知的着色方法GPU内核然后对相同颜色的元素执行元素刚度计算，然后将它们组装到全局刚度矩阵中。研究显示节点坐标的合并访问，但没有讨论数值积分的详细实现。该研究表明，分别在NVIDIA GeForce GT 430上使用四边形和六面体（8节点）单元。在Ohshima，Hayashi，Katagiri，andNakajima（2012）中发现了使用着色方法的类似组装策略。然而，作者提出了两种类型的实现的元素矩阵计算的8节点的三维元素。一个使用单个内核进行所有计算，而第二个使用三个内核。后一种实现被发现是更好的，减少了组装时间从2.44秒到0.65秒的512，000元素。然而，本文没有讨论在数值积分计算中线程之间的工作分配以及用于存储中间值的存储器类型。 在Bana's 、Paszew 'ski和Maciovich（2014）中，研究了可以扩展到高阶元素的数值积分策略。该方法使用一组线程（称为工作组），通过循环遍历积分点来完成单元刚度矩阵块的计算。在所有测试用例中，GPU策略的性能都优于CPU。在Komatitsch et al. （ 2009），基于GPU的地震模拟波传播是使用高阶谱元件完成的。为每个由125个节点组成的元素分配一个线程块，以计算元素贡献。用着色方法处理全局装配中的竞争条件问题Cecka等人（2011年）提出了另一种更通用FEM实现的全局装配详细分析。作者研究了使用不同GPU内存类型存储元素数据和不同内核设计进行汇编的许多方法一般结论建议使用两种不同的方法，低，高阶元素。对于低阶元素，元素数据存储在GPU的共享存储器中，而组装通过将线程与全局矩阵的非零（NZ）条目相关联来完成。对于高阶元素，最佳的方法是每个元素使用一个线程来计算元素数据并将其存储到GPU的全局内存中。通过并行归约完成方程组的组装基本子程序被视为一个黑盒子，这可能会妨碍使用高阶元素时的性能。虽然大多数文献都集中在非常具体的实现上，但在Reguly和Giles（2013）中，提出了一种可扩展到高阶有限元的FEM实现。元素矩阵计算是通过为每个元素分配一个线程来数值积分方案，它使用外循环形状函数和内循环高斯点。数值积分不使用任何存储器来存储中间数据，并且每次都通过进行冗余计算来计算它们其他的方法，使用本地存储器或全局存储器被发现是内存带宽的限制。它还提供了一个洞察的数据结构上的性能的影响除了全局矩阵组装方法，局部矩阵 方 法 （ LMA ） 和 无 矩 阵 方 法 也 在 文 献 中 讨 论 Markall 等 人（2013）已经在不同的众核架构上组合了几种不同的组装策略。作者提出了LMA方法的组装，并证明了这种方法的适用性，众核架构的二维网格。 在Fu et al. （2014），发现了一种不同的网格划分方法。将网格划分为称为面片的不相交节点集。元素分配是以这样一种方式完成的，即每个元素只属于一个面片那个...×联合Kiran等人 /Journal of Computational Design and Engineering 6（2019）705-718707通过以合并的方式读取输入数据并使用每个元件的单个线程来计算元件刚度矩阵（44）全局矩阵使用CSR存储格式，该格式具有预计算的列索引和行偏移。对于属于共享存储器中的补丁的所有元素进行组装，并且最终数据以合并的方式写入全局存储器中的全局矩阵。在Zayer，Steinberger和Seidel（2017）中发现了一种基于稀疏矩阵乘法的新有限元装配策略。该策略通过稀疏矩阵表示捕获网格的连通性信息，并使用它们来有效地组装元素矩阵，避免任何类型的预处理。另一种最近的有限元组装方法见于Sanfui和Sharma（2017）。在这项工作中，三维有限元的装配过程分为一个节点的节点符号部分和一个元素的元素的数字部分。上述研究中提出的策略在GPU上实现了显着的加速，然而，这些策略已经过时了。但这些方法仍然可行。在总结前人工作的基础上，我们发现数值积分方法不需要片外存储器，具有冗余计算的特点，更适合GPU实现。对于低阶单元，采用局部记忆或全局记忆的方法性能更好，但对于刚度矩阵尺寸较大的单元，记忆带宽受到在使用每个元素一个线程策略的情况下尤其如此。在这项工作中，我们提出了一种策略，使用多个线程做一个元素的计算这为我们提供了足够的片上存储空间，用于数值积分（无冗余计算）以及元素矩阵。现有的着色方法被用来避免在装配过程中的竞争条件。网格被划分成不相交的元素集表示不同的计算相应的每一种颜色是由GPU内核顺序进行。所提出的策略使用同一内核实现元素3. 预赛3.1. GPU架构和CUDAGPU是面向吞吐量的设备，设计用于处理数据GPU的硬件设计基于众核处理单元，与CPU等多核处理单元有很大不同众核设计包括大量的简单处理单元，非常小的高速缓存，和高内存带宽。另一方面，多核处理器由几个高度复杂和高效的核心以及大缓存和低内存带宽组成。本文采用了NVIDIA公司的GPU，它由多个流式多处理器（SM）组成每个SM包含多个流处理器（SP）和 片 上 存 储 器 。 它 还 有 一 个 被 称 为 全 局 内 存 的 片 外 存 储 器（DRAM），在GPU上具有最高的延迟。片外存储器也用于CPU和GPU之间的数据每个SM可用的片内存储器进一步分为寄存器、可配置共享存储器和只读数据缓存。寄存器被分配给一个线程，并且不能被任何其他线程访问.共享内存是线程块的公共内存，它可以被该块的所有线程访问。片上存储器比片外存储器快得多。本地内存是每个线程的私有内存空间。这种存储器通常占用寄存器中的空间，但可能溢出到片外存储器。常量内存对于GPU线程是只读的。 图 1显示NVIDIA Tesla K40流式多处理器以及内存层次结构。计算统一设备架构（CUDA）（NVIDIA，2016）由NVIDIA创建的并行编程平台。它允许软件开发人员使用支持CUDA的图形处理单元（GPU）进行通用计算。CUDA提供了使用多种语言（如C，C++，FORTRAN等）对硬件（GPU）进行编程的灵活性。感兴趣的读者可以参考（NVIDIACorporation，2016）了解更多详情。3.2. 线弹性有限元列式有限元法是求解偏微分方程的一种数值方法.控制偏微分方程，即线性弹性材料的强形式，由下式给出：$·rb¼0; 1其中r是柯西应力张量，b是单位体积的体积力强形式由下式给出的位移u和牵引力t边界条件补充：数据并行和数据吞吐量任务。数据并行应用是指在独立数据集上运行的应用，而数据吞吐任务需要处理大量的uu o在Cu上;tt Ct;ð2Þ图1.一、NVIDIA Tesla K40基于Kepler架构的SMX它有15个SMX和192个CUDA核心。该图还显示了线程的内存层次结构164Xð Þ642XuuuXzP PPðÞ×yzy4@/i=@5708U. Kiran等人/Journal of Computational Design and Engineering 6（2019）705-718式中，uo和t是在2@/i=@x0 03位移边界CU和牵引边界Ct。线弹性材料的本构方程由胡克定律给出为Bi¼60@/i=@y00 0@/i=@z0@/=@z@/=@y7ð12Þ我我6 7联系我们这里，D是材料本构张量，s是由下式给出的小应变张量：z0@/i=@x@/i=@y@/i=@x0本构矩阵D包含材料的弹性性质，并且对于各向同性材料具有以下形式：联系我们21 - mmm0003s¼2ruru：24小时m1-m m0 0 0E6m m1-m0 0 07对于任何任意的几何形状和边界条件，都很难得到（1）中给出的强形式的解因此，使用Galerkin加权残差方法（Reddy，1993）将强形式转换为弱形式。在其上的域直径1米直径1-2米直径60 0 01-2m000 0 0 0 1-2m00 0 0 0 01-2米ð13Þ所寻求的解决方案被离散化为称为元素的多个多边形或多面体。在一个典型的单元域Xe上的点x处的位移近似为：n联系我们ui/ix;51/1其中，n是每个元素的节点总数，ui是节点i的节点邻接元素向量，并且fix是节点i的形状函数。这可以重新排列为以下矩阵形式：8>ux9>=1/4年e其中E是杨氏采用基于高斯求积规则的数值积分方案来积分由（10）给出的表达式。 算法1给出了元素矩阵数值积分的基本步骤。在每个高斯点处，形状函数矩阵（f）和应变-位移矩阵（B）被一起Jacobian矩阵（J）的行列式（jJj）和逆还计算了将形状函数值从自然（参考）坐标变换到物理坐标所需的雅可比（J-1）一旦求出应变-位移矩阵B，就可以通过进行所需的矩阵乘法来生成单元刚度矩阵。乌乌>：uz>;1/4/D算法1.数值积分其中f是由下式给出的形状函数矩阵：/1个/¼ 00375;1071：将Ke、 Fe和Qe初始化为零2：对于q= 1到高斯点数，3：读取/（q）和d/（q）4：读取节点坐标5：计算J;j，j，j和J-1De是元素位移矢量，6：计算B7：Ke BT DBDe¼. u11 1 2y z x2 2.. . unn电话：+86-020-8：Fe/Tb这里，ui ui ui是9： Qe/T<$tx;y;z;第i个节点。将（6）代入弱形式，并在所有元素上进行组装，得到由下式给出的整体形式：KD F Q;9其中，K^e K^e; F^e F ^e; Q^e Q^e。元素矩阵的表达式如下：10：结束有限元计算需要两种数据集，一种是节点数据矩阵Cq，另一种是连通矩阵Pe;i。节点数据矩阵包含具有全局节点编号q的所有节点的坐标值和与该节点相关联的任何其它字段值。连接矩阵P∈e;i∈包含：每个元素e的第i个节点的全局节点号，即，它预先-Ke¼Fe¼Qe¼BT DBd x;e/Tbd x;/T<$tds：Ceð10Þ发送从元素的局部节点号到全局节点号的映射。整体矩阵是通过将网格所有元素的元素矩阵组合成一个方程组而得到的单元刚度矩阵Ke的阶数为3n×3n，单元体积力矢量Fe和单元外力矢量Qe的阶数为3n1，其中3n为每个单元的总自由度。矩阵B也称为应变-位移矩阵，定义为：B½ B1B2B3.. . ：Bn];113.3. 顺序有限元装配将各单元的单元刚度矩阵和单元载荷矢量每个元素对全局矩阵的非零值的贡献取决于网格的连通性模式。网格的连通性也影响全局矩阵 的 稀 疏 性 。 顺 序 球 组 装 的 最 常 见 方 式 被 称 为 Addto 算 法（Markall，Ham，Kelly，uuuR我75XR00/2个00... /n00/1个00/2个00.. . /n0Xe←¼ðÞ←ððÞðÞÞþ ¼ðÞ;Þ联合Kiran等人 /Journal of Computational Design and Engineering 6（2019）705-7187092010），其在算法2中呈现。所有单元的外循环使用单元子程序来计算单元刚度矩阵和载荷矢量。 一旦计算出这些值，则基于由连接矩阵Pdle; i给出的从局部DOF到全局DOF的映射来完成组装。算法2.顺序FEM装配或Addto算法1：将K和F初始化为零2：对于所有元素edo3：Ke;Fe←elementalSubroutinee4：对于所有本地节点n1，5：FPe;n1Fen1//地方部队向量6：对于所有本地节点n2，7：KP e;n1;P e;n2 Ken1;n2//局部矩阵8：结束9：结束10：结束3.4. 基于着色方法的着色是有限元装配中避免竞态的最简单、最常用的方法之一。它是健壮的，适用于各种问题和元素类型。如Cecka等人（2011）和Komatitsch等人（2009）所示，着色方法特别适合高阶元素。着色方法背后的主要思想是将网格划分为元素集合，使得属于同一集合的两个元素没有任何公共节点。在这里，每个集合都用唯一的颜色标识。现在可以同时对对应于一种颜色的元素进行组装，因为在不同元素上工作的线程将无法访问同一块内存。不同颜色的组件以顺序方式完成，从而防止任何竞争条件的可能性。算法中给出了基于着色方法的并行有限元装配3. CUDA内核按顺序为每种颜色启动，而在每个内核中，一个线程被分配给一个元素。Ek是地图-从第k个颜色集合中的局部元素号ping到全局元素号。输入数据由每个线程访问其全局元素编号。一旦单元刚度矩阵和载荷向量的计算结束，就由同一内核完成到全局矩阵的组装。算法3.基于GPU的着色法并行有限元装配1：对于所有颜色k做2：threadId←blockDim：xωthreadIdx：y块threadIdx：x3：ifthreadIdNk then//Nk：数组中元素的个数<第k色4：e Ek threadId//Ek：局部到全局元素映射5：计算Ke和Fe6：K累加（Ke）7：F累加（Fe）8：如果结束9：结束在采用逐单元装配策略时，必须确定用于存储输入和输出数据（如节点坐标、单元刚度矩阵和连接矩阵）的存储器类型。在GPU提供的片上和片外存储器之间进行选择。片外全局存储器速度慢且尺寸较大，而片上存储器速度快但尺寸较小。与输入和输出数据一起，数值积分中所需的几何参数（雅可比、行列式和雅可比的逆因此，为所有这些数据选择内存类型时，必须考虑其对性能的影响。此外，它还显示了根据内存类型的选择可以产生的不同实现的数量。在实验中，我们使用了基于Kepler架构的NVIDIA显卡 NVIDIATesla K40 GPU上可用的片上内存资源为每个SM 48 KB共享内存和每个SM 65536个寄存器。对于完全占用，共享内存由2048个线程共享（Kepler架构）。因此，对于每个线程，允许存储在共享内存中的值的数量由下式给出：48×10246 142048× 4考虑到每个值的大小等于四个字节。每个线程可用的寄存器数量应为32。这是在对比的要求，我们的测试问题，使用8节点hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hexahe-hold元素。每个元素需要空间来保存8个节点的坐标的24个值和每个节点的全局连通性的8个值。由于有八个节点，每个节点有三个自由度，因此每个元件的总自由度为24。因此，元件刚度矩阵的大小变为576。除此之外，还需要存储数值积分中的几何参数。我们可以很容易地观察到可用的快速内存资源的数量和元素需要处理的数据大小之间的差异。这表明，每个元素使用单个线程的组装策略对频繁需要的数据的快速片上存储器的使用施加了限制。我们解决这个问题，在我们提出的组装策略，通过分配一个扭曲的元素。对于我们的测试问题，我们尝试通过元素策略（每个元素使用一个线程）来优化组装。实现，我们稍后使用我们提出的装配策略进行比较，使用本地存储器来存储元素刚度矩阵。数值积分使用外循环高斯点，而两个形状函数的内部循环。每个高斯点计算一次几何参数，并在内部循环中用于计算所有刚度矩阵项。 形状函数在自然坐标系中的导数在CPU上预先计算。它存储在共享内存中。每次需要时，节点坐标都以合并方式读取。一旦单元刚度矩阵的计算结束，连接矩阵以合并的方式加载到本地存储器中。在连接矩阵和预先计算的CSR值数组索引的帮助下，以CSR格式完成对全局矩阵的组装。4. 建议的GPU组装策略我们提出的装配策略使用32个线程（经线）每个元素的计算元素刚度矩阵，×¼×←←ωþ←←ðÞ←ðÞ←ðÞ710U. Kiran等人/Journal of Computational Design and Engineering 6（2019）705-718集成到全球矩阵中。如第3.4节所讨论的，对于每个元素使用单个线程的策略，线程可以在每个大小为4字节的共享内存中具有6个值通过相同的逻辑，将一个warp分配给一个元素允许我们在共享内存中为每个元素保留32 6 192个值。这比保存节点数据和连通性矩阵所需的空间更多。我们提出的策略还使用共享内存来保持需要共享内存中的额外空间。由于基于warp shuffle的方法耗时最少，并且不需要额外的共享内存空间，因此本文选择它来寻找雅可比矩阵。算法4. 翘曲数值积分数值积分中的中间参数。更多的空间共享存储器意味着较低的寄存器压力、较高的数据重用以及很少或没有本地存储器访问。使用32个线程来计算单元刚度矩阵还需要减少每个线程上的计算强度。一个变形被分配给一个元素，用于计算元素矩阵以及组装。我们将计算平均分配给所有线程。因此，对于在单元刚度矩阵中具有576个条目的单元我们可以选择将总共576个线程中的任何18个条目分配给一个线程。但是考虑到它在实施和组装中的效果我们分别将这些策略称为“在以下各节中，介绍了这些战略的不同变体。本文中提出的组装策略提供了一个高效的GPU内核，可以与流行的着色方法一起使用4.1. 每排使用4根线通过经线组装图2中示出了用于在经线中工作的线程的18个条目的选择。每个节点具有三个自由度的8节点六面体单元的单元刚度矩阵（24 24）的条目被划分为32个线程。分配给每个线程的条目跨越三行六列。每排使用4根线的经线组装如图所示。3.第三章。节点坐标和连接矩阵以逐元素方式放置对于每个颜色集合完成输入数据的这种重新排序。这使得曲速能够读取以最有效的方式一次生成24个坐标值和8个连接值（合并）。在自然坐标系中的形状函数的导数是预先计算的CPU上的所有高斯点。一旦这些数据被加载到共享内存中，我们就通过涉及线程的许多线程来进行数值积分雅可比矩阵的计算方法是将节点坐标乘以形函数导数。雅可比矩阵计算的示意图如图所示。 四、24个线程从共享存储器读取形状函数导数的24个值，并与节点坐标的8个值CUDA的warp shuffle特性用于减少分布在八个线程上的八个值。在这里，一次计算三个雅可比它需要两个以上的计算通道，以找到雅可比矩阵的总共九个条目。然后使用图3所示的一个线程计算行列式。该决定因素是由一个传统的方法，涉及辅助因素的评价。总共九个线程用于找到雅可比矩阵的逆每个线程计算雅可比矩阵的一个条目的余因子并存储在其寄存器中。寄存器中的这些值除以共享内存中行列式的值最后，通过对雅可比矩阵进行反变换，将其存储在共享内存中。我们再次使用24个线程来找到物理坐标中形状函数的导数数值积分步骤如算法4所示。本文还研究了雅可比矩阵的其他计算表1显示了使用原子函数、共享内存和warp shuffle的方法的warp（4个线程）组装的运行时间。使用原子函数的方法比使用warp shuffle的方法大约多花34%的时间基于共享内存的方法比曲速重排多消耗5.8%的时间，1：threadIdblockDim：x threadIdx：ythreadIdx：x//块中线程的索引。2：LaneIdthreadId% 32//线程束中线程的索引。3：Derivative shapefn nat读取形状函数导数从全局内存到共享内存的自然协调4：对于所有高斯点，5：如果LaneId24，则<6：Jacobian Jacobian7：如果结束8：如果LaneId1，则<9：行列式计算行列式。10：如果结束11：如果LaneId9，则<12：雅可比矩阵 的反演算法13：如果结束14：如果LaneId24，则<15：求出物理坐标系中形状函数的导数16：如果结束17： 计算元件刚度矩阵18：结束元素矩阵的条目的计算由每个线程独立地完成。每个线程评估乘法-tion BT DB为他们的18个条目和存储在本地内存中。线程在物理上选取它们相应的形函数导数从共享存储器的坐标，并执行与材料属性矩阵D的乘法。通过适当的索引到形状函数导数数组中，我们可以在代码中没有任何分支的情况下进行此计算。每个线程现在负责将自己的18个条目组装成全局矩阵。在组装过程中，每个线程使用相同的索引访问存储在本地数组中的计算条目，这导致零重放开销。因此，值是以最有效的方式从本地存储器访问的。将其组装到具有CSR格式的全球矩阵将单元刚度矩阵的乘积累加到CSR格式的数值数组值数组中的条目的位置在CPU上预先计算我们称之为CSR指数。由于我们使用的是一个8节点的元素，每个元素在CSR索引中总共需要64个条目单元刚度矩阵的每一行需要八个CSR指标累积到整体矩阵中。正如本节开头所讨论的，四个线程负责将所有条目组装成一行。因此，它们分两次访问CSR索引。在第一遍中，访问位于偶数位置的索引，而在第二遍中，访问位于奇数位置的索引。CSR索引的存储模式如图所示。 五、对于每个元素，CSR索引基于它们在偶数和奇数位置上的位置被存储在两个如图所示，偶数位置索引由青色表示。数据存储在列优先的约定中，这样所有的线程束都可以读取连续的内存位置。每个线程通过使用连接矩阵和CSR索引将值累积到全局矩阵（CSR格式的值数组）接入联合Kiran等人 /Journal of Computational Design and Engineering 6（2019）705-718711图二. 对于每个节点具有3自由度的8节点六面体单元，每行4个线程的单元矩阵项的分布。图三. 使用4根线通过经线组装。每个线程独立计算其分配的条目并将其存储到本地内存中。因为线程束的不同线程将值累积到全局存储器中的跨越位置。4.2. 每排使用8根线通过经线组装在这种策略中，线程现在处理不同的条目集。每个线程被分配跨越元素矩阵的六行和三列的条目，如图6所示。该图以与第4.1节中讨论的类似方式表示元素矩阵项的分布。如使用图3所描述的那样读取输入数据。数值积分也遵循与算法4中给出的完全相同的策略。 由于每个线程现在计算不同的18个条目，BT DB的计算方法相应地改变。CSR索引像前面的策略一样分两次读取它现在使用不同的标准将CSR索引划分为两组。分配给一个线程的18个条目现在跨越了六行单元刚度矩阵。这些对应于2个节点的6个自由度。现在首先对第1、第3、第5和第7节点的所有自由这需要在第一次通过中对应CSR索引在第二遍中，读取第2、第4、第6和第8节点的CSR索引数据的存储和读取方式与前面的策略类似对程序集的全局内存的访问仍然没有合并，但与前一个策略相比，它更加本地化在每行使用4个线程的策略中，线程束必须一次在全局矩阵的8个不同行中累积值，而在此策略中，值仅在4个不同行中累积这减少了使用元素刚度矩阵的每行8个线程的该策略中所需的全局存储器事务的数量4.3. 使用PPNI的翘曲装配在这种装配策略中，我们开发了一个部分并行实现的数值积分（PPNI）的装配通过一个翘曲。712U. Kiran等人/Journal of Computational Design and Engineering 6（2019）705-718见图4。示出了用于雅可比矩阵计算的访问模式。线程将乘法值存储在它们的寄存器中，以用于扭曲混洗减少。表1整体运行时间的组件通过warp（4线程）策略，用于计算雅可比矩阵的三种不同方法。解决了650万个自由度的问题。雅可比计算法时间（秒）原子功能1.39共享存储器1.09曲速洗牌1.03单元刚度矩阵的计算需要在每个高斯点处对弱形式进行评估，对于每个高斯点，数值积分的评估需要计算许多几何参数（雅可比矩阵、行列式等）。之前的策略在高斯点上使用了一个循环来进行计算。在PPNI中，我们并行计算所有高斯点的几何参数，而每个高斯点的弱形式评估是在循环内完成图7中示出了通过翘曲使用PPNI的组件的示意图。中的节点坐标和形函数的导数自然坐标构成用于计算单元刚度矩阵的输入数据。形函数的导数在自然坐标系中，计算CPU上所有高斯点的。它存储在GPU的共享内存中，以便所有扭曲一个线程块可以访问。与前面的策略一样，一个扭曲被分配给一个元素。所有高斯点的雅可比矩阵同时计算并存储在共享内存中。 这需要每个元素72个条目的空间来存储对应于8个高斯点的所有雅可比矩阵。雅可比矩阵的计算遵循与使用4个线程的线程束装配中所述相同的过程，不同之处在于现在它使用线程束的所有线程。在非对称策略中，每次使用24根经线计算雅可比矩阵的三个条目。现在，它使用所有32个线程来计算雅可比矩阵的四个条目。重复该计算18次以计算所有72个条目。清单1显示了Jacobian计算的代码。雅可比矩阵计算之后是所有高斯点的雅可比矩阵的行列式和逆矩阵的计算。八个线程用于计算八个行列式。每个线程使用27个线程来计算三个雅可比矩阵的逆。线程以九的倍数进行，使得一个线程被分配以找到雅可比矩阵的一个条目的余因子。每个线程将其余因子除以相应的行列式。雅可比矩阵的逆矩阵被存储在共享存储器中，图五. CSR索引数据读取模式。联合Kiran等人 /Journal of Computational Design and Engineering 6（2019）705-718713见图6。 每行8个线程的元素条目分布。图7.第一次会议。装配采用部分并行数值积分（PPNI）策略。清单1.雅可比估计使用翘曲洗牌。该代码显示了24个雅可比矩阵的计算对应于8个高斯点，我们有8个雅可比矩阵，共有72个元素。714U. Kiran等人/Journal of Computational Design and Engineering 6（2019）705-718雅可比矩阵。该算法的其余部分，如计算形状函数在物理坐标和元素矩阵的导数保持与以前的策略相同BT DB的计算在每个高斯点上顺序进行通过计算形状函数在物理坐标系中的导数。单元刚度矩阵存储在所有经线的本地存储器中。一旦找到了单元刚度矩阵，连接矩阵就被加载到共享存储器中并完成装配。建议的PPNI组装策略是使用每行4个线程和每行8个线程分配来实现的。4.4. 两核在前面的三种策略中，每种颜色使用一个内核在该策略中，使用两个内核，使得单个复杂内核处的计算负载被分成两个内核。这两个内核经过优化，为它们分配的任务提供最佳结果第一内核计算将形状函数值及其导数从自然坐标变换到物理坐标所需的几何参数第二个核计算单元刚度矩阵并累积到全局矩阵中。这一战略如图所示。8.第八条。第一个内核从全局存储器读取输入数据（自然坐标中的节点坐标和形状函数导数）并存储到共享存储器中。雅可比矩阵、行列式和雅可比矩阵的逆矩阵的计算方法与使用PPNI的翘曲装配中所述的方法相同。一旦计算出逆，它与行列式一起存储在全局存储器中。然后在第二内核中以合并的方式读取逆。第二个内核将576个线程的线程块分配给一个元素。由于单元刚度矩阵中有576个条目，因此一个线程负责一个条目的计算和组装。同时计算所有高斯点的应变矩阵（B）并存储在共享存储器中。每个线程依次计算每个高斯点对应的单元刚度矩阵项，并将它们相加以获得最终值。与以前的策略一样，组装以非合并的方式完成CSR格式。元素刚度矩阵的条目以逐行方式分配给线程，使得前24个线程在第一行中被分配条目，而下一行中的下24个线程被分配条目。这确保了在组装期间，对于元件刚度矩阵中的连续条目，其对应的螺纹访问连续位置或具有较小刚度的位置。全球矩阵中的一大步。这在数据访问中提供了比每行使用4个线程和8个线程的策略更高级别的本地化。因此，存在更少的全局存储器事务和更高的高速缓存命中量。这比以前的策略更有效地减少了非协调性访问5. 结果和讨论5.1. 问题陈述建议的装配策略进行了比较，最常见的战略组件的三维弹性悬臂梁问题的元素。图9示出了一端具有点载荷而另一端固定的3D悬臂梁。为简单起见，网格划分采用8节点等边长砖形单元。然而，所提出的策略不使用任何几何属性的网格，并把它们作为非结构化。