支持GPU的分布式一致非结构均匀网格的可配置开源数据结构

3.1程序设计语言

数据结构是用C++语言实现的，这是一种用于现代高性能应用程序的最先进的编程语言。节点内CPU并行化基于OpenMP接口，使用CUDA框架实现对GPU使用的支持[34]MPI为分布式CPU和基于GPU的系统提供了高级并行化。我们广泛使用C++模板，而不是CUDA框架中支持有限的虚拟函数[34]. 一般来说，这种方法还可以产生更高效的代码[12,32]以增加编译时的工作量为代价。此外，它还鼓励在编译时检测编程错误。

3.2静态配置

根据所选择的编程范式，我们旨在满足高可配置性要求静态配置即在编译时提供和解析的配置。一个明显的缺点是，如果需要对配置进行运行时选择，那么它仅限于在编译时显式实例化的有限配置集。用户配置作为主类模板的模板参数提供网格这将在下一节中介绍。

配置的目的是允许用户通过仅选择必要的网格特征来针对特定目的微调网格表示。一个主要示例是消除应用程序无法访问的网格实体和关联矩阵。另一个例子是选择用于存储坐标和索引的内部数据类型，这可以进一步减少内存需求，从而提高计算效率。

附录中给出了一个完整的默认配置示例，从中可以轻松导出用户配置A类。配置允许更改以下参数：

—	单元的形状，表示所有其他网格实体的形状。它表示为单元格拓扑在附录中A类内置实体形状的示例如附录所示B类.
—	尺寸\网格嵌入空间的（N\ge D\）。注意N个表示为空格尺寸在附录中A类和D类是单元格的尺寸（表示为网格尺寸).
—	坐标数据类型（例如。，浮动,双重的)、全局和局部索引类型（例如。，整数和短整型). 类型表示为RealType（实际类型）,全局索引类型和本地索引类型在附录中A类.
—	尺寸对\（d_1，d_2），其中\（d_1\gt d_2），其中\（d_1）维实体将链接存储到\（d_2）维子实体。存储在网格中的每个实体都必须存储到其子曲面的链接，这些子曲面提供了它们的几何实现。查看函数子实体存储在附录中A类.
—	尺寸对\（d_1，d_2），其中\（d_1\lt d_2），其中\（d_1）维实体存储与相邻实体的链接\（d_2）维超实体。查看函数超级存储在附录中A类.
—	存储内部/边界标记的图元的尺寸。查看函数boundaryTags存储在附录中A类.
—	对偶图及其参数的存储\（n_｛common｝\）。查看函数双图形存储和属性双图形最小公共顶点在附录中A类.

3.3内部数据结构

一个重要的设计选择涉及用于存储非临时数据的内部数据结构。如第节所述2，可以使用可以以（二进制）稀疏矩阵格式存储的关联矩阵来表示整个网格拓扑。类似地，对偶图可以用它的邻接矩阵来表示\（A_D\），也可以以稀疏矩阵格式存储。其他数据（例如顶点坐标和图元标记）可以存储在由相应网格图元索引索引的普通数组中。与网格实体的面向对象表示相比，使用数组和合适的稀疏矩阵格式不仅降低了实现的复杂性，而且还提供了性能优势，因为可以有效缓存对数据元素的随机访问，并且可以使用合并内存访问模式[34]可以在GPU上使用。另一方面，在现有网格表示中插入新实体等操作的实现将非常复杂，并且由于数组和稀疏矩阵的重新分配，可能会导致性能下降。因此，不支持此类操作。

请注意，网格实体本身不必显式存储在计算机内存中。相反，他们可以实例化基于实体维度和索引的按需访问，因为所有拓扑信息以及顶点坐标等附加数据都可以从全局网格对象的内部数据结构中获取。在我们的实施中，网格实体是的泛型类模板d日-维度网格实体，它仅包含指向全局网格的指针和实体索引，而实体维度d日在类模板参数中编码。当通过网格指针访问全局数据结构时，各种成员函数提供实体的维度和索引。

表示网格实体之间连接的关联矩阵通常是稀疏的，因此我们考虑稀疏矩阵格式来组织内存中的数据。由于关联矩阵是一个二进制矩阵，因此通过省略矩阵元素值的存储并仅存储非零元素的列索引，对传统格式进行了修改。

许多稀疏矩阵存储格式的综合评估最近在[27]. 我们考虑了网格中使用的关联矩阵的以下常用格式：

(1)	这个压缩稀疏行（CSR）格式[9,40]易于实现，并且由于填充而没有任何开销。同一行中两个连续元素之间的步长始终为1。
(2)	这个Ellpack公司格式[24,40]在比最长行短的所有行的末尾添加填充零，以获得二维数组，该数组可以存储在row-mahor或column-major方向。因此，同一行中两个连续元素之间的步长为1或米，其中米是矩阵行数。如果有几个很长的行，填充导致的开销可能相当大。可以根据硬件体系结构选择行主方向或列主方向。例如，由于合并访问模式的较高内存吞吐量，列主方向对GPU有利[34]与赛艇主方位所暗示的跨步模式相比。
(3)	这个Ellpack切片格式[31,37]是Ellpack格式的块修改。存储阵列被拆分为多个块，其中包含以列主（或行主）方向存储的恒定行数。这在行和列中的连续元素之间提供了编译时常数步长。主要变体仍然允许合并访问，假设块大小是当前硬件上32的CUDA扭曲大小的倍数[34]. 此外，没有必要添加填充项来与全局最长的行对齐，但可以考虑每个块中最长的一行来减少内存开销。

