Spectral estimators and banding

[kslong]

The purpose of banding in photon generation was to improve ionization calculations by assuring that there were enough photons in wavelength regimes far from those that dominated the luminosity. A flexible mechanism was created to to this. However, when we decided to create spectral models of the spectrum in each cell, we did not couple the photons that were generated to the bands used for photon generation; instead we used hardwired bands. We are now encountering problems associated with this disconnect that need to be resolved.

(I'm creating this issue in order to allow us to determine a solution to this problem.)

[kslong]

Following a discussion of what we might do, I have implemented in bug798_bands a change which forces the bands used in the ionization calculation to be the same as used for photon generation. While the switch does not change the regression test models very much, I am hoping others will test this on models most relevant to the science they are doing to see if there are any obvious problems.

[kslong]

@jhmatthews

I have implemented a version python bug798_bands that calculates spectra in a cell. However, am am having trouble getting multi-threaded operations to work properly. The cell spectra are created as photons transit the wind. They then need to be renormalized to reflect the cell volume. This works properly for one thread, but not for multiple threads.

The cell spectra from the various threads are gathered together in para_update.c in the routine gather_spectra. However if you print out the values for one of the spectral elements from each thread before trying to sum them up and braodcasting them back to the various threads, you see that there is one thread has the correct value (about what you get with a single thread) and the rest have very much larger values. The larger value reflects what you would expect if the spectra had not been renormalized.

Here is an example. The two columns are the of the element 500 of the cell 15 spectra before gather spectra has done it's work. The second column is the first divided by the number of threads foo.pf has two cycles so you see two lines from each thread.

Note that the values for thread 4 are much lower than the rest of the thread. Thread 4 has a similar value to what you would get for a single thread. The larger numbers are all what you would expect if the spectra were not renormalized by the volume.

foo_0.diag:XXX1 2.966263e+48   3.707829e+47  
foo_0.diag:XXX1 3.982022e+48   4.977527e+47  
foo_0.diag:XXX1 2.966263e+48   3.707829e+47  
foo_0.diag:XXX1 3.982022e+48   4.977527e+47  
foo_1.diag:XXX1 5.116261e+48   6.395327e+47  
foo_1.diag:XXX1 3.585864e+48   4.482330e+47  
foo_2.diag:XXX1 4.531289e+48   5.664111e+47  
foo_2.diag:XXX1 3.569332e+48   4.461666e+47  
foo_3.diag:XXX1 4.105360e+48   5.131700e+47  
foo_3.diag:XXX1 3.344136e+48   4.180169e+47  
foo_4.diag:XXX1 3.038418e+09   3.798023e+08  
foo_4.diag:XXX1 3.359165e+09   4.198956e+08  
foo_5.diag:XXX1 3.293677e+48   4.117096e+47  
foo_5.diag:XXX1 2.304753e+48   2.880941e+47  
foo_6.diag:XXX1 3.533117e+48   4.416397e+47  
foo_6.diag:XXX1 4.406299e+48   5.507873e+47  
foo_7.diag:XXX1 3.962234e+48   4.952792e+47  
foo_7.diag:XXX1 4.959145e+48   6.198931e+47  

So the problem is somehow related to the fact that gather_spectra does not "receive" the renormlized cell spectra except for one of the threads. I though MPI_BARRIER was supposed to prevent that. Help is needed.

Note that there various XXX comments that print out diagnostics where the cell spectra are created and renomalized (which is in estimators_simple.c)

foo.pf.txt

[Higginbottom]

I have run a simple stellar wind model with a few cells - could be better to get a really spiky spectrum - but as a first look see Blue line is KSLs new detailed spectrum - orange line is the modelled cell spec...

If the detailed spectrum is intended to be J_nu, there is a normalisation issue... will investigate... (definately not volume - that is around 1e35...

OK, just needs dividing by the width of the frequency bin used to generate the fine spectrum....

[kslong]

That looks great @Higginbottom Do you think we should do this in the program or should it be part of post-processing.

Pros for in the program

• Reduces confusion and makes what is saved have units we expect them to have.

Cons for in the program

• The only con I can think of really would be if we decided to get clever and accumulate spectra over multiple cycles, and had to think about what we do on restarts. (I do keep wondering whether when we have too few photons in a given cell within a cycle, whether we would be better off keeping photons from multiple cycles, down-weighting those from previous cycles. This would probably generate a slower but more accurate path to convergence.)

We should make sure to document what these "spectra" actually are.

[Higginbottom]

I think how we deal with the spectrum depends on what we want to do with it.

In principle, we can use this spectrum (once normalised) to calculate the PI rates in a slightly more flexible manner than using the models. If it's not significantly slower than using the models (which it might be) then this could be great. It could be an extra ionization mode - but actually if it isn't any slower than using models, it should replace it.

Or, we can use the detailed spectrum to compute models - which was I think the initial idea. I suspect this will be more complicated (and perhaps slower) than just using the spectrum directly when needed!

The idea of collecting data across cycles is really interesting - I'm really not sure about what is best. Certainly right at the start it might bot be a great idea since the ionization structure is way off - but if it makes the code more stable it might still be a good idea...

Once scaled by bin width - I believe this spectrum is J_nu - erg/s/cm^3/Sr/Hz

[kslong]

@Higginbottom Can you just check that the cell spectra are properly normalized in the latest dev version before I close this issue (for now). We will leave further improvements to the future.

If so, I will change the version number to indicate we have finished coding for the year and merge dev into master.

[Higginbottom]

Here is a plot of the raw spectrum as output from windsave2table vs the model - normalising looks fine.

Just out of interest - here is the plot with 100x fewer photons...

In this case, the detailed spectrum certainly seems to do a better job at the edges - although the noise would require some care

And here it is for really tiny numbers - interesting...

The point I was trying to suggest (and I might be totally wrong here) is that there is a component of j_nu in a cell which is due to the cells 'own' radiation - recombination, free-free etc. In the past we have had problems when the number of wind photons being made is rather small, so the ionization state of cells that dont get many external photons is unhealthily dependant on very small numbers of 'self ionising' photons made in the cell. If we were to use these j_nu spectra we could include a contribution from the cells own emission spectrum, by generating that spectrum and adding it into the logged j_nu spectrum. This would essentially be the 'diffuse' continuum in the cell. Not sure if it makes sense as a concept - but it seems to be an idea worth discussing...

[kslong]

@Higginbottom Thanks for doing checking the normalization. This will allow me to close thing out for the year.

The idea of thinking about the difference between the "external" radiation and the "internal" radiation is interesting.

Do you recall if we have a good 1d example of the problem we would like to solve, or do we require a 2d wind so that shielding we need happens?

[jhmatthews]

There is some good stuff here. Nick makes an interesting point, that recording, in an "estimator" sense, the J_nu contribution from the cell's own radiation is in some sense inefficient, because in principle you know the radiation the cell should emit and you could more accurately characterise the "internal" J_nu. We probably don't have the appetite for implementing this, but a very interesting related test to do would be to take a slab (or thin shell) with a high tau, and see what emergent spectrum (and temperature structure) it has with number of cells N ranging from 1, to a few times tau. It will be interesting how high N needs to be, with different treatments of radiative transfer to get a reasonable answer or "converge".