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