PLUTO
Macros | Functions
cooling.h File Reference

Go to the source code of this file.

Macros

#define NIONS   1
 
#define X_HI   NFLX
 
#define frac_Z
 
#define frac_He
 

Functions

double GetMaxRate (double *, double *, double)
 
double H_MassFrac (void)
 
double CompEquil (double, double, double *)
 
void Radiat (double *, double *)
 

Macro Definition Documentation

#define frac_He
Value:
0.082 /* = N(He) / N(H), fractional number density of helium (He)
with respect to hydrogen (H) */

Definition at line 17 of file cooling.h.

#define frac_Z
Value:
1.e-3 /* = N(Z) / N(H), fractional number density of metals (Z)
with respect to hydrogen (H) */

Definition at line 15 of file cooling.h.

#define NIONS   1

Definition at line 10 of file cooling.h.

#define X_HI   NFLX

Definition at line 11 of file cooling.h.

Function Documentation

double CompEquil ( double  N,
double  T,
double *  v 
)
Parameters
[in]nthe particle number density (not needed, but kept for compatibility)
[in]Tthe temperature (in K) for which equilibrium must be found.
[in,out]van array of primitive variables. On input, only density needs to be defined. On output, fractions will be updated to the equilibrium values.

Compute the equilibrium ionization balance for (rho,T)

Parameters
[in]nthe particle number density (not needed, but kept for compatibility)
[in]Tthe temperature (in K) for which equilibrium must be found.
[in,out]van array of primitive variables. On input, only density needs to be defined. On output, fractions will be updated to the equilibrium values.

Compute the equilibrium ionization balance for (rho,T)

Definition at line 44 of file comp_equil.c.

56 {
57  double cr, ci, kr1, kr2, kr3, kr4;
58  double frac_Z, st, tev;
59  long double a, b, c, qc, scrh, az;
60  long double f, h, g;
61  double krv[9];
62 
63  H2RateTables(T, krv);
64  kr1 = krv[0];
65  kr2 = krv[1];
66  kr3 = krv[2];
67  kr4 = krv[3];
68  cr = krv[4];
69  ci = krv[5];
70 
71  kr2 /= kr1;
72  kr3 /= kr1;
73  kr4 /= kr1;
74  kr1 = 1.0;
75  frac_Z = ((1.0 - H_MASS_FRAC - He_MASS_FRAC)/CONST_AZ)*(CONST_AH/H_MASS_FRAC);
76 
77  az = 0.5*CONST_AZ*frac_Z;
78 
79 /* ---------------------------------------------------
80  Solve quadratic equation for fn or hn.
81  The solution depends only on temperature and the
82  physical solution is the one with - sign.
83  --------------------------------------------------- */
84 
85  if (T < 1.e16) { /* Solve for fn */
86  qc = ci/cr;
87  scrh = -0.5*(1.0 + qc);
88  a = (long double)(scrh*(kr2 + kr3*scrh + kr4*qc) - kr1);
89  b = (long double)( 0.5*(kr2 + kr3*scrh + kr4*qc + scrh*(kr3 + 2.0*kr4*az)));
90  c = (long double)(0.5*(0.5*kr3 + kr4*az));
91 
92  /* ------------------------------------------
93  Set fn = 0.0 below T = 300 K to prevent
94  machine accuracy problems
95  ------------------------------------------ */
96 
97  if (T > 300.0){
98  if (b >= 0.0) f = -(b + sqrtl(b*b - 4.0*a*c))/(2.0*a);
99  else f = 2.0*c/(sqrtl(b*b - 4.0*a*c) - b);
100  }else{
101  f = 0.0;
102  }
103  h = qc*f;
104  g = 0.5*(1.0 - f*(1.0 + qc));
105 
106  /* -- Set g = 0 if we are close to machine accuracy -- */
107 
108 /* if ( fabs(1.0 - h - f) < 1.e-15) g = 0.0; */
109 
110  }else{ /* Solve for hn */
111 
112  qc = cr/ci;
113  scrh = -0.5*(1.0 + qc);
114  a = scrh*(kr2*qc + kr3*scrh + kr4) - kr1*qc*qc;
115  b = 0.5*(kr2*qc + kr3*scrh + kr4 + scrh*(kr3 + 2.0*kr4*az));
116  c = 0.5*(0.5*kr3 + kr4*az);
117  h = -(b + sqrtl(b*b - 4.0*a*c))/(2.0*a);
118  f = qc*h;
119  g = 0.5*(1.0 - h*(1.0 + qc));
120  }
121 
122  v[X_HI] = f;
123  v[X_HII] = h;
124  v[X_H2] = g;
125 }
a
static double a
Definition: init.c:135
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
CONST_AH
#define CONST_AH
Definition: pluto.h:250
configure.scrh
tuple scrh
Definition: configure.py:200
frac_Z
#define frac_Z
Definition: cooling.h:15
He_MASS_FRAC
#define He_MASS_FRAC
Definition: pluto.h:551
X_HII
#define X_HII
Definition: cooling.h:31
CONST_AZ
#define CONST_AZ
Mean atomic weight of heavy elements.
Definition: pluto.h:252
X_HI
Definition: cooling.h:110
menu.c
tuple c
Definition: menu.py:375
X_H2
#define X_H2
Definition: cooling.h:30
H_MASS_FRAC
#define H_MASS_FRAC
Definition: pluto.h:542
H2RateTables
void H2RateTables(double, double *)
Definition: radiat.c:7

Here is the call graph for this function:

Here is the caller graph for this function:

double GetMaxRate ( real *  v0,
real *  k1,
real  T0 
)

Return an estimate of the maximum rate (dimension 1/time) in the chemical network. This will serve as a "stiffness" detector in the main ode integrator.

For integration to be carried explicitly all the time, return a small value (1.e-12).

Definition at line 4 of file maxrate.c.

17 {
18  return (1.e6);
19 }
double H_MassFrac ( void  )

Compute the mass fraction X of Hydrogen as function of the composition of the gas. 

        f_H A_H

X = -------------— f_k A_k

where

f_k : is the fractional abundance (by number), f_k = N_k/N_tot of atomic species (no matter ionization degree).

A_K : is the atomic weight

where N_{tot} is the total number density of particles

ARGUMENTS

none

Compute the mass fraction X of Hydrogen as function of the composition of the gas. 

      f_H A_H

X = -------------— f_k A_k

where

f_k : is the fractional abundance (by number), f_k = N_k/N_tot of atomic species (no matter ionization degree).

A_K : is the atomic weight

Note: In this module, f_H = 1.0

where N_{tot} is the total number density of particles

ARGUMENTS

none

Definition at line 721 of file radiat.c.

746 {
747  double mmw2;
748  int i, j;
749 
750  mmw2 = 0.0;
751  for (i = 0; i < (Fe_IONS==0?7:8); i++) {
752  mmw2 += elem_mass[i]*elem_ab[i]; /* Denominator part of X */
753  }
754  return ( (elem_mass[0]*elem_ab[0]) / mmw2);
755 }
elem_mass
const double elem_mass[]
Definition: radiat.c:14
elem_ab
const double elem_ab[]
Definition: radiat.c:9
j
int j
Definition: analysis.c:2
i
int i
Definition: analysis.c:2
Fe_IONS
#define Fe_IONS
Definition: cooling.h:18
void Radiat ( double *  v,
double *  rhs 
)

Cooling for optically thin plasma up to about 200,000 K Plasma composition: H, HeI-II, CI-V, NI-V, OI-V, NeI-V, SI-V Assumed abundances in elem_ab Uses S : Array = Variables vector x line points rhs : output for the system of ODE ibeg, iend : begin and end points of the current line

Cooling for neutral or singly ionized gas: good up to about 35,000 K in equilibrium or shocks in neutral gas up to about 80 km/s. Assumed abundances in ab Uses t : Kelvin dene : electron density cm*-3 fneut : hydrogen neutral fraction (adimensionale) ci,cr : H ionization and recombination rate coefficients

