cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*ConvBase main routine to drive ionization solution for all species, find total opacity 00004 * called by ConvIoniz */ 00005 /*lgConverg check whether ionization of element nelem has converged */ 00006 #include "cddefines.h" 00007 #include "dynamics.h" 00008 #include "trace.h" 00009 #include "elementnames.h" 00010 #include "punch.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 /*lgIonizConverg check whether ionization of element nelem has converged, called by ionize, 00033 * returns true if element is converged, false if not */ 00034 STATIC bool lgIonizConverg( 00035 /* atomic number on the C scale, 0 for H, 25 for Fe */ 00036 long int nelem , 00037 /* this is allowed error as a fractional change. Most are 0.15 */ 00038 double delta , 00039 /* option to print abundances */ 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 /* this is fractions [ion stage][nelem], ion stage = 0 for atom, nelem=0 for H*/ 00051 static realnum OldFracs[LIMELM+1][LIMELM+1]; 00052 00053 if( !dense.lgElmtOn[nelem] ) 00054 return true; 00055 00056 DEBUG_ENTRY( "lgIonizConverg()" ); 00057 00058 /* arguments are atomic number, ionization stage 00059 * OldFracs[nelem][0] is abundance of nelem (cm^-3) */ 00060 00061 /* this function returns true if ionization of element 00062 * with atomic number nelem has not changed by more than delta,*/ 00063 00064 /* check whether abundances exist yet, only do this for after first zone */ 00065 /*if( nzone > 0 )*/ 00066 /* >>chng 03 sep 02, check on changed ionization after first call through, 00067 * to insure converged constant eden / temperature models */ 00068 if( conv.nPres2Ioniz ) 00069 { 00070 /* >>chng 04 aug 31, this had been static, caused last hits on large changes 00071 * in ionization */ 00072 bigchange = 0.; 00073 change = 0.; 00074 Abund = dense.gas_phase[nelem]; 00075 00076 /* loop over all ionization stages, loop over all ions, not just active ones, 00077 * since this also sets old values, and so will zero out non-existant ones */ 00078 for( ion=0; ion <= (nelem+1); ++ion ) 00079 { 00080 /*lint -e727 OlsdFracs not initialized */ 00081 if( OldFracs[nelem][ion]/Abund > 1e-4 && 00082 dense.xIonDense[nelem][ion]/Abund > 1e-4 ) 00083 /*lint +e727 OlsdFracs not initialized */ 00084 { 00085 /* check change in old to new value */ 00086 one = fabs(dense.xIonDense[nelem][ion]-OldFracs[nelem][ion])/ 00087 OldFracs[nelem][ion]; 00088 change = MAX2(change, one ); 00089 /* remember abundances for largest change */ 00090 if( change>bigchange ) 00091 { 00092 bigchange = change; 00093 abundold = OldFracs[nelem][ion]/Abund; 00094 abundnew = dense.xIonDense[nelem][ion]/Abund; 00095 } 00096 } 00097 /* now set new value */ 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 /* option to print abundances */ 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 /*ConvBase main routine to drive ionization solution for all species, find total opacity 00137 * called by ConvIoniz 00138 * return 0 if ok, 1 if abort */ 00139 int ConvBase( 00140 /* this is zero the first time ConvBase is called by convIoniz, 00141 * counts number of call thereafter */ 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 int nelem, 00156 ion, 00157 loop_ion, 00158 i, 00159 ipISO; 00160 static double SumOTS=0. , OldSumOTS[2]={0.,0.}; 00161 double save_iso_grnd[NISO][LIMELM]; 00162 00163 DEBUG_ENTRY( "ConvBase()" ); 00164 00165 #if 0 00166 if( loopi == 0 ) 00167 { 00168 OldSumOTS[0] = 0.; 00169 OldSumOTS[1] = 0.; 00170 SumOTS = 0.; 00171 nzoneOTS = -1; 00172 } 00173 #endif 00174 00175 /* this is set to phycon.te in tfidle, is used to insure that all temp 00176 * vars are properly updated when conv_ionizeopacitydo is called 00177 * NB must be same type as phycon.te */ 00178 ASSERT( fp_equal( phycon.te, thermal.te_update ) ); 00179 00180 /* this allows zone number to be printed with slight increment as zone converged 00181 * conv.nPres2Ioniz is incremented at the bottom of this routine */ 00182 fnzone = (double)nzone + (double)conv.nPres2Ioniz/100.; 00183 00184 /* reevaluate pressure */ 00185 /* this sets values of pressure.PresTotlCurr, also calls tfidle, 00186 * and sets the total energy content of gas, which may be important for acvection */ 00187 PresTotCurrent(); 00188 00189 /* >>chng 04 sep 15, move EdenTrue_old up here, and will redo at bottom 00190 * to find change 00191 * find and save current true electron density - if this changes by more than the 00192 * tolerance then ionization solution is not converged */ 00193 /* >>chng 04 jul 27, move eden_sum here from after this loop, so that change in EdenTrue 00194 * can be monitored */ 00195 /* update EdenTrue, eden itself is actually changed in ConvEdenIoniz */ 00196 /* must not call eden_sum on very first time since for classic PDR total 00197 * ionization may still be zero on first call */ 00198 if( conv.nTotalIoniz ) 00199 { 00200 if( eden_sum() ) 00201 { 00202 /* non-zero return indicates abort condition */ 00203 ++conv.nTotalIoniz; 00204 return 1; 00205 } 00206 } 00207 00208 /* the following two were evaluated in eden_sum 00209 * will confirm that these are converged */ 00210 EdenTrue_old = dense.EdenTrue; 00211 EdenFromMolecOld = co.comole_eden; 00212 EdenFromGrainsOld = gv.TotalEden; 00213 HeatingOld = thermal.htot; 00214 CoolingOld = thermal.ctot; 00215 00216 /* remember current ground state population - will check if converged */ 00217 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00218 { 00219 for( nelem=ipISO; nelem<LIMELM;++nelem ) 00220 { 00221 if( dense.lgElmtOn[nelem] ) 00222 { 00223 /* save the ground state population */ 00224 save_iso_grnd[ipISO][nelem] = StatesElem[ipISO][nelem][0].Pop; 00225 } 00226 } 00227 } 00228 00229 for( i=0; i < mole.num_comole_calc; i++ ) 00230 { 00231 COmole[i]->comole_save = COmole[i]->hevmol; 00232 } 00233 00234 for( i=0; i < N_H_MOLEC; i++ ) 00235 { 00236 hmole_save[i] = hmi.Hmolec[i]; 00237 } 00238 00239 if( loopi==0 ) 00240 { 00241 /* these will be used to look for oscillating ots rates */ 00242 OldSumOTS[0] = 0.; 00243 OldSumOTS[1] = 0.; 00244 } 00245 00246 if( trace.lgTrace ) 00247 { 00248 fprintf( ioQQQ, 00249 " ConvBase called. %.2f Te:%.3e HI:%.3e HII:%.3e H2:%.3e Ne:%.3e htot:%.3e CSUP:%.2e Conv?%c\n", 00250 fnzone, 00251 phycon.te, 00252 dense.xIonDense[ipHYDROGEN][0], 00253 dense.xIonDense[ipHYDROGEN][1], 00254 hmi.H2_total, 00255 dense.eden, 00256 thermal.htot, 00257 secondaries.csupra[ipHYDROGEN][0] , 00258 TorF(conv.lgConvIoniz) ); 00259 } 00260 /* want this flag to be true when we exit, various problems will set falst */ 00261 conv.lgConvIoniz = true; 00262 00263 /* this routine is in heatsum.c, and zeros out the heating array */ 00264 HeatZero(); 00265 00266 /* if this is very first call, say not converged, since no meaningful way to 00267 * check on changes in quantities. this counter is false only on very first 00268 * call, when equal to zero */ 00269 if( !conv.nTotalIoniz ) 00270 { 00271 conv.lgConvIoniz = false; 00272 strcpy( conv.chConvIoniz, "first call" ); 00273 } 00274 00275 /* this will be flag to check whether any ionization stages 00276 * were trimmed down */ 00277 conv.lgIonStageTrimed = false; 00278 00279 /* must redo photoionization rates for all species on second and later tries */ 00280 /* always reevaluate the rates when . . . */ 00281 /* this flag says not to redo opac and photo rates, and following test will set 00282 * true if one of several tests not done*/ 00283 opac.lgRedoStatic = false; 00284 if( 00285 /* opac.lgOpacStatic (usually true), is set false with no static opacity command */ 00286 !opac.lgOpacStatic || 00287 /* we are in search mode */ 00288 conv.lgSearch || 00289 /* this is the first call to this zone */ 00290 conv.nPres2Ioniz == 0 ) 00291 { 00292 /* we need to redo ALL photoionization rates */ 00293 opac.lgRedoStatic = true; 00294 } 00295 00296 /* calculate radiative and dielectronic recombination rate coefficients */ 00297 ion_recom_calculate(); 00298 00299 /* this adjusts the lowest and highest stages of ionization we will consider, 00300 * only safe to call when lgRedoStatic is true since this could lower the 00301 * lowest stage of ionization, which needs all its photo rates */ 00302 00303 /* conv.nTotalIoniz is only 0 (false) on the very first call to here, 00304 * when the ionization distribution is not yet done */ 00305 if( conv.nTotalIoniz ) 00306 { 00307 bool lgIonizTrimCalled = false; 00308 static long int nZoneCalled = 0; 00309 00310 fixit(); // nZoneCalled should be reinitialized for each grid point? 00311 00312 /* ionization trimming only used for He and heavier, not H */ 00313 /* only do this one time per zone since up and down cycle can occur */ 00314 /* >>chng 05 jan 15, increasing temperature above default first conditions, also 00315 * no trim during search - this fixed major logical error when sim is 00316 * totally mechanically heated to coronal temperatures - 00317 * problem discovered by Ronnie Hoogerwerf */ 00318 /* do not keep changing the trim after the first call within 00319 * this zone once we are deep into layer - doing so introduced very 00320 * small level of noise as some stages 00321 * went up and down - O+2 - O+3 was good example, when small H+3 after He+ i-front 00322 * limit to one increase per element per zone */ 00323 if( conv.nTotalIoniz>2 && 00324 /* only call one time per zone except during search phase, 00325 * when only call after 20 times only if temperature has changed */ 00326 ( conv.lgSearch || nZoneCalled!=nzone) ) 00327 { 00328 lgIonizTrimCalled = true; 00329 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem ) 00330 { 00331 if( dense.lgElmtOn[nelem] ) 00332 { 00333 /* ion_trim will set conv.lgIonStageTrimed true is any ion has its 00334 * lowest stage of ionization dropped or trimmed */ 00335 ion_trim(nelem); 00336 } 00337 } 00338 nZoneCalled = nzone; 00339 } 00340 00341 /* following block only set of asserts */ 00342 # if !defined(NDEBUG) 00343 /* check that proper abundances are either positive or zero */ 00344 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem) 00345 { 00346 if( dense.lgElmtOn[nelem] ) 00347 { 00348 for( ion=0; ion<dense.IonLow[nelem]; ++ion ) 00349 { 00350 ASSERT( dense.xIonDense[nelem][ion] == 0. ); 00351 } 00352 /*if( nelem==5 ) fprintf(ioQQQ,"carbbb\t%li\n", dense.IonHigh[nelem]);*/ 00353 for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion ) 00354 { 00355 /* >>chng 02 feb 06, had been > o., chng to > SMALLFLOAT to 00356 * trip over VERY small floats that failed on alphas, but not 386 00357 * 00358 * in case where lower ionization stage was just lowered or 00359 * trimmed down the abundance 00360 * was set to SMALLFLOAT so test must be < SMALLFLOAT */ 00361 /* >>chng 02 feb 19, add check for search phase. During this search 00362 * models with extreme ionization (all neutral or all ionized) can 00363 * have extreme but non-zero abundances far from the ionization peak for 00364 * element with lots of electrons. These will go away once the model 00365 * becomes stable */ 00366 /* >>chng 03 dec 01, add check on whether ion trim was called 00367 * conserve.in threw assert when iontrim not called and abund grew small */ 00368 ASSERT( conv.lgSearch || !lgIonizTrimCalled || 00369 /* this can happen if all C is in the form of CO */ 00370 (ion==0 && dense.IonHigh[nelem]==0 ) || 00371 dense.xIonDense[nelem][ion] >= SMALLFLOAT || 00372 dense.xIonDense[nelem][ion]/dense.gas_phase[nelem] >= SMALLFLOAT ); 00373 } 00374 for( ion=dense.IonHigh[nelem]+1; ion<nelem+1; ++ion ) 00375 { 00376 ASSERT( ion >= 0 ); 00377 ASSERT( dense.xIonDense[nelem][ion] == 0. ); 00378 } 00379 } 00380 } 00381 # endif 00382 } 00383 00384 /* now check if anything trimmed down */ 00385 if( conv.lgIonStageTrimed ) 00386 { 00387 /* something was trimmed down, so say that ionization not yet stable */ 00388 /* say that ionization has not converged, secondaries changed by too much */ 00389 conv.lgConvIoniz = false; 00390 strcpy( conv.chConvIoniz, "IonTrimmed" ); 00391 } 00392 00393 /* reevaluate advective terms if turned on */ 00394 if( dynamics.lgAdvection ) 00395 DynaIonize(); 00396 00397 /* >>chng 04 feb 15, add loop over ionization until converged. Non-convergence 00398 * in this case means ionization ladder pop changed, probably because of way 00399 * that Auger ejection is treated - loop exits when conditions tested at end are met */ 00400 loop_ion = 0; 00401 do 00402 { 00403 00404 if( trace.nTrConvg>=4 ) 00405 { 00406 /* trace ionization */ 00407 fprintf( ioQQQ, 00408 " ConvBase4 ionization driver loop_ion %li converged? %c reason not converged %s\n" , 00409 loop_ion , 00410 TorF( conv.lgConvIoniz) , 00411 conv.chConvIoniz ); 00412 } 00413 00414 /* >>chng 04 sep 11, moved to before grains and other routines - must not clear this 00415 * flag before their calculations are done */ 00416 /* set the convergence flag to true, 00417 * if anything changes in ConvBase, it will be set false */ 00418 if( loop_ion ) 00419 conv.lgConvIoniz = true; 00420 00421 /* >>chng 01 apr 29, move charge transfer evaluation here, had been just before 00422 * HeLike, but needs to be here so that same rate coefficient used for H ion and other recomb */ 00423 /* fill in master charge transfer array, and zero out some arrays that track effects, 00424 * O_HIonRateis rate oxygen ionizes hydrogen - will need to converge this */ 00425 ChargTranEval( &O_HIonRate_Old ); 00426 00427 CO_update_species_cache(); /* Update densities of species controlled outside the chemical network */ 00428 00429 hmole_reactions(); 00430 00431 co.nitro_dissoc_rate = 0; 00432 00433 /* Update the rate constants -- generic version */ 00434 CO_update_rks(); 00435 00436 /* find grain temperature, charge, and gas heating rate */ 00437 /* >>chng 04 apr 15, moved here from after call to HeLike(), PvH */ 00438 GrainDrive(); 00439 00440 /* do all hydrogenic species, and fully do hydrogen itself, including molecules */ 00441 /*fprintf(ioQQQ,"DEBUG h2\t%.2f\t%.3e\t%.3e", fnzone,hmi.H2_total,hmi.Hmolec[ipMH2g]);*/ 00442 iso_solve( ipH_LIKE ); 00443 /*fprintf(ioQQQ,"\t%.3e\n", hmi.Hmolec[ipMH2g]);*/ 00444 00445 /* evaluate Compton heating, bound E Compton, cosmic rays */ 00446 highen(); 00447 00448 /* find corrections for three-body rec - collisional ionization */ 00449 atmdat_3body(); 00450 00451 /* deduce dielectronic suppression factors */ 00452 atmdat_DielSupres(); 00453 00454 /* >>chng 02 mar 08 move helike to before helium */ 00455 /* do all he-like species */ 00456 iso_solve( ipHE_LIKE ); 00457 00458 /*>>chng 04 may 09, add option to abort here, inspired by H2 pop failures 00459 * which cascaded downstream if we did not abort */ 00460 /* return if one of above solvers has declared abort condition */ 00461 if( lgAbort ) 00462 { 00463 ++conv.nTotalIoniz; 00464 return 1; 00465 } 00466 00467 /* find grain temperature, charge, and gas heating rate */ 00468 /*GrainDrive();*/ 00469 00470 /* only reevaluate this on first pass through on this zone */ 00471 if( opac.lgRedoStatic ) 00472 { 00473 /* inner shell ionization */ 00474 for( nelem=0; nelem< LIMELM; ++nelem ) 00475 { 00476 for( ion=0; ion<nelem+1; ++ion ) 00477 { 00478 ionbal.UTA_ionize_rate[nelem][ion] = 0.; 00479 ionbal.UTA_heat_rate[nelem][ion] = 0.; 00480 } 00481 } 00482 /* inner shell ionization by UTA lines */ 00483 /* this flag turned off with no UTA command */ 00484 if( ionbal.lgInnerShellLine_on ) 00485 { 00486 for( i=0; i<nUTA; ++i ) 00487 { 00488 if( UTALines[i].Emis->Aul > 0. ) 00489 { 00490 /* cs is actually the negative of the branching ratio for autoionization, 00491 * rateone is inverse lifetime of level against autoionization */ 00492 double rateone = UTALines[i].Emis->pump * -UTALines[i].Coll.cs; 00493 ionbal.UTA_ionize_rate[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1] += rateone; 00494 /* heating rate in erg atom-1 s-1 */ 00495 ionbal.UTA_heat_rate[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1] += rateone*UTALines[i].Coll.heat; 00496 { 00497 /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */ 00498 /*@-redef@*/ 00499 enum {DEBUG_LOC=false}; 00500 /*@+redef@*/ 00501 if( DEBUG_LOC /*&& UTALines[i].nelem==ipIRON+1 && (UTALines[i].IonStg==15||UTALines[i].IonStg==14)*/ ) 00502 { 00503 fprintf(ioQQQ,"DEBUG UTA %3i %3i %.3f %.2e %.2e %.2e\n", 00504 UTALines[i].Hi->nelem , UTALines[i].Hi->IonStg , UTALines[i].WLAng , 00505 rateone, UTALines[i].Coll.heat, 00506 UTALines[i].Coll.heat*dense.xIonDense[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1] ); 00507 } 00508 } 00509 } 00510 } 00511 } 00512 } 00513 00514 /* helium ionization balance */ 00515 IonHelium(); 00516 00517 /* do the ionization balance */ 00518 /* evaluate current opacities, OpacityAddTotal is called only here during actual calculation */ 00519 /* option to only evaluate opacity one time per zone, 00520 * rfield.lgOpacityReevaluate normally true, set false with no opacity reevaluate command */ 00521 /* >>chng 04 jul 16, move these out of block below, so always done 00522 * these all produce ions that are used in the CO network */ 00523 IonCarbo(); 00524 IonOxyge(); 00525 IonNitro(); 00526 IonSilic(); 00527 IonSulph(); 00528 IonChlor(); 00529 /* do carbon monoxide after oxygen */ 00530 CO_drive(); 00531 if( !conv.nPres2Ioniz || rfield.lgIonizReevaluate || 00532 conv.lgIonStageTrimed || conv.lgSearch ) 00533 { 00534 IonLithi(); 00535 IonBeryl(); 00536 IonBoron(); 00537 /*IonCarbo(); 00538 IonOxyge();*/ 00539 /* do carbon monoxide after oxygen */ 00540 /*fprintf(ioQQQ,"DEBUG co\t%.2f\t%.3e", fnzone,findspecies("CO")->hevmol);*/ 00541 /*CO_drive();*/ 00542 /*fprintf(ioQQQ,"\t%.3e\n", findspecies("CO")->hevmol);*/ 00543 IonFluor(); 00544 IonNeon(); 00545 IonSodiu(); 00546 IonMagne(); 00547 IonAlumi(); 00548 IonPhosi(); 00549 IonArgon(); 00550 IonPotas(); 00551 IonCalci(); 00552 IonScand(); 00553 IonTitan(); 00554 IonVanad(); 00555 IonChrom(); 00556 IonManga(); 00557 IonIron(); 00558 IonCobal(); 00559 IonNicke(); 00560 IonCoppe(); 00561 IonZinc(); 00562 } 00563 00564 /* all elements have now had ionization reevaluated - in some cases we may have upset 00565 * the ionization that was forced with an "element ionization" command - here we will 00566 * re-impose that set ionization */ 00567 /* >>chng 04 oct 13, add this logic */ 00568 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ ) 00569 { 00570 if( dense.lgSetIoniz[nelem] ) 00571 { 00572 dense.IonLow[nelem] = 0; 00573 dense.IonHigh[nelem] = nelem + 1; 00574 while( dense.SetIoniz[nelem][dense.IonLow[nelem]] == 0. ) 00575 ++dense.IonLow[nelem]; 00576 while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] == 0. ) 00577 --dense.IonHigh[nelem]; 00578 } 00579 } 00580 00581 /* redo Oxygen ct ion rate of H to see if it changed */ 00582 ChargTranEval( &O_HIonRate_New ); 00583 00584 /* lgIonizConverg is a function to check whether ionization has converged 00585 * check whether ionization changed by more than relative error 00586 * given by second number */ 00587 /* >>chng 04 feb 14, loop over all elements rather than just a few */ 00588 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem ) 00589 { 00590 if( !lgIonizConverg(nelem, 0.05 ,false ) ) 00591 { 00592 conv.lgConvIoniz = false; 00593 sprintf( conv.chConvIoniz , "%2s ion" , elementnames.chElementSym[nelem] ); 00594 } 00595 } 00596 00597 if( fabs(EdenTrue_old - dense.EdenTrue)/SDIV(dense.EdenTrue) > conv.EdenErrorAllowed/2. ) 00598 { 00599 conv.lgConvIoniz = false; 00600 sprintf( conv.chConvIoniz , "EdTr cng" ); 00601 conv.BadConvIoniz[0] = EdenTrue_old; 00602 conv.BadConvIoniz[1] = dense.EdenTrue; 00603 } 00604 ++loop_ion; 00605 } 00606 /* loop is not converged, less than loop limit, and we are reevaluating */ 00607 while( !conv.lgConvIoniz && loop_ion < LOOP_ION_LIMIT && rfield.lgIonizReevaluate); 00608 00609 /* >>chng 05 oct 29, move CT heating here from heat_sum since this sometimes contributes 00610 * cooling rather than heat and so needs to be sorted out before either heating or cooling 00611 * are derived first find net heating - */ 00612 thermal.char_tran_heat = ChargTranSumHeat(); 00613 /* net is cooling if negative */ 00614 thermal.char_tran_cool = MAX2(0. , -thermal.char_tran_heat ); 00615 thermal.char_tran_heat = MAX2(0. , thermal.char_tran_heat ); 00616 00617 /* get total cooling, thermal.ctot = does not occur since passes as pointer. This can add heat. 00618 * it calls coolSum at end to sum up the total cooling */ 00619 CoolEvaluate(&thermal.ctot ); 00620 00621 /* get total heating rate - first save old quantities to check how much it changes */ 00622 HeatOld = thermal.htot; 00623 00624 /* HeatSum will update ElecFrac, 00625 * secondary electron ionization and excitation efficiencies, 00626 * and sum up current secondary rates - remember old values before we enter */ 00627 SecondOld = secondaries.csupra[ipHYDROGEN][0]; 00628 00629 /* update the total heating - it was all set to zero in HeatZero at top of this routine, 00630 * occurs before secondaries bit below, since this updates electron fracs */ 00631 HeatSum(); 00632 00633 /* test whether we can just set the rate to the new one, or whether we should 00634 * take average and do it again. secondaries.sec2total was set in hydrogenic, and 00635 * is ratio of secondary to total hydrogen destruction rates */ 00636 /* >>chng 02 nov 20, add test on size of csupra - primal had very close to underflow */ 00637 if( (secondaries.csupra[ipHYDROGEN][0]>SMALLFLOAT) && (secondaries.sec2total>0.001) && 00638 fabs( 1. - SecondOld/SDIV(secondaries.csupra[ipHYDROGEN][0]) ) > 0.1 && 00639 SecondOld > 0. && secondaries.csupra[ipHYDROGEN][0] > 0.) 00640 { 00641 /* say that ionization has not converged, secondaries changed by too much */ 00642 conv.lgConvIoniz = false; 00643 strcpy( conv.chConvIoniz, "SecIonRate" ); 00644 conv.BadConvIoniz[0] = SecondOld; 00645 conv.BadConvIoniz[1] = secondaries.csupra[ipHYDROGEN][0]; 00646 } 00647 00648 # if 0 00649 static realnum hminus_old=0.; 00650 /* >>chng 04 apr 15, add this convergence test */ 00651 if( conv.nTotalIoniz ) 00652 { 00653 if( fabs( hminus_old-hmi.Hmolec[ipMHm])/SDIV(hmi.Hmolec[ipMHm]) > 00654 conv.EdenErrorAllowed/4. ) 00655 conv.lgConvIoniz = false; 00656 strcpy( conv.chConvIoniz, "Big H- chn" ); 00657 conv.BadConvIoniz[0] = hminus_old; 00658 conv.BadConvIoniz[1] = hmi.Hmolec[ipMHm]; 00659 } 00660 hminus_old = hmi.Hmolec[ipMHm]; 00661 # endif 00662 00663 if( HeatOld > 0. && !thermal.lgTemperatureConstant ) 00664 { 00665 /* check if heating has converged - tolerance is final match */ 00666 if( fabs(1.-thermal.htot/HeatOld) > conv.HeatCoolRelErrorAllowed*.5 ) 00667 { 00668 conv.lgConvIoniz = false; 00669 strcpy( conv.chConvIoniz, "Big d Heat" ); 00670 conv.BadConvIoniz[0] = HeatOld; 00671 conv.BadConvIoniz[1] = thermal.htot; 00672 } 00673 } 00674 00675 /* abort flag may have already been set - if so bail */ 00676 if( lgAbort ) 00677 { 00678 00679 return 1; 00680 } 00681 00682 /* evaluate current opacities, OpacityAddTotal is called only here during actual calculation */ 00683 /* option to only evaluate opacity one time per zone, 00684 * rfield.lgOpacityReevaluate normally true, set false with no opacity reevaluate command */ 00685 if( !conv.nPres2Ioniz || rfield.lgOpacityReevaluate || conv.lgIonStageTrimed ) 00686 OpacityAddTotal(); 00687 00688 /* >>chng 02 jun 11, call even first time that this routine is called - 00689 * this seemed to help convergence */ 00690 00691 /* do OTS rates for all lines and all continua since 00692 * we now have ionization balance of all species. Note that this is not 00693 * entirely self-consistent, since destruction probabilities here are not the same as 00694 * the ones used in the model atoms. Problems?? if near convergence 00695 * then should be nearly identical */ 00696 if( !conv.nPres2Ioniz || rfield.lgOpacityReevaluate || rfield.lgIonizReevaluate || 00697 conv.lgIonStageTrimed || conv.lgSearch || nzone!=nzoneOTS ) 00698 { 00699 RT_OTS(); 00700 nzoneOTS = nzone; 00701 00702 /* remember old ots rates */ 00703 OldSumOTS[0] = OldSumOTS[1]; 00704 OldSumOTS[1] = SumOTS; 00705 /*fprintf(ioQQQ," calling RT_OTS zone %.2f SumOTS is %.2e\n",fnzone,SumOTS);*/ 00706 00707 /* now update several components of the continuum, this only happens after 00708 * we have gone through entire solution for this zone at least one time. 00709 * there can be wild ots oscillation on first call */ 00710 /* the rel change of 0.2 was chosen by running hizqso - smaller increased 00711 * itrzn but larger did not further decrease it. */ 00712 RT_OTS_Update(&SumOTS , 0.20 ); 00713 /*fprintf(ioQQQ,"RT_OTS_Updateee\t%.3f\t%.2e\t%.2e\n", fnzone,SumOTS , OldSumOTS[1] );*/ 00714 } 00715 else 00716 SumOTS = 0.; 00717 00718 /* now check whether the ots rates changed */ 00719 if( SumOTS> 0. ) 00720 { 00721 /* the ots rate must be converged to the error in the electron density */ 00722 /* how fine should this be converged?? originally had above, 10%, but take 00723 * smaller ratio?? */ 00724 if( fabs(1.-OldSumOTS[1]/SumOTS) > conv.EdenErrorAllowed ) 00725 { 00726 /* this branch, ionization not converged due to large change in ots rates. 00727 * check whether ots rates are oscillating, if loopi > 1 so we have enough info*/ 00728 if( loopi > 1 ) 00729 { 00730 /* here we have three values, are they changing sign? */ 00731 if( (OldSumOTS[0]-OldSumOTS[1]) * ( OldSumOTS[1] - SumOTS ) < 0. ) 00732 { 00733 /* ots rates are oscillating, if so then use small fraction of new destruction rates 00734 * when updating Lyman line destruction probabilities in HydroPesc*/ 00735 conv.lgOscilOTS = true; 00736 } 00737 } 00738 00739 conv.lgConvIoniz = false; 00740 strcpy( conv.chConvIoniz, "OTSRatChng" ); 00741 conv.BadConvIoniz[0] = OldSumOTS[1]; 00742 conv.BadConvIoniz[1] = SumOTS; 00743 } 00744 /* produce info on the ots fields if either "trace ots" or 00745 * "trace convergence xxx ots " was entered */ 00746 if( ( trace.lgTrace && trace.lgOTSBug ) || 00747 ( trace.nTrConvg && trace.lgOTSBug ) ) 00748 { 00749 RT_OTS_PrtRate(SumOTS*0.05 , 'b' ); 00750 } 00751 /*fprintf(ioQQQ,"DEBUG opac\t%.2f\t%.3e\t%.3e\n",fnzone, 00752 dense.xIonDense[ipNICKEL][0] , 00753 dense.xIonDense[ipZINC][0] );*/ 00754 { 00755 /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */ 00756 /*@-redef@*/ 00757 enum {DEBUG_LOC=false}; 00758 /*@+redef@*/ 00759 if( DEBUG_LOC && (nzone>110) ) 00760 { 00761 # if 0 00762 # include "lines_service.h" 00763 DumpLine( &Transitions[ipH_LIKE][ipHYDROGEN][2][0] ); 00764 # endif 00765 /* last par 'l' for line, 'c' for continua, 'b' for both, 00766 * the numbers printed are: 00767 * cell i, [i], so 1 less than ipoint 00768 * anu[i], 00769 * otslin[i], 00770 * opacity_abs[i], 00771 * otslin[i]*opacity_abs[i], 00772 * rfield.chLineLabel[i] , 00773 * rfield.line_count[i] */ 00774 } 00775 } 00776 } 00777 { 00778 /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */ 00779 /*@-redef@*/ 00780 enum {DEBUG_LOC=false}; 00781 /*@+redef@*/ 00782 if( DEBUG_LOC && (nzone>200) ) 00783 { 00784 fprintf(ioQQQ,"debug otsss\t%li\t%.3e\t%.3e\t%.3e\n", 00785 nzone, 00786 Transitions[0][1][15][3].Emis->ots, 00787 TauLines[26].Emis->ots, 00788 opac.opacity_abs[2069]); 00789 } 00790 } 00791 00792 /* option to print OTS continuum with TRACE OTS */ 00793 if( trace.lgTrace && trace.lgOTSBug ) 00794 { 00795 /* find ots rates here, so we only print fraction of it, 00796 * SumOTS is both line and continuum contributing to ots, and is multiplied by opacity */ 00797 /* number is weakest rate to print */ 00798 RT_OTS_PrtRate(SumOTS*0.05 , 'b' ); 00799 } 00800 00801 /* now reevaluate only destruction probabilities if call is not first*/ 00802 /* >>chng 02 jun 13, had always been false, now reevaluate escape probabilities 00803 * and pumping rates on first call for each zone */ 00804 if( conv.nPres2Ioniz && !conv.lgSearch ) 00805 { 00806 RT_line_all(false , false ); 00807 } 00808 else 00809 { 00810 RT_line_all(true , false ); 00811 } 00812 00813 /* >>chng 01 mar 16, evaluate pressure here since changing and other values needed */ 00814 /* reevaluate pressure */ 00815 /* this sets values of pressure.PresTotlCurr, also calls tfidle */ 00816 PresTotCurrent(); 00817 00818 if( trace.lgTrace ) 00819 { 00820 fprintf( ioQQQ, 00821 " ConvBase return. %.2f Te:%.3e HI:%.3e HII:%.3e H2:%.3e Ne:%.3e htot:%.3e CSUP:%.2e Conv?%c\n", 00822 fnzone, 00823 phycon.te, 00824 dense.xIonDense[ipHYDROGEN][0], 00825 dense.xIonDense[ipHYDROGEN][1], 00826 hmi.H2_total, 00827 dense.eden, 00828 thermal.htot, 00829 secondaries.csupra[ipHYDROGEN][0] , 00830 TorF(conv.lgConvIoniz) ); 00831 } 00832 00833 /* this checks that all molecules and ions add up to gas phase abundance 00834 * but only if we have a converged solution */ 00835 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00836 { 00837 /* check that species add up if element is turned on, and 00838 * ionization is not set */ 00839 if( dense.lgElmtOn[nelem] && !dense.lgSetIoniz[nelem]) 00840 { 00841 /* this sum of over the molecules did not include the atom and first 00842 * ion that is calculated in the molecular solver */ 00843 double sum = dense.xMolecules[nelem]; 00844 for( ion=0; ion<nelem+2; ++ion ) 00845 { 00846 sum += dense.xIonDense[nelem][ion]; 00847 } 00848 ASSERT( sum > SMALLFLOAT ); 00849 /* check that sum agrees with total gas phase abundance */ 00853 /*if( fabs(sum-dense.gas_phase[nelem])/sum > 1e-2 ) 00854 { 00855 fprintf(ioQQQ,"DEBUG non conservation element %li sum %.3e gas phase %.3e rel error %.3f\n", 00856 nelem, sum , dense.gas_phase[nelem],(sum-dense.gas_phase[nelem])/sum); 00857 }*/ 00858 /* >>chng 04 jul 23, needed to increase to 2e-2 from 1e-2 to get conserve.in to work, 00859 * not clear what caused increase in error */ 00860 if( fabs(sum-dense.gas_phase[nelem])/SDIV(sum) > 2e-2 ) 00861 { 00862 /*fprintf(ioQQQ,"DEBUG non-conv nelem\t%li\tsum\t%e\ttot gas\t%e\trel err\t%.3e\n", 00863 nelem, 00864 sum, 00865 dense.gas_phase[nelem], 00866 (sum-dense.gas_phase[nelem])/sum);*/ 00867 conv.lgConvIoniz = false; 00868 strcpy( conv.chConvIoniz, "CO con4" ); 00869 sprintf( conv.chConvIoniz, "COnelem%2li",nelem ); 00870 conv.BadConvIoniz[0] = sum; 00871 conv.BadConvIoniz[1] = dense.gas_phase[nelem]; 00872 } 00873 } 00874 } 00875 00876 /* update some counters that keep track of how many times this routine 00877 * has been called */ 00878 00879 /* this counts number of times we are called by ConvPresTempEdenIoniz, 00880 * number of calls in this zone so first call is zero 00881 * reset to zero each time ConvPresTempEdenIoniz is called */ 00882 ++conv.nPres2Ioniz; 00883 00884 /* this is abort option set with SET PRESIONIZ command, 00885 * test on nzone since init can take many iterations 00886 * this is seldom used except in special cases */ 00887 if( conv.nPres2Ioniz > conv.limPres2Ioniz && nzone > 0) 00888 { 00889 fprintf(ioQQQ,"PROBLEM ConvBase sets lgAbort since nPres2Ioniz exceeds limPres2Ioniz.\n"); 00890 fprintf(ioQQQ,"Their values are %li and %li \n",conv.nPres2Ioniz , conv.limPres2Ioniz); 00891 lgAbort = true; 00892 00893 return 1; 00894 } 00895 00896 /* various checks on the convergence of the current solution */ 00897 /* >>chng 04 sep 15, add this check if the ionization so converged */ 00898 if( eden_sum() ) 00899 { 00900 /* non-zero return indicates abort condition */ 00901 return 1; 00902 } 00903 /* is electron density converged? */ 00905 if( fabs(EdenTrue_old - dense.EdenTrue)/dense.EdenTrue > conv.EdenErrorAllowed/2. ) 00906 { 00907 conv.lgConvIoniz = false; 00908 sprintf( conv.chConvIoniz, "eden chng" ); 00909 conv.BadConvIoniz[0] = EdenTrue_old; 00910 conv.BadConvIoniz[1] = dense.EdenTrue; 00911 } 00912 00913 /* check on molecular electron den */ 00914 if( fabs(EdenFromMolecOld - co.comole_eden)/dense.EdenTrue > conv.EdenErrorAllowed/2. ) 00915 { 00916 conv.lgConvIoniz = false; 00917 sprintf( conv.chConvIoniz, "edn chnCO" ); 00918 conv.BadConvIoniz[0] = EdenFromMolecOld; 00919 conv.BadConvIoniz[1] = dense.EdenTrue; 00920 } 00921 00922 /* check on molecular electron den */ 00923 if( fabs(EdenFromGrainsOld - gv.TotalEden)/dense.EdenTrue > conv.EdenErrorAllowed/2. ) 00924 { 00925 conv.lgConvIoniz = false; 00926 sprintf( conv.chConvIoniz, "edn grn e" ); 00927 conv.BadConvIoniz[0] = EdenFromGrainsOld; 00928 conv.BadConvIoniz[1] = gv.TotalEden; 00929 } 00930 00931 /* check on heating and cooling if vary temp model */ 00932 if( !thermal.lgTemperatureConstant ) 00933 { 00934 if( fabs(HeatingOld - thermal.htot)/thermal.htot > conv.HeatCoolRelErrorAllowed/2. ) 00935 { 00936 conv.lgConvIoniz = false; 00937 sprintf( conv.chConvIoniz, "heat chg" ); 00938 conv.BadConvIoniz[0] = HeatingOld; 00939 conv.BadConvIoniz[1] = thermal.htot; 00940 } 00941 00942 if( fabs(CoolingOld - thermal.ctot)/thermal.ctot > conv.HeatCoolRelErrorAllowed/2. ) 00943 { 00944 conv.lgConvIoniz = false; 00945 sprintf( conv.chConvIoniz, "cool chg" ); 00946 conv.BadConvIoniz[0] = CoolingOld; 00947 conv.BadConvIoniz[1] = thermal.ctot; 00948 } 00949 } 00950 00951 /* check on sum of grain and molecular electron den - often two large numbers that cancel */ 00952 if( fabs( (EdenFromMolecOld-EdenFromGrainsOld) - (co.comole_eden-gv.TotalEden) )/dense.eden > 00953 conv.EdenErrorAllowed/4. ) 00954 { 00955 conv.lgConvIoniz = false; 00956 sprintf( conv.chConvIoniz, "edn grn e" ); 00957 conv.BadConvIoniz[0] = (EdenFromMolecOld-EdenFromGrainsOld); 00958 conv.BadConvIoniz[1] = (co.comole_eden-gv.TotalEden); 00959 } 00960 00961 /* check whether molecular abundances are stable - these were evaluated in CO_drive */ 00962 for( i=0; i < mole.num_comole_calc; ++i ) 00963 { 00964 if( dense.lgElmtOn[COmole[i]->nelem_hevmol] && COmole[i]->hevmol>SMALLFLOAT ) 00965 { 00966 if( fabs(COmole[i]->hevmol-COmole[i]->comole_save)/dense.gas_phase[COmole[i]->nelem_hevmol]-1. > 00967 conv.EdenErrorAllowed/2.) 00968 { 00969 conv.lgConvIoniz = false; 00970 sprintf( conv.chConvIoniz, "ch %s",COmole[i]->label ); 00971 conv.BadConvIoniz[0] = COmole[i]->comole_save/dense.gas_phase[COmole[i]->nelem_hevmol]; 00972 conv.BadConvIoniz[1] = COmole[i]->hevmol/dense.gas_phase[COmole[i]->nelem_hevmol]; 00973 } 00974 } 00975 } 00976 00977 /* check whether H molecular abundances are stable */ 00978 for( i=0; i < N_H_MOLEC; ++i ) 00979 { 00980 if( fabs(hmi.Hmolec[i]-hmole_save[i])/dense.gas_phase[ipHYDROGEN]-1. > 00981 conv.EdenErrorAllowed/2.) 00982 { 00983 conv.lgConvIoniz = false; 00984 sprintf( conv.chConvIoniz, "ch %s",hmi.chLab[i] ); 00985 conv.BadConvIoniz[0] = hmole_save[i]/dense.gas_phase[ipHYDROGEN]; 00986 conv.BadConvIoniz[1] = hmi.Hmolec[i]/dense.gas_phase[ipHYDROGEN]; 00987 } 00988 } 00989 00990 /* >>chng 05 mar 26, add this convergence test - important for molecular or advective 00991 * sims since iso ion solver must sync up with chemistry */ 00992 /* remember current ground state population - will check if converged */ 00993 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00994 { 00995 for( nelem=ipISO; nelem<LIMELM;++nelem ) 00996 { 00997 if( dense.lgElmtOn[nelem] ) 00998 { 00999 /* only do this check for "significant" levels of ionization */ 01000 /*lint -e644 var possibly not init */ 01001 if( dense.xIonDense[nelem][nelem+1-ipISO]/dense.gas_phase[nelem] > 1e-5 ) 01002 { 01003 if( fabs(StatesElem[ipISO][nelem][0].Pop-save_iso_grnd[ipISO][nelem])/SDIV(StatesElem[ipISO][nelem][0].Pop)-1. > 01004 conv.EdenErrorAllowed) 01005 { 01006 conv.lgConvIoniz = false; 01007 sprintf( conv.chConvIoniz,"iso %3li%li",ipISO, nelem ); 01008 conv.BadConvIoniz[0] = save_iso_grnd[ipISO][nelem]/SDIV(StatesElem[ipISO][nelem][0].Pop); 01009 fixit(); // this looks like a bug, it Pop/Pop, unity for all Pop>SMALLFLOAT 01010 conv.BadConvIoniz[1] = StatesElem[ipISO][nelem][0].Pop/SDIV(StatesElem[ipISO][nelem][0].Pop); 01011 } 01012 } 01013 /*lint +e644 var possibly not init */ 01014 } 01015 } 01016 } 01017 01018 /* this counts how many times ConvBase has been called in this iteration, 01019 * located here at bottom of routine so that number is false on first 01020 * call, set to 0 in when iteration starts - used to create itr/zn 01021 * number in printout often used to tell whether this is the very 01022 * first attempt at solution in this iteration */ 01023 ++conv.nTotalIoniz; 01024 01025 /* this was set with STOP NTOTALIONIZ command - only a debugging aid 01026 * by default is zero so false */ 01027 if(StopCalc.nTotalIonizStop && conv.nTotalIoniz> StopCalc.nTotalIonizStop ) 01028 { 01029 /* non-zero return indicates abort condition */ 01030 fprintf(ioQQQ , "ABORT flag set since STOP nTotalIoniz was set and reached.\n"); 01031 return 1; 01032 } 01033 01034 if( punch.lgTraceConvergeBase ) 01035 { 01036 if( iteration > 1 && punch.lgTraceConvergeBaseHash ) 01037 { 01038 static int iter_punch=-1; 01039 if( iteration !=iter_punch ) 01040 fprintf( punch.ipDRout, "%s\n",punch.chHashString ); 01041 iter_punch = iteration; 01042 } 01043 01044 fprintf( punch.ipTraceConvergeBase, 01045 "%li\t%.4e\t%.4e\t%.4e\n", 01046 nzone , thermal.htot , thermal.ctot , dense.eden ); 01047 } 01048 01049 return 0; 01050 }