Merge branch 'fieldline_compute' into integrate_on_boundary

CHANGELOG.md

@@ -3,12 +3,18 @@ Changelog
 
 New Features
 
+- All vector variables are now computed in toroidal (R,phi,Z) coordinates by default.
+  Cartesian (X,Y,Z) coordinates can be requested with the compute keyword ``basis='xyz'``.
 - Add method ``from_values`` to ``FourierRZCurve`` to allow fitting of data points
 to a ``FourierRZCurve`` object, and ``to_FourierRZCurve`` methods to ``Curve`` class.
 - Adds the objective `CoilsetMinDistance`, which returns the minimum distance to another
 coil for each coil in a coilset.
 - Adds the objective `PlasmaCoilsetMinDistance`, which returns the minimum distance to the
 plasma surface for each coil in a coilset.
+- Add method ``is_self_intersecting`` to ``CoilSet``, which checks if any coils intersect eachother in the coilset.
+- Removes error in ``from_symmetry`` method of ``CoilSet`` when a coil crosses the symmetry plane,
+and instead adds a check for intersection, to allow for valid coilsets which may cross the
+symmetry plane but not be self-intersecting after rotation/reflection.
 
 v0.11.1
 -------

desc/coils.py

@@ -19,6 +19,7 @@
 from desc.compute import get_params, rpz2xyz, rpz2xyz_vec, xyz2rpz, xyz2rpz_vec
 from desc.compute.geom_utils import reflection_matrix
 from desc.compute.utils import _compute as compute_fun
