fourier bounce

desc/compute/_basis_vectors.py

@@ -3223,17 +3223,61 @@ def _e_sub_zeta_zz(params, transforms, profiles, data, **kwargs):
     transforms={},
     profiles=[],
     coordinates="rtz",
-    data=["e^rho", "e^theta", "e^zeta", "alpha_r", "alpha_t", "alpha_z"],
+    data=["periodic(grad(alpha))", "secular(grad(alpha))"],
 )
 def _grad_alpha(params, transforms, profiles, data, **kwargs):
-    data["grad(alpha)"] = (
-        data["alpha_r"] * data["e^rho"].T
+    data["grad(alpha)"] = data["periodic(grad(alpha))"] + data["secular(grad(alpha))"]
+    return data
+
+
+@register_compute_fun(
+    name="periodic(grad(alpha))",
+    label="\\mathrm{periodic}(\\nabla \\alpha)",
+    units="m^{-1}",
+    units_long="Inverse meters",
+    description=(
+        "Gradient of field line label, which is perpendicular to the field line, "
+        "periodic component"
+    ),
+    dim=3,
+    params=[],
+    transforms={},
+    profiles=[],
+    coordinates="rtz",
+    data=["e^rho", "e^theta", "e^zeta", "periodic(alpha_r)", "alpha_t", "alpha_z"],
+)
+def _periodic_grad_alpha(params, transforms, profiles, data, **kwargs):
+    data["periodic(grad(alpha))"] = (
+        data["periodic(alpha_r)"] * data["e^rho"].T
         + data["alpha_t"] * data["e^theta"].T
         + data["alpha_z"] * data["e^zeta"].T
     ).T
     return data
 
 
+@register_compute_fun(
+    name="secular(grad(alpha))",
+    label="\\mathrm{secular}(\\nabla \\alpha)",
+    units="m^{-1}",
+    units_long="Inverse meters",
+    description=(
+        "Gradient of field line label, which is perpendicular to the field line, "
+        "periodic component"
+    ),
+    dim=3,
+    params=[],
+    transforms={},
+    profiles=[],
+    coordinates="rtz",
+    data=["e^rho", "secular(alpha_r)"],
+)
+def _secular_grad_alpha(params, transforms, profiles, data, **kwargs):
+    data["secular(grad(alpha))"] = (
+        data["secular(alpha_r)"][:, jnp.newaxis] * data["e^rho"]
+    )
+    return data
+
+
 @register_compute_fun(
     name="grad(psi)",
     label="\\nabla\\psi",

desc/compute/_core.py

@@ -1518,14 +1518,48 @@ def _alpha(params, transforms, profiles, data, **kwargs):
     transforms={},
     profiles=[],
     coordinates="rtz",
-    data=["theta_PEST_r", "phi", "phi_r", "iota", "iota_r"],
+    data=["periodic(alpha_r)", "secular(alpha_r)"],
 )
 def _alpha_r(params, transforms, profiles, data, **kwargs):
-    data["alpha_r"] = (
-        data["theta_PEST_r"]
-        - data["iota_r"] * data["phi"]
-        - data["iota"] * data["phi_r"]
-    )
+    data["alpha_r"] = data["periodic(alpha_r)"] + data["secular(alpha_r)"]
+    return data
+
+
+@register_compute_fun(
+    name="periodic(alpha_r)",
+    label="\\mathrm{periodic}(\\partial_\\rho \\alpha)",
+    units="~",
+    units_long="None",
+    description="Field line label, derivative wrt radial coordinate, "
+    "periodic component",
+    dim=1,
+    params=[],
+    transforms={},
+    profiles=[],
+    coordinates="rtz",
+    data=["theta_PEST_r", "iota", "phi_r"],
+)
+def _periodic_alpha_r(params, transforms, profiles, data, **kwargs):
+    data["periodic(alpha_r)"] = data["theta_PEST_r"] - data["iota"] * data["phi_r"]
+    return data
+
+
+@register_compute_fun(
+    name="secular(alpha_r)",
+    label="\\mathrm{secular}(\\partial_\\rho \\alpha)",
+    units="~",
+    units_long="None",
+    description="Field line label, derivative wrt radial coordinate, "
+    "secular component",
+    dim=1,
+    params=[],
+    transforms={},
+    profiles=[],
+    coordinates="rtz",
+    data=["iota_r", "phi"],
+)
+def _secular_alpha_r(params, transforms, profiles, data, **kwargs):
+    data["secular(alpha_r)"] = -data["iota_r"] * data["phi"]
     return data
 
 

desc/compute/_metric.py

@@ -1941,8 +1941,8 @@ def _g_sup_ra(params, transforms, profiles, data, **kwargs):
     # Exact definition of the magnetic drifts taken from
     # eqn. 48 of Introduction to Quasisymmetry by Landreman
     # https://tinyurl.com/54udvaa4
-    label="\\mathrm{gbdrift} = 1/B^{2} (\\mathbf{b}\\times\\nabla B) \\cdot"
-    + "\\nabla \\alpha",
+    label="(\\nabla \\vert B \\vert)_{\\mathrm{drift}} = "
+    "(\\mathbf{b} \\times \\nabla B) \\cdot \\nabla \\alpha / \\vert B \\vert^{2}",
     units="1/(T-m^{2})",
     units_long="inverse Tesla meters^2",
     description="Binormal component of the geometric part of the gradB drift"
@@ -1952,13 +1952,61 @@ def _g_sup_ra(params, transforms, profiles, data, **kwargs):
     transforms={},
     profiles=[],
     coordinates="rtz",
-    data=["|B|^2", "b", "grad(alpha)", "grad(|B|)"],
+    data=["periodic(gbdrift)", "secular(gbdrift)"],
 )
 def _gbdrift(params, transforms, profiles, data, **kwargs):
