00001
00002
00003
00004
00005
00006 #include "cddefines.h"
00007 #include "physconst.h"
00008 #include "trace.h"
00009 #include "struc.h"
00010 #include "rfield.h"
00011 #include "mole.h"
00012 #include "dense.h"
00013 #include "stopcalc.h"
00014 #include "heavy.h"
00015 #include "wind.h"
00016 #include "geometry.h"
00017 #include "thermal.h"
00018 #include "radius.h"
00019 #include "phycon.h"
00020 #include "pressure.h"
00021 #include "conv.h"
00022 #include "hmi.h"
00023 #include "dynamics.h"
00024
00025
00026 static double dCoolNetDTOld = 0;
00027
00028 static double OxyInGrains , FracMoleMax;
00029
00030
00031
00032
00033
00034
00035 STATIC bool lgCoolHeatCheckConverge(
00036 double *CoolNet,
00037 bool lgReset )
00038 {
00039 static double HeatOld=-1. , CoolOld=-1.;
00040 DEBUG_ENTRY( "lgCoolHeatCheckConverge()" );
00041
00042 if( lgReset )
00043 {
00044 HeatOld = -1;
00045 CoolOld = -1;
00046 }
00047
00048 double HeatChange = thermal.htot - HeatOld;
00049 double CoolChange = thermal.ctot - CoolOld;
00050
00051
00052 double HeatCoolMax = MAX2( thermal.htot , thermal.ctot );
00053
00054
00055
00056 double HeatRelChange = fabs(HeatChange)/SDIV( HeatCoolMax );
00057 double CoolRelChange = fabs(CoolChange)/SDIV( HeatCoolMax );
00058 bool lgConverged = true;
00059 if( MAX2( HeatRelChange , CoolRelChange ) > conv.HeatCoolRelErrorAllowed )
00060 lgConverged = false;
00061
00062 *CoolNet = thermal.ctot - thermal.htot;
00063
00064 HeatOld = thermal.htot;
00065 CoolOld = thermal.ctot;
00066 return lgConverged;
00067 }
00068
00069
00070 STATIC bool lgCoolNetConverge(
00071 double *CoolNet ,
00072 double *dCoolNetDT,
00073 bool lgReset )
00074 {
00075 const long int LOOP_MAX=20;
00076 long int loop = 0;
00077 bool lgConvergeCoolHeat = false;
00078
00079 DEBUG_ENTRY( "lgCoolNetConverge()" );
00080
00081 while( loop < LOOP_MAX && !lgConvergeCoolHeat && !lgAbort )
00082 {
00083 if( ConvEdenIoniz() )
00084 {
00085
00086 lgAbort = true;
00087 }
00088 lgConvergeCoolHeat = lgCoolHeatCheckConverge( CoolNet, lgReset && loop==0 );
00089 if( trace.lgTrace || trace.nTrConvg>=3 )
00090 fprintf(ioQQQ," lgCoolNetConverge %li calls to ConvEdenIoniz, converged? %c\n",
00091 loop , TorF( lgConvergeCoolHeat ) );
00092 ++loop;
00093 }
00094
00095 if( thermal.ConstTemp > 0 )
00096 {
00097
00098
00099 *CoolNet = phycon.te - thermal.ConstTemp;
00100 *dCoolNetDT = 1.f;
00101 }
00102 else if( phycon.te < 4000.f )
00103 {
00104
00105
00106
00107
00108
00109 *dCoolNetDT = thermal.ctot / phycon.te;
00110 }
00111 else if( thermal.htot / thermal.ctot > 2.f )
00112 {
00113
00114
00115
00116 *dCoolNetDT = thermal.ctot / phycon.te;
00117 }
00118 else
00119 {
00120 *dCoolNetDT = thermal.dCooldT - thermal.dHeatdT;
00121 if( dCoolNetDTOld * *dCoolNetDT < 0. )
00122 {
00123
00124 fprintf(ioQQQ,"PROBLEM CoolNet derivative oscillating in initial solution "\
00125 "Te=%.3e dCoolNetDT=%.3e CoolNet=%.3e dCooldT=%.3e dHeatdT=%.3e\n",
00126 phycon.te , *dCoolNetDT , *CoolNet,
00127 thermal.dCooldT , thermal.dHeatdT);
00128 *dCoolNetDT = *CoolNet / phycon.te;
00129 }
00130 }
00131 return lgConvergeCoolHeat;
00132 }
00133
00134
00135 STATIC void ChemImportance( void )
00136 {
00137 long int i;
00138
00139 DEBUG_ENTRY( "ChemImportance()" );
00140
00141 FracMoleMax = 0.;
00142 long int ipMax = -1;
00143 for( i=0; i<mole.num_comole_calc; ++i )
00144 {
00145
00146 if( dense.lgElmtOn[COmole[i]->nelem_hevmol] && COmole[i]->n_nuclei>1 )
00147 {
00148 realnum frac = COmole[i]->hevmol / SDIV(dense.gas_phase[COmole[i]->nelem_hevmol]);
00149 frac *= (realnum) COmole[i]->nElem[COmole[i]->nelem_hevmol];
00150 if( frac > FracMoleMax )
00151 {
00152
00153 FracMoleMax = frac;
00154 ipMax = i;
00155 }
00156 }
00157 }
00158
00159
00160 OxyInGrains = 0;
00161 for(i=0;i<mole.num_comole_calc;++i)
00162 {
00163 OxyInGrains += (1 - COmole[i]->lgGas_Phase)*COmole[i]->hevmol*COmole[i]->nElem[ipOXYGEN];
00164 }
00165
00166 OxyInGrains /= SDIV(dense.gas_phase[ipOXYGEN]);
00167
00168 {
00169
00170 enum {DEBUG_LOC=false};
00171 if( DEBUG_LOC )
00172 {
00173 fprintf(ioQQQ,"DEBUG fraction molecular %.2e species %li O ices %.2e\n" ,
00174 FracMoleMax , ipMax ,OxyInGrains );
00175 }
00176 }
00177
00178 return;
00179 }
00180
00181 double FindTempChangeFactor(void)
00182 {
00183 double TeChangeFactor;
00184
00185 DEBUG_ENTRY( "FindTempChangeFactor()" );
00186
00187
00188
00189
00190
00191
00192
00193 ChemImportance();
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205 if( OxyInGrains > 0.99 )
00206 {
00207 TeChangeFactor = 0.999;
00208 }
00209 else if( OxyInGrains > 0.9 )
00210 {
00211 TeChangeFactor = 0.99;
00212 }
00213
00214 else if( phycon.te < 5. ||
00215
00216 OxyInGrains > 0.1 )
00217 {
00218 TeChangeFactor = 0.97;
00219 }
00220
00221 else if( phycon.te < 20. || FracMoleMax > 0.1 )
00222 {
00223
00224
00225 TeChangeFactor = 0.95;
00226 }
00227 else
00228 {
00229
00230 TeChangeFactor = 0.8;
00231 }
00232 return TeChangeFactor;
00233 }
00234
00235
00236 int ConvInitSolution(void)
00237 {
00238 long int i,
00239 ionlup,
00240 nelem ,
00241 nflux_old,
00242 nelem_reflux ,
00243 ion_reflux;
00244
00245 bool lgConvergeCoolHeat;
00246
00247 double
00248 TeChangeFactor,
00249 TeBoundLow,
00250 TeBoundHigh;
00251
00252 DEBUG_ENTRY( "ConvInitSolution()" );
00253
00254
00255 set_NaN( TeBoundLow );
00256 set_NaN( TeBoundHigh );
00257
00258
00259 conv.nPres2Ioniz = 0;
00260
00261
00262
00263 conv.nTotalIoniz = 0;
00264 conv.hist_pres_nzone = -1;
00265 conv.hist_temp_nzone = -1;
00266
00267 conv.lgOscilOTS = false;
00268
00269 lgAbort = false;
00270 dense.lgEdenBad = false;
00271 dense.nzEdenBad = 0;
00272
00273
00274 conv.BigEdenError = 0.;
00275 conv.AverEdenError = 0.;
00276 conv.BigHeatCoolError = 0.;
00277 conv.AverHeatCoolError = 0.;
00278 conv.BigPressError = 0.;
00279 conv.AverPressError = 0.;
00280
00281
00282 if( fp_equal( radius.sdrmin, radius.sdrmax ) )
00283 {
00284
00285 double rfac = radius.lgSdrmaxRel ? radius.Radius : 1.;
00286 radius.drad = MIN2( rfac*radius.sdrmax, radius.drad );
00287 radius.drad_mid_zone = radius.drad/2.;
00288 radius.drad_x_fillfac = radius.drad * geometry.FillFac;
00289 }
00290
00291
00292
00293
00294 ASSERT( dense.xIonDense[ipHYDROGEN][0] >=0 && dense.xIonDense[ipHYDROGEN][1]>= 0.);
00295
00296 if( trace.nTrConvg )
00297 fprintf( ioQQQ, "\nConvInitSolution entered \n" );
00298
00299
00300
00301
00302
00303
00304 if( iteration != 1 )
00305 {
00306
00307 if( trace.lgTrace || trace.nTrConvg )
00308 {
00309 fprintf( ioQQQ, " ConvInitSolution called, ITER=%2ld resetting Te to %10.2e\n",
00310 iteration, struc.testr[0] );
00311 }
00312
00313 if( trace.lgTrace || trace.nTrConvg )
00314 {
00315 fprintf( ioQQQ, " search set true\n" );
00316 }
00317
00318
00319
00320
00321 conv.lgSearch = true;
00322 conv.lgFirstSweepThisZone = true;
00323 conv.lgLastSweepThisZone = false;
00324
00325
00326 TempChange( struc.testr[0] , false);
00327
00328
00329 PresTotCurrent();
00330
00331
00332 if( ConvPresTempEdenIoniz() )
00333 {
00334
00335 lgAbort = true;
00336 return(1);
00337 }
00338
00339 if( trace.lgTrace || trace.nTrConvg )
00340 {
00341 fprintf( ioQQQ, " ========================================================================================================================\n");
00342 fprintf( ioQQQ, " ConvInitSolution: search set false 2 Te=%.3e========================================================================================\n" , phycon.te);
00343 fprintf( ioQQQ, " ========================================================================================================================\n");
00344 }
00345 conv.lgSearch = false;
00346
00347
00348 if( ConvTempEdenIoniz() )
00349 {
00350
00351 lgAbort = true;
00352 return(1);
00353 }
00354
00355
00356
00357 PresTotCurrent();
00358 }
00359
00360 else
00361 {
00362
00363
00364
00365
00366
00367
00368
00369 conv.lgSearch = true;
00370 conv.lgFirstSweepThisZone = true;
00371 conv.lgLastSweepThisZone = false;
00372
00373 if( trace.lgTrace )
00374 {
00375 fprintf( ioQQQ, " ConvInitSolution called, new temperature.\n" );
00376 }
00377
00378
00379 TeBoundLow = phycon.TEMP_LIMIT_LOW/1.001;
00380
00381
00382
00383
00384 TeBoundLow = MAX2( TeBoundLow , StopCalc.TeFloor );
00385
00386
00387
00388
00389 TeBoundHigh = phycon.TEMP_LIMIT_HIGH/1.2;
00390
00391
00392
00393 double TeFirst;
00394 if( thermal.ConstTemp > 0 )
00395 {
00396
00397
00398
00399 TeFirst = MAX2( 4000. , thermal.ConstTemp );
00400
00401
00402
00403 if( hmi.lgNoH2Mole )
00404 TeFirst = thermal.ConstTemp;
00405 }
00406 else if( thermal.lgTeHigh )
00407 {
00408
00409 TeFirst = MIN2( 1e6 , TeBoundHigh );
00410 }
00411 else
00412 {
00413
00414 TeFirst = MAX2( 4000., TeBoundLow );
00415 }
00416
00417
00418
00419
00420
00421 if( !thermal.lgTemperatureConstant )
00422 TeFirst = max( TeFirst , phycon.TEnerDen );
00423
00424 TempChange(TeFirst , false);
00425 if( trace.lgTrace || trace.nTrConvg>=2 )
00426 fprintf(ioQQQ,"ConvInitSolution doing initial solution with temp=%.2e\n",
00427 phycon.te);
00428
00429
00430 PresTotCurrent();
00431
00432 thermal.htot = 1.;
00433 thermal.ctot = 1.;
00434
00435
00436
00437
00438 double CoolNet=0, dCoolNetDT=0;
00439
00440 lgConvergeCoolHeat = false;
00441 for( ionlup=0; ionlup<2; ++ionlup )
00442 {
00443 if( trace.lgTrace || trace.nTrConvg>=2 )
00444 fprintf( ioQQQ, "ConvInitSolution calling lgCoolNetConverge "
00445 "initial setup %li with Te=%.3e C=%.2e H=%.2e d(C-H)/dT %.3e\n",
00446 ionlup,phycon.te , thermal.ctot , thermal.htot,dCoolNetDT );
00447
00448 dCoolNetDTOld = 0.;
00449 lgConvergeCoolHeat = lgCoolNetConverge( &CoolNet , &dCoolNetDT, true );
00450 if( lgAbort )
00451 break;
00452
00453
00454
00455
00456 radius_first();
00457 }
00458 if( !lgConvergeCoolHeat )
00459 fprintf(ioQQQ," PROBLEM ConvInitSolution did not converge the heating-cooling during initial search.\n");
00460
00461 if( lgAbort )
00462 {
00463
00464 fprintf( ioQQQ, " DISASTER PROBLEM ConvInitSolution: Search for "
00465 "initial conditions aborted - lgAbort set true.\n" );
00466 ShowMe();
00467
00468 return(1);
00469 }
00470
00471
00472
00473 if( thermal.ConstTemp<=0 )
00474 dCoolNetDT = thermal.ctot / phycon.te;
00475 bool lgConvergedLoop = false;
00476 long int LoopThermal = 0;
00477
00478
00479
00480 const long int LIMIT_THERMAL_LOOP=300;
00481 double CoolMHeatSave=-1. , TempSave=-1. , TeNew=-1.,CoolSave=-1;
00482 while( !lgConvergedLoop && LoopThermal < LIMIT_THERMAL_LOOP )
00483 {
00484
00485 CoolMHeatSave = CoolNet;
00486 dCoolNetDTOld = dCoolNetDT;
00487 CoolSave = thermal.ctot;
00488 TempSave = phycon.te;
00489
00490
00491 TeChangeFactor = FindTempChangeFactor();
00492 ASSERT( TeChangeFactor>0. && TeChangeFactor< 1. );
00493
00494 TeNew = phycon.te - CoolNet / dCoolNetDT;
00495
00496 TeNew = MAX2( phycon.te*TeChangeFactor , TeNew );
00497 TeNew = MIN2( phycon.te/TeChangeFactor , TeNew );
00498
00499 TempChange(TeNew , true);
00500 lgConvergeCoolHeat = lgCoolNetConverge( &CoolNet , & dCoolNetDT, false );
00501
00502 if( trace.lgTrace || trace.nTrConvg>=2 )
00503 fprintf(ioQQQ,"ConvInitSolution %4li TeNnew=%.3e "
00504 "TeChangeFactor=%.3e Cnet/H=%.2e dCnet/dT=%10.2e Ne=%.3e"
00505 " Ograins %.2e FracMoleMax %.2e\n",
00506 LoopThermal , TeNew , TeChangeFactor ,
00507 CoolNet/SDIV(thermal.htot) , dCoolNetDT,
00508 dense.eden , OxyInGrains , FracMoleMax );
00509
00510 if( lgAbort )
00511 return(1);
00512
00513
00514
00515
00516
00517 if( fabs(CoolNet)< SMALLDOUBLE )
00518
00519 lgConvergedLoop = true;
00520 else if( (CoolNet/fabs(CoolNet))*(CoolMHeatSave/fabs(CoolMHeatSave)) <= 0)
00521
00522 lgConvergedLoop = true;
00523 else if( fabs(CoolNet)/MAX2( thermal.ctot , thermal.htot ) <conv.HeatCoolRelErrorAllowed*10. )
00524 lgConvergedLoop = true;
00525 else if( phycon.te <= TeBoundLow || phycon.te>=TeBoundHigh)
00526 lgConvergedLoop = true;
00527 ++LoopThermal;
00528 }
00529
00530 if( !lgConvergedLoop )
00531 fprintf(ioQQQ,"PROBLEM ConvInitSolution: temperature solution not "
00532 "found in initial search, final Te=%.2e\n",
00533 phycon.te );
00534
00535
00536
00537
00538 if( thermal.ConstTemp<=0 &&
00539 ! fp_equal( TempSave , phycon.te ) )
00540 {
00541 if( trace.lgTrace || trace.nTrConvg>=2 )
00542 fprintf(ioQQQ," Change in sign of C-H found, Te brackets %.3e "
00543 "%.3e Cool brackets %.3e %.3e ",
00544 TempSave , phycon.te , CoolMHeatSave , CoolNet);
00545
00546
00547 double alpha = log(thermal.ctot/CoolSave) / log( phycon.te/TempSave);
00548 if( fabs(alpha) < SMALLFLOAT )
00549
00550 TeNew = phycon.te;
00551 else
00552 {
00553
00554 if( thermal.ctot<0. || thermal.htot<0. )
00555 TotalInsanity();
00556 double Alog = log( thermal.ctot ) - alpha * log( phycon.te );
00557
00558 double TeLn = (log( thermal.htot ) - Alog) / alpha ;
00559
00560
00561 if( TeLn < log( MIN2(phycon.te , TempSave )) )
00562 TeNew = MIN2(phycon.te , TempSave );
00563 else if( TeLn > log( MAX2(phycon.te , TempSave )) )
00564 TeNew = MAX2(phycon.te , TempSave );
00565 else
00566 TeNew = pow( EE , TeLn );
00567 }
00568
00569 ASSERT( TeNew>=MIN2 ( TempSave , phycon.te ) ||
00570 TeNew<=MAX2 ( TempSave , phycon.te ));
00571
00572 if( trace.lgTrace || trace.nTrConvg>=2 )
00573 fprintf(ioQQQ," interpolating to Te %.3e \n\n",
00574 TeNew);
00575
00576
00577 TempChange(TeNew , true);
00578 }
00579
00580 if( ConvTempEdenIoniz() )
00581 {
00582
00583 lgAbort = true;
00584 return(1);
00585 }
00586
00587
00588 PresTotCurrent();
00589
00590 if( trace.lgTrace || trace.nTrConvg )
00591 {
00592 fprintf( ioQQQ, "\nConvInitSolution: end 1st iteration search phase "
00593 "finding Te=%.3e\n\n" , phycon.te);
00594 }
00595 conv.lgSearch = false;
00596
00597 if( trace.lgTrace )
00598 {
00599 fprintf( ioQQQ, " ConvInitSolution return, TE:%10.2e==================\n",
00600 phycon.te );
00601 }
00602 }
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616 if( iteration==1 || dense.lgDenseInitConstant )
00617 {
00618 double PresNew = pressure.PresTotlCurr;
00619
00620
00621 if( pressure.lgPressureInitialSpecified )
00622
00623 PresNew = pressure.PressureInitialSpecified;
00624 if( trace.lgTrace )
00625 {
00626 fprintf( ioQQQ,
00627 " PresTotCurrent 1st zn old values of PresTotlInit:%.3e"
00628 " to %.3e \n",
00629 pressure.PresTotlInit,
00630 PresNew);
00631 }
00632
00633 pressure.PresTotlInit = PresNew;
00634 }
00635
00636 if( dynamics.lgTimeDependentStatic )
00637 {
00638
00639 static double PresTotalInitTime;
00640 if( iteration <= dynamics.n_initial_relax )
00641 {
00642 PresTotalInitTime = pressure.PresTotlInit;
00643 }
00644 else if (dense.lgPressureVaryTime)
00645 {
00646 pressure.PresTotlInit = PresTotalInitTime *
00647 pow( 1.+(dynamics.time_elapsed/dense.PressureVaryTimeTimescale) ,
00648 dense.PressureVaryTimeIndex);
00649 }
00650 fprintf(ioQQQ,"DEBUG conv_init_solution time dependent time=%.2e sets "
00651 "pressure to %.2e\n", dynamics.time_elapsed ,
00652 pressure.PresTotlInit);
00653 }
00654
00655
00656
00657 ConvPresTempEdenIoniz();
00658
00659
00660
00661
00662 conv.nPres2Ioniz = 0;
00663
00664 dense.lgEdenBad = false;
00665 dense.nzEdenBad = 0;
00666
00667
00668
00669 conv.nTotalIoniz_start = conv.nTotalIoniz;
00670
00671
00672
00673 conv.BigEdenError = 0.;
00674 conv.AverEdenError = 0.;
00675 conv.BigHeatCoolError = 0.;
00676 conv.AverHeatCoolError = 0.;
00677 conv.BigPressError = 0.;
00678 conv.AverPressError = 0.;
00679
00680
00681 if( iteration == 1 )
00682 {
00683
00684
00685 struc.testr[0] = (realnum)phycon.te;
00686
00687 struc.hden[0] = dense.gas_phase[ipHYDROGEN];
00688
00689 struc.DenParticles[0] = dense.pden;
00690 struc.heatstr[0] = thermal.htot;
00691 struc.coolstr[0] = thermal.ctot;
00692 struc.volstr[0] = (realnum)radius.dVeffAper;
00693 struc.drad_x_fillfac[0] = (realnum)radius.drad_x_fillfac;
00694 struc.histr[0] = dense.xIonDense[ipHYDROGEN][0];
00695 struc.hiistr[0] = dense.xIonDense[ipHYDROGEN][1];
00696 struc.ednstr[0] = (realnum)dense.eden;
00697 struc.o3str[0] = dense.xIonDense[ipOXYGEN][2];
00698 struc.DenMass[0] = dense.xMassDensity;
00699 struc.drad[0] = (realnum)radius.drad;
00700 }
00701
00702
00703
00704
00705
00706
00707
00708
00709 nflux_old = rfield.nflux;
00710 nelem_reflux = -1;
00711 ion_reflux = -1;
00712 for( nelem=2; nelem < LIMELM; nelem++ )
00713 {
00714
00715 for( i=0; i < nelem; i++ )
00716 {
00717 if( dense.xIonDense[nelem][i+1] > 0. )
00718 {
00719 if( Heavy.ipHeavy[nelem][i] > rfield.nflux )
00720 {
00721 rfield.nflux = Heavy.ipHeavy[nelem][i];
00722 nelem_reflux = nelem;
00723 ion_reflux = i;
00724 }
00725 }
00726 }
00727 }
00728
00729
00730
00731 if( nflux_old != rfield.nflux )
00732 {
00733
00734 rfield_opac_zero( nflux_old-1 , rfield.nflux );
00735
00736
00737
00738
00739
00740 PresTotCurrent();
00741
00742
00743 if( ConvBase(1) )
00744 {
00745
00746 lgAbort = true;
00747 return( 1 );
00748 }
00749
00750
00751 if( trace.lgTrace )
00752 {
00753 fprintf(ioQQQ," nflux updated from %li to %li, anu from %g to %g \n",
00754 nflux_old , rfield.nflux ,
00755 rfield.anu[nflux_old] , rfield.anu[rfield.nflux] );
00756 fprintf(ioQQQ," caused by element %li ion %li \n",
00757 nelem_reflux ,ion_reflux );
00758 }
00759 }
00760
00761
00762
00763 if( (strcmp(dense.chDenseLaw,"WIND") == 0) && !( wind.lgBallistic() || wind.lgStatic()) )
00764 {
00765
00766
00767
00768 # define PCHNG 0.98
00769
00770 if( ConvPresTempEdenIoniz() )
00771 {
00772
00773 lgAbort = true;
00774 return(1);
00775 }
00776
00777 pressure.PresTotlInit *= PCHNG;
00778 if( ConvPresTempEdenIoniz() )
00779 {
00780
00781 lgAbort = true;
00782 return(1);
00783 }
00784
00785 pressure.PresTotlInit /= PCHNG;
00786 if( ConvPresTempEdenIoniz() )
00787 {
00788
00789 lgAbort = true;
00790 return(1);
00791 }
00792
00793 # undef PCHNG
00794 }
00795
00796 return( lgAbort );
00797 }