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