Changeset 2031 for branches/newmole/source

Timestamp:

05/07/08 16:29:10 (7 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.

Location:

branches/newmole/source

Files:

• conv_init_solution.cpp (modified) (2 diffs)
• mole_eval_balance.cpp (modified) (1 diff)
• mole_solve.cpp (modified) (8 diffs)
• newton_step.cpp (modified) (1 diff)

branches/newmole/source/conv_init_solution.cpp

@@ -130 +130 @@
         for( i=0; i<mole.num_calc; ++i )
         {
-                for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++) 
-                {
-                        if( dense.lgElmtOn[nelem] && mole.list[i]->location == NULL )
+                if (mole.list[i]->location == NULL)
+                {
+                        for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++) 
                         {
-                                realnum dens_elemsp = (realnum) mole.list[i]->den*mole.list[i]->nElem[nelem];
-                                if (    FracMoleMax*dense.gas_phase[nelem] < dens_elemsp ) 
+                                if( dense.lgElmtOn[nelem] )
                                 {
-                                        FracMoleMax = dens_elemsp/dense.gas_phase[nelem];
-                                        ipMax = i;
+                                        realnum dens_elemsp = (realnum) mole.list[i]->den*mole.list[i]->nElem[nelem];
+                                        if ( FracMoleMax*dense.gas_phase[nelem] < dens_elemsp ) 
+                                        {
+                                                FracMoleMax = dens_elemsp/dense.gas_phase[nelem];
+                                                ipMax = i;
+                                        }
                                 }
                         }
@@ -148 +151 @@
         for(i=0;i<mole.num_calc;++i)
         {
-                OxyInGrains += (1 - mole.list[i]->lgGas_Phase)*mole.list[i]->den*mole.list[i]->nElem[ipOXYGEN]; 
+                if (! mole.list[i]->lgGas_Phase )
+                        OxyInGrains += mole.list[i]->den*mole.list[i]->nElem[ipOXYGEN]; 
         }
         /* this is now fraction of O in ices */

branches/newmole/source/mole_eval_balance.cpp

@@ -64 +64 @@
         }
 
-        /* Call all the routines that set up the matrix 
-         * CO_solve will call this routine, therefore all other matrix elements are
-         * included here so that, when CO_solve is called, everything is accounted for */       
-
-        /* All now cross-validated with new treatment, switching causes only v. minor
-         * differences in results */
-        /* Revised molecular network implementation */
-        /* Fills matrix elements for passive species -- these can be used to
-                 derive sources & sinks resulting from this part of the network */
         for(mole_reaction_i p=mole_priv.reactab.begin(); 
                         p != mole_priv.reactab.end(); p++) 

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++ )
                         {

branches/newmole/source/newton_step.cpp

@@ -54 +54 @@
         error0 = 1;
         error1 = grad = 0.;
-        fstep = flast = 1.;                     
+        fstep = 1.;                     
 
         for (loop=0;loop<20;loop++)