PLUTO
Functions
cooling.h File Reference

Go to the source code of this file.

Functions

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

Function Documentation

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 }
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

Here is the call graph for this function: