174 static double SecondOld;
175 static long int nzoneOTS=-1;
176 const int LOOP_ION_LIMIT = 10;
178 static double SumOTS=0. , OldSumOTS[2]={0.,0.};
181 IonizConverg lgIonizConverg;
208 if(
conv.nTotalIoniz )
220 EdenTrue_old =
dense.EdenTrue;
221 EdenFromMolecOld =
mole.elec;
222 EdenFromGrainsOld =
gv.TotalEden;
229 for(
long nelem=ipISO; nelem<
LIMELM;++nelem )
231 if(
dense.lgElmtOn[nelem] )
234 save_iso_grnd[ipISO][nelem] =
iso_sp[ipISO][nelem].st[0].Pop();
249 conv.lgOscilOTS =
false;
255 " ConvBase called. %.2f Te:%.3e HI:%.3e HII:%.3e H2:%.3e Ne:%.3e htot:%.3e CSUP:%.2e Conv?%c\n",
267 conv.resetConvIoniz();
275 if( !
conv.nTotalIoniz )
277 conv.setConvIonizFail(
"first call", 0.0, 0.0 );
282 conv.lgIonStageTrimed =
false;
288 opac.lgRedoStatic =
false;
291 !
opac.lgOpacStatic ||
295 conv.nPres2Ioniz == 0 )
298 opac.lgRedoStatic =
true;
310 if(
conv.nTotalIoniz )
312 bool lgIonizTrimCalled =
false;
313 static long int nZoneCalled = 0;
328 if(
conv.nTotalIoniz>2 &&
333 lgIonizTrimCalled =
true;
336 if(
dense.lgElmtOn[nelem] )
351 if(
dense.lgElmtOn[nelem] )
353 for(
long ion=0; ion<
dense.IonLow[nelem]; ++ion )
358 for(
long ion=
dense.IonLow[nelem]; ion<=
dense.IonHigh[nelem]; ++ion )
373 ASSERT(
conv.lgSearch || !lgIonizTrimCalled ||
375 (ion==0 &&
dense.IonHigh[nelem]==0 ) ||
379 for(
long ion=
dense.IonHigh[nelem]+1; ion<nelem+1; ++ion )
390 if(
conv.lgIonStageTrimed )
394 conv.setConvIonizFail(
"IonTrimmed", 0.0, 0.0 );
411 (*diatom)->CalcPhotoionizationRate();
417 if(
dense.lgSetIoniz[nelem] )
419 dense.IonLow[nelem] = 0;
420 dense.IonHigh[nelem] = nelem + 1;
421 while(
dense.SetIoniz[nelem][
dense.IonLow[nelem]] == 0. )
422 ++
dense.IonLow[nelem];
423 while(
dense.SetIoniz[nelem][
dense.IonHigh[nelem]] == 0. )
424 --
dense.IonHigh[nelem];
426 ionLowCache[nelem] =
dense.IonLow[nelem];
427 ionHighCache[nelem] =
dense.IonHigh[nelem];
443 conv.resetConvIoniz();
469 bool lgShortCircuit =
false;
470 for( ion_loop=0; ion_loop<nconv && !lgShortCircuit; ++ion_loop)
475 for (
long ion = 0; ion <= nelem+1; ++ion)
477 xIonDense0[ion_loop][nelem][ion] =
dense.xIonDense[nelem][ion];
484 if (!
dense.lgElmtOn[nelem])
488 ASSERT(
dense.IonHigh[nelem] <= ionHighCache[nelem]);
505 if (
conv.lgUpdateCouplings &&
conv.nTotalIoniz )
516 if (fabs(netion) > ion_cmp &&
520 conv.setConvIonizFail(
"OH CX inconsistency" , ion_cmp, netion);
525 bool lgCanShortCircuit = (ion_loop+1 < nconv);
528 if (!
dense.lgElmtOn[nelem])
530 for (
long ion =
dense.IonLow[nelem];
531 ion <=
dense.IonHigh[nelem]; ++ion)
533 double x0 = xIonDense0[ion_loop][nelem][ion];
534 double x1 =
dense.xIonDense[nelem][ion];
537 lgCanShortCircuit =
false;
542 lgShortCircuit = lgCanShortCircuit;
550 if (!
dense.lgElmtOn[nelem])
552 double tot0 = 0., tot1 = 0.;
554 for (
long ion = 0; ion <= nelem+1; ++ion)
556 double x0 = xIonDense0[nconv-2][nelem][ion];
557 double x1 = xIonDense0[nconv-1][nelem][ion];
558 double x2 =
dense.xIonDense[nelem][ion];
565 double extstep = 0.,predict=
x2,
566 step0 =
x1-
x0, step1 =
x2-
x1, abs1 = fabs(step1);
568 if ( abs1 > 1000.0*((
double)DBL_EPSILON)*
x2 )
570 double denom = fabs(step1-step0);
571 double sgn = (step1*step0 > 0)? 1.0 : -1.0;
574 const double MAXACC=100.0;
575 double extfac = 1.0/(denom/abs1 + 1.0/MAXACC);
576 extstep = sgn*extfac*step1;
578 predict =
x2+extstep;
580 xIonNew[ion] = predict;
584 if ( (nelem ==
ipNICKEL && ion <=2 ) )
585 fprintf(
ioQQQ,
"Extrap %3ld %3ld %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g\n",
587 x0,
x0-xIonDense0[nconv-3][nelem][ion],
x1-
x0,
x2-
x1,extstep,predict);
588 tot1 += xIonNew[ion];
592 double scal = tot0/tot1;
593 for (
long ion = 0; ion <= nelem+1; ++ion)
595 dense.xIonDense[nelem][ion] = scal*xIonNew[ion];
605 bool lgPostExtrapSolve =
true;
606 if (lgPostExtrapSolve)
610 if (!
dense.lgElmtOn[nelem])
640 if(
dense.lgSetIoniz[nelem] )
642 dense.IonLow[nelem] = 0;
643 dense.IonHigh[nelem] = nelem + 1;
644 while(
dense.SetIoniz[nelem][
dense.IonLow[nelem]] == 0. )
645 ++
dense.IonLow[nelem];
646 while(
dense.SetIoniz[nelem][
dense.IonHigh[nelem]] == 0. )
647 --
dense.IonHigh[nelem];
657 bool lgPopsConverged =
true;
658 double old_val, new_val;
659 (*diatom)->H2_LevelPops( lgPopsConverged, old_val, new_val );
660 if( !lgPopsConverged )
662 conv.setConvIonizFail(
"H2 pops", old_val, new_val);
672 lgIonizConverg(loop_ion,
conv.IonizErrorAllowed ,
false );
676 static double OldDeut[2] = {0., 0.};
677 for(
long ion=0; ion<2; ++ion )
679 if( fabs(
deut.xIonDense[ion] - OldDeut[ion] ) > 0.2*
conv.IonizErrorAllowed*fabs(OldDeut[ion]) )
681 conv.setConvIonizFail(
"D ion" , OldDeut[ion],
deut.xIonDense[ion]);
683 OldDeut[ion] =
deut.xIonDense[ion];
687 if( fabs(EdenTrue_old -
dense.EdenTrue) >
conv.EdenErrorAllowed/2.*fabs(
dense.EdenTrue) )
689 conv.setConvIonizFail(
"EdTr cng" , EdenTrue_old,
dense.EdenTrue);
694 for(
long nelem=ipISO; nelem<
LIMELM; ++nelem )
700 if(
trace.nTrConvg>=4 )
704 " ConvBase4 ionization driver loop_ion %li converged? %c reason not converged %s\n" ,
707 conv.chConvIoniz() );
713 while( !
conv.lgConvIoniz() && loop_ion < LOOP_ION_LIMIT &&
rfield.lgIonizReevaluate);
715 if(
conv.lgConvIoniz() )
751 conv.setConvIonizFail(
"SecIonRate", SecondOld,
758 if(
conv.nTotalIoniz )
761 if( fabs( hminus_old-hminus_den ) > fabs( hminus_den *
conv.EdenErrorAllowed ) )
763 conv.setConvIonizFail(
"Big H- chn", hminus_old, hminus_den );
765 hminus_old = hminus_den;
769 if( HeatOld > 0. && !
thermal.lgTemperatureConstant )
772 if( fabs(1.-
thermal.htot/HeatOld) >
conv.HeatCoolRelErrorAllowed*.5 )
774 conv.setConvIonizFail(
"Big d Heat", HeatOld,
thermal.htot);
789 if( !
conv.nPres2Ioniz ||
rfield.lgOpacityReevaluate ||
conv.lgIonStageTrimed )
800 if( !
conv.nPres2Ioniz ||
rfield.lgOpacityReevaluate ||
rfield.lgIonizReevaluate ||
807 OldSumOTS[0] = OldSumOTS[1];
808 OldSumOTS[1] = SumOTS;
828 if( fabs(1.-OldSumOTS[1]/SumOTS) >
conv.EdenErrorAllowed )
835 if( (OldSumOTS[0]-OldSumOTS[1]) * ( OldSumOTS[1] - SumOTS ) < 0. )
838 conv.lgOscilOTS =
true;
842 conv.setConvIonizFail(
"OTSRatChng" , OldSumOTS[1], SumOTS);
857 enum {DEBUG_LOC=
false};
858 if( DEBUG_LOC && (
nzone>110) )
878 enum {DEBUG_LOC=
false};
879 if( DEBUG_LOC && (
nzone>200) )
881 fprintf(
ioQQQ,
"debug otsss\t%li\t%.3e\t%.3e\t%.3e\n",
883 iso_sp[0][1].trans(15,3).Emis().ots(),
885 opac.opacity_abs[2069]);
914 " 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",
938 fprintf(
ioQQQ,
" PROBLEM ConvBase sets lgAbort since nPres2Ioniz exceeds limPres2Ioniz. ");
939 fprintf(
ioQQQ,
"Their values are %li and %li.\n",
conv.nPres2Ioniz ,
conv.limPres2Ioniz);
940 fprintf(
ioQQQ,
" Reset this limit with the SET PRES IONIZ command.\n");
954 if( fabs(EdenTrue_old -
dense.EdenTrue) > fabs(
dense.EdenTrue *
conv.EdenErrorAllowed/2.) )
956 conv.setConvIonizFail(
"eden chng", EdenTrue_old,
dense.EdenTrue);
960 if( fabs(EdenFromMolecOld -
mole.elec) > fabs(
dense.EdenTrue *
conv.EdenErrorAllowed/2.) )
962 conv.setConvIonizFail(
"edn chnCO", EdenFromMolecOld,
dense.EdenTrue);
965 if(
gv.lgGrainElectrons )
968 if( fabs(EdenFromGrainsOld -
gv.TotalEden) > fabs(
dense.EdenTrue *
conv.EdenErrorAllowed/2.) )
970 conv.setConvIonizFail(
"edn grn e", EdenFromGrainsOld,
gv.TotalEden);
974 if( fabs( (EdenFromMolecOld-EdenFromGrainsOld) - (
mole.elec-
gv.TotalEden) ) >
975 fabs(
dense.EdenTrue *
conv.EdenErrorAllowed/4.) )
977 conv.setConvIonizFail(
"edn mole-grn",
978 (EdenFromMolecOld-EdenFromGrainsOld),
979 (
mole.elec-
gv.TotalEden));
992 if( !
thermal.lgTemperatureConstant )
996 conv.setConvIonizFail(
"heat chg", HeatingOld,
thermal.htot);
1001 conv.setConvIonizFail(
"cool chg", CoolingOld,
thermal.ctot);
1009 if( fabs(
mole.species[i].den-mole_save[i])/
dense.xNucleiTotal-1. >
1010 conv.EdenErrorAllowed/2.)
1013 sprintf( chConvIoniz,
"ch %-4.4s",
mole_global.list[i]->label.c_str() );
1014 conv.setConvIonizFail( chConvIoniz,
1015 mole_save[i]/
dense.xNucleiTotal,
1016 mole.species[i].den/
dense.xNucleiTotal);
1025 for(
long nelem=ipISO; nelem<
LIMELM;++nelem )
1027 if(
dense.lgElmtOn[nelem] )
1031 if(
dense.xIonDense[nelem][nelem-ipISO]/
dense.gas_phase[nelem] > 1e-5 )
1033 if( fabs(
iso_sp[ipISO][nelem].st[0].Pop()-save_iso_grnd[ipISO][nelem])/
SDIV(
iso_sp[ipISO][nelem].st[0].Pop())-1. >
1034 conv.EdenErrorAllowed)
1037 sprintf( chConvIoniz,
"iso %2li %2li",ipISO, nelem );
1038 conv.setConvIonizFail(chConvIoniz,
1039 save_iso_grnd[ipISO][nelem],
1040 iso_sp[ipISO][nelem].st[0].Pop());
1056 if( !
conv.lgSearch )
1057 conv.lgFirstSweepThisZone =
false;
1064 fprintf(
ioQQQ ,
"ABORT flag set since STOP nTotalIoniz was set and reached.\n");
1068 if(
save.lgTraceConvergeBase )
1072 static int iter_punch=-1;
1074 fprintf(
save.ipTraceConvergeBase,
"%s\n",
save.chHashString );
1078 fprintf(
save.ipTraceConvergeBase,
1079 "%li\t%.4e\t%.4e\t%.4e\n",
1147 double ionsrc = 0., ionsnk = 0.;
1148 for(
long nelem=0; nelem <
LIMELM; ++nelem )
1150 if( !
dense.lgElmtOn[nelem] )
1152 ionsrc +=
ionbal.elecsrc[nelem];
1153 ionsnk +=
ionbal.elecsnk[nelem];
1154 for(
long ion_from = 0; ion_from <= nelem + 1; ++ion_from )
1156 for(
long ion_to = 0; ion_to <= nelem + 1; ++ion_to )
1158 if( ion_to-ion_from > 0 )
1160 ionsrc +=
gv.GrainChTrRate[nelem][ion_from][ion_to] *
1161 dense.xIonDense[nelem][ion_from] * (ion_to-ion_from);
1163 else if( ion_to-ion_from < 0 )
1165 ionsnk +=
gv.GrainChTrRate[nelem][ion_from][ion_to] *
1166 dense.xIonDense[nelem][ion_from] * (ion_from-ion_to);
1172 const double totsrc = ionsrc +
mole.species[ipMElec].src;
1173 const double totsnk = ionsnk +
mole.species[ipMElec].snk*
dense.EdenTrue;
1174 const double diff = (totsrc - totsnk);
1175 const double ave = ( fabs(totsrc) + fabs(totsnk) )/2.;
1177 const double error_allowed = 0.05 * ave;
1178 if( fabs(diff) > error_allowed )
1180 enum {DEBUG_LOC=
false};
1183 fprintf(
ioQQQ,
"PROBLEM large NetEdenSrc nzone %li\t%e\t%e\t%e\t%e\n",
1185 totsrc/
SDIV(totsnk),
1187 ionsrc +
mole.species[ipMElec].src,
1188 ionsnk +
mole.species[ipMElec].snk*
dense.EdenTrue);
1190 conv.setConvIonizFail(
"NetEdenSrc", diff, error_allowed);