macro atoms - one free free errors

[kslong]

We are seeing errors of the type

one_ff: Bad inputs f2 2.49313e+14 < f1 2.49827e+14 returning 0.0 t_e 119.656

in the C345 model with macro atoms.

This is being caused by the following code

and is called via through the following routiens

in macro_gov.

Note that the lower limit of the frequency is set by the bands, but the upper limit is set by the temperature of the grid cell. This is a logical problem.

We have tried to fix problems with this various times in the past, see #187 and #300, so obviously this is a touchy problem.

[kslong]

Part of the problem here is that electron temperatures vary widely, and some cells can get quite cold, but we still need photons in the hot (inner) part of the wind.

A couple of other things about this error. one_ff returns -1 as the frequency when it hits this problem, and an error code. I don't believe the calling routine checks for either of these cases.

The other question here is why is one_ff being called at all. If we are going to have the limit the way it was, then we should check that when total ff is calculated it returns 0, so one_ff is not called.

Finally, a some disconnected point, which may be irrelevant. I have been trying to set TMIN to 100 at various places in the code. We probably need to be consistent if possible.

[jhmatthews]

This error gets hit when the temperature is 50 Kelvin. Do we now why it gets quite that cold?

I can force freqmax to be higher than freqmin by some factor, but there could feasibly be numerical issues then (there's a reason why a maximum h nu / kT was put in). More generally it would seem to me that we shouldn't be allowing electron temperatures to get lower than (H * freqmin / k_B), or 100K would also do. freqmin being the global 1e14Hz limit.

[jhmatthews]

Sorry Knox, didn't read the initial statement correctly with the frequency and temperature quoted. Generally what I said above still applies though.

[kslong]

James, we should think about the location for some of this code. I assume the reason for the modifications made previously were to prevent some kind of exponent underflow (as your comments above suggest), but they ought to affect the simple case as well, suggesting that any kind of extra limitations on the frequencies should be made in one_ff, and not in the matom sections.

I also think that the reason this came up in the first place is associated with a disconnect between the way total_ff is calculated.

[kslong]

So James, I am a bit confused about the ff problem after looking at it a bit, and attempting a fix which failed.

In matom, around line 939, we call total_free, and use frequency limits between 0 and VERYBIG to calculate the ff emissivity, which is then used as part of the process to establish how often a kpkt is converted to a ff photon.

But later on (line 1175) we force the free free photon that is emitted to fall between between xband.f1[0] and a frequency determined by t_e (in the ionization cycles).

What do we want here? Should we be setting the total ff emissivity to 0, if the corrected upper limit is less than xban.f1[0], which would avoid the possibility of ever getting to one_ff, and should we be replacing VERYBIG with the "temperature corrected" upper limit? Or alternatively should we be setting one_ff so it has a lower limit of 0, thereby generating some very long wavelength photons.

[jhmatthews]

I think either of those would be ok except having very long wavelength photons. For the ionization cycles, I think that the call to total_free in kpkt should probably use freqmax based on hnu/kT, and if hnu/kT is less than zero then we should just set the free-free cooling channel to zeros.

For spectral cycles, I am unsure about whether what we have now is correct. i.e. is it right that the free-free cooling rate is calculated over the full frequency range, but the one_ff call is over the band limited range (Do you know the answer to this @ssim ?)

Again, I think one way to help is to have a minimum temperature and frequency that are related, i.e. have freqmin = H_OVER_K * NU_MIN, and ensure that the temperature can't drop below this. This would ensure one can't have NU = 100.0 * T_e / H_OVER_K being less than freqmin.

[kslong]

@jhmatthews @ssim I put a band aid on this issue (in the macro_diag branch), by checking whether the upper limit to the allowed frequency range was at least 1.1 freqmin in the ionization phase, even if this was larger than the limit derived from the temperature.

Although this eliminates the issue, I am not really convinced this is what should be done. I am concerned:

1. The emissivity which I assume is proportional to the fraction of the time a k-pkt decays by ff, is always calculated over the frequency range 0 - VERYBIG
2. For very low temperatures, I assume we decided wanted to make sure the photon did not disappear, and so, during ionization cycles, we give it a frequency that presumably is low enough that it will not interact again. But given that we remove photons due to adiabatic cooling, why don't we do that here.
3. And I, like James, don't really understand what is going on in the detailed spectrum case, and why there we always set it to the band limited range.

So Stuart, I would appreciate your comments on these areas, before we close the issue.

[kslong]

According Stuart, some of this historical, because originally we did not have banding. We could not in principle, put this in as a sink of energy, e.g to make it like a photon that was adiabatically eliminated. This would probably be more logically consistent. Alternatively, one could simply set the limit 0.

We talked about the detailed spectrum as well, and concluded what is done is wrong. We either need to calculate the emissivity in the band and create photons only in the band, or we have to use elimination.

[jhmatthews]

The biggest issue now in the code regarding kpkts and free-free concerns the frequency limits and some mistakes in my coding there.

The routine kpkt() should always be selecting the k-packet destruction channel according to the ratio of cooling rates in a cell. It does this correctly by calling total_free with wide frequency ranges. However, when it calls one_ff in the spectral cycles it currently forces the frequency to come out in the spectral range requested, which we think is wrong. It should allow a wide frequency range so that the rejection methods in get_matom_f (to get the right banded k-packet emissivity) and photo_gen_kpkt (to generate the photons with the right relative contributions from each process) work correctly.

There are a few ways one could do this, but I suggest the following:

• kpkt() always calls one_ff with a wide frequency range, say 0 to (100 * k T / h).
• if the photon is less than the global minimum frequency, we discard the photon with an istate set to P_LOFREQ_FF, in the same way we do for adiabatic cooling

An example proposal can be found here: agnwinds/python@dev...jhmatthews:kpkt_ff

[kslong]

I agree that this is a reasonable approach, and have looked at the code that you have written to implement it. I would vote to go ahead with this.

[jhmatthews]

The reason I haven't merged this yet is I'm not totally sure about my original proposed fix. I think there will be problems in the CDF if I adopt 0.0 as the lower limit for generated the photon frequency. Instead, I've proposed a fix that calculates the amount of free-free cooling that comes out above and below the global minimum frequency - one_ff is only called the right fraction of the time to account for these relative cooling rates. I think this is much more sensible.

The kpkt_ff branch has been updated with this: agnwinds/python@dev...jhmatthews:kpkt_ff

I'll do some tests then open a pull request.

[jhmatthews]

The tests generally look good - there are some small changes in temperature structure etc - but other than that ok and spectra identical. There are about 100 errors that occurs in the spectral cycles that I am trying to figure out - it says that it can't find a destruction channel in k-packet. Since I made some changes to this routine I'd like to figure out what the problem is there.

[jhmatthews]

I still don't understand these issues. I've documented the tests of a few different commits to branch jhmatthews/kpkt_ff on evernote under the page: "Python - 2019 - Testing of changes to free-free and macro-atom level pops". There are still multiple errors, some of which are expected and some of which appear critical. The spectra on the whole look OK.

I believe the main fundamental problem is still this error:
Error: matom.c: Failed to select a destruction process in kpkt. Abort.

Which should never happen, but happens about five times in the spectral cycles. Debugging the problem is difficult because the restarted runs don't display the same error. It looks like the restarts have different quantities in kpkt_abs, which implies yet another bug.

[jhmatthews]

The above error that occurs in my kpkt_ff branch disappears if you set kpkt_rates_known to -1 (i.e. force recalculation of the kpkt rates) at the start of the spectral cycles. So I could in principle put this in and merge this, but it's not clear to me why forcing the recalculation of kpacket cooling rates helps this problem go away. Hints at an underlying problem IMO.

[kslong]

@jhmatthews What is the status of this? I am now encountering errors that cause the Python to crash in attempting to run the following parameter file. What is happening (on dev) is that it is attempting to create a cdf associated with free-free when f1 = f2. This is happening at the beginning of detailed spectrum creation. If I make a small mod to make the cdf generation routine return an error (and not exit) when this condition is satisfied, I eventually fail with too many errors of the same type.

hot.pf.txt

[kslong]

Note that error associated with with hot.pf described above is not the same one that James was trying to address. Furthermore, I ran hot.pf with the version of python James made on the kpkt_ff branch, and hot.pf fails in exactly the same way as it did on the dev branch. So these may be unrelated issues, or they may both point to an underlying problem that we have not fully understood yet.

[kslong]

#511, which deals with the problem related to hot.pf is closed. In the end, that problem was only minimally associated with the one described here. @jhmatthews I still need you finish this one off.