ajwheeler / ajwheeler/Korg.jl

Faster `voigt`, `exp`

Ouverte
#80 6 commentaires 0 réactions 0 personnes assignées Voir sur GitHub
Langage dominant
Julia
Étoiles
67
Forks
15
Merge moyen
2 j 6 h
PR mergées (30 j)
6

Description

As promised, I'm posting my faster exponential function. I'm curious to hear if using this inside of `voigt` buys you any speed.

Tonight, I also put together a branchless version of `voigt`. Unfortunately, all of the extra work that is being done by the branchless voigt makes it slower than the unvectorized version. I suspect that if someone spent a lot of time refactoring, they might be able to get it to be slightly faster. But I don't think that's worth the time. As you had suggested a while back, I think a better use of time would be to replace the current approximation with a better approach.

I did some googling and the consensus seems to be that we should leverage the relationship of the Voigt function with the real function. According to [eqn 9 of Thompson (93)](https://aip.scitation.org/doi/abs/10.1063/1.4823236), `V(u,a) = Re[exp(z²) erfc(z)]`, where `z = a + i*u`. I suspect that nobody will have time to pursue this for a while, so I'll briefly document what I read. I've seen 2 ways to do this (there are probably more):

1. I think that [Turbospectrum's implementation](https://github.com/bertrandplez/Turbospectrum2019/blob/master/source-v19.1/newvoigt.f) is based on the same underlying approach as [Kuntz (1997)](https://www.sciencedirect.com/science/article/pii/S0022407396001628). (Interestingly, I think that turbospectrum and the article might be suggesting that for `V(a,v)`, `z = a - iv`...). Under this approach, I think that they directly approximate the whole expression `Re[exp(z²) erfc(z)]`, but it introduces at least some level of branching. I think this article talks a little about [variants of this approach](https://elib.dlr.de/120657/1/FSchreier_Voigt_1806.11560.pdf).

2. At a glance, I suspect this might potentially be slightly faster to adopt a strategy like the one used in [this article](https://opensky.ucar.edu/islandora/object/conference%3A828/datastream/PDF/view) (if we target a particular accuracy). They seemed to use an approximation for `exp(z²) erfc(z) ≈ (Σⱼᵖ aⱼ zʲ)/(z⁽ᵖ⁺¹⁾ + Σⱼᵖ bⱼ zʲ)` with coefficients from [here](https://ui.adsabs.harvard.edu/abs/1978JQSRT..19..509H/abstract), and directly computed `exp(z²)`. The benefit of this approach is that there is basically no-branching. However, plugging `Re(exp(-1*(x + i*y)^2) * (a + i b))` into wolfram alpha suggests that you also need to evaluate sin and cos (which I think boils down to 2 more polynomial evaluations)...

If either of these approaches are adopted, I believe some care is required to avoid using Julia's complex datatype (since I think ForwardDiff.jl is incompatible with complex numbers)

Guide de contribution

Aucun guide de contribution indexé pour ce dépôt

Évaluation

Cette issue n'a pas encore été évaluée.

Recevez les nouvelles issues par e-mail

Un résumé court des issues GitHub adaptées aux débutants.