-    data["gbdrift"] = (
-        1
+    data["gbdrift"] = data["periodic(gbdrift)"] + data["secular(gbdrift)"]
+    return data
+
+
+@register_compute_fun(
+    name="periodic(gbdrift)",
+    # Exact definition of the magnetic drifts taken from
+    # eqn. 48 of Introduction to Quasisymmetry by Landreman
+    # https://tinyurl.com/54udvaa4
+    label="\\mathrm{periodic}(\\nabla \\vert B \\vert)_{\\mathrm{drift}}",
+    units="1/(T-m^{2})",
+    units_long="inverse Tesla meters^2",
+    description="Binormal component of the geometric part of the gradB drift"
+    + " used for local stability analyses, Gamma_c, epsilon_eff etc."
+    " Periodic component.",
+    dim=1,
+    params=[],
+    transforms={},
+    profiles=[],
+    coordinates="rtz",
+    data=["|B|^2", "b", "periodic(grad(alpha))", "grad(|B|)"],
+)
+def _periodic_gbdrift(params, transforms, profiles, data, **kwargs):
+    data["periodic(gbdrift)"] = (
+        dot(data["b"], cross(data["grad(|B|)"], data["periodic(grad(alpha))"]))
+        / data["|B|^2"]
+    )
+    return data
+
+
+@register_compute_fun(
+    name="secular(gbdrift)",
+    # Exact definition of the magnetic drifts taken from
+    # eqn. 48 of Introduction to Quasisymmetry by Landreman
+    # https://tinyurl.com/54udvaa4
+    label="\\mathrm{secular}(\\nabla \\vert B \\vert)_{\\mathrm{drift}}",
+    units="1/(T-m^{2})",
+    units_long="inverse Tesla meters^2",
+    description="Binormal component of the geometric part of the gradB drift"
+    + " used for local stability analyses, Gamma_c, epsilon_eff etc. "
+    "Secular component.",
+    dim=1,
+    params=[],
+    transforms={},
+    profiles=[],
+    coordinates="rtz",
+    data=["|B|^2", "b", "secular(grad(alpha))", "grad(|B|)"],
+)
+def _secular_gbdrift(params, transforms, profiles, data, **kwargs):
+    data["secular(gbdrift)"] = (
+        dot(data["b"], cross(data["grad(|B|)"], data["secular(grad(alpha))"]))
         / data["|B|^2"]
-        * dot(data["b"], cross(data["grad(|B|)"], data["grad(alpha)"]))
     )
     return data
 
@@ -1979,11 +2027,35 @@ def _gbdrift(params, transforms, profiles, data, **kwargs):
     transforms={},
     profiles=[],
     coordinates="rtz",
-    data=["p_r", "psi_r", "|B|^2", "gbdrift"],
+    data=["periodic(cvdrift)", "secular(gbdrift)"],
 )
 def _cvdrift(params, transforms, profiles, data, **kwargs):
