Apparent mismatch between ballooning objective output and direct calculation of growth rates

[f0uriest]

setup:

eq0 = desc.examples.get("HELIOTRON")

nalpha = 8  # Number of field lines

# Field lines on which to evaluate ballooning stability
alpha = np.linspace(0, np.pi, nalpha, endpoint=False)

# Number of toroidal transits of the field line
nturns = 3

# Number of point along a field line per transit
nzetaperturn = 200

rho=0.8

obj2 = desc.objectives.BallooningStability(
            eq=eq0,
            rho=rho,
            alpha=alpha,
            nturns=nturns,
            nzetaperturn=nzetaperturn,
            weight=1,
            target=0,
        )
obj2.build()

We expect the following 3 to give the same answer since the objective has target=0 and weight=1

print(obj2.compute(eq0.params_dict, constants=obj2.constants))
print(obj2.compute_scaled_error(eq0.params_dict, constants=obj2.constants))
print(obj2.compute_scalar(eq0.params_dict, constants=obj2.constants))

which gives

221.8172472905935
[221.81724729]
221.81724729059223

which seems fine

But if we then compute the growth rates directly and apply the same reduction as in the objective function:

# range of the ballooning coordinate zeta
zeta = np.linspace(-np.pi * nturns, np.pi * nturns, nzetaperturn*nturns)

grid = desc.grid.Grid.create_meshgrid(
    [rho, alpha, zeta], coordinates="raz", period=(np.inf, 2 * np.pi, np.inf)
)

# make sure its on the same points
np.testing.assert_allclose(alpha, obj2.constants['alpha'], atol=1e-14)
np.testing.assert_allclose(zeta, obj2.constants['zeta'], atol=1e-14)

ball_data0 = eq0.compute(["ideal ballooning lambda"], grid=grid)
lam = ball_data0["ideal ballooning lambda"]

# same as in obj2.compute
lambda0, w0, w1 = obj2.constants["lambda0"], obj2.constants["w0"], obj2.constants["w1"]

# Shifted ReLU operation
data = (lam - lambda0) * (lam >= lambda0)

results = w0 * np.sum(data) + w1 * np.max(data)
print(results)

gives

14.325405153826495

which is way smaller. Am I missing something?

[rahulgaur104]

Different normalizations due to different values of data["a"]
For the first case it's 1.506 and
for the second case it's 0.95
which causes all the coefficients to change by a constant factor.
I think the second one might be wrong and we need to explicitly provide data["a"] to the compute function.

The easiest fix would be do exactly what the objective does: calculate "a" on a boundary grid and do
eq0.compute("ideal ballooning lambda", grid=grid, data = {"a":<a_from_boundary_grid>})
I have tested it by setting <a_from_boundary_grid> = 1.506 and it matches with the compute value.

I'll create a hotfix and add a test today.

[rahulgaur104]

Actually the data["a"] in the objective is incorrect and I can't figure out the reason.
When I do
eq.compute("a", grid=LinearGrid(rho=1.0, M=eq.M_grid, N=eq.N_grid, NFP=eq.NFP))[["a"]
I get the right answer(0.95) but when I print
compute_fun(eq, ["a"], params=params, transforms=constants["len_transforms"], profiles=constants["len_profiles"])["a"]
which is exactly what is happening inside the BallooningStability objective, I get a wrong value(1.506). I don't know why this is happening.

Have you seen this before? Could "constants" be causing issues?
I don't know why computing the same quantity inside objective gives a different answer.

[f0uriest]

I thought stuff like a should be computed on a quadrature grid, not a single surface grid? I think in eq.compute we override the user grid to do that, so that's likely why its giving a different answer. If that's the case the fix would just be to use a quadrature grid for len_transforms etc, but I thought there was some reason we didn't do that originally?

[rahulgaur104]

What I mean by "a" the effective minor radius which only depends on the boundary shape. Why would it need a quadrature(volume) grid? We just calculate the effective A(z) on the boundary and from that, calculate "a".
Also, a quadrature grid will never have rho points exactly on the boundary, right?

[f0uriest]

A(z) for an equilibrium is computed as a surface integral over the cross section, not a boundary integral:

DESC/desc/compute/_geometry.py

Lines 154 to 180 in 4fd41f1

	@register_compute_fun(
	name="A(z)",
	label="A(\\zeta)",
	units="m^{2}",
	units_long="square meters",
	description="Cross-sectional area as function of zeta",
	dim=1,
	params=[],
	transforms={"grid": []},
	profiles=[],
	coordinates="z",
	data=["|e_rho x e_theta|"],
	parameterization=[
	"desc.equilibrium.equilibrium.Equilibrium",
	"desc.geometry.surface.ZernikeRZToroidalSection",
	],
	resolution_requirement="rt",
	)
	def _A_of_z(params, transforms, profiles, data, **kwargs):
	data["A(z)"] = surface_integrals(
	transforms["grid"],
	data["|e_rho x e_theta|"],
	surface_label="zeta",
	expand_out=True,
	)
	return data

I think to only use a boundary grid needs #1094

[rahulgaur104]

Ok, then let me use quadrature grid for "a" inside the BallooningObjective. For calculating the eigenvalue using eq.compute, I'll pass "a" as data computed on a quadrature grid.

[unalmis]

The documentation for adding objectives would benefit from a section that tells the developer to make sure their dependencies are computed on the right grids, as it is not done automatically as is done through eq.compute.