*This (old) post originally appeared on the website of Infi, the company I co-founded in 2003.*

One month ago we (Infi) bought an iPad in France (iPads are not available in the Netherlands yet).Having an iPad and being a programmer, I wanted to program an “app” for this device. After spending a lot of hours setting up Xcode, getting the right provisioning keys and a bit on the programming itself, I ended up with realtime fluid dynamics running on the iPad (and iPhone).

Simulating fluid dynamics is a heavy proces. I was therefore interested in how the new iPad would compare to my old (3GS) iPhone; I will give you the numbers below.

But first watch the video to see Navier Stokes in action:

*Fluid, version 1.0.*

## iPhone vs iPad

The algorithm used is fast, but numerically heavy. Combined with the OpenGL rendering pipeline (see below), I was very interested in the performance of the application on different Apple touch-devices. So far, I have only tested it on an iPhone 3GS and an Ipad (both fullscreen). Here are the numbers:

FPS | cpu time per frame (in seconds) | |

iPhone 3GS | ~17 | ~0.05 |

iPad | ~25 | ~0.03 |

In the rest of this blog post I will explain the algorithm I used for this application.

## Navier Stokes equations

To model the behavior of particles (the colors) in a fluid I use two equations. The first equation describes the evolution of the fluid itself, the second equation tells us how the particles will move in this fluid.

A fluid can be described by a vector field describing the direction and magnitute of the velocity of the fluid for each point in space. Given a fluid described by a vectorfield *u*, the Navier Stokes equations tell us how the fluid will change in time:

Rougly speaking, the Navier Stokes equations tell us how the change in velocity of the fluid (the left side of the equation) depends on the three terms on the right site of the equation. If we try to interpret these three terms you can (again roughly) say that the first term is some kind of “transportation of velocity by the fluid itself”, the second term is some kind of “diffusion” and the third term represents the external forces* f* acting on the fluid.

The movement of the (color) particles in the fluid is described by a very simular equation:

Hereby the change of the density field of particles *ρ* is given in terms of the transportation of the density by the fluid (the first term on the right), a diffusion term and a last term of external added density *S* (from external sources).

## Algorithm

To solve both equations I use the algorithm presented by Jos Stam in his article Real-Time Fluid Dynamics for Games (2003)^{[1]}. In this article the vector field *u* and the density field *ρ* are represented by two matrices: a fluid-matrix and a density-matrix. The elements of this matrices have the value of the corresponding fields at the position of the center of the element. So element (1,0) in the fluid-matrix has the value of (1.5, 0.5) in the vectorfield *u*. The algorithm presented by Jos Stam gives us a nice (and simple) numerical solution of how this two matrices evolve in time, using the two equations given above.

- I made two changes in his algorithm. The first is that my density-field is 3-dimensional, corresponding with the colors red, green and blue.
- The second (big) change I made to the original algorithm is that I ignore the diffusion terms in both equations completely, by stating nu and kappa are equal to zero. This is of course a big ommisision which will make the solution of the algorithm probably not fysical correct (or believable) anymore. However: the result is still looking good and it simplifies the computations a lot which gives the application a nice speedup (and that was what I was looking for ;-)).

Having only a “transportation-step” (move: the first term of the right site of the equations) and the last term (adding external forces or density) the (simplified) algorithm of Jos Stam can be represented by the next figure:

Each cycle in the algorithm consist of two steps, evolving the matrix describing the fields in time:

- Add forces/velocity. When painting on the canvas using a single finger, density is “added” to the density-matrix. If you use two or more fingers, a force is added to the fluid-matrix in the direction of the movement of your fingers.
- Move. In the second step the density (or the velocity for the fluid-matrix) is moved to a new position, using the local velocity of the fluid. I will explain this step next:

## Algorithm: move

In the second step of the algorithm the density (or velocity) is moved according to the velocity of the fluid at each point. A naive approach would move the density at the center of each element of the matrix to a new position (see figure 1).

