PLUTO
Functions
init.c File Reference

Relativistic jet propagation. More...

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

Go to the source code of this file.

Functions

static double Profile (double r, int nv)
 
static void GetJetValues (double *vjet)
 
void Init (double *v, double x1, double x2, double x3)
 
void Analysis (const Data *d, Grid *grid)
 
void UserDefBoundary (const Data *d, RBox *box, int side, Grid *grid)
 

Detailed Description

Relativistic jet propagation.

This problem sets initial and boundary conditions for a jet propagating into a uniform medium with constant density and pressure. The computation can be carried using CYLINDRICAL or SPHERICAL coordinates. In the first case, the jet enters at the lower z-boundary with speed $ v_z = \beta$ while in the second case, a conical beam with a small aperture ( $\theta = 5^\circ$) is injected with the same speed from the lower radial boundary. At the border of the nozzle, jet values are smoothly joined with ambient values using the Profile() function (in the current setting the profile is a sharp transition). The jet is pressure-matched so that the beam and ambient pressure coincide.

The configuration is defined in terms of the following parameters:

  1. g_inputParam[BETA]: the jet velocity;
  2. g_inputParam[RHO_IN]: the jet density;
  3. g_inputParam[RHO_OUT]: the ambient density.
  4. g_inputParam[PRESS_IN]: the jet pressure (also equal to ambient pressure)

defined in pluto.ini. The TAUB equation of state is used.

  • Configurations #01 and #02 use CYLINDRICAL coordinates;
  • Configuration #03 employs SPHERICAL coordinates (see snapshot below)
rhd_jet.03.jpg
Density (log) for configuration #03 at t=200
Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
Sept 18, 2014

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

Compute volume-integrated magnetic pressure, Maxwell and Reynolds stresses. Save them to "averages.dat"

Definition at line 66 of file init.c.

71 {
72 
73 }
void GetJetValues ( double *  vjet)
static

Definition at line 114 of file init.c.