+from desc.compute.utils import safenorm
 from desc.geometry import (
     FourierPlanarCurve,
     FourierRZCurve,
@@ -797,13 +798,22 @@ class CoilSet(OptimizableCollection, _Coil, MutableSequence):
         assuming 'virtual' coils from the other half field period. Default = False.
     name : str
         Name of this CoilSet.
+    check_intersection: bool
+        Whether or not to check the coils in the coilset for intersections.
 
     """
 
     _io_attrs_ = _Coil._io_attrs_ + ["_coils", "_NFP", "_sym"]
     _io_attrs_.remove("_current")
 
-    def __init__(self, *coils, NFP=1, sym=False, name=""):
+    def __init__(
+        self,
+        *coils,
+        NFP=1,
+        sym=False,
+        name="",
+        check_intersection=True,
+    ):
         coils = flatten_list(coils, flatten_tuple=True)
         assert all([isinstance(coil, (_Coil)) for coil in coils])
         [_check_type(coil, coils[0]) for coil in coils]
@@ -812,6 +822,9 @@ def __init__(self, *coils, NFP=1, sym=False, name=""):
         self._sym = bool(sym)
         self._name = str(name)
 
+        if check_intersection:
+            self.is_self_intersecting()
+
     @property
     def name(self):
         """str: Name of the curve."""
@@ -899,7 +912,10 @@ def compute(
 
         """
         if params is None:
-            params = [get_params(names, coil) for coil in self]
+            params = [
+                get_params(names, coil, basis=kwargs.get("basis", "rpz"))
+                for coil in self
+            ]
         if data is None:
             data = [{}] * len(self)
 
@@ -940,9 +956,9 @@ def _compute_position(self, params=None, grid=None, **kwargs):
             Coil positions, in [R,phi,Z] or [X,Y,Z] coordinates.
 
         """
-        if params is None:
-            params = [get_params("x", coil) for coil in self]
         basis = kwargs.pop("basis", "xyz")
+        if params is None:
+            params = [get_params("x", coil, basis=basis) for coil in self]
         data = self.compute("x", grid=grid, params=params, basis=basis, **kwargs)
         data = tree_leaves(data, is_leaf=lambda x: isinstance(x, dict))
         x = jnp.dstack([d["x"].T for d in data]).T  # shape=(ncoils,num_nodes,3)
@@ -1009,7 +1025,9 @@ def compute_magnetic_field(
         assert basis.lower() in ["rpz", "xyz"]
         coords = jnp.atleast_2d(jnp.asarray(coords))
         if params is None:
-            params = [get_params(["x_s", "x", "s", "ds"], coil) for coil in self]
+            params = [
+                get_params(["x_s", "x", "s", "ds"], coil, basis=basis) for coil in self
+            ]
             for par, coil in zip(params, self):
                 par["current"] = coil.current
 
@@ -1169,21 +1187,6 @@ def from_symmetry(cls, coils, NFP=1, sym=False):
             # only need to check this for a CoilSet, not MixedCoilSet
             [_check_type(coil, coils[0]) for coil in coils]
 
-        # check toroidal extent of coils to be repeated
-        maxphi = 2 * np.pi / NFP / (sym + 1)
-        data = coils.compute("phi")
-        for i, cdata in enumerate(data):
-            errorif(
-                np.any(cdata["phi"] > maxphi),
-                ValueError,
-                f"coil {i} exceeds the toroidal extent for NFP={NFP} and sym={sym}",
-            )
-            warnif(
-                sym and np.any(cdata["phi"] < np.finfo(cdata["phi"].dtype).eps),
-                UserWarning,
-                f"coil {i} is on the symmetry plane phi=0",
-            )
-
         coilset = []
         if sym:
             # first reflect/flip original coilset
@@ -1498,6 +1501,100 @@ def to_SplineXYZ(self, knots=None, grid=None, method="cubic", name=""):
         coils = [coil.to_SplineXYZ(knots, grid, method) for coil in self]
         return self.__class__(*coils, NFP=self.NFP, sym=self.sym, name=name)
 
+    def is_self_intersecting(self, grid=None, tol=None):
+        """Check if any coils in the CoilSet intersect.
+
+        By default, checks intersection by checking that for each point on a given coil
+        the closest point in the coilset is on that same coil. If the closest point is
+        on another coil, that indicates that the coils may be close to intersecting.
+
+        If instead the ``tol`` argument is provided, then the function will
+        check the minimum distance from each coil to each other coil against
+        that tol and if it finds the minimum distance is less than the ``tol``,
+        it will take it as intersecting coils, returning True and raising a warning.
+
+        NOTE: If grid resolution used is too low, this function may fail to return
+        the correct answer.
+
+        Parameters
+        ----------
+        grid : Grid, optional
+            Collocation grid containing the nodes to evaluate the coil positions at.
+            If a list, must have the same structure as the coilset. Defaults to a
+            LinearGrid(N=100)
+        tol : float, optional
+            the tolerance (in meters) to check the intersections to, if points on any
+            two coils are closer than this tolerance, then the function will return
+            True and a warning will be raised. If not passed, then the method used
+            to determine coilset intersection will be based off of checking that
+            each point on a coil is closest to a point on the same coil, which does
+            not rely on a ``tol`` parameter.
+
+        Returns
+        -------
+        is_self_intersecting : bool
+            Whether or not any coils in the CoilSet come close enough to each other to
+            possibly be intersecting.
+
+        """
+        from desc.objectives._coils import CoilsetMinDistance
+
+        grid = grid if grid else LinearGrid(N=100)
+        obj = CoilsetMinDistance(self, grid=grid)
+        obj.build(verbose=0)
+        if tol:
+            min_dists = obj.compute(self.params_dict)
+            is_nearly_intersecting = np.any(min_dists < tol)
+            warnif(
+                is_nearly_intersecting,
+                UserWarning,
+                "Found coils which are nearly intersecting according to the given tol "
+                + "(min coil-coil distance = "
+                + f"{np.min(min_dists):1.3e} m < {tol:1.3e} m)"
+                + " in the coilset, it is recommended to check coils closely.",
+            )
+            return is_nearly_intersecting
+        else:
+
+            pts = obj._constants["coilset"]._compute_position(
+                params=self.params_dict, grid=obj._constants["grid"], basis="xyz"
+            )
+            pts = np.array(pts)
+            num_nodes = pts.shape[1]
+            bad_coil_inds = []
+            # We will raise the warning if the jth point on the
+            # kth coil is closer to a point on a different coil than
+            # it is to the neighboring points on itself
+            for k in range(self.num_coils):
+                # dist[i,j,n] is the distance from the jth point on the kth coil
+                # to the nth point on the ith coil
+                dist = np.asarray(
+                    safenorm(pts[k][None, :, None] - pts[:, None, :], axis=-1)
+                )
+                for j in range(num_nodes):
+                    dists_for_this_pt = dist[:, j, :].copy()
+                    dists_for_this_pt[k][
+                        j
+                    ] = np.inf  # Set the dist from the pt to itself to inf to ignore
+                    ind_min = np.argmin(dists_for_this_pt)
+                    # check if the index returned corresponds to a point on the same
+                    # coil. if it does not, then this jth pt on the kth coil is closer
+                    #  to a point on another coil than it is to pts on its own coil,
+                    # which means it may be intersecting it.
+                    if ind_min not in np.arange((num_nodes) * k, (num_nodes) * (k + 1)):
+                        bad_coil_inds.append(k)
+            bad_coil_inds = set(bad_coil_inds)
+            is_nearly_intersecting = True if bad_coil_inds else False
+            warnif(
+                is_nearly_intersecting,
+                UserWarning,
+                "Found coils which are nearly intersecting according to the given grid"
+                + " it is recommended to check coils closely or run function "
+                + "again with a higher resolution grid."
+                + f" Offending coil indices are {bad_coil_inds}.",
+            )
+            return is_nearly_intersecting
+
     def __add__(self, other):
         if isinstance(other, (CoilSet)):
             return CoilSet(*self.coils, *other.coils)
@@ -1548,18 +1645,22 @@ class MixedCoilSet(CoilSet):
         Collection of coils.
     name : str
         Name of this CoilSet.
+    check_intersection: bool
+        Whether or not to check the coils in the coilset for intersections.
 
     """
 
     _io_attrs_ = CoilSet._io_attrs_
 
-    def __init__(self, *coils, name=""):
+    def __init__(self, *coils, name="", check_intersection=True):
         coils = flatten_list(coils, flatten_tuple=True)
         assert all([isinstance(coil, (_Coil)) for coil in coils])
         self._coils = list(coils)
         self._NFP = 1
         self._sym = False
         self._name = str(name)
+        if check_intersection:
+            self.is_self_intersecting()
 
     @property
     def num_coils(self):

desc/compute/__init__.py

@@ -63,10 +63,10 @@ def _build_data_index():
     for p in data_index:
         for key in data_index[p]:
             full = {
-                "data": get_data_deps(key, p, has_axis=False),
-                "transforms": get_derivs(key, p, has_axis=False),
-                "params": get_params(key, p, has_axis=False),
-                "profiles": get_profiles(key, p, has_axis=False),
+                "data": get_data_deps(key, p, has_axis=False, basis="rpz"),
+                "transforms": get_derivs(key, p, has_axis=False, basis="rpz"),
+                "params": get_params(key, p, has_axis=False, basis="rpz"),
+                "profiles": get_profiles(key, p, has_axis=False, basis="rpz"),
             }
             data_index[p][key]["full_dependencies"] = full
 
@@ -81,9 +81,9 @@ def _build_data_index():
             else:
                 full_with_axis = {
                     "data": full_with_axis_data,
-                    "transforms": get_derivs(key, p, has_axis=True),
-                    "params": get_params(key, p, has_axis=True),
-                    "profiles": get_profiles(key, p, has_axis=True),
+                    "transforms": get_derivs(key, p, has_axis=True, basis="rpz"),
+                    "params": get_params(key, p, has_axis=True, basis="rpz"),
+                    "profiles": get_profiles(key, p, has_axis=True, basis="rpz"),
                 }
                 for _key, val in full_with_axis.items():
                     if full[_key] == val: