diff --git a/README.md b/README.md
index d63a6a1..7f936b2 100644
--- a/README.md
+++ b/README.md
@@ -1,11 +1,43 @@
**University of Pennsylvania, CIS 565: 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)
+* Jilin Liu
+ * [LinkedIn](https://www.linkedin.com/in/jilin-liu97/), [twitter](https://twitter.com/Jilin18043110).
+* Tested on: Windows 10, i7-8750H @ 2.20GHz, 16GB, GTX 1050Ti 4096MB (personal)
-### (TODO: Your README)
+## Results
-Include screenshots, analysis, etc. (Remember, this is public, so don't put
-anything here that you don't want to share with the world.)
+Demo of flocking with 50000 boids using coherent grid acceleration. Details about the uniform grid structure is discussed in the [instruction](./INSTRUCTION.md).
+
+
+
+## Performance Analysis
+
+### Number of boids
+**Q1:** For each implementation, increasing the number of boids will cause lower FPS, which is expected. The reason is that when the number of boids is greater than the number of threads that GPU can schedule to run in parallel, those threads will need additional cycles to get executed. And this additional cycle number is proportioanl to the number of boids so we can see a decrease in FPS as we increase the number of boids.
+
+
+
+When turning on the visualization, the performance of each implementation suffers since we are copying additional data to VBO to be used for OpenGL.
+
+
+
+### Block size
+**Q2:** Changing the block size does not have a noticable impact on the performance. My hppothesis is that CUDA will rearrange the threads in blocks so that in each cycle there will still be as many threads as there could to run in parallel. Even though such rearrangement has overhead, the performance is dominated by other factors like the number of boids.
+
+
+
+### Coherent uniform grid
+**Q3:** I do get a performance improvement with the more coherent uniform grid. The improvement is even greater with a larger number of boids. This outcome is expected because when the number of boids is greater, the particle velocity and position data are highly distributed and thus each boid needs to access a large portion of memory, increasing the count of cache miss. When we are reshuffling the particle data array, we also need to access different part of the memory but this happens only once.
+
+
+
+### Checking 27 vs 8 neighboring cells
+**Q4:** Checking 27 neighboring cells in my simulation is actually faster than checking 8 cells. This difference is more obvious when number of boids is larger. Suppose d is the maximum search distance of the 3 rules. When we check 27 neighbors, the side of each cell is d, so we are checking a space of volume (3d)3. On the other hand, when we check 8 neighbors, the side of each cell is 2d, so we are checking a space of volume (2\*2d)3. When there are a lot more boids in the scene, the data is dense in space and checking a larger volume means we have to discard more candidates, which causes worse performance.
+
+
+
+### Time measurement
+A general time measurement for these three methods is shown below.
+
+
diff --git a/images/blockSize.JPG b/images/blockSize.JPG
new file mode 100644
index 0000000..0a60fb1
Binary files /dev/null and b/images/blockSize.JPG differ
diff --git a/images/boidsCount.JPG b/images/boidsCount.JPG
new file mode 100644
index 0000000..399d13d
Binary files /dev/null and b/images/boidsCount.JPG differ
diff --git a/images/boidsCountVisualization.JPG b/images/boidsCountVisualization.JPG
new file mode 100644
index 0000000..82f12f2
Binary files /dev/null and b/images/boidsCountVisualization.JPG differ
diff --git a/images/checkingCell.JPG b/images/checkingCell.JPG
new file mode 100644
index 0000000..af5290d
Binary files /dev/null and b/images/checkingCell.JPG differ
diff --git a/images/demo.gif b/images/demo.gif
new file mode 100644
index 0000000..4c34bf6
Binary files /dev/null and b/images/demo.gif differ
diff --git a/images/time.JPG b/images/time.JPG
new file mode 100644
index 0000000..9e41546
Binary files /dev/null and b/images/time.JPG differ
diff --git a/images/uniformGrid.JPG b/images/uniformGrid.JPG
new file mode 100644
index 0000000..c936943
Binary files /dev/null and b/images/uniformGrid.JPG differ
diff --git a/src/kernel.cu b/src/kernel.cu
index 74dffcb..dcb56fa 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -86,6 +86,8 @@ 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_pos2;
+
// LOOK-2.1 - Grid parameters based on simulation parameters.
// These are automatically computed for you in Boids::initSimulation
int gridCellCount;
@@ -137,39 +139,61 @@ __global__ void kernGenerateRandomPosArray(int time, int N, glm::vec3 * arr, flo
* Initialize memory, update some globals
*/
void Boids::initSimulation(int N) {
- numObjects = N;
- dim3 fullBlocksPerGrid((N + blockSize - 1) / blockSize);
+ 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!");
+ // 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_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_vel2, N * sizeof(glm::vec3));
+ checkCUDAErrorWithLine("cudaMalloc dev_vel2 failed!");
- // LOOK-1.2 - This is a typical CUDA kernel invocation.
- kernGenerateRandomPosArray<<>>(1, numObjects,
+ // LOOK-1.2 - This is a typical CUDA kernel invocation.
+ kernGenerateRandomPosArray<<>>(1, numObjects,
dev_pos, scene_scale);
- checkCUDAErrorWithLine("kernGenerateRandomPosArray failed!");
+ 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;
+ // LOOK-2.1 computing grid params
- gridCellCount = gridSideCount * gridSideCount * gridSideCount;
- gridInverseCellWidth = 1.0f / gridCellWidth;
- float halfGridWidth = gridCellWidth * halfSideCount;
- gridMinimum.x -= halfGridWidth;
- gridMinimum.y -= halfGridWidth;
- gridMinimum.z -= halfGridWidth;
+#define testing27
+#ifdef testing27
+ gridCellWidth = std::max(std::max(rule1Distance, rule2Distance), rule3Distance);
+#else
+ gridCellWidth = 2.0f * std::max(std::max(rule1Distance, rule2Distance), rule3Distance);
+#endif
+
+ int halfSideCount = (int)(scene_scale / gridCellWidth) + 1;
+ gridSideCount = 2 * halfSideCount;
- // TODO-2.1 TODO-2.3 - Allocate additional buffers here.
- cudaDeviceSynchronize();
+ gridCellCount = gridSideCount * gridSideCount * gridSideCount;
+ gridInverseCellWidth = 1.0f / gridCellWidth;
+ float halfGridWidth = gridCellWidth * halfSideCount;
+ gridMinimum.x -= halfGridWidth;
+ gridMinimum.y -= halfGridWidth;
+ gridMinimum.z -= halfGridWidth;
+
+ // TODO-2.1 TODO-2.3 - Allocate additional buffers here.
+ cudaMalloc((void**)&dev_particleArrayIndices, N * sizeof(int));
+ checkCUDAErrorWithLine("cudaMalloc dev_particleArrayIndices failed!");
+
+ cudaMalloc((void**)&dev_particleGridIndices, N * sizeof(int));
+ checkCUDAErrorWithLine("cudaMalloc dev_particleGridIndices failed!");
+
+ cudaMalloc((void**)&dev_gridCellStartIndices, gridCellCount * sizeof(int));
+ checkCUDAErrorWithLine("cudaMalloc dev_gridCellStartIndices failed!");
+
+ cudaMalloc((void**)&dev_gridCellEndIndices, gridCellCount * sizeof(int));
+ checkCUDAErrorWithLine("cudaMalloc dev_gridCellEndIndices failed!");
+
+ cudaMalloc((void**)&dev_pos2, N * sizeof(glm::vec3));
+ checkCUDAErrorWithLine("cudaMalloc dev_pos2 failed!");
+
+ cudaDeviceSynchronize();
}
@@ -230,10 +254,51 @@ 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);
+ // 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
+
+
+ // Rule 1: Boids try to fly towards the center of mass of neighboring boids.
+ glm::vec3 perceivedCenter(0.f);
+ int numberOfNeighbors = 0;
+ for (int i = 0; i < N; i++) {
+ if (i != iSelf && glm::distance(pos[iSelf], pos[i]) < rule1Distance) {
+ perceivedCenter += pos[i];
+ numberOfNeighbors++;
+ }
+ }
+ glm::vec3 v1(0.f);
+ if (numberOfNeighbors > 0) {
+ perceivedCenter /= numberOfNeighbors;
+ v1 = (perceivedCenter - pos[iSelf]) * rule1Scale;
+ }
+
+ // Rule 2: Boids try to keep a small distance away from other objects.
+ glm::vec3 c(0.f);
+ for (int i = 0; i < N; i++) {
+ if (i != iSelf && glm::distance(pos[iSelf], pos[i]) < rule2Distance) {
+ c -= pos[i] - pos[iSelf];
+ }
+ }
+ glm::vec3 v2 = c * rule2Scale;
+
+ // Rule 3: Boids try to match velocity with near boids.
+ glm::vec3 perceivedVelocity(0.f);
+ numberOfNeighbors = 0;
+ for (int i = 0; i < N; i++) {
+ if (i != iSelf && glm::distance(pos[iSelf], pos[i]) < rule3Distance) {
+ perceivedVelocity += vel[i];
+ numberOfNeighbors++;
+ }
+ }
+ glm::vec3 v3(0.f);
+ if (numberOfNeighbors > 0) {
+ perceivedVelocity /= numberOfNeighbors;
+ v3 = perceivedVelocity * rule3Scale;
+ }
+
+ return v1 + v2 + v3;
}
/**
@@ -241,10 +306,24 @@ __device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *po
* 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?
+ 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;
+ }
+
+ glm::vec3 newVel = vel1[index] + computeVelocityChange(N, index, pos, vel1);
+ if (glm::length(newVel) > maxSpeed) {
+ newVel *= maxSpeed / glm::length(newVel);
+ }
+
+ // Because other threads will potentially access velocity of this boid,
+ // we cannot directly update vel1. Instead, we store the value in vel2.
+ vel2[index] = newVel;
}
/**
@@ -252,24 +331,24 @@ __global__ void kernUpdateVelocityBruteForce(int N, glm::vec3 *pos,
* 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;
+ // Update position by velocity
+ int index = threadIdx.x + (blockIdx.x * blockDim.x);
+ if (index >= N) {
+ return;
+ }
+ glm::vec3 thisPos = pos[index];
+ thisPos += vel[index] * dt;
- // Wrap the boids around so we don't lose them
- thisPos.x = thisPos.x < -scene_scale ? scene_scale : thisPos.x;
- thisPos.y = thisPos.y < -scene_scale ? scene_scale : thisPos.y;
- thisPos.z = thisPos.z < -scene_scale ? scene_scale : thisPos.z;
+ // Wrap the boids around so we don't lose them
+ thisPos.x = thisPos.x < -scene_scale ? scene_scale : thisPos.x;
+ thisPos.y = thisPos.y < -scene_scale ? scene_scale : thisPos.y;
+ thisPos.z = thisPos.z < -scene_scale ? scene_scale : thisPos.z;
- thisPos.x = thisPos.x > scene_scale ? -scene_scale : thisPos.x;
- thisPos.y = thisPos.y > scene_scale ? -scene_scale : thisPos.y;
- thisPos.z = thisPos.z > scene_scale ? -scene_scale : thisPos.z;
+ thisPos.x = thisPos.x > scene_scale ? -scene_scale : thisPos.x;
+ thisPos.y = thisPos.y > scene_scale ? -scene_scale : thisPos.y;
+ thisPos.z = thisPos.z > scene_scale ? -scene_scale : thisPos.z;
- pos[index] = thisPos;
+ pos[index] = thisPos;
}
// LOOK-2.1 Consider this method of computing a 1D index from a 3D grid index.
@@ -279,179 +358,478 @@ __global__ void kernUpdatePos(int N, float dt, glm::vec3 *pos, glm::vec3 *vel) {
// for(y)
// for(z)? Or some other order?
__device__ int gridIndex3Dto1D(int x, int y, int z, int gridResolution) {
- return x + y * gridResolution + z * gridResolution * gridResolution;
+ return x + y * gridResolution + z * gridResolution * gridResolution;
}
__global__ void kernComputeIndices(int N, int gridResolution,
- glm::vec3 gridMin, float inverseCellWidth,
- glm::vec3 *pos, int *indices, int *gridIndices) {
+ glm::vec3 gridMin, float inverseCellWidth,
+ glm::vec3 *pos, int *indices, int *gridIndices) {
// TODO-2.1
// - Label each boid with the index of its grid cell.
// - Set up a parallel array of integer indices as pointers to the actual
// boid data in pos and vel1/vel2
+
+ int index = threadIdx.x + blockIdx.x * blockDim.x;
+ if (index >= N) {
+ return;
+ }
+ glm::vec3 position = pos[index];
+ glm::vec3 offset = (position - gridMin) * inverseCellWidth;
+ glm::vec3 floor = glm::floor(offset);
+ int gi = gridIndex3Dto1D((int)floor.x, (int)floor.y, (int)floor.z,
+ gridResolution);
+ gridIndices[index] = gi;
+ indices[index] = index;
}
// 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;
- }
+ 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 *gridCellStartIndices, int *gridCellEndIndices) {
+ // TODO-2.1
+ // Identify the start point of each cell in the gridIndices array.
+ // This is basically a parallel unrolling of a loop that goes
+ // "this index doesn't match the one before it, must be a new cell!"
+
+ int index = threadIdx.x + blockIdx.x * blockDim.x;
+ if (index >= N) {
+ return;
+ }
+ else if (index == 0) {
+ int curCell = particleGridIndices[index];
+ gridCellStartIndices[curCell] = 0;
+ }
+ else {
+ int curCell = particleGridIndices[index];
+ int prevCell = particleGridIndices[index - 1];
+ if (curCell != prevCell) {
+ gridCellStartIndices[curCell] = index;
+ gridCellEndIndices[prevCell] = index;
+ }
+ if (index == N - 1) {
+ gridCellEndIndices[curCell] = N;
+ }
+ }
}
__global__ void kernUpdateVelNeighborSearchScattered(
- int N, int gridResolution, glm::vec3 gridMin,
- float inverseCellWidth, float cellWidth,
- int *gridCellStartIndices, int *gridCellEndIndices,
- int *particleArrayIndices,
- glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) {
- // TODO-2.1 - Update a boid's velocity using the uniform grid to reduce
- // the number of boids that need to be checked.
- // - Identify the grid cell that this particle is in
- // - Identify which cells may contain neighbors. This isn't always 8.
- // - For each cell, read the start/end indices in the boid pointer array.
- // - Access each boid in the cell and compute velocity change from
- // the boids rules, if this boid is within the neighborhood distance.
- // - Clamp the speed change before putting the new speed in vel2
+ int N, int gridResolution, glm::vec3 gridMin,
+ float inverseCellWidth, float cellWidth,
+ int *gridCellStartIndices, int *gridCellEndIndices,
+ int *particleArrayIndices,
+ glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) {
+ // TODO-2.1 - Update a boid's velocity using the uniform grid to reduce
+ // the number of boids that need to be checked.
+ // - Identify the grid cell that this particle is in
+ // - Identify which cells may contain neighbors. This isn't always 8.
+ // - For each cell, read the start/end indices in the boid pointer array.
+ // - Access each boid in the cell and compute velocity change from
+ // the boids rules, if this boid is within the neighborhood distance.
+ // - Clamp the speed change before putting the new speed in vel2
+
+ int index = threadIdx.x + blockIdx.x * blockDim.x;
+ if (index >= N) {
+ return;
+ }
+ int boid = particleArrayIndices[index];
+ glm::vec3 position = pos[boid];
+ glm::vec3 offset = (position - gridMin) * inverseCellWidth;
+ glm::vec3 floor = glm::floor(offset);
+
+#ifdef testing27
+ int sz = 3;
+ int3 start = make_int3((int)floor.x - 1, (int)floor.y - 1, (int)floor.z - 1);
+ int3 next = make_int3(1, 1, 1);
+#else
+ int sz = 2;
+ int3 start = make_int3((int)floor.x, (int)floor.y, (int)floor.z);
+ int3 next = make_int3(
+ (offset.x > floor.x + 0.5f) ? 1 : -1,
+ (offset.y > floor.y + 0.5f) ? 1 : -1,
+ (offset.z > floor.z + 0.5f) ? 1 : -1
+ );
+#endif
+ glm::vec3 v1, v2, v3;
+ int numberOfNeighbors1 = 0;
+ int numberOfNeighbors3 = 0;
+ glm::vec3 perceivedCenter(0.f);
+ glm::vec3 c(0.f);
+ glm::vec3 perceivedVelocity(0.f);
+ for (int k = 0; k < sz; k++) {
+ for (int j = 0; j < sz; j++) {
+ for (int i = 0; i < sz; i++) {
+ int3 cur = make_int3(
+ start.x + i * next.x,
+ start.y + j * next.y,
+ start.z + k * next.z
+ );
+ if (0 <= cur.x && cur.x < gridResolution &&
+ 0 <= cur.y && cur.y < gridResolution &&
+ 0 <= cur.z && cur.z < gridResolution) {
+ int gi = gridIndex3Dto1D(cur.x, cur.y, cur.z, gridResolution);
+ for (int nIdx = gridCellStartIndices[gi];
+ nIdx < gridCellEndIndices[gi];
+ nIdx++) {
+ int neighbor = particleArrayIndices[nIdx];
+ if (boid != neighbor) {
+ float dist = glm::distance(pos[neighbor], pos[boid]);
+ // Rule 1
+ if (dist < rule1Distance) {
+ perceivedCenter += pos[neighbor];
+ numberOfNeighbors1++;
+ }
+ // Rule 2
+ if (dist < rule2Distance) {
+ c -= (pos[neighbor] - pos[boid]);
+ }
+ // Rule 3
+ if (dist < rule3Distance) {
+ perceivedVelocity += vel1[neighbor];
+ numberOfNeighbors3++;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // Continue Rule 1
+ if (numberOfNeighbors1 > 0) {
+ perceivedCenter /= numberOfNeighbors1;
+ v1 = (perceivedCenter - pos[boid]) * rule1Scale;
+ }
+
+ // Continue Rule 2
+ v2 = c * rule2Scale;
+
+ // Continue Rule 3
+ if (numberOfNeighbors3 > 0) {
+ perceivedVelocity /= numberOfNeighbors3;
+ v3 = perceivedVelocity * rule3Scale;
+ }
+
+ glm::vec3 newVel = vel1[boid] + v1 + v2 + v3;
+ if (glm::length(newVel) > maxSpeed) {
+ newVel *= maxSpeed / glm::length(newVel);
+ }
+ vel2[boid] = newVel;
}
__global__ void kernUpdateVelNeighborSearchCoherent(
- int N, int gridResolution, glm::vec3 gridMin,
- float inverseCellWidth, float cellWidth,
- int *gridCellStartIndices, int *gridCellEndIndices,
- glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) {
- // TODO-2.3 - This should be very similar to kernUpdateVelNeighborSearchScattered,
- // except with one less level of indirection.
- // This should expect gridCellStartIndices and gridCellEndIndices to refer
- // directly to pos and vel1.
- // - Identify the grid cell that this particle is in
- // - Identify which cells may contain neighbors. This isn't always 8.
- // - For each cell, read the start/end indices in the boid pointer array.
- // DIFFERENCE: For best results, consider what order the cells should be
- // checked in to maximize the memory benefits of reordering the boids data.
- // - Access each boid in the cell and compute velocity change from
- // the boids rules, if this boid is within the neighborhood distance.
- // - Clamp the speed change before putting the new speed in vel2
+ int N, int gridResolution, glm::vec3 gridMin,
+ float inverseCellWidth, float cellWidth,
+ int *gridCellStartIndices, int *gridCellEndIndices,
+ glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) {
+ // TODO-2.3 - This should be very similar to kernUpdateVelNeighborSearchScattered,
+ // except with one less level of indirection.
+ // This should expect gridCellStartIndices and gridCellEndIndices to refer
+ // directly to pos and vel1.
+ // - Identify the grid cell that this particle is in
+ // - Identify which cells may contain neighbors. This isn't always 8.
+ // - For each cell, read the start/end indices in the boid pointer array.
+ // DIFFERENCE: For best results, consider what order the cells should be
+ // checked in to maximize the memory benefits of reordering the boids data.
+ // - Access each boid in the cell and compute velocity change from
+ // the boids rules, if this boid is within the neighborhood distance.
+ // - Clamp the speed change before putting the new speed in vel2
+
+ int index = threadIdx.x + blockIdx.x * blockDim.x;
+ if (index >= N) {
+ return;
+ }
+ glm::vec3 offset = (pos[index] - gridMin) * inverseCellWidth;
+ glm::vec3 floor = glm::floor(offset);
+
+#ifdef testing27
+ int sz = 3;
+ int3 start = make_int3((int)floor.x - 1, (int)floor.y - 1, (int)floor.z - 1);
+ int3 next = make_int3(1, 1, 1);
+#else
+ int sz = 2;
+ int3 start = make_int3((int)floor.x, (int)floor.y, (int)floor.z);
+ int3 next = make_int3(
+ (offset.x > floor.x + 0.5f) ? 1 : -1,
+ (offset.y > floor.y + 0.5f) ? 1 : -1,
+ (offset.z > floor.z + 0.5f) ? 1 : -1
+ );
+#endif
+
+ glm::vec3 v1, v2, v3;
+ int numberOfNeighbors1 = 0;
+ int numberOfNeighbors3 = 0;
+ glm::vec3 perceivedCenter(0.f);
+ glm::vec3 c(0.f);
+ glm::vec3 perceivedVelocity(0.f);
+ for (int k = 0; k < sz; k++) {
+ for (int j = 0; j < sz; j++) {
+ for (int i = 0; i < sz; i++) {
+ int3 cur = make_int3(
+ start.x + i * next.x,
+ start.y + j * next.y,
+ start.z + k * next.z
+ );
+ if (0 <= cur.x && cur.x < gridResolution &&
+ 0 <= cur.y && cur.y < gridResolution &&
+ 0 <= cur.z && cur.z < gridResolution) {
+ int gi = gridIndex3Dto1D(cur.x, cur.y, cur.z, gridResolution);
+ for (int nIdx = gridCellStartIndices[gi];
+ nIdx < gridCellEndIndices[gi];
+ nIdx++) {
+ if (index != nIdx) {
+ float dist = glm::distance(pos[nIdx], pos[index]);
+ // Rule 1
+ if (dist < rule1Distance) {
+ perceivedCenter += pos[nIdx];
+ numberOfNeighbors1++;
+ }
+ // Rule 2
+ if (dist < rule2Distance) {
+ c -= (pos[nIdx] - pos[index]);
+ }
+ // Rule 3
+ if (dist < rule3Distance) {
+ perceivedVelocity += vel1[nIdx];
+ numberOfNeighbors3++;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // Continue Rule 1
+ if (numberOfNeighbors1 > 0) {
+ perceivedCenter /= numberOfNeighbors1;
+ v1 = (perceivedCenter - pos[index]) * rule1Scale;
+ }
+
+ // Continue Rule 2
+ v2 = c * rule2Scale;
+
+ // Continue Rule 3
+ if (numberOfNeighbors3 > 0) {
+ perceivedVelocity /= numberOfNeighbors3;
+ v3 = perceivedVelocity * rule3Scale;
+ }
+
+ glm::vec3 newVel = vel1[index] + v1 + v2 + v3;
+ if (glm::length(newVel) > maxSpeed) {
+ newVel *= maxSpeed / glm::length(newVel);
+ }
+ vel2[index] = newVel;
+}
+
+__global__ void kernReshuffleParticle(
+ int N, int *particleArrayIndices,
+ glm::vec3* pos1, glm::vec3* pos2,
+ glm::vec3* vel1, glm::vec3* vel2) {
+ // Reshuffle the particle data
+ int index = threadIdx.x + blockIdx.x * blockDim.x;
+ if (index >= N) {
+ return;
+ }
+
+ pos2[index] = pos1[particleArrayIndices[index]];
+ vel2[index] = vel1[particleArrayIndices[index]];
}
/**
* 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
+
+ dim3 blocks(numObjects / threadsPerBlock.x + 1);
+ kernUpdateVelocityBruteForce<<>>(numObjects, dev_pos,
+ dev_vel1, dev_vel2);
+ kernUpdatePos<<>>(numObjects, dt, dev_pos, dev_vel2);
+ 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
+
+ // Be careful of different array length. Use different block numbers.
+ dim3 boidBlocks(numObjects / threadsPerBlock.x + 1);
+ dim3 cellBlocks(gridCellCount / threadsPerBlock.x + 1);
+
+ kernComputeIndices << > > (numObjects, gridSideCount,
+ gridMinimum, gridInverseCellWidth, dev_pos, dev_particleArrayIndices,
+ dev_particleGridIndices);
+
+ dev_thrust_particleArrayIndices = thrust::device_ptr(dev_particleArrayIndices);
+ dev_thrust_particleGridIndices = thrust::device_ptr(dev_particleGridIndices);
+ thrust::sort_by_key(dev_thrust_particleGridIndices,
+ dev_thrust_particleGridIndices + numObjects,
+ dev_thrust_particleArrayIndices);
+
+ kernResetIntBuffer << > > (gridCellCount,
+ dev_gridCellStartIndices, - 1);
+ kernResetIntBuffer << > > (gridCellCount,
+ dev_gridCellEndIndices, -1);
+ kernIdentifyCellStartEnd << > > (numObjects,
+ dev_particleGridIndices,
+ dev_gridCellStartIndices, dev_gridCellEndIndices);
+
+ kernUpdateVelNeighborSearchScattered << > > (
+ numObjects, gridSideCount, gridMinimum,
+ gridInverseCellWidth, gridCellWidth,
+ dev_gridCellStartIndices, dev_gridCellEndIndices,
+ dev_particleArrayIndices,
+ dev_pos, dev_vel1, dev_vel2);
+
+ kernUpdatePos << > > (numObjects, dt, dev_pos, dev_vel2);
+
+ std::swap(dev_vel1, dev_vel2);
}
void Boids::stepSimulationCoherentGrid(float dt) {
- // TODO-2.3 - start by copying Boids::stepSimulationNaiveGrid
- // Uniform Grid Neighbor search using Thrust sort on cell-coherent data.
- // In Parallel:
- // - Label each particle with its array index as well as its grid index.
- // Use 2x width grids
- // - Unstable key sort using Thrust. A stable sort isn't necessary, but you
- // are welcome to do a performance comparison.
- // - Naively unroll the loop for finding the start and end indices of each
- // cell's data pointers in the array of boid indices
- // - BIG DIFFERENCE: use the rearranged array index buffer to reshuffle all
- // the particle data in the simulation array.
- // CONSIDER WHAT ADDITIONAL BUFFERS YOU NEED
- // - Perform velocity updates using neighbor search
- // - Update positions
- // - Ping-pong buffers as needed. THIS MAY BE DIFFERENT FROM BEFORE.
+ // TODO-2.3 - start by copying Boids::stepSimulationNaiveGrid
+ // Uniform Grid Neighbor search using Thrust sort on cell-coherent data.
+ // In Parallel:
+ // - Label each particle with its array index as well as its grid index.
+ // Use 2x width grids
+ // - Unstable key sort using Thrust. A stable sort isn't necessary, but you
+ // are welcome to do a performance comparison.
+ // - Naively unroll the loop for finding the start and end indices of each
+ // cell's data pointers in the array of boid indices
+ // - BIG DIFFERENCE: use the rearranged array index buffer to reshuffle all
+ // the particle data in the simulation array.
+ // CONSIDER WHAT ADDITIONAL BUFFERS YOU NEED
+ // - Perform velocity updates using neighbor search
+ // - Update positions
+ // - Ping-pong buffers as needed. THIS MAY BE DIFFERENT FROM BEFORE.
+
+ dim3 boidBlocks(numObjects / threadsPerBlock.x + 1);
+ dim3 cellBlocks(gridCellCount / threadsPerBlock.x + 1);
+
+ kernComputeIndices << > > (numObjects, gridSideCount,
+ gridMinimum, gridInverseCellWidth, dev_pos, dev_particleArrayIndices,
+ dev_particleGridIndices);
+
+ dev_thrust_particleArrayIndices = thrust::device_ptr(dev_particleArrayIndices);
+ dev_thrust_particleGridIndices = thrust::device_ptr(dev_particleGridIndices);
+ thrust::sort_by_key(dev_thrust_particleGridIndices,
+ dev_thrust_particleGridIndices + numObjects,
+ dev_thrust_particleArrayIndices);
+
+ kernResetIntBuffer << > > (gridCellCount,
+ dev_gridCellStartIndices, -1);
+ kernResetIntBuffer << > > (gridCellCount,
+ dev_gridCellEndIndices, -1);
+ kernIdentifyCellStartEnd << > > (numObjects,
+ dev_particleGridIndices,
+ dev_gridCellStartIndices, dev_gridCellEndIndices);
+
+ kernReshuffleParticle << > > (numObjects,
+ dev_particleArrayIndices, dev_pos, dev_pos2, dev_vel1, dev_vel2);
+ std::swap(dev_pos, dev_pos2);
+ std::swap(dev_vel1, dev_vel2);
+
+ kernUpdateVelNeighborSearchCoherent << > > (
+ numObjects, gridSideCount, gridMinimum,
+ gridInverseCellWidth, gridCellWidth,
+ dev_gridCellStartIndices, dev_gridCellEndIndices,
+ dev_pos, dev_vel1, dev_vel2);
+
+ kernUpdatePos << > > (numObjects, dt, dev_pos, dev_vel2);
+ std::swap(dev_vel1, dev_vel2);
}
void Boids::endSimulation() {
- cudaFree(dev_vel1);
- cudaFree(dev_vel2);
- cudaFree(dev_pos);
-
- // TODO-2.1 TODO-2.3 - Free any additional buffers here.
+ cudaFree(dev_vel1);
+ cudaFree(dev_vel2);
+ cudaFree(dev_pos);
+
+ // TODO-2.1 TODO-2.3 - Free any additional buffers here.
+ cudaFree(dev_particleArrayIndices);
+ cudaFree(dev_particleGridIndices);
+ cudaFree(dev_gridCellStartIndices);
+ cudaFree(dev_gridCellEndIndices);
+ cudaFree(dev_pos2);
}
void Boids::unitTest() {
- // LOOK-1.2 Feel free to write additional tests here.
-
- // test unstable sort
- int *dev_intKeys;
- int *dev_intValues;
- int N = 10;
-
- std::unique_ptrintKeys{ new int[N] };
- std::unique_ptrintValues{ new int[N] };
-
- intKeys[0] = 0; intValues[0] = 0;
- intKeys[1] = 1; intValues[1] = 1;
- intKeys[2] = 0; intValues[2] = 2;
- intKeys[3] = 3; intValues[3] = 3;
- intKeys[4] = 0; intValues[4] = 4;
- intKeys[5] = 2; intValues[5] = 5;
- intKeys[6] = 2; intValues[6] = 6;
- intKeys[7] = 0; intValues[7] = 7;
- intKeys[8] = 5; intValues[8] = 8;
- intKeys[9] = 6; intValues[9] = 9;
-
- cudaMalloc((void**)&dev_intKeys, N * sizeof(int));
- checkCUDAErrorWithLine("cudaMalloc dev_intKeys failed!");
-
- cudaMalloc((void**)&dev_intValues, N * sizeof(int));
- checkCUDAErrorWithLine("cudaMalloc dev_intValues failed!");
-
- dim3 fullBlocksPerGrid((N + blockSize - 1) / blockSize);
-
- std::cout << "before unstable sort: " << std::endl;
- for (int i = 0; i < N; i++) {
- std::cout << " key: " << intKeys[i];
- std::cout << " value: " << intValues[i] << std::endl;
- }
-
- // How to copy data to the GPU
- cudaMemcpy(dev_intKeys, intKeys.get(), sizeof(int) * N, cudaMemcpyHostToDevice);
- cudaMemcpy(dev_intValues, intValues.get(), sizeof(int) * N, cudaMemcpyHostToDevice);
-
- // Wrap device vectors in thrust iterators for use with thrust.
- thrust::device_ptr dev_thrust_keys(dev_intKeys);
- thrust::device_ptr dev_thrust_values(dev_intValues);
- // LOOK-2.1 Example for using thrust::sort_by_key
- thrust::sort_by_key(dev_thrust_keys, dev_thrust_keys + N, dev_thrust_values);
-
- // How to copy data back to the CPU side from the GPU
- cudaMemcpy(intKeys.get(), dev_intKeys, sizeof(int) * N, cudaMemcpyDeviceToHost);
- cudaMemcpy(intValues.get(), dev_intValues, sizeof(int) * N, cudaMemcpyDeviceToHost);
- checkCUDAErrorWithLine("memcpy back failed!");
-
- std::cout << "after unstable sort: " << std::endl;
- for (int i = 0; i < N; i++) {
- std::cout << " key: " << intKeys[i];
- std::cout << " value: " << intValues[i] << std::endl;
- }
-
- // cleanup
- cudaFree(dev_intKeys);
- cudaFree(dev_intValues);
- checkCUDAErrorWithLine("cudaFree failed!");
- return;
+ // LOOK-1.2 Feel free to write additional tests here.
+
+ // test unstable sort
+ int *dev_intKeys;
+ int *dev_intValues;
+ int N = 10;
+
+ std::unique_ptrintKeys{ new int[N] };
+ std::unique_ptrintValues{ new int[N] };
+
+ intKeys[0] = 0; intValues[0] = 0;
+ intKeys[1] = 1; intValues[1] = 1;
+ intKeys[2] = 0; intValues[2] = 2;
+ intKeys[3] = 3; intValues[3] = 3;
+ intKeys[4] = 0; intValues[4] = 4;
+ intKeys[5] = 2; intValues[5] = 5;
+ intKeys[6] = 2; intValues[6] = 6;
+ intKeys[7] = 0; intValues[7] = 7;
+ intKeys[8] = 5; intValues[8] = 8;
+ intKeys[9] = 6; intValues[9] = 9;
+
+ cudaMalloc((void**)&dev_intKeys, N * sizeof(int));
+ checkCUDAErrorWithLine("cudaMalloc dev_intKeys failed!");
+
+ cudaMalloc((void**)&dev_intValues, N * sizeof(int));
+ checkCUDAErrorWithLine("cudaMalloc dev_intValues failed!");
+
+ dim3 fullBlocksPerGrid((N + blockSize - 1) / blockSize);
+
+ std::cout << "before unstable sort: " << std::endl;
+ for (int i = 0; i < N; i++) {
+ std::cout << " key: " << intKeys[i];
+ std::cout << " value: " << intValues[i] << std::endl;
+ }
+
+ // How to copy data to the GPU
+ cudaMemcpy(dev_intKeys, intKeys.get(), sizeof(int) * N, cudaMemcpyHostToDevice);
+ cudaMemcpy(dev_intValues, intValues.get(), sizeof(int) * N, cudaMemcpyHostToDevice);
+
+ // Wrap device vectors in thrust iterators for use with thrust.
+ thrust::device_ptr dev_thrust_keys(dev_intKeys);
+ thrust::device_ptr dev_thrust_values(dev_intValues);
+ // LOOK-2.1 Example for using thrust::sort_by_key
+ thrust::sort_by_key(dev_thrust_keys, dev_thrust_keys + N, dev_thrust_values);
+
+ // How to copy data back to the CPU side from the GPU
+ cudaMemcpy(intKeys.get(), dev_intKeys, sizeof(int) * N, cudaMemcpyDeviceToHost);
+ cudaMemcpy(intValues.get(), dev_intValues, sizeof(int) * N, cudaMemcpyDeviceToHost);
+ checkCUDAErrorWithLine("memcpy back failed!");
+
+ std::cout << "after unstable sort: " << std::endl;
+ for (int i = 0; i < N; i++) {
+ std::cout << " key: " << intKeys[i];
+ std::cout << " value: " << intValues[i] << std::endl;
+ }
+
+ // cleanup
+ cudaFree(dev_intKeys);
+ cudaFree(dev_intValues);
+ checkCUDAErrorWithLine("cudaFree failed!");
+ return;
}
diff --git a/src/main.cpp b/src/main.cpp
index b82c8c6..8ab2b9b 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -14,11 +14,11 @@
// LOOK-2.1 LOOK-2.3 - toggles for UNIFORM_GRID and COHERENT_GRID
#define VISUALIZE 1
-#define UNIFORM_GRID 0
-#define COHERENT_GRID 0
+#define UNIFORM_GRID 1
+#define COHERENT_GRID 1
// LOOK-1.2 - change this to adjust particle count in the simulation
-const int N_FOR_VIS = 5000;
+const int N_FOR_VIS = 50000;
const float DT = 0.2f;
/**
@@ -195,6 +195,12 @@ void initShaders(GLuint * program) {
cudaGLMapBufferObject((void**)&dptrVertPositions, boidVBO_positions);
cudaGLMapBufferObject((void**)&dptrVertVelocities, boidVBO_velocities);
+ // Timing using CUDA Events
+ cudaEvent_t start, stop;
+ cudaEventCreate(&start);
+ cudaEventCreate(&stop);
+ cudaEventRecord(start);
+
// execute the kernel
#if UNIFORM_GRID && COHERENT_GRID
Boids::stepSimulationCoherentGrid(DT);
@@ -204,6 +210,17 @@ void initShaders(GLuint * program) {
Boids::stepSimulationNaive(DT);
#endif
+ // Output execution time of each step
+ static float totalTime = 0;
+ static float frame = 0;
+ cudaEventRecord(stop);
+ cudaEventSynchronize(stop);
+ float milliseconds = 0;
+ cudaEventElapsedTime(&milliseconds, start, stop);
+ totalTime += milliseconds;
+ frame += 1;
+ //std::cout << (totalTime / frame) << std::endl;
+
#if VISUALIZE
Boids::copyBoidsToVBO(dptrVertPositions, dptrVertVelocities);
#endif