优化基于C++的2D空间粒子运动蒙特卡洛模拟程序

问题背景 #

我正在开发一款采用Monte Carlo方法模拟2D空间中粒子运动的C++程序，规则如下：

• 每个时间步长内，每个粒子有概率p移动到当前位置半径r范围内的随机新位置；
• 若粒子移动到已有粒子的位置，两者发生碰撞并移动到各自原位置半径s范围内的随机新位置；
• 粒子到达2D空间边缘时，需反弹并沿相反方向继续运动。

现有代码在粒子数量较少时可正常运行，但当粒子数达到1000级规模时，程序运行极慢甚至崩溃，代码如下：

void move_particles(vector<Particle> &particles, double p, double r, double s) {
  for (int i = 0; i < particles.size(); i++) {
    if (should_move(p)) {
      // choose new location within radius r
      double new_x = particles[i].x + r * (2 * rand() - 1);
      double new_y = particles[i].y + r * (2 * rand() - 1);

      // handle boundary conditions
      if (new_x < 0) {
        new_x = -new_x;
      }
      if (new_x > WIDTH) {
        new_x = 2 * WIDTH - new_x;
      }
      if (new_y < 0) {
        new_y = -new_y;
      }
      if (new_y > HEIGHT) {
        new_y = 2 * HEIGHT - new_y;
      }

      // check for collisions with other particles
      for (int j = 0; j < particles.size(); j++) {
        if (i == j) continue;
        if (distance(new_x, new_y, particles[j].x, particles[j].y) < 2 * s) {
          // choose new location within radius s
          new_x = particles[j].x + s * (2 * rand() - 1);
          new_y = particles[j].y + s * (2 * rand() - 1);
        }
      }

      particles[i].x = new_x;
      particles[i].y = new_y;
    }
  }
}

性能瓶颈与崩溃原因分析 #

1. O(n²)碰撞检测复杂度：每个粒子都遍历所有其他粒子做碰撞检查，n=1000时单次时间步的碰撞检查次数达1e6次，粒子数增长时计算量呈平方级爆炸，直接导致运行卡顿。
2. 随机位置生成不符合规则：当前用r*(2*rand()-1)生成的是正方形区域内的随机点，而非规则要求的圆形区域，既破坏了模拟准确性，还可能让粒子移动到超出预期的范围，增加不必要的边界处理和碰撞概率。
3. 碰撞处理逻辑错误：碰撞时仅修改移动粒子的位置，完全未处理被碰撞粒子的位置，且多次碰撞会重复覆盖新位置，导致粒子位置异常，极端情况下引发崩溃。
4. rand()函数低效且质量差：rand()是低效的随机数生成器，大量调用会拖慢速度，且不支持线程安全，不利于后续并行优化。
5. 边界处理冗余：多个if判断可以简化，减少分支预测开销。

优化与修正方案 #

1. 空间分区优化碰撞检测（降复杂度至O(n)） #

将2D空间划分为边长为max(r, 2*s)的网格，每个粒子仅需检查所在网格及相邻8个网格内的粒子，无需遍历所有粒子。示例实现思路：

// 定义网格哈希函数
struct PairHash {
    template <class T1, class T2>
    std::size_t operator () (const std::pair<T1,T2> &p) const {
        auto h1 = std::hash<T1>{}(p.first);
        auto h2 = std::hash<T2>{}(p.second);
        return h1 ^ (h2 << 1);
    }
};