The problem of this approach is that the new position of the density will in general not be the center of an element of the new matrix. Therefore you should, using this approach, somehow distribute the moved density at the end of this step over the elements surrounding the newly calculated positions. According to Jos Stam^{[1]} it is not trivial to find an algorithm for this, and proposed solutions often lead to unstable simulations.

A better way of “moving” the density (or velocity) is to work backwards: for the center of each element of the new matrix, you trace backward which point contains the density which would end up there. In general this point will also not be the center of an element (see figure 2), but you can calculate the correct value of the density at this point by lineair interpolation of the values of the surrounding elements. When using this “backward trace” method, the algorithm seems to be stable and behaves “well”.

## Algorithm: project

The method for moving the density and velocity presented above is good but not perfect. In general you will end up having just a little more (or less) total density after each step. Formally stated: the movement step is in general not perfectly mass conserving. For the movement of density this doesn’t really matter: it is hardly visible.

It is important however, that the velocity is mass convserving. The project-routine, presented in the article of Jos Stam, does that: it forces the velocity to be mass conserving. “Visually it forces the flow to have many vortices which produce realistic swirly-like flows. It is therefore an important part of the solver.”^{[1]}. The project-routine is called each step just before and after the movement step of the velocity.

I will roughly sketch why and how this routine works:

The first thing to note is that the fluid modeled is incompressible. There is also no fluid added (sources) or leaving (sinks) in the simulation. Every physicist (like me) will directly know that the field *u* therefore should be divergence free :). The project-routine is exactly trying to enforce this.

The first step in the project-routine is counter-intuitive: a scalar potential field is constructed (using Gauss-Seidel relaxation code) of which the gradient at each point is ‘as good as can be’ a approximation of the velocity field *u*. You can think of this potential field as a heightmap where sources correspond with peaks and sinks correspond with valleys. With this potential field a new vector field is constructed, whereby each point in this field is pointing in the direction of the ‘steepest descent’ of this “potential-heightmap”. This newly constructed field is a conservative vector field which has in fact a nonzero(!) divergence.

In the second step of the project-routine the newly constructed field is substracted from the original velocity field. Refering to the Helmholtz decomposition (which states that a field can be decomposed into an irrotational part and a source-free part) it is argued that the resulting velocity field is source-free (= divergence free) and thereby mass conserving, giving the nice curls we want. This is of course not a mathematical proof (and I doubt there is any available for this project routine) but it works :)

## Rendering

In the rendering step I render the fluid as a grid of triangles (see figure 3). The color of each corner of each triangle is given by the density of the corresponding element in the density-matrix. Because I use triangles to render the fluid, a lot of interpolation-artifacts can be seen (see for example the banner of this blog post). In a future improvement I will try to render the density to a texture which should remove these artifacts and (maybe) provide an extra speed improvement.

Because I use the density field *ρ* directly as RGB values while rendering, you get an additive color model: you draw on a black canvas and the colors red, green and blue will be added up to white when mixed.

By only changing the blending mode (in OpenGL) in the rendering pipeline, you can get a substractive color model. In that case you draw on a white canvas and all colors will be add up to black when mixed. I have tried this and the substractive color model is in fact as nice looking as the current additive color model. As a future improvement I will therefore add a configuration panel where you can switch between the two different color models.

## Future improvements

As noted before I have the intenion to make future improvements for this application. Some improvements which come to mind are:

- Give users the choice to use a substractive color model instead of the current additive color model.
- Rendering the fluid to a single texture. This gives less artifacts and possibly a speed increase.
- Wind tunnel: I want to add a ‘wind tunnel’ mode. In this mode you can draw shapes and watch the aerodynamics of your creations.

Meanwhile I’m trying to get this app in the AppStore, so you can try this need-to-have app for yourself :)

## References

- Jos Stam, Real-Time Fluid Dynamics for Games, 2003