PLUTO
Macros | Functions
pvte_law_H+.c File Reference

PVTE_LAW for a partially ionized gas. More...

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

Go to the source code of this file.

Macros

#define DIFF_ORDER
 
#define INTE_EXACT   1 /* Compute internal energy exactly in Gamma1() */
 

Functions

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

Detailed Description

PVTE_LAW for a partially ionized gas.

Compute the internal energy for a partially ionized hydrogen gas.

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

Reference

  • PLUTO Users' guide.

Definition in file pvte_law_H+.c.

Macro Definition Documentation

#define DIFF_ORDER
Value:
4 /* = 2 or 4 for 2nd or 4th accurate approximations
to derivatives in Gamma1() */

Definition at line 19 of file pvte_law_H+.c.

#define INTE_EXACT   1 /* Compute internal energy exactly in Gamma1() */

Definition at line 21 of file pvte_law_H+.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.

Parameters
[in]v1D array of primitive quantities
Returns
Value of first adiabatic index Gamma1.

Definition at line 113 of file pvte_law_H+.c.

144 {
145  double gmm1, T;
146  double cv, mu, chirho = 1.0, chiT = 1.0, rho, rhoe;
147  double epp, emm, ep, em, delta = 1.e-3;
148  double Tpp, Tmm, Tp, Tm, mupp, mumm, mup, mum, rhopp, rhomm, rhop, rhom;
149  double dmu_dT, dmu_drho, de_dT;
150 
151  return 1.666667;
152 
153 /* ---------------------------------------------
154  Obtain temperature and fractions.
155  --------------------------------------------- */
156 
157  #if NIONS == 0
158  GetPV_Temperature(v, &T);
159  GetMu(T, v[RHO], &mu);
160  #else
161  mu = MeanMolecularWeight(v);
162  T = v[PRS]/v[RHO]*KELVIN*mu;
163  #endif
164 
165 /* ---------------------------------------------------
166  Compute cV (Specific heat capacity for constant
167  volume) using centered derivative.
168  cV will be in code units. The corresponding cgs
169  units will be the same velocity^2/Kelvin.
170  --------------------------------------------------- */
171 
172  Tp = T*(1.0 + delta);
173  Tm = T*(1.0 - delta);
174  Tmm = T*(1.0 - 2.0*delta);
175  Tpp = T*(1.0 + 2.0*delta);
176 
177  #if INTE_EXACT
178  emm = InternalEnergyFunc(v, Tmm)/v[RHO]; /* in code units */
179  em = InternalEnergyFunc(v, Tm)/v[RHO]; /* in code units */
180  ep = InternalEnergyFunc(v, Tp)/v[RHO]; /* in code units */
181  epp = InternalEnergyFunc(v, Tpp)/v[RHO]; /* in code units */
182  #else
183  emm = InternalEnergy(v, Tmm)/v[RHO]; /* in code units */
184  em = InternalEnergy(v, Tm)/v[RHO]; /* in code units */
185  ep = InternalEnergy(v, Tp)/v[RHO]; /* in code units */
186  epp = InternalEnergy(v, Tpp)/v[RHO]; /* in code units */
187  #endif
188 
189  #if DIFF_ORDER == 2
190  de_dT = (ep - em)/(2.0*delta*T);
191  #elif DIFF_ORDER == 4
192  de_dT = (emm - 8.0*em + 8.0*ep - epp)/(12.0*delta*T);
193  #else
194  #error Order not allowed in Gamma1()
195  #endif
196 
197  cv = de_dT; /* this is code units. */
198 
199  #if NIONS == 0
200  rho = v[RHO];
201 
202  GetMu(Tp, rho, &mup);
203  GetMu(Tm, rho, &mum);
204  GetMu(Tpp, rho, &mupp);
205  GetMu(Tmm, rho, &mumm);
206 
207  #if DIFF_ORDER == 2
208  dmu_dT = (mup - mum)/(2.0*delta*T);
209  #elif DIFF_ORDER == 4
210  dmu_dT = (mumm - 8.0*mum + 8.0*mup - mupp)/(12.0*delta*T);
211  #else
212  #error Order not allowed in Gamma1()
213  #endif
214 
215  chiT = 1.0 - T/mu*dmu_dT;
216 
217  rhop = rho*(1.0 + delta);
218  rhom = rho*(1.0 - delta);
219  rhopp = rho*(1.0 + 2.0*delta);
220  rhomm = rho*(1.0 - 2.0*delta);
221 
222  GetMu(T, rhop, &mup);
223  GetMu(T, rhom, &mum);
224  GetMu(T, rhopp, &mupp);
225  GetMu(T, rhomm, &mumm);
226 
227  #if DIFF_ORDER == 2
228  dmu_drho = (mup - mum)/(2.0*delta*rho);
229  #elif DIFF_ORDER == 4
230  dmu_drho = (mumm - 8.0*mum + 8.0*mup - mupp)/(12.0*delta*rho);
231  #else
232  #error Order not allowed in Gamma1()
233  #endif
234  chirho = 1.0 - rho/mu*dmu_drho;
235  #endif
236 
237 /* --------------------------------------------
238  Compute first adiabatic index
239  -------------------------------------------- */
240 /*
241 printf ("gmm1: prs = %12.6e, T = %12.6e, chiT = %12.6e, chirho = %12.6e\n",
242  v[PRS], T, chiT, chirho);
243 */
244  gmm1 = v[PRS]/(cv*T*v[RHO])*chiT*chiT + chirho;
245 
246 //printf ("gmm1 = %f\n",gmm1);
247  return gmm1;
248 
249 }
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
GetMu
void GetMu(double T, double rho, double *mu)
Definition: pvte_law_H+.c:96
GetPV_Temperature
int GetPV_Temperature(double *v, double *T)
Definition: thermal_eos.c:221
InternalEnergyFunc
double InternalEnergyFunc(double *v, double T)
Definition: pvte_law_H+.c:24

Here is the call 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

Definition at line 96 of file pvte_law_H+.c.

106 {
107  double x;
108  x = SahaXFrac(T, rho);
109  *mu = 1.0/(1.0 + x);
110 }
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
SahaXFrac
static double SahaXFrac(double T, double rho)
Definition: pvte_law_H+.c:72

Here is the call graph for this function:

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 24 of file pvte_law_H+.c.

38 {
39  double x, rho, e, mu, kT;
40  double rhoe, chi = 13.6*CONST_eV;
41  double p0 = UNIT_DENSITY*UNIT_VELOCITY*UNIT_VELOCITY;
42 
43  kT = CONST_kB*T;
44  x = SahaXFrac(T, v[RHO]);
45  GetMu(T,v[RHO],&mu);
46  rho = v[RHO]*UNIT_DENSITY;
47  e = 1.5*kT/(mu*CONST_amu) + chi*x/CONST_mH;
48 
49 /*
50 FILE *fp;
51 double p, lnT;
52 fp = fopen("gamma.dat","w");
53 v[RHO] = 1.0;
54 for (lnT = 2; lnT < 6; lnT += 0.01){
55  T = pow(10.0, lnT);
56  kT = CONST_kB*T;
57  x = SahaXFrac(T, v[RHO]);
58  GetMu(T,v[RHO],&mu);
59  rho = v[RHO]*UNIT_DENSITY;
60  e = 1.5*kT/(mu*CONST_amu) + chi*x/CONST_mH;
61 p = kT/(mu*CONST_amu)*rho;
62 fprintf (fp,"%12.6e %12.6e %12.6e\n",T,p/(rho*e) + 1.0, x);
63 }
64 fclose(fp);
65 exit(1);
66 */
67 
68  return rho*e/p0; /* Convert to code units */
69 }
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
SahaXFrac
static double SahaXFrac(double T, double rho)
Definition: pvte_law_H+.c:72
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
RHO
#define RHO
Definition: mod_defs.h:19
CONST_amu
#define CONST_amu
Atomic mass unit.
Definition: pluto.h:253
func_param::rhoe
double rhoe
Definition: eos.h:107
func_param::v
double v[NVAR]
Definition: eos.h:106
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
GetMu
void GetMu(double T, double rho, double *mu)
Definition: pvte_law_H+.c:96

Here is the call graph for this function:

Here is the caller graph for this function:

double SahaXFrac ( double  T,
double  rho 
)
static

Use Saha equation to compute the degree of ionization.

Definition at line 72 of file pvte_law_H+.c.

78 {
79  double me, kT, h3, c, x, n;
80  double chi = 13.6*CONST_eV;
81 
82  rho *= UNIT_DENSITY;
83 
84  me = 2.0*CONST_PI*CONST_me;
85  kT = CONST_kB*T;
86  h3 = CONST_h*CONST_h*CONST_h;
87 
88  n = rho/CONST_mp; /* = n(protons) + n(neutrals) not n(total) */
89  c = me*kT*sqrt(me*kT)/(h3*n)*exp(-chi/kT);
90 
91  x = 2.0/(sqrt(1.0 + 4.0/c) + 1.0);
92  return x;
93 }
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
n
static int n
Definition: analysis.c:3
CONST_eV
#define CONST_eV
Electron Volt in erg.
Definition: pluto.h:256
CONST_mp
#define CONST_mp
Proton mass.
Definition: pluto.h:261
CONST_me
#define CONST_me
Electron mass.
Definition: pluto.h:263
CONST_kB
#define CONST_kB
Boltzmann constant.
Definition: pluto.h:259
menu.c
tuple c
Definition: menu.py:375
CONST_PI
#define CONST_PI
.
Definition: pluto.h:269

Here is the caller graph for this function: