Chapter 1: Introduction
4.2 The Volume Method Step
4.2.2 Transferring Particle Velocities to Grid
We explain the scenario of velocity transfer in the two dimensional case first.
In Figure 22, the black grid is the staggered MAC grid. The black circle represents the position of a particle. The orange arrow is the particle‟s velocity vector u⃗ = (u, v) (u⃗ = (u, v, w) in 3D). We are going to transfer the u-component of the velocity onto the center of vertical grid faces and the v-component onto the center of horizontal grid faces. We transfer velocities to the center of faces within the distance of one grid cell size. We refer to the center of faces where the u-component is stored as the grid‟s u-faces and the center of faces where the v-component is stored as the grid‟s v-faces in the rest of the paragraphs. We draw blue squares on the four u-faces that have a nonzero transferred value and purple (magenta) squares on the four v-faces with nonzero transferred value.
Figure 21: Two example configurations of the markers
Figure 22: A 2D example of transferring the velocity of one particle
The idea of transferring velocity u-component to the grid‟s u-faces is illustrated in Figure 23. The blue square resembles the u-faces. The arrows originated from the squares indicate the horizontal velocities on the u-faces.
In the following paragraphs, we describe how to traverse the particles in each grid cells parallelly. And we use this procedure to transfer the velocities on the particles to the grid. We follow the work of [Gre07] and modify their algorithm.
Algorithm overview: The procedure of visiting all the grid cells, VISITALLGRIDCELL, is listed in Figure 25. After the arrays are initialized, we build the data used for the traversal in the BUILDTRAVERSALDATA procedure (line 13). Then we visit each of the cell parallely from line 16 to line 19.
The BUILDTRAVERSALDATA procedure is listed in Figure 26. We sort the data (position, velocity) of the particles by its grid cell id (line 17). The cell id of a particle is the id of the cell where this particle is located. Then the data of the particles in the same cell are placed together in the data arrays. We record the array index of the first particle of each cell in an array start_index from line 20 to line 25. If there is no
Figure 23: An example of transferring u Figure 24: An example of transferring v
particle in a cell then the value of the array is the initial value Np, the number of particles.
The procedure of visiting a grid cell, VISITAGRIDCELL, is listed in Figure 27. This procedure takes the sorted data of the particles, position, velocity and cell id, and the lookup data start_index as input. We first get the array index of first particle by looking up in the start_index. Then we iterate through the particles in the for-loop (line 11). If the cell does not contain any particle, then the value of i_start is Np (the default value). In this case the termination condition of the for-loop is reached and the procedure ends. In the for-loop, there is another termination condition when the next cell id is different from the first id (line 15). It is because the data of the particles in a cell are placed together in the array. If the cell id is different from the first particle, then it means this particle belongs to anther cell. Consequently, the for-loop is exited.
Figure 25: The VISITALLGRIDCELL procedure
1 procedure VISITALLGRIDCELL(
16 for each (i,j,k) 1<=i<nx-1, 1<=j<ny-1, 1<=k<nz-1 in parallel 17 compute cell id : id=i+nx*(j+ ny*k)
18 VISITAGRIDCELL ( id, Np, sorted_cell_id, 19 p_sorted, v_ sorted, start_index ) 20 end
Figure 26: The BUILDTRAVERSALDATA procedure
1 procedure BUILDTRAVERSALDATA ( 2 in Np // number of particles 3 in Ncell // number of cells 4 in p[Np] // particle‟s position 5 in v[Np] // particle‟s velocity
6 out sorted_cell_id [Np] // particle‟s cell id after sorting 7 out p_sorted [Np] //
8 out v_sorted [Np] // particle‟s velocity after sorting 9 out start_index [Ncell]
10 // starting particle‟s index in p_sorted, v_sorted array ) 11 begin
12 // initialization
13 for each i, 0<=i<Ncell in parallel 14 start_index[i] = Np
15 //sorting
16 compute particle‟s cell id by its position in parallel 17 sort particle‟s position and velocity by its cell id 18 output: p_sorted[], v_sorted[], sorted_cell_id[]
19
20 for each i, 0<=i<Np in parallel 21 if i!= 0 then
22 if sorted_cell_id[i-1]!= sorted_cell_id[i] then 23 start_index[sorted_cell_id[i]]=i
24 else
25 start_index[sorted_cell_id[i]]=0 26 end
Figure 27: The VISITAGRIDCELL procedure
We give a two dimensional example in Figure 28 and the configuration of the data structures is shown in Figure 29. The configuration of the particles is in shown in Figure 28. Each square represents a cell and the black number in each cell is the id of this cell. The blue circles represent the particles. The particles actually are volume-less in our model, but we draw the circles big enough for visualization. The white numbers in the circles are the id of the particles. For the computing the id of a cell, we use the following equation:
cell(i, j) id = i + 𝑛𝑥 ∙ j (3)
nx is the number of grid cells in the x dimension, which is 4 in this case.
We start running the VISITALLGRIDCELL procedure, after initialization (line 8-11), we step into the BUILDTRAVERSALDATA procedure.
In line 13-14 of BUILDTRAVERSALDATA, each element of the start_index array is initialized the number of total particles, Np, which is 6.
In line 16, the particle ids are computed. We show the idea in Figure 28. We find where the particles are and see the cell ids, which are the numbers labeled on the cells.
For example, we see particle with label 0 is located at the cell label 9, so its cell id is 9.
We sort the id of cells and id of the particles according to the id of cells. The cell ids of the particles before sorting and after sorting are shown in the 2nd and 3rd column in Figure 29. The particle data (position and velocity) before sorting and after sorting are shown in the 4th and 5th column in Figure 29.
In line 20-25 of BUILDTRAVERSALDATA, we fill data into the start_index array.
When i=0, sorted_cell_id[0] is 4, then start_index[4] is 0.
When i=1, sorted_cell_id[0] == sorted_cell_id[1], nothing changed.
When i=2, sorted_cell_id[1] is 4 and sorted_cell_id[2] is 6, then start_index[6] is 2.
When i=3, sorted_cell_id[2] == sorted_cell_id[3], nothing changed.
When i=4, sorted_cell_id[3] == sorted_cell_id[4], nothing changed.
When i=5, sorted_cell_id[4] is 6 and sorted_cell_id[5] is 9, then start_index[9] is 5.
The result of star_index array is shown in the last column in Figure 29.
Figure 28: Configuration of the particles
Figure 29: An example of running the BUILDTRAVERSALDATA procedure
We step out of the BUILDTRAVERSALDATA procedure and go to line 16 of VISITALLGRIDCELL. We visit all the cells in the parallel for-loop and step into VISITAGRIDCELL in line 18.
We choose the cell with id 6 as the example to show the VISITAGRIDCELL
procedure. In line 10 of VISITAGRIDCELL, the i_start = start_index[ 6 ] is 2. Then the for-loop in line 11 becomes: for each i, 2 (i_start) <=i<6 (Np). The id_start is 6 When i=2, we visit particle with P1, V1. id_start is 6 and sorted_cell_id[3] is 6, the for-loop continues.
When i=3, we visit particle with P2, V2. id_start is 6 and sorted_cell_id[4] is 6, the for-loop continues.
When i=4, we visit particle with P4, V4. id_start is 6 and sorted_cell_id[5] is 9, the termination condition is reached and we exit the for-loop. Then the procedure ends.
The above operation is shown in Figure 30. In the figure, the first arrow indicates that we get i_start by looking-up start_index[6]. The particle data that are visited is colored with blue. The second arrow indicates that we reach the termination coniditon of the for-loop when id_start (6) and sorted_cell_id[ i+1 ] (9) is different.
We run the VISITAGRIDCELL procedure with the cell with id 9 as another example.
In line 10 of VISITAGRIDCELL, the i_start = start_index[ 9 ] is 5. Then the for-loop in line 11 becomes: for each i, 5(i_start) <=i<6 (Np). The id_start is 9
When i=5, we visit particle with P0, V0. i+1 is 6 >= Np, the for-loop terminates.
Let us choose the cell with id 7 (without any particle) and run the VISITAGRIDCELL
procedure. In line 10 of VISITAGRIDCELL, the i_start = start_index[ 7 ] is 6. Then the for-loop in line 11 becomes: for each i, 6(i_start) <=i<6 (Np). After that the procedure ends.
We apply the VISITALLGRIDCELL procedure to transfer one of the velocity components. In the body of the parallel for loop (line 17-19), we visit all the neighboring cells instead of visit just one cell. The VISITAGRIDCELL procedure is modified to output the weighting and the weighted velocity component of the visited particle. And then the velocity outputs from the VISITAGRIDCELL are accumulated and stored to the grid.
The three components (u, v, w) of the 𝑢⃗ = (𝑢, 𝑣, 𝑤) are stored on the staggered structure. When we transfer u, we shift the origin of the grid by half a cell size along the y, z direction. So that the location of where u becomes the grid points. Similarly, we shift the origin of the grid along the x, z direction by half a cell size when transferring v. We run the modified VISITALLGRIDCELL procedure three times to transfer the velocity components u, v and w, respectively.
We store the data used for traversal on the texture memory for acceleration because of the uncoalesced memory access pattern.
Figure 30: An example of running the VISITAGRIDCELL procedure