Disk-Planet interaction problem. More...

#include "pluto.h"

Include dependency graph for init.c:

Macros
#define	MIN_DENSITY 1e-8

#define	g_OmegaZ 0.0

Functions
static void	NormalizeDensity (const Data d, Grid g)

void	Init (double *us, 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

Disk-Planet interaction problem.

Simulate the interaction of a planet embedded in a disk as described in section 3.4 of Mignone et al., A&A (2012) 545, A152. This test is a nice benchmark for the FARGO module and the ROTATING_FRAME switch. For testing-purposes no viscosity is used here. The initial condition consists of a locally isothermal configuration with temperature profile $\propto T^{-1}$ yielding a disk vertical height to radius of 0.05. The gravitational potential due to the presence of the star and the planet is defined in BodyForcePotential() function.

The conventions used throught the implementation are the following:

r = spherical radius
R = cylindrical radius
z = cylindrical height
th = meridional angle

The test can be carried out in polar (2D or 3D) or spherical (3D) coordinates and the following parameters determine the initial configuration:

g_inputParam[Mstar]: controls the star mass (in solar masses)
g_inputParam[Mdisk]: controls the disk mass (in solar masses)
g_inputParam[Mplanet]: sets the planet mass (in earth masses)
g_inputParam[Viscosity]: sets the amount of viscosity

Computation can be carried in the rotating or in the observer's frame of reference (ROTATING_FRAME to YES or NO, respectively). In particular:

Configurations #01 and #02 are in 2D polar coordinates without and with FARGO, in the rotating frame.
Configurations #03, #04 and #05 are in spherical 3D coordinates with and without FARGO in the rotating frame
Configurations #06 and #07 are in 2D polar coordinates but in the observer's frame
Configuration #08 employs static AMR (grid levels are spatial dependent but not dependent on time) in the rotating frame.

Density map for configuration #08 using AMR

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

Date: Aug 16, 2012

References:

"A Conservative orbital advection scheme for simulations of magnetized shear flows with the PLUTO Code" Mignone et al, A&A (2012)

Definition in file init.c.

Macro Definition Documentation

#define g_OmegaZ 0.0

Definition at line 64 of file init.c.

#define MIN_DENSITY 1e-8

Definition at line 60 of 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 128 of file init.c.

 {
 
 }

void Init	(	double *	us,
		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 68 of file init.c.

 {
   double r, th, R, z, H, OmegaK, cs;
   double scrh;
 
   #if EOS == IDEAL
    g_gamma = 1.01;
   #endif
 
   #if ROTATING_FRAME == YES
    g_OmegaZ  = sqrt(1.0 + g_inputParam[Mplanet]/g_inputParam[Mstar]*CONST_Mearth/CONST_Msun);
    g_OmegaZ *= 2.0*CONST_PI;
   #endif
   
   #if GEOMETRY == POLAR
    R  = x1;
    #if DIMENSIONS == 2
     z  = 0.0;
     r  = R;
     th = 0.5*CONST_PI;
    #else
     z  = x3;
     r  = sqrt(R*R + z*z);
     th = atan2(R,z);
    #endif
   #elif GEOMETRY == SPHERICAL
    r  = x1;
    th = x2;
    R  = r*sin(th);
    z  = r*cos(th);
   #endif
   
   H      = 0.05*R;
   OmegaK = 2.0*CONST_PI/(R*sqrt(R));
   cs     = H*OmegaK;
   
   scrh   = (0.5*CONST_PI - th)*r/H;
   us[RHO] = 1.0/(R*sqrt(R))*exp(-0.5*scrh*scrh);
   us[VX1] = us[VX2] = us[VX3] = 0.0;
 
   us[iVPHI] = R*(OmegaK - g_OmegaZ);
   #if EOS == IDEAL
    us[PRS] = us[RHO]*cs*cs;
   #elif EOS == ISOTHERMAL
 //   g_isoSoundSpeed = cs;
    g_isoSoundSpeed = CONST_PI*0.1;
   #endif
 
 #if DUST == YES
   us[RHO_D] = 1.e-4;
   us[VX1_D] = us[VX2_D] = us[VX3_D] = 0.0;
   us[VX2_D] = us[iVPHI];
 #endif
 }

void NormalizeDensity	(	const Data *	d,
		Grid *	g
	)

static

Definition at line 199 of file init.c.

 {
   int   i, j, k;
   double mc;
         
   mc  = 0.5*g_inputParam[Mdisk]*CONST_Msun;
   mc /= UNIT_DENSITY*UNIT_LENGTH*UNIT_LENGTH*UNIT_LENGTH;
   DOM_LOOP(k,j,i){
     d->Vc[RHO][k][j][i] *= mc;
     #if EOS == IDEAL
      d->Vc[PRS][k][j][i] *= mc;
     #endif
   }
 }

Here is the call graph for this function:

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

Definition at line 137 of file init.c.

 {
   int   i, j, k, nv;
   double *x1, *x2, *x3, R, OmegaK, v[256];
   static int do_once = 1;
   
   x1 = grid[IDIR].x;
   x2 = grid[JDIR].x;
   x3 = grid[KDIR].x;
 
   #if DIMENSIONS == 3
   if (side == 0){
     if (do_once){
       NormalizeDensity(d, grid);
       do_once = 0;
     }
   }
   #endif
 
   if (side == X1_BEG){
     X1_BEG_LOOP(k,j,i){
       NVAR_LOOP(nv) d->Vc[nv][k][j][i] = d->Vc[nv][k][j][2*IBEG - i - 1];
       d->Vc[VX1][k][j][i] *= -1.0;
       #if GEOMETRY == POLAR
        R = x1[i];
       #elif GEOMETRY == SPHERICAL
        R = x1[i]*sin(x2[j]);
       #endif
       OmegaK = 2.0*CONST_PI/(R*sqrt(R));
       d->Vc[iVPHI][k][j][i] = R*(OmegaK - g_OmegaZ);
 #if DUST == YES      
 //      NDUST_LOOP(nv) d->Vc[nv][k][j][i] = 0.0;
       d->Vc[VX2_D][k][j][i] = d->Vc[iVPHI][k][j][i];
 #endif      
     }
   }
 
   if (side == X1_END){
     X1_END_LOOP(k,j,i){
       NVAR_LOOP(nv)  d->Vc[nv][k][j][i] = d->Vc[nv][k][j][IEND];
       #if GEOMETRY == POLAR
        R = x1[i];
 //       d->Vc[iVR][k][j][i] = 0.0;
       #elif GEOMETRY == SPHERICAL
        R = x1[i]*sin(x2[j]);
        d->Vc[iVR][k][j][i]  = 0.0;
        d->Vc[iVTH][k][j][i] = 0.0;
       #endif
       OmegaK = 2.0*CONST_PI/(R*sqrt(R));
       d->Vc[iVPHI][k][j][i] = R*(OmegaK - g_OmegaZ);
 #if DUST == YES      
       d->Vc[VX2_D][k][j][i] = d->Vc[iVPHI][k][j][i];
 #endif      
  
     }
   }
 }

Here is the call graph for this function:

Macros

Functions

Detailed Description

Macro Definition Documentation

Function Documentation