PLUTO
Functions
set_geometry.c File Reference

Initialize geometry-dependent grid quantities. More...

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

Go to the source code of this file.

Functions

void MakeGeometry (Grid *GXYZ)
 
double Length_1 (int i, int j, int k, Grid *grid)
 
double Length_2 (int i, int j, int k, Grid *grid)
 
double Length_3 (int i, int j, int k, Grid *grid)
 
double * GetInverse_dl (const Grid *grid)
 

Detailed Description

Initialize geometry-dependent grid quantities.

Compute grid quantities (such as interface areas, volumes, centroid of volume, etc..) that depend on the geometry.

Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
Dec 12, 2013

Definition in file set_geometry.c.

Function Documentation

double* GetInverse_dl ( const Grid *  grid)

Definition at line 205 of file set_geometry.c.

218 {
219 #if GEOMETRY == CARTESIAN || GEOMETRY == CYLINDRICAL
220 
221  return grid[g_dir].inv_dx;
222 
223 #elif GEOMETRY == POLAR
224 
225  if (g_dir != JDIR){
226  return grid[g_dir].inv_dx;
227  }else{
228  int j;
229  double r_1;
230  static double *inv_dl;
231 
232  if (inv_dl == NULL) {
233  #ifdef CHOMBO
234  inv_dl = ARRAY_1D(NX2_MAX, double);
235  #else
236  inv_dl = ARRAY_1D(NX2_TOT, double);
237  #endif
238  }
239  r_1 = grid[IDIR].r_1[g_i];
240  JTOT_LOOP(j) inv_dl[j] = grid[JDIR].inv_dx[j]*r_1;
241  return inv_dl;
242  }
243 
244 #elif GEOMETRY == SPHERICAL
245 
246  int j, k;
247  double r_1, s;
248  static double *inv_dl2, *inv_dl3;
249 
250  if (inv_dl2 == NULL) {
251  #ifdef CHOMBO
252  inv_dl2 = ARRAY_1D(NX2_MAX, double);
253  inv_dl3 = ARRAY_1D(NX3_MAX, double);
254  #else
255  inv_dl2 = ARRAY_1D(NX2_TOT, double);
256  inv_dl3 = ARRAY_1D(NX3_TOT, double);
257  #endif
258  }
259 
260  if (g_dir == IDIR){
261  return grid[IDIR].inv_dx;
262  }else if (g_dir == JDIR) {
263  r_1 = grid[IDIR].r_1[g_i];
264  JTOT_LOOP(j) inv_dl2[j] = grid[JDIR].inv_dx[j]*r_1;
265  return inv_dl2;
266  }else if (g_dir == KDIR){
267  r_1 = grid[IDIR].r_1[g_i];
268  s = grid[JDIR].x[g_j];
269  s = sin(s);
270  KTOT_LOOP(k) inv_dl3[k] = grid[KDIR].inv_dx[k]*r_1/s;
271  return inv_dl3;
272  }
273 
274 #endif
275 
276  return NULL;
277 }
NX2_TOT
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
Definition: globals.h:57
KTOT_LOOP
#define KTOT_LOOP(k)
Definition: macros.h:40
KDIR
#define KDIR
Definition: pluto.h:195
GRID::inv_dx
double * inv_dx
Definition: structs.h:90
NX3_MAX
#define NX3_MAX
Definition: pluto.h:743
g_i
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
Definition: globals.h:82
NX2_MAX
#define NX2_MAX
Definition: pluto.h:742
IDIR
#define IDIR
Definition: pluto.h:193
g_dir
int g_dir
Specifies the current sweep or direction of integration.
Definition: globals.h:86
g_j
int g_j
x2 grid index when sweeping along the x1 or x3 direction.
Definition: globals.h:83
NX3_TOT
long int NX3_TOT
Total number of zones in the X3 direction (boundaries included) for the local processor.
Definition: globals.h:59
j
int j
Definition: analysis.c:2
if
if(divB==NULL)
Definition: analysis.c:10
k
int k
Definition: analysis.c:2
GRID::x
double * x
Definition: structs.h:80
GRID::r_1
double * r_1
Geometrical factor 1/r.
Definition: structs.h:88
s
#define s
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
JTOT_LOOP
#define JTOT_LOOP(j)
Definition: macros.h:39
JDIR
#define JDIR
Definition: pluto.h:194

Here is the caller graph for this function:

double Length_1 ( int  i,
int  j,
int  k,
Grid *  grid 
)

Definition at line 156 of file set_geometry.c.

