Korg fails for some atmospheres

[ajwheeler]

edit: scroll down or click here for the current status

---

Details to come. The cool stars failures are OK for now, since we are targeting FGK in paper 1, but obviously we should figure out what's going on. I'm more concerned about the 3 in the middle.

[ajwheeler]

The few points in the middle were just parsing problems. #130 fixes them.

For the cool stars, the chemical equilibrium is either not converging, or returning negative number densities which causes things to crash later. Two possibilities:

• the solver fails, because it's starting point is the analytically-calculated no-molecules solution, and cool stars have lots of molecules
• somehow things don't work because MARCS and Korg don't use the same set of species

[segasai]

Hi,
I don't quite know if it's the same issue with cool stars but I also see failures for Teff=4000K MARCS
I.e. p4000_g+4.5_m0.0_t02_an_z+1.00_a-0.40_c+0.00_n+0.00_o-0.40_r+0.00_s+0.00.mod
with this abundance (1.0, {'O': 0.6, 'Ne': 0.6, 'Mg': 0.6, 'Si': 0.6, 'S': 0.6, 'Ar': 0.6, 'Ca': 0.6, 'Ti': 0.6})

ERROR: DomainError with -3.50386098169805e6:
log will only return a complex result if called with a complex argument. Try log(Complex(x)).
Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)
    @ Base.Math ./math.jl:33
  [2] _log(x::Float64, base::Val{:ℯ}, func::Symbol)
    @ Base.Math ./special/log.jl:301
  [3] log
    @ ./special/log.jl:267 [inlined]
  [4] metal_bf_absorption!(α::Vector{Float64}, νs::Vector{Float64}, T::Float64, number_densities::Dict{Korg.Species, Float64})
    @ Korg.ContinuumAbsorption ~/.julia/packages/Korg/Aas4k/src/ContinuumAbsorption/absorption_metals_bf.jl:39
  [5] total_continuum_absorption(νs::Vector{Float64}, T::Float64, nₑ::Float64, number_densities::Dict{Korg.Species, Float64}, partition_funcs::Dict{Korg.Species, Any}; error_oobounds::Bool)
    @ Korg.ContinuumAbsorption ~/.julia/packages/Korg/Aas4k/src/ContinuumAbsorption/ContinuumAbsorption.jl:72
  [6] total_continuum_absorption(νs::Vector{Float64}, T::Float64, nₑ::Float64, number_densities::Dict{Korg.Species, Float64}, partition_funcs::Dict{Korg.Species, Any})
    @ Korg.ContinuumAbsorption ~/.julia/packages/Korg/Aas4k/src/ContinuumAbsorption/ContinuumAbsorption.jl:46
  [7] (::Korg.var"#97#106"{Float64, Bool, Bool, Dict{Korg.Species, Any}, StepRangeLen{Float64, Float64, Float64, Int64}, Vector{Float64}, Matrix{Float64}, NamedTuple{(:atoms, :molecules, :equations, :absolute_abundances, :ionization_energies, :partition_fns, :equilibrium_constants), Tuple{UnitRange{UInt8}, Base.KeySet{Korg.Species, Dict{Korg.Species, Korg.CubicSplines.CubicSpline{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}}, Korg.var"#system#63"{Vector{Float64}, Dict{UInt8, Vector{Float64}}, Dict{Korg.Species, Any}, Dict{Korg.Species, Korg.CubicSplines.CubicSpline{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, Base.KeySet{Korg.Species, Dict{Korg.Species, Korg.CubicSplines.CubicSpline{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}}, UnitRange{UInt8}}, Vector{Float64}, Dict{UInt8, Vector{Float64}}, Dict{Korg.Species, Any}, Dict{Korg.Species, Korg.CubicSplines.CubicSpline{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}}}, Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}})(::Tuple{Int64, Korg.PlanarAtmosphereLayer{Float64}})
    @ Korg ~/.julia/packages/Korg/Aas4k/src/synthesize.jl:133
  [8] iterate
    @ ./generator.jl:47 [inlined]
  [9] collect(itr::Base.Generator{Base.Iterators.Enumerate{Vector{Korg.PlanarAtmosphereLayer{Float64}}}, Korg.var"#97#106"{Float64, Bool, Bool, Dict{Korg.Species, Any}, StepRangeLen{Float64, Float64, Float64, Int64}, Vector{Float64}, Matrix{Float64}, NamedTuple{(:atoms, :molecules, :equations, :absolute_abundances, :ionization_energies, :partition_fns, :equilibrium_constants), Tuple{UnitRange{UInt8}, Base.KeySet{Korg.Species, Dict{Korg.Species, Korg.CubicSplines.CubicSpline{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}}, Korg.var"#system#63"{Vector{Float64}, Dict{UInt8, Vector{Float64}}, Dict{Korg.Species, Any}, Dict{Korg.Species, Korg.CubicSplines.CubicSpline{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, Base.KeySet{Korg.Species, Dict{Korg.Species, Korg.CubicSplines.CubicSpline{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}}, UnitRange{UInt8}}, Vector{Float64}, Dict{UInt8, Vector{Float64}}, Dict{Korg.Species, Any}, Dict{Korg.Species, Korg.CubicSplines.CubicSpline{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}}}, Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}})
    @ Base ./array.jl:787
 [10] map
    @ ./abstractarray.jl:2961 [inlined]
 [11] synthesize(atm::Korg.PlanarAtmosphere{Float64}, linelist::Vector{Korg.Line{Float64}}, A_X::Vector{Float64}, λs::StepRangeLen{Float64, Int64, Float64, Int64}; vmic::Float64, line_buffer::Float64, cntm_step::Float64, hydrogen_lines::Bool, n_mu_points::Int64, line_cutoff_threshold::Float64, bezier_radiative_transfer::Bool, ionization_energies::Dict{UInt8, Vector{Float64}}, partition_funcs::Dict{Korg.Species, Any}, equilibrium_constants::Dict{Korg.Species, Korg.CubicSplines.CubicSpline{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}})
    @ Korg ~/.julia/packages/Korg/Aas4k/src/synthesize.jl:129
 [12] synthesize
    @ ~/.julia/packages/Korg/Aas4k/src/synthesize.jl:90 [inlined]
 [13] synthesize(atm::Korg.PlanarAtmosphere{Float64}, linelist::Vector{Korg.Line{Float64}}, A_X::Vector{Float64}, λ_start::Int64, λ_stop::Int64, λ_step::Float64; air_wavelengths::Bool, wavelength_conversion_warn_threshold::Float64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Korg ~/.julia/packages/Korg/Aas4k/src/synthesize.jl:88
 [14] synthesize (repeats 2 times)
    @ ~/.julia/packages/Korg/Aas4k/src/synthesize.jl:70 [inlined]
 [15] top-level scope
    @ REPL[10]:1

(using the v0.10.1 )

[ajwheeler]

Interesting, thank you for letting me know. This is indeed crashing in the same place, and almost certainly the same fundamental issue. I'm going to prioritize fixing this, or at least failing more gracefully.

[segasai]

Thanks.
It looks like putting an assert for positive densities at synthesize.jl or in molecular_equilibrium() may be a good idea to catch this early.

[ajwheeler]

#151 makes molecular equilibrium work slightly better. It works on your example, @segasai, and it emits an error which makes more sense when the molecular equilibrium fails.

Slightly improved:

[ajwheeler]

To elaborate on what I think is the fundamental cause here, it was definitely the case that the initial guess was throwing the solver to negative values for cool stars, because it was a no-molecules solution. #151 fixes that, I think, and the fundemental problem now is that Korg can't match the MARCS electron number density because it is using a different set of species (and different partition functions, molecular dissociation energies, solver, etc.). The solution going forward is going to be to self-consistently solve for the electron number density (probably using the value from the model atmosphere as a starting point).

This is inherently inconsistent, but there's no way around it except calculating our own model atmospheres. The current approach is also inconsistent, as is what everyone else does (except maybe turbospectrum, which might share moleq code with MARCS?)

[ajwheeler]

With the addition of 28 new molecules in #155, Korg's chemical equilibrium is much more likely to converge for cool stars.

[segasai]

Great!
Is it useful to take a look at the comparison of PHOENIX and korg there ?

[ajwheeler]

It seems very likely that achieving convergence everywhere is primarily a matter of using a list of species as similar to what MARCs uses as possible (or using different atmospheres). Of course, comparisons to PHOENIX might make it easier to identify missing features, which could point towards missing molecules. In general, comparisons with anything are interesting, especially for cool stars where models struggle more.

[ajwheeler]

#179 solves this for all atmospheres with solar metallicity (though I will open an issue soon on the places where the calculated and model-atmosphere electron densities dissagree). For lower-metallicity stars, things still sometimes don't work at the top of the red giant branch. I think this is fundamentally the same thing going on: if the set of molecules considered isn't quite the same, the combination of temperature and number density can't be matched.

[ajwheeler]

#221 eliminates any chemical equilibrium failures at metallicity = -1, so when that is merged this will be improved further.

metallicty = -2

metallicity = -2.5

The new failures for cool dwarfs at metallicity=-2.5 are interesting. Perhaps related to the discrepancies in #182?

[ajwheeler]

closed by #260