基于F=V.B的粒子偏转模拟：2D速度网格高效构建问询

你的核心问题在于Python层面的双重循环+布尔索引筛选，这在4亿多粒子的规模下效率极低（时间复杂度O(N*M)，N是粒子数，M是网格数）。下面给出两种基于NumPy向量化操作的优化方法，速度能提升几个数量级，同时保留原网格精度（不需要依赖降采样divider牺牲精度）。

方法一：用np.bincount实现快速统计（推荐） #

这个方法通过一次性计算所有粒子的网格索引，再用NumPy的高效统计函数bincount计算每个网格的动量总和与粒子数，最后得到平均速度。整个过程仅需几次数组遍历，时间复杂度为O(N)。

步骤说明： #

1. 计算每个粒子对应的网格索引：向量化计算，避免循环。
2. 处理边界索引：防止粒子坐标超出网格范围导致索引越界。
3. 统计每个网格的动量总和与粒子数：用bincount一次性完成。
4. 计算平均速度：总和除以粒子数再除以质子质量，最后reshape为2D数组。

完整代码： #

import numpy as np

# 原始参数
nx, ny = 1600, 900
xmin = -30.0e-6
xmax = 50.0e-6
ymin = -45.0e-6
ymax = 45.0e-6
m_p = 1.6726219e-27  # 质子质量

# 如果你需要降采样，调整这里的divider即可，不需要修改核心逻辑
divider = 1  # 用1就是原网格精度，完全不需要几天时间！
nx_new = nx // divider
ny_new = ny // divider

# 计算新网格的分辨率
res_x = (xmax - xmin) / nx_new
res_y = (ymax - ymin) / ny_new

# 1. 向量化计算每个粒子的网格索引（x和y方向）
# 注意：先减去xmin/ymin偏移，再除以分辨率，最后转成整数
grid_x_idx = ((grid_pro_x - xmin) / res_x).astype(int)
grid_y_idx = ((grid_pro_y - ymin) / res_y).astype(int)

# 2. 处理边界：将超出网格的粒子索引限制在0到nx_new-1/ny_new-1范围内
grid_x_idx = np.clip(grid_x_idx, 0, nx_new - 1)
grid_y_idx = np.clip(grid_y_idx, 0, ny_new - 1)

# 3. 计算线性索引（将2D索引转为1D，方便bincount处理）
linear_idx = grid_y_idx * nx_new + grid_x_idx

# 4. 统计每个网格的Px总和与粒子数
px_sum = np.bincount(linear_idx, weights=Px_pro, minlength=nx_new * ny_new)
particle_counts = np.bincount(linear_idx, minlength=nx_new * ny_new)

# 5. 处理粒子数为0的网格（避免除以0），这里设为0，你也可以根据需求调整
particle_counts[particle_counts == 0] = 1

# 6. 计算Vx并reshape为2D数组
Vx = (px_sum / particle_counts / m_p).reshape((ny_new, nx_new))

方法二：用np.histogram2d直接计算均值 #

NumPy的histogram2d也可以通过设置weights参数直接计算加权均值，代码更简洁：

import numpy as np

# 同样的参数设置...
divider = 1
nx_new = nx // divider
ny_new = ny // divider
res_x = (xmax - xmin)/nx_new
res_y = (ymax - ymin)/ny_new

# 生成网格的边界值
x_edges = np.linspace(xmin, xmax, nx_new + 1)
y_edges = np.linspace(ymin, ymax, ny_new + 1)

# 计算每个网格的Px总和（加权直方图）和粒子数（普通直方图）
px_sum, _, _ = np.histogram2d(grid_pro_x, grid_pro_y, bins=(x_edges, y_edges), weights=Px_pro)
particle_counts, _, _ = np.histogram2d(grid_pro_x, grid_pro_y, bins=(x_edges, y_edges))

# 处理0计数，计算Vx
particle_counts[particle_counts == 0] = 1
Vx = (px_sum / particle_counts / m_p).T  # 注意这里需要转置，因为histogram2d的输出是(nx_new, ny_new)，而你需要(ny_new, nx_new)

原代码的关键问题 #

除了效率低，你的原代码还有一个逻辑错误：
你计算网格边界时直接用xx = i*res_x、yy = j*res_y，完全忽略了xmin和ymin的偏移量！实际粒子的x坐标从-30e-6开始，而你的xx从0开始，导致筛选条件grid_pro_x >= xx几乎无法匹配到任何粒子，得到的结果完全错误。正确的边界计算应该是xx_min = xmin + j*res_x（j对应x方向）、yy_min = ymin + i*res_y（i对应y方向）。

性能对比 #

用原代码循环144万次（divider=1）需要数天，而用上述向量化方法，即使处理4亿多粒子，也只需要几分钟到几十分钟（取决于你的硬件），完全不需要依赖降采样牺牲精度。

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