Implement FourierRZCurveWindingSurface Class

[dpanici]

Resolves #838

TODO:

• Add Math
• Add figs,
• re-point PR to to Add ShareParameters Linear Objective #1320 so that can actually test optimization with multiple of these curves
• think of what happens when ratio of secular terms is constant but each individual secular term is increased (it increases the rate at which the curve is traversed I think, so increases dx/ds for instance, but the topology remains the same)
• connection to umbilic if secular term in s for zeta is larger than that for theta and is a rational number
• Fix compute quantiites that call surface.compute("R_tt") to use e_theta_t instead (etc...) requires Fix incorrect toroidal components of surface basis vectors #924 and then use the newly added FourierRZToroidalSurface parametrizations for R_tt etc, so that things like torsion of these curves can be calculated
• check that phi_tt etc are implemented as computable quantities from surfaces (or implement them as indexing into the basis vectors) Fix incorrect toroidal components of surface basis vectors #924
• Decide if surface should be passed to the compute as a Transform or as a param and be an optimizable_param -> i.e. do we want the surface on which the curves are defined to be optimized as well? (and how would it be an optimizable param, as the surface itself is an entire thing, but maybe we could just use R_lmn and Z_lmn of the surface as an optimizable param... would probably need to make new compute funs thuogh for this class if we wanted to say restrict the curvature of the surface or its aspect ratio...)
  • I implemented it by having the curve class have R_lmn and Z_lmn params, which are optimizable and inside the compute functions for this class get passed to the compute call when computing the underlying surface geometry. This lets the surface geometry be changed in order to optimize the coils, and for the info on the deriv of for example length to know about how the surface.R_lmn and Z_lmn affect its calculation. This will not let us, for example use PlasmaVesselDistance on the surface these things lie on, but if all we care about are the resulting coils, then this should work fine.
• Add a from_values method to allow fits (probably only need to pass in theta, zeta values and the surface, possibly allowing the secular terms to be specified instead of fit if the topology is already known)
• Deal with NFP when computing coil quantities
• Make secular terms be ints and NOT optimizable, as I think they must be integers for the curve to be closed
• if surface is allowed to optimize, needs the ShareParameters to keep R_lmn Z_lmn are shared in the CoilSet

MATH

idea is originally from a paper on COILOPT by Strickler in 2002

The idea is to create a curve class which is restricted to lie on a given winding surface. As it is implemented now (this may/should change to allow the underlying surface to be optimized as well, but unsure how best to do so, could give this class the R_lmn, Z_lmn and have them point to the surface classes R_lmn, Z_lmn, but for now is useful as a tool for use with the REGCOIL work)*]

To enforce the restriction to a surface, let's say the surface is known as $R_b(\theta,\zeta), Z_b(\theta,\zeta)$ in cylindrical coordinates (usually as a Fourier series but the exact implementation is irrelevant). Then, to define curves on this surface, we will consider a curve on the surface $R_c(s) = R_b((\theta(s),\zeta(s)), Z_c(s) = Z_b(\theta(s),\zeta(s))$, where $s$ is a curve-following parameter $s\in [0,2\pi]$ that parametrizes the curve. By defining the curve position as simply being the surface positions evaluated at some set of theta and zeta coordinates, we have ensured that the curve lies on the surface (i.e any point $R_c, Z_c$ by definition lies on the surface $R_b, Z_b$). The parametrization of the curve then falls to simply parametrizing the angle functions $\theta(s), \zeta(s)$ as functions of $s$. The parametrization is:

$$ \begin{align} \theta(s) &= \theta_{sec} s + \sum_n \theta_n \mathcal{F}^n(s)\\ \zeta(s) &= \zeta_{sec} s + \sum_n \zeta_n \mathcal{F}^n(s)\\ \end{align} $$

Where $\theta_{sec}, \zeta_{sec}$ are the secular terms in $s$ and $\theta_n,\zeta_n$ are Fourier coefficients associated with the periodic part. The secular terms are necessary to allow a curve to wrap around the torus either poloidally (if the $\theta_{sec}$ term is nonzero) or toroidally (if the $\zeta_{sec}$ term is nonzero) or both (if both secular terms are nonzero).

The simplest example is if we want a curve that is at a constant $\zeta$ but completes one poloidal turn around the surface, we simply set every coefficient to zero except the $\theta_sec$ term which we set to $1$, yielding a curve which completes one poloidal turn as $s$ goes from $0$ to $2\pi$

The toroidal analog is the same thing except $\zeta_{sec}=1$ instead, and this yields a curve of constant $\theta$ going around the surface toroidally.

Helicity of a curve is determined by the ratio of $\theta_{sec}$ to $\zeta_{sec}$, a helical curve that goes around once poloidally each toroidal turn is given by setting both of these terms to 1, for instance.

Curve Closure Conditions

For the curve to be closed, it must come back to itself after $2\pi$ in the curve parameter i.e.

$(\theta(s=0),\zeta(s=0)) = (mod(\theta(s=2/pi), mod(\zeta(s=2/pi))$

Looking at the above representations of the angles in $s$, we see after subtracting the LHS from the RHS (cancelling out the periodic fourier series) leaves us with just the secular terms

$(0,0) = (mod(2\pi\theta_{sec}, 2\pi) - \theta_{sec}, mod(2\pi\zeta_{sec}, 2\pi) - \zeta_{sec})$

therefore, the only way this is true is if the secular terms are integers.

[github-actions]

|             benchmark_name             |         dt(%)          |         dt(s)          |        t_new(s)        |        t_old(s)        | 
| -------------------------------------- | ---------------------- | ---------------------- | ---------------------- | ---------------------- |
 test_build_transform_fft_lowres         |     -0.72 +/- 3.79     | -3.88e-03 +/- 2.04e-02 |  5.33e-01 +/- 1.6e-02  |  5.37e-01 +/- 1.3e-02  |
 test_equilibrium_init_medres            |     -0.24 +/- 1.14     | -1.01e-02 +/- 4.81e-02 |  4.21e+00 +/- 4.2e-02  |  4.22e+00 +/- 2.4e-02  |
 test_equilibrium_init_highres           |     -0.10 +/- 1.00     | -5.41e-03 +/- 5.56e-02 |  5.56e+00 +/- 3.2e-02  |  5.57e+00 +/- 4.6e-02  |
 test_objective_compile_dshape_current   |     -0.21 +/- 0.75     | -8.09e-03 +/- 2.91e-02 |  3.87e+00 +/- 1.2e-02  |  3.87e+00 +/- 2.7e-02  |
 test_objective_compute_dshape_current   |     +0.61 +/- 3.47     | +2.11e-05 +/- 1.20e-04 |  3.49e-03 +/- 1.1e-04  |  3.47e-03 +/- 5.0e-05  |
 test_objective_jac_dshape_current       |     +0.22 +/- 4.32     | +9.05e-05 +/- 1.79e-03 |  4.15e-02 +/- 1.3e-03  |  4.14e-02 +/- 1.2e-03  |
 test_perturb_2                          |     -0.55 +/- 0.97     | -9.80e-02 +/- 1.72e-01 |  1.76e+01 +/- 9.1e-02  |  1.77e+01 +/- 1.5e-01  |
 test_proximal_freeb_jac                 |     -0.48 +/- 1.89     | -3.58e-02 +/- 1.40e-01 |  7.40e+00 +/- 1.2e-01  |  7.43e+00 +/- 6.8e-02  |
 test_solve_fixed_iter                   |     +1.45 +/- 59.01    | +7.20e-02 +/- 2.94e+00 |  5.05e+00 +/- 2.1e+00  |  4.98e+00 +/- 2.1e+00  |