PLUTO
Functions
tc.h File Reference

Thermal conduction (TC) module header file. More...

Go to the source code of this file.

Functions

void TC_Flux (double ***, const State_1D *, double **, int, int, Grid *)
 
void TC_kappa (double *, double, double, double, double *, double *, double *)
 
void GetGradient (double ***, double **, int, int, Grid *)
 

Detailed Description

Thermal conduction (TC) module header file.

Contains prototypes for the thermal conduction module.

Authors
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
T. Matsakos
Date
Sep 13, 2012

Definition in file tc.h.

Function Documentation

void GetGradient ( double ***  ,
double **  ,
int  ,
int  ,
Grid *   
)

Definition at line 4 of file tc_functions.c.

42 {
43  int i,j,k;
44  double *r, *rp;
45  double *inv_dx, *inv_dy, *inv_dz;
46  double *inv_dxi, *inv_dyi, *inv_dzi;
47  double dl1, dl2, dl3, theta, r_1, s_1;
48  double dx1, dx2, dx3;
49 
50  D_EXPAND(inv_dx = grid[IDIR].inv_dx; inv_dxi = grid[IDIR].inv_dxi; ,
51  inv_dy = grid[JDIR].inv_dx; inv_dyi = grid[JDIR].inv_dxi; ,
52  inv_dz = grid[KDIR].inv_dx; inv_dzi = grid[KDIR].inv_dxi;)
53 
54  r = grid[IDIR].x;
55  rp = grid[IDIR].xr;
56 
57  i = g_i; j = g_j; k = g_k;
58 
59  if (g_dir == IDIR) {
60 
61  #if GEOMETRY == SPHERICAL
62  theta = grid[JDIR].x[j];
63  s_1 = 1.0/sin(theta);
64  #endif
65  D_EXPAND( ,
66  dl2 = dx2 = inv_dy[j]; ,
67  dl3 = dx3 = inv_dz[k];
68  )
69  for (i = beg; i <= end; i++){
70  dl1 = inv_dxi[i];
71  #if GEOMETRY == POLAR && DIMENSIONS >= 2
72  dl2 = dx2/rp[i];
73  #elif GEOMETRY == SPHERICAL
74  D_EXPAND( ,
75  dl2 = dx2/rp[i]; ,
76  dl3 = dx3*s_1/rp[i];)
77  #endif
78  D_EXPAND(
79  gradT[i][0] = (T[k][j][i+1] - T[k][j][i])*dl1; ,
80  gradT[i][1] = 0.25*( T[k][j+1][i] + T[k][j+1][i+1]
81  - T[k][j-1][i] - T[k][j-1][i+1])*dl2; ,
82  gradT[i][2] = 0.25*( T[k+1][j][i] + T[k+1][j][i+1]
83  - T[k-1][j][i] - T[k-1][j][i+1])*dl3;
84  )
85  }
86 
87  }else if (g_dir == JDIR) {
88 
89  r_1 = 1.0/r[i];
90  D_EXPAND(
91  dl1 = dx1 = inv_dx[i]; ,
92  ,
93  dl3 = dx3 = inv_dz[k];
94  )
95  for (j = beg; j <= end; j++){
96  dl2 = inv_dyi[j];
97  #if GEOMETRY == POLAR
98  dl2 *= r_1;
99  #elif GEOMETRY == SPHERICAL
100  D_EXPAND( ,
101  dl2 *= r_1; ,
102  theta = grid[JDIR].xr[j];
103  dl3 = dx3*r_1/sin(theta);)
104  #endif
105  D_EXPAND(
106  gradT[j][0] = 0.25*( T[k][j][i+1] + T[k][j+1][i+1]
107  - T[k][j][i-1] - T[k][j+1][i-1])*dl1; ,
108  gradT[j][1] = (T[k][j+1][i] - T[k][j][i])*dl2; ,
109  gradT[j][2] = 0.25*( T[k+1][j][i] + T[k+1][j+1][i]
110  - T[k-1][j][i] - T[k-1][j+1][i])*dl3;
111  )
112  }
113 
114  }else if (g_dir == KDIR){
115 
116  dl1 = inv_dx[i];
117  dl2 = inv_dy[j];
118  #if GEOMETRY == POLAR || GEOMETRY == SPHERICAL
119  r_1 = 1.0/r[i];
120  dl2 *= r_1;
121  #endif
122  #if GEOMETRY == SPHERICAL
123  theta = grid[JDIR].x[j];
124  s_1 = 1.0/sin(theta);
125  #endif
126 
127  for (k = beg; k <= end; k++){
128  dl3 = inv_dzi[k];
129  #if GEOMETRY == SPHERICAL
130  dl3 *= r_1*s_1;
131  #endif
132  gradT[k][0] = 0.25*( T[k][j][i+1] + T[k+1][j][i+1]
133  - T[k][j][i-1] - T[k+1][j][i-1])*dl1;
134  gradT[k][1] = 0.25*( T[k][j+1][i] + T[k+1][j+1][i]
135  - T[k][j-1][i] - T[k+1][j-1][i])*dl2;
136  gradT[k][2] = (T[k+1][j][i] - T[k][j][i])*dl3;
137  }
138  }
139 }
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
GRID::end
int end
Global end index for the local array.
Definition: structs.h:116
KDIR
#define KDIR
Definition: pluto.h:195
GRID::inv_dx
double * inv_dx
Definition: structs.h:90
g_i
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
Definition: globals.h:82
IDIR
#define IDIR
Definition: pluto.h:193
g_dir
int g_dir
Specifies the current sweep or direction of integration.
Definition: globals.h:86
GRID::beg
int beg
Global start index for the local array.
Definition: structs.h:115
g_j
int g_j
x2 grid index when sweeping along the x1 or x3 direction.
Definition: globals.h:83
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
GRID::x
double * x
Definition: structs.h:80
g_k
int g_k
x3 grid index when sweeping along the x1 or x2 direction.
Definition: globals.h:84
GRID::r_1
double * r_1
Geometrical factor 1/r.
Definition: structs.h:88
D_EXPAND
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
i
int i
Definition: analysis.c:2
JDIR
#define JDIR
Definition: pluto.h:194
GRID::inv_dxi
double * inv_dxi
inverse of the distance between the center of two cells, inv_dxi = .
Definition: structs.h:91

