diff --git a/examples/3_Advanced/coil_forces.py b/examples/3_Advanced/coil_forces.py index 2ca07bb1e..dc09d47c3 100644 --- a/examples/3_Advanced/coil_forces.py +++ b/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)} diff --git a/src/simsopt/field/fieldalignment.py b/src/simsopt/field/fieldalignment.py index 7e4e9d87c..681ed3dde 100755 --- a/src/simsopt/field/fieldalignment.py +++ b/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) diff --git a/src/simsopt/field/selffield.py b/src/simsopt/field/selffield.py index 6873d8ae4..18245f1e5 100755 --- a/src/simsopt/field/selffield.py +++ b/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() \ No newline at end of file + return b_ext.B()