Interaction between macro_pops and saha/lucy_mazzali

[jhmatthews]

Originally I flagged this bug up as a concern about the lucy_mazzali routine interacting with macro_pops, but in the end we discovered that this was due to a problem to do with saha abundances and convergence criteria. This is a macro specific problem that is understood (to an extent) but unsolved, and requires some thought about the best fix -it will require a combination of stopping saha abundances being calculated for macro ions, and some thought about the best way to converge on macro level populations. There may also be questions that need to be asked about why we are even calling the saha function during the ionization cycles, as we should surely use the output from the last cycle as our initial guess and then calculate the corrected saha equation without overwriting everything.

The problem is broadly described as follows (see comments below for the progress on the bug over time and more discussion).

• in Lucy Mazzali ionisation mode (mode 3 in the pf file, and mode 2 in nebular concentrations), nebular_concentrations calls saha(), then macro_pops, and these are both inside a loop which has a convergence criteria of a fractional error of 0.03
• saha populates the xplasma->density arrays with saha ion abundances
• In the old version of the code, macro_pops was never called again, which means that the escape probabilities used to calculate the rate matrices (which are then inverted to get level pops) were incorrect as the densities used would be some weird saha/macro_pops hybrid
• Furthermore, in the old version of the code partition functions were different and lead
• in the new version, macro_pops was called again in lucy, which meant that it did another loop over ne, but the difference being
• An investigation of the convergence criteria implied that NEITHER version was actually converging on the correct value.

with FRACTIONAL_ERROR changed to 1e-12 in both versions, and saha not called after the first ionization cycle, we were able to get the following answers

58: !!wind_update: Wind luminosity  3.31e+36 (recomb 1.41e+35 ff 8.18e+33 lines 3.16e+36) after update
76: !!wind_update: Wind luminosity  3.23e+36 (recomb 1.40e+35 ff 8.31e+33 lines 3.08e+36) after update

which I am happy to declare this close to a success and proceed with a better way of getting the level populations in macro atom mode, provided that the spectra and level emissivities are fairly similar.

[jhmatthews]

We discussed this in the 130814 telecon with SS, NSH, JM and CK. By going through the code we came to a rather different conclusion than the one above.

in Lucy Mazzali ionisation mode (mode 3 in the pf file, and mode 2 in nebular concentrations), nebular_concentrations does the following

partition_functions (xplasma, mode);  // t_r with weights

m = concentrations (xplasma, 0);  // Saha equation using t_r

 m = lucy (xplasma);   // Main routine for running LucyMazzali

i.e., it works out the partition functions, then works out saha abundances, then applies the lucy mazzali correction factors ( in non macro mode). In macro mode, both of these functions also call

There are a number of problems with this

• in non-macro mode, there are no actual problems here, but it is worth noting that we don't need to iterate to find in ne in both the saha part of the code (concentrations()) and the lucy part (lucy()). We suspect the reason this is done is simply because there is an option to run just in LTE and so is somewhat historical.
• in macro atom mode there are a few problems. The first, is that before macro pops is called at all then we call saha(), which means that the initial guess at ion densities is always that of the Saha abundances- not what we should really be doing, as ideally we want to use the abundances from the last cycle.
• however, this shouldn't matter if the convergence loop is set correctly, because the loop should essentially forget what it's initial guess is. But, we don't think this is happening because of two quantities defined in python.h:

#define MAXITERATIONS   200 //The number of loops to do to try to converge in ne
#define FRACTIONAL_ERROR 0.03   //The change in n_e which causes a break out of the loop for ne

So these need to be set differently for macro pops to obtain the correct answer. Note also that MAXITERATIONS is only set to 20 in python 58.

Are we confident that they are good enough for the non macro atom mode?

So, what does this mean?

