merged latest changes from master

docs/source/example_vmec_only.rst

@@ -17,12 +17,13 @@ directory
 For this problem the two independent variables are ``RBC(1,1)`` and
 ``ZBS(1,1)``, which control the shape of the plasma boundary::
   
-  R(phi) = 1 + 0.1 * cos(theta) + RBC(1,1) * cos(theta - 5 * phi),
-  Z(phi) =     0.1 * sin(theta) + ZBS(1,1) * sin(theta - 5 * phi).
+  R(theta, phi) = 1 + 0.1 * cos(theta) + RBC(1,1) * cos(theta - 5 * phi),
+  Z(theta, phi) =     0.1 * sin(theta) + ZBS(1,1) * sin(theta - 5 * phi).
 
 Note that this boundary has five field periods. We consider the vacuum
 field inside this boundary magnetic surface, i.e. there is no plasma
-current or pressure.  The objective function is
+current or pressure. The toroidal flux enclosed by the boundary
+surface is held fixed. The objective function is
 
 .. code-block::
    

docs/source/simsopt.geo.rst

@@ -76,6 +76,14 @@ simsopt.geo.finitebuild module
    :undoc-members:
    :show-inheritance:
 
+simsopt.geo.framedcurve module
+------------------------------
+
+.. automodule:: simsopt.geo.framedcurve
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 simsopt.geo.jit module
 ----------------------
 
@@ -108,6 +116,14 @@ simsopt.geo.qfmsurface module
    :undoc-members:
    :show-inheritance:
 
+simsopt.geo.strain_optimization module
+--------------------------------------
+
+.. automodule:: simsopt.geo.strain_optimization
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 simsopt.geo.surface module
 --------------------------
 

examples/2_Intermediate/stage_two_optimization.py

@@ -25,7 +25,6 @@
 from pathlib import Path
 import numpy as np
 from scipy.optimize import minimize
-
 from simsopt.field import BiotSavart, Current, coils_via_symmetries
 from simsopt.geo import (SurfaceRZFourier, curves_to_vtk, create_equally_spaced_curves,
                          CurveLength, CurveCurveDistance, MeanSquaredCurvature,

examples/2_Intermediate/strain_optimization.py

@@ -0,0 +1,60 @@
+#!/usr/bin/env python
+
+"""
+This script performs an optimization of the HTS tape winding angle 
+with respect to binormal curvature and torsional strain cost functions as defined in
+
+    Paz Soldan, "Non-planar coil winding angle optimization for compatibility with 
+    non-insulated high-temperature superconducting magnets", Journal of Plasma Physics 
+    86 (2020), doi:10.1017/S0022377820001208.
+    
+The orientation of the tape is defined with respect to the Frenet-Serret Frame
+"""
+
+import numpy as np
+from scipy.optimize import minimize
+from simsopt.geo import CoilStrain, LPTorsionalStrainPenalty, LPBinormalCurvatureStrainPenalty
+from simsopt.geo import FrameRotation, FramedCurveFrenet, CurveXYZFourier
+from simsopt.configs import get_hsx_data
+from simsopt.util import in_github_actions
+
+MAXITER = 50 if in_github_actions else 400
+
+curves, currents, ma = get_hsx_data(Nt_coils=10, ppp=10)
+curve = curves[1]
+scale_factor = 0.1
+curve_scaled = CurveXYZFourier(curve.quadpoints, curve.order)
+curve_scaled.x = curve.x * scale_factor  # scale coil to magnify the strains
+rot_order = 10  # order of the Fourier expression for the rotation of the filament pack
+width = 1e-3  # tape width
+
+curve_scaled.fix_all()  # fix curve DOFs -> only optimize winding angle
+rotation = FrameRotation(curve_scaled.quadpoints, rot_order)
+
+framedcurve = FramedCurveFrenet(curve_scaled, rotation)
+
+tor_threshold = 0.02  # Threshold for strain parameters
+cur_threshold = 0.02
+
+Jtor = LPTorsionalStrainPenalty(framedcurve, p=2, threshold=tor_threshold)
+Jbin = LPBinormalCurvatureStrainPenalty(
+    framedcurve, p=2, threshold=cur_threshold)
+
+strain = CoilStrain(framedcurve, width)
+JF = Jtor + Jbin
+
+
+def fun(dofs):
+    JF.x = dofs
+    J = JF.J()
+    grad = JF.dJ()
+    outstr = f"Max torsional strain={np.max(strain.torsional_strain()):.1e}, Max curvature strain={np.max(strain.binormal_curvature_strain()):.1e}"
+    print(outstr)
+    return J, grad
+
+
+f = fun
+dofs = JF.x
+
+res = minimize(fun, dofs, jac=True, method='L-BFGS-B',
+               options={'maxiter': MAXITER, 'maxcor': 10, 'gtol': 1e-20, 'ftol': 1e-20}, tol=1e-20)

examples/run_serial_examples

@@ -18,6 +18,7 @@ set -ex
 ./2_Intermediate/stage_two_optimization.py
 ./2_Intermediate/stage_two_optimization_stochastic.py
 ./2_Intermediate/stage_two_optimization_finite_beta.py
+./2_Intermediate/strain_optimization.py
 ./2_Intermediate/permanent_magnet_MUSE.py
 ./2_Intermediate/permanent_magnet_QA.py
 ./2_Intermediate/permanent_magnet_PM4Stell.py

src/simsopt/field/normal_field.py

@@ -101,18 +101,18 @@ def from_spec(cls, filename):
         if ph['istellsym']:
             vnc = None
         else:
-            vnc = np.asarray(ph['vnc'][1:])
+            vnc = np.asarray(ph['vnc'])
 
-        nf = cls(
-            nfp=ph['nfp'], 
-            stellsym=ph['istellsym'], 
-            mpol=ph['Mpol'], 
+        normal_field = cls(
+            nfp=ph['nfp'],
+            stellsym=bool(ph['istellsym']),
+            mpol=ph['Mpol'],
             ntor=ph['Ntor'],
             vns=vns,
             vnc=vnc
         )
 
-        return nf
+        return normal_field
 
     def get_index_in_dofs(self, m, n, mpol=None, ntor=None, even=False):
         """

src/simsopt/field/tracing.py

@@ -799,7 +799,8 @@ class IterationStoppingCriterion(sopp.IterationStoppingCriterion):
     pass
 
 
-def plot_poincare_data(fieldlines_phi_hits, phis, filename, mark_lost=False, aspect='equal', dpi=300, xlims=None, ylims=None, surf=None):
+def plot_poincare_data(fieldlines_phi_hits, phis, filename, mark_lost=False, aspect='equal', dpi=300, xlims=None, 
+                       ylims=None, surf=None, s=2, marker='o'):
     """
     Create a poincare plot. Usage:
 
@@ -818,6 +819,8 @@ def plot_poincare_data(fieldlines_phi_hits, phis, filename, mark_lost=False, asp
     nrowcol = ceil(sqrt(len(phis)))
     plt.figure()
     fig, axs = plt.subplots(nrowcol, nrowcol, figsize=(8, 5))
+    for ax in axs.ravel():
+        ax.set_aspect(aspect)
     color = None
     for i in range(len(phis)):
         row = i//nrowcol
@@ -844,7 +847,7 @@ def plot_poincare_data(fieldlines_phi_hits, phis, filename, mark_lost=False, asp
             if data_this_phi.size == 0:
                 continue
             r = np.sqrt(data_this_phi[:, 2]**2+data_this_phi[:, 3]**2)
-            axs[row, col].scatter(r, data_this_phi[:, 4], marker='o', s=2, linewidths=0, c=color)
+            axs[row, col].scatter(r, data_this_phi[:, 4], marker=marker, s=s, linewidths=0, c=color)
 
         plt.rc('axes', axisbelow=True)
         axs[row, col].grid(True, linewidth=0.5)

src/simsopt/geo/__init__.py

@@ -8,7 +8,7 @@
 from .curvexyzfourier import *
 from .curveperturbed import *
 from .curveobjectives import *
-
+from .framedcurve import *
 from .finitebuild import *
 from .plotting import *
 
@@ -21,6 +21,7 @@
 from .surfacerzfourier import *
 from .surfacexyzfourier import *
 from .surfacexyztensorfourier import *
+from .strain_optimization import *
 
 from .permanent_magnet_grid import *
 from .windowpanecurve import *
@@ -35,3 +36,4 @@
            surfacerzfourier.__all__ + surfacexyzfourier.__all__ +
            surfacexyztensorfourier.__all__ + surfaceobjectives.__all__ +
            permanent_magnet_grid.__all__ + windowpanecurve.__all__)
+           strain_optimization.__all__ + framedcurve.__all__)