## Intro

Sometimes when you’re crafting visual effects for games or other graphics applications, particle systems just don’t cut it and you have to move into the realm of physical simulation. It can be pretty daunting. I understand the math quite well, but have not had much experience using numerical methods or computational fluid dynamics. Fortunately, with a few simple assumptions about the problem, you can reduce it to something that’s both manageable and very satisfying once you actually get it working!

Let me first just say that there are some great resources out there on this topic… the paper ‘Real-Time Fluid Dynamics for Games‘ by Jos Stam was monumentally helpful, and most of the code I will give here is based on it (he gives the full implementation in the paper).

Also, there’s a great chapter in GPUGems that gives a more in-depth explanation of the math at work in this simulation, and proposes a GPU implementation of this same simulation. This video shows my rendition of this fluid solver in action.

**The Math**

The Navier-Stokes equations can be used to accurately model most fluid flows that occur in nature, and are thus widely used as the mathematical foundation for simulating fluids in computer graphics applications. These are the two equations we will be using for the simulation:

### Upper Equation

The upper equation models the time evolution of the velocity vector field that moves the fluid. The astute reader will notice that all the terms on the right are just accelerations. The first term is called the ‘advection term’ because it models the self-advection of the velocity field, or perhaps more clearly, the transportation of velocity due to the velocity field itself. The second term is called the ‘diffusion term’ (nu is a constant that represents the viscosity of the fluid) and is pretty self explanatory… basically it models how the velocity diffuses through the field. The third term is for all external forces.

### Lower Equation

The lower equation models the time evolution of the density scalar field (basically how much liquid is where). Rho is the density in this case and kappa is the diffusion constant (both scalars). The first term is still the advection term which transports the density along the velocity vector field. The second term is still the diffusion term and models how the density diffuses throughout the field. The third term is for externally added density.

Technically there is a third equation that the velocity field must adhere to in order for the fluid to be both mass conserving and incompressible. This third equation states that the divergence of the velocity field must be zero at all times. This will be further explored later. First we need to set up the environment for the simulation.

**XNA and Simulation Surface Setup**

The simulation I describe from here on out will be pretty close to the simplest one possible using the equations above. The simulation surface will just be an N by N grid of cells with a vertex at the center of each cell (you should be using VertexPositionColor as the vertex type for your VertexBuffer). Each frame we will perform the necessary calculations and use the results to color each of the N^2 vertices with a color that reflects the value of the density field at that point. In the video above I just used the density field value at each vertex as the green and blue components of the Color that gets applied to the vertex such that high density results in what looks like a bloomed out cyan color.

### Visualization

The trickiest part of the simulation surface is turning this field of vertices into polygons that the graphics card can actually rasterize. Since XNA 4.0 no longer supports point sprites, the method of creating triangles from this grid of vertices seems to be the easiest way to go. We have a multitude of options for how we go about this, but the one that I chose was to simply keep the N^2 vertices in a single VertexBuffer and use an IndexBuffer to define the triangles. Setting up the index buffer only has to be done once (indices is a list of ushorts to save bandwidth to the GPU… note that if your grid gets much larger than 100 x 100 then you would need to use a list of ints):

The caveat of this visualization method is obvious from the video above… there are a significant amount of coloring artifacts due to the linear interpolation of the colors across the triangles. Such is the price we pay for taking the easy way out and doing the simulation entirely on the CPU!

On the highest level the simulation can be broken down into two phases:

- Update – This is where all the computation and preparation of the vertex data for drawing occurs.
- Draw – Here we draw the simulation surface and do any other post processing we might want.

Fortunately XNA provides this setup via the Game class, which contains an Update method and a Draw method that are called in an alternating fashion (most of the time). There is one rather cryptic flag in the Game class that needs to be set to prevent XNA from throttling the performance of your simulation. It’s called IsFixedTimeStep, and setting it to false will make it so that Update and Draw are called as rapidly as possible (by default Update and Draw are called such that the maximum framerate is 60 fps… if either Update or Draw is taking too long while this flag is set to true then some Update and Draw frames simply get skipped altogether).

