Incorrect computations with stellarator symmetry

[unalmis]

We have an inconsistent use of stellarator symmetry (the discrete one) that leads to computing quantities incorrectly.
Here is a toroid surface with an elliptic cross section and torsion.

surf = FourierRZToroidalSurface(
    R_lmn=[10, 1, 0.2],
    Z_lmn=[-2, -0.2],
    modes_R=[[0, 0], [1, 0], [0, 1]],
    modes_Z=[[-1, 0], [0, -1]],
)

DESC determines this is stellarator symmetric, and will also, by default, label equilibrium with such a boundary as a symmetric equilibrium. We use the same characteristic to construct computational grids for stellarator symmetry that truncate the poloidal domain from $\theta \in [0, 2 \pi]$ to $\theta \in [0, \pi]$. By doing so, we are assuming that, under stellarator symmetry, the following relation holds $\int_0^{2\pi} f(\theta) d \theta = 2 \int_0^{\pi} f(\theta) d \theta$, which is not true because f is not symmetric.

For example, this quantity is not computed correctly on symmetric objects such as the surface above when integrated along a boundary, so the Elongation and Aspect ratio objectives that target toroid surfaces introduced in #884 have never been correct on master for symmetric surfaces (the equilibrium ones are fine edit: fine for this surface).

One may expect that the same issue exists in the volume computation, and hence the magnetic well, but those also integrate over zeta, and even though their integrands, assuming symmetry, do not satisfy $\int_0^{2\pi} f(\theta) d \theta = 2 \int_0^{\pi} f(\theta) d \theta$, I think they do satisfy $\int_0^{2\pi} \int_0^{2\pi} f(\theta) d \theta d\zeta = 2 \int_0^{2\pi} \int_0^{\pi} f(\theta) d \theta d\zeta$ (one can show analytically they do for above surface), which is why all the other tests for correctness of this work. Still we should review code and review the symmetry condition.

(fyi these results were confirmed/tested with analytic tests using sympy).

[unalmis]

I read R.L. Dewar, S.R. Hudson, Stellarator symmetry, doi 10.1016/S0167-2789(97)00216-9.

Recall stellarator symmetry is a property of a curvilinear coordinate system, $(\rho, \theta, \zeta)$, such that $f(\rho, \theta, \zeta) = f(\rho, -\theta, -\zeta)$ Dewar, Hudson eq.8. The DESC coordinate system will be a stellarator symmetric coordinate system if the Fourier expansion of the flux surfaces have either the cosine or sine symmetry.

Now, assuming stellarator symmetry gives the first relation $$f(\rho, -\theta, -\zeta) = f(\rho, \theta, \zeta) \neq f(\rho, -\theta, \zeta \neq 0)$$ but the second relation does not follow (hence the $\neq$). So we should not expect any of our computations to be invariant to truncating the poloidal domain to above the midplane $\theta \in [0, \pi] \subset [0, 2 \pi)$.

If we are computing some function $g \colon \rho, \theta, \zeta \mapsto g(\rho, \theta, \zeta)$ that is just a pointwise evaluation of the basis functions, then we will of course still compute $g$ accurately above the midplane. However, if we are computing any function that is not a pointwise evaluation of the basis function, i.e. a function whose input takes multiple nodes as input and performs some type of reduction, e.g. $F \colon \rho, \theta, \zeta \mapsto \int f(\rho, \theta, \zeta) d S$, then in general $F$ will not be computed accurately if the computational domain is truncated to above the midplane.

In general, given

1. $f$ that evaluates the basis functions pointwise
2. $F$ that performs a reduction on $f$
3. stellarator symmetry: $f(\rho, \theta, \zeta) = f(\rho, -\theta, -\zeta)$

then $F$ is guaranteed to be able to be computed accurately on the truncated domain of computation $\theta \in [0, \pi] \subset [0, 2\pi)$ only¹ if $F$ is a linear reduction over $D \equiv [0, \pi] \times [0, 2 \pi) \ni (\theta, \zeta)$.

This means that if $F$ is a flux surface integral or volume integral of $f$, then it can be computed on grids that have nodes only above the midplane, i.e. grids such that grid.sym == True.

If $F$ is a nonlinear reduction or any reduction that is over a proper subset of $D$, then $F$ may not be computed accurately when the domain is truncated to above the midplane unless there is the additional symmetry
$$f(\rho, \theta, \zeta) = f(\rho, -\theta, \zeta)$$
Stellarator symmetry implies this relation holds for $\zeta = 0$. Therefore, stellarator symmetry and $\partial f / \partial \zeta = 0$ is sufficient, but not necessary, for this additional symmetry.

This means that if $F$ is a non-flux surface integral or line integral, then it cannot be computed accurately on grids that have nodes only above the midplane, i.e. grids such that grid.sym == True, unless the additional symmetry is satisfied.

For example, on the surface given at the top of this issue, A(z) cannot be computed on the boundary from a line integral over $\theta$ above the midplane. However, it can be computed from a surface integral of the cross section above the midplane. This is supported from the analytical tests done here. For a more general surface, e.g. one with $\partial_{\zeta} \Vert e_{\rho} \times e_{\theta}\Vert \neq 0$, the latter approach will also fail.

TODO:

• Add the explanations in this issue to Dev guide #353 .
• Use more complicated surfaces in tests.
• Add warnings to quantities that assume this additional symmetry to compute on symmetric grids in Integrate on boundary to compute length scale quantities #1094 .
• Update documentation of Grid class to reflect that grid.sym is unrelated to stellarator symmetry, rather it's a memory optimization to compute certain quantities under symmetry.
• Maybe someone can also check if this affects stuff like Galerkin and Spectral Condensation #382? If one was running a stellarator a symmetric stellarator / constraining the basis to be symmetric there, try the same but turn off the symmetry of all the computational grids, while preserving the symmetry of the basis.²

Footnotes

1. Since we have examples that show it won't work for subsets of $D$, I think I'm justified in claiming "only if reduction over $D$". I didn't bother to generate an example for the nonlinear claim, so if you doubt the linear qualifier in "only if linear reduction over $D$", you are justified, since I didn't show that. ↩

2. I'm thinking of the formulation of weighted residual pseudospectral methods in chapter 3 of Boyd's spectral book. If the weight function $w$ assigned to each collocation point does not possess stellarator symmetry, or if $w$ is a reduction like $F$ above that is not (linear and over $D$) then one should use the full grid to compute. In the Galerkin approach $w$ assigned to a collocation point is just the stellarator symmetric basis function evaluated there, so $w$ should be stellarator symmetric function. Then I think as long as the optimizer is doing some volume or flux average of w * f, it should be fine to use grids with nodes only above the midplane. ↩

[unalmis]

I won't get around to doing the other checkboxes, at least for foreseeable future.

[rahulgaur104]

I agree with you but the general definition of stellarator symmetry should also include an anti-symmetry relations. For example, a stellarator symmetric equilibrium satisfies R(theta, zeta)= R(-theta, -zeta) AND Z(theta, zeta) = -Z(-theta, -zeta).

In general a geometric quantity is function of R, Z and their derivatives.

[dpanici]

Mark as requires tz resolution and then do calculations (adding array flipped perhaps and dividing by 2) to then stitch together the full correct result from the sym grid calculation.