Magnetized accretion torus. More...

#include "pluto.h"

Include dependency graph for init.c:

Functions
static void	DipoleField (double x1, double x2, double x3, double Bx1, double Bx2, double *A)

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

void	Analysis (const Data d, Grid grid)

void	BackgroundField (double x1, double x2, double x3, double *B0)

void	UserDefBoundary (const Data d, RBox box, int side, Grid *grid)

Variables
static double	kappa

Detailed Description

Magnetized accretion torus.

Magnetized accretion torus in 2D spherical coordinates in a radially stratified atmosphere. The accretion torus is determined by the equilibrium condition

$-\frac{1}{R} + \frac{1}{2(1-a)}\frac{L_k^2}{r^{2-2a}} +(n+1)\frac{p}{\rho} = const$

where $p = K\rho^\gamma$ , R is the spherical radius, r is the cylindrical radius. Since $n = 1/(\gamma-1)$ then $1+n = \gamma/(\gamma-1)$ we can also write using $ a = 0 $ (constant angular momentum distribution, Kuwabara Eq [8]):

$-\frac{1}{R} + \frac{1}{2}\frac{L_k^2}{r^2} + \frac{\gamma}{\gamma-1}K\rho^{\gamma-1} = c = -\frac{1}{r_{\min}} + \frac{1}{2}\frac{L_k^2}{r_{\min}^2}$

where the constant on the right hand side is determined at the zero pressure surface. The previous equation is used to compute K at the point where the density is maximum and then again to obtain the density distribution:

$K = \frac{\gamma-1}{\gamma}\left[ c + \frac{1}{r_{\max}} - \frac{1}{2}\frac{L_k^2}{r_{\max}^2} \right]\frac{1}{\rho_{\max}^{\gamma-1}}\,;\qquad \rho = \left[\frac{\gamma-1}{K\gamma}\left(c + \frac{1}{R} - \frac{1}{2}\frac{L_k^2}{r^2}\right)\right]^{1/(\gamma-1)}$

Thus, specifying $r_{\min}$ , $r_{\max}$ and $\rho_{\max}$ determine the torus structure completely. The torus azimuthal velocity is compute from $v_\phi = L_k/r$ .

The atmosphere is assumed to be isothermal and it is given by the condition of hydrostatic balance:

$\rho_a(R) = \eta\rho_{\max} \exp\left[\left(\frac{1}{R}-\frac{1}{r_{\min}}\right) \frac{GM}{a^2}\right]$

where $\eta$ is the density contrast between the (maximum) torus density and the atmosphere density at $(r_{\min},0)$ .

The torus surface is defined by the condition $p_t > p_a$ .

The dimensionless form of the equation employs on the following units:

$[\rho] = \rho_{\max}$ : reference density;
$[L] = R_{*}$ : star radius;
: reference (squared) velocity.

The user defined parameters of this problem, as they appear in pluto.ini, allows for control in the Torus shape and the contrast of physical values:

g_inputParam[RMIN]: minimum cylindrical radius for the torus (inner rim)
g_inputParam[RMAX]: radius of the Torus where pressure is maximum
g_inputParam[RHO_CUT]: minimum density to define the last contour of the magnetic vec. pot.
g_inputParam[BETA]: plasma beta
g_inputParam[ETA] density contrast between atmosphere and Torus
g_inputParam[SCALE_HEIGHT]: atmosphere scale height

The magnetic field can be specified to be inside the torus or as a large-scale dipole field. In the first case, we give the vector potential:

$A_\phi = B_0(\rho_t - \rho_{\rm cut})$

In the second case, a dipole field is specified in the function DipoleField(). The user-defined constant USE_DIPOLE can be used to select between the two configurations.

Author: A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
T. Matsakos (titos.nosp@m.@odd.nosp@m.job.u.nosp@m.chic.nosp@m.ago.e.nosp@m.du)

Date: Aug 16, 2012

References

"Global Magnetohydrodynamical Simulations of Accretion Tori", Hawley, J.F., ApJ (2000) 528 462.
"The dynamical stability of differentially rotating discs with constant specific angular momentum" Papaloizou, J.C.B, Pringle, J.E, MNRAS (1984) 208 721
"The Acceleration Mechanism of Resistive Magnetohydrodynamic Jets Launched from Accretion Disks." Kuwabara, T., Shibata, K., Kudoh, T., Matsumoto, R. ApJ (2005) 621 921

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

208 {

209 }

void BackgroundField	(	double	x1,
		double	x2,
		double	x3,
		double *	B0
	)

Define the component of a static, curl-free background magnetic field.

Definition at line 212 of file init.c.

 {
    double A;
 
    B0[0] = 0.0;
    B0[1] = 0.0;
    B0[2] = 0.0;
    #if USE_DIPOLE == YES
     DipoleField (x1,x2,x3,B0, B0+1, &A);
    #endif
 }

Here is the call graph for this function:

void DipoleField	(	double	x1,
		double	x2,
		double	x3,
		double *	Bx1,
		double *	Bx2,
		double *	A
	)

static

Assign the polodial component of magnetic fields and vector potential for a dipole.

Definition at line 300 of file init.c.

 {
   double R, z, r, r2, r3, Bmag;
   double a;
 
   Bmag = sqrt(2.0*kappa/g_inputParam[BETA]);
 
   #if GEOMETRY == CYLINDRICAL
    R  = x1; z = x2; r = sqrt(x1*x1 + x2*x2);
    r2 = r*r;
    r3 = r2*r;
    *Bx1 = 3.0*Bmag*z*R/(r2*r3);
    *Bx2 =    -Bmag*(R*R - 3.0*z*z)/(r2*r3);
    *A   = Bmag*R/r3;
   #elif GEOMETRY == SPHERICAL
    r  = x1;
    r2 = r*r;
    r3 = r2*r;
    *Bx1 = 2.0*Bmag*cos(x2)/r3;
    *Bx2 =     Bmag*sin(x2)/r3;
    *A   = Bmag*sin(x2)/r2;
   #endif
 }

Here is the caller graph for this function:

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

 { 
   double rmin, rmax, beta;
   double r, R, z, A_phi, gmmr, H;
   double a,c, l_k, l_k2;
   double Bo, Vo;
   double prt, pra, phi;
   double rhot, rhoa, rhocut;
   
   gmmr = g_gamma/(g_gamma - 1.0);
   #if GEOMETRY == CYLINDRICAL
    r  = x1;   /* cylindrical radius */
    z  = x2;   
    R  = sqrt(x1*x1 + x2*x2);           /* spherical radius */
   #elif GEOMETRY == SPHERICAL
    r  = x1*sin(x2);   /* cylindrical radius */
    z  = x1*cos(x2);   
    R  = x1;           /* spherical radius */
   #endif
 
   phi  = 0.0; 
   
   rmax    = g_inputParam[RMAX];
   rmin    = g_inputParam[RMIN];
   beta    = g_inputParam[BETA];
   rhocut  = g_inputParam[RHO_CUT];  
   H       = g_inputParam[SCALE_HEIGHT];  
 
   l_k     = sqrt(rmax);
   l_k2    = l_k*l_k;
   Vo      = l_k;
 
 /* ----------------------------------------
    Solve for kappa and c assuming pra=0
    ---------------------------------------- */  
  
   c     =    - 1.0/rmin + 0.5*l_k2/(rmin*rmin);
   kappa = (c + 1.0/rmax - 0.5*l_k2/(rmax*rmax))/gmmr;
 
 /* ----------------------------
    Torus density + pressure
    ---------------------------- */  
 
   a    = c + 1.0/R - 0.5*l_k2/(r*r);
   a    = MAX(a, 0.0);
   rhot = pow(a/(gmmr*kappa), 1.0/(g_gamma-1));
   prt  = kappa*pow(rhot,g_gamma);
 
 /* -------------------------------------
     Atmospheric density + pressure 
     (force equilibrium Fg and \nabla P)
    ------------------------------------- */
 
   rhoa = g_inputParam[ETA]*exp((1./R - 1./rmin)/H); 
   pra  = rhoa*H;
 
   Bo = sqrt(2.0*kappa/beta);
 
   if (prt > pra && r > 2.0) {    /* Torus */
     v[RHO] = rhot;
     v[PRS] = prt;
     v[VX1] = 0.0;
     v[VX2] = 0.0;
     v[VX3] = Vo/r;
     v[TRC] = 1.0;
   }else{                          /* Atmosphere */
     v[RHO] = rhoa;
     v[PRS] = pra;
     v[VX1] = 0.0;
     v[VX2] = 0.0;
     v[VX3] = 0.0;
     v[TRC] = 0.0;
   }
   
   #if PHYSICS == MHD || PHYSICS == RMHD
    v[BX1] = 0.0;
    v[BX2] = 0.0;
    v[BX3] = 0.0;
 
    A_phi  = 0.0;
    #if USE_DIPOLE == NO
     if (rhot > rhocut && r > 2.0) A_phi = Bo*(rhot - rhocut);
    #endif
    
    v[AX1] = 0.0;
    v[AX2] = 0.0;
    v[AX3] = A_phi;
   #endif
 
   #if (USE_DIPOLE == YES) && (BACKGROUND_FIELD == NO)
    DipoleField (x1,x2,x3,v+BX1, v+BX2, v+AX3);
   #endif
 
   g_smallPressure = 1.e-8;
 }

Here is the call graph for this function:

void UserDefBoundary	(	const Data *	d,
		RBox *	box,
		int	side,
		Grid *	grid
	)

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

 {
   int     i, j, k, nv;
   double  *x1, *x2, *x3;
   static double vR1[NVAR];
   double R;
 
   x1 = grid[IDIR].x;
   x2 = grid[JDIR].x;
   x3 = grid[KDIR].x;
  
   if (side == X1_BEG){
     if (box->vpos == CENTER){
       BOX_LOOP(box,k,j,i){
         for (nv = 0; nv < NVAR; nv++){
           d->Vc[nv][k][j][i] = d->Vc[nv][k][j][IBEG];
         }
         d->Vc[VX1][k][j][i] = MIN(d->Vc[VX1][k][j][i],0.0); 
       }
     }else if (box->vpos == X2FACE){
       #ifdef STAGGERED_MHD
 //       BOX_LOOP(box,k,j,i) d->Vs[BX2s][k][j][i] = d->Vs[BX2s][k][j][IBEG];
       #endif
     }
   }
 }

Variable Documentation

double kappa

static

Definition at line 101 of file init.c.

Functions

Variables

Detailed Description

Function Documentation

Variable Documentation