recent work

examples/3_Advanced/coil_forces.py

@@ -26,8 +26,8 @@
 
 ncoils = 4
 R0 = 1
-R1 = 1
-order = 5
+R1 = 0.8
+order = 3
 
 LENGTH_WEIGHT = 1e-5
 
@@ -37,7 +37,7 @@
 CS_THRESHOLD = 0.4
 CS_WEIGHT = 10
 
-FORCE_WEIGHT = 1e-20
+FORCE_WEIGHT = 0  # 1e-13
 
 config_str = f"{ncoils}_coils_force_weight_{FORCE_WEIGHT}"
 #######################################################
@@ -53,7 +53,7 @@
 
 # Create the initial coils:
 base_curves = create_equally_spaced_curves(
-    ncoils, s.nfp, stellsym=True, R0=R0, R1=R1, order=order)
+    ncoils, s.nfp, stellsym=True, R0=R0, R1=R1, order=order, numquadpoints=50)
 base_currents = [Current(2.5e5) for i in range(ncoils)]
 # Since the target field is zero, one possible solution is just to set all
 # currents to 0. To avoid the minimizer finding that solution, we fix one
@@ -80,9 +80,6 @@
 
 
 JF = Jf \
-    + LENGTH_WEIGHT * sum(Jls) \
-    + CC_WEIGHT * Jccdist \
-    + CS_WEIGHT * Jcsdist \
     + Jforce * FORCE_WEIGHT
 
 MAXITER = 10
@@ -93,26 +90,23 @@ def fun(dofs):
     JF.x = dofs
     J = JF.J()
     grad = JF.dJ()
-    print(f"J={J:.3e}, ||∇J||={np.linalg.norm(grad):.3e}, J_force={Jforce.J():.3e}, Jflux={Jf.J():.3e}")
+    print(f"J={J:.3e}, ||∇J||={np.linalg.norm(grad):.3e}, J_force={FORCE_WEIGHT*Jforce.J():.3e}, Jflux={Jf.J():.3e}")
     return J, grad
 
 
 res = minimize(fun, dofs, jac=True, method='L-BFGS-B',
                options={'maxiter': MAXITER, 'maxcor': 300}, tol=1e-15)
-curves_to_vtk(curves, OUT_DIR + f"curves_opt_{config_str}")
-print("Coil force optimization")
+curves_to_vtk(base_curves, OUT_DIR + f"curves_opt_{config_str}")
 
-MAXITER = 100
 dofs = res.x
-FORCE_WEIGHT *= 1e6
+print("Force Optimization")
+FORCE_WEIGHT += 1e-13
 
 res = minimize(fun, dofs, jac=True, method='L-BFGS-B',
                options={'maxiter': MAXITER, 'maxcor': 300}, tol=1e-15)
 
+curves_to_vtk(base_curves, OUT_DIR + f"curves_opt_force_{config_str}")
 
-curves_to_vtk(curves, OUT_DIR + f"curves_opt_force_{config_str}")
-curves_to_vtk(base_curves, OUT_DIR +
-              f"curves_opt_hfp_{config_str}")
 pointData = {"B_N": np.sum(bs.B().reshape(
     (nphi, ntheta, 3)) * s.unitnormal(), axis=2)[:, :, None]/np.linalg.norm(bs.B().reshape(
         (nphi, ntheta, 3)), axis=2)}

src/simsopt/field/fieldalignment.py

@@ -1,13 +1,16 @@
 """Implements optimization of the angle between the magnetic field and the ReBCO c-axis."""
 import math
+import pickle
 from scipy import constants
 import numpy as np
 import jax.numpy as jnp
 from jax import grad
 from simsopt.field import BiotSavart, Coil
+from simsopt.field.selffield import field_on_coils
 from simsopt.geo.jit import jit
 from simsopt._core import Optimizable
 from simsopt._core.derivative import derivative_dec
+from simsopt.geo.finitebuild import rotated_centroid_frame
 
 Biot_savart_prefactor = constants.mu_0 / 4 / np.pi
 
@@ -59,10 +62,162 @@ def c_axis_pure(t, b):
 
 def c_axis_angle_pure(coil, B):
     """Angle between the magnetic field and the c-axis of the REBCO Tape"""
-    t, _, b = coil.curve.frenet_frame()
+    t, n, b = rotated_centroid_frame(coil.curve.gamma(
+    ), coil.curve.gammadash(), coil.curve.rotation.alpha(coil.curve.quadpoints))
     b_norm = jnp.einsum("ij, i->ij", B, 1/jnp.linalg.norm(B, axis=1))
 
     c_axis = c_axis_pure(t, b)
     angle = jnp.arccos(inner(c_axis, b_norm))
     return angle
 
