Cloudy & Associates

Commit 7a61cef8 authored by Ferland, Gary's avatar Ferland, Gary
Browse files

Protection against division by zero cooling.

These protect against division by zero.  Problems were exposed
by the lab sim posted in gtests / bugs / Mancini.  The
sim does get now much farther along but aborts when the line
stack becomes misaligned.

The sim is a lab plasma, so it is far from our home base, but we
should be able to compute it.  It is non-equilibrium.  It is
unclear why the cooling is zero - perhaps the total is net heating.

The changes seem safe since they are have never been needed and
will have no effect in most cases.
parent c00522d3
......@@ -39,14 +39,17 @@ void AgeCheck(void)
for( i=0; i < limit; i++ )
{
timesc.time_therm_long =
MAX2( timesc.time_therm_long ,
struc.DenParticles[i]*BOLTZMANN*1.5*struc.testr[i]/struc.coolstr[i]);
timesc.time_therm_short =
MIN2( timesc.time_therm_short ,
struc.DenParticles[i]*BOLTZMANN*1.5*struc.testr[i]/struc.coolstr[i]);
/*>>chng 99 feb 01, had div by heating, changed to cooling so constant
* temperature models are more realistic */
if( struc.coolstr[i]>SMALLFLOAT )
{
timesc.time_therm_long =
MAX2( timesc.time_therm_long ,
struc.DenParticles[i]*BOLTZMANN*1.5*struc.testr[i]/struc.coolstr[i]);
timesc.time_therm_short =
MIN2( timesc.time_therm_short ,
struc.DenParticles[i]*BOLTZMANN*1.5*struc.testr[i]/struc.coolstr[i]);
/*>>chng 99 feb 01, had div by heating, changed to cooling so constant
* temperature models are more realistic */
}
}
tlong = MAX2(tlong,timesc.time_therm_long);
......
......@@ -67,7 +67,10 @@ void CoolSave(FILE * io, const char chJob[])
if( strncmp( thermal.chClntLab[i], "adve", 4 ) == 0 )
continue;
csav[ip] = (realnum)( safe_div( MAX2(thermal.cooling[i],thermal.heatnt[i]), cool_total, 0. ));
if( cool_total>SMALLFLOAT )
csav[ip] = (realnum)( safe_div( MAX2(thermal.cooling[i],thermal.heatnt[i]), cool_total, 0. ));
else
csav[ip] = 0.;
/* save sign to remember if heating or cooling line */
if( thermal.heatnt[i] == 0. )
{
......
......@@ -41,12 +41,16 @@ void t_timesc::calc_therm_timesc( long int izone )
{
/* >>chng 99 feb 01, had div by heating, changed to cooling so constant
* temperature models are more realistic */
double dt = 1.5 * BOLTZMANN * struc.DenParticles[i] * struc.testr[i] / struc.coolstr[i];
time_therm_long = MAX2( time_therm_long , dt );
time_therm_short= MIN2( time_therm_short, dt );
// printf("dt = %g\t long = %g\t short = %g\n", dt, time_therm_long, time_therm_short);
}
if( struc.coolstr[i] > SMALLFLOAT )
{
double dt = 1.5 * BOLTZMANN * struc.DenParticles[i] * struc.testr[i] / struc.coolstr[i];
//double dt = 1.5 * BOLTZMANN * struc.DenParticles[i] * safe_div( (double)struc.testr[i] , struc.coolstr[i] );
time_therm_long = MAX2( time_therm_long , dt );
time_therm_short= MIN2( time_therm_short, dt );
// printf("dt = %g\t long = %g\t short = %g\n", dt, time_therm_long, time_therm_short);}
}
// printf( "*** long = %g\t short = %g\n", time_therm_long, time_therm_short );
}
return;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment