00001
00002
00003
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 "predcont.h"
00028 #include "pressure.h"
00029 #include "radius.h"
00030 #include "called.h"
00031 #include "wind.h"
00032 #include "hcmap.h"
00033
00034
00035 STATIC void dmpary(void);
00036
00037 int iter_end_check(void)
00038 {
00039 bool lgDone,
00040 lgEndFun_v,
00041 lgPrinted;
00042 long int i;
00043 double oxy_in_grains;
00044
00045 DEBUG_ENTRY( "iter_end_check()" );
00046
00047
00048
00049
00050 oxy_in_grains = 0.0f;
00051 count_ptr<chem_atom> elOxygen = unresolved_atom_list[ipOXYGEN];
00052 for(i=0;i<mole_global.num_calc;++i)
00053 {
00054
00055 if( mole_global.list[i]->nAtom.find(elOxygen) != mole_global.list[i]->nAtom.end() && mole_global.list[i]->parentLabel.empty() )
00056 oxy_in_grains += (1 - mole_global.list[i]->lgGas_Phase)*mole.species[i].den*mole_global.list[i]->nAtom[elOxygen];
00057 }
00058
00059
00060
00061
00062 if( trace.lgTrace )
00063 {
00064 fprintf( ioQQQ, " iter_end_check called, zone %li.\n" , nzone);
00065 }
00066 ASSERT( hcmap.MapZone >= 00 || !conv.lgSearch );
00067
00068
00069 if( nzone == 0 )
00070 {
00071 lgEndFun_v = false;
00072
00073 if( trace.lgTrace )
00074 {
00075 fprintf( ioQQQ, " iter_end_check returns, doing nothing since zone 0.\n" );
00076 }
00077 return( lgEndFun_v );
00078 }
00079
00080
00081 if( (nzone >= trace.nznbug && iteration >= trace.npsbug) && trace.lgTrOvrd )
00082 {
00083 if( trace.nTrConvg==0 )
00084 {
00085 geometry.nprint = 1;
00086 trace.lgTrace = true;
00087 }
00088 else
00089
00090
00091 trace.nTrConvg = abs( trace.nTrConvg );
00092 }
00093
00094
00095 if( prt.lgPrtStart && prt.nstart == nzone )
00096 {
00097 called.lgTalk = cpu.i().lgMPI_talk();
00098 }
00099
00100
00101 lgDone = false;
00102
00103
00104 conv.lgBadStop = false;
00105
00106
00107
00108
00109 if( phycon.te < StopCalc.TeFloor )
00110 {
00111 thermal.lgTemperatureConstant = true;
00112 thermal.ConstTemp = (realnum)StopCalc.TeFloor;
00113 phycon.te = thermal.ConstTemp;
00114 TempChange(thermal.ConstTemp , false);
00115 TotalInsanity();
00116 }
00117
00118
00119
00120 if( (pressure.lgPres_radiation_ON && pressure.pbeta > 1.0) &&
00121 (strcmp(dense.chDenseLaw ,"CPRE") == 0) &&
00122
00123 (iterations.lgLastIt||(pressure.pbeta>1000.)) &&
00124
00125
00126
00127 (pressure.lgRadPresAbortOK||(pressure.pbeta>1000.)) )
00128 {
00129
00130 if( called.lgTalk )
00131 {
00132 fprintf( ioQQQ, "\n STOP since P(rad)/P(gas)=%7.3f >1.0\n",
00133 pressure.pbeta );
00134
00135 fprintf( ioQQQ, " Line contributors to radiation pressure follows:\n" );
00136 PrtLinePres(ioQQQ);
00137 }
00138 lgDone = true;
00139 conv.lgBadStop = true;
00140 strncpy( StopCalc.chReasonStop, "of radiation pressure.", sizeof(StopCalc.chReasonStop) );
00141 }
00142
00143
00144 if( radius.drad_x_fillfac*radius.r1r0sq > BIGFLOAT/10.)
00145 {
00146
00147 lgDone = true;
00148 strncpy( StopCalc.chReasonStop, "volume too large for this cpu.", sizeof(StopCalc.chReasonStop) );
00149 }
00150
00151
00152 else if( !wind.lgVelPos && wind.lgBallistic() )
00153 {
00154
00155 lgDone = true;
00156 strncpy( StopCalc.chReasonStop, "wind veloc too small.", sizeof(StopCalc.chReasonStop) );
00157 }
00158
00159 else if( !wind.lgStatic() && fabs(wind.windv) < StopCalc.StopVelocity )
00160 {
00161
00162 lgDone = true;
00163 strncpy( StopCalc.chReasonStop, "wind V too small.", sizeof(StopCalc.chReasonStop) );
00164 }
00165
00166 else if( (StopCalc.nTotalIonizStop>0) && conv.nTotalIoniz> StopCalc.nTotalIonizStop )
00167 {
00168
00169
00170 lgDone = true;
00171 strncpy( StopCalc.chReasonStop, "nTotalIonizStop reached.", sizeof(StopCalc.chReasonStop) );
00172 }
00173
00174
00175 else if( StopCalc.lgStop21cm && (HFLines[0].Emis().TauCon() >= StopCalc.tauend) )
00176 {
00177 lgDone = true;
00178 strncpy( StopCalc.chReasonStop, "21 cm optical depth.", sizeof(StopCalc.chReasonStop) );
00179 }
00180
00181 else if( rfield.extin_mag_V_extended >= StopCalc.AV_extended )
00182 {
00183
00184 lgDone = true;
00185 strncpy( StopCalc.chReasonStop, "A_V reached.", sizeof(StopCalc.chReasonStop) );
00186 }
00187
00188 else if( rfield.extin_mag_V_point >= StopCalc.AV_point )
00189 {
00190
00191 lgDone = true;
00192 strncpy( StopCalc.chReasonStop, "A_V reached.", sizeof(StopCalc.chReasonStop) );
00193 }
00194
00195 else if( StopCalc.xMass!=0. &&
00196 log10(SDIV(dense.xMassTotal))+1.0992099+ 2.*log10(radius.rinner) >= StopCalc.xMass )
00197 {
00198
00199 lgDone = true;
00200 strncpy( StopCalc.chReasonStop, "mass reached.", sizeof(StopCalc.chReasonStop) );
00201 }
00202
00203
00204
00205
00206 else if( pressure.lgSonicPoint && pressure.lgSonicPointAbortOK )
00207 {
00208
00209 lgDone = true;
00210 strncpy( StopCalc.chReasonStop, "sonic point reached.", sizeof(StopCalc.chReasonStop) );
00211 }
00212
00213 else if( dense.EdenTrue==0 )
00214 {
00215
00216 conv.lgBadStop = true;
00217 lgDone = true;
00218 strncpy( StopCalc.chReasonStop, "zero electron density.", sizeof(StopCalc.chReasonStop) );
00219 }
00220
00221 else if( radius.lgdR2Small )
00222 {
00223 lgDone = true;
00224 conv.lgBadStop = true;
00225 strncpy( StopCalc.chReasonStop, "DR small rel to thick.", sizeof(StopCalc.chReasonStop) );
00226 }
00227
00228 else if( dense.eden < StopCalc.StopElecDensity )
00229 {
00230 lgDone = true;
00231 strncpy( StopCalc.chReasonStop, "lowest EDEN reached.", sizeof(StopCalc.chReasonStop) );
00232 }
00233
00234 else if( dense.eden/dense.gas_phase[ipHYDROGEN] < StopCalc.StopElecFrac )
00235 {
00236 lgDone = true;
00237 strncpy( StopCalc.chReasonStop, "low electron fraction.", sizeof(StopCalc.chReasonStop) );
00238 }
00239
00240
00241
00242
00243 else if( dense.lgElmtOn[ipOXYGEN] &&
00244 (oxy_in_grains/MAX2(SMALLFLOAT,dense.gas_phase[ipOXYGEN])) > StopCalc.StopDepleteFrac )
00245 {
00246 lgDone = true;
00247 strncpy( StopCalc.chReasonStop, "freeze out fraction.", sizeof(StopCalc.chReasonStop) );
00248 }
00249
00250
00251 else if( 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN] > StopCalc.StopH2MoleFrac )
00252 {
00253 lgDone = true;
00254 strncpy( StopCalc.chReasonStop, "large H_2/H fraction.", sizeof(StopCalc.chReasonStop) );
00255 }
00256
00257 else if( dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] <
00258 StopCalc.StopHPlusFrac )
00259 {
00260 lgDone = true;
00261 strncpy( StopCalc.chReasonStop, "low H_+/H fraction.", sizeof(StopCalc.chReasonStop) );
00262 }
00263
00264 else if( radius.lgDrMinUsed )
00265 {
00266
00267 conv.lgBadStop = true;
00268 lgDone = true;
00269 strncpy( StopCalc.chReasonStop, "DRAD small- set DRMIN.", sizeof(StopCalc.chReasonStop) );
00270 }
00271
00272 else if( radius.depth >= radius.StopThickness[iteration-1]/EPS )
00273 {
00274 lgDone = true;
00275 strncpy( StopCalc.chReasonStop, "outer radius reached.", sizeof(StopCalc.chReasonStop) );
00276 }
00277
00278 else if( StopCalc.iptnu >= 0 &&
00279 opac.TauAbsGeo[0][StopCalc.iptnu-1] >= StopCalc.tauend/EPS )
00280 {
00281 lgDone = true;
00282 strncpy( StopCalc.chReasonStop, "optical depth reached.", sizeof(StopCalc.chReasonStop) );
00283 }
00284
00285 else if( StopCalc.lgStopSpeciesColumn &&
00286 findspecieslocal( StopCalc.chSpeciesColumn )->column >= StopCalc.col_species/EPS )
00287 {
00288
00289 lgDone = true;
00290 sprintf( StopCalc.chReasonStop, "%s column dens reached.", StopCalc.chSpeciesColumn );
00291
00292 }
00293
00294 else if( colden.colden[ipCOL_HTOT] >= StopCalc.HColStop/EPS )
00295 {
00296
00297 lgDone = true;
00298 strncpy( StopCalc.chReasonStop, "H column dens reached.", sizeof(StopCalc.chReasonStop) );
00299 }
00300
00301 else if( colden.colden[ipCOL_Hp] >= StopCalc.colpls/EPS )
00302 {
00303 lgDone = true;
00304 strncpy( StopCalc.chReasonStop, "H+ column dens reached.", sizeof(StopCalc.chReasonStop) );
00305 }
00306
00307 else if( (colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]) >= StopCalc.col_h2/EPS )
00308 {
00309
00310 lgDone = true;
00311 strncpy( StopCalc.chReasonStop, "H2 column dens reached.", sizeof(StopCalc.chReasonStop) );
00312 }
00313
00314 else if( (2.*(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]) + colden.colden[ipCOL_H0]) >= StopCalc.col_h2_nut/EPS )
00315 {
00316
00317 lgDone = true;
00318 strncpy( StopCalc.chReasonStop, "H2+H0 column dens reached.", sizeof(StopCalc.chReasonStop) );
00319 }
00320
00321 else if( colden.H0_ov_Tspin >= StopCalc.col_H0_ov_Tspin/EPS )
00322 {
00323
00324 lgDone = true;
00325 strncpy( StopCalc.chReasonStop, "N(H0)/Tspin column dens reached.", sizeof(StopCalc.chReasonStop) );
00326 }
00327
00328 else if( findspecieslocal("CO")->column >= StopCalc.col_monoxco/EPS )
00329 {
00330
00331 lgDone = true;
00332 strncpy( StopCalc.chReasonStop, "CO column dens reached.", sizeof(StopCalc.chReasonStop) );
00333 }
00334
00335 else if( colden.colden[ipCOL_H0] >= StopCalc.colnut/EPS )
00336 {
00337 lgDone = true;
00338 strncpy( StopCalc.chReasonStop, "H0 column dens reached.", sizeof(StopCalc.chReasonStop) );
00339 }
00340
00341
00342
00343 else if(
00344 ( (colden.colden[ipCOL_H0]+2.*(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])) >= StopCalc.colnut/EPS) )
00345 {
00346 lgDone = true;
00347 strncpy( StopCalc.chReasonStop, "H0-H0+H2 column dens reached.", sizeof(StopCalc.chReasonStop) );
00348 }
00349
00350 else if( phycon.te > StopCalc.TempHiStopZone )
00351 {
00352 lgDone = true;
00353 strncpy( StopCalc.chReasonStop, "highest Te reached.", sizeof(StopCalc.chReasonStop) );
00354 }
00355
00356 else if( phycon.te < StopCalc.TempLoStopZone )
00357 {
00358 lgDone = true;
00359 strncpy( StopCalc.chReasonStop, "lowest Te reached.", sizeof(StopCalc.chReasonStop) );
00360 }
00361
00362 else if( nzone >= geometry.nend[iteration-1] )
00363 {
00364 lgDone = true;
00365 geometry.lgZoneTrp = true;
00366 strncpy( StopCalc.chReasonStop, "NZONE reached.", sizeof(StopCalc.chReasonStop) );
00367 }
00368
00369
00370
00371 else if( StopCalc.nstpl > 0 || StopCalc.ContIndex.size() > 0 )
00372 {
00373
00374
00375 for( i=0; i < StopCalc.nstpl; i++ )
00376 {
00377
00378 if( LineSv[StopCalc.ipStopLin2[i]].SumLine[StopCalc.nEmergent[i]] > 0. )
00379 {
00380 char chString[10];
00381 if( LineSv[StopCalc.ipStopLin1[i]].SumLine[StopCalc.nEmergent[i]]/
00382 LineSv[StopCalc.ipStopLin2[i]].SumLine[StopCalc.nEmergent[i]] >
00383 StopCalc.stpint[i] )
00384 {
00385 lgDone = true;
00386 sprt_wl( chString , StopCalc.StopLineWl1[i] );
00387 sprintf( StopCalc.chReasonStop, "line %s %s reached",
00388 StopCalc.chStopLabel1[i] , chString );
00389 }
00390 }
00391 }
00392
00393 for( size_t k=0; k < StopCalc.ContIndex.size(); ++k )
00394 {
00395
00396 long ind = t_PredCont::Inst().offset() + 4*StopCalc.ContIndex[k];
00397 double nFnu_model = LineSv[ind].SumLine[0]*pow(10.,radius.Conv2PrtInten);
00398 if( nFnu_model >= StopCalc.ContNFnu[k].get("erg/s/cm2") )
00399 {
00400 char chTemp[10];
00401 lgDone = true;
00402 sprt_wl( chTemp, LineSv[ind].wavelength );
00403 sprintf( StopCalc.chReasonStop, "flux %s %s reached",
00404 LineSv[ind].chALab, chTemp );
00405 }
00406 }
00407 }
00408
00409 else if( radius.drNext <= 0. )
00410 {
00411
00412 if( called.lgTalk )
00413 {
00414 fprintf( ioQQQ, " drNext=%10.2e STOP\n", radius.drNext );
00415 }
00416 lgDone = true;
00417 conv.lgBadStop = true;
00418 strncpy( StopCalc.chReasonStop, "internal error - DRAD.", sizeof(StopCalc.chReasonStop) );
00419 ShowMe();
00420 cdEXIT(EXIT_FAILURE);
00421 }
00422
00423 else if( lgAbort )
00424 {
00425
00426 conv.lgBadStop = true;
00427 lgDone = true;
00428 strncpy( StopCalc.chReasonStop, "calculation aborted.", sizeof(StopCalc.chReasonStop) );
00429 }
00430
00431 lgPrinted = false;
00432 if( lgDone )
00433 {
00434
00435 lgEndFun_v = true;
00436 PrtZone();
00437 lgPrinted = true;
00438 }
00439
00440 else
00441 {
00442
00443
00444 if( ((nzone/geometry.nprint)*geometry.nprint == nzone ||
00445 nzone == 1) || trace.nTrConvg )
00446 {
00447 PrtZone();
00448 lgPrinted = true;
00449 }
00450
00451 lgEndFun_v = false;
00452 }
00453
00454
00455 if( prt.nzdump == nzone || prt.nzdump == 0 )
00456 dmpary();
00457
00458
00459
00460
00461
00462 if( !hcmap.lgMapDone && (hcmap.MapZone == 0 || nzone == hcmap.MapZone) )
00463 {
00464
00465 if( !lgPrinted )
00466 {
00467 PrtZone();
00468 }
00469
00470
00471 hcmap.lgMapBeingDone = true;
00472
00473
00474
00475 if( ioMAP != NULL )
00476 map_do(ioMAP, " map");
00477 else
00478 map_do(ioQQQ, " map");
00479
00480
00481 lgEndFun_v = true;
00482 strncpy( StopCalc.chReasonStop, "MAP command used-stop.", sizeof(StopCalc.chReasonStop) );
00483
00484
00485
00486 iterations.itermx = 0;
00487
00488
00489 StopCalc.chReasonStop[sizeof(StopCalc.chReasonStop)-1] = '\0';
00490
00491 if( trace.lgTrace )
00492 {
00493 fprintf( ioQQQ, " iter_end_check returns after map.\n" );
00494 }
00495 return( lgEndFun_v );
00496 }
00497
00498 if( lgEndFun_v && prt.lgOnlyZone )
00499 {
00500 cdEXIT(EXIT_FAILURE);
00501 }
00502
00503
00504
00505
00506 if( strlen( StopCalc.chReasonStop ) >= nCHREASONSTOP-1 )
00507 TotalInsanity();
00508
00509 if( trace.lgTrace )
00510 {
00511 fprintf( ioQQQ, " iter_end_check bottom return.\n" );
00512 }
00513 return( lgEndFun_v );
00514 }
00515
00516 #ifdef EPS
00517 # undef EPS
00518 #endif
00519 #define EPS 0.005
00520
00521 STATIC void dmpary(void)
00522 {
00523 long int i;
00524 realnum ratio;
00525
00526 DEBUG_ENTRY( "dmpary()" );
00527
00528 fprintf( ioQQQ,
00529 " This is a print out of the cooling array for zone number %3ld\n",
00530 nzone );
00531
00532 fprintf( ioQQQ,
00533 " The total heating was HTOT=%10.2e erg/s/cm3, the total cooling was CTOT=%10.2e, and the temperature was%10.3eK.\n",
00534 thermal.htot, thermal.ctot, phycon.te );
00535
00536 fprintf( ioQQQ,
00537 " All coolants greater than%6.2f%% of the total will be printed.\n",
00538 EPS*100. );
00539
00540
00541 coolpr(ioQQQ,"ZERO",1,0.,"ZERO");
00542 for( i=0; i < thermal.ncltot; i++ )
00543 {
00544 ratio = (realnum)(thermal.cooling[i]/thermal.ctot);
00545 if( fabs(ratio) > EPS )
00546 {
00547 coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i],
00548 ratio,"DOIT");
00549 }
00550
00551 ratio = (realnum)(thermal.heatnt[i]/thermal.ctot);
00552 if( fabs(ratio) > EPS )
00553 {
00554 coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i],
00555 ratio,"DOIT");
00556 }
00557 }
00558 coolpr(ioQQQ,"DONE",1,0.,"DONE");
00559 return;
00560 }