diff --git a/CMakeLists.txt b/CMakeLists.txt index 2f16f40..3dce6c7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -56,6 +56,7 @@ set(headers src/kernel.h src/main.hpp src/utilityCore.hpp + src/boidImpl.h ) set(sources @@ -63,6 +64,7 @@ set(sources src/kernel.cu src/main.cpp src/utilityCore.cpp + src/boidImpl.cu ) list(SORT headers) diff --git a/README.md b/README.md index ee39093..9a2ffd5 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,58 @@ **University of Pennsylvania, CIS 5650: GPU Programming and Architecture, -Project 1 - Flocking** +Project 1 - Flocking - LATE DAY USED** -* (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) +* Henry Han + * https://github.com/sirenri2001 + * https://www.linkedin.com/in/henry-han-a832a6284/ +* Tested on: Windows 11 Pro 24H2, i7-9750H @ 2.60GHz 16GB, RTX 2070 Max-Q -### (TODO: Your README) +**I modified CMakeLists.txt** -Include screenshots, analysis, etc. (Remember, this is public, so don't put -anything here that you don't want to share with the world.) +## Screenshots + +![](/profiles/11.gif) + +Tested with coherent grid algorithm with particle spawn number of 1,000,000 + +Average FPS 130 + +## Performance Analysis + +### Analysis Conclusion + +This boid algorithm works well under a certain ratio of particles per cell grid. More particle in a grid, more performance impact on particle interacts with neighbors. Also, more grid means more cost on sorting particles along with grid indices. Therefore, this algorithm can demostrate best performance when average number of particle in each grid cell stays within 0.1 ~ 10. + +### Changing of Factors + +I tested the following factors that may have an impact on performance + +- Particle Count +- Neighbor Searching Strategy (8x or 27x grid method) +- Simulation Domain (larger domain contains more grid) + +Also I tested on blockSize, but it seems no impact on performance. + +Here is a diagram of my test result. + +![](/profiles/FrameRateFactors.png) + + +### Performance Analysis + +I use following setting as baseline: 1 million particles with initial grid size and domain size, as well as using coherent method and 27x neighbor grids. This yields an FPS of 142.599, as shown in bold text in the diagram below. + +![](/profiles/ProfileBaseLine.png) + +Also, note that yellow line shows when there are fewer grid cells, more particles will stay in the same cell, hence more cost will occur when computing neighbor particles (see picture below). This profile record shows `kernUpdateVelNeighborSearchCoherent` takes 92% of time per frame, which means a lot of cost when a particle finding its neighbor to update its velocity. + +![](/profiles/1MillionWith10KGrid.png) + +The grey line shows when there are too many grid cells. This leads to unnecessary cost for sorting the grid index. + +![](/profiles/1MillionWith515Million.png) + +The green line is more optimal setting compared to base line. It has a more reasonalbe particle / cell ratio which significantly improve FPS. + +### Future Works + +Future works may contains testing on more initial settings, with different grid sizes or particles. Also, interact ranges for particles might be another factors that affect the performance, since if particle pull other neighbors around them to themselves, there will be a dense area that is full of particles, which may drag down performance. \ No newline at end of file diff --git a/profiles/11.gif b/profiles/11.gif new file mode 100644 index 0000000..9b79de1 Binary files /dev/null and b/profiles/11.gif differ diff --git a/profiles/1MillionWith10KGrid.png b/profiles/1MillionWith10KGrid.png new file mode 100644 index 0000000..db7297c Binary files /dev/null and b/profiles/1MillionWith10KGrid.png differ diff --git a/profiles/1MillionWith515Million.png b/profiles/1MillionWith515Million.png new file mode 100644 index 0000000..fb197d3 Binary files /dev/null and b/profiles/1MillionWith515Million.png differ diff --git a/profiles/FrameRateFactors.png b/profiles/FrameRateFactors.png new file mode 100644 index 0000000..ac6bceb Binary files /dev/null and b/profiles/FrameRateFactors.png differ diff --git a/profiles/OtherFactors.png b/profiles/OtherFactors.png new file mode 100644 index 0000000..9062e36 Binary files /dev/null and b/profiles/OtherFactors.png differ diff --git a/profiles/ParticleCounts-FPS.png b/profiles/ParticleCounts-FPS.png new file mode 100644 index 0000000..112cab7 Binary files /dev/null and b/profiles/ParticleCounts-FPS.png differ diff --git a/profiles/ProfileBaseLine.png b/profiles/ProfileBaseLine.png new file mode 100644 index 0000000..a35b834 Binary files /dev/null and b/profiles/ProfileBaseLine.png differ diff --git a/profiles/ProfileSheet.png b/profiles/ProfileSheet.png new file mode 100644 index 0000000..c7c0f27 Binary files /dev/null and b/profiles/ProfileSheet.png differ diff --git a/profiles/profile_sheet.ods b/profiles/profile_sheet.ods new file mode 100644 index 0000000..bd4a860 Binary files /dev/null and b/profiles/profile_sheet.ods differ diff --git a/src/boidImpl.cu b/src/boidImpl.cu new file mode 100644 index 0000000..3fafb1e --- /dev/null +++ b/src/boidImpl.cu @@ -0,0 +1,397 @@ +#include "boidImpl.h" + +#include +#include +#include +#include + + + +// 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 +// order for iterating over neighboring grid cells? +// for(x) +// 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; +} + +// 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; + } +} + +/****************** +* stepSimulation * +******************/ + +/** +* LOOK-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 + glm::vec3 BoidPos = pos[iSelf]; + glm::vec3 AccumulatedPos = glm::vec3(0.0f, 0.0f, 0.0f); + glm::vec3 Propulsion = glm::vec3(0.0f, 0.0f, 0.0f); + glm::vec3 AccumulatedVelo = glm::vec3(0.0f, 0.0f, 0.0f); + int Counts1 = 0; + int Counts3 = 0; + for (int i = 0; i < N; i++) + { + glm::vec3 CurPos = pos[i]; + float Dist = glm::length(BoidPos - CurPos); + if (i != iSelf && Dist < rule1Distance) + { + AccumulatedPos += CurPos; + Counts1++; + } + if (i != iSelf && Dist < rule2Distance) + { + Propulsion += (BoidPos - CurPos); + } + if (i != iSelf && Dist < rule3Distance) + { + AccumulatedVelo += vel[i]; + Counts3++; + } + } + + glm::vec3 dvel; + if (Counts1>0) + { + dvel += rule1Scale * (AccumulatedPos / static_cast(Counts1) - BoidPos); + } + if (Counts3>0) + { + dvel += rule3Scale * AccumulatedVelo / static_cast(Counts3); + } + dvel += rule2Scale * Propulsion; + return dvel; +} + +/** +* 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? + + // Update position by velocity + int index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index >= N) { + return; + } + vel2[index] = vel1[index] + computeVelocityChange(N, index, pos, vel1); + if (glm::length(vel2[index]) > maxSpeed) + { + vel2[index] = glm::normalize(vel2[index]) * maxSpeed; + } +} + +__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 tid = threadIdx.x + (blockIdx.x * blockDim.x); + if (tid >= N) { + return; + } + if (tid==0) + { + gridCellStartIndices[0] = 0; + } + if(tid == N - 1) + { + gridCellEndIndices[particleGridIndices[tid]] = tid + 1; + } + if (tid > 0 && particleGridIndices[tid] != particleGridIndices[tid - 1]) + { + gridCellStartIndices[particleGridIndices[tid]] = tid; + gridCellEndIndices[particleGridIndices[tid - 1]] = tid; + } +} + +__device__ void accumulateOneCell( + int tid, int iX, int iY, int iZ, + int startIdx, int endIdx, + int* particleArrayIndices, + glm::vec3* pos, glm::vec3* vel1, glm::vec3* col, + glm::vec3& InOutAccuVelo1, glm::vec3& InOutAccuVelo2, glm::vec3& InOutAccuVelo3, + int& InOutCount1, int& InOutCount2, int& InOutCount3) +{ + glm::vec3 BoidPos = pos[tid]; + for (int i = startIdx; i < endIdx; i++) + { + int bufferIdx = particleArrayIndices[i]; + glm::vec3 CurPos = pos[bufferIdx]; + float Dist = glm::length(BoidPos - CurPos); + if (tid != bufferIdx && Dist < rule1Distance) + { + InOutAccuVelo1 += CurPos; + InOutCount1++; + } + if (tid != bufferIdx && Dist < rule2Distance) + { + InOutAccuVelo2 += (BoidPos - CurPos); + InOutCount2++; + } + if (tid != bufferIdx && Dist < rule3Distance) + { + InOutAccuVelo3 += vel1[i]; + InOutCount3++; + } + + col[bufferIdx].x = iX / 10 % 2; + col[bufferIdx].y = 0; + col[bufferIdx].z = 0; + } +} + +__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, glm::vec3* col) { + // TODO-2.1 - Update a boid's velocity using the uniform grid to reduce + // the number of boids that need to be checked. + + int tid = threadIdx.x + (blockIdx.x * blockDim.x); + if (tid >= N) { + return; + } + glm::vec3 BoidPos = pos[tid]; + glm::vec3 AccumulatedPos = glm::vec3(0.0f, 0.0f, 0.0f); + glm::vec3 Propulsion = glm::vec3(0.0f, 0.0f, 0.0f); + glm::vec3 AccumulatedVelo = glm::vec3(0.0f, 0.0f, 0.0f); + int Counts1 = 0; + int Counts2 = 0; + int Counts3 = 0; + int iX = floor((BoidPos.x - gridMin.x) * inverseCellWidth); + int iY = floor((BoidPos.y - gridMin.y) * inverseCellWidth); + int iZ = floor((BoidPos.z - gridMin.z) * inverseCellWidth); + + // - Identify the grid cell that this particle is in + int maxGridIndex = gridIndex3Dto1D(gridResolution, gridResolution, gridResolution, gridResolution) - 1; +#ifndef HALF_GRID_WIDTH + glm::vec3 posInCell; + posInCell.x = BoidPos.x - (iX * cellWidth + gridMin.x); + posInCell.y = BoidPos.y - (iY * cellWidth + gridMin.y); + posInCell.z = BoidPos.z - (iZ * cellWidth + gridMin.z); + int dX = (posInCell.x < cellWidth * 0.5f) ? -1 : 1; + int dY = (posInCell.y < cellWidth * 0.5f) ? -1 : 1; + int dZ = (posInCell.z < cellWidth * 0.5f) ? -1 : 1; + + int cellIndices[] = { + + gridIndex3Dto1D(iX, iY, iZ, gridResolution), + gridIndex3Dto1D(iX+dX, iY, iZ, gridResolution), + gridIndex3Dto1D(iX, iY+ dY, iZ, gridResolution), + gridIndex3Dto1D(iX+ dX,iY+ dY, iZ, gridResolution), + gridIndex3Dto1D(iX, iY, iZ+ dZ, gridResolution), + gridIndex3Dto1D(iX+ dX,iY, iZ+ dZ, gridResolution), + gridIndex3Dto1D(iX, iY+ dY, iZ+ dZ, gridResolution), + gridIndex3Dto1D(iX+ dX,iY+ dY, iZ+ dZ, gridResolution), + }; +#else + int cellIndices[27]; + int cellId = 0; + for (int dX : {-1, 0, 1}) + { + for (int dY : {-1, 0, 1}) + { + for (int dZ : {-1, 0, 1}) + { + cellIndices[cellId++] = gridIndex3Dto1D(iX + dX, iY + dY, iZ + dZ, gridResolution); + } + } + } +#endif + for (int i = 0;i<8;i++) + { + if (cellIndices[i] < 0 || cellIndices[i]>maxGridIndex) + { + continue; + } + accumulateOneCell(tid, iX, iY, iZ, + gridCellStartIndices[cellIndices[i]], + gridCellEndIndices[cellIndices[i]], particleArrayIndices, pos, vel1, col, + AccumulatedPos, Propulsion, AccumulatedVelo, + Counts1, Counts2, Counts3); + } + + glm::vec3 newVel = vel1[tid]; + if (Counts1>0) + { + newVel += rule1Scale * (AccumulatedPos / static_cast(Counts1) - BoidPos); + } + if (Counts2>0) + { + newVel += rule2Scale * Propulsion; + } + if (Counts3 > 0) + { + newVel += rule3Scale * AccumulatedVelo / static_cast(Counts3); + } + if (glm::length(newVel) > maxSpeed) + { + newVel = glm::normalize(newVel) * maxSpeed; + } + vel2[tid] = newVel; + // - 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 +} + +__device__ void accumulateOneCellCoherent( + int tid, int iX, int iY, int iZ, + int startIdx, int endIdx, + glm::vec3* pos, glm::vec3* vel1, + glm::vec3& InOutAccuVelo1, glm::vec3& InOutAccuVelo2, glm::vec3& InOutAccuVelo3, + int& InOutCount1, int& InOutCount2, int& InOutCount3) +{ + glm::vec3 BoidPos = pos[tid]; + for (int i = startIdx; i < endIdx; i++) + { + glm::vec3 CurPos = pos[i]; + float Dist = glm::length(BoidPos - CurPos); + if (tid != i && Dist < rule1Distance) + { + InOutAccuVelo1 += CurPos; + InOutCount1++; + } + if (tid != i && Dist < rule2Distance) + { + InOutAccuVelo2 -= (CurPos -BoidPos); + InOutCount2++; + } + if (tid != i && Dist < rule3Distance) + { + InOutAccuVelo3 += vel1[i]; + InOutCount3++; + } + } +} + +__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 tid = threadIdx.x + (blockIdx.x * blockDim.x); + if (tid >= N) { + return; + } + glm::vec3 BoidPos = pos[tid]; + glm::vec3 AccumulatedPos = glm::vec3(0.0f, 0.0f, 0.0f); + glm::vec3 Propulsion = glm::vec3(0.0f, 0.0f, 0.0f); + glm::vec3 AccumulatedVelo = glm::vec3(0.0f, 0.0f, 0.0f); + int Counts1 = 0; + int Counts2 = 0; + int Counts3 = 0; + int iX = floor((BoidPos.x - gridMin.x) * inverseCellWidth); + int iY = floor((BoidPos.y - gridMin.y) * inverseCellWidth); + int iZ = floor((BoidPos.z - gridMin.z) * inverseCellWidth); + + // - Identify the grid cell that this particle is in + // current cell + + glm::vec3 posInCell; + posInCell.x = BoidPos.x - (iX * cellWidth + gridMin.x); + posInCell.y = BoidPos.y - (iY * cellWidth + gridMin.y); + posInCell.z = BoidPos.z - (iZ * cellWidth + gridMin.z); + + int dX = (posInCell.x < cellWidth * 0.5f) ? -1 : 1; + int dY = (posInCell.y < cellWidth * 0.5f) ? -1 : 1; + int dZ = (posInCell.z < cellWidth * 0.5f) ? -1 : 1; + + int maxGridIndex = gridIndex3Dto1D(gridResolution, gridResolution, gridResolution, gridResolution) - 1; + + int cellIndices[] = { + gridIndex3Dto1D(iX, iY, iZ, gridResolution), + gridIndex3Dto1D(iX + dX, iY, iZ, gridResolution), + gridIndex3Dto1D(iX, iY + dY, iZ, gridResolution), + gridIndex3Dto1D(iX + dX, iY + dY, iZ, gridResolution), + gridIndex3Dto1D(iX, iY, iZ + dZ, gridResolution), + gridIndex3Dto1D(iX + dX, iY, iZ + dZ, gridResolution), + gridIndex3Dto1D(iX, iY + dY, iZ + dZ, gridResolution), + gridIndex3Dto1D(iX + dX, iY + dY, iZ + dZ, gridResolution), + }; + + for (int i = 0; i < 8; i++) + { + if (cellIndices[i] < 0 || cellIndices[i]>maxGridIndex) + { + continue; + } + accumulateOneCellCoherent(tid, iX, iY, iZ, + gridCellStartIndices[cellIndices[i]], + gridCellEndIndices[cellIndices[i]], pos, vel1, + AccumulatedPos, Propulsion, AccumulatedVelo, + Counts1, Counts2, Counts3); + } + glm::vec3 newVel; + newVel += vel1[tid]; + if (Counts1 > 0) + { + newVel += rule1Scale * (AccumulatedPos / static_cast(Counts1) - BoidPos); + } + newVel += rule2Scale * Propulsion; + if (Counts3 > 0) + { + newVel += rule3Scale * AccumulatedVelo / static_cast(Counts3); + } + if (glm::length(newVel) > maxSpeed) + { + newVel = glm::normalize(newVel) * maxSpeed; + } + vel2[tid] = newVel; + +} +__global__ void kernComputeIndices(int N, int gridResolution, + glm::vec3 gridMin, float inverseCellWidth, + glm::vec3* pos, int* particlePropertyIndex, int* particleGridIndex, glm::vec3* col) { + // TODO-2. + // - Label each boid with the index of its grid cell. + int tid = threadIdx.x + (blockIdx.x * blockDim.x); + if (tid >= N) { + return; + } + int bufferIndex = particlePropertyIndex[tid]; + int iX = floor((pos[bufferIndex].x - gridMin.x) * inverseCellWidth); + int iY = floor((pos[bufferIndex].y - gridMin.y) * inverseCellWidth); + int iZ = floor((pos[bufferIndex].z - gridMin.z) * inverseCellWidth); + // - Set up a parallel array of integer indices as pointers to the actual + // boid data in pos and vel1/vel2 + particleGridIndex[tid] = gridIndex3Dto1D(iX, iY, iZ, gridResolution); +} \ No newline at end of file diff --git a/src/boidImpl.h b/src/boidImpl.h new file mode 100644 index 0000000..4a4fd3b --- /dev/null +++ b/src/boidImpl.h @@ -0,0 +1,61 @@ +#include +#include <../device_launch_parameters.h> + +#define HALF_GRID_WIDTH + +// LOOK-2.1 potentially useful for doing grid-based neighbor search +#ifndef imax +#define imax( a, b ) ( ((a) > (b)) ? (a) : (b) ) +#endif + +#ifndef imin +#define imin( a, b ) ( ((a) < (b)) ? (a) : (b) ) +#endif + +#define checkCUDAErrorWithLine(msg) checkCUDAError(msg, __LINE__) + + +/***************** +* Configuration * +*****************/ + +/*! Block size used for CUDA kernel launch. */ +#define blockSize 1024 + +#define test_scale 2.f +// LOOK-1.2 Parameters for the boids algorithm. +// These worked well in our reference implementation. +#define rule1Distance (5.0f * test_scale) +#define rule2Distance (3.0f * test_scale) +#define rule3Distance (5.0f * test_scale) + +#define rule1Scale 0.01f +#define rule2Scale 0.1f +#define rule3Scale 0.1f + +#define maxSpeed 1.0f + +/*! Size of the starting area in simulation space. */ +#define scene_scale (100.0f * test_scale) + +#define GLM_FORCE_CUDA +#include +__global__ void kernUpdateVelocityBruteForce(int N, glm::vec3* pos, + glm::vec3* vel1, glm::vec3* vel2); + +__global__ void kernComputeIndices(int N, int gridResolution, + glm::vec3 gridMin, float inverseCellWidth, + glm::vec3* pos, int* particlePropertyIndex, int* particleGridIndex, glm::vec3* col); +__global__ void kernIdentifyCellStartEnd(int N, int* particleGridIndices, + int* gridCellStartIndices, int* gridCellEndIndices); +__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, glm::vec3* col); +__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); diff --git a/src/kernel.cu b/src/kernel.cu index 7149917..b8cd244 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -2,6 +2,7 @@ #include #include "kernel.h" +#include "boidImpl.h" #include "utilityCore.hpp" #include @@ -16,16 +17,6 @@ #include -// LOOK-2.1 potentially useful for doing grid-based neighbor search -#ifndef imax -#define imax( a, b ) ( ((a) > (b)) ? (a) : (b) ) -#endif - -#ifndef imin -#define imin( a, b ) ( ((a) < (b)) ? (a) : (b) ) -#endif - -#define checkCUDAErrorWithLine(msg) checkCUDAError(msg, __LINE__) /** * Check for CUDA errors; print and exit if there was a problem. @@ -42,28 +33,6 @@ void checkCUDAError(const char *msg, int line = -1) { } -/***************** -* Configuration * -*****************/ - -/*! Block size used for CUDA kernel launch. */ -#define blockSize 128 - -// LOOK-1.2 Parameters for the boids algorithm. -// These worked well in our reference implementation. -#define rule1Distance 5.0f -#define rule2Distance 3.0f -#define rule3Distance 5.0f - -#define rule1Scale 0.01f -#define rule2Scale 0.1f -#define rule3Scale 0.1f - -#define maxSpeed 1.0f - -/*! Size of the starting area in simulation space. */ -#define scene_scale 100.0f - /*********************************************** * Kernel state (pointers are device pointers) * ***********************************************/ @@ -79,6 +48,7 @@ dim3 threadsPerBlock(blockSize); glm::vec3 *dev_pos; glm::vec3 *dev_vel1; glm::vec3 *dev_vel2; +glm::vec3* dev_col; // LOOK-2.1 - these are NOT allocated for you. You'll have to set up the thrust // pointers on your own too. @@ -167,11 +137,17 @@ void Boids::initSimulation(int N) { checkCUDAErrorWithLine("kernGenerateRandomPosArray failed!"); // LOOK-2.1 computing grid params - gridCellWidth = 2.0f * std::max(std::max(rule1Distance, rule2Distance), rule3Distance); +#ifdef HALF_GRID_WIDTH + gridCellWidth = 5.f; +#else + gridCellWidth = 2.0f * 5.f; +#endif int halfSideCount = (int)(scene_scale / gridCellWidth) + 1; gridSideCount = 2 * halfSideCount; gridCellCount = gridSideCount * gridSideCount * gridSideCount; + std::cout << "[gridCellCount] " << gridCellCount << std::endl; + std::cout << "[particleSpawn] " << numObjects << std::endl; gridInverseCellWidth = 1.0f / gridCellWidth; float halfGridWidth = gridCellWidth * halfSideCount; gridMinimum.x -= halfGridWidth; @@ -179,6 +155,32 @@ void Boids::initSimulation(int N) { gridMinimum.z -= halfGridWidth; // TODO-2.1 TODO-2.3 - Allocate additional buffers here. + + + cudaMalloc((void**)&dev_col, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_col failed!"); + + + cudaMalloc((void**)&dev_particleGridIndices, N * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_particleGridIndices failed!"); + + cudaMalloc((void**)&dev_particleArrayIndices, N * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_particleArrayIndices 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!"); + + + std::vector particleArrayIndex; + particleArrayIndex.resize(numObjects); + for (int i = 0; i < numObjects; i++) + { + particleArrayIndex[i] = i; + } + cudaMemcpy(dev_particleArrayIndices, particleArrayIndex.data(), sizeof(int) * numObjects, cudaMemcpyHostToDevice); + cudaDeviceSynchronize(); } @@ -221,42 +223,13 @@ void Boids::copyBoidsToVBO(float *vbodptr_positions, float *vbodptr_velocities) dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); kernCopyPositionsToVBO << > >(numObjects, dev_pos, vbodptr_positions, scene_scale); - kernCopyVelocitiesToVBO << > >(numObjects, dev_vel1, vbodptr_velocities, scene_scale); + kernCopyVelocitiesToVBO << > >(numObjects, dev_col, vbodptr_velocities, scene_scale); checkCUDAErrorWithLine("copyBoidsToVBO failed!"); cudaDeviceSynchronize(); } - -/****************** -* stepSimulation * -******************/ - -/** -* LOOK-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); -} - -/** -* 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? -} - /** * LOOK-1.2 Since this is pretty trivial, we implemented it for you. * For each of the `N` bodies, update its position based on its current velocity. @@ -282,76 +255,6 @@ __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 -// order for iterating over neighboring grid cells? -// for(x) -// 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; -} - -__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 -} - -// LOOK-2.1 Consider how this could be useful for indicating that a cell -// does not enclose any boids -__global__ void kernResetIntBuffer(int N, int *intBuffer, int value) { - int index = (blockIdx.x * blockDim.x) + threadIdx.x; - if (index < N) { - intBuffer[index] = value; - } -} - -__global__ void 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!" -} - -__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 -} - -__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 -} /** * Step the entire N-body simulation by `dt` seconds. @@ -359,6 +262,15 @@ __global__ void kernUpdateVelNeighborSearchCoherent( void Boids::stepSimulationNaive(float dt) { // TODO-1.2 - use the kernels you wrote to step the simulation forward in time. // TODO-1.2 ping-pong the velocity buffers + glm::vec3* temp = dev_vel1; + dev_vel1 = dev_vel2; + dev_vel2 = temp; + dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); + kernUpdateVelocityBruteForce <<>> (numObjects, dev_pos, dev_vel1, dev_vel2); + kernUpdatePos <<>> (numObjects, dt, dev_pos, dev_vel2); + + cudaMemcpy(dev_col, dev_vel1, sizeof(glm::vec3) * numObjects, cudaMemcpyDeviceToDevice); + checkCUDAErrorWithLine("memcpy back failed!"); } void Boids::stepSimulationScatteredGrid(float dt) { @@ -369,11 +281,43 @@ void Boids::stepSimulationScatteredGrid(float dt) { // 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 + + dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); + kernComputeIndices <<>> (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, + dev_pos, dev_particleArrayIndices, dev_particleGridIndices, dev_col); + + checkCUDAErrorWithLine("kernComputeIndices failed!"); + thrust::device_ptr dev_thrust_keys(dev_particleGridIndices); + thrust::device_ptr dev_thrust_values(dev_particleArrayIndices); + thrust::sort_by_key(dev_thrust_keys, dev_thrust_keys + numObjects, dev_thrust_values); + checkCUDAErrorWithLine("thrust::sort_by_key failed!"); + + cudaMemset(dev_gridCellStartIndices, 0, sizeof(int) * gridCellCount); + cudaMemset(dev_gridCellEndIndices, 0, sizeof(int) * gridCellCount); + // - Ping-pong buffers as needed + + glm::vec3* temp = dev_vel1; + dev_vel1 = dev_vel2; + dev_vel2 = temp; + // - 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); + checkCUDAErrorWithLine("kernIdentifyCellStartEnd failed!"); + + // - Perform velocity updates using neighbor search + kernUpdateVelNeighborSearchScattered <<>>(numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, gridCellWidth, + dev_gridCellStartIndices, dev_gridCellEndIndices, dev_particleArrayIndices, + dev_pos, dev_vel1, dev_vel2, dev_col); + checkCUDAErrorWithLine("kernUpdateVelNeighborSearchScattered failed!"); + // - Update positions + kernUpdatePos <<>> (numObjects, dt, dev_pos, dev_vel1); + checkCUDAErrorWithLine("kernUpdatePos failed!"); + + cudaMemcpy(dev_col, dev_vel1, sizeof(glm::vec3) * numObjects, cudaMemcpyDeviceToDevice); + checkCUDAErrorWithLine("cudaMemcpy dev_col failed!"); + + } void Boids::stepSimulationCoherentGrid(float dt) { @@ -392,76 +336,118 @@ void Boids::stepSimulationCoherentGrid(float dt) { // - Perform velocity updates using neighbor search // - Update positions // - Ping-pong buffers as needed. THIS MAY BE DIFFERENT FROM BEFORE. + + dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); + kernComputeIndices << > > (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, + dev_pos, dev_particleArrayIndices, dev_particleGridIndices, dev_col); + + glm::vec3* temp = dev_vel1; + dev_vel1 = dev_vel2; + dev_vel2 = temp; + checkCUDAErrorWithLine("kernComputeIndices failed!"); + thrust::device_ptr dev_thrust_keys(dev_particleGridIndices); + thrust::device_ptr dev_thrust_vel1(dev_vel1); + thrust::device_ptr dev_thrust_pos(dev_pos); + thrust::sort_by_key(dev_thrust_keys, dev_thrust_keys + numObjects, thrust::make_zip_iterator(dev_thrust_vel1, dev_thrust_pos)); + checkCUDAErrorWithLine("thrust::sort_by_key failed!"); + + cudaMemset(dev_gridCellStartIndices, 0, sizeof(int) * gridCellCount); + cudaMemset(dev_gridCellEndIndices, 0, sizeof(int) * gridCellCount); + + // - Ping-pong buffers as needed + + // - 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); + checkCUDAErrorWithLine("kernIdentifyCellStartEnd failed!"); + + // - Perform velocity updates using neighbor search + kernUpdateVelNeighborSearchCoherent <<>> (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, gridCellWidth, + dev_gridCellStartIndices, dev_gridCellEndIndices, + dev_pos, dev_vel1, dev_vel2); + checkCUDAErrorWithLine("kernUpdateVelNeighborSearchScattered failed!"); + // - Update positions + kernUpdatePos << > > (numObjects, dt, dev_pos, dev_vel2); + checkCUDAErrorWithLine("kernUpdatePos failed!"); + + cudaMemcpy(dev_col, dev_vel1, sizeof(glm::vec3) * numObjects, cudaMemcpyDeviceToDevice); + checkCUDAErrorWithLine("cudaMemcpy dev_col failed!"); + } void Boids::endSimulation() { cudaFree(dev_vel1); cudaFree(dev_vel2); cudaFree(dev_pos); + cudaFree(dev_col); + cudaFree(dev_particleGridIndices); + cudaFree(dev_particleArrayIndices); + cudaFree(dev_gridCellStartIndices); + cudaFree(dev_gridCellEndIndices); // TODO-2.1 TODO-2.3 - Free any additional buffers here. } 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/kernel.h b/src/kernel.h index a38b64d..99bfa1b 100644 --- a/src/kernel.h +++ b/src/kernel.h @@ -1,5 +1,6 @@ #pragma once + namespace Boids { void initSimulation(int N); void stepSimulationNaive(float dt); diff --git a/src/main.cpp b/src/main.cpp index 9c917c0..797b29c 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -7,6 +7,9 @@ */ #include "main.hpp" + +#include + #include "kernel.h" #include @@ -15,6 +18,7 @@ #include #include +#include #include // ================ @@ -22,14 +26,14 @@ // ================ // 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 VISUALIZE 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; - /** * C main function. */ @@ -226,26 +230,29 @@ void initShaders(GLuint * program) { double timebase = 0; int frame = 0; + auto start = std::chrono::high_resolution_clock::now(); Boids::unitTest(); // LOOK-1.2 We run some basic example code to make sure // your CUDA development setup is ready to go. while (!glfwWindowShouldClose(window)) { glfwPollEvents(); - + //system("pause"); frame++; double time = glfwGetTime(); - if (time - timebase > 1.0) { + if (time - timebase > 3.0) { fps = frame / (time - timebase); timebase = time; frame = 0; } - + // auto duration = std::chrono::duration_cast(std::chrono::high_resolution_clock::now() - start); + //start = std::chrono::high_resolution_clock::now(); + //std::cout << duration.count() / 1000.0f << " ms" << std::endl; runCUDA(); std::ostringstream ss; ss << "["; - ss.precision(1); + ss.precision(3); ss << std::fixed << fps; ss << " fps] " << deviceName; glfwSetWindowTitle(window, ss.str().c_str());