cloudy  trunk
 All Data Structures Namespaces 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-2022 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 "newton_step.h"
12 #include "conv.h"
13 #include "secondaries.h"
14 #include "phycon.h"
15 #include "atmdat.h"
16 #include "heavy.h"
17 #include "elementnames.h"
18 #include "ionbal.h"
19 #include "trace.h"
20 #include "ion_trim.h"
21 #include "mole.h"
22 #include "freebound.h"
23 #include "dense.h"
24 #include "thirdparty.h"
25 
26 STATIC bool lgTrivialSolution( long nelem, double abund_total );
27 
28 STATIC void find_solution( long nelem, long ion_range, valarray<double> &xmat, valarray<double> &source);
29 
30 STATIC void fill_array( long int nelem, long ion_range, valarray<double> &xmat, valarray<double> &auger );
31 
32 STATIC void fill_ext_src_and_snk( long nelem, long ion_low, long ion_range, valarray<double> &xmat, valarray<double> &source );
33 
34 STATIC double get_total_abundance_ions( long int nelem );
35 
36 STATIC void HomogeneousSource( long nelem, long ion_low, long ion_range, valarray<double> &xmat, valarray<double> &source, double abund_total );
37 
38 STATIC void renorm_solution( long nelem, long ion_range, long ion_low, double *source, double abund_total, bool *lgNegPop );
39 
40 STATIC void store_new_densities( long nelem, long ion_range, long ion_low, double *source );
41 
42 STATIC void clean_up( long nelem, double abund_total );
43 
44 STATIC void PrintRates( long nelem, bool lgNegPop, double abund_total, valarray<double> &auger, bool lgPrintIt );
45 
46 void solveions(double *ion, double *rec, double *snk, double *src,
47  long int nlev, long int nmax);
48 
49  /* this will be used to address the 2d arrays */
50 # undef MAT
51 # define MAT(M_,I_,J_) ((M_)[(I_)*(ion_range)+(J_)])
52 
53 # undef MAT1
54 # define MAT1(M_,I_,J_) ((M_)[(I_)*(ion_range1)+(J_)])
55 
56 # undef MAT2
57 # define MAT2(M_,I_,J_) ((M_)[(I_)*(ion_range2)+(J_)])
58 
59 void ion_solver( long int nelem, bool lgPrintIt)
60 {
61  double abund_total = 0.0;
62  long ion_low = dense.IonLow[nelem];
63  bool lgNegPop = false;
64  double error = 0.0;
65 
66  DEBUG_ENTRY( "ion_solver()" );
67 
69 
70  long ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
71  valarray<double> xmat(ion_range*ion_range);
72  valarray<double> source(ion_range);
73  valarray<double> auger(LIMELM+1);
74 
75  // V/W cycle would do ion solution *first*, then iso, then
76  // ion again if iso made any changes
77  // Seem to need to break out after iso at the moment. Why?
78 
79  static ConvergenceCounter cctr=conv.register_("ION_SOLVES");
80  ++cctr;
81  for (long it=0; it<4; ++it)
82  {
83  static ConvergenceCounter cctrl=conv.register_("ISO_LOOPS");
84  ++cctrl;
85 
86  for( long ipISO=ipH_LIKE; ipISO<MIN2(NISO,nelem+1); ipISO++ )
87  {
88  if( (dense.IonHigh[nelem] >= nelem - ipISO) &&
89  (dense.IonLow[nelem] <= nelem - ipISO))
90  {
91  iso_set_ion_rates(ipISO, nelem);
92  }
93  }
94 
95  //if ( it > 1 && error > 0.0)
96  // fprintf(ioQQQ,"ion nelem %ld loop %ld error %g nzone %ld\n",nelem,it,error,nzone);
97 
98  abund_total = get_total_abundance_ions( nelem );
99 
100 
101  {
102  /* this sets up a fake ionization balance problem, with a trivial solution,
103  * for debugging the ionization solver */
104  enum {DEBUG_LOC=false};
105  if( DEBUG_LOC && nelem==ipCARBON )
106  {
107  broken();/* within optional debug print statement */
108  dense.IonLow[nelem] = 0;
109  dense.IonHigh[nelem] = 3;
110  abund_total = 1.;
111  long limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);
112  /* make up ionization and recombination rates */
113  for( long ion=dense.IonLow[nelem]; ion <= limit; ion++ )
114  {
115  double fac=1;
116  if(ion)
117  fac = 1e-10;
118  ionbal.RateRecomTot[nelem][ion] = 100.;
119  for( long ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
120  {
121  /* direct photoionization of this shell */
122  ionbal.PhotoRate_Shell[nelem][ion][ns][0] = fac;
123  }
124  }
125  }
126  }
127 
128  bool lgTrivial = lgTrivialSolution( nelem, abund_total );
129  if( !lgTrivial )
130  {
131  // fill xmat and source with appropriate terms
132  fill_array( nelem, ion_range, xmat, auger );
133 
134  // Fill source and sink terms from chemistry.
135  fill_ext_src_and_snk( nelem, ion_low, ion_range, xmat, source );
136 
137  // decide if matrix is homogeneous
138  HomogeneousSource( nelem, ion_low, ion_range, xmat, source, abund_total );
139 
140  // Now find the solution
141  find_solution( nelem, ion_range, xmat, source);
142 
143  renorm_solution( nelem, ion_range, ion_low, &source[0], abund_total, &lgNegPop );
144 
145  // save the results in the global density variables
146  store_new_densities( nelem, ion_range, ion_low, &source[0] );
147 
148  clean_up( nelem, abund_total );
149 
150  ASSERT(ion_range >= dense.IonHigh[nelem]-dense.IonLow[nelem]+1);
151  if (ion_range != dense.IonHigh[nelem]-dense.IonLow[nelem]+1)
152  {
153  ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
154  xmat.resize(ion_range*ion_range);
155  source.resize(ion_range);
156  }
157  }
158 
159  error = 0.0;
160  /* update and solve iso levels */
161  for( long ipISO=ipH_LIKE; ipISO<MIN2(NISO,nelem+1); ipISO++ )
162  {
163  double thiserror;
164  iso_solve( ipISO, nelem, thiserror );
165  if (thiserror > error)
166  error = thiserror;
167  }
168 
169  if ( error < 1e-4 )
170  break;
171  }
172 
173  iso_satellite_update(nelem);
174 
175  for( long ipISO=ipH_LIKE; ipISO<MIN2(NISO,nelem+1); ipISO++ )
176  {
177  /* now evaluate departure coefficients */
178  iso_departure_coefficients( ipISO, nelem );
179  }
180 
181  if( prt.lgPrtArry[nelem] || lgPrintIt )
182  PrintRates( nelem, lgNegPop, abund_total, auger, lgPrintIt );
183 
184  return;
185 }
186 
187 STATIC bool lgTrivialSolution( long nelem, double abund_total )
188 {
189  double renorm;
190  /* return if IonHigh is zero, since no ionization at all */
191  if( dense.IonHigh[nelem] == dense.IonLow[nelem] )
192  {
193  /* set the atom to the total gas phase abundance */
194  dense.xIonDense[nelem][dense.IonHigh[nelem]] = abund_total;
195  for ( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
196  iso_renorm(nelem, ipISO, renorm);
197  return true;
198  }
199  else if( dense.lgSetIoniz[nelem] )
200  {
201  /* option to force ionization distribution with element name ioniz */
202  for( long ion=0; ion<nelem+2; ++ion )
203  dense.xIonDense[nelem][ion] = dense.SetIoniz[nelem][ion]*abund_total;
204  for ( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
205  iso_renorm(nelem, ipISO, renorm);
206  return true;
207  }
208  else
209  return false;
210 }
211 
212 STATIC void find_solution( long nelem, long ion_range, valarray<double> &xmat, valarray<double> &source)
213 {
214  int32 nerror;
215  valarray<double> xmatsave(ion_range*ion_range);
216  valarray<double> sourcesave(ion_range);
217  valarray<int32> ipiv(ion_range);
218 
219  DEBUG_ENTRY( "find_solution()" );
220 
221  // save copy of xmat before continuing.
222  for( long i=0; i< ion_range; ++i )
223  {
224  sourcesave[i] = source[i];
225  for( long j=0; j< ion_range; ++j )
226  {
227  MAT( xmatsave, i, j ) = MAT( xmat, i, j );
228  }
229  }
230 
231  if (1)
232  {
233  nerror = solve_system(xmat,source,ion_range,NULL);
234 
235  if (nerror)
236  {
237  // solve_system doesn't return a valid solution, so just
238  // provide initial values to keep things running
239  fprintf(ioQQQ,"Error for nelem %ld, active ion range %ld--%ld\n",
240  nelem,dense.IonLow[nelem],dense.IonHigh[nelem]);
241  fprintf(ioQQQ,"Initial ion abundances:");
242  for( long j=0; j<nelem+2; ++j )
243  fprintf(ioQQQ," %g",dense.xIonDense[nelem][j]);
244  fprintf(ioQQQ,"\n");
245  for( long j=0; j<ion_range; ++j )
246  source[j] = dense.xIonDense[nelem][dense.IonLow[nelem]+j];
247  }
248  }
249  else
250  {
251  /* this is the default solver - now get new solution */
252  nerror = 0;
253  /* Use general matrix solver */
254  getrf_wrapper(ion_range, ion_range, &xmat[0], ion_range, &ipiv[0], &nerror);
255  if( nerror != 0 )
256  {
257  fprintf( ioQQQ,
258  " DISASTER ion_solver: dgetrf finds singular or ill-conditioned matrix nelem=%li %s ion_range=%li, nConv %li IonLow %li IonHi %li\n",
259  nelem ,
260  elementnames.chElementSym[nelem],
261  ion_range,
262  conv.nTotalIoniz ,
263  dense.IonLow[nelem], dense.IonHigh[nelem]);
264  fprintf( ioQQQ, " xmat follows\n");
265  for( long i=0; i<ion_range; ++i )
266  {
267  for( long j=0;j<ion_range;j++ )
268  {
269  fprintf(ioQQQ,"%e\t",MAT(xmatsave,j,i));
270  }
271  fprintf(ioQQQ,"\n");
272  }
273  fprintf(ioQQQ,"source follows\n");
274  for( long i=0; i<ion_range;i++ )
275  {
276  fprintf(ioQQQ,"%e\t",sourcesave[i]);
277  }
278  fprintf(ioQQQ,"\n");
280  }
281  getrs_wrapper('N', ion_range, 1, &xmat[0], ion_range, &ipiv[0], &source[0], ion_range, &nerror);
282  if( nerror != 0 )
283  {
284  fprintf( ioQQQ, " DISASTER ion_solver: dgetrs finds singular or ill-conditioned matrix nelem=%li ionrange=%li\n",
285  nelem , ion_range );
287  }
288  }
289 
290  {
291  /* this is to debug following failed assert */
292  enum {DEBUG_LOC=false};
293  if( DEBUG_LOC && nelem == ipHYDROGEN )
294  {
295  fprintf(ioQQQ,"debuggg ion_solver1 %ld\t%.2f\t%.4e\t%.4e\tIon\t%.3e\tRec\t%.3e\n",
296  nelem,
297  fnzone,
298  phycon.te,
299  dense.eden,
300  ionbal.RateIonizTot(nelem,0) ,
301  ionbal.RateRecomTot[nelem][0]);
302  fprintf(ioQQQ," Msrc %.3e %.3e\n", mole.source[nelem][0], mole.source[nelem][1]);
303  fprintf(ioQQQ," Msnk %.3e %.3e\n", mole.sink[nelem][0], mole.sink[nelem][1]);
304  fprintf(ioQQQ," Poprat %.3e nomol %.3e\n",source[1]/source[0],
305  ionbal.RateIonizTot(nelem,0)/ionbal.RateRecomTot[nelem][0]);
306  }
307  }
308 
309  for( long i=0; i<ion_range; i++ )
310  {
311  ASSERT( !isnan( source[i] ) );
312  ASSERT( source[i] < MAX_DENSITY );
313  }
314 
315  return;
316 }
317 
318 STATIC void renorm_solution( long nelem, long ion_range, long ion_low, double *source, double abund_total, bool *lgNegPop )
319 {
320  DEBUG_ENTRY( "renorm_solution()" );
321 
322  ASSERT( nelem >= 0 );
323  ASSERT( nelem < LIMELM );
324  ASSERT( ion_range <= nelem + 2 );
325  ASSERT( ion_low >= 0 );
326  ASSERT( ion_low <= nelem + 1 );
327 
328 #if 0
329 # define RJRW 0
330  if( RJRW && 0 )
331  {
332  /* verify that the rates are sensible */
333  double test;
334  for(long i=0; i<ion_range; i++) {
335  test = 0.;
336  for(long j=0; j<ion_range; j++) {
337  test = test+source[j]*MAT(xmatsave,j,i);
338  }
339  fprintf(ioQQQ,"%e\t",test);
340  }
341  fprintf(ioQQQ,"\n");
342 
343  test = 0.;
344  fprintf(ioQQQ," ion %li abundance %.3e\n",nelem,abund_total);
345  for( long ion=dense.IonLow[nelem]; ion < dense.IonHigh[nelem]; ion++ )
346  {
347  if( ionbal.RateRecomTot[nelem][ion] != 0 && source[ion-ion_low] != 0 )
348  fprintf(ioQQQ," %li %.3e %.3e : %.3e\n",ion,source[ion-ion_low],
349  source[ion-ion_low+1]/source[ion-ion_low],
350  ionbal.RateIonizTot(nelem,ion)/ionbal.RateRecomTot[nelem][ion]);
351  else
352  fprintf(ioQQQ," %li %.3e [One ratio infinity]\n",ion,source[ion-ion_low]);
353  test += source[ion-ion_low];
354  }
355  }
356 #endif
357 
358  /*
359  * >> chng 03 jan 15 rjrw:- terms are now included for
360  * molecular sources and sinks.
361  *
362  * When the network is not in equilibrium, this will lead to a
363  * change in the derived abundance of coupled ions after the matrix
364  * solution.
365  *
366  * We therefore renormalize to keep the total abundance of the
367  * states treated by the ionization ladder constant -- only the
368  * molecular network is allowed to change this.
369  *
370  * The difference between `renorm' and 1. is a measure of the
371  * quality of the solution (it will be 1. if the rate of transfer
372  * into the ionization ladder species balances the rate of transfer
373  * out, for the consistent relative abundances)
374  *
375  */
376 
377  /* source[i] contains new solution for ionization populations
378  * save resulting abundances into main ionization density array,
379  * while checking whether any negative level populations occurred */
380  *lgNegPop = false;
381  for( long i=0; i < ion_range; i++ )
382  {
383  long ion = i+ion_low;
384 
385  if( source[i] < 0. )
386  {
387  /* >>chng 04 dec 04, put test on neg abund here, don't print unless value is very -ve */
388  /* >>chng 06 feb 28, from -1e-10 to -1e-9, sim func_t10 had several negative
389  * during initial search, due to extremely high ionization */
390  /* >>chng 06 mar 11, from 1e-9 to 2e-9 make many struc elements floats from double */
391  if( source[i]<-2e-9 )
392  {
393  fprintf(ioQQQ,
394  " PROBLEM negative ion population in ion_solver, nelem=%li, %s ion=%li val=%.3e Search?%c zone=%li iteration=%li\n",
395  nelem ,
396  elementnames.chElementSym[nelem],
397  ion ,
398  source[i] ,
399  TorF(conv.lgSearch) ,
400  nzone ,
401  iteration );
402  *lgNegPop = true;
403  fixit("break PrintRates into one NegPop case and one trace? No auger defined here.");
404  //PrintRates( nelem, *lgNegPop, abund_total, auger );
405  }
406 
407  fixit("Kill this bit and force exit on negative populations.");
408 #if 1
409  source[i] = 0.;
410  /* if this is one of the iso seq model atoms then must also zero out pops */
411  //if( ion == nelem+1-NISO ) //newmole had this should have been next line?
412  if( ion > nelem-NISO && ion < nelem + 1 )
413  {
414  long int ipISO = nelem - ion;
415  ASSERT( ipISO>=0 && ipISO<NISO );
416  for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
417  iso_sp[ipISO][nelem].st[level].Pop() = 0.;
418  }
419 #endif
420  }
421  }
422 
423  double renormnew = 1.;
424  {
425  double dennew = 0.;
426 
427  /* find total population to renorm - also here check that negative pops did not occur */
428  for( long i=0;i < ion_range; i++ )
429  {
430  dennew += source[i];
431  }
432 
433  if(dennew > 0.)
434  {
435  renormnew = abund_total / dennew;
440  }
441  else
442  {
443  renormnew = 1.;
444  }
445 
446  for( long i=0;i < ion_range; i++ )
447  {
448  source[i] *= renormnew;
449  if( source[i] >= MAX_DENSITY )
450  {
451  long ion = i+ion_low;
452  fprintf( ioQQQ, "PROBLEM DISASTER Huge density in ion_solver, nelem %ld ion %ld source %e renormnew %e\n",
453  nelem, ion, source[i], renormnew );
454  }
455  }
456  }
457  /* check not negative, should be +ve, can be zero if species has become totally molecular.
458  * this happens for hydrogen if no cosmic rays, or cr ion set very low */
459  if( renormnew < 0)
460  {
461  fprintf(ioQQQ,"PROBLEM impossible value of renorm \n");
462  }
463  ASSERT( renormnew>=0 );
464 
465  return;
466 }
467 
468 STATIC void store_new_densities( long nelem, long ion_range, long ion_low, double *source )
469 {
470  DEBUG_ENTRY( "store_new_densities()" );
471 
472  for( long i=0; i < ion_range; i++ )
473  {
474  long ion = i+ion_low;
475  dense.xIonDense[nelem][ion] = source[i];
476  ASSERT( dense.xIonDense[nelem][ion] < MAX_DENSITY );
477  }
478 
479  return;
480 }
481 
482 STATIC void clean_up( long nelem, double abund_total )
483 {
484  DEBUG_ENTRY( "clean_up()" );
485 
486  fixit("this should only be done if trimming is not disabled?");
487 
488  /* Zero levels with abundances < 1e-25 which which will suffer numerical noise */
489  ion_trim_small(nelem, abund_total);
490 
491  for ( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
492  {
493  double renorm;
494  iso_renorm(nelem, ipISO, renorm);
495  }
496 
497  // sanity check, either offset stages of low and high ionization
498  ASSERT( (dense.IonLow[nelem] <= dense.IonHigh[nelem]) ||
499  // both totally neutral
500  (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) ||
501  // both fully stripped
502  ( dense.IonLow[nelem]==nelem+1 && dense.IonHigh[nelem]==nelem+1 ) );
503 
504  return;
505 }
506 
507 STATIC double get_total_abundance_ions( long int nelem )
508 {
509  DEBUG_ENTRY( "get_total_abundance_ions()" );
510 
511  ASSERT( nelem >= 0 );
512  ASSERT( nelem < LIMELM );
513 
514  ionbal.elecsnk[nelem] = 0.;
515  ionbal.elecsrc[nelem] = 0.;
516 
517  double abund_total = 0.;
518  for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
519  {
520  abund_total += dense.xIonDense[nelem][ion];
521  }
522 
523  realnum tot1 = dense.gas_phase[nelem];
524  realnum tot2 = (realnum)(dense.xMolecules(nelem)+abund_total);
525  if (0)
526  {
527  if (fabs(tot2-tot1) > conv.GasPhaseAbundErrorAllowed*tot1)
528  fprintf(ioQQQ,"%ld %13.8g %13.8g %13.8g %13.8g\n",nelem,tot1,tot2,abund_total,tot2/tot1 - 1.0);
529  }
530 
531  ASSERT( fp_equal_tol(tot1, tot2, realnum(conv.GasPhaseAbundErrorAllowed*tot1 + 100.f*FLT_MIN)) );
532 
533  ASSERT( abund_total < MAX_DENSITY );
534 
535  return abund_total;
536 }
537 
538 STATIC void fill_array( long int nelem, long ion_range, valarray<double> &xmat, valarray<double> &auger )
539 {
540  long int limit,
541  IonProduced;
542  double rateone;
543  long ion_low;
544 
545  valarray<double> sink(ion_range);
546  valarray<int32> ipiv(ion_range);
547 
548  DEBUG_ENTRY( "fill_array()" );
549 
550  /* this is on the c scale, so H is 0 */
551  ASSERT( nelem >= 0);
552  ASSERT( dense.IonLow[nelem] >= 0 );
553  ASSERT( dense.IonHigh[nelem] >= 0 );
554 
555  /* impossible for HIonFrac[nelem] to be zero if IonHigh(nelem)=nelem+1
556  * HIonFrac(nelem) is stripped to hydrogen */
557  /* >>chng 01 oct 30, to assert */
558  //ASSERT( (dense.IonHigh[nelem] < nelem + 1) || dense.xIonDense[nelem][nelem+1-ipH_LIKE] > 0. );
559 
560  /* zero out the ionization and recombination rates that we will modify here,
561  * but not the iso-electronic stages which are done elsewhere,
562  * the nelem stage of ionization is he-like,
563  * the nelem+1 stage of ionization is h-like */
564 
565  /* loop over stages of ionization that we solve for here,
566  * up through and including one less than nelem-NISO,
567  * never actually do highest NISO stages of ionization since they
568  * come from the ionization ratio from the next lower stage */
569  limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);
570 
571  /* the full range of ionization - this is number of ionization stages */
572  ASSERT( ion_range <= nelem+2 );
573 
574  ion_low = dense.IonLow[nelem];
575 
576  /* zero-out loop comes before main loop since there are off-diagonal
577  * elements in the main ionization loop, due to multi-electron processes,
578  * TotIonizRate and TotRecom were already set in h-like and he-like solvers
579  * other recombination rates were already set by routines responsible for them */
580  for( long ion_from=0; ion_from <= limit; ion_from++ )
581  {
582  for( long ion_to=0; ion_to < nelem+2; ion_to++ )
583  {
584  ionbal.RateIoniz[nelem][ion_from][ion_to] = 0.;
585  }
586  }
587 
588  /* auger is used only for debug printout - it is special because with multi-electron
589  * Auger ejection, very high stages of ionization can be produced, well beyond
590  * the highest stage considered here. so we malloc to the limit */
591  for( long i=0; i< LIMELM+1; ++i )
592  {
593  auger[i] = 0.;
594  }
595 
596  /* zero out xmat */
597  for( long i=0; i< ion_range; ++i )
598  {
599  for( long j=0; j< ion_range; ++j )
600  {
601  MAT( xmat, i, j ) = 0.;
602  }
603  }
604 
605  bool lgGrainsOn = gv.lgDustOn() && ionbal.lgGrainIonRecom && gv.lgGrainPhysicsOn;
606 
607  /* Now put in all recombination and ionization terms from CO_mole() that
608  * come from molecular reactions. this traces molecular process that
609  * change ionization stages with this ladder - but do not remove from
610  * the ladder */
611  for( long ion_to=dense.IonLow[nelem]; ion_to <= dense.IonHigh[nelem]; ion_to++ )
612  {
613  for( long ion_from=dense.IonLow[nelem]; ion_from <= dense.IonHigh[nelem]; ++ion_from )
614  {
615  /* do not do ion onto itself */
616  if( ion_to != ion_from )
617  {
618  /* CT with molecules */
619  rateone = mole.xMoleChTrRate[nelem][ion_from][ion_to] * atmdat.lgCTOn;
620  MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
621  MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
622  /* CT with grains */
623  rateone = gv.GrainChTrRate[nelem][ion_from][ion_to]*lgGrainsOn;
624  MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
625  MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
626  }
627  }
628  }
629 
630  for( long ion=dense.IonLow[nelem]; ion <= limit; ion++ )
631  {
632  /* thermal & secondary collisional ionization */
633  ionbal.RateIoniz[nelem][ion][ion+1] +=
634  ionbal.CollIonRate_Ground[nelem][ion][0] +
635  secondaries.csupra[nelem][ion] +
636  /* inner shell ionization by UTA lines */
637  ionbal.UTA_ionize_rate[nelem][ion];
638 
639  /* loop over all atomic sub-shells to include photoionization */
640  for( long ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
641  {
642  /* this is the primary ionization rate - add to diagonal element,
643  * test on ion stage is so that we don't include ionization from the very highest
644  * ionization stage to even higher - since those even higher stages are not considered
645  * this would appear as a sink - but populations of this highest level is ensured to
646  * be nearly trivial and neglecting it production of even higher ionization OK */
647  /* >>chng 04 nov 29 RJRW, include following in this branch so only
648  * evaluated when below ions done with iso-sequence */
649  if( ion+1-ion_low < ion_range )
650  {
651  /* this will be redistributed into charge states in following loop */
652 
653  /* t_yield::Inst().nelec_eject(nelem,ion,ns) is total number of electrons that can
654  * possibly be freed
655  * loop over nej, the number of electrons ejected including the primary,
656  * nej = 1 is primary, nej > 1 includes primary plus Auger
657  * t_yield::Inst().elec_eject_frac is prob of nej electrons */
658  for( long nej=1; nej <= t_yield::Inst().nelec_eject(nelem,ion,ns); nej++ )
659  {
660  /* this is the ion that is produced by this ejection,
661  * limited by highest possible stage of ionization -
662  * do not want to ignore ionization that go beyond this */
663  IonProduced = MIN2(ion+nej,dense.IonHigh[nelem]);
664  rateone = ionbal.PhotoRate_Shell[nelem][ion][ns][0]*
665  t_yield::Inst().elec_eject_frac(nelem,ion,ns,nej-1);
666 
667  /* direct photoionization of this shell */
668  ionbal.RateIoniz[nelem][ion][IonProduced] += rateone;
669 
670  /* only used for possible printout - multiple electron Auger rate -
671  * do not count one-electron as Auger */
672  if( nej>1 )
673  auger[IonProduced-1] += rateone;
674  }
675  }
676  }
677  }
678 
679  for( long ion=dense.IonLow[nelem]; ion < dense.IonHigh[nelem]; ion++ )
680  {
681  /* this is charge transfer recombination to H and He, ionization of other species */
682  double ction;
683  long ipISO=nelem-ion;
684 
685  if ( nelem < t_atmdat::NCX && nelem == ipISO )
686  {
687  ction = atmdat.CharExcIonTotal[nelem] * iso_sp[ipISO][nelem].st[0].Pop() / SDIV(dense.xIonDense[nelem][nelem-ipISO]);
688  }
689  else
690  {
691  ction=0;
692  for (long nelem1=0; nelem1 < t_atmdat::NCX; ++nelem1)
693  ction += atmdat.CharExcIonOf[nelem1][nelem][ion]*dense.xIonDense[nelem1][1];
694  }
695  /* depopulation processes enter with negative sign */
696  MAT( xmat, ion-ion_low, ion-ion_low ) -= ction;
697  MAT( xmat, ion-ion_low, ion+1-ion_low ) += ction;
698  }
699 
700  for( long ion=dense.IonLow[nelem]; ion < dense.IonHigh[nelem]; ion++ )
701  {
702  /* this is charge transfer ionization of H and He, recombination of other species */
703  double ctrec;
704  long ipISO=nelem-ion;
705 
706  // NCX number of ions with CX rates, so this is test on hydrogenic H, He
707  if ( nelem < t_atmdat::NCX && nelem == ipISO )
708  {
709  ctrec = atmdat.CharExcRecTotal[nelem];
710  }
711  else
712  {
713  // the remainder, so a loop over He^0 CX
714  ctrec = 0.;
715  for (long nelem1=0; nelem1<t_atmdat::NCX; ++nelem1)
716  {
717  long ipISO = nelem1;
718  ctrec +=
719  atmdat.CharExcRecTo[nelem1][nelem][ion]*iso_sp[ipISO][nelem1].st[0].Pop();
720  }
721  }
722 
723  MAT( xmat, ion+1-ion_low, ion+1-ion_low ) -= ctrec;
724  MAT( xmat, ion+1-ion_low, ion-ion_low ) += ctrec;
725  }
726 
727  for( long ion_from=dense.IonLow[nelem]; ion_from < dense.IonHigh[nelem]; ion_from++ )
728  {
729  for( long ion_to=ion_from+1; ion_to <= dense.IonHigh[nelem]; ion_to++ )
730  {
731  ionbal.elecsrc[nelem] += ionbal.RateIoniz[nelem][ion_from][ion_to]*dense.xIonDense[nelem][ion_from]*
732  (ion_to-ion_from);
733  /* depopulation processes enter with negative sign */
734  MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= ionbal.RateIoniz[nelem][ion_from][ion_to];
735  MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += ionbal.RateIoniz[nelem][ion_from][ion_to];
736  }
737  }
738 
739  for( long ion=dense.IonLow[nelem]; ion<dense.IonHigh[nelem]; ion++ )
740  {
741  /* loss of next higher ion due to recombination to this ion stage */
742  MAT( xmat, ion+1-ion_low, ion+1-ion_low ) -= ionbal.RateRecomTot[nelem][ion];
743  MAT( xmat, ion+1-ion_low, ion-ion_low ) += ionbal.RateRecomTot[nelem][ion];
744  ionbal.elecsnk[nelem] += ionbal.RateRecomTot[nelem][ion]*dense.xIonDense[nelem][ion+1];
745  }
746 
747  return;
748 }
749 
750 STATIC void fill_ext_src_and_snk( long nelem, long ion_low, long ion_range, valarray<double> &xmat, valarray<double> &source )
751 {
752  DEBUG_ENTRY( "fill_ext_src_and_snk()" );
753 
754  for( long i=0; i<ion_range;i++ )
755  {
756  source[i] = 0.;
757  }
758 
759  for( long i=0; i<ion_range;i++ )
760  {
761  long ion = i+ion_low;
762 
763  /* these are the external source and sink terms */
764  /* need negative sign to get positive pops */
765  source[i] -= mole.source[nelem][ion];
766  MAT( xmat, i, i ) -= mole.sink[nelem][ion];
767  }
768 
769  return;
770 }
771 
772 STATIC void HomogeneousSource( long nelem, long ion_low, long ion_range, valarray<double> &xmat, valarray<double> &source, double abund_total )
773 {
774  bool lgHomogeneous = true;
775  double totsrc,
776  totscl,
777  maxdiag;
778 
779  DEBUG_ENTRY( "lgHomogeneousSource()" );
780 
781  /* this will be sum of source and sink terms, will be used to decide if
782  * matrix is singular */
783  totsrc = totscl = maxdiag = 0.;
784  for( long i=0; i<ion_range;i++ )
785  {
786  long ion = i + ion_low;
787 
788  totsrc -= source[i];
789  fixit("need better way to determine total without molecular sinks");
790  // In old newmole, totscl was calculated before
791  // mole.sink terms are added into xmat, but those terms are now already done.
792  // the above hack takes care of that, but a cleaner solution is needed.
793  double diag = -(MAT( xmat, i, i )+mole.sink[nelem][ion]);
794  if (diag > maxdiag)
795  maxdiag = diag;
796  totscl += diag*dense.xIonDense[nelem][ion];
797  }
798 
799  // Largest condition number before we decide that conservation will get
800  // a better answer than the raw linear solver
801  const double CONDITION_SCALE = 1e8;
802  double conserve_scale = maxdiag/CONDITION_SCALE;
803  ASSERT( conserve_scale > 0.0 );
804 
805  /* matrix is not homogeneous if source is non-zero */
806  if( totsrc > 1e-10*totscl )
807  lgHomogeneous = false;
808 
809  fixit("dynamics rates need to be moved into fill_array?");
810  /* chng 03 jan 13 rjrw, add in dynamics if required here,
811  * - only do advection if we have not overrun the radius scale */
813  && !dynamics.lgEquilibrium && dynamics.Rate != 0. )
814  {
815  for( long i=0; i<ion_range;i++ )
816  {
817  long ion = i+ion_low;
818  MAT( xmat, i, i ) -= dynamics.Rate;
819  source[i] -= dynamics.Source[nelem][ion];
820  /* fprintf(ioQQQ," %li %li %.3e (%.3e %.3e)\n",i,i,MAT(*xmat,i,i),
821  dynamics.Rate, dynamics.Source[nelem][ion]);*/
822  }
823  lgHomogeneous = false;
824  }
825 
826  /* >>chng 06 nov 21, for very high ionization parameter sims the H molecular
827  * fraction can become so small that atom = atom + molecule. In this case we
828  * will not count system as an inhomogeneous case since linear algebra package
829  * will fail. If molecules are very small, count as homogeneous.
830  * Note that we have already added sink terms to the main matrix and they
831  * will not be removed, a possible source of error, but they must not
832  * have been significant, given that the molecular fraction is so small */
833  if( !lgHomogeneous && ion_range==2 )
834  {
835  /* solve algebraically */
836  double a = MAT( xmat, 0, 0 ),
837  b = MAT( xmat, 1, 0 ) ,
838  c = MAT( xmat, 0, 1 ) ,
839  d = MAT( xmat, 1, 1 );
840  b = SDIV(b);
841  d = SDIV(d);
842  double ratio1 = a/b , ratio2 = c/d , fratio1=fabs(a/b),fratio2=fabs(c/d);
843  if( fabs(ratio1-ratio2) <= DBL_EPSILON*1e4*MAX2(fratio1,fratio2) )
844  {
845  lgHomogeneous = true;
846  }
847  }
848 
849 #if 1
850  if(
851  fabs(dense.xMolecules(nelem)) <DBL_EPSILON*100.*dense.gas_phase[nelem] )
852  {
853  lgHomogeneous = true;
854  }
855 #endif
856 
857  /* this is true if no source terms
858  * we will use total population and species conservation to replace one
859  * set of balance equations since overdetermined */
860  if( 1 || lgHomogeneous )
861  {
862  double rate_ioniz=1., rate_recomb /*, scale = 0.*/;
863  /* Simple estimate of most abundant ion */
864  long jmax = 0;
865  for( long i=0; i<ion_range-1;i++)
866  {
867  long ion = i+ion_low;
868  rate_ioniz *= ionbal.RateIonizTot(nelem,ion);
869  rate_recomb = ionbal.RateRecomTot[nelem][ion];
870  /* trips if ion rate zero, so all the gas will be more neutral than this */
871  if( rate_ioniz == 0)
872  break;
873  /* rec rate is zero */
874  if( rate_recomb <= 0.)
875  break;
876 
877  rate_ioniz /= rate_recomb;
878  if( rate_ioniz > 1.)
879  {
880  /* this is peak ionization stage */
881  jmax = i;
882  rate_ioniz = 1.;
883  }
884  }
885 
886  /* replace its matrix elements with population sum */
887  for( long i=0; i<ion_range;i++ )
888  {
889  MAT(xmat,i,jmax) -= conserve_scale;
890  }
891  source[jmax] -= abund_total*conserve_scale;
892  }
893 
894 #if 0
895  if( false && nelem == ipHYDROGEN && (dynamics.lgAdvection || dynamics.lgTimeDependentStatic) && iteration>1 )
896  {
897  fprintf(ioQQQ,"DEBUGG Rate %.2f %.3e \n",fnzone,dynamics.Rate);
898  fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateIonizTot(nelem,0), ionbal.RateIonizTot(nelem,1) );
899  fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateRecomTot[nelem][0], ionbal.RateRecomTot[nelem][1]);
900  fprintf(ioQQQ," %.3e %.3e %.3e\n\n", dynamics.Source[nelem][0], dynamics.Source[nelem][1], dynamics.Source[nelem][2]);
901  }
902 
903  {
904  /* option to print matrix */
905  enum {DEBUG_LOC=false};
906  if( DEBUG_LOC && nzone==1 && lgPrintIt )
907  {
908  fprintf( ioQQQ,
909  " DEBUG ion_solver: nelem=%li ion_range=%li, limit=%li, nConv %li xmat follows\n",
910  nelem , ion_range,limit , conv.nTotalIoniz );
911  if( lgHomogeneous )
912  fprintf(ioQQQ , "Homogeneous \n");
913  for( long i=0; i<ion_range; ++i )
914  {
915  for( long j=0;j<ion_range;j++ )
916  {
917  fprintf(ioQQQ,"%e\t",MAT(xmat,i,j));
918  }
919  fprintf(ioQQQ,"\n");
920  }
921  fprintf(ioQQQ,"source follows\n");
922  for( long i=0; i<ion_range;i++ )
923  {
924  fprintf(ioQQQ,"%e\t",source[i]);
925  }
926  fprintf(ioQQQ,"\n");
927  fprintf(ioQQQ,"totsrc/totscl %g %g\n",totsrc,totscl);
928  }
929  }
930 #endif
931 
932 }
933 
934 STATIC void PrintRates( long nelem, bool lgNegPop, double abund_total, valarray<double> &auger, bool lgPrintIt )
935 {
936  DEBUG_ENTRY( "PrintRates()" );
937 
938  long ion;
939 
940  /* this should not happen */
941  if( lgNegPop )
942  {
943  fprintf( ioQQQ, " PROBLEM Negative population found for abundance of ionization stage of element %4.4s, ZONE=%4ld\n",
945 
946  fprintf( ioQQQ, " Populations were" );
947  for( ion=1; ion <= dense.IonHigh[nelem]+1; ion++ )
948  {
949  fprintf( ioQQQ, "%9.1e", dense.xIonDense[nelem][ion-1] );
950  }
951  fprintf( ioQQQ, "\n" );
952 
953  fprintf( ioQQQ, " destroy vector =" );
954  for( ion=1; ion <= dense.IonHigh[nelem]; ion++ )
955  {
956  fprintf( ioQQQ, "%9.1e", ionbal.RateIonizTot(nelem,ion-1) );
957  }
958  fprintf( ioQQQ, "\n" );
959 
960  fprintf( ioQQQ, " CTHeavy vector =" );
961  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
962  {
963  fprintf( ioQQQ, "%9.1e", atmdat.CharExcIonOf[ipHELIUM][nelem][ion] );
964  }
965  fprintf( ioQQQ, "\n" );
966 
967  fprintf( ioQQQ, " CharExcIonOf[ipHYDROGEN] vtr=" );
968  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
969  {
970  fprintf( ioQQQ, "%9.1e", atmdat.CharExcIonOf[ipHYDROGEN][nelem][ion] );
971  }
972  fprintf( ioQQQ, "\n" );
973 
974  fprintf( ioQQQ, " CollidRate vtr=" );
975  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
976  {
977  fprintf( ioQQQ, "%9.1e", ionbal.CollIonRate_Ground[nelem][ion][0] );
978  }
979  fprintf( ioQQQ, "\n" );
980 
981  /* photo rates per subshell */
982  fprintf( ioQQQ, " photo rates per subshell, ion\n" );
983  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
984  {
985  fprintf( ioQQQ, "%3ld", ion );
986  for( long ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
987  {
988  fprintf( ioQQQ, "%9.1e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
989  }
990  fprintf( ioQQQ, "\n" );
991  }
992 
993  /* now check out creation vector */
994  fprintf( ioQQQ, " create vector =" );
995  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
996  {
997  fprintf( ioQQQ, "%9.1e", ionbal.RateRecomTot[nelem][ion] );
998  }
999  fprintf( ioQQQ, "\n" );
1000  }
1001 
1002  /* option to print ionization and recombination arrays
1003  * prt flag set with print array print arrays command */
1004  if( lgPrintIt || prt.lgPrtArry[nelem] || lgNegPop )
1005  {
1006  const char* format = " %9.2e";
1007  /* say who we are, what we are doing .... */
1008  fprintf( ioQQQ,
1009  "\n %s ion_solver ion/rec rt [s-1] %s nz%.2f Te%.4e ne%.4e Tot abun:%.3e ion abun%.2e mole%.2e Search? %c\n",
1010  elementnames.chElementSym[nelem],
1011  elementnames.chElementName[nelem],
1012  fnzone,
1013  phycon.te ,
1014  dense.eden,
1015  dense.gas_phase[nelem],
1016  abund_total ,
1017  dense.xMolecules(nelem) , TorF(conv.lgSearch) );
1018  /* total ionization rate, all processes */
1019  fprintf( ioQQQ, " %s Ioniz total " ,elementnames.chElementSym[nelem]);
1020  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1021  {
1022  fprintf( ioQQQ, format, ionbal.RateIonizTot(nelem,ion) );
1023  }
1024  fprintf( ioQQQ, "\n" );
1025 
1026  /* sum of all creation processes */
1027  fprintf( ioQQQ, " %s Recom total " ,elementnames.chElementSym[nelem]);
1028  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1029  {
1030  fprintf( ioQQQ, format, ionbal.RateRecomTot[nelem][ion] );
1031  }
1032  fprintf( ioQQQ, "\n" );
1033 
1034  /* collisional ionization */
1035  fprintf( ioQQQ, " %s Coll ioniz " ,elementnames.chElementSym[nelem] );
1036  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1037  {
1038  double ColIoniz = ionbal.CollIonRate_Ground[nelem][ion][0];
1039  if( ion > nelem - NISO )
1040  {
1041  long ipISO = nelem-ion;
1042  ASSERT( ipISO >=0 && ipISO < NISO );
1043  if( dense.xIonDense[nelem][nelem-ipISO] > 0. )
1044  {
1045  ColIoniz *= iso_sp[ipISO][nelem].st[0].Pop();
1046  for( long ipLevel=1; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
1047  {
1048  ColIoniz += iso_sp[ipISO][nelem].fb[ipLevel].ColIoniz * dense.EdenHCorr * iso_sp[ipISO][nelem].st[ipLevel].Pop();
1049  }
1050  ColIoniz /= dense.xIonDense[nelem][nelem-ipISO];
1051  }
1052  }
1053  fprintf( ioQQQ, format, ColIoniz );
1054  }
1055  fprintf( ioQQQ, "\n" );
1056 
1057  /* UTA ionization */
1058  fprintf( ioQQQ, " %s UTA ioniz " ,elementnames.chElementSym[nelem] );
1059  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1060  {
1061  fprintf( ioQQQ, format, ionbal.UTA_ionize_rate[nelem][ion] );
1062  }
1063  fprintf( ioQQQ, "\n" );
1064 
1065  /* photo ionization */
1066  fprintf( ioQQQ, " %s Photoion snk" ,elementnames.chElementSym[nelem] );
1067  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1068  {
1069  double PhotIoniz = 0.;
1070  for( long ipShell = 0; ipShell < Heavy.nsShells[nelem][ion]; ipShell++ )
1071  PhotIoniz += ionbal.PhotoRate_Shell[nelem][ion][ipShell][0];
1072 
1073  // still don't have the total if stage is in one of the iso-sequences
1074  if( ion > nelem - NISO )
1075  {
1076  long ipISO = nelem-ion;
1077  ASSERT( ipISO >=0 && ipISO < NISO );
1078  if( dense.xIonDense[nelem][nelem-ipISO]>0 )
1079  {
1080  PhotIoniz *= iso_sp[ipISO][nelem].st[0].Pop();
1081  for( long ipLevel=1; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
1082  {
1083  PhotIoniz += iso_sp[ipISO][nelem].fb[ipLevel].gamnc * iso_sp[ipISO][nelem].st[ipLevel].Pop();
1084  }
1085  PhotIoniz /= dense.xIonDense[nelem][nelem-ipISO];
1086  }
1087  }
1088  fprintf( ioQQQ, format, PhotIoniz );
1089  }
1090  fprintf( ioQQQ, "\n" );
1091 
1092  /* photoionization source (source of this stage due to ionization of next stage) */
1093  fprintf( ioQQQ, " %s Photoion src" ,elementnames.chElementSym[nelem]);
1094  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1095  {
1096  double source = 0.;
1097  if( ion>0 )
1098  source = ionbal.RateIoniz[nelem][ion-1][ion];
1099 
1100  fprintf( ioQQQ, format, source );
1101  }
1102  fprintf( ioQQQ, "\n" );
1103 
1104  /* auger production (of this stage due to auger ionization of another) */
1105  fprintf( ioQQQ, " %s Auger src " ,elementnames.chElementSym[nelem]);
1106  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1107  {
1108  double source = 0.;
1109  if( ion>0 )
1110  source = auger[ion-1];
1111 
1112  fprintf( ioQQQ, format, source );
1113  }
1114  fprintf( ioQQQ, "\n" );
1115 
1116  /* secondary ionization */
1117  fprintf( ioQQQ, " %s Secon ioniz " ,elementnames.chElementSym[nelem]);
1118  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1119  {
1120  fprintf( ioQQQ, format,
1121  secondaries.csupra[nelem][ion] );
1122  }
1123  fprintf( ioQQQ, "\n" );
1124 
1125  /* grain ionization - not total rate but should be dominant process */
1126  fprintf( ioQQQ, " %s ion on grn " ,elementnames.chElementSym[nelem]);
1127  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1128  {
1129  fprintf( ioQQQ, format, gv.GrainChTrRate[nelem][ion][ion+1] );
1130  }
1131  fprintf( ioQQQ, "\n" );
1132 
1133  /* charge transfer ionization from chemistry */
1134  fprintf( ioQQQ, " %s ion xfr mol " ,elementnames.chElementSym[nelem]);
1135  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1136  {
1137  fprintf( ioQQQ, format, mole.xMoleChTrRate[nelem][ion][ion+1] );
1138  }
1139  fprintf( ioQQQ, "\n" );
1140 
1141  /* charge exchange ionization */
1142  fprintf( ioQQQ, " %s chr trn ion " ,elementnames.chElementSym[nelem] );
1143  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1144  {
1145  /* sum has units s-1 */
1146  double sum;
1147  long ipISO=nelem-ion;
1148 
1149  if( nelem < t_atmdat::NCX && nelem == ipISO )
1150  {
1151  sum = atmdat.CharExcIonTotal[nelem] * iso_sp[ipISO][nelem].st[0].Pop() / SDIV(dense.xIonDense[nelem][nelem-ipISO]);
1152  }
1153  else
1154  {
1155  sum = 0.0;
1156  for (long nelem1=0; nelem1 < t_atmdat::NCX; ++nelem1)
1157  sum += atmdat.CharExcIonOf[nelem1][nelem][ion] * dense.xIonDense[nelem1][1];
1158  }
1159 
1160  fprintf( ioQQQ, format, sum );
1161  }
1162  fprintf( ioQQQ, "\n" );
1163 
1164  /* radiative recombination */
1165  fprintf( ioQQQ, " %s radiati rec " ,elementnames.chElementSym[nelem]);
1166  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1167  {
1168  fprintf( ioQQQ, format, dense.eden*ionbal.RR_rate_coef_used[nelem][ion] );
1169  }
1170  fprintf( ioQQQ, "\n" );
1171 
1172  /* Badnell DR recombination */
1173  fprintf( ioQQQ, " %s drBadnel rec" ,elementnames.chElementSym[nelem]);
1174  for( ion=0; ion < min(nelem-1,dense.IonHigh[nelem]); ion++ )
1175  {
1176  fprintf( ioQQQ, format, dense.eden*ionbal.DR_Badnell_rate_coef[nelem][ion] );
1177  }
1178  fprintf( ioQQQ, "\n" );
1179 
1180  /* Cota rate */
1181  fprintf( ioQQQ, " %s CotaRate rec" ,elementnames.chElementSym[nelem]);
1182  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1183  {
1184  fprintf( ioQQQ, format, dense.eden*ionbal.CotaRate[ion] );
1185  }
1186  fprintf( ioQQQ, "\n" );
1187 
1188  /* grain recombination - not total but from next higher ion, should
1189  * be dominant */
1190  fprintf( ioQQQ, " %s rec on grn " ,elementnames.chElementSym[nelem]);
1191  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1192  {
1193  fprintf( ioQQQ, format, gv.GrainChTrRate[nelem][ion+1][ion] );
1194  }
1195  fprintf( ioQQQ, "\n" );
1196 
1197  /* charge transfer recombination from chemistry */
1198  fprintf( ioQQQ, " %s rec xfr mol " ,elementnames.chElementSym[nelem]);
1199  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1200  {
1201  fprintf( ioQQQ, format, mole.xMoleChTrRate[nelem][ion+1][ion] );
1202  }
1203  fprintf( ioQQQ, "\n" );
1204 
1205  /* charge exchange recombination */
1206  fprintf( ioQQQ, " %s chr trn rec " ,elementnames.chElementSym[nelem]);
1207  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1208  {
1209  double sum;
1210 
1211  if( nelem==ipHELIUM && ion==0 )
1212  {
1213  sum = atmdat.CharExcRecTotal[nelem];
1214  }
1215  else if( nelem==ipHYDROGEN && ion==0 )
1216  {
1217  sum = atmdat.CharExcRecTotal[nelem];
1218  }
1219  else
1220  {
1221  sum = 0.0;
1222  for (long nelem1=0; nelem1<t_atmdat::NCX; ++nelem1)
1223  {
1224  long ipISO=nelem1;
1225  sum += atmdat.CharExcRecTo[nelem1][nelem][ion] * iso_sp[ipISO][nelem1].st[0].Pop();
1226  }
1227  }
1228  fprintf( ioQQQ, format, sum );
1229  }
1230  fprintf( ioQQQ, "\n" );
1231 
1232  {
1233  /* sources from the chemistry network */
1234  fprintf(ioQQQ," %s molecule src",
1235  elementnames.chElementSym[nelem]);
1236  for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
1237  {
1238  fprintf(ioQQQ,format, mole.source[nelem][ion]/SDIV(dense.xIonDense[nelem][ion]) );
1239  }
1240  fprintf( ioQQQ, "\n" );
1241 
1242  /* sinks from the chemistry network */
1243  fprintf(ioQQQ," %s molecule snk",
1244  elementnames.chElementSym[nelem]);
1245  for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
1246  {
1247  fprintf(ioQQQ,format, mole.sink[nelem][ion] );
1248  }
1249  fprintf( ioQQQ, "\n" );
1250 
1252  {
1253  /* source from the dynamcs */
1254  fprintf(ioQQQ," %s dynamics src",
1255  elementnames.chElementSym[nelem]);
1256  for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
1257  {
1258  fprintf(ioQQQ,format, dynamics.Source[nelem][ion]/SDIV(dense.xIonDense[nelem][ion]) );
1259  }
1260  fprintf( ioQQQ, "\n" );
1261 
1262  /* sinks from the dynamics */
1263  fprintf(ioQQQ," %s dynamics snk",
1264  elementnames.chElementSym[nelem]);
1265  for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
1266  {
1267  fprintf(ioQQQ,format, dynamics.Rate );
1268  }
1269  fprintf( ioQQQ, "\n" );
1270  }
1271  }
1272 
1273  /* the "new" abundances the resulted from the previous ratio */
1274  fprintf( ioQQQ, " %s Abun [cm-3] " ,elementnames.chElementSym[nelem] );
1275  for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
1276  {
1277  fprintf( ioQQQ, format, dense.xIonDense[nelem][ion] );
1278  }
1279  fprintf( ioQQQ, "\n" );
1280  }
1281 
1282  if( lgNegPop )
1283  {
1284  ContNegative();
1285  ShowMe();
1287  }
1288 
1289  return;
1290 }
1291 
1292 /*
1293 
1294  Solve an ionization level system with specified ionization and
1295  recombination rates between neighboring ions, and additional sink
1296  and source terms. The sink array is overwritten, and the results
1297  appear in the source array.
1298 
1299  Written in matrix form, the algorithm is equivalent to the
1300  tridiagonal algorithm in Numerical Recipes applied to:
1301 
1302  / i_0+a_0 -r_0 . . . \ / x_0 \ / s_0 \
1303  | -i_0 i_1+a_1+r_0 -r_1 . . | | x_1 | | s_1 |
1304  | . -i_1 i_2+a_2+r_1 -r_2 . | | x_2 | | s_2 |
1305  | . . (etc....) | | ... | = | ... |
1306  \ . . . / \ / \ /
1307 
1308  where i, r are the ionization and recombination rates, s is the
1309  source rate and a is the sink rate.
1310 
1311  This matrix is diagonally dominant only when the sink terms are
1312  large -- the alternative method coded here prevents rounding error
1313  in the diagonal terms disturbing the solution when this is not the
1314  case.
1315 
1316 */
1317 
1318 /* solveions tridiagonal solver but optimized for structure of balance matrix */
1319 void solveions(double *ion, double *rec, double *snk, double *src,
1320  long int nlev, long int nmax)
1321 {
1322  double kap, bet;
1323  long int i;
1324 
1325  DEBUG_ENTRY( "solveions()" );
1326 
1327  if(nmax != -1)
1328  {
1329  /* Singular case */
1330  src[nmax] = 1.;
1331  for(i=nmax;i<nlev-1;i++)
1332  src[i+1] = src[i]*ion[i]/rec[i];
1333  for(i=nmax-1;i>=0;i--)
1334  src[i] = src[i+1]*rec[i]/ion[i];
1335  }
1336  else
1337  {
1338  kap = snk[0];
1339  for(i=0;i<nlev-1;i++)
1340  {
1341  bet = ion[i]+kap;
1342  if(bet == 0.)
1343  {
1344  fprintf(ioQQQ,"Ionization solver error\n");
1346  }
1347  bet = 1./bet;
1348  src[i] *= bet;
1349  src[i+1] += ion[i]*src[i];
1350  snk[i] = bet*rec[i];
1351  kap = kap*snk[i]+snk[i+1];
1352  }
1353  bet = kap;
1354  if(bet == 0.)
1355  {
1356  fprintf(ioQQQ,"Ionization solver error\n");
1358  }
1359  src[i] /= bet;
1360 
1361  for(i=nlev-2;i>=0;i--)
1362  {
1363  src[i] += snk[i]*src[i+1];
1364  }
1365  }
1366 }
1367 
1368 void ion_wrapper( long nelem )
1369 {
1370 
1371  DEBUG_ENTRY( "ion_wrapper()" );
1372 
1373  ASSERT( nelem >= ipHYDROGEN );
1374  ASSERT( nelem < LIMELM );
1375 
1376  switch( nelem )
1377  {
1378  case ipHYDROGEN:
1379  IonHydro();
1380  break;
1381  case ipHELIUM:
1382  IonHelium();
1383  break;
1384  default:
1385  IonNelem(false,nelem);
1386  break;
1387  }
1388 
1389  if( trace.lgTrace && dense.lgElmtOn[nelem] && trace.lgHeavyBug )
1390  {
1391  fprintf( ioQQQ, " ion_wrapper returns; %s fracs = ", elementnames.chElementSym[nelem] );
1392  for( long ion = 0; ion<=nelem+1; ion++ )
1393  fprintf( ioQQQ,"%10.3e ", dense.xIonDense[nelem][ion]/dense.gas_phase[nelem] );
1394  fprintf( ioQQQ, "\n" );
1395  }
1396 
1398 
1399  return;
1400 }
double ** DR_Badnell_rate_coef
Definition: ionbal.h:200
STATIC double get_total_abundance_ions(long int nelem)
Definition: ion_solver.cpp:507
t_atmdat atmdat
Definition: atmdat.cpp:6
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
realnum xMolecules(long nelem)
Definition: dense.h:88
qList st
Definition: iso.h:482
double EdenHCorr
Definition: dense.h:227
bool lgHeavyBug
Definition: trace.h:21
t_Heavy Heavy
Definition: heavy.cpp:5
realnum SetIoniz[LIMELM][LIMELM+1]
Definition: dense.h:172
const int NISO
Definition: cddefines.h:311
ConvergenceCounter register_(const string name)
Definition: conv.cpp:87
long int IonHigh[LIMELM+1]
Definition: dense.h:130
char TorF(bool l)
Definition: cddefines.h:749
bool lgAdvection
Definition: dynamics.h:66
bool lgTimeDependentStatic
Definition: dynamics.h:102
void ion_solver(long int nelem, bool lgPrintIt)
Definition: ion_solver.cpp:59
t_conv conv
Definition: conv.cpp:5
double elecsrc[LIMELM]
Definition: ionbal.h:248
t_phycon phycon
Definition: phycon.cpp:6
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition: cddefines.h:904
realnum GasPhaseAbundErrorAllowed
Definition: conv.h:281
bool lgSetIoniz[LIMELM]
Definition: dense.h:167
FILE * ioQQQ
Definition: cddefines.cpp:7
long int nzone
Definition: cddefines.cpp:14
t_dynamics dynamics
Definition: dynamics.cpp:42
vector< freeBound > fb
Definition: iso.h:481
#define MIN2(a, b)
Definition: cddefines.h:803
void IonNelem(bool lgPrintIt, long int nelem)
Definition: ion_nelem.cpp:9
t_dense dense
Definition: global.cpp:15
const double MAX_DENSITY
Definition: cddefines.h:319
static t_yield & Inst()
Definition: cddefines.h:209
STATIC void find_solution(long nelem, long ion_range, valarray< double > &xmat, valarray< double > &source)
Definition: ion_solver.cpp:212
t_elementnames elementnames
Definition: elementnames.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
double ** RateRecomTot
Definition: ionbal.h:194
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
double CharExcRecTotal[NCX]
Definition: atmdat.h:310
long int iteration
Definition: cddefines.cpp:16
double CharExcIonTotal[NCX]
Definition: atmdat.h:310
bool lgSearch
Definition: conv.h:168
t_trace trace
Definition: trace.cpp:5
double ** sink
Definition: mole.h:464
void broken(void)
Definition: service.cpp:1067
t_ionbal ionbal
Definition: ionbal.cpp:8
bool lgPrtArry[LIMELM]
Definition: prt.h:236
STATIC void store_new_densities(long nelem, long ion_range, long ion_low, double *source)
Definition: ion_solver.cpp:468
STATIC void PrintRates(long nelem, bool lgNegPop, double abund_total, valarray< double > &auger, bool lgPrintIt)
Definition: ion_solver.cpp:934
long int IonLow[LIMELM+1]
Definition: dense.h:129
long nelec_eject(long n, long i, long ns) const
Definition: yield.h:55
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
realnum elec_eject_frac(long n, long i, long ns, long ne) const
Definition: yield.h:48
void ContNegative(void)
STATIC void fill_ext_src_and_snk(long nelem, long ion_low, long ion_range, valarray< double > &xmat, valarray< double > &source)
Definition: ion_solver.cpp:750
#define STATIC
Definition: cddefines.h:118
bool lgTrace
Definition: trace.h:12
void iso_departure_coefficients(long ipISO, long nelem)
Definition: iso_solve.cpp:401
double ** source
Definition: mole.h:464
void iso_satellite_update(long nelem)
t_mole_local mole
Definition: mole.cpp:8
void iso_charge_transfer_update(long nelem)
double ** UTA_ionize_rate
Definition: ionbal.h:174
realnum *** xMoleChTrRate
Definition: mole.h:466
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
STATIC void HomogeneousSource(long nelem, long ion_low, long ion_range, valarray< double > &xmat, valarray< double > &source, double abund_total)
Definition: ion_solver.cpp:772
#define MAT(M_, I_, J_)
Definition: ion_solver.cpp:51
int32 solve_system(const valarray< double > &a, valarray< double > &b, long int n, error_print_t error_print)
bool lgEquilibrium
Definition: dynamics.h:180
#define cdEXIT(FAIL)
Definition: cddefines.h:482
long min(int a, long b)
Definition: cddefines.h:762
bool lgGrainPhysicsOn
Definition: grainvar.h:479
void iso_solve(long ipISO, long nelem, double &maxerr)
Definition: iso_solve.cpp:99
char chElementNameShort[LIMELM][CHARS_ELEMENT_NAME_SHORT]
Definition: elementnames.h:21
t_prt prt
Definition: prt.cpp:14
long int nTotalIoniz
Definition: conv.h:159
double **** PhotoRate_Shell
Definition: ionbal.h:112
bool lgElmtOn[LIMELM]
Definition: dense.h:160
STATIC void fill_array(long int nelem, long ion_range, valarray< double > &xmat, valarray< double > &auger)
Definition: ion_solver.cpp:538
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
realnum gas_phase[LIMELM]
Definition: dense.h:76
void ion_trim_small(long nelem, double abund_total)
Definition: ion_trim.cpp:115
double *** CollIonRate_Ground
Definition: ionbal.h:121
#define ASSERT(exp)
Definition: cddefines.h:613
STATIC bool lgTrivialSolution(long nelem, double abund_total)
Definition: ion_solver.cpp:187
int lgGrainIonRecom
Definition: ionbal.h:226
double Rate
Definition: dynamics.h:77
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:308
double elecsnk[LIMELM]
Definition: ionbal.h:248
double *** RateIoniz
Definition: ionbal.h:180
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
#define isnan
Definition: cddefines.h:659
const int ipHELIUM
Definition: cddefines.h:350
STATIC void clean_up(long nelem, double abund_total)
Definition: ion_solver.cpp:482
void solveions(double *ion, double *rec, double *snk, double *src, long int nlev, long int nmax)
double eden
Definition: dense.h:201
void iso_renorm(long nelem, long ipISO, double &renorm)
Definition: iso_solve.cpp:310
void IonHydro()
Definition: iso_solve.cpp:187
#define MAX2(a, b)
Definition: cddefines.h:824
double RateIonizTot(long nelem, long ion) const
Definition: ionbal.cpp:224
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
void IonHelium(void)
Definition: ion_helium.cpp:11
void ion_wrapper(long nelem)
bool lgCTOn
Definition: atmdat.h:325
realnum ** csupra
Definition: secondaries.h:33
double ** Source
Definition: dynamics.h:80
sys_float SDIV(sys_float x)
Definition: cddefines.h:1002
long int n_initial_relax
Definition: dynamics.h:132
const int ipCARBON
Definition: cddefines.h:354
long int numLevels_max
Definition: iso.h:524
char chElementName[LIMELM][CHARS_ELEMENT_NAME]
Definition: elementnames.h:17
#define fixit(a)
Definition: cddefines.h:417
bool lgElemsConserved(void)
Definition: dense.cpp:119
GrainVar gv
Definition: grainvar.cpp:5
t_secondaries secondaries
Definition: secondaries.cpp:5
double fnzone
Definition: cddefines.cpp:15
double CharExcIonOf[NCX][LIMELM][LIMELM+1]
Definition: atmdat.h:300
void ShowMe(void)
Definition: service.cpp:205
double te
Definition: phycon.h:21
bool lgDustOn() const
Definition: grainvar.h:475
const int ipHYDROGEN
Definition: cddefines.h:349
realnum GrainChTrRate[LIMELM][LIMELM+1][LIMELM+1]
Definition: grainvar.h:541
double CharExcRecTo[NCX][LIMELM][LIMELM+1]
Definition: atmdat.h:300
void iso_set_ion_rates(long ipISO, long nelem)
Definition: iso_level.cpp:811
long int numLevels_local
Definition: iso.h:529
STATIC void renorm_solution(long nelem, long ion_range, long ion_low, double *source, double abund_total, bool *lgNegPop)
Definition: ion_solver.cpp:318
realnum CotaRate[LIMELM]
Definition: ionbal.h:239
double ** RR_rate_coef_used
Definition: ionbal.h:210