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_solver solve the bi-diagonal matrix for ionization balance */ 00004 #include "cddefines.h" 00005 #include "yield.h" 00006 #include "prt.h" 00007 #include "continuum.h" 00008 #include "iso.h" 00009 #include "dynamics.h" 00010 #include "grainvar.h" 00011 #include "hmi.h" 00012 #include "mole.h" 00013 #include "thermal.h" 00014 #include "thirdparty.h" 00015 #include "conv.h" 00016 #include "secondaries.h" 00017 #include "phycon.h" 00018 #include "atmdat.h" 00019 #include "heavy.h" 00020 #include "elementnames.h" 00021 #include "dense.h" 00022 #include "radius.h" 00023 #include "ionbal.h" 00024 00025 void solveions(double *ion, double *rec, double *snk, double *src, 00026 long int nlev, long int nmax); 00027 00028 void ion_solver( 00029 /* this is element on the c scale, H is 0 */ 00030 long int nelem, 00031 /* option to print this element when called */ 00032 bool lgPrintIt) 00033 { 00034 /* use tridiag or general matrix solver? */ 00035 bool lgTriDiag=false; 00036 long int ion, 00037 limit, 00038 IonProduced, 00039 nej, 00040 ns, 00041 jmax=-1; 00042 double totsrc; 00043 double dennew, rateone; 00044 bool lgNegPop; 00045 static bool lgMustMalloc=true; 00046 static double *sink, *source , *sourcesave; 00047 int32 nerror; 00048 static int32 *ipiv; 00049 long ion_low, ion_range, i, j, ion_to , ion_from; 00050 static double sum_dense; 00051 /* only used for debug printout */ 00052 static double *auger; 00053 00054 double abund_total, renormnew; 00055 bool lgHomogeneous = true; 00056 static double *xmat , *xmatsave; 00057 00058 DEBUG_ENTRY( "ion_solver()" ); 00059 00060 /* this is on the c scale, so H is 0 */ 00061 ASSERT( nelem >= 0); 00062 ASSERT( dense.IonLow[nelem] >= 0 ); 00063 ASSERT( dense.IonHigh[nelem] >= 0 ); 00064 00065 /* H is special because its abundance spills into three routines - 00066 * the ion/atom solver (this routine), the H-mole solvers (hmole), and 00067 * the heavy molecule solver. xmolecules only includes the heavy mole 00068 * part for H only. So the difference between gas_phase and xmolecules 00069 * includes the H2 part of the molecular network. This branch does 00070 * this special H case, then the general case (heavy elements are 00071 * not part of the H2 molecular network) */ 00072 00073 /* >>chng 01 dec 07, define abund_total, total atom and ion abundance 00074 * here, removing molecules */ 00075 if( nelem == ipHYDROGEN ) 00076 { 00077 /* Hydrogen is a special case since hmole does not include the 00078 * atom/molecules - hmole collapses its network into H = H0 + H+ 00079 * and then forces the solution determined here for partitioning 00080 * between these two */ 00081 abund_total = dense.xIonDense[nelem][0] + dense.xIonDense[nelem][1]; 00082 } 00083 else 00084 { 00085 abund_total = SDIV( dense.gas_phase[nelem] - dense.xMolecules[nelem] ); 00086 } 00087 00088 /* >>chng 04 nov 22, 00089 * if gas nearly all atomic/ionic do not let source - sink terms from 00090 * molecular network force system of balance equations to become 00091 * inhomogeneous 00092 * what constitutes a source or sink IS DIFFERENT for hydrogen and the rest 00093 * the H solution must couple with hmole - and its defn of source and sink. For instance, oxygen charge 00094 * transfer goes into source and sink terms for hydrogen. So we never hose source and sink for H. 00095 * for the heavy elements, which couple onto comole, mole.source and sink represent terms that 00096 * remove atoms and ions from the ionization ladder. their presence makes the system of 00097 * equations inhomogeneous. we don't want to do this when comole has a trivial effect on 00098 * the ionization balance since the matrix becomes unstable */ 00099 /* >>chng 04 dec 06, limit from 0.01 to 1e-10 as per NA suggestion */ 00100 if( nelem>ipHYDROGEN && dense.xMolecules[nelem]/SDIV(dense.gas_phase[nelem]) < 1e-10 ) 00101 { 00102 for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion ) 00103 { 00104 mole.source[nelem][ion] = 0.; 00105 mole.sink[nelem][ion] = 0.; 00106 } 00107 } 00108 00109 /* protect against case where all gas phase abundances are in molecules, use the 00110 * atomic and first ion density from the molecule solver 00111 * >>chng 04 aug 15, NA change from 10 0000 to 10 pre-coef on 00112 * FLT_EPSILON for stability in PDR 00113 * the factor 10.*FLT_EPSILON also appears in mole_h_step in total H 00114 * conservation */ 00115 if( fabs( dense.gas_phase[nelem] - dense.xMolecules[nelem])/SDIV(dense.gas_phase[nelem] ) < 00116 10.*FLT_EPSILON ) 00117 { 00118 double renorm; 00119 /* >>chng 04 jul 31, add logic to conserve nuclei in fully molecular limit; 00120 * in first calls, when searching for solution, we may be VERY far 00121 * off, and sum of first ion and atom density may be far larger 00122 * than difference between total gas and molecular densities, 00123 * since they reflect the previous evaluation of the solution. Do 00124 * renorm to cover this case */ 00125 /* first form sum of all atoms and ions */ 00126 realnum sum = 0.; 00127 for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion ) 00128 sum += dense.xIonDense[nelem][ion]; 00129 /* now renorm to this sum - this should be unity, and is not if we have 00130 * now conserved particles, due to molecular fraction changing */ 00131 renorm = dense.gas_phase[nelem] / SDIV(sum + dense.xMolecules[nelem] ); 00132 00133 abund_total = renorm * sum; 00134 } 00135 00136 /* negative residual density */ 00137 if( abund_total < 0. ) 00138 { 00139 /* >>chng 05 dec 16, do not abort when net atomic/ ionic abundance 00140 * is negative due to chem net having too much of a species - this 00141 * happens naturally when ices become well formed, but the code will 00142 * converge away from it. Rather, we set the ionization 00143 * convergence flag to "not converge" and soldier on 00144 * if negative populations do not go away, we will eventually 00145 * terminate due to convergence failures */ 00146 /* print error if not search */ 00147 if(!conv.lgSearch ) 00148 fprintf(ioQQQ, 00149 "PROBLEM ion_solver - neg net atomic abundance zero for nelem= %li, rel val= %.2e conv.nTotalIoniz=%li, fixed\n", 00150 nelem, 00151 fabs(abund_total) / SDIV(dense.xMolecules[nelem]), 00152 conv.nTotalIoniz ); 00153 /* fix up is to use half the positive abundance, assuming chem is 00154 * trying to get too much of this species */ 00155 abund_total = -abund_total/2.; 00156 00157 /* say that ionization is not converged - do not abort - but if 00158 * cannot converge away from negative solution, this will become a 00159 * convergence failure abort */ 00160 conv.lgConvIoniz = false; 00161 strcpy( conv.chConvIoniz, "neg ion" ); 00162 } 00163 00164 /* return if IonHigh is zero, since no ionization at all */ 00165 if( dense.IonHigh[nelem] == 0 ) 00166 { 00167 /* set the atom to the total gas phase abundance */ 00168 dense.xIonDense[nelem][0] = (realnum)abund_total; 00169 return; 00170 } 00171 00172 /* >>chng 01 may 09, add option to force ionization distribution with 00173 * element name ioniz */ 00174 if( dense.lgSetIoniz[nelem] ) 00175 { 00176 for( ion=0; ion<nelem+2; ++ion ) 00177 { 00178 dense.xIonDense[nelem][ion] = dense.SetIoniz[nelem][ion]* 00179 (realnum)abund_total; 00180 } 00181 return; 00182 } 00183 00184 /* impossible for HIonFrac[nelem] to be zero if IonHigh(nelem)=nelem+1 00185 * HIonFrac(nelem) is stripped to hydrogen */ 00186 /* >>chng 01 oct 30, to assert */ 00187 ASSERT( (dense.IonHigh[nelem] < nelem + 1) || 00188 iso.pop_ion_ov_neut[ipH_LIKE][nelem] > 0. ); 00189 00190 /* zero out the ionization and recombination rates that we will modify here, 00191 * but not the iso-electronic stages which are done elsewhere, 00192 * the nelem stage of ionization is he-like, 00193 * the nelem+1 stage of ionization is h-like */ 00194 00195 /* loop over stages of ionization that we solve for here, 00196 * up through and including one less than nelem-NISO, 00197 * never actually do highest NISO stages of ionization since they 00198 * come from the ionization ratio from the next lower stage */ 00199 limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1); 00200 00201 /* the full range of ionization - this is number of ionization stages */ 00202 ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1; 00203 ASSERT( ion_range <= nelem+2 ); 00204 00205 ion_low = dense.IonLow[nelem]; 00206 00207 /* PARALLEL_MODE true if there can be several instances of this routine 00208 * running in the same memory pool - we must have fresh instances of the 00209 * arrays for each */ 00210 if( PARALLEL_MODE || lgMustMalloc ) 00211 { 00212 long int ncell=LIMELM+1; 00213 lgMustMalloc = false; 00214 if( PARALLEL_MODE ) 00215 ncell = ion_range; 00216 00217 /* this will be "new" matrix, with non-adjacent coupling included */ 00218 xmat=(double*)MALLOC( (sizeof(double)*(unsigned)(ncell*ncell) )); 00219 xmatsave=(double*)MALLOC( (sizeof(double)*(unsigned)(ncell*ncell) )); 00220 sink=(double*)MALLOC( (sizeof(double)*(unsigned)ncell )); 00221 source=(double*)MALLOC( (sizeof(double)*(unsigned)ncell )); 00222 sourcesave=(double*)MALLOC( (sizeof(double)*(unsigned)ncell )); 00223 ipiv=(int32*)MALLOC( (sizeof(int32)*(unsigned)ncell )); 00224 /* auger is used only for debug printout - it is special because with multi-electron 00225 * Auger ejection, very high stages of ionization can be produced, well beyond 00226 * the highest stage considered here. so we malloc to the limit */ 00227 auger=(double*)MALLOC( (sizeof(double)*(unsigned)(LIMELM+1) )); 00228 } 00229 00230 for( i=0; i<ion_range;i++ ) 00231 { 00232 source[i] = 0.; 00233 } 00234 00235 /* this will be used to address the 2d arrays */ 00236 # undef MAT 00237 # define MAT(M_,I_,J_) (*((M_)+(I_)*(ion_range)+(J_))) 00238 00239 /* zero-out loop comes before main loop since there are off-diagonal 00240 * elements in the main ionization loop, due to multi-electron processes, 00241 * TotIonizRate and TotRecom were already set in h-like and he-like solvers 00242 * other recombination rates were already set by routines responsible for them */ 00243 for( ion=0; ion <= limit; ion++ ) 00244 { 00245 ionbal.RateIonizTot[nelem][ion] = 0.; 00246 } 00247 00248 /* auger is used only for debug printout - it is special because with multi-electron 00249 * Auger ejection, very high stages of ionization can be produced, well beyond 00250 * the highest stage considered here. so we malloc to the limit */ 00251 for( i=0; i< LIMELM+1; ++i ) 00252 { 00253 auger[i] = 0.; 00254 } 00255 00256 /* zero out xmat */ 00257 for( i=0; i< ion_range; ++i ) 00258 { 00259 for( j=0; j< ion_range; ++j ) 00260 { 00261 MAT( xmat, i, j ) = 0.; 00262 } 00263 } 00264 00265 { 00266 /* this sets up a fake ionization balance problem, with a trivial solution, 00267 * for debugging the ionization solver */ 00268 enum {DEBUG_LOC=false}; 00269 if( DEBUG_LOC && nelem==ipCARBON ) 00270 { 00271 broken();/* within optional debug print statement */ 00272 dense.IonLow[nelem] = 0; 00273 dense.IonHigh[nelem] = 3; 00274 abund_total = 1.; 00275 ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1; 00276 /* make up ionization and recombination rates */ 00277 for( ion=dense.IonLow[nelem]; ion <= limit; ion++ ) 00278 { 00279 double fac=1; 00280 if(ion) 00281 fac = 1e-10; 00282 ionbal.RateRecomTot[nelem][ion] = 100.; 00283 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ ) 00284 { 00285 /* direct photoionization of this shell */ 00286 ionbal.PhotoRate_Shell[nelem][ion][ns][0] = fac; 00287 } 00288 } 00289 } 00290 } 00291 00292 /* Now put in all recombination and ionization terms from CO_mole() that 00293 * come from molecular reactions. this traces molecular process that 00294 * change ionization stages with this ladder - but do not remove from 00295 * the ladder */ 00296 for( ion_to=dense.IonLow[nelem]; ion_to <= limit; ion_to++ ) 00297 { 00298 for( ion_from=dense.IonLow[nelem]; ion_from <= dense.IonHigh[nelem]; ++ion_from ) 00299 { 00300 /* do not do ion onto itself */ 00301 if( ion_to != ion_from ) 00302 { 00303 /* this is the rate coefficient for charge transfer from ion to ion_to */ 00304 /*rateone = gv.GrainChTrRate[nelem][ion_from][ion_to];*/ 00305 rateone = mole.xMoleChTrRate[nelem][ion_from][ion_to]; 00306 ASSERT( rateone >= 0. ); 00307 MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone; 00308 MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone; 00309 } 00310 } 00311 } 00312 00313 /* now get actual arrays of ionization and recombination processes, 00314 * but only for the ions that are done as two-level systems */ 00315 /* in two-stage system, atom + first ion, limit is zero but must 00316 * include gv.GrainChTrRate[nelem][1][0] */ 00317 /* grain charge transfer */ 00318 if( gv.lgDustOn && ionbal.lgGrainIonRecom && gv.lgGrainPhysicsOn ) 00319 { 00320 long int low; 00321 /* do not double count this process for atoms that are in the co network - we use 00322 * a net recombination coefficient derived from the co solution, this includes grain ct */ 00323 /* >>chng 05 dec 23, add mole.lgElem_in_chemistry */ 00324 if( mole.lgElem_in_chemistry[nelem] ) 00325 /*if( nelem==ipHYDROGEN ||nelem==ipCARBON ||nelem== ipOXYGEN ||nelem==ipSILICON ||nelem==ipSULPHUR 00326 ||nelem==ipNITROGEN ||nelem==ipCHLORINE )*/ 00327 { 00328 low = MAX2(1, dense.IonLow[nelem] ); 00329 } 00330 else 00331 low = dense.IonLow[nelem]; 00332 00333 for( ion_to=low; ion_to <= dense.IonHigh[nelem]; ion_to++ ) 00334 { 00335 for( ion_from=dense.IonLow[nelem]; ion_from <= dense.IonHigh[nelem]; ++ion_from ) 00336 { 00337 /* do not do ion onto itself */ 00338 if( ion_to != ion_from ) 00339 { 00340 /* this is the rate coefficient for charge transfer from ion to ion_to 00341 * both grain charge transfer ionization and recombination */ 00342 rateone = gv.GrainChTrRate[nelem][ion_from][ion_to]; 00343 MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone; 00344 MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone; 00345 } 00346 } 00347 } 00348 } 00349 00350 for( ion=dense.IonLow[nelem]; ion <= limit; ion++ ) 00351 { 00352 /* thermal & secondary collisional ionization */ 00353 rateone = ionbal.CollIonRate_Ground[nelem][ion][0] + 00354 secondaries.csupra[nelem][ion] + 00355 /* inner shell ionization by UTA lines */ 00356 ionbal.UTA_ionize_rate[nelem][ion]; 00357 ionbal.RateIonizTot[nelem][ion] += rateone; 00358 00359 /* UTA ionization */ 00360 if( ion+1-ion_low < ion_range ) 00361 { 00362 /* depopulation processes enter with negative sign */ 00363 MAT( xmat, ion-ion_low, ion-ion_low ) -= rateone; 00364 MAT( xmat, ion-ion_low, ion+1-ion_low ) += rateone; 00365 } 00366 00367 /* total recombination rate */ 00368 if( ion-1-ion_low >= 0 ) 00369 { 00370 /* loss of this ion due to recombination to next lower ion stage */ 00371 MAT( xmat,ion-ion_low, ion-ion_low ) -= ionbal.RateRecomTot[nelem][ion-1]; 00372 MAT( xmat,ion-ion_low, ion-1-ion_low ) += ionbal.RateRecomTot[nelem][ion-1]; 00373 } 00374 00375 /* loop over all atomic sub-shells to include photoionization */ 00376 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ ) 00377 { 00378 /* direct photoionization of this shell */ 00379 ionbal.RateIonizTot[nelem][ion] += ionbal.PhotoRate_Shell[nelem][ion][ns][0]; 00380 00381 /* this is the primary ionization rate - add to diagonal element, 00382 * test on ion stage is so that we don't include ionization from the very highest 00383 * ionization stage to even higher - since those even higher stages are not considered 00384 * this would appear as a sink - but populations of this highest level is ensured to 00385 * be nearly trivial and neglecting it production of even higher ionization OK */ 00386 /* >>chng 04 nov 29 RJRW, include following in this branch so only 00387 * evaluated when below ions done with iso-sequence */ 00388 if( ion+1-ion_low < ion_range ) 00389 { 00390 /* this will be redistributed into charge states in following loop */ 00391 MAT( xmat, ion-ion_low, ion-ion_low ) -= ionbal.PhotoRate_Shell[nelem][ion][ns][0]; 00392 00393 /* t_yield::Inst().nelec_eject(nelem,ion,ns) is total number of electrons that can 00394 * possibly be freed 00395 * loop over nej, the number of electrons ejected including the primary, 00396 * nej = 1 is primary, nej > 1 includes primary plus Auger 00397 * t_yield::Inst().elec_eject_frac is prob of nej electrons */ 00398 for( nej=1; nej <= t_yield::Inst().nelec_eject(nelem,ion,ns); nej++ ) 00399 { 00400 /* this is the ion that is produced by this ejection, 00401 * limited by highest possible stage of ionization - 00402 * do not want to ignore ionization that go beyond this */ 00403 IonProduced = MIN2(ion+nej,dense.IonHigh[nelem]); 00404 rateone = ionbal.PhotoRate_Shell[nelem][ion][ns][0]* 00405 t_yield::Inst().elec_eject_frac(nelem,ion,ns,nej-1); 00406 /* >>chng 04 sep 06, above had included factor of nej to get rate, but 00407 * actually want events into particular ion */ 00408 /* number of electrons ejected 00409 *(double)nej; */ 00410 /* it goes into this charge state - recall upper cap due to ion stage trimming 00411 * note that compensating loss term on diagonal was done before this 00412 * loop, since frac_elec_eject adds to unity */ 00413 MAT( xmat, ion-ion_low, IonProduced-ion_low ) += rateone; 00414 00415 /* only used for possible printout - multiple electron Auger rate - 00416 * do not count one-electron as Auger */ 00417 if( nej>1 ) 00418 auger[IonProduced-1] += rateone; 00419 } 00420 } 00421 } 00422 00423 /* this is charge transfer ionization of this species by hydrogen and helium */ 00424 rateone = 00425 atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[ipHELIUM][1]+ 00426 atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[ipHYDROGEN][1]; 00427 ionbal.RateIonizTot[nelem][ion] += rateone; 00428 if( ion+1-ion_low < ion_range ) 00429 { 00430 MAT( xmat, ion-ion_low, ion-ion_low ) -= rateone; 00431 MAT( xmat, ion-ion_low, ion+1-ion_low ) += rateone; 00432 } 00433 } 00434 00435 /* after this loop, ionbal.RateIonizTot and ionbal.RateRecomTot have been defined for the 00436 * stages of ionization that are done with simple */ 00437 /* begin loop at first species that is treated with full model atom */ 00438 j = MAX2(0,limit+1); 00439 /* possible that lowest stage of ionization is higher than this */ 00440 j = MAX2( ion_low , j ); 00441 for( ion=j; ion<=dense.IonHigh[nelem]; ion++ ) 00442 { 00443 ASSERT( ion>=0 && ion<nelem+2 ); 00444 /* use total ionization/recombination rates for species done with ISO solver */ 00445 if( ion+1-ion_low < ion_range ) 00446 { 00447 /* depopulation processes enter with negative sign */ 00448 MAT( xmat, ion-ion_low, ion-ion_low ) -= ionbal.RateIonizTot[nelem][ion]; 00449 MAT( xmat, ion-ion_low, ion+1-ion_low ) += ionbal.RateIonizTot[nelem][ion]; 00450 } 00451 00452 if( ion-1-ion_low >= 0 ) 00453 { 00454 /* loss of this ion due to recombination to next lower ion stage */ 00455 MAT( xmat,ion-ion_low, ion-ion_low ) -= ionbal.RateRecomTot[nelem][ion-1]; 00456 MAT( xmat,ion-ion_low, ion-1-ion_low ) += ionbal.RateRecomTot[nelem][ion-1]; 00457 } 00458 } 00459 00460 /* this will be sum of source and sink terms, will be used to decide if 00461 * matrix is singular */ 00462 totsrc = 0.; 00463 /*>>chng 06 nov 28 only include source from molecules if we have an estimated first 00464 * solution - first test is that we have called mole at least twice, 00465 * second test is that we are on a later iteration */ 00466 if( conv.nTotalIoniz > 1 || iteration > 1 ) 00467 { 00468 00469 for( i=0; i<ion_range;i++ ) 00470 { 00471 ion = i+ion_low; 00472 00473 /* these are the external source and sink terms */ 00474 /* source first */ 00475 /* need negative sign to get positive pops */ 00476 source[i] -= mole.source[nelem][ion]; 00477 00478 totsrc += mole.source[nelem][ion]; 00479 /* sink next */ 00480 /*MAT( xmat, i, i ) += mole.sink[nelem][ion];*/ 00481 MAT( xmat, i, i ) -= mole.sink[nelem][ion]; 00482 } 00483 } 00484 00485 /* matrix is not homogeneous if source is non-zero */ 00486 if( totsrc != 0. ) 00487 lgHomogeneous = false; 00488 00489 /* chng 03 jan 13 rjrw, add in dynamics if required here, 00490 * last test - only do advection if we have not overrun the radius scale */ 00491 if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection 00492 /*&& radius.depth < dynamics.oldFullDepth*/ ) 00493 { 00494 for( i=0; i<ion_range;i++ ) 00495 { 00496 ion = i+ion_low; 00497 /*MAT( xmat, i, i ) += dynamics.Rate;*/ 00498 MAT( xmat, i, i ) -= dynamics.Rate; 00499 /*src[i] += dynamics.Source[nelem][ion];*/ 00500 source[i] -= dynamics.Source[nelem][ion]; 00501 /* fprintf(ioQQQ," %li %li %.3e (%.3e %.3e)\n",i,i,MAT(xmat,i,i), 00502 dynamics.Rate, dynamics.Source[nelem][ion]);*/ 00503 } 00504 lgHomogeneous = false; 00505 } 00506 00507 for( i=0; i< ion_range; ++i ) 00508 { 00509 for( j=0; j< ion_range; ++j ) 00510 { 00511 MAT( xmatsave, i, j ) = MAT( xmat, i, j ); 00512 } 00513 sourcesave[i] = source[i]; 00514 } 00515 00516 /* >>chng 06 nov 21, for very high ionization parameter sims the H molecular 00517 * fraction can become so small that atom = atom + molecule. In this case we 00518 * will not count system as an inhomogeneous case since linear algebra package 00519 * will fail. If molecules are very small, count as homogeneous. 00520 * Note that we have already added sink terms to the main matrix and they 00521 * will not be removed, a possible source of error, but they must not 00522 * have been significant, given that the molecular fraction is so small */ 00523 if( !lgHomogeneous && ion_range==2 ) 00524 { 00525 /* solve algebraically */ 00526 double a = MAT( xmatsave, 0, 0 ), 00527 b = MAT( xmatsave, 1, 0 ) , 00528 c = MAT( xmatsave, 0, 1 ) , 00529 d = MAT( xmatsave, 1, 1 ); 00530 b = SDIV(b); 00531 d = SDIV(d); 00532 double ratio1 = a/b , ratio2 = c/d , fratio1=fabs(a/b),fratio2=fabs(c/d); 00533 if( fabs(ratio1-ratio2)/MAX2(fratio1,fratio2) <DBL_EPSILON*1e4 ) 00534 { 00535 abund_total = SDIV( dense.gas_phase[nelem] - dense.xMolecules[nelem] ); 00536 lgHomogeneous = true; 00537 } 00538 } 00539 # if 0 00540 if( nelem==ipHYDROGEN && 00541 fabs(dense.xMolecules[nelem]) / SDIV(dense.gas_phase[ipHYDROGEN]) <DBL_EPSILON*100. ) 00542 { 00543 abund_total = SDIV( dense.gas_phase[nelem] - dense.xMolecules[nelem] ); 00544 lgHomogeneous = true; 00545 } 00546 # endif 00547 00548 /* this is true if no source terms 00549 * we will use total population and species conservation to replace one 00550 * set of balance equations since overdetermined */ 00551 if( lgHomogeneous ) 00552 { 00553 double rate_ioniz=1., rate_recomb /*, scale = 0.*/; 00554 /* Simple estimate of most abundant ion */ 00555 jmax = 0; 00556 for( i=0; i<ion_range-1;i++) 00557 { 00558 ion = i+ion_low; 00559 rate_ioniz *= ionbal.RateIonizTot[nelem][ion]; 00560 rate_recomb = ionbal.RateRecomTot[nelem][ion]; 00561 /* trips if ion rate zero, so ll the gas will be more neutral than this */ 00562 if( rate_ioniz == 0) 00563 break; 00564 /* rec rate is zero */ 00565 if( rate_recomb <= 0.) 00566 break; 00567 00568 rate_ioniz /= rate_recomb; 00569 if( rate_ioniz > 1.) 00570 { 00571 /* this is peak ionization stage */ 00572 jmax = i; 00573 rate_ioniz = 1.; 00574 } 00575 } 00576 00577 /* replace its matrix elements with population sum */ 00578 for( i=0; i<ion_range;i++ ) 00579 { 00580 MAT(xmat,i,jmax) = 1.; 00581 } 00582 source[jmax] = abund_total; 00583 } 00584 00585 for( i=0; i< ion_range; ++i ) 00586 { 00587 for( j=0; j< ion_range; ++j ) 00588 { 00589 MAT( xmatsave, i, j ) = MAT( xmat, i, j ); 00590 } 00591 sourcesave[i] = source[i]; 00592 } 00593 00594 if( false && nelem == ipHYDROGEN && dynamics.lgAdvection&& iteration>1 ) 00595 { 00596 fprintf(ioQQQ,"DEBUGG Rate %.2f %.3e \n",fnzone,dynamics.Rate); 00597 fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateIonizTot[nelem][0], ionbal.RateIonizTot[nelem][1]); 00598 fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateRecomTot[nelem][0], ionbal.RateRecomTot[nelem][1]); 00599 fprintf(ioQQQ," %.3e %.3e %.3e\n\n", dynamics.Source[nelem][0], dynamics.Source[nelem][1], dynamics.Source[nelem][2]); 00600 } 00601 00602 { 00603 /* option to print matrix */ 00604 /*@-redef@*/ 00605 enum {DEBUG_LOC=false}; 00606 /*@+redef@*/ 00607 if( DEBUG_LOC && nzone==1 && lgPrintIt ) 00608 { 00609 fprintf( ioQQQ, 00610 " DEBUG ion_solver: nelem=%li ion_range=%li, limit=%li, nConv %li xmat follows\n", 00611 nelem , ion_range,limit , conv.nTotalIoniz ); 00612 if( lgHomogeneous ) 00613 fprintf(ioQQQ , "Homogeneous \n"); 00614 for( i=0; i<ion_range; ++i ) 00615 { 00616 for( j=0;j<ion_range;j++ ) 00617 { 00618 fprintf(ioQQQ,"%e\t",MAT(xmatsave,i,j)); 00619 } 00620 fprintf(ioQQQ,"\n"); 00621 } 00622 fprintf(ioQQQ,"source follows\n"); 00623 for( i=0; i<ion_range;i++ ) 00624 { 00625 fprintf(ioQQQ,"%e\t",source[i]); 00626 } 00627 fprintf(ioQQQ,"\n"); 00628 } 00629 } 00630 00631 /* first branch is tridiag, following branch is just general matrix solver */ 00632 if( lgTriDiag ) 00633 { 00634 bool lgLapack=false , lgTriOptimized=true; 00635 /* there are two choices for the tridiag solver */ 00636 if(lgLapack) 00637 { 00638 /* this branch uses lapack tridiag solver, and should work 00639 * it is hardwired to not be used because not extensively tested 00640 * issues - is this slower than others, and is it robust in 00641 * singular cases? */ 00642 bool lgBad = false; 00643 /* Use Lapack tridiagonal solver */ 00644 double *dl, *d, *du; 00645 00646 d = (double *) MALLOC((unsigned)ion_range*sizeof(double)); 00647 du = (double *) MALLOC((unsigned)(ion_range-1)*sizeof(double)); 00648 dl = (double *) MALLOC((unsigned)(ion_range-1)*sizeof(double)); 00649 00650 for(i=0;i<ion_range-1;i++) 00651 { 00652 du[i] = MAT(xmat,i+1,i); 00653 d[i] = MAT(xmat,i,i); 00654 dl[i] = MAT(xmat,i,i+1); 00655 } 00656 d[i] = MAT(xmat,i,i); 00657 00658 # if 0 00659 int i1, i2; 00660 i1 = ion_range; 00661 i2 = 1; 00662 lgBad = dgtsv_wrapper(&i1, &i2, dl, d, du, source, &i2); 00663 # endif 00664 if( lgBad ) 00665 fprintf(ioQQQ," dgtsz error.\n"); 00666 00667 free(dl);free(du);free(d); 00668 } 00669 else if(lgTriOptimized) 00670 { 00671 /* Use tridiagonal solver re-coded to avoid rounding errors 00672 * on diagonal -- uses determination of jmax for the 00673 * singular case, but is otherwise independent of the array 00674 * filling code above */ 00675 00676 for(i=0;i<ion_range;i++) 00677 { 00678 source[i] = sink[i] = 0.; 00679 } 00680 if(nelem == ipHYDROGEN && !hmi.lgNoH2Mole) 00681 { 00682 for(i=0;i<ion_range;i++) 00683 { 00684 ion = i+ion_low; 00685 source[i] += mole.source[nelem][ion]; 00686 sink[i] += mole.sink[nelem][ion]; 00687 } 00688 } 00689 if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection && radius.depth < dynamics.oldFullDepth ) 00690 { 00691 for(i=0;i<ion_range;i++) 00692 { 00693 ion = i+ion_low; 00694 source[i] += dynamics.Source[nelem][ion]; 00695 sink[i] += dynamics.Rate; 00696 } 00697 } 00698 00699 solveions(ionbal.RateIonizTot[nelem]+ion_low,ionbal.RateRecomTot[nelem]+ion_low, 00700 sink,source,ion_range,jmax); 00701 } 00702 } 00703 else 00704 { 00705 /* this is the default solver - now get new solution */ 00706 nerror = 0; 00707 /* Use general matrix solver */ 00708 getrf_wrapper(ion_range, ion_range, xmat, ion_range, ipiv, &nerror); 00709 if( nerror != 0 ) 00710 { 00711 fprintf( ioQQQ, 00712 " DISASTER ion_solver: dgetrf finds singular or ill-conditioned matrix nelem=%li %s ion_range=%li, limit=%li, nConv %li xmat follows\n", 00713 nelem , 00714 elementnames.chElementSym[nelem], 00715 ion_range, 00716 limit , 00717 conv.nTotalIoniz ); 00718 for( i=0; i<ion_range; ++i ) 00719 { 00720 for( j=0;j<ion_range;j++ ) 00721 { 00722 fprintf(ioQQQ,"%e\t",MAT(xmatsave,j,i)); 00723 } 00724 fprintf(ioQQQ,"\n"); 00725 } 00726 fprintf(ioQQQ,"source follows\n"); 00727 for( i=0; i<ion_range;i++ ) 00728 { 00729 fprintf(ioQQQ,"%e\t",source[i]); 00730 } 00731 fprintf(ioQQQ,"\n"); 00732 cdEXIT(EXIT_FAILURE); 00733 } 00734 getrs_wrapper('N', ion_range, 1, xmat, ion_range, ipiv, source, ion_range, &nerror); 00735 if( nerror != 0 ) 00736 { 00737 fprintf( ioQQQ, " DISASTER ion_solver: dgetrs finds singular or ill-conditioned matrix nelem=%li ionrange=%li\n", 00738 nelem , ion_range ); 00739 cdEXIT(EXIT_FAILURE); 00740 } 00741 } 00742 00743 { 00744 /*@-redef@*/ 00745 /* this is to debug following failed assert */ 00746 enum {DEBUG_LOC=false}; 00747 /*@+redef@*/ 00748 if( DEBUG_LOC && nelem == ipHYDROGEN ) 00749 { 00750 fprintf(ioQQQ,"debuggg ion_solver1 \t%.2f\t%.4e\t%.4e\tIon\t%.3e\tRec\t%.3e\n", 00751 fnzone, 00752 phycon.te, 00753 dense.eden, 00754 ionbal.RateIonizTot[nelem][0] , 00755 ionbal.RateRecomTot[nelem][0]); 00756 fprintf(ioQQQ," Msrc %.3e %.3e\n", mole.source[ipHYDROGEN][0], mole.source[ipHYDROGEN][1]); 00757 fprintf(ioQQQ," Msnk %.3e %.3e\n", mole.sink[ipHYDROGEN][0], mole.sink[ipHYDROGEN][1]); 00758 fprintf(ioQQQ," Poprat %.3e nomol %.3e\n",source[1]/source[0], 00759 ionbal.RateIonizTot[nelem][0]/ionbal.RateRecomTot[nelem][0]); 00760 } 00761 } 00762 00763 for( i=0; i<ion_range; i++ ) 00764 { 00765 /* is this blows, source[i] is NaN */ 00766 ASSERT( !isnan( source[i] ) ); 00767 } 00768 00769 # define RJRW 0 00770 if( RJRW && 0 ) 00771 { 00772 /* verify that the rates are sensible */ 00773 double test; 00774 for(i=0; i<ion_range; i++) { 00775 test = 0.; 00776 for(j=0; j<ion_range; j++) { 00777 test = test+source[j]*MAT(xmatsave,j,i); 00778 } 00779 fprintf(ioQQQ,"%e\t",test); 00780 } 00781 fprintf(ioQQQ,"\n"); 00782 00783 test = 0.; 00784 fprintf(ioQQQ," ion %li abundance %.3e\n",nelem,abund_total); 00785 for( ion=dense.IonLow[nelem]; ion < dense.IonHigh[nelem]; ion++ ) 00786 { 00787 if( ionbal.RateRecomTot[nelem][ion] != 0 && source[ion-ion_low] != 0 ) 00788 fprintf(ioQQQ," %li %.3e %.3e : %.3e\n",ion,source[ion-ion_low], 00789 source[ion-ion_low+1]/source[ion-ion_low], 00790 ionbal.RateIonizTot[nelem][ion]/ionbal.RateRecomTot[nelem][ion]); 00791 else 00792 fprintf(ioQQQ," %li %.3e [One ratio infinity]\n",ion,source[ion-ion_low]); 00793 test += source[ion-ion_low]; 00794 } 00795 test += source[ion-ion_low]; 00796 } 00797 00798 /* 00799 * >> chng 03 jan 15 rjrw:- terms are now included for 00800 * molecular sources and sinks of H and H+. 00801 * 00802 * When the network is not in equilibrium, this will lead to a 00803 * change in the derived abundance of H and H+ when after the 00804 * matrix solution -- the difference between `renorm' and 1. is a 00805 * measure of the quality of the solution (it will be 1. if the 00806 * rate of transfer into Ho/H+ balances the rate of transfer 00807 * out, for the consistent relative abundances). 00808 * 00809 * We therefore renormalize to keep the total H abundance 00810 * correct -- only the molecular network is allowed to change 00811 * this. 00812 * 00813 * To do this, only the ion abundances are corrected, as the 00814 * molecular abundances may depend on several different 00815 * conserved species. 00816 * 00817 */ 00818 if( lgHomogeneous ) 00819 { 00820 dense.xIonDense[nelem][ion_low] = (realnum)abund_total; 00821 for( i=1;i < ion_range; i++ ) 00822 { 00823 dense.xIonDense[nelem][i+ion_low] = 0.; 00824 } 00825 } 00826 00827 renormnew = 1.; 00828 /* >>chng 06 mar 17, comment out test on old full depth - keep old solution if overrun scale */ 00829 if(iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection /*&& radius.depth < dynamics.oldFullDepth */ && 00830 nelem == ipHYDROGEN && hmi.lgNoH2Mole) 00831 { 00832 /* The normalization out of the matrix solution is correct and 00833 * should be retained if: dynamics is on and the total 00834 * abundance of HI & H+ isn't being controlled by the 00835 * molecular network */ 00836 renormnew = 1.; 00837 } 00838 else 00839 { 00840 dennew = 0.; 00841 sum_dense = 0.; 00842 00843 /* find total population to renorm - also here check that negative pops did not occur */ 00844 for( i=0;i < ion_range; i++ ) 00845 { 00846 ion = i+ion_low; 00847 sum_dense += dense.xIonDense[nelem][ion]; 00848 dennew += source[i]; 00849 } 00850 00851 if( dennew > 0.) 00852 { 00853 renormnew = sum_dense / dennew; 00858 } 00859 else 00860 { 00861 renormnew = 1.; 00862 } 00863 } 00864 /* check not negative, should be +ve, can be zero if species has become totally molecular. 00865 * this happens for hydrogen if no cosmic rays, or cr ion set very low */ 00866 if( renormnew < 0) 00867 { 00868 fprintf(ioQQQ,"PROBLEM impossible value of renorm \n"); 00869 } 00870 ASSERT( renormnew>=0 ); 00871 00872 /* source[i] contains new solution for ionization populations 00873 * save resulting abundances into main ionization density array, 00874 * while checking whether any negative level populations occurred */ 00875 lgNegPop = false; 00876 for( i=0; i < ion_range; i++ ) 00877 { 00878 ion = i+ion_low; 00879 00880 /*fprintf(ioQQQ," %li %li %.3e %.3e\n",nelem,ion,src[ion-ion_low+1],src[ion-ion_low]); 00881 pop_ion_ov_neut[ion] = src[ion-ion_low+1]/src[ion-ion_low]; */ 00882 if( source[i] < 0. ) 00883 { 00884 /* >>chng 04 dec 04, put test on neg abund here, don't print uncles value is very -ve */ 00885 /* >>chng 06 feb 28, from -1e-10 to -1e-9, sim func_t10 had several negative 00886 * during initial search, due to extremely high ionization */ 00887 /* >>chng 06 mar 11, from 1e-9 to 2e-9 make many struc elements floats from double */ 00888 if( source[i]<-2e-9 ) 00889 fprintf(ioQQQ, 00890 " PROBLEM negative ion population in ion_solver, nelem=%li, %s ion=%li val=%.3e Search?%c zone=%li iteration=%li\n", 00891 nelem , 00892 elementnames.chElementSym[nelem], 00893 ion , 00894 source[i] , 00895 TorF(conv.lgSearch) , 00896 nzone , 00897 iteration ); 00898 source[i] = 0.; 00899 /* if this is one of the iso seq model atoms then must also zero out Lya and Pop2Ion */ 00900 if( ion == nelem+1-NISO ) 00901 { 00902 long int ipISO = nelem - ion; 00903 ASSERT( ipISO>=0 && ipISO<NISO ); 00904 StatesElem[ipISO][nelem][0].Pop = 0.; 00905 } 00906 } 00907 00908 /* use new solution */ 00909 dense.xIonDense[nelem][ion] = (realnum)(source[i]*renormnew); 00910 } 00911 00912 /* Zero levels with abundances < 1e-25 which which will suffer numerical noise */ 00913 while( dense.IonHigh[nelem] > dense.IonLow[nelem] && 00914 dense.xIonDense[nelem][dense.IonHigh[nelem]] < 1e-25*abund_total ) 00915 { 00916 ASSERT( dense.xIonDense[nelem][dense.IonHigh[nelem]] >= 0. ); 00917 /* zero out abundance and heating due to stage of ionization we are about to zero out */ 00918 dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.; 00919 thermal.heating[nelem][dense.IonHigh[nelem]-1] = 0.; 00920 /* decrement counter */ 00921 --dense.IonHigh[nelem]; 00922 } 00923 00924 /* sanity check, either offset stages of low and high ionization, 00925 * or no ionization at all */ 00926 ASSERT( (dense.IonLow[nelem] < dense.IonHigh[nelem]) || 00927 (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) ); 00928 00929 /* this should not happen */ 00930 if( lgNegPop ) 00931 { 00932 fprintf( ioQQQ, " PROBLEM Negative population found for abundance of ionization stage of element %4.4s, ZONE=%4ld\n", 00933 elementnames.chElementNameShort[nelem], nzone ); 00934 00935 fprintf( ioQQQ, " Populations were" ); 00936 for( ion=1; ion <= dense.IonHigh[nelem]+1; ion++ ) 00937 { 00938 fprintf( ioQQQ, "%9.1e", dense.xIonDense[nelem][ion-1] ); 00939 } 00940 fprintf( ioQQQ, "\n" ); 00941 00942 fprintf( ioQQQ, " destroy vector =" ); 00943 for( ion=1; ion <= dense.IonHigh[nelem]; ion++ ) 00944 { 00945 fprintf( ioQQQ, "%9.1e", ionbal.RateIonizTot[nelem][ion-1] ); 00946 } 00947 fprintf( ioQQQ, "\n" ); 00948 00949 /* print some extra stuff if destroy was negative */ 00950 if( lgNegPop ) 00951 { 00952 fprintf( ioQQQ, " CTHeavy vector =" ); 00953 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 00954 { 00955 fprintf( ioQQQ, "%9.1e", atmdat.HeCharExcIonOf[nelem][ion] ); 00956 } 00957 fprintf( ioQQQ, "\n" ); 00958 00959 fprintf( ioQQQ, " HCharExcIonOf vtr=" ); 00960 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 00961 { 00962 fprintf( ioQQQ, "%9.1e", atmdat.HCharExcIonOf[nelem][ion] ); 00963 } 00964 fprintf( ioQQQ, "\n" ); 00965 00966 fprintf( ioQQQ, " CollidRate vtr=" ); 00967 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 00968 { 00969 fprintf( ioQQQ, "%9.1e", ionbal.CollIonRate_Ground[nelem][ion][0] ); 00970 } 00971 fprintf( ioQQQ, "\n" ); 00972 00973 /* photo rates per subshell */ 00974 fprintf( ioQQQ, " photo rates per subshell, ion\n" ); 00975 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 00976 { 00977 fprintf( ioQQQ, "%3ld", ion ); 00978 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ ) 00979 { 00980 fprintf( ioQQQ, "%9.1e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] ); 00981 } 00982 fprintf( ioQQQ, "\n" ); 00983 } 00984 } 00985 00986 /* now check out creation vector */ 00987 fprintf( ioQQQ, " create vector =" ); 00988 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 00989 { 00990 fprintf( ioQQQ, "%9.1e", ionbal.RateRecomTot[nelem][ion] ); 00991 } 00992 fprintf( ioQQQ, "\n" ); 00993 00994 ContNegative(); 00995 ShowMe(); 00996 cdEXIT(EXIT_FAILURE); 00997 } 00998 00999 /* option to print ionization and recombination arrays 01000 * prt flag set with print array print arrays command */ 01001 if( prt.lgPrtArry[nelem] || lgPrintIt ) 01002 { 01003 /* say who we are, what we are doing .... */ 01004 fprintf( ioQQQ, 01005 "\n %s ion_solver DEBUG ion/rec rt [s-1] %s nz%.2f Te%.4e ne%.4e Tot abun:%.3e ion abun%.2e mole%.2e\n", 01006 elementnames.chElementSym[nelem], 01007 elementnames.chElementName[nelem], 01008 fnzone, 01009 phycon.te , 01010 dense.eden, 01011 dense.gas_phase[nelem], 01012 abund_total , 01013 dense.xMolecules[nelem] ); 01014 /* total ionization rate, all processes */ 01015 fprintf( ioQQQ, " %s Ioniz total " ,elementnames.chElementSym[nelem]); 01016 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01017 { 01018 fprintf( ioQQQ, " %9.2e", ionbal.RateIonizTot[nelem][ion] ); 01019 } 01020 fprintf( ioQQQ, "\n" ); 01021 01022 /* sinks from the chemistry network */ 01023 fprintf(ioQQQ," %s mole sink ", 01024 elementnames.chElementSym[nelem]); 01025 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01026 { 01027 fprintf(ioQQQ," %9.2e", mole.sink[nelem][ion] ); 01028 } 01029 fprintf( ioQQQ, "\n" ); 01030 01031 if( dynamics.lgAdvection ) 01032 { 01033 /* sinks from the dynamics */ 01034 fprintf(ioQQQ," %s dynm sink ", 01035 elementnames.chElementSym[nelem]); 01036 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01037 { 01038 fprintf(ioQQQ," %9.2e", dynamics.Rate ); 01039 } 01040 fprintf( ioQQQ, "\n" ); 01041 } 01042 01043 /* sum of all creation processes */ 01044 fprintf( ioQQQ, " %s Recom total " ,elementnames.chElementSym[nelem]); 01045 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01046 { 01047 fprintf( ioQQQ, " %9.2e", ionbal.RateRecomTot[nelem][ion] ); 01048 } 01049 fprintf( ioQQQ, "\n" ); 01050 01051 /* sources from the chemistry network */ 01052 fprintf(ioQQQ," %s mole source ", 01053 elementnames.chElementSym[nelem]); 01054 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01055 { 01056 fprintf(ioQQQ," %9.2e", mole.source[nelem][ion] ); 01057 } 01058 fprintf( ioQQQ, "\n" ); 01059 01060 if( dynamics.lgAdvection ) 01061 { 01062 /* source from the dynamcs */ 01063 fprintf(ioQQQ," %s dynm sour ", 01064 elementnames.chElementSym[nelem]); 01065 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01066 { 01067 fprintf(ioQQQ," %9.2e", dynamics.Source[nelem][ion] ); 01068 } 01069 fprintf( ioQQQ, "\n" ); 01070 } 01071 01072 /* collisional ionization */ 01073 fprintf( ioQQQ, " %s Coll ioniz " ,elementnames.chElementSym[nelem] ); 01074 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01075 { 01076 fprintf( ioQQQ, " %9.2e", ionbal.CollIonRate_Ground[nelem][ion][0] ); 01077 } 01078 fprintf( ioQQQ, "\n" ); 01079 01080 /* UTA ionization */ 01081 fprintf( ioQQQ, " %s UTA ioniz " ,elementnames.chElementSym[nelem] ); 01082 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01083 { 01084 fprintf( ioQQQ, " %9.2e", ionbal.UTA_ionize_rate[nelem][ion] ); 01085 } 01086 fprintf( ioQQQ, "\n" ); 01087 01088 /* photo ionization */ 01089 fprintf( ioQQQ, " %s Phot ioniz " ,elementnames.chElementSym[nelem]); 01090 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01091 { 01092 fprintf( ioQQQ, " %9.2e", 01093 ionbal.PhotoRate_Shell[nelem][ion][Heavy.nsShells[nelem][ion]-1][0] ); 01094 } 01095 fprintf( ioQQQ, "\n" ); 01096 01097 /* auger ionization */ 01098 fprintf( ioQQQ, " %s Auger ioniz " ,elementnames.chElementSym[nelem]); 01099 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01100 { 01101 fprintf( ioQQQ, " %9.2e", 01102 auger[ion] ); 01103 } 01104 fprintf( ioQQQ, "\n" ); 01105 01106 /* secondary ionization */ 01107 fprintf( ioQQQ, " %s Secon ioniz " ,elementnames.chElementSym[nelem]); 01108 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01109 { 01110 fprintf( ioQQQ, " %9.2e", 01111 secondaries.csupra[nelem][ion] ); 01112 } 01113 fprintf( ioQQQ, "\n" ); 01114 01115 /* grain ionization - not total rate but should be dominant process */ 01116 fprintf( ioQQQ, " %s ion on grn " ,elementnames.chElementSym[nelem]); 01117 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01118 { 01119 fprintf( ioQQQ, " %9.2e", gv.GrainChTrRate[nelem][ion][ion+1] ); 01120 } 01121 fprintf( ioQQQ, "\n" ); 01122 01123 /* chemistry ionization - not total rate but should be dominant process 01124 * not source or sink but just increase ionization within ladder */ 01125 fprintf( ioQQQ, " %s ion in chem " ,elementnames.chElementSym[nelem]); 01126 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01127 { 01128 fprintf( ioQQQ, " %9.2e", mole.xMoleChTrRate[nelem][ion][ion+1] ); 01129 } 01130 fprintf( ioQQQ, "\n" ); 01131 01132 /* charge exchange ionization */ 01133 fprintf( ioQQQ, " %s chr trn ion " ,elementnames.chElementSym[nelem] ); 01134 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01135 { 01136 /* sum has units s-1 */ 01137 double sum = atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[ipHELIUM][1]+ 01138 atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[ipHYDROGEN][1]; 01139 01140 if( nelem==ipHELIUM && ion==0 ) 01141 { 01142 sum += atmdat.HeCharExcIonTotal; 01143 } 01144 else if( nelem==ipHYDROGEN && ion==0 ) 01145 { 01146 sum += atmdat.HCharExcIonTotal; 01147 } 01148 fprintf( ioQQQ, " %9.2e", sum ); 01149 } 01150 fprintf( ioQQQ, "\n" ); 01151 01152 /* radiative recombination */ 01153 fprintf( ioQQQ, " %s radiati rec " ,elementnames.chElementSym[nelem]); 01154 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01155 { 01156 fprintf( ioQQQ, " %9.2e", dense.eden*ionbal.RR_rate_coef_used[nelem][ion] ); 01157 } 01158 fprintf( ioQQQ, "\n" ); 01159 01160 /* old DR recombination */ 01161 fprintf( ioQQQ, " %s dr old rec " ,elementnames.chElementSym[nelem]); 01162 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01163 { 01164 fprintf( ioQQQ, " %9.2e", dense.eden*ionbal.DR_old_rate_coef[nelem][ion] ); 01165 } 01166 fprintf( ioQQQ, "\n" ); 01167 01168 /* Badnell DR recombination */ 01169 fprintf( ioQQQ, " %s drBadnel rec" ,elementnames.chElementSym[nelem]); 01170 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01171 { 01172 fprintf( ioQQQ, " %9.2e", dense.eden*ionbal.DR_Badnell_rate_coef[nelem][ion] ); 01173 } 01174 fprintf( ioQQQ, "\n" ); 01175 01176 /* grain recombination - not total but from next higher ion, should 01177 * be dominant */ 01178 fprintf( ioQQQ, " %s rec on grn " ,elementnames.chElementSym[nelem]); 01179 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01180 { 01181 fprintf( ioQQQ, " %9.2e", gv.GrainChTrRate[nelem][ion+1][ion] ); 01182 } 01183 fprintf( ioQQQ, "\n" ); 01184 01185 /* grain due to chemistry - not total but from next higher ion, should 01186 * be dominant - not source or sink, just moves up or down */ 01187 fprintf( ioQQQ, " %s rec in chem " ,elementnames.chElementSym[nelem]); 01188 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01189 { 01190 fprintf( ioQQQ, " %9.2e", mole.xMoleChTrRate[nelem][ion+1][ion] ); 01191 } 01192 fprintf( ioQQQ, "\n" ); 01193 01194 /* charge exchange recombination */ 01195 fprintf( ioQQQ, " %s chr trn rec " ,elementnames.chElementSym[nelem]); 01196 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 01197 { 01198 double sum = 01199 atmdat.HCharExcRecTo[nelem][ion]* 01200 StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop* 01201 dense.xIonDense[ipHYDROGEN][1] + 01202 01203 atmdat.HeCharExcRecTo[nelem][ion]* 01204 StatesElem[ipHE_LIKE][ipHELIUM][ipHe1s1S].Pop* 01205 dense.xIonDense[ipHELIUM][1]; 01206 01207 if( nelem==ipHELIUM && ion==0 ) 01208 { 01209 sum += atmdat.HeCharExcRecTotal; 01210 } 01211 else if( nelem==ipHYDROGEN && ion==0 ) 01212 { 01213 sum += atmdat.HCharExcRecTotal; 01214 } 01215 fprintf( ioQQQ, " %9.2e", sum ); 01216 } 01217 fprintf( ioQQQ, "\n" ); 01218 01219 /* the "new" abundances the resulted from the previous ratio */ 01220 fprintf( ioQQQ, " %s Abun [cm-3] " ,elementnames.chElementSym[nelem] ); 01221 for( ion=0; ion <= dense.IonHigh[nelem]; ion++ ) 01222 { 01223 fprintf( ioQQQ, " %9.2e", dense.xIonDense[nelem][ion] ); 01224 } 01225 fprintf( ioQQQ, "\n" ); 01226 } 01227 01228 if( PARALLEL_MODE ) 01229 { 01230 free(ipiv); 01231 free(source); 01232 free(sink); 01233 free(xmat); 01234 free(xmatsave); 01235 free(auger); 01236 } 01237 return; 01238 } 01239 01240 /* 01241 01242 Solve an ionization level system with specified ionization and 01243 recombination rates between neighboring ions, and additional sink 01244 and source terms. The sink array is overwritten, and the results 01245 appear in the source array. 01246 01247 Written in matrix form, the algorithm is equivalent to the 01248 tridiagonal algorithm in Numerical Recipes applied to: 01249 01250 / i_0+a_0 -r_0 . . . \ / x_0 \ / s_0 \ 01251 | -i_0 i_1+a_1+r_0 -r_1 . . | | x_1 | | s_1 | 01252 | . -i_1 i_2+a_2+r_1 -r_2 . | | x_2 | | s_2 | 01253 | . . (etc....) | | ... | = | ... | 01254 \ . . . / \ / \ / 01255 01256 where i, r are the ionization and recombination rates, s is the 01257 source rate and a is the sink rate. 01258 01259 This matrix is diagonally dominant only when the sink terms are 01260 large -- the alternative method coded here prevents rounding error 01261 in the diagonal terms disturbing the solution when this is not the 01262 case. 01263 01264 */ 01265 01266 /* solveions tridiagonal solver but optimized for structure of balance matrix */ 01267 void solveions(double *ion, double *rec, double *snk, double *src, 01268 long int nlev, long int nmax) 01269 { 01270 double kap, bet; 01271 long int i; 01272 01273 DEBUG_ENTRY( "solveions()" ); 01274 01275 if(nmax != -1) 01276 { 01277 /* Singular case */ 01278 src[nmax] = 1.; 01279 for(i=nmax;i<nlev-1;i++) 01280 src[i+1] = src[i]*ion[i]/rec[i]; 01281 for(i=nmax-1;i>=0;i--) 01282 src[i] = src[i+1]*rec[i]/ion[i]; 01283 } 01284 else 01285 { 01286 i = 0; 01287 kap = snk[0]; 01288 for(i=0;i<nlev-1;i++) 01289 { 01290 bet = ion[i]+kap; 01291 if(bet == 0.) 01292 { 01293 fprintf(ioQQQ,"Ionization solver error\n"); 01294 cdEXIT(EXIT_FAILURE); 01295 } 01296 bet = 1./bet; 01297 src[i] *= bet; 01298 src[i+1] += ion[i]*src[i]; 01299 snk[i] = bet*rec[i]; 01300 kap = kap*snk[i]+snk[i+1]; 01301 } 01302 bet = kap; 01303 if(bet == 0.) 01304 { 01305 fprintf(ioQQQ,"Ionization solver error\n"); 01306 cdEXIT(EXIT_FAILURE); 01307 } 01308 src[i] /= bet; 01309 01310 for(i=nlev-2;i>=0;i--) 01311 { 01312 src[i] += snk[i]*src[i+1]; 01313 } 01314 } 01315 }