MHD Shock Cloud interaction. More...

#include "pluto.h"

Include dependency graph for init.c:

Functions
void	Init (double *v, double x, double y, double z)

void	Analysis (const Data d, Grid grid)

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

Detailed Description

MHD Shock Cloud interaction.

Set the initial condition for the 2D or 3D shock-cloud problem. The shock propagates in the x-direction and it is initially located at x=0.6 with post- and -pre shock value equal to

$\left(\begin{array}{l} \rho \\ \noalign{\medskip} v_x \\ \noalign{\medskip} p \\ \noalign{\medskip} B_y \\ \noalign{\medskip} B_z \end{array}\right) = \left(\begin{array}{l} 3.86859 \\ \noalign{\medskip} 0 \\ \noalign{\medskip} 167.345 \\ \noalign{\medskip} B_{\rm post} \\ \noalign{\medskip} -B_{\rm post} \\ \noalign{\medskip} \end{array}\right) \quad\mathrm{for}\quad x < 0.6\,,\qquad\qquad \left(\begin{array}{l} \rho \\ \noalign{\medskip} v_x \\ \noalign{\medskip} p \\ \noalign{\medskip} B_y \\ \noalign{\medskip} B_z \end{array}\right) = \left(\begin{array}{l} 1 \\ \noalign{\medskip} -11.2536 \\ \noalign{\medskip} 1 \\ \noalign{\medskip} B_{\rm pre} \\ \noalign{\medskip} -B_{\rm pre} \\ \noalign{\medskip} \end{array}\right) \quad\mathrm{for}\quad x > 0.6\,,\qquad$

while the remaining vector components are 0. The cloud has radius $R$ centered at $r_0 = (x_0,y_0,z_0) = (0.8, 0.5, 0.5)$ and larger density $\rho = 10$ .

The interaction may be divided into two phases: 1) the collapse stage where the front of the cloud is strongly compressed and two fast shocks are generated and 2) the reexpansion phase which begins when the transmitted fast shock overtake the back of the cloud.

The runtime parameters that are read from pluto.ini are

g_inputParam[B_POST]: sets the post-shock magnetic field;
g_inputParam[B_PRE]: sets the pre-shock magnetic field;
g_inputParam[RADIUS]: sets the radius of the cloud;

A list of the available configurations is given in the following table:

Conf.	GEOMETRY	DIM	T. STEPPING	RECONSTRUCTION	divB	AMR
#01	CARTESIAN	2	RK2	LINEAR	CT	NO
#02	CARTESIAN	2	ChTr	LINEAR	CT	NO
#03	CARTESIAN	2	HANCOCK	LINEAR	8W	NO
#04	CARTESIAN	3	ChTr	LINEAR	CT	NO
#05	CARTESIAN	3	RK3	LINEAR	CT	NO
#06	CARTESIAN	3	ChTr	LINEAR	GLM	NO
#07	CARTESIAN	3	RK3	WENO3_FD	GLM	NO
#08	CARTESIAN	3	HANCOK	LINEAR	GLM	YES
#09	CARTESIAN	2	RK2	LINEAR	CT	NO

Density map with overplotted field lines for configuration #02

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

Date: Sept 17, 2014

References:

"Interactions between magnetohydrodynamical shocks and denser clouds" W. Dai an P.R. Woodward, ApJ (1994) 436, 776.
"The divB = 0 Constraint in Shock-Capturing MHD codes" G. Toth, JCP 161, 605-652 (2000)
"A central constrained transport scheme for ideal magnetohydrodynamics" U. Ziegler, JCP 196 (2004), 393

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

142 {

143 }

void Init	(	double *	v,
		double	x,
		double	y,
		double	z
	)

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

 {
  double r, x0, y0, z0;
 
   x0 = 0.8;
   y0 = 0.5;
   z0 = 0.5;
 
   if (x < 0.6) {
     v[RHO] = 3.86859;
     v[VX1] = 0.0;
     v[VX2] = 0.0;
     v[VX3] = 0.0;
     v[BX1] = 0.0;
     v[BX2] =  g_inputParam[B_POST];
     v[BX3] = -g_inputParam[B_POST];
     v[PRS] = 167.345;
   }else{
     v[RHO] = 1.0;
     v[VX1] = -11.2536;
     v[VX2] = 0.0;
     v[VX3] = 0.0;
     v[BX1] = 0.0;
     v[BX2] = g_inputParam[B_PRE];
     v[BX3] = g_inputParam[B_PRE];
     v[PRS] = 1.0;
   } 
 
   /*  ----  CLOUD  ----  */
 
   r = D_EXPAND((x-x0)*(x-x0), + (y-y0)*(y-y0) , + (z-z0)*(z-z0));
 
   if (sqrt(r) < g_inputParam[RADIUS]) v[RHO] = 10.0;
 
 /* no need for potential vector, 
    since CT_VEC_POT_INIT is set to NO
 
   v[AX1] = x3*v[BX2] - x2*v[BX3];
   v[AX2] = 0.0;
   v[AX3] = 0.0;
 */
 }

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.

Definition at line 146 of file init.c.

 {
   int   i, j, k;
 
   if (side == X1_END){          /* -- select the boundary side -- */
     if (box->vpos == CENTER){   /* -- select the variable position -- */
       BOX_LOOP(box,k,j,i){      /* -- Loop over boundary zones -- */
         d->Vc[RHO][k][j][i] = 1.0;
         EXPAND(d->Vc[VX1][k][j][i] = -11.2536;  ,
                d->Vc[VX2][k][j][i] = 0.0;       ,
                d->Vc[VX3][k][j][i] = 0.0;)
         d->Vc[PRS][k][j][i] = 1.0;
         EXPAND(d->Vc[BX1][k][j][i] = 0.0;        ,
                d->Vc[BX2][k][j][i] = g_inputParam[B_PRE]; ,
                d->Vc[BX3][k][j][i] = g_inputParam[B_PRE];)
       }
     }else if (box->vpos == X2FACE){  /* -- y staggered field -- */
       #ifdef STAGGERED_MHD
        BOX_LOOP(box,k,j,i) d->Vs[BX2s][k][j][i] = g_inputParam[B_PRE];
       #endif
     }else if (box->vpos == X3FACE){  /* -- z staggered field -- */
       #ifdef STAGGERED_MHD
        BOX_LOOP(box,k,j,i) d->Vs[BX3s][k][j][i] = g_inputParam[B_PRE];
       #endif
     }
   }
 }

Functions

Detailed Description

Function Documentation