PLUTO
Macros | Functions | Variables
thermal_eos.c File Reference

Thermodynamic relations for the thermal EoS, p=nkT and T=p/(nK) . More...

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

Go to the source code of this file.

Macros

#define PV_TEMPERATURE_TABLE_NX   768
 
#define PV_TEMPERATURE_TABLE_NY   768
 

Functions

static double TFunc (double T, void *par)
 
void MakePV_TemperatureTable ()
 
double Pressure (double *v, double T)
 
int GetPV_Temperature (double *v, double *T)
 

Variables

static Table2D Ttab
 A 2D table containing pre-computed values of temperature stored at equally spaced node values of Log(p/rho) and Log(rho) More...
 

Detailed Description

Thermodynamic relations for the thermal EoS, p=nkT and T=p/(nK) .

Collect miscellaneous functions involving computations with the thermal equation of state,

\[ p = nk_BT \quad {\rm(in\; cgs)} \qquad \Longrightarrow\qquad p = \frac{\rho T}{K\mu(\vec{X})} \quad {\rm(in\; code\quad units)} \]

where, in the second expression, rho and p are density and pressure (in code units), T is the temperature (in Kelvin), K is the KELVIN dimensionalization constant, mu(X) is the mean molecular weight and X is an array containing the ionization fractions of different elements. The previous relation is typically used to compute

  • p=p(T,rho,X) by calling Pressure();
  • T=T(p,rho,X) by calling GetPV_Temperature(). This relation can be either:
    • Explicit: ion fractions are evolved independently using a chemical reaction network (e.g. when H2_COOL or SNEq is enabled). Then one has $ T=(p/\rho) K\mu(\vec{X})$.
    • Implicit: ions are not evolved explicitly but are computed from equilibrium considerations (e.g. Saha of collisional ionization equilibrium) so that X=X(T,rho) and the mean molecular weight now depends on temperature. Then one has $T = (p/\rho) K \mu(T,\rho)$.

