HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
Functions
InterpolateGlobalnDVar.c File Reference

Functions to interpolate from one grid to another. More...

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <arrayfunctions.h>
#include <mathfunctions.h>

Go to the source code of this file.

Functions

static int isPowerOfTwo (int x)
 
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)
 
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)
 
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)
 

Detailed Description

Functions to interpolate from one grid to another.

Author
Debojyoti Ghosh

Definition in file InterpolateGlobalnDVar.c.

Function Documentation

◆ isPowerOfTwo()

static int isPowerOfTwo ( int  x)
static

Is the input number an integer power of 2?

Definition at line 13 of file InterpolateGlobalnDVar.c.

14 {
15  if (x == 0) return 0;
16 
17  while (x > 1) {
18  if (x%2 != 0) return 0;
19  x /= 2;
20  }
21  return 1;
22 }

◆ coarsen1D()

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 
)
static

Coarsen along a given dimension - the source and destination must have the same sizes along all other dimensions.

Note that the arrays must not have any ghost points!

Currently this function can only handle coarsening factors that are integer powers of 2.

Parameters
a_dim_srcGrid size of source data
a_dim_dstGrid size of destination data
a_u_srcSource solution
a_u_dstDestination solution
a_dirDimension along which to coarsen
a_nvarsNumber of vector components of the solution
a_ngptNumber of ghost points
a_ndimsNumber of spatial dimensions

Definition at line 32 of file InterpolateGlobalnDVar.c.

41 {
42  for (int d = 0; d < a_ndims; d++) {
43  if ((d != a_dir) && (a_dim_src[d] != a_dim_dst[d])) {
44  fprintf(stderr, "Error in coarsen1D() -\n");
45  fprintf(stderr, " a_dim_src[%d] != a_dim_dst[%d]\n", d, d);
46  return 1;
47  }
48  }
49 
50  int n_src = a_dim_src[a_dir];
51  int n_dst = a_dim_dst[a_dir];
52  if (n_dst > n_src) {
53  fprintf(stderr, "Error in coarsen1D() -\n");
54  fprintf(stderr, " destination grid is finer than source grid along a_dir!\n");
55  return 1;
56  }
57 
58  double fac = ((double) n_src) / ((double) n_dst);
59  int stride = (int) fac;
60  if (absolute(((double)stride)-fac) > _MACHINE_ZERO_) {
61  fprintf(stderr, "Error in coarsen1D() -\n");
62  fprintf(stderr, " non-integer coarsening factor!\n");
63  return 1;
64  }
65 
66  /* set interpolation coefficients depending on desired order */
67  double c0, c1, c2, c3, c4, c5;
68 // if (a_interp_order == 2) {
69 // c0 = c5 = 0.0;
70 // c1 = c4 = 0.0;
71 // c2 = c3 = 0.5;
72 // } else if (a_interp_order == 4) {
73 // c0 = c5 = 0.0;
74 // c1 = c4 = -1.0/16.0;
75 // c2 = c3 = 9.0/16.0;
76 // } else if (a_interp_order == 6) {
77  c0 = c5 = 3.0/256.0;
78  c1 = c4 = -25.0/256.0;
79  c2 = c3 = 150.0/256.0;
80 // } else {
81 // fprintf(stderr,"Invalid value of interpolation order!\n");
82 // return 1;
83 // }
84 
85  /* create bounds for the transverse loop, i.e., to loop over
86  * all 1D lines along dimension "a_dir" */
87  int bounds_transverse[a_ndims];
88  _ArrayCopy1D_(a_dim_src, bounds_transverse, a_ndims);
89  bounds_transverse[a_dir] = 1;
90 
91  int index_transverse[a_ndims], done = 0;
92  _ArraySetValue_(index_transverse, a_ndims, 0);
93  while (!done) {
94 
95  int index_dst[a_ndims], index_src[a_ndims];
96  _ArrayCopy1D_(index_transverse, index_dst, a_ndims);
97  _ArrayCopy1D_(index_transverse, index_src, a_ndims);
98 
99  for (int i_dst = 0; i_dst < n_dst; i_dst++) {
100 
101  int i_m1 = i_dst*stride + (stride/2-1);
102  int i_m3 = i_m1 - 2;
103  int i_m2 = i_m1 - 1;
104  int i_p1 = i_m1 + 1;
105  int i_p2 = i_m1 + 2;
106  int i_p3 = i_m1 + 3;
107 
108  int p;
109  index_dst[a_dir] = i_dst;
110  _ArrayIndex1D_(a_ndims, a_dim_dst, index_dst, a_ngpt, p);
111 
112  int p_m3;
113  index_src[a_dir] = i_m3;
114  _ArrayIndex1D_(a_ndims, a_dim_src, index_src, a_ngpt, p_m3);
115 
116  int p_m2;
117  index_src[a_dir] = i_m2;
118  _ArrayIndex1D_(a_ndims, a_dim_src, index_src, a_ngpt, p_m2);
119 
120  int p_m1;
121  index_src[a_dir] = i_m1;
122  _ArrayIndex1D_(a_ndims, a_dim_src, index_src, a_ngpt, p_m1);
123 
124  int p_p1;
125  index_src[a_dir] = i_p1;
126  _ArrayIndex1D_(a_ndims, a_dim_src, index_src, a_ngpt, p_p1);
127 
128  int p_p2;
129  index_src[a_dir] = i_p2;
130  _ArrayIndex1D_(a_ndims, a_dim_src, index_src, a_ngpt, p_p2);
131 
132  int p_p3;
133  index_src[a_dir] = i_p3;
134  _ArrayIndex1D_(a_ndims, a_dim_src, index_src, a_ngpt, p_p3);
135 
136  for (int v = 0; v < a_nvars; v++) {
137  double val = c0 * a_u_src[p_m3*a_nvars+v]
138  + c1 * a_u_src[p_m2*a_nvars+v]
139  + c2 * a_u_src[p_m1*a_nvars+v]
140  + c3 * a_u_src[p_p1*a_nvars+v]
141  + c4 * a_u_src[p_p2*a_nvars+v]
142  + c5 * a_u_src[p_p3*a_nvars+v];
143  a_u_dst[p*a_nvars+v] = val;
144  }
145 
146  }
147 
148  _ArrayIncrementIndex_(a_ndims, bounds_transverse, index_transverse, done);
149 
150  }
151 
152  return 0;
153 }
absolute
#define absolute(a)
Definition: math_ops.h:32
_MACHINE_ZERO_
#define _MACHINE_ZERO_
Definition: basic.h:26
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319

