surf

Module Contents

class surf.Surface(x0, external_dof_setter, names)

Surface is a base class for various representations of toroidal surfaces.

A Surface is modelled as a function \(\Gamma:[0, 1] \times [0, 1] \to R^3\) and is evaluated at quadrature points \(\{\phi_1, \ldots, \phi_{n_\phi}\}\times\{\theta_1, \ldots, \theta_{n_\theta}\}\).

RANGE_FULL_TORUS = 'full torus'
RANGE_FIELD_PERIOD = 'field period'
RANGE_HALF_PERIOD = 'half period'
classmethod from_nphi_ntheta(nphi=61, ntheta=62, range='full torus', nfp=1, **kwargs)

Initializes surface classes from the specified number of grid points along toroidal, \(\phi\), and poloidal, \(\theta\), directions. Additional parameters required for surface initialization could be supplied as keyword arguments.

Args:
nphi: Number of grid points \(\phi_j\) in the toroidal angle

\(\phi\).

ntheta: Number of grid points \(\theta_i\) in the poloidal angle

\(\theta\).

range: Toroidal extent of the \(\phi\) grid.

Set to "full torus" (or equivalently SurfaceRZFourier.RANGE_FULL_TORUS) to generate quadrature points up to 1 (with no point at 1). Set to "field period" (or equivalently SurfaceRZFourier.RANGE_FIELD_PERIOD) to generate points up to \(1/n_{fp}\) (with no point at \(1/n_{fp}\)). Set to "half period" (or equivalently SurfaceRZFourier.RANGE_HALF_PERIOD) to generate points up to \(1/(2 n_{fp})\), with all grid points shifted by half of the grid spacing in order to provide spectral convergence of integrals.

nfp: The number of field periods. kwargs: Additional arguments to initialize the surface classes. Look

at the docstrings of the specific class you are interested in.

static get_quadpoints(nphi=None, ntheta=None, range=None, nfp=1)

Sets the theta and phi grid points for Surface subclasses. It is typically called in when constructing Surface subclasses.

For more information about the arguments nphi, ntheta, range, quadpoints_phi, and quadpoints_theta, see the general documentation on surfaces.

Args:

nfp: The number of field periods. nphi: Number of grid points \(\phi_j\) in the toroidal angle \(\phi\). ntheta: Number of grid points \(\theta_j\) in the toroidal angle \(\theta\). range: Toroidal extent of the \(\phi\) grid.

Set to "full torus" (or equivalently Surface.RANGE_FULL_TORUS) to generate points up to 1 (with no point at 1). Set to "field period" (or equivalently Surface.RANGE_FIELD_PERIOD) to generate points up to \(1/n_{fp}\) (with no point at \(1/n_{fp}\)). Set to "half period" (or equivalently Surface.RANGE_HALF_PERIOD) to generate points up to \(1/(2 n_{fp})\), with all grid points shifted by half of the grid spacing in order to provide spectral convergence of integrals.

Returns:

Tuple containing

  • quadpoints_phi: List of grid points \(\phi_j\).

  • quadpoints_theta: List of grid points \(\theta_j\).

static get_theta_quadpoints(ntheta=None)

Sets the theta grid points for Surface subclasses.

Args:

ntheta: Number of grid points \(\theta_j\) in the toroidal angle \(\theta\).

Returns:

List of grid points \(\theta_j\).

static get_phi_quadpoints(nphi=None, range=None, nfp=1)

Sets the phi grid points for Surface subclasses.

Args:

nphi: Number of grid points \(\phi_j\) in the toroidal angle \(\phi\). range: Toroidal extent of the \(\phi\) grid.

Set to "full torus" (or equivalently Surface.RANGE_FULL_TORUS) to generate points up to 1 (with no point at 1). Set to "field period" (or equivalently Surface.RANGE_FIELD_PERIOD) to generate points up to \(1/n_{fp}\) (with no point at \(1/n_{fp}\)). Set to "half period" (or equivalently Surface.RANGE_HALF_PERIOD) to generate points up to \(1/(2 n_{fp})\), with all grid points shifted by half of the grid spacing in order to provide spectral convergence of integrals.

nfp: The number of field periods.

Returns:

List of grid points \(\phi_j\).

plot(engine='matplotlib', ax=None, show=True, close=False, axis_equal=True, plot_normal=False, plot_derivative=False, wireframe=True, **kwargs)

Plot the surface in 3D using matplotlib/mayavi/plotly.

Args:
engine: Selects the graphics engine. Currently supported options are "matplotlib" (default),

"mayavi", and "plotly".

ax: The figure/axis to be plotted on. This argument is useful when plotting multiple

objects on the same axes. If equal to the default None, a new axis will be created.

show: Whether to call the show() function of the graphics engine.

Should be set to False if more objects will be plotted on the same axes.

close: Whether to close the seams in the surface where the angles jump back to 0. axis_equal: For matplotlib, whether to adjust the scales of the x, y, and z axes so

distances in each direction appear equal.

plot_normal: Whether to plot the surface normal vectors. Only implemented for mayavi. plot_derivative: Whether to plot the surface derivatives. Only implemented for mayavi. wireframe: Whether to plot the wireframe in Mayavi. kwargs: Any additional arguments to pass to the plotting function, like color='r'.

Note: the ax and show parameters can be used to plot more than one surface:

ax = surface1.plot(show=False)
ax = surface2.plot(ax=ax, show=False)
surface3.plot(ax=ax, show=True)
Returns:

An axis which could be passed to a further call to the graphics engine so multiple objects are shown together.

abstractmethod to_RZFourier()

Return a SurfaceRZFourier instance corresponding to the shape of this surface. All subclasses should implement this abstract method.

cross_section(phi, thetas=None, tol=1e-13)

Computes an array of points on the cross section with cylindrical angle \(\phi\), using the bisection method. The poloidal angles of the points on the cross section are given by thetas. This function assumes that the surface intersection with the \(\phi\)-plane is a single curve: the surface should only go once around the z-axis, and it should not go back on itself.

Args:
phi (float):

The standard cylindrical angle (toroidal angle) normalized by \(2\pi\), i.e. \(\phi=0, 1\) corresponds to the standard cylindrical angle \(0, 2\pi\), respectively.

thetas (int, array, optional):

An optional argument indicating at which poloidal angle the cross section should be calculated. If thetas is an int, the cross section will be calculated at the poloidal angles in the array np.linspace(0, 1, thetas, endpoint=False). thetas can also be an array of poloidal angles between \(0\) and \(1\). If thetas is not provided, then the cross section will be calculated at the positions given by Surface.quadpoints_theta. Defaults to None.

tol (float):

The tolerance for the bisection root-finding. Defaults to 1e-13.

Returns:

(ntheta, 3) array of the Cartesian coordinates along the cross-section at each of the ntheta poloidal angles.

cross_section_phi(zetan=0.0, thetas=200)
aspect_ratio()

Note: cylindrical coordinates are \((R, \phi, Z)\), where \(\phi \in [-\pi,\pi)\) and the angles that parametrize the surface are \((\varphi, \theta) \in [0,1)^2\) For a given surface, this function computes its aspect ratio using the VMEC definition:

\[AR = R_{\text{major}} / R_{\text{minor}}\]

where

\[\begin{split}R_{\text{minor}} &= \sqrt{ \overline{A} / \pi } \\ R_{\text{major}} &= \frac{V}{2 \pi^2 R_{\text{minor}}^2}\end{split}\]

and \(V\) is the volume enclosed by the surface, and \(\overline{A}\) is the average cross sectional area.

daspect_ratio_by_dcoeff()

Return the derivative of the aspect ratio with respect to the surface coefficients

d2aspect_ratio_by_dcoeff_dcoeff()

Return the derivative of the aspect ratio with respect to the surface coefficients

minor_radius()

Return the minor radius of the surface using the formula

\[R_{\text{minor}} = \sqrt{ \overline{A} / \pi }\]

where \(\overline{A}\) is the average cross sectional area.

dminor_radius_by_dcoeff()

Return the derivative of the minor radius wrt surface coefficients

d2minor_radius_by_dcoeff_dcoeff()

Return the second derivative of the minor radius wrt surface coefficients

major_radius()

Return the major radius of the surface using the formula

\[R_{\text{major}} = \frac{V}{2 \pi^2 R_{\text{minor}}^2}\]

where \(\overline{A}\) is the average cross sectional area, and \(R_{\text{minor}}\) is the minor radius of the surface.

dmajor_radius_by_dcoeff()

Return the derivative of the major radius wrt surface coefficients

d2major_radius_by_dcoeff_dcoeff()

Return the second derivative of the major radius wrt surface coefficients

mean_cross_sectional_area()

