00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "cddefines.h"
00010
00011 #include "pressure_change.h"
00012
00013 #include "colden.h"
00014 #include "conv.h"
00015 #include "cosmology.h"
00016 #include "dark_matter.h"
00017 #include "dense.h"
00018 #include "dynamics.h"
00019 #include "geometry.h"
00020 #include "mole.h"
00021 #include "phycon.h"
00022 #include "pressure.h"
00023 #include "radius.h"
00024 #include "struc.h"
00025 #include "thermal.h"
00026 #include "trace.h"
00027 #include "wind.h"
00028
00029
00030 double zoneDensity()
00031 {
00032 double new_density;
00033
00034 DEBUG_ENTRY( "zoneDensity()" );
00035
00036 pressure.PresTotlError = 0.;
00037
00038
00039
00040 if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
00041 {
00042
00043 if( radius.glbdst < 0. )
00044 {
00045 fprintf( ioQQQ, " Globule distance is negative, internal overflow has occured, sorry.\n" );
00046 fprintf( ioQQQ, " This is routine lgConvPres, GLBDST is%10.2e\n",
00047 radius.glbdst );
00048 cdEXIT(EXIT_FAILURE);
00049 }
00050 new_density =
00051 (radius.glbden*pow(radius.glbrad/(radius.glbdst),radius.glbpow))*
00052 (scalingDensity()/dense.gas_phase[ipHYDROGEN]);
00053 }
00054 else if( cosmology.lgDo )
00055 {
00056
00057 double dnew = GetDensity( cosmology.redshift_current );
00058 new_density = dnew*(scalingDensity()/dense.gas_phase[ipHYDROGEN]);
00059 }
00060 else if( (strcmp(dense.chDenseLaw,"WIND") == 0) ) {
00061
00062 if ( wind.lgStatic() )
00063 {
00064 fprintf( ioQQQ, " PROBLEM WIND called with zero velocity - this is impossible.\n Sorry.\n" );
00065 TotalInsanity();
00066 }
00067
00068
00069 else if( wind.lgBallistic() )
00070 {
00071
00072
00073
00074 if( nzone > 1 )
00075 {
00076
00077
00078 if( trace.lgTrace && trace.lgWind )
00079 {
00080 fprintf(ioQQQ," lgConvPres sets AccelGravity %.3e lgDisk?%c\n",
00081 wind.AccelGravity ,
00082 TorF(wind.lgDisk) );
00083 }
00084
00085
00086
00087
00088
00089 double term = POW2(struc.windv[nzone-2]) + 2.*(wind.AccelTotalOutward - wind.AccelGravity)* radius.drad;
00090
00091
00092 fixit();
00093 if( term <= 1e3 )
00094 {
00095
00096
00097 wind.lgVelPos = false;
00098 }
00099 else
00100 {
00101
00102 term = sqrt(term);
00103 if (wind.windv > 0)
00104 {
00105 wind.windv = (realnum) term;
00106 }
00107 else
00108 {
00109 wind.windv = -(realnum) term;
00110 }
00111 wind.lgVelPos = true;
00112 }
00113
00114 if( trace.lgTrace && trace.lgWind )
00115 {
00116 fprintf(ioQQQ," lgConvPres new wind V zn%li %.3e AccelTotalOutward %.3e AccelGravity %.3e\n",
00117 nzone,wind.windv, wind.AccelTotalOutward, wind.AccelGravity );
00118 }
00119 }
00120
00121
00122 new_density = wind.emdot/(wind.windv*radius.r1r0sq) *
00123 (scalingDensity()/dense.gas_phase[ipHYDROGEN]);
00124 }
00125 else
00126 {
00127 fprintf(ioQQQ,"chDenseLaw==\"WIND\" must now be ballistic or static\n");
00128 TotalInsanity();
00129 }
00130 }
00131
00132 else if( strcmp(dense.chDenseLaw,"SINE") == 0 )
00133 {
00134
00135 if( dense.lgDenFlucRadius )
00136 {
00137 new_density =
00138 (dense.cfirst*cos(radius.depth*dense.flong+dense.flcPhase) +
00139 dense.csecnd)*(scalingDensity()/dense.gas_phase[ipHYDROGEN]);
00140 }
00141 else
00142 {
00143 new_density =
00144 (dense.cfirst*cos(colden.colden[ipCOL_HTOT]*dense.flong+dense.flcPhase) +
00145 dense.csecnd)*(scalingDensity()/dense.gas_phase[ipHYDROGEN]);
00146 }
00147 }
00148
00149 else if( strcmp(dense.chDenseLaw,"POWR") == 0 )
00150 {
00151
00152 double dnew = dense.den0*pow(radius.Radius/radius.rinner,(double)dense.DensityPower);
00153 new_density = dnew*(scalingDensity()/dense.gas_phase[ipHYDROGEN]);
00154 }
00155
00156 else if( strcmp(dense.chDenseLaw,"POWD") == 0 )
00157 {
00158
00159 double dnew = dense.den0*pow(1. + radius.depth/dense.rscale,(double)dense.DensityPower);
00160 new_density = dnew*(scalingDensity()/dense.gas_phase[ipHYDROGEN]);
00161 }
00162
00163 else if( strcmp(dense.chDenseLaw,"POWC") == 0 )
00164 {
00165
00166 double dnew = dense.den0*pow(1.f + colden.colden[ipCOL_HTOT]/
00167 dense.rscale,dense.DensityPower);
00168 new_density = dnew*(scalingDensity()/dense.gas_phase[ipHYDROGEN]);
00169 }
00170
00171
00172 else if( strncmp( dense.chDenseLaw ,"DLW" , 3) == 0 )
00173 {
00174 if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
00175 {
00176
00177 new_density =
00178 dense_fabden(radius.Radius,radius.depth)*
00179 (scalingDensity()/dense.gas_phase[ipHYDROGEN]);
00180 }
00181 else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
00182 {
00183
00184
00185 new_density = dense_tabden(radius.Radius,radius.depth) *
00186 (scalingDensity()/dense.gas_phase[ipHYDROGEN]);
00187 }
00188 else if( strcmp(dense.chDenseLaw,"DLW3") == 0 )
00189 {
00190
00191 new_density = dense_parametric_wind(radius.Radius) *
00192 (scalingDensity()/dense.gas_phase[ipHYDROGEN]);
00193 }
00194 else
00195 {
00196 fprintf( ioQQQ, " Insanity, lgConvPres gets chCPres=%4.4s\n",
00197 dense.chDenseLaw );
00198 cdEXIT(EXIT_FAILURE);
00199 }
00200
00201 if( new_density/scalingDensity() > 3. || new_density/scalingDensity()< 1./3 )
00202 {
00203 static bool lgWARN2BIG=false;
00204 if( !lgWARN2BIG )
00205 {
00206 lgWARN2BIG = true;
00207 fprintf(ioQQQ,"\n\n >========== Warning! The tabulated or functional change in density as a function of depth was VERY large. This is zone %li.\n",nzone);
00208 fprintf(ioQQQ," >========== Warning! This will cause convergence problems. \n");
00209 fprintf(ioQQQ," >========== Warning! The current radius is %.3e. \n",radius.Radius);
00210 fprintf(ioQQQ," >========== Warning! The current depth is %.3e. \n",radius.depth);
00211 fprintf(ioQQQ," >========== Warning! Consider using a more modest change in density vs radius. \n\n\n");
00212 }
00213 }
00214 }
00215
00216 else if( strcmp(dense.chDenseLaw,"CDEN") == 0 )
00217 {
00218
00219 new_density = scalingDensity();
00220 }
00221
00222 else
00223 {
00224 fprintf( ioQQQ, " Unknown density law=%s= This is a major internal error.\n",
00225 dense.chDenseLaw );
00226 ShowMe();
00227 cdEXIT(EXIT_FAILURE);
00228 }
00229
00230 return new_density;
00231 }
00232
00233 enum {
00234 CPRE,
00235 SUBSONIC,
00236 SUPERSONIC,
00237 STRONGD,
00238 ORIGINAL,
00239 SHOCK,
00240 ANTISHOCK,
00241 ANTISHOCK2
00242 };
00243
00245 STATIC double stepDensity(const PresMode & presmode, solverState &st);
00246
00247 STATIC void logPressureState()
00248 {
00249
00250 conv.hist_pres_density.push_back( dense.gas_phase[ipHYDROGEN] );
00251 conv.hist_pres_current.push_back( pressure.PresTotlCurr );
00252 conv.hist_pres_error.push_back( pressure.PresTotlError );
00253 }
00254
00255 STATIC bool lgTestPressureConvergence( double new_density)
00256 {
00257 double density_change_factor = new_density/scalingDensity();
00258
00259
00260 if( density_change_factor > 1. + conv.PressureErrorAllowed ||
00261 density_change_factor < 1. - conv.PressureErrorAllowed )
00262 {
00263 return false;
00264 }
00265 else
00266 {
00267 return true;
00268 }
00269 }
00270
00271 STATIC double limitedDensityScaling( double new_density, double dP_chng_factor )
00272 {
00273 double density_change_factor = new_density/scalingDensity();
00274
00275
00276 double pdelta = conv.MaxFractionalDensityStepPerIteration * dP_chng_factor;
00277
00278
00279 density_change_factor = MIN2(density_change_factor,1.+pdelta);
00280 density_change_factor = MAX2(density_change_factor,1.-pdelta);
00281 return density_change_factor;
00282 }
00283
00284
00285
00286
00287 bool PressureChange(
00288
00289 double dP_chng_factor, const PresMode & presmode, solverState &st )
00290 {
00291 DEBUG_ENTRY( "PressureChange()" );
00292
00293
00294
00295
00296 PresTotCurrent();
00297 logPressureState();
00298
00299 double new_density = stepDensity(presmode, st);
00300
00301 conv.lgConvPres = lgTestPressureConvergence(new_density);
00302
00303 double density_change_factor = limitedDensityScaling(new_density, dP_chng_factor);
00304
00305 {
00306
00307 enum{DEBUG_LOC=false};
00308 static long int nsave=-1;
00309
00310 if( DEBUG_LOC )
00311 {
00312 if( nsave-nzone ) fprintf(ioQQQ,"\n");
00313 nsave = nzone;
00314 fprintf(ioQQQ,"nnzzone\t%li\t%.2f%%\t%.3f\t%.2e\t%.2e\t%.2e\n",
00315 nzone,
00316 density_change_factor,
00317
00318 pressure.PresTotlError*100.,
00319 pressure.PresTotlCurr,
00320 pressure.PresGasCurr,
00321 pressure.pres_radiation_lines_curr );
00322 }
00323 }
00324
00325 if( trace.lgTrace )
00326 {
00327 fprintf( ioQQQ,
00328 " PressureChange called, changing HDEN from %10.3e to %10.3e Set fill fac to %10.3e\n",
00329 dense.gas_phase[ipHYDROGEN], density_change_factor*dense.gas_phase[ipHYDROGEN],
00330 geometry.FillFac );
00331 }
00332
00333 bool lgChange = ( density_change_factor != 1. );
00334
00335 if( lgChange )
00336 {
00337 ScaleAllDensities((realnum) density_change_factor);
00338
00339
00340
00341 TempChange(phycon.te , false);
00342 }
00343
00344 {
00345
00346 enum {DEBUG_LOC=false};
00347
00348 if( DEBUG_LOC && (nzone>215) )
00349 {
00350 fprintf( ioQQQ,
00351 "%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%c\n",
00352 radius.depth,
00353 pressure.PresTotlCurr,
00354 pressure.PresTotlInit + pressure.PresInteg,
00355 pressure.PresTotlInit,
00356 pressure.PresGasCurr,
00357 pressure.PresRamCurr,
00358 pressure.pres_radiation_lines_curr,
00359
00360 pressure.PresInteg - pressure.pinzon,
00361 wind.windv/1e5,
00362 sqrt(5.*pressure.PresGasCurr/3./dense.xMassDensity)/1e5,
00363 TorF(conv.lgConvPres) );
00364 }
00365 }
00366
00367 return lgChange;
00368 }
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388 STATIC double stepDensity(const PresMode &presmode, solverState &st)
00389 {
00390 DEBUG_ENTRY( "stepDensity()" );
00391
00392 double PresTotlCorrect = st.press;
00393
00394 double densityCurrent = scalingDensity();
00395
00396 double er = pressure.PresTotlCurr-PresTotlCorrect;
00397 pressure.PresTotlError = er/PresTotlCorrect;
00398
00399 if( trace.lgTrace && presmode.global != CPRE )
00400 {
00401 fprintf( ioQQQ,
00402 " DynaPresChngFactor set PresTotlCorrect=%.3e PresTotlInit=%.3e PresInteg=%.3e DivergePresInteg=%.3e\n",
00403 PresTotlCorrect , pressure.PresTotlInit , pressure.PresInteg*pressure.lgContRadPresOn,
00404 dynamics.DivergePresInteg );
00405 }
00406
00407 if( dynamics.lgTracePrint && presmode.global != CPRE)
00408 fprintf(ioQQQ,"Pressure: init %g +rad %g +diverge %g = %g cf %g\n",
00409 pressure.PresTotlInit, pressure.PresInteg*pressure.lgContRadPresOn,
00410 dynamics.DivergePresInteg, PresTotlCorrect, pressure.PresTotlCurr);
00411
00412
00413
00414 double rhohat;
00415 if( presmode.global == CPRE || presmode.global == ORIGINAL ||
00416 st.lastzone != nzone || fabs(er-st.erp) < SMALLFLOAT
00417 || fabs(densityCurrent-st.dp) < SMALLFLOAT)
00418 {
00419 double dpdrho;
00420 if ( presmode.zone == SUBSONIC )
00421 {
00422 dpdrho = pressure.PresTotlCurr/densityCurrent;
00423 }
00424 else
00425 {
00426 dpdrho = -pressure.PresTotlCurr/densityCurrent;
00427 }
00428 rhohat = densityCurrent - er/dpdrho;
00429 }
00430 else
00431 {
00432
00433
00434 double dpdrho = (st.erp-er)/(st.dp-densityCurrent);
00435 rhohat = densityCurrent - er/dpdrho;
00436
00437
00438
00439
00440
00441
00442 if(presmode.zone == SUBSONIC && dpdrho <= 0)
00443 {
00444 rhohat = 1.03*densityCurrent;
00445 }
00446 else if(presmode.zone == SUPERSONIC && dpdrho >= 0)
00447 {
00448 rhohat = 0.97*densityCurrent;
00449 }
00450 }
00451
00452 st.dp = densityCurrent;
00453 st.erp = er;
00454 st.lastzone = nzone;
00455
00456 if( presmode.global != CPRE && dynamics.lgTracePrint )
00457 fprintf(ioQQQ,"windv %li r %g v %g f %g\n",
00458 nzone,radius.depth,wind.windv,DynaFlux(radius.depth));
00459
00460
00461 if( presmode.global != CPRE && trace.nTrConvg>=1 )
00462 {
00463 fprintf( ioQQQ,
00464 " DynaPresChngFactor mode %s new density %.3f vel %.3e\n",
00465 dynamics.chPresMode , rhohat , wind.windv );
00466 }
00467
00468 return rhohat;
00469 }
00470
00471 void PresMode::set()
00472 {
00473 DEBUG_ENTRY("PresMode::set()");
00474
00475
00476
00477 if( !dynamics.lgSetPresMode )
00478 {
00479
00480
00481
00482 if(pressure.PresGasCurr < pressure.PresRamCurr)
00483 strcpy( dynamics.chPresMode , "supersonic" );
00484 else
00485 strcpy( dynamics.chPresMode , "subsonic" );
00486
00487 dynamics.lgSetPresMode = true;
00488 }
00489
00490
00491
00492
00493
00494
00495
00496 if (strcmp(dense.chDenseLaw,"CPRE") != 0)
00497 {
00498 ASSERT (strcmp(dense.chDenseLaw,"DYNA") == 0);
00499 }
00500
00501 if (strcmp(dense.chDenseLaw,"CPRE") == 0)
00502 {
00503 global = CPRE;
00504 pressure.lgSonicPointAbortOK = false;
00505 }
00506 else if( strcmp( dynamics.chPresMode , "original" ) == 0 )
00507 {
00508 global = ORIGINAL;
00509 pressure.lgSonicPointAbortOK = true;
00510 }
00511 else if( strcmp( dynamics.chPresMode , "subsonic" ) == 0 )
00512 {
00513 global = SUBSONIC;
00514 pressure.lgSonicPointAbortOK = true;
00515 }
00516
00517 else if( strcmp( dynamics.chPresMode , "supersonic" ) == 0 )
00518 {
00519 global = SUPERSONIC;
00520 pressure.lgSonicPointAbortOK = true;
00521 }
00522
00523 else if( strcmp( dynamics.chPresMode , "strongd" ) == 0 )
00524 {
00525 global = STRONGD;
00526 pressure.lgSonicPointAbortOK = false;
00527 }
00528 else if( strcmp( dynamics.chPresMode , "shock" ) == 0 )
00529 {
00530 global = SHOCK;
00531 pressure.lgSonicPointAbortOK = false;
00532 }
00533 else if( strcmp( dynamics.chPresMode , "antishock-by-mach" ) == 0 )
00534 {
00535 global = ANTISHOCK2;
00536 pressure.lgSonicPointAbortOK = false;
00537 }
00538 else if( strcmp( dynamics.chPresMode , "antishock" ) == 0 )
00539 {
00540 global = ANTISHOCK;
00541 pressure.lgSonicPointAbortOK = false;
00542 }
00543
00544
00545
00546 if(global == CPRE)
00547 {
00548 zone = SUBSONIC;
00549 }
00550 else if(global == ORIGINAL)
00551 {
00552
00553
00554 if(pressure.PresGasCurr > pressure.PresRamCurr)
00555 {
00556 zone = SUBSONIC;
00557 }
00558 else
00559 {
00560 zone = SUPERSONIC;
00561 }
00562 }
00563 else if(global == STRONGD)
00564 {
00565
00566 zone = SUPERSONIC;
00567 }
00568 else if(global == SUBSONIC)
00569 {
00570 zone = SUBSONIC;
00571 }
00572 else if(global == SUPERSONIC)
00573 {
00574 zone = SUPERSONIC;
00575 }
00576 else if(global == SHOCK)
00577 {
00578 if(radius.depth < dynamics.ShockDepth)
00579 {
00580 zone = SUBSONIC;
00581 }
00582 else
00583 {
00584 zone = SUPERSONIC;
00585 }
00586 }
00587 else if(global == ANTISHOCK)
00588 {
00589 if(radius.depth < dynamics.ShockDepth)
00590 {
00591 zone = SUPERSONIC;
00592 }
00593 else
00594 {
00595 zone = SUBSONIC;
00596 }
00597 }
00598 else if(global == ANTISHOCK2)
00599 {
00600
00601
00602
00603 if( fabs(wind.windv) > dynamics.ShockMach*sqrt(pressure.PresGasCurr/dense.xMassDensity))
00604 {
00605 zone = SUPERSONIC;
00606 }
00607 else
00608 {
00609 zone = SUBSONIC;
00610 }
00611 }
00612 else
00613 {
00614 printf("Need to set global pressure mode\n");
00615 cdEXIT( EXIT_FAILURE );
00616 }
00617 }
00618
00619 double pressureZone(const PresMode &presmode)
00620 {
00621 double press;
00622 if( presmode.global == CPRE )
00623 {
00624
00625 if( pressure.lgContRadPresOn )
00626 {
00627
00628
00629 press = pressure.PresTotlInit + pressure.PresInteg;
00630 }
00631 else
00632 {
00633
00634 press = pressure.PresTotlInit*
00635
00636 pow(radius.Radius/radius.rinner,(double)pressure.PresPowerlaw);
00637 }
00638
00639 if( dark.lgNFW_Set || pressure.gravity_symmetry >=0 )
00640 {
00641 fixit();
00642
00643
00644
00645
00646 press += pressure.IntegRhoGravity;
00647 }
00648 }
00649 else
00650 {
00651
00652 press = pressure.PresTotlInit + pressure.PresInteg*pressure.lgContRadPresOn
00653 + dynamics.DivergePresInteg;
00654 }
00655 return press;
00656 }
00657