Error when computing Poloidal flux and Boozer toroidal stream function

[IronLSY]

There exists errors when computing Boozer toroidal stream function and Poloidal flux in the example shown below, which uses the DSHAPE_output from DESC packge.

import numpy as np
import desc.io
from desc.grid import Grid
from desc.equilibrium import Equilibrium

#load the DSHAPE equilibrium
eq_fam = desc.io.load("DSHAPE_output.h5")

#node and node_list both contain the point (rho=0.5, theta=1.5PI, zeta=0)
node = Grid(np.array([0.5,1.5*np.pi,0]))
node_list = Grid(np.array([[0.5,1.5*np.pi,0], [0.5,1.7*np.pi,0]]))

#print their nodes to check if the points contained are correct
print('node:')
print(node.nodes)
print('node_list:')
print(node_list.nodes)

#computing the Boozer toroidal stream function and output
print('nu:')
print(eq_fam[0].compute('nu', node)['nu'])
print(eq_fam[0].compute('nu', node_list)['nu'])

#computing the Poloidal flux and output
print('chi:')
print(eq_fam[0].compute('chi', node)['chi'])

The output shown below hints some errors when computing nu and chi.
The nu should be the same at the same point (rho=0.5, theta=1.5PI, zeta=0), but it shows difference.
The chi can not be zero at the point (rho=0.5, theta=1.5PI, zeta=0) since rho ≠ 0.

node:
[[0.5        4.71238898 0.        ]]

node_list:
[[0.5        4.71238898 0.        ]
 [0.5        5.34070751 0.        ]]

nu:
[0.00400288]
[0.00665903 0.00363253]

chi:
[0.]

[f0uriest]

Hi IronLSY, thanks for the report!

I think both of these stem from the same cause, namely computing quantities at a point when the quantities require information at other points.

For example, we compute chi by integrating chi_r (radial derivative of the poloidal flux):

DESC/desc/compute/_profiles.py

Lines 143 to 146 in 8d2b96f

	chi_r = transforms["grid"].compress(data["chi_r"])
	chi = cumtrapz(chi_r, transforms["grid"].compress(data["rho"]), initial=0)
	data["chi"] = transforms["grid"].expand(chi)
	return data

So you would need to compute it from the core outwards to get the correct value at $\rho=0.5$. For most quantities we have logic in Equilibrium.compute to use the correct grids etc but I think not for this case, I'd have to think about how to make sure we always compute chi correctly.

The Boozer stream function is related I think, in that in order to compute it you need values over the entire flux surface. nu depends on the spectral coefficients w_Booozer_mn, which in turn depend on the spectral coefficients of B on the surface, so you need to include a full grid of points on that surface to properly compute the B_theta/zeta_mn harmonics (the Boozer stuff is also tricky because currently it only works on 1 surface at a time, though we're working on generalizing that in #369)

[IronLSY]

Thank you so much! This really helped me a lot. So now if I want to get the value of nu as accurate as possible on a specific magnetic surface, should I let theta take as many values as possible from (0,2PI)?

[f0uriest]

and zeta as well if its non-axisymmetric. You probably want to compute nu over a whole surface with a grid like LinearGrid and then interpolate:

grid = LinearGrid(rho=0.5, M=100, N=100, NFP=eq.NFP)
data = eq.compute("nu", grid=grid)
theta = grid.nodes[grid.unique_theta_idx, 1]
zeta = grid.nodes[grid.unique_zeta_idx, 2]
nu = data['nu'].reshape((grid.num_theta, grid.num_zeta))

theta_q = ... # where you want to interpolate
zeta_q = ... # where you want to interpolate

from desc.interpolate import interp2d
nu_q = interp2d(theta_q, zeta_q, theta, zeta, nu, period=(2*np.pi, 2*np.pi/eq.NFP)

[ddudt]

Rory beat me to it, but here are other examples of how you can accurately compute both values at the coordinates you wanted.

Poloidal flux:

grid = LinearGrid(L=100, theta=1.5*np.pi, zeta=0.0, NFP=eq.NFP)
data = eq.compute("chi", grid)
chi = data["chi"][grid.L // 2]  # chi at rho=0.5

Boozer toroidal stream function:

from desc.interpolate import interp1d

grid = LinearGrid(rho=0.5, M=2 * eq.M_grid, zeta=0.0, NFP=eq.NFP)  # ensure grid contains a lot of points
data = eq.compute(["theta", "nu"], grid, M_booz=2 * eq.M, N_booz=2 * eq.N)  # ensure Boozer spectrum is reasonably large
nu = interp1d(np.array([1.7 * np.pi]), data["theta"], data["nu"], period=2 * np.pi)  # nu at rho=0.5, theta=1.7PI, zeta=0

[IronLSY]

Thank you guys for your patient answers, the interpolation works well !

By the way, what's the definition of Jacobian determinant of Boozer coordinates: $\sqrt{g}_B$ ?
I think the Boozer poloidal angular coordinate $\theta_B$ is defined as $\theta_B=\theta+\lambda+\iota\nu$ and I check in the output that this relationship is correct.

However I thought: $\sqrt{g}_B=\sqrt{g}/\left[\partial\psi/\partial\rho\cdot\left(1+\partial\lambda/\partial\theta+\iota\partial\nu/\partial\theta\right)\right]$, where the variables above are consistent with your definition in the List of Variables. Or in the code style:sqrt(g)_B = sqrt(g)/[psi_r*(1+lambda_t+iota*nu_t)]. I find the output mismatch this relationship, then I guess replacing psi_r with chi_r may work but the output still mismatch(chi_r and nu_t is already integrated correctly).

So I'm confused about the definition of Jacobian determinant of Boozer coordinates: $\sqrt{g}_B$ .

Could you give me an analytical expression for it?

[f0uriest]

This is likely poor notation on our part.

Our "regular" sqrt(g) is the determinate of the coordinate transform between DESC coordinates $(\rho, \theta, \zeta)$ and the lab frame $(R, \phi, Z)$, and same for sqrt(g)_PEST = $(\rho, \theta_{PEST}, \zeta) -> (R, \phi, Z)$

However, what we call sqrt(g)_B is really the determinate of the coordinate change between boozer coordinates and DESC coordinates (not the lab frame): sqrt(g)_B = $(\rho, \theta, \zeta) -> (\rho, \theta_B, \zeta_B)$

The actual equation is $$\sqrt{g}_B = (1 + \partial \lambda /\partial \theta) * (1 + \partial \nu / \partial \zeta ) + (\iota - \partial \lambda / \partial \zeta) * \partial \nu / \partial \theta$$

Which in axisymmetry reduces down to what you have above but without the factor of $\sqrt{g}$

[IronLSY]

Thanks very much, and wish you all the best in your code development!

[dpanici]

@ddudt @f0uriest do we want to update this notation?

In the booz_xform style output PR #680 , when saving that file one of the fields is gmn_b which is the Boozer jacobian determinant $\sqrt{g}_B = \frac{\iota I + G}{B^2}$

to make this in our code I think we have to divide sqrt(g) by sqrt(g)_B, but then what do we call this new quantity? sqrt(g)_B is already taken... I vote we rename it to something jacobian_Boozer_DESC or something, so we can have sqrt(g)_B as its proper name

[f0uriest]

I'd be fine with that