Merging main

.github/workflows/publish.yml

@@ -0,0 +1,48 @@
+name: Publish to PyPI
+
+on:
+  push:
+    tags:
+      - "v*"
+
+jobs:
+  build-and-publish:
+    runs-on: ubuntu-latest
+    permissions:
+      contents: read
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Set up Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: "3.11"
+
+      - name: Install build tools
+        run: python -m pip install --upgrade build twine
+
+      - name: Clean old artifacts
+        run: rm -rf dist build src/*.egg-info
+
+      - name: Build distributions
+        run: python -m build
+
+      - name: Validate distributions
+        run: python -m twine check dist/*
+
+      - name: Upload to TestPyPI
+        env:
+          TWINE_USERNAME: __token__
+          TWINE_PASSWORD: ${{ secrets.TEST_PYPI_API_TOKEN }}
+        run: |
+          python -m twine upload \
+            --repository-url https://test.pypi.org/legacy/ \
+            --skip-existing \
+            dist/*
+
+      - name: Upload to PyPI
+        env:
+          TWINE_USERNAME: __token__
+          TWINE_PASSWORD: ${{ secrets.PYPI_API_TOKEN }}
+        run: python -m twine upload dist/*

.gitignore

@@ -1,11 +1,10 @@
 docs/_build/
+dist/
+build/
+*.egg-info/
 *.pyc
+__pycache__/
+.ipynb_checkpoints/
 *checkpoint.py
 *checkpoint.ipynb
-src/quadcoil.egg-info/
-__pycache__/
-src/quadcoil/.ipynb_checkpoints/surfacerzfourier_jax-checkpoint.py
-.gitignore
-src/quadcoil/io/.ipynb_checkpoints/simsopt-checkpoint.py
-.gitignore
-src/quadcoil/io/.ipynb_checkpoints/__init__-checkpoint.py
+.claude/

README.md

@@ -8,8 +8,7 @@ QUADCOIL is a global coil optimization code that approximates coils with a smoot
 In other words, it's a "winding surface" code. However, unlike other winding surface codes, QUADCOIL:
 
 - Supports constrained optimization.
-- Supports non-convex quadratic penalties/constraints, such as curvature 
-  $\mathbf{K} \cdot \nabla \mathbf{K}$.
+- Supports non-convex quadratic penalties/constraints, such as curvatures and Lorentz force.
 - Includes robust winding surface generators that do not produce self-intersections.
 - Calculates derivatives with respect to plasma shape, winding surface shape, objective weights, and constraint thresholds.
 

docs/conf.py

@@ -35,7 +35,7 @@
     'matplotlib',
 ]
 templates_path = ['_templates']
-exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
+exclude_patterns = ['_build', '.ipynb_checkpoints', 'Thumbs.db', '.DS_Store']
 
 language = 'en'
 

docs/index.rst

@@ -61,6 +61,8 @@ Publications
    tutorial_inputs
    tutorial_outputs
    tutorial_to_desc
+   tutorial_misc_outputs
    tutorial_to_simsopt
    quantity    
    quadcoil
+   version_history

docs/quadcoil.io.rst

@@ -6,24 +6,16 @@ This package handles file input/output, coil cutting, and provides interfaces fo
 Submodules
 ----------
 
-quadcoil.io.file module
------------------------
-
-.. automodule:: quadcoil.io.file
-   :members:
-   :show-inheritance:
-   :undoc-members:
-
 quadcoil.io.coil\_cutting module
----------------------------------
+--------------------------------
 
 .. automodule:: quadcoil.io.coil_cutting
    :members:
    :show-inheritance:
    :undoc-members:
 
 quadcoil.io.focus module
--------------------------
+------------------------
 
 .. automodule:: quadcoil.io.focus
    :members:

docs/quantity.rst

@@ -12,7 +12,7 @@ simply pass their names into ``quadcoil.quadcoil``.
 
 .. code-block:: python
 
-    phi_mn, out_dict, qp, status = quadcoil(
+    out_dict, qp, dofs_opt, solve_results = quadcoil(
         ...
         objective_name=('f_B',),
         constraint_name=('K_theta',),
@@ -29,12 +29,12 @@ directly imported as functions from ``quadcoil.quantity``:
 .. code-block:: python
 
     from quadcoil.quantity import K_theta
-    print(K_theta(qp, phi_mn))
+    print(K_theta(qp, dofs_opt))
 
 All members of ``quadcoil.quantity`` require the same inputs:
 
 - ``qp : QuadcoilParams`` - Stores the plasma and winding surface information.
-- ``phi_mn : ndarray`` - The Fourier Coefficients of :math:`\Phi_{sv}` produced by ``quadcoil.quadcoil``.
+- ``dofs_opt : dict`` - The optimized degrees of freedom produced by ``quadcoil.quadcoil``.
 
 Notation
 ---------
@@ -68,6 +68,10 @@ These objectives are related to the magnetic field on the plasma surface:
      - :math:`\mathbf{B}\cdot\hat{\mathbf{n}}`
      - :math:`(n_\phi^P, n_\theta^P)`
      - The normal magnetic field on the plasma surface.
+   * - ``'Bnormal2'``
+     - :math:`(\mathbf{B}\cdot\hat{\mathbf{n}})^2`
+     - :math:`(n_\phi^P, n_\theta^P)`
+     - The squared normal field error on the plasma surface.
    * - ``'f_B'``
      - :math:`f_B\equiv\frac{n_{FP}}{2}\oint_\text{plasma} da \|\mathbf{B}\cdot\hat{\mathbf{n}}\|^2`
      - Scalar
@@ -76,14 +80,14 @@ These objectives are related to the magnetic field on the plasma surface:
      - :math:`\frac{f_B}{f_B(\Phi_{sv}=0)}`
      - Scalar
      - :math:`f_B`, normalized by its value with only the net toroidal and poloidal currents.
-   * - ``'f_max_Bnormal_abs'``
+   * - ``'f_max_Bnormal'``
      - :math:`\max_\text{plasma surface} \|\mathbf{B}\cdot\hat{\mathbf{n}}\|`
      - Scalar
-     - The maximum normal magnetic field strength.
+     - The maximum normal field error.
    * - ``'f_max_Bnormal2'``
      - :math:`\max_\text{plasma surface} \|\mathbf{B}\cdot\hat{\mathbf{n}}\|^2`
      - Scalar
-     - The maximum normal magnetic field strength squared. A convex quadratic constraint may behave better than a linear constraint.
+     - The maximum normal field error squared. A convex quadratic constraint may behave better than a linear constraint.
 
 Current Magnitude and Sign
 --------------------------
@@ -104,19 +108,23 @@ These objectives are related to the magnitude and sign of the sheet current :mat
    * - ``'K2'``
      - :math:`\|\mathbf{K}\|^2`
      - :math:`(n_\phi^E, n_\theta^E)`
-     - The current strength on the winding surface.
+     - The current density squared on the winding surface.
    * - ``'K_theta'``
      - :math:`K_\theta`
      - :math:`(n_\phi^E, n_\theta^E)`
      - The poloidal current distribution on the winding surface.
    * - ``'f_K'``
      - :math:`\frac{n_{FP}}{2}\oint_\text{WS} da \|\mathbf{K}\|^2`
      - Scalar
-     - The integrated magnetic field strength on the winding surface. Also the REGCOIL regularization factor.
+     - The integrated current density squared on the winding surface. Also the REGCOIL regularization factor.
    * - ``'f_max_K2'``
      - :math:`\max_\text{WS}\|K\|_2^2`
      - Scalar
-     - The integrated magnetic field strength on the winding surface. Also the REGCOIL regularization factor.
+     - The maximum current density squared on the winding surface.
+   * - ``'f_huber_K'``
+     - Pseudo Huber penalty on :math:`\mathbf{K}`
+     - Scalar
+     - A smooth approximation of the surface L-1 norm of :math:`|\mathbf{K}|`. Promotes sparsity in :math:`\mathbf{K}`.
 
 Current Curvature
 -----------------
@@ -135,13 +143,13 @@ These objectives are related to the curvature of the sheet current:
      - :math:`(n_\phi^E, n_\theta^E, 3)`
      - The :math:`(x, y, z)` components of :math:`\mathbf{K}\cdot\nabla\mathbf{K}` on the winding surface.
    * - ``'K_dot_grad_K_cyl'``
-     - :math:`(\mathbf{K}\cdot\nabla\mathbf{K})_{(R, \Phi, Z)}`
+     - :math:`(\mathbf{K}\cdot\nabla\mathbf{K})_{(R, \phi, Z)}`
      - :math:`(n_\phi^E, n_\theta^E, 3)`
-     - The :math:`(R, \Phi, Z)` components of :math:`\mathbf{K}\cdot\nabla\mathbf{K}` on the winding surface.
+     - The :math:`(R, \phi, Z)` components of :math:`\mathbf{K}\cdot\nabla\mathbf{K}` on the winding surface.
    * - ``'f_max_K_dot_grad_K_cyl'``
-     - :math:`\max_\text{WS}\|(\mathbf{K}\cdot\nabla\mathbf{K})_{(R, \Phi, Z)}\|_\infty`
+     - :math:`\max_\text{WS}\|(\mathbf{K}\cdot\nabla\mathbf{K})_{(R, \phi, Z)}\|_\infty`
      - Scalar
-     - Maximum :math:`(R, \Phi, Z)` component of :math:`\mathbf{K}\cdot\nabla\mathbf{K}` over the winding surface.
+     - Maximum :math:`(R, \phi, Z)` component of :math:`\mathbf{K}\cdot\nabla\mathbf{K}` over the winding surface.
 
 Dipole
 ------
@@ -159,14 +167,14 @@ These objectives are related to dipole optimization:
      - :math:`\Phi_{sv}`
      - :math:`(n_\phi^E, n_\theta^E)`
      - The dipole density distribution on the winding surface. Also referred to as the single valued component of the current potential.
-   * - ``'Phi_abs'``
-     - :math:`\|\Phi_{sv}\|`
-     - :math:`(n_\phi^E, n_\theta^E)`
-     - The absolute value of the dipole density distribution on the winding surface.
    * - ``'Phi2'``
      - :math:`\|\Phi_{sv}\|^2`
      - :math:`(n_\phi^E, n_\theta^E)`
      - The squared dipole density distribution on the winding surface.
+   * - ``'f_Phi'``
+     - :math:`\frac{n_{FP}}{2}\oint_\text{WS} da\ \|\Phi_{sv}\|^2`
+     - Scalar
+     - The integrated squared dipole density on the winding surface.
    * - ``'Phi_with_net_current'``
      - :math:`\Phi = \Phi_{sv} + \frac{G\phi'}{2\pi} + \frac{I\theta'}{2\pi}`
      - :math:`(n_\phi^E, n_\theta^E)`
@@ -183,8 +191,32 @@ These objectives are related to dipole optimization:
      - :math:`\max_\text{WS}\|\Phi_{sv}\|^2`
      - Scalar
      - The maximum dipole density squared on the winding surface. A convex quadratic constraint may behave better than a linear constraint.
+   * - ``'f_max_Phi4'``
+     - :math:`\max_\text{WS}\|\Phi_{sv}\|^4`
+     - Scalar
+     - The maximum fourth power of the dipole density on the winding surface. Experimental. Added to test the convergence behavior of high-order convex terms.
 
 Lorentz Force
 -------------
 
-Lorentz force is not yet fully implemented.
+These objectives are related to the self-force of the sheet current. The force is reported in :math:`(R, \phi, Z)` components on the winding surface.
+
+.. list-table::
+   :header-rows: 1
+
+   * - Name
+     - Formula
+     - Output Shape
+     - Description
+   * - ``'f_max_force_cyl'``
+     - :math:`\max_\text{WS}\|\mathbf{F}_{(R, \phi, Z)}\|_\infty`
+     - Scalar
+     - Maximum :math:`(R, \phi, Z)` component of the sheet-current self-force over the winding surface.
+   * - ``'f_l1_force_cyl'``
+     - :math:`\int_\text{WS} dA\ \|\mathbf{F}_{(R, \phi, Z)}\|_1`
+     - Scalar
+     - L1 penalty on the :math:`(R, \phi, Z)` components of the sheet-current self-force.
+   * - ``'f_max_force2_cyl'``
+     - :math:`\max_\text{WS}\|\mathbf{F}_{(R, \phi, Z)}\|_\infty^2`
+     - Scalar
+     - Maximum squared :math:`(R, \phi, Z)` component of the sheet-current self-force over the winding surface. Experimental. Added to test the convergence behavior of high-order non-convex terms.

docs/tutorial_inputs.rst

@@ -13,14 +13,14 @@ A minimal example can be found in ``examples/simple_example.ipynb``:
 .. code-block:: python
 
     from quadcoil import quadcoil
-    from simsopt import load
+    from simsopt.mhd import Vmec
 
     # Loading an equilibrium's boundary using simsopt
     equil_qs = Vmec('wout_LandremanPaul2021_QA_lowres.nc', keep_all_files=True)
     plasma_surface = equil_qs.boundary
     net_poloidal_current_amperes = equil_qs.external_current()
 
-    nescoil_out_dict, nescoil_qp, nescoil_phi_mn, _ = quadcoil(
+    nescoil_out_dict, nescoil_qp, nescoil_dofs, _ = quadcoil(
         nfp=plasma_surface.nfp,
         stellsym=plasma_surface.stellsym,
         mpol=4, # 4 poloidal harmonics for the current potential
@@ -34,32 +34,33 @@ A minimal example can be found in ``examples/simple_example.ipynb``:
         # Set the objective to 
         # f_B
         objective_name='f_B',
-        objective_weight=None,
+        objective_weight=1.,
         objective_unit=None,
         # Set the output metrics to f_B and f_K
-        metric_name=('f_B', 'f_K')
+        metric_name=('f_B', 'f_K'),
     )
 
     # Plotting the solution
-    from quadcoil.objective import Phi_with_net_current
+    from quadcoil.quantity import Phi_with_net_current
     import matplotlib.pyplot as plt
 
     plt.contour(
         nescoil_qp.quadpoints_phi, 
         nescoil_qp.quadpoints_theta,
-        Phi_with_net_current(nescoil_qp, nescoil_phi_mn), 
+        Phi_with_net_current(nescoil_qp, nescoil_dofs), 
         levels=40
     )
 
-Here, we solved the NESCOIL problem (minimizing field error with no additional constraints) on the Landreman-Paul QS configuration. This tutorial will explain how to set up a more compelex coil optimizer/proxy with QUADCOIL using ``quadcoil.quadcoil()``, by going over all input parameters and their physical meaning. These parameters fall in 7 categories:
+Here, we solved the NESCOIL problem (minimizing field error with no additional constraints) on the Landreman-Paul QS configuration. This tutorial will explain how to set up a more complex coil optimizer/proxy with QUADCOIL using ``quadcoil.quadcoil()``, by going over all input parameters and their physical meaning. These parameters fall in 8 categories:
 
 1. Plasma boundary
 2. Sheet current properties (net current, resolution, ...)
 3. Coil-plasma distance or winding surface
 4. Objective functions for coil optimization. Encodes engineering requirements.
 5. Constraints for coil optimization. Encodes engineering requirements.
-6. Metrics for evaluating the coil set satisfying these requirements.
-7. (Optional) Augmented Lagrangial options.
+6. Important numerical settings.
+7. Metrics for evaluating the coil set satisfying these requirements.
+8. (Optional) Augmented Lagrangial options.
 
 For readability, we label:
 
@@ -253,7 +254,7 @@ As we will see below, every objective and constraint **must be accompanied** by
 Single-objective
 ~~~~~~~~~~~~~~~~
 
-In this mode, QUADDCOIL will minimize one quantity selected from the list. To select single-objective mode, pass a single ``str`` as the ``objective_name``.
+In this mode, QUADCOIL will minimize one quantity selected from the list. To select single-objective mode, pass a single ``str`` as the ``objective_name``.
 
 .. list-table::
    :header-rows: 1
@@ -378,13 +379,14 @@ convert the non-smooth problem to a smooth problem. The currently supported valu
 
 .. list-table::
    :header-rows: 1
+
    * - Value for ``smoothing``
      - Type
      - Advantages
      - Disadvantages
    * - ``'slack'``
-     - Exact conversion using slasck variables.
-     - More accurate optimimum.
+     - Exact conversion using slack variables.
+     - More accurate optimum.
      - Inaccurate adjoint differentiation. Slower, higher memory usage, and high constraint count.
    * - ``'approx'``
      - Approximate conversion by replacing maximum with LogSumExp functions.
@@ -437,12 +439,8 @@ The augmented Lagrangian solver can be fine-tuned for a specific problem if the
      - The *c* factor. Please see *Constrained Optimization and Lagrange* *Multiplier Methods*, Chapter 3.
    * - ``c_growth_rate``
      - ``float``, traced
-     - ``1.2``
+     - ``2.``
      - The growth rate of the *c* factor.
-   * - ``fstop_outer``
-     - ``float``, traced
-     - ``1e-6``
-     - :math:`f_{obj}(\Phi_{sv})` stopping criterion of the outer augmented Lagrangian loop. Terminates the convergence rate falls below this number.
    * - ``xstop_outer``
      - ``float``, traced
      - ``1e-6``
@@ -471,5 +469,21 @@ The augmented Lagrangian solver can be fine-tuned for a specific problem if the
      - ``int``, static
      - ``1000``
      - The maximum number of inner iterations permitted.
+   * - ``max_linesearch_steps``
+     - ``int``, static
+     - ``20``
+     - The maximum number of steps in the LBFGS line search.
+   * - ``svtol``
+     - ``float``, traced
+     - ``1e-7``
+     - Singular-value cut-off threshold during preconditioning.
+   * - ``merge_constraints``
+     - ``bool``, static
+     - ``False``
+     - When ``True``, combines compatible constraint evaluations before solving.
+   * - ``implicit_linear_solver``
+     - ``lineax.AbstractLinearSolver`` or ``None``, static
+     - ``None``
+     - Linear solver used for implicit differentiation.
 
 Thus far, we have successfully run an instance of QUADCOIL. The next section will explain how to interpret the outputs.