Streamer fluid modeling - An overview of ARCoS  1.0
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
interpol2.c
Go to the documentation of this file.
1 
28 #include <stdlib.h>
29 #include <stdio.h>
30 #include <stdarg.h>
31 #include <math.h>
32 
33 #include "grid.h"
34 #include "interpol2.h"
35 #include "parameters.h"
36 #include "proto.h"
37 #include "rz_array.h"
38 #include "species.h"
39 
40 #define anm(N_, M_) this->stencil[(N_) * this->method->q + (M_)]
41 
42 static double gen_apply (interpol_t *this, double r, double z);
43 void bilin_set_coeffs (interpol_t *this);
44 void quadlog_set_coeffs (interpol_t *this);
45 double quadlog_apply (interpol_t *this, double r, double z);
46 
48 interpol_t*
49 interpol_new_a (double Lr, double Lz, interpol_method_t *method)
50 {
51  interpol_t *interpol;
52 
53  interpol = xmalloc (sizeof (interpol_t));
54 
55  interpol->method = method;
56  interpol->Lr = Lr;
57  interpol->Lz = Lz;
58 
59  /* The number of pairs {i, j} (i, j >= 0) such that i + j < p
60  is p (p + 1) / 2 */
61  interpol->coeffs = xmalloc (sizeof(double) *
62  method->p * (method->p + 1) / 2);
63  interpol->stencil = xmalloc (sizeof(double) * method->q * method->q);
64 
65  return interpol;
66 }
67 
69 void
70 interpol_free (interpol_t *this)
71 {
72  free (this->coeffs);
73  free (this->stencil);
74  free (this);
75 }
76 
82 void
83 interpol_set_stencil (interpol_t *this, double r0, double z0, ...)
84 {
85  int i, j, indx;
86  va_list ap;
87 
88  va_start (ap, z0);
89 
90  for (indx = 0, i = 0; i < this->method->q; i++) {
91  for (j = 0; j < this->method->q; j++) {
92  this->stencil[indx++] = va_arg (ap, double);
93  }
94  }
95 
96  this->r0 = r0;
97  this->z0 = z0;
98 
99  va_end (ap);
100 
101  /* If the method has an optimized routine, use it. */
102  if (this->method->set_coeffs != NULL) {
103  this->method->set_coeffs (this);
104  }
105  else {
106  interpol_set_coeffs (this);
107  }
108 }
109 
115 void
116 interpol_set_stencil_at (grid_t *grid, interpol_t *this, double r0, double z0,
117  rz_array_t *ar,
118  int ir, int iz, int itheta)
119 {
120  int indx, i, j;
121  REAL *sp, *ap;
122 
123  assert (ar->ntheta > itheta);
124 
125  sp = this->stencil;
126  for (indx = 0, i = 0; i < this->method->q; i++) {
127  double r;
128  if (this->method->masses) {
129  r = cyl_r_at (ir + i - this->method->ir0, grid->level);
130  } else {
131  r = 1;
132  }
133 
134  ap = RZTP (ar, ir + i - this->method->ir0,
135  iz - this->method->iz0, itheta);
136 
137  for (j = 0; j < this->method->q; j++, ap += ar->strides[Z_INDX]) {
138  /* The case ir = -1 does not give any problem, since we set
139  up the correct boundaries. */
140  *sp = r * (*ap);
141  sp++;
142  }
143  }
144 
145  this->r0 = r0;
146  this->z0 = z0;
147 
148  /* If the method has an optimized routine, use it. */
149  if (this->method->set_coeffs != NULL) {
150  this->method->set_coeffs (this);
151  }
152  else {
153  interpol_set_coeffs (this);
154  }
155 }
156 
160 void
161 interpol_set_coeffs (interpol_t *this)
162 {
163  double Lrn[MAX_ORDER], Lzn[MAX_ORDER];
164  double *matrix, sum;
165  int n, ip, j;
166  int pindx, mindx;
167 
168  matrix = this->method->matrix;
169 
170  Lrn[0] = 1.0;
171  Lzn[0] = 1.0;
172 
173  for (ip = 0, pindx = 0, mindx = 0; ip < this->method->p; ip++) {
174  Lrn[ip + 1] = Lrn[ip] * this->Lr;
175  Lzn[ip + 1] = Lzn[ip] * this->Lz;
176  for (n = 0; n <= ip; n++) {
177  sum = 0;
178 
179  /* This is the coefficient of r^n z^m, with m + n = ip */
180  for (j = 0; j < this->method->q * this->method->q; j++)
181  sum += matrix[mindx++] * this->stencil[j];
182 
183  this->coeffs[pindx++] = sum / Lrn[n] / Lzn[ip - n];
184  }
185  }
186 }
187 
189 double
190 interpol_apply (interpol_t *this, double r, double z)
191 {
192  double result;
193 
194  /* If the method has an optimized routine, call it and return */
195  if (this->method->apply != NULL) {
196  result = this->method->apply (this, r, z);
197  } else {
198  /* Else, apply the general method. */
199  result = gen_apply (this, r, z);
200  }
201 
202  if (this->method->masses) {
203  return result / cyl_q(r);
204  } else {
205  return result;
206  }
207 }
208 
210 static double
211 gen_apply (interpol_t *this, double r, double z)
212 {
213  double deltar, deltaz, sum;
214 
215  deltar = r - this->r0;
216  deltaz = z - this->z0;
217 
218  int ipmax=this->method->p;
219  if (ipmax < 3) {
220  fprintf(stdout,"ipmax=%d\n",ipmax);
221  fatal("interpol2: order interpolation not equal to 3\n");
222  }
223 
224  if (deltar==deltaz) {
225  sum = this->coeffs[0] +
226  deltar * (this->coeffs[1] +
227  this->coeffs[2] +
228  deltar * (this->coeffs[3] +
229  this->coeffs[4] +
230  this->coeffs[5] ));
231 
232  } else
233  {
234  sum = this->coeffs[0] +
235  deltaz * (this->coeffs[1] + deltaz * this->coeffs[3]) +
236  deltar * (this->coeffs[2] +
237  deltaz * this->coeffs[4] +
238  deltar * this->coeffs[5] );
239  }
240  return sum;
241 }
242 
246 double matrix_zero[] = {1.0};
247 
248 interpol_method_t interpol_zero = {matrix_zero, 1, 1, 0, 0, FALSE,
249  interpol_set_coeffs,
250  gen_apply};
251 
252 interpol_method_t interpol_zero_masses = {matrix_zero, 1, 1, 0, 0, TRUE,
253  interpol_set_coeffs,
254  gen_apply};
255 
256 /* You are not supposed to understand these numbers.
257  Refer to InterpolArrays.nb */
258 
260 double matrix_quadratic[] =
261  {-0.1111111111111111, 0.2222222222222222, -0.1111111111111111,
262  0.2222222222222222, 0.5555555555555556, 0.2222222222222222,
263  -0.1111111111111111, 0.2222222222222222, -0.1111111111111111,
264  -0.1666666666666667, 0 , 0.1666666666666667,
265  -0.1666666666666667, 0 , 0.1666666666666667,
266  -0.1666666666666667, 0 , 0.1666666666666667,
267  -0.1666666666666667, -0.1666666666666667, -0.1666666666666667,
268  0 , 0 , 0 ,
269  0.1666666666666667, 0.1666666666666667, 0.1666666666666667,
270  0.1666666666666667, -0.3333333333333333, 0.1666666666666667,
271  0.1666666666666667, -0.3333333333333333, 0.1666666666666667,
272  0.1666666666666667, -0.3333333333333333, 0.1666666666666667,
273  0.25 , 0 , -0.25 ,
274  0.25 , 0 , -0.25 ,
275  0 , 0 , 0 ,
276  -0.25 , 0 , 0.25 ,
277  0.1666666666666667, 0.1666666666666667, 0.1666666666666667,
278  -0.3333333333333333, -0.3333333333333333, -0.3333333333333333,
279  0.1666666666666667, 0.1666666666666667, 0.1666666666666667
280  };
281 
282 interpol_method_t interpol_quadratic = {matrix_quadratic, 3, 3, 1, 1, FALSE,
283  interpol_set_coeffs,
284  gen_apply};
285 
286 interpol_method_t interpol_quadratic_masses = {matrix_quadratic, 3, 3, 1, 1,
287  TRUE,
288  interpol_set_coeffs,
289  gen_apply};
290 
293 double matrix_wackers[] =
294  { 0 , 0 , 0 ,
295  0 , 1.0 , 0 ,
296  0 , 0 , 0 ,
297  -0.1666666666666667, 0 , 0.1666666666666667,
298  -0.1666666666666667, 0 , 0.1666666666666667,
299  -0.1666666666666667, 0 , 0.1666666666666667,
300  -0.1666666666666667, -0.1666666666666667, -0.1666666666666667,
301  0 , 0 , 0 ,
302  0.1666666666666667, 0.1666666666666667, 0.1666666666666667,
303  0.1 , -0.2 , 0.1 ,
304  0.3 , -0.6 , 0.3 ,
305  0.1 , -0.2 , 0.1 ,
306  0.25 , 0 , -0.25 ,
307  0 , 0 , 0 ,
308  -0.25 , 0 , 0.25 ,
309  0.1 , 0.3 , 0.1 ,
310  -0.2 , -0.6 , -0.2 ,
311  0.1 , 0.3 , 0.1 };
312 
313 interpol_method_t interpol_wackers = {matrix_wackers, 3, 3, 1, 1, FALSE,
314  interpol_set_coeffs,
315  gen_apply};
316 
317 interpol_method_t interpol_wackers_masses = {matrix_wackers, 3, 3, 1, 1, TRUE,
318  interpol_set_coeffs,
319  gen_apply};
320 
321 
322 double matrix_averaged[] =
323 { -0.01880341880341880 , -0.008547008547008547, -0.01880341880341880 ,
324  -0.008547008547008547, 1.109401709401709 , -0.008547008547008547,
325  -0.01880341880341880 , -0.008547008547008547, -0.01880341880341880 ,
326  -0.1666666666666667 , 0 , 0.1666666666666667 ,
327  -0.1666666666666667 , 0 , 0.1666666666666667 ,
328  -0.1666666666666667 , 0 , 0.1666666666666667 ,
329  -0.1666666666666667 , -0.1666666666666667 , -0.1666666666666667 ,
330  0 , 0 , 0 ,
331  0.1666666666666667 , 0.1666666666666667 , 0.1666666666666667 ,
332  0.1128205128205128 , -0.1987179487179487 , 0.1128205128205128 ,
333  0.3012820512820513 , -0.6564102564102564 , 0.3012820512820513 ,
334  0.1128205128205128 , -0.1987179487179487 , 0.1128205128205128 ,
335  0.25 , 0 , -0.25 ,
336  0 , 0 , 0 ,
337  -0.25 , 0 , 0.25 ,
338  0.1128205128205128 , 0.3012820512820513 , 0.1128205128205128 ,
339  -0.1987179487179487 , -0.6564102564102564 , -0.1987179487179487 ,
340  0.1128205128205128 , 0.3012820512820513 , 0.1128205128205128 };
341 
342 interpol_method_t interpol_averaged = {matrix_averaged, 3, 3, 1, 1, FALSE,
343  interpol_set_coeffs,
344  gen_apply};
345 
346 interpol_method_t interpol_averaged_masses = {matrix_averaged, 3, 3, 1, 1, TRUE,
347  interpol_set_coeffs,
348  gen_apply};
349 
350 /* 4-point bilinear interpolation. */
351 double matrix_bilin[] =
352  { 1.0, 0.0, 0.0, 0.0,
353  -1.0, 0.0, 1.0, 0.0,
354  -1.0, 1.0, 0.0, 0.0,
355  0.0, 0.0, 0.0, 0.0,
356  1.0, -1.0, -1.0, 1.0,
357  0.0, 0.0, 0.0, 0.0};
358 
362 void
363 bilin_set_coeffs (interpol_t *this)
364 {
365  double Lr, Lz;
366  double *c;
367 
368  c = this->coeffs;
369 
370  Lr = this->Lr;
371  Lz = this->Lz;
372 
373  c[0] = anm(0, 0);
374  c[1] = (anm(0, 1) - anm(0, 0)) / Lz;
375  c[2] = (anm(1, 0) - anm(0, 0)) / Lr;
376  c[3] = 0; /* Never used anyway. */
377  c[4] = (anm(0, 0) - anm(1, 0) - anm(0, 1) + anm(1, 1)) / Lr / Lz;
378  c[5] = 0; /* Never used anyway. */
379 }
380 
381 double
382 bilin_apply (interpol_t *this, double r, double z)
383 {
384  double d_r, d_z;
385  double *c;
386 
387  c = this->coeffs;
388 
389  d_r = r - this->r0;
390  d_z = z - this->z0;
391 
392  return c[0] + d_z * (c[1] + c[4] * d_r) + c[2] * d_r;
393 }
394 
395 
396 interpol_method_t interpol_bilin = {matrix_bilin,
397  3, 2,
398  0, 0,
399  FALSE,
400  bilin_set_coeffs,
401  bilin_apply};
402 
403 interpol_method_t interpol_bilin_masses = {matrix_bilin,
404  3, 2,
405  0, 0, TRUE,
406  bilin_set_coeffs,
407  bilin_apply};
408 
409 void
410 quadlog_set_coeffs (interpol_t *this)
411 {
412  int indx, i, j;
413 
414  for (indx = 0, i = 0; i < this->method->q; i++) {
415  for (j = 0; j < this->method->q; j++) {
416  this->stencil[indx] = log (this->stencil[indx]);
417  indx++;
418  }
419  }
420  interpol_set_coeffs (this);
421 }
422 
423 double
424 quadlog_apply (interpol_t *this, double r, double z)
425 {
426  double rlog;
427 
428  rlog = gen_apply (this, r, z);
429  return exp (rlog);
430 }
431 
432 interpol_method_t interpol_quadlog = {matrix_quadratic, 3, 3, 1, 1,
433  FALSE,
434  quadlog_set_coeffs,
435  quadlog_apply};
436 
437 interpol_method_t *interpol_methods_index[MAX_INTERPOL_METHODS] = {
438  &interpol_zero_masses,
439  &interpol_quadratic_masses,
440  &interpol_wackers_masses,
441  &interpol_averaged_masses,
442  &interpol_zero,
443  &interpol_quadratic,
444  &interpol_wackers,
445  &interpol_averaged,
446  &interpol_quadlog };
Generated on Fri Jul 4 2014 15:40:19 for Streamer fluid modeling - An overview of ARCoS by   doxygen 1.8.2