39 static double HeatOld=-1. , CoolOld=-1.;
48 double HeatChange =
thermal.htot - HeatOld;
49 double CoolChange =
thermal.ctot - CoolOld;
56 double HeatRelChange = fabs(HeatChange)/
SDIV( HeatCoolMax );
57 double CoolRelChange = fabs(CoolChange)/
SDIV( HeatCoolMax );
58 bool lgConverged =
true;
59 if(
MAX2( HeatRelChange , CoolRelChange ) >
conv.HeatCoolRelErrorAllowed )
143 if (
mole.species[i].location == NULL && !
mole_global.list[i]->isMonatomic())
145 for( molecule::nAtomsMap::iterator atom=
mole_global.list[i]->nAtom.begin();
148 long nelem = atom->first->el->Z-1;
149 if(
dense.lgElmtOn[nelem] )
176 enum {DEBUG_LOC=
false};
179 fprintf(
ioQQQ,
"DEBUG fraction molecular %.2e species %li O ices %.2e\n" ,
251 bool lgConvergeCoolHeat;
265 conv.nPres2Ioniz = 0;
269 conv.nTotalIoniz = 0;
270 conv.hist_pres_nzone = -1;
271 conv.hist_temp_nzone = -1;
273 conv.lgOscilOTS =
false;
276 dense.lgEdenBad =
false;
280 conv.BigEdenError = 0.;
281 conv.AverEdenError = 0.;
282 conv.BigHeatCoolError = 0.;
283 conv.AverHeatCoolError = 0.;
284 conv.BigPressError = 0.;
285 conv.AverPressError = 0.;
303 fprintf(
ioQQQ,
"\nConvInitSolution entered \n" );
315 fprintf(
ioQQQ,
" ConvInitSolution called, ITER=%2ld resetting Te to %10.2e\n",
321 fprintf(
ioQQQ,
" search set true\n" );
327 conv.lgSearch =
true;
328 conv.lgFirstSweepThisZone =
true;
329 conv.lgLastSweepThisZone =
false;
347 fprintf(
ioQQQ,
" ========================================================================================================================\n");
348 fprintf(
ioQQQ,
" ConvInitSolution: search set false 2 Te=%.3e========================================================================================\n" ,
phycon.te);
349 fprintf(
ioQQQ,
" ========================================================================================================================\n");
351 conv.lgSearch =
false;
375 const bool lgDoInitConv =
true;
377 conv.lgSearch =
true;
378 conv.lgFirstSweepThisZone =
true;
379 conv.lgLastSweepThisZone =
false;
383 fprintf(
ioQQQ,
" ConvInitSolution called, new temperature.\n" );
387 TeBoundLow =
phycon.TEMP_LIMIT_LOW/1.001;
397 TeBoundHigh =
phycon.TEMP_LIMIT_HIGH/1.2;
410 TeFirst =
MAX2( 4000. , TeFirst );
419 TeFirst =
MIN2( 1e6 , TeBoundHigh );
424 TeFirst =
MAX2( 4000., TeBoundLow );
431 if( !
thermal.lgTemperatureConstant )
432 TeFirst =
max( TeFirst ,
phycon.TEnerDen );
436 fprintf(
ioQQQ,
"ConvInitSolution doing initial solution with temp=%.2e\n",
450 double CoolNet=0, dCoolNetDT=0;
452 lgConvergeCoolHeat =
false;
453 for( ionlup=0; ionlup<2; ++ionlup )
456 fprintf(
ioQQQ,
"ConvInitSolution calling lgCoolNetConverge "
457 "initial setup %li with Te=%.3e C=%.2e H=%.2e d(C-H)/dT %.3e\n",
470 if( !lgConvergeCoolHeat )
471 fprintf(
ioQQQ,
" PROBLEM ConvInitSolution did not converge the heating-cooling during initial search.\n");
476 fprintf(
ioQQQ,
" DISASTER PROBLEM ConvInitSolution: Search for "
477 "initial conditions aborted - lgAbort set true.\n" );
487 bool lgConvergedLoop =
false;
488 long int LoopThermal = 0;
492 const long int LIMIT_THERMAL_LOOP=300;
493 double CoolMHeatSave=-1. , TempSave=-1. , TeNew=-1.,
CoolSave=-1;
494 while( !lgConvergedLoop && LoopThermal < LIMIT_THERMAL_LOOP )
497 CoolMHeatSave = CoolNet;
504 ASSERT( TeChangeFactor>0. && TeChangeFactor< 1. );
506 TeNew =
phycon.te - CoolNet / dCoolNetDT;
508 TeNew =
MAX2(
phycon.te*TeChangeFactor , TeNew );
509 TeNew =
MIN2(
phycon.te/TeChangeFactor , TeNew );
515 fprintf(
ioQQQ,
"ConvInitSolution %4li TeNnew=%.3e "
516 "TeChangeFactor=%.3e Cnet/H=%.2e dCnet/dT=%10.2e Ne=%.3e"
517 " Ograins %.2e FracMoleMax %.2e\n",
518 LoopThermal , TeNew , TeChangeFactor ,
531 lgConvergedLoop =
true;
532 else if( (CoolNet/fabs(CoolNet))*(CoolMHeatSave/fabs(CoolMHeatSave)) <= 0)
534 lgConvergedLoop =
true;
535 else if(
phycon.te <= TeBoundLow ||
phycon.te>=TeBoundHigh)
536 lgConvergedLoop =
true;
540 if( !lgConvergedLoop )
541 fprintf(
ioQQQ,
"PROBLEM ConvInitSolution: temperature solution not "
542 "found in initial search, final Te=%.2e\n",
552 fprintf(
ioQQQ,
" Change in sign of C-H found, Te brackets %.3e "
553 "%.3e Cool brackets %.3e %.3e ",
554 TempSave ,
phycon.te , CoolMHeatSave , CoolNet);
566 double Alog = log(
thermal.ctot ) - alpha * log(
phycon.te );
568 double TeLn = (log(
thermal.htot ) - Alog) / alpha ;
571 if( TeLn < log(
MIN2(
phycon.te , TempSave )) )
573 else if( TeLn > log(
MAX2(
phycon.te , TempSave )) )
576 TeNew = pow(
EE , TeLn );
583 fprintf(
ioQQQ,
" interpolating to Te %.3e \n\n",
598 conv.lgSearch =
false;
614 fprintf(
ioQQQ,
"\nConvInitSolution: end 1st iteration search phase "
615 "finding Te=%.3e\n\n" ,
phycon.te);
620 fprintf(
ioQQQ,
" ConvInitSolution return, TE:%10.2e==================\n",
639 double PresNew =
pressure.PresTotlCurr;
642 if(
pressure.lgPressureInitialSpecified )
644 PresNew =
pressure.PressureInitialSpecified;
648 " PresTotCurrent 1st zn old values of PresTotlInit:%.3e"
657 if(
dynamics.lgTimeDependentStatic )
660 static double PresTotalInitTime;
663 PresTotalInitTime =
pressure.PresTotlInit;
665 else if (
dense.lgPressureVaryTime)
667 pressure.PresTotlInit = PresTotalInitTime *
668 pow( 1.+(
dynamics.time_elapsed/
dense.PressureVaryTimeTimescale) ,
669 dense.PressureVaryTimeIndex);
683 conv.nPres2Ioniz = 0;
685 dense.lgEdenBad =
false;
690 conv.nTotalIoniz_start =
conv.nTotalIoniz;
694 conv.BigEdenError = 0.;
695 conv.AverEdenError = 0.;
696 conv.BigHeatCoolError = 0.;
697 conv.AverHeatCoolError = 0.;
698 conv.BigPressError = 0.;
699 conv.AverPressError = 0.;
733 for( nelem=2; nelem <
LIMELM; nelem++ )
736 for( i=0; i < nelem; i++ )
738 if(
dense.xIonDense[nelem][i+1] > 0. )
743 nelem_reflux = nelem;
752 if( nflux_old !=
rfield.nflux )
774 fprintf(
ioQQQ,
" nflux updated from %li to %li, anu from %g to %g \n",
775 nflux_old ,
rfield.nflux ,
777 fprintf(
ioQQQ,
" caused by element %li ion %li \n",
778 nelem_reflux ,ion_reflux );
784 if( strcmp(
dense.chDenseLaw,
"DYNA") == 0 )
789 static const double PCHNG = 0.98;