diff --git a/.vscode/launch.json b/.vscode/launch.json new file mode 100644 index 0000000..8c41d6d --- /dev/null +++ b/.vscode/launch.json @@ -0,0 +1,14 @@ +{ + // Use IntelliSense to learn about possible attributes. + // Hover to view descriptions of existing attributes. + // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387 + "version": "0.2.0", + "configurations": [ + { + "name": "CUDA C++: Launch", + "type": "cuda-gdb", + "request": "launch", + "program": "${workspaceFolder}/build/bin/cis5650_boids" + } + ] +} \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..479969b --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,74 @@ +{ + "files.associations": { + "cctype": "cpp", + "clocale": "cpp", + "cmath": "cpp", + "csignal": "cpp", + "cstdarg": "cpp", + "cstddef": "cpp", + "cstdio": "cpp", + "cstdlib": "cpp", + "cstring": "cpp", + "ctime": "cpp", + "cwchar": "cpp", + "cwctype": "cpp", + "array": "cpp", + "atomic": "cpp", + "strstream": "cpp", + "bit": "cpp", + "*.tcc": "cpp", + "bitset": "cpp", + "chrono": "cpp", + "compare": "cpp", + "complex": "cpp", + "concepts": "cpp", + "condition_variable": "cpp", + "cstdint": "cpp", + "deque": "cpp", + "list": "cpp", + "map": "cpp", + "set": "cpp", + "string": "cpp", + "unordered_map": "cpp", + "vector": "cpp", + "exception": "cpp", + "algorithm": "cpp", + "functional": "cpp", + "iterator": "cpp", + "memory": "cpp", + "memory_resource": "cpp", + "numeric": "cpp", + "optional": "cpp", + "random": "cpp", + "ratio": "cpp", + "string_view": "cpp", + "system_error": "cpp", + "tuple": "cpp", + "type_traits": "cpp", + "utility": "cpp", + "fstream": "cpp", + "initializer_list": "cpp", + "iomanip": "cpp", + "iosfwd": "cpp", + "iostream": "cpp", + "istream": "cpp", + "limits": "cpp", + "mutex": "cpp", + "new": "cpp", + "numbers": "cpp", + "ostream": "cpp", + "semaphore": "cpp", + "sstream": "cpp", + "stdexcept": "cpp", + "stop_token": "cpp", + "streambuf": "cpp", + "thread": "cpp", + "cfenv": "cpp", + "cinttypes": "cpp", + "typeindex": "cpp", + "typeinfo": "cpp", + "valarray": "cpp", + "variant": "cpp", + "*.ipp": "cpp" + } +} \ No newline at end of file diff --git a/README.md b/README.md index ee39093..175f380 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,110 @@ -**University of Pennsylvania, CIS 5650: GPU Programming and Architecture, -Project 1 - Flocking** +# Project 1 - Flocking +**University of Pennsylvania, CIS 5650: GPU Programming and Architecture** -* (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) +| ![](images/uniform-coherent_100000.png) | +|-----------------------------------------------------------| +| *Uniform-Coherent Grid Implementation with 100,000 Boids* | -### (TODO: Your README) +* Xinran Tao + * [LinkedIn](https://www.linkedin.com/in/xinran-tao/), [Personal Website](https://www.xinrantao.com/), [GitHub](https://github.com/theBoilingPoint). +* Tested on: + ``` + - OS: Windows 11, + - CPU: AMD Ryzen 7 7800X3D @ 4.20GHz x 8, + - RAM: 64GB, + - GPU: NVIDIA GeForce RTX 2080 Ti + ``` -Include screenshots, analysis, etc. (Remember, this is public, so don't put -anything here that you don't want to share with the world.) +## Visualisation +Below are the GIFs of different implementations with **100,000** boids. The simulations are run using *NSight Graphics* with VSync off. + +### Naive Implementation +![](images/naive_100000.gif) + +### Scattered Grid Implementation +![](images/scattered_100000.gif) + +### Coherent Grid Implementation +![](images/coherent_100000.gif) + +### Coherent Grid with Shared Memory Implementation +![](images/shared_100000.gif) + +## Performance Analysis +All experiments are conducted with +- Project release mode; +- VSync disabled; +- visualisation turned off. + +### FPS Change with Number of Boids + +#### Method +The block size is fixed at 128 and the cell width is fixed at 2 units for all experiments. + +Given that calculating the average FPS per `n` number of iterations (and printing them out in the terminal) decreases the FPS, the FPS listed below are the mean of the maximum and minimum FPS after the algorithm has converged (i.e. after running for 60s). + +#### Experiment Results + +##### Raw Data +| Number of Boids | FPS-Naive | FPS-Scattered Grid | FPS-Coherent Grid | +| ------------------ | --------- | -------------------| ----------------- | +| 5000 | 1009.0 | 1155.0 | 1200.0 | +| 10000 | 554.0 | 720.0 | 739.6 | +| 20000 | 277.5 | 419.0 | 427.4 | +| 40000 | 125.5 | 210.0 | 215.9 | +| 80000 | 32.6 | 80.0 | 80.0 | +| 160000 | 9.0 | 22.1 | 22.2 | + +##### Graph +![](images/analysis/numBoids_FPS.png) + +#### Results Analysis +For each implementation, the FPS decreases as the number of boids increases. All three implementations share the same FPS change pattern. The drop rate sof each FPS roughly follows the increase rate of the number of boids. More specifically, if the number of boids doubles, the FPS will be halved. + +Overall, the coherent grid implementation has the highest FPS, followed by the scattered grid implementation, and the naive implementation has the lowest FPS. Although as the number of boids increases to above 40000, the advantage of the coherent grid implementation over the scattered grid implementation becomes less significant. + +The general pattern of the FPS change with the number of boids is as expected. The naive implementation has the lowest FPS because it has the highest time complexity. The scattered grid implementation has a higher FPS than the naive implementation because it reduces the number of boids that need to be checked for each boid. The coherent grid implementation has the highest FPS because it reduces the number of memory accesses. + +Nonetheless, the FPS of all implementations decreases as the number of boids increases. In the naive implementation, each boid checks every other boid, leading to a quadratic increase in computations, while the grid-based implementations reduce this to a more manageable linear complexity by limiting checks to local grid cells. However, even with these optimizations, the workload still doubles as the number of boids doubles, leading to proportional FPS reductions. Additionally, increased memory traffic and GPU resource limits, such as bandwidth and cache usage, further contribute to the performance drop, as the system struggles to handle the larger number of boids efficiently. Thus, the FPS is inversely proportional to the number of boids, resulting in a halving of FPS as the boid count doubles. + +### FPS Change with Block Size + +#### Method +The number of boids is fixed at 5000 and the cell width is fixed at 2 units for all experiments. + +Once again, the FPS listed below are the mean of the maximum and minimum FPS after the algorithm has converged (i.e. after running for 60s). + +#### Experiment Results +Note that initialising `blockSize` to 2048 and above will result in a CUDA error. + +##### Raw Data +| Block Size | FPS-Naive | FPS-Scattered Grid | FPS-Coherent Grid | +| --------------| --------- | -------------------| ----------------- | +| 32 | 1033.9 | 1141.8 | 1178.0 | +| 64 | 1024.2 | 1146.7 | 1194.8 | +| 128 | 992.6 | 1166.8 | 1183.7 | +| 256 | 976.6 | 1168.4 | 1190.0 | +| 512 | 888.0 | 1123.4 | 1155.6 | +| 1024 | 568.5 | 873.4 | 899.2 | + +##### Graph +![](images/analysis/blockSize_FPS.png) + +#### Results Analysis +Overall, the coherent grid implementation still has the highest FPS, followed by the scattered grid implementation, and the naive implementation has the lowest FPS. The reason for this is the same as the previous analysis. + +At smaller block sizes (32–64), the GPU can efficiently handle computations, as the workload is well distributed across threads, and there is relatively little overhead in managing thread blocks. The higher FPS values observed at small block sizes are due to: +- **Optimal resource utilization**: The smaller block sizes allow the GPU to fully load its processing cores, maximizing parallel execution. +- **Efficient memory access**: At small block sizes, threads within a block are generally accessing contiguous memory, reducing memory latency and maximizing cache utilization. + +FPS shows marginal increases or plateaus at medium block sizes (128–256), depending on the implementation. The initial benefit of increasing block size diminishes as the block size grows. While a medium block size still allows good parallelization, the performance benefit tapers off due to: +- **Saturation of GPU resources**: As the block size increases, the number of threads per block also increases, and eventually, the GPU’s available resources (e.g., registers, shared memory) are fully utilized. Any further increase in block size provides diminishing returns. +- **Thread divergence and memory contention**: Larger block sizes may introduce thread divergence and increased contention for shared resources (e.g., memory bandwidth), which can limit performance gains. + +FPS starts to decrease significantly for all implementations at larger block sizes (512–1024). This performance drop is primarily due to the fact that as block size increases, the overhead of managing threads and memory becomes a bottleneck: +- **Inefficient resource usage**: Larger blocks may result in fewer blocks being executed concurrently on the GPU, leading to underutilization of GPU cores. The GPU cannot launch enough thread blocks simultaneously to hide memory latency, reducing overall throughput. +- **Increased memory pressure**: As block size grows, the amount of memory accessed per block increases, leading to higher memory traffic and contention for global memory bandwidth. This causes a significant drop in FPS as memory access becomes the limiting factor. +- **Worse memory access patterns**: Larger blocks tend to access memory more randomly or in scattered patterns, which can degrade the memory coalescing benefits seen at smaller block sizes. This increases memory latency and reduces performance. + +### FPS Change with Cell Width +For all three implementations, increasing the cell width from 2 to 3 units decreases the FPS. This is expected because each boid will have to check more cells in the grid, leading to more memory accesses and computations. \ No newline at end of file diff --git a/images/analysis/blockSize_FPS.png b/images/analysis/blockSize_FPS.png new file mode 100644 index 0000000..412187f Binary files /dev/null and b/images/analysis/blockSize_FPS.png differ diff --git a/images/analysis/numBoids_FPS.png b/images/analysis/numBoids_FPS.png new file mode 100644 index 0000000..95314f2 Binary files /dev/null and b/images/analysis/numBoids_FPS.png differ diff --git a/images/coherent.gif b/images/coherent.gif new file mode 100644 index 0000000..c812382 Binary files /dev/null and b/images/coherent.gif differ diff --git a/images/coherent_100000.gif b/images/coherent_100000.gif new file mode 100644 index 0000000..404495a Binary files /dev/null and b/images/coherent_100000.gif differ diff --git a/images/naive.gif b/images/naive.gif new file mode 100644 index 0000000..8cc6f6f Binary files /dev/null and b/images/naive.gif differ diff --git a/images/naive_100000.gif b/images/naive_100000.gif new file mode 100644 index 0000000..edc0c8d Binary files /dev/null and b/images/naive_100000.gif differ diff --git a/images/scattered.gif b/images/scattered.gif new file mode 100644 index 0000000..55745b6 Binary files /dev/null and b/images/scattered.gif differ diff --git a/images/scattered_100000.gif b/images/scattered_100000.gif new file mode 100644 index 0000000..b36988e Binary files /dev/null and b/images/scattered_100000.gif differ diff --git a/images/shared_100000.gif b/images/shared_100000.gif new file mode 100644 index 0000000..2eb730e Binary files /dev/null and b/images/shared_100000.gif differ diff --git a/images/uniform-coherent_100000.png b/images/uniform-coherent_100000.png new file mode 100644 index 0000000..27b69d9 Binary files /dev/null and b/images/uniform-coherent_100000.png differ diff --git a/src/kernel.cu b/src/kernel.cu index 74dffcb..46e93a8 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -17,6 +17,11 @@ #define checkCUDAErrorWithLine(msg) checkCUDAError(msg, __LINE__) +// smaller cell width = the neighbourhood distance +// larger cell width = 2 * the neighbourhood distance +#define USE_LARGER_CELL_WIDTH true +#define USE_SHARED_MEMORY true + /** * Check for CUDA errors; print and exit if there was a problem. */ @@ -37,7 +42,7 @@ void checkCUDAError(const char *msg, int line = -1) { *****************/ /*! Block size used for CUDA kernel launch. */ -#define blockSize 128 +#define blockSize 128 // Default is 128 // LOOK-1.2 Parameters for the boids algorithm. // These worked well in our reference implementation. @@ -85,11 +90,15 @@ 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_pos_coherent; +glm::vec3 *dev_vel1_coherent; +glm::vec3 *dev_vel2_coherent; // LOOK-2.1 - Grid parameters based on simulation parameters. // These are automatically computed for you in Boids::initSimulation int gridCellCount; int gridSideCount; +float neighborDistance; float gridCellWidth; float gridInverseCellWidth; glm::vec3 gridMinimum; @@ -141,7 +150,11 @@ void Boids::initSimulation(int 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. + // Don't forget to cudaFree in Boids::endSimulation. + + // Note: We don't pass dev_pos because this means we write the result (i.e. the GPU mem address we've allocated) to + // where dev_pos is pointed to, which is garbage. Instead, we pass &dev_pos, which means we write the result to + // wherever &dev_pos is pointed to, which is an initialised address (because we've initialised the variable and it has an address). cudaMalloc((void**)&dev_pos, N * sizeof(glm::vec3)); checkCUDAErrorWithLine("cudaMalloc dev_pos failed!"); @@ -157,7 +170,13 @@ void Boids::initSimulation(int N) { checkCUDAErrorWithLine("kernGenerateRandomPosArray failed!"); // LOOK-2.1 computing grid params - gridCellWidth = 2.0f * std::max(std::max(rule1Distance, rule2Distance), rule3Distance); + neighborDistance = std::max(std::max(rule1Distance, rule2Distance), rule3Distance); +#if USE_LARGER_CELL_WIDTH + gridCellWidth = 2.0f * neighborDistance; +#else + gridCellWidth = neighborDistance; +#endif + int halfSideCount = (int)(scene_scale / gridCellWidth) + 1; gridSideCount = 2 * halfSideCount; @@ -169,6 +188,31 @@ 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!"); + + dev_thrust_particleArrayIndices = thrust::device_pointer_cast(dev_particleArrayIndices); + + cudaMalloc((void**)&dev_particleGridIndices, N * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_particleGridIndices failed!"); + + dev_thrust_particleGridIndices = thrust::device_pointer_cast(dev_particleGridIndices); + + 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_pos_coherent, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_pos_coherent failed!"); + + cudaMalloc((void**)&dev_vel1_coherent, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_vel1_coherent failed!"); + + cudaMalloc((void**)&dev_vel2_coherent, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_vel2_coherent failed!"); + cudaDeviceSynchronize(); } @@ -233,7 +277,45 @@ __device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *po // 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 perceived_center = glm::vec3(0.0f, 0.0f, 0.0f); + glm::vec3 c = glm::vec3(0.0f, 0.0f, 0.0f); + glm::vec3 perceived_velocity = glm::vec3(0.0f, 0.0f, 0.0f); + + int num_neighbors_r1 = 0; + int num_neighbors_r3 = 0; + + for (int i = 0; i < N; i++) { + if (i != iSelf) { + float distance = glm::distance(pos[i], pos[iSelf]); + if (distance < rule1Distance) { + perceived_center += pos[i]; + num_neighbors_r1++; + } + + if (distance < rule2Distance) { + c -= pos[i] - pos[iSelf]; + } + + if (distance < rule3Distance) { + perceived_velocity += vel[i]; + num_neighbors_r3++; + } + } + } + + if (num_neighbors_r1 > 0) { + perceived_center /= num_neighbors_r1; + } + + if (num_neighbors_r3 > 0) { + perceived_velocity /= num_neighbors_r3; + } + + glm::vec3 v1 = (perceived_center - pos[iSelf]) * rule1Scale; + glm::vec3 v2 = c * rule2Scale; + glm::vec3 v3 = perceived_velocity * rule3Scale; + + return vel[iSelf] + v1 + v2 + v3; } /** @@ -245,6 +327,12 @@ __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 = (blockIdx.x * blockDim.x) + threadIdx.x; + if (index < N) { + glm::vec3 vel = computeVelocityChange(N, index, pos, vel1); + vel = glm::length(vel) > maxSpeed ? glm::normalize(vel) * maxSpeed : vel; + vel2[index] = vel; + } } /** @@ -289,6 +377,19 @@ __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) { + return; + } + + indices[index] = index; + + gridIndices[index] = gridIndex3Dto1D( + (pos[index].x - gridMin.x) * inverseCellWidth, + (pos[index].y - gridMin.y) * inverseCellWidth, + (pos[index].z - gridMin.z) * inverseCellWidth, + gridResolution + ); } // LOOK-2.1 Consider how this could be useful for indicating that a cell @@ -300,17 +401,27 @@ __global__ void kernResetIntBuffer(int N, int *intBuffer, int value) { } } -__global__ void kernIdentifyCellStartEnd(int N, int *particleGridIndices, - int *gridCellStartIndices, int *gridCellEndIndices) { +__global__ void kernIdentifyCellStartEnd( + int N, const int* sortedGridIdx, + int* cellStart, int* cellEnd) +{ // 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!" + // At this point, the particleGridIndices has already been sorted. + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i >= N) return; + int cur = sortedGridIdx[i]; + int prev = (i > 0) ? sortedGridIdx[i - 1] : -1; + int next = (i < N - 1) ? sortedGridIdx[i + 1] : -1; + if (i == 0 || cur != prev) cellStart[cur] = i; + if (i == N - 1 || cur != next) cellEnd[cur] = i; } __global__ void kernUpdateVelNeighborSearchScattered( int N, int gridResolution, glm::vec3 gridMin, - float inverseCellWidth, float cellWidth, + float inverseCellWidth, float neighborDistance, int *gridCellStartIndices, int *gridCellEndIndices, int *particleArrayIndices, glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) { @@ -322,25 +433,289 @@ __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 tid = blockIdx.x * blockDim.x + threadIdx.x; + if (tid >= N) return; + + // Process boids in the same order as the sorted grid arrays (better locality) + const int self = particleArrayIndices[tid]; + + // ----- compute integer cell of this boid ----- + const glm::vec3 rel = (pos[self] - gridMin) * inverseCellWidth; + const int ix = max(0, min(gridResolution - 1, (int)floorf(rel.x))); + const int iy = max(0, min(gridResolution - 1, (int)floorf(rel.y))); + const int iz = max(0, min(gridResolution - 1, (int)floorf(rel.z))); + + // ----- map [pos - R, pos + R] to index ranges (inclusive) ----- + const float R = neighborDistance; + + int ixMin = (int)floorf((pos[self].x - R - gridMin.x) * inverseCellWidth); + int ixMax = (int)floorf((pos[self].x + R - gridMin.x) * inverseCellWidth); + int iyMin = (int)floorf((pos[self].y - R - gridMin.y) * inverseCellWidth); + int iyMax = (int)floorf((pos[self].y + R - gridMin.y) * inverseCellWidth); + int izMin = (int)floorf((pos[self].z - R - gridMin.z) * inverseCellWidth); + int izMax = (int)floorf((pos[self].z + R - gridMin.z) * inverseCellWidth); + + ixMin = max(0, min(gridResolution - 1, ixMin)); + ixMax = max(0, min(gridResolution - 1, ixMax)); + iyMin = max(0, min(gridResolution - 1, iyMin)); + iyMax = max(0, min(gridResolution - 1, iyMax)); + izMin = max(0, min(gridResolution - 1, izMin)); + izMax = max(0, min(gridResolution - 1, izMax)); + + // ----- accumulators ----- + glm::vec3 perceived_center(0.0f); + glm::vec3 c(0.0f); + glm::vec3 perceived_velocity(0.0f); + int num_r1 = 0, num_r3 = 0; + + const float r1sq = rule1Distance * rule1Distance; + const float r2sq = rule2Distance * rule2Distance; + const float r3sq = rule3Distance * rule3Distance; + + // ----- iterate only the needed neighbor cells (inclusive ranges) ----- + for (int z = izMin; z <= izMax; ++z) { + for (int y = iyMin; y <= iyMax; ++y) { + for (int x = ixMin; x <= ixMax; ++x) { + const int cell1D = gridIndex3Dto1D(x, y, z, gridResolution); + + const int start = gridCellStartIndices[cell1D]; + if (start == -1) continue; // empty cell + + const int end = gridCellEndIndices[cell1D]; + + // Indirection for scattered layout: turn sorted slot -> original boid index + for (int b = start; b <= end; ++b) { + const int j = particleArrayIndices[b]; + if (j == self) continue; + + const glm::vec3 d = pos[j] - pos[self]; + const float dist2 = glm::dot(d, d); + + if (dist2 < r1sq) { perceived_center += pos[j]; ++num_r1; } + if (dist2 < r2sq) { c -= d; } + if (dist2 < r3sq) { perceived_velocity += vel1[j]; ++num_r3; } + } + } + } + } + + if (num_r1 > 0) perceived_center /= (float)num_r1; + if (num_r3 > 0) perceived_velocity /= (float)num_r3; + + glm::vec3 v1 = (perceived_center - pos[self]) * rule1Scale; + glm::vec3 v2 = c * rule2Scale; + glm::vec3 v3 = perceived_velocity * rule3Scale; + + glm::vec3 v = vel1[self] + v1 + v2 + v3; + float speed = glm::length(v); + if (speed > maxSpeed) v *= (maxSpeed / speed); + + vel2[self] = v; +} + +__global__ void kernGenerateCoherentPosVal(int N, int *dev_particleArrayIndices, glm::vec3 *dev_pos, glm::vec3 *dev_vel1, glm::vec3 *dev_pos_coherent, glm::vec3 *dev_vel1_coherent) { + int index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index >= N) { + return; + } + + int array_index = dev_particleArrayIndices[index]; + dev_pos_coherent[index] = dev_pos[array_index]; + dev_vel1_coherent[index] = dev_vel1[array_index]; +} + +__global__ void kernRestorePosValFromCoherent(int N, int *dev_particleArrayIndices, glm::vec3 *dev_pos, glm::vec3 *dev_vel2, glm::vec3 *dev_pos_coherent, glm::vec3 *dev_vel2_coherent) { + int index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index >= N) { + return; + } + + int array_index = dev_particleArrayIndices[index]; + dev_pos[array_index] = dev_pos_coherent[index]; + dev_vel2[array_index] = dev_vel2_coherent[index]; } __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 + float inverseCellWidth, float neighborDistance, + const int *gridCellStartIndices, const int *gridCellEndIndices, + const glm::vec3 *pos, const glm::vec3 *vel1, glm::vec3 *vel2) { + + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i >= N) return; + + // Axis-wise min/max indices from [pos - R, pos + R] + auto toCell = [&](float x, float gmin) { + return (int)floorf((x - gmin) * inverseCellWidth); + }; + int ixMin = max(0, toCell(pos[i].x - neighborDistance, gridMin.x)); + int ixMax = min(gridResolution - 1, toCell(pos[i].x + neighborDistance, gridMin.x)); + int iyMin = max(0, toCell(pos[i].y - neighborDistance, gridMin.y)); + int iyMax = min(gridResolution - 1, toCell(pos[i].y + neighborDistance, gridMin.y)); + int izMin = max(0, toCell(pos[i].z - neighborDistance, gridMin.z)); + int izMax = min(gridResolution - 1, toCell(pos[i].z + neighborDistance, gridMin.z)); + + // ----- accumulators ----- + glm::vec3 perceived_center(0.0f); + glm::vec3 c(0.0f); + glm::vec3 perceived_velocity(0.0f); + int num_r1 = 0, num_r3 = 0; + + const float R1sq = rule1Distance * rule1Distance; + const float R2sq = rule2Distance * rule2Distance; + const float R3sq = rule3Distance * rule3Distance; + + // Iterate only the needed cells (1–8 when cellWidth=2R) + for (int z = izMin; z <= izMax; ++z) { + for (int y = iyMin; y <= iyMax; ++y) { + for (int x = ixMin; x <= ixMax; ++x) { + int cell1D = gridIndex3Dto1D(x, y, z, gridResolution); + int start = gridCellStartIndices[cell1D]; + if (start == -1) continue; + int end = gridCellEndIndices[cell1D]; + + // NOTE: coherent layout → direct indices b + for (int b = start; b <= end; ++b) { + if (b == i) continue; + + glm::vec3 d = pos[b] - pos[i]; + float dist2 = glm::dot(d, d); + + if (dist2 < R1sq) { perceived_center += pos[b]; ++num_r1; } + if (dist2 < R2sq) { c -= d; } + if (dist2 < R3sq) { perceived_velocity += vel1[b]; ++num_r3; } + } + } + } + } + + if (num_r1 > 0) perceived_center /= (float)num_r1; + if (num_r3 > 0) perceived_velocity /= (float)num_r3; + + glm::vec3 v1 = (perceived_center - pos[i]) * rule1Scale; + glm::vec3 v2 = c * rule2Scale; + glm::vec3 v3 = perceived_velocity * rule3Scale; + + glm::vec3 v = vel1[i] + v1 + v2 + v3; + float speed = glm::length(v); + if (speed > maxSpeed) v = v * (maxSpeed / speed); + + vel2[i] = v; +} + +__global__ void kernUpdateVelCoherentShared( + int gridResolution, + glm::vec3 gridMin, + float invCellWidth, + float neighborDistance, + const int* __restrict__ cellStart, + const int* __restrict__ cellEnd, + const glm::vec3* __restrict__ posCo, // coherent + const glm::vec3* __restrict__ velCo, // coherent + glm::vec3* __restrict__ outVelCo // coherent +) { + // We use float4 because it’s the most GPU-friendly way to store 3D vectors for CUDA kernels, especially with shared-memory tiling. + extern __shared__ float4 shmem[]; // 2*blockDim.x float4's + float4* shPos = shmem; + float4* shVel = shmem + blockDim.x; + + const int cell = blockIdx.x; + const int start = cellStart[cell]; + if (start == -1) return; + const int end = cellEnd[cell]; + + // This block processes boids in [start..end] + // Loop in case the cell has more boids than blockDim.x + // (increment by blockDim.x because boid with index < multiple of it will be processed by other threads) + for (int selfIdx = start + threadIdx.x; selfIdx <= end; selfIdx += blockDim.x) { + // Per-thread accumulators (registers) + glm::vec3 selfPos = posCo[selfIdx]; + glm::vec3 v1acc(0); + glm::vec3 v2acc(0); + glm::vec3 v3acc(0); + int n1 = 0, n3 = 0; + + // Compute neighbor-cell index ranges from [pos+-R] + const float R = neighborDistance; + auto toCell = [&](float x, float gmin) { + return (int)floorf((x - gmin) * invCellWidth); + }; + + int ixMin = max(0, toCell(selfPos.x - R, gridMin.x)); + int ixMax = min(gridResolution - 1, toCell(selfPos.x + R, gridMin.x)); + int iyMin = max(0, toCell(selfPos.y - R, gridMin.y)); + int iyMax = min(gridResolution - 1, toCell(selfPos.y + R, gridMin.y)); + int izMin = max(0, toCell(selfPos.z - R, gridMin.z)); + int izMax = min(gridResolution - 1, toCell(selfPos.z + R, gridMin.z)); + + const float r1sq = rule1Distance * rule1Distance; + const float r2sq = rule2Distance * rule2Distance; + const float r3sq = rule3Distance * rule3Distance; + + // Loop all neighbor cells (inclusive ranges) + for (int z = izMin; z <= izMax; ++z) { + for (int y = iyMin; y <= iyMax; ++y) { + for (int x = ixMin; x <= ixMax; ++x) { + int nCell = gridIndex3Dto1D(x, y, z, gridResolution); + int ns = cellStart[nCell]; + if (ns == -1) continue; + int ne = cellEnd[nCell]; + + // Tile the neighbor cell into shared memory + // Note that this is run for every neighbour cell. The loop makes sure that all neighboour boids are loaded. + for (int tile = ns; tile <= ne; tile += blockDim.x) { + int j = tile + threadIdx.x; + + // Cooperative load (coalesced) + if (j <= ne) { + shPos[threadIdx.x] = make_float4(posCo[j].x, posCo[j].y, posCo[j].z, 0); + shVel[threadIdx.x] = make_float4(velCo[j].x, velCo[j].y, velCo[j].z, 0); + } + __syncthreads(); + + int count = min(blockDim.x, ne - tile + 1); + // Consume the tile from shared memory + #pragma unroll + for (int t = 0; t < count; ++t) { + int idx = tile + t; + if (idx == selfIdx) continue; + + float3 d; + d.x = shPos[t].x - selfPos.x; + d.y = shPos[t].y - selfPos.y; + d.z = shPos[t].z - selfPos.z; + float dist2 = d.x*d.x + d.y*d.y + d.z*d.z; + + if (dist2 < r1sq) { v1acc.x += shPos[t].x; v1acc.y += shPos[t].y; v1acc.z += shPos[t].z; ++n1; } + if (dist2 < r2sq) { v2acc.x -= d.x; v2acc.y -= d.y; v2acc.z -= d.z; } + if (dist2 < r3sq) { v3acc.x += shVel[t].x; v3acc.y += shVel[t].y; v3acc.z += shVel[t].z; ++n3; } + } + __syncthreads(); + } + } + } + } + + if (n1 > 0) { v1acc /= (float)n1; v1acc -= selfPos; } + if (n3 > 0) { v3acc /= (float)n3; } + + // Scales + v1acc *= rule1Scale; + v2acc *= rule2Scale; + v3acc *= rule3Scale; + + glm::vec3 v = velCo[selfIdx]; + v += v1acc + v2acc + v3acc; + + // Clamp speed + float speed2 = glm::dot(v, v); + if (speed2 > maxSpeed * maxSpeed) { + float inv = maxSpeed * rsqrtf(speed2); + v *= inv; + } + + outVelCo[selfIdx] = glm::vec3(v.x, v.y, v.z); + } } /** @@ -349,6 +724,15 @@ __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 +748,22 @@ 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 fullBlocksPerCell((gridCellCount + blockSize - 1) / blockSize); + + kernComputeIndices << > >(numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, dev_pos, dev_particleArrayIndices, dev_particleGridIndices); + + thrust::sort_by_key(dev_thrust_particleGridIndices, dev_thrust_particleGridIndices + numObjects, dev_thrust_particleArrayIndices); + + kernResetIntBuffer << > >(gridCellCount, dev_gridCellStartIndices, -1); + + kernIdentifyCellStartEnd << > >(numObjects, dev_particleGridIndices, dev_gridCellStartIndices, dev_gridCellEndIndices); + + kernUpdateVelNeighborSearchScattered << > >(numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, neighborDistance, 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) { @@ -382,6 +782,55 @@ 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 fullBlocksPerCell((gridCellCount + blockSize - 1) / blockSize); + + kernComputeIndices << > >(numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, dev_pos, dev_particleArrayIndices, dev_particleGridIndices); + + thrust::sort_by_key(dev_thrust_particleGridIndices, dev_thrust_particleGridIndices + numObjects, dev_thrust_particleArrayIndices); + + kernResetIntBuffer << > >(gridCellCount, dev_gridCellStartIndices, -1); + + kernIdentifyCellStartEnd << > >(numObjects, dev_particleGridIndices, dev_gridCellStartIndices, dev_gridCellEndIndices); + + kernGenerateCoherentPosVal << > >(numObjects, dev_particleArrayIndices, dev_pos, dev_vel1, dev_pos_coherent, dev_vel1_coherent); + +#if USE_SHARED_MEMORY + // grid dim defines the number of blocks and block dim defines the number of threads + dim3 grid(gridCellCount); + int block = 32; + size_t shmem = 2 * block * sizeof(float4); + kernUpdateVelCoherentShared <<>> ( + gridSideCount, + gridMinimum, + gridInverseCellWidth, + neighborDistance, + dev_gridCellStartIndices, + dev_gridCellEndIndices, + dev_pos_coherent, + dev_vel1_coherent, + dev_vel2_coherent + ); +#else + kernUpdateVelNeighborSearchCoherent << > >( + numObjects, + gridSideCount, + gridMinimum, + gridInverseCellWidth, + neighborDistance, + dev_gridCellStartIndices, + dev_gridCellEndIndices, + dev_pos_coherent, + dev_vel1_coherent, + dev_vel2_coherent + ); +#endif + + kernUpdatePos << > >(numObjects, dt, dev_pos_coherent, dev_vel2_coherent); + + kernRestorePosValFromCoherent << > >(numObjects, dev_particleArrayIndices, dev_pos, dev_vel2, dev_pos_coherent, dev_vel2_coherent); + + std::swap(dev_vel1, dev_vel2); } void Boids::endSimulation() { @@ -390,6 +839,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_gridCellEndIndices); + cudaFree(dev_pos_coherent); + cudaFree(dev_vel1_coherent); + cudaFree(dev_vel2_coherent); } void Boids::unitTest() { diff --git a/src/main.cpp b/src/main.cpp index fe657ed..4cd62ac 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -16,12 +16,17 @@ // ================ // LOOK-2.1 LOOK-2.3 - toggles for UNIFORM_GRID and COHERENT_GRID +// NOTE: Apart from setting this to 0, remember to also turn off VSync in NVIDIA Control Panel. +// VSync caps your frame rate to the monitor’s refresh (e.g., 60 Hz), which hides the true performance of +// the CUDA boid simulation and also adds extra synchronization you don’t want when profiling. #define VISUALIZE 1 -#define UNIFORM_GRID 0 -#define COHERENT_GRID 0 +#define LIMIT_FRAMES 0 // turn this on for profiling if you want to limit the number of frames +#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 maxFrames = 10000; // the maximum number of frames to profile +const int N_FOR_VIS = 100000; const float DT = 0.2f; /** @@ -219,14 +224,21 @@ void initShaders(GLuint * program) { double fps = 0; double timebase = 0; int frame = 0; + int totalFrames = 0; Boids::unitTest(); // LOOK-1.2 We run some basic example code to make sure // your CUDA development setup is ready to go. while (!glfwWindowShouldClose(window)) { glfwPollEvents(); - frame++; + totalFrames++; +#if LIMIT_FRAMES + if (totalFrames > maxFrames) { + break; + } +#endif + double time = glfwGetTime(); if (time - timebase > 1.0) { @@ -259,6 +271,7 @@ void initShaders(GLuint * program) { glfwSwapBuffers(window); #endif } + glfwDestroyWindow(window); glfwTerminate(); }