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