Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Inconsistencies between Tardis versions #455

Closed
unoebauer opened this issue Dec 15, 2015 · 69 comments
Closed

Inconsistencies between Tardis versions #455

unoebauer opened this issue Dec 15, 2015 · 69 comments

Comments

@unoebauer
Copy link
Contributor

Main description:

TARDIS produces different output between 45eca24 and master with this YML file. Why? and what commit introduced this change.

The following YML files provide a minimum "working" example": issue_455.zip


The diagnosis of the problem at hand is still work in progress, but it seems the most recent Tardis version, c9e1b17, is inconsistent with a previous state of Tardis.

All this has been triggered by re-examining some detailed runs performed during the Abundance Tomography Project of Talytha. For that, we mainly used a Tardis version similar to aa7da17 (a working Open-MP implementation has just been incorporated, PR #338). Comparing the virtual spectra obtained with the recent Tardis version for the same setup, reveals huge discrepancies:

tardis_cmp

The new Tardis version finds a much higher inner temperature. This goes hand in hand with a much higher radiation temperature in the ejecta and a higher ionisation state (and in turn electron density).

ionisation_cmp

Setup:

  • tardis_00173_2.yml
atom_data: kurucz_cd23_chianti_H_He.h5
model:
  abundances:
    filename: abundances_00173_2.dat
    filetype: simple_ascii
    type: file
  structure:
    filename: densities_00173_2.dat
    filetype: simple_ascii
    type: file
    v_inner_boundary: 12000.000 km/s
    v_outer_boundary: 35000.000 km/s
montecarlo:
  black_body_sampling:
    num: 1000000
    start: 1 angstrom
    stop: 1000000 angstrom
  convergence_criteria:
    damping_constant: 1.0
    fraction: 0.8
    hold: 3
    t_inner:
      damping_constant: 1.0
    threshold: 0.05
    type: specific
  iterations: 20
  last_no_of_packets: 500000
  no_of_packets: 50000
  no_of_virtual_packets: 3
  nthreads: 16
  seed: 23111963
plasma:
  disable_electron_scattering: false
  excitation: dilute-lte
  ionization: nebular
  line_interaction_type: downbranch
  radiative_rates_type: dilute-blackbody
spectrum:
  num: 10000
  start: 500 angstrom
  stop: 20000 angstrom
supernova:
  luminosity_requested: 8.675 log_lsun
  luminosity_wavelength_end: 6500.000 angstrom
  luminosity_wavelength_start: 3500.000 angstrom
  time_explosion: 6.800 day
tardis_config_version: v1.0
  • densities_00173_2.dat
6.800000 day
#Index     Velocities [km/s] Densities [g/cm^3]
  0          1.2000e+04          9.652102e-13
  1          1.2408e+04          3.659182e-13
  2          1.2844e+04          3.932164e-13
  3          1.3312e+04          3.900119e-13
  4          1.3816e+04          3.090009e-13
  5          1.4359e+04          3.108902e-13
  6          1.4947e+04          2.463674e-13
  7          1.5584e+04          5.468102e-14
  8          1.6279e+04          3.228778e-14
  9          1.7039e+04          2.211246e-14
 10          1.7872e+04          1.542880e-14
 11          1.8792e+04          1.264987e-14
 12          1.9811e+04          8.297541e-15
 13          2.0948e+04          4.517051e-15
 14          2.2222e+04          5.480575e-16
 15          2.3662e+04          8.695515e-18
 16          2.5301e+04          4.354394e-18
 17          2.7184e+04          2.150320e-18
 18          2.9371e+04          1.040788e-18
 19          3.1939e+04          4.580299e-19
 20          3.5000e+04          1.893468e-19
  • abundances_00173_2.dat
#Index Z=1 - Z=32
 0 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 6.9280e-06 6.2377e-01 5.0460e-09 1.2570e-05 2.9230e-07 6.0000e-02 5.5630e-07 1.8000e-01 5.8250e-08 3.0000e-02 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-06 3.1720e-09 1.6600e-04 1.0820e-07 5.0000e-03 5.2010e-02 4.6180e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
 1 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 6.9280e-06 6.2377e-01 5.0460e-09 1.2570e-05 2.9230e-07 6.0000e-02 5.5630e-07 1.8000e-01 5.8250e-08 3.0000e-02 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-06 3.1720e-09 1.6600e-04 1.0820e-07 5.0000e-03 5.2010e-02 4.6180e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
 2 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 6.9280e-06 6.2377e-01 5.0460e-09 1.2570e-05 2.9230e-07 6.0000e-02 5.5630e-07 1.8000e-01 5.8250e-08 3.0000e-02 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-06 3.1720e-09 1.6600e-04 1.0820e-07 5.0000e-03 5.2010e-02 4.6180e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
 3 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 6.9280e-06 6.2377e-01 5.0460e-09 1.2570e-05 2.9230e-07 6.0000e-02 5.5630e-07 1.8000e-01 5.8250e-08 3.0000e-02 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-06 3.1720e-09 1.6600e-04 1.0820e-07 5.0000e-03 5.2010e-02 4.6180e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
 4 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 6.9280e-06 6.2377e-01 5.0460e-09 1.2570e-05 2.9230e-07 6.0000e-02 5.5630e-07 1.8000e-01 5.8250e-08 3.0000e-02 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-06 3.1720e-09 1.6600e-04 1.0820e-07 5.0000e-03 5.2010e-02 4.6180e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
 5 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 6.9280e-06 6.2377e-01 5.0460e-09 1.2570e-05 2.9230e-07 6.0000e-02 5.5630e-07 1.8000e-01 5.8250e-08 3.0000e-02 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-06 3.1720e-09 1.6600e-04 1.0820e-07 5.0000e-03 5.2010e-02 4.6180e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
 6 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 6.9280e-06 6.2377e-01 5.0460e-09 1.2570e-05 2.9230e-07 6.0000e-02 5.5630e-07 1.8000e-01 5.8250e-08 3.0000e-02 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-06 3.1720e-09 1.6600e-04 1.0820e-07 5.0000e-03 5.2010e-02 4.6180e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
 7 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 2.7585e-02 6.9280e-06 8.8141e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.3220e-02 5.5630e-07 2.3687e-02 5.8250e-08 5.1738e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 2.7997e-07 3.1720e-09 1.3517e-05 1.0820e-07 5.2687e-04 1.2043e-02 2.3331e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
 8 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
 9 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
10 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
11 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
12 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
13 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
14 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
15 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
16 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
17 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
18 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
19 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
20 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
@unoebauer
Copy link
Contributor Author

Tardis version aa7da17 produces compatible results:
tardis_cmp

@unoebauer
Copy link
Contributor Author

Unfortunately, this discrepancy also affects the simpler LTE treatment:

tardis_lte_cmp

@unoebauer unoebauer changed the title Iconsistentcies between Tardis versions Inconsistencies between Tardis versions Dec 15, 2015
@unoebauer
Copy link
Contributor Author

Lacking better ideas, I'll look through the commit history and will try to pin-point it to a specific change. Any other insights into the problem, @aoifeboyle, @wkerzendorf, @ssim, @orbitfold?

@unoebauer
Copy link
Contributor Author

Status of commit-bisection (works implies no discrepancies between :

still bisecting...

Note the following when reading through the above list:

  • works: no difference between virtual spectrum calculated with that version and the original calculation above the Monte Carlo noise level
  • fails: clear differences between the spectra which cannot be attributed to Monte Carlo noise
  • zeta-crash: calculation crashes before reaching the final spectral run because the radiative temperature is higher than the values in the table for the zeta values - see Note 2

Note 1:

---------------------------------------------------------------------------
PlasmaMissingModule                       Traceback (most recent call last)
<ipython-input-2-6282305a6c09> in <module>()
----> 1 mdl = tardis.run_tardis(config)

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/base.pyc in run_tardis(config, atom_data)
     38     tardis_config = config_reader.Configuration.from_config_dict(
     39         config_dict, atom_data=atom_data)
---> 40     radial1d_mdl = model.Radial1DModel(tardis_config)
     41 
     42     simulation.run_radial1d(radial1d_mdl)

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/model.pyc in __init__(self, tardis_config)
    134                                                          excitation_mode=tardis_config.plasma.excitation,
    135                                                          line_interaction_type=tardis_config.plasma.line_interaction_type,
--> 136                                                          link_t_rad_t_electron=0.9)
    137 
    138         self.spectrum = TARDISSpectrum(tardis_config.spectrum.frequency, tardis_config.supernova.distance)

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/plasma/standard_plasmas.pyc in __init__(self, number_densities, atomic_data, time_explosion, t_rad, delta_treatment, nlte_config, ionization_mode, excitation_mode, line_interaction_type, link_t_rad_t_electron)
    103             delta_input=delta_treatment, nlte_species=nlte_config.species,
    104             previous_electron_densities=initial_electron_densities,
--> 105             previous_beta_sobolevs=initial_beta_sobolevs)

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/plasma/base.pyc in __init__(self, plasma_properties, **kwargs)
     18         self.plasma_properties = self._init_properties(plasma_properties,
     19                                                        **kwargs)
---> 20         self._build_graph()
     21 #        self.write_to_tex('Plasma_Graph', 'Plasma_Formulae')
     22         self.update(**kwargs)

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/plasma/base.pyc in _build_graph(self)
     79                                               '{1} which has not been added'
     80                                               ' to this plasma'.format(
---> 81                         plasma_property.name, input))
     82                 try:
     83                     position = self.outputs_dict[input].outputs.index(input)

PlasmaMissingModule: Module PhiSahaNebular requires input t_electron which has not been added to this plasma

Note 2:

ValueError                                Traceback (most recent call last)
<ipython-input-2-6282305a6c09> in <module>()
----> 1 mdl = tardis.run_tardis(config)

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/base.pyc in run_tardis(config, atom_data)
     40     radial1d_mdl = model.Radial1DModel(tardis_config)
     41 
---> 42     simulation.run_radial1d(radial1d_mdl)
     43 
     44     return radial1d_mdl

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/simulation.pyc in run_radial1d(radial1d_model, history_fname)
     25         logger.info('Remaining run %d', radial1d_model.iterations_remaining)
     26         radial1d_model.simulate(update_radiation_field=update_radiation_field, enable_virtual=False, initialize_nlte=initialize_nlte,
---> 27                                 initialize_j_blues=initialize_j_blues)
     28         initialize_j_blues=False
     29         initialize_nlte=False

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/model.pyc in simulate(self, update_radiation_field, enable_virtual, initialize_j_blues, initialize_nlte)
    333 
    334         self.calculate_j_blues(init_detailed_j_blues=initialize_j_blues)
--> 335         self.update_plasmas(initialize_nlte=initialize_nlte)
    336 
    337 

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/model.pyc in update_plasmas(self, initialize_nlte)
    237 
    238         self.plasma_array.update_radiationfield(self.t_rads.value, self.ws, self.j_blues,
--> 239             self.tardis_config.plasma.nlte, initialize_nlte=initialize_nlte, n_e_convergence_threshold=0.05)
    240 
    241         if self.tardis_config.plasma.line_interaction_type in ('downbranch', 'macroatom'):

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/plasma/standard_plasmas.pyc in update_radiationfield(self, t_rad, ws, j_blues, nlte_config, t_electrons, n_e_convergence_threshold, initialize_nlte)
     50         if nlte_config.species:
     51             self.store_previous_properties()
---> 52         self.update(t_rad=t_rad, w=ws, j_blues=j_blues)
     53 
     54     def __init__(self, number_densities, atomic_data, time_explosion,

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/plasma/base.pyc in update(self, **kwargs)
    141 
    142         for module_name in self._resolve_update_list(kwargs.keys()):
--> 143             self.plasma_properties_dict[module_name].update()
    144 
    145     def _update_module_type_str(self):

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/plasma/properties/base.pyc in update(self)
     84         if len(self.outputs) == 1:
     85             setattr(self, self.outputs[0], self.calculate(
---> 86                 *self._get_input_values()))
     87         else:
     88             new_values = self.calculate(*self._get_input_values())

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/plasma/properties/ion_population.pyc in calculate(self, general_phi, t_rad, w, zeta_data, t_electrons, delta)
     65                              '- requested {2}'.format(
     66                 zeta_data.columns.values.min(), zeta_data.columns.values.max(),
---> 67                 t_rad))
     68         phis = general_phi * delta * w * (zeta + w * (1 - zeta)) * \
     69                (t_electrons/t_rad) ** .5