-    dp_dpsi = mu_0 * data["p_r"] / data["psi_r"]
-    data["cvdrift"] = 1 / data["|B|^2"] * dp_dpsi + data["gbdrift"]
+    data["cvdrift"] = data["periodic(cvdrift)"] + data["secular(gbdrift)"]
+    return data
+
+
+@register_compute_fun(
+    name="periodic(cvdrift)",
+    # Exact definition of the magnetic drifts taken from
+    # eqn. 48 of Introduction to Quasisymmetry by Landreman
+    # https://tinyurl.com/54udvaa4
+    label="\\mathrm{periodic(cvdrift)}",
+    units="1/(T-m^{2})",
+    units_long="inverse Tesla meters^2",
+    description="Binormal component of the geometric part of the curvature drift"
+    + " used for local stability analyses, Gamma_c, epsilon_eff etc. "
+    "Periodic component.",
+    dim=1,
+    params=[],
+    transforms={},
+    profiles=[],
+    coordinates="rtz",
+    data=["p_r", "psi_r", "|B|^2", "periodic(gbdrift)"],
+)
+def _periodic_cvdrift(params, transforms, profiles, data, **kwargs):
+    data["periodic(cvdrift)"] = (
+        mu_0 * data["p_r"] / data["psi_r"] / data["|B|^2"] + data["periodic(gbdrift)"]
+    )
     return data
 
 

desc/grid.py

@@ -671,6 +671,7 @@ class Grid(_Grid):
         Use np.inf to denote no periodicity.
     NFP : int
         Number of field periods (Default = 1).
+        Change this only if your nodes are placed within one field period.
     source_grid : Grid
         Grid from which coordinates were mapped from.
     sort : bool
@@ -794,7 +795,8 @@ def create_meshgrid(
         NFP : int
             Number of field periods (Default = 1).
             Only makes sense to change from 1 if last coordinate is periodic
-            with some constant divided by ``NFP``.
+            with some constant divided by ``NFP`` and the nodes are placed
+            within one field period.
 
         Returns
         -------
@@ -916,6 +918,8 @@ class LinearGrid(_Grid):
         Toroidal grid resolution.
     NFP : int
         Number of field periods (Default = 1).
+        Change this only if your nodes are placed within one field period
+        or should be interpreted as spanning one field period.
     sym : bool
         True for stellarator symmetry, False otherwise (Default = False).
     axis : bool
@@ -1011,6 +1015,8 @@ def _create_nodes(  # noqa: C901
             Toroidal grid resolution.
         NFP : int
             Number of field periods (Default = 1).
+            Only change this if your nodes are placed within one field period
+            or should be interpreted as spanning one field period.
         axis : bool
             True to include a point at rho=0 (default), False for rho[0] = rho[1]/2.
         endpoint : bool
@@ -1037,8 +1043,10 @@ def _create_nodes(  # noqa: C901
         """
         self._NFP = check_posint(NFP, "NFP", False)
         self._period = (np.inf, 2 * np.pi, 2 * np.pi / self._NFP)
-        # TODO:
+        # FIXME:
         #  https://github.com/PlasmaControl/DESC/pull/1204#pullrequestreview-2246771337
+        #  Quantities like alpha, grad(alpha), etc. are computed incorrectly at
+        #  phi > 2pi / NFP and theta > 2pi / NFP.
         axis = bool(axis)
         endpoint = bool(endpoint)
         theta_period = self.period[1]

desc/integrals/__init__.py

@@ -1,6 +1,6 @@
 """Classes for function integration."""
 
-from .bounce_integral import Bounce1D
+from .bounce_integral import Bounce1D, Bounce2D
 from .singularities import (
     DFTInterpolator,
     FFTInterpolator,