Note: cylindrical coordinates are \((R, \phi, Z)\), where \(\phi \in [-\pi,\pi)\) and the angles that parametrize the surface are \((\varphi, \theta) \in [0,1)^2\). The mean cross sectional area is given by the integral

\[\overline{A} = \frac{1}{2\pi} \int_{S_{\phi}} ~dS ~d\phi\]

where \(S_\phi\) is the cross section of the surface at the cylindrical angle \(\phi\). Note that \(\int_{S_\phi} ~dS\) can be rewritten as a line integral

\[\begin{split}\int_{S_\phi}~dS &= \int_{S_\phi} ~dR dZ \\ &= \int_{\partial S_\phi} [R,0] \cdot \mathbf n/\|\mathbf n\| ~dl \\ &= \int^1_{0} R \frac{\partial Z}{\partial \theta}~d\theta\end{split}\]

where \(\mathbf n = [n_R, n_Z] = [\partial Z/\partial \theta, -\partial R/\partial \theta]\) is the outward pointing normal. Consider the surface in cylindrical coordinates terms of its angles \([R(\varphi,\theta), \phi(\varphi,\theta), Z(\varphi,\theta)]\). The boundary of the cross section \(\partial S_\phi\) is given by the points \(\theta\rightarrow[R(\varphi(\phi,\theta),\theta),\phi, Z(\varphi(\phi,\theta),\theta)]\) for fixed \(\phi\). The cross sectional area of \(S_\phi\) becomes

\[\int^{1}_{0} R(\varphi(\phi,\theta),\theta) \frac{\partial}{\partial \theta}[Z(\varphi(\phi,\theta),\theta)] ~d\theta\]

Now, substituting this into the formula for the mean cross sectional area, we have

\[\overline{A} = \frac{1}{2\pi}\int^{\pi}_{-\pi}\int^{1}_{0} R(\varphi(\phi,\theta),\theta) \frac{\partial}{\partial \theta}[Z(\varphi(\phi,\theta),\theta)] ~d\theta ~d\phi\]

Instead of integrating over cylindrical \(\phi\), let’s complete the change of variables and integrate over \(\varphi\) using the mapping:

\[[\phi,\theta] \leftarrow [\text{atan2}(y(\varphi,\theta), x(\varphi,\theta)), \theta]\]

After the change of variables, the integral becomes:

\[\overline{A} = \frac{1}{2\pi}\int^{1}_{0}\int^{1}_{0} R(\varphi,\theta) \left[\frac{\partial Z}{\partial \varphi} \frac{\partial \varphi}{d \theta} + \frac{\partial Z}{\partial \theta} \right] \text{det} J ~d\theta ~d\varphi\]

where \(\text{det}J\) is the determinant of the mapping’s Jacobian.

dmean_cross_sectional_area_by_dcoeff()

Return the derivative of the mean cross sectional area wrt surface coefficients

d2mean_cross_sectional_area_by_dcoeff_dcoeff()

Return the second derivative of the mean cross sectional area wrt surface coefficients

arclength_poloidal_angle()

Computes a poloidal (angle) coordinate θ on a surface for which the arclength ∂|r|/∂θ is independent of θ in each φ plane. In other words, this function computes the uniform-arclength poloidal coordinate. The returned poloidal coordinate is in the range [0,1), and is used in methods evaluating the adjoint shape gradient.

Returns:

2d array of shape (numquadpoints_phi, numquadpoints_theta) containing the new poloidal angle

interpolate_on_arclength_grid(function, theta_evaluate)

Interpolate function onto the theta_evaluate grid in the arclength poloidal angle. This is required for evaluating the adjoint shape gradient for free-boundary calculations.

The theta_evaluate grid may have a different number of poloidal grid points compared to the surface’s numquadpoints_theta, but it must have the same number of toroidal grid points as the surface’s numquadpoints_phi. This is because we interpolate in theta but not phi.

Returns:
function_interpolated: 2d array (numquadpoints_phi,numquadpoints_theta)

defining interpolated function on arclength angle along curve at constant phi

make_theta_uniform_arclength()

Reparameterize the surface in terms of a uniform-arclength poloidal angle.

To do the conversion accurately, make sure the surface has both a sufficient number of quadrature points and a sufficiently large number of basis functions. More basis functions may be needed to represent the shape than for the original theta coordinate.

property deduced_range

The quadpoints of a surface can be anything, but are often set to ‘full torus’, ‘field period’ or ‘half period’. Since this is not stored in the object, but often useful to know this function deduces the range from the quadpoints