diff --git a/INSTRUCTION.md b/INSTRUCTION.md deleted file mode 100644 index 86e3553..0000000 --- a/INSTRUCTION.md +++ /dev/null @@ -1,371 +0,0 @@ -CUDA Introduction - Flocking -============================ - -This is due Sunday, September 13. - -**Summary:** In this project, you will get some real experience writing simple -CUDA kernels, using them, and analyzing their performance. You'll implement a -flocking simulation based on the Reynolds Boids algorithm, along with two levels -of optimization: a uniform grid, and a uniform grid with semi-coherent memory access. - -## Part 0: Nothing New - -This project (and all other CUDA projects in this course) requires an NVIDIA -graphics card with CUDA capability. Any card with Compute Capability 3.0 -(GTX 6xx and newer) or greater will work. Check your GPU on this -[compatibility table](https://developer.nvidia.com/cuda-gpus). -If you do not have a personal machine with these specs, you may use those -computers in the Moore 100B/C which have supported GPUs, but pay close -attention to the **setup** section below. - -**HOWEVER**: If you need to use the lab computer for your development, you will -not presently be able to do GPU performance profiling. This will be very -important for debugging performance bottlenecks in your program. If you do not -have administrative access to any CUDA-capable machine, please email the TA. - -## Part 1: Naive Boids Simulation - -### 1.0. Setup - The Usual - -See Project 0, Parts 1-3 for reference. - -If you are using the Nsight IDE (not Visual Studio) and started Project 0 -early, note that things have -changed slightly. Instead of creating a new project, use -*File->Import->General->Existing Projects Into Workspace*, and select the -project folder as the root directory. Under *Project->Build -Configurations->Set Active...*, you can now select various Release and Debug -builds. - -* `src/` contains the source code. -* `external/` contains the binaries and headers for GLEW, GLFW, and GLM. - -**CMake note:** Do not change any build settings or add any files to your -project directly (in Visual Studio, Nsight, etc.) Instead, edit the -`CMakeLists.txt` file. Any files you create must be added here. If you edit -it, just rebuild your VS/Nsight project to sync the changes into the IDE. - -Please run the project without modifications to ensure that everything works -correctly. We have provided a space in `kernel.cu` to run test code, and at the -moment it contains an example of how to use the `Thrust` library to perform -key-value sorting on the GPU. If everything is working right, you should see -some output in the console window and a cube of gray particles. The viewer is -equipped with camera controls: left-click and drag to move the camera view, -right-click and drag vertically to zoom in and out. - -**NOTE: Build the project in `release` mode for performance analysis and capturing.** - -### 1.1. Boids with Naive Neighbor Search - -In the Boids flocking simulation, particles representing birds or fish -(boids) move around the simulation space according to three rules: - -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 - -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 -and compute the velocity change contribution from each of the three rules. -Thus, a bare-bones boids implementation has each boid check every other boid in -the simulation. - -Here is some quick pseudocode to help you out: - -#### Rule 1: Boids try to fly towards the centre of mass of neighbouring boids - -``` -function rule1(Boid boid) - - Vector perceived_center - - foreach Boid b: - if b != boid and distance(b, boid) < rule1Distance then - perceived_center += b.position - endif - end - - perceived_center /= number_of_neighbors - - return (perceived_center - boid.position) * rule1Scale -end -``` - -#### Rule 2: Boids try to keep a small distance away from other objects (including other boids). - -``` -function rule2(Boid boid) - - Vector c = 0 - - foreach Boid b - if b != boid and distance(b, boid) < rule2Distance then - c -= (b.position - boid.position) - endif - end - - return c * rule2Scale -end -``` - -#### Rule 3: Boids try to match velocity with near boids. - -``` -function rule3(Boid boid) - - Vector perceived_velocity - - foreach Boid b - if b != boid and distance(b, boid) < rule3Distance then - perceived_velocity += b.velocity - endif - end - - perceived_velocity /= number_of_neighbors - - 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. - -We also have a simple [2D implementation in Processing](http://studio.sketchpad.cc/sp/pad/view/ro.9cbgCRcgbPOI6/rev.23) -that is very conceptually similar to what you will be writing. Feel free to use -this implementation as a math/code reference. - -For an idea of how the simulation "should" look in 3D, -[here's what our reference implementation looks like.](https://vimeo.com/181547860) - -**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 -``` - -This was for no reason other than that the TA who made the references missed this little detail when this project was originally created, but the results looked about right anyway :facepalm:. - -Following Conrad Parker's pseudocode exactly may lead to unexpected results using the default parameters that we ship with the project, which you may be able to resolve through parameter tuning. If you find good parameters for the "more correct" Boids implementation, feel free to share them on Piazza! - -However, since the purpose of this assignment is to introduce you to CUDA, we recommend at this time that you follow our pseudocode for your initial implementation of the algorithm so you don't have to tune parameters while debugging your implementation. - -### 1.2. Code walkthrough - -* `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. - -### 1.3. Play around - -The parameters we have provided result in a stable simulation using our -reference implementation, but your mileage may vary. Play around with the boid -count as well and see how the simulation responds. - -## Part 2: Let there be (better) flocking! - -### 2.0. A quick explanation of uniform grids - -Recall in part 1 that any two boids can only influence each other if they are -within some *neighborhood distance* of each other. -Based on this observation, we can see that having each boid check every -other boid is very inefficient, especially if (as in our standard parameters) -the number of boids is large and the neighborhood distance is much smaller than -the full simulation space. We can cull a lot of neighbor checks using a -datastructure called a **uniform spatial grid**. - -A uniform grid is made up of cells that are at least as wide as the neighborhood -distance and covers the entire simulation domain. -Before computing the new velocities of the boids, we "bin" them into the grid in -a preprocess step. -![a uniform grid in 2D](images/Boids%20Ugrid%20base.png) - -If the cell width is double the neighborhood distance, each boid only has to be -checked against other boids in 8 cells, or 4 in the 2D case. - -![a uniform grid in 2D with neighborhood and cells to search for some particles shown](images/Boids%20Ugrid%20neighbor%20search%20shown.png) - -You can build a uniform grid on the CPU by iterating over the boids, figuring out -its enclosing cell, and then keeping a pointer to the boid in a resizeable -array representing the cell. However, this doesn't transfer well to the GPU -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. - -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 -boids by these indices, we can ensure that pointers to boids in the same cells -are contiguous in memory. - -Then, we can walk over the array of sorted uniform grid indices and look at -every pair of values. If the values differ, we know that we are at the border -of the representation of two different cells. Storing these locations in a table -with an entry for each cell gives us a complete representation of the uniform -grid. This "table" can just be an array with as much space as there are cells. -This process is data parallel and can be naively parallelized. -![buffers for generating a uniform grid using index sort](images/Boids%20Ugrids%20buffers%20naive.png) - -### 2.1. Code walkthrough - -Instead of having you implement parallel sorting on the GPU for your first -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 -- `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` -- `dev_gridCellEndIndices` - buffer containing a pointer for each cell to the -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 -pointers. - -See the TODOs and LOOKs for Part 2.1 in the code for some pseudocode help. - -You can toggle between different timestep update modes using the defines in -`main.cpp`. - -### 2.2 Play around some more -Compare your uniform grid velocity update to your naive velocity update. -In the typical case, the uniform grid version should be considerably faster. -Try to push the limits of how many boids you can simulate. - -Change the cell width of the uniform grid to be the neighborhood distance, instead of twice the neighborhood distance. Now, 27 neighboring cells will need to be checked for intersection. Does this increase or decrease the efficiency of the flocking? - -### 2.3 Cutting out the middleman -Consider the uniform grid neighbor search outlined in 2.1: pointers to boids in -a single cell are contiguous in memory, but the boid data itself (velocities and -positions) is scattered all over the place. Try rearranging the boid data -itself so that all the velocities and positions of boids in one cell are also -contiguous in memory, so this data can be accessed directly using -`dev_gridCellStartIndices` and `dev_gridCellEndIndices` without -`dev_particleArrayIndices`. - -![buffers for generating a uniform grid using index sort, then making the boid data coherent](images/Boids%20Ugrids%20buffers%20data%20coherent.png) - -See the TODOs for Part 2.3. This should involve a slightly modified copy of -your code from 2.1. - -## Part 3: Performance Analysis -For this project, we will guide you through your performance analysis with some -basic questions. In the future, you will guide your own performance analysis - -but these simple questions will always be critical to answer. In general, we -want you to go above and beyond the suggested performance investigations and -explore how different aspects of your code impact performance as a whole. - -The provided framerate meter (in the window title) will be a useful base -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) -* 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. -* Describe the methodology you are using to benchmark. -* Performance plots are a good thing. - -### 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. - * To do this, change `#define VISUALIZE` to `0`. -* For tighter timing measurement, you can use CUDA events to measure just the - simulation CUDA kernel. Info on this can be found online easily. You will - probably have to average over several simulation steps, similar to the way - FPS is currently calculated. - -This section will not be graded for correctness, but please let us know your -hypotheses and insights. - -**Answer these:** - -* For each implementation, how does changing the number of boids affect -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? -* 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? -* 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! - -**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). -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 - -## Submit - -If you have modified any of the `CMakeLists.txt` files at all (aside from the -list of `SOURCE_FILES`), mention it explicitly. Beware of any build issues discussed on the Google Group. - -Open a GitHub pull request so that we can see that you have finished. -The title should be "Project 1: YOUR NAME". -The template of the comment section of your pull request is attached below, you can do some copy and paste: - -* [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 - * ... -* Feedback on the project itself, if any. - - -And you're done! - -## Tips -- If your simulation crashes before launch, use -`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. -- 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. -- 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) -- 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. If you choose to tackle this problem, have it done before Project 5. 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). diff --git a/README.md b/README.md index d63a6a1..cdc0cfe 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,56 @@ -**University of Pennsylvania, CIS 565: GPU Programming and Architecture, -Project 1 - Flocking** +**University of Pennsylvania, CIS 565: GPU Programming and Architecture** -* (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) +* Jiarui Yan + * [LinkedIn](https://www.linkedin.com/in/jiarui-yan-a06bb5197?lipi=urn%3Ali%3Apage%3Ad_flagship3_profile_view_base_contact_details%3BvRlITiOMSt%2B9Mgg6SZFKDQ%3D%3D), [personal website](https://jiaruiyan.pb.online/), [twitter](https://twitter.com/JerryYan1997), etc. +* Tested on: Windows 10 Home, i7-9700K @ 3.60GHz 16GB DDR4 RAM, RTX 2070 SUPER 8GB Dedicated GPU memory (Personal desktop) -### (TODO: Your README) +# Project 1 - Flocking -Include screenshots, analysis, etc. (Remember, this is public, so don't put -anything here that you don't want to share with the world.) +## Video demos + +![video screen shot](./images/Video_Capture.PNG) + +* [Youtube video link](https://youtu.be/Fz3OPSUzCes) + + +## Performance analysis +For this part, I take first 5000 frames for each test case and record there time in order to keep consistency. Then, I calculate the FPS by dividing frames to time duration. Besides, there are three experiments that I need to take care. In the first experiment, I compare the framerate change with increasing of boids for naive, scattered uniform grid and coherent uniform frid in the condition of visualization. The second experiment is same as the first one except for that it needs to be conducted without visualization. Finally, I conduct an experiment to reveal the effects of different block sizes. + +### Framerate change with increasing number for three methods with visualization + +![Experiment 1](./images/Boids_FPS_with_visualization.PNG) + +### Framerate change with increasing number for three methods without visualization + +![Experiment 2](./images/Boids_FPS_without_visualization.PNG) + +### Framerate change with increasing block size + +![Experiment 3](./images/Different_BlockSize.PNG) + +### Extra Credit -- Grid Optimization +For this part, I use the search radius of a particle to determine whether we need to check a particulour cell. In this way, we can exclude additional cells if a particle is lay exactly at the center line of a dimension. Here is the performance analysis by comparing the Uniform method and Coherent method before and after the optimization. You can find this part of code in velocity update section for these two methods. + +![Experiment 4](./images/Grid_Optimization.PNG) + + +## Questions and answers + +1. For each implementation, how does changing the number of boids affect performance? Why do you think this is? + + With the increasing of the number of boids, the performance would be decreasing. Generally speaking, I think it is caused by computational requirement increased by the number of boids. Specifically, there are several reasons. First, with the increase of number of boids, we need to update more particles and let GPU deal with more threads, which require more time. Second, with more particles, we have more neighbors for each particle. This means it signiciantly increases the number of accessing memory, which would also require more time. + +2. For each implementation, how does changing the block count and block size affect performance? Why do you think this is? + + For each implementation, the block size doesn't affect performance with an apparent rule. I think this is caused by the limit of cache size or file system size in each SM, because different implementation uses differnet amount of memory and there may be different sweet pots for different implementations. + +3. 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? + + Yes, but it depends. It gives less performance improvements by compare to uniform scattered method when the number of boids is small. However, it can give more improvement when the number of boids is large. I think this is a result that I expect, because coherent uniform grid method has more consumtption when it shuffle position and velocity. However, it is an o(1) consumption in parallel, which is less than the consumption of an additional level of indirect memory access when the number of boids is large. + +4. 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! + + Yes, it would affect the performance. According to my experiment conducted by using uniform scatterred grid with 50000 boids, checking 27 neighboring cells gives 941.553 FPS and checking 9 neighboring cells gives 650.395 FPS. Therefore, I think, with a finner grid, we can exclude lots of unnecessary particle checking, which can give us more performance than checking more particles with less checking cell number. + +## Acknowledgement +* [CIS565 2020 Fall Course link](https://cis565-fall-2020.github.io/) diff --git a/excel_data/AVG FPS for different boids.xlsx b/excel_data/AVG FPS for different boids.xlsx new file mode 100644 index 0000000..bab6768 Binary files /dev/null and b/excel_data/AVG FPS for different boids.xlsx differ diff --git a/images/.gitignore b/images/.gitignore deleted file mode 100644 index e69de29..0000000 diff --git a/images/Boids Ugrid base.png b/images/Boids Ugrid base.png deleted file mode 100644 index 7fc82dc..0000000 Binary files a/images/Boids Ugrid base.png and /dev/null differ diff --git a/images/Boids Ugrid neighbor search shown.png b/images/Boids Ugrid neighbor search shown.png deleted file mode 100644 index 1913caf..0000000 Binary files a/images/Boids Ugrid neighbor search shown.png and /dev/null differ diff --git a/images/Boids Ugrids buffers data coherent.png b/images/Boids Ugrids buffers data coherent.png deleted file mode 100644 index 5fc4f3a..0000000 Binary files a/images/Boids Ugrids buffers data coherent.png and /dev/null differ diff --git a/images/Boids Ugrids buffers naive.png b/images/Boids Ugrids buffers naive.png deleted file mode 100644 index ebb2c36..0000000 Binary files a/images/Boids Ugrids buffers naive.png and /dev/null differ diff --git a/images/Boids_FPS_with_visualization.PNG b/images/Boids_FPS_with_visualization.PNG new file mode 100644 index 0000000..5aa4cf7 Binary files /dev/null and b/images/Boids_FPS_with_visualization.PNG differ diff --git a/images/Boids_FPS_without_visualization.PNG b/images/Boids_FPS_without_visualization.PNG new file mode 100644 index 0000000..202814b Binary files /dev/null and b/images/Boids_FPS_without_visualization.PNG differ diff --git a/images/Different_BlockSize.PNG b/images/Different_BlockSize.PNG new file mode 100644 index 0000000..53a564b Binary files /dev/null and b/images/Different_BlockSize.PNG differ diff --git a/images/Grid_Optimization.PNG b/images/Grid_Optimization.PNG new file mode 100644 index 0000000..bb426c1 Binary files /dev/null and b/images/Grid_Optimization.PNG differ diff --git a/images/UnlockFPS.png b/images/UnlockFPS.png deleted file mode 100644 index 0b538d4..0000000 Binary files a/images/UnlockFPS.png and /dev/null differ diff --git a/images/Video_Capture.PNG b/images/Video_Capture.PNG new file mode 100644 index 0000000..dd28929 Binary files /dev/null and b/images/Video_Capture.PNG differ diff --git a/src/kernel.cu b/src/kernel.cu index 74dffcb..7930ac7 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -37,6 +37,7 @@ void checkCUDAError(const char *msg, int line = -1) { *****************/ /*! Block size used for CUDA kernel launch. */ +// #define blockSize 128 #define blockSize 128 // LOOK-1.2 Parameters for the boids algorithm. @@ -83,8 +84,12 @@ thrust::device_ptr dev_thrust_particleGridIndices; 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. +glm::vec3* dev_reshuffle_pos; +glm::vec3* dev_reshuffle_vel1; +glm::vec3* dev_reshuffle_vel2; // LOOK-2.1 - Grid parameters based on simulation parameters. // These are automatically computed for you in Boids::initSimulation @@ -158,6 +163,8 @@ void Boids::initSimulation(int N) { // LOOK-2.1 computing grid params gridCellWidth = 2.0f * std::max(std::max(rule1Distance, rule2Distance), rule3Distance); + // 2.2: + // gridCellWidth = std::max(std::max(rule1Distance, rule2Distance), rule3Distance); int halfSideCount = (int)(scene_scale / gridCellWidth) + 1; gridSideCount = 2 * halfSideCount; @@ -169,6 +176,34 @@ void Boids::initSimulation(int N) { gridMinimum.z -= halfGridWidth; // TODO-2.1 TODO-2.3 - Allocate additional buffers here. + // Buffer containing a pointer for each boid to its data in dev_pos and dev_vel1 and dev_vel2 + // Table 2 second row: + cudaMalloc((void**)&dev_particleArrayIndices, N * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_particleArrayIndices failed!"); + + // Buffer containing the grid index of each boid + // Table 2 first row: + cudaMalloc((void**)&dev_particleGridIndices, N * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_particleGridIndices failed!"); + + // Table 1: + // Buffer containing a pointer for each cell to the beginning of its data in dev_particleArrayIndices + cudaMalloc((void**)&dev_gridCellStartIndices, gridCellCount * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_gridCellStartIndices failed!"); + + cudaMalloc((void**)&dev_gridCellEndIndices, gridCellCount * sizeof(int)); + checkCUDAErrorWithLine("cudaMalloc dev_gridCellEndIndices failed!"); + + // TODO-2.3 additional buffers: + cudaMalloc((void**)&dev_reshuffle_pos, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_reshuffle_pos failed!"); + + cudaMalloc((void**)&dev_reshuffle_vel1, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_reshuffle_vel1 failed!"); + + cudaMalloc((void**)&dev_reshuffle_vel2, N * sizeof(glm::vec3)); + checkCUDAErrorWithLine("cudaMalloc dev_reshuffle_vel2 failed!"); + cudaDeviceSynchronize(); } @@ -230,21 +265,72 @@ void Boids::copyBoidsToVBO(float *vbodptr_positions, float *vbodptr_velocities) * in the `pos` and `vel` arrays. */ __device__ glm::vec3 computeVelocityChange(int N, int iSelf, const glm::vec3 *pos, const glm::vec3 *vel) { - // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves - // Rule 2: boids try to stay a distance d away from each other - // Rule 3: boids try to match the speed of surrounding boids - return glm::vec3(0.0f, 0.0f, 0.0f); + glm::vec3 temp_self_pos = pos[iSelf]; + glm::vec3 perceived_center(0.0, 0.0, 0.0); + glm::vec3 rule2_vel(0.0, 0.0, 0.0); + glm::vec3 rule3_vel(0.0, 0.0, 0.0); + glm::vec3 res_vel_change(0.0, 0.0, 0.0); + int rule1_neighbor_count = 0; + int rule3_neighbor_count = 0; + + for (int iTar = 0; iTar < N; ++iTar) + { + if (iTar == iSelf) continue; + glm::vec3 temp_tar_pos = pos[iTar]; + glm::vec3 temp_tar_vel = vel[iTar]; + float distance = (float) glm::distance(temp_self_pos, temp_tar_pos); + // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves. + if (distance < rule1Distance) { + perceived_center += temp_tar_pos; + rule1_neighbor_count ++; + } + // Rule 2: boids try to stay a distance d away from each other. + if (distance < rule2Distance) { + rule2_vel -= (temp_tar_pos - temp_self_pos); + } + // Rule 3: boids try to match the speed of surrounding boids. + if (distance < rule3Distance) { + rule3_vel += temp_tar_vel; + rule3_neighbor_count ++; + } + } + // Rule 1 result: + if (rule1_neighbor_count > 0) { + perceived_center /= rule1_neighbor_count; + glm::vec3 rule1_vel = (perceived_center - temp_self_pos) * rule1Scale; + res_vel_change += rule1_vel; + } + // Rule 2 result: + res_vel_change += (rule2_vel * rule2Scale); + // Rule 3 result: + if (rule3_neighbor_count > 0) { + rule3_vel /= rule3_neighbor_count; + res_vel_change += (rule3_vel * rule3Scale); + } + + return res_vel_change; } + + /** * 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) { + int index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index >= N) { + return; + } // Compute a new velocity based on pos and vel1 + glm::vec3 vel_change = computeVelocityChange(N, index, pos, vel1); + vel2[index] += vel_change; // Clamp the speed - // Record the new velocity into vel2. Question: why NOT vel1? + // Record the new velocity into vel2. Question: why NOT vel1? + if (glm::length(vel2[index]) > maxSpeed) { + vel2[index] = maxSpeed * glm::normalize(vel2[index]); + } } /** @@ -289,6 +375,16 @@ __global__ void kernComputeIndices(int N, int gridResolution, // - 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 particle_index = threadIdx.x + (blockIdx.x * blockDim.x); + if (particle_index >= N) { + return; + } + indices[particle_index] = particle_index; + glm::vec3 particle_relative_pos = pos[particle_index] - gridMin; + int grid_x_idx = glm::floor(particle_relative_pos[0] * inverseCellWidth); + int grid_y_idx = glm::floor(particle_relative_pos[1] * inverseCellWidth); + int grid_z_idx = glm::floor(particle_relative_pos[2] * inverseCellWidth); + gridIndices[particle_index] = gridIndex3Dto1D(grid_x_idx, grid_y_idx, grid_z_idx, gridResolution); } // LOOK-2.1 Consider how this could be useful for indicating that a cell @@ -302,10 +398,38 @@ __global__ void kernResetIntBuffer(int N, int *intBuffer, int value) { __global__ void kernIdentifyCellStartEnd(int N, int *particleGridIndices, int *gridCellStartIndices, int *gridCellEndIndices) { + // NOTE: N is the number of particles. + int table_index = threadIdx.x + (blockIdx.x * blockDim.x); + if (table_index >= N) { + return; + } // 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 curr_grid_idx = particleGridIndices[table_index]; + // First table element: + if (table_index == 0) { + gridCellStartIndices[curr_grid_idx] = 0; + } + // Last table element: + if (table_index == N - 1) { + gridCellEndIndices[curr_grid_idx] = N - 1; + } + // Look at left: + if (table_index != 0) { + int left_grid_index = particleGridIndices[table_index - 1]; + if (left_grid_index != curr_grid_idx) { + gridCellStartIndices[curr_grid_idx] = table_index; + } + } + // Look at right: + if (table_index != N - 1) { + int right_grid_index = particleGridIndices[table_index + 1]; + if (right_grid_index != curr_grid_idx) { + gridCellEndIndices[curr_grid_idx] = table_index; + } + } } __global__ void kernUpdateVelNeighborSearchScattered( @@ -314,14 +438,198 @@ __global__ void kernUpdateVelNeighborSearchScattered( int *gridCellStartIndices, int *gridCellEndIndices, int *particleArrayIndices, glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) { + int particle_index = threadIdx.x + (blockIdx.x * blockDim.x); + if (particle_index >= N) { + return; + } // TODO-2.1 - Update a boid's velocity using the uniform grid to reduce // the number of boids that need to be checked. // - Identify the grid cell that this particle is in - // - Identify which cells may contain neighbors. This isn't always 8. - // - For each cell, read the start/end indices in the boid pointer array. - // - Access each boid in the cell and compute velocity change from - // the boids rules, if this boid is within the neighborhood distance. + glm::vec3 particle_relative_pos = pos[particle_index] - gridMin; + + glm::vec3 temp_self_pos = pos[particle_index]; + glm::vec3 perceived_center(0.0, 0.0, 0.0); + glm::vec3 rule2_vel(0.0, 0.0, 0.0); + glm::vec3 rule3_vel(0.0, 0.0, 0.0); + glm::vec3 res_vel_change(0.0, 0.0, 0.0); + int rule1_neighbor_count = 0; + int rule3_neighbor_count = 0; + + // - Identify which cells may contain neighbors. This isn't always 8. --- Grid width : distance = 2 : 1. + glm::vec3 round_cell_relative_pos = glm::round(particle_relative_pos * inverseCellWidth); + for (int i = -1; i < 1; ++i) { + for (int j = -1; j < 1; ++j) { + for (int k = -1; k < 1; ++k) { + int neighbor_x_idx = i + (int)round_cell_relative_pos[0]; + int neighbor_y_idx = j + (int)round_cell_relative_pos[1]; + int neighbor_z_idx = k + (int)round_cell_relative_pos[2]; + // This neighbor cell is out of boundary: + if (neighbor_x_idx < 0 || neighbor_y_idx < 0 || neighbor_z_idx < 0 || neighbor_x_idx >= gridResolution || neighbor_y_idx >= gridResolution || neighbor_z_idx >= gridResolution) { + continue; + } + // Get the existing neighbor cell: + int curr_neighbor_idx = gridIndex3Dto1D(neighbor_x_idx, neighbor_y_idx, neighbor_z_idx, gridResolution); + // - For each cell, read the start/end indices in the boid pointer array. + int table_start_idx = gridCellStartIndices[curr_neighbor_idx]; + int table_end_idx = gridCellEndIndices[curr_neighbor_idx]; + if (table_start_idx < 0 || table_end_idx < 0) { + // This cell doesn't enclose a particle + continue; + } + // - Access each boid in the cell and compute velocity change from + // the boids rules, if this boid is within the neighborhood distance. + for (int iTar = table_start_idx; iTar <= table_end_idx; ++iTar) + { + int neighbor_particle_idx = particleArrayIndices[iTar]; + if (neighbor_particle_idx == particle_index) continue; + glm::vec3 temp_tar_pos = pos[neighbor_particle_idx]; + glm::vec3 temp_tar_vel = vel1[neighbor_particle_idx]; + float distance = (float)glm::distance(temp_self_pos, temp_tar_pos); + // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves. + if (distance < rule1Distance) { + perceived_center += temp_tar_pos; + rule1_neighbor_count++; + } + // Rule 2: boids try to stay a distance d away from each other. + if (distance < rule2Distance) { + rule2_vel -= (temp_tar_pos - temp_self_pos); + } + // Rule 3: boids try to match the speed of surrounding boids. + if (distance < rule3Distance) { + rule3_vel += temp_tar_vel; + rule3_neighbor_count++; + } + } + } + } + } + + + // Extra Credit -- Grid Optimization: + /* + float maxDistance = glm::max(glm::max(rule1Distance, rule2Distance), rule3Distance); + glm::vec3 radius_vector(maxDistance, maxDistance, maxDistance); + glm::vec3 upper_particle_search = temp_self_pos + radius_vector; + glm::vec3 lower_particle_search = temp_self_pos - radius_vector; + glm::ivec3 upper_cell_idx = glm::floor((upper_particle_search - gridMin) * inverseCellWidth); + glm::ivec3 lower_cell_idx = glm::floor((lower_particle_search - gridMin) * inverseCellWidth); + for (int neighbor_z_idx = lower_cell_idx[2]; neighbor_z_idx <= upper_cell_idx[2]; ++neighbor_z_idx) { + for (int neighbor_y_idx = lower_cell_idx[1]; neighbor_y_idx <= upper_cell_idx[1]; ++neighbor_y_idx) { + for (int neighbor_x_idx = lower_cell_idx[0]; neighbor_x_idx <= upper_cell_idx[0]; ++neighbor_x_idx) { + // This neighbor cell is out of boundary: + if (neighbor_x_idx < 0 || neighbor_y_idx < 0 || neighbor_z_idx < 0 || neighbor_x_idx >= gridResolution || neighbor_y_idx >= gridResolution || neighbor_z_idx >= gridResolution) { + continue; + } + // Get the existing neighbor cell: + int curr_neighbor_idx = gridIndex3Dto1D(neighbor_x_idx, neighbor_y_idx, neighbor_z_idx, gridResolution); + // - For each cell, read the start/end indices in the boid pointer array. + int table_start_idx = gridCellStartIndices[curr_neighbor_idx]; + int table_end_idx = gridCellEndIndices[curr_neighbor_idx]; + if (table_start_idx < 0 || table_end_idx < 0) { + // This cell doesn't enclose a particle + continue; + } + // - Access each boid in the cell and compute velocity change from + // the boids rules, if this boid is within the neighborhood distance. + for (int iTar = table_start_idx; iTar <= table_end_idx; ++iTar) + { + int neighbor_particle_idx = particleArrayIndices[iTar]; + if (neighbor_particle_idx == particle_index) continue; + glm::vec3 temp_tar_pos = pos[neighbor_particle_idx]; + glm::vec3 temp_tar_vel = vel1[neighbor_particle_idx]; + float distance = (float)glm::distance(temp_self_pos, temp_tar_pos); + // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves. + if (distance < rule1Distance) { + perceived_center += temp_tar_pos; + rule1_neighbor_count++; + } + // Rule 2: boids try to stay a distance d away from each other. + if (distance < rule2Distance) { + rule2_vel -= (temp_tar_pos - temp_self_pos); + } + // Rule 3: boids try to match the speed of surrounding boids. + if (distance < rule3Distance) { + rule3_vel += temp_tar_vel; + rule3_neighbor_count++; + } + } + } + } + } + */ + + /* + // TODO-2.2: Identify which cells may contain neighbors. This isn't always 27. --- Grid width : distance = 1 : 1. + int grid_x_idx = glm::floor(particle_relative_pos[0] * inverseCellWidth); + int grid_y_idx = glm::floor(particle_relative_pos[1] * inverseCellWidth); + int grid_z_idx = glm::floor(particle_relative_pos[2] * inverseCellWidth); + glm::vec3 self_cell_idx(grid_x_idx, grid_y_idx, grid_z_idx); + for (int i = -1; i < 2; ++i) { + for (int j = -1; j < 2; ++j) { + for (int k = -1; k < 2; ++k) { + int neighbor_x_idx = i + (int)self_cell_idx[0]; + int neighbor_y_idx = j + (int)self_cell_idx[1]; + int neighbor_z_idx = k + (int)self_cell_idx[2]; + // This neighbor cell is out of boundary: + if (neighbor_x_idx < 0 || neighbor_y_idx < 0 || neighbor_z_idx < 0 || neighbor_x_idx >= gridResolution || neighbor_y_idx >= gridResolution || neighbor_z_idx >= gridResolution) { + continue; + } + // Get the existing neighbor cell: + int curr_neighbor_idx = gridIndex3Dto1D(neighbor_x_idx, neighbor_y_idx, neighbor_z_idx, gridResolution); + // - For each cell, read the start/end indices in the boid pointer array. + int table_start_idx = gridCellStartIndices[curr_neighbor_idx]; + int table_end_idx = gridCellEndIndices[curr_neighbor_idx]; + if (table_start_idx < 0 || table_end_idx < 0) { + // This cell doesn't enclose a particle + continue; + } + // - Access each boid in the cell and compute velocity change from + // the boids rules, if this boid is within the neighborhood distance. + for (int iTar = table_start_idx; iTar <= table_end_idx; ++iTar) + { + int neighbor_particle_idx = particleArrayIndices[iTar]; + if (neighbor_particle_idx == particle_index) continue; + glm::vec3 temp_tar_pos = pos[neighbor_particle_idx]; + glm::vec3 temp_tar_vel = vel1[neighbor_particle_idx]; + float distance = (float)glm::distance(temp_self_pos, temp_tar_pos); + // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves. + if (distance < rule1Distance) { + perceived_center += temp_tar_pos; + rule1_neighbor_count++; + } + // Rule 2: boids try to stay a distance d away from each other. + if (distance < rule2Distance) { + rule2_vel -= (temp_tar_pos - temp_self_pos); + } + // Rule 3: boids try to match the speed of surrounding boids. + if (distance < rule3Distance) { + rule3_vel += temp_tar_vel; + rule3_neighbor_count++; + } + } + } + } + } + */ + + // Rule 1 result: + if (rule1_neighbor_count > 0) { + perceived_center /= rule1_neighbor_count; + glm::vec3 rule1_vel = (perceived_center - temp_self_pos) * rule1Scale; + res_vel_change += rule1_vel; + } + // Rule 2 result: + res_vel_change += (rule2_vel * rule2Scale); + // Rule 3 result: + if (rule3_neighbor_count > 0) { + rule3_vel /= rule3_neighbor_count; + res_vel_change += (rule3_vel * rule3Scale); + } // - Clamp the speed change before putting the new speed in vel2 + vel2[particle_index] += res_vel_change; + if (glm::length(vel2[particle_index]) > maxSpeed) { + vel2[particle_index] = maxSpeed * glm::normalize(vel2[particle_index]); + } } __global__ void kernUpdateVelNeighborSearchCoherent( @@ -329,59 +637,349 @@ __global__ void kernUpdateVelNeighborSearchCoherent( float inverseCellWidth, float cellWidth, int *gridCellStartIndices, int *gridCellEndIndices, glm::vec3 *pos, glm::vec3 *vel1, glm::vec3 *vel2) { + int particle_index = threadIdx.x + (blockIdx.x * blockDim.x); + if (particle_index >= N) { + return; + } // 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 + glm::vec3 particle_relative_pos = pos[particle_index] - gridMin; + + glm::vec3 temp_self_pos = pos[particle_index]; + glm::vec3 perceived_center(0.0, 0.0, 0.0); + glm::vec3 rule2_vel(0.0, 0.0, 0.0); + glm::vec3 rule3_vel(0.0, 0.0, 0.0); + glm::vec3 res_vel_change(0.0, 0.0, 0.0); + int rule1_neighbor_count = 0; + int rule3_neighbor_count = 0; + + // - Identify which cells may contain neighbors. This isn't always 8. --- Grid width : distance = 2 : 1. + + glm::vec3 round_cell_relative_pos = glm::round(particle_relative_pos * inverseCellWidth); + for (int i = -1; i < 1; ++i) { + for (int j = -1; j < 1; ++j) { + for (int k = -1; k < 1; ++k) { + int neighbor_x_idx = i + (int)round_cell_relative_pos[0]; + int neighbor_y_idx = j + (int)round_cell_relative_pos[1]; + int neighbor_z_idx = k + (int)round_cell_relative_pos[2]; + // This neighbor cell is out of boundary: + if (neighbor_x_idx < 0 || neighbor_y_idx < 0 || neighbor_z_idx < 0 || neighbor_x_idx >= gridResolution || neighbor_y_idx >= gridResolution || neighbor_z_idx >= gridResolution) { + continue; + } + // Get the existing neighbor cell: + int curr_neighbor_idx = gridIndex3Dto1D(neighbor_x_idx, neighbor_y_idx, neighbor_z_idx, gridResolution); + // - 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. + int table_start_idx = gridCellStartIndices[curr_neighbor_idx]; + int table_end_idx = gridCellEndIndices[curr_neighbor_idx]; + if (table_start_idx < 0 || table_end_idx < 0) { + // This cell doesn't enclose a particle + continue; + } + // - Access each boid in the cell and compute velocity change from + // the boids rules, if this boid is within the neighborhood distance. + for (int iTar = table_start_idx; iTar <= table_end_idx; ++iTar) + { + int neighbor_particle_idx = iTar; + if (neighbor_particle_idx == particle_index) continue; + glm::vec3 temp_tar_pos = pos[neighbor_particle_idx]; + glm::vec3 temp_tar_vel = vel1[neighbor_particle_idx]; + float distance = (float)glm::distance(temp_self_pos, temp_tar_pos); + // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves. + if (distance < rule1Distance) { + perceived_center += temp_tar_pos; + rule1_neighbor_count++; + } + // Rule 2: boids try to stay a distance d away from each other. + if (distance < rule2Distance) { + rule2_vel -= (temp_tar_pos - temp_self_pos); + } + // Rule 3: boids try to match the speed of surrounding boids. + if (distance < rule3Distance) { + rule3_vel += temp_tar_vel; + rule3_neighbor_count++; + } + } + } + } + } + + /* + // Extra Credit -- Grid Optimization: + float maxDistance = glm::max(glm::max(rule1Distance, rule2Distance), rule3Distance); + glm::vec3 radius_vector(maxDistance, maxDistance, maxDistance); + glm::vec3 upper_particle_search = temp_self_pos + radius_vector; + glm::vec3 lower_particle_search = temp_self_pos - radius_vector; + glm::ivec3 upper_cell_idx = glm::floor((upper_particle_search - gridMin) * inverseCellWidth); + glm::ivec3 lower_cell_idx = glm::floor((lower_particle_search - gridMin) * inverseCellWidth); + for (int neighbor_x_idx = lower_cell_idx[0]; neighbor_x_idx <= upper_cell_idx[0]; ++neighbor_x_idx) { + for (int neighbor_y_idx = lower_cell_idx[1]; neighbor_y_idx <= upper_cell_idx[1]; ++neighbor_y_idx) { + for (int neighbor_z_idx = lower_cell_idx[2]; neighbor_z_idx <= upper_cell_idx[2]; ++neighbor_z_idx) { + // This neighbor cell is out of boundary: + if (neighbor_x_idx < 0 || neighbor_y_idx < 0 || neighbor_z_idx < 0 || neighbor_x_idx >= gridResolution || neighbor_y_idx >= gridResolution || neighbor_z_idx >= gridResolution) { + continue; + } + // Get the existing neighbor cell: + int curr_neighbor_idx = gridIndex3Dto1D(neighbor_x_idx, neighbor_y_idx, neighbor_z_idx, gridResolution); + // - For each cell, read the start/end indices in the boid pointer array. + int table_start_idx = gridCellStartIndices[curr_neighbor_idx]; + int table_end_idx = gridCellEndIndices[curr_neighbor_idx]; + if (table_start_idx < 0 || table_end_idx < 0) { + // This cell doesn't enclose a particle + continue; + } + // - Access each boid in the cell and compute velocity change from + // the boids rules, if this boid is within the neighborhood distance. + for (int iTar = table_start_idx; iTar <= table_end_idx; ++iTar) + { + int neighbor_particle_idx = iTar; + if (neighbor_particle_idx == particle_index) continue; + glm::vec3 temp_tar_pos = pos[neighbor_particle_idx]; + glm::vec3 temp_tar_vel = vel1[neighbor_particle_idx]; + float distance = (float)glm::distance(temp_self_pos, temp_tar_pos); + // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves. + if (distance < rule1Distance) { + perceived_center += temp_tar_pos; + rule1_neighbor_count++; + } + // Rule 2: boids try to stay a distance d away from each other. + if (distance < rule2Distance) { + rule2_vel -= (temp_tar_pos - temp_self_pos); + } + // Rule 3: boids try to match the speed of surrounding boids. + if (distance < rule3Distance) { + rule3_vel += temp_tar_vel; + rule3_neighbor_count++; + } + } + } + } + } + */ + + // TODO-2.2: Identify which cells may contain neighbors. This isn't always 27. --- Grid width : distance = 1 : 1. + /* + int grid_x_idx = glm::floor(particle_relative_pos[0] * inverseCellWidth); + int grid_y_idx = glm::floor(particle_relative_pos[1] * inverseCellWidth); + int grid_z_idx = glm::floor(particle_relative_pos[2] * inverseCellWidth); + glm::vec3 self_cell_idx(grid_x_idx, grid_y_idx, grid_z_idx); + for (int i = -1; i < 2; ++i) { + for (int j = -1; j < 2; ++j) { + for (int k = -1; k < 2; ++k) { + int neighbor_x_idx = i + (int)self_cell_idx[0]; + int neighbor_y_idx = j + (int)self_cell_idx[1]; + int neighbor_z_idx = k + (int)self_cell_idx[2]; + // This neighbor cell is out of boundary: + if (neighbor_x_idx < 0 || neighbor_y_idx < 0 || neighbor_z_idx < 0 || neighbor_x_idx >= gridResolution || neighbor_y_idx >= gridResolution || neighbor_z_idx >= gridResolution) { + continue; + } + // Get the existing neighbor cell: + int curr_neighbor_idx = gridIndex3Dto1D(neighbor_x_idx, neighbor_y_idx, neighbor_z_idx, gridResolution); + // - For each cell, read the start/end indices in the boid pointer array. + int table_start_idx = gridCellStartIndices[curr_neighbor_idx]; + int table_end_idx = gridCellEndIndices[curr_neighbor_idx]; + if (table_start_idx < 0 || table_end_idx < 0) { + // This cell doesn't enclose a particle + continue; + } + // - Access each boid in the cell and compute velocity change from + // the boids rules, if this boid is within the neighborhood distance. + for (int iTar = table_start_idx; iTar <= table_end_idx; ++iTar) + { + int neighbor_particle_idx = iTar; + if (neighbor_particle_idx == particle_index) continue; + glm::vec3 temp_tar_pos = pos[neighbor_particle_idx]; + glm::vec3 temp_tar_vel = vel1[neighbor_particle_idx]; + float distance = (float)glm::distance(temp_self_pos, temp_tar_pos); + // Rule 1: boids fly towards their local perceived center of mass, which excludes themselves. + if (distance < rule1Distance) { + perceived_center += temp_tar_pos; + rule1_neighbor_count++; + } + // Rule 2: boids try to stay a distance d away from each other. + if (distance < rule2Distance) { + rule2_vel -= (temp_tar_pos - temp_self_pos); + } + // Rule 3: boids try to match the speed of surrounding boids. + if (distance < rule3Distance) { + rule3_vel += temp_tar_vel; + rule3_neighbor_count++; + } + } + } + } + } + */ + + // Rule 1 result: + if (rule1_neighbor_count > 0) { + perceived_center /= rule1_neighbor_count; + glm::vec3 rule1_vel = (perceived_center - temp_self_pos) * rule1Scale; + res_vel_change += rule1_vel; + } + // Rule 2 result: + res_vel_change += (rule2_vel * rule2Scale); + // Rule 3 result: + if (rule3_neighbor_count > 0) { + rule3_vel /= rule3_neighbor_count; + res_vel_change += (rule3_vel * rule3Scale); + } + // - Clamp the speed change before putting the new speed in vel2 + vel2[particle_index] = vel1[particle_index] + res_vel_change; + if (glm::length(vel2[particle_index]) > maxSpeed) { + vel2[particle_index] = maxSpeed * glm::normalize(vel2[particle_index]); + } } /** * Step the entire N-body simulation by `dt` seconds. */ void Boids::stepSimulationNaive(float dt) { - // TODO-1.2 - use the kernels you wrote to step the simulation forward in time. - // TODO-1.2 ping-pong the velocity buffers + dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); + // TODO-1.2 - use the kernels you wrote to step the simulation forward in time. + // Update velocity: + kernUpdateVelocityBruteForce <<< fullBlocksPerGrid, blockSize >>> (numObjects, dev_pos, dev_vel1, dev_vel2); + checkCUDAErrorWithLine("kernUpdateVelocityBruteForce failed!"); + // Update position: + kernUpdatePos <<< fullBlocksPerGrid, blockSize >>> (numObjects, dt, dev_pos, dev_vel2); + checkCUDAErrorWithLine("kernUpdateVelocityBruteForce failed!"); + // TODO-1.2 ping-pong the velocity buffers + dev_vel1 = dev_vel2; +} + +__global__ void particle_grid_idx_velocity_test(int N, int* particleGridIndices, glm::vec3* pos, glm::vec3* vel2) +{ + int index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index >= N) { + return; + } + vel2[index][0] = (float)particleGridIndices[index]; + vel2[index][1] = (float)particleGridIndices[index]; + vel2[index][2] = (float)particleGridIndices[index]; + if (glm::length(vel2[index]) > maxSpeed) { + vel2[index] = maxSpeed * glm::normalize(vel2[index]); + } } +__global__ void sorted_idx_velocity_test(int N, int* particleGridIndices, int* gridCellStartIndices, int* gridCellEndIndices, + int* particleArrayIndices, glm::vec3* pos, glm::vec3* vel2) +{ + int table_index = threadIdx.x + (blockIdx.x * blockDim.x); + if (table_index >= N) { + return; + } + int curr_particle_index = particleArrayIndices[table_index]; + vel2[curr_particle_index][0] = (float)particleGridIndices[table_index]; + vel2[curr_particle_index][1] = (float)particleGridIndices[table_index]; + vel2[curr_particle_index][2] = (float)particleGridIndices[table_index]; + if (glm::length(vel2[curr_particle_index]) > maxSpeed) { + vel2[curr_particle_index] = maxSpeed * glm::normalize(vel2[curr_particle_index]); + } +} + +// NOTE: dev_particleGridIndices[p_idx] gives the grid idx of the p_idx particle. void Boids::stepSimulationScatteredGrid(float dt) { + + dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); + dim3 cell_block_num((gridCellCount + blockSize - 1) / blockSize); + kernResetIntBuffer <<< cell_block_num, blockSize >>> (gridCellCount, dev_gridCellStartIndices, -1); + checkCUDAErrorWithLine("kernResetIntBuffer failed!"); + kernResetIntBuffer <<< cell_block_num, blockSize >>> (gridCellCount, dev_gridCellEndIndices, -1); + checkCUDAErrorWithLine("kernResetIntBuffer failed!"); // TODO-2.1 // Uniform Grid Neighbor search using Thrust sort. // In Parallel: // - label each particle with its array index as well as its grid index. // Use 2x width grids. + kernComputeIndices <<< fullBlocksPerGrid, blockSize >>> (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, dev_pos, dev_particleArrayIndices, dev_particleGridIndices); + checkCUDAErrorWithLine("kernComputeIndices failed!"); // - Unstable key sort using Thrust. A stable sort isn't necessary, but you // are welcome to do a performance comparison. + thrust::device_ptr dev_thrust_keys(dev_particleGridIndices); + thrust::device_ptr dev_thrust_values(dev_particleArrayIndices); + thrust::sort_by_key(dev_thrust_keys, dev_thrust_keys + numObjects, dev_thrust_values); // - Naively unroll the loop for finding the start and end indices of each // cell's data pointers in the array of boid indices + kernIdentifyCellStartEnd <<< fullBlocksPerGrid, blockSize >>> (numObjects, dev_particleGridIndices, dev_gridCellStartIndices, dev_gridCellEndIndices); + checkCUDAErrorWithLine("kernIdentifyCellStartEnd failed!"); + // - Perform velocity updates using neighbor search + kernUpdateVelNeighborSearchScattered <<< fullBlocksPerGrid, blockSize >>> (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, gridCellWidth, + dev_gridCellStartIndices, dev_gridCellEndIndices, dev_particleArrayIndices, dev_pos, dev_vel1, dev_vel2); + checkCUDAErrorWithLine("kernUpdateVelNeighborSearchScattered failed!"); + // - Update positions + kernUpdatePos <<< fullBlocksPerGrid, blockSize >>> (numObjects, dt, dev_pos, dev_vel2); + checkCUDAErrorWithLine("kernUpdatePos failed!"); // - Ping-pong buffers as needed + dev_vel1 = dev_vel2; + +} + +__global__ void reshuffle_particle_data(int N, int* particleArrayIndices, glm::vec3* pos, glm::vec3* vel1, glm::vec3* pos_reshuffle, glm::vec3* vel1_reshuffle) +{ + int table_index = threadIdx.x + (blockIdx.x * blockDim.x); + if (table_index >= N) { + return; + } + int particle_index = particleArrayIndices[table_index]; + pos_reshuffle[table_index] = pos[particle_index]; + vel1_reshuffle[table_index] = vel1[particle_index]; +} + + +__global__ void shuffle_back_particle_data(int N, int* particleArrayIndices, glm::vec3* vel2, glm::vec3* vel2_reshuffle) +{ + int table_index = threadIdx.x + (blockIdx.x * blockDim.x); + if (table_index >= N) { + return; + } + int particle_index = particleArrayIndices[table_index]; + vel2[particle_index] = vel2_reshuffle[table_index]; } void Boids::stepSimulationCoherentGrid(float dt) { + dim3 fullBlocksPerGrid((numObjects + blockSize - 1) / blockSize); + dim3 cell_block_num((gridCellCount + blockSize - 1) / blockSize); + kernResetIntBuffer << < cell_block_num, blockSize >> > (gridCellCount, dev_gridCellStartIndices, -1); + checkCUDAErrorWithLine("kernResetIntBuffer failed!"); + kernResetIntBuffer << < cell_block_num, blockSize >> > (gridCellCount, dev_gridCellEndIndices, -1); + checkCUDAErrorWithLine("kernResetIntBuffer failed!"); // 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 + kernComputeIndices << < fullBlocksPerGrid, blockSize >> > (numObjects, gridSideCount, gridMinimum, gridInverseCellWidth, dev_pos, dev_particleArrayIndices, dev_particleGridIndices); + checkCUDAErrorWithLine("kernComputeIndices failed!"); // - Unstable key sort using Thrust. A stable sort isn't necessary, but you // are welcome to do a performance comparison. + thrust::device_ptr dev_thrust_keys(dev_particleGridIndices); + thrust::device_ptr dev_thrust_values(dev_particleArrayIndices); + thrust::sort_by_key(dev_thrust_keys, dev_thrust_keys + numObjects, dev_thrust_values); // - Naively unroll the loop for finding the start and end indices of each // cell's data pointers in the array of boid indices + kernIdentifyCellStartEnd << < fullBlocksPerGrid, blockSize >> > (numObjects, dev_particleGridIndices, dev_gridCellStartIndices, dev_gridCellEndIndices); + checkCUDAErrorWithLine("kernIdentifyCellStartEnd failed!"); // - BIG DIFFERENCE: use the rearranged array index buffer to reshuffle all // the particle data in the simulation array. // CONSIDER WHAT ADDITIONAL BUFFERS YOU NEED + reshuffle_particle_data <<< fullBlocksPerGrid, blockSize >>> (numObjects, dev_particleArrayIndices, dev_pos, dev_vel1, dev_reshuffle_pos, dev_reshuffle_vel1); + checkCUDAErrorWithLine("reshuffle_particle_data failed!"); // - Perform velocity updates using neighbor search + kernUpdateVelNeighborSearchCoherent <<< fullBlocksPerGrid, blockSize >>> (numObjects, gridSideCount, gridMinimum, + gridInverseCellWidth, gridCellWidth, dev_gridCellStartIndices, dev_gridCellEndIndices, dev_reshuffle_pos, dev_reshuffle_vel1, dev_reshuffle_vel2); + checkCUDAErrorWithLine("kernUpdateVelNeighborSearchCoherent failed!"); + shuffle_back_particle_data <<< fullBlocksPerGrid, blockSize >>> (numObjects, dev_particleArrayIndices, dev_vel2, dev_reshuffle_vel2); // - Update positions + kernUpdatePos <<< fullBlocksPerGrid, blockSize >>> (numObjects, dt, dev_pos, dev_vel2); + checkCUDAErrorWithLine("kernUpdatePos failed!"); // - Ping-pong buffers as needed. THIS MAY BE DIFFERENT FROM BEFORE. + dev_vel1 = dev_vel2; } void Boids::endSimulation() { @@ -390,6 +988,13 @@ 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_reshuffle_pos); + cudaFree(dev_reshuffle_vel1); + cudaFree(dev_reshuffle_vel2); } void Boids::unitTest() { diff --git a/src/main.cpp b/src/main.cpp index b82c8c6..ad6a788 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -14,11 +14,14 @@ // LOOK-2.1 LOOK-2.3 - toggles for UNIFORM_GRID and COHERENT_GRID #define VISUALIZE 1 -#define UNIFORM_GRID 0 +#define UNIFORM_GRID 1 #define COHERENT_GRID 0 // LOOK-1.2 - change this to adjust particle count in the simulation -const int N_FOR_VIS = 5000; +// Uniform grid: const int N_FOR_VIS = 200000; ~= 70FPS +// Navie: const int N_FOR_VIS = 50000; ~= 40FPS +// TODO-2.2: It looks like the increased grid width can improve the efficiency. +const int N_FOR_VIS = 50000; const float DT = 0.2f; /** @@ -217,6 +220,10 @@ void initShaders(GLuint * program) { double timebase = 0; int frame = 0; + double total_fps = 0; + double start_time = glfwGetTime(); + int total_frame = 0; + Boids::unitTest(); // LOOK-1.2 We run some basic example code to make sure // your CUDA development setup is ready to go. @@ -224,6 +231,9 @@ void initShaders(GLuint * program) { glfwPollEvents(); frame++; + + total_frame++; + double time = glfwGetTime(); if (time - timebase > 1.0) { @@ -234,6 +244,12 @@ void initShaders(GLuint * program) { runCUDA(); + if (total_frame == 5000) { + double end_time = glfwGetTime(); + double avg_fps = total_frame / (end_time - start_time); + std::cout << "avg fps:" << avg_fps << std::endl; + } + std::ostringstream ss; ss << "["; ss.precision(1);