mathfunctions.h File Reference

Contains function definitions for common mathematical functions. More...

#include <math_ops.h>

Go to the source code of this file.

Functions
void	FindInterval (double, double, double *, int, int *, int *)
void	fillGhostCells (const int *const, const int, double *const, const int, const int, const int *const)
void	TrilinearInterpCoeffs (double, double, double, double, double, double, double, double, double, double *)
int	InterpolateGlobalnDVar (const int *const, double **const, const int *const, double *const, const int, const int, const int, const int *const)

Detailed Description

Contains function definitions for common mathematical functions.

Author

Debojyoti Ghosh

Definition in file mathfunctions.h.

Function Documentation

void FindInterval	(	double	a,
		double	b,
		double *	x,
		int	N,
		int *	imin,
		int *	imax
	)

Function to calculate the grid points corresponding to a given interval

Given an interval \(\left[a,b\right], a\leq b\), find grid indices imin and imax, such that

\begin{align} imin &= \min\ i\ {\rm satisfying}\ x_i \geq a\\ imax &= \max\ i\ {\rm satisfying}\ x_i \leq b \end{align}

where \(\left\{x_i; 0\leq i < N , x_i < x_{i+1} \forall i \right\}\) represents a 1-dimensional grid.

Note: This function handles 1-dimensional intervals and grids only.

Parameters

a	Lower bound of interval
b	Upper bound of interval
x	Array of spatial coordinates representing a grid
N	Number of grid points / size of x
imin	Lowest grid index within [a,b]
imax	Highest grid index within [a,b]

Definition at line 19 of file FindInterval.c.

{
  int i;
  *imax = -1;
  *imin =  N;

  double min_dx = x[1] - x[0];
  for (i = 2; i < N; i++) {
    double dx = x[i] - x[i-1];
    if (dx < min_dx) min_dx = dx;
  }
  double tol = 1e-10 * min_dx;

  for (i = 0; i < N; i++) {
    if (x[i] <= (b+tol)) *imax = i+1;
  }
  for (i = N-1; i > -1; i--) {
    if (x[i] >= (a-tol)) *imin = i;
  }

  return;
}

void fillGhostCells	(	const int *const	a_dim,
		const int	a_ngpt,
		double *const	a_u,
		const int	a_nvars,
		const int	a_ndims,
		const int *const	a_periodic
	)

Fill the ghost cells of a global n-dimensional array

Fill the ghost cells of a solution. Note that the solution array must be a global one (not one that is partitioned among MPI ranks). This is not a parallel operation (it will execute independently on multiple MPI ranks, if called by multiple processes).

For periodicity along any dimension, the ghost cells are filled appropriately from the other side of the domain. Otherwise, the interior data is extrapolated by a 4th order polynomial (assuming uniform grid spacing).

Parameters

a_dim	grid dimensions of solution
a_ngpt	number of ghost cells
a_u	solution array
a_nvars	number of vector components of the solution
a_ndims	number of spatial dimensions
a_periodic	periodic or not along each dimension

Definition at line 17 of file FillGhostCells.c.

{
  for (int d = 0; d < a_ndims; d++) {

    int bounds[a_ndims];
    _ArrayCopy1D_(a_dim, bounds, a_ndims);
    bounds[d] = a_ngpt;

    int index[a_ndims];
    _ArraySetValue_(index, a_ndims, 0);

    if (a_periodic[d]) {
      
      /* periodic case */

      int done = 0;
      while (!done) {

        {
          /* low end - face = 1 */

          int p_gpt = 0, 
              p_int = 0;

          int index_gpt[a_ndims];
          _ArrayCopy1D_(index, index_gpt, a_ndims);
          index_gpt[d] -= a_ngpt;
          _ArrayIndex1D_(a_ndims, a_dim, index_gpt, a_ngpt, p_gpt);

          int index_int[a_ndims];
          _ArrayCopy1D_(index, index_int, a_ndims);
          index_int[d] += (a_dim[d]-a_ngpt);
          _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int);

          _ArrayCopy1D_((a_u+a_nvars*p_int), (a_u+a_nvars*p_gpt), a_nvars);
        }

        {
          /* high end - face = -1 */

          int p_gpt = 0, 
              p_int = 0;

          int index_gpt[a_ndims];
          _ArrayCopy1D_(index, index_gpt, a_ndims);
          index_gpt[d] += a_dim[d];
          _ArrayIndex1D_(a_ndims, a_dim, index_gpt, a_ngpt, p_gpt);

          int index_int[a_ndims];
          _ArrayCopy1D_(index, index_int, a_ndims);
          _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int);

          _ArrayCopy1D_((a_u+a_nvars*p_int), (a_u+a_nvars*p_gpt), a_nvars);
        }

        _ArrayIncrementIndex_(a_ndims, bounds, index, done);

      }

    } else {

      /* not periodic - extrapolate */

      int done = 0;
      while (!done) {

        {
          /* low end - face = 1 */

          int p_gpt = 0, 
              p_int_0 = 0,
              p_int_1 = 0,
              p_int_2 = 0,
              p_int_3 = 0;

          int index_gpt[a_ndims];
          _ArrayCopy1D_(index, index_gpt, a_ndims);
          index_gpt[d] -= a_ngpt;
          _ArrayIndex1D_(a_ndims, a_dim, index_gpt, a_ngpt, p_gpt);

          int index_int[a_ndims];
          _ArrayCopy1D_(index, index_int, a_ndims);

          index_int[d] = 0;
          _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int_0);
          index_int[d]++;
          _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int_1);
          index_int[d]++;
          _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int_2);
          index_int[d]++;
          _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int_3);

          double alpha = - (double) (a_ngpt - index[d]);
          double c0 = -((-2.0 + alpha)*(-1.0 + alpha)*alpha)/6.0;
          double c1 = ((-2.0 + alpha)*(-1.0 + alpha)*(1.0 + alpha))/2.0;
          double c2 = (alpha*(2.0 + alpha - alpha*alpha))/2.0;
          double c3 = (alpha*(-1.0 + alpha*alpha))/6.0;

          for (int v = 0; v < a_nvars; v++) {
    
            a_u[p_gpt*a_nvars+v] =    c0 * a_u[p_int_0*a_nvars+v]
                                    + c1 * a_u[p_int_1*a_nvars+v]
                                    + c2 * a_u[p_int_2*a_nvars+v]
                                    + c3 * a_u[p_int_3*a_nvars+v];
    
          }

        }

        {
          /* high end - face = -1 */

          int p_gpt = 0, 
              p_int_0 = 0,
              p_int_1 = 0,
              p_int_2 = 0,
              p_int_3 = 0;

          int index_gpt[a_ndims];
          _ArrayCopy1D_(index, index_gpt, a_ndims);
          index_gpt[d] += a_dim[d];
          _ArrayIndex1D_(a_ndims, a_dim, index_gpt, a_ngpt, p_gpt);

          int index_int[a_ndims];
          _ArrayCopy1D_(index, index_int, a_ndims);

          index_int[d] = a_dim[d]-1;
          _ArrayIndex1D_(a_ndims, a_dim, index, a_ngpt, p_int_0);
          index_int[d]--;
          _ArrayIndex1D_(a_ndims, a_dim, index, a_ngpt, p_int_1);
          index_int[d]--;
          _ArrayIndex1D_(a_ndims, a_dim, index, a_ngpt, p_int_2);
          index_int[d]--;
          _ArrayIndex1D_(a_ndims, a_dim, index, a_ngpt, p_int_3);

          double alpha = - (double) (index[d]+1);
          double c0 = -((-2.0 + alpha)*(-1.0 + alpha)*alpha)/6.0;
          double c1 = ((-2.0 + alpha)*(-1.0 + alpha)*(1.0 + alpha))/2.0;
          double c2 = (alpha*(2.0 + alpha - alpha*alpha))/2.0;
          double c3 = (alpha*(-1.0 + alpha*alpha))/6.0;

          for (int v = 0; v < a_nvars; v++) {
    
            a_u[p_gpt*a_nvars+v] =    c0 * a_u[p_int_0*a_nvars+v]
                                    + c1 * a_u[p_int_1*a_nvars+v]
                                    + c2 * a_u[p_int_2*a_nvars+v]
                                    + c3 * a_u[p_int_3*a_nvars+v];
    
          }

        }

        _ArrayIncrementIndex_(a_ndims, bounds, index, done);

      }

    }

  }

  return;
}

void TrilinearInterpCoeffs	(	double	xmin,
		double	xmax,
		double	ymin,
		double	ymax,
		double	zmin,
		double	zmax,
		double	x,
		double	y,
		double	z,
		double *	coeffs
	)

Function to compute trilinear interpolation coefficients

This function computes the coefficients for a trilinear interpolation at a given point (x,y,z) inside a cube defined by [xmin,xmax] X [ymin,ymax] X [zmin,zmax]. The coefficients are stored in an array of size 8 with each element corresponding to a corner of the cube in the following order:
coeffs[0] => xmin,ymin,zmin
coeffs[1] => xmax,ymin,zmin
coeffs[2] => xmin,ymax,zmin
coeffs[3] => xmax,ymax,zmin
coeffs[4] => xmin,ymin,zmax
coeffs[5] => xmax,ymin,zmax
coeffs[6] => xmin,ymax,zmax
coeffs[7] => xmax,ymax,zmax

Parameters

xmin	x-coordinate of the lower-end
xmax	x-coordinate of the higher-end
ymin	y-coordinate of the lower-end
ymax	y-coordinate of the higher-end
zmin	z-coordinate of the lower-end
zmax	z-coordinate of the higher-end
x	x-coordinate of the point to interpolate at
y	y-coordinate of the point to interpolate at
z	z-coordinate of the point to interpolate at
coeffs	array of size 8 (pre-allocated) to store the coefficients in

Definition at line 22 of file TrilinearInterpolation.c.

{
    double vol_inv = 1 / ((xmax-xmin)*(ymax-ymin)*(zmax-zmin));
    double tldx1 = x - xmin;
    double tldx2 = xmax - x;
    double tldy1 = y - ymin;
    double tldy2 = ymax - y;
    double tldz1 = z - zmin;
    double tldz2 = zmax - z;

    coeffs[0] = tldz2 * tldy2 * tldx2 * vol_inv;
    coeffs[1] = tldz2 * tldy2 * tldx1 * vol_inv;
    coeffs[2] = tldz2 * tldy1 * tldx2 * vol_inv;
    coeffs[3] = tldz2 * tldy1 * tldx1 * vol_inv;
    coeffs[4] = tldz1 * tldy2 * tldx2 * vol_inv;
    coeffs[5] = tldz1 * tldy2 * tldx1 * vol_inv;
    coeffs[6] = tldz1 * tldy1 * tldx2 * vol_inv;
    coeffs[7] = tldz1 * tldy1 * tldx1 * vol_inv;

    return;
}

int InterpolateGlobalnDVar	(	const int *const	a_dim_dst,
		double **const	a_u_dst,
		const int *const	a_dim_src,
		double *const	a_u_src,
		const int	a_nvars,
		const int	a_ghosts,
		const int	a_ndims,
		const int *const	a_periodic
	)

Interpolate an n-dimensional grid variable from one grid to another

Interpolate n-dimensional data from one grid to another of a desired resolution. Note that along each dimension, the ratio of the number of grid points in the source grid and that in the destination grid must be an integer power of 2 (negative or positive).

The source data must be global. It will get deallocated at the end of this function. It must have the appropriate number of ghost points.

The incoming pointer for destination data must be NULL. After this function is executed, it will point to a chunk of memory with the interpolated solution. This is the global solution. It will have the specified number of ghost points appropriately filled.

Parameters

a_dim_dst	grid dimensions to interpolate to
a_u_dst	pointer to array containing interpolated data
a_dim_src	grid dimensions to interpolate from
a_u_src	pointer to array containing data to interpolate from
a_nvars	Number of vector components of the solution
a_ghosts	Number of ghost points
a_ndims	Number of spatial dimensions
a_periodic	Is the domain periodic or not along each dimension

Definition at line 317 of file InterpolateGlobalnDVar.c.

{
  if ((*a_u_dst) != NULL) {
    fprintf(stderr, "Error in InterpolateGlobalnDVar() - \n");
    fprintf(stderr, "  a_u_dst is not NULL!\n");
    return 1;
  }

  if (a_u_src == NULL) {
    fprintf(stderr, "Error in InterpolateGlobalnDVar() - \n");
    fprintf(stderr, "  a_u_src is NULL!\n");
    return 1;
  }

  int dim_to[a_ndims], dim_from[a_ndims];

  double *u_from;
  double *u_to;

  _ArrayCopy1D_(a_dim_src, dim_to, a_ndims);
  u_to = a_u_src;
  u_from = NULL;

  for (int dir = 0; dir < a_ndims; dir++) {

    _ArrayCopy1D_(dim_to, dim_from, a_ndims);
    dim_to[dir] = a_dim_dst[dir];

    if (dim_from[dir] == dim_to[dir]) continue;

    double fac = (dim_to[dir] > dim_from[dir] ? 
                      (double)dim_to[dir]/(double)dim_from[dir] 
                    : (double)dim_from[dir]/(double)dim_to[dir] );

    if (!isPowerOfTwo((int)fac)) {
      fprintf(stderr,"Error in interpolate() - \n");
      fprintf(stderr,"  refinement/coarsening factor not a power of 2!\n");
      return 1;
    }

    if (u_from != NULL) free(u_from);

    u_from = u_to;

    {
      long size = (long) a_nvars;
      for (int d = 0; d < a_ndims; d++) {
        size *= (long) (dim_to[d] + 2*a_ghosts);
      }
      u_to = (double*) calloc (size, sizeof(double));
    }

    fillGhostCells( dim_from, 
                    a_ghosts, 
                    u_from, 
                    a_nvars, 
                    a_ndims, 
                    a_periodic );

    if (dim_to[dir] < dim_from[dir]) {
      int retval = coarsen1D( dim_from, 
                              dim_to, 
                              u_from, 
                              u_to, 
                              dir, 
                              a_nvars, 
                              a_ghosts, 
                              a_ndims );
      if (retval) return retval;
    } else {
      int retval = refine1D(  dim_from, 
                              dim_to, 
                              u_from, 
                              u_to, 
                              dir, 
                              a_nvars, 
                              a_ghosts, 
                              a_ndims );
      if (retval) return retval;
    }

  }

  /* dim_to should be equal to a_dim_dst now */
  for (int d = 0; d < a_ndims; d++) {
    if (dim_to[d] != a_dim_dst[d]) {
      fprintf(stderr,"Error in InterpolateGlobalnDVar() - \n");
      fprintf(stderr,"  dim_to[%d] (%d) != a_dim_dst[%d] (%d)!\n", 
              d, dim_to[d], d, a_dim_dst[d]);
      return 1;
    }
  }

  if (u_from != NULL) free(u_from);
  (*a_u_dst) = u_to;

  return 0;
}