-
Notifications
You must be signed in to change notification settings - Fork 2
/
tc_kappa.c
107 lines (94 loc) · 3.71 KB
/
tc_kappa.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
#include "pluto.h"
#include "gamma_transp.h"
#include "current_table.h"
#include "transport_tables.h"
#include "capillary_wall.h"
#define KAPPAMAX 1e7
#define KAPPA_LOW 1e3
#define RHO_LOW (2.5e-9)
void TC_kappa(double *v, double x1, double x2, double x3,
double *kpar, double *knor, double *phi)
{
#if KAPPA_TABLE
static int tc_tab_not_done = 1;
#else
double mu=0.0, z=0.0;
#endif
double T=0.0;
double k=0.0;
// double unit_Mfield;
if (g_inputParam[KAPPA_GAU] > 0.0) {
// Fixed value from pluto.ini
*kpar = g_inputParam[KAPPA_GAU];
*knor = g_inputParam[KAPPA_GAU];
} else {
if (GetPV_Temperature(v, &(T) )!=0) {
#if WARN_ERR_COMP_TEMP
print1("\nTC_kappa:[Ema]Err.comp.temp");
#endif
}
#if KAPPA_TABLE
if (tc_tab_not_done) {
MakeThermConductivityTable();
tc_tab_not_done = 0;
}
if (GetThermConductivityFromTable(v[RHO]*UNIT_DENSITY, T, &k)!=0) {
print1("[TC_kappa] Error getting kappa from table\n");
print1("x1=%g, x2=%g, x3=%g\n", x1, x2, x3);
print1("rho=%g, prs=%g, T=%g",
v[RHO]*UNIT_DENSITY,
v[PRS]*UNIT_DENSITY*UNIT_VELOCITY*UNIT_VELOCITY,
T);
QUIT_PLUTO(1);
}
#else
GetMu(T, v[RHO], &mu);
z = fmax(1/mu - 1, IONIZMIN);
// unit_Mfield = COMPUTE_UNIT_MFIELD(UNIT_VELOCITY, UNIT_DENSITY);
// k = thermCond_norm(z, v[RHO]*UNIT_DENSITY, T*CONST_kB, 1, v[iBPHI]*unit_Mfield);
k = thermCond_norm_DD(z, v[RHO]*UNIT_DENSITY, T*CONST_kB);
// I increase the value of k to take into account the conductivity of a hydrogen gas at 2000K
// (according to Timrot value of H2 conductivity, cited in Mehl et al., "Ab initio transport...",2010 ),
// anyway, I am not sure this correction is ok, especially because I don't know what is the "low-density
// limit" where this value is applicable. Moreover, I think it might be that I should also add the
// conductivity of hydrogen ions (negligible at high T).
k = k + 8e4;
#endif
#ifdef KAPPAMAX
if (k > KAPPAMAX) {
k = KAPPAMAX;
}
#endif
#ifdef CONE_LOW_TCKAPPA
// Experimental: the smaller rho is, the smaller k is, with direct proportionality,
// k=alpha*rho (if rho=RHO_LOW, then k=KAPPA_LOW)
if (IsOutCone(CONE_LOW_TCKAPPA, x1, x2) && v[RHO]*UNIT_DENSITY < RHO_LOW)
k = v[RHO]*UNIT_DENSITY*KAPPA_LOW/RHO_LOW;
#endif
*knor = k;
*kpar = k; //This should be useless
// simplified formula present in the documentation
//*kpar = 5.6e-7*T*T*sqT;
//*knor = 5.6e-7*T*T*sqT;
}
/**************************************************
[Ema] Adimensionalization. It is slightly modified from
what is adviced in PLUTO's userguide:
I do not normalize by mu, since it is not constant over the domain,
I belive that it is not correct in general to normalize by mu.
Normalizing by mu could be a good idea when the EOS is IDEAL and no
change in mu over space/time is possible [at the time of writing it is
indeed necessary when EOS==IDEAL, due to the way temperature is computed in
the STS and EXPLICIT algorithms (T=p/rho). If one decides not to normalize by
mu then he should change the computation of T for the ideal case with
T=mu*p/rho]. For using ADI scheme I implemented, one must not normalize by mu.
*************************************************/
*kpar /= UNIT_KAPPA;
*knor /= UNIT_KAPPA;
/***************************************************/
*phi = 1; //[Ema] this should be useless with the saturation of the thermal conduction off
if (*kpar>1e15 || *knor>1e15) {
print1("\nDid you take the wrong convention for kappa?\n");
QUIT_PLUTO(1);
}
}