cloudy  trunk
 All Data Structures Namespaces 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-2017 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 /*lgConvTemp returns true if heating-cooling is converged */
6 /*CoolHeatError evaluate ionization, and difference in heating and cooling, for temperature temp */
7 /*DumpCoolStack helper routine to dump major coolants */
8 /*DumpHeatStack helper routine to dump major heating agents */
9 #include "cddefines.h"
10 #include "hmi.h"
11 #include "thermal.h"
12 #include "colden.h"
13 #include "pressure.h"
14 #include "dense.h"
15 #include "trace.h"
16 #include "phycon.h"
17 #include "conv.h"
18 #include "physconst.h"
19 #include "iter_track.h"
20 #include "radius.h"
21 
22 /*lgConvTemp returns true if heating-cooling is converged */
23 STATIC bool lgConvTemp(const iter_track& TeTrack);
24 /*CoolHeatError evaluate ionization, and difference in heating and cooling, for temperature temp */
25 STATIC double CoolHeatError( double temp );
26 
27 // debugging routines to print main sources of cooling and heating
28 STATIC void DumpCoolStack(double thres);
29 STATIC void DumpHeatStack(double thres);
30 
31 /*ConvTempEdenIoniz determine temperature, called by ConPresTempEdenIoniz,
32  * calls ConvEdenIoniz to get electron density and ionization
33  * returns 0 if ok, 1 if disaster */
34 int ConvTempEdenIoniz( void )
35 {
36  DEBUG_ENTRY( "ConvTempEdenIoniz()" );
37 
38  if( trace.lgTrace )
39  {
40  fprintf( ioQQQ, "\n ConvTempEdenIoniz called\n" );
41  }
42  if( trace.nTrConvg >= 2 )
43  {
44  fprintf( ioQQQ, " ConvTempEdenIoniz called, entering temp loop using solver %s.\n",
46  }
47 
48  // deal with special temperature laws first
50  {
51  if( thermal.lgTLaw )
52  {
53  double TeNew = phycon.te;
54  if( thermal.lgTeBD96 )
55  {
56  /* Bertoldi & Drain 96 temp law specified by TLAW BD96 command */
58  }
59  else if( thermal.lgTeSN99 )
60  {
61  /* Sternberg & Neufeld 99 temp law specified by TLAW SN99 command,
62  * this is equation 16 of
63  * >>refer H2 temp Sternberg, A., & Neufeld, D. A. 1999, ApJ, 516, 371-380 */
64  TeNew = thermal.T0SN99 /
65  (1. + 9.*POW4(2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN]) );
66  }
67  else if( thermal.lgTeTLaw )
68  {
69  /* Tabulate temperature tlaw */
71  }
72  else
73  TotalInsanity();
74 
75  TempChange( TeNew, false );
76  }
77 
78  if( ConvEdenIoniz() )
79  return 1;
81 
82  // convergence is automatic...
83  conv.lgConvTemp = true;
84 
85  if( trace.lgTrace || trace.nTrConvg >= 2 )
86  {
87  fprintf( ioQQQ, " ConvTempEdenIoniz: Te %e C %.4e H %.4e\n",
89  fprintf( ioQQQ, " ConvTempEdenIoniz returns ok.\n" );
90  }
91 
92  return 0;
93  }
94 
95  // this branch uses the van Wijngaarden-Dekker-Brent method
96  if( strcmp( conv.chSolverTemp , "vWDB" ) == 0 )
97  {
98  conv.lgConvTemp = false;
99 
100  // here starts the standard solver for variable temperature
101  iter_track TeTrack;
102  double t1=0, error1=0, t2, error2;
103 
104  t2 = phycon.te;
105  error2 = CoolHeatError( t2 );
106 
107  for( int n=0; n < 5 && !lgAbort; ++n )
108  {
109  const int DEF_ITER = 10;
110  const double DEF_FACTOR = 0.2;
111  double step, factor = DEF_FACTOR;
112 
113  TeTrack.clear();
114 
115  step = min( abs(safe_div( error2, conv.dCmHdT, 0. )), factor*t2 );
116 
117  // set up an initial guess for the bracket
118  // t2 was already initialized outside the main loop, or is copied from the
119  // previous iteration. don't record this evaluation, it may be poorly converged
120  for( int i=0; i < 100 && !lgAbort; ++i )
121  {
122  t1 = t2;
123  error1 = error2;
124 
125  double maxstep = factor*t1;
126  // limited testing on the auto test suite shows that sqrt(2)
127  // is close to the optimal value
128  step = SQRT2*step;
129  if( step == 0.0 || step > maxstep )
130  step = maxstep;
131  t2 = max( t1 + sign( step, -error1 ), phycon.TEMP_LIMIT_LOW );
132  error2 = CoolHeatError( t2 );
133  TeTrack.add( t2, error2 );
134 
135  // if n > 0, this indicates a previous failure to solve Te
136  // this could be due to hysteresis (e.g. O-H charge transfer)
137  // so ignore the first n steps, even if they seem to indicate
138  // that a bracket is found, to allow the code some time to settle
139  if( i >= n && error1*error2 <= 0. )
140  break;
141 
142  // test for i >= n here to give the code a chance to declare
143  // "bracket found" before aborting...
144  if( i >= n && fp_equal( t2, phycon.TEMP_LIMIT_LOW ) )
145  {
146  /* temp is too low */
147  fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature appears to be below the lower limit of the code,"
148  " %.3eK. It does not bracket thermal balance.\n",
150  fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n");
151  lgAbort = true;
152  return 1;
153  }
154  }
155 
156  if( trace.nTrConvg >= 2 && error1*error2 > 0. && !lgAbort )
157  {
158  fprintf( ioQQQ, " ConvTempEdenIoniz: bracket1 fails t1: %e %e t2: %e %e\n",
159  t1, error1, t2, error2 );
160  TeTrack.print_history();
161  }
162 
163  // keeping the history up until now has a bad effect on convergence
164  // so we wipe the slate clean....
165  TeTrack.clear();
166 
167  // the bracket should have been found, now set up the Brent solver
168  if( TeTrack.init_bracket( t1, error1, t2, error2 ) == 0 )
169  {
170  // The convergence criterion is based on the relative accuracy of Cool-Heat,
171  // combined with a relative accuracy on the temperature itself. We need to
172  // keep iterating until both accuracies are reached. Here we set tolerance on
173  // Te to 2 ulp. If bracket gets narrower than 3 ulp we declare a convergence
174  // failure to avoid changes getting lost in machine precision.
175  TeTrack.set_tol(2.*DBL_EPSILON*t2);
176 
177  if( error1 != 0.0 || error2 != 0.0 )
178  t2 = (t1*error2-t2*error1)/(error2-error1);
179  else
180  t2 = 0.5*(t1+t2);
181 
182  for( int i = 0; i < (1<<(n/2))*DEF_ITER && !lgAbort; i++ )
183  {
184  // check for convergence, as well as a pathologically narrow bracket
185  if( lgConvTemp(TeTrack) || TeTrack.bracket_width() < 3.*DBL_EPSILON*t2 )
186  break;
187 
188  error2 = CoolHeatError( t2 );
189  TeTrack.add( t2, error2 );
190  t2 = TeTrack.next_val(factor);
191  }
192 
193  if( conv.lgConvTemp )
194  break;
195 
196  if( trace.nTrConvg >= 2 && !conv.lgConvTemp && !lgAbort )
197  {
198  fprintf( ioQQQ, " ConvTempEdenIoniz: brent fails\n" );
199  TeTrack.print_history();
200  }
201  }
202  }
203 
204  if( lgAbort )
205  return 1;
206 
207  // only declare solution unstable if it is at least at the 2-sigma confidence level
208  thermal.lgUnstable = ( conv.dCmHdT + 2.*conv.sigma_dCmHdT < 0. );
209 
210  if( trace.lgTrace || trace.nTrConvg >= 2 )
211  {
212  fprintf( ioQQQ, " ConvTempEdenIoniz: Te %e C %.4e H %.4e (C-H)/H %.2f%%"
213  " d(C-H)/dT %.2e +/- %.2e\n",
215  (thermal.ctot/thermal.htot-1.)*100.,
217  fprintf( ioQQQ, " ConvTempEdenIoniz returns converged=%c\n", TorF(conv.lgConvTemp) );
218  }
219  }
220  else
221  {
222  fprintf( ioQQQ, "ConvTempEdenIoniz finds insane solver %s\n", conv.chSolverTemp );
223  ShowMe();
224  }
225 
226  return 0;
227 }
228 
229 
230 /* returns true if heating-cooling is converged */
231 STATIC bool lgConvTemp(const iter_track& TeTrack)
232 {
233  DEBUG_ENTRY( "lgConvTemp()" );
234 
235  if( lgAbort )
236  {
237  /* converging the temperature was aborted */
238  conv.lgConvTemp = false;
239  }
240  // The explicit test for H-C == 0. is needed since the requirement on the temperature
241  // bracket width may not be simultaneously satisfied. Since we have found the zero point
242  // exactly, we don't care about that... The temperature is as accurate as it is ever going
243  // to be. Not doing the explicit test for H-C == 0. is a bug. If the temparture bracket is
244  // too wide when we hit H-C == 0., this algorithm would never converge since vWDB requires
245  // that the endpoints of the bracket have non-zero function values, so they cannot get
246  // updated and the bracket width never gets smaller. So the requirement on the temperature
247  // bracket width would never be satisfied...
248  else if( thermal.htot - thermal.ctot == 0.
252  {
253  /* announce that temp is converged if relative heating - cooling mismatch
254  * is less than the relative heating cooling error allowed and the width of
255  * the temperature bracket is sufficiently small (this assures that the
256  * temperature is also well determined if H-C is a shallow function of T).
257  * If this is a constant temperature model, force convergence */
258  conv.lgConvTemp = true;
259  // remember numerical derivative to estimate initial stepsize on next call
260  conv.dCmHdT = TeTrack.deriv(conv.sigma_dCmHdT);
261  }
262  else
263  {
264  /* big mismatch, this has not converged */
265  conv.lgConvTemp = false;
266  }
267 
268  if( trace.nTrConvg >= 2 )
269  fprintf( ioQQQ, " lgConvTemp: C-H rel err %.4e Te rel err %.4e converged=%c\n",
271  TeTrack.bracket_width()/phycon.te,
272  TorF(conv.lgConvTemp) );
273 
274  return conv.lgConvTemp;
275 }
276 
277 /*CoolHeatError evaluate ionization, and difference in heating and cooling, for temperature temp */
278 STATIC double CoolHeatError( double temp )
279 {
280  DEBUG_ENTRY( "CoolHeatError()" );
281 
282  static ConvergenceCounter cctr=conv.register_("TEMP_CHANGES");
283  ++cctr;
284  TempChange( temp, false );
285 
286  /* converge the ionization and electron density;
287  * this calls ionize until lgIonDone is true */
288  /* NB should NOT set insanity - but rather return error condition */
289  if( ConvEdenIoniz() )
290  lgAbort = true;
291 
292  /* >>chng 01 mar 16, evaluate pressure here since changing and other values needed */
293  /* reevaluate pressure */
294  /* this sets values of pressure.PresTotlCurr */
295  PresTotCurrent();
296 
297  /* keep track of temperature solver in this zone
298  * conv.hist_temp_nzone is reset in ConvInitSolution */
299  if( nzone != conv.hist_temp_nzone )
300  {
301  /* first time in this zone - reset history */
303  conv.hist_temp_temp.clear();
304  conv.hist_temp_heat.clear();
305  conv.hist_temp_cool.clear();
306  }
307 
308  conv.hist_temp_temp.push_back( phycon.te );
309  conv.hist_temp_heat.push_back( thermal.htot );
310  conv.hist_temp_cool.push_back( thermal.ctot );
311 
312  // dump major contributors to heating and cooling - for debugging purposes
313  if( false )
314  {
317  }
318 
319  if( trace.nTrConvg >= 2 )
320  {
321 #if 0
322  species *spCO = findspecies("CO")->local()->dbase;
323 
324  fprintf( ioQQQ, "DEBUGGG CO %ld %ld\n",nzone,
325  spCO ? spCO->numLevels_local : -1);
326 #endif
327 
328  fprintf( ioQQQ, " CoolHeatError: Te: %.4e C: %.4e H: %.4e (C-H)/H: %.4e\n",
330  }
331 
332  double error = thermal.ctot - thermal.htot;
333 
334  // this can get set if temperature drops below floor temperature -> fake convergence
336  error = 0.;
337 
338  return error;
339 }
340 
341 STATIC void DumpCoolStack(double thres)
342 {
343  multimap<double,string> output;
344  char line[200];
345 
346  for( int i=0; i < thermal.ncltot; ++i )
347  {
348  double fraction;
349  if( abs(thermal.heatnt[i]) > thres )
350  {
351  fraction = thermal.heatnt[i]/thermal.ctot;
352  sprintf( line, "heat %s %e: %e %e\n",
353  thermal.chClntLab[i], thermal.collam[i], thermal.heatnt[i], fraction );
354  output.insert( pair<const double,string>( fraction, string(line) ) );
355  }
356  if( abs(thermal.cooling[i]) > thres )
357  {
358  fraction = thermal.cooling[i]/thermal.ctot;
359  sprintf( line, "cool %s %e: %e %e\n",
360  thermal.chClntLab[i], thermal.collam[i], thermal.cooling[i], fraction );
361  output.insert( pair<const double,string>( fraction, string(line) ) );
362  }
363  }
364 
365  dprintf( ioQQQ, " >>>>>>> STARTING COOLING DUMP <<<<<<\n" );
366  dprintf( ioQQQ, "total cooling %e\n", thermal.ctot );
367  // this will produce sorted output in reverse order (largest contributor first)
368  for( multimap<double,string>::reverse_iterator i=output.rbegin(); i != output.rend(); ++i )
369  dprintf( ioQQQ, "%s", i->second.c_str() );
370  dprintf( ioQQQ, " >>>>>>> FINISHED COOLING DUMP <<<<<<\n" );
371 }
372 
373 STATIC void DumpHeatStack(double thres)
374 {
375  multimap<double,string> output;
376  char line[200];
377 
378  for( int nelem=0; nelem < LIMELM; ++nelem )
379  {
380  for( int i=0; i < LIMELM; ++i )
381  {
382  double fraction = thermal.heating(nelem,i)/thermal.htot;
383  if( abs(thermal.heating(nelem,i)) > thres )
384  {
385  sprintf( line, "heating(%i,%i): %e %e\n",
386  nelem, i, thermal.heating(nelem,i), fraction );
387  output.insert( pair<const double,string>( fraction, string(line) ) );
388  }
389  }
390  }
391 
392  dprintf( ioQQQ, " >>>>>>> STARTING HEATING DUMP <<<<<<\n" );
393  dprintf( ioQQQ, "total heating %e\n", thermal.htot );
394  // this will produce sorted output in reverse order (largest contributor first)
395  for( multimap<double,string>::reverse_iterator i=output.rbegin(); i != output.rend(); ++i )
396  dprintf( ioQQQ, "%s", i->second.c_str() );
397  dprintf( ioQQQ, " >>>>>>> FINISHED HEATING DUMP <<<<<<\n" );
398 }
double Radius
Definition: radius.h:31
double depth
Definition: radius.h:31
double htot
Definition: thermal.h:169
t_thermal thermal
Definition: thermal.cpp:6
t_colden colden
Definition: colden.cpp:5
double bracket_width() const
Definition: iter_track.h:85
NORETURN void TotalInsanity(void)
Definition: service.cpp:1067
double next_val()
Definition: iter_track.cpp:57
void add(double x, double fx)
Definition: iter_track.h:120
void print_history() const
Definition: iter_track.h:186
ConvergenceCounter register_(const string name)
Definition: conv.cpp:88
STATIC bool lgConvTemp(const iter_track &TeTrack)
DepthTable tlaw
Definition: thermal.h:95
char TorF(bool l)
Definition: cddefines.h:753
STATIC void DumpHeatStack(double thres)
t_conv conv
Definition: conv.cpp:5
T sign(T x, T y)
Definition: cddefines.h:846
t_phycon phycon
Definition: phycon.cpp:6
realnum SigmaBD96
Definition: thermal.h:86
double cooling[NCOLNT]
Definition: thermal.h:104
bool lgTemperatureConstant
Definition: thermal.h:44
FILE * ioQQQ
Definition: cddefines.cpp:7
long int nzone
Definition: cddefines.cpp:14
STATIC void DumpCoolStack(double thres)
void TempChange(double TempNew, bool lgForceUpdate)
Definition: temp_change.cpp:31
int ConvTempEdenIoniz(void)
t_dense dense
Definition: global.cpp:15
void PresTotCurrent(void)
const double TEMP_LIMIT_LOW
Definition: phycon.h:121
vector< double > hist_temp_cool
Definition: conv.h:298
t_trace trace
Definition: trace.cpp:5
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:858
long int hist_temp_nzone
Definition: conv.h:299
bool lgUnstable
Definition: thermal.h:65
int init_bracket(double x1, double fx1, double x2, double fx2)
Definition: iter_track.h:104
int dprintf(FILE *fp, const char *format,...)
Definition: service.cpp:1198
char chClntLab[NCOLNT][NCOLNT_LAB_LEN+1]
Definition: thermal.h:108
#define STATIC
Definition: cddefines.h:118
bool lgTrace
Definition: trace.h:12
long numLevels_local
Definition: species.h:48
class molezone * local(void) const
Definition: mole.h:482
molecule * findspecies(const char buf[])
vector< double > hist_temp_temp
Definition: conv.h:298
bool lgTeTLaw
Definition: thermal.h:96
double sigma_dCmHdT
Definition: conv.h:289
double dCmHdT
Definition: conv.h:286
long max(int a, long b)
Definition: cddefines.h:821
double tabval(double r0, double depth) const
Definition: depth_table.cpp:8
bool lgTeBD96
Definition: thermal.h:84
long min(int a, long b)
Definition: cddefines.h:766
species * dbase
Definition: mole.h:416
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:1015
int nTrConvg
Definition: trace.h:27
realnum HeatCoolRelErrorAllowed
Definition: conv.h:276
t_radius radius
Definition: radius.cpp:5
double heating(long nelem, long ion)
Definition: thermal.h:186
realnum gas_phase[LIMELM]
Definition: dense.h:76
long int ncltot
Definition: thermal.h:106
realnum collam[NCOLNT]
Definition: thermal.h:103
bool lgConvTemp
Definition: conv.h:191
const int LIMELM
Definition: cddefines.h:307
void clear()
Definition: iter_track.h:76
bool lgTLaw
Definition: thermal.h:80
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
double H2_total
Definition: hmi.h:25
bool lgTeSN99
Definition: thermal.h:92
vector< double > hist_temp_heat
Definition: conv.h:298
realnum T0SN99
Definition: thermal.h:91
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
double deriv(int n, double &sigma) const
Definition: iter_track.cpp:184
t_hmi hmi
Definition: hmi.cpp:5
void ShowMe(void)
Definition: service.cpp:205
double te
Definition: phycon.h:21
const int ipHYDROGEN
Definition: cddefines.h:348
void set_tol(double tol)
Definition: iter_track.h:81
double heatnt[NCOLNT]
Definition: thermal.h:104
int ConvEdenIoniz(void)
realnum T0BD96
Definition: thermal.h:86
#define POW4
Definition: cddefines.h:997
realnum colden[NCOLD]
Definition: colden.h:32
char chSolverTemp[20]
Definition: conv.h:244
STATIC double CoolHeatError(double temp)
bool lgAbort
Definition: cddefines.cpp:10
double ctot
Definition: thermal.h:130