Changeset 2449 for branches/newmole/source

Timestamp:

11/09/08 12:46:54 (2 months ago)

Author:

rjrw

Message:

Merged from trunk r2439:2448

Location:

branches/newmole/source

Files:

• cddefines.h (modified) (2 diffs)
• dynamics.cpp (modified) (2 diffs)
• hydro_vs_rates.cpp (modified) (1 diff)
• TestFpCheck.cpp (copied) (copied from trunk/source/TestFpCheck.cpp)

branches/newmole/source/cddefines.h

@@ -788 +788 @@
 inline bool fp_equal( sys_float x, sys_float y, int n=3 )
 {
-#       ifdef _MSC_VER
-        /* disable warning that conditional expression is constant, true or false in if */
-#               pragma warning( disable : 4127 )
-#       endif
         ASSERT( n >= 1 );
         // mimic IEEE behavior
@@ -849 +845 @@
         return ( abs( x-y ) <= tol );
 }
+
+/** checks whether a number is within bounds */
+inline bool fp_bound( sys_float lo, sys_float x, sys_float hi, int n=3 )
+{
+        ASSERT( n >= 1 );
+        // mimic IEEE behavior
+        if( isnan(x) || isnan(lo) || isnan(hi) )
+                return false;
+        if( fp_equal(lo,hi,n) )
+                return fp_equal(0.5f*(lo+hi),x,n);
+        if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -((sys_float)n+0.1f)*FLT_EPSILON )
+                return false;
+        return true;
+}
+inline bool fp_bound( double lo, double x, double hi, int n=3 )
+{
+        ASSERT( n >= 1 );
+        // mimic IEEE behavior
+        if( isnan(x) || isnan(lo) || isnan(hi) )
+                return false;
+        if( fp_equal(lo,hi,n) )
+                return fp_equal(0.5*(lo+hi),x,n);
+        if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -((double)n+0.1)*DBL_EPSILON )
+                return false;
+        return true;
+}
+inline bool fp_bound_tol( sys_float lo, sys_float x, sys_float hi, sys_float tol )
+{
+        ASSERT( tol > 0.f );
+        // mimic IEEE behavior
+        if( isnan(x) || isnan(lo) || isnan(hi) )
+                return false;
+        if( fp_equal_tol(lo,hi,tol) )
+                return fp_equal_tol(0.5f*(lo+hi),x,tol);
+        if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -tol )
+                return false;
+        return true;
+}
+inline bool fp_bound_tol( double lo, double x, double hi, double tol )
+{
+        ASSERT( tol > 0. );
+        // mimic IEEE behavior
+        if( isnan(x) || isnan(lo) || isnan(hi) )
+                return false;
+        if( fp_equal_tol(lo,hi,tol) )
+                return fp_equal_tol(0.5*(lo+hi),x,tol);
+        if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -tol )
+                return false;
+        return true;
+}
+
 
 #undef POW2

branches/newmole/source/dynamics.cpp

@@ -827 +827 @@
 
         /* this is index of upstream cell in structure stored from previous iteration */
-        double upstream, dilution, dilutionleft, dilutionright, frac_next, testval;
+        double upstream, dilution, dilutionleft, dilutionright, frac_next;
 
         /* Properties for cell half as far upstream, used to converge timestep */
@@ -907 +907 @@
                         frac_next);
                 
-                testval = (Old_EnthalpyDensity[ipUpstream]*dilutionleft - AdvecSpecificEnthalpy) *
-                        (AdvecSpecificEnthalpy - Old_EnthalpyDensity[ipUpstream+1]*dilutionright);
-
                 ASSERT( Old_depth[ipUpstream] <= upstream && upstream <= Old_depth[ipUpstream+1] );
-                if( testval <= -SMALLFLOAT )
-                {
-                        fprintf(ioQQQ,"PROBLEM advected enthalpy density is not conserved %e\t%e\t%e\t%e\n",
-                                testval ,
-                                Old_EnthalpyDensity[ipUpstream]*dilutionleft,
-                                AdvecSpecificEnthalpy,
-                                Old_EnthalpyDensity[ipUpstream+1]*dilutionright);
-                }
-
-                ASSERT( testval > -SMALLFLOAT );
+
+                // we have a mix of realnum and double here, so make sure realnum is used for the test...
+                realnum lo = (realnum)(Old_EnthalpyDensity[ipUpstream]*dilutionleft);
+                realnum x = (realnum)AdvecSpecificEnthalpy;
+                realnum hi = (realnum)(Old_EnthalpyDensity[ipUpstream+1]*dilutionright);
+                if( ! fp_bound(lo,x,hi) )
+                {
+                        fprintf(ioQQQ,"PROBLEM advected enthalpy density is not conserved %.16e\t%.16e\t%.16e\t%e\t%e\n",
+                                lo, x, hi, (hi-x)/(hi-lo), (x-lo)/(hi-lo));
+                }
+
+                ASSERT( fp_bound(lo,x,hi) );
+
                 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
                 {

branches/newmole/source/hydro_vs_rates.cpp

@@ -36 +36 @@
         double coll_str;
 
-        global_ipISO = ipISO;
-        global_nelem = nelem;
-        global_ipHi = ipHi;
-        global_ipLo = ipLo;
-        global_temp = temp;
-        global_Collider = Collider;
-        global_Aul = Aul;
-
-        /* evaluate collision strength, two options,
-         * do thermal average (very slow) if set with
-         * SET COLLISION STRENGTH AVERAGE command,
-         * default just do point at kT
-         * tests show no impact on test suite, much faster */
-        if( iso.lgCollStrenThermAver )
-        {
-                /* do average over Maxwellian */
-                coll_str =  qg32( 0.0, 1.0, Therm_ave_coll_str_int_VS80);
-                coll_str += qg32( 1.0, 10.0, Therm_ave_coll_str_int_VS80);
-        }
-        else
-        {
-                /* the argument to the function is equivalent to evaluating the function at
-                 * hnu = kT */
-                coll_str = hydro_vs_coll_str( EVRYD*global_temp/TE1RYD );
+        if( Collider == ipELECTRON )
+        {
+                coll_str = hydro_vs_deexcit( ipISO, nelem, ipHi, ipLo, Aul);
+        }
+        else
+        {
+                global_ipISO = ipISO;
+                global_nelem = nelem;
+                global_ipHi = ipHi;
+                global_ipLo = ipLo;
+                global_temp = temp;
+                global_Collider = Collider;
+                global_Aul = Aul;
+
+                /* evaluate collision strength, two options,
+                 * do thermal average (very slow) if set with
+                 * SET COLLISION STRENGTH AVERAGE command,
+                 * default just do point at kT
+                 * tests show no impact on test suite, much faster */
+                if( iso.lgCollStrenThermAver )
+                {
+                        /* do average over Maxwellian */
+                        coll_str =  qg32( 0.0, 1.0, Therm_ave_coll_str_int_VS80);
+                        coll_str += qg32( 1.0, 10.0, Therm_ave_coll_str_int_VS80);
+                }
+                else
+                {
+                        /* the argument to the function is equivalent to evaluating the function at
+                         * hnu = kT */
+                        coll_str = hydro_vs_coll_str( EVRYD*global_temp/TE1RYD );
+                }
         }