Changeset 2110 for trunk/source/grains_qheat.cpp

Timestamp:

06/16/08 21:10:04 (7 months ago)

Author:

peter

Message:

source/grains_qheat.cpp:

Bug-fix: fix PR69 - failed assert in quantum heating code.

Files:

• trunk/source/grains_qheat.cpp (modified) (10 diffs)

trunk/source/grains_qheat.cpp

@@ -622 +622 @@
         if( trace.lgTrace && trace.lgDustBug )
         {
-                fprintf( ioQQQ, "   grain heating: %.4e, integral %.4e, total rate %.4e\n",
-                         gv.bin[nd]->GrainHeat,integral,Phi[0]);
+                double Rate2 = 0.;
+                for( int nz=0; nz < gv.bin[nd]->nChrg; nz++ )
+                        Rate2 += gv.bin[nd]->chrg[nz]->FracPop*gv.bin[nd]->chrg[nz]->HeatingRate2;
+
+                fprintf( ioQQQ, "   grain heating: %.4e, integral %.4e, total rate %.4e lgNegRate %c\n",
+                         gv.bin[nd]->GrainHeat,integral,Phi[0],TorF(lgNegRate));
                 fprintf( ioQQQ, "   av grain temp %.4e av grain enthalpy (Ryd) %.4e\n",
                          gv.bin[nd]->tedust,Umax);
                 fprintf( ioQQQ, "   fwhm^2/(4ln2*c2/c0): %.4e fwhm (Ryd) %.4e fwhm/Umax %.4e\n",
                          NumEvents,fwhm,FwhmRatio );
-                fprintf( ioQQQ, "   lgQHTooWide %d\n", gv.bin[nd]->lgQHTooWide );
+                fprintf( ioQQQ, "   HeatingRate1 %.4e HeatingRate2 %.4e lgQHTooWide %c\n",
+                         gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3, Rate2*gv.bin[nd]->cnv_H_pCM3,
+                         TorF(gv.bin[nd]->lgQHTooWide) );
         }
 
@@ -911 +917 @@
          * energies, the reshuffling for higher energies is done in the next loop
          * phiTilde has units events/H/s/cell at default depletion */
+
+        double NegHeatingRate = 0.;
 
         for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
@@ -983 +991 @@
                 *check += gptr->FracPop*check1*EN1RYD*gv.bin[nd]->cnv_H_pCM3;
 
+                sum += gptr->FracPop*check1*EN1RYD*gv.bin[nd]->cnv_H_pCM3;
+
                 if( DEBUG_LOC )
                 {
@@ -990 +1000 @@
                                 integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
                         }
-                        sum += gptr->FracPop*check1*EN1RYD*gv.bin[nd]->cnv_H_pCM3;
                         dprintf( ioQQQ, " integral test 1: integral %.6e %.6e\n", integral, sum );
                 }
@@ -1007 +1016 @@
 
                 /* >>chng 03 nov 06, check for extremely low HeatingRate and save CPU time, pah_crash.in, PvH */
-                if( fabs(gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3) > 0.05*CONSERV_TOL*gv.bin[nd]->GrainHeat ) 
+                if( gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3 > 0.05*CONSERV_TOL*gv.bin[nd]->GrainHeat ) 
                 {
                         double *RateArr,Sum,ESum,DSum,Delta,E_av2,Corr;
@@ -1076 +1085 @@
                         }
 
+                        sum += gptr->FracPop*gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3;
+
                         if( DEBUG_LOC )
                         {
@@ -1083 +1094 @@
                                         integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
                                 }
-                                sum += gptr->FracPop*gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3;
                                 dprintf( ioQQQ, " integral test 2: integral %.6e %.6e\n", integral, sum );
                         }
 
                         free( RateArr );
+                }
+                else
+                {
+                        NegHeatingRate += gptr->FracPop*gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3;
                 }
         }
@@ -1109 +1123 @@
 
         /* >>chng 03 nov 06, check for extremely low HeatingRate and save CPU time, PvH */
-        if( fabs(gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3) > 0.05*CONSERV_TOL*gv.bin[nd]->GrainHeat )
+        if( gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3 > 0.05*CONSERV_TOL*gv.bin[nd]->GrainHeat )
         {
                 /* limits for Taylor expansion of (1+x)*exp(-x) */
@@ -1153 +1167 @@
                         ylo = yhi;
                 }
-        }
-
-        if( DEBUG_LOC )
-        {
-                double integral = 0.;
+
+                sum += gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3;
+
+                if( DEBUG_LOC )
+                {
+                        double integral = 0.;
+                        for( i=0; i < gv.bin[nd]->qnflux; i++ )
+                        {
+                                integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
+                        }
+                        dprintf( ioQQQ, " integral test 3: integral %.6e %.6e\n", integral, sum );
+                }
+        }
+        else
+        {
+                NegHeatingRate += gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3;
+        }
+
+        /* here we account for the negative heating rates, we simply do that by scaling the entire
+         * phiTilde array down by a constant factor such that the total amount of energy is conserved
+         * This treatment assures that phiTilde never goes negative, which avoids problems further on */
+        if( NegHeatingRate < 0. )
+        {
+                double scale_fac = (sum+NegHeatingRate)/sum;
                 for( i=0; i < gv.bin[nd]->qnflux; i++ )
-                {
-                        integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
-                }
-                sum += gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3;
-                dprintf( ioQQQ, " integral test 3: integral %.6e %.6e\n", integral, sum );
-        }
+                        phiTilde[i] *= scale_fac;
+
+                sum += NegHeatingRate;
+
+                if( DEBUG_LOC )
+                {
+                        double integral = 0.;
+                        for( i=0; i < gv.bin[nd]->qnflux; i++ )
+                        {
+                                integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
+                        }
+                        dprintf( ioQQQ, " integral test 4: integral %.6e %.6e\n", integral, sum );
+                }
+        }
+
         return;
 }
@@ -1354 +1396 @@
                 RadCooling += dCool;
 
-/*              if( nzone >= 870 ) { */
-/*              printf(" k %ld T[k] %.6e U[k] %.6e p[k] %.6e dPdlnT[k] %.6e\n",k-1,qtemp[k-1],u1[k-1],p[k-1],dPdlnT[k-1]); */
-/*              printf(" k %ld T[k] %.6e U[k] %.6e p[k] %.6e dPdlnT[k] %.6e\n",k,qtemp[k],u1[k],p[k],dPdlnT[k]); */
-/*              } */
+//              if( nzone >= 24 && nd == 0 ) {
+//              printf(" k %ld T[k] %.6e U[k] %.6e p[k] %.6e dPdlnT[k] %.6e\n",k-1,qtemp[k-1],u1[k-1],p[k-1],dPdlnT[k-1]);
+//              printf(" k %ld T[k] %.6e U[k] %.6e p[k] %.6e dPdlnT[k] %.6e\n",k,qtemp[k],u1[k],p[k],dPdlnT[k]);
+//              }
 
                 /* if qtmin is far too low, p[k] can become extremely large, exceeding