我们实现的最重要标准是存储开销和GPU的联合访问模式。因此，我们使用基于Sliced Ellpack格式的稀疏矩阵，其中行主方向用于CPU计算，列主方向用于GPU计算。

上述稀疏矩阵格式没有显式存储每行非零元素的数量，因为它对于矩阵-向量乘法等常见算法来说是不必要的。然而，网格算法可能需要直接了解给定实体的子实体或超级实体的数量。在关联矩阵的情况下\（I_{d，d_1}\），\（d\gt d_1）表示实体及其子实体之间的关系，每行非零元素的数量（即\每个的（d_1）维子实体d日-维度实体）是先验的常量，因此可以将其作为编译时常量提供。另一方面\每个的（d_2）维超实体d日-维度实体不是常量，因此关联矩阵\（I_{d，d_2}\），\以稀疏矩阵格式存储的（d\ltd_2）必须由存储每行非零元素数的数组明确补充。邻接矩阵也是如此\对偶图的（A_D\）。

3.4内部存储层

数据结构需要包含用于网格表示的几个关联矩阵和数组，具体取决于配置。对于网格中实体的每个维度，内部数据结构可以组织在多个层中。可以使用C++功能组合各个存储层，例如递归继承和部分类模板专门化为了提供高效通用的D类-维度数据结构。递归继承是一种C++设计模式，其中类模板的实例继承自同一类模板的另一实例，但具有不同的模板参数。因此，在这种情况下，继承不是用来表达对象之间的关系，而是一般地将任何有限数量的不同类型的属性包含到单个对象中。

的继承图网格类模板如图所示4.模板参数配置，表示第节中描述的网格配置3.2，以及设备，可以是TNL:：设备：：主机或TNL:：设备：：Cuda并代表装置网格将在何处分配，将传播到层次结构中的所有基类模板。第一个基类模板是StorageLayer系列它没有任何数据成员，其目的是通过标签调度向所有单个存储层提供访问。递归继承在存储层由实体维度标记参数化的类模板，DimTag标签，并包含与d日-多维网格实体。这个StorageLayer系列类模板从维度0开始继承，每个存储层从同一个类模板继承，但带有递增的维度标记。最后，递归由空\（（D+1）\）维部分类模板专门化。

图4。

查看地物

通用的继承图d日-维度的存储层类模板如图所示5（实例之间的继承存储层使用参数DimTag标签和DimTag<d+1>如图所示4). 继承方案类似于图4，现在使用两个独立的存储层族来存储与子实体和超实体相关的数据d日-维度实体。我们介绍子实体存储层和超级实体存储层由实体维度参数化的类模板d日和\（d^{\prime}\）并分别存储一个特定的关联矩阵：子实体存储层商店\（I_{d，d^{\prime}}\），其中\（d\gt d^{\prime}\），以及超级存储层商店\（I_{d，d^{\prime}}\），其中\（d\ltd^{prime}\），以及包含每行非零元素数的附加数组（请参阅第节3.3). 请注意，的继承子实体存储层由启动\（d^{\prime}=0\）并继续递增\（d^{\prime}\），但对于超级存储层它是由\（d^{\prime}=d\）并继续递减\（d^{prime}\）。在这两种情况下，当实体维度与子实体或超级实体维度相同时（即。，\（d=d^{prime}\））。

图5。

查看地物

在下面的代码示例中，我们概述了图中方案的实现4使用递归继承和部分类模板专门化。递归继承可以简化为以下形式：

查看地物

假设DimTag标签参数最初提供了一个0维标记，该标记增加一个以继承下一个存储层。注意，由于C++语言的限制，我们必须使用DimTag标签类型模板参数而不是非类型模板参数，例如整数d表示维度，因为涉及其他类型模板参数的表达式不能用于部分专用化非类型模板参数。在使用了给定网格的所有存储层之后，将提供空的部分类模板专门化，以终止递归继承。空的专门化由以下内容标识配置：：meshDimension+1，其中配置：：meshDimension表示网格尺寸D类:

查看地物

在实践中，完全继承比图中所示的方案更复杂4和5，由于网格类模板有附加的基类，为了清楚起见，本文省略了这些基类。除了关联矩阵之外，各个存储层还存储更多数据。例如，有数组存储标签（例如边界实体,GhostEntity（鬼实体）或任意用户属性）。另一方面，一些不必通过实体维度参数化的数据成员（例如，顶点坐标数组）直接放置在网格类模板。最后，请注意，存储层层次结构中的所有数据成员始终存在，即使它们根据配置未使用。在配置中禁用存储只会跳过动态分配和统计_断言声明保护对禁用数据成员的访问。

3.5将用户数据映射到网格图元

