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

Definitions for the functions computing the first derivative. More...

Go to the source code of this file.

Macros

#define _FIRST_ORDER_   "1"
 
#define _SECOND_ORDER_CENTRAL_   "2"
 
#define _FOURTH_ORDER_CENTRAL_   "4"
 

Functions

int FirstDerivativeFirstOrder (double *, double *, int, int, void *, void *)
 
int FirstDerivativeSecondOrderCentral (double *, double *, int, int, void *, void *)
 
int FirstDerivativeFourthOrderCentral (double *, double *, int, int, void *, void *)
 
int FirstDerivativeSecondOrderCentralNoGhosts (double *, double *, int, int, int, int *, int, int, void *)
 
int gpuFirstDerivativeFourthOrderCentral (double *, double *, int, int, void *, void *)
 

Detailed Description

Definitions for the functions computing the first derivative.

Author
Debojyoti Ghosh

Definition in file firstderivative.h.

Macro Definition Documentation

#define _FIRST_ORDER_   "1"

First order scheme

Definition at line 10 of file firstderivative.h.

#define _SECOND_ORDER_CENTRAL_   "2"

Second order central scheme

Definition at line 12 of file firstderivative.h.

#define _FOURTH_ORDER_CENTRAL_   "4"

Fourth order central scheme

Definition at line 14 of file firstderivative.h.

Function Documentation

int FirstDerivativeFirstOrder ( double *  Df,
double *  f,
int  dir,
int  bias,
void *  s,
void *  m 
)

First order approximation to the first derivative (note: not divided by grid spacing)

Computes the first-order finite-difference approximation to the first derivative (Note: not divided by the grid spacing):

\begin{equation} \left(\partial f\right)_i = \left\{ \begin{array}{ll} f_{i+1} - f_i & {\rm bias} = 1 \\ f_i - f_{i-1} & {\rm bias} = -1 \end{array}\right. \end{equation}

where \(i\) is the grid index along the spatial dimension of the derivative.

Notes:

  • The first derivative is computed at the grid points or the cell centers.
  • The first derivative is computed at the ghost points too. Thus, biased schemes are used at and near the boundaries.
  • Df and f are 1D arrays containing the function and its computed derivatives on a multi- dimensional grid. The derivative along the specified dimension dir is computed by looping through all grid lines along dir.
Parameters
DfArray to hold the computed first derivative (with ghost points)
fArray containing the grid point function values whose first derivative is to be computed (with ghost points)
dirThe spatial dimension along which the derivative is computed
biasForward or backward differencing for non-central finite-difference schemes (-1: backward, 1: forward)
sSolver object of type SolverContext
mMPI object of type MPIContext

Definition at line 37 of file FirstDerivativeFirstOrder.c.

47 {
48  SolverContext *solver = (SolverContext*) s;
49  int i, j, v;
50 
51  int ghosts = solver->ghosts;
52  int ndims = solver->ndims;
53  int nvars = solver->nvars;
54  int *dim = solver->dim_local;
55 
56 
57  if ((!Df) || (!f)) {
58  fprintf(stderr, "Error in FirstDerivativeSecondOrder(): input arrays not allocated.\n");
59  return(1);
60  }
61 
62  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
63  dimension "dir" */
64  int indexC[ndims], index_outer[ndims], bounds_outer[ndims];
65  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
66  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
67 
68 #pragma omp parallel for schedule(auto) default(shared) private(i,j,v,index_outer,indexC)
69  for (j=0; j<N_outer; j++) {
70  _ArrayIndexnD_(ndims,j,bounds_outer,index_outer,0);
71  _ArrayCopy1D_(index_outer,indexC,ndims);
72  /* left boundary */
73  for (i = -ghosts; i < -ghosts+1; i++) {
74  int qC, qR;
75  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
76  indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR );
77  for (v=0; v<nvars; v++) Df[qC*nvars+v] = f[qR*nvars+v]-f[qC*nvars+v];
78  }
79  /* interior */
80  for (i = -ghosts+1; i < dim[dir]+ghosts-1; i++) {
81  int qC, qL, qR;
82  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
83  indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL);
84  indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR);
85  for (v=0; v<nvars; v++) Df[qC*nvars+v] = max(bias,0)*f[qR*nvars+v]-bias*f[qC*nvars+v]+min(bias,0)*f[qL*nvars+v];
86  }
87  /* right boundary */
88  for (i = dim[dir]+ghosts-1; i < dim[dir]+ghosts; i++) {
89  int qL, qC;
90  indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL );
91  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
92  for (v=0; v<nvars; v++) Df[qC*nvars+v] = f[qC*nvars+v]-f[qL*nvars+v];
93  }
94  }
95 
96  return(0);
97 }
_ArrayIndexnD_
#define _ArrayIndexnD_(N, index, imax, i, ghost)
Definition: arrayfunctions.h:21
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
HyPar::nvars
int nvars
Definition: hypar.h:29
HyPar::ndims
int ndims
Definition: hypar.h:26
_ArrayProduct1D_
#define _ArrayProduct1D_(x, size, p)
Definition: arrayfunctions.h:377
max
#define max(a, b)
Definition: math_ops.h:18
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
min
#define min(a, b)
Definition: math_ops.h:14
int FirstDerivativeSecondOrderCentral ( double *  Df,
double *  f,
int  dir,
int  bias,
void *  s,
void *  m 
)

