Debugging singular integrals

CHANGELOG.md

@@ -25,6 +25,7 @@ Bug Fixes
 - Fixed I/O bug when saving/loading ``_Profile`` classes that do not have a ``_params`` attribute.
 - Minor bugs described in [#1323](https://github.com/PlasmaControl/DESC/pull/1323).
 - Corrects basis vectors computations made on surface objects [#1175](https://github.com/PlasmaControl/DESC/pull/1175).
+- Fix issue with interpolator for singular integrals [#1522](https://github.com/PlasmaControl/DESC/issues/1522) and additional checks [1519](https://github.com/PlasmaControl/DESC/issues/1519).
 
 v0.13.0
 -------

desc/basis.py

@@ -261,9 +261,9 @@ def _get_modes(self, L):
             Each row is one basis function with modes (l,m,n).
 
         """
-        l = np.arange(L + 1).reshape((-1, 1))
-        z = np.zeros((L + 1, 2))
-        return np.hstack([l, z])
+        l = np.arange(L + 1)
+        z = np.zeros_like(l)
+        return np.array([l, z, z]).T
 
     def evaluate(
         self, nodes, derivatives=np.array([0, 0, 0]), modes=None, unique=False
@@ -374,10 +374,9 @@ def _get_modes(self, N):
             Each row is one basis function with modes (l,m,n).
 
         """
-        dim_tor = 2 * N + 1
-        n = np.arange(dim_tor).reshape((-1, 1)) - N
-        z = np.zeros((dim_tor, 2))
-        return np.hstack([z, n])
+        n = np.arange(-N, N + 1)
+        z = np.zeros_like(n)
+        return np.array([z, z, n]).T
 
     def evaluate(
         self, nodes, derivatives=np.array([0, 0, 0]), modes=None, unique=False
@@ -400,6 +399,9 @@ def evaluate(
         -------
         y : ndarray, shape(num_nodes,num_modes)
             Basis functions evaluated at nodes.
+            The Vandermonde matrix when ``modes is None`` is
+            given by ``y.reshape(-1,2*N+1)`` and is ordered
+            [sin(N𝛇), ..., sin(𝛇), 1, cos(𝛇), ..., cos(N𝛇)].
 
         """
         if modes is None:
@@ -477,9 +479,7 @@ def __init__(self, M, N, NFP=1, sym=False):
         self._NFP = check_posint(NFP, "NFP", False)
         self._sym = bool(sym) if not sym else str(sym)
         self._spectral_indexing = "linear"
-
         self._modes = self._get_modes(M=self.M, N=self.N)
-
         super().__init__()
 
     def _get_modes(self, M, N):
@@ -499,16 +499,13 @@ def _get_modes(self, M, N):
             Each row is one basis function with modes (l,m,n).
 
         """
-        dim_pol = 2 * M + 1
-        dim_tor = 2 * N + 1
-        m = np.arange(dim_pol) - M
-        n = np.arange(dim_tor) - N
-        mm, nn = np.meshgrid(m, n)
-        mm = mm.reshape((-1, 1), order="F")
-        nn = nn.reshape((-1, 1), order="F")
-        z = np.zeros_like(mm)
-        y = np.hstack([z, mm, nn])
-        return y
+        m = np.arange(-M, M + 1)
+        n = np.arange(-N, N + 1)
+        m, n = np.meshgrid(m, n, indexing="ij")
+        m = m.ravel()
+        n = n.ravel()
+        z = np.zeros_like(m)
+        return np.array([z, m, n]).T
 
     def evaluate(
         self, nodes, derivatives=np.array([0, 0, 0]), modes=None, unique=False
@@ -531,6 +528,11 @@ def evaluate(
         -------
         y : ndarray, shape(num_nodes,num_modes)
             Basis functions evaluated at nodes.
+            The Vandermonde matrix when ``modes is None`` is
+            given by ``y.reshape(-1,2*M+1,2*N+1)`` and
+            is an outer product of Fourier matrices with order
+            [sin(M𝛉), ..., sin(𝛉), 1, cos(𝛉), ..., cos(M𝛉)]
+            ⊗ [sin(N𝛇), ..., sin(𝛇), 1, cos(𝛇), ..., cos(N𝛇)].
 
         """
         if modes is None:
@@ -858,17 +860,14 @@ def _get_modes(self, L, M, N):
             Each row is one basis function with modes (l,m,n).
 
         """
-        dim_pol = 2 * M + 1
-        dim_tor = 2 * N + 1
         l = np.arange(L + 1)
-        m = np.arange(dim_pol) - M
-        n = np.arange(dim_tor) - N
-        ll, mm, nn = np.meshgrid(l, m, n)
-        ll = ll.reshape((-1, 1), order="F")
-        mm = mm.reshape((-1, 1), order="F")
-        nn = nn.reshape((-1, 1), order="F")
-        y = np.hstack([ll, mm, nn])
-        return y
+        m = np.arange(-M, M + 1)
+        n = np.arange(-N, N + 1)
+        l, m, n = np.meshgrid(l, m, n, indexing="ij")
+        l = l.ravel()
+        m = m.ravel()
+        n = n.ravel()
+        return np.array([l, m, n]).T
 
     def evaluate(
         self, nodes, derivatives=np.array([0, 0, 0]), modes=None, unique=False
@@ -891,6 +890,12 @@ def evaluate(
         -------
         y : ndarray, shape(num_nodes,num_modes)
             Basis functions evaluated at nodes.
+            The Vandermonde matrix when ``modes is None`` is given by
+            ``y.reshape(-1,L+1,2*M+1,2*N+1,3)`` and is
+            an outer product of Chebyshev and Fourier matrices with order
+            [T₀(𝛒), T₁(𝛒), ..., T_L(𝛒)]
+            ⊗ [sin(M𝛉), ..., sin(𝛉), 1, cos(𝛉), ..., cos(M𝛉)]
+            ⊗ [sin(N𝛇), ..., sin(𝛇), 1, cos(𝛇), ..., cos(N𝛇)].
 
         """
         if modes is None:
@@ -1212,9 +1217,9 @@ def _get_modes(self, L):
             Each row is one basis function with modes (l,m,n).
 
         """
-        l = np.arange(L + 1).reshape((-1, 1))
-        z = np.zeros((L + 1, 2))
-        return np.hstack([l, z])
+        l = np.arange(L + 1)
+        z = np.zeros_like(l)
+        return np.array([l, z, z]).T
 
     def evaluate(
         self, nodes, derivatives=np.array([0, 0, 0]), modes=None, unique=False
@@ -1627,7 +1632,7 @@ def chebyshev(r, l, dr=0):
 
     Parameters
     ----------
-    rho : ndarray, shape(N,)
+    r : ndarray, shape(N,)
         radial coordinates to evaluate basis
     l : ndarray of int, shape(K,)
         radial mode number(s)

desc/compute/_deprecated.py

@@ -351,13 +351,12 @@ def _Gamma_c_Velasco_1D(params, transforms, profiles, data, **kwargs):
 
     def Gamma_c(data):
         bounce = Bounce1D(grid, data, quad, is_reshaped=True)
-        points = bounce.points(data["pitch_inv"], num_well)
         v_tau, radial_drift, poloidal_drift = bounce.integrate(
             [_v_tau, _radial_drift, _poloidal_drift],
             data["pitch_inv"],
             data,
             ["cvdrift0", "gbdrift"],
-            points,
+            bounce.points(data["pitch_inv"], num_well),
         )
         # This is γ_c π/2.
         gamma_c = jnp.arctan(safediv(radial_drift, poloidal_drift))

desc/compute/_equil.py

@@ -824,12 +824,12 @@ def _beta_voltor(params, transforms, profiles, data, **kwargs):
 def _P_ISS04(params, transforms, profiles, data, **kwargs):
     rho = transforms["grid"].compress(data["rho"], surface_label="rho")
     iota = transforms["grid"].compress(data["iota"], surface_label="rho")
+    method = kwargs.get("method", "cubic")
     fx = {}
-    if "iota_r" in data:
-        fx["fx"] = transforms["grid"].compress(
-            data["iota_r"]
-        )  # noqa: unused dependency
-    iota_23 = interp1d(2 / 3, rho, iota, method=kwargs.get("method", "cubic"), **fx)
+    if "iota_r" in data and method == "cubic":
+        # noqa: unused dependency
+        fx["fx"] = transforms["grid"].compress(data["iota_r"])
+    iota_23 = interp1d(2 / 3, rho, iota, method=method, **fx)
     data["P_ISS04"] = 1e6 * (  # MW -> W
         jnp.abs(data["W_p"] / 1e6)  # J -> MJ
         / (

desc/compute/geom_utils.py

@@ -68,7 +68,7 @@ def xyz2rpz(pts):
 
     """
     x, y, z = pts.T
-    r = jnp.sqrt(x**2 + y**2)
+    r = jnp.hypot(x, y)
     p = jnp.arctan2(y, x)
     return jnp.array([r, p, z]).T
 

desc/grid.py

@@ -990,6 +990,7 @@ def __init__(
         self._node_pattern = "linear"
         self._coordinates = "rtz"
         self._is_meshgrid = True
+        self._can_fft2 = not sym and not endpoint
         self._period = (np.inf, 2 * np.pi, 2 * np.pi / self._NFP)
         self._nodes, self._spacing = self._create_nodes(
             L=L,
@@ -1205,6 +1206,11 @@ def _create_nodes(  # noqa: C901
         self._endpoint = (self._poloidal_endpoint or (t.size == 1 and z.size > 1)) and (
             self._toroidal_endpoint or (z.size == 1 and t.size > 1)
         )
+        self._can_fft2 = (
+            self._can_fft2
+            and not self._poloidal_endpoint
+            and not self._toroidal_endpoint
+        )
 
         r, t, z = map(np.ravel, np.meshgrid(r, t, z, indexing="ij"))
         dr, dt, dz = map(np.ravel, np.meshgrid(dr, dt, dz, indexing="ij"))