secondderivative.h File Reference

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

Go to the source code of this file.

Macros
#define	_SECOND_ORDER_CENTRAL_   "2"
#define	_FOURTH_ORDER_CENTRAL_   "4"

Functions
int	SecondDerivativeSecondOrderCentral (double *, double *, int, void *, void *)
int	SecondDerivativeFourthOrderCentral (double *, double *, int, void *, void *)
int	SecondDerivativeSecondOrderCentralNoGhosts (double *, double *, int, int, int *, int, int, void *)

Detailed Description

Definitions for the functions computing the second derivative.

Author

Debojyoti Ghosh

Definition in file secondderivative.h.

Macro Definition Documentation

_SECOND_ORDER_CENTRAL_

#define _SECOND_ORDER_CENTRAL_   "2"

Second order central scheme

Definition at line 10 of file secondderivative.h.

_FOURTH_ORDER_CENTRAL_

#define _FOURTH_ORDER_CENTRAL_   "4"

Fourth order central scheme

Definition at line 12 of file secondderivative.h.

Function Documentation

SecondDerivativeSecondOrderCentral()

int SecondDerivativeSecondOrderCentral	(	double *	D2f,
		double *	f,
		int	dir,
		void *	s,
		void *	m
	)

Second order approximation to the second derivative (note: not divided by square of grid spacing).

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

\begin{equation} \left(\partial^2 f\right)_i = f_{i-1} - 2f_i + f_{i+1} \end{equation}

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

Notes:

• The second derivative is computed at the grid points or the cell centers.
• Though the array D2f includes ghost points, the second derivative is not computed at these locations. Thus, array elements corresponding to the ghost points contain undefined values.
• D2f 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

D2f	Array to hold the computed second derivative (with ghost points) (same size and layout as f)
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
s	Solver object of type HyPar
m	MPI object of type MPIVariables

Definition at line 34 of file SecondDerivativeSecondOrder.c.

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

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


  if ((!D2f) || (!f)) {
    fprintf(stderr, "Error in SecondDerivativeSecondOrder(): input arrays not allocated.\n");
    return(1);
  }
  if (ghosts < _MINIMUM_GHOSTS_) {
    fprintf(stderr, "Error in SecondDerivativeSecondOrder(): 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 done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,indexC,ndims);
    for (i = 0; i < dim[dir]; i++) {
      int qL, qC, 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++)  D2f[qC*nvars+v] = f[qL*nvars+v]-2*f[qC*nvars+v]+f[qR*nvars+v];
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }
  
  return(0);
}

SecondDerivativeFourthOrderCentral()

int SecondDerivativeFourthOrderCentral	(	double *	D2f,
		double *	f,
		int	dir,
		void *	s,
		void *	m
	)

Fourth order approximation to the second derivative (note: not divided by square of grid spacing).

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

\begin{equation} \left(\partial^2 f\right)_i = -\frac{1}{12}f_{i-2} + \frac{4}{3}f_{i-1} - \frac{15}{6}f_i + \frac{4}{3}f_{i+1} - \frac{1}{12}f_{i+2} \end{equation}

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

Notes:

• The second derivative is computed at the grid points or the cell centers.
• Though the array D2f includes ghost points, the second derivative is not computed at these locations. Thus, array elements corresponding to the ghost points contain undefined values.
• D2f 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

D2f	Array to hold the computed second derivative (with ghost points) (same size and layout as f)
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
s	Solver object of type HyPar
m	MPI object of type MPIVariables

Definition at line 34 of file SecondDerivativeFourthOrder.c.

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

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

  static double one_twelve = 1.0/12.0;

  if ((!D2f) || (!f)) {
    fprintf(stderr, "Error in SecondDerivativeFourthOrder(): input arrays not allocated.\n");
    return(1);
  }
  if (ghosts < _MINIMUM_GHOSTS_) {
    fprintf(stderr, "Error in SecondDerivativeFourthOrder(): 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 done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,indexC,ndims);
    for (i = 0; i < dim[dir]; i++) {
      int qm2, qm1, qC, 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++)  
        D2f[qC*nvars+v] = (-f[qm2*nvars+v]+16*f[qm1*nvars+v]-30*f[qC*nvars+v]+16*f[qp1*nvars+v]-f[qp2*nvars+v])*one_twelve;
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }
  
  return(0);
}

SecondDerivativeSecondOrderCentralNoGhosts()

int SecondDerivativeSecondOrderCentralNoGhosts	(	double *	D2f,
		double *	f,
		int	dir,
		int	ndims,
		int *	dim,
		int	ghosts,
		int	nvars,
		void *	m
	)

Second order approximation to the second derivative (note: not divided by square of grid spacing).

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

\begin{equation} \left(\partial^2 f\right)_i = f_{i-1} - 2f_i + f_{i+1} \end{equation}

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

Notes:

• The second derivative is computed at the grid points or the cell centers.
• Though the array D2f must include the same number of ghost points as the array f (including 0), the second derivative is not computed at these locations. Thus, array elements corresponding to the ghost points contain undefined values.
• D2f 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.
• 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

D2f	Array to hold the computed second derivative (with ghost points) (same size and layout as f)
f	Array containing the grid point function values whose second derivative is to be computed (with ghost points)
dir	The spatial dimension along which the derivative is computed
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 39 of file SecondDerivativeSecondOrderNoGhosts.c.

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

  if ((!D2f) || (!f)) {
    fprintf(stderr, "Error in SecondDerivativeSecondOrder(): input arrays not allocated.\n");
    return(1);
  }
  if (ghosts < _MINIMUM_GHOSTS_) {
    fprintf(stderr, "Error in SecondDerivativeSecondOrderNoGhosts(): 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 done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,indexC,ndims);
    for (i = 0; i < dim[dir]; i++) {
      int qL, qC, 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++)  D2f[qC*nvars+v] = f[qL*nvars+v]-2*f[qC*nvars+v]+f[qR*nvars+v];
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  if (mpi->ip[dir] == 0) {
    /* left physical boundary: overwrite the leftmost value with biased finite-difference */
    int done = 0; _ArraySetValue_(index_outer,ndims,0);
    while (!done) {
      _ArrayCopy1D_(index_outer,indexC,ndims);
      i = 0;
      int qC, qR, qR2, qR3;
      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);
      indexC[dir] = i+3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR3);
      for (v=0; v<nvars; v++)  D2f[qC*nvars+v] = 2*f[qC*nvars+v]-5*f[qR*nvars+v]+4*f[qR2*nvars+v]-f[qR3*nvars+v];
      _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
    }
  }

  if (mpi->ip[dir] == (mpi->iproc[dir]-1)) {
    /* right physical boundary: overwrite the rightmost value with biased finite-difference */
    int done = 0; _ArraySetValue_(index_outer,ndims,0);
    while (!done) {
      _ArrayCopy1D_(index_outer,indexC,ndims);
      i = dim[dir]-1;
      int qL3, qL2, qL, qC;
      indexC[dir] = i-3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL3);
      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++)  D2f[qC*nvars+v] = 2*f[qC*nvars+v]-5*f[qL*nvars+v]+4*f[qL2*nvars+v]-f[qL3*nvars+v];
      _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
    }
  }

  
  return 0;
}