cloudy trunk
Loading...
Searching...
No Matches
mole_drive.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*mole_drive - public routine, calls mole_solve to converge molecules */
4#include "cddefines.h"
5#include "taulines.h"
6#include "dense.h"
7#include "hmi.h"
8#include "conv.h"
9#include "doppvel.h"
10#include "thermal.h"
11#include "gammas.h"
12#include "iso.h"
13#include "opacity.h"
14#include "phycon.h"
15#include "physconst.h"
16#include "radius.h"
17#include "rfield.h"
18#include "secondaries.h"
19#include "timesc.h"
20#include "trace.h"
21#include "thermal.h"
22#include "called.h"
23#include "hcmap.h"
24#include "coolheavy.h"
25#include "mole.h"
26#include "mole_priv.h"
27#include "dynamics.h"
28#include "h2.h"
29#include "co.h"
30#include "ionbal.h"
31
32/* this is the error tolerance for reporting non-convergence */
33static const double MOLETOLER = 0.10;
34
35STATIC void mole_effects(void);
37STATIC void mole_ion_trim(void);
39
40/*mole_drive main driver for heavy molecular equilibrium routines */
41void mole_drive(void)
42{
43 DEBUG_ENTRY( "mole_drive()" );
44
45 mole_update_species_cache(); /* Update densities of species controlled outside the chemical network */
46
48
50
51 double error = mole_solve();
52
53 bool lgConverged = (error < MOLETOLER);
54
55 if( !lgConverged )
56 {
57 // should something be done with this?
58 }
59
61
63
64 return;
65}
66
68{
69 for (ChemAtomList::iterator atom = unresolved_atom_list.begin();
70 atom != unresolved_atom_list.end(); ++atom)
71 {
72 const long int nelem=(*atom)->el->Z-1;
73 if( !dense.lgElmtOn[nelem] )
74 continue;
75
76 for (long int ion=0;ion<nelem+2;ion++)
77 {
78 if ((*atom)->ipMl[ion] != -1)
79 {
80 if (dense.IonLow[nelem] > ion )
81 {
82 if( dense.xIonDense[nelem][ion] > ( ionbal.trimlo * dense.gas_phase[nelem] ) )
83 dense.IonLow[nelem] = ion;
84 else
85 {
86 dense.xIonDense[nelem][dense.IonLow[nelem]] += dense.xIonDense[nelem][ion];
87 dense.xIonDense[nelem][ion] = 0.;
88 }
89 }
90
91 if (dense.IonHigh[nelem] < ion )
92 {
93 if( dense.xIonDense[nelem][ion] > ( ionbal.trimhi * dense.gas_phase[nelem] ) )
94 dense.IonHigh[nelem] = ion;
95 else
96 {
97 dense.xIonDense[nelem][dense.IonHigh[nelem]] += dense.xIonDense[nelem][ion];
98 dense.xIonDense[nelem][ion] = 0.;
99 }
100 }
101 }
102 }
103 }
104}
105
107{
108 DEBUG_ENTRY( "mole_update_sources()" );
109
112}
113
115{
116 double
117 plte;
118
119 DEBUG_ENTRY( "mole_effects()" );
120
122
123 co.CODissHeat = (realnum)(mole.findrate("PHOTON,CO=>C,O")*1e-12);
124
125 thermal.heating[0][9] = co.CODissHeat;
126
127 /* total H2 - all forms */
128 hmi.H2_total = findspecieslocal("H2*")->den + findspecieslocal("H2")->den;
129 hmi.HD_total = findspecieslocal("HD")->den;
130 /* first guess at ortho and para densities */
131 h2.ortho_density = 0.75*hmi.H2_total;
132 h2.para_density = 0.25*hmi.H2_total;
133 {
134 hmi.H2_total_f = (realnum)hmi.H2_total;
135 h2.ortho_density_f = (realnum)h2.ortho_density;
136 h2.para_density_f = (realnum)h2.para_density;
137 }
138
139 /* save rate H2 is destroyed units s-1 */
140 /* >>chng 05 mar 18, TE, add terms -
141 total destruction rate is: dest_tot = n_H2g/n_H2tot * dest_H2g + n_H2s/n_H2tot * dest_H2s */
142 /* as reactions that change H2s to H2g and vice versa are not counted destruction processes, the terms c[ipMH2g][ipMH2s] *
143 and c[ipMH2s][ipMH2g], which have a different sign than [ipMH2g][ipMH2g] and [ipMH2s][ipMH2s], have to be added */
144 hmi.H2_rate_destroy = mole.sink_rate_tot("H2");
145
146 /* establish local timescale for H2 molecule destruction */
147 if(hmi.H2_rate_destroy > SMALLFLOAT )
148 {
149 /* units are now seconds */
150 timesc.time_H2_Dest_here = 1./hmi.H2_rate_destroy;
151 }
152 else
153 {
154 timesc.time_H2_Dest_here = 0.;
155 }
156
157 /* local timescale for H2 formation
158 * both grains and radiative attachment
159 * >>chng 08 mar 15 rjrw -- ratio formation rate to density of _H2_ so comparable with destruction rate,
160 * use full rates not rate coefficients so have correct units */
161 timesc.time_H2_Form_here = (mole.findrate("H,H,grn=>H2*,grn")+mole.findrate("H,H,grn=>H2,grn"));
162 /* timescale is inverse of this rate */
163 if( timesc.time_H2_Form_here > SMALLFLOAT )
164 {
165 /* units are now seconds */
166 timesc.time_H2_Form_here = hmi.H2_total/timesc.time_H2_Form_here;
167 }
168 else
169 {
170 timesc.time_H2_Form_here = 0.;
171 }
172
173
174 /* heating due to H2 dissociation */
175
176 /* H2 photodissociation heating, eqn A9 of Tielens & Hollenbach 1985a */
177 /*hmi.HeatH2Dish_TH85 = (float)(1.36e-23*findspecieslocal("H2")->den*h2esc*hmi.UV_Cont_rel2_Habing_TH85_depth);*/
178 /* >>chng 04 feb 07, more general to express heating in terms of the assumed
179 * photo rates - the 0.25 was obtained by inverting A8 & A9 of TH85 to find that
180 * there are 0.25 eV per dissociative pumping, ie, 10% of total
181 * this includes both H2g and H2s - TH85 say just ground but they include
182 * process for both H2 and H2s - as we did above - both must be in
183 * heating term */
184 /* >>chng 05 mar 11, TE, old had used H2_Solomon_dissoc_rate_used, which was only
185 * H2g. in regions where Solomon process is fast, H2s has a large population
186 * and the heating rate was underestimated. */
187 /* >>chng 05 jun 23,
188 * >>chng 05 dec 05, TE, modified to approximate the heating better for the
189 * new approximation */
190 /* >>chng 00 nov 25, explicitly break out this heat agent */
191 /* 2.6 eV of heat per deexcitation, consider difference
192 * between deexcitation (heat) and excitation (cooling) */
193 /* >>chng 04 jan 27, code moved here and simplified */
194 /* >>chng 05 jul 10, GS*/
195 /* average collisional rate for H2* to H2g calculated from big H2, GS*/
196
197 /* TH85 dissociation heating - this is ALWAYS defined for reference
198 * may be output for comparison with other rates*/
199 hmi.HeatH2Dish_TH85 = 0.25 * EN1EV *
200 (hmi.H2_Solomon_dissoc_rate_used_H2g * findspecieslocal("H2")->den +
201 hmi.H2_Solomon_dissoc_rate_used_H2s * findspecieslocal("H2*")->den);
202
203 /* TH85 deexcitation heating */
204 hmi.HeatH2Dexc_TH85 = ((mole.findrate("H2*,H2=>H2,H2") + mole.findrate("H2*,H=>H2,H"))
205 - (mole.findrate("H2,H2=>H2*,H2") + mole.findrate("H2,H=>H2*,H"))) * 4.17e-12;
206 /* this is derivative wrt temperature, only if counted as a cooling source */
207 hmi.deriv_HeatH2Dexc_TH85 = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_TH85)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
208
209 if( hmi.chH2_small_model_type == 'H' )
210 {
211 /* Burton et al. 1990 */
212 hmi.HeatH2Dish_BHT90 = 0.25 * EN1EV *
213 (hmi.H2_Solomon_dissoc_rate_used_H2g * findspecieslocal("H2")->den +
214 hmi.H2_Solomon_dissoc_rate_used_H2s * findspecieslocal("H2*")->den);
215
216 /* Burton et al. 1990 heating */
217 hmi.HeatH2Dexc_BHT90 = ((mole.findrate("H2*,H2=>H2,H2") + mole.findrate("H2*,H=>H2,H"))
218 - (mole.findrate("H2,H2=>H2*,H2") + mole.findrate("H2,H=>H2*,H"))) * 4.17e-12;
219 /* this is derivative wrt temperature, only if counted as a cooling source */
220 hmi.deriv_HeatH2Dexc_BHT90 = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_BHT90)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
221 }
222 else if( hmi.chH2_small_model_type == 'B')
223 {
224 /* Bertoldi & Draine */
225 hmi.HeatH2Dish_BD96 = 0.25 * EN1EV *
226 (hmi.H2_Solomon_dissoc_rate_used_H2g * findspecieslocal("H2")->den +
227 hmi.H2_Solomon_dissoc_rate_used_H2s * findspecieslocal("H2*")->den);
228 /* Bertoldi & Draine heating, same as TH85 */
229 hmi.HeatH2Dexc_BD96 = ((mole.findrate("H2*,H2=>H2,H2") + mole.findrate("H2*,H=>H2,H"))
230 - (mole.findrate("H2,H2=>H2*,H2") + mole.findrate("H2,H=>H2*,H"))) * 4.17e-12;
231 /* this is derivative wrt temperature, only if counted as a cooling source */
232 hmi.deriv_HeatH2Dexc_BD96 = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_BD96)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
233 }
234 else if(hmi.chH2_small_model_type == 'E')
235 {
236 /* heating due to dissociation of H2
237 * >>chng 05 oct 19, TE, define new approximation for the heating due to the destruction of H2
238 * use this approximation for the specified cloud parameters, otherwise
239 * use limiting cases for 1 <= G0, G0 >= 1e7, n >= 1e7, n <= 1 */
240
241 double log_density,
242 f1, f2,f3, f4, f5;
243 static double log_G0_face = -1;
244 static double k_f4;
245
246
247 /* test for G0
248 * this is a constant so only do it in zone 0 */
249 if( !nzone )
250 {
251 if(hmi.UV_Cont_rel2_Draine_DB96_face <= 1.)
252 {
253 log_G0_face = 0.;
254 }
255 else if(hmi.UV_Cont_rel2_Draine_DB96_face >= 1e7)
256 {
257 log_G0_face = 7.;
258 }
259 else
260 {
261 log_G0_face = log10(hmi.UV_Cont_rel2_Draine_DB96_face);
262 }
263 /*>>chng 06 oct 24 TE change Go face for spherical geometry*/
264 log_G0_face /= radius.r1r0sq;
265 }
266
267 /* test for density */
268 if(dense.gas_phase[ipHYDROGEN] <= 1.)
269 {
270 log_density = 0.;
271 }
272 else if(dense.gas_phase[ipHYDROGEN] >= 1e7)
273 {
274 log_density = 7.;
275 }
276 else
277 {
278 log_density = log10(dense.gas_phase[ipHYDROGEN]);
279 }
280
281 f1 = 0.15 * log_density + 0.75;
282 f2 = -0.5 * log_density + 10.;
283
284 hmi.HeatH2Dish_ELWERT = 0.25 * EN1EV * f1 *
285 (hmi.H2_Solomon_dissoc_rate_used_H2g * findspecieslocal("H2")->den +
286 hmi.H2_Solomon_dissoc_rate_used_H2s * findspecieslocal("H2*")->den ) +
287 f2 * secondaries.x12tot * EN1EV * hmi.H2_total;
288
289 /*fprintf( ioQQQ, "f1: %.2e, f2: %.2e,heat Solomon: %.2e",f1,f2,hmi.HeatH2Dish_TH85);*/
290
291
292 /* heating due to collisional deexcitation within X of H2
293 * >>chng 05 oct 19, TE, define new approximation for the heating due to the destruction of H2
294 * use this approximation for the specified cloud parameters, otherwise
295 * use limiting cases for 1 <= G0, G0 >= 1e7, n >= 1e7, n <= 1 */
296
297 /* set value of k_f4 by testing on value of G0 */
298 if(hmi.UV_Cont_rel2_Draine_DB96_face <= 1.)
299 {
300 log_G0_face = 0.;
301 }
302 else if(hmi.UV_Cont_rel2_Draine_DB96_face >= 1e7)
303 {
304 log_G0_face = 7.;
305 }
306 else
307 {
308 log_G0_face = log10(hmi.UV_Cont_rel2_Draine_DB96_face);
309 }
310 /* 06 oct 24, TE introduce effects of spherical geometry */
311 log_G0_face /= radius.r1r0sq;
312
313 /* terms only dependent on G0_face */
314 k_f4 = -0.25 * log_G0_face + 1.25;
315
316 /* test for density */
317 if(dense.gas_phase[ipHYDROGEN] <= 1.)
318 {
319 log_density = 0.;
320 f4 = 0.;
321 }
322 else if(dense.gas_phase[ipHYDROGEN] >= 1e7)
323 {
324 log_density = 7.;
325 f4 = pow2(k_f4) * pow( 10. , 2.2211 * log_density - 29.8506);
326 }
327 else
328 {
329 log_density = log10(dense.gas_phase[ipHYDROGEN]);
330 f4 = pow2(k_f4) * pow( 10., 2.2211 * log_density - 29.8506);
331 }
332
333 f3 = MAX2(0.1, -4.75 * log_density + 24.25);
334 f5 = MAX2(1.,0.95 * log_density - 1.45) * 0.2 * log_G0_face;
335
336 hmi.HeatH2Dexc_ELWERT = ((mole.findrate("H2*,H2=>H2,H2") + mole.findrate("H2*,H=>H2,H"))
337 - (mole.findrate("H2,H2=>H2*,H2") + mole.findrate("H2,H=>H2*,H"))) * 4.17e-12 * f3 +
338 f4 * (mole.species[unresolved_atom_list[ipHYDROGEN]->ipMl[0]].den/dense.gas_phase[ipHYDROGEN]) +
339 f5 * secondaries.x12tot * EN1EV * hmi.H2_total;
340
341 if(log_G0_face == 0.&& dense.gas_phase[ipHYDROGEN] > 1.)
342 hmi.HeatH2Dexc_ELWERT *= 0.001 / dense.gas_phase[ipHYDROGEN];
343
344 /* >>chng 06 oct 24, TE introduce effects of spherical geometry */
345 /*if(radius.depth/radius.rinner >= 1.0) */
346 hmi.HeatH2Dexc_ELWERT /= POW2(radius.r1r0sq);
347
348 /* this is derivative wrt temperature, only if counted as a cooling source */
349 hmi.deriv_HeatH2Dexc_ELWERT = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_ELWERT)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
350
351 /*fprintf( ioQQQ, "\tf3: %.2e, f4: %.2e, f5: %.2e, heat coll dissoc: %.2e\n",f3,f4,f5,hmi.HeatH2Dexc_TH85);*/
352 }
353 /* end Elwert branch for photo rates */
354 else
356
357
358 if( findspecieslocal("H-")->den > 0. && hmi.rel_pop_LTE_Hmin > 0. )
359 {
360 hmi.hmidep = (double)findspecieslocal("H-")->den/ SDIV(
361 ((double)dense.xIonDense[ipHYDROGEN][0]*dense.eden*hmi.rel_pop_LTE_Hmin));
362 }
363 else
364 {
365 hmi.hmidep = 1.;
366 }
367
368 /* this will be net volume heating rate, photo heat - induc cool */
369 hmi.hmihet = hmi.HMinus_photo_heat*findspecieslocal("H-")->den - hmi.HMinus_induc_rec_cooling;
370
371 /* H2+ + HNU => H+ + H */
372 t_phoHeat photoHeat;
373 (void) GammaK(opac.ih2pnt[0],opac.ih2pnt[1],opac.ih2pof,1.,&photoHeat); /* Now used entirely for side-effect. Yuk */
374 hmi.h2plus_heat = photoHeat.HeatNet*findspecieslocal("H2+")->den;
375
376 /* departure coefficient for H2 defined rel to n(1s)**2
377 * (see Littes and Mihalas Solar Phys 93, 23) */
378 plte = (double)dense.xIonDense[ipHYDROGEN][0] * hmi.rel_pop_LTE_H2g * (double)dense.xIonDense[ipHYDROGEN][0];
379 if( plte > 0. )
380 {
381 hmi.h2dep = findspecieslocal("H2")->den/plte;
382 }
383 else
384 {
385 hmi.h2dep = 1.;
386 }
387
388 /* departure coefficient of H2+ defined rel to n(1s) n(p)
389 * sec den was HI before 85.34 */
390 plte = (double)dense.xIonDense[ipHYDROGEN][0]*hmi.rel_pop_LTE_H2p*(double)dense.xIonDense[ipHYDROGEN][1];
391 if( plte > 0. )
392 {
393 hmi.h2pdep = findspecieslocal("H2+")->den/plte;
394 }
395 else
396 {
397 hmi.h2pdep = 1.;
398 }
399
400 /* departure coefficient of H3+ defined rel to N(H2+) N(p) */
401 if( hmi.rel_pop_LTE_H3p > 0. )
402 {
403 hmi.h3pdep = findspecieslocal("H3+")->den/hmi.rel_pop_LTE_H3p;
404 }
405 else
406 {
407 hmi.h3pdep = 1.;
408 }
409
411
412 return;
413}
414#if defined(__HP_aCC)
415#pragma OPTIMIZE OFF
416#pragma OPTIMIZE ON
417#endif
419{
420 int i;
421 double
422 rate;
423
424 /* identify dominant H2 formation process */
425 {
426 /* following should be set true to identify H2 formation and destruction processes */
427 /*@-redef@*/
428 enum {DEBUG_LOC=false};
429 /*@+redef@*/
430 if( DEBUG_LOC && (nzone>50) )
431 {
432 double createsum ,create_from_Hn2 , create_3body_Ho, create_h2p,
433 create_h3p, create_h3pe, create_grains, create_hminus;
434 double destroysum, destroy_hm ,destroy_solomon ,destroy_2h ,destroy_hp,
435 destroy_h,destroy_hp2,destroy_h3p;
436
437 create_from_Hn2 = mole.findrate("H,H=>H2"); /* H(n=2) + H(n=1) -> H2 */
438 create_3body_Ho = mole.findrate("H,H,H=>H2,H");
439 create_h2p = mole.findrate("H,H2+=>H+,H2*");
440 create_h3p = mole.findrate("H,H3+=>H2,H2+");
441 create_h3pe = mole.findrate("e-,H3+=>H2,H");
442 create_grains = mole.findrate("H,H,grn=>H2*,grn")+mole.findrate("H,H,grn=>H2,grn");
443 create_hminus = (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-"));
444
445 createsum = create_from_Hn2 + create_3body_Ho + create_h2p +
446 create_h3p + create_h3pe + create_grains + create_hminus;
447
448 fprintf(ioQQQ,"H2 create zone\t%.2f \tsum\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
449 fnzone,
450 createsum,
451 create_hminus / createsum,
452 create_from_Hn2 / createsum,
453 create_3body_Ho / createsum,
454 create_h2p / createsum,
455 create_h3p / createsum,
456 create_h3pe / createsum,
457 create_grains / createsum );
458
459 /* all h2 molecular hydrogen destruction processes */
460 /* >>chng 04 jan 28, had wrong Boltzmann factor,
461 * caught by gargi Shaw */
462 destroy_hm = (mole.findrate("H2,e-=>H,H-")+mole.findrate("H2*,e-=>H,H-"));
463 /*destroy_hm2 = eh2hhm;*/
464 destroy_solomon = hmi.H2_Solomon_dissoc_rate_used_H2g * findspecieslocal("H2")->den;
465 destroy_2h = (mole.findrate("H2,e-=>H,H,e-")+mole.findrate("H2*,e-=>H,H,e-"));
466 destroy_hp = mole.findrate("H2,H+=>H3+");
467 destroy_h = mole.findrate("H,H2=>H,H,H");
468 destroy_hp2 = mole.findrate("H2,H+=>H2+,H");
469 destroy_h3p = mole.findrate("H2,H3+=>H2,H2+,H");
470 destroysum = destroy_hm + /*destroy_hm2 +*/ destroy_solomon + destroy_2h +
471 destroy_hp+ destroy_h+ destroy_hp2+ destroy_h3p;
472
473 fprintf(ioQQQ,"H2 destroy\t%.3f \t%.2e\tsum\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
474 fnzone,
475 destroysum,
476 destroy_hm / destroysum ,
477 destroy_solomon / destroysum ,
478 destroy_2h / destroysum ,
479 destroy_hp / destroysum ,
480 destroy_h / destroysum ,
481 destroy_hp2 / destroysum ,
482 destroy_h3p / destroysum );
483
484 }
485 }
486
487 {
488 /* following should be set true to identify H- formation and destruction processes */
489 /*@-redef@*/
490 enum {DEBUG_LOC=false};
491 /*@+redef@*/
492 if( DEBUG_LOC && (nzone>140)/**/ )
493 {
494 double create_from_Ho,create_3body_Ho,create_batach,destroy_photo,
495 destroy_coll_heavies,destroy_coll_electrons,destroy_Hattach,destroy_fhneut,
496 destsum , createsum;
497
498 create_from_Ho = mole.findrate("H,e-=>H-,PHOTON");
499 create_3body_Ho = mole.findrate("H,e-,e-=>H-,e-");
500 /* total formation is sum of g and s attachment */
501 create_batach = (mole.findrate("H2,e-=>H,H-") + mole.findrate("H2*,e-=>H,H-"));
502 destroy_photo = mole.findrate("H-,PHOTON=>H,e-");
503 destroy_coll_heavies = 0.;
504 for( i=ipLITHIUM; i < LIMELM; i++ )
505 {
506 char react[32];
507 if (dense.lgElmtOn[i])
508 {
509 sprintf(react,"H-,%s=>H,%s",mole_global.list[unresolved_atom_list[i]->ipMl[1]]->label.c_str(),unresolved_atom_list[i]->label().c_str() );
510 destroy_coll_heavies += mole.findrate(react);
511 }
512 }
513 destroy_coll_electrons = mole.findrate("H-,e-=>H-,e-,e-");
514 destroy_Hattach = (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-"));
515 destroy_fhneut = mole.findrate("H-,H+=>H,H");
516
517 destsum = destroy_photo + destroy_coll_heavies + destroy_coll_electrons +
518 destroy_Hattach + destroy_fhneut;
519 fprintf(ioQQQ,"H- destroy zone\t%.2f\tTe\t%.4e\tsum\t%.2e\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n",
520 fnzone,
521 phycon.te,
522 destsum,
523 destroy_photo/destsum ,
524 destroy_coll_heavies/destsum,
525 destroy_coll_electrons/destsum,
526 destroy_Hattach/destsum,
527 destroy_fhneut/destsum );
528
529 createsum = create_from_Ho+create_3body_Ho+create_batach;
530 fprintf(ioQQQ,"H- create\t%.2f\tTe\t%.4e\tsum\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
531 fnzone,
532 phycon.te,
533 createsum,
534 dense.eden,
535 create_from_Ho/createsum,
536 create_3body_Ho/createsum,
537 create_batach/createsum);
538 }
539 }
540 {
541 /* following should be set true to print populations */
542 /*@-redef@*/
543 enum {DEBUG_LOC=false};
544 /*@+redef@*/
545 if( DEBUG_LOC )
546 {
547 /* these are the raw results */
548 fprintf( ioQQQ, " MOLE raw; hi\t%.2e" , dense.xIonDense[ipHYDROGEN][0]);
549 for( i=0; i < mole_global.num_calc; i++ )
550 {
551 fprintf( ioQQQ, "\t%-4.4s\t%.2e", mole_global.list[i]->label.c_str(), mole.species[i].den );
552 }
553 fprintf( ioQQQ, " \n" );
554 }
555 }
556
557 if( trace.lgTrace && trace.lgTraceMole )
558 {
559 /* these are the raw results */
560 fprintf( ioQQQ, " raw; " );
561 for( i=0; i < mole_global.num_calc; i++ )
562 {
563 fprintf( ioQQQ, " %-4.4s:%.2e", mole_global.list[i]->label.c_str(), mole.species[i].den );
564 }
565 fprintf( ioQQQ, " \n" );
566
567 }
568
569 /* option to print rate H2 forms */
570 /* trace.lgTraceMole is trace molecules option,
571 * punch htwo */
572 if( (trace.lgTrace && trace.lgTraceMole) )
573 {
575 double H2_rate_create = mole.source_rate_tot("H2") + mole.source_rate_tot("H2*");
576
577 if( H2_rate_create > SMALLFLOAT )
578 {
579 fprintf( ioQQQ, " Create H2, rate=%10.2e grain;%5.3f hmin;%5.3f bhedis;%5.3f h2+;%5.3f hmi.radasc:%5.3f hmi.h3ph2p:%5.3f hmi.h3petc:%5.3f\n",
580 H2_rate_create,
581 (mole.findrate("H,H,grn=>H2*,grn")+mole.findrate("H,H,grn=>H2,grn"))/H2_rate_create,
582 (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-"))/H2_rate_create,
583 mole.findrate("H,H,H=>H2,H")/H2_rate_create,
584 mole.findrate("H,H2+=>H+,H2*")/H2_rate_create,
585 mole.findrate("H,H=>H2")/H2_rate_create,
586 mole.findrate("H,H3+=>H2,H2+")/H2_rate_create,
587 mole.findrate("H2,H3+=>H2,H2+,H")/H2_rate_create );
588 }
589 else
590 {
591 fprintf( ioQQQ, " Create H2, rate=0\n" );
592 }
593 }
594
595 /* this is H2+ */
596 if( trace.lgTrace && trace.lgTraceMole )
597 {
598 rate = mole.findrate("H2,H+=>H2+,H") + mole.findrate("H,H+,e-=>H2+,e-") +
599 mole.findrate("H,H3+=>H2,H2+") + mole.findrate("H2,H3+=>H2,H2+,H");
600 if( rate > 1e-25 )
601 {
602 fprintf( ioQQQ, " Create H2+, rate=%10.2e hmi.rh2h2p;%5.3f b2pcin;%5.3f hmi.h3ph2p;%5.3f hmi.h3petc+;%5.3f\n",
603 rate, mole.findrate("H2,H+=>H2+,H")/rate,
604 mole.findrate("H,H+,e-=>H2+,e-")/rate, mole.findrate("H,H3+=>H2,H2+")/rate, mole.findrate("H2,H3+=>H2,H2+,H")/rate );
605 }
606 else
607 {
608 fprintf( ioQQQ, " Create H2+, rate=0\n" );
609 }
610 }
611
612 if( trace.lgTrace && trace.lgTraceMole )
613 {
614
615 double destroy_coll_heavies = 0.;
616
617 for( i=ipLITHIUM; i < LIMELM; i++ )
618 {
619 char react[32];
620 if (dense.lgElmtOn[i])
621 {
622 sprintf(react,"H-,%s=>H,%s",mole_global.list[unresolved_atom_list[i]->ipMl[1]]->label.c_str(),unresolved_atom_list[i]->label().c_str() );
623 destroy_coll_heavies += mole.findrate(react);
624 }
625 }
626 fprintf( ioQQQ, " MOLE, Dep Coef, H-:%10.2e H2:%10.2e H2+:%10.2e\n",
627 hmi.hmidep, hmi.h2dep, hmi.h2pdep );
628 fprintf( ioQQQ, " H- creat: Rad atch%10.3e Induc%10.3e bHneut%10.2e 3bod%10.2e b=>H2%10.2e N(H-);%10.2e b(H-);%10.2e\n",
629 mole.findrate("H,e-=>H-,PHOTON"), hmi.HMinus_induc_rec_rate*findspecieslocal("H-")->den, mole.findrate("H,H=>H-,H+"),
630 mole.findrate("H,e-,e-=>H-,e-"),
631 mole.findrate("H2,e-=>H,H-"), findspecieslocal("H-")->den, hmi.hmidep );
632
633 fprintf( ioQQQ, " H- destr: Photo;%10.3e mut neut%10.2e e- coll ion%10.2e =>H2%10.2e x-ray%10.2e p+H-%10.2e\n",
634 mole.findrate("H-,PHOTON=>H,e-"), destroy_coll_heavies,
635 mole.findrate("H-,e-=>H-,e-,e-"),
636 (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-")), iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*findspecieslocal("H-")->den,
637 mole.findrate("H-,H+=>H,H") );
638 fprintf( ioQQQ, " H- heating:%10.3e Ind cooling%10.2e Spon cooling%10.2e\n",
639 hmi.hmihet, hmi.HMinus_induc_rec_cooling, hmi.hmicol );
640 }
641
642 /* identify creation and destruction processes for H2+ */
643 if( trace.lgTrace && trace.lgTraceMole )
644 {
645 rate = findspecieslocal("H2+")->snk;
646 if( rate != 0. )
647 {
648 rate *= findspecieslocal("H2+")->den;
649 if( rate > 0. )
650 {
651 fprintf( ioQQQ,
652 " Destroy H2+: rate=%10.2e e-;%5.3f phot;%5.3f hard gam;%5.3f H2col;%5.3f h2phhp;%5.3f pion;%5.3f bh2h2p:%5.3f\n",
653 rate, mole.findrate("H2+,e-=>H,H+,e-")/rate, mole.findrate("H2+,PHOTON=>H,H+")/rate,
654 mole.findrate("H2+,CRPHOT=>H,H+")/rate, mole.findrate("H2,H2+=>H,H3+")/rate, mole.findrate("H2+,H2=>H,H+,H2")/rate,
655 mole.findrate("H2+,H+=>H,H+,H+")/rate, mole.findrate("H,H2+=>H+,H2*")/rate );
656
657 fprintf( ioQQQ,
658 " Create H2+: rate=%.2e HII HI;%.3f Col H2;%.3f HII H2;%.3f HI HI;%.3f\n",
659 rate,
660 mole.findrate("H+,H=>H2+,PHOTON")/rate,
661 mole.findrate("H2,CRPHOT=>H2+,e-")/rate,
662 mole.findrate("H2,H+=>H2+,H")/rate,
663 mole.findrate("H,H+,e-=>H2+,e-")/rate );
664 }
665 else
666 {
667 fprintf( ioQQQ, " Create H2+: rate= is zero\n" );
668 }
669 }
670 }
671
672 {
673 /* following should be set true to print populations */
674 /*@-redef@*/
675 enum {DEBUG_LOC=false};
676 /*@+redef@*/
677 if( DEBUG_LOC )
678 {
679 fprintf(ioQQQ,"mole bugg\t%.3f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
680 fnzone,
681 iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].gamnc,
682 hmi.HMinus_photo_rate,
683 findspecieslocal("H2")->den ,
684 findspecieslocal("H-")->den ,
685 findspecieslocal("H+")->den);
686 }
687 }
688
689 {
690 /*@-redef@*/
691 /* this debug print statement compares H2 formation through grn vs H- */
692 enum {DEBUG_LOC=false};
693 /*@+redef@*/
694 if( DEBUG_LOC && nzone>140 )
695 {
696 fprintf(ioQQQ," debuggggrn grn\t%.2f\t%.3e\t%.3e\tfrac\t%.3e\tH-\t%.3e\t%.3e\tfrac\t%.3e\t%.3e\t%.3e\t%.3e\n",
697 fnzone ,
698 mole.findrate("H,H,grn=>H2*,grn")+mole.findrate("H,H,grn=>H2,grn") ,
699 hmi.H2_forms_grains+hmi.H2star_forms_grains ,
700 mole.findrate("H,H,grn=>H2*,grn")/SDIV(mole.findrate("H,H,grn=>H2*,grn")+mole.findrate("H,H,grn=>H2,grn")),
701 (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-")),
702 hmi.H2star_forms_hminus+hmi.H2_forms_hminus,
704 (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-")),findspecieslocal("H")->den,findspecieslocal("H-")->den
705 );
706 }
707 }
708
709 {
710 /*@-redef@*/
711 /* often the H- route is the most efficient formation mechanism for H2,
712 * will be through rate called mole_global.list[unresolved_atom_list[ipHYDROGEN]->ipMl]->den*hmi.assoc_detach
713 * this debug print statement is to trace h2 oscillations */
714 enum {DEBUG_LOC=false};
715 /*@+redef@*/
716 if( DEBUG_LOC && nzone>140/*&& iteration > 1*/)
717 {
718 /* rapid increase in H2 density caused by rapid increase in hmi.rel_pop_LTE_H2g */
719 fprintf(ioQQQ,
720 "hmi.assoc_detach_backwards_grnd\t%.2f\t%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
721 /* total forward rate */
722 fnzone,
723 phycon.te,
724 dense.eden,
725 (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-")),
726 mole.findrate("H2,e-=>H,H-"),
727 mole.findrate("H2*,e-=>H,H-"),
728 mole.findrate("H,e-=>H-,PHOTON"),
729 hmi.HMinus_induc_rec_rate*findspecieslocal("H-")->den,
730 mole.species[unresolved_atom_list[ipHYDROGEN]->ipMl[0]].den,
731 mole.species[unresolved_atom_list[ipHYDROGEN]->ipMl[1]].den,
732 findspecieslocal("H-")->den,
733 hmi.H2_total,
734 hmi.rel_pop_LTE_Hmin,
735 hmi.rel_pop_LTE_H2g,
736 hmi.rel_pop_LTE_H2s
737 );
738 }
739 }
740
741 if( (trace.lgTrace && trace.lgTraceMole) )
742 {
743 if( hmi.H2_rate_destroy > SMALLFLOAT )
744 {
745 fprintf( ioQQQ,
746 " H2 destroy rate=%.2e DIS;%.3f bat;%.3f h2dis;%.3f photoionize_rate;%.3f h2h2p;%.3f E-h;%.3f hmi.h2hph3p;%.3f sec;%.3f\n",
747 hmi.H2_rate_destroy,
748 hmi.H2_Solomon_dissoc_rate_used_H2g / hmi.H2_rate_destroy,
749 mole.findrate("H2,e-=>H,H-") / (hmi.H2_total*hmi.H2_rate_destroy),
750 mole.findrate("H,H2=>H,H,H") / (hmi.H2_total*hmi.H2_rate_destroy),
751 h2.photoionize_rate / hmi.H2_rate_destroy,
752 mole.findrate("H2,H+=>H2+,H") / (hmi.H2_total* hmi.H2_rate_destroy),
753 (mole.findrate("H2,e-=>H,H,e-")+mole.findrate("H2*,e-=>H,H,e-")) / (hmi.H2_total* hmi.H2_rate_destroy),
754 mole.findrate("H2,H+=>H3+") / (hmi.H2_total* hmi.H2_rate_destroy) ,
755 secondaries.csupra[ipHYDROGEN][0]*2.02 / hmi.H2_rate_destroy
756 );
757 }
758 else
759 {
760 fprintf( ioQQQ, " Destroy H2: rate=0\n" );
761 }
762 }
763
764 {
765 /* following should be set true to print populations */
766 /*@-redef@*/
767 enum {DEBUG_LOC=false};
768 /*@+redef@*/
769 if( DEBUG_LOC )
770 {
771 if( DEBUG_LOC && (nzone > 570) )
772 {
773 fprintf(ioQQQ,"Temperature %g\n",phycon.te);
774 fprintf(ioQQQ," Net mol ion rate [%g %g] %g\n",mole.source[ipHYDROGEN][1],mole.sink[ipHYDROGEN][1],
775 mole.source[ipHYDROGEN][1]-mole.sink[ipHYDROGEN][1]*mole.species[unresolved_atom_list[ipHYDROGEN]->ipMl[1]].den);
776 }
777 }
778 }
779}
780
782{
783 DEBUG_ENTRY( "mole_update_limiting_reactants()" );
784
785 /* largest fraction of atoms in molecules */
786 for( long i=0; i<mole_global.num_calc; ++i )
787 {
788 mole.species[i].xFracLim = 0.;
789 mole.species[i].nAtomLim = -1;
790 for (molecule::nAtomsMap::iterator atom = mole_global.list[i]->nAtom.begin();
791 atom != mole_global.list[i]->nAtom.end(); ++atom)
792 {
793 long nelem = atom->first->el->Z-1;
794 if( dense.lgElmtOn[nelem])
795 {
796 realnum dens_elemsp = (realnum)( mole.species[i].den * atom->second * atom->first->frac );
797 if (mole.species[i].xFracLim*dense.gas_phase[nelem] < dens_elemsp )
798 {
799 mole.species[i].xFracLim = dens_elemsp/dense.gas_phase[nelem];
800 mole.species[i].nAtomLim = nelem;
801 }
802 }
803 }
804 //fprintf( ioQQQ, "DEBUGGG mole_update_limiting_reactants %li\t%s\t%i\t%.12e\n", i, mole_global.list[i]->label.c_str(), mole.species[i].nAtomLim, mole.species[i].xFracLim );
805 }
806
807 return;
808}
809
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
double fnzone
Definition cddefines.cpp:15
#define STATIC
Definition cddefines.h:97
T pow2(T a)
Definition cddefines.h:931
const int LIMELM
Definition cddefines.h:258
#define POW2
Definition cddefines.h:929
const int ipLITHIUM
Definition cddefines.h:307
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
NORETURN void TotalInsanity(void)
Definition service.cpp:886
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
double snk
Definition mole.h:352
double den
Definition mole.h:358
t_co co
Definition co.cpp:5
double GammaK(long int ipLoEnr, long int ipHiEnr, long int ipOpac, double yield1, t_phoHeat *photoHeat)
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
t_hmi hmi
Definition hmi.cpp:5
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH1s
Definition iso.h:27
const int ipH_LIKE
Definition iso.h:62
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
molezone * findspecieslocal(const char buf[])
void mole_update_species_cache(void)
ChemAtomList unresolved_atom_list
STATIC void mole_h_rate_diagnostics(void)
void mole_drive(void)
STATIC void mole_ion_trim(void)
void mole_update_sources(void)
STATIC void mole_update_limiting_reactants()
STATIC void mole_effects(void)
static const double MOLETOLER
void mole_eval_sources(long int num_total)
double mole_solve(void)
double frac_H2star_hminus()
void mole_update_rks(void)
t_opac opac
Definition opacity.cpp:5
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double EN1EV
Definition physconst.h:192
t_radius radius
Definition radius.cpp:5
t_secondaries secondaries
double HeatNet
Definition thermal.h:172
t_thermal thermal
Definition thermal.cpp:5
t_timesc timesc
Definition timesc.cpp:5
t_trace trace
Definition trace.cpp:5