Changeset 2069 for branches/newmole

Timestamp:

05/15/08 13:33:22 (8 months ago)

Author:

rjrw

Message:

More optional debugging code,

Ensure we're using the exact Jacobian for constraints,

Remove exit(-1) from default behaviour when > 40 loops (probably still
fails elsewhere, though).

Location:

branches/newmole/source

Files:

• mole_solve.cpp (modified) (9 diffs)
• newton_step.cpp (modified) (7 diffs)

branches/newmole/source/mole_solve.cpp

@@ -49 +49 @@
 #endif
 
+static double molnum[LIMELM];
+
 double mole_solve()
 {
@@ -114 +116 @@
                                 fjac(mole.num_compacted*mole.num_compacted),
                                 wa(lwa);
-                        double tol = 1e-7;
+                        double tol = 1e-9;
                         for ( long i=0; i < mole.num_compacted; i++ )
                         {
@@ -121 +123 @@
                         hybrj1_(fcn, &mole.num_compacted, &newmols[0], &fvec[0], &fjac[0],
                                                         &mole.num_compacted, &tol, &info, &wa[0], &lwa);
-                        fprintf(ioQQQ,"hybrj1_ returns info %d\n",info);
+                        //fprintf(ioQQQ,"hybrj1_ returns info %d\n",info);
 #else
                         ;
@@ -325 +327 @@
         }
         if (*iflag == 1)
+        {
                 funjac(b2vec, fvec, &tmpvec[0], false);
+        }
         else if (*iflag == 2)
-                funjac(b2vec, &tmpvec[0], fjac, true);
+        {
+                funjac(b2vec, &tmpvec[0], &fjac[0], true);
+        }
 }
 #endif
@@ -495 +501 @@
                         // constraints, now look for common species
                         valarray<int> iamax(mole.num_compacted), ncons(LIMELM);
+                        double molnow[LIMELM];
+
+                        for ( nelem=ipHYDROGEN; nelem<LIMELM; nelem++)
+                        {
+                                molnow[nelem] = 0.;
+                        }
+                        for( i=0; i < mole.num_compacted; i++ )
+                        {
+                                for ( nelem=ipHYDROGEN; nelem<LIMELM; nelem++)
+                                {
+                                        molnow[nelem] += b2vec[i]*groupspecies[i]->nElem[nelem];
+                                }
+                        }
                         iamax = 0;
                         for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
@@ -500 +519 @@
                                 if (element_list[nelem]->ipMl[0] != -1)
                                 {
-                                        int imax = -1;
-                                        double amax = 0.;
-                                        for (i=0;i<mole.num_compacted;i++)
-                                        {
-                                                if (0 == iamax[i] && groupspecies[i]->nElem[nelem]*groupspecies[i]->den > amax)
+                                        if (0)
+                                        {
+                                                // determine common species
+                                                int imax = -1;
+                                                double amax = 0.;
+                                                for (i=0;i<mole.num_compacted;i++)
                                                 {
-                                                        amax = groupspecies[i]->nElem[nelem]*groupspecies[i]->den;
-                                                        imax = i;
+                                                        if (0 == iamax[i] && groupspecies[i]->nElem[nelem]*groupspecies[i]->den > amax)
+                                                        {
+                                                                amax = groupspecies[i]->nElem[nelem]*groupspecies[i]->den;
+                                                                imax = i;
+                                                        }
                                                 }
-                                        }
-                                        iamax[imax] = 1;
-                                        ncons[nelem] = groupspecies[imax]->index; // was element_list[nelem]->ipMl[0];
+                                                iamax[imax] = 1;
+                                                
+                                                ncons[nelem] = groupspecies[imax]->index;
+                                        }
+                                        else
+                                        {
+                                                ncons[nelem] = element_list[nelem]->ipMl[0];
+                                        }
                                 }
                         }
@@ -519 +547 @@
                                 if (element_list[nelem]->ipMl[0] != -1)
                                 {
+                                        double scale=1.;
                                         if (lgJac)
                                         {
-                                                double scale = -(ABSLIM+fabs(c[ncons[nelem]][ncons[nelem]]));
+                                                //scale = -(ABSLIM+fabs(c[ncons[nelem]][ncons[nelem]]));
                                                 for(i=0;i<mole.num_compacted;i++)
                                                 {
@@ -527 +556 @@
                                                 }
                                         }
-                                        
-                                        b[ncons[nelem]] = 0.;
+                                        b[ncons[nelem]] = (molnow[nelem]-molnum[nelem]);
+                                        //fprintf(ioQQQ,"Elem %ld err %g\n",nelem,b[ncons[nelem]]);
                                 }
                         }
@@ -773 +802 @@
         }       
 
+        for ( nelem=ipHYDROGEN; nelem<LIMELM; nelem++)
+        {
+                molnum[nelem] = 0.;
+        }
         for( i=0; i < mole.num_compacted; i++ )
         {
                 b0vec[i] = calcv[groupspecies[i]->index];
+                for ( nelem=ipHYDROGEN; nelem<LIMELM; nelem++)
+                {
+                        molnum[nelem] += b0vec[i]*groupspecies[i]->nElem[nelem];
+                }
         }
 }

branches/newmole/source/newton_step.cpp

@@ -8 +8 @@
 #include "cddefines.h"
 #include "thirdparty.h"
+#include "mole.h"
+#include "mole_priv.h"
 
 #define ABSLIM  1e-12
@@ -22 +24 @@
 /* mole_newton_step -- improve balance in chemical network along
  * descent direction, step limited to ensure improvement */
