cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*ConvTempEdenIoniz determine temperature, called by ConPresTempEdenIoniz, 00004 * calls ConvEdenIoniz to get electron density and ionization */ 00005 /*MakeDeriv derive numerical derivative of heating minus cooling */ 00006 /*PutHetCol save heating, cooling, and temperature in stack for numerical derivatives */ 00007 /*lgConvTemp returns true if heating-cooling is converged */ 00008 /*CoolHeatError evaluate ionization, and difference in heating and cooling, for temperature temp */ 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 "trace.h" 00019 #include "phycon.h" 00020 #include "conv.h" 00021 00022 /* >>chng 01 apr 06, from 10 to 20, also put damper at 0.5 below, 00023 * as part of strategy of letting code bracket solution */ 00024 /* >>chng 04 mar 17, from 20 to 40, now doing much finer changes in 00025 * many conditions, including smaller dT, so allow more loops */ 00026 /* >>chng 05 mar 25, from 40 to 60, thermal fronts */ 00027 static const int LIMDEF = 60; 00028 00029 /*CoolHeatError evaluate ionization, and difference in heating and cooling, for temperature temp */ 00030 STATIC double CoolHeatError( double temp ); 00031 00032 /* use brent's method to find heating-cooling match */ 00033 STATIC double TeBrent( 00034 double x1, 00035 double x2); 00036 00037 /*MakeDeriv derive numerical derivative of heating minus cooling */ 00038 STATIC void MakeDeriv( 00039 /* which job to do, either "zero" or "incr" */ 00040 const char *job, 00041 /* result, the numerical derivative */ 00042 double *DerivNumer); 00043 00044 /*PutHetCol save heating, cooling, and temperature in stack for numerical derivatives */ 00045 STATIC void PutHetCol(double te, 00046 double htot, 00047 double ctot); 00048 00049 /*ConvTempEdenIoniz determine temperature, called by ConPresTempEdenIoniz, 00050 * calls ConvEdenIoniz to get electron density and ionization 00051 * returns 0 if ok, 1 if disaster */ 00052 int ConvTempEdenIoniz( void ) 00053 { 00054 char chDEType; 00055 long int limtry, 00056 nTemperature_loop; 00057 /* following is zero because there are paths were not set following abort */ 00058 double CoolMHeat=0., 00059 Damper, 00060 dtl, 00061 fneut; 00062 int kase; 00063 double DerivNumer; 00064 static double OldCmHdT = 0.; 00065 long int nTemp_oscil; 00066 00067 /* derivative of cooling less heating wrt temperature */ 00068 double dCoolmHeatdT; 00069 00070 DEBUG_ENTRY( "ConvTempEdenIoniz()" ); 00071 00072 /* determine ionization and temperature, called by PIonte 00073 * nCall tells routine which call this is, 0 for first one, not changed here 00074 * nTemperature_loop counts how many loops performed */ 00075 00076 if( trace.lgTrace ) 00077 { 00078 fprintf( ioQQQ, "\n ConvTempEdenIoniz called\n" ); 00079 } 00080 if( trace.nTrConvg>=2 ) 00081 { 00082 fprintf( ioQQQ, " ConvTempEdenIoniz1 called. Te:%.4e Htot:%11.4e Ctot:%11.4e error:%9.1e%%, eden=%11.4e\n", 00083 phycon.te, thermal.htot, thermal.ctot, (thermal.htot - thermal.ctot)* 00084 100./MAX2(1e-36,thermal.htot), dense.eden ); 00085 } 00086 00087 if( strcmp( conv.chSolverTemp , "new" ) == 0 ) 00088 { 00089 /* this has never worked */ 00090 double t1 , t2 , 00091 error1 , error2; 00092 double TeNew; 00093 double factor = 1.02; 00094 00095 t1 = phycon.te; 00096 error1 = CoolHeatError( t1 ); 00097 t2 = t1 * factor; 00098 error2 = CoolHeatError( t2 ); 00099 TeNew = TeBrent( t1 , t2 ); 00100 fprintf(ioQQQ , "DEBUG new temp solver error1 %.2e error2 %.2e TeNew %.2e\n", 00101 error1 , error2 , TeNew ); 00102 } 00103 00104 else if( strcmp( conv.chSolverTemp , "simple") == 0 ) 00105 { 00106 00107 /* main routine to find ionization and temperature 00108 * following is flag to check for temp oscillations 00109 * it could be set true in main temp loop */ 00110 conv.lgTOscl = false; 00111 /* ots rates not oscillating */ 00112 conv.lgOscilOTS = false; 00113 00114 /* scale factor for changes in temp - make smaller if oscillations */ 00115 Damper = 1.; 00116 00117 /* flag for d(cool - heat)/dT changing sign, an internal error 00118 * >>chng 97 sep 03, only turn false one time, remains true if ever true */ 00119 if( nzone < 1 ) 00120 { 00121 conv.lgCmHOsc = false; 00122 } 00123 00124 /* option to make numerical derivatives of heating cooling */ 00125 MakeDeriv("zero",&DerivNumer); 00126 00127 /* this is will initial limit to number of tries at temp */ 00128 limtry = LIMDEF; 00129 00130 /* used to keep track of sign of dt, to check for oscillations */ 00131 dtl = 0.; 00132 00133 /* this counts how may thermal loops we go through 00134 * used in calling routine to make sure that at least 00135 * two solutions are performed */ 00136 nTemperature_loop = 0; 00137 nTemp_oscil = 0; 00138 00139 /* this is change in temp, which is used early in following loop, must be preset */ 00140 thermal.dTemper = 1e-3f*phycon.te; 00141 00142 /* this is the start of the heating-cooling match loop, 00143 * this exits by returning in the middle */ 00144 while ( true ) 00145 { 00146 00147 /* >>chng 02 dec 03 rjrw, insert these so dynamics cooling rate is calculated 00148 * before step is chosen -- but how much duplication results? 00149 * must have current estimate of heating and cooling before temperature is changed. 00150 * at bottom of following loop these two are called again to see whether heatin cooling 00151 * have converged */ 00152 if( ConvEdenIoniz() ) 00153 { 00154 return 1; 00155 } 00156 PresTotCurrent(); 00157 00158 /* we will try to make the following zero */ 00159 CoolMHeat = thermal.ctot - thermal.htot; 00160 /*fprintf(ioQQQ,"DEBUG grn 13 14\t%e\t%e\t%e\n", 00161 phycon.te, 00162 thermal.heating[0][13], 00163 thermal.heating[0][14] );*/ 00164 00165 if( thermal.lgTemperatureConstant ) 00166 { 00167 /* constant temp - declare this a success */ 00168 CoolMHeat = 0.; 00169 } 00170 00171 if( thermal.lgTLaw ) 00172 { 00173 double TeNew = phycon.te; 00174 /* temperature law has been specified */ 00175 if( thermal.lgTeBD96 ) 00176 { 00177 /* Bertoldi & Drain 96 temp law specified by TLAW BD96 command */ 00178 TeNew = thermal.T0BD96 / (1.f + thermal.SigmaBD96 * colden.colden[ipCOL_HTOT] ); 00179 } 00180 else if( thermal.lgTeSN99 ) 00181 { 00182 /* Sternberg & Neufeld 99 temp law specified by TLAW SN99 command, 00183 * this is equation 16 of 00184 * >>refer H2 temperature law Sternberg, A., & Neufeld, D.A. 1999, ApJ, 516, 371-380 */ 00185 TeNew = thermal.T0SN99 / 00186 /*(realnum)(1.f + 9.*pow(2.*hmi.Hmolec[ipMH2g]/dense.gas_phase[ipHYDROGEN], 4.) );*/ 00187 (realnum)(1.f + 9.*pow(2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN], 4.) ); 00188 } 00189 else 00190 TotalInsanity(); 00191 00192 TempChange(TeNew ,false); 00193 } 00194 00195 /* >>chng 02 dec 09, this block moved to just after loop entry above 00196 * to check on heating cooling and 00197 * exit if match is ok - prevents unnecessary reevaluate of ConvEdenIoniz and PresTotCurrent */ 00198 /* this is test for valid thermal solution, 00199 * lgConvTemp returns true if difference in heat cool is small */ 00200 /* >>chng 02 dec 04, add check on size of dTemper rel to te, 00201 * can become too small to increment a float */ 00202 /* >>chng 04 dec 13, for very cold gas the dTemper < 1e-6 Te limit was hit 00203 * during normal search. do not let this limit be tested here, rather use 00204 * fltepsilon to set limit to dTemper */ 00205 if( (lgConvTemp() /*|| fabs(thermal.dTemper) < 1e-6*phycon.te*/) 00206 && conv.lgConvIoniz ) 00207 { 00208 /* we have a good solution */ 00209 if( trace.lgTrace ) 00210 fprintf( ioQQQ, " ConvTempEdenIoniz returns ok.\n" ); 00211 00212 /* test for thermal stability, set flag if unstable, 00213 * later increment number of zones only once per zone */ 00214 /* this may become a check on thermal stability since neg slope is 00215 * unstable, but we may not have a valid solution, so only set flag */ 00216 if( thermal.dCooldT < 0. ) 00217 { 00218 thermal.lgUnstable = true; 00219 } 00220 else 00221 { 00222 thermal.lgUnstable = false; 00223 } 00224 00225 /* remember the highest and lowest temperatures */ 00226 thermal.thist = (realnum)MAX2(phycon.te,thermal.thist); 00227 thermal.tlowst = (realnum)MIN2(phycon.te,thermal.tlowst); 00228 00229 /* say that we have a valid solution */ 00230 conv.lgConvTemp = true; 00231 00232 if( trace.nTrConvg>=2 ) 00233 { 00234 fprintf( ioQQQ, " ConvTempEdenIoniz2 converged. Te:%11.4e Htot:%11.4e Ctot:%11.4e error:%9.1e%%, eden=%11.4e\n", 00235 phycon.te, thermal.htot, thermal.ctot, (thermal.htot - thermal.ctot)* 00236 100./MAX2(1e-36,thermal.htot), dense.eden ); 00237 } 00238 /******************************************************* 00239 * 00240 * this is return for valid solution 00241 * 00242 *******************************************************/ 00243 return 0; 00244 } 00245 else if(nTemperature_loop >= limtry || lgAbort ) 00246 { 00247 /* >>chng 02 dec 17, add test on ots fields oscillating 00248 * while ionization is not converged */ 00249 /* >>chng 04 may 25, do not test on lgOscilOTS since we do test 00250 * on lgConvIoniz */ 00251 if( nTemperature_loop == LIMDEF && !conv.lgTOscl && 00252 conv.lgConvIoniz /*&& !conv.lgOscilOTS*/ ) 00253 { 00254 /* LIMDEF is set to 60 above 00255 * >>chng 97 aug 10, from 2 to 4 to let keep going when not oscillation 00256 * this helped in getting over huge jumps in temperature at metal fronts */ 00257 limtry = LIMDEF*4; 00258 if( trace.nTrConvg>2 ) 00259 { 00260 fprintf( ioQQQ, " ConvTempEdenIoniz9 increases limtry to%3ld\n", 00261 limtry ); 00262 } 00263 } 00264 else 00265 { 00266 /* looped too long, so bail out anyway */ 00267 break; 00268 } 00269 } 00270 00271 /* get numerical heating-cooling derivative, but we may not use it */ 00272 MakeDeriv("incr",&DerivNumer); 00273 00274 /* remember this temp, heating, and cooling */ 00275 PutHetCol(phycon.te,thermal.htot,thermal.ctot); 00276 00277 /* the flag conv.lgConvEden was set during above call to ConvEdenIoniz, and 00278 * if is not true then eden failed, and we abort, 00279 * but don't abort on very first try, when flag has not been set 00280 * >> chng 02 dec 10 rjrw flag is now set on first try */ 00281 if( !conv.lgConvEden ) 00282 { 00283 /* set this to zero so that will not flag temperature failure, since we are 00284 * aborting due to eden failure */ 00285 CoolMHeat = 0.; 00286 if( trace.nTrConvg>=2 ) 00287 { 00288 fprintf( ioQQQ, 00289 " ConvTempEdenIoniz3 breaks since lgConvEden returns failed electron density.\n"); 00290 } 00291 break; 00292 } 00293 00294 /* save old derivative to check for oscillations */ 00295 if( nzone < 1 ) 00296 { 00297 /* >>chng 96 jun 07, was set to dCoolmHeatdT, is not defined initially */ 00298 OldCmHdT = thermal.ctot/phycon.te; 00299 } 00300 00301 # define USENUMER false 00302 if( USENUMER ) 00303 { 00304 dCoolmHeatdT = DerivNumer; 00305 chDEType = 'n'; 00306 } 00307 else if( phycon.te < 1e4 && thermal.dCooldT < 0. ) 00308 { 00309 /* use numerical guess, HTOT nearly equal to CTOT, drive Te lower 00310 * this is to overcome thermal inflection points */ 00311 dCoolmHeatdT = (double)(thermal.htot)/phycon.te*2.; 00312 chDEType = 'd'; 00313 } 00314 else 00315 { 00316 /* >>chng 97 mar 18, all below growth code just did the following 00317 * use analytical estimate of change in heat - cooling wrt temp */ 00318 dCoolmHeatdT = thermal.dCooldT - thermal.dHeatdT; 00319 chDEType = 'a'; 00320 } 00321 00322 /* check whether derivative changing sign - an internal error or 00323 * numerical instability */ 00324 if( OldCmHdT*dCoolmHeatdT < 0. ) 00325 { 00326 conv.lgCmHOsc = true; 00327 00328 /* derivative can change sign when heating and cooling balance, 00329 * and processes have exactly same temp dependence 00330 * happened in mods of ism.in test 00331 * >>chng 97 sep 02, added following test for oscillating derivs */ 00332 dCoolmHeatdT = thermal.ctot/phycon.te; 00333 chDEType = 'o'; 00334 } 00335 OldCmHdT = dCoolmHeatdT; 00336 00337 /* if not first pass through then numerical derivative should 00338 * not be zero required change in temperature, this may be too large */ 00339 thermal.dTemper = (realnum)(CoolMHeat/dCoolmHeatdT); 00340 kase = 0; 00341 00342 /* try to stop numerical oscillation if dTemper is changing sign */ 00343 /* when first time through and far from solution (when zones a bit too thick) first 00344 * step can be in wrong direction, and all subsequent are in right direction. 00345 * do not call this an oscillation - only declare oscillation when second occurs, 00346 * nTemperature_loop is zero first time through, (dtl is also zero), 00347 * is 1 first time through with memory of this zone, (when false oscillation occurs) 00348 * and is >1 thereafter */ 00349 /* >>chng 05 mar 25, had counted as oscillations when loop count greater than 1, change to 3 00350 * to allow things to settle down. Also, for first three evaluations, do not let te 00351 * change by too much, until lower solvers settle down */ 00352 if( dtl*thermal.dTemper < 0. && nTemperature_loop > 3) 00353 { 00354 conv.lgTOscl = true; 00355 /* make DT smaller and smaller */ 00356 /* >>chng 01 mar 16, from 0.5 to 0.75 00357 Damper *= 0.75;*/ 00358 Damper *= 0.5; 00359 nTemp_oscil = nTemperature_loop; 00360 } 00361 else if( Damper < 1. && nTemperature_loop-nTemp_oscil > 10 ) 00362 { 00363 /* >>chng 05 mar 26, some cases there is noise in heating-cooling 00364 * in first few evaluations and this causes te oscillations. but it may 00365 * settle down and become well behaved. this allows that to happen */ 00366 conv.lgTOscl = false; 00367 if( !((nTemperature_loop-nTemp_oscil)%5) ) 00368 Damper = MIN2( 1., Damper*1.5 ); 00369 } 00370 00371 /* SIGN: sign of sec arg and abs val of first */ 00372 /* >>chng 96 nov 08, from 0.1 to 0.08 in fneut */ 00373 /*fneut = (dense.xIonDense[ipHYDROGEN][0] + 2.*hmi.Hmolec[ipMH2g])/dense.gas_phase[ipHYDROGEN];*/ 00374 fneut = (dense.xIonDense[ipHYDROGEN][0] + 2.*hmi.H2_total)/dense.gas_phase[ipHYDROGEN]; 00375 /*>>chng 04 jan 11, when big H2 is on and heating is important, 00376 * only allow small changes in temp since must converge pops */ 00377 /* >>chng 04 jan 15, chng ratio of heating in h2 to total as the criteria, 00378 * to the error. temp changes will become much larger when this is 00379 * not tripped, and h2 can start to oscillate wildly. 00380 * NB this limit should be kept parallel with similar test for convergence 00381 * of heating in h2 */ 00382 if( 0 && h2.lgH2ON && fabs(hmi.HeatH2Dexc_BigH2)/thermal.ctot > conv.HeatCoolRelErrorAllowed ) 00383 { 00384 /* heating due to collisional deexcitation of the large H2 molecule 00385 * is very important. this rate depends strongly on the temperature 00386 * due to Boltzmann factors, and large changes in temperature 00387 * inhibit convergence of the level pops in the big H2. That 00388 * atom has tests on converging the heating, so want only small 00389 * temperature changes when H2 heating is important. */ 00390 double absdt = fabs(thermal.dTemper); 00391 thermal.dTemper = (sign(MIN2(absdt,0.0025*phycon.te), 00392 thermal.dTemper)); 00393 kase = 1; 00394 } 00395 else if( fneut > 0.08 && hydro.lgHColionImp ) 00396 { 00397 /* lgHColionImp is true if model collisionally ionized 00398 * do not let TE change by too much if lgHColionImp is set by HLEVEL */ 00399 double absdt = fabs(thermal.dTemper); 00400 thermal.dTemper = (sign(MIN2(absdt,0.005*phycon.te), 00401 thermal.dTemper)); 00402 kase = 2; 00403 } 00404 # if 0 00405 /* >>chng 05 mar 25, for first few steps through solvers in a new zone, don't let temp 00406 * change by too much initially, to avoid a te oscil which sets the oscil flag */ 00407 else if( 0 && nTemperature_loop <= 3) 00408 { 00409 /* >>chng 05 mar 25, add this branch, an early pass through these solvers, 00410 * do not let te change by too much to avoid te oscillations */ 00411 double absdt = fabs(thermal.dTemper); 00412 thermal.dTemper = (sign(MIN2(absdt,0.002*phycon.te), 00413 thermal.dTemper)); 00414 kase = 3; 00415 } 00416 # endif 00417 else 00418 { 00419 /* this branch, we will use the change in temperature that came from the deriv, 00420 * but do not let temp change by more than 2 percent */ 00421 double absdt = fabs(thermal.dTemper); 00422 thermal.dTemper = (sign(MIN2(absdt,0.02*phycon.te), 00423 thermal.dTemper)); 00424 kase = 4; 00425 } 00426 00427 /* >>chng 04 mar 14, add this test, cooling flow clouds had temp failure 00428 * when dt as large as 4% were allowed with the big H2 */ 00429 /* in no case do we want the temperature to change by more than 1% when the big 00430 * H2 molecule is turned on - it needs small temp changes to converge */ 00431 if( 0 && h2.lgH2ON && fabs(thermal.dTemper)/ phycon.te > 0.01 ) 00432 { 00433 double absdt = fabs(thermal.dTemper); 00434 thermal.dTemper = (sign(MIN2(absdt,0.01*phycon.te), 00435 thermal.dTemper)); 00436 kase = 5; 00437 } 00438 00439 /* >>chng 03 mar 16, to following logic, previous had been impossible 00440 * to hit - only do fine changes in Te when at half-way point in i-front 00441 * or if model is in collisional recombination equilibrium 00442 * if so then small change in temp causes big change in equilibrium 00443 * lgHColRec set in hlevel, only true is totl rec >> rad rec */ 00444 /* are most recombinations due to 3-body? */ 00445 if( dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] < 0.9 && 00446 /* is hydrogen more than 30% ionized */ 00447 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] > 0.01 ) 00448 { 00449 if( iso.RecomCollisFrac[ipH_LIKE][ipHYDROGEN] > 0.8 ) 00450 /* >>chng 04 feb 07, add test for H mostly collisionally ionized */ 00451 /* is hydrogen less than 70% ionized? */ 00452 /* >>chng 04 feb 07, from 70% to 90% ionized - oscill in hard 00453 * x-ray continua (laser at 80 ryd) when i front hit */ 00454 { 00455 double absdt = fabs(thermal.dTemper); 00456 /* >>chng 03 nov 03 , fac from 0.001 to 0.004 */ 00457 /* >>chng 04 feb 24 , fac back from 0.004 to 0.001, 00458 * this is needed for col den 25, high density (14) blr models 00459 * with flux of 22 - small changes in temperature drive large 00460 * changes in electron density due to three body recomb 00461 * this small change is ok due to test above, which makes sure 00462 * that we are in collisional recombination limit 00463 * sign returns sign of sec par times value of first */ 00464 thermal.dTemper = (sign(MIN2(absdt,0.001* 00465 phycon.te),thermal.dTemper)); 00466 kase = 6; 00467 } 00468 } 00469 00470 /* make steps smaller if oscillations, caused by overshots, occur */ 00471 thermal.dTemper *= Damper; 00472 00473 /* do not let dTemper/Te grow smaller than 10 * flt epsilon, the smaller number that 00474 * can be added to a float and still work */ 00475 if( fabs(thermal.dTemper)/phycon.te <10.*DBL_EPSILON ) 00476 { 00477 kase = 7; 00478 thermal.dTemper = (sign( 10.*DBL_EPSILON*phycon.te , thermal.dTemper)); 00479 } 00480 00481 /* save last change in temperature */ 00482 dtl = thermal.dTemper; 00483 00484 /* THIS IS IT !!! this is it !!! this is where te changes. */ 00485 if( !thermal.lgTLaw ) 00486 { 00487 /* usual case - valid thermal solution */ 00488 TempChange(phycon.te - thermal.dTemper , false); 00489 } 00490 00491 if( trace.nTrConvg>=2 ) 00492 { 00493 fprintf( ioQQQ, 00494 " ConvTempEdenIoniz4 a %4li T:%.4e dt:%10.2e%7.2f%% C:%10.2e H:%10.2e" 00495 " CoolMHeat/h:%7.2f%% dC-H/dT:%.2e kase:%i%c damp%.2f\n", 00496 nTemperature_loop, 00497 phycon.te, 00498 thermal.dTemper, 00499 thermal.dTemper/(phycon.te+thermal.dTemper)*100, 00500 thermal.ctot, 00501 thermal.htot, 00502 CoolMHeat/MIN2(thermal.ctot , thermal.htot )*100. , 00503 dCoolmHeatdT, 00504 kase, 00505 /* which type of deriv we have */ 00506 chDEType, 00507 Damper); 00508 /* dCoolmHeatdT is the actual number used for getting dT 00509 * thermal.dCooldT is just cooling */ 00510 } 00511 00512 if( trace.lgTrace ) 00513 { 00514 fprintf( ioQQQ, 00515 " ConvTempEdenIoniz TE:%.5e dT%.3e HTOT:%.5e CTOT:%.5e dCoolmHeatdT:%c%.4e C-H%.4e HDEN%.2e kase %i\n", 00516 phycon.te, 00517 thermal.dTemper, 00518 thermal.htot, 00519 thermal.ctot, 00520 chDEType, 00521 dCoolmHeatdT, 00522 CoolMHeat, 00523 dense.gas_phase[ipHYDROGEN], 00524 kase ); 00525 } 00526 00527 /* increment loop counter */ 00528 ++nTemperature_loop; 00529 } 00530 00531 /* fall through, exceeded limit to number of solutions, 00532 * >>chng 01 mar 14, or no electron density convergence 00533 * no temperature convergence */ 00534 if( fabs(CoolMHeat/thermal.htot ) > conv.HeatCoolRelErrorAllowed ) 00535 { 00536 conv.lgConvTemp = false; 00537 if( trace.nTrConvg>=2 ) 00538 { 00539 fprintf( ioQQQ, " ConvTempEdenIoniz7 fell through loop, heating cooling not converged, setting conv.lgConvTemp = false\n"); 00540 } 00541 } 00542 else 00543 { 00544 /* possible that ionization has failed, and we fell through to here 00545 * with a valid thermal solution */ 00546 conv.lgConvTemp = true; 00547 if( trace.nTrConvg>=2 ) 00548 { 00549 fprintf( ioQQQ, " ConvTempEdenIoniz8 fell through loop, heating cooling is converged (ion not?), setting conv.lgConvTemp = true\n"); 00550 } 00551 } 00552 } 00553 else 00554 { 00555 fprintf( ioQQQ, "ConvTempEdenIoniz finds insane solver%s \n" , conv.chSolverTemp ); 00556 ShowMe(); 00557 cdEXIT(EXIT_FAILURE); 00558 } 00559 return 0; 00560 } 00561 00562 /*MakeDeriv derive numerical derivative of heating minus cooling */ 00563 STATIC void MakeDeriv( 00564 const char *job, 00565 double *DerivNumer) 00566 { 00567 static long int nstore; 00568 static double OldTe , OldCool , OldHeat; 00569 00570 DEBUG_ENTRY( "MakeDeriv()" ); 00571 00572 if( strcmp(job,"zero") == 0 ) 00573 { 00574 nstore = 0; 00575 OldTe = 0.; 00576 OldCool = 0.; 00577 OldHeat = 0.; 00578 } 00579 00580 else if( strcmp(job,"incr") == 0 ) 00581 { 00582 00583 if( nstore > 0 && !thermal.lgTemperatureConstant && fabs(phycon.te - OldTe)>SMALLFLOAT ) 00584 { 00585 /* get numerical derivatives */ 00586 double dCdT = (thermal.ctot - OldCool)/(phycon.te - OldTe); 00587 double dHdT = (thermal.htot - OldHeat)/(phycon.te - OldTe); 00588 *DerivNumer = dCdT - dHdT; 00589 # if 0 00590 fprintf(ioQQQ,"derivderiv\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 00591 nzone, dCdT , thermal.dCooldT , 00592 thermal.ctot , OldCool,phycon.te , OldTe 00593 /*(thermal.ctot - OldCool),(phycon.te - OldTe),dHdT, thermal.dHeatdT*/); 00594 # endif 00595 } 00596 00597 else 00598 { 00599 /* if early in soln, or constant temperature, return zero */ 00600 *DerivNumer = 0.; 00601 } 00602 00603 OldTe = phycon.te; 00604 OldCool = thermal.ctot; 00605 OldHeat = thermal.htot; 00606 ++ nstore; 00607 } 00608 00609 else 00610 { 00611 fprintf( ioQQQ, " MakeDeriv called with insane job=%4.4s\n", 00612 job ); 00613 cdEXIT(EXIT_FAILURE); 00614 } 00615 return; 00616 } 00617 00618 /*PutHetCol save heating, cooling, and temperature in stack for numerical derivatives */ 00619 STATIC void PutHetCol(double te, 00620 double htot, 00621 double ctot) 00622 { 00623 00624 DEBUG_ENTRY( "PutHetCol()" ); 00625 00626 if( nzone == 0 ) 00627 { 00628 /* code is searching for first temp now - do not remember these numbers */ 00629 thermal.ipGrid = 0; 00630 } 00631 else 00632 { 00633 if( thermal.ipGrid >= NGRID ) 00634 thermal.ipGrid = 0; 00635 00636 ASSERT( thermal.ipGrid >= 0 ); 00637 ASSERT( thermal.ipGrid < NGRID ); 00638 ASSERT( nzone>=0 ); 00639 00640 thermal.TeGrid[thermal.ipGrid] = (realnum)te; 00641 thermal.HtGrid[thermal.ipGrid] = (realnum)htot; 00642 thermal.ClGrid[thermal.ipGrid] = (realnum)ctot; 00643 thermal.nZonGrid[thermal.ipGrid] = nzone; 00644 00645 ++thermal.ipGrid; 00646 } 00647 return; 00648 } 00649 00650 /*CoolHeatError evaluate ionization, and difference in heating and cooling, for temperature temp */ 00651 STATIC double CoolHeatError( double temp ) 00652 { 00653 DEBUG_ENTRY( "CoolHeatError()" ); 00654 00655 TempChange(temp , false); 00656 00657 /* converge the ionization and electron density; 00658 * this calls ionize until lgIonDone is true */ 00659 /* NB should NOT set insanity - but rather return error condition */ 00660 if( ConvEdenIoniz() ) 00661 lgAbort = true; 00662 00663 /* >>chng 01 mar 16, evaluate pressure here since changing and other values needed */ 00664 /* reevaluate pressure */ 00665 /* this sets values of pressure.PresTotlCurr */ 00666 PresTotCurrent(); 00667 00668 /* remember this temp, heating, and cooling */ 00669 PutHetCol(phycon.te,thermal.htot,thermal.ctot); 00670 return thermal.ctot - thermal.htot; 00671 } 00672 00673 /*#define EPS 3.e-8*/ 00674 #define ITERMAX 100 00675 /* use brent's method to find heating-cooling match */ 00676 STATIC double TeBrent( 00677 double x1, 00678 double x2) 00679 { 00680 long int iter; 00681 double c, 00682 d, 00683 e, 00684 fx1, 00685 fx2, 00686 fc, 00687 p, 00688 q, 00689 r, 00690 s, 00691 xm; 00692 00693 DEBUG_ENTRY( "TeBrent()" ); 00694 00695 /* first evaluate function at two boundaries, and check 00696 * that we have bracketed the target */ 00697 /*NB big logical error - this routine ignores return value of 00698 * routine it calls, so will continue despite disaster */ 00699 fx1 = CoolHeatError(x1); 00700 fx2 = CoolHeatError(x2); 00701 00702 if( fx1*fx2 > 0. ) 00703 { 00704 /* both values either positive or negative, this is insane */ 00705 fprintf( ioQQQ, " TeBrent called without proper bracket - this is impossible\n" ); 00706 fprintf( ioQQQ, " Sorry.\n" ); 00707 cdEXIT(EXIT_FAILURE); 00708 } 00709 00710 /* we have bracketed the target */ 00711 c = x2; 00712 fc = fx2; 00713 iter = 0; 00714 /* these are sanity checks, to get code past lint */ 00715 e = DBL_MAX; 00716 d = DBL_MAX; 00717 while( iter < ITERMAX ) 00718 { 00719 if( (fx2 > 0. && fc > 0.) || (fx2 < 0. && fc < 0.) ) 00720 { 00721 c = x1; 00722 fc = fx1; 00723 d = x2 - x1; 00724 e = d; 00725 } 00726 if( fabs(fc) < fabs(fx2) ) 00727 { 00728 x1 = x2; 00729 x2 = c; 00730 c = x1; 00731 fx1 = fx2; 00732 fx2 = fc; 00733 fc = fx1; 00734 } 00735 00736 xm = .5*(c - x2); 00737 00738 /* solution converged if residual less than tolerance, or hit zero */ 00739 if( fabs(xm) <= thermal.htot*conv.HeatCoolRelErrorAllowed || fx2 == 0. ) 00740 break; 00741 00742 if( fabs(e) >= thermal.htot*conv.HeatCoolRelErrorAllowed && fabs(fx1) > fabs(fx2) ) 00743 { 00744 double aa , bb; 00745 s = fx2/fx1; 00746 if( fp_equal( x1, c ) ) 00747 { 00748 p = 2.*xm*s; 00749 q = 1. - s; 00750 } 00751 else 00752 { 00753 q = fx1/fc; 00754 r = fx2/fc; 00755 p = s*(2.*xm*q*(q - r) - (x2 - x1)*(r - 1.)); 00756 q = (q - 1.)*(r - 1.)*(s - 1.); 00757 } 00758 if( p > 0. ) 00759 q = -q; 00760 00761 p = fabs(p); 00762 aa = fabs(thermal.htot*conv.HeatCoolRelErrorAllowed*q); 00763 bb = fabs(e*q); 00764 if( 2.*p < MIN2(3.*xm*q-aa,bb) ) 00765 { 00766 e = d; 00767 d = p/q; 00768 } 00769 else 00770 { 00771 d = xm; 00772 e = d; 00773 } 00774 } 00775 else 00776 { 00777 d = xm; 00778 e = d; 00779 } 00780 x1 = x2; 00781 fx1 = fx2; 00782 if( fabs(d) > thermal.htot*conv.HeatCoolRelErrorAllowed ) 00783 { 00784 x2 += d; 00785 } 00786 else 00787 { 00788 x2 += sign(thermal.htot*conv.HeatCoolRelErrorAllowed,xm); 00789 } 00790 fx2 = CoolHeatError(x2); 00791 00792 ++iter; 00793 } 00794 00795 /* check whether we fell through (hit limit to number of iterations) 00796 * of broke out when complete */ 00797 if( iter == ITERMAX ) 00798 { 00799 /* both values either positive or negative, this is insane */ 00800 fprintf( ioQQQ, " TeBrent did not converge in %i iterations\n",ITERMAX ); 00801 fprintf( ioQQQ, " Sorry.\n" ); 00802 cdEXIT(EXIT_FAILURE); 00803 } 00804 return( x2 ); 00805 00806 } 00807 00808 /* returns true if heating-cooling is converged */ 00809 bool lgConvTemp(void) 00810 { 00811 bool lgRet; 00812 DEBUG_ENTRY( "lgConvTemp()" ); 00813 00814 /* >>chng 03 sep 19, add check on phycon.TEMP_LIMIT_LOW */ 00815 if( fabs(thermal.htot - thermal.ctot)/thermal.htot <= conv.HeatCoolRelErrorAllowed || 00816 thermal.lgTemperatureConstant || phycon.te <= phycon.TEMP_LIMIT_LOW ) 00817 { 00818 /* announce that temp is converged if relative heating - cooling mismatch 00819 * is less than the relative heating cooling error allowed, or 00820 * if this is a constant temperature model */ 00821 lgRet = true; 00822 conv.lgConvTemp = true; 00823 } 00824 else 00825 { 00826 /* big mismatch, this has not converged */ 00827 lgRet = false; 00828 conv.lgConvTemp = false; 00829 } 00830 return lgRet; 00831 }