PLUTO
Macros | Functions
tc_flux.c File Reference

Compute the thermal conduction flux. More...

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

Go to the source code of this file.

Macros

#define TC_SATURATED_FLUX   YES
 
#define HYPERBOLIC_SAT_FLUX   YES
 

Functions

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

Detailed Description

Compute the thermal conduction flux.

Compute the thermal conduction flux along one row of computational zones for the HD and MHD modules according to Spitzer (1962):

\[ \vec{F}_c = \frac{q}{|\vec{F}_{\rm class}| + q}\vec{F}_{\rm class} \]

where $ \vec{F}_{\rm class} = \kappa_\parallel \nabla T$ is the classical (hydrodynamic) thermal conduction flux, $ q = F_{\rm sat} = 5\phi \rho c_{\rm iso}^3$ is the saturated flux. Temperature is, at present, simply computed as $p/\rho$.

The classical flux is purely parabolic, and it is discretized using standard finite differences. The saturated flux is included only when the macro TC_SATURATED_FLUX is enabled and it is treated in an upwind manner following the guidelines given in the appendix of Mignone et al. 2012 (see also Balsara (2008) for alternative discretization methods).

In MHD, the classical flux further splits into 2 components, along and across the magnetic field lines (see Eq. [7] of Mignone et al.).

This function also computes the inverse of the time step and return its maximum over the current sweep.

References

  • "Simulating anisotropic thermal conduction in supernova remnants - I. Numerical methods", Balsara, Tilley, Howk, MNRAS (2008) 386, 627
  • "The PLUTO Code for Adaptive Mesh Computations in Astrophysical Fluid Dynamics"
    Mignone et al, ApJS (2012) 198, 7M
Authors
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
T. Matsakos
Date
Aug 24, 2015

Definition in file tc_flux.c.

Macro Definition Documentation

#define HYPERBOLIC_SAT_FLUX   YES

When set to YES, saturated flux is computed using an upwind selection rule. When set to NO, staurated flux is treated in the same manner as the conduction flux.

Definition at line 54 of file tc_flux.c.

#define TC_SATURATED_FLUX   YES

When set to YES (default), include saturation effects when computing the conduction flux.

Definition at line 47 of file tc_flux.c.

Function Documentation

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: