PLUTO
Macros | Functions
pvte_law.c File Reference

Basic interface definition for the PVTE_LAW EoS. More...

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

Go to the source code of this file.

Macros

#define HELIUM_IONIZATION   NO
 
#define DEG_x   0 /* hydrogen ionization degree */
 
#define DEG_y   1 /* molecular hydrogen dissociation degree */
 
#define DEG_z1   2 /* helium single ionization degree */
 
#define DEG_z2   3 /* helium double ionization degree */
 

Functions

void GetFuncDum (double, double *)
 
static void GetSahaHFracs (double, double, double *)
 
double InternalEnergyFunc (double *v, double T)
 
void GetMu (double T, double rho, double *mu)
 
double Gamma1 (double *v)
 

Detailed Description

Basic interface definition for the PVTE_LAW EoS.

Collect the basic set of functions required by the PVTE_LAW equation of state:

  • InternalEnergyFunc() defines the gas internal energy as a function temperature and ionization fractions (for non-equilibrium chemistry) or temperature and density (in Local Thermodynamic Equilibrium - LTE - or Collisional Ionization Equilibrium - CIE).
  • GetMu() computes the mean molecular weight (in LTE or CIE).

Auxiliary functions (needed by the the previous ones) are also included here.

The current version implements the EoS given by D'Angelo (2013).

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

Reference

  • D'Angelo et. al ApJ 778, 2013 (Eq. [26-27])

Definition in file pvte_law.c.

Macro Definition Documentation

#define DEG_x   0 /* hydrogen ionization degree */

Definition at line 33 of file pvte_law.c.

#define DEG_y   1 /* molecular hydrogen dissociation degree */

Definition at line 34 of file pvte_law.c.

#define DEG_z1   2 /* helium single ionization degree */

Definition at line 35 of file pvte_law.c.

#define DEG_z2   3 /* helium double ionization degree */

Definition at line 36 of file pvte_law.c.

#define HELIUM_IONIZATION   NO

Definition at line 31 of file pvte_law.c.

Function Documentation

double Gamma1 ( double *  v)

Calculate the value of the first adiabatic index:

\[ \Gamma_1 = \frac{1}{c_V}\frac{p}{\rho T} \chi_T^2 + \chi_\rho^2 \qquad{\rm where}\quad \chi_T = \left(\pd{\log p}{\log T}\right)_{\rho} = 1 - \pd{\log{\mu}}{\log T} \,;\quad \chi_\rho = \left(\pd{\log p}{\log\rho}\right)_{T} = 1 - \pd{\log{\mu}}{\log\rho} \,;\quad \]

where p and rho are in c.g.s units. Note that if species are evolved explicitly (non-equilibrium chemistry), we set chi=1.

The heat capacity at constant volume, cV, is defined as the derivative of specific internal energy with respect to temperature:

\[ c_V = \left.\pd{e}{T}\right|_V \]

and it is computed numerically using a centered derivative.

This function is needed (at present) only when computing the sound speed in the Riemann solver. Since this is only needed for an approximated value, 5/3 (upper bound) should be ok.

Reference

  • D'Angelo et. al ApJ 778, 2013 (Eq. [26-27])
Parameters
[in]v1D array of primitive quantities
Returns
Value of first adiabatic index Gamma1.

Definition at line 231 of file pvte_law.c.

265 {
266  double gmm1, T;
267  double cv, mu, chirho = 1.0, chiT = 1.0, rho, rhoe;
268  double ep, em, delta = 1.e-2;
269  double Tp, Tm, mup, mum, rhop, rhom;
270  double dmu_dT, dmu_drho, de_dT;
271 
272  return 1.6667; /* !! For the current implementation, we always return
273  5/3 rather then computing the actual Gamma1 index
274  (which is quite expensive).
275  We do not think this is a major problem since
276  Gamma1 only appears in the estimate of the sound speed
277  needed by the Riemann solver. An upper bound will do. */
278 /* ---------------------------------------------
279  Obtain temperature and fractions.
280  --------------------------------------------- */
281 
282  #if NIONS == 0
283  GetPV_Temperature(v, &T);
284  GetMu(T, v[RHO], &mu);
285  #else
286  mu = MeanMolecularWeight(v);
287  T = v[PRS]/v[RHO]*KELVIN*mu;
288  #endif
289 
290 /* ---------------------------------------------------
291  Compute cV (Specific heat capacity for constant
292  volume) using centered derivative.
293  cV will be in code units. The corresponding cgs
294  units will be the same velocity^2/Kelvin.
295  --------------------------------------------------- */
296 
297  Tp = T*(1.0 + delta);
298  Tm = T*(1.0 - delta);
299  em = InternalEnergy(v, Tm)/v[RHO]; /* in code units */
300  ep = InternalEnergy(v, Tp)/v[RHO]; /* in code units */
301 
302  de_dT = (ep - em)/(2.0*delta*T);
303  cv = de_dT; /* this is code units. */
304 
305  #if NIONS == 0
306  rho = v[RHO];
307 
308  GetMu(Tp, rho, &mup);
309  GetMu(Tm, rho, &mum);
310  dmu_dT = (mup - mum)/(2.0*delta*T);
311  chiT = 1.0 - T/mu*dmu_dT;
312 
313  rhop = rho*(1.0 + delta);
314  rhom = rho*(1.0 - delta);
315  GetMu(T, rhop, &mup);
316  GetMu(T, rhom, &mum);
317  dmu_drho = (mup - mum)/(2.0*delta*rho);
318  chirho = 1.0 - rho/mu*dmu_drho;
319  #endif
320 
321 /* --------------------------------------------
322  Compute first adiabatic index
323  -------------------------------------------- */
324 
325  gmm1 = v[PRS]/(cv*T*v[RHO])*chiT*chiT + chirho;
326 
327  return gmm1;
328 }
GetMu
void GetMu(double T, double rho, double *mu)
Definition: pvte_law.c:115
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
InternalEnergy
double InternalEnergy(double *v, double T)
Definition: internal_energy.c:255
RHO
#define RHO
Definition: mod_defs.h:19
func_param::rhoe
double rhoe
Definition: eos.h:107
func_param::v
double v[NVAR]
Definition: eos.h:106
KELVIN
#define KELVIN
Definition: pluto.h:401
MeanMolecularWeight
double MeanMolecularWeight(double *V)
Definition: radiat.modified_by_Bhargav.c:178
GetPV_Temperature
int GetPV_Temperature(double *v, double *T)
Definition: thermal_eos.c:221

Here is the call graph for this function:

Here is the caller graph for this function:

void GetFuncDum ( double  T,
double *  funcdum_val 
)

Interpolate the value of a function of zetaR from the table that is used to estimate the value of EH2.

Parameters
[in]TValue of temperature in kelvin.
[out]*funcdum_valPointer to the value of function of zetaR
Returns
This function has no return value.

Reference: D'Angelo, G. et al, ApJ 778 2013.

Definition at line 16 of file zeta_tables.c.

29 {
30  int klo, khi, kmid;
31  int nsteps = 5000;
32  double mu, Tmid, dT;
33  static double *lnT, *funcdum;
34  double y, dy;
35  int indx;
36 
37 /* -------------------------------------------
38  Make table on first call
39  ------------------------------------------- */
40 
41  if (lnT == NULL){
42  lnT = ARRAY_1D(nsteps, double);
43  funcdum = ARRAY_1D(nsteps, double);
44  MakeZetaTables(lnT, funcdum, nsteps);
45  }
46  y = log(T);
47 
48 /* -------------------------------------------------
49  Since the table has regular spacing in log T,
50  we divide by the increment to find the nearest
51  node in the table.
52  ------------------------------------------------- */
53 
54  if (y > lnT[nsteps-2]) {
55  *funcdum_val = funcdum[nsteps-2];
56  } else if (y < lnT[0]) {
57  *funcdum_val = funcdum[0];
58  } else{
59  dy = lnT[1] - lnT[0];
60  indx = floor((y - lnT[0])/dy);
61 
62  if (indx >= nsteps || indx < 0){
63  print1 ("! GetFuncDum: indx out of range, indx = %d\n",indx);
64  print1 ("! T = %12.6e\n",T);
65  QUIT_PLUTO(1);
66  }
67  *funcdum_val = (funcdum[indx]*(lnT[indx+1] - y)
68  + funcdum[indx+1]*(y - lnT[indx]))/dy;
69  }
70 }
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
MakeZetaTables
static void MakeZetaTables(double *, double *, int)
Definition: zeta_tables.c:73
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the caller graph for this function:

void GetMu ( double  T,
double  rho,
double *  mu 
)

Calculate the mean molecular weight for the case in which hydrogen fractions are estimated using Saha Equations.

Parameters
[in]TGas temperature in Kelvin.
[in]rhoGas density (code units)
[out]muMean molecular weight

Reference

  • Black and Bodenheimer, ApJ 199, 1975 (Eq. 8)
  • D'Angelo, G. et al ApJ 778, 2013 (Eqs. 15, 16.)

Definition at line 115 of file pvte_law.c.

