Compute A(r)

[ddudt]

Presently we have a compute function for cross-sectional area A, which is a scalar quantity for the area through the $\rho=1$ surface. I would like to compute the same quantity but as a flux function, $A(\rho)$.

[f0uriest]

This should be possible using Stokes theorem: $\int \intD (\nabla \times \mathbf{A}) \cdot d\mathbf{S} = \int{\partial D} \mathbf{A} \cdot d\mathbf{l}$

$d\mathbf{S} = \mathbf{n} |e\rho \times e\theta| d\rho d\theta$

So we want $\nabla \times \mathbf{A} = \mathbf{n}$

Any ideas?

[unalmis]

Let's pretend we renamed the cross sectional area as $B(r)$ to not get confused with existing notation on wikipedia. We want to generate a vector potential whose curl points along $e^{\zeta}$. Or at least one where $(\nabla \times A) \cdot e^{\zeta}$ is a constant. $e^{\zeta}$ is divergence free, so $A$ exists. The task becomes to equate these components of $e^{\zeta}$: with the components of $\nabla \times A$ given here: https://en.wikipedia.org/wiki/Del_in_cylindrical_and_spherical_coordinates#Del_formula.

I haven't done the math for this, but it is a standard calculus problem. It's possible there would need to be some integrals of $R(\rho, \theta, \zeta)$ etc. though...

Are we sure we want to do this integral along its boundary? Using this technique to compute $V(\rho)$ was a good idea because

1. The integration becomes more tractable. We only need to consider n^2 rather than n^3 points, so this also leads to higher accuracy
2. Transforming the integration to one over periodic domains makes it more accurate, especially at low resolution, (midpoint rule integral on periodic domains converge exponentially to true integral).
3. We knew a cheap to compute transformation with $\nabla \cdot Z = 1$.

For $B(\rho)$, point 2 is still valid, but maybe less important as I'd expect an integral over less dimensions to have less issues with accuracy. For point 1 -- if we still want to take the mean of the cross sectional area over all $\zeta$ surfaces, then while integrating along the boundary reduces the dimension of each integral it doesn't reduce the number of integrals that need to be done. We would still need to do grid.num_zeta * grid.num_rho integrals before taking the mean across the cross sections.

If the vector potential ends up being complicated, I think we could fall back to computing $\partial_{\rho} B(\rho, \zeta)=$ line_integrals(grid, q=q, line_label="theta", fix_surface=(rho, [grid.unique_rho_idx])). Then integrating over $\rho$ with np.cumtrapz like we did with chi, then average over $\zeta$ for the mean.

[unalmis]

Equation 2.6.37 in D'haeseleer's book actually gives the vector potential for the covariant basis vector. We want the contravariant one though, so maybe it can come out of that.

[unalmis]

There's also this https://en.wikipedia.org/wiki/Helmholtz_decomposition

[unalmis]

Instead of finding a vector potential it's easier to integrate along the boundary like:

$\iintD dS (\nabla \cdot A) = \int{\partial D} d\ell (n \cdot A)$ where $n$ is the normal vector to the coordinate curve. Since we want $d \ell$ to parametrize the boundary of a $\rho$ surface at constant $\zeta$, the tangent vector to the curve points in direction $e{\theta}$, so the normal vector points in direction $n = \partial{\theta} e_{\theta}$.

def other_way_to_compute_area(eq, grid):
    data = eq.compute(["g_tt", "Z", "e_theta_t", "rho"], grid=grid)
    # line length Jacobian, dl or tangent vector magnitude
    dl = np.sqrt(data["g_tt"])
    # unit normal vector to coordinate curve
    n = -(data["e_theta_t"].T / np.linalg.norm(data["e_theta_t"], axis=-1)).T
    # 2d cylindrical divergence (fix zeta) of Z over zeta surface is 1
    q = dl * n[:, 2] * data["Z"]  # Z is the A in divergence theorem in github comment above

    # In theory this part can be made into just 1 integral with some changes to line_integrals function
    A = np.empty(grid.num_rho)  # this is the cross-sectional area A(r) as function of rho
    for i, r in enumerate(grid.compress(data["rho"])):
        A[i] = np.mean(
            line_integrals(
                grid=grid,
                q=q,
                line_label="theta",
                fix_surface=("rho", r),
                expand_out=False,
            )
        )
    # axis limit is 0
    return A

[unalmis]

