Skip to content

Commit

Permalink
Got the henneberg example finalized. Working on a QH example now. Sav…
Browse files Browse the repository at this point in the history
…ing so I can switch to desktop.
  • Loading branch information
akaptano committed Dec 30, 2024
1 parent f7907d4 commit 90ea154
Show file tree
Hide file tree
Showing 10 changed files with 2,034 additions and 54 deletions.
514 changes: 514 additions & 0 deletions examples/3_Advanced/passive_coils/henneberg_experiment.py

Large diffs are not rendered by default.

496 changes: 496 additions & 0 deletions examples/3_Advanced/passive_coils/passive_coils_QH.py

Large diffs are not rendered by default.

406 changes: 406 additions & 0 deletions examples/3_Advanced/passive_coils/passive_coils_QH_continuation.py

Large diffs are not rendered by default.

17 changes: 9 additions & 8 deletions examples/3_Advanced/passive_coils/passive_coils_henneberg.py
Original file line number Diff line number Diff line change
Expand Up @@ -80,8 +80,8 @@ def initialize_coils_QA(TEST_DIR, s):
# generate planar TF coils
ncoils = 2
R0 = s.get_rc(0, 0) * 1.4
R1 = s.get_rc(1, 0) * 4
order = 16
R1 = s.get_rc(1, 0) * 5
order = 8

from simsopt.mhd.vmec import Vmec
vmec_file = 'wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc'
Expand Down Expand Up @@ -122,12 +122,12 @@ def initialize_coils_QA(TEST_DIR, s):
aa = 0.05
bb = 0.05

Nx = 6
Nx = 5
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=False,
s, s_inner, s_outer, Nx, Ny, Nz, order=order, coil_coil_flag=False, jax_flag=False,
)
import warnings

Expand Down Expand Up @@ -158,7 +158,7 @@ def initialize_coils_QA(TEST_DIR, s):
for ii in range(len(base_curves)):
counter = 0
for i in range(base_curves[0].gamma().shape[0]):
eps = 1e-2
eps = 0.2
dij = np.sqrt(np.sum((base_curves[ii].gamma()[i, :]) ** 2))
conflict_bool = (dij < (1.0 + eps) * s.get_rc(0, 0))
if conflict_bool:
Expand Down Expand Up @@ -208,9 +208,9 @@ def initialize_coils_QA(TEST_DIR, s):
for j in range(2 * order + 1):
base_curves[i].fix('x' + str(j))
# Fix center points of each coil
base_curves[i].fix('x' + str(2 * order + 5))
base_curves[i].fix('x' + str(2 * order + 6))
base_curves[i].fix('x' + str(2 * order + 7))
# base_curves[i].fix('x' + str(2 * order + 5))
# base_curves[i].fix('x' + str(2 * order + 6))
# base_curves[i].fix('x' + str(2 * order + 7))

