PLUTO
cooling.h File Reference

Go to the source code of this file.

Macros

#define NIONS   3
 
#define X_HI   (NFLX)
 
#define X_H2   (NFLX + 1)
 
#define X_HII   (NFLX + 2)
 

Functions

void CompEquil (double n, double T, double *v)
 
double GetMaxRate (double *, double *, double)
 
void Radiat (double *, double *)
 
void NormalizeIons (double *)
 
void H2RateTables (double, double *)
 

Macro Definition Documentation

#define NIONS   3

Definition at line 28 of file cooling.h.

#define X_H2   (NFLX + 1)

Definition at line 30 of file cooling.h.

#define X_HI   (NFLX)

Definition at line 29 of file cooling.h.

#define X_HII   (NFLX + 2)

Definition at line 31 of file cooling.h.

Function Documentation

void CompEquil ( double  N,
double  T,
double *  v0 
)
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 }
static double a
Definition: init.c:135
tuple T
Definition: Sph_disk.py:33
#define CONST_AH
Definition: pluto.h:250
tuple scrh
Definition: configure.py:200
#define frac_Z
Definition: cooling.h:15
#define He_MASS_FRAC
Definition: pluto.h:551
#define X_HII
Definition: cooling.h:31
#define CONST_AZ
Mean atomic weight of heavy elements.
Definition: pluto.h:252
Definition: cooling.h:110
tuple c
Definition: menu.py:375
#define X_H2
Definition: cooling.h:30
#define H_MASS_FRAC
Definition: pluto.h:542
void H2RateTables(double, double *)
Definition: radiat.c:7
double GetMaxRate ( double *  v0,
double *  k1,
double  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 }

Here is the caller graph for this function:

void H2RateTables ( double  T,
double *  krvals 
)

On first call, compute and store rate coefficient arrays needed for the H2_COOL module. On subsequent calls, use linear interpolation between adjacent tabulated values to obtain the coefficients at the desired temperature T.

Definition at line 7 of file radiat.c.

15 {
16  int i, indx_lo, indx_hi;
17  double st, t3, tev, lnT, lnTmin, lnTmax, dlnT;
18  double scrh, dum_lo, dum_hi;
19  double Tmin, Tmax;
20  static double *lnTarr, *Rate_arr1, *Rate_arr2, *Rate_arr3, *Rate_arr4;
21  static double *Rate_arr5, *Rate_arr6;
22  static double *Emiss_arr1, *Emiss_arr2, *Emiss_arr3;
23  double Tdust = 15.0;
24 
25  Tmin = 1.0e-2;
26  Tmax = 1.0e8;
27  lnTmin = log10(Tmin);
28  lnTmax = log10(Tmax);
29  dlnT = (lnTmax - lnTmin)/((double)TABLE_NPT - 1.0);
30 
31  T = MAX(T, Tmin); /* We constrain local T to lie between [Tmin, Tmax] so */
32  T = MIN(T, Tmax); /* that coefficients will be constant for values of
33  temperature overflow / underflow. */
34 
35  if (lnTarr == NULL){
36  lnTarr = ARRAY_1D(TABLE_NPT, double);
37  Rate_arr1 = ARRAY_1D(TABLE_NPT, double);
38  Rate_arr2 = ARRAY_1D(TABLE_NPT, double);
39  Rate_arr3 = ARRAY_1D(TABLE_NPT, double);
40  Rate_arr4 = ARRAY_1D(TABLE_NPT, double);
41  Rate_arr5 = ARRAY_1D(TABLE_NPT, double);
42  Rate_arr6 = ARRAY_1D(TABLE_NPT, double);
43  Emiss_arr1 = ARRAY_1D(TABLE_NPT, double);
44  Emiss_arr2 = ARRAY_1D(TABLE_NPT, double);
45  Emiss_arr3 = ARRAY_1D(TABLE_NPT, double);
46 
47  for(i = 0 ; i < TABLE_NPT; i++){
48  lnTarr[i] = lnTmin + i*dlnT;
49  T = POW10(lnTarr[i]);
50  st = sqrt(T);
51  t3 = T/1.e3;
52  tev = T*CONST_kB/CONST_eV;
53 
54  Rate_arr1[i] = 3.e-18*st/(1.0 + 0.04*st + 2.e-3*T + 8.e-6*T*T); /* fa ~ 1 and Tg << T */
55  Rate_arr2[i] = 1.067e-10*pow(tev,2.012)*exp(-(4.463/tev)*pow((1.0 + 0.2472),3.512));
56  Rate_arr3[i] = 1.e-8*exp(-84100.0/T);
57  Rate_arr4[i] = 4.4e-10*pow(T,0.35)*exp(-102000.0/T);
58  Rate_arr5[i] = 2.6e-11/st;
59  Rate_arr6[i] = 1.08e-8*st*exp(-157890.0/T)/(13.6*13.6);
60  Emiss_arr1[i] = 6.7e-19*exp(-5.86/t3) + 1.6e-18*exp(-11.7/t3);
61  Emiss_arr2[i] = (9.5e-22*pow(t3,3.76)/(1. + 0.12*pow(t3,2.1)))*exp(pow((-0.13/t3),3.0))
62  + 3.0e-24*exp(-0.51/t3);
63  Emiss_arr3[i] = 3.8e-33*st*(T-Tdust)*(1.0 - 0.8*exp(-75.0/T));
64  }
65  }else{ /* Perform linear interpolation */
66  lnT = log10(T);
67 
68  scrh = (lnT - lnTmin)/dlnT;
69  indx_lo = INT_FLOOR(scrh);
70  indx_hi = indx_lo + 1;
71 
72  dum_lo = (lnTarr[indx_hi] - lnT)/dlnT;
73  dum_hi = (lnT - lnTarr[indx_lo])/dlnT;
74 
75  if (indx_lo < 0){
76  print ("! Make_RateTable: Cannot Locate Index for T = %12.6e\n",T);
77  QUIT_PLUTO(1);
78  }else{
79  krvals[0] = Rate_arr1[indx_lo]*dum_lo + Rate_arr1[indx_hi]*dum_hi; /* kr1 */
80  krvals[1] = Rate_arr2[indx_lo]*dum_lo + Rate_arr2[indx_hi]*dum_hi; /* kr2 */
81  krvals[2] = Rate_arr3[indx_lo]*dum_lo + Rate_arr3[indx_hi]*dum_hi; /* kr3 */
82  krvals[3] = Rate_arr4[indx_lo]*dum_lo + Rate_arr4[indx_hi]*dum_hi; /* kr4 */
83  krvals[4] = Rate_arr5[indx_lo]*dum_lo + Rate_arr5[indx_hi]*dum_hi; /* cr */
84  krvals[5] = Rate_arr6[indx_lo]*dum_lo + Rate_arr6[indx_hi]*dum_hi; /* ci */
85  krvals[6] = Emiss_arr1[indx_lo]*dum_lo + Emiss_arr1[indx_hi]*dum_hi; /* em[1] */
86  krvals[7] = Emiss_arr2[indx_lo]*dum_lo + Emiss_arr2[indx_hi]*dum_hi; /* em[2] */
87  krvals[8] = Emiss_arr3[indx_lo]*dum_lo + Emiss_arr3[indx_hi]*dum_hi; /* em[11] */
88  }
89  }
90  return;
91 }
#define MAX(a, b)
Definition: macros.h:101
tuple T
Definition: Sph_disk.py:33
#define CONST_eV
Electron Volt in erg.
Definition: pluto.h:256
tuple scrh
Definition: configure.py:200
#define MIN(a, b)
Definition: macros.h:104
#define TABLE_NPT
Definition: radiat.c:4
#define POW10(x)
Definition: radiat.c:3
#define INT_FLOOR(z)
Definition: macros.h:98
#define CONST_kB
Boltzmann constant.
Definition: pluto.h:259
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
int i
Definition: analysis.c:2
#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 NormalizeIons ( double *  )

