Skip to content

Use Arb library for bessel function evaluation

Eskil Vik requested to merge bugfix/move-bessel-calculations-to-arb-library into master

Overview

The modified Bessel functions have until now been manually implemented in IW2D, based on Taylor series and asymptotic series approximations from https://dlmf.nist.gov. This approach led to large differences with well-regarded implementations such as that by Amos (which is used by e.g. SciPy and Matlab), especially in the intermediate region between the area of validity for the Taylor series and the area of validity for the asymptotic series. The worst cases saw a result of Bessel function evaluations many times larger than it should be.

This MR aims to fix this issue, and other potential future issues, by using the third party library Arb (https://arblib.org) to calculate the Bessel functions instead of a custom local implementation.

Fixes issues #3 (closed) and #4 (closed), and the large error present in the discussion in !13 (merged).

Implementation details

The Bessel functions are computed using the Arb library's acb_hypgeom_bessel_i and acb_hypgeom_bessel_k functions. See docs for implementation details.

The Arb library uses so called "ball arithmetic" to keep track of possible floating point errors and give a range of confidence for the result. The Arb Bessel functions are run iteratively with higher and higher precision (doubled in each iteration) until the range of confidence becomes more precise than the original machine precision. If this does not happen during the first 5 iterations, the process is stopped, the current best result is returned, and a warning is printed.

Since Arb uses its own multiprecision types, conversion functions to convert from amp::campf (used in the rest of IW2D) to acb_t and back, were implemented (acb_from_campf and campf_from_acb, respectively).

Other notes

Removed support for the exponentially normalized Bessel functions from bessel, as they were not currently used by the IW2D code, and they were no longer a "free" byproduct of computing the non-normalized Bessel functions.

The Makefiles, README and pytests were updated to account for the changes.

The reference output files in examples/output_files/ have been modified accordingly. Also, the wake in front of the source, for the round chamber wake case, being subject to numerical instabilities (i.e. machine-dependent values), most of it was removed from the test.


Documentation of Performance to follow in the comments.

Edited by Nicolas Mounet

Merge request reports