Changeset 1714 for branches/newmole/source

Timestamp:

12/22/07 09:53:10 (11 months ago)

Author:

rjrw

Message:

Revise coding of grain photoevaporation rates.

Switch back error scaling in mole_newton_step to get more decks working.

Location:

branches/newmole/source

Files:

• ion_solver.cpp (modified) (1 diff)
• mole.h (modified) (1 diff)
• mole_newton_step.cpp (modified) (4 diffs)
• mole_reactions.cpp (modified) (1 diff)
• mole_species.cpp (modified) (3 diffs)
• punch_do.cpp (modified) (1 diff)

branches/newmole/source/ion_solver.cpp

@@ -68 +68 @@
                 if (abund_total < -10.0*FLT_EPSILON*dense.gas_phase[nelem])
                 {
+                        double sum;
                         fprintf(stdout," DISASTER ion_solver: Bad density %li %g %g %g\n",nelem,dense.gas_phase[nelem],dense.xMolecules[nelem],
                                                         dense.gas_phase[nelem]-dense.xMolecules[nelem]);
+                        sum = 0;
+                        fprintf(stdout," Levels:");
+                        for (ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
+                        {
+                                fprintf(stdout," %ld %g;",ion,dense.xIonDense[nelem][ion]);
+                                sum += dense.xIonDense[nelem][ion];
+                        }
+                        fprintf(stdout," tot=%g\n",sum);
                         puts( "[Stop in ion_solver]" );
                         cdEXIT(EXIT_FAILURE);

branches/newmole/source/mole.h

@@ -107 +107 @@
         int num_total, num_calc, num_compacted;
         
-        realnum eden_f, grain_area;
+        realnum eden_f, grain_area, grain_avg_area;
 
         /** total charge in molecules */

branches/newmole/source/mole_newton_step.cpp

@@ -26 +26 @@
 
 #define ABSLIM  1e-12
+#define ERRLIM  1e-12
 #define SMALLABUND 1e-24
 #       ifdef MAT
@@ -291 +292 @@
                 {
                         /* Smooth the error mode tailoff */
-                        etmp = b[groupspecies[i]->index];
-                        //((ABSLIM+b0vec[i])*(-c[groupspecies[i]->index][groupspecies[i]->index]));
+                        etmp = b[groupspecies[i]->index]/
+                                ((ERRLIM+fabs(b0vec[i]))*(ERRLIM+fabs(c[groupspecies[i]->index][groupspecies[i]->index])));
                         etmp *= etmp;
                         if (etmp > emax)
@@ -477 +478 @@
                                         iworst = i;
                                 }
-                                if (b0vec[i] > SMALLABUND) // || b2vec[i] < -SMALLABUND)
+                                if (b0vec[i] > SMALLABUND || b2vec[i] < -SMALLABUND)
                                         ++*nBad;
                         }
@@ -488 +489 @@
                 updateMolecules( b2vec, fion );         
         }
+        //fprintf(ioQQQ,"Final fstep %g\n",fstep);
 
         *error = (realnum) MIN2(error0,1e30);

branches/newmole/source/mole_reactions.cpp

@@ -990 +990 @@
          * The following treatment uses their equation 3 */
 
-        double grain_area;
-        int nd;
-
         DEBUG_ENTRY("grn_photo()");
 
-        grain_area = 0.;
-        for( nd=0; nd < gv.nBin; nd++ )
-        {
-                grain_area += gv.bin[nd]->AvArea;
-        }
-
-        /* This is the number of molecules removed from the grain per incident photon
-         * use same value of 1e-4 as used in Bergin */
-
-        return grain_area* hmi.UV_Cont_rel2_Draine_DB96_depth;
+        /* This is the number of molecules removed from the grain per
+         * incident photon use same value of 1e-4 as used in Bergin, include
+         * conversion to radiation field units to get the dimensions right
+         * (cf the d'Hendercourt reference in Bergin et al) */
+
+        fixit(); // Should this be 4* (i.e. the total not the cs area)?  Depends on definition of radiation field.
+        return mole.grain_avg_area * hmi.UV_Cont_rel2_Draine_DB96_depth *(1.232e7f * 1.71f);
 }
 

branches/newmole/source/mole_species.cpp

@@ -212 +212 @@
                 sp = newspecies("H2Ogrn",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, -238.9f);
                 sp = newspecies("OHgrn ",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, 38.4f);
-                if (0)
+                if (1)
                 {
                         sp = newspecies("Cgrn  ",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, 711.2f);
@@ -596 +596 @@
         int i;
         int nd;
-        double den_times_area;
+        double den_times_area, dengrains;
 
         DEBUG_ENTRY("mole_update_species_cache()");
@@ -611 +611 @@
         /* >>chng 06 feb 28, turn off this rate when no grain molecules */
         /* >>chng 06 dec 05 rjrw: do this in newreact rather than rate */
-        den_times_area = 0.0;
-        for( nd=0; nd < gv.nBin; nd++ )
-        {
-                /* >>chng 06 mar 04, update expression for projected grain surface area, PvH */
-                den_times_area += gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
-        }
-
-        mole.grain_area = den_times_area;
+        if(gv.lgDustOn)
+        {
+                den_times_area = dengrains = 0.0;
+                for( nd=0; nd < gv.nBin; nd++ )
+                {
+                        /* >>chng 06 mar 04, update expression for projected grain surface area, PvH */
+                        den_times_area += gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
+                        dengrains += gv.bin[nd]->cnv_GR_pCM3;
+                }
+                
+                mole.grain_area = den_times_area;
+                mole.grain_avg_area = den_times_area/dengrains;
+        }
+        else
+        {
+                mole.grain_area = mole.grain_avg_area = 0.0;
+        }
 
         for (i=0;i<mole.num_total;i++) 

branches/newmole/source/punch_do.cpp

@@ -2530 +2530 @@
                                         for(i=0; i<mole.num_calc; ++i )
                                         {
-                                                fprintf( punch.ipPnunit[ipPun], "\t%-4.4s", COmole[i]->label );
+                                                fprintf( punch.ipPnunit[ipPun], "\t%-7.7s", COmole[i]->label );
                                         }
                                         fprintf ( punch.ipPnunit[ipPun], "\n");