CN111353261A - 浸没式边界高精度统计动力学流体仿真的gpu优化方法 - Google Patents

Abstract

本发明提供了一种浸没式边界高精度统计动力学流体仿真的GPU优化方法。本发明应用前景广阔，主要有以下几个方面：第一)在较低分辨率场景下，做到相对较高精度的实时仿真，可以实现交互式预览设计。第二)在较高分辨率的大场景下，充分利用多节点多GPU的可扩展性，使得这一类仿真也能非常高效地完成，这对该方法在多领域实际应用的推广帮助极大，例如搭建大规模高性能工业设计平台等。

Description

浸没式边界高精度统计动力学流体仿真的GPU优化方法

技术领域

本发明涉及一种基于格子玻尔兹曼方程的浸没式边界高阶自适应统计动力学方法的单GPU和多节点多GPU大规模并行计算加速优化方法，以实现兼具精度及效率的大规模湍流仿真和流固耦合仿真，属于流体动力学仿真、计算机图形学和大规模并行计算跨学科技术领域。

背景技术

流体仿真对于工业产品设计、医学诊断、可视动画、环境保护、能源等领域都是非常重要的。传统的流体仿真通常用不同类型的离散化方案求解不可压缩的纳维叶-斯托克斯方程(INSE)，如有限差分、有限体积和有限元法，并在不同方面进行了改进。虽然已经有大量可用的离散求解器，其中一些对于湍流的仿真也具有相对较高的精度，但它们都不可避免地需要求解一个较大的全局稀疏线性系统(例如，压力方程，或由基于稳定性原因的隐式公式产生的方程)，这使得这些求解器很难在GPU上并行化进行高效计算，特别是在仿真分辨率非常高的情况下。因此，不需要全局方程的局部求解器的优势逐渐凸显，利用格子玻尔兹曼方程(LBE)建模的统计动力学求解器被认为是高性能流体仿真的一种选择。这一方法可以达到与目前市面上较为成熟的基于三角网格的仿真软件如ANSYS同样的精度(二阶)，但仿真效率远高于后者。

由于其局部动力学特性，这种新的方案非常简单地求解了底层的模型方程，并且很容易在并行系统上进行大规模场景的高分辨率仿真。然而，在该方法的发展过程和过去的多次尝试中(它们主要依赖于BGK-SRT模型和RM-MRT模型)，由于计算精度较低和稳定性较差，这一方法在流体仿真领域中很少被使用。自适应中心矩多弛豫时间(ACM-MRT)模型的发展使统计动力学方法进行高精度稳定流体仿真成为可能。结果表明，这种新的模型对于高效的流体仿真是非常有前途的，不仅可以生成比最先进的纳维叶-斯托克斯求解器更好的可视化效果，还能达到更高的计算精度。此外，它还可以很容易地与浸没式边界(IB)方法相结合，更加灵活地处理浸入流体区域内的静态和动态刚体，甚至复杂几何形状的流固耦合仿真。

然而，一套基于ACM-MRT模型的统计动力学求解器在GPU上的直接实现并没有达到最佳性能。在之前的一些工作中，有一些提升统计动力学求解器的并行性能的尝试，优化主要集中于对内存数据布局的重新排列和计算顺序的调整，但这些优化方案都针对非常旧的BGK-SRT模型或传统的RM-MRT模型，难以达到较高流体仿真的精度和稳定性。此外，在统计动力学方法中，没有GPU优化工作来针对性地提高浸没式边界-格子玻尔兹曼方法的性能，以支持流体内部与各种复杂静态或动态固体的耦合仿真。因此，提高基于ACM-MRT模型的浸没式边界高阶自适应统计动力学方法在GPU平台上的仿真性能有着重要的意义，将为实现单GPU和多节点多GPU大规模场景的高性能仿真以及小规模多场景的实时仿真提供可能性。

发明内容

本发明的目的是：提高基于ACM-MRT模型的浸没式边界高阶自适应统计动力学方法在GPU平台上的仿真性能，加速多场景下流体仿真的过程，以实现单GPU和多节点多GPU大规模场景的高分辨高效仿真，以及较低分辨率下的实时仿真。

为了达到上述目的，本发明的技术方案是提供了一种浸没式边界高精度统计动力学流体仿真的GPU优化方法，其特征在于，包括以下步骤：

步骤1、在GPU平台上进行算法基本实现

根据浸没式边界-统计动力学方法算法，按步骤流程进行基本GPU并行实现，其中，浸没式边界-统计动力学方法算法包括以下步骤：

步骤101、初始化

将微观速度c离散成多个格子速度，c_i表示第i个格子速度，每个格子速度对应一个速度分布函数，则第i个格子速度对应格子离散速度分布方程f_i，在仿真之前，分配内存来存储格子离散速度分布方程f_i，以及密度ρ、速度u和外力g，对于向量场，被存储为向量的线性阵列，表示为一个结构体，所有这些结构体数组都存储在GPU全局内存中；

步骤102、流处理

从时间t时刻开始，将流体网格节点作为并行对象，通过访问GPU全局内存将每个流体网格节点不同方向的格子离散速度分布方程f_i从附近对应流体网格节点并行平流输送到当前流体网格节点；

步骤103、计算物理宏观量

经过流处理后，将流体网格节点作为并行对象，利用各流体网格节点流处理后的格子离散速度分布方程f_i计算得到该流体网格节点的宏观物理量——密度ρ和速度u：

ρ＝∑_if_i，

步骤104、刚体边界处理

采用浸没式边界方法进行刚体边界条件处理，当计算刚体采样点的受力时，并行对象是这些刚体采样点，每个线程分别访问该刚体采样点附近一定区域内的流体网格节点；当传播力时，并行对象是刚体的包围盒中的流体网格节点，每个线程分别计算该流体网格节点所受反作用力，访问该流体网格节点周围一定区域内的全部刚体采样点；

步骤105、碰撞处理

经过刚体边界处理后，用流处理后的格子离散速度分布方程f_i计算得到的宏观物理量密度ρ和速度u和刚体边界处理计算得到合力项G_i进行碰撞过程Ω_i计算；

步骤106、更新过程

根据计算得到的碰撞过程Ω_i来更新得到每个流体网格节点在时间t+1时刻的格子离散速度分布方程f_i，通过并行流体网格节点来实现；

步骤107、区域边界处理

采用回弹方法对区域边界条件进行处理

步骤2、在单GPU系统上，针对算法本身通过单GPU优化技术，以显著提高湍流及流固耦合仿真的性能。

本发明首先提出一种基于内存访问连续性原则的新的参数化内存数据布局模型，以提高浸没式边界-统计动力学求解器的内存访问性能。其次，提出一种新的浸没式边界并行实现方法，同时在区域边界上进行分位置计算实现，以改善数据负载不平衡和判断带来的多支路发散。同时，进一步简化自适应中心矩多弛豫时间模型(ACM-MRT模型)，以减少每个线程的寄存器使用量，从而提高GPU占用率以提高运算性能。

具体而言，该单GPU优化技术包括：基于内存访问连续性原则的参数化内存数据布局模型，以提高浸没式边界-统计动力学求解器的内存访问性能。减少边界处理时各线程负荷失衡和线程发散。在处理流固耦合的边界问题处理时，使用浸没式边界方法来进行刚体边界条件处理，并且以刚体采样点作为并行对象进行实现，而不是在统计动力学方法中得到了广泛应用的回弹方法。这不仅为处理静态和动态刚体的流固耦合创造了灵活性，而且还避免了代码分支，减少了线程束的发散。但在处理整个仿真区域的边界时，依然采用传统的回弹方法，因为区域边界相较于刚体几何形状更加规则，位置更加固定。进一步简化ACM-MRT模型，以减少每个线程的寄存器使用量，从而提高GPU占用率以提高运算性能；

在多节点多GPU的大规模扩展平台系统上，首先将单GPU优化应用于多节点多GPU系统，同时根据多节点多GPU特征，进行多节点多GPU平台实现优化。由于其很好的可扩展性，可以很好地为不同场景要求，不同分辨率要求，不同精度要求的应用服务。

优选地，步骤105中，进行碰撞过程Ω_i计算所用的松弛模型基于自适应中心矩多弛豫时间模型。这一模型相较于传统的巴特纳加尔-格罗斯-克鲁克单弛豫时间(BGK-SRT)模型和基于原始矩的多弛豫时间(RM-MRT)模型，能获得更高的精度和仿真稳定性，但计算过程更为复杂，计算量远远大于传统方法。

优选地，步骤2中，所述单GPU优化技术进一步包括优化GPU线程分配：当流体网格节点被并行化处理时，在GPU块中分配三维线程组，其中线程组的x维度应该是线程束(warp)大小的倍数(线程束大小通常是32)，y和z维度可根据正在使用的GPU的类型进行调节的，以实现对应GPU最佳并行性能。然而，当刚体采样点被并行化时，总是使用一维线程组，进而每个线程通过初始化时刚体采样点存储顺序的索引分配访问。通过合理分配各线程块的线程数目，能有效平衡资源利用，提升GPU运行效率。

内存布局在统计动力学方法的GPU优化中起着至关重要的作用，其目的是最大限度地合并多个线程对全局内存的访问，利用GPU缓存来提高内存访问性能，减少访问时间。则优选地，步骤2中，所述基于内存访问连续性原则的参数化内存数据布局模型，以提高浸没式边界-统计动力学求解器的内存访问性能包括：

对于流体网格节点的内存布局优化，采用带有参数α的CSoA内存布局并应用于流体网格节点的所有向量场，α表示CSoA内存布局中穿插重复向量分量的个数；

构造一个围绕刚体的包围盒，并将包围盒细分为多个边长为l的正方体。对于刚体采样点的内存布局优化，由于边界表面上的刚体采样点可能不能以合理的顺序存储在GPU内存中(取决于如何构建三角网格，以及每个三角形中的采样点如何排列)，即采样点的空间局部性和内存连续性不匹配，这可能会打破内存合并访问的原则，降低缓存利用率。本发明提出了一个带有参数l的优化内存布局模型，抽象表示为H^t(l,α)，通过选取合适的l，使得刚体样本点在空间上更连续，更好地共享附近的流体网格节点，提高缓存利用率：

式中，t_a表示预仿真开始时间，t_c表示预仿真结束时间，

表示t_s到t_c时间段中的迭代次数。通过预仿真，根据l和α的不同配置，计算多个迭代所需的总时间，最终找到最优的参数l和α，以最小化运行时间。

