diff --git a/README.md b/README.md
index ee39093..bd2e537 100644
--- a/README.md
+++ b/README.md
@@ -1,11 +1,85 @@
**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)
+* Matt Schwartz
+ * [LinkedIn](https://www.linkedin.com/in/matthew-schwartz-37019016b/)
+ * [Personal website](https://mattzschwartz.web.app/)
+* Tested on: Windows 10 22H2, Intel(R) Core(TM) i7-10750H CPU @ 2.60GHz, NVIDIA GeForce RTX 2060
-### (TODO: Your README)
+# Reynold Boids Flocking
-Include screenshots, analysis, etc. (Remember, this is public, so don't put
-anything here that you don't want to share with the world.)
+
+
+
+
+## Background
+
+This repository contains several implementations of Reynolds' famous Boids algorithm, in which multitudes of particles are simulated interacting with each other via three simple rules:
+
+1. Cohesion - each Boid steers itself towards the perceived center of mass of the group.
+1. Separation - each Boid maintains some small distance away from other Boids (to prevent crashes).
+1. Alignment - each Boid attempts to align its velocity with the perceived group velocity.
+
+Following only these rules, and tuning weights for each rule carefully, we can achieve a group motion that resembles flocking (of birds, fish, etc.), as seen above.
+
+Because each Boid can be simulated at a given timestep independently of all other Boids, this algorithm lends itself well to parallelization via the GPU. In this project, I have implemented the three different GPU-based approaches as follows:
+
+
+### Naive (Brute-force)
+
+In this approach, the kernel which simulates each Boid does so by comparing a given Boid against every other Boid (more or less - a target Boid could be out of range, and its impact excluded, but is still considered nonetheless).
+
+This approach has the advantage of simplicity; no extra data structures beyond the positions and velocities of all Boids are necessary to optimize this approach. The downside, of course, is the time complexity required. As the number of Boids in the simulation grows, the simulation time grows with its square.
+
+### Scattered Uniform Grid
+
+Since each of our aforementioned rules has an associated distance of effect, we don't *need* to compare every Boid to every other. We only need to compare Boids within a given vicinity. Thus, the second approach buckets Boids into cells of a uniform spatial grid (see image below). With a well-chosen grid cell width, we can limit our neighboring Boid checks to only those within the surrounding cells.
+
+In theory, this optimization should increase performance, as we have fewer Boids to check against. In practice, we'll see below that the performance gains significant, but that there is still much room for improvement. The data structures used in the approach, combined with a need for a preprocessing sorting step, lead to scattered global memory accesses, whose latency offsets some of the gains from the use of a uniform grid.
+
+
+
+### Coherent Uniform Grid
+
+We can improve on the scattered grid approach by changing the underlying data structures involved. By sorting the Boid position and velocity arrays on each timestep, and cutting out a middleman look-up table into these arrays, each Boid's information can be stored and accessed in contiguous global memory.
+
+The downsides to this approach are the code complexity involved and extra memory usage for book-keeping buffers, but the resulting performance gains are notable.
+
+
+## Performance Analysis
+
+In this section, I will analyze how each of the above approaches performs with respect to two variables: number of Boids simulated, and number of threads per GPU block. The metric of performance I use will be average framerate; the greater the framerate, the better the performance of the implementation. During measurements, visualation will be turned off so that the framerate is only a factor of the simulation speed. V-sync will also be turned off.
+
+### Performance as a function of number of Boids simulated
+
+
+
+
+
+Unsurpsingly, increasing the number of Boids in the simulation decreases the framerate across all implementations. This makes sense, as the GPU has more work to do as more Boids are simulated, which will slow it down. It also appears that the naive implementation suffers most drastically - as is evident from the concavity of its curve. This also makes sense, since the Naive algorithm has O(N^2) time complexity, so it will suffer the most of all implementations.
+
+Another trend of interest: the coherent uniform grid outperforms the scattered uniform grid across all numbers of Boids. However, the difference in performance is vanishing as the number of Boids decreases. I have three potential theories to explain this (none of which are mutually exclusive):
+
+1. With a small enough number of Boids, a larger portion of Boid data can fit in cache memory / fewer cache misses, negating the impact of random global memory access.
+1. With a small enough number of Boids, a larger portion of Boids can be simulated in parallel, so the latency of global memory access is less noticeable.
+1. The coherent grid requires an extra sort step, which entails overhead. The cost of this extra sort is more noticeable with a smaller number of boids.
+
+### Performance as a function of number of threads per block
+
+
+
+
+
+Little to talk about here - (perhaps) surprisingly, the number of threads per block has very little impact on the performance of the simulation. The most apparent answer here, is that the GPU is already saturated at the minimum tested 64 threads/block, and so increasing the number of threads per block any further does not help performance.
+
+### Bonus: variable cell width
+
+
+
+
+
+This graph analyzes the impact of the uniform grid's cell width on performance. In the above plot, cell width is a ratio where the denominator is the maximum distance by which a Boid will be influenced by another Boid. At a ratio of 1x or greater, a Boid will have to check at most 27 neighboring cells for other Boids (usually less, though). At a ratio of < 1x, a Boid has to check more and more cells (for instance, 0.5x ratio -> up to 125 cell checks).
+
+On the one hand, more cell checks means more work. On the other hand, smaller cells means fewer boids per cell (and thus fewer boid checks). *Too many* cells to check is not good for performance, but *too few, big* cells to check is also bad (in the limit, it approaches one big cell, i.e. no grid at all!). The balance point, intuitively, lies around where the Boid search range equals the cell width. At this point, we don't have to perform too many extraneous cell checks, and we also don't have to perform too many extraneous Boid checks.
+
+This is confirmed by the above plot, where the optimal cell width is shown to be around 1-2x the maximum Boid search distance.
\ No newline at end of file
diff --git a/images/BoidsDemo.gif b/images/BoidsDemo.gif
new file mode 100644
index 0000000..70b8606
Binary files /dev/null and b/images/BoidsDemo.gif differ
diff --git a/images/BoidsDemo2.gif b/images/BoidsDemo2.gif
new file mode 100644
index 0000000..0f6d469
Binary files /dev/null and b/images/BoidsDemo2.gif differ
diff --git a/images/cellWidth.svg b/images/cellWidth.svg
new file mode 100644
index 0000000..4b0e404
--- /dev/null
+++ b/images/cellWidth.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/images/numBoids.svg b/images/numBoids.svg
new file mode 100644
index 0000000..3e320d3
--- /dev/null
+++ b/images/numBoids.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/images/threadsPerBlock.svg b/images/threadsPerBlock.svg
new file mode 100644
index 0000000..3b7a2af
--- /dev/null
+++ b/images/threadsPerBlock.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/src/kernel.cu b/src/kernel.cu
index 74dffcb..7534609 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -1,10 +1,12 @@
#define GLM_FORCE_CUDA
+#include
#include
#include
#include
#include
#include "utilityCore.hpp"
#include "kernel.h"
+#include
// LOOK-2.1 potentially useful for doing grid-based neighbor search
#ifndef imax
@@ -54,6 +56,10 @@ void checkCUDAError(const char *msg, int line = -1) {
/*! Size of the starting area in simulation space. */
#define scene_scale 100.0f
+/*! Defines how big a cell is compared to the maximum rule distance */
+#define cellWidthToSearchDistanceRatio 2.0f
+#define inverseCellWidthToSearchDistanceRatio (1.0f / cellWidthToSearchDistanceRatio)
+
/***********************************************
* Kernel state (pointers are device pointers) *
***********************************************/
@@ -67,21 +73,27 @@ dim3 threadsPerBlock(blockSize);
// boid cares about its neighbors' velocities.
// These are called ping-pong buffers.
glm::vec3 *dev_pos;
+glm::vec3 *dev_pos_sorted;
glm::vec3 *dev_vel1;
glm::vec3 *dev_vel2;
+glm::vec3* dev_vel1_sorted;
// 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; // This array answers the question: What index in dev_pos and dev_velX represents this particle?
+int *dev_particleGridIndices; // This array answers the question: What grid cell is this particle in?
// needed for use with thrust
thrust::device_ptr dev_thrust_particleArrayIndices;
thrust::device_ptr dev_thrust_particleGridIndices;
+thrust::device_ptr dev_thrust_pos;
+thrust::device_ptr dev_thrust_vel;
+thrust::device_ptr dev_thrust_pos_sorted;
+thrust::device_ptr dev_thrust_vel1_sorted;
-int *dev_gridCellStartIndices; // What part of dev_particleArrayIndices belongs
-int *dev_gridCellEndIndices; // to this cell?
+int *dev_gridCellStartIndices; // Together, these two arrays answer the question:
+int *dev_gridCellEndIndices; // What part of dev_particleArrayIndices belongs 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.
@@ -145,19 +157,25 @@ void Boids::initSimulation(int N) {
cudaMalloc((void**)&dev_pos, N * sizeof(glm::vec3));
checkCUDAErrorWithLine("cudaMalloc dev_pos failed!");
+ cudaMalloc((void**)&dev_pos_sorted, N * sizeof(glm::vec3));
+ checkCUDAErrorWithLine("cudaMalloc dev_pos_sorted 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!");
+ cudaMalloc((void**)&dev_vel1_sorted, N * sizeof(glm::vec3));
+ checkCUDAErrorWithLine("cudaMalloc dev_vel1_sorted 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);
+ gridCellWidth = cellWidthToSearchDistanceRatio * std::max(std::max(rule1Distance, rule2Distance), rule3Distance);
int halfSideCount = (int)(scene_scale / gridCellWidth) + 1;
gridSideCount = 2 * halfSideCount;
@@ -168,7 +186,25 @@ void Boids::initSimulation(int N) {
gridMinimum.y -= halfGridWidth;
gridMinimum.z -= halfGridWidth;
- // TODO-2.1 TODO-2.3 - Allocate additional buffers here.
+ cudaMalloc((void**)&dev_particleArrayIndices, N * sizeof(int));
+ checkCUDAErrorWithLine("cudaMalloc dev_particleArrayIndices failed!");
+
+ cudaMalloc((void**)&dev_particleGridIndices, N * sizeof(int));
+ checkCUDAErrorWithLine("cudaMalloc dev_particleGridIndicies 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!");
+
+ dev_thrust_particleArrayIndices = thrust::device_ptr(dev_particleArrayIndices);
+ dev_thrust_particleGridIndices = thrust::device_ptr(dev_particleGridIndices);
+ dev_thrust_pos = thrust::device_ptr(dev_pos);
+ dev_thrust_vel = thrust::device_ptr(dev_vel1);
+ dev_thrust_pos_sorted = thrust::device_ptr(dev_pos_sorted);
+ dev_thrust_vel1_sorted = thrust::device_ptr(dev_vel1_sorted);
+
cudaDeviceSynchronize();
}
@@ -223,6 +259,72 @@ void Boids::copyBoidsToVBO(float *vbodptr_positions, float *vbodptr_velocities)
* stepSimulation *
******************/
+// Rule 1: boids fly towards their local perceived center of mass, which excludes themselves
+__device__ glm::vec3 computeCenterOfMassVelChange(int N, int iSelf, const glm::vec3* pos) {
+ glm::vec3 perceived_center(0, 0, 0);
+ const glm::vec3 pos_self = pos[iSelf];
+ int neighborCount = 0;
+
+ for (int i = 0; i < N; ++i) {
+ if (i == iSelf) continue;
+
+ const glm::vec3& pos_i = *(pos + i);
+ if (glm::length(pos_i - pos_self) > rule1Distance) continue;
+
+ ++neighborCount;
+ perceived_center += pos_i;
+ }
+
+ if (neighborCount > 0) {
+ perceived_center /= neighborCount;
+ return (perceived_center - pos_self) * rule1Scale;
+ }
+ return glm::vec3(0, 0, 0);
+}
+
+// Rule 2: boids try to stay a distance d away from each other
+__device__ glm::vec3 computeMaintainDistanceVelChange(int N, int iSelf, const glm::vec3* pos) {
+ glm::vec3 adjustment_velocity(0, 0, 0);
+ const glm::vec3 pos_self = pos[iSelf];
+
+ for (int i = 0; i < N; ++i) {
+ if (i == iSelf) continue;
+
+ const glm::vec3 pos_i = *(pos + i);
+ const glm::vec3 pos_diff = pos_i - pos_self;
+
+ if (glm::length(pos_diff) > rule2Distance) continue;
+
+ adjustment_velocity -= pos_diff;
+ }
+
+ return adjustment_velocity * rule2Scale;
+}
+
+// Rule 3: boids try to match the speed of surrounding boids
+__device__ glm::vec3 computeVelocityMatchVelChange(int N, int iSelf, const glm::vec3* pos, const glm::vec3* vel) {
+ glm::vec3 perceived_velocity(0, 0, 0);
+ const glm::vec3 pos_self = pos[iSelf];
+ int neighborCount = 0;
+
+ for (int i = 0; i < N; ++i) {
+ if (i == iSelf) continue;
+ const glm::vec3 pos_i = *(pos + i);
+ if (glm::length(pos_i - pos_self) > rule3Distance) continue;
+
+ const glm::vec3 vel_i = *(vel + i);
+ ++neighborCount;
+ perceived_velocity += vel_i;
+ }
+
+ if (neighborCount > 0) {
+ perceived_velocity /= neighborCount;
+ return perceived_velocity * rule3Scale;
+ }
+
+ return glm::vec3(0, 0, 0);
+}
+
/**
* LOOK-1.2 You can use this as a helper for kernUpdateVelocityBruteForce.
* __device__ code can be called from a __global__ context
@@ -230,21 +332,32 @@ void Boids::copyBoidsToVBO(float *vbodptr_positions, float *vbodptr_velocities)
* in the `pos` and `vel` arrays.
*/
__device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *pos, const glm::vec3 *vel) {
- // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves
- // Rule 2: boids try to stay a distance d away from each other
- // Rule 3: boids try to match the speed of surrounding boids
- return glm::vec3(0.0f, 0.0f, 0.0f);
+ glm::vec3 deltaV1 = computeCenterOfMassVelChange(N, iSelf, pos);
+ glm::vec3 deltaV2 = computeMaintainDistanceVelChange(N, iSelf, pos);
+ glm::vec3 deltaV3 = computeVelocityMatchVelChange(N, iSelf, pos, vel);
+
+ return deltaV1 + deltaV2 + deltaV3;
}
/**
-* 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?
+ int index = threadIdx.x + (blockIdx.x * blockDim.x);
+ if (index >= N) {
+ return;
+ }
+
+ const glm::vec3 deltaV = computeVelocityChange(N, index, pos, vel1);
+ glm::vec3 updatedVelocity = vel1[index] + deltaV;
+ const float updatedSpeed = glm::length(updatedVelocity);
+ if (updatedSpeed > maxSpeed) {
+ updatedVelocity *= (maxSpeed / updatedSpeed);
+ }
+
+ // We record the new velocity into vel2 because we don't want to affect threads running in parallel that need the pre-update velocity of each boid.
+ vel2[index] = updatedVelocity;
}
/**
@@ -284,11 +397,20 @@ __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
- // - 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
+ const glm::vec3 *pos, int *indices, int *gridIndices) {
+ int index = (blockIdx.x * blockDim.x) + threadIdx.x;
+ if (index >= N) return;
+
+ indices[index] = index; // This seems useless right now, but it's necessary for sorting later since we don't
+ // have map data structures in CUDA.
+
+ // Determine what grid cell we're in.
+ const glm::vec3 pos_i = pos[index];
+ const glm::vec3 gridCell3D_i = glm::floor((pos_i - gridMin) * inverseCellWidth);
+ int gridCell1D_i = gridIndex3Dto1D(gridCell3D_i.x, gridCell3D_i.y, gridCell3D_i.z, gridResolution);
+
+ gridIndices[index] = gridCell1D_i;
+
}
// LOOK-2.1 Consider how this could be useful for indicating that a cell
@@ -300,12 +422,39 @@ __global__ void kernResetIntBuffer(int N, int *intBuffer, int 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;
+ if (index >= N) return;
+
+ int gridCell_i = particleGridIndices[index];
+
+ int previousGridCell = (index - 1 < 0) ? -1 : particleGridIndices[index - 1];
+ int nextGridCell = (index + 1 >= N) ? -1 : particleGridIndices[index + 1];
+
+ if (previousGridCell != gridCell_i) {
+ gridCellStartIndices[gridCell_i] = index;
+ }
+
+ if (nextGridCell != gridCell_i) {
+ gridCellEndIndices[gridCell_i] = index;
+ }
+}
+
+__device__ glm::vec3 wrapIndices(const glm::vec3& index, int size) {
+ return glm::vec3(
+ (static_cast(index.x) % size + size) % size,
+ (static_cast(index.y) % size + size) % size,
+ (static_cast(index.z) % size + size) % size
+ );
+}
+
+__device__ bool isBoidWithinRadiusOfGridCell(const glm::vec3& pos, const glm::vec3& gridCellMinPoint, float cellWidth, float radius) {
+ const glm::vec3 maxBounds = gridCellMinPoint + glm::vec3(cellWidth);
+ glm::vec3 closestPoint = glm::clamp(pos, gridCellMinPoint, maxBounds);
+
+ return glm::distance(pos, closestPoint) <= radius;
}
__global__ void kernUpdateVelNeighborSearchScattered(
@@ -313,15 +462,94 @@ __global__ void kernUpdateVelNeighborSearchScattered(
float inverseCellWidth, float cellWidth,
int *gridCellStartIndices, int *gridCellEndIndices,
int *particleArrayIndices,
- glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) {
+ const glm::vec3 *pos, const glm::vec3 *vel1, glm::vec3 *vel2) {
+
+ int index = (blockIdx.x * blockDim.x) + threadIdx.x;
+ if (index >= N) return;
+
// 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
+ const glm::vec3 pos_i = pos[index];
+ const glm::vec3 gridCell3D_i = glm::floor((pos_i - gridMin) * inverseCellWidth);
+
+ // Initialize tracking variables for different rules
+ glm::vec3 perceived_center(0, 0, 0);
+ glm::vec3 separation_velocity_adjustment(0, 0, 0);
+ glm::vec3 perceived_velocity(0, 0, 0);
+ glm::vec3 rule1DeltaV(0, 0, 0);
+ glm::vec3 rule2DeltaV(0, 0, 0);
+ glm::vec3 rule3DeltaV(0, 0, 0);
+ int rule1NeighborCount = 0;
+ int rule3NeighborCount = 0;
+
+ // Search over neighboring cells within the neighborhood distance of the boid
+ int searchRange = ceil(imax(inverseCellWidthToSearchDistanceRatio, 1.0f));
+ for (int x = -searchRange; x <= searchRange; ++x) {
+ for (int y = -searchRange; y <= searchRange; ++y) {
+ for (int z = -searchRange; z <= searchRange; ++z) {
+ const glm::vec3 neighborCell3DUnwrapped = gridCell3D_i + glm::vec3(x, y, z);
+ const glm::vec3 neighborCell3DWrapped = wrapIndices(neighborCell3DUnwrapped, gridResolution);
+ const glm::vec3 neighborGridCellMinPoint = gridMin + (neighborCell3DUnwrapped * cellWidth); // Note the use of the *unwrapped* neighbor here
+
+ // We can quickly eliminate some neighboring cells whose nearest point is not in range
+ if (!isBoidWithinRadiusOfGridCell(pos_i, neighborGridCellMinPoint, cellWidth, inverseCellWidthToSearchDistanceRatio * cellWidth)) continue;
+
+ // - For each cell, read the start/end indices in the boid pointer array.
+ int neighborCell1D = gridIndex3Dto1D(neighborCell3DWrapped.x, neighborCell3DWrapped.y, neighborCell3DWrapped.z, gridResolution);
+ int cellStart = gridCellStartIndices[neighborCell1D];
+ int cellEnd = gridCellEndIndices[neighborCell1D];
+ if (cellStart == -1 || cellEnd == -1) continue;
+
+ // - Collect pos and vel information from boids in neighboring cells and apply rules
+ for (int i = cellStart; i <= cellEnd; ++i) {
+ int neighborIdx = particleArrayIndices[i];
+ if (neighborIdx == index) continue;
+
+ const glm::vec3 neighborPos = pos[neighborIdx];
+ const glm::vec3 neighborVel = vel1[neighborIdx];
+ float neighborDist = glm::distance(neighborPos, pos_i);
+
+ if (neighborDist < rule1Distance) {
+ perceived_center += neighborPos;
+ ++rule1NeighborCount;
+ }
+
+ if (neighborDist < rule2Distance) {
+ rule2DeltaV -= (neighborPos - pos_i);
+ }
+
+ if (neighborDist < rule3Distance) {
+ perceived_velocity += neighborVel;
+ ++rule3NeighborCount;
+ }
+ }
+
+ }
+ }
+ }
+
+ if (rule1NeighborCount > 0) {
+ perceived_center /= rule1NeighborCount;
+ rule1DeltaV = (perceived_center - pos_i) * rule1Scale;
+ }
+
+ rule2DeltaV *= rule2Scale;
+
+ if (rule3NeighborCount > 0) {
+ perceived_velocity /= rule3NeighborCount;
+ rule3DeltaV = (perceived_velocity * rule3Scale);
+ }
+
+ const glm::vec3 deltaV = rule1DeltaV + rule2DeltaV + rule3DeltaV;
+ glm::vec3 updatedVelocity = vel1[index] + deltaV;
+ const float updatedSpeed = glm::length(updatedVelocity);
+ if (updatedSpeed > maxSpeed) { // - Clamp the speed change before putting the new speed in vel2
+ updatedVelocity *= (maxSpeed / updatedSpeed);
+ }
+
+ // We record the new velocity into vel2 because we don't want to affect threads running in parallel that need the pre-update velocity of each boid.
+ vel2[index] = updatedVelocity;
}
__global__ void kernUpdateVelNeighborSearchCoherent(
@@ -329,67 +557,210 @@ __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,
- // except with one less level of indirection.
- // This should expect gridCellStartIndices and gridCellEndIndices to refer
- // directly to pos and vel1.
- // - Identify the grid cell that this particle is in
- // - Identify which cells may contain neighbors. This isn't always 8.
- // - For each cell, read the start/end indices in the boid pointer array.
- // DIFFERENCE: For best results, consider what order the cells should be
- // checked in to maximize the memory benefits of reordering the boids data.
- // - Access each boid in the cell and compute velocity change from
- // the boids rules, if this boid is within the neighborhood distance.
- // - Clamp the speed change before putting the new speed in vel2
+ int index = (blockIdx.x * blockDim.x) + threadIdx.x;
+ if (index >= N) return;
+
+ // 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
+ const glm::vec3 pos_i = pos[index];
+ const glm::vec3 gridCell3D_i = glm::floor((pos_i - gridMin) * inverseCellWidth);
+
+ // Initialize tracking variables for different rules
+ glm::vec3 perceived_center(0, 0, 0);
+ glm::vec3 separation_velocity_adjustment(0, 0, 0);
+ glm::vec3 perceived_velocity(0, 0, 0);
+ glm::vec3 rule1DeltaV(0, 0, 0);
+ glm::vec3 rule2DeltaV(0, 0, 0);
+ glm::vec3 rule3DeltaV(0, 0, 0);
+ int rule1NeighborCount = 0;
+ int rule3NeighborCount = 0;
+
+ // Search over neighboring cells within the neighborhood distance of the boid
+ // Note that the outermost loop is over the z-direction. This aligns our computation pattern with our boid layout in data, for optimal coherency.
+ int searchRange = ceil(imax(inverseCellWidthToSearchDistanceRatio, 1.0f));
+ for (int z = -searchRange; z <= searchRange; ++z) {
+ for (int y = -searchRange; y <= searchRange; ++y) {
+ for (int x = -searchRange; x <= searchRange; ++x) {
+ const glm::vec3 neighborCell3DUnwrapped = gridCell3D_i + glm::vec3(x, y, z);
+ const glm::vec3 neighborCell3DWrapped = wrapIndices(neighborCell3DUnwrapped, gridResolution);
+ const glm::vec3 neighborGridCellMinPoint = gridMin + (neighborCell3DUnwrapped * cellWidth); // Note the use of the *unwrapped* neighbor here
+
+ // We can quickly eliminate some neighboring cells whose nearest point is not in range
+ if (!isBoidWithinRadiusOfGridCell(pos_i, neighborGridCellMinPoint, cellWidth, inverseCellWidthToSearchDistanceRatio * cellWidth)) continue;
+
+ // - For each cell, read the start/end indices in the boid pointer array.
+ int neighborCell1D = gridIndex3Dto1D(neighborCell3DWrapped.x, neighborCell3DWrapped.y, neighborCell3DWrapped.z, gridResolution);
+ int cellStart = gridCellStartIndices[neighborCell1D];
+ int cellEnd = gridCellEndIndices[neighborCell1D];
+ if (cellStart == -1 || cellEnd == -1) continue;
+
+ // - Collect pos and vel information from boids in neighboring cells and apply rules
+ for (int i = cellStart; i <= cellEnd; ++i) {
+ if (i == index) continue;
+
+ const glm::vec3 neighborPos = pos[i];
+ const glm::vec3 neighborVel = vel1[i];
+ float neighborDist = glm::distance(neighborPos, pos_i);
+
+ if (neighborDist < rule1Distance) {
+ perceived_center += neighborPos;
+ ++rule1NeighborCount;
+ }
+
+ if (neighborDist < rule2Distance) {
+ rule2DeltaV -= (neighborPos - pos_i);
+ }
+
+ if (neighborDist < rule3Distance) {
+ perceived_velocity += neighborVel;
+ ++rule3NeighborCount;
+ }
+ }
+
+ }
+ }
+ }
+
+ if (rule1NeighborCount > 0) {
+ perceived_center /= rule1NeighborCount;
+ rule1DeltaV = (perceived_center - pos_i) * rule1Scale;
+ }
+
+ rule2DeltaV *= rule2Scale;
+
+ if (rule3NeighborCount > 0) {
+ perceived_velocity /= rule3NeighborCount;
+ rule3DeltaV = (perceived_velocity * rule3Scale);
+ }
+
+ const glm::vec3 deltaV = rule1DeltaV + rule2DeltaV + rule3DeltaV;
+ glm::vec3 updatedVelocity = vel1[index] + deltaV;
+ const float updatedSpeed = glm::length(updatedVelocity);
+ if (updatedSpeed > maxSpeed) { // - Clamp the speed change before putting the new speed in vel2
+ updatedVelocity *= (maxSpeed / updatedSpeed);
+ }
+
+ // We record the new velocity into vel2 because we don't want to affect threads running in parallel that need the pre-update velocity of each boid.
+ vel2[index] = updatedVelocity;
}
/**
* 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
+ dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize);
+ kernUpdateVelocityBruteForce<<>>(numObjects, dev_pos, dev_vel1, dev_vel2);
+
+ cudaDeviceSynchronize();
+ checkCUDAErrorWithLine("Error during execution of velocity update kernel!");
+
+ kernUpdatePos<<>>(numObjects, dt, dev_pos, dev_vel2);
+
+ cudaDeviceSynchronize();
+ checkCUDAErrorWithLine("Error during execution of position update kernel!");
+
+ glm::vec3* tmp_dev_vel = dev_vel1;
+ dev_vel1 = dev_vel2;
+ dev_vel2 = tmp_dev_vel;
+
}
void Boids::stepSimulationScatteredGrid(float dt) {
+ dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize);
+
+ // First, reset the cell start and end indices to be populated with -1, so we know when a cell is empty.
+ dim3 gridIndicesResetBlocks((gridCellCount + blockSize - 1) / blockSize);
+ kernResetIntBuffer<<>>(gridCellCount, dev_gridCellStartIndices, -1);
+ kernResetIntBuffer<<>>(gridCellCount, dev_gridCellEndIndices, -1);
+
// 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.
+ 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
+ kernIdentifyCellStartEnd<<>>(numObjects, dev_particleGridIndices,
+ dev_gridCellStartIndices, dev_gridCellEndIndices);
+
// - Perform velocity updates using neighbor search
+ kernUpdateVelNeighborSearchScattered<<>>(numObjects, gridSideCount, gridMinimum,
+ gridInverseCellWidth, gridCellWidth,
+ dev_gridCellStartIndices, dev_gridCellEndIndices,
+ dev_particleArrayIndices, dev_pos, dev_vel1, dev_vel2);
+
// - Update positions
+ kernUpdatePos<<>>(numObjects, dt, dev_pos, dev_vel2);
+
// - Ping-pong buffers as needed
+ glm::vec3* tmp_dev_vel = dev_vel1;
+ dev_vel1 = dev_vel2;
+ dev_vel2 = tmp_dev_vel;
}
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.
+ dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize);
+
+ // First, reset the cell start and end indices to be populated with -1, so we know when a cell is empty.
+ dim3 gridIndicesResetBlocks((gridCellCount + blockSize - 1) / blockSize);
+ kernResetIntBuffer<<>> (gridCellCount, dev_gridCellStartIndices, -1);
+ kernResetIntBuffer<<>> (gridCellCount, dev_gridCellEndIndices, -1);
+
+ 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);
+ thrust::gather(dev_thrust_particleArrayIndices, dev_thrust_particleArrayIndices + numObjects, dev_thrust_pos, dev_thrust_pos_sorted);
+ thrust::gather(dev_thrust_particleArrayIndices, dev_thrust_particleArrayIndices + numObjects, dev_thrust_vel, dev_thrust_vel1_sorted);
+
+ // - Naively unroll the loop for finding the start and end indices of each
+ // cell's data pointers in the array of boid indices
+ kernIdentifyCellStartEnd<<>>(numObjects, dev_particleGridIndices,
+ dev_gridCellStartIndices, dev_gridCellEndIndices);
+
+ // - Perform velocity updates using neighbor search
+ 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);
+
+ // - Ping-pong buffers as needed
+ glm::vec3* tmp_dev_vel = dev_vel1;
+ dev_vel1 = dev_vel2;
+ dev_vel2 = tmp_dev_vel;
+
+ glm::vec3* tmp_dev_pos = dev_pos;
+ dev_pos = dev_pos_sorted;
+ dev_pos_sorted = tmp_dev_pos;
+
+ // And update thrust pointers
+ dev_thrust_pos = thrust::device_ptr(dev_pos);
+ dev_thrust_pos_sorted = thrust::device_ptr(dev_pos_sorted);
+ dev_thrust_vel = thrust::device_ptr(dev_vel1);
}
void Boids::endSimulation() {
cudaFree(dev_vel1);
cudaFree(dev_vel2);
+ cudaFree(dev_vel1_sorted);
cudaFree(dev_pos);
+ cudaFree(dev_pos_sorted);
- // TODO-2.1 TODO-2.3 - Free any additional buffers here.
+ cudaFree(dev_particleArrayIndices);
+ cudaFree(dev_particleGridIndices);
+ cudaFree(dev_gridCellEndIndices);
+ cudaFree(dev_gridCellStartIndices);
}
void Boids::unitTest() {
diff --git a/src/main.cpp b/src/main.cpp
index fe657ed..bc00324 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -17,8 +17,8 @@
// LOOK-2.1 LOOK-2.3 - toggles for UNIFORM_GRID and COHERENT_GRID
#define VISUALIZE 1
-#define UNIFORM_GRID 0
-#define COHERENT_GRID 0
+#define UNIFORM_GRID 1
+#define COHERENT_GRID 1
// LOOK-1.2 - change this to adjust particle count in the simulation
const int N_FOR_VIS = 5000;
+
+
+