diff --git a/README.md b/README.md index ee39093..46873dc 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,63 @@ **University of Pennsylvania, CIS 5650: 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) +- Dineth Meegoda + - [LinkedIn](https://www.linkedin.com/in/dinethmeegoda/), [personal website](https://www.dinethmeegoda.com). +- Tested on: Windows 10 Pro, Ryzen 9 5900X 12 Core @ 3.7GHz 32GB, RTX 3070 8GB -### (TODO: Your README) +## Summary -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/demo.gif) + +This project implements a flocking simulation based on Craig Reynold's Boids Rules and Algorithm. This project was optimized from the original naive method to include uniform grid based neighbor searching and coherent memory access. + +This was done with the following approaches: + +- **Naive**: Reynold's three rules are followed, and each boid checks every other boid to see how it influences itself (even if it is beyond its search distance). + +- **Uniform**: Each boid is placed into a grid of uniform cell widths. As a result, boids only search the nearest cells that its search distance overlaps. This increases performance by a significant amount. + +- **Coherent**: The memory for each position and velocity for boids is placed contigiously in memory as opposed to scattered. This allows for a considerable margin of performance to be gained. + +## Performance Analysis + +_Notes: The Performance based on FPS was calculated by averaging the reported Frame Count over 12 seconds, taking measurements every second. The first two counts of frames were removed to allow the simulation to stablize._ + +### For each implementation, how does changing the number of boids affect performance? Why do you think this is? + +- **Naive**: Changing the number of boids heavily impacts performance since the algorithm has O(n^2) complexity. Each boid must check every other boid, which causes poor scaling. + +- **Uniform**: This approach has a significant increase in performance over the naive approach since the number of checks is reduced with using the grid. As a result, the falloff of performance as the boids count increases is more linear, rather than exponential. + +- **Coherent** This approach is a bit faster since it performs an additional optimization on top of the existing grid method. The increasing of boids still causes a linear decrease in FPS, however the performance is better and causes a less drastic drop off. The algorithm is still tracking a similar amount of neighbors, so algorithmic effiency has not changed much, but the memory management optimization causes some performance boost. + +
+ + + + + +
+ Image 1 +

Performance Graph with Visualization

+ 		+ Image 2 +

Performance Graph with No Visualization

+
+
+ +### For each implementation, how does changing the block count and block size affect performance? Why do you think this is? + +For all the approaches, the performance increases as the block size grows to about 64, and then plataeus and slightly falls off from there. The difference in performance between the implementation comes from the algorithmic discrepencies of computing the simulation. When the block size is under 64, the GPU seems to be underutilized. However, when the block size grows larger, some performance does drop off since there is a trade off between increasing the warps run on a block and resource availability (memory, scheduling with stalls, etc.). This hardware limit causes the limit on performance that is able to be gained by increasing the block size. + +![](images/blockSize.png) + +### 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? + +I did expect the coherent uniform grid to give a performance boost since the positions and velocities of each boid was contigiously in memory as opposed to being scattered in memory with the previous implementation. The largest performance inhibitor is reading and writing to memory, and improving memory latency by making them easier to access does give a significant performance boost. + +### Did changing cell width and checking 27 vs 8 neighboring cells affect performance? Why or why not? + +More cell checks generally resulted in more work, and less performance, but smaller cells means that there are less boids in each cell, which means less checks. However, there still must be a good amount of cells in order to be able to effectively use the cell/grid method as an optimization. + +Generally, I experienced the greatest amount of performance when the cell width was a little less than twice the maximum search distance. This means that there were more than 8 cell checks, but the smaller amount of cells to check means less boids had to be checked than necessary. The real amount that we must limit is the number of boids that we check that are not within the search distance. In fact, as the amount of boids increased, and the simulation became more dense, performance was saved more by the greater amount of cell checks (of smaller cells). diff --git a/images/blockSize.png b/images/blockSize.png new file mode 100644 index 0000000..9375549 Binary files /dev/null and b/images/blockSize.png differ diff --git a/images/boidsNoVis.png b/images/boidsNoVis.png new file mode 100644 index 0000000..2d507c0 Binary files /dev/null and b/images/boidsNoVis.png differ diff --git a/images/boidsVis.png b/images/boidsVis.png new file mode 100644 index 0000000..d0260c6 Binary files /dev/null and b/images/boidsVis.png differ diff --git a/images/demo.gif b/images/demo.gif new file mode 100644 index 0000000..fe12fe5 Binary files /dev/null and b/images/demo.gif differ diff --git a/src/kernel.cu b/src/kernel.cu index 74dffcb..9974fac 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -86,6 +86,9 @@ int *dev_gridCellEndIndices; // to this cell? // TODO-2.3 - consider what additional buffers you might need to reshuffle // the position and velocity data to be coherent within cells. +glm::vec3* dev_coherentPos; +glm::vec3* dev_coherentVel; + // LOOK-2.1 - Grid parameters based on simulation parameters. // These are automatically computed for you in Boids::initSimulation int gridCellCount; @@ -169,6 +172,26 @@ void Boids::initSimulation(int N) { gridMinimum.z -= halfGridWidth; // TODO-2.1 TODO-2.3 - Allocate additional buffers here. + + cudaMalloc((void**)&dev_particleArrayIndices, N * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_particleArrayIndices failed!"); + + cudaMalloc((void**)&dev_particleGridIndices, N * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_particleGridIndices failed!"); + + cudaMalloc((void**)&dev_gridCellStartIndices, gridCellCount * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_gridCellStartIndices failed!"); + + cudaMalloc((void**)&dev_gridCellEndIndices, gridCellCount * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_gridCellEndIndices failed!"); + + cudaMalloc((void**)&dev_coherentPos, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_coherentPos failed!"); + + cudaMalloc((void**)&dev_coherentVel, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_coherentVel failed!"); + + cudaDeviceSynchronize(); } @@ -230,10 +253,49 @@ 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); + + glm::vec3 result(0.0f); + glm::vec3 perceived_center(0.0f), c( 0.0f), perceived_velocity(0.0f); + + int r1_num_neighbors = 0, r3_num_neighbors = 0; + + for (int i = 0; i < N; i++) { + + float distance = glm::distance(pos[i], pos[iSelf]); + + // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves + if (i != iSelf && distance < rule1Distance) { + perceived_center += pos[i]; + r1_num_neighbors++; + } + + // Rule 2: boids try to stay a distance d away from each other + if (i != iSelf && distance < rule2Distance) { + c -= (pos[i] - pos[iSelf]); + } + + // Rule 3: boids try to match the speed of surrounding boids + if (i != iSelf && distance < rule3Distance) { + perceived_velocity += vel[i]; + r3_num_neighbors++; + } + } + + // Finalize Rule 1 and Rule 3 + if (r1_num_neighbors > 0) { + perceived_center /= r1_num_neighbors; + result += (perceived_center - pos[iSelf]) * rule1Scale; + } + + if (r3_num_neighbors > 0) { + perceived_velocity /= r3_num_neighbors; + result += perceived_velocity * rule3Scale; + } + + // Finalize Rule 2 + result += c * rule2Scale; + + return result; } /** @@ -245,6 +307,15 @@ __global__ void kernUpdateVelocityBruteForce(int N, glm::vec3 *pos, // 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) { + glm::vec3 velocity = vel1[index] + computeVelocityChange(N, index, pos, vel1); + float speed = glm::length(velocity); + if (speed > maxSpeed) { + velocity = glm::normalize(velocity) * maxSpeed; + } + vel2[index] = velocity; + } } /** @@ -289,6 +360,12 @@ __global__ void kernComputeIndices(int N, int gridResolution, // - 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 + int index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index < N) { + indices[index] = index; + glm::vec3 gridIndex3D = glm::floor((pos[index] - gridMin) * inverseCellWidth); + gridIndices[index] = gridIndex3Dto1D(gridIndex3D.x, gridIndex3D.y, gridIndex3D.z, gridResolution); + } } // LOOK-2.1 Consider how this could be useful for indicating that a cell @@ -306,6 +383,38 @@ __global__ void kernIdentifyCellStartEnd(int N, int *particleGridIndices, // 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 index = threadIdx.x + (blockIdx.x * blockDim.x); + + if (index < N) { + // If at the start, set the start index as current + if (index == 0) { + gridCellStartIndices[particleGridIndices[index]] = index; + } + // If at the end, set the end index as current + else if (index == N - 1) { + gridCellEndIndices[particleGridIndices[index]] = index; + } + else { + // If the current grid index is different from the previous one, set the end index + // of the previous grid index as the current index - 1 and the start index of the + // current grid index as the current index + if (particleGridIndices[index] != particleGridIndices[index - 1]) { + gridCellEndIndices[particleGridIndices[index - 1]] = index - 1; + gridCellStartIndices[particleGridIndices[index]] = index; + } + } + } + +} + +__global__ void kernReshuffleBuffer(int N, int* particleArrayIndices, glm::vec3* pos, glm::vec3* vel, glm::vec3* coherentPos, glm::vec3* coherentVel) { + int index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index < N) { + int coherentIndex = particleArrayIndices[index]; + coherentPos[index] = pos[coherentIndex]; + coherentVel[index] = vel[coherentIndex]; + } } __global__ void kernUpdateVelNeighborSearchScattered( @@ -322,6 +431,80 @@ __global__ void kernUpdateVelNeighborSearchScattered( // - 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 index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index < N) { + glm::vec3 perceived_center(0.0f), c(0.0f), perceived_velocity(0.0f), result(0.0f); + int num_neighbors_rule1 = 0, num_neighbors_rule3 = 0; + + glm::vec3 position = pos[index], gridPos = glm::floor((position - gridMin) * inverseCellWidth); + + // Iterate over the 3x3x3 grid of cells around the current cell (upto 8 cells) + for (int x = imax(gridPos.x - 1, 0); x < imin(gridPos.x + 1, gridResolution - 1); x++) { + for (int y = imax(gridPos.y - 1, 0); y < imin(gridPos.y + 1, gridResolution - 1); y++) { + for (int z = imax(gridPos.z - 1, 0); z < imin(gridPos.z + 1, gridResolution - 1); z++) { + + int currGridIndex = gridIndex3Dto1D(x, y, z, gridResolution); + int currStartIndex = gridCellStartIndices[currGridIndex]; + int currEndIndex = gridCellEndIndices[currGridIndex]; + + // If no boids are in the current cell, continue + if (currStartIndex == -1 || currEndIndex == -1) { + continue; + } + + + for (int i = currStartIndex; i <= currEndIndex; i++) { + + int boidIndex = particleArrayIndices[i]; + + if (boidIndex != index) { + float distance = glm::distance(pos[boidIndex], position); + + // Rule 1 + if (distance < rule1Distance) { + perceived_center += pos[boidIndex]; + num_neighbors_rule1++; + } + // Rule 2 + if (distance < rule2Distance) { + c -= (pos[boidIndex] - position); + } + // Rule 3 + if (distance < rule3Distance) { + perceived_velocity += vel1[boidIndex]; + num_neighbors_rule3++; + } + } + } + } + } + } + + // Finalize Rule 1 and Rule 3 + + if (num_neighbors_rule1 > 0) { + perceived_center /= num_neighbors_rule1; + result += (perceived_center - position) * rule1Scale; + } + + if (num_neighbors_rule3 > 0) { + perceived_velocity /= num_neighbors_rule3; + result += perceived_velocity * rule3Scale; + } + + // Finalize Rule 2 + result += c * rule2Scale; + + result += vel1[index]; + + // Clamp the speed + float speed = glm::length(result); + if (speed > maxSpeed) { + result = glm::normalize(result) * maxSpeed; + } + + vel2[index] = result; + } } __global__ void kernUpdateVelNeighborSearchCoherent( @@ -341,6 +524,79 @@ __global__ void kernUpdateVelNeighborSearchCoherent( // - 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 index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index < N) { + glm::vec3 perceived_center(0.0f), c(0.0f), perceived_velocity(0.0f), result(0.0f); + int num_neighbors_rule1 = 0, num_neighbors_rule3 = 0; + + glm::vec3 position = pos[index], gridPos = glm::floor((position - gridMin) * inverseCellWidth); + + // Iterate over the 3x3x3 grid of cells around the current cell (upto 8 cells) + for (int x = imax(gridPos.x - 1, 0); x < imin(gridPos.x + 1, gridResolution - 1); x++) { + for (int y = imax(gridPos.y - 1, 0); y < imin(gridPos.y + 1, gridResolution - 1); y++) { + for (int z = imax(gridPos.z - 1, 0); z < imin(gridPos.z + 1, gridResolution - 1); z++) { + + int currGridIndex = gridIndex3Dto1D(x, y, z, gridResolution); + int currStartIndex = gridCellStartIndices[currGridIndex]; + int currEndIndex = gridCellEndIndices[currGridIndex]; + + // If no boids are in the current cell, continue + if (currStartIndex == -1 || currEndIndex == -1) { + continue; + } + + + for (int i = currStartIndex; i <= currEndIndex; i++) { + + if (i != index) { + float distance = glm::distance(pos[i], position); + + // Rule 1 + if (distance < rule1Distance) { + perceived_center += pos[i]; + num_neighbors_rule1++; + } + // Rule 2 + if (distance < rule2Distance) { + c -= (pos[i] - position); + } + // Rule 3 + if (distance < rule3Distance) { + perceived_velocity += vel1[i]; + num_neighbors_rule3++; + } + } + } + } + } + } + + // Finalize Rule 1 and Rule 3 + + if (num_neighbors_rule1 > 0) { + perceived_center /= num_neighbors_rule1; + result += (perceived_center - position) * rule1Scale; + } + + if (num_neighbors_rule3 > 0) { + perceived_velocity /= num_neighbors_rule3; + result += perceived_velocity * rule3Scale; + } + + // Finalize Rule 2 + result += c * rule2Scale; + + result += vel1[index]; + + // Clamp the speed + float speed = glm::length(result); + if (speed > maxSpeed) { + result = glm::normalize(result) * maxSpeed; + } + + vel2[index] = result; + } } /** @@ -349,6 +605,14 @@ __global__ void kernUpdateVelNeighborSearchCoherent( 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 + dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); + kernUpdateVelocityBruteForce <<>> (numObjects, dev_pos, dev_vel1, dev_vel2); + checkCUDAErrorWithLine("kernUpdateVelocityBruteForce failed!"); + + kernUpdatePos <<>> (numObjects, dt, dev_pos, dev_vel2); + checkCUDAErrorWithLine("kernUpdatePos failed!"); + + std::swap(dev_vel1, dev_vel2); } void Boids::stepSimulationScatteredGrid(float dt) { @@ -364,6 +628,38 @@ void Boids::stepSimulationScatteredGrid(float dt) { // - Perform velocity updates using neighbor search // - Update positions // - Ping-pong buffers as needed + + dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); + dim3 fullCellsPerGrid((gridCellCount + blockSize - 1) / blockSize); + + // Compute indices + kernComputeIndices <<>> (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, dev_pos, dev_particleArrayIndices, dev_particleGridIndices); + checkCUDAErrorWithLine("kernComputeIndices failed!"); + + // Unstable key sort + thrust::device_ptr dev_thrust_keys(dev_particleGridIndices); + thrust::device_ptr dev_thrust_values(dev_particleArrayIndices); + thrust::sort_by_key(dev_thrust_keys, dev_thrust_keys + numObjects, dev_thrust_values); + + // Reset buffers to -1 + kernResetIntBuffer <<>> (gridCellCount, dev_gridCellStartIndices, -1); + checkCUDAErrorWithLine("kernResetIntBuffer for Start failed!"); + kernResetIntBuffer <<>> (gridCellCount, dev_gridCellEndIndices, -1); + checkCUDAErrorWithLine("kernResetIntBuffer for End failed!"); + + // Identify cell start and end indices + kernIdentifyCellStartEnd <<>> (numObjects, dev_particleGridIndices, dev_gridCellStartIndices, dev_gridCellEndIndices); + checkCUDAErrorWithLine("kernIdentifyCellStartEnd failed!"); + + // Update velocity using neighbor search + kernUpdateVelNeighborSearchScattered <<>> (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, gridCellWidth, dev_gridCellStartIndices, dev_gridCellEndIndices, dev_particleArrayIndices, dev_pos, dev_vel1, dev_vel2); + checkCUDAErrorWithLine("kernUpdateVelNeighborSearchScattered failed!"); + + // Update positions + kernUpdatePos <<>> (numObjects, dt, dev_pos, dev_vel2); + checkCUDAErrorWithLine("kernUpdatePos failed!"); + // Ping Pong buffers + std::swap(dev_vel1, dev_vel2); } void Boids::stepSimulationCoherentGrid(float dt) { @@ -382,6 +678,48 @@ void Boids::stepSimulationCoherentGrid(float dt) { // - 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 fullCellsPerGrid((gridCellCount + blockSize - 1) / blockSize); + + // Compute indices + kernComputeIndices << > > (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, dev_pos, dev_particleArrayIndices, dev_particleGridIndices); + checkCUDAErrorWithLine("kernComputeIndices failed!"); + + // Unstable key sort + thrust::device_ptr dev_thrust_keys(dev_particleGridIndices); + thrust::device_ptr dev_thrust_values(dev_particleArrayIndices); + thrust::sort_by_key(dev_thrust_keys, dev_thrust_keys + numObjects, dev_thrust_values); + + // Reset buffers to -1 + kernResetIntBuffer <<>> (gridCellCount, dev_gridCellStartIndices, -1); + checkCUDAErrorWithLine("kernResetIntBuffer for Start failed!"); + kernResetIntBuffer <<>> (gridCellCount, dev_gridCellEndIndices, -1); + checkCUDAErrorWithLine("kernResetIntBuffer for End failed!"); + + // Identify cell start and end indices + kernIdentifyCellStartEnd <<>> (numObjects, dev_particleGridIndices, dev_gridCellStartIndices, dev_gridCellEndIndices); + checkCUDAErrorWithLine("kernIdentifyCellStartEnd failed!"); + + // Reshuffle buffers + kernReshuffleBuffer <<>> (numObjects, dev_particleArrayIndices, dev_pos, dev_vel1, dev_coherentPos, dev_coherentVel); + checkCUDAErrorWithLine("kernReshuffleBuffer failed!"); + + // Update velocity using neighbor search + kernUpdateVelNeighborSearchCoherent <<>> (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, gridCellWidth, dev_gridCellStartIndices, dev_gridCellEndIndices, dev_coherentPos, dev_coherentVel, dev_vel2); + checkCUDAErrorWithLine("kernUpdateVelNeighborSearchCoherent failed!"); + + // Update positions + kernUpdatePos <<>> (numObjects, dt, dev_coherentPos, dev_vel2); + checkCUDAErrorWithLine("kernUpdatePos failed!"); + + // Ping Pong buffers + std::swap(dev_pos, dev_coherentPos); + std::swap(dev_vel1, dev_coherentVel); + std::swap(dev_vel1, dev_vel2); + + + } void Boids::endSimulation() { @@ -390,6 +728,13 @@ void Boids::endSimulation() { 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_gridCellStartIndices); + + cudaFree(dev_coherentPos); + cudaFree(dev_coherentVel); } void Boids::unitTest() { diff --git a/src/main.cpp b/src/main.cpp index fe657ed..6c80628 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -17,11 +17,11 @@ // 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 int N_FOR_VIS = 40000; const float DT = 0.2f; /**