ValueError: t_rads outside of zeta factor interpolation zeta_min=2000.00 zeta_max=40000.00 - requested [ 24829.47853128  29255.41589326  33609.08295046  37345.62521055
  40019.28004965  40953.81576507  38383.82271919  34434.31244872
  32665.86581665  30041.40250129  27808.87307025  25354.81563128
  23264.0689412   21349.97604201  20834.53491125  20280.93368913
  20109.54322745  19707.24624127  19069.69093371  18789.89685196]

@aoifeboyle
Copy link
Contributor

@unoebauer I have no idea. Pinpointing the commit that first created the problem is probably the best thing to do for now. I'm surprised Travis didn't catch this?

@unoebauer
Copy link
Contributor Author

@aoifeboyle: no, unfortunately, Travis did never catch this. However, the problem considered here is much harder than our final spectral test...

@aoifeboyle
Copy link
Contributor

@unoebauer Ah I see. We should try to create some kind of faster test to determine if the code is going wrong without having to run 20 iterations with 50k packets, I think? So we need to find the simplest possible case where the divergence is observed.

@unoebauer
Copy link
Contributor Author

@aoifeboyle: I think that I've pinpointed the regime in the commit history during which the problem was introduced - see list above.

It happens somewhere between PR #360 and PR #365:

Any thoughts, @aoifeboyle?

I'll try to devise a less expensive test which triggers the same problem.

@unoebauer
Copy link
Contributor Author

I've created a "low-resolution" version of the test problem:

  • tardis_test.yml
atom_data: kurucz_cd23_chianti_H_He.h5
model:
  abundances:
    filename: abundances_test.dat
    filetype: simple_ascii
    type: file
  structure:
    filename: densities_test.dat
    filetype: simple_ascii
    type: file
    v_inner_boundary: 12000.000 km/s
    v_outer_boundary: 35000.000 km/s
montecarlo:
  black_body_sampling:
    num: 1000000
    start: 1 angstrom
    stop: 1000000 angstrom
  convergence_criteria:
    damping_constant: 1.0
    fraction: 0.8
    hold: 3
    t_inner:
      damping_constant: 1.0
    threshold: 0.05
    type: specific
  iterations: 15
  last_no_of_packets: 20000
  no_of_packets: 10000
  no_of_virtual_packets: 3
  nthreads: 16
  seed: 23111963
plasma:
  disable_electron_scattering: false
  excitation: dilute-lte
  ionization: nebular
  line_interaction_type: downbranch
  radiative_rates_type: dilute-blackbody
spectrum:
  num: 10000
  start: 500 angstrom
  stop: 20000 angstrom
supernova:
  luminosity_requested: 8.675 log_lsun
  luminosity_wavelength_end: 6500.000 angstrom
  luminosity_wavelength_start: 3500.000 angstrom
  time_explosion: 6.800 day
tardis_config_version: v1.0
  • densities_test.dat
6.800000 day
#Index     Velocities [km/s] Densities [g/cm^3]
  0      1.2000e+04      9.652102e-13
  1      1.2844e+04      3.804829e-13
  2      1.3816e+04      3.465327e-13
  3      1.4947e+04      2.760799e-13
  4      1.6279e+04      4.252074e-14
  5      1.7872e+04      1.846216e-14
  6      1.9811e+04      1.025111e-14
  7      2.2222e+04      2.306574e-15
  8      2.5301e+04      6.245017e-18
  9      2.9371e+04      1.513200e-18
 10      3.5000e+04      3.003242e-19
  • abundances_test.dat
#Index Z=1 - Z=32
  0 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 6.9280e-06 6.2377e-01 5.0460e-09 1.2570e-05 2.9230e-07 6.0000e-02 5.5630e-07 1.8000e-01 5.8250e-08 3.0000e-02 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-06 3.1720e-09 1.6600e-04 1.0820e-07 5.0000e-03 5.2010e-02 4.6180e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
  1 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 6.9280e-06 6.2377e-01 5.0460e-09 1.2570e-05 2.9230e-07 6.0000e-02 5.5630e-07 1.8000e-01 5.8250e-08 3.0000e-02 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-06 3.1720e-09 1.6600e-04 1.0820e-07 5.0000e-03 5.2010e-02 4.6180e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
  2 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 6.9280e-06 6.2377e-01 5.0460e-09 1.2570e-05 2.9230e-07 6.0000e-02 5.5630e-07 1.8000e-01 5.8250e-08 3.0000e-02 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-06 3.1720e-09 1.6600e-04 1.0820e-07 5.0000e-03 5.2010e-02 4.6180e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
  3 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 6.9280e-06 6.2377e-01 5.0460e-09 1.2570e-05 2.9230e-07 6.0000e-02 5.5630e-07 1.8000e-01 5.8250e-08 3.0000e-02 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-06 3.1720e-09 1.6600e-04 1.0820e-07 5.0000e-03 5.2010e-02 4.6180e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
  4 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 2.8896e-02 6.9280e-06 8.9366e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.1471e-02 5.5630e-07 1.6254e-02 5.8250e-08 3.9934e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 1.4488e-07 3.1720e-09 6.2670e-06 1.0820e-07 3.1418e-04 1.0143e-02 2.2244e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
  5 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
  6 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
  7 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
  8 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
  9 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00
 10 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.0000e-02 6.9280e-06 9.0397e-01 5.0460e-09 1.2570e-05 2.9230e-07 2.0000e-02 5.5630e-07 1.0000e-02 5.8250e-08 3.0000e-03 8.2020e-08 7.3410e-07 3.0650e-08 3.0000e-03 4.6460e-10 3.1210e-08 3.1720e-09 1.6600e-07 1.0820e-07 1.3520e-04 8.5440e-03 2.1330e-02 7.2000e-09 1.7370e-08 0.0000e+00 0.0000e+00

It triggers the same problem as described above, but completes 10 times faster than the original one.

tardis_cmp

@aoifeboyle does this help?

@aoifeboyle
Copy link
Contributor

Yes it should do! Was at a conference today but I have time to investigate now. I'll keep you updated.

@ssim
Copy link
Contributor

ssim commented Dec 15, 2015

this all sounds rather worrying!
do we know if this is specific to using readin files for density/composition, or does a problem manifest even for simple input models?

@unoebauer
Copy link
Contributor Author

@ssim: currently trying to isolate the problem. The tardis_example problem doesn't seem to be affected.
tardis_cmp

@ssim
Copy link
Contributor

ssim commented Dec 16, 2015

Ok - perhaps we can take a look in Wurzburg?

In the meantime, good luck Sherlock!

@unoebauer
Copy link
Contributor Author

Issue is not related to OMP: I get the same results when running the tardis_test setup with 16threads or serially, without activating OMP.

@unoebauer
Copy link
Contributor Author

Maybe some progress towards narrowing down the exact cause of the problem:

Three more test setups, derived from tardis_test.yml (detailed setup files will be provided in next comment):

  • tardis_test_unifab:
    • replace stratified abundances structure by uniform composition (reflecting the overall composition of the original test problem)
    • Strong differences persist

tardis_cmp

  • tardis_test_unifab_branch:
    • same as tardis_test_unifab but now also with branch85_w7 density profile (i.e. no structure files are used)
    • Differences are weaker

tardis_cmp

  • tardis_test_unifab_branch_late:
    • same as tardis_test_unifab_branch but with a longer time since explosion
    • Differences are very weak - results almost identical

tardis_cmp

@aoifeboyle
Copy link
Contributor

@unoebauer That's interesting. I was working through the yaml file and changing parameters to see what eliminated the problem, but looks like you beat me to it. Can we safely say the issue is related to the structure file then? Or do you think there's something else involved?

@unoebauer
Copy link
Contributor Author

@aoifeboyle: not completely sure that it's the structure file alone. I have the suspicion that it is related to high densities...but I need to confirm that. I'll update you as soon as I have the results.

@unoebauer
Copy link
Contributor Author

Here are the yaml files for the different runs:

  • tardis_test_unifab.yml
atom_data: kurucz_cd23_chianti_H_He.h5
model:
  abundances:
    type: uniform
    C: 0.01899
    N: 0.00001
    O: 0.80114
    Ne: 0.00001
    Mg: 0.03468
    Si: 0.07239
    S: 0.01291
    Ca: 0.00300
    Cr: 0.00006
    Fe: 0.00192
    Co: 0.02450
    Ni: 0.03045
  structure:
    filename: densities_test.dat
    filetype: simple_ascii
    type: file
    v_inner_boundary: 12000.000 km/s
    v_outer_boundary: 35000.000 km/s
montecarlo:
  black_body_sampling:
    num: 1000000
    start: 1 angstrom
    stop: 1000000 angstrom
  convergence_criteria:
    damping_constant: 1.0
    fraction: 0.8
    hold: 3
    t_inner:
      damping_constant: 1.0
    threshold: 0.05
    type: specific
  iterations: 15
  last_no_of_packets: 20000
  no_of_packets: 10000
  no_of_virtual_packets: 3
  nthreads: 16
  seed: 23111963
plasma:
  disable_electron_scattering: false
  excitation: dilute-lte
  ionization: nebular
  line_interaction_type: downbranch
  radiative_rates_type: dilute-blackbody
spectrum:
  num: 10000
  start: 500 angstrom
  stop: 20000 angstrom
supernova:
  luminosity_requested: 8.675 log_lsun
  luminosity_wavelength_end: 6500.000 angstrom
  luminosity_wavelength_start: 3500.000 angstrom
  time_explosion: 6.800 day
tardis_config_version: v1.0
  • tardis_test_unifab_branch.yml
atom_data: kurucz_cd23_chianti_H_He.h5
model:
  abundances:
    type: uniform
    C: 0.01899
    N: 0.00001
    O: 0.80114
    Ne: 0.00001
    Mg: 0.03468
    Si: 0.07239
    S: 0.01291
    Ca: 0.00300
    Cr: 0.00006
    Fe: 0.00192
    Co: 0.02450
    Ni: 0.03045
  structure:
    type: specific
    velocity:
      start: 12000.000 km/s
      stop: 35000.000 km/s
      num: 10
    density:
      type: branch85_w7
montecarlo:
  black_body_sampling:
    num: 1000000
    start: 1 angstrom
    stop: 1000000 angstrom
  convergence_criteria:
    damping_constant: 1.0
    fraction: 0.8
    hold: 3
    t_inner:
      damping_constant: 1.0
    threshold: 0.05
    type: specific
  iterations: 15
  last_no_of_packets: 20000
  no_of_packets: 10000
  no_of_virtual_packets: 3
  nthreads: 16
  seed: 23111963
plasma:
  disable_electron_scattering: false
  excitation: dilute-lte
  ionization: nebular
  line_interaction_type: downbranch
  radiative_rates_type: dilute-blackbody
spectrum:
  num: 10000
  start: 500 angstrom
  stop: 20000 angstrom
supernova:
  luminosity_requested: 8.675 log_lsun
  luminosity_wavelength_end: 6500.000 angstrom
  luminosity_wavelength_start: 3500.000 angstrom
  time_explosion: 6.800 day
tardis_config_version: v1.0
  • tardis_test_unifab_branch_late.yml
atom_data: kurucz_cd23_chianti_H_He.h5
model:
  abundances:
    type: uniform
    C: 0.01899
    N: 0.00001
    O: 0.80114
    Ne: 0.00001
    Mg: 0.03468
    Si: 0.07239
    S: 0.01291
    Ca: 0.00300
    Cr: 0.00006
    Fe: 0.00192
    Co: 0.02450
    Ni: 0.03045
  structure:
    type: specific
    velocity:
      start: 12000.000 km/s
      stop: 35000.000 km/s
      num: 10
    density:
      type: branch85_w7
montecarlo:
  black_body_sampling:
    num: 1000000
    start: 1 angstrom
    stop: 1000000 angstrom
  convergence_criteria:
    damping_constant: 1.0
    fraction: 0.8
    hold: 3
    t_inner:
      damping_constant: 1.0
    threshold: 0.05
    type: specific
  iterations: 15
  last_no_of_packets: 20000
  no_of_packets: 10000
  no_of_virtual_packets: 3
  nthreads: 16
  seed: 23111963
