Adding bounce-averaging functionality

[rahulgaur104]

Bounce=averaging is used in many different calculations in stellarator optimization. Some of the examples are the $\epsilon_{eff}$ proxy, the energetic particle proxy $\Gamma_c$, neoclassical solvers, and gyrokinetic solvers (such as GS2).

Hence, it would be useful to add that feature to DESC. Below is an idea that I have used in the past to calculate bounce-averaged drift frequency of electrons $\omega_{De}$. I think we can use the same logic for all bounce averages.

Below I have briefly explained how to apply the Gauss-Chebyshev quadrature to calculate the bounce-averaged integrals of the form $$\langle X \rangle = \left(\int_{-\theta_b}^{\thetab} \frac{1}{\mathbf{B}\cdot\mathbf{\nabla}\theta}\frac{d\theta}{\sqrt{1-\lambda B}}\right)^{-1} \int{-\theta_b}^{\theta_b} \frac{1}{\mathbf{B}\cdot\mathbf{\nabla}\theta}\frac{d\theta\, f(\theta)}{\sqrt{1-\lambda B}}$$ where $\lambda$ is the pitch angle, $B$ is the magnetic field strength along a field line on a given flux surface, and $f(\theta)$ is a function dependent on the equilibrium geometry. Making the use of the transformation $$x = \sin\left(-\frac{\theta}{\thetab} \frac{\pi}{2}\right)$$ we can write the bounce average as $$\langle X \rangle = \left(\int{-1}^{1} \frac{dx\, X2(x)}{\sqrt{1-x^2}}\right)^{-1} \int{-1}^{1} \frac{dx\, X_1(x)}{\sqrt{1-x^2}}$$ where $$X_1 = \frac{1}{\mathbf{B}\cdot\mathbf{\nabla}\theta}\frac{f(x)}{\sqrt{1-\lambda B}}$$ $$X2 = \frac{1}{\mathbf{B}\cdot\mathbf{\nabla}\theta}\frac{1}{\sqrt{1-\lambda B}}$$ Using the Gauss-Chebyshev quadrature rule, we can write $$\langle X \rangle = \left(\sum{i=1}^{n} w_i X_2(xi) \right)^{-1} \sum{i=1}^{n} w_i X_1(x_i)$$ where $w_i$ are fixed weights associated with the Gauss-Chebyshev quadrature and $x_i$ are set of points where we should evaluate the function $X_1(xi)$. For a up-down asymmetric or a general magnetic field profile with bounce points that lie between $\theta{b1}$ and $\theta{b2}$, we may first use the transformation $$x = \sin\left[\left(1 + 2\frac{(\theta - \theta{b2})}{\theta{b2}-\theta{b1}}\right)\frac{\pi}{2}\right]$$ This will transform an integral of the type $$\langle X \rangle = \left(\int{\theta{b1}}^{\theta{b2}} \frac{1}{\mathbf{B}\cdot\mathbf{\nabla}\theta}\frac{d\theta}{\sqrt{1-\lambda B}}\right)^{-1} \int{\theta{b1}}^{\theta{b2}} \frac{1}{\mathbf{B}\cdot\mathbf{\nabla}\theta}\frac{d\theta\, f(\theta)}{\sqrt{1-\lambda B}}$$ into the following integral $$\langle X \rangle = \left(\int_{-1}^{1} \frac{dx\, X2(x)}{\sqrt{1-x^2}}\right)^{-1} \int{-1}^{1} \frac{dx\, X_1(x)}{\sqrt{1-x^2}}$$ where $$X_1 = \frac{1}{\mathbf{B}\cdot\mathbf{\nabla}\theta}\frac{f(x)}{\sqrt{1-\lambda B}}$$ $$X_2 = \frac{1}{\mathbf{B}\cdot\mathbf{\nabla}\theta}\frac{1}{\sqrt{1-\lambda B}}$$ and we can use the same quadrature rule to obtain the required integral. For a magnetic field, one would have to find all the intersection points of the pitch angle lambda with the function $B(\theta)$. This has been done before in branch jc/gammaC. I have a different script for this. A typical set of intersection points are shown below

For each horizontal line, we find all the intersection points and calculate the integral between each pair of bounce points. This should be enough to calculate the various metrics mentioned above and hopefully do some optimization.

[f0uriest]

I'm a bit confused by the use of gauss-chebyshev above. You factor out the $1/\sqrt(1-x^2)$ but the leftover part still has a $1/\sqrt(1-\lambda B)$ which is still singular at the endpoints. Even if you avoid evaluating it at the endpoints, its still a singular integrand which likely won't be well approximated by a polynomial, unless there's something I'm missing?

[rahulgaur104]

Yes, the term 1/sqrt(1-lambda B) is singular at the endpoints but by virtue of Gauss-Chebyshev quadrature, none of the quadrature points will lie at the singular points. You can get arbitrarily close without actually being at the singular point.

The objective is not to transform the singularity into some nice, non-singular integral but to use a robust quadrature that excludes the singular points without excluding the contributions of the integrand near the integrable singularities.

[f0uriest]

hmmm ok. Once we get the basics working it would be good to compare gauss-chebyshev vs tanh-sinh or modified trapezoid.

[rahulgaur104]

https://cptc.wisc.edu/wp-content/uploads/sites/327/2017/09/UW-CPTC_15-4.pdf

Paper with analytical expressions for bounce-averaged drifts.

[rahulgaur104]

Here's a detailed comparison of these quadratures for a simple shifted circle case.

• At lower resolutions Gauss-Chebyshev is more accurate, especially for the integral of the type 1/sqrt(1-lambda B).
• At each resolution, sinh-tanh quadrature takes around 40% longer than Gauss-Chebyshev. Most of the time is spent to calculate the weights and quadrature points so that can just be done once for the whole calculation.
• Sinh-tanh is clearly a better method but I would still like to keep both quadratures since there may be differences when we do this for a more complicated geometry.
• The transformation from [-theta_b, theta_b] -> [-1, 1] can have many different forms, I chose a sin function but you can choose any odd function of the form x^{n} or x^{1/n} where n is an odd number or any invertible function. One can use arcsin if they want.

[unalmis]

I think the better transformation would be arcsin. The Chebyshev and tanh sinh quadratures are built to integrate singularities at boundaries, but obviously the less singular the better for faster convergence.

Using the sin transformation as above introduces an additional singular factor of 1 / sqrt(1 - x^2). As a higher order effect, it also slightly pushes the quadrature points x_k further toward the singularities, which seems undesirable as the closer the evaluation points are to the singularities the more punishing the round off error will be when evaluating 1 / sqrt(1 - x^2) and the bounce integrand.

An arcsin transformation on the other hand introduces a factor proportional to cosine that will compete with the singularity of the bounce integral, making it easier for the quadrature to approximate. And it contracts the quadrature points toward regions that induce less round off error.

Basically we should pick a change of variable that pulls the mass of the integral away from the singularity rather than pushing more toward it.

[unalmis]

So I generalized the bounce integral PR to allow for any transformation. As I mentioned above the sin automorphism works well to suppress singularities. In fact, the transformation with arcsin plus Gauss-Legendre quadrature beats the transformation with sin and tanh-sinh quadrature (16 nodes to 18 nodes)

[rahulgaur104]

The analytical bounce averaging turned out to be more tedious than I expected but here are the 8 terms we need to calculate and all the integrals here are called incomplete elliptic functions that can be calculated using scipy.

bavg_notes.pdf

@unalmis Sorry for the delay. I'll upload the modified test today but it will still take this week to correct signs, factors of 2 in these analytical terms. But that's it, after that there is no more analytical work needed for this test other than debugging.