Nonlinear equation can be solved with a root finder (typically Brent's method) or a tabulated approach (default). In both case, a 2D table (Ttab) of the temperature T(i,j) is pre-calculated for fixed values of p/rho and rho by MakePV_TemperatureTable() and stored into memory. The table is then used to bracket the root (in the case of a root finder) or to replace the runtime computation with a simpler array indexing operation followed by a combination of lookup table and bilinear (direct or inverse) interpolation.

Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
B. Vaidya
Date
31 Aug, 2014

Definition in file thermal_eos.c.

Macro Definition Documentation

#define PV_TEMPERATURE_TABLE_NX   768

Definition at line 50 of file thermal_eos.c.

#define PV_TEMPERATURE_TABLE_NY   768

Definition at line 53 of file thermal_eos.c.

Function Documentation

int GetPV_Temperature ( double *  v,
double *  T 
)

Compute gas temperature (in Kelvin) from density, pressure and fractions:

\[ T = \left\{\begin{array}{ll} \DS \frac{p}{\rho} K \mu(\vec{X}) &\quad{(\rm with\; chemistry)} \\ \noalign{\medskip} \DS \frac{p}{\rho}K\mu(\vec{X}(T,\rho)) & \quad{(\rm without\; chemistry)} \end{array}\right. \]

where K is the KELVIN macro.

The first relation is used when ion fractions are explicitly evolved by the code (NIONS>0) while the second one is inverted numerically using root finder or lookup table when either LTE or CIE is assumed.

Parameters
[in]v1D array of primitive quantities.
[out]TTemperature in K
Returns
0 on success, 1 on failure.

Definition at line 221 of file thermal_eos.c.

243 {
244  int status;
245 
246 #if NIONS == 0 && PV_TEMPERATURE_TABLE == YES
247 
248  double rho, T1;
249 
250  if (v[PRS] < 0.0) return 1;
251 
252  rho = v[RHO];
253  T1 = v[PRS]/rho;
254  status = Table2DInterpolate(&Ttab, T1, rho, T);
255 
256  if (status != 0){
257  print ("! GetPV_Temperature: table interpolation failure ");
258  print ("[T1 = %8.3e, rho = %8.3e]\n",T1,rho);
259 
260  if (status == -1) return 1; /* hit lower x range (easy fix) */
261  else QUIT_PLUTO(1); /* not so easy to fix */
262  }
263 
264  if (*T < T_CUT_RHOE) return 1;
265  else return 0;
266 
267 #elif NIONS == 0 && PV_TEMPERATURE_TABLE == NO
268 
269  int i,j;
270  double Tmin, Tmax, lnx, lny, rho, prs;
271  struct func_param par;
272 
273  if (v[PRS] < 0.0) return 1;
274 
275  par.v[RHO] = v[RHO];
276  par.T1 = v[PRS]/v[RHO]*KELVIN;
277 
278 /* -- Locate indices (i,j) in the table -- */
279 
280  rho = v[RHO];
281  prs = v[PRS];
282  lnx = log10(prs/rho);
283  lny = log10(rho);
284  i = INT_FLOOR((lnx - Ttab.lnxmin)*Ttab.dlnx_1);
285  j = INT_FLOOR((lny - Ttab.lnymin)*Ttab.dlny_1);
286 
287 /* -------------------------------------------------------------------
288  Since Ttab[j][i] increases with i and j, the min and max between
289  four adjacent node values are the two opposite points on the main
290  diagonal given, respectively, by (i,j) and (i+1,j+1).
291  ------------------------------------------------------------------- */
292 
293  Tmin = MIN(Ttab.f[j][i], Ttab.f[j+1][i+1]);
294  Tmax = MAX(Ttab.f[j][i], Ttab.f[j+1][i+1]);
295 
296 /* status = Ridder(TFunc, &p, Tmin, Tmax, -1, 1.0e-7, &T); */
297  status = Brent(TFunc, &par, Tmin, Tmax, -1, 1.0e-7, T);
298  if (status != 0){
299  print ("! PrimitiveTemperature: could not find root, ");
300  print ("! [Tmin, Tmax] = [%12.6e, %12.6e]\n",Tmin, Tmax);
301  QUIT_PLUTO(1);
302  }
303 
304  if (*T < T_CUT_RHOE || status != 0) return 1;
305  return 0;
306 
307 #else
308 
309  *T = (v[PRS]/v[RHO])*KELVIN*MeanMolecularWeight(v);
310  if (*T < T_CUT_RHOE) return 1;
311  return 0;
312 
313 #endif
314 
315 }
MAX
#define MAX(a, b)
Definition: macros.h:101
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
T_CUT_RHOE
#define T_CUT_RHOE
Sets the lowest cut-off temperature used in the PVTE_LAW equation of state.
Definition: eos.h:69
RHO
#define RHO
Definition: mod_defs.h:19
func_param::T1
double T1
Definition: eos.h:108
TABLE2D::dlny_1
double dlny_1
Definition: structs.h:205
func_param::v
double v[NVAR]
Definition: eos.h:106
TABLE2D::f
double ** f
Definition: structs.h:186
KELVIN
#define KELVIN
Definition: pluto.h:401
MIN
#define MIN(a, b)
Definition: macros.h:104
Brent
int Brent(double(*Func)(double, void *), void *param, double x1, double x2, double abs_acc, double rel_acc, double *xroot)
Definition: math_root_finders.c:14
TABLE2D::dlnx_1
double dlnx_1
Definition: structs.h:204
INT_FLOOR
#define INT_FLOOR(z)
Definition: macros.h:98
j
int j
Definition: analysis.c:2
Table2DInterpolate
int Table2DInterpolate(Table2D *tab, double x, double y, double *f)
Definition: math_table2D.c:430
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
TFunc
static double TFunc(double T, void *par)
Definition: thermal_eos.c:319
func_param
Definition: eos.h:105
i
int i
Definition: analysis.c:2
MeanMolecularWeight
double MeanMolecularWeight(double *V)
Definition: radiat.modified_by_Bhargav.c:178
Ttab
static Table2D Ttab
A 2D table containing pre-computed values of temperature stored at equally spaced node values of Log...
Definition: thermal_eos.c:56
TABLE2D::lnxmin
double lnxmin
lower limit (in log10) in the x-direction
Definition: structs.h:198
TABLE2D::lnymin
double lnymin
lower limit (in log10) in the y-direction
Definition: structs.h:200
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the call graph for this function:

Here is the caller graph for this function:

void MakePV_TemperatureTable ( )

Pre-calculate a 2D table of temperature T(p,rho) by inverting, at specified values of p and rho, the nonlinear equation (only in LTE or CIE)

\[ f(T) = T - \frac{p}{\rho} K\mu(T,\rho) = 0 \qquad\Longrightarrow\qquad T_{ij} = T(x_i, y_j) \qquad \left(x=\log_{10}\frac{p}{\rho}\,,\quad y=\log_{10}\rho\right) \]

For convenience, the table is constructed using equally spaced bins of x=log10(p/rho) and y=log10(rho) where p and rho are in code units. This function must be called only once to initialize the 2D table Ttab.

Definition at line 63 of file thermal_eos.c.

81 {
82  int i,j, status;
83  double prs, rho, T1, Tlo, Thi, T, logK = log10(KELVIN);
84  double mu_lo, mu_hi;
85  struct func_param par;
86 
87  print1 ("> MakePV_TemperatureTable: Generating table...\n");
88 
89 /* --------------------------------------------------------------
90  Initialize table. The two table axis are given by
91  ln(x) = log(p/rho) and ln(y) = log(rho).
92  The lower and upper x-bounds correspond to T/mu = 1 K and
93  T/mu = 10^7 K, respectively.
94  -------------------------------------------------------------- */
95 
96  InitializeTable2D(&Ttab,1.0/KELVIN, 1.e7/KELVIN, PV_TEMPERATURE_TABLE_NX,
97  1.e-7 , 1.e7, PV_TEMPERATURE_TABLE_NY);
98 
99 /* -----------------------------------------------------------------
100  Guess the smallest and largest value of \mu.
101  This is likely to happen at very low and very high temperatures
102  ----------------------------------------------------------------- */
103 
104  GetMu(1.0, 1.0, &mu_lo);
105  GetMu(1.e12, 1.0, &mu_hi);
106 
107  for (j = 0; j < Ttab.ny; j++){
108  for (i = 0; i < Ttab.nx; i++){
109 
110  T1 = Ttab.x[i]*KELVIN; /* = p/rho*KELVIN */
111  rho = Ttab.y[j];
112 
113  par.T1 = T1;
114  par.v[RHO] = rho;
115 
116  /* -- set interval endpoints for root finder -- */
117 
118  Tlo = par.T1*mu_lo; /* T > T1*mu_\min */
119  Thi = par.T1*mu_hi; /* T < T1*mu_\max */
120 
121 /* status = Ridder(TFunc, &par, Tlo, Thi, -1, 1.0e-12, &T); */
122  status = Brent(TFunc, &par, Tlo, Thi, -1, 1.0e-12, &T);
123  if (status != 0) {
124  print1 ("! MakePV_TemperatureTable: ");
125  print1 ("root could not be found [%d]\n",status);
126  QUIT_PLUTO(1);
127  }
128  Ttab.f[j][i] = T;
129  }}
130 
131  FinalizeTable2D(&Ttab);
132  if (prank == 0) WriteBinaryTable2D("T_tab.bin",&Ttab);
133 
134 #if 0 /* Table-to-Table conversion: attempt to interpolate from another table */
135 {
136  int status;
137  double T, p, mu = 1.2;
138  Table2D rhotab; /* test table rho = rho(T,p); */
139 
140  InitializeTable2D(&rhotab, 1.e2, 1.e7, 64, 1.e-2, 1.e4, 80);
141 
142  for (j = 0; j < rhotab.ny; j++){
143  for (i = 0; i < rhotab.nx; i++){
144  T = rhotab.x[i];
145  p = rhotab.y[j];
146 mu = (T/5.e3);
147 mu = 1.2 + 0.0/cosh(mu*mu);
148  rhotab.f[j][i] = p*KELVIN*mu/T; /* rho */
149  }}
150 
151  FinalizeTable2D(&rhotab);
152  if (prank == 0) WriteBinaryTable2D("rhotab.bin",&rhotab);
153 
154 /* -- Build T(p/rho, rho) table -- */
155 
156  status = 0;
157  for (j = 0; j < Ttab.ny; j++){
158  for (i = 0; i < Ttab.nx; i++){
159  rho = Ttab.y[j];
160  p = Ttab.x[i]*rho;
161  status = InverseLookupTable2D(&rhotab, p, rho, &T);
162  if (status != 0){
163  printf ("! Error in inverting table\n");
164  T = -1.0;
165  }
166  Ttab.f[j][i] = T;
167  }}
168 
169  if (prank == 0) WriteBinaryTable2D("T_tab2.bin",&Ttab);
170 
171  printf ("KELVIN = %12.6e\n",KELVIN);
172  exit(1);
173 }
174 #endif /* 0 */
175 }
GetMu
void GetMu(double T, double rho, double *mu)
Definition: pvte_law.c:115
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
TABLE2D::y
double * y
array of y-values (not uniform)
Definition: structs.h:181
PV_TEMPERATURE_TABLE_NY
#define PV_TEMPERATURE_TABLE_NY
Definition: thermal_eos.c:53
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
TABLE2D::nx
int nx
Number of columns or points in the x direction.
Definition: structs.h:174
RHO
#define RHO
Definition: mod_defs.h:19
func_param::T1
double T1
Definition: eos.h:108
PV_TEMPERATURE_TABLE_NX
#define PV_TEMPERATURE_TABLE_NX
Definition: thermal_eos.c:50
TABLE2D::f
double ** f
Definition: structs.h:186
KELVIN
#define KELVIN
Definition: pluto.h:401
prank
int prank
Processor rank.
Definition: globals.h:33
WriteBinaryTable2D
void WriteBinaryTable2D(char *fname, Table2D *tab)
Definition: math_table2D.c:528
Brent
int Brent(double(*Func)(double, void *), void *param, double x1, double x2, double abs_acc, double rel_acc, double *xroot)
Definition: math_root_finders.c:14
InverseLookupTable2D
int InverseLookupTable2D(Table2D *tab, double y, double f, double *x)
Definition: math_table2D.c:196
j
int j
Definition: analysis.c:2
TABLE2D::ny
int ny
Number of rows or points in the y direction.
Definition: structs.h:175
TFunc
static double TFunc(double T, void *par)
Definition: thermal_eos.c:319
func_param
Definition: eos.h:105
InitializeTable2D
void InitializeTable2D(Table2D *tab, double xmin, double xmax, int nx, double ymin, double ymax, int ny)
Definition: math_table2D.c:15
i
int i
Definition: analysis.c:2
Ttab
static Table2D Ttab
A 2D table containing pre-computed values of temperature stored at equally spaced node values of Log...
Definition: thermal_eos.c:56
FinalizeTable2D
void FinalizeTable2D(Table2D *tab)
Definition: math_table2D.c:121
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
TABLE2D::x
double * x
array of x-values (not uniform)
Definition: structs.h:180

Here is the call graph for this function:

Here is the caller graph for this function:

double Pressure ( double *  v,
double  T 
)

Return pressure as a function of temperature, density and fractions. If T(p/rho,rho) has been tabulated, this requires inverting the table by calling InverseLookupTable2D().

Parameters
[in]va pointer to a vector of primitive quantities (only density and fractions will be actually used).
[in]TThe temperature (in Kelvin);
Returns
Pressure in code units.

Definition at line 179 of file thermal_eos.c.

191 {
192 #if PV_TEMPERATURE_TABLE == YES
193 
194  int status;
195  double rho, T1;
196 
197  rho = v[RHO];
198  status = InverseLookupTable2D(&Ttab, rho, T, &T1);
199  if (status != 0){
200  print ("! Pressure: table interpolation failure [rho = %12.6e]\n", rho);
201  QUIT_PLUTO(1);
202  }
203 
204  return T1*rho;
205 
206 #else
207 
208  double mu;
209 
210  #if NIONS == 0
211  GetMu(T, v[RHO], &mu);
212  #else
213  mu = MeanMolecularWeight(v);
214  #endif
215  return v[RHO]*T/(KELVIN*mu);
216 
217 #endif
218 }
GetMu
void GetMu(double T, double rho, double *mu)
Definition: pvte_law.c:115
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
RHO
#define RHO
Definition: mod_defs.h:19
func_param::T1
double T1
Definition: eos.h:108
func_param::v
double v[NVAR]
Definition: eos.h:106
KELVIN
#define KELVIN
Definition: pluto.h:401
InverseLookupTable2D
int InverseLookupTable2D(Table2D *tab, double y, double f, double *x)
Definition: math_table2D.c:196
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
MeanMolecularWeight
double MeanMolecularWeight(double *V)
Definition: radiat.modified_by_Bhargav.c:178
Ttab
static Table2D Ttab
A 2D table containing pre-computed values of temperature stored at equally spaced node values of Log...
Definition: thermal_eos.c:56
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the call graph for this function:

Here is the caller graph for this function:

double TFunc ( double  T,
void *  par 
)
static

Return the nonlinear function f(T) = T - T1 mu(T,rho) used by the root finder to find temperature T.

Parameters
[in]TThe temperature in Kelvin
[in]parA pointer to a func_param structure containing density required for computing internal energy.
Returns
f(T)

Definition at line 319 of file thermal_eos.c.

329 {
330  double T1, rho, mu;
331  struct func_param *p = (struct func_param *) par;
332 
333  T1 = p->T1;
334  rho = p->v[RHO];
335  GetMu(T, rho, &mu);
336  return T - T1*mu;
337 }
GetMu
void GetMu(double T, double rho, double *mu)
Definition: pvte_law.c:115
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
RHO
#define RHO
Definition: mod_defs.h:19
func_param::T1
double T1
Definition: eos.h:108
func_param::v
double v[NVAR]
Definition: eos.h:106
func_param
Definition: eos.h:105

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

Table2D Ttab
static

A 2D table containing pre-computed values of temperature stored at equally spaced node values of Log(p/rho) and Log(rho)

Definition at line 56 of file thermal_eos.c.