HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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
aLower bound of interval
bUpper bound of interval
xArray of spatial coordinates representing a grid
NNumber of grid points / size of x
iminLowest grid index within [a,b]
imaxHighest grid index within [a,b]

Definition at line 19 of file FindInterval.c.

27 {
28  int i;
29  *imax = -1;
30  *imin = N;
31 
32  double min_dx = x[1] - x[0];
33  for (i = 2; i < N; i++) {
34  double dx = x[i] - x[i-1];
35  if (dx < min_dx) min_dx = dx;
36  }
37  double tol = 1e-10 * min_dx;
38 
39  for (i = 0; i < N; i++) {
40  if (x[i] <= (b+tol)) *imax = i+1;
41  }
42  for (i = N-1; i > -1; i--) {
43  if (x[i] >= (a-tol)) *imin = i;
44  }
45 
46  return;
47 }
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_dimgrid dimensions of solution
a_ngptnumber of ghost cells
a_usolution array
a_nvarsnumber of vector components of the solution
a_ndimsnumber of spatial dimensions
a_periodicperiodic or not along each dimension

Definition at line 17 of file FillGhostCells.c.

24 {
25  for (int d = 0; d < a_ndims; d++) {
26 
27  int bounds[a_ndims];
28  _ArrayCopy1D_(a_dim, bounds, a_ndims);
29  bounds[d] = a_ngpt;
30 
31  int index[a_ndims];
32  _ArraySetValue_(index, a_ndims, 0);
33 
34  if (a_periodic[d]) {
35 
36  /* periodic case */
37 
38  int done = 0;
39  while (!done) {
40 
41  {
42  /* low end - face = 1 */
43 
44  int p_gpt = 0,
45  p_int = 0;
46 
47  int index_gpt[a_ndims];
48  _ArrayCopy1D_(index, index_gpt, a_ndims);
49  index_gpt[d] -= a_ngpt;
50  _ArrayIndex1D_(a_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
51 
52  int index_int[a_ndims];
53  _ArrayCopy1D_(index, index_int, a_ndims);
54  index_int[d] += (a_dim[d]-a_ngpt);
55  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int);
56 
57  _ArrayCopy1D_((a_u+a_nvars*p_int), (a_u+a_nvars*p_gpt), a_nvars);
58  }
59 
60  {
61  /* high end - face = -1 */
62 
63  int p_gpt = 0,
64  p_int = 0;
65 
66  int index_gpt[a_ndims];
67  _ArrayCopy1D_(index, index_gpt, a_ndims);
68  index_gpt[d] += a_dim[d];
69  _ArrayIndex1D_(a_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
70 
71  int index_int[a_ndims];
72  _ArrayCopy1D_(index, index_int, a_ndims);
73  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int);
74 
75  _ArrayCopy1D_((a_u+a_nvars*p_int), (a_u+a_nvars*p_gpt), a_nvars);
76  }
77 
78  _ArrayIncrementIndex_(a_ndims, bounds, index, done);
79 
80  }
81 
82  } else {
83 
84  /* not periodic - extrapolate */
85 
86  int done = 0;
87  while (!done) {
88 
89  {
90  /* low end - face = 1 */
91 
92  int p_gpt = 0,
93  p_int_0 = 0,
94  p_int_1 = 0,
95  p_int_2 = 0,
96  p_int_3 = 0;
97 
98  int index_gpt[a_ndims];
99  _ArrayCopy1D_(index, index_gpt, a_ndims);
100  index_gpt[d] -= a_ngpt;
101  _ArrayIndex1D_(a_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
102 
103  int index_int[a_ndims];
104  _ArrayCopy1D_(index, index_int, a_ndims);
105 
106  index_int[d] = 0;
107  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int_0);
108  index_int[d]++;
109  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int_1);
110  index_int[d]++;
111  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int_2);
112  index_int[d]++;
113  _ArrayIndex1D_(a_ndims, a_dim, index_int, a_ngpt, p_int_3);
114 
115  double alpha = - (double) (a_ngpt - index[d]);
116  double c0 = -((-2.0 + alpha)*(-1.0 + alpha)*alpha)/6.0;
117  double c1 = ((-2.0 + alpha)*(-1.0 + alpha)*(1.0 + alpha))/2.0;
118  double c2 = (alpha*(2.0 + alpha - alpha*alpha))/2.0;
119  double c3 = (alpha*(-1.0 + alpha*alpha))/6.0;
120 
121  for (int v = 0; v < a_nvars; v++) {
122 
123  a_u[p_gpt*a_nvars+v] = c0 * a_u[p_int_0*a_nvars+v]
124  + c1 * a_u[p_int_1*a_nvars+v]
125  + c2 * a_u[p_int_2*a_nvars+v]
126  + c3 * a_u[p_int_3*a_nvars+v];
127 
128  }
129 
130  }
131 
132  {
133  /* high end - face = -1 */
134 
135  int p_gpt = 0,
136  p_int_0 = 0,
137  p_int_1 = 0,
138  p_int_2 = 0,
139  p_int_3 = 0;
140 
141  int index_gpt[a_ndims];
142  _ArrayCopy1D_(index, index_gpt, a_ndims);
143  index_gpt[d] += a_dim[d];
144  _ArrayIndex1D_(a_ndims, a_dim, index_gpt, a_ngpt, p_gpt);
145 
146  int index_int[a_ndims];
147  _ArrayCopy1D_(index, index_int, a_ndims);
148 
149  index_int[d] = a_dim[d]-1;
150  _ArrayIndex1D_(a_ndims, a_dim, index, a_ngpt, p_int_0);
151  index_int[d]--;
152  _ArrayIndex1D_(a_ndims, a_dim, index, a_ngpt, p_int_1);
153  index_int[d]--;
154  _ArrayIndex1D_(a_ndims, a_dim, index, a_ngpt, p_int_2);
155  index_int[d]--;
156  _ArrayIndex1D_(a_ndims, a_dim, index, a_ngpt, p_int_3);
157 
158  double alpha = - (double) (index[d]+1);
159  double c0 = -((-2.0 + alpha)*(-1.0 + alpha)*alpha)/6.0;
160  double c1 = ((-2.0 + alpha)*(-1.0 + alpha)*(1.0 + alpha))/2.0;
161  double c2 = (alpha*(2.0 + alpha - alpha*alpha))/2.0;
162  double c3 = (alpha*(-1.0 + alpha*alpha))/6.0;
163 
164  for (int v = 0; v < a_nvars; v++) {
165 
166  a_u[p_gpt*a_nvars+v] = c0 * a_u[p_int_0*a_nvars+v]
167  + c1 * a_u[p_int_1*a_nvars+v]
168  + c2 * a_u[p_int_2*a_nvars+v]
169  + c3 * a_u[p_int_3*a_nvars+v];
170 
171  }
172 
173  }
174 
175  _ArrayIncrementIndex_(a_ndims, bounds, index, done);
176 
177  }
178 
179  }
180 
181  }
182 
183  return;
184 }
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayCopy1D_(x, y, size)
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
xminx-coordinate of the lower-end
xmaxx-coordinate of the higher-end
yminy-coordinate of the lower-end
ymaxy-coordinate of the higher-end
zminz-coordinate of the lower-end
zmaxz-coordinate of the higher-end
xx-coordinate of the point to interpolate at
yy-coordinate of the point to interpolate at
zz-coordinate of the point to interpolate at
coeffsarray of size 8 (pre-allocated) to store the coefficients in

