Negative jbar contribution in macro-atom run

[jhmatthews]

This is an error report documenting one of the (many) issues that have recently come up with macro-atoms (see #425), whereby the code tries to add a negative contribution to the jbar estimator in bb_estimators_increment().

To reproduce:
Commit: agnwinds/python@39f1517
Pf file: h10_standard80.pf.txt
I find this fails in the first cycle if I run with mpirun -n 3 py h10_standard80, with the error

Negative contribution to jbar. Abort

My current understanding of this problem is:

• the negative contribution to jbar arises because of a negative photon weight
• The negative photon weight arises because of our new "upweighting scheme", which multiplies the photon weight by (1- nu_0/nu) around line 1210 in resonate.c. For some reason, nu is actually less than nu_0, i.e. a photon with a frequency less than
• this results in the photon weight being multiplied by a -ve number.

The reason nu is less than nu_0 appears to be due to Doppler shifts. When we calculate the bound-free opacity in calculate_ds, we use "freq_inner" which I think is the photon frequency shifted to the frame of the plasma at the start of the path element. This is used to get the opacity and is higher than the frequency that we use

Question: is p->freq always the frequency of photon in the observer frame? @kslong @Higginbottom @ssim

If so, there's an error here, and there possibly always has been, in how we calculate bound-free opacities and heating in macro-atom mode. It seems we came across a similar issue before in simple-atoms (see #73 ) - I'm betting this is the issue.

[jhmatthews]

Proposed fix for this individual issue opened as a pull request in #437

I've carried out a quick review of the various places that estimators get incremented and so on to see if this issue could occur anywhere else. This file contains all the instances I can find of where photon frequencies are used in a context where a possible error could happen:
frequencies.txt - this only accounts for when the "freq" is in a photon pointer so there could be others where we pass around a double.

Having had a look through, it seems we do things more or less correctly in radiation.c (in that we choose an average comoving frame frequency) and similarly we do things right when we calculate opacities. I'm assuming that all the emission and matom deactivations bits are done ok - it should choose a comoving frequency then shift it to observer with doppler() or something. I think there is one other main error:
In indivisible packet mode, there is no shift in bf_estimators_increment, which is kind of the analogue of radiation.c, so I think the bound-free estimators will be a little wrong, as will xplasma->ave_freq and xplasma->j, and the ionization parameter diagnostics.

We could also improve some stuff codewise:

• we have a bit of code duplication in how we calculate doppler shift
• probably more importantly: I think it would probably make sense to have two separate variables in the photon structure: observer and comoving frame frequency. One could then do away with most of the issues here and minimise human error. One could either update comoving frame frequency every time a photon is moved, but if this is too expensive one could have a function called update_comoving_frequency() that does it at less frequent intervals when needed. However, this is a change we could easily screw stuff up so caution needed as usual...

[jhmatthews]

I note that having coded up the solution to this (#437), there are now a lot of my errors being reported, whereby the frequency is just below threshold. I think this is probably because when one calculates the opacity one uses some average comoving frequency along the path element, and this can be different to the actual comoving frequency when the scatter event is identified.

Probably this doesn't matter at all, because the heating effect of these photons is by definition small - so we just need to make a decision on what is done.

[kslong]

Since we know what the frequency is at the interaction point, why can’t we fix this properly?

[jhmatthews]

The point is that although the frequency of the photon is known at the interaction point, the frequency used to calculate the opacity in the edge is an average frequency along some small path element. This means it is possible for the frequency of the photon at the point of interaction to be lower than the edge frequency. So a decision need to be made about what to do.

[kslong]

@jhmatthews Is there anything more to do here, or should this be closed, with perhaps additional documentation.

[jhmatthews]

I've added the final fix here, which is only report an error if the variable prob_kpkt is above 1e-3 in magnitude. This could be improved by using information about the velocity gradient and path length, but I'm closing the issue since we've made the major fixes here and in any case prob_kpkt gets set to 0 if negative.