cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
conv_base.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 /*ConvBase main routine to drive ionization solution for all species, find total opacity
4  * called by ConvIoniz */
5 /*lgConverg check whether ionization of element nelem has converged */
6 #include "cddefines.h"
7 #include "dynamics.h"
8 #include "trace.h"
9 #include "elementnames.h"
10 #include "punch.h"
11 #include "phycon.h"
12 #include "secondaries.h"
13 #include "stopcalc.h"
14 #include "grainvar.h"
15 #include "highen.h"
16 #include "dense.h"
17 #include "hmi.h"
18 #include "rfield.h"
19 #include "pressure.h"
20 #include "taulines.h"
21 #include "rt.h"
22 #include "grains.h"
23 #include "atmdat.h"
24 #include "ionbal.h"
25 #include "opacity.h"
26 #include "cooling.h"
27 #include "thermal.h"
28 #include "mole.h"
29 #include "iso.h"
30 #include "conv.h"
31 
32 /*lgIonizConverg check whether ionization of element nelem has converged, called by ionize,
33  * returns true if element is converged, false if not */
35  /* atomic number on the C scale, 0 for H, 25 for Fe */
36  long int nelem ,
37  /* this is allowed error as a fractional change. Most are 0.15 */
38  double delta ,
39  /* option to print abundances */
40  bool lgPrint )
41 {
42  bool lgConverg_v;
43  long int ion;
44  double Abund,
45  bigchange ,
46  change ,
47  one;
48  double abundold=0. , abundnew=0.;
49 
50  /* this is fractions [ion stage][nelem], ion stage = 0 for atom, nelem=0 for H*/
51  static realnum OldFracs[LIMELM+1][LIMELM+1];
52 
53  if( !dense.lgElmtOn[nelem] )
54  return true;
55 
56  DEBUG_ENTRY( "lgIonizConverg()" );
57 
58  /* arguments are atomic number, ionization stage
59  * OldFracs[nelem][0] is abundance of nelem (cm^-3) */
60 
61  /* this function returns true if ionization of element
62  * with atomic number nelem has not changed by more than delta,*/
63 
64  /* check whether abundances exist yet, only do this for after first zone */
65  /*if( nzone > 0 )*/
66  /* >>chng 03 sep 02, check on changed ionization after first call through,
67  * to insure converged constant eden / temperature models */
68  if( conv.nPres2Ioniz )
69  {
70  /* >>chng 04 aug 31, this had been static, caused last hits on large changes
71  * in ionization */
72  bigchange = 0.;
73  change = 0.;
74  Abund = dense.gas_phase[nelem];
75 
76  /* loop over all ionization stages, loop over all ions, not just active ones,
77  * since this also sets old values, and so will zero out non-existant ones */
78  for( ion=0; ion <= (nelem+1); ++ion )
79  {
80  /*lint -e727 OlsdFracs not initialized */
81  if( OldFracs[nelem][ion]/Abund > 1e-4 &&
82  dense.xIonDense[nelem][ion]/Abund > 1e-4 )
83  /*lint +e727 OlsdFracs not initialized */
84  {
85  /* check change in old to new value */
86  one = fabs(dense.xIonDense[nelem][ion]-OldFracs[nelem][ion])/
87  OldFracs[nelem][ion];
88  change = MAX2(change, one );
89  /* remember abundances for largest change */
90  if( change>bigchange )
91  {
92  bigchange = change;
93  abundold = OldFracs[nelem][ion]/Abund;
94  abundnew = dense.xIonDense[nelem][ion]/Abund;
95  }
96  }
97  /* now set new value */
98  OldFracs[nelem][ion] = dense.xIonDense[nelem][ion];
99  }
100 
101  if( change < delta )
102  {
103  lgConverg_v = true;
104  }
105  else
106  {
107  lgConverg_v = false;
108  conv.BadConvIoniz[0] = abundold;
109  conv.BadConvIoniz[1] = abundnew;
110  ASSERT( abundold>0. && abundnew>0. );
111  }
112  }
113  else
114  {
115  for( ion=0; ion <= (nelem+1); ++ion )
116  {
117  OldFracs[nelem][ion] = dense.xIonDense[nelem][ion];
118  }
119 
120  lgConverg_v = true;
121  }
122 
123  /* option to print abundances */
124  if( lgPrint )
125  {
126  fprintf(ioQQQ," element %li converged? %c ",nelem, TorF(lgConverg_v));
127  for( ion=0; ion<(nelem+1); ++ion )
128  {
129  fprintf(ioQQQ,"\t%.2e", dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]);
130  }
131  fprintf(ioQQQ,"\n");
132  }
133  return( lgConverg_v );
134 }
135 
136 /*ConvBase main routine to drive ionization solution for all species, find total opacity
137  * called by ConvIoniz
138  * return 0 if ok, 1 if abort */
139 int ConvBase(
140  /* this is zero the first time ConvBase is called by convIoniz,
141  * counts number of call thereafter */
142  long loopi )
143 {
144  double O_HIonRate_New , O_HIonRate_Old;
145  double HeatOld,
146  EdenTrue_old,
147  EdenFromMolecOld,
148  EdenFromGrainsOld,
149  HeatingOld ,
150  CoolingOld;
151  realnum hmole_save[N_H_MOLEC];
152  static double SecondOld;
153  static long int nzoneOTS=-1;
154 # define LOOP_ION_LIMIT 10
155  long int nelem,
156  ion,
157  loop_ion,
158  i,
159  ipISO;
160  static double SumOTS=0. , OldSumOTS[2]={0.,0.};
161  double save_iso_grnd[NISO][LIMELM];
162 
163  DEBUG_ENTRY( "ConvBase()" );
164 
165 #if 0
166  if( loopi == 0 )
167  {
168  OldSumOTS[0] = 0.;
169  OldSumOTS[1] = 0.;
170  SumOTS = 0.;
171  nzoneOTS = -1;
172  }
173 #endif
174 
175  /* this is set to phycon.te in tfidle, is used to insure that all temp
176  * vars are properly updated when conv_ionizeopacitydo is called
177  * NB must be same type as phycon.te */
179 
180  /* this allows zone number to be printed with slight increment as zone converged
181  * conv.nPres2Ioniz is incremented at the bottom of this routine */
182  fnzone = (double)nzone + (double)conv.nPres2Ioniz/100.;
183 
184  /* reevaluate pressure */
185  /* this sets values of pressure.PresTotlCurr, also calls tfidle,
186  * and sets the total energy content of gas, which may be important for acvection */
187  PresTotCurrent();
188 
189  /* >>chng 04 sep 15, move EdenTrue_old up here, and will redo at bottom
190  * to find change
191  * find and save current true electron density - if this changes by more than the
192  * tolerance then ionization solution is not converged */
193  /* >>chng 04 jul 27, move eden_sum here from after this loop, so that change in EdenTrue
194  * can be monitored */
195  /* update EdenTrue, eden itself is actually changed in ConvEdenIoniz */
196  /* must not call eden_sum on very first time since for classic PDR total
197  * ionization may still be zero on first call */
198  if( conv.nTotalIoniz )
199  {
200  if( eden_sum() )
201  {
202  /* non-zero return indicates abort condition */
203  ++conv.nTotalIoniz;
204  return 1;
205  }
206  }
207 
208  /* the following two were evaluated in eden_sum
209  * will confirm that these are converged */
210  EdenTrue_old = dense.EdenTrue;
211  EdenFromMolecOld = co.comole_eden;
212  EdenFromGrainsOld = gv.TotalEden;
213  HeatingOld = thermal.htot;
214  CoolingOld = thermal.ctot;
215 
216  /* remember current ground state population - will check if converged */
217  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
218  {
219  for( nelem=ipISO; nelem<LIMELM;++nelem )
220  {
221  if( dense.lgElmtOn[nelem] )
222  {
223  /* save the ground state population */
224  save_iso_grnd[ipISO][nelem] = StatesElem[ipISO][nelem][0].Pop;
225  }
226  }
227  }
228 
229  for( i=0; i < mole.num_comole_calc; i++ )
230  {
231  COmole[i]->comole_save = COmole[i]->hevmol;
232  }
233 
234  for( i=0; i < N_H_MOLEC; i++ )
235  {
236  hmole_save[i] = hmi.Hmolec[i];
237  }
238 
239  if( loopi==0 )
240  {
241  /* these will be used to look for oscillating ots rates */
242  OldSumOTS[0] = 0.;
243  OldSumOTS[1] = 0.;
244  }
245 
246  if( trace.lgTrace )
247  {
248  fprintf( ioQQQ,
249  " ConvBase called. %.2f Te:%.3e HI:%.3e HII:%.3e H2:%.3e Ne:%.3e htot:%.3e CSUP:%.2e Conv?%c\n",
250  fnzone,
251  phycon.te,
254  hmi.H2_total,
255  dense.eden,
256  thermal.htot,
258  TorF(conv.lgConvIoniz) );
259  }
260  /* want this flag to be true when we exit, various problems will set falst */
261  conv.lgConvIoniz = true;
262 
263  /* this routine is in heatsum.c, and zeros out the heating array */
264  HeatZero();
265 
266  /* if this is very first call, say not converged, since no meaningful way to
267  * check on changes in quantities. this counter is false only on very first
268  * call, when equal to zero */
269  if( !conv.nTotalIoniz )
270  {
271  conv.lgConvIoniz = false;
272  strcpy( conv.chConvIoniz, "first call" );
273  }
274 
275  /* this will be flag to check whether any ionization stages
276  * were trimmed down */
277  conv.lgIonStageTrimed = false;
278 
279  /* must redo photoionization rates for all species on second and later tries */
280  /* always reevaluate the rates when . . . */
281  /* this flag says not to redo opac and photo rates, and following test will set
282  * true if one of several tests not done*/
283  opac.lgRedoStatic = false;
284  if(
285  /* opac.lgOpacStatic (usually true), is set false with no static opacity command */
286  !opac.lgOpacStatic ||
287  /* we are in search mode */
288  conv.lgSearch ||
289  /* this is the first call to this zone */
290  conv.nPres2Ioniz == 0 )
291  {
292  /* we need to redo ALL photoionization rates */
293  opac.lgRedoStatic = true;
294  }
295 
296  /* calculate radiative and dielectronic recombination rate coefficients */
298 
299  /* this adjusts the lowest and highest stages of ionization we will consider,
300  * only safe to call when lgRedoStatic is true since this could lower the
301  * lowest stage of ionization, which needs all its photo rates */
302 
303  /* conv.nTotalIoniz is only 0 (false) on the very first call to here,
304  * when the ionization distribution is not yet done */
305  if( conv.nTotalIoniz )
306  {
307  bool lgIonizTrimCalled = false;
308  static long int nZoneCalled = 0;
309 
310  fixit(); // nZoneCalled should be reinitialized for each grid point?
311 
312  /* ionization trimming only used for He and heavier, not H */
313  /* only do this one time per zone since up and down cycle can occur */
314  /* >>chng 05 jan 15, increasing temperature above default first conditions, also
315  * no trim during search - this fixed major logical error when sim is
316  * totally mechanically heated to coronal temperatures -
317  * problem discovered by Ronnie Hoogerwerf */
318  /* do not keep changing the trim after the first call within
319  * this zone once we are deep into layer - doing so introduced very
320  * small level of noise as some stages
321  * went up and down - O+2 - O+3 was good example, when small H+3 after He+ i-front
322  * limit to one increase per element per zone */
323  if( conv.nTotalIoniz>2 &&
324  /* only call one time per zone except during search phase,
325  * when only call after 20 times only if temperature has changed */
326  ( conv.lgSearch || nZoneCalled!=nzone) )
327  {
328  lgIonizTrimCalled = true;
329  for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
330  {
331  if( dense.lgElmtOn[nelem] )
332  {
333  /* ion_trim will set conv.lgIonStageTrimed true is any ion has its
334  * lowest stage of ionization dropped or trimmed */
335  ion_trim(nelem);
336  }
337  }
338  nZoneCalled = nzone;
339  }
340 
341  /* following block only set of asserts */
342 # if !defined(NDEBUG)
343  /* check that proper abundances are either positive or zero */
344  for( nelem=ipHELIUM; nelem<LIMELM; ++nelem)
345  {
346  if( dense.lgElmtOn[nelem] )
347  {
348  for( ion=0; ion<dense.IonLow[nelem]; ++ion )
349  {
350  ASSERT( dense.xIonDense[nelem][ion] == 0. );
351  }
352  /*if( nelem==5 ) fprintf(ioQQQ,"carbbb\t%li\n", dense.IonHigh[nelem]);*/
353  for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
354  {
355  /* >>chng 02 feb 06, had been > o., chng to > SMALLFLOAT to
356  * trip over VERY small floats that failed on alphas, but not 386
357  *
358  * in case where lower ionization stage was just lowered or
359  * trimmed down the abundance
360  * was set to SMALLFLOAT so test must be < SMALLFLOAT */
361  /* >>chng 02 feb 19, add check for search phase. During this search
362  * models with extreme ionization (all neutral or all ionized) can
363  * have extreme but non-zero abundances far from the ionization peak for
364  * element with lots of electrons. These will go away once the model
365  * becomes stable */
366  /* >>chng 03 dec 01, add check on whether ion trim was called
367  * conserve.in threw assert when iontrim not called and abund grew small */
368  ASSERT( conv.lgSearch || !lgIonizTrimCalled ||
369  /* this can happen if all C is in the form of CO */
370  (ion==0 && dense.IonHigh[nelem]==0 ) ||
371  dense.xIonDense[nelem][ion] >= SMALLFLOAT ||
372  dense.xIonDense[nelem][ion]/dense.gas_phase[nelem] >= SMALLFLOAT );
373  }
374  for( ion=dense.IonHigh[nelem]+1; ion<nelem+1; ++ion )
375  {
376  ASSERT( ion >= 0 );
377  ASSERT( dense.xIonDense[nelem][ion] == 0. );
378  }
379  }
380  }
381 # endif
382  }
383 
384  /* now check if anything trimmed down */
385  if( conv.lgIonStageTrimed )
386  {
387  /* something was trimmed down, so say that ionization not yet stable */
388  /* say that ionization has not converged, secondaries changed by too much */
389  conv.lgConvIoniz = false;
390  strcpy( conv.chConvIoniz, "IonTrimmed" );
391  }
392 
393  /* reevaluate advective terms if turned on */
394  if( dynamics.lgAdvection )
395  DynaIonize();
396 
397  /* >>chng 04 feb 15, add loop over ionization until converged. Non-convergence
398  * in this case means ionization ladder pop changed, probably because of way
399  * that Auger ejection is treated - loop exits when conditions tested at end are met */
400  loop_ion = 0;
401  do
402  {
403 
404  if( trace.nTrConvg>=4 )
405  {
406  /* trace ionization */
407  fprintf( ioQQQ,
408  " ConvBase4 ionization driver loop_ion %li converged? %c reason not converged %s\n" ,
409  loop_ion ,
410  TorF( conv.lgConvIoniz) ,
411  conv.chConvIoniz );
412  }
413 
414  /* >>chng 04 sep 11, moved to before grains and other routines - must not clear this
415  * flag before their calculations are done */
416  /* set the convergence flag to true,
417  * if anything changes in ConvBase, it will be set false */
418  if( loop_ion )
419  conv.lgConvIoniz = true;
420 
421  /* >>chng 01 apr 29, move charge transfer evaluation here, had been just before
422  * HeLike, but needs to be here so that same rate coefficient used for H ion and other recomb */
423  /* fill in master charge transfer array, and zero out some arrays that track effects,
424  * O_HIonRateis rate oxygen ionizes hydrogen - will need to converge this */
425  ChargTranEval( &O_HIonRate_Old );
426 
427  CO_update_species_cache(); /* Update densities of species controlled outside the chemical network */
428 
429  hmole_reactions();
430 
431  co.nitro_dissoc_rate = 0;
432 
433  /* Update the rate constants -- generic version */
434  CO_update_rks();
435 
436  /* find grain temperature, charge, and gas heating rate */
437  /* >>chng 04 apr 15, moved here from after call to HeLike(), PvH */
438  GrainDrive();
439 
440  /* do all hydrogenic species, and fully do hydrogen itself, including molecules */
441  /*fprintf(ioQQQ,"DEBUG h2\t%.2f\t%.3e\t%.3e", fnzone,hmi.H2_total,hmi.Hmolec[ipMH2g]);*/
442  iso_solve( ipH_LIKE );
443  /*fprintf(ioQQQ,"\t%.3e\n", hmi.Hmolec[ipMH2g]);*/
444 
445  /* evaluate Compton heating, bound E Compton, cosmic rays */
446  highen();
447 
448  /* find corrections for three-body rec - collisional ionization */
449  atmdat_3body();
450 
451  /* deduce dielectronic suppression factors */
453 
454  /* >>chng 02 mar 08 move helike to before helium */
455  /* do all he-like species */
456  iso_solve( ipHE_LIKE );
457 
458  /*>>chng 04 may 09, add option to abort here, inspired by H2 pop failures
459  * which cascaded downstream if we did not abort */
460  /* return if one of above solvers has declared abort condition */
461  if( lgAbort )
462  {
463  ++conv.nTotalIoniz;
464  return 1;
465  }
466 
467  /* find grain temperature, charge, and gas heating rate */
468  /*GrainDrive();*/
469 
470  /* only reevaluate this on first pass through on this zone */
471  if( opac.lgRedoStatic )
472  {
473  /* inner shell ionization */
474  for( nelem=0; nelem< LIMELM; ++nelem )
475  {
476  for( ion=0; ion<nelem+1; ++ion )
477  {
478  ionbal.UTA_ionize_rate[nelem][ion] = 0.;
479  ionbal.UTA_heat_rate[nelem][ion] = 0.;
480  }
481  }
482  /* inner shell ionization by UTA lines */
483  /* this flag turned off with no UTA command */
485  {
486  for( i=0; i<nUTA; ++i )
487  {
488  if( UTALines[i].Emis->Aul > 0. )
489  {
490  /* cs is actually the negative of the branching ratio for autoionization,
491  * rateone is inverse lifetime of level against autoionization */
492  double rateone = UTALines[i].Emis->pump * -UTALines[i].Coll.cs;
493  ionbal.UTA_ionize_rate[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1] += rateone;
494  /* heating rate in erg atom-1 s-1 */
495  ionbal.UTA_heat_rate[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1] += rateone*UTALines[i].Coll.heat;
496  {
497  /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */
498  /*@-redef@*/
499  enum {DEBUG_LOC=false};
500  /*@+redef@*/
501  if( DEBUG_LOC /*&& UTALines[i].nelem==ipIRON+1 && (UTALines[i].IonStg==15||UTALines[i].IonStg==14)*/ )
502  {
503  fprintf(ioQQQ,"DEBUG UTA %3i %3i %.3f %.2e %.2e %.2e\n",
504  UTALines[i].Hi->nelem , UTALines[i].Hi->IonStg , UTALines[i].WLAng ,
505  rateone, UTALines[i].Coll.heat,
506  UTALines[i].Coll.heat*dense.xIonDense[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1] );
507  }
508  }
509  }
510  }
511  }
512  }
513 
514  /* helium ionization balance */
515  IonHelium();
516 
517  /* do the ionization balance */
518  /* evaluate current opacities, OpacityAddTotal is called only here during actual calculation */
519  /* option to only evaluate opacity one time per zone,
520  * rfield.lgOpacityReevaluate normally true, set false with no opacity reevaluate command */
521  /* >>chng 04 jul 16, move these out of block below, so always done
522  * these all produce ions that are used in the CO network */
523  IonCarbo();
524  IonOxyge();
525  IonNitro();
526  IonSilic();
527  IonSulph();
528  IonChlor();
529  /* do carbon monoxide after oxygen */
530  CO_drive();
533  {
534  IonLithi();
535  IonBeryl();
536  IonBoron();
537  /*IonCarbo();
538  IonOxyge();*/
539  /* do carbon monoxide after oxygen */
540  /*fprintf(ioQQQ,"DEBUG co\t%.2f\t%.3e", fnzone,findspecies("CO")->hevmol);*/
541  /*CO_drive();*/
542  /*fprintf(ioQQQ,"\t%.3e\n", findspecies("CO")->hevmol);*/
543  IonFluor();
544  IonNeon();
545  IonSodiu();
546  IonMagne();
547  IonAlumi();
548  IonPhosi();
549  IonArgon();
550  IonPotas();
551  IonCalci();
552  IonScand();
553  IonTitan();
554  IonVanad();
555  IonChrom();
556  IonManga();
557  IonIron();
558  IonCobal();
559  IonNicke();
560  IonCoppe();
561  IonZinc();
562  }
563 
564  /* all elements have now had ionization reevaluated - in some cases we may have upset
565  * the ionization that was forced with an "element ionization" command - here we will
566  * re-impose that set ionization */
567  /* >>chng 04 oct 13, add this logic */
568  for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
569  {
570  if( dense.lgSetIoniz[nelem] )
571  {
572  dense.IonLow[nelem] = 0;
573  dense.IonHigh[nelem] = nelem + 1;
574  while( dense.SetIoniz[nelem][dense.IonLow[nelem]] == 0. )
575  ++dense.IonLow[nelem];
576  while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] == 0. )
577  --dense.IonHigh[nelem];
578  }
579  }
580 
581  /* redo Oxygen ct ion rate of H to see if it changed */
582  ChargTranEval( &O_HIonRate_New );
583 
584  /* lgIonizConverg is a function to check whether ionization has converged
585  * check whether ionization changed by more than relative error
586  * given by second number */
587  /* >>chng 04 feb 14, loop over all elements rather than just a few */
588  for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
589  {
590  if( !lgIonizConverg(nelem, 0.05 ,false ) )
591  {
592  conv.lgConvIoniz = false;
593  sprintf( conv.chConvIoniz , "%2s ion" , elementnames.chElementSym[nelem] );
594  }
595  }
596 
597  if( fabs(EdenTrue_old - dense.EdenTrue)/SDIV(dense.EdenTrue) > conv.EdenErrorAllowed/2. )
598  {
599  conv.lgConvIoniz = false;
600  sprintf( conv.chConvIoniz , "EdTr cng" );
601  conv.BadConvIoniz[0] = EdenTrue_old;
603  }
604  ++loop_ion;
605  }
606  /* loop is not converged, less than loop limit, and we are reevaluating */
607  while( !conv.lgConvIoniz && loop_ion < LOOP_ION_LIMIT && rfield.lgIonizReevaluate);
608 
609  /* >>chng 05 oct 29, move CT heating here from heat_sum since this sometimes contributes
610  * cooling rather than heat and so needs to be sorted out before either heating or cooling
611  * are derived first find net heating - */
613  /* net is cooling if negative */
616 
617  /* get total cooling, thermal.ctot = does not occur since passes as pointer. This can add heat.
618  * it calls coolSum at end to sum up the total cooling */
620 
621  /* get total heating rate - first save old quantities to check how much it changes */
622  HeatOld = thermal.htot;
623 
624  /* HeatSum will update ElecFrac,
625  * secondary electron ionization and excitation efficiencies,
626  * and sum up current secondary rates - remember old values before we enter */
627  SecondOld = secondaries.csupra[ipHYDROGEN][0];
628 
629  /* update the total heating - it was all set to zero in HeatZero at top of this routine,
630  * occurs before secondaries bit below, since this updates electron fracs */
631  HeatSum();
632 
633  /* test whether we can just set the rate to the new one, or whether we should
634  * take average and do it again. secondaries.sec2total was set in hydrogenic, and
635  * is ratio of secondary to total hydrogen destruction rates */
636  /* >>chng 02 nov 20, add test on size of csupra - primal had very close to underflow */
638  fabs( 1. - SecondOld/SDIV(secondaries.csupra[ipHYDROGEN][0]) ) > 0.1 &&
639  SecondOld > 0. && secondaries.csupra[ipHYDROGEN][0] > 0.)
640  {
641  /* say that ionization has not converged, secondaries changed by too much */
642  conv.lgConvIoniz = false;
643  strcpy( conv.chConvIoniz, "SecIonRate" );
644  conv.BadConvIoniz[0] = SecondOld;
646  }
647 
648 # if 0
649  static realnum hminus_old=0.;
650  /* >>chng 04 apr 15, add this convergence test */
651  if( conv.nTotalIoniz )
652  {
653  if( fabs( hminus_old-hmi.Hmolec[ipMHm])/SDIV(hmi.Hmolec[ipMHm]) >
655  conv.lgConvIoniz = false;
656  strcpy( conv.chConvIoniz, "Big H- chn" );
657  conv.BadConvIoniz[0] = hminus_old;
659  }
660  hminus_old = hmi.Hmolec[ipMHm];
661 # endif
662 
663  if( HeatOld > 0. && !thermal.lgTemperatureConstant )
664  {
665  /* check if heating has converged - tolerance is final match */
666  if( fabs(1.-thermal.htot/HeatOld) > conv.HeatCoolRelErrorAllowed*.5 )
667  {
668  conv.lgConvIoniz = false;
669  strcpy( conv.chConvIoniz, "Big d Heat" );
670  conv.BadConvIoniz[0] = HeatOld;
672  }
673  }
674 
675  /* abort flag may have already been set - if so bail */
676  if( lgAbort )
677  {
678 
679  return 1;
680  }
681 
682  /* evaluate current opacities, OpacityAddTotal is called only here during actual calculation */
683  /* option to only evaluate opacity one time per zone,
684  * rfield.lgOpacityReevaluate normally true, set false with no opacity reevaluate command */
686  OpacityAddTotal();
687 
688  /* >>chng 02 jun 11, call even first time that this routine is called -
689  * this seemed to help convergence */
690 
691  /* do OTS rates for all lines and all continua since
692  * we now have ionization balance of all species. Note that this is not
693  * entirely self-consistent, since destruction probabilities here are not the same as
694  * the ones used in the model atoms. Problems?? if near convergence
695  * then should be nearly identical */
697  conv.lgIonStageTrimed || conv.lgSearch || nzone!=nzoneOTS )
698  {
699  RT_OTS();
700  nzoneOTS = nzone;
701 
702  /* remember old ots rates */
703  OldSumOTS[0] = OldSumOTS[1];
704  OldSumOTS[1] = SumOTS;
705  /*fprintf(ioQQQ," calling RT_OTS zone %.2f SumOTS is %.2e\n",fnzone,SumOTS);*/
706 
707  /* now update several components of the continuum, this only happens after
708  * we have gone through entire solution for this zone at least one time.
709  * there can be wild ots oscillation on first call */
710  /* the rel change of 0.2 was chosen by running hizqso - smaller increased
711  * itrzn but larger did not further decrease it. */
712  RT_OTS_Update(&SumOTS , 0.20 );
713  /*fprintf(ioQQQ,"RT_OTS_Updateee\t%.3f\t%.2e\t%.2e\n", fnzone,SumOTS , OldSumOTS[1] );*/
714  }
715  else
716  SumOTS = 0.;
717 
718  /* now check whether the ots rates changed */
719  if( SumOTS> 0. )
720  {
721  /* the ots rate must be converged to the error in the electron density */
722  /* how fine should this be converged?? originally had above, 10%, but take
723  * smaller ratio?? */
724  if( fabs(1.-OldSumOTS[1]/SumOTS) > conv.EdenErrorAllowed )
725  {
726  /* this branch, ionization not converged due to large change in ots rates.
727  * check whether ots rates are oscillating, if loopi > 1 so we have enough info*/
728  if( loopi > 1 )
729  {
730  /* here we have three values, are they changing sign? */
731  if( (OldSumOTS[0]-OldSumOTS[1]) * ( OldSumOTS[1] - SumOTS ) < 0. )
732  {
733  /* ots rates are oscillating, if so then use small fraction of new destruction rates
734  * when updating Lyman line destruction probabilities in HydroPesc*/
735  conv.lgOscilOTS = true;
736  }
737  }
738 
739  conv.lgConvIoniz = false;
740  strcpy( conv.chConvIoniz, "OTSRatChng" );
741  conv.BadConvIoniz[0] = OldSumOTS[1];
742  conv.BadConvIoniz[1] = SumOTS;
743  }
744  /* produce info on the ots fields if either "trace ots" or
745  * "trace convergence xxx ots " was entered */
746  if( ( trace.lgTrace && trace.lgOTSBug ) ||
747  ( trace.nTrConvg && trace.lgOTSBug ) )
748  {
749  RT_OTS_PrtRate(SumOTS*0.05 , 'b' );
750  }
751  /*fprintf(ioQQQ,"DEBUG opac\t%.2f\t%.3e\t%.3e\n",fnzone,
752  dense.xIonDense[ipNICKEL][0] ,
753  dense.xIonDense[ipZINC][0] );*/
754  {
755  /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */
756  /*@-redef@*/
757  enum {DEBUG_LOC=false};
758  /*@+redef@*/
759  if( DEBUG_LOC && (nzone>110) )
760  {
761 # if 0
762 # include "lines_service.h"
764 # endif
765  /* last par 'l' for line, 'c' for continua, 'b' for both,
766  * the numbers printed are:
767  * cell i, [i], so 1 less than ipoint
768  * anu[i],
769  * otslin[i],
770  * opacity_abs[i],
771  * otslin[i]*opacity_abs[i],
772  * rfield.chLineLabel[i] ,
773  * rfield.line_count[i] */
774  }
775  }
776  }
777  {
778  /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */
779  /*@-redef@*/
780  enum {DEBUG_LOC=false};
781  /*@+redef@*/
782  if( DEBUG_LOC && (nzone>200) )
783  {
784  fprintf(ioQQQ,"debug otsss\t%li\t%.3e\t%.3e\t%.3e\n",
785  nzone,
786  Transitions[0][1][15][3].Emis->ots,
787  TauLines[26].Emis->ots,
788  opac.opacity_abs[2069]);
789  }
790  }
791 
792  /* option to print OTS continuum with TRACE OTS */
793  if( trace.lgTrace && trace.lgOTSBug )
794  {
795  /* find ots rates here, so we only print fraction of it,
796  * SumOTS is both line and continuum contributing to ots, and is multiplied by opacity */
797  /* number is weakest rate to print */
798  RT_OTS_PrtRate(SumOTS*0.05 , 'b' );
799  }
800 
801  /* now reevaluate only destruction probabilities if call is not first*/
802  /* >>chng 02 jun 13, had always been false, now reevaluate escape probabilities
803  * and pumping rates on first call for each zone */
804  if( conv.nPres2Ioniz && !conv.lgSearch )
805  {
806  RT_line_all(false , false );
807  }
808  else
809  {
810  RT_line_all(true , false );
811  }
812 
813  /* >>chng 01 mar 16, evaluate pressure here since changing and other values needed */
814  /* reevaluate pressure */
815  /* this sets values of pressure.PresTotlCurr, also calls tfidle */
816  PresTotCurrent();
817 
818  if( trace.lgTrace )
819  {
820  fprintf( ioQQQ,
821  " ConvBase return. %.2f Te:%.3e HI:%.3e HII:%.3e H2:%.3e Ne:%.3e htot:%.3e CSUP:%.2e Conv?%c\n",
822  fnzone,
823  phycon.te,
826  hmi.H2_total,
827  dense.eden,
828  thermal.htot,
830  TorF(conv.lgConvIoniz) );
831  }
832 
833  /* this checks that all molecules and ions add up to gas phase abundance
834  * but only if we have a converged solution */
835  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
836  {
837  /* check that species add up if element is turned on, and
838  * ionization is not set */
839  if( dense.lgElmtOn[nelem] && !dense.lgSetIoniz[nelem])
840  {
841  /* this sum of over the molecules did not include the atom and first
842  * ion that is calculated in the molecular solver */
843  double sum = dense.xMolecules[nelem];
844  for( ion=0; ion<nelem+2; ++ion )
845  {
846  sum += dense.xIonDense[nelem][ion];
847  }
848  ASSERT( sum > SMALLFLOAT );
849  /* check that sum agrees with total gas phase abundance */
853  /*if( fabs(sum-dense.gas_phase[nelem])/sum > 1e-2 )
854  {
855  fprintf(ioQQQ,"DEBUG non conservation element %li sum %.3e gas phase %.3e rel error %.3f\n",
856  nelem, sum , dense.gas_phase[nelem],(sum-dense.gas_phase[nelem])/sum);
857  }*/
858  /* >>chng 04 jul 23, needed to increase to 2e-2 from 1e-2 to get conserve.in to work,
859  * not clear what caused increase in error */
860  if( fabs(sum-dense.gas_phase[nelem])/SDIV(sum) > 2e-2 )
861  {
862  /*fprintf(ioQQQ,"DEBUG non-conv nelem\t%li\tsum\t%e\ttot gas\t%e\trel err\t%.3e\n",
863  nelem,
864  sum,
865  dense.gas_phase[nelem],
866  (sum-dense.gas_phase[nelem])/sum);*/
867  conv.lgConvIoniz = false;
868  strcpy( conv.chConvIoniz, "CO con4" );
869  sprintf( conv.chConvIoniz, "COnelem%2li",nelem );
870  conv.BadConvIoniz[0] = sum;
871  conv.BadConvIoniz[1] = dense.gas_phase[nelem];
872  }
873  }
874  }
875 
876  /* update some counters that keep track of how many times this routine
877  * has been called */
878 
879  /* this counts number of times we are called by ConvPresTempEdenIoniz,
880  * number of calls in this zone so first call is zero
881  * reset to zero each time ConvPresTempEdenIoniz is called */
882  ++conv.nPres2Ioniz;
883 
884  /* this is abort option set with SET PRESIONIZ command,
885  * test on nzone since init can take many iterations
886  * this is seldom used except in special cases */
887  if( conv.nPres2Ioniz > conv.limPres2Ioniz && nzone > 0)
888  {
889  fprintf(ioQQQ,"PROBLEM ConvBase sets lgAbort since nPres2Ioniz exceeds limPres2Ioniz.\n");
890  fprintf(ioQQQ,"Their values are %li and %li \n",conv.nPres2Ioniz , conv.limPres2Ioniz);
891  lgAbort = true;
892 
893  return 1;
894  }
895 
896  /* various checks on the convergence of the current solution */
897  /* >>chng 04 sep 15, add this check if the ionization so converged */
898  if( eden_sum() )
899  {
900  /* non-zero return indicates abort condition */
901  return 1;
902  }
903  /* is electron density converged? */
905  if( fabs(EdenTrue_old - dense.EdenTrue)/dense.EdenTrue > conv.EdenErrorAllowed/2. )
906  {
907  conv.lgConvIoniz = false;
908  sprintf( conv.chConvIoniz, "eden chng" );
909  conv.BadConvIoniz[0] = EdenTrue_old;
911  }
912 
913  /* check on molecular electron den */
914  if( fabs(EdenFromMolecOld - co.comole_eden)/dense.EdenTrue > conv.EdenErrorAllowed/2. )
915  {
916  conv.lgConvIoniz = false;
917  sprintf( conv.chConvIoniz, "edn chnCO" );
918  conv.BadConvIoniz[0] = EdenFromMolecOld;
920  }
921 
922  /* check on molecular electron den */
923  if( fabs(EdenFromGrainsOld - gv.TotalEden)/dense.EdenTrue > conv.EdenErrorAllowed/2. )
924  {
925  conv.lgConvIoniz = false;
926  sprintf( conv.chConvIoniz, "edn grn e" );
927  conv.BadConvIoniz[0] = EdenFromGrainsOld;
929  }
930 
931  /* check on heating and cooling if vary temp model */
933  {
934  if( fabs(HeatingOld - thermal.htot)/thermal.htot > conv.HeatCoolRelErrorAllowed/2. )
935  {
936  conv.lgConvIoniz = false;
937  sprintf( conv.chConvIoniz, "heat chg" );
938  conv.BadConvIoniz[0] = HeatingOld;
940  }
941 
942  if( fabs(CoolingOld - thermal.ctot)/thermal.ctot > conv.HeatCoolRelErrorAllowed/2. )
943  {
944  conv.lgConvIoniz = false;
945  sprintf( conv.chConvIoniz, "cool chg" );
946  conv.BadConvIoniz[0] = CoolingOld;
948  }
949  }
950 
951  /* check on sum of grain and molecular electron den - often two large numbers that cancel */
952  if( fabs( (EdenFromMolecOld-EdenFromGrainsOld) - (co.comole_eden-gv.TotalEden) )/dense.eden >
954  {
955  conv.lgConvIoniz = false;
956  sprintf( conv.chConvIoniz, "edn grn e" );
957  conv.BadConvIoniz[0] = (EdenFromMolecOld-EdenFromGrainsOld);
959  }
960 
961  /* check whether molecular abundances are stable - these were evaluated in CO_drive */
962  for( i=0; i < mole.num_comole_calc; ++i )
963  {
965  {
966  if( fabs(COmole[i]->hevmol-COmole[i]->comole_save)/dense.gas_phase[COmole[i]->nelem_hevmol]-1. >
968  {
969  conv.lgConvIoniz = false;
970  sprintf( conv.chConvIoniz, "ch %s",COmole[i]->label );
973  }
974  }
975  }
976 
977  /* check whether H molecular abundances are stable */
978  for( i=0; i < N_H_MOLEC; ++i )
979  {
980  if( fabs(hmi.Hmolec[i]-hmole_save[i])/dense.gas_phase[ipHYDROGEN]-1. >
982  {
983  conv.lgConvIoniz = false;
984  sprintf( conv.chConvIoniz, "ch %s",hmi.chLab[i] );
985  conv.BadConvIoniz[0] = hmole_save[i]/dense.gas_phase[ipHYDROGEN];
987  }
988  }
989 
990  /* >>chng 05 mar 26, add this convergence test - important for molecular or advective
991  * sims since iso ion solver must sync up with chemistry */
992  /* remember current ground state population - will check if converged */
993  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
994  {
995  for( nelem=ipISO; nelem<LIMELM;++nelem )
996  {
997  if( dense.lgElmtOn[nelem] )
998  {
999  /* only do this check for "significant" levels of ionization */
1000  /*lint -e644 var possibly not init */
1001  if( dense.xIonDense[nelem][nelem+1-ipISO]/dense.gas_phase[nelem] > 1e-5 )
1002  {
1003  if( fabs(StatesElem[ipISO][nelem][0].Pop-save_iso_grnd[ipISO][nelem])/SDIV(StatesElem[ipISO][nelem][0].Pop)-1. >
1005  {
1006  conv.lgConvIoniz = false;
1007  sprintf( conv.chConvIoniz,"iso %3li%li",ipISO, nelem );
1008  conv.BadConvIoniz[0] = save_iso_grnd[ipISO][nelem]/SDIV(StatesElem[ipISO][nelem][0].Pop);
1009  fixit(); // this looks like a bug, it Pop/Pop, unity for all Pop>SMALLFLOAT
1010  conv.BadConvIoniz[1] = StatesElem[ipISO][nelem][0].Pop/SDIV(StatesElem[ipISO][nelem][0].Pop);
1011  }
1012  }
1013  /*lint +e644 var possibly not init */
1014  }
1015  }
1016  }
1017 
1018  /* this counts how many times ConvBase has been called in this iteration,
1019  * located here at bottom of routine so that number is false on first
1020  * call, set to 0 in when iteration starts - used to create itr/zn
1021  * number in printout often used to tell whether this is the very
1022  * first attempt at solution in this iteration */
1023  ++conv.nTotalIoniz;
1024 
1025  /* this was set with STOP NTOTALIONIZ command - only a debugging aid
1026  * by default is zero so false */
1028  {
1029  /* non-zero return indicates abort condition */
1030  fprintf(ioQQQ , "ABORT flag set since STOP nTotalIoniz was set and reached.\n");
1031  return 1;
1032  }
1033 
1035  {
1037  {
1038  static int iter_punch=-1;
1039  if( iteration !=iter_punch )
1040  fprintf( punch.ipDRout, "%s\n",punch.chHashString );
1041  iter_punch = iteration;
1042  }
1043 
1044  fprintf( punch.ipTraceConvergeBase,
1045  "%li\t%.4e\t%.4e\t%.4e\n",
1047  }
1048 
1049  return 0;
1050 }

Generated for cloudy by doxygen 1.8.1.1