cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mole_co_drive.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*lgMolecAver average old and new molecular equilibrium balance from CO_solve */
4 /*CO_drive - public routine, calls CO_solve to converge molecules */
5 #include "cddefines.h"
6 #include "taulines.h"
7 #include "dense.h"
8 #include "hmi.h"
9 #include "conv.h"
10 #include "phycon.h"
11 #include "trace.h"
12 #include "thermal.h"
13 #include "called.h"
14 #include "hcmap.h"
15 #include "coolheavy.h"
16 #include "mole.h"
17 #include "mole_co_priv.h"
18 
19 /* says whether the co network is currently set to zero */
20 static bool lgMoleZeroed=true;
21 
22 /* the limit to H2/H where we will solve for CO */
23 static double h2lim;
24 
25 /* this is the relative change in OH, which is used as a convergence
26  * criteria in lgMolecAver */
27 static const double COTOLER_MOLAV = 0.10;
28 
29 /* flag to control debug statement, if 0 then none, just loops when 1, more when 2 */
30 static const int CODEBUG = 0;
31 /* this is the maximum number of loops CO_drive will go over while
32  * trying to converge the molecular network */
33 static const int LUPMAX_CODRIV = 100;
34 
35 /*lgMolecAver average old and new molecular equilibrium balance from CO_solve,
36  * returns true if converged */
37 STATIC bool lgMolecAver(
38  /* job == "SETUP", set things up during first call,
39  * job == "AVER" take average of arrays */
40  const char chJob[10] ,
41  char *chReason );
42 
43 /*CO_drive main driver for heavy molecular equilibrium routines */
44 void CO_drive(void)
45 {
46  bool lgNegPop,
47  lgZerPop,
48  lgFirstTry;
49  long int i;
50  bool lgConverged;
51  /* count how many times through loop */
52  long int loop;
53  long int nelem , ion;
54 
55 
56  /* this will give reason CO not converged */
57  char chReason[100];
58  /*static bool lgMoleZeroed=true;*/
59 
60  DEBUG_ENTRY( "CO_drive()" );
61 
62  /* h2lim is smallest ratio of H2/HDEN for which we will
63  * even try to invert heavy element network */
64  if( hmi.lgNoH2Mole )
65  {
66  /* H2 network is turned off, ignore this limit */
67  h2lim = 0.;
68  }
69  else
70  {
71  /* this limit is important since the co molecular network is first turned
72  * on when this is hit. It's conservation law will only ever include the
73  * initial O, O+, and C, C+, so must not start before nearly all
74  * C and O is in these forms */
75  h2lim = 1e-15;
76  /* >>chng 05 jul 15, from 1e-15 to 1e-12, see whether results are stable,
77  * this does include CO in the H+ zone in orion_hii_pdr
78  * a problem is that, during search phase, the first temp is 4000K and the
79  * H2 abundance is larger than it will be at the illuminated face. try to
80  * avoid turning H2 on too soon */
81  h2lim = 1e-12;
82  }
83 
84  /* do we want to try to calculate the CO?
85  * >>chng 03 sep 07, add first test, so that we do not turn off
86  * heavy element network if it has ever been turned on
87  * >>chng 04 apr 23, change logic so never done when temp > 20000K */
88  /* make series of tests setting co.lgCODoCalc for simplicity */
89  co.lgCODoCalc = true;
90  /* too hot - don't bother Te > 2e4K */
91  if( phycon.te > 20000. )
92  co.lgCODoCalc = false;
93 
94  /* molecules have been turned off */
95  if( co.lgNoCOMole )
96  co.lgCODoCalc = false;
97 
98  /* a critical element has been turned off */
100  co.lgCODoCalc = false;
101 
102  /* >>chng 04 jul 20, do not attempt co if no O or C atoms present,
103  * since co net needs their atomic data */
104  if( dense.IonLow[ipCARBON]>0 || dense.IonLow[ipOXYGEN]>0 )
105  co.lgCODoCalc = false;
106 
107  /* do not do if H2/Htot < h2lim which is set to 1e-12 by default above */
108  if( iteration!=co.iteration_co &&
110  co.lgCODoCalc = false;
111 
112  /* >>chng 06 sep 10, do not call CO network if H+/Htot is greater than 0.5 -
113  * this change suggested by Nick Abel after cooling plasma developed CO
114  * convergence problems when H+/H = 0.97 */
116  co.lgCODoCalc = false;
117 
118  if( dense.lgElmtOn[ipSILICON] )
119  {
120  /* >>chng 04 jun 28, add test on atomic Si, some H II region
121  * models crashed due to matrix failure when co network turned
122  * on with atomic Si at very small abundance */
123  if( dense.xIonDense[ipSILICON][0]/dense.gas_phase[ipSILICON] < 1e-15 )
124  co.lgCODoCalc = false;
125  }
126 
127  /* this branch we will not do the CO equilibrium, set some variables that
128  * would be calculated by the chemistry, zero others, and return */
129  if( !co.lgCODoCalc )
130  {
131  /* these are heavy - heavy charge transfer rate that will still be needed
132  * for recombination of Si+, S+, and C+ */
133  struct COmole_rate_s *rate;
134 
135  mole.xMoleChTrRate[ipSILICON][1][0] = 0.;
136  rate = CO_findrate_s("Si+,Fe=>Si,Fe+");
137  mole.xMoleChTrRate[ipSILICON][1][0] += (realnum)(rate->rk*rate->rate_species[1]->hevmol);
138  rate = CO_findrate_s("Si+,Mg=>Si,Mg+");
139  mole.xMoleChTrRate[ipSILICON][1][0] += (realnum)(rate->rk*rate->rate_species[1]->hevmol);
140 
141  mole.xMoleChTrRate[ipSULPHUR][1][0] = 0.;
142  rate = CO_findrate_s("S+,Fe=>S,Fe+");
143  mole.xMoleChTrRate[ipSULPHUR][1][0] += (realnum)(rate->rk*rate->rate_species[1]->hevmol);
144  rate = CO_findrate_s("S+,Mg=>S,Mg+");
145  mole.xMoleChTrRate[ipSULPHUR][1][0] += (realnum)(rate->rk*rate->rate_species[1]->hevmol);
146 
147  mole.xMoleChTrRate[ipCARBON][1][0] = 0.;
148  rate = CO_findrate_s("C+,Fe=>C,Fe+");
149  mole.xMoleChTrRate[ipCARBON][1][0] += (realnum)(rate->rk*rate->rate_species[1]->hevmol);
150  rate = CO_findrate_s("C+,Mg=>C,Mg+");
151  mole.xMoleChTrRate[ipCARBON][1][0] += (realnum)(rate->rk*rate->rate_species[1]->hevmol);
152 
153 #if 0
155  COmole_rate[ipR_SiP_Mg_Si_MgP].rk +
156  dense.xIonDense[ipIRON][0]*COmole_rate[ipR_SiP_Fe_Si_FeP].rk);
157 
159  COmole_rate[ipR_SP_Mg_S_MgP].rk +
160  dense.xIonDense[ipIRON][0]*COmole_rate[ipR_SP_Fe_S_FeP].rk);
161 
163  COmole_rate[ipR_CP_Mg_C_MgP].rk +
164  dense.xIonDense[ipIRON][0]*COmole_rate[ipR_CP_Fe_C_FeP].rk);
165 #endif
166 
167  /* do we need to zero out the arrays and vars? */
168  if( !lgMoleZeroed )
169  {
170  lgMoleZeroed = true;
171  for( i=0; i<mole.num_comole_calc; ++i )
172  {
173  COmole[i]->hevmol = 0.;
174  }
175  /* heating due to CO photodissociation */
176  thermal.heating[0][9] = 0.;
177  /* CO cooling */
178  CoolHeavy.C12O16Rot = 0.;
179  CoolHeavy.dC12O16Rot = 0.;
180  CoolHeavy.C13O16Rot = 0.;
181  CoolHeavy.dC13O16Rot = 0.;
182  co.CODissHeat = 0.;
183 
185  for( i=0; i < nCORotate; i++ )
186  {
187  C12O16Rotate[i].Lo->Pop = 0.;
188  C13O16Rotate[i].Lo->Pop = 0.;
189  /* population of upper level */
190  C12O16Rotate[i].Hi->Pop = 0.;
191  C13O16Rotate[i].Hi->Pop = 0.;
192  /* population of lower level with correction for stimulated emission */
193  C12O16Rotate[i].Emis->PopOpc = 0.;
194  C13O16Rotate[i].Emis->PopOpc = 0.;
195  /* following two heat exchange excitation, deexcitation */
196  C12O16Rotate[i].Coll.cool = 0.;
197  C12O16Rotate[i].Coll.heat = 0.;
198  C13O16Rotate[i].Coll.cool = 0.;
199  C13O16Rotate[i].Coll.heat = 0.;
200  /* intensity of line */
201  C12O16Rotate[i].Emis->xIntensity = 0.;
202  C13O16Rotate[i].Emis->xIntensity = 0.;
203  /* number of photons emitted in line */
204  C12O16Rotate[i].Emis->phots = 0.;
205  C13O16Rotate[i].Emis->phots = 0.;
206  }
207  for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem )
208  {
209  for( ion=0; ion<nelem+2; ++ion )
210  {
211  mole.source[nelem][ion] = 0.;
212  mole.sink[nelem][ion] = 0.;
213  }
214  }
215  }
216 
217  if( trace.nTrConvg>=4 )
218  {
219  /* trace ionization */
220  fprintf( ioQQQ,
221  " CO_drive4 do not evaluate CO chemistry.\n");
222  }
223  return;
224  }
225 
226  CO_update_species_cache(); /* Update densities of species controlled outside the network */
227 
228  /* Update the rate constants -- final generic version */
229  CO_update_rks();
230 
231  /* this flag is used to stop calculation due to negative abundances */
232  loop = 0;
233  lgNegPop = false;
234  lgConverged = lgMolecAver("SETUP",chReason);
235  do
236  {
237 
238  if( trace.nTrConvg>=5 )
239  {
240  /* trace ionization */
241  fprintf( ioQQQ,
242  " CO_drive5 CO chemistry loop %li chReason %s.\n" , loop, chReason );
243  }
244 
245  /*oh_h_o_h2ld = findspecies("OH")->hevmol;*/
246  CO_solve(&lgNegPop,&lgZerPop );
247  /* this one takes average first, then updates molecular densities in co.hevmol,
248  * returns true if change in OH density is less than macro COTOLER_MOLAV set above*/
249  lgConverged = lgMolecAver("AVER",chReason);
250  if( CODEBUG )
251  fprintf(ioQQQ,"codrivbug %.2f %li Neg?:%c\tZero?:%c\tOH new\t%.3e\tCO\t%.3e\tTe\t%.3e\tH2/H\t%.2e\n",
252  fnzone ,
253  loop ,
254  TorF(lgNegPop) ,
255  TorF(lgZerPop) ,
256  findspecies("OH")->hevmol ,
257  findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]),
258  phycon.te,
260 
261  ++loop;
262  } while (
263  /* these are a series of conditions to stop this loop -
264  * has the limit to the number of loops been reached? */
265  (loop < LUPMAX_CODRIV) && !lgConverged );
266 
267  /* say that we have found a solution before */
268  lgMoleZeroed = false;
269 
270  /* check whether we have converged */
271  if( loop == LUPMAX_CODRIV)
272  {
273  conv.lgConvIoniz = false;
274  strcpy( conv.chConvIoniz, "CO con1" );
275  if( CODEBUG )
276  fprintf(ioQQQ,"CO_drive not converged\n");
277  }
278 
279  /* this flag says whether this is the first zone we are trying
280  * to converge the molecules - there will be problems during initial
281  * search so do not print anything in this case */
282  lgFirstTry = (nzone==co.co_nzone && iteration==co.iteration_co);
283 
284  /* did the molecule network have negative pops? */
285  if( lgNegPop )
286  {
288  (iteration==1) )
289  {
290  /* we are in search phase during the first iteration,
291  * and the H2/H ratio has fallen
292  * substantially below the threshold for solving the CO network.
293  * Turn it off */
294  /* >> chng 07 jan 10 rjrw: this was CO_Init(), but the comment suggests
295  * it should really be CO_zero */
296  CO_zero();
297  }
298  else
299  {
300  if( called.lgTalk && !lgFirstTry )
301  {
302  conv.lgConvPops = false;
303  fprintf( ioQQQ,
304  " PROBLEM CO_drive failed to converge1 due to negative pops, zone=%.2f, CO/C=%.2e OH/C=%.2e H2/H=%.2e\n",
305  fnzone,
306  findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]),
307  findspecies("OH")->hevmol/SDIV(dense.gas_phase[ipCARBON]),
309  ConvFail( "pops" , "CO" );
310  }
311  }
312  }
313 
314  /* this test, hit upper limit to number of iterations */
315  else if( loop == LUPMAX_CODRIV )
316  {
317  /* do not print this if we are in the first zone where molecules are turned on */
318  if( called.lgTalk && !lgFirstTry )
319  {
320  conv.lgConvPops = false;
321  fprintf( ioQQQ,
322  " PROBLEM CO_drive failed to converge2 in %li iter, zone=%.2f, CO/C=%.2e negpop?%1c reason:%s\n",
323  loop,
324  fnzone,
325  findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]),
326  TorF(lgNegPop) ,
327  chReason);
328  ConvFail( "pops" , "CO" );
329  }
330  }
331 
332  /* make sure we have not used more than all the CO */
333  ASSERT(conv.lgSearch || findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) <= 1.01 ||
334  /* this is true when loop did not converge co */
335  loop == LUPMAX_CODRIV );
336  /*fprintf(ioQQQ,"ratioo o\t%c\t%.2f\t%f\n", TorF(conv.lgSearch),
337  fnzone , co.hevmol[ipCO]/dense.gas_phase[ipCARBON] );*/
338 
339  /* these are the elements whose converge we check */
340  /* this is a self consistency check, for converged solution */
341  /* >>chng 04 dec 02, this test is turn back on - don't know why it was turned off */
342  if(0) {
343  double sum_ion , sum_mol;
344  sum_ion = dense.xIonDense[ipCARBON][0] + dense.xIonDense[ipCARBON][1];
345  sum_mol = findspecies("C")->hevmol + findspecies("C+")->hevmol;
346  if( fabs(sum_ion-sum_mol)/SDIV(sum_mol) > 1e-2 )
347  {
348  /*fprintf(ioQQQ,
349  "non conservation element %li zone %.2f sum_ion %.3e sum_mol %.3e\n",
350  5, fnzone, sum_ion , sum_mol);*/
351  conv.lgConvIoniz = false;
352  strcpy( conv.chConvIoniz, "CO con2" );
353  conv.BadConvIoniz[0] = sum_ion;
354  conv.BadConvIoniz[1] = sum_mol;
355  if( CODEBUG )
356  fprintf(ioQQQ,"CO_drive not converged\n");
357  }
358  sum_ion = dense.xIonDense[ipOXYGEN][0] + dense.xIonDense[ipOXYGEN][1];
359  sum_mol = findspecies("O")->hevmol + findspecies("O+")->hevmol;
360  if( fabs(sum_ion-sum_mol)/SDIV(sum_mol) > 1e-2 )
361  {
362  /*fprintf(ioQQQ,
363  "non conservation element %li zone %.2f sum_ion %.3e sum_mol %.3e\n",
364  7, fnzone, sum_ion , sum_mol);*/
365  conv.lgConvIoniz = false;
366  strcpy( conv.chConvIoniz, "CO con3" );
367  conv.BadConvIoniz[0] = sum_ion;
368  conv.BadConvIoniz[1] = sum_mol;
369  if( CODEBUG )
370  fprintf(ioQQQ,"CO_drive not converged\n");
371  }
372  }
373  return;
374 }
375 
377  /* job == "SETUP", set things up during first call,
378  * job == "AVER" take average of arrays */
379  const char chJob[10] ,
380  char *chReason )
381 {
382  long int i;
383  bool lgReturn;
384  /* this will be used to rescale old saved abundances in constant pressure model */
385  static realnum hden_save_prev_call;
386  double conv_fac;
387  struct chem_element_s *element;
388 
389  DEBUG_ENTRY( "lgMolecAver()" );
390  const bool DEBUG_MOLECAVER = false;
391 
392  /* this will become reason not converged */
393  strcpy( chReason , "none given" );
394 
395  /* when called with SETUP, just about to enter solver loop */
396  if( strcmp( "SETUP" , chJob ) == 0 )
397  {
398  static realnum hden_save_prev_iter;
399 
400  /* >>chng 04 apr 29, use PDR solution, scaled to density of neutral carbon, as first guess*/
401  /* CO_Init set this to -2 when code initialized, so negative
402  * number shows very first call in this model */
403  /* >>chng 05 jul 15, do not init if we have never evaluated CO in this iteration
404  * and we are below limit where it should be evaluated */
405  if( iteration!=co.iteration_co &&
407  {
408 
409  lgReturn = true;
410  return lgReturn;
411  }
412 
413  else if( co.iteration_co < 0 || lgMoleZeroed )
414  {
415 
416  /* very first attempt at ever obtaining a solution -
417  * called one time per model since co.iteration_co set negative
418  * when cdInit called */
419 
420  /* >>chng 05 jun 24, during map phase do not reset molecules to zero
421  * since we probably have a better estimate right now */
423  {
424  for( i=0; i < mole.num_comole_calc; i++ )
425  {
427  COmole[i]->co_save = COmole[i]->hevmol;
428  }
429  }
430 
431  /* we should have a neutral carbon solution at this point */
433 
434  /* first guess of the elements and charges come from ionization solvers */
435  for(i=0;i<mole.num_elements;i++) {
436  element = chem_element[i];
437  COmole[element->ipMl]->hevmol = dense.xIonDense[element->ipCl][0];
438  COmole[element->ipMlP]->hevmol = dense.xIonDense[element->ipCl][1];
439  }
440 
441  /* set iteration flag */
443  co.co_nzone = nzone;
444 
445  /* for constant pressure models when molecules are reset on second
446  * and higher iterations, total density will be different, so
447  * must take this into account when abundances are reset */
448  hden_save_prev_iter = dense.gas_phase[ipHYDROGEN];
449  hden_save_prev_call = dense.gas_phase[ipHYDROGEN];
450 
451  if( DEBUG_MOLECAVER )
452  fprintf(ioQQQ," DEBUG lgMolecAver iteration %li zone %li zeroing since very first call. H2/H=%.2e\n",
454  }
455  else if( iteration > co.iteration_co )
456  {
457  realnum ratio = dense.gas_phase[ipHYDROGEN] / hden_save_prev_iter;
458 
459  /* this is first call on new iteration, reset
460  * to first initial abundances from previous iteration */
461  for( i=0; i < mole.num_comole_calc; i++ )
462  {
463  COmole[i]->hevmol = COmole[i]->hev_reinit*ratio;
464  COmole[i]->co_save = COmole[i]->hev_reinit*ratio;
465  }
466 
468  co.co_nzone = nzone;
469 
470  if( DEBUG_MOLECAVER )
471  fprintf(ioQQQ," DEBUG lgMolecAver iteration %li zone %li resetting since first call on new iteration. H2/H=%.2e\n",
472  iteration,
473  nzone,
475  }
476  else if( iteration == co.iteration_co && nzone==co.co_nzone+1 )
477  {
478  /* this branch, second zone with solution, so save previous
479  * zone's solution to reset things in next iteration */
480  for( i=0; i < mole.num_comole_calc; i++ )
481  {
482  COmole[i]->hev_reinit = COmole[i]->hevmol;
483  }
484 
485  co.co_nzone = -2;
486  hden_save_prev_iter = dense.gas_phase[ipHYDROGEN];
487  hden_save_prev_call = dense.gas_phase[ipHYDROGEN];
488 
489  if( DEBUG_MOLECAVER )
490  fprintf(ioQQQ,"DEBUG lgMolecAver iteration %li zone %li second zone on new iteration, saving reset.\n", iteration,nzone);
491  }
492 
493  /* didn't do anything, but say converged */
494  lgReturn = true;
495  }
496 
497  else if( strcmp( "AVER" , chJob ) == 0 )
498  {
499  /* >>chng 04 jan 22, co.hevmol is current value, we want to save old
500  * value, which is in co.co_save */
501  /*realnum oldoh = findspecies("OH")->hevmol;*/
502  realnum oldoh = findspecies("OH")->co_save;
503  lgReturn = true;
504 
505  /* get new numbers - take average of old and new */
506  /* >>chng 03 sep 16, only use damper for molecular species,
507  * not ion/atoms */
508  for( i=0; i < mole.num_comole_calc; i++ )
509  {
510  /* parameter to mix old and new,
511  * damper is fraction of old solution to take */
512  realnum damper = 0.2f;
513 
514  /* >>chng 04 sep 11, correct for den chng in var den models */
515  /* the ratio of densities is unity for constant density model, but in const pres or other
516  * models, account for change in entire density scale */
517  COmole[i]->co_save *= dense.gas_phase[ipHYDROGEN] / hden_save_prev_call;
518 
519  /* if co.co_save is zero then don't take average of it and
520  * new solution */
521  if( COmole[i]->co_save< SMALLFLOAT )
522  {
523  COmole[i]->co_save = COmole[i]->hevmol;
524  }
525  else
526  {
527  /* >>chng 04 jan 24, add logic to check how different old and
528  * new solutions are. if quite different then take average
529  * in the log */
530  double ratio = (double)COmole[i]->hevmol/(double)COmole[i]->co_save;
531  ASSERT( COmole[i]->co_save>0. && COmole[i]->hevmol>=0. );
532  if( fabs(ratio-1.) < 1.5 ||
533  fabs(ratio-1.) > 0.66 )
534  {
535  /* this branch moderate changes */
536  COmole[i]->hevmol = COmole[i]->hevmol*(1.f - damper) +
537  COmole[i]->co_save*damper;
538  }
539  else
540  {
541  /* this branch - large changes take average in the log */
542  COmole[i]->hevmol = (realnum)pow(10. ,
543  (log10(COmole[i]->hevmol) + log10(COmole[i]->co_save))/2. );
544  }
545  COmole[i]->co_save = COmole[i]->hevmol;
546  }
547  }
548  /* remember the current H density */
549  hden_save_prev_call = dense.gas_phase[ipHYDROGEN];
550 
551  /* >>chng 05 jul 14, add logic to not be as fine a tolerance when
552  * only trace amounts of molecules are present */
553  if( oldoh/SDIV(dense.gas_phase[ipOXYGEN]) < 1e-10 )
554  conv_fac = 3.;
555  else
556  conv_fac = 1.;
557 
558  if(fabs(oldoh/SDIV(findspecies("OH")->hevmol)-1.) < COTOLER_MOLAV*conv_fac )
559  {
560  lgReturn = true;
561  }
562  else
563  {
564  /* say we are not converged */
565  lgReturn = false;
566  sprintf( chReason , "OH change from %.3e to %.3e",
567  oldoh ,
568  findspecies("OH")->hevmol);
569  }
570  }
571  else
572  TotalInsanity();
573 
574  /* also not converged if co exceeds c */
575  if( findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) > 1.+FLT_EPSILON )
576  {
577  lgReturn = false;
578  sprintf( chReason , "CO/C >1, value is %e",
579  findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) );
580  }
581 
582  if( (trace.lgTrace && trace.lgTr_CO_Mole) || DEBUG_MOLECAVER )
583  {
584  fprintf( ioQQQ, " DEBUG lgMolecAver converged? %c" ,TorF(lgReturn) );
585  fprintf( ioQQQ, " CO/C=%9.1e", findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) );
586  fprintf( ioQQQ, " OH/O=%9.1e", findspecies("OH")->hevmol/SDIV(dense.gas_phase[ipOXYGEN]) );
587  fprintf( ioQQQ, "\n" );
588  }
589 
590  return lgReturn;
591 }

Generated for cloudy by doxygen 1.8.1.1