优选地，步骤2中，所述简化自适应中心矩多弛豫时间模型包括：由于受到ACM-MRT模型复杂度的限制，每个线程计算需要过多的寄存器，可同时运行的线程数目大大减少，这使得GPU的占用率大大降低。本发明采用的优化策略之一是简化碰撞过程的代数表达式，从而减少计算所需的局部变量数量。同时，为了进一步减少每个线程所需的变量数，可以通过将碰撞过程的计算分离成n个不同的内核函数来实现。通过选择适当的n，这样的计算分离可以确保每个核函数的所有变量都存储在寄存器中，从而提高了仿真的性能。

优选地，步骤2中，多节点多GPU系统优化由于统计动力学方法是局部的，只需要一个环形区域内的流体网格节点，所以很容易被扩展到多节点多GPU系统。假设有m个GPU(在单个集群节点上，通常m＝4)，流体网格节点将沿一个维度被细分为多个子块(例如，z维，因为流体网格节点按x-y-z顺序存储)，确保每个子块包含几乎相同数量的流体网格节点。对于每个子块，沿着细分维度扩展一个单位格子作为共享面，这个共享面的数据将从存储在相邻GPU中的附近子域中进行更新。同时由于硬件限制，一些多节点多GPU系统不完全支持PCI-E接口，因此需要依赖于CPU内存的通用性。对于刚体的处理，当刚体不移动时，将属于不同子块的刚体采样点缓存到相应的GPU内存中。然而，当刚体在流体场中移动时，由于每个子块的刚体采样点可能会随着时间的推移而变化，因此我们将所有的刚体采样点都分别缓存到每个GPU以提高效率。除此以外，通过使用OpenMPI进行多节点之间高效的数据传输和通信(需要拷贝的数据与单节点多GPU情况相同，即各节点间的共享面)，以进一步提升计算效率和可仿真场景大小。

本发明应用前景广阔，主要有以下几个方面：第一)在较低分辨率场景下，做到相对较高精度的实时仿真，可以实现交互式预览设计。第二)在较高分辨率的大场景下，充分利用多节点多GPU的可扩展性，使得这一类仿真也能非常高效地完成，这对该方法在多领域实际应用的推广帮助极大，例如搭建大规模高性能工业设计平台等。

本发明与现有技术相比的优点在于：

1.本发明所使用的ACM-MRT模型的浸没式边界高阶自适应统计动力学方法在可视化质量和计算精度上能比其他仿真方法达到更高的水平，特别是对于复杂刚体几何的流固耦合仿真和湍流仿真，在此基础上能实现不同场景下的高效仿真。

2.本发明提出一种基于内存访问连续性原则的新的参数化内存布局模型，以提高浸没式边界-统计动力学求解器的内存访问性能。

3.本发明提出一种新的浸没式边界并行实现方法，同时在区域边界上进行分位置计算实现，以改善数据负载不平衡和判断带来的多支路发散。

4.本发明进一步简化ACM-MRT模型GPU实现，以减少每个线程的寄存器使用量，从而提高GPU占用率以提高运算性能。

5.本发明的方法能够在多节点多GPU的大规模扩展平台进行实现，可以很好地为不同场景要求，不同分辨率要求，不同精度要求的应用服务。

附图说明

图1：D3Q27格子结构示意图；

图2：流传输优化前后性能对比示意图；

图3：宏观量计算优化前后性能对比示意图；

图4：浸没式边界方法优化前后性能对比示意图；

图5：碰撞过程优化前后性能对比示意图；

图6：区域边界处理优化前后性能对比示意图；

图7：参数模型选取不同参数性能对比示意图；

图8：高分辨率仿真结果示意图：流固耦合城市污染扩散预测(仿真分辨率：1200×250×840)；

图9：高分辨率仿真结果示意图：流固耦合飞机尾流仿真(仿真分辨率：2100×240×600)；

图10：高分辨率仿真结果示意图：流固耦合汽车尾流仿真(仿真分辨率：2240×320×480)；

图11：低分辨率仿真结果示意图：流固耦合圆环旋转湍流仿真(仿真分辨率：400×200×400)；

图12：低分辨率仿真结果示意图：流固耦合流体对撞仿真(仿真分辨率：400×200×400)；

图13：与传统纳维叶-斯托克斯方法仿真性能对比示意图；

图14：与传统统计动力学方法湍流速度场结果对比示意图(图左(a)为基于BGK-SRT模型统计动力学方法结果，图中(b)基于传统RM-MRT模型统计动力学方法结果，图右(c)基于ACM-MRT模型统计动力学方法结果)；

图15：与传统统计动力学方法碰撞过程仿真性能对比示意图；

图16：与ACM-MRT模型统计动力学方法基本实现仿真性能对比示意图；

图17：三种不同数据内存结构示意图；

图18：刚体采样点与流体网格节点内存排布结构参数模型示意图；

图19：浸没式边界方法实现及优化示意图；

图20：回弹边界方法实现示意图；

图21：多GPU系统实现示意图；

图22：不同数目GPU系统仿真性能对比示意图(图(a)仿真分辨率200×400×200，图(b)仿真分辨率400×400×400)。

具体实施方式

下面结合具体实施例，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。此外应理解，在阅读了本发明讲授的内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。

