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

2nd order discretization of the second derivative. More...

#include <stdio.h>
#include <stdlib.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <secondderivative.h>
#include <mpivars.h>

Go to the source code of this file.

Macros

#define _MINIMUM_GHOSTS_   1
 

Typedefs

typedef MPIVariables MPIContext
 

Functions

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

Detailed Description

2nd order discretization of the second derivative.

Author
Ping-Hsuan Tsai, Debojyoti Ghosh

Definition in file SecondDerivativeSecondOrderNoGhosts.c.

Macro Definition Documentation

#define _MINIMUM_GHOSTS_   1

Minimum number of ghost points required.

Definition at line 19 of file SecondDerivativeSecondOrderNoGhosts.c.

Typedef Documentation

Definition at line 13 of file SecondDerivativeSecondOrderNoGhosts.c.

Function Documentation

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

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

49 {
50  MPIContext* mpi = (MPIContext*) m;
51  int i, v;
52 
53  if ((!D2f) || (!f)) {
54  fprintf(stderr, "Error in SecondDerivativeSecondOrder(): input arrays not allocated.\n");
55  return(1);
56  }
57  if (ghosts < _MINIMUM_GHOSTS_) {
58  fprintf(stderr, "Error in SecondDerivativeSecondOrderNoGhosts(): insufficient number of ghosts.\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 
67  int done = 0; _ArraySetValue_(index_outer,ndims,0);
68  while (!done) {
69  _ArrayCopy1D_(index_outer,indexC,ndims);
70  for (i = 0; i < dim[dir]; i++) {
71  int qL, qC, qR;
72  indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL);
73  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC);
74  indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR);
75  for (v=0; v<nvars; v++) D2f[qC*nvars+v] = f[qL*nvars+v]-2*f[qC*nvars+v]+f[qR*nvars+v];
76  }
77  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
78  }
79 
80  if (mpi->ip[dir] == 0) {
81  /* left physical boundary: overwrite the leftmost value with biased finite-difference */
82  int done = 0; _ArraySetValue_(index_outer,ndims,0);
83  while (!done) {
84  _ArrayCopy1D_(index_outer,indexC,ndims);
85  i = 0;
86  int qC, qR, qR2, qR3;
87  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC);
88  indexC[dir] = i+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR);
89  indexC[dir] = i+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR2);
90  indexC[dir] = i+3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR3);
91  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];
92  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
93  }
94  }
95 
96  if (mpi->ip[dir] == (mpi->iproc[dir]-1)) {
97  /* right physical boundary: overwrite the rightmost value with biased finite-difference */
98  int done = 0; _ArraySetValue_(index_outer,ndims,0);
99  while (!done) {
100  _ArrayCopy1D_(index_outer,indexC,ndims);
101  i = dim[dir]-1;
102  int qL3, qL2, qL, qC;
103  indexC[dir] = i-3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL3);
104  indexC[dir] = i-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL2);
105  indexC[dir] = i-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL);
106  indexC[dir] = i ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qC);
107  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];
108  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
109  }
110  }
111 
112 
113  return 0;
114 }
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayCopy1D_(x, y, size)
Structure of MPI-related variables.