Constant FLUID_ADVECT_GLSL 

pub const FLUID_ADVECT_GLSL: &str = r#"
#version 430 core

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

layout(std430, binding = 0) buffer VelocityX  { float vel_x[]; };
layout(std430, binding = 1) buffer VelocityY  { float vel_y[]; };
layout(std430, binding = 2) buffer VelocityXn { float vel_xn[]; };
layout(std430, binding = 3) buffer VelocityYn { float vel_yn[]; };
layout(std430, binding = 4) readonly buffer Density { float density[]; };

uniform int grid_w;
uniform int grid_h;
uniform float dt;
uniform float dissipation;

int idx(int x, int y) { return clamp(x, 0, grid_w-1) + clamp(y, 0, grid_h-1) * grid_w; }

float sample_x(float px, float py) {
    int x0 = int(floor(px)); int x1 = x0 + 1;
    int y0 = int(floor(py)); int y1 = y0 + 1;
    float tx = fract(px); float ty = fract(py);
    return mix(mix(vel_x[idx(x0,y0)], vel_x[idx(x1,y0)], tx),
               mix(vel_x[idx(x0,y1)], vel_x[idx(x1,y1)], tx), ty);
}

float sample_y(float px, float py) {
    int x0 = int(floor(px)); int x1 = x0 + 1;
    int y0 = int(floor(py)); int y1 = y0 + 1;
    float tx = fract(px); float ty = fract(py);
    return mix(mix(vel_y[idx(x0,y0)], vel_y[idx(x1,y0)], tx),
               mix(vel_y[idx(x0,y1)], vel_y[idx(x1,y1)], tx), ty);
}

void main() {
    int x = int(gl_GlobalInvocationID.x);
    int y = int(gl_GlobalInvocationID.y);
    if (x >= grid_w || y >= grid_h) return;

    int i = idx(x, y);
    float vx = vel_x[i];
    float vy = vel_y[i];

    // Backtrace
    float px = float(x) - vx * dt;
    float py = float(y) - vy * dt;

    vel_xn[i] = sample_x(px, py) * dissipation;
    vel_yn[i] = sample_y(px, py) * dissipation;
}
"#;

Fluid simulation velocity advection compute pass.