The core of most real-time fluid simulators, like the one in EmberGen, are based on the "Stable Fluids" algorithm by Jos Stam, which to my knowledge was first presented at SIGGRAPH '99. This is a post about one part of this algorithm that's often underestimated: Projection
MG4_F32.mp4
The Stable Fluids algorithm solves a subset of the famous "Navier Stokes equations", which describe how fluids interact and move. In particular, it typically solves what's called the "incompressible Euler equations", where viscous forces are often ignored.
The two fundamental parts that make up this kind of fluid simulation algorithm are Advection and Projection. Each are a huge topic of their own. For an excellent visual introduction to the general algorithm itself check out Jamie Wong's article here.
Advection is the process of moving "stuff" around based on some underlying velocity field. This "stuff" can be smoke densities, colors, temperatures, and even the velocity itself.
Projection is what makes this velocity "well behaved" and "nice": It enforces incompressibility.
Incompressibility is the most important property at play here. By enforcing incompressibility we can make sure that concentrations (like smoke density) and momentum/velocities remains conserved. In other words: Things don't appear from or disappear into nothing.
If you disturb a pool of water it will start to swirl. The characteristic swirly motions of fluids is a direct side effect of incompressibility. This is not a dynamic (time dependent) effect, but happens immediately in the presence of any disturbance.
You can imagine pushing on some water surface at a single point (e.g. with a stick): The water you push wants to go somewhere, and this water will simultaneously try to push on the nearby water. At the same time there's a "void" where the water you pushed on resided. The easiest way for the water to fill in the void from where you pushed on it is for it to curl back on itself, generating a coiling motion, just like a magnetic field from a bar magnet (for similar reasons). This is the path of least resistance.
merged.mp4
The most common way to enforce incompressibility is a process called "Helmholtz decomposition". By applying this technique we end up with a "Poisson equation" for "pressure", and this pressure field is used to enforce incompressibility. The exact details does not matter for now.
The Poisson equation can be "solved" by different means, including various techniques from numerical linear algebra. We will focus on two techniques: the "Jacobi method" and the "Multigrid method". Again, the exact details here does not matter for now.
The important thing to keep in mind is that the Jacobi method is simple to understand and implement, but has severe drawbacks in terms of efficiency. The Multigrid algorithm on the other hand is much more complex to understand and to implement, but it is incredibly efficient. As a result almost every resource on the Stable Fluids algorithm opts into using the Jacobi method for the "pressure solve", which has resulted in the Jacobi method being the (I expect involuntary) default method of choice, even in "production code", e.g. here in Shadow of the Tomb Raider.
As a consequence of picking a poor solver the results of Advection gets even poorer, resulting in boring smooth fluid motion that do not swirl much at all, and any concentration of "stuff" gets smeared out, disappers, appears or remain static. Computing resources are also wasted.
Instead of focusing on improving the Projection first people instead resort to "hacks" as a crutch to overcome some of these issues, often without much care or understanding why they need them in the first place. E.g. things like "MacCormack/BFECC" and "Vorticity Confinement". While techniques like MacCormack or BFECC advection and Vorticity Confinement are brilliant and very useful in their own right, they are often not the right solution until we have a good Projection, and effectively sweeps the primary problem under the rug, hiding it.
In the remaining parts we will focus on some of the practical consequences of picking the simple Jacobi method over a more sophisticated choice like the Multigrid algorithm. Our reference will be a "simple" simulation made in a work-in-progress sparse version of EmberGen.
Here we are using a Multigrid solver with 4 V-cycles, but the details are not important. Keep in mind the amount of details, in both the motion and the smoke itself.
mg4.mp4
This is a screenshot of the last frame. Take note of the ~240ms spent on the Multigrid part, which is a little bit more than half of the entire simulation time (the other half being Advection). There's 134 million active voxels simulated, so it's a reasonably sized simulation.
Next we'll compare this with 1, 10 and 20 Jacobi iterations (10, 20 or 40 iterations are often cited as the threshold of choice). On the upper left we keep the Multigrid results, then we move to the right and down from there.
output_mg4_1_10_20.mp4
Here's a screenshot of the last frame of each case. Note how the total time spent on 20 iteration Jacobi version is getting close to the Multigrid reference already, but the results are still quite far off in terms of quality, and it's not even close to 134 million voxels yet.
For 40, 80 and 160 Jacobi iterations we're starting to see some reasonable motion, but clearly the smoke smears out and actually grow in size over time, so it's definitely not conserving any concentrations yet.
output_mg4_40_80_160.mp4
Screenshots again. Every single case is now much more expensive than the Multigrid method (which is already a bit over the top), but it is finally starting to fill up all the 134 million voxels now.
Finally we have the case of 320, 640 and 1280 iterations.
output_mg4_320_640_1280.mp4
And screenshots again. Only at the last couple setups, with close to or over 1000 iterations (!), it's starting to reasonably look like the reference Multigrid implementation, but at an immense cost.
Here's 20, 160 and 1280 Jacobi iterations at the same time. These are the same as the ones presented before, but juxtaposed to see a larger range of behaviors.
output_mg4_20_160_1280.mp4
In practice a single V-cycle of the Multigrid algorithm is equivalent to about a thousand Jacobi iterations, at the equivalent cost as around 10 Jacobi iterations. Even a single Multigrid V-cycle gives decent results. Here I'm comparing 1, 2, 4 and 8 V-cycles.
output_mg1_2_4_8.mp4
The individual screenshots of the final frame of the Multigrid examples.
Hopefully this serves as an illustration of the importance of having a solid Projection solver for a realtime fluid simulation, and hopefully serves as motivation to learn about more complex algorithms like Multigrid. Unfortunately the literature and resources on the Multigrid algorithm are quite sparse and often convoluted and riddled in academic noise. I hope that I will one day be able to write an approachable introduction to it.
I have glossed over a lot in this post, and there's so much to cover for completeness. The Multigrid algorithm is not perfect, and not sufficient in isolation to get great results. Different interpolation schemes for advection and grid types are really important, too. There are also a lot of good alternatives to Multigrid, like the (Preconditioned) Conjugate Gradient algorithm, which in some cases might be preferable (e.g. for liquids). Multigrid can even serve as the preconditioner here. But all this is a topic for another day.
One shouldn't take this as a signal to dismiss the Jacobi method, as it is indeed a fundamental building block of the Multigrid method, and some of its properties (the very property that makes it slow) is what makes the Multigrid method so efficient!
The videos on this page are not functioning, Just a heads up. Ive supplied screenshots of a few of them, however it seems to be all of them
