00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "cddefines.h"
00010 #include "hmi.h"
00011 #include "thermal.h"
00012 #include "iso.h"
00013 #include "hydrogenic.h"
00014 #include "colden.h"
00015 #include "h2.h"
00016 #include "pressure.h"
00017 #include "dense.h"
00018 #include "struc.h"
00019 #include "thirdparty.h"
00020 #include "trace.h"
00021 #include "phycon.h"
00022 #include "conv.h"
00023
00024
00025 STATIC bool lgConvTemp(const iter_track& TeTrack);
00026
00027 STATIC double CoolHeatError( double temp );
00028
00029
00030 STATIC void DumpCoolStack(double thres);
00031 STATIC void DumpHeatStack(double thres);
00032
00033
00034
00035
00036 int ConvTempEdenIoniz( void )
00037 {
00038 DEBUG_ENTRY( "ConvTempEdenIoniz()" );
00039
00040 if( trace.lgTrace )
00041 {
00042 fprintf( ioQQQ, "\n ConvTempEdenIoniz called\n" );
00043 }
00044 if( trace.nTrConvg >= 2 )
00045 {
00046 fprintf( ioQQQ, " ConvTempEdenIoniz called, entering temp loop using solver %s.\n",
00047 conv.chSolverTemp );
00048 }
00049
00050 if( strcmp( conv.chSolverTemp , "vWDB" ) == 0 )
00051 {
00052 conv.lgConvTemp = false;
00053
00054
00055 if( thermal.lgTLaw )
00056 {
00057 double TeNew = phycon.te;
00058 if( thermal.lgTeBD96 )
00059 {
00060
00061 TeNew = thermal.T0BD96 / (1. + thermal.SigmaBD96 * colden.colden[ipCOL_HTOT]);
00062 }
00063 else if( thermal.lgTeSN99 )
00064 {
00065
00066
00067
00068 TeNew = thermal.T0SN99 /
00069 (1. + 9.*POW4(2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN]) );
00070 }
00071 else
00072 TotalInsanity();
00073
00074 TempChange( TeNew, false );
00075 }
00076
00077 if( thermal.lgTemperatureConstant || thermal.lgTLaw )
00078 {
00079 if( ConvEdenIoniz() )
00080 return 1;
00081 PresTotCurrent();
00082
00083
00084 conv.lgConvTemp = true;
00085
00086 if( trace.lgTrace || trace.nTrConvg >= 2 )
00087 {
00088 fprintf( ioQQQ, " ConvTempEdenIoniz: Te %e C %.4e H %.4e\n",
00089 phycon.te, thermal.ctot, thermal.htot );
00090 fprintf( ioQQQ, " ConvTempEdenIoniz returns ok.\n" );
00091 }
00092
00093 return 0;
00094 }
00095
00096
00097 iter_track TeTrack;
00098 double t1=0, error1=0, t2, error2;
00099
00100 t2 = phycon.te;
00101 error2 = CoolHeatError( t2 );
00102
00103 for( int n=0; n < 5 && !lgAbort; ++n )
00104 {
00105 const int DEF_ITER = 10;
00106 const double DEF_FACTOR = 0.2;
00107 double step, factor = DEF_FACTOR;
00108
00109 TeTrack.clear();
00110
00111 step = min( abs(safe_div( error2, conv.dCmHdT, 0. )), factor*t2 );
00112
00113
00114
00115
00116 for( int i=0; i < 100 && !lgAbort; ++i )
00117 {
00118 t1 = t2;
00119 error1 = error2;
00120
00121 double maxstep = factor*t1;
00122
00123
00124 step = SQRT2*step;
00125 if( step == 0.0 || step > maxstep )
00126 step = maxstep;
00127 t2 = max( t1 + sign( step, -error1 ), phycon.TEMP_LIMIT_LOW );
00128 error2 = CoolHeatError( t2 );
00129 TeTrack.add( t2, error2 );
00130
00131
00132
00133
00134
00135 if( i >= n && error1*error2 <= 0. )
00136 break;
00137
00138
00139
00140 if( i >= n && fp_equal( t2, phycon.TEMP_LIMIT_LOW ) )
00141 {
00142
00143 fprintf(ioQQQ," PROBLEM DISASTER - the lower limit of the code,"
00144 " %.3eK, does not bracket thermal balance.\n",
00145 phycon.TEMP_LIMIT_LOW );
00146 fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n");
00147 lgAbort = true;
00148 return 1;
00149 }
00150 }
00151
00152 if( trace.nTrConvg >= 2 && error1*error2 > 0. && !lgAbort )
00153 {
00154 fprintf( ioQQQ, " ConvTempEdenIoniz: bracket1 fails t1: %e %e t2: %e %e\n",
00155 t1, error1, t2, error2 );
00156 TeTrack.print_history();
00157 }
00158
00159
00160
00161 TeTrack.clear();
00162
00163
00164 if( TeTrack.init_bracket( t1, error1, t2, error2 ) == 0 )
00165 {
00166
00167
00168
00169
00170
00171 TeTrack.set_tol(2.*DBL_EPSILON*t2);
00172
00173 if( error1 != 0.0 || error2 != 0.0 )
00174 t2 = (t1*error2-t2*error1)/(error2-error1);
00175 else
00176 t2 = 0.5*(t1+t2);
00177
00178 for( int i = 0; i < (1<<(n/2))*DEF_ITER && !lgAbort; i++ )
00179 {
00180
00181 if( lgConvTemp(TeTrack) || TeTrack.bracket_width() < 3.*DBL_EPSILON*t2 )
00182 break;
00183
00184 error2 = CoolHeatError( t2 );
00185 TeTrack.add( t2, error2 );
00186 t2 = TeTrack.next_val(factor);
00187 }
00188
00189 if( conv.lgConvTemp )
00190 break;
00191
00192 if( trace.nTrConvg >= 2 && !conv.lgConvTemp && !lgAbort )
00193 {
00194 fprintf( ioQQQ, " ConvTempEdenIoniz: brent fails\n" );
00195 TeTrack.print_history();
00196 }
00197 }
00198 }
00199
00200 if( lgAbort )
00201 return 1;
00202
00203
00204 thermal.lgUnstable = ( conv.dCmHdT + 2.*conv.sigma_dCmHdT < 0. );
00205
00206 if( trace.lgTrace || trace.nTrConvg >= 2 )
00207 {
00208 fprintf( ioQQQ, " ConvTempEdenIoniz: Te %e C %.4e H %.4e (C-H)/H %.2f%%"
00209 " d(C-H)/dT %.2e +/- %.2e\n",
00210 phycon.te, thermal.ctot, thermal.htot,
00211 (thermal.ctot/thermal.htot-1.)*100.,
00212 conv.dCmHdT, conv.sigma_dCmHdT );
00213 fprintf( ioQQQ, " ConvTempEdenIoniz returns converged=%c\n", TorF(conv.lgConvTemp) );
00214 }
00215 }
00216 else
00217 {
00218 fprintf( ioQQQ, "ConvTempEdenIoniz finds insane solver %s\n", conv.chSolverTemp );
00219 ShowMe();
00220 }
00221
00222 return 0;
00223 }
00224
00225
00226
00227 STATIC bool lgConvTemp(const iter_track& TeTrack)
00228 {
00229 DEBUG_ENTRY( "lgConvTemp()" );
00230
00231 if( lgAbort )
00232 {
00233
00234 conv.lgConvTemp = false;
00235 }
00236
00237
00238
00239
00240
00241
00242
00243
00244 else if( thermal.htot - thermal.ctot == 0.
00245 || ( abs(thermal.htot - thermal.ctot)/thermal.htot <= conv.HeatCoolRelErrorAllowed &&
00246 TeTrack.bracket_width()/phycon.te <= conv.HeatCoolRelErrorAllowed/3. )
00247 || thermal.lgTemperatureConstant )
00248 {
00249
00250
00251
00252
00253
00254 conv.lgConvTemp = true;
00255
00256 conv.dCmHdT = TeTrack.deriv(conv.sigma_dCmHdT);
00257 }
00258 else
00259 {
00260
00261 conv.lgConvTemp = false;
00262 }
00263
00264 if( trace.nTrConvg >= 2 )
00265 fprintf( ioQQQ, " lgConvTemp: C-H rel err %.4e Te rel err %.4e converged=%c\n",
00266 abs(thermal.htot - thermal.ctot)/thermal.htot,
00267 TeTrack.bracket_width()/phycon.te,
00268 TorF(conv.lgConvTemp) );
00269
00270 return conv.lgConvTemp;
00271 }
00272
00273
00274 STATIC double CoolHeatError( double temp )
00275 {
00276 DEBUG_ENTRY( "CoolHeatError()" );
00277
00278 conv.incrementCounter(TEMP_CHANGES);
00279 TempChange( temp, false );
00280
00281
00282
00283
00284 if( ConvEdenIoniz() )
00285 lgAbort = true;
00286
00287
00288
00289
00290 PresTotCurrent();
00291
00292
00293
00294 if( nzone != conv.hist_temp_nzone )
00295 {
00296
00297 conv.hist_temp_nzone = nzone;
00298 conv.hist_temp_temp.clear();
00299 conv.hist_temp_heat.clear();
00300 conv.hist_temp_cool.clear();
00301 }
00302
00303 conv.hist_temp_temp.push_back( phycon.te );
00304 conv.hist_temp_heat.push_back( thermal.htot );
00305 conv.hist_temp_cool.push_back( thermal.ctot );
00306
00307
00308 if( false )
00309 {
00310 DumpCoolStack( conv.HeatCoolRelErrorAllowed/5.*thermal.ctot );
00311 DumpHeatStack( conv.HeatCoolRelErrorAllowed/5.*thermal.htot );
00312 }
00313
00314 if( trace.nTrConvg >= 2 )
00315 fprintf( ioQQQ, " CoolHeatError: Te: %.4e C: %.4e H: %.4e (C-H)/H: %.4e\n",
00316 temp, thermal.ctot, thermal.htot, thermal.ctot/thermal.htot-1. );
00317
00318 double error = thermal.ctot - thermal.htot;
00319
00320
00321 if( thermal.lgTemperatureConstant )
00322 error = 0.;
00323
00324 return error;
00325 }
00326
00327 STATIC void DumpCoolStack(double thres)
00328 {
00329 multimap<double,string> output;
00330 char line[200];
00331
00332 for( int i=0; i < thermal.ncltot; ++i )
00333 {
00334 double fraction;
00335 if( abs(thermal.heatnt[i]) > thres )
00336 {
00337 fraction = thermal.heatnt[i]/thermal.ctot;
00338 sprintf( line, "heat %s %e: %e %e\n",
00339 thermal.chClntLab[i], thermal.collam[i], thermal.heatnt[i], fraction );
00340 output.insert( pair<const double,string>( fraction, string(line) ) );
00341 }
00342 if( abs(thermal.cooling[i]) > thres )
00343 {
00344 fraction = thermal.cooling[i]/thermal.ctot;
00345 sprintf( line, "cool %s %e: %e %e\n",
00346 thermal.chClntLab[i], thermal.collam[i], thermal.cooling[i], fraction );
00347 output.insert( pair<const double,string>( fraction, string(line) ) );
00348 }
00349 }
00350
00351 dprintf( ioQQQ, " >>>>>>> STARTING COOLING DUMP <<<<<<\n" );
00352 dprintf( ioQQQ, "total cooling %e\n", thermal.ctot );
00353
00354 for( multimap<double,string>::reverse_iterator i=output.rbegin(); i != output.rend(); ++i )
00355 dprintf( ioQQQ, "%s", i->second.c_str() );
00356 dprintf( ioQQQ, " >>>>>>> FINISHED COOLING DUMP <<<<<<\n" );
00357 }
00358
00359 STATIC void DumpHeatStack(double thres)
00360 {
00361 multimap<double,string> output;
00362 char line[200];
00363
00364 for( int nelem=0; nelem < LIMELM; ++nelem )
00365 {
00366 for( int i=0; i < LIMELM; ++i )
00367 {
00368 double fraction = thermal.heating[nelem][i]/thermal.htot;
00369 if( abs(thermal.heating[nelem][i]) > thres )
00370 {
00371 sprintf( line, "heating[%i][%i]: %e %e\n",
00372 nelem, i, thermal.heating[nelem][i], fraction );
00373 output.insert( pair<const double,string>( fraction, string(line) ) );
00374 }
00375 }
00376 }
00377
00378 dprintf( ioQQQ, " >>>>>>> STARTING HEATING DUMP <<<<<<\n" );
00379 dprintf( ioQQQ, "total heating %e\n", thermal.htot );
00380
00381 for( multimap<double,string>::reverse_iterator i=output.rbegin(); i != output.rend(); ++i )
00382 dprintf( ioQQQ, "%s", i->second.c_str() );
00383 dprintf( ioQQQ, " >>>>>>> FINISHED HEATING DUMP <<<<<<\n" );
00384 }