Changeset 2426 for branches/newmole

Timestamp:

10/27/08 15:50:47 (2 months ago)

Author:

rjrw

Message:

Fix for limit_lte_he1_ste.in -- increase size of initial guess on
rlimit in newton_step.cpp. This may well need further tuning to
balance robustness and speed (or generally an improved approach).

Also add template to deal with problems on VS/Sun compilers getting
pointer-to-const from const valarray.

Location:

branches/newmole/source

Files:

• mole.h (modified) (1 diff)
• mole_solve.cpp (modified) (11 diffs)
• newton_step.cpp (modified) (1 diff)

branches/newmole/source/mole.h

@@ -192 +192 @@
 extern struct chem_element_s *element_list[LIMELM];
 
+template<class T> const T* get_const_address(const valarray<T> v)
+{
+        return const_cast<const T*>(&const_cast<valarray<T>&>(v)[0]);
+}
+
 #endif /* _MOLE_H_ */

branches/newmole/source/mole_solve.cpp

@@ -46 +46 @@
 #define SMALLABUND 1e-24
 
-static double molnum[LIMELM];
+static double molElems[LIMELM],molElemsNew[LIMELM];
 
 double mole_solve()
@@ -106 +106 @@
 
                 MoleMap.setup(&oldmols[0]);
+
+                grouped_elems(&oldmols[0],molElems);
                 
                 long j;
@@ -113 +115 @@
                         newton_step(oldmols, newmols, &error, mole.num_compacted, &rlimit, escale, funjac); 
                         
+                        if (0)
+                        {
+                                grouped_elems(&newmols[0],molElemsNew);
+                                for (long nelem=0; nelem <LIMELM; nelem++)
+                                {
+                                        if (molElems[nelem] > 0.0)
+                                        {
+                                                fprintf(ioQQQ,"StepChg %ld %g %g\n",nelem,molElems[nelem],molElemsNew[nelem]);
+                                        }
+                                }
+                        }
+                
                         /* check for negative populations */
                         nBad = 0; 
@@ -418 +432 @@
                 }
         }
+
+        if (0)
+        {
+                for (nelem=ipHYDROGEN;nelem<=ipHELIUM;nelem++)
+                {
+                        double sum1 = 0., sum2 = 0.;
+                        for(i=0;i<mole.num_calc;i++)
+                        {
+                                sum1 += b[i]*mole.list[i]->nElem[nelem];
+                                sum2 += fabs( b[i]*mole.list[i]->nElem[nelem] );
+                        }
+                        fprintf(ioQQQ,"Elem ungroup %ld cons %g %g %g\n",nelem,sum1,dense.gas_phase[nelem],sum2);
+                }
+        }
         
         for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
@@ -436 +464 @@
         }
         
+        if (0)
+        {
+                for (nelem=ipHYDROGEN;nelem<=ipHELIUM;nelem++)
+                {
+                        double sum1 = 0.;
+                        for(i=0;i<mole.num_calc;i++)
+                        {
+                                sum1 += b[i]*mole.list[i]->nElem[nelem];
+                        }
+                        fprintf(ioQQQ,"Elem group %ld cons %g %g\n",nelem,sum1,dense.gas_phase[nelem]);
+                }
+        }
+
         /* Species are now grouped -- only mole.num_compacted elements now active */
 
         /* For Newton-Raphson method, want the change in populations to be zero,
          * so the conserved component must also be zero */
-        if (lgConserve)
-        {
-                if (0)
-                {
-                        for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
-                        {
-                                double sum1 = 0.;
+        if (0)
+        {
+                for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
+                {
+                        double sum1 = 0.;
+                        for(i=0;i<mole.num_compacted;i++)
+                        {
+                                sum1 += b[groupspecies[i]->index]*groupspecies[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_compacted;i++)
                                 {
-                                        sum1 += b[groupspecies[i]->index]*groupspecies[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_compacted;j++)
-                                                {
-                                                        sum1 += c[groupspecies[i]->index][groupspecies[j]->index]*
-                                                                groupspecies[j]->nElem[nelem];
-                                                }
-                                                fprintf(ioQQQ,"Eqn %ld error %g\n",i,sum1);
-                                        }
+                                        sum1 = 0.;
+                                        for (long j=0;j<mole.num_compacted;j++)
+                                        {
+                                                sum1 += c[groupspecies[i]->index][groupspecies[j]->index]*
+                                                        groupspecies[j]->nElem[nelem];
+                                        }
+                                        fprintf(ioQQQ,"Eqn %ld error %g\n",i,sum1);
                                 }
                         }
@@ -480 +518 @@
                         double molnow[LIMELM];
 
-                        grouped_elems(&b2vec[0], molnow);
-                        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];
-                                }
-                        }
+                        grouped_elems(get_const_address(b2vec), molnow);
+
                         iamax = 0;
                         for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
@@ -534 +562 @@
                                                 }
                                         }
-                                        b[ncons[nelem]] = (molnow[nelem]-molnum[nelem]);
+                                        b[ncons[nelem]] = (molnow[nelem]-molElems[nelem]);
                                         //fprintf(ioQQQ,"Elem %ld err %g\n",nelem,b[ncons[nelem]]);
                                 }
@@ -568 +596 @@
                                 for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
                                 {
-                                        u[i][nelem] = mole.list[groupspecies[i]->index]->nElem[nelem];
+                                        u[i][nelem] = groupspecies[i]->nElem[nelem];
                                 }
                         }
@@ -698 +726 @@
                                 ervals[i] = b[groupspecies[i]->index];
                         }
-                        
+
                         if (lgJac) 
                         {
@@ -711 +739 @@
                 }
         }
+
+        if (0) 
+        {
+                grouped_elems(&ervals[0],molElemsNew);
+                for ( long nelem=0; nelem<LIMELM; nelem++)
+                {
+                        if (molElemsNew[nelem] != 0.0)
+                                fprintf(ioQQQ,"Newmole errs nelem %ld %g\n",nelem,molElemsNew[nelem]);
+                }               
+        }
+
 }
 
@@ -799 +838 @@
                 b0vec[i] = calcv[groupspecies[i]->index];
         }
-        grouped_elems(&b0vec[0],molnum);
 }
 

branches/newmole/source/newton_step.cpp

@@ -86 +86 @@
                                         }
                                 }
-                                *rlimit *= 1e-11;
+                                *rlimit *= 1e-8;
                         }
                         for( i=0; i < n; i++ )