cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
conv_temp_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 /*ConvTempEdenIoniz determine temperature, called by ConPresTempEdenIoniz,
4  * calls ConvEdenIoniz to get electron density and ionization */
5 /*MakeDeriv derive numerical derivative of heating minus cooling */
6 /*PutHetCol save heating, cooling, and temperature in stack for numerical derivatives */
7 /*lgConvTemp returns true if heating-cooling is converged */
8 /*CoolHeatError evaluate ionization, and difference in heating and cooling, for temperature temp */
9 #include "cddefines.h"
10 #include "hmi.h"
11 #include "thermal.h"
12 #include "iso.h"
13 #include "hydrogenic.h"
14 #include "colden.h"
15 #include "h2.h"
16 #include "pressure.h"
17 #include "dense.h"
18 #include "trace.h"
19 #include "phycon.h"
20 #include "conv.h"
21 
22 /* >>chng 01 apr 06, from 10 to 20, also put damper at 0.5 below,
23  * as part of strategy of letting code bracket solution */
24 /* >>chng 04 mar 17, from 20 to 40, now doing much finer changes in
25  * many conditions, including smaller dT, so allow more loops */
26 /* >>chng 05 mar 25, from 40 to 60, thermal fronts */
27 static const int LIMDEF = 60;
28 
29 /*CoolHeatError evaluate ionization, and difference in heating and cooling, for temperature temp */
30 STATIC double CoolHeatError( double temp );
31 
32 /* use brent's method to find heating-cooling match */
33 STATIC double TeBrent(
34  double x1,
35  double x2);
36 
37 /*MakeDeriv derive numerical derivative of heating minus cooling */
38 STATIC void MakeDeriv(
39  /* which job to do, either "zero" or "incr" */
40  const char *job,
41  /* result, the numerical derivative */
42  double *DerivNumer);
43 
44 /*PutHetCol save heating, cooling, and temperature in stack for numerical derivatives */
45 STATIC void PutHetCol(double te,
46  double htot,
47  double ctot);
48 
49 /*ConvTempEdenIoniz determine temperature, called by ConPresTempEdenIoniz,
50  * calls ConvEdenIoniz to get electron density and ionization
51  * returns 0 if ok, 1 if disaster */
52 int ConvTempEdenIoniz( void )
53 {
54  char chDEType;
55  long int limtry,
56  nTemperature_loop;
57  /* following is zero because there are paths were not set following abort */
58  double CoolMHeat=0.,
59  Damper,
60  dtl,
61  fneut;
62  int kase;
63  double DerivNumer;
64  static double OldCmHdT = 0.;
65  long int nTemp_oscil;
66 
67  /* derivative of cooling less heating wrt temperature */
68  double dCoolmHeatdT;
69 
70  DEBUG_ENTRY( "ConvTempEdenIoniz()" );
71 
72  /* determine ionization and temperature, called by PIonte
73  * nCall tells routine which call this is, 0 for first one, not changed here
74  * nTemperature_loop counts how many loops performed */
75 
76  if( trace.lgTrace )
77  {
78  fprintf( ioQQQ, "\n ConvTempEdenIoniz called\n" );
79  }
80  if( trace.nTrConvg>=2 )
81  {
82  fprintf( ioQQQ, " ConvTempEdenIoniz1 called. Te:%.4e Htot:%11.4e Ctot:%11.4e error:%9.1e%%, eden=%11.4e\n",
84  100./MAX2(1e-36,thermal.htot), dense.eden );
85  }
86 
87  if( strcmp( conv.chSolverTemp , "new" ) == 0 )
88  {
89  /* this has never worked */
90  double t1 , t2 ,
91  error1 , error2;
92  double TeNew;
93  double factor = 1.02;
94 
95  t1 = phycon.te;
96  error1 = CoolHeatError( t1 );
97  t2 = t1 * factor;
98  error2 = CoolHeatError( t2 );
99  TeNew = TeBrent( t1 , t2 );
100  fprintf(ioQQQ , "DEBUG new temp solver error1 %.2e error2 %.2e TeNew %.2e\n",
101  error1 , error2 , TeNew );
102  }
103 
104  else if( strcmp( conv.chSolverTemp , "simple") == 0 )
105  {
106 
107  /* main routine to find ionization and temperature
108  * following is flag to check for temp oscillations
109  * it could be set true in main temp loop */
110  conv.lgTOscl = false;
111  /* ots rates not oscillating */
112  conv.lgOscilOTS = false;
113 
114  /* scale factor for changes in temp - make smaller if oscillations */
115  Damper = 1.;
116 
117  /* flag for d(cool - heat)/dT changing sign, an internal error
118  * >>chng 97 sep 03, only turn false one time, remains true if ever true */
119  if( nzone < 1 )
120  {
121  conv.lgCmHOsc = false;
122  }
123 
124  /* option to make numerical derivatives of heating cooling */
125  MakeDeriv("zero",&DerivNumer);
126 
127  /* this is will initial limit to number of tries at temp */
128  limtry = LIMDEF;
129 
130  /* used to keep track of sign of dt, to check for oscillations */
131  dtl = 0.;
132 
133  /* this counts how may thermal loops we go through
134  * used in calling routine to make sure that at least
135  * two solutions are performed */
136  nTemperature_loop = 0;
137  nTemp_oscil = 0;
138 
139  /* this is change in temp, which is used early in following loop, must be preset */
140  thermal.dTemper = 1e-3f*phycon.te;
141 
142  /* this is the start of the heating-cooling match loop,
143  * this exits by returning in the middle */
144  while ( true )
145  {
146 
147  /* >>chng 02 dec 03 rjrw, insert these so dynamics cooling rate is calculated
148  * before step is chosen -- but how much duplication results?
149  * must have current estimate of heating and cooling before temperature is changed.
150  * at bottom of following loop these two are called again to see whether heatin cooling
151  * have converged */
152  if( ConvEdenIoniz() )
153  {
154  return 1;
155  }
156  PresTotCurrent();
157 
158  /* we will try to make the following zero */
159  CoolMHeat = thermal.ctot - thermal.htot;
160  /*fprintf(ioQQQ,"DEBUG grn 13 14\t%e\t%e\t%e\n",
161  phycon.te,
162  thermal.heating[0][13],
163  thermal.heating[0][14] );*/
164 
166  {
167  /* constant temp - declare this a success */
168  CoolMHeat = 0.;
169  }
170 
171  if( thermal.lgTLaw )
172  {
173  double TeNew = phycon.te;
174  /* temperature law has been specified */
175  if( thermal.lgTeBD96 )
176  {
177  /* Bertoldi & Drain 96 temp law specified by TLAW BD96 command */
178  TeNew = thermal.T0BD96 / (1.f + thermal.SigmaBD96 * colden.colden[ipCOL_HTOT] );
179  }
180  else if( thermal.lgTeSN99 )
181  {
182  /* Sternberg & Neufeld 99 temp law specified by TLAW SN99 command,
183  * this is equation 16 of
184  * >>refer H2 temperature law Sternberg, A., & Neufeld, D.A. 1999, ApJ, 516, 371-380 */
185  TeNew = thermal.T0SN99 /
186  /*(realnum)(1.f + 9.*pow(2.*hmi.Hmolec[ipMH2g]/dense.gas_phase[ipHYDROGEN], 4.) );*/
187  (realnum)(1.f + 9.*pow(2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN], 4.) );
188  }
189  else
190  TotalInsanity();
191 
192  TempChange(TeNew ,false);
193  }
194 
195  /* >>chng 02 dec 09, this block moved to just after loop entry above
196  * to check on heating cooling and
197  * exit if match is ok - prevents unnecessary reevaluate of ConvEdenIoniz and PresTotCurrent */
198  /* this is test for valid thermal solution,
199  * lgConvTemp returns true if difference in heat cool is small */
200  /* >>chng 02 dec 04, add check on size of dTemper rel to te,
201  * can become too small to increment a float */
202  /* >>chng 04 dec 13, for very cold gas the dTemper < 1e-6 Te limit was hit
203  * during normal search. do not let this limit be tested here, rather use
204  * fltepsilon to set limit to dTemper */
205  if( (lgConvTemp() /*|| fabs(thermal.dTemper) < 1e-6*phycon.te*/)
206  && conv.lgConvIoniz )
207  {
208  /* we have a good solution */
209  if( trace.lgTrace )
210  fprintf( ioQQQ, " ConvTempEdenIoniz returns ok.\n" );
211 
212  /* test for thermal stability, set flag if unstable,
213  * later increment number of zones only once per zone */
214  /* this may become a check on thermal stability since neg slope is
215  * unstable, but we may not have a valid solution, so only set flag */
216  if( thermal.dCooldT < 0. )
217  {
218  thermal.lgUnstable = true;
219  }
220  else
221  {
222  thermal.lgUnstable = false;
223  }
224 
225  /* remember the highest and lowest temperatures */
228 
229  /* say that we have a valid solution */
230  conv.lgConvTemp = true;
231 
232  if( trace.nTrConvg>=2 )
233  {
234  fprintf( ioQQQ, " ConvTempEdenIoniz2 converged. Te:%11.4e Htot:%11.4e Ctot:%11.4e error:%9.1e%%, eden=%11.4e\n",
236  100./MAX2(1e-36,thermal.htot), dense.eden );
237  }
238  /*******************************************************
239  *
240  * this is return for valid solution
241  *
242  *******************************************************/
243  return 0;
244  }
245  else if(nTemperature_loop >= limtry || lgAbort )
246  {
247  /* >>chng 02 dec 17, add test on ots fields oscillating
248  * while ionization is not converged */
249  /* >>chng 04 may 25, do not test on lgOscilOTS since we do test
250  * on lgConvIoniz */
251  if( nTemperature_loop == LIMDEF && !conv.lgTOscl &&
252  conv.lgConvIoniz /*&& !conv.lgOscilOTS*/ )
253  {
254  /* LIMDEF is set to 60 above
255  * >>chng 97 aug 10, from 2 to 4 to let keep going when not oscillation
256  * this helped in getting over huge jumps in temperature at metal fronts */
257  limtry = LIMDEF*4;
258  if( trace.nTrConvg>2 )
259  {
260  fprintf( ioQQQ, " ConvTempEdenIoniz9 increases limtry to%3ld\n",
261  limtry );
262  }
263  }
264  else
265  {
266  /* looped too long, so bail out anyway */
267  break;
268  }
269  }
270 
271  /* get numerical heating-cooling derivative, but we may not use it */
272  MakeDeriv("incr",&DerivNumer);
273 
274  /* remember this temp, heating, and cooling */
276 
277  /* the flag conv.lgConvEden was set during above call to ConvEdenIoniz, and
278  * if is not true then eden failed, and we abort,
279  * but don't abort on very first try, when flag has not been set
280  * >> chng 02 dec 10 rjrw flag is now set on first try */
281  if( !conv.lgConvEden )
282  {
283  /* set this to zero so that will not flag temperature failure, since we are
284  * aborting due to eden failure */
285  CoolMHeat = 0.;
286  if( trace.nTrConvg>=2 )
287  {
288  fprintf( ioQQQ,
289  " ConvTempEdenIoniz3 breaks since lgConvEden returns failed electron density.\n");
290  }
291  break;
292  }
293 
294  /* save old derivative to check for oscillations */
295  if( nzone < 1 )
296  {
297  /* >>chng 96 jun 07, was set to dCoolmHeatdT, is not defined initially */
298  OldCmHdT = thermal.ctot/phycon.te;
299  }
300 
301 # define USENUMER false
302  if( USENUMER )
303  {
304  dCoolmHeatdT = DerivNumer;
305  chDEType = 'n';
306  }
307  else if( phycon.te < 1e4 && thermal.dCooldT < 0. )
308  {
309  /* use numerical guess, HTOT nearly equal to CTOT, drive Te lower
310  * this is to overcome thermal inflection points */
311  dCoolmHeatdT = (double)(thermal.htot)/phycon.te*2.;
312  chDEType = 'd';
313  }
314  else
315  {
316  /* >>chng 97 mar 18, all below growth code just did the following
317  * use analytical estimate of change in heat - cooling wrt temp */
318  dCoolmHeatdT = thermal.dCooldT - thermal.dHeatdT;
319  chDEType = 'a';
320  }
321 
322  /* check whether derivative changing sign - an internal error or
323  * numerical instability */
324  if( OldCmHdT*dCoolmHeatdT < 0. )
325  {
326  conv.lgCmHOsc = true;
327 
328  /* derivative can change sign when heating and cooling balance,
329  * and processes have exactly same temp dependence
330  * happened in mods of ism.in test
331  * >>chng 97 sep 02, added following test for oscillating derivs */
332  dCoolmHeatdT = thermal.ctot/phycon.te;
333  chDEType = 'o';
334  }
335  OldCmHdT = dCoolmHeatdT;
336 
337  /* if not first pass through then numerical derivative should
338  * not be zero required change in temperature, this may be too large */
339  thermal.dTemper = (realnum)(CoolMHeat/dCoolmHeatdT);
340  kase = 0;
341 
342  /* try to stop numerical oscillation if dTemper is changing sign */
343  /* when first time through and far from solution (when zones a bit too thick) first
344  * step can be in wrong direction, and all subsequent are in right direction.
345  * do not call this an oscillation - only declare oscillation when second occurs,
346  * nTemperature_loop is zero first time through, (dtl is also zero),
347  * is 1 first time through with memory of this zone, (when false oscillation occurs)
348  * and is >1 thereafter */
349  /* >>chng 05 mar 25, had counted as oscillations when loop count greater than 1, change to 3
350  * to allow things to settle down. Also, for first three evaluations, do not let te
351  * change by too much, until lower solvers settle down */
352  if( dtl*thermal.dTemper < 0. && nTemperature_loop > 3)
353  {
354  conv.lgTOscl = true;
355  /* make DT smaller and smaller */
356  /* >>chng 01 mar 16, from 0.5 to 0.75
357  Damper *= 0.75;*/
358  Damper *= 0.5;
359  nTemp_oscil = nTemperature_loop;
360  }
361  else if( Damper < 1. && nTemperature_loop-nTemp_oscil > 10 )
362  {
363  /* >>chng 05 mar 26, some cases there is noise in heating-cooling
364  * in first few evaluations and this causes te oscillations. but it may
365  * settle down and become well behaved. this allows that to happen */
366  conv.lgTOscl = false;
367  if( !((nTemperature_loop-nTemp_oscil)%5) )
368  Damper = MIN2( 1., Damper*1.5 );
369  }
370 
371  /* SIGN: sign of sec arg and abs val of first */
372  /* >>chng 96 nov 08, from 0.1 to 0.08 in fneut */
373  /*fneut = (dense.xIonDense[ipHYDROGEN][0] + 2.*hmi.Hmolec[ipMH2g])/dense.gas_phase[ipHYDROGEN];*/
375  /*>>chng 04 jan 11, when big H2 is on and heating is important,
376  * only allow small changes in temp since must converge pops */
377  /* >>chng 04 jan 15, chng ratio of heating in h2 to total as the criteria,
378  * to the error. temp changes will become much larger when this is
379  * not tripped, and h2 can start to oscillate wildly.
380  * NB this limit should be kept parallel with similar test for convergence
381  * of heating in h2 */
383  {
384  /* heating due to collisional deexcitation of the large H2 molecule
385  * is very important. this rate depends strongly on the temperature
386  * due to Boltzmann factors, and large changes in temperature
387  * inhibit convergence of the level pops in the big H2. That
388  * atom has tests on converging the heating, so want only small
389  * temperature changes when H2 heating is important. */
390  double absdt = fabs(thermal.dTemper);
391  thermal.dTemper = (sign(MIN2(absdt,0.0025*phycon.te),
392  thermal.dTemper));
393  kase = 1;
394  }
395  else if( fneut > 0.08 && hydro.lgHColionImp )
396  {
397  /* lgHColionImp is true if model collisionally ionized
398  * do not let TE change by too much if lgHColionImp is set by HLEVEL */
399  double absdt = fabs(thermal.dTemper);
400  thermal.dTemper = (sign(MIN2(absdt,0.005*phycon.te),
401  thermal.dTemper));
402  kase = 2;
403  }
404 # if 0
405  /* >>chng 05 mar 25, for first few steps through solvers in a new zone, don't let temp
406  * change by too much initially, to avoid a te oscil which sets the oscil flag */
407  else if( 0 && nTemperature_loop <= 3)
408  {
409  /* >>chng 05 mar 25, add this branch, an early pass through these solvers,
410  * do not let te change by too much to avoid te oscillations */
411  double absdt = fabs(thermal.dTemper);
412  thermal.dTemper = (sign(MIN2(absdt,0.002*phycon.te),
413  thermal.dTemper));
414  kase = 3;
415  }
416 # endif
417  else
418  {
419  /* this branch, we will use the change in temperature that came from the deriv,
420  * but do not let temp change by more than 2 percent */
421  double absdt = fabs(thermal.dTemper);
422  thermal.dTemper = (sign(MIN2(absdt,0.02*phycon.te),
423  thermal.dTemper));
424  kase = 4;
425  }
426 
427  /* >>chng 04 mar 14, add this test, cooling flow clouds had temp failure
428  * when dt as large as 4% were allowed with the big H2 */
429  /* in no case do we want the temperature to change by more than 1% when the big
430  * H2 molecule is turned on - it needs small temp changes to converge */
431  if( 0 && h2.lgH2ON && fabs(thermal.dTemper)/ phycon.te > 0.01 )
432  {
433  double absdt = fabs(thermal.dTemper);
434  thermal.dTemper = (sign(MIN2(absdt,0.01*phycon.te),
435  thermal.dTemper));
436  kase = 5;
437  }
438 
439  /* >>chng 03 mar 16, to following logic, previous had been impossible
440  * to hit - only do fine changes in Te when at half-way point in i-front
441  * or if model is in collisional recombination equilibrium
442  * if so then small change in temp causes big change in equilibrium
443  * lgHColRec set in hlevel, only true is totl rec >> rad rec */
444  /* are most recombinations due to 3-body? */
446  /* is hydrogen more than 30% ionized */
448  {
449  if( iso.RecomCollisFrac[ipH_LIKE][ipHYDROGEN] > 0.8 )
450  /* >>chng 04 feb 07, add test for H mostly collisionally ionized */
451  /* is hydrogen less than 70% ionized? */
452  /* >>chng 04 feb 07, from 70% to 90% ionized - oscill in hard
453  * x-ray continua (laser at 80 ryd) when i front hit */
454  {
455  double absdt = fabs(thermal.dTemper);
456  /* >>chng 03 nov 03 , fac from 0.001 to 0.004 */
457  /* >>chng 04 feb 24 , fac back from 0.004 to 0.001,
458  * this is needed for col den 25, high density (14) blr models
459  * with flux of 22 - small changes in temperature drive large
460  * changes in electron density due to three body recomb
461  * this small change is ok due to test above, which makes sure
462  * that we are in collisional recombination limit
463  * sign returns sign of sec par times value of first */
464  thermal.dTemper = (sign(MIN2(absdt,0.001*
466  kase = 6;
467  }
468  }
469 
470  /* make steps smaller if oscillations, caused by overshots, occur */
471  thermal.dTemper *= Damper;
472 
473  /* do not let dTemper/Te grow smaller than 10 * flt epsilon, the smaller number that
474  * can be added to a float and still work */
475  if( fabs(thermal.dTemper)/phycon.te <10.*DBL_EPSILON )
476  {
477  kase = 7;
478  thermal.dTemper = (sign( 10.*DBL_EPSILON*phycon.te , thermal.dTemper));
479  }
480 
481  /* save last change in temperature */
482  dtl = thermal.dTemper;
483 
484  /* THIS IS IT !!! this is it !!! this is where te changes. */
485  if( !thermal.lgTLaw )
486  {
487  /* usual case - valid thermal solution */
488  TempChange(phycon.te - thermal.dTemper , false);
489  }
490 
491  if( trace.nTrConvg>=2 )
492  {
493  fprintf( ioQQQ,
494  " ConvTempEdenIoniz4 a %4li T:%.4e dt:%10.2e%7.2f%% C:%10.2e H:%10.2e"
495  " CoolMHeat/h:%7.2f%% dC-H/dT:%.2e kase:%i%c damp%.2f\n",
496  nTemperature_loop,
497  phycon.te,
498  thermal.dTemper,
500  thermal.ctot,
501  thermal.htot,
502  CoolMHeat/MIN2(thermal.ctot , thermal.htot )*100. ,
503  dCoolmHeatdT,
504  kase,
505  /* which type of deriv we have */
506  chDEType,
507  Damper);
508  /* dCoolmHeatdT is the actual number used for getting dT
509  * thermal.dCooldT is just cooling */
510  }
511 
512  if( trace.lgTrace )
513  {
514  fprintf( ioQQQ,
515  " ConvTempEdenIoniz TE:%.5e dT%.3e HTOT:%.5e CTOT:%.5e dCoolmHeatdT:%c%.4e C-H%.4e HDEN%.2e kase %i\n",
516  phycon.te,
517  thermal.dTemper,
518  thermal.htot,
519  thermal.ctot,
520  chDEType,
521  dCoolmHeatdT,
522  CoolMHeat,
524  kase );
525  }
526 
527  /* increment loop counter */
528  ++nTemperature_loop;
529  }
530 
531  /* fall through, exceeded limit to number of solutions,
532  * >>chng 01 mar 14, or no electron density convergence
533  * no temperature convergence */
534  if( fabs(CoolMHeat/thermal.htot ) > conv.HeatCoolRelErrorAllowed )
535  {
536  conv.lgConvTemp = false;
537  if( trace.nTrConvg>=2 )
538  {
539  fprintf( ioQQQ, " ConvTempEdenIoniz7 fell through loop, heating cooling not converged, setting conv.lgConvTemp = false\n");
540  }
541  }
542  else
543  {
544  /* possible that ionization has failed, and we fell through to here
545  * with a valid thermal solution */
546  conv.lgConvTemp = true;
547  if( trace.nTrConvg>=2 )
548  {
549  fprintf( ioQQQ, " ConvTempEdenIoniz8 fell through loop, heating cooling is converged (ion not?), setting conv.lgConvTemp = true\n");
550  }
551  }
552  }
553  else
554  {
555  fprintf( ioQQQ, "ConvTempEdenIoniz finds insane solver%s \n" , conv.chSolverTemp );
556  ShowMe();
557  cdEXIT(EXIT_FAILURE);
558  }
559  return 0;
560 }
561 
562 /*MakeDeriv derive numerical derivative of heating minus cooling */
564  const char *job,
565  double *DerivNumer)
566 {
567  static long int nstore;
568  static double OldTe , OldCool , OldHeat;
569 
570  DEBUG_ENTRY( "MakeDeriv()" );
571 
572  if( strcmp(job,"zero") == 0 )
573  {
574  nstore = 0;
575  OldTe = 0.;
576  OldCool = 0.;
577  OldHeat = 0.;
578  }
579 
580  else if( strcmp(job,"incr") == 0 )
581  {
582 
583  if( nstore > 0 && !thermal.lgTemperatureConstant && fabs(phycon.te - OldTe)>SMALLFLOAT )
584  {
585  /* get numerical derivatives */
586  double dCdT = (thermal.ctot - OldCool)/(phycon.te - OldTe);
587  double dHdT = (thermal.htot - OldHeat)/(phycon.te - OldTe);
588  *DerivNumer = dCdT - dHdT;
589 # if 0
590  fprintf(ioQQQ,"derivderiv\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
591  nzone, dCdT , thermal.dCooldT ,
592  thermal.ctot , OldCool,phycon.te , OldTe
593  /*(thermal.ctot - OldCool),(phycon.te - OldTe),dHdT, thermal.dHeatdT*/);
594 # endif
595  }
596 
597  else
598  {
599  /* if early in soln, or constant temperature, return zero */
600  *DerivNumer = 0.;
601  }
602 
603  OldTe = phycon.te;
604  OldCool = thermal.ctot;
605  OldHeat = thermal.htot;
606  ++ nstore;
607  }
608 
609  else
610  {
611  fprintf( ioQQQ, " MakeDeriv called with insane job=%4.4s\n",
612  job );
613  cdEXIT(EXIT_FAILURE);
614  }
615  return;
616 }
617 
618 /*PutHetCol save heating, cooling, and temperature in stack for numerical derivatives */
619 STATIC void PutHetCol(double te,
620  double htot,
621  double ctot)
622 {
623 
624  DEBUG_ENTRY( "PutHetCol()" );
625 
626  if( nzone == 0 )
627  {
628  /* code is searching for first temp now - do not remember these numbers */
629  thermal.ipGrid = 0;
630  }
631  else
632  {
633  if( thermal.ipGrid >= NGRID )
634  thermal.ipGrid = 0;
635 
636  ASSERT( thermal.ipGrid >= 0 );
637  ASSERT( thermal.ipGrid < NGRID );
638  ASSERT( nzone>=0 );
639 
644 
645  ++thermal.ipGrid;
646  }
647  return;
648 }
649 
650 /*CoolHeatError evaluate ionization, and difference in heating and cooling, for temperature temp */
651 STATIC double CoolHeatError( double temp )
652 {
653  DEBUG_ENTRY( "CoolHeatError()" );
654 
655  TempChange(temp , false);
656 
657  /* converge the ionization and electron density;
658  * this calls ionize until lgIonDone is true */
659  /* NB should NOT set insanity - but rather return error condition */
660  if( ConvEdenIoniz() )
661  lgAbort = true;
662 
663  /* >>chng 01 mar 16, evaluate pressure here since changing and other values needed */
664  /* reevaluate pressure */
665  /* this sets values of pressure.PresTotlCurr */
666  PresTotCurrent();
667 
668  /* remember this temp, heating, and cooling */
670  return thermal.ctot - thermal.htot;
671 }
672 
673 /*#define EPS 3.e-8*/
674 #define ITERMAX 100
675 /* use brent's method to find heating-cooling match */
677  double x1,
678  double x2)
679 {
680  long int iter;
681  double c,
682  d,
683  e,
684  fx1,
685  fx2,
686  fc,
687  p,
688  q,
689  r,
690  s,
691  xm;
692 
693  DEBUG_ENTRY( "TeBrent()" );
694 
695  /* first evaluate function at two boundaries, and check
696  * that we have bracketed the target */
697  /*NB big logical error - this routine ignores return value of
698  * routine it calls, so will continue despite disaster */
699  fx1 = CoolHeatError(x1);
700  fx2 = CoolHeatError(x2);
701 
702  if( fx1*fx2 > 0. )
703  {
704  /* both values either positive or negative, this is insane */
705  fprintf( ioQQQ, " TeBrent called without proper bracket - this is impossible\n" );
706  fprintf( ioQQQ, " Sorry.\n" );
707  cdEXIT(EXIT_FAILURE);
708  }
709 
710  /* we have bracketed the target */
711  c = x2;
712  fc = fx2;
713  iter = 0;
714  /* these are sanity checks, to get code past lint */
715  e = DBL_MAX;
716  d = DBL_MAX;
717  while( iter < ITERMAX )
718  {
719  if( (fx2 > 0. && fc > 0.) || (fx2 < 0. && fc < 0.) )
720  {
721  c = x1;
722  fc = fx1;
723  d = x2 - x1;
724  e = d;
725  }
726  if( fabs(fc) < fabs(fx2) )
727  {
728  x1 = x2;
729  x2 = c;
730  c = x1;
731  fx1 = fx2;
732  fx2 = fc;
733  fc = fx1;
734  }
735 
736  xm = .5*(c - x2);
737 
738  /* solution converged if residual less than tolerance, or hit zero */
739  if( fabs(xm) <= thermal.htot*conv.HeatCoolRelErrorAllowed || fx2 == 0. )
740  break;
741 
742  if( fabs(e) >= thermal.htot*conv.HeatCoolRelErrorAllowed && fabs(fx1) > fabs(fx2) )
743  {
744  double aa , bb;
745  s = fx2/fx1;
746  if( fp_equal( x1, c ) )
747  {
748  p = 2.*xm*s;
749  q = 1. - s;
750  }
751  else
752  {
753  q = fx1/fc;
754  r = fx2/fc;
755  p = s*(2.*xm*q*(q - r) - (x2 - x1)*(r - 1.));
756  q = (q - 1.)*(r - 1.)*(s - 1.);
757  }
758  if( p > 0. )
759  q = -q;
760 
761  p = fabs(p);
763  bb = fabs(e*q);
764  if( 2.*p < MIN2(3.*xm*q-aa,bb) )
765  {
766  e = d;
767  d = p/q;
768  }
769  else
770  {
771  d = xm;
772  e = d;
773  }
774  }
775  else
776  {
777  d = xm;
778  e = d;
779  }
780  x1 = x2;
781  fx1 = fx2;
782  if( fabs(d) > thermal.htot*conv.HeatCoolRelErrorAllowed )
783  {
784  x2 += d;
785  }
786  else
787  {
789  }
790  fx2 = CoolHeatError(x2);
791 
792  ++iter;
793  }
794 
795  /* check whether we fell through (hit limit to number of iterations)
796  * of broke out when complete */
797  if( iter == ITERMAX )
798  {
799  /* both values either positive or negative, this is insane */
800  fprintf( ioQQQ, " TeBrent did not converge in %i iterations\n",ITERMAX );
801  fprintf( ioQQQ, " Sorry.\n" );
802  cdEXIT(EXIT_FAILURE);
803  }
804  return( x2 );
805 
806 }
807 
808 /* returns true if heating-cooling is converged */
809 bool lgConvTemp(void)
810 {
811  bool lgRet;
812  DEBUG_ENTRY( "lgConvTemp()" );
813 
814  /* >>chng 03 sep 19, add check on phycon.TEMP_LIMIT_LOW */
817  {
818  /* announce that temp is converged if relative heating - cooling mismatch
819  * is less than the relative heating cooling error allowed, or
820  * if this is a constant temperature model */
821  lgRet = true;
822  conv.lgConvTemp = true;
823  }
824  else
825  {
826  /* big mismatch, this has not converged */
827  lgRet = false;
828  conv.lgConvTemp = false;
829  }
830  return lgRet;
831 }

Generated for cloudy by doxygen 1.8.4