HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
NavierStokes2DSource_GPU.cu
Go to the documentation of this file.
1 /*! @file NavierStokes2DSource_GPU.cu
2  @author Youngdae Kim
3  @brief Compute the gravitational source term for the 2D Navier Stokes system
4 */
5 #include <stdlib.h>
6 #include <basic_gpu.h>
7 #include <arrayfunctions.h>
8 #include <physicalmodels/navierstokes2d.h>
9 #include <mpivars.h>
10 #include <hypar.h>
11 
12 #ifdef CUDA_VAR_ORDERDING_AOS
13 
14 static int gpuNavierStokes2DSourceFunction (double*,const double*,const double*,void*,void*,double,int);
15 static int gpuNavierStokes2DSourceUpwind (double*,const double*,const double*,const double*,int,void*,double);
16 
17 __global__
18 void NavierStokes2DSource_Xdir_kernel(
19  int ngrid_points,
20  int ghosts,
21  int ndims,
22  double gamma,
23  double RT,
24  const int *dim,
25  const double *dxinv,
26  const double *grav_field_f,
27  const double *SourceI,
28  const double *u,
29  double *source
30 )
31 {
32  int tx = threadIdx.x + (blockDim.x * blockIdx.x);
33  if (tx < ngrid_points) {
34  int v,p,p1,p2;
35  int index[GPU_MAX_NDIMS],index1[GPU_MAX_NDIMS],index2[GPU_MAX_NDIMS],dim_interface[GPU_MAX_NDIMS];
36 
37  _ArrayIndexnD_(ndims,tx,dim,index,0);
38  _ArrayCopy1D_(dim,dim_interface,ndims); dim_interface[_XDIR_]++;
39  _ArrayCopy1D_(index,index1,ndims);
40  _ArrayCopy1D_(index,index2,ndims); index2[_XDIR_]++;
41  _ArrayIndex1D_(ndims,dim ,index ,ghosts,p );
42  _ArrayIndex1D_(ndims,dim_interface,index1,0 ,p1);
43  _ArrayIndex1D_(ndims,dim_interface,index2,0 ,p2);
44 
45  double dx_inverse; _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,dxinv,dx_inverse);
46  double rho, uvel, vvel, e, P; _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*p),rho,uvel,vvel,e,P,gamma);
47  double term[_MODEL_NVARS_] = {0.0, rho*RT, 0.0, rho*RT*vvel};
48  for (v=0; v<_MODEL_NVARS_; v++) {
49  source[_MODEL_NVARS_*p+v] += ( (term[v]*grav_field_f[p])
50  * (SourceI[_MODEL_NVARS_*p2+v]-SourceI[_MODEL_NVARS_*p1+v])*dx_inverse );
51  }
52  uvel = P; /* useless statement to avoid compiler warnings */
53  }
54  return;
55 }
56 
57 __global__
58 void NavierStokes2DSource_Ydir_kernel(
59  int ngrid_points,
60  int ghosts,
61  int ndims,
62  double gamma,
63  double RT,
64  const int *dim,
65  const double *dxinv,
66  const double *grav_field_f,
67  const double *SourceI,
68  const double *u,
69  double *source
70 )
71 {
72  int tx = threadIdx.x + (blockDim.x * blockIdx.x);
73  if (tx < ngrid_points) {
74  int v,p,p1,p2;
75  int index[GPU_MAX_NDIMS],index1[GPU_MAX_NDIMS],index2[GPU_MAX_NDIMS],dim_interface[GPU_MAX_NDIMS];
76 
77  _ArrayIndexnD_(ndims,tx,dim,index,0);
78  _ArrayCopy1D_(dim,dim_interface,ndims); dim_interface[_YDIR_]++;
79  _ArrayCopy1D_(index,index1,ndims);
80  _ArrayCopy1D_(index,index2,ndims); index2[_YDIR_]++;
81  _ArrayIndex1D_(ndims,dim ,index ,ghosts,p );
82  _ArrayIndex1D_(ndims,dim_interface,index1,0 ,p1);
83  _ArrayIndex1D_(ndims,dim_interface,index2,0 ,p2);
84 
85  double dy_inverse; _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,dxinv,dy_inverse);
86  double rho, uvel, vvel, e, P; _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*p),rho,uvel,vvel,e,P,gamma);
87  double term[_MODEL_NVARS_] = {0.0, 0.0, rho*RT, rho*RT*vvel};
88  for (v=0; v<_MODEL_NVARS_; v++) {
89  source[_MODEL_NVARS_*p+v] += ( (term[v]*grav_field_f[p])
90  * (SourceI[_MODEL_NVARS_*p2+v]-SourceI[_MODEL_NVARS_*p1+v])*dy_inverse );
91  }
92  uvel = P; /* useless statement to avoid compiler warnings */
93  }
94  return;
95 }
96 
97 __global__
98 void NavierStokes2DSourceUpwind_kernel(
99  int ngrid_points,
100  const double *fL,
101  const double *fR,
102  double *fI
103 )
104 {
105  int p = threadIdx.x + (blockDim.x * blockIdx.x);
106  if (p < ngrid_points) {
107  for (int k = 0; k < _MODEL_NVARS_; k++)
108  (fI+_MODEL_NVARS_*p)[k] = 0.5 * ((fL+_MODEL_NVARS_*p)[k] + (fR+_MODEL_NVARS_*p)[k]);
109  }
110  return;
111 }
112 
113 __global__
114 void NavierStokes2DSourceFunction_kernel(
115  int ngrid_points,
116  int dir,
117  const double *grav_field_g,
118  double *f
119 )
120 {
121  int p = threadIdx.x + (blockDim.x * blockIdx.x);
122  if (p < ngrid_points) {
123  (f+_MODEL_NVARS_*p)[0] = 0.0;
124  (f+_MODEL_NVARS_*p)[1] = grav_field_g[p] * (dir == _XDIR_);
125  (f+_MODEL_NVARS_*p)[2] = grav_field_g[p] * (dir == _YDIR_);
126  (f+_MODEL_NVARS_*p)[3] = grav_field_g[p];
127  }
128  return;
129 }
130 
131 /*! Computes the gravitational source term using a well-balanced formulation: the source term
132  is rewritten as follows:
133  \f{equation}{
134  \left[\begin{array}{c} 0 \\ -\rho {\bf g}\cdot{\bf \hat{i}} \\ -\rho {\bf g}\cdot{\bf \hat{j}} \\ -\rho u {\bf g}\cdot{\bf \hat{i}} - \rho v {\bf g}\cdot{\bf \hat{j}} \end{array}\right]
135  =
136  \left[\begin{array}{cccc} 0 & p_0 \varrho^{-1} & 0 & p_0 u \varrho^{-1} \end{array}\right] \cdot \frac{\partial}{\partial x}\left[\begin{array}{c} 0 \\ \varphi \\ 0 \\ \varphi \end{array}\right]
137  +
138  \left[\begin{array}{cccc} 0 & 0 & p_0 \varrho^{-1} & p_0 v \varrho^{-1} \end{array}\right] \cdot \frac{\partial}{\partial y}\left[\begin{array}{c} 0 \\ 0 \\ \varphi \\ \varphi \end{array}\right]
139  \f}
140  where \f$\varphi = \varphi\left(x,y\right)\f$ and \f$\varrho = \varrho\left(x,y\right)\f$ are computed in
141  NavierStokes2DGravityField(). The derivatives are computed in an exactly identical manner as the hyperbolic
142  flux (i.e., using a conservative finite-difference formulation) (see HyperbolicFunction()).
143  \n\n
144  References:
145  + Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source
146  Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889,
147  7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX,
148  http://dx.doi.org/10.2514/6.2015-2889
149  + Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm
150  for Atmospheric Flows, AIAA Journal, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580.
151 */
152 extern "C" int gpuNavierStokes2DSource(
153  double *source, /*!< Array to hold the computed source */
154  double *gpu_u, /*!< Solution vector array */
155  void *s, /*!< Solver object of type #HyPar */
156  void *m, /*!< MPI object of type #MPIVariables */
157  double t /*!< Current simulation time */
158 )
159 {
160  HyPar *solver = (HyPar* ) s;
161  MPIVariables *mpi = (MPIVariables*) m;
162  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
163 
164  if ((param->grav_x == 0.0) && (param->grav_y == 0.0))
165  return(0); /* no gravitational forces */
166 
167  double *SourceI = solver->fluxI; /* interace source term */
168  double *SourceC = solver->fluxC; /* cell-centered source term */
169  double *SourceL = solver->fL;
170  double *SourceR = solver->fR;
171 
172  int ndims = solver->ndims;
173  int ghosts = solver->ghosts;
174  int *dim = solver->dim_local;
175  double *gpu_x = solver->gpu_x;
176  double *gpu_dxinv = solver->gpu_dxinv;
177  double RT = param->p0 / param->rho0;
178 
179  /* Along X-direction */
180  if (param->grav_x != 0.0) {
181  printf("ALONG X-DIRECTION not tested yet\n");
182  exit(0);
183 
184  /* calculate the split source function exp(-phi/RT) */
185  gpuNavierStokes2DSourceFunction(SourceC,gpu_u,gpu_x,solver,mpi,t,_XDIR_);
186  /* calculate the left and right interface source terms */
187  solver->InterpolateInterfacesHyp(SourceL,SourceC,gpu_u,gpu_x, 1,_XDIR_,solver,mpi,0);
188  solver->InterpolateInterfacesHyp(SourceR,SourceC,gpu_u,gpu_x,-1,_XDIR_,solver,mpi,0);
189  /* calculate the final interface source term */
190  gpuNavierStokes2DSourceUpwind(SourceI,SourceL,SourceR,gpu_u,_XDIR_,solver,t);
191  /* calculate the final cell-centered source term */
192 
193  int ngrid_points = 1; for (int i = 0; i < ndims; i++) ngrid_points *= dim[i];
194  int nblocks = (ngrid_points - 1) / GPU_THREADS_PER_BLOCK + 1;
195  NavierStokes2DSource_Xdir_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
196  ngrid_points, ghosts, ndims, param->gamma, RT,
197  solver->gpu_dim_local, gpu_dxinv, param->gpu_grav_field_f, SourceI, gpu_u, source
198  );
199  }
200 
201  /* Along Y-direction */
202  if (param->grav_y != 0.0) {
203  printf("ALONG Y-DIRECTION\n");
204  /* set interface dimensions */
205  /*_ArrayCopy1D_(dim,dim_interface,ndims); dim_interface[_YDIR_]++;*/
206  /* calculate the split source function exp(-phi/RT) */
207  gpuNavierStokes2DSourceFunction(SourceC,gpu_u,gpu_x,solver,mpi,t,_YDIR_);
208  /* calculate the left and right interface source terms */
209  solver->InterpolateInterfacesHyp(SourceL,SourceC,gpu_u,gpu_x, 1,_YDIR_,solver,mpi,0);
210  solver->InterpolateInterfacesHyp(SourceR,SourceC,gpu_u,gpu_x,-1,_YDIR_,solver,mpi,0);
211  /* calculate the final interface source term */
212  gpuNavierStokes2DSourceUpwind(SourceI,SourceL,SourceR,gpu_u,_YDIR_,solver,t);
213 
214  int ngrid_points = 1; for (int i = 0; i < ndims; i++) ngrid_points *= dim[i];
215  int nblocks = (ngrid_points - 1) / GPU_THREADS_PER_BLOCK + 1;
216  NavierStokes2DSource_Ydir_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
217  ngrid_points, ghosts, ndims, param->gamma, RT,
218  solver->gpu_dim_local, gpu_dxinv, param->gpu_grav_field_f, SourceI, gpu_u, source
219  );
220  }
221  return(0);
222 }
223 
224 /*! Compute the source function in the well-balanced treatment of the source terms. The source
225  function is:
226  \f{equation}{
227  dir = x \rightarrow \left[\begin{array}{c}0 \\ \varphi\left(x,y\right) \\ 0 \\ \varphi\left(x,y\right) \end{array}\right],
228  \ dir = y \rightarrow \left[\begin{array}{c}0 \\ 0 \\ \varphi\left(x,y\right) \\ \varphi\left(x,y\right) \end{array}\right]
229  \f}
230  where \f$\varphi\f$ (#NavierStokes2D::grav_field_g) is computed in NavierStokes2D::GravityField().
231  \n\n
232  References:
233  + Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source
234  Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889,
235  7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX,
236  http://dx.doi.org/10.2514/6.2015-2889
237  + Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm
238  for Atmospheric Flows, Submitted
239 */
240 int gpuNavierStokes2DSourceFunction(
241  double *gpu_f, /*!< Array to hold the computed source function */
242  const double *gpu_u, /*!< Solution vector array */
243  const double *gpu_x, /*!< Array of spatial coordinates (grid) */
244  void *s, /*!< Solver object of type #HyPar */
245  void *m, /*!< MPI object of type #MPIVariables */
246  double t, /*!< Current simulation time */
247  int dir /*!< Spatial dimension (x or y) */
248 )
249 {
250  HyPar *solver = (HyPar* ) s;
251  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
252 
253  int ghosts = solver->ghosts;
254  int *dim = solver->dim_local;
255  int ndims = solver->ndims;
256 
257  double cpu_time = 0.0;
258  clock_t cpu_start, cpu_end;
259 
260  int ngrid_points = 1; for (int i = 0; i < ndims; i++) ngrid_points *= (dim[i] + 2*ghosts);
261  int nblocks = (ngrid_points - 1) / GPU_THREADS_PER_BLOCK + 1;
262 
263  cpu_start = clock();
264  NavierStokes2DSourceFunction_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
265  ngrid_points, dir, param->gpu_grav_field_g, gpu_f);
266  cudaDeviceSynchronize();
267  cpu_end = clock();
268  cpu_time += (double)(cpu_end - cpu_start) / CLOCKS_PER_SEC;
269 
270  return(0);
271 }
272 
273 /*! Compute the "upwind" source function value at the interface: the upwinding is just the
274  arithmetic average of the left and right biased interpolated values of the source function
275  at the interface.
276  \n\n
277  References:
278  + Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source
279  Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889,
280  7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX,
281  http://dx.doi.org/10.2514/6.2015-2889
282  + Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm
283  for Atmospheric Flows, Submitted
284 */
285 int gpuNavierStokes2DSourceUpwind(
286  double *gpu_fI, /*!< Array to hold the computed "upwind" interface source function */
287  const double *gpu_fL, /*!< Interface source function value computed using left-biased interpolation */
288  const double *gpu_fR, /*!< Interface source function value computed using right-biased interpolation */
289  const double *gpu_u, /*!< Solution vector array */
290  int dir, /*!< Spatial dimension (x or y) */
291  void *s, /*!< Solver object of type #HyPar */
292  double t /*!< Current simulation time */
293 )
294 {
295  HyPar *solver = (HyPar*) s;
296  _DECLARE_IERR_;
297 
298  int ndims = solver->ndims;
299  int *dim = solver->dim_local;
300  int bounds_inter[ndims];
301 
302  double cpu_time = 0.0;
303  clock_t cpu_start, cpu_end;
304 
305  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
306 
307  int ngrid_points; _ArrayProduct1D_(bounds_inter,ndims,ngrid_points);
308  int nblocks = (ngrid_points - 1) / GPU_THREADS_PER_BLOCK + 1;
309 
310  cpu_start = clock();
311  NavierStokes2DSourceUpwind_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
312  ngrid_points, gpu_fL, gpu_fR, gpu_fI);
313  cudaDeviceSynchronize();
314  cpu_end = clock();
315  cpu_time += (double)(cpu_end - cpu_start) / CLOCKS_PER_SEC;
316 
317  return(0);
318 }
319 
320 #endif