Projection - sahilshahpatel/fluid-sim GitHub Wiki

This chapter will be the most math-heavy of them all. To understand the goal of projection, we'll need to understand divergence and curl in vector fields.

Before we dive in, let's think about the expected behavior when the dye slides into the boundary. Right now it sort of leaks out like this

If that dye didn't leak out, where would it go? Well, it would certainly move to the sides, but it should also sort of bounce off the wall and curl backwards. By the end of this chapter we'll get something like this

You'll notice that we used the word curl to describe the mushroom-like effect that's happening here. That's because curl is also the name of a fundamental property of vector fields.

Our vector field has curl in a certain spot when a disk floating on that spot would naturally spin. This happens, for example, when the top part of the disk is moving right and the bottom is moving left.

It may be more helpful to see this in the fluid simulation itself

This render uses the techniques described in Appendix B to help show the vector arrows. If you're struggling with your own projection implementation, this can be a powerful debugging tool!

The sign of curl (positive or negative) indicates the direction of rotation (clockwise or counter-clockwise). In calculus terms we often abbreviate the curl of a vector field F as . The formal definition on Wikipedia is pretty good.

Divergence is the second important property of a vector field. A field has divergence in a certain spot when the fluid is moving inward our outward from that spot. Another descriptive names for these spots are sources and sinks.

The sign of divergence indicates if the spot is a source or sink. In calculus terms we often abbreviate the divergence of a vector field F as . The formal definition on Wikipedia is pretty good.

It turns out that divergence and curl are the fundamental aspects of all vector fields. Every vector field can be decomposed into a field that has no curl and another with no divergence that sum to the original field.

This process is called Helmholtz Decomposition, and will be integral to the projection step.

One of the stipulations on our fluid simulation is incompressibility. This is a property of fluids that says that you cannot make them more or less dense, you can only move them around. Imagine, for example, holding water in your hands. If you squeeze your hands together, the water will squirt out between your palms before they touch! This is because water is incompressible, so rather than squeeze into the smaller space the water must move elsewhere.

In our simulation, incompressibility translates to a no-divergence condition. Our user interactivity, advection, and diffusion forces can introduce divergence in our velocity field. This is why our last step is projection, where we will use Helmholtz decomposition to find the divergence-free component of our velocity field and discard the rest.

Unfortunately, there isn't a direct way to compute the divergence-free portion of a vector field. Instead we'll compute the curl-free portion and subtract it from the original.

To do this we'll first need to calculate the pressure at every cell. The equation for pressure fits the bill for solution by Jacobi iteration just like diffusion did:

We also need to use our boundary condition shader here. This is where the pure-Neumann condition comes into play. In a pure-Neumann fluid, the derivative of pressure must be 0 at the boundary. This means that the pressure of boundary cells must be the same as their inner neighbor.

So to calculate pressure we need to:

1. Calculate divergence and store it in a texture (requires a new one-line GLSL shader)
2. Use the Jacobi program to solve for the pressure texture, enforcing the boundary conditions after each iteration

The gradient of the pressure texture is the curl-free portion of our vector field (and so the part we will subtract out). We can create one last GLSL shader which calculates the gradient and subtracts it out. This shader will have the original velocity and pressure textures as inputs.

For the sake of being complete, Wikipedia's formal definition of the gradient may also be useful.

Ok, that was a lot. Let's quickly summarize:

1. Calculate the divergence of our velocity field using a new shader program and store it in a texture
2. Calculate the pressure of our velocity field using Jacobi Iteration and store it in a texture
3. Subtract the gradient of the pressure field from the original velocity field using a new shader program to remove the divergence

If you're struggling with how to implement this, check out our supporting sources in Appendix C. In particular, the GPU Gems chapter is useful here

Once you have done this the fluid will no longer be able to leak. It will be incompressible, and so it will be squished out to the sides when forced against the wall. In the end you should get something like this:

Reusing our jacobi shader from diffusion

As mentioned above, Jacobi iteration is a general technique for solving systems of equations. We are using it for both diffusion and projection, so it will be useful to generalize our shader. For both cases it is possible to distill our equation to

Here x and y are textures (they may be the same texture), and alpha and beta are constants.

We can rearrange our diffusion and projection equations to fit this form in order to be able to use the same shader for both stages. We leave this rearrangment as an exercise for the reader.

At this point you have implemented all of Jos Stam's original algorithm! If you haven't looked at Appendix A, you should definitely do that now to make your render more visually appealing.

There is also one more feature which we think is critical to creating a more interesting fluid flow: vorticity confinement. This procedure encourages the fluid to create more vorticies or whirlpools -- a feature of fluid flow that is diminished by some of our other techniques. That will be the focus of the next chapter.