PLUTO
Functions | Variables
init.c File Reference

Axisymmetric jet propagation. More...

#include "pluto.h"
Include dependency graph for init.c:

Go to the source code of this file.

Functions

void GetJetValues (double, double *)
 
static double Profile (double, double)
 
void Init (double *v, double x1, double x2, double x3)
 
void Analysis (const Data *d, Grid *grid)
 
void UserDefBoundary (const Data *d, RBox *box, int side, Grid *grid)
 

Variables

static double a = 0.8
 
static double Ta
 
static double pa
 
static double mua
 

Detailed Description

Axisymmetric jet propagation.

Description.
The jet problem is set up in axisymmetric cylindrical coordinates $ (R,z) $. We label ambient and jet values with the suffix "a" and "j", respectively. The ambient medium is at rest and is characterized by uniform density and pressure, $ \rho_a $ and $ p_a $. The beam enters from the lower-z boundary from a circular nozzle of radius $ R_j $ and carries a constant poloidal field $ B_z $ and a (radially-varying) toroidal component $ B_\phi(R) $. Flow variables are prescribed as

\[ \left\{\begin{array}{lcl} \rho(R) &=& \rho_j \\ \noalign{\medskip} v_z(R) &=& v_j \\ \noalign{\medskip} B_\phi(R) &=& \DS \left\{\begin{array}{ll} -B_m R/a & \quad{\rm for}\quad R<a \\ \noalign{\medskip} -B_m a/R & \quad{\rm for}\quad R>a \end{array}\right. \\ \noalign{\medskip} B_z(R) &=& B_{z0}\quad\mathrm{(const)} \end{array}\right. \]

Here a is the magnetization radius while $ v_R = B_R = 0 $. These profiles are similar to the ones used by Tesileanu et al. (2008). The pressure distribution is found by solving the radial momentum balance between thermal, centrifugal and magnetic forces:

\[ \frac{dp}{dR} = \frac{\rho v_\phi^2}{R} - \frac{1}{2}\left[\frac{1}{R^2}\frac{d(R^2B^2_\phi)}{dR} +\frac{dB_z^2}{dR}\right] \]

Neglecting rotation and assuming Bz to be constant the solution of radial momentum balance becomes

\[ p(R) = p_a + B_m^2\left[1-\min\left(\frac{R^2}{a^2},1\right)\right] \]

and the jet on-axis pressure increases for increasing toroidal field:

\[ p(R=0) \equiv p_j = p_a + B_m^2 \]

where $p_a$ is the ambient pressure.

Normalization.
Since the MHD equations are scale-invariant we have the freedom to specify a reference length, density and velocity. Here we choose

  • Length: jet radius $R_j=1$;
  • Density: ambient density $\rho_a=1$;
  • Velocity:
    • adiabtic setup (no cooling): ambient sound speed: $c_a = \Gamma p_a/\rho_a = 1$ (from which it follows $p_a = 1/\Gamma$).
    • radiative setup (with cooling): 1 Km/s (set by UNIT_VELOCITY). In this case, the ambient pressure is computed from the ambient temperature Ta = 2500 K .

In this way the number of parameters is reduced to 4:

  1. g_inputParam[ETA]: density contrast $ \eta = \rho_j/\rho_a \qquad\Longrightarrow\qquad \rho_j = \eta $
  2. g_inputParam[JET_VEL]: Jet velocity $ v_j $. This is also the Mach number of the beam with respect to the ambient in the adiabatic setup;
  3. g_inpurParam[SIGMA_Z]: poloidal magnetization;
  4. g_inpurParam[SIGMA_PHI]: toroidal magnetization; Magnetization are defined as follows:

    \[ \left\{\begin{array}{ll} \sigma_z &= \DS \frac{B_z^2}{2p_a} \\ \noalign{\medskip} \sigma_\phi &= \DS \frac{<B_\phi^2>}{2p_a} \end{array}\right. \qquad\Longrightarrow\qquad \left\{\begin{array}{ll} B_z^2 &= \DS 2\sigma_zp_a \\ \noalign{\medskip} B_m^2 &= \DS \frac{2\sigma_\phi}{a^2(1/2-2\log a)}p_a \end{array}\right. \]

    Here the average value of $B^2_\phi$ is simply found as

    \[ <B^2_\phi> = \frac{\int_0^1 B^2_\phi R\,dR}{\int_0^1 R\,dR} = B_m^2a^2\left(\frac{1}{2} - 2\log a\right) \]

The following MAPLE code can be used to verify the solution:

restart; # Solve radial equilibrium balance equation
assume (a < 1, a > 0);
B[phi] := -B[m]*R/a;
ode := diff(p(R),R) = -1/(2*R^2)*diff(R^2*B[phi]^2,R);
dsolve({ode,p(a)=pa}); # only for R < a

When cooling is enabled, two additional parameters controlling the amplitude and frequency of perturbation can be used:

  1. g_inputParam[PERT_AMPLITUDE]: perturbation amplitude (in units of jet velocity);
  2. g_inputParam[PERT_PERIOD]: perturbation period (in years).

The jet problem has been tested using a variety of configurations, namely:

Conf.GEOMETRY DIMT. STEPPINGRECONSTRUCTIONdivBEOS COOLING
#01 CYLINDRICAL2 RK2 LINEAR CT IDEAL NO
#02 CYLINDRICAL2 HANCOCK LINEAR CT IDEAL NO
#03 CYLINDRICAL2 ChTr PARABOLIC CT IDEAL NO
#04 CYLINDRICAL2.5RK2 LINEAR NONEIDEAL NO
#05 CYLINDRICAL2.5RK2 LINEAR CT IDEAL NO
#06 CYLINDRICAL2.5RK2 LINEAR CT PVTE_LAWNO
#07 CYLINDRICAL2.5ChTr PARABOLIC NONEIDEAL SNEq
#08 CYLINDRICAL2.5ChTr PARABOLIC NONEIDEAL MINEq
#09 CYLINDRICAL2.5ChTr PARABOLIC NONEIDEAL H2_COOL

Here 2.5 is a short-cut for to 2 dimensions and 3 components.

The following image show the (log) density map at the end of simulation for setup #01.

mhd_jet.01.jpg
Density map at the end of the computation for configuration #01
Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
Sept 12, 2014

References

  • "Simulating radiative astrophysical flows with the PLUTO code: a non-equilibrium, multi-species cooling function" Tesileanu, Mignone and Massaglia, A&A (2008) 488,429.

Definition in file init.c.

Function Documentation

void Analysis ( const Data *  d,
Grid *  grid 
)

Perform runtime data analysis.

Parameters
[in]dthe PLUTO Data structure
[in]gridpointer to array of Grid structures

Definition at line 187 of file init.c.

192 {
193 
194 }
void GetJetValues ( double  R,
double *  vj 
)

Define jet value as functions of the (cylindrical) radial coordinate. For a periodic jet configuration, density and verical velocity smoothly join the ambient values. This has no effect on equilibrium as long as rotation profile depends on the chosen density profile.

Definition at line 264 of file init.c.

273 {
274  int nv;
275  static int first_call = 1;
276  static double veq[NVAR];
277  double Bz, Bm = 0.0, x;
278 
279  if (fabs(R) < 1.e-9) R = 1.e-9;
280  x = R/a;
281 
282  vj[RHO] = 1.0 + (g_inputParam[ETA] - 1.0);
283  Bz = sqrt(2.0*g_inputParam[SIGMA_Z]*pa);
284  Bm = sqrt(2.0*g_inputParam[SIGMA_PHI]*pa/(a*a*(0.5 - 2.0*log(a))));
285 
286  EXPAND( vj[BX1] = 0.0; ,
287  vj[BX2] = Bz; ,
288  vj[BX3] = -Bm*(x < 1.0 ? x:1.0/x);)
289  vj[AX1] = vj[AX2] = 0.0;
290  vj[AX3] = 0.5*R*Bz;
291 
292  EXPAND( vj[VX1] = 0.0; ,
293  vj[VX2] = g_inputParam[JET_VEL]; ,
294  vj[VX3] = 0.0; )
295 
296  vj[PRS] = pa + Bm*Bm*(1.0 - MIN(x*x,1.0));
297  vj[TRC] = 1.0;
298 
299  #if COOLING != NO
300  static double Tj,muj;
301  #if COOLING == H2_COOL
302  Tj = Ta*(1.0 + Bm*Bm/pa)/g_inputParam[ETA]; /* Guess Tj assuming muj = mua */
303  if (first_call){
304  CompEquil(200.0*g_inputParam[ETA], Tj, veq);
305  muj = MeanMolecularWeight(veq);
306  Tj *= muj/mua; /* Improve guess */
307  }
308  #else
309  if (first_call) CompEquil(200.0*g_inputParam[ETA], Ta, veq);
310  #endif
311  for (nv = NFLX; nv < NFLX+NIONS; nv++) vj[nv] = veq[nv];
312  #endif
313  first_call = 0;
314 
315 }
a
static double a
Definition: init.c:135
VX2
#define VX2
Definition: mod_defs.h:29
RHO
#define RHO
Definition: mod_defs.h:19
CompEquil
void CompEquil(double n, double T, double *v)
Definition: comp_equil.c:44
AX2
#define AX2
Definition: mod_defs.h:86
JET_VEL
#define JET_VEL
Definition: definitions_01.h:27
TRC
#define TRC
Definition: pluto.h:581
pa
static double pa
Definition: init.c:136
VX1
#define VX1
Definition: mod_defs.h:28
Ta
static double Ta
Definition: init.c:136
MIN
#define MIN(a, b)
Definition: macros.h:104
SIGMA_PHI
#define SIGMA_PHI
Definition: definitions_01.h:29
NFLX
#define NFLX
Definition: mod_defs.h:32
NIONS
#define NIONS
Definition: cooling.h:28
AX3
#define AX3
Definition: mod_defs.h:87
g_inputParam
double g_inputParam[32]
Array containing the user-defined parameters.
Definition: globals.h:131
GRID::x
double * x
Definition: structs.h:80
BX3
#define BX3
Definition: mod_defs.h:27
MeanMolecularWeight
double MeanMolecularWeight(double *V)
Definition: radiat.modified_by_Bhargav.c:178
ETA
#define ETA
Definition: definitions_01.h:23
AX1
#define AX1
Definition: mod_defs.h:85
BX1
#define BX1
Definition: mod_defs.h:25
mua
static double mua
Definition: init.c:136
VX3
#define VX3
Definition: mod_defs.h:30
BX2
#define BX2
Definition: mod_defs.h:26
SIGMA_Z
#define SIGMA_Z
Definition: definitions_01.h:28
NVAR
#define NVAR
Definition: pluto.h:609

Here is the call graph for this function:

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]va pointer to a vector of primitive variables
[in]x1coordinate point in the 1st dimension
[in]x2coordinate point in the 2nd dimension
[in]x3coordinate 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 139 of file init.c.

143 {
144  int nv;
145  static int first_call = 1;
146  double nH = 200.0;
147  static double veq[NVAR];
148 
149  v[RHO] = 1.0;
150  v[VX1] = v[VX2] = v[VX3] = 0.0;
151 #if COOLING == NO
152  v[PRS] = pa = 0.6;
153  /* v[PRS] = Pressure(v, Ta); */
154 #else
155  #if COOLING == H2_COOL
156  Ta = 100.0; /* Ambient temperature for molecular cooling */
157  #else
158  Ta = 2500.0; /* Ambient temperature for atomic cooling */
159  #endif
160  if (first_call) CompEquil(nH, Ta, veq);
161  for (nv = NFLX; nv < NFLX + NIONS; nv++) v[nv] = veq[nv];
162  mua = MeanMolecularWeight(v);
163  v[PRS] = pa = Ta*v[RHO]/(KELVIN*mua);
164 #endif
165 
166  EXPAND(v[BX1] = 0.0; ,
167  v[BX2] = sqrt(2.0*g_inputParam[SIGMA_Z]*pa); ,
168  v[BX3] = 0.0;)
169 
170  v[AX1] = v[AX2] = 0.0;
171  v[AX3] = 0.5*x1*v[BX2];
172 
173  v[TRC] = 0.0;
174 
175 #if COOLING == H2_COOL
176  g_minCoolingTemp = 10.0;
177  g_maxCoolingRate = 0.1;
178 #else
179  g_minCoolingTemp = 2500.0;
180  g_maxCoolingRate = 0.2;
181 #endif
182  g_smallPressure = 1.e-3;
183 
184  first_call = 0;
185 }
VX2
#define VX2
Definition: mod_defs.h:29
RHO
#define RHO
Definition: mod_defs.h:19
CompEquil
void CompEquil(double n, double T, double *v)
Definition: comp_equil.c:44
AX2
#define AX2
Definition: mod_defs.h:86
TRC
#define TRC
Definition: pluto.h:581
g_smallPressure
double g_smallPressure
Small value for pressure fix.
Definition: globals.h:110
KELVIN
#define KELVIN
Definition: pluto.h:401
pa
static double pa
Definition: init.c:136
VX1
#define VX1
Definition: mod_defs.h:28
Ta
static double Ta
Definition: init.c:136
NFLX
#define NFLX
Definition: mod_defs.h:32
NIONS
#define NIONS
Definition: cooling.h:28
AX3
#define AX3
Definition: mod_defs.h:87
g_inputParam
double g_inputParam[32]
Array containing the user-defined parameters.
Definition: globals.h:131
BX3
#define BX3
Definition: mod_defs.h:27
g_minCoolingTemp
double g_minCoolingTemp
The minimum temperature (in K) below which cooling is suppressed.
Definition: globals.h:106
MeanMolecularWeight
double MeanMolecularWeight(double *V)
Definition: radiat.modified_by_Bhargav.c:178
AX1
#define AX1
Definition: mod_defs.h:85
BX1
#define BX1
Definition: mod_defs.h:25
mua
static double mua
Definition: init.c:136
VX3
#define VX3
Definition: mod_defs.h:30
BX2
#define BX2
Definition: mod_defs.h:26
g_maxCoolingRate
double g_maxCoolingRate
The maximum fractional variation due to cooling from one step to the next.
Definition: globals.h:104
SIGMA_Z
#define SIGMA_Z
Definition: definitions_01.h:28
NVAR
#define NVAR
Definition: pluto.h:609

Here is the call graph for this function:

double Profile ( double  R,
double  nv 
)
static

Definition at line 318 of file init.c.

323 {
324  double R4 = R*R*R*R, R8 = R4*R4;
325 
326  #if PHYSICS == MHD && COMPONENTS == 3 /* Start with a smoother profile */
327  if (g_time < 0.1) return 1.0/cosh(R4); /* with toroidal magnetic fields */
328  #endif
329  return 1.0/cosh(R8);
330 }
g_time
double g_time
The current integration time.
Definition: globals.h:117

Here is the caller graph for this function:

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

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

Definition at line 196 of file init.c.

208 {
209  int i, j, k, nv;
210  double *x1, *x2, *x3;
211  double r, vjet[256], vout[NVAR], q;
212  double t, omega; /* pulsed jet frequency */
213 
214  x1 = grid[IDIR].xgc;
215  x2 = grid[JDIR].xgc;
216  x3 = grid[KDIR].xgc;
217 
218  if (side == 0) { /* -- check solution inside domain -- */
219  DOM_LOOP(k,j,i){}
220  }
221 
222  if (side == X2_BEG){ /* -- X2_BEG boundary -- */
223  if (box->vpos == CENTER){ /* -- cell-centered boundary conditions -- */
224  BOX_LOOP(box,k,j,i){
225  GetJetValues (x1[i], vjet);
226 
227  /* -- add pulsed perturbation -- */
228 
229  #if COOLING != NO
230  t = g_time*UNIT_LENGTH/UNIT_VELOCITY; /* time in seconds */
231  omega = 2.0*CONST_PI/(86400.*365.*g_inputParam[PERT_PERIOD]);
232  vjet[VX2] *= 1.0 + g_inputParam[PERT_AMPLITUDE]*sin(omega*t);
233  #endif
234 
235  /* -- copy and reflect ambient medium values -- */
236 
237  VAR_LOOP(nv) vout[nv] = d->Vc[nv][k][2*JBEG - j - 1][i];
238  vout[VX2] *= -1.0;
239  EXPAND(vout[BX1] *= -1.0; ,
240  ; ,
241  vout[BX3] *= -1.0;)
242 
243  VAR_LOOP(nv){
244  d->Vc[nv][k][j][i] = vout[nv]-(vout[nv]-vjet[nv])*Profile(x1[i], nv);
245  }
246  }
247 
248  }else if (box->vpos == X1FACE){ /* -- staggered fields -- */
249  #ifdef STAGGERED_MHD
250  x1 = grid[IDIR].A;
251  BOX_LOOP(box,k,j,i){
252  vout[BX1] = -d->Vs[BX1s][k][2*JBEG - j - 1][i];
253  d->Vs[BX1s][k][j][i] = vout[BX1]
254  - (vout[BX1] - vjet[BX1])*Profile(x1[i], BX1);
255  }
256  #endif
257  }else if (box->vpos == X3FACE){
258 
259  }
260  }
261 }
PERT_PERIOD
#define PERT_PERIOD
Definition: definitions_07.h:31
DATA::Vs
double **** Vs
The main four-index data array used for face-centered staggered magnetic fields.
Definition: structs.h:43
DOM_LOOP
DOM_LOOP(k, j, i)
Definition: analysis.c:22
VX2
#define VX2
Definition: mod_defs.h:29
PERT_AMPLITUDE
#define PERT_AMPLITUDE
Definition: definitions_07.h:30
CENTER
#define CENTER
Definition: pluto.h:200
RBOX::vpos
int vpos
Location of the variable inside the cell.
Definition: structs.h:359
X3FACE
#define X3FACE
Definition: pluto.h:203
BOX_LOOP
#define BOX_LOOP(B, k, j, i)
Definition: macros.h:70
DATA::Vc
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
KDIR
#define KDIR
Definition: pluto.h:195
X1FACE
#define X1FACE
Definition: pluto.h:201
UNIT_VELOCITY
#define UNIT_VELOCITY
Unit velocity in cm/sec.
Definition: pluto.h:377
VAR_LOOP
#define VAR_LOOP(n)
Definition: macros.h:226
IDIR
#define IDIR
Definition: pluto.h:193
UNIT_LENGTH
#define UNIT_LENGTH
Unit Length in cm.
Definition: pluto.h:373
g_inputParam
double g_inputParam[32]
Array containing the user-defined parameters.
Definition: globals.h:131
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
BX1s
#define BX1s
Definition: ct.h:27
GRID::xgc
double * xgc
Cell volumetric centroid (!= x when geometry != CARTESIAN).
Definition: structs.h:84
BX3
#define BX3
Definition: mod_defs.h:27
Profile
static double Profile(double, double)
Definition: init.c:318
i
int i
Definition: analysis.c:2
GetJetValues
void GetJetValues(double, double *)
Definition: init.c:264
g_time
double g_time
The current integration time.
Definition: globals.h:117
BX1
#define BX1
Definition: mod_defs.h:25
CONST_PI
#define CONST_PI
.
Definition: pluto.h:269
X2_BEG
#define X2_BEG
Boundary region at X2 beg.
Definition: pluto.h:148
JDIR
#define JDIR
Definition: pluto.h:194
JBEG
long int JBEG
Lower grid index of the computational domain in the the X2 direction for the local processor...
Definition: globals.h:39
NVAR
#define NVAR
Definition: pluto.h:609
q
static Runtime q
Definition: runtime_setup.c:490
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87

Here is the call graph for this function:

Variable Documentation

double a = 0.8
static

Definition at line 135 of file init.c.

double mua
static

Definition at line 136 of file init.c.

double pa
static

Definition at line 136 of file init.c.

double Ta
static

Definition at line 136 of file init.c.