cloudy trunk
Loading...
Searching...
No Matches
radius_first.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/*radius_first derive thickness of first zone, called after conditions in first zone
4 * are established, sets
5radius.drad_x_fillfac
6radius.drad
7 */
8#include "cddefines.h"
9#include "wind.h"
10#include "stopcalc.h"
11#include "thermal.h"
12#include "dynamics.h"
13#include "trace.h"
14#include "save.h"
15#include "pressure.h"
16#include "iso.h"
17#include "h2.h"
18#include "rfield.h"
19#include "dense.h"
20#include "hmi.h"
21#include "geometry.h"
22#include "opacity.h"
23#include "ipoint.h"
24#include "radius.h"
25
26void radius_first(void)
27{
28 long int i ,
29 ip;
30
31 bool lgDoPun;
32
33 int indexOfSmallest = 0;
34
35 const double Z = 1.0001;
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.lgBallistic() )
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.lgTimeDependentStatic && 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.StopThickness[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_sp[ipH_LIKE][ipHYDROGEN].RadRec_caseB/
156 POW2((double)dense.gas_phase[ipHYDROGEN]);
157
158 // Hubble radius is largest value
159 drStromgren = MIN2(1e28 , drStromgren );
160
161 /* different logic if this is a sphere */
162 if( drStromgren/radius.rinner > 1. )
163 {
164 /* >>chng 03 mar 15, to double to avoid FP overflow, PvH */
165 drStromgren = (double)rfield.qhtot*3./(radius.rinner*
166 iso_sp[ipH_LIKE][ipHYDROGEN].RadRec_caseB*POW2((double)dense.gas_phase[ipHYDROGEN]) );
167 drStromgren += 1.;
168 /* this results in r_out / r_in */
169 drStromgren = pow( drStromgren , 0.33333);
170 /* make it a physics thickness in cm */
171 drStromgren *= radius.rinner;
172 }
173
174 /* remember the Stromgren thickness */
175 radius.thickness_stromgren = (realnum)drStromgren;
176
177 /* take one hundredth of this, do nothing if near underflow */
178 if( drStromgren > SMALLFLOAT *100.)
179 drStromgren /= 100.;
180 }
181 else
182 {
183 drStromgren = FLT_MAX;
184 radius.thickness_stromgren = FLT_MAX;
185 }
186
187 /***********************************************************************
188 *
189 * find largest opacity, to keep the first zone optical depth 1
190 * this is usually the physics that sets the first zone thickness
191 *
192 ***********************************************************************/
193
194 /* >>>chng 99 jun 25, this is to simulate behavior of code before extension
195 * of continuum array to 1e-8 Ryd */
196 ip = ipoint(1e-5);
197
198 /* find largest opacity */
199 BigOpacity = 0.;
200 for( i=ip; i < rfield.nflux; i++ )
201 {
202 /* remember largest opacity, and energy where this happened,
203 * make sure flux at energy is gt 0, can be zero for case where
204 * nflux increased to include emission from some ions */
205 if( rfield.flux[0][i]>0. && opac.opacity_abs[i] > BigOpacity )
206 {
207 BigOpacity = opac.opacity_abs[i];
208 }
209 }
210 /* BigOpacity may be zero on very first call */
211
212 /* drChange set with set didz command, is only number set with this command,
213 * default in zerologic is 0.15
214 * set drad to small part of*/
215 if( BigOpacity > SMALLFLOAT )
216 {
217 drOpacity = (radius.drChange/100.)/BigOpacity/geometry.FillFac;
218 }
219 else
220 {
221 drOpacity = 1e30;
222 }
223
224 /***********************************************************************
225 *
226 * thermalization length of typical lines
227 *
228 ***********************************************************************/
229
230 drthrm = 1.5e31/MAX2(1.,POW2((double)dense.gas_phase[ipHYDROGEN]));
231 /* thermalization length is irrelevant above critical density, which we
232 * regard as 1e16 cm-3. Put a floor at the corresponding drthrm. */
233 drthrm = MAX2( 0.15, drthrm );
234
235 /***********************************************************************
236 *
237 * make sure we resolve initial structure in dense_tabden command
238 * if interpolated table we need to make sure that we resolve the
239 * initial changes in the structure
240 *
241 ***********************************************************************/
242
243 if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
244 {
245 drTabDen = 1.;
246 i = 1;
247 factor = 0.;
248 while( i < 100 && factor < 0.05 && radius.Radius+drTabDen*2.<radius.StopThickness[0] )
249 {
250 /* check densities at ever larger dr's, until factor becomes more than 5% */
251 factor = dense.gas_phase[ipHYDROGEN]/
252 dense_tabden(radius.Radius+drTabDen, drTabDen );
253 /* density change can be positive or negative sign */
254 factor = fabs(factor-1.);
255 drTabDen *= 2.;
256 i += 1;
257 }
258 drTabDen /= 2.;
259 }
260 else
261 {
262 drTabDen = 1e30;
263 }
264
265 /* >>chng 03 mar 20, add check on lyman band optical depth - want first zone
266 * to be thin in H2 bands */
267 /* some tests are fully molecular with solomon process turned off,
268 * do not sense this when already almost fully molecular */
269 if( hmi.H2_total/dense.gas_phase[ipHYDROGEN] < 0.1 )
270 {
271 change = 0.1;
272 }
273 else
274 {
275 /* >>chng 04 mar 14, this branch, H is quite molecular,
276 * still do not want large changes in solomon rate since linearization
277 * would not work in hmole network, bu do not need such fine steps */
278 change = 1.;
279 }
280
281 /* >>chng 04 mar 14 go back to original logic since molecular
282 * pdr's had big jump in conditions from
283 * first to second zon even when most H in H2
284 change = 0.1; */
285 /* >>chng 04 apr 18, change from 0.1 to 0.001, inital zones too large in
286 * leiden test case f1 */
287 change = 0.001;
288 /* >>chng 04 mar 13, not too large when big H2 is on */
289 if( h2.lgEnabled && h2.lgEvaluated )
290 {
291 if( fabs(h2.HeatDexc)/thermal.ctot > 0.05 )
292 {
293 /* changes in H2 heating caused by changes in solomon rate
294 * would drive temperature failures */
295 /* >>chng 04 apr 18, change from 0.001 to 0.0001, inital zones too large in
296 * leiden test case f1 */
297 change = 0.0001;
298 }
299 else
300 {
301 /* >>chng 04 apr 18, change from 0.01 to 0.001, inital zones too large in
302 * leiden test case f1 */
303 change = 0.001;
304 }
305 }
306 drH2 = change / SDIV(
307 hmi.H2_total * geometry.FillFac * hmi.H2Opacity );
308
309 /* >>chng 06 feb 01, very high U ulirg models had dramatic increase in
310 * cont pre in first few zones,
311 * in constant total pressure case, don't want acceleration across first zone to
312 * be large compared with current gas pressure */
313 if( (strcmp( dense.chDenseLaw, "CPRE" )==0) && pressure.lgContRadPresOn )
314 {
315 /* radiative acceleration was evaluated in PressureTotal */
316 drContPres = 0.05 * pressure.PresTotlCurr /
317 ((double)wind.AccelTotalOutward*dense.xMassDensity*geometry.FillFac*geometry.DirectionalCosin);
318 }
319 else if( !wind.lgStatic() )
320 {
321 /* acceleration and change in v in wind */
322 double g = fabs(wind.AccelTotalOutward-wind.AccelGravity);
323 /* wind - do not let velocity change by too much */
324 drContPres = 0.05*POW2(wind.windv)/(2.*SDIV(g));
325 }
326 else
327 drContPres = 1e30;
328
329 drValues[0].dr = drOpacity;
330 drValues[1].dr = radius.Radius/20.;
331 drValues[2].dr = drStromgren;
332 drValues[3].dr = radius.StopThickness[iteration-1]/10.;
333 drValues[4].dr = drcol;
334 drValues[5].dr = dr912;
335 drValues[6].dr = drthrm;
336 drValues[7].dr = winddr;
337 drValues[8].dr = dradf;
338 drValues[9].dr = drTabDen;
339 drValues[10].dr = drH2;
340 drValues[11].dr = drContPres;
341 drValues[12].dr = dr_time_dep;
342
343 strcpy( drValues[0].whatToSay, "drOpacity" );
344 strcpy( drValues[1].whatToSay, "radius.Radius/20.");
345 strcpy( drValues[2].whatToSay, "drStromgren");
346 strcpy( drValues[3].whatToSay, "radius.StopThickness[iteration-1]/10.");
347 strcpy( drValues[4].whatToSay, "drcol");
348 strcpy( drValues[5].whatToSay, "dr912");
349 strcpy( drValues[6].whatToSay, "drthrm");
350 strcpy( drValues[7].whatToSay, "winddr");
351 strcpy( drValues[8].whatToSay, "dradf");
352 strcpy( drValues[9].whatToSay, "drTabDen");
353 strcpy( drValues[10].whatToSay, "drH2");
354 strcpy( drValues[11].whatToSay, "drContPres");
355 strcpy( drValues[12].whatToSay, "dr_time_dep");
356
357 for( i=0; i<NUM_DR_TYPES; i++ )
358 {
359 if( drValues[i].dr < drValues[indexOfSmallest].dr )
360 {
361 indexOfSmallest = i;
362 }
363 }
364
365 radius.drad = drValues[indexOfSmallest].dr;
366
367 double rfacmin = radius.lgSdrminRel ? radius.Radius : 1.;
368 /* reset if radius.drad is less than radius.sdrmin */
369 if( rfacmin*radius.sdrmin >= radius.drad )
370 {
371 radius.drad = rfacmin*radius.sdrmin;
372 /* set flag for comment if the previous line forced a larger dr than
373 * would otherwise have been chosen. will cause comment to be generated
374 * in PrtComment if set true*/
375 radius.lgDR2Big = true;
376 }
377 else
378 {
379 radius.lgDR2Big = false;
380 }
381
382 /* this min had been in the big min set above, but caused a false alarm
383 * on the lgDR2Big test above since the set dr command sets both in and max */
384 // lgSdrmaxRel true if sdrmax is relative to current radius, false if limit in cm
385 double rfacmax = radius.lgSdrmaxRel ? radius.Radius : 1.;
386 radius.drad = MIN2( rfacmax*radius.sdrmax, radius.drad );
387 radius.drad_mid_zone = radius.drad/2.;
388
389#if 0
390 /***********************************************************************
391 *
392 * we have now generated range of estimates of first thickness,
393 * now choose smallest of the group
394 *
395 ***********************************************************************/
396
397 /* radius div by 20, to prevent big change in iron ionization for high ioniz gas,
398 * this is also the ONLY place that sphericity comes in */
399 radius.drad = MIN4( MIN3( drOpacity, radius.Radius/20., drStromgren ),
400 MIN3( radius.StopThickness[iteration-1]/10., drcol, dr912 ),
401 MIN4( drthrm, winddr, dradf, drTabDen ),
402 MIN3( drH2, drContPres, dr_time_dep ) );
403
404 /* option to set lower limit to zone thickness, with set drmin command*/
405 radius.drad = MAX2( radius.drad, rfacmin*radius.sdrmin );
406
407 /* set flag for comment if the previous line forced a larger dr than
408 * would otherwise have been chosen. will cause comment to be generated
409 * in PrtComment if set true*/
410 if( fp_equal( radius.drad, rfacmin*radius.sdrmin ) )
411 {
412 radius.lgDR2Big = true;
413 }
414 else
415 {
416 radius.lgDR2Big = false;
417 }
418
419 /* this min had been in the big min set above, but caused a false alarm
420 * on the lgDR2Big test above since the set dr command sets both in and max */
421 radius.drad = MIN2( rfacmax*radius.sdrmax, radius.drad );
422#endif
423
424 radius.drad_x_fillfac = radius.drad * geometry.FillFac;
425
426 /* save dr for this iteration */
427 drad_last_iteration = radius.drad;
428
429 /* drMinimum is smallest acceptable DRAD, and is 1/100 OF DRAD(1) */
430 /* this can be turned off by GLOB command */
431 if( radius.lgDrMnOn )
432 {
433 /* >>chng 05 mar 05, drMinimum is now drad * hden, to make propro to optical depth
434 * avoid false trigger across thermal fronts
435 * add * dense.gas_phase */
436 /* NB - drMinimum not used in code - delete? */
437 radius.drMinimum = (realnum)(radius.drad * dense.gas_phase[ipHYDROGEN]/1e7);
438 }
439 else
440 {
441 radius.drMinimum = 0.;
442 }
443
444 /* if set drmin is used, make sure drMinimum (which will cause an abort) is
445 * smaller than drmin */
446 if( radius.lgSMinON )
447 {
448 /* >>chng 05 mar 05, drMinimum is now drad * hden, to make propro to optical depth
449 * avoid false trigger across thermal fronts
450 * add * dense.gas_phase */
451 /* NB - drMinimum not used in code - delete? */
452 radius.drMinimum = MIN2(radius.drMinimum * dense.gas_phase[ipHYDROGEN],
453 (realnum)(rfacmin*radius.sdrmin/10.f) );
454 }
455
456 if( trace.lgTrace )
457 {
458 fprintf( ioQQQ,
459 " radius_first called, finds dr=%13.5e drMinimum=%12.3e sdrmin=%10.2e sdrmax=%10.2e\n",
460 radius.drad, radius.drMinimum/ dense.gas_phase[ipHYDROGEN],
461 rfacmin*radius.sdrmin, rfacmax*radius.sdrmax );
462 }
463
464 if( radius.drad < SMALLFLOAT*1.1 )
465 {
466 fprintf( ioQQQ,
467 " PROBLEM radius_first detected likely insanity, found dr=%13.5e \n", radius.drad);
468 fprintf( ioQQQ,
469 " radius_first: calculation continuing but crash is likely. \n");
470 /* this sets flag that insanity has occurred */
472 }
473
474 /* all this is to only save on last iteration
475 * the save dr command is not really a save command, making this necessary
476 * lgDRon is set true if "save dr" entered */
477 if( save.lgDROn )
478 {
479 lgDoPun = true;
480 }
481 else
482 {
483 lgDoPun = false;
484 }
485
486 /* save what we decided up? */
487 if( lgDoPun )
488 {
489 /* create hash marks on second and later iterations */
490 if( iteration > 1 && save.lgDRHash )
491 {
492 static int iter_punch=-1;
493 if( iteration !=iter_punch )
494 fprintf( save.ipDRout, "%s\n",save.chHashString );
495 iter_punch = iteration;
496 }
497 /* this is common part of each line, the zone count, depth, chosen dr, and depth2go */
498 /* >>chng 05 aug 15, had printed drNext, always zero, rather the drad, which is set here */
499 fprintf( save.ipDRout , "%ld\t%.5e\t%.3e\t%.3e\t", nzone, radius.depth, radius.drad, radius.Depth2Go );
500
501 if( radius.lgDR2Big )
502 {
503 fprintf( save.ipDRout,
504 "radius_first keys from radius.sdrmin\n");
505
506 }
507 else if( fp_equal( radius.drad, rfacmax*radius.sdrmax ) )
508 {
509
510 fprintf( save.ipDRout,
511 "radius_first keys from radius.sdrmax\n");
512 }
513 else
514 {
515 ASSERT( indexOfSmallest < NUM_DR_TYPES - 1 );
516 fprintf( save.ipDRout, "radius_first keys from %s\n",
517 drValues[indexOfSmallest].whatToSay);
518 }
519
520 /* \todo 1 improve this printout and the drValues treatment above. */
521
522#if 0
523 if( fp_equal( radius.drad, drOpacity ) )
524 {
525 fprintf( save.ipDRout,
526 "radius_first keys from drOpacity, opac was %.2e at %.2e Ryd\n",
527 BigOpacity , BigOpacityAnu );
528 }
529 else if( fp_equal( radius.drad, radius.Radius/20. ) )
530 {
531 fprintf( save.ipDRout,
532 "radius_first keys from radius.Radius\n" );
533 }
534 else if( fp_equal( radius.drad, drStromgren ) )
535 {
536 fprintf( save.ipDRout,
537 "radius_first keys from drStromgren\n");
538 }
539 else if( fp_equal( radius.drad, dr_time_dep ) )
540 {
541 fprintf( save.ipDRout,
542 "radius_first keys from time dependent\n");
543 }
544 else if( fp_equal( radius.drad, radius.StopThickness[iteration-1]/10. ) )
545 {
546 fprintf( save.ipDRout,
547 "radius_first keys from radius.StopThickness[iteration-1]\n");
548 }
549 else if( fp_equal( radius.drad, drcol ) )
550 {
551 fprintf( save.ipDRout,
552 "radius_first keys from drcol\n");
553 }
554 else if( fp_equal( radius.drad, rfacmin*radius.sdrmin ) )
555 {
556 fprintf( save.ipDRout,
557 "radius_first keys from radius.sdrmin\n");
558 }
559 else if( fp_equal( radius.drad, dr912 ) )
560 {
561 fprintf( save.ipDRout,
562 "radius_first keys from dr912\n");
563 }
564 else if( fp_equal( radius.drad, rfacmax*radius.sdrmax ) )
565 {
566 fprintf( save.ipDRout,
567 "radius_first keys from radius.sdrmax\n");
568 }
569 else if( fp_equal( radius.drad, drthrm ) )
570 {
571 fprintf( save.ipDRout,
572 "radius_first keys from drthrm\n");
573 }
574 else if( fp_equal( radius.drad, winddr ) )
575 {
576 fprintf( save.ipDRout,
577 "radius_first keys from winddr\n");
578 }
579 else if( fp_equal( radius.drad, drH2 ) )
580 {
581 fprintf( save.ipDRout,
582 "radius_first keys from H2 lyman lines\n");
583 }
584 else if( fp_equal( radius.drad, dradf ) )
585 {
586 fprintf( save.ipDRout,
587 "radius_first keys from dradf\n");
588 }
589 else if( fp_equal( radius.drad, drTabDen ) )
590 {
591 fprintf( save.ipDRout,
592 "radius_first keys from drTabDen\n");
593 }
594 else if( fp_equal( radius.drad, drContPres ) )
595 {
596 fprintf( save.ipDRout,
597 "radius_first keys from radiative acceleration across zone\n");
598 }
599 else
600 {
601 fprintf( save.ipDRout, "radius_first insanity\n" );
602 fprintf( ioQQQ, "radius_first insanity, radius is %e\n" ,
603 radius.drad);
604 ShowMe();
605 }
606#endif
607
608 }
609 return;
610}
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
long int iteration
Definition cddefines.cpp:16
#define ASSERT(exp)
Definition cddefines.h:578
#define MIN4(a, b, c, d)
Definition cddefines.h:771
#define MIN2
Definition cddefines.h:761
#define MIN3(a, b, c)
Definition cddefines.h:766
#define POW2
Definition cddefines.h:929
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
NORETURN void TotalInsanity(void)
Definition service.cpp:886
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
long ipoint(double energy_ryd)
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
double dense_tabden(double r0, double depth)
t_dynamics dynamics
Definition dynamics.cpp:44
t_geometry geometry
Definition geometry.cpp:5
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
t_hmi hmi
Definition hmi.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH_LIKE
Definition iso.h:62
t_opac opac
Definition opacity.cpp:5
t_pressure pressure
Definition pressure.cpp:5
t_radius radius
Definition radius.cpp:5
void radius_first(void)
t_rfield rfield
Definition rfield.cpp:8
t_save save
Definition save.cpp:5
static double * g
Definition species2.cpp:28
t_StopCalc StopCalc
Definition stopcalc.cpp:5
t_thermal thermal
Definition thermal.cpp:5
t_trace trace
Definition trace.cpp:5
Wind wind
Definition wind.cpp:5