Here is the call graph for this function:

Here is the caller graph for this function:

void TC_Flux ( double ***  T,
const State_1D *  state,
double **  dcoeff,
int  beg,
int  end,
Grid *  grid 
)

Compute the thermal conduction flux, state->par_flx.

Parameters
[in]T3D array containing the dimensionless temperature
[in,out]statepointer to a State_1D structure
[out]dcoeffthe diffusion coefficient needed for computing the time step.
[in]beginitial index of computation
[in]endfinal index of computation
[in]gridpointer to an array of Grid structures
Returns
This function has no return value.

Definition at line 57 of file tc_flux.c.

73 {
74  int i, j, k, nv;
75  double bgradT, Bmag, dTmag;
76  double Fc, Fcmag, Fsat, g1 = g_gamma - 1.0;
77  double x1, x2, x3;
78  double alpha, uL, uR, suL, suR, bn;
79  double vi[NVAR], kpar=0.0, knor=0.0, phi;
80  double bck_fld[3];
81  static double **gradT;
82 
83 /* -----------------------------------------------------------
84  1. Allocate memory, compute temperature gradient in the
85  required direction and set 2 out of the 3 coordinates.
86  ----------------------------------------------------------- */
87 
88  if (gradT == NULL) {
89  gradT = ARRAY_2D(NMAX_POINT, 3, double);
90  }
91 
92  GetGradient (T, gradT, beg, end, grid);
93  if (g_dir == JDIR || g_dir == KDIR) x1 = grid[IDIR].x[g_i];
94  if (g_dir == IDIR || g_dir == KDIR) x2 = grid[JDIR].x[g_j];
95  if (g_dir == IDIR || g_dir == JDIR) x3 = grid[KDIR].x[g_k];
96 
97 /* -----------------------------------------------
98  2. Compute Thermal Conduction Flux (tcflx).
99  ----------------------------------------------- */
100 
101  for (i = beg; i <= end; i++){
102 
103  for (nv = 0; nv < NVAR; nv++) {
104  vi[nv] = 0.5*(state->vh[i][nv] + state->vh[i+1][nv]);
105  }
106 
107  /* ------------------------------------------------
108  2a. Obtain the thermal conduction coefficients
109  along (kpar) and across (knor) the field
110  lines. This is done at the cell interface.
111  ------------------------------------------------ */
112 
113  if (g_dir == IDIR) x1 = grid[IDIR].xr[i];
114  if (g_dir == JDIR) x2 = grid[JDIR].xr[i];
115  if (g_dir == KDIR) x3 = grid[KDIR].xr[i];
116 
117 #if PHYSICS == MHD && BACKGROUND_FIELD == YES
118  BackgroundField(x1,x2,x3, bck_fld);
119  EXPAND(vi[BX1] += bck_fld[0]; ,
120  vi[BX2] += bck_fld[1]; ,
121  vi[BX3] += bck_fld[2];)
122 #endif
123 
124  TC_kappa(vi, x1, x2, x3, &kpar, &knor, &phi);
125  dTmag = D_EXPAND( gradT[i][0]*gradT[i][0],
126  + gradT[i][1]*gradT[i][1],
127  + gradT[i][2]*gradT[i][2]);
128  dTmag = sqrt(dTmag) + 1.e-12;
129 
130  /* ---------------------------------------------------
131  2b. Compute magnitude of saturated flux using
132  a Roe-like average
133  --------------------------------------------------- */
134 
135 #if TC_SATURATED_FLUX == YES
136  #if HYPERBOLIC_SAT_FLUX == YES
137 /*
138  uL = state->vL[i][PRS]/g1; suL = sqrt(uL);
139  uR = state->vR[i][PRS]/g1; suR = sqrt(uR);
140  scrh = 5.0*phi*g1*sqrt(g1/vi[RHO]);
141  csat[i] = scrh*(uR + uL + suR*suL)/(suL + suR);
142  Fsat = 0.5*scrh*(uL*suL + uR*suR);
143  if (gradT[i][g_dir] > 0.0) Fsat += 0.5*csat[i]*(uR - uL);
144  else if (gradT[i][g_dir] < 0.0) Fsat -= 0.5*csat[i]*(uR - uL);
145 
146 */
147  uL = state->vL[i][PRS];
148  uR = state->vR[i][PRS];
149  if (gradT[i][g_dir] > 0.0) Fsat = 5.0*phi/sqrt(vi[RHO])*uR*sqrt(uR);
150  else if (gradT[i][g_dir] < 0.0) Fsat = 5.0*phi/sqrt(vi[RHO])*uL*sqrt(uL);
151  else Fsat = 5.0*phi*vi[PRS]*sqrt(vi[PRS]/vi[RHO]);
152  #else
153  Fsat = 5.0*phi*vi[PRS]*sqrt(vi[PRS]/vi[RHO]);
154  #endif
155 #endif
156 
157  /* ------------------------------------------------------
158  2c. Compute the classical MHD thermal conduction
159  flux Fc and the diffusion coefficient dcoeff.
160  ------------------------------------------------------ */
161 
162 #if PHYSICS == HD
163 
164  Fc = kpar*gradT[i][g_dir]; /* -- classical thermal conduction flux -- */
165 /*
166 {
167  double x = kpar*dTmag/Fsat;
168  double flx_lim = exp(-x*x);
169  state->par_flx[i][ENG] = flx_lim*Fc
170  + (1.0 - flx_lim)*Fsat*gradT[i][g_dir]/dTmag;
171  dcoeff[i][ENG] = flx_lim*kpar/vi[RHO]*(g_gamma - 1.0);
172 }
173 */
174  #if TC_SATURATED_FLUX == YES
175  alpha = Fsat/(Fsat + kpar*dTmag);
176  #else
177  alpha = 1.0;
178  #endif
179  state->tc_flux[i][ENG] = alpha*Fc;
180  dcoeff[i][ENG] = fabs(alpha*kpar/vi[RHO])*g1;
181 
182 #elif PHYSICS == MHD
183 
184  Bmag = EXPAND(vi[BX1]*vi[BX1], + vi[BX2]*vi[BX2], + vi[BX3]*vi[BX3]);
185  Bmag = sqrt(Bmag) + 1.e-12;
186 
187  bgradT = D_EXPAND( vi[BX1]*gradT[i][0],
188  + vi[BX2]*gradT[i][1],
189  + vi[BX3]*gradT[i][2]);
190  bgradT /= Bmag;
191 
192  bn = vi[BX1 + g_dir]/Bmag; /* -- unit vector component -- */
193  Fc = kpar*bgradT*bn + knor*(gradT[i][g_dir] - bn*bgradT);
194  Fcmag = sqrt((kpar*kpar - knor*knor)*bgradT*bgradT +
195  knor*knor*dTmag*dTmag);
196 
197  #if TC_SATURATED_FLUX == YES
198  alpha = Fsat/(Fsat + Fcmag);
199  #else
200  alpha = 1.0;
201  #endif
202  state->tc_flux[i][ENG] = alpha*Fc;
203  dcoeff[i][ENG] = fabs(Fcmag/dTmag*bn*alpha/vi[RHO])*g1;
204 #endif
205  }
206 }
alpha
static double alpha
Definition: init.c:3
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
g_gamma
double g_gamma
Definition: globals.h:112
GRID::end
int end
Global end index for the local array.
Definition: structs.h:116
STATE_1D::vh
double ** vh
Primitive state at n+1/2 (only for one step method)
Definition: structs.h:162
GRID::xr
double * xr
Definition: structs.h:81
BackgroundField
void BackgroundField(double x1, double x2, double x3, double *B0)
Definition: init.c:79
RHO
#define RHO
Definition: mod_defs.h:19
STATE_1D::vR
double ** vR
Primitive variables to the right of the interface, .
Definition: structs.h:139
KDIR
#define KDIR
Definition: pluto.h:195
TC_kappa
void TC_kappa(double *, double, double, double, double *, double *, double *)
Definition: tc_kappa.c:20
g_i
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
Definition: globals.h:82
IDIR
#define IDIR
Definition: pluto.h:193
g_dir
int g_dir
Specifies the current sweep or direction of integration.
Definition: globals.h:86
GRID::beg
int beg
Global start index for the local array.
Definition: structs.h:115
g_j
int g_j
x2 grid index when sweeping along the x1 or x3 direction.
Definition: globals.h:83
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
GRID::x
double * x
Definition: structs.h:80
STATE_1D::tc_flux
double ** tc_flux
Thermal conduction flux.
Definition: structs.h:152
g_k
int g_k
x3 grid index when sweeping along the x1 or x2 direction.
Definition: globals.h:84
GetGradient
void GetGradient(double ***, double **, int, int, Grid *)
Definition: tc_functions.c:4
BX3
#define BX3
Definition: mod_defs.h:27
D_EXPAND
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
NMAX_POINT
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
i
int i
Definition: analysis.c:2
BX1
#define BX1
Definition: mod_defs.h:25
STATE_1D::vL
double ** vL
Primitive variables to the left of the interface, .
Definition: structs.h:136
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
BX2
#define BX2
Definition: mod_defs.h:26
JDIR
#define JDIR
Definition: pluto.h:194
NVAR
#define NVAR
Definition: pluto.h:609

Here is the call graph for this function:

Here is the caller graph for this function:

void TC_kappa ( double *  v,
double  x1,
double  x2,
double  x3,
double *  kpar,
double *  knor,
double *  phi 
)

Compute thermal conduction coefficients.

Parameters
[in]varray of primitive variables
[in]x1coordinate in the X1 direction
[in]x2coordinate in the X2 direction
[in]x3coordinate in the X3 direction
[out]kparpointer to the conduction coefficient $ \kappa_\parallel $ in the direction of magnetic field
[out]knorpointer to the conduction coefficient $ \kappa_\perp $ perpendicular to magnetic field
[out]phipointer to the parameter $ \phi $ controlling the magnitude of the saturated flux.

Definition at line 20 of file tc_kappa.c.

38 {
39  double mu, T, B2_cgs, nH;
40 
41  mu = 0.5;
42  T = v[PRS]/v[RHO]*mu*KELVIN;
43 
44  *kpar = 5.6e-7*T*T*sqrt(T);
45  #if PHYSICS == MHD
46  B2_cgs = EXPAND(v[BX1]*v[BX1], + v[BX2]*v[BX2], + v[BX3]*v[BX3]) + 1.e-12;
47  B2_cgs *= 4.0*CONST_PI*UNIT_DENSITY*UNIT_VELOCITY*UNIT_VELOCITY;
48  nH = v[RHO]*UNIT_DENSITY/CONST_mp;
49  *knor = 3.3e-16*nH*nH/(sqrt(T)*B2_cgs);
50  #endif
51 
52  *phi = 0.3;
53 }
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
UNIT_DENSITY
#define UNIT_DENSITY
Unit density in gr/cm^3.
Definition: pluto.h:369
RHO
#define RHO
Definition: mod_defs.h:19
KELVIN
#define KELVIN
Definition: pluto.h:401
CONST_mp
#define CONST_mp
Proton mass.
Definition: pluto.h:261
UNIT_VELOCITY
#define UNIT_VELOCITY
Unit velocity in cm/sec.
Definition: pluto.h:377
BX3
#define BX3
Definition: mod_defs.h:27
BX1
#define BX1
Definition: mod_defs.h:25
CONST_PI
#define CONST_PI
.
Definition: pluto.h:269
BX2
#define BX2
Definition: mod_defs.h:26

Here is the caller graph for this function: