1.简介
格子Boltzmann方法（LBM）是计算流体动力学中用于求解Navier–Stokes方程的一种简单算法。格子Boltzmann方法的有效实现受到了计算机科学家群体的极大关注[1,2,三]. LBM有几个特性，从面向性能的算法角度来看，这些特性使它很有趣。浮点操作与内存访问的比率相对较低，因此该方法的性能通常受到带宽限制。LBM使用的变量超过了求解Navier–Stokes方程所需的数学变量，这增加了相对较大的数据流量。尽管如此，如果实施正确，LBM通常被认为是有效的。该算法非常适合大规模并行实现[4]. 时间积分是二阶精度的，尽管它只取决于前一个时间步长[5]. 与有限差分格式不同，LBM的一个特殊特征是局部（即节点）时间积分格式的输入变量数等于输出变量数。这一点很重要，因为它意味着每个输入数据只需要一次。在有限差分格式中，函数值通常必须收集在网格节点的有限邻域中，更新多个网格节点需要相同的数据。近年来，大多数算法和硬件优化都侧重于数据的重用。空间分块和扫描线算法是为加速有限差分方法开发的突出示例[6]. 缓存是一种通用的硬件优化，旨在提高重复使用数据的性能。值得注意的是，LBM并没有从这种优化中受益，因为它不会在单个时间步骤中重用数据。LBM只能间接地从硬件缓存中获益，或者当缓存太大时，它可以保存状态变量，直到在下一个时间步骤中重用它们。这种零数据冗余不一定是LBM的一个缺点。例如，Asinari等人[7]LBM效率低下，因为使用了分布函数，这增加了内存需求，并且应该直观地增加CPU和主内存之间的通信量。Asinari等人[8]建议用宏观变量替换分布函数，这减少了内存占用，但由于计算有限差分时必须重复读取数据，因此反而增加了内存流量[7]. 此外，节省的内存比天真的假设要小：LBM可以用输出数据覆盖每个输入数据，而有限差分方法需要一个源数组和一个目标数组，以避免覆盖相邻网格节点仍然需要的数据。除了这些纯粹的算法差异外，有限差分和LBM实现之间的精度也存在差异。根据模型和验证示例的细节，可以使用有限差分法[9]或LBM[10]胜过其他人。在当前的论文中，我们将不讨论LBM背后的数值建模。相反，我们将完全关注其有效实现。通过这样做，我们利用了该方法的特点，得出了一个与高效实现有限差分方法截然不同的结果。 2.格子Boltzmann方法
在本节中，我们从纯算法的角度描述LBM，而不考虑建模的细节。我们指导读者阅读各种教科书，以了解该方法[11,12,13,14]. LBM是一种求解晶格Boltzmann方程的计算方法，该方程可以在三维中写成： 这里是索引我,j个和k个是表示分布方向的整数值在笛卡尔晶格上移动，节点位于相应位置x个,年和z（z）时间步长表示为t吨.所谓的碰撞算子是所有局部分布的函数但既不依赖于任何其他晶格节点的状态，也不依赖于同一节点在除t吨因此可以看出，晶格节点之间的通信仅发生在(1)所谓的流媒体。右侧，即所谓的碰撞，完全是局部的。还可以看出，晶格玻尔兹曼方程是从变量。每个基准仅使用一次。 LBM可以用于不同的速度集，其中索引我,j个和k个取不同的值。在不失一般性的情况下，我们将我们限制在技术上最相关的情况，即27个速度和我,j个,.
3.格子Boltzmann方法的实现
LBM在算法上与时域有限差分技术有着表面上的相似性，可以以类似的方式实现它。简单的方法将利用两个数据阵列和并重写(1)作为： 在完成一个时间步骤之后，碰撞后分布可以被重新定义为新的碰撞前分布： 此更新算法表示为AB模式[15,16]. 在时域有限差分法中，为了避免相邻节点的数据被覆盖，必须进行这种更新。在LBM中，这是不必要的，因为冲突将一组输入映射到相同数量的输出。因此，允许通过输出直接覆盖输入，并且不需要二级数组。这个概念的简单实现需要在单独的步骤中移动分布： 我们希望实现LBM的方式是，只需要一个数组，而流和冲突在一个步骤中组合在一起，即我们正在寻找一种只需要单个数组来保存分布并在每个时间步中只读取和写入每个数据一次的算法。这种技术称为就地流，文献中提出了几种方法。就地流算法的早期示例是Mattila等人的交换算法[17]和Latt[18]. 交换算法在一个步骤中将流和冲突结合在一起，但它不会通过输出覆盖本地冲突的输入。因此，交换算法需要按一定的顺序处理晶格节点。这是一个严重的缺点，因为它限制了算法在指令顺序明确的串行计算中的应用。该方法不适用于大规模并行计算。 就地流的另一种替代方法是压缩网格，它大约节省一半的内存[19]. 压缩网格背后的思想是将矩形模拟域向各个方向扩展一行节点，并在网格节点的处理过程中保持一定的顺序，在更新第一行节点后，网格节点被写入备用的节点行。处理完第一行中的节点后，就不再需要它们了，可以将第二行节点写入第一行节点，以此类推，直到处理完所有行。压缩网格方法仅适用于矩形域，并且需要固定的节点更新顺序。 贝利[16]研究了图形卡上LBM的实现，其中网格节点由数千个线程并行以随机顺序更新。他提出了一种称为AA-pattern的算法，该算法通过区分奇偶时间步长来用输出覆盖输入（参见图1). 在奇数时间步长中，分布不会被流式传输，本地分布会被直接覆盖。在偶数时间步中，从相邻节点获取分布，并在碰撞后以相反的方向写回相邻节点。以这种方式，分布每秒钟移动两个晶格间距。AA-pattern要求代码实现两次，一次用于偶数时间步，一次用来奇数时间步。 AA-pattern解决了交换算法的并发问题，但它是为在完整的矩阵实现上运行而设计的。计算流体动力学的实际应用通常需要较高的几何灵活性和局部网格细化。为此，通常希望将LBM实现为稀疏矩阵方法。AA-pattern允许在稀疏矩阵上实现，但它需要指向所有相邻节点的指针，这意味着内存占用相对较大，我们介绍并讨论了Esoteric Twist数据结构的特性，该数据结构解决了与AA-pattern类似的并发问题，但在使用间接寻址时会大大节省内存。
4.神秘扭曲
Esoteric Twist（简称EsoTwist）算法结合了以下所需功能：
交换算法实现了除最后两个功能外的所有功能，而AA-模式实现了除最新功能外的全部功能。AA模式必须在每二个时间步长中访问所有相邻节点，这意味着如果使用间接寻址，它必须首先从索引数组中读取指向所有相邻节点的指针。下面我们展示了EsoTwist方法在内存占用和间接寻址的内存读取方面更经济。在本节中，我们首先介绍一般情况下的EsoTwist方法。
Esoteric Twist这个名称来源于一个非直观（深奥）的观察结果，即如果与冲突前的读取相比，以相反（扭曲）的顺序写回分布，则可以消除流步骤。与AA-pattern的情况一样，线程安全的就地流式处理的关键是，冲突操作符只向它从中提取输入的内存写入。我们假设现代计算机可以访问一个缓慢但巨大的内存，该内存容纳所有节点的所有分布，以及一个快速但较小的内存，在该内存中可以执行与单个节点冲突相关的所有操作。我们将前者称为主内存，后者称为寄存器。我们假设，与AA-模式一样，属于一个节点的分布在冲突之前从主内存传输到寄存器，在冲突之后写回主内存。我们假设主存和寄存器之间的传输开销很大，应该最小化。由于分布应该移动，因此它们要么需要在冲突前从相邻节点（拉模式）中提取，要么必须在冲突后发送到相邻节点（推模式）。在EsoTwist方法中，所有沿负方向移动的分布我,j个,当所有分布沿正向移动时我,j个,被按下。该过程如所示图2对于二维情况。中的黄色节点图2就是我们考虑一个更新步骤的那个。节点上的碰撞取决于问传入分发及其返回问传出分发。只有传入的分布正向移动x个-和年-方向是从黄色节点本身读取的。从图中所示的正方向上的相应邻域读取沿负方向移动的分布。碰撞后，将分布写回读取反向移动分布的位置。正如AA模式中的情况一样，每个分布都会交换位置，并且在碰撞后分布会朝相反的方向移动。为了使分布发生移动，该方法必须区分奇数和偶数时间步长，如图2这可以通过不同的方式实现。与AA-pattern不同，EsoTwist对于奇数和偶数时间步长不一定需要两种方法实现。通过将分布存储在阵列结构中，即通过根据分布的移动方向将分布分组在阵列中，可以获得该方法的更简单的实现。完成一个时间步后，所有指向数组的指针都将与指向相反方向移动的数组的指针交换。然后，如前所述执行下一个时间步骤。然而，不一定需要数组的结构。EsoTwist方法也可以像AA-pattern一样通过显式区分奇数和偶数时间步长来实现。算法1在一维中显示了EsoTwist算法的概念。三维的完整算法如附录A. 这里值得注意的是，在中提出了一个与EsoTwist几乎相同的算法[20]. 然而，作者显然忽略了这样一个事实，即他们可以简单地将碰撞和交换方向结合起来，形成一个单一的步骤。相反，他们在单独的步骤中交换了分布，这不必要地使数据流量加倍。| 算法1：用于阵列格式结构中具有三种速度的一维LBM的基本EsoTwist算法。星号标记碰撞后的分布。首都F类表示从主存中的分布读取分布的寄存器（f）因为碰撞。分布正向运动与平稳分布在索引处读取n个而分布向负方向移动是从正方向的邻居处读取的。碰撞后，分布正向移动在最初读取反向分布的时隙处，以正方向写入邻居。负向分布现在已写入位置n个这样，沿负方向的分布和沿正方向的分布在时间步长内都移动了一个节点的距离。在所有节点发生冲突后，交换指向数组的指针就足够了和执行所有节点的流。 |
![计算05 00019 i001]() |
4.1. 间接寻址
在几乎所有计算流体动力学的实际应用中，使用矩形笛卡尔网格离散计算域是不现实的。要离散任意形状的域，稀疏矩阵是最灵活的选择。在稀疏矩阵实现中，每个节点都必须通过指针或等效的索引访问其邻居，索引必须存储在额外的数组中。例如，如果使用AA-pattern在稀疏矩阵上实现具有27个速度的三维LBM，则每个节点必须在偶数时间步长中找到节点的26个邻居。AA-pattern是高效的，因为它只能每隔一个时间点找到邻居。不过，它必须存储指向这些邻居的指针。指针通常占用与标量数据相同的内存量，因此26个指针占用的内存量基本上与分布本身相同。这种内存可以通过应用指针追逐来显著减少。例如，使用AA-pattern，就足以存储六个轴向的指针，并通过查询（例如，正方向的邻居）来访问对角邻居x个-其相邻方向为负年-方向。应用这种技术可以将每个节点的指针数量从26个减少到6个。如果节点沿阵列的一个基本方向对齐，则此数字可以进一步减少到四个。然而，这些优化都降低了方法的灵活性。与间接寻址相关的网格沿其中一个主轴对齐通常意味着网格在生成后无法修改，这是自适应模拟的一个障碍。也就是说，一旦节点沿一维对齐，为了消除此方向的指针，稀疏矩阵网格上的节点可以随意插入和删除的功能就会丢失。因此，在稀疏矩阵网格的分析中，我们不会考虑沿着一个基本方向对齐数据。
使用指针跟踪，AA-pattern必须在每个节点存储到六个邻居的链接(的数据值N个节点）。对于具有27个离散速度的LBM，在没有进一步元数据的情况下，总内存需求为数据值。由于AA-pattern仅每隔一个时间步长传输数据，因此它需要读取和使用读奇数和以偶数时间步长读取。有趣的是，虽然指针跟踪减少了内存占用，但它并没有减少读取次数。所有指向所有邻居的26个指针必须每隔一个时间点读取一次。
与AA-pattern相比，EsoTwist减少了内存占用和读取次数。这是因为EsoTwist在空间上是不对称的。只能访问两个维度中正向的邻居。在三维空间中，使用27个速度，一个节点访问的所有分布都存储在该节点的索引和周围立方体上八分位的七个邻居的索引中。因此，通过应用指针跟踪，必须存储链接，并且每个时间步骤都需要读取和写入。因此，总内存需求是（比AA-图案少9%）。数据读取量甚至比AA-pattern少25%，并且数据写入量相同。在他们的分析中，Wittmann等人[三]得出的结论是，AA-pattern和EsoTwist的数据流量在已知的流算法中最低，但他们的分析忽略了EsoTwost需要的指针和读取次数比AA-patter少。 除了直接节省内存外，由于使用间接寻址时需要较少的虚节点，因此与AA-模式相比，还可以间接节省内存（请参阅图3). EsoTwist和AA-pattern都将与每个节点相关的部分分布存储在一些相邻节点上。即使节点本身不属于模拟域，也必须为这些节点分配内存。在AA-pattern中，必须分配现有节点的所有26个邻居，而在EsoTwist中，只需分配周围多维数据集上八分位中的七个邻居。根据几何体的复杂性，这可以通过较少数量的重影节点来实现可观的节约。 应该注意的是，上述分析忽略了缓存对读写次数的影响，因此在此阶段无法得出关于性能的一般结论。现代CPU通常一次读取完整的缓存线，这样可以在数据对齐时提高性能。然而，由于在当前小节中，我们考虑了稀疏矩阵网格，因此我们不能保证节点按任何有利的顺序排列，也不能保证我们会从缓存中受益。上述分析适用于数据完全不一致的最坏情况。
4.2. 变体
尽管指针是机器级所有软件中不可或缺的元素，但一些流行的编程语言（例如Fortran）并不明确支持指针。这是非常不幸的，因为硬件在纳秒内高效执行的简单操作必须由更复杂、当然效率更低的操作来模拟。尽管如此，没有指针也可以实现EsoTwist。一种可能性是将沿相反方向移动的分布阵列组合成一个长度为两倍的阵列。以D3Q27速度集为例，我们将有一个长度数组N个用于分配和十三个长度数组对于其他分布。例如，让是分布的数组和现在，让我们和对于奇数时间步长和和对于偶数时间步。可以在不使用指针交换通道的情况下访问分布： 因此，也可以在不支持指针的编程语言中实现EsoTwist。
另一个考虑因素是数组数据类型的结构。在大规模并行硬件（如GPGPU）上，应用单指令多数据操作，可以根据方向将数据存储在不同的阵列中。然而，这并不是串行计算的最佳数据布局，因为所有27个分布都存储在主内存的不同位置，并且没有一个分布共享一条缓存线。在这些机器上，最好将属于某个节点的所有数据存储在一起。盖勒[21]提出了一个折衷方案，增强了EsoTwist在串行计算机上的性能。他建议使用两个容器，一个用于本地分发，另一个用于相邻节点的分发。这样，碰撞所需的一半数据始终可以一起访问。只有静态分布被单独存储。 EsoTwist也已成功应用于块结构网格[22,23]. 块结构可用于增加数据的局部性，并有助于负载平衡并行代码。在块内，不需要指向相邻节点的指针。与EsoTwist连接时，每个块只需访问其周围立方体上八分位中的七个相邻块即可进行数据交换。因此，如果使用指针跟踪，则每个块只需要三个到相邻块的链接。块结构网格可以通过高数据局部性来提高性能，它们大大简化了方法的并行化，但与EsoTwist在稀疏矩阵网格上的应用相比，它们也有一些缺点。块之间接口的数据传输不会自动进行，需要明确的数据交换步骤。通常，还需要至少一个额外的节点层来促进数据交换。就三维立体块的内存占用而言节点，通常需要ghost节点。使用指针追逐，每个块只需要三个到相邻块的链接。相反，稀疏矩阵要求每个节点有三个链路，并且仅在域边界上有虚单元。假设每个数据（分布和指针）在内存中占用相同的大小，当块结构网格占用的内存少于稀疏矩阵表示时，我们可以确定内存消耗中的盈亏平衡点。对于速度为27的方法，我们必须将三个指针等同于块结构网格所需的平均额外数据量： 这解决了。因此，我们看到块大小至少为与稀疏矩阵相比，块结构网格上的节点需要占用更少的内存。此分析没有考虑到块结构网格通常包含不属于域的节点。原则上，只要块中至少有一个节点位于流体域中，就必须完全分配块。因此可以看出，内存经济通常是使用块结构网格的不良动机，至少在与EsoTwist结合时是如此。
EsoTwist数据结构也已在一种变体中实现，该变体可确保所有发行版的内存对齐[24]. 这对于在老一代Nvidia GPU上实施格子Boltzmann方法具有历史重要性[25]其中数据总是以十六个相邻的四字节值的形式读取。EsoStripe中称为变体的技巧是，相邻点的分布以16的距离存储，这样，只要向量处理器同时处理16个节点，数据访问就完全对齐。虽然这导致GeForce GTX 465的性能提高了25%以上，但与非对齐读写相比，在较新的硬件上没有观察到性能提高[24]. 在这种情况下，发现数据对齐对于较新的硬件来说并不重要。 4.3。隐式反弹
EsoTwist有一个有趣的特性，它涉及从计算中排除的晶格节点。如果未接触晶格节点，即如果未读取任何分布且未写入任何分布，则此节点上的填充移动方向将通过指针交换反转。在LBM中，分布方向的反转用于在所谓的简单反弹方案中模拟实体墙。使用EsoTwist反弹将隐式应用于更新期间跳过的所有节点，因此这些节点的行为类似于实心墙。在EsoTwist中应用简单的反弹完全是免费的，因为各个节点甚至不需要被触摸。这种通过省略节点来强制反弹的方式称为隐式反弹。该属性由其他交换方法共享，如Mattila等人提出的方法[17]和Latt[18]和AA-图案[16]. 5.结果
EsoTwist数据结构已在并行CPU系统上成功实现[22,23]和在GPU上[26,27,28]. 这里，我们重点介绍与间接寻址相关的GPU的实现。 5.1. 性能模型
我们的性能模型基于这样的观察，即LBM通常是一种内存带宽有限的算法[三,29,30,31,32]. 27速晶格的EsoTwist方法需要27个分布，必须在每个时间步长中读取和写入这些分布。即使每个节点只存储三个指向邻居的指针，我们仍然需要在每个节点和时间步长中读取七个这样的指针。另一个数字用于指示节点的类型。由于此指示器还用作指向边界条件数组的指针（每个节点可以有其单独的边界条件），因此此数字还必须至少有四个字节。总的来说，我们的方法需要每个节点和时间步长27+7+1=35次读取和27次写入，其总和为=62×4字节=248字节。我们代码的性能以每秒百万节点更新（MNUPS）为单位进行衡量。带宽占用率P（P）计算公式为： 当比较EsoTwist的结果与完整矩阵代码的结果时，我们必须考虑到完整矩阵代码不需要任何指针，因此每个节点的每个时间步中只有28次读取和27次写入。
5.2. 与全矩阵实现的比较
在这里，我们比较了EsoTwist在稀疏矩阵和完整矩阵实现上的性能。为了使设置具有可比性，我们选择了一个具有节点。在第一个示例中，我们使用隐式反弹边界条件模拟管道流动，并使用Bhatnagar Gross Krook碰撞模型[33]. 此示例在旧的Nvidia Tesla C1060 GPU上运行，带宽为102 GB/s。完整矩阵实现获得205 MNUPS（41.2%带宽），稀疏矩阵实现获得191 MNUSP（43.2%带宽）。 在我们的第二个示例中，我们比较了具有相同网格大小但具有级联碰撞内核的模拟的性能[34]在GeForce GTX Titan X（Maxwell GM200）上，这是比特斯拉C1060更新的GPU。GeForce GTX Titan X的理论峰值带宽为GBytes/s。代码的完整矩阵版本使用AB格式[15,16]. 尽管级联内核的计算量更大，但我们在这个新硬件上获得了更好的性能数据。全矩阵版本使用999.425 MNUPS（60.8%带宽）运行，稀疏矩阵版本使用993.514 MNUPS，68.2%带宽。 在这两种情况下都可以观察到，虽然全矩阵在晶格更新方面的性能更高，但稀疏矩阵版本可以更好地利用带宽。执行速度的差异很小。稀疏矩阵版本在过时的特斯拉C1060上获得了全矩阵版本93.2%的执行速度，在更新的GeForce GTX Titan X上获得了99.4%的执行速度。这一结果进一步证实了Linxweiler的观察结果[24]硬件方面的进步降低了高度优化的全矩阵码相对于更灵活的码的优势。 5.3. 各向同性
我们通过比较同一几何体在三个不同方向下的模拟执行速度来探讨EsoTwist方法的各向同性。之所以这样做，是因为此栅格上的节点首先沿一个基本方向填充，这取决于几何体相邻节点的方向是靠近还是远离。因此，该问题探讨了该方法如何对不同级别的数据局部性做出反应。几何图形如所示图4它由三个具有不同横截面的矩形管组成。每条管道有128个网格节点长。管道的横截面为,和网格节点。模拟在Nvidia GeForce GTX Titan GPU上执行，该GPU运行的设备驱动程序版本为331.82，CUDA（计算统一设备架构）工具包版本为6.0。此GPU的理论带宽为GB/秒[35]. 使用的计算模型是单精度的BGK碰撞核。无论域在空间中的方向如何，我们都观察到类似的性能。对于沿x个-轴，我们获得738.5MNUPS（59.1%峰值带宽）；对于沿年-轴，我们获得729.8 MNUPS（58.4%峰值带宽），对于z（z）方向，我们获得727.3 MNUPS（58.2%峰值带宽）。我们的性能数据与Tomczak和Szafran报告的GeForce GTX Titan 407到684 MNUPS之间的性能数据进行了比较[36]使用19速格子Boltzmann方法（D3Q19格子），在平铺布局（类似于块结构）和两个数据集（AB-pattern）中实现双精度，对应于48%到72.6%的带宽（数字直接取自[36]因为我们的性能模型不直接应用于他们的代码）。相对较宽的性能范围[36]是由于不同测试用例中的数据局部性的差异。 5.4. 多孔材料
我们使用EsoTwist在GPU上间接寻址实现LBM，参与了不同计算流体动力学（CFD）方法之间的比较研究。仿真结果如所示[27]. 与求解多孔介质中Navier–Stokes方程的竞争方法相比，我们的方法是最有效的。这里，我们只提供性能值。 我们离散了半径为6864个单分散球体（珠）的随机堆积mm，如所示图5。使用两种不同的分辨率，使用笛卡尔网格离散珠包：40米(图6)和20米(图7). 网格的完整矩阵表示将包含和40个网格节点m和20m种情况。通过使用间接寻址，节点数量减少到（52%）和（46.5%）。因此，在这种情况下，我们可以看到间接寻址消除了大约一半的网格节点。值得注意的是，在没有指针追踪的AA-pattern的情况下，稀疏数据结构与完整矩阵在这种情况下几乎没有任何优势，因为间接寻址所需的指针数量几乎与存储分布所需的内存一样大。在EsoTwist与指针跟踪相结合的情况下，每个网格节点只需存储三个指针。20人m案例实时模拟72s，转化为1600000个时间步长。与之前的模拟不同，当前测试用例应用了二阶精确边界条件来捕获珠子的曲率，并且模拟域被分解并分布在多个GPU上。在两台Nvidia K40c GPU上，我们获得了525.2 MNUPS的性能，在六台特斯拉C1060上，我们得到了448.6 MNUPS。40人在单个K40c GPU上模拟了m个测试用例，并运行了400.000个时间步长。它获得了233.3 MNUPS的性能。然而，这些性能数据并不包括以下事实：多孔介质的模拟以边界为主，二阶精确插值边界条件的应用[26]在单独的步骤中完成。此外，还必须考虑到复杂累积量碰撞核[26]已使用。应用二阶边界条件产生的数据流量与一次碰撞产生的流量大致相同。如果是40m模拟，我们大约边界节点的带宽占用（包括边界条件的应用）计算为27.5%。因此，这些性能数据表明，即使在具有复杂边界条件的实际应用程序中从未获得最佳性能，但性能仍与最佳情况相当。 6.结论
我们提出了在大规模并行硬件上高效线程安全执行LBM的Esoteric Twist数据结构。EsoTwist只需要对每个时间步中的每个数据进行一次读写操作，并且只需要在主内存中的一个位置。当与间接寻址相结合时，EsoTwist要求使用我们已知的稀疏矩阵的任何LBM算法中指向相邻节点的指针数量最少。据我们所知，EsoTowist在任何稀疏矩阵LBM实现中占用的内存最少，并且它需要使用稀疏矩阵的任何已知LBM算法的最小读取次数。通过与Wittmann等人最近的分析进行比较，也证实了这些说法[三]. 事实上，忽略EsoTwist所需的指针和读取数比AA-pattern少，Wittmann等人得出结论，AA-patter和EsoTwost的数据流量是已知流算法中最低的。EsoTwist不仅效率高，而且非常灵活，可以轻松应用于具有复杂边界的实际应用程序。它还与一些出版物中显示的局部网格细化技术结合得非常好[22,23,26,28]. 我们的结果显示，在内存带宽方面（高达～68.2%的峰值带宽），与在最近的Nvidia GPU上使用两种发行版的块结构代码相比，性能相当[36]和全矩阵代码。在执行时间方面，发现稀疏矩阵EsoTwist方法在最近的GPU上获得了99.4%的AB模式性能。值得注意的是，旧硬件上的比率较低（特斯拉C1060为93.2%）。因此可以看出，硬件方面的进步允许代码具有更大的灵活性，而不会造成严重的性能损失。
