00001
00002
00003
00004 #include "cddefines.h"
00005 #include "yield.h"
00006 #include "prt.h"
00007 #include "continuum.h"
00008 #include "iso.h"
00009 #include "dynamics.h"
00010 #include "grainvar.h"
00011 #include "hmi.h"
00012 #include "mole.h"
00013 #include "thermal.h"
00014 #include "thirdparty.h"
00015 #include "conv.h"
00016 #include "secondaries.h"
00017 #include "phycon.h"
00018 #include "atmdat.h"
00019 #include "heavy.h"
00020 #include "elementnames.h"
00021 #include "dense.h"
00022 #include "radius.h"
00023 #include "ionbal.h"
00024
00025 void solveions(double *ion, double *rec, double *snk, double *src,
00026 long int nlev, long int nmax);
00027
00028 void ion_solver(
00029
00030 long int nelem,
00031
00032 bool lgPrintIt)
00033 {
00034
00035 bool lgTriDiag=false;
00036 long int ion,
00037 limit,
00038 IonProduced,
00039 nej,
00040 ns,
00041 jmax=-1;
00042 double totsrc;
00043 double dennew, rateone;
00044 bool lgNegPop;
00045 static bool lgMustMalloc=true;
00046 static double *sink, *source , *sourcesave;
00047 int32 nerror;
00048 static int32 *ipiv;
00049 long ion_low, ion_range, i, j, ion_to , ion_from;
00050 static double sum_dense;
00051
00052 static double *auger;
00053
00054 double abund_total, renormnew;
00055 bool lgHomogeneous = true;
00056 static double *xmat , *xmatsave;
00057
00058 DEBUG_ENTRY( "ion_solver()" );
00059
00060
00061 ASSERT( nelem >= 0);
00062 ASSERT( dense.IonLow[nelem] >= 0 );
00063 ASSERT( dense.IonHigh[nelem] >= 0 );
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074 if( nelem == ipHYDROGEN )
00075 {
00076
00077
00078
00079
00080 abund_total = dense.xIonDense[nelem][0] + dense.xIonDense[nelem][1];
00081 }
00082 else
00083 {
00084 abund_total = SDIV( dense.gas_phase[nelem] - dense.xMolecules[nelem] );
00085 }
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098 if( nelem>ipHYDROGEN && dense.xMolecules[nelem]/SDIV(dense.gas_phase[nelem]) < 1e-10 )
00099 {
00100 for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
00101 {
00102 mole.source[nelem][ion] = 0.;
00103 mole.sink[nelem][ion] = 0.;
00104 }
00105 }
00106
00107
00108
00109
00110
00111
00112 if( fabs( dense.gas_phase[nelem] - dense.xMolecules[nelem])/SDIV(dense.gas_phase[nelem] ) <
00113 10.*FLT_EPSILON )
00114 {
00115 double renorm;
00116
00117
00118
00119
00120
00121 float sum = 0.;
00122 for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
00123 sum += dense.xIonDense[nelem][ion];
00124
00125
00126 renorm = dense.gas_phase[nelem] / SDIV(sum + dense.xMolecules[nelem] );
00127
00128 abund_total = renorm * sum;
00129 }
00130
00131
00132 if( abund_total < 0. )
00133 {
00134
00135
00136
00137
00138
00139
00140
00141 if(!conv.lgSearch )
00142 fprintf(ioQQQ,
00143 "PROBLEM ion_solver - neg net atomic abundance zero for nelem= %li, rel val= %.2e conv.nTotalIoniz=%li, fixed\n",
00144 nelem,
00145 fabs(abund_total) / SDIV(dense.xMolecules[nelem]),
00146 conv.nTotalIoniz );
00147
00148
00149 abund_total = -abund_total/2.;
00150
00151
00152
00153 conv.lgConvIoniz = false;
00154 strcpy( conv.chConvIoniz, "neg ion" );
00155 }
00156
00157
00158 if( dense.IonHigh[nelem] == 0 )
00159 {
00160
00161 dense.xIonDense[nelem][0] = (float)abund_total;
00162 DEBUG_EXIT( "ion_solver()" );
00163 return;
00164 }
00165
00166
00167 if( dense.lgSetIoniz[nelem] )
00168 {
00169 for( ion=0; ion<nelem+2; ++ion )
00170 {
00171 dense.xIonDense[nelem][ion] = dense.SetIoniz[nelem][ion]*(float)abund_total;
00172 }
00173 DEBUG_EXIT( "ion_solver()" );
00174 return;
00175 }
00176
00177
00178
00179
00180 ASSERT( (dense.IonHigh[nelem] < nelem + 1) || iso.pop_ion_ov_neut[ipH_LIKE][nelem] > 0. );
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191 limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);
00192
00193
00194 ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
00195 ASSERT( ion_range <= nelem+2 );
00196
00197 ion_low = dense.IonLow[nelem];
00198
00199
00200
00201 if( PARALLEL_MODE || lgMustMalloc )
00202 {
00203 long int ncell=LIMELM+1;
00204 lgMustMalloc = false;
00205 if( PARALLEL_MODE )
00206 ncell = ion_range;
00207
00208
00209 xmat=(double*)MALLOC( (sizeof(double)*(unsigned)(ncell*ncell) ));
00210 xmatsave=(double*)MALLOC( (sizeof(double)*(unsigned)(ncell*ncell) ));
00211 sink=(double*)MALLOC( (sizeof(double)*(unsigned)ncell ));
00212 source=(double*)MALLOC( (sizeof(double)*(unsigned)ncell ));
00213 sourcesave=(double*)MALLOC( (sizeof(double)*(unsigned)ncell ));
00214 ipiv=(int32*)MALLOC( (sizeof(int32)*(unsigned)ncell ));
00215
00216
00217
00218 auger=(double*)MALLOC( (sizeof(double)*(unsigned)(LIMELM+1) ));
00219 }
00220
00221 for( i=0; i<ion_range;i++ )
00222 {
00223 source[i] = 0.;
00224 }
00225
00226
00227 # undef MAT
00228 # define MAT(M_,I_,J_) (*((M_)+(I_)*(ion_range)+(J_)))
00229
00230
00231
00232
00233
00234 for( ion=0; ion <= limit; ion++ )
00235 {
00236 ionbal.RateIonizTot[nelem][ion] = 0.;
00237 }
00238
00239
00240
00241
00242 for( i=0; i< LIMELM+1; ++i )
00243 {
00244 auger[i] = 0.;
00245 }
00246
00247
00248 for( i=0; i< ion_range; ++i )
00249 {
00250 for( j=0; j< ion_range; ++j )
00251 {
00252 MAT( xmat, i, j ) = 0.;
00253 }
00254 }
00255
00256 {
00257
00258
00259
00260 enum {DEBUG_LOC=false};
00261
00262 if( DEBUG_LOC && nelem==ipCARBON )
00263 {
00264 broken();
00265 dense.IonLow[nelem] = 0;
00266 dense.IonHigh[nelem] = 3;
00267 abund_total = 1.;
00268 ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
00269
00270 for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
00271 {
00272 double fac=1;
00273 if(ion)
00274 fac = 1e-10;
00275 ionbal.RateRecomTot[nelem][ion] = 100.;
00276 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
00277 {
00278
00279 ionbal.PhotoRate_Shell[nelem][ion][ns][0] = fac;
00280 }
00281 }
00282 }
00283 }
00284
00285
00286
00287
00288 for( ion_to=dense.IonLow[nelem]; ion_to <= limit; ion_to++ )
00289 {
00290 for( ion_from=dense.IonLow[nelem]; ion_from <= dense.IonHigh[nelem]; ++ion_from )
00291 {
00292
00293 if( ion_to != ion_from )
00294 {
00295
00296
00297 rateone = mole.xMoleChTrRate[nelem][ion_from][ion_to];
00298 ASSERT( rateone >= 0. );
00299 MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
00300 MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
00301 }
00302 }
00303 }
00304
00305
00306
00307
00308
00309
00310 if( gv.lgDustOn && ionbal.lgGrainIonRecom && gv.lgGrainPhysicsOn )
00311 {
00312 long int low;
00313
00314
00315
00316 if( mole.lgElem_in_chemistry[nelem] )
00317
00318
00319 {
00320 low = MAX2(1, dense.IonLow[nelem] );
00321 }
00322 else
00323 low = dense.IonLow[nelem];
00324
00325 for( ion_to=low; ion_to <= limit; ion_to++ )
00326 {
00327 for( ion_from=dense.IonLow[nelem]; ion_from <= dense.IonHigh[nelem]; ++ion_from )
00328 {
00329
00330 if( ion_to != ion_from )
00331 {
00332
00333
00334 rateone = gv.GrainChTrRate[nelem][ion_from][ion_to];
00335 MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
00336 MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
00337 }
00338 }
00339 }
00340 }
00341
00342 for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
00343 {
00344
00345 rateone = ionbal.CollIonRate_Ground[nelem][ion][0] +
00346 secondaries.csupra[nelem][ion] +
00347
00348 ionbal.UTA_ionize_rate[nelem][ion];
00349 ionbal.RateIonizTot[nelem][ion] += rateone;
00350
00351
00352 if( ion+1-ion_low < ion_range )
00353 {
00354
00355 MAT( xmat, ion-ion_low, ion-ion_low ) -= rateone;
00356 MAT( xmat, ion-ion_low, ion+1-ion_low ) += rateone;
00357 }
00358
00359
00360 if( ion-1-ion_low >= 0 )
00361 {
00362
00363 MAT( xmat,ion-ion_low, ion-ion_low ) -= ionbal.RateRecomTot[nelem][ion-1];
00364 MAT( xmat,ion-ion_low, ion-1-ion_low ) += ionbal.RateRecomTot[nelem][ion-1];
00365 }
00366
00367
00368 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
00369 {
00370
00371 ionbal.RateIonizTot[nelem][ion] += ionbal.PhotoRate_Shell[nelem][ion][ns][0];
00372
00373
00374
00375
00376
00377
00378
00379
00380 if( ion+1-ion_low < ion_range )
00381 {
00382
00383 MAT( xmat, ion-ion_low, ion-ion_low ) -= ionbal.PhotoRate_Shell[nelem][ion][ns][0];
00384
00385
00386
00387
00388
00389
00390 for( nej=1; nej <= t_yield::Inst().nelec_eject(nelem,ion,ns); nej++ )
00391 {
00392
00393
00394
00395 IonProduced = MIN2(ion+nej,dense.IonHigh[nelem]);
00396 rateone = ionbal.PhotoRate_Shell[nelem][ion][ns][0]*
00397 t_yield::Inst().elec_eject_frac(nelem,ion,ns,nej-1);
00398
00399
00400
00401
00402
00403
00404
00405 MAT( xmat, ion-ion_low, IonProduced-ion_low ) += rateone;
00406
00407
00408
00409
00410 if( nej>1 )
00411 auger[IonProduced-1] += rateone;
00412
00413 }
00414 }
00415 }
00416
00417
00418 rateone =
00419 atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[ipHELIUM][1]+
00420 atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[ipHYDROGEN][1];
00421 ionbal.RateIonizTot[nelem][ion] += rateone;
00422 if( ion+1-ion_low < ion_range )
00423 {
00424 MAT( xmat, ion-ion_low, ion-ion_low ) -= rateone;
00425 MAT( xmat, ion-ion_low, ion+1-ion_low ) += rateone;
00426 }
00427 }
00428
00429
00430
00431
00432 j = MAX2(0,limit+1);
00433
00434 j = MAX2( ion_low , j );
00435 for( ion=j; ion<=dense.IonHigh[nelem]; ion++ )
00436 {
00437 ASSERT( ion>=0 && ion<nelem+2 );
00438
00439 if( ion+1-ion_low < ion_range )
00440 {
00441
00442 MAT( xmat, ion-ion_low, ion-ion_low ) -= ionbal.RateIonizTot[nelem][ion];
00443 MAT( xmat, ion-ion_low, ion+1-ion_low ) += ionbal.RateIonizTot[nelem][ion];
00444 }
00445
00446 if( ion-1-ion_low >= 0 )
00447 {
00448
00449 MAT( xmat,ion-ion_low, ion-ion_low ) -= ionbal.RateRecomTot[nelem][ion-1];
00450 MAT( xmat,ion-ion_low, ion-1-ion_low ) += ionbal.RateRecomTot[nelem][ion-1];
00451 }
00452 }
00453
00454
00455
00456 totsrc = 0.;
00457
00458
00459
00460 if( conv.nTotalIoniz > 1 || iteration > 1 )
00461 {
00462
00463 for( i=0; i<ion_range;i++ )
00464 {
00465 ion = i+ion_low;
00466
00467
00468
00469
00470 source[i] -= mole.source[nelem][ion];
00471
00472 totsrc += mole.source[nelem][ion];
00473
00474
00475 MAT( xmat, i, i ) -= mole.sink[nelem][ion];
00476 }
00477 }
00478
00479
00480 if( totsrc != 0. )
00481 lgHomogeneous = false;
00482
00483
00484
00485 if( iteration >= 2 && dynamics.lgAdvection )
00486 {
00487 for( i=0; i<ion_range;i++ )
00488 {
00489 ion = i+ion_low;
00490
00491 MAT( xmat, i, i ) -= dynamics.Rate;
00492
00493 source[i] -= dynamics.Source[nelem][ion];
00494
00495
00496 }
00497 lgHomogeneous = false;
00498 }
00499
00500 for( i=0; i< ion_range; ++i )
00501 {
00502 for( j=0; j< ion_range; ++j )
00503 {
00504 MAT( xmatsave, i, j ) = MAT( xmat, i, j );
00505 }
00506 sourcesave[i] = source[i];
00507 }
00508
00509
00510
00511
00512
00513
00514
00515
00516 if( !lgHomogeneous && ion_range==2 )
00517 {
00518
00519 double a = MAT( xmatsave, 0, 0 ),
00520 b = MAT( xmatsave, 1, 0 ) ,
00521 c = MAT( xmatsave, 0, 1 ) ,
00522 d = MAT( xmatsave, 1, 1 );
00523 b = SDIV(b);
00524 d = SDIV(d);
00525 double ratio1 = a/b , ratio2 = c/d , fratio1=fabs(a/b),fratio2=fabs(c/d);
00526 if( fabs(ratio1-ratio2)/MAX2(fratio1,fratio2) <DBL_EPSILON*1e4 )
00527 {
00528 abund_total = SDIV( dense.gas_phase[nelem] - dense.xMolecules[nelem] );
00529 lgHomogeneous = true;
00530 }
00531 }
00532 # if 0
00533 if( nelem==ipHYDROGEN &&
00534 fabs(dense.xMolecules[nelem]) / SDIV(dense.gas_phase[ipHYDROGEN]) <DBL_EPSILON*100. )
00535 {
00536 abund_total = SDIV( dense.gas_phase[nelem] - dense.xMolecules[nelem] );
00537 lgHomogeneous = true;
00538 }
00539 # endif
00540
00541
00542
00543
00544 if( lgHomogeneous )
00545 {
00546 double rate_ioniz=1., rate_recomb ;
00547
00548 jmax = 0;
00549 for( i=0; i<ion_range-1;i++)
00550 {
00551 ion = i+ion_low;
00552 rate_ioniz *= ionbal.RateIonizTot[nelem][ion];
00553 rate_recomb = ionbal.RateRecomTot[nelem][ion];
00554
00555 if( rate_ioniz == 0)
00556 break;
00557
00558 if( rate_recomb <= 0.)
00559 break;
00560
00561 rate_ioniz /= rate_recomb;
00562 if( rate_ioniz > 1.)
00563 {
00564
00565 jmax = i;
00566 rate_ioniz = 1.;
00567 }
00568 }
00569
00570
00571 for( i=0; i<ion_range;i++ )
00572 {
00573 MAT(xmat,i,jmax) = 1.;
00574 }
00575 source[jmax] = abund_total;
00576 }
00577
00578 for( i=0; i< ion_range; ++i )
00579 {
00580 for( j=0; j< ion_range; ++j )
00581 {
00582 MAT( xmatsave, i, j ) = MAT( xmat, i, j );
00583 }
00584 sourcesave[i] = source[i];
00585 }
00586
00587 if( false && nelem == ipHYDROGEN && dynamics.lgAdvection&& iteration>1 )
00588 {
00589 fprintf(ioQQQ,"DEBUGG Rate %.2f %.3e \n",fnzone,dynamics.Rate);
00590 fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateIonizTot[nelem][0], ionbal.RateIonizTot[nelem][1]);
00591 fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateRecomTot[nelem][0], ionbal.RateRecomTot[nelem][1]);
00592 fprintf(ioQQQ," %.3e %.3e %.3e\n\n", dynamics.Source[nelem][0], dynamics.Source[nelem][1], dynamics.Source[nelem][2]);
00593 }
00594
00595 {
00596
00597
00598 enum {DEBUG_LOC=false};
00599
00600 if( DEBUG_LOC && nzone==1 && lgPrintIt )
00601 {
00602 fprintf( ioQQQ,
00603 " DEBUG ion_solver: nelem=%li ion_range=%li, limit=%li, nConv %li xmat follows\n",
00604 nelem , ion_range,limit , conv.nTotalIoniz );
00605 if( lgHomogeneous )
00606 fprintf(ioQQQ , "Homogeneous \n");
00607 for( i=0; i<ion_range; ++i )
00608 {
00609 for( j=0;j<ion_range;j++ )
00610 {
00611 fprintf(ioQQQ,"%e\t",MAT(xmatsave,i,j));
00612 }
00613 fprintf(ioQQQ,"\n");
00614 }
00615 fprintf(ioQQQ,"source follows\n");
00616 for( i=0; i<ion_range;i++ )
00617 {
00618 fprintf(ioQQQ,"%e\t",source[i]);
00619 }
00620 fprintf(ioQQQ,"\n");
00621 }
00622 }
00623
00624
00625 if( lgTriDiag )
00626 {
00627 bool lgLapack=false , lgTriOptimized=true;
00628
00629 if(lgLapack)
00630 {
00631
00632
00633
00634
00635 bool lgBad = false;
00636
00637 double *dl, *d, *du;
00638
00639 d = (double *) MALLOC((unsigned)ion_range*sizeof(double));
00640 du = (double *) MALLOC((unsigned)(ion_range-1)*sizeof(double));
00641 dl = (double *) MALLOC((unsigned)(ion_range-1)*sizeof(double));
00642
00643 for(i=0;i<ion_range-1;i++)
00644 {
00645 du[i] = MAT(xmat,i+1,i);
00646 d[i] = MAT(xmat,i,i);
00647 dl[i] = MAT(xmat,i,i+1);
00648 }
00649 d[i] = MAT(xmat,i,i);
00650
00651 # if 0
00652 int i1, i2;
00653 i1 = ion_range;
00654 i2 = 1;
00655 lgBad = dgtsv_wrapper(&i1, &i2, dl, d, du, source, &i2);
00656 # endif
00657 if( lgBad )
00658 fprintf(ioQQQ," dgtsz error.\n");
00659
00660 free(dl);free(du);free(d);
00661 }
00662 else if(lgTriOptimized)
00663 {
00664
00665
00666
00667
00668
00669 for(i=0;i<ion_range;i++)
00670 {
00671 source[i] = sink[i] = 0.;
00672 }
00673 if(nelem == ipHYDROGEN && !hmi.lgNoH2Mole)
00674 {
00675 for(i=0;i<ion_range;i++)
00676 {
00677 ion = i+ion_low;
00678 source[i] += mole.source[nelem][ion];
00679 sink[i] += mole.sink[nelem][ion];
00680 }
00681 }
00682 if( iteration >= 2 && dynamics.lgAdvection && radius.depth < dynamics.oldFullDepth )
00683 {
00684 for(i=0;i<ion_range;i++)
00685 {
00686 ion = i+ion_low;
00687 source[i] += dynamics.Source[nelem][ion];
00688 sink[i] += dynamics.Rate;
00689 }
00690 }
00691
00692 solveions(ionbal.RateIonizTot[nelem]+ion_low,ionbal.RateRecomTot[nelem]+ion_low,
00693 sink,source,ion_range,jmax);
00694 }
00695 }
00696 else
00697 {
00698
00699 nerror = 0;
00700
00701 getrf_wrapper(ion_range, ion_range, xmat, ion_range, ipiv, &nerror);
00702 if( nerror != 0 )
00703 {
00704 fprintf( ioQQQ,
00705 " DISASTER ion_solver: dgetrf finds singular or ill-conditioned matrix nelem=%li %s ion_range=%li, limit=%li, nConv %li xmat follows\n",
00706 nelem ,
00707 elementnames.chElementSym[nelem],
00708 ion_range,
00709 limit ,
00710 conv.nTotalIoniz );
00711 for( i=0; i<ion_range; ++i )
00712 {
00713 for( j=0;j<ion_range;j++ )
00714 {
00715 fprintf(ioQQQ,"%e\t",MAT(xmatsave,j,i));
00716 }
00717 fprintf(ioQQQ,"\n");
00718 }
00719 fprintf(ioQQQ,"source follows\n");
00720 for( i=0; i<ion_range;i++ )
00721 {
00722 fprintf(ioQQQ,"%e\t",source[i]);
00723 }
00724 fprintf(ioQQQ,"\n");
00725 puts( "[Stop in ion_solver]" );
00726 cdEXIT(EXIT_FAILURE);
00727 }
00728 getrs_wrapper('N', ion_range, 1, xmat, ion_range, ipiv, source, ion_range, &nerror);
00729 if( nerror != 0 )
00730 {
00731 fprintf( ioQQQ, " DISASTER ion_solver: dgetrs finds singular or ill-conditioned matrix nelem=%li ionrange=%li\n",
00732 nelem , ion_range );
00733 puts( "[Stop in ion_solver]" );
00734 cdEXIT(EXIT_FAILURE);
00735 }
00736 }
00737
00738 {
00739
00740
00741 enum {DEBUG_LOC=false};
00742
00743 if( DEBUG_LOC && nelem == ipHYDROGEN )
00744 {
00745 fprintf(ioQQQ,"debuggg ion_solver1 \t%.2f\t%.4e\t%.4e\tIon\t%.3e\tRec\t%.3e\n",
00746 fnzone,
00747 phycon.te,
00748 dense.eden,
00749 ionbal.RateIonizTot[nelem][0] ,
00750 ionbal.RateRecomTot[nelem][0]);
00751 fprintf(ioQQQ," Msrc %.3e %.3e\n", mole.source[ipHYDROGEN][0], mole.source[ipHYDROGEN][1]);
00752 fprintf(ioQQQ," Msnk %.3e %.3e\n", mole.sink[ipHYDROGEN][0], mole.sink[ipHYDROGEN][1]);
00753 fprintf(ioQQQ," Poprat %.3e nomol %.3e\n",source[1]/source[0],
00754 ionbal.RateIonizTot[nelem][0]/ionbal.RateRecomTot[nelem][0]);
00755 }
00756 }
00757
00758 for( i=0; i<ion_range; i++ )
00759 {
00760
00761 ASSERT( !isnan( source[i] ) );
00762 }
00763
00764 # define RJRW 0
00765 if( RJRW && 0 )
00766 {
00767
00768 double test;
00769 for(i=0; i<ion_range; i++) {
00770 test = 0.;
00771 for(j=0; j<ion_range; j++) {
00772 test = test+source[j]*MAT(xmatsave,j,i);
00773 }
00774 fprintf(ioQQQ,"%e\t",test);
00775 }
00776 fprintf(ioQQQ,"\n");
00777
00778 test = 0.;
00779 fprintf(ioQQQ," ion %li abundance %.3e\n",nelem,abund_total);
00780 for( ion=dense.IonLow[nelem]; ion < dense.IonHigh[nelem]; ion++ )
00781 {
00782 if( ionbal.RateRecomTot[nelem][ion] != 0 && source[ion-ion_low] != 0 )
00783 fprintf(ioQQQ," %li %.3e %.3e : %.3e\n",ion,source[ion-ion_low],
00784 source[ion-ion_low+1]/source[ion-ion_low],
00785 ionbal.RateIonizTot[nelem][ion]/ionbal.RateRecomTot[nelem][ion]);
00786 else
00787 fprintf(ioQQQ," %li %.3e [One ratio infinity]\n",ion,source[ion-ion_low]);
00788 test += source[ion-ion_low];
00789 }
00790 test += source[ion-ion_low];
00791 }
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813 if( lgHomogeneous )
00814 {
00815 dense.xIonDense[nelem][ion_low] = (float)abund_total;
00816 for( i=1;i < ion_range; i++ )
00817 {
00818 dense.xIonDense[nelem][i+ion_low] = 0.;
00819 }
00820 }
00821
00822 renormnew = 1.;
00823
00824 if(iteration >= 2 && dynamics.lgAdvection &&
00825 nelem == ipHYDROGEN && hmi.lgNoH2Mole)
00826 {
00827
00828
00829
00830
00831 renormnew = 1.;
00832 }
00833 else
00834 {
00835 dennew = 0.;
00836 sum_dense = 0.;
00837
00838
00839 for( i=0;i < ion_range; i++ )
00840 {
00841 ion = i+ion_low;
00842 sum_dense += dense.xIonDense[nelem][ion];
00843 dennew += source[i];
00844 }
00845
00846 if( dennew > 0.)
00847 {
00848 renormnew = sum_dense / dennew;
00853 }
00854 else
00855 {
00856 renormnew = 1.;
00857 }
00858 }
00859
00860
00861 if( renormnew < 0)
00862 {
00863 fprintf(ioQQQ,"PROBLEM impossible value of renorm \n");
00864 }
00865 ASSERT( renormnew>=0 );
00866
00867
00868
00869
00870 lgNegPop = false;
00871 for( i=0; i < ion_range; i++ )
00872 {
00873 ion = i+ion_low;
00874
00875
00876
00877 if( source[i] < 0. )
00878 {
00879
00880
00881
00882
00883 if( source[i]<-2e-9 )
00884 fprintf(ioQQQ,
00885 " PROBLEM negative ion population in ion_solver, nelem=%li, %s ion=%li val=%.3e Search?%c zone=%li iteration=%li\n",
00886 nelem ,
00887 elementnames.chElementSym[nelem],
00888 ion ,
00889 source[i] ,
00890 TorF(conv.lgSearch) ,
00891 nzone ,
00892 iteration );
00893 source[i] = 0.;
00894
00895 if( ion == nelem+1-NISO )
00896 {
00897 long int ipISO = nelem - ion;
00898 ASSERT( ipISO>=0 && ipISO<NISO );
00899 iso.Pop2Ion[ipISO][nelem][0] = 0.;
00900 }
00901 }
00902
00903
00904 dense.xIonDense[nelem][ion] = (float)(source[i]*renormnew);
00905 }
00906
00907
00908 while( dense.IonHigh[nelem] > dense.IonLow[nelem] &&
00909 dense.xIonDense[nelem][dense.IonHigh[nelem]] < 1e-25*abund_total )
00910 {
00911 ASSERT( dense.xIonDense[nelem][dense.IonHigh[nelem]] >= 0. );
00912
00913 dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
00914 thermal.heating[nelem][dense.IonHigh[nelem]-1] = 0.;
00915
00916 --dense.IonHigh[nelem];
00917 }
00918
00919
00920
00921 ASSERT( (dense.IonLow[nelem] < dense.IonHigh[nelem]) ||
00922 (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) );
00923
00924
00925 if( lgNegPop )
00926 {
00927 fprintf( ioQQQ, " PROBLEM Negative population found for abundance of ionization stage of element %4.4s, ZONE=%4ld\n",
00928 elementnames.chElementNameShort[nelem], nzone );
00929
00930 fprintf( ioQQQ, " Populations were" );
00931 for( ion=1; ion <= dense.IonHigh[nelem]+1; ion++ )
00932 {
00933 fprintf( ioQQQ, "%9.1e", dense.xIonDense[nelem][ion-1] );
00934 }
00935 fprintf( ioQQQ, "\n" );
00936
00937 fprintf( ioQQQ, " destroy vector =" );
00938 for( ion=1; ion <= dense.IonHigh[nelem]; ion++ )
00939 {
00940 fprintf( ioQQQ, "%9.1e", ionbal.RateIonizTot[nelem][ion-1] );
00941 }
00942 fprintf( ioQQQ, "\n" );
00943
00944
00945 if( lgNegPop )
00946 {
00947 fprintf( ioQQQ, " CTHeavy vector =" );
00948 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
00949 {
00950 fprintf( ioQQQ, "%9.1e", atmdat.HeCharExcIonOf[nelem][ion] );
00951 }
00952 fprintf( ioQQQ, "\n" );
00953
00954 fprintf( ioQQQ, " HCharExcIonOf vtr=" );
00955 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
00956 {
00957 fprintf( ioQQQ, "%9.1e", atmdat.HCharExcIonOf[nelem][ion] );
00958 }
00959 fprintf( ioQQQ, "\n" );
00960
00961 fprintf( ioQQQ, " CollidRate vtr=" );
00962 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
00963 {
00964 fprintf( ioQQQ, "%9.1e", ionbal.CollIonRate_Ground[nelem][ion][0] );
00965 }
00966 fprintf( ioQQQ, "\n" );
00967
00968
00969 fprintf( ioQQQ, " photo rates per subshell, ion\n" );
00970 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
00971 {
00972 fprintf( ioQQQ, "%3ld", ion );
00973 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
00974 {
00975 fprintf( ioQQQ, "%9.1e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
00976 }
00977 fprintf( ioQQQ, "\n" );
00978 }
00979 }
00980
00981
00982 fprintf( ioQQQ, " create vector =" );
00983 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
00984 {
00985 fprintf( ioQQQ, "%9.1e", ionbal.RateRecomTot[nelem][ion] );
00986 }
00987 fprintf( ioQQQ, "\n" );
00988
00989 ContNegative();
00990 ShowMe();
00991 puts( "[Stop in ion_solver]" );
00992 cdEXIT(EXIT_FAILURE);
00993 }
00994
00995
00996
00997 if( prt.lgPrtArry[nelem] || lgPrintIt )
00998 {
00999
01000 fprintf( ioQQQ,
01001 "\n %s ion_solver DEBUG ion/rec rt [s-1] %s nz%.2f Te%.4e ne%.4e Tot abun:%.3e ion abun%.2e mole%.2e\n",
01002 elementnames.chElementSym[nelem],
01003 elementnames.chElementName[nelem],
01004 fnzone,
01005 phycon.te ,
01006 dense.eden,
01007 dense.gas_phase[nelem],
01008 abund_total ,
01009 dense.xMolecules[nelem] );
01010
01011 fprintf( ioQQQ, " %s Ioniz total " ,elementnames.chElementSym[nelem]);
01012 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01013 {
01014 fprintf( ioQQQ, " %9.2e", ionbal.RateIonizTot[nelem][ion] );
01015 }
01016 fprintf( ioQQQ, "\n" );
01017
01018
01019 fprintf(ioQQQ," %s mole sink ",
01020 elementnames.chElementSym[nelem]);
01021 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01022 {
01023 fprintf(ioQQQ," %9.2e", mole.sink[nelem][ion] );
01024 }
01025 fprintf( ioQQQ, "\n" );
01026
01027 if( dynamics.lgAdvection )
01028 {
01029
01030 fprintf(ioQQQ," %s dynm sink ",
01031 elementnames.chElementSym[nelem]);
01032 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01033 {
01034 fprintf(ioQQQ," %9.2e", dynamics.Rate );
01035 }
01036 fprintf( ioQQQ, "\n" );
01037 }
01038
01039
01040 fprintf( ioQQQ, " %s Recom total " ,elementnames.chElementSym[nelem]);
01041 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01042 {
01043 fprintf( ioQQQ, " %9.2e", ionbal.RateRecomTot[nelem][ion] );
01044 }
01045 fprintf( ioQQQ, "\n" );
01046
01047
01048 fprintf(ioQQQ," %s mole source ",
01049 elementnames.chElementSym[nelem]);
01050 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01051 {
01052 fprintf(ioQQQ," %9.2e", mole.source[nelem][ion] );
01053 }
01054 fprintf( ioQQQ, "\n" );
01055
01056 if( dynamics.lgAdvection )
01057 {
01058
01059 fprintf(ioQQQ," %s dynm sour ",
01060 elementnames.chElementSym[nelem]);
01061 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01062 {
01063 fprintf(ioQQQ," %9.2e", dynamics.Source[nelem][ion] );
01064 }
01065 fprintf( ioQQQ, "\n" );
01066 }
01067
01068
01069 fprintf( ioQQQ, " %s Coll ioniz " ,elementnames.chElementSym[nelem] );
01070 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01071 {
01072 fprintf( ioQQQ, " %9.2e", ionbal.CollIonRate_Ground[nelem][ion][0] );
01073 }
01074 fprintf( ioQQQ, "\n" );
01075
01076
01077 fprintf( ioQQQ, " %s UTA ioniz " ,elementnames.chElementSym[nelem] );
01078 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01079 {
01080 fprintf( ioQQQ, " %9.2e", ionbal.UTA_ionize_rate[nelem][ion] );
01081 }
01082 fprintf( ioQQQ, "\n" );
01083
01084
01085 fprintf( ioQQQ, " %s Phot ioniz " ,elementnames.chElementSym[nelem]);
01086 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01087 {
01088 fprintf( ioQQQ, " %9.2e",
01089 ionbal.PhotoRate_Shell[nelem][ion][Heavy.nsShells[nelem][ion]-1][0] );
01090 }
01091 fprintf( ioQQQ, "\n" );
01092
01093
01094 fprintf( ioQQQ, " %s Auger ioniz " ,elementnames.chElementSym[nelem]);
01095 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01096 {
01097 fprintf( ioQQQ, " %9.2e",
01098 auger[ion] );
01099 }
01100 fprintf( ioQQQ, "\n" );
01101
01102
01103 fprintf( ioQQQ, " %s Secon ioniz " ,elementnames.chElementSym[nelem]);
01104 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01105 {
01106 fprintf( ioQQQ, " %9.2e",
01107 secondaries.csupra[nelem][ion] );
01108 }
01109 fprintf( ioQQQ, "\n" );
01110
01111
01112 fprintf( ioQQQ, " %s ion on grn " ,elementnames.chElementSym[nelem]);
01113 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01114 {
01115 fprintf( ioQQQ, " %9.2e", gv.GrainChTrRate[nelem][ion][ion+1] );
01116 }
01117 fprintf( ioQQQ, "\n" );
01118
01119
01120
01121 fprintf( ioQQQ, " %s ion in chem " ,elementnames.chElementSym[nelem]);
01122 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01123 {
01124 fprintf( ioQQQ, " %9.2e", mole.xMoleChTrRate[nelem][ion][ion+1] );
01125 }
01126 fprintf( ioQQQ, "\n" );
01127
01128
01129 fprintf( ioQQQ, " %s chr trn ion " ,elementnames.chElementSym[nelem] );
01130 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01131 {
01132 double sum = atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[ipHELIUM][1]+
01133 atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[ipHYDROGEN][1];
01134
01135 if( nelem==ipHELIUM && ion==0 )
01136 {
01137 sum += atmdat.HeCharExcIonTotal;
01138 }
01139 else if( nelem==ipHYDROGEN && ion==0 )
01140 {
01141 sum += atmdat.HCharExcIonTotal;
01142 }
01143 fprintf( ioQQQ, " %9.2e", sum );
01144 }
01145 fprintf( ioQQQ, "\n" );
01146
01147
01148 fprintf( ioQQQ, " %s radiati rec " ,elementnames.chElementSym[nelem]);
01149 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01150 {
01151 fprintf( ioQQQ, " %9.2e", dense.eden*ionbal.RR_rate_coef_used[nelem][ion] );
01152 }
01153 fprintf( ioQQQ, "\n" );
01154
01155
01156 fprintf( ioQQQ, " %s dr old rec " ,elementnames.chElementSym[nelem]);
01157 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01158 {
01159 fprintf( ioQQQ, " %9.2e", dense.eden*ionbal.DR_old_rate_coef[nelem][ion] );
01160 }
01161 fprintf( ioQQQ, "\n" );
01162
01163
01164 fprintf( ioQQQ, " %s drBadnel rec" ,elementnames.chElementSym[nelem]);
01165 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01166 {
01167 fprintf( ioQQQ, " %9.2e", dense.eden*ionbal.DR_Badnell_rate_coef[nelem][ion] );
01168 }
01169 fprintf( ioQQQ, "\n" );
01170
01171
01172
01173 fprintf( ioQQQ, " %s rec on grn " ,elementnames.chElementSym[nelem]);
01174 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01175 {
01176 fprintf( ioQQQ, " %9.2e", gv.GrainChTrRate[nelem][ion+1][ion] );
01177 }
01178 fprintf( ioQQQ, "\n" );
01179
01180
01181
01182 fprintf( ioQQQ, " %s rec in chem " ,elementnames.chElementSym[nelem]);
01183 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01184 {
01185 fprintf( ioQQQ, " %9.2e", mole.xMoleChTrRate[nelem][ion+1][ion] );
01186 }
01187 fprintf( ioQQQ, "\n" );
01188
01189
01190 fprintf( ioQQQ, " %s chr trn rec " ,elementnames.chElementSym[nelem]);
01191 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01192 {
01193 double sum =
01194 atmdat.HCharExcRecTo[nelem][ion]*
01195 iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]*
01196 dense.xIonDense[ipHYDROGEN][1] +
01197
01198 atmdat.HeCharExcRecTo[nelem][ion]*
01199 iso.Pop2Ion[ipHE_LIKE][ipHELIUM][ipHe1s1S]*
01200 dense.xIonDense[ipHELIUM][1];
01201
01202 if( nelem==ipHELIUM && ion==0 )
01203 {
01204 sum += atmdat.HeCharExcRecTotal;
01205 }
01206 else if( nelem==ipHYDROGEN && ion==0 )
01207 {
01208 sum += atmdat.HCharExcRecTotal;
01209 }
01210 fprintf( ioQQQ, " %9.2e", sum );
01211 }
01212 fprintf( ioQQQ, "\n" );
01213
01214
01215 fprintf( ioQQQ, " %s Abun [cm-3] " ,elementnames.chElementSym[nelem] );
01216 for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
01217 {
01218 fprintf( ioQQQ, " %9.2e", dense.xIonDense[nelem][ion] );
01219 }
01220 fprintf( ioQQQ, "\n" );
01221 }
01222
01223 if( PARALLEL_MODE )
01224 {
01225 free(ipiv);
01226 free(source);
01227 free(sink);
01228 free(xmat);
01229 free(xmatsave);
01230 free(auger);
01231 }
01232
01233 DEBUG_EXIT( "ion_solver()" );
01234 return;
01235 }
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264 void solveions(double *ion, double *rec, double *snk, double *src,
01265 long int nlev, long int nmax)
01266 {
01267 double kap, bet;
01268 long int i;
01269
01270 if(nmax != -1)
01271 {
01272
01273 src[nmax] = 1.;
01274 for(i=nmax;i<nlev-1;i++)
01275 src[i+1] = src[i]*ion[i]/rec[i];
01276 for(i=nmax-1;i>=0;i--)
01277 src[i] = src[i+1]*rec[i]/ion[i];
01278 }
01279 else
01280 {
01281 i = 0;
01282 kap = snk[0];
01283 for(i=0;i<nlev-1;i++)
01284 {
01285 bet = ion[i]+kap;
01286 if(bet == 0.)
01287 {
01288 fprintf(ioQQQ,"Ionization solver error\n");
01289 puts("[Stop in solveions]");
01290 cdEXIT(EXIT_FAILURE);
01291 }
01292 bet = 1./bet;
01293 src[i] *= bet;
01294 src[i+1] += ion[i]*src[i];
01295 snk[i] = bet*rec[i];
01296 kap = kap*snk[i]+snk[i+1];
01297 }
01298 bet = kap;
01299 if(bet == 0.)
01300 {
01301 fprintf(ioQQQ,"Ionization solver error\n");
01302 puts("[Stop in solveions]");
01303 cdEXIT(EXIT_FAILURE);
01304 }
01305 src[i] /= bet;
01306
01307 for(i=nlev-2;i>=0;i--)
01308 {
01309 src[i] += snk[i]*src[i+1];
01310 }
01311 }
01312 }