### Setting up the Vector and Scalar Fields

Implicitly a vector or scalar field has a value at every point in space, but for our purposes we only need the values at the center of the grid cells. The velocity and density fields will thus be represented by arrays of size (N + 2) * (N + 2). We add 2 to each dimension so that there is an extra row of cells around the edges of the grid for defining boundary conditions so the fluid doesn’t escape. Due to the advection method that we will be using I’ll just say that it’s necessary for us to maintain two arrays for each field… one that holds the current values and one that holds the values from the previous time step. This gives us 6 arrays in all if you break the velocity field into two separate arrays for the x and y components of the velocity. Now we’re ready to get to the good stuff.

## Diffusion

The intuitive way to approach the problem of advancing time varying physical quantities such as velocity according to an equation is a method known as Forward Euler Integration. So for instance to advance the position of a projectile for a given time step you would use p_current = p_old + v_current * dt + 0.5 * a_current * dt^2. Unfortunately this method is often unconditionally unstable for physical simulations… meaning that your simulation * will *numerically blow up. To get around this problem (which exists for both the diffusion and advection) we’ll use a backwards integration method. We’ll determine which values for each neighboring cell, when propagated backward in time, diffuse to yield the value that we started with. This involves solving a rather large linear system of equations, but luckily the matrix is sparse so this can be accomplished easily using a technique called Gauss-Seidel relaxation. This method of backwards diffusion is stable for any time step or diffusion constant values you throw at it. Here’s the code:

The more iterations you do for the outermost loop, the more accurate the approximation of the solution to the system will be… albeit at the cost of greatly taxing your CPU and possibly slowing the framerate. SetBoundary will be discussed later on.

## Advection

Just like the diffusion, we will use backward advection to move the density along the velocity field. This time we propagate the center of each grid cell back in time along the velocity field, obtain the density value at that point, and linearly interpolate that value with the density values at the centers of the 4 neighboring cells to get the final density for the original cell. This is the code:

## Density Step

This is what the final Density Step routine looks like:

## Velocity

The velocity step reuses pretty much the same Diffusion and Advect methods shown above to diffuse and advect the velocity field. Since the advection of the velocity field is non-linear, extra steps must be taken to ensure that the divergence of the velocity field is always zero at the end of each time step in order to retain mass conservation. The following method uses the Hemholtz-Hodge Decomposition Theorem to ensure that the velocity field always has zero divergence.

Once again, the number of iterations in the Gauss-Seidel relaxation will greatly affect the accuracy and speed of your simulation. Since Project is called multiple times per velocity step the number of iterations here has a much more profound effect on framerate. The final velocity step looks like this:

## Set Boundary

The boundary conditions must be set such that the net horizontal velocity along each cell of the vertical edges of the grid is zero and the net vertical velocity along each cell of the horizontal edges of the grid is zero. The SetBoundary method goes along each edge and adds an equal magnitude velocity reflected over the axis of the edge (such that the net velocity perpendicular to the edge is zero) to each cell. This effectively gives the appearance that the fluid is confined in the grid. Arbitrary boundaries can be set this way to simulate flow around objects (perhaps a blog for another time). Here’s the code:

## Conclusion and Extensions

Once again I highly recommend the paper by Jos Stam and the Nvidia article if you are genuinely interesting in implementing this or something similar (or if you were thoroughly confused by this abrupt summary of a difficult problem). There are certainly a number of improvements that can be made and I do intend to visit these in a future blog post… namely porting this method to the GPU. I’ll leave you with the few improvements I have made (adding buoyancy and vorticity confinement to more accurately model the smoke that might be emitted from the tip of a cigarette… ignore the awful banding):

If you have any questions about building your software product or need help with your software project then feel free to contact us.

BlogSlayer is the official blog of FrogSlayer, a custom software product development shop in Bryan/College Station, Texas. Our specialty is getting your product to market in 90 days or less. If you would like a free consultation for your project or big idea, email us at howdy@frogslayer.com. You can also connect with us on Twitter, Facebook, or LinkedIn.