cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ion_solver.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*ion_solver solve the bi-diagonal matrix for ionization balance */
4 #include "cddefines.h"
5 #include "yield.h"
6 #include "prt.h"
7 #include "continuum.h"
8 #include "iso.h"
9 #include "dynamics.h"
10 #include "grainvar.h"
11 #include "hmi.h"
12 #include "mole.h"
13 #include "thermal.h"
14 #include "thirdparty.h"
15 #include "conv.h"
16 #include "secondaries.h"
17 #include "phycon.h"
18 #include "atmdat.h"
19 #include "heavy.h"
20 #include "elementnames.h"
21 #include "dense.h"
22 #include "radius.h"
23 #include "ionbal.h"
24 
25 void solveions(double *ion, double *rec, double *snk, double *src,
26  long int nlev, long int nmax);
27 
29  /* this is element on the c scale, H is 0 */
30  long int nelem,
31  /* option to print this element when called */
32  bool lgPrintIt)
33 {
34  /* use tridiag or general matrix solver? */
35  bool lgTriDiag=false;
36  long int ion,
37  limit,
38  IonProduced,
39  nej,
40  ns,
41  jmax=-1;
42  double totsrc;
43  double dennew, rateone;
44  bool lgNegPop;
45  static bool lgMustMalloc=true;
46  static double *sink, *source , *sourcesave;
47  int32 nerror;
48  static int32 *ipiv;
49  long ion_low, ion_range, i, j, ion_to , ion_from;
50  static double sum_dense;
51  /* only used for debug printout */
52  static double *auger;
53 
54  double abund_total, renormnew;
55  bool lgHomogeneous = true;
56  static double *xmat , *xmatsave;
57 
58  DEBUG_ENTRY( "ion_solver()" );
59 
60  /* this is on the c scale, so H is 0 */
61  ASSERT( nelem >= 0);
62  ASSERT( dense.IonLow[nelem] >= 0 );
63  ASSERT( dense.IonHigh[nelem] >= 0 );
64 
65  /* H is special because its abundance spills into three routines -
66  * the ion/atom solver (this routine), the H-mole solvers (hmole), and
67  * the heavy molecule solver. xmolecules only includes the heavy mole
68  * part for H only. So the difference between gas_phase and xmolecules
69  * includes the H2 part of the molecular network. This branch does
70  * this special H case, then the general case (heavy elements are
71  * not part of the H2 molecular network) */
72 
73  /* >>chng 01 dec 07, define abund_total, total atom and ion abundance
74  * here, removing molecules */
75  if( nelem == ipHYDROGEN )
76  {
77  /* Hydrogen is a special case since hmole does not include the
78  * atom/molecules - hmole collapses its network into H = H0 + H+
79  * and then forces the solution determined here for partitioning
80  * between these two */
81  abund_total = dense.xIonDense[nelem][0] + dense.xIonDense[nelem][1];
82  }
83  else
84  {
85  abund_total = SDIV( dense.gas_phase[nelem] - dense.xMolecules[nelem] );
86  }
87 
88  /* >>chng 04 nov 22,
89  * if gas nearly all atomic/ionic do not let source - sink terms from
90  * molecular network force system of balance equations to become
91  * inhomogeneous
92  * what constitutes a source or sink IS DIFFERENT for hydrogen and the rest
93  * the H solution must couple with hmole - and its defn of source and sink. For instance, oxygen charge
94  * transfer goes into source and sink terms for hydrogen. So we never hose source and sink for H.
95  * for the heavy elements, which couple onto comole, mole.source and sink represent terms that
96  * remove atoms and ions from the ionization ladder. their presence makes the system of
97  * equations inhomogeneous. we don't want to do this when comole has a trivial effect on
98  * the ionization balance since the matrix becomes unstable */
99  /* >>chng 04 dec 06, limit from 0.01 to 1e-10 as per NA suggestion */
100  if( nelem>ipHYDROGEN && dense.xMolecules[nelem]/SDIV(dense.gas_phase[nelem]) < 1e-10 )
101  {
102  for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
103  {
104  mole.source[nelem][ion] = 0.;
105  mole.sink[nelem][ion] = 0.;
106  }
107  }
108 
109  /* protect against case where all gas phase abundances are in molecules, use the
110  * atomic and first ion density from the molecule solver
111  * >>chng 04 aug 15, NA change from 10 0000 to 10 pre-coef on
112  * FLT_EPSILON for stability in PDR
113  * the factor 10.*FLT_EPSILON also appears in mole_h_step in total H
114  * conservation */
115  if( fabs( dense.gas_phase[nelem] - dense.xMolecules[nelem])/SDIV(dense.gas_phase[nelem] ) <
116  10.*FLT_EPSILON )
117  {
118  double renorm;
119  /* >>chng 04 jul 31, add logic to conserve nuclei in fully molecular limit;
120  * in first calls, when searching for solution, we may be VERY far
121  * off, and sum of first ion and atom density may be far larger
122  * than difference between total gas and molecular densities,
123  * since they reflect the previous evaluation of the solution. Do
124  * renorm to cover this case */
125  /* first form sum of all atoms and ions */
126  realnum sum = 0.;
127  for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
128  sum += dense.xIonDense[nelem][ion];
129  /* now renorm to this sum - this should be unity, and is not if we have
130  * now conserved particles, due to molecular fraction changing */
131  renorm = dense.gas_phase[nelem] / SDIV(sum + dense.xMolecules[nelem] );
132 
133  abund_total = renorm * sum;
134  }
135 
136  /* negative residual density */
137  if( abund_total < 0. )
138  {
139  /* >>chng 05 dec 16, do not abort when net atomic/ ionic abundance
140  * is negative due to chem net having too much of a species - this
141  * happens naturally when ices become well formed, but the code will
142  * converge away from it. Rather, we set the ionization
143  * convergence flag to "not converge" and soldier on
144  * if negative populations do not go away, we will eventually
145  * terminate due to convergence failures */
146  /* print error if not search */
147  if(!conv.lgSearch )
148  fprintf(ioQQQ,
149  "PROBLEM ion_solver - neg net atomic abundance zero for nelem= %li, rel val= %.2e conv.nTotalIoniz=%li, fixed\n",
150  nelem,
151  fabs(abund_total) / SDIV(dense.xMolecules[nelem]),
152  conv.nTotalIoniz );
153  /* fix up is to use half the positive abundance, assuming chem is
154  * trying to get too much of this species */
155  abund_total = -abund_total/2.;
156 
157  /* say that ionization is not converged - do not abort - but if
158  * cannot converge away from negative solution, this will become a
159  * convergence failure abort */
160  conv.lgConvIoniz = false;
161  strcpy( conv.chConvIoniz, "neg ion" );
162  }
163 
164  /* return if IonHigh is zero, since no ionization at all */
165  if( dense.IonHigh[nelem] == 0 )
166  {
167  /* set the atom to the total gas phase abundance */
168  dense.xIonDense[nelem][0] = (realnum)abund_total;
169  return;
170  }
171 
172  /* >>chng 01 may 09, add option to force ionization distribution with
173  * element name ioniz */
174  if( dense.lgSetIoniz[nelem] )
175  {
176  for( ion=0; ion<nelem+2; ++ion )
177  {
178  dense.xIonDense[nelem][ion] = dense.SetIoniz[nelem][ion]*
179  (realnum)abund_total;
180  }
181  return;
182  }
183 
184  /* impossible for HIonFrac[nelem] to be zero if IonHigh(nelem)=nelem+1
185  * HIonFrac(nelem) is stripped to hydrogen */
186  /* >>chng 01 oct 30, to assert */
187  ASSERT( (dense.IonHigh[nelem] < nelem + 1) ||
188  iso.pop_ion_ov_neut[ipH_LIKE][nelem] > 0. );
189 
190  /* zero out the ionization and recombination rates that we will modify here,
191  * but not the iso-electronic stages which are done elsewhere,
192  * the nelem stage of ionization is he-like,
193  * the nelem+1 stage of ionization is h-like */
194 
195  /* loop over stages of ionization that we solve for here,
196  * up through and including one less than nelem-NISO,
197  * never actually do highest NISO stages of ionization since they
198  * come from the ionization ratio from the next lower stage */
199  limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);
200 
201  /* the full range of ionization - this is number of ionization stages */
202  ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
203  ASSERT( ion_range <= nelem+2 );
204 
205  ion_low = dense.IonLow[nelem];
206 
207  /* PARALLEL_MODE true if there can be several instances of this routine
208  * running in the same memory pool - we must have fresh instances of the
209  * arrays for each */
210  if( PARALLEL_MODE || lgMustMalloc )
211  {
212  long int ncell=LIMELM+1;
213  lgMustMalloc = false;
214  if( PARALLEL_MODE )
215  ncell = ion_range;
216 
217  /* this will be "new" matrix, with non-adjacent coupling included */
218  xmat=(double*)MALLOC( (sizeof(double)*(unsigned)(ncell*ncell) ));
219  xmatsave=(double*)MALLOC( (sizeof(double)*(unsigned)(ncell*ncell) ));
220  sink=(double*)MALLOC( (sizeof(double)*(unsigned)ncell ));
221  source=(double*)MALLOC( (sizeof(double)*(unsigned)ncell ));
222  sourcesave=(double*)MALLOC( (sizeof(double)*(unsigned)ncell ));
223  ipiv=(int32*)MALLOC( (sizeof(int32)*(unsigned)ncell ));
224  /* auger is used only for debug printout - it is special because with multi-electron
225  * Auger ejection, very high stages of ionization can be produced, well beyond
226  * the highest stage considered here. so we malloc to the limit */
227  auger=(double*)MALLOC( (sizeof(double)*(unsigned)(LIMELM+1) ));
228  }
229 
230  for( i=0; i<ion_range;i++ )
231  {
232  source[i] = 0.;
233  }
234 
235  /* this will be used to address the 2d arrays */
236 # undef MAT
237 # define MAT(M_,I_,J_) (*((M_)+(I_)*(ion_range)+(J_)))
238 
239  /* zero-out loop comes before main loop since there are off-diagonal
240  * elements in the main ionization loop, due to multi-electron processes,
241  * TotIonizRate and TotRecom were already set in h-like and he-like solvers
242  * other recombination rates were already set by routines responsible for them */
243  for( ion=0; ion <= limit; ion++ )
244  {
245  ionbal.RateIonizTot[nelem][ion] = 0.;
246  }
247 
248  /* auger is used only for debug printout - it is special because with multi-electron
249  * Auger ejection, very high stages of ionization can be produced, well beyond
250  * the highest stage considered here. so we malloc to the limit */
251  for( i=0; i< LIMELM+1; ++i )
252  {
253  auger[i] = 0.;
254  }
255 
256  /* zero out xmat */
257  for( i=0; i< ion_range; ++i )
258  {
259  for( j=0; j< ion_range; ++j )
260  {
261  MAT( xmat, i, j ) = 0.;
262  }
263  }
264 
265  {
266  /* this sets up a fake ionization balance problem, with a trivial solution,
267  * for debugging the ionization solver */
268  enum {DEBUG_LOC=false};
269  if( DEBUG_LOC && nelem==ipCARBON )
270  {
271  broken();/* within optional debug print statement */
272  dense.IonLow[nelem] = 0;
273  dense.IonHigh[nelem] = 3;
274  abund_total = 1.;
275  ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
276  /* make up ionization and recombination rates */
277  for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
278  {
279  double fac=1;
280  if(ion)
281  fac = 1e-10;
282  ionbal.RateRecomTot[nelem][ion] = 100.;
283  for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
284  {
285  /* direct photoionization of this shell */
286  ionbal.PhotoRate_Shell[nelem][ion][ns][0] = fac;
287  }
288  }
289  }
290  }
291 
292  /* Now put in all recombination and ionization terms from CO_mole() that
293  * come from molecular reactions. this traces molecular process that
294  * change ionization stages with this ladder - but do not remove from
295  * the ladder */
296  for( ion_to=dense.IonLow[nelem]; ion_to <= limit; ion_to++ )
297  {
298  for( ion_from=dense.IonLow[nelem]; ion_from <= dense.IonHigh[nelem]; ++ion_from )
299  {
300  /* do not do ion onto itself */
301  if( ion_to != ion_from )
302  {
303  /* this is the rate coefficient for charge transfer from ion to ion_to */
304  /*rateone = gv.GrainChTrRate[nelem][ion_from][ion_to];*/
305  rateone = mole.xMoleChTrRate[nelem][ion_from][ion_to];
306  ASSERT( rateone >= 0. );
307  MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
308  MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
309  }
310  }
311  }
312 
313  /* now get actual arrays of ionization and recombination processes,
314  * but only for the ions that are done as two-level systems */
315  /* in two-stage system, atom + first ion, limit is zero but must
316  * include gv.GrainChTrRate[nelem][1][0] */
317  /* grain charge transfer */
319  {
320  long int low;
321  /* do not double count this process for atoms that are in the co network - we use
322  * a net recombination coefficient derived from the co solution, this includes grain ct */
323  /* >>chng 05 dec 23, add mole.lgElem_in_chemistry */
324  if( mole.lgElem_in_chemistry[nelem] )
325  /*if( nelem==ipHYDROGEN ||nelem==ipCARBON ||nelem== ipOXYGEN ||nelem==ipSILICON ||nelem==ipSULPHUR
326  ||nelem==ipNITROGEN ||nelem==ipCHLORINE )*/
327  {
328  low = MAX2(1, dense.IonLow[nelem] );
329  }
330  else
331  low = dense.IonLow[nelem];
332 
333  for( ion_to=low; ion_to <= dense.IonHigh[nelem]; ion_to++ )
334  {
335  for( ion_from=dense.IonLow[nelem]; ion_from <= dense.IonHigh[nelem]; ++ion_from )
336  {
337  /* do not do ion onto itself */
338  if( ion_to != ion_from )
339  {
340  /* this is the rate coefficient for charge transfer from ion to ion_to
341  * both grain charge transfer ionization and recombination */
342  rateone = gv.GrainChTrRate[nelem][ion_from][ion_to];
343  MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
344  MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
345  }
346  }
347  }
348  }
349 
350  for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
351  {
352  /* thermal & secondary collisional ionization */
353  rateone = ionbal.CollIonRate_Ground[nelem][ion][0] +
354  secondaries.csupra[nelem][ion] +
355  /* inner shell ionization by UTA lines */
356  ionbal.UTA_ionize_rate[nelem][ion];
357  ionbal.RateIonizTot[nelem][ion] += rateone;
358 
359  /* UTA ionization */
360  if( ion+1-ion_low < ion_range )
361  {
362  /* depopulation processes enter with negative sign */
363  MAT( xmat, ion-ion_low, ion-ion_low ) -= rateone;
364  MAT( xmat, ion-ion_low, ion+1-ion_low ) += rateone;
365  }
366 
367  /* total recombination rate */
368  if( ion-1-ion_low >= 0 )
369  {
370  /* loss of this ion due to recombination to next lower ion stage */
371  MAT( xmat,ion-ion_low, ion-ion_low ) -= ionbal.RateRecomTot[nelem][ion-1];
372  MAT( xmat,ion-ion_low, ion-1-ion_low ) += ionbal.RateRecomTot[nelem][ion-1];
373  }
374 
375  /* loop over all atomic sub-shells to include photoionization */
376  for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
377  {
378  /* direct photoionization of this shell */
379  ionbal.RateIonizTot[nelem][ion] += ionbal.PhotoRate_Shell[nelem][ion][ns][0];
380 
381  /* this is the primary ionization rate - add to diagonal element,
382  * test on ion stage is so that we don't include ionization from the very highest
383  * ionization stage to even higher - since those even higher stages are not considered
384  * this would appear as a sink - but populations of this highest level is ensured to
385  * be nearly trivial and neglecting it production of even higher ionization OK */
386  /* >>chng 04 nov 29 RJRW, include following in this branch so only
387  * evaluated when below ions done with iso-sequence */
388  if( ion+1-ion_low < ion_range )
389  {
390  /* this will be redistributed into charge states in following loop */
391  MAT( xmat, ion-ion_low, ion-ion_low ) -= ionbal.PhotoRate_Shell[nelem][ion][ns][0];
392 
393  /* t_yield::Inst().nelec_eject(nelem,ion,ns) is total number of electrons that can
394  * possibly be freed
395  * loop over nej, the number of electrons ejected including the primary,
396  * nej = 1 is primary, nej > 1 includes primary plus Auger
397  * t_yield::Inst().elec_eject_frac is prob of nej electrons */
398  for( nej=1; nej <= t_yield::Inst().nelec_eject(nelem,ion,ns); nej++ )
399  {
400  /* this is the ion that is produced by this ejection,
401  * limited by highest possible stage of ionization -
402  * do not want to ignore ionization that go beyond this */
403  IonProduced = MIN2(ion+nej,dense.IonHigh[nelem]);
404  rateone = ionbal.PhotoRate_Shell[nelem][ion][ns][0]*
405  t_yield::Inst().elec_eject_frac(nelem,ion,ns,nej-1);
406  /* >>chng 04 sep 06, above had included factor of nej to get rate, but
407  * actually want events into particular ion */
408  /* number of electrons ejected
409  *(double)nej; */
410  /* it goes into this charge state - recall upper cap due to ion stage trimming
411  * note that compensating loss term on diagonal was done before this
412  * loop, since frac_elec_eject adds to unity */
413  MAT( xmat, ion-ion_low, IonProduced-ion_low ) += rateone;
414 
415  /* only used for possible printout - multiple electron Auger rate -
416  * do not count one-electron as Auger */
417  if( nej>1 )
418  auger[IonProduced-1] += rateone;
419  }
420  }
421  }
422 
423  /* this is charge transfer ionization of this species by hydrogen and helium */
424  rateone =
425  atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[ipHELIUM][1]+
427  ionbal.RateIonizTot[nelem][ion] += rateone;
428  if( ion+1-ion_low < ion_range )
429  {
430  MAT( xmat, ion-ion_low, ion-ion_low ) -= rateone;
431  MAT( xmat, ion-ion_low, ion+1-ion_low ) += rateone;
432  }
433  }
434 
435  /* after this loop, ionbal.RateIonizTot and ionbal.RateRecomTot have been defined for the
436  * stages of ionization that are done with simple */
437  /* begin loop at first species that is treated with full model atom */
438  j = MAX2(0,limit+1);
439  /* possible that lowest stage of ionization is higher than this */
440  j = MAX2( ion_low , j );
441  for( ion=j; ion<=dense.IonHigh[nelem]; ion++ )
442  {
443  ASSERT( ion>=0 && ion<nelem+2 );
444  /* use total ionization/recombination rates for species done with ISO solver */
445  if( ion+1-ion_low < ion_range )
446  {
447  /* depopulation processes enter with negative sign */
448  MAT( xmat, ion-ion_low, ion-ion_low ) -= ionbal.RateIonizTot[nelem][ion];
449  MAT( xmat, ion-ion_low, ion+1-ion_low ) += ionbal.RateIonizTot[nelem][ion];
450  }
451 
452  if( ion-1-ion_low >= 0 )
453  {
454  /* loss of this ion due to recombination to next lower ion stage */
455  MAT( xmat,ion-ion_low, ion-ion_low ) -= ionbal.RateRecomTot[nelem][ion-1];
456  MAT( xmat,ion-ion_low, ion-1-ion_low ) += ionbal.RateRecomTot[nelem][ion-1];
457  }
458  }
459 
460  /* this will be sum of source and sink terms, will be used to decide if
461  * matrix is singular */
462  totsrc = 0.;
463  /*>>chng 06 nov 28 only include source from molecules if we have an estimated first
464  * solution - first test is that we have called mole at least twice,
465  * second test is that we are on a later iteration */
466  if( conv.nTotalIoniz > 1 || iteration > 1 )
467  {
468 
469  for( i=0; i<ion_range;i++ )
470  {
471  ion = i+ion_low;
472 
473  /* these are the external source and sink terms */
474  /* source first */
475  /* need negative sign to get positive pops */
476  source[i] -= mole.source[nelem][ion];
477 
478  totsrc += mole.source[nelem][ion];
479  /* sink next */
480  /*MAT( xmat, i, i ) += mole.sink[nelem][ion];*/
481  MAT( xmat, i, i ) -= mole.sink[nelem][ion];
482  }
483  }
484 
485  /* matrix is not homogeneous if source is non-zero */
486  if( totsrc != 0. )
487  lgHomogeneous = false;
488 
489  /* chng 03 jan 13 rjrw, add in dynamics if required here,
490  * last test - only do advection if we have not overrun the radius scale */
492  /*&& radius.depth < dynamics.oldFullDepth*/ )
493  {
494  for( i=0; i<ion_range;i++ )
495  {
496  ion = i+ion_low;
497  /*MAT( xmat, i, i ) += dynamics.Rate;*/
498  MAT( xmat, i, i ) -= dynamics.Rate;
499  /*src[i] += dynamics.Source[nelem][ion];*/
500  source[i] -= dynamics.Source[nelem][ion];
501  /* fprintf(ioQQQ," %li %li %.3e (%.3e %.3e)\n",i,i,MAT(xmat,i,i),
502  dynamics.Rate, dynamics.Source[nelem][ion]);*/
503  }
504  lgHomogeneous = false;
505  }
506 
507  for( i=0; i< ion_range; ++i )
508  {
509  for( j=0; j< ion_range; ++j )
510  {
511  MAT( xmatsave, i, j ) = MAT( xmat, i, j );
512  }
513  sourcesave[i] = source[i];
514  }
515 
516  /* >>chng 06 nov 21, for very high ionization parameter sims the H molecular
517  * fraction can become so small that atom = atom + molecule. In this case we
518  * will not count system as an inhomogeneous case since linear algebra package
519  * will fail. If molecules are very small, count as homogeneous.
520  * Note that we have already added sink terms to the main matrix and they
521  * will not be removed, a possible source of error, but they must not
522  * have been significant, given that the molecular fraction is so small */
523  if( !lgHomogeneous && ion_range==2 )
524  {
525  /* solve algebraically */
526  double a = MAT( xmatsave, 0, 0 ),
527  b = MAT( xmatsave, 1, 0 ) ,
528  c = MAT( xmatsave, 0, 1 ) ,
529  d = MAT( xmatsave, 1, 1 );
530  b = SDIV(b);
531  d = SDIV(d);
532  double ratio1 = a/b , ratio2 = c/d , fratio1=fabs(a/b),fratio2=fabs(c/d);
533  if( fabs(ratio1-ratio2)/MAX2(fratio1,fratio2) <DBL_EPSILON*1e4 )
534  {
535  abund_total = SDIV( dense.gas_phase[nelem] - dense.xMolecules[nelem] );
536  lgHomogeneous = true;
537  }
538  }
539 # if 0
540  if( nelem==ipHYDROGEN &&
541  fabs(dense.xMolecules[nelem]) / SDIV(dense.gas_phase[ipHYDROGEN]) <DBL_EPSILON*100. )
542  {
543  abund_total = SDIV( dense.gas_phase[nelem] - dense.xMolecules[nelem] );
544  lgHomogeneous = true;
545  }
546 # endif
547 
548  /* this is true if no source terms
549  * we will use total population and species conservation to replace one
550  * set of balance equations since overdetermined */
551  if( lgHomogeneous )
552  {
553  double rate_ioniz=1., rate_recomb /*, scale = 0.*/;
554  /* Simple estimate of most abundant ion */
555  jmax = 0;
556  for( i=0; i<ion_range-1;i++)
557  {
558  ion = i+ion_low;
559  rate_ioniz *= ionbal.RateIonizTot[nelem][ion];
560  rate_recomb = ionbal.RateRecomTot[nelem][ion];
561  /* trips if ion rate zero, so ll the gas will be more neutral than this */
562  if( rate_ioniz == 0)
563  break;
564  /* rec rate is zero */
565  if( rate_recomb <= 0.)
566  break;
567 
568  rate_ioniz /= rate_recomb;
569  if( rate_ioniz > 1.)
570  {
571  /* this is peak ionization stage */
572  jmax = i;
573  rate_ioniz = 1.;
574  }
575  }
576 
577  /* replace its matrix elements with population sum */
578  for( i=0; i<ion_range;i++ )
579  {
580  MAT(xmat,i,jmax) = 1.;
581  }
582  source[jmax] = abund_total;
583  }
584 
585  for( i=0; i< ion_range; ++i )
586  {
587  for( j=0; j< ion_range; ++j )
588  {
589  MAT( xmatsave, i, j ) = MAT( xmat, i, j );
590  }
591  sourcesave[i] = source[i];
592  }
593 
594  if( false && nelem == ipHYDROGEN && dynamics.lgAdvection&& iteration>1 )
595  {
596  fprintf(ioQQQ,"DEBUGG Rate %.2f %.3e \n",fnzone,dynamics.Rate);
597  fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateIonizTot[nelem][0], ionbal.RateIonizTot[nelem][1]);
598  fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateRecomTot[nelem][0], ionbal.RateRecomTot[nelem][1]);
599  fprintf(ioQQQ," %.3e %.3e %.3e\n\n", dynamics.Source[nelem][0], dynamics.Source[nelem][1], dynamics.Source[nelem][2]);
600  }
601 
602  {
603  /* option to print matrix */
604  /*@-redef@*/
605  enum {DEBUG_LOC=false};
606  /*@+redef@*/
607  if( DEBUG_LOC && nzone==1 && lgPrintIt )
608  {
609  fprintf( ioQQQ,
610  " DEBUG ion_solver: nelem=%li ion_range=%li, limit=%li, nConv %li xmat follows\n",
611  nelem , ion_range,limit , conv.nTotalIoniz );
612  if( lgHomogeneous )
613  fprintf(ioQQQ , "Homogeneous \n");
614  for( i=0; i<ion_range; ++i )
615  {
616  for( j=0;j<ion_range;j++ )
617  {
618  fprintf(ioQQQ,"%e\t",MAT(xmatsave,i,j));
619  }
620  fprintf(ioQQQ,"\n");
621  }
622  fprintf(ioQQQ,"source follows\n");
623  for( i=0; i<ion_range;i++ )
624  {
625  fprintf(ioQQQ,"%e\t",source[i]);
626  }
627  fprintf(ioQQQ,"\n");
628  }
629  }
630 
631  /* first branch is tridiag, following branch is just general matrix solver */
632  if( lgTriDiag )
633  {
634  bool lgLapack=false , lgTriOptimized=true;
635  /* there are two choices for the tridiag solver */
636  if(lgLapack)
637  {
638  /* this branch uses lapack tridiag solver, and should work
639  * it is hardwired to not be used because not extensively tested
640  * issues - is this slower than others, and is it robust in
641  * singular cases? */
642  bool lgBad = false;
643  /* Use Lapack tridiagonal solver */
644  double *dl, *d, *du;
645 
646  d = (double *) MALLOC((unsigned)ion_range*sizeof(double));
647  du = (double *) MALLOC((unsigned)(ion_range-1)*sizeof(double));
648  dl = (double *) MALLOC((unsigned)(ion_range-1)*sizeof(double));
649 
650  for(i=0;i<ion_range-1;i++)
651  {
652  du[i] = MAT(xmat,i+1,i);
653  d[i] = MAT(xmat,i,i);
654  dl[i] = MAT(xmat,i,i+1);
655  }
656  d[i] = MAT(xmat,i,i);
657 
658 # if 0
659  int i1, i2;
660  i1 = ion_range;
661  i2 = 1;
662  lgBad = dgtsv_wrapper(&i1, &i2, dl, d, du, source, &i2);
663 # endif
664  if( lgBad )
665  fprintf(ioQQQ," dgtsz error.\n");
666 
667  free(dl);free(du);free(d);
668  }
669  else if(lgTriOptimized)
670  {
671  /* Use tridiagonal solver re-coded to avoid rounding errors
672  * on diagonal -- uses determination of jmax for the
673  * singular case, but is otherwise independent of the array
674  * filling code above */
675 
676  for(i=0;i<ion_range;i++)
677  {
678  source[i] = sink[i] = 0.;
679  }
680  if(nelem == ipHYDROGEN && !hmi.lgNoH2Mole)
681  {
682  for(i=0;i<ion_range;i++)
683  {
684  ion = i+ion_low;
685  source[i] += mole.source[nelem][ion];
686  sink[i] += mole.sink[nelem][ion];
687  }
688  }
690  {
691  for(i=0;i<ion_range;i++)
692  {
693  ion = i+ion_low;
694  source[i] += dynamics.Source[nelem][ion];
695  sink[i] += dynamics.Rate;
696  }
697  }
698 
699  solveions(ionbal.RateIonizTot[nelem]+ion_low,ionbal.RateRecomTot[nelem]+ion_low,
700  sink,source,ion_range,jmax);
701  }
702  }
703  else
704  {
705  /* this is the default solver - now get new solution */
706  nerror = 0;
707  /* Use general matrix solver */
708  getrf_wrapper(ion_range, ion_range, xmat, ion_range, ipiv, &nerror);
709  if( nerror != 0 )
710  {
711  fprintf( ioQQQ,
712  " DISASTER ion_solver: dgetrf finds singular or ill-conditioned matrix nelem=%li %s ion_range=%li, limit=%li, nConv %li xmat follows\n",
713  nelem ,
714  elementnames.chElementSym[nelem],
715  ion_range,
716  limit ,
717  conv.nTotalIoniz );
718  for( i=0; i<ion_range; ++i )
719  {
720  for( j=0;j<ion_range;j++ )
721  {
722  fprintf(ioQQQ,"%e\t",MAT(xmatsave,j,i));
723  }
724  fprintf(ioQQQ,"\n");
725  }
726  fprintf(ioQQQ,"source follows\n");
727  for( i=0; i<ion_range;i++ )
728  {
729  fprintf(ioQQQ,"%e\t",source[i]);
730  }
731  fprintf(ioQQQ,"\n");
732  cdEXIT(EXIT_FAILURE);
733  }
734  getrs_wrapper('N', ion_range, 1, xmat, ion_range, ipiv, source, ion_range, &nerror);
735  if( nerror != 0 )
736  {
737  fprintf( ioQQQ, " DISASTER ion_solver: dgetrs finds singular or ill-conditioned matrix nelem=%li ionrange=%li\n",
738  nelem , ion_range );
739  cdEXIT(EXIT_FAILURE);
740  }
741  }
742 
743  {
744  /*@-redef@*/
745  /* this is to debug following failed assert */
746  enum {DEBUG_LOC=false};
747  /*@+redef@*/
748  if( DEBUG_LOC && nelem == ipHYDROGEN )
749  {
750  fprintf(ioQQQ,"debuggg ion_solver1 \t%.2f\t%.4e\t%.4e\tIon\t%.3e\tRec\t%.3e\n",
751  fnzone,
752  phycon.te,
753  dense.eden,
754  ionbal.RateIonizTot[nelem][0] ,
755  ionbal.RateRecomTot[nelem][0]);
756  fprintf(ioQQQ," Msrc %.3e %.3e\n", mole.source[ipHYDROGEN][0], mole.source[ipHYDROGEN][1]);
757  fprintf(ioQQQ," Msnk %.3e %.3e\n", mole.sink[ipHYDROGEN][0], mole.sink[ipHYDROGEN][1]);
758  fprintf(ioQQQ," Poprat %.3e nomol %.3e\n",source[1]/source[0],
759  ionbal.RateIonizTot[nelem][0]/ionbal.RateRecomTot[nelem][0]);
760  }
761  }
762 
763  for( i=0; i<ion_range; i++ )
764  {
765  /* is this blows, source[i] is NaN */
766  ASSERT( !isnan( source[i] ) );
767  }
768 
769 # define RJRW 0
770  if( RJRW && 0 )
771  {
772  /* verify that the rates are sensible */
773  double test;
774  for(i=0; i<ion_range; i++) {
775  test = 0.;
776  for(j=0; j<ion_range; j++) {
777  test = test+source[j]*MAT(xmatsave,j,i);
778  }
779  fprintf(ioQQQ,"%e\t",test);
780  }
781  fprintf(ioQQQ,"\n");
782 
783  test = 0.;
784  fprintf(ioQQQ," ion %li abundance %.3e\n",nelem,abund_total);
785  for( ion=dense.IonLow[nelem]; ion < dense.IonHigh[nelem]; ion++ )
786  {
787  if( ionbal.RateRecomTot[nelem][ion] != 0 && source[ion-ion_low] != 0 )
788  fprintf(ioQQQ," %li %.3e %.3e : %.3e\n",ion,source[ion-ion_low],
789  source[ion-ion_low+1]/source[ion-ion_low],
790  ionbal.RateIonizTot[nelem][ion]/ionbal.RateRecomTot[nelem][ion]);
791  else
792  fprintf(ioQQQ," %li %.3e [One ratio infinity]\n",ion,source[ion-ion_low]);
793  test += source[ion-ion_low];
794  }
795  test += source[ion-ion_low];
796  }
797 
798  /*
799  * >> chng 03 jan 15 rjrw:- terms are now included for
800  * molecular sources and sinks of H and H+.
801  *
802  * When the network is not in equilibrium, this will lead to a
803  * change in the derived abundance of H and H+ when after the
804  * matrix solution -- the difference between `renorm' and 1. is a
805  * measure of the quality of the solution (it will be 1. if the
806  * rate of transfer into Ho/H+ balances the rate of transfer
807  * out, for the consistent relative abundances).
808  *
809  * We therefore renormalize to keep the total H abundance
810  * correct -- only the molecular network is allowed to change
811  * this.
812  *
813  * To do this, only the ion abundances are corrected, as the
814  * molecular abundances may depend on several different
815  * conserved species.
816  *
817  */
818  if( lgHomogeneous )
819  {
820  dense.xIonDense[nelem][ion_low] = (realnum)abund_total;
821  for( i=1;i < ion_range; i++ )
822  {
823  dense.xIonDense[nelem][i+ion_low] = 0.;
824  }
825  }
826 
827  renormnew = 1.;
828  /* >>chng 06 mar 17, comment out test on old full depth - keep old solution if overrun scale */
829  if(iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection /*&& radius.depth < dynamics.oldFullDepth */ &&
830  nelem == ipHYDROGEN && hmi.lgNoH2Mole)
831  {
832  /* The normalization out of the matrix solution is correct and
833  * should be retained if: dynamics is on and the total
834  * abundance of HI & H+ isn't being controlled by the
835  * molecular network */
836  renormnew = 1.;
837  }
838  else
839  {
840  dennew = 0.;
841  sum_dense = 0.;
842 
843  /* find total population to renorm - also here check that negative pops did not occur */
844  for( i=0;i < ion_range; i++ )
845  {
846  ion = i+ion_low;
847  sum_dense += dense.xIonDense[nelem][ion];
848  dennew += source[i];
849  }
850 
851  if( dennew > 0.)
852  {
853  renormnew = sum_dense / dennew;
858  }
859  else
860  {
861  renormnew = 1.;
862  }
863  }
864  /* check not negative, should be +ve, can be zero if species has become totally molecular.
865  * this happens for hydrogen if no cosmic rays, or cr ion set very low */
866  if( renormnew < 0)
867  {
868  fprintf(ioQQQ,"PROBLEM impossible value of renorm \n");
869  }
870  ASSERT( renormnew>=0 );
871 
872  /* source[i] contains new solution for ionization populations
873  * save resulting abundances into main ionization density array,
874  * while checking whether any negative level populations occurred */
875  lgNegPop = false;
876  for( i=0; i < ion_range; i++ )
877  {
878  ion = i+ion_low;
879 
880  /*fprintf(ioQQQ," %li %li %.3e %.3e\n",nelem,ion,src[ion-ion_low+1],src[ion-ion_low]);
881  pop_ion_ov_neut[ion] = src[ion-ion_low+1]/src[ion-ion_low]; */
882  if( source[i] < 0. )
883  {
884  /* >>chng 04 dec 04, put test on neg abund here, don't print uncles value is very -ve */
885  /* >>chng 06 feb 28, from -1e-10 to -1e-9, sim func_t10 had several negative
886  * during initial search, due to extremely high ionization */
887  /* >>chng 06 mar 11, from 1e-9 to 2e-9 make many struc elements floats from double */
888  if( source[i]<-2e-9 )
889  fprintf(ioQQQ,
890  " PROBLEM negative ion population in ion_solver, nelem=%li, %s ion=%li val=%.3e Search?%c zone=%li iteration=%li\n",
891  nelem ,
892  elementnames.chElementSym[nelem],
893  ion ,
894  source[i] ,
895  TorF(conv.lgSearch) ,
896  nzone ,
897  iteration );
898  source[i] = 0.;
899  /* if this is one of the iso seq model atoms then must also zero out Lya and Pop2Ion */
900  if( ion == nelem+1-NISO )
901  {
902  long int ipISO = nelem - ion;
903  ASSERT( ipISO>=0 && ipISO<NISO );
904  StatesElem[ipISO][nelem][0].Pop = 0.;
905  }
906  }
907 
908  /* use new solution */
909  dense.xIonDense[nelem][ion] = (realnum)(source[i]*renormnew);
910  }
911 
912  /* Zero levels with abundances < 1e-25 which which will suffer numerical noise */
913  while( dense.IonHigh[nelem] > dense.IonLow[nelem] &&
914  dense.xIonDense[nelem][dense.IonHigh[nelem]] < 1e-25*abund_total )
915  {
916  ASSERT( dense.xIonDense[nelem][dense.IonHigh[nelem]] >= 0. );
917  /* zero out abundance and heating due to stage of ionization we are about to zero out */
918  dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
919  thermal.heating[nelem][dense.IonHigh[nelem]-1] = 0.;
920  /* decrement counter */
921  --dense.IonHigh[nelem];
922  }
923 
924  /* sanity check, either offset stages of low and high ionization,
925  * or no ionization at all */
926  ASSERT( (dense.IonLow[nelem] < dense.IonHigh[nelem]) ||
927  (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) );
928 
929  /* this should not happen */
930  if( lgNegPop )
931  {
932  fprintf( ioQQQ, " PROBLEM Negative population found for abundance of ionization stage of element %4.4s, ZONE=%4ld\n",
934 
935  fprintf( ioQQQ, " Populations were" );
936  for( ion=1; ion <= dense.IonHigh[nelem]+1; ion++ )
937  {
938  fprintf( ioQQQ, "%9.1e", dense.xIonDense[nelem][ion-1] );
939  }
940  fprintf( ioQQQ, "\n" );
941 
942  fprintf( ioQQQ, " destroy vector =" );
943  for( ion=1; ion <= dense.IonHigh[nelem]; ion++ )
944  {
945  fprintf( ioQQQ, "%9.1e", ionbal.RateIonizTot[nelem][ion-1] );
946  }
947  fprintf( ioQQQ, "\n" );
948 
949  /* print some extra stuff if destroy was negative */
950  if( lgNegPop )
951  {
952  fprintf( ioQQQ, " CTHeavy vector =" );
953  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
954  {
955  fprintf( ioQQQ, "%9.1e", atmdat.HeCharExcIonOf[nelem][ion] );
956  }
957  fprintf( ioQQQ, "\n" );
958 
959  fprintf( ioQQQ, " HCharExcIonOf vtr=" );
960  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
961  {
962  fprintf( ioQQQ, "%9.1e", atmdat.HCharExcIonOf[nelem][ion] );
963  }
964  fprintf( ioQQQ, "\n" );
965 
966  fprintf( ioQQQ, " CollidRate vtr=" );
967  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
968  {
969  fprintf( ioQQQ, "%9.1e", ionbal.CollIonRate_Ground[nelem][ion][0] );
970  }
971  fprintf( ioQQQ, "\n" );
972 
973  /* photo rates per subshell */
974  fprintf( ioQQQ, " photo rates per subshell, ion\n" );
975  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
976  {
977  fprintf( ioQQQ, "%3ld", ion );
978  for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
979  {
980  fprintf( ioQQQ, "%9.1e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
981  }
982  fprintf( ioQQQ, "\n" );
983  }
984  }
985 
986  /* now check out creation vector */
987  fprintf( ioQQQ, " create vector =" );
988  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
989  {
990  fprintf( ioQQQ, "%9.1e", ionbal.RateRecomTot[nelem][ion] );
991  }
992  fprintf( ioQQQ, "\n" );
993 
994  ContNegative();
995  ShowMe();
996  cdEXIT(EXIT_FAILURE);
997  }
998 
999  /* option to print ionization and recombination arrays
1000  * prt flag set with print array print arrays command */
1001  if( prt.lgPrtArry[nelem] || lgPrintIt )
1002  {
1003  /* say who we are, what we are doing .... */
1004  fprintf( ioQQQ,
1005  "\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",
1006  elementnames.chElementSym[nelem],
1007  elementnames.chElementName[nelem],
1008  fnzone,
1009  phycon.te ,
1010  dense.eden,
1011  dense.gas_phase[nelem],
1012  abund_total ,
1013  dense.xMolecules[nelem] );
1014  /* total ionization rate, all processes */
1015  fprintf( ioQQQ, " %s Ioniz total " ,elementnames.chElementSym[nelem]);
1016  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1017  {
1018  fprintf( ioQQQ, " %9.2e", ionbal.RateIonizTot[nelem][ion] );
1019  }
1020  fprintf( ioQQQ, "\n" );
1021 
1022  /* sinks from the chemistry network */
1023  fprintf(ioQQQ," %s mole sink ",
1024  elementnames.chElementSym[nelem]);
1025  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1026  {
1027  fprintf(ioQQQ," %9.2e", mole.sink[nelem][ion] );
1028  }
1029  fprintf( ioQQQ, "\n" );
1030 
1031  if( dynamics.lgAdvection )
1032  {
1033  /* sinks from the dynamics */
1034  fprintf(ioQQQ," %s dynm sink ",
1035  elementnames.chElementSym[nelem]);
1036  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1037  {
1038  fprintf(ioQQQ," %9.2e", dynamics.Rate );
1039  }
1040  fprintf( ioQQQ, "\n" );
1041  }
1042 
1043  /* sum of all creation processes */
1044  fprintf( ioQQQ, " %s Recom total " ,elementnames.chElementSym[nelem]);
1045  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1046  {
1047  fprintf( ioQQQ, " %9.2e", ionbal.RateRecomTot[nelem][ion] );
1048  }
1049  fprintf( ioQQQ, "\n" );
1050 
1051  /* sources from the chemistry network */
1052  fprintf(ioQQQ," %s mole source ",
1053  elementnames.chElementSym[nelem]);
1054  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1055  {
1056  fprintf(ioQQQ," %9.2e", mole.source[nelem][ion] );
1057  }
1058  fprintf( ioQQQ, "\n" );
1059 
1060  if( dynamics.lgAdvection )
1061  {
1062  /* source from the dynamcs */
1063  fprintf(ioQQQ," %s dynm sour ",
1064  elementnames.chElementSym[nelem]);
1065  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1066  {
1067  fprintf(ioQQQ," %9.2e", dynamics.Source[nelem][ion] );
1068  }
1069  fprintf( ioQQQ, "\n" );
1070  }
1071 
1072  /* collisional ionization */
1073  fprintf( ioQQQ, " %s Coll ioniz " ,elementnames.chElementSym[nelem] );
1074  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1075  {
1076  fprintf( ioQQQ, " %9.2e", ionbal.CollIonRate_Ground[nelem][ion][0] );
1077  }
1078  fprintf( ioQQQ, "\n" );
1079 
1080  /* UTA ionization */
1081  fprintf( ioQQQ, " %s UTA ioniz " ,elementnames.chElementSym[nelem] );
1082  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1083  {
1084  fprintf( ioQQQ, " %9.2e", ionbal.UTA_ionize_rate[nelem][ion] );
1085  }
1086  fprintf( ioQQQ, "\n" );
1087 
1088  /* photo ionization */
1089  fprintf( ioQQQ, " %s Phot ioniz " ,elementnames.chElementSym[nelem]);
1090  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1091  {
1092  fprintf( ioQQQ, " %9.2e",
1093  ionbal.PhotoRate_Shell[nelem][ion][Heavy.nsShells[nelem][ion]-1][0] );
1094  }
1095  fprintf( ioQQQ, "\n" );
1096 
1097  /* auger ionization */
1098  fprintf( ioQQQ, " %s Auger ioniz " ,elementnames.chElementSym[nelem]);
1099  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1100  {
1101  fprintf( ioQQQ, " %9.2e",
1102  auger[ion] );
1103  }
1104  fprintf( ioQQQ, "\n" );
1105 
1106  /* secondary ionization */
1107  fprintf( ioQQQ, " %s Secon ioniz " ,elementnames.chElementSym[nelem]);
1108  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1109  {
1110  fprintf( ioQQQ, " %9.2e",
1111  secondaries.csupra[nelem][ion] );
1112  }
1113  fprintf( ioQQQ, "\n" );
1114 
1115  /* grain ionization - not total rate but should be dominant process */
1116  fprintf( ioQQQ, " %s ion on grn " ,elementnames.chElementSym[nelem]);
1117  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1118  {
1119  fprintf( ioQQQ, " %9.2e", gv.GrainChTrRate[nelem][ion][ion+1] );
1120  }
1121  fprintf( ioQQQ, "\n" );
1122 
1123  /* chemistry ionization - not total rate but should be dominant process
1124  * not source or sink but just increase ionization within ladder */
1125  fprintf( ioQQQ, " %s ion in chem " ,elementnames.chElementSym[nelem]);
1126  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1127  {
1128  fprintf( ioQQQ, " %9.2e", mole.xMoleChTrRate[nelem][ion][ion+1] );
1129  }
1130  fprintf( ioQQQ, "\n" );
1131 
1132  /* charge exchange ionization */
1133  fprintf( ioQQQ, " %s chr trn ion " ,elementnames.chElementSym[nelem] );
1134  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1135  {
1136  /* sum has units s-1 */
1137  double sum = atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[ipHELIUM][1]+
1138  atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[ipHYDROGEN][1];
1139 
1140  if( nelem==ipHELIUM && ion==0 )
1141  {
1142  sum += atmdat.HeCharExcIonTotal;
1143  }
1144  else if( nelem==ipHYDROGEN && ion==0 )
1145  {
1146  sum += atmdat.HCharExcIonTotal;
1147  }
1148  fprintf( ioQQQ, " %9.2e", sum );
1149  }
1150  fprintf( ioQQQ, "\n" );
1151 
1152  /* radiative recombination */
1153  fprintf( ioQQQ, " %s radiati rec " ,elementnames.chElementSym[nelem]);
1154  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1155  {
1156  fprintf( ioQQQ, " %9.2e", dense.eden*ionbal.RR_rate_coef_used[nelem][ion] );
1157  }
1158  fprintf( ioQQQ, "\n" );
1159 
1160  /* old DR recombination */
1161  fprintf( ioQQQ, " %s dr old rec " ,elementnames.chElementSym[nelem]);
1162  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1163  {
1164  fprintf( ioQQQ, " %9.2e", dense.eden*ionbal.DR_old_rate_coef[nelem][ion] );
1165  }
1166  fprintf( ioQQQ, "\n" );
1167 
1168  /* Badnell DR recombination */
1169  fprintf( ioQQQ, " %s drBadnel rec" ,elementnames.chElementSym[nelem]);
1170  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1171  {
1172  fprintf( ioQQQ, " %9.2e", dense.eden*ionbal.DR_Badnell_rate_coef[nelem][ion] );
1173  }
1174  fprintf( ioQQQ, "\n" );
1175 
1176  /* grain recombination - not total but from next higher ion, should
1177  * be dominant */
1178  fprintf( ioQQQ, " %s rec on grn " ,elementnames.chElementSym[nelem]);
1179  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1180  {
1181  fprintf( ioQQQ, " %9.2e", gv.GrainChTrRate[nelem][ion+1][ion] );
1182  }
1183  fprintf( ioQQQ, "\n" );
1184 
1185  /* grain due to chemistry - not total but from next higher ion, should
1186  * be dominant - not source or sink, just moves up or down */
1187  fprintf( ioQQQ, " %s rec in chem " ,elementnames.chElementSym[nelem]);
1188  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1189  {
1190  fprintf( ioQQQ, " %9.2e", mole.xMoleChTrRate[nelem][ion+1][ion] );
1191  }
1192  fprintf( ioQQQ, "\n" );
1193 
1194  /* charge exchange recombination */
1195  fprintf( ioQQQ, " %s chr trn rec " ,elementnames.chElementSym[nelem]);
1196  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1197  {
1198  double sum =
1199  atmdat.HCharExcRecTo[nelem][ion]*
1201  dense.xIonDense[ipHYDROGEN][1] +
1202 
1203  atmdat.HeCharExcRecTo[nelem][ion]*
1205  dense.xIonDense[ipHELIUM][1];
1206 
1207  if( nelem==ipHELIUM && ion==0 )
1208  {
1209  sum += atmdat.HeCharExcRecTotal;
1210  }
1211  else if( nelem==ipHYDROGEN && ion==0 )
1212  {
1213  sum += atmdat.HCharExcRecTotal;
1214  }
1215  fprintf( ioQQQ, " %9.2e", sum );
1216  }
1217  fprintf( ioQQQ, "\n" );
1218 
1219  /* the "new" abundances the resulted from the previous ratio */
1220  fprintf( ioQQQ, " %s Abun [cm-3] " ,elementnames.chElementSym[nelem] );
1221  for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
1222  {
1223  fprintf( ioQQQ, " %9.2e", dense.xIonDense[nelem][ion] );
1224  }
1225  fprintf( ioQQQ, "\n" );
1226  }
1227 
1228  if( PARALLEL_MODE )
1229  {
1230  free(ipiv);
1231  free(source);
1232  free(sink);
1233  free(xmat);
1234  free(xmatsave);
1235  free(auger);
1236  }
1237  return;
1238 }
1239 
1240 /*
1241 
1242  Solve an ionization level system with specified ionization and
1243  recombination rates between neighboring ions, and additional sink
1244  and source terms. The sink array is overwritten, and the results
1245  appear in the source array.
1246 
1247  Written in matrix form, the algorithm is equivalent to the
1248  tridiagonal algorithm in Numerical Recipes applied to:
1249 
1250  / i_0+a_0 -r_0 . . . \ / x_0 \ / s_0 \
1251  | -i_0 i_1+a_1+r_0 -r_1 . . | | x_1 | | s_1 |
1252  | . -i_1 i_2+a_2+r_1 -r_2 . | | x_2 | | s_2 |
1253  | . . (etc....) | | ... | = | ... |
1254  \ . . . / \ / \ /
1255 
1256  where i, r are the ionization and recombination rates, s is the
1257  source rate and a is the sink rate.
1258 
1259  This matrix is diagonally dominant only when the sink terms are
1260  large -- the alternative method coded here prevents rounding error
1261  in the diagonal terms disturbing the solution when this is not the
1262  case.
1263 
1264 */
1265 
1266 /* solveions tridiagonal solver but optimized for structure of balance matrix */
1267 void solveions(double *ion, double *rec, double *snk, double *src,
1268  long int nlev, long int nmax)
1269 {
1270  double kap, bet;
1271  long int i;
1272 
1273  DEBUG_ENTRY( "solveions()" );
1274 
1275  if(nmax != -1)
1276  {
1277  /* Singular case */
1278  src[nmax] = 1.;
1279  for(i=nmax;i<nlev-1;i++)
1280  src[i+1] = src[i]*ion[i]/rec[i];
1281  for(i=nmax-1;i>=0;i--)
1282  src[i] = src[i+1]*rec[i]/ion[i];
1283  }
1284  else
1285  {
1286  i = 0;
1287  kap = snk[0];
1288  for(i=0;i<nlev-1;i++)
1289  {
1290  bet = ion[i]+kap;
1291  if(bet == 0.)
1292  {
1293  fprintf(ioQQQ,"Ionization solver error\n");
1294  cdEXIT(EXIT_FAILURE);
1295  }
1296  bet = 1./bet;
1297  src[i] *= bet;
1298  src[i+1] += ion[i]*src[i];
1299  snk[i] = bet*rec[i];
1300  kap = kap*snk[i]+snk[i+1];
1301  }
1302  bet = kap;
1303  if(bet == 0.)
1304  {
1305  fprintf(ioQQQ,"Ionization solver error\n");
1306  cdEXIT(EXIT_FAILURE);
1307  }
1308  src[i] /= bet;
1309 
1310  for(i=nlev-2;i>=0;i--)
1311  {
1312  src[i] += snk[i]*src[i+1];
1313  }
1314  }
1315 }

Generated for cloudy by doxygen 1.8.4