Changeset 2185

Timestamp:

07/09/08 12:05:44 (6 months ago)

Author:

peter

Message:

source/conv_itercheck.cpp:

Use safe way to check for convergence (see PR73).

Files:

• trunk/source/conv_itercheck.cpp (modified) (5 diffs)

trunk/source/conv_itercheck.cpp

@@ -70 +70 @@
                                         {
                                                 /* test if Lya converged, nLyaLevel is upper level of Lya for iso seq */
-                                                if( fabs(Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauTot/
-                                                        (Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn*rt.DoubleTau)-1.) > 
-                                                        conv.autocv )
+                                                if( fabs(Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauTot -
+                                                         Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn*rt.DoubleTau) > 
+                                                    conv.autocv*fabs(Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn*rt.DoubleTau) )
                                                 {
                                                         /* not converged to within AUTOCV, normally 15 percent */
@@ -93 +93 @@
                                                 }
 
-                                                if( fabs(Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot/
-                                                        (Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn*rt.DoubleTau)-1.) > conv.autocv )
+                                                if( fabs(Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot -
+                                                         Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn*rt.DoubleTau) >
+                                                    conv.autocv*fabs(Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn*rt.DoubleTau) )
                                                 {
                                                         /* not converged to within AUTOCV, normally 15 percent */
@@ -122 +123 @@
                 for( i=0; i<NCOLD; ++i )
                 {
-                        double differ = fabs(colden.colden_old[i]-colden.colden[i]) /
-                                SDIV(colden.colden[i]);
                         /* was the species column density significant relative to
                          * the total H column density, and was its abundance changing? */
-                        if( (colden.colden[i]/colden.colden[ipCOL_HTOT] > 1e-5) &&
-                                (differ > conv.autocv) )
+                        if( colden.colden[i]/colden.colden[ipCOL_HTOT] > 1e-5 &&
+                            fabs(colden.colden_old[i]-colden.colden[i]) > conv.autocv*colden.colden[i] )
                         {
                                 /* not converged to within conv.autocv, normally 20 percent */
@@ -151 +150 @@
                 for( i=0; i<mole.num_comole_calc; ++i )
                 {
-                        double differ;
                         if(COmole[i]->n_nuclei == 1)
                                 continue;
 
                         /* was the species abundance and changing? */
-                        differ = (double)fabs(COmole[i]->hevcol_old-COmole[i]->hevcol) /
-                                (double)SDIV(COmole[i]->hevcol);
-                        if( (COmole[i]->hevcol/colden.colden[ipCOL_HTOT] > 1e-5) &&
-                                (differ > conv.autocv) )
+                        if( COmole[i]->hevcol/colden.colden[ipCOL_HTOT] > 1e-5 &&
+                            fabs(COmole[i]->hevcol_old-COmole[i]->hevcol) > conv.autocv*COmole[i]->hevcol )
                         {
                                 /* not converged to within conv.autocv, normally 20 percent */
@@ -182 +178 @@
                 if( dynamics.lgAdvection )
                 {
-                        double error1 = dynamics.convergence_error / SDIV(dynamics.error_scale2);
-                        double error2 = dynamics.discretization_error / SDIV(dynamics.error_scale2);
                         /* >>chng 02 nov 29, as per Will Henney email */
-                        if( MAX2(error1/dynamics.convergence_tolerance, error2 ) > conv.autocv )
-                        /*if( MAX2(error1, error2 ) > conv.autocv )*/
+                        if( dynamics.convergence_error > conv.autocv*dynamics.error_scale2*dynamics.convergence_tolerance ||
+                            dynamics.discretization_error > conv.autocv*dynamics.error_scale2 )
                         {
                                 lgConverged = false;