cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ion_trim.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 /*ion_trim raise or lower most extreme stages of ionization considered,
4  * called by ConvBase - ion limits were originally set by */
5 #include "cddefines.h"
6 #include "ion_trim.h"
7 
8 #include "heavy.h"
9 #include "conv.h"
10 #include "rfield.h"
11 #include "phycon.h"
12 #include "mole.h"
13 #include "thermal.h"
14 #include "iso.h"
15 #include "struc.h"
16 #include "ionbal.h"
17 #include "dense.h"
18 
19 void ion_trim_untrim( long nelem )
20 {
21  DEBUG_ENTRY( "ion_trim_untrim()" );
22  dense.IonLow[nelem] = 0;
23  dense.IonHigh[nelem] = nelem+1;
24 }
25 void ion_trim_invalidate( long nelem )
26 {
27  DEBUG_ENTRY( "ion_trim_invalidate()" );
28  dense.IonLow[nelem] = -1;
29  dense.IonHigh[nelem] = -1;
30 }
31 
32 STATIC void ion_trim_from_set( long nelem);
33 
35 {
36  DEBUG_ENTRY( "ion_trim_init()" );
37  for(long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
38  {
39  if( !dense.lgElmtOn[nelem] )
40  {
41  ion_trim_invalidate(nelem);
42  continue;
43  }
44 
45  if( dense.lgSetIoniz[nelem] )
46  {
47  /* >>chng 04 jan 13, add this test, caught by Orly Gnat */
48  /* check on actual zero abundances of lower stages - this will only
49  * happen when ionization is set with element ionization command */
50  ion_trim_from_set(nelem);
51  }
52  else
53  {
54  // IonHigh[n] is the highest stage of ionization present
55  // the IonHigh array index is on the C scale, so [0] is hydrogen
56  // the value is also on the C scale, so element [nelem] can range
57  // from 0 to nelem+1
58  dense.IonHigh[nelem] = nelem + 1;
59 
60  dense.IonLow[nelem] = 0;
61  // for very intense radiation fields very heavy elements, N>Fe, will fail
62  // in ion_solver due to ill conditioned matrix. all populations are in
63  // fully stripped state. Start will fully stripped ion distribution in this case.
64  if( rfield.uh > 1e15 )
65  {
66  //trim down highest stage to be within incident radiation field
67  while ( ionbal.lgTrimhiOn &&
68  rfield.anu(Heavy.ipHeavy[nelem][dense.IonHigh[nelem]-1]) >
69  rfield.anu(rfield.nflux) && dense.IonHigh[nelem]>1 )
70  --dense.IonHigh[nelem];
71 
72  if (ionbal.lgTrimloOn)
73  dense.IonLow[nelem] = max( 0 , dense.IonHigh[nelem]-1 );
74  }
75  }
76 
77  // make low-stage populations zero
78  for( long ion=0; ion<dense.IonLow[nelem]; ++ion )
79  {
80  dense.xIonDense[nelem][dense.IonLow[nelem]] += dense.xIonDense[nelem][ion];
81  dense.xIonDense[nelem][ion] = 0.;
82  }
83  for( long ion=nelem+1; ion>dense.IonHigh[nelem]; --ion )
84  {
85  dense.xIonDense[nelem][dense.IonHigh[nelem]] += dense.xIonDense[nelem][ion];
86  dense.xIonDense[nelem][ion] = 0.;
87  }
88  }
89 }
90 
92 {
93  DEBUG_ENTRY( "ion_trim_from_set()" );
94  for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
95  {
96  if( dense.lgSetIoniz[nelem] )
97  {
98  ion_trim_from_set(nelem);
99  }
100  }
101 }
102 
103 STATIC void ion_trim_from_set( long nelem )
104 {
105  DEBUG_ENTRY( "ion_trim_from_set()" );
106  dense.IonLow[nelem] = 0;
107  dense.IonHigh[nelem] = nelem + 1;
108  while( ionbal.lgTrimloOn && dense.SetIoniz[nelem][dense.IonLow[nelem]] < dense.density_low_limit )
109  ++dense.IonLow[nelem];
110  while( ionbal.lgTrimhiOn && dense.SetIoniz[nelem][dense.IonHigh[nelem]] < dense.density_low_limit )
111  --dense.IonHigh[nelem];
112 }
113 
114 void ion_trim_small (long nelem, double abund_total)
115 {
116  DEBUG_ENTRY( "ion_trim_small()" );
117  while( ionbal.lgTrimhiOn && dense.IonHigh[nelem] > dense.IonLow[nelem] &&
118  dense.xIonDense[nelem][dense.IonHigh[nelem]] < 1e-25*abund_total &&
119  dense.IonHigh[nelem] > 1)
120  {
121  ASSERT( dense.xIonDense[nelem][dense.IonHigh[nelem]] >= 0. );
122  /* zero out abundance and heating due to stage of ionization we are about to zero out */
123  dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
124  thermal.setHeating(nelem,dense.IonHigh[nelem]-1,0.);
125  /* decrement counter */
126  --dense.IonHigh[nelem];
127  }
128 }
129 
130 void mole_ion_trim(void)
131 {
132  DEBUG_ENTRY( "mole_ion_trim()" );
133  for (ChemNuclideList::iterator atom = nuclide_list.begin();
134  atom != nuclide_list.end(); ++atom)
135  {
136  if( !(*atom)->lgHasLinkedIon())
137  continue;
138  const long int nelem=(*atom)->el->Z-1;
139  if( !dense.lgElmtOn[nelem] )
140  continue;
141 
142  for (long int ion=0;ion<nelem+2;ion++)
143  {
144  if ((*atom)->ipMl[ion] != -1)
145  {
146  if ( ionbal.lgTrimloOn && dense.IonLow[nelem] > ion )
147  {
148  if( dense.xIonDense[nelem][ion] > ( ionbal.trimlo * dense.gas_phase[nelem] ) )
149  {
150  dense.IonLow[nelem] = ion;
151  }
152  else
153  {
154  dense.xIonDense[nelem][dense.IonLow[nelem]] += dense.xIonDense[nelem][ion];
155  dense.xIonDense[nelem][ion] = 0.;
156  }
157  }
158 
159  if ( ionbal.lgTrimhiOn && dense.IonHigh[nelem] < ion )
160  {
161  if( dense.xIonDense[nelem][ion] > ( ionbal.trimhi * dense.gas_phase[nelem] ) )
162  {
163  dense.IonHigh[nelem] = ion;
164  }
165  else
166  {
167  dense.xIonDense[nelem][dense.IonHigh[nelem]] += dense.xIonDense[nelem][ion];
168  dense.xIonDense[nelem][ion] = 0.;
169  }
170  }
171  }
172  }
173  }
174 }
175 
176 void ion_trim(
177  /* nelem is on the C scale, 0 for H, 5 for C */
178  long int nelem )
179 {
180 
181  /* this will remember that higher stages trimed up or down */
182  bool lgHi_Down = false;
183  bool lgHi_Up = false;
184  bool lgHi_Up_enable;
185  /* this will remember that lower stages trimmed up or own*/
186  bool lgLo_Up = false;
187  bool lgLo_Down = false;
188  long int ion_lo_save = dense.IonLow[nelem],
189  ion_hi_save = dense.IonHigh[nelem];
190  long int ion;
191  realnum trimhi , trimlo;
192 
193  /* this is debugging code that turns on print after certain number of calls */
194  /*realnum xlimit = SMALLFLOAT *10.;*/
195  /*static int ncall=0;
196  if( nelem==5 )
197  ++ncall;*/
198 
199  DEBUG_ENTRY( "ion_trim()" );
200 
201  /* this routine should not be called early in the simulation, before
202  * things have settled down */
203  ASSERT( conv.nTotalIoniz>2 );
204 
205  /*confirm all ionization stages are within their range of validity */
206  ASSERT( nelem >= ipHELIUM && nelem < LIMELM );
207  ASSERT( dense.IonLow[nelem] >= 0 );
208  ASSERT( dense.IonHigh[nelem] <= nelem+1 );
209  /* IonHigh can be equal to IonLow */
210  ASSERT( dense.IonLow[nelem] <= dense.IonHigh[nelem] );
211 
212  /* during search phase of mostly neutral matter the electron density
213  * can be vastly too large, and the ionization suppressed as a result.
214  * during search do not trim down or up as much */
215  if( conv.lgSearch )
216  {
217  trimhi = (realnum)ionbal.trimhi * 1e-4f;
218  trimlo = (realnum)ionbal.trimlo * 1e-4f;
219  }
220  else
221  {
222  trimhi = (realnum)ionbal.trimhi;
223  trimlo = (realnum)ionbal.trimlo;
224  }
225 
226  /* helium is special case since abundance is so high, and He+ CT with
227  * CO is the dominant CO destruction process in molecular regions */
228  if( nelem == ipHELIUM )
229  {
230  /* never want to trip up a lower stage of ionization */
231  trimlo = SMALLFLOAT;
232 
233  /* if He+ is highest stage of ionization, probably want to keep it
234  * since important for CO chemistry in molecular clouds */
235  if( dense.IonHigh[ipHELIUM] == 1 )
236  {
237  trimhi = MIN2( trimhi , 1e-20f );
238  }
239  else if( dense.IonHigh[ipHELIUM] == 2 )
240  {
241  if( conv.lgSearch )
242  {
243  /* during search phase we may be quite far from solution, in the
244  * case of a PDR sim the electron density may be vastly higher
245  * than it will be with stable solution. He++ can be very important
246  * for the chemistry and we want to consider it. Make the
247  * threshold for ignoring He++ much higher than normal to prevent
248  * a premature removal of the ion */
249  trimhi = MIN2( trimhi , 1e-17f );
250  }
251  else
252  {
253  /* similar smaller upper limit for ion*/
254  trimhi = MIN2( trimhi , 1e-12f );
255  }
256  }
257  }
258 
259  /* logic for PDRs, for elements included in chemistry, need stable solutions,
260  * keep 3 ion stages in most cases - added by NA to do HII/PDR sims */
261  if( !mole_global.lgNoMole )
262  {
263  trimlo = SMALLFLOAT;
264  if(dense.IonHigh[nelem] ==2)
265  {
266  trimhi = MIN2(trimhi, 1e-20f);
267  }
268  }
269 
270  /* raise or lower highest and lowest stages of ionization to
271  * consider only significant stages
272  * IonHigh[nelem] stage of ionization */
273 
274  /* this is a special block for initialization only - it checks on absurd abundances
275  * and can trim multiple stages of ionization at one time. */
276  if( conv.lgSearch )
277  {
278  /* special - trim down higher stages if they have essentially zero abundance */
279  while( ionbal.lgTrimhiOn &&
280  (dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] < dense.density_low_limit ||
281  dense.xIonDense[nelem][dense.IonHigh[nelem]] < dense.density_low_limit ) &&
282  /* >>chng 02 may 12, rm +1 since this had effect of not allowing fully atomic */
283  dense.IonHigh[nelem] > dense.IonLow[nelem] )
284  {
285  /* dense.xIonDense[nelem][i] is density of i+1 th ionization stage (cm^-3)
286  * the -1 is correct for the heating, -1 since heating is caused by ionization of stage below it */
287  dense.xIonDense[nelem][dense.IonHigh[nelem]-1] +=
288  dense.xIonDense[nelem][dense.IonHigh[nelem]];
289  dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
290  thermal.setHeating(nelem,dense.IonHigh[nelem]-1,0.);
291  if( dense.IonHigh[nelem] == nelem+1-NISO )
292  {
293  long int ipISO = nelem - dense.IonHigh[nelem];
294  ASSERT( ipISO>=0 && ipISO<NISO );
295  for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
296  iso_sp[ipISO][nelem].st[level].Pop() = 0.;
297  }
298 
299  /* decrement high stage limit */
300  --dense.IonHigh[nelem];
301  ASSERT( dense.IonHigh[nelem] >= 0);
302  /* remember this happened */
303  lgHi_Down = true;
304  }
305 
306  /* special - trim up lower stages trim if they have essentially zero abundance */
307  while( ionbal.lgTrimloOn &&
308  (dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] < dense.density_low_limit ||
309  dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit ) &&
310  dense.IonLow[nelem] < dense.IonHigh[nelem] - 1 )
311  {
312  /* dense.xIonDense[nelem][i] is density of ith ionization stage (cm^-3)
313  * there is no-1 since we are removing the agent that heats */
314  dense.xIonDense[nelem][dense.IonLow[nelem]+1] +=
315  dense.xIonDense[nelem][dense.IonLow[nelem]];
316  dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.;
317  /* >>chng 01 aug 04, remove -1 which clobbers thermal.heating when IonLow == 0 */
318  thermal.setHeating(nelem,dense.IonLow[nelem],0.);
319  if( dense.IonLow[nelem] == nelem+1-NISO )
320  {
321  long int ipISO = nelem - dense.IonLow[nelem];
322  ASSERT( ipISO>=0 && ipISO<NISO );
323  for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
324  iso_sp[ipISO][nelem].st[level].Pop() = 0.;
325  }
326 
327  /* increment low stage limit */
328  ++dense.IonLow[nelem];
329  lgLo_Up = true;
330  }
331  }
332 
333  /* sanity check */
334  ASSERT( dense.IonLow[nelem] <= dense.IonHigh[nelem] );
335 
336  /* this checks on error condition where the gas is stupendously highly ionized,
337  * the low limit is already one less than high, and we need to raise the low
338  * limit further */
339  if( dense.IonHigh[nelem] > 1 &&
340  (dense.IonLow[nelem]==dense.IonHigh[nelem]-1) &&
341  (dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit) )
342  {
343 # if 0
344  fprintf(ioQQQ,
345  "PROBLEM DISASTER the density of ion %li of element %s is too small to be computed on this cpu.\n",
346  dense.IonLow[nelem]+1,
347  elementnames.chElementName[nelem]);
348  fprintf(ioQQQ,
349  "Turn off the element with the command ELEMENT XXX OFF or do not consider "
350  "gas with low density, the current hydrogen density is %.2e cm-3.\n",
352  fprintf(ioQQQ,
353  "The outer radius of the cloud is %.2e cm - should a stopping "
354  "radius be set?\n",
355  radius.Radius );
356  fprintf(ioQQQ,
357  "The abort flag is being set - the calculation is stopping.\n");
358 # endif
359  //lgAbort = true;
361  return;
362  }
363 
364  /* trim down high stages that have too small an abundance */
365  while( ionbal.lgTrimhiOn &&
366  dense.IonHigh[nelem] > dense.IonLow[nelem] &&
367  ( (dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] <=
368  trimhi ) ||
369  (dense.xIonDense[nelem][dense.IonHigh[nelem]] <= dense.density_low_limit ) )
370  )
371  {
372  /* >>chng 03 sep 30, the atom and its first ion are a special case
373  * since we want to compute even trivial ions in molecular clouds */
374  if( dense.IonHigh[nelem]>1 ||
375  (dense.IonHigh[nelem]==1&&dense.xIonDense[nelem][1]<100.*dense.density_low_limit) )
376  {
377  dense.xIonDense[nelem][dense.IonHigh[nelem]-1] +=
378  dense.xIonDense[nelem][dense.IonHigh[nelem]];
379  dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
380  thermal.setHeating(nelem,dense.IonHigh[nelem]-1,0.);
381  if( dense.IonHigh[nelem] == nelem+1-NISO )
382  {
383  long int ipISO = nelem - dense.IonHigh[nelem];
384  ASSERT( ipISO>=0 && ipISO<NISO );
385  for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
386  iso_sp[ipISO][nelem].st[level].Pop() = 0.;
387  }
388  --dense.IonHigh[nelem];
389  lgHi_Down = true;
390  }
391  else
392  {
393  break;
394  }
395  }
396 
397  /* trim up highest stages - will this be done? */
398  lgHi_Up_enable = true;
399  /* trimming can be turned off */
400  if( !ionbal.lgTrimhiOn )
401  lgHi_Up_enable = false;
402  /* high stage is already fully stripped */
403  if( dense.IonHigh[nelem]>=nelem+1 )
404  lgHi_Up_enable = false;
405  /* we have previously trimmed this stage down */
406  if( lgHi_Down )
407  lgHi_Up_enable = false;
408  /* we are in neutral gas */
410  lgHi_Up_enable = false;
411 
412  if( lgHi_Up_enable )
413  {
414  realnum abundold = 0;
415  if( nzone>2 )
416  abundold = struc.xIonDense[nzone-3][nelem][dense.IonHigh[nelem]]/
417  SDIV( struc.gas_phase[nzone-3][nelem]);
418  realnum abundnew = dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem];
419  /* only raise highest stage if ionization potential of next highest stage is within
420  * continuum array and the abundance of the highest stage is significant */
421  if( (Heavy.Valence_IP_Ryd[nelem][dense.IonHigh[nelem]] < rfield.anu(rfield.nflux-1) ||
422  Heavy.Valence_IP_Ryd[nelem][dense.IonHigh[nelem]] < phycon.te_ryd*10.) &&
423  dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] > trimhi &&
424  /* this checks that abundance of highest stage is increasing */
425  abundnew > abundold*1.01 )
426  {
427  /*fprintf(ioQQQ,"uuppp %li %li \n", nelem, dense.IonHigh[nelem] );*/
428  /* raise highest level of ionization */
429  ++dense.IonHigh[nelem];
430  lgHi_Up = true;
431  /* set abundance to almost zero so that sanity check elsewhere will be ok */
432  dense.xIonDense[nelem][dense.IonHigh[nelem]] = (dense.density_low_limit*10.);
433  dense.xIonDense[nelem][dense.IonHigh[nelem]-1] -=
434  dense.xIonDense[nelem][dense.IonHigh[nelem]];
435  long int ipISO = nelem - dense.IonHigh[nelem];
436  if (ipISO >= 0 && ipISO < NISO)
437  iso_sp[ipISO][nelem].st[0].Pop() = dense.xIonDense[nelem][dense.IonHigh[nelem]];
438  }
439  }
440 
441  /* sanity check */
442  ASSERT( dense.IonLow[nelem] <= dense.IonHigh[nelem] );
443 
444  /* lower lowest stage of ionization if we have significant abundance at current lowest */
445  if( dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] > trimlo &&
446  dense.IonLow[nelem] > 0 )
447  {
448  /* lower lowest level of ionization */
449  --dense.IonLow[nelem];
450  lgLo_Down = true;
451  /* >>chng 01 aug 02, set this to zero so that sanity check elsewhere will be ok */
453  dense.xIonDense[nelem][dense.IonLow[nelem]+1] -=
454  dense.xIonDense[nelem][dense.IonLow[nelem]];
455  long int ipISO = nelem - dense.IonLow[nelem];
456  if (ipISO >= 0 && ipISO < NISO)
457  iso_sp[ipISO][nelem].st[0].Pop() = dense.xIonDense[nelem][dense.IonLow[nelem]];
458  }
459 
460  /* raise lowest stage of ionization, but only if we are near illuminated face of cloud*/
461  else if( ionbal.lgTrimloOn && nzone < 10 &&
462  (dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] <= (realnum)trimlo) &&
463  (dense.IonLow[nelem] < (dense.IonHigh[nelem] - 2) ) )
464  {
465  /* raise lowest level of ionization */
466  dense.xIonDense[nelem][dense.IonLow[nelem]+1] +=
467  dense.xIonDense[nelem][dense.IonLow[nelem]];
468  dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.;
469  /* no minus one in below since this is low bound, already bounds at atom */
470  thermal.setHeating(nelem,dense.IonLow[nelem],0.);
471  if( dense.IonLow[nelem] == nelem+1-NISO )
472  {
473  long int ipISO = nelem - dense.IonLow[nelem];
474  ASSERT( ipISO>=0 && ipISO<NISO );
475  for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
476  iso_sp[ipISO][nelem].st[level].Pop() = 0.;
477  }
478  ++dense.IonLow[nelem];
479  lgLo_Up = true;
480  }
481  /* test on zero */
482  else if( ionbal.lgTrimloOn && (dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit) &&
483  /*>>chng 06 may 24, from IonLow < IonHigh to <IonHigh-1 because
484  * old test would allow low to be raised too high, which then
485  * leads to insanity */
486  (dense.IonLow[nelem] < dense.IonHigh[nelem]-1) )
487  {
488  while(dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit &&
489  /* >>chng 07 feb 27 from < IonHigh to < IonHigh-1 to prevent raising
490  * low to high */
491  (dense.IonLow[nelem] < dense.IonHigh[nelem]-1) )
492  {
493  /* raise lowest level of ionization */
494  dense.xIonDense[nelem][dense.IonLow[nelem]+1] +=
495  dense.xIonDense[nelem][dense.IonLow[nelem]];
496  dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.;
497  /* no minus one in below since this is low bound, already bounds at atom */
498  thermal.setHeating(nelem,dense.IonLow[nelem],0.);
499  if( dense.IonLow[nelem] == nelem+1-NISO )
500  {
501  long int ipISO = nelem - dense.IonLow[nelem];
502  ASSERT( ipISO>=0 && ipISO<NISO );
503  for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
504  iso_sp[ipISO][nelem].st[level].Pop() = 0.;
505  }
506  ++dense.IonLow[nelem];
507  lgLo_Up = true;
508  }
509  }
510 
511  /* sanity check */
512  ASSERT( dense.IonLow[nelem] <= dense.IonHigh[nelem] );
513 
514  /* these are standard bounds checks that appear throughout this routine
515  * dense.xIonDense[IonLow] is first one with positive abundances
516  *
517  * in case where lower ionization stage was just lowered the abundance
518  * was set to dense.density_low_limit so test must be < dense.density_low_limit */
519  ASSERT( dense.IonLow[nelem] <= 1 ||
520  dense.xIonDense[nelem][dense.IonLow[nelem]-1] == 0. );
521 
522  ASSERT( (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) || lgLo_Up ||
523  ! ionbal.lgTrimloOn ||
524  dense.xIonDense[nelem][dense.IonLow[nelem]] >= dense.density_low_limit ||
525  dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] >= dense.density_low_limit ||
526  /*>>chng 06 may 24, include this to cover case where cap is set by not allowing
527  * lower limit to come up to upper limit */
528  (dense.IonLow[nelem]==dense.IonHigh[nelem]-1 ));
529 
530  {
531  /* option to print out what has happened so far .... */
532  enum {DEBUG_LOC=false};
533  if( DEBUG_LOC && nelem == ipHELIUM )
534  {
535  if( lgHi_Down ||lgHi_Up ||lgLo_Up ||lgLo_Down )
536  {
537  fprintf(ioQQQ,"DEBUG TrimZone\t%li\t",nzone );
538  if( lgHi_Down )
539  {
540  fprintf(ioQQQ,"high dn %li to %li",
541  ion_hi_save ,
542  dense.IonHigh[nelem] );
543  }
544  if( lgHi_Up )
545  {
546  fprintf(ioQQQ,"high up %li to %li",
547  ion_hi_save ,
548  dense.IonHigh[nelem] );
549  }
550  if( lgLo_Up )
551  {
552  fprintf(ioQQQ,"low up %li to %li",
553  ion_lo_save ,
554  dense.IonLow[nelem] );
555  }
556  if( lgLo_Down )
557  {
558  fprintf(ioQQQ,"low dn %li to %li",
559  ion_lo_save ,
560  dense.IonLow[nelem] );
561  }
562  fprintf(ioQQQ,"\n" );
563  }
564  }
565  }
566 
567  /* option to print if any trimming occurred */
568  if( lgHi_Down || lgHi_Up || lgLo_Up || lgLo_Down )
569  {
570  conv.lgIonStageTrimed = true;
571  {
572  /* option to print out what has happened so far .... */
573  enum {DEBUG_LOC=false};
574  if( DEBUG_LOC && nelem==ipHELIUM )
575  {
576  fprintf(ioQQQ,"DEBUG ion_trim zone\t%.2f \t%li\t", fnzone, nelem);
577  if( lgHi_Down )
578  fprintf(ioQQQ,"\tHi_Down");
579  if( lgHi_Up )
580  fprintf(ioQQQ,"\tHi_Up");
581  if( lgLo_Up )
582  fprintf(ioQQQ,"\tLo_Up");
583  if( lgLo_Down )
584  fprintf(ioQQQ,"\tLo_Down");
585  fprintf(ioQQQ,"\n");
586  }
587  }
588  }
589 
590  /* asserts that that densities of ion stages that are now turned
591  * off are zero */
592  for( ion=0; ion<dense.IonLow[nelem]; ++ion )
593  {
594  ASSERT( dense.xIonDense[nelem][ion] == 0. );
595  }
596  for( ion=dense.IonHigh[nelem]+1; ion<nelem+1; ++ion )
597  {
598  ASSERT( dense.xIonDense[nelem][ion] == 0. );
599  }
600 
601  for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
602  {
603  /* in case where ionization stage was changed the
604  * density is zero, we set it to SMALLFLOAT since a density
605  * of zero is not possible */
606  dense.xIonDense[nelem][ion] = MAX2(dense.xIonDense[nelem][ion], SMALLFLOAT);
607  }
608  return;
609 }
610 
611 // Try winnowing down ion_trim
613  /* nelem is on the C scale, 0 for H, 5 for C */
614  long int nelem )
615 {
616 
617  /* this will remember that higher stages trimed up or down */
618  bool lgHi_Down = false;
619  /* this will remember that lower stages trimmed up or own*/
620  bool lgLo_Up = false;
621  long int ion_lo_save = dense.IonLow[nelem],
622  ion_hi_save = dense.IonHigh[nelem];
623  long int ion;
624  realnum trimhi , trimlo, trimcharge;
625 
626  /* this is debugging code that turns on print after certain number of calls */
627  /*realnum xlimit = SMALLFLOAT *10.;*/
628  /*static int ncall=0;
629  if( nelem==5 )
630  ++ncall;*/
631 
632  DEBUG_ENTRY( "ion_trim2()" );
633 
634  /*confirm all ionization stages are within their range of validity */
635  ASSERT( nelem >= ipHYDROGEN && nelem < LIMELM );
636  ASSERT( dense.IonLow[nelem] >= 0 );
637  ASSERT( dense.IonHigh[nelem] <= nelem+1 );
638  /* IonHigh can be equal to IonLow */
639  ASSERT( dense.IonLow[nelem] <= dense.IonHigh[nelem] );
640 
641  // Trimming limit on fraction of total charge in system
642  trimcharge = 1e-6;
643 
644  /* during search phase of mostly neutral matter the electron density
645  * can be vastly too large, and the ionization suppressed as a result.
646  * during search do not trim down or up as much */
647  if( conv.lgSearch )
648  {
649  trimhi = (realnum)ionbal.trimhi * 1e-4f;
650  trimlo = (realnum)ionbal.trimlo * 1e-4f;
651  }
652  else
653  {
654  trimhi = (realnum)ionbal.trimhi;
655  trimlo = (realnum)ionbal.trimlo;
656  }
657 
658  /* helium is special case since abundance is so high, and He+ CT with
659  * CO is the dominant CO destruction process in molecular regions */
660  if( nelem == ipHELIUM )
661  {
662  /* never want to trip up a lower stage of ionization */
663  trimlo = SMALLFLOAT;
664 
665  /* if He+ is highest stage of ionization, probably want to keep it
666  * since important for CO chemistry in molecular clouds */
667  if( dense.IonHigh[ipHELIUM] == 1 )
668  {
669  trimhi = MIN2( trimhi , 1e-20f );
670  }
671  else if( dense.IonHigh[ipHELIUM] == 2 )
672  {
673  if( conv.lgSearch )
674  {
675  /* during search phase we may be quite far from solution, in the
676  * case of a PDR sim the electron density may be vastly higher
677  * than it will be with stable solution. He++ can be very important
678  * for the chemistry and we want to consider it. Make the
679  * threshold for ignoring He++ much higher than normal to prevent
680  * a premature removal of the ion */
681  trimhi = MIN2( trimhi , 1e-17f );
682  }
683  else
684  {
685  /* similar smaller upper limit for ion*/
686  trimhi = MIN2( trimhi , 1e-12f );
687  }
688  }
689  }
690 
691  /* logic for PDRs, for elements included in chemistry, need stable solutions,
692  * keep 3 ion stages in most cases - added by NA to do HII/PDR sims */
693  if( !mole_global.lgNoMole )
694  {
695  trimlo = SMALLFLOAT;
696  if(dense.IonHigh[nelem] ==2)
697  {
698  trimhi = MIN2(trimhi, 1e-20f);
699  }
700  }
701 
702  /* raise or lower highest and lowest stages of ionization to
703  * consider only significant stages
704  * IonHigh[nelem] stage of ionization */
705 
706  /* sanity check */
707  ASSERT( dense.IonLow[nelem] <= dense.IonHigh[nelem] );
708 
709  /* this checks on error condition where the gas is stupendously highly ionized,
710  * the low limit is already one less than high, and we need to raise the low
711  * limit further */
712  if( dense.IonHigh[nelem] > 1 &&
713  (dense.IonLow[nelem]==dense.IonHigh[nelem]-1) &&
714  (dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit) )
715  {
716 # if 0
717  fprintf(ioQQQ,
718  "PROBLEM DISASTER the density of ion %li of element %s is too small to be computed on this cpu.\n",
719  dense.IonLow[nelem]+1,
720  elementnames.chElementName[nelem]);
721  fprintf(ioQQQ,
722  "Turn off the element with the command ELEMENT XXX OFF or do not consider "
723  "gas with low density, the current hydrogen density is %.2e cm-3.\n",
725  fprintf(ioQQQ,
726  "The outer radius of the cloud is %.2e cm - should a stopping "
727  "radius be set?\n",
728  radius.Radius );
729  fprintf(ioQQQ,
730  "The abort flag is being set - the calculation is stopping.\n");
731 # endif
732  //lgAbort = true;
734  return;
735  }
736 
737  double itrim = 1.0;
738  /* trim down high stages that have too small an abundance */
739  while ( ionbal.lgTrimhiOn &&
740  dense.IonHigh[nelem] > dense.IonLow[nelem] &&
741  ( ( dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] <= trimhi*itrim &&
742  dense.xIonDense[nelem][dense.IonHigh[nelem]] <= trimcharge*dense.EdenTrue )
743  ||
744  dense.xIonDense[nelem][dense.IonHigh[nelem]] <= dense.density_low_limit )
745  )
746  {
747  itrim = 0.125;
748  dense.xIonDense[nelem][dense.IonHigh[nelem]-1] +=
749  dense.xIonDense[nelem][dense.IonHigh[nelem]];
750  dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
751  thermal.setHeating(nelem,dense.IonHigh[nelem]-1,0.);
752  long int ipISO = nelem - dense.IonHigh[nelem];
753  if ( ipISO>=0 && ipISO<NISO )
754  {
755  for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
756  iso_sp[ipISO][nelem].st[level].Pop() = 0.;
757  }
758  --dense.IonHigh[nelem];
759  lgHi_Down = true;
760  }
761 
762  /* sanity check */
763  ASSERT( dense.IonLow[nelem] <= dense.IonHigh[nelem] );
764 
765  itrim = 1.0;
766  /* raise lowest stage of ionization, but only if we are near illuminated face of cloud*/
767  while( ionbal.lgTrimloOn &&
768  (dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] <= (realnum)trimlo*itrim) &&
769  (dense.IonLow[nelem] < (dense.IonHigh[nelem] - 2) ) )
770  {
771  itrim = 0.125;
772  /* raise lowest level of ionization */
773  dense.xIonDense[nelem][dense.IonLow[nelem]+1] +=
774  dense.xIonDense[nelem][dense.IonLow[nelem]];
775  dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.;
776  /* no minus one in below since this is low bound, already bounds at atom */
777  thermal.setHeating(nelem,dense.IonLow[nelem],0.);
778  long int ipISO = nelem - dense.IonLow[nelem];
779  if ( ipISO>=0 && ipISO<NISO )
780  {
781  for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
782  iso_sp[ipISO][nelem].st[level].Pop() = 0.;
783  }
784  ++dense.IonLow[nelem];
785  lgLo_Up = true;
786  }
787 
788  /* sanity check */
789  ASSERT( dense.IonLow[nelem] <= dense.IonHigh[nelem] );
790 
791  /* these are standard bounds checks that appear throughout this routine
792  * dense.xIonDense[IonLow] is first one with positive abundances
793  *
794  * in case where lower ionization stage was just lowered the abundance
795  * was set to dense.density_low_limit so test must be < dense.density_low_limit */
796  ASSERT( dense.IonLow[nelem] <= 1 ||
797  dense.xIonDense[nelem][dense.IonLow[nelem]-1] == 0. );
798 
799  {
800  /* option to print out what has happened so far .... */
801  enum {DEBUG_LOC=false};
802  if( DEBUG_LOC && nelem == ipHELIUM )
803  {
804  if( lgHi_Down ||lgLo_Up )
805  {
806  fprintf(ioQQQ,"DEBUG TrimZone\t%li\t",nzone );
807  if( lgHi_Down )
808  {
809  fprintf(ioQQQ,"high dn %li to %li",
810  ion_hi_save ,
811  dense.IonHigh[nelem] );
812  }
813  if( lgLo_Up )
814  {
815  fprintf(ioQQQ,"low up %li to %li",
816  ion_lo_save ,
817  dense.IonLow[nelem] );
818  }
819  fprintf(ioQQQ,"\n" );
820  }
821  }
822  }
823 
824  /* option to print if any trimming occurred */
825  if( lgHi_Down || lgLo_Up )
826  {
827  conv.lgIonStageTrimed = true;
828  {
829  /* option to print out what has happened so far .... */
830  enum {DEBUG_LOC=false};
831  if( DEBUG_LOC && nelem==ipHELIUM )
832  {
833  fprintf(ioQQQ,"DEBUG ion_trim zone\t%.2f \t%li\t", fnzone, nelem);
834  if( lgHi_Down )
835  fprintf(ioQQQ,"\tHi_Down");
836  if( lgLo_Up )
837  fprintf(ioQQQ,"\tLo_Up");
838  fprintf(ioQQQ,"\n");
839  }
840  }
841  }
842 
843  /* asserts that that densities of ion stages that are now turned
844  * off are zero */
845  for( ion=0; ion<dense.IonLow[nelem]; ++ion )
846  {
847  ASSERT( dense.xIonDense[nelem][ion] == 0. );
848  }
849  for( ion=dense.IonHigh[nelem]+1; ion<nelem+1; ++ion )
850  {
851  ASSERT( dense.xIonDense[nelem][ion] == 0. );
852  }
853 
854  for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
855  {
856  /* in case where ionization stage was changed the
857  * density is zero, we set it to SMALLFLOAT since a density
858  * of zero is not possible */
859  dense.xIonDense[nelem][ion] = MAX2(dense.xIonDense[nelem][ion], SMALLFLOAT);
860  }
861  return;
862 }
863 
864 void ion_widen(long nelem)
865 {
866  DEBUG_ENTRY( "ion_widen()" );
867  if ( dense.lgSetIoniz[nelem] ||
868  ( dense.IonLow[nelem] == 0 && dense.IonHigh[nelem] == nelem+1 ) )
869  {
870  return;
871  }
872  double abund = 0.0;
873  for (long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion)
874  {
875  abund += dense.xIonDense[nelem][ion];
876  }
877  double abundnew = abund;
878  if (dense.IonLow[nelem] > 0 )
879  {
880  --dense.IonLow[nelem];
881  dense.xIonDense[nelem][dense.IonLow[nelem]] = (dense.density_low_limit*10.);
882  abundnew += (dense.density_low_limit*10.);
883  long ipISO = nelem - dense.IonLow[nelem];
884  if (ipISO >= 0 && ipISO < NISO)
885  iso_sp[ipISO][nelem].st[0].Pop() = dense.xIonDense[nelem][dense.IonLow[nelem]];
886  }
887  if ( dense.IonHigh[nelem] < nelem+1)
888  {
889  ++dense.IonHigh[nelem];
890  dense.xIonDense[nelem][dense.IonHigh[nelem]] = (dense.density_low_limit*10.);
891  abundnew += (dense.density_low_limit*10.);
892  long ipISO = nelem - dense.IonHigh[nelem];
893  if (ipISO >= 0 && ipISO < NISO)
894  iso_sp[ipISO][nelem].st[0].Pop() = dense.xIonDense[nelem][dense.IonHigh[nelem]];
895  }
896  double frac = abund/abundnew;
897  for (long ion=dense.IonLow[nelem];ion<=dense.IonHigh[nelem];++ion)
898  {
899  dense.xIonDense[nelem][ion] *= frac;
900  }
901 }
902 
903 #ifdef NDEBUG
904 void ion_trim_validate(long, bool)
905 {
906  (void)0;
907 }
908 #else
909 void ion_trim_validate(long nelem, bool lgIonizTrimCalled)
910 {
911  DEBUG_ENTRY( "ion_trim_validate()" );
912  for( long ion=0; ion<dense.IonLow[nelem]; ++ion )
913  {
914  ASSERT( dense.xIonDense[nelem][ion] == 0. );
915  }
916  /*if( nelem==5 ) fprintf(ioQQQ,"carbbb\t%li\n", dense.IonHigh[nelem]);*/
917  for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
918  {
919  /* >>chng 02 feb 06, had been > o., chng to > SMALLFLOAT to
920  * trip over VERY small floats that failed on alphas, but not 386
921  *
922  * in case where lower ionization stage was just lowered or
923  * trimmed down the abundance
924  * was set to SMALLFLOAT so test must be < SMALLFLOAT */
925  /* >>chng 02 feb 19, add check for search phase. During this search
926  * models with extreme ionization (all neutral or all ionized) can
927  * have extreme but non-zero abundances far from the ionization peak for
928  * element with lots of electrons. These will go away once the model
929  * becomes stable */
930  /* >>chng 03 dec 01, add check on whether ion trim was called
931  * conserve.in threw assert when iontrim not called and abund grew small */
932  ASSERT( conv.lgSearch || !lgIonizTrimCalled ||
933  /* this can happen if all C is in the form of CO --
934  generalize test to being that at least one ion must survive */
935  dense.IonLow[nelem] == dense.IonHigh[nelem] ||
936  dense.xIonDense[nelem][ion] >= SMALLFLOAT ||
937  dense.xIonDense[nelem][ion]/dense.gas_phase[nelem] >= SMALLFLOAT );
938  }
939  for( long ion=dense.IonHigh[nelem]+1; ion<=nelem+1; ++ion )
940  {
941  ASSERT( ion >= 0 );
942  ASSERT( dense.xIonDense[nelem][ion] == 0. );
943  }
944 }
945 #endif
void ion_trim(long int nelem)
Definition: ion_trim.cpp:176
t_mole_global mole_global
Definition: mole.cpp:7
double Radius
Definition: radius.h:31
bool lgTrimhiOn
Definition: ionbal.h:90
t_thermal thermal
Definition: thermal.cpp:6
qList st
Definition: iso.h:482
t_struc struc
Definition: struc.cpp:6
t_Heavy Heavy
Definition: heavy.cpp:5
const realnum SMALLFLOAT
Definition: cpu.h:246
realnum SetIoniz[LIMELM][LIMELM+1]
Definition: dense.h:172
const int NISO
Definition: cddefines.h:310
void ion_widen(long nelem)
Definition: ion_trim.cpp:864
long int IonHigh[LIMELM+1]
Definition: dense.h:130
void ion_trim_untrim(long nelem)
Definition: ion_trim.cpp:19
void ion_trim_invalidate(long nelem)
Definition: ion_trim.cpp:25
t_conv conv
Definition: conv.cpp:5
void mole_ion_trim(void)
Definition: ion_trim.cpp:130
t_phycon phycon
Definition: phycon.cpp:6
double trimhi
Definition: ionbal.h:83
void ion_trim_init()
Definition: ion_trim.cpp:34
bool lgSetIoniz[LIMELM]
Definition: dense.h:167
FILE * ioQQQ
Definition: cddefines.cpp:7
long int nzone
Definition: cddefines.cpp:14
ChemNuclideList nuclide_list
#define MIN2(a, b)
Definition: cddefines.h:807
double anu(size_t i) const
Definition: mesh.h:111
t_dense dense
Definition: global.cpp:15
t_elementnames elementnames
Definition: elementnames.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
bool lgIonStageTrimed
Definition: conv.h:184
STATIC void ion_trim_from_set(long nelem)
Definition: ion_trim.cpp:103
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
bool lgSearch
Definition: conv.h:168
t_ionbal ionbal
Definition: ionbal.cpp:8
t_abund abund
Definition: abund.cpp:5
long int IonLow[LIMELM+1]
Definition: dense.h:129
void ion_trim2(long int nelem)
Definition: ion_trim.cpp:612
#define STATIC
Definition: cddefines.h:118
double EdenTrue
Definition: dense.h:232
t_rfield rfield
Definition: rfield.cpp:9
float realnum
Definition: cddefines.h:124
realnum uh
Definition: rfield.h:347
long max(int a, long b)
Definition: cddefines.h:821
t_radius radius
Definition: radius.cpp:5
long int nTotalIoniz
Definition: conv.h:159
bool lgElmtOn[LIMELM]
Definition: dense.h:160
void setHeating(long nelem, long ion, double heating)
Definition: thermal.h:190
realnum gas_phase[LIMELM]
Definition: dense.h:76
void ion_trim_small(long nelem, double abund_total)
Definition: ion_trim.cpp:114
#define ASSERT(exp)
Definition: cddefines.h:617
const int LIMELM
Definition: cddefines.h:307
double Valence_IP_Ryd[LIMELM][LIMELM]
Definition: heavy.h:24
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
const int ipHELIUM
Definition: cddefines.h:349
double te_ryd
Definition: phycon.h:27
realnum ** gas_phase
Definition: struc.h:75
#define MAX2(a, b)
Definition: cddefines.h:828
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
sys_float SDIV(sys_float x)
Definition: cddefines.h:1006
double trimlo
Definition: ionbal.h:83
long int numLevels_max
Definition: iso.h:524
char chElementName[LIMELM][CHARS_ELEMENT_NAME]
Definition: elementnames.h:17
bool lgTrimloOn
Definition: ionbal.h:93
double fnzone
Definition: cddefines.cpp:15
double frac(double d)
bool lgNoMole
Definition: mole.h:325
void ion_trim_validate(long nelem, bool lgIonizTrimCalled)
Definition: ion_trim.cpp:909
const int ipHYDROGEN
Definition: cddefines.h:348
long int nflux
Definition: rfield.h:48
realnum *** xIonDense
Definition: struc.h:64
long int ipHeavy[LIMELM][LIMELM]
Definition: heavy.h:11
double density_low_limit
Definition: dense.h:208