Changeset 1804 for branches/newmole/source

Timestamp:

02/19/08 15:23:09 (9 months ago)

Author:

rjrw

Message:

Merge from trunk r1800:1803

Correct iso_level handling of molecular sources, along lines of Ryan's
fix to the grain charge transfer.

Location:

branches/newmole/source

Files:

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

branches/newmole/source/ion_solver.cpp

@@ -246 +246 @@
                 low = dense.IonLow[nelem];
 
-                for( ion_to=low; ion_to <= limit; ion_to++ )
+                for( ion_to=low; ion_to <= dense.IonHigh[nelem]; ion_to++ )
                 {
                         for( ion_from=dense.IonLow[nelem]; ion_from <= dense.IonHigh[nelem]; ++ion_from )

branches/newmole/source/iso_ionize_recombine.cpp

@@ -7 +7 @@
 #include "atmdat.h"
 #include "dense.h"
-#include "grainvar.h"
 #include "hmi.h"
 #include "mole.h"
@@ -115 +114 @@
         long int nelem )
 {
-        long int level,
-                ion;
-        double Recom3Body,
-                sum;
+        long int level;
+        double Recom3Body;
 
         DEBUG_ENTRY( "iso_ionize_recombine()" );
@@ -167 +164 @@
         /* all following go into or out of ground state */
         level = 0;
-
-        /* grain charge transfer recombination and ionization, 
-         * assume goes into and comes from ground state */
-        /* first grain surface losses - both to lower and higher stages of ionization */
-        sum = 0.;
-        for(ion=0; ion<nelem+1; ++ion )
-                if( ion!=nelem-ipISO )
-                        sum += gv.GrainChTrRate[nelem][nelem-ipISO][ion];
-        /*>>chng 04 sep 08, replaced grain ct rate from old sum with current sum */
-        iso.RateLevel2Cont[ipISO][nelem][level] += sum;
-
-        /* >>chng 04 apr 10, add formally correct grain surface ionization and recombination - 
-         * this HAD NOT been included for H-like species */
-        /* >>chng 04 sep 08, use explicit term from higher stage of ionization - note that added
-         * GrainCreat[nelem][nelem-ipISO] for case of atomic He incorrectly added two-stage recombination
-         * to existing matrix logic */
-        /*iso.RateCont2Level[ipISO][nelem][level] += ionbal.GrainCreat[nelem][nelem-ipISO];*/
-        iso.RateCont2Level[ipISO][nelem][level] += gv.GrainChTrRate[nelem][nelem+1-ipISO][nelem-ipISO];
-
-        for(ion=0; ion<nelem+1; ++ion)
-                if ( ion < nelem-ipISO )
-                        iso.RateLevel2Cont[ipISO][nelem][level] += mole.xMoleChTrRate[nelem][nelem-ipISO][ion];
-                else if ( ion > nelem-ipISO )
-                        iso.RateCont2Level[ipISO][nelem][level] += mole.xMoleChTrRate[nelem][ion][nelem-ipISO];
 
         /* now charge transfer - all into/from ground, two cases, H and not H */

branches/newmole/source/iso_level.cpp

@@ -9 +9 @@
 #include "dynamics.h"
 #include "elementnames.h"
+#include "grainvar.h"
 #include "he.h"
 #include "helike.h"
@@ -246 +247 @@
                 {
                         /* all process depopulating level and placing into the continuum
-                         * this includes grain charge transfer ionization */
+                         * this does NOT include grain charge transfer ionization, added below */
                         z[level][level] = iso.RateLevel2Cont[ipISO][nelem][level];
-
+                        
                         if( ipISO == ipHE_LIKE && level == ipHe2s3S )
                                 /* >>chng 05 dec 21, rm eden to make into rate coefficient */
@@ -353 +354 @@
                 z[1+ipISO][1+ipISO] += iso.TwoNu_induc_dn[ipISO][nelem]*iso.lgInd2nu_On;
 
+                /* grain charge transfer recombination and ionization to ALL other stages */
+                for( long ion=0; ion<=nelem+1; ++ion )
+                {
+                        if( ion!=nelem-ipISO )
+                        {
+                                /* recombination must be multiplied by a ratio of densities to get proper rate. */
+                                creation[0] += gv.GrainChTrRate[nelem][ion][nelem-ipISO] * 
+                                        dense.xIonDense[nelem][ion] / SDIV(dense.xIonDense[nelem][nelem+1-ipISO]);
+                                z[0][0] += gv.GrainChTrRate[nelem][nelem-ipISO][ion];
+                        }
+                }
+
                 /* >>chng 02 Sep 06 rjrw -- all elements have these terms */
                 /*>>>chng 02 oct 01, only include if lgAdvection is set */
@@ -372 +385 @@
 
                 /* add in source and sink terms from molecular network. */
-                if( nelem == ipISO && conv.nTotalIoniz )
-                {
-                        /* these are the external source and sink terms */
-                        /* source first */
+                if( conv.nTotalIoniz && nelem-ipISO < N_MOLE_ION)
+                {
                         creation[0] += mole.source[nelem][nelem-ipISO]/SDIV(dense.xIonDense[nelem][nelem+1-ipISO]);
-
                         for( ipLo=0; ipLo<numlevels_local; ++ipLo )
                         {
                                 z[ipLo][ipLo] += mole.sink[nelem][nelem-ipISO];
+                        }
+                        
+                        for( long ion=0; ion<N_MOLE_ION; ++ion )
+                        {
+                                if( ion!=nelem-ipISO && ion < nelem+1 )
+                                {
+                                        /* recombination must be multiplied by a ratio of densities to get proper rate. */
+                                        creation[0] += mole.xMoleChTrRate[nelem][ion][nelem-ipISO] * 
+                                                dense.xIonDense[nelem][ion] / SDIV(dense.xIonDense[nelem][nelem+1-ipISO]);
+                                        for( ipLo=0; ipLo<numlevels_local; ++ipLo )
+                                        {
+                                                z[ipLo][ipLo] += mole.xMoleChTrRate[nelem][nelem-ipISO][ion];
+                                        }
+                                }
                         }
                 }