Misaligned target and eq fields in omnigenity objective

[rahulgaur104]

trafficstars

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.

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]

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

[dpanici]

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]

Let me know if #1166 resolves this

[rahulgaur104]

Unlikely, but I can try.

[rahulgaur104]

:rofl:

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

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

[dpanici]

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]

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

[rahulgaur104]

The alignment can be fixed using the flip_theta function from desc.compat.

[dpanici]

Does that fix this issue? what if the theta=0 point was not on the midplane?

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

[rahulgaur104]

Yes, flip_theta function works in a way that after using it you get min(B) at theta = zeta = 0 instead of max(B) at theta=zeta=0. For the two or three OH cases I looked at, it works. It is possible this will fail for non-stellarator symmetric cases but I am not working on that right now.