网格数据结构主要用于计算算法，如有限元或有限体积方法。网格计算的一个重要部分是定义用户数据和网格实体之间的映射。与OpenMesh等其他库不同[8]在mesh对象中可以直接存储任意用户数据的地方，本文中介绍的数据结构并没有解决将用户数据映射到mesh的问题。网格实体可以通过其C++类型（即拓扑）和全局索引进行标识，并且可以使用局部索引访问相邻实体。此信息可用于定义网格实体与适当的数据容器之间的自然映射，该数据容器独立于网格存储用户数据。根据具体应用，合适的容器可能是（一维）数组、多维数组、结构数组或数组结构。TNL提供一维和通用多维数组的实现(阵列和NDArray公司)包括全面的GPU支持。后者的存储布局可通过索引排列进行配置，这允许轻松切换内存布局（例如，如果是二维数组，则为row-mahor和column-major）。

3.6网格图元的方向

对于有限元或有限体积法等计算算法，通常需要定义方向网格实体的数量。不幸的是，这个词方向通常用于两种不同的上下文：

(1)

在有限体积法的情况下，人们对唯一分配给网格每个面的法向量的方向感兴趣。每个面最多可以有两个相邻单元，其法向量可以两种方向：相对于第一个相邻单元向外或向内（即相对于另一个单元向内或向外）。这些信息对于一些低阶有限元离散化也足够了，例如，我们在中使用的最低阶Raviart–Thomas元素[15]. 基于此上下文的面方向隐式存储在显示的数据结构中。虽然它没有显式存储每个面的法向量，但可以根据面及其相邻单元的拓扑推断其方向。核心实现细节是，关联矩阵每行中的列索引没有排序。因此，对于每个网格实体，我们可以定义其子实体和超级实体的局部索引，这些索引对应于与该实体相关的行中的位置。例如，关联矩阵\（I_{D-1，D}\）每行包含一个或两个条目，对应于与特定行索引给定的面相邻的单元格。第一个条目被分配一个局部索引0，其列索引对应于所谓的所有者单元格而第二个条目（如果存在）被分配一个局部索引1，其列索引与所谓的相邻小区面部表情。面方向可以从其所有者单元的拓扑推断出来（参见附录B类)，因为所有者单元格是初始化人脸的单元格。 在适当的情况下，需要使用诸如以下的适当算法从数据结构中手动提取人脸的方向。例如，给定一个用索引表示的单元格K（K）其相邻面用索引表示F类，我们感兴趣的是K（K）是脸的主人或邻居F类。首先，我们可以使用模板化函数获取超级实体计数获取与面相邻的单元格数F类.如果我们有一个，那一定是牢房K（K）这也是所有者单元格。否则，我们可以将局部索引0传递给模板化函数获取超级索引获取全局索引\可与索引进行比较的脸所有者单元格的（K_0\）K（K）.如果K（K）等于\（K_0\），然后K（K）是所有者单元格，否则是相邻单元格和的方向F类需要根据\（K_0\）。

(2)

对于高阶有限元方法，通常需要在每个单元上定义局部坐标系，以便从每个相邻单元查看面或边时，可以一致地定位与相邻单元共享的面和边相关联的自由度。在这种情况下，实体的方向与它们的子类相对于相邻单元的局部顺序有关。这种情况比前者更为普遍，因为考虑了面以外的实体（例如，3D中的边），甚至面上法线向量的方向也可以基于面子曲面的一致局部排序。 然而，当前所示数据结构的实现并没有提供根据此上下文处理一般网格实体方向的任何方法。因此，实现需要此信息的计算方法（例如高阶有限元方法）可能需要使用额外的数据结构来存储必要的网格实体方向，或者采用不同的方法来避免对额外存储的需要。有关处理这一一般情况的各种策略，请参阅[1].

4.1网格初始化

每个网格对象，即网格<配置，设备>类型，必须在首次使用之前进行初始化。前一个类模板参数，配置，表示特定的配置类（参见第节3.2)和后一个类模板参数，设备，表示装置其中网格将被分配。TNL目前支持两种设备类型：TNL:：设备：：主机用于CPU计算和系统内存分配，以及TNL:：设备：：Cuda用于GPU计算和分配。在我们的实现中，初始化过程是连续的，完全在CPU上完成，然后，可以使用重载赋值运算符将整个网格传输到GPU。初始化是通过公共成员函数完成的初始化的网格类模板，仅当TNL:：设备：：主机被指定为类模板参数设备.

成员函数初始化采用两个参数：点和细胞种子。前一个参数是一个数组点类型对象，后一个参数是一个数组细胞种子类型物体。这个我-中的第个对象点数组表示的笛卡尔坐标我欧氏空间中的第个顶点。这个细胞种子类型类型是泛型的实例种子类模板，它是仅在初始化期间使用的辅助类模板。它旨在根据特定顺序的顶点列表唯一地标识网格中的任意非顶点实体（参见附录B类). 因此j个-中的第个对象细胞种子数组表示定义j个第th个细胞。每个单元的其他子元素是从单元种子中顶点索引的特定顺序推导出来的。

初始化D类-在以下阶段进行三维网格划分。首先，从点然后释放的数组。接下来，通过从单元种子复制其subvertex索引来初始化网格单元细胞种子数组已释放。那么，对于所有人来说\（d=d-1，d-2，\ldots，1）d日-通过以下步骤根据网格配置初始化维度实体：

(1)	为创建中间种子d日-维度实体并将其存储在中间索引集容器表示为标准：：unordered_map<种子，索引>在我们的实施中。请注意，要确定两个种子是否表示同一实体，必须首先对subvertex索引进行排序，但通过散列可以避免大多数比较。为了避免对subversitices进行排序，subversitex索引的散列必须与交换操作（例如整数加法）相结合。
(2)	初始化的subditesd日-维度实体（即关联矩阵\（I_{d，0}）以及相反的关系（即关联矩阵\（I_{0，d}\））。
(3)	对于所有人\（s=d+1，d+2，\ldots，d\），初始化链接秒-标注实体到其d日-维子实体（即关联矩阵\（I_{s，d}\））。相反方向的链接（来自d日-标注实体到其秒-维度超实体）记录在过程和关联矩阵中\（I_{d，s}\）也被初始化。
(4)	取消分配所有中间数据结构。

最后，初始化补充数据成员，如对偶图或边界实体标记。请注意，如果网格配置中禁用了一些矩阵，则会省略上述步骤，以减少初始化的内存和计算成本。

初始化过程会导致系统错误或填充的网格对象，其中包含配置指定的所有实体及其之间的链接。网格表示的正确性取决于输入数据，在某种意义上，网格实体的拓扑属性不由初始值设定项本身检查。为了方便起见，一个名为网格生成器还提供了基本的输入数据验证和常见格式的文本文件读取器。

4.2从文件输入

实际上，科学求解器通常通过读取另一个专用工具准备的输入文件来初始化网格。TNL库目前为传统VTK、基于XML的VTK和Netgen文件格式的输入和输出提供了便利。分布式网格也支持基于XML的并行格式PVTU，这将在第节中介绍5.

将网格从文件导入到网格对象分为三个阶段。首先，通过检查文件确定网格的特征。接下来网格类模板在运行时解析。最后，将网格表示从文件导入到适当类型的网格对象中。

可以使用名为网格类型解析器来自TNL库。网格类型解析器在运行时为默认网格配置类模板选择模板参数，类似于[28]. 实施网格类型解析器相当复杂，但其用法可以用以下伪代码概括：

查看地物

用户提供的函子对象通常是C++lambda表达式它使用完全初始化的网格继续计算。

4.3实体的重新排序

第节中描述的稀疏矩阵格式3.3提供对网格实体的子实体和超实体索引的有效访问，但效率不会传播到使用获得的索引对其他实体的访问。例如，如果我们访问关联矩阵\（I_{2,1}）获取与二维网格实体相关的边的索引，然后使用获得的索引查询关联矩阵\（I{1,0}）要获得边缘子空间的索引，不能保证访问后一个关联矩阵的效率，这取决于获得的边缘索引的局部性。对于类似场景，一种常见的优化策略是通过应用合适的重排序算法来改善数据局部性[43]. 常见方法通常包括顺序遍历考虑网格实体的几何局部性的四叉树或八叉树，或图排序，如Cuthill–McKee算法[13]它可以减少图的邻接矩阵的带宽，或图的着色技术，例如[18].

这个网格类模板提供了一个接口，用于应用网格实体的任何自定义顺序。这个d日-通过调用模板化成员函数，可以对网格的维度实体进行重新排序重新排序实体并指定d日作为模板参数：

查看地物

这个全局索引矢量type代表的向量全球指数在与网格相同的设备上分配的元素。向量烫发和吐温表示从0到\（N_d-1\），其中\（N_d\）是d日-网格中的标注实体。重新排序是根据METIS公司图书馆[23]. 让\（\mathcal{M}\）是原始网格\（\mathcal{M}^{prime}\）置换网格。然后，具有索引的实体我在里面\（\mathcal{M}^{prime}\）是具有索引的实体烫发[\（\，i\，\）]在里面\（\mathcal{M}\）和带有索引的实体j个在里面\（\mathcal{M}\）是具有索引的实体吐温[\（\，j\，\）]在里面\（\mathcal{M}^{prime}\）。

通过将排列应用于所有相关关联矩阵和数组来执行重新排序。第一，烫发用于排列所有关联矩阵的行\（I_{d，d^{\prime}}\），其中\（d^{prime}\ne d\）代表的子实体或超实体的维度d日-维度实体。接下来，吐温用于排列所有关联矩阵的列\（I_{d^{prime}，d}\），其中\（d^{prime}\ned\）具有相同的含义。两者都有烫发和吐温也用同样的方法排列对偶图\（I_{d，d}\）。最后，烫发用于排列存储在普通数组中的数据（例如，如果\（d=0\））和其他属性（如边界标记）被重新初始化。

5.1分布式网格生成

在非结构化网格上进行分布式计算的明显要求是生成分布式网格，这通常是通过按照区域分解方法将现有网格划分为几个子域来执行的。分解网格时，可以将每个子域分配给不同的MPI秩，以提供并行计算。