• macro pops is not really converging correctly in both versions of the code
• the differences in line luminosity and bf luminosity between 58 and 76 (see matom line + bf luminosity bug #37 ) appear to be at least in part due to the fact that the Saha equation gives different answers because it calculates different partition functions, and because of this and the fact the convergence is not robust the error propagates through
• an initial test of changing the fractional error to a smaller fraction gives a substantially different line luminosity (factor of 2 ish)...although not necessarily one that agrees between the versions...testing...
• need to also test the sensitivity to the saha abundances by using some kind of macro atom modde if statement

[kslong]

So here are a couple of comments about all of this.

First, it is definitely true that much of what is is contained in the ionization calculation part of the routine is "historical" and developed as we were trying to understand how to do the ionization calculation. We also were only really testing situations where the plasma was essentially fully ionized, so that ne did not change much. Indeed, my understanding at the time was, that the only reason we had to iterate was to get a better estimate of ne.

I presume, when James talks about needing more iterations, he is talking about cases where the temperature is fairly low, and this is the reason ne keeps changing.

Even so, generally speaking, it seems odd to me that we should under normal situations have to iterated for 200 times and still just get a fractional error of 0.03, and if we are getting close to this, we ought to try to figure out if there is a better approach. If memory serves me correctly, it is possible to write down the set of equations that includes all of the ion abundances and ne in a way that avoids the iteration. I thought that was what Stuart did in his program.

[ssim]

Hi Knox - hope you're enjoying the beach!

The problem isn't that ne doesn't converge - probably ne converges immediately (without iteration).

The problem is that the macro atom level populations are currently being overwritten by LTE populations, which are being inserted as part of the process of computing LTE ion populations prior to going into the Lucy/Mazzali correction stuff. That shouldn't have been allowed - the "density" quantities associated with the macro atom levels should not have been allowed to be modified by anything other than "macro pops" - the effect of this is that the escape probabilities used to solve the rate equations will not be consistent with the level populations themselves, which will have a knock on effect.

If there had been iterations for ne occurring then the situation would have been somewhat better (because that part of the code does not interfere with the macro atom populations - so during the iteration on ne it should have effectively been iterating on the macro atom populations too, which would be fine). However, that's not really the elegant way to fix it - the correct thing would be to somehow forbid "saha" from overwriting the macro atom densities...or avoid the "saha" call completely.

Cheers from a wet and grey Belfast!
Stuart

[jhmatthews]

Yes, sorry I probably didn't explain that clearly enough.

A few interesting tests- as yet I've just been testing the sensitivity to FRACTIONAL_ERROR and N-ITERATIONS in both 58 and 76...the results are a little odd. Remember I'm still doing these diagnostics with one cycle runs for ease and simplicity, and it is not about whether we converge on ne but rather whether we run macro_pops enough times after calling saha that we forget about our initial guess, which is first of all not the guess we should be using, and secondly is a different guess between the 2 versions.

• First of all, here's what happens to the wind luminosity in 58 (default values in 58 are 20 and 0.03):

N_ITERATIONS	FRACTIONAL_ERROR	wind luminosity
20	0.03	3.76e+36 (recomb 1.93e+35 ff 8.98e+33 lines 3.56e+36)
200	0.03	3.76e+36 (recomb 1.93e+35 ff 8.98e+33 lines 3.56e+36)
20	1e-8	1.11e+36 (recomb 7.46e+34 ff 8.98e+33 lines 1.02e+36)
200	1e-8	1.11e+36 (recomb 7.46e+34 ff 8.98e+33 lines 1.02e+36)
20000	1e-10	1.10e+36 (recomb 7.44e+34 ff 8.98e+33 lines 1.02e+36)

• So it does have a pretty sizeable effect, which is encouraging (well, encouraging in that it helps with the bug at least, not so much in terms of accuracy of old models). The fractional error here is the limiting factor, rather than N_ITERATIONS, and changing it 1e-8 seems to be a good enough convergence criteria (note less may also work).
• So, naively, we might hope that changing the same things in 76 might lead to the same values coming out...here's the same table for 76 (default values in 76 are 200 and 0.03)::

N_ITERATIONS	FRACTIONAL_ERROR	wind luminosity
20	0.03	1.99e+36 (recomb 1.15e+35 ff 8.99e+33 lines 1.86e+36)
200	0.03	1.99e+36 (recomb 1.15e+35 ff 8.99e+33 lines 1.86e+36)
20	1e-8	1.98e+36 (recomb 1.15e+35 ff 8.99e+33 lines 1.86e+36)
200	1e-8	1.98e+36 (recomb 1.15e+35 ff 8.99e+33 lines 1.86e+36)
20000	1e-10	1.98e+36 (recomb 1.15e+35 ff 8.99e+33 lines 1.86e+36)

..but of course, it wouldn't be that simple! You'll see that in 76, there is almost no sensitivity to these two variables...

Mysteries:

• if this is the only bug, we should get the same value!
• regardless of whether it's the only bug, if it changes this much in 58, why does it barely change at all in 76?

Ideas

• Are LTE level populations actually a fairly decent guess in 76, but not in 58 due to the partition function difference, which means that 0.03 is actually good enough to converge on the correct populations in our most recent version?
• is there another bug which means that in 58 there is a great sensitivity to this but there isn't in 76, and that bug could then also explain why we don't get the same value out when we use more robust convergence criteria?

Conclusions:

• changes in level populations can have a significant effect on line and fb luminosities
• changes in ion densities used as inputs to escape probability calculations can have a significant impact on the output level populations, and in cases where the convergence criteria for ne are not sufficient for the macro_pops routine to forget it's initial guess this will result in an incorrect calculation
• there must be another bug, or something linked to this which goes wrong somewhere else

Things to do:

• first step, code a temporary solution such that saha only gets called before the ionization cycles- to follow.

[ssim]

Hi James,

Sounds like you are making some progress - I agree that it sounds odd that v76 doesn't seem to care about this though.

I guess I'm not up to speed - why are the partition functions different, and what values do they have in the two cases (I would have expected for H you always get "2" for HI and "1" for HII, pretty much regardless of what one does!? Are the partition functions even computed for the macro atoms?

Stuart

[jhmatthews]

Ok, this is getting really odd. When I change the code so that the only time saha gets called for macro ions is in the wind initialisation, and also keep FRACTIONAL_ERROR and MAX_ITERATIONS at 1e-8 and 20000, then I get 58 to agree ALMOST PERFECTLY with 76. Here's the actual figures when you don't call saha in either and both have the same improved convergence conditions:

58: Wind lum 1.99e+36 (recomb 1.16e+35 ff 8.99e+33 lines 1.87e+36)
76: Wind lum 1.98e+36 (recomb 1.15e+35 ff 8.99e+33 lines 1.86e+36)

What is going on?? This is bizarre, but without understanding it yet I think it may be almost solved!

I'll take a look at the partition functions in a moment- the only relevance they have here is that they are used in the saha equation, but they shouldn't be used at all in a macro atom scheme if it's done properly.

[jhmatthews]

Right, just talking with Christian and I think we've cracked it. The reason the above happens is because, actually, the lucy_mazzali1 routine in 58 doesn't call macro_pops at all (something I'd missed because we'd been debugging 76 in the telecon, woops), but it does in 76.

So here's what happens

• in nebular_concentrations it starts out with a guess of the saha abundances, and every time it goes round the loop in saha it also calls the saha equation.
• in 58, the saha equation has incorrect, or at least worse, partition functions, which means that your initial guess is wrong, but also that the density arrays are populated with incorrect saha abundances directly before every call to macro_pops
• this means that as long as you still call saha in that loop, not matter what you change FRACTIONAL_ERROR to you will never be calculating the correct level populations
• the fact the the initial saha guess is wrong also means that, even if you remove saha from the loop, the FRACTIONAL_ERROR is not low enough that you forget about your initial guess
• but if you remove saha from the loop AND change FRACTIONAL_ERROR then you get the right answer
• in 76, the reason this doesn't matter i two fold. 1, it also calls macro pops again in lucy mazzali so has a chance to 'forget' it's bad guess from the saha loop, and 2, the saha guess is already much closer to the correct answer than it is in 58

Hope that makes some sort of sense, it's still a bit fuzzy in my head, but I think we might've finally cracked it!

The bad news of course, is that this could mean the line strengths in Sim et al. 2005 are out by a factor of ~2- I'm currently running a test with 58 to see how this changes the final spectrum, and if that spectrum will even agree with 76. I'll also probably need to get hold of 55 to do a test of Sim 2005 properly, as that is the version stuart used.

[kslong]

So do I infer from all of this, that the problem is isolated to the macro-atom case? I hope so. Or has the fact that we set the convergence criterion to 0.03 also produced some inaccurate results.

Also, since you have investigated a lot of time in this, please add a lot of comments to the various code bits to explain better why we/I did what we did. We ought to make sure the code is clean as possible as well, eliminating any options that are unused, and explaining what situations the various ones are used for. Ideally some of the modes, which appear to be transferred as numbers, should be replaced by #defined variables, so that they are more transparent. We probably should use this as an opportunity to explain in the "user's manual" section of the google site, how each of the various modes is calculated at least in general terms.

[jhmatthews]

Yes, I think it is unique to macro atoms, although we may want to think about how that value is set.

Agreed that we'll need some cleaning up, and substantial commenting of this section.

[jhmatthews]

Ok, so I think the best solution to this is going to be as follows. saha is called in concentrations, and, if we are going to preserve the current way of doing things then this still needs to be done so that python can treat some ions and simple ions and some as macro.

The function saha loops over all ions and populates the xplasma structure with saha ion densities. Here's the relevant loop:

for (nion = first + 1; nion < last; nion++)
    {
      b = xsaha * partition[nion]
        * exp (-ion[nion - 1].ip / (BOLTZMANN * t)) / (ne *
                               partition[nion -
                                     1]);
//        if (b > big && nh < 1e5) this is a line to only modify things if the density is high enough to matter
      if (b > big)
        b = big;        //limit step so there is no chance of overflow
      a = density[nion - 1] * b;
      sum += density[nion] = a;
      if (density[nion] < 0.0)
        mytrap ();
      if (sane_check (sum))
        Error ("saha:sane_check failed for density summation\n");
    }

      a = nh * ele[nelem].abun / sum;
      for (nion = first; nion < last; nion++)
    {
      density[nion] *= a;
      if (sane_check (density[nion]))
        Error ("saha:sane_check failed for density of ion %i\n", nion);
    }
    }

