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