TNL库提供了工具tnl-分解网格它使用METIS公司图书馆[23]. 类似于分解ParOpenFOAM框架中的工具[19],tnl-分解网格在主计算之前作为单独的进程顺序运行。它从VTK、VTU或Netgen格式的文件中读取输入网格表示，并将分解的网格写入基于XML的并行格式PVTU中，该格式每个子域包含多个文件。

除了网格划分之外，tnl-分解网格处理我们的分布式网格数据结构所需的几个特性。首先，它为顶点和单元指定全局索引，以便主人甚至对于子域之间接口上的公共实体，也可以轻松识别每个实体的。全局索引在每个维度中分别进行分配，分配方式为每两个等级我和j个,\（i \lt j），全部d日-等级拥有的维度实体我全球指数比任何一个都小d日-等级拥有的维度实体j个。接下来，它可以生成重叠（也称为重影区域)跨子域。这将导致从其他子域向局部网格（即表示给定子域的空间离散化的网格）添加靠近子域边界的几层顶点和单元。这对于需要交换与网格实体关联的数据的MPI列组之间的通信模式来说是必要的。最后，生成每个子域上网格实体的适当重新编号。在每个维度中，对网格实体进行排序，以便首先对当前MPI列所拥有的实体进行编号，然后对第一个邻居所拥有的重影实体进行编号以及对第二个邻居所具有的重影图元进行编号，依此类推。这种排序简化了实现并提高了虚区域上数据同步的效率，这将在后面进行描述。

5.2分布式网格实现

由于存储在PVTU文件中的分布式网格包含与原始网格相比的附加数据（例如，全局实体索引），因此需要提供一个附加数据结构，用于加载PVTU文档。TNL库提供了DistributedMesh<类型名Mesh>包含以下部分的类模板：

—	类型的局部网格网格表示特定子域。
—	重影区域级别的数量，存储为整数(整数).
—	局部网格中每个实体的全局索引数组。如果是顶点和单元，则从PVTU文件中读取这些数组。对于其他维度，如果需要，必须生成它们。
—	顶点和单元类型数组(vtk重影类型)以VTK文件格式输出数据所需的。

本地网格包含MPI列组所拥有的所有实体，该列组被分配了相应的子域，以及其他列组所具有的重影实体。可以通过查询实体标记数组来获得每个实体的类型（参见第节3.4). 具体而言，每个实体要么地方的或鬼，从几何角度来看，每个实体都是内部的或边界.

5.3Ghost区域上的数据同步

遵循域分解方法的计算通常涉及同步，以确保虚区域上的数据一致性。这是通过将由其子域上的MPI秩计算的数据发送到其邻居的重影区域来实现的，并从另一侧将邻居计算的数据接收到子域周围的重影区中。

数据同步过程在分布式网格同步器类模板作为单独的数据结构，可以为特定的分布式网格和实体维度实例化d日（即被同步实体的维度）。同步器可用于同步任何位置分配的数据阵列，例如，支持主机到主机、主机到GPU以及GPU到GPU通信。注意，由于所谓的CUDA-aware MPI[26]，在接收或发送GPU上分配的数据时，我们不必处理自定义缓冲区。

同步器需要在使用前进行初始化。初始化阶段包括以下步骤：

(1)	不对称\（P\乘以P\），其中P（P）表示MPI列数，通信矩阵G公司是这样创建的：\（G_{i，j}\）表示d日-列上的多维虚实体我按等级拥有的j个.
(2)	每个排名收集每个邻居拥有的第一个实体的索引，并将其存储在大小数组中P（P）.
(3)	每个等级收集本地实体的索引，这些实体是不同子域上的重影。这可以用二进制稀疏矩阵表示P（P）排。

在同步期间，每对列组根据通信矩阵交换数据。要发送的数据可能是非连续的，必须在传输之前将其复制到相邻的缓冲区中。另一方面，由于网格实体的适当编号由tnl-分解网格，因此在这种情况下不需要缓冲。

上述初始化过程足以同步与已经在PVTU文件中分配了其全局索引的顶点或单元相关联的数据。另一方面，输入文件中不存在诸如边或面之类的剩余实体，因此它们还没有指定全局索引。在这种情况下，必须在运行时通过MPI初始化其余实体，然后才能将同步器用于这些实体。该算法从每个MPI等级开始，确定并计算所有地方的未初始化的实体，与其他MPI列组交换计数，并为所有本地实体分配全局索引。在第一阶段无法确定全局索引的所有剩余实体都位于子域之间的接口上或虚区域中，并且它们的全局索引必须从拥有该实体的列组接收。可以观察到，未初始化的实体可以通过其subdite明确标识，subdites已经分配了全局索引。因此，相邻列交换全局subvertex索引的元组，在其本地网格中识别所请求的实体，并将全局索引发送回相邻列。最后，每个秩基于全局索引对其局部网格上新初始化的实体进行重新排序。在附加初始化之后，数据同步本身通过分布式网格同步器与单元和顶点的方法相同。

6.1计算单元尺寸、曲线长度和表面积

