26 January 2014 by Woody Chow on fluid simulation | cublas | CUDA | navier-stokes | blas | ray marching | cusparse | berkeley

Eulerian Fluid Simulation in CUDA


This was the final project for both the graduate level Computer Graphics class and parallel computing class at Berkeley. A high performance Eulerian fluid simulation kernel is implemented in CUDA. The core of this kernel is a highly parallelized Conjugate Gradient matrix solver, which heavily utilized CUDA's BLAS library. The smoke simulation data is then rendered in OpenGL using ray marching.


A 128x128x128 smoke simulation in a closed box

Implementation Details

Simulation Overview

The simulation method follows [2] closely. The incompressible flow Navier-Stokes equations can be divided into 4 main parts. These 4 parts are all written in CUDA to improve performance. Please refer to my final report for the parallel computing class, Eulerian Fluid Simulation for Computer Graphics in CUDA for details on parallelizing and optimizing the simulation.

Figure 1: Flow of the Simulation Kernel


Figure 2: Advection Formula
where q is the quantity to be advected, and x is the position.

Semi-Lagrangian Advection is used here, which is, updating quantities of each grid cell by tracking a particle centered at the grid cell back a time step. In our case of smoke, velocity, smoke concentration and temperature are advected.

Body Force

Figure 3: Buoyancy Formula
where a and b are arbitrary non negative constants.

Body force on the fluid includes gravitational force in most cases, and buoyancy in our case of smoke. Body force is applied to each grid cell individually at each time step.


One way to fake diffusion is to simply run a few particles per grid cell, starting from random positions inside the grid cell, through the velocity field.

Another approach, which is the one implemented in this final assignment, is to model diffuse using a Laplacian term in PDE as described in [2]. In short, finite differences of smoke and concentration are taken with the neighbouring cells are added to the cell, adjusted by an arbitary parameter for each respectively.

The video on the top of the page contains simulations with diffuse on and diffuse off.

Pressure Solve

Solve for pressure for each grid cell such that the result velocity field satisfies the incompressibility condition of the incompressible flow Navier-Stokes equations inside the fluid.

This involves solving a matrix since the pressure of each cell is implicitly defined and depends on the neighbouring pressure values, which happens to be a discretized 3D Poisson equation in the case of Eulerian simulation. The matrix is giantic since the dimension is number of grids by the number of grids, which is width * height * depth by width * height * depth. Luckily, the matrix happens to be symmetric positive-definite, which one can take advantage of. Conjugate Gradient method is used to solve the matrix iteratively.


The smoke is rendered by using concentration value stored in each grid cell. The concentration array is passed into OpenGL as a 3D texture.

Then, shader ray marching is performed to visualize the smoke in 3D. Ray marching march rays into the scene from the camera, and take samples of some values along the ray to calculate color, which it will be smoke concentration in our case. Ray marching consists of four steps:

  1. Position the 3D smoke texture(essentially a box) in the scene.
  2. Render the front face of the box to determine entry point of the ray.
  3. Render the back face to determine the leaving point of the ray.
  4. Take samples between the entry point and the leaving point of the ray.

Figure 4, 5: World position coordinates of the front face and the back face of the box

The color calculation loosely follows the ray marching tutorial in [1], which the transmittance of a ray is calculated by multiplying all the transmittance of the samples along the ray. Final color of each ray is calculated by summing up the contributions of each sample.

Ray marching can absolutely be performed in real time. The performance of the simulation is only hindered by the simulation engine. If performance is not a concern, path tracer that is based on radiative transfer can generate more realistic results.


  1. Magnus W., Nafees B. Z.: Volumetric Methods in Visual Effects. Siggraph 2010 Course Notes.
  2. Robert B.: Fluid Simulation for Computer Graphics. AK Peters (2008).

First published at http://woodychow.com/berkeley/cs283/fluid/ in May 2013.