From 1c49be10958f7c14874c70b2af60089923eddaed Mon Sep 17 00:00:00 2001 From: Alan Kaptanoglu Date: Mon, 20 May 2024 19:51:10 -0400 Subject: [PATCH] Added a grid for window pane coils, mostly to simplify the optimization problem and figure out why BFGS is having issues. Can now optimize window panes. When only optimizing coil currents, seems to work fine. Issue is clearly with A_deriv function somehow (the issue is that finite differences outperform BFGS on small problems because it is not getting the downhill direction quite right). Rescaled things so that mu0 factors only get multiplied at the very end. Hopefully reduces numerical accumulation. May need to remove in future and just optimized the unnormalized quantities. --- examples/2_Intermediate/QA_psc_example.py | 329 +++++++++++++++++ examples/2_Intermediate/QA_wp_example.py | 344 ++++++++++++++++++ examples/2_Intermediate/QH_psc_example.py | 20 +- src/simsopt/geo/psc_grid.py | 36 +- .../util/permanent_magnet_helper_functions.py | 17 +- src/simsoptpp/psc.cpp | 8 +- src/simsoptpp/psc.h | 1 - tests/geo/test_psc_grid.py | 37 +- 8 files changed, 737 insertions(+), 55 deletions(-) create mode 100644 examples/2_Intermediate/QA_psc_example.py create mode 100644 examples/2_Intermediate/QA_wp_example.py diff --git a/examples/2_Intermediate/QA_psc_example.py b/examples/2_Intermediate/QA_psc_example.py new file mode 100644 index 000000000..137471ba8 --- /dev/null +++ b/examples/2_Intermediate/QA_psc_example.py @@ -0,0 +1,329 @@ +#!/usr/bin/env python +r""" +This example script optimizes a set of relatively simple toroidal field coils +and passive superconducting coils (PSCs) +for an ARIES-CS reactor-scale version of the precise-QH stellarator from +Landreman and Paul. + +The script should be run as: + mpirun -n 1 python QH_psc_example.py +on a cluster machine but + python QH_psc_example.py +is sufficient on other machines. Note that this code does not use MPI, but is +parallelized via OpenMP, so will run substantially +faster on multi-core machines (make sure that all the cores +are available to OpenMP, e.g. through setting OMP_NUM_THREADS). + +""" +from pathlib import Path + +import numpy as np +from scipy.optimize import minimize + +from simsopt.field import BiotSavart, Current, coils_via_symmetries +from simsopt.geo import SurfaceRZFourier, curves_to_vtk +from simsopt.geo.psc_grid import PSCgrid +from simsopt.objectives import SquaredFlux +from simsopt.util import in_github_actions +from simsopt.util.permanent_magnet_helper_functions import * +import time + +def coil_optimization_QA(s, bs, base_curves, curves, out_dir=''): + """ + Optimize the coils for the QA, QH, or other configurations. + + Args: + s: plasma boundary. + bs: Biot Savart class object, presumably representing the + magnetic fields generated by the coils. + base_curves: List of CurveXYZFourier class objects. + curves: List of Curve class objects. + out_dir: Path or string for the output directory for saved files. + + Returns: + bs: Biot Savart class object, presumably representing the + OPTIMIZED magnetic fields generated by the coils. + """ + + from simsopt.geo import CurveLength, CurveCurveDistance, \ + MeanSquaredCurvature, LpCurveCurvature, CurveSurfaceDistance + from simsopt.objectives import QuadraticPenalty + from simsopt.geo import curves_to_vtk + from simsopt.objectives import SquaredFlux + + out_dir = Path(out_dir) + nphi = len(s.quadpoints_phi) + ntheta = len(s.quadpoints_theta) + ncoils = len(base_curves) + + # Weight on the curve lengths in the objective function: + LENGTH_WEIGHT = 1e-4 + + # Threshold and weight for the coil-to-coil distance penalty in the objective function: + CC_THRESHOLD = 0.2 + CC_WEIGHT = 1e1 + + # Threshold and weight for the coil-to-surface distance penalty in the objective function: + CS_THRESHOLD = 0.5 + CS_WEIGHT = 1 + + # Threshold and weight for the curvature penalty in the objective function: + CURVATURE_THRESHOLD = 10 + CURVATURE_WEIGHT = 1e-10 + + # Threshold and weight for the mean squared curvature penalty in the objective function: + MSC_THRESHOLD = 10 + MSC_WEIGHT = 1e-10 + + MAXITER = 500 # number of iterations for minimize + + # Define the objective function: + Jf = SquaredFlux(s, bs) + Jls = [CurveLength(c) for c in base_curves] + Jccdist = CurveCurveDistance(curves, CC_THRESHOLD, num_basecurves=ncoils) + Jcsdist = CurveSurfaceDistance(curves, s, CS_THRESHOLD) + Jcs = [LpCurveCurvature(c, 2, CURVATURE_THRESHOLD) for c in base_curves] + Jmscs = [MeanSquaredCurvature(c) for c in base_curves] + + # Form the total objective function. + JF = Jf \ + + LENGTH_WEIGHT * sum(Jls) \ + + CC_WEIGHT * Jccdist \ + + CS_WEIGHT * Jcsdist \ + + CURVATURE_WEIGHT * sum(Jcs) \ + + MSC_WEIGHT * sum(QuadraticPenalty(J, MSC_THRESHOLD) for J in Jmscs) + + def fun(dofs): + """ Function for coil optimization grabbed from stage_two_optimization.py """ + JF.x = dofs + J = JF.J() + grad = JF.dJ() + jf = Jf.J() + BdotN = np.mean(np.abs(np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2))) + outstr = f"J={J:.1e}, Jf={jf:.1e}, ⟨B·n⟩={BdotN:.1e}" + cl_string = ", ".join([f"{J.J():.1f}" for J in Jls]) + kap_string = ", ".join(f"{np.max(c.kappa()):.1f}" for c in base_curves) + msc_string = ", ".join(f"{J.J():.1f}" for J in Jmscs) + outstr += f", Len=sum([{cl_string}])={sum(J.J() for J in Jls):.1f}, ϰ=[{kap_string}], ∫ϰ²/L=[{msc_string}]" + outstr += f", C-C-Sep={Jccdist.shortest_distance():.2f}, C-S-Sep={Jcsdist.shortest_distance():.2f}" + outstr += f", ║∇J║={np.linalg.norm(grad):.1e}" + print(outstr) + return J, grad + + print(""" + ################################################################################ + ### Perform a Taylor test ###################################################### + ################################################################################ + """) + f = fun + dofs = JF.x + np.random.seed(1) + h = np.random.uniform(size=dofs.shape) + + J0, dJ0 = f(dofs) + dJh = sum(dJ0 * h) + for eps in [1e-3, 1e-4, 1e-5, 1e-6, 1e-7]: + J1, _ = f(dofs + eps*h) + J2, _ = f(dofs - eps*h) + print("err", (J1-J2)/(2*eps) - dJh) + + print(""" + ################################################################################ + ### Run the optimisation ####################################################### + ################################################################################ + """) + minimize(fun, dofs, jac=True, method='L-BFGS-B', options={'maxiter': MAXITER, 'maxcor': 300}, tol=1e-15) + curves_to_vtk(curves, out_dir / "curves_opt") + bs.set_points(s.gamma().reshape((-1, 3))) + return bs + + +np.random.seed(1) # set a seed so that the same PSCs are initialized each time + +# Set some parameters -- if doing CI, lower the resolution +if in_github_actions: + nphi = 4 # nphi = ntheta >= 64 needed for accurate full-resolution runs + ntheta = nphi +else: + # Resolution needs to be reasonably high if you are doing permanent magnets + # or small coils because the fields are quite local + nphi = 32 # nphi = ntheta >= 64 needed for accurate full-resolution runs + ntheta = nphi + # Make higher resolution surface for plotting Bnormal + qphi = nphi * 4 + quadpoints_phi = np.linspace(0, 1, qphi, endpoint=True) + quadpoints_theta = np.linspace(0, 1, ntheta * 4, endpoint=True) + +poff = 0.25 # PSC grid will be offset 'poff' meters from the plasma surface +coff = 0.3 # PSC grid will be initialized between 1 m and 2 m from plasma + +# Read in the plasma equilibrium file +input_name = 'input.LandremanPaul2021_QA_lowres' +TEST_DIR = (Path(__file__).parent / ".." / ".." / "tests" / "test_files").resolve() +surface_filename = TEST_DIR / input_name +range_param = 'half period' +s = SurfaceRZFourier.from_vmec_input( + surface_filename, range=range_param, nphi=nphi, ntheta=ntheta +) +# Print major and minor radius +print('s.R = ', s.get_rc(0, 0)) +print('s.r = ', s.get_rc(1, 0)) + +# Make inner and outer toroidal surfaces very high resolution, +# which helps to initialize coils precisely between the surfaces. +s_inner = SurfaceRZFourier.from_vmec_input( + surface_filename, range=range_param, nphi=nphi * 4, ntheta=ntheta * 4 +) +s_outer = SurfaceRZFourier.from_vmec_input( + surface_filename, range=range_param, nphi=nphi * 4, ntheta=ntheta * 4 +) + +# Make the inner and outer surfaces by extending the plasma surface +s_inner.extend_via_normal(poff) +s_outer.extend_via_normal(poff + coff) + +# Make the output directory +out_str = "QA_psc_output/" +out_dir = Path("QA_psc_output") +out_dir.mkdir(parents=True, exist_ok=True) + +# Save the inner and outer surfaces for debugging purposes +s_inner.to_vtk(out_str + 'inner_surf') +s_outer.to_vtk(out_str + 'outer_surf') + +# initialize the coils +base_curves, curves, coils = initialize_coils('qa_psc', TEST_DIR, s, out_dir) +currents = np.array([coil.current.get_value() for coil in coils]) +print('Currents = ', currents) + +# Set up BiotSavart fields +bs = BiotSavart(coils) + +# Calculate average, approximate on-axis B field strength +calculate_on_axis_B(bs, s) + +# Make high resolution, full torus version of the plasma boundary for plotting +s_plot = SurfaceRZFourier.from_vmec_input( + surface_filename, + quadpoints_phi=quadpoints_phi, + quadpoints_theta=quadpoints_theta +) + +# Plot initial Bnormal on plasma surface from un-optimized BiotSavart coils +make_Bnormal_plots(bs, s_plot, out_dir, "biot_savart_initial") + +# optimize the currents in the TF coils and plot results +# fix all the coil shapes so only the currents are optimized +# for i in range(ncoils): +# base_curves[i].fix_all() +bs = coil_optimization_QA(s, bs, base_curves, curves, out_dir) +curves_to_vtk(curves, out_dir / "TF_coils", close=True) +bs.set_points(s.gamma().reshape((-1, 3))) +Bnormal = np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2) +make_Bnormal_plots(bs, s_plot, out_dir, "biot_savart_TF_optimized") + +# check after-optimization average on-axis magnetic field strength +B_axis = calculate_on_axis_B(bs, s) +make_Bnormal_plots(bs, s_plot, out_dir, "biot_savart_TF_optimized", B_axis) + +# Finally, initialize the psc class +kwargs_geo = {"Nx": 5, "out_dir": out_str, + "random_initialization": True, "poff": poff,} + # "interpolated_field": True} +psc_array = PSCgrid.geo_setup_between_toroidal_surfaces( + s, coils, s_inner, s_outer, **kwargs_geo +) +print('Number of PSC locations = ', len(psc_array.grid_xyz)) + +currents = [] +for i in range(psc_array.num_psc): + currents.append(Current(psc_array.I[i])) +all_coils = coils_via_symmetries( + psc_array.curves, currents, nfp=psc_array.nfp, stellsym=psc_array.stellsym +) +B_PSC = BiotSavart(all_coils) +# Plot initial errors from only the PSCs, and then together with the TF coils +make_Bnormal_plots(B_PSC, s_plot, out_dir, "biot_savart_PSC_initial", B_axis) +make_Bnormal_plots(bs + B_PSC, s_plot, out_dir, "PSC_and_TF_initial", B_axis) + +# Check SquaredFlux values using different ways to calculate it +x0 = np.ravel(np.array([psc_array.alphas, psc_array.deltas])) +fB = SquaredFlux(s, bs, np.zeros((nphi, ntheta))).J() +print('fB only TF coils = ', fB / (B_axis ** 2 * s.area())) +psc_array.least_squares(np.zeros(x0.shape)) +bs.set_points(s.gamma().reshape(-1, 3)) +Bnormal = np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2) +print('fB only TF direct = ', np.sum(Bnormal.reshape(-1) ** 2 * psc_array.grid_normalization ** 2 + ) / (2 * B_axis ** 2 * s.area())) +make_Bnormal_plots(B_PSC, s_plot, out_dir, "biot_savart_PSC_initial", B_axis) +fB = SquaredFlux(s, B_PSC, np.zeros((nphi, ntheta))).J() +print(fB/ (B_axis ** 2 * s.area())) +fB = SquaredFlux(s, B_PSC + bs, np.zeros((nphi, ntheta))).J() +print('fB with both, before opt = ', fB / (B_axis ** 2 * s.area())) +fB = SquaredFlux(s, B_PSC, -Bnormal).J() +print('fB with both (minus sign), before opt = ', fB / (B_axis ** 2 * s.area())) + +# Actually do the minimization now +print('beginning optimization: ') +opt_bounds = tuple([(0, 2 * np.pi) for i in range(psc_array.num_psc * 2)]) +options = {"disp": True, "maxiter": 100} +# print(opt_bounds) +# x0 = np.random.rand(2 * psc_array.num_psc) * 2 * np.pi +verbose = True + +# Run STLSQ with BFGS in the loop +kwargs_manual = { + "out_dir": out_str, + "plasma_boundary" : s, + "coils_TF" : coils + } +I_threshold = 0.0 +STLSQ_max_iters = 10 +for k in range(STLSQ_max_iters): + x0 = np.ravel(np.array([psc_array.alphas, psc_array.deltas])) + print('Number of PSCs = ', len(x0) // 2, ' in iteration ', k) + x_opt = minimize(psc_array.least_squares, x0, args=(verbose,), + method='L-BFGS-B', + # bounds=opt_bounds, + jac=psc_array.least_squares_jacobian, + options=options, + tol=1e-10, + ) + I = psc_array.I + small_I_inds = np.ravel(np.where(np.abs(I) < I_threshold)) + grid_xyz = psc_array.grid_xyz + alphas = psc_array.alphas + deltas = psc_array.deltas + if len(small_I_inds) > 0: + grid_xyz = grid_xyz[small_I_inds, :] + alphas = alphas[small_I_inds] + deltas = deltas[small_I_inds] + else: + print('STLSQ converged, breaking out of loop') + break + kwargs_manual["alphas"] = alphas + kwargs_manual["deltas"] = deltas + # Initialize new PSC array with coils only at the remaining locations + # with initial orientations from the solve using BFGS + psc_array = PSCgrid.geo_setup_manual( + grid_xyz, psc_array.R, **kwargs_manual + ) + +# psc_array.setup_orientations(x_opt.x[:len(x_opt) // 2], x_opt.x[len(x_opt) // 2:]) +psc_array.setup_curves() +psc_array.plot_curves('final_') +currents = [] +for i in range(psc_array.num_psc): + currents.append(Current(psc_array.I[i])) +all_coils = coils_via_symmetries( + psc_array.curves, currents, nfp=psc_array.nfp, stellsym=psc_array.stellsym +) +B_PSC = BiotSavart(all_coils) + +# Check that direct Bn calculation agrees with optimization calculation +fB = SquaredFlux(s, B_PSC + bs, np.zeros((nphi, ntheta))).J() +print('fB with both, after opt = ', fB / (B_axis ** 2 * s.area())) +make_Bnormal_plots(B_PSC, s_plot, out_dir, "PSC_final", B_axis) +make_Bnormal_plots(bs + B_PSC, s_plot, out_dir, "PSC_and_TF_final", B_axis) +print('end') + diff --git a/examples/2_Intermediate/QA_wp_example.py b/examples/2_Intermediate/QA_wp_example.py new file mode 100644 index 000000000..36ec0a107 --- /dev/null +++ b/examples/2_Intermediate/QA_wp_example.py @@ -0,0 +1,344 @@ +#!/usr/bin/env python +r""" +This example script optimizes a set of relatively simple toroidal field coils +and passive superconducting coils (WPs) +for an ARIES-CS reactor-scale version of the precise-QH stellarator from +Landreman and Paul. + +The script should be run as: + mpirun -n 1 python QH_wp_example.py +on a cluster machine but + python QH_wp_example.py +is sufficient on other machines. Note that this code does not use MPI, but is +parallelized via OpenMP, so will run substantially +faster on multi-core machines (make sure that all the cores +are available to OpenMP, e.g. through setting OMP_NUM_THREADS). + +""" +from pathlib import Path + +import numpy as np +from scipy.optimize import minimize + +from simsopt.field import BiotSavart, Current, coils_via_symmetries +from simsopt.geo import SurfaceRZFourier, curves_to_vtk +from simsopt.geo.wp_grid import WPgrid +from simsopt.objectives import SquaredFlux +from simsopt.util import in_github_actions +from simsopt.util.permanent_magnet_helper_functions import * +import time + +def coil_optimization_QA(s, bs, base_curves, curves, out_dir=''): + """ + Optimize the coils for the QA, QH, or other configurations. + + Args: + s: plasma boundary. + bs: Biot Savart class object, presumably representing the + magnetic fields generated by the coils. + base_curves: List of CurveXYZFourier class objects. + curves: List of Curve class objects. + out_dir: Path or string for the output directory for saved files. + + Returns: + bs: Biot Savart class object, presumably representing the + OPTIMIZED magnetic fields generated by the coils. + """ + + from simsopt.geo import CurveLength, CurveCurveDistance, \ + MeanSquaredCurvature, LpCurveCurvature, CurveSurfaceDistance + from simsopt.objectives import QuadraticPenalty + from simsopt.geo import curves_to_vtk + from simsopt.objectives import SquaredFlux + + out_dir = Path(out_dir) + nphi = len(s.quadpoints_phi) + ntheta = len(s.quadpoints_theta) + ncoils = len(base_curves) + + # Weight on the curve lengths in the objective function: + LENGTH_WEIGHT = 1e-4 + + # Threshold and weight for the coil-to-coil distance penalty in the objective function: + CC_THRESHOLD = 0.2 + CC_WEIGHT = 1e1 + + # Threshold and weight for the coil-to-surface distance penalty in the objective function: + CS_THRESHOLD = 0.5 + CS_WEIGHT = 1 + + # Threshold and weight for the curvature penalty in the objective function: + CURVATURE_THRESHOLD = 10 + CURVATURE_WEIGHT = 1e-10 + + # Threshold and weight for the mean squared curvature penalty in the objective function: + MSC_THRESHOLD = 10 + MSC_WEIGHT = 1e-10 + + MAXITER = 500 # number of iterations for minimize + + # Define the objective function: + Jf = SquaredFlux(s, bs) + Jls = [CurveLength(c) for c in base_curves] + Jccdist = CurveCurveDistance(curves, CC_THRESHOLD, num_basecurves=ncoils) + Jcsdist = CurveSurfaceDistance(curves, s, CS_THRESHOLD) + Jcs = [LpCurveCurvature(c, 2, CURVATURE_THRESHOLD) for c in base_curves] + Jmscs = [MeanSquaredCurvature(c) for c in base_curves] + + # Form the total objective function. + JF = Jf \ + + LENGTH_WEIGHT * sum(Jls) \ + + CC_WEIGHT * Jccdist \ + + CS_WEIGHT * Jcsdist \ + + CURVATURE_WEIGHT * sum(Jcs) \ + + MSC_WEIGHT * sum(QuadraticPenalty(J, MSC_THRESHOLD) for J in Jmscs) + + def fun(dofs): + """ Function for coil optimization grabbed from stage_two_optimization.py """ + JF.x = dofs + J = JF.J() + grad = JF.dJ() + jf = Jf.J() + BdotN = np.mean(np.abs(np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2))) + outstr = f"J={J:.1e}, Jf={jf:.1e}, ⟨B·n⟩={BdotN:.1e}" + cl_string = ", ".join([f"{J.J():.1f}" for J in Jls]) + kap_string = ", ".join(f"{np.max(c.kappa()):.1f}" for c in base_curves) + msc_string = ", ".join(f"{J.J():.1f}" for J in Jmscs) + outstr += f", Len=sum([{cl_string}])={sum(J.J() for J in Jls):.1f}, ϰ=[{kap_string}], ∫ϰ²/L=[{msc_string}]" + outstr += f", C-C-Sep={Jccdist.shortest_distance():.2f}, C-S-Sep={Jcsdist.shortest_distance():.2f}" + outstr += f", ║∇J║={np.linalg.norm(grad):.1e}" + print(outstr) + return J, grad + + print(""" + ################################################################################ + ### Perform a Taylor test ###################################################### + ################################################################################ + """) + f = fun + dofs = JF.x + np.random.seed(1) + h = np.random.uniform(size=dofs.shape) + + J0, dJ0 = f(dofs) + dJh = sum(dJ0 * h) + for eps in [1e-3, 1e-4, 1e-5, 1e-6, 1e-7]: + J1, _ = f(dofs + eps*h) + J2, _ = f(dofs - eps*h) + print("err", (J1-J2)/(2*eps) - dJh) + + print(""" + ################################################################################ + ### Run the optimisation ####################################################### + ################################################################################ + """) + minimize(fun, dofs, jac=True, method='L-BFGS-B', options={'maxiter': MAXITER, 'maxcor': 300}, tol=1e-15) + curves_to_vtk(curves, out_dir / "curves_opt") + bs.set_points(s.gamma().reshape((-1, 3))) + return bs + + +np.random.seed(1) # set a seed so that the same WPs are initialized each time + +# Set some parameters -- if doing CI, lower the resolution +if in_github_actions: + nphi = 4 # nphi = ntheta >= 64 needed for accurate full-resolution runs + ntheta = nphi +else: + # Resolution needs to be reasonably high if you are doing permanent magnets + # or small coils because the fields are quite local + nphi = 32 # nphi = ntheta >= 64 needed for accurate full-resolution runs + ntheta = nphi + # Make higher resolution surface for plotting Bnormal + qphi = nphi * 2 + quadpoints_phi = np.linspace(0, 1, qphi, endpoint=True) + quadpoints_theta = np.linspace(0, 1, ntheta * 2, endpoint=True) + +poff = 0.2 # WP grid will be offset 'poff' meters from the plasma surface +coff = 0.3 # WP grid will be initialized between 1 m and 2 m from plasma + +# Read in the plasma equilibrium file +input_name = 'input.LandremanPaul2021_QA_lowres' +TEST_DIR = (Path(__file__).parent / ".." / ".." / "tests" / "test_files").resolve() +surface_filename = TEST_DIR / input_name +range_param = 'half period' +s = SurfaceRZFourier.from_vmec_input( + surface_filename, range=range_param, nphi=nphi, ntheta=ntheta +) +# Print major and minor radius +print('s.R = ', s.get_rc(0, 0)) +print('s.r = ', s.get_rc(1, 0)) + +# Make inner and outer toroidal surfaces very high resolution, +# which helps to initialize coils precisely between the surfaces. +s_inner = SurfaceRZFourier.from_vmec_input( + surface_filename, range=range_param, nphi=nphi * 4, ntheta=ntheta * 4 +) +s_outer = SurfaceRZFourier.from_vmec_input( + surface_filename, range=range_param, nphi=nphi * 4, ntheta=ntheta * 4 +) + +# Make the inner and outer surfaces by extending the plasma surface +s_inner.extend_via_normal(poff) +s_outer.extend_via_normal(poff + coff) + +# Make the output directory +out_str = "QA_wp_output/" +out_dir = Path("QA_wp_output") +out_dir.mkdir(parents=True, exist_ok=True) + +# Save the inner and outer surfaces for debugging purposes +s_inner.to_vtk(out_str + 'inner_surf') +s_outer.to_vtk(out_str + 'outer_surf') + +# initialize the coils +base_curves, curves, coils = initialize_coils('qa_psc', TEST_DIR, s, out_dir) +currents = np.array([coil.current.get_value() for coil in coils]) +print('Currents = ', currents) + +# Set up BiotSavart fields +bs = BiotSavart(coils) + +# Calculate average, approximate on-axis B field strength +calculate_on_axis_B(bs, s) + +# Make high resolution, full torus version of the plasma boundary for plotting +s_plot = SurfaceRZFourier.from_vmec_input( + surface_filename, + quadpoints_phi=quadpoints_phi, + quadpoints_theta=quadpoints_theta +) + +# Plot initial Bnormal on plasma surface from un-optimized BiotSavart coils +make_Bnormal_plots(bs, s_plot, out_dir, "biot_savart_initial") + +# optimize the currents in the TF coils and plot results +# fix all the coil shapes so only the currents are optimized +# for i in range(ncoils): +# base_curves[i].fix_all() +# bs = coil_optimization_QA(s, bs, base_curves, curves, out_dir) +curves_to_vtk(curves, out_dir / "TF_coils", close=True) +bs.set_points(s.gamma().reshape((-1, 3))) +Bnormal = np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2) +make_Bnormal_plots(bs, s_plot, out_dir, "biot_savart_TF_optimized") + +# check after-optimization average on-axis magnetic field strength +B_axis = calculate_on_axis_B(bs, s) +make_Bnormal_plots(bs, s_plot, out_dir, "biot_savart_TF_optimized", B_axis) + +# Finally, initialize the wp class +kwargs_geo = {"Nx": 4, "out_dir": out_str, + "random_initialization": True, "poff": poff,} + # "interpolated_field": True} +wp_array = WPgrid.geo_setup_between_toroidal_surfaces( + s, coils, s_inner, s_outer, **kwargs_geo +) +print('Number of WP locations = ', len(wp_array.grid_xyz)) + +currents = [] +for i in range(wp_array.num_wp): + currents.append(Current(wp_array.I[i])) +all_coils = coils_via_symmetries( + wp_array.curves, currents, nfp=wp_array.nfp, stellsym=wp_array.stellsym +) +B_WP = BiotSavart(all_coils) +# Plot initial errors from only the WPs, and then together with the TF coils +make_Bnormal_plots(B_WP, s_plot, out_dir, "biot_savart_WP_initial", B_axis) +make_Bnormal_plots(bs + B_WP, s_plot, out_dir, "WP_and_TF_initial", B_axis) + +# Check SquaredFlux values using different ways to calculate it +x0 = np.ravel(np.array([wp_array.alphas, wp_array.deltas, wp_array.I])) +fB = SquaredFlux(s, bs, np.zeros((nphi, ntheta))).J() +print('fB only TF coils = ', fB / (B_axis ** 2 * s.area())) +bs.set_points(s.gamma().reshape(-1, 3)) +Bnormal = np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2) +print('fB only TF direct = ', np.sum(Bnormal.reshape(-1) ** 2 * wp_array.grid_normalization ** 2 + ) / (2 * B_axis ** 2 * s.area())) +make_Bnormal_plots(B_WP, s_plot, out_dir, "biot_savart_WP_initial", B_axis) +fB = SquaredFlux(s, B_WP, np.zeros((nphi, ntheta))).J() +print(fB/ (B_axis ** 2 * s.area())) +fB = SquaredFlux(s, B_WP + bs, np.zeros((nphi, ntheta))).J() +print('fB with both, before opt = ', fB / (B_axis ** 2 * s.area())) +fB = SquaredFlux(s, B_WP, -Bnormal).J() +print('fB with both (minus sign), before opt = ', fB / (B_axis ** 2 * s.area())) +Bnormal_WP = np.sum(B_WP.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2) +fB = SquaredFlux(s, bs, -Bnormal_WP).J() +print('fB with both, before opt = ', fB / (B_axis ** 2 * s.area())) + +# Actually do the minimization now +print('beginning optimization: ') +opt_bounds = [(0, 2 * np.pi) for i in range(wp_array.num_wp * 2)] +for i in range(wp_array.num_wp): + opt_bounds.append((None, None)) +opt_bounds = tuple(opt_bounds) +options = {"disp": True, "maxiter": 200, "eps": 1e-4} +# print(opt_bounds) +# x0 = np.random.rand(2 * wp_array.num_wp) * 2 * np.pi +verbose = True + +# Run STLSQ with BFGS in the loop +kwargs_manual = { + "out_dir": out_str, + "plasma_boundary" : s, + "coils_TF" : coils + } +# I_threshold = 0.0 +# STLSQ_max_iters = 10 +# for k in range(STLSQ_max_iters): +# x0 = np.ravel(np.array([wp_array.alphas, wp_array.deltas, wp_array.I])) +# print('Number of WPs = ', len(x0) // 4, ' in iteration ', k) + +# x0 = np.ravel(np.array([wp_array.I])) +x0 = np.ravel(np.array([wp_array.alphas, wp_array.I])) +x_opt = minimize(wp_array.least_squares, x0, args=(verbose,), + method='L-BFGS-B', + # method='BFGS', + # bounds=opt_bounds, + jac=wp_array.least_squares_jacobian, + options=options, + tol=1e-20, + ) +# print(x_opt.message) +# print(x_opt.jac) +# print(x_opt.hess_inv) + # I = wp_array.I + # small_I_inds = np.ravel(np.where(np.abs(I) < I_threshold)) + # grid_xyz = wp_array.grid_xyz + # alphas = wp_array.alphas + # deltas = wp_array.deltas + # if len(small_I_inds) > 0: + # grid_xyz = grid_xyz[small_I_inds, :] + # alphas = alphas[small_I_inds] + # deltas = deltas[small_I_inds] + # else: + # print('STLSQ converged, breaking out of loop') + # break + # kwargs_manual["alphas"] = alphas + # kwargs_manual["deltas"] = deltas + # Initialize new WP array with coils only at the remaining locations + # with initial orientations from the solve using BFGS + # wp_array = WPgrid.geo_setup_manual( + # grid_xyz, wp_array.R, **kwargs_manual + # ) + +wp_array.I = x_opt.x +# ind3 = len(x_opt.x) // 3 +# wp_array.I = x_opt.x[2 * ind3:] +# wp_array.setup_orientations(x_opt.x[:ind3], x_opt.x[ind3:2 * ind3]) +wp_array.setup_curves() +wp_array.plot_curves('final_') +currents = [] +for i in range(wp_array.num_wp): + currents.append(Current(wp_array.I[i])) +all_coils = coils_via_symmetries( + wp_array.curves, currents, nfp=wp_array.nfp, stellsym=wp_array.stellsym +) +B_WP = BiotSavart(all_coils) + +# Check that direct Bn calculation agrees with optimization calculation +fB = SquaredFlux(s, B_WP + bs, np.zeros((nphi, ntheta))).J() +print('fB with both, after opt = ', fB / (B_axis ** 2 * s.area())) +make_Bnormal_plots(B_WP, s_plot, out_dir, "WP_final", B_axis) +make_Bnormal_plots(bs + B_WP, s_plot, out_dir, "WP_and_TF_final", B_axis) +print('end') + diff --git a/examples/2_Intermediate/QH_psc_example.py b/examples/2_Intermediate/QH_psc_example.py index 05104ecda..1fb6308e9 100644 --- a/examples/2_Intermediate/QH_psc_example.py +++ b/examples/2_Intermediate/QH_psc_example.py @@ -44,7 +44,7 @@ quadpoints_theta = np.linspace(0, 1, ntheta * 4, endpoint=True) poff = 1.0 # PSC grid will be offset 'poff' meters from the plasma surface -coff = 1.2 # PSC grid will be initialized between 1 m and 2 m from plasma +coff = 0.5 # PSC grid will be initialized between 1 m and 2 m from plasma # Read in the plasma equilibrium file input_name = 'input.LandremanPaul2021_QH_reactorScale_lowres' @@ -115,7 +115,7 @@ make_Bnormal_plots(bs, s_plot, out_dir, "biot_savart_TF_optimized", B_axis) # Finally, initialize the psc class -kwargs_geo = {"Nx": 10, "out_dir": out_str, "random_initialization": True, "poff": poff} +kwargs_geo = {"Nx": 10, "out_dir": out_str, "random_initialization": True} psc_array = PSCgrid.geo_setup_between_toroidal_surfaces( s, coils, s_inner, s_outer, **kwargs_geo ) @@ -136,7 +136,6 @@ x0 = np.ravel(np.array([psc_array.alphas, psc_array.deltas])) fB = SquaredFlux(s, bs, np.zeros((nphi, ntheta))).J() print('fB only TF coils = ', fB / (B_axis ** 2 * s.area())) -psc_array.least_squares(np.zeros(x0.shape)) bs.set_points(s.gamma().reshape(-1, 3)) Bnormal = np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2) print('fB only TF direct = ', np.sum(Bnormal.reshape(-1) ** 2 * psc_array.grid_normalization ** 2 @@ -152,8 +151,8 @@ # Actually do the minimization now from scipy.optimize import minimize print('beginning optimization: ') -# opt_bounds = tuple([(0, 2 * np.pi) for i in range(psc_array.num_psc * 2)]) -# options = {"disp": True} #, "maxiter": 100} +opt_bounds = tuple([(0, 2 * np.pi) for i in range(psc_array.num_psc * 2)]) +options = {"disp": True, "maxiter": 100} # print(opt_bounds) # x0 = np.random.rand(2 * psc_array.num_psc) * 2 * np.pi verbose = True @@ -164,15 +163,15 @@ "plasma_boundary" : s, "coils_TF" : coils } -I_threshold = 1e4 -STLSQ_max_iters = 2 +I_threshold = 0.0 +STLSQ_max_iters = 10 for k in range(STLSQ_max_iters): + # x0 = np.ravel(np.array([psc_array.alphas, psc_array.deltas])) print('Number of PSCs = ', len(x0) // 2, ' in iteration ', k) - x0 = np.ravel(np.array([psc_array.alphas, psc_array.deltas])) x_opt = minimize(psc_array.least_squares, x0, args=(verbose,), method='L-BFGS-B', # bounds=opt_bounds, - jac=psc_array.least_squares_jacobian, + # jac=psc_array.least_squares_jacobian, options=options, tol=1e-10, ) @@ -185,6 +184,9 @@ grid_xyz = grid_xyz[small_I_inds, :] alphas = alphas[small_I_inds] deltas = deltas[small_I_inds] + else: + print('STLSQ converged, breaking out of loop') + break kwargs_manual["alphas"] = alphas kwargs_manual["deltas"] = deltas # Initialize new PSC array with coils only at the remaining locations diff --git a/src/simsopt/geo/psc_grid.py b/src/simsopt/geo/psc_grid.py index efdd30274..2972f1c51 100644 --- a/src/simsopt/geo/psc_grid.py +++ b/src/simsopt/geo/psc_grid.py @@ -28,6 +28,7 @@ class PSCgrid: def __init__(self): self.mu0 = 4 * np.pi * 1e-7 + self.fac = 1e-7 # Define a set of quadrature points and weights for the N point # Gaussian quadrature rule num_quad = 10 @@ -69,8 +70,8 @@ def _setup_uniform_grid(self): # This is not a guarantee that coils will not touch but inductance # matrix blows up if they do so it is easy to tell when they do - self.R = Nmin / 4.0 - self.a = self.R / 100.0 # Hard-coded aspect ratio of 100 right now + self.R = Nmin / 4.0 # , self.poff / 4.0) + self.a = self.R / 10.0 # Hard-coded aspect ratio of 100 right now print('Major radius of the coils is R = ', self.R) print('Coils are spaced so that every coil of radius R ' ' is at least 2R away from the next coil' @@ -106,10 +107,15 @@ def _setup_uniform_grid(self): for j in range(Ny): for k in range(Nz): phi = np.arctan2(Y[i, j, k], X[i, j, k]) - phi2 = np.arctan2(self.R, X[i, j, k]) + if self.nfp == 4: + phi2 = np.arctan2(self.R, X[i, j, k]) + elif self.nfp == 2: + phi2 = np.arctan2(self.R, self.plasma_boundary.get_rc(0, 0)) # Add a little factor to avoid phi = pi / n_p degrees # exactly, which can intersect with a symmetrized # coil if not careful + # print(phi2) + # exit() if phi >= (np.pi / self.nfp - phi2) or phi < 0.0: inds.append(int(i * Ny * Nz + j * Nz + k)) good_inds = np.setdiff1d(np.arange(Nx * Ny * Nz), inds) @@ -622,7 +628,7 @@ def setup_currents_and_fields(self): for stell in self.stell_list: xyz = self.grid_xyz_all[q * nn: (q + 1) * nn, :] t1 = time.time() - A_matrix += sopp.A_matrix( + A_matrix += 2 * sopp.A_matrix( contig(xyz), contig(self.plasma_points), contig(self.alphas_total[q * nn: (q + 1) * nn]), @@ -640,7 +646,7 @@ def setup_currents_and_fields(self): # ) q = q + 1 self.A_matrix = A_matrix - self.Bn_PSC = (A_matrix @ self.I).reshape(-1) * self.grid_normalization + self.Bn_PSC = (A_matrix @ self.I).reshape(-1) # self.B_PSC = B_PSC # currents = [] # for i in range(self.num_psc): @@ -780,7 +786,7 @@ def b_vector(self): Bn_TF = np.sum( self.B_TF.B().reshape(-1, 3) * self.plasma_unitnormals, axis=-1 ) - self.b_opt = (Bn_TF + Bn_plasma) * self.grid_normalization + self.b_opt = (Bn_TF + Bn_plasma) / self.fac def least_squares(self, kappas, verbose=False): """ @@ -806,14 +812,14 @@ def least_squares(self, kappas, verbose=False): alphas = kappas[:len(kappas) // 2] deltas = kappas[len(kappas) // 2:] self.setup_orientations(alphas, deltas) - Ax_b = self.Bn_PSC + self.b_opt - BdotN2 = 0.5 * Ax_b.T @ Ax_b / self.normalization - outstr = f"Normalized f_B = {BdotN2:.3e} " + Ax_b = (self.Bn_PSC + self.b_opt) * self.grid_normalization + BdotN2 = 0.5 * Ax_b.T @ Ax_b / self.normalization * self.fac ** 2 + # outstr = f"Normalized f_B = {BdotN2:.3e} " # for i in range(len(kappas)): - # outstr += f"kappas[{i:d}] = {kappas[i]:.2e} " + # outstr += f"kappas[{i:d}] = {kappas[i]:.2e} " # outstr += "\n" - if verbose: - print(outstr) + # if verbose: + # print(outstr) return BdotN2 def least_squares_jacobian(self, kappas, verbose=False): @@ -839,9 +845,7 @@ def least_squares_jacobian(self, kappas, verbose=False): deltas = kappas[len(kappas) // 2:] self.setup_orientations(alphas, deltas) # Two factors of grid normalization since it is not added to the gradients - # A has shape (num_plasma_points) - Ax_b = (self.Bn_PSC + self.b_opt) / self.normalization - # Ax_b = (self.Bn_PSC * self.grid_normalization + self.b_opt) / self.normalization + Ax_b = (self.Bn_PSC + self.b_opt) / self.normalization * self.grid_normalization # t1 = time.time() A_deriv = self.A_deriv() Linv = self.L_inv[:self.num_psc, :self.num_psc] @@ -915,7 +919,7 @@ def L_deriv(self): contig(self.deltas), contig(self.quad_points_phi), contig(self.quad_weights_phi) - ) / (4 * np.pi) # * self.quad_dphi ** 2 / (np.pi ** 2) + ) / (4 * np.pi) # symmetrize it L_deriv = (L_deriv + np.transpose(L_deriv, axes=[0, 2, 1])) diff --git a/src/simsopt/util/permanent_magnet_helper_functions.py b/src/simsopt/util/permanent_magnet_helper_functions.py index c06c383de..10bb0df51 100644 --- a/src/simsopt/util/permanent_magnet_helper_functions.py +++ b/src/simsopt/util/permanent_magnet_helper_functions.py @@ -117,7 +117,7 @@ def coil_optimization(s, bs, base_curves, curves, out_dir=''): CC_WEIGHT = 1 # Threshold and weight for the coil-to-surface distance penalty in the objective function: - CS_THRESHOLD = 0.1 + CS_THRESHOLD = 1.5 CS_WEIGHT = 1e-2 # Threshold and weight for the curvature penalty in the objective function: @@ -305,7 +305,7 @@ def initialize_coils(config_flag, TEST_DIR, s, out_dir=''): base_curves[i].fix_all() elif config_flag == 'qh': # generate planar TF coils - ncoils = 2 + ncoils = 4 R0 = s.get_rc(0, 0) R1 = s.get_rc(1, 0) * 4 order = 2 @@ -317,7 +317,10 @@ def initialize_coils(config_flag, TEST_DIR, s, out_dir=''): print('Total current = ', total_current) base_curves = create_equally_spaced_curves(ncoils, s.nfp, stellsym=True, R0=R0, R1=R1, order=order, numquadpoints=64) - base_currents = [(Current(total_current / ncoils * 1e-5) * 1e5) for _ in range(ncoils-1)] + base_currents = [(Current(total_current / ncoils * 1e-5) * 1e5) for _ in range(ncoils - 1)] + # base_currents = [(Current(total_current / ncoils * 1e-5) * 1e5) for _ in range(ncoils)] + # base_currents[0].fix_all() + total_current = Current(total_current) total_current.fix_all() base_currents += [total_current - sum(base_currents)] @@ -346,9 +349,9 @@ def initialize_coils(config_flag, TEST_DIR, s, out_dir=''): elif config_flag == 'qa_psc': # generate planar TF coils - ncoils = 6 + ncoils = 2 R0 = 1.0 - R1 = 0.65 + R1 = 0.9 order = 5 # qa needs to be scaled to 0.1 T on-axis magnetic field strength @@ -362,8 +365,8 @@ def initialize_coils(config_flag, TEST_DIR, s, out_dir=''): base_currents += [total_current - sum(base_currents)] coils = coils_via_symmetries(base_curves, base_currents, s.nfp, True) # fix all the coil shapes so only the currents are optimized - for i in range(ncoils): - base_curves[i].fix_all() + # for i in range(ncoils): + # base_curves[i].fix_all() # Initialize the coil curves and save the data to vtk curves = [c.curve for c in coils] diff --git a/src/simsoptpp/psc.cpp b/src/simsoptpp/psc.cpp index 08295ed68..4bb836d4f 100644 --- a/src/simsoptpp/psc.cpp +++ b/src/simsoptpp/psc.cpp @@ -867,7 +867,7 @@ Array A_matrix(Array& points, Array& plasma_points, Array& alphas, Array& deltas A(i, j) = Bx_rot * nx + By_rot * ny + Bz_rot * nz; } } - return A * fac; + return A; } Array B_PSC(Array& points, Array& plasma_points, Array& alphas, Array& deltas, Array& psc_currents, double R) @@ -983,7 +983,6 @@ Array A_matrix_direct(Array& points, Array& plasma_points, Array& alphas, Array& // this variable is the A matrix in the least-squares term so A * I = Bn Array A = xt::zeros({num_plasma_points, num_coils}); - double fac = 1.0e-7; using namespace boost::math; #pragma omp parallel for schedule(static) @@ -1041,7 +1040,7 @@ Array A_matrix_direct(Array& points, Array& plasma_points, Array& alphas, Array& A(i, j) = (Rxx * Bx + Rxy * By + Rxz * Bz) * nx + (Ryx * Bx + Ryy * By + Ryz * Bz) * ny + (Rzx * Bx + Rzy * By + Rzz * Bz) * nz; } } - return A * fac; + return A; } @@ -1069,7 +1068,6 @@ Array dA_dkappa(Array& points, Array& plasma_points, Array& alphas, Array& delta // this variable is the A matrix in the least-squares term so A * I = Bn Array dA = xt::zeros({num_coils * 2, num_plasma_points}); - double fac = 1.0e-7; using namespace boost::math; #pragma omp parallel for schedule(static) @@ -1197,7 +1195,7 @@ Array dA_dkappa(Array& points, Array& plasma_points, Array& alphas, Array& delta dA(j + num_coils, i) += (Rzx * (Bx4 + Bx5) + Rzy * (By4 + By5) + Rzz * (Bz4 + Bz5)) * nz; } } - return dA * fac; + return dA; } diff --git a/src/simsoptpp/psc.h b/src/simsoptpp/psc.h index de4de07ae..4a4904d5f 100644 --- a/src/simsoptpp/psc.h +++ b/src/simsoptpp/psc.h @@ -10,7 +10,6 @@ #include // c++ tuples #include // for string class #include -#include #include "xtensor-python/pyarray.hpp" // Numpy bindings typedef xt::pyarray Array; diff --git a/tests/geo/test_psc_grid.py b/tests/geo/test_psc_grid.py index aac7dc8e9..083e315b2 100644 --- a/tests/geo/test_psc_grid.py +++ b/tests/geo/test_psc_grid.py @@ -31,7 +31,7 @@ def test_dpsi_analytic_derivatives(self): ) A = psc_array.A_matrix psi = psc_array.psi - b = psc_array.b_opt / psc_array.grid_normalization + b = psc_array.b_opt # / psc_array.grid_normalization L = psc_array.L Linv = psc_array.L_inv Bn_objective = 0.5 * (A @ Linv @ psi + b).T @ (A @ Linv @ psi + b) @@ -61,7 +61,7 @@ def test_dpsi_analytic_derivatives(self): ) A = psc_array.A_matrix psi = psc_array.psi - b = psc_array.b_opt / psc_array.grid_normalization + b = psc_array.b_opt # / psc_array.grid_normalization Linv = psc_array.L_inv Bn_objective = 0.5 * (A @ Linv @ psi + b).T @ (A @ Linv @ psi + b) psi_deriv = psc_array.psi_deriv() @@ -90,7 +90,7 @@ def test_dpsi_analytic_derivatives(self): ) A = psc_array.A_matrix psi = psc_array.psi - b = psc_array.b_opt / psc_array.grid_normalization + b = psc_array.b_opt # / psc_array.grid_normalization Linv = psc_array.L_inv Bn_objective = 0.5 * (A @ Linv @ psi + b).T @ (A @ Linv @ psi + b) psi_deriv = psc_array.psi_deriv() @@ -119,7 +119,7 @@ def test_dpsi_analytic_derivatives(self): ) A = psc_array.A_matrix psi = psc_array.psi - b = psc_array.b_opt / psc_array.grid_normalization + b = psc_array.b_opt # / psc_array.grid_normalization Linv = psc_array.L_inv Bn_objective = 0.5 * (A @ Linv @ psi + b).T @ (A @ Linv @ psi + b) psi_deriv = psc_array.psi_deriv() @@ -160,7 +160,7 @@ def test_L_analytic_derivatives(self): ) A = psc_array.A_matrix psi = psc_array.psi - b = psc_array.b_opt / psc_array.grid_normalization + b = psc_array.b_opt # / psc_array.grid_normalization L = psc_array.L Linv = psc_array.L_inv Bn_objective = 0.5 * (A @ Linv @ psi + b).T @ (A @ Linv @ psi + b) @@ -199,7 +199,7 @@ def test_L_analytic_derivatives(self): ) A = psc_array.A_matrix psi = psc_array.psi - b = psc_array.b_opt / psc_array.grid_normalization + b = psc_array.b_opt # / psc_array.grid_normalization L = psc_array.L Linv = psc_array.L_inv Bn_objective = 0.5 * (A @ Linv @ psi + b).T @ (A @ Linv @ psi + b) @@ -240,7 +240,7 @@ def test_L_analytic_derivatives(self): ) A = psc_array.A_matrix psi = psc_array.psi - b = psc_array.b_opt / psc_array.grid_normalization + b = psc_array.b_opt # / psc_array.grid_normalization L = psc_array.L Linv = psc_array.L_inv Bn_objective = 0.5 * (A @ Linv @ psi + b).T @ (A @ Linv @ psi + b) @@ -281,7 +281,7 @@ def test_L_analytic_derivatives(self): ) A = psc_array.A_matrix psi = psc_array.psi - b = psc_array.b_opt / psc_array.grid_normalization + b = psc_array.b_opt # / psc_array.grid_normalization L = psc_array.L Linv = psc_array.L_inv Bn_objective = 0.5 * (A @ Linv @ psi + b).T @ (A @ Linv @ psi + b) @@ -335,9 +335,9 @@ def test_dA_analytic_derivatives(self): Linv = psc_array.L_inv[:ncoils, :ncoils] psi = psc_array.psi Linv_psi = Linv @ psi - b = psc_array.b_opt / psc_array.grid_normalization + b = psc_array.b_opt # / psc_array.grid_normalization Bn_objective = 0.5 * (A @ Linv_psi + b).T @ (A @ Linv_psi + b) - epsilon = 1e-4 + epsilon = 1e-6 alphas_new = alphas alphas_new[0] += epsilon deltas_new = deltas @@ -352,7 +352,7 @@ def test_dA_analytic_derivatives(self): assert np.isclose(dBn_objective, dBn_analytic[0], rtol=1e-2) dA_dalpha = (A_new - A) / epsilon dA_dkappa_analytic = A_deriv - assert np.allclose(dA_dalpha[:, 0], dA_dkappa_analytic[:, 0]) + assert np.allclose(dA_dalpha[:, 0], dA_dkappa_analytic[:, 0], rtol=1e-2) # Repeat but change coil 3 points = (np.random.rand(ncoils, 3) - 0.5) * 5 @@ -365,7 +365,7 @@ def test_dA_analytic_derivatives(self): Linv = psc_array.L_inv[:ncoils, :ncoils] psi = psc_array.psi Linv_psi = Linv @ psi - b = psc_array.b_opt / psc_array.grid_normalization + b = psc_array.b_opt # / psc_array.grid_normalization Bn_objective = 0.5 * (A @ Linv_psi + b).T @ (A @ Linv_psi + b) epsilon = 1e-4 alphas_new = alphas @@ -382,7 +382,7 @@ def test_dA_analytic_derivatives(self): assert np.isclose(dBn_objective, dBn_analytic[4], rtol=1e-2) dA_dalpha = (A_new - A) / epsilon dA_dkappa_analytic = A_deriv - assert np.allclose(dA_dalpha[:, 4], dA_dkappa_analytic[:, 4]) + assert np.allclose(dA_dalpha[:, 4], dA_dkappa_analytic[:, 4], rtol=1e-2) # Repeat but changing delta instead of alpha deltas = (np.random.rand(ncoils) - 0.5) * 2 * np.pi @@ -394,7 +394,7 @@ def test_dA_analytic_derivatives(self): Linv = psc_array.L_inv[:ncoils, :ncoils] psi = psc_array.psi Linv_psi = Linv @ psi - b = psc_array.b_opt / psc_array.grid_normalization + b = psc_array.b_opt # / psc_array.grid_normalization Bn_objective = 0.5 * (A @ Linv_psi + b).T @ (A @ Linv_psi + b) epsilon = 1e-4 alphas_new = alphas @@ -409,6 +409,7 @@ def test_dA_analytic_derivatives(self): A_deriv = psc_array.A_deriv() # notice ncoils: instead of :ncoils below dBn_analytic = (A @ Linv_psi + b).T @ (A_deriv[:, ncoils:] * Linv_psi) + print(dBn_objective, dBn_analytic[0]) assert np.isclose(dBn_objective, dBn_analytic[0], rtol=1e-2) dA_ddelta = (A_new - A) / epsilon dA_dkappa_analytic = A_deriv @@ -424,7 +425,7 @@ def test_dA_analytic_derivatives(self): Linv = psc_array.L_inv[:ncoils, :ncoils] psi = psc_array.psi Linv_psi = Linv @ psi - b = psc_array.b_opt / psc_array.grid_normalization + b = psc_array.b_opt # / psc_array.grid_normalization Bn_objective = 0.5 * (A @ Linv_psi + b).T @ (A @ Linv_psi + b) epsilon = 1e-4 alphas_new = alphas @@ -589,7 +590,7 @@ def test_L(self): contig(psc_array.deltas_total[q * nn: (q + 1) * nn]), contig(psc_array.plasma_boundary.unitnormal().reshape(-1, 3)), psc_array.R, - ) @ psc_array.I + ) @ psc_array.I * 2e-7 B_PSC += sopp.B_PSC( contig(xyz), contig(psc_array.plasma_boundary.gamma().reshape(-1, 3)), @@ -636,6 +637,7 @@ def test_L(self): Bn_direct_all = np.sum(B_direct_all.B().reshape( -1, 3) * psc_array.plasma_boundary.unitnormal().reshape(-1, 3), axis=-1) + psc_array.Bn_PSC = psc_array.Bn_PSC # / psc_array.grid_normalization # Robust test of all the B and Bn calculations from circular coils assert(np.allclose(psc_array.Bn_PSC * 1e10, Bn_PSC * 1e10, atol=1e3)) assert(np.allclose(psc_array.Bn_PSC * 1e10, Bn_circular_coils * 1e10, atol=1e3)) @@ -681,7 +683,7 @@ def test_L(self): contig(psc_array.deltas_total[q * nn: (q + 1) * nn]), contig(psc_array.plasma_boundary.unitnormal().reshape(-1, 3)), psc_array.R, - ) @ (psc_array.I) + ) @ (psc_array.I) * 2e-7 B_PSC += sopp.B_PSC( contig(xyz), contig(psc_array.plasma_points), @@ -737,6 +739,7 @@ def test_L(self): -1, 3) * psc_array.plasma_boundary.unitnormal().reshape(-1, 3), axis=-1) # Robust test of all the B and Bn calculations from circular coils + psc_array.Bn_PSC = psc_array.Bn_PSC #/ psc_array.grid_normalization assert(np.allclose(psc_array.Bn_PSC * 1e10, Bn_PSC * 1e10, rtol=1e-2, atol=1e3)) assert(np.allclose(psc_array.Bn_PSC * 1e10, Bn_circular_coils * 1e10, rtol=1e-2, atol=1e3)) assert(np.allclose(psc_array.Bn_PSC * 1e10, Bn_direct * 1e10, rtol=1e-2, atol=1e3))