164 {
165  return (grid[0].dx[i]);
166 }
GRID::dx
double * dx
Definition: structs.h:83
i
int i
Definition: analysis.c:2

Here is the caller graph for this function:

double Length_2 ( int  i,
int  j,
int  k,
Grid *  grid 
)

Definition at line 169 of file set_geometry.c.

177 {
178  #if GEOMETRY == CARTESIAN || GEOMETRY == CYLINDRICAL
179  return (grid[1].dx[j]);
180  #elif GEOMETRY == POLAR || GEOMETRY == SPHERICAL
181  return (fabs(grid[0].xgc[i])*grid[1].dx[j]);
182  #endif
183 }
GRID::dx
double * dx
Definition: structs.h:83
j
int j
Definition: analysis.c:2
GRID::xgc
double * xgc
Cell volumetric centroid (!= x when geometry != CARTESIAN).
Definition: structs.h:84
i
int i
Definition: analysis.c:2

Here is the caller graph for this function:

double Length_3 ( int  i,
int  j,
int  k,
Grid *  grid 
)

Definition at line 186 of file set_geometry.c.

194 {
195  #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
196  return (grid[2].dx[k]);
197  #elif GEOMETRY == CYLINDRICAL
198  return (fabs(grid[0].xgc[i])*grid[2].dx[k]);
199  #elif GEOMETRY == SPHERICAL
200  return (fabs(grid[0].xgc[i]*sin(grid[1].xgc[j]))*grid[2].dx[k]);
201  #endif
202 }
GRID::dx
double * dx
Definition: structs.h:83
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
GRID::xgc
double * xgc
Cell volumetric centroid (!= x when geometry != CARTESIAN).
Definition: structs.h:84
i
int i
Definition: analysis.c:2

Here is the caller graph for this function:

void MakeGeometry ( Grid *  GXYZ)
Parameters
[in,out]GXYZPointer to an array of Grid structures;

Definition at line 16 of file set_geometry.c.

