Fluid.cu

Go to the documentation of this file.

/*
 * Copyright 2009 StingNET.  All rights reserved.
 *
 * Any use of this source code in individual and commercial software must
 * include, in the user documentation and internal comments to the code,
 * credit to the author Dirk Hekhuis.
 *
 * The code in this program is based on code from 
 * Jos Stam's paper 'Real-Time Fluid Dynamics for Games'.
 */

#include <stdio.h>
#include <stdlib.h>
#include <cutil.h>
#include "glut.h"

#define IX(i,j) ((i)+(N+2)*(j)) 
#define SWAP(x0,x) {float *tmp = x0; x0 = x; x = tmp;} 
#define DIM 128 
#define SIZE DIM*DIM 
#define max(x,y) x>y?x:y 
#define BLOCK_SIZE 32 


void printArray (int N, float *x);


void init (void);


void display (void);


void idle (void);


void reshape (int x, int y);


void keyboard (unsigned char key, int x, int y);


void click (int button, int state, int x, int y);


void motion (int x, int y);


void timer (int value);


void updateFluid (int button, int mouseMove, int state);


__global__ void gpuKernelAdd_Source (int N, float *gpuX, float *gpuS, float dt);


__global__ void gpuKernelDiffuse(int N,int b, float *gpuX, float *gpuX0, float diff,float dt,int *gpuBoundary);


__global__ void gpuKernelAdvect(int N, int b, float *gpuD, float *gpuD0, float *gpuU, float *gpuV, float dt, int *gpuBoundary);


__global__ void gpuKernelProjectStep1(int N, float *gpuU, float *gpuV, float *gpuP, float *gpuDiv, int *gpuBoundary);


__global__ void gpuKernelProjectStep2(int N, float *gpuP, float *gpuDiv, int *gpuBoundary);


__global__ void gpuKernelProjectStep3(int N, float *gpuU, float *gpuV, float *gpuP, int *gpuBoundary);


void gpuCalcFluid(void);


float u[SIZE]; 
float v[SIZE]; 
float u_prev[SIZE]; 
float v_prev[SIZE]; 
float dens[SIZE]; 
float dens_prev[SIZE]; 
int boundary[SIZE]; 
int wWidth = max(512,DIM); 
int wHeight = max(512,DIM); 
float cellSize; 
float diff = 0.5f; 
float dt = .2f; 
float visc = 0.0f; 
int stepToggle = 0; 
int runToggle = 1; 
int consoleToggle = 0; 
int fluidToggle = 0; 
int densityToggle = 1; 
int velocityToggle = 0; 
int gridToggle = 0; 
int boundaryToggle = 1; 
int drawBoundaryToggle = 0; 
int id; 
int lastX = 0; 
int lastY = 0; 
int currentX = 0; 
int currentY = 0; 
int clicked = 0; 
int leftClick = 0; 
int N = DIM-2; 
int emitterLoc = DIM/2; 
int frameRate = 60; 
int updown = 0; 
int lastCellX; 
int lastCellY; 
float totalDensity = 0; 


int main(int argc, char *argv[]) {
        int i;

        // Initialize arrays to 0
        for (i=0; i<SIZE; i++)
        {
                u[i] = 0.0f;
                u_prev[i] = 0.0f;
                v[i] = 0.0f;
                v_prev[i] = 0.0f;
                dens[i] = 0.0f;
                dens_prev[i] = 0.0f;
                boundary[i] = 0;
        }
        
        // Initialize GLUT
        glutInit(&argc, argv);
        glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE);
        glutInitWindowSize(wWidth, wHeight);
        glutCreateWindow("Compute Stable Fluids");
        init();
        glutDisplayFunc(display);
        glutReshapeFunc(reshape);
        glutIdleFunc(idle);
        glutKeyboardFunc(keyboard);
        glutMouseFunc(click);
        glutMotionFunc(motion);
        glutTimerFunc(100, timer, frameRate);

        CUDA_SAFE_CALL(cudaSetDevice(2));

        // Start GLUTs main looop
        glutMainLoop();

        return 0;
}

void printArray (int N, float *x)
{
        int i, j;

        // Outer loop for the row
        for (i=N+1; i>=0; i--)
        {
                // Inner loop for the column
                for (j=0; j<N+2; j++)
                {
                        printf("%f ", x[IX(i,j)]);
                }

                // Star printing on the next line
                printf("\n");
        }
}

void init (void)
{
        // Set the size of the each cell to be drawn
        cellSize = wWidth/N;

        // Create a list for OpenGL
        id = glGenLists(1);

        // Create a display list to hold the cell to be drawn
        glNewList(id, GL_COMPILE);
        {
                glBegin(GL_QUADS);
                {
                        glVertex3f(-cellSize/2.0f, cellSize/2.0f, 0);
                        glVertex3f(-cellSize/2.0f, -cellSize/2.0f, 0);
                        glVertex3f(cellSize/2.0f, -cellSize/2.0f, 0);
                        glVertex3f(cellSize/2.0f, cellSize/2.0f, 0);
                }
                glEnd();
        }
        glEndList();

        // Set the clear color to black
        glClearColor(0x00, 0x00, 0x00, 1);
}

