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