PLUTO
ppm_coeffs.h File Reference

Go to the source code of this file.

Classes

struct  PPM_COEFFS
 

Macros

#define PPM_ORDER   4
 
#define POLY_2(a0, a1, a2, x)    ( a0 + x*(a1 + x*a2) )
 
#define POLY_3(a0, a1, a2, a3, x)   ( a0 + x*(a1 + x*(a2 + x*a3)) )
 
#define POLY_4(a0, a1, a2, a3, a4, x)   ( a0 + x*(a1 + x*(a2 + x*(a3 + x*a4))) )
 
#define POLY_5(a0, a1, a2, a3, a4, a5, x)   ( a0 + x*(a1 + x*(a2 + x*(a3 + x*(a4 + x*a5)))) )
 
#define POLY_6(a0, a1, a2, a3, a4, a5, a6, x)   ( a0 + x*(a1 + x*(a2 + x*(a3 + x*(a4 + x*(a5 + x*a6))))) )
 
#define POLY_7(a0, a1, a2, a3, a4, a5, a6, a7, x)   ( a0 + x*(a1 + x*(a2 + x*(a3 + x*(a4 + x*(a5 + x*(a6 + x*a7)))))) )
 
#define POLY_8(a0, a1, a2, a3, a4, a5, a6, a7, a8, x)   ( a0 + x*(a1 + x*(a2 + x*(a3 + x*POLY_4(a4, a5, a6, a7, a8, x)))) )
 
#define POLY_10(a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, x)   ( a0 + x*(a1 + x*(a2 + x*(a3 + x*(a4 + x*POLY_5(a5, a6, a7, a8, a9, a10, x))))) )
 

Typedefs

typedef struct PPM_COEFFS PPM_Coeffs
 

Functions

void PPM_CoefficientsSet (Grid *grid)
 
void PPM_CoefficientsGet (PPM_Coeffs *, int)
 
PPM_CoeffsPPM_Coefficients (int action, Grid *grid)
 
void PPM_Q6_Coeffs (double *hp, double *hm, int dir, Grid *grid)
 

Macro Definition Documentation

#define POLY_10 (   a0,
  a1,
  a2,
  a3,
  a4,
  a5,
  a6,
  a7,
  a8,
  a9,
  a10,
 
)    ( a0 + x*(a1 + x*(a2 + x*(a3 + x*(a4 + x*POLY_5(a5, a6, a7, a8, a9, a10, x))))) )

Definition at line 21 of file ppm_coeffs.h.

#define POLY_2 (   a0,
  a1,
  a2,
 
)    ( a0 + x*(a1 + x*a2) )

Definition at line 5 of file ppm_coeffs.h.

#define POLY_3 (   a0,
  a1,
  a2,
  a3,
 
)    ( a0 + x*(a1 + x*(a2 + x*a3)) )

Definition at line 7 of file ppm_coeffs.h.

#define POLY_4 (   a0,
  a1,
  a2,
  a3,
  a4,
 
)    ( a0 + x*(a1 + x*(a2 + x*(a3 + x*a4))) )

Definition at line 9 of file ppm_coeffs.h.

#define POLY_5 (   a0,
  a1,
  a2,
  a3,
  a4,
  a5,
 
)    ( a0 + x*(a1 + x*(a2 + x*(a3 + x*(a4 + x*a5)))) )

Definition at line 11 of file ppm_coeffs.h.

#define POLY_6 (   a0,
  a1,
  a2,
  a3,
  a4,
  a5,
  a6,
 
)    ( a0 + x*(a1 + x*(a2 + x*(a3 + x*(a4 + x*(a5 + x*a6))))) )

Definition at line 13 of file ppm_coeffs.h.

#define POLY_7 (   a0,
  a1,
  a2,
  a3,
  a4,
  a5,
  a6,
  a7,
 
)    ( a0 + x*(a1 + x*(a2 + x*(a3 + x*(a4 + x*(a5 + x*(a6 + x*a7)))))) )

Definition at line 15 of file ppm_coeffs.h.

#define POLY_8 (   a0,
  a1,
  a2,
  a3,
  a4,
  a5,
  a6,
  a7,
  a8,
 
)    ( a0 + x*(a1 + x*(a2 + x*(a3 + x*POLY_4(a4, a5, a6, a7, a8, x)))) )

Definition at line 18 of file ppm_coeffs.h.

#define PPM_ORDER   4

Definition at line 2 of file ppm_coeffs.h.