void display (void)
{
        int i, j;
        float c;

        // Clear the screen for the next drawing
        glClear (GL_COLOR_BUFFER_BIT);

        // The program is running or being stepped thru
        if(runToggle || stepToggle)
        {
                // Calculations
                gpuCalcFluid();

                // Calculate the total density
                totalDensity = 0;
                for (i=N; i>=1; i--)
                {
                        for (j=1; j<=N; j++)
                        {
                                totalDensity += dens[IX(i,j)];
                        }
                }

                /* Print to the console */
                if (consoleToggle)
                {
                        // Print out the densisty field
                        printf("\n\nDensity Field\n");
                        printArray(N, dens);

                        // Print out the velocity field in the x direction
                        printf("\n\nu Field\n");
                        printArray(N, u);

                        // Print out the velocity field in the y direction
                        printf("\n\nv Field\n");
                        printArray(N, v);

                        // Print out the total density in the simulation
                        printf("\n\nTotal Density: %f\n", totalDensity);
                }

                // Reset the step toggle
                stepToggle = 0;
        }

        // Loop thru all the elements of the arrays
        for (i=N; i>=1; i--)
        {
                for (j=1; j<=N; j++)
                {
                        /* Draws the boundary */
                        if(boundaryToggle)
                        {
                                // If there is a boundary make the cell blue
                                if(boundary[IX(i,j)] == 1)
                                        glColor3ub(0x00, 0x00, 0xFF);
                                // If there is no boundary make the cell black
                                else
                                        glColor3ub(0x00, 0x00, 0x00);

                                glPushMatrix();
                                {
                                        // Move the cell into the correct position on the screen
                                        glTranslatef(((j)*cellSize)+(cellSize/2), ((i)*cellSize)+(cellSize/2), 0);

                                        // Draw the cell
                                        glCallList(id);
                                }
                                glPopMatrix();
                        }

                        /* Draws the density field */
                        if(densityToggle)
                        {
                                if(dens[IX(i,j)] > 0)
                                {
                                        // Set the color of the cell based on the amount of density in the cell
                                        c = (dens[IX(i,j)]);
                                        if(c == 0)
                                                glColor3f(0,0,0);
                                        else if(c < .33)
                                                glColor3f(0,0,c);
                                        else if(c < .66)
                                                glColor3f(0,c,0);
                                        else
                                                glColor3f(c,0,0); 
                                        
                                        glPushMatrix();
                                        {
                                                // Move the cell into the correct position on the screen
                                                glTranslatef(((j)*cellSize)+(cellSize/2), ((i)*cellSize)+(cellSize/2), 0);

                                                // Draw the cell
                                                glCallList(id);
                                        }
                                        glPopMatrix();
                                }
                        }

                        /* Draws the velocity field */
                        if(velocityToggle)
                        {
                                // Set the color of the velocity lines to red
                                glColor3ub(0xFF, 0x00, 0x00);
                                
                                glPushMatrix();
                                {
                                        float x, y;

                                        // The x position of the end of the line
                                        x = u[IX(i,j)] * cellSize * 10;

                                        // The y position of the end of the line
                                        y = v[IX(i,j)] * cellSize * 10;

                                        // Move the line to the center of the cell
                                        glTranslatef(((j)*cellSize)+(cellSize/2), ((i)*cellSize)+(cellSize/2), 0);

                                        // Draw the line
                                        glBegin(GL_LINES);
                                        {
                                                glVertex3f(0, 0, 0);
                                                glVertex3f(x, y, 0);
                                        }
                                        glEnd();
                                }
                                glPopMatrix();
                        }

                        /* Draws the grid */
                        if(gridToggle)
                        {
                                // Set the color of the grid to green
                                glColor3ub(0x00, 0xFF, 0x00);

                                glPushMatrix();
                                {
                                        glPushAttrib(GL_POLYGON_BIT);
                                        {
                                                // Show the outline of the cell instead of filling it in
                                                glPolygonMode (GL_FRONT, GL_LINE);

                                                // Move the cell into the correct position on the screen
                                                glTranslatef(((j)*cellSize)+(cellSize/2), ((i)*cellSize)+(cellSize/2), 0);

                                                // Draw the cell
                                                glCallList(id);
                                        }
                                        glPopAttrib();
                                }
                                glPopMatrix();
                        }
                }
        }

        // Swap the back buffer with the top buffer to draw the new image
        glutSwapBuffers();
}

void idle (void)
{
        // Calls display()
        glutPostRedisplay();
}


void reshape (int x, int y)
{
        // Set the width to the new window width
        wWidth = x;

        // Set the height to the new window height
        wHeight = y;

        // Delete the current display list
        glDeleteLists(id,0);

        // Called to create a new display list with the new cell size
        init();

        // Set up the camera
        glViewport(0, 0, x, y);
        glMatrixMode(GL_PROJECTION);
        glLoadIdentity();
        glOrtho(0, wWidth, 0, wHeight, 0, 1);
        glMatrixMode(GL_MODELVIEW);
        glLoadIdentity();

        // Calls display()
        glutPostRedisplay();
}

void keyboard (unsigned char key,int x, int y)
{
        int i;

        switch(key)
        {
                // Esc key was pressed
                case 27:
                        exit(0);

                // c key was pressed
                case 'c':
                        consoleToggle = !consoleToggle;
                        break;

                // s key was pressed
                case 's':
                        runToggle = !runToggle;
                        break;

                // f key was pressed
                case 'f':
                        dens_prev[IX(1,emitterLoc)] = N * cellSize;
                        v_prev[IX(1,emitterLoc)] = cellSize;
                        break;

                // d key was pressed
                case 'd':
                        densityToggle = !densityToggle;
                        break;

                // v key was pressed
                case 'v':
                        velocityToggle = !velocityToggle;
                        break;

                // g key was pressed
                case 'g':
                        gridToggle = !gridToggle;
                        break;

                // b key was pressed
                case 'b':
                        boundaryToggle = !boundaryToggle;
                        break;

                // shift b key was pressed
                case 'B':
                        drawBoundaryToggle = !drawBoundaryToggle;
                        break;

                // r key was pressed
                case 'r':
                        // Initialize arrays to 0
                        for (i=0; i<SIZE; i++)
                        {
                                u[i] = 0.0f;
                                u_prev[i] = 0.0f;
                                v[i] = 0.0f;
                                v_prev[i] = 0.0f;
                                dens[i] = 0.0f;
                                dens_prev[i] = 0.0f;
                                boundary[i] = 0;
                        }
                        break;

                // Spacebar was pressed
                case 32:
                        stepToggle = 1;
                        break;
        }

        // Calls display()
        glutPostRedisplay();
}

void click (int button, int state, int x, int y)
{
        updown = state;

        // The left mouse button is pushed down
        if(state == 0)
        {
                lastCellX = -1;
                lastCellY = -1;
        }

        // Set the last position of the mouse
        lastX = currentX;
        lastY = currentY;

        // Set the current position of the mouse
        currentX = x;
        currentY = y;

        // Call to interact with gui and fluid
        updateFluid(button, 0, updown);
}

void motion (int x, int y)
{
        // The left mouse button is pushed down
        if(updown == 0)
        {
                lastX = currentX;
                lastY = currentY;

                currentX = x;
                currentY = y;

                // Call to interact with gui and fluid
                updateFluid(0, 1, updown);
        }
}

void timer(int value) 
{
        // Calls display()
        glutPostRedisplay();

        // Create a new timer to replace the current timer
        glutTimerFunc(1000/frameRate, timer, value);
}

void updateFluid (int button, int mouseMove, int state)
{
        float fx, fy;
        int nx, ny;

        // Find the cell that was clicked on
        fx = currentX / (float)wWidth;
        fy = (wHeight - currentY) / (float)wHeight;
        nx = (int)(fx * DIM);
        ny = (int)(fy * DIM);

        // If the mouse goes out of fluid box set it to the edge of the box
        if(nx < 1)
                nx = 1;
        else if(nx > N)
                nx = N;
        if(ny < 1)
                ny = 1;
        else if(ny > N)
                ny = N;

        // Right mouse button was clicked
        if(button == 2)
        {
                // Add density to the cell that was clicked
                dens_prev[IX(ny,nx)] = N * cellSize;

                // Add velocity to the cells around the cell that
                // density was just added to, to move the inserted density
                u_prev[IX(ny,nx-1)] = -cellSize;
                u_prev[IX(ny,nx+1)] = cellSize;
                v_prev[IX(ny-1,nx)] = -cellSize;
                v_prev[IX(ny+1,nx)] = cellSize;
        }

        // The left mouse button is pushed down, the mouse is moving and the user
        // isn't drawing boundaries
        if(button == 0 && mouseMove && drawBoundaryToggle == 0 && updown == 0)
        {
                // Add velocity to the velocity fields of the cell that was clicked on
                // The amount of velocity to add is based on the how far the mouse moves
                u_prev[IX(ny,nx)] += (currentX-lastX) * cellSize;
                v_prev[IX(ny,nx)] += (lastY-currentY) * cellSize;
        }

        // The left mouse button is pushed down and the user is drawing boundaries
        if(button == 0 && drawBoundaryToggle && updown == 0)
        {
                // There is not a boundary at the cell that is selected
                if(boundary[IX(ny,nx)] == 0)
                {
                        // The selected cell hasn't been modified in the last step
                        if(nx != lastCellX || ny != lastCellY)
                        {
                                // Make the selected cell a boundary
                                boundary[IX(ny,nx)] = 1;

                                // Set the selected cell as the last cell to be modified
                                lastCellX = nx;
                                lastCellY = ny;
                        }
                }
                
                // There is a boundary at the cell that is selected
                else
                {
                        // The selected cell hasn't been modified in the last step
                        if(lastCellX != nx || lastCellY != ny)
                        {
                                // Make the selected cell not a boundary
                                boundary[IX(ny,nx)] = 0;

                                // Set the selected cell as the last cell to be modified
                                lastCellX = nx;
                                lastCellY = ny;
                        }
                }
        }

        // Calls display()
        glutPostRedisplay();
}

__global__ void gpuKernelAdd_Source (int N, float *gpuX, float *gpuS, float dt)
{
        int i, size = (N+2)*(N+2);
        int thisEntry;

        thisEntry = blockIdx.x*BLOCK_SIZE + threadIdx.x;

        i = thisEntry % size;

        // Add the value from each position of the array gpuS multiplied by
        // the time step to the corresponding position in array gpuX
        gpuX[i] += dt*gpuS[i];
}