Second order approximation to the first derivative (note: not divided by grid spacing)

Computes the second-order finite-difference approximation to the first derivative (Note: not divided by the grid spacing):

\begin{equation} \left(\partial f\right)_i = \left\{ \begin{array}{ll} \frac{1}{2}\left(-3f_i+4f_{i+1}-f_{i+2}\right) & i = -g \\ \frac{1}{2}\left( f_{i+1} - f_{i-1} \right) & -g+1 \leq i \leq N+g-2 \\ \frac{1}{2}\left( f_{i-2} -4f_{i-1}+3f_i \right) & i = N+g-1 \end{array}\right. \end{equation}

where \(i\) is the grid index along the spatial dimension of the derivative, \(g\) is the number of ghost points, and \(N\) is the number of grid points (not including the ghost points) in the spatial dimension of the derivative.

Notes:

  • The first derivative is computed at the grid points or the cell centers.
  • The first derivative is computed at the ghost points too. Thus, biased schemes are used at and near the boundaries.
  • Df and f are 1D arrays containing the function and its computed derivatives on a multi- dimensional grid. The derivative along the specified dimension dir is computed by looping through all grid lines along dir.
Parameters
DfArray to hold the computed first derivative (with ghost points)
fArray containing the grid point function values whose first derivative is to be computed (with ghost points)
dirThe spatial dimension along which the derivative is computed
biasForward or backward differencing for non-central finite-difference schemes (-1: backward, 1: forward)
sSolver object of type SolverContext
mMPI object of type MPIContext

Definition at line 36 of file FirstDerivativeSecondOrder.c.

46 {
47  SolverContext *solver = (SolverContext*) s;
48  int i, j, v;
49 
50  int ghosts = solver->ghosts;
51  int ndims = solver->ndims;
52  int nvars = solver->nvars;
53  int *dim = solver->dim_local;
54 
55 
56  if ((!Df) || (!f)) {
57  fprintf(stderr, "Error in FirstDerivativeSecondOrder(): input arrays not allocated.\n");
58  return(1);
59  }
60 
61  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
62  dimension "dir" */
63  int indexC[ndims], index_outer[ndims], bounds_outer[ndims];
64  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
65  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
66 
67 #pragma omp parallel for schedule(auto) default(shared) private(i,j,v,index_outer,indexC)
68  for (j=0; j<N_outer; j++) {
69  _ArrayIndexnD_(ndims,j,bounds_outer,index_outer,0);
70  _ArrayCopy1D_(index_outer,indexC,ndims);
71  /* left boundary */
72  for (i = -ghosts; i < -ghosts+1; i++) {
73  int qC, qR, qRR;
74  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
75  indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR );
76  indexC[dir] = i+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qRR);
77  for (v=0; v<nvars; v++) Df[qC*nvars+v] = 0.5 * (-3*f[qC*nvars+v]+4*f[qR*nvars+v]-f[qRR*nvars+v]);
78  }
79  /* interior */
80  for (i = -ghosts+1; i < dim[dir]+ghosts-1; i++) {
81  int qC, qL, qR;
82  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
83  indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL);
84  indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR);
85  for (v=0; v<nvars; v++) Df[qC*nvars+v] = 0.5 * (f[qR*nvars+v]-f[qL*nvars+v]);
86  }
87  /* right boundary */
88  for (i = dim[dir]+ghosts-1; i < dim[dir]+ghosts; i++) {
89  int qLL, qL, qC;
90  indexC[dir] = i-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qLL);
91  indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL );
92  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
93  for (v=0; v<nvars; v++) Df[qC*nvars+v] = 0.5 * (3*f[qC*nvars+v]-4*f[qL*nvars+v]+f[qLL*nvars+v]);
94  }
95  }
96 
97  return(0);
98 }
_ArrayIndexnD_
#define _ArrayIndexnD_(N, index, imax, i, ghost)
Definition: arrayfunctions.h:21
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
HyPar::nvars
int nvars
Definition: hypar.h:29
HyPar::ndims
int ndims
Definition: hypar.h:26
_ArrayProduct1D_
#define _ArrayProduct1D_(x, size, p)
Definition: arrayfunctions.h:377
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int FirstDerivativeFourthOrderCentral ( double *  Df,
double *  f,
int  dir,
int  bias,
void *  s,
void *  m 
)