-void newton_step(const valarray<double> &b0vec, valarray<double> &b2vec, realnum *error, const long n,
-                        void (*jacobn)(const valarray<double> &b2vec, double * const ervals, double * const amat, const bool lgJac))
+void newton_step(const valarray<double> &b0vec, valarray<double> &b2vec, 
+                                                                 realnum *error, const long n,
+                                                                 void (*jacobn)(const valarray<double> &b2vec, 
+                                                                                                                                double * const ervals, double * const amat, 
+                                                                                                                                const bool lgJac))
 {
         long int i, 
@@ -55 +60 @@
         error0 = error2 = grad = f1 = f2 = 0.;
 
-        for (loop=0;loop<20;loop++) 
+        for (loop=0;loop<40;loop++) 
         {
                 jacobn(b2vec,&ervals[0],&amat[0],loop==0);
@@ -92 +97 @@
                         grad = -2*error0;
                         f1 = 1.0;
+
+                        if (0)
+                        {
+                                for (i=0; i<LIMELM; ++i)
+                                {
+                                        double felem = 0., fproj = 0.;
+                                        for (long j=0; j<n; ++j)
+                                        {
+                                                felem += b1vec[j]*groupspecies[j]->nElem[i];
+                                                fproj += groupspecies[j]->nElem[i];
+                                        }
+                                        //fprintf(ioQQQ,"Elem %ld cons %g\n",i,felem);
+                                        if (felem != 0.)
+                                        {
+                                                felem /= fproj;
+                                                for (long j=0; j<n; ++j)
+                                                {
+                                                        b1vec[j] -= groupspecies[j]->nElem[i]*felem;
+                                                }
+                                        }
+                                }
+                        }
+                        if (0)
+                        {
+                                double rmax = 0., rval;
+                                int imax=0;
+                                for( i=0; i < n; ++i )
+                                {
+                                        if (b0vec[i] != 0.)
+                                                rval = b1vec[i]/b0vec[i];
+                                        else
+                                                rval = 0.;
+                                        if (fabs(rval) > rmax) 
+                                        {
+                                                rmax = fabs(rval);
+                                                imax = i;
+                                        }
+                                }
+                                fprintf(ioQQQ,"Biggest is %s %g\n",groupspecies[imax]->label,
+                                                                b1vec[imax]/b0vec[imax]);
+                                /* if (rmax > 0.5)
+                                         f1 = 0.5/rmax;
+                                         if (f1 < 1e-3)
+                                         f1 = 1e-3; */
+                        }
                 } else {
+                        //fprintf(ioQQQ,"Newt %ld stp %11.4g err %11.4g grd %11.4g fdg %11.4g e0 %11.4g de %11.4g\n",
+                        //                              loop,f1,error1,grad,(error1-error0)/f1,error0,error1-error0);
                         if (error1 < (1-2e-4*f1)*error0 || error1 < 1e-20)
                         {
@@ -100 +152 @@
                         if (loop == 1) 
                         {
-                                f2 = 1.;
-                                f1 = -0.5*grad/SDIV(error1-error0-grad);
+                                f2 = f1;
+                                f1 = -0.5*f1*grad/SDIV(error1-error0-f1*grad);
                         }
                         else
@@ -117 +169 @@
                                 f1 = 0.5*f2;
                         else if (f1 < 0.1*f2)
-                                f1 = 0.1*f2;
-                }
-
+                                f1 = 0.1*f2;                    
+                }
+                
                 /*
                  * b1vec[] is descent direction
@@ -129 +181 @@
                         b2vec[i] = b0vec[i]-b1vec[i]*f1;                                
                 }
+        }
+        if (0 && 40 == loop)
+        {
+                double rmax = 0., rval;
+                int imax=0;
+                for( i=0; i < n; ++i )
+                {
+                        if (b0vec[i] != 0.)
+                                rval = b1vec[i]/b0vec[i];
+                        else
+                                rval = 0.;
+                        fprintf(ioQQQ,"%7s %11.4g %11.4g\n",
+                                                        groupspecies[i]->label,rval,b0vec[i]);
+                        if (fabs(rval) > rmax) 
+                        {
+                                rmax = fabs(rval);
+                                imax = i;
+                        }
+                }
+                fprintf(ioQQQ,"Biggest is %s\n",groupspecies[imax]->label);
+                if (0)
+                { // Verify Jacobian
+                        long j;
+                        for( j=0; j < n; j++ )
+                        {
+                                b1vec[j] = b0vec[j];
+                        }
+                        jacobn(b1vec,&ervals[0],&amat[0],false);
+                        for( i=0; i < n; i++ )
+                        {
+                                for( j=0; j < n; j++ )
+                                {
+                                        b1vec[j] = b0vec[j];
+                                }
+                                double db = 1e-3*fabs(b0vec[i])+1e-9;
+                                b1vec[i] += db;
+                                db = b1vec[i]-b0vec[i];
+                                jacobn(b1vec,&escale[0],&amat[0],false);
+                                for( j=0; j < n; j++ )
+                                {
+                                        double e1 = MAT(amat,i,j);
+                                        double e2 = (escale[j]-ervals[j])/db;
+                                        if (fabs(e1-e2) > 0.01*fabs(e1+e2))
+                                                fprintf(ioQQQ,"%7s %7s %11.4g %11.4g %11.4g\n",
+                                                                                groupspecies[i]->label,groupspecies[j]->label,e1,e2,
+                                                                                ervals[j]/db);
+                                }
+                        }
+                }
+                exit(-1);
         }