diff --git a/README.md b/README.md
index ee39093..500f41b 100644
--- a/README.md
+++ b/README.md
@@ -1,11 +1,141 @@
-**University of Pennsylvania, CIS 5650: GPU Programming and Architecture,
-Project 1 - Flocking**
+# **Project 1 - CUDA 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)
+* Jefferson Koumba Moussadji Lu
+ * [LinkedIn](https://www.linkedin.com/in/-jeff-koumba-0b356721b/)
+* Tested on: Personal Laptop, Windows 11 Home, Intel(R) Core(TM) i9-14900HX @ 2.22GHz @ 24 Cores @ 32GB RAM, Nvidia GeForce RTX 4090 @ 16 GB @ SM 8.9
-### (TODO: Your README)
+
+ 
+
Welcome to the world of Boids!
+
-Include screenshots, analysis, etc. (Remember, this is public, so don't put
-anything here that you don't want to share with the world.)
+## Overview
+
+This project implements a massively parallel flocking simulation based on Reynolds' Boids algorithm using CUDA and the GPU. The implementation features three progressive optimization levels:
+
+- Naive brute-force neighbor search: Each boid compares against all others with time complexity of (O(N²)).
+- Uniform Grid spatial partitioning with scattered memory access: Bin boids into cells, sort by cell index, and only check nearby cells’ occupants.
+- Coherent Grid with optimized memory layout: Like scattered, but reorders boid data in memory so neighbors are contiguous; improves memory coalescing and cache locality.
+- Grid-Looping Optimization (Extra Credit) with dynamic search boundaries: removes the hard‑coded “8 or 27 neighbor cells” assumption and loops only over cells that actually intersect a boid’s interaction radius.
+
+
+ 
+
Visualization of 10,000 boids, 100,000 boids, and 1,000,000 boids from left to right and top to bottom
+
+
+## Techniques
+
+### Naive
+
+One CUDA thread per boid; each thread loops over all other boids and accumulates cohesion, separation, and alignment contributions per rule radius/scale
+
+### Uniform grid scattered
+Compute each boid’s grid cell index and sort (key = cell, value = boid index) using Thrust. Build gridCellStartIndices/gridCellEndIndices so each cell maps to a contiguous range in a sorted index buffer. Neighbor checks are limited to nearby cells (8 or 27, depending on grid width).
+
+### Uniform grid (coherent)
+After sorting, reorder positions/velocities into sorted order so cell members are also contiguous in the position/velocity arrays. This removes the extra indirection and yields coalesced memory access during the neighbor search kernel.
+
+### Grid‑Looping Optimization (extra credit)
+Instead of hard‑coding a list of neighbor cells (8 or 27), compute (per boid) the min/max cell indices along each axis that actually intersect the sphere of influence and loop only those cells. That is, bound your loop by (minCellX…maxCellX) × (minCellY…maxCellY) × (minCellZ…maxCellZ) for the current boid. It’s flexible and avoids unnecessary checks of far corner cells.
+
+
+## Performance Analysis
+
+### Impact of the Number of Boids on Performance
+
+
+ 
+
Framerate change with increasing number of boids for naive, scattered uniform grid, and coherent uniform grid with visualization
+
+
+ 
+
Framerate change with increasing number of boids for naive, scattered uniform grid, and coherent uniform grid without visualization
+
+| # Boids | Naive (Viz) | Naive (No Viz) | Uniform (Viz) | Uniform (No Viz) | Coherent (Viz) | **Coherent (No Viz)** |
+| --------: | ----------: | -------------: | ------------: | ---------------: | -------------: | --------------------: |
+| 10,000 | 240.1 | 559.9 | 240.0 | 1684.4 | 240.0 | **2516.6** |
+| 50,000 | 60.5 | 77.7 | 240.0 | 802.4 | 240.0 | **1631.6** |
+| 100,000 | 19.1 | 20.0 | 240.0 | 333.2 | 240.0 | **1594.7** |
+| 500,000 | 0.9 | 1.1 | 28.1 | 30.8 | **238.0** | **289.5** |
+| 1,000,000 | 0.2 | 0.2 | 6.8 | 7.7 | **75.6** | **88.1** |
+
+From the table, we can deduct that for performance, the lower number Boids is, the higher the FPS with or without visualization.
+- Naive declines sharply with N (quadratic work). At 10k, several hundred FPS; at 100k, tens; at 1M, ~0.2 FPS.
+- Uniform (scattered) scales ~linearly with N (sort + limited neighbor checks). Hundreds to thousands of FPS for 10k–100k; tens for 1M.
+- Coherent is the fastest at all N because coalesced memory reduces global memory traffic in the neighbor kernel.
+
+With the rendering enabled, low counts are often vsync‑limited ~240 FPS on my setup. At higher counts, the compute becomes the bottleneck. We notice that coherent uniform remains much faster than the naive, and significantly faster than scattered uniform at large N.
+
+Coherent uniform achieves the highest simulation FPS at every size while the naive collapses past ~50k boids. Turning viz off reveals true simulation throughput.
+
+This makes sense since neighbor candidates are pruned by spatial binning (vs. O(N²)), and coherent access makes each candidate check cheaper.
+
+### Impact of the Block Size and Block Count on Performance
+
+
+ 
+
Framerate change with increasing block size
+
+| Block size | 16 | 32 | 64 | 128 | **256** | 512 | 1024 |
+| ---------: | -----: | -----: | -----: | -----: | ---------: | -----: | -----: |
+| FPS | 1171.4 | 1508.7 | 1600.2 | 1556.8 | **1638.4** | 1515.9 | 1440.9 |
+
+We can notice that the FPS peak around 256 threads/block on my RTX 4090 Laptop GPU and the framerate drops when the block size are too small or too big.
+
+This can be explained by the fact that very small blocks under‑utilize warps and increase scheduling overhead while very large blocks reduce SM occupancy and may exhaust registers/shared memory, hurting latency hiding.
+
+### Coherent uniform grid performance
+
+Coherent uniform grid improved performance at medium and large N, which is what I expected. For example, we can see that at 10,000 boids, we have 2516.6 vs 1684.4 FPS, and at 100,000 boids we have 1594.7 vs 333.2 FPS.
+
+By reordering positions and velocities so that boids in the same cell are contiguous, the neighbor kernel reads memory in a coalesced way. This reduces the number of global memory transactions per warp. The reorder costs extra work, so at very small N the benefit can be small or negative, but at scale the improvement dominates.
+
+### 27 vs 8 neighboring cells
+
+Using neighborhood‑sized cells (27 neighbors) did not slow things, often it was comparable or slightly faster. Although there are more cells to visit, each cell holds far fewer boids, so the total candidates can be lower than with 8 larger cells. The trade‑off depends on spatial distribution. With relatively uniform flocks, 27‑cell searches can reduce wasted distance checks while with highly clustered flocks, benefits shrink.
+
+Performance depends on how many neighbors you actually evaluate, not just how many cells you loop over.
+
+## Grid-Looping Optimization (Extra-credit)
+
+### Goal
+The goal is to dynamically determine (minCell…maxCell) per axis that intersect the influence radius and iterate only those cells—no hard‑coded “8 or 27” list.
+
+See ```Boids::stepSimulationGridLoopOptimized``` and ```__global__ void kernUpdateVelNeighborSearchCoherentOptimized``` in kernel.cu.
+
+To enable it, in main.cpp set:
+```
+#define UNIFORM_GRID 1
+#define COHERENT_GRID 1
+#define GRID_LOOP_OPTIMIZATION 1
+```
+
+
+ 
+
Framerate change with increasing number of boids for naive, scattered uniform grid, and coherent uniform grid without visualization
+
+| # Boids | **Grid‑Looping Optimized (FPS)** |
+| --------: | -------------------------------: |
+| 10,000 | **2562.9** |
+| 50,000 | **1723.8** |
+| 100,000 | **2163.9** |
+| 500,000 | **525.5** |
+| 1,000,000 | **168.4** |
+
+
+When the radius intersects fewer than the full 3×3×3 set, we skip empty/out‑of‑range cells, reducing per‑boid work while keeping correctness
+
+We can observe that Grid-Looping outperforms standard uniform coherent.
+
+At 10k–100k, Grid‑Looping does better than standard coherent grid by a healthy margin. Since fewer cells actually intersect, we have fewer candidate checks.
+
+At 500k–1M, the gains narrow. Most nearby cells contain boids, sort/reorder costs dominate, and total neighbor counts are high regardless. It still performs better than uniform coherent grid at 1M (~168 FPS vs ~88 FPS) with a smaller gap.
+
+The grid-looping optimization demonstrated that dynamic neighborhood bounds can prune unnecessary checks and outperform the standard coherent approach at moderate scales, while still remaining competitive at very large counts.
+
+## Conclusion
+
+These experiments highlight a central theme in GPU programming: algorithmic complexity alone is not enough. Data layout, memory access patterns, and occupancy tuning are just as critical for unlocking the GPU’s full potential. By carefully combining spatial partitioning with memory coherence, I achieved simulations of up to 1,000,000 boids at interactive framerates on modern hardware.
+
+This project not only deepened my understanding of CUDA and GPU architecture but also reinforced the importance of profiling, experimentation, and iterative optimization when designing high-performance systems. It was rewarding to see abstract concepts like warp divergence, coalescing, and occupancy translate directly into measurable FPS gains and scalable, visually engaging simulations.
\ No newline at end of file
diff --git a/images/Bloid_1000000.gif b/images/Bloid_1000000.gif
new file mode 100644
index 0000000..58c1fd6
Binary files /dev/null and b/images/Bloid_1000000.gif differ
diff --git a/images/Framerate_change_block_size.png b/images/Framerate_change_block_size.png
new file mode 100644
index 0000000..01b18f9
Binary files /dev/null and b/images/Framerate_change_block_size.png differ
diff --git a/images/Framerate_change_visualization.png b/images/Framerate_change_visualization.png
new file mode 100644
index 0000000..13461ec
Binary files /dev/null and b/images/Framerate_change_visualization.png differ
diff --git a/images/Framerate_change_without_visualization.png b/images/Framerate_change_without_visualization.png
new file mode 100644
index 0000000..82f089e
Binary files /dev/null and b/images/Framerate_change_without_visualization.png differ
diff --git a/images/Grid_Beginning_2.gif b/images/Grid_Beginning_2.gif
new file mode 100644
index 0000000..5b783d7
Binary files /dev/null and b/images/Grid_Beginning_2.gif differ
diff --git a/images/Nvidia Nsight System Analysis.png b/images/Nvidia Nsight System Analysis.png
new file mode 100644
index 0000000..a25743a
Binary files /dev/null and b/images/Nvidia Nsight System Analysis.png differ
diff --git a/images/Nvidis Nsight system Timeline.png b/images/Nvidis Nsight system Timeline.png
new file mode 100644
index 0000000..78994b9
Binary files /dev/null and b/images/Nvidis Nsight system Timeline.png differ
diff --git a/images/Presentation.gif b/images/Presentation.gif
new file mode 100644
index 0000000..78aa8bd
Binary files /dev/null and b/images/Presentation.gif differ
diff --git a/images/bloid_10000.gif b/images/bloid_10000.gif
new file mode 100644
index 0000000..c69d67b
Binary files /dev/null and b/images/bloid_10000.gif differ
diff --git a/images/grid_loop_optimization.png b/images/grid_loop_optimization.png
new file mode 100644
index 0000000..ccbb1b1
Binary files /dev/null and b/images/grid_loop_optimization.png differ
diff --git a/src/kernel.cu b/src/kernel.cu
index 7149917..1737ab1 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -1,4 +1,4 @@
-#define GLM_FORCE_CUDA
+#define GLM_FORCE_CUDA
#include
#include "kernel.h"
@@ -30,15 +30,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);
- }
}
@@ -47,7 +47,7 @@ void checkCUDAError(const char *msg, int line = -1) {
*****************/
/*! Block size used for CUDA kernel launch. */
-#define blockSize 128
+#define blockSize 128 //Originally 128
// LOOK-1.2 Parameters for the boids algorithm.
// These worked well in our reference implementation.
@@ -76,25 +76,28 @@ 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.
+glm::vec3* dev_pos_coherent; // positions reordered to match sorted-by-cell order
+glm::vec3* dev_vel1_coherent; // "current" velocities in coherent order (input to neighbor search)
+glm::vec3* dev_vel2_coherent; // "next" velocities in coherent order (output of neighbor search)
// LOOK-2.1 - Grid parameters based on simulation parameters.
// These are automatically computed for you in Boids::initSimulation
@@ -109,13 +112,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;
}
/**
@@ -123,63 +126,97 @@ __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.
+
+
+ ///////////////////////////////////////////////////////////////////////////
+ //TODO-2.1 - Allocate memory for uniform grid data structures
+
+ // Allocate memory for the particle array indices and grid indices
+ cudaMalloc((void**)&dev_particleArrayIndices, N * sizeof(int));
+ checkCUDAErrorWithLine("cudaMalloc dev_particleArrayIndices failed!");
+
+ // Allocate memory for the particle grid indices
+ cudaMalloc((void**)&dev_particleGridIndices, N * sizeof(int));
+ checkCUDAErrorWithLine("cudaMalloc dev_particleGridIndices failed!");
+
+ // Set up the thrust pointers
+ cudaMalloc((void**)&dev_gridCellStartIndices, gridCellCount * sizeof(int));
+ checkCUDAErrorWithLine("cudaMalloc dev_gridCellStartIndices failed!");
+
+ // Allocate memory for the grid cell end indices
+ cudaMalloc((void**)&dev_gridCellEndIndices, gridCellCount * sizeof(int));
+ checkCUDAErrorWithLine("cudaMalloc dev_gridCellEndIndices failed!");
+
+ ///////////////////////////////////////////////////////////////////////////
+ // TODO-2.3 - Allocate additional buffers here.
+ // Allocate memory for the coherent position and velocity buffers
+ 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();
}
@@ -190,42 +227,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();
}
@@ -239,47 +276,123 @@ 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
+ glm::vec3 perceived_center(0.0f, 0.0f, 0.0f);
+ int rule1_neighbors = 0; //number of neighbors within a certain distance for rule 1
+
+ // Rule 2: boids try to stay a distance d away from each other
+ glm::vec3 c(0.0f, 0.0f, 0.0f); // the center of mass of the local neighborhood
+
+ // Rule 3: boids try to match the speed of surrounding boids
+ glm::vec3 perceived_velocity(0.0f, 0.0f, 0.0f);
+ int rule3_neighbors = 0; //number of neighbors within a certain distance for rule 3
+
+ // Implementation of the 3 boids rules
+
+ // Loop over all boids
+ for (int i = 0; i < N; i++) {
+ if (i != iSelf) {
+ float distance = glm::distance(pos[i], pos[iSelf]); //distance between the two boids
+
+ // Rule1: Cohesion - Move towards perceived center
+ if (distance < rule1Distance) {
+ perceived_center += pos[i];
+ rule1_neighbors++;
+ }
+
+ // Rule2: Separation - Avoid crowding neighbors
+ if (distance < rule2Distance) {
+ c -= (pos[i] - pos[iSelf]);
+ }
+
+ // Rule 3: Alignment - Match velocity with nearby boids
+ if (distance < rule3Distance) {
+ perceived_velocity += vel[i];
+ rule3_neighbors++;
+ }
+ }
+ }
+
+ glm::vec3 result(0.0f, 0.0f, 0.0f); // The change in velocity
+
+ // Apply Rule 1
+ if (rule1_neighbors > 0) {
+ perceived_center /= rule1_neighbors; // Average position of neighbors => perceived_center /= number_of_neighbors
+ result += (perceived_center - pos[iSelf]) * rule1Scale; // Move towards the perceived center => (perceived_center - boid.position) * rule1Scale
+ }
+
+ // Apply Rule 2
+ result += c * rule2Scale; // Move away from neighbors => c * rule2Scale
+
+ // Apply Rule 3
+ if (rule3_neighbors > 0) {
+ perceived_velocity /= rule3_neighbors; // Average velocity of neighbors => perceived_velocity /= number_of_neighbors
+ result += perceived_velocity * rule3Scale; // Match the perceived velocity => perceived_velocity * rule3Scale
+ }
+
+ return result; // Return the total change in velocity
}
/**
* 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?
+ /*
+ We’re doing a read‑old and write‑new update
+ Every thread must see the same old velocities (vel1) while computing
+ if we write back into vel1 we will read partially‑updated data from other threads in the same step causing race behavior
+ */
+
+ int index = threadIdx.x + (blockIdx.x * blockDim.x); // Current boid index
+
+ // Ensure we don't go out of bounds
+ if (index >= N) {
+ return;
+ }
+
+ // Compute velocity change based on three rules
+ glm::vec3 velocityChange = computeVelocityChange(N, index, pos, vel1);
+
+ // Update velocity and clamp to max speed
+ glm::vec3 newVel = vel1[index] + velocityChange; // New velocity after applying rules
+ float speed = glm::length(newVel); // Calculate the speed (magnitude of velocity)
+
+ if (speed > maxSpeed) {
+ newVel = (newVel / speed) * maxSpeed; // Clamp to max speed
+ }
+
+ vel2[index] = newVel; // Store new velocity in vel2 for ping-ponging
}
/**
* 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.
@@ -289,179 +402,790 @@ __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
+
+ // Compute grid indices for each particle
+ int index = (blockIdx.x * blockDim.x) + threadIdx.x;
+
+ // Ensure we don't go out of bounds
+ if (index >= N) {
+ return;
+ }
+
+ // Compute 3D grid index for this particle
+ glm::vec3 gridPos = (pos[index] - gridMin) * inverseCellWidth; // Convert position to grid coordinates
+ int gridX = (int)floor(gridPos.x);
+ int gridY = (int)floor(gridPos.y);
+ int gridZ = (int)floor(gridPos.z);
+
+ // Clamp to grid bounds
+ gridX = imax(0, imin(gridX, gridResolution - 1));
+ gridY = imax(0, imin(gridY, gridResolution - 1));
+ gridZ = imax(0, imin(gridZ, gridResolution - 1));
+
+ // Store particle array index and grid index
+ indices[index] = index;
+ gridIndices[index] = gridIndex3Dto1D(gridX, gridY, gridZ, gridResolution);
}
// 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!"
+
+ int index = (blockIdx.x * blockDim.x) + threadIdx.x; // Current particle index
+
+ // Ensure we don't go out of bounds
+ if (index >= N) {
+ return;
+ }
+
+ int currentGridIndex = particleGridIndices[index]; // Grid index of the current particle
+
+ // If this is the first particle, it marks the start of its grid cell
+ if (index == 0) {
+ gridCellStartIndices[currentGridIndex] = index; // Start of current cell
+
+ }
+ else {
+ int previousGridIndex = particleGridIndices[index - 1]; // Grid index of the previous particle
+
+ // If the grid index changes, mark the end of the previous cell and start of the new cell
+ if (currentGridIndex != previousGridIndex) {
+ gridCellEndIndices[previousGridIndex] = index; // End of previous cell
+ gridCellStartIndices[currentGridIndex] = index; // Start of current cell
+ }
+ }
+
+ // If this is the last particle, it marks the end of its grid cell
+ if (index == N - 1) {
+ gridCellEndIndices[currentGridIndex] = index + 1; // End is exclusive
+ }
}
-__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!"
+////////////////////////////////////////////////////////////////////////
+// Helper functions
+
+// Reorder position/velocity into coherent arrays based on the sorted index mapping.
+// sortedIdx => originalIdx mapping lives in particleArrayIndices[sortedIdx].
+__global__ void kernReorderDataToCoherent(int N,
+ const int* __restrict__ particleArrayIndices,
+ const glm::vec3* __restrict__ pos_in,
+ const glm::vec3* __restrict__ vel1_in,
+ glm::vec3* __restrict__ pos_out,
+ glm::vec3* __restrict__ vel1_out) {
+
+ int sortedIdx = blockIdx.x * blockDim.x + threadIdx.x;
+
+ if (sortedIdx >= N) return;
+
+ const int originalIdx = particleArrayIndices[sortedIdx];
+ pos_out[sortedIdx] = pos_in[originalIdx];
+ vel1_out[sortedIdx] = vel1_in[originalIdx];
}
+// Scatter newly computed velocities from coherent order back to original order.
+// particleArrayIndices[sortedIdx] tells us where that coherent item came from.
+__global__ void kernScatterCoherentVelToUnsorted(int N,
+ const int* __restrict__ particleArrayIndices,
+ const glm::vec3* __restrict__ vel2_coherent,
+ glm::vec3* __restrict__ vel2_unsorted) {
+
+ int sortedIdx = blockIdx.x * blockDim.x + threadIdx.x;
+
+ if (sortedIdx >= N) return;
+
+ const int originalIdx = particleArrayIndices[sortedIdx];
+ vel2_unsorted[originalIdx] = vel2_coherent[sortedIdx];
+}
+
+////////////////////////////////////////////////////////////////////////
+
__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
+
+ // Grid-based neighbor search with scattered memory access
+ int index = (blockIdx.x * blockDim.x) + threadIdx.x;
+
+ // Ensure we don't go out of bounds
+ if (index >= N) {
+ return;
+ }
+
+ glm::vec3 thisPos = pos[index]; // Current particle position
+ glm::vec3 thisVel = vel1[index]; // Current particle velocity
+
+ // Get grid position
+ glm::vec3 gridPos = (thisPos - gridMin) * inverseCellWidth;
+ int gridX = (int)floor(gridPos.x);
+ int gridY = (int)floor(gridPos.y);
+ int gridZ = (int)floor(gridPos.z);
+
+ // Rule accumulators
+ glm::vec3 perceived_center(0.0f); // Rule 1
+ int rule1_neighbors = 0; //number of neighbors within a certain distance for rule 1
+
+ glm::vec3 c(0.0f); // Rule 2
+
+ glm::vec3 perceived_velocity(0.0f); // Rule 3
+ int rule3_neighbors = 0; //number of neighbors within a certain distance for rule 3
+
+ // Check neighboring cells (8 cells for cellWidth = 2 * searchRadius)
+
+ // Iterate over neighboring cells in a 3x3x3 cube around the current cell
+ for (int offsetX = -1; offsetX <= 1; offsetX++) {
+ for (int offsetY = -1; offsetY <= 1; offsetY++) {
+ for (int offsetZ = -1; offsetZ <= 1; offsetZ++) {
+
+ // Compute neighbor cell coordinates
+ int neighborX = gridX + offsetX;
+ int neighborY = gridY + offsetY;
+ int neighborZ = gridZ + offsetZ;
+
+ // Skip if out of bounds
+ if (neighborX < 0 || neighborX >= gridResolution ||
+ neighborY < 0 || neighborY >= gridResolution ||
+ neighborZ < 0 || neighborZ >= gridResolution) {
+ continue;
+ }
+
+ // Get start/end indices of boids in this neighbor cell
+ int neighborGridIndex = gridIndex3Dto1D(neighborX, neighborY, neighborZ, gridResolution);
+ int startIndex = gridCellStartIndices[neighborGridIndex];
+ int endIndex = gridCellEndIndices[neighborGridIndex];
+
+ // Skip empty cells
+ if (startIndex == -1) {
+ continue;
+ }
+
+ // Check all particles in this grid cell
+ for (int i = startIndex; i <= endIndex; i++) {
+ int neighborIndex = particleArrayIndices[i];
+
+ // Skip self
+ if (neighborIndex != index) {
+ glm::vec3 neighborPos = pos[neighborIndex];
+ float distance = glm::distance(neighborPos, thisPos);
+
+ // Rule 1: Cohesion - Move towards perceived center
+ if (distance < rule1Distance) {
+ perceived_center += neighborPos;
+ rule1_neighbors++;
+ }
+
+ // Rule 2: Separation - Avoid crowding neighbors
+ if (distance < rule2Distance) {
+ c -= (neighborPos - thisPos);
+ }
+
+ // Rule 3: Alignment - Match velocity with nearby boids
+ if (distance < rule3Distance) {
+ perceived_velocity += vel1[neighborIndex];
+ rule3_neighbors++;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // Compute final velocity change
+ glm::vec3 result(0.0f);
+
+ // Apply Rule 1
+ if (rule1_neighbors > 0) {
+ perceived_center /= rule1_neighbors;
+ result += (perceived_center - thisPos) * rule1Scale;
+ }
+
+ // Apply Rule 2
+ result += c * rule2Scale;
+
+ // Apply Rule 3
+ if (rule3_neighbors > 0) {
+ perceived_velocity /= rule3_neighbors;
+ result += perceived_velocity * rule3Scale;
+ }
+
+ // Update velocity and clamp
+ glm::vec3 newVel = thisVel + result; // New velocity after applying rules
+ float speed = glm::length(newVel); // Calculate the speed (magnitude of velocity)
+ if (speed > maxSpeed) {
+ newVel = (newVel / speed) * maxSpeed; // Clamp to max speed
+ }
+
+ // Write new velocity to vel2
+ vel2[index] = 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
+
+
+ const int index = blockIdx.x * blockDim.x + threadIdx.x; // Current boid index
+
+ if (index >= N) return; // Ensure we don't go out of bounds
+
+ // These arrays are already in cell-sorted (coherent) order.
+ const glm::vec3 currentPos = pos[index]; // Current boid position
+ const glm::vec3 currentVel = vel1[index]; // Current boid velocity
+
+ // Figure out which grid cell this boid lives in
+ const glm::vec3 gridPos = (currentPos - gridMin) * inverseCellWidth; // Convert position to grid coordinates
+ const int gridX = (int)floor(gridPos.x);
+ const int gridY = (int)floor(gridPos.y);
+ const int gridZ = (int)floor(gridPos.z);
+
+ // Accumulators for Boids rules
+
+ glm::vec3 perceived_center(0.0f); int rule1_neighbors = 0; // Rule 1 (cohesion)
+ glm::vec3 c(0.0f); // Rule 2 (separation)
+ glm::vec3 perceived_velocity(0.0f); int rule3_neighbors = 0; // Rule 3 (alignment)
+
+ // Visit neighboring cells in a 3x3x3 block
+ for (int offsetZ = -1; offsetZ <= 1; ++offsetZ) {
+ for (int offsetY = -1; offsetY <= 1; ++offsetY) {
+ for (int offsetX = -1; offsetX <= 1; ++offsetX) {
+
+ // Compute neighbor cell coordinates
+ int neighborX = gridX + offsetX;
+ int neighborY = gridY + offsetY;
+ int neighborZ = gridZ + offsetZ;
+
+ // Skip if out of bounds
+ if (neighborX < 0 || neighborX >= gridResolution ||
+ neighborY < 0 || neighborY >= gridResolution ||
+ neighborZ < 0 || neighborZ >= gridResolution) {
+ continue;
+ }
+
+ // Compute the 1D index of the neighbor cell
+ const int neighborCell1D = gridIndex3Dto1D(neighborX, neighborY, neighborZ, gridResolution);
+
+ // Look up the contiguous run [start, end) of boids in this cell
+ const int startIndex = gridCellStartIndices[neighborCell1D]; // INCLUSIVE
+ const int endIndex = gridCellEndIndices[neighborCell1D]; // EXCLUSIVE
+
+ if (startIndex == -1) continue; // empty cell
+
+ // Sequentially scan neighbors in this cell (coalesced reads)
+ for (int j = startIndex; j < endIndex; ++j) {
+
+ if (j == index) continue; // skip self (same coherent index)
+
+ const glm::vec3 neighborPos = pos[j];
+ const float distanceBoid = glm::distance(neighborPos, currentPos);
+
+ // Rule 1: Cohesion – move toward center of mass of neighbors within rule1Distance
+ if (distanceBoid < rule1Distance) {
+ perceived_center += neighborPos;
+ ++rule1_neighbors;
+ }
+
+ // Rule 2: Separation – avoid crowding neighbors within rule2Distance
+ if (distanceBoid < rule2Distance) {
+ c -= (neighborPos - currentPos);
+ }
+
+ // Rule 3: Alignment – align velocity with neighbors within rule3Distance
+ if (distanceBoid < rule3Distance) {
+ perceived_velocity += vel1[j];
+ ++rule3_neighbors;
+ }
+ }
+ }
+ }
+ }
+
+ // Combine rule contributions
+ glm::vec3 delta_v(0.0f); // change in velocity (delta v)
+
+ //rule 1
+ if (rule1_neighbors > 0) {
+ perceived_center /= rule1_neighbors; // Average position of neighbors
+ delta_v += (perceived_center - currentPos) * rule1Scale; // Move towards the perceived center
+ }
+
+ //rule 2
+ delta_v += c * rule2Scale; // Move away from neighbors
+
+ //rule 3
+ if (rule3_neighbors > 0) {
+ perceived_velocity /= rule3_neighbors; // Average velocity of neighbors
+ delta_v += perceived_velocity * rule3Scale; // Match the perceived velocity
+ }
+
+ // Clamp speed
+ glm::vec3 outVel = currentVel + delta_v; // New velocity after applying rules
+ const float speed = glm::length(outVel); // Calculate the speed (magnitude of velocity)
+
+ if (speed > maxSpeed) {
+ outVel = (outVel / speed) * maxSpeed; // Clamp to max speeds
+ }
+
+ // Write new velocity in coherent order
+ vel2[index] = outVel;
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////////
+// Extra credit - grid loop optimization
+
+// For this kernel, the uniform grid is constructed such that the cell width is
+__global__ void kernUpdateVelNeighborSearchCoherentOptimized(
+ int N, int gridResolution, glm::vec3 gridMin,
+ float inverseCellWidth, float cellWidth,
+ int* gridCellStartIndices, int* gridCellEndIndices,
+ glm::vec3* pos, glm::vec3* vel1, glm::vec3* vel2)
+{
+
+ int index = (blockIdx.x * blockDim.x) + threadIdx.x; // Current boid index
+
+ if (index >= N) return; // Ensure we don't go out of bounds
+
+ // Current boid’s position and velocity (coherent arrays)
+ glm::vec3 currentPos = pos[index]; // Current boid position
+ glm::vec3 currentVel = vel1[index]; // Current boid velocity
+
+ // Determine this boid’s grid cell coordinates
+ glm::vec3 gridPos = (currentPos - gridMin) * inverseCellWidth; // Convert position to grid coordinates
+ const int gridX = (int)floor(gridPos.x);
+ const int gridY = (int)floor(gridPos.y);
+ const int gridZ = (int)floor(gridPos.z);
+
+ // Compute search radius (max neighbor influence distance)
+ float maxDistance = fmaxf(rule1Distance, fmaxf(rule2Distance, rule3Distance));
+
+ // Determine grid cell index range in each dimension that lies within `maxDistance`
+ int minCellX = static_cast(floorf((currentPos.x - gridMin.x - maxDistance) * inverseCellWidth));
+ int maxCellX = static_cast(floorf((currentPos.x - gridMin.x + maxDistance) * inverseCellWidth));
+ int minCellY = static_cast(floorf((currentPos.y - gridMin.y - maxDistance) * inverseCellWidth));
+ int maxCellY = static_cast(floorf((currentPos.y - gridMin.y + maxDistance) * inverseCellWidth));
+ int minCellZ = static_cast(floorf((currentPos.z - gridMin.z - maxDistance) * inverseCellWidth));
+ int maxCellZ = static_cast(floorf((currentPos.z - gridMin.z + maxDistance) * inverseCellWidth));
+
+ // Clamp the cell index ranges to the grid bounds
+ if (minCellX < 0) minCellX = 0;
+ if (minCellY < 0) minCellY = 0;
+ if (minCellZ < 0) minCellZ = 0;
+ if (maxCellX >= gridResolution) maxCellX = gridResolution - 1;
+ if (maxCellY >= gridResolution) maxCellY = gridResolution - 1;
+ if (maxCellZ >= gridResolution) maxCellZ = gridResolution - 1;
+
+ // Accumulators for the three Boids rules
+ glm::vec3 perceived_center(0.0f); int rule1_neighbors = 0; // Cohesion
+ glm::vec3 c(0.0f); // Separation
+ glm::vec3 perceived_velocity(0.0f); int rule3_neighbors = 0; // Alignment
+
+ // Loop over all candidate neighbor cells in the computed range
+ for (int z = minCellZ; z <= maxCellZ; ++z) {
+ for (int y = minCellY; y <= maxCellY; ++y) {
+ for (int x = minCellX; x <= maxCellX; ++x) {
+
+
+ int cellIndex = gridIndex3Dto1D(x, y, z, gridResolution); // Get the 1D index of this neighbor cell
+ int startIndex = gridCellStartIndices[cellIndex]; // INCLUSIVE
+ int endIndex = gridCellEndIndices[cellIndex];// EXCLUSIVE
+
+ if (startIndex == -1) continue; // skip empty cells
+
+ // Check all boids in this cell
+ for (int j = startIndex; j < endIndex; ++j) {
+ if (j == index) continue; // skip itself
+
+ // Compute distance to this neighbor
+ glm::vec3 neighborPos = pos[j]; // Neighbor boid position
+ float dist = glm::distance(neighborPos, currentPos); // Distance to neighbor
+
+
+ if (dist < rule1Distance) { //rules1 - Cohesion
+ perceived_center += neighborPos;
+ rule1_neighbors++;
+ }
+
+ if (dist < rule2Distance) { // rule2 - Separation
+ c -= (neighborPos - currentPos);
+ }
+
+ if (dist < rule3Distance) { // rule3 -Alignment
+ perceived_velocity += vel1[j];
+ rule3_neighbors++;
+ }
+ }
+ }
+ }
+ }
+
+ // Apply the three rules to compute velocity change
+ glm::vec3 deltaV(0.0f); // Change in velocity (delta v)
+
+ // Rule 1: Cohesion
+ if (rule1_neighbors > 0) {
+ perceived_center /= rule1_neighbors;
+ deltaV += (perceived_center - currentPos) * rule1Scale;
+ }
+
+ // Rule 2: Separation
+ deltaV += c * rule2Scale;
+
+ // Rule 3: Alignment
+ if (rule3_neighbors > 0) {
+ perceived_velocity /= rule3_neighbors;
+ deltaV += perceived_velocity * rule3Scale;
+ }
+
+ // Clamp the new speed to maxSpeed
+ glm::vec3 newVel = currentVel + deltaV;
+ float speed = glm::length(newVel);
+ if (speed > maxSpeed) {
+ newVel = (newVel / speed) * maxSpeed;
+ }
+
+ // Write back the new velocity (coherent order)
+ vel2[index] = newVel;
}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
/**
* Step the entire N-body simulation by `dt` seconds.
*/
void Boids::stepSimulationNaive(float dt) {
- // TODO-1.2 - use the kernels you wrote to step the simulation forward in time.
- // TODO-1.2 ping-pong the velocity buffers
+ // TODO-1.2 - use the kernels you wrote to step the simulation forward in time.
+ // TODO-1.2 ping-pong the velocity buffers
+ //Naive method
+
+ dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); // Number of blocks needed
+
+ // Update velocity (naive O(N^2) neighbor search)
+ kernUpdateVelocityBruteForce << > > (numObjects, dev_pos, dev_vel1, dev_vel2);
+ checkCUDAErrorWithLine("kernUpdateVelocityBruteForce failed!");
+
+ // Integrate positions using the NEW velocities in dev_vel2s
+ kernUpdatePos << > > (numObjects, dt, dev_pos, dev_vel2);
+ checkCUDAErrorWithLine("kernUpdatePos failed!");
+
+ // Swap 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.
+ // In Parallel:
+ // - label each particle with its array index as well as its grid index.
+ // Use 2x width grids.
+ // - Unstable key sort using Thrust. A stable sort isn't necessary, but you
+ // are welcome to do a performance comparison.
+ // - Naively unroll the loop for finding the start and end indices of each
+ // cell's data pointers in the array of boid indices
+ // - Perform velocity updates using neighbor search
+ // - Update positions
+ // - Ping-pong buffers as needed
+
+ dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); // Number of blocks needed
+ dim3 gridBlocksPerGrid((gridCellCount + blockSize - 1) / blockSize); // Number of blocks for grid cells
+
+ // Reset grid indices
+ kernResetIntBuffer << > > (gridCellCount, dev_gridCellStartIndices, -1); // Reset start indices to -1
+ kernResetIntBuffer << > > (gridCellCount, dev_gridCellEndIndices, -1); // Reset end indices to -1
+ checkCUDAErrorWithLine("kernResetIntBuffer failed!");
+
+ // Compute grid indices for each particle
+ kernComputeIndices << > > (numObjects, gridSideCount,
+ gridMinimum, gridInverseCellWidth, dev_pos, dev_particleArrayIndices, dev_particleGridIndices); // Label each particle with its array index and grid index
+ checkCUDAErrorWithLine("kernComputeIndices failed!");
+
+ // Sort particles by grid index
+ dev_thrust_particleArrayIndices = thrust::device_pointer_cast(dev_particleArrayIndices); // Set up the thrust pointers
+ dev_thrust_particleGridIndices = thrust::device_pointer_cast(dev_particleGridIndices); // Set up the thrust pointers
+ thrust::sort_by_key(dev_thrust_particleGridIndices, dev_thrust_particleGridIndices + numObjects,
+ dev_thrust_particleArrayIndices); // Sort based on grid indices
+
+ // Identify cell start and end
+ kernIdentifyCellStartEnd << > > (numObjects,
+ dev_particleGridIndices, dev_gridCellStartIndices, dev_gridCellEndIndices); // Identify start/end of each cell
+ checkCUDAErrorWithLine("kernIdentifyCellStartEnd failed!");
+
+ // Update velocities using grid
+ kernUpdateVelNeighborSearchScattered << > > (
+ numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, gridCellWidth,
+ dev_gridCellStartIndices, dev_gridCellEndIndices, dev_particleArrayIndices,
+ dev_pos, dev_vel1, dev_vel2); // Update velocities based on neighbor search
+ checkCUDAErrorWithLine("kernUpdateVelNeighborSearchScattered failed!");
+
+ // Update positions
+ kernUpdatePos << > > (numObjects, dt, dev_pos, dev_vel2); // Update positions
+ checkCUDAErrorWithLine("kernUpdatePos failed!");
+
+ // Swap velocity buffers
+ std::swap(dev_vel1, dev_vel2);
+
}
void Boids::stepSimulationCoherentGrid(float dt) {
- // TODO-2.3 - start by copying Boids::stepSimulationNaiveGrid
- // Uniform Grid Neighbor search using Thrust sort on cell-coherent data.
- // In Parallel:
- // - Label each particle with its array index as well as its grid index.
- // Use 2x width grids
- // - Unstable key sort using Thrust. A stable sort isn't necessary, but you
- // are welcome to do a performance comparison.
- // - Naively unroll the loop for finding the start and end indices of each
- // cell's data pointers in the array of boid indices
- // - BIG DIFFERENCE: use the rearranged array index buffer to reshuffle all
- // the particle data in the simulation array.
- // CONSIDER WHAT ADDITIONAL BUFFERS YOU NEED
- // - Perform velocity updates using neighbor search
- // - Update positions
- // - Ping-pong buffers as needed. THIS MAY BE DIFFERENT FROM BEFORE.
+ // TODO-2.3 - start by copying Boids::stepSimulationNaiveGrid
+ // Uniform Grid Neighbor search using Thrust sort on cell-coherent data.
+ // In Parallel:
+ // - Label each particle with its array index as well as its grid index.
+ // Use 2x width grids
+ // - Unstable key sort using Thrust. A stable sort isn't necessary, but you
+ // are welcome to do a performance comparison.
+ // - Naively unroll the loop for finding the start and end indices of each
+ // cell's data pointers in the array of boid indices
+ // - BIG DIFFERENCE: use the rearranged array index buffer to reshuffle all
+ // the particle data in the simulation array.
+ // CONSIDER WHAT ADDITIONAL BUFFERS YOU NEED
+ // - Perform velocity updates using neighbor search
+ // - Update positions
+ // - Ping-pong buffers as needed. THIS MAY BE DIFFERENT FROM BEFORE.
+
+ // Blocks for boids vs grid
+ dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize);
+ dim3 gridBlocksPerGrid((gridCellCount + blockSize - 1) / blockSize);
+
+ // 1) Clear cell start/end to empty (-1)
+ kernResetIntBuffer << > > (gridCellCount, dev_gridCellStartIndices, -1);
+ kernResetIntBuffer << > > (gridCellCount, dev_gridCellEndIndices, -1);
+ checkCUDAErrorWithLine("kernResetIntBuffer failed!");
+
+ // 2) Label each boid with (a) its original index and (b) its grid cell id
+ kernComputeIndices << > > (
+ numObjects, gridSideCount, gridMinimum, gridInverseCellWidth,
+ dev_pos, dev_particleArrayIndices, dev_particleGridIndices);
+ checkCUDAErrorWithLine("kernComputeIndices failed!");
+
+ // 3) Sort by grid cell (keys: cell ids, values: boid indices)
+ dev_thrust_particleArrayIndices = thrust::device_pointer_cast(dev_particleArrayIndices);
+ dev_thrust_particleGridIndices = thrust::device_pointer_cast(dev_particleGridIndices);
+ thrust::sort_by_key(dev_thrust_particleGridIndices,
+ dev_thrust_particleGridIndices + numObjects,
+ dev_thrust_particleArrayIndices);
+
+ // 4) Build cell start/end tables over the sorted cell-id array
+ kernIdentifyCellStartEnd << > > (
+ numObjects, dev_particleGridIndices,
+ dev_gridCellStartIndices, dev_gridCellEndIndices);
+ checkCUDAErrorWithLine("kernIdentifyCellStartEnd failed!");
+
+ // 5) Reorder position/velocity into coherent arrays (contiguous per cell)
+ kernReorderDataToCoherent << > > (
+ numObjects, dev_particleArrayIndices,
+ dev_pos, dev_vel1,
+ dev_pos_coherent, dev_vel1_coherent);
+ checkCUDAErrorWithLine("kernReorderDataToCoherent failed!");
+
+ // 6) Neighbor search directly over coherent arrays, write vel2_coherent
+ kernUpdateVelNeighborSearchCoherent << > > (
+ numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, gridCellWidth,
+ dev_gridCellStartIndices, dev_gridCellEndIndices,
+ dev_pos_coherent, dev_vel1_coherent, dev_vel2_coherent);
+ checkCUDAErrorWithLine("kernUpdateVelNeighborSearchCoherent failed!");
+
+ // 7) Scatter new velocities back to original order for integration
+ kernScatterCoherentVelToUnsorted << > > (
+ numObjects, dev_particleArrayIndices, dev_vel2_coherent, dev_vel2);
+ checkCUDAErrorWithLine("kernScatterCoherentVelToUnsorted failed!");
+
+ // 8) Integrate positions with the new velocities (unsorted arrays)
+ kernUpdatePos << > > (numObjects, dt, dev_pos, dev_vel2);
+ checkCUDAErrorWithLine("kernUpdatePos failed!");
+
+ // 9) Ping-pong velocity buffers (dev_vel1 holds newest after swap)
+ std::swap(dev_vel1, dev_vel2);
}
-void Boids::endSimulation() {
- cudaFree(dev_vel1);
- cudaFree(dev_vel2);
- cudaFree(dev_pos);
+///////////////////////////////////////////////////////////////////////
+// Extra Credit - grid loop optimization
+
+
+void Boids::stepSimulationGridLoopOptimized(float dt) {
+
+ // Launch configurations
+ dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize);
+ dim3 gridBlocksPerGrid((gridCellCount + blockSize - 1) / blockSize);
+
+ // 1) Reset grid cell start and end indices to -1 (empty)
+ kernResetIntBuffer << > > (gridCellCount, dev_gridCellStartIndices, -1);
+ kernResetIntBuffer << > > (gridCellCount, dev_gridCellEndIndices, -1);
+ checkCUDAErrorWithLine("kernResetIntBuffer failed!");
+
+ // 2) Compute grid indices for each boid (map each particle to its grid cell)
+ kernComputeIndices << > > (
+ numObjects, gridSideCount, gridMinimum, gridInverseCellWidth,
+ dev_pos, dev_particleArrayIndices, dev_particleGridIndices);
+ checkCUDAErrorWithLine("kernComputeIndices failed!");
+
+ // 3) Sort boids by grid index using Thrust (cell indices as keys, boid indices as values)
+ dev_thrust_particleArrayIndices = thrust::device_pointer_cast(dev_particleArrayIndices);
+ dev_thrust_particleGridIndices = thrust::device_pointer_cast(dev_particleGridIndices);
+ thrust::sort_by_key(dev_thrust_particleGridIndices,
+ dev_thrust_particleGridIndices + numObjects,
+ dev_thrust_particleArrayIndices);
+
+ // 4) Identify the start and end index of each cell’s boid list
+ kernIdentifyCellStartEnd << > > (
+ numObjects, dev_particleGridIndices,
+ dev_gridCellStartIndices, dev_gridCellEndIndices);
+ checkCUDAErrorWithLine("kernIdentifyCellStartEnd failed!");
+
+ // 5) Reorder boid data into contiguous memory by cell (coherent arrays)
+ kernReorderDataToCoherent << > > (
+ numObjects, dev_particleArrayIndices,
+ dev_pos, dev_vel1,
+ dev_pos_coherent, dev_vel1_coherent);
+ checkCUDAErrorWithLine("kernReorderDataToCoherent failed!");
+
+ // 6) Update velocities using dynamic grid-loop neighbor search (optimized kernel)
+ kernUpdateVelNeighborSearchCoherentOptimized << > > (
+ numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, gridCellWidth,
+ dev_gridCellStartIndices, dev_gridCellEndIndices,
+ dev_pos_coherent, dev_vel1_coherent, dev_vel2_coherent);
+ checkCUDAErrorWithLine("kernUpdateVelNeighborSearchCoherentOptimized failed!");
+
+ // 7) Scatter the new velocities from coherent order back to original order
+ kernScatterCoherentVelToUnsorted << > > (
+ numObjects, dev_particleArrayIndices, dev_vel2_coherent, dev_vel2);
+ checkCUDAErrorWithLine("kernScatterCoherentVelToUnsorted failed!");
+
+ // 8) Integrate positions using the updated velocities (dev_vel2)
+ kernUpdatePos << > > (numObjects, dt, dev_pos, dev_vel2);
+ checkCUDAErrorWithLine("kernUpdatePos failed!");
+
+ // 9) Ping-pong the velocity buffers for next iteration (make dev_vel1 current)
+ std::swap(dev_vel1, dev_vel2);
+}
+
+
+////////////////////////////////////////////////////////////////////////////
- // TODO-2.1 TODO-2.3 - Free any additional buffers here.
+void Boids::endSimulation() {
+ cudaFree(dev_vel1);
+ cudaFree(dev_vel2);
+ cudaFree(dev_pos);
+
+ // TODO-2.1 TODO-2.3 - Free any additional buffers here.
+
+ //TODO-2.1 - Free the memory allocated for the uniform grid data structures
+ cudaFree(dev_particleArrayIndices); // Free particle array indices
+ cudaFree(dev_particleGridIndices); // Free particle grid indices
+ cudaFree(dev_gridCellStartIndices); // Free grid cell start indices
+ cudaFree(dev_gridCellEndIndices); // Free grid cell end indices
+
+ //TODO-2.3 - Free the additional buffers for coherent grid
+ cudaFree(dev_pos_coherent); // Free coherent position buffer
+ cudaFree(dev_vel1_coherent); // Free coherent velocity buffer
+ cudaFree(dev_vel2_coherent); // Free coherent velocity buffer
}
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;
+}
\ No newline at end of file
diff --git a/src/kernel.h b/src/kernel.h
index a38b64d..cf9a078 100644
--- a/src/kernel.h
+++ b/src/kernel.h
@@ -5,6 +5,8 @@ namespace Boids {
void stepSimulationNaive(float dt);
void stepSimulationScatteredGrid(float dt);
void stepSimulationCoherentGrid(float dt);
+ void stepSimulationGridLoopOptimized(float dt); // Extra credit
+
void copyBoidsToVBO(float *vbodptr_positions, float *vbodptr_velocities);
void endSimulation();
diff --git a/src/main.cpp b/src/main.cpp
index 9c917c0..5c954aa 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -23,26 +23,30 @@
// 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
+
+#define GRID_LOOP_OPTIMIZATION 0 // Enable Grid-Looping Optimization mode - extra credit
+
// LOOK-1.2 - change this to adjust particle count in the simulation
-const int N_FOR_VIS = 5000;
+const int N_FOR_VIS = 75000; //original 5000
const float DT = 0.2f;
/**
* C main function.
*/
int main(int argc, char* argv[]) {
- projectName = "5650 CUDA Intro: Boids";
-
- if (init(argc, argv)) {
- mainLoop();
- Boids::endSimulation();
- return 0;
- } else {
- return 1;
- }
+ projectName = "5650 CUDA Intro: Boids";
+
+ if (init(argc, argv)) {
+ mainLoop();
+ Boids::endSimulation();
+ return 0;
+ }
+ else {
+ return 1;
+ }
}
//-------------------------------
@@ -50,260 +54,274 @@ int main(int argc, char* argv[]) {
//-------------------------------
std::string deviceName;
-GLFWwindow *window;
+GLFWwindow* window;
/**
* Initialization of CUDA and GLFW.
*/
-bool init(int argc, char **argv) {
- // Set window title to "Student Name: [SM 2.0] GPU Name"
- cudaDeviceProp deviceProp;
- int gpuDevice = 0;
- int device_count = 0;
- cudaGetDeviceCount(&device_count);
- if (gpuDevice > device_count) {
- std::cout
- << "Error: GPU device number is greater than the number of devices!"
- << " Perhaps a CUDA-capable GPU is not installed?"
- << std::endl;
- return false;
- }
- cudaGetDeviceProperties(&deviceProp, gpuDevice);
- int major = deviceProp.major;
- int minor = deviceProp.minor;
-
- std::ostringstream ss;
- ss << projectName << " [SM " << major << "." << minor << " " << deviceProp.name << "]";
- deviceName = ss.str();
-
- // Window setup stuff
- glfwSetErrorCallback(errorCallback);
-
- if (!glfwInit()) {
- std::cout
- << "Error: Could not initialize GLFW!"
- << " Perhaps OpenGL 3.3 isn't available?"
- << std::endl;
- return false;
- }
-
- glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3);
- glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 3);
- glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE);
- glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
-
- window = glfwCreateWindow(width, height, deviceName.c_str(), NULL, NULL);
- if (!window) {
- glfwTerminate();
- return false;
- }
- glfwMakeContextCurrent(window);
- glfwSetKeyCallback(window, keyCallback);
- glfwSetCursorPosCallback(window, mousePositionCallback);
- glfwSetMouseButtonCallback(window, mouseButtonCallback);
+bool init(int argc, char** argv) {
+ // Set window title to "Student Name: [SM 2.0] GPU Name"
+ cudaDeviceProp deviceProp;
+ int gpuDevice = 0;
+ int device_count = 0;
+ cudaGetDeviceCount(&device_count);
+ if (gpuDevice > device_count) {
+ std::cout
+ << "Error: GPU device number is greater than the number of devices!"
+ << " Perhaps a CUDA-capable GPU is not installed?"
+ << std::endl;
+ return false;
+ }
+ cudaGetDeviceProperties(&deviceProp, gpuDevice);
+ int major = deviceProp.major;
+ int minor = deviceProp.minor;
+
+ std::ostringstream ss;
+ ss << projectName << " [SM " << major << "." << minor << " " << deviceProp.name << "]";
+ deviceName = ss.str();
+
+ // Window setup stuff
+ glfwSetErrorCallback(errorCallback);
+
+ if (!glfwInit()) {
+ std::cout
+ << "Error: Could not initialize GLFW!"
+ << " Perhaps OpenGL 3.3 isn't available?"
+ << std::endl;
+ return false;
+ }
+
+ glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3);
+ glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 3);
+ glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE);
+ glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
- glewExperimental = GL_TRUE;
- if (glewInit() != GLEW_OK) {
- return false;
- }
+ window = glfwCreateWindow(width, height, deviceName.c_str(), NULL, NULL);
+ if (!window) {
+ glfwTerminate();
+ return false;
+ }
+ glfwMakeContextCurrent(window);
+ glfwSetKeyCallback(window, keyCallback);
+ glfwSetCursorPosCallback(window, mousePositionCallback);
+ glfwSetMouseButtonCallback(window, mouseButtonCallback);
+
+ glewExperimental = GL_TRUE;
+ if (glewInit() != GLEW_OK) {
+ return false;
+ }
- // Initialize drawing state
- initVAO();
+ // Initialize drawing state
+ initVAO();
- // Default to device ID 0. If you have more than one GPU and want to test a non-default one,
- // change the device ID.
- cudaGLSetGLDevice(0);
+ // Default to device ID 0. If you have more than one GPU and want to test a non-default one,
+ // change the device ID.
+ cudaGLSetGLDevice(0);
- cudaGLRegisterBufferObject(boidVBO_positions);
- cudaGLRegisterBufferObject(boidVBO_velocities);
+ cudaGLRegisterBufferObject(boidVBO_positions);
+ cudaGLRegisterBufferObject(boidVBO_velocities);
- // Initialize N-body simulation
- Boids::initSimulation(N_FOR_VIS);
+ // Initialize N-body simulation
+ Boids::initSimulation(N_FOR_VIS);
- updateCamera();
+ updateCamera();
- initShaders(program);
+ initShaders(program);
- glEnable(GL_DEPTH_TEST);
+ glEnable(GL_DEPTH_TEST);
- return true;
+ return true;
}
void initVAO() {
- std::unique_ptr bodies{ new GLfloat[4 * (N_FOR_VIS)] };
- std::unique_ptr bindices{ new GLuint[N_FOR_VIS] };
+ std::unique_ptr bodies{ new GLfloat[4 * (N_FOR_VIS)] };
+ std::unique_ptr bindices{ new GLuint[N_FOR_VIS] };
- glm::vec4 ul(-1.0, -1.0, 1.0, 1.0);
- glm::vec4 lr(1.0, 1.0, 0.0, 0.0);
+ glm::vec4 ul(-1.0, -1.0, 1.0, 1.0);
+ glm::vec4 lr(1.0, 1.0, 0.0, 0.0);
- for (int i = 0; i < N_FOR_VIS; i++) {
- bodies[4 * i + 0] = 0.0f;
- bodies[4 * i + 1] = 0.0f;
- bodies[4 * i + 2] = 0.0f;
- bodies[4 * i + 3] = 1.0f;
- bindices[i] = i;
- }
+ for (int i = 0; i < N_FOR_VIS; i++) {
+ bodies[4 * i + 0] = 0.0f;
+ bodies[4 * i + 1] = 0.0f;
+ bodies[4 * i + 2] = 0.0f;
+ bodies[4 * i + 3] = 1.0f;
+ bindices[i] = i;
+ }
- glGenVertexArrays(1, &boidVAO); // Attach everything needed to draw a particle to this
- glGenBuffers(1, &boidVBO_positions);
- glGenBuffers(1, &boidVBO_velocities);
- glGenBuffers(1, &boidIBO);
+ glGenVertexArrays(1, &boidVAO); // Attach everything needed to draw a particle to this
+ glGenBuffers(1, &boidVBO_positions);
+ glGenBuffers(1, &boidVBO_velocities);
+ glGenBuffers(1, &boidIBO);
- glBindVertexArray(boidVAO);
+ glBindVertexArray(boidVAO);
- // Bind the positions array to the boidVAO by way of the boidVBO_positions
- glBindBuffer(GL_ARRAY_BUFFER, boidVBO_positions); // bind the buffer
- glBufferData(GL_ARRAY_BUFFER, 4 * (N_FOR_VIS) * sizeof(GLfloat), bodies.get(), GL_DYNAMIC_DRAW); // transfer data
+ // Bind the positions array to the boidVAO by way of the boidVBO_positions
+ glBindBuffer(GL_ARRAY_BUFFER, boidVBO_positions); // bind the buffer
+ glBufferData(GL_ARRAY_BUFFER, 4 * (N_FOR_VIS) * sizeof(GLfloat), bodies.get(), GL_DYNAMIC_DRAW); // transfer data
- glEnableVertexAttribArray(positionLocation);
- glVertexAttribPointer((GLuint)positionLocation, 4, GL_FLOAT, GL_FALSE, 0, 0);
+ glEnableVertexAttribArray(positionLocation);
+ glVertexAttribPointer((GLuint)positionLocation, 4, GL_FLOAT, GL_FALSE, 0, 0);
- // Bind the velocities array to the boidVAO by way of the boidVBO_velocities
- glBindBuffer(GL_ARRAY_BUFFER, boidVBO_velocities);
- glBufferData(GL_ARRAY_BUFFER, 4 * (N_FOR_VIS) * sizeof(GLfloat), bodies.get(), GL_DYNAMIC_DRAW);
- glEnableVertexAttribArray(velocitiesLocation);
- glVertexAttribPointer((GLuint)velocitiesLocation, 4, GL_FLOAT, GL_FALSE, 0, 0);
+ // Bind the velocities array to the boidVAO by way of the boidVBO_velocities
+ glBindBuffer(GL_ARRAY_BUFFER, boidVBO_velocities);
+ glBufferData(GL_ARRAY_BUFFER, 4 * (N_FOR_VIS) * sizeof(GLfloat), bodies.get(), GL_DYNAMIC_DRAW);
+ glEnableVertexAttribArray(velocitiesLocation);
+ glVertexAttribPointer((GLuint)velocitiesLocation, 4, GL_FLOAT, GL_FALSE, 0, 0);
- glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, boidIBO);
- glBufferData(GL_ELEMENT_ARRAY_BUFFER, (N_FOR_VIS) * sizeof(GLuint), bindices.get(), GL_STATIC_DRAW);
+ glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, boidIBO);
+ glBufferData(GL_ELEMENT_ARRAY_BUFFER, (N_FOR_VIS) * sizeof(GLuint), bindices.get(), GL_STATIC_DRAW);
- glBindVertexArray(0);
+ glBindVertexArray(0);
}
-void initShaders(GLuint * program) {
- GLint location;
+void initShaders(GLuint* program) {
+ GLint location;
- program[PROG_BOID] = glslUtility::createProgram(
- "shaders/boid.vert.glsl",
- "shaders/boid.geom.glsl",
- "shaders/boid.frag.glsl", attributeLocations, 2);
+ program[PROG_BOID] = glslUtility::createProgram(
+ "shaders/boid.vert.glsl",
+ "shaders/boid.geom.glsl",
+ "shaders/boid.frag.glsl", attributeLocations, 2);
glUseProgram(program[PROG_BOID]);
if ((location = glGetUniformLocation(program[PROG_BOID], "u_projMatrix")) != -1) {
- glUniformMatrix4fv(location, 1, GL_FALSE, &projection[0][0]);
+ glUniformMatrix4fv(location, 1, GL_FALSE, &projection[0][0]);
}
if ((location = glGetUniformLocation(program[PROG_BOID], "u_cameraPos")) != -1) {
- glUniform3fv(location, 1, &cameraPosition[0]);
+ glUniform3fv(location, 1, &cameraPosition[0]);
}
- }
+}
- //====================================
- // Main loop
- //====================================
- void runCUDA() {
+//====================================
+// Main loop
+//====================================
+void runCUDA() {
// Map OpenGL buffer object for writing from CUDA on a single GPU
// No data is moved (Win & Linux). When mapped to CUDA, OpenGL should not
// use this buffer
- float4 *dptr = NULL;
- float *dptrVertPositions = NULL;
- float *dptrVertVelocities = NULL;
+ float4* dptr = NULL;
+ float* dptrVertPositions = NULL;
+ float* dptrVertVelocities = NULL;
cudaGLMapBufferObject((void**)&dptrVertPositions, boidVBO_positions);
cudaGLMapBufferObject((void**)&dptrVertVelocities, boidVBO_velocities);
// execute the kernel
- #if UNIFORM_GRID && COHERENT_GRID
+/*
+* // Original code
+#if UNIFORM_GRID && COHERENT_GRID
+ Boids::stepSimulationCoherentGrid(DT);
+#elif UNIFORM_GRID
+ Boids::stepSimulationScatteredGrid(DT);
+#else
+ Boids::stepSimulationNaive(DT);
+#endif
+*/
+
+// Extra Credit: Grid-Looping Optimization
+#if UNIFORM_GRID && COHERENT_GRID && GRID_LOOP_OPTIMIZATION
+ Boids::stepSimulationGridLoopOptimized(DT);
+#elif UNIFORM_GRID && COHERENT_GRID
Boids::stepSimulationCoherentGrid(DT);
- #elif UNIFORM_GRID
+#elif UNIFORM_GRID
Boids::stepSimulationScatteredGrid(DT);
- #else
+#else
Boids::stepSimulationNaive(DT);
- #endif
+#endif
- #if VISUALIZE
+#if VISUALIZE
Boids::copyBoidsToVBO(dptrVertPositions, dptrVertVelocities);
- #endif
+#endif
// unmap buffer object
cudaGLUnmapBufferObject(boidVBO_positions);
cudaGLUnmapBufferObject(boidVBO_velocities);
- }
+}
- void mainLoop() {
+void mainLoop() {
double fps = 0;
double timebase = 0;
int frame = 0;
Boids::unitTest(); // LOOK-1.2 We run some basic example code to make sure
- // your CUDA development setup is ready to go.
+ // your CUDA development setup is ready to go.
while (!glfwWindowShouldClose(window)) {
- glfwPollEvents();
+ glfwPollEvents();
- frame++;
- double time = glfwGetTime();
+ frame++;
+ double time = glfwGetTime();
- if (time - timebase > 1.0) {
- fps = frame / (time - timebase);
- timebase = time;
- frame = 0;
- }
+ if (time - timebase > 1.0) {
+ fps = frame / (time - timebase);
+ timebase = time;
+ frame = 0;
+ }
- runCUDA();
+ runCUDA();
- std::ostringstream ss;
- ss << "[";
- ss.precision(1);
- ss << std::fixed << fps;
- ss << " fps] " << deviceName;
- glfwSetWindowTitle(window, ss.str().c_str());
+ std::ostringstream ss;
+ ss << "[";
+ ss.precision(1);
+ ss << std::fixed << fps;
+ ss << " fps] " << deviceName;
+ glfwSetWindowTitle(window, ss.str().c_str());
- glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+ glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
- #if VISUALIZE
- glUseProgram(program[PROG_BOID]);
- glBindVertexArray(boidVAO);
- glPointSize((GLfloat)pointSize);
- glDrawElements(GL_POINTS, N_FOR_VIS + 1, GL_UNSIGNED_INT, 0);
- glPointSize(1.0f);
+#if VISUALIZE
+ glUseProgram(program[PROG_BOID]);
+ glBindVertexArray(boidVAO);
+ glPointSize((GLfloat)pointSize);
+ glDrawElements(GL_POINTS, N_FOR_VIS + 1, GL_UNSIGNED_INT, 0);
+ glPointSize(1.0f);
- glUseProgram(0);
- glBindVertexArray(0);
+ glUseProgram(0);
+ glBindVertexArray(0);
- glfwSwapBuffers(window);
- #endif
+ glfwSwapBuffers(window);
+#endif
}
glfwDestroyWindow(window);
glfwTerminate();
- }
+}
- void errorCallback(int error, const char *description) {
+void errorCallback(int error, const char* description) {
fprintf(stderr, "error %d: %s\n", error, description);
- }
+}
- void keyCallback(GLFWwindow* window, int key, int scancode, int action, int mods) {
+void keyCallback(GLFWwindow* window, int key, int scancode, int action, int mods) {
if (key == GLFW_KEY_ESCAPE && action == GLFW_PRESS) {
- glfwSetWindowShouldClose(window, GL_TRUE);
+ glfwSetWindowShouldClose(window, GL_TRUE);
}
- }
+}
- void mouseButtonCallback(GLFWwindow* window, int button, int action, int mods) {
+void mouseButtonCallback(GLFWwindow* window, int button, int action, int mods) {
leftMousePressed = (button == GLFW_MOUSE_BUTTON_LEFT && action == GLFW_PRESS);
rightMousePressed = (button == GLFW_MOUSE_BUTTON_RIGHT && action == GLFW_PRESS);
- }
+}
- void mousePositionCallback(GLFWwindow* window, double xpos, double ypos) {
+void mousePositionCallback(GLFWwindow* window, double xpos, double ypos) {
if (leftMousePressed) {
- // compute new camera parameters
- phi += (xpos - lastX) / width;
- theta -= (ypos - lastY) / height;
- theta = std::fmax(0.01f, std::fmin(theta, 3.14f));
- updateCamera();
+ // compute new camera parameters
+ phi += (xpos - lastX) / width;
+ theta -= (ypos - lastY) / height;
+ theta = std::fmax(0.01f, std::fmin(theta, 3.14f));
+ updateCamera();
}
else if (rightMousePressed) {
- zoom += (ypos - lastY) / height;
- zoom = std::fmax(0.1f, std::fmin(zoom, 5.0f));
- updateCamera();
+ zoom += (ypos - lastY) / height;
+ zoom = std::fmax(0.1f, std::fmin(zoom, 5.0f));
+ updateCamera();
}
- lastX = xpos;
- lastY = ypos;
- }
+ lastX = xpos;
+ lastY = ypos;
+}
- void updateCamera() {
+void updateCamera() {
cameraPosition.x = zoom * sin(phi) * sin(theta);
cameraPosition.z = zoom * cos(theta);
cameraPosition.y = zoom * cos(phi) * sin(theta);
@@ -317,6 +335,6 @@ void initShaders(GLuint * program) {
glUseProgram(program[PROG_BOID]);
if ((location = glGetUniformLocation(program[PROG_BOID], "u_projMatrix")) != -1) {
- glUniformMatrix4fv(location, 1, GL_FALSE, &projection[0][0]);
+ glUniformMatrix4fv(location, 1, GL_FALSE, &projection[0][0]);
}
- }
+}
\ No newline at end of file