22 {
23  int i, j, k, idim, ngh, ileft;
24  int iright;
25  int iend, jend, kend;
26  double xiL, xiR, dxi, dvol;
27  double x, dx, xr, xl;
28  double y, dy, yr, yl;
29  struct GRID *GG;
30 
31  iend = GXYZ[0].lend + GXYZ[0].nghost;
32  jend = GXYZ[1].lend + GXYZ[1].nghost;
33  kend = GXYZ[2].lend + GXYZ[2].nghost;
34 
35 /* --------------------------------------------------------------
36  Memory allocation. All values are defined at the cell center
37  with the exception of the area element which is defined on a
38  staggered mesh and therefore starts at [-1].
39  ----------------------------------------------------------- */
40 
41  for (idim = 0; idim < 3; idim++) {
42  (GXYZ + idim)->A = ARRAY_1D (GXYZ[idim].np_tot+1, double)+1;
43  (GXYZ + idim)->xgc = ARRAY_1D (GXYZ[idim].np_tot, double);
44  (GXYZ + idim)->dV = ARRAY_1D (GXYZ[idim].np_tot, double);
45  (GXYZ + idim)->r_1 = ARRAY_1D (GXYZ[idim].np_tot, double);
46  (GXYZ + idim)->ct = ARRAY_1D (GXYZ[idim].np_tot, double);
47  (GXYZ + idim)->inv_dx = ARRAY_1D (GXYZ[idim].np_tot, double);
48  (GXYZ + idim)->inv_dxi = ARRAY_1D (GXYZ[idim].np_tot, double);
49  }
50 
51 /* ------------------------------------------------------------
52  Define area (A), volume element (dV) and cell geometrical
53  centers for each direction.
54 
55  Conventions:
56  - G->dx: spacing (always > 0)
57  - G->x: cell-center. Can be > 0 or < 0
58  - G->xgc: geometrical cell center.
59  ----------------------------------------------------------- */
60 
61 /* ------------------------------------------------------------
62  X1 (IDIR) Direction
63  ------------------------------------------------------------ */
64 
65  GG = GXYZ;
66  for (i = 0; i <= iend; i++) {
67 
68  dx = GG->dx[i];
69  x = GG->x[i];
70  xr = x + 0.5*dx;
71  xl = x - 0.5*dx;
72 
73  #if GEOMETRY == CARTESIAN
74  GG->A[i] = 1.0;
75  if (i == 0) GG->A[-1] = 1.0;
76  GG->dV[i] = dx;
77  GG->xgc[i] = x;
78  #elif GEOMETRY == CYLINDRICAL || GEOMETRY == POLAR
79  GG->A[i] = fabs(xr);
80  if (i == 0) GG->A[-1] = fabs(xl);
81  GG->dV[i] = fabs(x)*dx;
82  GG->xgc[i] = x + dx*dx/(12.0*x);
83  GG->r_1[i] = 1.0/x;
84  #elif GEOMETRY == SPHERICAL
85  GG->A[i] = xr*xr;
86  if (i == 0) GG->A[-1] = xl*xl;
87 
88  GG->dV[i] = fabs(xr*xr*xr - xl*xl*xl)/3.0;
89  GG->xgc[i] = x + 2.0*x*dx*dx/(12.0*x*x + dx*dx);
90  GG->r_1[i] = 1.0/x;
91  #endif
92  }
93 
94 /* ------------------------------------------------------------
95  X2 (JDIR) Direction
96  ------------------------------------------------------------ */
97 
98  GG = GXYZ + 1;
99  for (j = 0; j <= jend; j++) {
100 
101  dx = GG->dx[j];
102  x = GG->x[j];
103  xr = x + 0.5*dx;
104  xl = x - 0.5*dx;
105  #if GEOMETRY != SPHERICAL
106  GG->A[j] = 1.0;
107  if (j == 0) GG->A[-1] = 1.0;
108  GG->dV[j] = dx;
109  GG->xgc[j] = x;
110 
111  #else
112  GG->A[j] = fabs(sin(xr));
113  if (j == 0) GG->A[-1] = fabs(sin(xl));
114  GG->dV[j] = fabs(cos(xl) - cos(xr));
115 
116  GG->xgc[j] = (sin(xr) - sin(xl) + xl*cos(xl) - xr*cos(xr));
117  GG->xgc[j] /= cos(xl) - cos(xr);
118  GG->ct[j] = 1.0/tan(x); /* (sin(xr) - sin(xl))/(cos(xl) - cos(xr)); */
119  #endif
120  }
121 
122 /* ------------------------------------------------------------
123  X3 (KDIR) Direction
124  ------------------------------------------------------------ */
125 
126  GG = GXYZ + 2;
127  for (k = 0; k <= kend; k++) {
128  dx = GG->dx[k];
129  x = GG->x[k];
130  xr = x + 0.5*dx;
131  xl = x - 0.5*dx;
132 
133  GG->A[k] = 1.0;
134  if (k == 0) GG->A[-1] = 1.0;
135  GG->dV[k] = dx;
136  GG->xgc[k] = x;
137  }
138 
139 /* ---------------------------------------------------------
140  compute and store the reciprocal of cell spacing
141  between interfaces (inv_dx) and cell centers (inv_dxi)
142  --------------------------------------------------------- */
143 
144  for (idim = 0; idim < DIMENSIONS; idim++){
145  for (i = 0; i < GXYZ[idim].np_tot; i++) {
146  GXYZ[idim].inv_dx[i] = 1.0/(GXYZ[idim].dx[i]);
147  }
148 
149  for (i = 0; i < GXYZ[idim].np_tot-1; i++) {
150  GXYZ[idim].inv_dxi[i] = 2.0/(GXYZ[idim].dx[i] + GXYZ[idim].dx[i+1]);
151  }
152  }
153 }
GRID::xr
double * xr
Definition: structs.h:81
GRID::lend
int lend
Local end index for the local array.
Definition: structs.h:118
GRID::dV
double * dV
Cell volume.
Definition: structs.h:86
GRID::dx
double * dx
Definition: structs.h:83
GRID::inv_dx
double * inv_dx
Definition: structs.h:90
GRID::xl
double * xl
Definition: structs.h:82
GRID
Definition: structs.h:78
j
int j
Definition: analysis.c:2
GRID::nghost
int nghost
Number of ghost zones.
Definition: structs.h:104
k
int k
Definition: analysis.c:2
GRID::x
double * x
Definition: structs.h:80
GRID::xgc
double * xgc
Cell volumetric centroid (!= x when geometry != CARTESIAN).
Definition: structs.h:84
GRID::r_1
double * r_1
Geometrical factor 1/r.
Definition: structs.h:88
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
i
int i
Definition: analysis.c:2
GRID::ct
double * ct
Geometrical factor cot(theta).
Definition: structs.h:89
GRID::np_tot
int np_tot
Total number of points in the local domain (boundaries included).
Definition: structs.h:100
GRID::inv_dxi
double * inv_dxi
inverse of the distance between the center of two cells, inv_dxi = .
Definition: structs.h:91
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87
DIMENSIONS
#define DIMENSIONS
Definition: definitions_01.h:2

Here is the caller graph for this function: