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 /*ion_trim raise or lower most extreme stages of ionization considered, 00004 * called by ConvBase - ion limits were originally set by */ 00005 #include "cddefines.h" 00006 #include "elementnames.h" 00007 #include "radius.h" 00008 #include "heavy.h" 00009 #include "conv.h" 00010 #include "rfield.h" 00011 #include "phycon.h" 00012 #include "mole.h" 00013 #include "thermal.h" 00014 #include "iso.h" 00015 #include "dense.h" 00016 #include "conv.h" 00017 #include "struc.h" 00018 #include "ionbal.h" 00019 00020 void ion_trim( 00021 /* nelem is on the C scale, 0 for H, 5 for C */ 00022 long int nelem ) 00023 { 00024 00025 /* this will remember that higher stages trimed up or down */ 00026 bool lgHi_Down = false; 00027 bool lgHi_Up = false; 00028 bool lgHi_Up_enable; 00029 /* this will remember that lower stages trimmed up or own*/ 00030 bool lgLo_Up = false; 00031 bool lgLo_Down = false; 00032 long int ion_lo_save = dense.IonLow[nelem], 00033 ion_hi_save = dense.IonHigh[nelem]; 00034 long int ion; 00035 realnum trimhi , trimlo; 00036 00037 /* this is debugging code that turns on print after certain number of calls */ 00038 /*realnum xlimit = SMALLFLOAT *10.;*/ 00039 /*static int ncall=0; 00040 if( nelem==5 ) 00041 ++ncall;*/ 00042 00043 DEBUG_ENTRY( "ion_trim()" ); 00044 00045 /* this routine should not be called early in the simulation, before 00046 * things have settled down */ 00047 ASSERT( conv.nTotalIoniz>2 ); 00048 00049 /*confirm all ionization stages are within their range of validity */ 00050 ASSERT( nelem >= ipHELIUM && nelem < LIMELM ); 00051 ASSERT( dense.IonLow[nelem] >= 0 ); 00052 ASSERT( dense.IonHigh[nelem] <= nelem+1 ); 00053 /* IonHigh can be equal to IonLow if both are zero - so no ionization, 00054 * or is ionization has been set with "element ionization" command */ 00055 ASSERT( dense.IonLow[nelem] < dense.IonHigh[nelem] || 00056 (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) || 00057 dense.lgSetIoniz[nelem] ); 00058 00059 /* during search phase of mostly neutral matter the electron density 00060 * can be vastly too large, and the ionization suppressed as a result. 00061 * during search do not trim down or up as much */ 00062 if( conv.lgSearch ) 00063 { 00064 trimhi = (realnum)ionbal.trimhi * 1e-4f; 00065 trimlo = (realnum)ionbal.trimlo * 1e-4f; 00066 } 00067 else 00068 { 00069 trimhi = (realnum)ionbal.trimhi; 00070 trimlo = (realnum)ionbal.trimlo; 00071 } 00072 00073 /* helium is special case since abundance is so high, and He+ CT with 00074 * CO is the dominant CO destruction process in molecular regions */ 00075 if( nelem == ipHELIUM ) 00076 { 00077 /* never want to trip up a lower stage of ionization */ 00078 trimlo = SMALLFLOAT; 00079 00080 /* if He+ is highest stage of ionization, probably want to keep it 00081 * since important for CO chemistry in molecular clouds */ 00082 if( dense.IonHigh[ipHELIUM] == 1 ) 00083 { 00084 trimhi = MIN2( trimhi , 1e-20f ); 00085 } 00086 else if( dense.IonHigh[ipHELIUM] == 2 ) 00087 { 00088 if( conv.lgSearch ) 00089 { 00090 /* during search phase we may be quite far from solution, in the 00091 * case of a PDR sim the electron density may be vastly higher 00092 * than it will be with stable solution. He++ can be very important 00093 * for the chemistry and we want to consider it. Make the 00094 * threshold for ignoring He++ much higher than normal to prevent 00095 * a premature removal of the ion */ 00096 trimhi = MIN2( trimhi , 1e-17f ); 00097 } 00098 else 00099 { 00100 /* similar smaller upper limit for ion*/ 00101 trimhi = MIN2( trimhi , 1e-12f ); 00102 } 00103 } 00104 } 00105 00106 /* logic for PDRs, for elements included in chemistry, need stable solutions, 00107 * keep 3 ion stages in most cases - added by NA to do HII/PDR sims */ 00108 if( mole.lgElem_in_chemistry[nelem] ) 00109 { 00110 trimlo = SMALLFLOAT; 00111 if(dense.IonHigh[nelem] ==2) 00112 { 00113 trimhi = MIN2(trimhi, 1e-20f); 00114 } 00115 } 00116 00117 /* raise or lower highest and lowest stages of ionization to 00118 * consider only significant stages 00119 * IonHigh[nelem] stage of ionization */ 00120 00121 /* this is a special block for initialization only - it checks on absurd abundances 00122 * and can trim multiple stages of ionization at one time. */ 00123 if( conv.lgSearch ) 00124 { 00125 /* special - trim down higher stages if they have essentially zero abundance */ 00126 while( 00127 (dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] < dense.density_low_limit || 00128 dense.xIonDense[nelem][dense.IonHigh[nelem]] < dense.density_low_limit ) && 00129 /* >>chng 02 may 12, rm +1 since this had effect of not allowing fully atomic */ 00130 dense.IonHigh[nelem] > dense.IonLow[nelem] ) 00131 { 00132 /* dense.xIonDense[nelem][i] is density of i+1 th ionization stage (cm^-3) 00133 * the -1 is correct for the heating, -1 since heating is caused by ionization of stage below it */ 00134 dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.; 00135 thermal.heating[nelem][dense.IonHigh[nelem]-1] = 0.; 00136 if( dense.IonHigh[nelem] == nelem+1-NISO ) 00137 { 00138 long int ipISO = nelem - dense.IonHigh[nelem]; 00139 ASSERT( ipISO>=0 && ipISO<NISO ); 00140 StatesElem[ipISO][nelem][0].Pop = 0.; 00141 } 00142 00143 /* decrement high stage limit */ 00144 --dense.IonHigh[nelem]; 00145 ASSERT( dense.IonHigh[nelem] >= 0); 00146 /* remember this happened */ 00147 lgHi_Down = true; 00148 } 00149 00150 /* special - trim up lower stages trim if they have essentially zero abundance */ 00151 while( 00152 (dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] < dense.density_low_limit || 00153 dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit ) && 00154 dense.IonLow[nelem] < dense.IonHigh[nelem] - 1 ) 00155 { 00156 /* dense.xIonDense[nelem][i] is density of ith ionization stage (cm^-3) 00157 * there is no-1 since we are removing the agent that heats */ 00158 dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.; 00159 /* >>chng 01 aug 04, remove -1 which clobbers thermal.heating when IonLow == 0 */ 00160 /*thermal.heating[nelem][dense.IonLow[nelem]-1] = 0.;*/ 00161 thermal.heating[nelem][dense.IonLow[nelem]] = 0.; 00162 if( dense.IonLow[nelem] == nelem+1-NISO ) 00163 { 00164 long int ipISO = nelem - dense.IonLow[nelem]; 00165 ASSERT( ipISO>=0 && ipISO<NISO ); 00166 StatesElem[ipISO][nelem][0].Pop = 0.; 00167 } 00168 00169 /* increment low stage limit */ 00170 ++dense.IonLow[nelem]; 00171 lgLo_Up = true; 00172 } 00173 } 00174 00175 /* sanity check */ 00176 /* IonHigh can be equal to IonLow if both are zero - no ionization*/ 00177 ASSERT( dense.IonLow[nelem] < dense.IonHigh[nelem] || 00178 (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) || 00179 dense.lgSetIoniz[nelem] ); 00180 00181 /* this checks on error condition where the gas is stupendously highly ionized, 00182 * the low limit is already one less than high, and we need to raise the low 00183 * limit further */ 00184 if( dense.IonHigh[nelem] > 1 && 00185 (dense.IonLow[nelem]==dense.IonHigh[nelem]-1) && 00186 (dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit) ) 00187 { 00188 fprintf(ioQQQ, 00189 "PROBLEM DISASTER the density of ion %li of element %s is too small to be computed on this cpu.\n", 00190 dense.IonLow[nelem]+1, 00191 elementnames.chElementName[nelem]); 00192 fprintf(ioQQQ, 00193 "Turn off the element with the command ELEMENT XXX OFF or do not consider " 00194 "gas with low density, the current hydrogen density is %.2e cm-3.\n", 00195 dense.gas_phase[ipHYDROGEN]); 00196 fprintf(ioQQQ, 00197 "The outer radius of the cloud is %.2e cm - should a stopping " 00198 "radius be set?\n", 00199 radius.Radius ); 00200 fprintf(ioQQQ, 00201 "The abort flag is being set - the calculation is stopping.\n"); 00202 lgAbort = true; 00203 return; 00204 } 00205 00206 /* trim down high stages that have too small an abundance */ 00207 if( 00208 dense.IonHigh[nelem] > dense.IonLow[nelem] && 00209 ( (dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] <= 00210 trimhi ) || 00211 (dense.xIonDense[nelem][dense.IonHigh[nelem]] <= dense.density_low_limit ) ) 00212 ) 00213 { 00214 /* >>chng 03 sep 30, the atom and its first ion are a special case 00215 * since we want to compute even trivial ions in molecular clouds */ 00216 if( dense.IonHigh[nelem]>1 || 00217 (dense.IonHigh[nelem]==1&&dense.xIonDense[nelem][1]<100.*dense.density_low_limit) ) 00218 { 00219 dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.; 00220 thermal.heating[nelem][dense.IonHigh[nelem]-1] = 0.; 00221 if( dense.IonHigh[nelem] == nelem+1-NISO ) 00222 { 00223 long int ipISO = nelem - dense.IonHigh[nelem]; 00224 ASSERT( ipISO>=0 && ipISO<NISO ); 00225 StatesElem[ipISO][nelem][0].Pop = 0.; 00226 } 00227 --dense.IonHigh[nelem]; 00228 lgHi_Down = true; 00229 } 00230 } 00231 00232 /* trim up highest stages - will this be done? */ 00233 lgHi_Up_enable = true; 00234 /* trimming can be turned off */ 00235 if( !ionbal.lgTrimhiOn ) 00236 lgHi_Up_enable = false; 00237 /* high stage is already fully stripped */ 00238 if( dense.IonHigh[nelem]>=nelem+1 ) 00239 lgHi_Up_enable = false; 00240 /* we have previously trimmed this stage down */ 00241 if( lgHi_Down ) 00242 lgHi_Up_enable = false; 00243 /* we are in neutral gas */ 00244 if( dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN]<0.9 ) 00245 lgHi_Up_enable = false; 00246 00247 if( lgHi_Up_enable ) 00248 { 00249 realnum abundold = 0; 00250 if( nzone>2 ) 00251 abundold = struc.xIonDense[nelem][dense.IonHigh[nelem]][nzone-3]/ 00252 SDIV( struc.gas_phase[nelem][nzone-3]); 00253 realnum abundnew = dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem]; 00254 /* only raise highest stage if ionization potential of next highest stage is within 00255 * continuum array and the abundance of the highest stage is significant */ 00256 if( (Heavy.Valence_IP_Ryd[nelem][dense.IonHigh[nelem]] < rfield.anu[rfield.nflux-1] || 00257 Heavy.Valence_IP_Ryd[nelem][dense.IonHigh[nelem]] < phycon.te_ryd*10.) && 00258 dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] > 1e-4f && 00259 /* this checks that abundance of highest stage is increasing */ 00260 abundnew > abundold*1.01 ) 00261 { 00262 /*fprintf(ioQQQ,"uuppp %li %li \n", nelem, dense.IonHigh[nelem] );*/ 00263 /* raise highest level of ionization */ 00264 ++dense.IonHigh[nelem]; 00265 lgHi_Up = true; 00266 /* set abundance to almost zero so that sanity check elsewhere will be ok */ 00267 dense.xIonDense[nelem][dense.IonHigh[nelem]] = (realnum)(dense.density_low_limit*10.); 00268 } 00269 } 00270 00271 /* sanity check 00272 * IonHigh can be equal to IonLow if both are zero - no ionization*/ 00273 ASSERT( dense.IonLow[nelem] < dense.IonHigh[nelem] || 00274 (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) || 00275 dense.lgSetIoniz[nelem] ); 00276 00277 /* lower lowest stage of ionization if we have significant abundance at current lowest */ 00278 if( dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] > 1e-3f && 00279 dense.IonLow[nelem] > 0 ) 00280 { 00281 /* lower lowest level of ionization */ 00282 --dense.IonLow[nelem]; 00283 lgLo_Down = true; 00284 /* >>chng 01 aug 02, set this to zero so that sanity check elsewhere will be ok */ 00285 dense.xIonDense[nelem][dense.IonLow[nelem]] = (realnum)dense.density_low_limit; 00286 } 00287 00288 /* raise lowest stage of ionization, but only if we are near illuminated face of cloud*/ 00289 else if( nzone < 10 && 00290 (dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] <= (realnum)trimlo) && 00291 (dense.IonLow[nelem] < (dense.IonHigh[nelem] - 2) ) ) 00292 { 00293 /* raise lowest level of ionization */ 00294 dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.; 00295 /* no minus one in below since this is low bound, already bounds at atom */ 00296 thermal.heating[nelem][dense.IonLow[nelem]] = 0.; 00297 if( dense.IonLow[nelem] == nelem+1-NISO ) 00298 { 00299 long int ipISO = nelem - dense.IonLow[nelem]; 00300 ASSERT( ipISO>=0 && ipISO<NISO ); 00301 StatesElem[ipISO][nelem][0].Pop = 0.; 00302 } 00303 ++dense.IonLow[nelem]; 00304 lgLo_Up = true; 00305 } 00306 /* test on zero */ 00307 else if( (dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit) && 00308 /*>>chng 06 may 24, from IonLow < IonHigh to <IonHigh-1 because 00309 * old test would allow low to be raised too high, which then 00310 * leads to insanity */ 00311 (dense.IonLow[nelem] < dense.IonHigh[nelem]-1) ) 00312 { 00313 while(dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit && 00314 /* >>chng 07 feb 27 from < IonHigh to < IonHigh-1 to prevent raising 00315 * low to high */ 00316 (dense.IonLow[nelem] < dense.IonHigh[nelem]-1) ) 00317 { 00318 /* raise lowest level of ionization */ 00319 dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.; 00320 /* no minus one in below since this is low bound, already bounds at atom */ 00321 thermal.heating[nelem][dense.IonLow[nelem]] = 0.; 00322 if( dense.IonLow[nelem] == nelem+1-NISO ) 00323 { 00324 long int ipISO = nelem - dense.IonLow[nelem]; 00325 ASSERT( ipISO>=0 && ipISO<NISO ); 00326 StatesElem[ipISO][nelem][0].Pop = 0.; 00327 } 00328 ++dense.IonLow[nelem]; 00329 lgLo_Up = true; 00330 } 00331 } 00332 00333 /* sanity check */ 00334 /* IonHigh can be equal to IonLow if both are zero - no ionization*/ 00335 ASSERT( dense.IonLow[nelem] < dense.IonHigh[nelem] || 00336 (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) || 00337 dense.lgSetIoniz[nelem] ); 00338 00339 /* these are standard bounds checks that appear throughout this routine 00340 * dense.xIonDense[IonLow] is first one with positive abundances 00341 * 00342 * in case where lower ionization stage was just lowered the abundance 00343 * was set to dense.density_low_limit so test must be < dense.density_low_limit */ 00344 ASSERT( dense.IonLow[nelem] <= 1 || 00345 dense.xIonDense[nelem][dense.IonLow[nelem]-1] == 0. ); 00346 00347 ASSERT( (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) || lgLo_Up || 00348 dense.xIonDense[nelem][dense.IonLow[nelem]] >= dense.density_low_limit || 00349 dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] >= dense.density_low_limit || 00350 /*>>chng 06 may 24, include this to cover case where cap is set by not allowing 00351 * lower limit to come up to upper limit */ 00352 (dense.IonLow[nelem]==dense.IonHigh[nelem]-1 )); 00353 00354 { 00355 /* option to print out what has happened so far .... */ 00356 enum {DEBUG_LOC=false}; 00357 if( DEBUG_LOC && nelem == ipHELIUM ) 00358 { 00359 if( lgHi_Down ||lgHi_Up ||lgLo_Up ||lgLo_Down ) 00360 { 00361 fprintf(ioQQQ,"DEBUG TrimZone\t%li\t",nzone ); 00362 if( lgHi_Down ) 00363 { 00364 fprintf(ioQQQ,"high dn %li to %li", 00365 ion_hi_save , 00366 dense.IonHigh[nelem] ); 00367 } 00368 if( lgHi_Up ) 00369 { 00370 fprintf(ioQQQ,"high up %li to %li", 00371 ion_hi_save , 00372 dense.IonHigh[nelem] ); 00373 } 00374 if( lgLo_Up ) 00375 { 00376 fprintf(ioQQQ,"low up %li to %li", 00377 ion_lo_save , 00378 dense.IonLow[nelem] ); 00379 } 00380 if( lgLo_Down ) 00381 { 00382 fprintf(ioQQQ,"low dn %li to %li", 00383 ion_lo_save , 00384 dense.IonLow[nelem] ); 00385 } 00386 fprintf(ioQQQ,"\n" ); 00387 } 00388 } 00389 } 00390 00391 /* only assert that the stage above the highest has a population of 00392 * zero if that stage exists */ 00393 if(dense.IonHigh[nelem] < nelem+1 ) 00394 ASSERT( dense.xIonDense[nelem][dense.IonHigh[nelem]+1] == 0. ); 00395 00396 /* option to print if any trimming occurred */ 00397 if( lgHi_Down || lgHi_Up || lgLo_Up || lgLo_Down ) 00398 { 00399 conv.lgIonStageTrimed = true; 00400 { 00401 /* option to print out what has happened so far .... */ 00402 enum {DEBUG_LOC=false}; 00403 if( DEBUG_LOC && nelem==ipHELIUM ) 00404 { 00405 fprintf(ioQQQ,"DEBUG ion_trim zone\t%.2f \t%li\t", fnzone, nelem); 00406 if( lgHi_Down ) 00407 fprintf(ioQQQ,"\tHi_Down"); 00408 if( lgHi_Up ) 00409 fprintf(ioQQQ,"\tHi_Up"); 00410 if( lgLo_Up ) 00411 fprintf(ioQQQ,"\tLo_Up"); 00412 if( lgLo_Down ) 00413 fprintf(ioQQQ,"\tLo_Down"); 00414 fprintf(ioQQQ,"\n"); 00415 } 00416 } 00417 } 00418 00419 /* asserts that that densities of ion stages that are now turned 00420 * off are zero */ 00421 for( ion=0; ion<dense.IonLow[nelem]; ++ion ) 00422 { 00423 ASSERT( dense.xIonDense[nelem][ion] == 0. ); 00424 } 00425 for( ion=dense.IonHigh[nelem]+1; ion<nelem+1; ++ion ) 00426 { 00427 ASSERT( dense.xIonDense[nelem][ion] == 0. ); 00428 } 00429 00430 for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion ) 00431 { 00432 /* in case where ionization stage was changed the 00433 * density is zero, we set it to SMALLFLOAT since a density 00434 * of zero is not possible */ 00435 dense.xIonDense[nelem][ion] = MAX2(dense.xIonDense[nelem][ion], SMALLFLOAT); 00436 } 00437 return; 00438 }