diff --git a/INSTRUCTION.md b/INSTRUCTION.md index ab40a66..fbc88bb 100644 --- a/INSTRUCTION.md +++ b/INSTRUCTION.md @@ -63,7 +63,7 @@ In the Boids flocking simulation, particles representing birds or fish 1. cohesion - boids move towards the perceived center of mass of their neighbors 2. separation - boids avoid getting to close to their neighbors 3. alignment - boids generally try to move with the same direction and speed as -their neighbors + their neighbors These three rules specify a boid's velocity change in a timestep. At every timestep, a boid thus has to look at each of its neighboring boids @@ -127,6 +127,7 @@ function rule3(Boid boid) return perceived_velocity * rule3Scale end ``` + Based on [Conard Parker's notes](http://www.vergenet.net/~conrad/boids/pseudocode.html) with slight adaptations. For the purposes of an interesting simulation, we will say that two boids only influence each other according if they are within a certain **neighborhood distance** of each other. @@ -141,10 +142,13 @@ For an idea of how the simulation "should" look in 3D, **Please Note** that our pseudocode, our 2D implementation, and our reference code (from which we derived the parameters that ship with the basecode) differ from Conrad Parker's notes in Rule 3 - our references do not subtract the boid's own velocity from the perceived velocity: Our pseuodocode: + ``` return perceived_velocity * rule3Scale ``` + Conrad Parker's notes: + ``` RETURN (pvJ - bJ.velocity) / 8 ``` @@ -159,12 +163,12 @@ However, since the purpose of this assignment is to introduce you to CUDA, we re * `src/main.cpp`: Performs all of the CUDA/OpenGL setup and OpenGL visualization. + * `src/kernel.cu`: CUDA device functions, state, kernels, and CPU functions for kernel invocations. In place of a unit testing/sandbox framework, there is space in here for individually running your kernels and getting the output back from the GPU before running the actual simulation. PLEASE make use of this in Part 2 to individually test your kernels. - 1. Search the code for `TODO-1.2` and `LOOK-1.2`. * `src/kernel.cu`: Use what you learned in the first lectures to figure out how to resolve these X Part 1 TODOs. @@ -205,7 +209,7 @@ because: 1. We don't have resizeable arrays on the GPU 2. Naively parallelizing the iteration may lead to race conditions, where two -particles need to be written into the same bucket on the same clock cycle. + particles need to be written into the same bucket on the same clock cycle. Instead, we will construct the uniform grid by sorting. If we label each boid with an index representing its enclosing cell and then sort the list of @@ -227,13 +231,14 @@ homework, we will use the value/key sort built into **Thrust**. See `Boids::unitTest` in `kernel.cu` for an example of how to use this. Your uniform grid will probably look something like this in GPU memory: + - `dev_particleArrayIndices` - buffer containing a pointer for each boid to its -data in dev_pos and dev_vel1 and dev_vel2 + data in dev_pos and dev_vel1 and dev_vel2 - `dev_particleGridIndices` - buffer containing the grid index of each boid - `dev_gridCellStartIndices` - buffer containing a pointer for each cell to the -beginning of its data in `dev_particleArrayIndices` + beginning of its data in `dev_particleArrayIndices` - `dev_gridCellEndIndices` - buffer containing a pointer for each cell to the -end of its data in `dev_particleArrayIndices`. + end of its data in `dev_particleArrayIndices`. Here the term `pointer` when used with buffers is largely interchangeable with the term `index`, however, you will effectively be using array indices as @@ -280,9 +285,10 @@ metric, but adding your own `cudaTimer`s, etc., will allow you to do more fine-grained benchmarking of various parts of your code. REMEMBER: + * Do your performance testing in `Release` mode! * Turn off Vertical Sync in Nvidia Control Panel: -![Unlock FPS](images/UnlockFPS.png) + ![Unlock FPS](images/UnlockFPS.png) * Performance should always be measured relative to some baseline when possible. A GPU can make your program faster - but by how much? * If a change impacts performance, show a comparison. Describe your changes. @@ -292,6 +298,7 @@ REMEMBER: ### Questions There are two ways to measure performance: + * Disable visualization so that the framerate reported will be for the the simulation only, and not be limited to 60 fps. This way, the framerate reported in the window title will be useful. @@ -307,27 +314,26 @@ hypotheses and insights. **Answer these:** * For each implementation, how does changing the number of boids affect -performance? Why do you think this is? + performance? Why do you think this is? * For each implementation, how does changing the block count and block size -affect performance? Why do you think this is? + affect performance? Why do you think this is? * For the coherent uniform grid: did you experience any performance improvements -with the more coherent uniform grid? Was this the outcome you expected? -Why or why not? + with the more coherent uniform grid? Was this the outcome you expected? + Why or why not? * Did changing cell width and checking 27 vs 8 neighboring cells affect performance? -Why or why not? Be careful: it is insufficient (and possibly incorrect) to say -that 27-cell is slower simply because there are more cells to check! + Why or why not? Be careful: it is insufficient (and possibly incorrect) to say + that 27-cell is slower simply because there are more cells to check! **NOTE: Nsight performance analysis tools *cannot* presently be used on the lab computers, as they require administrative access.** If you do not have access to a CUDA-capable computer, the lab computers still allow you to do timing mesasurements! However, the tools are very useful for performance debugging. - ## Part 4: Write-up 1. Take a screenshot of the boids **and** use a gif tool like [licecap](http://www.cockos.com/licecap/) to record an animations of the boids with a fixed camera. -Put this at the top of your README.md. Take a look at [How to make an attractive -GitHub repo](https://github.com/pjcozzi/Articles/blob/master/CIS565/GitHubRepo/README.md). + Put this at the top of your README.md. Take a look at [How to make an attractive + GitHub repo](https://github.com/pjcozzi/Articles/blob/master/CIS565/GitHubRepo/README.md). 2. Add your performance analysis. Graphs to include: - Framerate change with increasing # of boids for naive, scattered uniform grid, and coherent uniform grid (with and without visualization) - Framerate change with increasing block size @@ -343,9 +349,9 @@ The template of the comment section of your pull request is attached below, you * [Repo Link](https://link-to-your-repo) * (Briefly) Mentions features that you've completed. Especially those bells and whistles you want to highlight - * Feature 0 - * Feature 1 - * ... + * Feature 0 + * Feature 1 + * ... * Feedback on the project itself, if any. And you're done! @@ -353,21 +359,21 @@ And you're done! ## Tips - If your simulation crashes before launch, use -`checkCUDAErrorWithLine("message")` after CUDA invocations + `checkCUDAErrorWithLine("message")` after CUDA invocations - `ctrl + f5` in Visual Studio will launch the program but won't let the window -close if the program crashes. This way you can see any `checkCUDAErrorWithLine` -output. + close if the program crashes. This way you can see any `checkCUDAErrorWithLine` + output. - For debugging purposes, you can transfer data to and from the GPU. -See `Boids::unitTest` in `kernel.cu` for an example of how to use this. + See `Boids::unitTest` in `kernel.cu` for an example of how to use this. - For high DPI displays like 4K monitors or the Macbook Pro with Retina Display, you might want to double the rendering resolution and point size. See `main.hpp`. - Your README.md will be done in github markdown. You can find a [cheatsheet here](https://guides.github.com/pdfs/markdown-cheatsheet-online.pdf). There is -also a [live preview plugin](https://atom.io/packages/markdown-preview) for the -[atom text editor](https://atom.io/) from github. The same for [VS Code](https://www.visualstudio.com/en-us/products/code-vs.aspx) + also a [live preview plugin](https://atom.io/packages/markdown-preview) for the + [atom text editor](https://atom.io/) from github. The same for [VS Code](https://www.visualstudio.com/en-us/products/code-vs.aspx) - If your framerate is capped at 60fps, [disable V-sync](http://support.enmasse.com/tera/enable-v-sync-to-fix-graphics-issues-screen-tearing) ## Optional Extra Credit * Shared-Memory Optimization: - * Add fast nearest neighbor search using shared memory and the uniform grid. Include additional graphs and performance analysis, showing clearly how much better the program performed using shared memory. + * Add fast nearest neighbor search using shared memory and the uniform grid. Include additional graphs and performance analysis, showing clearly how much better the program performed using shared memory. * Grid-Looping Optimization: - * Instead of hard-coding a search of the designated area, limit the search area based on the grid cells that have any aspect of them within the max_distance. This prevents the excessive positional comparisons with the corner points of each grid cell, while at the same time also allowing a more flexible approach (since we're just defining a min cell index and max cell index in all three cardinal directions). That is, there is no longer a manual check for a hard-coded specific number of surrounding cells depending on the implementation (such as the 8 surrounding cells, 27 surrounding cells, etc). + * Instead of hard-coding a search of the designated area, limit the search area based on the grid cells that have any aspect of them within the max_distance. This prevents the excessive positional comparisons with the corner points of each grid cell, while at the same time also allowing a more flexible approach (since we're just defining a min cell index and max cell index in all three cardinal directions). That is, there is no longer a manual check for a hard-coded specific number of surrounding cells depending on the implementation (such as the 8 surrounding cells, 27 surrounding cells, etc). diff --git a/README.md b/README.md index ee39093..89d159e 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,39 @@ **University of Pennsylvania, CIS 5650: GPU Programming and Architecture, Project 1 - Flocking** -* (TODO) YOUR NAME HERE - * (TODO) [LinkedIn](), [personal website](), [twitter](), etc. -* Tested on: (TODO) Windows 22, i7-2222 @ 2.22GHz 22GB, GTX 222 222MB (Moore 2222 Lab) +* Zhanbo Lin + * [LinkedIn](https://www.linkedin.com/in/zhanbo-lin) +* Tested on: Windows 10, i5-10400F @ 2.90GHz 48GB, RTX-3080 10GB (Personal) +* GPU Compute Capability: 8.6 -### (TODO: Your README) +### Boids Simulation -Include screenshots, analysis, etc. (Remember, this is public, so don't put -anything here that you don't want to share with the world.) +![](images/simulation.gif) + +boids count = 100000 + +### Performance Analysis + +![](images/FPS-Boids.png) + +Figure-1 + +![](images/FPS-BlockSize.png) + +Figure-2 + +##### For each implementation, how does changing the number of boids affect performance? Why do you think this is? + +The framerate drops as the number of boids increases (Figure-1). This is expected and intuitive, that it takes more compute power to simulate more boids. + +##### For each implementation, how does changing the block count and block size affect performance? Why do you think this is? + +How the block count and block size affects performance depends on the actual implementation (Figure-2). The search algorithm reaches peak performance at blockSize=256. I think it is more of a hardware-dependent thing, that this blockSize aligns better with number of cores in each SM, resulting in more efficient scheduling. + +##### For the coherent uniform grid: did you experience any performance improvements with the more coherent uniform grid? Was this the outcome you expected? Why or why not? + +With the number of boids less than 100000, the performance difference varies between 8-neighbors and 27-neighbors search algorithm (Figure-1). For the 8-neighbors search, the coherent version performs better. In contrast, for the 27-neighbors serch, the scattered version performs better. This is expected because the coherent version has better data locality, which reduces cache misses. + +##### Did changing cell width and checking 27 vs 8 neighboring cells affect performance? Why or why not? + +The 27-neighbor version performs better for the scattered version but not the coherent version (Figure-1). I suspect this is related to data locality, possibly because the 27-neighbor search results in fewer boids per cell, though I am not certain of the exact cause.” diff --git a/images/FPS-BlockSize.png b/images/FPS-BlockSize.png new file mode 100644 index 0000000..bdd6462 Binary files /dev/null and b/images/FPS-BlockSize.png differ diff --git a/images/FPS-Boids.png b/images/FPS-Boids.png new file mode 100644 index 0000000..7338ae4 Binary files /dev/null and b/images/FPS-Boids.png differ diff --git a/images/simulation.gif b/images/simulation.gif new file mode 100644 index 0000000..8f17941 Binary files /dev/null and b/images/simulation.gif differ diff --git a/src/kernel.cu b/src/kernel.cu index 7149917..ac054d6 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -30,15 +30,18 @@ /** * Check for CUDA errors; print and exit if there was a problem. */ -void checkCUDAError(const char *msg, int line = -1) { - cudaError_t err = cudaGetLastError(); - if (cudaSuccess != err) { - if (line >= 0) { - fprintf(stderr, "Line %d: ", line); +void checkCUDAError(const char* msg, int line = -1) +{ + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err) + { + if (line >= 0) + { + fprintf(stderr, "Line %d: ", line); + } + fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString(err)); + exit(EXIT_FAILURE); } - fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString(err)); - exit(EXIT_FAILURE); - } } @@ -64,6 +67,7 @@ void checkCUDAError(const char *msg, int line = -1) { /*! Size of the starting area in simulation space. */ #define scene_scale 100.0f +#define DOUBLE_DISTANCE 0 /*********************************************** * Kernel state (pointers are device pointers) * ***********************************************/ @@ -76,22 +80,26 @@ 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_vel1; -glm::vec3 *dev_vel2; +glm::vec3* dev_pos; +glm::vec3* dev_vel1; +glm::vec3* dev_vel2; // LOOK-2.1 - these are NOT allocated for you. You'll have to set up the thrust // pointers on your own too. // For efficient sorting and the uniform grid. These should always be parallel. -int *dev_particleArrayIndices; // What index in dev_pos and dev_velX represents this particle? -int *dev_particleGridIndices; // What grid cell is this particle in? +int* dev_particleArrayIndices; // What index in dev_pos and dev_velX represents this particle? +int* dev_particleGridIndices; // What grid cell is this particle in? // needed for use with thrust thrust::device_ptr dev_thrust_particleArrayIndices; thrust::device_ptr dev_thrust_particleGridIndices; +thrust::device_ptr dev_thrust_pos; +thrust::device_ptr dev_thrust_vel1; +thrust::device_ptr dev_thrust_vel2; -int *dev_gridCellStartIndices; // What part of dev_particleArrayIndices belongs -int *dev_gridCellEndIndices; // to this cell? + +int* dev_gridCellStartIndices; // What part of dev_particleArrayIndices belongs +int* dev_gridCellEndIndices; // to this cell? // TODO-2.3 - consider what additional buffers you might need to reshuffle // the position and velocity data to be coherent within cells. @@ -104,82 +112,128 @@ float gridCellWidth; float gridInverseCellWidth; glm::vec3 gridMinimum; +glm::ivec3 gridCoordMinumum; + +/********************* + * Helper Functions * + ********************/ +__device__ glm::ivec3 getGridCoord(glm::vec3 pos, float invGridWidth) +{ + glm::ivec3 grid = glm::floor(pos * invGridWidth); + return grid; +} + +__device__ glm::ivec3 getBottomLeftGridCoord(glm::vec3 pos, float halfGridWidth, float invGridWidth) +{ + glm::ivec3 bottomLeft = glm::floor((pos - halfGridWidth) * invGridWidth); + return bottomLeft; +} + + + /****************** * initSimulation * ******************/ -__host__ __device__ unsigned int hash(unsigned int a) { - a = (a + 0x7ed55d16) + (a << 12); - a = (a ^ 0xc761c23c) ^ (a >> 19); - a = (a + 0x165667b1) + (a << 5); - a = (a + 0xd3a2646c) ^ (a << 9); - a = (a + 0xfd7046c5) + (a << 3); - a = (a ^ 0xb55a4f09) ^ (a >> 16); - return a; +__host__ __device__ unsigned int hash(unsigned int a) +{ + a = (a + 0x7ed55d16) + (a << 12); + a = (a ^ 0xc761c23c) ^ (a >> 19); + a = (a + 0x165667b1) + (a << 5); + a = (a + 0xd3a2646c) ^ (a << 9); + a = (a + 0xfd7046c5) + (a << 3); + a = (a ^ 0xb55a4f09) ^ (a >> 16); + return a; } /** * LOOK-1.2 - this is a typical helper function for a CUDA kernel. * Function for generating a random vec3. */ -__host__ __device__ glm::vec3 generateRandomVec3(float time, int index) { - thrust::default_random_engine rng(hash((int)(index * time))); - thrust::uniform_real_distribution unitDistrib(-1, 1); +__host__ __device__ glm::vec3 generateRandomVec3(float time, int index) +{ + thrust::default_random_engine rng(hash((int)(index * time))); + thrust::uniform_real_distribution unitDistrib(-1, 1); - return glm::vec3((float)unitDistrib(rng), (float)unitDistrib(rng), (float)unitDistrib(rng)); + 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. */ -__global__ void kernGenerateRandomPosArray(int time, int N, glm::vec3 * arr, float scale) { - int index = (blockIdx.x * blockDim.x) + threadIdx.x; - if (index < N) { - glm::vec3 rand = generateRandomVec3(time, index); - arr[index].x = scale * rand.x; - arr[index].y = scale * rand.y; - arr[index].z = scale * rand.z; - } +__global__ void kernGenerateRandomPosArray(int time, int N, glm::vec3* arr, float scale) +{ + int index = (blockIdx.x * blockDim.x) + threadIdx.x; + if (index < N) + { + glm::vec3 rand = generateRandomVec3(time, index); + arr[index].x = scale * rand.x; + arr[index].y = scale * rand.y; + arr[index].z = scale * rand.z; + } } /** * Initialize memory, update some globals */ -void Boids::initSimulation(int N) { - numObjects = N; - dim3 fullBlocksPerGrid((N + blockSize - 1) / blockSize); - - // LOOK-1.2 - This is basic CUDA memory management and error checking. - // 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_vel2, N * sizeof(glm::vec3)); - checkCUDAErrorWithLine("cudaMalloc dev_vel2 failed!"); - - // LOOK-1.2 - This is a typical CUDA kernel invocation. - kernGenerateRandomPosArray<<>>(1, numObjects, - dev_pos, scene_scale); - checkCUDAErrorWithLine("kernGenerateRandomPosArray failed!"); - - // LOOK-2.1 computing grid params - gridCellWidth = 2.0f * std::max(std::max(rule1Distance, rule2Distance), rule3Distance); - int halfSideCount = (int)(scene_scale / gridCellWidth) + 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. - cudaDeviceSynchronize(); +void Boids::initSimulation(int N) +{ + numObjects = N; + dim3 fullBlocksPerGrid((N + blockSize - 1) / blockSize); + + // LOOK-1.2 - This is basic CUDA memory management and error checking. + // 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_vel2, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_vel2 failed!"); + + // LOOK-1.2 - This is a typical CUDA kernel invocation. + kernGenerateRandomPosArray << > > (1, numObjects, + dev_pos, scene_scale); + checkCUDAErrorWithLine("kernGenerateRandomPosArray failed!"); + + // LOOK-2.1 computing grid params +#if DOUBLE_DISTANCE + gridCellWidth = 2.0f * std::max(std::max(rule1Distance, rule2Distance), rule3Distance); +#else + gridCellWidth = std::max(std::max(rule1Distance, rule2Distance), rule3Distance); +#endif + int halfSideCount = (int)(scene_scale / gridCellWidth) + 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. + // 2.1 buffers + cudaMalloc((void**)&dev_particleArrayIndices, numObjects * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_particleArrayIndices failed"); + cudaMalloc((void**)&dev_particleGridIndices, numObjects * 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"); + + dev_thrust_particleArrayIndices = thrust::device_ptr(dev_particleArrayIndices); + dev_thrust_particleGridIndices = thrust::device_ptr(dev_particleGridIndices); + dev_thrust_pos = thrust::device_ptr(dev_pos); + dev_thrust_vel1 = thrust::device_ptr(dev_vel1); + dev_thrust_vel2 = thrust::device_ptr(dev_vel2); + checkCUDAErrorWithLine("assign thrust device pointers failed"); + + + cudaDeviceSynchronize(); } @@ -190,42 +244,47 @@ void Boids::initSimulation(int N) { /** * Copy the boid positions into the VBO so that they can be drawn by OpenGL. */ -__global__ void kernCopyPositionsToVBO(int N, glm::vec3 *pos, float *vbo, float s_scale) { - int index = threadIdx.x + (blockIdx.x * blockDim.x); - - float c_scale = -1.0f / s_scale; - - if (index < N) { - vbo[4 * index + 0] = pos[index].x * c_scale; - vbo[4 * index + 1] = pos[index].y * c_scale; - vbo[4 * index + 2] = pos[index].z * c_scale; - vbo[4 * index + 3] = 1.0f; - } +__global__ void kernCopyPositionsToVBO(int N, glm::vec3* pos, float* vbo, float s_scale) +{ + int index = threadIdx.x + (blockIdx.x * blockDim.x); + + float c_scale = -1.0f / s_scale; + + if (index < N) + { + vbo[4 * index + 0] = pos[index].x * c_scale; + vbo[4 * index + 1] = pos[index].y * c_scale; + vbo[4 * index + 2] = pos[index].z * c_scale; + vbo[4 * index + 3] = 1.0f; + } } -__global__ void kernCopyVelocitiesToVBO(int N, glm::vec3 *vel, float *vbo, float s_scale) { - int index = threadIdx.x + (blockIdx.x * blockDim.x); +__global__ void kernCopyVelocitiesToVBO(int N, glm::vec3* vel, float* vbo, float s_scale) +{ + int index = threadIdx.x + (blockIdx.x * blockDim.x); - if (index < N) { - vbo[4 * index + 0] = vel[index].x + 0.3f; - vbo[4 * index + 1] = vel[index].y + 0.3f; - vbo[4 * index + 2] = vel[index].z + 0.3f; - vbo[4 * index + 3] = 1.0f; - } + if (index < N) + { + vbo[4 * index + 0] = vel[index].x + 0.3f; + vbo[4 * index + 1] = vel[index].y + 0.3f; + vbo[4 * index + 2] = vel[index].z + 0.3f; + vbo[4 * index + 3] = 1.0f; + } } /** * Wrapper for call to the kernCopyboidsToVBO CUDA kernel. */ -void Boids::copyBoidsToVBO(float *vbodptr_positions, float *vbodptr_velocities) { - dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); +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); + kernCopyPositionsToVBO << > > (numObjects, dev_pos, vbodptr_positions, scene_scale); + kernCopyVelocitiesToVBO << > > (numObjects, dev_vel1, vbodptr_velocities, scene_scale); - checkCUDAErrorWithLine("copyBoidsToVBO failed!"); + checkCUDAErrorWithLine("copyBoidsToVBO failed!"); - cudaDeviceSynchronize(); + cudaDeviceSynchronize(); } @@ -239,47 +298,274 @@ 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 - // Rule 3: boids try to match the speed of surrounding boids - return glm::vec3(0.0f, 0.0f, 0.0f); +__device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3* pos, const glm::vec3* vel) +{ + using vec3 = glm::vec3; + + vec3 selfPos = pos[iSelf]; + vec3 selfVelocity = vel[iSelf]; + + // Rule 1 + vec3 center {0}; + int r1Count = 0; + + // Rule 2 + vec3 c{0}; + + // Rule 3 + vec3 v{ 0 }; + int r3Count = 0; + + for (int i = 0; i < N; i++) + { + vec3 bPos = pos[i]; + vec3 bVel = vel[i]; + float dist = glm::length(bPos - selfPos); + + if (i == iSelf) { continue; } + + // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves + if (dist < rule1Distance) + { + center += bPos; + r1Count++; + } + + // Rule 2: boids try to stay a distance d away from each other + if (dist < rule2Distance) + { + c -= (bPos - selfPos); + } + + // Rule 3: boids try to match the speed of surrounding boids + if (dist < rule3Distance) + { + v += bVel; + r3Count++; + } + } + + // Rule 1 + vec3 v1 {0}; + if (r1Count) + { + center /= imax(r1Count, 1); + v1 = (center - selfPos) * rule1Scale; + } + + // Rule 2 + vec3 v2 = c * rule2Scale; + + // Rule 3 + vec3 v3 {0}; + if (r3Count) + { + v /= imax(r3Count, 1); + v3 = v * rule3Scale; + } + + return selfVelocity + v1 + v2 + v3; +} + +__device__ glm::vec3 computeVelocityChangeGrid(int N, int* neighbors, int iSelf, const glm::vec3* pos, const glm::vec3* vel, + const int* particleArrayIndices, const int* gridCellStartIndices, const int* gridCellEndIndices) +{ + using glm::vec3; + + vec3 selfPos = pos[iSelf]; + vec3 selfVelocity = vel[iSelf]; + + // Rule 1 + vec3 center{ 0 }; + int r1Count = 0; + + // Rule 2 + vec3 c{ 0 }; + + // Rule 3 + vec3 v{ 0 }; + int r3Count = 0; + + for (int neighborIdx = 0; neighborIdx < N; neighborIdx++) + { + int cellID = neighbors[neighborIdx]; + int start = gridCellStartIndices[cellID]; + int end = gridCellEndIndices[cellID]; + for (int bIdx = start; bIdx < end; bIdx++) + { + int bID = particleArrayIndices[bIdx]; + vec3 bPos = pos[bID]; + vec3 bVel = vel[bID]; + float dist = glm::length(bPos - selfPos); + + if (bID == iSelf) { continue; } + + // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves + if (dist < rule1Distance) + { + center += bPos; + r1Count++; + } + + // Rule 2: boids try to stay a distance d away from each other + if (dist < rule2Distance) + { + c -= (bPos - selfPos); + } + + // Rule 3: boids try to match the speed of surrounding boids + if (dist < rule3Distance) + { + v += bVel; + r3Count++; + } + } + } + + // Rule 1 + vec3 v1{ 0 }; + if (r1Count) + { + center /= imax(r1Count, 1); + v1 = (center - selfPos) * rule1Scale; + } + + // Rule 2 + vec3 v2 = c * rule2Scale; + + // Rule 3 + vec3 v3{ 0 }; + if (r3Count) + { + v /= imax(r3Count, 1); + v3 = v * rule3Scale; + } + + return selfVelocity + v1 + v2 + v3; +} + +__device__ glm::vec3 computeVelocityChangeGridCoherent(int N, int* neighbors, int iSelf, const glm::vec3* pos, const glm::vec3* vel, + const int* gridCellStartIndices, const int* gridCellEndIndices) +{ + using glm::vec3; + + vec3 selfPos = pos[iSelf]; + vec3 selfVelocity = vel[iSelf]; + + // Rule 1 + vec3 center{ 0 }; + int r1Count = 0; + + // Rule 2 + vec3 c{ 0 }; + + // Rule 3 + vec3 v{ 0 }; + int r3Count = 0; + + for (int neighborIdx = 0; neighborIdx < N; neighborIdx++) + { + int cellID = neighbors[neighborIdx]; + int start = gridCellStartIndices[cellID]; + int end = gridCellEndIndices[cellID]; + for (int bID = start; bID < end; bID++) + { + vec3 bPos = pos[bID]; + vec3 bVel = vel[bID]; + float dist = glm::length(bPos - selfPos); + + if (bID == iSelf) { continue; } + + // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves + if (dist < rule1Distance) + { + center += bPos; + r1Count++; + } + + // Rule 2: boids try to stay a distance d away from each other + if (dist < rule2Distance) + { + c -= (bPos - selfPos); + } + + // Rule 3: boids try to match the speed of surrounding boids + if (dist < rule3Distance) + { + v += bVel; + r3Count++; + } + } + } + + // Rule 1 + vec3 v1{ 0 }; + if (r1Count) + { + center /= imax(r1Count, 1); + v1 = (center - selfPos) * rule1Scale; + } + + // Rule 2 + vec3 v2 = c * rule2Scale; + + // Rule 3 + vec3 v3{ 0 }; + if (r3Count) + { + v /= imax(r3Count, 1); + v3 = v * rule3Scale; + } + + return selfVelocity + v1 + v2 + v3; } /** * 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? +__global__ void kernUpdateVelocityBruteForce(int N, glm::vec3* pos, + glm::vec3* vel1, glm::vec3* vel2) +{ + int idx = threadIdx.x + (blockIdx.x * blockDim.x); + if (idx >= N) { return; } + + // Compute a new velocity based on pos and vel1 + auto newVel = computeVelocityChange(N, idx, pos, vel1); + + // Clamp the speed + if (glm::length(newVel) > maxSpeed) + newVel = glm::normalize(newVel) * maxSpeed; + + // Record the new velocity into vel2. Question: why NOT vel1? + vel2[idx] = newVel; } /** * 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. */ -__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; +__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; - // 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. @@ -288,180 +574,436 @@ __global__ void kernUpdatePos(int N, float dt, glm::vec3 *pos, glm::vec3 *vel) { // 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; +__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) { + 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 idx = threadIdx.x + (blockIdx.x * blockDim.x); + if (idx >= N) { return; } + + // get grid index + glm::vec3 boidPos = pos[idx]; + glm::vec3 boidCoord = boidPos - gridMin; + glm::ivec3 gridCoord = getGridCoord(boidCoord, inverseCellWidth); + int gridCellIdx = gridIndex3Dto1D(gridCoord.x, gridCoord.y, gridCoord.z, gridResolution); + + indices[idx] = idx; + gridIndices[idx] = gridCellIdx; } // 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 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 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 idx = threadIdx.x + (blockIdx.x * blockDim.x); + if (idx >= N) { return; } + + int prevIdx = idx - 1; + int nextIdx = idx + 1; + int gridIdx = particleGridIndices[idx]; + + // Start + if (prevIdx < 0) + { + gridCellStartIndices[gridIdx] = idx; + } + else + { + int prevGridIdx = particleGridIndices[prevIdx]; + if (gridIdx != prevGridIdx) + { + gridCellStartIndices[gridIdx] = idx; + } + } + + // End + if (nextIdx >= N) + { + gridCellEndIndices[gridIdx] = nextIdx; + } + else + { + int nextGridIdx = particleGridIndices[nextIdx]; + if (gridIdx != nextGridIdx) + { + gridCellEndIndices[gridIdx] = nextIdx; + } + } } __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 + using glm::vec3; + using glm::ivec3; + + int selfIdx = threadIdx.x + (blockIdx.x * blockDim.x); + if (selfIdx >= N) { return; } + + vec3 selfPos = pos[selfIdx]; + + // get neighbor search locations + int neighborCount = 0; + +#if DOUBLE_DISTANCE + int neighborCell[8]; + ivec3 searchGridBottomLeft = getBottomLeftGridCoord(selfPos - gridMin, cellWidth * 0.5f, inverseCellWidth); + ivec3 searchGridTopRight = searchGridBottomLeft + ivec3(1); + + searchGridBottomLeft = glm::max(searchGridBottomLeft, ivec3(0)); + searchGridTopRight = glm::min(searchGridTopRight, ivec3(gridResolution - 1)); +#else + int neighborCell[27]; + ivec3 searchCellCenter = getGridCoord(selfPos - gridMin, inverseCellWidth); + ivec3 searchGridBottomLeft = searchCellCenter + ivec3(-1); + ivec3 searchGridTopRight = searchCellCenter + ivec3(1); + + searchGridBottomLeft = glm::max(searchGridBottomLeft, ivec3(0)); + searchGridTopRight = glm::min(searchGridTopRight, ivec3(gridResolution - 1)); +#endif + + for (int z = searchGridBottomLeft.z; z <= searchGridTopRight.z; z++) + { + for (int y = searchGridBottomLeft.y; y <= searchGridTopRight.y; y++) + { + for (int x = searchGridBottomLeft.x; x <= searchGridTopRight.x; x++) + { + int gridIdx = gridIndex3Dto1D(x, y, z, gridResolution); + if (gridCellStartIndices[gridIdx] == -1) + { + continue; + } + neighborCell[neighborCount++] = gridIdx; + } + } + } + + // compute new velocity from neighbor cells + vec3 newVel = computeVelocityChangeGrid(neighborCount, neighborCell, selfIdx, pos, vel1, particleArrayIndices, gridCellStartIndices, gridCellEndIndices); + if (glm::length(newVel) > maxSpeed) + { + newVel = glm::normalize(newVel) * maxSpeed; + } + vel2[selfIdx] = 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 + using glm::vec3; + using glm::ivec3; + + int selfIdx = threadIdx.x + (blockIdx.x * blockDim.x); + if (selfIdx >= N) { return; } + + vec3 selfPos = pos[selfIdx]; + + // get neighbor search locations + int neighborCount = 0; + +#if DOUBLE_DISTANCE + int neighborCell[8]; + ivec3 searchGridBottomLeft = getBottomLeftGridCoord(selfPos - gridMin, cellWidth * 0.5f, inverseCellWidth); + ivec3 searchGridTopRight = searchGridBottomLeft + ivec3(1); + + searchGridBottomLeft = glm::max(searchGridBottomLeft, ivec3(0)); + searchGridTopRight = glm::min(searchGridTopRight, ivec3(gridResolution - 1)); +#else + int neighborCell[27]; + ivec3 searchCellCenter = getGridCoord(selfPos - gridMin, inverseCellWidth); + ivec3 searchGridBottomLeft = searchCellCenter + ivec3(-1); + ivec3 searchGridTopRight = searchCellCenter + ivec3(1); + + searchGridBottomLeft = glm::max(searchGridBottomLeft, ivec3(0)); + searchGridTopRight = glm::min(searchGridTopRight, ivec3(gridResolution - 1)); +#endif + + for (int z = searchGridBottomLeft.z; z <= searchGridTopRight.z; z++) + { + for (int y = searchGridBottomLeft.y; y <= searchGridTopRight.y; y++) + { + for (int x = searchGridBottomLeft.x; x <= searchGridTopRight.x; x++) + { + int gridIdx = gridIndex3Dto1D(x, y, z, gridResolution); + if (gridCellStartIndices[gridIdx] == -1) + { + continue; + } + neighborCell[neighborCount++] = gridIdx; + } + } + } + + // compute new velocity from neighbor cells + vec3 newVel = computeVelocityChangeGridCoherent(neighborCount, neighborCell, selfIdx, pos, vel1, gridCellStartIndices, gridCellEndIndices); + if (glm::length(newVel) > maxSpeed) + { + newVel = glm::normalize(newVel) * maxSpeed; + } + vel2[selfIdx] = newVel; } + /** * 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 +void Boids::stepSimulationNaive(float dt) +{ + dim3 numBlocks(utilityCore::divup(numObjects, blockSize)); + + // TODO-1.2 - use the kernels you wrote to step the simulation forward in time. + kernUpdateVelocityBruteForce<<>>(numObjects, dev_pos, dev_vel1, dev_vel2); + kernUpdatePos<<>>(numObjects, dt, dev_pos, dev_vel2); + + // TODO-1.2 ping-pong the velocity buffers + 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 +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. + // ZB: thrust::stable_sort_by_key() + + // - 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 + + // set up + dim3 numBlocks(utilityCore::divup(numObjects, blockSize)); + dim3 gridNumBlocks(utilityCore::divup(gridCellCount, blockSize)); + + // identify the locating grids + kernComputeIndices <<< numBlocks, threadsPerBlock >>> (numObjects, gridSideCount, gridMinimum, + gridInverseCellWidth, dev_pos, dev_particleArrayIndices, dev_particleGridIndices); + + // sort by grid + thrust::sort_by_key(dev_thrust_particleGridIndices, dev_thrust_particleGridIndices + numObjects, + dev_thrust_particleArrayIndices); + + // identify starts and ends + kernResetIntBuffer <<< gridNumBlocks, threadsPerBlock >>> (gridCellCount, dev_gridCellStartIndices, -1); + kernResetIntBuffer <<< gridNumBlocks, threadsPerBlock >>> (gridCellCount, dev_gridCellEndIndices, -1); + kernIdentifyCellStartEnd <<< numBlocks, threadsPerBlock >>> (numObjects, dev_particleGridIndices, dev_gridCellStartIndices, dev_gridCellEndIndices); + + // update velocity and position + kernUpdateVelNeighborSearchScattered <<< numBlocks, threadsPerBlock >>> (numObjects, gridSideCount, gridMinimum, + gridInverseCellWidth, gridCellWidth, + dev_gridCellStartIndices, dev_gridCellEndIndices, + dev_particleArrayIndices, + dev_pos, dev_vel1, dev_vel2); + + kernUpdatePos <<< numBlocks, threadsPerBlock >>>(numObjects, dt, dev_pos, dev_vel2); + + // swap + 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. +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. + // set up + dim3 numBlocks(utilityCore::divup(numObjects, blockSize)); + dim3 gridNumBlocks(utilityCore::divup(gridCellCount, blockSize)); + + + // identify the locating grids + kernComputeIndices << < numBlocks, threadsPerBlock >> > (numObjects, gridSideCount, gridMinimum, + gridInverseCellWidth, dev_pos, dev_particleArrayIndices, dev_particleGridIndices); + + // sort by grid + dev_thrust_vel1 = thrust::device_ptr(dev_vel1); + auto zipItr = thrust::make_zip_iterator(dev_thrust_pos, dev_thrust_vel1); + + thrust::sort_by_key(dev_thrust_particleGridIndices, dev_thrust_particleGridIndices + numObjects, zipItr); + + // identify starts and ends + kernResetIntBuffer <<< gridNumBlocks, threadsPerBlock >>> (gridCellCount, dev_gridCellStartIndices, -1); + kernResetIntBuffer <<< gridNumBlocks, threadsPerBlock >>> (gridCellCount, dev_gridCellEndIndices, -1); + kernIdentifyCellStartEnd <<< numBlocks, threadsPerBlock >>> (numObjects, dev_particleGridIndices, dev_gridCellStartIndices, dev_gridCellEndIndices); + + // update velocity and position + kernUpdateVelNeighborSearchCoherent <<< numBlocks, threadsPerBlock >>> (numObjects, gridSideCount, gridMinimum, + gridInverseCellWidth, gridCellWidth, + dev_gridCellStartIndices, dev_gridCellEndIndices, + dev_pos, dev_vel1, dev_vel2); + + kernUpdatePos <<< numBlocks, threadsPerBlock >>> (numObjects, dt, dev_pos, dev_vel2); + + // swap + std::swap(dev_vel1, dev_vel2); } -void Boids::endSimulation() { - cudaFree(dev_vel1); - cudaFree(dev_vel2); - cudaFree(dev_pos); +void Boids::endSimulation() +{ + cudaFree(dev_vel1); + cudaFree(dev_vel2); + cudaFree(dev_pos); - // TODO-2.1 TODO-2.3 - Free any additional buffers here. + // TODO-2.1 TODO-2.3 - Free any additional buffers here. + cudaFree(dev_particleArrayIndices); + cudaFree(dev_particleGridIndices); + cudaFree(dev_gridCellStartIndices); + cudaFree(dev_gridCellEndIndices); + + checkCUDAErrorWithLine("cudaFree failed!"); } -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; +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; + } + + /*std::cout << "gridCellCount: " << gridCellCount << std::endl; + std::cout << "gridSideCount: " << gridSideCount << std::endl; + std::cout << "gridCellWidth: " << gridCellWidth << std::endl; + std::cout << "gridInverseCellWidth: " << gridInverseCellWidth << std::endl; + std::cout << "gridMinimum: " << gridMinimum.x << "," << gridMinimum.y << "," << gridMinimum.z << std::endl; + + for (int i = 0; i <= 4; i++) + { + for (int j = 0; j <= 4; j++) + { + for (int k = 0; k <= 4; k++) + { + glm::vec3 pos = gridMinimum; + pos += scene_scale * 2 * glm::vec3(i, j, k) / 4.0f; + + glm::ivec3 gridPos = glm::floor((pos - gridCellWidth*0.5f) * gridInverseCellWidth); + + std::cout << "pos: " << pos.x << "," << pos.y << "," << pos.z << " grid: " << gridPos.x << "," << gridPos.y << "," << gridPos.z << std::endl; + } + } + }*/ + // cleanup + cudaFree(dev_intKeys); + cudaFree(dev_intValues); + checkCUDAErrorWithLine("cudaFree failed!"); + return; } diff --git a/src/main.cpp b/src/main.cpp index 9c917c0..5713c58 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -22,12 +22,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 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 = 100000; // 58241;? const float DT = 0.2f; /** diff --git a/src/utilityCore.hpp b/src/utilityCore.hpp index fd8d313..484f7e6 100644 --- a/src/utilityCore.hpp +++ b/src/utilityCore.hpp @@ -33,6 +33,12 @@ extern void printCudaMat4(const cudaMat4 &m); extern std::string convertIntToString(int number); extern std::istream& safeGetline(std::istream& is, std::string& t); //Thanks to http://stackoverflow.com/a/6089413 + template + T divup(T size, T div) + { + return (size + div - 1) / div; + } + //----------------------------- //-------GLM Printers---------- //-----------------------------