119 {
120  vjet[RHO] = g_inputParam[RHO_IN];
121  #if GEOMETRY == CYLINDRICAL || GEOMETRY == CARTESIAN
122  vjet[VX1] = 0.0;
123  vjet[VX2] = g_inputParam[BETA];
124  #elif GEOMETRY == SPHERICAL
125  vjet[VX1] = g_inputParam[BETA];
126  vjet[VX2] = 0.0;
127  #endif
128  vjet[PRS] = g_inputParam[PRESS_IN];
129 }
VX2
#define VX2
Definition: mod_defs.h:29
BETA
#define BETA
Definition: definitions_01.h:26
RHO_IN
#define RHO_IN
Definition: definitions_01.h:23
RHO
#define RHO
Definition: mod_defs.h:19
VX1
#define VX1
Definition: mod_defs.h:28
g_inputParam
double g_inputParam[32]
Array containing the user-defined parameters.
Definition: globals.h:131
PRESS_IN
#define PRESS_IN
Definition: definitions_01.h:23
void Init ( double *  v,
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 46 of file init.c.

52 {
53  double scrh;
54  #if EOS == IDEAL
55  g_gamma = 5.0/3.0;
56  #endif
57 
58  v[RHO] = g_inputParam[RHO_OUT];
59  v[VX1] = 0.0;
60  v[VX2] = 0.0;
61  v[PRS] = g_inputParam[PRESS_IN];
62 
63  g_smallPressure = v[PRS]/500.0;
64 }
g_gamma
double g_gamma
Definition: globals.h:112
VX2
#define VX2
Definition: mod_defs.h:29
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
g_smallPressure
double g_smallPressure
Small value for pressure fix.
Definition: globals.h:110
RHO_OUT
#define RHO_OUT
Definition: definitions_01.h:28
VX1
#define VX1
Definition: mod_defs.h:28
g_inputParam
double g_inputParam[32]
Array containing the user-defined parameters.
Definition: globals.h:131
PRESS_IN
#define PRESS_IN
Definition: definitions_01.h:23
double Profile ( double  r,
int  nv 
)
static

Definition at line 132 of file init.c.

137 {
138  int xn = 14;
139  double r0 = 1.0;
140 
141  if (nv == RHO) r0 = 1.1;
142 
143  #if GEOMETRY == SPHERICAL
144  r0 = 5.0/180.0*CONST_PI;
145  #endif
146  return 1.0/cosh(pow(r/r0,xn));
147 }
RHO
#define RHO
Definition: mod_defs.h:19
CONST_PI
#define CONST_PI
.
Definition: pluto.h:269

Here is the caller 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.

Assign user-defined boundary conditions.

Parameters
[in/out]d pointer 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 on which side boundary conditions need to be assigned. side 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.

Set the injection boundary condition at the lower z-boundary (X2-beg must be set to userdef in pluto.ini). For $ R \le 1 $ we set constant input values (given by the GetJetValues() function while for $ R > 1 $ the solution has equatorial symmetry with respect to the z=0 plane. To avoid numerical problems with a "top-hat" discontinuous jump, we smoothly merge the inlet and reflected value using a profile function Profile().

Provide inner radial boundary condition in polar geometry. Zero gradient is prescribed on density, pressure and magnetic field. For the velocity, zero gradient is imposed on v/r (v = vr, vphi).

The user-defined boundary is used to impose stress-free boundary and purely vertical magnetic field at the top and bottom boundaries, as done in Bodo et al. (2012). In addition, constant temperature and hydrostatic balance are imposed. For instance, at the bottom boundary, one has:

\[ \left\{\begin{array}{lcl} \DS \frac{dp}{dz} &=& \rho g_z \\ \noalign{\medskip} p &=& \rho c_s^2 \end{array}\right. \qquad\Longrightarrow\qquad \frac{p_{k+1}-p_k}{\Delta z} = \frac{p_{k} + p_{k+1}}{2c_s^2}g_z \]

where $g_z$ is the value of gravity at the lower boundary. Solving for $p_k$ at the bottom boundary where $k=k_b-1$ gives:

\[ \left\{\begin{array}{lcl} p_k &=& \DS p_{k+1} \frac{1-a}{1+a} \\ \noalign{\medskip} \rho_k &=& \DS \frac{p_k}{c_s^2} \end{array}\right. \qquad\mathrm{where}\qquad a = \frac{\Delta z g_z}{2c_s^2} > 0 \]

where, for simplicity, we keep constant temperature in the ghost zones rather than at the boundary interface (this seems to give a more stable behavior and avoids negative densities). A similar treatment holds at the top boundary.

Assign user-defined boundary conditions. At the inner boundary we use outflow conditions, except for velocity which we reset to zero when there's an inflow

Definition at line 75 of file init.c.

80 {
81  int nv, i, j, k;
82  double r, vjet[NVAR], vout[NVAR];
83 
84  #if GEOMETRY == SPHERICAL
85  if (side == X1_BEG){
86  GetJetValues(vjet);
87  X1_BEG_LOOP(k,j,i){
88  VAR_LOOP(nv) vout[nv] = d->Vc[nv][k][j][2*IBEG-i-1];
89  vout[VX2] *= -1.0;
90 
91  r = grid[JDIR].x[j];
92  VAR_LOOP(nv)
93  d->Vc[nv][k][j][i] = vout[nv] + (vjet[nv] - vout[nv])*Profile(r,nv);
94  }
95  }
96  #endif
97 
98  #if GEOMETRY == CYLINDRICAL || GEOMETRY == CARTESIAN
99  if (side == X2_BEG){
100  GetJetValues(vjet);
101  X2_BEG_LOOP(k,j,i){
102  VAR_LOOP(nv) vout[nv] = d->Vc[nv][k][2*JBEG-j-1][i];
103  vout[VX2] *= -1.0;
104 
105  r = grid[IDIR].x[i];
106  for (nv = 0; nv < NVAR; nv++)
107  d->Vc[nv][k][j][i] = vout[nv] + (vjet[nv] - vout[nv])*Profile(r,nv);
108  }
109  }
110  #endif
111 }
X1_BEG
#define X1_BEG
Boundary region at X1 beg.
Definition: pluto.h:146
VX2
#define VX2
Definition: mod_defs.h:29
X2_BEG_LOOP
#define X2_BEG_LOOP(k, j, i)
Definition: macros.h:47
VAR_LOOP
#define VAR_LOOP(n)
Definition: macros.h:226
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
X1_BEG_LOOP
#define X1_BEG_LOOP(k, j, i)
Definition: macros.h:46
i
int i
Definition: analysis.c:2
GetJetValues
void GetJetValues(double, double *)
Definition: init.c:264
X2_BEG
#define X2_BEG
Boundary region at X2 beg.
Definition: pluto.h:148
JDIR
#define JDIR
Definition: pluto.h:194
JBEG
long int JBEG
Lower grid index of the computational domain in the the X2 direction for the local processor...
Definition: globals.h:39
NVAR
#define NVAR
Definition: pluto.h:609
IBEG
long int IBEG
Lower grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:35
Profile
static double Profile(double r, int nv)
Definition: init.c:132

Here is the call graph for this function: