From 1b85d0fb165fbf827697ab9f91f76cbbf0f64538 Mon Sep 17 00:00:00 2001 From: Paul Huslage Date: Wed, 18 Oct 2023 16:00:40 +0200 Subject: [PATCH] field alignment with analyical model for critical current and FramedCurve --- src/simsopt/field/fieldalignment.py | 125 ++++++++-------------------- 1 file changed, 37 insertions(+), 88 deletions(-) diff --git a/src/simsopt/field/fieldalignment.py b/src/simsopt/field/fieldalignment.py index 681ed3dde..2eda9efbc 100755 --- a/src/simsopt/field/fieldalignment.py +++ b/src/simsopt/field/fieldalignment.py @@ -6,11 +6,12 @@ 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.field.selffield import B_regularized, regularization_rect, regularization_circ, B_regularized_pure +from simsopt.geo import FramedCurve +from simsopt.geo.framedcurve import rotated_frenet_frame 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 @@ -62,82 +63,61 @@ 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, n, b = rotated_centroid_frame(coil.curve.gamma( - ), coil.curve.gammadash(), coil.curve.rotation.alpha(coil.curve.quadpoints)) + t, _, b = coil.curve.rotated_frame() 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 angle - return np.min(Ic) +def critical_current_pure(gamma, gammadash, gammadashdash, alpha, quadpoints, current, a, b): -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) + regularization = regularization_rect(a, b) + tangent, normal, binormal = rotated_frenet_frame( + gamma, gammadash, gammadashdash, alpha) - with open(file_path, "rb") as file: - # Load the contents of the .pck file - data = pickle.load(file) + # Fit parameters for reduced Kim-like model of the critical current (doi:10.1088/0953-2048/24/6/065005) + xi = -0.3 + k = 0.7 + B_0 = 42.6e-3 + Ic_0 = 1 - ic_interpolation = data['Ic_interpolation'] - ic_interpolation_reversed = data['Ic_interpolation_reversed'] + field = B_regularized_pure( + gamma, gammadash, gammadashdash, quadpoints, current, regularization) + B_proj = field - inner(field, tangent)[:, None] * tangent - 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)] + B_perp = inner(B_proj, normal) + B_par = inner(B_proj, binormal) + Ic = Ic_0*(jnp.sqrt(k**2 * B_par**2 + B_perp**2) / B_0)**xi 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) +def critical_current(framedcoil, a, b, JANUS=True): + """Critical current on a coil assuming a homogeneous current distribution. Replace with analytical model""" + Ic = critical_current_pure(framedcoil.curve.curve.gamma(), framedcoil.curve.curve.gammadash(), framedcoil.curve.curve.gammadashdash( + ), framedcoil.curve.rotation.alpha(framedcoil.curve.quadpoints), framedcoil.curve.quadpoints, framedcoil.current.current, a, b) + return Ic - 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'] +def critical_current_obj_pure(gamma, gammadash, gammadashdash, alpha, quadpoints, current, a, b): + Ic = critical_current_pure( + gamma, gammadash, gammadashdash, alpha, quadpoints, current, a, b) + obj = np.min(Ic) + return obj - 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_obj(framedcoil, a, b, JANUS=True): + """Objective for field alignement optimization: Target minimum of the critical current along the coil""" + return np.min(critical_current(framedcoil, a, b)) class CrtitcalCurrentOpt(Optimizable): """Optimizable class to optimize the critical on a ReBCO coil""" - def __init__(self, coil, coils, a=0.05): + def __init__(self, framedcoil, coils, a=0.05): self.coil = coil self.curve = coil.curve self.coils = coils @@ -146,9 +126,9 @@ def __init__(self, coil, coils, a=0.05): 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.quadpoints = coil.curve.quadpoints + self.J_jax = jit(lambda gamma, gammadash, gammadashdash, alpha, quadpoints, current, a, + b: critical_current_obj_pure(gamma, gammadash, gammadashdash, alpha, quadpoints, current, a, b)) 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)) @@ -190,34 +170,3 @@ def dJ(self): + 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)