Changeset 1676 for branches/newmole/source

Timestamp:

12/15/07 08:24:13 (11 months ago)

Author:

rjrw

Message:

Back out changes to mole_newton_step from r1645:1652 which seemed to
cause more problems than they solved.

Some finesse of the handling of output buffering.

Location:

branches/newmole/source

Files:

• cddrive.cpp (modified) (1 diff)
• mole.h (modified) (1 diff)
• mole_co_solve.cpp (modified) (7 diffs)
• mole_newton_step.cpp (modified) (18 diffs)
• prt_final.cpp (modified) (1 diff)

branches/newmole/source/cddrive.cpp

@@ -2102 +2102 @@
                         fflush( ioQQQ );
                         /* stderr is always unbuffered */
-                        ioQQQ = stderr;
+                        //ioQQQ = stderr;
+                        setvbuf( ioQQQ, NULL, _IONBF, 0); 
                         /* will be used to generate comment at end */
-                        input.lgSetNoBuffering = false;
+                        input.lgSetNoBuffering = true;
                 }
         }

branches/newmole/source/mole.h

@@ -42 +42 @@
 
 /** Take one Newton step of the chemical network 
-\param *nFixup
+\param *nBad
 \param *error
 */
-void mole_newton_step(int *nFixup, realnum *error, realnum numelem[]);
+void mole_newton_step(int *nBad, realnum *error, realnum numelem[]);
 
 /** \verbatim >>chng 03 feb 09, rm ipH3P_hev, since not used, and decrement NUM_HEAVY_MOLEC to 17 

branches/newmole/source/mole_co_solve.cpp

@@ -36 +36 @@
 static void mole_h_fixup(void);
 /**hmole determine populations of hydrogen molecules */
-STATIC double hmole(bool *lgNegPop, bool *lgZerPop);
+STATIC double hmole( void );
 
 
@@ -62 +62 @@
         ASSERT(lgElemsConserved());
 
-        error = hmole(lgNegPop, lgZerPop);              
+        error = hmole();                
 
         ASSERT(lgElemsConserved());
@@ -196 +196 @@
  * For conditions of distribution and use see copyright notice in license.txt */
 /*hmole determine populations of hydrogen molecules */
-STATIC double hmole( bool *lgNegPop, bool *lgZerPop )
+STATIC double hmole( void )
 {
-        int nFixup, i, nelem;
+        int nBad, i, nelem;
         realnum error,
                 numelem[LIMELM],
@@ -227 +227 @@
         /* loop until run to limit, or no fix ups needed and error is small */
         /* will be used to keep track of neg solns */
-        nFixup = 0;
+        nBad = 0;
         error = -1.;
         for(i=0; i < LIM_LOOP;i++) 
@@ -247 +247 @@
                         }
                 }
-                /* nFixup is number of times negative abundances were fixed, should be zero 
-                 * at return for valid soln */
-                nFixup = 0;
-
-                mole_newton_step(&nFixup, &error, numelem);
+                /* nBad returns number of times negative abundances were fixed
+                 * for latest step, should be zero at return for valid soln */
+
+                mole_newton_step(&nBad, &error, numelem);
                 
-                //fprintf(stdout,"Hmole zone %ld -- %7.2f loop %d error %15.8g fixup %d\n",nzone,fnzone,i,error,nFixup);
+                //fprintf(stdout,"Hmole zone %ld -- %7.2f loop %d error %15.8g fixup %d\n",nzone,fnzone,i,error,nBad);
 
                 {
@@ -284 +283 @@
                         }
                 }
-                if (error < BIGERROR && !nFixup) /* Stop early if good enough */
+
+                if (error < BIGERROR && nBad == 0) /* Stop early if good enough */
                         break;
         }
 
-        if( (i == LIM_LOOP && error > BIGERROR)  || nFixup != 0 )
+        if( (i == LIM_LOOP && error > BIGERROR)  || nBad != 0 )
         {
                 /* most of these failures occur just one time during search phase -
@@ -294 +294 @@
                 conv.lgConvPops = false;
 
-                if( !conv.lgSearch && called.lgTalk )
+                if( 1||!conv.lgSearch && called.lgTalk )
                 {
                         fprintf(ioQQQ," PROBLEM  hmole, zone %li: %d iters, %d bad; final error: %g lgSearch %i\n",
                         nzone,
                         i,
-                        nFixup,
+                        nBad,
                         error,
                         conv.lgSearch);

branches/newmole/source/mole_newton_step.cpp

@@ -40 +40 @@
 /* mole_newton_step -- improve balance in chemical network along
  * descent direction, step limited to ensure improvement */
-void mole_newton_step(int *nFixup, realnum *error, realnum numelem[])
+void mole_newton_step(int *nBad, realnum *error, realnum numelem[])
 {
         enum {PRINTSOL = false};
@@ -76 +76 @@
                 *calcv=NULL,
                 *b=NULL,
-                **c=NULL,
-                *scalev;
+                **c=NULL;
         int32 merror;
         double maxrate = 0.,
@@ -96 +95 @@
                 b1vec = ((double*)MALLOC( (size_t)(mole.num_compacted)*sizeof(double) ));
                 b2vec = ((double*)MALLOC( (size_t)(mole.num_compacted)*sizeof(double) ));
-                scalev = ((double*)MALLOC( (size_t)(mole.num_compacted)*sizeof(double) ));
                 if (mole.num_calc != 0)
                 {
@@ -180 +178 @@
         error1 = grad = 0.;
         fstep = flast = 1.;                     
-        *nFixup = 0; 
+        *nBad = 0; 
 
         lgNegPop = false;
@@ -196 +194 @@
                 
                 /* add positive ions and neutral atoms: ratios are set by ion_solver,
-                 * we determine abundance of the group as a whole here */
+                 * we determine sum of two here */
                 
                 if (loop == 0) {
@@ -237 +235 @@
                                 }
                         }
-                        maxrate = 0.0;
-                        for(i=0;i<mole.num_compacted;i++)
-                        {
-                                if (-c[groupspecies[i]->index][groupspecies[i]->index] > maxrate)
-                                        maxrate = -c[groupspecies[i]->index][groupspecies[i]->index];
-                        }
-                        for(i=0;i<mole.num_compacted;i++) 
-                        {
-                                //ASSERT (c[COmole[i]->index][COmole[i]->index] < 0);
-                                scalev[groupspecies[i]->index] = 1./(1e-6*maxrate-c[groupspecies[i]->index][groupspecies[i]->index]);
-                                for (j=0;j<mole.num_compacted;j++) 
-                                {
-                                        c[groupspecies[j]->index][groupspecies[i]->index] *= scalev[groupspecies[i]->index];
-                                }
-                        }
                 }
 
@@ -271 +254 @@
                 }
 
-                for (i=0;i<mole.num_compacted;++i) 
-                {
-                        b[groupspecies[i]->index] *= scalev[groupspecies[i]->index];
-                }
-
                 /* For Newton-Raphson method, want the change in populations to be zero,
                  * so the conserved component must also be zero */
@@ -308 +286 @@
                 for( i=0; i < mole.num_compacted; i++ )
                 {
-                        etmp = (b[groupspecies[i]->index]-scalerate*(b2vec[i]-b0vec[i]));
+                        /* Smooth the error mode tailoff */
+                        etmp = (b[groupspecies[i]->index]-scalerate*(b2vec[i]-b0vec[i]))/
+                                ((ABSLIM+b0vec[i])*(-c[groupspecies[i]->index][groupspecies[i]->index]));
                         etmp *= etmp;
                         if (etmp > emax)
@@ -321 +301 @@
                 }
 
-                //fprintf(stdout,"Loop %ld error %15.8g fstep %15.8g; worst species %s abund %g error %g\n",
-                //                              loop,error0,fstep,groupspecies[iworst]->label,groupspecies[iworst]->hevmol,emax);
+                // fprintf(stdout,"Loop %ld error %15.8g fstep %15.8g; worst species %s abund %g  error %g\n",
+                //                              loop,error0,fstep,groupspecies[iworst]->label,groupspecies[iworst]->hevmol,,emax);
 
                 /*------------------------------------------------------------------ */
@@ -397 +377 @@
                                                  groupspecies[i]->hevmol,
                                                  -b1vec[i]); */
+                                        if(b0vec[i] > 1e-12)
+                                        {
+                                                fixit(); /* Nadger low populations -- magic number */
                                                 
-                                        if (b0vec[i]-b1vec[i] < 0) 
-                                                fprintf(stdout,"Bad NR step for %s wants %g step %g orig %g\n",groupspecies[i]->label,b0vec[i]-b1vec[i],b1vec[i],b0vec[i]);
-                                        if (fstep*(1e-24+b1vec[i]) > 0.9*(1e-24+b0vec[i])) 
-                                        {
-                                                fstep = 0.9*(1e-24+b0vec[i])/(1e-24+b1vec[i]);
-                                                iworst = i;
-                                                fprintf(stdout,"Limiting %s %g %g %g\n",groupspecies[i]->label,b0vec[i],b1vec[i],fstep);
-                                        }                               
-                                }
-                                // if (fstep < 0.1*flast)
-                                //      fstep = 0.1*flast;                      
+                                                if (b0vec[i]-b1vec[i] < 0) 
+                                                        fprintf(stdout,"Bad NR step for %s %g %g\n",groupspecies[i]->label,b0vec[i],b1vec[i]);
+                                                if (fstep*b1vec[i] > 0.9*(ABSLIM+b0vec[i])) 
+                                                {
+                                                        fstep = 0.9*(ABSLIM+b0vec[i])/(b1vec[i]);
+                                                        iworst = i;
+                                                }
+                                        }
+                                
+                                }
+                                if (fstep < 0.1*flast)
+                                        fstep = 0.1*flast;                      
                         }
                         //fprintf(ioQQQ,"Step factor %g bad species %s abund %g step %g\n",
@@ -420 +404 @@
                         //                              iworst!=-1?-b1vec[iworst]:0.); 
                         //fprintf(stdout,"Grads %g %g (%g %g)\n",grad,(error0-error1)/fstep,error0,error1);
-                        if ((error0 < (1-2e-4*fstep)*error1 || error0 < 1e-20))
+                        if (!lgNegPop && (error0 < (1-2e-4*fstep)*error1 || error0 < 1e-20))
                         {
                                 // fprintf(stdout,"OK %g %g\n",error0,error1);
                                 break;
-                        }
-                        else
-                        {
-                                // fprintf(ioQQQ,"Continuing: tests %c %g %g\n",TorF(lgNegPop),error0,(1-2e-4*fstep)*error1-error0);
                         }
                         fstep = -0.5*grad*fstep*fstep/SDIV(error0-error1-grad*fstep);
@@ -518 +498 @@
                 /* check for negative populations */
                 lgNegPop = false;
-                *nFixup = 0; 
+                *nBad = 0; 
                 fracneg = 0.;
                 iworst = -1;
@@ -525 +505 @@
                         if( b2vec[i] < 0. ) 
                         {
-                                //fprintf(stdout,"Species %s negative at %g from %g\n",groupspecies[i]->label,b2vec[i],b0vec[i]);
+                                // fprintf(stdout,"Species %s negative at %g from %g\n",groupspecies[i]->label,b2vec[i],b0vec[i]);
                                 /* largest relative change in solution for neg soln */
                                 /* this can only occur for negative solutions since fracneg starts
@@ -534 +514 @@
                                         iworst = i;
                                 }
-                                if (b2vec[i] < -SMALLABUND)
-                                        lgNegPop = true;
-                                b2vec[i] = 0.;
-                        }
-                }
+                                //if (b2vec[i] < -SMALLABUND)
+                                //      lgNegPop = true;
+                                //b2vec[i] = 0.;
+                        }
+                }
+                
                 
 #       if 0
@@ -566 +547 @@
                         {
                                 if (b0vec[i] > SMALLABUND)
-                                        ++*nFixup;
+                                        ++*nBad;
                         }
                 }
@@ -578 +559 @@
                                 {
                                         if (b0vec[i] > SMALLABUND)
-                                                ++*nFixup;
+                                                ++*nBad;
                                         if (0)
                                                 b2vec[i] = 0.;
@@ -586 +567 @@
                         /* very common to obtain negative solution on very first try - 
                          * don't print in this case */
-                        if(0 && *nFixup > 0)
+                        if(*nBad > 0)
                         {
                                 fprintf( ioQQQ, " PROBLEM  hmole_step finds negative H molecule, in zone %.2f, loop %ld step %g.\n",fnzone,loop,fstep );
@@ -882 +863 @@
         if ( merror != 0 )
         {
-                fprintf( ioQQQ, " CO_solve getrf_wrapper error %ld\n",merror );
+                fprintf( ioQQQ, " CO_solve getrf_wrapper error %ld\n",(long int) merror );
                 puts( "[Stop in CO_solve]" );
                 cdEXIT(EXIT_FAILURE);

branches/newmole/source/prt_final.cpp

@@ -1940 +1940 @@
         /* print input title for model */
         fprintf( ioQQQ, "%s\n\n", input.chTitle );
+        fflush(ioQQQ);
         return;
 }