我们使用两个基准问题来评估数据结构在不同配置下的效率。第一个基准是计算网格中每个单元格的度量值，这涉及启动网格中单元格数量相同的独立并行任务，计算每个任务一个单元格的度量，并将结果写入全局数组。单元度量的计算涉及读取单元的子单元，并在单元度量的公式中使用它们的空间坐标。例如，三角网格的计算内核可以使用扩展lambda表达式[34]如下：

查看地物

请注意，lambda表达式内核_度量传递给TNL平行于在指定设备上并行执行它的函数。

第二个基准是分别计算围绕顶点及其相邻三角形或四面体的最小曲线或曲面的长度或面积。与第一个基准不同，每个面都执行独立的并行任务，计算包括计算面部测量值、读取面部相邻单元格、查找与面部相对的顶点并将面部测量值添加到输出数组中。由于可能有多个并行任务同时将其值写入输出数组的同一元素，因此必须自动执行加法。简单网格的计算内核可以实现如下：

查看地物

lambda表达式内核_曲面再次被移交给TNL平行于在指定设备上并行执行它的函数。

功能获取实体度量值第6行中使用的是计算给定网格中给定实体的度量，类似于前一个内核中三角形的显式计算。辅助功能获取本地面索引第11行中使用检查给定单元格的哪一面具有给定的全局索引，并返回其相对于单元格的局部索引。由于简单网格中相对的顶点和面具有相同的局部索引（参见附录B类)，变量v（v）标识相对于给定单元格与面相对的顶点。

要比较网格在TNL中使用另一个库实现了数据结构，我们还使用MOAB实现了这两个基准测试[44]. 然而，存在一些不可避免的差异：

—	MOAB中的所有网格都是隐式三维的，因此除了检查顶点坐标的第三个分量外，无法检查三角形是否来自二维或三维空间。可以使用\（2乘2）行列式，但在3D中计算三角形面积的公式涉及昂贵的平方根函数。为了实现涉及二维三角形网格的第一个基准，我们使用了基于\使用TNL网格匹配基准的（2乘2）行列式。但是生成的函数不是通用的，不能在实践中使用，因为它给出了一般三维三角形的错误结果。
—	第二个基准无法使用MOAB等效地实现，因为它依赖于单纯网格中相对顶点和面的局部索引约定。此约定在MOAB中不可用，MOAB返回给定网格实体的子实体，这些子实体按其全局索引排序，因此无法通过局部索引进行寻址。因此，得到的MOAB算法的复杂度稍高，因为我们需要得到给定单纯形的所有子曲面，并找到与给定面相对的顶点。

本小节中描述的所有基准的代码在以下项目中作为开源软件提供给社区：

—	https://mmg-gitlab.fjfi.cvut.cz/gitlab/tnl/tnl-benchmark-mesh（TNL实施），
—	https://mmg-gitlab.fjfi.cvut.cz/gitlab/tnl/moab-bequinchmark（MOAB实现）。

第一个基准测试问题的结果如表所示4第二个基准问题的结果如表所示5在所有情况下，计算时间\（CT\；[\text{ms}]\）是每次计算100次迭代的平均值。有效带宽\（EBW\；[\text{GB/s}]\）通过将有用数据大小对于给定的计算，按实现的计算时间。在第一个基准中，数据大小是由所有顶点坐标的大小加上每个单元格的所有subvertex索引的大小得出的，即。，\（\begin{equation*}\text{data-size}_1=\text{sizeof}（\texttt{Real}）\cdot D\cdot\#{\mathcal{V}}+\text{sizeof}（\textt{GlobalIndex}）\ cdot（D+1）\cdot\#{\mathcal{C}}，\end{equation*}\）其中函数\（\text{sizeof}\）返回给定类型的大小（以字节为单位），D类表示网格尺寸和实体数\（\#\mathcal{V}\），\（#\mathcal{C}）如表所示2基于第二个基准点的算法，相应的数据大小由所有顶点坐标的大小加上每个面的所有subvertex索引的大小加2个计数和每个面的相邻单元的全局索引加上每个面相对顶点的2个全局索引得出，即。，\（开始{align*}\text{data-size}_2&=（\text{sizeof}{\texttt{Real}）}\cdot D\cdot\#\mathcal{V}+\text{sizeof}（{\textt{GlobalIndex}；\text{sizeof}（\texttt{LocalIndex}）\right）\cdot\#\mathcal{F}，\end{align*}\）其中符号\（\text{sizeof}\）和D类具有相同的含义和实体计数\（\#\mathcal{V}\），\（\#\mathcal{F}\）可以在表中找到2。请注意，计算的带宽是有效的，因为我们在数据大小计算中假设所有数据只读取一次，但实际上，像顶点坐标这样的实体会从周围的实体多次访问，这些访问是否能有效缓存取决于其他因素。对表中所示的所有网格进行基准测试2，但为了简洁起见，仅显示2D和3D（2D）中最精细网格的结果\（^\三角形_5\）和3D\（三角形_5）分别如表所示4和5.