void move_particles(vector<Particle> &particles, double p, double r, double s) {
    const double GRID_SIZE = std::max(r, 2*s);
    // 构建粒子到网格的映射
    std::unordered_map<std::pair<int, int>, std::vector<int>, PairHash> grid;
    for (int i = 0; i < particles.size(); ++i) {
        int grid_x = static_cast<int>(particles[i].x / GRID_SIZE);
        int grid_y = static_cast<int>(particles[i].y / GRID_SIZE);
        grid[{grid_x, grid_y}].push_back(i);
    }

    // 其余逻辑...
    // 碰撞检查时，遍历当前网格及相邻网格：
    int new_grid_x = static_cast<int>(new_x / GRID_SIZE);
    int new_grid_y = static_cast<int>(new_y / GRID_SIZE);
    for (int dx = -1; dx <= 1; ++dx) {
        for (int dy = -1; dy <= 1; ++dy) {
            auto it = grid.find({new_grid_x + dx, new_grid_y + dy});
            if (it == grid.end()) continue;
            for (int j : it->second) {
                if (i == j) continue;
                // 用距离平方代替距离，避免开根号
                double dx_dist = new_x - particles[j].x;
                double dy_dist = new_y - particles[j].y;
                if (dx_dist*dx_dist + dy_dist*dy_dist < (2*s)*(2*s)) {
                    // 碰撞处理逻辑
                }
            }
        }
    }
}

2. 修正随机位置生成（圆形区域） #

用极坐标生成圆形区域内的均匀随机点：

// 初始化随机数生成器（全局或类内初始化一次）
std::mt19937 rng(std::random_device{}());
std::uniform_real_distribution<double> dist_theta(0.0, 2*M_PI);
std::uniform_real_distribution<double> dist_r(0.0, 1.0);

// 生成半径r内的随机点
double theta = dist_theta(rng);
double rand_r = r * std::sqrt(dist_r(rng)); // 确保均匀分布
double new_x = particles[i].x + rand_r * std::cos(theta);
double new_y = particles[i].y + rand_r * std::sin(theta);

3. 修复碰撞处理逻辑 #

碰撞时需同时更新两个粒子的位置，且避免重复处理同一对碰撞：

// 保存所有粒子的原位置
std::vector<std::pair<double, double>> original_pos(particles.size());
for (int i = 0; i < particles.size(); ++i) {
    original_pos[i] = {particles[i].x, particles[i].y};
}
// 存储最终位置
std::vector<std::pair<double, double>> final_pos = original_pos;
// 标记已处理碰撞的粒子
std::vector<bool> processed(particles.size(), false);

for (int i = 0; i < particles.size(); ++i) {
    if (processed[i]) continue;
    if (!should_move(p)) continue;

    // 生成新位置（圆形区域+边界处理）
    double new_x, new_y;
    // ...生成逻辑...

    // 碰撞检查（空间分区方式）
    int collision_j = -1;
    // ...查找碰撞的j...

    if (collision_j != -1 && !processed[collision_j]) {
        // 生成i的碰撞后位置（基于原位置）
        double theta_i = dist_theta(rng);
        double rand_r_i = s * std::sqrt(dist_r(rng));
        final_pos[i] = {original_pos[i].first + rand_r_i * std::cos(theta_i),
                        original_pos[i].second + rand_r_i * std::sin(theta_i)};
        // 生成j的碰撞后位置（基于原位置）
        double theta_j = dist_theta(rng);
        double rand_r_j = s * std::sqrt(dist_r(rng));
        final_pos[collision_j] = {original_pos[collision_j].first + rand_r_j * std::cos(theta_j),
                                original_pos[collision_j].second + rand_r_j * std::sin(theta_j)};
        processed[i] = true;
        processed[collision_j] = true;
    } else {
        final_pos[i] = {new_x, new_y};
    }
}

// 统一更新粒子位置
for (int i = 0; i < particles.size(); ++i) {
    particles[i].x = final_pos[i].first;
    particles[i].y = final_pos[i].second;
}

4. 简化边界处理 #

用对称性简化反弹计算：

// x方向反弹
new_x = std::abs(new_x);
new_x = WIDTH - std::abs(new_x - WIDTH);
// y方向反弹
new_y = std::abs(new_y);
new_y = HEIGHT - std::abs(new_y - HEIGHT);

5. 其他优化点 #

• 避免开根号：碰撞检测时用距离平方代替距离，减少浮点运算开销；
• 并行化：在粒子遍历阶段用OpenMP并行处理（需确保碰撞处理逻辑无竞态）；
• 减少拷贝：仅保存粒子的坐标而非整个对象，降低内存开销。

内容的提问来源于stack exchange，提问作者user3607022