Pseudo-spectral improvements for bounce integrals over surfaces

[unalmis]

Basically make effective ripple /gamma c type optimization faster/more memory efficient.

some thoughts: To evaluate a quantity along a field line, one can generate a grid in the Clebsh-Type field line coordinate system $\rho, \alpha, \zeta$, then vary $\zeta$ at fixed $\rho, \alpha$. We can then transform into the DESC computational coordinates to evaluate Recall $\alpha$ is multivalued in this coordinate system; after a toroidal transit, the field line is still specified by $\alpha$, but unless it closes on itself, it is at a different physical poloidal coordinate. It would be useful to have a utility that maps a field line onto a rectangular grid, cutting the field line after each toroidal transit.

Applications

1. bounce integration/ neoclassical stuff

Currently, we generate one-dimensional splines of ‖𝐵‖ along field lines. In stellarators, there can be fine ripples in ‖𝐵‖(𝜁) which are only captured by splines at high knot density. Here fine ripple means the distance between bounce points is small compared to the number of toroidal transits. We are integrating over a large domain 𝜁 ∈ 𝐼 ⊂ ℝ with more accurate surface averages obtained as 𝐼 → ℝ. In that limit, it becomes infeasible to maintain the knot density required to capture all ripples, as the total number of required knots grows with |𝐼|. Since the field line wraps around a surface, the knots that support reconstruction of ‖𝐵‖ at small 𝜁 should also support ‖𝐵‖ at large 𝜁. Indeed, this is the job of a 2D spline in $\alpha, \zeta$, and switching to such an interpolation will detach the required number of knots from the length of the field line integration.

Note that we don’t want a 3D spline; quantities that rely on bounce integrals to compute usually reduce to flux surface functions. An optimization will involve selecting a finite number of flux surfaces over which to optimize the quantity. Unless the radial resolution is also high, a 3D spline would hinder reconstruction of ‖𝐵‖ over a surface.

Implementation suggestions:

1. Can compute the jump in $\alpha$ at the boundary of the domain at $\zeta = 0 = 2 \pi$ given the rotational transform.
2. So can generate a set $A_1$ of alpha values and the corresponding field lines from $\zeta \in [0, 2\pi]$. Can generate another set $A_2 = A_1 + \alpha_{\text{jump}}$ where $\alpha_{\text{jump}}$ is computed from the rotational transform, and so on.
3. Use interpax to compute a 2D spline of ‖𝐵‖($\alpha$, 𝜁). For a fixed $\alpha$, the roots of this spline can be computed analytically, in the same manner currently implemented on bounce branch. For bounce points across the "branch" cut at $\zeta = 2 \pi$, treat the boundary $\zeta = 2\pi$ as a right bounce point for all field lines in $A_1$ and $\zeta = 0$ as the left bounce point for all field lines in $A_2$. This is fine because our integration technique works fine for singular and non-singular integrals.

Might need to stagger the knots in $\zeta$ across $\alpha$ to save memory. So if the $\alpha = 0$ field line has knots at every 1/10 increment in $\zeta$ then the $\alpha = 0 + \Delta \alpha $ field line should have knots shifted by some amount to lie in between the knots of the other field line.

[f0uriest]

A few comments:

1. to your point 3 above, interpax can do 2d splines just fine, as long as its over a tensor product grid of $\alpha, \zeta$, so your note about staggering grids might not work.
2. Another option might be to compute $f(\theta, \zeta)$ and $\alpha(\theta, \zeta)$ and spline both, then when evaluating $f(\alpha, \zeta)$ first do root-finding on $\alpha(\theta, \zeta)$ to find $\theta$ then evaluate the $f(\theta, \zeta)$ spline. This would have the advantage of only needing to evaluate things on a fixed grid in regular DESC coordinates, and the only root-finding is done on a spline which is much cheaper. Additionally, since both $f$ and $\alpha$ are periodic in DESC coordinates, you could use FFTs to upsample before cubic interpolation to get much higher accuracy (or similarly use NUFFT for the interpolation)

[unalmis]

That's a good idea, thanks Rory. I've been doing some reading form Boyd's spectral book and thinking for the optimal method to compute these quantities over the last week. my plan is a little different than your suggestion, mainly ~because we can do this (epsilon effective) directly in DESC coordinates~. Much of the logic from PR #854 will be reused, with smarter interpolation methods.