00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "dense.h"
00007 #include "trace.h"
00008 #include "thermal.h"
00009 #include "thirdparty.h"
00010 #include "phycon.h"
00011 #include "conv.h"
00012
00013
00014 static const int LOOPMAX = 35;
00015
00016
00017
00018 int ConvEdenIoniz(void)
00019 {
00020 long int LoopEden ,
00021
00022 LoopLimit ,
00023 LoopBig;
00024 bool lgFailLastTry;
00025 int kase = -1;
00026 static bool lgSlope_oscil=false;
00027
00028 static bool lgEden_ever_oscil = false;
00029 static double eden_assumed_old ,
00030 eden_assumed_minus_true_old,
00031 eden_assumed_minus_true,
00032 slope_ls=0., slope_ls_lastzone=0.;
00033 static bool lgEvalSlop_ls=0;
00034 double EdenEntry;
00035 static bool lgMustMalloc_temp_history = true;
00036 static double slope=-1.;
00037 double
00038
00039 EdenUpperLimit,
00040 EdenLowerLimit ,
00041 change_eden_rel2_tolerance ,
00042 PreviousChange;
00043
00044 DEBUG_ENTRY( "ConvEdenIoniz()" );
00045
00046
00047
00048
00049
00050 EdenEntry = dense.eden;
00051
00052 if( trace.lgTrace )
00053 {
00054 fprintf( ioQQQ, " \n" );
00055 fprintf( ioQQQ, " ConvEdenIoniz entered \n" );
00056 }
00057
00058 if( trace.nTrConvg>=3 )
00059 {
00060 fprintf( ioQQQ,
00061 " ConvEdenIoniz3 called, entering eden loop using solver %sw.\n",
00062 conv.chSolverEden);
00063 }
00064
00065
00066
00067 if( !conv.nTotalIoniz )
00068 {
00069
00070 conv.hist_temp_ntemp = -1;
00071 conv.hist_temp_nzone = -1;
00072
00073
00074
00075 }
00076
00077
00078 if( nzone!=conv.hist_temp_nzone )
00079 {
00080
00081 conv.hist_temp_nzone = nzone;
00082
00083
00084 conv.hist_temp_ntemp = 0;
00085 }
00086
00087
00088
00089
00090
00091
00092
00093
00094 if( lgMustMalloc_temp_history )
00095 {
00096 lgMustMalloc_temp_history = false;
00097 conv.hist_temp_limit = 200;
00098 conv.hist_temp_temp = (double *)MALLOC( (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double) ) );
00099 conv.hist_temp_heat = (double *)MALLOC( (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double) ) );
00100 conv.hist_temp_cool = (double *)MALLOC( (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double) ) );
00101 conv.chNotConverged = (char *)MALLOC( (size_t)((unsigned long)INPUT_LINE_LENGTH*sizeof(char) ) );
00102 conv.chConvEden = (char *)MALLOC( (size_t)((unsigned long)INPUT_LINE_LENGTH*sizeof(char) ) );
00103 conv.chConvIoniz = (char *)MALLOC( (size_t)((unsigned long)INPUT_LINE_LENGTH*sizeof(char) ) );
00104 strcpy( conv.chNotConverged, "none" );
00105 strcpy( conv.chConvEden, "none" );
00106 strcpy( conv.chConvIoniz, "none" );
00107 }
00108
00109 if( (conv.hist_temp_ntemp+1) >= conv.hist_temp_limit )
00110 {
00111 conv.hist_temp_limit *= 3;
00112 conv.hist_temp_temp = (double *)REALLOC( conv.hist_temp_temp, (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double)));
00113 conv.hist_temp_heat = (double *)REALLOC( conv.hist_temp_heat, (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double)));
00114 conv.hist_temp_cool = (double *)REALLOC( conv.hist_temp_cool, (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double)));
00115 }
00116
00117
00118 conv.hist_temp_temp[conv.hist_temp_ntemp] = phycon.te;
00119 conv.hist_temp_heat[conv.hist_temp_ntemp] = thermal.htot;
00120 conv.hist_temp_cool[conv.hist_temp_ntemp] = thermal.ctot;
00121 ++conv.hist_temp_ntemp;
00122
00123 if( conv.nPres2Ioniz < 5 )
00124 lgEden_ever_oscil = false;
00125
00126
00127
00128 # define PRT_FAIL_LAST_TRY false
00129 lgFailLastTry = false;
00130 if( PRT_FAIL_LAST_TRY ) fprintf(ioQQQ,"lgFailLastTry set false\n");
00131
00132
00133
00134
00135 LoopBig = 0;
00136 while( LoopBig==0 || (lgFailLastTry && LoopBig==1 ) )
00137 {
00138
00139 if( strcmp( conv.chSolverEden , "simple" )== 0 )
00140 {
00141 static double damp;
00142 long int nLoopOscil;
00143 LoopEden = 0;
00144 conv.nConvIonizFails = 0;
00145
00146
00147
00148 LoopLimit = LOOPMAX;
00149
00150
00151 conv.lgEdenOscl = false;
00152 nLoopOscil = 0;
00153 # define PRT_EDEN_OSCIL false
00154 if( PRT_EDEN_OSCIL ) fprintf(ioQQQ," PRT_EDEN_OSCIIL setting lgEdenOscl FALSE2\n");
00155 lgFailLastTry = false;
00156 if( PRT_FAIL_LAST_TRY ) fprintf(ioQQQ,"lgFailLastTry set false\n");
00157
00158
00159 PreviousChange = 0.;
00160
00161 strcpy( conv.chConvIoniz, " NONE!!!!!" );
00162
00163
00164 eden_assumed_old = 0.;
00165 eden_assumed_minus_true_old = 0.;
00166
00167 damp = 1.;
00168
00169 if( !nzone )
00170 slope = 0.;
00171
00172
00173 do
00174 {
00175 double slopenew , slopeold;
00176
00177
00178 conv.lgConvEden = true;
00179
00180
00181 if( ConvIoniz() || lgAbort )
00182 {
00183 DEBUG_EXIT( "ConvEdenIoniz()" );
00184 return 1;
00185 }
00186
00187
00188 eden_assumed_minus_true = dense.eden - dense.EdenTrue;
00189
00190 {
00191 static long int limEdenHistory=1000;
00192 static bool lgMustMalloc_eden_error_history = true;
00193 bool lgHistUpdate=false;
00194 static double *eden_error_history=NULL, *eden_assumed_history=NULL ,
00195 *stdarray;
00196 static long int nEval=-1 , nzoneUsed=-1, iterUsed=-1;
00197 if( lgMustMalloc_eden_error_history )
00198 {
00199 lgMustMalloc_eden_error_history = false;
00200 lgSlope_oscil = false;
00201 eden_error_history =
00202 (double*)MALLOC(sizeof(double)*(unsigned)limEdenHistory );
00203 eden_assumed_history =
00204 (double*)MALLOC(sizeof(double)*(unsigned)limEdenHistory );
00205 stdarray =
00206 (double*)MALLOC(sizeof(double)*(unsigned)limEdenHistory );
00207 }
00208 if( nzoneUsed!=nzone || iterUsed!=iteration )
00209 {
00210
00211 if( nzone==0 )
00212 {
00213 lgEvalSlop_ls = true;
00214 slope_ls = 1.;
00215 }
00216
00217 nEval = 0;
00218 nzoneUsed = nzone;
00219 iterUsed = iteration;
00220 lgSlope_oscil = false;
00221 slope_ls_lastzone = slope_ls;
00222 }
00223
00224 else if( !conv.lgSearch && !lgSlope_oscil )
00225 {
00226
00227 int ip = MAX2(0,nEval-1);
00228 if(
00229
00230 fabs(dense.eden - dense.EdenTrue)/dense.eden > conv.EdenErrorAllowed &&
00231
00232 (!nEval ||
00233
00234
00235 fabs((dense.eden - dense.EdenTrue)-eden_error_history[ip]*dense.gas_phase[ipHYDROGEN] )/dense.EdenTrue > 1e-5 ||
00236 fabs(dense.eden-eden_assumed_history[ip]*dense.gas_phase[ipHYDROGEN])/ dense.EdenTrue > 1e-5 ))
00237 {
00238
00239 eden_error_history[nEval] = (dense.eden - dense.EdenTrue)/dense.gas_phase[ipHYDROGEN];
00240 eden_assumed_history[nEval] = dense.eden/dense.gas_phase[ipHYDROGEN];
00241 stdarray[nEval] = 0.;
00242
00243
00244 if( eden_error_history[nEval] != 0. )
00245 {
00246 ++nEval;
00247 lgHistUpdate = true;
00248 }
00249
00250
00251
00252
00253
00254 }
00255 if( nEval>=limEdenHistory )
00256 {
00257
00258 limEdenHistory *=10;
00259 eden_error_history = (double *)REALLOC( eden_error_history, (unsigned)(limEdenHistory)*sizeof(double));
00260 eden_assumed_history = (double *)REALLOC( eden_error_history, (unsigned)(limEdenHistory)*sizeof(double));
00261 stdarray = (double *)REALLOC( eden_error_history, (unsigned)(limEdenHistory)*sizeof(double));
00262 }
00263 if( nEval>3 && lgHistUpdate )
00264 {
00265 double y_intercept, stdx, stdy, slope_save;
00266 slope_save = slope_ls;
00267
00268 if( linfit(nEval,
00269 eden_assumed_history,
00270 eden_error_history,
00271 y_intercept,
00272 stdy,
00273 slope_ls,
00274 stdx) )
00275 {
00276 int i;
00277 fprintf(ioQQQ," ConvEdenIoniz: linfit returns impossible condition, history follows:\n");
00278 fprintf(ioQQQ,"eden_error_history\tstdarray\n");
00279 for( i=0; i<nEval; ++i )
00280 {
00281 fprintf(ioQQQ,"%.3e\t%.4e\n",
00282 eden_error_history[i] ,
00283 stdarray[i] );
00284 }
00285 fprintf(ioQQQ,"\n");
00286 ShowMe();
00287 puts( "[Stop in ConvEdenIoniz]" );
00288 cdEXIT(EXIT_FAILURE);
00289 }
00290 lgEvalSlop_ls = true;
00291
00292 if( slope_ls*slope_save<0. )
00293 {
00294 slope_ls = slope_ls_lastzone;
00295 lgSlope_oscil = true;
00296 }
00297
00298
00299
00300
00301
00302
00303 }
00304 }
00305 }
00306 slopeold = slope;
00307
00308 slopenew = 0.;
00309 if( fabs(eden_assumed_minus_true_old) > 0. )
00310 {
00311 if( fabs( eden_assumed_old-dense.eden ) > SMALLFLOAT )
00312 {
00313 # define OLDFAC 0.75
00314 slopenew = (eden_assumed_minus_true_old - eden_assumed_minus_true) / (eden_assumed_old-dense.eden );
00315
00316 # if 0
00317 if( slope !=0. && slope*slopenew>=0.)
00318 slope = OLDFAC*slope + (1.-OLDFAC)*slopenew;
00319 else if( slope==0.)
00320 slope = slopenew;
00321 # endif
00322 if( slope !=0. )
00323 slope = OLDFAC*slope + (1.-OLDFAC)*slopenew;
00324 else
00325 slope = slopenew;
00326 # undef OLDFAC
00327 }
00328
00329
00330 }
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341 if( fabs(dense.eden-dense.EdenTrue)/SDIV(dense.eden) >=
00342 conv.EdenErrorAllowed )
00343 {
00344 double eden_proposed;
00345
00346
00347 conv.lgConvEden = false;
00348 strcpy( conv.chConvIoniz, "Ne big chg" );
00349
00350
00351 ASSERT( dense.eden > SMALLFLOAT );
00352
00353
00354
00355
00356 if( conv.lgSearch && dense.EdenTrue > SMALLFLOAT )
00357 {
00358 dense.eden = pow(10.,
00359 (log10(dense.eden) + log10(dense.EdenTrue))/ 2.);
00360 }
00361 else
00362 {
00363
00364
00365
00366
00367
00368
00369
00370 if( 0 && dense.eden_from_metals > 0.5 )
00371 {
00372 static long int nzone_eval=-1;
00373 static bool lgOscilkase10 = false;
00374
00375
00376
00377 # define KASE_EDEN_NOT_ION 10
00378 if( nzone != nzone_eval )
00379 {
00380 nzone_eval = nzone;
00381 lgOscilkase10 = false;
00382 }
00383
00384
00385 if( PreviousChange*(dense.EdenTrue-dense.eden) < 0. && nzone_eval > 6 )
00386 lgOscilkase10 = true;
00387
00388
00389
00390 change_eden_rel2_tolerance = 0.2;
00391 if( lgOscilkase10 )
00392 {
00393
00394 change_eden_rel2_tolerance = 0.05;
00395 }
00396 kase = KASE_EDEN_NOT_ION;
00397 }
00398
00399
00400
00401
00402
00403 else if( PreviousChange*(dense.EdenTrue-dense.eden) < 0. )
00404 {
00405
00406 conv.lgEdenOscl = true;
00407 nLoopOscil = LoopEden;
00408
00409
00410 damp = MAX2( 0.05, damp * 0.75 );
00411 if( PRT_EDEN_OSCIL )
00412 fprintf(ioQQQ," PRT_EDEN_OSCIIL setting lgEdenOscl true, previous change %e new change %e\n",
00413 PreviousChange,dense.EdenTrue-dense.eden);
00414
00415
00416
00417
00418 change_eden_rel2_tolerance = 0.5;
00419
00420 change_eden_rel2_tolerance = 0.25;
00421 kase = 0;
00422 }
00423 else if( lgFailLastTry )
00424 {
00425
00426
00427 change_eden_rel2_tolerance = 0.5;
00428 kase = 1;
00429 }
00430
00431
00432 else if( LoopEden > 4 && !lgEden_ever_oscil &&
00433 fabs( (dense.EdenTrue-dense.eden)/dense.eden) > 3.*conv.EdenErrorAllowed )
00434 {
00435
00436
00437
00438 change_eden_rel2_tolerance = 3.;
00439 kase = 2;
00440 }
00441 else if( conv.lgEdenOscl )
00442 {
00443
00444
00445 change_eden_rel2_tolerance = 0.25;
00446 kase = 4;
00447 }
00448 else
00449 {
00450
00451
00452
00453
00454 change_eden_rel2_tolerance = 1.;
00455 kase = 3;
00456 }
00457
00458
00459 PreviousChange = dense.EdenTrue - dense.eden;
00460
00461
00462 eden_assumed_minus_true_old = eden_assumed_minus_true;
00463
00464
00465 eden_assumed_old = dense.eden;
00466
00467
00468
00469 EdenUpperLimit = dense.eden * (1. + damp*conv.EdenErrorAllowed*change_eden_rel2_tolerance);
00470 EdenLowerLimit = dense.eden / (1. + damp*conv.EdenErrorAllowed*change_eden_rel2_tolerance);
00471
00472 # define USE_NEW true
00473 # define PRT_NEW false
00474 if( PRT_NEW && USE_NEW )
00475 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",
00476 nzone,
00477 conv.nPres2Ioniz,
00478 slope ,
00479 slopeold ,
00480 slopenew ,
00481 dense.eden,
00482 dense.EdenTrue,
00483 (dense.eden-dense.EdenTrue)/SDIV(dense.eden) ,
00484 change_eden_rel2_tolerance ,
00485 EdenLowerLimit);
00486 if( lgEvalSlop_ls )
00487 slope = slope_ls;
00488 if( USE_NEW && slope !=0. )
00489 {
00490
00491
00492 eden_proposed = dense.eden + (dense.EdenTrue - dense.eden)/slope_ls * damp;
00493 }
00494 else
00495 {
00496 eden_proposed = dense.EdenTrue;
00497 }
00498 # if 0
00499 {
00500 #include "grainvar.h"
00501 fprintf(ioQQQ,"DEBUG eden grn\t%e\t%e\t%e\t%e\n",
00502 dense.eden, eden_proposed,dense.EdenTrue, -gv.TotalEden);
00503 }
00504 # endif
00505
00506
00507
00508
00509 if( eden_proposed > EdenUpperLimit )
00510 {
00511
00512 dense.eden = EdenUpperLimit;
00513 }
00514
00515 else if( eden_proposed < EdenLowerLimit )
00516 {
00517
00518 dense.eden = EdenLowerLimit;
00519 }
00520 else
00521 {
00522
00523
00524 dense.eden = eden_proposed;
00525 }
00526
00527 if( trace.lgTrace && trace.lgNeBug )
00528 {
00529 fprintf( ioQQQ,
00530 " ConvEdenIoniz, chg ne to %.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e\n",
00531 dense.eden, eden_assumed_old, dense.eden/eden_assumed_old, dense.EdenTrue );
00532 }
00533 }
00534
00535 }
00536
00537
00538
00539 else
00540 {
00541
00542
00543 if( fabs(dense.eden-dense.EdenTrue)/dense.EdenTrue > conv.EdenErrorAllowed/10. &&
00544
00545
00546 kase !=KASE_EDEN_NOT_ION )
00547 {
00548
00549
00550
00551
00552
00553 dense.eden = dense.eden + (dense.EdenTrue-dense.eden)/MAX2(1.,slope_ls)*damp;
00554
00555 if( trace.nTrConvg>=3 )
00556 {
00557 fprintf( ioQQQ,
00558 " ConvEdenIoniz3 converged eden, re-calling ConvIoniz with EdenTrue=%.4e (was %.4e).\n",
00559 dense.EdenTrue ,
00560 dense.eden );
00561 }
00562
00563
00564
00565
00566
00567
00568 if( ConvIoniz() )
00569 {
00570 DEBUG_EXIT( "ConvEdenIoniz()" );
00571 return 1;
00572 }
00573 }
00574 else if( trace.nTrConvg>=3 )
00575 {
00576 fprintf( ioQQQ,
00577 " ConvEdenIoniz3 no need to re-call ConvIoniz since eden did not change much.\n");
00578 }
00579
00580
00581
00582
00583
00584 if(
00585 fabs(dense.eden-dense.EdenTrue)/SDIV(dense.eden) >=
00586 conv.EdenErrorAllowed )
00587 {
00588
00589
00590 conv.lgConvEden = false;
00591
00592
00593
00594 lgFailLastTry = true;
00595 if( PRT_FAIL_LAST_TRY ) fprintf(ioQQQ,"lgFailLastTry set true\n");
00596
00597 lgEden_ever_oscil = true;
00598 }
00599 }
00600 {
00601
00602
00603 enum {DEBUG_LOC=false};
00604
00605 if( DEBUG_LOC )
00606 {
00607 fprintf(ioQQQ,"edendebugg %li\t%.2e\t%.2e\t%.2e\t%.2e\n",
00608 nzone,dense.eden ,eden_assumed_old, (dense.eden-eden_assumed_old)/dense.eden, dense.EdenTrue);
00609 }
00610 }
00611 if( trace.lgTrace && trace.lgNeBug )
00612 {
00613 fprintf( ioQQQ,
00614 " ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e converge=%c\n",
00615 dense.eden, eden_assumed_old, dense.eden/SDIV(eden_assumed_old), dense.EdenTrue ,TorF( conv.lgConvEden));
00616 }
00617
00618 if( trace.nTrConvg>=3 )
00619 {
00620
00621 fprintf( ioQQQ, " ConvEdenIoniz3 loop:%4ld ",
00622 LoopEden );
00623
00624
00625 if( conv.lgConvEden )
00626 {
00627 fprintf( ioQQQ, " converged, eden rel error:%g ",
00628 (dense.EdenTrue-dense.eden)/dense.EdenTrue );
00629 }
00630 else if( eden_assumed_old > SMALLFLOAT )
00631 {
00632
00633 fprintf( ioQQQ, " NOT Converged! corr:%8.4f prop:%9.5f ",
00634 (dense.EdenTrue-eden_assumed_old)/eden_assumed_old ,
00635 (dense.eden-eden_assumed_old)/eden_assumed_old );
00636 }
00637
00638 fprintf( ioQQQ, " new ne:%.6e true:%.6e kase:%i slope:%.3e osc:%c",
00639 dense.eden ,
00640 dense.EdenTrue ,
00641 kase ,
00642 slope_ls,
00643 TorF(lgSlope_oscil));
00644
00645 if( conv.lgEdenOscl )
00646 {
00647 fprintf( ioQQQ, " EDEN OSCIL1 damp %.4f\n", damp);
00648 }
00649 else
00650 {
00651 fprintf( ioQQQ, "\n");
00652 }
00653 }
00654
00655
00656 ++LoopEden;
00657
00658
00659
00660
00661 if( LoopEden - nLoopOscil >4 && fabs(slope_ls) < 3. )
00662 {
00663 conv.lgEdenOscl = false;
00664
00665 damp = 1.;
00666 }
00667
00668
00669
00670
00671 if( LoopEden == (LOOPMAX-1) && !conv.lgEdenOscl && conv.nConvIonizFails==0)
00672
00673
00674 LoopLimit *= 4;
00675
00676 } while( !conv.lgConvEden && LoopEden < LoopLimit && !lgAbort );
00677 }
00678
00679 else if( strcmp( conv.chSolverEden , "new" )== 0 )
00680 {
00681
00682 double PreviousEdenError = 0.;
00683 double CurrentEdenError = 0.;
00684 double CurrentEden = 0.;
00685 double ProposedEden,
00686 FractionChangeEden = 0.;
00687
00688
00689 LoopEden = 0;
00690 conv.nConvIonizFails = 0;
00691
00692
00693 LoopLimit = LOOPMAX;
00694
00695
00696 conv.lgEdenOscl = false;
00697 if( PRT_EDEN_OSCIL ) fprintf(ioQQQ," PRT_EDEN_OSCIIL setting lgEdenOscl FALSE1\n");
00698
00699
00700
00701 change_eden_rel2_tolerance = 0.02;
00702
00703 strcpy( conv.chConvIoniz, " NONE!!!!!" );
00704
00705
00706
00707 do
00708 {
00709
00710
00711
00712 conv.lgConvEden = true;
00713
00714 if( trace.nTrConvg>=3 )
00715 fprintf( ioQQQ, " ConvEdenIoniz3 calling ConvIoniz with eden= %.4e\n",dense.eden);
00716
00717
00718 if( ConvIoniz() )
00719 {
00720 DEBUG_EXIT( "ConvEdenIoniz()" );
00721 return 1;
00722 }
00723
00724
00725
00726
00727
00728
00729
00730
00731 CurrentEden = dense.eden;
00732
00733
00734 CurrentEdenError = dense.EdenTrue - dense.eden;
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746 if(
00747 (LoopEden==0 && fabs(CurrentEdenError)/dense.EdenTrue <= conv.EdenErrorAllowed/2.) ||
00748 (LoopEden>0 && fabs(CurrentEdenError)/dense.EdenTrue <= conv.EdenErrorAllowed
00749 && FractionChangeEden < conv.EdenErrorAllowed / 2.) )
00750 break;
00751
00752
00753
00754 conv.lgConvEden = false;
00755 strcpy( conv.chConvIoniz, "Ne big chg" );
00756
00757
00758 if( conv.lgSearch )
00759 {
00760 dense.eden = pow(10.,
00761 (log10(dense.eden) + log10(dense.EdenTrue))/ 2.);
00762 }
00763 else
00764 {
00765
00766
00767
00768 if( PreviousEdenError*CurrentEdenError < 0. )
00769 {
00770
00771
00772 slope = (PreviousEdenError - CurrentEdenError ) /
00773 ( eden_assumed_old - CurrentEden );
00774
00775 ProposedEden = eden_assumed_old - PreviousEdenError / slope;
00776
00777
00778 ASSERT( ProposedEden >= MIN2( eden_assumed_old , CurrentEden ) );
00779 ASSERT( ProposedEden <= MAX2( eden_assumed_old , CurrentEden ) );
00780
00781 change_eden_rel2_tolerance *= 0.25;
00782 if( trace.nTrConvg>=3 )
00783 fprintf( ioQQQ,
00784 " ConvEdenIoniz3 bracketed electron density factor now %.3e error is %.4f proposed ne %.4e\n",
00785 change_eden_rel2_tolerance,
00786 (dense.EdenTrue-dense.eden)/dense.EdenTrue, ProposedEden);
00787 }
00788 else
00789 {
00790
00791
00792
00793 EdenUpperLimit = dense.eden * (1. + change_eden_rel2_tolerance);
00794 EdenLowerLimit = dense.eden / (1. + change_eden_rel2_tolerance);
00795
00796
00797 if( dense.EdenTrue > EdenUpperLimit )
00798 {
00799
00800 ProposedEden = EdenUpperLimit;
00801 }
00802 else if( dense.EdenTrue < EdenLowerLimit )
00803 {
00804
00805 ProposedEden = EdenLowerLimit;
00806 }
00807 else
00808 {
00809
00810 ProposedEden = dense.EdenTrue;
00811 }
00812 if( trace.nTrConvg>=3 )
00813 fprintf( ioQQQ,
00814 " ConvEdenIoniz3 elec dens fac %.3e err is %.4f prop ne %.4e cor ne %.4e slope=%.2e\n",
00815 change_eden_rel2_tolerance,
00816 (dense.EdenTrue-dense.eden)/dense.EdenTrue ,
00817 ProposedEden,
00818 dense.EdenTrue,slope);
00819 }
00820
00821
00822 PreviousEdenError = CurrentEdenError;
00823 eden_assumed_old = CurrentEden;
00824 FractionChangeEden = fabs( dense.eden - ProposedEden ) / dense.EdenTrue;
00825 dense.eden = ProposedEden;
00826
00827 if( trace.lgTrace && trace.lgNeBug )
00828 {
00829 fprintf( ioQQQ,
00830 " ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e\n",
00831 dense.eden, eden_assumed_old, dense.eden/eden_assumed_old, dense.EdenTrue );
00832 }
00833 }
00834
00835 ++LoopEden;
00836 } while( LoopEden < LoopLimit && !lgAbort );
00837
00838
00839
00840
00841
00842 if( trace.nTrConvg>=3 )
00843 fprintf( ioQQQ,
00844 " ConvEdenIoniz3 converged eden, current error is %.4f, calling ConvIoniz with final density=EdenTrue=%.4e\n",
00845 (dense.EdenTrue - dense.eden)/dense.EdenTrue,dense.EdenTrue);
00846 dense.eden = dense.EdenTrue;
00847
00848
00849 if( ConvIoniz() )
00850 {
00851 if( trace.nTrConvg>=3 )
00852 fprintf( ioQQQ,
00853 " ConvEdenIoniz3 eden no longer converged eden, current error is %.4f\n",
00854 (dense.EdenTrue - dense.eden)/dense.EdenTrue);
00855 }
00856
00857
00858
00859 if(
00860 fabs(dense.eden-dense.EdenTrue)/dense.EdenTrue >
00861 conv.EdenErrorAllowed )
00862 {
00863
00864
00865 conv.lgConvEden = false;
00866 if( trace.nTrConvg>=3 )
00867 fprintf( ioQQQ,
00868 " ConvEdenIoniz3 setting eden not converged, error %.4f, exiting\n",
00869 (dense.eden-dense.EdenTrue)/dense.EdenTrue );
00870 }
00871 else
00872 {
00873 conv.lgConvEden = true;
00874 }
00875
00876 if( trace.lgTrace && trace.lgNeBug )
00877 {
00878 fprintf( ioQQQ,
00879 " ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e converge=%c\n",
00880 dense.eden, eden_assumed_old, dense.eden/eden_assumed_old, dense.EdenTrue ,TorF( conv.lgConvEden));
00881
00882
00883 }
00884
00885 if( trace.nTrConvg>=3 )
00886 {
00887 fprintf( ioQQQ, " ConvEdenIoniz3 %4ld new ne=%12.4e true=%12.4e",
00888 LoopEden, dense.eden , dense.EdenTrue );
00889
00890
00891 if( conv.lgConvEden )
00892 {
00893 fprintf( ioQQQ, " converged, eden rel error is %g ",
00894 (dense.EdenTrue-dense.eden)/dense.EdenTrue );
00895 }
00896 else
00897 {
00898 fprintf( ioQQQ, " NOCONV corr:%7.3f prop:%7.3f ",
00899 (dense.EdenTrue-eden_assumed_old)/eden_assumed_old ,
00900 (dense.eden-eden_assumed_old)/eden_assumed_old );
00901 }
00902 if( conv.lgEdenOscl )
00903 fprintf( ioQQQ, " EDEN OSCIL2 \n");
00904 }
00905
00906 }
00907 else
00908 {
00909 fprintf( ioQQQ, "ConvEdenIoniz finds insane solver %s \n" , conv.chSolverEden );
00910 ShowMe();
00911 }
00912 ++LoopBig;
00913 }
00914
00915 if( trace.nTrConvg>=3 )
00916 {
00917 fprintf( ioQQQ, " ConvEdenIoniz3 exits, lgConvEden=%1c entry eden %.4e -> %.4e rel change %.3f\n" ,
00918 TorF(conv.lgConvEden ),
00919 EdenEntry ,
00920 dense.eden ,
00921 (EdenEntry-dense.eden)/SDIV( dense.eden ) );
00922 }
00923 else if( trace.lgTrace )
00924 {
00925 fprintf( ioQQQ, " ConvEdenIoniz exits zn %.2f lgConvEden=%1c entry eden %.4e -> %.4e rel change %.3f\n" ,
00926 fnzone,
00927 TorF(conv.lgConvEden ),
00928 EdenEntry ,
00929 dense.eden ,
00930 (EdenEntry-dense.eden)/SDIV( dense.eden ) );
00931 }
00932
00933 DEBUG_EXIT( "ConvEdenIoniz()" );
00934
00935 return 0;
00936 }
00937
00938
00939 bool lgConvEden(void)
00940 {
00941 bool lgRet;
00942
00943
00944 if( fabs(dense.eden-dense.EdenTrue)/SDIV(dense.eden) >=
00945 conv.EdenErrorAllowed )
00946 {
00947 lgRet = false;
00948 conv.lgConvEden = false;
00949 }
00950 else
00951 {
00952 lgRet = true;
00953 conv.lgConvEden = true;
00954 }
00955 return lgRet;
00956 }