PLUTO
ppm_coeffs.v4.1.c
Go to the documentation of this file.
1 #include "pluto.h"
2 
3 #if INTERPOLATION == WENO3
4  #undef PPM_VERSION
5  #define PPM_VERSION PPM3
6 #endif
7 #define PPM_VERSION PPM_ORDER
8 #define PPM3 3
9 #define PPM4 4
10 #define PPM5 5
11 
12 /* -----------------------------------------------------
13  Static local variables
14  ----------------------------------------------------- */
15 
16 static double ***Wp3D, ***Wm3D, **hp3D, **hm3D;
17 
18 
19 static void PPM_CartCoeff(double **, double **, int, int);
20 static void PPM_FindWeights(double **, double **, int, Grid *);
21 
22 /* ********************************************************************* */
23 void PPM_CoefficientsSet(Grid *grid)
24 /*!
25  * Compute interpolation coefficients for PPM in every coordinate system.
26  *
27  *********************************************************************** */
28 {
29  int i, iL, iR, n, beg, end, d;
30  int uniform_grid[3];
31  double i1, i2, rp, dr, den;
32  double **wp, **wm;
33 
34  if (Wp3D == NULL){
35  Wp3D = ArrayBox(0, DIMENSIONS-1, 0, NMAX_POINT-1, -2, 2);
36  Wm3D = ArrayBox(0, DIMENSIONS-1, 0, NMAX_POINT-1, -2, 2);
37  hp3D = ARRAY_2D(DIMENSIONS, NMAX_POINT, double);
38  hm3D = ARRAY_2D(DIMENSIONS, NMAX_POINT, double);
39  }
40 
41 /* -------------------------------------------------------
42  Check whether the grid is regularly spaced.
43  This is always true when chombo is used, but not
44  necessarily so with the static grid version.
45  ------------------------------------------------------- */
46 
47  #ifdef CHOMBO
48  for (d = 0; d < DIMENSIONS; d++) uniform_grid[d] = 1;
49  #else
50  for (d = 0; d < DIMENSIONS; d++) {
51  uniform_grid[d] = 1;
52  for (i = 0; i < grid[d].np_tot-1; i++){
53  if (fabs(grid[d].dx[i+1]/grid[d].dx[i]-1.0) > 1.e-9){
54  uniform_grid[d] = 0;
55  break;
56  }
57  }
58  }
59  #endif
60 
61 /* -----------------------------------------------------
62  set the local (iL, iR) and global (beg, end)
63  reconstruction stencils and the order of accuracy n
64  ----------------------------------------------------- */
65 
66  #if PPM_VERSION == PPM3
67  iL = 1; iR = 1;
68  #elif PPM_VERSION == PPM4
69  iL = 1; iR = 2;
70  #elif PPM_VERSION == PPM5
71  iL = 2; iR = 2;
72  #endif
73  n = iR + iL + 1;
74 
75 /* ----------------------------------------------------------
76  Loop over dimensions
77  ---------------------------------------------------------- */
78 
79  for (d = 0; d < DIMENSIONS; d++){
80 
81  wp = Wp3D[d];
82  wm = Wm3D[d];
83 
84  beg = iL;
85  end = grid[d].np_tot - 1 - iR;
86 
87  /* -------------------------------------------------------------
88  Interpolation weights for Cartesian geometry (called anyway
89  to initialize cofficients to zero where necessary).
90  ------------------------------------------------------------- */
91 
92  PPM_CartCoeff(wp, wm, beg, end);
93 
94  /* -------------------------------------------------------------
95  Interpolation weights for cylindrical/polar geometries
96  ------------------------------------------------------------- */
97 
98  #if GEOMETRY == CYLINDRICAL || GEOMETRY == POLAR
99  if (d == IDIR) for (i = beg; i <= end; i++){
100  rp = grid[IDIR].xr[i];
101  dr = grid[IDIR].dx[i];
102  i1 = rp/dr;
103 // i1 = (double)(i + (grid[dir].beg-IBEG) - IBEG + 1);
104  i2 = i1*i1;
105  #if PPM_VERSION == PPM3
106  den = 12.0*POLY_3(1.0, -1.0, -3.0, 2.0, i1);
107 
108  wp[i][-2] = 0.0;
109  wp[i][-1] = POLY_3( -3.0, 2.0, 6.0, -4.0, i1)/den;
110  wp[i][ 0] = POLY_3( 11.0, -13.0, -28.0, 20.0, i1)/den;
111  wp[i][ 1] = POLY_3( 4.0, -1.0, -14.0, 8.0, i1)/den;
112  wp[i][ 2] = 0.0;
113 
114  wm[i][-2] = 0.0;
115  wm[i][-1] = POLY_3( 3.0, -5.0, -10.0, 8.0, i1)/den;
116  wm[i][ 0] = POLY_3( 10.0, -9.0, -32.0, 20.0, i1)/den;
117  wm[i][ 1] = POLY_3( -1.0, 2.0, 6.0, -4.0, i1)/den;
118  wm[i][ 2] = 0.0;
119  #elif PPM_VERSION == PPM4
120  den = 24.0*POLY_2(4.0, -15.0, 5.0, i2);
121 
122  wp[i][-2] = 0.0;
123  wp[i][-1] = POLY_4(-12.0, -1.0, 30.0, -1.0, -10.0, i1)/den;
124  wp[i][ 0] = POLY_4( 60.0, -27.0, -210.0, 13.0, 70.0, i1)/den;
125  wp[i][ 1] = POLY_4( 60.0, 27.0, -210.0, -13.0, 70.0, i1)/den;
126  wp[i][ 2] = POLY_4(-12.0, 1.0, 30.0, 1.0, -10.0, i1)/den;
127 
128  #elif PPM_VERSION == PPM5
129  den = 120.0*(2.0*i1-1.0)*POLY_4(12.0, 16.0, -13.0, -6.0, 3.0, i1);
130 
131  wp[i][-2] = POLY_5( -80, 32, 200, -80, -60, 24, i1)/den;
132  wp[i][-1] = POLY_5( 492, -193, -1230, 535, 384, -156, i1)/den;
133  wp[i][ 0] = POLY_5(-1276, 1157, 4090, -2075, -1332, 564, i1)/den;
134  wp[i][ 1] = POLY_5( -684, 27, 2610, -885, -888, 324, i1)/den;
135  wp[i][ 2] = POLY_5( 108, -63, -270, 105, 96, -36, i1)/den;
136 
137  wm[i][-2] = POLY_5( 60, -84, -261, 129, 84, -36, i1)/den;
138  wm[i][-1] = POLY_5( -504, 660, 2133, -1197, -732, 324, i1)/den;
139  wm[i][ 0] = POLY_5(-1128, 604, 4487, -1763, -1488, 564, i1)/den;
140  wm[i][ 1] = POLY_5( 168, -292, -1119, 511, 396, -156, i1)/den;
141  wm[i][ 2] = POLY_5( -36, 72, 160, -80, -60, 24, i1)/den;
142  #endif
143  }
144  #endif
145 
146  /* -----------------------------------------------------
147  Interpolation weights for spherical geometry
148  ------------------------------------------------------- */
149 
150  #if GEOMETRY == SPHERICAL
151  if (d == IDIR) for (i = beg; i <= end; i++){
152  rp = grid[IDIR].xr[i];
153  dr = grid[IDIR].dx[i];
154  i1 = fabs(rp/dr);
155 /* i1 = (double)(i + (grid[dir].beg-IBEG) - IBEG + 1); */
156  i2 = i1*i1;
157  #if PPM_VERSION == PPM3
158  den = 18.0*POLY_6(4.0, -6.0, -9.0, 20.0, 15.0, -30.0, 10.0, i1);
159 
160  wp[i][-2] = 0.0;
161  wp[i][-1] = -POLY_2(7.0, -9.0, 3.0, i1)
162  *POLY_2(3.0, -9.0, 10.0, i2)/den;
163  wp[i][ 0] = POLY_2(1.0, -3.0, 3.0, i1)
164  *POLY_4(69.0, 96.0, -63.0, -90.0, 50.0, i1)/den;
165  wp[i][ 1] = 2.0*POLY_2( 1.0, 3.0, 3.0, i1)
166  *POLY_4(12.0, -48.0, 72.0, -45.0, 10.0, i1)/den;
167  wp[i][ 2] = 0.0;
168 
169  wm[i][-2] = 0.0;
170  wm[i][-1] = 2.0*POLY_2(7.0, -9.0, 3.0, i1)
171  *POLY_4(1.0, -1.0, -3.0, 5.0, 10.0, i1)/den;
172  wm[i][ 0] = POLY_2(1.0, -3.0, 3.0, i1)
173  *POLY_4(62.0, 100.0, -33.0, -110.0, 50.0, i1)/den;
174  wm[i][ 1] = -POLY_2(1.0, 3.0, 3.0, i1)
175  *POLY_4(4.0, -22.0, 51.0, -40.0, 10.0, i1)/den;
176  wm[i][ 2] = 0.0;
177  #elif PPM_VERSION == PPM4
178  den = 36.0*POLY_4(16.0, -60.0, 150.0, -85.0, 15.0, i2);
179 
180  wp[i][-2] = 0.0;
181  wp[i][-1] = -POLY_2( 7, -9, 3, i1)/den
182  *POLY_6(12, 16, -30, -48.0, 23, 48, 15, i1);
183  wp[i][ 0] = POLY_2( 1, -3, 3, i1)/den
184  *POLY_6(372, 1008.0, 510, -720, -487, 144, 105, i1);
185  wp[i][ 1] = POLY_2( 1, 3, 3.0, i1)/den
186  *POLY_6(372, -1008, 510, 720, -487, -144, 105, i1);
187  wp[i][ 2] = -POLY_2( 7, 9, 3, i1)/den
188  *POLY_6(12, -16, -30, 48, 23, -48, 15, i1);
189  #elif PPM_VERSION == PPM5
190  den = POLY_10( 48.0, -48.0, -164.0, 200.0, 390.0,
191  -399.0, -161.0, 210.0, 0.0, -35.0, 7.0, i1);
192  den *= 180.0;
193 
194  wp[i][-2] = 2.0*POLY_2(19., -15., 3., i1)/den
195  *POLY_4(16., -60., 94., -45., 7., i2);
196 
197  wp[i][-1] = -POLY_2(7.0, -9.0, 3.0, i1)/den;
198  wp[i][-1] *= POLY_8(508., 240., -1740., -795., 2417.,
199  930., -780., -175., 91., i1);
200 
201  wp[i][0] = POLY_2(1.0, -3.0, 3.0, i1)/den;
202  wp[i][0] *= POLY_8(8132., 15120., -5700., -20325., 3863.,
203  8670., -1800., -1225., 329., i1);
204 
205  wp[i][1] = POLY_2(1.0, 3.0, 3.0, i1)/den;
206  wp[i][1] *= POLY_8(4212., -15120., 16560., 1275., -11517.,
207  4350., 1620., -1225., 189., i1);
208 
209  wp[i][2] = -POLY_2(7.0, 9.0, 3.0, i1)/den;
210  wp[i][2] *= POLY_8( 108., -240., -120., 645., -223.,
211  -510., 510., -175., 21., i1);
212 
213  wm[i][-2] = -POLY_2( 19., -15., 3., i1)/den;
214  wm[i][-2] *= POLY_8( 16., -16., -60., 96., 222.,
215  -51., -127., 7., 21., i1);
216 
217  wm[i][-1] = POLY_2( 7., -9., 3., i1)/den;
218  wm[i][-1] *= POLY_8( 344., -164., -1350., 1184., 4888.,
219  1071., -1663., -287., 189., i1);
220 
221  wm[i][0] = POLY_2(1.0, -3.0, 3.0, i1)/den;
222  wm[i][0] *= POLY_8(7064., 15196., -310., -21376., 368.,
223  9431., -1163., -1407., 329., i1);
224 
225  wm[i][1] = -POLY_2(1.0, 3.0, 3.0, i1)/den;
226  wm[i][1] *= POLY_8( 696., -3516., 6850., -1544., -4388.,
227  2329., 543., -553., 91., i1);
228 
229  wm[i][2] = 2.0*POLY_2(7.0, 9.0, 3.0, i1)/den;
230  wm[i][2] *= POLY_8( 12., -42., 25., 132., -91.,
231  -122., 151., -56., 7., i1);
232  #endif /* PPM_VERSION */
233  } else if (d == JDIR) {
234 
235  PPM_FindWeights(wp, wm, d, grid);
236 
237  } else if (d == KDIR){
238 
239  PPM_CartCoeff(wp, wm, beg, end);
240 
241  }
242  #endif /* GEOMETRY == SPHERICAL */
243 
244  /* ------------------------------------------------------
245  Compute Q6 coefficients
246  ------------------------------------------------------ */
247 
248  PPM_Q6_Coeffs(hp3D[d], hm3D[d], d, grid);
249 
250  } /* End main loop on dimensions */
251 
252 }
253 
254 double BetaTheta(double x);
255 static double gx0;
256 static int gk;
257 /* ********************************************************************* */
258 void PPM_FindWeights(double **wp, double **wm, int dir, Grid *grid)
259 /*!
260  * Find the coefficients numerically by inverting the beta matrix
261  * We define the matrix beta[0...n][0...j] so we need to
262  * index it as beta[k][j-jb]
263  *
264  *********************************************************************** */
265 {
266  int i, j, k, n;
267  int jb, je, iL, iR, beg, end;
268  int indx[16];
269  double rp, rc, rm, vol, d, a[16];
270  static double **beta;
271 
272 /* ----------------------------------------------------
273  set the reconstruction stencil & order of accuracy
274  ---------------------------------------------------- */
275 
276  #if PPM_VERSION == PPM3
277  iL = 1; iR = 1;
278  #elif PPM_VERSION == PPM4
279  iL = 1; iR = 2;
280  #elif PPM_VERSION == PPM5
281  iL = 2; iR = 2;
282  #endif
283  n = iR + iL + 1;
284 
285  beg = iL;
286  end = grid[dir].np_tot - 1 - iR;
287 
288  if (beta == NULL) beta = ARRAY_2D(8, 8, double);
289 
290 /* -----------------------------------------------------------------
291  Start main loop
292  ----------------------------------------------------------------- */
293 
294  for (i = beg; i <= end; i++){
295 
296  rc = grid[dir].x[i];
297  jb = i - iL; je = i + iR;
298 
299  if (dir == IDIR){
300  double m;
301  #if GEOMETRY == CARTESIAN
302  m = 0;
303  #elif GEOMETRY == CYLINDRICAL || GEOMETRY == POLAR
304  m = 1;
305  #elif GEOMETRY == SPHERICAL
306  m = 2;
307  #endif
308 
309  for (j = jb; j <= je; j++){ /* -- stencil loop -- */
310  rp = grid[dir].xr[j];
311  rm = grid[dir].xl[j];
312  vol = (pow(rp,m+1) - pow(rm,m+1))/(m+1.0);
313  for (k = 0; k < n; k++) { /* -- order loop -- */
314 /* beta[k][j-jb] = (pow(rp-ri,k+m+1) - pow(rm-ri,k+m+1))/(k+m+1.0)/vol; */
315 
316  #if GEOMETRY == CARTESIAN
317  beta[k][j-jb] = (pow(rp-rc,k+1) - pow(rm-rc,k+1))/(k+1.0)/vol;
318  #elif GEOMETRY == CYLINDRICAL || GEOMETRY == POLAR
319  beta[k][j-jb] = pow(rp-rc,k+1)*((k+1.0)*rp + rc)
320  -pow(rm-rc,k+1)*((k+1.0)*rm + rc);
321  beta[k][j-jb] /= (k+2.0)*(k+1.0)*vol;
322  #elif GEOMETRY == SPHERICAL
323  beta[k][j-jb] = pow(rp-rc,k+1)*((k*k + 3.0*k + 2.0)*rp*rp
324  + 2.0*rc*(k+1.0)*rp + 2.0*rc*rc)
325  -pow(rm-rc,k+1)*((k*k + 3.0*k + 2.0)*rm*rm
326  + 2.0*rc*(k+1.0)*rm + 2.0*rc*rc);
327 
328  beta[k][j-jb] /= (k+3.0)*(k+2.0)*(k+1.0)*vol;
329  #endif
330  }
331  }
332  } else if (dir == JDIR){
333  if (n > 5) {
334  print1 ("! Order of accuracy too large for jdir in PPM_FindWeights\n");
335  QUIT_PLUTO(1);
336  }
337  for (j = jb; j <= je; j++){ /* -- stencil loop -- */
338  #if GEOMETRY == SPHERICAL
339  double zp, zm;
340  double zp2, zm2;
341 
342  rp = grid[dir].xr[j];
343  rm = grid[dir].xl[j];
344  vol = cos(rm) - cos(rp);
345 
346  zp = rp - rc; zp2 = zp*zp;
347  zm = rm - rc; zm2 = zm*zm;
348 
349  /* Table of integral, page 216, equation 2.635 */
350 
351  beta[0][j-jb] = vol;
352  beta[1][j-jb] = (-zp*cos(rp) + sin(rp))
353  -(-zm*cos(rm) + sin(rm));
354  beta[2][j-jb] = ((2.0-zp2)*cos(rp) + 2.0*zp*sin(rp))
355  -((2.0-zm2)*cos(rm) + 2.0*zm*sin(rm));
356  beta[3][j-jb] = (zp*(6.0-zp2)*cos(rp) + 3.0*(zp2-2.0)*sin(rp))
357  -(zm*(6.0-zm2)*cos(rm) + 3.0*(zm2-2.0)*sin(rm));
358  beta[4][j-jb] = ((-24.0 + zp2*(12.0-zp2))*cos(rp)
359  + 4.0*zp*(zp2-6.0)*sin(rp))
360  -((-24.0 + zm2*(12.0 - zm2))*cos(rm)
361  + 4.0*zm*(zm2-6.0)*sin(rm));
362 
363 gx0 = rc;
364 for (gk = 0; gk < n; gk++){
365 double scrh;
366 scrh = GaussQuadrature(&BetaTheta, rm, rp, 1, 5);
367 /*
368 if (fabs(scrh/beta[gk][j-jb]-1.0) > 1.e-4){
369  printf ("gk = %d, I = %12.6e, beta = %12.6e, err = %12.6e\n",
370  gk, scrh, beta[gk][j-jb], fabs(scrh-beta[gk][j-jb]));
371  QUIT_PLUTO(1);
372 }
373 */
374 beta[gk][j-jb] = scrh;
375 
376 }
377 
378  beta[0][j-jb] /= vol;
379  beta[1][j-jb] /= vol;
380  beta[2][j-jb] /= vol;
381  beta[3][j-jb] /= vol;
382  beta[4][j-jb] /= vol;
383 
384  #else
385  #endif
386  }
387  } else {
388  print1 ("! PPM_FindWeights: not defined for x3 direction\n");
389  QUIT_PLUTO(1);
390  }
391 
392  LUDecompose(beta, n, indx, &d);
393 
394  rp = grid[dir].xr[i];
395  rm = grid[dir].xl[i];
396 
397  /* ---------------------------------------------------------
398  Solve B.w = xi^k
399  --------------------------------------------------------- */
400 
401  a[0] = 1.0;
402  for (k = 1; k < n; k++) a[k] = a[k-1]*(rp - rc);
403  LUBackSubst (beta, n, indx, a);
404  for (j = 0; j < n; j++) {
405 /*
406  if (fabs(a[j]/wp[i][-iL+j]-1.0) > 1.e-9){
407  printf ("! Err, i = %d, j = %d, Wp = %18.10e, a = %18.10e\n",
408  i, j, wp[i][-iL+j], a[j]);
409  QUIT_PLUTO(1);
410  }
411 */
412  wp[i][-iL+j] = a[j];
413  }
414 
415  #if PPM_VERSION != PPM4
416  a[0] = 1.0;
417  for (k = 1; k < n; k++) a[k] = a[k-1]*(rm - rc);
418  LUBackSubst (beta, n, indx, a);
419  for (j = 0; j < n; j++) {
420 /*
421  if (fabs(a[j]/wm[i][-iL+j]-1.0) > 1.e-9){
422  printf ("! Err, i = %d, j = %d, Wm = %12.6e, a = %12.6e\n",
423  i, j, wm[i][-iL+j], a[j]);
424  QUIT_PLUTO(1);
425  }
426 */
427  wm[i][-iL+j] = a[j];
428  }
429  #endif
430  } /* -- end loop on zones -- */
431 }
432 
433 /* ********************************************************************* */
434 void PPM_CartCoeff(double **wp, double **wm, int beg, int end)
435 /*!
436  * Compute the standard PPM weigth coefficients for a uniformly
437  * spaced Cartesian grid.
438  *
439  *********************************************************************** */
440 {
441  int i;
442 
443  for (i = beg; i <= end; i++){
444  #if PPM_VERSION == PPM3
445  wp[i][-2] = 0.0;
446  wp[i][-1] = -1.0/6.0;
447  wp[i][ 0] = 5.0/6.0;
448  wp[i][ 1] = 1.0/3.0;
449  wp[i][ 2] = 0.0;
450 
451  wm[i][-2] = 0.0;
452  wm[i][-1] = 1.0/3.0;
453  wm[i][ 0] = 5.0/6.0;
454  wm[i][ 1] = -1.0/6.0;
455  wm[i][ 2] = 0.0;
456  #elif PPM_VERSION == PPM4
457  wp[i][-2] = 0.0;
458  wp[i][-1] = -1.0/12.0;
459  wp[i][ 0] = 7.0/12.0;
460  wp[i][ 1] = 7.0/12.0;
461  wp[i][ 2] = -1.0/12.0;
462  #elif PPM_VERSION == PPM5
463  wp[i][-2] = 1.0/30.0;
464  wp[i][-1] = -13.0/60.0;
465  wp[i][ 0] = 47.0/60.0;
466  wp[i][ 1] = 9.0/20.0;
467  wp[i][ 2] = -1.0/20.0;
468 
469  wm[i][-2] = -1.0/20.0;
470  wm[i][-1] = 9.0/20.0;
471  wm[i][ 0] = 47.0/60.0;
472  wm[i][ 1] = -13.0/60.0;
473  wm[i][ 2] = 1.0/30.0;
474  #endif
475  }
476 }
477 
478 double BetaTheta(double x)
479 {
480  return pow(x - gx0, gk)*sin(x);
481 }
482 
483 /* ********************************************************************* */
484 void PPM_Q6_Coeffs(double *hp, double *hm, int dir, Grid *grid)
485 /*
486  *
487  *
488  *********************************************************************** */
489 {
490  int i, beg, end;
491 
492  beg = 0;
493  end = grid[dir].np_tot-1;
494 
495  if (dir == IDIR){
496  double den, *r, *dr;
497 
498  r = grid[IDIR].x;
499  dr = grid[IDIR].dx;
500  for (i = beg; i <= end; i++){
501  #if GEOMETRY == CARTESIAN
502  hp[i] = 3.0;
503  hm[i] = 3.0;
504  #elif GEOMETRY == CYLINDRICAL
505  hp[i] = 3.0 + 0.5*dr[i]/r[i];
506  hm[i] = 3.0 - 0.5*dr[i]/r[i];
507  #elif GEOMETRY == SPHERICAL
508  den = 20.0*r[i]*r[i] + dr[i]*dr[i];
509  hp[i] = 3.0 + 2.0*dr[i]*(10.0*r[i] + dr[i])/den;
510  hm[i] = 3.0 - 2.0*dr[i]*(10.0*r[i] - dr[i])/den;
511  #endif
512  }
513  }else if (dir == JDIR){
514  double *th, *dth, *thp, cp, cm, sp, sm;
515  double dmu, dmu_t;
516 
517  th = grid[JDIR].x;
518  thp = grid[JDIR].xr;
519  dth = grid[JDIR].dx;
520  for (i = beg; i <= end; i++){
521  #if GEOMETRY != SPHERICAL
522  hp[i] = 3.0;
523  hm[i] = 3.0;
524  #else
525  cp = cos(thp[i]); sp = sin(thp[i]);
526  cm = cos(thp[i-1]); sm = sin(thp[i-1]);
527  dmu = cm - cp;
528  dmu_t = sm - sp;
529  hp[i] = dth[i]*(dmu_t + dth[i]*cp)/(dth[i]*(sp + sm) - 2.0*dmu);
530  hm[i] = -dth[i]*(dmu_t + dth[i]*cm)/(dth[i]*(sp + sm) - 2.0*dmu);
531  #endif
532  }
533  }else if (dir == KDIR){
534  for (i = beg; i <= end; i++){
535  hp[i] = 3.0;
536  hm[i] = 3.0;
537  }
538  }
539 }
540 
541 /* ********************************************************************* */
542 void PPM_CoefficientsGet(PPM_Coeffs *ppm_coeffs, int dir)
543 /*!
544  * Retrieve interpolation coefficients for parabolic reconstruction.
545  * This function can be called only if the previous one has been
546  * completed already.
547  *
548  * \param[in] grid pointer to array of grid structure
549  *
550  * \return a pointer to a static structure containing the coefficients
551  * (actually pointers to array of coefficients) in the
552  * direction given by ::g_dir.
553  *********************************************************************** */
554 {
555  if (Wp3D == NULL) {
556  print1 ("! PPM_Coefficients: coefficients not set.\n");
557  QUIT_PLUTO(1);
558  }
559 
560  ppm_coeffs->wp = Wp3D[dir];
561  ppm_coeffs->wm = Wm3D[dir];
562 
563  ppm_coeffs->hp = hp3D[dir];
564  ppm_coeffs->hm = hm3D[dir];
565 
566 }
PPM_Q6_Coeffs
void PPM_Q6_Coeffs(double *hp, double *hm, int dir, Grid *grid)
Definition: ppm_coeffs.v4.1.c:484
a
static double a
Definition: init.c:135
PPM_CoefficientsGet
void PPM_CoefficientsGet(PPM_Coeffs *ppm_coeffs, int dir)
Definition: ppm_coeffs.v4.1.c:542
n
static int n
Definition: analysis.c:3
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
POLY_6
#define POLY_6(a0, a1, a2, a3, a4, a5, a6, x)
Definition: ppm_coeffs.h:13
gk
static int gk
Definition: ppm_coeffs.v4.1.c:256
GRID::xr
double * xr
Definition: structs.h:81
configure.scrh
tuple scrh
Definition: configure.py:200
Wp3D
static double *** Wp3D
Definition: ppm_coeffs.v4.1.c:16
Wm3D
static double *** Wm3D
Definition: ppm_coeffs.v4.1.c:16
ArrayBox
double *** ArrayBox(long int nrl, long int nrh, long int ncl, long int nch, long int ndl, long int ndh)
Definition: arrays.c:341
PPM_COEFFS
Definition: ppm_coeffs.h:29
GRID::dx
double * dx
Definition: structs.h:83
KDIR
#define KDIR
Definition: pluto.h:195
POLY_5
#define POLY_5(a0, a1, a2, a3, a4, a5, x)
Definition: ppm_coeffs.h:11
hp3D
static double ** hp3D
Definition: ppm_coeffs.v4.1.c:16
hm3D
static double ** hm3D
Definition: ppm_coeffs.v4.1.c:16
PPM_CoefficientsSet
void PPM_CoefficientsSet(Grid *grid)
Definition: ppm_coeffs.v4.1.c:23
PPM_COEFFS::hm
double * hm
Definition: ppm_coeffs.h:33
PPM_COEFFS::wp
double ** wp
Definition: ppm_coeffs.h:30
GRID::xl
double * xl
Definition: structs.h:82
LUBackSubst
void LUBackSubst(double **a, int n, int *indx, double b[])
Definition: math_lu_decomp.c:84
IDIR
#define IDIR
Definition: pluto.h:193
GRID
Definition: structs.h:78
GaussQuadrature
double GaussQuadrature(double(*func)(double), double xb, double xe, int nstep, int order)
Definition: math_quadrature.c:20
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
GRID::x
double * x
Definition: structs.h:80
PPM_FindWeights
static void PPM_FindWeights(double **, double **, int, Grid *)
Definition: ppm_coeffs.v4.1.c:258
pluto.h
PLUTO main header file.
gx0
static double gx0
Definition: ppm_coeffs.v4.1.c:255
NMAX_POINT
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
i
int i
Definition: analysis.c:2
PPM_COEFFS::hp
double * hp
Definition: ppm_coeffs.h:32
PPM_COEFFS::wm
double ** wm
Definition: ppm_coeffs.h:31
POLY_3
#define POLY_3(a0, a1, a2, a3, x)
Definition: ppm_coeffs.h:7
POLY_2
#define POLY_2(a0, a1, a2, x)
Definition: ppm_coeffs.h:5
POLY_4
#define POLY_4(a0, a1, a2, a3, a4, x)
Definition: ppm_coeffs.h:9
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
GRID::np_tot
int np_tot
Total number of points in the local domain (boundaries included).
Definition: structs.h:100
LUDecompose
int LUDecompose(double **a, int n, int *indx, double *d)
Definition: math_lu_decomp.c:12
BetaTheta
double BetaTheta(double x)
Definition: ppm_coeffs.v4.1.c:478
JDIR
#define JDIR
Definition: pluto.h:194
POLY_10
#define POLY_10(a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, x)
Definition: ppm_coeffs.h:21
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
DIMENSIONS
#define DIMENSIONS
Definition: definitions_01.h:2
PPM_CartCoeff
static void PPM_CartCoeff(double **, double **, int, int)
Definition: ppm_coeffs.v4.1.c:434
POLY_8
#define POLY_8(a0, a1, a2, a3, a4, a5, a6, a7, a8, x)
Definition: ppm_coeffs.h:18