Constant KERNEL_FLUID_DIFFUSE 

pub const KERNEL_FLUID_DIFFUSE: &str = r#"
// Jacobi iteration for 2D fluid diffusion
// Reads from one grid, writes to another (ping-pong).
// d(x)/dt = k * laplacian(x)
// Jacobi: x_new[i,j] = (x_old[i,j] + alpha * (x[i-1,j] + x[i+1,j] + x[i,j-1] + x[i,j+1])) / (1 + 4*alpha)
// where alpha = k * dt / (dx * dx)

layout(local_size_x = 16, local_size_y = 16) in;

layout(std430, binding = 0) readonly buffer GridIn {
    float grid_in[];
};

layout(std430, binding = 1) writeonly buffer GridOut {
    float grid_out[];
};

uniform uint u_width;
uniform uint u_height;
uniform float u_alpha;  // diffusion_rate * dt / (dx * dx)
uniform float u_r_beta;  // 1.0 / (1.0 + 4.0 * alpha)

uint idx2d(uint x, uint y) {
    return y * u_width + x;
}

void main() {
    uint x = gl_GlobalInvocationID.x;
    uint y = gl_GlobalInvocationID.y;

    if (x >= u_width || y >= u_height) return;

    // Boundary: clamp to edge
    uint x0 = max(x, 1u) - 1u;
    uint x1 = min(x + 1u, u_width - 1u);
    uint y0 = max(y, 1u) - 1u;
    uint y1 = min(y + 1u, u_height - 1u);

    float center = grid_in[idx2d(x, y)];
    float left   = grid_in[idx2d(x0, y)];
    float right  = grid_in[idx2d(x1, y)];
    float down   = grid_in[idx2d(x, y0)];
    float up     = grid_in[idx2d(x, y1)];

    float result = (center + u_alpha * (left + right + down + up)) * u_r_beta;
    grid_out[idx2d(x, y)] = result;
}
"#;

Fluid diffusion kernel using Jacobi iteration.