diff --git a/README.md b/README.md index d63a6a1..2c9437a 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,46 @@ **University of Pennsylvania, CIS 565: GPU Programming and Architecture, Project 1 - Flocking** -* (TODO) YOUR NAME HERE - * (TODO) [LinkedIn](), [personal website](), [twitter](), etc. -* Tested on: (TODO) Windows 22, i7-2222 @ 2.22GHz 22GB, GTX 222 222MB (Moore 2222 Lab) +* Emily Vo + * [LinkedIn](linkedin.com/in/emilyvo), [personal website](emilyhvo.com) +* Tested on: Windows 10, i7-7700HQ @ 2.8GHz 16GB, GTX 1060 6GB (Personal Computer) +Updated the CMakeLists.txt to sm_61. -### (TODO: Your README) +### README -Include screenshots, analysis, etc. (Remember, this is public, so don't put -anything here that you don't want to share with the world.) +![](images/boids.gif) + +This simulation has thousands of boids that exhibit flocking behavior. Their velocity is computed based on the positions of their neighbors, and the velocity computations follow three rules. + +The three rules are: +1) cohesion, where each boid flocks towards a perceived center of mass. +2) separation, where each boid maintains a minimum distance from others +3) where nearby boids will try to match each others velocities. + +Over time, you can see that the flock behavior becomes more orderly. + +This simulation was first implemented with naive neighbor searching, where each boid iterated through every other boid in the simulation to find the nearest neighbors. It was then optimized through a uniform grid data structure, that allowed the boids to look for neighbors in nearby grid cells. Coherence was also implemented to reduce pointer indirection when searching the uniform grid data structure. + +### PERFORMANCE ANALYSIS AND CONCEPT QUESTIONS +* For each implementation, how does changing the number of boids affect performance? Why do you think this is? +![](images/FPS_vs_numBoids.PNG) +![](images/FPS_vs_numBoids_visOFF.PNG) +Generally, increasing the number of boids decreases performance. However, when we use the uniform data structure and coherent data structure, having too few boids decreased performance because of spacial overhead (there is time spent preparing the data structures, and with so few boids, the performance benefit in the search step is not significant). + +* For each implementation, how does changing the block count and block size affect performance? Why do you think this is? +![](images/FPS_vs_threadsPerBlocks.PNG) +Increasing the block size increased performance until the threads per block exceeded 32. There is no performance increase because my GPU's warp count is 32. + +* For the coherent uniform grid: did you experience any performance improvements with the more coherent uniform grid? Was this the outcome you expected? Why or why not? + * For a smaller number of boids, the coherent uniform grid actually detracted from performance because we had spacial overhead. For a higher number of boids, the coherent uniform grid performed significantly better because there was less pointer indirection. The results were as expected, because conceptually, speed is significantly affected by cache. Furthermore, the pointer might point to something that isn't cached. The value will then have to be fetched. The program can stall twice if the target pointer's value is also not in cache. + +![](images/Uniform_vs_Coherent.PNG) + +* Did changing cell width and checking 27 vs 8 neighboring cells affect performance? Why or why not? Be careful: it is insufficient (and possibly incorrect) to say that 27-cell is slower simply because there are more cells to check! + + * *Changing the cell width*: To test how changing the cell width affected performanced, I increased a dividing factor by the power of two for each run. As the cell width decreased, the grid resolution increased. This led to a decrease in performance because it would increase the amount of memory used. The grid memory had to be reset and updated each time, so with a larger array used for the grid, there would have to be more time spent pre-processing per frame. + + ![](images/Changing_Cell_Width.PNG) + + + * *Changing the number of cell checks*: When I wanted to change my code to check 27 cells, I simply had to take away consideration of the octant and increase the range of the offset from the grid cell I was currently in. At 10,000 boids, the 27-cell check had a frame rate of 798.8 FPS. The 8-cell check had a frame rate of 975.0 FPS. At 100,000 boids, the 27-cell check had a frame rate of 180.0 FPS. The 8-Cell Check had a frame rate of 279.4 FPS. When deconstructing our algorithm, we know that the grid cell size is two times the rule's radius, so when extending that to 3 dimensions, we only need to check 8 neighboring cells. Checking 27-cells is unnecessary because we iterate over many more boids that are not even in our rule's range. \ No newline at end of file diff --git a/images/Changing_Cell_Width.PNG b/images/Changing_Cell_Width.PNG new file mode 100644 index 0000000..749bea1 Binary files /dev/null and b/images/Changing_Cell_Width.PNG differ diff --git a/images/FPS_vs_numBoids.PNG b/images/FPS_vs_numBoids.PNG new file mode 100644 index 0000000..bf38eb5 Binary files /dev/null and b/images/FPS_vs_numBoids.PNG differ diff --git a/images/FPS_vs_numBoids_visOFF.PNG b/images/FPS_vs_numBoids_visOFF.PNG new file mode 100644 index 0000000..5f4e1fb Binary files /dev/null and b/images/FPS_vs_numBoids_visOFF.PNG differ diff --git a/images/FPS_vs_threadsPerBlocks.PNG b/images/FPS_vs_threadsPerBlocks.PNG new file mode 100644 index 0000000..5728952 Binary files /dev/null and b/images/FPS_vs_threadsPerBlocks.PNG differ diff --git a/images/Uniform_vs_Coherent.PNG b/images/Uniform_vs_Coherent.PNG new file mode 100644 index 0000000..db148a1 Binary files /dev/null and b/images/Uniform_vs_Coherent.PNG differ diff --git a/images/boids.gif b/images/boids.gif new file mode 100644 index 0000000..6437e95 Binary files /dev/null and b/images/boids.gif differ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index fdd636d..b737097 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -10,5 +10,5 @@ set(SOURCE_FILES cuda_add_library(src ${SOURCE_FILES} - OPTIONS -arch=sm_20 + OPTIONS -arch=sm_61 ) diff --git a/src/kernel.cu b/src/kernel.cu index 74dffcb..cee0f4d 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -21,14 +21,14 @@ * Check for CUDA errors; print and exit if there was a problem. */ void checkCUDAError(const char *msg, int line = -1) { - cudaError_t err = cudaGetLastError(); - if (cudaSuccess != err) { - if (line >= 0) { - fprintf(stderr, "Line %d: ", line); - } - fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString(err)); - exit(EXIT_FAILURE); - } + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err) { + if (line >= 0) { + fprintf(stderr, "Line %d: ", line); + } + fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString(err)); + exit(EXIT_FAILURE); + } } @@ -70,6 +70,9 @@ glm::vec3 *dev_pos; glm::vec3 *dev_vel1; glm::vec3 *dev_vel2; +glm::vec3 *coherent_pos; +glm::vec3 *coherent_vel; + // LOOK-2.1 - these are NOT allocated for you. You'll have to set up the thrust // pointers on your own too. @@ -99,13 +102,13 @@ glm::vec3 gridMinimum; ******************/ __host__ __device__ unsigned int hash(unsigned int a) { - a = (a + 0x7ed55d16) + (a << 12); - a = (a ^ 0xc761c23c) ^ (a >> 19); - a = (a + 0x165667b1) + (a << 5); - a = (a + 0xd3a2646c) ^ (a << 9); - a = (a + 0xfd7046c5) + (a << 3); - a = (a ^ 0xb55a4f09) ^ (a >> 16); - return a; + a = (a + 0x7ed55d16) + (a << 12); + a = (a ^ 0xc761c23c) ^ (a >> 19); + a = (a + 0x165667b1) + (a << 5); + a = (a + 0xd3a2646c) ^ (a << 9); + a = (a + 0xfd7046c5) + (a << 3); + a = (a ^ 0xb55a4f09) ^ (a >> 16); + return a; } /** @@ -113,10 +116,10 @@ __host__ __device__ unsigned int hash(unsigned int a) { * Function for generating a random vec3. */ __host__ __device__ glm::vec3 generateRandomVec3(float time, int index) { - thrust::default_random_engine rng(hash((int)(index * time))); - thrust::uniform_real_distribution unitDistrib(-1, 1); + thrust::default_random_engine rng(hash((int)(index * time))); + thrust::uniform_real_distribution unitDistrib(-1, 1); - return glm::vec3((float)unitDistrib(rng), (float)unitDistrib(rng), (float)unitDistrib(rng)); + return glm::vec3((float)unitDistrib(rng), (float)unitDistrib(rng), (float)unitDistrib(rng)); } /** @@ -124,52 +127,72 @@ __host__ __device__ glm::vec3 generateRandomVec3(float time, int index) { * CUDA kernel for generating boids with a specified mass randomly around the star. */ __global__ void kernGenerateRandomPosArray(int time, int N, glm::vec3 * arr, float scale) { - int index = (blockIdx.x * blockDim.x) + threadIdx.x; - if (index < N) { - glm::vec3 rand = generateRandomVec3(time, index); - arr[index].x = scale * rand.x; - arr[index].y = scale * rand.y; - arr[index].z = scale * rand.z; - } + int index = (blockIdx.x * blockDim.x) + threadIdx.x; + if (index < N) { + glm::vec3 rand = generateRandomVec3(time, index); + arr[index].x = scale * rand.x; + arr[index].y = scale * rand.y; + arr[index].z = scale * rand.z; + } } /** * Initialize memory, update some globals */ void Boids::initSimulation(int N) { - numObjects = N; - dim3 fullBlocksPerGrid((N + blockSize - 1) / blockSize); - - // LOOK-1.2 - This is basic CUDA memory management and error checking. - // Don't forget to cudaFree in Boids::endSimulation. - cudaMalloc((void**)&dev_pos, N * sizeof(glm::vec3)); - checkCUDAErrorWithLine("cudaMalloc dev_pos failed!"); - - cudaMalloc((void**)&dev_vel1, N * sizeof(glm::vec3)); - checkCUDAErrorWithLine("cudaMalloc dev_vel1 failed!"); - - cudaMalloc((void**)&dev_vel2, N * sizeof(glm::vec3)); - checkCUDAErrorWithLine("cudaMalloc dev_vel2 failed!"); - - // LOOK-1.2 - This is a typical CUDA kernel invocation. - kernGenerateRandomPosArray<<>>(1, numObjects, - dev_pos, scene_scale); - checkCUDAErrorWithLine("kernGenerateRandomPosArray failed!"); - - // LOOK-2.1 computing grid params - gridCellWidth = 2.0f * std::max(std::max(rule1Distance, rule2Distance), rule3Distance); - int halfSideCount = (int)(scene_scale / gridCellWidth) + 1; - gridSideCount = 2 * halfSideCount; - - gridCellCount = gridSideCount * gridSideCount * gridSideCount; - gridInverseCellWidth = 1.0f / gridCellWidth; - float halfGridWidth = gridCellWidth * halfSideCount; - gridMinimum.x -= halfGridWidth; - gridMinimum.y -= halfGridWidth; - gridMinimum.z -= halfGridWidth; - - // TODO-2.1 TODO-2.3 - Allocate additional buffers here. - cudaDeviceSynchronize(); + numObjects = N; + dim3 fullBlocksPerGrid((N + blockSize - 1) / blockSize); + + // LOOK-1.2 - This is basic CUDA memory management and error checking. + // Don't forget to cudaFree in Boids::endSimulation. + cudaMalloc((void**)&dev_pos, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_pos failed!"); + + cudaMalloc((void**)&dev_vel1, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_vel1 failed!"); + + cudaMalloc((void**)&coherent_pos, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_pos failed!"); + + cudaMalloc((void**)&coherent_vel, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_vel1 failed!"); + + cudaMalloc((void**)&dev_vel2, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_vel2 failed!"); + + // LOOK-1.2 - This is a typical CUDA kernel invocation. + kernGenerateRandomPosArray << > > (1, numObjects, + dev_pos, scene_scale); + checkCUDAErrorWithLine("kernGenerateRandomPosArray failed!"); + + // LOOK-2.1 computing grid params + float cellWidthDivider = 1.0; + gridCellWidth = (1.0 / cellWidthDivider) * 2.0f * std::max(std::max(rule1Distance, rule2Distance), rule3Distance); + int halfSideCount = (int)(scene_scale / gridCellWidth) + 1; + gridSideCount = 2 * halfSideCount; + + gridCellCount = gridSideCount * gridSideCount * gridSideCount; + gridInverseCellWidth = 1.0f / gridCellWidth; + float halfGridWidth = gridCellWidth * halfSideCount; + gridMinimum.x -= halfGridWidth; + gridMinimum.y -= halfGridWidth; + gridMinimum.z -= halfGridWidth; + + // TODO-2.1 TODO-2.3 - Allocate additional buffers here. + + cudaMalloc((void**)&dev_particleArrayIndices, N * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_pos failed!"); + + cudaMalloc((void**)&dev_particleGridIndices, N * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_vel1 failed!"); + + cudaMalloc((void**)&dev_gridCellStartIndices, gridCellCount * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_pos failed!"); + + cudaMalloc((void**)&dev_gridCellEndIndices, gridCellCount * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_vel1 failed!"); + + cudaDeviceSynchronize(); } @@ -181,41 +204,41 @@ void Boids::initSimulation(int N) { * Copy the boid positions into the VBO so that they can be drawn by OpenGL. */ __global__ void kernCopyPositionsToVBO(int N, glm::vec3 *pos, float *vbo, float s_scale) { - int index = threadIdx.x + (blockIdx.x * blockDim.x); + int index = threadIdx.x + (blockIdx.x * blockDim.x); - float c_scale = -1.0f / s_scale; + float c_scale = -1.0f / s_scale; - if (index < N) { - vbo[4 * index + 0] = pos[index].x * c_scale; - vbo[4 * index + 1] = pos[index].y * c_scale; - vbo[4 * index + 2] = pos[index].z * c_scale; - vbo[4 * index + 3] = 1.0f; - } + if (index < N) { + vbo[4 * index + 0] = pos[index].x * c_scale; + vbo[4 * index + 1] = pos[index].y * c_scale; + vbo[4 * index + 2] = pos[index].z * c_scale; + vbo[4 * index + 3] = 1.0f; + } } __global__ void kernCopyVelocitiesToVBO(int N, glm::vec3 *vel, float *vbo, float s_scale) { - int index = threadIdx.x + (blockIdx.x * blockDim.x); - - if (index < N) { - vbo[4 * index + 0] = vel[index].x + 0.3f; - vbo[4 * index + 1] = vel[index].y + 0.3f; - vbo[4 * index + 2] = vel[index].z + 0.3f; - vbo[4 * index + 3] = 1.0f; - } + int index = threadIdx.x + (blockIdx.x * blockDim.x); + + if (index < N) { + vbo[4 * index + 0] = vel[index].x + 0.3f; + vbo[4 * index + 1] = vel[index].y + 0.3f; + vbo[4 * index + 2] = vel[index].z + 0.3f; + vbo[4 * index + 3] = 1.0f; + } } /** * Wrapper for call to the kernCopyboidsToVBO CUDA kernel. */ void Boids::copyBoidsToVBO(float *vbodptr_positions, float *vbodptr_velocities) { - dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); + dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); - kernCopyPositionsToVBO << > >(numObjects, dev_pos, vbodptr_positions, scene_scale); - kernCopyVelocitiesToVBO << > >(numObjects, dev_vel1, vbodptr_velocities, scene_scale); + kernCopyPositionsToVBO << > > (numObjects, dev_pos, vbodptr_positions, scene_scale); + kernCopyVelocitiesToVBO << > > (numObjects, dev_vel1, vbodptr_velocities, scene_scale); - checkCUDAErrorWithLine("copyBoidsToVBO failed!"); + checkCUDAErrorWithLine("copyBoidsToVBO failed!"); - cudaDeviceSynchronize(); + cudaDeviceSynchronize(); } @@ -230,10 +253,55 @@ void Boids::copyBoidsToVBO(float *vbodptr_positions, float *vbodptr_velocities) * in the `pos` and `vel` arrays. */ __device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *pos, const glm::vec3 *vel) { - // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves - // Rule 2: boids try to stay a distance d away from each other - // Rule 3: boids try to match the speed of surrounding boids - return glm::vec3(0.0f, 0.0f, 0.0f); + // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves + // Rule 2: boids try to stay a distance d away from each other + // Rule 3: boids try to match the speed of surrounding boids + // Update position by velocity + glm::vec3 thisPos = pos[iSelf]; + glm::vec3 thisVel = vel[iSelf]; + + int rule1N = 0; + int rule3N = 0; + + glm::vec3 center(0.0f, 0.0f, 0.0f); + glm::vec3 separate(0.0f, 0.0f, 0.0f); + glm::vec3 cohesion(0.0f, 0.0f, 0.0f); + + glm::vec3 newVel(0.0f, 0.0f, 0.0f); + + for (int i = 0; i < N; i++) { + if (i != iSelf) { + glm::vec3 thatPos = pos[i]; + glm::vec3 thatVel = vel[i]; + float distance = glm::length(thisPos - thatPos); + if (distance < rule1Distance) { + rule1N++; + center += thatPos; + } + + if (distance < rule2Distance) { + separate -= (thatPos - thisPos); + } + + if (distance < rule3Distance) { + rule3N++; + cohesion += thatVel; + } + } + } + + if (rule1N > 0) { + center /= rule1N; + newVel += (center - thisPos) * rule1Scale; + } + if (rule3N > 0) { + cohesion /= rule3N; + newVel += cohesion * rule3Scale; + } + + newVel += separate * rule2Scale; + + return newVel; } /** @@ -241,10 +309,20 @@ __device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *po * For each of the `N` bodies, update its position based on its current velocity. */ __global__ void kernUpdateVelocityBruteForce(int N, glm::vec3 *pos, - glm::vec3 *vel1, glm::vec3 *vel2) { - // Compute a new velocity based on pos and vel1 - // Clamp the speed - // Record the new velocity into vel2. Question: why NOT vel1? + glm::vec3 *vel1, glm::vec3 *vel2) { + // Compute a new velocity based on pos and vel1 + // Clamp the speed + // Record the new velocity into vel2. Question: why NOT vel1? + int index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index >= N) { + return; + } + + vel2[index] = vel1[index] + computeVelocityChange(N, index, pos, vel1); + float speed = glm::length(vel2[index]); + if (speed > maxSpeed) { + vel2[index] = glm::normalize(vel2[index]) * maxSpeed; + } } /** @@ -252,24 +330,24 @@ __global__ void kernUpdateVelocityBruteForce(int N, glm::vec3 *pos, * For each of the `N` bodies, update its position based on its current velocity. */ __global__ void kernUpdatePos(int N, float dt, glm::vec3 *pos, glm::vec3 *vel) { - // Update position by velocity - int index = threadIdx.x + (blockIdx.x * blockDim.x); - if (index >= N) { - return; - } - glm::vec3 thisPos = pos[index]; - thisPos += vel[index] * dt; - - // Wrap the boids around so we don't lose them - thisPos.x = thisPos.x < -scene_scale ? scene_scale : thisPos.x; - thisPos.y = thisPos.y < -scene_scale ? scene_scale : thisPos.y; - thisPos.z = thisPos.z < -scene_scale ? scene_scale : thisPos.z; - - thisPos.x = thisPos.x > scene_scale ? -scene_scale : thisPos.x; - thisPos.y = thisPos.y > scene_scale ? -scene_scale : thisPos.y; - thisPos.z = thisPos.z > scene_scale ? -scene_scale : thisPos.z; - - pos[index] = thisPos; + // Update position by velocity + int index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index >= N) { + return; + } + glm::vec3 thisPos = pos[index]; + thisPos += vel[index] * dt; + + // Wrap the boids around so we don't lose them + thisPos.x = thisPos.x < -scene_scale ? scene_scale : thisPos.x; + thisPos.y = thisPos.y < -scene_scale ? scene_scale : thisPos.y; + thisPos.z = thisPos.z < -scene_scale ? scene_scale : thisPos.z; + + thisPos.x = thisPos.x > scene_scale ? -scene_scale : thisPos.x; + thisPos.y = thisPos.y > scene_scale ? -scene_scale : thisPos.y; + thisPos.z = thisPos.z > scene_scale ? -scene_scale : thisPos.z; + + pos[index] = thisPos; } // LOOK-2.1 Consider this method of computing a 1D index from a 3D grid index. @@ -279,179 +357,435 @@ __global__ void kernUpdatePos(int N, float dt, glm::vec3 *pos, glm::vec3 *vel) { // for(y) // for(z)? Or some other order? __device__ int gridIndex3Dto1D(int x, int y, int z, int gridResolution) { - return x + y * gridResolution + z * gridResolution * gridResolution; + return x + y * gridResolution + z * gridResolution * gridResolution; } __global__ void kernComputeIndices(int N, int gridResolution, - glm::vec3 gridMin, float inverseCellWidth, - glm::vec3 *pos, int *indices, int *gridIndices) { - // TODO-2.1 - // - Label each boid with the index of its grid cell. - // - Set up a parallel array of integer indices as pointers to the actual - // boid data in pos and vel1/vel2 + glm::vec3 gridMin, float inverseCellWidth, + glm::vec3 *pos, int *indices, int *gridIndices) { + + int index = (blockIdx.x * blockDim.x) + threadIdx.x; + if (index >= N) { + return; + } + + glm::vec3 thisPos = pos[index]; + glm::vec3 cellPos = glm::floor((thisPos - gridMin) * inverseCellWidth); + int cellIdx = gridIndex3Dto1D(cellPos[0], cellPos[1], cellPos[2], gridResolution); + + // TODO-2.1 + // - Label each boid with the index of its grid cell. + gridIndices[index] = cellIdx; + + // - Set up a parallel array of integer indices as pointers to the actual + // boid data in pos and vel1/vel2 + indices[index] = index; + } // LOOK-2.1 Consider how this could be useful for indicating that a cell // does not enclose any boids __global__ void kernResetIntBuffer(int N, int *intBuffer, int value) { - int index = (blockIdx.x * blockDim.x) + threadIdx.x; - if (index < N) { - intBuffer[index] = value; - } + int index = (blockIdx.x * blockDim.x) + threadIdx.x; + if (index < N) { + intBuffer[index] = value; + } } __global__ void kernIdentifyCellStartEnd(int N, int *particleGridIndices, - int *gridCellStartIndices, int *gridCellEndIndices) { - // TODO-2.1 - // Identify the start point of each cell in the gridIndices array. - // This is basically a parallel unrolling of a loop that goes - // "this index doesn't match the one before it, must be a new cell!" + int *gridCellStartIndices, int *gridCellEndIndices) { + // TODO-2.1 + // Identify the start point of each cell in the gridIndices array. + // This is basically a parallel unrolling of a loop that goes + // "this index doesn't match the one before it, must be a new cell!" + int arrIdx = (blockIdx.x * blockDim.x) + threadIdx.x; + if (arrIdx >= N) { + return; + } + + int gridIdx = particleGridIndices[arrIdx]; + int nextGridIdx = particleGridIndices[arrIdx + 1]; + + if (arrIdx == N - 1) { + gridCellEndIndices[gridIdx] = arrIdx; + return; + } + if (arrIdx == 0) { + gridCellStartIndices[gridIdx] = arrIdx; + } + if (gridIdx != nextGridIdx) { + gridCellEndIndices[gridIdx] = arrIdx; + gridCellStartIndices[nextGridIdx] = arrIdx + 1; + } + +} + +__global__ void kernFlattenPosVel(int N, glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *coherent_pos, glm::vec3 *coherent_vel, int *particleArrayIndices) { + int idx = (blockIdx.x * blockDim.x) + threadIdx.x; + if (idx >= N) { + return; + } + + int particleIdx = particleArrayIndices[idx]; + coherent_pos[idx] = pos[particleIdx]; + coherent_vel[idx] = vel1[particleIdx]; } __global__ void kernUpdateVelNeighborSearchScattered( - int N, int gridResolution, glm::vec3 gridMin, - float inverseCellWidth, float cellWidth, - int *gridCellStartIndices, int *gridCellEndIndices, - int *particleArrayIndices, - glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) { - // TODO-2.1 - Update a boid's velocity using the uniform grid to reduce - // the number of boids that need to be checked. - // - Identify the grid cell that this particle is in - // - Identify which cells may contain neighbors. This isn't always 8. - // - For each cell, read the start/end indices in the boid pointer array. - // - Access each boid in the cell and compute velocity change from - // the boids rules, if this boid is within the neighborhood distance. - // - Clamp the speed change before putting the new speed in vel2 + int N, int gridResolution, glm::vec3 gridMin, + float inverseCellWidth, float cellWidth, + int *gridCellStartIndices, int *gridCellEndIndices, + int *particleArrayIndices, + glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) { + int idx = (blockIdx.x * blockDim.x) + threadIdx.x; + if (idx >= N) { + return; + } + // TODO-2.1 - Update a boid's velocity using the uniform grid to reduce + // the number of boids that need to be checked. + // - Identify the grid cell that this particle is in + glm::vec3 thisPos = pos[idx]; + glm::vec3 thisVel = vel1[idx]; + + int rule1N = 0; + int rule3N = 0; + + glm::vec3 center(0.0f, 0.0f, 0.0f); + glm::vec3 separate(0.0f, 0.0f, 0.0f); + glm::vec3 cohesion(0.0f, 0.0f, 0.0f); + + glm::vec3 newVel(0.0f, 0.0f, 0.0f); + + glm::vec3 gridPos = (thisPos - gridMin) * inverseCellWidth; + glm::vec3 gridIdx = glm::floor(gridPos); + int flatIdx = gridIndex3Dto1D(gridIdx[0], gridIdx[1], gridIdx[2], gridResolution); + + glm::vec3 octant = glm::sign((thisPos - gridMin) - (gridIdx + glm::vec3(0.5)) * cellWidth); + + // - Identify which cells may contain neighbors. This isn't always 8. + + for (int z = 0; z <= 1; z++) { + for (int y = 0; y <= 1; y++) { + for (int x = 0; x <= 1; x++) { + glm::vec3 offset = glm::vec3(x, y, z) * octant; + glm::vec3 neighborIdx = gridIdx + offset; + if (neighborIdx[0] < 0 || neighborIdx[0] >= gridResolution || + neighborIdx[1] < 0 || neighborIdx[1] >= gridResolution || + neighborIdx[2] < 0 || neighborIdx[2] >= gridResolution) { + continue; + } + + + int flatNeighborIdx = gridIndex3Dto1D(neighborIdx[0], neighborIdx[1], neighborIdx[2], gridResolution); + // - For each cell, read the start/end indices in the boid pointer array. + for (int i = gridCellStartIndices[flatNeighborIdx]; + i <= gridCellEndIndices[flatNeighborIdx]; i++) { + if (i < 0) break; + int boidIdx = particleArrayIndices[i]; + if (boidIdx != idx) { + glm::vec3 thatPos = pos[boidIdx]; + glm::vec3 thatVel = vel1[boidIdx]; + float distance = glm::length(thisPos - thatPos); + if (distance < rule1Distance) { + rule1N++; + center += thatPos; + } + + if (distance < rule2Distance) { + separate -= (thatPos - thisPos); + } + + if (distance < rule3Distance) { + rule3N++; + cohesion += thatVel; + } + } + } + + } + } + } + + if (rule1N > 0) { + center /= rule1N; + newVel += (center - thisPos) * rule1Scale; + } + if (rule3N > 0) { + cohesion /= rule3N; + newVel += cohesion * rule3Scale; + } + + newVel += thisVel + separate * rule2Scale; + float speed = glm::length(newVel); + if (speed > maxSpeed) { + newVel = glm::normalize(newVel) * maxSpeed; + } + + vel2[idx] = newVel; } __global__ void kernUpdateVelNeighborSearchCoherent( - int N, int gridResolution, glm::vec3 gridMin, - float inverseCellWidth, float cellWidth, - int *gridCellStartIndices, int *gridCellEndIndices, - glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) { - // TODO-2.3 - This should be very similar to kernUpdateVelNeighborSearchScattered, - // except with one less level of indirection. - // This should expect gridCellStartIndices and gridCellEndIndices to refer - // directly to pos and vel1. - // - Identify the grid cell that this particle is in - // - Identify which cells may contain neighbors. This isn't always 8. - // - For each cell, read the start/end indices in the boid pointer array. - // DIFFERENCE: For best results, consider what order the cells should be - // checked in to maximize the memory benefits of reordering the boids data. - // - Access each boid in the cell and compute velocity change from - // the boids rules, if this boid is within the neighborhood distance. - // - Clamp the speed change before putting the new speed in vel2 + int N, int gridResolution, glm::vec3 gridMin, + float inverseCellWidth, float cellWidth, + int *gridCellStartIndices, int *gridCellEndIndices, + glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) { + // TODO-2.3 - This should be very similar to kernUpdateVelNeighborSearchScattered, + // except with one less level of indirection. + // This should expect gridCellStartIndices and gridCellEndIndices to refer + // directly to pos and vel1. + // - Identify the grid cell that this particle is in + // - Identify which cells may contain neighbors. This isn't always 8. + // - For each cell, read the start/end indices in the boid pointer array. + // DIFFERENCE: For best results, consider what order the cells should be + // checked in to maximize the memory benefits of reordering the boids data. + // - Access each boid in the cell and compute velocity change from + // the boids rules, if this boid is within the neighborhood distance. + // - Clamp the speed change before putting the new speed in vel2 + int idx = (blockIdx.x * blockDim.x) + threadIdx.x; + if (idx >= N) { + return; + } + // TODO-2.1 - Update a boid's velocity using the uniform grid to reduce + // the number of boids that need to be checked. + // - Identify the grid cell that this particle is in + glm::vec3 thisPos = pos[idx]; + glm::vec3 thisVel = vel1[idx]; + + int rule1N = 0; + int rule3N = 0; + + glm::vec3 center(0.0f, 0.0f, 0.0f); + glm::vec3 separate(0.0f, 0.0f, 0.0f); + glm::vec3 cohesion(0.0f, 0.0f, 0.0f); + + glm::vec3 newVel(0.0f, 0.0f, 0.0f); + + glm::vec3 gridPos = (thisPos - gridMin) * inverseCellWidth; + glm::vec3 gridIdx = glm::floor(gridPos); + int flatIdx = gridIndex3Dto1D(gridIdx[0], gridIdx[1], gridIdx[2], gridResolution); + + glm::vec3 octant = glm::sign((thisPos - gridMin) - (gridIdx + glm::vec3(0.5)) * cellWidth); + + // - Identify which cells may contain neighbors. This isn't always 8. + + for (int z = 0; z <= 1; z++) { + for (int y = 0; y <= 1; y++) { + for (int x = 0; x <= 1; x++) { + glm::vec3 offset = glm::vec3(x, y, z) * octant; + glm::vec3 neighborIdx = gridIdx + offset; + if (neighborIdx[0] < 0 || neighborIdx[0] >= gridResolution || + neighborIdx[1] < 0 || neighborIdx[1] >= gridResolution || + neighborIdx[2] < 0 || neighborIdx[2] >= gridResolution) { + continue; + } + + + int flatNeighborIdx = gridIndex3Dto1D(neighborIdx[0], neighborIdx[1], neighborIdx[2], gridResolution); + // - For each cell, read the start/end indices in the boid pointer array. + for (int i = gridCellStartIndices[flatNeighborIdx]; + i <= gridCellEndIndices[flatNeighborIdx]; i++) { + if (i < 0) break; + int boidIdx = i; + if (boidIdx != idx) { + glm::vec3 thatPos = pos[boidIdx]; + glm::vec3 thatVel = vel1[boidIdx]; + float distance = glm::length(thisPos - thatPos); + if (distance < rule1Distance) { + rule1N++; + center += thatPos; + } + + if (distance < rule2Distance) { + separate -= (thatPos - thisPos); + } + + if (distance < rule3Distance) { + rule3N++; + cohesion += thatVel; + } + } + } + + } + } + } + + if (rule1N > 0) { + center /= rule1N; + newVel += (center - thisPos) * rule1Scale; + } + if (rule3N > 0) { + cohesion /= rule3N; + newVel += cohesion * rule3Scale; + } + + newVel += thisVel + separate * rule2Scale; + float speed = glm::length(newVel); + if (speed > maxSpeed) { + newVel = glm::normalize(newVel) * maxSpeed; + } + + vel2[idx] = newVel; } /** * Step the entire N-body simulation by `dt` seconds. */ void Boids::stepSimulationNaive(float dt) { - // TODO-1.2 - use the kernels you wrote to step the simulation forward in time. - // TODO-1.2 ping-pong the velocity buffers + // TODO-1.2 - use the kernels you wrote to step the simulation forward in time. + // TODO-1.2 ping-pong the velocity buffers + dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); + kernUpdatePos << > > (numObjects, dt, dev_pos, dev_vel1); + kernUpdateVelocityBruteForce << > > (numObjects, dev_pos, dev_vel1, dev_vel2); + std::swap(dev_vel1, dev_vel2); } void Boids::stepSimulationScatteredGrid(float dt) { - // TODO-2.1 - // Uniform Grid Neighbor search using Thrust sort. - // In Parallel: - // - label each particle with its array index as well as its grid index. - // Use 2x width grids. - // - Unstable key sort using Thrust. A stable sort isn't necessary, but you - // are welcome to do a performance comparison. - // - Naively unroll the loop for finding the start and end indices of each - // cell's data pointers in the array of boid indices - // - Perform velocity updates using neighbor search - // - Update positions - // - Ping-pong buffers as needed + // TODO-2.1 + // Uniform Grid Neighbor search using Thrust sort. + // In Parallel: + // - label each particle with its array index as well as its grid index. + // Use 2x width grids. + // - Unstable key sort using Thrust. A stable sort isn't necessary, but you + // are welcome to do a performance comparison. + // - Naively unroll the loop for finding the start and end indices of each + // cell's data pointers in the array of boid indices + // - Perform velocity updates using neighbor search + // - Update positions + // - Ping-pong buffers as needed + dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); + dim3 fullBlocksPerGridCells((gridCellCount + blockSize - 1) / blockSize); + + kernResetIntBuffer << > > (gridCellCount, dev_gridCellStartIndices, -1); + kernResetIntBuffer << > > (gridCellCount, dev_gridCellEndIndices, -1); + + kernComputeIndices << > > (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, + dev_pos, dev_particleArrayIndices, dev_particleGridIndices); + thrust::device_ptr dev_thrust_grid_indices(dev_particleGridIndices); + thrust::device_ptr dev_thrust_array_indices(dev_particleArrayIndices); + thrust::sort_by_key(dev_thrust_grid_indices, dev_thrust_grid_indices + numObjects, dev_thrust_array_indices); + + kernIdentifyCellStartEnd << > > (numObjects, dev_particleGridIndices, dev_gridCellStartIndices, + dev_gridCellEndIndices); + kernUpdateVelNeighborSearchScattered << > > (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, + gridCellWidth, dev_gridCellStartIndices, dev_gridCellEndIndices, + dev_particleArrayIndices, dev_pos, dev_vel1, dev_vel2); + kernUpdatePos << > > (numObjects, dt, dev_pos, dev_vel2); + std::swap(dev_vel1, dev_vel2); } void Boids::stepSimulationCoherentGrid(float dt) { - // TODO-2.3 - start by copying Boids::stepSimulationNaiveGrid - // Uniform Grid Neighbor search using Thrust sort on cell-coherent data. - // In Parallel: - // - Label each particle with its array index as well as its grid index. - // Use 2x width grids - // - Unstable key sort using Thrust. A stable sort isn't necessary, but you - // are welcome to do a performance comparison. - // - Naively unroll the loop for finding the start and end indices of each - // cell's data pointers in the array of boid indices - // - BIG DIFFERENCE: use the rearranged array index buffer to reshuffle all - // the particle data in the simulation array. - // CONSIDER WHAT ADDITIONAL BUFFERS YOU NEED - // - Perform velocity updates using neighbor search - // - Update positions - // - Ping-pong buffers as needed. THIS MAY BE DIFFERENT FROM BEFORE. + // TODO-2.3 - start by copying Boids::stepSimulationNaiveGrid + // Uniform Grid Neighbor search using Thrust sort on cell-coherent data. + // In Parallel: + // - Label each particle with its array index as well as its grid index. + // Use 2x width grids + // - Unstable key sort using Thrust. A stable sort isn't necessary, but you + // are welcome to do a performance comparison. + // - Naively unroll the loop for finding the start and end indices of each + // cell's data pointers in the array of boid indices + // - BIG DIFFERENCE: use the rearranged array index buffer to reshuffle all + // the particle data in the simulation array. + // CONSIDER WHAT ADDITIONAL BUFFERS YOU NEED + // - Perform velocity updates using neighbor search + // - Update positions + // - Ping-pong buffers as needed. THIS MAY BE DIFFERENT FROM BEFORE. + dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); + dim3 fullBlocksPerGridCells((gridCellCount + blockSize - 1) / blockSize); + + kernResetIntBuffer<<>>(gridCellCount, dev_gridCellStartIndices, -1); + kernResetIntBuffer<<>>(gridCellCount, dev_gridCellEndIndices, -1); + + kernComputeIndices<<>>(numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, + dev_pos, dev_particleArrayIndices, dev_particleGridIndices); + thrust::device_ptr dev_thrust_grid_indices(dev_particleGridIndices); + thrust::device_ptr dev_thrust_array_indices(dev_particleArrayIndices); + thrust::sort_by_key(dev_thrust_grid_indices, dev_thrust_grid_indices + numObjects, dev_thrust_array_indices); + + kernIdentifyCellStartEnd<<>>(numObjects, dev_particleGridIndices, dev_gridCellStartIndices, + dev_gridCellEndIndices); + //glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *coherent_pos, glm::vec3 *coherent_vel, int *particleArrayIndices + kernFlattenPosVel << > > (numObjects, dev_pos, dev_vel1, coherent_pos, coherent_vel, dev_particleArrayIndices); + kernUpdateVelNeighborSearchCoherent<<>>(numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, + gridCellWidth, dev_gridCellStartIndices, dev_gridCellEndIndices, coherent_pos, coherent_vel, dev_vel2); + std::swap(dev_pos, coherent_pos); + kernUpdatePos<<>>(numObjects, dt, dev_pos, dev_vel2); + std::swap(dev_vel1, dev_vel2); + } void Boids::endSimulation() { - cudaFree(dev_vel1); - cudaFree(dev_vel2); - cudaFree(dev_pos); - - // TODO-2.1 TODO-2.3 - Free any additional buffers here. + cudaFree(dev_vel1); + cudaFree(dev_vel2); + cudaFree(dev_pos); + + // TODO-2.1 TODO-2.3 - Free any additional buffers here. + cudaFree(dev_particleArrayIndices); + cudaFree(dev_particleGridIndices); + cudaFree(dev_gridCellStartIndices); + cudaFree(dev_gridCellEndIndices); } void Boids::unitTest() { - // LOOK-1.2 Feel free to write additional tests here. - - // test unstable sort - int *dev_intKeys; - int *dev_intValues; - int N = 10; - - std::unique_ptrintKeys{ new int[N] }; - std::unique_ptrintValues{ new int[N] }; - - intKeys[0] = 0; intValues[0] = 0; - intKeys[1] = 1; intValues[1] = 1; - intKeys[2] = 0; intValues[2] = 2; - intKeys[3] = 3; intValues[3] = 3; - intKeys[4] = 0; intValues[4] = 4; - intKeys[5] = 2; intValues[5] = 5; - intKeys[6] = 2; intValues[6] = 6; - intKeys[7] = 0; intValues[7] = 7; - intKeys[8] = 5; intValues[8] = 8; - intKeys[9] = 6; intValues[9] = 9; - - cudaMalloc((void**)&dev_intKeys, N * sizeof(int)); - checkCUDAErrorWithLine("cudaMalloc dev_intKeys failed!"); - - cudaMalloc((void**)&dev_intValues, N * sizeof(int)); - checkCUDAErrorWithLine("cudaMalloc dev_intValues failed!"); - - dim3 fullBlocksPerGrid((N + blockSize - 1) / blockSize); - - std::cout << "before unstable sort: " << std::endl; - for (int i = 0; i < N; i++) { - std::cout << " key: " << intKeys[i]; - std::cout << " value: " << intValues[i] << std::endl; - } - - // How to copy data to the GPU - cudaMemcpy(dev_intKeys, intKeys.get(), sizeof(int) * N, cudaMemcpyHostToDevice); - cudaMemcpy(dev_intValues, intValues.get(), sizeof(int) * N, cudaMemcpyHostToDevice); - - // Wrap device vectors in thrust iterators for use with thrust. - thrust::device_ptr dev_thrust_keys(dev_intKeys); - thrust::device_ptr dev_thrust_values(dev_intValues); - // LOOK-2.1 Example for using thrust::sort_by_key - thrust::sort_by_key(dev_thrust_keys, dev_thrust_keys + N, dev_thrust_values); - - // How to copy data back to the CPU side from the GPU - cudaMemcpy(intKeys.get(), dev_intKeys, sizeof(int) * N, cudaMemcpyDeviceToHost); - cudaMemcpy(intValues.get(), dev_intValues, sizeof(int) * N, cudaMemcpyDeviceToHost); - checkCUDAErrorWithLine("memcpy back failed!"); - - std::cout << "after unstable sort: " << std::endl; - for (int i = 0; i < N; i++) { - std::cout << " key: " << intKeys[i]; - std::cout << " value: " << intValues[i] << std::endl; - } - - // cleanup - cudaFree(dev_intKeys); - cudaFree(dev_intValues); - checkCUDAErrorWithLine("cudaFree failed!"); - return; + // LOOK-1.2 Feel free to write additional tests here. + + // test unstable sort + int *dev_intKeys; + int *dev_intValues; + int N = 10; + + std::unique_ptrintKeys{ new int[N] }; + std::unique_ptrintValues{ new int[N] }; + + intKeys[0] = 0; intValues[0] = 0; + intKeys[1] = 1; intValues[1] = 1; + intKeys[2] = 0; intValues[2] = 2; + intKeys[3] = 3; intValues[3] = 3; + intKeys[4] = 0; intValues[4] = 4; + intKeys[5] = 2; intValues[5] = 5; + intKeys[6] = 2; intValues[6] = 6; + intKeys[7] = 0; intValues[7] = 7; + intKeys[8] = 5; intValues[8] = 8; + intKeys[9] = 6; intValues[9] = 9; + + cudaMalloc((void**)&dev_intKeys, N * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_intKeys failed!"); + + cudaMalloc((void**)&dev_intValues, N * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_intValues failed!"); + + dim3 fullBlocksPerGrid((N + blockSize - 1) / blockSize); + + std::cout << "before unstable sort: " << std::endl; + for (int i = 0; i < N; i++) { + std::cout << " key: " << intKeys[i]; + std::cout << " value: " << intValues[i] << std::endl; + } + + // How to copy data to the GPU + cudaMemcpy(dev_intKeys, intKeys.get(), sizeof(int) * N, cudaMemcpyHostToDevice); + cudaMemcpy(dev_intValues, intValues.get(), sizeof(int) * N, cudaMemcpyHostToDevice); + + // Wrap device vectors in thrust iterators for use with thrust. + thrust::device_ptr dev_thrust_keys(dev_intKeys); + thrust::device_ptr dev_thrust_values(dev_intValues); + // LOOK-2.1 Example for using thrust::sort_by_key + thrust::sort_by_key(dev_thrust_keys, dev_thrust_keys + N, dev_thrust_values); + + // How to copy data back to the CPU side from the GPU + cudaMemcpy(intKeys.get(), dev_intKeys, sizeof(int) * N, cudaMemcpyDeviceToHost); + cudaMemcpy(intValues.get(), dev_intValues, sizeof(int) * N, cudaMemcpyDeviceToHost); + checkCUDAErrorWithLine("memcpy back failed!"); + + std::cout << "after unstable sort: " << std::endl; + for (int i = 0; i < N; i++) { + std::cout << " key: " << intKeys[i]; + std::cout << " value: " << intValues[i] << std::endl; + } + + // cleanup + cudaFree(dev_intKeys); + cudaFree(dev_intValues); + checkCUDAErrorWithLine("cudaFree failed!"); + return; } diff --git a/src/main.cpp b/src/main.cpp index b82c8c6..0c3b2e4 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -14,12 +14,12 @@ // LOOK-2.1 LOOK-2.3 - toggles for UNIFORM_GRID and COHERENT_GRID #define VISUALIZE 1 -#define UNIFORM_GRID 0 -#define COHERENT_GRID 0 +#define UNIFORM_GRID 1 +#define COHERENT_GRID 1 // LOOK-1.2 - change this to adjust particle count in the simulation -const int N_FOR_VIS = 5000; -const float DT = 0.2f; +const int N_FOR_VIS = 100000; +const float DT = 0.1f; /** * C main function. @@ -226,7 +226,7 @@ void initShaders(GLuint * program) { frame++; double time = glfwGetTime(); - if (time - timebase > 1.0) { + if (time - timebase > 10.0) { fps = frame / (time - timebase); timebase = time; frame = 0;