CN112989683A - 一种sph的向量化并行计算方法及装置 - Google Patents

一种sph的向量化并行计算方法及装置 Download PDF

Info

Publication number
CN112989683A
CN112989683A CN202110418253.0A CN202110418253A CN112989683A CN 112989683 A CN112989683 A CN 112989683A CN 202110418253 A CN202110418253 A CN 202110418253A CN 112989683 A CN112989683 A CN 112989683A
Authority
CN
China
Prior art keywords
particle
data
particles
target
adjacent
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202110418253.0A
Other languages
English (en)
Inventor
杨灿群
龙思凡
郭晓威
范小康
李超
崔英博
黄春
唐滔
彭林
方建滨
张鹏
夏泽宇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202110418253.0A priority Critical patent/CN112989683A/zh
Publication of CN112989683A publication Critical patent/CN112989683A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本申请涉及一种SPH的向量化并行计算方法、装置、计算机设备和存储介质。所述方法包括：通过将邻近粒子搜索范围内粒子的原始AoS数据重新组织为SoA数据，使得SoA数据符合所使用的SIMD指令集要求，通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将所述粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算，确定邻近粒子；由SIMD指令集封装的应力计算函数，向量化并行计算目标粒子所受的应力信息，根据应力信息更新目标粒子在下一时间步长的粒子信息。本发明充分利用了SIMD指令集向量化计算的性能，提高了SPH程序在CPU内的运行效率。

Description

