MHD Rotor test problem. More...

#include "pluto.h"

Include dependency graph for init.c:

Functions
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)

Detailed Description

MHD Rotor test problem.

The rotor problem consists of a rapidly spinning cylinder embedded in a static background medium with uniform density and pressure and a constant magnetic field along the x direction $B_x = 5/\sqrt{4\pi}$ . The cylinder rotates uniformly with constant angular velocity $\omega$ and has larger density:

$(\rho, v_\phi) = \left\{\begin{array}{ll} (10, \omega r) & \qquad {\rm for}\quad r < r_0 \\ \noalign{\medskip} (1+9f,f\omega r_0) & \qquad {\rm for}\quad r_0 \le r \le r_1 \\ \noalign{\medskip} (1,0) & \qquad {\rm otherwise} \end{array}\right.$

Here $ r_0=0.1$ and $ = (r_1-r)/(r_1-r_0)$ is a taper function. The ideal equation of state with $\Gamma = 1.4$ is used. As the disk rotates, strong torsional Alfven waves form and propagate outward carrying angular momentum from the disk to the ambient.

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

Conf.	GEOMETRY	divB	BCK_FIELD	AMR
#01	CARTESIAN	CT	NO	NO
#02	POLAR	CT	NO	NO
#03	POLAR	8W	NO	NO
#04	POLAR	CT	YES	NO
#05	POLAR	GLM	NO	NO
#06	CARTESIAN	GLM	NO	NO
#07	CARTESIAN	GLM	NO	YES
#08	CARTESIAN	8W	NO	YES
#09	POLAR	CT	YES	NO
#10	POLAR	CT	YES	NO
#11	POLAR	GLM	YES	YES

A snapshot of the solution using static and AMR grid is given below.

Density map (in log scale) at t = 0.15 in Cartesian coordinates using 400 x 400 grid zones (configuration #01). Field lines are super-imposed.

Density map (in log scale) at t = 0.15 in Cartesian coordinates using a base grid of 64x64 and 4 levels of refinement (conf. #08).

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

Date: July 08, 2014

Reference:

"On the divergence-free condition in Godunov-type schemes for ideal MHD: the upwind constrained transport method" P. Londrillo, L. Del Zanna, JCP (2004), 195, 17
"PLUTO: A numerical code for computational astrophysics" Mignone et al., ApJS (2007) 170, 228 (see Sect. 5.4)

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

Definition at line 132 of file init.c.

 {
 
 }

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

Definition at line 143 of file init.c.

 {
   double Bx;
   Bx = 5.0/sqrt(4.0*CONST_PI);
 
   #if GEOMETRY == CARTESIAN
    B0[0] = Bx;
    B0[1] = 0.0;
    B0[2] = 0.0;
   #elif GEOMETRY == POLAR
    B0[0] =   Bx*cos(x2);
    B0[1] = - Bx*sin(x2);
    B0[2] = 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 63 of file init.c.

 {
   double r, r0, r1, Bx, f, omega;
 
   g_gamma   = 1.4;
   omega = 20.0;
   Bx    = 5.0/sqrt(4.0*CONST_PI);
   r0    = 0.1;
   r1    = 0.115;
 
   #if GEOMETRY == CARTESIAN
    r = sqrt(x1*x1 + x2*x2);
 
    v[PRS] = 1.0;
    v[BX1] = Bx;
    v[BX2] = 0.0;
 
    f = (r1 - r)/(r1 - r0);
    if (r <= r0) {
      v[RHO] = 10.0;
      v[VX1] = -omega*x2;
      v[VX2] =  omega*x1;
    }else if (r < r1) {
      v[RHO] = 1.0 + 9.0*f;
      v[VX1] = -f*omega*x2*r0/r;
      v[VX2] =  f*omega*x1*r0/r;
    } else {
      v[RHO] = 1.0;
      v[VX1] = 0.0;
      v[VX2] = 0.0;
    }
 
    v[AX1] = v[AX2] = 0.0;
    v[AX3] = Bx*x2;
   #elif GEOMETRY == POLAR
    r = x1;
 
    v[BX1] =  Bx*cos(x2);
    v[BX2] = -Bx*sin(x2);
    v[VX1] = 0.0;
    v[PRS] = 1.0;
 
   f = (r1 - r)/(r1 - r0);
    if (r <= r0) {
      v[RHO] = 10.0;
      v[VX2] = omega*r;
    }else if (r < r1) {
      v[RHO] = 1.0 + 9.0*f;
      v[VX2] = f*omega*r0;
    } else {
      v[RHO] = 1.0;
      v[VX2] = 0.0;
    }
   
    v[AX1] = v[AX2] = 0.0;
    v[AX3] = - r*sin(x2)*Bx;
   #endif
 
   #if BACKGROUND_FIELD == YES  /* With background splitting, the initial 
                                   condition should contain deviation only */
    v[BX1] = v[BX2] = v[BX3] =
    v[AX1] = v[AX2] = v[AX3] = 0.0;
   #endif
 }

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

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).

Definition at line 164 of file init.c.

 {
   int   i, j, k, nv;
   double  *r, slp;
 
   if (side == X1_BEG){
     r = grid[IDIR].x;
     if (box->vpos == CENTER) {
       BOX_LOOP(box,k,j,i){
         slp = r[i]/r[IBEG];
         d->Vc[RHO][k][j][i] = d->Vc[RHO][k][j][IBEG];
         d->Vc[VX1][k][j][i] = slp*d->Vc[VX1][k][j][IBEG];
         d->Vc[VX2][k][j][i] = slp*d->Vc[VX2][k][j][IBEG]; 
         d->Vc[PRS][k][j][i] = d->Vc[PRS][k][j][IBEG];
         d->Vc[BX1][k][j][i] = d->Vc[BX1][k][j][IBEG];
         d->Vc[BX2][k][j][i] = d->Vc[BX2][k][j][IBEG];
       }
     }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
     }
   }
 }

Functions

Detailed Description

Function Documentation