Constant SPH_FORCE_WGSL 

pub const SPH_FORCE_WGSL: &str = r#"
pub(super) struct SphForceParams {
    n_particles: u32,
    h: f32,
    mu: f32,
    _pad: u32,
}

@group(0) @binding(0) var<storage, read>       positions:  array<vec4<f32>>;
@group(0) @binding(1) var<storage, read>       velocities: array<vec4<f32>>;
@group(0) @binding(2) var<storage, read>       densities:  array<f32>;
@group(0) @binding(3) var<storage, read>       pressures:  array<f32>;
@group(0) @binding(4) var<storage, read>       masses:     array<f32>;
@group(0) @binding(5) var<storage, read_write> forces:     array<vec4<f32>>;
@group(0) @binding(6) var<uniform>             params:     SphForceParams;

@compute @workgroup_size(64)
pub(super) fn main(@builtin(global_invocation_id) id: vec3<u32>) {
    let i = id.x;
    if (i >= params.n_particles) { return; }

    let pi = positions[i].xyz;
    let vi = velocities[i].xyz;
    var force: vec3<f32> = vec3<f32>(0.0, 0.0, 0.0);
    let h2 = params.h * params.h;

    for (var j: u32 = 0u; j < params.n_particles; j++) {
        if (j == i) { continue; }
        let r_vec = pi - positions[j].xyz;
        let r2 = dot(r_vec, r_vec);
        let r = sqrt(r2);
        if (r < 0.0001 || r >= params.h) { continue; }

        // Spiky gradient
        let h_r = params.h - r;
        let grad_mag = -45.0 / (3.14159265 * pow(params.h, 6.0)) * h_r * h_r;
        let grad = r_vec / r * grad_mag;

        // Pressure force
        let p_term = -masses[j] * (pressures[i] + pressures[j]) / (2.0 * densities[j]);
        force += p_term * grad;

        // Viscosity
        let v_lap = 45.0 / (3.14159265 * pow(params.h, 6.0)) * h_r;
        force += params.mu * masses[j] * (velocities[j].xyz - vi) / densities[j] * v_lap;
    }

    forces[i] = vec4<f32>(force, 0.0);
}
"#;

WGSL compute shader for SPH pressure force computation.