cloudy  trunk
 All Data Structures 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-2008 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 "elementnames.h"
7 #include "radius.h"
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 "dense.h"
16 #include "conv.h"
17 #include "struc.h"
18 #include "ionbal.h"
19 
20 void ion_trim(
21  /* nelem is on the C scale, 0 for H, 5 for C */
22  long int nelem )
23 {
24 
25  /* this will remember that higher stages trimed up or down */
26  bool lgHi_Down = false;
27  bool lgHi_Up = false;
28  bool lgHi_Up_enable;
29  /* this will remember that lower stages trimmed up or own*/
30  bool lgLo_Up = false;
31  bool lgLo_Down = false;
32  long int ion_lo_save = dense.IonLow[nelem],
33  ion_hi_save = dense.IonHigh[nelem];
34  long int ion;
35  realnum trimhi , trimlo;
36 
37  /* this is debugging code that turns on print after certain number of calls */
38  /*realnum xlimit = SMALLFLOAT *10.;*/
39  /*static int ncall=0;
40  if( nelem==5 )
41  ++ncall;*/
42 
43  DEBUG_ENTRY( "ion_trim()" );
44 
45  /* this routine should not be called early in the simulation, before
46  * things have settled down */
47  ASSERT( conv.nTotalIoniz>2 );
48 
49  /*confirm all ionization stages are within their range of validity */
50  ASSERT( nelem >= ipHELIUM && nelem < LIMELM );
51  ASSERT( dense.IonLow[nelem] >= 0 );
52  ASSERT( dense.IonHigh[nelem] <= nelem+1 );
53  /* IonHigh can be equal to IonLow if both are zero - so no ionization,
54  * or is ionization has been set with "element ionization" command */
55  ASSERT( dense.IonLow[nelem] < dense.IonHigh[nelem] ||
56  (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) ||
57  dense.lgSetIoniz[nelem] );
58 
59  /* during search phase of mostly neutral matter the electron density
60  * can be vastly too large, and the ionization suppressed as a result.
61  * during search do not trim down or up as much */
62  if( conv.lgSearch )
63  {
64  trimhi = (realnum)ionbal.trimhi * 1e-4f;
65  trimlo = (realnum)ionbal.trimlo * 1e-4f;
66  }
67  else
68  {
69  trimhi = (realnum)ionbal.trimhi;
70  trimlo = (realnum)ionbal.trimlo;
71  }
72 
73  /* helium is special case since abundance is so high, and He+ CT with
74  * CO is the dominant CO destruction process in molecular regions */
75  if( nelem == ipHELIUM )
76  {
77  /* never want to trip up a lower stage of ionization */
78  trimlo = SMALLFLOAT;
79 
80  /* if He+ is highest stage of ionization, probably want to keep it
81  * since important for CO chemistry in molecular clouds */
82  if( dense.IonHigh[ipHELIUM] == 1 )
83  {
84  trimhi = MIN2( trimhi , 1e-20f );
85  }
86  else if( dense.IonHigh[ipHELIUM] == 2 )
87  {
88  if( conv.lgSearch )
89  {
90  /* during search phase we may be quite far from solution, in the
91  * case of a PDR sim the electron density may be vastly higher
92  * than it will be with stable solution. He++ can be very important
93  * for the chemistry and we want to consider it. Make the
94  * threshold for ignoring He++ much higher than normal to prevent
95  * a premature removal of the ion */
96  trimhi = MIN2( trimhi , 1e-17f );
97  }
98  else
99  {
100  /* similar smaller upper limit for ion*/
101  trimhi = MIN2( trimhi , 1e-12f );
102  }
103  }
104  }
105 
106  /* logic for PDRs, for elements included in chemistry, need stable solutions,
107  * keep 3 ion stages in most cases - added by NA to do HII/PDR sims */
108  if( mole.lgElem_in_chemistry[nelem] )
109  {
110  trimlo = SMALLFLOAT;
111  if(dense.IonHigh[nelem] ==2)
112  {
113  trimhi = MIN2(trimhi, 1e-20f);
114  }
115  }
116 
117  /* raise or lower highest and lowest stages of ionization to
118  * consider only significant stages
119  * IonHigh[nelem] stage of ionization */
120 
121  /* this is a special block for initialization only - it checks on absurd abundances
122  * and can trim multiple stages of ionization at one time. */
123  if( conv.lgSearch )
124  {
125  /* special - trim down higher stages if they have essentially zero abundance */
126  while(
127  (dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] < dense.density_low_limit ||
128  dense.xIonDense[nelem][dense.IonHigh[nelem]] < dense.density_low_limit ) &&
129  /* >>chng 02 may 12, rm +1 since this had effect of not allowing fully atomic */
130  dense.IonHigh[nelem] > dense.IonLow[nelem] )
131  {
132  /* dense.xIonDense[nelem][i] is density of i+1 th ionization stage (cm^-3)
133  * the -1 is correct for the heating, -1 since heating is caused by ionization of stage below it */
134  dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
135  thermal.heating[nelem][dense.IonHigh[nelem]-1] = 0.;
136  if( dense.IonHigh[nelem] == nelem+1-NISO )
137  {
138  long int ipISO = nelem - dense.IonHigh[nelem];
139  ASSERT( ipISO>=0 && ipISO<NISO );
140  StatesElem[ipISO][nelem][0].Pop = 0.;
141  }
142 
143  /* decrement high stage limit */
144  --dense.IonHigh[nelem];
145  ASSERT( dense.IonHigh[nelem] >= 0);
146  /* remember this happened */
147  lgHi_Down = true;
148  }
149 
150  /* special - trim up lower stages trim if they have essentially zero abundance */
151  while(
152  (dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] < dense.density_low_limit ||
153  dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit ) &&
154  dense.IonLow[nelem] < dense.IonHigh[nelem] - 1 )
155  {
156  /* dense.xIonDense[nelem][i] is density of ith ionization stage (cm^-3)
157  * there is no-1 since we are removing the agent that heats */
158  dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.;
159  /* >>chng 01 aug 04, remove -1 which clobbers thermal.heating when IonLow == 0 */
160  /*thermal.heating[nelem][dense.IonLow[nelem]-1] = 0.;*/
161  thermal.heating[nelem][dense.IonLow[nelem]] = 0.;
162  if( dense.IonLow[nelem] == nelem+1-NISO )
163  {
164  long int ipISO = nelem - dense.IonLow[nelem];
165  ASSERT( ipISO>=0 && ipISO<NISO );
166  StatesElem[ipISO][nelem][0].Pop = 0.;
167  }
168 
169  /* increment low stage limit */
170  ++dense.IonLow[nelem];
171  lgLo_Up = true;
172  }
173  }
174 
175  /* sanity check */
176  /* IonHigh can be equal to IonLow if both are zero - no ionization*/
177  ASSERT( dense.IonLow[nelem] < dense.IonHigh[nelem] ||
178  (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) ||
179  dense.lgSetIoniz[nelem] );
180 
181  /* this checks on error condition where the gas is stupendously highly ionized,
182  * the low limit is already one less than high, and we need to raise the low
183  * limit further */
184  if( dense.IonHigh[nelem] > 1 &&
185  (dense.IonLow[nelem]==dense.IonHigh[nelem]-1) &&
186  (dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit) )
187  {
188  fprintf(ioQQQ,
189  "PROBLEM DISASTER the density of ion %li of element %s is too small to be computed on this cpu.\n",
190  dense.IonLow[nelem]+1,
191  elementnames.chElementName[nelem]);
192  fprintf(ioQQQ,
193  "Turn off the element with the command ELEMENT XXX OFF or do not consider "
194  "gas with low density, the current hydrogen density is %.2e cm-3.\n",
196  fprintf(ioQQQ,
197  "The outer radius of the cloud is %.2e cm - should a stopping "
198  "radius be set?\n",
199  radius.Radius );
200  fprintf(ioQQQ,
201  "The abort flag is being set - the calculation is stopping.\n");
202  lgAbort = true;
203  return;
204  }
205 
206  /* trim down high stages that have too small an abundance */
207  if(
208  dense.IonHigh[nelem] > dense.IonLow[nelem] &&
209  ( (dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] <=
210  trimhi ) ||
211  (dense.xIonDense[nelem][dense.IonHigh[nelem]] <= dense.density_low_limit ) )
212  )
213  {
214  /* >>chng 03 sep 30, the atom and its first ion are a special case
215  * since we want to compute even trivial ions in molecular clouds */
216  if( dense.IonHigh[nelem]>1 ||
217  (dense.IonHigh[nelem]==1&&dense.xIonDense[nelem][1]<100.*dense.density_low_limit) )
218  {
219  dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
220  thermal.heating[nelem][dense.IonHigh[nelem]-1] = 0.;
221  if( dense.IonHigh[nelem] == nelem+1-NISO )
222  {
223  long int ipISO = nelem - dense.IonHigh[nelem];
224  ASSERT( ipISO>=0 && ipISO<NISO );
225  StatesElem[ipISO][nelem][0].Pop = 0.;
226  }
227  --dense.IonHigh[nelem];
228  lgHi_Down = true;
229  }
230  }
231 
232  /* trim up highest stages - will this be done? */
233  lgHi_Up_enable = true;
234  /* trimming can be turned off */
235  if( !ionbal.lgTrimhiOn )
236  lgHi_Up_enable = false;
237  /* high stage is already fully stripped */
238  if( dense.IonHigh[nelem]>=nelem+1 )
239  lgHi_Up_enable = false;
240  /* we have previously trimmed this stage down */
241  if( lgHi_Down )
242  lgHi_Up_enable = false;
243  /* we are in neutral gas */
245  lgHi_Up_enable = false;
246 
247  if( lgHi_Up_enable )
248  {
249  realnum abundold = 0;
250  if( nzone>2 )
251  abundold = struc.xIonDense[nelem][dense.IonHigh[nelem]][nzone-3]/
252  SDIV( struc.gas_phase[nelem][nzone-3]);
253  realnum abundnew = dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem];
254  /* only raise highest stage if ionization potential of next highest stage is within
255  * continuum array and the abundance of the highest stage is significant */
256  if( (Heavy.Valence_IP_Ryd[nelem][dense.IonHigh[nelem]] < rfield.anu[rfield.nflux-1] ||
257  Heavy.Valence_IP_Ryd[nelem][dense.IonHigh[nelem]] < phycon.te_ryd*10.) &&
258  dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] > 1e-4f &&
259  /* this checks that abundance of highest stage is increasing */
260  abundnew > abundold*1.01 )
261  {
262  /*fprintf(ioQQQ,"uuppp %li %li \n", nelem, dense.IonHigh[nelem] );*/
263  /* raise highest level of ionization */
264  ++dense.IonHigh[nelem];
265  lgHi_Up = true;
266  /* set abundance to almost zero so that sanity check elsewhere will be ok */
267  dense.xIonDense[nelem][dense.IonHigh[nelem]] = (realnum)(dense.density_low_limit*10.);
268  }
269  }
270 
271  /* sanity check
272  * IonHigh can be equal to IonLow if both are zero - no ionization*/
273  ASSERT( dense.IonLow[nelem] < dense.IonHigh[nelem] ||
274  (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) ||
275  dense.lgSetIoniz[nelem] );
276 
277  /* lower lowest stage of ionization if we have significant abundance at current lowest */
278  if( dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] > 1e-3f &&
279  dense.IonLow[nelem] > 0 )
280  {
281  /* lower lowest level of ionization */
282  --dense.IonLow[nelem];
283  lgLo_Down = true;
284  /* >>chng 01 aug 02, set this to zero so that sanity check elsewhere will be ok */
286  }
287 
288  /* raise lowest stage of ionization, but only if we are near illuminated face of cloud*/
289  else if( nzone < 10 &&
290  (dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] <= (realnum)trimlo) &&
291  (dense.IonLow[nelem] < (dense.IonHigh[nelem] - 2) ) )
292  {
293  /* raise lowest level of ionization */
294  dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.;
295  /* no minus one in below since this is low bound, already bounds at atom */
296  thermal.heating[nelem][dense.IonLow[nelem]] = 0.;
297  if( dense.IonLow[nelem] == nelem+1-NISO )
298  {
299  long int ipISO = nelem - dense.IonLow[nelem];
300  ASSERT( ipISO>=0 && ipISO<NISO );
301  StatesElem[ipISO][nelem][0].Pop = 0.;
302  }
303  ++dense.IonLow[nelem];
304  lgLo_Up = true;
305  }
306  /* test on zero */
307  else if( (dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit) &&
308  /*>>chng 06 may 24, from IonLow < IonHigh to <IonHigh-1 because
309  * old test would allow low to be raised too high, which then
310  * leads to insanity */
311  (dense.IonLow[nelem] < dense.IonHigh[nelem]-1) )
312  {
313  while(dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit &&
314  /* >>chng 07 feb 27 from < IonHigh to < IonHigh-1 to prevent raising
315  * low to high */
316  (dense.IonLow[nelem] < dense.IonHigh[nelem]-1) )
317  {
318  /* raise lowest level of ionization */
319  dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.;
320  /* no minus one in below since this is low bound, already bounds at atom */
321  thermal.heating[nelem][dense.IonLow[nelem]] = 0.;
322  if( dense.IonLow[nelem] == nelem+1-NISO )
323  {
324  long int ipISO = nelem - dense.IonLow[nelem];
325  ASSERT( ipISO>=0 && ipISO<NISO );
326  StatesElem[ipISO][nelem][0].Pop = 0.;
327  }
328  ++dense.IonLow[nelem];
329  lgLo_Up = true;
330  }
331  }
332 
333  /* sanity check */
334  /* IonHigh can be equal to IonLow if both are zero - no ionization*/
335  ASSERT( dense.IonLow[nelem] < dense.IonHigh[nelem] ||
336  (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) ||
337  dense.lgSetIoniz[nelem] );
338 
339  /* these are standard bounds checks that appear throughout this routine
340  * dense.xIonDense[IonLow] is first one with positive abundances
341  *
342  * in case where lower ionization stage was just lowered the abundance
343  * was set to dense.density_low_limit so test must be < dense.density_low_limit */
344  ASSERT( dense.IonLow[nelem] <= 1 ||
345  dense.xIonDense[nelem][dense.IonLow[nelem]-1] == 0. );
346 
347  ASSERT( (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) || lgLo_Up ||
348  dense.xIonDense[nelem][dense.IonLow[nelem]] >= dense.density_low_limit ||
349  dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] >= dense.density_low_limit ||
350  /*>>chng 06 may 24, include this to cover case where cap is set by not allowing
351  * lower limit to come up to upper limit */
352  (dense.IonLow[nelem]==dense.IonHigh[nelem]-1 ));
353 
354  {
355  /* option to print out what has happened so far .... */
356  enum {DEBUG_LOC=false};
357  if( DEBUG_LOC && nelem == ipHELIUM )
358  {
359  if( lgHi_Down ||lgHi_Up ||lgLo_Up ||lgLo_Down )
360  {
361  fprintf(ioQQQ,"DEBUG TrimZone\t%li\t",nzone );
362  if( lgHi_Down )
363  {
364  fprintf(ioQQQ,"high dn %li to %li",
365  ion_hi_save ,
366  dense.IonHigh[nelem] );
367  }
368  if( lgHi_Up )
369  {
370  fprintf(ioQQQ,"high up %li to %li",
371  ion_hi_save ,
372  dense.IonHigh[nelem] );
373  }
374  if( lgLo_Up )
375  {
376  fprintf(ioQQQ,"low up %li to %li",
377  ion_lo_save ,
378  dense.IonLow[nelem] );
379  }
380  if( lgLo_Down )
381  {
382  fprintf(ioQQQ,"low dn %li to %li",
383  ion_lo_save ,
384  dense.IonLow[nelem] );
385  }
386  fprintf(ioQQQ,"\n" );
387  }
388  }
389  }
390 
391  /* only assert that the stage above the highest has a population of
392  * zero if that stage exists */
393  if(dense.IonHigh[nelem] < nelem+1 )
394  ASSERT( dense.xIonDense[nelem][dense.IonHigh[nelem]+1] == 0. );
395 
396  /* option to print if any trimming occurred */
397  if( lgHi_Down || lgHi_Up || lgLo_Up || lgLo_Down )
398  {
399  conv.lgIonStageTrimed = true;
400  {
401  /* option to print out what has happened so far .... */
402  enum {DEBUG_LOC=false};
403  if( DEBUG_LOC && nelem==ipHELIUM )
404  {
405  fprintf(ioQQQ,"DEBUG ion_trim zone\t%.2f \t%li\t", fnzone, nelem);
406  if( lgHi_Down )
407  fprintf(ioQQQ,"\tHi_Down");
408  if( lgHi_Up )
409  fprintf(ioQQQ,"\tHi_Up");
410  if( lgLo_Up )
411  fprintf(ioQQQ,"\tLo_Up");
412  if( lgLo_Down )
413  fprintf(ioQQQ,"\tLo_Down");
414  fprintf(ioQQQ,"\n");
415  }
416  }
417  }
418 
419  /* asserts that that densities of ion stages that are now turned
420  * off are zero */
421  for( ion=0; ion<dense.IonLow[nelem]; ++ion )
422  {
423  ASSERT( dense.xIonDense[nelem][ion] == 0. );
424  }
425  for( ion=dense.IonHigh[nelem]+1; ion<nelem+1; ++ion )
426  {
427  ASSERT( dense.xIonDense[nelem][ion] == 0. );
428  }
429 
430  for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
431  {
432  /* in case where ionization stage was changed the
433  * density is zero, we set it to SMALLFLOAT since a density
434  * of zero is not possible */
435  dense.xIonDense[nelem][ion] = MAX2(dense.xIonDense[nelem][ion], SMALLFLOAT);
436  }
437  return;
438 }

Generated for cloudy by doxygen 1.8.1.1