PlasmaControl / DESC

Stellarator Equilibrium and Optimization Suite
MIT License
87 stars 22 forks source link

Misaligned target and eq fields in omnigenity objective #1143

Open rahulgaur104 opened 1 month ago

rahulgaur104 commented 1 month ago

For the attached data.zip OH equilibrium, the |B| in Boozer coordinates using the equilibrium object, and a field object have max and min(|B|) at different points.

eq_modB field_modB

Since the definition of Boozer coordinates is the same, this implies that the field object's |B| is unphysical and incorrect.

To reproduce the result, run the script below with the data inside the attached zip file.

import pdb
import os
import numpy as np
import matplotlib.pyplot as plt
from desc.equilibrium import Equilibrium
from desc.grid import LinearGrid
from desc.plotting import *

fname_path1 = "eq_final.h5"
fname_path2 = "field_final.h5"

eq0 = Equilibrium.load(f"{fname_path1}")
field = Equilibrium.load(f"{fname_path2}")

N = int(200)
grid = LinearGrid(L = N)
rho = np.linspace(0, 1, N+1)

data_keys = ["iota", "D_Mercier"]
data1 = eq0.compute(data_keys, grid = grid)

rho0 = 1.0
iota = data1["iota"]

fig, ax, Boozer_data0 = plot_boozer_surface(eq0, rho=rho0, return_data=True)
plt.show()
plt.close()

fig, ax, Boozer_data1 = plot_boozer_surface(field, rho=rho0, iota=np.interp(rho0, rho, iota), return_data=True)
plt.show()
plt.close()
f0uriest commented 1 month ago

Also looks like |B| in the 2nd plot goes negative? and are the big white areas where its NaN or something?

dpanici commented 1 month ago

I think the point here is, even for a well behaved omni field, the max |B| point in $(\theta_B, \zeta_B)$ space might be misaligned with that of the equilibrium when its |B| is converted to Boozer coordinates. This could cause the optimizer to have to work harder to first align the |B| contours by re-parametrizing the surface poloidal angle, and then it can do the actual work of making the field omnigenous. At least is my interpretation of this

unalmis commented 1 month ago

Let me know if #1166 resolves this

rahulgaur104 commented 1 month ago

Unlikely, but I can try.

rahulgaur104 commented 1 month ago

:rofl:

After using #1166 it actually looks worse. Run the test script with the data shared in this issue and you shall see this.

target

If you don't want to use the field from the zip file shared at the top of issue, you can create your own using the following script.

import pdb
import os
import numpy as np
import matplotlib.pyplot as plt

from desc.equilibrium import Equilibrium
from desc.grid import LinearGrid

from desc.magnetic_fields import OmnigenousField

from desc.plotting import *

keyword_arr = ["OH"]

fname_path1 = "eq_final.h5"
fname_path2 = "field_final.h5"

eq0 = Equilibrium.load(f"{fname_path1}")

L_shift = int(2)
M_shift = int(1)
N_shift = int(1)
L_well = int(2)
M_well = int(3)
mirror_ratio = 0.4
well_width = 0.03
radial_shift = 0.002

field = OmnigenousField(
        L_B=L_well,  # radial resolution of B_lm parameters
        M_B=M_well,  # number of spline knots on each flux surface
        L_x=L_shift,  # radial resolution of x_lmn parameters
        M_x=M_shift,  # eta resolution of x_lmn parameters
        N_x=N_shift,  # alpha resolution of x_lmn parameters
        NFP=eq0.NFP,  # number of field periods; should always be equal to Equilibrium.NFP
        helicity=(1, eq0.NFP),  # helicity for toroidally closed |B| contours
        B_lm = np.array([[1-mirror_ratio, 1+well_width, 1+mirror_ratio], [radial_shift, radial_shift, radial_shift], [0., 0., 0,]]).flatten()
        )

N = int(200)
grid = LinearGrid(L = N)
rho = np.linspace(0, 1, N+1)

data_keys = ["iota", "D_Mercier"]

data1 = eq0.compute(data_keys, grid = grid)

rho0 = 1.0
iota = data1["iota"]

fig, ax, Boozer_data0 = plot_boozer_surface(eq0, rho=rho0, return_data=True)
plt.show()
plt.close()

fig, ax, Boozer_data1 = plot_boozer_surface(field, rho=rho0, iota=np.interp(rho0, rho, iota), return_data=True)
plt.show()
plt.close()

and so the target field becomes

target2
dpanici commented 3 weeks ago

Could we add a check to the Omnigenity objective to check if the initial eq's B max is at theta_booz=0? @ddudt

dpanici commented 1 week ago

Add ability to shift B values after transformation to match where eq's Bmax is