HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions | Variables
NavierStokes3DUpwind_GPU.cu File Reference

Contains functions to compute the upwind flux at grid interfaces for the 3D Navier Stokes equations. More...

#include <math.h>
#include <basic_gpu.h>
#include <arrayfunctions_gpu.h>
#include <mathfunctions.h>
#include <matmult_native.h>
#include <physicalmodels/navierstokes3d.h>
#include <hypar.h>

Go to the source code of this file.

Functions

__global__ void gpuNavierStokes3DUpwindRusanov_kernel (int npoints_grid, int npoints_local_wghosts, int dir, int ghosts, double gamma, const int *__restrict__ dim, const double *__restrict__ grav_field_g, const double *__restrict__ fL, const double *__restrict__ fR, const double *__restrict__ uL, const double *__restrict__ uR, const double *__restrict__ u, double *__restrict__ fI)
 
int gpuNavierStokes3DUpwindRusanov (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
__global__ void gpuNavierStokes3DUpwindRoe_kernel (int npoints_grid, int npoints_local_wghosts, int dir, int ghosts, double gamma, const int *__restrict__ dim, const double *__restrict__ grav_field_g, const double *__restrict__ fL, const double *__restrict__ fR, const double *__restrict__ uL, const double *__restrict__ uR, const double *__restrict__ u, double *__restrict__ fI)
 
int gpuNavierStokes3DUpwindRoe (double *__restrict__ fI, double *__restrict__ fL, double *__restrict__ fR, double *__restrict__ uL, double *__restrict__ uR, double *__restrict__ u, int dir, void *s, double t)
 

Variables

static const int dummy = 1
 

Detailed Description

Contains functions to compute the upwind flux at grid interfaces for the 3D Navier Stokes equations.

Author
Youngdae Kim

Definition in file NavierStokes3DUpwind_GPU.cu.

Function Documentation

__global__ void gpuNavierStokes3DUpwindRusanov_kernel ( int  npoints_grid,
int  npoints_local_wghosts,
int  dir,
int  ghosts,
double  gamma,
const int *__restrict__  dim,
const double *__restrict__  grav_field_g,
const double *__restrict__  fL,
const double *__restrict__  fR,
const double *__restrict__  uL,
const double *__restrict__  uR,
const double *__restrict__  u,
double *__restrict__  fI 
)

Definition at line 257 of file NavierStokes3DUpwind_GPU.cu.

272 {
273  int p = blockDim.x * blockIdx.x + threadIdx.x;
274 
275  if (p < npoints_grid) {
276  double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_];
277  int bounds_inter[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_],
278  indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_];
279 
280  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir] += 1;
281  _ArrayIndexnD_(_MODEL_NDIMS_,p,bounds_inter,index_inter,0);
282  _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
283  _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
284  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,ghosts,pL);
285  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,ghosts,pR);
286 
287  /* Modified Rusanov's upwinding scheme */
288 
289  udiff[0] = 0.5 * (uR[p+ 0] - uL[p+ 0]);
290  udiff[1] = 0.5 * (uR[p+ npoints_grid] - uL[p+ npoints_grid]);
291  udiff[2] = 0.5 * (uR[p+2*npoints_grid] - uL[p+2*npoints_grid]);
292  udiff[3] = 0.5 * (uR[p+3*npoints_grid] - uL[p+3*npoints_grid]);
293  udiff[4] = 0.5 * (uR[p+4*npoints_grid] - uL[p+4*npoints_grid]);
294 
295  _NavierStokes3DRoeAverage_(uavg,npoints_local_wghosts,(u+pL),(u+pR),gamma);
296 
297  double c, vel[_MODEL_NDIMS_], rho,E,P;
298  _NavierStokes3DGetFlowVar_((u+pL),npoints_local_wghosts,rho,vel[0],vel[1],vel[2],E,P,gamma);
299  c = sqrt(gamma*P/rho);
300  double alphaL = c + absolute(vel[dir]);
301  _NavierStokes3DGetFlowVar_((u+pR),npoints_local_wghosts,rho,vel[0],vel[1],vel[2],E,P,gamma);
302  c = sqrt(gamma*P/rho);
303  double alphaR = c + absolute(vel[dir]);
304  _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,gamma);
305  c = sqrt(gamma*P/rho);
306  double alphaavg = c + absolute(vel[dir]);
307 
308  double kappa = max(grav_field_g[pL],grav_field_g[pR]);
309  double alpha = kappa*max3(alphaL,alphaR,alphaavg);
310 
311  fI[p+0] = 0.5*(fL[p+ 0]+fR[p+ 0])-alpha*udiff[0];
312  fI[p+ npoints_grid] = 0.5*(fL[p+ npoints_grid]+fR[p+ npoints_grid])-alpha*udiff[1];
313  fI[p+2*npoints_grid] = 0.5*(fL[p+2*npoints_grid]+fR[p+2*npoints_grid])-alpha*udiff[2];
314  fI[p+3*npoints_grid] = 0.5*(fL[p+3*npoints_grid]+fR[p+3*npoints_grid])-alpha*udiff[3];
315  fI[p+4*npoints_grid] = 0.5*(fL[p+4*npoints_grid]+fR[p+4*npoints_grid])-alpha*udiff[4];
316  }
317  return;
318 }
_NavierStokes3DGetFlowVar_
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)
Definition: navierstokes3d.h:98
max3
#define max3(a, b, c)
Definition: math_ops.h:27
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayIndexnD_
#define _ArrayIndexnD_(N, index, imax, i, ghost)
Definition: arrayfunctions.h:21
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
absolute
#define absolute(a)
Definition: math_ops.h:32
_NavierStokes3DRoeAverage_
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
Definition: navierstokes3d.h:144
max
#define max(a, b)
Definition: math_ops.h:18
dummy
static const int dummy
Definition: NavierStokes3DUpwind_GPU.cu:13
int gpuNavierStokes3DUpwindRusanov ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

Rusanov's upwinding scheme.

\begin{equation} {\bf f}_{j+1/2} = \frac{1}{2}\left[ {\bf f}_{j+1/2}^L + {\bf f}_{j+1/2}^R - \max_{j,j+1} \nu_j \left( {\bf u}_{j+1/2}^R - {\bf u}_{j+1/2}^L \right)\right] \end{equation}

where \(\nu = c + \left|u\right|\).

  • Rusanov, V. V., "The calculation of the interaction of non-stationary shock waves and obstacles," USSR Computational Mathematics and Mathematical Physics, Vol. 1, No. 2, 1962, pp. 304–320

This upwinding scheme is modified for the balanced discretization of the 3D Navier Stokes equations when there is a non-zero gravitational force. See the reference below. For flows without any gravitational forces, it reduces to its original form.

  • Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889, 7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX, http://dx.doi.org/10.2514/6.2015-2889
  • Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm for Atmospheric Flows, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x,y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 340 of file NavierStokes3DUpwind_GPU.cu.

351 {
352  HyPar *solver = (HyPar*) s;
353  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
354  int *dim = solver->dim_local;
355 
356  int bounds_inter[_MODEL_NDIMS_];
357  _ArrayCopy1D_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
358  int npoints_grid; _ArrayProduct1D_(bounds_inter,_MODEL_NDIMS_,npoints_grid);
359  int nblocks = (npoints_grid-1)/GPU_THREADS_PER_BLOCK + 1;
360 
361  gpuNavierStokes3DUpwindRusanov_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
362  npoints_grid, solver->npoints_local_wghosts, dir, solver->ghosts, param->gamma, solver->gpu_dim_local,
363  param->gpu_grav_field_g, fL, fR, uL, uR, u, fI
364  );
365  cudaDeviceSynchronize();
366 
367  return 0;
368 }
HyPar::npoints_local_wghosts
int npoints_local_wghosts
Definition: hypar.h:42
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
GPU_THREADS_PER_BLOCK
#define GPU_THREADS_PER_BLOCK
Definition: basic_gpu.h:7
HyPar::ghosts
int ghosts
Definition: hypar.h:52
HyPar::gpu_dim_local
int * gpu_dim_local
Definition: hypar.h:455
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
NavierStokes3D
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
Definition: navierstokes3d.h:482
_ArrayProduct1D_
#define _ArrayProduct1D_(x, size, p)
Definition: arrayfunctions.h:377
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
NavierStokes3D::gpu_grav_field_g
double * gpu_grav_field_g
Definition: navierstokes3d.h:546
__global__ void gpuNavierStokes3DUpwindRoe_kernel ( int  npoints_grid,
int  npoints_local_wghosts,
int  dir,
int  ghosts,
double  gamma,
const int *__restrict__  dim,
const double *__restrict__  grav_field_g,
const double *__restrict__  fL,
const double *__restrict__  fR,
const double *__restrict__  uL,
const double *__restrict__  uR,
const double *__restrict__  u,
double *__restrict__  fI 
)

