Changeset 1834 for branches/newmole/source

Timestamp:

03/09/08 07:08:13 (9 months ago)

Author:

rjrw

Message:

Merged from trunk r1830:1833

Files:

• branches/newmole/source/iso_level.cpp (modified) (10 diffs)

branches/newmole/source/iso_level.cpp

@@ -50 +50 @@
         bool lgNegPop=false;
         static double SaveHe23S_photorate=0.;
-        static int32 *ipiv; /* malloc out to [numlevels_local] */
+        valarray<int32> ipiv(numlevels_local); 
         /* this block of variables will be obtained and freed here */
-        static double
-                *creation,
-                *error/*[numlevels_local]*/,
-                *work/*[numlevels_local]*/,
-                *Save_creation;
-        static bool lgSpaceMalloc=false;
-        double *PopPerN;
+        valarray<double>
+                creation(numlevels_local),
+                error(numlevels_local),
+                work(numlevels_local),
+                Save_creation(numlevels_local);
+        double source=0.,
+                sink=0.;
+        valarray<double> PopPerN(iso.n_HighestResolved_local[ipISO][nelem]+1);
 
         multi_arr<double,2,C_TYPE> z, SaveZ;
-
-        /* this flag true means to malloc and free all space each time */
-#       define PARALLEL false
 
         DEBUG_ENTRY( "iso_level()" );
@@ -70 +68 @@
         ASSERT( nelem >= ipISO );
         ASSERT( nelem < LIMELM );
-
-        /* malloc some scratch space */
-        if( PARALLEL )
-        {
-                ipiv = (int32 *)MALLOC(sizeof(int32)*(unsigned)(numlevels_local) );
-                creation = (double *)MALLOC(sizeof(double)*(unsigned)(numlevels_local) );
-                error = (double *)MALLOC(sizeof(double)*(unsigned)(numlevels_local) );
-                work = (double *)MALLOC(sizeof(double)*(unsigned)(numlevels_local) );
-                Save_creation = (double *)MALLOC(sizeof(double)*(unsigned)(numlevels_local) );
-        }
-        else if( !lgSpaceMalloc )
-        {
-                /* space has not been malloced yet, but will only do it one time */
-                lgSpaceMalloc = true;
-
-                ipiv = (int32 *)MALLOC(sizeof(int32)*(unsigned)(max_num_levels) );
-                creation = (double *)MALLOC(sizeof(double)*(unsigned)(max_num_levels) );
-                error = (double *)MALLOC(sizeof(double)*(unsigned)(max_num_levels) );
-                work = (double *)MALLOC(sizeof(double)*(unsigned)(max_num_levels) );
-                Save_creation = (double *)MALLOC(sizeof(double)*(unsigned)(max_num_levels) );
-        }
 
         /* now do the 2D array */
@@ -370 +347 @@
                         {
                                 /* 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]);
+                                source += gv.GrainChTrRate[nelem][ion][nelem-ipISO] * 
+                                        dense.xIonDense[nelem][ion] ;
                                 /* take ionization out of every level. */
-                                for( level=0; level < numlevels_local; level++ )
-                                {
-                                        z[level][level] += gv.GrainChTrRate[nelem][nelem-ipISO][ion];
-                                }
+                                sink += gv.GrainChTrRate[nelem][nelem-ipISO][ion];
                         }
                 }
@@ -382 +356 @@
                 /* >>chng 02 Sep 06 rjrw -- all elements have these terms */
                 /*>>>chng 02 oct 01, only include if lgAdvection is set */
-                if( dynamics.lgAdvection )
+                if( dynamics.lgAdvection && dynamics.lgISO[ipISO])
                 {
                         /* add in advection - these terms normally zero */
                         /* assume for present that all advection is into ground state */
-                        creation[0] += 
-                                dynamics.Source[nelem][nelem-ipISO]/
-                                SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*
-                                dynamics.lgISO[ipISO];
+                        source += dynamics.Source[nelem][nelem-ipISO];
                         /* >>chng 02 Sep 06 rjrw -- advective term not recombination */
                         /* can sink from all components (must do, for conservation) */
-                        for( ipLo=0; ipLo<numlevels_local; ++ipLo )
-                        {
-                                z[ipLo][ipLo] += dynamics.Rate*dynamics.lgISO[ipISO];
-                        }
+                        sink += dynamics.Rate;
                 }
 
@@ -401 +369 @@
                 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];
-                        }
+                        source += mole.source[nelem][nelem-ipISO];
+                        sink += mole.sink[nelem][nelem-ipISO];
                         
                         for( long ion=0; ion<N_MOLE_ION; ++ion )
@@ -412 +377 @@
                                 {
                                         /* 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];
-                                        }
+                                        source += mole.xMoleChTrRate[nelem][ion][nelem-ipISO] * 
+                                                dense.xIonDense[nelem][ion] ;
+                                        sink += mole.xMoleChTrRate[nelem][nelem-ipISO][ion];
                                 }
                         }
+                }
+
+                creation[0] += source/SDIV(dense.xIonDense[nelem][nelem+1-ipISO]);
+                for( level=0; level < numlevels_local; level++ )
+                {
+                        z[level][level] += sink;
                 }
 
@@ -494 +462 @@
 
                 getrf_wrapper(numlevels_local,numlevels_local,
-                              z.data(),numlevels_local,ipiv,&nerror);
-
-                getrs_wrapper('N',numlevels_local,1,z.data(),numlevels_local,ipiv,
-                              creation,numlevels_local,&nerror);
+                              z.data(),numlevels_local,&ipiv[0],&nerror);
+
+                getrs_wrapper('N',numlevels_local,1,z.data(),numlevels_local,&ipiv[0],
+                              &creation[0],numlevels_local,&nerror);
 
                 if( nerror != 0 )
@@ -744 +712 @@
          * of highest to next highest stage of ionization */
         {
-                PopPerN = (double *)MALLOC(sizeof(double)*(unsigned)(iso.n_HighestResolved_local[ipISO][nelem]+1) );
 
                 for( long n=0; n <= iso.n_HighestResolved_local[ipISO][nelem]; n++ )
@@ -809 +776 @@
                 }
         }
-
-        free( PopPerN );
 
         /************************************************/
@@ -896 +861 @@
         }
 #endif
-
-        /* only release this if option for parallel processing */
-        if( PARALLEL )
-        {
-                free( Save_creation );
-                free( work );
-                free( error );
-                free( creation );
-                free( ipiv );
-        }
 
         return;