43 static double MIN_CONV_FRAC=1e-4;
56 for (
long ion = 0; ion <= nelem+1; ++ion)
70 bool lgConverg_v =
false;
83 double abundold=0. , abundnew=0.;
106 for(
long ion=0; ion <= (nelem+1); ++ion )
109 if( OldFracs[nelem][ion]/Abund > MIN_CONV_FRAC &&
115 OldFracs[nelem][ion];
116 change =
MAX2(change, one );
118 if( change>bigchange )
121 abundold = OldFracs[nelem][ion]/Abund;
130 if( change >= delta )
135 ASSERT( abundold>0. && abundnew>0. );
143 for(
long ion=0; ion <= (nelem+1); ++ion )
153 fprintf(
ioQQQ,
" nz %ld loop %ld element %li converged? %c worst %ld change %g\n",
154 nzone, loop_ion, nelem,
TorF(lgConverg_v),ionchg,bigchange);
155 for(
long ion=0; ion<(nelem+1); ++ion )
165 double dground(
long nelem,
long ion)
168 if( ! ( nelem >=
ipHYDROGEN && nelem < LIMELM &&
169 ion >= 0 && ion <= nelem+1 ) )
172 "ERROR: Invalid ionization stage in dground()."
173 " nelem= %ld\t ion= %ld\n",
177 long ipISO = nelem - ion;
178 if (ipISO >= 0 && ipISO <
NISO)
179 return iso_sp[ipISO][nelem].
st[0].Pop();
199 static double SecondOld;
200 static long int nzoneOTS=-1;
201 const int LOOP_ION_LIMIT = 10;
203 static double SumOTS=0. , OldSumOTS[2]={0.,0.};
206 IonizConverg lgIonizConverg;
210 bool lgNewTrim =
false;
211 bool lgIonizWidened =
false;
215 static double wide_te_used = -1.;
216 static long int nZoneIonizWidenCalled = 0;
218 bool lgWiden = ( nZoneIonizWidenCalled !=
nzone ||
223 nZoneIonizWidenCalled =
nzone;
225 lgIonizWidened =
true;
239 # if !defined(NDEBUG)
296 for(
long nelem=ipISO; nelem<
LIMELM;++nelem )
301 save_iso_grnd[ipISO][nelem] =
iso_sp[ipISO][nelem].
st[0].Pop();
317 " ConvBase called. %.2f Te:%.3e HI:%.3e HII:%.3e H2:%.3e Ne:%.3e htot:%.3e CSUP:%.2e Conv?%c\n",
369 static long int nZoneIonizTrimCalled = 0;
372 nZoneIonizTrimCalled = 0;
380 bool lgIonizTrimCalled =
false;
400 lgIonizTrimCalled =
true;
411 nZoneIonizTrimCalled =
nzone;
415 # if !defined(NDEBUG)
469 (*diatom)->CalcPhotoionizationRate();
480 double xIonDense0[nconv][
LIMELM][LIMELM+1];
481 double xIonDense_prev[
LIMELM][LIMELM+1];
483 for (
long nelem=0; nelem <
LIMELM; ++nelem)
487 lgOscill[nelem] =
false;
534 for (
long ion = 0; ion <= nelem+1; ++ion)
536 xIonDense_prev[nelem][ion] = xIonDense0[0][nelem][ion];
542 bool lgShortCircuit =
false;
543 for( ion_loop=0; ion_loop<nconv && !lgShortCircuit; ++ion_loop)
551 for (
long ion = 0; ion <= nelem+1; ++ion)
560 for (
long nhi=0; nhi<
LIMELM; ++nhi)
561 netion[nlo][nhi] = 0.0;
575 long nlo =
MIN2(nelem, n2);
576 long nhi =
MAX2(nelem, n2);
577 if (nlo >= t_atmdat::NCX)
580 double dl1 = dground(nlo,0);
581 double ul1 = dground(nlo,1);
582 double dl2 = dground(nhi,0);
583 double ul2 = dground(nhi,1);
589 netion[nlo][nhi] += hion;
593 netion[nlo][nhi] -= hion;
607 for (
long nhi = nlo+1; nhi <
LIMELM; ++nhi)
611 double dl1 = dground(nlo,0);
612 double ul1 = dground(nlo,1);
613 double dl2 = dground(nhi,0);
614 double ul2 = dground(nhi,1);
620 bool lgCheckAll =
false;
622 fabs(netion[nlo][nhi]) > ion_cmp && ul1*ul2 > 1e-8*dl1*dl2 )
626 <<
" CX inconsistency";
634 bool lgCanShortCircuit = (ion_loop+1 < nconv);
635 for(
long nelem=
ipHYDROGEN; nelem<LIMELM && lgCanShortCircuit; ++nelem )
642 double x0 = xIonDense0[ion_loop][nelem][ion];
646 lgCanShortCircuit =
false;
651 lgShortCircuit = lgCanShortCircuit;
661 double tot0 = 0., tot1 = 0.;
662 double xIonNew[LIMELM+1];
663 for (
long ion = 0; ion <= nelem+1; ++ion)
665 double x0 = xIonDense0[nconv-2][nelem][ion];
666 double x1 = xIonDense0[nconv-1][nelem][ion];
674 double step0 = x1-
x0, step1 = x2-
x1, abs1 = fabs(step1);
676 if ( abs1 > 1000.0*((
double)DBL_EPSILON)*x2 )
678 double denom = fabs(step1-step0);
679 double sgn = (step1*step0 > 0)? 1.0 : -1.0;
682 const double MAXACC=100.0;
683 double extfac = 1.0/(denom/abs1 + 1.0/MAXACC);
684 double extstep = sgn*extfac*step1;
686 double predict = x2+extstep;
688 xIonNew[ion] = predict;
691 fprintf(
ioQQQ,
"Extrap %3ld %3ld %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g\n",
693 x0,x0-xIonDense0[nconv-3][nelem][ion],x1-x0,x2-x1,extstep,predict);
695 tot1 += xIonNew[ion];
699 double scal = tot0/tot1;
700 for (
long ion = 0; ion <= nelem+1; ++ion)
706 long ion = nelem-ipISO;
724 bool lgPostExtrapSolve =
true;
725 if (lgPostExtrapSolve)
744 for (
long ion = 0; ion <= nelem+1; ++ion)
750 long ion = nelem-ipISO;
781 vector<const molecule*>debug_list;
792 for (
long ion = 0; ion <= nelem+1; ++ion)
795 (xIonDense0[0][nelem][ion]-xIonDense_prev[nelem][ion]) < 0)
797 lgOscill[nelem] =
true;
818 bool lgPopsConverged =
true;
819 double old_val, new_val;
820 (*diatom)->H2_LevelPops( lgPopsConverged, old_val, new_val );
821 if( !lgPopsConverged )
837 static double OldDeut[2] = {0., 0.};
838 for(
long ion=0; ion<2; ++ion )
858 for(
long nelem=ipISO; nelem<
LIMELM; ++nelem )
889 if ( loop_ion == 0 && lgIonizWidened )
905 " ConvBase4 ionization driver loop_ion %li converged? %c reason not converged %s\n" ,
949 OldSumOTS[0] = OldSumOTS[1];
950 OldSumOTS[1] = SumOTS;
977 if( (OldSumOTS[0]-OldSumOTS[1]) * ( OldSumOTS[1] - SumOTS ) < 0. )
999 enum {DEBUG_LOC=
false};
1000 if( DEBUG_LOC && (
nzone>110) )
1030 bool lgCheckEsc =
false;
1081 hminus_old = hminus_den;
1114 " ConvBase return. fnzone %.2f nPres2Ioniz %li Te:%.3e HI:%.3e HII:%.3e H2:%.3e Ne:%.3e htot:%.3e CSUP:%.2e Conv?%c reason:%s\n",
1138 fprintf(
ioQQQ,
" PROBLEM ConvBase sets lgAbort since nPres2Ioniz exceeds limPres2Ioniz. ");
1141 fprintf(
ioQQQ,
" Reset this limit with the SET PRES IONIZ command.\n");
1179 (EdenFromMolecOld-EdenFromGrainsOld),
1217 for(
long nelem=ipISO; nelem<
LIMELM;++nelem )
1225 if( fabs(
iso_sp[ipISO][nelem].st[0].Pop()-save_iso_grnd[ipISO][nelem])/
SDIV(
iso_sp[ipISO][nelem].st[0].Pop())-1. >
1229 sprintf( chConvIoniz,
"iso %2li %2li",ipISO, nelem );
1231 save_iso_grnd[ipISO][nelem],
1232 iso_sp[ipISO][nelem].st[0].Pop());
1256 fprintf(
ioQQQ ,
"ABORT flag set since STOP nTotalIoniz was set and reached.\n");
1264 static int iter_punch=-1;
1271 "%li\t%.4e\t%.4e\t%.4e\n",
1286 for(
long nelem=0; nelem<
LIMELM; ++nelem )
1288 for(
long ion=0; ion<nelem+1; ++ion )
1302 ionbal.
UTA_heat_rate[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1] += rateone*UTALines[i].Coll().heat();
1306 enum {DEBUG_LOC=
false};
1310 fprintf(
ioQQQ,
"DEBUG UTA %3i %3i %.3f %.2e %.2e %.2e\n",
1311 (*UTALines[i].Hi()).nelem() , (*UTALines[i].Hi()).IonStg() , UTALines[i].WLAng() ,
1312 rateone, UTALines[i].Coll().heat(),
1313 UTALines[i].Coll().heat()*
dense.
xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1] );
1326 fixit(
"lgNetEdenSrcSmall needs to be enabled");
1331 fixit(
"grain rates not well tested below");
1337 double ionsrc = 0., ionsnk = 0.;
1338 for(
long nelem=0; nelem <
LIMELM; ++nelem )
1344 for(
long ion_from = 0; ion_from <= nelem + 1; ++ion_from )
1346 for(
long ion_to = 0; ion_to <= nelem + 1; ++ion_to )
1348 if( ion_to-ion_from > 0 )
1353 else if( ion_to-ion_from < 0 )
1362 const double totsrc = ionsrc +
mole.
species[ipMElec].src;
1364 const double diff = (totsrc - totsnk);
1365 const double ave = ( fabs(totsrc) + fabs(totsnk) )/2.;
1366 fixit(
"Need to tighten up e- population convergence criterion");
1367 const double error_allowed = 0.05 * ave;
1368 if( fabs(diff) > error_allowed )
1370 enum {DEBUG_LOC=
false};
1373 fprintf(
ioQQQ,
"PROBLEM large NetEdenSrc nzone %li\t%e\t%e\t%e\t%e\n",
1375 totsrc/
SDIV(totsnk),
void ion_trim(long int nelem)
t_mole_global mole_global
void RT_OTS_Update(double *SumOTS)
void DumpLine(const TransitionProxy &t)
void CoolSave(FILE *io, const char chJob[])
TransitionList UTALines("UTALines",&AnonStates)
bool lgFirstSweepThisZone
ConvergenceCounter register_(const string name)
void ion_widen(long nelem)
long int IonHigh[LIMELM+1]
void RT_OTS_PrtRate(double weak, int chFlag)
bool lgTimeDependentStatic
STATIC bool lgNetEdenSrcSmall(void)
double ChargTranSumHeat(void)
void CoolEvaluate(double *tot)
char chHashString[INPUT_LINE_LENGTH]
bool lgTemperatureConstant
molezone * findspecieslocal(const char buf[])
void RT_line_all_escape(realnum *error)
void lgStatesConserved(long nelem, long ionStage, qList states, long numStates, realnum err_tol, long loop_ion)
void mole_save(FILE *punit, const char speciesname[], const char args[], bool lgHeader, bool lgData, bool lgCoef, double depth)
t_elementnames elementnames
t_iso_sp iso_sp[NISO][LIMELM]
void PresTotCurrent(void)
STATIC void ion_trim_from_set(long nelem)
double xIonDense[LIMELM][LIMELM+1]
bool fp_equal(sys_float x, sys_float y, int n=3)
void iso_multiplet_opacities(void)
void iso_collapsed_update(void)
const char * chConvIoniz() const
FILE * ipTraceConvergeBase
long int IonLow[LIMELM+1]
void ion_trim2(long int nelem)
molecule * findspecies(const char buf[])
double ** UTA_ionize_rate
void OpacityAddTotal(void)
valarray< class molezone > species
realnum IonizErrorAllowed
const int INPUT_LINE_LENGTH
vector< diatomics * > diatoms
void ion_recom_calculate(void)
void iso_solve(long ipISO, long nelem, double &maxerr)
realnum HeatCoolRelErrorAllowed
void iso_update_rates(void)
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
realnum gas_phase[LIMELM]
void mole_update_sources(void)
void SetDeuteriumIonization(const double &xNeutral, const double &xIonized)
void mole_dominant_rates(const vector< const molecule * > &debug_list, FILE *ioOut, bool lgPrintReagents, size_t NPRINT, double fprint)
#define DEBUG_ENTRY(funcname)
void iso_renorm(long nelem, long ipISO, double &renorm)
double RateIonizTot(long nelem, long ion) const
int fprintf(const Output &stream, const char *format,...)
void ion_wrapper(long nelem)
sys_float SDIV(sys_float x)
void setConvIonizFail(const char *reason, double oldval, double newval)
bool lgElemsConserved(void)
t_secondaries secondaries
bool lgTraceConvergeBaseHash
double CharExcIonOf[NCX][LIMELM][LIMELM+1]
void ion_trim_validate(long nelem, bool lgIonizTrimCalled)
realnum GrainChTrRate[LIMELM][LIMELM+1][LIMELM+1]
double CharExcRecTo[NCX][LIMELM][LIMELM+1]
vector< diatomics * >::iterator diatom_iter
const bool lgConvBaseHeatTest