cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*iter_end_check after each zone by Cloudy, determines whether model is complete */ 00004 #include "cddefines.h" 00005 /* */ 00006 #ifdef EPS 00007 # undef EPS 00008 #endif 00009 #define EPS 1.00001 00010 #include "lines.h" 00011 #include "mole.h" 00012 #include "conv.h" 00013 #include "rfield.h" 00014 #include "iterations.h" 00015 #include "trace.h" 00016 #include "dense.h" 00017 #include "colden.h" 00018 #include "taulines.h" 00019 #include "hmi.h" 00020 #include "prt.h" 00021 #include "phycon.h" 00022 #include "geometry.h" 00023 #include "stopcalc.h" 00024 #include "opacity.h" 00025 #include "thermal.h" 00026 #include "cooling.h" 00027 #include "pressure.h" 00028 #include "radius.h" 00029 #include "called.h" 00030 #include "wind.h" 00031 #include "hcmap.h" 00032 00033 /*dmpary print all coolants for some zone, as from print cooling command */ 00034 STATIC void dmpary(void); 00035 00036 int iter_end_check(void) 00037 { 00038 bool lgDone, 00039 lgEndFun_v, 00040 lgPrinted; 00041 long int i; 00042 double oxy_in_grains; 00043 00044 DEBUG_ENTRY( "iter_end_check()" ); 00045 00046 /* >>chng 05 nov 22 - NPA. Stop calculation when fraction of oxygen frozen 00047 * out on grains gets too high - 00048 * NB this test is not used since StopCalc.StopDepleteFrac is set to > 1 */ 00049 oxy_in_grains = 0.0f; 00050 for(i=0;i<mole.num_comole_calc;++i) 00051 { 00052 /* define the abundance of oxygen frozen out on grain surfaces */ 00053 oxy_in_grains += (1 - COmole[i]->lgGas_Phase)*COmole[i]->hevmol*COmole[i]->nElem[ipOXYGEN]; 00054 } 00055 /*fprintf(ioQQQ, "DEBUG oxy in grains %.2e %e %e\n", 00056 oxy_in_grains , 00057 oxy_in_grains/MAX2(SMALLFLOAT,dense.gas_phase[ipOXYGEN]) , StopCalc.StopDepleteFrac );*/ 00058 00059 if( trace.lgTrace ) 00060 { 00061 fprintf( ioQQQ, " iter_end_check called, zone %li.\n" , nzone); 00062 } 00063 ASSERT( hcmap.MapZone >= 00 || !conv.lgSearch ); 00064 00065 /* >>chng 97 jun 09, now called before first zone with nzone 0 */ 00066 if( nzone == 0 ) 00067 { 00068 lgEndFun_v = false; 00069 00070 if( trace.lgTrace ) 00071 { 00072 fprintf( ioQQQ, " iter_end_check returns, doing nothing since zone 0.\n" ); 00073 } 00074 return( lgEndFun_v ); 00075 } 00076 00077 /* check whether trace is needed for this zone and iteration */ 00078 if( (nzone >= trace.nznbug && iteration >= trace.npsbug) && trace.lgTrOvrd ) 00079 { 00080 if( trace.nTrConvg==0 ) 00081 { 00082 geometry.nprint = 1; 00083 trace.lgTrace = true; 00084 } 00085 else 00086 /* trace convergence entered = but with negative flag = make positive, 00087 * abs and not mult by -1 since may trigger more than one time */ 00088 trace.nTrConvg = abs( trace.nTrConvg ); 00089 } 00090 00091 /* option to turn printout on after certain zone */ 00092 if( prt.lgPrtStart && prt.nstart == nzone ) 00093 { 00094 called.lgTalk = true; 00095 } 00096 00097 /* check whether model is done, various criteria used. */ 00098 lgDone = false; 00099 00100 /* this is flag to check whether stopping reason was bad */ 00101 conv.lgBadStop = false; 00102 00103 /* set temperature floor option - 00104 * go to constant temperature calculation if temperature 00105 * falls below floor */ 00106 if( phycon.te < StopCalc.TeFloor ) 00107 { 00108 thermal.lgTemperatureConstant = true; 00109 thermal.ConstTemp = (realnum)StopCalc.TeFloor; 00110 phycon.te = thermal.ConstTemp; 00111 TempChange(thermal.ConstTemp , false); 00112 TotalInsanity(); 00113 } 00114 00115 /* check on radiation pressure - constant pressure unstable if dominated 00116 * by radiation pressure */ 00117 if( (pressure.lgPres_radiation_ON && pressure.pbeta > 1.0) && 00118 (strcmp(dense.chDenseLaw ,"CPRE") == 0) && 00119 /* >>chng 03 aug 20, check on extreme values of pbeta, and abort if true */ 00120 (iterations.lgLastIt||(pressure.pbeta>1000.)) && 00121 /* >>chng 03 aug 19, add check on pbeta, and abort even if "no abort" 00122 * was specified, since rad pres dominated limit can lead to VERY 00123 * small H densities and crash due to underflow */ 00124 (pressure.lgRadPresAbortOK||(pressure.pbeta>1000.)) ) 00125 { 00126 /* const total pres model; if RadPres>PGAS, then unstable, stop */ 00127 if( called.lgTalk ) 00128 { 00129 fprintf( ioQQQ, "\n STOP since P(rad)/P(gas)=%7.3f >1.0\n", 00130 pressure.pbeta ); 00131 00132 fprintf( ioQQQ, " Line contributors to radiation pressure follows:\n" ); 00133 PrtLinePres(); 00134 } 00135 lgDone = true; 00136 conv.lgBadStop = true; 00137 strncpy( StopCalc.chReasonStop, "of radiation pressure.", sizeof(StopCalc.chReasonStop) ); 00138 } 00139 00140 /* radius and resulting volume too large for this cpu */ 00141 if( radius.drad_x_fillfac*radius.r1r0sq > BIGFLOAT/10.) 00142 { 00143 /* too big */ 00144 lgDone = true; 00145 strncpy( StopCalc.chReasonStop, "volume too large for this cpu.", sizeof(StopCalc.chReasonStop) ); 00146 } 00147 /* supersonic outflowing wind, initial velocity, windv0, was > 0, 00148 * but it has coasted to a stop - lgVelPos is false */ 00149 else if( !wind.lgVelPos && wind.windv0 > 0.) 00150 { 00151 /* wind solution with negative velocity */ 00152 lgDone = true; 00153 strncpy( StopCalc.chReasonStop, "wind veloc too small.", sizeof(StopCalc.chReasonStop) ); 00154 } 00155 00156 else if( wind.windv!=0. && fabs(wind.windv) < StopCalc.StopVelocity ) 00157 { 00158 /* stop if absolute value of velocity falls below value */ 00159 lgDone = true; 00160 strncpy( StopCalc.chReasonStop, "wind V too small.", sizeof(StopCalc.chReasonStop) ); 00161 } 00162 00163 else if( (StopCalc.nTotalIonizStop>0) && conv.nTotalIoniz> StopCalc.nTotalIonizStop ) 00164 { 00165 /* stop if exceed number of calls to conv base set with 00166 * stop nTotalIonizStop command */ 00167 lgDone = true; 00168 strncpy( StopCalc.chReasonStop, "nTotalIonizStop reached.", sizeof(StopCalc.chReasonStop) ); 00169 } 00170 00171 /* this flag says that 21cm line optical depth is the stop quantity */ 00172 else if( StopCalc.lgStop21cm && (HFLines[0].Emis->TauCon >= StopCalc.tauend) ) 00173 { 00174 lgDone = true; 00175 strncpy( StopCalc.chReasonStop, "21 cm optical depth.", sizeof(StopCalc.chReasonStop) ); 00176 } 00177 00178 else if( rfield.extin_mag_V_extended >= StopCalc.AV_extended ) 00179 { 00180 /* stop at specified AV for (1-g) in scattering opacity */ 00181 lgDone = true; 00182 strncpy( StopCalc.chReasonStop, "A_V reached.", sizeof(StopCalc.chReasonStop) ); 00183 } 00184 00185 else if( rfield.extin_mag_V_point >= StopCalc.AV_point ) 00186 { 00187 /* stop at specified AV without (1-g) in scattering opacity */ 00188 lgDone = true; 00189 strncpy( StopCalc.chReasonStop, "A_V reached.", sizeof(StopCalc.chReasonStop) ); 00190 } 00191 00192 else if( StopCalc.xMass!=0. && 00193 log10(SDIV(dense.xMassTotal))+1.0992099+ 2.*log10(radius.rinner) >= StopCalc.xMass ) 00194 { 00195 /* stop at specified AV without (1-g) in scattering opacity */ 00196 lgDone = true; 00197 strncpy( StopCalc.chReasonStop, "mass reached.", sizeof(StopCalc.chReasonStop) ); 00198 } 00199 00200 /* >>chg 02 may 31, added pressure.lgSonicPoint logic */ 00201 /* WJH 19 May 2004: for some models, we want to get through the 00202 * sonic point and out the other side */ 00203 else if( pressure.lgSonicPoint && pressure.lgSonicPointAbortOK ) 00204 { 00205 /* D-critical solution reached sonic point */ 00206 lgDone = true; 00207 strncpy( StopCalc.chReasonStop, "sonic point reached.", sizeof(StopCalc.chReasonStop) ); 00208 } 00209 00210 else if( dense.EdenTrue==0 ) 00211 { 00212 /* calculation failed */ 00213 conv.lgBadStop = true; 00214 lgDone = true; 00215 strncpy( StopCalc.chReasonStop, "zero electron density.", sizeof(StopCalc.chReasonStop) ); 00216 } 00217 00218 else if( radius.lgdR2Small ) 00219 { 00220 lgDone = true; 00221 conv.lgBadStop = true; 00222 strncpy( StopCalc.chReasonStop, "DR small rel to thick.", sizeof(StopCalc.chReasonStop) ); 00223 } 00224 00225 else if( dense.eden < StopCalc.StopElecDensity ) 00226 { 00227 lgDone = true; 00228 strncpy( StopCalc.chReasonStop, "lowest EDEN reached.", sizeof(StopCalc.chReasonStop) ); 00229 } 00230 00231 else if( dense.eden/dense.gas_phase[ipHYDROGEN] < StopCalc.StopElecFrac ) 00232 { 00233 lgDone = true; 00234 strncpy( StopCalc.chReasonStop, "low electron fraction.", sizeof(StopCalc.chReasonStop) ); 00235 } 00236 00237 /* >>chng 05 nov 22, NA add this stop condition - stop when too many molecules 00238 * are ices or solids on grains - the limit StopCalc.StopDepleteFrac is changed 00239 * with the stop molecular depletion command */ 00240 else if( dense.lgElmtOn[ipOXYGEN] && 00241 (oxy_in_grains/MAX2(SMALLFLOAT,dense.gas_phase[ipOXYGEN])) > StopCalc.StopDepleteFrac ) 00242 { 00243 lgDone = true; 00244 strncpy( StopCalc.chReasonStop, "freeze out fraction.", sizeof(StopCalc.chReasonStop) ); 00245 } 00246 00247 /*else if( 2.*hmi.Hmolec[ipMH2g]/dense.gas_phase[ipHYDROGEN] < StopCalc.StopH2MoleFrac )*/ 00248 else if( 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN] > StopCalc.StopH2MoleFrac ) 00249 { 00250 lgDone = true; 00251 strncpy( StopCalc.chReasonStop, "large H_2/H fraction.", sizeof(StopCalc.chReasonStop) ); 00252 } 00253 00254 else if( dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] < 00255 StopCalc.StopHPlusFrac ) 00256 { 00257 lgDone = true; 00258 strncpy( StopCalc.chReasonStop, "low H_+/H fraction.", sizeof(StopCalc.chReasonStop) ); 00259 } 00260 00261 else if( radius.lgDrMinUsed ) 00262 { 00263 /* dr became too small */ 00264 conv.lgBadStop = true; 00265 lgDone = true; 00266 strncpy( StopCalc.chReasonStop, "DRAD small- set DRMIN.", sizeof(StopCalc.chReasonStop) ); 00267 } 00268 00269 else if( radius.depth >= radius.router[iteration-1]/EPS || radius.lgDrNeg ) 00270 { 00271 lgDone = true; 00272 strncpy( StopCalc.chReasonStop, "outer radius reached.", sizeof(StopCalc.chReasonStop) ); 00273 } 00274 00275 else if( StopCalc.iptnu >= 0 && 00276 opac.TauAbsGeo[0][StopCalc.iptnu-1] >= StopCalc.tauend/EPS ) 00277 { 00278 lgDone = true; 00279 strncpy( StopCalc.chReasonStop, "optical depth reached.", sizeof(StopCalc.chReasonStop) ); 00280 } 00281 00282 else if( colden.colden[ipCOL_HTOT] >= StopCalc.HColStop/EPS ) 00283 { 00284 /* StopCalc.HColStop default set to COLUMN_INIT == 1e30 */ 00285 lgDone = true; 00286 strncpy( StopCalc.chReasonStop, "H column dens reached.", sizeof(StopCalc.chReasonStop) ); 00287 } 00288 00289 else if( colden.colden[ipCOL_Hp] >= StopCalc.colpls/EPS ) 00290 { 00291 lgDone = true; 00292 strncpy( StopCalc.chReasonStop, "H+ column dens reached.", sizeof(StopCalc.chReasonStop) ); 00293 } 00294 00295 else if( (colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]) >= StopCalc.col_h2/EPS ) 00296 { 00297 /* >>chng 03 apr 15, add molecular hydrogen */ 00298 lgDone = true; 00299 strncpy( StopCalc.chReasonStop, "H2 column dens reached.", sizeof(StopCalc.chReasonStop) ); 00300 } 00301 00302 else if( (2.*(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]) + colden.colden[ipCOL_H0]) >= StopCalc.col_h2_nut/EPS ) 00303 { 00304 /* >>chng 04 feb 10, stopping command for H2 + H I */ 00305 lgDone = true; 00306 strncpy( StopCalc.chReasonStop, "H2+H0 column dens reached.", sizeof(StopCalc.chReasonStop) ); 00307 } 00308 00309 else if( colden.H0_ov_Tspin >= StopCalc.col_H0_ov_Tspin/EPS ) 00310 { 00311 /* >>chng 05 jan 09, stopping command for N(H0) / Tspin */ 00312 lgDone = true; 00313 strncpy( StopCalc.chReasonStop, "N(H0)/Tspin column dens reached.", sizeof(StopCalc.chReasonStop) ); 00314 } 00315 00316 else if( findspecies("CO")->hevcol >= StopCalc.col_monoxco/EPS ) 00317 { 00318 /* >>chng 03 oct 27--Nick Abel, add carbon monoxide */ 00319 lgDone = true; 00320 strncpy( StopCalc.chReasonStop, "CO column dens reached.", sizeof(StopCalc.chReasonStop) ); 00321 } 00322 00323 else if( colden.colden[ipCOL_H0] >= StopCalc.colnut/EPS ) 00324 { 00325 lgDone = true; 00326 strncpy( StopCalc.chReasonStop, "H0 column dens reached.", sizeof(StopCalc.chReasonStop) ); 00327 } 00328 00329 /* >>chng 99 jul 7, when no h2 molecules, include h2 as neutral atomic hydrogen, 00330 * hmi.lgNoH2Mole is set false in zero, true with command no hydrogen molecules */ 00331 else if( hmi.lgNoH2Mole && 00332 ( (colden.colden[ipCOL_H0]+2.*(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])) >= StopCalc.colnut/EPS) ) 00333 { 00334 lgDone = true; 00335 strncpy( StopCalc.chReasonStop, "H0-H0+H2 column dens reached.", sizeof(StopCalc.chReasonStop) ); 00336 } 00337 00338 else if( phycon.te > StopCalc.T2High ) 00339 { 00340 lgDone = true; 00341 strncpy( StopCalc.chReasonStop, "highest Te reached.", sizeof(StopCalc.chReasonStop) ); 00342 } 00343 00344 else if( phycon.te < StopCalc.tend ) 00345 { 00346 lgDone = true; 00347 strncpy( StopCalc.chReasonStop, "lowest Te reached.", sizeof(StopCalc.chReasonStop) ); 00348 } 00349 00350 else if( nzone >= geometry.nend[iteration-1] ) 00351 { 00352 lgDone = true; 00353 geometry.lgZoneTrp = true; 00354 strncpy( StopCalc.chReasonStop, "NZONE reached.", sizeof(StopCalc.chReasonStop) ); 00355 } 00356 00357 /* option to stop calculation when line intensity ratio reaches certain value, 00358 * nstpl is number of stop line commands entered */ 00359 else if( StopCalc.nstpl > 0 ) 00360 { 00361 /* line ratio exceeded maximum permitted value 00362 * do not consider case where norm line has zero intensity */ 00363 for( i=0; i < StopCalc.nstpl; i++ ) 00364 { 00365 /* the second line is always set to something, default is H beta */ 00366 if( LineSv[StopCalc.ipStopLin2[i]].sumlin[LineSave.lgLineEmergent] > 0. ) 00367 { 00368 char chString[10]; 00369 if( LineSv[StopCalc.ipStopLin1[i]].sumlin[LineSave.lgLineEmergent]/ 00370 LineSv[StopCalc.ipStopLin2[i]].sumlin[LineSave.lgLineEmergent] > StopCalc.stpint[i] ) 00371 { 00372 lgDone = true; 00373 sprt_wl( chString , StopCalc.StopLineWl1[i] ); 00374 sprintf( StopCalc.chReasonStop, "line %s %s reached", 00375 StopCalc.chStopLabel1[i] , chString ); 00376 } 00377 } 00378 } 00379 } 00380 00381 else if( radius.drNext <= 0. ) 00382 { 00383 /* this cant happen */ 00384 if( called.lgTalk ) 00385 { 00386 fprintf( ioQQQ, " drNext=%10.2e STOP\n", radius.drNext ); 00387 } 00388 lgDone = true; 00389 conv.lgBadStop = true; 00390 strncpy( StopCalc.chReasonStop, "internal error - DRAD.", sizeof(StopCalc.chReasonStop) ); 00391 ShowMe(); 00392 } 00393 00394 else if( lgAbort ) 00395 { 00396 /* calculation failed */ 00397 conv.lgBadStop = true; 00398 lgDone = true; 00399 strncpy( StopCalc.chReasonStop, "calculation aborted.", sizeof(StopCalc.chReasonStop) ); 00400 } 00401 00402 lgPrinted = false; 00403 if( lgDone ) 00404 { 00405 /* flag to call it quits */ 00406 lgEndFun_v = true; 00407 PrtZone(); 00408 lgPrinted = true; 00409 } 00410 00411 else 00412 { 00413 /* passed all the tests, keep going */ 00414 /* check whether this zone should be printed */ 00415 if( ((nzone/geometry.nprint)*geometry.nprint == nzone || 00416 nzone == 1) || trace.nTrConvg ) 00417 { 00418 PrtZone(); 00419 lgPrinted = true; 00420 } 00421 /* flag to keep going */ 00422 lgEndFun_v = false; 00423 } 00424 00425 /* dump cooling arrays for this zone? */ 00426 if( prt.nzdump == nzone || prt.nzdump == 0 ) 00427 dmpary(); 00428 00429 /* do map of cooling function if desired, and not yet done */ 00430 /* >>chng 02 may 29, MapZone < = to <= 0 - map 0 did not work */ 00431 /* >>chng 04 jun 16, change to MapZone = 0 for map of first zone then quit, 00432 * -1 is not set, positive, do map of that zone */ 00433 if( !hcmap.lgMapDone && (hcmap.MapZone == 0 || nzone == hcmap.MapZone) ) 00434 { 00435 /* print last zone if not already done */ 00436 if( !lgPrinted ) 00437 { 00438 PrtZone(); 00439 } 00440 00441 /* say that we are doing a map */ 00442 hcmap.lgMapBeingDone = true; 00443 00444 /* save old output file then redirect to map file */ 00445 /* >>chng 01 mar 28, ioMAP may not be initialized, PvH */ 00446 if( ioMAP != NULL ) 00447 map_do(ioMAP, " map"); 00448 else 00449 map_do(ioQQQ, " map"); 00450 00451 /* stop after doing map */ 00452 lgEndFun_v = true; 00453 strncpy( StopCalc.chReasonStop, "MAP command used-stop.", sizeof(StopCalc.chReasonStop) ); 00454 00455 /* >>chng 03 jun 06, reset iterations since we want to stop even if 00456 * iterate xx is specified, bug caught by Joop Schaye */ 00457 iterations.itermx = 0; 00458 00459 /* make really sure that the string contained in StopCalc.chReasonStop is properly terminated */ 00460 StopCalc.chReasonStop[sizeof(StopCalc.chReasonStop)-1] = '\0'; 00461 00462 if( trace.lgTrace ) 00463 { 00464 fprintf( ioQQQ, " iter_end_check returns after map.\n" ); 00465 } 00466 return( lgEndFun_v ); 00467 } 00468 00469 if( lgEndFun_v && prt.lgOnlyZone ) 00470 { 00471 cdEXIT(EXIT_FAILURE); 00472 } 00473 00474 /* the string contained in StopCalc.chReasonStop must be properly 00475 * terminated -this can't fail - strlen returns the number of characters 00476 * in str, excluding the terminal NULL.*/ 00477 if( strlen( StopCalc.chReasonStop ) >= nCHREASONSTOP-1 ) 00478 TotalInsanity(); 00479 00480 if( trace.lgTrace ) 00481 { 00482 fprintf( ioQQQ, " iter_end_check bottom return.\n" ); 00483 } 00484 return( lgEndFun_v ); 00485 } 00486 00487 #ifdef EPS 00488 # undef EPS 00489 #endif 00490 #define EPS 0.005 00491 /*dmpary print all coolants for some zone, as from print cooling command */ 00492 STATIC void dmpary(void) 00493 { 00494 long int i; 00495 realnum ratio; 00496 00497 DEBUG_ENTRY( "dmpary()" ); 00498 00499 fprintf( ioQQQ, 00500 " This is a print out of the cooling array for zone number %3ld\n", 00501 nzone ); 00502 00503 fprintf( ioQQQ, 00504 " The total heating was HTOT=%10.2e erg/s/cm3, the total cooling was CTOT=%10.2e, and the temperature was%10.3eK.\n", 00505 thermal.htot, thermal.ctot, phycon.te ); 00506 00507 fprintf( ioQQQ, 00508 " All coolants greater than%6.2f%% of the total will be printed.\n", 00509 EPS*100. ); 00510 00511 /* flag all significant coolants */ 00512 coolpr(ioQQQ,"ZERO",1,0.,"ZERO"); 00513 for( i=0; i < thermal.ncltot; i++ ) 00514 { 00515 ratio = (realnum)(thermal.cooling[i]/thermal.ctot); 00516 if( fabs(ratio) > EPS ) 00517 { 00518 coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i], 00519 ratio,"DOIT"); 00520 } 00521 00522 ratio = (realnum)(thermal.heatnt[i]/thermal.ctot); 00523 if( fabs(ratio) > EPS ) 00524 { 00525 coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i], 00526 ratio,"DOIT"); 00527 } 00528 } 00529 coolpr(ioQQQ,"DONE",1,0.,"DONE"); 00530 return; 00531 }