Shooting method solver failing for Fe

[kslong]

Nick showed fairly convincing evidence that our shooting method for ionization balance of the Lucy Mazalli method is failing for Fe. He thinks, and I agree, that the most likely problem is due to the fact that we have a large number of lower level ions, with very low abundances, and this causes a problem when you get to high ionization states.

If this is the case, then the simple way out of this problem would be to determine in advance where the most abundant ions are for a given set of conditions, and work out from there in determining the ratio.

It would be desirable to sort out this problem, especially if we are going to leave the shooting method solver in the code when the code becomes public, but also to verify the matrix inversion techique.

[jhmatthews]

Personally, I'm not worried about actually 'fixing' this / I believe this was the motivation for putting in the variable temperature scheme.

I would prefer to have only the matrix method in the code come public release, or at least gradually phase out the shooting method. Even if we do have the shooting method in the code then I don't think this is a particularly big issue anyway. The user would have at least 2 or 3 methods available to them which don't have this problem and make as good, if not better, implicit assumptions.

The past week or so's discussion have shown that all this can be a significant time sink, as understanding the subtleties of the schemes is not necessarily easy. I think we just need to

• be able to verify the matrix scheme, which we can do well enough without this I think
• test that this is not affecting my CV models in any significant way (which it won't be) by e.g. removing iron, or running it with different ionization schemes.

[kslong]

I don't agree with this statement

Personally, I'm not worried about actually 'fixing' this / I believe this was the motivation for putting in the variable temperature scheme.

... or at least the last portion of it. The variable temperature scheme was put into place because the spectra in portions of the wind deviated significantly from a a blackbody, and this allowed us to make a better estimate of the actual ionization. For a pure BB spectrum, the two should give identical results..

I fully support the idea that we want to move to the matrix approach, as a better numerical way to solve the ionization for any set of assumptions about how the ionization should be calculated, since I think it is likely to be more stable than the shooting method, which tends to build up errors when you have lots of ions in an atom.

[kslong]

I've made a certain amount of progress in making the shooting method more robust for the lucy_mazzali part of saha.c. Instead of working from the lowest ion; the new method works outward from what is nominally the most abundant ion. As far as I can tell the results are identical for low z ions like oxygen, but produce the following for a Fe

Shooting method

This is to be compared with the rate matrix approach for the same set of conditions

The new shooting method approach avoids the jagged results we were seeing previously, but there are still some differences between the two approaches. It is not clear why any differences reamin. Note that the shooting method is not converging at the lowest temperatures. The new version of saha is in the matrix_ksl branch

[jhmatthews]

Ok, so by variable temperature I mean solving the saha equation at a temperature which gives ratios of around 1, then applying a correction factor. This was put in for practical/numerical reasons. The other issue regarding ionization relates to the modelling of the SED (i.e. dilute BB/spectral model) and is more motivated by a physical problem.

''For a pure BB spectrum, the two should give identical results..'' - this is what I've been trying to get at with my emails and what Stuart pointed out in the telecon. This is only true if you have the right temperature dependence in your recombination coefficients and in general, I think especially for more complex ions, this is not true. So one can say they should probably agree relatively well but I don't think it's really possible to know exactly how well they should do without fudging the cross sections.

I guess it's good that this has been implemented, though. I'm a little confused as to the shape of the rate matrix graph off to the left- why does the shooting method have a sharp cut off and always seems to sum to the same ratio, whereas that doesn't happen in the matrix method? Which ion is taking up the slack when Fe I is not 100%??

[jhmatthews]

Knox, did you merge this / do you want to? And do we understand my question about the rate matrix graph with Fe I?

[kslong]

No, I did not merge it. Chicken I guess. And yes please do so.

And no, I don’t remember the question

Knox

[jhmatthews]

I mean the one in my previous comment:

"I'm a little confused as to the shape of the rate matrix graph off to the left- why does the shooting method have a sharp cut off and always seems to sum to the same ratio, whereas that doesn't happen in the matrix method? Which ion is taking up the slack when Fe I is not 100%??"

[jhmatthews]

@kslong - This work is in the matrix_ksl branch. You should be able to see what would happen if we merged matrix_ksl into dev by going to the following link:

https://github.com/agnwinds/python/compare/dev...matrix_ksl

[kslong]

James will merge and check that my update works. We may want to ultimately delete this method. James said one could not quite reproduce the Lucy and Mazzlli with the rate matrix. James said it was because of the assumption about the recombination.

[jhmatthews]

I've tested this for thin shells and the results look sensible as Knox found above. However, it seems to produce errors in examples/core/cv_standard.pf. Namely, lots of errors of this type:

ionization_on_the_spot: nebular_concentrations failed to converge

that are not produced in the current dev version. This is the error that is produced when the code fails to produce a change in ne of less than FRACTIONAL_ERROR inside the maximum number of iterations (MAXITERATIONS). Increasing MAXITERATIONS to the dev value of 200 (Knox had set it to 100) does not help.

Unsure why this is happening but it seems this method is not more stable in current form.

[kslong]

@jhmatthews I am not sure what the path forward is here. Do you have a .pf file that illustrates any remaining problems. Or should we bite the bullet, and remove some of our ionization approaches from python, since we don't use them and would not recommend others using them. We really need to either include .pf files for all our problems, or put them somewhere and reference them from github.

[jhmatthews]

@kslong I always try and use a standard pf file if possible, and simply give the location in the Python directory as I did above.

I think I would advocate the following way forward here, which is all linked to ionization modes too (see #215):

• Keep the source files for the old method in the code
• Make the shooting method (modes 3,6,7 currently?) only accessible via advanced mode or a compiler flag
• Add in some error logic to stop a user using these modes
• This allows us to test but means we use matrix methods where possible
• Consider changing ionization mode to a string for clarity

Questions to answer include how we want the LTE modes etc to work.

[kslong]

@jhmatthews I think we/I dropped the ball on this one in terms of getting it into dev But perhaps we should try to recover.

[kslong]

The matrix_ksl branch still exists, but this version of the routine is so old it requires the fitsio library to make the code.

The problem we though we were having was that when you had a lot of ions for an element, and the saha factors between adjacent ions was too great that one could end up with numbers whose exponents were out of range for the computer. We originally limited this to 10**20 which is too large. We set the maximum difference between adjacent ions to

big=pow(10.,250.)/(last-first)

where last and first refer to the first and last ion number. In principle this should address the original problem. It is is different from what was there originally ), although ksl does not know if it actually solved the problem.

The version of saha that exists in the matrix_ksl branch is different. For any element, it first identifies the ion that is likely to be dominant (or near dominant), given the temperature, and then works in two directions (one toward lower ionization states and the other toward higher ionization states. Once you get below a certain density it stops. This avoids the problem caused by the accumulation of many powers of 10, and is somewhat more "elegant".

However, what is needed is a robust test of the code for this problem. Before we decide there is really anything that needs to be done, since these modes really are mostly out of date at this point.

[jhmatthews]

This is a very seldom used mode for calculating ionization states and should only be used in test cases (for example, when comparing the SN model to tardis). As a result this can be closed with caveat, but this may pop up if someone cares for some reason about Fe in, e.g., the Lucy Mazzali scheme.