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