Changeset 2056

Timestamp:

05/13/08 17:51:39 (4 months ago)

Author:

rjrw

Message:

Move to overwriting common species with conservation constraints
rather than always the atomic species, to try to improve robustness
using standard scheme (doesn't appear to have much effect).

Files:

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

branches/newmole/source/mole_solve.cpp

@@ -464 +464 @@
                         {
                                 double sum1 = 0.;
-                                for(i=0;i<mole.num_calc;i++)
-                                {
-                                        sum1 += b[i]*mole.list[i]->nElem[nelem];
+                                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]);
@@ -474 +474 @@
                                         {
                                                 sum1 = 0.;
-                                                for (long j=0;j<mole.num_calc;j++)
+                                                for (long j=0;j<mole.num_compacted;j++)
                                                 {
-                                                        sum1 += c[i][j]*mole.list[j]->nElem[nelem];
+                                                        sum1 += c[groupspecies[i]->index][groupspecies[j]->index]*
+                                                                groupspecies[j]->nElem[nelem];
                                                 }
                                                 fprintf(ioQQQ,"Eqn %ld error %g\n",i,sum1);
@@ -491 +492 @@
                 if (lgConserve)
                 {
-                        // Original scheme -- replace atom rows with conservation constraints
+                        // Original scheme -- replace atom rows with conservation
+                        // constraints, now look for common species
+                        valarray<int> iamax(mole.num_compacted), ncons(LIMELM);
+                        iamax = 0;
+                        for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
+                        {
+                                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)
+                                                {
+                                                        amax = groupspecies[i]->nElem[nelem]*groupspecies[i]->den;
+                                                        imax = i;
+                                                }
+                                        }
+                                        iamax[imax] = 1;
+                                        ncons[nelem] = groupspecies[imax]->index; // was element_list[nelem]->ipMl[0];
+                                }
+                        }
+
                         double scale;
                         for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
@@ -499 +522 @@
                                         if (lgJac)
                                         {
-                                                scale = -(ABSLIM+fabs(c[element_list[nelem]->ipMl[0]][element_list[nelem]->ipMl[0]]));
-                                                for(i=0;i<mole.num_calc;i++)
+                                                scale = -(ABSLIM+fabs(c[ncons[nelem]][ncons[nelem]]));
+                                                for(i=0;i<mole.num_compacted;i++)
                                                 {
-                                                        c[i][element_list[nelem]->ipMl[0]] = mole.list[i]->nElem[nelem]*scale;
+                                                        c[groupspecies[i]->index][ncons[nelem]] = groupspecies[i]->nElem[nelem]*scale;
                                                 }
                                         }
-                                        else
-                                        {
-                                                scale = c[element_list[nelem]->ipMl[0]][element_list[nelem]->ipMl[0]];
-                                        }                               
                                         
-                                        b[element_list[nelem]->ipMl[0]] = 0.;
+                                        b[ncons[nelem]] = 0.;
                                 }
                         }