本发明的目的是提高基于ACM-MRT模型的浸没式边界高阶自适应统计动力学方法在GPU平台上的仿真性能，加速多场景下流体仿真的过程，以实现单GPU和多节点多GPU大规模场景的高分辨高效仿真，以及较低分辨率下的实时仿真。该发明应用前景广阔，主要有以下几个方面。第一，在较低分辨率场景下，做到相对较高精度的实时仿真，可以实现交互式预览设计。第二，在较高分辨率的大场景下，充分利用多节点多GPU的可扩展性，使得这一类仿真也能非常高效地完成，这对该方法在多领域实际应用的推广帮助极大，例如搭建大规模高性能工业设计平台等。

基于ACM-MRT模型的浸没式边界高阶自适应统计动力学方法的基本算法和计算框架介绍如下：

作为INSE的一种替代选择，玻尔兹曼方程(BE)起源于统计力学，通过演化速度分布函数f，通过以更一般的方式描述流体流动，它指定了在位置x和时间t处用微观速度c找到粒子的概率。BE可以写为以下基本形式：

式(1)中，

表示f的梯度，g表示外力，

表示f关于c的梯度。通过使用巴特纳加尔-格罗斯-克鲁克(BGK)模型建模的Ω，BE可以通过零阶矩和一阶矩恢复弱可压缩的纳维叶-斯托克斯方程，该方程可用于在相对较低的速度(远离声速)条件下近似INSE。为了求解这一方程，通常可以将其离散为以下格子玻尔兹曼方程(LBE)的形式，它们具有局部，密度守恒且动量守恒的特点：

f_i(x+c_i,t+1)-f_i(x,t)＝Ω_i(f_i)+G_i (2)

式(2)中，ACM-MRT模型在三维仿真中使用D3Q27晶格结构将微观速度c离散成27个格子速度，c_i表示第i个格子速度；27个格子速度对应27个速度分布函数，第i个格子速度c_i对应的速度分布函数为f_i(x,t)，即f_i；Ω_i和G_i是第i个格子速度的离散碰撞和外力项。物理宏观量密度和速度可以通过以下关于c_i的公式计算得到：

式(3)中，ρ表示密度，u表示速度。由于其简单的代数形式，LBE的求解过程通常可以分为一个流传输过程：

f_i ^*(x,t)＝f_i(x-c_i,t) (4)

在这之后，进行边界条件处理，以及带有外力项G_i的碰撞过程：

f_i(x,t+1)＝f_i ^*(x,t)+Ω_i(f_i ^*)+G_i (5)

式(5)中，f_i ^*为f_i ^*(x,t)；

Ω＝-M^-1DM(f-f^eq)＝-M^-1D(m-m^eq)，其中，Ω包含每个流体网格节点上的所有碰撞结果Ω_i；f是包含分布函数f_i的集合向量，f^eq是包含平衡分布函数f^eq的集合向量，m是f在中心矩空间中的相应分布，m^eq是f^eq在中心矩空间中的相应分布；M(27×27)是由非正交中心矩构造的中心矩空间投影矩阵；D(27×27)是包含不同弛豫速率的对角线矩阵，可分为低阶驰豫速率(D_i≤8)和高阶驰豫速率(D_i>8)，D_i表示不同弛豫速率的对角线矩阵的各个分量。低阶驰豫速率完全由运动粘度等物理参数决定，高阶驰豫速率是根据流体速度反比关系自适应确定的。在计算外力时，通过加入近似到LBE的二阶离散化中，G_i的计算公式如下：

其中,τ是由运动粘度确定的弛豫时间，w_i是各流体网格对应的网格权重，c_s是声速，

表示f_i关于c_i的梯度，

表示f_i ^eq关于c_i的梯度，f_i ^eq表示f_i对应的平衡分布函数。

以上这样的一套算法流程能很好地在GPU系统上被实现。

基于ACM-MRT模型的浸没式边界高阶自适应统计动力学方法在GPU平台上的基本实现介绍如下：

基于以上浸没式边界-统计动力学方法算法，按步骤流程进行基本GPU并行实现：

步骤(1)初始化。在仿真之前，需要分配内存来存储格子离散速度分布方程f_i，以及密度ρ、速度u和外力g。对于向量场(例如f_i)，它们被存储为向量的线性阵列(表示为一个结构体)，所有这些结构体数组都存储在GPU全局内存中。

步骤(2)流处理。从时间t时刻开始，将流体网格节点作为并行对象，通过访问全局内存将每个流体点不同方向的f_i从附近对应节点并行平流输送到当前节点。

步骤(3)计算物理宏观量。经过流处理后，将流体网格节点作为并行对象，利用各节点流处理后的f_i及相关方程计算得到该节点的宏观物理量密度ρ和速度u。

步骤(4)刚体边界处理。采用浸没式边界方法进行刚体边界条件处理。当计算刚体采样点的受力时，并行对象是这些刚体采样点，每个线程分别访问该刚体采样点附近一定区域内的流体网格节点；当传播力时，并行对象是刚体的包围盒中的流体网格节点，每个线程分别计算该流体网格节点所受反作用力，需要访问该流体点周围一定区域内的全部刚体采样点。

步骤(5)碰撞处理。经过刚体边界处理后，用流处理后的f_i，计算得到的宏观物理量密度ρ和速度u和边界处理计算得到合力项G_i进行碰撞过程Ω_i计算。

步骤(6)更新过程。根据计算得到的Ω_i来更新得到每个流体网格节点在时间t+1时刻的f_i，这是通过并行流体网格节点来实现的。

