Gamma_c (Nemov and Velasco et. al) Γ_c

[unalmis]

Resolves #522 . Automatically differentiable

[github-actions]

|             benchmark_name             |         dt(%)          |         dt(s)          |        t_new(s)        |        t_old(s)        | 
| -------------------------------------- | ---------------------- | ---------------------- | ---------------------- | ---------------------- |
 test_build_transform_fft_lowres         |     +0.38 +/- 6.30     | +2.06e-03 +/- 3.40e-02 |  5.42e-01 +/- 2.1e-02  |  5.40e-01 +/- 2.7e-02  |
 test_build_transform_fft_midres         |     -0.95 +/- 6.07     | -5.71e-03 +/- 3.66e-02 |  5.98e-01 +/- 2.6e-02  |  6.04e-01 +/- 2.5e-02  |
 test_build_transform_fft_highres        |     +1.18 +/- 4.91     | +1.17e-02 +/- 4.85e-02 |  1.00e+00 +/- 3.1e-02  |  9.88e-01 +/- 3.7e-02  |
 test_equilibrium_init_lowres            |     +2.01 +/- 3.27     | +7.54e-02 +/- 1.22e-01 |  3.82e+00 +/- 9.7e-02  |  3.75e+00 +/- 7.5e-02  |
 test_equilibrium_init_medres            |     -0.85 +/- 9.48     | -3.82e-02 +/- 4.28e-01 |  4.47e+00 +/- 2.7e-01  |  4.51e+00 +/- 3.3e-01  |
 test_equilibrium_init_highres           |     -0.95 +/- 2.57     | -5.84e-02 +/- 1.57e-01 |  6.07e+00 +/- 1.3e-01  |  6.13e+00 +/- 8.1e-02  |
 test_objective_compile_dshape_current   |     -0.33 +/- 6.32     | -1.21e-02 +/- 2.28e-01 |  3.60e+00 +/- 1.6e-01  |  3.61e+00 +/- 1.6e-01  |
 test_objective_compile_atf              |     +1.50 +/- 2.09     | +1.08e-01 +/- 1.50e-01 |  7.28e+00 +/- 1.2e-01  |  7.18e+00 +/- 8.9e-02  |
 test_objective_compute_dshape_current   |     +2.25 +/- 2.28     | +8.13e-05 +/- 8.23e-05 |  3.69e-03 +/- 7.3e-05  |  3.61e-03 +/- 3.7e-05  |
 test_objective_compute_atf              |     -0.05 +/- 2.99     | -8.94e-06 +/- 4.86e-04 |  1.63e-02 +/- 2.7e-04  |  1.63e-02 +/- 4.0e-04  |
 test_objective_jac_dshape_current       |     -0.95 +/- 8.16     | -3.94e-04 +/- 3.37e-03 |  4.09e-02 +/- 1.9e-03  |  4.13e-02 +/- 2.8e-03  |
 test_objective_jac_atf                  |     -1.67 +/- 2.60     | -3.29e-02 +/- 5.11e-02 |  1.93e+00 +/- 4.0e-02  |  1.96e+00 +/- 3.2e-02  |
 test_perturb_1                          |     -4.11 +/- 3.79     | -5.25e-01 +/- 4.84e-01 |  1.22e+01 +/- 4.6e-02  |  1.28e+01 +/- 4.8e-01  |
 test_perturb_2                          |     -2.46 +/- 2.88     | -4.41e-01 +/- 5.16e-01 |  1.75e+01 +/- 4.3e-01  |  1.79e+01 +/- 2.9e-01  |
 test_proximal_jac_atf                   |     -0.40 +/- 1.38     | -3.01e-02 +/- 1.04e-01 |  7.49e+00 +/- 5.4e-02  |  7.52e+00 +/- 8.9e-02  |
 test_proximal_freeb_compute             |     -1.48 +/- 1.43     | -2.91e-03 +/- 2.81e-03 |  1.93e-01 +/- 2.0e-03  |  1.96e-01 +/- 1.9e-03  |
 test_proximal_freeb_jac                 |    -12.25 +/- 14.41    | -1.05e+00 +/- 1.24e+00 |  7.54e+00 +/- 2.6e-01  |  8.59e+00 +/- 1.2e+00  |

[f0uriest]

