bug in iso_level.cpp?
Marios Chatzikos
Going over the matrix system in iso_level.cpp in the context of level
populations in cooling sims, I think I have spotted a conceptual error.
The code snippet in question is the following:
374 /* ionization from/recombination from lower
ionization stages */
375 for(long ion_from=dense.IonLow[nelem]; ion_from <
MIN2( dense.IonHigh[nelem], nelem-ipISO ) ; ion_from++ )
376 {
377 if(
ionbal.RateIoniz[nelem][ion_from][nelem-ipISO] >= 0. )
378 {
379 /* ionization from lower ionization
stages, cm-3 s-1 */
380 source +=
ionbal.RateIoniz[nelem][ion_from][nelem-ipISO] *
dense.xIonDense[nelem][ion_from];
381 }
382 /* recombination to next lower ionization
stage, s-1 */
383 if( ion_from == nelem-1-ipISO )
384 sink +=
ionbal.RateRecomTot[nelem][ion_from];
385 }
386
387 ASSERT( source >= 0.f );
388 if (0)
389 {
390 /*
391 * Collisional ionization and
photoionization can only be to the ground state in H iso-sequences
392 * to conserve energy.
393 */
394 creation[0] += source;
395 /*for( level=0; level < numlevels_local;
level++ )
396 {
397 z[level][level] += sink;
398 }*/
399 /*recombination is only done from the ground
state */
400 z[0][0] += sink;
401 }
402 else
403 {
404 // Try Boltzmann weighting to capture LTE
limit correctly
405 t_iso_sp* sp = &iso_sp[ipISO][nelem];
406 double partfun=0.0;
407 for ( level = 0; level < numlevels_local;
level++ )
408 {
409 partfun +=
sp->st[level].Boltzmann()*sp->st[level].g();
410 }
411 source /= partfun;
412 for( level=0; level < numlevels_local; level++ )
413 {
414 creation[level] += source*
415 sp->st[level].Boltzmann()*sp->st[level].g();
416 z[level][level] += sink;
417 }
418 }
419 creation[0] += ctsource;
420 z[0][0] += ctsink;
Lines 404-417 carry out the partition of the source terms among the
various levels. The sink term, however, in line 416 is used for all
levels without partitioning. As shown in lines 375-385, it contains the
total recombination rate into the ion.
This is what seems like a conceptual error to me. The total
recombination should either be partitioned among the levels (although
this is probably inaccurate), or added to the ground state as in line
400 (although this is certainly inaccurate). Level-resolved
recombination ought to be the answer here, but that being years away, I
think partitioning the total recombination rate should be favored.
Yet, this is probably not an oversight, so if you cared to share your
thoughts, that'd be great.
Thanks,
Marios
Gary J. Ferland
--
--
http://groups.google.com/group/cloudy-dev
---
You received this message because you are subscribed to the Google Groups "cloudy-dev" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cloudy-dev+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cloudy-dev/43f83ab6-1bb4-6fbc-aeb3-0cc7dec898c9%40gmail.com.
Physics, Univ of Kentucky
Lexington KY 40506 USA
Tel: 859 257-8795
https://pa.as.uky.edu/users/gary
Marios Chatzikos
$ svn log -r 6023
------------------------------------------------------------------------
r6023 | rporter | 2012-03-17 13:57:04 -0400 (Sat, 17 Mar
2012) | 1 line
trunk - merge newmole branch onto trunk
------------------------------------------------------------------------
To view this discussion on the web visit https://groups.google.com/d/msgid/cloudy-dev/CAOxTjFSJGq5aTxuZrSpFh76dt07Efqc1HZedEevJ1%2Bj7afH6xQ%40mail.gmail.com.
Gary J. Ferland
To view this discussion on the web visit https://groups.google.com/d/msgid/cloudy-dev/1eef49f0-656d-2a00-5eed-c96062be2ab8%40gmail.com.
Marios Chatzikos
Sure.
To view this discussion on the web visit https://groups.google.com/d/msgid/cloudy-dev/CAOxTjFSfGtPtf93cXinTXSBQK4DJGuAZo0G1Y2c8DSKgQ1gRrg%40mail.gmail.com.
Francisco Guzman Fulgencio
Sent: Monday, November 7, 2016 6:04:52 PM To: cloudy-dev Subject: Re: [cloudy-dev] RydDep botches We met to discuss the status of RydDep - the issue of the partition function is serious but not something that can be fixed easily. Currently we put all source / sink terms into iso_levels with weighting by the partition function. That is not correct since each atomic process deposites the product in a particular state. For instance, our valence photoionization cross sections populate the ground state of the parent ion. Some charge exchange reactions populate particular excited states. Most chemical / CX reactions don't have the state specified, although it probably could be figured out. Not a job for this release. So our idea is to leave iso_level as it is now, and has been for quite some time, and mark it as a priority for future development. Leaving the partition function the old way, RydDep gets the following serious botches in auto. The biggest change was a bugfix - the iso sequences had used Vriem & Smeets for all ions. Their theory was only for neutrals. We now use Percival & Richards for ions. That changes strong lines in the dense BLR sims because, probably, the UV / XUV lines of He II change - we can see that where they are directly monitored, and they produce ionizing radiation which affects the thermal balance. We have seen changes like this in the BLR sims before. So this seems fine to me. Thoughts? thanks, Gary
Sent: Thursday, February 18, 2021 9:19 AM
To: cloud...@googlegroups.com <cloud...@googlegroups.com>
Subject: Re: [cloudy-dev] bug in iso_level.cpp?
Francisco Guzman Fulgencio
Sent: Thursday, February 18, 2021 9:26 AM
Marios Chatzikos
Hi Fran,
I should clarify that I am not taking issue with the partition
function per se. I recognize that this would be a major
development effort.
My point is that the total recombination rate should not be added to each and every level. Either partition it among the levels, or add it to ground only.
Thanks,
Marios
To view this discussion on the web visit https://groups.google.com/d/msgid/cloudy-dev/BN7PR07MB5284099585C102C346FF05B3E0859%40BN7PR07MB5284.namprd07.prod.outlook.com.
Robin Williams
For the sinks, you'd expect the rate to be the product of phase space factors, Boltzmann suppression max[1,exp(-Delta E/kT)] (like Metropolis-Hastings MCMC, NB interactions with partners with excitations may be exothermic), and possibly quantum forbidden rules, averaged over every process that causes coupling. That's a lot of detail to cover. Spreading the loss across all levels is an "uninformative prior," which might even be a decent guess at least when kT >> Delta E.
To view this discussion on the web visit https://groups.google.com/d/msgid/cloudy-dev/45bd7f38-d424-5771-7744-f3988d42b1fc%40gmail.com.
Marios Chatzikos
I am not arguing that the recombination should come out of the ground state. Your experience clearly shows otherwise. But the 'uninformative prior' might be preferable for conservation.
The following is from iso_level.cpp on branch newmole at r2500:
372 /* derive source/sink terms for ion/rec
terms coupling resolved atoms
373 * to lower ionization stages */
374 if( nelem-ipISO >= 1 &&
ionbal.RateIonizTot[nelem][nelem-ipISO-1] > 0.)
375 {
376 double sum_recom =0. ,
sum_ioniz = 0.;
377 sum_popn_ov_ion = 0.;
378 /* if lower stage of
ionization is itself resolved then we must
379 * sum over rates between
levels and its continuum, the ion we
380 * are now solving */
381 if( ipISO+1<NISO )
382 {
383 for( level=0; level
< iso.numLevels_local[ipISO+1][nelem]; level++ )
384 {
385
sum_popn_ov_ion += StatesElem[ipISO+1][nelem][level].Pop;
386
387 /* sum of all
ionization processes from levels to ion we
388 * now
consider - use difference in ionization / recombination
389 * since at
high densities each of these terms are very
390 * large but
nearly cancel */
391 double one =
iso.RateCont2Level[ipISO+1][nelem][level] -
392
StatesElem[ipISO+1][nelem][level].Pop *
393
iso.RateLevel2Cont[ipISO+1][nelem][level];
394
395 /* sign
determines whether this is net ionization or
396 *
recombination process. Ground & metastable states have
397 * large
departure coefficients so are net ionization,
398 * excited
states have departure coefficients that
399 *
asymptotically approach unity from below, so are net
400 *
recombination */
401 if( one >
0 )
402
sum_recom += one;
403 else
404
sum_ioniz -= one;
405 }
406 ASSERT(
sum_recom>=0. && sum_ioniz>=0. );
407 sum_ioniz /=
MAX2(SMALLDOUBLE , sum_popn_ov_ion);
408 }
409 else
410 {
411 /* lower stage is not
resolved */
412 sum_ioniz =
ionbal.RateIonizTot[nelem][nelem-ipISO-1];
413 sum_recom =
ionbal.RateRecomTot[nelem][nelem-ipISO-1];
414 }
415 /* now add ionization from
next lower stage */
416 source += sum_ioniz*
417
dense.xIonDense[nelem][nelem-ipISO-1];
418
419 /* now add recombination to
next lower stage, units are s-1 */
420 sink += sum_recom;
421 }
422
423 creation[0] +=
source/SDIV(dense.xIonDense[nelem][nelem+1-ipISO]);
424 for( level=0; level <
numlevels_local; level++ )
425 {
426 z[level][level] += sink;
427 }
It looks like that the total recombination has always been used as a sink for every single level, and I simply don't see the reason for that.
Line 413 above is all that survived to the current version of the file (line 384). It was used for unresolved species, but even for the resolved ones the total recombination was just lumped together.
Marios
To view this discussion on the web visit https://groups.google.com/d/msgid/cloudy-dev/CAJBtPxeFO8%2BPKFeWcrad4MS4%3D6yQH8zfQ%2BYPZON7NPyckGejxw%40mail.gmail.com.
Robin Williams
To view this discussion on the web visit https://groups.google.com/d/msgid/cloudy-dev/988e0870-7971-6e21-e669-6b910dde37d5%40gmail.com.
Marios Chatzikos
No doubt, but the sum over all levels should be equal to the
total recombination rate. At the moment we are effectively
inflating the recombination rate by a factor equal to the number
of levels.
To view this discussion on the web visit https://groups.google.com/d/msgid/cloudy-dev/CAJBtPxfkSOnfj8gKknKQqPOX%3DhyY%3DDNLmt15ZfVYi8%2BmSnT9Kw%40mail.gmail.com.
Robin Williams
To view this discussion on the web visit https://groups.google.com/d/msgid/cloudy-dev/67bf536a-ed88-d538-c618-9cae036a6ae6%40gmail.com.
Marios Chatzikos
To view this discussion on the web visit https://groups.google.com/d/msgid/cloudy-dev/CAJBtPxfW_2vJMD9gjOJB5NDS3MdT6vtWq3yX7puZrX4USSZKRA%40mail.gmail.com.
Francisco Guzman Fulgencio
RateRecomTot is the sum over all recombination processes (RR, DR, etc) times the electron density. The RR data, say, come from Badnell 2006:
https://ui.adsabs.harvard.edu/abs/2006ApJ...651L..73B
https://ui.adsabs.harvard.edu/abs/2006ApJS..167..334B
which are summed over final states (see end of section 2.2 in the second reference).
That is, the rates in RateRecomTot refer to the entire ion. It seems odd to me that we use the total rate for each level separately.
To: 'Robin Williams' via cloudy-dev <cloud...@googlegroups.com>
Robin Williams
To view this discussion on the web visit https://groups.google.com/d/msgid/cloudy-dev/BN7PR07MB52842D09AEFF0FD4120C0C5DE0849%40BN7PR07MB5284.namprd07.prod.outlook.com.