cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
iso_level.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*iso_level solve for iso-sequence level populations */
4 #include "cddefines.h"
5 #include "atmdat.h"
6 #include "continuum.h"
7 #include "conv.h"
8 #include "doppvel.h"
9 #include "dynamics.h"
10 #include "elementnames.h"
11 #include "grainvar.h"
12 #include "he.h"
13 #include "ionbal.h"
14 #include "iso.h"
15 #include "opacity.h"
16 #include "phycon.h"
17 #include "rfield.h"
18 #include "trace.h"
19 #include "mole.h"
20 #include "freebound.h"
21 #include "two_photon.h"
22 #include "dense.h"
23 #include "vectorize.h"
24 #include "prt.h"
25 #include "iterations.h"
26 
27 
29  const long int ipISO, const long int nelem);
30 
31 /* solve for level populations */
32 void iso_level( const long int ipISO, const long int nelem, double &renorm,
33  bool lgPrtMatrix )
34 {
35  long int ipHi,
36  ipLo,
37  i,
38  level,
39  level_error;
40 
41  t_iso_sp* sp = &iso_sp[ipISO][nelem];
42  const long int numlevels_local = sp->numLevels_local;
43 
44  double BigError;
45 
46  int32 nerror;
47  double HighestPColOut = 0.;
48  bool lgNegPop=false;
49  static vector<int32> ipiv;
50  ipiv.resize(numlevels_local);
51  /* this block of variables will be obtained and freed here */
52  double source=0., sink=0.;
53  static vector<double> PopPerN;
54  PopPerN.resize(sp->n_HighestResolved_local+1);
55 
56  DEBUG_ENTRY( "iso_level()" );
57  renorm = 1.;
58 
59  /* check that we were called with valid charge */
60  ASSERT( nelem >= ipISO );
61  ASSERT( nelem < LIMELM );
62 
63  /* these two collision rates must be the same or we are in big trouble,
64  * since used interchangeably */
65  ASSERT( ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0]< SMALLFLOAT ||
66  fabs( (sp->fb[0].ColIoniz* dense.EdenHCorr) /
67  SDIV(ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0] ) - 1.) < 0.001 );
68 
69  if( (dense.IonHigh[nelem] >= nelem - ipISO) &&
70  (dense.IonLow[nelem] <= nelem - ipISO) &&
71  ionbal.RateRecomIso[nelem][ipISO] > 0. )
72  {
73  /* get simple estimate of atom to ion ratio */
74  sp->xIonSimple = ionbal.RateIonizTot(nelem,nelem-ipISO)/ionbal.RateRecomIso[nelem][ipISO];
75  }
76  else
77  {
78  sp->xIonSimple = 0.;
79  }
80 
81  /* which case atom to solve??? */
82  if( dense.xIonDense[nelem][nelem+1-ipISO] < SMALLFLOAT || sp->xIonSimple < 1e-35 )
83  {
84  /* don't bother if no ionizing radiation */
85  strcpy( sp->chTypeAtomUsed, "zero " );
86  if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) )
87  {
88  fprintf( ioQQQ, " iso_level iso=%2ld nelem=%2ld simple II/I=%10.2e so not doing equilibrium, doing %s.\n",
89  ipISO, nelem, sp->xIonSimple, sp->chTypeAtomUsed );
90  }
91 
92  /* total ionization will just be the ground state */
93  lgNegPop = false;
94  sp->st[0].Pop() = dense.xIonDense[nelem][nelem-ipISO];
95  for( long n=1; n < numlevels_local; n++ )
96  {
97  sp->st[n].Pop() = 0.;
98  }
99  sp->qTot2S = 0.;
100  }
101  else
102  {
103  /* fill in recombination vector - values were set in iso_ionize_recombine.cpp */
104  static vector<double> creation, error, work, Save_creation;
105  creation.resize(numlevels_local);
106  error.resize(numlevels_local);
107  work.resize(numlevels_local);
108  Save_creation.resize(numlevels_local);
109 
110  ASSERT( dense.xIonDense[nelem][nelem+1-ipISO] >= 0.f );
111  for( level=0; level < numlevels_local; level++ )
112  {
113  /* total recombination from once more ionized [cm-3 s-1] */
114  creation[level] = sp->fb[level].RateCont2Level * dense.xIonDense[nelem][nelem+1-ipISO];
115  }
116 
117  double ctsource=0.0, ctsink=0.0, ctrec=0.0;
118  /* now charge transfer - all into/from ground, two cases, H and not H */
119  if( nelem==ipHYDROGEN )
120  {
121  /* charge transfer, hydrogen onto everything else */
122  /* charge exchange recombination */
123  ctrec += atmdat.CharExcRecTotal[nelem];
124  ctsink += atmdat.CharExcIonTotal[nelem];
125  }
126  else if( nelem==ipHELIUM && ipISO==ipHE_LIKE )
127  {
128  /* this is recom of He due to ct with all other gas constituents */
129  ctrec += atmdat.CharExcRecTotal[nelem];
130  ctsink += atmdat.CharExcIonTotal[nelem];
131  }
132  else
133  {
134  for (long nelem1=ipHYDROGEN; nelem1 < t_atmdat::NCX; ++nelem1)
135  {
136  long ipISO=nelem1;
137  ctrec += atmdat.CharExcRecTo[nelem1][nelem][nelem-ipISO]*iso_sp[ipISO][nelem1].st[0].Pop();
138  ctsink += atmdat.CharExcIonOf[nelem1][nelem][nelem-ipISO]*dense.xIonDense[nelem1][1];
139  }
140  }
141  ctsource += ctrec*dense.xIonDense[nelem][nelem+1-ipISO];
142 
143  if ( nelem > ipISO )
144  {
145  double ction=0.0;
146  if( nelem==ipHELIUM )
147  {
148  ctsink += atmdat.CharExcRecTotal[nelem];
149  ction += atmdat.CharExcIonTotal[nelem];
150  }
151  else
152  {
153  for (long nelem1=ipHYDROGEN; nelem1 < t_atmdat::NCX; ++nelem1)
154  {
155  long ipISO1=nelem1;
156  ctsink += atmdat.CharExcRecTo[nelem1][nelem][nelem-(ipISO+1)]*iso_sp[ipISO1][nelem1].st[0].Pop();
157  ction += atmdat.CharExcIonOf[nelem1][nelem][nelem-(ipISO+1)]*dense.xIonDense[nelem1][1];
158  }
159  }
160  ctsource += ction * dense.xIonDense[nelem][nelem-(ipISO+1)];
161  }
162 
163  /* now do the 2D array */
165  z.alloc(numlevels_local,numlevels_local);
166  z.zero();
167 
168  /* this branch is main solver, full level populations
169  * assert since this code must change if NISO ever increased */
170  ASSERT( NISO == 2 );
171  long SpinUsed[NISO] = {2,1};
172  long indexNmaxP =
173  sp->QuantumNumbers2Index[ sp->n_HighestResolved_local ][1][SpinUsed[ipISO]];
174 
175  strcpy( sp->chTypeAtomUsed, "Popul" );
176  if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) )
177  {
178  fprintf( ioQQQ, " iso_level iso=%2ld nelem=%2ld doing regular matrix inversion, %s\n",
179  ipISO, nelem, sp->chTypeAtomUsed );
180  }
181 
182  //qList::const_iterator StElm = StatesElemNEW[nelem][nelem-ipISO].begin();
183  qList::const_iterator StElm = sp->st.begin();
184  static vector<double> Boltzmann_overg;
185  Boltzmann_overg.resize(numlevels_local-1);
186  for (unsigned lev = 0; lev < Boltzmann_overg.size(); ++lev)
187  Boltzmann_overg[lev] = 1.0/(double)StElm[lev].g();
188 
189  /* >>chng 05 dec 21, rm eden to make into rate coefficient */
190  sp->qTot2S = sp->fb[1].ColIoniz;
191 
192  static vector<double> coll_down, RadDecay, pump;
193  coll_down.resize(numlevels_local);
194  RadDecay.resize(numlevels_local);
195  pump.resize(numlevels_local);
196  avx_ptr<double> arg(1,numlevels_local), bstep(1,numlevels_local);
197 
198  for( level=1; level < numlevels_local; level++ )
199  {
200  fixit("Wouldn't need to mask this out if levels were in order");
201  arg[level] = -max((StElm[level].energy().K()-StElm[level-1].energy().K()) / phycon.te, -38.);
202  }
203  vexp( arg.ptr0(), bstep.ptr0(), 1, numlevels_local );
204 
206 
207  enum { DEBUG_RATES = false };
208 
209  if( DEBUG_RATES && iterations.lgLastIt )
210  fprintf( stdout, "# ipISO\tnelem\tlevel\tcollDown\tcollIonz\tradRecom\tRadDecay\n" );
211 
212  /* master balance equation, use when significant population */
213  for( level=0; level < numlevels_local; level++ )
214  {
215  double coll_down_total = 0.;
216 
217  /* all process depopulating level and placing into the continuum
218  * this does NOT include grain charge transfer ionization, added below */
219  z[level][level] += sp->fb[level].RateLevel2Cont;
220 
221  if (level != 0)
222  {
223  for ( ipLo = 0; ipLo < level; ++ipLo )
224  Boltzmann_overg[ipLo] *= bstep[level];
225  }
226 
227  /* all processes populating level from below */
228  for( ipLo=0; ipLo < level; ipLo++ )
229  {
230  coll_down[ipLo] = sp->trans(level,ipLo).Coll().ColUL( colld );
231  if( DEBUG_RATES )
232  {
233  coll_down_total += coll_down[ipLo];
234  }
235 
236  if ( rfield.lgPlasNu && sp->trans(level,ipLo).EnergyRyd()<rfield.plsfrq )
237  {
238  RadDecay[ipLo] = iso_ctrl.SmallA;
239  pump[ipLo] = iso_ctrl.SmallA;
240  }
241  else
242  {
243  RadDecay[ipLo] = MAX2( iso_ctrl.SmallA, sp->trans(level,ipLo).Emis().Aul()*
244  (sp->trans(level,ipLo).Emis().Ploss()) );
245  pump[ipLo] = MAX2( iso_ctrl.SmallA, sp->trans(level,ipLo).Emis().pump() );
246  }
247  }
248 
249  if( iso_ctrl.lgRandErrGen[ipISO] )
250  {
251  for( ipLo=0; ipLo < level; ipLo++ )
252  {
253  coll_down[ipLo] *= sp->ex[level][ipLo].ErrorFactor[IPCOLLIS];
254  RadDecay[ipLo] *= sp->ex[level][ipLo].ErrorFactor[IPRAD];
255  pump[ipLo] *= sp->ex[level][ipLo].ErrorFactor[IPRAD];
256  }
257  }
258 
259  double glev = (double)StElm[level].g(), rglev = 1.0/glev;
260  for( ipLo=0; ipLo < level; ipLo++ )
261  {
262  double coll_up = coll_down[ipLo] * glev *
263  Boltzmann_overg[ipLo];
264 
265  z[ipLo][ipLo] += coll_up + pump[ipLo] ;
266  z[ipLo][level] -= coll_up + pump[ipLo] ;
267 
268  double pump_down = pump[ipLo] *
269  (double)StElm[ipLo].g() * rglev;
270 
271  z[level][level] += RadDecay[ipLo] + pump_down + coll_down[ipLo];
272  z[level][ipLo] -= RadDecay[ipLo] + pump_down + coll_down[ipLo];
273  if( level == indexNmaxP )
274  {
275  HighestPColOut += coll_down[ipLo];
276  }
277  if( ipLo == indexNmaxP )
278  {
279  HighestPColOut += coll_up;
280  }
281 
282  /* find total collisions out of 2^3S to singlets. */
283  if( (level == 1) && (ipLo==0) )
284  {
285  sp->qTot2S += coll_down[ipLo]/dense.EdenHCorr;
286  }
287  if( (ipLo == 1) && (ipISO==ipH_LIKE || StElm[level].S()==1) )
288  {
289  sp->qTot2S += coll_up/dense.EdenHCorr;
290  }
291  }
292 
293  if( DEBUG_RATES && iterations.lgLastIt )
294  {
295  fprintf( stdout,
296  "%2ld\t%2ld\t%2ld\t"
297  "%.4e\t"
298  "%.4e\t"
299  "%.4e\t"
300  "%.4e\n",
301  ipISO, nelem, level,
302  coll_down_total,
303  iso_sp[ipISO][nelem].fb[level].ColIoniz * dense.eden,
304  iso_sp[ipISO][nelem].fb[level].RadRecomb[0] * dense.eden,
305  1. / iso_sp[ipISO][nelem].st[level].lifetime() );
306  }
307  }
308 
309  if( DEBUG_RATES && iterations.lgLastIt )
310  fprintf( stdout, "\n\n" );
311 
312  if( ipISO == nelem )
313  {
314  /* iso_ctrl.lgCritDensLMix[ipISO] is a flag used to print warning if density is
315  * too low for first collapsed level to be l-mixed. Check is if l-mixing collisions
316  * out of highest resolved singlet P are greater than sum of transition probs out. */
317  if( HighestPColOut < 1./sp->st[indexNmaxP].lifetime() )
318  {
319  iso_ctrl.lgCritDensLMix[ipISO] = false;
320  }
321  }
322 
324  ASSERT( ipISO <= ipHE_LIKE );
325  for( vector<two_photon>::iterator tnu = sp->TwoNu.begin(); tnu != sp->TwoNu.end(); ++tnu )
326  {
327  // induced two photon emission - special because upward and downward are
328  // not related by ratio of statistical weights
329  // iso.lgInd2nu_On is controlled with SET IND2 ON/OFF command
330 
331  fixit("need Pesc or the equivalent to multiply AulTotal?");
332  // downward rate
333  z[tnu->ipHi][tnu->ipLo] -= tnu->AulTotal;
334  z[tnu->ipHi][tnu->ipHi] += tnu->AulTotal;
335 
336  }
337 
338  if (iso_ctrl.lgInd2nu_On)
339  {
340  for( vector<two_photon>::iterator tnu = sp->TwoNu.begin(); tnu != sp->TwoNu.end(); ++tnu )
341  {
342  // induced two photon emission - special because upward and downward are
343  // not related by ratio of statistical weights
344  // iso.lgInd2nu_On is controlled with SET IND2 ON/OFF command
345 
346  z[tnu->ipHi][tnu->ipLo] -= tnu->induc_dn;
347  z[tnu->ipHi][tnu->ipHi] += tnu->induc_dn;
348 
349  // upward rate
350  z[tnu->ipLo][tnu->ipHi] -= tnu->induc_up;
351  z[tnu->ipLo][tnu->ipLo] += tnu->induc_up;
352  }
353  }
354 
355  /* grain charge transfer recombination and ionization to ALL other stages */
356  for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
357  {
358  if( ion!=nelem-ipISO )
359  {
360  source += gv.GrainChTrRate[nelem][ion][nelem-ipISO] *
361  dense.xIonDense[nelem][ion];
362  sink += gv.GrainChTrRate[nelem][nelem-ipISO][ion];
363 
364  source += mole.xMoleChTrRate[nelem][ion][nelem-ipISO] *
365  dense.xIonDense[nelem][ion] * atmdat.lgCTOn;
366  sink += mole.xMoleChTrRate[nelem][nelem-ipISO][ion] * atmdat.lgCTOn;
367 
368  }
369  }
370 
371  /* add in source and sink terms from molecular network. */
372  source += mole.source[nelem][nelem-ipISO];
373  sink += mole.sink[nelem][nelem-ipISO];
374 
375 #if 1
376  /* >>chng 02 Sep 06 rjrw -- all elements have these terms */
377  /*>>>chng 02 oct 01, only include if lgAdvection or lgTimeDependentStatic is set */
380  dynamics.Rate != 0.0 &&
381  !dynamics.lgEquilibrium && dynamics.lgISO[ipISO] )
382  {
383  /* add in advection - these terms normally zero */
384  source += dynamics.Source[nelem][nelem-ipISO];
385  /* >>chng 02 Sep 06 rjrw -- advective term not recombination */
386  sink += dynamics.Rate;
387  }
388 #else
389  /*>>>chng 02 oct 01, only include if lgAdvection or lgTimeDependentStatic is set */
392  dynamics.Rate != 0.0 &&
394  {
395  for( level=0; level < numlevels_local; level++ )
396  {
397  creation[level] += dynamics.StatesElem[nelem][nelem-ipISO][level];
398  z[level][level] += dynamics.Rate;
399  }
400  }
401 #endif
402 
403  /* ionization from/recombination from lower ionization stages */
404  for(long ion_from=dense.IonLow[nelem]; ion_from < MIN2( dense.IonHigh[nelem], nelem-ipISO ) ; ion_from++ )
405  {
406  if( ionbal.RateIoniz[nelem][ion_from][nelem-ipISO] >= 0. )
407  {
408  /* ionization from lower ionization stages, cm-3 s-1 */
409  source += ionbal.RateIoniz[nelem][ion_from][nelem-ipISO] * dense.xIonDense[nelem][ion_from];
410  }
411  /* recombination to next lower ionization stage, s-1 */
412  if( ion_from == nelem-1-ipISO )
413  sink += ionbal.RateRecomTot[nelem][ion_from];
414  }
415 
416  ASSERT( source >= 0.f );
417  if (0)
418  {
419  /*
420  * Collisional ionization and photoionization can only be to the ground state in H iso-sequences
421  * to conserve energy.
422  */
423  creation[0] += source;
424  /*for( level=0; level < numlevels_local; level++ )
425  {
426  z[level][level] += sink;
427  }*/
428  /*recombination is only done from the ground state */
429  z[0][0] += sink;
430  }
431  else
432  {
433  // Try Boltzmann weighting to capture LTE limit correctly
434  t_iso_sp* sp = &iso_sp[ipISO][nelem];
435  double partfun=0.0;
436  for ( level = 0; level < numlevels_local; level++ )
437  {
438  partfun += sp->st[level].Boltzmann()*sp->st[level].g();
439  }
440  source /= partfun;
441  for( level=0; level < numlevels_local; level++ )
442  {
443  creation[level] += source*
444  sp->st[level].Boltzmann()*sp->st[level].g();
445  z[level][level] += sink;
446  }
447  }
448  creation[0] += ctsource;
449  z[0][0] += ctsink;
450 
451  /* >>chng 04 nov 30, atom XX-like collisions off turns this off too */
452  if( sp->trans(iso_ctrl.nLyaLevel[ipISO],0).Coll().rate_lu_nontherm() * iso_ctrl.lgColl_excite[ipISO] > 0. )
453  {
454  /* now add on supra thermal excitation */
455  for( level=1; level < numlevels_local; level++ )
456  {
457  double RateUp , RateDown;
458 
459  RateUp = sp->trans(level,0).Coll().rate_lu_nontherm();
460  RateDown = RateUp * (double)sp->st[ipH1s].g() /
461  (double)sp->st[level].g();
462 
463  /* total rate out of lower level */
464  z[ipH1s][ipH1s] += RateUp;
465 
466  /* rate from the upper level to ground */
467  z[level][ipH1s] -= RateDown;
468 
469  /* rate from ground to upper level */
470  z[ipH1s][level] -= RateUp;
471 
472  z[level][level] += RateDown;
473  }
474  }
475 
476  /* ===================================================================
477  *
478  * at this point all matrix elements have been established
479  *
480  * ==================================================================== */
481  /* save matrix, this allocates SaveZ */
482  SaveZ = z;
483 
484  for( ipLo=0; ipLo < numlevels_local; ipLo++ )
485  Save_creation[ipLo] = creation[ipLo];
486 
487  if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) )
488  {
489  const long numlevels_print = numlevels_local;
490  fprintf( ioQQQ, " pop level others => (iso_level)\n" );
491  for( long n=0; n < numlevels_print; n++ )
492  {
493  fprintf( ioQQQ, " %s %s %2ld", iso_ctrl.chISO[ipISO], elementnames.chElementNameShort[nelem], n );
494  for( long j=0; j < numlevels_print; j++ )
495  {
496  fprintf( ioQQQ,"\t%.9e", z[j][n] );
497  }
498  fprintf( ioQQQ, "\n" );
499  }
500  fprintf(ioQQQ," recomb ct %.2e co %.2e hectr %.2e hctr %.2e\n",
502  findspecieslocal("CO")->den ,
503  atmdat.CharExcRecTo[ipHELIUM][nelem][nelem-ipISO]*iso_sp[ipHE_LIKE][ipHELIUM].st[0].Pop() ,
504  atmdat.CharExcRecTo[ipHYDROGEN][nelem][nelem-ipISO]*iso_sp[ipH_LIKE][ipHYDROGEN].st[0].Pop() );
505  fprintf(ioQQQ," recomb ");
506  for( long n=0; n < numlevels_print; n++ )
507  {
508  fprintf( ioQQQ,"\t%.9e", creation[n] );
509  }
510  fprintf( ioQQQ, "\n" );
511  }
512 
513 
514  if( lgPrtMatrix )
515  {
516  valarray<double> c( get_ptr(creation), creation.size() );
517  prt.matrix.prtRates( numlevels_local, z, c );
518  }
519 
520  nerror = 0;
521 
522  getrf_wrapper(numlevels_local,numlevels_local,
523  z.data(),numlevels_local,&ipiv[0],&nerror);
524 
525  getrs_wrapper('N',numlevels_local,1,z.data(),numlevels_local,&ipiv[0],
526  &creation[0],numlevels_local,&nerror);
527 
528  if( nerror != 0 )
529  {
530  fprintf( ioQQQ, " iso_level: dgetrs finds singular or ill-conditioned matrix\n" );
532  }
533 
534  /* check whether solution is valid */
535  /* >>chng 06 aug 28, both of these from numLevels_max to _local. */
536  for( level=ipH1s; level < numlevels_local; level++ )
537  {
538  double qn = 0., qx = 0.;
539  error[level] = 0.;
540  for( long n=ipH1s; n < numlevels_local; n++ )
541  {
542  double q = SaveZ[n][level]*creation[n];
543 
544  /* remember the largest size of element in sum to div by below */
545  if ( q > qx )
546  qx = q;
547  else if (q < qn)
548  qn = q;
549 
550  error[level] += q;
551  }
552 
553  if (-qn > qx)
554  qx = -qn;
555 
556  if( qx > 0. )
557  {
558  error[level] = (error[level] - Save_creation[level])/qx;
559  }
560  else
561  {
562  error[level] = 0.;
563  }
564  }
565 
566  /* remember largest residual in matrix inversion */
567  BigError = -1.;
568  level_error = -1;
569  /* >>chng 06 aug 28, from numLevels_max to _local. */
570  for( level=ipH1s; level < numlevels_local; level++ )
571  {
572  double abserror;
573  abserror = fabs( error[level]);
574  /* this will be the largest residual in the matrix inversion */
575  if( abserror > BigError )
576  {
577  BigError = abserror;
578  level_error = level;
579  }
580  }
581 
582  /* matrix inversion should be nearly as good as the accuracy of a double,
583  * but demand that it is better than epsilon for a float */
584  if( BigError > FLT_EPSILON )
585  {
586  if( !conv.lgSearch )
587  fprintf(ioQQQ," PROBLEM" );
588 
589  fprintf(ioQQQ,
590  " iso_level zone %.2f - largest residual in iso=%li %s nelem=%li matrix inversion is %g "
591  "level was %li Search?%c \n",
592  fnzone,
593  ipISO,
595  nelem ,
596  BigError ,
597  level_error,
598  TorF(conv.lgSearch) );
599  }
600 
601  // Force level balance to LTE
602  if ( iso_ctrl.lgLTE_levels[ipISO] )
603  {
604  t_iso_sp* sp = &iso_sp[ipISO][nelem];
605  double partfun=0.0;
606  for ( level = 0; level < numlevels_local; level++ )
607  {
608  partfun += sp->st[level].Boltzmann()*sp->st[level].g();
609  }
610  double scale = dense.xIonDense[nelem][nelem-ipISO]/partfun;
611  for ( level = 0; level < numlevels_local; level++ )
612  {
613  creation[level] = sp->st[level].Boltzmann()*sp->st[level].g()*scale;
614  }
615  }
616 
617  for( level = numlevels_local-1; level > 0; --level )
618  {
619  /* check for negative populations */
620  if( creation[level] < 0. )
621  lgNegPop = true;
622  }
623 
624  if( lgNegPop && dense.lgSetIoniz[nelem] )
625  {
626  // simulation can become unphysical if ionization is fixed.
627  // in this case, just put everything in ground.
628  // It's really the best we can do.
629  for( level = 1; level < numlevels_local; ++level )
630  creation[level] = 0.;
631  creation[0] = dense.xIonDense[nelem][nelem-ipISO];
632  // flip flag back
633  lgNegPop = false;
634  }
635 
636  /* put level populations into master array */
637  for( level=0; level < numlevels_local; level++ )
638  {
639  ASSERT( creation[level] >= 0. );
640  sp->st[level].Pop() = creation[level];
641 
642  if( sp->st[level].Pop() <= 0 && !conv.lgSearch )
643  {
644  fprintf(ioQQQ,
645  "PROBLEM non-positive level pop for iso = %li, nelem = "
646  "%li = %s, level=%li val=%.3e nTotalIoniz %li\n",
647  ipISO,
648  nelem ,
649  elementnames.chElementSym[nelem],
650  level,
651  sp->st[level].Pop() ,
652  conv.nTotalIoniz);
653  }
654  }
655 
656  /* zero populations of unused levels. */
657  for( level=numlevels_local; level < sp->numLevels_max; level++ )
658  {
659  sp->st[level].Pop() = 0.;
660  /* >>chng 06 jul 25, no need to zero this out, fix limit to 3-body heating elsewhere. */
661  /* sp->st[level].PopLTE = 0.; */
662  }
663 
664  /* TotalPopExcited is sum of excited level pops */
665  /* renormalize the populations to agree with ion solver */
666  iso_renorm( nelem, ipISO, renorm );
667 
668  double TotalPopExcited = 0.;
669  /* create sum of populations */
670  for( level=1; level < numlevels_local; level++ )
671  TotalPopExcited += sp->st[level].Pop();
672  ASSERT( TotalPopExcited >= 0. );
673  double TotalPop = TotalPopExcited + sp->st[0].Pop();
674 
675  /* option to force ionization */
676  if( dense.lgSetIoniz[nelem] )
677  {
678  if( !fp_equal( TotalPop, (double)dense.xIonDense[nelem][nelem-ipISO] ) )
679  {
680  if( TotalPopExcited >= dense.xIonDense[nelem][nelem-ipISO] )
681  {
682  for( level=0; level < numlevels_local; level++ )
683  sp->st[level].Pop() *=
684  dense.xIonDense[nelem][nelem-ipISO] / TotalPop;
685  }
686  else
687  {
688  sp->st[0].Pop() =
689  MAX2( 1e-30 * dense.xIonDense[nelem][nelem-ipISO],
690  dense.xIonDense[nelem][nelem-ipISO] - TotalPopExcited );
691  }
692  sp->lgPopsRescaled = true;
693  }
694  ASSERT( sp->st[0].Pop() >= 0. );
695  }
696  }
697  /* all solvers end up here */
698 
699  /* check on the sum of the populations */
700  if( lgNegPop || dense.xIonDense[nelem][nelem-ipISO] < 0. )
701  {
702  fprintf( ioQQQ,
703  " DISASTER iso_level finds negative ion fraction for iso=%2ld nelem=%2ld "
704  "%s using solver %s, IonFrac=%.3e simple=%.3e TE=%.3e ZONE=%4ld\n",
705  ipISO,
706  nelem,
707  elementnames.chElementSym[nelem],
708  sp->chTypeAtomUsed,
709  dense.xIonDense[nelem][nelem+1-ipISO]/SDIV(dense.xIonDense[nelem][nelem-ipISO]),
710  sp->xIonSimple,
711  phycon.te,
712  nzone );
713 
714  fprintf( ioQQQ, " level pop are:\n" );
715  for( i=0; i < numlevels_local; i++ )
716  {
717  fprintf( ioQQQ,PrintEfmt("%8.1e", sp->st[i].Pop() ));
718  if( (i!=0) && !(i%10) ) fprintf( ioQQQ,"\n" );
719  }
720  fprintf( ioQQQ, "\n" );
721  ContNegative();
722  ShowMe();
724  }
725 
726  for( ipHi=1; ipHi<numlevels_local; ++ipHi )
727  {
728  for( ipLo=0; ipLo<ipHi; ++ipLo )
729  {
730  if( sp->trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
731  continue;
732 
733  /* population of lower level, corrected for stimulated emission */
734  sp->trans(ipHi,ipLo).Emis().PopOpc() =
735  sp->st[ipLo].Pop() - sp->st[ipHi].Pop()*
736  sp->st[ipLo].g()/sp->st[ipHi].g();
737 
738  // don't allow masers from collapsed levels or in Case B
739  if( N_(ipHi) > sp->n_HighestResolved_local || opac.lgCaseB )
740  sp->trans(ipHi,ipLo).Emis().PopOpc() = MAX2( 0., sp->trans(ipHi,ipLo).Emis().PopOpc() );
741  }
742  }
743 
744  // Zero PopOpc of inactive transitions.
745  for( ipHi=numlevels_local; ipHi < sp->numLevels_max; ++ipHi )
746  {
747  for( ipLo=0; ipLo<ipHi; ++ipLo )
748  {
749  sp->trans(ipHi,ipLo).Emis().PopOpc() = 0.;
750  }
751  }
752  return;
753 }
754 
757 {
758  for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
759  {
760  if( ! dense.lgElmtOn[nelem] )
761  continue;
762  for ( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
763  {
764  if( (dense.IonHigh[nelem] >= nelem - ipISO) &&
765  (dense.IonLow[nelem] <= nelem - ipISO) )
766  {
767  iso_multiplet_opacities_one(ipISO, nelem);
768  }
769  }
770  }
771 }
772 
774  const long int ipISO, const long int nelem)
775 {
776  DEBUG_ENTRY( "iso_multiplet_opacities_one()" );
777  t_iso_sp* sp = &iso_sp[ipISO][nelem];
778  // Allocate array for calculating n-n' multiplet opacities.
779  multi_arr<double,2> MultOpacs;
780  long nMax = sp->n_HighestResolved_max + sp->nCollapsed_max;
781  MultOpacs.reserve( nMax+1 );
782  for( long n=2; n <= nMax; ++n )
783  MultOpacs.reserve( n, n+1 );
784  MultOpacs.alloc();
785  MultOpacs.zero();
786 
787  double rDopplerWidth = 1.0/GetDopplerWidth(dense.AtomicWeight[nelem]);
788 
789  // Sum n-n' multiplet opacities.
790  for( long ipHi=1; ipHi < sp->numLevels_max; ++ipHi )
791  {
792  for( long ipLo=0; ipLo < ipHi; ++ipLo )
793  {
794  const TransitionProxy& tr = sp->trans(ipHi,ipLo);
795  MultOpacs[ sp->st[ipHi].n() ][ sp->st[ipLo].n() ] += tr.Emis().PopOpc() *
796  tr.Emis().opacity() * rDopplerWidth;
797  }
798  }
799 
800  // Now store n-n' multiplet opacities.
801  for( long ipHi=1; ipHi < sp->numLevels_max; ++ipHi )
802  {
803  for( long ipLo=0; ipLo < ipHi; ++ipLo )
804  {
805  const TransitionProxy& tr = sp->trans(ipHi,ipLo);
806  tr.Emis().mult_opac() = MultOpacs[ sp->st[ipHi].n() ][ sp->st[ipLo].n() ];
807  }
808  }
809 }
810 
811 void iso_set_ion_rates( long ipISO, long nelem)
812 {
813  DEBUG_ENTRY( "iso_set_ion_rates()" );
814  t_iso_sp* sp = &iso_sp[ipISO][nelem];
815  const long int numlevels_local = sp->numLevels_local;
816  /* this is total ionization rate, s-1, of this species referenced to
817  * the total abundance */
818  double TotalPop = 0.;
819  ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] = 0.;
820  for( long level=0; level < numlevels_local; level++ )
821  {
822  /* sum of all ionization processes from this atom to ion, cm-3 s-1 now,
823  * but is divided by TotalPop below to become s-1 */
824  ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] +=
825  sp->st[level].Pop() * sp->fb[level].RateLevel2Cont;
826  TotalPop += sp->st[level].Pop();
827  }
828 
829  if( ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] > BIGDOUBLE )
830  {
831  fprintf( ioQQQ, "DISASTER RateIonizTot for Z=%li, ion %li is larger than BIGDOUBLE. This is a big problem.",
832  nelem+1, nelem-ipISO);
834  }
835 
836  if (TotalPop <= SMALLFLOAT)
837  ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] = sp->fb[0].RateLevel2Cont;
838  else
839  ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] /= TotalPop;
840 
841  if( ionbal.RateRecomIso[nelem][ipISO] > 0. )
842  sp->xIonSimple = ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1]/ionbal.RateRecomIso[nelem][ipISO];
843  else
844  sp->xIonSimple = 0.;
845 
846  ASSERT( ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] >= 0. );
847 
848  if( ipISO == ipHE_LIKE && nelem==ipHELIUM && nzone>0 )
849  {
850  /* find fraction of He0 destructions due to photoionization of 2^3S */
851  double ratio;
852  double rateOutOf2TripS = sp->st[ipHe2s3S].Pop() * sp->fb[ipHe2s3S].RateLevel2Cont;
853  if( rateOutOf2TripS > SMALLFLOAT )
854  {
855  ratio = rateOutOf2TripS /
856  ( rateOutOf2TripS + sp->st[ipHe1s1S].Pop()*ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] );
857  }
858  else
859  {
860  ratio = 0.;
861  }
862  if( ratio > he.frac_he0dest_23S )
863  {
864  /* remember zone where this happended and fraction, and frac due to photoionization */
865  he.nzone = nzone;
866  he.frac_he0dest_23S = ratio;
867  rateOutOf2TripS = sp->st[ipHe2s3S].Pop() *sp->fb[ipHe2s3S].gamnc;
868  if( rateOutOf2TripS > SMALLFLOAT )
869  {
870  he.frac_he0dest_23S_photo = rateOutOf2TripS /
871  ( rateOutOf2TripS + sp->st[ipHe1s1S].Pop()*ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] );
872  }
873  else
874  {
876  }
877  }
878  }
879 }
t_atmdat atmdat
Definition: atmdat.cpp:6
realnum & opacity() const
Definition: emission.h:638
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
qList st
Definition: iso.h:482
T * ptr0()
Definition: vectorize.h:221
const int ipHE_LIKE
Definition: iso.h:65
T * get_ptr(T *v)
Definition: cddefines.h:1113
double EdenHCorr
Definition: dense.h:227
const double BIGDOUBLE
Definition: cpu.h:249
t_opac opac
Definition: opacity.cpp:5
const realnum SMALLFLOAT
Definition: cpu.h:246
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
const int NISO
Definition: cddefines.h:310
bool lgIsoTraceFull[NISO]
Definition: trace.h:85
void prtRates(const long nlevels_local, const multi_arr< double, 2, C_TYPE > &a, valarray< double > &b)
Definition: prt.cpp:232
long int IonHigh[LIMELM+1]
Definition: dense.h:130
char TorF(bool l)
Definition: cddefines.h:753
const int ipHe2s3S
Definition: iso.h:46
#define PrintEfmt(F, V)
Definition: cddefines.h:1364
bool lgAdvection
Definition: dynamics.h:66
bool lgTimeDependentStatic
Definition: dynamics.h:102
long int nCollapsed_max
Definition: iso.h:518
t_conv conv
Definition: conv.cpp:5
t_phycon phycon
Definition: phycon.cpp:6
long int ipIsoTrace[NISO]
Definition: trace.h:88
bool lgSetIoniz[LIMELM]
Definition: dense.h:167
FILE * ioQQQ
Definition: cddefines.cpp:7
molezone * findspecieslocal(const char buf[])
long int nzone
Definition: cddefines.cpp:14
bool lgRandErrGen[NISO]
Definition: iso.h:430
const int ipHe1s1S
Definition: iso.h:43
bool lgLTE_levels[NISO]
Definition: iso.h:367
t_dynamics dynamics
Definition: dynamics.cpp:42
vector< freeBound > fb
Definition: iso.h:481
#define MIN2(a, b)
Definition: cddefines.h:807
bool lgISO[NISO]
Definition: dynamics.h:89
#define IPRAD
Definition: iso.h:88
t_dense dense
Definition: global.cpp:15
double ** RateRecomIso
Definition: ionbal.h:194
t_elementnames elementnames
Definition: elementnames.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
void vexp(const double x[], double y[], long nlo, long nhi)
realnum SmallA
Definition: iso.h:391
double ** RateRecomTot
Definition: ionbal.h:191
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
long int n_HighestResolved_local
Definition: iso.h:538
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
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:858
void iso_multiplet_opacities(void)
Definition: iso_level.cpp:756
double ** sink
Definition: mole.h:464
bool lgColl_excite[NISO]
Definition: iso.h:362
t_ionbal ionbal
Definition: ionbal.cpp:8
ColliderList colliders(dense)
vector< two_photon > TwoNu
Definition: iso.h:598
long int n_HighestResolved_max
Definition: iso.h:536
long nzone
Definition: he.h:29
long int IonLow[LIMELM+1]
Definition: dense.h:129
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
double * Boltzmann()
Definition: quantumstate.h:159
void ContNegative(void)
double energy(const genericState &gs)
const int ipH1s
Definition: iso.h:29
#define STATIC
Definition: cddefines.h:118
bool lgTrace
Definition: trace.h:12
double ** source
Definition: mole.h:464
double frac_he0dest_23S_photo
Definition: he.h:25
realnum Ploss() const
Definition: emission.h:130
EmissionList::reference Emis() const
Definition: transition.h:447
#define N_(A_)
Definition: iso.h:22
const char * chISO[NISO]
Definition: iso.h:348
t_mole_local mole
Definition: mole.cpp:8
realnum rate_lu_nontherm() const
Definition: collision.h:215
t_rfield rfield
Definition: rfield.cpp:9
#define IPCOLLIS
Definition: iso.h:89
bool lgCaseB
Definition: opacity.h:173
realnum *** xMoleChTrRate
Definition: mole.h:466
#define EXIT_FAILURE
Definition: cddefines.h:168
bool lgEquilibrium
Definition: dynamics.h:180
long max(int a, long b)
Definition: cddefines.h:821
#define cdEXIT(FAIL)
Definition: cddefines.h:484
realnum GetDopplerWidth(realnum massAMU)
t_iterations iterations
Definition: iterations.cpp:6
char chElementNameShort[LIMELM][CHARS_ELEMENT_NAME_SHORT]
Definition: elementnames.h:21
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:490
t_prt prt
Definition: prt.cpp:14
bool lgCritDensLMix[NISO]
Definition: iso.h:422
long int nTotalIoniz
Definition: conv.h:159
bool lgElmtOn[LIMELM]
Definition: dense.h:160
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
Definition: iso.h:470
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
int nLyaLevel[NISO]
Definition: iso.h:397
bool lgInd2nu_On
Definition: iso.h:375
bool lgPopsRescaled
Definition: iso.h:515
double *** CollIonRate_Ground
Definition: ionbal.h:118
realnum AtomicWeight[LIMELM]
Definition: dense.h:80
multi_arr< extra_tr, 2 > ex
Definition: iso.h:478
#define ASSERT(exp)
Definition: cddefines.h:617
bool lgLastIt
Definition: iterations.h:47
double qTot2S
Definition: iso.h:594
STATIC void iso_multiplet_opacities_one(const long int ipISO, const long int nelem)
Definition: iso_level.cpp:773
double Rate
Definition: dynamics.h:77
void reserve(size_type i1)
t_he he
Definition: he.cpp:5
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:307
t_prt_matrix matrix
Definition: prt.h:238
double *** RateIoniz
Definition: ionbal.h:177
CollisionProxy Coll() const
Definition: transition.h:463
double & mult_opac() const
Definition: emission.h:648
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
double & PopOpc() const
Definition: emission.h:658
iterator begin()
Definition: quantumstate.h:406
const int ipHELIUM
Definition: cddefines.h:349
double frac_he0dest_23S
Definition: he.h:25
double *** StatesElem
Definition: dynamics.h:83
double eden
Definition: dense.h:201
pointer data()
double EnergyRyd() const
Definition: transition.h:95
void iso_renorm(long nelem, long ipISO, double &renorm)
Definition: iso_solve.cpp:310
#define MAX2(a, b)
Definition: cddefines.h:828
double RateIonizTot(long nelem, long ion) const
Definition: ionbal.cpp:223
double xIonSimple
Definition: iso.h:496
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
bool lgCTOn
Definition: atmdat.h:325
double ** Source
Definition: dynamics.h:80
sys_float SDIV(sys_float x)
Definition: cddefines.h:1006
long int n_initial_relax
Definition: dynamics.h:132
long int numLevels_max
Definition: iso.h:524
char chElementName[LIMELM][CHARS_ELEMENT_NAME]
Definition: elementnames.h:17
#define fixit(a)
Definition: cddefines.h:416
GrainVar gv
Definition: grainvar.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
const int ipHYDROGEN
Definition: cddefines.h:348
realnum GrainChTrRate[LIMELM][LIMELM+1][LIMELM+1]
Definition: grainvar.h:543
realnum plsfrq
Definition: rfield.h:428
double CharExcRecTo[NCX][LIMELM][LIMELM+1]
Definition: atmdat.h:300
realnum & Aul() const
Definition: emission.h:668
bool lgPlasNu
Definition: rfield.h:426
void iso_set_ion_rates(long ipISO, long nelem)
Definition: iso_level.cpp:811
long int numLevels_local
Definition: iso.h:529
double ColUL(const ColliderList &colls) const
Definition: collision.h:106
void iso_level(const long ipISO, const long nelem, double &renorm, bool lgPrtMatrix)
char chTypeAtomUsed[10]
Definition: iso.h:587
double & pump() const
Definition: emission.h:518