ncoils = len(base_curves)
a_list = np.ones(len(base_curves)) * aa
Expand Down Expand Up @@ -462,6 +462,7 @@ def fun(dofs):
bpsc = btot.Bfields[0]
bpsc.set_points(s_plot.gamma().reshape((-1, 3)))
dipole_currents = [c.current.get_value() for c in bpsc.coils]
psc_array.recompute_currents()
print(dipole_currents)
curves_to_vtk(
[c.curve for c in bpsc.coils],
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,12 @@ def pointData_forces_torques(coils, allcoils, aprimes, bprimes, nturns_list):
return point_data

from simsopt import load
input_dir = "passive_coils_henneberg_ndofs7_TForder16_n23_p1.50e+00_c2.00e+00_lw1.00e-02_lt9.00e+01_lkw1.00e+04_cct8.00e-01_ccw1.00e+02_cst1.30e+00_csw1.00e+02_fw0.00e+00_fww0.000000e+00_tw0.00e+00_tww0.000000e+00/"
input_dir = "passive_coils_henneberg_continuation_ndofs3_TForder16_n7_lw1.00e-02_lt8.00e+01_lkw1.00e+04_cct8.00e-01_ccw1.00e+02_cst1.30e+00_csw1.00e+02_fw0.00e+00_fww0.000000e+00_tw0.00e+00_tww0.000000e+00/"
#"passive_coils_henneberg_continuation_ndofs3_TForder16_n7_lw1.00e-02_lt8.00e+01_lkw1.00e+04_cct5.00e-01_ccw1.00e+02_cst1.30e+00_csw1.00e+02_fw0.00e+00_fww0.000000e+00_tw0.00e+00_tww0.000000e+00/"
#"passive_coils_henneberg_ndofs3_TForder16_n7_p1.30e+00_c4.00e+00_lw1.00e-02_lt9.00e+01_lkw1.00e+03_cct8.00e-01_ccw1.00e+00_cst1.30e+00_csw1.00e+02_fw0.00e+00_fww0.000000e+00_tw0.00e+00_tww0.000000e+00/"
#"passive_coils_henneberg_ndofs3_TForder16_n8_p1.30e+00_c4.00e+00_lw1.00e-02_lt9.00e+01_lkw1.00e+04_cct7.00e-01_ccw1.00e+02_cst1.30e+00_csw1.00e+02_fw0.00e+00_fww0.000000e+00_tw0.00e+00_tww0.000000e+00/"
#"passive_coils_henneberg_ndofs0_TForder20_n37_p1.50e+00_c2.00e+00_lw1.00e-02_lt9.00e+01_lkw1.00e+04_cct8.00e-01_ccw1.00e+02_cst1.30e+00_csw1.00e+02_fw0.00e+00_fww0.000000e+00_tw0.00e+00_tww0.000000e+00/"
#"passive_coils_henneberg_ndofs7_TForder16_n23_p1.50e+00_c2.00e+00_lw1.00e-02_lt9.00e+01_lkw1.00e+04_cct8.00e-01_ccw1.00e+02_cst1.30e+00_csw1.00e+02_fw0.00e+00_fww0.000000e+00_tw0.00e+00_tww0.000000e+00/"
coils = load(input_dir + "psc_coils.json")
coils_TF = load(input_dir + "TF_coils.json")
curves = [c.curve for c in coils]
Expand All @@ -92,6 +97,47 @@ def pointData_forces_torques(coils, allcoils, aprimes, bprimes, nturns_list):
# Initialize the PSCArray object
eval_points = s.gamma().reshape(-1, 3)
psc_array = PSCArray(base_curves, coils_TF, eval_points, a_list, b_list, nfp=s.nfp, stellsym=s.stellsym)
coils = psc_array.coils
coils_TF = psc_array.coils_TF
print([c.current.get_value() for c in coils])
# Remove the coil with small current in it
coils = [c for c in coils if np.abs(c.current.get_value()) > 1e5]
curves = [c.curve for c in coils]
base_curves = curves[:len(curves) // 4]
base_coils = coils[:len(coils) // 4]
curves_TF = [c.curve for c in coils_TF]
base_curves_TF = curves_TF[:len(curves_TF) // 4]
base_coils_TF = coils_TF[:len(coils_TF) // 4]
# Reinitialize the PSCArray object with one less coil
psc_array = PSCArray(base_curves, coils_TF, eval_points, a_list, b_list, nfp=s.nfp, stellsym=s.stellsym)
coils = psc_array.coils
coils_TF = psc_array.coils_TF
curves = [c.curve for c in coils]

order = 0
base_curves = curves[:len(curves) // 4]
for i in range(len(base_curves)):
# unfix orientations of each coil
base_curves[i].unfix('x' + str(2 * order + 1))
base_curves[i].unfix('x' + str(2 * order + 2))
base_curves[i].unfix('x' + str(2 * order + 3))
base_curves[i].unfix('x' + str(2 * order + 4))

# unfix shape of each coil
# for j in range(2 * order + 1):
# base_curves[i].fix('x' + str(j))
# unfix center points of each coil
# base_curves[i].fix('x' + str(2 * order + 5))
# base_curves[i].fix('x' + str(2 * order + 6))
# base_curves[i].fix('x' + str(2 * order + 7))

base_coils = coils[:len(coils) // 4]
curves_TF = [c.curve for c in coils_TF]
base_curves_TF = curves_TF[:len(curves_TF) // 4]
base_coils_TF = coils_TF[:len(coils_TF) // 4]
ncoils = len(base_curves)
print([c.current.get_value() for c in coils])
print(ncoils)

# Calculate average, approximate on-axis B field strength
calculate_on_axis_B(psc_array.biot_savart_TF, s)
Expand All @@ -111,7 +157,7 @@ def pointData_forces_torques(coils, allcoils, aprimes, bprimes, nturns_list):
LENGTH_WEIGHT = Weight(0.01)
LENGTH_TARGET = 80
LINK_WEIGHT = 1e4
CC_THRESHOLD = 0.8
CC_THRESHOLD = 1.0
CC_WEIGHT = 1e2
CS_THRESHOLD = 1.3
CS_WEIGHT = 1e2
Expand All @@ -121,7 +167,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) # 1e-22 Forces are in Newtons, and typical values are ~10^5, 10^6 Newtons
# Directory for output
OUT_DIR = ("./passive_coils_henneberg_continuation_ndofs{:d}_TForder{:d}_n{:d}_lw{:.2e}_lt{:.2e}_lkw{:.2e}" +
OUT_DIR = ("./passive_coils_henneberg_continuation2_ndofs{:d}_TForder{:d}_n{:d}_lw{:.2e}_lt{:.2e}_lkw{:.2e}" +
"_cct{:.2e}_ccw{:.2e}_cst{:.2e}_csw{:.2e}_fw{:.2e}_fww{:2e}_tw{:.2e}_tww{:2e}/").format(
len(base_curves[0].x), base_curves_TF[0].order, ncoils, LENGTH_WEIGHT.value, LENGTH_TARGET, LINK_WEIGHT,
CC_THRESHOLD, CC_WEIGHT, CS_THRESHOLD, CS_WEIGHT, FORCE_WEIGHT.value,
Expand Down Expand Up @@ -200,7 +246,7 @@ def pointData_forces_torques(coils, allcoils, aprimes, bprimes, nturns_list):

CURVATURE_THRESHOLD = 0.5
MSC_THRESHOLD = 0.05
CURVATURE_WEIGHT = 1e-4
CURVATURE_WEIGHT = 1e-1
MSC_WEIGHT = 1e-5
Jcs = [LpCurveCurvature(c.curve, 2, CURVATURE_THRESHOLD) for c in base_coils_TF]
Jmscs = [MeanSquaredCurvature(c.curve) for c in base_coils_TF]
Expand Down Expand Up @@ -303,11 +349,11 @@ def fun(dofs):
""")

n_saves = 1
MAXITER = 1000
MAXITER = 2000
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', # bounds=opt_bounds,
options={'maxiter': MAXITER, 'maxcor': 300}, tol=1e-15)
options={'maxiter': MAXITER, 'maxcor': 500}, tol=1e-15)
# dofs = res.x

bpsc = btot.Bfields[0]
Expand Down
Loading

0 comments on commit 90ea154

Please sign in to comment.