cloudy  trunk
conv_base.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 /*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 "save.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 #include "h2.h"
32 #include "deuterium.h"
33 
34 STATIC bool lgNetEdenSrcSmall( void );
35 
36 void UpdateUTAs( void );
37 
38 /*lgIonizConverg check whether ionization of element nelem has converged, called by ionize,
39  * returns true if element is converged, false if not */
40 namespace
41 {
42  class IonizConverg
43  {
44  /* this is fractions [ion stage][nelem], ion stage = 0 for atom, nelem=0 for H*/
45  double OldFracs[LIMELM+1][LIMELM+1];
46  public:
47  IonizConverg()
48  {
49  for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
50  {
51  for (long ion = 0; ion <= nelem+1; ++ion)
52  {
53  OldFracs[nelem][ion] = dense.xIonDense[nelem][ion];
54  }
55  }
56  }
57  void operator()(
58  /* atomic number on the C scale, 0 for H, 25 for Fe */
59  long loop_ion ,
60  /* this is allowed error as a fractional change. Most are 0.15 */
61  double delta ,
62  /* option to print abundances */
63  bool lgPrint )
64  {
65  bool lgConverg_v = false;
66  double Abund,
67  bigchange ,
68  change ,
69  one;
70 
71  DEBUG_ENTRY( "IonizConverg::operator()" );
72 
73  for (long nelem =ipHYDROGEN; nelem<LIMELM; ++nelem )
74  {
75  if( !dense.lgElmtOn[nelem] )
76  continue;
77 
78  double abundold=0. , abundnew=0.;
79  long ionchg=-1;
80 
81  /* arguments are atomic number, ionization stage
82  * OldFracs[nelem][0] is abundance of nelem (cm^-3) */
83 
84  /* this function returns true if ionization of element
85  * with atomic number nelem has not changed by more than delta,*/
86 
87  /* check whether abundances exist yet, only do this for after first zone */
88  /*if( nzone > 0 )*/
89  /* >>chng 03 sep 02, check on changed ionization after first call through,
90  * to insure converged constant eden / temperature models */
91  if( conv.nPres2Ioniz )
92  {
93  /* >>chng 04 aug 31, this had been static, caused last hits on large changes
94  * in ionization */
95  bigchange = 0.;
96  change = 0.;
97  Abund = dense.gas_phase[nelem];
98 
99  /* loop over all ionization stages, loop over all ions, not just active ones,
100  * since this also sets old values, and so will zero out non-existant ones */
101  for( long ion=0; ion <= (nelem+1); ++ion )
102  {
103  /*lint -e727 OlsdFracs not initialized */
104  if( OldFracs[nelem][ion]/Abund > 1e-4 &&
105  dense.xIonDense[nelem][ion]/Abund > 1e-4 )
106  /*lint +e727 OlsdFracs not initialized */
107  {
108  /* check change in old to new value */
109  one = fabs(dense.xIonDense[nelem][ion]-OldFracs[nelem][ion])/
110  OldFracs[nelem][ion];
111  change = MAX2(change, one );
112  /* remember abundances for largest change */
113  if( change>bigchange )
114  {
115  bigchange = change;
116  abundold = OldFracs[nelem][ion]/Abund;
117  abundnew = dense.xIonDense[nelem][ion]/Abund;
118  ionchg = ion;
119  }
120  }
121  /* now set new value */
122  OldFracs[nelem][ion] = dense.xIonDense[nelem][ion];
123  }
124 
125  if( change >= delta )
126  {
127  char chConvIoniz[INPUT_LINE_LENGTH];
128  sprintf( chConvIoniz , "%2s ion" , elementnames.chElementSym[nelem] );
129  conv.setConvIonizFail(chConvIoniz,abundold,abundnew);
130  ASSERT( abundold>0. && abundnew>0. );
131  }
132  else
133  lgConverg_v = true;
134  }
135  else
136  {
137  for( long ion=0; ion <= (nelem+1); ++ion )
138  {
139  OldFracs[nelem][ion] = dense.xIonDense[nelem][ion];
140  }
141 
142  }
143 
144  /* option to print abundances */
145  if( lgPrint )
146  {
147  fprintf(ioQQQ," nz %ld loop %ld element %li converged? %c worst %ld change %g\n",
148  nzone, loop_ion, nelem, TorF(lgConverg_v),ionchg,bigchange);
149  for( long ion=0; ion<(nelem+1); ++ion )
150  {
151  fprintf(ioQQQ,"\t%.5e", dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]);
152  }
153  fprintf(ioQQQ,"\n");
154  }
155  }
156  }
157  };
158 }
159 
160 /*ConvBase main routine to drive ionization solution for all species, find total opacity
161  * called by ConvIoniz
162  * return 0 if ok, 1 if abort */
163 int ConvBase(
164  /* this is zero the first time ConvBase is called by convIoniz,
165  * counts number of call thereafter */
166  long loopi )
167 {
168  double HeatOld,
169  EdenTrue_old,
170  EdenFromMolecOld,
171  EdenFromGrainsOld,
172  HeatingOld ,
173  CoolingOld;
174  static double SecondOld;
175  static long int nzoneOTS=-1;
176  const int LOOP_ION_LIMIT = 10;
177  long int loop_ion;
178  static double SumOTS=0. , OldSumOTS[2]={0.,0.};
179  double save_iso_grnd[NISO][LIMELM];
180  valarray<realnum> mole_save(mole_global.num_calc);
181  IonizConverg lgIonizConverg;
182 
183  DEBUG_ENTRY( "ConvBase()" );
184 
185  /* this is set to phycon.te in tfidle, is used to insure that all temp
186  * vars are properly updated when conv_ionizeopacitydo is called
187  * NB must be same type as phycon.te */
188  ASSERT( fp_equal( phycon.te, thermal.te_update ) );
189 
190  /* this allows zone number to be printed with slight increment as zone converged
191  * conv.nPres2Ioniz is incremented at the bottom of this routine */
192  fnzone = (double)nzone + (double)conv.nPres2Ioniz/100.;
193 
194  /* reevaluate pressure */
195  /* this sets values of pressure.PresTotlCurr, also calls tfidle,
196  * and sets the total energy content of gas, which may be important for acvection */
197  PresTotCurrent();
198 
199  /* >>chng 04 sep 15, move EdenTrue_old up here, and will redo at bottom
200  * to find change
201  * find and save current true electron density - if this changes by more than the
202  * tolerance then ionization solution is not converged */
203  /* >>chng 04 jul 27, move eden_sum here from after this loop, so that change in EdenTrue
204  * can be monitored */
205  /* update EdenTrue, eden itself is actually changed in ConvEdenIoniz */
206  /* must not call eden_sum on very first time since for classic PDR total
207  * ionization may still be zero on first call */
208  if( conv.nTotalIoniz )
209  {
210  if( eden_sum() )
211  {
212  /* non-zero return indicates abort condition */
213  ++conv.nTotalIoniz;
214  return 1;
215  }
216  }
217 
218  /* the following two were evaluated in eden_sum
219  * will confirm that these are converged */
220  EdenTrue_old = dense.EdenTrue;
221  EdenFromMolecOld = mole.elec;
222  EdenFromGrainsOld = gv.TotalEden;
223  HeatingOld = thermal.htot;
224  CoolingOld = thermal.ctot;
225 
226  /* remember current ground state population - will check if converged */
227  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
228  {
229  for( long nelem=ipISO; nelem<LIMELM;++nelem )
230  {
231  if( dense.lgElmtOn[nelem] )
232  {
233  /* save the ground state population */
234  save_iso_grnd[ipISO][nelem] = iso_sp[ipISO][nelem].st[0].Pop();
235  }
236  }
237  }
238 
239  for( long i=0; i < mole_global.num_calc; i++ )
240  {
241  mole_save[i] = (realnum) mole.species[i].den;
242  }
243 
244  if( loopi==0 )
245  {
246  /* these will be used to look for oscillating ots rates */
247  OldSumOTS[0] = 0.;
248  OldSumOTS[1] = 0.;
249  conv.lgOscilOTS = false;
250  }
251 
252  if( trace.lgTrace )
253  {
254  fprintf( ioQQQ,
255  " ConvBase called. %.2f Te:%.3e HI:%.3e HII:%.3e H2:%.3e Ne:%.3e htot:%.3e CSUP:%.2e Conv?%c\n",
256  fnzone,
257  phycon.te,
258  dense.xIonDense[ipHYDROGEN][0],
259  dense.xIonDense[ipHYDROGEN][1],
260  hmi.H2_total,
261  dense.eden,
262  thermal.htot,
263  secondaries.csupra[ipHYDROGEN][0] ,
264  TorF(conv.lgConvIoniz()) );
265  }
266  /* want this flag to be true when we exit, various problems will set falst */
267  conv.resetConvIoniz();
268 
269  /* this routine is in heatsum.c, and zeros out the heating array */
270  HeatZero();
271 
272  /* if this is very first call, say not converged, since no meaningful way to
273  * check on changes in quantities. this counter is false only on very first
274  * call, when equal to zero */
275  if( !conv.nTotalIoniz )
276  {
277  conv.setConvIonizFail( "first call", 0.0, 0.0 );
278  }
279 
280  /* this will be flag to check whether any ionization stages
281  * were trimmed down */
282  conv.lgIonStageTrimed = false;
283 
284  /* must redo photoionization rates for all species on second and later tries */
285  /* always reevaluate the rates when . . . */
286  /* this flag says not to redo opac and photo rates, and following test will set
287  * true if one of several tests not done*/
288  opac.lgRedoStatic = false;
289  if(
290  /* opac.lgOpacStatic (usually true), is set false with no static opacity command */
291  !opac.lgOpacStatic ||
292  /* we are in search mode */
293  conv.lgSearch ||
294  /* this is the first call to this zone */
295  conv.nPres2Ioniz == 0 )
296  {
297  /* we need to redo ALL photoionization rates */
298  opac.lgRedoStatic = true;
299  }
300 
301  /* calculate radiative and dielectronic recombination rate coefficients */
302  ion_recom_calculate();
303 
304  /* this adjusts the lowest and highest stages of ionization we will consider,
305  * only safe to call when lgRedoStatic is true since this could lower the
306  * lowest stage of ionization, which needs all its photo rates */
307 
308  /* conv.nTotalIoniz is only 0 (false) on the very first call to here,
309  * when the ionization distribution is not yet done */
310  if( conv.nTotalIoniz )
311  {
312  bool lgIonizTrimCalled = false;
313  static long int nZoneCalled = 0;
314 
315  fixit(); // nZoneCalled should be reinitialized for each grid point?
316 
317  /* ionization trimming only used for He and heavier, not H */
318  /* only do this one time per zone since up and down cycle can occur */
319  /* >>chng 05 jan 15, increasing temperature above default first conditions, also
320  * no trim during search - this fixed major logical error when sim is
321  * totally mechanically heated to coronal temperatures -
322  * problem discovered by Ronnie Hoogerwerf */
323  /* do not keep changing the trim after the first call within
324  * this zone once we are deep into layer - doing so introduced very
325  * small level of noise as some stages
326  * went up and down - O+2 - O+3 was good example, when small H+3 after He+ i-front
327  * limit to one increase per element per zone */
328  if( conv.nTotalIoniz>2 &&
329  /* only call one time per zone except during search phase,
330  * when only call after 20 times only if temperature has changed */
331  ( conv.lgSearch || nZoneCalled!=nzone) )
332  {
333  lgIonizTrimCalled = true;
334  for( long nelem=ipHELIUM; nelem<LIMELM; ++nelem )
335  {
336  if( dense.lgElmtOn[nelem] )
337  {
338  /* ion_trim will set conv.lgIonStageTrimed true is any ion has its
339  * lowest stage of ionization dropped or trimmed */
340  ion_trim(nelem);
341  }
342  }
343  nZoneCalled = nzone;
344  }
345 
346  /* following block only set of asserts */
347 # if !defined(NDEBUG)
348  /* check that proper abundances are either positive or zero */
349  for( long nelem=ipHELIUM; nelem<LIMELM; ++nelem)
350  {
351  if( dense.lgElmtOn[nelem] )
352  {
353  for( long ion=0; ion<dense.IonLow[nelem]; ++ion )
354  {
355  ASSERT( dense.xIonDense[nelem][ion] == 0. );
356  }
357  /*if( nelem==5 ) fprintf(ioQQQ,"carbbb\t%li\n", dense.IonHigh[nelem]);*/
358  for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
359  {
360  /* >>chng 02 feb 06, had been > o., chng to > SMALLFLOAT to
361  * trip over VERY small floats that failed on alphas, but not 386
362  *
363  * in case where lower ionization stage was just lowered or
364  * trimmed down the abundance
365  * was set to SMALLFLOAT so test must be < SMALLFLOAT */
366  /* >>chng 02 feb 19, add check for search phase. During this search
367  * models with extreme ionization (all neutral or all ionized) can
368  * have extreme but non-zero abundances far from the ionization peak for
369  * element with lots of electrons. These will go away once the model
370  * becomes stable */
371  /* >>chng 03 dec 01, add check on whether ion trim was called
372  * conserve.in threw assert when iontrim not called and abund grew small */
373  ASSERT( conv.lgSearch || !lgIonizTrimCalled ||
374  /* this can happen if all C is in the form of CO */
375  (ion==0 && dense.IonHigh[nelem]==0 ) ||
376  dense.xIonDense[nelem][ion] >= SMALLFLOAT ||
377  dense.xIonDense[nelem][ion]/dense.gas_phase[nelem] >= SMALLFLOAT );
378  }
379  for( long ion=dense.IonHigh[nelem]+1; ion<nelem+1; ++ion )
380  {
381  ASSERT( ion >= 0 );
382  ASSERT( dense.xIonDense[nelem][ion] == 0. );
383  }
384  }
385  }
386 # endif
387  }
388 
389  /* now check if anything trimmed down */
390  if( conv.lgIonStageTrimed )
391  {
392  /* something was trimmed down, so say that ionization not yet stable */
393  /* say that ionization has not converged, secondaries changed by too much */
394  conv.setConvIonizFail( "IonTrimmed", 0.0, 0.0 );
395  }
396 
397  /* reevaluate advective terms if turned on */
398  if( dynamics.lgAdvection )
399  DynaIonize();
400 
401  /* evaluate Compton heating, bound E Compton, cosmic rays */
402  highen();
403 
404  // Depends on dense.eden, phycon.te and trimming level of H, He
405  iso_collapsed_update();
406 
407  iso_update_rates( );
408 
409  for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
410  {
411  (*diatom)->CalcPhotoionizationRate();
412  }
413 
414  long ionLowCache[LIMELM], ionHighCache[LIMELM];
415  for( long nelem=ipHYDROGEN; nelem<LIMELM; nelem++ )
416  {
417  if( dense.lgSetIoniz[nelem] )
418  {
419  dense.IonLow[nelem] = 0;
420  dense.IonHigh[nelem] = nelem + 1;
421  while( dense.SetIoniz[nelem][dense.IonLow[nelem]] == 0. )
422  ++dense.IonLow[nelem];
423  while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] == 0. )
424  --dense.IonHigh[nelem];
425  }
426  ionLowCache[nelem] = dense.IonLow[nelem];
427  ionHighCache[nelem] = dense.IonHigh[nelem];
428  }
429 
430  /* >>chng 04 feb 15, add loop over ionization until converged. Non-convergence
431  * in this case means ionization ladder pop changed, probably because of way
432  * that Auger ejection is treated - loop exits when conditions tested at end are met */
433  conv.incrementCounter(CONV_BASE_CALLS);
434  loop_ion = 0;
435  do
436  {
437  conv.incrementCounter(CONV_BASE_LOOPS);
438  ASSERT(lgElemsConserved());
439 
440  /* set the convergence flag to true,
441  * if anything changes in ConvBase, it will be set false */
442  if( loop_ion )
443  conv.resetConvIoniz();
444 
445  /* charge transfer evaluation needs to be here so that same rate
446  * coefficient used for H ion and other recombination */
447  /* fill in master charge transfer array, and zero out some arrays that track effects */
448 
449  ChargTranEval();
450 
451  /* find grain temperature, charge, and gas heating rate */
452  /* >>chng 04 apr 15, moved here from after call to HeLike(), PvH */
453  GrainDrive();
454 
455  /* evaluate Compton heating, bound E Compton, cosmic rays */
456  highen();
457 
458  /* find corrections for three-body rec - collisional ionization */
459  atmdat_3body();
460 
461  /* update UTA inner-shell ionization rates */
462  UpdateUTAs();
463 
464  ASSERT(lgElemsConserved());
465 
466  long ion_loop = 0;
467  const int nconv = 3;
468  double xIonDense0[nconv][LIMELM][LIMELM+1];
469  bool lgShortCircuit = false;
470  for( ion_loop=0; ion_loop<nconv && !lgShortCircuit; ++ion_loop)
471  {
472  conv.incrementCounter(CONV_BASE_ACCELS);
473  for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
474  {
475  for (long ion = 0; ion <= nelem+1; ++ion)
476  {
477  xIonDense0[ion_loop][nelem][ion] = dense.xIonDense[nelem][ion];
478  }
479  }
480 
481  double netion = 0.0;
482  for( long nelem=ipHYDROGEN; nelem<LIMELM; nelem++ )
483  {
484  if (!dense.lgElmtOn[nelem])
485  continue;
486 
487  ASSERT(dense.IonLow[nelem] >= ionLowCache[nelem]);
488  ASSERT(dense.IonHigh[nelem] <= ionHighCache[nelem]);
489  ion_wrapper( nelem );
490  if (nelem == ipHYDROGEN || nelem == ipOXYGEN)
491  {
492  double hion =
493  atmdat.CharExcRecTo[ipHYDROGEN][ipOXYGEN][0]*
494  iso_sp[ipH_LIKE][ipHYDROGEN].st[0].Pop()*
495  dense.xIonDense[ipOXYGEN][1]-
496  atmdat.CharExcIonOf[ipHYDROGEN][ipOXYGEN][0]*
497  dense.xIonDense[ipHYDROGEN][1]*
498  dense.xIonDense[ipOXYGEN][0];
499  if (nelem == ipHYDROGEN)
500  netion += hion;
501  else
502  netion -= hion;
503  }
504 
505  if ( conv.lgUpdateCouplings && conv.nTotalIoniz )
506  mole_update_sources();
507  }
508  if (dense.lgElmtOn[ipHYDROGEN] && dense.lgElmtOn[ipOXYGEN] && ion_loop == 1)
509  {
510  // Comparison rate is whether error is either a small fraction of the minimum ionization rate,
511  // or a tighter check on whether the fluxes are balanced as well as is feasible numerically
512  double ion_cmp = MAX2(0.01*MIN2(ionbal.RateIonizTot(ipHYDROGEN,0), ionbal.RateIonizTot(ipOXYGEN,0)),
513  1e-4*atmdat.CharExcRecTo[ipHYDROGEN][ipOXYGEN][0]*
514  iso_sp[ipH_LIKE][ipHYDROGEN].st[0].Pop()*dense.xIonDense[ipOXYGEN][1]);
515 
516  if (fabs(netion) > ion_cmp &&
517  dense.xIonDense[ipOXYGEN][1]*dense.xIonDense[ipHYDROGEN][1] >
518  1e-8*iso_sp[ipH_LIKE][ipHYDROGEN].st[0].Pop()*dense.xIonDense[ipOXYGEN][0] )
519  {
520  conv.setConvIonizFail( "OH CX inconsistency" , ion_cmp, netion);
521  }
522  }
523 
524  // If not going to end anyhow, check whether changes were small
525  bool lgCanShortCircuit = (ion_loop+1 < nconv);
526  for( long nelem=ipHYDROGEN; nelem<LIMELM && lgCanShortCircuit; ++nelem )
527  {
528  if (!dense.lgElmtOn[nelem])
529  continue;
530  for (long ion = dense.IonLow[nelem];
531  ion <= dense.IonHigh[nelem]; ++ion)
532  {
533  double x0 = xIonDense0[ion_loop][nelem][ion];
534  double x1 = dense.xIonDense[nelem][ion];
535  if (fabs(x0-x1) > 1e-6*(x0+x1))
536  {
537  lgCanShortCircuit = false;
538  break;
539  }
540  }
541  }
542  lgShortCircuit = lgCanShortCircuit;
543  }
544 
545  if (!lgShortCircuit)
546  {
547  // Apply convergence acceleration
548  for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
549  {
550  if (!dense.lgElmtOn[nelem])
551  continue;
552  double tot0 = 0., tot1 = 0.;
553  double xIonNew[LIMELM+1];
554  for (long ion = 0; ion <= nelem+1; ++ion)
555  {
556  double x0 = xIonDense0[nconv-2][nelem][ion];
557  double x1 = xIonDense0[nconv-1][nelem][ion];
558  double x2 = dense.xIonDense[nelem][ion];
559  xIonNew[ion] = x2;
560  tot0 += x2;
561  // Richardson extrapolation formula to accelerate convergence
562  // Assumes convergence to x^ is geometric, i.e.
563  // x_i = x^ + a d^i
564 
565  double extstep = 0.,predict=x2,
566  step0 = x1-x0, step1 = x2-x1, abs1 = fabs(step1);
567  // Protect against roundoff/noise in inner solver
568  if ( abs1 > 1000.0*((double)DBL_EPSILON)*x2 )
569  {
570  double denom = fabs(step1-step0);
571  double sgn = (step1*step0 > 0)? 1.0 : -1.0;
572  // Greatest acceleration allowed is MAXACC*latest step length
573  // Can we do better than static tuning this parameter?
574  const double MAXACC=100.0;
575  double extfac = 1.0/(denom/abs1 + 1.0/MAXACC);
576  extstep = sgn*extfac*step1;
577  //extstep = sgn*MIN2(extfac*step1,0.01*x2);
578  predict = x2+extstep;
579  if (predict > 0.0)
580  xIonNew[ion] = predict;
581  }
582  if ( 0 )
583  //if ( nelem == ipHYDROGEN || (nelem == ipIRON && ion <=2 ) )
584  if ( (nelem == ipNICKEL && ion <=2 ) )
585  fprintf(ioQQQ,"Extrap %3ld %3ld %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g\n",
586  nelem,ion,
587  x0,x0-xIonDense0[nconv-3][nelem][ion],x1-x0,x2-x1,extstep,predict);
588  tot1 += xIonNew[ion];
589  }
590  if ( tot1 > SMALLFLOAT )
591  {
592  double scal = tot0/tot1;
593  for (long ion = 0; ion <= nelem+1; ++ion)
594  {
595  dense.xIonDense[nelem][ion] = scal*xIonNew[ion];
596  }
597  for ( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
598  {
599  double renorm;
600  iso_renorm(nelem, ipISO, renorm);
601  }
602  }
603  }
604 
605  bool lgPostExtrapSolve = true;
606  if (lgPostExtrapSolve)
607  {
608  for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
609  {
610  if (!dense.lgElmtOn[nelem])
611  continue;
612  ion_wrapper( nelem );
613  }
614  }
615  }
616  SetDeuteriumIonization( dense.xIonDense[ipHYDROGEN][0], dense.xIonDense[ipHYDROGEN][1] );
617 
618  /*>>chng 04 may 09, add option to abort here, inspired by H2 pop failures
619  * which cascaded downstream if we did not abort */
620  /* return if one of above solvers has declared abort condition */
621  if( lgAbort )
622  {
623  ++conv.nTotalIoniz;
624  return 1;
625  }
626 
627  ASSERT(lgElemsConserved());
628 
629  /* drive chemistry network*/
630  mole_drive();
631 
632  ASSERT(lgElemsConserved());
633 
634  /* all elements have now had ionization reevaluated - in some cases we may have upset
635  * the ionization that was forced with an "element ionization" command - here we will
636  * re-impose that set ionization */
637  /* >>chng 04 oct 13, add this logic */
638  for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
639  {
640  if( dense.lgSetIoniz[nelem] )
641  {
642  dense.IonLow[nelem] = 0;
643  dense.IonHigh[nelem] = nelem + 1;
644  while( dense.SetIoniz[nelem][dense.IonLow[nelem]] == 0. )
645  ++dense.IonLow[nelem];
646  while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] == 0. )
647  --dense.IonHigh[nelem];
648  }
649  }
650 
651  /* redo ct ion rate for reasons now unclear */
652  ChargTranEval();
653 
654  /* evaluate molecular hydrogen level populations */
655  for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
656  {
657  bool lgPopsConverged = true;
658  double old_val, new_val;
659  (*diatom)->H2_LevelPops( lgPopsConverged, old_val, new_val );
660  if( !lgPopsConverged )
661  {
662  conv.setConvIonizFail( "H2 pops", old_val, new_val);
663  }
664  }
665 
666  ASSERT(lgElemsConserved());
667 
668  /* lgIonizConverg is a function to check whether ionization has converged
669  * check whether ionization changed by more than relative error
670  * given by second number */
671  /* >>chng 04 feb 14, loop over all elements rather than just a few */
672  lgIonizConverg(loop_ion, conv.IonizErrorAllowed , false );
673 
674  if( deut.lgElmtOn )
675  {
676  static double OldDeut[2] = {0., 0.};
677  for( long ion=0; ion<2; ++ion )
678  {
679  if( fabs(deut.xIonDense[ion] - OldDeut[ion] ) > 0.2*conv.IonizErrorAllowed*fabs(OldDeut[ion]) )
680  {
681  conv.setConvIonizFail( "D ion" , OldDeut[ion], deut.xIonDense[ion]);
682  }
683  OldDeut[ion] = deut.xIonDense[ion];
684  }
685  }
686 
687  if( fabs(EdenTrue_old - dense.EdenTrue) > conv.EdenErrorAllowed/2.*fabs(dense.EdenTrue) )
688  {
689  conv.setConvIonizFail( "EdTr cng" , EdenTrue_old, dense.EdenTrue);
690  }
691 
692  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
693  {
694  for( long nelem=ipISO; nelem<LIMELM; ++nelem )
695  {
696  lgStatesConserved( nelem, nelem-ipISO, iso_sp[ipISO][nelem].st, iso_sp[ipISO][nelem].numLevels_local, 1e-3f, loop_ion );
697  }
698  }
699 
700  if( trace.nTrConvg>=4 )
701  {
702  /* trace ionization */
703  fprintf( ioQQQ,
704  " ConvBase4 ionization driver loop_ion %li converged? %c reason not converged %s\n" ,
705  loop_ion ,
706  TorF( conv.lgConvIoniz()) ,
707  conv.chConvIoniz() );
708  }
709 
710  ++loop_ion;
711  }
712  /* loop is not converged, less than loop limit, and we are reevaluating */
713  while( !conv.lgConvIoniz() && loop_ion < LOOP_ION_LIMIT && rfield.lgIonizReevaluate);
714 
715  if( conv.lgConvIoniz() )
716  lgNetEdenSrcSmall();
717 
718  /* >>chng 05 oct 29, move CT heating here from heat_sum since this sometimes contributes
719  * cooling rather than heat and so needs to be sorted out before either heating or cooling
720  * are derived first find net heating - */
721  thermal.char_tran_heat = ChargTranSumHeat();
722  /* net is cooling if negative */
723  thermal.char_tran_cool = MAX2(0. , -thermal.char_tran_heat );
724  thermal.char_tran_heat = MAX2(0. , thermal.char_tran_heat );
725 
726  /* get total cooling, thermal.ctot = does not occur since passes as pointer. This can add heat.
727  * it calls coolSum at end to sum up the total cooling */
728  CoolEvaluate(&thermal.ctot );
729 
730  /* get total heating rate - first save old quantities to check how much it changes */
731  HeatOld = thermal.htot;
732 
733  /* HeatSum will update ElecFrac,
734  * secondary electron ionization and excitation efficiencies,
735  * and sum up current secondary rates - remember old values before we enter */
736  SecondOld = secondaries.csupra[ipHYDROGEN][0];
737 
738  /* update the total heating - it was all set to zero in HeatZero at top of this routine,
739  * occurs before secondaries bit below, since this updates electron fracs */
740  HeatSum();
741 
742  /* test whether we can just set the rate to the new one, or whether we should
743  * take average and do it again. secondaries.sec2total was set in hydrogenic, and
744  * is ratio of secondary to total hydrogen destruction rates */
745  /* >>chng 02 nov 20, add test on size of csupra - primal had very close to underflow */
746  if( (secondaries.csupra[ipHYDROGEN][0]>SMALLFLOAT) && (secondaries.sec2total>0.001) &&
747  fabs( 1. - SecondOld/SDIV(secondaries.csupra[ipHYDROGEN][0]) ) > 0.1 &&
748  SecondOld > 0. && secondaries.csupra[ipHYDROGEN][0] > 0.)
749  {
750  /* say that ionization has not converged, secondaries changed by too much */
751  conv.setConvIonizFail( "SecIonRate", SecondOld,
752  secondaries.csupra[ipHYDROGEN][0] );
753  }
754 
755 # if 0
756  static realnum hminus_old=0.;
757  /* >>chng 04 apr 15, add this convergence test */
758  if( conv.nTotalIoniz )
759  {
760  realnum hminus_den = findspecieslocal("H-")->den;
761  if( fabs( hminus_old-hminus_den ) > fabs( hminus_den * conv.EdenErrorAllowed ) )
762  {
763  conv.setConvIonizFail( "Big H- chn", hminus_old, hminus_den );
764  }
765  hminus_old = hminus_den;
766  }
767 # endif
768 
769  if( HeatOld > 0. && !thermal.lgTemperatureConstant )
770  {
771  /* check if heating has converged - tolerance is final match */
772  if( fabs(1.-thermal.htot/HeatOld) > conv.HeatCoolRelErrorAllowed*.5 )
773  {
774  conv.setConvIonizFail( "Big d Heat", HeatOld, thermal.htot);
775  }
776  }
777 
778  /* abort flag may have already been set - if so bail */
779  if( lgAbort )
780  {
781 
782  return 1;
783  }
784 
785  /* evaluate current opacities
786  * rfield.lgOpacityReevaluate normally true,
787  * set false with no opacity reevaluate command, an option to only
788  * evaluate opacity one time per zone */
789  if( !conv.nPres2Ioniz || rfield.lgOpacityReevaluate || conv.lgIonStageTrimed )
790  OpacityAddTotal();
791 
792  /* >>chng 02 jun 11, call even first time that this routine is called -
793  * this seemed to help convergence */
794 
795  /* do OTS rates for all lines and all continua since
796  * we now have ionization balance of all species. Note that this is not
797  * entirely self-consistent, since destruction probabilities here are not the same as
798  * the ones used in the model atoms. Problems?? if near convergence
799  * then should be nearly identical */
800  if( !conv.nPres2Ioniz || rfield.lgOpacityReevaluate || rfield.lgIonizReevaluate ||
801  conv.lgIonStageTrimed || conv.lgSearch || nzone!=nzoneOTS )
802  {
803  RT_OTS();
804  nzoneOTS = nzone;
805 
806  /* remember old ots rates */
807  OldSumOTS[0] = OldSumOTS[1];
808  OldSumOTS[1] = SumOTS;
809  /*fprintf(ioQQQ," calling RT_OTS zone %.2f SumOTS is %.2e\n",fnzone,SumOTS);*/
810 
811  /* now update several components of the continuum, this only happens after
812  * we have gone through entire solution for this zone at least one time.
813  * there can be wild ots oscillation on first call */
814  /* the rel change of 0.2 was chosen by running hizqso - smaller increased
815  * itrzn but larger did not further decrease it. */
816  RT_OTS_Update(&SumOTS);
817  /*fprintf(ioQQQ,"RT_OTS_Updateee\t%.3f\t%.2e\t%.2e\n", fnzone,SumOTS , OldSumOTS[1] );*/
818  }
819  else
820  SumOTS = 0.;
821 
822  /* now check whether the ots rates changed */
823  if( SumOTS> 0. )
824  {
825  /* the ots rate must be converged to the error in the electron density */
826  /* how fine should this be converged?? originally had above, 10%, but take
827  * smaller ratio?? */
828  if( fabs(1.-OldSumOTS[1]/SumOTS) > conv.EdenErrorAllowed )
829  {
830  /* this branch, ionization not converged due to large change in ots rates.
831  * check whether ots rates are oscillating, if loopi > 1 so we have enough info*/
832  if( loopi > 1 )
833  {
834  /* here we have three values, are they changing sign? */
835  if( (OldSumOTS[0]-OldSumOTS[1]) * ( OldSumOTS[1] - SumOTS ) < 0. )
836  {
837  /* ots rates are oscillating */
838  conv.lgOscilOTS = true;
839  }
840  }
841 
842  conv.setConvIonizFail( "OTSRatChng" , OldSumOTS[1], SumOTS);
843  }
844 
845  /* produce info on the ots fields if either "trace ots" or
846  * "trace convergence xxx ots " was entered */
847  if( ( trace.lgTrace && trace.lgOTSBug ) ||
848  ( trace.nTrConvg && trace.lgOTSBug ) )
849  {
850  RT_OTS_PrtRate(SumOTS*0.05 , 'b' );
851  }
852  /*fprintf(ioQQQ,"DEBUG opac\t%.2f\t%.3e\t%.3e\n",fnzone,
853  dense.xIonDense[ipNICKEL][0] ,
854  dense.xIonDense[ipZINC][0] );*/
855  {
856  /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */
857  enum {DEBUG_LOC=false};
858  if( DEBUG_LOC && (nzone>110) )
859  {
860 # if 0
861 # include "lines_service.h"
862  DumpLine( &iso_sp[ipH_LIKE][ipHYDROGEN].trans(2,0) );
863 # endif
864  /* last par 'l' for line, 'c' for continua, 'b' for both,
865  * the numbers printed are:
866  * cell i, [i], so 1 less than ipoint
867  * anu[i],
868  * otslin[i],
869  * opacity_abs[i],
870  * otslin[i]*opacity_abs[i],
871  * rfield.chLineLabel[i] ,
872  * rfield.line_count[i] */
873  }
874  }
875  }
876  {
877  /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */
878  enum {DEBUG_LOC=false};
879  if( DEBUG_LOC && (nzone>200) )
880  {
881  fprintf(ioQQQ,"debug otsss\t%li\t%.3e\t%.3e\t%.3e\n",
882  nzone,
883  iso_sp[0][1].trans(15,3).Emis().ots(),
884  TauLines[26].Emis().ots(),
885  opac.opacity_abs[2069]);
886  }
887  }
888 
889  /* option to print OTS continuum with TRACE OTS */
890  if( trace.lgTrace && trace.lgOTSBug )
891  {
892  /* find ots rates here, so we only print fraction of it,
893  * SumOTS is both line and continuum contributing to ots, and is multiplied by opacity */
894  /* number is weakest rate to print */
895  RT_OTS_PrtRate(SumOTS*0.05 , 'b' );
896  }
897 
898  // all RT routines called
899  RT_line_all( );
900 
901  /* >>chng 01 mar 16, evaluate pressure here since changing and other values needed */
902  /* reevaluate pressure */
903  /* this sets values of pressure.PresTotlCurr, also calls tfidle */
904  PresTotCurrent();
905 
906  ASSERT(lgElemsConserved());
907 
908  /* update some counters that keep track of how many times this routine
909  * has been called */
910 
911  if( trace.lgTrace )
912  {
913  fprintf( ioQQQ,
914  " ConvBase return. fnzone %.2f nPres2Ioniz %li Te:%.3e HI:%.3e HII:%.3e H2:%.3e Ne:%.3e htot:%.3e CSUP:%.2e Conv?%c reason:%s\n",
915  fnzone,
916  conv.nPres2Ioniz,
917  phycon.te,
918  dense.xIonDense[ipHYDROGEN][0],
919  dense.xIonDense[ipHYDROGEN][1],
920  hmi.H2_total,
921  dense.eden,
922  thermal.htot,
923  secondaries.csupra[ipHYDROGEN][0] ,
924  TorF(conv.lgConvIoniz()) ,
925  conv.chConvIoniz());
926  }
927 
928  /* this counts number of times we are called by ConvPresTempEdenIoniz,
929  * number of calls in this zone so first call is zero
930  * reset to zero each time ConvPresTempEdenIoniz is called */
931  ++conv.nPres2Ioniz;
932 
933  /* this is abort option set with SET PRES IONIZ command,
934  * test on nzone since init can take many iterations
935  * this is seldom used except in special cases */
936  if( conv.nPres2Ioniz > conv.limPres2Ioniz && nzone > 0)
937  {
938  fprintf(ioQQQ," PROBLEM ConvBase sets lgAbort since nPres2Ioniz exceeds limPres2Ioniz. ");
939  fprintf(ioQQQ,"Their values are %li and %li.\n",conv.nPres2Ioniz , conv.limPres2Ioniz);
940  fprintf(ioQQQ," Reset this limit with the SET PRES IONIZ command.\n");
941  lgAbort = true;
942  return 1;
943  }
944 
945  /* various checks on the convergence of the current solution */
946  if( eden_sum() )
947  {
948  /* non-zero return indicates abort condition */
949  return 1;
950  }
951 
952  /* is electron density converged? */
954  if( fabs(EdenTrue_old - dense.EdenTrue) > fabs(dense.EdenTrue * conv.EdenErrorAllowed/2.) )
955  {
956  conv.setConvIonizFail( "eden chng", EdenTrue_old, dense.EdenTrue);
957  }
958 
959  /* check on molecular electron den */
960  if( fabs(EdenFromMolecOld - mole.elec) > fabs(dense.EdenTrue * conv.EdenErrorAllowed/2.) )
961  {
962  conv.setConvIonizFail( "edn chnCO", EdenFromMolecOld, dense.EdenTrue);
963  }
964 
965  if( gv.lgGrainElectrons )
966  {
967  /* check on grain electron den */
968  if( fabs(EdenFromGrainsOld - gv.TotalEden) > fabs(dense.EdenTrue * conv.EdenErrorAllowed/2.) )
969  {
970  conv.setConvIonizFail( "edn grn e", EdenFromGrainsOld, gv.TotalEden);
971  }
972 
973  /* check on sum of grain and molecular electron den - often two large numbers that cancel */
974  if( fabs( (EdenFromMolecOld-EdenFromGrainsOld) - (mole.elec-gv.TotalEden) ) >
975  fabs(dense.EdenTrue * conv.EdenErrorAllowed/4.) )
976  {
977  conv.setConvIonizFail( "edn mole-grn",
978  (EdenFromMolecOld-EdenFromGrainsOld),
979  (mole.elec-gv.TotalEden));
980  }
981  }
982 
983  /* check on heating and cooling if vary temp model
984  * >>chng 08 jul 01, over the code's entire history it had tested whether
985  * this is a constant temperature simulation and did not do this test if
986  * the thermal solution was not done. There are some cases where we do
987  * want to specify the temperature and then find the heating or cooling -
988  * this is done in calculations of cooling curves for instance. With this
989  * change the heating/cooling are converged even in a constant temperature
990  * sim. this does make CT sims run more slowly but with greater accuracy
991  * if heating or cooling is reported */
992  if( !thermal.lgTemperatureConstant )
993  {
994  if( fabs(HeatingOld - thermal.htot)/thermal.htot > conv.HeatCoolRelErrorAllowed/2. )
995  {
996  conv.setConvIonizFail( "heat chg", HeatingOld, thermal.htot);
997  }
998 
999  if( fabs(CoolingOld - thermal.ctot)/thermal.ctot > conv.HeatCoolRelErrorAllowed/2. )
1000  {
1001  conv.setConvIonizFail( "cool chg", CoolingOld, thermal.ctot);
1002  }
1003  }
1004 
1005  /* check whether molecular abundances are stable */
1006  for( long i=0; i < mole_global.num_calc; ++i )
1007  {
1008  // done relative to total nuclei density so that we can do pure metal plasmas.
1009  if( fabs(mole.species[i].den-mole_save[i])/dense.xNucleiTotal-1. >
1010  conv.EdenErrorAllowed/2.)
1011  {
1012  char chConvIoniz[INPUT_LINE_LENGTH];
1013  sprintf( chConvIoniz, "ch %-4.4s",mole_global.list[i]->label.c_str() );
1014  conv.setConvIonizFail( chConvIoniz,
1015  mole_save[i]/dense.xNucleiTotal,
1016  mole.species[i].den/dense.xNucleiTotal);
1017  }
1018  }
1019 
1020  /* >>chng 05 mar 26, add this convergence test - important for molecular or advective
1021  * sims since iso ion solver must sync up with chemistry */
1022  /* remember current ground state population - will check if converged */
1023  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
1024  {
1025  for( long nelem=ipISO; nelem<LIMELM;++nelem )
1026  {
1027  if( dense.lgElmtOn[nelem] )
1028  {
1029  /* only do this check for "significant" levels of ionization */
1030  /*lint -e644 var possibly not init */
1031  if( dense.xIonDense[nelem][nelem-ipISO]/dense.gas_phase[nelem] > 1e-5 )
1032  {
1033  if( fabs(iso_sp[ipISO][nelem].st[0].Pop()-save_iso_grnd[ipISO][nelem])/SDIV(iso_sp[ipISO][nelem].st[0].Pop())-1. >
1034  conv.EdenErrorAllowed)
1035  {
1036  char chConvIoniz[INPUT_LINE_LENGTH];
1037  sprintf( chConvIoniz,"iso %2li %2li",ipISO, nelem );
1038  conv.setConvIonizFail(chConvIoniz,
1039  save_iso_grnd[ipISO][nelem],
1040  iso_sp[ipISO][nelem].st[0].Pop());
1041  }
1042  }
1043  /*lint +e644 var possibly not init */
1044  }
1045  }
1046  }
1047 
1048  /* this counts how many times ConvBase has been called in this iteration,
1049  * located here at bottom of routine so that number is false on first
1050  * call, set to 0 in when iteration starts - used to create itr/zn
1051  * number in printout often used to tell whether this is the very
1052  * first attempt at solution in this iteration */
1053  ++conv.nTotalIoniz;
1054 
1055  /* first sweep is over if this is not search phase */
1056  if( !conv.lgSearch )
1057  conv.lgFirstSweepThisZone = false;
1058 
1059  /* this was set with STOP NTOTALIONIZ command - only a debugging aid
1060  * by default is zero so false */
1061  if(StopCalc.nTotalIonizStop && conv.nTotalIoniz> StopCalc.nTotalIonizStop )
1062  {
1063  /* non-zero return indicates abort condition */
1064  fprintf(ioQQQ , "ABORT flag set since STOP nTotalIoniz was set and reached.\n");
1065  return 1;
1066  }
1067 
1068  if( save.lgTraceConvergeBase )
1069  {
1070  if( iteration > 1 && save.lgTraceConvergeBaseHash )
1071  {
1072  static int iter_punch=-1;
1073  if( iteration !=iter_punch )
1074  fprintf( save.ipTraceConvergeBase, "%s\n",save.chHashString );
1075  iter_punch = iteration;
1076  }
1077 
1078  fprintf( save.ipTraceConvergeBase,
1079  "%li\t%.4e\t%.4e\t%.4e\n",
1080  nzone , thermal.htot , thermal.ctot , dense.eden );
1081  }
1082 
1083  return 0;
1084 }
1085 
1086 void UpdateUTAs( void )
1087 {
1088  DEBUG_ENTRY( "UpdateUTAs()" );
1089 
1090  /* only reevaluate this on first pass through on this zone */
1091  if( conv.lgFirstSweepThisZone )
1092  {
1093  /* inner shell ionization */
1094  for( long nelem=0; nelem< LIMELM; ++nelem )
1095  {
1096  for( long ion=0; ion<nelem+1; ++ion )
1097  {
1098  ionbal.UTA_ionize_rate[nelem][ion] = 0.;
1099  ionbal.UTA_heat_rate[nelem][ion] = 0.;
1100  }
1101  }
1102  /* inner shell ionization by UTA lines */
1103  /* this flag turned off with no UTA command */
1104  if( ionbal.lgInnerShellLine_on )
1105  {
1106  for( long i=0; i < nUTA; ++i )
1107  {
1108  /* rateone is inverse lifetime of level against autoionization */
1109  double rateone = UTALines[i].Emis().pump() * UTALines[i].Emis().AutoIonizFrac();
1110  ionbal.UTA_ionize_rate[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1] += rateone;
1111  /* heating rate in erg atom-1 s-1 */
1112  ionbal.UTA_heat_rate[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1] += rateone*UTALines[i].Coll().heat();
1113  {
1114  /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */
1115  /*@-redef@*/
1116  enum {DEBUG_LOC=false};
1117  /*@+redef@*/
1118  if( DEBUG_LOC /*&& UTALines[i].nelem==ipIRON+1 && (UTALines[i].IonStg==15||UTALines[i].IonStg==14)*/ )
1119  {
1120  fprintf(ioQQQ,"DEBUG UTA %3i %3i %.3f %.2e %.2e %.2e\n",
1121  (*UTALines[i].Hi()).nelem() , (*UTALines[i].Hi()).IonStg() , UTALines[i].WLAng() ,
1122  rateone, UTALines[i].Coll().heat(),
1123  UTALines[i].Coll().heat()*dense.xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1] );
1124  }
1125  }
1126  }
1127  }
1128  }
1129 
1130  return;
1131 }
1132 
1133 STATIC bool lgNetEdenSrcSmall( void )
1134 {
1135  fixit(); // this routine needs to be enabled
1136  return true;
1137  DEBUG_ENTRY( "lgNetEdenSrcSmall()" );
1138 
1139  if( conv.lgSearch )
1140  return true;
1141  fixit(); // grain rates not well tested below
1142  if( gv.lgDustOn() )
1143  return true;
1144 
1145  // Check for consistency of explicit source and sink rates for
1146  // electrons with population derived from neutrality
1147  double ionsrc = 0., ionsnk = 0.;
1148  for( long nelem=0; nelem < LIMELM; ++nelem )
1149  {
1150  if( !dense.lgElmtOn[nelem] )
1151  continue;
1152  ionsrc += ionbal.elecsrc[nelem];
1153  ionsnk += ionbal.elecsnk[nelem];
1154  for( long ion_from = 0; ion_from <= nelem + 1; ++ion_from )
1155  {
1156  for( long ion_to = 0; ion_to <= nelem + 1; ++ion_to )
1157  {
1158  if( ion_to-ion_from > 0 )
1159  {
1160  ionsrc += gv.GrainChTrRate[nelem][ion_from][ion_to] *
1161  dense.xIonDense[nelem][ion_from] * (ion_to-ion_from);
1162  }
1163  else if( ion_to-ion_from < 0 )
1164  {
1165  ionsnk += gv.GrainChTrRate[nelem][ion_from][ion_to] *
1166  dense.xIonDense[nelem][ion_from] * (ion_from-ion_to);
1167  }
1168  }
1169  }
1170  }
1171  long ipMElec = findspecies("e-")->index;
1172  const double totsrc = ionsrc + mole.species[ipMElec].src;
1173  const double totsnk = ionsnk + mole.species[ipMElec].snk*dense.EdenTrue;
1174  const double diff = (totsrc - totsnk);
1175  const double ave = ( fabs(totsrc) + fabs(totsnk) )/2.;
1176  fixit(); // Need to tighten up e- population convergence criterion
1177  const double error_allowed = 0.05 * ave; //conv.EdenErrorAllowed * ave;
1178  if( fabs(diff) > error_allowed )
1179  {
1180  enum {DEBUG_LOC=false};
1181  if( DEBUG_LOC )
1182  {
1183  fprintf(ioQQQ,"PROBLEM large NetEdenSrc nzone %li\t%e\t%e\t%e\t%e\n",
1184  nzone,
1185  totsrc/SDIV(totsnk),
1186  dense.EdenTrue,
1187  ionsrc + mole.species[ipMElec].src,
1188  ionsnk + mole.species[ipMElec].snk*dense.EdenTrue);
1189  }
1190  conv.setConvIonizFail( "NetEdenSrc", diff, error_allowed);
1191  return false;
1192  }
1193  else
1194  return true;
1195 }
RT_line_all
void RT_line_all(void)
Definition: rt_line_all.cpp:26
thermal.h
t_atmdat::CharExcIonOf
double CharExcIonOf[NCX][LIMELM][LIMELM+1]
Definition: atmdat.h:152
findspecieslocal
molezone * findspecieslocal(const char buf[])
Definition: mole_species.cpp:833
TorF
char TorF(bool l)
Definition: cddefines.h:710
lgAbort
bool lgAbort
Definition: cddefines.cpp:10
ipOXYGEN
const int ipOXYGEN
Definition: cddefines.h:312
atmdat_3body
void atmdat_3body(void)
Definition: atmdat_3body.cpp:52
StopCalc
t_StopCalc StopCalc
Definition: stopcalc.cpp:5
t_rfield::lgOpacityReevaluate
bool lgOpacityReevaluate
Definition: rfield.h:121
t_save::lgTraceConvergeBase
bool lgTraceConvergeBase
Definition: save.h:340
t_dense::eden
double eden
Definition: dense.h:190
DumpLine
void DumpLine(const TransitionProxy &t)
Definition: transition.cpp:100
dense
t_dense dense
Definition: dense.cpp:24
elementnames.h
secondaries.h
t_dynamics::lgAdvection
bool lgAdvection
Definition: dynamics.h:60
ion_wrapper
void ion_wrapper(long nelem)
Definition: ion_solver.cpp:1498
rfield
t_rfield rfield
Definition: rfield.cpp:8
t_deuterium::lgElmtOn
bool lgElmtOn
Definition: deuterium.h:19
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
ConvBase
int ConvBase(long loopi)
Definition: conv_base.cpp:163
diatom_iter
vector< diatomics * >::iterator diatom_iter
Definition: h2.h:13
realnum
float realnum
Definition: cddefines.h:103
conv.h
rfield.h
t_conv::EdenErrorAllowed
double EdenErrorAllowed
Definition: conv.h:267
STATIC
#define STATIC
Definition: cddefines.h:97
GrainVar::TotalEden
double TotalEden
Definition: grainvar.h:528
t_conv::lgFirstSweepThisZone
bool lgFirstSweepThisZone
Definition: conv.h:155
nUTA
long int nUTA
Definition: taulines.cpp:26
mole.h
iso_collapsed_update
void iso_collapsed_update(void)
Definition: iso_solve.cpp:27
t_conv::IonizErrorAllowed
realnum IonizErrorAllowed
Definition: conv.h:280
t_mole_local::elec
double elec
Definition: mole.h:390
molecule::index
int index
Definition: mole.h:169
ChargTranSumHeat
double ChargTranSumHeat(void)
Definition: atmdat_char_tran.cpp:569
t_conv::setConvIonizFail
void setConvIonizFail(const char *reason, double oldval, double newval)
Definition: conv.h:107
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
PresTotCurrent
void PresTotCurrent(void)
Definition: pressure_total.cpp:34
t_thermal::char_tran_heat
double char_tran_heat
Definition: thermal.h:146
diatoms
vector< diatomics * > diatoms
Definition: h2.cpp:8
highen.h
phycon
t_phycon phycon
Definition: phycon.cpp:6
trace.h
t_mole_global::list
MoleculeList list
Definition: mole.h:317
t_ionbal::UTA_heat_rate
double ** UTA_heat_rate
Definition: ionbal.h:172
t_conv::lgOscilOTS
bool lgOscilOTS
Definition: conv.h:193
t_opac::lgOpacStatic
bool lgOpacStatic
Definition: opacity.h:140
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:952
t_dense::gas_phase
realnum gas_phase[LIMELM]
Definition: dense.h:71
lines_service.h
t_dense::EdenTrue
double EdenTrue
Definition: dense.h:221
CONV_BASE_CALLS
@ CONV_BASE_CALLS
Definition: conv.h:75
dynamics.h
t_conv::chConvIoniz
const char * chConvIoniz() const
Definition: conv.h:119
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
t_hmi::H2_total
double H2_total
Definition: hmi.h:16
t_ionbal::UTA_ionize_rate
double ** UTA_ionize_rate
Definition: ionbal.h:170
iso.h
HeatSum
void HeatSum(void)
Definition: heat_sum.cpp:33
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
t_thermal::ctot
double ctot
Definition: thermal.h:112
opac
t_opac opac
Definition: opacity.cpp:5
atmdat.h
t_trace::lgOTSBug
bool lgOTSBug
Definition: trace.h:103
lgNetEdenSrcSmall
STATIC bool lgNetEdenSrcSmall(void)
Definition: conv_base.cpp:1133
x1
static double x1[83]
Definition: atmdat_3body.cpp:28
t_thermal::lgTemperatureConstant
bool lgTemperatureConstant
Definition: thermal.h:32
MIN2
#define MIN2
Definition: cddefines.h:761
SetDeuteriumIonization
void SetDeuteriumIonization(const double &xNeutral, const double &xIonized)
Definition: deuterium.cpp:26
t_secondaries::sec2total
realnum sec2total
Definition: secondaries.h:27
t_save::lgTraceConvergeBaseHash
bool lgTraceConvergeBaseHash
Definition: save.h:341
nzone
long int nzone
Definition: cddefines.cpp:14
GrainVar::lgGrainElectrons
bool lgGrainElectrons
Definition: grainvar.h:494
highen
void highen(void)
Definition: highen.cpp:21
ion_trim
void ion_trim(long int nelem)
Definition: ion_trim.cpp:21
t_thermal::char_tran_cool
double char_tran_cool
Definition: thermal.h:146
t_conv::lgUpdateCouplings
bool lgUpdateCouplings
Definition: conv.h:259
x0
static double x0[83]
Definition: atmdat_3body.cpp:23
iso_update_rates
void iso_update_rates(void)
Definition: iso_solve.cpp:51
mole_update_sources
void mole_update_sources(void)
Definition: mole_drive.cpp:106
t_opac::lgRedoStatic
bool lgRedoStatic
Definition: opacity.h:147
t_conv::lgConvIoniz
bool lgConvIoniz() const
Definition: conv.h:115
mole_drive
void mole_drive(void)
Definition: mole_drive.cpp:41
dense.h
t_ionbal::elecsrc
double elecsrc[LIMELM]
Definition: ionbal.h:251
RT_OTS_Update
void RT_OTS_Update(double *SumOTS)
Definition: rt_ots.cpp:488
t_conv::resetConvIoniz
void resetConvIoniz()
Definition: conv.h:100
mole
t_mole_local mole
Definition: mole.cpp:7
cooling.h
trace
t_trace trace
Definition: trace.cpp:5
cddefines.h
lgStatesConserved
void lgStatesConserved(long nelem, long ionStage, qList states, long numStates, realnum err_tol, long loop_ion)
Definition: dense.cpp:160
UpdateUTAs
void UpdateUTAs(void)
Definition: conv_base.cpp:1086
GrainVar::GrainChTrRate
realnum GrainChTrRate[LIMELM][LIMELM+1][LIMELM+1]
Definition: grainvar.h:541
thermal
t_thermal thermal
Definition: thermal.cpp:5
t_trace::nTrConvg
int nTrConvg
Definition: trace.h:27
t_rfield::lgIonizReevaluate
bool lgIonizReevaluate
Definition: rfield.h:128
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
t_save::chHashString
char chHashString[INPUT_LINE_LENGTH]
Definition: save.h:295
t_conv::lgIonStageTrimed
bool lgIonStageTrimed
Definition: conv.h:189
CONV_BASE_LOOPS
@ CONV_BASE_LOOPS
Definition: conv.h:76
grains.h
deut
t_deuterium deut
Definition: deuterium.cpp:8
t_opac::opacity_abs
double * opacity_abs
Definition: opacity.h:95
hmi.h
t_thermal::te_update
double te_update
Definition: thermal.h:128
MAX2
#define MAX2
Definition: cddefines.h:782
ionbal
t_ionbal ionbal
Definition: ionbal.cpp:5
pressure.h
iso_renorm
void iso_renorm(long nelem, long ipISO, double &renorm)
Definition: iso_solve.cpp:272
RT_OTS_PrtRate
void RT_OTS_PrtRate(double weak, int chFlag)
Definition: rt_ots.cpp:712
LIMELM
const int LIMELM
Definition: cddefines.h:258
t_conv::limPres2Ioniz
long int limPres2Ioniz
Definition: conv.h:161
t_conv::nPres2Ioniz
long int nPres2Ioniz
Definition: conv.h:152
t_dense::IonLow
long int IonLow[LIMELM+1]
Definition: dense.h:119
t_iso_sp::st
qList st
Definition: iso.h:453
UTALines
TransitionList UTALines("UTALines", &AnonStates)
t_deuterium::xIonDense
double xIonDense[2]
Definition: deuterium.h:21
save.h
x2
static double x2[63]
Definition: atmdat_3body.cpp:20
t_dense::IonHigh
long int IonHigh[LIMELM+1]
Definition: dense.h:120
iteration
long int iteration
Definition: cddefines.cpp:16
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
TauLines
TransitionList TauLines("TauLines", &AnonStates)
DynaIonize
void DynaIonize(void)
Definition: dynamics.cpp:186
t_thermal::htot
double htot
Definition: thermal.h:149
INPUT_LINE_LENGTH
const int INPUT_LINE_LENGTH
Definition: cddefines.h:254
grainvar.h
OpacityAddTotal
void OpacityAddTotal(void)
Definition: opacity_addtotal.cpp:27
t_secondaries::csupra
realnum ** csupra
Definition: secondaries.h:21
rt.h
t_conv::nTotalIoniz
long int nTotalIoniz
Definition: conv.h:166
ionbal.h
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
GrainDrive
void GrainDrive(void)
Definition: grains.cpp:1591
hmi
t_hmi hmi
Definition: hmi.cpp:5
secondaries
t_secondaries secondaries
Definition: secondaries.cpp:5
dynamics
t_dynamics dynamics
Definition: dynamics.cpp:44
t_dense::xNucleiTotal
realnum xNucleiTotal
Definition: dense.h:104
t_ionbal::lgInnerShellLine_on
bool lgInnerShellLine_on
Definition: ionbal.h:175
conv
t_conv conv
Definition: conv.cpp:5
CoolEvaluate
void CoolEvaluate(double *tot)
Definition: cool_eval.cpp:45
gv
GrainVar gv
Definition: grainvar.cpp:5
findspecies
molecule * findspecies(const char buf[])
Definition: mole_species.cpp:814
t_ionbal::RateIonizTot
double RateIonizTot(long nelem, long ion)
Definition: ionbal.h:254
t_conv::incrementCounter
void incrementCounter(const counter_type type)
Definition: conv.h:308
t_StopCalc::nTotalIonizStop
long int nTotalIonizStop
Definition: stopcalc.h:127
RT_OTS
void RT_OTS(void)
Definition: rt_ots.cpp:39
fixit
void fixit(void)
Definition: service.cpp:991
CONV_BASE_ACCELS
@ CONV_BASE_ACCELS
Definition: conv.h:77
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
taulines.h
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:398
t_dense::lgSetIoniz
bool lgSetIoniz[LIMELM]
Definition: dense.h:149
lgElemsConserved
bool lgElemsConserved(void)
Definition: dense.cpp:99
t_conv::HeatCoolRelErrorAllowed
realnum HeatCoolRelErrorAllowed
Definition: conv.h:278
eden_sum
int eden_sum(void)
Definition: eden_sum.cpp:18
phycon.h
atmdat
t_atmdat atmdat
Definition: atmdat.cpp:6
t_dense::SetIoniz
realnum SetIoniz[LIMELM][LIMELM+1]
Definition: dense.h:154
opacity.h
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
t_atmdat::CharExcRecTo
double CharExcRecTo[NCX][LIMELM][LIMELM+1]
Definition: atmdat.h:153
h2.h
t_save::ipTraceConvergeBase
FILE * ipTraceConvergeBase
Definition: save.h:342
fnzone
double fnzone
Definition: cddefines.cpp:15
molezone::den
double den
Definition: mole.h:358
ChargTranEval
void ChargTranEval(void)
Definition: atmdat_char_tran.cpp:44
t_conv::lgSearch
bool lgSearch
Definition: conv.h:175
stopcalc.h
t_phycon::te
double te
Definition: phycon.h:11
NISO
const int NISO
Definition: cddefines.h:261
ion_recom_calculate
void ion_recom_calculate(void)
Definition: ion_recomb_Badnell.cpp:1202
GrainVar::lgDustOn
bool lgDustOn() const
Definition: grainvar.h:471
mole_global
t_mole_global mole_global
Definition: mole.cpp:6
t_ionbal::elecsnk
double elecsnk[LIMELM]
Definition: ionbal.h:251
t_mole_global::num_calc
int num_calc
Definition: mole.h:314
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
save
t_save save
Definition: save.cpp:5
ipNICKEL
const int ipNICKEL
Definition: cddefines.h:332
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
deuterium.h
HeatZero
void HeatZero(void)
Definition: heat_sum.cpp:928
TransitionList::Emis
EmissionList & Emis()
Definition: transition.h:329
Generated by   doxygen 1.8.17