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:
[3|1|2|1]
If we advect this grid according to the backwards lookup, at the next timestamp, we will have the following:
[0|3|3|2]
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:
And here's the shader I am using:
#define WORKGROUP_SIZE 16
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) {
return;
}
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?
