cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
iter_end_chk.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 /*iter_end_check after each zone by Cloudy, determines whether model is complete */
4 #include "cddefines.h"
5 /* */
6 #ifdef EPS
7 # undef EPS
8 #endif
9 #define EPS 1.00001
10 #include "lines.h"
11 #include "mole.h"
12 #include "conv.h"
13 #include "rfield.h"
14 #include "iterations.h"
15 #include "trace.h"
16 #include "dense.h"
17 #include "colden.h"
18 #include "taulines.h"
19 #include "hmi.h"
20 #include "prt.h"
21 #include "phycon.h"
22 #include "geometry.h"
23 #include "stopcalc.h"
24 #include "opacity.h"
25 #include "thermal.h"
26 #include "cooling.h"
27 #include "pressure.h"
28 #include "radius.h"
29 #include "called.h"
30 #include "wind.h"
31 #include "hcmap.h"
32 
33 /*dmpary print all coolants for some zone, as from print cooling command */
34 STATIC void dmpary(void);
35 
36 int iter_end_check(void)
37 {
38  bool lgDone,
39  lgEndFun_v,
40  lgPrinted;
41  long int i;
42  double oxy_in_grains;
43 
44  DEBUG_ENTRY( "iter_end_check()" );
45 
46  /* >>chng 05 nov 22 - NPA. Stop calculation when fraction of oxygen frozen
47  * out on grains gets too high -
48  * NB this test is not used since StopCalc.StopDepleteFrac is set to > 1 */
49  oxy_in_grains = 0.0f;
50  for(i=0;i<mole.num_comole_calc;++i)
51  {
52  /* define the abundance of oxygen frozen out on grain surfaces */
53  oxy_in_grains += (1 - COmole[i]->lgGas_Phase)*COmole[i]->hevmol*COmole[i]->nElem[ipOXYGEN];
54  }
55  /*fprintf(ioQQQ, "DEBUG oxy in grains %.2e %e %e\n",
56  oxy_in_grains ,
57  oxy_in_grains/MAX2(SMALLFLOAT,dense.gas_phase[ipOXYGEN]) , StopCalc.StopDepleteFrac );*/
58 
59  if( trace.lgTrace )
60  {
61  fprintf( ioQQQ, " iter_end_check called, zone %li.\n" , nzone);
62  }
63  ASSERT( hcmap.MapZone >= 00 || !conv.lgSearch );
64 
65  /* >>chng 97 jun 09, now called before first zone with nzone 0 */
66  if( nzone == 0 )
67  {
68  lgEndFun_v = false;
69 
70  if( trace.lgTrace )
71  {
72  fprintf( ioQQQ, " iter_end_check returns, doing nothing since zone 0.\n" );
73  }
74  return( lgEndFun_v );
75  }
76 
77  /* check whether trace is needed for this zone and iteration */
78  if( (nzone >= trace.nznbug && iteration >= trace.npsbug) && trace.lgTrOvrd )
79  {
80  if( trace.nTrConvg==0 )
81  {
82  geometry.nprint = 1;
83  trace.lgTrace = true;
84  }
85  else
86  /* trace convergence entered = but with negative flag = make positive,
87  * abs and not mult by -1 since may trigger more than one time */
88  trace.nTrConvg = abs( trace.nTrConvg );
89  }
90 
91  /* option to turn printout on after certain zone */
92  if( prt.lgPrtStart && prt.nstart == nzone )
93  {
94  called.lgTalk = true;
95  }
96 
97  /* check whether model is done, various criteria used. */
98  lgDone = false;
99 
100  /* this is flag to check whether stopping reason was bad */
101  conv.lgBadStop = false;
102 
103  /* set temperature floor option -
104  * go to constant temperature calculation if temperature
105  * falls below floor */
106  if( phycon.te < StopCalc.TeFloor )
107  {
111  TempChange(thermal.ConstTemp , false);
112  TotalInsanity();
113  }
114 
115  /* check on radiation pressure - constant pressure unstable if dominated
116  * by radiation pressure */
117  if( (pressure.lgPres_radiation_ON && pressure.pbeta > 1.0) &&
118  (strcmp(dense.chDenseLaw ,"CPRE") == 0) &&
119  /* >>chng 03 aug 20, check on extreme values of pbeta, and abort if true */
120  (iterations.lgLastIt||(pressure.pbeta>1000.)) &&
121  /* >>chng 03 aug 19, add check on pbeta, and abort even if "no abort"
122  * was specified, since rad pres dominated limit can lead to VERY
123  * small H densities and crash due to underflow */
125  {
126  /* const total pres model; if RadPres>PGAS, then unstable, stop */
127  if( called.lgTalk )
128  {
129  fprintf( ioQQQ, "\n STOP since P(rad)/P(gas)=%7.3f >1.0\n",
130  pressure.pbeta );
131 
132  fprintf( ioQQQ, " Line contributors to radiation pressure follows:\n" );
133  PrtLinePres();
134  }
135  lgDone = true;
136  conv.lgBadStop = true;
137  strncpy( StopCalc.chReasonStop, "of radiation pressure.", sizeof(StopCalc.chReasonStop) );
138  }
139 
140  /* radius and resulting volume too large for this cpu */
142  {
143  /* too big */
144  lgDone = true;
145  strncpy( StopCalc.chReasonStop, "volume too large for this cpu.", sizeof(StopCalc.chReasonStop) );
146  }
147  /* supersonic outflowing wind, initial velocity, windv0, was > 0,
148  * but it has coasted to a stop - lgVelPos is false */
149  else if( !wind.lgVelPos && wind.windv0 > 0.)
150  {
151  /* wind solution with negative velocity */
152  lgDone = true;
153  strncpy( StopCalc.chReasonStop, "wind veloc too small.", sizeof(StopCalc.chReasonStop) );
154  }
155 
156  else if( wind.windv!=0. && fabs(wind.windv) < StopCalc.StopVelocity )
157  {
158  /* stop if absolute value of velocity falls below value */
159  lgDone = true;
160  strncpy( StopCalc.chReasonStop, "wind V too small.", sizeof(StopCalc.chReasonStop) );
161  }
162 
164  {
165  /* stop if exceed number of calls to conv base set with
166  * stop nTotalIonizStop command */
167  lgDone = true;
168  strncpy( StopCalc.chReasonStop, "nTotalIonizStop reached.", sizeof(StopCalc.chReasonStop) );
169  }
170 
171  /* this flag says that 21cm line optical depth is the stop quantity */
172  else if( StopCalc.lgStop21cm && (HFLines[0].Emis->TauCon >= StopCalc.tauend) )
173  {
174  lgDone = true;
175  strncpy( StopCalc.chReasonStop, "21 cm optical depth.", sizeof(StopCalc.chReasonStop) );
176  }
177 
179  {
180  /* stop at specified AV for (1-g) in scattering opacity */
181  lgDone = true;
182  strncpy( StopCalc.chReasonStop, "A_V reached.", sizeof(StopCalc.chReasonStop) );
183  }
184 
186  {
187  /* stop at specified AV without (1-g) in scattering opacity */
188  lgDone = true;
189  strncpy( StopCalc.chReasonStop, "A_V reached.", sizeof(StopCalc.chReasonStop) );
190  }
191 
192  else if( StopCalc.xMass!=0. &&
193  log10(SDIV(dense.xMassTotal))+1.0992099+ 2.*log10(radius.rinner) >= StopCalc.xMass )
194  {
195  /* stop at specified AV without (1-g) in scattering opacity */
196  lgDone = true;
197  strncpy( StopCalc.chReasonStop, "mass reached.", sizeof(StopCalc.chReasonStop) );
198  }
199 
200  /* >>chg 02 may 31, added pressure.lgSonicPoint logic */
201  /* WJH 19 May 2004: for some models, we want to get through the
202  * sonic point and out the other side */
204  {
205  /* D-critical solution reached sonic point */
206  lgDone = true;
207  strncpy( StopCalc.chReasonStop, "sonic point reached.", sizeof(StopCalc.chReasonStop) );
208  }
209 
210  else if( dense.EdenTrue==0 )
211  {
212  /* calculation failed */
213  conv.lgBadStop = true;
214  lgDone = true;
215  strncpy( StopCalc.chReasonStop, "zero electron density.", sizeof(StopCalc.chReasonStop) );
216  }
217 
218  else if( radius.lgdR2Small )
219  {
220  lgDone = true;
221  conv.lgBadStop = true;
222  strncpy( StopCalc.chReasonStop, "DR small rel to thick.", sizeof(StopCalc.chReasonStop) );
223  }
224 
225  else if( dense.eden < StopCalc.StopElecDensity )
226  {
227  lgDone = true;
228  strncpy( StopCalc.chReasonStop, "lowest EDEN reached.", sizeof(StopCalc.chReasonStop) );
229  }
230 
232  {
233  lgDone = true;
234  strncpy( StopCalc.chReasonStop, "low electron fraction.", sizeof(StopCalc.chReasonStop) );
235  }
236 
237  /* >>chng 05 nov 22, NA add this stop condition - stop when too many molecules
238  * are ices or solids on grains - the limit StopCalc.StopDepleteFrac is changed
239  * with the stop molecular depletion command */
240  else if( dense.lgElmtOn[ipOXYGEN] &&
242  {
243  lgDone = true;
244  strncpy( StopCalc.chReasonStop, "freeze out fraction.", sizeof(StopCalc.chReasonStop) );
245  }
246 
247  /*else if( 2.*hmi.Hmolec[ipMH2g]/dense.gas_phase[ipHYDROGEN] < StopCalc.StopH2MoleFrac )*/
249  {
250  lgDone = true;
251  strncpy( StopCalc.chReasonStop, "large H_2/H fraction.", sizeof(StopCalc.chReasonStop) );
252  }
253 
256  {
257  lgDone = true;
258  strncpy( StopCalc.chReasonStop, "low H_+/H fraction.", sizeof(StopCalc.chReasonStop) );
259  }
260 
261  else if( radius.lgDrMinUsed )
262  {
263  /* dr became too small */
264  conv.lgBadStop = true;
265  lgDone = true;
266  strncpy( StopCalc.chReasonStop, "DRAD small- set DRMIN.", sizeof(StopCalc.chReasonStop) );
267  }
268 
269  else if( radius.depth >= radius.router[iteration-1]/EPS || radius.lgDrNeg )
270  {
271  lgDone = true;
272  strncpy( StopCalc.chReasonStop, "outer radius reached.", sizeof(StopCalc.chReasonStop) );
273  }
274 
275  else if( StopCalc.iptnu >= 0 &&
277  {
278  lgDone = true;
279  strncpy( StopCalc.chReasonStop, "optical depth reached.", sizeof(StopCalc.chReasonStop) );
280  }
281 
283  {
284  /* StopCalc.HColStop default set to COLUMN_INIT == 1e30 */
285  lgDone = true;
286  strncpy( StopCalc.chReasonStop, "H column dens reached.", sizeof(StopCalc.chReasonStop) );
287  }
288 
289  else if( colden.colden[ipCOL_Hp] >= StopCalc.colpls/EPS )
290  {
291  lgDone = true;
292  strncpy( StopCalc.chReasonStop, "H+ column dens reached.", sizeof(StopCalc.chReasonStop) );
293  }
294 
296  {
297  /* >>chng 03 apr 15, add molecular hydrogen */
298  lgDone = true;
299  strncpy( StopCalc.chReasonStop, "H2 column dens reached.", sizeof(StopCalc.chReasonStop) );
300  }
301 
303  {
304  /* >>chng 04 feb 10, stopping command for H2 + H I */
305  lgDone = true;
306  strncpy( StopCalc.chReasonStop, "H2+H0 column dens reached.", sizeof(StopCalc.chReasonStop) );
307  }
308 
310  {
311  /* >>chng 05 jan 09, stopping command for N(H0) / Tspin */
312  lgDone = true;
313  strncpy( StopCalc.chReasonStop, "N(H0)/Tspin column dens reached.", sizeof(StopCalc.chReasonStop) );
314  }
315 
316  else if( findspecies("CO")->hevcol >= StopCalc.col_monoxco/EPS )
317  {
318  /* >>chng 03 oct 27--Nick Abel, add carbon monoxide */
319  lgDone = true;
320  strncpy( StopCalc.chReasonStop, "CO column dens reached.", sizeof(StopCalc.chReasonStop) );
321  }
322 
323  else if( colden.colden[ipCOL_H0] >= StopCalc.colnut/EPS )
324  {
325  lgDone = true;
326  strncpy( StopCalc.chReasonStop, "H0 column dens reached.", sizeof(StopCalc.chReasonStop) );
327  }
328 
329  /* >>chng 99 jul 7, when no h2 molecules, include h2 as neutral atomic hydrogen,
330  * hmi.lgNoH2Mole is set false in zero, true with command no hydrogen molecules */
331  else if( hmi.lgNoH2Mole &&
333  {
334  lgDone = true;
335  strncpy( StopCalc.chReasonStop, "H0-H0+H2 column dens reached.", sizeof(StopCalc.chReasonStop) );
336  }
337 
338  else if( phycon.te > StopCalc.T2High )
339  {
340  lgDone = true;
341  strncpy( StopCalc.chReasonStop, "highest Te reached.", sizeof(StopCalc.chReasonStop) );
342  }
343 
344  else if( phycon.te < StopCalc.tend )
345  {
346  lgDone = true;
347  strncpy( StopCalc.chReasonStop, "lowest Te reached.", sizeof(StopCalc.chReasonStop) );
348  }
349 
350  else if( nzone >= geometry.nend[iteration-1] )
351  {
352  lgDone = true;
353  geometry.lgZoneTrp = true;
354  strncpy( StopCalc.chReasonStop, "NZONE reached.", sizeof(StopCalc.chReasonStop) );
355  }
356 
357  /* option to stop calculation when line intensity ratio reaches certain value,
358  * nstpl is number of stop line commands entered */
359  else if( StopCalc.nstpl > 0 )
360  {
361  /* line ratio exceeded maximum permitted value
362  * do not consider case where norm line has zero intensity */
363  for( i=0; i < StopCalc.nstpl; i++ )
364  {
365  /* the second line is always set to something, default is H beta */
367  {
368  char chString[10];
371  {
372  lgDone = true;
373  sprt_wl( chString , StopCalc.StopLineWl1[i] );
374  sprintf( StopCalc.chReasonStop, "line %s %s reached",
375  StopCalc.chStopLabel1[i] , chString );
376  }
377  }
378  }
379  }
380 
381  else if( radius.drNext <= 0. )
382  {
383  /* this cant happen */
384  if( called.lgTalk )
385  {
386  fprintf( ioQQQ, " drNext=%10.2e STOP\n", radius.drNext );
387  }
388  lgDone = true;
389  conv.lgBadStop = true;
390  strncpy( StopCalc.chReasonStop, "internal error - DRAD.", sizeof(StopCalc.chReasonStop) );
391  ShowMe();
392  }
393 
394  else if( lgAbort )
395  {
396  /* calculation failed */
397  conv.lgBadStop = true;
398  lgDone = true;
399  strncpy( StopCalc.chReasonStop, "calculation aborted.", sizeof(StopCalc.chReasonStop) );
400  }
401 
402  lgPrinted = false;
403  if( lgDone )
404  {
405  /* flag to call it quits */
406  lgEndFun_v = true;
407  PrtZone();
408  lgPrinted = true;
409  }
410 
411  else
412  {
413  /* passed all the tests, keep going */
414  /* check whether this zone should be printed */
415  if( ((nzone/geometry.nprint)*geometry.nprint == nzone ||
416  nzone == 1) || trace.nTrConvg )
417  {
418  PrtZone();
419  lgPrinted = true;
420  }
421  /* flag to keep going */
422  lgEndFun_v = false;
423  }
424 
425  /* dump cooling arrays for this zone? */
426  if( prt.nzdump == nzone || prt.nzdump == 0 )
427  dmpary();
428 
429  /* do map of cooling function if desired, and not yet done */
430  /* >>chng 02 may 29, MapZone < = to <= 0 - map 0 did not work */
431  /* >>chng 04 jun 16, change to MapZone = 0 for map of first zone then quit,
432  * -1 is not set, positive, do map of that zone */
433  if( !hcmap.lgMapDone && (hcmap.MapZone == 0 || nzone == hcmap.MapZone) )
434  {
435  /* print last zone if not already done */
436  if( !lgPrinted )
437  {
438  PrtZone();
439  }
440 
441  /* say that we are doing a map */
442  hcmap.lgMapBeingDone = true;
443 
444  /* save old output file then redirect to map file */
445  /* >>chng 01 mar 28, ioMAP may not be initialized, PvH */
446  if( ioMAP != NULL )
447  map_do(ioMAP, " map");
448  else
449  map_do(ioQQQ, " map");
450 
451  /* stop after doing map */
452  lgEndFun_v = true;
453  strncpy( StopCalc.chReasonStop, "MAP command used-stop.", sizeof(StopCalc.chReasonStop) );
454 
455  /* >>chng 03 jun 06, reset iterations since we want to stop even if
456  * iterate xx is specified, bug caught by Joop Schaye */
457  iterations.itermx = 0;
458 
459  /* make really sure that the string contained in StopCalc.chReasonStop is properly terminated */
460  StopCalc.chReasonStop[sizeof(StopCalc.chReasonStop)-1] = '\0';
461 
462  if( trace.lgTrace )
463  {
464  fprintf( ioQQQ, " iter_end_check returns after map.\n" );
465  }
466  return( lgEndFun_v );
467  }
468 
469  if( lgEndFun_v && prt.lgOnlyZone )
470  {
471  cdEXIT(EXIT_FAILURE);
472  }
473 
474  /* the string contained in StopCalc.chReasonStop must be properly
475  * terminated -this can't fail - strlen returns the number of characters
476  * in str, excluding the terminal NULL.*/
477  if( strlen( StopCalc.chReasonStop ) >= nCHREASONSTOP-1 )
478  TotalInsanity();
479 
480  if( trace.lgTrace )
481  {
482  fprintf( ioQQQ, " iter_end_check bottom return.\n" );
483  }
484  return( lgEndFun_v );
485 }
486 
487 #ifdef EPS
488 # undef EPS
489 #endif
490 #define EPS 0.005
491 /*dmpary print all coolants for some zone, as from print cooling command */
492 STATIC void dmpary(void)
493 {
494  long int i;
495  realnum ratio;
496 
497  DEBUG_ENTRY( "dmpary()" );
498 
499  fprintf( ioQQQ,
500  " This is a print out of the cooling array for zone number %3ld\n",
501  nzone );
502 
503  fprintf( ioQQQ,
504  " The total heating was HTOT=%10.2e erg/s/cm3, the total cooling was CTOT=%10.2e, and the temperature was%10.3eK.\n",
506 
507  fprintf( ioQQQ,
508  " All coolants greater than%6.2f%% of the total will be printed.\n",
509  EPS*100. );
510 
511  /* flag all significant coolants */
512  coolpr(ioQQQ,"ZERO",1,0.,"ZERO");
513  for( i=0; i < thermal.ncltot; i++ )
514  {
515  ratio = (realnum)(thermal.cooling[i]/thermal.ctot);
516  if( fabs(ratio) > EPS )
517  {
519  ratio,"DOIT");
520  }
521 
522  ratio = (realnum)(thermal.heatnt[i]/thermal.ctot);
523  if( fabs(ratio) > EPS )
524  {
526  ratio,"DOIT");
527  }
528  }
529  coolpr(ioQQQ,"DONE",1,0.,"DONE");
530  return;
531 }

Generated for cloudy by doxygen 1.8.3.1