So I think two if statements and substantial commenting here should do the trick

          if ((ion[nion].macro_info == 0) || (geo.macro_ioniz_mode == 0))
          {

inside each loop of nion. This means the loops will only be entered if we are before the ionization cycles and want our initial saha guess at densities ( geo.macro_ioniz_mode == 0 ) or if the ion in question is being treated as a simple ion (ion[nion].macro_info == 0).

We will also need to do something similar once the macro atom mode is implemented in ionization mode 7 ( #43 ), and we probably don't need to call macro_pops in both concentrations and lucy. And then we may need to think about FRACTIONAL_ERROR, if it's correct, and if we even want macro_pops to be called in the same place with the ne criteria.

[jhmatthews]

We now understand this problem, but a few questions still remain

• should the macro atom convergence be done over something else, i.e. not a loop over ne
• if so
  • is the fractional error appropriate just for converging on ne anyway? 0.03 is not exactly stringent
• If not
  • what do we set FRACTIONAL_ERROR to to get a matom level populations and ion densities convergence?
  • does setting it to a significantly lower value cause a slow down in either matom mode or non matom mode?

I have done some tests with a fiducial style model, but with mode 3 ionization so we call saha. It seems to suggest that the routine doesn't spend much time in this routine at all and so decreasing fractional error should be fine in that regard. Thus, an easy fix may be just to set this to a low value, although I'm not sure how reliable this is. If we set it to be 1/ne_max, where ne_max is the highest thing we care about then that should be stringent enough in most situations.

We also have an error report for 'does not converge on ne' ... with max iterations set to 200, this was never called in my test model.

Slightly unsure as to how to proceed with this one in a reliable way. Perhaps the only solution is to rewrite the procedure in matom mode, but it could be tricky given that we want simple ions to be treated in the normal way (which was presumably the motivation for doing it this way in the first place).

[jhmatthews]

I'm closing this issue, as the FRACTIONAL_ERROR problem is now identified elsewhere.