三木的创意空间 - 3D交互体验展示

数万粒子同屏运动，帧率稳定在 60 FPS，运行在浏览器里，不依赖任何原生插件。这是 jeantimex 开源的 WebGPU 流体仿真项目 fluid 的运行状态。

代码库里 TypeScript 占 71%，WGSL 占 28%。物理计算全部跑在 Compute Shader 里，没有任何 CPU 端的粒子循环。

粒子邻域搜索：8 Pass GPU 排序管线

SPH 的核心操作是每帧为每个粒子找到邻居粒子。朴素做法两层循环 O(n²)，几万粒子的量级在 CPU 上根本维持不了实时帧率。

fluid 用一套 8 Pass 的 GPU 排序管线解决这个问题，把粒子按空间位置排进连续内存，邻域查询退化成读几个连续格子。8 个 Pass 依次是：hash → clearOffsets → countOffsets → prefixScan → prefixCombine → scatter → reorder → copyBack。

第一步：给每个粒子算哈希 key。

fn hashCell3D(cellX: i32, cellY: i32, cellZ: i32) -> u32 {
    let blockSize = 50u;
    let ucell = vec3<u32>(
        u32(cellX + i32(blockSize / 2u)),
        u32(cellY + i32(blockSize / 2u)),
        u32(cellZ + i32(blockSize / 2u))
    );
    let localCell = ucell % blockSize;
    let blockID   = ucell / blockSize;
    let blockHash = blockID.x * 15823u
                  + blockID.y * 9737333u
                  + blockID.z * 440817757u;
    return localCell.x
         + blockSize * (localCell.y + blockSize * localCell.z)
         + blockHash;
}

这里用的是 Unity 的 block-based 空间哈希：把 3D 空间切成 50³ 的块，块内用线性下标（最多 125000），块之间用大素数偏移。相比纯素数哈希，同一块内的格子 key 值连续，排序后内存局部性更好，邻域查询时 cache miss 更少。

第二步：Blelloch 并行前缀和，把 count 数组变成 offset 数组。

这步是整条管线的核心。输入是每个 bucket 里有多少粒子（计数），输出是每个 bucket 的起始下标（偏移），只有这样才能把粒子 scatter 到正确位置。

Blelloch 算法分两个相：

// UP-SWEEP：从叶往根累加，根节点存总和
for (var d = 256u; d > 0u; d = d >> 1u) {
    workgroupBarrier();
    if (tid < d) {
        let ai = offset * (2u * tid + 1u) - 1u;
        let bi = offset * (2u * tid + 2u) - 1u;
        temp[bi] = temp[bi] + temp[ai];   // 右孩子 = 左 + 右
    }
    offset = offset * 2u;
}
 
// 根节点清零，开始 DOWN-SWEEP
if (tid == 0u) { temp[511u] = 0u; }
 
// DOWN-SWEEP：从根往叶分发前缀值
for (var d = 1u; d < 512u; d = d * 2u) {
    offset = offset >> 1u;
    workgroupBarrier();
    if (tid < d) {
        let ai = offset * (2u * tid + 1u) - 1u;
        let bi = offset * (2u * tid + 2u) - 1u;
        let t    = temp[ai];
        temp[ai] = temp[bi];          // 左孩子 = 父
        temp[bi] = temp[bi] + t;      // 右孩子 = 父 + 旧左
    }
}

每个 workgroup 256 个线程处理 512 个元素，全程在 shared memory 里跑，不碰全局内存。超过 512 的数组用 L0/L1/L2 三级层次处理，最大支持 134M 元素。

P2G 传递：原子操作处理写冲突

FLIP 方法需要把粒子速度"喷"到网格上（P2G，Particle-to-Grid），再从网格读回来（G2P）。问题是多个粒子可能同时写同一个格子，Compute Shader 里没有锁，需要用原子操作。

fluid 的做法是先把浮点速度乘以一个 SCALE 系数转成整数，再用 atomicAdd 累加：

@compute @workgroup_size(PARTICLE_WORKGROUP_SIZE)
fn transferToGrid(@builtin(global_invocation_id) id: vec3<u32>) {
    let pos = positions[pIdx].xyz;
    let vel = velocities[pIdx].xyz;
    let g   = worldToGrid(pos);
 
    for (var di = 0; di <= 1; di++) {
        for (var dj = 0; dj <= 1; dj++) {
            for (var dk = 0; dk <= 1; dk++) {
                // MAC 交错布局：三个速度分量各用不同的采样偏移
                let xPos = vec3<f32>(f32(cellX),       f32(cellY)+0.5, f32(cellZ)+0.5);
                let yPos = vec3<f32>(f32(cellX)+0.5,   f32(cellY),     f32(cellZ)+0.5);
                let zPos = vec3<f32>(f32(cellX)+0.5,   f32(cellY)+0.5, f32(cellZ)    );
 
                let xWeight = kernel(g - xPos);
                let yWeight = kernel(g - yPos);
                let zWeight = kernel(g - zPos);
 
                // 固定小数点原子加，避免浮点竞争写
                atomicAdd(&gridVelAtomic[cellIdx].x, i32(vel.x * xWeight * SCALE));
                atomicAdd(&gridVelAtomic[cellIdx].y, i32(vel.y * yWeight * SCALE));
                atomicAdd(&gridVelAtomic[cellIdx].z, i32(vel.z * zWeight * SCALE));
            }
        }
    }
}

注意三个速度分量的采样位置不同：Vx 在 (i, j+0.5, k+0.5)，Vy 在 (i+0.5, j, k+0.5)，Vz 在 (i+0.5, j+0.5, k)。这是 MAC 交错网格（Staggered MAC Grid）的布局，速度定义在格子面的中心而不是格子中心，能避免经典的棋盘格不稳定问题。

G2P 与 PIC/FLIP 混合

G2P 阶段的核心逻辑只有三行：

let vFlip = velOld + (vGridNew - vGridOld);  // FLIP：旧粒子速度 + 网格增量
let vPic  = vGridNew;                         // PIC：直接用网格速度
let vNew  = mix(vPic, vFlip, uniforms.fluidity);

PIC 直接把网格速度赋给粒子，数值耗散大，但稳定。FLIP 只把网格速度的变化量叠加到粒子上，保留了更多动能，但积累误差后会出现噪声。fluidity 参数（通常 0.95–0.99）控制两者的混合比例：偏 FLIP 更真实，偏 PIC 更稳定。

压力求解：红黑 vs Jacobi

FLIP 每帧需要解一个压力泊松方程，保证流体不可压缩。fluid 同时实现了两种迭代法，差别只在更新顺序。

Jacobi 每次迭代用上一步的全部值，所有格子并行更新，天然适合 GPU 但收敛慢。红黑 Gauss-Seidel 把格子按棋盘格染色，先更新红格再更新黑格，更新红格时黑格已经是本轮最新值，每次迭代就能用到更多新信息，收敛速度约快 2 倍：

@compute @workgroup_size(8, 4, 4)
fn jacobiRed(@builtin(global_invocation_id) id: vec3<u32>) {
    // 奇偶判断：(x+y+z) 为偶数的是红格，否则跳过
    let parity = (id.x + id.y + id.z) % 2u;
    if (parity != 0u) { return; }
 
    let si  = scalarIdx(id.x, id.y, id.z);
    if (marker[si] == 0u) { return; }  // 跳过气体格
 
    let div = divergence[si];
    var pL = 0.0; var pR = 0.0;
    var pB = 0.0; var pT = 0.0;
    var pBk = 0.0; var pFr = 0.0;
 
    if (id.x > 0u)              { pL  = pressure[scalarIdx(id.x-1u, id.y, id.z)]; }
    if (id.x < uniforms.nx-1u) { pR  = pressure[scalarIdx(id.x+1u, id.y, id.z)]; }
    if (id.y > 0u)              { pB  = pressure[scalarIdx(id.x, id.y-1u, id.z)]; }
    if (id.y < uniforms.ny-1u) { pT  = pressure[scalarIdx(id.x, id.y+1u, id.z)]; }
    if (id.z > 0u)              { pBk = pressure[scalarIdx(id.x, id.y, id.z-1u)]; }
    if (id.z < uniforms.nz-1u) { pFr = pressure[scalarIdx(id.x, id.y, id.z+1u)]; }
 
    pressure[si] = (invDx2*(pL+pR) + invDy2*(pB+pT) + invDz2*(pBk+pFr) - div)
                 * uniforms.precomputeJacobi;
}

jacobiBlack 的代码结构完全相同，只把 parity != 0u 换成 parity != 1u。每帧交替调用这两个 Pass 50 次即完成压力迭代。

6 条渲染管线，共享同一套物理状态

物理求解之外，fluid 提供 6 种可视化路径，并且支持在不中断仿真的情况下热切换：

Billboard Particles：相机对齐矩形片，带动态阴影和速度颜色映射
Marching Cubes：每帧从密度场实时提取等值面，生成 Metaball 风格网格
Volumetric Raymarching：对密度场做光线步进，支持光衰减和折射
Screen-Space Fluid：多 Pass 管线，包含深度平滑、厚度计算和泡沫模拟
FLIP Fluid：专为 PIC/FLIP 求解器配套，支持动态边界响应
Unified Dashboard：统一入口，热切换上述任意渲染器，仿真状态不丢失

写在最后

fluid 是研究和演示性质的项目，尚无生产级的工程配套。运行需要 Chrome/Edge 113 以上或 Firefox Nightly，WebGPU 在 Safari 上仍处于实验阶段。

代码注释密度很高，hash.wgsl 里每个算法步骤都有配图解释，prefix_sum.wgsl 里 Blelloch 算法的两个阶段有完整的树形示意。作为 WebGPU Compute Shader 的学习材料，可读性在同类项目里算好的。对 GPU 流体仿真感兴趣、想看完整 WGSL 实现的，可以从 spatial_grid.ts 和 flip_simulation.wgsl 两个文件入手。

GitHub：https://github.com/jeantimex/fluid (opens in a new tab)

ditherwave kiln-render