Over the last couple of decades, the lattice Boltzmann method (LBM) has grown into a potent alternative methodology for computational fluid dynamics (CFD). As such, CFD has been traditionally developed using the Navier-Stokes equations as a starting point. But there are several advantages to using the LBM, chiefly the extreme simplicity of the approach (easy coding) and the inherently parallel nature of the computations (scales easily to large machines). The LBM is akin to an explicit numerical scheme, in that the variable at the new time level can be calculated based on the (local) values of the variable at the previous time level.

Unlike regular CFD, where the primary dependent variables are the velocity components, the density and the pressure, the LBM works with particle velocity distribution functions, or PDFs. A simplified way of understanding the LBM is to imagine that there are a set of discrete points in space (lattice points) and there is a set of finite directions at each lattice point along which fluid particles can move or "stream". During each time step, particles "stream" or hop from one lattice point to another along these finite directions and then "collide" with each other, conserving mass and momentum. This collision event leads to a new distribution of particles, which then streams again to neighboring lattice points. For more information about the LBM, you can refer to several excellent review articles available on the internet.

A simple application of the LBM is to solve for two-dimensional flow inside a square cavity. In this problem, the square cavity is completely filled with a constant viscosity fluid. The three side walls of the cavity are fixed and the top wall suddenly starts moving to the right with a constant velocity, dragging the adjacent fluid along with it and causing a re-circulating flow to develop inside the cavity. In this blog, we will provide a canonical LBM implementation of this problem (for low Reynolds numbers) and provide detailed CPU and GPU versions of the code so interested users can see for themselves how easy it is to migrate from the C++ (CPU version)  to the CUDA (GPU version) and the amazing speed-up than can result from doing so.

<strong>The CPU approach</strong>

A simple LBM code that can be implemented on the CPU is copied below. This 2D implementation uses the D2Q9 model, which means we have a 2 dimensional model with 9 discrete velocity directions.
<pre class="lang:default decode:true">/*

Simulation of flow inside a 2D square cavity
using the lattice Boltzmann method (LBM)

Written by:       Abhijit Joshi (abhijit@accelereyes.com)

Last modified on: Saturday, March 23 @ 6:41 pm

Build instructions: g++ cpu_lbm.cpp -o lbmCPU.x

Run instructions: ./lbmCPU.x

*/

#include&lt;iostream&gt;

// problem parameters
const int N = 512;            // number of node points along X and Y (cavity length in lattice units)
const int TIME_STEPS = 1000;  // number of time steps for which the simulation is run

const int NDIR = 9;           // number of discrete velocity directions used in the D2Q9 model

const double DENSITY = 2.7;          // fluid density in lattice units
const double LID_VELOCITY = 0.05;    // lid velocity in lattice units
const double REYNOLDS_NUMBER = 100;  // Re = 

int main(int argc, char *argv[])
{
    // the base vectors and associated weight coefficients
    double *ex = new double[NDIR];
    double *ey = new double[NDIR];
    double *alpha = new double[NDIR];

    ex[0] =  0.0;   ey[0] =  0.0;   alpha[0] = 4.0 /  9.0;
    ex[1] =  1.0;   ey[1] =  0.0;   alpha[1] = 1.0 /  9.0;
    ex[2] =  0.0;   ey[2] =  1.0;   alpha[2] = 1.0 /  9.0;
    ex[3] = -1.0;   ey[3] =  0.0;   alpha[3] = 1.0 /  9.0;
    ex[4] =  0.0;   ey[4] = -1.0;   alpha[4] = 1.0 /  9.0;
    ex[5] =  1.0;   ey[5] =  1.0;   alpha[5] = 1.0 / 36.0;
    ex[6] = -1.0;   ey[6] =  1.0;   alpha[6] = 1.0 / 36.0;
    ex[7] = -1.0;   ey[7] = -1.0;   alpha[7] = 1.0 / 36.0;
    ex[8] =  1.0;   ey[8] = -1.0;   alpha[8] = 1.0 / 36.0;

    // anti-vector direction for each base vector direction
    // (useful for implementing the bounce-back scheme)
    int *ant = new int[NDIR];

    ant[0] = 0;      //      6        2        5
    ant[1] = 3;      //               ^       
    ant[2] = 4;      //               |
    ant[3] = 1;      //               |  
    ant[4] = 2;      //      3 &lt;----- 0 -----&gt; 1
    ant[5] = 7;      //               |
    ant[6] = 8;      //               |
    ant[7] = 5;      //               v
    ant[8] = 6;      //      7        4        8

    // distribution functions
    double *f = new double[N*N*NDIR];
    double *feq = new double[N*N*NDIR];
    double *f_new = new double[N*N*NDIR];

    // density
    double *rh = new double[N*N];

    // velocity components 
    double *ux = new double[N*N];
    double *uy = new double[N*N];

    // initialize density and velocity fields inside the cavity
    for(int i=0;i&lt;N;i++) {
        for(int j=0;j&lt;N;j++) {
            rh[i+N*j] = DENSITY;
            ux[i+N*j] = 0;
            uy[i+N*j] = 0;
        }
    }
    for(int i=0;i&lt;N;i++) ux[i+N*(N-1)] = LID_VELOCITY;

    // assign initial values for distribution functions
    for(int i=0;i&lt;N;i++) {
        for(int j=0;j&lt;N;j++) {
            int ixy = i+N*j;
            for(int dir=0;dir&lt;NDIR;dir++) {
                int index = i+j*N+dir*N*N;
                double edotu = ex[dir]*ux[ixy] + ey[dir]*uy[ixy];
                double udotu = ux[ixy]*ux[ixy] + uy[ixy]*uy[ixy];
                feq[index] = rh[ixy] * alpha[dir] * (1 + 3*edotu + 4.5*edotu*edotu - 1.5*udotu);
                f[index] = feq[index];
                f_new[index] = feq[index];
            }
        }
    }

    // calculate fluid viscosity based on the Reynolds number
    double kinematicViscosity = LID_VELOCITY * (double)N / REYNOLDS_NUMBER;

    // calculate relaxation time tau
    double tau =  0.5 + 3.0 * kinematicViscosity;

    // time integration
    int time=0;
    while(time&lt;TIME_STEPS) {

        time++;

        // collision
        for(int i=1;i&lt;N-1;i++) {
            for(int j=1;j&lt;N-1;j++) {
                int ixy = i+N*j;
                for(int dir=0;dir&lt;NDIR;dir++) {
                    int index = i+j*N+dir*N*N;
                    double edotu = ex[dir]*ux[ixy] + ey[dir]*uy[ixy];
                    double udotu = ux[ixy]*ux[ixy] + uy[ixy]*uy[ixy];
                    feq[index] = rh[ixy] * alpha[dir] * (1 + 3*edotu + 4.5*edotu*edotu - 1.5*udotu);
                }
            }
        }

        // streaming from interior node points
        for(int i=1;i&lt;N-1;i++) {
            for(int j=1;j&lt;N-1;j++) {
                for(int dir=0;dir&lt;NDIR;dir++) {

                    int index = i+j*N+dir*N*N;   // (i,j,dir)
                    int index_new = (i+ex[dir]) + (j+ey[dir])*N + dir*N*N;
                    int index_ant = i + j*N + ant[dir]*N*N;

                    // post-collision distribution at (i,j) along "dir"
                    double f_plus = f[index] - (f[index] - feq[index])/tau;

                    if((i+ex[dir]==0) || (i+ex[dir]==N-1) || (j+ey[dir]==0) || (j+ey[dir]==N-1)) {
                        // bounce back
                        int ixy = i+ex[dir] + N*(j+ey[dir]);
                        double ubdote = ux[ixy]*ex[dir] + uy[ixy]*ey[dir];
                        f_new[index_ant] = f_plus - 6.0 * DENSITY * alpha[dir] * ubdote;
                    }
                    else {
                        // stream to neighbor
                        f_new[index_new] = f_plus;
                    }
                }
            }
        }

        // push f_new into f
        for(int index=0;index&lt;N*N*NDIR;index++) f[index] = f_new[index];

        // update density at interior nodes
        for(int i=1;i&lt;N-1;i++) {
            for(int j=1;j&lt;N-1;j++) {
                double rho=0;
                for(int dir=0;dir&lt;NDIR;dir++) {
                    int index = i+j*N+dir*N*N;
                    rho+=f_new[index];
                }
                rh[i+N*j] = rho;
            }
        }

        // update velocity at interior nodes
        for(int i=1;i&lt;N-1;i++) {
            for(int j=1;j&lt;N-1;j++) {
                double velx=0;
                double vely=0;
                for(int dir=0;dir&lt;NDIR;dir++) {
                    int index = i+j*N+dir*N*N;
                    velx+=f_new[index]*ex[dir];
                    vely+=f_new[index]*ey[dir];
                }
                ux[i+N*j] = velx/rh[i+N*j];
                uy[i+N*j] = vely/rh[i+N*j];
            }
        }

    } // while loop ends

/*
    // centerline velocity profile
    for(int i=1;i&lt;N-1;i++) {
        std::cout &lt;&lt; 0.5 + (double)(i-1) &lt;&lt; "  " &lt;&lt; uy[i + N*(N/2)] &lt;&lt; std::endl;
    }
*/
    return 0;</pre>
Like the Laplace equation discussed in a related blog, the CPU version sweeps through each lattice point in a sequential manner and finds the new distribution functions after streaming and collision at that lattice point. However, all these updates are completely independent of each other and thus, they can be performed in parallel.

<strong>The GPU approach</strong>

The basic idea behind this particular GPU implementation is very simple. We allocate and initialize all the relevant arrays directly on the GPU. Once this is done (we will assume that we have enough memory to accomplish this task), the actual updating of the PDFs because of streaming and collision is performed in parallel using thousands of threads that each work on independent lattice points. The complete CUDA code for implementing this idea is provided below for reference. As we said, this is just a demonstration of how to port the CPU based C++ code to the GPU and we do not wish to imply that this is the best and most efficient implementation of the LBM on the GPU. But a code is worth a thousand words and here is the code:
<pre class="lang:default decode:true">  
Simulation of flow inside a 2D square cavity
using the lattice Boltzmann method (LBM)

Written by:       Abhijit Joshi (abhijit@accelereyes.com)

Last modified on: Sunday, March 24 @ 9:00 am

Build instructions: nvcc -arch=sm_13 gpu_lbm.cu -o lbmGPU.x

Run instructions: ./lbmGPU.x

*/

#include&lt;iostream&gt;

// problem parameters
const int N = 512;            // number of node points along X and Y (cavity length in lattice units)
const int TIME_STEPS = 1000;  // number of time steps for which the simulation is run

const int NDIR = 9;           // number of discrete velocity directions used in the D2Q9 model

const double DENSITY = 2.7;          // fluid density in lattice units
const double LID_VELOCITY = 0.05;    // lid velocity in lattice units
const double REYNOLDS_NUMBER = 100;  // Re = 

__global__ void initialize(double *_ex, double *_ey, double *_alpha, int *_ant, double *_rh, double *_ux, double *_uy, double *_f, double *_feq, double *_f_new)
{
    // these could be initialized directly inside the GPU kernel and moved out of the CPU part...
    _ex[0] =  0.0;   _ey[0] =  0.0;   _alpha[0] = 4.0 /  9.0;
    _ex[1] =  1.0;   _ey[1] =  0.0;   _alpha[1] = 1.0 /  9.0;
    _ex[2] =  0.0;   _ey[2] =  1.0;   _alpha[2] = 1.0 /  9.0;
    _ex[3] = -1.0;   _ey[3] =  0.0;   _alpha[3] = 1.0 /  9.0;
    _ex[4] =  0.0;   _ey[4] = -1.0;   _alpha[4] = 1.0 /  9.0;
    _ex[5] =  1.0;   _ey[5] =  1.0;   _alpha[5] = 1.0 / 36.0;
    _ex[6] = -1.0;   _ey[6] =  1.0;   _alpha[6] = 1.0 / 36.0;
    _ex[7] = -1.0;   _ey[7] = -1.0;   _alpha[7] = 1.0 / 36.0;
    _ex[8] =  1.0;   _ey[8] = -1.0;   _alpha[8] = 1.0 / 36.0;

    // these could be initialized directly inside the GPU kernel and moved out of the CPU part...
    _ant[0] = 0;      //      6        2        5
    _ant[1] = 3;      //               ^       
    _ant[2] = 4;      //               |
    _ant[3] = 1;      //               |  
    _ant[4] = 2;      //      3 &lt;----- 0 -----&gt; 1
    _ant[5] = 7;      //               |
    _ant[6] = 8;      //               |
    _ant[7] = 5;      //               v
    _ant[8] = 6;      //      7        4        8

    // compute the "i" and "j" location and the "dir"
    // handled by this thread

    int i = blockIdx.x * blockDim.x + threadIdx.x ;
    int j = blockIdx.y * blockDim.y + threadIdx.y ;

    // initialize density and velocity fields inside the cavity
    _rh[i+N*j] = DENSITY;
    _ux[i+N*j] = 0;
    _uy[i+N*j] = 0;

    if(j==N-1) _ux[i+N*(N-1)] = LID_VELOCITY;

    // assign initial values for distribution functions
    int ixy = i+N*j;
    for(int dir=0;dir&lt;NDIR;dir++) {
        int index = i+j*N+dir*N*N;
        double edotu = _ex[dir]*_ux[ixy] + _ey[dir]*_uy[ixy];
        double udotu = _ux[ixy]*_ux[ixy] + _uy[ixy]*_uy[ixy];
        _feq[index] = _rh[ixy] * _alpha[dir] * (1 + 3*edotu + 4.5*edotu*edotu - 1.5*udotu);
        _f[index] = _feq[index];
        _f_new[index] = _feq[index];
    }
}

__global__ void timeIntegration(double *_ex, double *_ey, double *_alpha, int *_ant, double *_rh, double *_ux, double *_uy, double *_f, double *_feq, double *_f_new)
{
    // calculate fluid viscosity based on the Reynolds number
    double kinematicViscosity = LID_VELOCITY * (double)N / REYNOLDS_NUMBER;

    // calculate relaxation time tau
    double tau =  0.5 + 3.0 * kinematicViscosity;

    // compute the "i" and "j" location and the "dir"
    // handled by this thread

    int i = blockIdx.x * blockDim.x + threadIdx.x ;
    int j = blockIdx.y * blockDim.y + threadIdx.y ;

    {
        // collision
        if((i&gt;0) &amp;&amp; (i&lt;N-1) &amp;&amp; (j&gt;0) &amp;&amp; (j&lt;N-1)) {
                int ixy = i+N*j;
                for(int dir=0;dir&lt;NDIR;dir++) {
                    int index = i+j*N+dir*N*N;
                    double edotu = _ex[dir]*_ux[ixy] + _ey[dir]*_uy[ixy];
                    double udotu = _ux[ixy]*_ux[ixy] + _uy[ixy]*_uy[ixy];
                    _feq[index] = _rh[ixy] * _alpha[dir] * (1 + 3*edotu + 4.5*edotu*edotu - 1.5*udotu);
                }
        }

        // streaming from interior node points
        if((i&gt;0) &amp;&amp; (i&lt;N-1) &amp;&amp; (j&gt;0) &amp;&amp; (j&lt;N-1)) {
                for(int dir=0;dir&lt;NDIR;dir++) {

                    int index = i+j*N+dir*N*N;   // (i,j,dir)
                    int index_new = (i+_ex[dir]) + (j+_ey[dir])*N + dir*N*N;
                    int index_ant = i + j*N + _ant[dir]*N*N;

                    // post-collision distribution at (i,j) along "dir"
                    double f_plus = _f[index] - (_f[index] - _feq[index])/tau;

                    if((i+_ex[dir]==0) || (i+_ex[dir]==N-1) || (j+_ey[dir]==0) || (j+_ey[dir]==N-1)) {
                        // bounce back
                        int ixy = i+_ex[dir] + N*(j+_ey[dir]);
                        double ubdote = _ux[ixy]*_ex[dir] + _uy[ixy]*_ey[dir];
                        _f_new[index_ant] = f_plus - 6.0 * DENSITY * _alpha[dir] * ubdote;
                    }
                    else {
                        // stream to neighbor
                        _f_new[index_new] = f_plus;
                    }
                }
        }

        // push f_new into f
        if((i&gt;0) &amp;&amp; (i&lt;N-1) &amp;&amp; (j&gt;0) &amp;&amp; (j&lt;N-1)) {
                for(int dir=0;dir&lt;NDIR;dir++) {
                    int index = i+j*N+dir*N*N;   // (i,j,dir)
                    _f[index] = _f_new[index];
                }
        }

        // update density at interior nodes
        if((i&gt;0) &amp;&amp; (i&lt;N-1) &amp;&amp; (j&gt;0) &amp;&amp; (j&lt;N-1)) {
                double rho=0;
                for(int dir=0;dir&lt;NDIR;dir++) {
                    int index = i+j*N+dir*N*N;
                    rho+=_f_new[index];
                }
                _rh[i+N*j] = rho;
        }

        // update velocity at interior nodes
        if((i&gt;0) &amp;&amp; (i&lt;N-1) &amp;&amp; (j&gt;0) &amp;&amp; (j&lt;N-1)) {
                double velx=0;
                double vely=0;
                for(int dir=0;dir&lt;NDIR;dir++) {
                    int index = i+j*N+dir*N*N;
                    velx+=_f_new[index]*_ex[dir];
                    vely+=_f_new[index]*_ey[dir];
                }
                _ux[i+N*j] = velx/_rh[i+N*j];
                _uy[i+N*j] = vely/_rh[i+N*j];
        }
    }

}

int main(int argc, char *argv[])
{
    // the base vectors and associated weight coefficients (GPU)
    double *_ex, *_ey, *_alpha;  // pointers to device (GPU) memory
    cudaMalloc((void **)&amp;_ex,NDIR*sizeof(double));
    cudaMalloc((void **)&amp;_ey,NDIR*sizeof(double));
    cudaMalloc((void **)&amp;_alpha,NDIR*sizeof(double));

    // ant vector (GPU)
    int *_ant;  // gpu memory
    cudaMalloc((void **)&amp;_ant,NDIR*sizeof(int));

    // allocate memory on the GPU
    double *_f, *_feq, *_f_new;
    cudaMalloc((void **)&amp;_f,N*N*NDIR*sizeof(double));
    cudaMalloc((void **)&amp;_feq,N*N*NDIR*sizeof(double));
    cudaMalloc((void **)&amp;_f_new,N*N*NDIR*sizeof(double));

    double *_rh, *_ux, *_uy;
    cudaMalloc((void **)&amp;_rh,N*N*sizeof(double));
    cudaMalloc((void **)&amp;_ux,N*N*sizeof(double));
    cudaMalloc((void **)&amp;_uy,N*N*sizeof(double));

    // assign a 3D distribution of CUDA "threads" within each CUDA "block"    
    int threadsAlongX=16, threadsAlongY=16;
    dim3 dimBlock(threadsAlongX, threadsAlongY, 1);

    // calculate number of blocks along X and Y in a 2D CUDA "grid"
    dim3 dimGrid( ceil(float(N)/float(dimBlock.x)), ceil(float(N)/float(dimBlock.y)), 1 );

    initialize&lt;&lt;&lt;dimGrid,dimBlock&gt;&gt;&gt;(_ex, _ey, _alpha, _ant, _rh, _ux, _uy, _f, _feq, _f_new);

    // time integration
    int time=0;
    while(time&lt;TIME_STEPS) {

        time++;

        timeIntegration&lt;&lt;&lt;dimGrid,dimBlock &gt;&gt;&gt;(_ex, _ey, _alpha, _ant, _rh, _ux, _uy, _f, _feq, _f_new);

    }

/*
    double *uy = new double [N*N];
    cudaMemcpy(uy,_uy,N*N*sizeof(double),cudaMemcpyDeviceToHost);
    // centerline velocity profiles
    for(int i=1;i&lt;N-1;i++) {
        std::cout &lt;&lt; 0.5 + (double)(i-1) &lt;&lt; "  " &lt;&lt; uy[i + N*(N/2)] &lt;&lt; std::endl;
    }
*/
    return 0;
}</pre>
The initial part of the code starts off on the CPU and allocates space to store the PDF arrays on the GPU using cudaMalloc. After this, the main portion of the algorithm runs completely on the GPU. The only instance where we would be required to move data from the GPU to the CPU is when we want to write it to output files for post-processing purposes. This is not done in this example for clarity and conciseness.

<strong>RESULTS</strong>

We ran this problem using both the CPU version and the GPU version for a range of grid sizes, starting from a 128 x 128 lattice and going up to a 512 x 512 lattice. The results for total simulation time are summarized in the graph below. It can be observed that the GPU version completely outperforms the CPU version, even for the 128 x 128 grid.
<p style="text-align: center;"><a href="http://blog.accelereyes.com/blog/wp-content/uploads/2013/03/Screen-Shot-2013-03-25-at-2.11.49-PM.png"><img class="aligncenter  wp-image-3921" alt="Screen Shot 2013-03-25 at 2.11.49 PM" src="http://blog.accelereyes.com/blog/wp-content/uploads/2013/03/Screen-Shot-2013-03-25-at-2.11.49-PM-1024x767.png" width="614" height="460" /></a></p>
<p style="text-align: center;"></p>