So I am trying to implement a fluid simulation, as described in Chapter 38 of GPU Gems.

I am also drawing heavily from the excellent implementation by Pavel.

Where I'm getting a bit stuck is with advection.

So as I understand it, advection is the process that moves quantities along with the movement of the fluid. So for example, if at one timestamp, you have a velocity v at position p, then in the next timestamp, you have velocity V at position p + v * dt, at least according to advection.

So as it says in GPU gems, there is a difficulty implementing this logic on the GPU:

We implement our simulation in fragment programs, which cannot change the locations of the fragments they are writing. This forward-integration method requires the ability to "move" the particles, so it cannot be implemented on current GPUs.

The solution is to invert the problem and use an implicit method (Stam 1999). Rather than advecting quantities by computing where a particle moves over the current time step, we trace the trajectory of the particle from each grid cell back in time to its former position, and we copy the quantities at that position to the starting grid cell.

So in other words, instead of projecting values forward in the grid, based on the velocity value in the current timestamp, we use the current velocity to trace backwards in the grid, to see what the next value should be which is coming in the direction of our current grid cell.

The problem with this approach appears to occur when velocity is zero. So for instance, consider that we have a 1 dimensional grid of velocities:


If we advect this grid according to the backwards lookup, at the next timestamp, we will have the following:


So this is working as expected.

However this breaks down, and if we have a zero in any position, there's no way to get another value to replace the zero:

t+0: [3|1|2|0]

t+1: [0|3|3|0]

t+2: [0|0|0|0]

So basically, if the value of the cell is zero, it will never look up any value other than itself, so there is no way for another quantity to advect into it.

And indeed this is exactly the behavior I am seeing in my simulation:

enter image description here

And here's the shader I am using:

layout (local_size_x = WORKGROUP_SIZE, local_size_y = WORKGROUP_SIZE, local_size_z = 1 ) in;

precision highp float;
precision highp sampler2D;

layout(push_constant) uniform PushConstant {
    vec2 viewport;
    vec2 texelSize;
    vec2 dyeTexelSize;
    float dt;
    float dissipation;
} pushConstants;

layout(binding = 0, rgba8) uniform writeonly image2D outputImage;

layout(binding = 1) uniform sampler2D uVelocity;
layout(binding = 2) uniform sampler2D uSource;

vec4 bilerp(sampler2D sam, vec2 uv, vec2 tsize) {
    vec2 st = uv / tsize - 0.5;
    vec2 iuv = floor(st);
    vec2 fuv = fract(st);
    vec4 a = texture(sam, (iuv + vec2(0.5, 0.5)) * tsize);
    vec4 b = texture(sam, (iuv + vec2(1.5, 0.5)) * tsize);
    vec4 c = texture(sam, (iuv + vec2(0.5, 1.5)) * tsize);
    vec4 d = texture(sam, (iuv + vec2(1.5, 1.5)) * tsize);
    return mix(mix(a, b, fuv.x), mix(c, d, fuv.x), fuv.y);

void main() {
    vec2 viewport = pushConstants.viewport;
    vec2 texelSize = pushConstants.texelSize;
    vec2 dyeTexelSize = pushConstants.dyeTexelSize;
    float dt = pushConstants.dt;
    float dissipation = pushConstants.dissipation;
    vec2 position = vec2(gl_GlobalInvocationID.x / viewport.x, gl_GlobalInvocationID.y / viewport.y);
    if (position.x > 1.0 || position.y > 1.0) {
    vec2 uv = position + texelSize/2;
    vec2 vUv = params.uv;
    vec2 coord = vUv - dt * bilerp(uVelocity, uv, texelSize).xy * texelSize;
    vec4 result = bilerp(uSource, coord, dyeTexelSize);
    float decay = 1.0 + dissipation * dt;
    imageStore(outputImage, ivec2(gl_GlobalInvocationID.xy), result / decay);

So my question is, how are quantities supposed to advect forward into a void? Am I missing part of the process here?


1 Answer 1


The cells with zero velocity are boundaries (walls) of the simulation domain. You will need to handle the boundary conditions explicitly. The advection itself cannot do much for you.

From the context it seems you are simulating incompressible flow. Intuitively, if the fluid flows toward a wall then the fluid is expected to be "pushed" to the sides of the wall. That is handled by pressure projection. Specifically it is around Equation 16 on the page you linked.

That Equation 16 is a discretization of the pressure term. That will lead you to a linear system. You can solve it iteratively (Stam99). The solution of that system will give you the divergence-free velocity, and it will advect your fluid to the sides of walls as expected.

Regarding to your question: if the boundary is a wall then nothing will advect into it. The incompressibility (the pressure here) will enforce the fluid to flow somewhere else.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.