首页 / 资讯中心 / 文章详情

LAMMPS 并行算法深度解析:域分解、Ghost 原子通信、邻居表与并行 FFT 的实现原理

LAMMPS 并行算法深度解析:域分解、Ghost 原子通信、邻居表与并行 FFT 的实现原理 ★ FEATURED ARTICLE
科研科学计算高性能计算【免费下载链接】lammpsPublic development project of the LAMMPS MD software package项目地址https://gitcode.com/gh_mirrors/la/lammps点击查看免费下载LAMMPS 采用基于 MPI 的分布式内存并行架构通过空间域分解domain decomposition ghost 原子通信 Verlet 邻居表三级机制实现对大规模原子体系的强扩展与弱扩展并行。本文以 LAMMPS 开发者手册中的并行算法章节为核心结合仓库中src/下的实际源码comm_brick、comm_tiled、邻居表、KSPACE/FFT 与 OPENMP 包系统讲解域分解方案brick/tiled、forward/reverse 通信模式、邻居表与 bin/stencil 构建算法、PPPM 长程静电力的并行 FFT以及 MPI 内层的 OpenMP 线程并行设计帮助开发者完整理解 LAMMPS 并行的底层原理并指导负载优化实践。一、整体并行架构强扩展与弱扩展LAMMPS 的设计目标是支持 MPI 标准下的分布式内存并行通过域分解把模拟盒拆分为若干互不重叠的子域subdomain每个 MPI 进程拥有own一个子域及其中的原子。其并行化追求两类扩展性强扩展strong scaling同一体系下增加进程数能近似按比例加速弱扩展weak scaling体系规模扩大时计算成本与体系规模成正比。在此之上还可以在同一 MPI 子域内部使用 GPU 或 OpenMP 做第二层并行。手册章节为清晰起见多用 2d 示意图但所述算法在 2d 和 3d 情形完全通用。注意该章节文字与配图改编自 LAMMPS 正式发表论文中并行算法一节见文档中的 论文引用。二、域分解Partitioningbrick 与 tiled 两种模式2.1 模拟盒与两种分解模式LAMMPS 的模拟盒是一个 3d或 2d体积可以是**正交orthogonal或单斜/三斜triclinic**形状正交盒盒边沿x、y、z笛卡尔轴对齐各面均为矩形三斜盒盒边沿三个任意向量对齐各面为平行四边形。每个维度上的盒面可以是周期边界的也可以是非周期的fixed 固定边界原子移出即删除shrink-wrapped盒面随原子包络连续收缩调整。分解方式由comm_style命令选定只有两种comm_style特点适用场景brick默认规则网格所有子域大小形状相同原子密度大致均匀tiled不规则瓦片式划分支持递归二分RCB密度不均、需要粒子数均衡正交盒与三斜盒都可以在模拟中连续变形例如压缩固体、剪切液体处理器子域会随之变形。源码层面src/comm.h 中的Comm基类明确定义了BRICK6 向 stencil 通信与TILED不规则瓦片通信两种 style以及LAYOUT_UNIFORM / LAYOUT_NONUNIFORM / LAYOUT_TILED三种布局src/comm_brick.cpp 与 src/comm_tiled.cpp 分别实现两种分解。处理器网格的划分逻辑则位于 src/domain.cppDomain类负责依据体系密度比与进程数优化procgrid[3]。2.2 非均匀密度下的负载不均衡与三种缓解策略密度不均匀的模型中默认 brick 分解会导致各进程粒子数不均衡整体模拟速率受最慢负载最大进程限制。手册给出了三种官方缓解手段调整处理器网格processor grid 默认基于模拟盒体积、并尝试优化子域的体积/表面积比可用processors命令手动指定Nx Ny Nz平移分解平面balance命令可以平移定义子域大小的并行平面常与更优的处理器网格组合使用递归二分RCB tiled 通信balance的 RCB 模式结合comm_style tiled可产生每个子域粒子数相等的划分。这些实现在源码中一一对应src/balance.cpp 及其imbalance_*系列文件实现了各种平衡算法包括 RCBsrc/comm_tiled.cpp 中包含递归遍历 RCB 切割树rcbinfo数组记录每个进程的切割维度与切割比例cutfrac来判定原子归属、计算子域间重叠的逻辑。手册用一个含真空的 2d 体系12 个 MPI rank演示了策略差异默认分解下多个 rank 分不到原子严重失衡强制1x12x1网格后每个 rank 都有计算量但子域内原子数仍不均且薄片形状增大通信量2x6x1网格并平移分解平面后负载与通信进一步改善RCB 递归二分给出最均衡的划分原子颜色从绿到红表示子域负载优劣。三、Ghost 原子通信forward 与 reverse 两种模式3.1 为什么需要 ghost 原子按选定的分解方案所有 per-atom 数据都分布在 MPI 进程之间——这正是 LAMMPS 能用足够多的 MPI 进程处理超大体系的前提。但计算短程作用力时每个进程不仅需要访问自己拥有的原子还需要邻域子域中的原子信息即 LAMMPS 所说的ghost 原子邻居进程原子的副本保存最多到通信截断距离communication cutoff内的 per-atom 数据。这一机制同时也是周期性边界条件的实现方式位于周期边界另一侧、距离在截断之内的原子同样存为 ghost 原子取自子域的周期复制串行运行时即取自本进程自己。由此带来的一个重要性质是LAMMPS 的力计算不受最小镜像约定minimum image convention限制截断距离可以大于模拟域尺寸的一半。3.2 两阶段3d 为三阶段交换模式由于 ghost 数据必须在每个 MD 时间步或最小化步更新通信模式必须高效。2d 情形下通信分两阶段3d 分三阶段进行规则与不规则分解均适用规则brick情形先沿x方向交换再沿y方向处理器子域网格中每个进程有 4 个邻居。x阶段中rank 1 和 2 把红色阴影区域中的 owned 原子发给 rank 0反之亦然y阶段中rank 3 和 4 把蓝色阴影区域的原子发给 rank 0——这些原子包含它们在x阶段收到的 ghost 原子。rank 0 由此集齐全部 ghost 原子图中实心蓝色角区里的原子会被通信两次后才到达 rank 0。不规则tiled情形两阶段逻辑相似但一个处理器在每个方向上可能有多个邻居。x阶段中 rank 1、2、3 把红色区域原子发给 rank 0且只包含位于 rank 0 子域y上下边界之间的原子y阶段中 rank 4、5、6 把蓝色区域原子发给 rank 0其中可能包含x阶段收到的 ghost 原子但仅限 rank 0 填充其 ±y方向扩展 ghost 区域所必需的。因此 rank 5、6 之间互相收到的 ghost 原子不会转发给 rank 0。要点是尽管不规则分解的通信模式更复杂仍然只需两阶段3d 三阶段、仅与相邻处理器交换即可完成。3.3 forward / reverse 通信的定义与实现细节forward 通信owned 原子的属性发给邻居进程、成为其 ghost 原子的属性。在每个原子迁移、邻居表重建的时间步每个处理器会列出自己是哪些邻居进程 ghost 原子的 owned 原子清单此后直到下次 reneighboring都用这些清单把 per-atom 坐标等数据打包进消息缓冲。reverse 通信计算出的 ghost 原子属性送回原子所属进程典型用途是把 ghost 原子上的部分力累加回 owned 原子的完整力。两阶段顺序与 forward 相反使用同一套原子清单打包/解包力数据缓冲解包时 ghost 力直接加到 owned 原子力上。图中四个蓝色角区原子的力在两个阶段各发送/接收/累加一次owned 处理器最终才拿到完整力。这两个操作在 LAMMPS 内部被广泛复用坐标与力交换之外manybody 势用它共享 per-atom 中间量、刚体积分器用它让 body 内各原子访问 body 属性。手册还列出若干底层细节均能在源码中印证与不同处理器交换时使用MPI_Send()/MPI_Irecv()若处理器与自己交换如单进程跨周期边界造 ghost只做 pack/unpack。src/comm_brick.cpp 的forward_comm()与reverse_comm()正是这种先MPI_Irecv发异步接收、再MPI_Send发数据的经典实现。跨周期边界发送 owned 原子坐标时收发过程中直接加减周期盒长因此建邻居表与算力时无需再做最小镜像。ghost 交换截断通常等于邻居截断但可设得更大如多原子刚体的半直径、扭结作用力的 3 倍以上拉伸键长甚至可以超过周期盒尺寸。规则模式下若截断超出某方向邻居子域范围会沿同一方向做多次交换第一次只发 owned 原子后续各次发送之前收到的 ghost不规则模式则检测扩展 ghost 子域与其他处理器在每个方向上的全部重叠。Comm基类中的cutghost[3]、cutghostuser以及SINGLE/MULTI多集合截断模式src/comm.h第 31–42 行对应了单截断 / 多集合截断两种通信模式forward_comm()/reverse_comm()为纯虚函数并被Pair、Bond、Fix、Compute、Dump等重载src/comm.h供各组件复用同一通信基础设施。四、邻居表Neighbor ListsVerlet 列表、bin 与 stencil4.1 Verlet 列表与 skin 距离为高效算力每个处理器构建 Verlet 风格邻居表枚举所有i,j对i ownedj owned 或 ghost中距离小于邻居截断的原子对。实现要点邻居表存储在多页multiple-page数据结构中每页是一段连续内存为多个i原子存放邻居j的向量可整块分配/释放。源码对应 src/my_page.cpp 中的MyPage页分配器与 src/neigh_list.cpp 中的NeighList类。邻居表通常是 LAMMPS 中内存占用最大的数据结构。邻居表每隔若干步从零重建一次其间每步重复使用。邻居截断距离为 $R_n R_f \Delta_s$$R_f$ 是势函数定义的最大力截断$\Delta_s$ 是skin 距离允许假设原子相邻步间不会移动太远而多步复用列表。通常当任一原子自上次重建以来移动了半 skin 距离时触发重建相关选项由neigh_modify命令调整原子排序则由atom_modify控制。4.2 重建步原子迁移、周期性包裹与空间排序执行 reneighboring 的时间步上移出所属子域的原子先经通信迁移到新进程周期性边界条件也仅在这些步执行保证每个原子被分配到正确进程迁移后owned 原子存入连续向量每隔若干步对 owned 原子做空间排序spatial sort先用 bin 分箱再重排使同 bin 原子在向量中相邻以改善算力与建表时的缓存效率。4.3 bin 网格与 stencil 搜索为线性时间构建本地邻居表模拟域上叠加规则 3d或 2dbin 网格。每个处理器存储覆盖子域 邻居截断 $R_n$的 bin 集合。bin无需对齐处理器边界每个维度的整数 bin 数是按整个模拟盒尺寸整除适配的。最常见的构建是半列表half list每对i,j只存一次。为此使用围绕中心 bin 的预计算stencil——一组x,y,z整数偏移列出足够近、可能包含任何距中心 bin 内原子i在 $R_n$ 范围内的邻居j的 bin。半列表的 stencil 可以是不对称的。手册给出的 stencil 数值bin 尺寸为 $\frac{1}{2}R_n$ 时维度正交三斜2d13 个 stencil bin15 个3d63 个含原 bin 平面内 13 个z方向上下两层各 25 个75 个三斜 stencil 多出的 bin 是因为 bin 平铺的是三斜域的包围盒本身相对于模拟盒并不周期故判定哪些i,j对入表的逻辑略有调整。构建算法先遍历 owned ghost 原子并分配 bin用整型向量在每个 bin 内建原子索引链表再三重嵌套循环——owned 原子i→i所在 bin 的 stencil → 每个 stencil bin 内的j原子含 ghost——若 $r_{ij} R_n$ 则把j加入i的邻居向量。4.4 邻居表构建的扩展选项手册列出的构建选项与源码/命令对应关系bin 尺寸可调$\frac{1}{2}R_n$ 在多数典型情形最优——bin 更小则循环开销增大bin 更大则距离计算变多。bin 更小时 2d stencil 更接近半圆3d 半球方形角上的 bin 因距离过远被剔除多列表与多种 stencil不同截断算力 vs 如 RDF 诊断、r-RESPA 按距离分层的子集、full 列表i,j存两次用更大的对称 stencil、部分枚举 ghost 邻居的列表——后两者被各类 manybody 势使用入表判据除距离外还可依赖两粒子直径如多分散颗粒pair_style hybrid需为主邻居表按子势生成子列表每个子列表只含相应原子子集间的i,j对异构截断体系如溶胶-胶体体系colloid/colloid、colloid/solvent、solvent/solvent 截断差异巨大、粒径比可达 10–100 倍的多分散颗粒。LAMMPS 提供高效算法对不同截断使用不同 stencil 并裁剪出仅跨越截断球面的 bin还有更新的多 stencil 多 bin 尺寸方法可高效处理任意粒径比体系neighbor nsq小/稀疏体系或兜底方法无分箱、子域内全 $O(N^2)$ 双重循环newton pair开关的影响newton pair off时半列表包含所有j为 ghost 的对newton pair on时仅当j的z或z相等时y、再相等时x坐标更大才入表。对均匀致密体系这使 owned 原子固定从同一方向扇区取邻居邻居表长度更一致、降低负载不均衡风险。五、长程相互作用的并行实现PPPM 与 3d 并行 FFT5.1 PPPM 的五阶段流程带电体系的长程库仑作用由 FFT 基 PPPM 方法kspace_style pppm及其变体计算库仑作用拆为短程部分实空间、邻居表内成对循环与长程部分倒易空间、kspace 风格计算。PPPM 在模拟盒上叠加 3d 规则 FFT 网格流程为每个原子的点电荷插值到邻近网格点前向 3d FFT倒易空间卷积一个或多个逆 3d FFT由原子附近网格点的电场值插值回原子得到力。5.2 网格的域分解与 grid ghost无论哪种空间分解每个处理器都拥有其子域内的砖形brick网格点集合。两次插值都用到围绕每个原子的网格点 stencil为此每个处理器额外保存包围其砖块的若干层ghost 网格点。网格点值的 forward/reverse 通信与原子数据通信方式类似owned 网格点上的电场发给邻居成为其 ghost 点值ghost 点上的电荷值发送并累加回 owned 点。三斜盒情形FFT 网格平面平行于盒面但电荷/电场的网格映射在约化坐标中进行把倾斜盒概念上视为单位立方体stencil 与 FFT 操作不变不过达到给定精度所需的 FFT 网格比正交盒更大。5.3 -pencil 分解的并行 FFT并行 3d FFT 的通信量相对其计算量很大。3d FFT 实现为沿x、y、z三个方向的 1d FFT 序列因此网格不能像原子那样按 3 维分解只能按 1 维平面或 2 维铅笔分解步骤之间需要转置使每个 MPI 进程拥有的网格部分恰好完整覆盖其待做 1d FFT 的方向。LAMMPS 采用图中所示的pencil 分解初始最左每处理器拥有子域内的同色网格砖brick-to-pencil 通信转为x方向 1d 铅笔布局左二随后用配置的 FFT 库对每条完整拥有的铅笔做 1d FFTpencil-to-pencil 通信转为y铅笔相当于转置x、y维做y方向 1d FFT最后一次y→z铅笔转置 z方向 1d FFT完成前向 FFT数据保持z-pencil 布局供卷积逆 FFT 按相反顺序执行 1d FFT 与通信最终网格值布局与初始砖布局相同。每次 FFT 内部通信brick-to-pencil、pencil-to-pencil、pencil-to-brick都是把一种tiling把 3d 网格的小砖形子集分配给各处理器的覆盖方案转换为另一种。LAMMPS 配套的fftMPI库支持任意 tiling 定义因此不规则域分解可直接使用tiling 转换用点对点通信实现——每个处理器只向少数处理器发送数据初始 tiling 中每个 tile 只与目标 tiling 的少数几个 tile 重叠。虽然理论上可用单次MPI_Alltoall()做全体 $P$ 个进程集合通信但通常慢得多。对图中 brick/pencil tiling甚至不需要全体集合通信$8^3$ 网格512 网格单元分到 64 处理器、每处理器拥有一块 2×2×2 砖时brick-to-pencil 只需在4 个处理器的小组内做集合通信图中 4 种颜色一般地可把 $P$ 个处理器分为 $P^{2/3}$ 个小组、每组 $P^{1/3}$ 个处理器各小组内部做集合通信。pencil-to-pencil 则分为 $P^{1/2}$ 个小组、每组 $P^{1/2}$ 个处理器如图示y⇒z情形中 8 处理器小组拥有前部yz平面并在组内转置。默认使用点对点通信kspace_modify collective yes可切换为分区集合通信模式此时代码检测不相交子组大小把单个 $P$ 大小的 communicator 划分为多个小 communicator 分别调用集合通信。文档引用了 Argonne 国家实验室 IBM Blue Gene/Q 上的测试大进程数下分区集合通信的 FFT 性能显著优于纯点对点或全体 $P$ 进程集合通信。相关 FFT 实现位于 KSPACE 包如 src/KSPACE/fft3d.cpp、src/KSPACE/fft3d_wrap.cppKOKKOS 包另有 src/KOKKOS/fft3d_kokkos.cpp。手册补充的 FFT 细节fftMPI允许每个网格维度是小素因子2、3、5的倍数支持任意处理器数brick/pencil 分解不一定对齐但两种通信的通信子组规模仍按 $O(P^{1/3})$ 与 $O(P^{1/2})$ 缩放为 1d FFT 效率转置操作同时重排 3d 数据使不同方向在内存中连续重排可并入 MPI 缓冲的 pack/unpack大体系、多进程下并行 FFT 的主导成本常是通信而非 1d FFT 计算后者按每方向网格点数 $N$ 计为 $N\log N$因为网格只能 2d 铅笔分解而原子数据短程力/能计算可高效 3d 分解。两种对策MPI OpenMP 混合并行package omp减少参与 FFT 通信的 MPI 进程数MPI 域内用 OpenMP 保持算力。部分计算并行效率略降可能小于 FFT 通信开销多分区 run_style verlet/split启动多分区multi-partition计算用整数 1:p 映射把p个原子子域映射到 1 个 FFT 网格子域非键合/键合力与能在较大分区计算PPPM kspace 计算并发地在较小专用分区执行。此外LAMMPS 还有基于 PPPM 的偶极、色散Lennard-Jones长程作用求解器可与点电荷库仑长程作用联用GridComm类在模拟域上叠加规则网格并按与处理器子域一致的方式分区提供 owned/ghost 网格点值的 forward/reverse 通信方法被 PPPMFFT 网格与 MSM从细到粗的网格级联库仑力算法使用也适用于连续场与粒子耦合的模型例如两温模型TTM中fix ttm/grid用GridComm在分布式网格上做有限差分热方程替代使用复制网格的fix ttm。六、第二层并行OPENMP 包的设计与数据竞争规避INTEL、KOKKOS、OPENMP 三个包均提供 OpenMP 线程并行主要把循环分布在本地数据上与 MPI 空间域分解构成正交的两层并行策略。手册以结构最简单的OPENMP 包为例讲解INTEL/KOKKOS 更复杂支持向量化、加速器与双精度/单精度/混合精度选择其写法示例见开发者手册中为 OPENMP 包编写 style一节。6.1 引入 OpenMP 的理由资源利用效率多核机器上 MPI 任务常受通信带宽限制节点间通信共享同一网络链路减少 MPI 任务、增加 OpenMP 线程可降低通信开销与内存争用。代价是节点内内存带宽需求更高——域分解天然利于 CPU 缓存局部性而多线程实现为规避数据竞争使用 per-thread 复制存储需要额外的规约reduction线程数少时效率最高PPPM 扩展性长程静电求解器在高节点数下受 3d FFT all-to-all 通信限制子域变小后计算/通信比恶化用 OpenMP 后求解器可在 MPI 任务子集上运行改善整体扩展并行粒度OpenMP 提供节点内按粒子/邻居的更细粒度并行负载均衡域分解隐式假设体系均匀非均匀体系中域分解的负载均衡手段也有局限。改用 OpenMP 后每个 MPI rank 的子域更大负载均衡通常更有效高性能计算实践混合 MPIOpenMP 常是大型 HPC 集群/超算获得最优性能的关键。6.2 包的设计继承 最小修改 辅助类OPENMP 包的关键设计决策是保持源码改动最小便于维护并与非线程版标准实现保持同步继承OPENMP 版本从标准版本派生如PairLJCutOMP继承自PairLJCut只覆写多线程化或需修改的方法沿用 OPT 包的做法。仓库中 src/OPENMP/ 目录包含上百个此类派生类如 src/OPENMP/pair_lj_cut_omp.h 等最小修改修改后的代码结构与标准版尽量保持一致便于逐行对照辅助类多线程支持功能拆分到三个类——ThrOMP提供_thr 变体的 tally 等多线程感知功能多重继承见 src/OPENMP/thr_omp.hThrData管理 per-thread 数据结构避免 false sharing 拖慢见 src/OPENMP/thr_data.hFixOMP管理由package omp命令激活的全局多线程状态与设置见 src/OPENMP/fix_omp.h。6.3 规避数据竞争的三条路线与 LAMMPS 的选择MD 代码线程并行的核心问题是更新力、能、应力、力矩等 per-atom 累积量时的数据竞争交互计算总是涉及多个原子多线程会同时更新同一原子的数据。不做大幅重写的前提下有三条路线锁或原子操作保证同一时刻只有一个线程更新适合冲突很少的场景每个拥有元组中原子的线程都计算全部 per-atom 属性只存 owned 原子的数据线程多时最有效GPU 加速中常用per-thread 数据副本各线程独立更新自己的副本最后用 reduction 合并回常规 per-atom 存储需要额外的初始化/清零/合并步骤。OPENMP 包选择复制的 per-thread 数据结构因为单线程时性能无损区别于另两方案代码可维护且与非线程版相似且对 2–8 线程这种 OpenMP 的典型场景效率最高。6.4 循环调度与邻居表并行循环调度多线程并行通过把外层循环静态分配到各线程实现通常是处理邻居表i,j对或 bond/angle/dihedral/improper 列表时对本地原子i的循环。由于均匀体系中每原子邻居数相近线程间显著负载不均衡不常见。邻居表构建并行化除力计算外邻居表生成本身也被并行化——每个线程处理一块不同的 owned 原子并用MyPage页分配器实例管理自己的一组邻居页。七、实践要点小结场景推荐配置均匀密度体系默认comm_style brick 默认处理器网格密度不均真空、界面、沉积processors手动网格 balance平移或 RCB comm_style tiled大截断 / 刚体 / 扭结依赖多轮 ghost 交换机制必要时显式设置 ghost 截断大节点数下 PPPM 成为瓶颈package omp混合并行或 multi-partition run_style verlet/split专用 FFT 分区小/稀疏体系、兜底neighbor nsq多子势pair_style hybrid自动生成子邻居表从源码结构看上述机制的落点分别是分解与通信在 src/comm_brick.cpp 和 src/comm_tiled.cpp进程网格在 src/domain.cpp负载均衡在 src/balance.cpp邻居表与页分配在 src/neigh_list.cpp 与 src/my_page.cpp并行 FFT 在 src/KSPACE/fft3d.cpp 系列文件OpenMP 层在 src/OPENMP/ 目录。结合这些文件阅读本文各节可以完整掌握 LAMMPS 并行算法从域分解、原子通信、邻居表到长程作用与多线程实现的全链路原理。赞分享科研科学计算高性能计算【免费下载链接】lammpsPublic development project of the LAMMPS MD software package项目地址https://gitcode.com/gh_mirrors/la/lammps点击查看免费下载相关推荐LAMMPS 并行通信机制详解Ghost 原子、forward/reverse 通信与 MPI 实现LAMMPS 并行通信机制详解Ghost 原子、forward/reverse 通信与 MPI 实现 LAMMPSLarge scale Atomic/Mo科研科学计算高性能计算LAMMPS 长程相互作用与并行 FFT 实现PPPM 网格分解与通信机制深度解析LAMMPS 长程相互作用与并行 FFT 实现PPPM 网格分解与通信机制深度解析 导读 对于带电体系LAMMPS 通过基于 FFT 的 particle科研科学计算高性能计算解锁多语言输入困境Fcitx5安卓输入法框架的完整解决方案解锁多语言输入困境Fcitx5安卓输入法框架的完整解决方案 有没有想过为什么在移动设备上切换不同语言的输入体验总是如此割裂为什么我们不得不为每种语言安装独移动开发工具创作声明:本文部分内容由AI辅助生成(AIGC),仅供参考
阅读完成 · 觉得有帮助?
咨询建站