◆ refine1D()

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 
)
static

Refine along a given dimension - the source and destination must have the same sizes along all other dimensions.

Note that the arrays must not have any ghost points!

Currently this function can only handle refinement factors that are integer powers of 2.

Parameters
a_dim_srcGrid size of source data
a_dim_dstGrid size of destination data
a_u_srcSource solution
a_u_dstDestination solution
a_dirDimension along which to coarsen
a_nvarsNumber of vector components of the solution
a_ngptNumber of ghost points
a_ndimsNumber of spatial dimensions

Definition at line 163 of file InterpolateGlobalnDVar.c.

172 {
173  for (int d = 0; d < a_ndims; d++) {
174  if ((d != a_dir) && (a_dim_src[d] != a_dim_dst[d])) {
175  fprintf(stderr, "Error in refine1D() -\n");
176  fprintf(stderr, " a_dim_src[%d] != a_dim_dst[%d]\n", d, d);
177  return 1;
178  }
179  }
180 
181  int n_src = a_dim_src[a_dir];
182  int n_dst = a_dim_dst[a_dir];
183  if (n_dst < n_src) {
184  fprintf(stderr, "Error in refine1D() -\n");
185  fprintf(stderr, " destination grid is coarser than source grid along a_dir!\n");
186  return 1;
187  }
188 
189  double fac = ((double) n_dst) / ((double) n_src);
190  int stride = (int) fac;
191  if (absolute(((double)stride)-fac) > _MACHINE_ZERO_) {
192  fprintf(stderr, "Error in refine1D() -\n");
193  fprintf(stderr, " non-integer refinement factor!\n");
194  return 1;
195  }
196 
197  /* create bounds for the transverse loop, i.e., to loop over
198  * all 1D lines along dimension "a_dir" */
199  int bounds_transverse[a_ndims];
200  _ArrayCopy1D_(a_dim_src, bounds_transverse, a_ndims);
201  bounds_transverse[a_dir] = 1;
202 
203  int index_transverse[a_ndims], done = 0;
204  _ArraySetValue_(index_transverse, a_ndims, 0);
205  while (!done) {
206 
207  int index_dst [a_ndims],
208  index_src0[a_ndims],
209  index_src1[a_ndims],
210  index_src2[a_ndims],
211  index_src3[a_ndims],
212  index_src4[a_ndims],
213  index_src5[a_ndims];
214  _ArrayCopy1D_(index_transverse, index_dst , a_ndims);
215  _ArrayCopy1D_(index_transverse, index_src0, a_ndims);
216  _ArrayCopy1D_(index_transverse, index_src1, a_ndims);
217  _ArrayCopy1D_(index_transverse, index_src2, a_ndims);
218  _ArrayCopy1D_(index_transverse, index_src3, a_ndims);
219  _ArrayCopy1D_(index_transverse, index_src4, a_ndims);
220  _ArrayCopy1D_(index_transverse, index_src5, a_ndims);
221 
222  for (int i_dst = 0; i_dst < n_dst; i_dst++) {
223 
224  double xi_dst = ((double) i_dst + 0.5) / ((double) stride) - 0.5;
225 
226  int i_src_2 = floor(xi_dst);
227  int i_src_3 = ceil(xi_dst);
228  int i_src_0 = i_src_2 - 2;
229  int i_src_1 = i_src_2 - 1;
230  int i_src_4 = i_src_3 + 1;
231  int i_src_5 = i_src_3 + 2;
232 
233  double alpha = (xi_dst - (double)i_src_2) / ((double)i_src_3 - (double)i_src_2);
234 
235  index_dst[a_dir] = i_dst;
236  int p; _ArrayIndex1D_(a_ndims, a_dim_dst, index_dst, a_ngpt, p);
237 
238  index_src0[a_dir] = i_src_0;
239  int q0; _ArrayIndex1D_(a_ndims, a_dim_src, index_src0, a_ngpt, q0);
240 
241  index_src1[a_dir] = i_src_1;
242  int q1; _ArrayIndex1D_(a_ndims, a_dim_src, index_src1, a_ngpt, q1);
243 
244  index_src2[a_dir] = i_src_2;
245  int q2; _ArrayIndex1D_(a_ndims, a_dim_src, index_src2, a_ngpt, q2);
246 
247  index_src3[a_dir] = i_src_3;
248  int q3; _ArrayIndex1D_(a_ndims, a_dim_src, index_src3, a_ngpt, q3);
249 
250  index_src4[a_dir] = i_src_4;
251  int q4; _ArrayIndex1D_(a_ndims, a_dim_src, index_src4, a_ngpt, q4);
252 
253  index_src5[a_dir] = i_src_5;
254  int q5; _ArrayIndex1D_(a_ndims, a_dim_src, index_src5, a_ngpt, q5);
255 
256  /* set interpolation coefficients depending on desired order */
257  double c0, c1, c2, c3, c4, c5;
258 // if (a_interp_order == 2) {
259 // c0 = 0.0;
260 // c1 = 0.0;
261 // c2 = (1.0-alpha);
262 // c3 = alpha;
263 // c4 = 0.0;
264 // c5 = 0.0;
265 // } else if (a_interp_order == 4) {
266 // c0 = 0.0;
267 // c1 = -((-2.0 + alpha)*(-1.0 + alpha)*alpha)/6.0;
268 // c2 = ((-2.0 + alpha)*(-1.0 + alpha)*(1.0 + alpha))/2.0;
269 // c3 = (alpha*(2.0 + alpha - alpha*alpha))/2.0;
270 // c4 = (alpha*(-1.0 + alpha*alpha))/6.0;
271 // c5 = 0.0;
272 // } else if (a_interp_order == 6) {
273  c0 = -((-3.0 + alpha)*(-2.0 + alpha)*(-1.0 + alpha)*alpha*(1.0 + alpha))/120.0;
274  c1 = ((-3.0 + alpha)*(-2.0 + alpha)*(-1.0 + alpha)*alpha*(2.0 + alpha))/24.0;
275  c2 = -((-3.0 + alpha)*(-2.0 + alpha)*(-1.0 + alpha)*(1.0 + alpha)*(2.0 + alpha))/12.0;
276  c3 = ((-3.0 + alpha)*(-2.0 + alpha)*alpha*(1.0 + alpha)*(2.0 + alpha))/12.0;
277  c4 = -((-3.0 + alpha)*(-1.0 + alpha)*alpha*(1.0 + alpha)*(2.0 + alpha))/24.0;
278  c5 = (alpha*(4.0 - 5.0*alpha*alpha + alpha*alpha*alpha*alpha))/120.0;
279 // } else {
280 // fprintf(stderr,"Invalid value of interpolation order!\n");
281 // return 1;
282 // }
283 
284  for (int v = 0; v < a_nvars; v++) {
285 
286  a_u_dst[p*a_nvars+v] = c0 * a_u_src[q0*a_nvars+v]
287  + c1 * a_u_src[q1*a_nvars+v]
288  + c2 * a_u_src[q2*a_nvars+v]
289  + c3 * a_u_src[q3*a_nvars+v]
290  + c4 * a_u_src[q4*a_nvars+v]
291  + c5 * a_u_src[q5*a_nvars+v];
292 
293  }
294 
295  }
296 
297  _ArrayIncrementIndex_(a_ndims, bounds_transverse, index_transverse, done);
298 
299  }
300 
301  return 0;
302 }
absolute
#define absolute(a)
Definition: math_ops.h:32
_MACHINE_ZERO_
#define _MACHINE_ZERO_
Definition: basic.h:26
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319

◆ InterpolateGlobalnDVar()

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 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 }
fillGhostCells
void fillGhostCells(const int *const, const int, double *const, const int, const int, const int *const)
Definition: FillGhostCells.c:17
refine1D
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)
Definition: InterpolateGlobalnDVar.c:163
coarsen1D
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)
Definition: InterpolateGlobalnDVar.c:32
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
isPowerOfTwo
static int isPowerOfTwo(int x)
Definition: InterpolateGlobalnDVar.c:13