pyoculus.geo.interpolate_coordinates

class pyoculus.geo.interpolate_coordinates.SurfacesToroidal(nsurfaces=2, mpol=10, ntor=10, Nfp=1, stellar_sym=True)

Toroidal Surface respresented as a Fourier Harmonic Serie.

Define surfaces \(s(\vartheta, \zeta)\) and the angle transformations \(\theta (\vartheta, \zeta)\) written in terms of Fourier harmonics.

..math:

s(\vartheta, \zeta) &= \sum_{m=0}^{m_\text{pol}}\sum_{n=-n_\text{tor}}^{n_\text{tor}}
s_{c,m,n} \cos(m \vartheta - n n_\text{fp} \zeta) + s_{s,m,n} \sin(m \vartheta - n n_\text{fp} \zeta),
\theta(\vartheta, \zeta) = \vartheta + \sum_{m=0}^{m_\text{pol}}\sum_{n=-n_\text{tor}}^{n_\text{tor}}
\theta_{c,m,n} \cos(m \vartheta - n n_\text{fp} \zeta) +        heta_{s,m,n} \sin(m \vartheta - n n_\text{fp} \zeta).

If stellarator symmetry is assumed, then \(s\) will only have cosine components and \(\theta\) will only have sine components.

Each surface is assigned a radial label \(\rho\).

Now a new set of coordinates

$( ho, artheta, zeta) $ will be constructed

given the relationship

$s( ho, artheta, zeta), heta( ho, artheta, zeta) $.

If

$ ho $ is not on one of the known surfaces, an interpolation will be performed on all the Fourier harmonics radially using the method of the user’s choice.

remove_surface(i: int)

! Remove a surface @param i the id of the surface

add_surface(rho: float, scn, tsn, ssn=None, tcn=None)

! Adding a surface into the system with radial label rho

@param rho the new coordinate

$ ho $ for this new surface

@param scn the cosine components of

$s( ho, artheta, zeta) $

@param tsn the sine components of

$ heta( ho, artheta, zeta) $

@param ssn the sine componets of

$s( ho, artheta, zeta) $

@param tcn the cosine components of

$ heta( ho, artheta, zeta) $

replace_surface(idx: int, rho: float = None, scn=None, tsn=None, ssn=None, tcn=None)

! Replacing a surface by the new one

@param idx the index of the surface to be replaced @param rho the new coordinate

$ ho $ for this new surface. If this is None then keep same.

@param scn the cosine components of

$s( ho, artheta, zeta) $. If this is None then keep same.

@param tsn the sine components of

$ heta( ho, artheta, zeta) $. If this is None then keep same.

@param ssn the sine componets of

$s( ho, artheta, zeta) $. If this is None then keep same.

@param tcn the cosine components of

$ heta( ho, artheta, zeta) $. If this is None then keep same.

get_coords(sarr=[0], tarr=[0], zarr=[0], derivative=0, input1D=True)

! Compute the coordinates and their derivatives given $ ho, artheta, zeta $

@param sarr

$ ho $, an array

@param tarr

$ artheta $, an array

@param zarr

$zeta $, an array

@param derivative 0:only return

$s, heta $, 1: also return the derivatives, 2: also return the second derivatives

@param input1D if False, create a meshgrid with sarr, tarr and zarr, if True, treat them as a list of points @returns coords coords.s, coords.t contains the coordinates

$s, heta $. coords.ds and dt contains the first derivatives. coords.dds and ddt contains second derivates. coords.jacobi contains the jacobi matrix $J = partial(s, heta,zeta)/partial( ho, artheta,zeta) $, coords.djacobi contains the derivative.

jacobi_transform(jacobi, coords: CoordsOutput)

! Compute the derivatives wrt the new coordinate given the derivatives in the old coordinates

@param jacobi the Jacobi matrix

$partial(f_1, cdots, f_N) / partial (s, heta,zeta) $

@param coords the output of get_coords, should match the dimension of J @returns jacobi_output the transformed Jacobi matrix and the derivatives wrt

$( ho, artheta,zeta) $

Given

$N $ variables $f_1, cdots, f_N $, the jacobi matrix (derivatives wrt the new coordinate is given by the chain rule

[

rac{partial(f_1, cdots, f_N)}{partial( ho, artheta, zeta)} = rac{partial(f_1, cdots, f_N)}{partial(s, heta, zeta)}

cdot

rac{partial(s, heta, zeta)}{partial( ho, artheta, zeta)}

]

metric_transform(g, coords: CoordsOutput)

! Compute the coordinate transformation for a metric tensor $g_{ij} $ (lower!) (derivative computation under construction)

@param g the metric

$g_{ij} $ of the old coordinate $(s, heta,zeta) $ to be transformed, should have the dimension (3,3,…)

@param coords the output of get_coords, should match the dimension of g @returns goutput, [dgoutput], the transformed metric and the derivative wrt

$( ho, artheta,zeta) $

Let the input metric for

$(s, heta, zeta) $ being $mathbf{g} $,

the metric for

$( ho, artheta, zeta) $ is given by

[

mathbf{G} = mathbf{J}^T mathbf{g} mathbf{J},

]

where

[

mathbf{J} =

rac{partial(s, heta,zeta)}{partial( ho, artheta,zeta)}

]

is the Jacobi matrix.

contra_vector_transform(v, coords: CoordsOutput, has_jacobian=False, derivative=False, dv=None)

! Compute the coordinate transformation for a contravariant vector $v^i $ (upper!)

@param v the vector

$v^i $ to be transformed, should have the dimension (3,…)

@param coords the output of get_coords, should match the dimension of v @param has_jacobian if the given vector already contains the jacobian @param derivative if the derivatives are needed or not. If True, dv is required. @param dv the derivative of v wrt

$(s, heta,zeta) $, have the dimension (3 derivatives, 3 component,…)

@returns voutput, dvoutput, the transformed metric and the derivative wrt

$( ho, artheta,zeta) $

The vector transformation is given by

[

left(egin{array}{c} v^{

ho} v^{ artheta} v^{zeta}

end{array}

ight)

= mathbf{J}^{-1} left(egin{array}{c} v^{{s}} v^{ heta} v^{zeta} end{array}

ight),

]

where

[

mathbf{J} =

rac{partial(s, heta,zeta)}{partial( ho, artheta,zeta)}

]

is the Jacobi matrix. For computing the result we use matrix solve instead of matrix inversion.

When the derivative is needed, we note that

[

rac{partial }{partial ho}mathbf{J}

left(egin{array}{c} v^{

ho} v^{ artheta} v^{zeta}

end{array}

ight)

+mathbf{J}

rac{partial }{partial ho}left(egin{array}{c}

v^{

ho} v^{ artheta} v^{zeta}

end{array}

ight)

=

rac{partial }{partial ho}left(egin{array}{c}

v^{{s}} v^{ heta} v^{zeta} end{array}

ight),

]

therefore

[

rac{partial}{partial ho}left(egin{array}{c}

v^{

ho} v^{ artheta} v^{zeta}

end{array}

ight)

= mathbf{J}^{-1} left[

rac{partial }{partial ho}left(egin{array}{c}

v^{{s}} v^{ heta} v^{zeta} end{array}

ight)

rac{partial mathbf{J}}{partial ho}

left(egin{array}{c} v^{

ho} v^{ artheta} v^{zeta}

end{array}

ight) ight].

]

Similarly we can compute the

$ artheta $ and $zeta $ derivatives.

Now,

[

rac{partial (v^{{s}}, v^{ heta}, v^{zeta})}{partial( ho, artheta, zeta)} =

rac{partial (v^{{s}}, v^{ heta}, v^{zeta})}{partial(s, heta, zeta)}

imes

rac{partial(s, heta, zeta)}{partial( ho, artheta, zeta)}

]

Summing up

[

rac{partial ( v^{ ho}, v^{ artheta}, v^{zeta})}{partial( ho, artheta, zeta)} =

mathbf{J}^{-1} cdot

rac{partial (v^{{s}}, v^{ heta}, v^{zeta})}{partial(s, heta, zeta)} cdot mathbf{J}

• mathbf{J}^{-1} cdot (mathbf{J}_

ho, mathbf{J}_ artheta, mathbf{J}_zeta) cdot left(egin{array}{c}

v^{

ho} v^{ artheta} v^{zeta}

end{array}

ight)

]

construct_interpolant(rhosurfs=None, method='cubic_spline', **kwargs)

! Construct the interpolant between surfaces for the given method

@param rhosurfs The

$ ho $ coordinates for each surface. If input is None, the default or internally stored will be used.

@param method The method being used, can be one of [‘cubic_spline’, ‘cubic_hermite’, ‘pchip’]. The corresponding scipy class is being used. @param **kwargs Other parameters passing onto the interpolant constructor

plot(zeta=0, npoints=129, **kwargs)

!Plot the surfaces @param zeta the toroidal plane @param npoints the number of points in the plot @param **kwargs all other parameters going into plt.plot

read_surfaces_from_file(filename='data.npz', Nfp=None)

! Read the surfaces into a numpy array on disk @param filename the filename to load from @param Nfp the toroidal periodicity, if known.

write_surfaces_to_file(filename=None)

! Save the surfaces into a numpy array on disk @param filename the filename to save to