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

_FIRST_ORDER_

#define _FIRST_ORDER_   "1"

First order scheme

Definition at line 10 of file firstderivative.h.

_SECOND_ORDER_CENTRAL_

#define _SECOND_ORDER_CENTRAL_   "2"

Second order central scheme

Definition at line 12 of file firstderivative.h.

_FOURTH_ORDER_CENTRAL_

#define _FOURTH_ORDER_CENTRAL_   "4"

Fourth order central scheme

Definition at line 14 of file firstderivative.h.

Function Documentation

FirstDerivativeFirstOrder()

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

Df	Array to hold the computed first derivative (with ghost points)
f	Array containing the grid point function values whose first derivative is to be computed (with ghost points)
dir	The spatial dimension along which the derivative is computed
bias	Forward or backward differencing for non-central finite-difference schemes (-1: backward, 1: forward)
s	Solver object of type SolverContext
m	MPI object of type MPIContext

Definition at line 37 of file FirstDerivativeFirstOrder.c.

{
  SolverContext *solver = (SolverContext*) s;
  int           i, j, v;

  int ghosts = solver->ghosts;
  int ndims  = solver->ndims;
  int nvars  = solver->nvars;
  int *dim   = solver->dim_local;


  if ((!Df) || (!f)) {
    fprintf(stderr, "Error in FirstDerivativeSecondOrder(): input arrays not allocated.\n");
    return(1);
  }

  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
     dimension "dir"                                                                    */
  int indexC[ndims], index_outer[ndims], bounds_outer[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);

#pragma omp parallel for schedule(auto) default(shared) private(i,j,v,index_outer,indexC)
  for (j=0; j<N_outer; j++) {
    _ArrayIndexnD_(ndims,j,bounds_outer,index_outer,0);
    _ArrayCopy1D_(index_outer,indexC,ndims);
    /* left boundary */
    for (i = -ghosts; i < -ghosts+1; i++) {
      int qC, qR;
      indexC[dir] = i  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
      indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR );
      for (v=0; v<nvars; v++)  Df[qC*nvars+v] = f[qR*nvars+v]-f[qC*nvars+v];
    }
    /* interior */
    for (i = -ghosts+1; i < dim[dir]+ghosts-1; i++) {
      int qC, qL, qR;
      indexC[dir] = i  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
      indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL);
      indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR);
      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];
    }
    /* right boundary */
    for (i = dim[dir]+ghosts-1; i < dim[dir]+ghosts; i++) {
      int qL, qC;
      indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL );
      indexC[dir] = i  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
      for (v=0; v<nvars; v++)  Df[qC*nvars+v] = f[qC*nvars+v]-f[qL*nvars+v];
    }
  }
  
  return(0);
}

FirstDerivativeSecondOrderCentral()

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

Df	Array to hold the computed first derivative (with ghost points)
f	Array containing the grid point function values whose first derivative is to be computed (with ghost points)
dir	The spatial dimension along which the derivative is computed
bias	Forward or backward differencing for non-central finite-difference schemes (-1: backward, 1: forward)
s	Solver object of type SolverContext
m	MPI object of type MPIContext

Definition at line 36 of file FirstDerivativeSecondOrder.c.