Definition at line 372 of file NavierStokes3DUpwind_GPU.cu.

387 {
388  int p = blockDim.x * blockIdx.x + threadIdx.x;
389 
390  if (p < npoints_grid) {
391  double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
392  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
393  modA[_MODEL_NVARS_*_MODEL_NVARS_];
394  double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
395  int bounds_inter[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_],
396  indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_];
397 
398  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir] += 1;
399  _ArrayIndexnD_(_MODEL_NDIMS_,p,bounds_inter,index_inter,0);
400  _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
401  _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
402  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,ghosts,pL);
403  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,ghosts,pR);
404 
405  /* Modified Roe's upwinding scheme */
406 
407  udiff[0] = 0.5 * (uR[p+0] - uL[p+0]);
408  udiff[1] = 0.5 * (uR[p+ npoints_grid] - uL[p+ npoints_grid]);
409  udiff[2] = 0.5 * (uR[p+2*npoints_grid] - uL[p+2*npoints_grid]);
410  udiff[3] = 0.5 * (uR[p+3*npoints_grid] - uL[p+3*npoints_grid]);
411  udiff[4] = 0.5 * (uR[p+4*npoints_grid] - uL[p+4*npoints_grid]);
412 
413  _NavierStokes3DRoeAverage_ (uavg,npoints_local_wghosts,(u+pL),(u+pR),gamma);
414  _NavierStokes3DEigenvalues_ (uavg,dummy,D,gamma,dir);
415  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,gamma,dir);
416  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,gamma,dir);
417 
418  /* Harten's Entropy Fix - Page 362 of Leveque */
419  int k;
420  double delta = 0.000001, delta2 = delta*delta;
421  k=0; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
422  k=6; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
423  k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
424  k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
425  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
426 
427  MatMult5(_MODEL_NVARS_,DL,D,L);
428  MatMult5(_MODEL_NVARS_,modA,R,DL);
429  MatVecMult5(_MODEL_NVARS_,udiss,modA,udiff);
430 
431  fI[p+0] = 0.5 * (fL[p+0]+fR[p+0]) - udiss[0];
432  fI[p+ npoints_grid] = 0.5 * (fL[p+ npoints_grid]+fR[p+ npoints_grid]) - udiss[1];
433  fI[p+2*npoints_grid] = 0.5 * (fL[p+2*npoints_grid]+fR[p+2*npoints_grid]) - udiss[2];
434  fI[p+3*npoints_grid] = 0.5 * (fL[p+3*npoints_grid]+fR[p+3*npoints_grid]) - udiss[3];
435  fI[p+4*npoints_grid] = 0.5 * (fL[p+4*npoints_grid]+fR[p+4*npoints_grid]) - udiss[4];
436  }
437  return;
438 }
_NavierStokes3DLeftEigenvectors_
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
Definition: navierstokes3d.h:288
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayIndexnD_
#define _ArrayIndexnD_(N, index, imax, i, ghost)
Definition: arrayfunctions.h:21
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
MatVecMult5
#define MatVecMult5(N, y, A, x)
Definition: matmult_native.h:158
_NavierStokes3DEigenvalues_
#define _NavierStokes3DEigenvalues_(u, stride, D, gamma, dir)
Definition: navierstokes3d.h:254
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
MatMult5
#define MatMult5(N, A, X, Y)
Definition: matmult_native.h:82
_NavierStokes3DRightEigenvectors_
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
Definition: navierstokes3d.h:384
absolute
#define absolute(a)
Definition: math_ops.h:32
_NavierStokes3DRoeAverage_
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
Definition: navierstokes3d.h:144
dummy
static const int dummy
Definition: NavierStokes3DUpwind_GPU.cu:13
int gpuNavierStokes3DUpwindRoe ( double *__restrict__  fI,
double *__restrict__  fL,
double *__restrict__  fR,
double *__restrict__  uL,
double *__restrict__  uR,
double *__restrict__  u,
int  dir,
void *  s,
double  t 
)

