void Simulate() {

// CLEAR GRID
std::size_t grid_size = grid.size();

// 确保 grid_size 不超出 int 的范围

pragma omp parallel for

for (int i = 0; i < static_cast<int>(grid_size); ++i) {
    grid[i].vel = glm::vec3(0.0f);
    grid[i].mass = 0.0f;
}

// P2G_1
//质量与动量转移
 #pragma omp parallel for
for (int i = 0; i < particles.size(); i++) {
    auto &p = particles[i];
    glm::uvec3 cell_idx = glm::uvec3(p.pos);
    glm::vec3 cell_diff = (p.pos - glm::vec3(cell_idx)) - 0.5f;

    glm::vec3 weights_[3];
    weights_[0] = 0.5f * glm::pow(0.5f - cell_diff, glm::vec3(2.0f));
    weights_[1] = 0.75f - glm::pow(cell_diff, glm::vec3(2.0f));
    weights_[2] = 0.5f * glm::pow(0.5f + cell_diff, glm::vec3(2.0f));

    for (uint32_t gx = 0; gx < 3; ++gx) {
        for (uint32_t gy = 0; gy < 3; ++gy) {
            for (uint32_t gz = 0; gz < 3; ++gz) {
                float weight = weights_[gx].x * weights_[gy].y * weights_[gz].z;

                glm::uvec3 cell_pos = glm::uvec3(cell_idx.x + gx - 1,
                                                 cell_idx.y + gy - 1,
                                                 cell_idx.z + gz - 1);
                glm::vec3 cell_dist = (glm::vec3(cell_pos) - p.pos) + 0.5f;
                glm::vec3 Q = p.C * cell_dist;

                int cell_index = (int) cell_pos.x +
                                 (int) cell_pos.y * grid_res +
                                 (int) cell_pos.z * grid_res * grid_res;

                float mass_contrib = weight * particle_mass;
                grid[cell_index].mass += mass_contrib;
                grid[cell_index].vel += mass_contrib * (p.vel + Q);
            }

        }
    }
}

// P2G_2
//密度计算和压力计算
 #pragma omp parallel for
for (int i = 0; i < particles.size(); i++) {
    auto &p = particles[i];

    glm::uvec3 cell_idx = glm::uvec3(p.pos);
    glm::vec3 cell_diff = (p.pos - glm::vec3(cell_idx)) - 0.5f;

    glm::vec3 weights_[3];
    weights_[0] = 0.5f * glm::pow(0.5f - cell_diff, glm::vec3(2.0f));
    weights_[1] = 0.75f - glm::pow(cell_diff, glm::vec3(2.0f));
    weights_[2] = 0.5f * glm::pow(0.5f + cell_diff, glm::vec3(2.0f));

    float density = 0.0f;
    for (uint32_t gx = 0; gx < 3; ++gx) {
        for (uint32_t gy = 0; gy < 3; ++gy) {
            for (uint32_t gz = 0; gz < 3; ++gz) {
                float weight = weights_[gx].x * weights_[gy].y * weights_[gz].z;

                glm::uvec3 cell_pos = glm::uvec3(cell_idx.x + gx - 1,
                                                 cell_idx.y + gy - 1,
                                                 cell_idx.z + gz - 1);

                int cell_index = (int) cell_pos.x +
                                 (int) cell_pos.y * grid_res +
                                 (int) cell_pos.z * grid_res * grid_res;

                density += grid[cell_index].mass * weight;
            }
        }
    }

    float volume = particle_mass / density;
    float pressure = std::max(-0.1f, eos_stiffness *
                                     (std::pow(density / rest_density, eos_power) - 1.0f));

    glm::mat3 stress = glm::mat3(
            -pressure, 0, 0,
            0, -pressure, 0,
            0, 0, -pressure
    );

     glm::mat3 strain = p.C;

     //float trace = strain[0][0] + strain[1][0] + strain[2][0]; // DEBUG
     float trace = glm::determinant(strain);
     strain[0][0] = strain[1][0] = strain[2][0] = trace;

     glm::mat3 viscosity_term = dynamic_viscosity * strain;
     stress += viscosity_term;

    auto eq_16_term_0 = -volume * 4 * stress * dt;

    for (uint32_t gx = 0; gx < 3; ++gx) {
        for (uint32_t gy = 0; gy < 3; ++gy) {
            for (uint32_t gz = 0; gz < 3; ++gz) {
                float weight = weights_[gx].x * weights_[gy].y * weights_[gz].z;

                glm::uvec3 cell_pos = glm::uvec3(cell_idx.x + gx - 1,
                                                 cell_idx.y + gy - 1,
                                                 cell_idx.z + gz - 1);

                glm::vec3 cell_dist = (glm::vec3(cell_pos) - p.pos) + 0.5f;

                int cell_index = (int) cell_pos.x +
                                 (int) cell_pos.y * grid_res +
                                 (int) cell_pos.z * grid_res * grid_res;

                glm::vec3 momentum = (eq_16_term_0 * weight) * cell_dist;
                grid[cell_index].vel += momentum;
            }
        }
    }
}

// GRID UPDATE
// 获取 grid 的大小
grid_size = grid.size();

// GRID UPDATE

pragma omp parallel for

for (int i = 0; i < static_cast<int>(grid_size); ++i) {
    auto& cell = grid[static_cast<std::size_t>(i)]; // 使用 static_cast 将索引转换为 std::size_t
    if (cell.mass > 0) {
        cell.vel /= cell.mass;
        cell.vel += dt * glm::vec3(0.0f, gravity, 0.0f);

        int index = cell.index;
        int x = index / (grid_res * grid_res);
        index /= grid_res;
        int y = (index / grid_res) % grid_res;
        index /= grid_res;
        int z = index;

        if (x < 1 || x > grid_res - 2) { cell.vel.x = 0.0f; }
        if (y < 1 || y > grid_res - 2) { cell.vel.y = 0.0f; }
        if (z < 1 || z > grid_res - 2) { cell.vel.z = 0.0f; }
    }
}

// G2P
 #pragma omp parallel for
for (int i = 0; i < static_cast<int>(particles.size()); ++i) {
    auto& p = particles[static_cast<std::size_t>(i)]; // 使用索引访问粒子

    p.vel = glm::vec3(0.0f);

    glm::uvec3 cell_idx = glm::uvec3(p.pos);
    glm::vec3 cell_diff = (p.pos - glm::vec3(cell_idx)) - 0.5f;

    glm::vec3 weights_[3];

    weights_[0] = 0.5f * glm::pow(0.5f - cell_diff, glm::vec3(2.0f));
    weights_[1] = 0.75f - glm::pow(cell_diff, glm::vec3(2.0f));
    weights_[2] = 0.5f * glm::pow(0.5f + cell_diff, glm::vec3(2.0f));

    glm::mat3 B = glm::mat3(0.0f);
    for (uint32_t gx = 0; gx < 3; ++gx) {
        for (uint32_t gy = 0; gy < 3; ++gy) {
            for (uint32_t gz = 0; gz < 3; ++gz) {
                float weight = weights_[gx].x * weights_[gy].y * weights_[gz].z;
                // std::cout << weight << std::endl;
                glm::uvec3 cell_pos = glm::uvec3(cell_idx.x + gx - 1,
                    cell_idx.y + gy - 1,
                    cell_idx.z + gz - 1);

                glm::vec3 cell_dist = (glm::vec3(cell_pos) - p.pos) + 0.5f;

                int cell_index = (int)cell_pos.x +
                    (int)cell_pos.y * grid_res +
                    (int)cell_pos.z * grid_res * grid_res;

                glm::vec3 weighted_velocity = grid[cell_index].vel * weight;

                B += glm::mat3(weighted_velocity * cell_dist.x,
                    weighted_velocity * cell_dist.y,
                    weighted_velocity * cell_dist.z);

                p.vel += weighted_velocity;
            }
        }
    }

    p.C = B * 4.0f;
    p.vel *= damping;
    p.pos += p.vel * dt;
    p.pos = glm::clamp(p.pos, 1.0f, grid_res - 2.0f);

    glm::vec3 x_n = p.pos + p.vel;
    const float wall_min = 3.0f;
    const float wall_max = grid_res - 4.0f;
    if (x_n.x < wall_min) p.vel.x += (wall_min - x_n.x);
    if (x_n.x > wall_max) p.vel.x += (wall_max - x_n.x);
    if (x_n.y < wall_min) p.vel.y += (wall_min - x_n.y);
    if (x_n.y > wall_max) p.vel.y += (wall_max - x_n.y);
    if (x_n.z < wall_min) p.vel.z += (wall_min - x_n.z);
    if (x_n.z > wall_max) p.vel.z += (wall_max - x_n.z);
}

}这是使用openmp的进行了加速，但是我有在youtube上看到可以对其用simd指令优化到单线程的也能流畅运行的程度，但本人对此一窍不通（悲），所以来探讨一下是否有可行性。（本人尝试过使用simd进行改写，但是debug的太痛苦了，模板元什么的最讨厌了）https://github.com/Hab5/mls-mpm这是对应的完整代码仓库，https://www.youtube.com/watch?v=4Y58Pg9tpSo这是对应的YouTube视频介绍（梦开始的地方）

这是你的仓库吗？看到这个仓库使用了OpenGL计算着色器，为什么说是单CPU的，可以参考一下我的：https://github.com/archibate/opengl_mpm

P2G这一步里的+=需要是atomic的吧？不然累计结果会不对，导致仿真结果爆炸。

你可以用chrono测一测每一步的时间，你应该会发现P2G这一步花的时间是最长的，因为P2G涉及了scatter操作，会出现多个线程同时访问同一个grid的情况，需要低效的atomic。为了避免多线程竞争导致的性能损失，有两种方案：

先把粒子按照莫顿码排个序，排序后因为openmp会把值接近的i放在一个线程里执行，从而只有边缘存在多线程同时atomic访问同一个grid的情况。
如果有4个cpu，那就把粒子按位置分成4块，分别在各自的区域里，那就不需要atomic了。

看到你把omp都注释了，是没加atomic炸了然后debug了半天是吧

                    float mass_contrib = weight * particle_mass;
                    grid[cell_index].mass += mass_contrib;
                    grid[cell_index].vel += mass_contrib * (p.vel + Q);

需要改成：

                    float mass_contrib = weight * particle_mass;
#pragma omp atomic
                    grid[cell_index].mass += mass_contrib;
#pragma omp atomic
                    grid[cell_index].vel.x += mass_contrib * (p.vel.x + Q);
#pragma omp atomic
                    grid[cell_index].vel.y += mass_contrib * (p.vel.y + Q);
#pragma omp atomic
                    grid[cell_index].vel.z += mass_contrib * (p.vel.z + Q);

所以，并行并不是儿童玩具，并不是说加上#pragma omp parallel for就能免费提升性能。并行是有成本的，99%的情况下，无脑加#pragma omp parallel for要么崩溃，要么没有任何加速效果。你必须改变编程思路，适应并行编程的特点，而不是继续按照串行的思路来写代码，这是一套截然不同的思路，并不是免费的。CPU的并行编程和串行时有很大不同，再加上SIMD优化也会完全不同，GPU更是完全不同，在GPU上并行编程若想获得理想的加速效果，有时甚至需要加大计算量来提升并行度，用空间换时间等。所以，实际上很多算法，因为本身不适合GPU的并行编程方式，而没有办法用GPU并行或者收效甚微。好消息是，MPM还算是一个GPU可以加速的算法，这里的P2G俗称scatter操作，thrust里是有稳定的方案的，当然CPU也有相应方案，无论如何，这些并行的方案所写的代码都会和你现在截然不同。

串行写法

例如，在串行编程中过滤出一个vector所有的偶数：

std::vector<int> in = {1, 2, 3, 4};
std::vector<int> out;
for (size_t i = 0; i < in.size(); ++i) {
  if (in[i] % 2 == 0)
    out.push_back(in[i]);
}

很明显，如果你无脑加上#pragma omp parallel for，这个std::vector的push_back会直接奔溃，因为他并不是线程安全的。

串行思维的并行

坚持使用“串行思维”写并行的方法是使用并发安全的vector容器tbb::concurrent_vector，其push_back是线程安全的：

std::vector<int> in = {1, 2, 3, 4};
tbb::concurrent_vector<int> out;
#pragma omp parallel for
for (size_t i = 0; i < in.size(); ++i) {
  if (in[i] % 2 == 0)
    out.push_back(in[i]);
}

但是你会发现性能提升不大，甚至可能反而倒退，这就是错误用“串行思维”写并行代码的后果（你的代码中需要atomic地+=一系列随机位置的grid成员，就属于这种错误思想）。

真正的大并行

真正的并行方案可能令你大开眼界，需要把上面的单个for拆成三步：

std::vector<int> in = {1, 2, 3, 4};
std::vector<int> cnt(in.size());
#pragma omp parallel for
for (size_t i = 0; i < in.size(); ++i) {
  cnt[i] = in[i] % 2 == 0;
}
tbb::parallel_exclusive_scan(cnt.begin(), cnt.end());
std::vector<int> out(cnt.back() + in.back());
#pragma omp parallel for
for (size_t i = 0; i < in.size(); ++i) {
  if (in[i] % 2 == 0)
    out[cnt[i]] = in[i];
}

连续判断了两遍，很低效？不好意思。并行编程就是需要经常牺牲“计算量”来提升并行度，最终在高并行环境中反而更高效。

Building-Blocks

其中parallel_exclusive_scan是“前缀和”算法，是可以高度并行的，在GPU和CPU上均有高效的实现。此类常用的并行算法被称为“基本算法”，也是Thread-Building-Blocks的Building-Blocks得名由来。通常写并行算法时，实际上会先拆分成多个“基本算法”，然后就可以直接调用TBB里的基本块来实现。TBB就像乐高积木一样，你需要自己根据你想要拼出模型的形状，选择适合的积木块。 parallel_for，parallel_exclusive_scan，parallel_reduce，parallel_merge_sort，parallel_bitonic_sort，parallel_radix_sort，parallel_scatter，parallel_gather等都是这样的基本算法。而在你原始的头脑中，只有parallel_for这一种最无脑的基本算法，因为其他的算法你在串行编程中可能根本没见过。不去把自己的算法拆解成这些可以高度并行的“基本算法”，是你并行编程这么困难的根本原因。顺便一提，上面这种过滤出指定条件（例如“偶数”）元素的算法被称为 filter，虽然可以进一步拆成 for+scan+for，但因为非常常用，也被当作是一种基本算法了，thrust中就有 thrust::copy_if。

CPU优化可以走邪道

当然，上面这种基于Building-Blocks的“三步走”战略是并行得最彻底的，好处是：不论是CPU并行，CPU SIMD，GPU并行，多GPU并行，多机联网并行，都可以轻松移植过去。但是因为牺牲了一点“计算量”，如果CPU核心数量不多（例如只有2核），可能还不如串行来得快。因此一般CPU并行和CPU SIMD都会有单独的优化方案，例如上面这种过滤出偶数的，用x86 SIMD的AVX2指令集有一种基于_mm_maskstore_epi32的黑科技可以单核加速20倍，如果用AVX512又会有所不同，CPU并行也有黑科技可以加速，不必采用最泛用的“三步走”。但是如果你已经写出并行的最彻底的GPU方案（例如上面这种“三步走”），那么也是可以无缝降级到CPU并行和CPU SIMD的，只不过在CPU上会不如专为SIMD优化的高效。当然，GPU也有自己的邪道，例如利用shared memory在单个block内局部地三步走，而不必整个在全局范围内跑最通用的“三步走”浪费计算量和缓存局域性，但是如果你要方便移植的话基于Building-Blocks是最方便的。所以，你考虑好你的目标用户到底是GPU还是CPU。据我估算，你这种算法上GPU可以获得免费的20倍提升，专门优化后可以获得100倍。CPU如果无脑OpenMP（搞定atomic问题后）可以获得免费的2倍提升，但是深度优化（SIMD+并行）可以继续获得20倍提升。

好的，感谢小彭老师的指导，我使用#pragma omp parallel for能得到一点效率上的提升，debug的话主要是想用_mm_maskstore_epi32类似这种simd、avx指令来重写一下（虽然我也感觉需要一定的同步机制代码，可能编译器自带的锁机制满足了已经，没遇到崩溃的情况），但是现在我觉得走gpu加速会更好一点，所以是将这部分代码修改为在.cu代码文件，在cpu和gpu之间传递数据吗

parallel101 / simdtutor

小彭老师大大，有关mpm的流体模拟函数优化 #10

pragma omp parallel for

pragma omp parallel for

串行写法

串行思维的并行

真正的大并行

Building-Blocks

CPU优化可以走邪道