shonDy 是一款基于无网格粒子方法的三维多物理场数值仿真软件。该软件通过将流体属性存储于粒子上，消除了仿真过程中复杂网格带来的约束。传统的网格间相互作用转化为粒子与粒子之间的相互作用。这些相互作用具有确定的影响范围，称为影响半径。以某粒子为球心、半径为 r 定义一个球体，位于该球体内的粒子即视为邻近粒子。
为便于说明，下方以二维情形为例进行展示：
以粒子6为例，其邻近粒子的编号为3、5、7和10。
问题分析
为计算粒子间相互作用，首先需要确定每个粒子的邻近粒子。假设共有 n 个粒子，暴力遍历方法的时间复杂度为 O(n²)。在大规模项目中，粒子数量可达数百万乃至数千万，如此高的计算时间是不可接受的。由于该算法本身时间复杂度较高，基于其进行并行优化的收益十分有限。因此，需要一种更高效的方法。
问题的关键在于，每个粒子在搜索邻近粒子时必须遍历所有其他粒子，而其中绝大多数遍历是不必要的。因此，可通过减少每个粒子在查找邻近粒子时需遍历的粒子数量来降低时间复杂度，并在此优化算法的基础上进行并行化。
具体实现
通过将三维空间中的每个粒子映射到三维网格上，将粒子与粒子之间的位置关系转化为网格单元与网格单元之间的关系，从而将粒子遍历转换为网格遍历。
如图所示，每个粒子编号对应其所在网格单元的编号。
为提高计算效率，网格边长通常设置为粒子的影响半径。在计算某粒子的邻近粒子时，只需计算目标粒子与周围26个网格单元内各粒子的距离，并将每个距离与影响半径进行比较，即可确定所有邻近粒子。
在搜索粒子6的邻近粒子时，只需在网格单元07周围的网格单元中的粒子范围内进行搜索，这大幅缩小了搜索范围，从而降低了时间复杂度。
常用并行方法简介
1. 数据并行
该方法本质上是将程序中的数据进行划分，分配给不同的计算单元处理。数据可包括输入和输出。在进行数据划分时，以下三点考量至关重要：
- 负载均衡：每个计算单元处理的数据量和计算负载应大致相等。否则，部分单元可能处于空闲状态，导致计算资源浪费。
- 信息交换：并行计算中的内存配置通常分为两类：共享内存和分布式内存。共享内存通常用于多核处理器，多个处理单元共享同一内存空间，需采取措施防止内存冲突。分布式内存常用于由多处理器组成的系统，每个处理器拥有独立的内存空间，计算过程中需要进行通信。
2. 任务并行
该方法将问题分解为若干任务，分配给不同的处理器执行。任务之间必须相互独立，不存在依赖关系，或依赖关系可以解决。
邻近粒子搜索的并行化
我们采用数据并行方法解决此问题。（为便于说明，以下示例使用两个进程：进程0和进程1）
将粒子均匀分配到不同进程
并行粒子分配 计算每个粒子的网格编号，并将粒子编号与网格编号对应。
粒子网格编号 边界网格中的粒子信息传递：部分网格单元位于两个进程的边界处，例如网格06。在查找邻近粒子时，网格03、07和11中的粒子也可能是网格06中粒子的邻近粒子。由于分布式内存访问模型，进程0（承载网格06）无法获取其邻近网格和粒子的信息，需要与进程1进行通信。进程1的边界网格亦同理。
遍历当前进程中的每个网格，获取其邻近网格，从而获取邻近网格中的粒子信息。依次计算当前网格中粒子与邻近网格中粒子之间的距离。该过程涉及大量计算，粒子间的逻辑操作相同且无需信息交换，非常适合并行处理。
获取每个粒子的邻近粒子编号。
扩展
GPU 加速
当前主流 GPU 具有强大的计算能力，性能瓶颈主要出现在数据读写操作上。为实现高性能，应关注程序的数据结构是否符合 GPU 特性，通常采用连续存储方式。此外，在程序设计时应尽量减少对全局内存的频繁读写。Hilbert 曲线映射
该算法利用 Hilbert 曲线将三维网格映射到一维编号。二维定义：Hilbert 曲线是一条连续参数曲线。通过选取适当的函数绘制一条连续参数曲线，当参数 t 在区间 [0, 1] 内取值时，曲线遍历单位正方形内的每一个点，从而形成一条填满整个空间的曲线。Hilbert 曲线是连续但不可微的。
二维 Hilbert 曲线示意 其本质是一种降维方式：将二维空间转化为一维空间的曲线，可从二维推广至三维。
三维 Hilbert 曲线示意