Magnetic axis limits and new variables

[unalmis]

Math to compute the rotational transform in DESC.pdf

• Resolves Implement correct axis limits #486
• Update desc.plotting.plot_fsa to be able to plot at the axis
• Remove desc.compute.compress and desc.compute.expand functions in favor grid.compress and grid.expand object methods
• Miscellaneous formatting changes caught when doing above tasks
• Update test durations to prevent GitHub CI timeout failure

If the limits were unable to be computed with the version of l’Hôpital’s below, then I probably added some comments to outline whatever alternative technique I used. Please note that this version of l’Hôpital’s rule I have used can only resolve the $\mathbf{0}/0$ indeterminate form not the (something / $\infty$) form.

An $ℝ^n → ℝ^m$ l’Hôpital’s rule for $\mathbf{0}/0$

Part of my work for this PR included proving non-existence of the limit of some of the more complicated functions. (Recall that for a function $a/b$, if $da/db$ is proven to not exist, we still know nothing about the limit of $a/b$). The functions whose limits do not exist are listed at the top of tests/test_axis_limits.py. If you don't want to confirm this non-existence with math then a heuristic alternative is just plotting (script below) and seeing if they blow up. For the limits that do exist, I used the rule below. I didn't think it's widely known that you can use a l’Hôpital’s rule for our functions, so I wrote this explanation.

l’Hôpital’s rule relies on Cauchy’s generalized mean value theorem, which is a property proven unique to functions that map $ℝ → ℝ$ (see paper 1). The quantities we compute map $ℝ^3 → ℝ^m$, so it would be very wrong to evaluate limits by comparing the end behavior of their derivatives. Our polar coordinates do allow us to write the limit we desire as $\lim_{\rho \to 0} f(\rho, \theta, \zeta)$, but this does not let us reduce our limits to single variable limits. See here for an already written explanation reminding this.

Fortunately, some parts of l’Hôpital’s rule have (recently!) been generalized to multivariable functions that map $ℝ^n → ℝ$ where $𝑛 ∈ ℕ$. Since all of our singularities are non-isolated (every ball that contains a point of the magnetic axis contains more than one point on the magnetic axis) and, when the quantity can be written as a fraction, the zero set of the numerator contains the zero set of the denominator inside some small tube surrounding the magnetic axis, we are able to apply theorem 4, claim 1 of paper 2: multivariable form of l’Hôpital’s rule for our scalar valued functions that have the $0/0$ indeterminate form. The linked paper does not generalizes this result for vector-valued multivariable functions from $ℝ^n → ℝ^m$ -- which is a result we require; more on this below.

The fundamental reason why generalizing a l’Hôpital rule to a vector-valued function is not possible with a mean value theorem is because for
$$\frac{X}{Y}(\rho) = \left( \frac{x_1}{y_1}(\rho), \frac{x_2}{y_2}(\rho), \frac{x_3}{y_3}(\rho) \right)$$, while you may be able to argue something like
$$\left\lvert \lim \frac{X}{Y}(\rho) \right\rvert = \left\lvert \left(\lim \frac{x_1}{y_1}, \lim \frac{x_2}{y_2}, \lim \frac{x_3}{y_3} \right) \right\rvert = \left\lvert (\lim \frac{dx_1/d\rho}{dy_1/d\rho}, \lim \frac{dx_2/d\rho}{dy_2/d\rho}, \lim \frac{dx_3/d\rho}{dy_3/d\rho}) \right\rvert$$
assuming l’Hôpital’s rule applies to each component, the direction of the limit of the vector is unresolved because the mean value theorem cannot guarantee that the three independent limits are to be evaluated at the same point. So all you know is $$(\lim \frac{x_1}{y_1}(\rho), \lim \frac{x_2}{y_2}(\rho), \lim \frac{x_3}{y_3}(\rho)) = (\frac{dx_1/d\rho}{dy_1/d\rho}(\rho=\xi_1),\frac{dx_2/d\rho}{dy_2/d\rho}(\rho=\xi_2), \frac{dx_3/d\rho}{dy_3/d\rho}(\rho=\xi_3)) $$
The right hand side is equal to $\frac{d X / d \rho}{d Y / d\rho}$ if somehow $\xi_1=\xi_2=\xi_3$.
Edit: The wikipedia page explains similar.

The earlier linked paper 1 takes a more topological approach to sidestep this issue and recovers l’Hôpital's rule for maps from $ℝ → ℝ^m$, where the additional constraint that the sequence in the denominator be the same function from $ℝ → ℝ$ for all components. However, their proof relies on a property that is unique to the domain $ℝ$. It does not generalize to our multivariable setting $ℝ^3 → ℝ^m$.

Generalizing l’Hôpital’s rule to the multivariable real-vector-valued setting

