26 static bool lgSlope_oscil=
false;
28 static bool lgEden_ever_oscil =
false;
29 static double eden_assumed_old ,
30 eden_assumed_minus_true_old,
31 eden_assumed_minus_true,
32 slope_ls=0., slope_ls_lastzone=0.;
33 static bool lgEvalSlop_ls=0;
35 static bool lgMustMalloc_temp_history =
true;
36 static double slope=-1.;
41 change_eden_rel2_tolerance ,
54 fprintf(
ioQQQ,
" \n" );
55 fprintf(
ioQQQ,
" ConvEdenIoniz entered \n" );
61 " ConvEdenIoniz3 called, entering eden loop using solver %sw.\n",
94 if( lgMustMalloc_temp_history )
96 lgMustMalloc_temp_history =
false;
124 lgEden_ever_oscil =
false;
128 # define PRT_FAIL_LAST_TRY false
129 lgFailLastTry =
false;
136 while( LoopBig==0 || (lgFailLastTry && LoopBig==1 ) )
153 # define PRT_EDEN_OSCIL false
155 lgFailLastTry =
false;
164 eden_assumed_old = 0.;
165 eden_assumed_minus_true_old = 0.;
175 double slopenew , slopeold;
190 static long int limEdenHistory=1000;
191 static bool lgMustMalloc_eden_error_history =
true;
192 bool lgHistUpdate=
false;
193 static double *eden_error_history=NULL, *eden_assumed_history=NULL ,
195 static long int nEval=-1 , nzoneUsed=-1, iterUsed=-1;
196 if( lgMustMalloc_eden_error_history )
198 lgMustMalloc_eden_error_history =
false;
199 lgSlope_oscil =
false;
201 (
double*)
MALLOC(
sizeof(
double)*(unsigned)limEdenHistory );
202 eden_assumed_history =
203 (
double*)
MALLOC(
sizeof(
double)*(unsigned)limEdenHistory );
205 (
double*)
MALLOC(
sizeof(
double)*(unsigned)limEdenHistory );
212 lgEvalSlop_ls =
true;
219 lgSlope_oscil =
false;
220 slope_ls_lastzone = slope_ls;
226 int ip =
MAX2(0,nEval-1);
240 stdarray[nEval] = 0.;
243 if( eden_error_history[nEval] != 0. )
254 if( nEval>=limEdenHistory )
258 eden_error_history = (
double *)
REALLOC( eden_error_history, (
unsigned)(limEdenHistory)*
sizeof(
double));
259 eden_assumed_history = (
double *)
REALLOC( eden_error_history, (
unsigned)(limEdenHistory)*
sizeof(
double));
260 stdarray = (
double *)
REALLOC( eden_error_history, (
unsigned)(limEdenHistory)*
sizeof(
double));
262 if( nEval>3 && lgHistUpdate )
264 double y_intercept, stdx, stdy, slope_save;
265 slope_save = slope_ls;
268 eden_assumed_history,
276 fprintf(
ioQQQ,
" ConvEdenIoniz: linfit returns impossible condition, history follows:\n");
277 fprintf(
ioQQQ,
"eden_error_history\tstdarray\n");
278 for( i=0; i<nEval; ++i )
280 fprintf(
ioQQQ,
"%.3e\t%.4e\n",
281 eden_error_history[i] ,
288 lgEvalSlop_ls =
true;
290 if( slope_ls*slope_save<0. )
292 slope_ls = slope_ls_lastzone;
293 lgSlope_oscil =
true;
307 if( fabs(eden_assumed_minus_true_old) > 0. )
312 slopenew = (eden_assumed_minus_true_old - eden_assumed_minus_true) / (eden_assumed_old-
dense.
eden );
315 if( slope !=0. && slope*slopenew>=0.)
342 double eden_proposed;
370 static long int nzone_eval=-1;
371 static bool lgOscilkase10 =
false;
375 # define KASE_EDEN_NOT_ION 10
376 if(
nzone != nzone_eval )
379 lgOscilkase10 =
false;
384 lgOscilkase10 =
true;
388 change_eden_rel2_tolerance = 0.2;
392 change_eden_rel2_tolerance = 0.05;
405 nLoopOscil = LoopEden;
408 damp =
MAX2( 0.05, damp * 0.75 );
410 fprintf(
ioQQQ,
" PRT_EDEN_OSCIIL setting lgEdenOscl true, previous change %e new change %e\n",
416 change_eden_rel2_tolerance = 0.5;
418 change_eden_rel2_tolerance = 0.25;
421 else if( lgFailLastTry )
425 change_eden_rel2_tolerance = 0.5;
430 else if( LoopEden > 4 && !lgEden_ever_oscil &&
436 change_eden_rel2_tolerance = 3.;
443 change_eden_rel2_tolerance = 0.25;
452 change_eden_rel2_tolerance = 1.;
460 eden_assumed_minus_true_old = eden_assumed_minus_true;
470 # define USE_NEW true
471 # define PRT_NEW false
473 fprintf(
ioQQQ,
"DEBUG slope\t%li\t%li\t%.3f\t%.3f\t%.3f\t%.4e\t%.4e\t%.4f\t%.3f\t%.3e\n",
482 change_eden_rel2_tolerance ,
499 fprintf(
ioQQQ,
"DEBUG eden grn\t%e\t%e\t%e\t%e\n",
507 if( eden_proposed > EdenUpperLimit )
513 else if( eden_proposed < EdenLowerLimit )
528 " ConvEdenIoniz, chg ne to %.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e\n",
556 " ConvEdenIoniz3 converged eden, re-calling ConvIoniz with EdenTrue=%.4e (was %.4e).\n",
574 " ConvEdenIoniz3 no need to re-call ConvIoniz since eden did not change much.\n");
591 lgFailLastTry =
true;
594 lgEden_ever_oscil =
true;
600 enum {DEBUG_LOC=
false};
604 fprintf(
ioQQQ,
"edendebugg %li\t%.2e\t%.2e\t%.2e\t%.2e\n",
611 " ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e converge=%c\n",
618 fprintf(
ioQQQ,
" ConvEdenIoniz3 loop:%4ld ",
624 fprintf(
ioQQQ,
" converged, eden rel error:%g ",
630 fprintf(
ioQQQ,
" NOT Converged! corr:%8.4f prop:%9.5f ",
632 (
dense.
eden-eden_assumed_old)/eden_assumed_old );
635 fprintf(
ioQQQ,
" new ne:%.6e true:%.6e kase:%i slope:%.3e osc:%c",
640 TorF(lgSlope_oscil));
644 fprintf(
ioQQQ,
" EDEN OSCIL1 damp %.4f\n", damp);
648 fprintf(
ioQQQ,
"\n");
658 if( LoopEden - nLoopOscil >4 && fabs(slope_ls) < 3. )
679 double PreviousEdenError = 0.;
680 double CurrentEdenError = 0.;
681 double CurrentEden = 0.;
683 FractionChangeEden = 0.;
698 change_eden_rel2_tolerance = 0.02;
712 fprintf(
ioQQQ,
" ConvEdenIoniz3 calling ConvIoniz with eden= %.4e\n",
dense.
eden);
764 if( PreviousEdenError*CurrentEdenError < 0. )
768 slope = (PreviousEdenError - CurrentEdenError ) /
769 ( eden_assumed_old - CurrentEden );
771 ProposedEden = eden_assumed_old - PreviousEdenError / slope;
774 ASSERT( ProposedEden >=
MIN2( eden_assumed_old , CurrentEden ) );
775 ASSERT( ProposedEden <=
MAX2( eden_assumed_old , CurrentEden ) );
777 change_eden_rel2_tolerance *= 0.25;
780 " ConvEdenIoniz3 bracketed electron density factor now %.3e error is %.4f proposed ne %.4e\n",
781 change_eden_rel2_tolerance,
789 EdenUpperLimit =
dense.
eden * (1. + change_eden_rel2_tolerance);
790 EdenLowerLimit =
dense.
eden / (1. + change_eden_rel2_tolerance);
796 ProposedEden = EdenUpperLimit;
801 ProposedEden = EdenLowerLimit;
810 " ConvEdenIoniz3 elec dens fac %.3e err is %.4f prop ne %.4e cor ne %.4e slope=%.2e\n",
811 change_eden_rel2_tolerance,
818 PreviousEdenError = CurrentEdenError;
819 eden_assumed_old = CurrentEden;
826 " ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e\n",
832 }
while( LoopEden < LoopLimit && !
lgAbort );
840 " ConvEdenIoniz3 converged eden, current error is %.4f, calling ConvIoniz with final density=EdenTrue=%.4e\n",
849 " ConvEdenIoniz3 eden no longer converged eden, current error is %.4f\n",
864 " ConvEdenIoniz3 setting eden not converged, error %.4f, exiting\n",
875 " ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e converge=%c\n",
881 fprintf(
ioQQQ,
" ConvEdenIoniz3 %4ld new ne=%12.4e true=%12.4e",
887 fprintf(
ioQQQ,
" converged, eden rel error is %g ",
892 fprintf(
ioQQQ,
" NOCONV corr:%7.3f prop:%7.3f ",
894 (
dense.
eden-eden_assumed_old)/eden_assumed_old );
897 fprintf(
ioQQQ,
" EDEN OSCIL2 \n");
911 fprintf(
ioQQQ,
" ConvEdenIoniz3 exits, lgConvEden=%1c entry eden %.4e -> %.4e rel change %.3f\n" ,
919 fprintf(
ioQQQ,
" ConvEdenIoniz exits zn %.2f lgConvEden=%1c entry eden %.4e -> %.4e rel change %.3f\n" ,