plasma:
  disable_electron_scattering: false
  excitation: dilute-lte
  ionization: nebular
  line_interaction_type: downbranch
  radiative_rates_type: dilute-blackbody
spectrum:
  num: 10000
  start: 500 angstrom
  stop: 20000 angstrom
supernova:
  luminosity_requested: 8.675 log_lsun
  luminosity_wavelength_end: 6500.000 angstrom
  luminosity_wavelength_start: 3500.000 angstrom
  time_explosion: 13.0 day
tardis_config_version: v1.0

@unoebauer
Copy link
Contributor Author

@aoifeboyle, @ssim: From the comparison calculations I've performed, it looks like very high densities are triggering the problem. The branch85_w7 has lower densities in the innermost cells compared to what has been set up in densities_test.dat. When choosing a longer time since explosion (homologous expansion leads to lower densities), the differences are very weak, comparable to the MC noise level.

This could also explain, why Travis didn't catch this problem and why it doesn't show up in the tardis_example setup: we typically use times since explosions around 15d. The resulting densities are comparable to the last test run.

@aoifeboyle
Copy link
Contributor

@unoebauer I find that the differences are small when using the original yaml file, but with scatter instead of downbranch. This is a comparison for commit 45eca24 (blue) vs. the current version (green).
scatter_comparison

@unoebauer
Copy link
Contributor Author

@aoifeboyle: Very interesting! - maybe that's an indication that something is wrong with the downbranch scheme....How does it look if you use macroatom instead?

@aoifeboyle
Copy link
Contributor

The issue also occurs when macroatom is used. That narrows the issue down a lot. So something to do with high densities and downbranch/macroatom.

@unoebauer
Copy link
Contributor Author

Very interesting, @aoifeboyle. It may not be so simple though...I've tried to trigger the issue with the tardis_example setup as well. I've changed a few parameters to mimic the tardis_test_unifab_branch setup, but so far the problem didn't show up (in all of the runs macroatom has been used):

  • only changing the time to explosion to 6days: no differences
  • time since explosion 6days and lower luminosity, 8.6 log Lsun: no difference
  • time since explosion 6days, 8.6 log Lsun, wavelength cut for requested luminosity 3500 - 6000A: no difference

Update

No difference always refers to "above the noise level" when comparing the current version with commit 45eca24. The only test in which the differences might be slightly larger than the noise level is the second one...

@aoifeboyle
Copy link
Contributor

I've narrowed the issue down to be something related to the j_blues values. Hoping to find the exact issue soon.

@unoebauer
Copy link
Contributor Author

@aoifeboyle: Thanks - hopefully that's it. We will discuss this issue today during the workshop in Wuerzburg.

@aoifeboyle
Copy link
Contributor

@unoebauer I've been trying to work this out all day today and I'm not getting very far. I think it might be something to do with the montecarlo stuff rather than the plasma. I fixed the zeta problem in d864eb3, since it's the first one that fails, and I've been comparing it with 45eca24, the last one that works. The two plasmas seem identical after they are first initialised and first updated, and then after the first montecarlo runs the returned packet nus average out to a different value. For some reason the montecarlo runner is included in 45eca24 but not in d864eb3. I'll try to make some more progress tomorrow. @wkerzendorf do you have any ideas?

@unoebauer
Copy link
Contributor Author

@aoifeboyle Thanks for the detailed analysis. This seems to be consistent with some of my findings: as a test I looked at two different setups (always comparing current version with 45eca24):

  • fix initial black body temperature and radiation temperature to 10000 K and perform only one ionisation iteration and the final spectral run. Only having one plasma update seems to be enough to introduce significant differences in the final spectra
  • same as above, now with no ionisation iterations but only the spectral calculation: the spectra are almost identical - only in two line features the differences may be above the noise level

I think that these two tests point into the same direction, namely that something goes wrong when updating the plasma. In my opinion there are two possibilities:

  • something got screwed up in C part of the Monte Carlo routine and the nu averages and/or related determinations go wrong
  • something was changed in the convergence scheme, i.e. how to determine new values for T_R, T_BB and W from their old values and the radiation field state reconstructed from the Monte Carlo step

@wkerzendorf performed a git-bisect bug search and identified a commit in which only the C part was changed. I was sceptical at the beginning that we performed the correct search, but thinking about it in light of @aoifeboyle's findings and these tests, git-bisect may have been right. @wkerzendorf - can you elaborate on the git-bisect results?

After X-mas of course ;)

@unoebauer
Copy link
Contributor Author

FYI @aoifeboyle, @wkerzendorf, @ssim: I've had a closer look:

UPDATE

See my follow-up comment: discrepancies in tau_sobolevs are mostly likely not real. The virtual spectrum comparison plot and the printout of the nu_bar_estimators are still relevant, though

ORIGINAL POST

I reinvestigated the problem with the fixed plasma and boundary state, running only the spectral calculation (the yaml file is posted below). Comparing the results obtained with 45eca24 and the current branch revealed the following differences in the plasma state:

  • maximum relative difference in the electron densities: 0
  • maximum relative difference in the ion populations: ~1e-15
  • maximum relative difference in the level populations: ~1e-15
  • maximum relative difference in the tau sobolevs: ~1e+65

There are a number of line transitions (more accurately 13519), for which the relative difference in the tau Sobolevs is above 1e-10. Some of these lines are quite strong, as the following plot of the relative vs. the absolute difference in tau Sobolev shows:

00173_2_fix_tau_sob_diff_cmp

These differences could also explain the (small but visible) discrepancies in the virtual spectrum:

00173_2_fix_spectrum_cmp

It could be possible that these differences in the tau Sobolevs are the reason for the different nu_bar estimates after the Monte Carlo calculation. Having looked into the source code, I didn't spot any differences in the algorithm with which these estimates are calculated between the two Tardis versions. Moreover, the relative differences between the nubar_estimators are on the few per cent level:

In [1]: config = yaml.safe_load(open("tardis_00173_2_fix.yml"))                                                                                                

In [2]: mdl_old = tardis.run_tardis(config)
tardis.io.config_reader - INFO - Reading Atomic Data from kurucz_cd23_chianti_H_He.h5
tardis.atomic - INFO - Read Atom Data with UUID=5ca3035ca8b311e3bb684437e69d75d7 and MD5=21095dd25faa1683f4c90c911a00c3f8
tardis.io.model_reader - WARNING - v_inner_boundary requested too small for readin file. Boundary shifted to match file.
tardis.io.model_reader - WARNING - v_outer_boundary requested too large for readin file. Boundary shifted to match file.
tardis.io.config_reader - WARNING - Abundances have not been normalized to 1. - normalizing
tardis.io.config_reader - WARNING - No "species" given - ignoring other NLTE options given:
{   'classical_nebular': False, 'coronal_approximation': False}
tardis.io.config_reader - WARNING - No convergence criteria selected - just damping by 0.5 for w, t_rad and t_inner
tardis.simulation - INFO - Doing last run
tardis.model - INFO - Calculating J_blues for radiative_rates_type=dilute-blackbody
tardis.packet_source - INFO - Calculating 500000 packets for t_inner=10000.00
Running with OpenMP - 16 threadstardis.simulation - INFO - Finished in 1 iterations and took 498.82 s

In [3]: mdl_old.runner.nu_bar_estimator
Out[3]: 
array([  5.86773400e+28,   5.02563144e+28,   4.34765138e+28,
         3.70558614e+28,   2.96558255e+28,   2.02856011e+28,
         1.37465171e+28,   1.25736437e+28,   1.19699029e+28,
         1.14731265e+28,   1.10559996e+28,   1.04484887e+28,
         9.92324586e+27,   9.82713255e+27,   1.04365404e+28,
         1.14142807e+28,   1.26513448e+28,   1.42094539e+28,
         1.61340125e+28,   1.85833208e+28])
  • current
In [1]: config = yaml.safe_load(open("tardis_00173_2_fix.yml"))

In [2]: mdl_new = tardis.run_tardis(config)
tardis.io.config_reader - INFO - Reading Atomic Data from kurucz_cd23_chianti_H_He.h5
tardis.atomic - INFO - Read Atom Data with UUID=5ca3035ca8b311e3bb684437e69d75d7 and MD5=21095dd25faa1683f4c90c911a00c3f8
tardis.io.model_reader - WARNING - v_inner_boundary requested too small for readin file. Boundary shifted to match file.
tardis.io.model_reader - WARNING - v_outer_boundary requested too large for readin file. Boundary shifted to match file.
tardis.io.config_reader - WARNING - Abundances have not been normalized to 1. - normalizing
tardis.io.config_reader - WARNING - No "species" given - ignoring other NLTE options given:
{   'classical_nebular': False, 'coronal_approximation': False}
tardis.io.config_reader - WARNING - No convergence criteria selected - just damping by 0.5 for w, t_rad and t_inner
tardis.plasma.properties.atomic - WARNING - Zeta_data missing - replaced with 1s. Missing ions: [(9, 10), (11, 12), (12, 13), (13, 14), (14, 15), (15, 16), (16, 17), (17, 18), (18, 19), (19, 20), (20, 21), (21, 22), (22, 23), (23, 24), (24, 25), (25, 26), (26, 27), (27, 28), (28, 29), (29, 1), (29, 2), (29, 3), (29, 4), (29, 5), (29, 6), (29, 7), (29, 8), (29, 9), (29, 10), (29, 11), (29, 12), (29, 13), (29, 14), (29, 15), (29, 16), (29, 17), (29, 18), (29, 19), (29, 20), (29, 21), (29, 22), (29, 23), (29, 24), (29, 25), (29, 26), (29, 27), (29, 28), (29, 29), (29, 30), (30, 1), (30, 2), (30, 3), (30, 4), (30, 5), (30, 6), (30, 7), (30, 8), (30, 9), (30, 10), (30, 11), (30, 12), (30, 13), (30, 14), (30, 15), (30, 16), (30, 17), (30, 18), (30, 19), (30, 20), (30, 21), (30, 22), (30, 23), (30, 24), (30, 25), (30, 26), (30, 27), (30, 28), (30, 29), (30, 30), (30, 31)]
tardis.model - INFO - Calculating J_blues for radiative_rates_type=dilute-blackbody
tardis.simulation.base - INFO - Doing last run
Running with OpenMP - 1 threads
tardis.simulation.base - INFO - Finished in 0 iterations and took 149.83 s

In [3]: mdl_new.runner.nu_bar_estimator
Out[3]: 
array([  6.01178210e+28,   5.13461157e+28,   4.43216973e+28,
         3.77940227e+28,   3.04304078e+28,   2.10994277e+28,
         1.43991981e+28,   1.31111096e+28,   1.22309477e+28,
         1.16231763e+28,   1.11305261e+28,   1.05664707e+28,
         1.00899515e+28,   1.00840642e+28,   1.07745764e+28,
         1.17769597e+28,   1.30562338e+28,   1.46607607e+28,
         1.66597070e+28,   1.91696842e+28])

What we see in the original setup (posted at the beginning of this issue) is a runaway process during many (20) ionization iterations.

Config file: tardis_00173_2_fix.yml

atom_data: kurucz_cd23_chianti_H_He.h5
model:
  abundances:
    filename: abundances_00173_2.dat
    filetype: simple_ascii
    type: file
  structure:
    filename: densities_00173_2.dat
    filetype: simple_ascii
    type: file
    v_inner_boundary: 12000.000 km/s
    v_outer_boundary: 35000.000 km/s
montecarlo:
  black_body_sampling:
    num: 1000000
    start: 1 angstrom
    stop: 1000000 angstrom
  convergence_criteria:
    damping_constant: 1.0
    fraction: 0.8
    hold: 3
    t_inner:
      damping_constant: 1.0
    threshold: 0.05
    type: specific
  iterations: 1
  last_no_of_packets: 500000
  no_of_packets: 50000
  no_of_virtual_packets: 3
  nthreads: 16
  seed: 23111963
plasma:
  initial_t_rad: 10000 K
  initial_t_inner: 10000 K
  disable_electron_scattering: false
  excitation: dilute-lte
  ionization: nebular
  line_interaction_type: downbranch
  radiative_rates_type: dilute-blackbody
spectrum:
  num: 10000
  start: 500 angstrom
  stop: 20000 angstrom
supernova:
  luminosity_requested: 8.675 log_lsun
  luminosity_wavelength_end: 6500.000 angstrom
  luminosity_wavelength_start: 3500.000 angstrom
  time_explosion: 6.800 day