Roe's upwinding scheme.

\begin{equation} {\bf f}_{j+1/2} = \frac{1}{2}\left[ {\bf f}_{j+1/2}^L + {\bf f}_{j+1/2}^R - \left| A\left({\bf u}_{j+1/2}^L,{\bf u}_{j+1/2}^R\right) \right| \left( {\bf u}_{j+1/2}^R - {\bf u}_{j+1/2}^L \right)\right] \end{equation}

This upwinding scheme is modified for the balanced discretization of the 3D Navier Stokes equations when there is a non-zero gravitational force. See the reference below. For flows without any gravitational forces, it reduces to its original form.

  • Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889, 7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX, http://dx.doi.org/10.2514/6.2015-2889
  • Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm for Atmospheric Flows, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x,y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 460 of file NavierStokes3DUpwind_GPU.cu.

471 {
472  HyPar *solver = (HyPar*) s;
473  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
474  int *dim = solver->dim_local;
475 
476  int bounds_inter[_MODEL_NDIMS_];
477  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
478  int npoints_grid; _ArrayProduct1D_(bounds_inter,_MODEL_NDIMS_,npoints_grid);
479  int nblocks = (npoints_grid-1)/GPU_THREADS_PER_BLOCK + 1;
480 
481 #if defined(GPU_STAT)
482  cudaEvent_t startEvent, stopEvent;
483  float milliseconds = 0;
484 
485  checkCuda( cudaEventCreate(&startEvent) );
486  checkCuda( cudaEventCreate(&stopEvent) );
487  checkCuda( cudaEventRecord(startEvent, 0) );
488 #endif
489 
490  gpuNavierStokes3DUpwindRoe_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
491  npoints_grid, solver->npoints_local_wghosts, dir, solver->ghosts, param->gamma, solver->gpu_dim_local,
492  param->gpu_grav_field_g, fL, fR, uL, uR, u, fI
493  );
494 
495 #if defined(GPU_STAT)
496  checkCuda( cudaEventRecord(stopEvent, 0) );
497  checkCuda( cudaEventSynchronize(stopEvent) );
498 #endif
499 
500  cudaDeviceSynchronize();
501 
502 #if defined(GPU_STAT)
503  checkCuda( cudaEventElapsedTime(&milliseconds, startEvent, stopEvent) );
504  printf("%-50s GPU time (secs) = %.6f dir = %d\n",
505  "NavierStokes3DUpwindRoe2", milliseconds*1e-3, dir);
506  checkCuda( cudaEventDestroy(startEvent) );
507  checkCuda( cudaEventDestroy(stopEvent) );
508 #endif
509 
510 
511 
512 
513  return 0;
514 }
HyPar::npoints_local_wghosts
int npoints_local_wghosts
Definition: hypar.h:42
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayCopy1D3_
#define _ArrayCopy1D3_(x, y, size)
Definition: arrayfunctions.h:343
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
GPU_THREADS_PER_BLOCK
#define GPU_THREADS_PER_BLOCK
Definition: basic_gpu.h:7
HyPar::ghosts
int ghosts
Definition: hypar.h:52
HyPar::gpu_dim_local
int * gpu_dim_local
Definition: hypar.h:455
NavierStokes3D
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
Definition: navierstokes3d.h:482
_ArrayProduct1D_
#define _ArrayProduct1D_(x, size, p)
Definition: arrayfunctions.h:377
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
NavierStokes3D::gpu_grav_field_g
double * gpu_grav_field_g
Definition: navierstokes3d.h:546

Variable Documentation

const int dummy = 1
static

Definition at line 13 of file NavierStokes3DUpwind_GPU.cu.