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