一种SPH的向量化并行计算方法及装置
技术领域
本申请涉及计算流体力学技术领域，特别是涉及一种SPH的向量化并行计算方法、装置、计算机设备和存储介质。
背景技术
随着计算技术的飞速发展，越来越多的实验可以通过计算机模拟得出具有实际意义的结果。因此，计算机模拟技术已经深入的应用到各个领域，比如爆炸与冲击、水文治理、侵彻等等需要耗费大量资源的领域，并且取得了很好的效果。但是通过计算机模拟会存在一些问题，其中比较普遍的是要解决模拟时间的问题。由于计算机模拟是通过数学离散的方法将一系列的偏微分方程离散为大量能够进行近似计算的数学方程组，因此在求解近似方程组则会耗费大量的计算资源。
光滑粒子动力学(Smoothed Particle Hydrodynamics,SPH)是一种拉格朗日型无网格粒子方法，已经成功地应用到了工程和科学的众多领域。SPH使用粒子离散及代表所模拟的介质，并且基于粒子体系估算和近似介质运动的控制方程。但是该算法需要计算邻近粒子之间的相互作用，因此当计算涉及到数以万计的粒子时会呈现出计算效率低，无法在规定的时间内满足计算的精度要求。
因此，现有技术存在计算效率低的问题。
发明内容
基于此，有必要针对上述技术问题，提供一种能够提高SPH程序计算效率的SPH的向量化并行计算方法、装置、计算机设备和存储介质。
一种SPH的向量化并行计算方法，所述方法包括：
根据实验要求设定SPH计算的求解区域并设置边界条件；
根据目标粒子的位置，确定所述目标粒子在所述求解区域中的邻近粒子搜索范围；
读取所述邻近粒子搜索范围内粒子的原始AoS数据，将所述原始AoS数据重新组织为SoA数据；其中，所述原始AoS数据通过结构体链表形式以粒子为单位存储所述粒子的多个属性变量；所述SoA数据通过多个数组分别存储粒子的每个属性变量；
通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将所述粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算所述目标粒子和所读取数据的粒子之间的距离，根据所述距离确定所述目标粒子的邻近粒子；所述粒子数据为所述SoA数据中的粒子变量；
根据所述邻近粒子，由SIMD指令集封装的应力计算函数，向量化并行计算所述目标粒子所受的应力信息，根据所述应力信息更新所述目标粒子在下一时间步长的粒子信息。
在其中一个实施例中，还包括：将所述求解区域分为多个元胞；
在一维度以所述目标粒子所在元胞和前后相邻的2个元胞为邻近粒子搜索范围，在二维度以所述目标粒子所在元胞和相邻8个元胞为邻近粒子搜索范围，在三维度以所述目标粒子所在元胞和相邻26个元胞为邻近粒子搜索范围。
在其中一个实施例中，还包括：读取所述邻近粒子搜索范围内粒子的原始AoS数据为：
[x1,y1,z1…m1],[x2,y2,z2…m2]…[xn,yn,zn…mn]
其中，x,y,z…m为所述粒子的多种变量，n为所述邻近粒子的总数；
将所述原始AoS数据重新组织为SoA数据，所述SoA数据中的数组元素为：
X[n]＝x1,x2,x3…xn
Y[n]＝y1,y2,y3…yn
Z[n]＝z1,z2,z3…zn
…
M[n]＝m1,m2,m3…mn
其中，X[n]为粒子1…n的变量X组成的数组；Y[n]为粒子1…n的变量Y组成的数组；Z[n]为粒子1…n的变量Z组成的数组；M[n]为粒子1…n的变量M组成的数组。
在其中一个实施例中，还包括：读取所述邻近粒子搜索范围内粒子的原始AoS数据，将所述原始AoS数据重新组织为SoA数据之后，根据SMID指令集对应的寄存器宽度和粒子数据所占比特数，得到每读取一次数据可以得到的向量化并行计算的数据组数；
判断所述邻近粒子的总数是否能被所述数据组数整除，当不能整除时，将余数个粒子按照原始AoS数据形式进行常规串行计算。
在其中一个实施例中，还包括：通过SIMD技术的向量化部件，一次读取所述邻近粒子的一组粒子数据，由SIMD指令集封装的应力计算函数向量化并行计算所述目标粒子受到的应力；
将所述一组粒子数据的计算结果进行归约操作；
根据归约后的信息得到所述一组粒子对所述目标粒子施加的应力的累积结果；比如需要对所有邻近粒子的插值结果进行累加；
根据所有邻近粒子的所述累积结果得到应力的最终累积结果；
根据所述最终累积结果更新所述目标粒子在下一时间步长的粒子信息。
在其中一个实施例中，还包括：通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将所述粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算所述目标粒子和所读取数据的粒子之间的距离，当所述距离小于预先设定的光滑长度2h时，标记所述粒子为邻近粒子；所述粒子数据为所述SoA数据中的粒子变量。
在其中一个实施例中，还包括：所述粒子的多个属性变量包括：粒子坐标、粒子速度、粒子密度、粒子所受内力、粒子所受外力、粒子所受人工粘滞力、粒子能量、粒子能量变化率。
一种SPH的向量化并行计算装置，所述装置包括：
初始化模块，用于根据实验要求设定SPH计算的求解区域并设置边界条件；
邻近粒子搜索范围确定模块，用于根据目标粒子的位置，确定所述目标粒子的邻近粒子搜索范围；
SoA数据获取模块，用于读取所述邻近粒子搜索范围内粒子的原始AoS数据，将所述原始AoS数据重新组织为SoA数据；其中，所述原始AoS数据通过结构体链表形式以粒子为单位存储所述粒子的多个属性变量；所述SoA数据通过多个数组分别存储粒子的每个属性变量；
邻近粒子确定模块，用于通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将所述粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算所述目标粒子和所读取数据的粒子之间的距离，根据所述距离确定所述目标粒子的邻近粒子；所述粒子数据为所述SoA数据中的粒子变量
粒子信息更新模块，用于根据所述邻近粒子，由SIMD指令集封装的应力计算函数，向量化并行计算所述目标粒子所受的应力信息，根据所述应力信息更新所述目标粒子在下一时间步长的粒子信息。
一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现以下步骤：
根据实验要求设定SPH计算的求解区域并设置边界条件；
根据目标粒子的位置，确定所述目标粒子在所述求解区域中的邻近粒子搜索范围；
读取所述邻近粒子搜索范围内粒子的原始AoS数据，将所述原始AoS数据重新组织为SoA数据；其中，所述原始AoS数据通过结构体链表形式以粒子为单位存储所述粒子的多个属性变量；所述SoA数据通过多个数组分别存储粒子的每个属性变量；
通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将所述粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算所述目标粒子和所读取数据的粒子之间的距离，根据所述距离确定所述目标粒子的邻近粒子；所述粒子数据为所述SoA数据中的粒子变量；
根据所述邻近粒子，由SIMD指令集封装的应力计算函数，向量化并行计算所述目标粒子所受的应力信息，根据所述应力信息更新所述目标粒子在下一时间步长的粒子信息。
一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：
根据实验要求设定SPH计算的求解区域并设置边界条件；
根据目标粒子的位置，确定所述目标粒子在所述求解区域中的邻近粒子搜索范围；
读取所述邻近粒子搜索范围内粒子的原始AoS数据，将所述原始AoS数据重新组织为SoA数据；其中，所述原始AoS数据通过结构体链表形式以粒子为单位存储所述粒子的多个属性变量；所述SoA数据通过多个数组分别存储粒子的每个属性变量；
通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将所述粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算所述目标粒子和所读取数据的粒子之间的距离，根据所述距离确定所述目标粒子的邻近粒子；所述粒子数据为所述SoA数据中的粒子变量；
根据所述邻近粒子，由SIMD指令集封装的应力计算函数，向量化并行计算所述目标粒子所受的应力信息，根据所述应力信息更新所述目标粒子在下一时间步长的粒子信息。
上述SPH的向量化并行计算方法、装置、计算机设备和存储介质，通过将邻近粒子搜索范围内粒子的原始AoS数据重新组织为SoA数据，使得SoA数据符合所使用的SIMD指令集要求，通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将所述粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算目标粒子和所读取数据的粒子之间的距离，根据距离确定目标粒子的邻近粒子；根据邻近粒子的信息，由SIMD指令集封装的应力计算函数，向量化并行计算目标粒子所受的应力信息，根据应力信息更新目标粒子在下一时间步长的粒子信息。本发明通过将SPH程序开源框架中的数据重新组织为SoA数据形式，将邻近粒子搜索模块和粒子应力计算模块分别封装，使得SPH程序兼容SIMD的指令集，充分利用了SIMD指令集向量化计算的性能，提高了SPH程序的运行效率。
附图说明
图1为一个实施例中SPH的向量化并行计算方法的流程示意图；
图2为一个实施例中SPH算法基本流程图；
图3为一个实施例中固壁边界虚粒子生成示意图；
图4为一个实施例中小元胞划分求解区域示意图；
图5为一个实施例中原始AoS数据和SoA数据存储内存示意图；
图6为一个实施例中SIMD多数据流计算示意图；
图7为一个实施例中粒子数据组织为SoA形式后的示意图；
图8为另一个实施例中SPH的向量化并行计算方法的流程示意图；
图9为一个实施例中中SPH的向量化并行计算装置的示意图；
图10为一个实施例中计算机设备的内部结构图。
具体实施方式
为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。
本申请提供的SPH的向量化并行计算方法，可以应用于如下应用环境中。根据实验要求设定SPH计算的求解区域并设置边界条件；根据目标粒子的位置，确定目标粒子在求解区域中的邻近粒子搜索范围，将邻近粒子搜索范围内粒子的原始AoS数据重新组织为SoA数据，使得SoA数据符合所使用的SIMD指令集要求，通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将所述粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算目标粒子和所读取数据的粒子之间的距离，根据距离确定目标粒子的邻近粒子；根据邻近粒子的信息，由SIMD指令集封装的应力计算函数，向量化并行计算目标粒子所受的应力信息，根据应力信息更新目标粒子在下一时间步长的粒子信息。
在一个实施例中，如图1所示，提供了一种SPH的向量化并行计算方法，包括以下步骤：
步骤102，根据实验要求设定SPH计算的求解区域并设置边界条件。
光滑粒子动力学(Smoothed Particle Hydrodynamics,SPH)是一种拉格朗日型无网格粒子方法，SPH使用粒子离散及代表所模拟的介质，基于粒子体系估算和近似介质运动的控制方程。SPH算法基本流程图如图2所示，开始后设置几何参数和初始条件，进行邻近粒子搜索后，计算粒子速度和坐标，计算粒子应力后，更新粒子信息，直到超过最大迭代次数，结束循环。
SPH程序中的边界条件设置是通过数值差分或者拟合沿着边界粒子继续向外扩展生产虚粒子补足边界的缺陷问题；同样，对于固壁边界，使用镜面虚粒子法，也就是以固壁为镜面对称生成虚粒子，并且确保虚粒子与实粒子携带有相同的物理量。如图3所示，在问题域内(空心圆圈所在区域)，坐标轴为固壁边界，在外侧生成对称的虚粒子(图中的黑点表示)。这样就会防止因为物理穿透，造成计算上的错误。
步骤104，根据目标粒子的位置，确定目标粒子在求解区域中的邻近粒子搜索范围。
将求解区域划分为若干个小元胞，如图4所示，通过元胞将求解区域分为非邻近粒子搜索范围和邻近粒子搜索范围，在邻近粒子搜索范围内根据距离阈值2h确定目标粒子P1的邻近粒子。一维度只需搜索目标粒子所在元胞和2个相邻元胞内的邻近粒子，二维度只需搜索目标粒子所在元胞和8个相邻元胞内的邻近粒子，三维度只需搜索目标粒子所在元胞和26个相邻元胞内的邻近粒子，所要搜索的元胞范围即为邻近粒子搜索范围。
步骤106，读取邻近粒子搜索范围内粒子的原始AoS数据，将原始AoS数据重新组织为SoA数据。
在SPH开源框架中，粒子数据以原始AoS数据形式存储，原始AoS数据通过结构体链表形式以粒子为单位存储粒子的多个属性变量，如图5所示，若粒子由X，Y，Z，M四种变量，原始AoS数据在存储粒子数据时，先存储粒子1的变量值x1，y1，z1，m1，接着存储粒子2的变量值x2，y2，z2，m2，接着是粒子3、粒子4……，一直到粒子n的变量xn，yn，zn，mn。
将原始AoS数据重新组织为SoA数据，SoA数据通过多个数组分别存储粒子的每个属性变量，使得多个粒子的某种变量的数据内存相邻，如图5所示，若粒子由X，Y，Z，M四种变量，SoA数据则有四个数组，分别为X[n]，Y[n]，Z[n]，M[n]，其中，X[n]＝x1,x2,x3...xn，Y[n]＝y1,y2,y3...yn，Z[n]＝z1,z2,z3..zn，M[n]＝m1,m2,m3...mn，每个数组中的元素内存相邻。
步骤108，通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算目标粒子和所读取数据的粒子之间的距离，根据距离确定目标粒子的邻近粒子。
SIMD全称Single Instruction Multiple Data，单指令多数据流，能够复制多个操作数，并把它们打包在大型寄存器的一组指令集，如图6所示，区别于标量计算每次执行两个数据的运算，SIMD每次执行多组数据的运算。其历代版本也一直在更新换代，比较具有代表性如Intel推出的版本为MMX，SSE，AVX以及AVX-512等，直到后续的SSE2、SSE3、AVX2。其中，AVX2指令集的寄存器宽度为256位，AVX-512指令集的寄存器宽度为512位。
粒子数据为SoA数据中的粒子变量，通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将粒子数据批量加载到内存，若使用的是AVX2指令集，粒子数据为单精度浮点数，由于单精度浮点数所占内存为4byte，即32位，AVX2指令集256位的寄存器可以存储8个单精度浮点的粒子数据，如图7所示，8个粒子的粒子数据为一组，因此，每读取一次数据，可以读到8个粒子的粒子数据并行计算，理论上，通过使用SIMD技术则可以将计算次数减少到[n/8]次；若使用的是AVX-512指令集，粒子数据为单精度浮点数，由于单精度浮点数所占内存为4byte，即32位，AVX-512指令集512位的寄存器可以存储16个单精度浮点的粒子数据，因此，每读取一次数据，可以读到16个粒子的粒子数据做计算，理论上，通过使用SIMD技术则可以将计算次数减少到[n/16]次。
距离计算函数的公式为：
Figure BDA0003026777300000081
其中dij为粒子i到粒子j的距离，x,y,z依次为粒子的三维坐标。通过由SIMD指令集对距离计算函数进行重新封装，使得距离计算函数可以并行处理多个粒子数据。当所计算的粒子与目标粒子的距离dij小于某一预先设定的光滑长度2h时，确认所计算的粒子为邻近粒子。
步骤110，根据邻近粒子，由SIMD指令集封装的应力计算函数，向量化并行计算目标粒子所受的应力信息，根据应力信息更新目标粒子在下一时间步长的粒子信息。
应力计算函数也是由SIMD指令集重新封装的，使得目标粒子所受应力的计算也是向量化并行的，进一步提高算法的计算效率。根据应力信息更新目标粒子在下一时间步长的粒子信息。
上述SPH的向量化并行计算方法中，通过将邻近粒子搜索范围内粒子的原始AoS数据重新组织为SoA数据，使得SoA数据符合所使用的SIMD指令集要求，通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算目标粒子和所读取数据的粒子之间的距离，根据距离确定目标粒子的邻近粒子；根据邻近粒子的信息，由SIMD指令集封装的应力计算函数，向量化并行计算目标粒子所受的应力信息，根据应力信息更新目标粒子在下一时间步长的粒子信息。本发明通过将SPH程序开源框架中的数据重新组织为SoA数据形式，将邻近粒子搜索模块和粒子应力计算模块分别封装，使得SPH程序兼容SIMD的指令集，充分利用了SIMD指令集向量化计算的性能，提高了SPH程序的运行效率。
在其中一个实施例中，还包括：将求解区域分为多个元胞；
在一维度以目标粒子所在元胞和前后相邻的2个元胞为邻近粒子搜索范围，在二维度以目标粒子所在元胞和相邻8个元胞为邻近粒子搜索范围，在三维度以目标粒子所在元胞和相邻26个元胞为邻近粒子搜索范围。
在其中一个实施例中，还包括：读取邻近粒子搜索范围内粒子的原始AoS数据为：
[x1,y1,z1…m1],[x2,y2,z2…m2]…[xn,yn,zn…mn]
其中，x,y,z…m为粒子的多种变量，n为邻近粒子的总数；
将原始AoS数据重新组织为SoA数据，SoA数据中的数组元素为：
X[n]＝x1,x2,x3…xn
Y[n]＝y1,y2,y3…yn
Z[n]＝z1,z2,z3…zn
…
M[n]＝m1,m2,m3…mn
其中，X[n]为粒子1…n的变量X组成的数组；Y[n]为粒子1…n的变量Y组成的数组；Z[n]为粒子1…n的变量Z组成的数组；M[n]为粒子1…n的变量M组成的数组。
在其中一个实施例中，还包括：读取邻近粒子搜索范围内粒子的原始AoS数据，将原始AoS数据重新组织为SoA数据之后，根据SMID指令集对应的寄存器宽度和粒子数据所占比特数，得到每读取一次数据可以得到的向量化并行计算的数据组数；
判断邻近粒子的总数是否能被数据组数整除，当不能整除时，将余数个粒子按照原始AoS数据形式进行常规串行计算。
为了确保进行向量化并行计算的数据符合SIMD指令集要求，必须确保SIMD指令集所读取的寄存器中存储的是完整的数据，以便能完全兼容数据对齐原则。对于尾部的数据，不需要做向量化并行处理，而是在原始串行数据的基础上直接进行串行计算即可。例如寄存器宽度为256位，粒子数据为单精度浮点数，每个数据占32位，得到向量化并行计算的数据组数为8，若粒子总数为405，405除以8的余数为5，那么最后5组粒子不进行向量化并行，而是按传统方式串行计算。由于尾部数据的量很少，即使是串行计算也可以很快完成，对整体算法的计算效率影响不大。
在其中一个实施例中，还包括：通过SIMD技术的向量化部件，一次读取邻近粒子的一组粒子数据，由SIMD指令集封装的应力计算函数向量化并行计算目标粒子受到的应力；将一组粒子数据的计算结果进行归约操作；根据归约后的信息得到一组粒子对目标粒子施加的应力的累积结果；根据所有邻近粒子的累积结果得到应力的最终累积结果；根据最终累积结果更新目标粒子在下一时间步长的粒子信息。
归约操作是向量化并行计算的必要操作，使得向量化元素的计算结果可以进行进一步的整合计算。
在其中一个实施例中，还包括：通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算目标粒子和所读取数据的粒子之间的距离，当距离小于预先设定的光滑长度2h时，标记粒子为邻近粒子；粒子数据为SoA数据中的粒子变量。
在其中一个实施例中，还包括：粒子的多个属性变量包括：粒子坐标、粒子速度、粒子密度、粒子所受内力、粒子所受外力、粒子所受人工粘滞力、粒子能量、粒子能量变化率。
应该理解的是，虽然图1的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。
在一个具体实施例中，如图8所示，提供了一种SPH的向量化并行计算方法，包括：对SPH算法进行改进，将邻近粒子打包、邻近粒子封装和邻近粒子批量计算作SIMD并行处理，提高算法运行效率。
在一个实施例中，如图9所示，提供了一种SPH的向量化并行计算装置，包括：初始化模块902、邻近粒子搜索范围确定模块904、SoA数据获取模块906、邻近粒子确定模块908和粒子信息更新模块910，其中：
初始化模块902，用于根据实验要求设定SPH计算的求解区域并设置边界条件；
邻近粒子搜索范围确定模块904，用于根据目标粒子的位置，确定目标粒子的邻近粒子搜索范围；
SoA数据获取模块906，用于读取邻近粒子搜索范围内粒子的原始AoS数据，将原始AoS数据重新组织为SoA数据；其中，原始AoS数据通过结构体链表形式以粒子为单位存储粒子的多个属性变量；SoA数据通过多个数组分别存储粒子的每个属性变量；
邻近粒子确定模块908，用于通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算目标粒子和所读取数据的粒子之间的距离，根据距离确定目标粒子的邻近粒子；粒子数据为SoA数据中的粒子变量；
粒子信息更新模块910，用于根据邻近粒子，由SIMD指令集封装的应力计算函数，向量化并行计算目标粒子所受的应力信息，根据应力信息更新目标粒子在下一时间步长的粒子信息。
邻近粒子搜索范围确定模块904还用于将求解区域分为多个元胞；在一维度以目标粒子所在元胞和前后相邻的2个元胞为邻近粒子搜索范围，在二维度以目标粒子所在元胞和相邻8个元胞为邻近粒子搜索范围，在三维度以目标粒子所在元胞和相邻26个元胞为邻近粒子搜索范围。
SoA数据获取模块906还用于读取邻近粒子搜索范围内粒子的原始AoS数据为：
[x1,y1,z1…m1],[x2,y2,z2…m2]…[xn,yn,zn…mn]
其中，x,y,z…m为粒子的多种变量，n为邻近粒子的总数；
将原始AoS数据重新组织为SoA数据，SoA数据中的数组元素为：
X[n]＝x1,x2,x3…xn
Y[n]＝y1,y2,y3…yn
Z[n]＝z1,z2,z3…zn
…
M[n]＝m1,m2,m3…mn
其中，X[n]为粒子1…n的变量X组成的数组；Y[n]为粒子1…n的变量Y组成的数组；Z[n]为粒子1…n的变量Z组成的数组；M[n]为粒子1…n的变量M组成的数组。
邻近粒子确定模块908还用于根据SMID指令集对应的寄存器宽度和粒子数据所占比特数，得到每读取一次数据可以得到的向量化并行计算的数据组数；判断邻近粒子的总数是否能被数据组数整除，当不能整除时，将余数个粒子按照原始AoS数据形式进行常规串行计算。
粒子信息更新模块910还用于通过SIMD技术的向量化部件，一次读取邻近粒子的一组粒子数据，由SIMD指令集封装的应力计算函数向量化并行计算目标粒子受到的应力；将一组粒子数据的计算结果进行归约操作；根据归约后的信息得到一组粒子对目标粒子施加的应力的累积结果；根据所有邻近粒子的累积结果得到应力的最终累积结果；根据最终累积结果更新目标粒子在下一时间步长的粒子信息。
邻近粒子确定模块908还用于通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算目标粒子和所读取数据的粒子之间的距离，当距离小于预先设定的光滑长度时，标记粒子为邻近粒子；粒子数据为SoA数据中的粒子变量。
关于SPH的向量化并行计算装置的具体限定可以参见上文中对于SPH的向量化并行计算方法的限定，在此不再赘述。上述SPH的向量化并行计算装置中的各个模块可全部或部分通过软件、硬件及其组合来实现。上述各模块可以硬件形式内嵌于或独立于计算机设备中的处理器中，也可以以软件形式存储于计算机设备中的存储器中，以便于处理器调用执行以上各个模块对应的操作。
在一个实施例中，提供了一种计算机设备，该计算机设备可以是终端，其内部结构图可以如图10所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口、显示屏和输入装置。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种SPH的向量化并行计算方法。该计算机设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该计算机设备的输入装置可以是显示屏上覆盖的触摸层，也可以是计算机设备外壳上设置的按键、轨迹球或触控板，还可以是外接的键盘、触控板或鼠标等。
本领域技术人员可以理解，图10中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
在一个实施例中，提供了一种计算机设备，包括存储器和处理器，该存储器存储有计算机程序，该处理器执行计算机程序时实现上述方法实施例中的步骤。
在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述方法实施例中的步骤。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(ROM)、可编程ROM(PROM)、电可编程ROM(EPROM)、电可擦除可编程ROM(EEPROM)或闪存。易失性存储器可包括随机存取存储器(RAM)或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM(SRAM)、动态RAM(DRAM)、同步DRAM(SDRAM)、双数据率SDRAM(DDRSDRAM)、增强型SDRAM(ESDRAM)、同步链路(Synchlink)DRAM(SLDRAM)、存储器总线(Rambus)直接RAM(RDRAM)、直接存储器总线动态RAM(DRDRAM)、以及存储器总线动态RAM(RDRAM)等。
以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。

