Jos Stam Advection Scheme

We can now finally start with writing the first-step of our fluid simulator, the advection scheme. In order to implement this on the GPU as shaders, we need to review some information about GPU architecture so we can come up with an appropriate way to store and operate on our data.

GPU Data Structures and Operations

Whether you choose to implement the fluid simulator using compute shaders or fragment shaders is entirely up to you.

Eulerian Data on the GPU

Since we are doing an Eulerian simulation (which just means that it is a grid-based simulation), it makes sense to store our data in textures. For a fragment program, the position of the fragment would represent the position in the grid you are working on, and for a compute program, the invocation ((A.K.A. thread or GPU work unit) The smallest unit of GPU execution within a kernel function. For GPGPU programming, invocations are organized into 3D grids and are dispatched in batches called “workgroups”.
) index would tell you what position in the grid you are working on.

Alternatively, for a compute program you could also store the data in a 1D buffer and calculate the buffer index from the grid position. The benefit of doing this is the freedom to store more data. In a texture you only get up to 4 channels to store data to for each cell, however, a compute buffer can store contains structs of larger sizes making it easier to store data.

You can get around the problem of texture storage space being limited to 4-channels by simply adding more textures to your shader of course. In this book we assume the use of textures as they make our explanations easier to understand.

Slab Operations

Every step that we are going to carry out in our fluid simulation (advection, diffusion, external forces, projection, etc.) is going to be carried out as a slab operation (often abbreviated as slabop).

A slab operation involves reading all of the old simulation data in parallel, computing the new values, writing the new values to a separate temporary texture, and finally copying the data from this new texture into the original texture (which is called a texture update). This can be summarized as:

1. Read old values
2. Calculate new values
3. Write new values to temporary storage
4. After all new values are written to temporary storage, copy this data to original buffer

Classical Advection Scheme

As we discussed earlier, advection is the process by which some fluid property is moved by the velocity of the fluid. Even the velocities themselves are moved each simulation step by their own velocity.

Describing the classical advection scheme using this information because very simple. Suppose we want to advect the velocities of the fluid (that means we want to move the position of the velocities, by the velocity itself):

WARNING

The following graphic might appear incorrectly on Safari / Webkit browsers.

Issues with Classical Advection Scheme

This approach is unsuitable for our needs for a few reasons:

1. This scheme is highly unsuitable for GPU implementation because it means that multiple GPU invocations ((A.K.A. thread or GPU work unit) The smallest unit of GPU execution within a kernel function. For GPGPU programming, invocations are organized into 3D grids and are dispatched in batches called “workgroups”.
   ) can write to the same pixel. We want each pixel to be updated by one GPU invocation ((A.K.A. thread or GPU work unit) The smallest unit of GPU execution within a kernel function. For GPGPU programming, invocations are organized into 3D grids and are dispatched in batches called “workgroups”.
   ) at most.
2. The new lookup position $\vec{p}_{\text{new}}$ will very likely not lie exactly on a pixel center, so the location you are trying to write to is actually an interpolation of nearby pixels, which means you will have to update multiple pixels at once.

Stam’s Advection Scheme

A staple of writing your first real-time fluid simulator is reading the amazing paper by Stam (2003). Stam describes a different approach to advection that is much more suited to be implemented on the GPU:

1. Read the velocity $\vec{v}$ at the current position $\vec{p}$
2. Look “backwards” using this velocity: $\vec{p}_{\text{new}} = \vec{p} - \vec{v} \cdot \Delta t$
3. Look up the velocity at $\vec{p}_{\text{new}}$
4. Copy that velocity and use it as the new velocity for that point

It is very likely that $\vec{p}_{\text{new}}$ will not be located exactly on a pixel center, so instead, look up the 4 nearest pixels and perform a bilinear interpolation* to calculate the new velocity value.

If that’s difficult to visualize, I have made an animation to hopefully illustrate this advection scheme better. You can also see how part of the slab operation is performed, after this the “new velocity texture” data is copied to the “old velocity texture” so the calculation can be re-run with the new data on the next time step:

*You only need to do bilinear interpolation if you are looking up the values of the velocity texture using texel positions. If you are using a sampler2D and looking up values by UV coordinates instead, then all you need to do is ensure that your sampler2D is configured with bilinear filtering enabled and it will be done automatically for you!

// If you are using any of these methods to look up values
// then you need to do bilinear filtering yourself:
texelFetch(VELOCITY_TEXTURE, texel, 0);
imageLoad(VELOCITY_IMAGE, texel);

// If you are using a sampler, simply set up the sampler with bilinear filtering:
texture(VELOCITY_TEXTURE, uv);
texture(VELOCITY_TEXTURE, texel / (resolution - 1.0));

This type of advection scheme is called a semi-Lagrangian advection scheme. Lagrangian simulations are fluid simulations where the fluid is composed of simulated particles. On the other hand, what we are doing is called an Eulerian simulation, where our simulation data exists on a fixed grid. So the term semi-Lagrangian means that we pretend / imagine that there are particles in our fluid and evaluate their behavior, but we never actually store these particles as data anywhere. We are simply using the “concept” of particles.

Stam’s advection scheme is a better approach for two reasons:

1. It is very well suited for implementation on the GPU
2. It is unconditionally stable, which means that you can use it even when the time-step is large.

Issues with Stam’s Advection Scheme

Consider running Stam’s advection scheme on the following case. Here is an example current velocity texture:

WARNING

The following graphic might appear incorrectly on Safari / Webkit browsers.

Only the center grid position has a velocity pointing to the right, the rest of the grid positions don’t have any velocity (represented in the diagram as “0”).

For all the positions with zero velocity, their lookup position is identical to their current position, which means their output velocity will also still be zero.

However, for the center grid position, it inverts it’s velocity, looks up the values to the left (which are all zero), and copies that velocity, resulting in the following output:

WARNING

The following graphic might appear incorrectly on Safari / Webkit browsers.

This is obviously incorrect, as advection tells us that the velocity must move by the velocity itself. So the correct result should be something like this:

WARNING

The following graphic might appear incorrectly on Safari / Webkit browsers.

So what’s going on here? Why does Stam’s advection scheme work for simple fluid simulations if it gives us an obviously incorrect answer here? Take a moment to think about it.

The devil lies in the example problem itself. The example I have given you here is not a valid simulation state in the first place . Recall our incompressibility condition we discussed with the Navier-Stokes equations:

It states that the divergence of the velocity field must always be zero. Divergence is a measure of how much the vectors of a vector field are “flowing out” of a certain point. Look back at this chapter for a review. A non-zero divergence means that matter is either being created or destroyed, which violates the law of conservation of mass.

If you look at the example above, it is not possible for the center grid cell to have a non-zero velocity value when all of it’s neighbors are zero. That means that there is fluid flowing out of that cell without any fluid entering that cell, which means matter is being mysteriously created there. A valid example might look something like this instead:

WARNING

The following graphic might appear incorrectly on Safari / Webkit browsers.

Notice how the velocities flow into each other so that the velocity field has no divergence.

When we get to the later parts of our simulator where we employ the Helmholtz-Hodge decomposition to obtain a divergence-free velocity field, Stam’s advection scheme works much better and proves to be a simple yet efficient advection model.

Compute Shader Implementation

We’re going to set up quite a few things here because they are going to be used in future compute stages as well. In order to minimize and avoid unnecessary CPU-GPU communication in between our compute stages, it is better for the CPU to upload all of the required data from the first compute stage itself, and have that data common between all the compute stages.

Even if our first stage (advection) doesn’t require the pressure textures, we define it here anyways. Otherwise, the CPU has to update the storage buffers and texture buffers between every compute stage and reupload them to the GPU which wastes a lot of time as memory bandwidth is limited.

Texture Buffers Setup

As discussed earlier, our slab operation requires temporary storage for us to write our output values to. So, I have created two image textures for the velocity field on the compute device, however, at some point we will need to sample the textures themselves (for bilinear sampling of the input velocity, or for sampling the output texture to copy its data back to the input texture for the next simulation step).

So for both textures we attach the image itself and a sampler, giving us the following texture buffers:

#version 450

layout(set = 0, binding = 0, rg32f) uniform image2D VELOCITY_INPUT;
layout(set = 0, binding = 1) uniform sampler2D VELOCITY_INPUT_SAMPLER;
layout(set = 0, binding = 2, rg32f) uniform image2D VELOCITY_OUTPUT;
layout(set = 0, binding = 3) uniform sampler2D VELOCITY_OUTPUT_SAMPLER;

Notice that the image textures are of type rg32f which means they only have a red and green channel, each storing a 32-bit signed floating point number. We only need two channels because we are storing 2D velocities.

We do the same for the pressure field which is a scalar field, so the image textures here are of type r32f because we only need one channel to store scalar values:

layout(set = 0, binding = 4, r32f) uniform image2D PRESSURE_INPUT;
layout(set = 0, binding = 5) uniform sampler2D PRESSURE_INPUT_SAMPLER;
layout(set = 0, binding = 6, r32f) uniform image2D PRESSURE_OUTPUT;
layout(set = 0, binding = 7) uniform sampler2D PRESSURE_OUTPUT_SAMPLER;

Finally, later we are going to need a way to visualize our fluid. So we create appropriate textures for us to update later:

layout(set = 0, binding = 8, rgba32f) uniform image2D UPDATE_TEXTURE;
layout(set = 0, binding = 9, rgba32f) uniform image2D DISPLAY_TEXTURE;

Storage Buffers Setup

Here we want to define all the user controlled parameters and any extra data we want to send from the CPU such as the mouse position and whether the mouse button is clicked.

layout(set = 1, binding = 0, std430) restrict buffer Parameters
{
	float dt;
	float dx;
	float viscosity;
	float click_radius;
	float click_force;
}
parameters;

Finally, since the mouse data needs to be updated every frame, we send it using push constants which are extremely fast. However, the amount of data you can upload with the push constant is limited.

layout(push_constant, std430) uniform MouseData {
	vec2 mouse_pos;
};

Advection Compute

First we define a helper function to help us look up the interpolated value from the input velocity texture:

vec4 Tv(vec2 _texel)
{
	// Bilinear Texture Lookup for Velocity using sampler's filtering
	// Make sure RDSamplerState mag_filter, min_filter, mip_filter are all set to linear!!!
	vec2 resolution = vec2(textureSize(VELOCITY_INPUT_SAMPLER,0));
	vec2 uv = _texel / (resolution - 1.0);
	return texture(VELOCITY_INPUT_SAMPLER, uv);
}

Now, simply apply Stam’s advection scheme:

layout(local_size_x = 8, local_size_y = 8, local_size_z = 1) in;
void main()
{
	ivec2 texel = ivec2(gl_GlobalInvocationID.xy);
	vec2 ftexel = vec2(texel); // Convert from int to float

	// Follow velocity field "back in time"
	float dt = parameters.dt;
	vec2 lookup_pos = ftexel - Tv(ftexel).xy * dt;

	// Copy value from lookup location
	vec2 new_vel = Tv(lookup_pos).xy;

	// Write out data
	vec4 color = vec4(new_vel, 0.0, 0.0);
	imageStore(VELOCITY_OUTPUT, texel, color);
}

Finally, and arguably the most important step. Is to copy the data from VELOCITY_OUTPUT back to VELOCITY_INPUT after this pass has run to complete the slab operation. In order to do this we need will develop a second “texture copy” pass that we will re-use every time we need to carry out a slabop:

#version 450

layout(set = 0, binding = 0, rg32f) uniform image2D VELOCITY_INPUT;
layout(set = 0, binding = 1) uniform sampler2D VELOCITY_INPUT_SAMPLER;
layout(set = 0, binding = 2, rg32f) uniform image2D VELOCITY_OUTPUT;
layout(set = 0, binding = 3) uniform sampler2D VELOCITY_OUTPUT_SAMPLER;

layout(set = 0, binding = 4, r32f) uniform image2D PRESSURE_INPUT;
layout(set = 0, binding = 5) uniform sampler2D PRESSURE_INPUT_SAMPLER;
layout(set = 0, binding = 6, r32f) uniform image2D PRESSURE_OUTPUT;
layout(set = 0, binding = 7) uniform sampler2D PRESSURE_OUTPUT_SAMPLER;

layout(set = 0, binding = 8, rgba32f) uniform image2D UPDATE_TEXTURE;
layout(set = 0, binding = 9, rgba32f) uniform image2D DISPLAY_TEXTURE;

layout(set = 1, binding = 0, std430) restrict buffer Parameters
{
	float dt;
	float dx;
	float viscosity;
	float click_radius;
	float click_force;
}
parameters;

layout(push_constant, std430) uniform MouseData {
	vec2 mouse_pos;
};

layout(local_size_x = 8, local_size_y = 8, local_size_z = 1) in;
void main()
{
	// Copy data from VELOCITY_OUTPUT into VELOCITY_INPUT
	ivec2 texel = ivec2(gl_GlobalInvocationID.xy);
	imageStore(VELOCITY_INPUT, texel, imageLoad(VELOCITY_OUTPUT, texel));
}

Live Demo

I have placed initial velocities in the center that radiate outwards from the center. You can see how this advection scheme is unconditionally stable, however, because the velocities outside this circle are all zero, nothing moves.

Don’t worry about this. As we discussed earlier, once we apply the other steps of our fluid simulator, this advection scheme performs much better.

Here the absolute values of the velocity vector field are visualized. The red channel holds the absolute values of the horizontal velocity and the green channel holds the absolute values of the vertical velocity.

A Note on Alternatives

Keep in mind that other advection schemes also exist. Many are a lot more accurate than Stam’s advection scheme, however, Stam’s is one of the simplest to understand and easiest to implement.

After we have completed building a basic fluid simulator we will explore a different, more accurate advection scheme that you can choose to replace Stam’s advection scheme with.