diff --git a/README.md b/README.md index 98dd9a8..a806190 100644 --- a/README.md +++ b/README.md @@ -1,10 +1,47 @@ **University of Pennsylvania, CIS 565: GPU Programming and Architecture, Project 1 - Flocking** -* (TODO) YOUR NAME HERE -* Tested on: (TODO) Windows 22, i7-2222 @ 2.22GHz 22GB, GTX 222 222MB (Moore 2222 Lab) +* (TODO) Anton Khabbaz +* Pennkey: akhabbaz +* Tested on: Windows 10 , i7-6600U @ 2.6 GHz 16 GB, GTX 965M 10000MB (personal Windows SUrface book) ### (TODO: Your README) +// update src/CMakeLists.txt with '-arch=sm_52' since that is the cuda version supported on my graphics card. +p +I added a line to increase the frame rate by turning off vsync with +glfwSwapInterval(false) (Josh helped me). + + +Got 1.2 and the flocking to work in the naive case. I changed parameters +and found that with 1000 boids I got a frame rate of 26 fps with my gpu. + +I wrote the code for 2.1 the grid search. + +I determined that the gridCells and the start and endIndices are correct. This is all the input needed to update the velocity. +I am sure that the indices. I ran the code with 100 boids and a 20 sized max radius. These parameters worked for the naive grid and had the boids converge. I used those to get the uniform grid scatter to work. With those, I wrote a Boid::testGRidArray that transfeerred the intArrays to the host and printed them out nicely. + +Those indicated that the Start, End grid indices were correct. Those algorithms that I wrote were all parallel, and had no loops over all the Boids. +For example to find the grid Start and stop I compared each gridCellindex and its predecessor. There was no loop over all the Boids. + +THere was a tricky routine to getting the neigboring cells. I used the maximum of rule1 -3 distances, called maxdist and I added (maxdist, maxdist, maxdist) to the current Boid position. This was the upper limit of the gridcell that was needed. I then realized that the 4 highest cells are next to (-1 in x, -1 in y, (-1, -1) in x, y) this cell. Similarly I subtracted the same offset. Since the cell width is 2 * maxdist, this measure will get the 8 cells and no grid cells outside of the core 8. This created an array of possible grid Cells Each thread pulled the boids in order of the grid Cells. This way all boids in a cell are processed together and this seemed to help the speed. + +A debugging issue that really slowed me down by about a day was that I could not break in any kernel inucluding the update velocity kernels using Nsight. More specifically, I could break prior to running thrust, but afterwards the break points caused Visual Studio to stop with an error. I checked the arrays with print statements and all the values were fine. The debugging failure also occured in the naive implementation and that for sure was working. +Many people helped today including a TA, but we could not find the bug. + +By selectively getting rule1 only to work and comparing with the naive solution, and by studying the code line by line, I traced down many bugs and got the grid Scattering to work. + +The coherent Grid was straightforward. I wrote a parallel kernel that would rearrange the position and velocity buffers according to the particle array buffer. One extra position buffer was needed to prevent writing over values. + +Results: +![](images/boids.png) +![](images/boids2.png) +![](images/resultsProfile.pdf) +See the images and also the results of several time tests. +By considering all the boids in nearest neighbor, the rate was about 26 fps with 10000 boids (7 = maxdist). By using the uniform grid with Scattering and taking processing nearest neighbors at a time, I got a rate of about 205 fps. By rearranging so that the boids were coherent and removing the indices, The reate improved to about 207 fps (slight improvement). The uniform grid sped up the processing by about a factor of 8. The coherent grid did better because boids in one grid were all in contiguous memory. The scatter was nearly as good. Boids there were still processed in grid order and I think for that reason the data was cached and so the processing time was nearly as good. + +Increasing the number of boids sightly reduced processing time: at 50000 the rate was about 156, about a 25% reduction in frame rate. I think the effect is low because of the parallel architecture. tripling the grid size (from 7 to 20) reduced the rate by over a factor of 9; now there are about 27 times as many neighbors, so it is expected that the rate would be reduced by a huge factor. Since there is only 8 grid cells that could possibly contribute, increasing the grid count would only slow down the frame rate further. Include screenshots, analysis, etc. (Remember, this is public, so don't put anything here that you don't want to share with the world.) + +: diff --git a/Recitation 1.pdf b/Recitation 1.pdf new file mode 100644 index 0000000..d0c5856 Binary files /dev/null and b/Recitation 1.pdf differ diff --git a/images/boids.png b/images/boids.png new file mode 100644 index 0000000..5193faf Binary files /dev/null and b/images/boids.png differ diff --git a/images/boids2.png b/images/boids2.png new file mode 100644 index 0000000..20cf333 Binary files /dev/null and b/images/boids2.png differ diff --git a/images/profileResults.pdf b/images/profileResults.pdf new file mode 100644 index 0000000..036ccfc Binary files /dev/null and b/images/profileResults.pdf differ diff --git a/profileResults.xlsx b/profileResults.xlsx new file mode 100644 index 0000000..e4eac88 Binary files /dev/null and b/profileResults.xlsx differ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index fdd636d..eeaabd4 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -10,5 +10,5 @@ set(SOURCE_FILES cuda_add_library(src ${SOURCE_FILES} - OPTIONS -arch=sm_20 + OPTIONS -arch=sm_52 ) diff --git a/src/kernel.cu b/src/kernel.cu index aaf0fbf..cbd99ff 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -2,9 +2,13 @@ #include #include #include +#include #include #include "utilityCore.hpp" #include "kernel.h" +#include +#include +#include // LOOK-2.1 potentially useful for doing grid-based neighbor search #ifndef imax @@ -41,9 +45,9 @@ void checkCUDAError(const char *msg, int line = -1) { // 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 rule1Distance 20.0f +#define rule2Distance 4.0f +#define rule3Distance 10.0f #define rule1Scale 0.01f #define rule2Scale 0.1f @@ -53,6 +57,9 @@ void checkCUDAError(const char *msg, int line = -1) { /*! Size of the starting area in simulation space. */ #define scene_scale 100.0f +// this is the number of cells that could be neighbors +// add 1 to indicate the end of the array. +#define maxGridCells 9 /*********************************************** * Kernel state (pointers are device pointers) * @@ -66,7 +73,7 @@ dim3 threadsPerBlock(blockSize); // Consider why you would need two velocity buffers in a simulation where each // boid cares about its neighbors' velocities. // These are called ping-pong buffers. -glm::vec3 *dev_pos; +glm::vec3 *dev_pos; glm::vec3 *dev_vel1; glm::vec3 *dev_vel2; @@ -85,7 +92,10 @@ 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. - +// during the swap, we need to rearrange the positions. This involves copying +// positions from one location to another. A second buffer is needed, not just a single +// glm::vec3 position because it is not a simple swap. +glm::vec3 *dev_posCopy; // LOOK-2.1 - Grid parameters based on simulation parameters. // These are automatically computed for you in Boids::initSimulation int gridCellCount; @@ -94,6 +104,14 @@ float gridCellWidth; float gridInverseCellWidth; glm::vec3 gridMinimum; + +// print index Array matrices + +static const int numPerRow{20}; +static const int spacePerInt{5}; +static const int numPerRowVec3{3}; +static const int spacePerFloat{10}; +static const int precision{2}; /****************** * initSimulation * ******************/ @@ -119,6 +137,7 @@ __host__ __device__ glm::vec3 generateRandomVec3(float time, int index) { return glm::vec3((float)unitDistrib(rng), (float)unitDistrib(rng), (float)unitDistrib(rng)); } + /** * LOOK-1.2 - This is a basic CUDA kernel. * CUDA kernel for generating boids with a specified mass randomly around the star. @@ -147,9 +166,13 @@ void Boids::initSimulation(int N) { cudaMalloc((void**)&dev_vel1, N * sizeof(glm::vec3)); checkCUDAErrorWithLine("cudaMalloc dev_vel1 failed!"); - + cudaMemset(dev_vel1, 0.0, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("Initializing dev_vel1 failed"); + cudaMalloc((void**)&dev_vel2, N * sizeof(glm::vec3)); checkCUDAErrorWithLine("cudaMalloc dev_vel2 failed!"); + cudaMemset(dev_vel2, 0.0, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("Initializing dev_vel2 failed"); // LOOK-1.2 - This is a typical CUDA kernel invocation. kernGenerateRandomPosArray<<>>(1, numObjects, @@ -157,18 +180,40 @@ void Boids::initSimulation(int N) { 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; + gridCellWidth = 2.0f * std::max({rule1Distance, rule2Distance, + rule3Distance}); + gridInverseCellWidth = 1.0f / gridCellWidth; + int halfSideCount = static_cast(scene_scale * gridInverseCellWidth) + 1; gridSideCount = 2 * halfSideCount; gridCellCount = gridSideCount * gridSideCount * gridSideCount; - gridInverseCellWidth = 1.0f / gridCellWidth; float halfGridWidth = gridCellWidth * halfSideCount; gridMinimum.x -= halfGridWidth; gridMinimum.y -= halfGridWidth; gridMinimum.z -= halfGridWidth; // TODO-2.1 TODO-2.3 - Allocate additional buffers here. + // allocation for 2.1 arrays + cudaMalloc(reinterpret_cast(&dev_particleArrayIndices), numObjects * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_particleArrayIndices failed!"); + + cudaMalloc(reinterpret_cast(&dev_particleGridIndices), numObjects * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_particleGridIndices failed!"); + + cudaMalloc(reinterpret_cast(&dev_gridCellStartIndices), gridCellCount * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc devgridCellStartIndices failed!"); + + cudaMalloc(reinterpret_cast(&dev_gridCellEndIndices), gridCellCount * sizeof(int)); + + + dev_thrust_particleArrayIndices = thrust::device_ptr(dev_particleArrayIndices); + dev_thrust_particleGridIndices = thrust::device_ptr(dev_particleGridIndices); + checkCUDAErrorWithLine("Assign dev_thrust_pointers failed!"); + + // allocate for the position copy + cudaMalloc((void**)&dev_posCopy, numObjects * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_posCopy failed!"); + cudaThreadSynchronize(); } @@ -229,11 +274,272 @@ void Boids::copyBoidsToVBO(float *vbodptr_positions, float *vbodptr_velocities) * Compute the new velocity on the body with index `iSelf` due to the `N` boids * in the `pos` and `vel` arrays. */ -__device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *pos, const glm::vec3 *vel) { - // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves - // Rule 2: boids try to stay a distance d away from each other +__device__ glm::vec3 rule1(int N, int iSelf, const glm::vec3 * pos, float scale) +{ + glm::vec3 pc(0.f, 0.0f, 0.0f); + int n {0}; + for ( int i {0}; i < iSelf; ++i){ + if ( glm::length(pos[i] - pos[iSelf]) < rule1Distance) { + pc += pos[i]; + ++n; + } + } + for( int i {iSelf + 1} ; i < N; ++i){ + if ( glm::length(pos[i] - pos[iSelf]) < rule1Distance) { + pc += pos[i]; + ++n; + } + } + if ( n == 0) + { + return glm::vec3(0.f); + } + pc /= n; + return (pc - pos[iSelf]) * scale; +} +__device__ glm::vec3 rule1Scatter(int iSelf, int const * gridCellStartIndices, + int const * gridCellEndIndices, int const * particleArrayIndices, + glm::vec3 const * pos, int const usedGridIndices[maxGridCells], + float scale) +{ + glm::vec3 pc(0.0f, 0.0f, 0.0f); + int n {0}; + for (int k{0}; usedGridIndices[k] != -1; ++k) + { + int gridC{ usedGridIndices[k] }; + //printf(" %d\n", gridC); + for ( int i {gridCellStartIndices[gridC]}; i < gridCellEndIndices[gridC]; ++i){ + int neighbor = particleArrayIndices[i]; + //printf("neighbor: %d\n", neighbor); + //int neighbor = i; + if ( neighbor == iSelf) { + continue; + } + if ( glm::length(pos[neighbor] - pos[iSelf]) < rule1Distance) { + pc += pos[neighbor]; + ++n; + } + } + } + if ( n == 0) + { + return glm::vec3(0.f); + } + pc /= n; + return (pc - pos[iSelf]) * scale; +} +__device__ glm::vec3 rule1Coherent(int iSelf, int const * gridCellStartIndices, + int const * gridCellEndIndices, + glm::vec3 const * pos, int const usedGridIndices[maxGridCells], + float scale) +{ + glm::vec3 pc(0.0f, 0.0f, 0.0f); + int n {0}; + for (int k{0}; usedGridIndices[k] != -1; ++k) + { + int gridC{ usedGridIndices[k] }; + //printf(" %d\n", gridC); + for ( int neighbor {gridCellStartIndices[gridC]}; neighbor < gridCellEndIndices[gridC]; ++neighbor){ + if ( neighbor == iSelf) { + continue; + } + if ( glm::length(pos[neighbor] - pos[iSelf]) < rule1Distance) { + pc += pos[neighbor]; + ++n; + } + } + } + if ( n == 0) + { + return glm::vec3(0.f); + } + pc /= n; + return (pc - pos[iSelf]) * scale; +} +__device__ glm::vec3 rule2(int N, int iSelf, const glm::vec3 * pos, float scale) +{ + glm::vec3 c{ 0.f }; + for ( int i {0}; i < iSelf; ++i){ + if ( glm::length(pos[i] - pos[iSelf]) < rule2Distance) { + c -= pos[i] - pos[iSelf]; + } + } + for( int i {iSelf + 1}; i < N; ++i){ + if ( glm::length(pos[i] - pos[iSelf]) < rule2Distance) { + c -= pos[i] - pos[iSelf]; + } + } + return c * scale; +} +__device__ glm::vec3 rule2Scatter(int iSelf, int const * gridCellStartIndices, + int const * gridCellEndIndices, int const * particleArrayIndices, + glm::vec3 const * pos, int const usedGridIndices[maxGridCells], + float scale) +{ + glm::vec3 c{ 0.f, 0.f, 0.f }; + for (int k{0}; usedGridIndices[k] != -1; ++k) + { + int gridC{ usedGridIndices[k] }; + for ( int i {gridCellStartIndices[gridC]}; i < gridCellEndIndices[gridC]; ++i){ + int neighbor = particleArrayIndices[i]; + if ( neighbor == iSelf) { + continue; + } + if ( glm::length(pos[neighbor] - pos[iSelf]) < rule2Distance) { + c -= pos[neighbor] - pos[iSelf]; + } + } + } + return c * scale; +} +__device__ glm::vec3 rule2Coherent(int iSelf, int const * gridCellStartIndices, + int const * gridCellEndIndices, + glm::vec3 const * pos, int const usedGridIndices[maxGridCells], + float scale) +{ + glm::vec3 c{ 0.f, 0.f, 0.f }; + for (int k{0}; usedGridIndices[k] != -1; ++k) + { + int gridC{ usedGridIndices[k] }; + for ( int neighbor {gridCellStartIndices[gridC]}; neighbor < gridCellEndIndices[gridC]; + ++neighbor){ + if ( neighbor == iSelf) { + continue; + } + if ( glm::length(pos[neighbor] - pos[iSelf]) < rule2Distance) { + c -= pos[neighbor] - pos[iSelf]; + } + } + } + return c * scale; +} +__device__ glm::vec3 rule3(int N, int iSelf, const glm::vec3 * pos, + const glm::vec3 *vel, float scale) +{ + glm::vec3 vsum{ 0.0f, 0.0f, 0.0f }; + int n{0}; + for ( int i {0}; i < iSelf; ++i){ + if ( glm::length(pos[i] - pos[iSelf]) < rule3Distance) { + vsum += vel[i]; + ++n; + } + } + for( int i {iSelf + 1} ; i < N; ++i){ + if ( glm::length(pos[i] - pos[iSelf]) < rule3Distance) { + vsum += vel[i]; + ++n; + } + } + if (n != 0) { + vsum *= scale / n; + } + return vsum; +} +__device__ glm::vec3 rule3Scatter(int iSelf, int const * gridCellStartIndices, + int const * gridCellEndIndices, int const * particleArrayIndices, + glm::vec3 const * pos, glm::vec3 const * vel, + int const usedGridIndices[maxGridCells], + float scale) +{ + glm::vec3 vsum{ 0.0f }; + int n{0}; + for (int k{0}; usedGridIndices[k] != -1; ++k) + { + int gridC{ usedGridIndices[k] }; + for ( int i {gridCellStartIndices[gridC]}; i < gridCellEndIndices[gridC]; ++i){ + int neighbor = particleArrayIndices[i]; + if ( neighbor == iSelf) { + continue; + } + if ( glm::length(pos[neighbor] - pos[iSelf]) < rule3Distance) { + vsum += vel[neighbor]; + ++n; + } + } + } + if (n != 0) { + vsum *= scale / n; + } + return vsum; +} +__device__ glm::vec3 rule3Coherent(int iSelf, int const * gridCellStartIndices, + int const * gridCellEndIndices, + glm::vec3 const * pos, glm::vec3 const * vel, + int const usedGridIndices[maxGridCells], + float scale) +{ + glm::vec3 vsum{ 0.0f }; + int n{0}; + for (int k{0}; usedGridIndices[k] != -1; ++k) + { + int gridC{ usedGridIndices[k] }; + for ( int neighbor {gridCellStartIndices[gridC]}; neighbor < gridCellEndIndices[gridC]; + ++neighbor){ + if ( neighbor == iSelf) { + continue; + } + if ( glm::length(pos[neighbor] - pos[iSelf]) < rule3Distance) { + vsum += vel[neighbor]; + ++n; + } + } + } + if (n != 0) { + vsum *= scale / n; + } + return vsum; +} + +// compute Velocity Change Scattered will get the velocity change with the uniform grid +__device__ glm::vec3 computeVelocityChangeScattered(int iSelf, const int * gridCellStartIndices, + const int * gridCellEndIndices, const int * particleArrayIndices, const glm::vec3 *pos, + const glm::vec3 *vel, int usedGridIndices[maxGridCells] ) +{ + // Rule 1: boids fly towards their local pe:rceived center of mass, which excludes themselves + // move towards the average + glm::vec3 v { rule1Scatter(iSelf, gridCellStartIndices, gridCellEndIndices, + particleArrayIndices, pos, usedGridIndices, rule1Scale)}; + //printf("%f %f %f\n", v.x, v.y, v.z); + // keep a boids apart + //Rule 2: boids try to stay a distance d away from each other + v += rule2Scatter(iSelf, gridCellStartIndices, gridCellEndIndices, + particleArrayIndices, pos, usedGridIndices, rule2Scale); + // Rule 3: boids try to match the speed of surrounding boids + v += rule3Scatter(iSelf, gridCellStartIndices, gridCellEndIndices, + particleArrayIndices, pos, vel, usedGridIndices, rule3Scale); + return v; +} +// compute Velocity Change Scattered will get the velocity change with the uniform grid +__device__ glm::vec3 computeVelocityChangeCoherent(int iSelf, const int * gridCellStartIndices, + const int * gridCellEndIndices, const glm::vec3 *pos, + const glm::vec3 *vel, int usedGridIndices[maxGridCells] ) +{ + // Rule 1: boids fly towards their local pe:rceived center of mass, which excludes themselves + // move towards the average + glm::vec3 v { rule1Coherent(iSelf, gridCellStartIndices, gridCellEndIndices, + pos, usedGridIndices, rule1Scale)}; + //printf("%f %f %f\n", v.x, v.y, v.z); + // keep a boids apart + //Rule 2: boids try to stay a distance d away from each other + v += rule2Coherent(iSelf, gridCellStartIndices, gridCellEndIndices, + pos, usedGridIndices, rule2Scale); + // Rule 3: boids try to match the speed of surrounding boids + v += rule3Coherent(iSelf, gridCellStartIndices, gridCellEndIndices, + pos, vel, usedGridIndices, rule3Scale); + return v; +} +__device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *pos, + const glm::vec3 *vel) +{ + // Rule 1: boids fly towards their local pe:rceived center of mass, which excludes themselves + // move towards the average + glm::vec3 v { rule1(N, iSelf, pos, rule1Scale)}; + // keep a boids apart + //Rule 2: boids try to stay a distance d away from each other + v += rule2(N, iSelf, pos, rule2Scale); // Rule 3: boids try to match the speed of surrounding boids - return glm::vec3(0.0f, 0.0f, 0.0f); + v += rule3(N, iSelf, pos, vel, rule3Scale); + return v; } /** @@ -242,7 +548,16 @@ __device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *po */ __global__ void kernUpdateVelocityBruteForce(int N, glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) { - // Compute a new velocity based on pos and vel1 + int boid = threadIdx.x + blockIdx.x * blockDim.x; + if (boid >= N) { + return; + } + glm::vec3 v3{ computeVelocityChange(N, boid, pos, vel1)}; + vel2[boid] = vel1[boid] + v3; + float length = glm::length(vel2[boid]); + if (length > maxSpeed){ + vel2[boid] = maxSpeed/length * vel2[boid]; + } // Clamp the speed // Record the new velocity into vel2. Question: why NOT vel1? } @@ -281,18 +596,41 @@ __global__ void kernUpdatePos(int N, float dt, glm::vec3 *pos, glm::vec3 *vel) { __device__ int gridIndex3Dto1D(int x, int y, int z, int gridResolution) { return x + y * gridResolution + z * gridResolution * gridResolution; } - +// returns true if the grid location is in range +__device__ bool gridInRange(int3 val, int gridResolution) { + return (val.x >= 0 && val.x <= gridResolution && + val.y >= 0 && val.y <= gridResolution && + val.z >= 0 && val.z <= gridResolution); +} + +__device__ int3 boidGrid(const glm::vec3 location, const glm::vec3 gridMin, + float inverseCellWidth) +{ + glm::vec3 boidCell = (location - gridMin) * inverseCellWidth; + return make_int3( (int) boidCell.x, (int) boidCell.y, (int) boidCell.z); +} +/// compute the grid and index arrays __global__ void kernComputeIndices(int N, int gridResolution, glm::vec3 gridMin, float inverseCellWidth, glm::vec3 *pos, int *indices, int *gridIndices) { + + int boid = (blockIdx.x * blockDim.x) + threadIdx.x; + if (boid >= N) { + return; + } + int3 gridCoord = boidGrid(pos[boid] , gridMin, inverseCellWidth); + int gridIndex = gridIndex3Dto1D(gridCoord.x, gridCoord.y, gridCoord.z, gridResolution); + gridIndices[boid] = gridIndex; + indices[boid] = boid; + + // 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 +// : 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) { @@ -300,12 +638,83 @@ __global__ void kernResetIntBuffer(int N, int *intBuffer, int value) { } } + +// findGridCells will find the neigboring grid cells relative to +// a boid. usedGridCells is the array of 9 gridCells that may hold +// the neighboring gridCells +// the idea here is that boidoffset is longer than the max nearest neighbor +// distance but less than a cellWidth. the two extremes of a 8 grid cell +// cube are given here and all neighbors have to be within these 8 cells. +__device__ void findGridCells(int iSelf, int gridResolution, + const glm::vec3 gridMin, float inverseCellWidth, + float cellWidth,const int *gridCellStartIndices, + const glm::vec3 *pos, int usedGridCells[maxGridCells]) +{ + // includes the grid cells that are used here. last indicates end of array + int xoffsets[4] = { 0, 1, 0, 1}; + int yoffsets[4] = { 0, 0, 1, 1}; + int currLoc = 0; + glm::vec3 boidOffset = glm::vec3(cellWidth * 0.5f, cellWidth * 0.5f, cellWidth * 0.5f); + glm::vec3 extremeBoid = pos[iSelf] + boidOffset; + int3 cellindices = boidGrid( extremeBoid, gridMin, inverseCellWidth); + for (int i = 0; i < 4; ++i) + { + int3 cpy = make_int3(cellindices.x - xoffsets[i], + cellindices.y - yoffsets[i], cellindices.z); + int gridCell = gridIndex3Dto1D(cpy.x, cpy.y, cpy.z, gridResolution); + if (gridInRange(cpy, gridResolution) && + gridCellStartIndices[gridCell] != -1) { + usedGridCells[currLoc++] = gridCell; + } + } + extremeBoid = pos[iSelf] - boidOffset; + cellindices = boidGrid( extremeBoid, gridMin, inverseCellWidth); + for (int i = 0; i < 4; ++i) + { + int3 cpy = make_int3(cellindices.x + xoffsets[i], + cellindices.y + yoffsets[i], cellindices.z); + int gridCell = gridIndex3Dto1D(cpy.x, cpy.y, cpy.z, gridResolution); + if (gridInRange(cpy, gridResolution) && + gridCellStartIndices[gridCell] != -1) { + usedGridCells[currLoc++] = gridCell; + } + } + usedGridCells[currLoc] = -1; +} + __global__ void kernIdentifyCellStartEnd(int N, int *particleGridIndices, - int *gridCellStartIndices, int *gridCellEndIndices) { + 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!" + // indicates the last cell that was found + // use the parallel structure to unroll always comparing with the prior + // The last grid occupied must be terminated. + // the ordering is [gridCellStartIndex gridCellEndIndex) (inclusive, exclusive). + int index = blockIdx.x * blockDim.x + threadIdx.x; + if (index >= N) { + return; + } + int currentGridLoc = particleGridIndices[index]; + if ( index == 0) { + gridCellStartIndices[currentGridLoc] = index; + return; + } + int priorGridLoc = particleGridIndices[index - 1]; + // no change indicated here + if ( priorGridLoc == currentGridLoc) { + return; + } + // there is a change compared to the prior--update Gridindices + gridCellEndIndices[priorGridLoc] = index; + gridCellStartIndices[currentGridLoc] = index; + // add the end of the currentNow + if (index == N - 1) { + gridCellEndIndices[currentGridLoc] = N; + } } __global__ void kernUpdateVelNeighborSearchScattered( @@ -314,6 +723,7 @@ __global__ void kernUpdateVelNeighborSearchScattered( int *gridCellStartIndices, int *gridCellEndIndices, int *particleArrayIndices, glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) { + // TODO-2.1 - Update a boid's velocity using the uniform grid to reduce // the number of boids that need to be checked. // - Identify the grid cell that this particle is in @@ -322,8 +732,43 @@ __global__ void kernUpdateVelNeighborSearchScattered( // - 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 iSelf = blockIdx.x * blockDim.x + threadIdx.x; + if (iSelf >= N) { + return; + } + int usedGridCells[maxGridCells]; + iSelf = particleArrayIndices[iSelf]; + findGridCells(iSelf, gridResolution, gridMin, inverseCellWidth, + cellWidth, gridCellStartIndices, pos, usedGridCells); + + glm::vec3 vchange {computeVelocityChangeScattered(iSelf, gridCellStartIndices, + gridCellEndIndices, particleArrayIndices, pos, + vel1, usedGridCells) }; + vel2[iSelf] = vel1[iSelf] + vchange; + float length = glm::length(vel2[iSelf]); + if (length > maxSpeed){ + vel2[iSelf] = maxSpeed/length * vel2[iSelf]; + } } +// kernMakeCoherent will take the particleArrayIndices and rewrite the glm::vec3 array +// (could be position or velocity) and put it into the to location based on the pointers in +// the particleArrayIndices. if particleArrayIndex[10] says 5 that means that +// posNew[10] = posOld[5]; +__global__ void kernMakeCoherent(int N, int const * particleArrayIndices, + glm::vec3 * to1, glm::vec3 const * from1, glm::vec3 * to2, + glm::vec3 const * from2) +{ + int index = threadIdx.x + blockIdx.x * blockDim.x; + if (index >= N) + { + return; + } + int oldLocation = particleArrayIndices[index]; + to1[index] = from1[oldLocation]; + to2[index] = from2[oldLocation]; +} +/// from vel1 to vel2 where vel2 is the new one. __global__ void kernUpdateVelNeighborSearchCoherent( int N, int gridResolution, glm::vec3 gridMin, float inverseCellWidth, float cellWidth, @@ -341,17 +786,45 @@ __global__ void kernUpdateVelNeighborSearchCoherent( // - Access each boid in the cell and compute velocity change from // the boids rules, if this boid is within the neighborhood distance. // - Clamp the speed change before putting the new speed in vel2 + int iSelf = blockIdx.x * blockDim.x + threadIdx.x; + if (iSelf >= N) { + return; + } + int usedGridCells[maxGridCells]; + findGridCells(iSelf, gridResolution, gridMin, inverseCellWidth, + cellWidth, gridCellStartIndices, pos, usedGridCells); + + glm::vec3 vchange {computeVelocityChangeCoherent(iSelf, gridCellStartIndices, + gridCellEndIndices, pos, + vel1, usedGridCells) }; + vel2[iSelf] = vel1[iSelf] + vchange; + float length = glm::length(vel2[iSelf]); + if (length > maxSpeed){ + vel2[iSelf] = maxSpeed/length * vel2[iSelf]; + } } /** * Step the entire N-body simulation by `dt` seconds. */ void Boids::stepSimulationNaive(float dt) { + + dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); + kernUpdateVelocityBruteForce<<>>(numObjects, dev_pos, + dev_vel1, dev_vel2); + checkCUDAErrorWithLine("brute force failed"); + kernUpdatePos<<>>(numObjects, dt, dev_pos, dev_vel2); + checkCUDAErrorWithLine("update Position Function Failed"); + + + std::swap(dev_vel1, dev_vel2); + cudaDeviceSynchronize(); // TODO-1.2 - use the kernels you wrote to step the simulation forward in time. // TODO-1.2 ping-pong the velocity buffers } void Boids::stepSimulationScatteredGrid(float dt) { + // TODO-2.1 // Uniform Grid Neighbor search using Thrust sort. // In Parallel: @@ -364,6 +837,40 @@ void Boids::stepSimulationScatteredGrid(float dt) { // - Perform velocity updates using neighbor search // - Update positions // - Ping-pong buffers as needed + + // compute the grid index of each boid + dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); + kernComputeIndices<<>> + (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, + dev_pos, dev_particleArrayIndices, dev_particleGridIndices); + + dim3 gridBlocksPerGrid((gridCellCount + blockSize -1) / blockSize); + // initialize to -1 all the gridCells + kernResetIntBuffer<<>>(gridCellCount, + dev_gridCellStartIndices, -1); + kernResetIntBuffer<<>>(gridCellCount, + dev_gridCellEndIndices, -1); + //Boids::testGridArrays("Before sorting Grid Cells"); + thrust::sort_by_key(dev_thrust_particleGridIndices, dev_thrust_particleGridIndices + numObjects, + dev_thrust_particleArrayIndices); + //Boids::testGridArrays("AfterGrid Cell Sorting"); + // update the start and end indices + kernIdentifyCellStartEnd<<>>(numObjects, + dev_particleGridIndices, dev_gridCellStartIndices, + dev_gridCellEndIndices); + cudaDeviceSynchronize(); + //Boids::testGridArrays("After StartEnd Assignment"); + // update the velocity + kernUpdateVelNeighborSearchScattered<<>> + (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, + gridCellWidth, dev_gridCellStartIndices, dev_gridCellEndIndices, + dev_particleArrayIndices, dev_pos, dev_vel1, dev_vel2); + + // update the position + kernUpdatePos<<>>(numObjects, dt, dev_pos, dev_vel2); + std::swap(dev_vel1, dev_vel2); + cudaDeviceSynchronize(); + } void Boids::stepSimulationCoherentGrid(float dt) { @@ -382,6 +889,48 @@ 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); + + dim3 gridBlocksPerGrid((gridCellCount + blockSize -1) / blockSize); + // initialize to -1 all the gridCells + kernResetIntBuffer<<>>(gridCellCount, + dev_gridCellStartIndices, -1); + kernResetIntBuffer<<>>(gridCellCount, + dev_gridCellEndIndices, -1); + //Boids::testGridArrays("Before sorting Grid Cells"); + thrust::sort_by_key(dev_thrust_particleGridIndices, dev_thrust_particleGridIndices + numObjects, + dev_thrust_particleArrayIndices); + + // make the position and velocity coherent-- when a grid cell uses + // particle indices 2 to 4 that should also be the boid number; + // Boids::seeFloatArray("positions before Sort", numObjects, dev_pos); + //now dev_vel2 has the coherently ordered velocity + // dev_posCopy also is coherently ordered + kernMakeCoherent<<>>(numObjects, + dev_particleArrayIndices, dev_posCopy, dev_pos, + dev_vel2, dev_vel1); + // swap puts the coherent position back into dev_pos + std::swap(dev_posCopy, dev_pos); + //Boids::seeFloatArray("Positions after Sort", numObjects, dev_posCopy); + //Boids::testGridArrays("AfterGrid Cell Sorting"); + // update the start and end indices + kernIdentifyCellStartEnd<<>>(numObjects, + dev_particleGridIndices, dev_gridCellStartIndices, + dev_gridCellEndIndices); + //cudaDeviceSynchronize(); + // Boids::testGridArrays("After StartEnd Assignment"); + // update the velocity + // dev_vel1 now has the updated velocity + kernUpdateVelNeighborSearchCoherent<<>> + (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, + gridCellWidth, dev_gridCellStartIndices, dev_gridCellEndIndices, + dev_pos, dev_vel2, dev_vel1); + // update the position + kernUpdatePos<<>>(numObjects, dt, dev_pos, dev_vel1); + cudaDeviceSynchronize(); } void Boids::endSimulation() { @@ -390,6 +939,11 @@ void Boids::endSimulation() { 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_posCopy); } void Boids::unitTest() { @@ -457,3 +1011,92 @@ void Boids::unitTest() { checkCUDAErrorWithLine("cudaFree failed!"); return; } +void printIndexPair(std::string name, int const * arr1, int const * arr2, int N) +{ + + std::cout << name << std::endl; + for (int nleft {0}; nleft < N; nleft+=numPerRow) + { + for (int i{nleft}; i < N && i < nleft + numPerRow; ++i) + { + std::cout << std::setw(spacePerInt) << i; + } + std::cout << std::endl; + for (int i{nleft}; i < N && i < nleft + numPerRow; ++i) + { + std::cout << std::setw(spacePerInt) << arr1[i]; + } + std::cout << std::endl; + for (int i{nleft}; i < N && i < nleft + numPerRow; ++i) + { + std::cout << std::setw(spacePerInt) << arr2[i]; + } + std::cout << std::endl << std::endl; + } +} +void printFloatArray(std::string name, int N, glm::vec3 const * arr) +{ + + std::cout << name << std::endl; + std::setprecision(precision); + + for (int nleft {0}; nleft < N; nleft+=numPerRowVec3) + { + for (int i{nleft}; i < N && i < nleft + numPerRowVec3; ++i) + { + std::cout << std::setw(spacePerInt) << i << ':'; + std::cout << std::fixed << std::setprecision(precision) << + std::setw(spacePerFloat) << arr[i].x; + std::cout << std::fixed << std::setw(spacePerFloat) << + std::setprecision(precision) << arr[i].y; + std::cout << std::fixed << std::setprecision(precision) << + std::setw(spacePerFloat) << arr[i].z; + std::cout << ';'; + } + std::cout << std::endl; + } + std::cout << std::endl; +} + +void Boids::testGridArrays(std::string msg) +{ + + int *host_gridCellStartIndices { new int[gridCellCount]}; + int *host_gridCellEndIndices { new int[gridCellCount]}; + int *host_particleArrayIndices { new int[numObjects]}; + int *host_particleGridIndices { new int[numObjects]}; + + + // How to copy data back to the CPU side from the GPU + cudaMemcpy(host_gridCellStartIndices, dev_gridCellStartIndices, sizeof(int) * gridCellCount, + cudaMemcpyDeviceToHost); + cudaMemcpy(host_gridCellEndIndices, dev_gridCellEndIndices, sizeof(int) * gridCellCount, + cudaMemcpyDeviceToHost); + cudaMemcpy(host_particleArrayIndices, dev_particleArrayIndices, sizeof(int) * numObjects, + cudaMemcpyDeviceToHost); + cudaMemcpy(host_particleGridIndices, dev_particleGridIndices, sizeof(int) * numObjects, + cudaMemcpyDeviceToHost); + checkCUDAErrorWithLine("memcpy back of grid arrays failed!"); + std::cout << msg << std::endl; + printIndexPair("GridIndex/ParticleIndex", host_particleGridIndices, host_particleArrayIndices, + numObjects); + printIndexPair("GridStartEndArray", host_gridCellStartIndices, host_gridCellEndIndices, + gridCellCount); + delete[] host_gridCellStartIndices; + delete[] host_gridCellEndIndices; + delete[] host_particleArrayIndices; + delete[] host_particleGridIndices; +} +void Boids::seeFloatArray(std::string msg, int N, glm::vec3 const * dev_array) +{ + + glm::vec3 *host_array { new glm::vec3[N]}; + + + // How to copy data back to the CPU side from the GPU + cudaMemcpy(host_array, dev_array, sizeof(glm::vec3) * N, cudaMemcpyDeviceToHost); + checkCUDAErrorWithLine("memcpy back of grid arrays failed!"); + std::cout << msg << std::endl; + printFloatArray("float Array", N, host_array); + delete[] host_array; +} diff --git a/src/kernel.h b/src/kernel.h index 3d3da72..88a8459 100644 --- a/src/kernel.h +++ b/src/kernel.h @@ -10,12 +10,15 @@ #include namespace Boids { - void initSimulation(int N); - void stepSimulationNaive(float dt); - void stepSimulationScatteredGrid(float dt); - void stepSimulationCoherentGrid(float dt); - void copyBoidsToVBO(float *vbodptr_positions, float *vbodptr_velocities); + void initSimulation(int N); + void stepSimulationNaive(float dt); + void stepSimulationScatteredGrid(float dt); + void stepSimulationCoherentGrid(float dt); + void copyBoidsToVBO(float *vbodptr_positions, float *vbodptr_velocities); - void endSimulation(); - void unitTest(); + void endSimulation(); + void unitTest(); + void testGridArrays(std::string); + // see the values in a float array + void seeFloatArray(std::string, int N, glm::vec3 const * dev_array); } diff --git a/src/main.cpp b/src/main.cpp index a29471d..447e286 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -13,12 +13,12 @@ // ================ // 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 1 +#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 = 10000; const float DT = 0.2f; /** @@ -92,7 +92,8 @@ bool init(int argc, char **argv) { glfwSetKeyCallback(window, keyCallback); glfwSetCursorPosCallback(window, mousePositionCallback); glfwSetMouseButtonCallback(window, mouseButtonCallback); - + // turn off vsync to increase the frame rate. + glfwSwapInterval(false); glewExperimental = GL_TRUE; if (glewInit() != GLEW_OK) { return false; @@ -220,7 +221,7 @@ void initShaders(GLuint * program) { double timebase = 0; int frame = 0; - Boids::unitTest(); // LOOK-1.2 We run some basic example code to make sure + Boids::unitTest(); // LOOK-1.2 We run some basic example code to make sure // your CUDA development setup is ready to go. while (!glfwWindowShouldClose(window)) { diff --git a/src/main.hpp b/src/main.hpp index 6cdaa93..a99c0a2 100644 --- a/src/main.hpp +++ b/src/main.hpp @@ -36,8 +36,9 @@ const float fovy = (float) (PI / 4); const float zNear = 0.10f; const float zFar = 10.0f; // LOOK-1.2: for high DPI displays, you may want to double these settings. -int width = 1280; -int height = 720; +// doubled +int width = 2560; +int height = 1440; int pointSize = 2.0f; // For camera controls