em(1) = TOTAL EMISSIVITY : (ergs cm**3 s**-1) em(2) = Ly alpha + two photon continuum: Aggarwal MNRAS 202, 10**4.3 K em(3) = H alpha: Aggarwal, Case B em(4) = He I 584 + two photon + 623 (all n=2 excitations): Berrington &Kingston,JPB 20 em(5) = C I 9850 + 9823: Mendoza, IAU 103, 5000 K em(6) = C II, 156 micron: Mendoza, 10,000 K em(7) = C II] 2325 A: Mendoza, 15,000 K em(8) = N I 5200 A: Mendoza, 7500 K em(9) = N II 6584 + 6548 A: Mendoza em(10) = O I 63 micron: Mendoza,2500 K em(11) = O I 6300 A + 6363 A: Mendoza, 7500 K em(12) = O II 3727: Mendoza em(13) = Mg II 2800: Mendoza em(14) = Si II 35 micron: Dufton&Kingston, MNRAS 248 em(15) = S II 6717+6727: Mendoza em(16) = Fe II 25 micron: Nussbaumer&Storey em(17) = Fe II 1.6 micron em(18) = thermal energy lost by ionization em(19) = thermal energy lost by recombination (2/3 kT per recombination. The ionization energy lost is not included here.

Provide r.h.s. for tabulated cooling.

Definition at line 94 of file radiat.c.

106 {
107  int ii, k, nv, status;
108  double T, mu, t3, rho, prs, fn, gn, hn, x;
109  double N_H, n_el, cr, ci, rlosst, src_pr, crold,crnew;
110  double kr1, kr2, kr3, kr4, em[20], krvalues[9];
111 
112  static int first_call = 1;
113  static real E_cost, Unit_Time, N_H_rho, frac_He, frac_Z, N_tot;
114 
115  if (first_call) {
116  E_cost = UNIT_LENGTH/UNIT_DENSITY/pow(UNIT_VELOCITY, 3.0);
117  Unit_Time = UNIT_LENGTH/UNIT_VELOCITY;
118 
119  N_H_rho = (UNIT_DENSITY/CONST_amu)*(H_MASS_FRAC/CONST_AH);
120  frac_He = (He_MASS_FRAC/CONST_AHe)*(CONST_AH/H_MASS_FRAC);
121  frac_Z = ((1 - H_MASS_FRAC - He_MASS_FRAC)/CONST_AZ)*(CONST_AH/H_MASS_FRAC);
122  N_tot = N_H_rho*(1.0 + frac_He + frac_Z);
123  H2RateTables(100.0, krvalues);
124  first_call = 0;
125  }
126 
127 /* ---------------------------------------------
128  Force fneut and fmol to stay between [0,1]
129  --------------------------------------------- */
130 
131  NIONS_LOOP(nv){
132  v[nv] = MAX(v[nv], 0.0);
133  v[nv] = MIN(v[nv], 1.0);
134  }
135  v[X_H2] = MIN(v[X_H2], 0.5);
136 
137 /* ---------------------------------------------------
138  Compute temperature from rho, rhoe and fractions
139  --------------------------------------------------- */
140 
141  rho = v[RHO];
142  mu = MeanMolecularWeight(v);
143  if (mu < 0.0){
144  print1 ("! Radiat: mu = %f < 0 \n",mu);
145  QUIT_PLUTO(1);
146  }
147  #if EOS == IDEAL
148  if (v[RHOE] < 0.0) v[RHOE] = g_smallPressure/(g_gamma-1.0);
149  prs = v[RHOE]*(g_gamma-1.0);
150  T = prs/rho*KELVIN*mu;
151  #else
152  status = GetEV_Temperature(v[RHOE], v, &T);
153  if (status != 0) {
154  T = T_CUT_RHOE;
155  v[RHOE] = InternalEnergy(v, T);
156  }
157  #endif
158 
159  N_H = N_H_rho*rho; /* Hydrogen number density
160  N_H = N(X_HI) + 2N(X_H2)+ N(X_HII) */
161  fn = v[X_HI];
162  gn = v[X_H2];
163  hn = v[X_HII];
164 
165 /* Recombination and ionization coefficients */
166 
167  n_el = N_H*(hn + 0.5*CONST_AZ*frac_Z); /* -- electron number density, in cm^{-3} -- */
168 
169  if (n_el < 0){
170  print1 ("! Radiat: negative electron density\n");
171  QUIT_PLUTO(1);
172  }
173 
174  t3 = T/1.e3;
175 
176 /* The Units of kr1, kr2, kr3, kr4, cr, ci = cm^{3}s^{-1}*/
177 
178  H2RateTables(T, krvalues);
179  kr1 = krvalues[0];
180  kr2 = krvalues[1];
181  kr3 = krvalues[2];
182  kr4 = krvalues[3];
183  cr = krvalues[4];
184  ci = krvalues[5];
185 
186  /*
187  kr1 = 3.e-18*st/(1.0 + 0.04*st + 2.e-3*T + 8.e-6*T*T);
188  kr2 = 1.067e-10*pow(tev,2.012)*exp(-(4.463/tev)*pow((1.0 + 0.2472),3.512));
189  kr3 = 1.e-8*exp(-84100.0/T);
190  kr4 = 4.4e-10*pow(T,0.35)*exp(-102000.0/T);
191  */
192 /* ************ ENSURING SUM IS UNITY **********************
193  * This is done to ensure that the sum of speicies
194  * equal to 1. Total Hydrogen number N_H composed of
195  * N_H = nHI + 2.*nH2 + nHII.
196  * while fraction gn = nH2/N_H, thus the rhs[X_H2] has
197  * a pre-factor of 0.5.
198  * And the fact the sum = 1 is ensured if the
199  * corresponding sum of RHS = 0.
200  * i.e., rhs[X_H1] + 2.0*rhs[X_H2] + rhs[X_HII] = 0.0.
201  ********************************************************** */
202 
203  x = n_el/N_H;
204  rhs[X_HI] = Unit_Time*N_H*( gn*(kr2*fn + kr3*gn + kr4*x) + cr*hn*x
205  -fn*(ci*x + kr1*fn));
206  rhs[X_H2] = 0.5*Unit_Time*N_H*(kr1*fn*fn - gn*(kr2*fn + kr3*gn + kr4*x));
207  rhs[X_HII] = Unit_Time*N_H*(ci*fn - cr*hn)*x;
208 /*
209  double sum;
210  sum = 0.0;
211  sum=fabs(rhs[X_HI] + 2.0*rhs[X_H2] + rhs[X_HII]);
212  if (sum > 1.e-9){
213  print("%le > Sum(RHS) != 0 for H!!\n",sum);
214  QUIT_PLUTO(1);
215  }
216 */
217 
218  double logt3 = log(t3);
219 
220 /* The unit of cooling function Lambda in erg cm^{-3} s^{-1}*/
221 
222  /* em[1] = 6.7e-19*exp(-5.86/t3) + 1.6e-18*exp(-11.7/t3);
223  em[2] = (9.5e-22*pow(t3,3.76)/(1. + 0.12*pow(t3,2.1)))*exp(pow((-0.13/t3),3.0))
224  + 3.0e-24*exp(-0.51/t3);*/
225 
226 
227  em[1] = krvalues[6];
228  em[2] = krvalues[7];
229  em[3] = em[1] + em[2];
230 
231  if (t3 < 0.1){ /* Following the Table 8 from Glover & Abel, MNRAS (2008) */
232  em[4] = -16.818342 + logt3*(37.383713 + logt3*(58.145166 + logt3*(48.656103 + logt3*(20.159831 + logt3*(3.8479610)))));
233  }else if(t3 < 1.0){
234  em[4] = -24.311209 + logt3*(3.5692468 + logt3*(-11.332860 + logt3*(-27.850082 + logt3*(-21.328264 + logt3*(-4.2519023)))));
235 
236  }else{
237  em[4] = -24.311209 + logt3*(4.6450521 + logt3*(-3.7209846 + logt3*(5.9369081 + logt3*(-5.5108047 + logt3*(1.5538288)))));
238 
239  }
240  em[4] = POW10(em[4]);
241 
242  em[5] = -23.962112 + logt3*(2.09433740 + logt3*(-0.77151436 + logt3*(0.43693353 + logt3*(-0.14913216 + logt3*(-0.033638326)))));
243 
244  em[5] = POW10(em[5]);
245  em[6] = (fn*em[4] + gn*em[5])*N_H;
246 
247 /* Cooling due to X_H2 rotational vibration. */
248 
249  em[7] = gn*N_H*em[3]/(1.0 + (em[3]/em[6]));
250 
251 /* Cooling due to ionization */
252 
253  em[8] = ci*13.6*1.6e-12*fn;
254 
255 /* Cooling due to radiative recombination */
256 
257  em[9] = cr*0.67*1.6e-12*hn*T/11590.0;
258 
259 /* Cooling due to H2 disassociation */
260 
261  em[10] = 7.18e-12*(kr2*fn*gn + kr3*gn*gn + kr4*gn*(n_el/N_H))*N_H*N_H;
262 
263 /* Cooling due to gas grain cooling */
264 
265  em[11] = krvalues[8]*N_H*N_H*fn*fn; //3.8e-33*st*(T-Tdust)*(1.0 - 0.8*exp(-75.0/T))*N_H*N_H;
266 
267  em[13] = em[11] + em[10] + em[7] + (em[8] + em[9])*n_el*N_H;
268 
269 /* ---------------------------------------------------
270  rlosst is the energy loss in units of erg/cm^3/s;
271  it must be multiplied by cost_E in order to match
272  non-dimensional units.
273  Source term for the neutral fraction scales with
274  UNIT_TIME
275  --------------------------------------------------- */
276 
277  rlosst = em[13];
278  rhs[RHOE] = -E_cost*rlosst;
279  rhs[RHOE] *= 1.0/(1.0 + exp(-(T - g_minCoolingTemp)/100.0)); /* -- lower cutoff -- */
280 }
MAX
#define MAX(a, b)
Definition: macros.h:101
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
NIONS_LOOP
#define NIONS_LOOP(n)
Definition: pluto.h:614
g_gamma
double g_gamma
Definition: globals.h:112
UNIT_DENSITY
#define UNIT_DENSITY
Unit density in gr/cm^3.
Definition: pluto.h:369
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
InternalEnergy
double InternalEnergy(double *v, double T)
Definition: internal_energy.c:255
real
double real
Definition: pluto.h:488
CONST_AH
#define CONST_AH
Definition: pluto.h:250
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
RHO
#define RHO
Definition: mod_defs.h:19
CONST_amu
#define CONST_amu
Atomic mass unit.
Definition: pluto.h:253
frac_He
#define frac_He
Definition: cooling.h:17
frac_Z
#define frac_Z
Definition: cooling.h:15
g_smallPressure
double g_smallPressure
Small value for pressure fix.
Definition: globals.h:110
KELVIN
#define KELVIN
Definition: pluto.h:401
He_MASS_FRAC
#define He_MASS_FRAC
Definition: pluto.h:551
MIN
#define MIN(a, b)
Definition: macros.h:104
UNIT_VELOCITY
#define UNIT_VELOCITY
Unit velocity in cm/sec.
Definition: pluto.h:377
X_HII
#define X_HII
Definition: cooling.h:31
POW10
#define POW10(x)
Definition: radiat.c:3
GetEV_Temperature
int GetEV_Temperature(double rhoe, double *v, double *T)
Definition: internal_energy.c:290
UNIT_LENGTH
#define UNIT_LENGTH
Unit Length in cm.
Definition: pluto.h:373
CONST_AZ
#define CONST_AZ
Mean atomic weight of heavy elements.
Definition: pluto.h:252
X_HI
Definition: cooling.h:110
k
int k
Definition: analysis.c:2
g_minCoolingTemp
double g_minCoolingTemp
The minimum temperature (in K) below which cooling is suppressed.
Definition: globals.h:106
MeanMolecularWeight
double MeanMolecularWeight(double *V)
Definition: radiat.modified_by_Bhargav.c:178
X_H2
#define X_H2
Definition: cooling.h:30
CONST_AHe
#define CONST_AHe
Atomic weight of Helium.
Definition: pluto.h:251
H2RateTables
void H2RateTables(double T, double *krvals)
Definition: radiat.c:7
H_MASS_FRAC
#define H_MASS_FRAC
Definition: pluto.h:542
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125