PLUTO
Macros | Functions | Variables
ppm_coeffs.c File Reference

Compute coefficients for high-order reconstruction methods. More...

#include "pluto.h"
Include dependency graph for ppm_coeffs.c:

Go to the source code of this file.

Macros

#define PPM_ORDER   3
 

Functions

static void PPM_CartCoeff (double **, double **, int, int)
 
static void PPM_FindWeights (double **, double **, int, Grid *)
 
static double BetaTheta (double)
 
void PPM_CoefficientsSet (Grid *grid)
 
void PPM_Q6_Coeffs (double *hp, double *hm, int dir, Grid *grid)
 
void PPM_CoefficientsGet (PPM_Coeffs *ppm_coeffs, int dir)
 

Variables

static double *** s_Wp3D
 
static double *** s_Wm3D
 
static double ** s_hp3D
 
static double ** s_hm3D
 
static double s_x0
 
static int s_k
 

Detailed Description

Compute coefficients for high-order reconstruction methods.

Compute the interpolation coefficients needed by high-order (3rd, 4th and 5th) reconstruction methods such as PPM or WENO3. The function PPM_CoefficientsSet() must be called to initialize arrays and compute the coefficients after the grid has been generated. The function PPM_CoefficientsGet() can be used at anytime to retrieve the coefficients in a particular direction.

Reconstruction coefficients are computed in different ways depending on the geometry and grid uniformity:

  • Uniform Cartesian grids: reconstruction coefficients are computed by PPM_CartCoeff().
  • Uniform curvilinear grids: coefficients are computed in PPM_CoefficientsSet() for radial grids (polar/cyindrical and spherical), by PPM_FindWeights() for meridional spherical coordinate.
  • Non uniform grids: reconstruction coeffients are always computed by inverting the Vandermonde-like equation (see Eq. [21] in Mignone, JCP 2014) in PPM_FindWeights()
Authors
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
Dec 29, 2014

Reference

  • "High-order conservative reconstruction schemes for finite volume methods in cylindrical and spherical coordinates", A. Mignone, JCP (2014), 270, 784.

Definition in file ppm_coeffs.c.

Macro Definition Documentation

#define PPM_ORDER   3

Definition at line 42 of file ppm_coeffs.c.

Function Documentation

double BetaTheta ( double  x)
static

Definition at line 527 of file ppm_coeffs.c.

528 {
529  return pow(x - s_x0, s_k)*sin(x);
530 }
GRID::x
double * x
Definition: structs.h:80
s_x0
static double s_x0
Definition: ppm_coeffs.c:50
s_k
static int s_k
Definition: ppm_coeffs.c:51

Here is the caller graph for this function:

void PPM_CartCoeff ( double **  wp,
double **  wm,
int  beg,
int  end 
)
static

Compute the standard PPM weigth coefficients for a uniformly spaced Cartesian grid.

Definition at line 483 of file ppm_coeffs.c.

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 }
GRID::end
int end
Global end index for the local array.
Definition: structs.h:116
GRID::beg
int beg
Global start index for the local array.
Definition: structs.h:115
i
int i
Definition: analysis.c:2

Here is the caller graph for this function:

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

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);
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 }
GRID::end
int end
Global end index for the local array.
Definition: structs.h:116
n
static int n
Definition: analysis.c:3
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
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
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_FindWeights
static void PPM_FindWeights(double **, double **, int, Grid *)
Definition: ppm_coeffs.c:327
IDIR
#define IDIR
Definition: pluto.h:193
GRID::beg
int beg
Global start index for the local array.
Definition: structs.h:115
PPM_Q6_Coeffs
void PPM_Q6_Coeffs(double *hp, double *hm, int dir, Grid *grid)
Definition: ppm_coeffs.c:533
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
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
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
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
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

Here is the caller graph for this function:

void PPM_FindWeights ( double **  wp,
double **  wm,
int  dir,
Grid *  grid 
)
static

Find the coefficients numerically by inverting the beta matrix We define the matrix beta[0...n][0...j] so we need to index it as beta[k][j-jb]

Definition at line 327 of file ppm_coeffs.c.

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 }
a
static double a
Definition: init.c:135
GRID::end
int end
Global end index for the local array.
Definition: structs.h:116
n
static int n
Definition: analysis.c:3
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
GRID::xr
double * xr
Definition: structs.h:81
configure.scrh
tuple scrh
Definition: configure.py:200
BetaTheta
static double BetaTheta(double)
Definition: ppm_coeffs.c:527
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::beg
int beg
Global start index for the local array.
Definition: structs.h:115
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
s_x0
static double s_x0
Definition: ppm_coeffs.c:50
i
int i
Definition: analysis.c:2
s_k
static int s_k
Definition: ppm_coeffs.c:51
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
JDIR
#define JDIR
Definition: pluto.h:194

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

Here is the caller graph for this function:

Variable Documentation

double ** s_hm3D
static

Definition at line 49 of file ppm_coeffs.c.

double ** s_hp3D
static

Definition at line 49 of file ppm_coeffs.c.

int s_k
static

Definition at line 51 of file ppm_coeffs.c.

double *** s_Wm3D
static

Definition at line 49 of file ppm_coeffs.c.

double*** s_Wp3D
static

Definition at line 49 of file ppm_coeffs.c.

double s_x0
static

Definition at line 50 of file ppm_coeffs.c.