Strain optimization

CMakeLists.txt

@@ -51,7 +51,7 @@ else()
     unset(COMPILER_SUPPORTS_MARCH_NATIVE CACHE)
     CHECK_CXX_COMPILER_FLAG(-march=native COMPILER_SUPPORTS_MARCH_NATIVE)
     if(COMPILER_SUPPORTS_MARCH_NATIVE)
-    	set(CMAKE_CXX_FLAGS "-O3 -march=native -mfma -ffp-contract=fast")
+    	set(CMAKE_CXX_FLAGS "-O3 -mfma -ffp-contract=fast")
     elseif(${CMAKE_HOST_SYSTEM_PROCESSOR} STREQUAL "arm64")  
     	set(CMAKE_CXX_FLAGS "-O3 -mcpu=apple-a14 -mfma -ffp-contract=fast")
     else()

examples/2_Intermediate/strain_optimization_simple.py

@@ -0,0 +1,61 @@
+"""
+This script performs an optimization of the HTS tape winding angle 
+with respect to binormal curvature and torsional strain cost functions. 
+The orientation of the tape is defined wrt the Frener-Serret Frame
+"""
+
+import numpy as np
+import os
+import numpy as np
+from scipy.optimize import minimize
+from simsopt.geo import StrainOpt, LPTorsionalStrainPenalty, LPBinormalCurvatureStrainPenalty
+from simsopt.geo import FrameRotation, FramedCurveFrenet, CurveXYZFourier
+from simsopt.field import load_coils_from_makegrid_file
+from simsopt.configs import get_hsx_data
+import matplotlib.pyplot as plt
+
+ci = "CI" in os.environ and os.environ['CI'].lower() in ['1', 'true']
+MAXITER = 50 if ci 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 = StrainOpt(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
+
+print("""
+################################################################################
+### Run the optimisation #######################################################
+################################################################################
+""")
+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/3_Advanced/strain_simple.py

@@ -0,0 +1,52 @@
+"""
+This script performs a stage two coil optimization
+with additional terms for torsion and binormal curvature in the cost function.
+We use a multifilament approach that initializes the coils using the Frener-Serret Frame
+"""
+
+import os
+from pathlib import Path
+import numpy as np
+from scipy.optimize import minimize
+from simsopt.field import BiotSavart
+from simsopt.field import Current, Coil, apply_symmetries_to_curves, apply_symmetries_to_currents
+from simsopt.geo import curves_to_vtk, create_equally_spaced_curves
+from simsopt.geo import CurveLength, CurveCurveDistance, CurveSurfaceDistance
+from simsopt.geo import SurfaceRZFourier
+from simsopt.geo import ArclengthVariation
+from simsopt.objectives import SquaredFlux
+from simsopt.objectives import QuadraticPenalty
+from simsopt.geo.strain_optimization_classes import StrainOpt
+from simsopt.geo import create_multifilament_grid, ZeroRotation, FramedCurveFrenet, FrameRotation
+# from exportcoils import export_coils, import_coils, import_coils_fb, export_coils_fb
+from simsopt.configs import get_ncsx_data
+import matplotlib.pyplot as plt
+
+curves, currents, ma = get_ncsx_data()
+curve = curves[0]
+
+# Set up the winding pack
+rot_order = 5  # order of the Fourier expression for the rotation of the filament pack
+
+width = 12
+
+rotation = FrameRotation(curve.quadpoints, rot_order)
+framedcurve = FramedCurveFrenet(curve, rotation)
+
+strain = StrainOpt(framedcurve, width=width)
+
+tor = strain.binormal_curvature_strain()
+
+plt.figure()
+plt.plot(tor)
+plt.show()
+
+# tor = strain.torsional_strain()
+tor_frame = framedcurve.frame_torsion()
+tor_curve = curve.torsion()
+dl = curve.incremental_arclength()
+
+plt.figure()
+plt.plot(tor_frame)
+plt.plot(tor_curve)
+plt.show()

src/simsopt/field/magneticfieldclasses.py

@@ -257,8 +257,8 @@ def __init__(self, phi_str):
         self.phi_str = phi_str
         self.phi_parsed = parse_expr(phi_str)
         R, Z, Phi = sp.symbols('R Z phi')
-        self.Blambdify = sp.lambdify((R, Z, Phi), [self.phi_parsed.diff(R)+1e-30*Phi*R*Z,\
-                                                   self.phi_parsed.diff(Phi)/R+1e-30*Phi*R*Z,\
+        self.Blambdify = sp.lambdify((R, Z, Phi), [self.phi_parsed.diff(R)+1e-30*Phi*R*Z, \
+                                                   self.phi_parsed.diff(Phi)/R+1e-30*Phi*R*Z, \
                                                    self.phi_parsed.diff(Z)+1e-30*Phi*R*Z])
         self.dBlambdify_by_dX = sp.lambdify(
             (R, Z, Phi),

src/simsopt/field/self_field.py

@@ -0,0 +1 @@
+# self field calculation

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 *
 
 __all__ = (curve.__all__ + curvehelical.__all__ +
            curverzfourier.__all__ + curvexyzfourier.__all__ +