Monday 10 May 2010

Barely Functional Fluid Dynamics (2)

Continuing on from last times adventures on Barely Functional Fluid Dynamics, it's time to start actually visualizing the results. All of the OpenGL is heavily based on the original C code by Jos Stam available here (.zip).

GLUT is the simplest way to get started with OpenGL. It's a window system toolkit designed specifically for writing OpenGL programs. The Haskell Platform comes complete with a set of Haskell bindings to GLUT. The documentation doesn't give a whole lot of help regarding how to get started, but thankfully there's a huge bundle of examples within the Cabal package that help you get started.

To draw the fluid dynamics results, I needed to have some representation of the various state. State encapsulates the various bits necessary to draw the fluid dynamics simulation. This includes the data (density and velocity) together with state about the mouse position and so on. This is stored in an IORef that provides a mutable reference within the IO monad. IORef seems very simliar to Clojure's references, with similar options for changing the state inside the reference via a supplied function.

OpenGL supports several primitive types, including lines, points, quads and triangles. Each of these primitives is specified using a sequence of vertices. For example, lines are specified as a series of pairs, and a point is drawn between each pair. Similarly quads are specified using groups of four points. For drawing the fluid we use the GL_QUAD to specify a series of squares. The colour of each vertex is the density at that particular point. In order to make it look pretty the color is based on not just the density but the x and y co-ordinates. This gives green in the top-left corner down to purple in the bottom right. It's not based on any clever principles, I just thought it looked nice!

Example colours that look nice!

To draw primitive types the Haskell GLUT library exposes a function called renderPrimitive which takes the PrimitiveMode (e.g. Lines / Quads / Polygon) and an IO action. The IO action represents the drawing of the list of points, so all we need to do is loop through each point and construct an IO action to choose a color and create the vertex.

Similarly, to draw the velocity field we just need to draw a series of lines.

Putting this all together with a few event handlers and timers gives a little OpenGL application (compile with ghc -O2 --make -lglut Main.hs MFluid.hs - I appear to need to specify -lglut manually on my system at least?).

I was quite suprised about how easy it was to get something together in Haskell to do this. Reading the web, you'd guess that it would require writing a monad tutorial and getting a deep understanding of category theory. In reality, you just download the Haskell platform and get going. There's no need to understand all the details about what's going on under the hood straight away (though no doubt it helps!).

The full code is available here. Any suggestions for improvements greatly appreciated!

Wednesday 5 May 2010

Barely Functional Fluid Dynamics

Functional Fluid Dynamics in Clojure is an implementation of a fluid dynamics simulator based on a paper by Jos Stam, Real Time Fluid Dynamics for Games. Demonstration code is available in OpenGL / C (on a side note, more papers should do this as it makes research easier to extend and far more approachable). I thought it would be fun to reimplement this in Haskell as this would give me a chance to refresh my knowledge of OpenGL (last seen many years ago in Principles of Computer Graphics.)

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.)