{
  SolverContext *solver = (SolverContext*) s;
  int           i, j, v;

  int ghosts = solver->ghosts;
  int ndims  = solver->ndims;
  int nvars  = solver->nvars;
  int *dim   = solver->dim_local;


  if ((!Df) || (!f)) {
    fprintf(stderr, "Error in FirstDerivativeSecondOrder(): input arrays not allocated.\n");
    return(1);
  }

  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
     dimension "dir"                                                                    */
  int indexC[ndims], index_outer[ndims], bounds_outer[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);

#pragma omp parallel for schedule(auto) default(shared) private(i,j,v,index_outer,indexC)
  for (j=0; j<N_outer; j++) {
    _ArrayIndexnD_(ndims,j,bounds_outer,index_outer,0);
    _ArrayCopy1D_(index_outer,indexC,ndims);
    /* left boundary */
    for (i = -ghosts; i < -ghosts+1; i++) {
      int qC, qR, qRR;
      indexC[dir] = i  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
      indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR );
      indexC[dir] = i+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qRR);
      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]);
    }
    /* interior */
    for (i = -ghosts+1; i < dim[dir]+ghosts-1; i++) {
      int qC, qL, qR;
      indexC[dir] = i  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
      indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL);
      indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR);
      for (v=0; v<nvars; v++)  Df[qC*nvars+v] = 0.5 * (f[qR*nvars+v]-f[qL*nvars+v]);
    }
    /* right boundary */
    for (i = dim[dir]+ghosts-1; i < dim[dir]+ghosts; i++) {
      int qLL, qL, qC;
      indexC[dir] = i-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qLL);
      indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL );
      indexC[dir] = i  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
      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]);
    }
  }
  
  return(0);
}

FirstDerivativeFourthOrderCentral()

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

Df	Array to hold the computed first derivative (with ghost points)
f	Array containing the grid point function values whose first derivative is to be computed (with ghost points)
dir	The spatial dimension along which the derivative is computed
bias	Forward or backward differencing for non-central finite-difference schemes (-1: backward, 1: forward)
s	Solver object of type SolverContext
m	MPI object of type MPIContext

Definition at line 36 of file FirstDerivativeFourthOrder.c.

{
  SolverContext *solver = (SolverContext*) s;
  int           i, j, v;

  int ghosts = solver->ghosts;
  int ndims  = solver->ndims;
  int nvars  = solver->nvars;
  int *dim   = solver->dim_local;


  if ((!Df) || (!f)) {
    fprintf(stderr, "Error in FirstDerivativeFourthOrder(): input arrays not allocated.\n");
    return(1);
  }

  static double one_twelve = 1.0/12.0;

  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
     dimension "dir"                                                                    */
  int indexC[ndims], index_outer[ndims], bounds_outer[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);

#pragma omp parallel for schedule(auto) default(shared) private(i,j,v,index_outer,indexC)
  for (j=0; j<N_outer; j++) {
    _ArrayIndexnD_(ndims,j,bounds_outer,index_outer,0);
    _ArrayCopy1D_(index_outer,indexC,ndims);
    /* left boundary */
    for (i = -ghosts; i < -ghosts+1; i++) {
      int     qC, qp1, qp2, qp3, qp4;
      indexC[dir] = i  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
      indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
      indexC[dir] = i+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
      indexC[dir] = i+3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp3);
      indexC[dir] = i+4; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp4);
      for (v=0; v<nvars; v++)
        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;
    }
    for (i = -ghosts+1; i < -ghosts+2; i++) {
      int qC, qm1, qp1, qp2, qp3;
      indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
      indexC[dir] = i  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
      indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
      indexC[dir] = i+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
      indexC[dir] = i+3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp3);
      for (v=0; v<nvars; v++)
        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;
    }
    /* interior */
    for (i = -ghosts+2; i < dim[dir]+ghosts-2; i++) {
      int qC, qm1, qm2, qp1, qp2;
      indexC[dir] = i-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
      indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
      indexC[dir] = i  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
      indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
      indexC[dir] = i+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
      for (v=0; v<nvars; v++)
        Df[qC*nvars+v] = (f[qm2*nvars+v]-8*f[qm1*nvars+v]+8*f[qp1*nvars+v]-f[qp2*nvars+v])*one_twelve;
    }
    /* right boundary */
    for (i = dim[dir]+ghosts-2; i < dim[dir]+ghosts-1; i++) {
      int qC, qm3, qm2, qm1, qp1;
      indexC[dir] = i-3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
      indexC[dir] = i-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
      indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
      indexC[dir] = i  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
      indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
      for (v=0; v<nvars; v++)
        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;
    }
    for (i = dim[dir]+ghosts-1; i < dim[dir]+ghosts; i++) {
      int qC, qm4, qm3, qm2, qm1;
      indexC[dir] = i-4; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm4);
      indexC[dir] = i-3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
      indexC[dir] = i-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
      indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
      indexC[dir] = i  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
      for (v=0; v<nvars; v++)
        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;
    }
  }

  return(0);
}

FirstDerivativeSecondOrderCentralNoGhosts()

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

Df	Array to hold the computed first derivative (with ghost points)
f	Array containing the grid point function values whose first derivative is to be computed (with ghost points)
dir	The spatial dimension along which the derivative is computed
bias	Forward or backward differencing for non-central finite-difference schemes (-1: backward, 1: forward)
ndims	Number of spatial/coordinate dimensions
dim	Local dimensions
ghosts	Number of ghost points
nvars	Number of vector components at each grid points
m	MPI object of type MPIContext

Definition at line 42 of file FirstDerivativeSecondOrderNoGhosts.c.

{
  MPIContext* mpi = (MPIContext*) m;
  int  i, j, v;

  if ((!Df) || (!f)) {
    fprintf(stderr, "Error in FirstDerivativeSecondOrder(): input arrays not allocated.\n");
    return(1);
  }
  if (ghosts < _MINIMUM_GHOSTS_) {
    fprintf(stderr, "Error in FirstDerivativeSecondOrderCentralNoGhosts(): insufficient number of ghosts.\n");
    return(1);
  }

  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
     dimension "dir"                                                                    */
  int indexC[ndims], index_outer[ndims], bounds_outer[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  int N_outer;  _ArrayProduct1D_(bounds_outer,ndims,N_outer);

#pragma omp parallel for schedule(auto) default(shared) private(i,j,v,index_outer,indexC)
  for (j=0; j<N_outer; j++) {
    _ArrayIndexnD_(ndims,j,bounds_outer,index_outer,0);
    _ArrayCopy1D_(index_outer,indexC,ndims);
    for (i = 0; i < dim[dir]; i++) {
      int qC, qL, qR;
      indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL);
      indexC[dir] = i  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
      indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR);
      for (v=0; v<nvars; v++)  Df[qC*nvars+v] = 0.5 * (f[qR*nvars+v]-f[qL*nvars+v]);
    }
  }

  if (mpi->ip[dir] == 0) {
    /* left physical boundary: overwrite the leftmost value with biased finite-difference */
#pragma omp parallel for schedule(auto) default(shared) private(i,j,v,index_outer,indexC)
    for (j=0; j<N_outer; j++) {
      _ArrayIndexnD_(ndims,j,bounds_outer,index_outer,0);
      _ArrayCopy1D_(index_outer,indexC,ndims);
      i = 0;
      int qC, qR, qR2;
      indexC[dir] = i  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
      indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR);
      indexC[dir] = i+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR2);
      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]);
    }
  }

  if (mpi->ip[dir] == (mpi->iproc[dir]-1)) {
    /* right physical boundary: overwrite the rightmost value with biased finite-difference */
#pragma omp parallel for schedule(auto) default(shared) private(i,j,v,index_outer,indexC)
    for (j=0; j<N_outer; j++) {
      _ArrayIndexnD_(ndims,j,bounds_outer,index_outer,0);
      _ArrayCopy1D_(index_outer,indexC,ndims);
      i = dim[dir] - 1;
      int qC, qL, qL2;
      indexC[dir] = i-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL2);
      indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL);
      indexC[dir] = i  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC );
      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]);
    }
  }

  return 0;
}

gpuFirstDerivativeFourthOrderCentral()

int gpuFirstDerivativeFourthOrderCentral	(	double *	,
		double *	,
		int	,
		int	,
		void *	,
		void *
	)

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