gvec.gframe

pyGVEC G-frame

This module provides functions for computing the G-Frame from a surface, cutting a surface at a given plane, change from xyz to xyz_hat coordinates and read and write the G-Frame from/to a file.

gvec.gframe.check_field_periodicity(xyz: ndarray, nfp: int, atol=1e-12)

checks if all xyz positions of the surface on a full turn, xyz[0:nz*nfp,0:nt,0:2], have the field periodicity with nfp. returns the sign of the rotation +2pi/nfp or -2pi/nfp

gvec.gframe.construct_gframe_from_surface(xyz_in: ndarray, nfp: int, name: str, tolerance_output: float = 1e-08, format: Literal['yaml', 'toml'] = 'yaml', tolerance_clean_surface: float = 0.0, impose_stell_symmetry: bool = False, theta0: float = 0.0, zeta0: float = 0.0, cutoff_gframe: int = -1, atol_field_periodicity: float = 1e-12, boundary_coefficients: bool = False, writeFiles: bool = True, logger: Logger | None = None)

Construct a G-Frame from a surface given by its cartesian points.

Parameters:

• xyz_in (ndarray, shape (nz*nfp,nt,3)) – Cartesian points xyz[0:nz*nfp,0:nt,0:2] of the surface.

• nfp (int) – Number of field periods.

• name (str) – Name of the output file.

• tolerance_output (float, optional) – Tolerance for the output surface, computes the necessary modes in X1,X2 without changing the output data. Defaults to 1e-8.

• format (str, optional) – Output parameter file format. Defaults to “yaml”.

• tolerance_clean_surface (float, optional) – Tolerance to reduce input surface resolution, computes the necessary modes for one field period, and recomputes the input surface with these modes. Defaults to 0.0.

• impose_stell_symmetry (bool, optional) – If set, imposes stellarator symmetry for the input surface. Use this with great care! Defaults to False.

• theta0 (float, optional) – First point in logical theta direction where xyz was sampled. Defaults to 0.

• zeta0 (float, optional) – First point in logical zeta direction where xyz was sampled. Defaults to 0.

• cutoff_gframe (int, optional) – Maximum mode number (>=0) to be used along the toroidal direction to construct the G-frame. Default -1 means no cutoff.

• atol_field_periodicity (float, optional) – Absolute tolerance for field-periodicity check. Default is 1e-12

• boundary_coefficients (bool, optional) – Write the boundary data as fourier coefficients in the parameter file instead of pointing to the netCDF file.

• writeFiles (bool, optional) – If True, write the GVEC parameters to file and the G-frame data to netcdf file.

• logger (logging.Logger, optional) – Logger for logging messages. If None, a new logger is created.

Returns:

• parameters (dict) – Dictionary containing the GVEC parameters.

• dict_out (dict) – Dictionary containing the G-frame data.

gvec.gframe.cut_surf(xyz, nfp, xyz0_in, N_in, B_in)

Given a surface \(\mathbf{x}(\zeta,\theta)\) on the full torus, find intersection point of lines of \(\theta=\) const with all cutting planes defined by an origin curve xyz0 and two base vectors N_in and B_in (normal of the cutting plane is then \(N \times B\)). The intersection points are projected on the base vectors, yielding the \(X^1,X^2\) coordinates in all cutting planes. Cutting planes must be given on the full torus and must be a multiple of the number of field periods.

The input surface is assumed to have the same number of field periods, and must be given on the full torus. It is described by an array of cartesian point positions \(\mathbf{x}(\zeta_i,\vartheta_j)\), assumed to be evaluated at a meshgrid of a toroidal and poloidal parameterization \((\zeta,\theta)\), excluding the periodic endpoints:

\(\zeta_i=2\pi \frac{i}{n_\zeta n_{FP}},i=0\dots,(n_\zeta n_{FP})-1,\quad \theta_j=2\pi\frac{j}{n_\theta},j=0,\dots,n_\theta-1\)

The surface can thus be evaluated using a Fourier transform.

Note: If the number of cutting planes does not match the number of toroidal points of the surface, the cutting planes are resampled to match the number of toroidal points, before cutting.

Parameters:

• xyz (cartesian coordinates of the surface on the full torus, excluding the periodic endpoints: shape is (nzeta,ntheta,3))

• nfp (number of field periods)

• xyz0_in (origin of the cutting plane, shape is (nzeta_cut*nfp,3).)

• N_in (first base vector to span the cutting plane, shape is (nzeta_cut*nfp,3))

• B_in (second base vector to span the cutting plane, shape is (nzeta_cut*nfp,3).)

Returns:

coordinates of the intersection points for cutting planes in one field period. shape is (nzeta,ntheta).

Return type:

x1_cut,x2_cut

gvec.gframe.eval_curve(zeta_in, xyz, dft_dict)

evaluates the curve at a single point zeta_in given by cartesian positions of a periodic curve xyz[0:len(zeta1d)+1,0:2], evaluated zeta1d[0:2pi[,

gvec.gframe.eval_distance_to_curve(zeta_in, origin, normal, xyz, dft_dict)

gvec.gframe.eval_distance_to_plane(xyz_in, origin, normal)

gvec.gframe.find_zeta_cuts(zeta_in, origin, normal, xyz, dft_dict, zeta_bracket: float)

gvec.gframe.frenet_frame(X0: ndarray, nzeta: int | None = None) → dict

Compute the Frenet frame of a closed (!) curve defined by \(X_0(\zeta)\), by computing the tangent vector \(T\), normal vector \(N\) and binormal vector \(B\).

Derivatives of \(X_0\) are computed via dft, so the curve is assumed to be given on an equispaced grid in zeta=zeta0+np.linspace(0,2*np.pi,npoints,endpoint=False), excluding the endpoint.

Warning

The curvature kappa scales with the second derivative. If it is zero at any point, normal and bi-normal vectors are not defined. N and B and torsion tau will return as None.

Parameters:

• X0 (np.ndarray) – closed curve positions, in cartesian coordinates, shape (npoints,3), excluding periodic endpoint!

• nzeta (int, optional) – number of points for the output, default is None, which then uses npoints.

Returns:

dict_frame – dictionary with:

• X0: input curve positions, shape (npoints,3)

• T: tangent vector, shape (npoints,3)

• N: normal vector, shape (npoints,3). If curvature is zero at any point, None is returned

• B: binormal vector, shape (npoints,3). If curvature is zero at any point, None is returned

• lp: arc-length element of the curve \(|X_0^{\prime}(\zeta)|\), shape (npoints,)

• kappa: curvature, shape (npoints,)

• tau: torsion, shape (npoints,). If curvature is zero at any point, None is returned

Return type:

dict

gvec.gframe.frenet_frame_evaluate(X0, X0p, X0pp, X0ppp) → dict

Compute the Frenet frame from the first three derivatives of the curve \(X_0(\zeta)\) in cartesian coordinates.

The tangent vector, binormal vector and normal vector are computed as

\[T = \frac{X_0^\prime}{|X_0^\prime|}, \quad B = \frac{X_0^\prime \times X_0^{\prime\prime}}{|X_0^\prime \times X_0^{\prime\prime}|}, \quad N = B \times T\]

The arc-length element is \(\ell^\prime = |X_0^\prime|\), and curvature and torsion are computed as

\[\kappa = \frac{|X_0^\prime \times X_0^{\prime\prime}|}{(\ell^\prime)^3},\quad \tau = \frac{((X_0^\prime \times X_0^{\prime\prime}) \cdot X_0^{\prime\prime\prime})}{|X_0^\prime \times X_0^{\prime\prime}|^2}\,.\]

Warning

The curvature kappa scales with the second derivative. If it is zero at any point, normal and bi-normal vectors are not defined. N and B and torsion tau will return as None.

Parameters:

• X0 (np.ndarray) – closed curve positions, in cartesian coordinates, shape (npoints,3), excluding periodic endpoint!

• X0p (np.ndarray) – first derivative of the curve, shape (npoints, 3)

• X0pp (np.ndarray) – second derivative of the curve, shape (npoints, 3)

• X0ppp (np.ndarray) – third derivative of the curve, shape (npoints, 3)

Returns:

dict_frame – dictionary with:

• T: tangent vector, shape (npoints,3)

• N: normal vector, shape (npoints,3). If any kappa*lp <1e-8`, ``None is returned

• B: binormal vector, shape (npoints,3). If any kappa*lp <1e-8`, ``None is returned

• lp: arc-length element of the curve \(|X_0^\prime(\zeta)|\), shape (npoints,)

• kappa: curvature, shape (npoints,)

• tau: torsion, shape (npoints,). If any kappa*lp <1e-8`, ``None is returned

Return type:

dict

gvec.gframe.get_X0_N_B(xyz)

Get guiding curve and two guiding vectors from the cartesian coordinates of a surface.

gvec.gframe.get_xyz_cut(zeta_start, origins, normals, xyz_in, dft_dict, nfp)

Find the xyz positions of the cut of a surface with a given number of planes.

Parameters:

• zeta_start (float, shape (nz_out,)) – Initial guess for zeta position for each cut.

• origins (float, shape (nz_out, 3)) – Origin of the cutting plane.

• normals (float, shape (nz_out, 3)) – Normal of the cutting plane.

• dft_dict (dict) – Dictionary containing the dft from the zeta points of xyz_in, to be used for evaluation of theta=const curves at arbitrary zeta.

Returns:

xyz_cut – XYZ positions from the cuts.

Return type:

float, shape (nz_out, nt, 3)

gvec.gframe.minimal_modes(X, Y, Z=None, tolerance=1e-08)

Find the minimal maximum mode numbers (M, N) such that the error is below the tolerance. First dimension of X and Y is assumed to be theta (starting at 0., without endpoint), second dimension is assumed to be zeta (starting at 0., without endpoint).

gvec.gframe.plot_cross_section_comparison(dict_surf, dict_RZ, step=1, halfperiod=True)

Plot cross-sections from two dictionaries

Parameters:

• dict_surf (dict) – Dictionary from to_surface function

• dict_RZ (dict) – Dictionary from to_RZ function

• step (int, optional) – Step in zeta array. Defaults to 1.

• halfperiod (bool, optional) – If True, only half of the zeta array is plotted. Defaults to True.

Returns:

fig – Figure object

Return type:

matplotlib.figure.Figure

gvec.gframe.read_Gframe_ncfile(ncfile: str | Path)

Read G-frame netcdf file and store data in a dictionary.

Parameters:

ncfile (str or Path) – Name/path to netcdf file

Returns:

dict_out – Dictionary with the data (with ‘axis’ and ‘boundary’ groups, if they exist in the netcdf file)

Return type:

dict

gvec.gframe.rodrigues(xyz: ndarray, angle, origin=array([0., 0., 0.]), rot_axis=array([0., 0., 1.]))

Rodrigues rotation function.

Parameters:

• xyz (ndarray) – cartesian point positions. if multidimensional, the LAST dimension must contain the cartesian components.

• angle (float) – The rotation angle.

Returns:

pos_rot – The rotated position vector (cartesian components).

Return type:

ndarray

gvec.gframe.to_RZ(xyz: ndarray, nfp: int, nzeta=81, ntheta=81, zeta0: float = 0.0, theta0: float = 0.0, tolerance: float = 1e-08)

Cut a xyz surface to yield a R,Z positions on one field period.

Parameters:

• xyz (ndarray) – Boundary surface in cartesian coordinates, with shape (nzeta*nfp,ntheta,3)

• nfp (int) – Number of field periods

• zeta0 (float, optional) – First point in logical zeta direction where xyz was sampled. Defaults to 0.

• theta0 (float, optional) – First point in logical theta direction where xyz was sampled. Defaults to 0.

• nzeta (int, optional) – Number of zeta positions (=geometric angle -phi) for the output, to sample on one field period. Defaults to 81.

• ntheta (int, optional) – Number of theta positions for the output. Defaults to 81.

• tolerance (float, optional) – Tolerance for finding minimal mode numbers. Defaults to 1e-8.

Returns:

Dictionary with: - zeta : zeta positions on one field period, length nzeta - theta : theta positions, length ntheta - R : R positions on one field period, with shape (ntheta,nzeta) - Z : Z positions on one field period, with shape (ntheta,nzeta) - nfp : number of field periods - lasym : logical for asymmetry, =False if stellarator symmetry is found - Mmax,``Nmax`` : maximum mode numbers needed for the given tolerance - Rc,``Rs``,``Zc``,``Zs`` : R and Z cosine and sine Fourier mode coefficients, shape is (Mmax+1,2*Nmax+1) - m_modes: poloidal mode numbers (m) in first dimension of Rc, Rs,Zc, Zs - n_modes: toroidal mode numbers (n) in second dimension of Rc, Rs,Zc, Zs - tolerance: input tolerance

Return type:

dict

gvec.gframe.to_axis(dict_in: dict, nzeta: int = 81)

Convert the “axis” of a gframe file to a different resolution.

Parameters:

• dict_in (dict) – dictionary of the Gframe netcdf file, from gvec.gframe.read_Gframe_ncfile(filename)

• nzeta (int, optional) – number of zeta positions for the output, to sample on one field period (default: 81)

Returns:

dict_axis – dictionary containing

• nfp: number of field periods

• axis: dictionary for the G-Frame data:

  • nzeta: number of zeta positions for one field-period

  • zetafull: equidistant positions along the curve parameter zeta, without endpoint, length is nfp*nzeta

  • xyz: cartesian positions along the curve, shape is (3,nfp*nzeta)

  • Nxyz: cartesian components of the ‘normal’ vector, shape is (3,nfp*nzeta)

  • Bxyz: cartesian components of the ‘bi-normal’ vector, shape is (3,nfp*nzeta)

Return type:

dict

gvec.gframe.to_surface(dict_in: dict, nzeta: int = 81, ntheta: int = 81, tolerance: float = 1e-08)

Convert a gframe file with axis+boundary to boundary surface in cartesian coordinates.

Parameters:

• dict_in (dict) – dictionary of the Gframe netcdf file, from gvec.gframe.read_Gframe_ncfile(filename)

• nzeta (int, optional) – number of zeta positions for the output, to sample on one field period (default: 81)

• ntheta (int, optional) – number of theta positions for the output (default: 81)

Returns:

dictionary with:

• xyz : boundary surface in cartesian coordinates, with shape (nzeta*nfp,ntheta,3)

• X1, X2 : boundary in G-Frame, shape (ntheta,nzeta)

• zetafull : zeta values of the boundary surface

• theta : theta values of the boundary surface

• lasym : logical for asymmetry, =False if stellarator symmetry is found

• nfp : number of field periods

• Mmax, Nmax : maximum mode numbers needed for the given tolerance

• X1c, X1s, X2c, X2s : boundary modes in G-Frame, up to Mmax, Nmax

• m_modes: poloidal mode numbers (m) in first dimension of Rc, Rs,Zc, Zs

• n_modes: toroidal mode numbers (n) in second dimension of Rc, Rs,Zc, Zs

Return type:

dict

gvec.gframe.twist_of_ribbon(X0: ndarray, N: ndarray, nint: int | None = None, return_integrand: bool = False) → float

Compute the twist of a closed (!) ribbon defined by a centerline curve \(X_0(\zeta)\) and non-vanishing “normal” vector \(N(\zeta)\). The vector N only needs to be linearly independent of the tangent vector \(|X_0^\prime(\zeta) \times N| >0\)

The twist is computed as

\[ \text{Tw} = \frac{1}{2\pi}\int_0^{2\pi}\frac{\left ({N\times \left [{N^\prime\left|{X_0^\prime}\right|^2 - \left({N \cdot X_0^\prime}\right) X_0^{\prime\prime}}\right]}\right) \cdot X_0^\prime}{\left|{N\left|{X_0^\prime}\right|^2-\left({N \cdot X_0^\prime}\right) X_0^\prime}\right|^2}\left|{X_0^\prime}\right| d\zeta \]

Derivatives of \(X_0\) and \(N\) are computed via fft, so the curve is assumed to be given on an equispaced grid in zeta=np.linspace(0,2*np.pi,npoints,endpoint=False), excluding the endpoint.

Parameters:

• X0 (np.ndarray) – positions along closed curve, in cartesian coordinates, shape (npoints,3), excluding periodic endpoint!

• N (np.ndarray) – linearly independent, “normal” vector at curve positions, in cartesian coordinates, same shape as X0. Does not need to be normalized.

• nint (int, optional) – number of itnegration points, None sets default =npoints. Only used if nint > npoints, to upsample X0 and N for more accurate integration.

• return_integrand (bool, optional) – If True, the point values of the integrand of the twist integral is returned. If False, the integral result, np.average(integrand) is returned. Default is False.

Returns:

Tw – twist of the ribbon defined by X0 and N

Return type:

float

gvec.gframe.write_Gframe_ncfile(filename: str | Path, dict_in)

Write the G-Frame [and boundary] to a GVEC-compatible netCDF file.

Parameters:

• filename (str | Path) – name of the file to be written, should include .nc ending, overwrites existing file

• dict_in (dict) –

  dictionary containing

  • nfp: number of field periods

  • axis: dictionary for the G-Frame data:

    • nzeta: number of zeta positions for one field-period

    • zetafull: equidistant positions along the curve parameter zeta, without endpoint, length is nfp*nzeta

    • xyz: cartesian positions along the curve, shape is (3,nfp*nzeta)

    • Nxyz: cartesian components of the ‘normal’ vector, shape is (3,nfp*nzeta)

    • Bxyz: cartesian components of the ‘bi-normal’ vector, shape is (3,nfp*nzeta)

  • boundary: dictionary for the boundary data (written to file if exists):

    • nzeta: number of (toroidal ) zeta positions for one field-period

    • zeta: zeta positions where boundary was sampled, length is nzeta

    • n_max: maximum toroidal mode number ``>=(nzeta-1)//2

    • ntheta: number of (poloidal) theta positions for one field-period

    • theta: theta positions where boundary was sampled, length is ntheta

    • m_max: maximum poloidal mode number ``>=(ntheta-1)//2

    • X1: position of the boundary in the G-Frame, first component, in direction of N, shape is (ntheta,nzeta)

    • X2: position of the boundary in the G-Frame, second component, in direction of B, shape is (ntheta,nzeta)

    • lasym: logical for asymmetry, =False if boundary is stellarator-symmetric

Return type:

None

gvec.gframe.writhe(X0: ndarray, N: ndarray | None = None, width: float = 1e-06, nint: int | None = None)

Compute the writhe of a closed (!) curve defined by \(X_0(\zeta)\), by attaching a ribbon defined along a non-vanishing “normal” vector \(N(\zeta)\), and computing writhe as the difference of the linking number of the ribbon and the twist of the ribbon, Wr = Lk - Tw`. The vector ``N only needs to be linearly independent of the tangent vector \(|X_0^\prime(\zeta) \times N| >0\).

If N is not provided, the centroid frame of the curve is used.

Note that this method is more accurate for a smooth curve than using a polygonal approximation writhe_from_polygon

Parameters:

• X0 (np.ndarray) – positions along closed curve, in cartesian coordinates, shape (npoints,3), excluding periodic endpoint!

• N (np.ndarray, optional) – linearly independent, non-normalized “normal” vector at curve positions, in cartesian coordinates, shape (npoints,3). If not provided, N is computed from the centroid frame of X0

• width (float, optional) – In order to compute the linking number, the width of the ribbon must be specified, generating a second curve X0+N*width. Default width is 1e-6

• nint (int, optional) – number of points for integration for twist, None sets default =npoints

Returns:

• Wr (float) – Write of the closed curve X0

• Lk (float) – Linking number of the ribbon defined by X0 and N

• Tw (float) – twist of the ribbon defined by X0 and N

gvec.gframe.xyz_hat_to_xyz(xhat: ndarray, yhat: ndarray, zhat: ndarray, zeta: ndarray, sign_rot: float)

Change from xyz ‘hat’ coordinates to cartesian xyz positions.

Parameters:

• xhat (ndarray) – Hat coordinates, size [0:nz,0:nt,0:2].

• yhat (ndarray) – Hat coordinates, size [0:nz,0:nt,0:2].

• zhat (ndarray) – Hat coordinates, size [0:nz,0:nt,0:2].

• zeta (ndarray) – 1d array of zeta positions corresponding to xyz positions (without endpoint), size [0:nz].

• sign_rot (float) – Direction of zeta for the rotation into hat coordinates, +1 or -1.

Returns:

xyz – Cartesian xyz positions of the surface, size [0:nz,0:nt,0:2].

Return type:

ndarray

gvec.gframe.xyz_to_xyz_hat(xyz_in: ndarray, zeta: ndarray, sign_rot: float)

Change from cartesian xyz positions of the surface (can be a full torus or a single field period) on to ‘hat’ coordinates, which are periodic on field period.

Parameters:

• xyz_in (ndarray) – xyz positions of the surface, xyz[0:nz,0:nt,0:2], sampled at zeta positions, must exclude the endpoint

• zeta (ndarray) – 1d array of zeta positions belonging to the surface (without endpoint), size [0:nz]

• sign_rot (float) – direction of zeta for the rotation into hat coordinates, +1 or -1

Returns:

xhat, yhat, zhat – hat coordinates, periodic on one field period, size [0:nz,0:nt,0:2]

Return type:

ndarray