To find the source of the instabilities, I pulled my code apart into more independent steps, that could individually be turned on or off. This did result in a speed hit, but allowed me to quickly trace the source of the problems to the advection routine. This is the part that moves the water along at its own speed.
To get rid of the problems, I tried the backward advection method described by Jos Stam in his 2003 article Real-Time Fluid Dynamics for Games. I needed to adjust the method a bit, because I am using a staggered grid: horizontal velocities are specified on vertical cell boundaries and vice versa. The backward advection method did solve the problem, but also causes too much smoothing. Presumably this happens because the method is not momentum-conserving. I went back to the naïve forward scheme I had been using.
Eventually, I traced part of the problem back to my treatment of boundary conditions. I had been setting the velocity of the fluid along the top wall, but this caused a discontinuous velocity field. I should have been setting the velocity inside the wall and the problem disappeared. The other instability problem was solved easily by decreasing the time step; with methods of this kind, the so-called CFL condition dictates, roughly, that the fluid must not move faster than one grid cell per time step. This makes a lot of sense, considering that interaction only takes place between neighbouring cells.
Now it was time to add a distinction between water and air. This comes down mostly to bookkeeping: full cells, empty cells and surface (partly filled) cells are all treated differently. Approximating the actual location of the surface from these values is another job. By the end of the week, I had it mostly working.
My test case is a so-called “broken dam problem”, where all of the water starts out in a block on the left side, and gravity is applied to eventually make it horizontal. The screenshot shows the situation after the water splashed up to the right side of the box, displaying a beautiful breaking wave. The little blue lines represent the surface orientation; although the orientation looks right their position is not yet correct. Other problems are that single partially-filled cells will remain suspended in the air, and that the water level will slowly sink until only a puddle is left. However, I’m sure that these are solveable problems.
Also, large velocities will exceed the CFL criterion and cause the simulation to blow up. In a game, lowering the time step is not an option, so I’m thinking to put an artificial cap on velocity to avoid this problem.
With full compiler optimization, this test can be run in real time on a single core of my Intel i7 CPU for a grid of 64 by 64. The code itself is still not optimized, and I expect to be able to increase speed at least fivefold by doing so. That would allow me to run the same thing in real time on grids of 128 by 128, which is approximately the minimum size needed for a decent full-screen simulation. The point has clearly been proven: what I want is possible. And I’m not even using the GPU yet!
Now that this first hurdle is out of the way, I will start working on the second: how do we turn a water simulation into an actual game?