PLUTO
init.c File Reference

Propagation of a conduction front. More...

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

Go to the source code of this file.

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

Propagation of a conduction front.

The problem sets the initial condition for a propagating conduction front for which an analytical solution exists (Reale 1995).

The equation

\[ \frac{\partial T}{\partial t} = a\frac{\partial}{\partial s} \left(T^n\frac{\partial T}{\partial s}\right) \]

has the solution

\[ T = T_c\left(1 - \frac{s^2}{s_f^2}\right)^{1/n} \]

where

\[ T_c = \left(\frac{Q^2}{at}\right)^{1/(n+2)} \left(\frac{n}{2(n+2)}\zeta_0^2\right)^{1/n}\,, \quad\quad\quad s_f = (aQ^nt)^{1/(n+2)}\zeta_0\,, \quad\quad\quad \zeta_0 = \left[\frac{(n+2)^{1+n}2^{1-n}}{n\pi^{n/2}}\right]^{1/(n+2)} \left[\frac{\Gamma(1/2 + 1/n)}{\Gamma(1/n)}\right]^{n/(n+2)}\,, \]

$Q$ is the integral over the whole space, and $\Gamma$ is the gamma function.

The setup is built to compare the numerical solution with the analytical one.

1D-hd-exp.jpg
Evolution of the TC front in the 1D case; numerical solution (points) and analytical (lines)

.

Configurations (EXPLICIT and STS):

  • #01-#04: 1D cartesian
  • #05-#06: 2D cartesian
  • #07-#09: 3D cartesian
  • #10-#11: 3D polar
  • #12-#13: 3D spherical
Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
July 09, 2014

References:

  • Section 4 of Reale, F. Computer Physics Communications (1995), 86, 13

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

Compute volume-integrated magnetic pressure, Maxwell and Reynolds stresses. Save them to "averages.dat"

Definition at line 139 of file init.c.

144 {
145 }
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]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 54 of file init.c.

57 {
58 #if DIMENSIONS == 1
59  double n0, k0, Q, t0, j0, xf, Tc, T0;
60  double G1_2p1_n, G1_n;
61 
62  n0 = 5./2.;
63  k0 = 4.412;
64  Q = 1.2e15;
65  t0 = 0.1;
66 
67  G1_2p1_n = 1.068628702119319;
68  G1_n = 2.218159543757688;
69 
70  j0 = pow((pow(n0+2.,1.+n0)*pow(2.,1.-n0))/(n0*pow(CONST_PI,n0/2.)),1./(n0+2.))*pow(G1_2p1_n/G1_n,n0/(n0+2.));
71  xf = pow(k0*pow(Q,n0)*t0, 1./(n0+2.))*j0/1.0e8;
72  Tc = pow(Q*Q/(k0*t0), 1./(n0+2.))*pow(n0*j0*j0/(2*(n0+2.)), 1./n0);
73  T0 = 6.05736872274111e7;
74 
75  us[RHO] = 1.0;
76  us[VX1] = 0.0;
77  us[VX2] = 0.0;
78  us[VX3] = 0.0;
79 
80  if ((1.0 - x1*x1/xf/xf) < 1.0e-20){
81  us[PRS] = us[RHO]*Tc/T0*pow(1.0e-20,0.4);
82  }else{
83  us[PRS] = us[RHO]*Tc/T0*pow( (1.0 - x1*x1/xf/xf) ,0.4);
84  }
85  us[TRC] = 0.0;
86 
87  #if PHYSICS == MHD
88  us[BX1] = 1.0;
89  us[BX2] = 0.0;
90  us[BX3] = 0.0;
91  #endif
92 
93 #else
94 
95  double n0, k0, Q, t0, j0, xf, Tc, T0;
96  double G5_2p1_n, G1p1_n, G3_2, x;
97 
98  n0 = 5./2.;
99  k0 = 4.412;
100  Q = 1.2e30;
101  t0 = 0.1;
102 
103  G5_2p1_n = 1.827355080624036;
104  G1p1_n = 0.887263817503075;
105  G3_2 = 0.886226925452758;
106 
107  j0 = pow((3.*n0+2.)/(pow(2., n0-1.)*n0*pow(CONST_PI,n0)),1./(3.*n0+2.))*pow(G5_2p1_n/G1p1_n/G3_2,n0/(3.*n0+2.));
108  xf = pow(k0*pow(Q,n0)*t0, 1./(3.*n0+2.))*j0/1.e8;
109  Tc = j0*j0*j0/xf/xf/xf*1.e-24*Q*pow(n0*j0*j0/(2.*(3.*n0+2.)), 1./n0);
110  T0 = 6.05736872274111e7;
111 
112  #if GEOMETRY == CARTESIAN || GEOMETRY == CYLINDRICAL
113  x = sqrt(D_EXPAND(x1*x1, + x2*x2, + x3*x3));
114  #elif GEOMETRY == POLAR
115  x = sqrt(D_EXPAND(x1*x1, + 0.0, + x3*x3));
116  #elif GEOMETRY == SPHERICAL
117  x = x1;
118  #endif
119 
120 
121  us[RHO] = 1.0;
122  us[VX1] = 0.0;
123  us[VX2] = 0.0;
124  us[VX3] = 0.0;
125  if ((1.0 - x*x/xf/xf) < 1.0e-20){
126  #if GEOMETRY == SPHERICAL
127  us[PRS] = us[RHO]*Tc/T0*pow(1.0e-8,0.4);
128  #else
129  us[PRS] = us[RHO]*Tc/T0*pow(1.0e-20,0.4);
130  #endif
131  }else{
132  us[PRS] = us[RHO]*Tc/T0*pow( (1.0 - x*x/xf/xf) ,0.4);
133  }
134  us[TRC] = 0.0;
135 #endif
136 }
#define VX2
Definition: mod_defs.h:29
#define RHO
Definition: mod_defs.h:19
#define TRC
Definition: pluto.h:581
double xf
Leftmost and rightmost point in the local domain.
Definition: structs.h:79
#define VX1
Definition: mod_defs.h:28
double * x
Definition: structs.h:80
#define BX3
Definition: mod_defs.h:27
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
Definition: analysis.c:27
#define BX1
Definition: mod_defs.h:25
#define VX3
Definition: mod_defs.h:30
#define CONST_PI
.
Definition: pluto.h:269
#define BX2
Definition: mod_defs.h:26
static double t0
Definition: failsafe.c:5

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]dpointer to the PLUTO data structure containing cell-centered primitive quantities (d->Vc) and staggered magnetic fields (d->Vs, when used) to be filled.
[in]boxpointer 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]sidespecifies 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]gridpointer 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]boxpointer 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]sidespecifies 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]gridpointer 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 147 of file init.c.

