40 if( strcmp(
dense.chDenseLaw,
"GLOB") == 0 )
45 fprintf(
ioQQQ,
" Globule distance is negative, internal overflow has occured, sorry.\n" );
46 fprintf(
ioQQQ,
" This is routine lgConvPres, GLBDST is%10.2e\n",
60 else if( (strcmp(
dense.chDenseLaw,
"WIND") == 0) ) {
62 if (
wind.lgStatic() )
64 fprintf(
ioQQQ,
" PROBLEM WIND called with zero velocity - this is impossible.\n Sorry.\n" );
69 else if(
wind.lgBallistic() )
80 fprintf(
ioQQQ,
" lgConvPres sets AccelGravity %.3e lgDisk?%c\n",
97 wind.lgVelPos =
false;
111 wind.lgVelPos =
true;
116 fprintf(
ioQQQ,
" lgConvPres new wind V zn%li %.3e AccelTotalOutward %.3e AccelGravity %.3e\n",
127 fprintf(
ioQQQ,
"chDenseLaw==\"WIND\" must now be ballistic or static\n");
132 else if( strcmp(
dense.chDenseLaw,
"SINE") == 0 )
135 if(
dense.lgDenFlucRadius )
149 else if( strcmp(
dense.chDenseLaw,
"POWR") == 0 )
156 else if( strcmp(
dense.chDenseLaw,
"POWD") == 0 )
163 else if( strcmp(
dense.chDenseLaw,
"POWC") == 0 )
172 else if( strncmp(
dense.chDenseLaw ,
"DLW" , 3) == 0 )
174 if( strcmp(
dense.chDenseLaw,
"DLW1") == 0 )
181 else if( strcmp(
dense.chDenseLaw,
"DLW2") == 0 )
188 else if( strcmp(
dense.chDenseLaw,
"DLW3") == 0 )
196 fprintf(
ioQQQ,
" Insanity, lgConvPres gets chCPres=%4.4s\n",
203 static bool lgWARN2BIG=
false;
207 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);
208 fprintf(
ioQQQ,
" >========== Warning! This will cause convergence problems. \n");
209 fprintf(
ioQQQ,
" >========== Warning! The current radius is %.3e. \n",
radius.Radius);
210 fprintf(
ioQQQ,
" >========== Warning! The current depth is %.3e. \n",
radius.depth);
211 fprintf(
ioQQQ,
" >========== Warning! Consider using a more modest change in density vs radius. \n\n\n");
216 else if( strcmp(
dense.chDenseLaw,
"CDEN") == 0 )
224 fprintf(
ioQQQ,
" Unknown density law=%s= This is a major internal error.\n",
251 conv.hist_pres_current.push_back(
pressure.PresTotlCurr );
252 conv.hist_pres_error.push_back(
pressure.PresTotlError );
260 if( density_change_factor > 1. +
conv.PressureErrorAllowed ||
261 density_change_factor < 1. -
conv.PressureErrorAllowed )
276 double pdelta =
conv.MaxFractionalDensityStepPerIteration * dP_chng_factor;
279 density_change_factor =
MIN2(density_change_factor,1.+pdelta);
280 density_change_factor =
MAX2(density_change_factor,1.-pdelta);
281 return density_change_factor;
307 enum{DEBUG_LOC=
false};
308 static long int nsave=-1;
314 fprintf(
ioQQQ,
"nnzzone\t%li\t%.2f%%\t%.3f\t%.2e\t%.2e\t%.2e\n",
316 density_change_factor,
321 pressure.pres_radiation_lines_curr );
328 " PressureChange called, changing HDEN from %10.3e to %10.3e Set fill fac to %10.3e\n",
333 bool lgChange = ( density_change_factor != 1. );
346 enum {DEBUG_LOC=
false};
348 if( DEBUG_LOC && (
nzone>215) )
351 "%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%c\n",
392 double PresTotlCorrect = st.
press;
396 double er =
pressure.PresTotlCurr-PresTotlCorrect;
397 pressure.PresTotlError = er/PresTotlCorrect;
402 " DynaPresChngFactor set PresTotlCorrect=%.3e PresTotlInit=%.3e PresInteg=%.3e DivergePresInteg=%.3e\n",
408 fprintf(
ioQQQ,
"Pressure: init %g +rad %g +diverge %g = %g cf %g\n",
422 dpdrho =
pressure.PresTotlCurr/densityCurrent;
426 dpdrho = -
pressure.PresTotlCurr/densityCurrent;
428 rhohat = densityCurrent - er/dpdrho;
434 double dpdrho = (st.
erp-er)/(st.
dp-densityCurrent);
435 rhohat = densityCurrent - er/dpdrho;
444 rhohat = 1.03*densityCurrent;
448 rhohat = 0.97*densityCurrent;
452 st.
dp = densityCurrent;
457 fprintf(
ioQQQ,
"windv %li r %g v %g f %g\n",
464 " DynaPresChngFactor mode %s new density %.3f vel %.3e\n",
483 strcpy(
dynamics.chPresMode ,
"supersonic" );
485 strcpy(
dynamics.chPresMode ,
"subsonic" );
496 if (strcmp(
dense.chDenseLaw,
"CPRE") != 0)
501 if (strcmp(
dense.chDenseLaw,
"CPRE") == 0)
504 pressure.lgSonicPointAbortOK =
false;
506 else if( strcmp(
dynamics.chPresMode ,
"original" ) == 0 )
509 pressure.lgSonicPointAbortOK =
true;
511 else if( strcmp(
dynamics.chPresMode ,
"subsonic" ) == 0 )
514 pressure.lgSonicPointAbortOK =
true;
517 else if( strcmp(
dynamics.chPresMode ,
"supersonic" ) == 0 )
520 pressure.lgSonicPointAbortOK =
true;
523 else if( strcmp(
dynamics.chPresMode ,
"strongd" ) == 0 )
526 pressure.lgSonicPointAbortOK =
false;
528 else if( strcmp(
dynamics.chPresMode ,
"shock" ) == 0 )
531 pressure.lgSonicPointAbortOK =
false;
533 else if( strcmp(
dynamics.chPresMode ,
"antishock-by-mach" ) == 0 )
536 pressure.lgSonicPointAbortOK =
false;
538 else if( strcmp(
dynamics.chPresMode ,
"antishock" ) == 0 )
541 pressure.lgSonicPointAbortOK =
false;
614 printf(
"Need to set global pressure mode\n");
NORETURN void TotalInsanity(void)
#define DEBUG_ENTRY(funcname)
realnum GetDensity(realnum z)
void ScaleAllDensities(realnum factor)
realnum scalingDensity(void)
double dense_tabden(double r0, double depth)
double dense_fabden(double radius, double depth)
double dense_parametric_wind(double rad)
realnum DynaFlux(double depth)
void PresTotCurrent(void)
STATIC bool lgTestPressureConvergence(double new_density)
STATIC void logPressureState()
double pressureZone(const PresMode &presmode)
STATIC double stepDensity(const PresMode &presmode, solverState &st)
STATIC double limitedDensityScaling(double new_density, double dP_chng_factor)
bool PressureChange(double dP_chng_factor, const PresMode &presmode, solverState &st)
void TempChange(double TempNew, bool lgForceUpdate)