PLUTO
Functions
init.c File Reference
#include "pluto.h"
Include dependency graph for init.c:

Go to the source code of this file.

Functions

void Init (double *us, double x1, double x2, double x3)
 
void Analysis (const Data *d, Grid *grid)
 
void UserDefBoundary (const Data *d, RBox *box, int side, Grid *grid)
 

Function Documentation

void Analysis ( const Data *  d,
Grid *  grid 
)

Perform runtime data analysis.

Parameters
[in]dthe PLUTO Data structure
[in]gridpointer to array of Grid structures

Definition at line 61 of file init.c.

66 {
67  int i,j,k,nv;
68  static int first_call = 1;
69  double v[NVAR];
70  double *R, *phi, *dVR, *dVphi, *dR, *dphi, vol;
71  double **vR, **vphi;
72  double E=0.0, w, pw, wmin, pwmin;
73  double Ltot, Etot, wtot, pwtot;
74  static double voltot;
75  FILE *fp;
76 
77  R = grid[IDIR].x; phi = grid[JDIR].x;
78  dVR = grid[IDIR].dV; dVphi = grid[JDIR].dV;
79  dR = grid[IDIR].dx; dphi = grid[JDIR].dx;
80 
81 /* ---------------------------------------------------------
82  compute total volume once at the beginning
83  --------------------------------------------------------- */
84 
85  if (first_call){
86  voltot = 0.0;
87  DOM_LOOP(k,j,i) voltot += dVR[i]*dVphi[j];
88  #ifdef PARALLEL
89  MPI_Allreduce (&voltot, &vol, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
90  MPI_Barrier (MPI_COMM_WORLD);
91  voltot = vol;
92  #endif
93  }
94 
95 /* --------------------------------------------------------
96  compute volume averages
97  -------------------------------------------------------- */
98 
99  vR = d->Vc[VX1][0]; vphi = d->Vc[VX2][0];
100 
101  Etot = Ltot = wtot = pwtot = 0.0;
102  wmin = pwmin = 1.e20;
103  DOM_LOOP(k,j,i){
104  vol = dVR[i]*dVphi[j]/voltot;
105 
106  for (nv = NVAR; nv--; ) v[nv] = d->Vc[nv][k][j][i];
107 
108  w = ( (vphi[j][i+1] - 1.0/sqrt(R[i+1]))*R[i+1]
109  - (vphi[j][i-1] - 1.0/sqrt(R[i-1]))*R[i-1])/(2.0*R[i]*dR[i]);
110  w -= (vR[j+1][i] - vR[j-1][i])/(2.0*R[i]*dphi[j]);
111 
112  pw = w/v[RHO];
113  #if EOS == IDEAL
114  E = v[PRS]/(g_gamma-1.0)
115  + 0.5*v[RHO]*(v[VX1]*v[VX1] + v[VX2]*v[VX2]) - v[RHO]/R[i];
116  #if PHYSICS == MHD
117  E += 0.5*(v[BX1]*v[BX1] + v[BX2]*v[BX2]);
118  #endif
119  #endif
120 
121  Ltot += R[i]*v[RHO]*v[VX2]*vol;
122  Etot += E*vol;
123  wtot += w*vol;
124  pwtot += pw*vol;
125 
126  wmin = MIN(wmin, w);
127  pwmin = MIN(pwmin, pw);
128  }
129 
130  #ifdef PARALLEL
131  MPI_Allreduce (&Ltot, &vol, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
132  Ltot = vol;
133 
134  MPI_Allreduce (&Etot, &vol, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
135  Etot = vol;
136 
137  MPI_Allreduce (&wtot, &vol, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
138  wtot = vol;
139 
140  MPI_Allreduce (&pwtot, &vol, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
141  pwtot = vol;
142 
143  MPI_Allreduce (&wmin, &vol, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
144  wmin = vol;
145 
146  MPI_Allreduce (&pwmin, &vol, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
147  pwmin = vol;
148 
149  MPI_Barrier (MPI_COMM_WORLD);
150  #endif
151 /*
152 printf ("<pm> = %12.6e, <pr> = %12.6e, <pr>/<pm> = %12.6e, <pr*f>/<pm*f> = %12.6e\n",
153  pmtot, prtot, prtot/pmtot, beta_num/beta_den);
154 exit(1);
155 */
156 /* -----------------------------------------------------
157  processor #0 does the actual writing
158  ----------------------------------------------------- */
159 
160  if (prank == 0){
161  static double tpos = -1.0;
162  if (g_stepNumber == 0){
163  fp = fopen("averages.dat","w");
164  fprintf (fp,"#%s %12s %13s %13s %13s %14s %12s %12s\n",
165  " t","dt","<L>","<E>","<w>","<w/rho>","min(w)","min(pw)");
166  }else{
167  if (tpos < 0.0){ /* obtain time coordinate of to last entry */
168  char sline[512];
169  fp = fopen("averages.dat","r");
170  while (fgets(sline, 512, fp)) {}
171  sscanf(sline, "%lf\n",&tpos);
172  fclose(fp);
173  }
174  fp = fopen("averages.dat","a");
175  }
176  if (g_time > tpos){
177  fprintf (fp, "%12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n",
178  g_time, g_dt, Ltot, Etot, wtot, pwtot, wmin, pwmin);
179  }
180  fclose(fp);
181  }
182  first_call = 0;
183 }
g_gamma
double g_gamma
Definition: globals.h:112
DOM_LOOP
DOM_LOOP(k, j, i)
Definition: analysis.c:22
VX2
#define VX2
Definition: mod_defs.h:29
RHO
#define RHO
Definition: mod_defs.h:19
DATA::Vc
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
g_dt
double g_dt
The current integration time step.
Definition: globals.h:118
GRID::dV
double * dV
Cell volume.
Definition: structs.h:86
GRID::dx
double * dx
Definition: structs.h:83
prank
int prank
Processor rank.
Definition: globals.h:33
VX1
#define VX1
Definition: mod_defs.h:28
MIN
#define MIN(a, b)
Definition: macros.h:104
g_stepNumber
long int g_stepNumber
Gives the current integration step number.
Definition: globals.h:97
IDIR
#define IDIR
Definition: pluto.h:193
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
GRID::x
double * x
Definition: structs.h:80
i
int i
Definition: analysis.c:2
g_time
double g_time
The current integration time.
Definition: globals.h:117
fp
FILE * fp
Definition: analysis.c:7
BX1
#define BX1
Definition: mod_defs.h:25
BX2
#define BX2
Definition: mod_defs.h:26
JDIR
#define JDIR
Definition: pluto.h:194
NVAR
#define NVAR
Definition: pluto.h:609

Here is the call graph for this function:

void Init ( double *  us,
double  x1,
double  x2,
double  x3 
)

The Init() function can be used to assign initial conditions as as a function of spatial position.

Parameters
[out]va pointer to a vector of primitive variables
[in]x1coordinate point in the 1st dimension
[in]x2coordinate point in the 2nd dimension
[in]x3coordinate point in the 3rdt dimension

The meaning of x1, x2 and x3 depends on the geometry:

\[ \begin{array}{cccl} x_1 & x_2 & x_3 & \mathrm{Geometry} \\ \noalign{\medskip} \hline x & y & z & \mathrm{Cartesian} \\ \noalign{\medskip} R & z & - & \mathrm{cylindrical} \\ \noalign{\medskip} R & \phi & z & \mathrm{polar} \\ \noalign{\medskip} r & \theta & \phi & \mathrm{spherical} \end{array} \]

Variable names are accessed by means of an index v[nv], where nv = RHO is density, nv = PRS is pressure, nv = (VX1, VX2, VX3) are the three components of velocity, and so forth.

Definition at line 4 of file init.c.

10 {
11  static int first_call = 1;
12  double x, y, R, Rc, phi, c, s;
13  double arg, kp, h;
14  double dvx, dvy;
15 
16  kp = g_inputParam[KPAR];
17  h = g_inputParam[HSCALE]; /* -- the scale of the vortex -- */
18  Rc = 1.0; /* -- the radial position of the vortex -- */
19 
20 /* --------------------------------------------
21  set coordinates
22  -------------------------------------------- */
23 
24  R = x1;
25  phi = x2;
26  c = cos(phi);
27  s = sin(phi);
28 
29  x = R*c - Rc/sqrt(2.0);
30  y = R*s - Rc/sqrt(2.0);
31 
32  us[RHO] = 1.0;
33  #if EOS == IDEAL
34  us[PRS] = 1.0/(g_inputParam[MACH]*g_inputParam[MACH]*g_gamma);
35  #else
36  g_isoSoundSpeed = 1.0/g_inputParam[MACH];
37  #endif
38  us[VX1] = 0.0;
39  us[VX2] = 1.0/sqrt(R);
40 
41  arg = - (x*x + y*y)/(h*h);
42 
43  dvx = -y*kp*exp(arg);
44  dvy = x*kp*exp(arg);
45 
46  us[VX1] += dvx*c + dvy*s;
47  us[VX2] += -dvx*s + dvy*c;
48 
49  #if ROTATING_FRAME
50  g_OmegaZ = 1.0;
51  us[VX2] -= g_OmegaZ*R;
52  #endif
53 
54  #if NTRACER == 1
55  us[TRC] = (exp(arg) > 1.e-2);
56  #endif
57 
58 }
g_gamma
double g_gamma
Definition: globals.h:112
MACH
#define MACH
Definition: definitions_01.h:25
VX2
#define VX2
Definition: mod_defs.h:29
RHO
#define RHO
Definition: mod_defs.h:19
g_OmegaZ
#define g_OmegaZ
Definition: init.c:64
TRC
#define TRC
Definition: pluto.h:581
VX1
#define VX1
Definition: mod_defs.h:28
KPAR
#define KPAR
Definition: definitions_01.h:23
g_inputParam
double g_inputParam[32]
Array containing the user-defined parameters.
Definition: globals.h:131
GRID::x
double * x
Definition: structs.h:80
menu.c
tuple c
Definition: menu.py:375
s
#define s
HSCALE
#define HSCALE
Definition: definitions_01.h:24
void UserDefBoundary ( const Data *  d,
RBox *  box,
int  side,
Grid *  grid 
)

Assign user-defined boundary conditions.

Parameters
[in,out]dpointer to the PLUTO data structure containing cell-centered primitive quantities (d->Vc) and staggered magnetic fields (d->Vs, when used) to be filled.
[in]boxpointer to a RBox structure containing the lower and upper indices of the ghost zone-centers/nodes or edges at which data values should be assigned.
[in]sidespecifies the boundary side where ghost zones need to be filled. It can assume the following pre-definite values: X1_BEG, X1_END, X2_BEG, X2_END, X3_BEG, X3_END. The special value side == 0 is used to control a region inside the computational domain.
[in]gridpointer to an array of Grid structures.

Definition at line 185 of file init.c.

206 {
207  int i, j, k, nv;
208  double *R, v[256];
209 
210  R = grid[IDIR].x;
211 
212  if (side == X1_BEG){
213 
214  if (box->vpos == CENTER){
215  BOX_LOOP(box,k,j,i){
216  d->Vc[RHO][k][j][i] = d->Vc[RHO][k][j][IBEG];
217  #if EOS == IDEAL
218  d->Vc[PRS][k][j][i] = d->Vc[PRS][k][j][IBEG];
219  #endif
220  d->Vc[VX1][k][j][i] = d->Vc[VX1][k][j][IBEG];
221  d->Vc[VX2][k][j][i] = 1.0/sqrt(R[i]) + d->Vc[VX2][k][j][IBEG]
222  - 1.0/sqrt(R[IBEG]);
223  #if ROTATING_FRAME == YES
224  d->Vc[VX2][k][j][i] += g_OmegaZ*(R[i] - R[IBEG]);
225  #endif
226  #if NTRACER == 1
227  d->Vc[TRC][k][j][i] = 0.0;
228  #endif
229  }
230  }
231 
232  } else if (side == X1_END) {
233 
234  if (box->vpos == CENTER){
235  BOX_LOOP(box,k,j,i){
236  d->Vc[RHO][k][j][i] = d->Vc[RHO][k][j][IEND];
237  #if EOS == IDEAL
238  d->Vc[PRS][k][j][i] = d->Vc[PRS][k][j][IEND];
239  #endif
240  d->Vc[VX1][k][j][i] = d->Vc[VX1][k][j][IEND];
241  d->Vc[VX2][k][j][i] = 1.0/sqrt(R[i]) + d->Vc[VX2][k][j][IEND]
242  - 1.0/sqrt(R[IEND]);
243  #if ROTATING_FRAME == YES
244  d->Vc[VX2][k][j][i] += g_OmegaZ*(R[i] - R[IEND]);
245  #endif
246  }
247  }
248  }
249 }
X1_BEG
#define X1_BEG
Boundary region at X1 beg.
Definition: pluto.h:146
VX2
#define VX2
Definition: mod_defs.h:29
CENTER
#define CENTER
Definition: pluto.h:200
RBOX::vpos
int vpos
Location of the variable inside the cell.
Definition: structs.h:359
RHO
#define RHO
Definition: mod_defs.h:19
BOX_LOOP
#define BOX_LOOP(B, k, j, i)
Definition: macros.h:70
DATA::Vc
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
g_OmegaZ
#define g_OmegaZ
Definition: init.c:64
TRC
#define TRC
Definition: pluto.h:581
X1_END
#define X1_END
Boundary region at X1 end.
Definition: pluto.h:147
VX1
#define VX1
Definition: mod_defs.h:28
IDIR
#define IDIR
Definition: pluto.h:193
j
int j
Definition: analysis.c:2
IEND
long int IEND
Upper grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:37
k
int k
Definition: analysis.c:2
GRID::x
double * x
Definition: structs.h:80
i
int i
Definition: analysis.c:2
IBEG
long int IBEG
Lower grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:35