__global__ void gpuKernelDiffuse(int N,int b, float *gpuX, float *gpuX0, float diff,float dt,int *gpuBoundary)
{
        int i,j,k;
        float a = dt*diff*N*N;
        int thisEntry; 

        // Perform this step multiple times for increased accuracy
        for(k = 0;k < 20;k++) 
        {
                thisEntry = blockIdx.x*BLOCK_SIZE + threadIdx.x;

                j = thisEntry / (N+2);
                i = thisEntry % (N+2);

                // If the element in array x is a boundary then no fluid
                // should be at that spot so set the density at that spot to 0
                // Also don't calculate the values for the invisible boundaries,
                // just set them to 0
                if (i == 0 || i == N+1 || j == 0 || j == N+1 || gpuBoundary[thisEntry]) 
                {
                        gpuX[thisEntry] = 0.0f;
                }

                // Set the density at the current position to the density of
                // the temporary density field at the current position plus
                // the density of the cells adjacent to the current cell.
                else
                {
                        gpuX[thisEntry] = 
                        (gpuX0[thisEntry] + a*(gpuX[IX(i-1,j)] +
                                                                        gpuX[IX(i+1,j)] +
                                                                        gpuX[IX(i,j-1)] +
                                                                        gpuX[IX(i,j+1)])) / (1+4*a);
                }

                // Set the values in the invisible boundary
                // for the velocity field
                if(b == 0)
                {
                        if (i == 0) 
                        {
                                gpuX[thisEntry] = gpuX[IX(1,j)];
                        }
                        if (i == (N+1))
                        { 
                                gpuX[thisEntry] = gpuX[IX(N,i)];
                        }
                        if (j == 0) 
                        {
                                gpuX[thisEntry] = gpuX[IX(i,1)];
                        }
                        if (j == (N+1))
                        { 
                                gpuX[thisEntry] = gpuX[IX(i,N)];
                        }
        
                        if ((i == 0) && (j == 0))
                        {
                                gpuX[thisEntry] = 0.5f*(gpuX[IX(1,0  )] + gpuX[IX(0  ,1)]);
                        }
                        if ((i == 0) && (j == (N+1)))
                        {
                                gpuX[thisEntry] = 0.5f*(gpuX[IX(1,N+1)] + gpuX[IX(0  ,N)]);
                        }
                        if ((i == (N+1)) && (j == 0))
                        {
                                gpuX[thisEntry] = 0.5f*(gpuX[IX(N,0  )] + gpuX[IX(N+1,1)]);
                        }
                        if ((i == (N+1)) && (j == (N+1)))
                        {
                                gpuX[thisEntry] = 0.5f*(gpuX[IX(N,N+1)] + gpuX[IX(N+1,N)]);
                        }
                }

                // Set the values in the invisible boundary
                // for the density field
                else
                {
                        if (i == 0) 
                        {
                                gpuX[thisEntry] = 0;
                }
                        if (i == (N+1))
                { 
                                gpuX[thisEntry] = 0;
                        }
                        if (j == 0) 
                        {
                                gpuX[thisEntry] = 0;
                        }               
                        if (j == (N+1))
                        { 
                                gpuX[thisEntry] = 0;
                        }
                }
        }
}

__global__ void gpuKernelAdvect(int N, int b, float *gpuD, float *gpuD0, float *gpuU, float *gpuV, float dt, int *gpuBoundary)
{
        int i, j, i0, j0, i1, j1;
        float x, y, s0, t0, s1, t1, dt0, temp;
        int thisEntry;

        // The time step multiplied by the size of the array
        dt0 = dt*N;

        thisEntry = blockIdx.x*BLOCK_SIZE + threadIdx.x;
        j = thisEntry / (N+2);
        i = thisEntry % (N+2);

        x = i-dt0*gpuU[IX(i,j)];
        y = j-dt0*gpuV[IX(i,j)];

        if (x<0.5f)
                x = 0.5f;
        if (x>N+0.5f)
                x = N+0.5f;
        i0 = (int)x;
        i1 = i0+1;

        if (y<0.5f)
                y = 0.5f;
        if (y>N+0.5f)
                y = N+0.5f;
        j0 = (int)y;
        j1 = j0+1;

        s1 = x-i0;
        s0 = 1-s1;

        t1 = y-j0;
        t0 = 1-t1;

        // If the element in array x is a boundary then no fluid
        // should be at that spot so set the density at that spot to 0
        // Also don't calculate the values for the invisible boundaries,
        // just set them to 0
        if (i == 0 || i == N+1 || j == 0 || j == N+1 || gpuBoundary[thisEntry])
                gpuD[thisEntry] = 0.0f;

        // Set the current cells new velocity value
        else
        {
                // Find the new velocity value for the current cell
                temp = s0 * (t0*gpuD0[IX(i0,j0)] + t1*gpuD0[IX(i0,j1)]) +
                                        s1 * (t0*gpuD0[IX(i1,j0)] + t1*gpuD0[IX(i1,j1)]);

                // If velocity is moving in the horizontal direction and is moving left
                // and there is a boundary to the left of the current cell
                if(b == 1 && temp < 0 && gpuBoundary[IX(i,j-1)])
                        // Make the direction of the velocity go in the opposite direction
                        gpuD[IX(i,j)] = -temp;

                // If the velocity is moving in the horizontal direction and is moving right
                // and there is a boundary to the right of the current cell
                else if(b == 1 && temp > 0 && gpuBoundary[IX(i,j+1)])
                        // Make the direction of the velocity go in the opposite direction
                        gpuD[IX(i,j)] = -temp;

                // If the velocity is moving in the vertical direction and is moving down
                // and there is a boundary below the current cell
                else if(b == 2 && temp < 0 && gpuBoundary[IX(i-1,j)])
                        // Make the direction of the velocity go in the opposite direction
                        gpuD[IX(i,j)] = -temp;

                // If the velocity is moving in the verticle direction and is moving down
                // and there is a boundary above the current cell
                else if(b == 2 && temp > 0 && gpuBoundary[IX(i+1,j)])
                        // Make the direction of the velocity go in the opposite direction
                        gpuD[IX(i,j)] = -temp;

                // Set the new velocity value to the current cell
                else
                        gpuD[IX(i,j)] = temp;   
        }

        // Set the values in the invisible boundary
        // for the velocity field
        if(b == 0)
        {
                if (i == 0) 
                {
                        gpuD[thisEntry] = gpuD[IX(1,j)];
        }
                if (i == (N+1)) 
                { 
                        gpuD[thisEntry] = gpuD[IX(N,i)];
        }
                if (j == 0) 
                {
                        gpuD[thisEntry] = gpuD[IX(i,1)];
        }
                if (j == (N+1))
        { 
                        gpuD[thisEntry] = gpuD[IX(i,N)];
        }

        if ((i == 0) && (j == 0)) 
                {
                        gpuD[thisEntry] = 0.5f*(gpuD[IX(1,0  )] + gpuD[IX(0  ,1)]);
                }
        if ((i == 0) && (j == (N+1))) 
                {
                        gpuD[thisEntry] = 0.5f*(gpuD[IX(1,N+1)] + gpuD[IX(0  ,N)]);
                }
        if ((i == (N+1)) && (j == 0)) 
                {
                        gpuD[thisEntry] = 0.5f*(gpuD[IX(N,0  )] + gpuD[IX(N+1,1)]);
                }
        if ((i == (N+1)) && (j == (N+1))) 
                {
                        gpuD[thisEntry] = 0.5f*(gpuD[IX(N,N+1)] + gpuD[IX(N+1,N)]);
                }
        }

        // Set the values in the invisible boundary
        // for the density field
        else
        {
                if (i == 0) 
                {
                        gpuD[thisEntry] = 0;
        }
                if (i == (N+1))
        { 
                        gpuD[thisEntry] = 0;
        }
                if (j == 0) 
                {
                        gpuD[thisEntry] = 0;
        }
                if (j == (N+1))
        { 
                        gpuD[thisEntry] = 0;
        }
        }
}

__global__ void gpuKernelProjectStep1(int N, float *gpuU, float *gpuV, float *gpuP, float *gpuDiv, int *gpuBoundary)
{
        int i, j;
        float h, temp;
        int thisEntry;

        // 1 divided by the size of the array
        h = 1.0f/N;

        thisEntry = blockIdx.x*BLOCK_SIZE + threadIdx.x;
        i = thisEntry / (N+2);
        j = thisEntry % (N+2);

        // If the element in array x is a boundary then no fluid
        // should be at that spot so set the density at that spot to 0
        // Also don't calculate the values for the invisible boundaries,
        // just set them to 0
        if (i == 0 || i == N+1 || j == 0 || j == N+1 || gpuBoundary[thisEntry])
                gpuDiv[thisEntry] = 0.0f;
        else
        {       
                // Find the new velocity value of the current cell
                temp = -0.5f*h*(gpuU[IX(i+1,j)] - gpuU[IX(i-1,j)] + gpuV[IX(i,j+1)] - gpuV[IX(i,j-1)]);

                // If the velocity is moving in the vertical direction and is moving down
                // and there is a boundary below the current cell
                if(temp < 0 && gpuBoundary[IX(i-1,j)])
                        // Make the direction of the velocity go in the opposite direction
                        gpuDiv[thisEntry] = -temp;

                // If the velocity is moving in the verticle direction and is moving down
                // and there is a boundary above the current cell
                else if(temp > 0 && gpuBoundary[IX(i+1,j)])
                        // Make the direction of the velocity go in the opposite direction
                        gpuDiv[thisEntry] = -temp;
                        
                // Set the new velocity value to the current cell
                else
                        gpuDiv[thisEntry] = temp;
        }

        // Set the horizontal velocity of the current cell to 0
        gpuP[thisEntry] = 0.0f;

        __syncthreads();
        

        // Set the values in the invisible boundary
        if (i == 0) 
        {
                gpuDiv[thisEntry] = gpuDiv[IX(1,j)];
    }
        if (i == (N+1)) 
        { 
                gpuDiv[thisEntry] = gpuDiv[IX(N,i)];
    }
        if (j == 0) 
        {
                gpuDiv[thisEntry] = gpuDiv[IX(i,1)];
    }
        if (j == (N+1))
    { 
                gpuDiv[thisEntry] = gpuDiv[IX(i,N)];
    }

    if ((i == 0) && (j == 0)) 
        {
                gpuDiv[thisEntry] = 0.5f*(gpuDiv[IX(1,0  )] + gpuDiv[IX(0  ,1)]);
        }
    if ((i == 0) && (j == (N+1))) 
        {
                gpuDiv[thisEntry] = 0.5f*(gpuDiv[IX(1,N+1)] + gpuDiv[IX(0  ,N)]);
        }
    if ((i == (N+1)) && (j == 0)) 
        {
                gpuDiv[thisEntry] = 0.5f*(gpuDiv[IX(N,0  )] + gpuDiv[IX(N+1,1)]);
        }
    if ((i == (N+1)) && (j == (N+1))) 
        {
                gpuDiv[thisEntry] = 0.5f*(gpuDiv[IX(N,N+1)] + gpuDiv[IX(N+1,N)]);
        }


        // Set the values in the invisible boundary
        if (i == 0) 
        {
                gpuP[thisEntry] = gpuP[IX(1,j)];
    }
        if (i == (N+1)) 
        { 
                gpuP[thisEntry] = gpuP[IX(N,i)];
    }
        if (j == 0) 
        {
                gpuP[thisEntry] = gpuP[IX(i,1)];
    }
        if (j == (N+1))
    { 
                gpuP[thisEntry] = gpuP[IX(i,N)];
    }

    if ((i == 0) && (j == 0)) 
        {
                gpuP[thisEntry] = 0.5f*(gpuP[IX(1,0  )] + gpuP[IX(0  ,1)]);
        }
    if ((i == 0) && (j == (N+1))) 
        {
                gpuP[thisEntry] = 0.5f*(gpuP[IX(1,N+1)] + gpuP[IX(0  ,N)]);
        }
    if ((i == (N+1)) && (j == 0)) 
        {
                gpuP[thisEntry] = 0.5f*(gpuP[IX(N,0  )] + gpuP[IX(N+1,1)]);
        }
    if ((i == (N+1)) && (j == (N+1))) 
        {
                gpuP[thisEntry] = 0.5f*(gpuP[IX(N,N+1)] + gpuP[IX(N+1,N)]);
        }
}

