00001
00002
00003
00004
00005
00006 #include "cddefines.h"
00007 #include "dense.h"
00008 #include "trace.h"
00009 #include "thermal.h"
00010 #include "thirdparty.h"
00011 #include "phycon.h"
00012 #include "conv.h"
00013
00014
00015 STATIC bool lgConvEden(void);
00016
00017 STATIC double EdenError(double eden);
00018
00019
00020
00021 int ConvEdenIoniz(void)
00022 {
00023 DEBUG_ENTRY( "ConvEdenIoniz()" );
00024
00025
00026
00027
00028 if( trace.lgTrace )
00029 {
00030 fprintf( ioQQQ, "\n" );
00031 fprintf( ioQQQ, " ConvEdenIoniz entered\n" );
00032 }
00033 if( trace.nTrConvg>=3 )
00034 {
00035 fprintf( ioQQQ,
00036 " ConvEdenIoniz called, entering eden loop using solver %s.\n",
00037 conv.chSolverEden);
00038 }
00039
00040
00041 double EdenEntry = dense.eden;
00042
00043
00044 if( strcmp( conv.chSolverEden , "vWDB" )== 0 )
00045 {
00046 conv.lgConvEden = false;
00047
00048 iter_track NeTrack;
00049 double n1, error1, n2, error2;
00050
00051 for( int n=0; n < 3; ++n )
00052 {
00053 const int DEF_ITER = 10;
00054
00055 const double factor = 0.02;
00056
00057 NeTrack.clear();
00058
00059
00060
00061
00062
00063
00064 n1 = dense.eden;
00065 error1 = EdenError( n1 );
00066 NeTrack.add( n1, error1 );
00067
00068 if( conv.lgSearch && dense.EdenTrue > SMALLFLOAT )
00069 n2 = sqrt(dense.eden*dense.EdenTrue);
00070 else if( abs(safe_div( error1, n1 )) < factor )
00071 n2 = dense.EdenTrue;
00072 else
00073 n2 = ( error1 > 0. ) ? n1*(1.-factor) : n1*(1.+factor);
00074
00075
00076 if( !fp_equal( n1, n2 ) )
00077 error2 = EdenError( n2 );
00078 else
00079 error2 = error1;
00080 NeTrack.add( n2, error2 );
00081
00082 int j = 0;
00083
00084
00085 while( error1*error2 > 0. && j++ < DEF_ITER )
00086 {
00087 n1 = n2;
00088 error1 = error2;
00089 double deriv = NeTrack.deriv(5);
00090
00091 double step = safe_div( -1.2*error1, deriv, 0. );
00092 step = sign( min( abs(step), factor*n1 ), step );
00093 n2 = n1 + step;
00094 error2 = EdenError( n2 );
00095 NeTrack.add( n2, error2 );
00096 }
00097
00098 if( error1*error2 > 0. && trace.nTrConvg >= 3 )
00099 {
00100 fprintf( ioQQQ, " ConvEdenIoniz: bracket failure 1 n1: %e %e n2: %e %e\n",
00101 n1, error1, n2, error2 );
00102 NeTrack.print_history();
00103 }
00104
00105
00106
00107 while( error1*error2 > 0. && j++ < 20*DEF_ITER )
00108 {
00109 n1 = n2;
00110 error1 = error2;
00111 n2 = ( error1 > 0. ) ? n1*(1.-factor) : n1*(1.+factor);
00112 error2 = EdenError( n2 );
00113 NeTrack.add( n2, error2 );
00114 }
00115
00116 if( error1*error2 > 0. && trace.nTrConvg >= 3 )
00117 {
00118 fprintf( ioQQQ, " ConvEdenIoniz: bracket failure 2 n1: %e %e n2: %e %e\n",
00119 n1, error1, n2, error2 );
00120 NeTrack.print_history();
00121 }
00122
00123 NeTrack.clear();
00124
00125
00126 if( NeTrack.init_bracket( n1, error1, n2, error2 ) == 0 )
00127 {
00128
00129
00130 NeTrack.set_tol(2.*DBL_EPSILON*n2);
00131
00132 double NeNew = 0.5*(n1+n2);
00133 for( int i = 0; i < (1<<(n/2))*DEF_ITER; i++ )
00134 {
00135
00136 if( lgConvEden() || NeTrack.bracket_width() < 3.*DBL_EPSILON*n2 )
00137 break;
00138
00139 NeTrack.add( NeNew, EdenError( NeNew ) );
00140 NeNew = NeTrack.next_val(factor);
00141 }
00142 }
00143
00144 if( conv.lgConvEden )
00145 break;
00146
00147 if( trace.nTrConvg >= 3 )
00148 {
00149 fprintf( ioQQQ, " ConvEdenIoniz: brent fails\n" );
00150 NeTrack.print_history();
00151 }
00152 }
00153
00154 if( trace.lgTrace || trace.nTrConvg >= 3 )
00155 {
00156 fprintf( ioQQQ, " ConvEdenIoniz: entry eden %.4e -> %.4e rel chng %.2f%% accuracy %.2f%%\n",
00157 EdenEntry, dense.eden, (safe_div(dense.eden,EdenEntry,1.)-1.)*100.,
00158 (safe_div(dense.eden,dense.EdenTrue,1.)-1.)*100. );
00159 fprintf( ioQQQ, " ConvEdenIoniz returns converged=%c reason %s\n",
00160 TorF(conv.lgConvEden), conv.chConvIoniz() );
00161 }
00162 }
00163 else
00164 {
00165 fprintf( ioQQQ, "ConvEdenIoniz finds insane solver %s\n", conv.chSolverEden );
00166 ShowMe();
00167 }
00168
00169 return 0;
00170 }
00171
00172
00173 STATIC bool lgConvEden(void)
00174 {
00175 conv.lgConvEden = ( abs(dense.eden-dense.EdenTrue) < abs(dense.eden)*conv.EdenErrorAllowed );
00176 if( !conv.lgConvEden )
00177 {
00178 conv.setConvIonizFail( "Ne big chg" , dense.EdenTrue, dense.eden);
00179 }
00180 return conv.lgConvEden;
00181 }
00182
00183
00184 STATIC double EdenError(double eden)
00185 {
00186
00187 ASSERT( eden > 0. );
00188
00189
00190 conv.incrementCounter(EDEN_CHANGES);
00191 EdenChange( eden );
00192
00193 for( int i=0; i < 5; ++i )
00194 {
00195 if( ConvIoniz() )
00196 lgAbort = true;
00197
00198 if( conv.lgConvIoniz() )
00199 break;
00200 }
00201
00202 double error = dense.eden - dense.EdenTrue;
00203
00204 if( trace.nTrConvg >= 3 )
00205 fprintf( ioQQQ, " EdenError: eden %.4e EdenTrue %.4e rel. err. %.4e\n",
00206 dense.eden, dense.EdenTrue, safe_div(dense.eden,dense.EdenTrue,1.)-1. );
00207
00208 return error;
00209 }
00210