PLUTO
Functions | Variables
init.c File Reference

Circularly polarized Alfven waves. More...

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

Go to the source code of this file.

Functions

static double AX_VEC (double x, double y, double z)
 
static double AY_VEC (double x, double y, double z)
 
static double AZ_VEC (double x, double y, double z)
 
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)
 

Variables

static double ta = 0.0
 
static double tb = 0.0
 
static double ca
 
static double sa
 
static double cg
 
static double sg
 

Detailed Description

Circularly polarized Alfven waves.

This is a rotated version of a 1D setup where

\[ \rho = 1 \,,\quad v_x = V_0 \,,\quad v_y = |\epsilon|\cos(\phi) \,,\quad v_z = |\epsilon|\sin(\phi) \,,\quad B_x = 1 \,,\quad B_y = \epsilon\cos(\phi) \,,\quad B_z = \epsilon\cos(\phi) \,,\quad P = P_0 \,, \]

where $V_0$ is the translation velocity in the 1D $x$ direction, $\psi$ is the phase ( $k_xx$ in 1D and ( $k_xx + k_yy$) in 2D), and $\epsilon$ is the wave amplitude ( $\epsilon > 0$ implies right going waves; $\epsilon < 0$ implies left going waves).

With this normalization, the Alfven velocity $V_A = B_x/\sqrt{\rho}$ is always one.

Rotations are specified by $t_a = k_y/k_x$ and $t_b = k_z/k_x$ expressing the ratios between the $y$- and $z$- components of the wave vector with the $x$ component. We choose the wavelength in $x$-direction to be always 1, so that $k_x = 2\pi/l_x = 2\pi$, $k_y = 2\pi/l_y$, $k_z = 2\pi/l_z$ so that $t_z = 1/l_y$, $t_y = 1/l_z$. By choosing the domain in the $x$-direction to be 1, the extent in $y$ and $z$ should be adjusted so that $L_y/l_y = 1,2,3,4...$ (and similarly for $L_z/l_z$) becomes an integer number to ensure periodicity. If one wants 1 wave length in both directions than $L_x = 1$, $L_y = 1/t_a$, $L_z = 1/t_b$.

If $N$ wavelengths are specified in the $y$-direction than $L_x = 1$, $L_y = N/t_a$, $L_z = N/t_b$

1D is recovered by specifying $t_z = t_y = 0$.

The final time step is one period and is found from

$V_A = \omega/|k|$ –> $1/T = \sqrt{1 + t_a^2 + t_b^2}$

The runtime parameters that are read from pluto.ini are

  • g_inputParam[EPS]: sets the wave amplitude $\epsilon$;
  • g_inputParam[VEL0]: sets $V_0$;
  • g_inputParam[PR0]: sets the pressure of the 1D solution;
  • g_inputParam[ALPHA_GLM]: ;

Configurations:

  • #01-05, 08-09: 2D cartesian;
  • #06-07: 3D cartesian;
Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
July 09, 2014

References:

  • "High-order conservative finite difference GLM-MHD schemes for cell-centered MHD" Mignone, Tzeferacos & Bodo JCP (2010)

Definition in file init.c.

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 151 of file init.c.

156 {
157  int i,j,k;
158  double bmag, err[3], gerr[3], Vtot, dV;
159  double ***Bx, ***By, ***Bz;
160  double *dx, *dy, *dz;
161  static double ***Bx0, ***By0, ***Bz0;
162  FILE *fp;
163 
164  Bx = d->Vc[BX1]; By = d->Vc[BX2]; Bz = d->Vc[BX3];
165  dx = grid[IDIR].dx; dy = grid[JDIR].dx; dz = grid[KDIR].dx;
166 
167  if (Bx0 == NULL){
168  Bx0 = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
169  By0 = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
170  Bz0 = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
171  print1 ("ANALYSIS: Initializing bmag\n");
172  DOM_LOOP(k,j,i){
173  Bx0[k][j][i] = Bx[k][j][i];
174  By0[k][j][i] = By[k][j][i];
175  Bz0[k][j][i] = Bz[k][j][i];
176  }
177  return;
178  }
179 
180  for (i = 0; i < 3; i++) err[i] = 0.0;
181  DOM_LOOP(k,j,i){
182  dV = D_EXPAND(dx[i], *dy[j], *dz[k]);
183  err[0] += fabs(Bx[k][j][i] - Bx0[k][j][i])*dV;
184  err[1] += fabs(By[k][j][i] - By0[k][j][i])*dV;
185  err[2] += fabs(Bz[k][j][i] - Bz0[k][j][i])*dV;
186  }
187 
188  Vtot = D_EXPAND(
189  (g_domEnd[IDIR] - g_domBeg[IDIR]),
190  *(g_domEnd[JDIR] - g_domBeg[JDIR]),
191  *(g_domEnd[KDIR] - g_domBeg[KDIR]));
192 
193  for (i = 0; i < 3; i++) err[i] /= Vtot;
194 
195  #ifdef PARALLEL
196  MPI_Allreduce (err, gerr, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
197  for (i = 0; i < 3; i++) err[i] = gerr[i];
198  MPI_Barrier (MPI_COMM_WORLD);
199  #endif
200 
201  if (prank == 0){
202  double errtot;
203  errtot = sqrt(err[0]*err[0] + err[1]*err[1] + err[2]*err[2]);
204  fp = fopen("bmag.dat","w");
205  printf ("analysis: writing to disk...\n");
206  fprintf (fp,"%d %10.3e\n", (int)(1.0/dx[IBEG]), errtot);
207  fclose (fp);
208  }
209 }
DOM_LOOP
DOM_LOOP(k, j, i)
Definition: analysis.c:22
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
NX2_TOT
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
Definition: globals.h:57
DATA::Vc
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
GRID::dV
double * dV
Cell volume.
Definition: structs.h:86
ARRAY_3D
#define ARRAY_3D(nx, ny, nz, type)
Definition: prototypes.h:172
GRID::dx
double * dx
Definition: structs.h:83
prank
int prank
Processor rank.
Definition: globals.h:33
KDIR
#define KDIR
Definition: pluto.h:195
IDIR
#define IDIR
Definition: pluto.h:193
Bx
static double Bx
Definition: hlld.c:62
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
k
int k
Definition: analysis.c:2
g_domBeg
double g_domBeg[3]
Lower limits of the computational domain.
Definition: globals.h:125
BX3
#define BX3
Definition: mod_defs.h:27
D_EXPAND
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
Definition: analysis.c:27
i
int i
Definition: analysis.c:2
fp
FILE * fp
Definition: analysis.c:7
BX1
#define BX1
Definition: mod_defs.h:25
BX2
#define BX2
Definition: mod_defs.h:26
g_domEnd
double g_domEnd[3]
Upper limits of the computational domain.
Definition: globals.h:126
JDIR
#define JDIR
Definition: pluto.h:194
NX1_TOT
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.
Definition: globals.h:55
IBEG
long int IBEG
Lower grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:35

Here is the call graph for this function:

double AX_VEC ( double  x,
double  y,
double  z 
)
static

Definition at line 218 of file init.c.

219 {
220  double xn, yn, zn, eps;
221  double kx, ky, kz, phi;
222  double kxn, kyn, kzn;
223  double Axn, Ayn, Azn;
224 
225  eps = g_inputParam[EPS];
226  kx = 2.0*CONST_PI/1.0;
227  ky = ta*kx;
228  kz = tb*kx;
229 
230  phi = D_EXPAND(kx*x, + ky*y, + kz*z);
231 
232  xn = ca*cg*x + sa*cg*y + sg*z;
233  yn = -sa*x + ca*y;
234  zn = -ca*sg*x - sa*sg*y + cg*z;
235 
236  kxn = ca*cg*kx + sa*cg*ky + sg*kz;
237  kyn = -sa*kx + ca*ky;
238  kzn = -ca*sg*kx - sa*sg*ky + cg*kz;
239 
240  Axn = 0.0;
241  Ayn = eps*sin(phi)/kxn;
242  Azn = eps*cos(phi)/kxn + yn;
243 
244  return(Axn*cg*ca - Ayn*sa - Azn*sg*ca);
245 }
EPS
#define EPS
Definition: definitions_01.h:26
cg
static double cg
Definition: init.c:75
ca
static double ca
Definition: init.c:75
ta
static double ta
Definition: init.c:74
tb
static double tb
Definition: init.c:74
eps
#define eps
g_inputParam
double g_inputParam[32]
Array containing the user-defined parameters.
Definition: globals.h:131
GRID::x
double * x
Definition: structs.h:80
sg
static double sg
Definition: init.c:75
D_EXPAND
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
Definition: analysis.c:27
CONST_PI
#define CONST_PI
.
Definition: pluto.h:269
sa
static double sa
Definition: init.c:75

Here is the call graph for this function:

Here is the caller graph for this function:

double AY_VEC ( double  x,
double  y,
double  z 
)
static

Definition at line 247 of file init.c.

248 {
249  double xn, yn, zn, eps;
250  double kx, ky, kz, phi;
251  double kxn, kyn, kzn;
252  double Axn, Ayn, Azn;
253 
254  eps = g_inputParam[EPS];
255  kx = 2.0*CONST_PI/1.0;
256  ky = ta*kx;
257  kz = tb*kx;
258 
259  phi = D_EXPAND(kx*x, + ky*y, + kz*z);
260 
261  xn = ca*cg*x + sa*cg*y + sg*z;
262  yn = -sa*x + ca*y;
263  zn = -ca*sg*x - sa*sg*y + cg*z;
264 
265  kxn = ca*cg*kx + sa*cg*ky + sg*kz;
266  kyn = -sa*kx + ca*ky;
267  kzn = -ca*sg*kx - sa*sg*ky + cg*kz;
268 
269  Axn = 0.0;
270  Ayn = eps*sin(phi)/kxn;
271  Azn = eps*cos(phi)/kxn + yn;
272 
273  return(Axn*cg*sa + Ayn*ca - Azn*sg*sa);
274 }
EPS
#define EPS
Definition: definitions_01.h:26
cg
static double cg
Definition: init.c:75
ca
static double ca
Definition: init.c:75
ta
static double ta
Definition: init.c:74
tb
static double tb
Definition: init.c:74
eps
#define eps
g_inputParam
double g_inputParam[32]
Array containing the user-defined parameters.
Definition: globals.h:131
GRID::x
double * x
Definition: structs.h:80
sg
static double sg
Definition: init.c:75
D_EXPAND
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
Definition: analysis.c:27
CONST_PI
#define CONST_PI
.
Definition: pluto.h:269
sa
static double sa
Definition: init.c:75

Here is the call graph for this function:

Here is the caller graph for this function:

double AZ_VEC ( double  x,
double  y,
double  z 
)
static

Definition at line 275 of file init.c.

276 {
277  double xn, yn, zn, eps;
278  double kx, ky, kz, phi;
279  double kxn, kyn, kzn;
280  double Axn, Ayn, Azn;
281 
282  eps = g_inputParam[EPS];
283  kx = 2.0*CONST_PI/1.0;
284  ky = ta*kx;
285  kz = tb*kx;
286 
287  phi = D_EXPAND(kx*x, + ky*y, + kz*z);
288 
289  xn = ca*cg*x + sa*cg*y + sg*z;
290  yn = -sa*x + ca*y;
291  zn = -ca*sg*x - sa*sg*y + cg*z;
292 
293  kxn = ca*cg*kx + sa*cg*ky + sg*kz;
294  kyn = -sa*kx + ca*ky;
295  kzn = -ca*sg*kx - sa*sg*ky + cg*kz;
296 
297  Axn = 0.0;
298  Ayn = eps*sin(phi)/kxn;
299  Azn = eps*cos(phi)/kxn + yn;
300 
301  return(Axn*sg + Azn*cg);
302 }
EPS
#define EPS
Definition: definitions_01.h:26
cg
static double cg
Definition: init.c:75
ca
static double ca
Definition: init.c:75
ta
static double ta
Definition: init.c:74
tb
static double tb
Definition: init.c:74
eps
#define eps
g_inputParam
double g_inputParam[32]
Array containing the user-defined parameters.
Definition: globals.h:131
GRID::x
double * x
Definition: structs.h:80
sg
static double sg
Definition: init.c:75
D_EXPAND
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
Definition: analysis.c:27
CONST_PI
#define CONST_PI
.
Definition: pluto.h:269
sa
static double sa
Definition: init.c:75

Here is the call graph for this function:

Here is the caller 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 78 of file init.c.

82 {
83  static int first_call = 1;
84  double tg, eps;
85  double x, y, z, kx, ky,kz;
86  double phi, cb, sb;
87  double vx, vy, vz, Bx, By, Bz, rho, pr;
88  double Lx, Ly, Lz;
89 
90  Lx = g_domEnd[IDIR] - g_domBeg[IDIR];
91  Ly = g_domEnd[JDIR] - g_domBeg[JDIR];
92  Lz = g_domEnd[KDIR] - g_domBeg[KDIR];
93 
94  if (first_call == 1){
95  D_EXPAND( ,
96  ta = Lx/Ly; ,
97  tb = Lx/Lz;)
98 
99  if (prank == 0) printf ("ta = %f; tb = %f\n",ta, tb);
100  first_call = 0;
101  }
102  x = x1; y = x2; z = x3;
103 
104  eps = g_inputParam[EPS]; /* wave amplitude */
105  kx = 2.0*CONST_PI/1.0;
106  ky = ta*kx;
107  kz = tb*kx;
108 
109  ca = 1.0/sqrt(ta*ta + 1.0);
110  cb = 1.0/sqrt(tb*tb + 1.0);
111  sa = ta*ca;
112  sb = tb*cb;
113 
114  phi = D_EXPAND(kx*x, + ky*y, + kz*z);
115 
116 /* -- define the 1-D solution -- */
117 
118  rho = 1.0;
119  pr = g_inputParam[PR0];
120  vx = g_inputParam[VEL0]; /* translational velocity */
121  vy = fabs(eps)*sin(phi);
122  vz = fabs(eps)*cos(phi);
123  Bx = 1.0;
124  By = eps*sqrt(rho)*sin(phi);
125  Bz = eps*sqrt(rho)*cos(phi);
126 
127 /* -- rotate solution -- */
128 
129  us[RHO] = rho;
130  us[PRS] = pr;
131 
132  tg = tb*ca;
133  cg = 1.0/sqrt(tg*tg + 1.0); sg = cg*tg;
134 
135  us[VX1] = vx*cg*ca - vy*sa - vz*sg*ca;
136  us[VX2] = vx*cg*sa + vy*ca - vz*sg*sa;
137  us[VX3] = vx*sg + vz*cg;
138  us[TRC] = 0.0;
139 
140  us[BX1] = Bx*cg*ca - By*sa - Bz*sg*ca;
141  us[BX2] = Bx*cg*sa + By*ca - Bz*sg*sa;
142  us[BX3] = Bx*sg + Bz*cg;
143 
144  #ifdef STAGGERED_MHD
145  us[AX1] = AX_VEC(x,y,z);
146  us[AX2] = AY_VEC(x,y,z);
147  us[AX3] = AZ_VEC(x,y,z);
148  #endif
149 }
VX2
#define VX2
Definition: mod_defs.h:29
n
static int n
Definition: analysis.c:3
VEL0
#define VEL0
Definition: definitions_01.h:27
EPS
#define EPS
Definition: definitions_01.h:26
RHO
#define RHO
Definition: mod_defs.h:19
cg
static double cg
Definition: init.c:75
ca
static double ca
Definition: init.c:75
PR0
#define PR0
Definition: definitions_01.h:28
ta
static double ta
Definition: init.c:74
AX2
#define AX2
Definition: mod_defs.h:86
TRC
#define TRC
Definition: pluto.h:581
AX_VEC
static double AX_VEC(double x, double y, double z)
Definition: init.c:218
prank
int prank
Processor rank.
Definition: globals.h:33
VX1
#define VX1
Definition: mod_defs.h:28
KDIR
#define KDIR
Definition: pluto.h:195
tb
static double tb
Definition: init.c:74
eps
#define eps
IDIR
#define IDIR
Definition: pluto.h:193
AX3
#define AX3
Definition: mod_defs.h:87
Bx
static double Bx
Definition: hlld.c:62
g_inputParam
double g_inputParam[32]
Array containing the user-defined parameters.
Definition: globals.h:131
AZ_VEC
static double AZ_VEC(double x, double y, double z)
Definition: init.c:275
g_domBeg
double g_domBeg[3]
Lower limits of the computational domain.
Definition: globals.h:125
GRID::x
double * x
Definition: structs.h:80
sg
static double sg
Definition: init.c:75
BX3
#define BX3
Definition: mod_defs.h:27
D_EXPAND
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
Definition: analysis.c:27
AX1
#define AX1
Definition: mod_defs.h:85
BX1
#define BX1
Definition: mod_defs.h:25
VX3
#define VX3
Definition: mod_defs.h:30
CONST_PI
#define CONST_PI
.
Definition: pluto.h:269
sa
static double sa
Definition: init.c:75
BX2
#define BX2
Definition: mod_defs.h:26
g_domEnd
double g_domEnd[3]
Upper limits of the computational domain.
Definition: globals.h:126
JDIR
#define JDIR
Definition: pluto.h:194
AY_VEC
static double AY_VEC(double x, double y, double z)
Definition: init.c:247

Here is the call graph for this function:

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.

Assign user-defined boundary conditions in the lower boundary ghost zones. The profile is top-hat:

\[ V_{ij} = \left\{\begin{array}{ll} V_{\rm jet} & \quad\mathrm{for}\quad r_i < 1 \\ \noalign{\medskip} \mathrm{Reflect}(V) & \quad\mathrm{otherwise} \end{array}\right. \]

where $ V_{\rm jet} = (\rho,v,p)_{\rm jet} = (1,M,1/\Gamma)$ and M is the flow Mach number (the unit velocity is the jet sound speed, so $ v = M$).

Assign user-defined boundary conditions:

  • left side (x beg): constant shocked values
  • bottom side (y beg): constant shocked values for x < 1/6 and reflective boundary otherwise.
  • top side (y end): time-dependent boundary: for $ x < x_s(t) = 10 t/\sin\alpha + 1/6 + 1.0/\tan\alpha $ we use fixed (post-shock) values. Unperturbed values otherwise.

Assign user-defined boundary conditions at inner and outer radial boundaries. Reflective conditions are applied except for the azimuthal velocity which is fixed.

Definition at line 211 of file init.c.

215 { }

Variable Documentation

double ca
static

Definition at line 75 of file init.c.

double cg
static

Definition at line 75 of file init.c.

double sa
static

Definition at line 75 of file init.c.

double sg
static

Definition at line 75 of file init.c.

double ta = 0.0
static

Definition at line 74 of file init.c.

double tb = 0.0
static

Definition at line 74 of file init.c.