Merge branch 'rc/coils_init' of github.com:PlasmaControl/DESC into rc…

…/coils_init

.github/workflows/changelog_update.yml

@@ -0,0 +1,35 @@
+name: Check changelog updated
+
+on:
+  pull_request:
+    branches:
+      - master
+    types: [opened, synchronize, labeled, unlabeled]
+
+jobs:
+  check_changelog_updated:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Filter changes
+        id: changes
+        uses: dorny/paths-filter@v3
+        with:
+          filters: |
+            has_changes:
+              - 'desc/**'
+              - 'requirements.txt'
+              - 'requirements_conda.yml'
+              - '.github/workflows/changelog_update.yml'
+
+      - name: Check for relevant changes
+        id: check_changes
+        run: echo "has_changes=${{ steps.changes.outputs.has_changes }}" >> $GITHUB_ENV
+
+      - uses: actions/checkout@v4
+
+      - if: ${{ !contains(github.event.pull_request.labels.*.name, 'skip_changelog') && env.has_changes == 'true'}}
+        uses: danieljimeneznz/[email protected]
+        with:
+          require-changes-to: |
+            CHANGELOG.md
+          token: ${{ secrets.GITHUB_TOKEN }}

CHANGELOG.md

@@ -2,7 +2,8 @@ Changelog
 =========
 
 New Features
-- Add ``from_input_file`` method to ``Equilibrium`` class to generate an ``Equilibrium`` object with boundary, profiles, resolution and flux specified in a given DESC or VMEC input file
+- Add ``from_input_file`` method to ``Equilibrium`` class to generate an ``Equilibrium`` object with boundary, profiles, resolution and flux specified in a given DESC or VMEC input file.
+- Add ``desc.coils.initialize_modular_coils`` and ``desc.coils.initialize_saddle_coils`` for creating an initial guess for stage 2 optimization.
 
 
 Bug Fixes

desc/coils.py

@@ -2039,7 +2039,7 @@ def to_FourierPlanar(
     def to_FourierRZ(
         self, N=10, grid=None, NFP=None, sym=False, name="", check_intersection=True
     ):
-        """Convert all coils to FourierRZCoil representaion.
+        """Convert all coils to FourierRZCoil representation.
 
         Note that some types of coils may not be representable in this basis.
 
@@ -2058,7 +2058,7 @@ def to_FourierRZ(
         name : str
             Name for this coilset.
         check_intersection: bool
-            Whether or not to check the coils in the new coiilset for intersections.
+            Whether or not to check the coils in the new coilset for intersections.
 
         Returns
         -------
@@ -2091,7 +2091,7 @@ def to_FourierXYZ(self, N=10, grid=None, s=None, name="", check_intersection=Tru
         name : str
             Name for the new CoilSet.
         check_intersection: bool
-            Whether or not to check the coils in the new coiilset for intersections.
+            Whether or not to check the coils in the new coilset for intersections.
 
         Returns
         -------
@@ -2136,7 +2136,7 @@ def to_SplineXYZ(
         name : str
             Name for the new CoilSet.
         check_intersection: bool
-            Whether or not to check the coils in the new coiilset for intersections.
+            Whether or not to check the coils in the new coilset for intersections.
 
         Returns
         -------
@@ -2549,7 +2549,7 @@ def to_FourierPlanar(
         name : str
             Name for the new MixedCoilSet.
         check_intersection: bool
-            Whether or not to check the coils in the new coiilset for intersections.
+            Whether or not to check the coils in the new coilset for intersections.
 
         Returns
         -------
@@ -2588,7 +2588,7 @@ def to_FourierRZ(
         name : str
             Name for the new MixedCoilSet.
         check_intersection: bool
-            Whether or not to check the coils in the new coiilset for intersections.
+            Whether or not to check the coils in the new coilset for intersections.
 
         Returns
         -------
@@ -2620,7 +2620,7 @@ def to_FourierXYZ(self, N=10, grid=None, s=None, name="", check_intersection=Tru
         name : str
             Name for the new MixedCoilSet.
         check_intersection: bool
-            Whether or not to check the coils in the new coiilset for intersections.
+            Whether or not to check the coils in the new coilset for intersections.
 
         Returns
         -------
@@ -2662,7 +2662,7 @@ def to_SplineXYZ(
         name : str
             Name for the new MixedCoilSet.
         check_intersection: bool
-            Whether or not to check the coils in the new coiilset for intersections.
+            Whether or not to check the coils in the new coilset for intersections.
 
         Returns
         -------
@@ -2881,15 +2881,15 @@ def _linking_number(x1, x2, x1_s, x2_s, dx1, dx2):
     return ratio.sum()
 
 
-def initialize_modular_coils(eq, num_coils, r_over_a=1.5):
-    """Initilize a CoilSet of modular coils for stage 2 optimization.
+def initialize_modular_coils(eq, num_coils, r_over_a=2.0):
+    """Initialize a CoilSet of modular coils for stage 2 optimization.
 
-    The coils will be planar, circular coils centered on the equilibrium magnetic axis, and aligned such that
-    the normal to the coil points along the axis. The currents will be set to match
-    the equilibrium required poloidal linking current.
+    The coils will be planar, circular coils centered on the equilibrium magnetic axis,
+    and aligned such that the normal to the coil points along the axis. The currents
+    will be set to match the equilibrium required poloidal linking current.
 
-    The coils will be FourierPlanarCoils, if another type is desired use
-    ``coilset.to_FourierXYZ()``, ``coilset.to_SplineXYZ()`` etc.
+    The coils will be ``FourierPlanarCoil`` with N=0, if another type is desired use
+    ``coilset.to_FourierXYZ(N=10)``, ``coilset.to_SplineXYZ()`` etc.
 
     Parameters
     ----------
@@ -2899,7 +2899,9 @@ def initialize_modular_coils(eq, num_coils, r_over_a=1.5):
         Number of coils to create per field period. For stellarator symmetric
         equilibria, this will be the number of coils per half-period.
     r_over_a : float
-        Minor radius of the coils, in units of equilibrium minor radius.
+        Minor radius of the coils, in units of equilibrium minor radius. Note that for
+        strongly shaped equilibria this may need to be large to avoid having the coils
+        intersect the plasma.
 
     Returns
     -------
@@ -2932,14 +2934,14 @@ def initialize_modular_coils(eq, num_coils, r_over_a=1.5):
 
 
 def initialize_saddle_coils(eq, num_coils, r_over_a=0.5, offset=2.0, position="outer"):
-    """Initilize a CoilSet of saddle coils for stage 2 optimization.
+    """Initialize a CoilSet of saddle coils for stage 2 optimization.
 
-    The coils will be positioned around the plasma without linking it, and aligned
-    such that the normal to the coil points towards the magnetic axis. The currents
-    will be initialized to zero.
+    The coils will be planar, circular coils positioned around the plasma without
+    linking it, and aligned such that the normal to the coil points towards the
+    magnetic axis. The currents will be initialized to zero.
 
-    The coils will be FourierPlanarCoils, if another type is desired use
-    ``coilset.to_FourierXYZ()``, ``coilset.to_SplineXYZ()`` etc.
+    The coils will be ``FourierPlanarCoil`` with N=0, if another type is desired use
+    ``coilset.to_FourierXYZ(N=10)``, ``coilset.to_SplineXYZ()`` etc.
 
     Parameters
     ----------
@@ -2952,6 +2954,8 @@ def initialize_saddle_coils(eq, num_coils, r_over_a=0.5, offset=2.0, position="o
         Minor radius of the coils, in units of equilibrium minor radius.
     offset : float
         Distance from coil to magnetic axis, in units of equilibrium minor radius.
+        Note that for strongly shaped equilibria this may need to be large to avoid
+        having the coils intersect the plasma.
     position : {"outer", "inner", "top", "bottom"}
         Placement of coils relative to plasma. "outer" will place coils on the outboard
         side, "inner" on the inboard side, "top" will place coils above the plasma,

tests/conftest.py

@@ -284,13 +284,13 @@ def DummyMixedCoilSet(tmpdir_factory):
         vf_coil, displacement=[0, 0, 2], n=3, endpoint=True
     )
     xyz_coil = FourierXYZCoil(current=2)
-    phi = 2 * np.pi * np.linspace(0, 1, 20, endpoint=True) ** 2
+    phi = 2 * np.pi * np.linspace(0, 1, 20, endpoint=True)
     spline_coil = SplineXYZCoil(
         current=1,
         X=np.cos(phi),
         Y=np.sin(phi),
         Z=np.zeros_like(phi),
-        knots=np.linspace(0, 2 * np.pi, len(phi)),
+        knots=phi,
     )
     full_coilset = MixedCoilSet(
         (tf_coilset, vf_coilset, xyz_coil, spline_coil), check_intersection=False

tests/test_examples.py

@@ -1324,7 +1324,20 @@ def test_second_stage_optimization_CoilSet():
 
 @pytest.mark.slow
 @pytest.mark.unit
-def test_optimize_with_all_coil_types(DummyCoilSet, DummyMixedCoilSet):
+@pytest.mark.parametrize(
+    "coil_type",
+    [
+        "FourierPlanarCoil",
+        "FourierRZCoil",
+        "FourierXYZCoil",
+        "SplineXYZCoil",
+        "CoilSet sym",
+        "CoilSet asym",
+        "MixedCoilSet",
+        "nested CoilSet",
+    ],
+)
+def test_optimize_with_all_coil_types(DummyCoilSet, DummyMixedCoilSet, coil_type):
     """Test optimizing for every type of coil and dummy coil sets."""
     sym_coils = load(load_from=str(DummyCoilSet["output_path_sym"]), file_format="hdf5")
     asym_coils = load(
@@ -1340,66 +1353,63 @@ def test_optimize_with_all_coil_types(DummyCoilSet, DummyMixedCoilSet):
     quad_eval_grid = LinearGrid(M=2, sym=True)
     quad_field_grid = LinearGrid(N=2)
 
-    def test(c, method):
-        target = 11
-        rtol = 1e-3
-        # first just check that quad flux works for a couple iterations
-        # as this is an expensive objective to compute
-        obj = ObjectiveFunction(
-            QuadraticFlux(
-                eq=eq,
-                field=c,
-                vacuum=True,
-                weight=1e-4,
-                eval_grid=quad_eval_grid,
-                field_grid=quad_field_grid,
-            )
-        )
-        optimizer = Optimizer(method)
-        (c,), _ = optimizer.optimize(c, obj, maxiter=2, ftol=0, xtol=1e-15)
-
-        # now check with optimizing geometry and actually check result
-        objs = [
-            CoilLength(c, target=target),
-        ]
-        extra_msg = ""
-        if isinstance(c, MixedCoilSet):
-            # just to check they work without error
-            objs.extend(
-                [
-                    CoilCurvature(c, target=0.5, weight=1e-2),
-                    CoilTorsion(c, target=0, weight=1e-2),
-                ]
-            )
-            rtol = 3e-2
-            extra_msg = " with curvature and torsion obj"
-
-        obj = ObjectiveFunction(objs)
-
-        (c,), _ = optimizer.optimize(c, obj, maxiter=25, ftol=5e-3, xtol=1e-15)
-        flattened_coils = tree_leaves(
-            c, is_leaf=lambda x: isinstance(x, _Coil) and not isinstance(x, CoilSet)
-        )
-        lengths = [coil.compute("length")["length"] for coil in flattened_coils]
-        np.testing.assert_allclose(
-            lengths, target, rtol=rtol, err_msg=f"lengths {c}" + extra_msg
-        )
-
     spline_coil = mixed_coils.coils[-1].copy()
 
-    # single coil
-    test(FourierPlanarCoil(), "fmintr")
-    test(FourierRZCoil(), "fmintr")
-    test(FourierXYZCoil(), "fmintr")
-    test(spline_coil, "fmintr")
+    types = {
+        "FourierPlanarCoil": (FourierPlanarCoil(), "fmintr"),
+        "FourierRZCoil": (FourierRZCoil(), "fmintr"),
+        "FourierXYZCoil": (FourierXYZCoil(), "fmintr"),
+        "SplineXYZCoil": (spline_coil, "fmintr"),
+        "CoilSet sym": (sym_coils, "lsq-exact"),
+        "CoilSet asym": (asym_coils, "lsq-exact"),
+        "MixedCoilSet": (mixed_coils, "lsq-exact"),
+        "nested CoilSet": (nested_coils, "lsq-exact"),
+    }
+    c, method = types[coil_type]
+
+    target = 11
+    rtol = 1e-3
+    # first just check that quad flux works for a couple iterations
+    # as this is an expensive objective to compute
+    obj = ObjectiveFunction(
+        QuadraticFlux(
+            eq=eq,
+            field=c,
+            vacuum=True,
+            weight=1e-4,
+            eval_grid=quad_eval_grid,
+            field_grid=quad_field_grid,
+        )
+    )
+    optimizer = Optimizer(method)
+    (cc,), _ = optimizer.optimize(c, obj, maxiter=2, ftol=0, xtol=1e-8, copy=True)
+
+    # now check with optimizing geometry and actually check result
+    objs = [
+        CoilLength(c, target=target),
+    ]
+    extra_msg = ""
+    if isinstance(c, MixedCoilSet):
+        # just to check they work without error
+        objs.extend(
+            [
+                CoilCurvature(c, target=0.5, weight=1e-2),
+                CoilTorsion(c, target=0, weight=1e-2),
+            ]
+        )
+        rtol = 3e-2
+        extra_msg = " with curvature and torsion obj"
 
-    # CoilSet
-    test(sym_coils, "lsq-exact")
-    test(asym_coils, "lsq-exact")
+    obj = ObjectiveFunction(objs)
 
-    # MixedCoilSet
-    test(mixed_coils, "lsq-exact")
-    test(nested_coils, "lsq-exact")
+    (c,), _ = optimizer.optimize(c, obj, maxiter=25, ftol=5e-3, xtol=1e-8)
+    flattened_coils = tree_leaves(
+        c, is_leaf=lambda x: isinstance(x, _Coil) and not isinstance(x, CoilSet)
+    )
+    lengths = [coil.compute("length")["length"] for coil in flattened_coils]
+    np.testing.assert_allclose(
+        lengths, target, rtol=rtol, err_msg=f"lengths {c}" + extra_msg
+    )
 
 
 @pytest.mark.unit