cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*ConvEdenIoniz called by ConvTempIonz, calls ConvIoniz solving for eden */ 00004 /*lgConvEden returns true if electron density is converged */ 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 /* limit to how many loops we will do */ 00013 /* >>chng 04 sep 27, from 25 to 35 */ 00014 static const int LOOPMAX = 35; 00015 00016 /*ConvEdenIoniz called by ConvTempIonz, calls ConvIoniz solving for eden 00017 * returns 0 if ok, 1 if abort */ 00018 int ConvEdenIoniz(void) 00019 { 00020 long int LoopEden , 00021 /* this will be LOOPMAX at first, then doubled if no oscillations occur*/ 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 /* upper and lower limits to the range of the change we want to consider */ 00039 EdenUpperLimit, 00040 EdenLowerLimit , 00041 change_eden_rel2_tolerance , 00042 PreviousChange; 00043 00044 DEBUG_ENTRY( "ConvEdenIoniz()" ); 00045 00046 /* this routine is called by ConvTempIonz, it calls ConvIonize 00047 * and changes the electron density until it converges */ 00048 00049 /* save entry value of eden */ 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 /* keep track of temperature solver in this zone - 00066 * conv.nTotalIoniz is zero in first call of new iteration */ 00067 if( !conv.nTotalIoniz ) 00068 { 00069 /* one time initialization of variables */ 00070 conv.hist_temp_ntemp = -1; 00071 conv.hist_temp_nzone = -1; 00072 /*>>chng 06 jun 25, do not reset upper limit to number of 00073 * points saved - only set this when malloc or remalloc */ 00074 /*conv.hist_temp_limit = 0;*/ 00075 } 00076 /* this will save the history of density - pressure relationship 00077 * for the current zone */ 00078 if( nzone!=conv.hist_temp_nzone ) 00079 { 00080 /* first time in this zone - reset counters */ 00081 conv.hist_temp_nzone = nzone; 00082 /* this counter will be updated after vars are saved so will be 00083 * total number of saved values */ 00084 conv.hist_temp_ntemp = 0; 00085 } 00086 /* do we need to allocate, or reallocate, memory for the vectors 00087 * check initial allocation first */ 00088 /*>>chng 06 jun 25, do not test on this limit for need to malloc 00089 * rather set static flag in this routine - fix of memory leak 00090 * discovered by PvH - routine entered with conv.nTotalIoniz zero 00091 * during first call of each iteration so memory was malloced without 00092 * freeing previous memory - memory leak resulted */ 00093 /*if( conv.hist_temp_limit==0 )*/ 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 /* ran out of space - allocate more */ 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 /* >>chng 04 feb 11, add option to remember current density and pressure */ 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 /* this flag says wether we converged the density but failed after final tweek, 00127 * need to use smaller steps */ 00128 # define PRT_FAIL_LAST_TRY false 00129 lgFailLastTry = false; 00130 if( PRT_FAIL_LAST_TRY ) fprintf(ioQQQ,"lgFailLastTry set false\n"); 00131 00132 /* >>chng 03 mar 17, add this big loop, only go through one time if fully 00133 * converged at end, but if eden jumps off in final touchup, 00134 * do whole thing again */ 00135 LoopBig = 0; 00136 while( LoopBig==0 || (lgFailLastTry && LoopBig==1 ) ) 00137 { 00138 /* the old default solver */ 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 /* this will be increased by 2x if, at the end, no oscillations have occurred */ 00147 /* >>chng 04 aug 04, from 2x to 4x, to following through on eden front */ 00148 LoopLimit = LOOPMAX; 00149 00150 /* will be set true if sign of change, ever changes */ 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 /* will be used to look for oscillations in the electron density */ 00159 PreviousChange = 0.; 00160 00161 strcpy( conv.chConvIoniz, " NONE!!!!!" ); 00162 00163 /* we have not yet */ 00164 eden_assumed_old = 0.; 00165 eden_assumed_minus_true_old = 0.; 00166 00167 damp = 1.; 00168 /* if working with zone 0 (nzone=0) reset slope to zero */ 00169 if( !nzone ) 00170 slope = 0.; 00171 00172 /* this is electron density convergence loop */ 00173 do 00174 { 00175 double slopenew , slopeold; 00176 /* this flag will be set false below if electron density not within tolerance 00177 * after ionization is recomputed */ 00178 conv.lgConvEden = true; 00179 00180 /* converge the current level of ionization, this calls eden_sum which updates EdenTrue */ 00181 if( ConvIoniz() || lgAbort ) 00182 { 00183 return 1; 00184 } 00185 /* this is the new error in the electron density, the difference between 00186 * the assumed and true electron density - we want this to be zero */ 00187 eden_assumed_minus_true = dense.eden - dense.EdenTrue; 00188 00189 { 00190 static long int limEdenHistory=1000; 00191 static bool lgMustMalloc_eden_error_history = true; 00192 bool lgHistUpdate=false; 00193 static double *eden_error_history=NULL, *eden_assumed_history=NULL , 00194 *stdarray; 00195 static long int nEval=-1 , nzoneUsed=-1, iterUsed=-1; 00196 if( lgMustMalloc_eden_error_history ) 00197 { 00198 lgMustMalloc_eden_error_history = false; 00199 lgSlope_oscil = false; 00200 eden_error_history = 00201 (double*)MALLOC(sizeof(double)*(unsigned)limEdenHistory ); 00202 eden_assumed_history = 00203 (double*)MALLOC(sizeof(double)*(unsigned)limEdenHistory ); 00204 stdarray = 00205 (double*)MALLOC(sizeof(double)*(unsigned)limEdenHistory ); 00206 } 00207 if( nzoneUsed!=nzone || iterUsed!=iteration ) 00208 { 00209 /* first evaluation this zone */ 00210 if( nzone==0 ) 00211 { 00212 lgEvalSlop_ls = true; 00213 slope_ls = 1.; 00214 } 00215 /* reset some vars at start of new zone */ 00216 nEval = 0; 00217 nzoneUsed = nzone; 00218 iterUsed = iteration; 00219 lgSlope_oscil = false; 00220 slope_ls_lastzone = slope_ls; 00221 } 00222 /* do not evaluate this during search phase or if slope is oscillating */ 00223 else if( !conv.lgSearch && !lgSlope_oscil ) 00224 { 00225 /* may cycle through here with exactly same edensity before tripping if */ 00226 int ip = MAX2(0,nEval-1); 00227 if( 00228 /* no need to do this if already converged */ 00229 fabs(dense.eden - dense.EdenTrue)/dense.eden > conv.EdenErrorAllowed && 00230 /* do it if never evaluated */ 00231 (!nEval || 00232 /* this test - don't want tiny corrections, within the tolerance, which 00233 * can have noise since at the level code is converging, affecting the slope */ 00234 fabs((dense.eden - dense.EdenTrue)-eden_error_history[ip]*dense.gas_phase[ipHYDROGEN] )/dense.EdenTrue > 1e-5 || 00235 fabs(dense.eden-eden_assumed_history[ip]*dense.gas_phase[ipHYDROGEN])/ dense.EdenTrue > 1e-5 )) 00236 { 00237 /* use relative fraction for constant pressure case */ 00238 eden_error_history[nEval] = (dense.eden - dense.EdenTrue)/dense.gas_phase[ipHYDROGEN]; 00239 eden_assumed_history[nEval] = dense.eden/dense.gas_phase[ipHYDROGEN]; 00240 stdarray[nEval] = 0.; 00241 /* >>chng 05 feb 22, only update nEval if residuals are non-zero - 00242 * this happens in constant temperature models - happened in optimization run */ 00243 if( eden_error_history[nEval] != 0. ) 00244 { 00245 ++nEval; 00246 lgHistUpdate = true; 00247 } 00248 /*fprintf(ioQQQ,"DEBUG history\t%li\t%e\t%e\t%e\n", 00249 nEval, 00250 dense.gas_phase[ipHYDROGEN], 00251 eden_assumed_history[nEval], 00252 eden_error_history[nEval]);*/ 00253 } 00254 if( nEval>=limEdenHistory ) 00255 { 00256 /* used a lot of points - must create more space */ 00257 limEdenHistory *=10; 00258 eden_error_history = (double *)REALLOC( eden_error_history, (unsigned)(limEdenHistory)*sizeof(double)); 00259 eden_assumed_history = (double *)REALLOC( eden_error_history, (unsigned)(limEdenHistory)*sizeof(double)); 00260 stdarray = (double *)REALLOC( eden_error_history, (unsigned)(limEdenHistory)*sizeof(double)); 00261 } 00262 if( nEval>3 && lgHistUpdate ) 00263 { 00264 double y_intercept, stdx, stdy, slope_save; 00265 slope_save = slope_ls; 00266 /* enough data to do least squares */ 00267 if( linfit(nEval, 00268 eden_assumed_history, 00269 eden_error_history, 00270 y_intercept, 00271 stdy, 00272 slope_ls, 00273 stdx) ) 00274 { 00275 int i; 00276 fprintf(ioQQQ," ConvEdenIoniz: linfit returns impossible condition, history follows:\n"); 00277 fprintf(ioQQQ,"eden_error_history\tstdarray\n"); 00278 for( i=0; i<nEval; ++i ) 00279 { 00280 fprintf(ioQQQ,"%.3e\t%.4e\n", 00281 eden_error_history[i] , 00282 stdarray[i] ); 00283 } 00284 fprintf(ioQQQ,"\n"); 00285 ShowMe(); 00286 cdEXIT(EXIT_FAILURE); 00287 } 00288 lgEvalSlop_ls = true; 00289 /* check for slope changing sign - if does, use last zone's slope */ 00290 if( slope_ls*slope_save<0. ) 00291 { 00292 slope_ls = slope_ls_lastzone; 00293 lgSlope_oscil = true; 00294 } 00295 /*fprintf(ioQQQ,"DEBUG LS\t%.2f\t%li\t%.3e\t%.3e\t%.3e\n", 00296 fnzone , 00297 nEval , 00298 slope_ls, 00299 y_intercept , 00300 slope);*/ 00301 } 00302 } 00303 } 00304 slopeold = slope; 00305 /* may not be set below, but could print it. set to zero as flag for this case */ 00306 slopenew = 0.; 00307 if( fabs(eden_assumed_minus_true_old) > 0. ) 00308 { 00309 if( fabs( eden_assumed_old-dense.eden ) > SMALLFLOAT ) 00310 { 00311 # define OLDFAC 0.75 00312 slopenew = (eden_assumed_minus_true_old - eden_assumed_minus_true) / (eden_assumed_old-dense.eden ); 00313 /* >>chng 04 sep 15 do not update slope if changing sign */ 00314 # if 0 00315 if( slope !=0. && slope*slopenew>=0.) 00316 slope = OLDFAC*slope + (1.-OLDFAC)*slopenew; 00317 else if( slope==0.) 00318 slope = slopenew; 00319 # endif 00320 if( slope !=0. ) 00321 slope = OLDFAC*slope + (1.-OLDFAC)*slopenew; 00322 else 00323 slope = slopenew; 00324 # undef OLDFAC 00325 } 00326 /*else 00327 slope = 0.;*/ 00328 } 00329 /* following block is to not let electron density change by 00330 * too much 00331 * conv.EdenErrorAllowed is allowed error 00332 * the default value of conv.EdenErrorAllowed is 0.01 in zerologic 00333 * is changed with the SET Eden Error command 00334 * eden_assumed_old was incoming value of eden 00335 * EdenTrue is correct value from eden_sum for current ionization 00336 * new eden will be set using these, trying to prevent oscillations */ 00337 00338 /* is error larger than the tolerance we are trying to beat? */ 00339 if( fabs(dense.eden-dense.EdenTrue)/SDIV(dense.eden) >= 00340 conv.EdenErrorAllowed ) 00341 { 00342 double eden_proposed; 00343 /* difference is assumed and true electron density is larger 00344 * than tolerance, we are not yet converged, default is 0.01 */ 00345 conv.lgConvEden = false; 00346 strcpy( conv.chConvIoniz, "Ne big chg" ); 00347 00348 /* these are needed since will take log below */ 00349 ASSERT( dense.eden > SMALLFLOAT ); 00350 00351 /* SEARCH is true if this is only search for first solution */ 00352 /* EdenTrue is allowed to go negative during search for soln - not physical, 00353 * but can happen in math search for soln */ 00354 if( conv.lgSearch && dense.EdenTrue > SMALLFLOAT ) 00355 { 00356 dense.eden = pow(10., 00357 (log10(dense.eden) + log10(dense.EdenTrue))/ 2.); 00358 } 00359 else 00360 { 00361 /* >>chng 03 jul 07, add case for grains contain significant 00362 * fraction of electrons */ 00363 /* >>chng 03 nov 28, add this test on fraction of electrons from ions, 00364 * this branch is to identify limit where molecules and grains have 00365 * most of the electrons */ 00366 /* >>chng 04 sep 14, change from all ions to just metals */ 00367 /* this was patch to fix steep slope */ 00368 if( 0 && dense.eden_from_metals > 0.5 ) 00369 { 00370 static long int nzone_eval=-1; 00371 static bool lgOscilkase10 = false; 00372 /* this flag says be careful, if kase is this then 00373 * small changes in assumed eden result in large changes 00374 * in returned eden */ 00375 # define KASE_EDEN_NOT_ION 10 00376 if( nzone != nzone_eval ) 00377 { 00378 nzone_eval = nzone; 00379 lgOscilkase10 = false; 00380 } 00381 00382 /* >>chng 03 dec 25, check whether oscillations have occurred */ 00383 if( PreviousChange*(dense.EdenTrue-dense.eden) < 0. && nzone_eval > 6 ) 00384 lgOscilkase10 = true; 00385 00386 /* few of the electrons are due to ions, most are grains and 00387 * or molecules, be very careful */ 00388 change_eden_rel2_tolerance = 0.2; 00389 if( lgOscilkase10 ) 00390 { 00391 /* eden is oscillating, make changes small */ 00392 change_eden_rel2_tolerance = 0.05; 00393 } 00394 kase = KASE_EDEN_NOT_ION; 00395 } 00396 /* was the sign of the change in the electron density changing, 00397 * or has it ever changed? Also use if we are close to converged? */ 00398 /* >>chng 04 aug 04, rm test on eden converged, since fooled solver when 00399 * initially very close to soln, but near eden front, stopped aggressive 00400 * pursuit of true eden, pdr_leiden_v2 */ 00401 else if( PreviousChange*(dense.EdenTrue-dense.eden) < 0. ) 00402 { 00403 /* remember that oscillations are happening */ 00404 conv.lgEdenOscl = true; 00405 nLoopOscil = LoopEden; 00406 /* turn this back on 04 ep 26 */ 00407 /* >>chng 04 sep 27, from * 0.5 to * 0.75 */ 00408 damp = MAX2( 0.05, damp * 0.75 ); 00409 if( PRT_EDEN_OSCIL ) 00410 fprintf(ioQQQ," PRT_EDEN_OSCIIL setting lgEdenOscl true, previous change %e new change %e\n", 00411 PreviousChange,dense.EdenTrue-dense.eden); 00412 00413 /* changes in electron density are changing sign - be careful, 00414 * change_eden_rel2_tolerance multiplies the error tolerance on the electron density, 00415 * largest change allowed is related to this */ 00416 change_eden_rel2_tolerance = 0.5; 00417 /* >>chng 04 aug 09, from 0.5 to 0.25 */ 00418 change_eden_rel2_tolerance = 0.25; 00419 kase = 0; 00420 } 00421 else if( lgFailLastTry ) 00422 { 00423 /* this is case where last evaluation of eden, after solution, 00424 * caused failure */ 00425 change_eden_rel2_tolerance = 0.5; 00426 kase = 1; 00427 } 00428 /* >>chng 03 apr 22, add this branch 00429 * >>chng 04 feb 23, replace conv.lgEdenOscl with lgEden_ever_oscil */ 00430 else if( LoopEden > 4 && !lgEden_ever_oscil && 00431 fabs( (dense.EdenTrue-dense.eden)/dense.eden) > 3.*conv.EdenErrorAllowed ) 00432 { 00433 /* we have gone a few times around this loop, the electron density is 00434 * not oscillating, and we have a long ways to go. Open up the 00435 * allowed change */ 00436 change_eden_rel2_tolerance = 3.; 00437 kase = 2; 00438 } 00439 else if( conv.lgEdenOscl ) 00440 { 00441 /* >>chng 04 aug 14, add this to stop oscil in primal.in */ 00442 /* oscillations have occurred */ 00443 change_eden_rel2_tolerance = 0.25; 00444 kase = 4; 00445 } 00446 else 00447 { 00448 /* stable so far . .. */ 00449 /* change_eden_rel2_tolerance = 2.;*/ 00450 /* had been 2, change to 1 stopped oscillations from developing in 00451 * loc blr grid */ 00452 change_eden_rel2_tolerance = 1.; 00453 kase = 3; 00454 } 00455 00456 /* now remember this change for the next time through the loop */ 00457 PreviousChange = dense.EdenTrue - dense.eden; 00458 00459 /* old difference between assumed and trye */ 00460 eden_assumed_minus_true_old = eden_assumed_minus_true; 00461 00462 /* remember electron density before we update it */ 00463 eden_assumed_old = dense.eden; 00464 00465 /* is the correct electron density - is it too different? 00466 * default on conv.EdenErrorAllowed is 0.01, dyanmics is 0.001 */ 00467 EdenUpperLimit = dense.eden * (1. + damp*conv.EdenErrorAllowed*change_eden_rel2_tolerance); 00468 EdenLowerLimit = dense.eden / (1. + damp*conv.EdenErrorAllowed*change_eden_rel2_tolerance); 00469 00470 # define USE_NEW true 00471 # define PRT_NEW false 00472 if( PRT_NEW && USE_NEW ) 00473 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", 00474 nzone, 00475 conv.nPres2Ioniz, 00476 slope , 00477 slopeold , 00478 slopenew , 00479 dense.eden, 00480 dense.EdenTrue, 00481 (dense.eden-dense.EdenTrue)/SDIV(dense.eden) , 00482 change_eden_rel2_tolerance , 00483 EdenLowerLimit); 00484 if( lgEvalSlop_ls ) 00485 slope = slope_ls; 00486 if( USE_NEW && slope !=0. ) 00487 { 00488 /* slope is d(ne_true) / d(ne_assumed) */ 00489 /* >>chng 04 sep 26, add damp here too */ 00490 eden_proposed = dense.eden + (dense.EdenTrue - dense.eden)/slope_ls * damp; 00491 } 00492 else 00493 { 00494 eden_proposed = dense.EdenTrue; 00495 } 00496 # if 0 00497 { 00498 #include "grainvar.h" 00499 fprintf(ioQQQ,"DEBUG eden grn\t%e\t%e\t%e\t%e\n", 00500 dense.eden, eden_proposed,dense.EdenTrue, -gv.TotalEden); 00501 } 00502 # endif 00503 00504 /* THIS IS IT !!! this is it !!! this is where eden changes. */ 00505 /* get the new electron density */ 00506 /*if( dense.EdenTrue > EdenUpperLimit )*/ 00507 if( eden_proposed > EdenUpperLimit ) 00508 { 00509 /* this branch, proposed eden too big */ 00510 dense.eden = EdenUpperLimit; 00511 } 00512 /*else if( dense.EdenTrue < EdenLowerLimit )*/ 00513 else if( eden_proposed < EdenLowerLimit ) 00514 { 00515 /* proposed eden too small */ 00516 dense.eden = EdenLowerLimit; 00517 } 00518 else 00519 { 00520 /* eden was within fac of the correct value, use it */ 00521 /*dense.eden = dense.EdenTrue;*/ 00522 dense.eden = eden_proposed; 00523 } 00524 00525 if( trace.lgTrace && trace.lgNeBug ) 00526 { 00527 fprintf( ioQQQ, 00528 " ConvEdenIoniz, chg ne to %.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e\n", 00529 dense.eden, eden_assumed_old, dense.eden/eden_assumed_old, dense.EdenTrue ); 00530 } 00531 } 00532 /* we now have the proposed new electron density */ 00533 } 00534 /* >>chng 00 oct 21, did not have this branch before so did not make small changes to 00535 * electron density (smaller than size that determined eden not converged */ 00536 /* this branch electron density was pretty close, we just need to make a small update */ 00537 else 00538 { 00539 /* this is test for whether will call ConvIoniz again - check how close old and correct 00540 * electron densities are before updating EdenTrue */ 00541 if( fabs(dense.eden-dense.EdenTrue)/dense.EdenTrue > conv.EdenErrorAllowed/10. && 00542 /* >>chng 03 nov 29, do not update if we are in not-ion limit of eden, 00543 * this means that molecules and grains have most of the electrons */ 00544 kase !=KASE_EDEN_NOT_ION ) 00545 { 00546 00547 /* now update eden to EdenTrue, since we are so close to is */ 00548 /* >>chng 04 sep 20, from just update to taking into account slope */ 00549 /*dense.eden = dense.EdenTrue;*/ 00550 /* max is to avoid overshoots - we don't want a large correction this late */ 00551 dense.eden = dense.eden + (dense.EdenTrue-dense.eden)/MAX2(1.,slope_ls)*damp; 00552 00553 if( trace.nTrConvg>=3 ) 00554 { 00555 fprintf( ioQQQ, 00556 " ConvEdenIoniz3 converged eden, re-calling ConvIoniz with EdenTrue=%.4e (was %.4e).\n", 00557 dense.EdenTrue , 00558 dense.eden ); 00559 } 00560 00561 /* we have changed the density slightly, so now must reevaluate ionization with this new value */ 00562 /* converge the current level of ionization 00563 * but only do this if change was significant */ 00564 /* >>chng 02 may 31, always call ConvIoniz (basically did so anyway) and remove eden_sum from here*/ 00565 /* , this calls eden_sum which updates EdenTrue */ 00566 if( ConvIoniz() ) 00567 { 00568 return 1; 00569 } 00570 } 00571 else if( trace.nTrConvg>=3 ) 00572 { 00573 fprintf( ioQQQ, 00574 " ConvEdenIoniz3 no need to re-call ConvIoniz since eden did not change much.\n"); 00575 } 00576 00577 /* >>chng 01 mar 14, check whether electron density is no longer converged 00578 * after reevaluating ionization */ 00579 /* >>chng 04 sep 26, had div by min of eden or edentrue, use eden since 00580 * EdenTrue can be negaive when not converged */ 00581 if( 00582 fabs(dense.eden-dense.EdenTrue)/SDIV(dense.eden) >= 00583 conv.EdenErrorAllowed ) 00584 { 00585 /* difference in assumed and true electron density is larger 00586 * than tolerance, we are not yet converged, default is 0.01 */ 00587 conv.lgConvEden = false; 00588 /* >>chng 04 aug 04, do not set oscillations flag if not oscillating */ 00589 /* make steps smaller by setting this flag 00590 conv.lgEdenOscl = true; */ 00591 lgFailLastTry = true; 00592 if( PRT_FAIL_LAST_TRY ) fprintf(ioQQQ,"lgFailLastTry set true\n"); 00593 /* this will stay true through this zone */ 00594 lgEden_ever_oscil = true; 00595 } 00596 } 00597 { 00598 /*@-redef@*/ 00599 /* debug print statement for change in electron density */ 00600 enum {DEBUG_LOC=false}; 00601 /*@+redef@*/ 00602 if( DEBUG_LOC ) 00603 { 00604 fprintf(ioQQQ,"edendebugg %li\t%.2e\t%.2e\t%.2e\t%.2e\n", 00605 nzone,dense.eden ,eden_assumed_old, (dense.eden-eden_assumed_old)/dense.eden, dense.EdenTrue); 00606 } 00607 } 00608 if( trace.lgTrace && trace.lgNeBug ) 00609 { 00610 fprintf( ioQQQ, 00611 " ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e converge=%c\n", 00612 dense.eden, eden_assumed_old, dense.eden/SDIV(eden_assumed_old), dense.EdenTrue ,TorF( conv.lgConvEden)); 00613 } 00614 00615 if( trace.nTrConvg>=3 ) 00616 { 00617 /* these prints all have : delimiter to parse into excel */ 00618 fprintf( ioQQQ, " ConvEdenIoniz3 loop:%4ld ", 00619 LoopEden ); 00620 00621 /* this is flag says whether or not eden/eden has converged */ 00622 if( conv.lgConvEden ) 00623 { 00624 fprintf( ioQQQ, " converged, eden rel error:%g ", 00625 (dense.EdenTrue-dense.eden)/dense.EdenTrue ); 00626 } 00627 else if( eden_assumed_old > SMALLFLOAT ) 00628 { 00629 /* NB - use 8.4 fmt so chng sign not misaleign chngs*/ 00630 fprintf( ioQQQ, " NOT Converged! corr:%8.4f prop:%9.5f ", 00631 (dense.EdenTrue-eden_assumed_old)/eden_assumed_old , 00632 (dense.eden-eden_assumed_old)/eden_assumed_old ); 00633 } 00634 00635 fprintf( ioQQQ, " new ne:%.6e true:%.6e kase:%i slope:%.3e osc:%c", 00636 dense.eden , 00637 dense.EdenTrue , 00638 kase , 00639 slope_ls, 00640 TorF(lgSlope_oscil)); 00641 00642 if( conv.lgEdenOscl ) 00643 { 00644 fprintf( ioQQQ, " EDEN OSCIL1 damp %.4f\n", damp); 00645 } 00646 else 00647 { 00648 fprintf( ioQQQ, "\n"); 00649 } 00650 } 00651 00652 /* loop until converged, or we give up */ 00653 ++LoopEden; 00654 00655 /* this is case where overshoots did occur, but no longer */ 00656 /* >>chng 04 sep 23, introduce this reset */ 00657 /* disable resetting loop unless slope is small */ 00658 if( LoopEden - nLoopOscil >4 && fabs(slope_ls) < 3. ) 00659 { 00660 conv.lgEdenOscl = false; 00661 /* turn this back on 04 ep 26 */ 00662 damp = 1.; 00663 } 00664 00665 /* if last iteration through here and not converged, and no oscillations, 00666 * and no ionization failures , 00667 * then increase limit by 2x */ 00668 if( LoopEden == (LOOPMAX-1) && !conv.lgEdenOscl && conv.nConvIonizFails==0) 00669 /* double the limit, but only one time, and only if no oscillations */ 00670 /* >>chng 04 aug 04, from 2x to 4x, to follow eden changes through eden front */ 00671 LoopLimit *= 4; 00672 00673 } while( !conv.lgConvEden && LoopEden < LoopLimit && !lgAbort ); 00674 } 00675 /* turned on with set eden solver new */ 00676 else if( strcmp( conv.chSolverEden , "new" )== 0 ) 00677 { 00678 /* will be used to look for bracketing in the electron density */ 00679 double PreviousEdenError = 0.; 00680 double CurrentEdenError = 0.; 00681 double CurrentEden = 0.; 00682 double ProposedEden, 00683 FractionChangeEden = 0.; 00684 00685 /* new method of converging electron density */ 00686 LoopEden = 0; 00687 conv.nConvIonizFails = 0; 00688 00689 /* this will be increased by 2x if, at the end, no oscillations have occurred */ 00690 LoopLimit = LOOPMAX; 00691 00692 /* will be set true if sign of change, ever changes */ 00693 conv.lgEdenOscl = false; 00694 if( PRT_EDEN_OSCIL ) fprintf(ioQQQ," PRT_EDEN_OSCIIL setting lgEdenOscl FALSE1\n"); 00695 00696 /* this is relative change in electron density that we will permit - it will become 00697 * smaller each time we bracket the true electron density */ 00698 change_eden_rel2_tolerance = 0.02; 00699 00700 strcpy( conv.chConvIoniz, " NONE!!!!!" ); 00701 00702 /* this is ionization/electron density convergence loop 00703 * keep calling ionize until lgIonDone is true */ 00704 do 00705 { 00706 00707 /* this flag will be set false below if electron density not within tolerance 00708 * after ionization is recomputed */ 00709 conv.lgConvEden = true; 00710 00711 if( trace.nTrConvg>=3 ) 00712 fprintf( ioQQQ, " ConvEdenIoniz3 calling ConvIoniz with eden= %.4e\n",dense.eden); 00713 00714 /* converge the current level of ioinization, this calls eden_sum which updates EdenTrue */ 00715 if( ConvIoniz() ) 00716 { 00717 return 1; 00718 } 00719 00720 /* the history of how we got here 00721 EdenAssumed3 = EdenAssumed2, 00722 EdenObtained3 = EdenObtained2; 00723 EdenAssumed2 = EdenAssumed1, 00724 EdenObtained2 = EdenObtained1,*/ 00725 00726 /* remember the old electron density for possible printout */ 00727 CurrentEden = dense.eden; 00728 00729 /* this is the current error in the electron density */ 00730 CurrentEdenError = dense.EdenTrue - dense.eden; 00731 00732 /* conv.EdenErrorAllowed is allowed error 00733 * the default value of conv.EdenErrorAllowed is 0.01 in zerologic 00734 * is changed with the SET Eden Error command 00735 * eden_assumed_old was incoming value of eden 00736 * EdenTrue is correct value from eden_sum for current ionization 00737 * new eden will be set using these, trying to prevent oscillations */ 00738 00739 /* is error larger than the tolerance we are trying to beat? 00740 * first check is whether error is very small after very first check, since 00741 * ionization soln may not have settled down yet */ 00742 if( 00743 (LoopEden==0 && fabs(CurrentEdenError)/dense.EdenTrue <= conv.EdenErrorAllowed/2.) || 00744 (LoopEden>0 && fabs(CurrentEdenError)/dense.EdenTrue <= conv.EdenErrorAllowed 00745 && FractionChangeEden < conv.EdenErrorAllowed / 2.) ) 00746 break; 00747 00748 /* diference is assumed and true electron density is larger 00749 * than tolerance, we are not yet converged, default is 0.01 */ 00750 conv.lgConvEden = false; 00751 strcpy( conv.chConvIoniz, "Ne big chg" ); 00752 00753 /* SEARCH is true if this is only search for first solution */ 00754 if( conv.lgSearch ) 00755 { 00756 dense.eden = pow(10., 00757 (log10(dense.eden) + log10(dense.EdenTrue))/ 2.); 00758 } 00759 else 00760 { 00761 00762 /* was the sign of the change in the electron density changing, 00763 * if so then we have bracked the target */ 00764 if( PreviousEdenError*CurrentEdenError < 0. ) 00765 { 00766 /* we have bracketed the correct electron density, 00767 * make changes smaller, also solve linear equation for zero error */ 00768 slope = (PreviousEdenError - CurrentEdenError ) / 00769 ( eden_assumed_old - CurrentEden ); 00770 00771 ProposedEden = eden_assumed_old - PreviousEdenError / slope; 00772 /* as sanity check, this must be between the two limits we examined, 00773 * since zero was between them */ 00774 ASSERT( ProposedEden >= MIN2( eden_assumed_old , CurrentEden ) ); 00775 ASSERT( ProposedEden <= MAX2( eden_assumed_old , CurrentEden ) ); 00776 00777 change_eden_rel2_tolerance *= 0.25; 00778 if( trace.nTrConvg>=3 ) 00779 fprintf( ioQQQ, 00780 " ConvEdenIoniz3 bracketed electron density factor now %.3e error is %.4f proposed ne %.4e\n", 00781 change_eden_rel2_tolerance, 00782 (dense.EdenTrue-dense.eden)/dense.EdenTrue, ProposedEden); 00783 } 00784 else 00785 { 00786 /* did not bracket the target, keep moving in its direction */ 00787 00788 /* the correct electron density - is it too different? */ 00789 EdenUpperLimit = dense.eden * (1. + change_eden_rel2_tolerance); 00790 EdenLowerLimit = dense.eden / (1. + change_eden_rel2_tolerance); 00791 00792 /* get the new electron density */ 00793 if( dense.EdenTrue > EdenUpperLimit ) 00794 { 00795 /* this branch, proposed eden too big */ 00796 ProposedEden = EdenUpperLimit; 00797 } 00798 else if( dense.EdenTrue < EdenLowerLimit ) 00799 { 00800 /* proposed eden too small */ 00801 ProposedEden = EdenLowerLimit; 00802 } 00803 else 00804 { 00805 /* eden was within fac of the correct value, use it */ 00806 ProposedEden = dense.EdenTrue; 00807 } 00808 if( trace.nTrConvg>=3 ) 00809 fprintf( ioQQQ, 00810 " ConvEdenIoniz3 elec dens fac %.3e err is %.4f prop ne %.4e cor ne %.4e slope=%.2e\n", 00811 change_eden_rel2_tolerance, 00812 (dense.EdenTrue-dense.eden)/dense.EdenTrue , 00813 ProposedEden, 00814 dense.EdenTrue,slope); 00815 } 00816 00817 /* now remember this change for the next time through the loop */ 00818 PreviousEdenError = CurrentEdenError; 00819 eden_assumed_old = CurrentEden; 00820 FractionChangeEden = fabs( dense.eden - ProposedEden ) / dense.EdenTrue; 00821 dense.eden = ProposedEden; 00822 00823 if( trace.lgTrace && trace.lgNeBug ) 00824 { 00825 fprintf( ioQQQ, 00826 " ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e\n", 00827 dense.eden, eden_assumed_old, dense.eden/eden_assumed_old, dense.EdenTrue ); 00828 } 00829 } 00830 /* loop until we give up -- normal exit is break in center of loop */ 00831 ++LoopEden; 00832 } while( LoopEden < LoopLimit && !lgAbort ); 00833 00834 /* we now have the proposed new electron density */ 00835 /* we have changed the density slightly, so now must reevaluate ionization with this new value */ 00836 /* converge the current level of ioinization 00837 * but only do this if change was significant */ 00838 if( trace.nTrConvg>=3 ) 00839 fprintf( ioQQQ, 00840 " ConvEdenIoniz3 converged eden, current error is %.4f, calling ConvIoniz with final density=EdenTrue=%.4e\n", 00841 (dense.EdenTrue - dense.eden)/dense.EdenTrue,dense.EdenTrue); 00842 dense.eden = dense.EdenTrue; 00843 /* , this calls eden_sum which updates EdenTrue */ 00844 /* >>chng 02 may 31, always call ConvIoniz (basically did so anyway) and remove eden_sum from here*/ 00845 if( ConvIoniz() ) 00846 { 00847 if( trace.nTrConvg>=3 ) 00848 fprintf( ioQQQ, 00849 " ConvEdenIoniz3 eden no longer converged eden, current error is %.4f\n", 00850 (dense.EdenTrue - dense.eden)/dense.EdenTrue); 00851 } 00852 00853 /* >>chng 01 mar 14, check whether electron density is no longer converged 00854 * after reevaluating ionization */ 00855 if( 00856 fabs(dense.eden-dense.EdenTrue)/dense.EdenTrue > 00857 conv.EdenErrorAllowed ) 00858 { 00859 /* diference is assumed and true electron density is larger 00860 * than tolerance, we are not yet converged, default is 0.01 */ 00861 conv.lgConvEden = false; 00862 if( trace.nTrConvg>=3 ) 00863 fprintf( ioQQQ, 00864 " ConvEdenIoniz3 setting eden not converged, error %.4f, exiting\n", 00865 (dense.eden-dense.EdenTrue)/dense.EdenTrue ); 00866 } 00867 else 00868 { 00869 conv.lgConvEden = true; 00870 } 00871 00872 if( trace.lgTrace && trace.lgNeBug ) 00873 { 00874 fprintf( ioQQQ, 00875 " ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e converge=%c\n", 00876 dense.eden, eden_assumed_old, dense.eden/eden_assumed_old, dense.EdenTrue ,TorF( conv.lgConvEden)); 00877 } 00878 00879 if( trace.nTrConvg>=3 ) 00880 { 00881 fprintf( ioQQQ, " ConvEdenIoniz3 %4ld new ne=%12.4e true=%12.4e", 00882 LoopEden, dense.eden , dense.EdenTrue ); 00883 00884 /* this is flag says whether or not eden/eden has converged */ 00885 if( conv.lgConvEden ) 00886 { 00887 fprintf( ioQQQ, " converged, eden rel error is %g ", 00888 (dense.EdenTrue-dense.eden)/dense.EdenTrue ); 00889 } 00890 else 00891 { 00892 fprintf( ioQQQ, " NOCONV corr:%7.3f prop:%7.3f ", 00893 (dense.EdenTrue-eden_assumed_old)/eden_assumed_old , 00894 (dense.eden-eden_assumed_old)/eden_assumed_old ); 00895 } 00896 if( conv.lgEdenOscl ) 00897 fprintf( ioQQQ, " EDEN OSCIL2 \n"); 00898 } 00899 00900 } 00901 else 00902 { 00903 fprintf( ioQQQ, "ConvEdenIoniz finds insane solver %s \n" , conv.chSolverEden ); 00904 ShowMe(); 00905 } 00906 ++LoopBig; 00907 } 00908 00909 if( trace.nTrConvg>=3 ) 00910 { 00911 fprintf( ioQQQ, " ConvEdenIoniz3 exits, lgConvEden=%1c entry eden %.4e -> %.4e rel change %.3f\n" , 00912 TorF(conv.lgConvEden ), 00913 EdenEntry , 00914 dense.eden , 00915 (EdenEntry-dense.eden)/SDIV( dense.eden ) ); 00916 } 00917 else if( trace.lgTrace ) 00918 { 00919 fprintf( ioQQQ, " ConvEdenIoniz exits zn %.2f lgConvEden=%1c entry eden %.4e -> %.4e rel change %.3f\n" , 00920 fnzone, 00921 TorF(conv.lgConvEden ), 00922 EdenEntry , 00923 dense.eden , 00924 (EdenEntry-dense.eden)/SDIV( dense.eden ) ); 00925 } 00926 00927 return 0; 00928 } 00929 00930 /* returns true if electron density is converged */ 00931 bool lgConvEden(void) 00932 { 00933 bool lgRet; 00934 00935 /* >>chng 04 sep 26, div by eden, not min of eden and edentrue since latter now poss neg */ 00936 if( fabs(dense.eden-dense.EdenTrue)/SDIV(dense.eden) >= 00937 conv.EdenErrorAllowed ) 00938 { 00939 lgRet = false; 00940 conv.lgConvEden = false; 00941 } 00942 else 00943 { 00944 lgRet = true; 00945 conv.lgConvEden = true; 00946 } 00947 return lgRet; 00948 }