Navier Stokes: iPhone vs iPad - reindernijhoff.net

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 real-time fluid dynamics running on the iPad (and iPhone).

Simulating fluid dynamics is a heavy process. 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:

Roughly 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 side 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 centre of the element. So element (1,0) in the fluid-matrix has the value of (1.5, 0.5) in the vector field u. The algorithm presented by Jos Stam gives us a nice (and simple) numerical solution of how these two matrices evolve over 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 colours 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 omission which will make the solution of the algorithm probably not physical 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 side 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 centre of each element of the matrix to a new position (see figure 1).

The problem with this approach is that the new position of the density will in general not be the centre 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 centre of each element of the new matrix, you trace backwards which point contains the density which would end up there. In general, this point will also not be the centre of an element (see figure 2), but you can calculate the correct value of the density at this point by linear 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 conserving. 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 modelled 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 subtracted from the original velocity field. Referring 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 colour 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-artefacts 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 artefacts and (maybe) provide an extra speed improvement.

Because I use the density field ρ directly as RGB values while rendering, you get an additive colour model: you draw on a black canvas and the colours 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 subtractive colour model. In that case, you draw on a white canvas and all colours will add up to black when mixed. I have tried this and the subtractive colour model is in fact as nice looking as the current additive colour model. As a future improvement, I will, therefore, add a configuration panel where you can switch between the two different colour models.

Future improvements

As noted before I have the intention 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