Changeset 2094 for branches/newmole

Timestamp:

05/28/08 14:24:08 (7 months ago)

Author:

rjrw

Message:

Fix for problems in

pdr_coolbb.in -- step control in newton_step.cpp

blr_nf84.in -- have to switch off NeH+ chemistry

Location:

branches/newmole/source

Files:

• mole_species.cpp (modified) (1 diff)
• newton_step.cpp (modified) (7 diffs)

branches/newmole/source/mole_species.cpp

@@ -229 +229 @@
         sp = newspecies("ClO+  ",MOLECULE,MOLE_ACTIVE,NULL, 1158.0f);
         sp = newspecies("ArH+  ",MOLECULE,MOLE_ACTIVE,NULL, 0.0f)/* Binding Energy ?? */;
-        sp = newspecies("NeH+  ",MOLECULE,MOLE_ACTIVE,NULL, 0.0f)/* Binding Energy ?? */;
+        if (0) // Including NeH+ breaks blr_nf84.in ???
+                sp = newspecies("NeH+  ",MOLECULE,MOLE_ACTIVE,NULL, 0.0f)/* Binding Energy ?? */;
         if(gv.lgDustOn && mole.lgGrain_mole_deplete )
         {               

branches/newmole/source/newton_step.cpp

@@ -42 +42 @@
                 error1,
                 error0,
-                grad;
+                grad,
+                pred;
 
         valarray<double> 
@@ -58 +59 @@
         }
 
-        error0 = error2 = grad = f1 = f2 = 0.;
+        pred = error0 = error2 = grad = f1 = f2 = 0.;
 
         for (loop=0;loop<40;loop++) 
@@ -77 +78 @@
                 {
                         /* Smooth the error mode tailoff */
-                        etmp = ervals[i]/((ERRLIM+fabs(b0vec[i]))*(ERRLIM+fabs(escale[i])));
+                        etmp = ervals[i]/(1e-24+fabs(b0vec[i]*escale[i]));
                         etmp *= etmp;
                         if (etmp > emax)
@@ -93 +94 @@
                                 b1vec[i] = ervals[i];
                         }
-                        //fprintf(ioQQQ,"Chem ");
                         merror = solve_system(amat,b1vec,n);
                         error0 = error1;
                         grad = -2*error0;
                         f1 = 1.0;
-
+                        if (error0 > 1e20) // Limit step for really bad initial conditions
+                                f1 = 1e-3;
+                        if (0)
+                                for ( i=0; i < n; i++ )
+                                {
+                                        if (b2vec[i] > 0. && f1*fabs(b1vec[i]) > 0.1*b2vec[i])
+                                        {
+                                                f1 = 0.1*b2vec[i]/fabs(b1vec[i]);
+                                        }
+                                }
+                        
                         if (0)
                         {
@@ -144 +154 @@
                         }
                 } else {
-                        //fprintf(ioQQQ,"Newt %ld stp %11.4g err %11.4g grd %11.4g fdg %11.4g e0 %11.4g de %11.4g\n",
-                        //                              loop,f1,error1,grad,(error1-error0)/f1,error0,error1-error0);
+                        //fprintf(ioQQQ,"Newt %ld stp %11.4g err %11.4g grd %11.4g fdg %11.4g e0 %11.4g de %11.4g pred %11.4g\n",
+                        //                              loop,f1,error1,grad,(error1-error0)/f1,error0,error1-error0,pred);
                         if (error1 < (1-2e-4*f1)*error0 || error1 < 1e-20)
                         {
@@ -154 +164 @@
                         {
                                 f2 = f1;
-                                f1 = -0.5*f1*grad/SDIV(error1-error0-f1*grad);
+                                f1 *= -0.5*f1*grad/(error1-error0-f1*grad);
+                                pred = error0+0.5*grad*f1;
                         }
                         else
@@ -163 +174 @@
                                 double b = (-f2*a1/(f1*f1)+f1*b1/(f2*f2));
                                 double ft = f1;
-                                f1 = (-b+sqrt(b*b-3.*a*grad*(f1-f2)))/(3.*a);
+                                //fprintf(ioQQQ,"Pred 1 %g %g\n",error1,error0+ft*(grad+ft/(ft-f2)*(b+a*ft)));
+                                //fprintf(ioQQQ,"Pred 2 %g %g\n",error2,error0+f2*(grad+f2/(ft-f2)*(b+a*f2)));                          
+                                //f1 = (-b+sqrt(b*b-3.*a*grad*(ft-f2)))/(3.*a);
+                                //fprintf(ioQQQ,"Pred x %g %g\n",f1,error0+f1*(grad+f1/(ft-f2)*(b+a*f1)));
+                                f1 = 1.-3.*(a/b)*(grad/b)*(ft-f2);
+                                if (f1 > 0.) 
+                                        f1 = b/(3.*a)*(-1.-sqrt(f1));
+                                else
+                                        f1 = -b/(3.*a);
+                                //fprintf(ioQQQ,"Pred n %g %g\n",f1,error0+f1*(grad+f1/(ft-f2)*(b+a*f1)));                              
+                                pred = error0+f1*(grad+f1/(ft-f2)*(b+a*f1));
                                 f2 = ft;
                         }