__global__ void gpuKernelProjectStep2(int N, float *gpuP, float *gpuDiv, int *gpuBoundary)
{
        int i, j, k;
        float temp;
        int thisEntry;

        // Perform this step multiple times for increased accuracy
        for (k=0; k<20; k++)
        {                       
                thisEntry = blockIdx.x*BLOCK_SIZE + threadIdx.x;
                i = thisEntry / (N+2);
                j = thisEntry % (N+2);

                // If the element in array x is a boundary then no fluid
                // should be at that spot so set the density at that spot to 0
                // Also don't calculate the values for the invisible boundaries,
                // just set them to 0
                if (i == 0 || i == N+1 || j == 0 || j == N+1 || gpuBoundary[thisEntry])
                        gpuP[thisEntry] = 0.0f;
                else
                {
                        // Find the new velocity value of the current cell
                        temp = (gpuDiv[thisEntry] + gpuP[IX(i-1,j)] + gpuP[IX(i+1,j)] + gpuP[IX(i,j-1)] + gpuP[IX(i,j+1)]) / 4;

                        // If velocity is moving in the horizontal direction and is moving left
                        // and there is a boundary to the left of the current cell
                        if(temp < 0 && gpuBoundary[IX(i,j-1)])
                                // Make the direction of the velocity go in the opposite direction
                                gpuP[thisEntry] = -temp;
                        
                        // If the velocity is moving in the horizontal direction and is moving right
                        // and there is a boundary to the right of the current cell
                        else if(temp > 0 && gpuBoundary[IX(i,j+1)])
                                // Make the direction of the velocity go in the opposite direction
                                gpuP[thisEntry] = -temp;

                        // Set the new velocity value to the current cell
                        else
                                gpuP[thisEntry] = temp;
                }
                __syncthreads();
                
                // Set the values in the invisible boundary
                if (i == 0) 
                {
                        gpuP[thisEntry] = gpuP[IX(1,j)];
                }
                if (i == (N+1)) 
                { 
                        gpuP[thisEntry] = gpuP[IX(N,i)];
                }
                if (j == 0) 
                {
                        gpuP[thisEntry] = gpuP[IX(i,1)];
                }
                if (j == (N+1))
                { 
                        gpuP[thisEntry] = gpuP[IX(i,N)];
                }

                if ((i == 0) && (j == 0)) 
                {
                        gpuP[thisEntry] = 0.5f*(gpuP[IX(1,0  )] + gpuP[IX(0  ,1)]);
                }
                if ((i == 0) && (j == (N+1))) 
                {
                        gpuP[thisEntry] = 0.5f*(gpuP[IX(1,N+1)] + gpuP[IX(0  ,N)]);
                }
                if ((i == (N+1)) && (j == 0)) 
                {
                        gpuP[thisEntry] = 0.5f*(gpuP[IX(N,0  )] + gpuP[IX(N+1,1)]);
                }
                if ((i == (N+1)) && (j == (N+1))) 
                {
                        gpuP[thisEntry] = 0.5f*(gpuP[IX(N,N+1)] + gpuP[IX(N+1,N)]);
                }               
        }
}

