HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
BCSlipWall_GPU.cu
Go to the documentation of this file.
1 /*! @file BCSlipWall_GPU.cu
2  @author Youngdae Kim
3  @brief GPU implmentation of slip-wall boundary conditions
4 */
5 #include <stdlib.h>
6 #include <basic_gpu.h>
7 #include <arrayfunctions_gpu.h>
8 #include <boundaryconditions.h>
9 
10 #include <physicalmodels/navierstokes2d.h>
11 #include <physicalmodels/navierstokes3d.h>
12 
13 #ifdef CUDA_VAR_ORDERDING_AOS
14 
15 /*! Kernel for 2D implementation of gpuBCSlipWallU() */
16 __global__
17 void BCSlipWallU_dim2_kernel(
18  int grid_size,
19  int dim,
20  int face,
21  int ghosts,
22  int ndims,
23  int nvars,
24  double gamma,
25  const int *size,
26  const int *ie,
27  const int *is,
28  const double *FlowVelocity,
29  double *phi
30 )
31 {
32  int index = threadIdx.x + (blockDim.x * blockIdx.x);
33  if (index < grid_size) {
34  const int max_ndims = 3;
35  int p1, p2;
36  int bounds[max_ndims], indexb[max_ndims], indexi[max_ndims];
37  double inv_gamma_m1 = 1.0/(gamma-1.0);
38 
39  _ArraySubtract1D_(bounds,ie,is,ndims);
40  _ArrayIndexnD_(ndims,index,bounds,indexb,0);
41  _ArrayCopy1D_(indexb,indexi,ndims);
42  _ArrayAdd1D_(indexi,indexi,is,ndims);
43  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
44  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
45  else return;
46  _ArrayIndex1DWO_(ndims,size,indexb,is,ghosts,p1);
47  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
48 
49  /* flow variables in the interior */
50  double rho, uvel, vvel, energy, pressure;
51  double rho_gpt, uvel_gpt, vvel_gpt, energy_gpt, pressure_gpt;
52  _NavierStokes2DGetFlowVar_((phi+nvars*p2),rho,uvel,vvel,energy,pressure,gamma);
53  /* set the ghost point values */
54  rho_gpt = rho;
55  pressure_gpt = pressure;
56  if (dim == _XDIR_) {
57  uvel_gpt = 2.0*FlowVelocity[_XDIR_] - uvel;
58  vvel_gpt = vvel;
59  } else if (dim == _YDIR_) {
60  uvel_gpt = uvel;
61  vvel_gpt = 2.0*FlowVelocity[_YDIR_] - vvel;
62  } else {
63  uvel_gpt = 0.0;
64  vvel_gpt = 0.0;
65  }
66  energy_gpt = inv_gamma_m1*pressure_gpt
67  + 0.5 * rho_gpt * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt);
68 
69  phi[nvars*p1+0] = rho_gpt;
70  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
71  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
72  phi[nvars*p1+3] = energy_gpt;
73  }
74  return;
75 }
76 
77 /*! Kernel for 3D implementation of gpuBCSlipWallU() */
78 __global__
79 void BCSlipWallU_dim3_kernel(
80  int grid_size,
81  int dim,
82  int face,
83  int ghosts,
84  int ndims,
85  int nvars,
86  double gamma,
87  const int *size,
88  const int *ie,
89  const int *is,
90  const double *FlowVelocity,
91  double *phi
92 )
93 {
94  int index = threadIdx.x + (blockDim.x * blockIdx.x);
95  if (index < grid_size) {
96  const int max_ndims = 3;
97  int p1, p2;
98  int bounds[max_ndims], indexb[max_ndims], indexi[max_ndims];
99  double inv_gamma_m1 = 1.0/(gamma-1.0);
100 
101  _ArraySubtract1D_(bounds,ie,is,ndims);
102  _ArrayIndexnD_(ndims,index,bounds,indexb,0);
103  _ArrayCopy1D_(indexb,indexi,ndims);
104  _ArrayAdd1D_(indexi,indexi,is,ndims);
105  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
106  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
107  else return;
108  _ArrayIndex1DWO_(ndims,size,indexb,is,ghosts,p1);
109  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
110 
111  /* flow variables in the interior */
112  double rho, uvel, vvel, wvel, energy, pressure;
113  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
114  _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
115  /* set the ghost point values */
116  rho_gpt = rho;
117  pressure_gpt = pressure;
118  if (dim == _XDIR_) {
119  uvel_gpt = 2.0*FlowVelocity[_XDIR_] - uvel;
120  vvel_gpt = vvel;
121  wvel_gpt = wvel;
122  } else if (dim == _YDIR_) {
123  uvel_gpt = uvel;
124  vvel_gpt = 2.0*FlowVelocity[_YDIR_] - vvel;
125  wvel_gpt = wvel;
126  } else if (dim == _ZDIR_) {
127  uvel_gpt = uvel;
128  vvel_gpt = vvel;
129  wvel_gpt = 2.0*FlowVelocity[_ZDIR_] - wvel;
130  } else {
131  uvel_gpt = 0.0;
132  vvel_gpt = 0.0;
133  wvel_gpt = 0.0;
134  }
135  energy_gpt = inv_gamma_m1*pressure_gpt
136  + 0.5 * rho_gpt
137  * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);
138 
139  phi[nvars*p1+0] = rho_gpt;
140  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
141  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
142  phi[nvars*p1+3] = rho_gpt * wvel_gpt;
143  phi[nvars*p1+4] = energy_gpt;
144  }
145  return;
146 }
147 
148 /*! Applies the slip-wall boundary condition: This is specific to the two and three
149  dimensional Euler and Navier-Stokes systems (#Euler2D, #NavierStokes2D, #NavierStokes3D).
150  It is used for simulating inviscid walls or symmetric boundaries. The pressure, density,
151  and tangential velocity at the ghost points are extrapolated from the interior, while the
152  normal velocity at the ghost points is set such that the interpolated value at the boundary
153  face is equal to the specified wall velocity.
154 
155  \sa BCSlipWallU()
156 */
157 extern "C" int gpuBCSlipWallU(
158  void *b, /*!< Boundary object of type #DomainBoundary */
159  void *m, /*!< MPI object of type #MPIVariables */
160  int ndims, /*!< Number of spatial dimensions */
161  int nvars, /*!< Number of variables/DoFs per grid point */
162  int *gpu_size, /*!< Integer array with the number of grid points in each spatial dimension */
163  int ghosts, /*!< Number of ghost points */
164  double *gpu_phi, /*!< The solution array on which to apply the boundary condition */
165  double waqt /*!< Current solution time */
166  )
167 {
168  DomainBoundary *boundary = (DomainBoundary*) b;
169 
170  int dim = boundary->dim;
171  int face = boundary->face;
172 
173  if (ndims == 2) {
174 
175  if (boundary->on_this_proc) {
176  int bounds[ndims];
177  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
178  int ngrid_points = 1; for(int i = 0; i < ndims; i++) ngrid_points *= bounds[i];
179  int nblocks = (ngrid_points - 1) / GPU_THREADS_PER_BLOCK + 1;
180 
181  gpuMemcpy( boundary->gpu_ie, boundary->ie, ndims*sizeof(int), gpuMemcpyHostToDevice );
182  gpuMemcpy( boundary->gpu_is, boundary->is, ndims*sizeof(int), gpuMemcpyHostToDevice );
183 
184  BCSlipWallU_dim2_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
185  ngrid_points, dim, face, ghosts, ndims, nvars, boundary->gamma,
186  gpu_size, boundary->gpu_ie, boundary->gpu_is, boundary->gpu_FlowVelocity, gpu_phi);
187  }
188 
189  } else if (ndims == 3) {
190 
191  if (boundary->on_this_proc) {
192  int bounds[ndims];
193  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
194  int ngrid_points = 1; for(int i = 0; i < ndims; i++) ngrid_points *= bounds[i];
195  int nblocks = (ngrid_points - 1) / GPU_THREADS_PER_BLOCK + 1;
196 
197  BCSlipWallU_dim3_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
198  ngrid_points, dim, face, ghosts, ndims, nvars, boundary->gamma,
199  gpu_size, boundary->gpu_ie, boundary->gpu_is, boundary->gpu_FlowVelocity, gpu_phi);
200  }
201 
202  } else {
203 
204  fprintf(stderr,"Error in gpuBCSlipWallU(): not implemented for ndims=%d.\n", ndims);
205  return 1;
206 
207  }
208 
209  return 0;
210 }
211 
212 #else
213 
214 /*! Kernel for 3D implementation of gpuBCSlipWallU() */
215 __global__
216 void BCSlipWallU_dim3_kernel(
217  int npoints_bounds,
218  int npoints_local_wghosts,
219  int face,
220  int ndims,
221  int dim,
222  int ghosts,
223  int nvars,
224  double gamma,
225  const int * __restrict__ bounds,
226  const int * __restrict__ size,
227  const int * __restrict__ boundary_is,
228  const double * __restrict__ FlowVelocity,
229  double * __restrict__ phi
230 )
231 {
232  int index = threadIdx.x + (blockDim.x * blockIdx.x);
233  if (index < npoints_bounds) {
234  const int max_ndims = 3;
235  int p1, p2;
236  int indexb[max_ndims], indexi[max_ndims];
237  double inv_gamma_m1 = 1.0/(gamma-1.0);
238 
239  _ArrayIndexnD_(ndims,index,bounds,indexb,0);
240  _ArrayCopy1D_(indexb,indexi,ndims);
241  _ArrayAdd1D_(indexi,indexi,boundary_is,ndims);
242  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
243  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
244  else return;
245  _ArrayIndex1DWO_(ndims,size,indexb,boundary_is,ghosts,p1);
246  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
247 
248  /* flow variables in the interior */
249  double rho, uvel, vvel, wvel, energy, pressure;
250  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
251  _NavierStokes3DGetFlowVar_((phi+p2),npoints_local_wghosts,rho,uvel,vvel,wvel,energy,pressure,gamma);
252  /* set the ghost point values */
253  rho_gpt = rho;
254  pressure_gpt = pressure;
255  if (dim == _XDIR_) {
256  uvel_gpt = 2.0*FlowVelocity[_XDIR_] - uvel;
257  vvel_gpt = vvel;
258  wvel_gpt = wvel;
259  } else if (dim == _YDIR_) {
260  uvel_gpt = uvel;
261  vvel_gpt = 2.0*FlowVelocity[_YDIR_] - vvel;
262  wvel_gpt = wvel;
263  } else if (dim == _ZDIR_) {
264  uvel_gpt = uvel;
265  vvel_gpt = vvel;
266  wvel_gpt = 2.0*FlowVelocity[_ZDIR_] - wvel;
267  } else {
268  uvel_gpt = 0.0;
269  vvel_gpt = 0.0;
270  wvel_gpt = 0.0;
271  }
272  energy_gpt = inv_gamma_m1*pressure_gpt
273  + 0.5 * rho_gpt
274  * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);
275 
276  phi[p1+0] = rho_gpt;
277  phi[p1+1*npoints_local_wghosts] = rho_gpt * uvel_gpt;
278  phi[p1+2*npoints_local_wghosts] = rho_gpt * vvel_gpt;
279  phi[p1+3*npoints_local_wghosts] = rho_gpt * wvel_gpt;
280  phi[p1+4*npoints_local_wghosts] = energy_gpt;
281  }
282  return;
283 }
284 
285 /*! Applies the slip-wall boundary condition: This is specific to the two and three
286  dimensional Euler and Navier-Stokes systems (#Euler2D, #NavierStokes2D, #NavierStokes3D).
287  It is used for simulating inviscid walls or symmetric boundaries. The pressure, density,
288  and tangential velocity at the ghost points are extrapolated from the interior, while the
289  normal velocity at the ghost points is set such that the interpolated value at the boundary
290  face is equal to the specified wall velocity.
291 
292  \sa BCSlipWallU()
293 */
294 extern "C" int gpuBCSlipWallU(
295  void * __restrict__ b, /*!< Boundary object of type #DomainBoundary */
296  void * __restrict__ m, /*!< MPI object of type #MPIVariables */
297  int ndims, /*!< Number of spatial dimensions */
298  int nvars, /*!< Number of variables/DoFs per grid point */
299  int * __restrict__ size, /*!< Integer array with the number of grid points in each spatial dimension */
300  int ghosts, /*!< Number of ghost points */
301  double * __restrict__ phi, /*!< The solution array on which to apply the boundary condition */
302  double waqt /*!< Current solution time */
303  )
304 {
305  DomainBoundary *boundary = (DomainBoundary*) b;
306 
307  int dim = boundary->dim;
308  int face = boundary->face;
309 
310  if (ndims == 3) {
311 
312  if (boundary->on_this_proc) {
313  int nblocks = (boundary->gpu_npoints_bounds - 1) / GPU_THREADS_PER_BLOCK + 1;
314 
315  BCSlipWallU_dim3_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
316  boundary->gpu_npoints_bounds, boundary->gpu_npoints_local_wghosts,
317  face, ndims, dim, ghosts, nvars, boundary->gamma,
318  boundary->gpu_bounds, size, boundary->gpu_is, boundary->gpu_FlowVelocity, phi);
319  }
320 
321  } else {
322 
323  fprintf(stderr,"Error in gpuBCSlipWallU(): not implemented for ndims=%d.\n", ndims);
324  return 1;
325 
326  }
327 
328 
329  return 0;
330 }
331 
332 #endif