136 static double te_old=-1;
137 double timestep_Hp_temp , timestep_return;
141 timestep_return =
dynamics.timestep;
159 double dTdStep = fabs(te_new-te_old)/te_new;
161 double dT_factor = 0.04/
SDIV(dTdStep);
162 dT_factor =
MIN2( 2.0 , dT_factor );
163 dT_factor =
MAX2( 0.01 , dT_factor );
164 timestep_Hp_temp =
dynamics.timestep * dT_factor;
168 timestep_Hp_temp = -1.;
171 if( timestep_Hp_temp > 0. )
172 timestep_return = timestep_Hp_temp;
176 timestep_return =
dynamics.timestep_init;
178 fprintf(
ioQQQ,
"DEBUG timestep_next returns %.3e, old temp %.2e\n" , timestep_return, te_old );
180 return( timestep_return );
192 if( !
dynamics.lgTimeDependentStatic )
211 double depth =
radius.depth;
213 ( !
dynamics.lgTimeDependentStatic &&
214 ( depth < 0 || depth >
dynamics.oldFullDepth ) ) )
225 for(
long ion=0; ion<nelem+2; ++ion )
232 for(
long nelem=ipISO; nelem<
LIMELM; ++nelem)
234 if(
dense.lgElmtOn[nelem] )
236 for(
long level=0; level <
iso_sp[ipISO][nelem].numLevels_local; ++level )
238 dynamics.StatesElem[nelem][nelem-ipISO][level] = 0.;
252 fprintf(
ioQQQ,
"workwork\t%li\t%.3e\t%.3e\t%.3e\n",
272 "dynamics cool-heat\t%li\t%.3e\t%i\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
308 if(
dense.lgElmtOn[nelem] )
314 "PROBLEM conservation error: zn %li elem %li upstream %.8e abund %.8e (up-ab)/up %.2e\n",
318 dense.gas_phase[nelem] ,
323 for(
long ion=0; ion<
dense.IonLow[nelem]; ++ion )
327 for(
long ion=
dense.IonLow[nelem]; ion<=
dense.IonHigh[nelem]; ++ion )
340 for(
long ion=
dense.IonHigh[nelem]+1;ion<nelem+2; ++ion )
351 for(
long nelem=ipISO; nelem<
LIMELM; ++nelem)
353 if(
dense.lgElmtOn[nelem] )
355 for(
long level=0; level <
iso_sp[ipISO][nelem].numLevels_local; ++level )
357 dynamics.StatesElem[nelem][nelem-ipISO][level] =
365 fprintf(
ioQQQ,
"dynamiccc\t%li\t%.2e\t%.2e\t%.2e\t%.2e\n",
376 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",
381 dense.xIonDense[nelem][ion],
384 dense.xIonDense[nelem][ion+1],
393 fprintf(
ioQQQ,
" DynaIonize, %4li photo=%.2e , H recom= %.2e \n",
415 double upstream, dilution, dilutionleft, dilutionright, frac_next;
418 double hupstream, hnextfrac=-
BIGFLOAT, hion, hmol, hiso;
421 double ynextfrac=-
BIGFLOAT, yion, ymol, yiso;
423 long int nelem , ion, mol;
434 for( ion=0; ion<nelem+2; ++ion )
441 for( nelem=ipISO; nelem<
LIMELM; ++nelem)
443 if(
dense.lgElmtOn[nelem] )
445 for(
long level=0; level <
iso_sp[ipISO][nelem].numLevels_max; ++level )
474 hupstream = 0.5*(
radius.depth + upstream);
501 dilution = 1./
dynamics.Upstream_density;
516 fprintf(
ioQQQ,
"PROBLEM interpolated enthalpy density is not within range %.16e\t%.16e\t%.16e\t%e\t%e\n",
517 lo, x, hi, (hi-x)/(hi-lo), (x-lo)/(hi-lo));
525 for( ion=0; ion<nelem+2; ++ion )
546 for( nelem=ipISO; nelem<
LIMELM; ++nelem)
548 if(
dense.lgElmtOn[nelem] )
550 for(
long level=0; level <
iso_sp[ipISO][nelem].numLevels_max; ++level )
571 if(
mole.species[mol].location == NULL &&
mole_global.list[mol]->parentLabel.empty())
573 for(molecule::nAtomsMap::iterator atom=
mole_global.list[mol]->nAtom.begin();
574 atom !=
mole_global.list[mol]->nAtom.end(); ++atom)
590 dilution = 1./
dynamics.Upstream_density;
596 for( ion=0; ion<nelem+2; ++ion )
607 for( nelem=ipISO; nelem<
LIMELM; ++nelem)
609 if(
dense.lgElmtOn[nelem] )
611 for(
long level=0; level <
iso_sp[ipISO][nelem].numLevels_max; ++level )
625 if(
mole.species[mol].location == NULL &&
mole_global.list[mol]->parentLabel.empty())
627 for(molecule::nAtomsMap::iterator atom=
mole_global.list[mol]->nAtom.begin();
628 atom !=
mole_global.list[mol]->nAtom.end(); ++atom)
675 const double STEP_FACTOR=0.05;
679 for( ion=0; ion<nelem+2; ++ion )
734 for( nelem=ipISO; nelem<
LIMELM; ++nelem)
736 if(
dense.lgElmtOn[nelem] )
738 for(
long level=0; level <
iso_sp[ipISO][nelem].numLevels_max; ++level )
846 fprintf(
ioQQQ,
" DynaStartZone, %4li photo=%.2e , H recom= %.2e dil %.2e \n",
862 fprintf(
ioQQQ,
"Check dp: %g %g mom %g %g mas %g\n",
881 static long int nTime_dt_array_element = 0;
893 fprintf(
ioQQQ,
"DYNAMICS DynaIterEnd sets stop radius to %.2e after "
894 "dynamics.n_initial_relax=%li iterations.\n",
912 if( !
dynamics.lgTimeDependentStatic )
938 fprintf(
ioQQQ,
" DynaIterEnd, dr=%.2e \n",
951 static double HeatInitial=-1. , HeatRadiated=-1. ,
957 "DEBUG times enter dynamics.timestep %.2e elapsed_time %.2e iteration %li relax %li \n",
970 thermal.lgTemperatureConstant =
false;
997 HeatInitial = 1.5 *
pressure.PresGasCurr;
1000 fprintf(
ioQQQ,
"DEBUG relaxing times requested %li this is step %li\n",
1003 fprintf(
ioQQQ,
"DEBUG heat conser HeatInitial=%.2e HeatRadiated=%.2e\n",
1004 HeatInitial , HeatRadiated );
1009 fprintf(
ioQQQ,
"DEBUG lgtimes -- stop time reached.\n" );
1017 ++nTime_dt_array_element;
1025 fprintf(
ioQQQ,
"DEBUG dynamics turn on recombination logic\n");
1032 radius.lgSdrmaxRel =
false;
1038 fprintf(
ioQQQ,
"DEBUG lgtimes increment Time to %li %.2e\n" ,nTime_dt_array_element,
1050 fprintf(
ioQQQ,
"DEBUG lgtimes increment Timeby dynamics.timestep_factor to %li %.2e\n" ,
1051 nTime_dt_array_element,
1055 fprintf(
ioQQQ,
"DEBUG lgtimes but reset to %.2e\n" ,
dynamics.timestep );
1072 fprintf(
ioQQQ,
"DEBUG times exit dynamics.timestep %.2e elapsed_time %.2e scale %.2e ",
1075 rfield.time_continuum_scale);
1082 fprintf(
ioQQQ,
"\n" );
1091 dynamics.discretization_error = 0.;
1121 for(i=0;i<
nzone;++i)
1147 for( ion=0; ion<nelem+2; ++ion )
1169 for( nelem=ipISO; nelem<
LIMELM; ++nelem)
1171 if(
dense.lgElmtOn[nelem] )
1173 for(
long level=0; level <
iso_sp[ipISO][nelem].numLevels_local; ++level )
1185 dynamics.convergence_error +=
POW2(Oldi_iso/Oldi_density-
struc.StatesElem[nelem][nelem-ipISO][level][i]/
struc.hden[i]) ;
1220 fprintf(
ioQQQ,
"DYNAMICS DynaNewStep: Dyn_dr %.2e convergence_error %.2e discretization_error %.2e error_scale1 %.2e error_scale2 %.2e\n",
1245 for( i=0; i<
nzone; ++i )
1269 for( ion=0; ion<nelem+2; ++ion )
1276 for( nelem=ipISO; nelem<
LIMELM; ++nelem)
1278 if(
dense.lgElmtOn[nelem] )
1280 for(
long level=0; level <
iso_sp[ipISO][nelem].numLevels_max; ++level )
1282 Old_StatesElem[i][nelem][nelem-ipISO][level] =
struc.StatesElem[nelem][nelem-ipISO][level][i];
1313 flux *=
dense.xMassDensity0;
1348 dynamics.lgStatic_completed =
false;
1351 dynamics.lgTimeDependentStatic =
false;
1366 dynamics.convergence_tolerance = 0.1;
1391 dynamics.discretization_error = 0.;
1426 dynamics.Source[nelem] = ((
double*)
MALLOC( (
size_t)(nelem+2)*
sizeof(
double) ));
1428 for( ion=0; ion<nelem+2; ++ion )
1437 if(
dense.lgElmtOn[nelem] )
1440 for(
long ion=0; ion<nelem+1; ion++ )
1442 long ipISO = nelem-ion;
1460 if(
dense.lgElmtOn[nelem] )
1462 dynamics.StatesElem[nelem] = (
double**)
MALLOC(
sizeof(
double*)*(
unsigned)(nelem+1) );
1463 for(
long ion=0; ion<nelem+1; ion++ )
1465 long ipISO = nelem-ion;
1468 dynamics.StatesElem[nelem][nelem-ipISO] = (
double*)
MALLOC(
sizeof(
double)*(
unsigned)
iso_sp[ipISO][nelem].numLevels_max);
1473 dynamics.StatesElem[nelem][nelem-ipISO] = NULL;
1510 for( ns=0; ns <
struc.nzlim; ++ns )
1531 for( nelem=0; nelem<
LIMELM; ++nelem )
1537 for( nelem=0; nelem<
LIMELM; ++nelem )
1539 if(
dense.lgElmtOn[nelem] )
1543 for( ion=0; ion<nelem+1; ion++ )
1545 long ipISO = nelem-ion;
1561 for( i=0; i <
struc.nzlim; i++ )
1577 for( nelem=0; nelem<
LIMELM; ++nelem )
1579 for( ion=0; ion<
LIMELM+1; ++ion )
1586 for( nelem=ipISO; nelem<
LIMELM; ++nelem)
1588 if(
dense.lgElmtOn[nelem] )
1590 for(
long level=0; level <
iso_sp[ipISO][nelem].numLevels_max; ++level )
1626 pressure.lgPres_radiation_ON =
true;
1627 pressure.lgPres_magnetic_ON =
true;
1634 conv.EdenErrorAllowed = 1e-3;
1640 conv.HeatCoolRelErrorAllowed = 3e-4f;
1641 conv.PressureErrorAllowed = 1e-3f;
1645 conv.EdenErrorAllowed = 1e-5;
1646 conv.PressureErrorAllowed = 1e-5f;
1659 dynamics.lgTimeDependentStatic =
true;
1692 " Hit EOF while reading time-continuum list; use END to end list.\n" );
1700 while( p.
strcmp(
"END") != 0 )
1705 " Too many time points have been entered; the limit is %d. Increase variable NTIME in dynamics.c.\n",
1724 fprintf(
ioQQQ,
" Adding cycles with period = %e s.\n", period );
1730 fprintf(
ioQQQ,
" Hit EOF while reading line list; use END to end list.\n" );
1751 "scale factor to increase time",-1.);
1754 if( p.
nMatch(
"RECOMBIN") )
1772 fprintf(
ioQQQ,
" Hit EOF while reading line list; use END to end list.\n" );
1780 fprintf(
ioQQQ,
" At least two instances of time must be specified. There is an implicit instance at t=0.\n" \
1781 " The user must specify at least one additional time. Sorry.\n" );
1787 fprintf(
ioQQQ,
"DEBUG time dep %.2e %.2e %.2e %.2e\n",
1793 fprintf(
ioQQQ,
"\n" );
1801 bool lgModeSet=
false;
1828 wind.setBallistic();
1840 if ( 1 == iVelocity_Type && !lgModeSet)
1842 if (
wind.windv0 > 0.)
1844 fprintf(
ioQQQ,
"Warning, BALListic option needed to switch off pressure gradient terms\n");
1846 else if (
wind.windv0 == 0.)
1848 fprintf(
ioQQQ,
"Warning, STATic option needed for zero speed solutions\n");
1864 "centre of mass flux distribution");
1872 "power law index of mass flux distribution");
1876 if(iVelocity_Type == 1)
1904 else if(iVelocity_Type == 2)
1908 fprintf(
ioQQQ,
"Can't specify gradient when flux is constant!\n");
1926 wind.windv0 = -0.01f;
1935 if (
wind.windv0 < 0.)
1939 else if (
wind.windv0 > 0.)
1941 wind.setBallistic();
1960 fprintf(
ioQQQ,
"Scale %g (*%c) Index %g Center %g\n",
1970 strcpy(
dense.chDenseLaw,
"DYNA" );
1976 if(
wind.windv0 <= 1.e6 )
1979 fprintf(
ioQQQ,
" >>>>Initial wind velocity should be greater than speed of sound; calculation only valid above sonic point.\n" );
1980 wind.lgWindOK =
false;
1988 wind.lgDisk =
false;
1992 strcpy(
dense.chDenseLaw,
"WIND" );
2018 fprintf(
ioQQQ,
" DYNAMICS Advection: Uad %.2f Uwd%.2e FRCcool: %4.2f Heat %4.2f\n",
2019 timesc.sound_speed_adiabatic/1e5 ,
2027 fprintf(
ioQQQ,
" DYNAMICS Eexcit:%.4e Eion:%.4e Ebin:%.4e Ekin:%.4e ET+pdv %.4e EnthalpyDensity/rho%.4e AdvSpWork%.4e\n",
2044 if( strcmp( chJob ,
"END" ) == 0 )
2126 "%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\n" ,
2129 rfield.time_continuum_scale ,
2155 fprintf( ipPnunit ,
"%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
2213 rfield.time_continuum_scale = 1.;
2227 fprintf(
ioQQQ,
" PROBLEM - DynaIterStart - I need the continuum at time %.2e but the table ends at %.2e.\n" ,
2245 fprintf(
ioQQQ,
"DEBUG time dep reset continuum iter %ld dynamics.timestep %.2e elapsed time %.2e scale %.2e",
2249 rfield.time_continuum_scale);
2252 fprintf(
ioQQQ,
" recom");
2254 fprintf(
ioQQQ,
"\n");
2257 long int nTimeVary = 0;
2258 for(
long int i=0; i <
rfield.nShape; i++ )
2262 if(
rfield.lgTimeVary[i] )
2266 if(
hextra.lgTurbHeatVaryTime )
2270 fprintf(
ioQQQ,
"DEBUG TurbHeat vary new heat %.2e\n",
2273 else if( !nTimeVary )
2275 fprintf(
ioQQQ,
" DISASTER - there were no variable continua "
2276 "or heat sources - put TIME option on at least one "
2277 "luminosity or hextra command.\n");
bool fp_bound(sys_float lo, sys_float x, sys_float hi, int n=3)
NORETURN void TotalInsanity(void)
sys_float SDIV(sys_float x)
#define DEBUG_ENTRY(funcname)
int cdIonFrac(const char *chLabel, long int IonStage, double *fracin, const char *chWeight, bool lgDensity)
int cdTemp(const char *chLabel, long int IonStage, double *TeMean, const char *chWeight)
bool nMatch(const char *chKey) const
double getNumberCheckAlwaysLog(const char *chDesc)
int strcmp(const char *s2)
double getNumberPlain(const char *chDesc)
double getNumberCheck(const char *chDesc)
double getNumberCheckAlwaysLogLim(const char *chDesc, double flim)
double getNumberDefault(const char *chDesc, double fdef)
double getNumberDefaultAlwaysLog(const char *chDesc, double fdef)
realnum GetHubbleFactor(realnum z)
UNUSED const realnum BIGFLOAT
realnum scalingDensity(void)
realnum scalingZoneDensity(long i)
static double ** UpstreamIon
static realnum * Old_ednstr
void ParseDynaWind(Parser &p)
STATIC void DynaSaveLast(void)
static double AdvecSpecificEnthalpy
void ParseDynaTime(Parser &p)
static realnum * Old_histr
static double * time_flux_ratio
static realnum * Old_hiistr
static double * time_dt_scale_factor
static realnum **** Old_StatesElem
static realnum ** Old_gas_phase
static long int nOld_zone
static double *** UpstreamStatesElem
STATIC void advection_set_default(bool lgWind)
void DynaSave(FILE *ipPnunit, char chJob)
static realnum ** Old_molecules
static realnum * Old_xLyman_depth
static double * Upstream_molecules
STATIC double timestep_next(void)
static realnum * Old_pressure
static realnum * Old_depth
static realnum * Old_density
static realnum *** Old_xIonDense
void DynaPunchTimeDep(FILE *ipPnunit, const char *chJob)
STATIC void DynaNewStep(void)
static double * UpstreamElem
realnum DynaFlux(double depth)
static realnum * Old_DenMass
static double * time_elapsed_time
static long int nTime_flux
void DynaCreateArrays(void)
static realnum * EnthalpyDensity
static realnum * Old_EnthalpyDensity
t_iso_sp iso_sp[NISO][LIMELM]
t_mole_global mole_global
molezone * findspecieslocal(const char buf[])
double linint(const double x[], const double y[], long n, double xval)