Typedef Documentation

typedef struct PPM_COEFFS PPM_Coeffs

Simple structure used to retrieve 1D reconstruction weights (c, w, d) used by piecewise linear interpolation (see states_plm.c)

Function Documentation

PPM_Coeffs* PPM_Coefficients ( int  action,
Grid grid 
)
void PPM_CoefficientsGet ( PPM_Coeffs ppm_coeffs,
int  dir 
)

Retrieve interpolation coefficients for parabolic reconstruction. This function can be called only if the previous one has been completed already.

Parameters
[in]gridpointer to array of grid structure
Returns
a pointer to a static structure containing the coefficients (actually pointers to array of coefficients) in the direction given by g_dir.

Definition at line 591 of file ppm_coeffs.c.

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 }
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
static double *** s_Wm3D
Definition: ppm_coeffs.c:49
double * hm
Definition: ppm_coeffs.h:33
double ** wp
Definition: ppm_coeffs.h:30
static double ** s_hp3D
Definition: ppm_coeffs.c:49
static double *** s_Wp3D
Definition: ppm_coeffs.c:49
double * hp
Definition: ppm_coeffs.h:32
double ** wm
Definition: ppm_coeffs.h:31
static double ** s_hm3D
Definition: ppm_coeffs.c:49
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the call graph for this function:

Here is the caller graph for this function:

void PPM_CoefficientsSet ( Grid grid)

Compute interpolation coefficients for PPM in every coordinate system.

Definition at line 58 of file ppm_coeffs.c.

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);
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 }
int end
Global end index for the local array.
Definition: structs.h:116
static int n
Definition: analysis.c:3
#define POLY_6(a0, a1, a2, a3, a4, a5, a6, x)
Definition: ppm_coeffs.h:13
double * xr
Definition: structs.h:81
double *** ArrayBox(long int nrl, long int nrh, long int ncl, long int nch, long int ndl, long int ndh)
Definition: arrays.c:341
double * dx
Definition: structs.h:83
#define KDIR
Definition: pluto.h:195
#define POLY_5(a0, a1, a2, a3, a4, a5, x)
Definition: ppm_coeffs.h:11
static double *** s_Wm3D
Definition: ppm_coeffs.c:49
static void PPM_FindWeights(double **, double **, int, Grid *)
Definition: ppm_coeffs.c:327
#define IDIR
Definition: pluto.h:193
int beg
Global start index for the local array.
Definition: structs.h:115
void PPM_Q6_Coeffs(double *hp, double *hm, int dir, Grid *grid)
Definition: ppm_coeffs.c:533
static double ** s_hp3D
Definition: ppm_coeffs.c:49
static void PPM_CartCoeff(double **, double **, int, int)
Definition: ppm_coeffs.c:483
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
static double *** s_Wp3D
Definition: ppm_coeffs.c:49
int i
Definition: analysis.c:2
#define POLY_3(a0, a1, a2, a3, x)
Definition: ppm_coeffs.h:7
#define POLY_2(a0, a1, a2, x)
Definition: ppm_coeffs.h:5
#define POLY_4(a0, a1, a2, a3, a4, x)
Definition: ppm_coeffs.h:9
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
int np_tot
Total number of points in the local domain (boundaries included).
Definition: structs.h:100
static double ** s_hm3D
Definition: ppm_coeffs.c:49
#define JDIR
Definition: pluto.h:194
#define POLY_10(a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, x)
Definition: ppm_coeffs.h:21
#define DIMENSIONS
Definition: definitions_01.h:2
#define POLY_8(a0, a1, a2, a3, a4, a5, a6, a7, a8, x)
Definition: ppm_coeffs.h:18

Here is the call graph for this function:

Here is the caller graph for this function:

void PPM_Q6_Coeffs ( double *  hp,
double *  hm,
int  dir,
Grid grid 
)

Definition at line 533 of file ppm_coeffs.c.

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 }
int end
Global end index for the local array.
Definition: structs.h:116
double * xr
Definition: structs.h:81
double * dx
Definition: structs.h:83
#define KDIR
Definition: pluto.h:195
#define IDIR
Definition: pluto.h:193
int beg
Global start index for the local array.
Definition: structs.h:115
double * x
Definition: structs.h:80
int i
Definition: analysis.c:2
int np_tot
Total number of points in the local domain (boundaries included).
Definition: structs.h:100
#define JDIR
Definition: pluto.h:194

Here is the caller graph for this function: