Adding computation of complete elliptic integrals into amrex::Math so… #4151
Conversation
… that it can be executed on accelerator devices.
Src/Base/AMReX_Math.H
Outdated
| { | ||
| // Computing K based on DLMF | ||
| // https://dlmf.nist.gov/19.8 | ||
| T tol = 1e-12; |
There was a problem hiding this comment.
I will add a check on it to adjust the tolerance for a float to 1e-6.
There was a problem hiding this comment.
This should be fixed now.
| g0 = g; | ||
| } | ||
|
|
||
| return 0.5*pi<T>()/a; |
There was a problem hiding this comment.
If m is 0, a is uninitialized.
There was a problem hiding this comment.
For comp_ellint_1(0.5), I get 1.685750355 with std::comp_ellint_1, but 1.854074677 with amrex::Math::comp_ellint_1. The error seems to be too big.
There was a problem hiding this comment.
This is a definition issue. In most cases the integrals are defined with m as the argument. In this case with gcc it is k = sqrt(m). I will update it to be consistent and I will additionally add some test figures doing comparisons.
There was a problem hiding this comment.
@WeiqunZhang in general though do you think this is the right spot for this calculation?
There was a problem hiding this comment.
It seems that we could use T tol = std::numeric_limits<T>::epsilon();. If we want to be conservative, we could make it 10x bigger than epsilon. In my testing, it only increases the number of iteration slightly.
There was a problem hiding this comment.
I think the latest push should resolve these concerns. I have also added a comparison of the AGM algorithm to what is done in Python's ellipm and ellipk routines.
There was a problem hiding this comment.
The above figure is the difference between the two cases for the Z-components of the analytical solution to a current loop.
…to mention that it is an argument of k.
… that it can be executed on accelerator devices.
Summary
Currently the parser recognizes comp_ellint_1 and comp_ellint_2, but only executes on CPU when compiled with gcc. An implementation of the calculations have been added to amrex::Math to compute these quantities using the method described here:
https://dlmf.nist.gov/19.8
It is a quadratically converging iterative method that is compatible with device execution.
Gauss's Arithmetic-Geometric Mean
Additional background
Checklist
The proposed changes: