Trapped photon in macro atom mode

[kslong]

I am encountering a trapped photon problem with the current dev version of the code for the following .pf file:

warm_macro.pf.txt

The problem arises during the detailed spectral stage, on the very first photon one is trying to generate (namely 66637). The trapping is occurring because this loop is never being escaped:

in photo_gen_matom.c

The reason the loop is never escaped is because in the routine kpkt, there are a set of options to determine what kind of photon is being produced. The choice puts the photon into this region of the code:

and *escape sets to 0, and you return. At that point, i is always 0 (nres 576) which is likely Ly-alpha. This triggers a repeat.

mplasma is as above.

If I take this at face value, it looks like we want to generate a line photons, but we cannot because all the line cooling is due to a macro atom, and we are not allowed to generate a line photon for a macro atp,

Note that the logic of the kpkt routine is fairly obscure. There seems to have been a choice made to avoid generating more than one random number. There are various confusing places where one subtracts from the destruction_choice and then adds back to it.

Note: It is possible that this is a layered problem in the sense that there are other errors that are putting us into this situation, but it still seems to me one should be able to extract oneself from this problem, given what is here.

@jhmatthews I think you are probably going to have to solve (or at least comment on this).

[kslong]

The problem also arises in the same model with cylindrical coordinates. Like the routine with spherical coordinates this routine also produces other errors, that may be related specifically, the routine produces errors, like
cdf_gen_from_array: all y's were zero or xmin xmax out of range of array x-- returning uniform distribution 0

That is coming from a call from one_fb, namely

fb (xplasma, xplasma->t_e, fb_x[nnn], nions, FB_FULL)

where this is always returning 0. It's hard to know how related this is to our bug, as this is associated with fb, whereas the problem looks to be due to bb. It's possible though that fb_tot should be a lot larger than it appears in mplasma above.

[jhmatthews]

Ok, this is a tricky one. I agree that the procedure here is a little convoluted. photo_gen_kpkt wants to produce photons from the k->r process, but not photons where a kpkt excites a macro-atom (i.e. k->A*->r or similar in Lucy's notation). Hence we keep creating kpkts and just wait until a k->r process happens, as Knox says.

It's not clear to me that this is just due to the section of kpkt regarding lines - the bf cooling rate is comparable to the bb cooling rate so just because it was a line process this time doesn't mean it always will be.

I think fundamentally the issue is that the kpkt emissivity is non-zero, but the chances of a k->r packet lying in the wavelength range 6500->6600 are absolutely tiny. But I've a feeling there might be other issues here too. E.g., should the emissivity be non-zero? Regardless, clearly the code should be able to deal with this situation in a nicer way.

[jhmatthews]

First observation is that if you increase the wavelength range to 3000-6500 this problem goes away, which implies my above guess is right.

[kslong]

I have run ito another example of this in trying to look at #290

original.pf.txt

This .pf file takes a good while to run, unfortunately

@jhmatthews What to we want to happen in this case?

[kslong]

This error no longer occurs with the two examples we have used in the past, and so I am closing this issue.

It should be noted that for original.pf, the creation of the detailed spectra can much longer than the the ionization cycles.

Specifically, for original.pf the ionization cycles running with 16 processors take 45 minutes, while the detailed spectral runs minutes 132 minutes. It would be interesting to try to speed up the spectral runs.

Here is the updated .pf file for the longer of the two runs
original.pf.txt

This is an agn type model where we are trying to get the detailed spectrum from 3800 to 5500 A running 10**8 photons per cycle, which undoubtedly is part of the reason the detailed spectral cycles are slow.