tardis_config_version: v1.0

@yeganer
Copy link
Contributor

yeganer commented Jan 15, 2016

I did another bisect run and there seem to be several bugs indeed.
when using

plasma:
  ionization: lte

bisect reports e9fb778 as first bad commit.
But using

plasma:
  ionization: nebular

gives commits e2eac6d - ec8fcda as possible first bad commits.

Maybe one of the issues was fixed by @aoifeboyle

I didn't look into the case with ionization: lte further but for ionization: nebular i did several binary dumps of the data passed to the storage model and electron_density, inverse_electron_density, line_lists_tau_sobolevs and transition_probabilities have relative errors > 1e-5 with max errors >1

@aoifeboyle
Copy link
Contributor

@wkerzendorf I accidentally deleted my comment from last night, could you please paste it into this conversation somewhere? I think you should have it as an email?

The bug that I found in the plasma/restructure commit is something that I had already corrected at a later time, but it was not corrected at the time of the merge of the plasma/restructure. So I think that explains what you found @yeganer. My aim for now is to get consistent results for our test between the last working version (45eca24) and the plasma restructure commit (41d8fcd). There are too many differences between 45eca24 and the current version to easily debug it. But once I get this right, it should either solve all the problems or at least considerably reduce the number of problems.

@yeganer
Copy link
Contributor

yeganer commented Jan 15, 2016

@aoifeboyle in my opinion it is better to compare 30d1cb9 with ec8fcda because that's the smalles difference in code where the bug appears.

I added your comment at the bottom of my comment, if that's okay

@unoebauer
Copy link
Contributor Author

I am re-posting @aoifeboyle's originial comment:

The bug that I found in the plasma/restructure commit is something that I had already corrected at a later time, but it was not corrected at the time of the merge of the plasma/restructure. So I think that explains what you found. My aim for now is to get consistent results for our test between the last working version (45eca24) and the plasma restructure commit (41d8fcd). There are too many differences between 45eca24 and the current version to easily debug it. But once I get this right, it should either solve all the problems or at least considerably reduce the number of problems.

@aoifeboyle
Copy link
Contributor

Also @yeganer, what I'm comparing is basically the same thing that you're suggesting. 30d1cb9 is effectively 45eca24 because the model was still calling the old plasma structure. And ec8fcda is when the nebular ionisation was added to the new plasma, but this was on the old plasma/restructure branch rather than on the main Tardis branch, so it's less confusing for me to work with 41d8fcd.

I think I'm making good progress anyway. The only discrepancy I can find in the plasma now is j_blues being an array in one version and a dataframe in the other, so I'm going to correct that.

@aoifeboyle
Copy link
Contributor

@unoebauer @yeganer
PR #462 is an updated version of 41d8fcd (the point of new plasma structure merging). It now produces the exact same plasma, transition probabilities, taus, etc. as 45eca24. However, the montecarlo code still returns the wrong nubar_estimators, even though the montecarlo code is the same and the inputs from the plasma (appear) the same. I think this is possibly a view/copy issue. @yeganer you mentioned before that you know how to examine the binary data being passed from the python to C. Could you checkout and perform the same analysis to compare 2e7a85f (the fix I uploaded) and 45eca24?

@unoebauer
Copy link
Contributor Author

@aoifeboyle: I just looked at your PR #462. Unfortunately, it crashes when running the plasma_fail_test.yml problem:

In [2]: config = yaml.safe_load(open("plasma_fail_test.yml"))

In [3]: mdl = tardis.run_tardis(config)
tardis.io.config_reader - INFO - Reading Atomic Data from kurucz_cd23_chianti_H_He.h5
tardis.atomic - INFO - Read Atom Data with UUID=5ca3035ca8b311e3bb684437e69d75d7 and MD5=21095dd25faa1683f4c90c911a00c3f8
tardis.io.model_reader - WARNING - v_inner_boundary requested too small for readin file. Boundary shifted to match file.
tardis.io.model_reader - WARNING - v_outer_boundary requested too large for readin file. Boundary shifted to match file.
tardis.io.config_reader - WARNING - Abundances have not been normalized to 1. - normalizing
tardis.io.config_reader - WARNING - No "species" given - ignoring other NLTE options given:
{   'classical_nebular': False, 'coronal_approximation': False}
tardis.io.config_reader - WARNING - No convergence criteria selected - just damping by 0.5 for w, t_rad and t_inner
tardis.plasma.properties.atomic - WARNING - Zeta_data missing - replaced with 1s
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-3-6282305a6c09> in <module>()
----> 1 mdl = tardis.run_tardis(config)

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/base.pyc in run_tardis(config, atom_data)
     38     tardis_config = config_reader.Configuration.from_config_dict(
     39         config_dict, atom_data=atom_data)
---> 40     radial1d_mdl = model.Radial1DModel(tardis_config)
     41 
     42     simulation.run_radial1d(radial1d_mdl)

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/model.pyc in __init__(self, tardis_config)
    130                                                          excitation_mode=tardis_config.plasma.excitation,
    131                                                          line_interaction_type=tardis_config.plasma.line_interaction_type,
--> 132                                                          link_t_rad_t_electron=0.9)
    133 
    134         self.spectrum = TARDISSpectrum(tardis_config.spectrum.frequency, tardis_config.supernova.distance)

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/plasma/standard_plasmas.pyc in __init__(self, number_densities, atomic_data, time_explosion, t_rad, delta_treatment, nlte_config, ionization_mode, excitation_mode, line_interaction_type, link_t_rad_t_electron)
    103             delta_input=delta_treatment, nlte_species=nlte_config.species,
    104             previous_electron_densities=initial_electron_densities,
--> 105             previous_beta_sobolevs=initial_beta_sobolevs)

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/plasma/base.pyc in __init__(self, plasma_properties, **kwargs)
     20         self._build_graph()
     21 #        self.write_to_tex('Plasma_Graph', 'Plasma_Formulae')
---> 22         self.update(**kwargs)
     23 
     24     def __getattr__(self, item):

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/plasma/base.pyc in update(self, **kwargs)
    141 
    142         for module_name in self._resolve_update_list(kwargs.keys()):
--> 143             self.plasma_properties_dict[module_name].update()
    144 
    145     def _update_module_type_str(self):

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/plasma/properties/base.pyc in update(self)
     84         if len(self.outputs) == 1:
     85             setattr(self, self.outputs[0], self.calculate(
---> 86                 *self._get_input_values()))
     87         else:
     88             new_values = self.calculate(*self._get_input_values())

/afs/mpa/home/ulnoe/local/lib/python2.7/site-packages/tardis_sn-1.0.1-py2.7-linux-x86_64.egg/tardis/plasma/properties/radiative_properties.pyc in calculate(lines, j_blues_array, beta_rad)
    159     @staticmethod
    160     def calculate(lines, j_blues_array, beta_rad):
--> 161         return pd.DataFrame(j_blues_array, index=lines.index)
    162 
    163 class LTEJBlues(ProcessingPlasmaProperty):

/opt/python-2.7.9/lib/python2.7/site-packages/pandas/core/frame.pyc in __init__(self, data, index, columns, dtype, copy)
    237             else:
    238                 mgr = self._init_ndarray(data, index, columns, dtype=dtype,
--> 239                                          copy=copy)
    240         elif isinstance(data, (list, types.GeneratorType)):
    241             if isinstance(data, types.GeneratorType):

/opt/python-2.7.9/lib/python2.7/site-packages/pandas/core/frame.pyc in _init_ndarray(self, values, index, columns, dtype, copy)
    405             values = _possibly_infer_to_datetimelike(values)
    406 
--> 407         return create_block_manager_from_blocks([values], [columns, index])
    408 
    409     @property

/opt/python-2.7.9/lib/python2.7/site-packages/pandas/core/internals.pyc in create_block_manager_from_blocks(blocks, axes)
   3526         blocks = [getattr(b, 'values', b) for b in blocks]
   3527         tot_items = sum(b.shape[0] for b in blocks)
-> 3528         construction_error(tot_items, blocks[0].shape[1:], axes, e)
   3529 
   3530 

/opt/python-2.7.9/lib/python2.7/site-packages/pandas/core/internals.pyc in construction_error(tot_items, block_shape, axes, e)
   3507         raise e
   3508     raise ValueError("Shape of passed values is {0}, indices imply {1}".format(
-> 3509         passed,implied))
   3510 
   3511 

ValueError: Shape of passed values is (0, 0), indices imply (0, 270394)

In [4]: 

@aoifeboyle
Copy link
Contributor

@unoebauer Strange - it worked fine for me. Let me look and see where I went wrong.

@aoifeboyle
Copy link
Contributor

@unoebauer I tried checking it out directly from the PR rather than using my own branch and it still works for me. Is this the test you're talking about?

atom_data: kurucz_cd23_chianti_H_He.h5
model:
abundances:
filename: abundances_test.dat
filetype: simple_ascii
type: file
structure:
filename: densities_test.dat
filetype: simple_ascii
type: file
v_inner_boundary: 12000.000 km/s
v_outer_boundary: 35000.000 km/s
montecarlo:
black_body_sampling:
num: 1000000
start: 1 angstrom
stop: 1000000 angstrom
iterations: 1
last_no_of_packets: 1000
no_of_packets: 1000
no_of_virtual_packets: 3
nthreads: 16
seed: 23111963
plasma:
initial_t_rad: 10000 K
initial_t_inner: 10000 K
disable_electron_scattering: false
excitation: dilute-lte
ionization: nebular
line_interaction_type: downbranch
radiative_rates_type: dilute-blackbody
spectrum:
num: 10000
start: 500 angstrom
stop: 20000 angstrom
supernova:
luminosity_requested: 8.675 log_lsun
luminosity_wavelength_end: 6500.000 angstrom
luminosity_wavelength_start: 3500.000 angstrom
time_explosion: 6.800 day
tardis_config_version: v1.0

@aoifeboyle
Copy link
Contributor

In case anyone's following this, the problem is apparently just our different pandas versions. I'm going to tweak it a bit so it works with both versions now.

@aoifeboyle
Copy link
Contributor

@unoebauer I've fixed it for pandas-0.16.0. Could you please check that it's working for you now?

@unoebauer
Copy link
Contributor Author

@aoifeboyle: Looks good now - thanks. The plasma_fail_test.yml setup runs without problems and the Travis tests pass.

@unoebauer
Copy link
Contributor Author

The discussion is getting out of hand. Too many loose and dead ends have accumulated here (and I am as much to blame for this as anybody) and are muddying the picture. I think it's time to take a step back, review the current status and take a fresh, re-focused start. For this purpose, I am closing this issue and will open a new one in the hope of simplifying the problem solving process. Please don't post anything here any more but use the new issue, #464.

unoebauer added a commit to unoebauer/tardis that referenced this issue Jan 16, 2016
yeganer pushed a commit to yeganer/tardis that referenced this issue Jan 19, 2016
yeganer pushed a commit to yeganer/tardis that referenced this issue Jan 20, 2016
yeganer pushed a commit to yeganer/tardis that referenced this issue Jan 20, 2016
unoebauer added a commit to unoebauer/tardis that referenced this issue Jan 22, 2016
* added config files for the original  setup triggering tardis-sn#455
* renamed it from tardis_00173_2 to abn_tom_test (for abundance
  tomography test)
yeganer pushed a commit to yeganer/tardis that referenced this issue Feb 1, 2016
Deleted atom_data._lines and atom_data._levels
This change ensures that the BasePlasma always gets prepared data.
A proper fix will take more time and will be done later.
yeganer pushed a commit to yeganer/tardis that referenced this issue Feb 2, 2016
The old files contained wrong values (see issue tardis-sn#455)
With this change the following tests will run again:

* test_packet_output
* test_plasma_estimates
* TestPlasma.test_lte_plasma

One should consider updating tardis/simulation/tests/test_simulation.py
since it creates a new model which essentially tests the whole codebase.
These tests will always fail if something in the codebase changes ( see comment
in the file ).
yeganer pushed a commit to yeganer/tardis that referenced this issue Feb 3, 2016
Changed behavior to raise a custom exception instead of a generic one.
This is done to better reflect the error.

Related Issues: tardis-sn#455 tardis-sn#464
Pull Request tardis-sn#471
unoebauer added a commit that referenced this issue Feb 3, 2016
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

7 participants