Changeset 1851 for branches/newmole/source/ion_solver.cpp

Timestamp:

03/16/08 03:11:06 (10 months ago)

Author:

rjrw

Message:

Deal with FPE in dense.cpp

Use iteratively improved matrix solver in ion_solver

Change criteria for effectively homogeneous matrix in ion_solver

Files:

• branches/newmole/source/ion_solver.cpp (modified) (13 diffs)

branches/newmole/source/ion_solver.cpp

@@ -32 +32 @@
         bool lgPrintIt)
 {
-        /* use tridiag or general matrix solver? */
-        bool lgTriDiag=false;
         long int ion, 
           limit, 
@@ -40 +38 @@
           ns,
           jmax=-1;
-        double totsrc;
+        double totsrc,
+                totscl;
         double rateone;
         bool lgNegPop;
-        static bool lgMustMalloc=true;
-        static double *sink, *source , *sourcesave;
         int32 nerror;
-        static int32 *ipiv;
-        long ion_low, ion_range, i, j, ion_to , ion_from;
+        long ion_low, i, j, ion_to , ion_from;
         /* only used for debug printout */
-        static double *auger;
+        const long ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
 
         double abund_total, renormnew;
         bool lgHomogeneous = true;
-        static double *xmat , *xmatsave;
+        valarray<double> xmat(ion_range*ion_range),xmatsave(ion_range*ion_range);
+        valarray<double> sink(ion_range), source(ion_range) , sourcesave(ion_range);
+        valarray<int32> ipiv(ion_range);
+        valarray<double> auger(LIMELM+1);
 
         DEBUG_ENTRY( "ion_solver()" );
@@ -125 +124 @@
 
         /* the full range of ionization - this is number of ionization stages */
-        ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
         ASSERT( ion_range <= nelem+2 );
 
@@ -133 +131 @@
          * running in the same memory pool - we must have fresh instances of the 
          * arrays for each */
-        if( PARALLEL_MODE || lgMustMalloc )
-        {
-                long int ncell=LIMELM+1;
-                lgMustMalloc = false;
-                if( PARALLEL_MODE )
-                        ncell = ion_range;
-
-                /* this will be "new" matrix, with non-adjacent coupling included */
-                xmat=(double*)MALLOC( (sizeof(double)*(unsigned)(ncell*ncell) ));
-                xmatsave=(double*)MALLOC( (sizeof(double)*(unsigned)(ncell*ncell) ));
-                sink=(double*)MALLOC( (sizeof(double)*(unsigned)ncell ));
-                source=(double*)MALLOC( (sizeof(double)*(unsigned)ncell ));
-                sourcesave=(double*)MALLOC( (sizeof(double)*(unsigned)ncell ));
-                ipiv=(int32*)MALLOC( (sizeof(int32)*(unsigned)ncell ));
-                /* auger is used only for debug printout - it is special because with multi-electron
-                 * Auger ejection, very high stages of ionization can be produced, well beyond
-                 * the highest stage considered here.  so we malloc to the limit */
-                auger=(double*)MALLOC( (sizeof(double)*(unsigned)(LIMELM+1) ));
-        }
+
+        /* auger is used only for debug printout - it is special because with multi-electron
+         * Auger ejection, very high stages of ionization can be produced, well beyond
+         * the highest stage considered here.  so we malloc to the limit */
 
         for( i=0; i<ion_range;i++ )
@@ -160 +143 @@
         /* this will be used to address the 2d arrays */
 #       undef MAT
-#       define MAT(M_,I_,J_)    (*((M_)+(I_)*(ion_range)+(J_)))
+#       define MAT(M_,I_,J_)    ((M_)[(I_)*(ion_range)+(J_)])
 
         /* zero-out loop comes before main loop since there are off-diagonal
@@ -198 +181 @@
                         dense.IonHigh[nelem] = 3;
                         abund_total = 1.;
-                        ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
                         /* make up ionization and recombination rates */
                         for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
@@ -377 +359 @@
         /* this will be sum of source and sink terms, will be used to decide if
          * matrix is singular */
-        totsrc = 0.;
+        totsrc = totscl = 0.;
 
         for( i=0; i<ion_range;i++ )
@@ -387 +369 @@
                 source[i] -= mole.source[nelem][ion];
                 totsrc += mole.source[nelem][ion];
+                totscl += -MAT( xmat, i, i )*dense.xIonDense[nelem][ion];
                 MAT( xmat, i, i ) -= mole.sink[nelem][ion];
         }
 
         /* matrix is not homogeneous if source is non-zero */
-        if( totsrc != 0. )
+        if( totsrc > 1e-10*totscl )
                 lgHomogeneous = false;
 
@@ -530 +513 @@
                         }
                         fprintf(ioQQQ,"\n");
-                }
-        }
-
-        /* first branch is tridiag, following branch is just general matrix solver */
-        if( lgTriDiag ) 
-        {
-                bool lgLapack=false , lgTriOptimized=true;
-                /* there are two choices for the tridiag solver */
-                if(lgLapack) 
-                {
-                        /* this branch uses lapack tridiag solver, and should work 
-                         * it is hardwired to not be used because not extensively tested
-                         * issues - is this slower than others, and is it robust in
-                         * singular cases? */
-                        bool lgBad = false;
-                        /* Use Lapack tridiagonal solver */
-                        double *dl, *d, *du;
-
-                        d = (double *) MALLOC((unsigned)ion_range*sizeof(double));
-                        du = (double *) MALLOC((unsigned)(ion_range-1)*sizeof(double));
-                        dl = (double *) MALLOC((unsigned)(ion_range-1)*sizeof(double));
-
-                        for(i=0;i<ion_range-1;i++) 
-                        {
-                                du[i] = MAT(xmat,i+1,i);
-                                d[i] = MAT(xmat,i,i);
-                                dl[i] = MAT(xmat,i,i+1);
-                        }
-                        d[i] = MAT(xmat,i,i);
-
-#                       if 0
-                        int i1, i2;
-                        i1 = ion_range;
-                        i2 = 1;
-                        lgBad = dgtsv_wrapper(&i1, &i2, dl, d, du, source, &i2);
-#                       endif
-                        if( lgBad )
-                                fprintf(ioQQQ," dgtsz error.\n");
-
-                        free(dl);free(du);free(d);
-                } 
-                else if(lgTriOptimized)
-                {
-                        /* Use tridiagonal solver re-coded to avoid rounding errors
-                         * on diagonal -- uses determination of jmax for the
-                         * singular case, but is otherwise independent of the array
-                         * filling code above */
-
-                        for(i=0;i<ion_range;i++) 
-                        {
-                                source[i] = sink[i] = 0.;
-                        }
-                        for(i=0;i<ion_range;i++) 
-                        {
-                                ion = i+ion_low;
-                                source[i] += mole.source[nelem][ion];
-                                sink[i] += mole.sink[nelem][ion];
-                        }
-                        if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection && radius.depth < dynamics.oldFullDepth )
-                        {
-                                for(i=0;i<ion_range;i++) 
-                                {
-                                        ion = i+ion_low;
-                                        source[i] += dynamics.Source[nelem][ion];
-                                        sink[i] += dynamics.Rate;
-                                }
-                        }
-
-                        solveions(ionbal.RateIonizTot[nelem]+ion_low,ionbal.RateRecomTot[nelem]+ion_low,
-                                                                sink,source,ion_range,jmax);
-                }
-        } 
+                        fprintf(ioQQQ,"totsrc/totscl %g %g\n",totsrc,totscl);
+                }
+        }
+
+        if (1)
+        {
+                int32 solve_system(const valarray<double> &a, valarray<double> &b, const long int n);
+                nerror = solve_system(xmat,source,ion_range);
+        }
         else
         {
@@ -607 +527 @@
                 nerror = 0;
                 /* Use general matrix solver */
-                getrf_wrapper(ion_range, ion_range, xmat, ion_range, ipiv, &nerror);
+                getrf_wrapper(ion_range, ion_range, &xmat[0], ion_range, &ipiv[0], &nerror);
                 if( nerror != 0 )
                 {
                         fprintf( ioQQQ, 
-                                " DISASTER ion_solver: dgetrf finds singular or ill-conditioned matrix nelem=%li %s ion_range=%li, limit=%li, nConv %li xmat follows\n",
-                                nelem , 
-                                elementnames.chElementSym[nelem],
-                                ion_range,
-                                limit , 
-                                conv.nTotalIoniz );
+                                                         " DISASTER ion_solver: dgetrf finds singular or ill-conditioned matrix nelem=%li %s ion_range=%li, limit=%li, nConv %li xmat follows\n",
+                                                         nelem , 
+                                                         elementnames.chElementSym[nelem],
+                                                         ion_range,
+                                                         limit , 
+                                                         conv.nTotalIoniz );
                         for( i=0; i<ion_range; ++i )
                         {
@@ -633 +553 @@
                         cdEXIT(EXIT_FAILURE);
                 }
-                getrs_wrapper('N', ion_range, 1, xmat, ion_range, ipiv, source, ion_range, &nerror);
+                getrs_wrapper('N', ion_range, 1, &xmat[0], ion_range, &ipiv[0], &source[0], ion_range, &nerror);
                 if( nerror != 0 )
                 {
                         fprintf( ioQQQ, " DISASTER ion_solver: dgetrs finds singular or ill-conditioned matrix nelem=%li ionrange=%li\n",
-                                nelem , ion_range );
+                                                         nelem , ion_range );
                         cdEXIT(EXIT_FAILURE);
                 }
@@ -694 +614 @@
                         test += source[ion-ion_low];
                 }
-                test += source[ion-ion_low];
         }
 
@@ -1126 +1045 @@
         }
 
-        if( PARALLEL_MODE )
-        {
-                free(ipiv);
-                free(source);
-                free(sink);
-                free(xmat);
-                free(xmatsave);
-                free(auger);
-        }
         return;
 }