diff --git a/README.md b/README.md index ee39093..8306d08 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,44 @@ -**University of Pennsylvania, CIS 5650: GPU Programming and Architecture, -Project 1 - Flocking** +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) +**University of Pennsylvania, CIS 5650: GPU Programming and Architecture, Project 0** -### (TODO: Your README) +* ADITHYA RAJEEV + * [LinkedIn](https://www.linkedin.com/in/adithyar262/) +* Tested on: Windows 11, i7 13th Gen @ 2.40GHz 16GB, GeForce RTX 4050 8GB (Personal) -Include screenshots, analysis, etc. (Remember, this is public, so don't put -anything here that you don't want to share with the world.) +# Boid Simulation + +![](images/Flocking_GIF.gif) + +# Performance Analysis + +## 1. FPS vs No. of Boids + +![](images/Boid_FPS_Table.png) + +![](images/Boid_FPS_Plot.png) + +## 2. FPS vs Block Size + +![](images/Block_Size_FPS_Table.png) + +![](images/Block_Size_FPS_Plot.png) + +# Questions + +## 1. For each implementation, how does changing the number of boids affect performance? Why do you think this is? +The increase in the number of boids leads to a decrease in frames per second (FPS). The inverse relationship between the number of boids and FPS is a natural consequence of the increased computational load and resource requirements as the simulation scales up. This effect is observed across both naive and optimized implementations, although optimized versions may scale better with increasing boid counts. + +## 2. For each implementation, how does changing the block count and block size affect performance? Why do you think this is? +Changing the block size and number of blocks does not significantly affect the FPS (Frames Per Second) in the boids simulation. Changing the block size and number of blocks is more about optimizing how the work is distributed on the GPU rather than changing the fundamental amount or speed of parallel operations being performed. As a result, these changes typically don't lead to significant variations in FPS for the boids simulation. + +## 3. 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? +The coherent uniform grid implementation demonstrated a significant performance improvement over the scattered grid approach. In the coherent grid method, we eliminate the need for an additional matrix to redirect to the position and velocity array indices for boids. This removal of indirection means one less memory lookup operation (or warp) to complete for each data access. By directly accessing the boid data, we reduce the overall number of memory operations, leading to faster execution times. The coherent grid approach also organizes boid data in a way that results in more contiguous and sequential memory accesses when retrieving position and velocity information. +The combination of reduced indirection and more coherent memory access patterns results in more efficient use of the GPU's memory subsystem. This leads to reduced latency, increased throughput, and ultimately, the observed performance improvement in the coherent uniform grid implementation. + +## 4. 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 number of neighboring cells checked to 27 (3x3x3 grid) improves performance in the boid simulation. This quantization of the search space leads to checking smaller volumes around each boid, often resulting in fewer overall boid comparisons. The more compact search area enhances spatial locality and cache usage, leading to more efficient memory access patterns. Despite checking more cells, the total volume searched is smaller, reducing computational overhead. This demonstrates that a finer-grained approach to neighbor searching can yield significant performance benefits by decreasing workload and enhancing memory efficiency. + + + diff --git a/images/Block_Size_FPS_Plot.png b/images/Block_Size_FPS_Plot.png new file mode 100644 index 0000000..bb2fdb6 Binary files /dev/null and b/images/Block_Size_FPS_Plot.png differ diff --git a/images/Block_Size_FPS_Table.png b/images/Block_Size_FPS_Table.png new file mode 100644 index 0000000..8969768 Binary files /dev/null and b/images/Block_Size_FPS_Table.png differ diff --git a/images/Boid_FPS_Plot.png b/images/Boid_FPS_Plot.png new file mode 100644 index 0000000..3009d49 Binary files /dev/null and b/images/Boid_FPS_Plot.png differ diff --git a/images/Boid_FPS_Table.png b/images/Boid_FPS_Table.png new file mode 100644 index 0000000..cd95cd2 Binary files /dev/null and b/images/Boid_FPS_Table.png differ diff --git a/images/Flocking_GIF.gif b/images/Flocking_GIF.gif new file mode 100644 index 0000000..3568c24 Binary files /dev/null and b/images/Flocking_GIF.gif differ diff --git a/images/Flocking_Video.mp4 b/images/Flocking_Video.mp4 new file mode 100644 index 0000000..816b8b1 Binary files /dev/null and b/images/Flocking_Video.mp4 differ diff --git a/src/kernel.cu b/src/kernel.cu index 74dffcb..55a31dd 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -20,15 +20,15 @@ /** * 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); +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); } - fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString(err)); - exit(EXIT_FAILURE); - } } @@ -66,22 +66,22 @@ dim3 threadsPerBlock(blockSize); // Consider why you would need two velocity buffers in a simulation where each // boid cares about its neighbors' velocities. // These are called ping-pong buffers. -glm::vec3 *dev_pos; -glm::vec3 *dev_vel1; -glm::vec3 *dev_vel2; +glm::vec3* dev_pos; +glm::vec3* dev_vel1; +glm::vec3* dev_vel2; // LOOK-2.1 - these are NOT allocated for you. You'll have to set up the thrust // pointers on your own too. // For efficient sorting and the uniform grid. These should always be parallel. -int *dev_particleArrayIndices; // What index in dev_pos and dev_velX represents this particle? -int *dev_particleGridIndices; // What grid cell is this particle in? +int* dev_particleArrayIndices; // What index in dev_pos and dev_velX represents this particle? +int* dev_particleGridIndices; // What grid cell is this particle in? // needed for use with thrust thrust::device_ptr dev_thrust_particleArrayIndices; thrust::device_ptr dev_thrust_particleGridIndices; -int *dev_gridCellStartIndices; // What part of dev_particleArrayIndices belongs -int *dev_gridCellEndIndices; // to this cell? +int* dev_gridCellStartIndices; // What part of dev_particleArrayIndices belongs +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. @@ -99,13 +99,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,63 +113,76 @@ __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)); } /** * LOOK-1.2 - This is a basic CUDA kernel. * 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; - } +__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; + } } /** * 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**)&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. + 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!"); + + // Initialize Thrust device pointers + dev_thrust_particleArrayIndices = thrust::device_pointer_cast(dev_particleArrayIndices); + dev_thrust_particleGridIndices = thrust::device_pointer_cast(dev_particleGridIndices); + + cudaDeviceSynchronize(); } @@ -180,42 +193,42 @@ 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); +__global__ void kernCopyPositionsToVBO(int N, glm::vec3* pos, float* vbo, float s_scale) { + 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); +__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; - } + 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); +void Boids::copyBoidsToVBO(float* vbodptr_positions, float* vbodptr_velocities) { + 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(); } @@ -229,47 +242,84 @@ void Boids::copyBoidsToVBO(float *vbodptr_positions, float *vbodptr_velocities) * Compute the new velocity on the body with index `iSelf` due to the `N` boids * 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); +__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 + glm::vec3 selfPos = pos[iSelf]; + glm::vec3 selfVel = vel[iSelf]; + + glm::vec3 cohesion(0.0f); + glm::vec3 separation(0.0f); + glm::vec3 alignment(0.0f); + int neighborCount = 0; + + for (int i = 0; i < N; i++) { + if (i == iSelf) continue; + + float dist = glm::distance(pos[i], selfPos); + + if (dist < rule1Distance) { + cohesion += pos[i]; + alignment += vel[i]; + neighborCount++; + + if (dist < rule2Distance) { + separation -= (pos[i] - selfPos); + } + } + } + + glm::vec3 velocityChange(0.0f); + + if (neighborCount > 0) { + cohesion = (cohesion / float(neighborCount) - selfPos) * rule1Scale; + alignment = (alignment / float(neighborCount) - selfVel) * rule3Scale; + velocityChange = cohesion + alignment + (separation * rule2Scale); + } + + return selfVel + velocityChange; } /** * TODO-1.2 implement basic flocking * 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? +__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? + int i = threadIdx.x + (blockIdx.x * blockDim.x); + if (i >= N) { + return; + } + vel2[i] = glm::clamp(computeVelocityChange(N, i, pos, vel1), -maxSpeed, maxSpeed); } /** * LOOK-1.2 Since this is pretty trivial, we implemented it for you. * 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; +__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; + // 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; + 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; + pos[index] = thisPos; } // LOOK-2.1 Consider this method of computing a 1D index from a 3D grid index. @@ -279,179 +329,482 @@ __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) { + 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 + int index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index < N) { + // Compute the grid cell coordinates for this boid + glm::vec3 relativePos = pos[index] - gridMin; + int gridX = static_cast(relativePos.x * inverseCellWidth); + int gridY = static_cast(relativePos.y * inverseCellWidth); + int gridZ = static_cast(relativePos.z * inverseCellWidth); + + // Clamp the grid coordinates to ensure they're within the grid + gridX = glm::clamp(gridX, 0, gridResolution - 1); + gridY = glm::clamp(gridY, 0, gridResolution - 1); + gridZ = glm::clamp(gridZ, 0, gridResolution - 1); + + // - Label each boid with the index of its grid cell. + int gridIndex = gridIndex3Dto1D(gridX, gridY, gridZ, gridResolution); + gridIndices[index] = gridIndex; + + // - Set up a parallel array of integer indices as pointers to the actual + // boid data in pos and vel1/vel2 + indices[index] = index; + } + return; } // 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; - } +__global__ void kernResetIntBuffer(int N, int* intBuffer, int 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!" +__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 index = threadIdx.x + (blockIdx.x * blockDim.x); + + if (index < N) { + int currentCellIndex = particleGridIndices[index]; + + // Check if this is the start of a new cell + if (index == 0 || currentCellIndex != particleGridIndices[index - 1]) { + gridCellStartIndices[currentCellIndex] = index; + } + + // Check if this is the end of a cell + if (index == N - 1 || currentCellIndex != particleGridIndices[index + 1]) { + gridCellEndIndices[currentCellIndex] = index + 1; + } + } + + return; + } __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) { + // 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 index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index >= N) return; + + int boidIndex = particleArrayIndices[index]; + glm::vec3 boidPos = pos[boidIndex]; + glm::vec3 cellCoord = (boidPos - gridMin) * inverseCellWidth; + glm::ivec3 cellIndex = glm::floor(cellCoord); + + glm::ivec3 searchStart = glm::max(cellIndex - 1, glm::ivec3(0)); + glm::ivec3 searchEnd = glm::min(cellIndex + 1, glm::ivec3(gridResolution - 1)); + + for (int i = 0; i < 3; ++i) { + if (cellCoord[i] > cellIndex[i] + 0.5f) { + searchStart[i] = cellIndex[i]; + } + else { + searchEnd[i] = cellIndex[i]; + } + } + + glm::vec3 velChange(0.0f); + glm::vec3 perceivedCenter(0.0f); + glm::vec3 separation(0.0f); + glm::vec3 perceivedVelocity(0.0f); + int numNeighbors1 = 0, numNeighbors3 = 0; + + for (int x = searchStart.x; x <= searchEnd.x; ++x) { + for (int y = searchStart.y; y <= searchEnd.y; ++y) { + for (int z = searchStart.z; z <= searchEnd.z; ++z) { + int cellIdx = gridIndex3Dto1D(x, y, z, gridResolution); + if (gridCellStartIndices[cellIdx] == -1) continue; + + for (int i = gridCellStartIndices[cellIdx]; i <= gridCellEndIndices[cellIdx]; ++i) { + int neighborIndex = particleArrayIndices[i]; + if (neighborIndex == boidIndex) continue; + + glm::vec3 neighborPos = pos[neighborIndex]; + float dist = glm::distance(neighborPos, boidPos); + + if (dist < rule1Distance) { + perceivedCenter += neighborPos; + ++numNeighbors1; + } + if (dist < rule2Distance) { + separation -= (neighborPos - boidPos); + } + if (dist < rule3Distance) { + perceivedVelocity += vel1[neighborIndex]; + ++numNeighbors3; + } + } + } + } + } + + if (numNeighbors1 > 0) { + perceivedCenter /= float(numNeighbors1); + velChange += (perceivedCenter - boidPos) * rule1Scale; + } + velChange += separation * rule2Scale; + if (numNeighbors3 > 0) { + perceivedVelocity /= float(numNeighbors3); + velChange += perceivedVelocity * rule3Scale; + } + + glm::vec3 newVel = vel1[boidIndex] + velChange; + vel2[boidIndex] = glm::length(newVel) > maxSpeed ? maxSpeed * glm::normalize(newVel) : 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 index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index >= N) return; + + glm::vec3 boidPos = pos[index]; + glm::vec3 cellCoord = (boidPos - gridMin) * inverseCellWidth; + glm::ivec3 cellIndex = glm::floor(cellCoord); + + glm::ivec3 searchStart = glm::max(cellIndex - 1, glm::ivec3(0)); + glm::ivec3 searchEnd = glm::min(cellIndex + 1, glm::ivec3(gridResolution - 1)); + + for (int i = 0; i < 3; ++i) { + if (cellCoord[i] > cellIndex[i] + 0.5f) { + searchStart[i] = cellIndex[i]; + } + else { + searchEnd[i] = cellIndex[i]; + } + } + + glm::vec3 velChange(0.0f); + glm::vec3 perceivedCenter(0.0f); + glm::vec3 separation(0.0f); + glm::vec3 perceivedVelocity(0.0f); + int numNeighbors1 = 0, numNeighbors3 = 0; + + // Reorder the cell checking to maximize memory coherence + for (int z = searchStart.z; z <= searchEnd.z; ++z) { + for (int y = searchStart.y; y <= searchEnd.y; ++y) { + for (int x = searchStart.x; x <= searchEnd.x; ++x) { + int cellIdx = gridIndex3Dto1D(x, y, z, gridResolution); + if (gridCellStartIndices[cellIdx] == -1) continue; + + for (int i = gridCellStartIndices[cellIdx]; i <= gridCellEndIndices[cellIdx]; ++i) { + if (i == index) continue; + + glm::vec3 neighborPos = pos[i]; + float dist = glm::distance(neighborPos, boidPos); + + if (dist < rule1Distance) { + perceivedCenter += neighborPos; + ++numNeighbors1; + } + if (dist < rule2Distance) { + separation -= (neighborPos - boidPos); + } + if (dist < rule3Distance) { + perceivedVelocity += vel1[i]; + ++numNeighbors3; + } + } + } + } + } + + if (numNeighbors1 > 0) { + perceivedCenter /= float(numNeighbors1); + velChange += (perceivedCenter - boidPos) * rule1Scale; + } + velChange += separation * rule2Scale; + if (numNeighbors3 > 0) { + perceivedVelocity /= float(numNeighbors3); + velChange += perceivedVelocity * rule3Scale; + } + + glm::vec3 newVel = vel1[index] + velChange; + vel2[index] = glm::length(newVel) > maxSpeed ? maxSpeed * glm::normalize(newVel) : 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 + int n = (numObjects + blockSize - 1) / blockSize; + // TODO-1.2 - use the kernels you wrote to step the simulation forward in time. + kernUpdateVelocityBruteForce << > > (numObjects, dev_pos, dev_vel1, dev_vel2); + kernUpdatePos << > > (numObjects, dt, dev_pos, dev_vel2); + // TODO-1.2 ping-pong the velocity buffers + 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. + + dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); + + // - label each particle with its array index as well as its grid index. + // Use 2x width grids. + kernComputeIndices << > > ( + numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, + dev_pos, dev_particleArrayIndices, dev_particleGridIndices); + + // - Unstable key sort using Thrust. A stable sort isn't necessary, but you + // are welcome to do a performance comparison. + thrust::sort_by_key(dev_thrust_particleGridIndices, + dev_thrust_particleGridIndices + numObjects, + dev_thrust_particleArrayIndices); + + kernResetIntBuffer << > > ( + gridCellCount, dev_gridCellStartIndices, -1); + kernResetIntBuffer << > > ( + gridCellCount, dev_gridCellEndIndices, -1); + + // - Naively unroll the loop for finding the start and end indices of each + // cell's data pointers in the array of boid indices + kernIdentifyCellStartEnd << > > ( + numObjects, dev_particleGridIndices, dev_gridCellStartIndices, dev_gridCellEndIndices); + // - Perform velocity updates using neighbor search + kernUpdateVelNeighborSearchScattered << > > ( + numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, gridCellWidth, + dev_gridCellStartIndices, dev_gridCellEndIndices, dev_particleArrayIndices, + dev_pos, dev_vel1, dev_vel2); + // - Update positions + kernUpdatePos << > > ( + numObjects, dt, dev_pos, dev_vel2); + + // - Ping-pong buffers as needed + std::swap(dev_vel1, dev_vel2); + +} + +__global__ void kernReorderBoidsData( + int N, int* particleArrayIndices, + glm::vec3* pos, glm::vec3* vel1, + glm::vec3* coherent_pos, glm::vec3* coherent_vel1) { + + int index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index >= N) return; + + int sortedIndex = particleArrayIndices[index]; + coherent_pos[index] = pos[sortedIndex]; + coherent_vel1[index] = vel1[sortedIndex]; +} + +__global__ void kernRestoreBoidsOrder( + int N, int* particleArrayIndices, + glm::vec3* coherent_pos, glm::vec3* coherent_vel2, + glm::vec3* pos, glm::vec3* vel1) { + + int index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index >= N) return; + + int originalIndex = particleArrayIndices[index]; + pos[originalIndex] = coherent_pos[index]; + vel1[originalIndex] = coherent_vel2[index]; +} + +__global__ void kernReshuffleData( + int N, int* particleArrayIndices, + glm::vec3* pos_in, glm::vec3* vel_in, + glm::vec3* pos_out, glm::vec3* vel_out) { + + int index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index >= N) return; + + int sortedIndex = particleArrayIndices[index]; + pos_out[index] = pos_in[sortedIndex]; + vel_out[index] = vel_in[sortedIndex]; } 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); + + // Label each particle with its array index and grid index + kernComputeIndices << > > ( + numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, + dev_pos, dev_particleArrayIndices, dev_particleGridIndices); + + // Sort particles by grid index + thrust::sort_by_key(dev_thrust_particleGridIndices, + dev_thrust_particleGridIndices + numObjects, + dev_thrust_particleArrayIndices); + + // Reset grid cell start and end indices + kernResetIntBuffer << > > ( + gridCellCount, dev_gridCellStartIndices, -1); + kernResetIntBuffer << > > ( + gridCellCount, dev_gridCellEndIndices, -1); + + // Identify start and end indices for each grid cell + kernIdentifyCellStartEnd << > > ( + numObjects, dev_particleGridIndices, dev_gridCellStartIndices, dev_gridCellEndIndices); + + // BIG DIFFERENCE: Reshuffle particle data based on sorted indices + // Allocate temporary buffers for coherent data + glm::vec3* dev_pos_coherent; + glm::vec3* dev_vel1_coherent; + cudaMalloc(&dev_pos_coherent, numObjects * sizeof(glm::vec3)); + cudaMalloc(&dev_vel1_coherent, numObjects * sizeof(glm::vec3)); + + // Kernel to reshuffle data + kernReshuffleData << > > ( + numObjects, dev_particleArrayIndices, + dev_pos, dev_vel1, + dev_pos_coherent, dev_vel1_coherent); + + // Update velocities using coherent data + kernUpdateVelNeighborSearchCoherent << > > ( + numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, gridCellWidth, + dev_gridCellStartIndices, dev_gridCellEndIndices, + dev_pos_coherent, dev_vel1_coherent, dev_vel2); + + // Update positions + kernUpdatePos << > > ( + numObjects, dt, dev_pos_coherent, dev_vel2); + + // Ping-pong buffers + std::swap(dev_vel1_coherent, dev_vel2); + + // Copy coherent data back to original buffers + cudaMemcpy(dev_pos, dev_pos_coherent, numObjects * sizeof(glm::vec3), cudaMemcpyDeviceToDevice); + cudaMemcpy(dev_vel1, dev_vel1_coherent, numObjects * sizeof(glm::vec3), cudaMemcpyDeviceToDevice); + + // Free temporary buffers + cudaFree(dev_pos_coherent); + cudaFree(dev_vel1_coherent); } void Boids::endSimulation() { - cudaFree(dev_vel1); - cudaFree(dev_vel2); - cudaFree(dev_pos); + 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); - // TODO-2.1 TODO-2.3 - Free any additional buffers here. } 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 fe657ed..5ca608a 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -17,11 +17,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 int N_FOR_VIS = 100000; +//const int N_FOR_VIS = 100; const float DT = 0.2f; /**