Changeset 1838

Timestamp:

03/09/08 16:59:05 (10 months ago)

Author:

rjrw

Message:

Minor changes to resulting from trying to chase down differences on
pdf_leiden_f1 --

ion_solver: allow molecule effects in first iteration,
iso_level: defensive programming for z[][]
mole_eval_balance: move some definitions for sp
mole_co_atom: add SDIV (aka fixit...) to prevent FPE in some circumstances

Location:

branches/newmole/source

Files:

• ion_solver.cpp (modified) (3 diffs)
• iso_level.cpp (modified) (4 diffs)
• mole_co_atom.cpp (modified) (1 diff)
• mole_eval_balance.cpp (modified) (3 diffs)

branches/newmole/source/ion_solver.cpp

@@ -378 +378 @@
          * matrix is singular */
         totsrc = 0.;
-        /*>>chng 06 nov 28 only include source from molecules if we have an estimated first 
-         * solution - first test is that we have called mole at least twice,
-         * second test is that we are on a later iteration */
-        if( conv.nTotalIoniz > 1 || iteration > 1 )
-        {
-
-                for( i=0; i<ion_range;i++ )
-                {
-                        ion = i+ion_low;
-
-                        /* these are the external source and sink terms */
-                        /* need negative sign to get positive pops */
-                        source[i] -= mole.source[nelem][ion];
-                        totsrc += mole.source[nelem][ion];
-                        MAT( xmat, i, i ) -= mole.sink[nelem][ion];
-                }
+
+        for( i=0; i<ion_range;i++ )
+        {
+                ion = i+ion_low;
+                
+                /* these are the external source and sink terms */
+                /* need negative sign to get positive pops */
+                source[i] -= mole.source[nelem][ion];
+                totsrc += mole.source[nelem][ion];
+                MAT( xmat, i, i ) -= mole.sink[nelem][ion];
         }
 
@@ -1030 +1024 @@
                 fprintf( ioQQQ, "\n" );
 
+                /* grain ionization - not total rate but should be dominant process */
+                fprintf( ioQQQ, " %s ion xfr mol "  ,elementnames.chElementSym[nelem]);
+                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
+                {
+                        fprintf( ioQQQ, " %9.2e", mole.xMoleChTrRate[nelem][ion][ion+1] );
+                }
+                fprintf( ioQQQ, "\n" );
+
                 /* charge exchange ionization */
                 fprintf( ioQQQ, " %s chr trn ion " ,elementnames.chElementSym[nelem] );
@@ -1080 +1082 @@
                 {
                         fprintf( ioQQQ, " %9.2e", gv.GrainChTrRate[nelem][ion+1][ion] );
+                }
+                fprintf( ioQQQ, "\n" );
+
+                fprintf( ioQQQ, " %s rec xfr mol "  ,elementnames.chElementSym[nelem]);
+                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
+                {
+                        fprintf( ioQQQ, " %9.2e", mole.xMoleChTrRate[nelem][ion+1][ion] );
                 }
                 fprintf( ioQQQ, "\n" );

branches/newmole/source/iso_level.cpp

@@ -71 +71 @@
         /* now do the 2D array */
         z.alloc(numlevels_local,numlevels_local);
+        z.zero();
 
         /* >>chng 03 may 01, move defn of these two up here, before entering matrix */
@@ -235 +236 @@
                         /* all process depopulating level and placing into the continuum
                          * this does NOT include grain charge transfer ionization, added below */
-                        z[level][level] = iso.RateLevel2Cont[ipISO][nelem][level];
+                        z[level][level] += iso.RateLevel2Cont[ipISO][nelem][level];
                         
                         if( ipISO == ipHE_LIKE && level == ipHe2s3S )
@@ -271 +272 @@
 
                                 z[ipLo][ipLo] += coll_up + pump ;
-                                z[ipLo][level] = - ( coll_up + pump );
+                                z[ipLo][level] -= coll_up + pump ;
 
                                 double pump_down = pump *
@@ -278 +279 @@
 
                                 z[level][level] += RadDecay + pump_down + coll_down;
-                                z[level][ipLo] = - (RadDecay + pump_down + coll_down);
+                                z[level][ipLo] -= RadDecay + pump_down + coll_down;
 
                                 if( ipISO == ipHE_LIKE )

branches/newmole/source/mole_co_atom.cpp

@@ -333 +333 @@
                         (*Rotate)[j].Emis->ColOvTot = 0.0;
                 else
-                        (*Rotate)[j].Emis->ColOvTot = (realnum)(CollRate[j][j+1]/(CollRate[j][j+1]+(*Rotate)[j].Emis->pump) );
+                        (*Rotate)[j].Emis->ColOvTot = (realnum)(CollRate[j][j+1]/SDIV(CollRate[j][j+1]+(*Rotate)[j].Emis->pump) );
 
                 /* two cases - collisionally excited (usual case) or 

branches/newmole/source/mole_eval_balance.cpp

@@ -94 +94 @@
                 
                 for(i=0;i<rate->nreactants;i++)
-                {
+                {       
+                        sp = rate->reactants[i];
                         if (rate->rvector[i] == NULL)
                         {
-                                sp = rate->reactants[i];
                                 mole.zone[sp->index].snk += rate_deriv[i];
                                 if (sp == debug_species) 
@@ -112 +112 @@
                                 for( nelem=0; nelem< LIMELM; ++nelem )
                                 {
-                                        if (rate->reactants[i]->nElem[nelem] != 0)
+                                        if (sp->nElem[nelem] != 0)
                                         {
-                                                mole.xMoleChTrRate[nelem][rate->reactants[i]->charge][rate->rvector[i]->charge] += rate_deriv[i];
+                                                mole.xMoleChTrRate[nelem][sp->charge][rate->rvector[i]->charge] += rate_deriv[i];
                                                 break;
                                         }
@@ -123 +123 @@
                 for(i=0;i<rate->nproducts;i++)
                 {
+                        sp = rate->products[i];
                         if (rate->pvector[i] == NULL)
                         {
-                                sp = rate->products[i];
                                 mole.zone[sp->index].src += rate_tot;
                                 if (sp == debug_species)