Change to CoilSetLinkingNumber

CHANGELOG.md

@@ -10,6 +10,7 @@ Performance Improvements
 Bug Fixes
 
 - Fixes bug in ``auglag`` optimizers which prevented them from accepting solver hyperparameters.
+- Fixes computation of ``CoilSetLinkingNumber`` to exclude coil writhe.
 - Fixes bug in modified Cholesky factorization used by the trust-region
   subproblems when the Gershgorin lower bound of the Hessian was exactly zero
   (e.g. a Hessian with an all-zero row), producing NaN steps in ``fmintr`` and

desc/coils.py

@@ -1888,6 +1888,8 @@ def _compute_linking_number(self, params=None, grid=None):
         link : ndarray, shape(num_coils, num_coils)
             Linking number of each coil with each other coil. link=0 means they are not
             linked, +/- 1 means the coils link each other in one direction or another.
+            Diagonal entries represent the writhe of a coil, and can be non-zero for
+            non-planar coils.
 
         """
         if grid is None:

desc/objectives/_coils.py

@@ -2634,7 +2634,9 @@ def compute(self, params, constants=None):
             params=params, grid=constants["grid"]
         )
 
-        return jnp.abs(link).sum(axis=0)
+        # the diagonal entries of "link" should be excluded
+        mask = ~jnp.eye(self._dim_f, dtype=bool)
+        return jnp.abs(link).sum(axis=0, where=mask)
 
 
 class SurfaceCurrentRegularization(_Objective):

tests/test_coils.py

@@ -1463,9 +1463,12 @@ def test_linking_number():
     link = coilset2._compute_linking_number(grid=grid)
 
     # modular coils don't link each other
+    # note we remove diagonal entries explicitly, as they do not
+    # represent linking numbers (they are zero here anyway).
+    mask = ~jnp.eye(link.shape[0], dtype=bool)
+    link = jnp.where(mask, link, 0)
     np.testing.assert_allclose(link[:-1, :-1], 0, atol=1e-14)
-    # axis coil doesn't link itself
-    np.testing.assert_allclose(link[-1, -1], 0, atol=1e-14)
+
     # we expect the axis coil to link all the modular coils, with alternating sign
     # due to alternating orientation of the coils due to symmetry.
     expected = [1, -1] * 5

tests/test_objective_funs.py

@@ -1758,6 +1758,7 @@ def test_coil_linking_number(self):
         # one way and half going the other way
         coilset = CoilSet.from_symmetry(coil, NFP=5, sym=True, check_intersection=False)
         coil2 = FourierRZCoil()
+
         # add a coil along the axis that links all the other coils
         coilset2 = MixedCoilSet(coilset, coil2, check_intersection=False)
 
@@ -1769,6 +1770,17 @@ def test_coil_linking_number(self):
         expected = np.array([1] * 10 + [10])
         np.testing.assert_allclose(out, expected, rtol=1e-3)
 
+        # produce a coilset of non-planar coils with nonzero writhe
+        coil3 = FourierXYZCoil(
+            X_n=[0, 0.5, 5, 1, 0], Y_n=[0, 0, 0, 0, 0.5], Z_n=[0, -1, 0, 0, 1]
+        )
+        coilset3 = CoilSet.from_symmetry(coil3, NFP=5, check_intersection=False)
+        obj = CoilSetLinkingNumber(coilset3)
+        obj.build()
+        out = obj.compute_scaled_error(coilset3.params_dict)
+        expected = np.array([0] * 5)
+        np.testing.assert_allclose(out, expected, atol=1e-12)
+
     @pytest.mark.unit
     def test_signed_plasma_vessel_distance(self):
         """Test calculation of signed distance from plasma to vessel."""