Bug in Tardis atomic / standard data set?

[ssim]

It looks to me like the set of atomic data we currently are providing has missing basic ionizating data for Z=4 (Be). This data is not actually important for current applications (abundances are low); but if running with full ARTIS model this leads to a crash since the (new?) plasma insists that the full data is available (and although negligibly small, the Be abundance is not = 0 in the models).

Data appears to be missing the top ion stage (expected).

@wkerzendorf can you easily regenerate our standard atomic dataset will all ions for Be and replace the existing one?

[aoifeboyle]

Hi @ssim. Yeah, the missing Be data has been a problem for ages now. I think we agreed that rebuilding the atomic dataset was complicated.

In my own version of Tardis I have a temporary solution that manually adds the Be IV ionisation data. It's not a long-term solution but it does work. I'll make a PR with it (not to be merged ever, just as a temporary solution).

[unoebauer]

Just realised that this issue is quite severe because it breaks a key feature of Tardis, namely setting up a stratified abundance structure (using the abundances config parameters type: file, filetype: simple_ascii). Pushing up the priority.

[unoebauer]

@ssim, @aoifeboyle, @wkerzendorf - maybe we should, for the time being, merge the temporary Be fix (PR #439) to restore the possibility of calculating stratified models. If desired, we could add a warning message, notifying the user, once Be is present in the model, that the Be treatment is incomplete in the current version and that this issue will be resolved in future releases.

The other possibility is of course, as @ssim suggested, to re-generate the atomic data file, now containing the missing data...(However, I have the feeling that this may turn out to be a bit tricky, given the obscurity of the involved scripts).

[unoebauer]

@ssim, @aoifeboyle, @wkerzendorf: I have to go back on my statement from before: the Be issue does not break the generic Tardis capability of handling stratified abundances, as long as all Be abundances are set to zero.

[chvogl]

@unoebauer: As a temporary solution it is possible to manually add the missing Be data to the existing atomic data file:

with h5py.File('kurucz_cd23_chianti_H_He.h5') as atom_data:
    new_ion_data = np.insert(atom_data['ionization_data'][:],9,(4,4,217.71865))
    del atom_data['ionization_data’]
    atom_data.create_dataset(u'ionization_data', data=new_ion_data)

[unoebauer]

Would be good if we could do something about this for v1.5 - maybe just adopt @chvogl's quickfix and update the atomic database shipped with tardis_example?

[unoebauer]

Important Update: This fix will work only for LTE setups.

The quickfix only adds the missing ionization data. The modified atomic data set still lacks levels and transitions data for the ions (4, 3), (4,4), i.e. Be IV and Be V.

This leads to a shape mismatch problem in the plasma calculation in the nebular ionisation case. The partition function uses the available levels to set its shape and the nebular deltas rely on the dimensions of the ionization data. In my stratified W7 test case this leads to the following situation:

/home/ulrich/python-virtualenv/tardis/lib/python2.7/site-packages/tardis_sn-1.5.dev2074-py2.7-linux-x86_64.egg/tardis/plasma/properties/ion_population.pyc in calculate(t_rad, w, zeta_data, t_electrons, delta, g_electron, beta_rad, partition_function, ionization_data)
     91             partition_function, ionization_data)
     92         zeta = PhiSahaNebular.get_zeta_values(zeta_data, phi_lte.index, t_rad)
---> 93         phis = phi_lte * w * ((zeta * delta) + w * (1 - zeta)) * \
     94                (t_electrons/t_rad) ** .5
     95         return phis

ValueError: operands could not be broadcast together with shapes (462,27) (464,27) 

In [18]: %debug
> /home/ulrich/python-virtualenv/tardis/lib/python2.7/site-packages/tardis_sn-1.5.dev2074-py2.7-linux-x86_64.egg/tardis/plasma/properties/ion_population.py(93)calculate()
     92         zeta = PhiSahaNebular.get_zeta_values(zeta_data, phi_lte.index, t_rad)
---> 93         phis = phi_lte * w * ((zeta * delta) + w * (1 - zeta)) * \
     94                (t_electrons/t_rad) ** .5

ipdb> p delta.shape
(464, 27)
ipdb> p phi_lte.shape
(462, 27)
ipdb> 

Bottom line: unfortunately, the quickfix does not solve the problem. Any suggestions, @ssim, @aoifeboyle, @wkerzendorf, @chvogl?

[ssim]

Very quick thoughtette: probably the way I would try to fix this is just to have a catch such that if the data do not exist it assumes 1.0 for zeta and carries on / logs a warning rather than throws a wobbly. Can that be done easily?

[unoebauer]

@ssim - don't think this will work. The problem is not with missing zeta data (which is actually already done the way you are suggesting) but in the design of the plasma calculation. Well that sounds as if something is wrong with the plasma part - which is not true. Rather the problem comes from the quickfix, as suggested above, only introducing part of the missing data. The plasma calculation just assumes that for each of the ionisation stages which appear in the ionization_data part of the atomic_data set, corresponding levels_data is also available - which is not the case. This then leads to different shapes for different quantities in the plasma calculation and the run crashing.

A viable quickfix would be, however, to check whether Be is present in the model. In case it is, Tardis logs a warning and just sets the abundance to zero everywhere.

[aoifeboyle]

@unoebauer You mean automatically setting the abundance of Be to zero? This problem is difficult to deal with properly in any other way. The problem with the quick fixes that I have tried is that it's hard to insert data into multiindex dataframes, and the method is not consistent across the recent pandas versions. So a fix that may work on my laptop in that anaconda environment won't work on my desktop, or stops working if a package is updated. @ssim I'm doing the final run of my paper simulations now, and the low mass model results gives the exact same results as I got from the old plasma (which is a reassuring result for the new plasma), but the high mass model does not (I'm pretty sure this is because of the Be, which the high mass model has and the low mass model does not). Since the high mass model only has negligible quantities of Be, I've just removed it for now. Is that ok with you?

[ssim]

@aoifeboyle @unoebauer As a matter of pragmatism, yes just switching Be out (i.e. setting abundance to zero) is a fine solution, just not very elegant. For purposes of the paper, however, it should be ok.

[unoebauer]

After some back-and-forth discussion with @wkerzendorf, the consensus is now that we should avoid a dirty quickfix (c.f. PR #502) but rather migrate the full Be data into our standard atomic_data file.

@wkerzendorf will update the data base shipped with tardis_example until the v1.5 release.

[unoebauer]

This should be fixed with the carsus project - something for @yeganer to explicitly check

[unoebauer]

Finally, we can safely close this issue. The new standard atomic data set, produced with carsus and located in the tardis-refdata repository has the full Be data included.