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

Timestamp:

05/06/08 16:35:36 (8 months ago)

Author:

rjrw

Message:

Further refactor newton_step to more purely numerical function.

Remove more dead or leftover code and data.

Files:

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

branches/newmole/source/mole_solve.cpp

@@ -35 +35 @@
 void check_co_ion_converge(void);
 
+STATIC void funjac(const valarray<double> &b2vec, double * const ervals, double * const amat, const bool lgJac);
+
 STATIC void mole_h_fixup(void);
+
+#define SMALLABUND 1e-24
 
 double mole_solve()
@@ -42 +46 @@
         realnum
                 error;
-        valarray<double> b2vec(mole.num_compacted);
+        valarray<double> oldmols(mole.num_compacted),
+                newmols(mole.num_compacted);
 
         DEBUG_ENTRY( "mole_solve()" );
@@ -89 +94 @@
                 nPrevBad = nBad;
 
-                MoleMap.setup(&b2vec[0]);
-                newton_step(b2vec, &nBad, &error, mole.num_compacted);
-                MoleMap.updateMolecules( &b2vec[0] );
+                MoleMap.setup(&oldmols[0]);
+                newton_step(oldmols, newmols, &error, mole.num_compacted, funjac);
+
+                /* check for negative populations */
+                nBad = 0; 
+                long fracneg = 0.;
+                long iworst = -1;
+                for( long i=0; i < mole.num_compacted; i++ )
+                {
+                        if( newmols[i] < 0. ) 
+                        {
+                                if(-newmols[i] > fracneg*SDIV(oldmols[i]) && oldmols[i] > SMALLABUND) 
+                                {
+                                        fracneg = -newmols[i]/oldmols[i];
+                                        iworst = i;
+                                }
+                                if (oldmols[i] > SMALLABUND || newmols[i] < -SMALLABUND)
+                                        ++nBad;
+                        }
+                }
+                
+                MoleMap.updateMolecules( &newmols[0] );
                 
                 //fprintf(stdout,"Mole zone %ld -- %7.2f loop %d error %15.8g nBad %d\n",nzone,fnzone,i,error,nBad);
@@ -248 +272 @@
 #define ABSLIM  1e-12
 #define ERRLIM  1e-12
-#define SMALLABUND 1e-24
 #       ifdef MAT
 #               undef MAT
@@ -260 +283 @@
 GroupMap MoleMap;
 
-void funjac(valarray<double> &b2vec, double *ervals, double *amat, bool lgJac)
+STATIC void funjac(const valarray<double> &b2vec, double * const ervals, double * const amat, const bool lgJac)
 {
         long
@@ -520 +543 @@
 }
 
-void GroupMap::updateMolecules(double *b2)
+void GroupMap::updateMolecules(const double * const b2)
 {
         long int