sorting.cuh - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
sorting.cuh (6233B)
---
1 #ifndef SORTING_CUH_
2 #define SORTING_CUH_
3
4 // Returns the cellID containing the particle, based cubic grid
5 // See Bayraktar et al. 2009
6 // Kernel is executed on the device, and is callable from the device only
7 __device__ unsigned int calcCellID(Float3 x)
8 {
9 unsigned int i_x, i_y, i_z;
10
11 // Calculate integral coordinates:
12 i_x = floor((x.x - devC_grid.origo[0]) / (devC_grid.L[0]/devC_grid.num[0]));
13 i_y = floor((x.y - devC_grid.origo[1]) / (devC_grid.L[1]/devC_grid.num[1]));
14 i_z = floor((x.z - devC_grid.origo[2]) / (devC_grid.L[2]/devC_grid.num[2]));
15
16 // Integral coordinates are converted to 1D coordinate:
17 return (i_z * devC_grid.num[1])
18 * devC_grid.num[0] + i_y * devC_grid.num[0] + i_x;
19
20 } // End of calcCellID(...)
21
22
23 // Calculate hash value for each particle, based on position in grid.
24 // Kernel executed on device, and callable from host only.
25 __global__ void calcParticleCellID(unsigned int* dev_gridParticleCellID,
26 unsigned int* dev_gridParticleIndex,
27 Float4* dev_x)
28 {
29 unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;
30
31 if (idx < devC_np) { // Condition prevents block size error
32
33 //volatile Float4 x = dev_x[idx]; // Ensure coalesced read
34 Float4 x = dev_x[idx]; // Ensure coalesced read
35
36 unsigned int cellID = calcCellID(MAKE_FLOAT3(x.x, x.y, x.z));
37
38 // Check for NaN
39 if (x.x != x.x || x.y != x.y || x.z != x.z)
40 printf("\ncalcParticleCellID: Error! NaN encountered. "
41 "idx = %d: cellID = %d, "
42 "x = %f,%f,%f\n",
43 idx, cellID, x.x, x.y, x.z);
44
45 // Store values
46 __syncthreads();
47 dev_gridParticleCellID[idx] = cellID;
48 dev_gridParticleIndex[idx] = idx;
49
50 }
51 } // End of calcParticleCellID(...)
52
53
54 // Reorder particle data into sorted order, and find the start and end particle
55 // indexes of each cell in the sorted hash array.
56 #ifdef SPHERE_GPU
57 __global__ void reorderArrays(unsigned int* dev_cellStart,
58 unsigned int* dev_cellEnd,
59 unsigned int* dev_gridParticleCellID,
60 unsigned int* dev_gridParticleIndex,
61 Float4* dev_x,
62 Float4* dev_vel,
63 Float4* dev_angvel,
64 Float4* dev_x_sorted,
65 Float4* dev_vel_sorted,
66 Float4* dev_angvel_sorted)
67 {
68
69 // Create hash array in shared on-chip memory. The size of the array
70 // (threadsPerBlock + 1) is determined at launch time (extern notation).
71 extern __shared__ unsigned int shared_data[];
72
73 // Thread index in block
74 unsigned int tidx = threadIdx.x;
75
76 // Thread index in grid
77 unsigned int idx = threadIdx.x + blockIdx.x * blockDim.x;
78
79 // CellID hash value of particle idx
80 unsigned int cellID;
81
82 // Read cellID data and store it in shared memory (shared_data)
83 if (idx < devC_np) { // Condition prevents block size error
84 cellID = dev_gridParticleCellID[idx];
85
86 // Load hash data into shared memory, allowing access to neighbor
87 // particle cellID values
88 shared_data[tidx+1] = cellID;
89
90 if (idx > 0 && tidx == 0) {
91 // First thread in block must load neighbor particle hash
92 shared_data[0] = dev_gridParticleCellID[idx-1];
93 }
94 }
95 //if (cellID != 0)
96 //printf("reorderArrays: %d,%d\tcellID = %d\n", tidx, idx, cellID);
97
98 // Pause completed threads in this block, until all
99 // threads are done loading data into shared memory
100 __syncthreads();
101
102 // Find lowest and highest particle index in each cell
103 if (idx < devC_np) { // Condition prevents block size error
104 // If this particle has a different cell index to the previous particle,
105 // it's the first particle in the cell -> Store the index of this
106 // particle in the cell. The previous particle must be the last particle
107 // in the previous cell.
108 if (idx == 0 || cellID != shared_data[tidx]) {
109 dev_cellStart[cellID] = idx;
110 if (idx > 0)
111 dev_cellEnd[shared_data[tidx]] = idx;
112 }
113
114 // Check wether the thread is the last one
115 if (idx == (devC_np - 1))
116 dev_cellEnd[cellID] = idx + 1;
117
118
119 // Use the sorted index to reorder the position and velocity data
120 unsigned int sortedIndex = dev_gridParticleIndex[idx];
121
122 // Fetch from global read
123 Float4 x = dev_x[sortedIndex];
124 Float4 vel = dev_vel[sortedIndex];
125 Float4 angvel = dev_angvel[sortedIndex];
126
127 __syncthreads();
128 // Write sorted data to global memory
129 dev_x_sorted[idx] = x;
130 dev_vel_sorted[idx] = vel;
131 dev_angvel_sorted[idx] = angvel;
132 }
133 } // End of reorderArrays(...)
134
135 #else
136 // CPU variant: identical result, but reads the neighbor hash from global
137 // memory instead of the shared-memory staging used on the GPU.
138 __global__ void reorderArrays(unsigned int* dev_cellStart,
139 unsigned int* dev_cellEnd,
140 unsigned int* dev_gridParticleCellID,
141 unsigned int* dev_gridParticleIndex,
142 Float4* dev_x,
143 Float4* dev_vel,
144 Float4* dev_angvel,
145 Float4* dev_x_sorted,
146 Float4* dev_vel_sorted,
147 Float4* dev_angvel_sorted)
148 {
149 unsigned int idx = threadIdx.x + blockIdx.x * blockDim.x;
150
151 if (idx < devC_np) { // Condition prevents block size error
152 const unsigned int cellID = dev_gridParticleCellID[idx];
153
154 // Find lowest and highest particle index in each cell
155 if (idx == 0 || cellID != dev_gridParticleCellID[idx-1]) {
156 dev_cellStart[cellID] = idx;
157 if (idx > 0)
158 dev_cellEnd[dev_gridParticleCellID[idx-1]] = idx;
159 }
160 if (idx == devC_np - 1)
161 dev_cellEnd[cellID] = idx + 1;
162
163 // Use the sorted index to reorder the position and velocity data
164 const unsigned int sortedIndex = dev_gridParticleIndex[idx];
165 dev_x_sorted[idx] = dev_x[sortedIndex];
166 dev_vel_sorted[idx] = dev_vel[sortedIndex];
167 dev_angvel_sorted[idx] = dev_angvel[sortedIndex];
168 }
169 } // End of reorderArrays(...)
170 #endif // SPHERE_GPU
171
172 #endif
173 // vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4