Claims (10)

1.一种SPH的向量化并行计算方法，其特征在于，所述方法包括：
根据实验要求设定SPH计算的求解区域并设置边界条件；
根据目标粒子的位置，确定所述目标粒子在所述求解区域中的邻近粒子搜索范围；
读取所述邻近粒子搜索范围内粒子的原始AoS数据，将所述原始AoS数据重新组织为SoA数据；其中，所述原始AoS数据通过结构体链表形式以粒子为单位存储所述粒子的多个属性变量；所述SoA数据通过多个数组分别存储粒子的每个属性变量；
通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将所述粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算所述目标粒子和所读取数据的粒子之间的距离，根据所述距离确定所述目标粒子的邻近粒子；所述粒子数据为所述SoA数据中的粒子变量；
根据所述邻近粒子，由SIMD指令集封装的应力计算函数，向量化并行计算所述目标粒子所受的应力信息，根据所述应力信息更新所述目标粒子在下一时间步长的粒子信息。
2.根据权利要求1所述的方法，其特征在于，所述根据目标粒子的位置，确定所述目标粒子在所述求解区域中的邻近粒子搜索范围包括：
将所述求解区域分为多个元胞；
在一维度以所述目标粒子所在元胞和前后相邻的2个元胞为邻近粒子搜索范围，在二维度以所述目标粒子所在元胞和相邻8个元胞为邻近粒子搜索范围，在三维度以所述目标粒子所在元胞和相邻26个元胞为邻近粒子搜索范围。
3.根据权利要求2所述的方法，其特征在于，读取所述邻近粒子搜索范围内粒子的原始AoS数据，将所述原始AoS数据重新组织为SoA数据包括：
读取所述邻近粒子搜索范围内粒子的原始AoS数据为：
[x1,y1,z1···m1],[x2,y2,z2···m2]···[xn,yn,zn···mn]
其中，x,y,z···m为所述粒子的多种变量，n为所述邻近粒子的总数；
将所述原始AoS数据重新组织为SoA数据，所述SoA数据中的数组元素为：
X[n]＝x1,x2,x3···xn
Y[n]＝y1,y2,y3···yn
Z[n]＝z1,z2,z3···zn
···
M[n]＝m1,m2,m3···mn
其中，X[n]为粒子1···n的变量X组成的数组；Y[n]为粒子1···n的变量Y组成的数组；Z[n]为粒子1···n的变量Z组成的数组；M[n]为粒子1···n的变量M组成的数组。
4.根据权利要求3任意一项所述的方法，其特征在于，读取所述邻近粒子搜索范围内粒子的原始AoS数据，将所述原始AoS数据重新组织为SoA数据之后，包括：
根据SMID指令集对应的寄存器宽度和粒子数据所占比特数，得到每读取一次数据可以得到的向量化并行计算的数据组数；
判断所述邻近粒子的总数是否能被所述数据组数整除，当不能整除时，将余数个粒子按照原始AoS数据形式进行常规串行计算。
5.根据权利要求4任意一项所述的方法，其特征在于，根据所述邻近粒子，由SIMD指令集封装的应力计算函数，向量化并行计算所述目标粒子所受的应力，根据所述应力信息更新所述目标粒子在下一时间步长的粒子信息，包括：
通过SIMD技术的向量化部件，一次读取所述邻近粒子的一组粒子数据，由SIMD指令集封装的应力计算函数向量化并行计算所述目标粒子受到的应力；
将所述一组粒子数据的计算结果进行归约操作；
根据归约后的信息得到所述一组粒子对所述目标粒子施加的应力的累积结果；
根据所有邻近粒子的所述累积结果得到应力的最终累积结果；
根据所述最终累积结果更新所述目标粒子在下一时间步长的粒子信息。
6.根据权利要求5所述的方法，其特征在于，通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将所述粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算所述目标粒子和所读取数据的粒子之间的距离，根据所述距离确定所述目标粒子的邻近粒子；所述粒子数据为所述SoA数据中的粒子变量，包括：
通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将所述粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算所述目标粒子和所读取数据的粒子之间的距离，当所述距离小于预先设定的光滑长度时，标记所述粒子为邻近粒子；所述粒子数据为所述SoA数据中的粒子变量。
7.根据权利要求1至6任一项所述的方法，其特征在于，所述粒子的多个属性变量包括：粒子坐标、粒子速度、粒子密度、粒子所受内力、粒子所受外力、粒子所受人工粘滞力、粒子能量、粒子能量变化率。
8.一种SPH的向量化并行计算装置，其特征在于，所述装置包括：
初始化模块，用于根据实验要求设定SPH计算的求解区域并设置边界条件；
邻近粒子搜索范围确定模块，用于根据目标粒子的位置，确定所述目标粒子的邻近粒子搜索范围；
SoA数据获取模块，用于读取所述邻近粒子搜索范围内粒子的原始AoS数据，将所述原始AoS数据重新组织为SoA数据；其中，所述原始AoS数据通过结构体链表形式以粒子为单位存储所述粒子的多个属性变量；所述SoA数据通过多个数组分别存储粒子的每个属性变量；
邻近粒子确定模块，用于通过SIMD技术的向量化部件，一次读取多个粒子的粒子数据，将所述粒子数据批量加载到内存，通过由SIMD指令集封装的距离计算函数向量化并行计算所述目标粒子和所读取数据的粒子之间的距离，根据所述距离确定所述目标粒子的邻近粒子；所述粒子数据为所述SoA数据中的粒子变量；
粒子信息更新模块，用于根据所述邻近粒子，由SIMD指令集封装的应力计算函数，向量化并行计算所述目标粒子所受的应力信息，根据所述应力信息更新所述目标粒子在下一时间步长的粒子信息。
9.一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至7中任一项所述方法的步骤。
10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至7中任一项所述的方法的步骤。
CN202110418253.0A 2021-04-19 2021-04-19 一种sph的向量化并行计算方法及装置 Pending CN112989683A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110418253.0A CN112989683A (zh) 2021-04-19 2021-04-19 一种sph的向量化并行计算方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110418253.0A CN112989683A (zh) 2021-04-19 2021-04-19 一种sph的向量化并行计算方法及装置

