cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
conv_init_solution.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 /*ConvInitSolution drive search for initial temperature, for illuminated face */
4 /*lgTempTooHigh returns true if temperature is too high */
5 /*lgCoolHeatCheckConverge return true if converged, false if change in heating cooling larger than set tolerance */
6 #include "cddefines.h"
7 #include "physconst.h"
8 #include "trace.h"
9 #include "struc.h"
10 #include "rfield.h"
11 #include "mole.h"
12 #include "dense.h"
13 #include "stopcalc.h"
14 #include "heavy.h"
15 #include "wind.h"
16 #include "geometry.h"
17 #include "thermal.h"
18 #include "radius.h"
19 #include "phycon.h"
20 #include "pressure.h"
21 #include "conv.h"
22 
23 /* derivative of net cooling wrt temperature to check on sign oscillations */
24 static double dCoolNetDTOld = 0;
25 
26 static double OxyInGrains , FracMoleMax;
27 /* determine whether chemistry is important - what is the largest
28  * fraction of an atom in molecules? Also is ice formation
29  * important
30  * sets above two file static variables */
31 
32 /*lgCoolHeatCheckConverge return true if converged, false if change in heating cooling larger than set tolerance */
34  double *CoolNet )
35 {
36  static double HeatOld=-1. , CoolOld=-1.;
37  DEBUG_ENTRY( "lgCoolHeatCheckConverge()" );
38 
39  double HeatChange = thermal.htot - HeatOld;
40  double CoolChange = thermal.ctot - CoolOld;
41  /* will use larger of heating or cooling in tests for relative convergence
42  * because in constant temperature models one may have trivially small value */
43  double HeatCoolMax = MAX2( thermal.htot , thermal.ctot );
44 
45  /* is the heating / cooling converged?
46  * get max relative change, determines whether converged heating/cooling */
47  double HeatRelChange = fabs(HeatChange)/SDIV( HeatCoolMax );
48  double CoolRelChange = fabs(CoolChange)/SDIV( HeatCoolMax );
49  bool lgConverged = true;
50  if( MAX2( HeatRelChange , CoolRelChange ) > conv.HeatCoolRelErrorAllowed )
51  lgConverged = false;
52 
53  *CoolNet = thermal.ctot - thermal.htot;
54 
55  HeatOld = thermal.htot;
56  CoolOld = thermal.ctot;
57  return lgConverged;
58 }
59 
60 /* call lgCoolHeatCheckConverge until cooling and heating are converged */
62  double *CoolNet ,
63  double *dCoolNetDT )
64 {
65  const long int LOOP_MAX=20;
66  long int loop = 0;
67  bool lgConvergeCoolHeat = false;
68 
69  DEBUG_ENTRY( "lgCoolNetConverge()" );
70 
71  while( loop < LOOP_MAX && !lgConvergeCoolHeat && !lgAbort )
72  {
73  if( ConvEdenIoniz() )
74  {
75  /* this is an error return, calculation will immediately stop */
76  lgAbort = true;
77  }
78  lgConvergeCoolHeat = lgCoolHeatCheckConverge( CoolNet );
79  if( trace.lgTrace || trace.nTrConvg>=3 )
80  fprintf(ioQQQ," lgCoolNetConverge %li calls to ConvEdenIoniz, converged? %c\n",
81  loop , TorF( lgConvergeCoolHeat ) );
82  ++loop;
83  }
84 
85  if( thermal.ConstTemp > 0 )
86  {
87  /* constant temperature model - this is trick so that temperature
88  * updater uses derivative to go to the set constant temperature */
89  *CoolNet = phycon.te - thermal.ConstTemp;
90  *dCoolNetDT = 1.f;
91  }
92  else if( phycon.te < 4000.f )
93  {
94  /* at low temperatures the analytical cooling-heating derivative
95  * is often of no value - the species densities change when the
96  * temperature changes. Use simple dCnet/dT=(C-H)/T - this is
97  * usually too large and results is smaller than optical dT, but
98  * we do eventually converge */
99  *dCoolNetDT = thermal.ctot / phycon.te;
100  }
101  else if( thermal.htot / thermal.ctot > 2.f )
102  {
103  /* we are far from the solution and the temperature is much lower than
104  * equilibrium - analytical derivative from cooling evaluation is probably
105  * wrong since coolants are not correct */
106  *dCoolNetDT = thermal.ctot / phycon.te;
107  }
108  else
109  {
110  *dCoolNetDT = thermal.dCooldT - thermal.dHeatdT;
111  if( dCoolNetDTOld * *dCoolNetDT < 0. )
112  {
113  /* derivative is oscillating */
114  *dCoolNetDT = *CoolNet / phycon.te;
115  fprintf(ioQQQ,"PROBLEM CoolNet derivative oscillating in initial solution \n");
116  }
117  }
118  return lgConvergeCoolHeat;
119 }
120 
121 /*ChemImportance find fraction of heavy elements in molecules, O in ices */
122 STATIC void ChemImportance( void )
123 {
124  long int i;
125 
126  DEBUG_ENTRY( "ChemImportance()" );
127 
128  FracMoleMax = 0.;
129  long int ipMax = -1;
130  for( i=0; i<mole.num_comole_calc; ++i )
131  {
132  /* is element included in the chemistry, and is it a molecule? */
133  if( dense.lgElmtOn[COmole[i]->nelem_hevmol] && COmole[i]->n_nuclei>1 )
134  {
136  frac *= (realnum) COmole[i]->nElem[COmole[i]->nelem_hevmol];
137  if( frac > FracMoleMax )
138  {
139  /* remember largest fraction */
140  FracMoleMax = frac;
141  ipMax = i;
142  }
143  }
144  }
145 
146  /* fraction of all oxygen in ices on grains */
147  OxyInGrains = 0;
148  for(i=0;i<mole.num_comole_calc;++i)
149  {
150  OxyInGrains += (1 - COmole[i]->lgGas_Phase)*COmole[i]->hevmol*COmole[i]->nElem[ipOXYGEN];
151  }
152  /* this is now fraction of O in ices */
154 
155  {
156  /* option to print out entire matrix */
157  enum {DEBUG_LOC=false};
158  if( DEBUG_LOC )
159  {
160  fprintf(ioQQQ,"DEBUG fraction molecular %.2e species %li O ices %.2e\n" ,
161  FracMoleMax , ipMax ,OxyInGrains );
162  }
163  }
164 
165  return;
166 }
167 
169 {
170  double TeChangeFactor;
171 
172  DEBUG_ENTRY( "FindTempChangeFactor()" );
173 
174  /* find fraction of atoms in molecules - need small changes
175  * in temperature if fully molecular since chemistry
176  * network is complex and large changes in solution would
177  * cause linearization to fail */
178  /* this evaluates the global variables FracMoleMax and
179  * OxyInGrains */
180  ChemImportance();
181 
182  /* Following are series of chemical
183  * tests that determine the factor by which
184  * which can change the temperature. must be VERY small when
185  * ice formation on grains is important due to exponential sublimation
186  * rate. Also small when molecules are important, to keep chemistry
187  * within limits of linearized solver
188  * the complete linearization that is implicit in all these solvers
189  * will not work when large Te jumps occur when molecules/ices
190  * are important - solvers can't return to solution if too far away
191  * keep temp steps small when mole/ice is important */
192  if( OxyInGrains > 0.99 )
193  {
194  TeChangeFactor = 0.999;
195  }
196  else if( OxyInGrains > 0.9 )
197  {
198  TeChangeFactor = 0.99;
199  }
200  /* >>chng 97 mar 3, make TeChangeFactor small when close to lower bound */
201  else if( phycon.te < 5. ||
202  /*>>chng 06 jul 30, or if molecules/ices are important */
203  OxyInGrains > 0.1 )
204  {
205  TeChangeFactor = 0.97;
206  }
207  /*>>chng 07 feb 23, add test of chemistry being very important */
208  else if( phycon.te < 20. || FracMoleMax > 0.1 )
209  {
210  /* >>chng 07 feb 24, changed from 0,9 to 0.95 to converge
211  * pdr_coolbb.in test */
212  TeChangeFactor = 0.95;
213  }
214  else
215  {
216  /* this is the default largest step */
217  TeChangeFactor = 0.8;
218  }
219  return TeChangeFactor;
220 }
221 
222 /*ConvInitSolution drive search for initial temperature, for illuminated face */
224 {
225  long int i,
226  ionlup,
227  nelem ,
228  nflux_old,
229  nelem_reflux ,
230  ion_reflux;
231  /* current value of Cooling - Heating */
232  bool lgConvergeCoolHeat;
233 
234  double
235  TeChangeFactor,
236  TeBoundLow,
237  TeBoundHigh;
238 
239  DEBUG_ENTRY( "ConvInitSolution()" );
240 
241  /* set NaN for safety */
242  set_NaN( TeBoundLow );
243  set_NaN( TeBoundHigh );
244 
245  /* this counts number of times ConvBase is called by PressureChange, in current zone */
246  conv.nPres2Ioniz = 0;
247  /* this counts how many times ConvBase is called in this iteration
248  * and is flag used by various routines to understand they are
249  * being called the first time*/
250  conv.nTotalIoniz = 0;
251  /* ots rates not oscillating */
252  conv.lgOscilOTS = false;
253 
254  lgAbort = false;
255  dense.lgEdenBad = false;
256  dense.nzEdenBad = 0;
257  /* these are variables to remember the biggest error in the
258  * electron density, and the zone in which it occurred */
259  conv.BigEdenError = 0.;
260  conv.AverEdenError = 0.;
261  conv.BigHeatCoolError = 0.;
262  conv.AverHeatCoolError = 0.;
263  conv.BigPressError = 0.;
264  conv.AverPressError = 0.;
265 
266  /* these are equal if set dr was used, and we know the first dr */
268  {
271  }
272 
273  /* the H+ density is zero in sims with no H-ionizing radiation and no cosmic rays,
274  * the code does work in this limit. But it is unphysical for the H0 density
275  * to be zero. This could only happen in the fully molecular limit, and the
276  * code does not go to this limit at this stage, due to simple approximations. */
278 
279  if( trace.nTrConvg )
280  {
281  fprintf( ioQQQ, "\nConvInitSolution entered \n" );
282  }
283 
284  /********************************************************************
285  *
286  * this is second or higher iteration, reestablish original temperature
287  *
288  *********************************************************************/
289  if( iteration != 1 )
290  {
291  /* this is second or higher iteration on multi-iteration model */
292  if( trace.lgTrace || trace.nTrConvg )
293  {
294  fprintf( ioQQQ, " ConvInitSolution called, ITER=%2ld resetting Te to %10.2e\n",
295  iteration, struc.testr[0] );
296  }
297 
298  if( trace.lgTrace || trace.nTrConvg )
299  {
300  fprintf( ioQQQ, " search set true\n" );
301  }
302 
303  /* search phase must be turned on so that variables such as the ots rates,
304  * secondary ionization, and auger yields, can converge more quickly to
305  * proper values */
306  conv.lgSearch = true;
307 
308  /* this is the temperature and pressure from the previous iteration */
309  /*phycon.te = TempPrevIteration;*/
310  /*>>chng 06 jun 09, reset to this rather than init value set in this routine */
311  TempChange( struc.testr[0] , false);
312 
313  /* the initial pressure should now be valid */
314  /* this sets values of pressure.PresTotlCurr */
315  PresTotCurrent();
316 
317  /* now get final temperature */
318  if( ConvPresTempEdenIoniz() )
319  {
320  /* this is an error return, calculation will immediately stop */
321  lgAbort = true;
322  return(1);
323  }
324 
325  if( trace.lgTrace || trace.nTrConvg )
326  {
327  fprintf( ioQQQ, " ========================================================================================================================\n");
328  fprintf( ioQQQ, " ConvInitSolution: search set false 2 Te=%.3e========================================================================================\n" , phycon.te);
329  fprintf( ioQQQ, " ========================================================================================================================\n");
330  }
331  conv.lgSearch = false;
332 
333  /* now get final temperature */
334  if( ConvTempEdenIoniz() )
335  {
336  /* this is an error return, calculation will immediately stop */
337  lgAbort = true;
338  return(1);
339  }
340 
341  /* the initial pressure should now be valid
342  * this sets values of pressure.PresTotlCurr */
343  PresTotCurrent();
344  }
345 
346  else
347  {
348  /********************************************************************
349  *
350  * do first te from scratch
351  *
352  *********************************************************************/
353 
354  /* say that we are in search phase */
355  conv.lgSearch = true;
356 
357  if( trace.lgTrace )
358  {
359  fprintf( ioQQQ, " ConvInitSolution called, new temperature.\n" );
360  }
361 
362  /* coming up to here Te is either 4,000 K (usually) or 10^6 */
363  TeBoundLow = phycon.TEMP_LIMIT_LOW/1.001;
364 
365  /* set temperature floor option - StopCalc.TeFloor is usually
366  * zero but reset this this command - will go over to constant
367  * temperature calculation if temperature falls below floor */
368  TeBoundLow = MAX2( TeBoundLow , StopCalc.TeFloor );
369 
370  /* highest possible temperature - must not get this high since
371  * parts of code will abort if value is reached.
372  * divide by 1.2 to prevent checking on temp > 1e10 */
373  TeBoundHigh = phycon.TEMP_LIMIT_HIGH/1.2;
374 
375  /* set initial temperature, options are constant temperature,
376  * or approach equilibrium from either high or low TE */
377  double TeFirst;
378  if( thermal.ConstTemp > 0 )
379  {
380  /* constant temperature , set 4000 K then relax to desired value
381  * for cold temperatures, to slowly approach fully molecular solution.
382  * if constant temperature is higher than 4000., go right to it */
383  TeFirst = MAX2( 4000. , thermal.ConstTemp );
384  /* true allow ionization stage trimming, false blocks it */
385  }
386  else if( thermal.lgTeHigh )
387  {
388  /* approach from high TE */
389  TeFirst = 1e6;
390  }
391  else
392  {
393  /* approach from low TE */
394  TeFirst = 4000.;
395  }
396  TempChange(TeFirst , false);
397  if( trace.lgTrace || trace.nTrConvg>=2 )
398  fprintf(ioQQQ,"ConvInitSolution doing initial solution with temp=%.2e\n",
399  phycon.te);
400 
401  /* this sets values of pressure.PresTotlCurr */
402  PresTotCurrent();
403 
404  thermal.htot = 1.;
405  thermal.ctot = 1.;
406 
407  /* call cooling, heating, opacity, loop to convergence
408  * this is very first call to it, by default is at 4000K */
409 
410  double CoolNet=0, dCoolNetDT=0;
411  /* do ionization twice to get stable solution, evaluating initial dR each time */
412  lgConvergeCoolHeat = false;
413  for( ionlup=0; ionlup<2; ++ionlup )
414  {
415  if( trace.lgTrace || trace.nTrConvg>=2 )
416  fprintf( ioQQQ, "ConvInitSolution calling lgCoolNetConverge "
417  "initial setup %li with Te=%.3e C=%.2e H=%.2e d(C-H)/dT %.3e\n",
418  ionlup,phycon.te , thermal.ctot , thermal.htot,dCoolNetDT );
419  /* do not flag oscillating d(C-H)/dT until stable */
420  dCoolNetDTOld = 0.;
421  lgConvergeCoolHeat = lgCoolNetConverge( &CoolNet , &dCoolNetDT );
422  if( lgAbort )
423  break;
424  /* set thickness of first zone, this affects the pumping rates due
425  * to correction for attenuation across zone, so must be known
426  * for ConvEdenIoniz to get valid solution - this is why it
427  * is in this loop */
428  radius_first();
429  }
430  if( !lgConvergeCoolHeat )
431  fprintf(ioQQQ," PROBLEM ConvInitSolution did not converge the heating-cooling during initial search.\n");
432 
433  if( lgAbort )
434  {
435  /* we hit an abort */
436  fprintf( ioQQQ, " DISASTER PROBLEM ConvInitSolution: Search for "
437  "initial conditions aborted - lgAbort set true.\n" );
438  ShowMe();
439  /* this is an error return, calculation will immediately stop */
440  return(1);
441  }
442 
443  /* we now have error in C-H and its derivative - following is better
444  * derivative for global case where we may be quite far from C==H */
445  if( thermal.ConstTemp<=0 )
446  dCoolNetDT = thermal.ctot / phycon.te;
447  bool lgConvergedLoop = false;
448  long int LoopThermal = 0;
449  /* initial temperature is usually 4000K, if the dTe scale factor is 0.95
450  * it will take 140 integrations to lower temperature to 3K,
451  * and many more if ices are important. */
452  const long int LIMIT_THERMAL_LOOP=300;
453  double CoolMHeatSave=-1. , TempSave=-1. , TeNew=-1.,CoolSave=-1;
454  while( !lgConvergedLoop && LoopThermal < LIMIT_THERMAL_LOOP )
455  {
456  /* change temperature until we sign of C-H changes */
457  CoolMHeatSave = CoolNet;
458  dCoolNetDTOld = dCoolNetDT;
459  CoolSave = thermal.ctot;
460  TempSave = phycon.te;
461 
462  /* find temperature change factor, a number less than 1*/
463  TeChangeFactor = FindTempChangeFactor();
464  ASSERT( TeChangeFactor>0. && TeChangeFactor< 1. );
465 
466  TeNew = phycon.te - CoolNet / dCoolNetDT;
467 
468  TeNew = MAX2( phycon.te*TeChangeFactor , TeNew );
469  TeNew = MIN2( phycon.te/TeChangeFactor , TeNew );
470  /* change temperature */
471  TempChange(TeNew , true);
472  lgConvergeCoolHeat = lgCoolNetConverge( &CoolNet , & dCoolNetDT );
473 
474  if( trace.lgTrace || trace.nTrConvg>=2 )
475  fprintf(ioQQQ,"ConvInitSolution %4li TeNnew=%.3e "
476  "TeChangeFactor=%.3e Cnet/H=%.2e dCnet/dT=%10.2e Ne=%.3e"
477  " Ograins %.2e FracMoleMax %.2e\n",
478  LoopThermal , TeNew , TeChangeFactor ,
479  CoolNet/SDIV(thermal.htot) , dCoolNetDT,
481  if( lgAbort )
482  {
483  /* we hit an abort */
484  fprintf( ioQQQ, " Search for initial conditions aborted - "
485  "lgAbort set true, lgConvergeCoolHeat=%c.\n",
486  TorF(lgConvergeCoolHeat) );
487  /* this is an error return, calculation will immediately stop */
488  return(1);
489  }
490 
491  /* keep changing temperature until sign of CoolNet changes
492  * in constant temperature case CoolNet is delta Temperature
493  * so is zero for last pass in this loop
494  * this is for constant temperature case */
495  if( fabs(CoolNet)< SMALLDOUBLE )
496  /* CoolNet is zero */
497  lgConvergedLoop = true;
498  else if( (CoolNet/fabs(CoolNet))*(CoolMHeatSave/fabs(CoolMHeatSave)) <= 0)
499  /* change in sign of net cooling */
500  lgConvergedLoop = true;
501  else if( fabs(CoolNet)/MAX2( thermal.ctot , thermal.htot ) <conv.HeatCoolRelErrorAllowed*10. )
502  lgConvergedLoop = true;
503  else if( phycon.te <= TeBoundLow || phycon.te>=TeBoundHigh)
504  lgConvergedLoop = true;
505  ++LoopThermal;
506  }
507 
508  if( !lgConvergedLoop )
509  fprintf(ioQQQ,"PROBLEM ConvInitSolution: temperature solution not "
510  "found in initial search, final Te=%.2e\n",
511  phycon.te );
512 
513  /* interpolate temperature where C=H if not constant temperature sim
514  * will have set constant temperature mode above if we hit temperature
515  * floor */
516  if( thermal.ConstTemp<=0 &&
517  ! fp_equal( TempSave , phycon.te ) )
518  {
519  if( trace.lgTrace || trace.nTrConvg>=2 )
520  fprintf(ioQQQ," Change in sign of C-H found, Te brackets %.3e "
521  "%.3e Cool brackets %.3e %.3e ",
522  TempSave , phycon.te , CoolMHeatSave , CoolNet);
523  /* interpolate new temperature assuming Cool = a T^alpha,
524  * first find alpha */
525  double alpha = log(thermal.ctot/CoolSave) / log( phycon.te/TempSave);
526  if( fabs(alpha) < SMALLFLOAT )
527  /* alpha close to 0 means constant temperature */
528  TeNew = phycon.te;
529  else
530  {
531  /* next find log a - work in logs to avoid infinities */
532  if( thermal.ctot<0. || thermal.htot<0. )
533  TotalInsanity();
534  double Alog = log( thermal.ctot ) - alpha * log( phycon.te );
535  /* the interpolated temperature where heating equals cooling */
536  double TeLn = (log( thermal.htot ) - Alog) / alpha ;
537 
538  /* a sanity check */
539  if( TeLn < log( MIN2(phycon.te , TempSave )) )
540  TeNew = MIN2(phycon.te , TempSave );
541  else if( TeLn > log( MAX2(phycon.te , TempSave )) )
542  TeNew = MAX2(phycon.te , TempSave );
543  else
544  TeNew = pow( EE , TeLn );
545  }
546 
547  ASSERT( TeNew>=MIN2 ( TempSave , phycon.te ) ||
548  TeNew<=MAX2 ( TempSave , phycon.te ));
549 
550  if( trace.lgTrace || trace.nTrConvg>=2 )
551  fprintf(ioQQQ," interpolating to Te %.3e \n\n",
552  TeNew);
553 
554  /* change temperature */
555  TempChange(TeNew , true);
556  }
557 
558  if( ConvTempEdenIoniz() )
559  {
560  /* this is an error return, calculation will immediately stop */
561  lgAbort = true;
562  return(1);
563  }
564 
565  /* this sets values of pressure.PresTotlCurr */
566  PresTotCurrent();
567 
568  if( trace.lgTrace || trace.nTrConvg )
569  {
570  fprintf( ioQQQ, "\nConvInitSolution: end 1st iteration search phase "
571  "finding Te=%.3e\n\n" , phycon.te);
572  }
573  conv.lgSearch = false;
574 
575  if( trace.lgTrace )
576  {
577  fprintf( ioQQQ, " ConvInitSolution return, TE:%10.2e==================\n",
578  phycon.te );
579  }
580  }
581 
582  /* we now have a fairly good temperature and ionization
583  * iteration is 1 on first iteration - remember current pressure
584  * on first iteration so we can hold this constant in constant
585  * pressure simulation.
586  *
587  * flag dense.lgDenseInitConstant false if
588  * *constant pressure reset* is used -
589  * default is true, after first iteration initial density is used for
590  * first zone no matter what pressure results. Small changes occur due
591  * to radiative transfer convergence,
592  * when set false with *reset* option the density on later iterations
593  * can change to keep the pressure constant */
595  {
596  double PresNew = pressure.PresTotlCurr;
597 
598  /* option to specify the pressure as option on constant pressure command */
600  /* this is log of nT product - if not present then set zero */
602  if( trace.lgTrace )
603  {
604  fprintf( ioQQQ,
605  " PresTotCurrent 1st zn old values of PresTotlInit:%.3e"
606  " to %.3e \n",
608  PresNew);
609  }
610 
611  pressure.PresTotlInit = PresNew;
612  }
613 
614  /* now bring current pressure to the correct pressure - must do this
615  * before initial values are saved in iter start/end */
617 
618  /* this counts number of times ConvBase is called by PressureChange, in
619  * current zone these are reset here, so that we count from first zone
620  * not search */
621  conv.nPres2Ioniz = 0;
622 
623  dense.lgEdenBad = false;
624  dense.nzEdenBad = 0;
625 
626  /* save counter upon exit so we can see how many iterations were
627  * needed to do true zones */
629 
630  /* these are variables to remember the biggest error in the
631  * electron density, and the zone in which it occurred */
632  conv.BigEdenError = 0.;
633  conv.AverEdenError = 0.;
634  conv.BigHeatCoolError = 0.;
635  conv.AverHeatCoolError = 0.;
636  conv.BigPressError = 0.;
637  conv.AverPressError = 0.;
638 
639  /*>>chng 06 jun 09, only reset on first try - otherwise have stable solution? */
640  if( iteration == 1 )
641  {
642  /* now remember some things we may need even in first zone, these
643  * are normally set towards end of zone calculation in RT_tau_inc */
644  struc.testr[0] = (realnum)phycon.te;
645  /* >>chng 02 May 2001 rjrw: add hden for dilution */
647  /* pden is the total number of particles per unit vol */
649  struc.heatstr[0] = thermal.htot;
650  struc.coolstr[0] = thermal.ctot;
655  struc.ednstr[0] = (realnum)dense.eden;
658  struc.drad[0] = (realnum)radius.drad;
659  }
660 
661  /* check that nflux extends above IP of highest ionization species present.
662  * for collisional case it is possible for species to exist that are higher
663  * IP than the limit to the continuum. Need continuum to encompass all
664  * possible emission - to account for diffuse emission
665  * NB
666  * on second iteration of multi-iteration model this may result in rfield.nflux increasing
667  * which can change the solution */
668  nflux_old = rfield.nflux;
669  nelem_reflux = -1;
670  ion_reflux = -1;
671  for( nelem=2; nelem < LIMELM; nelem++ )
672  {
673  /* do not include hydrogenic species in following */
674  for( i=0; i < nelem; i++ )
675  {
676  if( dense.xIonDense[nelem][i+1] > 0. )
677  {
678  if( Heavy.ipHeavy[nelem][i] > rfield.nflux )
679  {
680  rfield.nflux = Heavy.ipHeavy[nelem][i];
681  nelem_reflux = nelem;
682  ion_reflux = i;
683  }
684  }
685  }
686  }
687 
688  /* was the upper limit to the continuum updated? if so need to define
689  * continuum variables to this new range */
690  if( nflux_old != rfield.nflux )
691  {
692  /* zero out parts of rfield arrays that were previously undefined */
693  rfield_opac_zero( nflux_old-1 , rfield.nflux );
694 
695  /* if continuum reset up, we need to define gaunt factors through high end */
696  /*tfidle(false);*/
697  /* this calls tfidle, among other things */
698  /* this sets values of pressure.PresTotlCurr */
699  PresTotCurrent();
700 
701  /* redo ionization and update opacities */
702  if( ConvBase(1) )
703  {
704  /* this is catastrophic failure */
705  lgAbort = true;
706  return( 1 );
707  }
708 
709  /* need to update continuum opacities */
710  if( trace.lgTrace )
711  {
712  fprintf(ioQQQ," nflux updated from %li to %li, anu from %g to %g \n",
713  nflux_old , rfield.nflux ,
714  rfield.anu[nflux_old] , rfield.anu[rfield.nflux] );
715  fprintf(ioQQQ," caused by element %li ion %li \n",
716  nelem_reflux ,ion_reflux );
717  }
718  }
719 
720  /* wind with negative velocity is flow from PDR - change density slightly
721  * and call pressure solver to see if it returns to original density */
722  if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv < 0. )
723  {
724  /* >>chng 04 apr 23, change pressure and let solver come back to correct
725  * pressure. This trick makes sure we are correctly either super or sub sonic
726  * since solver will jump from one branch to the other if necessary */
727 # define PCHNG 0.98
728  /* this ignores return condition, assume must be ok if got this far */
729  if( ConvPresTempEdenIoniz() )
730  {
731  /* this is an error return, calculation will immediately stop */
732  lgAbort = true;
733  return(1);
734  }
735 
737  if( ConvPresTempEdenIoniz() )
738  {
739  /* this is an error return, calculation will immediately stop */
740  lgAbort = true;
741  return(1);
742  }
743 
745  if( ConvPresTempEdenIoniz() )
746  {
747  /* this is an error return, calculation will immediately stop */
748  lgAbort = true;
749  return(1);
750  }
751 
752 # undef PCHNG
753  }
754  /* this is the only valid return and lgAbort should be false*/
755  return( lgAbort );
756 }

Generated for cloudy by doxygen 1.8.1.2