+
+def critical_current_obj(coil, JANUS=True):
+    field = field_on_coils(coil)
+    angle = c_axis_angle_pure(coil, field)
+    file_path = "/Users/paulhuslage/PhD/epos-simsopt/src/simsopt/examples/3_Advanced/inputs/Ic_angle_field_interpolation_with units.pck"
+    B_norm = np.linalg.norm(field, axis=1)
+
+    with open(file_path, "rb") as file:
+        # Load the contents of the .pck file
+        data = pickle.load(file)
+
+    ic_interpolation = data['Ic_interpolation']
+    ic_interpolation_reversed = data['Ic_interpolation_reversed']
+
+    if JANUS:
+        Ic = [ic_interpolation_reversed([a, f]) for a, f in zip(angle, B_norm)]
+    else:
+        Ic = [ic_interpolation([a, f]) for a, f in zip(angle, B_norm)]
+
+    return np.min(Ic)
+
+
+def critical_current(coil, JANUS=True):
+    field = field_on_coils(coil)
+    angle = c_axis_angle_pure(coil, field)
+    file_path = "/Users/paulhuslage/PhD/epos-simsopt/src/simsopt/examples/3_Advanced/inputs/Ic_angle_field_interpolation_with units.pck"
+    B_norm = np.linalg.norm(field, axis=1)
+
+    with open(file_path, "rb") as file:
+        # Load the contents of the .pck file
+        data = pickle.load(file)
+
+    ic_interpolation = data['Ic_interpolation']
+    ic_interpolation_reversed = data['Ic_interpolation_reversed']
+
+    if JANUS:
+        Ic = [ic_interpolation_reversed([a, f]) for a, f in zip(angle, B_norm)]
+    else:
+        Ic = [ic_interpolation([a, f]) for a, f in zip(angle, B_norm)]
+
+    return Ic
+
+
+def critical_current_pure(coil, JANUS=True):
+    field = field_on_coils(coil)
+    angle = c_axis_angle_pure(coil, field)
+    file_path = "/Users/paulhuslage/PhD/epos-simsopt/src/simsopt/examples/3_Advanced/inputs/Ic_angle_field_interpolation_with units.pck"
+    B_norm = np.linalg.norm(field, axis=1)
+
+    with open(file_path, "rb") as file:
+        # Load the contents of the .pck file
+        data = pickle.load(file)
+
+    ic_interpolation = data['Ic_interpolation']
+    ic_interpolation_reversed = data['Ic_interpolation_reversed']
+
+    if JANUS:
+        Ic = [ic_interpolation_reversed([a, f]) for a, f in zip(angle, B_norm)]
+    else:
+        Ic = [ic_interpolation([a, f]) for a, f in zip(angle, B_norm)]
+
+    return Ic
+
+
+class CrtitcalCurrentOpt(Optimizable):
+    """Optimizable class to optimize the critical on a ReBCO coil"""
+
+    def __init__(self, coil, coils, a=0.05):
+        self.coil = coil
+        self.curve = coil.curve
+        self.coils = coils
+        self.a = a
+        self.B_ext = BiotSavart(coils).set_points(self.curve.gamma()).B()
+        self.B_self = 0
+        self.B = 0
+        self.alpha = coil.curve.rotation.alpha(coil.curve.quadpoints)
+        self.J_jax = jit(lambda gamma, gammadash, gammadashdash,
+                         current, phi, phidash, B_ext: force_opt_pure(gamma, gammadash, gammadashdash,
+                                                                      current, phi, phidash, B_ext))
+
+        self.thisgrad0 = jit(lambda gamma, gammadash, gammadashdash, current, phi, phidash, B_ext: grad(
+            self.J_jax, argnums=0)(gamma, gammadash, gammadashdash, current, phi, phidash, B_ext))
+        self.thisgrad1 = jit(lambda gamma, gammadash, gammadashdash, current, phi, phidash, B_ext: grad(
+            self.J_jax, argnums=1)(gamma, gammadash, gammadashdash, current, phi, phidash, B_ext))
+        self.thisgrad2 = jit(lambda gamma, gammadash, gammadashdash, current, phi, phidash, B_ext: grad(
+            self.J_jax, argnums=2)(gamma, gammadash, gammadashdash, current, phi, phidash, B_ext))
+
+        super().__init__(depends_on=[coil])
+
+    def J(self):
+        gamma = self.coil.curve.gamma()
+        d1gamma = self.coil.curve.gammadash()
+        d2gamma = self.coil.curve.gammadashdash()
+        current = self.coil.current.get_value()
+        phi = self.coil.curve.quadpoints
+        phidash = self.coil.curve.quadpoints
+        B_ext = self.B_ext
+        return self.J_jax(gamma, d1gamma, d2gamma, current, phi, phidash, B_ext)
+
+    @derivative_dec
+    def dJ(self):
+        gamma = self.coil.curve.gamma()
+        d1gamma = self.coil.curve.gammadash()
+        d2gamma = self.coil.curve.gammadashdash()
+        current = self.coil.current.get_value()
+        phi = self.coil.curve.quadpoints
+        phidash = self.coil.curve.quadpoints
+        B_ext = self.B_ext
+
+        grad0 = self.thisgrad0(gamma, d1gamma, d2gamma,
+                               current, phi, phidash, B_ext)
+        grad1 = self.thisgrad0(gamma, d1gamma, d2gamma,
+                               current, phi, phidash, B_ext)
+        grad2 = self.thisgrad0(gamma, d1gamma, d2gamma,
+                               current, phi, phidash, B_ext)
+
+        return self.coil.curve.dgamma_by_dcoeff_vjp(grad0) + self.coil.curve.dgammadash_by_dcoeff_vjp(grad1) \
+            + self.coil.curve.dgammadashdash_by_dcoeff_vjp(grad2)
+
+    return_fn_map = {'J': J, 'dJ': dJ}
+
+
+def field_alignment_pure(curve, field):
+    angle = c_axis_angle_pure(curve, field)
+    Ic = 0
+    return Ic
+
+
+def grad_field_alignment_pure(curve, field):
+    angle = c_axis_angle_pure(curve, field)
+    return 0
+
+
+class FieldAlignment(Optimizable):
+
+    def init(self, curve, field):
+        Optimizable.__init__(self, depends_on=[curve])
+        self.curve = curve
+        self.field = field
+        # new_dofs = np.zeros(len(dofs))
+        # for i in range(len(dofs)):
+        # new_dofs[i] = dofs[i]
+
+    def J(self):
+
+        return field_alignment_pure(self.curve, self.field)
+
+    @derivative_dec
+    def dJ(self):
+        # return Derivative({self: lambda x0: gradfunction(x0, self.order, self.curves, self.n)})
+        return grad_field_alignment_pure(self.curve, self.field)

src/simsopt/field/selffield.py

@@ -117,7 +117,7 @@ def field_on_coils(coil, a=0.05):
     n_quad = phidash.shape[0]
     gamma = coil.curve.gamma()
     gammadash = coil.curve.gammadash() / 2 / np.pi
-    gammadashdash = coil.curve.gammadashdash() / 4 / np.pi**2
+    gammadashdash = coil.curve.curve.gammadashdash() / 4 / np.pi**2
     integral_term = np.zeros((n_quad, 3))
 
     A = singularity_term_circ(gamma, gammadash, gammadashdash, a)
@@ -157,6 +157,7 @@ def G(x, y):
 
 
 def K(u, v, kappa_1, kappa_2, p, q, a, b):
+    """Auxiliary function for the calculation of the internal field of a rectangular crosssction"""
     K = - 2 * u * v * (kappa_1 * q - kappa_2 * p) * jnp.log(a * u**2 / b + b * v**2 / a) + \
         (kappa_2 * q - kappa_1 * p) * (a * u**2 / b + b * v**2 / a) * jnp.log(a * u**2 / b + b * v**2 / a) + \
         4 * a * u**2 * kappa_2 * p / b * \
@@ -198,15 +199,9 @@ def local_field_rect(coil, u, v, a, b):
     b_b = jnp.zeros((N_phi, 3))
     b_0 = jnp.zeros((N_phi, 3))
 
-    k = (4 * b) / (3 * a) * jnp.arctan(a/b) + (4*a)/(3*b)*jnp.arctan(b/a) + \
-        (b**2)/(6*a**2)*jnp.log(b/a) + (a**2)/(6*b**2)*jnp.log(a/b) - \
-        (a**4 - 6*a**2*b**2 + b**4)/(6*a**2*b**2)*jnp.log(a/b+b/a)
-
-    delta = jnp.exp(-25/6 + k)
-
     for i in range(N_phi):
         b_b.at[i].set(kappa[i] * b[i] / 2 *
-                      (4 + 2*jnp.log(2) + jnp.log(delta)))
+                      (4 + 2*jnp.log(2) + jnp.log(delta(a, b))))
         b_kappa.at[i].set(1 / 16 * (K(u - 1, v - 1, kappa_1[i], kappa_2[i], p[i], q[i], a, b) + K(u + 1, v + 1, kappa_1[i], kappa_2[i], p[i], q[i], a, b)
                                     - K(u - 1, v + 1, kappa_1[i], kappa_2[i], p[i], q[i], a, b) - K(u + 1, v - 1, kappa_1[i], kappa_2[i], p[i], q[i], a, b)))
         b_0.at[i].set(1 / (a * b) * ((G(b * (v - 1), a * (u - 1)) + G(b * (v + 1), a * (u + 1)) - G(b * (v + 1), a * (u - 1)) - G(b * (v - 1), a * (u + 1))) * q -
@@ -278,4 +273,4 @@ def field_from_other_coils_pure(gamma, curves, currents, a=0.05):
     coils = [Coil(curve, current) for curve, current in zip(curves, currents)]
     b_ext = BiotSavart(coils)
     b_ext.set_points(gamma)
-    return b_ext.B()
+    return b_ext.B()