Changeset 2107

Timestamp:

06/11/08 02:41:14 (3 months ago)

Author:

gary

Message:

iso_level.cpp - use difference in ionization / recombination rates in source/sink for exchange with next resolved atoms in next lower ionization stage. At high densities excited states have very large ion/rec rates which cancel one another when level is in LTE. Using the differences is more stable fur such densities.

trivial changes in comments
ion_solver.cpp
iso.h
iso_ionize_recombine.cpp

with this change in iso_level.cpp test suite has following issues:
blr_hizqso.out: crashed
blr_n09_p20_Z20.out: crashed
limit_conserve.out: crashed
all crashes are due to radius_increment finds inconsistent populations
radius_increment.cpp at line number 152

following have asserts thrown;
may be serious, should look at:
orion_hii_pdr - molecular ions changes, esp HeH+

trivial changes;
pdr_coolbb
blr_n13_p18_Z20 - took more iterations per zone
blr_n12_p19 slight change in 1640

known problem unrelated to this:
grains_conserve

Location:

trunk

Files:

• source/ion_solver.cpp (modified) (2 diffs)
• source/iso.h (modified) (1 diff)
• source/iso_ionize_recombine.cpp (modified) (1 diff)
• source/iso_level.cpp (modified) (5 diffs)
• tsuite/auto/grid_line_ratios.in (modified) (1 diff)

trunk/source/ion_solver.cpp

@@ -600 +600 @@
         {
                 /* option to print matrix */
-                /*@-redef@*/
                 enum {DEBUG_LOC=false};
-                /*@+redef@*/
                 if( DEBUG_LOC && nzone==1 && lgPrintIt )
                 {
@@ -740 +738 @@
 
         {
-                /*@-redef@*/
                 /* this is to debug following failed assert */
                 enum {DEBUG_LOC=false};
-                /*@+redef@*/
                 if( DEBUG_LOC && nelem == ipHYDROGEN )
                 {

trunk/source/iso.h

@@ -254 +254 @@
         multi_arr<double,3> xIsoLevNIonRyd;
 
-        /** the ratio fully stripped ion to atom for all hydrogenic species
-         * HIonSimple is simple estimate, should agree at low density */
+        /** the ratio of ion to atom for all iso species
+         * xIonSimple is simple estimate, should agree at low density */
         double xIonSimple[NISO][LIMELM];
 

trunk/source/iso_ionize_recombine.cpp

@@ -241 +241 @@
         }
 
-        /* fraction of total recombs due to three body - when most are due to this
-         * small changes in temperature can produce large changes in rec coef,
+        /* fraction of total recombinations due to three body - when most are due to this
+         * small changes in temperature can produce large changes in recombination coefficient,
          * and so in ionization */
         iso.RecomCollisFrac[ipISO][nelem] = Recom3Body* dense.eden / ionbal.RateRecomTot[nelem][nelem-ipISO];
 
-        /* get simple estimate of level of ionization */
+        /* get simple estimate of atom to ion ratio */
         if( ionbal.RateRecomTot[nelem][nelem-ipISO] > 0. )
         {

trunk/source/iso_level.cpp

@@ -11 +11 @@
 #include "grainvar.h"
 #include "he.h"
-#include "heavy.h"
 #include "helike.h"
 #include "hydrogenic.h"
@@ -24 +23 @@
 #include "thirdparty.h"
 #include "trace.h"
-#include "yield.h"
 
 /* solve for level populations  */
@@ -90 +88 @@
 
         /* these two collision rates must be the same or we are in big trouble,
-         * since used interchangeably */
+         * since used interchangably */
         ASSERT( ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0]< SMALLFLOAT ||
                 fabs( iso.ColIoniz[ipISO][nelem][0]* collider /
@@ -375 +373 @@
 #endif
 
+                /* derive source/sink terms for ion/rec terms coupling resolved atoms 
+                 * to lower ionization stages */
                 if( nelem-ipISO >= 1 && ionbal.RateIonizTot[nelem][nelem-ipISO-1] > 0.)
                 {
+                        double sum_recom =0. , sum_ioniz = 0.;
+                        sum_popn_ov_ion = 0.;
+                        /* if lower stage of ionization is itself resolved then we must
+                         * sum over rates between levels and its continuum, the ion we
+                         * are now solving */
+                        if( ipISO+1<NISO )
+                        {
+                                for( level=0; level < iso.numLevels_local[ipISO+1][nelem]; level++ )
+                                {
+                                        sum_popn_ov_ion += StatesElem[ipISO+1][nelem][level].Pop;
+
+                                        /* sum of all ionization processes from levels to ion we
+                                         * now consider - use difference in ionization / recombination 
+                                         * since at high densities each of these terms are very
+                                         * large but nearly cancel */
+                                        double one = iso.RateCont2Level[ipISO+1][nelem][level] -
+                                                StatesElem[ipISO+1][nelem][level].Pop * 
+                                                iso.RateLevel2Cont[ipISO+1][nelem][level];
+
+                                        /* sign determines whether this is net ionization or 
+                                         * recombination process.  Ground & metastable states have 
+                                         * large departure coefficients so are net ionization,
+                                         * excited states have departure coefficients that
+                                         * asymptotically approach unity from below, so are net
+                                         * recombination */
+                                        if( one > 0 )
+                                                sum_recom += one;
+                                        else
+                                                sum_ioniz -= one;
+                                }
+                                ASSERT( sum_recom>=0. && sum_ioniz>=0. );
+                                sum_ioniz  /= MAX2(SMALLDOUBLE , sum_popn_ov_ion);
+                        }
+                        else
+                        {
+                                /* lower stage is not resolved */
+                                sum_ioniz = ionbal.RateIonizTot[nelem][nelem-ipISO-1];
+                                sum_recom = ionbal.RateRecomTot[nelem][nelem-ipISO-1];
+                        }
                         /* now add ionization from next lower stage */
-                        source += ionbal.RateIonizTot[nelem][nelem-ipISO-1]*
+                        source += sum_ioniz*
                                 dense.xIonDense[nelem][nelem-ipISO-1];  
 
-                        /* now add recombination to next lower stage ground */
-                        sink += ionbal.RateRecomTot[nelem][nelem-ipISO-1];
-                }
-                
+                        /* now add recombination to next lower stage, units are s-1 */
+                        sink += sum_recom;
+                }
+
                 creation[0] += source/SDIV(dense.xIonDense[nelem][nelem+1-ipISO]);
-
                 for( level=0; level < numlevels_local; level++ )
                 {
@@ -577 +615 @@
 
         /* all solvers end up here */
-        /* sum_popn_ov_ion will become the ratio of iso to parent
+
+        /* sum_popn_ov_ion will become the ratio of this atom to parent
          * species, create sum of level pops per ion first */
         sum_popn_ov_ion = 0.;
 
-        /* this is total ionization of this species referenced to the ground state */
+        /* this is total ionization rate, s-1, of this species referenced to 
+         * the total abundance */
         ionbal.RateIonizTot[nelem][nelem-ipISO] = 0.;
-
         for( level=0; level < numlevels_local; level++ )
         {

trunk/tsuite/auto/grid_line_ratios.in

@@ -2 +2 @@
 c 
 c commands controlling continuum =====
-blackbody 4e4K
+blackbody 40,000K
 ionization parameter -2
 c