From looking through STELLOPTs version of Gamma_C here I see a few possible issues with their implementation that might lead to differences:

• Ntransits and Nzetapertransit hard coded to 400 each, and npitch hard coded to 80, so will need to recompile if we want to scan resolution (from a comment it looks like ROSE uses nzetapertransit=290?)
• They don't seem to do any root finding to get bounce points, or even interpolation. They take the bounce points as the first and last zeta points in the fixed zeta grid where B < 1/lambda. I think this means they probably miss a lot of the contribution to Gamma_C near the tops of wells
• The integration (over both lambda and zeta) seems to only be a first order Reimann sum so may need super high resolution.

In contrast, the WISTELL implementation here seems a lot better

• They spline qtys along a field line using cubic splines
• bounce points are found through rootfinding on these splines
• bounce integrals are done using adaptive tanh-sinh quadrature over splines
• integral over lambda also seems to be done with adaptive gauss-kronrod, by splining gamma_c(lambda)
• seems like they also can calculate eps_eff

Might be worth reaching out to Aaron and Ben to talk about benchmarking vs WISTELL

[unalmis]

From looking through STELLOPTs version of Gamma_C here

Thanks! I know code is not always a leisure read. I agree it seems problematic. The integrals are strong functions of the bounce points; so it's necessary to obtain them accurately. A specialized quadrature is also recommended. Prior to writing the DESC bounce integrals, I tested this (because if we could just run along the field line using a reiman sum, less expensive algorithms could be made to compute this). Assuming the bounce points are known with accuracy, a non-adaptive trapezoidal scheme to evaluate these bounce integrals required orders of magnitude more quadrature points. Might become a race with the floating point plateau.

Note that there are two STELLOPT implementations; the second one is better. They mention use of tanh-sinh, but it's been a while since I looked at it, and I don't remember where it used now.

[f0uriest]

Yeah I looked at that as well but it seems the same? There are more comments explaining some of the coordinate transformations but the bounce points are still rounded to the nearest grid point (there is some stuff about root finding but it's inside a "if 0 > 1" block), and the following stuff doesn't seem to use it. Also looks like it's still just doing first order integration over both zeta and pitch, there's a lone comment about tanh sinh but I didn't see anything that actually looked like it was being used

[unalmis]

the WISTELL implementation

I'll look over this more thoroughly later, but from a quick glance, some points to consider

integral over lambda also seems to be done with adaptive gauss-kronrod, by splining gamma_c(lambda)

I see that a cubic spline is fit to the integrand, i.e. onto $f$ in $\int f(\lambda) d\lambda$, and then an adaptive quadrature is used to integrate this spline. Simpson's rule integrates a cubic spline exactly, so I don't understand that strategy.

Likewise the integrand $f(\lambda)$ is not smooth. Each derivative of $f(\lambda)$ is less smooth than its primitive, and for $\Gamma_c$, we have $f(\lambda) \to \infty$ if it intersects a local maximum of |B| (unless there is some magical cancellation between the diverging integrals among the different ripples), so any usual function approximation on $f$ will underestimate the integral.

bounce integrals are done using adaptive tanh-sinh quadrature over splines

FYI, some reasons I chose not to use an adaptive quadrature for the bounce integrals

1. A fixed quadrature makes the computation homogenous, and this will improve performance on GPUs/numpy/jax etc.
2. The goal is computing a double integral (1. over pitch $\lambda$ and 2. some bounce integrals). For $\Gamma_c$, in theory, my hunch is an adaptive quadrature on the bounce integrals will over-fit on the diverging integrals and require more pitch values to capture the true double integral
3. it's more expensive

[unalmis]

Unrelated, but since we're discussing quadrature, I'm not a fan of regular tanh-sinh because it packs nodes too densely near bounce points. This could make floating point errors in the bounce points punishing. Nemov's effective ripple paper mentions how there's additional phenomena to capture by considering particles trapped in wells with ripples, so intuitively that means when we do these bounce integrals we should try to capture more than just the features near the bounce points (hence reduce clustering of quadrature near bounce points as much as possible once the singularity is resolved).

[unalmis]

In previous commits, I had an kwarg option to perform the integral over $\lambda$ with any adaptive quadrature from quadax (without splining the integrand, so the functions behavior is truly captured). Do we want this option back?