HLSL中利用内存连续性优化SPH流体模拟性能的疑问

我在Unity中使用并行HLSL后端开发3D平滑粒子流体动力学（SPH）模拟，这是一种拉格朗日流体模拟方法，通过粒子离散连续场并利用核函数平滑属性。目前功能正常，正进行性能优化以提升支持的粒子数量，已实现某论文中的高效空间分区算法，将同网格单元的粒子在Position等缓冲区中连续存储，具体步骤如下：

• 将粒子位置离散到网格单元（如3D下的(1,1,1)）；
• 用论文中的哈希函数将网格单元映射到表索引（如(1,1,1)→42）；
• 重排Particles缓冲区，使同单元的粒子连续存储；
• 邻域搜索时仅遍历目标粒子所在单元的27个相邻单元对应的缓冲区连续区间（单元大小等于核半径）；
• 每个单元的startIndex和endIndex存储在Offsets缓冲区中，startIndex=Offsets[hash]，endIndex=Offsets[hash+1]。

相关缓冲区（Position、Velocity、Density等）已按此规则排序，邻域搜索代码如下：

void CalculatePCISPHComponents(uint i, out float3 viscosity, out float3 surfaceTension, out float3 pressure, out float3 xsph) {
    viscosity = 0;
    surfaceTension = 0;
    xsph = 0;

    float3 posI = Positions[i];
    int3 gridPosI = GetGridPos(posI);

    float3 velI = Velocities[i];

    for (int x = -1; x < 2; x++) {
        for (int y = -1; y < 2; y++) {
            for (int z = -1; z < 2; z++) {
                int3 gridPosJ = gridPosI + int3(x, y, z);
                if (!IsInBounds(gridPosJ)) continue;

                uint hash = CalculateHashFromGrid(gridPosJ);

                uint startIndex = Offsets[hash];
                uint endIndex = Offsets[hash + 1];

                for (uint j = startIndex; j < endIndex; j++) {
                    if (i == j) continue;

                    float3 posOffset = posI - Positions[j];
                    float rSq = dot(posOffset, posOffset);

                    if (rSq < Epsilon * Epsilon || rSq > 4 * smoothingRadius * smoothingRadius) continue;

                    float r = length(posOffset);
                    float3 velOffset = velI - Velocities[j];
                    
                    float kernel = CubicSplineKernel(r);
                    float3 grad = CubicSplineGrad(posOffset, r);

                    viscosity += CalculateViscosityContribution(posOffset, velOffset, grad, i, j);
                    surfaceTension += -surfaceTensionMultiplier * kernel * (posOffset / r);
                    pressure += CalculatePressureContribution(posOffset, grad, i, j);

                    float massOverDensity = particleMass / Densities[j];
                    xsph += velocitySmoothing * massOverDensity * kernel * -velOffset;
                }
            }
        }
    }
}

该函数以单线程单粒子的方式调用，调度代码如下：

[numthreads(NumThreads, 1, 1)]
void CalculateNonPressureAcceleration(uint3 id : SV_DispatchThreadID) {
    if (id.x >= instanceCount) return;

    IntermediateAccelerations[id.x] = CalculatePCISPHAcceleration(id.x);
}

尽管遍历的是内存连续区域，但粒子数超过约15K后性能急剧下降，而同类实现可支持数十万粒子。我怀疑存在GPU内存访问的理解误区，已尝试以下方案：

• 基于单元的调度模式，先将邻域粒子复制到共享内存再调度线程；
• 更换哈希函数（尝试过Morton编码和3D线性化）；
• 帧捕获分析显示，float3 posOffset = posI - Positions[j];对应的编译代码占某内核16.5%的计算时间，推测内存访问是瓶颈。

运行环境：M4 Metal 10核GPU，每组256线程，8192个bin。希望获得优化思路或排查方向。

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