Stellar wind test problem. More...

#include "pluto.h"

Include dependency graph for init.c:

Functions
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

Stellar wind test problem.

Sets initial condition for a spherically symmetric radial wind blowing from the origin of coordinates, see [Mig14] The initial condition consists of a constant-density ambient medium with

$\rho = \rho_a \,,\qquad p = p_a\,,\qquad \vec{v} = v_{csm}\hvec{j}$

where v_csm is the velocity of the star with respect to the background. The wind is injeted using the INTERNAL_BOUNDARY where flow quantities are kept constant in time and equal to

$r^2v_r\rho = \rho_0 V_0^2r_0^2 \,,\qquad v_r = V_0\hvec{r} \,,\qquad p = \frac{c_s^2}{\Gamma}\rho^\Gamma$

These value are defined through the UserDefBoundary() function when side is equal to 0.

Dimensions are chosen so that the spherical wind shell has radius 1, density 1 and velocity 1 ( $r_0 = V_0 = \rho_0 = 1$ ).

The input parameters that control the problem dynamics are

g_inputParam[CS_WIND]: sets the sound speed in the wind region;
g_inputParam[RHO_AMB]: sets the ambient density
g_inputParam[CS_AMB]: sets the ambient sound speed;
g_inputParam[V_CSM]: sets the velocity of the star with respect to the background;

Configurations #01-05 and #07-08 work in 2D cylindrical axisymmetric coordinates while conf. #06 is 3D Cartesian. An AMR setup is available with configuration #04.

Density map at the end of computation for configuration #8

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

Date: Sept 25, 2012

References

[Mig14]: "High-order conservative reconstruction schemes for finite volume methods in cylindrical and spherical coordinates", Mignone, JCP (2014) 270, 784

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

 {
   double rho, cs, R;
 
   rho = g_inputParam[RHO_AMB];
   cs  = g_inputParam[CS_AMB];
   us[RHO] = rho;
   us[PRS] = cs*cs*rho/g_gamma;
 
   #if GEOMETRY == CARTESIAN
    us[VX1] = 0.0;         
    us[VX2] = 0.0;         
    us[VX3] = -g_inputParam[V_CSM];       
   #elif GEOMETRY == CYLINDRICAL
    us[VX1] = 0.0;         
    us[VX2] = -g_inputParam[V_CSM];          
   #elif GEOMETRY == SPHERICAL
    us[VX1] =  g_inputParam[V_CSM]*cos(x2);
    us[VX2] = -g_inputParam[V_CSM]*sin(x2);
   #endif
 
   R = sqrt(x1*x1 + x2*x2);
   us[TRC] = (R <= 1.0 ? 1.0:0.0);
   g_smallPressure = 1.e-5;
 }

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.

Definition at line 90 of file init.c.

 {
   int   i, j, k, nv;
   double *x1, *x2, *x3;
   double  r, r0, cs;
   double  Vwind  = 1.0, rho, vr;
 
   x1 = grid[IDIR].xgc;
   x2 = grid[JDIR].xgc;
   x3 = grid[KDIR].xgc;
 
   if (side == 0){
     r0 = 1.0;
     cs = g_inputParam[CS_WIND];
     TOT_LOOP(k,j,i){ 
       #if GEOMETRY == CARTESIAN
        r  = sqrt(x1[i]*x1[i] + x2[j]*x2[j] + x3[k]*x3[k]);
        if (r <= r0){  
          vr    = tanh(r/r0/0.1)*Vwind;
          rho   = Vwind*r0*r0/(vr*r*r);
          d->Vc[RHO][k][j][i] = rho;
          d->Vc[VX1][k][j][i] = Vwind*x1[i]/r;
          d->Vc[VX2][k][j][i] = Vwind*x2[j]/r;
          d->Vc[VX3][k][j][i] = Vwind*x3[k]/r;
          d->Vc[PRS][k][j][i] = cs*cs/g_gamma*pow(rho,g_gamma);
          d->flag[k][j][i]   |= FLAG_INTERNAL_BOUNDARY;
        }    
       #elif GEOMETRY == CYLINDRICAL
        r  = sqrt(x1[i]*x1[i] + x2[j]*x2[j]);
        if (r <= r0){  
          vr    = tanh(r/r0/0.1)*Vwind;
          rho   = Vwind*r0*r0/(vr*r*r);
          d->Vc[RHO][k][j][i] = rho;
          d->Vc[VX1][k][j][i] = Vwind*x1[i]/r;
          d->Vc[VX2][k][j][i] = Vwind*x2[j]/r;
          d->Vc[PRS][k][j][i] = cs*cs/g_gamma*pow(rho,g_gamma);
          d->flag[k][j][i]   |= FLAG_INTERNAL_BOUNDARY;
        }    
       #endif
     }
   }
 
   if (side == X1_BEG){  /* -- X1_BEG boundary -- */
   }
 
   if (side == X1_END){  /* -- X1_END boundary -- */
   }
 
   if (side == X2_BEG){  /* -- X2_BEG boundary -- */
    X2_BEG_LOOP(k,j,i){ }
   }
 
   if (side == X2_END){  /* -- X2_END boundary -- */
     
     cs  = g_inputParam[CS_AMB];
     rho = g_inputParam[RHO_AMB];
     X2_END_LOOP(k,j,i){
       #if GEOMETRY == CYLINDRICAL
        d->Vc[VX1][k][j][i] = 0.0;
        d->Vc[VX2][k][j][i] = -g_inputParam[V_CSM];  
        d->Vc[RHO][k][j][i] =  rho;       
        d->Vc[PRS][k][j][i] =  cs*cs*rho/g_gamma; 
       #endif
     }
   }
 
   if (side == X3_BEG){  /* -- X3_BEG boundary -- */
     X3_BEG_LOOP(k,j,i){}
   }
  
   if (side == X3_END){  /* -- X3_END boundary -- */
     
     cs  = g_inputParam[CS_AMB];
     rho = g_inputParam[RHO_AMB];
     X3_END_LOOP(k,j,i){
       #if GEOMETRY == CARTESIAN
        d->Vc[VX1][k][j][i] = 0.0;
        d->Vc[VX2][k][j][i] = 0.0;
        d->Vc[VX3][k][j][i] = -g_inputParam[V_CSM];  
        d->Vc[RHO][k][j][i] =  rho;       
        d->Vc[PRS][k][j][i] =  cs*cs*rho/g_gamma; 
       #endif
     }
   }
 }

Functions

Detailed Description

Function Documentation