步骤(7)区域边界处理。采用传统回弹方法对区域边界条件进行处理。

单GPU系统优化方法介绍如下：

上述简单的GPU实现相较于传统的纳维叶-斯托克斯方法还是更快的，因为整个算法均是在较为局部的范围内进行的。然而，这样的直接实现仍然存在几个比较严重的问题，通过对GPU实现进行优化，可以达到更高的性能。最终，优化实现结果比直接的GPU实现快5-10倍。具体优化方法如下：

方法(1)：优化GPU线程分配。首先描述GPU线程通常是如何在GPU上分配给流体网格节点和刚体采样点的。通常，当流体网格节点被并行化时，每个线程都会按照空间排列x-y-z顺序分配一个索引，有序访问流体网格节点。对于并行化流体网格节点，在GPU块中分配三维线程组，其中：线程组的x维度应该是线程束(warp)大小的倍数(线程束大小通常是32)；y和z维度是根据正在使用的GPU的类型进行调节的，以实现对应GPU最佳并行性能。然而，当刚体采样点被并行化时，总是使用一维线程组，进而每个线程通过初始化时刚体采样点存储顺序的索引分配访问。通过合理分配各线程块的线程数目，能有效平衡资源利用，提升GPU运行效率。

方法(2)：提高全局内存访问性能。正如之前的许多工作所指出的，由于受到内存总量的限制，统计动力学方法的数据大部分存储在全局内存中，访问全局内存又是极为耗时的。因此，内存数据布局在统计动力学方法的GPU优化中起着至关重要的作用，其目的是最大限度地合并多个线程对全局内存的访问，利用GPU缓存来提高内存访问性能，减少内存访问时间。在GPU体系结构设计中，访问全局内存是以半线程束高度并行的方式完成的，理想情况下，这需要同一个半线程束中的不同线程不仅使用相同的代码指令(SIMD)操作，而且还访问处理连续存储的数据，从而最小化内存访问次数和耗时。然而，在实践中，这种要求很难被严格满足，但往往可以进行优化。因此，GPU上的内存数据布局设计背后的一个关键原则是最大限度地提高内存访问的连续性，这意味着给定任何线程，其在同一线程束中的附近线程应该尽可能连续地访问内存数据。

以下，分别讨论流体网格节点和刚体采样点的内存优化：

(1)流体网格节点的内存数据布局优化。首先讨论如何为流体网格节点设计内存数据布局。在统计动力学方法中，向量场的内存数据布局(例如f_i、u和g)通常分为两种类型：结构体数组(AOS)和数组结构体(SOA)。在AoS布局里，数组中的每个元素是存储某一流体网格节点所有向量组件的结构体，然而在SoA布局里，整个仿真区域内的流体网格节点的每个向量中的其中一个分量首先被集中放在一个一维数组中，然后将这多个一维数组连续存储在内存中。如果根据先前提出的原则，即线程束中的不同线程应该尽可能连续地访问内存数据，可以发现SoA完全满足了对流体网格节点的这种要求，这就是为什么许多现有的工作建议对流体网格节点的向量场进行SoA内存数据布局。为了进一步提高性能，最近的一些工作建议使用具有较好缓存性能的CSoA结构，通过将SoA结构放入到AoS结构中进行组合得到。因此，在流体网格节点的GPU优化中保留了CSoA内存数据布局，用于统计动力学方法中流体网格节点的所有向量场。如果使用α来表示CSoA中穿插重复向量分量的个数，并且β表示存储在原始数组中的向量的维数，那么使用CSoA存储的向量的每个分量都可以通过以下索引方式在线程中访问：

j＝βind(x)+2^αi+ind(x)％2^α (6)

式(6)中，ind(x)表示流体网格节点x在整个空间中的排列序号，j表示所访问数据在内存中的实际位置。

(2)刚体采样点的内存数据布局优化。当刚体物体浸入流体中时，在刚体与流体表面应满足边界条件。如上所述，由于灵活性，采用浸没式边界方法在刚体边界周围施加额外惩罚力来使之满足边界条件，而不是传统的回弹方法，因为回弹方法需要标记刚体采样点和流体网格节点，并将影响流体网格节点的内存布局的空间连续性。在浸没式边界方法中，由于刚体采样点和流体网格节点存储相互独立，不会干扰内存中流体网格节点排列的连续性。因此，在浸没式边界-统计动力学方法中，CSoA内存数据布局仍然是流体网格节点向量场内存的理想选择。但边界表面上的刚体采样点可能不能以合理的顺序存储在GPU内存中(取决于如何构建三角网格，以及每个三角形中的采样点如何排列)，即采样点的空间局部性和内存连续性不匹配，这可能会打破内存合并访问的原则，降低缓存利用率。在每一个时刻的计算中，每个刚体采样点将插值速度，并将其作用力扩展到附近具有相同核函数的流体网格节点。由于刚体采样点的排序可能是任意的，不一定能确保在插值和计算反作用力过程中同一线程束中的线程访问的流体网格节点被连续存储，这违反了数据访问的连续性，这可能会降低内存缓存利用率并引入延迟。

为了优化浸没式边界-统计动力学方法的内存访问合并率和缓存利用率，将刚体采样点沿流体网格节点空间存储顺序(x-y-z序列)重新排序存储在内存中，以便在为每个刚体采样点访问附近的流体网格节点时，使流体网格节点能够更好地被刚体采样点共享。因此，提出了一个刚体采样点的参数内存数据布局。首先，构造一个围绕刚体的包围盒，并将包围盒细分为多个边长为l的正方体，其中l可以是可根据刚体几何结构调节的。由于我们的目标是保持刚体与流体数据访问的连续性，我们试图将刚体块与周围流体网格节点的顺序对齐(即2D中的x-y顺序和3D中的x-y-z顺序)。然后，对于每个正方体块内的刚体采样点，计算每个采样点空间位置的摩顿(Morton)码，并对同一个正方体块里每个Morton码进行z排序，这样将使访问流体网格节点这一操作具有更高的缓存利用率。事实上，l对性能有不同的影响，因为它产生了不同的缓存利用率。当将l从1增加到2，刚体样本点在空间上更连续，这将更好地共享附近的流体网格节点，提高缓存利用率，从而提高内存访问性能。然而，如果进一步增加l，流体数据共享可能会被破坏，从而降低内存访问性能。所以，这个过程最终将根据如何选择切分长度l来决定刚体采样点在内存中如何存储，同时决定最终的性能表现。

方法(3)：减少边界处理时各线程负荷失衡和线程发散。在处理流固耦合的边界问题处理时，使用浸没式边界方法来执行刚体边界条件，而不是在统计动力学方法中得到了广泛应用的回弹方法。这不仅为处理静态和动态刚体与流体耦合创造了灵活性，而且还避免了代码分支，减少了线程束的发散度。但在处理整个仿真区域的边界时，依然采用传统的回弹方法，因为区域边界相对规则固定。以下分别介绍两种边界处理的具体优化方法：

(1)浸没式边界方法：在统计动力学求解器中直接实现浸没式边界方法可能会导致负载不平衡。在作用力传播过程中，如果将流体网格节点作为线程并行对象，流体网格节点上的作用力是通过将附近一定区域内的刚体采样点作用力加权求和来计算的。虽然泊松圆盘采样确保了刚体表面的基本实现均匀采样，但由于局部几何形状的变化，每个2×2×2局部区域可能包含不同数量的刚体采样点，从而导致在靠近刚体边界的不同流体网格节点的邻域内的求和样本数不同。这就造成了负载不平衡，特别是当刚体几何是非常复杂时。基于GPU的浸没式边界方法中的负载不平衡可以通过并行化刚体采样点来大大减少。在相互作用过程中，每个GPU线程只处理一个采样点的计算，每个采样点从附近一定区域的流体网格节点插值速度。在扩散过程中，每个刚体采样点将其作用力重新分配到附近的流体网格节点，并将其添加到相应的速度中。虽然每个刚体采样点的流体网格节点的数量仍然不同，但不同情况的数量却显著减少，求和的数量也显著减少，使负载不平衡问题变得不那么严重。请注意，这样的并行加法可能有冲突，导致不正确的求和。因此，需要在GPU代码中添加原子求和以确保正确性，这不会导致明显的性能降低。

(2)传统回弹方法：虽然上述GPU实现的浸没式边界方法可以显着地减少处理刚体边界时线程负载不平衡的问题，但这一方法通常很难应用于仿真区域边界，所以对于区域边界我们还是需要保持传统的回弹方法进行特殊处理，这会导致由于条件判断语句而产生的线程束发散。由于区域边界形状相对更加规则，同时位置固定，因此判断条件较少且确定，可以通过多个独立的内核函数来分别处理每个边界条件。更具体地说，在三维仿真中，有六个区域面，使用六个核调用来分别处理每个面。对于任何一个面，哪些分布函数需要回弹处理已经预先知道(由于简单和固定的平面几何)，它删除了算法中的条件语句，大大避免了在最初直接的实现中发生的线程束发散。

方法(4)：减少局部资源使用量，提升GPU占用率。GPU的占用率通常定义为线程束实际执行的线程总数除以初始化分配运行的线程总数。有两种情况可能会降低GPU的占用率：当多个线程同时访问内存数据和/或每个线程消耗太多的本地资源(寄存器和本地内存)时。由于内存布局设计，前者可以显著减少，使得GPU占用率增加。然而，由于我们用于处理碰撞过程的ACM-MRT模型包含非常复杂的代数表达式，导致每个线程计算需要过多的寄存器。受到线程块总资源量的限制，可同时运行的线程数目大大减少，这使得GPU的占用率大大降低。这里需要注意的是，内核线程中的局部变量首先存储在寄存器中；如果所有寄存器数量不足，其余的局部变量则会被存储在局部内存中，而局部内存的访问效率将会明显低于寄存器访问效率。因此，我们采用的优化策略之一是简化碰撞的代数表达式，从而减少计算所需的局部变量数。同时，需要进一步减少每个线程所需的变量数，这可以通过将碰撞过程的计算分离成n个不同的内核函数来实现(每个内核启动只计算Ω_i的一部分)。通过选择适当的n，这样的计算分离可以确保每个核函数的所有变量都存储在寄存器中，从而提高了仿真的性能。

方法(5)：构建并优化仿真性能参数模型。以上不同方面的GPU优化有效地提高基于ACM-MRT模型的浸没式边界-统计动力学方法的性能。这些优化似乎更独立，但实际上是相互关联的。在整个算法的许多部分，内存访问是不可避免的。因此，内存数据布局应通过考虑整体仿真性能来确定，这是由之前提到的刚体采样点内存参数l和流体网格节点内存参数α控制的。因此，上述GPU优化可以转化为(非线性)参数模型，抽象表示为H^t(l,α)，优化模型如下：

式(7)中，t_s表示预仿真开始时间，t_c表示预仿真结束时间，

表示t_s到t_c时间段中的迭代次数。通过预仿真，根据l和α的不同配置，计算多个迭代所需的总时间，最终找到最优的参数l和α，以最小化运行时间。

多节点多GPU系统实现及优化方法介绍如下：

由于GPU内存大小的限制，仿真只能在一定的分辨率范围内进行。为了支持大规模场景的高分辨率仿真，需要在多节点多GPU上运行和优化统计动力学方法。注意，由于统计动力学方法是局部的，只需要一个环形区域内的流体网格节点，所以很容易将以上优化方法扩展到多节点多GPU系统。假设有m个GPU(在单个集群节点上，通常m＝4)，首先将流体网格节点沿一个维度细分为多个子块(例如，z维，如果流体网格节点按x-y-z顺序存储)，确保每个子块包含几乎相同数量的流体网格节点。对于每个子块，沿着细分维度扩展一个单位格子作为共享面，这个共享面的数据将从存储在相邻GPU中的附近子域中进行更新。为了复制数据，我们依赖于CPU内存的通用性(由于硬件限制，一些多节点多GPU系统不完全支持PCI-E接口)。之前所提到的单GPU系统优化可以独立地应用于每个块。当刚体被浸入但不移动时，将属于不同子块的刚体采样点缓存到相应的GPU内存中。然而，当刚体在流体场中移动时，由于每个子块的刚体采样点可能会随着时间的推移而变化，因此我们将所有的刚体采样点都分别缓存到每个GPU以提高效率。除此以外，通过使用OpenMPI进行多节点之间高效的数据传输和通信(需要拷贝的数据与单节点多GPU情况相同，即各节点间的共享面)，以进一步提升计算效率和可仿真场景大小。

应用以上单GPU，多节点多GPU系统优化方法，基于ACM-MRT模型的统计动力学方法以及浸没式边界方法在GPU平台上的仿真性能将会大大提高，从而加速多场景下流体仿真的过程，以实现不同场景下的实际应用。

本发明公开了一种基于格子玻尔兹曼方程的统计动力学方法以及浸没式边界方法的单GPU和多节点多GPU大规模并行计算加速优化方法，以实现兼具精度及效率的大规模湍流仿真和流固耦合仿真。流体仿真对于工业产品设计、医学诊断、可视动画、环境、能源等领域都有非常重要的应用。目前，常见的流体仿真通常是通过求解不可压缩的纳维叶-斯托克斯方程(INSE)来实现的，这需要全局求解压力方程等稀疏线性系统，其求解过程很难在GPU系统上高效并行实现。最近，基于统计力学的统计动力学方法已经展现出仿真湍流的精确度和稳定性，具有比众多纳维叶-斯托克斯求解方法更好的可视化质量和更高的计算精度。由于其局部求解的特点，统计动力学求解器具有非常高的并行性。然而，当需要高精度计算或采用浸没边界法处理复杂刚体边界时，统计动力学方法的直接并行化实现的效率仍然不够高。因此，需要针对算法本身提出GPU并行计算优化技术，以显著提高湍流仿真的性能。为了实现这一目标，本发明首先提出一种基于内存访问连续性原则的新的参数化内存数据布局模型，以提高浸没式边界-统计动力学求解器的内存访问性能。其次，提出一种新的浸没式边界并行实现方法，同时在区域边界上进行分位置计算实现，以改善数据负载不平衡和判断带来的多支路发散。同时，进一步简化自适应中心矩多弛豫时间(ACM-MRT)模型，以减少每个线程的寄存器使用量，从而提高GPU占用率以提高运算性能。优化后的仿真方法能够在多节点多GPU的多尺度扩展平台进行实现，可以很好地为不同场景要求，不同分辨率要求，不同精度要求的应用服务。该发明应用前景广阔，主要有以下几个方面。第一，在较低分辨率场景下，做到相对较高精度的实时仿真，可以实现交互式预览设计。第二，在较高分辨率的大场景下，充分利用多节点多GPU的可扩展性，使得这一类仿真也能非常高效地完成，这对该方法在多领域实际应用的推广帮助极大，例如搭建大规模高性能工业设计平台等。

Claims (6)

1.一种浸没式边界高精度统计动力学流体仿真的GPU优化方法，其特征在于，包括以下步骤：

步骤1、在GPU平台上进行算法基本实现

根据浸没式边界-统计动力学方法算法，按步骤流程进行GPU并行实现，其中，浸没式边界-统计动力学方法算法包括以下步骤：

步骤101、初始化

将微观速度c离散成多个格子速度，c_i表示第i个格子速度，每个格子速度对应一个速度分布函数，则第i个格子速度对应格子离散速度分布方程f_i，在仿真之前，分配内存来存储格子离散速度分布方程f_i，以及密度ρ、速度u和外力g，对于向量场，被存储为向量的线性阵列，表示为一个结构体，所有这些结构体数组都存储在GPU全局内存中；

步骤102、流处理

从时间t时刻开始，将流体网格节点作为并行对象，通过访问GPU全局内存将每个流体网格节点不同方向的格子离散速度分布方程f_i从附近对应流体网格节点并行平流输送到当前流体网格节点；

步骤103、计算物理宏观量

经过流处理后，将流体网格节点作为并行对象，利用各流体网格节点流处理后的格子离散速度分布方程f_i计算得到该流体网格节点的宏观物理量——密度ρ和速度u：

ρ＝∑_if_i，

步骤104、刚体边界处理

采用浸没式边界方法进行刚体边界条件处理，当计算刚体采样点的受力时，并行对象是这些刚体采样点，每个线程分别访问该刚体采样点附近一定区域内的流体网格节点；当传播力时，并行对象是刚体的包围盒中的流体网格节点，每个线程分别计算该流体网格节点所受反作用力，访问该流体网格节点周围一定区域内的全部刚体采样点；

步骤105、碰撞处理

经过刚体边界处理后，用流处理后的格子离散速度分布方程f_i计算得到的宏观物理量密度ρ和速度u和刚体边界处理计算得到合力项G_i进行碰撞过程Ω_i计算；

步骤106、更新过程

根据计算得到的碰撞过程Ω_i来更新得到每个流体网格节点在时间t+1时刻的格子离散速度分布方程f_i，通过并行流体网格节点来实现；

步骤107、区域边界处理

采用回弹方法对区域边界条件进行处理

步骤2、在单GPU系统上，针对算法本身通过单GPU优化技术，以显著提高湍流及流固耦合仿真的性能，其中，该单GPU优化技术包括：基于内存访问连续性原则的参数化内存数据布局模型，以提高浸没式边界-统计动力学求解器的内存访问性能；使用浸没式边界方法来进行刚体边界条件处理，并且以刚体采样点作为并行对象进行实现，在处理整个仿真区域的边界时，采用回弹方法，以改善数据负载不平衡和判断带来的多支路发散；进一步简化自适应中心矩多弛豫时间模型，以减少每个线程的寄存器使用量，从而提高GPU占用率以提高运算性能；

在多节点多GPU的大规模扩展平台系统上，首先将单GPU优化技术应用于多GPU系统，同时根据多GPU特征，进行多节点多GPU平台实现优化。

2.如权利要求1所述的一种浸没式边界高精度统计动力学流体仿真的GPU优化方法，其特征在于，步骤105中，进行碰撞过程Ω_i计算所用的松弛模型基于自适应中心矩多弛豫时间模型。

3.如权利要求1所述的一种浸没式边界高精度统计动力学流体仿真的GPU优化方法，其特征在于，步骤2中，所述单GPU优化技术进一步包括优化GPU线程分配：当流体网格节点被并行化处理时，在GPU块中分配三维线程组，其中线程组的x维度是线程束warp大小的倍数，y和z维度根据正在使用的GPU的类型进行调节的，以实现对应GPU最佳并行性能；当刚体采样点被并行化时，使用一维线程组，进而每个线程通过初始化时刚体采样点存储顺序的索引分配访问，通过合理分配各线程块的线程数目，有效平衡资源利用，提升GPU运行效率。

4.如权利要求1所述的一种浸没式边界高精度统计动力学流体仿真的GPU优化方法，其特征在于，步骤2中，所述基于内存访问连续性原则的参数化内存数据布局模型，以提高浸没式边界-统计动力学求解器的内存访问性能包括：

对于流体网格节点的内存布局优化，采用带有参数α的CSoA内存布局并应用于流体网格节点的所有向量场，α表示CSoA内存布局中穿插重复向量分量的个数；

构造一个围绕刚体的包围盒，并将包围盒细分为多个边长为l的正方体，对于刚体采样点的内存布局优化，提出了一个带有参数l的优化内存布局模型，抽象表示为H^t(l，α)，通过选取合适的l，使得刚体样本点在空间上更连续，更好地共享附近的流体网格节点，提高缓存利用率：

式中，t_s表示预仿真开始时间，t_c表示预仿真结束时间，

表示t_s到t_c时间段中的迭代次数，通过预仿真，根据l和α的不同配置，计算多个迭代所需的总时间，最终找到最优的参数l和α，以最小化运行时间。

5.如权利要求1所述的一种浸没式边界高精度统计动力学流体仿真的GPU优化方法，其特征在于，步骤2中，所述简化自适应中心矩多弛豫时间模型包括：简化碰撞过程的代数表达式，从而减少计算所需的局部变量数量；通过将碰撞过程的计算分离成n个不同的内核函数来实现，通过选择适当的n，这样的计算分离确保每个核函数的所有变量都存储在寄存器中，从而提高了仿真的性能。

6.如权利要求1所述的一种浸没式边界高精度统计动力学流体仿真的GPU优化方法，其特征在于，步骤2中，在多节点多GPU的大规模扩展平台系统上，假设有m个GPU，流体网格节点沿一个维度被细分为多个子块，确保每个子块包含几乎相同数量的流体网格节点；对于每个子块，沿着细分维度扩展一个单位格子作为共享面，这个共享面的数据从存储在相邻GPU中的附近子域中进行更新；

对于刚体的处理：当刚体不移动时，将属于不同子块的刚体采样点缓存到相应的GPU内存中；当刚体在流体场中移动时，将所有的刚体采样点都分别缓存到每个GPU以提高效率。

Images