Skip to content

Bug in contribution function? #447

@jwreep

Description

@jwreep

I was running the tutorial on calculating contribution functions, and hit a snag. Not sure I understand the root of the issue, but this is with fiasco version 0.8.1 and dbase version 8.0.7. It looks like the issue arises during the level population calculation, so possibly a change relating to the updates there with newer versions of the database is causing a problem?

I ran the following lines (copied verbatim from the how-to):

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np

from astropy.visualization import quantity_support

from fiasco import Ion

quantity_support()

Te = np.geomspace(0.1, 10, 26) * u.MK
ne = 1e8 * u.cm**-3

ion = Ion('O 5+', Te)
print(ion)

contribution_func = ion.contribution_function(ne)

which gives the following traceback:

WARNING: No autoionization or level-resolved radiative recombination data available for O 6. Using single-ion model for level populations calculation. To silence this warning, set use_two_ion_model=False. [fiasco.ions]
WARNING: No proton data available for O 6. Not including proton excitation and de-excitation in level populations calculation. [fiasco.ions]

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[4], line 1
----> 1 contribution_func = ion.contribution_function(ne)

File ~/anaconda3/lib/python3.12/site-packages/fiasco/util/decorators.py:32, in needs_dataset.<locals>.decorator.<locals>.func_wrapper(*args, **kwargs)
     27     except KeyError:
     28         raise MissingDatasetException(
     29             f"{n} dataset missing for {getattr(obj, 'ion_name', obj)}."
     30         )
---> 32 return func(*args, **kwargs)

File ~/anaconda3/lib/python3.12/site-packages/astropy/units/decorators.py:316, in QuantityInput.__call__.<locals>.wrapper(*func_args, **func_kwargs)
    314 # Call the original function with any equivalencies in force.
    315 with equiv_context:
--> 316     return_ = wrapped_function(*func_args, **func_kwargs)
    318 # Return
    319 ra = wrapped_signature.return_annotation

File ~/anaconda3/lib/python3.12/site-packages/fiasco/ions.py:1135, in Ion.contribution_function(self, density, **kwargs)
   1083 r"""
   1084 Contribution function :math:`G(n_e,T)` for all transitions.
   1085 
   (...)   1132 level_populations
   1133 """
   1134 couple_density_to_temperature = kwargs.get('couple_density_to_temperature', False)
-> 1135 populations = self.level_populations(density, **kwargs)
   1136 if couple_density_to_temperature:
   1137     term = self.ionization_fraction / density

File ~/anaconda3/lib/python3.12/site-packages/astropy/units/decorators.py:316, in QuantityInput.__call__.<locals>.wrapper(*func_args, **func_kwargs)
    314 # Call the original function with any equivalencies in force.
    315 with equiv_context:
--> 316     return_ = wrapped_function(*func_args, **func_kwargs)
    318 # Return
    319 ra = wrapped_signature.return_annotation

File ~/anaconda3/lib/python3.12/site-packages/fiasco/ions.py:639, in Ion.level_populations(self, density, include_protons, include_level_resolved_rate_correction, couple_density_to_temperature, use_two_ion_model)
    637 b = np.zeros(c_matrix.shape[2:])
    638 b[-1] = 1.0
--> 639 pop = np.linalg.solve(c_matrix.value, b)
    640 pop[pop<0] = 0.0
    641 np.divide(pop, pop.sum(axis=1)[:, np.newaxis], out=pop)

File ~/anaconda3/lib/python3.12/site-packages/numpy/linalg/linalg.py:409, in solve(a, b)
    407 signature = 'DD->D' if isComplexType(t) else 'dd->d'
    408 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 409 r = gufunc(a, b, signature=signature, extobj=extobj)
    411 return wrap(r.astype(result_t, copy=False))

ValueError: solve: Input operand 1 does not have enough dimensions (has 1, gufunc core with signature (m,m),(m,n)->(m,n) requires 2)

Metadata

Metadata

Assignees

No one assigned

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions