00001
00002
00003
00004
00005
00006 #include "cddefines.h"
00007 #include "dynamics.h"
00008 #include "trace.h"
00009 #include "elementnames.h"
00010 #include "save.h"
00011 #include "phycon.h"
00012 #include "secondaries.h"
00013 #include "stopcalc.h"
00014 #include "grainvar.h"
00015 #include "highen.h"
00016 #include "dense.h"
00017 #include "hmi.h"
00018 #include "rfield.h"
00019 #include "pressure.h"
00020 #include "taulines.h"
00021 #include "rt.h"
00022 #include "grains.h"
00023 #include "atmdat.h"
00024 #include "ionbal.h"
00025 #include "opacity.h"
00026 #include "cooling.h"
00027 #include "thermal.h"
00028 #include "mole.h"
00029 #include "iso.h"
00030 #include "conv.h"
00031 #include "h2.h"
00032 #include "deuterium.h"
00033
00034 STATIC bool lgNetEdenSrcSmall( void );
00035
00036 void UpdateUTAs( void );
00037
00038
00039
00040 namespace
00041 {
00042 class IonizConverg
00043 {
00044
00045 double OldFracs[LIMELM+1][LIMELM+1];
00046 public:
00047 IonizConverg()
00048 {
00049 for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
00050 {
00051 for (long ion = 0; ion <= nelem+1; ++ion)
00052 {
00053 OldFracs[nelem][ion] = dense.xIonDense[nelem][ion];
00054 }
00055 }
00056 }
00057 void operator()(
00058
00059 long loop_ion ,
00060
00061 double delta ,
00062
00063 bool lgPrint )
00064 {
00065 bool lgConverg_v = false;
00066 double Abund,
00067 bigchange ,
00068 change ,
00069 one;
00070
00071 DEBUG_ENTRY( "IonizConverg::operator()" );
00072
00073 for (long nelem =ipHYDROGEN; nelem<LIMELM; ++nelem )
00074 {
00075 if( !dense.lgElmtOn[nelem] )
00076 continue;
00077
00078 double abundold=0. , abundnew=0.;
00079 long ionchg=-1;
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091 if( conv.nPres2Ioniz )
00092 {
00093
00094
00095 bigchange = 0.;
00096 change = 0.;
00097 Abund = dense.gas_phase[nelem];
00098
00099
00100
00101 for( long ion=0; ion <= (nelem+1); ++ion )
00102 {
00103
00104 if( OldFracs[nelem][ion]/Abund > 1e-4 &&
00105 dense.xIonDense[nelem][ion]/Abund > 1e-4 )
00106
00107 {
00108
00109 one = fabs(dense.xIonDense[nelem][ion]-OldFracs[nelem][ion])/
00110 OldFracs[nelem][ion];
00111 change = MAX2(change, one );
00112
00113 if( change>bigchange )
00114 {
00115 bigchange = change;
00116 abundold = OldFracs[nelem][ion]/Abund;
00117 abundnew = dense.xIonDense[nelem][ion]/Abund;
00118 ionchg = ion;
00119 }
00120 }
00121
00122 OldFracs[nelem][ion] = dense.xIonDense[nelem][ion];
00123 }
00124
00125 if( change >= delta )
00126 {
00127 char chConvIoniz[INPUT_LINE_LENGTH];
00128 sprintf( chConvIoniz , "%2s ion" , elementnames.chElementSym[nelem] );
00129 conv.setConvIonizFail(chConvIoniz,abundold,abundnew);
00130 ASSERT( abundold>0. && abundnew>0. );
00131 }
00132 else
00133 lgConverg_v = true;
00134 }
00135 else
00136 {
00137 for( long ion=0; ion <= (nelem+1); ++ion )
00138 {
00139 OldFracs[nelem][ion] = dense.xIonDense[nelem][ion];
00140 }
00141
00142 }
00143
00144
00145 if( lgPrint )
00146 {
00147 fprintf(ioQQQ," nz %ld loop %ld element %li converged? %c worst %ld change %g\n",
00148 nzone, loop_ion, nelem, TorF(lgConverg_v),ionchg,bigchange);
00149 for( long ion=0; ion<(nelem+1); ++ion )
00150 {
00151 fprintf(ioQQQ,"\t%.5e", dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]);
00152 }
00153 fprintf(ioQQQ,"\n");
00154 }
00155 }
00156 }
00157 };
00158 }
00159
00160
00161
00162
00163 int ConvBase(
00164
00165
00166 long loopi )
00167 {
00168 double HeatOld,
00169 EdenTrue_old,
00170 EdenFromMolecOld,
00171 EdenFromGrainsOld,
00172 HeatingOld ,
00173 CoolingOld;
00174 static double SecondOld;
00175 static long int nzoneOTS=-1;
00176 # define LOOP_ION_LIMIT 10
00177 long int loop_ion;
00178 static double SumOTS=0. , OldSumOTS[2]={0.,0.};
00179 double save_iso_grnd[NISO][LIMELM];
00180 valarray<realnum> mole_save(mole_global.num_calc);
00181 IonizConverg lgIonizConverg;
00182
00183 DEBUG_ENTRY( "ConvBase()" );
00184
00185
00186
00187
00188 ASSERT( fp_equal( phycon.te, thermal.te_update ) );
00189
00190
00191
00192 fnzone = (double)nzone + (double)conv.nPres2Ioniz/100.;
00193
00194
00195
00196
00197 PresTotCurrent();
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208 if( conv.nTotalIoniz )
00209 {
00210 if( eden_sum() )
00211 {
00212
00213 ++conv.nTotalIoniz;
00214 return 1;
00215 }
00216 }
00217
00218
00219
00220 EdenTrue_old = dense.EdenTrue;
00221 EdenFromMolecOld = mole.elec;
00222 EdenFromGrainsOld = gv.TotalEden;
00223 HeatingOld = thermal.htot;
00224 CoolingOld = thermal.ctot;
00225
00226
00227 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00228 {
00229 for( long nelem=ipISO; nelem<LIMELM;++nelem )
00230 {
00231 if( dense.lgElmtOn[nelem] )
00232 {
00233
00234 save_iso_grnd[ipISO][nelem] = iso_sp[ipISO][nelem].st[0].Pop();
00235 }
00236 }
00237 }
00238
00239 for( long i=0; i < mole_global.num_calc; i++ )
00240 {
00241 mole_save[i] = (realnum) mole.species[i].den;
00242 }
00243
00244 if( loopi==0 )
00245 {
00246
00247 OldSumOTS[0] = 0.;
00248 OldSumOTS[1] = 0.;
00249 conv.lgOscilOTS = false;
00250 }
00251
00252 if( trace.lgTrace )
00253 {
00254 fprintf( ioQQQ,
00255 " ConvBase called. %.2f Te:%.3e HI:%.3e HII:%.3e H2:%.3e Ne:%.3e htot:%.3e CSUP:%.2e Conv?%c\n",
00256 fnzone,
00257 phycon.te,
00258 dense.xIonDense[ipHYDROGEN][0],
00259 dense.xIonDense[ipHYDROGEN][1],
00260 hmi.H2_total,
00261 dense.eden,
00262 thermal.htot,
00263 secondaries.csupra[ipHYDROGEN][0] ,
00264 TorF(conv.lgConvIoniz()) );
00265 }
00266
00267 conv.resetConvIoniz();
00268
00269
00270 HeatZero();
00271
00272
00273
00274
00275 if( !conv.nTotalIoniz )
00276 {
00277 conv.setConvIonizFail( "first call", 0.0, 0.0 );
00278 }
00279
00280
00281
00282 conv.lgIonStageTrimed = false;
00283
00284
00285
00286
00287
00288 opac.lgRedoStatic = false;
00289 if(
00290
00291 !opac.lgOpacStatic ||
00292
00293 conv.lgSearch ||
00294
00295 conv.nPres2Ioniz == 0 )
00296 {
00297
00298 opac.lgRedoStatic = true;
00299 }
00300
00301
00302 ion_recom_calculate();
00303
00304
00305
00306
00307
00308
00309
00310 if( conv.nTotalIoniz )
00311 {
00312 bool lgIonizTrimCalled = false;
00313 static long int nZoneCalled = 0;
00314
00315 fixit();
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328 if( conv.nTotalIoniz>2 &&
00329
00330
00331 ( conv.lgSearch || nZoneCalled!=nzone) )
00332 {
00333 lgIonizTrimCalled = true;
00334 for( long nelem=ipHELIUM; nelem<LIMELM; ++nelem )
00335 {
00336 if( dense.lgElmtOn[nelem] )
00337 {
00338
00339
00340 ion_trim(nelem);
00341 }
00342 }
00343 nZoneCalled = nzone;
00344 }
00345
00346
00347 # if !defined(NDEBUG)
00348
00349 for( long nelem=ipHELIUM; nelem<LIMELM; ++nelem)
00350 {
00351 if( dense.lgElmtOn[nelem] )
00352 {
00353 for( long ion=0; ion<dense.IonLow[nelem]; ++ion )
00354 {
00355 ASSERT( dense.xIonDense[nelem][ion] == 0. );
00356 }
00357
00358 for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
00359 {
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373 ASSERT( conv.lgSearch || !lgIonizTrimCalled ||
00374
00375 (ion==0 && dense.IonHigh[nelem]==0 ) ||
00376 dense.xIonDense[nelem][ion] >= SMALLFLOAT ||
00377 dense.xIonDense[nelem][ion]/dense.gas_phase[nelem] >= SMALLFLOAT );
00378 }
00379 for( long ion=dense.IonHigh[nelem]+1; ion<nelem+1; ++ion )
00380 {
00381 ASSERT( ion >= 0 );
00382 ASSERT( dense.xIonDense[nelem][ion] == 0. );
00383 }
00384 }
00385 }
00386 # endif
00387 }
00388
00389
00390 if( conv.lgIonStageTrimed )
00391 {
00392
00393
00394 conv.setConvIonizFail( "IonTrimmed", 0.0, 0.0 );
00395 }
00396
00397
00398 if( dynamics.lgAdvection )
00399 DynaIonize();
00400
00401
00402 highen();
00403
00404
00405 iso_collapsed_update();
00406
00407 iso_update_rates( );
00408
00409 long ionLowCache[LIMELM], ionHighCache[LIMELM];
00410 for( long nelem=ipHYDROGEN; nelem<LIMELM; nelem++ )
00411 {
00412 if( dense.lgSetIoniz[nelem] )
00413 {
00414 dense.IonLow[nelem] = 0;
00415 dense.IonHigh[nelem] = nelem + 1;
00416 while( dense.SetIoniz[nelem][dense.IonLow[nelem]] == 0. )
00417 ++dense.IonLow[nelem];
00418 while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] == 0. )
00419 --dense.IonHigh[nelem];
00420 }
00421 ionLowCache[nelem] = dense.IonLow[nelem];
00422 ionHighCache[nelem] = dense.IonHigh[nelem];
00423 }
00424
00425
00426
00427
00428 conv.incrementCounter(CONV_BASE_CALLS);
00429 loop_ion = 0;
00430 do
00431 {
00432 conv.incrementCounter(CONV_BASE_LOOPS);
00433 ASSERT(lgElemsConserved());
00434
00435
00436
00437 if( loop_ion )
00438 conv.resetConvIoniz();
00439
00440
00441
00442
00443
00444 ChargTranEval();
00445
00446
00447
00448 GrainDrive();
00449
00450
00451 highen();
00452
00453
00454 atmdat_3body();
00455
00456
00457 UpdateUTAs();
00458
00459 ASSERT(lgElemsConserved());
00460
00461
00462 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
00463 {
00464 bool lgPopsConverged = true;
00465 double old_val, new_val;
00466 (*diatom)->H2_LevelPops( lgPopsConverged, old_val, new_val );
00467 if( !lgPopsConverged )
00468 {
00469 conv.setConvIonizFail( "H2 pops", old_val, new_val);
00470 }
00471 }
00472
00473 ASSERT(lgElemsConserved());
00474
00475 long ion_loop = 0;
00476 const int nconv = 3;
00477 double xIonDense0[nconv][LIMELM][LIMELM+1];
00478 bool lgShortCircuit = false;
00479 for( ion_loop=0; ion_loop<nconv && !lgShortCircuit; ++ion_loop)
00480 {
00481 conv.incrementCounter(CONV_BASE_ACCELS);
00482 for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
00483 {
00484 for (long ion = 0; ion <= nelem+1; ++ion)
00485 {
00486 xIonDense0[ion_loop][nelem][ion] = dense.xIonDense[nelem][ion];
00487 }
00488 }
00489 for( long nelem=ipHYDROGEN; nelem<LIMELM; nelem++ )
00490 {
00491 if (!dense.lgElmtOn[nelem])
00492 continue;
00493
00494 ASSERT(dense.IonLow[nelem] >= ionLowCache[nelem]);
00495 ASSERT(dense.IonHigh[nelem] <= ionHighCache[nelem]);
00496 ion_wrapper( nelem );
00497
00498 if ( conv.lgUpdateCouplings )
00499 mole_update_sources();
00500
00501 }
00502
00503 bool lgCanShortCircuit = (ion_loop+1 < nconv);
00504 for( long nelem=ipHYDROGEN; nelem<LIMELM && lgCanShortCircuit; ++nelem )
00505 {
00506 if (!dense.lgElmtOn[nelem])
00507 continue;
00508 for (long ion = dense.IonLow[nelem];
00509 ion <= dense.IonHigh[nelem]; ++ion)
00510 {
00511 double x0 = xIonDense0[ion_loop][nelem][ion];
00512 double x1 = dense.xIonDense[nelem][ion];
00513 if (fabs(x0-x1) > 1e-6*(x0+x1))
00514 {
00515 lgCanShortCircuit = false;
00516 break;
00517 }
00518 }
00519 }
00520 lgShortCircuit = lgCanShortCircuit;
00521 }
00522
00523 if (!lgShortCircuit)
00524 {
00525
00526 for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
00527 {
00528 if (!dense.lgElmtOn[nelem])
00529 continue;
00530 double tot0 = 0., tot1 = 0.;
00531 double xIonNew[LIMELM+1];
00532 for (long ion = 0; ion <= nelem+1; ++ion)
00533 {
00534 double x0 = xIonDense0[nconv-2][nelem][ion];
00535 double x1 = xIonDense0[nconv-1][nelem][ion];
00536 double x2 = dense.xIonDense[nelem][ion];
00537 xIonNew[ion] = x2;
00538 tot0 += x2;
00539
00540
00541
00542
00543 double extstep = 0.,predict=x2,
00544 step0 = x1-x0, step1 = x2-x1, abs1 = fabs(step1);
00545
00546 if ( abs1 > 1000.0*((double)DBL_EPSILON)*x2 )
00547 {
00548 double denom = fabs(step1-step0);
00549 double sgn = (step1*step0 > 0)? 1.0 : -1.0;
00550
00551
00552 const double MAXACC=100.0;
00553 double extfac = 1.0/(denom/abs1 + 1.0/MAXACC);
00554 extstep = sgn*extfac*step1;
00555
00556 predict = x2+extstep;
00557 if (predict > 0.0)
00558 xIonNew[ion] = predict;
00559 }
00560 if ( 0 )
00561
00562 if ( (nelem == ipNICKEL && ion <=2 ) )
00563 fprintf(ioQQQ,"Extrap %3ld %3ld %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g\n",
00564 nelem,ion,
00565 x0,x0-xIonDense0[nconv-3][nelem][ion],x1-x0,x2-x1,extstep,predict);
00566 tot1 += xIonNew[ion];
00567 }
00568 if ( tot1 > SMALLFLOAT )
00569 {
00570 double scal = tot0/tot1;
00571 for (long ion = 0; ion <= nelem+1; ++ion)
00572 {
00573 dense.xIonDense[nelem][ion] = scal*xIonNew[ion];
00574 }
00575 for ( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00576 {
00577 double renorm;
00578 iso_renorm(nelem, ipISO, renorm);
00579 }
00580 }
00581 }
00582
00583 bool lgPostExtrapSolve = true;
00584 if (lgPostExtrapSolve)
00585 {
00586 for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
00587 {
00588 if (!dense.lgElmtOn[nelem])
00589 continue;
00590 ion_wrapper( nelem );
00591 }
00592 }
00593 }
00594 SetDeuteriumIonization( dense.xIonDense[ipHYDROGEN][0], dense.xIonDense[ipHYDROGEN][1] );
00595
00596
00597
00598
00599 if( lgAbort )
00600 {
00601 ++conv.nTotalIoniz;
00602 return 1;
00603 }
00604
00605 ASSERT(lgElemsConserved());
00606
00607
00608 mole_drive();
00609
00610 ASSERT(lgElemsConserved());
00611
00612
00613
00614
00615
00616 for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
00617 {
00618 if( dense.lgSetIoniz[nelem] )
00619 {
00620 dense.IonLow[nelem] = 0;
00621 dense.IonHigh[nelem] = nelem + 1;
00622 while( dense.SetIoniz[nelem][dense.IonLow[nelem]] == 0. )
00623 ++dense.IonLow[nelem];
00624 while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] == 0. )
00625 --dense.IonHigh[nelem];
00626 }
00627 }
00628
00629
00630 ChargTranEval();
00631
00632
00633
00634
00635
00636 lgIonizConverg(loop_ion, conv.IonizErrorAllowed , false );
00637
00638 if( deut.lgElmtOn )
00639 {
00640 static double OldDeut[2] = {0., 0.};
00641 for( long ion=0; ion<2; ++ion )
00642 {
00643 if( fabs(deut.xIonDense[ion] - OldDeut[ion] ) > 0.2*conv.IonizErrorAllowed*fabs(OldDeut[ion]) )
00644 {
00645 conv.setConvIonizFail( "D ion" , OldDeut[ion], deut.xIonDense[ion]);
00646 }
00647 OldDeut[ion] = deut.xIonDense[ion];
00648 }
00649 }
00650
00651 if( fabs(EdenTrue_old - dense.EdenTrue) > conv.EdenErrorAllowed/2.*fabs(dense.EdenTrue) )
00652 {
00653 conv.setConvIonizFail( "EdTr cng" , EdenTrue_old, dense.EdenTrue);
00654 }
00655
00656 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00657 {
00658 for( long nelem=ipISO; nelem<LIMELM; ++nelem )
00659 {
00660 lgStatesConserved( nelem, nelem-ipISO, iso_sp[ipISO][nelem].st, iso_sp[ipISO][nelem].numLevels_local, 1e-3f, loop_ion );
00661 }
00662 }
00663
00664 if( trace.nTrConvg>=4 )
00665 {
00666
00667 fprintf( ioQQQ,
00668 " ConvBase4 ionization driver loop_ion %li converged? %c reason not converged %s\n" ,
00669 loop_ion ,
00670 TorF( conv.lgConvIoniz()) ,
00671 conv.chConvIoniz() );
00672 }
00673
00674 ++loop_ion;
00675 }
00676
00677 while( !conv.lgConvIoniz() && loop_ion < LOOP_ION_LIMIT && rfield.lgIonizReevaluate);
00678
00679 if( conv.lgConvIoniz() )
00680 lgNetEdenSrcSmall();
00681
00682
00683
00684
00685 thermal.char_tran_heat = ChargTranSumHeat();
00686
00687 thermal.char_tran_cool = MAX2(0. , -thermal.char_tran_heat );
00688 thermal.char_tran_heat = MAX2(0. , thermal.char_tran_heat );
00689
00690
00691
00692 CoolEvaluate(&thermal.ctot );
00693
00694
00695 HeatOld = thermal.htot;
00696
00697
00698
00699
00700 SecondOld = secondaries.csupra[ipHYDROGEN][0];
00701
00702
00703
00704 HeatSum();
00705
00706
00707
00708
00709
00710 if( (secondaries.csupra[ipHYDROGEN][0]>SMALLFLOAT) && (secondaries.sec2total>0.001) &&
00711 fabs( 1. - SecondOld/SDIV(secondaries.csupra[ipHYDROGEN][0]) ) > 0.1 &&
00712 SecondOld > 0. && secondaries.csupra[ipHYDROGEN][0] > 0.)
00713 {
00714
00715 conv.setConvIonizFail( "SecIonRate", SecondOld,
00716 secondaries.csupra[ipHYDROGEN][0] );
00717 }
00718
00719 # if 0
00720 static realnum hminus_old=0.;
00721
00722 if( conv.nTotalIoniz )
00723 {
00724 realnum hminus_den = findspecieslocal("H-")->den;
00725 if( fabs( hminus_old-hminus_den ) > fabs( hminus_den * conv.EdenErrorAllowed ) )
00726 {
00727 conv.setConvIonizFail( "Big H- chn", hminus_old, hminus_den );
00728 }
00729 hminus_old = hminus_den;
00730 }
00731 # endif
00732
00733 if( HeatOld > 0. && !thermal.lgTemperatureConstant )
00734 {
00735
00736 if( fabs(1.-thermal.htot/HeatOld) > conv.HeatCoolRelErrorAllowed*.5 )
00737 {
00738 conv.setConvIonizFail( "Big d Heat", HeatOld, thermal.htot);
00739 }
00740 }
00741
00742
00743 if( lgAbort )
00744 {
00745
00746 return 1;
00747 }
00748
00749
00750
00751
00752
00753 if( !conv.nPres2Ioniz || rfield.lgOpacityReevaluate || conv.lgIonStageTrimed )
00754 OpacityAddTotal();
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764 if( !conv.nPres2Ioniz || rfield.lgOpacityReevaluate || rfield.lgIonizReevaluate ||
00765 conv.lgIonStageTrimed || conv.lgSearch || nzone!=nzoneOTS )
00766 {
00767 RT_OTS();
00768 nzoneOTS = nzone;
00769
00770
00771 OldSumOTS[0] = OldSumOTS[1];
00772 OldSumOTS[1] = SumOTS;
00773
00774
00775
00776
00777
00778
00779
00780 RT_OTS_Update(&SumOTS);
00781
00782 }
00783 else
00784 SumOTS = 0.;
00785
00786
00787 if( SumOTS> 0. )
00788 {
00789
00790
00791
00792 if( fabs(1.-OldSumOTS[1]/SumOTS) > conv.EdenErrorAllowed )
00793 {
00794
00795
00796 if( loopi > 1 )
00797 {
00798
00799 if( (OldSumOTS[0]-OldSumOTS[1]) * ( OldSumOTS[1] - SumOTS ) < 0. )
00800 {
00801
00802 conv.lgOscilOTS = true;
00803 }
00804 }
00805
00806 conv.setConvIonizFail( "OTSRatChng" , OldSumOTS[1], SumOTS);
00807 }
00808
00809
00810
00811 if( ( trace.lgTrace && trace.lgOTSBug ) ||
00812 ( trace.nTrConvg && trace.lgOTSBug ) )
00813 {
00814 RT_OTS_PrtRate(SumOTS*0.05 , 'b' );
00815 }
00816
00817
00818
00819 {
00820
00821 enum {DEBUG_LOC=false};
00822 if( DEBUG_LOC && (nzone>110) )
00823 {
00824 # if 0
00825 # include "lines_service.h"
00826 DumpLine( &iso_sp[ipH_LIKE][ipHYDROGEN].trans(2,0) );
00827 # endif
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837 }
00838 }
00839 }
00840 {
00841
00842 enum {DEBUG_LOC=false};
00843 if( DEBUG_LOC && (nzone>200) )
00844 {
00845 fprintf(ioQQQ,"debug otsss\t%li\t%.3e\t%.3e\t%.3e\n",
00846 nzone,
00847 iso_sp[0][1].trans(15,3).Emis().ots(),
00848 TauLines[26].Emis().ots(),
00849 opac.opacity_abs[2069]);
00850 }
00851 }
00852
00853
00854 if( trace.lgTrace && trace.lgOTSBug )
00855 {
00856
00857
00858
00859 RT_OTS_PrtRate(SumOTS*0.05 , 'b' );
00860 }
00861
00862
00863 RT_line_all( );
00864
00865
00866
00867
00868 PresTotCurrent();
00869
00870 ASSERT(lgElemsConserved());
00871
00872
00873
00874
00875 if( trace.lgTrace )
00876 {
00877 fprintf( ioQQQ,
00878 " 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",
00879 fnzone,
00880 conv.nPres2Ioniz,
00881 phycon.te,
00882 dense.xIonDense[ipHYDROGEN][0],
00883 dense.xIonDense[ipHYDROGEN][1],
00884 hmi.H2_total,
00885 dense.eden,
00886 thermal.htot,
00887 secondaries.csupra[ipHYDROGEN][0] ,
00888 TorF(conv.lgConvIoniz()) ,
00889 conv.chConvIoniz());
00890 }
00891
00892
00893
00894
00895 ++conv.nPres2Ioniz;
00896
00897
00898
00899
00900 if( conv.nPres2Ioniz > conv.limPres2Ioniz && nzone > 0)
00901 {
00902 fprintf(ioQQQ,"PROBLEM ConvBase sets lgAbort since nPres2Ioniz exceeds limPres2Ioniz. ");
00903 fprintf(ioQQQ,"Their values are %li and %li.\n",conv.nPres2Ioniz , conv.limPres2Ioniz);
00904 lgAbort = true;
00905 return 1;
00906 }
00907
00908
00909 if( eden_sum() )
00910 {
00911
00912 return 1;
00913 }
00914
00915
00917 if( fabs(EdenTrue_old - dense.EdenTrue) > fabs(dense.EdenTrue * conv.EdenErrorAllowed/2.) )
00918 {
00919 conv.setConvIonizFail( "eden chng", EdenTrue_old, dense.EdenTrue);
00920 }
00921
00922
00923 if( fabs(EdenFromMolecOld - mole.elec) > fabs(dense.EdenTrue * conv.EdenErrorAllowed/2.) )
00924 {
00925 conv.setConvIonizFail( "edn chnCO", EdenFromMolecOld, dense.EdenTrue);
00926 }
00927
00928 if( gv.lgGrainElectrons )
00929 {
00930
00931 if( fabs(EdenFromGrainsOld - gv.TotalEden) > fabs(dense.EdenTrue * conv.EdenErrorAllowed/2.) )
00932 {
00933 conv.setConvIonizFail( "edn grn e", EdenFromGrainsOld, gv.TotalEden);
00934 }
00935
00936
00937 if( fabs( (EdenFromMolecOld-EdenFromGrainsOld) - (mole.elec-gv.TotalEden) ) >
00938 fabs(dense.EdenTrue * conv.EdenErrorAllowed/4.) )
00939 {
00940 conv.setConvIonizFail( "edn mole-grn",
00941 (EdenFromMolecOld-EdenFromGrainsOld),
00942 (mole.elec-gv.TotalEden));
00943 }
00944 }
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955 if( !thermal.lgTemperatureConstant )
00956 {
00957 if( fabs(HeatingOld - thermal.htot)/thermal.htot > conv.HeatCoolRelErrorAllowed/2. )
00958 {
00959 conv.setConvIonizFail( "heat chg", HeatingOld, thermal.htot);
00960 }
00961
00962 if( fabs(CoolingOld - thermal.ctot)/thermal.ctot > conv.HeatCoolRelErrorAllowed/2. )
00963 {
00964 conv.setConvIonizFail( "cool chg", CoolingOld, thermal.ctot);
00965 }
00966 }
00967
00968
00969 for( long i=0; i < mole_global.num_calc; ++i )
00970 {
00971
00972 if( fabs(mole.species[i].den-mole_save[i])/dense.xNucleiTotal-1. >
00973 conv.EdenErrorAllowed/2.)
00974 {
00975 char chConvIoniz[INPUT_LINE_LENGTH];
00976 sprintf( chConvIoniz, "ch %-4.4s",mole_global.list[i]->label.c_str() );
00977 conv.setConvIonizFail( chConvIoniz,
00978 mole_save[i]/dense.xNucleiTotal,
00979 mole.species[i].den/dense.xNucleiTotal);
00980 }
00981 }
00982
00983
00984
00985
00986 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00987 {
00988 for( long nelem=ipISO; nelem<LIMELM;++nelem )
00989 {
00990 if( dense.lgElmtOn[nelem] )
00991 {
00992
00993
00994 if( dense.xIonDense[nelem][nelem-ipISO]/dense.gas_phase[nelem] > 1e-5 )
00995 {
00996 if( fabs(iso_sp[ipISO][nelem].st[0].Pop()-save_iso_grnd[ipISO][nelem])/SDIV(iso_sp[ipISO][nelem].st[0].Pop())-1. >
00997 conv.EdenErrorAllowed)
00998 {
00999 char chConvIoniz[INPUT_LINE_LENGTH];
01000 sprintf( chConvIoniz,"iso %2li %2li",ipISO, nelem );
01001 conv.setConvIonizFail(chConvIoniz,
01002 save_iso_grnd[ipISO][nelem],
01003 iso_sp[ipISO][nelem].st[0].Pop());
01004 }
01005 }
01006
01007 }
01008 }
01009 }
01010
01011
01012
01013
01014
01015
01016 ++conv.nTotalIoniz;
01017
01018
01019 if( !conv.lgSearch )
01020 conv.lgFirstSweepThisZone = false;
01021
01022
01023
01024 if(StopCalc.nTotalIonizStop && conv.nTotalIoniz> StopCalc.nTotalIonizStop )
01025 {
01026
01027 fprintf(ioQQQ , "ABORT flag set since STOP nTotalIoniz was set and reached.\n");
01028 return 1;
01029 }
01030
01031 if( save.lgTraceConvergeBase )
01032 {
01033 if( iteration > 1 && save.lgTraceConvergeBaseHash )
01034 {
01035 static int iter_punch=-1;
01036 if( iteration !=iter_punch )
01037 fprintf( save.ipTraceConvergeBase, "%s\n",save.chHashString );
01038 iter_punch = iteration;
01039 }
01040
01041 fprintf( save.ipTraceConvergeBase,
01042 "%li\t%.4e\t%.4e\t%.4e\n",
01043 nzone , thermal.htot , thermal.ctot , dense.eden );
01044 }
01045
01046 return 0;
01047 }
01048
01049 void UpdateUTAs( void )
01050 {
01051 DEBUG_ENTRY( "UpdateUTAs()" );
01052
01053
01054 if( conv.lgFirstSweepThisZone )
01055 {
01056
01057 for( long nelem=0; nelem< LIMELM; ++nelem )
01058 {
01059 for( long ion=0; ion<nelem+1; ++ion )
01060 {
01061 ionbal.UTA_ionize_rate[nelem][ion] = 0.;
01062 ionbal.UTA_heat_rate[nelem][ion] = 0.;
01063 }
01064 }
01065
01066
01067 if( ionbal.lgInnerShellLine_on )
01068 {
01069 for( long i=0; i < nUTA; ++i )
01070 {
01071
01072 double rateone = UTALines[i].Emis().pump() * UTALines[i].Emis().AutoIonizFrac();
01073 ionbal.UTA_ionize_rate[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1] += rateone;
01074
01075 ionbal.UTA_heat_rate[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1] += rateone*UTALines[i].Coll().heat();
01076 {
01077
01078
01079 enum {DEBUG_LOC=false};
01080
01081 if( DEBUG_LOC )
01082 {
01083 fprintf(ioQQQ,"DEBUG UTA %3i %3i %.3f %.2e %.2e %.2e\n",
01084 (*UTALines[i].Hi()).nelem() , (*UTALines[i].Hi()).IonStg() , UTALines[i].WLAng() ,
01085 rateone, UTALines[i].Coll().heat(),
01086 UTALines[i].Coll().heat()*dense.xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1] );
01087 }
01088 }
01089 }
01090 }
01091 }
01092
01093 return;
01094 }
01095
01096 STATIC bool lgNetEdenSrcSmall( void )
01097 {
01098 fixit();
01099 return true;
01100 DEBUG_ENTRY( "lgNetEdenSrcSmall()" );
01101
01102 if( conv.lgSearch )
01103 return true;
01104 fixit();
01105 if( gv.lgDustOn() )
01106 return true;
01107
01108
01109
01110 double ionsrc = 0., ionsnk = 0.;
01111 for( long nelem=0; nelem < LIMELM; ++nelem )
01112 {
01113 if( !dense.lgElmtOn[nelem] )
01114 continue;
01115 ionsrc += ionbal.elecsrc[nelem];
01116 ionsnk += ionbal.elecsnk[nelem];
01117 for( long ion_from = 0; ion_from <= nelem + 1; ++ion_from )
01118 {
01119 for( long ion_to = 0; ion_to <= nelem + 1; ++ion_to )
01120 {
01121 if( ion_to-ion_from > 0 )
01122 {
01123 ionsrc += gv.GrainChTrRate[nelem][ion_from][ion_to] *
01124 dense.xIonDense[nelem][ion_from] * (ion_to-ion_from);
01125 }
01126 else if( ion_to-ion_from < 0 )
01127 {
01128 ionsnk += gv.GrainChTrRate[nelem][ion_from][ion_to] *
01129 dense.xIonDense[nelem][ion_from] * (ion_from-ion_to);
01130 }
01131 }
01132 }
01133 }
01134 long ipMElec = findspecies("e-")->index;
01135 const double totsrc = ionsrc + mole.species[ipMElec].src;
01136 const double totsnk = ionsnk + mole.species[ipMElec].snk*dense.EdenTrue;
01137 const double diff = (totsrc - totsnk);
01138 const double ave = ( fabs(totsrc) + fabs(totsnk) )/2.;
01139 fixit();
01140 const double error_allowed = 0.05 * ave;
01141 if( fabs(diff) > error_allowed )
01142 {
01143 enum {DEBUG_LOC=false};
01144 if( DEBUG_LOC )
01145 {
01146 fprintf(ioQQQ,"PROBLEM large NetEdenSrc nzone %li\t%e\t%e\t%e\t%e\n",
01147 nzone,
01148 totsrc/SDIV(totsnk),
01149 dense.EdenTrue,
01150 ionsrc + mole.species[ipMElec].src,
01151 ionsnk + mole.species[ipMElec].snk*dense.EdenTrue);
01152 }
01153 conv.setConvIonizFail( "NetEdenSrc", diff, error_allowed);
01154 return false;
01155 }
01156 else
01157 return true;
01158 }