Definition at line 22 of file TrilinearInterpolation.c.

34 {
35  double vol_inv = 1 / ((xmax-xmin)*(ymax-ymin)*(zmax-zmin));
36  double tldx1 = x - xmin;
37  double tldx2 = xmax - x;
38  double tldy1 = y - ymin;
39  double tldy2 = ymax - y;
40  double tldz1 = z - zmin;
41  double tldz2 = zmax - z;
42 
43  coeffs[0] = tldz2 * tldy2 * tldx2 * vol_inv;
44  coeffs[1] = tldz2 * tldy2 * tldx1 * vol_inv;
45  coeffs[2] = tldz2 * tldy1 * tldx2 * vol_inv;
46  coeffs[3] = tldz2 * tldy1 * tldx1 * vol_inv;
47  coeffs[4] = tldz1 * tldy2 * tldx2 * vol_inv;
48  coeffs[5] = tldz1 * tldy2 * tldx1 * vol_inv;
49  coeffs[6] = tldz1 * tldy1 * tldx2 * vol_inv;
50  coeffs[7] = tldz1 * tldy1 * tldx1 * vol_inv;
51 
52  return;
53 }
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_dstgrid dimensions to interpolate to
a_u_dstpointer to array containing interpolated data
a_dim_srcgrid dimensions to interpolate from
a_u_srcpointer to array containing data to interpolate from
a_nvarsNumber of vector components of the solution
a_ghostsNumber of ghost points
a_ndimsNumber of spatial dimensions
a_periodicIs the domain periodic or not along each dimension

Definition at line 317 of file InterpolateGlobalnDVar.c.

326 {
327  if ((*a_u_dst) != NULL) {
328  fprintf(stderr, "Error in InterpolateGlobalnDVar() - \n");
329  fprintf(stderr, " a_u_dst is not NULL!\n");
330  return 1;
331  }
332 
333  if (a_u_src == NULL) {
334  fprintf(stderr, "Error in InterpolateGlobalnDVar() - \n");
335  fprintf(stderr, " a_u_src is NULL!\n");
336  return 1;
337  }
338 
339  int dim_to[a_ndims], dim_from[a_ndims];
340 
341  double *u_from;
342  double *u_to;
343 
344  _ArrayCopy1D_(a_dim_src, dim_to, a_ndims);
345  u_to = a_u_src;
346  u_from = NULL;
347 
348  for (int dir = 0; dir < a_ndims; dir++) {
349 
350  _ArrayCopy1D_(dim_to, dim_from, a_ndims);
351  dim_to[dir] = a_dim_dst[dir];
352 
353  if (dim_from[dir] == dim_to[dir]) continue;
354 
355  double fac = (dim_to[dir] > dim_from[dir] ?
356  (double)dim_to[dir]/(double)dim_from[dir]
357  : (double)dim_from[dir]/(double)dim_to[dir] );
358 
359  if (!isPowerOfTwo((int)fac)) {
360  fprintf(stderr,"Error in interpolate() - \n");
361  fprintf(stderr," refinement/coarsening factor not a power of 2!\n");
362  return 1;
363  }
364 
365  if (u_from != NULL) free(u_from);
366 
367  u_from = u_to;
368 
369  {
370  long size = (long) a_nvars;
371  for (int d = 0; d < a_ndims; d++) {
372  size *= (long) (dim_to[d] + 2*a_ghosts);
373  }
374  u_to = (double*) calloc (size, sizeof(double));
375  }
376 
377  fillGhostCells( dim_from,
378  a_ghosts,
379  u_from,
380  a_nvars,
381  a_ndims,
382  a_periodic );
383 
384  if (dim_to[dir] < dim_from[dir]) {
385  int retval = coarsen1D( dim_from,
386  dim_to,
387  u_from,
388  u_to,
389  dir,
390  a_nvars,
391  a_ghosts,
392  a_ndims );
393  if (retval) return retval;
394  } else {
395  int retval = refine1D( dim_from,
396  dim_to,
397  u_from,
398  u_to,
399  dir,
400  a_nvars,
401  a_ghosts,
402  a_ndims );
403  if (retval) return retval;
404  }
405 
406  }
407 
408  /* dim_to should be equal to a_dim_dst now */
409  for (int d = 0; d < a_ndims; d++) {
410  if (dim_to[d] != a_dim_dst[d]) {
411  fprintf(stderr,"Error in InterpolateGlobalnDVar() - \n");
412  fprintf(stderr," dim_to[%d] (%d) != a_dim_dst[%d] (%d)!\n",
413  d, dim_to[d], d, a_dim_dst[d]);
414  return 1;
415  }
416  }
417 
418  if (u_from != NULL) free(u_from);
419  (*a_u_dst) = u_to;
420 
421  return 0;
422 }
static int refine1D(const int *const a_dim_src, const int *const a_dim_dst, const double *const a_u_src, double *const a_u_dst, const int a_dir, const int a_nvars, const int a_ngpt, const int a_ndims)
void fillGhostCells(const int *const, const int, double *const, const int, const int, const int *const)
static int coarsen1D(const int *const a_dim_src, const int *const a_dim_dst, const double *const a_u_src, double *const a_u_dst, const int a_dir, const int a_nvars, const int a_ngpt, const int a_ndims)
#define _ArrayCopy1D_(x, y, size)
static int isPowerOfTwo(int x)