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