cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
radius_first.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 /*radius_first derive thickness of first zone, called after conditions in first zone
4  * are established, sets
5 radius.drad_x_fillfac
6 radius.drad
7  */
8 #include "cddefines.h"
9 #define Z 1.0001
10 #include "wind.h"
11 #include "stopcalc.h"
12 #include "thermal.h"
13 #include "dynamics.h"
14 #include "trace.h"
15 #include "punch.h"
16 #include "pressure.h"
17 #include "iso.h"
18 #include "h2.h"
19 #include "rfield.h"
20 #include "dense.h"
21 #include "hmi.h"
22 #include "geometry.h"
23 #include "opacity.h"
24 #include "ipoint.h"
25 #include "radius.h"
26 
27 void radius_first(void)
28 {
29  long int i ,
30  ip;
31 
32  bool lgDoPun;
33 
34  int indexOfSmallest = 0;
35 
36  const int NUM_DR_TYPES = 13;
37 
38  struct t_drValues{
39  double dr;
40  char whatToSay[40];
41  } drValues[NUM_DR_TYPES];
42 
43  double accel,
44  BigOpacity,
45  change,
46  dr912,
47  drH2 ,
48  drContPres ,
49  drOpacity ,
50  drStromgren, /* used for Stromgren length */
51  drTabDen,
52  dradf,
53  drcol,
54  dr_time_dep,
55  drthrm,
56  factor,
57  winddr;
58  static double drad_last_iteration=-1.;
59 
60  DEBUG_ENTRY( "radius_first()" );
61 
62  /***********************************************************************
63  *
64  * wind model, use acceleration length
65  *
66  ***********************************************************************/
67 
68  if( wind.windv > 0. )
69  {
70  /* evaluate total pressure, although value not used (so stuffed into dr912) */
71  /* >>chng 01 nov 02, remove call, confirm vals defined with assert */
72  ASSERT( dense.pden > 0. && dense.wmole > 0. );
73  accel = 1.3e-10*dense.pden*dense.wmole;
74  winddr = POW2(wind.windv)/25./accel;
75  }
76  else
77  {
78  winddr = 1e30;
79  }
80 
81  /* key off of Lyman continuum optical depth */
82  if( StopCalc.taunu > 0.99 && StopCalc.taunu < 3. )
83  {
84  dr912 = StopCalc.tauend/6.3e-18/(dense.xIonDense[ipHYDROGEN][0]*geometry.FillFac)*Z/50.;
85  }
86  else
87  {
88  dr912 = 1e30;
89  }
90 
91  if( dynamics.lgStatic && iteration > 2 )
92  {
93  /* when time dependent case on do not let dr change since current continuum
94  * is not good indicator of conditions */
95  dr_time_dep = drad_last_iteration;
96  }
97  else
98  {
99  dr_time_dep = 1e30;
100  }
101 
102  /***********************************************************************
103  *
104  * key off of column density; total, neutral, or ionized
105  *
106  ***********************************************************************/
107 
108  if( StopCalc.HColStop < 5e29 )
109  {
110  /* this is useful for very thin columns, normally larger than 1st DR */
111  drcol = log10(StopCalc.HColStop) - log10(dense.gas_phase[ipHYDROGEN]*geometry.FillFac* 20.);
112  }
113  else if( StopCalc.colpls < 5e29 )
114  {
115  /* ionized column density */
116  drcol = log10(StopCalc.colpls) - log10(dense.xIonDense[ipHYDROGEN][1]*geometry.FillFac* 20.);
117  }
118  else if( StopCalc.colnut < 5e29 )
119  {
120  /* neutral column denisty */
121  drcol = log10(StopCalc.colnut) - log10(dense.xIonDense[ipHYDROGEN][0]*geometry.FillFac*50.);
122  }
123  else
124  {
125  /* not used */
126  drcol = 30.;
127  }
128  /* finally convert the drived column density to linear scale */
129  drcol = pow(10.,MIN2(35.,drcol));
130 
131  /***********************************************************************
132  *
133  * key off of density or abundance fluctuations, must be small part of wavelength
134  *
135  ***********************************************************************/
136 
137  if( dense.flong != 0. )
138  {
139  /* flong set => density fluctuations */
140  dradf = 6.283/dense.flong/10.;
141  dradf = MIN4(dradf,radius.router[iteration-1]*Z,drcol,dr912);
142  }
143  else
144  {
145  dradf = FLT_MAX;
146  }
147 
148  /* >>>chng 99 nov 18, add check on stromgren length */
149  /* estimate Stromgren length, but only if there are ionizing photons
150  * and not constant temperature model */
151  if( (rfield.qhtot>0.) && (rfield.qhtot> rfield.qbal*0.01) && (rfield.uh>1e-10) )
152  {
153  /* >>chng 99 dec 23, double to allow lte.in to work on alphas */
154  /* >>chng 03 mar 15, density to double to avoid overflow, PvH */
155  drStromgren = (double)(rfield.qhtot)/iso.RadRec_caseB[ipH_LIKE][ipHYDROGEN]/
156  POW2((double)dense.gas_phase[ipHYDROGEN]);
157 
158  /* different logic if this is a sphere */
159  if( drStromgren/radius.rinner > 1. )
160  {
161  /* >>chng 03 mar 15, to double to avoid FP overflow, PvH */
162  drStromgren = (double)rfield.qhtot*3./(radius.rinner*
163  iso.RadRec_caseB[ipH_LIKE][ipHYDROGEN]*POW2((double)dense.gas_phase[ipHYDROGEN]) );
164  drStromgren += 1.;
165  /* this results in r_out / r_in */
166  drStromgren = pow( drStromgren , 0.33333);
167  /* make it a physics thickness in cm */
168  drStromgren *= radius.rinner;
169  }
170 
171  /* remember the Stromgren thickness */
172  radius.thickness_stromgren = (realnum)drStromgren;
173 
174  /* take one hundredth of this */
175  drStromgren /= 100.;
176  }
177  else
178  {
179  drStromgren = FLT_MAX;
180  radius.thickness_stromgren = FLT_MAX;
181  }
182 
183  /***********************************************************************
184  *
185  * find largest opacity, to keep the first zone optical depth 1
186  * this is usually the physics that sets the first zone thickness
187  *
188  ***********************************************************************/
189 
190  /* >>>chng 99 jun 25, this is to simulate behavior of code before extension
191  * of continuum array to 1e-8 Ryd */
192  ip = ipoint(1e-5);
193 
194  /* find largest opacity */
195  BigOpacity = 0.;
196  for( i=ip; i < rfield.nflux; i++ )
197  {
198  /* remember largest opacity, and energy where this happened,
199  * make sure flux at energy is gt 0, can be zero for case where
200  * nflux increased to include emission from some ions */
201  if( rfield.flux[i]>0. && opac.opacity_abs[i] > BigOpacity )
202  {
203  BigOpacity = opac.opacity_abs[i];
204  }
205  }
206  /* BigOpacity may be zero on very first call */
207 
208  /* drChange set with set didz command, is only number set with this command,
209  * default in zerologic is 0.15
210  * set drad to small part of*/
211  if( BigOpacity > SMALLFLOAT )
212  {
213  drOpacity = (radius.drChange/100.)/BigOpacity/geometry.FillFac;
214  }
215  else
216  {
217  drOpacity = 1e30;
218  }
219 
220  /***********************************************************************
221  *
222  * thermalization length of typical lines
223  *
224  ***********************************************************************/
225 
226  drthrm = 1.5e31/MAX2(1.,POW2((double)dense.gas_phase[ipHYDROGEN]));
227 
228  /***********************************************************************
229  *
230  * make sure we resolve initial structure in dense_tabden command
231  * if interpolated table we need to make sure that we resolve the
232  * initial changes in the structure
233  *
234  ***********************************************************************/
235 
236  if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
237  {
238  drTabDen = 1.;
239  i = 1;
240  factor = 0.;
241  while( i < 100 && factor < 0.05 )
242  {
243  /* check densities at ever larger dr's, until factor becomes more than 5% */
244  factor = dense.gas_phase[ipHYDROGEN]/
245  dense_tabden(radius.Radius+drTabDen, drTabDen );
246  /* density change can be positive or negative sign */
247  factor = fabs(factor-1.);
248  drTabDen *= 2.;
249  i += 1;
250  }
251  drTabDen /= 2.;
252  }
253  else
254  {
255  drTabDen = 1e30;
256  }
257 
258  /* >>chng 03 mar 20, add check on lyman band optical depth - want first zone
259  * to be thin in H2 bands */
260  /* some tests are fully molecular with solomon process turned off,
261  * do not sense this when already almost fully molecular */
262  if( hmi.H2_total/dense.gas_phase[ipHYDROGEN] < 0.1 )
263  {
264  change = 0.1;
265  }
266  else
267  {
268  /* >>chng 04 mar 14, this branch, H is quite molecular,
269  * still do not want large changes in solomon rate since linearization
270  * would not work in hmole network, bu do not need such fine steps */
271  change = 1.;
272  }
273 
274  /* >>chng 04 mar 14 go back to original logic since molecular
275  * pdr's had big jump in conditions from
276  * first to second zon even when most H in H2
277  change = 0.1; */
278  /* >>chng 04 apr 18, change from 0.1 to 0.001, inital zones too large in
279  * leiden test case f1 */
280  change = 0.001;
281  /* >>chng 04 mar 13, not too large when big H2 is on */
282  if( h2.lgH2ON && hmi.lgBigH2_evaluated )
283  {
284  if( fabs(hmi.HeatH2Dexc_BigH2)/thermal.ctot > 0.05 )
285  {
286  /* changes in H2 heating caused by changes in solomon rate
287  * would drive temperature failures */
288  /* >>chng 04 apr 18, change from 0.001 to 0.0001, inital zones too large in
289  * leiden test case f1 */
290  change = 0.0001;
291  }
292  else
293  {
294  /* >>chng 04 apr 18, change from 0.01 to 0.001, inital zones too large in
295  * leiden test case f1 */
296  change = 0.001;
297  }
298  }
299  drH2 = change / SDIV(
301 
302  /* >>chng 06 feb 01, very high U ulirg models had dramatic increase in
303  * cont pre in first few zones,
304  * in constant total pressure case, don't want acceleration across first zone to
305  * be large compared with current gas pressure */
306  if( (strcmp( dense.chDenseLaw, "CPRE" )==0) && pressure.lgContRadPresOn )
307  {
308  /* radiative acceleration was evaluated in PressureTotal */
309  drContPres = 0.05 * pressure.PresTotlCurr /
311  }
312  else if( wind.windv != 0. )
313  {
314  /* acceleration and change in v in wind */
315  double g = fabs(wind.AccelTot-wind.AccelGravity);
316  /* wind - do not let velocity change by too much */
317  drContPres = 0.05*POW2(wind.windv)/(2.*SDIV(g));
318  }
319  else
320  drContPres = 1e30;
321 
322  drValues[0].dr = drOpacity;
323  drValues[1].dr = radius.Radius/20.;
324  drValues[2].dr = drStromgren;
325  drValues[3].dr = radius.router[iteration-1]/10.;
326  drValues[4].dr = drcol;
327  drValues[5].dr = dr912;
328  drValues[6].dr = drthrm;
329  drValues[7].dr = winddr;
330  drValues[8].dr = dradf;
331  drValues[9].dr = drTabDen;
332  drValues[10].dr = drH2;
333  drValues[11].dr = drContPres;
334  drValues[12].dr = dr_time_dep;
335 
336  strcpy( drValues[0].whatToSay, "drOpacity" );
337  strcpy( drValues[1].whatToSay, "radius.Radius/20.");
338  strcpy( drValues[2].whatToSay, "drStromgren");
339  strcpy( drValues[3].whatToSay, "radius.router[iteration-1]/10.");
340  strcpy( drValues[4].whatToSay, "drcol");
341  strcpy( drValues[5].whatToSay, "dr912");
342  strcpy( drValues[6].whatToSay, "drthrm");
343  strcpy( drValues[7].whatToSay, "winddr");
344  strcpy( drValues[8].whatToSay, "dradf");
345  strcpy( drValues[9].whatToSay, "drTabDen");
346  strcpy( drValues[10].whatToSay, "drH2");
347  strcpy( drValues[11].whatToSay, "drContPres");
348  strcpy( drValues[12].whatToSay, "dr_time_dep");
349 
350  for( i=0; i<NUM_DR_TYPES; i++ )
351  {
352  if( drValues[i].dr < drValues[indexOfSmallest].dr )
353  {
354  indexOfSmallest = i;
355  }
356  }
357 
358  radius.drad = drValues[indexOfSmallest].dr;
359 
360  /* reset if radius.drad is less than radius.sdrmin */
361  if( radius.sdrmin >= radius.drad )
362  {
364  /* set flag for comment if the previous line forced a larger dr than
365  * would otherwise have been chosen. will cause comment to be generated
366  * in PrtComment if set true*/
367  radius.lgDR2Big = true;
368  }
369  else
370  {
371  radius.lgDR2Big = false;
372  }
373 
374  /* this min had been in the big min set above, but caused a false alarm
375  * on the lgDR2Big test above since the set dr command sets both in and max */
377 
378 #if 0
379  /***********************************************************************
380  *
381  * we have now generated range of estimates of first thickness,
382  * now choose smallest of the group
383  *
384  ***********************************************************************/
385 
386  /* radius div by 20, to prevent big change in iron ionization for high ioniz gas,
387  * this is also the ONLY place that sphericity comes in */
388  radius.drad = MIN4( MIN3( drOpacity, radius.Radius/20., drStromgren ),
389  MIN3( radius.router[iteration-1]/10., drcol, dr912 ),
390  MIN4( drthrm, winddr, dradf, drTabDen ),
391  MIN3( drH2, drContPres, dr_time_dep ) );
392 
393  /* option to set lower limit to zone thickness, with set drmin command*/
395 
396  /* set flag for comment if the previous line forced a larger dr than
397  * would otherwise have been chosen. will cause comment to be generated
398  * in PrtComment if set true*/
399  if( fp_equal( radius.drad, radius.sdrmin ) )
400  {
401  radius.lgDR2Big = true;
402  }
403  else
404  {
405  radius.lgDR2Big = false;
406  }
407 
408  /* this min had been in the big min set above, but caused a false alarm
409  * on the lgDR2Big test above since the set dr command sets both in and max */
411 #endif
412 
414 
415  /* save dr for this iteration */
416  drad_last_iteration = radius.drad;
417 
418  /* drMinimum is smallest acceptable DRAD, and is 1/100 OF DRAD(1) */
419  /* this can be turned off by GLOB command */
420  if( radius.lgDrMnOn )
421  {
422  /* >>chng 05 mar 05, drMinimum is now drad * hden, to make propro to optical depth
423  * avoid false trigger across thermal fronts
424  * add * dense.gas_phase */
425  /* NB - drMinimum not used in code - delete? */
426  radius.drMinimum = (realnum)(radius.drad * dense.gas_phase[ipHYDROGEN]/1e7);
427  }
428  else
429  {
430  radius.drMinimum = 0.;
431  }
432 
433  /* if set drmin is used, make sure drMinimum (which will cause an abort) is
434  * smaller than drmin */
435  if( radius.lgSMinON )
436  {
437  /* >>chng 05 mar 05, drMinimum is now drad * hden, to make propro to optical depth
438  * avoid false trigger across thermal fronts
439  * add * dense.gas_phase */
440  /* NB - drMinimum not used in code - delete? */
442  }
443 
444  if( trace.lgTrace )
445  {
446  fprintf( ioQQQ,
447  " radius_first called, finds dr=%13.5e drMinimum=%12.3e sdrmin=%10.2e sdrmax=%10.2e\n",
449  }
450 
451  if( radius.drad < SMALLFLOAT*1.1 )
452  {
453  fprintf( ioQQQ,
454  " PROBLEM radius_first detected likely insanity, found dr=%13.5e \n", radius.drad);
455  fprintf( ioQQQ,
456  " radius_first: calculation continuing but crash is likely. \n");
457  /* this sets flag that insanity has occurred */
458  TotalInsanity();
459  }
460 
461  /* all this is to only punch on last iteration
462  * the punch dr command is not really a punch command, making this necessary
463  * lgDRon is set true if "punch dr" entered */
464  if( punch.lgDROn )
465  {
466  lgDoPun = true;
467  }
468  else
469  {
470  lgDoPun = false;
471  }
472 
473  /* punch what we decided up? */
474  if( lgDoPun )
475  {
476  /* create hash marks on second and later iterations */
477  if( iteration > 1 && punch.lgDRHash )
478  {
479  static int iter_punch=-1;
480  if( iteration !=iter_punch )
481  fprintf( punch.ipDRout, "%s\n",punch.chHashString );
482  iter_punch = iteration;
483  }
484  /* this is common part of each line, the zone count, depth, chosen dr, and depth2go */
485  /* >>chng 05 aug 15, had printed drNext, always zero, rather the drad, which is set here */
486  fprintf( punch.ipDRout , "%ld\t%.5e\t%.3e\t%.3e\t", nzone, radius.depth, radius.drad, radius.Depth2Go );
487 
488  if( radius.lgDR2Big )
489  {
490  fprintf( punch.ipDRout,
491  "radius_first keys from radius.sdrmin\n");
492 
493  }
494  else if( fp_equal( radius.drad, radius.sdrmax ) )
495  {
496 
497  fprintf( punch.ipDRout,
498  "radius_first keys from radius.sdrmax\n");
499  }
500  else
501  {
502  ASSERT( indexOfSmallest < NUM_DR_TYPES - 1 );
503  fprintf( punch.ipDRout, "radius_first keys from %s\n",
504  drValues[indexOfSmallest].whatToSay);
505  }
506 
507  /* \todo 1 improve this printout and the drValues treatment above. */
508 
509 #if 0
510  if( fp_equal( radius.drad, drOpacity ) )
511  {
512  fprintf( punch.ipDRout,
513  "radius_first keys from drOpacity, opac was %.2e at %.2e Ryd\n",
514  BigOpacity , BigOpacityAnu );
515  }
516  else if( fp_equal( radius.drad, radius.Radius/20. ) )
517  {
518  fprintf( punch.ipDRout,
519  "radius_first keys from radius.Radius\n" );
520  }
521  else if( fp_equal( radius.drad, drStromgren ) )
522  {
523  fprintf( punch.ipDRout,
524  "radius_first keys from drStromgren\n");
525  }
526  else if( fp_equal( radius.drad, dr_time_dep ) )
527  {
528  fprintf( punch.ipDRout,
529  "radius_first keys from time dependent\n");
530  }
531  else if( fp_equal( radius.drad, radius.router[iteration-1]/10. ) )
532  {
533  fprintf( punch.ipDRout,
534  "radius_first keys from radius.router[iteration-1]\n");
535  }
536  else if( fp_equal( radius.drad, drcol ) )
537  {
538  fprintf( punch.ipDRout,
539  "radius_first keys from drcol\n");
540  }
541  else if( fp_equal( radius.drad, radius.sdrmin ) )
542  {
543  fprintf( punch.ipDRout,
544  "radius_first keys from radius.sdrmin\n");
545  }
546  else if( fp_equal( radius.drad, dr912 ) )
547  {
548  fprintf( punch.ipDRout,
549  "radius_first keys from dr912\n");
550  }
551  else if( fp_equal( radius.drad, radius.sdrmax ) )
552  {
553  fprintf( punch.ipDRout,
554  "radius_first keys from radius.sdrmax\n");
555  }
556  else if( fp_equal( radius.drad, drthrm ) )
557  {
558  fprintf( punch.ipDRout,
559  "radius_first keys from drthrm\n");
560  }
561  else if( fp_equal( radius.drad, winddr ) )
562  {
563  fprintf( punch.ipDRout,
564  "radius_first keys from winddr\n");
565  }
566  else if( fp_equal( radius.drad, drH2 ) )
567  {
568  fprintf( punch.ipDRout,
569  "radius_first keys from H2 lyman lines\n");
570  }
571  else if( fp_equal( radius.drad, dradf ) )
572  {
573  fprintf( punch.ipDRout,
574  "radius_first keys from dradf\n");
575  }
576  else if( fp_equal( radius.drad, drTabDen ) )
577  {
578  fprintf( punch.ipDRout,
579  "radius_first keys from drTabDen\n");
580  }
581  else if( fp_equal( radius.drad, drContPres ) )
582  {
583  fprintf( punch.ipDRout,
584  "radius_first keys from radiative acceleration across zone\n");
585  }
586  else
587  {
588  fprintf( punch.ipDRout, "radius_first insanity\n" );
589  fprintf( ioQQQ, "radius_first insanity, radius is %e\n" ,
590  radius.drad);
591  ShowMe();
592  }
593 #endif
594 
595  }
596  return;
597 }

Generated for cloudy by doxygen 1.8.1.1