HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
InterpolateGlobalnDVar.c
Go to the documentation of this file.
1 
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <math.h>
9 #include <arrayfunctions.h>
10 #include <mathfunctions.h>
11 
13 static int isPowerOfTwo(int x)
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 }
23 
32 static int coarsen1D( const int* const a_dim_src,
33  const int* const a_dim_dst,
34  const double* const a_u_src,
35  double* const a_u_dst,
36  const int a_dir,
37  const int a_nvars,
38  const int a_ngpt,
39  const int a_ndims
40  )
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 }
154 
163 static int refine1D(const int* const a_dim_src,
164  const int* const a_dim_dst,
165  const double* const a_u_src,
166  double* const a_u_dst,
167  const int a_dir,
168  const int a_nvars,
169  const int a_ngpt,
170  const int a_ndims
171  )
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 }
303 
317 int InterpolateGlobalnDVar( const int* const a_dim_dst,
318  double** const a_u_dst,
319  const int* const a_dim_src,
320  double* const a_u_src,
321  const int a_nvars,
322  const int a_ghosts,
323  const int a_ndims,
324  const int* const a_periodic
325  )
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 }
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
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 _ArrayIndex1D_(N, imax, i, ghost, index)
int InterpolateGlobalnDVar(const int *const, double **const, const int *const, double *const, const int, const int, const int, const int *const)
Contains function definitions for common mathematical functions.
#define _ArrayCopy1D_(x, y, size)
#define _MACHINE_ZERO_
Definition: basic.h:26
Contains macros and function definitions for common array operations.
#define absolute(a)
Definition: math_ops.h:32
static int isPowerOfTwo(int x)