From bf3dab0647d6dee0e23bfd0500628b9f1d1c5042 Mon Sep 17 00:00:00 2001 From: Alan Kaptanoglu Date: Wed, 23 Oct 2024 17:51:13 -0400 Subject: [PATCH] Trying to get the fixed surface dipole example working better. --- .../QA_reactorScale_dipoleArrays.py | 10 +-- .../QA_reactorscale_fixed_dipoles.py | 69 +++++++++++-------- 2 files changed, 45 insertions(+), 34 deletions(-) diff --git a/examples/3_Advanced/QA_reactorScale_dipoleArrays.py b/examples/3_Advanced/QA_reactorScale_dipoleArrays.py index 478d4cf08..07ffa9b7e 100644 --- a/examples/3_Advanced/QA_reactorScale_dipoleArrays.py +++ b/examples/3_Advanced/QA_reactorScale_dipoleArrays.py @@ -41,8 +41,8 @@ range_param = "half period" nphi = 32 ntheta = 32 -poff = 1.5 -coff = 3.0 +poff = 2.25 +coff = 3.5 s = SurfaceRZFourier.from_vmec_input(filename, range=range_param, nphi=nphi, ntheta=ntheta) s_inner = SurfaceRZFourier.from_vmec_input(filename, range=range_param, nphi=nphi * 4, ntheta=ntheta * 4) s_outer = SurfaceRZFourier.from_vmec_input(filename, range=range_param, nphi=nphi * 4, ntheta=ntheta * 4) @@ -307,8 +307,10 @@ def pointData_forces_torques(coils, allcoils, aprimes, bprimes, nturns_list): Jlength = QuadraticPenalty(sum(Jls_TF), LENGTH_TARGET, "max") # coil-coil and coil-plasma distances should be between all coils -Jccdist = CurveCurveDistance(curves + curves_TF, CC_THRESHOLD, num_basecurves=len(coils + coils_TF)) -Jcsdist = CurveSurfaceDistance(curves + curves_TF, s, CS_THRESHOLD) + +### Jcc below removed the dipoles! +Jccdist = CurveCurveDistance(curves_TF, CC_THRESHOLD, num_basecurves=len(coils_TF)) +Jcsdist = CurveSurfaceDistance(curves_TF, s, CS_THRESHOLD) # While the coil array is not moving around, they cannot # interlink. diff --git a/examples/3_Advanced/QA_reactorscale_fixed_dipoles.py b/examples/3_Advanced/QA_reactorscale_fixed_dipoles.py index a22152406..3b5cef953 100644 --- a/examples/3_Advanced/QA_reactorscale_fixed_dipoles.py +++ b/examples/3_Advanced/QA_reactorscale_fixed_dipoles.py @@ -41,15 +41,21 @@ range_param = "half period" nphi = 32 ntheta = 32 -poff = 1.5 -coff = 0.5 +poff = 2.5 +coff = 3.0 s = SurfaceRZFourier.from_vmec_input(filename, range=range_param, nphi=nphi, ntheta=ntheta) s_inner = SurfaceRZFourier.from_vmec_input(TEST_DIR / 'input.circular_tokamak', range=range_param, nphi=nphi * 4, ntheta=ntheta * 4) s_outer = SurfaceRZFourier.from_vmec_input(TEST_DIR / 'input.circular_tokamak', range=range_param, nphi=nphi * 4, ntheta=ntheta * 4) +s_inner.set_rc(0, 0, s.get_rc(0, 0)) +s_inner.set_rc(1, 0, s.get_rc(1, 0) * poff) +s_inner.set_zs(1, 0, s.get_zs(1, 0) * poff) +s_outer.set_rc(0, 0, s.get_rc(0, 0)) +s_outer.set_rc(1, 0, s.get_rc(1, 0) * coff) +s_outer.set_zs(1, 0, s.get_zs(1, 0) * coff) # Make the inner and outer surfaces by extending the plasma surface -s_inner.extend_via_normal(poff) -s_outer.extend_via_normal(poff + coff) +# s_inner.extend_via_normal(poff) +# s_outer.extend_via_normal(poff + coff) qphi = nphi * 2 qtheta = ntheta * 2 @@ -142,12 +148,12 @@ def initialize_coils_QA(TEST_DIR, s): aa = 0.05 bb = 0.05 -Nx = 9 +Nx = 6 Ny = Nx Nz = Nx # Create the initial coils: base_curves, all_curves = create_planar_curves_between_two_toroidal_surfaces( - s, s_inner, s_outer, Nx, Ny, Nz, order=order, coil_coil_flag=True, jax_flag=True, + s, s_inner, s_outer, Nx, Ny, Nz, order=order, coil_coil_flag=False, jax_flag=True, # numquadpoints=10 # Defaults is (order + 1) * 40 so this halves it ) import warnings @@ -183,29 +189,29 @@ def initialize_coils_QA(TEST_DIR, s): point = (base_curves[i].get_dofs()[-3:]) dists = np.sum((point - plasma_points) ** 2, axis=-1) min_ind = np.argmin(dists) - coil_normals[i, :] = plasma_unitnormals[min_ind, :] - # coil_normals[i, :] = -(plasma_points[min_ind, :] - point) -alphas = np.arctan2( - -coil_normals[:, 1], - np.sqrt(coil_normals[:, 0] ** 2 + coil_normals[:, 2] ** 2)) -deltas = np.arcsin(coil_normals[:, 0] / \ - np.sqrt(coil_normals[:, 0] ** 2 + coil_normals[:, 2] ** 2)) + # coil_normals[i, :] = plasma_unitnormals[min_ind, :] + coil_normals[i, :] = (plasma_points[min_ind, :] - point) +# alphas = np.arctan2( +# -coil_normals[:, 1], +# np.sqrt(coil_normals[:, 0] ** 2 + coil_normals[:, 2] ** 2)) +# deltas = np.arcsin(coil_normals[:, 0] / \ +# np.sqrt(coil_normals[:, 0] ** 2 + coil_normals[:, 2] ** 2)) for i in range(len(base_curves)): - alpha2 = alphas[i] / 2.0 - delta2 = deltas[i] / 2.0 - calpha2 = np.cos(alpha2) - salpha2 = np.sin(alpha2) - cdelta2 = np.cos(delta2) - sdelta2 = np.sin(delta2) - base_curves[i].set('x' + str(2 * order + 1), calpha2 * cdelta2) - base_curves[i].set('x' + str(2 * order + 2), salpha2 * cdelta2) - base_curves[i].set('x' + str(2 * order + 3), calpha2 * sdelta2) - base_curves[i].set('x' + str(2 * order + 4), -salpha2 * sdelta2) + # alpha2 = alphas[i] / 2.0 + # delta2 = deltas[i] / 2.0 + # calpha2 = np.cos(alpha2) + # salpha2 = np.sin(alpha2) + # cdelta2 = np.cos(delta2) + # sdelta2 = np.sin(delta2) + # base_curves[i].set('x' + str(2 * order + 1), calpha2 * cdelta2) + # base_curves[i].set('x' + str(2 * order + 2), salpha2 * cdelta2) + # base_curves[i].set('x' + str(2 * order + 3), calpha2 * sdelta2) + # base_curves[i].set('x' + str(2 * order + 4), -salpha2 * sdelta2) # Fix orientations of each coil - base_curves[i].fix('x' + str(2 * order + 1)) - base_curves[i].fix('x' + str(2 * order + 2)) - base_curves[i].fix('x' + str(2 * order + 3)) - base_curves[i].fix('x' + str(2 * order + 4)) + # base_curves[i].fix('x' + str(2 * order + 1)) + # base_curves[i].fix('x' + str(2 * order + 2)) + # base_curves[i].fix('x' + str(2 * order + 3)) + # base_curves[i].fix('x' + str(2 * order + 4)) # Fix shape of each coil for j in range(2 * order + 1): @@ -273,7 +279,7 @@ def pointData_forces_torques(coils, allcoils, aprimes, bprimes, nturns_list): TORQUE_WEIGHT = Weight(0.0) # Forces are in Newtons, and typical values are ~10^5, 10^6 Newtons TORQUE_WEIGHT2 = Weight(0.0) # Forces are in Newtons, and typical values are ~10^5, 10^6 Newtons # Directory for output -OUT_DIR = ("./QA_n{:d}_p{:.2e}_c{:.2e}_lw{:.2e}_lt{:.2e}_lkw{:.2e}" + \ +OUT_DIR = ("./QA_fixeddipoles_n{:d}_p{:.2e}_c{:.2e}_lw{:.2e}_lt{:.2e}_lkw{:.2e}" + \ "_cct{:.2e}_ccw{:.2e}_cst{:.2e}_csw{:.2e}_fw{:.2e}_fww{:2e}_tw{:.2e}_tww{:2e}/").format( ncoils, poff, coff, LENGTH_WEIGHT.value, LENGTH_TARGET, LINK_WEIGHT, CC_THRESHOLD, CC_WEIGHT, CS_THRESHOLD, CS_WEIGHT, FORCE_WEIGHT.value, @@ -284,6 +290,9 @@ def pointData_forces_torques(coils, allcoils, aprimes, bprimes, nturns_list): shutil.rmtree(OUT_DIR) os.makedirs(OUT_DIR, exist_ok=True) +s_inner.to_vtk(OUT_DIR + "s_inner") +s_outer.to_vtk(OUT_DIR + "s_outer") + curves_to_vtk( curves_TF, OUT_DIR + "curves_TF_0", @@ -333,7 +342,7 @@ def pointData_forces_torques(coils, allcoils, aprimes, bprimes, nturns_list): # coil-coil and coil-plasma distances should be between all coils Jccdist = CurveCurveDistance(curves_TF, CC_THRESHOLD, num_basecurves=len(coils + coils_TF)) -Jcsdist = CurveSurfaceDistance(curves_TF, s, CS_THRESHOLD) +Jcsdist = CurveSurfaceDistance(curves + curves_TF, s, CS_THRESHOLD) # While the coil array is not moving around, they cannot # interlink. @@ -522,7 +531,7 @@ def fun(dofs): # print('dJtorques time = ', t2 - t1, ' s') n_saves = 1 -MAXITER = 60 +MAXITER = 200 for i in range(1, n_saves + 1): print('Iteration ' + str(i) + ' / ' + str(n_saves)) res = minimize(fun, dofs, jac=True, method='L-BFGS-B',