Here is the caller graph for this function:

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 
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 }
#define MAX(a, b)
Definition: macros.h:101
tuple T
Definition: Sph_disk.py:33
#define NIONS_LOOP(n)
Definition: pluto.h:614
double g_gamma
Definition: globals.h:112
#define UNIT_DENSITY
Unit density in gr/cm^3.
Definition: pluto.h:369
#define T_CUT_RHOE
Sets the lowest cut-off temperature used in the PVTE_LAW equation of state.
Definition: eos.h:69
double InternalEnergy(double *v, double T)
double real
Definition: pluto.h:488
#define CONST_AH
Definition: pluto.h:250
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
#define RHO
Definition: mod_defs.h:19
#define CONST_amu
Atomic mass unit.
Definition: pluto.h:253
#define frac_He
Definition: cooling.h:17
#define frac_Z
Definition: cooling.h:15
double g_smallPressure
Small value for pressure fix.
Definition: globals.h:110
#define KELVIN
Definition: pluto.h:401
#define He_MASS_FRAC
Definition: pluto.h:551
#define MIN(a, b)
Definition: macros.h:104
#define UNIT_VELOCITY
Unit velocity in cm/sec.
Definition: pluto.h:377
#define X_HII
Definition: cooling.h:31
#define POW10(x)
Definition: radiat.c:3
int GetEV_Temperature(double rhoe, double *v, double *T)
#define UNIT_LENGTH
Unit Length in cm.
Definition: pluto.h:373
#define CONST_AZ
Mean atomic weight of heavy elements.
Definition: pluto.h:252
Definition: cooling.h:110
int k
Definition: analysis.c:2
double g_minCoolingTemp
The minimum temperature (in K) below which cooling is suppressed.
Definition: globals.h:106
double MeanMolecularWeight(double *V)
#define X_H2
Definition: cooling.h:30
#define CONST_AHe
Atomic weight of Helium.
Definition: pluto.h:251
void H2RateTables(double T, double *krvals)
Definition: radiat.c:7
#define H_MASS_FRAC
Definition: pluto.h:542
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the caller graph for this function: