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