diff --git a/README.md b/README.md
index ee39093..a8e11e5 100644
--- a/README.md
+++ b/README.md
@@ -1,11 +1,69 @@
**University of Pennsylvania, CIS 5650: GPU Programming and Architecture,
Project 1 - Flocking**
-* (TODO) YOUR NAME HERE
- * (TODO) [LinkedIn](), [personal website](), [twitter](), etc.
-* Tested on: (TODO) Windows 22, i7-2222 @ 2.22GHz 22GB, GTX 222 222MB (Moore 2222 Lab)
+* Nadine Adnane
+ * [LinkedIn](https://www.linkedin.com/in/nadnane/)
+* Tested on my personal laptop (ASUS ROG Zephyrus M16):
+* **OS:** Windows 11
+* **Processor:** 12th Gen Intel(R) Core(TM) i9-12900H, 2500 Mhz, 14 Core(s), 20 Logical Processor(s)
+* **GPU:** NVIDIA GeForce RTX 3070 Ti Laptop GPU
-### (TODO: Your README)
+Note: I useed a late day for this assignment.
-Include screenshots, analysis, etc. (Remember, this is public, so don't put
-anything here that you don't want to share with the world.)
+### Results
+
+## Naive Implementation
+5000 boids
+
+
+50,000 boids
+
+
+100,000 boids
+
+
+## Uniform Grid
+5000 boids
+
+
+50,000 boids
+
+
+100,000 boids
+
+
+## Coherent Grid
+5000 boids
+
+
+50,000 boids
+
+
+100,000 boids
+
+
+## Graphs
+
+
+
+
+
+
+
+## Performance Analysis
+
+# How does changing the number of boids affect performance? Why do you think this is?
+
+Across all of the implementation methods, the performance decreases as the number of boids increases (as one would expect!). This makes sense since for each boid, we have to consider the impact of its neighboring boids on its velocity. More boids in the simulation means more neighbors to consider per boid, and thus more calculations to be done, especially for the naive implementation where we are not culling the pool of potential neighbors! For the naive implementation, there was a significant drop in the frame rate as the number of particles in the simulation increased. The decreased performance also makes sense for the grid-based methods, as even though we are culling the pool of potential neighbors, there are still more neighbors to consider per cell as the total number of boids increases (and thus more calculations to be performed!).
+
+# How does changing the block count and block size affect performance? Why do you think this is?
+
+Based off of my tests, there appeared to be a slight increase in performance when the block size is set to around 32, but otherwise it was somewhat difficult to tell if the difference in performance was particularly significant, as there wasn't much of a change in FPS during the tests that I ran. From there on, there seemed to be a point of diminishing returns. I think that increasing the block size causes an increase in performance at first because there are more threads available to be used.
+
+# For the coherent uniform grid: did you experience any performance improvements with the more coherent uniform grid? Was this the outcome you expected? Why or why not?
+
+Yes! There was a significant improvement in performance with the coherent uniform grid as compared to the scattered grid method, which was as I had expected. This makes sense, since the coherent grid method involves sorting the boid data (velocities and positions) to be contiguous in memory (and match the ordering of the grid index), so the data could be accessed more directly without the need for dev_particleArrayIndices.
+
+# Did changing cell width and checking 27 vs 8 neighboring cells affect performance? Why or why not? Be careful: it is insufficient (and possibly incorrect) to say that 27-cell is slower simply because there are more cells to check!
+
+Yes, changing the cell width and checking 27 cells instead of 8 did seem to affect the performance - I actually saw a very slight increase in performance when the cell width was changed from 1X to 2X and we were checking 27 neighboring cells instead of only 8. At first this was surprising - I assumed that checking 27 cells would surely be slower. However, I think the performance increase is due to the fact that even if we are only checking 8 cells, we still have to go through the process of figuring out which of the 8 cells to check. Also, since I only observed a slight increase in performance (at least from the FPS), I think more testing would be needed to actually be sure that the performance increase was significant.
\ No newline at end of file
diff --git a/images/coherent-100000.gif b/images/coherent-100000.gif
new file mode 100644
index 0000000..2a407d0
Binary files /dev/null and b/images/coherent-100000.gif differ
diff --git a/images/coherent-5000.gif b/images/coherent-5000.gif
new file mode 100644
index 0000000..1b9bd45
Binary files /dev/null and b/images/coherent-5000.gif differ
diff --git a/images/coherent-50000.gif b/images/coherent-50000.gif
new file mode 100644
index 0000000..6455a9a
Binary files /dev/null and b/images/coherent-50000.gif differ
diff --git a/images/naive-100000.gif b/images/naive-100000.gif
new file mode 100644
index 0000000..2574cab
Binary files /dev/null and b/images/naive-100000.gif differ
diff --git a/images/naive-5000.gif b/images/naive-5000.gif
new file mode 100644
index 0000000..6cc600e
Binary files /dev/null and b/images/naive-5000.gif differ
diff --git a/images/naive-50000.gif b/images/naive-50000.gif
new file mode 100644
index 0000000..eea3aa6
Binary files /dev/null and b/images/naive-50000.gif differ
diff --git a/images/uniform-100000.gif b/images/uniform-100000.gif
new file mode 100644
index 0000000..01fcf6c
Binary files /dev/null and b/images/uniform-100000.gif differ
diff --git a/images/uniform-5000.gif b/images/uniform-5000.gif
new file mode 100644
index 0000000..b8bbe8f
Binary files /dev/null and b/images/uniform-5000.gif differ
diff --git a/images/uniform-50000.gif b/images/uniform-50000.gif
new file mode 100644
index 0000000..cd213ea
Binary files /dev/null and b/images/uniform-50000.gif differ
diff --git a/src/kernel.cu b/src/kernel.cu
index 74dffcb..467f0cb 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -6,7 +6,7 @@
#include "utilityCore.hpp"
#include "kernel.h"
-// LOOK-2.1 potentially useful for doing grid-based neighbor search
+// LOOKED-2.1 potentially useful for doing grid-based neighbor search
#ifndef imax
#define imax( a, b ) ( ((a) > (b)) ? (a) : (b) )
#endif
@@ -39,7 +39,7 @@ void checkCUDAError(const char *msg, int line = -1) {
/*! Block size used for CUDA kernel launch. */
#define blockSize 128
-// LOOK-1.2 Parameters for the boids algorithm.
+// LOOKED-1.2 Parameters for the boids algorithm.
// These worked well in our reference implementation.
#define rule1Distance 5.0f
#define rule2Distance 3.0f
@@ -61,7 +61,7 @@ void checkCUDAError(const char *msg, int line = -1) {
int numObjects;
dim3 threadsPerBlock(blockSize);
-// LOOK-1.2 - These buffers are here to hold all your boid information.
+// LOOKED-1.2 - These buffers are here to hold all your boid information.
// These get allocated for you in Boids::initSimulation.
// Consider why you would need two velocity buffers in a simulation where each
// boid cares about its neighbors' velocities.
@@ -70,7 +70,7 @@ 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
+// LOOKED-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.
@@ -83,10 +83,12 @@ thrust::device_ptr dev_thrust_particleGridIndices;
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
+// DONE-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_sorted;
+glm::vec3* dev_vel1_sorted;
-// LOOK-2.1 - Grid parameters based on simulation parameters.
+// LOOKED-2.1 - Grid parameters based on simulation parameters.
// These are automatically computed for you in Boids::initSimulation
int gridCellCount;
int gridSideCount;
@@ -109,7 +111,7 @@ __host__ __device__ unsigned int hash(unsigned int a) {
}
/**
-* LOOK-1.2 - this is a typical helper function for a CUDA kernel.
+* LOOKED-1.2 - this is a typical helper function for a CUDA kernel.
* Function for generating a random vec3.
*/
__host__ __device__ glm::vec3 generateRandomVec3(float time, int index) {
@@ -120,7 +122,7 @@ __host__ __device__ glm::vec3 generateRandomVec3(float time, int index) {
}
/**
-* LOOK-1.2 - This is a basic CUDA kernel.
+* LOOKED-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) {
@@ -140,7 +142,7 @@ 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.
+ // LOOKED-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!");
@@ -151,12 +153,12 @@ void Boids::initSimulation(int N) {
cudaMalloc((void**)&dev_vel2, N * sizeof(glm::vec3));
checkCUDAErrorWithLine("cudaMalloc dev_vel2 failed!");
- // LOOK-1.2 - This is a typical CUDA kernel invocation.
+ // LOOKED-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
+ // LOOKED-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;
@@ -168,7 +170,28 @@ void Boids::initSimulation(int N) {
gridMinimum.y -= halfGridWidth;
gridMinimum.z -= halfGridWidth;
- // TODO-2.1 TODO-2.3 - Allocate additional buffers here.
+ // DONE-2.1 DONE-2.3 - Allocate additional buffers here.
+ cudaMalloc((void**)&dev_particleArrayIndices, N * sizeof(int));
+ checkCUDAErrorWithLine("cudaMalloc dev_particleArrayIndices failed!");
+
+ cudaMalloc((void**)&dev_particleGridIndices, N * sizeof(int));
+ checkCUDAErrorWithLine("cudaMalloc dev_particleGridIndices failed!");
+
+ cudaMalloc((void**)&dev_gridCellStartIndices, gridCellCount * sizeof(int));
+ checkCUDAErrorWithLine("cudaMalloc dev_gridCellStartIndices failed!");
+
+ cudaMalloc((void**)&dev_gridCellEndIndices, gridCellCount * sizeof(int));
+ checkCUDAErrorWithLine("cudaMalloc dev_gridCellEndIndices failed!");
+
+ cudaMalloc((void**)&dev_pos_sorted, N * sizeof(glm::vec3));
+ checkCUDAErrorWithLine("cudaMalloc dev_pos_sorted failed!");
+
+ cudaMalloc((void**)&dev_vel1_sorted, N * sizeof(glm::vec3));
+ checkCUDAErrorWithLine("cudaMalloc dev_vel1_sorted failed!");
+
+ dev_thrust_particleArrayIndices = thrust::device_ptr(dev_particleArrayIndices);
+ dev_thrust_particleGridIndices = thrust::device_ptr(dev_particleGridIndices);
+
cudaDeviceSynchronize();
}
@@ -224,31 +247,111 @@ void Boids::copyBoidsToVBO(float *vbodptr_positions, float *vbodptr_velocities)
******************/
/**
-* LOOK-1.2 You can use this as a helper for kernUpdateVelocityBruteForce.
+* LOOKED-1.2 You can use this as a helper for kernUpdateVelocityBruteForce.
* __device__ code can be called from a __global__ context
* 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);
+
+ glm::vec3 perceivedCenter(0.f, 0.f, 0.f);
+ glm::vec3 perceivedVelocity(0.f, 0.f, 0.f);
+
+ glm::vec3 delta1(0.f, 0.f, 0.f), delta2(0.f, 0.f, 0.f), delta3(0.f, 0.f, 0.f);
+
+ // Keep track of the number of neighbors for rules 1 and 3
+ int rule1Neighbors = 0;
+ int rule3Neighbors = 0;
+
+ // Final delta V
+ glm::vec3 result;
+
+ // Loop through all of the boids and exclude current boid
+ glm::vec3 c = glm::vec3(0, 0, 0);
+ for(int i = 0; i < N; i++)
+ {
+ if (i == iSelf)
+ {
+ // Exclude the current boid
+ continue;
+ }
+
+ glm::vec3 neighborBoidPos = pos[i];
+ float distance = glm::length(neighborBoidPos - pos[iSelf]);
+
+ // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves
+ if (distance < rule1Distance)
+ {
+ perceivedCenter += neighborBoidPos;
+ rule1Neighbors++;
+ }
+
+ // Rule 2: boids try to stay a distance d away from each other
+ if (distance < rule2Distance)
+ {
+ c -= (pos[i] - pos[iSelf]);
+ }
+
+ // Rule 3: boids try to match the speed of surrounding boids
+ if (distance < rule3Distance)
+ {
+ perceivedVelocity += vel[i];
+ rule3Neighbors++;
+ }
+
+ }
+
+ // Check if any of the neighbors should influence the current boid
+ // For rules 1 and 3
+ if (rule1Neighbors > 0)
+ {
+ delta1 = (perceivedCenter / (float) rule1Neighbors - pos[iSelf]) * rule1Scale;
+ }
+
+ delta2 = c * rule2Scale;
+
+ if (rule3Neighbors > 0)
+ {
+ delta3 = (perceivedVelocity / (float) rule3Neighbors) * rule3Scale;
+ }
+
+ result += delta1 + delta2 + delta3;
+
+ return result;
}
/**
-* TODO-1.2 implement basic flocking
+* DONE-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) {
+
+ // First figure out which boid we are working with
+ int index = threadIdx.x + (blockIdx.x * blockDim.x);
+ if (index >= N)
+ {
+ return;
+ }
+
// Compute a new velocity based on pos and vel1
+ // vel1[index] represents current velocity
+ glm::vec3 deltaV = computeVelocityChange(N, index, pos, vel1);
+ glm::vec3 newVel = vel1[index] + deltaV;
+
// Clamp the speed
+ if (glm::length(newVel) > maxSpeed)
+ {
+ newVel = glm::normalize(newVel) * maxSpeed;
+ }
+
// Record the new velocity into vel2. Question: why NOT vel1?
+ // Answer: Because we don't want to write over the data we're reading from (I.e. we need neighbor's data)
+ vel2[index] = newVel;
}
/**
-* LOOK-1.2 Since this is pretty trivial, we implemented it for you.
+* LOOKED-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) {
@@ -272,8 +375,8 @@ __global__ void kernUpdatePos(int N, float dt, glm::vec3 *pos, glm::vec3 *vel) {
pos[index] = thisPos;
}
-// LOOK-2.1 Consider this method of computing a 1D index from a 3D grid index.
-// LOOK-2.3 Looking at this method, what would be the most memory efficient
+// LOOKED-2.1 Consider this method of computing a 1D index from a 3D grid index.
+// LOOKED-2.3 Looking at this method, what would be the most memory efficient
// order for iterating over neighboring grid cells?
// for(x)
// for(y)
@@ -285,13 +388,27 @@ __device__ int gridIndex3Dto1D(int x, int y, int z, int gridResolution) {
__global__ void kernComputeIndices(int N, int gridResolution,
glm::vec3 gridMin, float inverseCellWidth,
glm::vec3 *pos, int *indices, int *gridIndices) {
- // TODO-2.1
+ // DONE-2.1
+ // First figure out which boid we are working with
+ int index = threadIdx.x + (blockIdx.x * blockDim.x);
+ if (index >= N)
+ {
+ return;
+ }
+
// - Label each boid with the index of its grid cell.
+ // Calculate the index of the grid cell using the boid's position
+ glm::vec3 offset = pos[index] - gridMin;
+ glm::vec3 cellIndex = glm::floor(inverseCellWidth * offset);
+ // Label this boid in the grid index table
+ gridIndices[index] = gridIndex3Dto1D(cellIndex.x, cellIndex.y, cellIndex.z, gridResolution);
+
// - Set up a parallel array of integer indices as pointers to the actual
// boid data in pos and vel1/vel2
+ indices[index] = index;
}
-// LOOK-2.1 Consider how this could be useful for indicating that a cell
+// LOOKED-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;
@@ -302,10 +419,35 @@ __global__ void kernResetIntBuffer(int N, int *intBuffer, int value) {
__global__ void kernIdentifyCellStartEnd(int N, int *particleGridIndices,
int *gridCellStartIndices, int *gridCellEndIndices) {
- // TODO-2.1
+ // DONE-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!"
+
+ // First figure out which particle we are working with
+ int index = threadIdx.x + (blockIdx.x * blockDim.x);
+ if (index >= N) {
+ return;
+ }
+
+ // Find the grid index for this particle
+ int particleGridIndex = particleGridIndices[index];
+
+ // Find the start point for each cell in the gridIndices array
+ // If this is the first index or if this index doesn't match the one before it
+ // Then set this as a new cell!
+ if ((index == 0) || (particleGridIndex != particleGridIndices[index - 1]))
+ {
+ gridCellStartIndices[particleGridIndex] = index;
+ }
+
+ // Find the end point for each cell in the gridIndices array
+ // If this is the last index, or if this index doesn't match the next one,
+ // Then set this as the end of the current cell!
+ if ((index == N - 1) || (particleGridIndex != particleGridIndices[index + 1]))
+ {
+ gridCellEndIndices[particleGridIndex] = index;
+ }
}
__global__ void kernUpdateVelNeighborSearchScattered(
@@ -314,14 +456,129 @@ __global__ void kernUpdateVelNeighborSearchScattered(
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
+ // DONE-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
+
+ // First find the current boid's index
+ int index = threadIdx.x + (blockIdx.x * blockDim.x);
+ if (index >= N)
+ {
+ return;
+ }
+
+ // Find the max neighborhood distance
+ float neighborhoodDistance = fmax(fmax(rule1Distance, rule2Distance), rule3Distance);
+
+ // Figure out which grid cell this boid is in
+ glm::vec3 boidPos = pos[index];
+ glm::vec3 gridPos = boidPos - gridMin;
+ int gridIdx = gridIndex3Dto1D(
+ gridPos.x * inverseCellWidth,
+ gridPos.y * inverseCellWidth,
+ gridPos.z * inverseCellWidth,
+ gridResolution);
+
// - Identify which cells may contain neighbors. This isn't always 8.
+ int minX = imax(0, static_cast((gridPos.x - neighborhoodDistance) * inverseCellWidth));
+ int maxX = imin(gridResolution - 1, static_cast((gridPos.x + neighborhoodDistance) * inverseCellWidth));
+
+ int minY = imax(0, static_cast((gridPos.y - neighborhoodDistance) * inverseCellWidth));
+ int maxY = imin(gridResolution - 1, static_cast((gridPos.y + neighborhoodDistance) * inverseCellWidth));
+
+ int minZ = imax(0, static_cast((gridPos.z - neighborhoodDistance) * inverseCellWidth));
+ int maxZ = imin(gridResolution - 1, static_cast((gridPos.z + neighborhoodDistance) * inverseCellWidth));
+
// - For each cell, read the start/end indices in the boid pointer array.
+ glm::vec3 c(0.f, 0.f, 0.f);
+ glm::vec3 deltaV(0.f, 0.f, 0.f);
+ glm::vec3 perceivedCenter(0.f, 0.f, 0.f);
+ glm::vec3 perceivedVelocity(0.f, 0.f, 0.f);
+
+ // Keep track of the number of neighbors for rules 1 and 3
+ int rule1Neighbors = 0;
+ int rule3Neighbors = 0;
+
// - Access each boid in the cell and compute velocity change from
// the boids rules, if this boid is within the neighborhood distance.
+ for (int z = minZ; z <= maxZ; z++)
+ {
+ for (int y = minY; y <= maxY; y++)
+ {
+ for (int x = minX; x <= maxX; x++)
+ {
+ glm::vec3 cell(x * cellWidth, y * cellWidth, z * cellWidth);
+
+ // Find the nearest point on current cell to our boid
+ glm::vec3 nearest;
+ nearest.x = fmax(cell.x, fmin(gridPos.x, cell.x + cellWidth));
+ nearest.y = fmax(cell.y, fmin(gridPos.y, cell.y + cellWidth));
+ nearest.z = fmax(cell.z, fmin(gridPos.z, cell.z + cellWidth));
+
+ // If the nearest point is within neighborhood distance, consider the neighb boids in that cell
+ if (glm::distance(gridPos, nearest) <= neighborhoodDistance) {
+ int startEndIdx = gridIndex3Dto1D(x, y, z, gridResolution);
+
+ int gridCellStart = gridCellStartIndices[startEndIdx];
+ int gridCellEnd = gridCellEndIndices[startEndIdx];
+
+ if (gridCellStart == -1) {
+ continue;
+ }
+
+ for (int cellIndex = gridCellStart; cellIndex <= gridCellEnd; cellIndex++) {
+ int boidIdx = particleArrayIndices[cellIndex];
+ glm::vec3 bPos = pos[boidIdx];
+ glm::vec3 bVel = vel1[boidIdx];
+
+ if (boidIdx != index) {
+
+ // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves
+ if (glm::distance(bPos, boidPos) < rule1Distance) {
+ rule1Neighbors++;
+ perceivedCenter += bPos;
+ }
+
+ // Rule 2: boids try to stay a distance d away from each other
+ if (glm::distance(bPos, boidPos) < rule2Distance) {
+ c -= (bPos - boidPos);
+ }
+
+ // Rule 3: boids try to match the speed of surrounding boids
+ if (glm::distance(bPos, boidPos) < rule3Distance) {
+ rule3Neighbors++;
+ perceivedVelocity += bVel;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ glm::vec3 rule1(0.f);
+ glm::vec3 rule2(0.f);
+ glm::vec3 rule3(0.f);
+
+ if (rule1Neighbors > 0) {
+ perceivedCenter /= glm::vec3(rule1Neighbors);
+ rule1 = (perceivedCenter - boidPos) * rule1Scale;
+ }
+
+ rule2 = c * rule2Scale;
+
+ if (rule3Neighbors > 0) {
+ perceivedVelocity /= glm::vec3(rule3Neighbors);
+ rule3 = perceivedVelocity * rule3Scale;
+ }
+
// - Clamp the speed change before putting the new speed in vel2
+ glm::vec3 newV = vel1[index] + (rule1 + rule2 + rule3);
+ if (glm::length(newV) > maxSpeed) {
+ newV = glm::normalize(newV) * maxSpeed;
+ }
+
+ vel2[index] = newV;
+
}
__global__ void kernUpdateVelNeighborSearchCoherent(
@@ -329,7 +586,7 @@ __global__ void kernUpdateVelNeighborSearchCoherent(
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,
+ // DONE-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.
@@ -341,33 +598,248 @@ __global__ void kernUpdateVelNeighborSearchCoherent(
// - Access each boid in the cell and compute velocity change from
// the boids rules, if this boid is within the neighborhood distance.
// - Clamp the speed change before putting the new speed in vel2
+
+ // First find the current boid's index
+ int index = threadIdx.x + (blockIdx.x * blockDim.x);
+ if (index >= N)
+ {
+ return;
+ }
+
+ // Figure out which grid cell this boid is in
+ glm::vec3 boidPos = pos[index];
+ glm::vec3 gridPos = boidPos - gridMin;
+ int gridIdx = gridIndex3Dto1D(
+ gridPos.x * inverseCellWidth,
+ gridPos.y * inverseCellWidth,
+ gridPos.z * inverseCellWidth,
+ gridResolution);
+
+ // Find the max neighborhood distance
+ float neighborhoodDistance = fmax(fmax(rule1Distance, rule2Distance), rule3Distance);
+
+ // - Identify which cells may contain neighbors. This isn't always 8.
+ // Calculate the mins in the x, y, and z directions
+ int minX = imax(0, static_cast((gridPos.x - neighborhoodDistance) * inverseCellWidth));
+ int minY = imax(0, static_cast((gridPos.y - neighborhoodDistance) * inverseCellWidth));
+ int minZ = imax(0, static_cast((gridPos.z - neighborhoodDistance) * inverseCellWidth));
+ // Calculate the maxes in the x, y, and z directions
+ int maxX = imin(gridResolution - 1, static_cast((gridPos.x + neighborhoodDistance) * inverseCellWidth));
+ int maxY = imin(gridResolution - 1, static_cast((gridPos.y + neighborhoodDistance) * inverseCellWidth));
+ int maxZ = imin(gridResolution - 1, static_cast((gridPos.z + neighborhoodDistance) * inverseCellWidth));
+
+ // - For each cell, read the start/end indices in the boid pointer array.
+ glm::vec3 c(0.f, 0.f, 0.f);
+ glm::vec3 deltaV(0.f, 0.f, 0.f);
+ glm::vec3 perceivedCenter(0.f, 0.f, 0.f);
+ glm::vec3 perceivedVelocity(0.f, 0.f, 0.f);
+
+ // Keep track of the number of neighbors for rules 1 and 3
+ int rule1Neighbors = 0;
+ int rule3Neighbors = 0;
+
+ // - Access each boid in the cell and compute velocity change from
+ // the boids rules, if this boid is within the neighborhood distance.
+ for (int z = minZ; z <= maxZ; z++)
+ {
+ for (int y = minY; y <= maxY; y++)
+ {
+ for (int x = minX; x <= maxX; x++)
+ {
+ glm::vec3 cell(x * cellWidth, y * cellWidth, z * cellWidth);
+
+ // Find the nearest point on current cell to our boid
+ glm::vec3 nearest;
+ nearest.x = fmax(cell.x, fmin(gridPos.x, cell.x + cellWidth));
+ nearest.y = fmax(cell.y, fmin(gridPos.y, cell.y + cellWidth));
+ nearest.z = fmax(cell.z, fmin(gridPos.z, cell.z + cellWidth));
+
+ if (glm::distance(gridPos, nearest) <= neighborhoodDistance)
+ {
+ // Find the grid start and end indices
+ int startEndIndex = gridIndex3Dto1D(x, y, z, gridResolution);
+
+ // Only proceed with a valid start
+ int gridCellStart = gridCellStartIndices[startEndIndex];
+ if (gridCellStart == -1) {
+ continue;
+ }
+
+ // Find the grid cell end
+ int gridCellEnd = gridCellEndIndices[startEndIndex];
+
+ // Iterate through the neighboids (get it? neighboring boids? :D)
+ for (int neighboidIndex = gridCellStart; neighboidIndex <= gridCellEnd; neighboidIndex++) {
+
+ // Get the position and velocity for the current boid
+ glm::vec3 neighboidPos = pos[neighboidIndex];
+ glm::vec3 neighboidVel = vel1[neighboidIndex];
+
+ // Exclude the current boid, we only care about neighbors
+ if (neighboidIndex == index)
+ {
+ continue;
+ }
+
+ // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves
+ if (glm::distance(boidPos, boidPos) < rule1Distance) {
+ rule1Neighbors++;
+ perceivedCenter += boidPos;
+ }
+
+ // Rule 2: boids try to stay a distance d away from each other
+ if (glm::distance(boidPos, neighboidPos) < rule2Distance) {
+ c -= (neighboidPos - boidPos);
+ }
+
+ // Rule 3: boids try to match the speed of surrounding boids
+ if (glm::distance(neighboidPos, boidPos) < rule3Distance) {
+ rule3Neighbors++;
+ perceivedVelocity += neighboidVel;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ glm::vec3 delta1(0.f), delta2(0.f), delta3(0.f);
+
+ // Rule 1
+ if (rule1Neighbors > 0) {
+ perceivedCenter /= glm::vec3(rule1Neighbors);
+ delta1 = (perceivedCenter - boidPos) * rule1Scale;
+ }
+
+ // Rule 2
+ delta2 = c * rule2Scale;
+
+ // Rule3
+ if (rule3Neighbors > 0) {
+ perceivedVelocity /= glm::vec3(rule3Neighbors);
+ delta3 = perceivedVelocity * rule3Scale;
+ }
+
+ // - Clamp the speed change before putting the new speed in vel2
+ glm::vec3 newV = vel1[index] + (delta1 + delta2 + delta3);
+ if (glm::length(newV) > maxSpeed)
+ {
+ newV = glm::normalize(newV) * maxSpeed;
+ }
+
+ vel2[index] = newV;
}
/**
* 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
+ // DONE-1.2 - use the kernels you wrote to step the simulation forward in time.
+ dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize);
+
+ // Update the boid's velocity
+ kernUpdateVelocityBruteForce << > > (numObjects, dev_pos, dev_vel1, dev_vel2);
+
+ // Update the boid's position
+ kernUpdatePos << > > (numObjects, dt, dev_pos, dev_vel2);
+
+ // DONE-1.2 ping-pong the velocity buffers
+ glm::vec3* temp = dev_vel2;
+ dev_vel2 = dev_vel1;
+ dev_vel1 = temp;
}
void Boids::stepSimulationScatteredGrid(float dt) {
- // TODO-2.1
+ // DONE-2.1
// Uniform Grid Neighbor search using Thrust sort.
+ // Calculate blocks per grid
+ int blocksPerGrid = (numObjects + threadsPerBlock.x - 1) / threadsPerBlock.x;
+
// In Parallel:
// - label each particle with its array index as well as its grid index.
// Use 2x width grids.
+ kernComputeIndices << > > (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, dev_pos, dev_particleArrayIndices, dev_particleGridIndices);
+
// - Unstable key sort using Thrust. A stable sort isn't necessary, but you
// are welcome to do a performance comparison.
+ thrust::sort_by_key(dev_thrust_particleGridIndices, dev_thrust_particleGridIndices + numObjects, dev_thrust_particleArrayIndices);
+
// - Naively unroll the loop for finding the start and end indices of each
// cell's data pointers in the array of boid indices
+ dim3 gridCellBlocksPerGrid((gridCellCount + threadsPerBlock.x - 1) / threadsPerBlock.x);
+
+ kernResetIntBuffer << > > (gridCellCount, dev_gridCellStartIndices, -1);
+ checkCUDAErrorWithLine("kernResetIntBuffer failed!");
+ kernResetIntBuffer << > > (gridCellCount, dev_gridCellEndIndices, -1);
+ checkCUDAErrorWithLine("kernResetIntBuffer failed!");
+
+ kernIdentifyCellStartEnd << > > (numObjects, dev_particleGridIndices, dev_gridCellStartIndices, dev_gridCellEndIndices);
+ checkCUDAErrorWithLine("kernIdentifyCellStartEnd failed!");
+
+#if 0
+ std::vector host_particleGridIndices;
+ std::vector host_particleArrayIndices;
+ std::vector host_start;
+ std::vector host_end;
+
+ host_particleGridIndices.resize(numObjects);
+ host_particleArrayIndices.resize(numObjects);
+ host_start.resize(gridCellCount);
+ host_end.resize(gridCellCount);
+
+ cudaMemcpy(host_particleGridIndices.data(), dev_particleGridIndices, numObjects * sizeof(int), cudaMemcpyDeviceToHost);
+ checkCUDAErrorWithLine("Memcpy host_particleGridIndices.data failed!");
+
+ cudaMemcpy(host_particleArrayIndices.data(), dev_particleArrayIndices, numObjects * sizeof(int), cudaMemcpyDeviceToHost);
+ checkCUDAErrorWithLine("Memcpy host_particleArrayIndices.data failed!");
+
+ cudaMemcpy(host_start.data(), dev_gridCellStartIndices, gridCellCount * sizeof(int), cudaMemcpyDeviceToHost);
+ checkCUDAErrorWithLine("Memcpy start.data failed!");
+
+ cudaMemcpy(host_end.data(), dev_gridCellEndIndices, gridCellCount * sizeof(int), cudaMemcpyDeviceToHost);
+ checkCUDAErrorWithLine("Memcpy end.data failed!");
+#endif
+
// - Perform velocity updates using neighbor search
+ kernUpdateVelNeighborSearchScattered << > > (
+ numObjects,
+ gridSideCount,
+ gridMinimum,
+ gridInverseCellWidth,
+ gridCellWidth,
+ dev_gridCellStartIndices,
+ dev_gridCellEndIndices,
+ dev_particleArrayIndices,
+ dev_pos,
+ dev_vel1,
+ dev_vel2
+ );
+ checkCUDAErrorWithLine("kernUpdateVelNeighborSearchScattered failed!");
+
// - Update positions
+ kernUpdatePos << > > (numObjects, dt, dev_pos, dev_vel2);
+ checkCUDAErrorWithLine("kernUpdatePos failed!");
+
// - Ping-pong buffers as needed
+ cudaMemcpy(dev_vel1, dev_vel2, numObjects * sizeof(glm::vec3), cudaMemcpyDeviceToDevice);
+ checkCUDAErrorWithLine("Memcpy (dev_vel2 to dev_vel1) failed!");
+}
+
+// Helper function to copy the position and velocity into new buffers with sorted indices based on grid index
+__global__ void kernReshufflePosAndVel(
+ int N, glm::vec3* pos, glm::vec3* vel, glm::vec3* pos_sorted,
+ glm::vec3* vel_sorted, int* particleArrayIndices) {
+ int index = (blockIdx.x * blockDim.x) + threadIdx.x;
+ if (index >= N) {
+ return;
+ }
+
+ pos_sorted[index] = pos[particleArrayIndices[index]];
+ vel_sorted[index] = vel[particleArrayIndices[index]];
}
void Boids::stepSimulationCoherentGrid(float dt) {
- // TODO-2.3 - start by copying Boids::stepSimulationNaiveGrid
+ // DONE-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.
@@ -382,6 +854,58 @@ void Boids::stepSimulationCoherentGrid(float dt) {
// - Perform velocity updates using neighbor search
// - Update positions
// - Ping-pong buffers as needed. THIS MAY BE DIFFERENT FROM BEFORE.
+ //==================================================================================
+
+ // First calculate the number of blocks per grid
+ int blocksPerGrid = (numObjects + threadsPerBlock.x - 1) / threadsPerBlock.x;
+
+ // - Label each particle with its array index as well as its grid index. Use 2x width grids
+ kernComputeIndices << > > (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, dev_pos, dev_particleArrayIndices, dev_particleGridIndices);
+ checkCUDAErrorWithLine("kernComputeIndices failed!");
+
+ // - Unstable key sort using Thrust. A stable sort isn't necessary, but you are welcome to do a performance comparison.
+ thrust::sort_by_key(dev_thrust_particleGridIndices, dev_thrust_particleGridIndices + numObjects, dev_thrust_particleArrayIndices);
+ kernReshufflePosAndVel << > > (numObjects, dev_pos, dev_vel1, dev_pos_sorted, dev_vel1_sorted, dev_particleArrayIndices);
+
+ // - Naively unroll the loop for finding the start and end indices of each cell's data pointers in the array of boid indices
+ dim3 gridCellBlocksPerGrid((gridCellCount + threadsPerBlock.x - 1) / threadsPerBlock.x);
+ // Reset start index buffer
+ kernResetIntBuffer << > > (gridCellCount, dev_gridCellStartIndices, -1);
+ checkCUDAErrorWithLine("kernResetIntBuffer failed!");
+ // Reset ending index buffer
+ kernResetIntBuffer << > > (gridCellCount, dev_gridCellEndIndices, -1);
+ checkCUDAErrorWithLine("kernResetIntBuffer failed!");
+ // Run kernel to try and identify cell start and end indices
+ kernIdentifyCellStartEnd << > > (numObjects, dev_particleGridIndices, dev_gridCellStartIndices, dev_gridCellEndIndices);
+ checkCUDAErrorWithLine("kernIdentifyCellStartEnd failed!");
+
+ // - 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
+
+ kernUpdateVelNeighborSearchCoherent << > > (
+ numObjects,
+ gridSideCount,
+ gridMinimum,
+ gridInverseCellWidth,
+ gridCellWidth,
+ dev_gridCellStartIndices,
+ dev_gridCellEndIndices,
+ dev_pos_sorted,
+ dev_vel1_sorted,
+ dev_vel2
+ );
+
+ // - Update positions
+ kernUpdatePos << > > (numObjects, dt, dev_pos_sorted, dev_vel2);
+ checkCUDAErrorWithLine("kernUpdatePos failed!");
+
+ // - Ping-pong buffers as needed. THIS MAY BE DIFFERENT FROM BEFORE.
+ cudaMemcpy(dev_pos, dev_pos_sorted, numObjects * sizeof(glm::vec3), cudaMemcpyDeviceToDevice);
+ checkCUDAErrorWithLine("Memcpy (dev_pos_sorted to dev_pos) failed!");
+
+ cudaMemcpy(dev_vel1, dev_vel2, numObjects * sizeof(glm::vec3), cudaMemcpyDeviceToDevice);
+ checkCUDAErrorWithLine("Memcpy (dev_vel2 to dev_vel1) failed!");
}
void Boids::endSimulation() {
@@ -389,11 +913,16 @@ void Boids::endSimulation() {
cudaFree(dev_vel2);
cudaFree(dev_pos);
- // TODO-2.1 TODO-2.3 - Free any additional buffers here.
+ // DONE-2.1 DONE-2.3 - Free any additional buffers here.
+ cudaFree(dev_gridCellStartIndices);
+ cudaFree(dev_gridCellEndIndices);
+ cudaFree(dev_particleArrayIndices);
+ cudaFree(dev_particleGridIndices);
}
void Boids::unitTest() {
- // LOOK-1.2 Feel free to write additional tests here.
+ // LOOKED-1.2 Feel free to write additional tests here.
+
// test unstable sort
int *dev_intKeys;
@@ -435,7 +964,7 @@ void Boids::unitTest() {
// 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
+ // LOOKED-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
@@ -449,7 +978,7 @@ void Boids::unitTest() {
std::cout << " value: " << intValues[i] << std::endl;
}
- // cleanup
+ // Cleanup
cudaFree(dev_intKeys);
cudaFree(dev_intValues);
checkCUDAErrorWithLine("cudaFree failed!");
diff --git a/src/main.cpp b/src/main.cpp
index fe657ed..10b54f5 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -15,13 +15,13 @@
// Configuration
// ================
-// LOOK-2.1 LOOK-2.3 - toggles for UNIFORM_GRID and COHERENT_GRID
+// LOOKED-2.1 LOOKED-2.3 - toggles for UNIFORM_GRID and COHERENT_GRID
#define VISUALIZE 1
-#define UNIFORM_GRID 0
-#define COHERENT_GRID 0
+#define UNIFORM_GRID 1
+#define COHERENT_GRID 1
-// LOOK-1.2 - change this to adjust particle count in the simulation
-const int N_FOR_VIS = 5000;
+// LOOKED-1.2 - change this to adjust particle count in the simulation
+const int N_FOR_VIS = 100000;
const float DT = 0.2f;
/**
@@ -220,7 +220,7 @@ void initShaders(GLuint * program) {
double timebase = 0;
int frame = 0;
- Boids::unitTest(); // LOOK-1.2 We run some basic example code to make sure
+ Boids::unitTest(); // LOOKED-1.2 We run some basic example code to make sure
// your CUDA development setup is ready to go.
while (!glfwWindowShouldClose(window)) {
diff --git a/src/main.hpp b/src/main.hpp
index faf87cf..11a3b8c 100644
--- a/src/main.hpp
+++ b/src/main.hpp
@@ -32,7 +32,7 @@ const unsigned int PROG_BOID = 0;
const float fovy = (float) (PI / 4);
const float zNear = 0.10f;
const float zFar = 10.0f;
-// LOOK-1.2: for high DPI displays, you may want to double these settings.
+// LOOKED-1.2: for high DPI displays, you may want to double these settings.
int width = 1280;
int height = 720;
int pointSize = 2;