wtbarnes / fiasco

Python interface to the CHIANTI atomic database
http://fiasco.readthedocs.io/
BSD 3-Clause "New" or "Revised" License
21 stars 17 forks source link

Differing behavior on two different systems #304

Open STBadman opened 2 weeks ago

STBadman commented 2 weeks ago

With fiasco.version = 0.2.3, numpy.version = 1.24.3, astropy.version=6.0.1, python 3.11 i get the following output for creating an oxygen element for 1e5K which seems reasonable:

from fiasco import Element
import astropy.units as u
Element('oxygen',1e5*u.K).equilibrium_ionization
<Quantity [[3.44576031e-05, 4.01115956e-02, 6.47236538e-01,
            3.09443217e-01, 3.17408912e-03, 1.03448308e-07,
            1.02775598e-12, 2.62084865e-16, 1.96052787e-16]]>

but my student running exactly the same thing with the same set of versions above gets a fully neutral result for the same input. Is there anything I'm not checking to ensure things are set up identically?


from fiasco import Element
import astropy.units as u
Element('oxygen',1e5*u.K).equilibrium_ionization
<Quantity [[1, 0, 0, 0, 0, 0, 0, 0, 0]]>
wtbarnes commented 2 weeks ago

That's really odd. I just built a fresh environment using the versions you list above (on an Apple M1 Macbook Pro) and I get an answer that is the same (up to numerical precision) with your first result. This is also the answer I would expect given the formation temperatures of the O ions.

A few things to try:

  1. What version of the CHIANTI database are they using? Is it possible this file is missing? Have a look in $HOME/.fiasco.
  2. Send your student the exact database file you're using. If you didn't specify otherwise, this should live in $HOME/.fiasco/chianti_dbase.h5. They can either copy it to that same location on their machine or point to a custom location using the hdf5_dbase_root= keyword argument to Element. If they then get the same answer as you do, then something probably went wrong when they built the HDF5 version of the database.
  3. Have a look at the ionization and recombination rates of each O ion,
    oxygen = fiasco.Element('O',1e5*u.K)
    for ion in oxygen:
    print(ion.recombination_rate)
    print(ion.ionization_rate)

    If both of you still get the same answer, then the actual solution to the ionization equilibria is behaving differently which is more worrying.

STBadman commented 2 weeks ago

We'll try using an identical database file between our systems, will let you know.

Another puzzling thing is we get identical results to each other (to floating point accuracy) if we change oxygen to carbon or iron.

wtbarnes commented 2 weeks ago

That is even more puzzling!

If you loop over all elements (fiasco.list_elements()), is it only oxygen that is the problem?