128 {
129  double f[4];
130 
131  GetSahaHFracs(T, rho, f);
132  #if HELIUM_IONIZATION == YES
133  *mu = 4.0/( 2.*H_MASS_FRAC*(1.0 + f[DEG_y] + 2.0*f[DEG_x]*f[DEG_y])
134  + He_MASS_FRAC*(1.0 + f[DEG_z1]*(1.0 + f[DEG_z2])));
135  #else
136  *mu = 4.0/(2.*H_MASS_FRAC*(1.0 + f[DEG_y] + 2.0*f[DEG_x]*f[DEG_y])
137  + He_MASS_FRAC);
138  #endif
139 }
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
GetSahaHFracs
static void GetSahaHFracs(double, double, double *)
Definition: pvte_law.c:142
DEG_y
#define DEG_y
Definition: pvte_law.c:34
He_MASS_FRAC
#define He_MASS_FRAC
Definition: pluto.h:551
DEG_z2
#define DEG_z2
Definition: pvte_law.c:36
DEG_x
#define DEG_x
Definition: pvte_law.c:33
DEG_z1
#define DEG_z1
Definition: pvte_law.c:35
H_MASS_FRAC
#define H_MASS_FRAC
Definition: pluto.h:542

Here is the call graph for this function:

Here is the caller graph for this function:

void GetSahaHFracs ( double  T,
double  rho,
double *  fdeg 
)
static

Compute degree of ionization and dissociation using Saha Equations. The quadratic equation ay^2 + by + c = 0 is solved using a = 1, c = -b (for hydrogen). A similar expression is used for Helium.

Parameters
[in]TGas temperature in Kelvin.
[in]rhoGas density (in code units)
[out]fdegarray of ionization/dissociation degrees.

Reference

  • D'Angelo, G. et al ApJ 778, 2013 (Eqs. 15, 16.)

Definition at line 142 of file pvte_law.c.

156 {
157  double rhs1, rhs2, rhs3;
158  double b,c, scrh;
159  double kT = CONST_kB*T;
160  double hbar = CONST_h/(2.0*CONST_PI);
161 
162  rho *= UNIT_DENSITY;
163 
164  rhs1 = CONST_mH/(2.0*H_MASS_FRAC*rho);
165  rhs2 = CONST_mH*kT/(4.0*CONST_PI*hbar*hbar);
166  rhs3 = exp(-4.48*CONST_eV/kT);
167  b = rhs1*rhs2*sqrt(rhs2)*rhs3;
168 
169  fdeg[DEG_y] = 2.0/(1.0 + sqrt(1.0 + 4.0/b)); /* solution of quadratic equation */
170 
171  rhs1 = CONST_mH/(H_MASS_FRAC*rho);
172  rhs2 = CONST_me*kT/(2.0*CONST_PI*hbar*hbar);
173  rhs3 = exp(-13.60*CONST_eV/kT);
174  b = rhs1*rhs2*sqrt(rhs2)*rhs3;
175 
176  fdeg[DEG_x] = 2.0/(1.0 + sqrt(1.0 + 4.0/b)); /* solution of quadratic equation */
177 
178  #if HELIUM_IONIZATION == YES
179  rhs3 = 4.0*CONST_mH/rho*rhs2*sqrt(rhs2)*exp(-24.59*CONST_eV/kT);
180  b = 4.0/He_MASS_FRAC*(H_MASS_FRAC + rhs3);
181  c = -rhs3*4.0/He_MASS_FRAC;
182  fdeg[DEG_z1] = -2.0*c/(b + sqrt(b*b - 4.0*c));
183 
184  rhs3 = CONST_mH/rho*rhs2*sqrt(rhs2)*exp(-54.42*CONST_eV/kT);
185  b = 4.0/He_MASS_FRAC*(H_MASS_FRAC + 0.25*He_MASS_FRAC + rhs3);
186  c = -rhs3*4.0/He_MASS_FRAC;
187  fdeg[DEG_z2] = -2.0*c/(b + sqrt(b*b - 4.0*c));
188  #endif
189 }
CONST_h
#define CONST_h
Planck Constant.
Definition: pluto.h:258
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
CONST_eV
#define CONST_eV
Electron Volt in erg.
Definition: pluto.h:256
configure.scrh
tuple scrh
Definition: configure.py:200
DEG_y
#define DEG_y
Definition: pvte_law.c:34
He_MASS_FRAC
#define He_MASS_FRAC
Definition: pluto.h:551
CONST_me
#define CONST_me
Electron mass.
Definition: pluto.h:263
DEG_z2
#define DEG_z2
Definition: pvte_law.c:36
CONST_kB
#define CONST_kB
Boltzmann constant.
Definition: pluto.h:259
CONST_mH
#define CONST_mH
Hydrogen atom mass.
Definition: pluto.h:264
menu.c
tuple c
Definition: menu.py:375
DEG_x
#define DEG_x
Definition: pvte_law.c:33
CONST_PI
#define CONST_PI
.
Definition: pluto.h:269
DEG_z1
#define DEG_z1
Definition: pvte_law.c:35
H_MASS_FRAC
#define H_MASS_FRAC
Definition: pluto.h:542

