The paper describes a simplified model for fluid dynamics based on the Navier-Stokes equations. The fluid is modelled as points in a 2D plane where each point has a density and a velocity. Given an initial density and velocity, the equations calculate a delta to the next density and velocity. The C implementation uses one dimensional arrays to store the points. I decided to try and model the C code as closly as possible in Haskell (and then hopefully contrast it with a purely functional solution). The main data structure is a
Grid, instances of which are used to store the density and velocity components. The grid consists of a vector and a size. Just like the C implementation, there's an area of padding around the outside to make calculating the edge cases simpler.
Most of the operations involve simply looping over the pixels and doing some calculations. The use of mutable vectors means almost every function is used only for its side-effects. This is really ugly, and I assume that a better Haskell programmer would separate the two more! The following functions remove a bit of boiler-plate from the code.
Diffusion models the spread of fluid between each point. Each point will lose density to its neighbours, but also increase by getting density from its own neighbours. For any given cell the net difference will be the sum of fluid leaving via the neighbours minus the amount coming in which leads to the following equation. An iterative technique, Gauss-Seidel relaxation, is used to solve the corresponding equations.
Advection models the changes due to the fluid's bulk motion. Haskell wise it follows a similar pattern as before (loop over each pixel, then write). Projection models the change in velocity.
Before this is plugging into a UI, we should profile it to find out how fast it is. I used the same board size as the Clojure example. Any comparison between the two is completely meaningless. There's different CPU's, different methods (see Clojure's diffusion which only uses a single iteration of the relaxation) and a completely different platform. (that should be enough of a disclaimer!). I used the awesome benchmarking framework Criterion and ended up with:
Running after compilation gives pretty good performance. The "set bounds" function (which isn't shown above) takes about 12 micro-seconds per iteration. Projection is a little slower taking 10 ms (this is actually slower than the Clojure implementation; the reason is that I'm doing 20 iterations, just like the C code - if I switch to a single iteration it takes 1.5 ms). Even with 20 iterations, it should still be possible to get a frame rate of around 100 / second. Now all that remains is plugging it into OpenGL and getting something like this:
(it looks much better when its moving honest! At the moment it looks like some kind of crazy echocardiogram)
The full code is available on my git hub repository. As always any suggestions on how I can improve this code much appreciated!
(update, OpenGL visualization code in next entry.)