151 {
152  int i, j, k, nv;
153  double x1, x2, x3;
154 
155  if (side == 0) { /* -- check solution inside domain -- */
156  DOM_LOOP(k,j,i){
157  EXPAND(d->Vc[VX1][k][j][i] = 0.0; ,
158  d->Vc[VX2][k][j][i] = 0.0; ,
159  d->Vc[VX3][k][j][i] = 0.0;)
160  }
161  #if ENTROPY_SWITCH == YES
162  TOT_LOOP(k,j,i) d->flag[k][j][i] |= FLAG_ENTROPY;
163  #endif
164  }
165 
166  if (side == X1_BEG){ /* -- X1_BEG boundary -- */
167  if (box->vpos == CENTER){
168  BOX_LOOP(box,k,j,i){
169  x1 = grid[IDIR].x[i];
170  d->Vc[RHO][k][j][i] = d->Vc[RHO][k][j][2*IBEG - i - 1];
171  EXPAND(
172  d->Vc[VX1][k][j][i] = -d->Vc[VX1][k][j][2*IBEG - i - 1]; ,
173  d->Vc[VX2][k][j][i] = d->Vc[VX2][k][j][2*IBEG - i - 1]; ,
174  d->Vc[VX3][k][j][i] = d->Vc[VX3][k][j][2*IBEG - i - 1];
175  )
176  d->Vc[PRS][k][j][i] = d->Vc[PRS][k][j][2*IBEG - i - 1];
177  #if PHYSICS == MHD
178  EXPAND(
179  d->Vc[BX1][k][j][i] = d->Vc[BX1][k][j][2*IBEG - i - 1]; ,
180  d->Vc[BX2][k][j][i] = d->Vc[BX2][k][j][2*IBEG - i - 1]; ,
181  d->Vc[BX3][k][j][i] = d->Vc[BX3][k][j][2*IBEG - i - 1];
182  )
183  #endif
184  }
185  }
186  }
187 
188  if (side == X1_END){ /* -- X1_END boundary -- */
189  X1_END_LOOP(k,j,i){}
190  }
191 
192  if (side == X2_BEG){ /* -- X2_BEG boundary -- */
193  X2_BEG_LOOP(k,j,i){}
194  }
195 
196  if (side == X2_END){ /* -- X2_END boundary -- */
197  X2_END_LOOP(k,j,i){}
198  }
199 
200  if (side == X3_BEG){ /* -- X3_BEG boundary -- */
201  X3_BEG_LOOP(k,j,i){}
202  }
203 
204  if (side == X3_END) { /* -- X3_END boundary -- */
205  X3_END_LOOP(k,j,i){}
206  }
207 }
#define X3_BEG
Boundary region at X3 beg.
Definition: pluto.h:150
#define X1_BEG
Boundary region at X1 beg.
Definition: pluto.h:146
DOM_LOOP(k, j, i)
Definition: analysis.c:22
#define VX2
Definition: mod_defs.h:29
#define CENTER
Definition: pluto.h:200
int vpos
Location of the variable inside the cell.
Definition: structs.h:359
#define RHO
Definition: mod_defs.h:19
#define BOX_LOOP(B, k, j, i)
Definition: macros.h:70
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
#define X3_END_LOOP(k, j, i)
Definition: macros.h:52
#define X1_END
Boundary region at X1 end.
Definition: pluto.h:147
#define VX1
Definition: mod_defs.h:28
#define X2_BEG_LOOP(k, j, i)
Definition: macros.h:47
#define X2_END_LOOP(k, j, i)
Definition: macros.h:51
#define TOT_LOOP(k, j, i)
Definition: macros.h:44
#define FLAG_ENTROPY
Update pressure using entropy equation.
Definition: pluto.h:182
#define IDIR
Definition: pluto.h:193
#define X2_END
Boundary region at X2 end.
Definition: pluto.h:149
int j
Definition: analysis.c:2
int k
Definition: analysis.c:2
double * x
Definition: structs.h:80
#define X1_END_LOOP(k, j, i)
Definition: macros.h:50
#define BX3
Definition: mod_defs.h:27
#define X3_END
Boundary region at X3 end.
Definition: pluto.h:151
int i
Definition: analysis.c:2
#define BX1
Definition: mod_defs.h:25
#define VX3
Definition: mod_defs.h:30
#define X2_BEG
Boundary region at X2 beg.
Definition: pluto.h:148
#define BX2
Definition: mod_defs.h:26
#define X3_BEG_LOOP(k, j, i)
Definition: macros.h:48
long int IBEG
Lower grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:35

Here is the call graph for this function: