53 #define DIAG_PRINT false
54 #define MAINPRINT false
147 static double te_old=-1;
148 double timestep_Hp_temp , timestep_return;
170 double dTdStep = fabs(te_new-te_old)/te_new;
172 double dT_factor = 0.04/
SDIV(dTdStep);
173 dT_factor =
MIN2( 2.0 , dT_factor );
174 dT_factor =
MAX2( 0.01 , dT_factor );
175 timestep_Hp_temp =
timestep * dT_factor;
179 timestep_Hp_temp = -1.;
182 if( timestep_Hp_temp > 0. )
183 timestep_return = timestep_Hp_temp;
189 fprintf(
ioQQQ,
"DEBUG timestep_next returns %.3e, old temp %.2e\n" , timestep_return, te_old );
191 return( timestep_return );
224 static double dp = -1.,
230 static int lastzone = -1,
246 " DynaPresChngFactor set PresTotlCorrect=%.3e PresTotlInit=%.3e PresInteg=%.3e DivergePresInteg=%.3e\n",
252 fprintf(stdout,
"Pressure: init %g +rad %g +diverge %g = %g cf %g\n",
303 globalpresmode =
SHOCK;
332 else if(globalpresmode ==
STRONGD)
345 else if(globalpresmode ==
SHOCK)
383 printf(
"Need to set global pressure mode\n");
407 if(fabs(factor-1.) > 0.01)
409 factor = 1.+
sign(0.01,factor-1.);
448 if(fabs(factor-1.) > 3*width)
450 factor = 1.+
sign(3*width,factor-1.);
466 double a, b, c, q, dmin, emin, dsol, f1 , f2;
471 b = (erp-erpp)/(dp-dpp) - a * (dp+dpp);
472 c = erp - dp*(a*dp+b);
474 emin = (a*dmin+b)*dmin+c;
483 fprintf(
ioQQQ,
"Check 2: %g %g\n",erp,(a*dp+b)*dp+c);
484 fprintf(
ioQQQ,
"Check 3: %g %g\n",erpp,(a*dpp+b)*dpp+c);
496 if(globalpresmode ==
STRONGD && -q > 1e-3*b*b)
519 dsol = -(b+sqrt(q))/(2*a);
523 dsol = 2*c/(-b+sqrt(q));
530 dsol = (-b+sqrt(q))/(2*a);
534 dsol = -2*c/(b+sqrt(q));
539 fprintf(
ioQQQ,
"Check 4: %g %g %d %g %g\n",dsol,(a*dsol+b)*dsol+c,iroot,dmin,emin);
548 if(fabs(factor-1.) > 3*width)
550 factor = 1.+
sign(3*width,factor-1.);
569 fprintf(
ioQQQ,
"windv %li r %g v %g f %g\n",
576 " DynaPresChngFactor mode %s scale factor %.3f vel %.3e\n",
582 enum{DEBUG_LOC=
false};
586 char chPRE[][4] = {
"gas" ,
"ram",
"sec",
"par" };
588 "pre %s\tfac\t%.5f\n",
644 for( ion=0; ion<nelem+2; ++ion )
663 fprintf(
ioQQQ,
"workwork\t%li\t%.3e\t%.3e\t%.3e\n",
701 "dynamics cool-heat\t%li\t%.3e\t%.3e\t%.3e\t%.3e\t %.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
754 "PROBLEM conservation error: zn %li elem %li upstream %.8e abund %.8e (up-ab)/up %.2e f1(H n CO) %.2e f2(H n CO) %.2e\n",
789 fprintf(
ioQQQ,
"dynamiccc\t%li\t%.2e\t%.2e\t%.2e\t%.2e\n",
800 fprintf(
ioQQQ,
"dynaionizeeee\t%li\t%i\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
817 fprintf(
ioQQQ,
" DynaIonize, %4li photo=%.2e , H recom= %.2e \n",
839 double upstream, dilution, dilutionleft, dilutionright, frac_next;
842 double hupstream, hnextfrac=-
BIGFLOAT, hion, hmol;
845 double ynextfrac=-
BIGFLOAT, yion, ymol;
847 long int nelem , ion, mol;
858 for( ion=0; ion<nelem+2; ++ion )
911 dilutionright = 1./
Old_hden[ipUpstream+1];
927 fprintf(
ioQQQ,
"PROBLEM advected enthalpy density is not conserved %e\t%e\t%e\t%e\n",
940 for( ion=0; ion<nelem+2; ++ion )
977 if(
COmole[mol]->n_nuclei > 1)
1001 for( ion=0; ion<nelem+2; ++ion )
1019 if(
COmole[mol]->n_nuclei > 1)
1087 for( ion=0; ion<nelem+2; ++ion )
1229 fprintf(
ioQQQ,
" DynaStartZone, %4li photo=%.2e , H recom= %.2e dil %.2e \n",
1244 fprintf(
ioQQQ,
"Check dp: %g %g mom %g %g mas %g\n",
1263 static long int nTime_dt_array_element = 0;
1275 fprintf(
ioQQQ,
"DYNAMICS DynaEndIter sets stop radius to %.2e after"
1276 "dynamics.n_initial_relax=%li iterations.\n",
1317 fprintf(
ioQQQ,
" DynaEndIter, dr=%.2e \n",
1330 static double HeatInitial=-1. , HeatRadiated=-1. ,
1336 "DEBUG times enter timestep %.2e elapsed_time %.2e iteration %li relax %li \n",
1374 " PROBLEM - DynaEnditer - I need the continuum at time %.2e but the table ends at %.2e.\n" ,
1392 fprintf(
ioQQQ,
"DEBUG time dep reset continuum zero timestep %.2e elapsed time %.2e scale %.2e",
1398 fprintf(
ioQQQ,
" recom");
1400 fprintf(
ioQQQ,
"\n");
1416 fprintf(
ioQQQ,
"DEBUG TurbHeat vary new heat %.2e\n",
1423 else if( nTimeVary )
1435 else if( !nTimeVary )
1437 fprintf(
ioQQQ,
" DISASTER - there were no variable continua "
1438 "or heat sources - put TIME option on at least one "
1439 "luminosity or hextra command.\n");
1453 fprintf(
ioQQQ,
"DEBUG relaxing times requested %li this is step %li\n",
1456 fprintf(
ioQQQ,
"DEBUG heat conser HeatInitial=%.2e HeatRadiated=%.2e\n",
1457 HeatInitial , HeatRadiated );
1465 ++nTime_dt_array_element;
1473 fprintf(
ioQQQ,
"DEBUG dynamics turn on recombination logic\n");
1485 fprintf(
ioQQQ,
"DEBUG lgtimes increment Time to %li %.2e\n" ,nTime_dt_array_element,
1497 fprintf(
ioQQQ,
"DEBUG lgtimes increment Timeby timestep_factor to %li %.2e\n" ,
1498 nTime_dt_array_element,
1502 fprintf(
ioQQQ,
"DEBUG lgtimes but reset to %.2e\n" ,
timestep );
1511 fprintf(
ioQQQ,
"DEBUG times exit timestep %.2e elapsed_time %.2e scale %.2e \n",
1553 for(i=0;i<
nzone;++i)
1579 for( ion=0; ion<nelem+2; ++ion )
1641 fprintf(
ioQQQ,
"DYNAMICS DynaNewStep: Dyn_dr %.2e convergence_error %.2e discretization_error %.2e error_scale1 %.2e error_scale2 %.2e\n",
1666 for( i=0; i<
nzone; ++i )
1695 for( ion=0; ion<nelem+2; ++ion )
1834 for( ion=0; ion<nelem+2; ++ion )
1884 for( nelem=0; nelem< (LIMELM+3);++nelem )
1907 for( nelem=0; nelem< (LIMELM+3);++nelem )
1909 for( ion=0; ion<LIMELM+1; ++ion )
2001 fprintf(
ioQQQ,
"WARNING - the log of the initial time step is too "
2002 "large, I shall probably crash. The value was log t = %.2e\n",
2044 " Hit EOF while reading time-continuum list; use END to end list.\n" );
2049 strcpy( chCap, chCard );
2056 while( strncmp(chCap,
"END" ,3 ) != 0 )
2061 " Too many time points have been entered; the limit is %d. Increase variable NTIME in dynamics.c.\n",
2072 " each line must have two numbers, log of time (s0 and flux ratio.\n");
2089 if(
nMatch(
"RECO" , chCap ) &&
nMatch(
"MBIN" , chCap ) )
2107 fprintf(
ioQQQ,
" Hit EOF while reading line list; use END to end list.\n" );
2112 strcpy( chCap, chCard );
2118 fprintf(
ioQQQ,
"DEBUG time dep %.2e %.2e %.2e %.2e\n",
2124 fprintf(
ioQQQ,
"\n" );
2147 if( (i =
nMatch(
"VELO",chCard))>0 )
2155 if( (i =
nMatch(
"DFDR",chCard))>0 )
2164 if( (i =
nMatch(
"CENT",chCard))>0 )
2172 if( (i =
nMatch(
"INDE",chCard))>0 )
2180 if(iVelocity_Type == 1)
2208 else if(iVelocity_Type == 2)
2212 fprintf(
ioQQQ,
"Can't specify gradient when flux is constant!\n");
2255 fprintf(
ioQQQ,
"Scale %g (*%c) Index %g Center %g\n",
2261 if(
nMatch(
"ADVE" , chCard ) )
2273 fprintf(
ioQQQ,
" >>>>Initial wind velocity should be greater than speed of sound; calculation only valid above sonic point.\n" );
2285 if(
nMatch(
"DISK" , chCard ) )
2294 if(
nMatch(
"O CO",chCard) )
2315 fprintf(
ioQQQ,
" DYNAMICS Advection: Uad %.2f Uwd%.2e FRCcool: %4.2f Heat %4.2f\n",
2324 fprintf(
ioQQQ,
" DYNAMICS Eexcit:%.4e Eion:%.4e Ebin:%.4e Ekin:%.4e ET+pdv %.4e EnthalpyDensity/rho%.4e AdvSpWork%.4e\n",
2341 if( strcmp( chJob ,
"END" ) == 0 )
2423 "%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\n" ,
2449 fprintf( ipPnunit ,
"%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",