Fourth order approximation to the first derivative (note: not divided by grid spacing)

Computes the fourth-order finite-difference approximation to the first derivative (Note: not divided by the grid spacing):

\begin{equation} \left(\partial f\right)_i = \left\{ \begin{array}{ll} \frac{1}{12}\left(-25f_i+48f_{i+1}-36f_{i+2}+16f_{i+3}-3f_{i+4}\right) & i=-g \\ \frac{1}{12}\left(-3f_{i-1}-10f_i+18f_{i+1}-6f_{i+2}+f_{i+3}\right) & i = -g+1 \\ \frac{1}{2}\left( f_{i-2}-8f_{i-1}+8f_{i+1}-f_{i+2} \right) & -g+2 \leq i \leq N+g-3 \\ \frac{1}{12}\left( -f_{i-3}+6f_{i-2}-18f_{i-1}+10f_i+3f_{i+1}\right) & i = N+g-2 \\ \frac{1}{12}\left( 3f_{i-4}-16f_{i-3}+36f_{i-2}-48f_{i-1}+25f_i \right) & i = N+g-1 \end{array}\right. \end{equation}

where \(i\) is the grid index along the spatial dimension of the derivative, \(g\) is the number of ghost points, and \(N\) is the number of grid points (not including the ghost points) in the spatial dimension of the derivative.

Notes:

  • The first derivative is computed at the grid points or the cell centers.
  • The first derivative is computed at the ghost points too. Thus, biased schemes are used at and near the boundaries.
  • Df and f are 1D arrays containing the function and its computed derivatives on a multi- dimensional grid. The derivative along the specified dimension dir is computed by looping through all grid lines along dir.
Parameters
DfArray to hold the computed first derivative (with ghost points)
fArray containing the grid point function values whose first derivative is to be computed (with ghost points)
dirThe spatial dimension along which the derivative is computed
biasForward or backward differencing for non-central finite-difference schemes (-1: backward, 1: forward)
sSolver object of type SolverContext
mMPI object of type MPIContext

Definition at line 36 of file FirstDerivativeFourthOrder.c.

46 {
47  SolverContext *solver = (SolverContext*) s;
48  int i, j, v;
49 
50  int ghosts = solver->ghosts;
51  int ndims = solver->ndims;
52  int nvars = solver->nvars;
53  int *dim = solver->dim_local;
54 
55 
56  if ((!Df) || (!f)) {
57  fprintf(stderr, "Error in FirstDerivativeFourthOrder(): input arrays not allocated.\n");
58  return(1);
59  }
60 
61  static double one_twelve = 1.0/12.0;
62 
63  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
64  dimension "dir" */
65  int indexC[ndims], index_outer[ndims], bounds_outer[ndims];
66  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
67  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
68 
69 #pragma omp parallel for schedule(auto) default(shared) private(i,j,v,index_outer,indexC)
70  for (j=0; j<N_outer; j++) {
71  _ArrayIndexnD_(ndims,j,bounds_outer,index_outer,0);
72  _ArrayCopy1D_(index_outer,indexC,ndims);
73  /* left boundary */
74  for (i = -ghosts; i < -ghosts+1; i++) {
75  int qC, qp1, qp2, qp3, qp4;
76  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
77  indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
78  indexC[dir] = i+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
79  indexC[dir] = i+3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp3);
80  indexC[dir] = i+4; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp4);
81  for (v=0; v<nvars; v++)
82  Df[qC*nvars+v] = (-25*f[qC*nvars+v]+48*f[qp1*nvars+v]-36*f[qp2*nvars+v]+16*f[qp3*nvars+v]-3*f[qp4*nvars+v])*one_twelve;
83  }
84  for (i = -ghosts+1; i < -ghosts+2; i++) {
85  int qC, qm1, qp1, qp2, qp3;
86  indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
87  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
88  indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
89  indexC[dir] = i+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
90  indexC[dir] = i+3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp3);
91  for (v=0; v<nvars; v++)
92  Df[qC*nvars+v] = (-3*f[qm1*nvars+v]-10*f[qC*nvars+v]+18*f[qp1*nvars+v]-6*f[qp2*nvars+v]+f[qp3*nvars+v])*one_twelve;
93  }
94  /* interior */
95  for (i = -ghosts+2; i < dim[dir]+ghosts-2; i++) {
96  int qC, qm1, qm2, qp1, qp2;
97  indexC[dir] = i-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
98  indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
99  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
100  indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
101  indexC[dir] = i+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
102  for (v=0; v<nvars; v++)
103  Df[qC*nvars+v] = (f[qm2*nvars+v]-8*f[qm1*nvars+v]+8*f[qp1*nvars+v]-f[qp2*nvars+v])*one_twelve;
104  }
105  /* right boundary */
106  for (i = dim[dir]+ghosts-2; i < dim[dir]+ghosts-1; i++) {
107  int qC, qm3, qm2, qm1, qp1;
108  indexC[dir] = i-3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
109  indexC[dir] = i-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
110  indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
111  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
112  indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
113  for (v=0; v<nvars; v++)
114  Df[qC*nvars+v] = (-f[qm3*nvars+v]+6*f[qm2*nvars+v]-18*f[qm1*nvars+v]+10*f[qC*nvars+v]+3*f[qp1*nvars+v])*one_twelve;
115  }
116  for (i = dim[dir]+ghosts-1; i < dim[dir]+ghosts; i++) {
117  int qC, qm4, qm3, qm2, qm1;
118  indexC[dir] = i-4; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm4);
119  indexC[dir] = i-3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
120  indexC[dir] = i-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
121  indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
122  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
123  for (v=0; v<nvars; v++)
124  Df[qC*nvars+v] = (3*f[qm4*nvars+v]-16*f[qm3*nvars+v]+36*f[qm2*nvars+v]-48*f[qm1*nvars+v]+25*f[qC*nvars+v])*one_twelve;
125  }
126  }
127 
128  return(0);
129 }
_ArrayIndexnD_
#define _ArrayIndexnD_(N, index, imax, i, ghost)
Definition: arrayfunctions.h:21
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
HyPar::nvars
int nvars
Definition: hypar.h:29
HyPar::ndims
int ndims
Definition: hypar.h:26
_ArrayProduct1D_
#define _ArrayProduct1D_(x, size, p)
Definition: arrayfunctions.h:377
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int FirstDerivativeSecondOrderCentralNoGhosts ( double *  Df,
double *  f,
int  dir,
int  bias,
int  ndims,
int *  dim,
int  ghosts,
int  nvars,
void *  m 
)

Second order approximation to the first derivative (note: not divided by grid spacing)

Computes the second-order finite-difference approximation to the first derivative (Note: not divided by the grid spacing):

\begin{equation} \left(\partial f\right)_i = \left\{ \begin{array}{ll} \frac{1}{2}\left(-3f_i+4f_{i+1}-f_{i+2}\right) & i = -g \\ \frac{1}{2}\left( f_{i+1} - f_{i-1} \right) & -g+1 \leq i \leq N+g-2 \\ \frac{1}{2}\left( f_{i-2} -4f_{i-1}+3f_i \right) & i = N+g-1 \end{array}\right. \end{equation}

where \(i\) is the grid index along the spatial dimension of the derivative, \(g\) is the number of ghost points, and \(N\) is the number of grid points (not including the ghost points) in the spatial dimension of the derivative.

Notes:

  • The first derivative is computed at the grid points or the cell centers.
  • Df and f are serialized 1D arrays containing the function and its computed derivatives on a multi-dimensional grid. The derivative along the specified dimension dir is computed by looping through all grid lines along dir.
  • Df and f must have the same number of ghost points.
  • Biased stencils are used at the physical boundaries and ghost point values of f are not used. So this function can be used when the ghost values at physical boundaries are not filled. The ghost values at internal (MPI) boundaries are still needed.
Parameters
DfArray to hold the computed first derivative (with ghost points)
fArray containing the grid point function values whose first derivative is to be computed (with ghost points)
dirThe spatial dimension along which the derivative is computed
biasForward or backward differencing for non-central finite-difference schemes (-1: backward, 1: forward)
ndimsNumber of spatial/coordinate dimensions
dimLocal dimensions
ghostsNumber of ghost points
nvarsNumber of vector components at each grid points
mMPI object of type MPIContext

Definition at line 42 of file FirstDerivativeSecondOrderNoGhosts.c.

53 {
54  MPIContext* mpi = (MPIContext*) m;
55  int i, j, v;
56 
57  if ((!Df) || (!f)) {
58  fprintf(stderr, "Error in FirstDerivativeSecondOrder(): input arrays not allocated.\n");
59  return(1);
60  }
61  if (ghosts < _MINIMUM_GHOSTS_) {
62  fprintf(stderr, "Error in FirstDerivativeSecondOrderCentralNoGhosts(): insufficient number of ghosts.\n");
63  return(1);
64  }
65 
66  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
67  dimension "dir" */
68  int indexC[ndims], index_outer[ndims], bounds_outer[ndims];
69  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
70  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
71 
72 #pragma omp parallel for schedule(auto) default(shared) private(i,j,v,index_outer,indexC)
73  for (j=0; j<N_outer; j++) {
74  _ArrayIndexnD_(ndims,j,bounds_outer,index_outer,0);
75  _ArrayCopy1D_(index_outer,indexC,ndims);
76  for (i = 0; i < dim[dir]; i++) {
77  int qC, qL, qR;
78  indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL);
79  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
80  indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR);
81  for (v=0; v<nvars; v++) Df[qC*nvars+v] = 0.5 * (f[qR*nvars+v]-f[qL*nvars+v]);
82  }
83  }
84 
85  if (mpi->ip[dir] == 0) {
86  /* left physical boundary: overwrite the leftmost value with biased finite-difference */
87 #pragma omp parallel for schedule(auto) default(shared) private(i,j,v,index_outer,indexC)
88  for (j=0; j<N_outer; j++) {
89  _ArrayIndexnD_(ndims,j,bounds_outer,index_outer,0);
90  _ArrayCopy1D_(index_outer,indexC,ndims);
91  i = 0;
92  int qC, qR, qR2;
93  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
94  indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR);
95  indexC[dir] = i+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR2);
96  for (v=0; v<nvars; v++) Df[qC*nvars+v] = (-0.5*f[qR2*nvars+v]+2*f[qR*nvars+v]-1.5*f[qC*nvars+v]);
97  }
98  }
99 
100  if (mpi->ip[dir] == (mpi->iproc[dir]-1)) {
101  /* right physical boundary: overwrite the rightmost value with biased finite-difference */
102 #pragma omp parallel for schedule(auto) default(shared) private(i,j,v,index_outer,indexC)
103  for (j=0; j<N_outer; j++) {
104  _ArrayIndexnD_(ndims,j,bounds_outer,index_outer,0);
105  _ArrayCopy1D_(index_outer,indexC,ndims);
106  i = dim[dir] - 1;
107  int qC, qL, qL2;
108  indexC[dir] = i-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL2);
109  indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL);
110  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
111  for (v=0; v<nvars; v++) Df[qC*nvars+v] = (0.5*f[qL2*nvars+v]-2*f[qL*nvars+v]+1.5*f[qC*nvars+v]);
112  }
113  }
114 
115  return 0;
116 }
_ArrayIndexnD_
#define _ArrayIndexnD_(N, index, imax, i, ghost)
Definition: arrayfunctions.h:21
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
_ArrayProduct1D_
#define _ArrayProduct1D_(x, size, p)
Definition: arrayfunctions.h:377
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
int gpuFirstDerivativeFourthOrderCentral ( double *  Df,
double *  f,
int  dir,
int  bias,
void *  s,
void *  m 
)

Fourth order approximation to the first derivative (note: not divided by grid spacing)

Computes the fourth-order finite-difference approximation to the first derivative (Note: not divided by the grid spacing):

\begin{equation} \left(\partial f\right)_i = \left\{ \begin{array}{ll} \frac{1}{12}\left(-25f_i+48f_{i+1}-36f_{i+2}+16f_{i+3}-3f_{i+4}\right) & i=-g \\ \frac{1}{12}\left(-3f_{i-1}-10f_i+18f_{i+1}-6f_{i+2}+f_{i+3}\right) & i = -g+1 \\ \frac{1}{2}\left( f_{i-2}-8f_{i-1}+8f_{i+1}-f_{i+2} \right) & -g+2 \leq i \leq N+g-3 \\ \frac{1}{12}\left( -f_{i-3}+6f_{i-2}-18f_{i-1}+10f_i+3f_{i+1}\right) & i = N+g-2 \\ \frac{1}{12}\left( 3f_{i-4}-16f_{i-3}+36f_{i-2}-48f_{i-1}+25f_i \right) & i = N+g-1 \end{array}\right. \end{equation}

where \(i\) is the grid index along the spatial dimension of the derivative, \(g\) is the number of ghost points, and \(N\) is the number of grid points (not including the ghost points) in the spatial dimension of the derivative.

Notes:

  • The first derivative is computed at the grid points or the cell centers.
  • The first derivative is computed at the ghost points too. Thus, biased schemes are used at and near the boundaries.
  • Df and f are 1D arrays containing the function and its computed derivatives on a multi- dimensional grid. The derivative along the specified dimension dir is computed by looping through all grid lines along dir.
See Also
FirstDerivativeFourthOrderCentral(), gpuFirstDerivativeFourthOrderCentral()
Parameters
DfArray to hold the computed first derivative (with ghost points)
fArray containing the grid point function values whose first derivative is to be computed (with ghost points)
dirThe spatial dimension along which the derivative is computed
biasForward or backward differencing for non-central finite-difference schemes (-1: backward, 1: forward)
sSolver object of type SolverContext
mMPI object of type MPIContext

Definition at line 351 of file FirstDerivativeFourthOrder_GPU.cu.

361 {
362  SolverContext *solver = (SolverContext*) s;
363 
364  int ghosts = solver->ghosts;
365  int ndims = solver->ndims;
366  int nvars = solver->nvars;
367  int *dim = solver->dim_local;
368 
369  if ((!Df) || (!f)) {
370  fprintf(stderr, "Error in FirstDerivativeFourthOrder(): input arrays not allocated.\n");
371  return(1);
372  }
373 
374  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
375  dimension "dir" */
376  int bounds_outer[ndims];
377  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
378  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
379  int nblocks = (N_outer - 1) / GPU_THREADS_PER_BLOCK + 1;
380 
381 #if defined(GPU_STAT)
382  cudaEvent_t startEvent, stopEvent;
383  float milliseconds;
384 
385  checkCuda( cudaEventCreate(&startEvent) );
386  checkCuda( cudaEventCreate(&stopEvent) );
387 
388  checkCuda( cudaEventRecord(startEvent, 0) );
389 #endif
390 
391  FirstDerivativeFourthOrderCentral_boundary_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
392  N_outer, solver->npoints_local_wghosts, ghosts, ndims, nvars, dir, solver->gpu_dim_local, f, Df
393  );
394 
395 #if defined(GPU_STAT)
396  checkCuda( cudaEventRecord(stopEvent, 0) );
397  checkCuda( cudaEventSynchronize(stopEvent) );
398  checkCuda( cudaEventElapsedTime(&milliseconds, startEvent, stopEvent) );
399 
400  printf("%-50s GPU time (secs) = %.6f\n",
401  "FirstDerivativeFourthOrderCentral_boundary", milliseconds*1e-3);
402 #endif
403 
404  int npoints_grid = N_outer*(dim[dir] + 2*ghosts - 4);
405  nblocks = (npoints_grid - 1) / GPU_THREADS_PER_BLOCK + 1;
406 
407 #if defined(GPU_STAT)
408  checkCuda( cudaEventRecord(startEvent, 0) );
409 #endif
410 
411  FirstDerivativeFourthOrderCentral_interior_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
412  npoints_grid, solver->npoints_local_wghosts, ghosts, ndims, nvars, dir, solver->gpu_dim_local, f, Df
413  );
414  cudaDeviceSynchronize();
415 
416 #if defined(GPU_STAT)
417  checkCuda( cudaEventRecord(stopEvent, 0) );
418  checkCuda( cudaEventSynchronize(stopEvent) );
419  checkCuda( cudaEventElapsedTime(&milliseconds, startEvent, stopEvent) );
420 
421  printf("%-50s GPU time (secs) = %.6f\n",
422  "FirstDerivativeFourthOrderCentral_interior", milliseconds*1e-3);
423 #endif
424 
425  return(0);
426 }
HyPar::npoints_local_wghosts
int npoints_local_wghosts
Definition: hypar.h:42
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
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
HyPar::nvars
int nvars
Definition: hypar.h:29
HyPar::ndims
int ndims
Definition: hypar.h:26
_ArrayProduct1D_
#define _ArrayProduct1D_(x, size, p)
Definition: arrayfunctions.h:377
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23