Diffusion - sahilshahpatel/fluid-sim GitHub Wiki

In this chapter we will be implementing diffusion. Diffusion is the process that makes properties of the fluid spread to their surroundings over time.

Just like with advection, both velocity and the dye should diffuse over time. The rate of that diffusion is connected to the viscosity of the fluid we are simulating. For example, if you were to pour water on a flat surface it would quickly spread out into a thin layer. Honey, on the other hand, would take a lot longer to diffuse.

Now that we have a good idea of what diffusion is, let's take a look at the math.

Section 1: Diffusion equations

When we imagine diffusion we think of the boundary of the fluid spreading out, but in reality every part of the fluid is mixing with everything around it. This means every cell will spread some of it's fluid into each of its neighboring cells, and each of its neighbors will be spreading their fluid into it.

With that in mind, we might come up with the following equation (where superscript n means "next" time step and c means "current"):

$x^n_{i,j} = x^c_{i,j} + k*(x^c_{i-1,j} + x^c_{i+1,j} + x^c_{i,j-1} + x^c_{i,j+1} - 4x^c_{i,j})$

Here our next cell value is the current value plus what it receives from its neighbors, minus what it gives to its neighbors, where the amount transferred between cells is controlled by the variable k. k is itself equal to the diffusion rate times the time step.

Above we used the example of the fluid spreading to its neighbors, but the velocity of a cell should also spread to its neighbors! How much velocity diffuses in a liquid is a very important property, one which you have probably heard called viscosity.

Unfortunately, while this method makes intuitive sense, it has a problem! When the value of k gets too high, the output of these equations won't behave well. In other words, these equations are unstable.

You may be wondering why we can't just cap the value of k. There's a few reasons for this:

k is impacted by the time step which we can't exactly control
k is also formed from the viscosity, which shouldn't have bounds
We could use other methods to contain k while allowing the viscosity to be boundless, but this would ruin the unit relationship between the two and would cause changes in viscosity to be represented inaccurately

Section 2: Stable diffusion

What we need is a stable method of simulation. This means that our equations should work with any parameter values.

Just like with advection, it can be helpful to try it reversed. With that approach we might get:

$x^c_{i,j} = x^n_{i,j} - k*(x^n_{i-1,j} + x^n_{i+1,j} + x^n_{i,j-1} + x^n_{i,j+1} - 4x^n_{i,j})$

Note that these two versions are not equivalent, but both represent a similar concept of tracking the in-flow and out-flow between neighboring cells.

We can solve our new equation for the next cell value: $x^n_{i,j} = x^c_{i,j} + \frac{k*(x^n_{i-1,j} + x^n_{i+1,j} + x^n_{i,j-1} + x^n_{i,j+1})}{1+4k}$

Section 3: Jacobi iteration

But wait... we're still using the neighbors' next values to calculate our next value. So how does that work?

Well, since all the next-step values are co-dependent, what we've really created is a huge system of equations. The solution to that system of equations will be different each frame, so we need a quick method of solving it. Using matrix math we could find an exact solution, but that would be quite slow. Instead we will use a very simple algorithm to solve this system iteratively. That algorithm? Jacobi iteration.

We'll quickly note here that both Jos Stam's paper and Inspecto's video uses the term Gauss-Seidel iteration instead. You can read about the difference here, but with our GPU-based approach we have to use the less-accurate Jacobi iteration. (The GPU Gems chapter on this topic linked at the home page has some other suggestions for more accurate algorithms if you want to try those instead.)

So how does this algorithm work? Basically we set our initial guess for each next-step value to the previous value. We then use our equation to update our texture of guesses and repeat that process many times to increase accuracy. For diffusion, a good number of iterations is 20, but you can choose any value you like.

After implementing diffusion we get an output like this

Diffusion GIF

The above example used a diffusion strength of 10.

Section 5: Implementation hints

Since Jacobi iteration is used in both diffusion and projection, we recommend making a shader at glsl/jacobi.glsl for this chapter. At this point you know enough to implement it yourself, but we're including a few extra hints in case you need some help.

Shader uniform outline

Your shader uniforms should look like this glsl uniform sampler2D x; uniform sampler2D y; uniform float alpha; uniform float beta; uniform vec2 res;

How to implement multiple jacobi iterations

For jacobi iteration our loop needs to be in JavaScript, not GLSL. This is because we need to update the whole texture of guesses before making new guesses.

Code organization of jacobi iterations

In order to keep your JavaScript wrapper function adaptable it is best to keep your iterateration loop in update rather than doing it from within the function. This is because we can't swap this.outputTexture and our x texture from inside jacobi. Another option here is to create a texture this.jacobiTexture which is only used from within the jacobi function.

Now that our velocities are able to move and spread over time, we need to start enforcing our boundary conditions. That will be the focus of our next chapter.

Now is also a great time to check out Appendix A: Beautifying Your Render which shows you how to make your render shader produce a more visually appealing output.