I've tested the above code matches what we expect analytically for a torus (on eq = Equilibrium()). When I compare B_of_r[-1] to data["A"] on eq = examples.get("W7-X") I get different but reasonable answers (1 vs 0.8, respectively). If we agree that my math above is correct, then B_of_r[-1] is in theory more accurate. @ddudt can you check sometime that the above code gives a reasonable result for your use case?

[ddudt]

@unalmis this does not work for elliptical cross-sections. I tested with both axisymmetric ellipses and with rotating ellipses. In all cases, your code underestimates the correct value (but converges to the correct answer for circular cross-sections).

The "correct" implementation I was testing with was adapted from some ASCOT source code to compute the cross-sectional area of VMEC equilibria in the phi=0 plane:

def area(eq, grid):
    data = eq.compute(["R", "Z"], grid=grid)
    # reshape to (num_rho, num_theta, num_zeta)
    r = np.transpose(
        data["R"].reshape((grid.num_theta, grid.num_rho, grid.num_zeta), order="F"),
        (1, 0, 2),
    )
    z = np.transpose(
        data["Z"].reshape((grid.num_theta, grid.num_rho, grid.num_zeta), order="F"),
        (1, 0, 2),
    )
    A = np.zeros(r.shape[0])
    for i in range(A.size):
        # only computing for the zeta=0 cross-section
        x = r[i, :, 0]
        y = z[i, :, 0]
        # can this be replaced with a line integral?
        A[i] = 0.5 * np.abs(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1)))
    return A

This function agrees with our existing implementation eq.compute("A")["A"] for the boundary surface. I think it is a crude trapazoidal integration and we might be able to do something similar.

[unalmis]

I don't know why it wouldn't work. The curl line integral probably won't work either then. Here's the basic trapezoidal integration.

def area_trapezoidal_looped(eq, grid):
    data = eq.compute(["|e_rho x e_theta|", "rho"], grid=grid)
    rho = grid.compress(data["rho"])
    dA_dr = np.empty((grid.num_rho, grid.num_zeta))
    for i, r in enumerate(rho):
        dA_dr[i] = line_integrals(
            grid=grid,
            q=data["|e_rho x e_theta|"],
            line_label="theta",
            fix_surface=("rho", r),
            expand_out=False,
        )
    A = np.mean(cumtrapz(dA_dr, rho, axis=0, initial=0), axis=-1)
    return A

[f0uriest]

What about doing something like V(r)/R0?

[ddudt]

The trapezoidal integration does work, but so does V(r) / (2*pi*R0) (they both give the "correct" results). The later is probably the most elegant solution.

[unalmis]

@ddudt @f0uriest I figured out why the original line integral failed for ellipses. The original code I provided assumes the curve is parameterized with unit speed (e.g. in the Fenet Frame or parameterized by arc length). For eliptical cross sections, a curve parameterized by $\theta$ does not travel at unit speed.

For a parameterization with non-constant speed, the unit normal vector cannot be computed as the derivative of the tangent vector of the curve (i.e. the acceleration) because some component of the acceleration vector must point along the tangent curve to change the speed of parameterization. In general for curves with arbitrary parameterizations, we can compute the unit normal vector as the acceleration vector minus the component of the acceleration vector along the tangent vector.

The integration along the boundary is provided below is correct and avoids the numerical errors that come with a trapezoidal integration. If V(r) / (2*pi*R0) always works then yea that is probably nicer.

def area_boundary_line_integral(eq, grid):
    data = eq.compute(["g_tt", "Z", "e_theta", "e_theta_t", "rho"], grid=grid)
    # line length Jacobian (tangent vector magnitude)
    dl = np.sqrt(data["g_tt"])
    # unit tangent vector
    tangent = (data["e_theta"].T / dl).T
    # outward unit normal vector
    normal = -(data["e_theta_t"].T - dot(data["e_theta_t"], tangent) * tangent.T).T
    normal = (normal.T / np.linalg.norm(normal, axis=-1)).T
    # 2d cylindrical divergence (fix zeta) of Z over zeta surface is 1
    integrand = dl * normal[:, 2] * data["Z"]

    # with some changes to line_integrals function this for loop can be removed
    A = np.array(
        [
            line_integrals(
                grid,
                integrand,
                line_label="theta",
                fix_surface=("rho", r),
                expand_out=False,
            ).mean()
            for r in grid.compress(data["rho"])
        ]
    )
    # axis limit is 0
    return A