cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
conv_eden_ioniz.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*ConvEdenIoniz called by ConvTempIonz, calls ConvIoniz solving for eden */
4 /*lgConvEden returns true if electron density is converged */
5 #include "cddefines.h"
6 #include "dense.h"
7 #include "trace.h"
8 #include "thermal.h"
9 #include "thirdparty.h"
10 #include "phycon.h"
11 #include "conv.h"
12 /* limit to how many loops we will do */
13 /* >>chng 04 sep 27, from 25 to 35 */
14 static const int LOOPMAX = 35;
15 
16 /*ConvEdenIoniz called by ConvTempIonz, calls ConvIoniz solving for eden
17  * returns 0 if ok, 1 if abort */
18 int ConvEdenIoniz(void)
19 {
20  long int LoopEden ,
21  /* this will be LOOPMAX at first, then doubled if no oscillations occur*/
22  LoopLimit ,
23  LoopBig;
24  bool lgFailLastTry;
25  int kase = -1;
26  static bool lgSlope_oscil=false;
27 
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;
34  double EdenEntry;
35  static bool lgMustMalloc_temp_history = true;
36  static double slope=-1.;
37  double
38  /* upper and lower limits to the range of the change we want to consider */
39  EdenUpperLimit,
40  EdenLowerLimit ,
41  change_eden_rel2_tolerance ,
42  PreviousChange;
43 
44  DEBUG_ENTRY( "ConvEdenIoniz()" );
45 
46  /* this routine is called by ConvTempIonz, it calls ConvIonize
47  * and changes the electron density until it converges */
48 
49  /* save entry value of eden */
50  EdenEntry = dense.eden;
51 
52  if( trace.lgTrace )
53  {
54  fprintf( ioQQQ, " \n" );
55  fprintf( ioQQQ, " ConvEdenIoniz entered \n" );
56  }
57 
58  if( trace.nTrConvg>=3 )
59  {
60  fprintf( ioQQQ,
61  " ConvEdenIoniz3 called, entering eden loop using solver %sw.\n",
63  }
64 
65  /* keep track of temperature solver in this zone -
66  * conv.nTotalIoniz is zero in first call of new iteration */
67  if( !conv.nTotalIoniz )
68  {
69  /* one time initialization of variables */
70  conv.hist_temp_ntemp = -1;
71  conv.hist_temp_nzone = -1;
72  /*>>chng 06 jun 25, do not reset upper limit to number of
73  * points saved - only set this when malloc or remalloc */
74  /*conv.hist_temp_limit = 0;*/
75  }
76  /* this will save the history of density - pressure relationship
77  * for the current zone */
79  {
80  /* first time in this zone - reset counters */
82  /* this counter will be updated after vars are saved so will be
83  * total number of saved values */
85  }
86  /* do we need to allocate, or reallocate, memory for the vectors
87  * check initial allocation first */
88  /*>>chng 06 jun 25, do not test on this limit for need to malloc
89  * rather set static flag in this routine - fix of memory leak
90  * discovered by PvH - routine entered with conv.nTotalIoniz zero
91  * during first call of each iteration so memory was malloced without
92  * freeing previous memory - memory leak resulted */
93  /*if( conv.hist_temp_limit==0 )*/
94  if( lgMustMalloc_temp_history )
95  {
96  lgMustMalloc_temp_history = false;
97  conv.hist_temp_limit = 200;
98  conv.hist_temp_temp = (double *)MALLOC( (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double) ) );
99  conv.hist_temp_heat = (double *)MALLOC( (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double) ) );
100  conv.hist_temp_cool = (double *)MALLOC( (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double) ) );
101  conv.chNotConverged = (char *)MALLOC( (size_t)((unsigned long)INPUT_LINE_LENGTH*sizeof(char) ) );
102  conv.chConvEden = (char *)MALLOC( (size_t)((unsigned long)INPUT_LINE_LENGTH*sizeof(char) ) );
103  conv.chConvIoniz = (char *)MALLOC( (size_t)((unsigned long)INPUT_LINE_LENGTH*sizeof(char) ) );
104  strcpy( conv.chNotConverged, "none" );
105  strcpy( conv.chConvEden, "none" );
106  strcpy( conv.chConvIoniz, "none" );
107  }
108  /* ran out of space - allocate more */
110  {
111  conv.hist_temp_limit *= 3;
112  conv.hist_temp_temp = (double *)REALLOC( conv.hist_temp_temp, (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double)));
113  conv.hist_temp_heat = (double *)REALLOC( conv.hist_temp_heat, (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double)));
114  conv.hist_temp_cool = (double *)REALLOC( conv.hist_temp_cool, (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double)));
115  }
116 
117  /* >>chng 04 feb 11, add option to remember current density and pressure */
122 
123  if( conv.nPres2Ioniz < 5 )
124  lgEden_ever_oscil = false;
125 
126  /* this flag says wether we converged the density but failed after final tweek,
127  * need to use smaller steps */
128 # define PRT_FAIL_LAST_TRY false
129  lgFailLastTry = false;
130  if( PRT_FAIL_LAST_TRY ) fprintf(ioQQQ,"lgFailLastTry set false\n");
131 
132  /* >>chng 03 mar 17, add this big loop, only go through one time if fully
133  * converged at end, but if eden jumps off in final touchup,
134  * do whole thing again */
135  LoopBig = 0;
136  while( LoopBig==0 || (lgFailLastTry && LoopBig==1 ) )
137  {
138  /* the old default solver */
139  if( strcmp( conv.chSolverEden , "simple" )== 0 )
140  {
141  static double damp;
142  long int nLoopOscil;
143  LoopEden = 0;
144  conv.nConvIonizFails = 0;
145 
146  /* this will be increased by 2x if, at the end, no oscillations have occurred */
147  /* >>chng 04 aug 04, from 2x to 4x, to following through on eden front */
148  LoopLimit = LOOPMAX;
149 
150  /* will be set true if sign of change, ever changes */
151  conv.lgEdenOscl = false;
152  nLoopOscil = 0;
153 # define PRT_EDEN_OSCIL false
154  if( PRT_EDEN_OSCIL ) fprintf(ioQQQ," PRT_EDEN_OSCIIL setting lgEdenOscl FALSE2\n");
155  lgFailLastTry = false;
156  if( PRT_FAIL_LAST_TRY ) fprintf(ioQQQ,"lgFailLastTry set false\n");
157 
158  /* will be used to look for oscillations in the electron density */
159  PreviousChange = 0.;
160 
161  strcpy( conv.chConvIoniz, " NONE!!!!!" );
162 
163  /* we have not yet */
164  eden_assumed_old = 0.;
165  eden_assumed_minus_true_old = 0.;
166 
167  damp = 1.;
168  /* if working with zone 0 (nzone=0) reset slope to zero */
169  if( !nzone )
170  slope = 0.;
171 
172  /* this is electron density convergence loop */
173  do
174  {
175  double slopenew , slopeold;
176  /* this flag will be set false below if electron density not within tolerance
177  * after ionization is recomputed */
178  conv.lgConvEden = true;
179 
180  /* converge the current level of ionization, this calls eden_sum which updates EdenTrue */
181  if( ConvIoniz() || lgAbort )
182  {
183  return 1;
184  }
185  /* this is the new error in the electron density, the difference between
186  * the assumed and true electron density - we want this to be zero */
187  eden_assumed_minus_true = dense.eden - dense.EdenTrue;
188 
189  {
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 ,
194  *stdarray;
195  static long int nEval=-1 , nzoneUsed=-1, iterUsed=-1;
196  if( lgMustMalloc_eden_error_history )
197  {
198  lgMustMalloc_eden_error_history = false;
199  lgSlope_oscil = false;
200  eden_error_history =
201  (double*)MALLOC(sizeof(double)*(unsigned)limEdenHistory );
202  eden_assumed_history =
203  (double*)MALLOC(sizeof(double)*(unsigned)limEdenHistory );
204  stdarray =
205  (double*)MALLOC(sizeof(double)*(unsigned)limEdenHistory );
206  }
207  if( nzoneUsed!=nzone || iterUsed!=iteration )
208  {
209  /* first evaluation this zone */
210  if( nzone==0 )
211  {
212  lgEvalSlop_ls = true;
213  slope_ls = 1.;
214  }
215  /* reset some vars at start of new zone */
216  nEval = 0;
217  nzoneUsed = nzone;
218  iterUsed = iteration;
219  lgSlope_oscil = false;
220  slope_ls_lastzone = slope_ls;
221  }
222  /* do not evaluate this during search phase or if slope is oscillating */
223  else if( !conv.lgSearch && !lgSlope_oscil )
224  {
225  /* may cycle through here with exactly same edensity before tripping if */
226  int ip = MAX2(0,nEval-1);
227  if(
228  /* no need to do this if already converged */
230  /* do it if never evaluated */
231  (!nEval ||
232  /* this test - don't want tiny corrections, within the tolerance, which
233  * can have noise since at the level code is converging, affecting the slope */
234  fabs((dense.eden - dense.EdenTrue)-eden_error_history[ip]*dense.gas_phase[ipHYDROGEN] )/dense.EdenTrue > 1e-5 ||
235  fabs(dense.eden-eden_assumed_history[ip]*dense.gas_phase[ipHYDROGEN])/ dense.EdenTrue > 1e-5 ))
236  {
237  /* use relative fraction for constant pressure case */
238  eden_error_history[nEval] = (dense.eden - dense.EdenTrue)/dense.gas_phase[ipHYDROGEN];
239  eden_assumed_history[nEval] = dense.eden/dense.gas_phase[ipHYDROGEN];
240  stdarray[nEval] = 0.;
241  /* >>chng 05 feb 22, only update nEval if residuals are non-zero -
242  * this happens in constant temperature models - happened in optimization run */
243  if( eden_error_history[nEval] != 0. )
244  {
245  ++nEval;
246  lgHistUpdate = true;
247  }
248  /*fprintf(ioQQQ,"DEBUG history\t%li\t%e\t%e\t%e\n",
249  nEval,
250  dense.gas_phase[ipHYDROGEN],
251  eden_assumed_history[nEval],
252  eden_error_history[nEval]);*/
253  }
254  if( nEval>=limEdenHistory )
255  {
256  /* used a lot of points - must create more space */
257  limEdenHistory *=10;
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));
261  }
262  if( nEval>3 && lgHistUpdate )
263  {
264  double y_intercept, stdx, stdy, slope_save;
265  slope_save = slope_ls;
266  /* enough data to do least squares */
267  if( linfit(nEval,
268  eden_assumed_history,
269  eden_error_history,
270  y_intercept,
271  stdy,
272  slope_ls,
273  stdx) )
274  {
275  int i;
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 )
279  {
280  fprintf(ioQQQ,"%.3e\t%.4e\n",
281  eden_error_history[i] ,
282  stdarray[i] );
283  }
284  fprintf(ioQQQ,"\n");
285  ShowMe();
286  cdEXIT(EXIT_FAILURE);
287  }
288  lgEvalSlop_ls = true;
289  /* check for slope changing sign - if does, use last zone's slope */
290  if( slope_ls*slope_save<0. )
291  {
292  slope_ls = slope_ls_lastzone;
293  lgSlope_oscil = true;
294  }
295  /*fprintf(ioQQQ,"DEBUG LS\t%.2f\t%li\t%.3e\t%.3e\t%.3e\n",
296  fnzone ,
297  nEval ,
298  slope_ls,
299  y_intercept ,
300  slope);*/
301  }
302  }
303  }
304  slopeold = slope;
305  /* may not be set below, but could print it. set to zero as flag for this case */
306  slopenew = 0.;
307  if( fabs(eden_assumed_minus_true_old) > 0. )
308  {
309  if( fabs( eden_assumed_old-dense.eden ) > SMALLFLOAT )
310  {
311 # define OLDFAC 0.75
312  slopenew = (eden_assumed_minus_true_old - eden_assumed_minus_true) / (eden_assumed_old-dense.eden );
313  /* >>chng 04 sep 15 do not update slope if changing sign */
314 # if 0
315  if( slope !=0. && slope*slopenew>=0.)
316  slope = OLDFAC*slope + (1.-OLDFAC)*slopenew;
317  else if( slope==0.)
318  slope = slopenew;
319 # endif
320  if( slope !=0. )
321  slope = OLDFAC*slope + (1.-OLDFAC)*slopenew;
322  else
323  slope = slopenew;
324 # undef OLDFAC
325  }
326  /*else
327  slope = 0.;*/
328  }
329  /* following block is to not let electron density change by
330  * too much
331  * conv.EdenErrorAllowed is allowed error
332  * the default value of conv.EdenErrorAllowed is 0.01 in zerologic
333  * is changed with the SET Eden Error command
334  * eden_assumed_old was incoming value of eden
335  * EdenTrue is correct value from eden_sum for current ionization
336  * new eden will be set using these, trying to prevent oscillations */
337 
338  /* is error larger than the tolerance we are trying to beat? */
339  if( fabs(dense.eden-dense.EdenTrue)/SDIV(dense.eden) >=
341  {
342  double eden_proposed;
343  /* difference is assumed and true electron density is larger
344  * than tolerance, we are not yet converged, default is 0.01 */
345  conv.lgConvEden = false;
346  strcpy( conv.chConvIoniz, "Ne big chg" );
347 
348  /* these are needed since will take log below */
350 
351  /* SEARCH is true if this is only search for first solution */
352  /* EdenTrue is allowed to go negative during search for soln - not physical,
353  * but can happen in math search for soln */
355  {
356  dense.eden = pow(10.,
357  (log10(dense.eden) + log10(dense.EdenTrue))/ 2.);
358  }
359  else
360  {
361  /* >>chng 03 jul 07, add case for grains contain significant
362  * fraction of electrons */
363  /* >>chng 03 nov 28, add this test on fraction of electrons from ions,
364  * this branch is to identify limit where molecules and grains have
365  * most of the electrons */
366  /* >>chng 04 sep 14, change from all ions to just metals */
367  /* this was patch to fix steep slope */
368  if( 0 && dense.eden_from_metals > 0.5 )
369  {
370  static long int nzone_eval=-1;
371  static bool lgOscilkase10 = false;
372  /* this flag says be careful, if kase is this then
373  * small changes in assumed eden result in large changes
374  * in returned eden */
375 # define KASE_EDEN_NOT_ION 10
376  if( nzone != nzone_eval )
377  {
378  nzone_eval = nzone;
379  lgOscilkase10 = false;
380  }
381 
382  /* >>chng 03 dec 25, check whether oscillations have occurred */
383  if( PreviousChange*(dense.EdenTrue-dense.eden) < 0. && nzone_eval > 6 )
384  lgOscilkase10 = true;
385 
386  /* few of the electrons are due to ions, most are grains and
387  * or molecules, be very careful */
388  change_eden_rel2_tolerance = 0.2;
389  if( lgOscilkase10 )
390  {
391  /* eden is oscillating, make changes small */
392  change_eden_rel2_tolerance = 0.05;
393  }
394  kase = KASE_EDEN_NOT_ION;
395  }
396  /* was the sign of the change in the electron density changing,
397  * or has it ever changed? Also use if we are close to converged? */
398  /* >>chng 04 aug 04, rm test on eden converged, since fooled solver when
399  * initially very close to soln, but near eden front, stopped aggressive
400  * pursuit of true eden, pdr_leiden_v2 */
401  else if( PreviousChange*(dense.EdenTrue-dense.eden) < 0. )
402  {
403  /* remember that oscillations are happening */
404  conv.lgEdenOscl = true;
405  nLoopOscil = LoopEden;
406  /* turn this back on 04 ep 26 */
407  /* >>chng 04 sep 27, from * 0.5 to * 0.75 */
408  damp = MAX2( 0.05, damp * 0.75 );
409  if( PRT_EDEN_OSCIL )
410  fprintf(ioQQQ," PRT_EDEN_OSCIIL setting lgEdenOscl true, previous change %e new change %e\n",
411  PreviousChange,dense.EdenTrue-dense.eden);
412 
413  /* changes in electron density are changing sign - be careful,
414  * change_eden_rel2_tolerance multiplies the error tolerance on the electron density,
415  * largest change allowed is related to this */
416  change_eden_rel2_tolerance = 0.5;
417  /* >>chng 04 aug 09, from 0.5 to 0.25 */
418  change_eden_rel2_tolerance = 0.25;
419  kase = 0;
420  }
421  else if( lgFailLastTry )
422  {
423  /* this is case where last evaluation of eden, after solution,
424  * caused failure */
425  change_eden_rel2_tolerance = 0.5;
426  kase = 1;
427  }
428  /* >>chng 03 apr 22, add this branch
429  * >>chng 04 feb 23, replace conv.lgEdenOscl with lgEden_ever_oscil */
430  else if( LoopEden > 4 && !lgEden_ever_oscil &&
432  {
433  /* we have gone a few times around this loop, the electron density is
434  * not oscillating, and we have a long ways to go. Open up the
435  * allowed change */
436  change_eden_rel2_tolerance = 3.;
437  kase = 2;
438  }
439  else if( conv.lgEdenOscl )
440  {
441  /* >>chng 04 aug 14, add this to stop oscil in primal.in */
442  /* oscillations have occurred */
443  change_eden_rel2_tolerance = 0.25;
444  kase = 4;
445  }
446  else
447  {
448  /* stable so far . .. */
449  /* change_eden_rel2_tolerance = 2.;*/
450  /* had been 2, change to 1 stopped oscillations from developing in
451  * loc blr grid */
452  change_eden_rel2_tolerance = 1.;
453  kase = 3;
454  }
455 
456  /* now remember this change for the next time through the loop */
457  PreviousChange = dense.EdenTrue - dense.eden;
458 
459  /* old difference between assumed and trye */
460  eden_assumed_minus_true_old = eden_assumed_minus_true;
461 
462  /* remember electron density before we update it */
463  eden_assumed_old = dense.eden;
464 
465  /* is the correct electron density - is it too different?
466  * default on conv.EdenErrorAllowed is 0.01, dyanmics is 0.001 */
467  EdenUpperLimit = dense.eden * (1. + damp*conv.EdenErrorAllowed*change_eden_rel2_tolerance);
468  EdenLowerLimit = dense.eden / (1. + damp*conv.EdenErrorAllowed*change_eden_rel2_tolerance);
469 
470 # define USE_NEW true
471 # define PRT_NEW false
472  if( PRT_NEW && USE_NEW )
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",
474  nzone,
476  slope ,
477  slopeold ,
478  slopenew ,
479  dense.eden,
480  dense.EdenTrue,
482  change_eden_rel2_tolerance ,
483  EdenLowerLimit);
484  if( lgEvalSlop_ls )
485  slope = slope_ls;
486  if( USE_NEW && slope !=0. )
487  {
488  /* slope is d(ne_true) / d(ne_assumed) */
489  /* >>chng 04 sep 26, add damp here too */
490  eden_proposed = dense.eden + (dense.EdenTrue - dense.eden)/slope_ls * damp;
491  }
492  else
493  {
494  eden_proposed = dense.EdenTrue;
495  }
496 # if 0
497  {
498 #include "grainvar.h"
499  fprintf(ioQQQ,"DEBUG eden grn\t%e\t%e\t%e\t%e\n",
500  dense.eden, eden_proposed,dense.EdenTrue, -gv.TotalEden);
501  }
502 # endif
503 
504  /* THIS IS IT !!! this is it !!! this is where eden changes. */
505  /* get the new electron density */
506  /*if( dense.EdenTrue > EdenUpperLimit )*/
507  if( eden_proposed > EdenUpperLimit )
508  {
509  /* this branch, proposed eden too big */
510  dense.eden = EdenUpperLimit;
511  }
512  /*else if( dense.EdenTrue < EdenLowerLimit )*/
513  else if( eden_proposed < EdenLowerLimit )
514  {
515  /* proposed eden too small */
516  dense.eden = EdenLowerLimit;
517  }
518  else
519  {
520  /* eden was within fac of the correct value, use it */
521  /*dense.eden = dense.EdenTrue;*/
522  dense.eden = eden_proposed;
523  }
524 
525  if( trace.lgTrace && trace.lgNeBug )
526  {
527  fprintf( ioQQQ,
528  " ConvEdenIoniz, chg ne to %.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e\n",
529  dense.eden, eden_assumed_old, dense.eden/eden_assumed_old, dense.EdenTrue );
530  }
531  }
532  /* we now have the proposed new electron density */
533  }
534  /* >>chng 00 oct 21, did not have this branch before so did not make small changes to
535  * electron density (smaller than size that determined eden not converged */
536  /* this branch electron density was pretty close, we just need to make a small update */
537  else
538  {
539  /* this is test for whether will call ConvIoniz again - check how close old and correct
540  * electron densities are before updating EdenTrue */
542  /* >>chng 03 nov 29, do not update if we are in not-ion limit of eden,
543  * this means that molecules and grains have most of the electrons */
544  kase !=KASE_EDEN_NOT_ION )
545  {
546 
547  /* now update eden to EdenTrue, since we are so close to is */
548  /* >>chng 04 sep 20, from just update to taking into account slope */
549  /*dense.eden = dense.EdenTrue;*/
550  /* max is to avoid overshoots - we don't want a large correction this late */
551  dense.eden = dense.eden + (dense.EdenTrue-dense.eden)/MAX2(1.,slope_ls)*damp;
552 
553  if( trace.nTrConvg>=3 )
554  {
555  fprintf( ioQQQ,
556  " ConvEdenIoniz3 converged eden, re-calling ConvIoniz with EdenTrue=%.4e (was %.4e).\n",
557  dense.EdenTrue ,
558  dense.eden );
559  }
560 
561  /* we have changed the density slightly, so now must reevaluate ionization with this new value */
562  /* converge the current level of ionization
563  * but only do this if change was significant */
564  /* >>chng 02 may 31, always call ConvIoniz (basically did so anyway) and remove eden_sum from here*/
565  /* , this calls eden_sum which updates EdenTrue */
566  if( ConvIoniz() )
567  {
568  return 1;
569  }
570  }
571  else if( trace.nTrConvg>=3 )
572  {
573  fprintf( ioQQQ,
574  " ConvEdenIoniz3 no need to re-call ConvIoniz since eden did not change much.\n");
575  }
576 
577  /* >>chng 01 mar 14, check whether electron density is no longer converged
578  * after reevaluating ionization */
579  /* >>chng 04 sep 26, had div by min of eden or edentrue, use eden since
580  * EdenTrue can be negaive when not converged */
581  if(
582  fabs(dense.eden-dense.EdenTrue)/SDIV(dense.eden) >=
584  {
585  /* difference in assumed and true electron density is larger
586  * than tolerance, we are not yet converged, default is 0.01 */
587  conv.lgConvEden = false;
588  /* >>chng 04 aug 04, do not set oscillations flag if not oscillating */
589  /* make steps smaller by setting this flag
590  conv.lgEdenOscl = true; */
591  lgFailLastTry = true;
592  if( PRT_FAIL_LAST_TRY ) fprintf(ioQQQ,"lgFailLastTry set true\n");
593  /* this will stay true through this zone */
594  lgEden_ever_oscil = true;
595  }
596  }
597  {
598  /*@-redef@*/
599  /* debug print statement for change in electron density */
600  enum {DEBUG_LOC=false};
601  /*@+redef@*/
602  if( DEBUG_LOC )
603  {
604  fprintf(ioQQQ,"edendebugg %li\t%.2e\t%.2e\t%.2e\t%.2e\n",
605  nzone,dense.eden ,eden_assumed_old, (dense.eden-eden_assumed_old)/dense.eden, dense.EdenTrue);
606  }
607  }
608  if( trace.lgTrace && trace.lgNeBug )
609  {
610  fprintf( ioQQQ,
611  " ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e converge=%c\n",
612  dense.eden, eden_assumed_old, dense.eden/SDIV(eden_assumed_old), dense.EdenTrue ,TorF( conv.lgConvEden));
613  }
614 
615  if( trace.nTrConvg>=3 )
616  {
617  /* these prints all have : delimiter to parse into excel */
618  fprintf( ioQQQ, " ConvEdenIoniz3 loop:%4ld ",
619  LoopEden );
620 
621  /* this is flag says whether or not eden/eden has converged */
622  if( conv.lgConvEden )
623  {
624  fprintf( ioQQQ, " converged, eden rel error:%g ",
626  }
627  else if( eden_assumed_old > SMALLFLOAT )
628  {
629  /* NB - use 8.4 fmt so chng sign not misaleign chngs*/
630  fprintf( ioQQQ, " NOT Converged! corr:%8.4f prop:%9.5f ",
631  (dense.EdenTrue-eden_assumed_old)/eden_assumed_old ,
632  (dense.eden-eden_assumed_old)/eden_assumed_old );
633  }
634 
635  fprintf( ioQQQ, " new ne:%.6e true:%.6e kase:%i slope:%.3e osc:%c",
636  dense.eden ,
637  dense.EdenTrue ,
638  kase ,
639  slope_ls,
640  TorF(lgSlope_oscil));
641 
642  if( conv.lgEdenOscl )
643  {
644  fprintf( ioQQQ, " EDEN OSCIL1 damp %.4f\n", damp);
645  }
646  else
647  {
648  fprintf( ioQQQ, "\n");
649  }
650  }
651 
652  /* loop until converged, or we give up */
653  ++LoopEden;
654 
655  /* this is case where overshoots did occur, but no longer */
656  /* >>chng 04 sep 23, introduce this reset */
657  /* disable resetting loop unless slope is small */
658  if( LoopEden - nLoopOscil >4 && fabs(slope_ls) < 3. )
659  {
660  conv.lgEdenOscl = false;
661  /* turn this back on 04 ep 26 */
662  damp = 1.;
663  }
664 
665  /* if last iteration through here and not converged, and no oscillations,
666  * and no ionization failures ,
667  * then increase limit by 2x */
668  if( LoopEden == (LOOPMAX-1) && !conv.lgEdenOscl && conv.nConvIonizFails==0)
669  /* double the limit, but only one time, and only if no oscillations */
670  /* >>chng 04 aug 04, from 2x to 4x, to follow eden changes through eden front */
671  LoopLimit *= 4;
672 
673  } while( !conv.lgConvEden && LoopEden < LoopLimit && !lgAbort );
674  }
675  /* turned on with set eden solver new */
676  else if( strcmp( conv.chSolverEden , "new" )== 0 )
677  {
678  /* will be used to look for bracketing in the electron density */
679  double PreviousEdenError = 0.;
680  double CurrentEdenError = 0.;
681  double CurrentEden = 0.;
682  double ProposedEden,
683  FractionChangeEden = 0.;
684 
685  /* new method of converging electron density */
686  LoopEden = 0;
687  conv.nConvIonizFails = 0;
688 
689  /* this will be increased by 2x if, at the end, no oscillations have occurred */
690  LoopLimit = LOOPMAX;
691 
692  /* will be set true if sign of change, ever changes */
693  conv.lgEdenOscl = false;
694  if( PRT_EDEN_OSCIL ) fprintf(ioQQQ," PRT_EDEN_OSCIIL setting lgEdenOscl FALSE1\n");
695 
696  /* this is relative change in electron density that we will permit - it will become
697  * smaller each time we bracket the true electron density */
698  change_eden_rel2_tolerance = 0.02;
699 
700  strcpy( conv.chConvIoniz, " NONE!!!!!" );
701 
702  /* this is ionization/electron density convergence loop
703  * keep calling ionize until lgIonDone is true */
704  do
705  {
706 
707  /* this flag will be set false below if electron density not within tolerance
708  * after ionization is recomputed */
709  conv.lgConvEden = true;
710 
711  if( trace.nTrConvg>=3 )
712  fprintf( ioQQQ, " ConvEdenIoniz3 calling ConvIoniz with eden= %.4e\n",dense.eden);
713 
714  /* converge the current level of ioinization, this calls eden_sum which updates EdenTrue */
715  if( ConvIoniz() )
716  {
717  return 1;
718  }
719 
720  /* the history of how we got here
721  EdenAssumed3 = EdenAssumed2,
722  EdenObtained3 = EdenObtained2;
723  EdenAssumed2 = EdenAssumed1,
724  EdenObtained2 = EdenObtained1,*/
725 
726  /* remember the old electron density for possible printout */
727  CurrentEden = dense.eden;
728 
729  /* this is the current error in the electron density */
730  CurrentEdenError = dense.EdenTrue - dense.eden;
731 
732  /* conv.EdenErrorAllowed is allowed error
733  * the default value of conv.EdenErrorAllowed is 0.01 in zerologic
734  * is changed with the SET Eden Error command
735  * eden_assumed_old was incoming value of eden
736  * EdenTrue is correct value from eden_sum for current ionization
737  * new eden will be set using these, trying to prevent oscillations */
738 
739  /* is error larger than the tolerance we are trying to beat?
740  * first check is whether error is very small after very first check, since
741  * ionization soln may not have settled down yet */
742  if(
743  (LoopEden==0 && fabs(CurrentEdenError)/dense.EdenTrue <= conv.EdenErrorAllowed/2.) ||
744  (LoopEden>0 && fabs(CurrentEdenError)/dense.EdenTrue <= conv.EdenErrorAllowed
745  && FractionChangeEden < conv.EdenErrorAllowed / 2.) )
746  break;
747 
748  /* diference is assumed and true electron density is larger
749  * than tolerance, we are not yet converged, default is 0.01 */
750  conv.lgConvEden = false;
751  strcpy( conv.chConvIoniz, "Ne big chg" );
752 
753  /* SEARCH is true if this is only search for first solution */
754  if( conv.lgSearch )
755  {
756  dense.eden = pow(10.,
757  (log10(dense.eden) + log10(dense.EdenTrue))/ 2.);
758  }
759  else
760  {
761 
762  /* was the sign of the change in the electron density changing,
763  * if so then we have bracked the target */
764  if( PreviousEdenError*CurrentEdenError < 0. )
765  {
766  /* we have bracketed the correct electron density,
767  * make changes smaller, also solve linear equation for zero error */
768  slope = (PreviousEdenError - CurrentEdenError ) /
769  ( eden_assumed_old - CurrentEden );
770 
771  ProposedEden = eden_assumed_old - PreviousEdenError / slope;
772  /* as sanity check, this must be between the two limits we examined,
773  * since zero was between them */
774  ASSERT( ProposedEden >= MIN2( eden_assumed_old , CurrentEden ) );
775  ASSERT( ProposedEden <= MAX2( eden_assumed_old , CurrentEden ) );
776 
777  change_eden_rel2_tolerance *= 0.25;
778  if( trace.nTrConvg>=3 )
779  fprintf( ioQQQ,
780  " ConvEdenIoniz3 bracketed electron density factor now %.3e error is %.4f proposed ne %.4e\n",
781  change_eden_rel2_tolerance,
782  (dense.EdenTrue-dense.eden)/dense.EdenTrue, ProposedEden);
783  }
784  else
785  {
786  /* did not bracket the target, keep moving in its direction */
787 
788  /* the correct electron density - is it too different? */
789  EdenUpperLimit = dense.eden * (1. + change_eden_rel2_tolerance);
790  EdenLowerLimit = dense.eden / (1. + change_eden_rel2_tolerance);
791 
792  /* get the new electron density */
793  if( dense.EdenTrue > EdenUpperLimit )
794  {
795  /* this branch, proposed eden too big */
796  ProposedEden = EdenUpperLimit;
797  }
798  else if( dense.EdenTrue < EdenLowerLimit )
799  {
800  /* proposed eden too small */
801  ProposedEden = EdenLowerLimit;
802  }
803  else
804  {
805  /* eden was within fac of the correct value, use it */
806  ProposedEden = dense.EdenTrue;
807  }
808  if( trace.nTrConvg>=3 )
809  fprintf( ioQQQ,
810  " ConvEdenIoniz3 elec dens fac %.3e err is %.4f prop ne %.4e cor ne %.4e slope=%.2e\n",
811  change_eden_rel2_tolerance,
813  ProposedEden,
814  dense.EdenTrue,slope);
815  }
816 
817  /* now remember this change for the next time through the loop */
818  PreviousEdenError = CurrentEdenError;
819  eden_assumed_old = CurrentEden;
820  FractionChangeEden = fabs( dense.eden - ProposedEden ) / dense.EdenTrue;
821  dense.eden = ProposedEden;
822 
823  if( trace.lgTrace && trace.lgNeBug )
824  {
825  fprintf( ioQQQ,
826  " ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e\n",
827  dense.eden, eden_assumed_old, dense.eden/eden_assumed_old, dense.EdenTrue );
828  }
829  }
830  /* loop until we give up -- normal exit is break in center of loop */
831  ++LoopEden;
832  } while( LoopEden < LoopLimit && !lgAbort );
833 
834  /* we now have the proposed new electron density */
835  /* we have changed the density slightly, so now must reevaluate ionization with this new value */
836  /* converge the current level of ioinization
837  * but only do this if change was significant */
838  if( trace.nTrConvg>=3 )
839  fprintf( ioQQQ,
840  " ConvEdenIoniz3 converged eden, current error is %.4f, calling ConvIoniz with final density=EdenTrue=%.4e\n",
843  /* , this calls eden_sum which updates EdenTrue */
844  /* >>chng 02 may 31, always call ConvIoniz (basically did so anyway) and remove eden_sum from here*/
845  if( ConvIoniz() )
846  {
847  if( trace.nTrConvg>=3 )
848  fprintf( ioQQQ,
849  " ConvEdenIoniz3 eden no longer converged eden, current error is %.4f\n",
851  }
852 
853  /* >>chng 01 mar 14, check whether electron density is no longer converged
854  * after reevaluating ionization */
855  if(
858  {
859  /* diference is assumed and true electron density is larger
860  * than tolerance, we are not yet converged, default is 0.01 */
861  conv.lgConvEden = false;
862  if( trace.nTrConvg>=3 )
863  fprintf( ioQQQ,
864  " ConvEdenIoniz3 setting eden not converged, error %.4f, exiting\n",
866  }
867  else
868  {
869  conv.lgConvEden = true;
870  }
871 
872  if( trace.lgTrace && trace.lgNeBug )
873  {
874  fprintf( ioQQQ,
875  " ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e converge=%c\n",
876  dense.eden, eden_assumed_old, dense.eden/eden_assumed_old, dense.EdenTrue ,TorF( conv.lgConvEden));
877  }
878 
879  if( trace.nTrConvg>=3 )
880  {
881  fprintf( ioQQQ, " ConvEdenIoniz3 %4ld new ne=%12.4e true=%12.4e",
882  LoopEden, dense.eden , dense.EdenTrue );
883 
884  /* this is flag says whether or not eden/eden has converged */
885  if( conv.lgConvEden )
886  {
887  fprintf( ioQQQ, " converged, eden rel error is %g ",
889  }
890  else
891  {
892  fprintf( ioQQQ, " NOCONV corr:%7.3f prop:%7.3f ",
893  (dense.EdenTrue-eden_assumed_old)/eden_assumed_old ,
894  (dense.eden-eden_assumed_old)/eden_assumed_old );
895  }
896  if( conv.lgEdenOscl )
897  fprintf( ioQQQ, " EDEN OSCIL2 \n");
898  }
899 
900  }
901  else
902  {
903  fprintf( ioQQQ, "ConvEdenIoniz finds insane solver %s \n" , conv.chSolverEden );
904  ShowMe();
905  }
906  ++LoopBig;
907  }
908 
909  if( trace.nTrConvg>=3 )
910  {
911  fprintf( ioQQQ, " ConvEdenIoniz3 exits, lgConvEden=%1c entry eden %.4e -> %.4e rel change %.3f\n" ,
912  TorF(conv.lgConvEden ),
913  EdenEntry ,
914  dense.eden ,
915  (EdenEntry-dense.eden)/SDIV( dense.eden ) );
916  }
917  else if( trace.lgTrace )
918  {
919  fprintf( ioQQQ, " ConvEdenIoniz exits zn %.2f lgConvEden=%1c entry eden %.4e -> %.4e rel change %.3f\n" ,
920  fnzone,
921  TorF(conv.lgConvEden ),
922  EdenEntry ,
923  dense.eden ,
924  (EdenEntry-dense.eden)/SDIV( dense.eden ) );
925  }
926 
927  return 0;
928 }
929 
930 /* returns true if electron density is converged */
931 bool lgConvEden(void)
932 {
933  bool lgRet;
934 
935  /* >>chng 04 sep 26, div by eden, not min of eden and edentrue since latter now poss neg */
936  if( fabs(dense.eden-dense.EdenTrue)/SDIV(dense.eden) >=
938  {
939  lgRet = false;
940  conv.lgConvEden = false;
941  }
942  else
943  {
944  lgRet = true;
945  conv.lgConvEden = true;
946  }
947  return lgRet;
948 }

Generated for cloudy by doxygen 1.8.1.1