Relativistic magnetized jet in axisymmetric coordinates. More...

#include "pluto.h"

Include dependency graph for init.c:

Functions
static void	GetJetValues (double x1, double x2, double x3, double *vj)

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)

Variables
static double	pj

static double	g_gamma

Detailed Description

Relativistic magnetized jet in axisymmetric coordinates.

The jet problem is set up in axisymmetric cylindrical coordinates $ (r,z) $ as described in section 4.2.3 of Mignone, Ugliano & Bodo (2009, [MUB09] hereafter). We label ambient and jet values with the suffix "a" and "j", respectively. The ambient medium is at rest and is characterized by uniform density and pressure, $\rho_a$ and $ p_a $ . The beam enters from the lower-z boundary from a circular nozzle of radius 1, carries a constant poloidal field $ B_z $ and a (radially-varying) toroidal component $B_\phi(R)$ . Flow variables are prescribed as

$\left\{\begin{array}{lcl} \rho(R) &=& \rho_j \\ \noalign{\medskip} v_z(R) &=& v_j \\ \noalign{\medskip} B_\phi(R) &=& \DS \left\{\begin{array}{ll} -\gamma_j b_m R/a & \quad{\rm for}\quad r<a \\ \noalign{\medskip} -\gamma_j b_m a/R & \quad{\rm for}\quad r>a \end{array}\right. \\ \noalign{\medskip} B_z(R) &=& B_{z0}\quad\mathrm{(const)} \end{array}\right.$

Here a is the magnetization radius while $ v_r = B_r = 0 $ . The pressure distribution is found by solving the radial momentum balance between thermal, centrifugal and magnetic forces. Neglecting rotation and assuming Bz to be constant the solution is given by

$p(R) = p_j + b_m^2\left[1-\min\left(\frac{r^2}{a^2},1\right)\right]$

Here $p_j=p_a$ is the jet/ambient pressure at r=1 and its value depends on the Mach number, see [MUB09].

The parameters controlling the problem is

g_inputParam[MACH]: the jet Mach number where $c_s=\sqrt{\Gamma p_j/(\rho h)}$ is the sound speed and $\rho h = \rho + \Gamma p_j/(\Gamma-1)$ . This is used to recover ;
g_inputParam[LORENTZ]: the jet Lorentz factor;
g_inpurParam[RHOJ]: the jet density;
g_inpurParam[SIGMA_POL]: magnetization strength for poloidal magnetic field component (see [MUB09], Eq [65]);
g_inpurParam[SIGMA_TOR]: magnetization strength for toroidal magnetic field component (see [MUB09], Eq [65]);

The different configurations are:

configurations #01-02 employ a purely poloidal field and are similar to section 4.3.3 of Mignone & Bodo (2006);
configuration #03

The following image show the (log) density map at the end of simulation for setup #01.

Density map at the end of the computation for configuration #01

References

"A five-wave Harten-Lax-van Leer Riemann solver for relativistic magnetohydrodynamics", Mignone, Ugliano & Bodo, MNARS (200) 393, 1141
"An HLLC Rieman solver for relativistic flows - II. Magnetohydrodynamics" Mignone & Bodo, MNRAS (2006) 368, 1040

Author: A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)

Date: Sept 17, 2014

Definition in file init.c.

Function Documentation

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

Perform runtime data analysis.

Parameters

[in]	d	the PLUTO Data structure
[in]	grid	pointer to array of Grid structures

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

Definition at line 118 of file init.c.

 {
 
 }

void GetJetValues	(	double	x1,
		double	x2,
		double	x3,
		double *	vj
	)

static

Definition at line 183 of file init.c.

 {
   static int  first_call = 1;
   double bm, b2_av, bphi, lor, a = 0.5;
   double r, x, scrh, sig_z, sig_phi;
 
   r   = x1;
   lor = g_inputParam[LORENTZ];
   sig_z   = g_inputParam[SIGMA_POL];
   sig_phi = g_inputParam[SIGMA_TOR];
   if (fabs(r) < 1.e-9) r = 1.e-9;
 
   x = r/a;
   vj[RHO] = g_inputParam[RHOJ];
 
   EXPAND(vj[VX1] = 0.0;                        ,
          vj[VX2] = sqrt(1.0 - 1.0/(lor*lor));  , /* 3-vel */
          vj[VX3] = 0.0;)
 
   scrh = g_inputParam[MACH]/vj[VX2];
   pj   = vj[RHO]/(g_gamma*scrh*scrh - g_gamma/(g_gamma - 1.0));
 
   bm  = 4.0*pj*sig_phi;
   bm /= a*a*(1.0 - 4.0*log(a) - 2.0*sig_phi);
   bm  = sqrt(bm);
 
   scrh    = MIN(x*x, 1.0);
   vj[PRS] = pj + bm*bm*(1.0 - scrh);
 
   bphi  = bm*(fabs(x) < 1.0 ? x: 1.0/x);
 
   #if GEOMETRY == CYLINDRICAL
    EXPAND(vj[iBR]   = 0.0;                                 ,
           vj[iBZ]   = sqrt(sig_z*(bm*bm*a*a + 2.0*pj));  ,
           vj[iBPHI] = lor*bphi;)
    vj[AX1] = 0.0;
    vj[AX2] = 0.0;
    vj[AX3] = 0.5*r*vj[iBZ];
   #endif
 
 /* -- set tracer value -- */
 
   #if NTRACER > 0
    vj[TRC] = 1.0;
   #endif
   #ifdef PSI_GLM 
    vj[PSI_GLM] = 0.0;
   #endif
 
 }

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]	v	a pointer to a vector of primitive variables
[in]	x1	coordinate point in the 1st dimension
[in]	x2	coordinate point in the 2nd dimension
[in]	x3	coordinate 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 80 of file init.c.

 {
   int    nv;
   double r, z, p0, vjet[256], vamb[256];
 
   g_gamma = 5./3.;
   GetJetValues (x1, x2, x3, vjet);
 
 /* -- set ambient values -- */
 
   v[RHO] = g_inputParam[RHOA];
   EXPAND(v[VX1] = 0.0;   ,
          v[VX2] = 0.0;   ,
          v[VX3] = 0.0;)
   v[PRS] = pj;
   #if PHYSICS == RMHD
    EXPAND(v[BX1] = 0.0;        ,
           v[BX2] = vjet[BX2];  ,
           v[BX3] = 0.0;)
    v[AX1] = 0.0;
    v[AX2] = 0.0;
    v[AX3] = vjet[AX3];
   #endif
   #if NTRACER > 0
    v[TRC] = 0.0;
   #endif
   #ifdef PSI_GLM 
    v[PSI_GLM] = 0.0;
   #endif
 
   g_smallPressure = 1.e-5;
 }

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]	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]	box	pointer 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]	side	specifies 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]	grid	pointer 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]	box	pointer 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]	side	specifies 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]	grid	pointer 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 128 of file init.c.

 {
   int   i, j, k, nv;
   double *x1, *x2, *x3;
   double ***bxs, ***bys, ***bzs;
   double prof, vjet[256], vout[NVAR];
 
   #ifdef STAGGERED_MHD
    D_EXPAND(bxs = d->Vs[BX1s];  , 
             bys = d->Vs[BX2s];  , 
             bzs = d->Vs[BX3s];)
   #endif
 
   x1  = grid[IDIR].xgc;  
   x2  = grid[JDIR].xgc;  
   x3  = grid[KDIR].xgc;
 
   if (side == X2_BEG){  /* -- X2_BEG boundary -- */
     if (box->vpos == CENTER){    /* -- cell-centered boundary conditions -- */
       BOX_LOOP(box,k,j,i){
         GetJetValues (x1[i], x2[j], x3[k], vjet); /* Jet Values */
         VAR_LOOP(nv) vout[nv] = d->Vc[nv][k][2*JBEG - j - 1][i]; /* Ambient */
 
         vout[VX2] *= -1.0; 
         EXPAND(vout[BX1] *= -1.0;  ,  
                                 ;  , 
                vout[BX3] *= -1.0;)
         #ifdef PSI_GLM
          vjet[PSI_GLM]  = d->Vc[PSI_GLM][k][JBEG][i] - d->Vc[BX2][k][JBEG][i];
          vout[PSI_GLM] *= -1.0;
         #endif
 
         prof = (fabs(x1[i]) <= 1.0); 
         VAR_LOOP(nv) d->Vc[nv][k][j][i] = vout[nv] - (vout[nv] - vjet[nv])*prof;
       }
 
     }else if (box->vpos == X1FACE){  /* -- staggered fields -- */
       #ifdef STAGGERED_MHD
        x1 = grid[IDIR].xr;
        vjet[BX1] = 0.0; 
        BOX_LOOP(box,k,j,i){
           vout[BX1] = -bxs[k][2*JBEG - j - 1][i];
           prof = (fabs(x1[i]) <= 1.0);
           bxs[k][j][i] = vout[BX1] - (vout[BX1] - vjet[BX1])*prof;
        }
       #endif
     }
   }
 }

Here is the call graph for this function:

Variable Documentation

double g_gamma

static

Definition at line 76 of file init.c.

double pj

static

Definition at line 73 of file init.c.

Functions

Variables

Detailed Description

Function Documentation

Variable Documentation