Publications (1)

Publication Number Publication Date
CN112989683A true CN112989683A (zh) 2021-06-18

Family

ID=76341104

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110418253.0A Pending CN112989683A (zh) 2021-04-19 2021-04-19 一种sph的向量化并行计算方法及装置

Country Status (1)

Country Link
CN (1) CN112989683A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114357907A (zh) * 2022-01-07 2022-04-15 中国空气动力研究与发展中心计算空气动力研究所 一种适用于拉格朗日型粒子类数值模拟的并行方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120249560A1 (en) * 2011-04-01 2012-10-04 Paul Frederick Cilgrim Dickenson Parallel computation of matrix problems
CN104991999A (zh) * 2015-06-17 2015-10-21 大连理工大学 一种基于二维sph的溃坝洪水演进模拟方法
CN106484532A (zh) * 2016-09-19 2017-03-08 华东师范大学 面向sph流体模拟的gpgpu并行计算方法
CN110321161A (zh) * 2019-06-26 2019-10-11 中国人民解放军国防科技大学 使用simd指令的向量函数快速查表法、系统及介质
CN110673877A (zh) * 2019-08-22 2020-01-10 成都信息工程大学 一种基于手动向量化的并行计算方法
EP3751444A1 (en) * 2019-06-11 2020-12-16 Dassault Systemes Simulia Corp. Computer simulation of physical fluids on irregular spatial grids with stabilized explicit numerical diffusion

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120249560A1 (en) * 2011-04-01 2012-10-04 Paul Frederick Cilgrim Dickenson Parallel computation of matrix problems
CN104991999A (zh) * 2015-06-17 2015-10-21 大连理工大学 一种基于二维sph的溃坝洪水演进模拟方法
CN106484532A (zh) * 2016-09-19 2017-03-08 华东师范大学 面向sph流体模拟的gpgpu并行计算方法
EP3751444A1 (en) * 2019-06-11 2020-12-16 Dassault Systemes Simulia Corp. Computer simulation of physical fluids on irregular spatial grids with stabilized explicit numerical diffusion
CN110321161A (zh) * 2019-06-26 2019-10-11 中国人民解放军国防科技大学 使用simd指令的向量函数快速查表法、系统及介质
CN110673877A (zh) * 2019-08-22 2020-01-10 成都信息工程大学 一种基于手动向量化的并行计算方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
周煜坤等: "基于CUDA的大规模流体实时模拟", 《计算机应用与软件》 *
范小康等: "基于ARM SVE的光滑粒子流体动力学SIMD加速方法" *
范小康等: "基于ARM SVE的光滑粒子流体动力学SIMD加速方法", 《计算机工程与科学》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114357907A (zh) * 2022-01-07 2022-04-15 中国空气动力研究与发展中心计算空气动力研究所 一种适用于拉格朗日型粒子类数值模拟的并行方法
CN114357907B (zh) * 2022-01-07 2023-03-21 中国空气动力研究与发展中心计算空气动力研究所 一种适用于拉格朗日型粒子类数值模拟的并行方法

Similar Documents

Publication Publication Date Title
Zhou et al. An active learning radial basis function modeling method based on self-organization maps for simulation-based design problems
Ramamurti et al. A parallel implicit incompressible flow solver using unstructured meshes
CN115016951B (zh) 流场数值模拟方法、装置、计算机设备和存储介质
CN113821878B (zh) 一种改善高超声速气动热热流分布异常的计算方法和装置
CN110750933A (zh) 一种耦合Lagrange质点和Euler方法的精确界面追踪处理方法
Bruno et al. Implementing the beam and warming method on the hypercube
Bocharov et al. Implicit method for the solution of supersonic and hypersonic 3D flow problems with Lower-Upper Symmetric-Gauss-Seidel preconditioner on multiple graphics processing units
CN111930932A (zh) 网络空间安全领域知识图谱表示学习方法和装置
CN112989683A (zh) 一种sph的向量化并行计算方法及装置
CN112765871B (zh) 一种基于曲线坐标的并行粒子追踪方法和装置
CN111930491B (zh) 一种全局通信优化加速方法、装置和计算机设备
Khimich et al. Numerical study of the stability of composite materials on computers of hybrid architecture
Fujita et al. Acceleration of element-by-element kernel in unstructured implicit low-order finite-element earthquake simulation using openacc on pascal gpus
CN115630559B (zh) 一种基于粒子网格适配算法的流固耦合方法以及装置
Chen et al. Developing an advanced neural network and physics solver coupled framework for accelerating flow field simulations
US9600446B2 (en) Parallel multicolor incomplete LU factorization preconditioning processor and method of use thereof
Nie et al. Development of an object-oriented finite element program with adaptive mesh refinement for multi-physics applications
Marshall et al. Performance evaluation and enhancements of a flood simulator application for heterogeneous hpc environments
Sanfui et al. Symbolic and numeric kernel division for graphics processing unit-based finite element analysis assembly of regular meshes with modified sparse storage formats
CN115048877A (zh) 环境流体分布信息生成方法、装置、电子设备及存储介质
Zapata et al. A GPU parallel finite volume method for a 3D Poisson equation on arbitrary geometries
Thompson et al. Simest: Technique for Model Aggregation with Considerations of Chaos
CN112287622B (zh) 基于链接方向人工压缩的快速湍流数值模拟方法和装置
Tomczak et al. Parallel accelerator for discontinuous Galerkin method for Navier-Stokes equations
Tumeo et al. A flexible CUDA LU-based solver for small, batched linear systems

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication

Application publication date: 20210618

RJ01 Rejection of invention patent application after publication