HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
NavierStokes3DModifiedSolution_GPU.cu
Go to the documentation of this file.
1 /*! @file NavierStokes3DModifiedSolution_GPU.cu
2  @author Youngdae Kim
3  @brief Compute the modified solution for the well-balanced treatment of gravitational source terms.
4 */
5 #include <basic_gpu.h>
6 #include <arrayfunctions_gpu.h>
7 #include <physicalmodels/navierstokes3d.h>
8 #include <mpivars.h>
9 #include <hypar.h>
10 
11 #ifdef CUDA_VAR_ORDERDING_AOS
12 
13 /*! Kernel for gpuNavierStokes3DModifiedSolution() */
14 __global__
15 void gpuNavierStokes3DModifiedSolution_kernel(
16  int npoints_grid,
17  double gamma,
18  double inv_gamma_m1,
19  const double * __restrict__ grav_field_f,
20  const double * __restrict__ grav_field_g,
21  const double * __restrict__ u,
22  double * __restrict__ uC
23 )
24 {
25  int p = blockDim.x * blockIdx.x + threadIdx.x;
26 
27  if (p < npoints_grid) {
28  double rho, uvel, vvel, wvel, E, P;
29 
30  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*p),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,E,P,gamma);
31  uC[_MODEL_NVARS_*p+0] = u[_MODEL_NVARS_*p+0] * grav_field_f[p];
32  uC[_MODEL_NVARS_*p+1] = u[_MODEL_NVARS_*p+1] * grav_field_f[p];
33  uC[_MODEL_NVARS_*p+2] = u[_MODEL_NVARS_*p+2] * grav_field_f[p];
34  uC[_MODEL_NVARS_*p+3] = u[_MODEL_NVARS_*p+3] * grav_field_f[p];
35  uC[_MODEL_NVARS_*p+4] = (P*inv_gamma_m1)*(1.0/grav_field_g[p]) + (0.5*rho*(uvel*uvel+vvel*vvel+wvel*wvel))*grav_field_f[p];
36  }
37  return;
38 }
39 
40 /*!
41  This function computes the modified solution for the well-balanced treatment of the
42  gravitational source terms. The modified solution vector is given by
43  \f{equation}{
44  {\bf u}^* = \left[\begin{array}{c} \rho \varrho^{-1}\left(x,y\right) \\ \rho u \varrho^{-1}\left(x,y\right) \\ \rho v \varrho^{-1}\left(x,y\right) \\ \rho w \varrho^{-1}\left(x,y\right) \\ e^* \end{array}\right]
45  \f}
46  where
47  \f{equation}{
48  e^* = \frac{p \varphi^{-1}\left(x,y\right)}{\gamma-1} + \frac{1}{2}\rho \varrho^{-1}\left(x,y\right) \left(u^2+v^2+w^2\right)
49  \f}
50  \f$\varrho\f$ and \f$\varphi\f$ are computed in #NavierStokes3DGravityField(). For flows without gravity, \f${\bf u}^* = {\bf u}\f$.
51 
52  References:
53  + Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source
54  Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889,
55  7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX,
56  http://dx.doi.org/10.2514/6.2015-2889
57  + Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm
58  for Atmospheric Flows, AIAA Journal, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580.
59 */
60 extern "C" int gpuNavierStokes3DModifiedSolution(
61  double * __restrict__ uC, /*!< Array to hold the computed modified solution */
62  double * __restrict__ u, /*!< Solution vector array */
63  int d, /*!< spatial dimension (not used) */
64  void * __restrict__ s, /*!< Solver object of type #HyPar */
65  void * __restrict__ m, /*!< MPI object of time #MPIVariables */
66  double waqt /*!< Current simulation time */
67 )
68 {
69  HyPar *solver = (HyPar*) s;
70  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
71 
72  double inv_gamma_m1 = 1.0 / (param->gamma-1.0);
73  int nblocks = (solver->npoints_local_wghosts-1)/GPU_THREADS_PER_BLOCK + 1;
74 
75  gpuNavierStokes3DModifiedSolution_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
76  solver->npoints_local_wghosts, param->gamma, inv_gamma_m1,
77  param->gpu_grav_field_f, param->gpu_grav_field_g, u, uC
78  );
79  cudaDeviceSynchronize();
80 
81  return 0;
82 }
83 
84 #else
85 
86 /*! Kernel for gpuNavierStokes3DModifiedSolution() */
87 __global__
88 void gpuNavierStokes3DModifiedSolution_kernel(
89  int npoints_grid,
90  double gamma,
91  double inv_gamma_m1,
92  const double * __restrict__ grav_field_f,
93  const double * __restrict__ grav_field_g,
94  const double * __restrict__ u,
95  double * __restrict__ uC
96 )
97 {
98  int p = blockDim.x * blockIdx.x + threadIdx.x;
99 
100  if (p < npoints_grid) {
101  double rho, uvel, vvel, wvel, E, P;
102 
103  _NavierStokes3DGetFlowVar_((u+p),npoints_grid,rho,uvel,vvel,wvel,E,P,gamma);
104  uC[p+0] = u[p+0] * grav_field_f[p];
105  uC[p+ npoints_grid] = u[p+ npoints_grid] * grav_field_f[p];
106  uC[p+2*npoints_grid] = u[p+2*npoints_grid] * grav_field_f[p];
107  uC[p+3*npoints_grid] = u[p+3*npoints_grid] * grav_field_f[p];
108  uC[p+4*npoints_grid] = (P*inv_gamma_m1)*(1.0/grav_field_g[p]) + (0.5*rho*(uvel*uvel+vvel*vvel+wvel*wvel))*grav_field_f[p];
109  }
110  return;
111 }
112 
113 /*!
114  This function computes the modified solution for the well-balanced treatment of the
115  gravitational source terms. The modified solution vector is given by
116  \f{equation}{
117  {\bf u}^* = \left[\begin{array}{c} \rho \varrho^{-1}\left(x,y\right) \\ \rho u \varrho^{-1}\left(x,y\right) \\ \rho v \varrho^{-1}\left(x,y\right) \\ \rho w \varrho^{-1}\left(x,y\right) \\ e^* \end{array}\right]
118  \f}
119  where
120  \f{equation}{
121  e^* = \frac{p \varphi^{-1}\left(x,y\right)}{\gamma-1} + \frac{1}{2}\rho \varrho^{-1}\left(x,y\right) \left(u^2+v^2+w^2\right)
122  \f}
123  \f$\varrho\f$ and \f$\varphi\f$ are computed in #NavierStokes3DGravityField(). For flows without gravity, \f${\bf u}^* = {\bf u}\f$.
124 
125  References:
126  + Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source
127  Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889,
128  7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX,
129  http://dx.doi.org/10.2514/6.2015-2889
130  + Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm
131  for Atmospheric Flows, AIAA Journal, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580.
132 */
133 extern "C" int gpuNavierStokes3DModifiedSolution(
134  double * __restrict__ uC, /*!< Array to hold the computed modified solution */
135  double * __restrict__ u, /*!< Solution vector array */
136  int d, /*!< spatial dimension (not used) */
137  void * __restrict__ s, /*!< Solver object of type #HyPar */
138  void * __restrict__ m, /*!< MPI object of time #MPIVariables */
139  double waqt /*!< Current simulation time */
140 )
141 {
142  HyPar *solver = (HyPar*) s;
143  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
144 
145  double inv_gamma_m1 = 1.0 / (param->gamma-1.0);
146  int nblocks = (solver->npoints_local_wghosts-1)/GPU_THREADS_PER_BLOCK + 1;
147 
148  gpuNavierStokes3DModifiedSolution_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
149  solver->npoints_local_wghosts, param->gamma, inv_gamma_m1,
150  param->gpu_grav_field_f, param->gpu_grav_field_g, u, uC
151  );
152  cudaDeviceSynchronize();
153 
154  return 0;
155 }
156 #endif
157