这两个基准测试都是使用不同配置的网格在TNL中为模板参数使用不同数据类型的类：浮动(（f）)以及双重的(d日)的真实的,整数(我)以及长整型(我)的全球指数，以及短整型(秒)以及整数(我)的本地索引为了简单起见，表中只显示了四组具有代表性的类型。对于MOAB，我们只使用了一组类型(双重的对于顶点坐标，长整型对于全球指数和整数相邻实体的计数），这是在安装MOAB时指定的。可以在表中看到4和5为真实的,全球指数，以及本地索引与较大类型的模板相比，模板参数通常会导致CPU和GPU上的计算速度更快。最后，我们想指出的是，使用TNL库执行的所有基准测试的性能并不依赖于其他网格配置参数，例如子实体存储或超级实体存储，可用于禁用不必要关联矩阵的存储。表中未包含此行为4和5，因为这将导致具有相同值的重复表。因此，这些配置参数会影响网格表示所需的动态尺寸，如表所示三而不是使用数据结构进行计算的性能。

基准测试中包含的操作涉及许多独立的任务，因此可以轻松地进行并行化。在这个基准测试中，CPU计算只考虑OpenMP并行化。对于表中所示的电池测量计算4，CPU并行化的效率取决于计算操作和内存访问之间的比率。该基准的2D情况涉及更简单的计算（a\计算三角形面积的（2乘以2）行列式，1.3 ops/mem）比3D情况（a\（3乘以3）行列式计算四面体体积，2.4 ops/mem）。因此，有效带宽电子束焊接用于在2D上使用12个CPU线程进行计算\（三角形5）网格接近硬件峰值性能（见表1)和加速\（Sp_{12}\）远低于12，但对于3D\（^\三角形5\）加速\（Sp{12}\）更高，有效带宽电子束焊接要低得多。另一方面，GPU并行化的效率主要取决于正在处理的数据元素的数量。最精细的二维三角形网格\（^\三角形_5\）包含大约\（60\cdot 10^3）个单元，而最精细的四面体网格3D\（^\三角形_5\）包含大约\（2\cdot 10^6\）单元格（见表2). 因此，该基准无法充分利用GPU在三角形网格上的潜力，在三角形网格中，GPU仅实现了与CPU和速度相比略高的有效带宽\（GSp{12}\）仅略大于1。另一方面，四面体网格允许基准充分利用GPU，从而使有效带宽比2D情况高出约4倍。观察GPU加速\（GSp_1）和\（GSp_{12}\）对于数据类型最小的网格配置，GPU计算大于100\（\倍\）比顺序CPU计算快，几乎10倍\比使用12个CPU线程的计算速度快（\倍\）。

由于算法差异，第二个基准与第一个基准具有不同的限制因素。根据空间维度，基准计算2D线段长度或3D三角形区域。在这两种情况下，都有必要评估平方英尺每个网格面运行一次。因此，我们希望这个基准是计算绑定的，而不是内存绑定的。实际上，表中CPU计算的有效带宽5明显低于硬件峰值性能。另一方面，GPU计算的有效带宽仍然很高，可以与第一个基准中实现的最大有效带宽相媲美\（^\三角形_5\），GPU加速\（GSp_1）上升到100以上\（GSp{12}）高于10表明，即使在网格上具有相对较少的单元，GPU也可以得到高度利用。对于四面体网格3D\（^\三角形_5\），GPU加速\（GSp_1）和\（GSp{12}）分别高于300和30。

最后，我们想评论一下TNL和MOAB中实现的数据结构之间的性能差异[44]. 对于表中所示的第一个基准4，使用MOAB的计算大约为15\2D速度较慢，约为6倍\3D中速度较慢。对于表中所示的第二个基准5，使用MOAB的计算大约为25\2D速度较慢，超过50次\3D中速度较慢。这令人惊讶，因为第二个基准应该是计算绑定而不是内存绑定，并且是最昂贵的计算，即平方英尺函数，无论使用哪种数据结构都需要完成。我们的研究表明，大约80-90%的计算时间用于MOAB函数获取相邻项（_A）它负责获取与给定网格实体相邻的实体。在第一个基准测试中，我们只使用了函数获取连接性，它在MOAB中有效地实现，以获得给定网格单元的颠覆。在第二个基准中，我们需要访问顶点以外的相邻实体，这只能通过函数实现获取相邻项（_A）这在MOAB中没有有效实现。如前所述，当使用MOAB数据结构时，第二个基准的算法不能依赖于单纯形上相对顶点和面的局部索引约定，但仅算法复杂性的增加并不能证明表中观察到的低效率是合理的5另一方面，MOAB被设计为一种动态数据结构，并提供了许多高级操作，如非结构化网格的自适应细化，这在TNL中当前实现的数据结构中是不可能的。

6.2广义McWhorter–Sunada问题

为了证明该数据结构适用于高级数值格式，我们实现了一个基于MHFEM的方案[15]求解多孔介质中的两相流。在本节中，我们给出了广义McWhorter–Sunada问题的结果，这是多孔介质中两相流的一个特例，具有已知的半解析解[16,30]已用于[15]验证数值格式。在本文中，计算域、边界条件和物理参数的选择与[15]，以及表中列出的网格2与原始文章中用于CPU计算的三角形和四面体网格相同。

与以前的基准问题不同，广义McWhorter–Sunada问题的计算时间取决于稀疏线性系统的分辨率，而不是涉及非结构化网格的操作。利用Jacobi预条件，通过BiCGstab求解器求解稀疏线性系统。解算器可以在CPU或GPU上执行。CPU执行可以使用OpenMP（共享内存，基于线程）或MPI（分布式内存，基于进程）进行并行化。GPU执行使用CUDA[34]用于在单个GPU上并行化，并且可以与MPI结合使用[26]用于多GPU执行。此基准测试中只包括使用MPI和/或CUDA进行并行化的结果。

计算时间计算机断层扫描，并行效率\（E！f！f）和GPU加速\此基准测试的（GSp_\ell）如表所示6和7GPU加速\（GSp_\ell）随着网格的细化而增加，在最精细的网格上，GPU计算速度比连续计算快20倍（2D）或50倍（3D），比使用12个内核的CPU计算快2倍（二维）或9倍（三维）。当表中最精细的2D网格上的CPU执行从6核变为12核时6，超线性加速表现为效率的提高\可以注意到（E \！f \！f）。这种现象可能是由缓存大小增加引起的，因为每个内核都有自己的大小为1 MiB的二级缓存（参见表1)其他内核无法使用的。如果内存需求与缓存大小相当，例如网格2D上提到的计算出现问题\（^\三角形5\），增加核心数可能会比核心数成比例地提高性能，因为缓存中可能有更多的数据可用，并且可以避免对系统内存的访问。

计算时间的比较计算机断层扫描和加速\涉及多GPU计算的基准测试的（Sp_\ell）如表所示8和9。请注意，每个MPI列组管理其专用GPU，因此每次计算使用的GPU数量与MPI列的数量相同。可以在表中观察到8使用多个GPU不会导致更快的计算（加速\（Sp_\ell\）小于1）。这可以归因于多个GPU之间的通信延迟，这与实际计算相比是很高的，因为基准中使用的三角形网格没有提供足够的自由度来充分利用GPU的计算能力。另一方面，基准中使用的四面体网格提供了更多的自由度，因此加快了速度\表中的（Sp_\ell\）9高于1。在最精细的四面体网格上，4阶的加速比大约为3.3，这与理想情况相比相当小(\（Sp_4=4\））。这是由于基准测试中使用的强缩放方法造成的：每次计算都使用不同数量的MPI列在同一网格上重复，随后，随着MPI列数的增加，每个GPU的利用率越来越低。通过优化多GPU计算的线性系统解算器（BiCGstab）可以提高速度，但这超出了本文的范围。

表中比较了使用MPI分布在多个节点上的CPU计算10。此强缩放性能研究仅在最精细的四面体网格3D上进行\（^\三角形_5\）。用于计算的计算集群由20个双处理器节点组成，其中包含表中列出的CPU1，但一次MPI计算最多只能使用16个节点。加速服务提供商使用12个岩芯进行计算，包括在表中7已经。可以注意到，加速通常高于CPU的数量，这可能再次归因于缓存总大小的增加以及在强规模研究中问题大小保持不变。比较网格3D上的计算时间\表中的（^\三角形_5\）9带表10可以看出，与使用8个CPU（4个节点）相比，使用1个GPU的计算时间更快，与使用4个GPU相比，至少需要32个CPU（16个节点）才能获得更快的计算时间。

本节中显示的最终结果是表中总计算时间的细分10在最精细的四面体网格3D上\（^\三角形_5\）。在表中11可以观察到，主要部分对应于稀疏线性系统解算器（BiCGstab），它在每次迭代中涉及两个稀疏矩阵-向量乘法、几个点积和其他BLAS-1操作[40]. 前两种类型是集体操作，在使用MPI时会引入通信延迟。在每个稀疏矩阵向量相乘之前，输入向量的虚区域需要与同步分布式网格同步器表的第四行中包括了对总计算时间的贡献11.使用MPI_全部减少函数及其对计算时间的贡献包含在表的第五行中11剩余的操作，如稀疏矩阵汇编和MHFEM离散化的各种操作，只占总计算时间的一小部分，甚至小于线性系统求解器所需的通信时间。

7.1未来的工作

所呈现的数据结构仅限于符合和同种类的非结构化网格。虽然前者是一个核心限制，很难在实现中推广，但后者可以更容易地放宽。

为了支持任意二维多边形网格，必须定义一个实体拓扑，其中子实体的数量是在运行时而不是在编译时指定的。这将对内部产生影响子实体存储层类，它必须显式存储运行时属性，而不是依赖编译时属性。相关超级存储层类已经可以处理运行时属性，因为在编译时不知道超实体的数量，所以子实体存储层课程可以很容易地进行调整。

因此，三维弱均匀网格（例如具有三角形或四边形面的棱柱网格）可以用具有所有编译时属性的单元拓扑表示，但其面是一般多边形。一般的三维多面体网格也可以由一组关联矩阵表示，假设多面体及其面之间的链接是强制性的。多面体的拓扑可以像多边形一样实现。

当然，必须解决多边形、多面体和弱均匀网格的初始化问题，并且必须适当调整使用广义网格实体的所有算法。还有一些与性能和效率有关的问题，特别是在处理GPU上的一般多面体网格时，但只有在完成一般实现后才能解决这些问题。