Here is the caller graph for this function:

double InternalEnergyFunc ( double *  v,
double  T 
)

Compute the gas internal energy as a function of temperature and fractions (or density):

  • rhoe = rhoe(T,rho) in LTE or CIE;
  • rhoe = rhoe(T,X) in non-equilibrium chemistry.
Parameters
[in]v1D Array of primitive variables containing density and species. Other variables are ignored.
[in]TGas temperature
Returns
The gas internal energy (rhoe) in code units.

Definition at line 47 of file pvte_law.c.

61 {
62  double eH, eHe, eHpH, eHplus, eH2, rhoe, eHeplus, eHeplusplus;
63  double func_zetaR;
64  double f[4];
65 
66  #if NIONS == 0
67  GetSahaHFracs(T, v[RHO], f);
68  #else
69  GetHFracs (v, f);
70  if (f[DEG_x] > 1.001 || f[DEG_y] > 1.001){
71  printf ("! InternalEnergy: x or y are too large %f, %f\n",
72  f[DEG_x],f[DEG_y]);
73  QUIT_PLUTO(1);
74  }
75  #endif
76 
77  GetFuncDum(T, &func_zetaR); /* = (1.5 + e(rot) + e(vib)) */
78 /* func_zetaR = 1.5; to recover ideal EoS */
79 
80 /* -- Estimate contributions to Egas -- */
81 
82  eH = 1.5*H_MASS_FRAC*(1.0 + f[DEG_x])*f[DEG_y];
83  eHe = 3.0*He_MASS_FRAC/8.0;
84  eH2 = 0.5*H_MASS_FRAC*(1.0 - f[DEG_y])*func_zetaR;
85 
86 /* -- constant terms -- */
87 
88  #if COOLING == NO
89  eHpH = 4.48*CONST_eV*H_MASS_FRAC*f[DEG_y]/(2.0*CONST_kB*T);
90  eHplus = 13.60*CONST_eV*H_MASS_FRAC*f[DEG_x]*f[DEG_y]/(CONST_kB*T);
91 
92  #if HELIUM_IONIZATION == YES
93  eHeplus = 24.59*CONST_eV*He_MASS_FRAC*f[DEG_z1]*(1.0 - f[DEG_z2])/(4.0*CONST_kB*T);
94  eHeplusplus = 54.42*CONST_eV*He_MASS_FRAC*f[DEG_z1]*f[DEG_z2]/(4.0*CONST_kB*T);
95  #else
96  eHeplus = eHeplusplus = 0.0;
97  #endif
98 
99  #else
100  eHpH = eHplus = eHeplus = eHeplusplus = 0.0;
101  #endif
102 
103 /* ----------------------------------------------------------------
104  Compute rhoe in cgs except for density which is left in
105  code units (to avoid extra multiplication and division later
106  ---------------------------------------------------------------- */
107 
108  rhoe = (eH2 + eH + eHe + eHpH + eHplus + eHeplus + eHeplusplus)*(CONST_kB*T*v[RHO]/CONST_mH);
109 
110  return rhoe/(UNIT_VELOCITY*UNIT_VELOCITY); /* convert to code units */
111 }
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
GetSahaHFracs
static void GetSahaHFracs(double, double, double *)
Definition: pvte_law.c:142
CONST_eV
#define CONST_eV
Electron Volt in erg.
Definition: pluto.h:256
GetFuncDum
void GetFuncDum(double, double *)
Definition: zeta_tables.c:16
RHO
#define RHO
Definition: mod_defs.h:19
func_param::rhoe
double rhoe
Definition: eos.h:107
func_param::v
double v[NVAR]
Definition: eos.h:106
DEG_y
#define DEG_y
Definition: pvte_law.c:34
He_MASS_FRAC
#define He_MASS_FRAC
Definition: pluto.h:551
DEG_z2
#define DEG_z2
Definition: pvte_law.c:36
UNIT_VELOCITY
#define UNIT_VELOCITY
Unit velocity in cm/sec.
Definition: pluto.h:377
CONST_kB
#define CONST_kB
Boltzmann constant.
Definition: pluto.h:259
CONST_mH
#define CONST_mH
Hydrogen atom mass.
Definition: pluto.h:264
DEG_x
#define DEG_x
Definition: pvte_law.c:33
DEG_z1
#define DEG_z1
Definition: pvte_law.c:35
H_MASS_FRAC
#define H_MASS_FRAC
Definition: pluto.h:542
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: