如何将轨道模拟的串行for循环改造为并行运行模式

我编写了一个模拟程序，用于计算流浪行星闯入太阳系后地球轨道的变化情况，当前串行版本的核心代码如下：

import numpy as np 

total_time = 86400*30
dt = 720
position_Rogue = np.zeros((int(total_time / dt), 3))
velocity_Rogue = np.zeros((int(total_time / dt), 3))
acceleration_Rogue = np.zeros((int(total_time / dt), 3))
position_earth = np.zeros((int(total_time / dt), 3))
velocity_earth = np.zeros((int(total_time / dt), 3))
acceleration_earth = np.zeros((int(total_time / dt), 3))

position_earth[0] = np.array([-8.720402184155567E+07 , -1.249081082441206E+08 , 3.737986027865112E+04]) 
velocity_earth[0] = np.array([2.401992974509115E+01 , -1.707462009431704E+01 ,1.433502989955038E-03]) 

position_Rogue[0] = np.array([9.720402184155567E+07 , -1.549081082441206E+08 , 3.737986027865112E+04]) 
velocity_Rogue[0] = np.array([3.401992974509115E+01 , 1.707462009431704E+01 ,2.433502989955038E-03])
 
G = 6.674e-20
M_SUN= float(1.989 * 10**30)
M_earth = float(5.972 * 10**24)
M_Rogue = float(M_earth*300)#increased from m_earth to30*m_earth .... till 300

for n in range(1, int(total_time/dt)):
    
        acceleration_Rogue[n] = (- G * M_SUN* position_Rogue[n-1]/np.linalg.norm(position_Rogue[n-1])**3- G * M_earth * (position_Rogue[n-1]-position_earth[n-1])/np.linalg.norm(position_Rogue[n-1]-position_earth[n-1]**3)
        
        )
        velocity_Rogue[n] = velocity_Rogue[n-1] + acceleration_Rogue[n-1]*dt
        position_Rogue[n] = position_Rogue[n-1] + velocity_Rogue[n-1]*dt + 0.5*dt**2*acceleration_Rogue[n-1]

        
        acceleration_earth[n] = (
        - G * M_SUN * position_earth[n-1]/np.linalg.norm(position_earth[n-1])**3   
         
        - G * M_Rogue* (position_earth[n-1] - position_Rogue[n-1])/np.linalg.norm(position_earth[n-1] - position_Rogue[n-1])**3
            )

当前循环为串行运行，每次迭代计算天体新位置时依赖前一次的数据。我希望将其改造为并行运行以适配超级计算机，目标是模拟总时长约100年、时间步长为1秒，以下是具体实现思路：

并行化核心思路 #

1. 时间分块+MPI节点间并行 #

由于模拟总时长100年对应约3.15e9个时间步，直接串行完全不可行。可以将整个时间轴拆分为多个连续的块，每个MPI进程负责处理一个块的计算：

• 每个进程只存储自己负责块的位置、速度、加速度数组，无需加载整个时间序列，大幅降低内存占用
• 块与块之间仅需传递边界数据：前一个块的最后一步位置、速度作为当前块的初始条件
• 适合超级计算机的多节点分布式架构，并行度等于分块数量（可根据节点数调整）

2. 单时间步内的SIMD/OpenMP并行 #

每个时间步内的加速度、速度、位置计算包含大量向量运算，可利用CPU的SIMD指令集或OpenMP实现节点内并行：

• 对地球和流浪行星的加速度计算可同时启动线程并行处理
• 将向量运算（如np.linalg.norm、数组乘法）拆解为单个元素的操作，用OpenMP对元素级操作并行化
• 结合NumPy的向量化特性，避免显式循环，让底层自动利用SIMD并行

3. GPU加速并行 #

利用GPU的大规模并行处理能力，将每个时间步的计算转移到GPU执行：

• 用CuPy替代NumPy，自动将数组运算映射到GPU线程，实现单时间步内的向量级并行
• 对于循环迭代，可利用CUDA核函数手动实现时间步的并行计算（需注意迭代依赖，仅能并行单步内的独立运算）
• 借助专业天体模拟GPU库（如Rebound-GPU），直接调用优化后的并行积分器，无需手动实现循环

4. MPI+OpenMP混合并行 #

结合以上两种方式，充分利用超级计算机的节点间和节点内资源：

• MPI负责跨节点的时间分块并行
• 每个节点内部用OpenMP对单时间步内的向量运算做并行加速
• 这种模式能最大化超算的硬件利用率，适合大规模长时间尺度的模拟

关键注意事项 #

• 固定1秒时间步会导致总步数极多，建议改用自适应时间步长：当天体运动平缓时增大步长，接近相遇时缩小步长，减少计算量的同时保证精度
• 显式积分的迭代依赖无法完全消除，因此无法实现全时间步的完全并行，上述思路都是在依赖约束下最大化并行度
• 优先考虑使用成熟的并行天体模拟库（如Rebound、AMUSE），这些库已针对超算做了深度优化，比手动实现并行效率更高

内容的提问来源于stack exchange，提问作者M.M. CAN