__global__ void gpuKernelProjectStep3(int N, float *gpuU, float *gpuV, float *gpuP, int *gpuBoundary)
{
        int i, j;
        float h, temp;
        int thisEntry;

        // 1 divided by the size of the array
        h = 1.0f/N;

        thisEntry = blockIdx.x*BLOCK_SIZE + threadIdx.x;
        i = thisEntry / (N+2);
        j = thisEntry % (N+2);

        // If the element in array x is a boundary then no fluid
        // should be at that spot so set the density at that spot to 0
        // Also don't calculate the values for the invisible boundaries,
        // just set them to 0
        if (i == 0 || i == N+1 || j == 0 || j == N+1 || gpuBoundary[thisEntry])
        {
                gpuU[thisEntry] = 0.0f;
                gpuV[thisEntry] = 0.0f;
        }

        // The current cell is not a boundary
        else
        {       
                // Find the new horizontal velocity value of the current cell
                temp = 0.5f * (gpuP[IX(i+1,j)] - gpuP[IX(i-1,j)]) / h;

                // If velocity is moving in the horizontal direction and is moving left
                // and there is a boundary to the left of the current cell
                if(temp < 0 && gpuBoundary[IX(i,j-1)])
                        // Make the direction of the velocity go in the opposite direction
                        gpuU[thisEntry] -= -temp;

                // If the velocity is moving in the horizontal direction and is moving right
                // and there is a boundary to the right of the current cell
                else if(temp > 0 && gpuBoundary[IX(i,j+1)])
                        // Make the direction of the velocity go in the opposite direction
                        gpuU[thisEntry] -= -temp;

                // Set the new horizontal velocity value to the current cell
                else
                        gpuU[thisEntry] -= temp;

                // Find the new vertical velocity value of the current cell
                temp = 0.5f * (gpuP[IX(i,j+1)] - gpuP[IX(i,j-1)]) / h;

                // If the velocity is moving in the vertical direction and is moving down
                // and there is a boundary below the current cell
                if(temp < 0 && gpuBoundary[IX(i-1,j)])
                        // Make the direction of the velocity go in the opposite direction
                        gpuV[thisEntry] -= -temp;

                // If the velocity is moving in the verticle direction and is moving down
                // and there is a boundary above the current cell
                else if(temp > 0 && gpuBoundary[IX(i+1,j)])
                        // Make the direction of the velocity go in the opposite direction
                        gpuV[thisEntry] -= -temp;

                // Set the new vertical velocity value to the current cell
                else
                        gpuV[thisEntry] -= temp;
        }

        __syncthreads();

        
        // Set the values in the invisible boundary
        if (i == 0) 
        {
                gpuU[thisEntry] = gpuU[IX(1,j)];
        }
        if (i == (N+1)) 
        { 
                gpuU[thisEntry] = gpuU[IX(N,i)];
        }
        if (j == 0) 
        {
                gpuU[thisEntry] = gpuU[IX(i,1)];
        }
        if (j == (N+1))
        { 
                gpuU[thisEntry] = gpuU[IX(i,N)];
        }

        if ((i == 0) && (j == 0)) 
        {
                gpuU[thisEntry] = 0.5f*(gpuU[IX(1,0  )] + gpuU[IX(0  ,1)]);
        }
        if ((i == 0) && (j == (N+1))) 
        {
                gpuU[thisEntry] = 0.5f*(gpuU[IX(1,N+1)] + gpuU[IX(0  ,N)]);
        }
        if ((i == (N+1)) && (j == 0)) 
        {
                gpuU[thisEntry] = 0.5f*(gpuU[IX(N,0  )] + gpuU[IX(N+1,1)]);
        }
        if ((i == (N+1)) && (j == (N+1))) 
        {
                gpuU[thisEntry] = 0.5f*(gpuU[IX(N,N+1)] + gpuU[IX(N+1,N)]);
        }


        // Set the values in the invisible boundary
        if (i == 0) 
        {
                gpuV[thisEntry] = gpuV[IX(1,j)];
        }
        if (i == (N+1)) 
        { 
                gpuV[thisEntry] = gpuV[IX(N,i)];
        }
        if (j == 0) 
        {
                gpuV[thisEntry] = gpuV[IX(i,1)];
        }
        if (j == (N+1))
        { 
                gpuV[thisEntry] = gpuV[IX(i,N)];
        }

        if ((i == 0) && (j == 0)) 
        {
                gpuV[thisEntry] = 0.5f*(gpuV[IX(1,0  )] + gpuV[IX(0  ,1)]);
        }
        if ((i == 0) && (j == (N+1))) 
        {
                gpuV[thisEntry] = 0.5f*(gpuV[IX(1,N+1)] + gpuV[IX(0  ,N)]);
        }
        if ((i == (N+1)) && (j == 0)) 
        {
                gpuV[thisEntry] = 0.5f*(gpuV[IX(N,0  )] + gpuV[IX(N+1,1)]);
        }
        if ((i == (N+1)) && (j == (N+1))) 
        {
                gpuV[thisEntry] = 0.5f*(gpuV[IX(N,N+1)] + gpuV[IX(N+1,N)]);
        }
}

void gpuCalcFluid()
{
        float *gpuU;
        float *gpuV;
        float *gpuUPrev;
        float *gpuVPrev;
        float *gpuDens;
        float *gpuDensPrev;
        int *gpuBoundary;
        int result;

        // Allocate the arrays in the GPU memory
        result = cudaMalloc((void**)&gpuU, sizeof(float)*SIZE);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMalloc for gpuU failed.\n");
                exit(1);
        }
        result = cudaMalloc((void**)&gpuV, sizeof(float)*SIZE);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMalloc for gpuV failed.\n");
                exit(1);
        }
        result = cudaMalloc((void**)&gpuUPrev, sizeof(float)*SIZE);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMalloc for gpuUPrev failed.\n");
                exit(1);
        }
        result = cudaMalloc((void**)&gpuVPrev, sizeof(float)*SIZE);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMalloc for gpuVPrev failed.\n");
                exit(1);
        }
        result = cudaMalloc((void**)&gpuDens, sizeof(float)*SIZE);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMalloc for gpuDens failed.\n");
                exit(1);
        }
        result = cudaMalloc((void**)&gpuDensPrev, sizeof(float)*SIZE);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMalloc for gpuDensPrev failed.\n");
                exit(1);
        }
        result = cudaMalloc((void**)&gpuBoundary, sizeof(int)*SIZE);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMalloc for gpuBoundary failed.\n");
                exit(1);
        }

        // Copy the arrays from main memory to the memory of the GPU
        result = cudaMemcpy(gpuU, u, sizeof(float)*SIZE, cudaMemcpyHostToDevice);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMemcpy for gpuU failed.\n");
                exit(1);
        }
        result = cudaMemcpy(gpuV, v, sizeof(float)*SIZE, cudaMemcpyHostToDevice);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMemcpy for gpuV failed.\n");
                exit(1);
        }
        result = cudaMemcpy(gpuUPrev, u_prev, sizeof(float)*SIZE, cudaMemcpyHostToDevice);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMemcpy for gpuUPrev failed.\n");
                exit(1);
        }
        result = cudaMemcpy(gpuVPrev, v_prev, sizeof(float)*SIZE, cudaMemcpyHostToDevice);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMemcpy for gpuVPrev failed.\n");
                exit(1);
        }
        result = cudaMemcpy(gpuDens, dens, sizeof(float)*SIZE, cudaMemcpyHostToDevice);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMemcpy for gpuDens failed.\n");
                exit(1);
        }
        result = cudaMemcpy(gpuDensPrev, dens_prev, sizeof(float)*SIZE, cudaMemcpyHostToDevice);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMemcpy for gpuDensPrev failed.\n");
                exit(1);
        }
        result = cudaMemcpy(gpuBoundary, boundary, sizeof(int)*SIZE, cudaMemcpyHostToDevice);
        if (result != cudaSuccess) {
                printf("In gpuProject: cudaMemcpy for gpuBoundary failed.\n");
                exit(1);
        }

        // Indicate the dimension of the block
        dim3 dimblock(BLOCK_SIZE);
        // Indicate the dimension of the grid measured in blocks
        dim3 dimgrid(SIZE/BLOCK_SIZE);

        // Densisty Step
        gpuKernelAdd_Source<<<dimgrid,dimblock>>>(N,gpuDens,gpuDensPrev,dt);
        SWAP (gpuDensPrev, gpuDens);
        gpuKernelDiffuse<<<dimgrid,dimblock>>>(N,1,gpuDens,gpuDensPrev,diff,dt,gpuBoundary);
        SWAP (gpuDensPrev, gpuDens);
        gpuKernelAdvect<<<dimgrid,dimblock>>>(N,1,gpuDens,gpuDensPrev,gpuU,gpuV,dt,gpuBoundary);

        // Velocity Step
        gpuKernelAdd_Source<<<dimgrid,dimblock>>>(N,gpuU,gpuUPrev,dt);
        gpuKernelAdd_Source<<<dimgrid,dimblock>>>(N,gpuV,gpuVPrev,dt);
        SWAP (gpuUPrev,gpuU);
        gpuKernelDiffuse<<<dimgrid,dimblock>>>(N,0,gpuU,gpuUPrev,visc,dt,gpuBoundary);
        SWAP (gpuVPrev,gpuV);
        gpuKernelDiffuse<<<dimgrid,dimblock>>>(N,0,gpuV,gpuVPrev,visc,dt,gpuBoundary);
        gpuKernelProjectStep1<<<dimgrid,dimblock>>>(N, gpuU, gpuV, gpuUPrev, gpuVPrev, gpuBoundary);
        gpuKernelProjectStep2<<<dimgrid,dimblock>>>(N, gpuUPrev, gpuVPrev, gpuBoundary);
        gpuKernelProjectStep3<<<dimgrid,dimblock>>>(N, gpuU, gpuV, gpuUPrev, gpuBoundary);
        SWAP (gpuUPrev,gpuU);
        SWAP (gpuVPrev,gpuV);
        gpuKernelAdvect<<<dimgrid,dimblock>>>(N,0,gpuU,gpuUPrev,gpuUPrev,gpuVPrev,dt,gpuBoundary);
        gpuKernelAdvect<<<dimgrid,dimblock>>>(N,0,gpuV,gpuVPrev,gpuUPrev,gpuVPrev,dt,gpuBoundary);
        gpuKernelProjectStep1<<<dimgrid,dimblock>>>(N, gpuU, gpuV, gpuUPrev, gpuVPrev, gpuBoundary);
        gpuKernelProjectStep2<<<dimgrid,dimblock>>>(N, gpuUPrev, gpuVPrev, gpuBoundary);
        gpuKernelProjectStep3<<<dimgrid,dimblock>>>(N, gpuU, gpuV, gpuUPrev, gpuBoundary);
        
        // Copy the results back from the GPU memory to the CPU memory
        result = cudaMemcpy(u, gpuU, sizeof(float)*SIZE, cudaMemcpyDeviceToHost);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMemcpy back to host for u failed.\n");
                exit(1);
        }
        result = cudaMemcpy(v, gpuV, sizeof(float)*SIZE, cudaMemcpyDeviceToHost);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMemcpy back to host for v failed.\n");
                exit(1);
        }
        result = cudaMemcpy(u_prev, gpuUPrev, sizeof(float)*SIZE, cudaMemcpyDeviceToHost);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMemcpy back to host for u_prev failed.\n");
                exit(1);
        }
        result = cudaMemcpy(v_prev, gpuVPrev, sizeof(float)*SIZE, cudaMemcpyDeviceToHost);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMemcpy back to host for v_prev failed.\n");
                exit(1);
        }
                result = cudaMemcpy(dens, gpuDens, sizeof(float)*SIZE, cudaMemcpyDeviceToHost);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMemcpy back to host for dens failed.\n");
                exit(1);
        }
        result = cudaMemcpy(dens_prev, gpuDensPrev, sizeof(float)*SIZE, cudaMemcpyDeviceToHost);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaMemcpy back to host for dens_prev failed.\n");
                exit(1);
        }

        // Free the arrays
        result = cudaFree(gpuU);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaFree gpuU failed.\n");
                exit(1);
        }
        result = cudaFree(gpuV);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaFree gpuV failed.\n");
                exit(1);
        }
        result = cudaFree(gpuUPrev);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaFree gpuUPrev failed.\n");
                exit(1);
        }
        result = cudaFree(gpuVPrev);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaFree gpuVPrev failed.\n");
                exit(1);
        }
        result = cudaFree(gpuDens);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaFree gpuDens failed.\n");
                exit(1);
        }
        result = cudaFree(gpuDensPrev);
        if (result != cudaSuccess) {
                printf("In gpuCalcFluid: cudaFree gpuDensPrev failed.\n");
                exit(1);
        }
        result = cudaFree(gpuBoundary);
        if (result != cudaSuccess) {
                printf("In gpuProject: cudaFree gpuBoundary failed.\n");
                exit(1);
        }
}

---