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 bool lgHighTeSearch = false;
00100
00101 t2 = phycon.te;
00102 error2 = CoolHeatError( t2 );
00103
00104 for( int n=0; n < 5 && !lgAbort; ++n )
00105 {
00106 const int DEF_ITER = 10;
00107 const double DEF_FACTOR = 0.02;
00108 double step, factor = DEF_FACTOR;
00109
00110 TeTrack.clear();
00111
00112
00113
00114
00115 for( int i=0; i < 2 && !lgAbort; ++i )
00116 {
00117 t1 = t2;
00118 error1 = error2;
00119
00120
00121 step = safe_div( -1.2*error1, conv.dCmHdT, 0. );
00122 step = sign( min( abs(step), factor*t1 ), step );
00123 t2 = t1 + step;
00124 error2 = CoolHeatError( t2 );
00125 TeTrack.add( t2, error2 );
00126 }
00127
00128 int j = 0;
00129
00130
00131 while( error1*error2 > 0. && j++ < 2*DEF_ITER && !lgAbort )
00132 {
00133 t1 = t2;
00134 error1 = error2;
00135 double deriv = TeTrack.deriv(7);
00136
00137 step = safe_div( -1.2*error1, deriv, 0. );
00138 step = sign( min( abs(step), factor*t1 ), step );
00139 t2 = t1 + step;
00140 error2 = CoolHeatError( t2 );
00141 TeTrack.add( t2, error2 );
00142 }
00143
00144 if( trace.nTrConvg >= 2 && error1*error2 > 0. && !lgAbort )
00145 {
00146 fprintf( ioQQQ, " ConvTempEdenIoniz: bracket1 fails t1: %e %e t2: %e %e\n",
00147 t1, error1, t2, error2 );
00148 TeTrack.print_history();
00149 }
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159 int nUP = 2;
00160
00161 if( error2 < 0. && !lgHighTeSearch && !lgAbort )
00162 {
00163 lgHighTeSearch = true;
00164 nUP = 25;
00165 }
00166
00167
00168
00169
00170
00171 while( error1*error2 > 0. && j++ < (2+nUP)*DEF_ITER && !lgAbort &&
00172 t2*(1.+factor) <= 0.99999*phycon.TEMP_LIMIT_HIGH )
00173 {
00174 t1 = t2;
00175 error1 = error2;
00176 if( lgHighTeSearch && t2 > 4000. )
00177 factor = 0.05;
00178 t2 = t1*(1.+factor);
00179 error2 = CoolHeatError( t2 );
00180 TeTrack.add( t2, error2 );
00181 }
00182
00183
00184
00185
00186
00187
00188
00189 while( error1*error2 > 0. && j++ < (2+2*nUP+9)*DEF_ITER && !lgAbort )
00190 {
00191 t1 = t2;
00192 error1 = error2;
00193 if( lgHighTeSearch && t2 <= 4000. )
00194 factor = DEF_FACTOR;
00195 t2 = t1*(1.-factor);
00196 error2 = CoolHeatError( t2 );
00197 TeTrack.add( t2, error2 );
00198 }
00199
00200 factor = DEF_FACTOR;
00201
00202 if( trace.nTrConvg >= 2 && error1*error2 > 0. && !lgAbort )
00203 {
00204 fprintf( ioQQQ, " ConvTempEdenIoniz: bracket2 fails t1: %e %e t2: %e %e\n",
00205 t1, error1, t2, error2 );
00206 TeTrack.print_history();
00207 }
00208
00209
00210
00211 TeTrack.clear();
00212
00213
00214 if( TeTrack.init_bracket( t1, error1, t2, error2 ) == 0 )
00215 {
00216
00217
00218
00219
00220
00221 TeTrack.set_tol(2.*DBL_EPSILON*t2);
00222
00223 t2 = 0.5*(t1+t2);
00224 for( int i = 0; i < (1<<(n/2))*DEF_ITER && !lgAbort; i++ )
00225 {
00226
00227 if( lgConvTemp(TeTrack) || TeTrack.bracket_width() < 3.*DBL_EPSILON*t2 )
00228 break;
00229
00230 error2 = CoolHeatError( t2 );
00231 TeTrack.add( t2, error2 );
00232 t2 = TeTrack.next_val(factor);
00233 }
00234
00235 if( conv.lgConvTemp )
00236 break;
00237
00238 if( trace.nTrConvg >= 2 && !lgAbort )
00239 {
00240 fprintf( ioQQQ, " ConvTempEdenIoniz: brent fails\n" );
00241 TeTrack.print_history();
00242 }
00243 }
00244 }
00245
00246 if( lgAbort )
00247 return 1;
00248
00249
00250 thermal.lgUnstable = ( conv.dCmHdT + 2.*conv.sigma_dCmHdT < 0. );
00251
00252 if( trace.lgTrace || trace.nTrConvg >= 2 )
00253 {
00254 fprintf( ioQQQ, " ConvTempEdenIoniz: Te %e C %.4e H %.4e (C-H)/H %.2f%%"
00255 " d(C-H)/dT %.2e +/- %.2e\n",
00256 phycon.te, thermal.ctot, thermal.htot,
00257 (thermal.ctot/thermal.htot-1.)*100.,
00258 conv.dCmHdT, conv.sigma_dCmHdT );
00259 fprintf( ioQQQ, " ConvTempEdenIoniz returns converged=%c\n", TorF(conv.lgConvTemp) );
00260 }
00261 }
00262 else
00263 {
00264 fprintf( ioQQQ, "ConvTempEdenIoniz finds insane solver %s\n", conv.chSolverTemp );
00265 ShowMe();
00266 }
00267
00268 return 0;
00269 }
00270
00271
00272
00273 STATIC bool lgConvTemp(const iter_track& TeTrack)
00274 {
00275 DEBUG_ENTRY( "lgConvTemp()" );
00276
00277 if( lgAbort )
00278 {
00279
00280 conv.lgConvTemp = false;
00281 }
00282 else if( ( abs(thermal.htot - thermal.ctot)/thermal.htot <= conv.HeatCoolRelErrorAllowed &&
00283 TeTrack.bracket_width()/phycon.te <= conv.HeatCoolRelErrorAllowed/3. ) ||
00284 thermal.lgTemperatureConstant || phycon.te <= phycon.TEMP_LIMIT_LOW )
00285 {
00286
00287
00288
00289 conv.lgConvTemp = true;
00290
00291 conv.dCmHdT = TeTrack.deriv(conv.sigma_dCmHdT);
00292 }
00293 else
00294 {
00295
00296 conv.lgConvTemp = false;
00297 }
00298
00299 if( trace.nTrConvg >= 2 )
00300 fprintf( ioQQQ, " lgConvTemp: C-H abs err %.4e Te err %.4e converged=%c\n",
00301 abs(thermal.htot - thermal.ctot)/thermal.htot,
00302 TeTrack.bracket_width()/phycon.te,
00303 TorF(conv.lgConvTemp) );
00304
00305 return conv.lgConvTemp;
00306 }
00307
00308
00309 STATIC double CoolHeatError( double temp )
00310 {
00311 DEBUG_ENTRY( "CoolHeatError()" );
00312
00313 TempChange( temp, false );
00314
00315
00316
00317
00318 if( ConvEdenIoniz() )
00319 lgAbort = true;
00320
00321
00322
00323
00324 PresTotCurrent();
00325
00326
00327
00328 if( nzone != conv.hist_temp_nzone )
00329 {
00330
00331 conv.hist_temp_nzone = nzone;
00332 conv.hist_temp_temp.clear();
00333 conv.hist_temp_heat.clear();
00334 conv.hist_temp_cool.clear();
00335 }
00336
00337 conv.hist_temp_temp.push_back( phycon.te );
00338 conv.hist_temp_heat.push_back( thermal.htot );
00339 conv.hist_temp_cool.push_back( thermal.ctot );
00340
00341
00342 if( false )
00343 {
00344 DumpCoolStack( conv.HeatCoolRelErrorAllowed/5.*thermal.ctot );
00345 DumpHeatStack( conv.HeatCoolRelErrorAllowed/5.*thermal.htot );
00346 }
00347
00348 if( trace.nTrConvg >= 2 )
00349 fprintf( ioQQQ, " CoolHeatError: Te: %.4e C: %.4e H: %.4e (C-H)/H: %.4e\n",
00350 temp, thermal.ctot, thermal.htot, thermal.ctot/thermal.htot-1. );
00351
00352 double error = thermal.ctot - thermal.htot;
00353
00354
00355 if( thermal.lgTemperatureConstant )
00356 error = 0.;
00357
00358 return error;
00359 }
00360
00361 STATIC void DumpCoolStack(double thres)
00362 {
00363 multimap<double,string> output;
00364 char line[200];
00365
00366 for( int i=0; i < thermal.ncltot; ++i )
00367 {
00368 double fraction;
00369 if( abs(thermal.heatnt[i]) > thres )
00370 {
00371 fraction = thermal.heatnt[i]/thermal.ctot;
00372 sprintf( line, "heat %s %e: %e %e\n",
00373 thermal.chClntLab[i], thermal.collam[i], thermal.heatnt[i], fraction );
00374 output.insert( pair<const double,string>( fraction, string(line) ) );
00375 }
00376 if( abs(thermal.cooling[i]) > thres )
00377 {
00378 fraction = thermal.cooling[i]/thermal.ctot;
00379 sprintf( line, "cool %s %e: %e %e\n",
00380 thermal.chClntLab[i], thermal.collam[i], thermal.cooling[i], fraction );
00381 output.insert( pair<const double,string>( fraction, string(line) ) );
00382 }
00383 }
00384
00385 dprintf( ioQQQ, " >>>>>>> STARTING COOLING DUMP <<<<<<\n" );
00386 dprintf( ioQQQ, "total cooling %e\n", thermal.ctot );
00387
00388 for( multimap<double,string>::reverse_iterator i=output.rbegin(); i != output.rend(); ++i )
00389 dprintf( ioQQQ, "%s", i->second.c_str() );
00390 dprintf( ioQQQ, " >>>>>>> FINISHED COOLING DUMP <<<<<<\n" );
00391 }
00392
00393 STATIC void DumpHeatStack(double thres)
00394 {
00395 multimap<double,string> output;
00396 char line[200];
00397
00398 for( int nelem=0; nelem < LIMELM; ++nelem )
00399 {
00400 for( int i=0; i < LIMELM; ++i )
00401 {
00402 double fraction = thermal.heating[nelem][i]/thermal.htot;
00403 if( abs(thermal.heating[nelem][i]) > thres )
00404 {
00405 sprintf( line, "heating[%i][%i]: %e %e\n",
00406 nelem, i, thermal.heating[nelem][i], fraction );
00407 output.insert( pair<const double,string>( fraction, string(line) ) );
00408 }
00409 }
00410 }
00411
00412 dprintf( ioQQQ, " >>>>>>> STARTING HEATING DUMP <<<<<<\n" );
00413 dprintf( ioQQQ, "total heating %e\n", thermal.htot );
00414
00415 for( multimap<double,string>::reverse_iterator i=output.rbegin(); i != output.rend(); ++i )
00416 dprintf( ioQQQ, "%s", i->second.c_str() );
00417 dprintf( ioQQQ, " >>>>>>> FINISHED HEATING DUMP <<<<<<\n" );
00418 }