Skip to content

Commit

Permalink
debugging objective calculation
Browse files Browse the repository at this point in the history
  • Loading branch information
@WillhHoffman
WillhHoffman committed Jan 31, 2025
1 parent 268ded8 commit 1ba3079
Showing 1 changed file with 70 additions and 36 deletions.
106 changes: 70 additions & 36 deletions examples/2_Intermediate/dipole_QA.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
ntheta = nphi
dr = 0.05 # cylindrical bricks with radial extent 5 cm
else:
nphi = 32 # nphi = ntheta >= 64 needed for accurate full-resolution runs
nphi = 64 # nphi = ntheta >= 64 needed for accurate full-resolution runs
ntheta = nphi
# dr = 0.02 # cylindrical bricks with radial extent 2 cm
Nx = 64
Expand All @@ -30,7 +30,7 @@
poff = 0.03 # PM grid end offset ~ 15 cm from the plasma surface
input_name = 'input.LandremanPaul2021_QA_lowres'

nIter_max = 5000
nIter_max = 1000

# Read in the plas/ma equilibrium file
TEST_DIR = (Path(__file__).parent / ".." / ".." / "tests" / "test_files").resolve()
Expand Down Expand Up @@ -127,8 +127,8 @@
# Print optimized metrics
print("Total fB = ",
0.5 * np.sum((pm_opt.A_obj @ pm_opt.m - pm_opt.b_obj) ** 2))
print("Total fB (sparse) = ",
0.5 * np.sum((pm_opt.A_obj @ pm_opt.m_proxy - pm_opt.b_obj) ** 2))
# print("Total fB (sparse) = ",a
# 0.5 * np.sum((pm_opt.A_obj @ pm_opt.m_proxy - pm_opt.b_obj) ** 2))

bs.set_points(s_plot.gamma().reshape((-1, 3)))
Bnormal = np.sum(bs.B().reshape((qphi, ntheta, 3)) * s_plot.unitnormal(), axis=2)
Expand All @@ -151,13 +151,8 @@

# Print optimized f_B and other metrics
print('B_field shape = ',b_dipole.B().shape)
# print('B_field = ',b_dipole.B())
f_B_sf = SquaredFlux(s_plot, b_dipole, -Bnormal).J()

# print('s_plot = ',s_plot)
# print('b_dipole = ',b_dipole)
# print('- Bnorm = ',-Bnormal)

print('f_B = ', f_B_sf)
total_volume = np.sum(np.sqrt(np.sum(pm_opt.m.reshape(pm_opt.ndipoles, 3) ** 2, axis=-1))) * s.nfp * 2 * mu0 / B_max
print('Total volume = ', total_volume)
Expand All @@ -176,40 +171,79 @@
m_maxima = pm_opt.m_maxima
)

b_comp.set_points(s_plot.gamma().reshape((-1, 3)))
b_comp._toVTK(out_dir / "comp_Fields", pm_opt.dx, pm_opt.dy, pm_opt.dz)
assert pm_comp.dx == pm_opt.dx
assert pm_comp.dy == pm_opt.dy
assert pm_comp.dz == pm_opt.dz

print('exact total fB = ', 0.5 * np.sum((pm_comp.A_obj @ pm_comp.m - pm_comp.b_obj) ** 2))
b_comp.set_points(s_plot.gamma().reshape((-1, 3)))
b_comp._toVTK(out_dir / "magnet_fields", pm_comp.dx, pm_comp.dy, pm_comp.dz)

# Print optimized metrics
# print("Total test fB = ",
# 0.5 * np.sum((pm_opt.A_obj @ pm_opt.m - pm_opt.b_obj) ** 2))
# in order to get a correct fB, have to have an exact grid that compows the dipole optimzation, placing magnets
# in all the same positions. Actually, do I need this even for the correct field?
print("comp fB = ",
0.5 * np.sum((pm_comp.A_obj @ pm_opt.m - pm_opt.b_obj) ** 2))

# For plotting Bn on the full torus surface at the end with just the dipole fields
make_Bnormal_plots(b_comp, s_plot, out_dir, "only_m__comp_optimized")
bs.set_points(s_plot.gamma().reshape((-1, 3)))
Bcnormal = np.sum(bs.B().reshape((qphi, ntheta, 3)) * s_plot.unitnormal(), axis=2)
make_Bnormal_plots(bs, s_plot, out_dir, "biot_savart_optimized")
Bcnormal_magnets = np.sum(b_comp.B().reshape((qphi, ntheta, 3)) * s_plot.unitnormal(), axis=-1)
Bcnormal_total = Bcnormal + Bcnormal_magnets

# Compute metrics with permanent magnet results
magnets_m = pm_opt.m.reshape(pm_opt.ndipoles, 3)
num_nonzero = np.count_nonzero(np.sum(magnets_m ** 2, axis=-1)) / pm_opt.ndipoles * 100
print("Number of possible magnets = ", pm_opt.ndipoles)
print("% of magnets that are nonzero = ", num_nonzero)

# For plotting Bn on the full torus surface at the end with just the magnet fields
make_Bnormal_plots(b_comp, s_plot, out_dir, "only_m_comp_optimized")
pointData = {"B_N": Bnormal_total[:, :, None]}
s_plot.to_vtk(out_dir / "m_comp_optimized", extra_data=pointData)

# Print optimized f_B and other metrics
print('B_Field_comp shape = ',b_comp.B().shape)
# print('B_Field_comp = ',b_comp.B())
f_Btest_sf = SquaredFlux(s_plot, b_comp, -Bnormal).J()

print('b_comp = ',b_comp)

print('f_B_comp = ', f_Btest_sf)

# b_dipole = DipoleField(
# pm_opt.dipole_grid_xyz,
# pm_opt.m,
# nfp=s.nfp,
# coordinate_flag=pm_opt.coordinate_flag,
# m_maxima=pm_opt.m_maxima
# )
# b_dipole.set_points(s_plot.gamma().reshape((-1, 3)))
# dipoles = pm_opt.m.reshape(pm_opt.ndipoles, 3)
# num_nonzero_sparse = np.count_nonzero(np.sum(dipoles ** 2, axis=-1)) / pm_opt.ndipoles * 100
b_comp.set_points(s.gamma().reshape((-1, 3)))
bs.set_points(s.gamma().reshape((-1, 3)))
Bnormal = np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2)
f_B_sf = SquaredFlux(s, b_comp, -Bnormal).J()

print('f_B = ', f_B_sf)
total_volume = np.sum(np.sqrt(np.sum(pm_opt.m.reshape(pm_opt.ndipoles, 3) ** 2, axis=-1))) * s.nfp * 2 * mu0 / B_max
print('Total volume = ', total_volume)

t_end = time.time()
print('Total time = ', t_end - t_start)
plt.show()b_comp.set_points(s_plot.gamma().reshape((-1, 3)))
b_comp._toVTK(out_dir / "magnet_fields", pm_opt.dx, pm_opt.dy, pm_opt.dz)

# Print optimized metrics
print("Total fB = ",
0.5 * np.sum((pm_opt.A_obj @ pm_opt.m - pm_opt.b_obj) ** 2))

bs.set_points(s_plot.gamma().reshape((-1, 3)))
Bnormal = np.sum(bs.B().reshape((qphi, ntheta, 3)) * s_plot.unitnormal(), axis=2)
make_Bnormal_plots(bs, s_plot, out_dir, "biot_savart_optimized")
Bnormal_magnets = np.sum(b_comp.B().reshape((qphi, ntheta, 3)) * s_plot.unitnormal(), axis=-1)
Bnormal_total = Bnormal + Bnormal_magnets

# Compute metrics with permanent magnet results
magnets_m = pm_opt.m.reshape(pm_opt.ndipoles, 3)
num_nonzero = np.count_nonzero(np.sum(magnets_m ** 2, axis=-1)) / pm_opt.ndipoles * 100
print("Number of possible magnets = ", pm_opt.ndipoles)
print("% of magnets that are nonzero = ", num_nonzero)

# For plotting Bn on the full torus surface at the end with just the magnet fields
make_Bnormal_plots(b_comp, s_plot, out_dir, "only_m_optimized")
pointData = {"B_N": Bnormal_total[:, :, None]}
s_plot.to_vtk(out_dir / "m_optimized", extra_data=pointData)

# Print optimized f_B and other metrics
b_comp.set_points(s.gamma().reshape((-1, 3)))
bs.set_points(s.gamma().reshape((-1, 3)))
Bcnormal = np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2)
f_B_sf = SquaredFlux(s, b_comp, -Bcnormal).J()

print('f_B_comp = ', f_B_sf)
total_volume = np.sum(np.sqrt(np.sum(pm_opt.m.reshape(pm_opt.ndipoles, 3) ** 2, axis=-1))) * s.nfp * 2 * mu0 / B_max
print('Total volume = ', total_volume)

t_end = time.time()
print('Total time = ', t_end - t_start)
Expand Down

0 comments on commit 1ba3079

Please sign in to comment.