If we can (I did) recast all our vector-valued indeterminacy into the form $\mathbf{0}/0$, then we can generalize the multivariable l’Hôpital’s rule to evaluate our limit. My claim for this is as follows. We know that a vector valued sequence $X$ converges to the vector that is the limit of each component if and only if the limit of each component exists. Now if $(x_1, y_1, x_2, y_2, x_3, y_3)$ converges to a tuple of zeros as $\rho \to 0$, then we can write
$$\lim \frac{X}{Y} = (\lim \frac{x_1}{y_1} , \lim \frac{x_2}{y_2} , \lim \frac{x_3}{y_3} ) = (\frac{\lim \frac{x_1 - 0}{\lvert \rho - 0 \rvert}}{\lim \frac{y_1-0}{\lvert \rho - 0 \rvert}}, \frac{\lim \frac{x_2 - 0}{\lvert \rho - 0 \rvert}}{\lim \frac{y_2-0}{\lvert \rho - 0 \rvert}}, \frac{\lim \frac{x_3 - 0}{\lvert \rho - 0 \rvert}}{\lim \frac{y_3-0}{\lvert \rho - 0 \rvert}} ) = \left(\frac{\partial x_1/ \partial\rho}{\partial y_1/ \partial \rho} , \frac{\partial x_2/ \partial \rho}{\partial y_2/ \partial \rho} , \frac{\partial x_3/ \partial\rho}{\partial y_3/ \partial\rho} \right) = \frac{\partial X / \partial \rho}{\partial Y / \partial \rho} $$
The second to last equality follows if we take the limit along the path that holds all other variables except $\rho$ as constant. Now if we know the limits of each of the components exist (not just along the path that only varies $\rho$), then, for any component, the path taken varying only $\rho$ above must be the limit.

follow up comment here: #556 (comment)

(Details...
If you're wondering, this argument is "necessary" in the sense that the result does not follow from just applying the multivariable l’Hôpital’s rule to each component. The latter approach tells you that the limit is the vector of quotients of partial derivatives where each component need not be differentiated with respect to the same variable. The above shows that existence of the limit of each component implies the limit of the vector has all components as quotients of partial derivatives with respect to the same variable.
)

scripts may be useful for reviewers

from desc.grid import LinearGrid
from desc.examples import get
from desc.backend import jnp
from desc.compute.utils import surface_integrals_map
from desc.compute import data_index
from tests.test_axis_limits import not_finite_limits
import matplotlib.pyplot as plt


eq = get("W7-X")  # should also try QAS or other fixed-current
rho = jnp.linspace(0, 0.5, 30) ** 1.5
grid = LinearGrid(rho=rho, M=8, N=8, NFP=eq.NFP)
data = eq.compute(list(not_finite_limits), grid=grid)
integrate = surface_integrals_map(grid, expand_out=False)

for name in not_finite_limits:
    # make single variable function of rho
    if data_index["desc.equilibrium.equilibrium.Equilibrium"][name]["coordinates"] == "r":
        # already function of rho
        profile = grid.compress(data[name])
    else:
        rho_surface_has_nan = integrate(jnp.isnan(data[name])).astype(bool)
        profile = jnp.where(
            rho_surface_has_nan,
            # integration below replaced nans with 0, put them back
            jnp.nan,
            # Norms and integrals are continuous functions, so their composition
            # cannot disrupt existing continuity. Note that the absolute value
            # before the integration ensures that a discontinuous integrand does
            # not become continuous once integrated.
            integrate(jnp.abs(data[name])),
        )
    fig, ax = plt.subplots()
    ax.plot(rho, profile, marker='o')
    ax.set_title(name)
    plt.pause(0.1)  # show the plot before requesting user input
    input("Press Enter to continue...")
    plt.close(fig)

mathematica notebook to help do limit computations:

lhop limits.zip

[dpanici]

I will check after this and the coils are merged in, or before I merge coils in

[dpanici]

I can't seem to see the diff, or the file in the tests/baseline folder, for the following plot tests:
They must have been inadvertently deleted?

• test_plot_vac_normF_2d.png
• test_plot_vac_normF_section.png

[f0uriest]

I can't seem to see the diff, or the file in the tests/baseline folder, for the following plot tests:
They must have been inadvertently deleted?

• test_plot_vac_normF_2d.png
• test_plot_vac_normF_section.png

They're both there, but for some reason GitHub sometimes doesn't show image diffs

[unalmis]

I can't seem to see the diff, or the file in the tests/baseline folder, for the following plot tests: They must have been inadvertently deleted?

• test_plot_vac_normF_2d.png
• test_plot_vac_normF_section.png

They were deleted because the tests that generated those outputs were renamed in earlier PRs. I think the new outputs don't have the vac in their name.

[unalmis]

They were deleted because the tests that generated those outputs were renamed in earlier PRs. I think the new outputs don't have the vac in their name.

My guess is something went wrong with a git merge because Dario added these in 6e7161b, but the git history on master never shows the tests in test/plotting.py.

[dpanici]

Hm I see, yea that is fine, the tests dont actually exist so no need for the pngs.