Virtual Spectrum Zero

[unoebauer]

I noticed this while investigating issue #191: using the current upstream/master 34f667a, Tardis does not produce a virtual spectrum any more. The following model attributes are all arrays populated exclusively with zeros:

• mdl.montecarlo_virtual_luminosity
• mdl.spectrum_virtual.lumosity_density_nu
• mdl.spectrum_virtual.lumosity_density_lambda

This is a serious issue. It is also worrisome that none of our automated tests picked this one up. @wkerzendorf: we should include a test_virtual_spectrum routine into the test_tardis_full script.

[unoebauer]

Seems to be related to the virtual packet logging scheme introduced in PR #433. When running with Tardis complied with the flag --with-vpacket-logging, the virtual spectrum is properly and correctly calculated. Once this flag is deactivated, I get zeros.

This is clearly not what we intended with PR #433. Here, only the additional memory intensive virtual_packet diagnostic arrays should be activated or de-activated. The virtual spectrum, however, should be calculated in any case (correct me if I am wrong, @wkerzendorf, @orbitfold)

I'll have a closer look at this shortly and will try to fix it.

[orbitfold]

If I understand it correctly then this is "correct" behavior. It should not be stored when --with-vpacket-logging is not enabled. Virtual packets are still used, just the spectrum is not stored. I could have gotten this wrong though.

[unoebauer]

The only reason to use virtual packets is to increase the signal-to-noise in the spectrum. So, there would be not much point to use virtual packets, which slow down the calculation, and then not benefit from the only advantage that they offer, namely decreasing the MC noise in the spectrum.

So, whenever virtual packets are active in the simulation, we definitely want to calculate/store the spectrum. However, we not necessary need always to keep track of all the interaction details of the virtual packets. That said, your switch, which is very useful, should only activate and deactivate the storage of arrays such as

• runner.virt_packet_last_line_interaction_in_id
• runner.virt_packet_last_line_interaction_out_id
  ...etc

but should not interfere with calculating and storing the virtual spectrum. Again, correct me if I am wrong, @ssim, @wkerzendorf.

[orbitfold]

Well the spectrum is still calculated and stored. The reason for storing the spectrum for virtual packets alone was to use the information in uncertainty calculations. It wasn't stored before we needed it this summer, I believe. Which is also why tests didn't fail. The full spectrum test uses virtual packets. It would have failed if they stopped working.

[unoebauer]

Seems to be not the case, as the following comparisons show:

Both plots were produced with the following instructions:

config = yaml.safe_load(open('tardis_example.yml'))
config["montecarlo"]["iterations"] = 1 #only to speed-up the calculation
mdl = run_tardis(config)

plt.plot(mdl.spectrum.wavelength, mdl.spectrum.luminosity_density_lambda)
plt.plot(mdl.spectrum_virtual.wavelength, mdl.spectrum_virtual.luminosity_density_lambda

• results obtained with Tardis compiled with vpacket-logging (real packet spectrum blue, virtual packet spectrum green):
  
• results obtained with Tardis compiled without vpacket-logging (real packet spectrum blue, virtual packet spectrum green):
  

[unoebauer]

Found the problem, I think - see the inline note I just added to the PR #433

[unoebauer]

The virtual spectrum is calculated in simulation/base.py lines 318-323, using the runner attribute legacy_montecarlo_virtual_luminosity. This is linked to storage.spectrum_virt_nu (see montecarlo.pyx, lines 197,198).

Whether this array is populated (cmontecarlo.c, lines 451-456) is currently subject to a large ifdef block and is thus only performed if vpacket loggin is active. I'll open a PR to fix this shortly.