Changeset 2118 for branches/newmole/source/newton_step.cpp

Timestamp:

06/22/08 10:51:55 (5 months ago)

Author:

rjrw

Message:

Use internal timestep limiting (...in an attempt...) to improve
robustness of molecular solver.

Files:

• branches/newmole/source/newton_step.cpp (modified) (10 diffs)

branches/newmole/source/newton_step.cpp

@@ -36 +36 @@
         double 
                 etmp,
+                etmpeq,
                 emax,
                 f1,
@@ -42 +43 @@
                 error1,
                 error0,
+                erroreq,
                 grad,
-                pred;
+                pred,
+                errstp;
 
         valarray<double> 
@@ -49 +52 @@
                 b1vec(n),
                 escale(n),
-                ervals(n);
+                ervals(n),
+                ervals1(n);
         int32 merror;
 
@@ -59 +63 @@
         }
 
-        pred = error0 = error2 = grad = f1 = f2 = 0.;
+        pred = error0 = error2 = erroreq = grad = f1 = f2 = 0.;
+        errstp = 0.0;
 
         for (loop=0;loop<40;loop++) 
@@ -67 +72 @@
                 if (loop == 0) 
                 {
+                        // Limit timestep to value where linearizations are reasonable
+                        iworst = -1;
                         for( i=0; i < n; i++ )
                         {
                                 escale[i] = MAT(amat,i,i);
-                        }
-                }
-
+                                if (b0vec[i] > 0 && errstp*b0vec[i] < fabs(ervals[i]))
+                                {
+                                        errstp = fabs(ervals[i])/(1e-24+b0vec[i]);
+                                        iworst = i;
+                                }
+                        }
+                        errstp *= 3.;
+                        for (i=0; i < n; i++)
+                        { 
+                                MAT(amat,i,i) -= errstp;
+                        }
+                        //if (iworst != -1)
+                        //      fprintf(stdout,"Timescale %g rate %g worst %s %g %g %g\n",
+                        //                                      1./SDIV(errstp),errstp,groupspecies[iworst]->label,
+                        //                                      ervals[iworst],b0vec[iworst],escale[iworst]);
+                }
+
+                // Error from fully relaxed equilibrium solution
+                erroreq = 0.;
+                for( i=0; i < n; i++ )
+                {
+                        etmpeq = ervals[i];//(1e-24+fabs(b0vec[i]*escale[i]));
+                        erroreq += etmpeq*etmpeq;
+                }
+
+                // Error from timestep-limited solution
+                for (i=0; i < n; i++)
+                {
+                        ervals1[i] = ervals[i]-errstp*(b2vec[i]-b0vec[i]);
+                }
                 error1 = emax = 0.;
                 iworst = -1;
@@ -78 +112 @@
                 {
                         /* Smooth the error mode tailoff */
-                        etmp = ervals[i]/(1e-24+fabs(b0vec[i]*escale[i]));
+                        etmp = ervals1[i];//(1e-24+fabs(b0vec[i]*escale[i]));
                         etmp *= etmp;
                         if (etmp > emax)
@@ -92 +126 @@
                         for( i=0; i < n; i++ )
                         {
-                                b1vec[i] = ervals[i];
+                                b1vec[i] = ervals1[i];
                         }
                         merror = solve_system(amat,b1vec,n);
@@ -230 +264 @@
                                 b1vec[j] = b0vec[j];
                         }
-                        jacobn(b1vec,&ervals[0],&amat[0],false);
+                        jacobn(b1vec,&ervals1[0],&amat[0],false);
                         for( i=0; i < n; i++ )
                         {
@@ -244 +278 @@
                                 {
                                         double e1 = MAT(amat,i,j);
-                                        double e2 = (escale[j]-ervals[j])/db;
+                                        double e2 = (escale[j]-ervals1[j])/db;
                                         if (fabs(e1-e2) > 0.01*fabs(e1+e2))
                                                 fprintf(ioQQQ,"%7s %7s %11.4g %11.4g %11.4g\n",
                                                                                 groupspecies[i]->label,groupspecies[j]->label,e1,e2,
-                                                                                ervals[j]/db);
+                                                                                ervals1[j]/db);
                                 }
                         }
@@ -255 +289 @@
         }
 
-        *error = (realnum) MIN2(error1,1e30);
+        *error = (realnum) MIN2(erroreq,1e30);
 
         return;