Strain optimization

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/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 *
 
@@ -33,4 +34,5 @@
            surfacegarabedian.__all__ + surfacehenneberg.__all__ +
            surfacerzfourier.__all__ + surfacexyzfourier.__all__ +
            surfacexyztensorfourier.__all__ + surfaceobjectives.__all__ +
+           strain_optimization.__all__ + framedcurve.__all__ + 
            permanent_magnet_grid.__all__)