Changeset 2031 for branches/newmole/source/mole_solve.cpp

Timestamp:

05/07/08 16:29:10 (8 months ago)

Author:

rjrw

Message:

Stubs for using minpack routines (they don't improve things, possibly
due to the Jacobian fix for species conservation).

Re-order a loop in conv_init_solution.

Remove some out of date comments.

Files:

• branches/newmole/source/mole_solve.cpp (modified) (8 diffs)

branches/newmole/source/mole_solve.cpp

@@ -40 +40 @@
 
 #define SMALLABUND 1e-24
+#ifdef MINPACK
+#include <minpack.h>
+STATIC void fcn(const int *n, const double *x, double *fvec, double *fjac,
+                                                                const int *ldfjac, int *iflag );
+#endif
 
 double mole_solve()
@@ -95 +100 @@
 
                 MoleMap.setup(&oldmols[0]);
-                newton_step(oldmols, newmols, &error, mole.num_compacted, funjac);
+                if (1)
+                {
+                        newton_step(oldmols, newmols, &error, mole.num_compacted, funjac); 
+                }
+                else
+                {
+#ifdef MINPACK
+                        int info, lwa = (mole.num_compacted*(mole.num_compacted+13))/2;
+                        valarray<double> fvec(mole.num_compacted),
+                                fjac(mole.num_compacted*mole.num_compacted),
+                                wa(lwa);
+                        double tol = 1e-7;
+                        for ( long i=0; i < mole.num_compacted; i++ )
+                        {
+                                newmols[i] = oldmols[i];
+                        }
+                        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);
+#else
+                        ;
+#endif
+                }
 
                 /* check for negative populations */
@@ -117 +144 @@
                 MoleMap.updateMolecules( &newmols[0] );
                 
-                //fprintf(stdout,"Mole zone %ld -- %7.2f loop %d error %15.8g nBad %d\n",nzone,fnzone,i,error,nBad);
+                // fprintf(stdout,"Mole zone %ld -- %7.2f error %15.8g nBad %d\n",nzone,fnzone,error,nBad);
         
                 if (error < BIGERROR && nBad == 0 && nPrevBad == 0) /* Stop early if good enough */
@@ -283 +310 @@
 GroupMap MoleMap;
 
+#ifdef MINPACK
+STATIC void fcn(const int *n, const double *x, double *fvec, double *fjac,
+                                                                const int *ldfjac, int *iflag )
+{
+        /* Function to interface funjac with MINPACK routines */
+        ASSERT (*ldfjac == *n);
+        valarray<double> b2vec(*n), tmpvec(*n);
+        for (long i=0; i < *n; i++) 
+        {
+                b2vec[i] = x[i];
+        }
+        if (*iflag == 1)
+                funjac(b2vec, fvec, &tmpvec[0], false);
+        else if (*iflag == 2)
+                funjac(b2vec, &tmpvec[0], fjac, true);
+}
+#endif
+
 STATIC void funjac(const valarray<double> &b2vec, double * const ervals, double * const amat, const bool lgJac)
 {
@@ -334 +379 @@
         
         if (lgJac) {
-                // c[i][i] can be positive for auto-catalytic reactions
-                //for(i=0;i<mole.num_calc;i++) 
-                //{
-                //      ASSERT (c[i][i] <= 0.);
-                //}
                 for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
                 {
@@ -401 +441 @@
         if (lgConserve)
         {
+                if (0)
+                {
+                        double sum1;
+                        for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
+                        {
+                                sum1 = 0.;
+                                for(i=0;i<mole.num_calc;i++)
+                                {
+                                        sum1 += b[i]*mole.list[i]->nElem[nelem];
+                                }
+                                fprintf(ioQQQ,"Elem %ld cons %g %g\n",nelem,sum1,dense.gas_phase[nelem]);
+                                if (lgJac)
+                                {
+                                        for(i=0;i<mole.num_calc;i++)
+                                        {
+                                                sum1 = 0.;
+                                                for (long j=0;j<mole.num_calc;j++)
+                                                {
+                                                        sum1 += c[i][j]*mole.list[j]->nElem[nelem];
+                                                }
+                                                fprintf(ioQQQ,"Eqn %ld error %g\n",i,sum1);
+                                        }
+                                }
+                        }
+                }
+
                 double scale;
                 for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
@@ -415 +481 @@
                                 }
                                 else
+                                {
                                         scale = c[element_list[nelem]->ipMl[0]][element_list[nelem]->ipMl[0]];
-                                
+                                }                               
+
                                 b[element_list[nelem]->ipMl[0]] = 0.;
                         }
@@ -459 +527 @@
                 for( i=0; i < mole.num_compacted; i++ )
                 {
-                        // c[i][i] can be +ve for (effectively) autocatalytic reactions,
-                        // e.g. H-,H+=>H,H
-                        //ASSERT (c[groupspecies[i]->index][groupspecies[i]->index] <= 0.);
                         for( j=0; j < mole.num_compacted; j++ )
                         {