simsopt.geo package

class simsopt.geo.ArclengthVariation(curve, nintervals='full')

Bases: Optimizable

J()

__init__(curve, nintervals='full')

This class penalizes variation of the arclength along a curve. The idea of this class is to avoid ill-posedness of curve objectives due to non-uniqueness of the underlying parametrization. Essentially we want to achieve constant arclength along the curve. Since we can not expect perfectly constant arclength along the entire curve, this class has some support to relax this notion. Consider a partition of the \([0, 1]\) interval into intervals \(\{I_i\}_{i=1}^L\), and tenote the average incremental arclength on interval \(I_i\) by \(\ell_i\). This objective then penalises the variance

\[J = \mathrm{Var}(\ell_i)\]

it remains to choose the number of intervals \(L\) that \([0, 1]\) is split into. If nintervals="full", then the number of intervals \(L\) is equal to the number of quadrature points of the curve. If nintervals="partial", then the argument is as follows:

A curve in 3d space is defined uniquely by an initial point, an initial direction, and the arclength, curvature, and torsion along the curve. For a simsopt.geo.curvexyzfourier.CurveXYZFourier, the intuition is now as follows: assuming that the curve has order \(p\), that means we have \(3*(2p+1)\) degrees of freedom in total. Assuming that three each are required for both the initial position and direction, \(6p-3\) are left over for curvature, torsion, and arclength. We want to fix the arclength, so we can afford \(2p-1\) constraints, which corresponds to \(L=2p\).

Finally, the user can also provide an integer value for nintervals and thus specify the number of intervals directly.

dJ(*args, partials=False, **kwargs)

return_fn_map: Dict[str, Callable] = {'J': <function ArclengthVariation.J>, 'dJ': <function derivative_dec.<locals>._derivative_dec>}

class simsopt.geo.Area(surface, range=None, nphi=None, ntheta=None)

Bases: Optimizable

Wrapper class for surface area label.

J()

Compute the area of a surface.

d2J_by_dsurfacecoefficientsdsurfacecoefficients()

Calculate the second partial derivatives with respect to the surface coefficients.

dJ(*args, partials=False, **kwargs)

dJ_by_dsurfacecoefficients()

Calculate the partial derivatives with respect to the surface coefficients.

class simsopt.geo.AspectRatio(surface, range=None, nphi=None, ntheta=None)

Bases: Optimizable

Wrapper class for surface aspect ratio.

J()

Compute the aspect ratio of a surface.

d2J_by_dsurfacecoefficientsdsurfacecoefficients()

Calculate the second partial derivatives with respect to the surface coefficients.

dJ(*args, partials=False, **kwargs)

dJ_by_dsurfacecoefficients()

Calculate the partial derivatives with respect to the surface coefficients.

class simsopt.geo.BoozerResidual(boozer_surface, bs)

Bases: Optimizable

This term returns the Boozer residual penalty term

\[J = \int_0^{1/n_{\text{fp}}} \int_0^1 \| \mathbf r \|^2 ~d\theta ~d\varphi + w (\text{label.J()-boozer_surface.constraint_weight})^2.\]

where

\[\mathbf r = \frac{1}{\|\mathbf B\|}[G\mathbf B_\text{BS}(\mathbf x) - ||\mathbf B_\text{BS}(\mathbf x)||^2 (\mathbf x_\varphi + \iota \mathbf x_\theta)]\]

J()

Return the value of the penalty function.

compute()

dJ(*args, partials=False, **kwargs)

dJ_by_dB()

Return the partial derivative of the objective with respect to the magnetic field

recompute_bell(parent=None)

Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.

This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.

Need to be implemented by classes that provide a dof_setter for external handling of DOFs.

class simsopt.geo.BoozerSurface(biotsavart, surface, label, targetlabel, constraint_weight=None, options=None)

Bases: Optimizable

The BoozerSurface class computes a flux surface of a BiotSavart magnetic field where the angles of the surface are Boozer angles [1,2]. The class takes as input a Surface representation (SurfaceXYZFourier or SurfaceXYZTensorFourier), a BiotSavart magnetic field, a flux surface label evaluator, and a target value of the label.

The Boozer angles are computed by solving a constrained least squares problem,

\[\min_x J(x) = \frac{1}{2} \mathbf r^T(x) \mathbf r(x)\]

subject to

\[ \begin{align}\begin{aligned}l(x) = l_0\\z(\varphi=0,\theta=0) = 0\end{aligned}\end{align} \]

where \(\mathbf r\) is a vector of residuals computed by boozer_surface_residual, \(l\) is a surface label function with target value \(l_0\). The degrees of freedom are the surface coefficients, the rotational transform, \(\iota\), and the value of Boozer’s \(G\) on the surface. This objective is zero when the surface corresponds to a magnetic surface of the field, \((\phi,\theta)\) that parametrize the surface correspond to Boozer angles, and the constraints are satisfied.

The recommended approach to finding the Boozer angles is to use the run_code method,

run_code(iota_guess, G=G_guess).

Depending on how the class is initialized, run_code, will use either the BoozerLS [2] or BoozerExact [1] approach to finding the flux surface. The BoozerLS approach finds the flux surface by solving the constrained least squares problem mentioned above. The methods

1. minimize_boozer_penalty_constraints_LBFGS

2. minimize_boozer_penalty_constraints_newton

3. minimize_boozer_penalty_constraints_ls

scalarize the constrained problem using a quadratic penalty method, and apply L-BFGS, Newton, or scipy.optimize.least_squares to solve the penalty problem. Alternatively, the constraints can be enforced exactly (not with a penalty) using,

minimize_boozer_exact_constraints_newton

In this approach, Newton’s method is used to solve the first order necessary conditions for optimality. Note that this differs from the BoozerExact approach.The BoozerExact approach solves the residual equations directly at a specific set of colocation points on the surface,

\[ \begin{align}\begin{aligned}\mathbf r(x) = 0\\l(x) = l_0\\z(\varphi=0,\theta=0) = 0\end{aligned}\end{align} \]

The colocation points are chosen such that the number of colocation points is equal to the number of unknowns in on the surface, so that the resulting nonlinear system of equations can be solved using Newton’s method. The BoozerExact approach is implemented in

solve_residual_equation_exactly_newton

Generally, the BoozerExact approach is faster than the BoozerLS approach, but it is less robust. Note that there are specific requirements on the set of colocation points, i.e. surface.quadpoints_phi and surface.quadpoints_theta, for stellarator symmetric BoozerExact surfaces. See the class method solve_residual_equation_exactly_newton and get_stellsym_mask() for more information.

[1]: Giuliani A, Wechsung F, Stadler G, Cerfon A, Landreman M. Direct computation of magnetic surfaces in Boozer coordinates and coil optimization for quasisymmetry. Journal of Plasma Physics. 2022;88(4):905880401. doi:10.1017/S0022377822000563

[2]: Giuliani, A., Wechsung, F., Cerfon, A., Landreman, M., & Stadler, G. (2023). Direct stellarator coil optimization for nested magnetic surfaces with precise quasi-symmetry. Physics of Plasmas, 30(4).

__init__(biotsavart, surface, label, targetlabel, constraint_weight=None, options=None)

Parameters:

• biotsavart (BiotSavart) – BiotSavart object.

• surface (SurfaceXYZFourier, SurfaceXYZTensorFourier) – Surface object.

• label (Optimizable) – A method that computes a flux surface label for the surface, such as Volume, Area, or ToroidalFlux.

• targetlabel (float) – The target value of the label on the surface.

• constraint_weight (float, Optional) – The weight of the label constraint used when solving Boozer least squares. If None, then Boozer Exact is used in the run_code method.

• options (dict, Optional) –

  A dictionary of solver options. If a keyword is not specified, then a default value is used. Possible keywords are:

  • verbose (bool): display convergence information. Defaults to True.

  • newton_tol (float): tolerance for newton solver. Defaults to 1e-13 for BoozerExact and 1e-11 for BoozerLS.

  • bfgs_tol (float): tolerance for bfgs solver. Defaults to 1e-10.

  • newton_maxiter (int): maximum number of iterations for Newton solver. Defaults to 40.

  • bfgs_maxiter (int): maximum number of iterations for BFGS solver. Defaults to 1500.

  • limited_memory (bool): True if L-BFGS solver is desired, False if the BFGS solver otherwise. Defaults to False.

  • weight_inv_modB (float): for BoozerLS surfaces, weight the residual by modB so that it does not scale with coil currents. Defaults to True.

_get_residual_vector_and_jacobian(x, constraint_weight, optimize_G, weight_inv_modB)

Helper function to get residual vector and Jacobian for least_squares

boozer_exact_constraints(xl, derivatives=0, optimize_G=True)

This function returns the optimality conditions corresponding to the minimization problem

\[\text{min}_x ~J(x)\]

subject to

\[\begin{split}l - l_0 &= 0 \\ z(\varphi=0,\theta=0) - 0 &= 0\end{split}\]

The function can additionally return the first derivatives of these optimality conditions.

Parameters:

• xl (ndarray) – The degrees of freedom of the Surface object, followed by the value of iota and G. e.g. [surface.x, iota, G] or [surface.x, iota] if optimize_G=False.

• derivatives (int, Optional) – 0 if no derivatives are requested, 1 if first derivatives are requested.

• optimize_G (bool, Optional) – True if G is a variable in the optimization problem, False otherwise.

Returns:

If derivatives=0, return res the residual of the optimization problem. If derivatives=1, return (res, dres) the residual and the Jacobian of the optimization problem.

boozer_penalty_constraints_vectorized(dofs, derivatives=0, constraint_weight=1.0, optimize_G=False, weight_inv_modB=True)

Replacement for the previous boozer_penalty_constraints function, which has issues on ubuntu. It is much faster and uses less memory since it calls a vectorized implementation in cpp. This is especially true when derivatives=2, i.e., when the Hessian is requested.

Parameters:

• dofs (ndarray) – The degrees of freedom of the Surface object, followed by the value of iota and G. e.g. [surface.x, iota, G] or [surface.x, iota] if optimize_G=False.

• derivatives (int, Optional) – 0 if no derivatives are requested, 1 if first derivatives are requested.

• constraint_weight (float, Optional) – The weight of the label constraint used when solving Boozer least squares.

• optimize_G (bool, Optional) – True if G is a variable in the optimization problem, False otherwise.

• weight_inv_modB (bool, Optional) – If True, weight the residual by modB so that it does not scale with coil currents. Defaults to True.

Returns:

(r, J, H) The residual vector, the Jacobian of the optimization problem, and the Hessian of the optimization problem. If derivatives=0, then J and H are None. If derivatives=1, then H is None.

Return type:

tuple

minimize_boozer_exact_constraints_newton(tol=1e-12, maxiter=10, iota=0.0, G=None, lm=[0.0, 0.0])

This function solves the constrained optimization problem

\[\text{min}_x ~ J(x)\]

subject to

\[\begin{split}l - l_0 &= 0 \\ z(\varphi=0,\theta=0) - 0 &= 0\end{split}\]

using the method of Lagrange multipliers and applying Newton’s method. In the above, \(J(x) = \frac{1}{2}\mathbf r(x)^T \mathbf r(x)\), and \(\mathbf r(x)\) contains the Boozer residuals at quadrature points \(1,\dots,n\).

The final constraint is not necessary for stellarator symmetric surfaces as it is automatically satisfied by the stellarator symmetric surface parametrization.

Parameters:

• tol (float, Optional) – The tolerance for the optimization. Defaults to 1e-12.

• maxiter (int, Optional) – The maximum number of iterations for the optimization. Defaults to 10.

• iota (float, Optional) – The initial guess for the value of the rotational transform on the surface. Defaults to 0.

• G (float, Optional) – The initial guess for the value of G on the surface. Defaults to None.

• lm (list, Optional) – The initial guesses for the Lagrange multipliers. Defaults to [0., 0.].

Returns:

A dictionary containing the results of the optimization. The dictionary contains the following keys in addition to others:

• ’residual’: the value of the residual at the solution

• ’jacobian’: the value of the jacobian at the solution

• ’iter’: the number of iterations taken to converge

• ’success’: True if the optimization converged, False otherwise

• ’G’: the value of G on the surface

• ’lm’: the value of the Lagrange multipliers

Return type:

dict

minimize_boozer_penalty_constraints_LBFGS(tol=0.001, maxiter=1000, constraint_weight=1.0, iota=0.0, G=None, limited_memory=True, weight_inv_modB=True, verbose=False)

This function uses L-BFGS to find the surface that approximately solves

\[\text{min}_x ~J(x) + \frac{1}{2} w_c (l - l_0)^2 + \frac{1}{2} w_c (z(\varphi=0, \theta=0) - 0)^2\]

where \(J(x) = \frac{1}{2}\mathbf r(x)^T \mathbf r(x)\), and \(\mathbf r(x)\) contains the Boozer residuals at quadrature points \(1,\dots,n\).

Parameters:

• tol (float, Optional) – The tolerance for the optimization. Defaults to 1e-3.

• maxiter (int, Optional) – The maximum number of iterations for the optimization. Defaults to 1000.

• constraint_weight (float, Optional) – The weight of the label constraint used when solving Boozer least squares.

• iota (float, Optional) – The initial guess for the value of the rotational transform on the surface. Defaults to 0.

• G (float, Optional) – The initial guess for the value of G on the surface. Defaults to None.

• limited_memory (bool, Optional) – If True, use the limited memory version of L-BFGS. Defaults to True.

• weight_inv_modB (bool, Optional) – If True, weight the residual by modB so that it does not scale with coil currents. Defaults to True.

• verbose (bool, Optional) – If True, print the optimization progress. Defaults to False.

Returns:

A dictionary containing the results of the optimization. The dictionary contains the following keys in addition to others:

• ’fun’: the value of the objective function at the solution

• ’gradient’: the gradient of the objective function at the solution

• ’iter’: the number of iterations taken to converge

• ’info’: the optimization result

• ’success’: True if the optimization converged, False otherwise

• ’G’: the value of G on the surface

• ’s’: the surface object

• ’iota’: the value of iota on the surface

• ’weight_inv_modB’: the value of weight_inv_modB used in the optimization

• ’type’: the type of optimization used

Return type:

res (dict)

minimize_boozer_penalty_constraints_ls(tol=1e-12, maxiter=10, constraint_weight=1.0, iota=0.0, G=None, method='lm', weight_inv_modB=True)

This function does the same as minimize_boozer_penalty_constraints_LBFGS, but instead of LBFGS it uses a nonlinear least squares algorithm when method='lm'. Options for the method are the same as for scipy.optimize.least_squares. If method='manual', then a damped Gauss-Newton method is used.

Parameters:

• tol (float, Optional) – The tolerance for the optimization. Defaults to 1e-12.

• maxiter (int, Optional) – The maximum number of iterations for the optimization. Defaults to 10.

• constraint_weight (float, Optional) – The weight of the label constraint used when solving Boozer least squares.

• iota (float, Optional) – The initial guess for the value of the rotational transform on the surface. Defaults to 0.

• G (float, Optional) – The initial guess for the value of G on the surface. Defaults to None.

• method (str, Optional) – The method to use for the optimization. Defaults to ‘lm’.

• weight_inv_modB (bool, Optional) – If True, weight the residual by modB so that it does not scale with coil currents. Defaults to True.

Returns:

A dictionary containing the results of the optimization. The dictionary contains the following keys in addition to others:

• ’residual’: the value of the residual at the solution

• ’gradient’: the value of the gradient at the solution

• ’jacobian’: the value of the jacobian at the solution

• ’success’: True if the optimization converged, False otherwise

• ’G’: the value of G on the surface

• ’s’: the surface object

• ’iota’: the value of iota on the surface

Return type:

res (dict)

minimize_boozer_penalty_constraints_newton(tol=1e-12, maxiter=10, constraint_weight=1.0, iota=0.0, G=None, stab=0.0, weight_inv_modB=True, verbose=False)

This function does the same as minimize_boozer_penalty_constraints_LBFGS, but instead of LBFGS it uses Newton’s method.

Parameters:

• tol (float, Optional) – The tolerance for the optimization. Defaults to 1e-12.

• maxiter (int, Optional) – The maximum number of iterations for the optimization. Defaults to 10.

• constraint_weight (float, Optional) – The weight of the label constraint used when solving Boozer least squares.

• iota (float, Optional) – The initial guess for the value of the rotational transform on the surface. Defaults to 0.

• G (float, Optional) – The initial guess for the value of G on the surface. Defaults to None.

• stab (float, Optional) – The stabilization parameter for the Newton method. Defaults to 0.

• weight_inv_modB (bool, Optional) – If True, weight the residual by modB so that it does not scale with coil currents. Defaults to True.

• verbose (bool, Optional) – If True, print the optimization progress. Defaults to False.

Returns:

A dictionary containing the results of the optimization. The dictionary contains the following keys in addition to others:

• ’residual’: the value of the residual at the solution

• ’jacobian’: the value of the Jacobian at the solution

• ’hessian’: the value of the Hessian at the solution

• ’iter’: the number of iterations taken to converge

• ’success’: True if the optimization converged, False otherwise

• ’G’: the value of G on the surface

• ’iota’: the value of iota on the surface

• ’PLU’: the LU decomposition of the hessian

• ’type’: ‘ls’.

• ’weight_inv_modB’: the value of weight_inv_modB used in the optimization

Return type:

dict

recompute_bell(parent=None)

Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.

This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.

Need to be implemented by classes that provide a dof_setter for external handling of DOFs.

run_code(iota, G=None)

Run the default solvers, i.e., run Newton’s method directly if you are computing a BoozerExact surface, and run BFGS followed by Newton if you are computing a BoozerLS surface.

Parameters:

• iota (float) – Guess for value of rotational transform on the surface.

• G (float, Optional) – Guess for value of G on surface, defaults to None. Note that if None is used, then the coil currents must be fixed.

Returns:

A dictionary containing the results of the optimization. The dictionary contains the following keys in addition to others:

• ”residual”: the residual of the optimization problem

• ”iter”: the number of iterations taken to converge

• ”success”: True if the optimization converged, False otherwise

• ”G”: the value of G on the surface

• ”s”: the surface object

• ”iota”: the value of iota on the surface

• ”PLU”: the LU decomposition of the hessian

Return type:

dict

solve_residual_equation_exactly_newton(tol=1e-10, maxiter=10, iota=0.0, G=None, verbose=False)

The function implements the BoozerExact approach by solving residual equation exactly using Newtons method.

For Newton’s method to be applied, we need the right balance of quadrature points, degrees of freedom and constraints. For this reason, this function is only implemented for surfaces of type SurfaceXYZTensorFourier right now.

Given ntor, mpol, nfp and stellsym, the surface is expected to be created in the following way:

phis = np.linspace(0, 1/nfp, 2*ntor+1, endpoint=False)
thetas = np.linspace(0, 1, 2*mpol+1, endpoint=False)
s = SurfaceXYZTensorFourier(
    mpol=mpol, ntor=ntor, stellsym=stellsym, nfp=nfp,
    quadpoints_phi=phis, quadpoints_theta=thetas)

Or the following two are also possible in the stellsym case:

phis = np.linspace(0, 1/nfp, 2*ntor+1, endpoint=False)
thetas = np.linspace(0, 0.5, mpol+1, endpoint=False)

or:

phis = np.linspace(0, 1/(2*nfp), ntor+1, endpoint=False)
thetas = np.linspace(0, 1, 2*mpol+1, endpoint=False)

and then:

s = SurfaceXYZTensorFourier(
    mpol=mpol, ntor=ntor, stellsym=stellsym, nfp=nfp,
    quadpoints_phi=phis, quadpoints_theta=thetas)

For the stellarator symmetric case, there is some redundancy between DOFs. This is taken care of inside this function.

In the non-stellarator-symmetric case, the surface has (2*ntor+1)*(2*mpol+1) many quadrature points and 3*(2*ntor+1)*(2*mpol+1) many dofs.

Equations:

• Boozer residual in x, y, and z at all quadrature points

• z(0, 0) = 0

• label constraint (e.g. volume or flux)

Unknowns:

• Surface dofs

• iota

• G

So we end up having 3*(2*ntor+1)*(2*mpol+1) + 2 equations and the same number of unknowns.

In the stellarator-symmetric case, we have D = (ntor+1)*(mpol+1)+ ntor*mpol + 2*(ntor+1)*mpol + 2*ntor*(mpol+1) = 6*ntor*mpol + 3*ntor + 3*mpol + 1 many dofs in the surface. After calling surface.get_stellsym_mask() we have kicked out 2*ntor*mpol + ntor + mpol quadrature points, i.e. we have 2*ntor*mpol + ntor + mpol + 1 quadrature points remaining. In addition we know that the x coordinate of the residual at phi=0=theta is also always satisfied. In total this leaves us with 3*(2*ntor*mpol + ntor + mpol) + 2 equations for the boozer residual, plus 1 equation for the label, which is the same as the number of surface dofs + 2 extra unknowns given by iota and G.

Parameters:

• tol (float, Optional) – The tolerance for the optimization. Defaults to 1e-10.

• maxiter (int, Optional) – The maximum number of iterations for the optimization. Defaults to 10.

• iota (float, Optional) – The initial guess for the value of the rotational transform on the surface. Defaults to 0.

• G (float, Optional) – The initial guess for the value of G on the surface. Defaults to None.

• verbose (bool, Optional) – If True, print the optimization progress. Defaults to False.

Returns:

A dictionary containing the results of the optimization. The dictionary contains the following keys in addition to others:

• ’residual’: the value of the residual at the solution

• ’jacobian’: the value of the jacobian at the solution

• ’iter’: the number of iterations taken to converge

• ’success’: True if the optimization converged, False otherwise

• ’G’: the value of G on the surface

• ’s’: the surface object

• ’iota’: the value of iota on the surface

• ’PLU’: the LU decomposition of the jacobian

• ’mask’: a mask for the residuals that are not used in the optimization

• ’type’: ‘exact’.

• ’vjp’: the vector-Jacobian product for the optimization

Return type:

dict

class simsopt.geo.CircularPort(ox=1.0, oy=0.0, oz=0.0, ax=1.0, ay=0.0, az=0.0, ir=0.5, thick=0.01, l0=0.0, l1=0.5)

Bases: Port

Class representing a port with cylindrical geometry; specifically, with a circular cross section.

__init__(ox=1.0, oy=0.0, oz=0.0, ax=1.0, ay=0.0, az=0.0, ir=0.5, thick=0.01, l0=0.0, l1=0.5)

Initializes the CircularPort class according to the geometric port parameters.

Parameters:

• ox (floats) – Cartesian x, y, and z coordinates of the origin point of the port in 3d space

• oy (floats) – Cartesian x, y, and z coordinates of the origin point of the port in 3d space

• oz (floats) – Cartesian x, y, and z coordinates of the origin point of the port in 3d space

• ax (floats) – Cartesian x, y, and z components of a vector in the direction of the port’s axis

• ay (floats) – Cartesian x, y, and z components of a vector in the direction of the port’s axis

• az (floats) – Cartesian x, y, and z components of a vector in the direction of the port’s axis

• ir (float) – Inner radius of the port

• thick (float) – Wall thickness of the port

• l0 (double) – Locations of the ends of the port, specified by the signed distance along the port axis from the origin point

• l1 (double) – Locations of the ends of the port, specified by the signed distance along the port axis from the origin point

collides(x, y, z, gap=0.0)

Determines if the user input point(s) collide with the port

Parameters:

• x (float arrays) – Cartesian x, y, and z coordinates of the point(s) to be assessed. Dimensions must be the same.

• y (float arrays) – Cartesian x, y, and z coordinates of the point(s) to be assessed. Dimensions must be the same.

• z (float arrays) – Cartesian x, y, and z coordinates of the point(s) to be assessed. Dimensions must be the same.

• gap (float (optional)) – Minimum gap spacing required between the point and the exterior of the port. Default is 0.

Returns:

colliding – Array of logical values indicating whether each of the input points collide with the port.

Return type:

logical array

mesh_representation(n_edges=100)

Constructs a triangular mesh representation of the port for plotting and visualization.

Parameters:

n_edges (integer (optional)) – Number of edges for the polygon used to approximate the circular cross-section. Default is 100.

Returns:

• x, y, z (1D arrays) – Cartesian x, y, and z coordinates of the vertices of the mesh.

• triangles (2D array) – Indices of the vertices (as listed in the x, y, and z arrays) surrounding each triangular face within the mesh. Each row (dimension 1) represents a face.

plot(n_edges=100, **kwargs)

Places a representation of the port on a 3D plot. Currently only works with the mayavi plotting package.

Parameters:

• n_edges (integer (optional)) – Number of edges for the polygon used to approximate the circular cross-section. Default is 100.

• **kwargs – Keyword arguments to be passed to the mayavi surface module for the port

Returns:

surf – Reference to the plotted port

Return type:

mayavi.modules.surface.Surface class instance

repeat_via_symmetries(nfp, stell_sym, include_self=True)

Generates a set of equivalent ports to uphold toroidal and/or stellarator symmetry.

Parameters:

• nfp (integer) – Number of toroidal field periods

• stell_sym (logical) – If true, stellarator symmetry will be assumed and equivalent ports will be created for every half-period

• include_self (logical (optional)) – If true, the port instance calling this method will be included in the returned set of ports (see below). Default is True.

Returns:

ports – A set of all the symmetric ports, including the one represented by the calling Port class instance if include_self==True.

Return type:

PortSet class instance

to_vtk(filename, n_edges=100)

Export a mesh representation of the port to a VTK file, which can be read with ParaView.

Note: This function requires the pyevtk python package, which can be installed using pip install pyevtk.

Parameters:

• filename (string) – Name of the VTK file, without extension, to create.

• n_edges (integer (optional)) – Number of edges for the polygon used to approximate the circular cross-section. Default is 100.

class simsopt.geo.CoilStrain(framedcurve, width=0.001)

Bases: Optimizable

This class evaluates the torsional and binormal curvature strains on HTS, based on a filamentary model of the coil and the orientation of the HTS tape.

As defined in,

Paz Soldan, “Non-planar coil winding angle optimization for compatibility with non-insulated high-temperature superconducting magnets”, Journal of Plasma Physics 86 (2020), doi:10.1017/S0022377820001208,

the expressions for the strains are:

\[ \begin{align}\begin{aligned}\epsilon_{\text{tor}} = \frac{\tau^2 w^2}{12}\\\epsilon_{\text{bend}} = \frac{w |\hat{\textbf{b}} \cdot \boldsymbol{\kappa}|}{2},\end{aligned}\end{align} \]

where \(\tau\) is the torsion of the tape frame, \(\hat{\textbf{b}}\) is the frame binormal vector, \(\boldsymbol{\kappa}\) is the curvature vector of the filamentary coil, and \(w\) is the width of the tape.

This class is not intended to be used as an objective function inside optimization. For that purpose you should instead use LPBinormalCurvatureStrainPenalty or LPTorsionalStrainPenalty. Those classes also compute gradients whereas this class does not.

binormal_curvature_strain()

Returns the value of the torsional strain, \(\epsilon_{\text{bend}}\), along the quadpoints defining the filamentary coil.

torsional_strain()

Returns the value of the torsional strain, \(\epsilon_{\text{tor}}\), along the quadpoints defining the filamentary coil.

class simsopt.geo.Curve(**kwargs)

Bases: Optimizable

Curve is a base class for various representations of curves in SIMSOPT using the graph based Optimizable framework with external handling of DOFS as well.

centroid()

Compute the centroid of the curve

\[\mathbf{c} = \frac{1}{L} \int_0^L \mathbf{\gamma}(l) dl\]

where \(\gamma\) is the position on the curve. Note that this function was once called center but this conflicts with the center property of the C++ CurvePlanarFourier implementation.

dfrenet_frame_by_dcoeff()

This function returns the derivative of the curve’s Frenet frame,

\[\left(\frac{\partial \mathbf{t}}{\partial \mathbf{c}}, \frac{\partial \mathbf{n}}{\partial \mathbf{c}}, \frac{\partial \mathbf{b}}{\partial \mathbf{c}}\right),\]

with respect to the curve dofs, where \((\mathbf t, \mathbf n, \mathbf b)\) correspond to the Frenet frame, and \(\mathbf c\) are the curve dofs.

dgamma_by_dcoeff_vjp(v)

dgammadash_by_dcoeff_vjp(v)

dgammadashdash_by_dcoeff_vjp(v)

dgammadashdashdash_by_dcoeff_vjp(v)

dincremental_arclength_by_dcoeff_vjp(v)

This function returns the vector Jacobian product

\[v^T \frac{\partial \|\Gamma'\|}{\partial \mathbf{c}}\]

where \(\|\Gamma'\|\) is the incremental arclength, \(\Gamma'\) is the tangent to the curve and \(\mathbf{c}\) are the curve dofs.

dkappa_by_dcoeff_impl(dkappa_by_dcoeff)

This function computes the derivative of the curvature with respect to the curve coefficients.

\[\frac{\partial \kappa}{\partial \mathbf c}\]

where \(\mathbf c\) are the curve dofs, \(\kappa\) is the curvature.

dkappa_by_dcoeff_vjp(v)

This function returns the vector Jacobian product

\[v^T \frac{\partial \kappa}{\partial \mathbf{c}}\]

where \(\mathbf c\) are the curve dofs and \(\kappa\) is the curvature.

dkappadash_by_dcoeff()

This function returns

\[\frac{\partial \kappa'(\phi)}{\partial \mathbf{c}}.\]

where \(\mathbf{c}\) are the curve dofs, and \(\kappa\) is the curvature.

dtorsion_by_dcoeff_impl(dtorsion_by_dcoeff)

This function returns the derivative of torsion with respect to the curve dofs.

\[\frac{\partial \tau}{\partial \mathbf c}\]

where \(\mathbf c\) are the curve dofs, and \(\tau\) is the torsion.

dtorsion_by_dcoeff_vjp(v)

This function returns the vector Jacobian product

\[v^T \frac{\partial \tau}{\partial \mathbf{c}}\]

where \(\mathbf c\) are the curve dofs, and \(\tau\) is the torsion.

frenet_frame()

This function returns the Frenet frame, \((\mathbf{t}, \mathbf{n}, \mathbf{b})\), associated to the curve.

kappa_impl(kappa)

This function implements the curvature, \(\kappa(\varphi)\).

kappadash()

This function returns \(\kappa'(\phi)\), where \(\kappa\) is the curvature.

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

Plot the curve in 3D using matplotlib.pyplot, mayavi, or plotly.

Parameters:

• engine – The graphics engine to use. Available settings are "matplotlib", "mayavi", and "plotly".

• ax – The axis object on which to plot. 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.

• plot_derivative – Whether to plot the tangent of the curve too. Not implemented for plotly.

• close – Whether to connect the first and last point on the curve. Can lead to surprising results when only quadrature points on a part of the curve are considered, e.g. when exploting rotational symmetry.

• axis_equal – For matplotlib, whether all three dimensions should be scaled equally.

• kwargs – Any additional arguments to pass to the plotting function, like color='r'.

Returns:

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

recompute_bell(parent=None)

For derivative classes of Curve, all of which also subclass from C++ Curve class, call invalidate_cache which is implemented in C++ side.

torsion_impl(torsion)

This function returns the torsion, \(\tau\), of a curve.

class simsopt.geo.CurveCurveDistance(curves, minimum_distance, num_basecurves=None, downsample=1)

Bases: Optimizable

CurveCurveDistance is a class that computes

\[J = \sum_{i = 1}^{\text{num_coils}} \sum_{j = 1}^{i-1} d_{i,j}\]

where

\[\begin{split}d_{i,j} = \int_{\text{curve}_i} \int_{\text{curve}_j} \max(0, d_{\min} - \| \mathbf{r}_i - \mathbf{r}_j \|_2)^2 ~dl_j ~dl_i\\\end{split}\]

and \(\mathbf{r}_i\), \(\mathbf{r}_j\) are points on coils \(i\) and \(j\), respectively. \(d_\min\) is a desired threshold minimum intercoil distance. This penalty term is zero when the points on coil \(i\) and coil \(j\) lie more than \(d_\min\) away from one another, for \(i, j \in \{1, \cdots, \text{num_coils}\}\)

If num_basecurves is passed, then the code only computes the distance to the first num_basecurves many curves, which is useful when the coils satisfy symmetries that can be exploited.

J()

This returns the value of the quantity.

compute_candidates()

dJ(*args, partials=False, **kwargs)

recompute_bell(parent=None)

Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.

This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.

Need to be implemented by classes that provide a dof_setter for external handling of DOFs.

return_fn_map: Dict[str, Callable] = {'J': <function CurveCurveDistance.J>, 'dJ': <function derivative_dec.<locals>._derivative_dec>}

shortest_distance()

shortest_distance_among_candidates()

class simsopt.geo.CurveFilament(framedcurve, dn, db)

Bases: FramedCurve

__init__(framedcurve, dn, db)

Given a FramedCurve, defining a normal and binormal vector, create a grid of curves by shifting along the normal and binormal vector.

The idea is explained well in Figure 1 in the reference:

Singh et al, “Optimization of finite-build stellarator coils”, Journal of Plasma Physics 86 (2020), doi:10.1017/S0022377820000756.

Parameters:

• curve – the underlying curve

• dn – how far to move in normal direction

• db – how far to move in binormal direction

• rotation – angle along the curve to rotate the frame.

dgamma_by_dcoeff_vjp(v)

dgammadash_by_dcoeff_vjp(v)

gamma_impl(self: simsoptpp.Curve, arg0: numpy.ndarray[numpy.float64], arg1: numpy.ndarray[numpy.float64]) → None

gammadash_impl(gammadash)

recompute_bell(parent=None)

For derivative classes of Curve, all of which also subclass from C++ Curve class, call invalidate_cache which is implemented in C++ side.

class simsopt.geo.CurveHelical(quadpoints, order, m=5, ell=2, R0=1.0, r=0.3, **kwargs)

Bases: JaxCurve

A helical curve lying on a circular-cross-section axisymmetric torus.

This curve type can describe the coils from the Hanson & Cary (1984) and Cary & Hanson (1986) papers on optimization for reduced stochasticity, as well as the Compact Toroidal Hybrid (CTH) device at Auburn University.

The helical curve lies on a circular-cross-section axisymmetric torus with major radius \(R_0\) and minor radius \(r\). Following Hanson and Cary, the position along the curve is written in terms of a poloidal angle on the winding surface \(\eta\) and the standard toroidal angle \(\phi\):

\[\begin{split}\begin{align*} x &= [R_0 + r \cos(\eta)] \cos(\phi), \\ y &= [R_0 + r \cos(\eta)] \sin(\phi), \\ z &= -r \sin(\eta). \end{align*}\end{split}\]

Along the curve, the poloidal angle depends on the toroidal angle both through a secular linear term and periodic Fourier terms:

\[\eta = \frac{m \phi}{l} + A_0 + \sum_{k=1}^N \left[ A_k \cos(k \phi m / l) + B_k \sin(k \phi m / l) \right]\]

When the toroidal angle \(\phi\) increases by \(2 \pi l\), the poloidal angle \(\eta\) increases by \(2 \pi m\), so the “rotational transform” of the curve is \(m / l\). The shape of the helical curve on the surface can be adjusted through the \(A_k\) and \(B_k\) Fourier coefficients, which are the optimizable degrees of freedom for this class. Although \(\phi\) in the above formulas is periodic in the range \([0, 2 \pi)\), the quadpoints argument should be in the range \([0, 1)\) as usual for simsopt curves.

The 1986 Cary-Hanson paper has a more general parameterization than the 1984 Hanson-Cary paper, relaxing the constraint that the curves lie on a circle in the R-z plane, although no results in the paper use this freedom. We do not consider this more general parameterization here in this class. For a more general parameterization of helical curves and coils in simsopt, see simsopt.geo.CurveXYZFourierSymmetries.

The order in which the degrees of freedom are stored in the state vector x is

\[A_0, A_1, \ldots, A_N, B_1, \ldots, B_N.\]

The default values of m, ell, R0, and r correspond to the coils in the Hanson-Cary and Cary-Hanson papers.

Parameters:

• quadpoints – (int or array-like) Grid points (or number thereof) along the curve.

• order – Maximum (inclusive) Fourier mode number \(N\) for \(A_k\) and \(B_k\).

• m – Integer \(m\) describing helicity of the coil.

• ell – Integer \(l\) describing helicity of the coil.

• R0 – Major radius of the toroidal surface on which the coil lies.

• r – Minor radius of the toroidal surface on which the coil lies.

get_dofs(self: simsoptpp.Curve) → list[float]

num_dofs(self: simsoptpp.Curve) → int

set_dofs_impl(self: simsoptpp.Curve, arg0: list[float]) → None

class simsopt.geo.CurveLength(curve)

Bases: Optimizable

CurveLength is a class that computes the length of a curve, i.e.

\[J = \int_{\text{curve}}~dl.\]

J()

This returns the value of the quantity.

dJ(*args, partials=False, **kwargs)

return_fn_map: Dict[str, Callable] = {'J': <function CurveLength.J>, 'dJ': <function derivative_dec.<locals>._derivative_dec>}

class simsopt.geo.CurvePerturbed(curve, sample)

Bases: Curve, Curve

A perturbed curve.

__init__(curve, sample)

Perturb a underlying simsopt.geo.curve.Curve object by drawing a perturbation from a GaussianSampler.

Comment: Doing anything involving randomness in a reproducible way requires care. Even more so, when doing things in parallel. Let’s say we have a list of simsopt.geo.curve.Curve objects curves that represent a stellarator, and now we want to consider N perturbed stellarators. Let’s also say we have multiple MPI ranks. To avoid the same thing happening on the different MPI ranks, we could pick a different seed on each rank. However, then we get different results depending on the number of MPI ranks that we run on. Not ideal. Instead, we should pick a new seed for each \(1\le i\le N\). e.g.

from np.random import SeedSequence, PCG64DXSM, Generator
import numpy as np
curves = ...
sigma = 0.01
length_scale = 0.2
sampler = GaussianSampler(curves[0].quadpoints, sigma, length_scale, n_derivs=1)
globalseed = 1
N = 10 # number of perturbed stellarators
seeds = SeedSequence(globalseed).spawn(N)
idx_start, idx_end = split_range_between_mpi_rank(N) # e.g. [0, 5) on rank 0, [5, 10) on rank 1
perturbed_curves = [] # this will be a List[List[Curve]], with perturbed_curves[i] containing the perturbed curves for the i-th stellarator
for i in range(idx_start, idx_end):
    rg = Generator(PCG64DXSM(seeds_sys[j]))
    stell = []
    for c in curves:
        pert = PerturbationSample(sampler_systematic, randomgen=rg)
        stell.append(CurvePerturbed(c, pert))
    perturbed_curves.append(stell)

dgamma_by_dcoeff_vjp(v)

dgammadash_by_dcoeff_vjp(v)

dgammadashdash_by_dcoeff_vjp(v)

dgammadashdashdash_by_dcoeff_vjp(v)

gamma_impl(self: simsoptpp.Curve, arg0: numpy.ndarray[numpy.float64], arg1: numpy.ndarray[numpy.float64]) → None

gammadash_impl(gammadash)

gammadashdash_impl(gammadashdash)

gammadashdashdash_impl(gammadashdashdash)

recompute_bell(parent=None)

For derivative classes of Curve, all of which also subclass from C++ Curve class, call invalidate_cache which is implemented in C++ side.

resample()

class simsopt.geo.CurvePlanarFourier(quadpoints, order, dofs=None)

Bases: CurvePlanarFourier, Curve

CurvePlanarFourier is a curve that is restricted to lie in a plane. The shape of the curve within the plane is represented by a Fourier series in polar coordinates centered at the center of curve. The resulting planar curve is then rotated in three dimensions using a quaternion, and finally a translation is applied by the center point (X, Y, Z). The Fourier series in polar coordinates is

\[r(\phi) = \sum_{m=0}^{\text{order}} r_{c,m}\cos(m \phi) + \sum_{m=1}^{\text{order}} r_{s,m}\sin(m \phi).\]

The rotation quaternion is

\[ \begin{align}\begin{aligned}\bf{q} &= [q0,qi,qj,qk]\\&= [\cos(\theta / 2), \hat{x}\sin(\theta / 2), \hat{y}\sin(\theta / 2), \hat{z}\sin(\theta / 2)]\end{aligned}\end{align} \]

where \(\theta\) is the counterclockwise rotation angle about a unit axis \((\hat{x},\hat{y},\hat{z})\). Details of the quaternion rotation can be found for example in pages 575-576 of https://www.cis.upenn.edu/~cis5150/ws-book-Ib.pdf.

A quaternion is used for rotation rather than other methods for rotation to prevent gimbal locking during optimization. The quaternion is normalized before being applied to prevent scaling of the curve. The dofs themselves are not normalized. This results in a redundancy in the optimization, where several different sets of dofs may correspond to the same normalized quaternion. Normalizing the dofs directly would create a dependence between the quaternion dofs, which may cause issues during optimization.

The dofs are stored in the order

\[[r_{c,0}, \cdots, r_{c,\text{order}}, r_{s,1}, \cdots, r_{s,\text{order}}, q0, qi, qj, qk, X, Y, Z]\]

Parameters:

• quadpoints (array) – Array of quadrature points.

• order (int) – Order of the Fourier series.

• dofs (array, optional) – Array of dofs.

_make_names(order)

This function returns the names of the dofs associated to this object.

Parameters:

order (int) – Order of the Fourier series.

Returns:

List of dof names.

get_dofs()

This function returns the dofs associated to this object.

set_dofs(dofs)

This function sets the dofs associated to this object.

class simsopt.geo.CurveRZFourier(quadpoints, order, nfp, stellsym, dofs=None)

Bases: CurveRZFourier, Curve

CurveRZFourier is a curve that is represented in cylindrical coordinates using the following Fourier series:

\[\begin{split}r(\phi) &= \sum_{m=0}^{\text{order}} r_{c,m}\cos(n_{\text{fp}} m \phi) + \sum_{m=1}^{\text{order}} r_{s,m}\sin(n_{\text{fp}} m \phi) \\ z(\phi) &= \sum_{m=0}^{\text{order}} z_{c,m}\cos(n_{\text{fp}} m \phi) + \sum_{m=1}^{\text{order}} z_{s,m}\sin(n_{\text{fp}} m \phi)\end{split}\]

If stellsym = True, then the \(\sin\) terms for \(r\) and the \(\cos\) terms for \(z\) are zero.

For the stellsym = False case, the dofs are stored in the order

\[[r_{c,0}, \cdots, r_{c,\text{order}}, r_{s,1}, \cdots, r_{s,\text{order}}, z_{c,0},....]\]

or in the stellsym = True case they are stored

\[[r_{c,0},...,r_{c,order},z_{s,1},...,z_{s,order}]\]

get_dofs()

This function returns the dofs associated to this object.

set_dofs(dofs)

This function sets the dofs associated to this object.

class simsopt.geo.CurveSurfaceDistance(curves, surface, minimum_distance)

Bases: Optimizable

CurveSurfaceDistance is a class that computes

\[J = \sum_{i = 1}^{\text{num_coils}} d_{i}\]

where

\[\begin{split}d_{i} = \int_{\text{curve}_i} \int_{surface} \max(0, d_{\min} - \| \mathbf{r}_i - \mathbf{s} \|_2)^2 ~dl_i ~ds\\\end{split}\]

and \(\mathbf{r}_i\), \(\mathbf{s}\) are points on coil \(i\) and the surface, respectively. \(d_\min\) is a desired threshold minimum coil-to-surface distance. This penalty term is zero when the points on all coils \(i\) and on the surface lie more than \(d_\min\) away from one another.

J()

This returns the value of the quantity.

compute_candidates()

dJ(*args, partials=False, **kwargs)

recompute_bell(parent=None)

Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.

This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.

Need to be implemented by classes that provide a dof_setter for external handling of DOFs.

return_fn_map: Dict[str, Callable] = {'J': <function CurveSurfaceDistance.J>, 'dJ': <function derivative_dec.<locals>._derivative_dec>}

shortest_distance()

shortest_distance_among_candidates()

class simsopt.geo.CurveXYZFourier(quadpoints, order, dofs=None)

Bases: CurveXYZFourier, Curve

CurveXYZFourier is a curve that is represented in Cartesian coordinates using the following Fourier series:

\[\begin{split}x(\theta) &= \sum_{m=0}^{\text{order}} x_{c,m}\cos(m\theta) + \sum_{m=1}^{\text{order}} x_{s,m}\sin(m\theta) \\ y(\theta) &= \sum_{m=0}^{\text{order}} y_{c,m}\cos(m\theta) + \sum_{m=1}^{\text{order}} y_{s,m}\sin(m\theta) \\ z(\theta) &= \sum_{m=0}^{\text{order}} z_{c,m}\cos(m\theta) + \sum_{m=1}^{\text{order}} z_{s,m}\sin(m\theta)\end{split}\]

The dofs are stored in the order

\[[x_{c,0}, x_{s,1}, x_{c,1},\cdots x_{s,\text{order}}, x_{c,\text{order}},y_{c,0},y_{s,1},y_{c,1},\cdots]\]

_make_names(order)

This function returns the names of the dofs associated to this object. :param order: Order of the Fourier series. :type order: int

Returns:

List of dof names.

get_dofs()

This function returns the dofs associated to this object.

static load_curves_from_file(filename, order=None, ppp=20, delimiter=',')

This function loads a file containing Fourier coefficients for several coils. The file is expected to have 6*num_coils many columns, and order+1 many rows. The columns are in the following order,

sin_x_coil1, cos_x_coil1, sin_y_coil1, cos_y_coil1, sin_z_coil1, cos_z_coil1, sin_x_coil2, cos_x_coil2, sin_y_coil2, cos_y_coil2, sin_z_coil2, cos_z_coil2, …

static load_curves_from_makegrid_file(filename: str, order: int, ppp=20, group_names=None)

This function loads a Makegrid input file containing the Cartesian coordinates for several coils and finds the corresponding Fourier coefficients through an fft. The format is described at https://princetonuniversity.github.io/STELLOPT/MAKEGRID

Parameters:

• filename – file to load.

• order – maximum mode number in the Fourier series.

• ppp – points-per-period: number of quadrature points per period.

• group_names – List of coil group names (str). If not ‘None’, only get coils in coil groups that are in the list.

Returns:

A list of CurveXYZFourier objects.

set_dofs(dofs)

This function sets the dofs associated to this object.

class simsopt.geo.CurveXYZFourierSymmetries(quadpoints, order, nfp, stellsym, ntor=1, **kwargs)

Bases: JaxCurve

A curve representation that allows for stellarator and discrete rotational symmetries. This class can be used to represent a helical coil that does not lie on a torus. The coordinates of the curve are given by:

\[\begin{split}x(\theta) &= \hat x(\theta) \cos(2 \pi \theta n_{\text{tor}}) - \hat y(\theta) \sin(2 \pi \theta n_{\text{tor}})\\ y(\theta) &= \hat x(\theta) \sin(2 \pi \theta n_{\text{tor}}) + \hat y(\theta) \cos(2 \pi \theta n_{\text{tor}})\\ z(\theta) &= \sum_{m=1}^{\text{order}} z_{s,m} \sin(2 \pi n_{\text{fp}} m \theta)\end{split}\]

where

\[\begin{split}\hat x(\theta) &= x_{c, 0} + \sum_{m=1}^{\text{order}} x_{c,m} \cos(2 \pi n_{\text{fp}} m \theta)\\ \hat y(\theta) &= \sum_{m=1}^{\text{order}} y_{s,m} \sin(2 \pi n_{\text{fp}} m \theta)\\\end{split}\]

if the coil is stellarator symmetric. When the coil is not stellarator symmetric, the formulas above become

\[\begin{split}x(\theta) &= \hat x(\theta) \cos(2 \pi \theta n_{\text{tor}}) - \hat y(\theta) \sin(2 \pi \theta n_{\text{tor}})\\ y(\theta) &= \hat x(\theta) \sin(2 \pi \theta n_{\text{tor}}) + \hat y(\theta) \cos(2 \pi \theta n_{\text{tor}})\\ z(\theta) &= z_{c, 0} + \sum_{m=1}^{\text{order}} \left[ z_{c, m} \cos(2 \pi n_{\text{fp}} m \theta) + z_{s, m} \sin(2 \pi n_{\text{fp}} m \theta) \right]\end{split}\]

where

\[\begin{split}\hat x(\theta) &= x_{c, 0} + \sum_{m=1}^{\text{order}} \left[ x_{c, m} \cos(2 \pi n_{\text{fp}} m \theta) + x_{s, m} \sin(2 \pi n_{\text{fp}} m \theta) \right] \\ \hat y(\theta) &= y_{c, 0} + \sum_{m=1}^{\text{order}} \left[ y_{c, m} \cos(2 \pi n_{\text{fp}} m \theta) + y_{s, m} \sin(2 \pi n_{\text{fp}} m \theta) \right] \\\end{split}\]

Parameters:

• quadpoints – number of grid points/resolution along the curve,

• order – how many Fourier harmonics to include in the Fourier representation,

• nfp – discrete rotational symmetry number,

• stellsym – stellaratory symmetry if True, not stellarator symmetric otherwise,

• ntor – the number of times the curve wraps toroidally before biting its tail. Note, it is assumed that nfp and ntor are coprime. If they are not coprime, then then the curve actually has nfp_new:=nfp // gcd(nfp, ntor), and ntor_new:=ntor // gcd(nfp, ntor). The operator // is integer division. To avoid confusion, we assert that ntor and nfp are coprime at instantiation.

get_dofs(self: simsoptpp.Curve) → list[float]

num_dofs(self: simsoptpp.Curve) → int

set_dofs_impl(self: simsoptpp.Curve, arg0: list[float]) → None

class simsopt.geo.FrameRotation(quadpoints, order, scale=1.0, dofs=None)

Bases: Optimizable

__init__(quadpoints, order, scale=1.0, dofs=None)

Defines the rotation angle with respect to a reference orthonormal frame (either frenet or centroid). For example, can be used to define the rotation of a multifilament pack; alpha in Figure 1 of Singh et al, “Optimization of finite-build stellarator coils”, Journal of Plasma Physics 86 (2020), doi:10.1017/S0022377820000756

alpha(quadpoints)

alphadash(quadpoints)

dalpha_by_dcoeff_vjp(quadpoints, v)

dalphadash_by_dcoeff_vjp(quadpoints, v)

class simsopt.geo.FramedCurve(curve, rotation=None)

Bases: Curve, Curve

__init__(curve, rotation=None)

A FramedCurve defines an orthonormal basis around a Curve, where one basis is taken to be the tangent along the Curve. The frame is defined with respect to a reference frame, either centroid or frenet. A rotation angle defines the rotation with respect to this reference frame.

class simsopt.geo.FramedCurveCentroid(curve, rotation=None)

Bases: FramedCurve

Implementation of the centroid frame introduced in Singh et al, “Optimization of finite-build stellarator coils”, Journal of Plasma Physics 86 (2020), doi:10.1017/S0022377820000756. Given a curve, one defines a reference frame using the normal and binormal vector based on the centoid of the coil. In addition, we specify an angle along the curve that defines the rotation with respect to this reference frame.

The idea is explained well in Figure 1 in the reference above.

dframe_binormal_curvature_by_dcoeff_vjp(v)

dframe_torsion_by_dcoeff_vjp(v)

frame_binormal_curvature()

frame_torsion()

Exports frame torsion along a curve

rotated_frame()

rotated_frame_dash()

rotated_frame_dash_dcoeff_vjp(v, dn, db, arg=0)

rotated_frame_dcoeff_vjp(v, dn, db, arg=0)

class simsopt.geo.FramedCurveFrenet(curve, rotation=None)

Bases: FramedCurve

Given a curve, one defines a reference frame using the Frenet normal and binormal vectors:

tangent = dr/dl

normal = (dtangent/dl)/||dtangent/dl||

binormal = tangent x normal

In addition, we specify an angle along the curve that defines the rotation with respect to this reference frame.

dframe_binormal_curvature_by_dcoeff_vjp(v)

dframe_torsion_by_dcoeff_vjp(v)

frame_binormal_curvature()

frame_torsion()

Exports frame torsion along a curve

rotated_frame()

rotated_frame_dash()

rotated_frame_dash_dcoeff_vjp(v, dn, db, arg=0)

rotated_frame_dcoeff_vjp(v, dn, db, arg=0)

class simsopt.geo.GaussianSampler(points: Sequence[Real] | ndarray[tuple[Any, ...], dtype[float64]], sigma: float, length_scale: float, n_derivs: int = 1)

Bases: GSONable

Generate a periodic gaussian process on the interval [0, 1] on a given list of quadrature points. The process has standard deviation sigma a correlation length scale length_scale. Large values of length_scale correspond to smooth processes, small values result in highly oscillatory functions. Also has the ability to sample the derivatives of the function.

We consider the kernel

\[\kappa(d) = \sigma^2 \exp(-d^2/l^2)\]

and then consider a Gaussian process with covariance

\[Cov(X(s), X(t)) = \sum_{i=-\infty}^\infty \sigma^2 \exp(-(s-t+i)^2/l^2)\]

the sum is used to make the kernel periodic and in practice the infinite sum is truncated.

Parameters:

• points – the quadrature points along which the perturbation should be computed.

• sigma – standard deviation of the underlying gaussian process (measure for the magnitude of the perturbation).

• length_scale – length scale of the underlying gaussian process (measure for the smoothness of the perturbation).

• n_derivs – number of derivatives of the gaussian process to sample.

draw_sample(randomgen=None)

Returns a list of n_derivs+1 arrays of size (len(points), 3), containing the perturbation and the derivatives.

length_scale: float

n_derivs: int = 1

points: Sequence[Real] | ndarray[tuple[Any, ...], dtype[float64]]

sigma: float

class simsopt.geo.Iotas(boozer_surface)

Bases: Optimizable

This term returns the rotational transform on a boozer surface as well as its derivative with respect to the coil degrees of freedom.

J()

compute()

dJ(*args, partials=False, **kwargs)

recompute_bell(parent=None)

Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.

This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.

Need to be implemented by classes that provide a dof_setter for external handling of DOFs.

class simsopt.geo.JaxCurve(quadpoints, gamma_pure, **kwargs)

Bases: Curve, Curve

A class for curves defined by a pure function.

Parameters:

• quadpoints (array) – Array of quadrature points.

• gamma_pure (function) – Pure function for the curve.

• **kwargs – Additional keyword arguments.

dgamma_by_dcoeff_impl(dgamma_by_dcoeff)

This function returns

\[\frac{\partial \Gamma}{\partial \mathbf c}\]

where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.

dgamma_by_dcoeff_vjp_impl(v)

This function returns the vector Jacobian product

\[v^T \frac{\partial \Gamma}{\partial \mathbf c}\]

where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.

dgammadash_by_dcoeff_impl(dgammadash_by_dcoeff)

This function returns

\[\frac{\partial \Gamma'}{\partial \mathbf c}\]

where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.

dgammadash_by_dcoeff_vjp_impl(v)

This function returns

\[\mathbf v^T \frac{\partial \Gamma'}{\partial \mathbf c}\]

where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.

dgammadashdash_by_dcoeff_impl(dgammadashdash_by_dcoeff)

This function returns

\[\frac{\partial \Gamma''}{\partial \mathbf c}\]

where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.

dgammadashdash_by_dcoeff_vjp_impl(v)

This function returns the vector Jacobian product

\[v^T \frac{\partial \Gamma''}{\partial \mathbf c}\]

where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.

dgammadashdashdash_by_dcoeff_impl(dgammadashdashdash_by_dcoeff)

This function returns

\[\frac{\partial \Gamma'''}{\partial \mathbf c}\]

where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.

dgammadashdashdash_by_dcoeff_vjp_impl(v)

This function returns the vector Jacobian product

\[v^T \frac{\partial \Gamma'''}{\partial \mathbf c}\]

where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.

dkappa_by_dcoeff_vjp(v)

This function returns the vector Jacobian product

\[v^T \frac{\partial \kappa}{\partial \mathbf{c}}\]

where \(\mathbf{c}\) are the curve dofs and \(\kappa\) is the curvature.

dtorsion_by_dcoeff_vjp(v)

This function returns the vector Jacobian product

\[v^T \frac{\partial \tau}{\partial \mathbf{c}}\]

where \(\mathbf{c}\) are the curve dofs, and \(\tau\) is the torsion.

gamma_impl(gamma, quadpoints)

This function returns the x,y,z coordinates of the curve \(\Gamma\).

gammadash_impl(gammadash)

This function returns \(\Gamma'(\varphi)\), where \(\Gamma\) are the x, y, z coordinates of the curve.

gammadashdash_impl(gammadashdash)

This function returns \(\Gamma''(\varphi)\), and \(\Gamma\) are the x, y, z coordinates of the curve.

gammadashdashdash_impl(gammadashdashdash)

This function returns \(\Gamma'''(\varphi)\), and \(\Gamma\) are the x, y, z coordinates of the curve.

incremental_arclength_impl()

This function returns the incremental arclength of the curve.

set_dofs(dofs)

This function sets the dofs of the curve.

class simsopt.geo.JaxCurvePlanarFourier(quadpoints, order, dofs=None)

Bases: JaxCurve

A Python+Jax implementation of the CurvePlanarFourier class. This is an autodiff compatible version of the same CurvePlanarFourier class in the C++ implementation in simsoptpp. The point of this class is to illustrate how jax can be used to define a geometric object class and calculate all the derivatives (both with respect to dofs and with respect to the angle \(\theta\)) automatically.

[r_{c,0}, cdots, r_{c,text{order}}, r_{s,1}, cdots, r_{s,text{order}}, q_0, q_i, q_j, q_k, x_{text{center}}, y_{text{center}}, z_{text{center}}]

Parameters:

• quadpoints (array) – Array of quadrature points.

• order (int) – Order of the Fourier series.

• dofs (array) – Array of dofs.

_make_names(order)

This function returns the names of the dofs associated to this object.

Parameters:

order (int) – Order of the Fourier series.

Returns:

List of dof names.

get_dofs()

This function returns the dofs associated to this object.

num_dofs()

This function returns the number of dofs associated to this object.

set_dofs_impl(dofs)

This function sets the dofs associated to this object.

class simsopt.geo.JaxCurveXYZFourier(quadpoints, order, dofs=None)

Bases: JaxCurve

A Python+Jax implementation of the CurveXYZFourier class. This is an autodiff compatible version of the same CurveXYZFourier class in the C++ implementation in simsoptpp. The point of this class is to illustrate how jax can be used to define a geometric object class and calculate all the derivatives (both with respect to dofs and with respect to the angle :math:` heta`) automatically.

Parameters:

• quadpoints (array) – Array of quadrature points.

• order (int) – Order of the Fourier series.

• dofs (array) – Array of dofs.

_make_names(order)

This function returns the names of the dofs associated to this object.

Parameters:

order (int) – Order of the Fourier series.

Returns:

List of dof names.

get_dofs()

This function returns the dofs associated to this object.

num_dofs()

This function returns the number of dofs associated to this object.

set_dofs_impl(dofs)

This function sets the dofs associated to this object.

class simsopt.geo.LPBinormalCurvatureStrainPenalty(framedcurve, width=0.001, p=2, threshold=0)

Bases: Optimizable

This class computes a penalty term based on the \(L_p\) norm of the binormal curvature strain, and penalizes where the local strain exceeds a threshold

\[J = \frac{1}{p} \int_{\text{curve}} \text{max}(\epsilon_{\text{bend}} - \epsilon_0, 0)^p ~dl,\]

where

\[\epsilon_{\text{bend}} = \frac{w |\hat{\textbf{b}} \cdot \boldsymbol{\kappa}|}{2},\]

\(w\) is the width of the tape, \(\hat{\textbf{b}}\) is the frame binormal vector, \(\boldsymbol{\kappa}\) is the curvature vector of the filamentary coil, and \(\epsilon_0\) is a threshold strain, given by the argument threshold.

J()

This returns the value of the quantity.

dJ(*args, partials=False, **kwargs)

return_fn_map: Dict[str, Callable] = {'J': <function LPBinormalCurvatureStrainPenalty.J>, 'dJ': <function derivative_dec.<locals>._derivative_dec>}

class simsopt.geo.LPTorsionalStrainPenalty(framedcurve, width=0.001, p=2, threshold=0)

Bases: Optimizable

This class computes a penalty term based on the \(L_p\) norm of the torsional strain, and penalizes where the local strain exceeds a threshold

\[J = \frac{1}{p} \int_{\text{curve}} \text{max}(\epsilon_{\text{tor}} - \epsilon_0, 0)^p ~dl\]

where

\[\epsilon_{\text{tor}} = \frac{\tau^2 w^2}{12},\]

\(\tau\) is the torsion of the tape frame, \(w\) is the width of the tape, and \(\epsilon_0\) is a threshold strain, given by the argument threshold.

J()

This returns the value of the quantity.

dJ(*args, partials=False, **kwargs)

return_fn_map: Dict[str, Callable] = {'J': <function LPTorsionalStrainPenalty.J>, 'dJ': <function derivative_dec.<locals>._derivative_dec>}

class simsopt.geo.LinkingNumber(curves, downsample=1)

Bases: Optimizable

J()

dJ(*args, partials=False, **kwargs)

dphis

Compute the Gauss linking number of a set of curves, i.e. whether the curves are interlocked or not.

The value is an integer, >= 1 if the curves are interlocked, 0 if not. For each pair of curves, the contribution to the linking number is

\[Link(c_1, c_2) = \frac{1}{4\pi} \left| \oint_{c_1}\oint_{c_2}\frac{\textbf{r}_1 - \textbf{r}_2}{|\textbf{r}_1 - \textbf{r}_2|^3} (d\textbf{r}_1 \times d\textbf{r}_2) \right|\]

where \(c_1\) is the first curve, \(c_2\) is the second curve, \(\textbf{r}_1\) is the position vector along the first curve, and \(\textbf{r}_2\) is the position vector along the second curve.

Parameters:

• curves (list of Curve, shape (n_curves)) – The set of curves for which the linking number should be computed.

• downsample (int, default=1) – Factor by which to downsample the quadrature points by skipping through the array by a factor of downsample, e.g. curve.gamma()[::downsample, :]. Setting this parameter to a value larger than 1 will speed up the calculation, which may be useful if the set of coils is large, though it may introduce inaccuracy if downsample is set too large, or not a multiple of the total number of quadrature points (since this will produce a nonuniform set of points). This parameter is used to speed up expensive calculations during optimization, while retaining higher accuracy for the other objectives.

class simsopt.geo.LpCurveCurvature(curve, p, threshold=0.0)

Bases: Optimizable

This class computes a penalty term based on the \(L_p\) norm of the curve’s curvature, and penalizes where the local curve curvature exceeds a threshold

\[J = \frac{1}{p} \int_{\text{curve}} \text{max}(\kappa - \kappa_0, 0)^p ~dl\]

where \(\kappa_0\) is a threshold curvature, given by the argument threshold.

J()

This returns the value of the quantity.

dJ(*args, partials=False, **kwargs)

return_fn_map: Dict[str, Callable] = {'J': <function LpCurveCurvature.J>, 'dJ': <function derivative_dec.<locals>._derivative_dec>}

class simsopt.geo.LpCurveTorsion(curve, p, threshold=0.0)

Bases: Optimizable

LpCurveTorsion is a class that computes a penalty term based on the \(L_p\) norm of the curve’s torsion:

\[J = \frac{1}{p} \int_{\text{curve}} \max(|\tau|-\tau_0, 0)^p ~dl.\]

J()

This returns the value of the quantity.

dJ(*args, partials=False, **kwargs)

return_fn_map: Dict[str, Callable] = {'J': <function LpCurveTorsion.J>, 'dJ': <function derivative_dec.<locals>._derivative_dec>}

class simsopt.geo.MajorRadius(boozer_surface)

Bases: Optimizable

This wrapper objective computes the major radius of a toroidal Boozer surface and supplies its derivative with respect to coils

Parameters:

boozer_surface – The surface to use for the computation

J()

compute()

dJ(*args, partials=False, **kwargs)

recompute_bell(parent=None)

Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.

This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.

Need to be implemented by classes that provide a dof_setter for external handling of DOFs.

class simsopt.geo.MeanSquaredCurvature(curve)

Bases: Optimizable

J()

__init__(curve)

Compute the mean of the squared curvature of a curve.

\[J = (1/L) \int_{\text{curve}} \kappa^2 ~dl\]

where \(L\) is the curve length, \(\ell\) is the incremental arclength, and \(\kappa\) is the curvature.

Parameters:

curve – the curve of which the curvature should be computed.

dJ(*args, partials=False, **kwargs)

class simsopt.geo.NonQuasiSymmetricRatio(boozer_surface, bs, sDIM=20, quasi_poloidal=False)

Bases: Optimizable

This objective decomposes the field magnitude \(B(\varphi,\theta)\) into quasisymmetric and non-quasisymmetric components. For quasi-axisymmetry, we compute

\[\begin{split}B_{\text{QS}} &= \frac{\int_0^1 B \|\mathbf n\| ~d\varphi}{\int_0^1 \|\mathbf n\| ~d\varphi} \\ B_{\text{non-QS}} &= B - B_{\text{QS}}\end{split}\]

where \(B = \| \mathbf B(\varphi,\theta) \|_2\). For quasi-poloidal symmetry, an analagous formula is used, but the integral is computed in the \(\theta\) direction. The objective computed by this penalty is

\[\begin{split}J &= \frac{\int_{\Gamma_{s}} B_{\text{non-QS}}^2~dS}{\int_{\Gamma_{s}} B_{\text{QS}}^2~dS} \\\end{split}\]

When \(J\) is zero, then there is perfect QS on the given boozer surface. The ratio of the QS and non-QS components of the field is returned to avoid dependence on the magnitude of the field strength. Note that this penalty is computed on an auxilliary surface with quadrature points that are different from those on the input Boozer surface. This is to allow for a spectrally accurate evaluation of the above integrals. Note that if boozer_surface.surface.stellsym == True, computing this term on the half-period with shifted quadrature points is ~not~ equivalent to computing on the full-period with unshifted points. This is why we compute on an auxilliary surface with quadrature points on the full period.

Parameters:

• boozer_surface – input boozer surface on which the penalty term is evaluated,

• biotsavart – biotsavart object (not necessarily the same as the one used on the Boozer surface).

• sDIM – integer that determines the resolution of the quadrature points placed on the auxilliary surface.

• quasi_poloidal – False for quasiaxisymmetry and True for quasipoloidal symmetry

J()

compute()

dJ(*args, partials=False, **kwargs)

dJ_by_dB()

Return the partial derivative of the objective with respect to the magnetic field

dJ_by_dsurfacecoefficients()

Return the partial derivative of the objective with respect to the surface coefficients

recompute_bell(parent=None)

Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.

This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.

Need to be implemented by classes that provide a dof_setter for external handling of DOFs.

class simsopt.geo.PermanentMagnetGrid(plasma_boundary: Surface, Bn, coordinate_flag='cartesian')

Bases: object

PermanentMagnetGrid is a class for setting up the grid, plasma surface, and other objects needed to perform permanent magnet optimization for stellarators. The class takes as input two toroidal surfaces specified as SurfaceRZFourier objects, and initializes a set of points (in cylindrical coordinates) between these surfaces. If an existing FAMUS grid file called from FAMUS is desired, use from_famus() instead. It finishes initialization by pre-computing a number of quantities required for the optimization such as the geometric factor in the dipole part of the magnetic field and the target Bfield that is a sum of the coil and plasma magnetic fields.

Parameters:

• plasma_boundary – Surface class object Representing the plasma boundary surface. Gets converted into SurfaceRZFourier object for ease of use.

• Bn – 2D numpy array, shape (ntheta_quadpoints, nphi_quadpoints) Magnetic field (coils and plasma) at the plasma boundary. Typically this will be the optimized plasma magnetic field from a stage-1 optimization, and the optimized coils from a basic stage-2 optimization. This variable must be specified to run the permanent magnet optimization.

• coordinate_flag – string Flag to specify the coordinate system used for the grid and optimization. This is primarily used to tell the optimizer which coordinate directions should be considered, “grid-aligned”. Therefore, you can do weird stuff like cartesian grid-alignment on a simple toroidal grid, which might be self-defeating. However, if coordinate_flag=’cartesian’ a rectangular Cartesian grid is initialized using the Nx, Ny, and Nz parameters. If the coordinate_flag=’cylindrical’, a uniform cylindrical grid is initialized using the dr and dz parameters.

_print_initial_opt()

Print out initial errors and the bulk optimization parameters before the permanent magnets are optimized.

_setup_uniform_grid()

Initializes a uniform grid in cartesian coordinates and sets some important grid variables for later.

coordinate_flag

Validates that the input is among the specified strings.

Parameters:

*args – All the acceptable values for the string

classmethod geo_setup_between_toroidal_surfaces(plasma_boundary, Bn, inner_toroidal_surface: Surface, outer_toroidal_surface: Surface, **kwargs)

Function to initialize a SIMSOPT PermanentMagnetGrid from a volume defined by two toroidal surfaces. These must be specified directly. Often a good choice is made by extending the plasma boundary by its normal vectors.

Parameters:

• inner_toroidal_surface (Surface class object) – Representing the inner toroidal surface of the volume. Gets converted into SurfaceRZFourier object for ease of use.

• outer_toroidal_surface (Surface object representing) – the outer toroidal surface of the volume. Typically want this to have same quadrature points as the inner surface for a functional grid setup. Gets converted into SurfaceRZFourier object for ease of use.

• kwargs (The following are valid keyword arguments.) –

  m_maxima: float or 1D numpy array, shape (Ndipoles)

  Optional array of maximal dipole magnitudes for the permanent magnets. If not provided, defaults to the common magnets used in the MUSE design, with strengths ~ 1 Tesla.

  pol_vectors: 3D numpy array, shape (Ndipoles, Ncoords, 3)

  Optional set of local coordinate systems for each dipole, which specifies which directions should be considered grid-aligned. Ncoords can be > 3, as in the PM4Stell design.

  dr: double

  Radial grid spacing in the permanent magnet manifold. Used only if coordinate_flag = cylindrical, then dr is the radial size of the cylindrical bricks in the grid.

  dz: double

  Axial grid spacing in the permanent magnet manifold. Used only if coordinate_flag = cylindrical, then dz is the axial size of the cylindrical bricks in the grid.

  Nx: int

  Number of points in x to use in a cartesian grid, taken between the inner and outer toroidal surfaces. Used only if the coordinate_flag = cartesian, then Nx is the x-size of the rectangular cubes in the grid.

  Ny: int

  Number of points in y to use in a cartesian grid, taken between the inner and outer toroidal surfaces. Used only if the coordinate_flag = cartesian, then Ny is the y-size of the rectangular cubes in the grid.

  Nz: int

  Number of points in z to use in a cartesian grid, taken between the inner and outer toroidal surfaces. Used only if the coordinate_flag = cartesian, then Nz is the z-size of the rectangular cubes in the grid.

  coordinate_flag: string

  Flag to specify the coordinate system used for the grid and optimization. This is primarily used to tell the optimizer which coordinate directions should be considered, “grid-aligned”. Therefore, you can do weird stuff like cartesian grid-alignment on a simple toroidal grid, which might be self-defeating. However, if coordinate_flag=’cartesian’ a rectangular Cartesian grid is initialized using the Nx, Ny, and Nz parameters. If the coordinate_flag=’cylindrical’, a uniform cylindrical grid is initialized using the dr and dz parameters.

Returns:

pm_grid

Return type:

An initialized PermanentMagnetGrid class object.

classmethod geo_setup_from_famus(plasma_boundary, Bn, famus_filename, **kwargs)

Function to initialize a SIMSOPT PermanentMagnetGrid from a pre-existing FAMUS file defining a grid of magnets. For example, this function is used for the half-Tesla NCSX configuration (c09r000) and MUSE grids.

Parameters:

• famus_filename (string) – Filename of a FAMUS grid file (a pre-made dipole grid).

• kwargs (The following are valid keyword arguments.) –

  m_maxima: float or 1D numpy array, shape (Ndipoles)

  Optional array of maximal dipole magnitudes for the permanent magnets. If not provided, defaults to the common magnets used in the MUSE design, with strengths ~ 1 Tesla.

  pol_vectors: 3D numpy array, shape (Ndipoles, Ncoords, 3)

  Optional set of local coordinate systems for each dipole, which specifies which directions should be considered grid-aligned. Ncoords can be > 3, as in the PM4Stell design.

  coordinate_flag: string

  Flag to specify the coordinate system used for the grid and optimization. This is primarily used to tell the optimizer which coordinate directions should be considered, “grid-aligned”. Therefore, you can do weird stuff like cartesian grid-alignment on a simple toroidal grid, which might be self-defeating. However, if coordinate_flag=’cartesian’ a rectangular Cartesian grid is initialized using the Nx, Ny, and Nz parameters. If the coordinate_flag=’cylindrical’, a uniform cylindrical grid is initialized using the dr and dz parameters.

  downsample: int

  Optional integer for downsampling the FAMUS grid, since the MUSE and other grids can be very high resolution and this makes CI take a long while.

Returns:

pm_grid

Return type:

An initialized PermanentMagnetGrid class object.

rescale_for_opt(reg_l0, reg_l1, reg_l2, nu)

Scale regularizers to the largest scale of ATA (~1e-6) to avoid regularization >> ||Am - b|| ** 2 term in the optimization. The prox operator uses reg_l0 * nu for the threshold so normalization below allows reg_l0 and reg_l1 values to be exactly the thresholds used in calculation of the prox. Then add contributions to ATA and ATb coming from extra loss terms such as L2 regularization and relax-and-split. Currently does not rescale the L2 and nu hyperparameters, but users may want to play around with this.

Parameters:

• reg_l0 – L0 regularization.

• reg_l1 – L1 regularization.

• reg_l2 – L2 regularization.

• nu – nu hyperparameter in relax-and-split optimization.

Returns:

Rescaled L0 regularization. reg_l1: Rescaled L1 regularization. reg_l2: Rescaled L2 regularization. nu: Rescaled nu hyperparameter in relax-and-split optimization.

Return type:

reg_l0

write_to_famus(out_dir='')

Takes a PermanentMagnetGrid object and saves the geometry and optimization solution into a FAMUS input file.

Parameters:

out_dir – Path object for the output directory for saved files.

class simsopt.geo.PerturbationSample(sampler, randomgen=None, sample=None)

Bases: GSONable

This class represents a single sample of a perturbation. The point of having a dedicated class for this is so that we can apply the same perturbation to multipe curves (e.g. in the case of multifilament approximations to finite build coils). The main way to interact with this class is via the overloaded __getitem__ (i.e. [ ] indexing). For example:

sample = PerturbationSample(...)
g = sample[0] # get the values of the perturbation
gd = sample[1] # get the first derivative of the perturbation

__getitem__(deriv)

Get the perturbation (if deriv=0) or its deriv-th derivative.

resample()

class simsopt.geo.Port

Bases: ABC

Abstract base class for ports

abstract collides()

Determines if the user input point(s) collide with the port, with an optional gap spacing enforced around the port’s exterior.

abstract mesh_representation()

Constructs a triangular mesh representation of the port for plotting and visualization.

Returns:

• x, y, z (1D arrays) – Cartesian x, y, and z coordinates of the vertices of the mesh.

• triangles (2D array) – Indices of the vertices (as listed in the x, y, and z arrays) surrounding each triangular face within the mesh. Each row (dimension 1) represents a face.

abstract plot()

Returns handle to a three-dimensional plot with a visual depiction of the port.

abstract repeat_via_symmetries()

Returns a PortSet class instance containing the calling Port class instances as well ports with parameters transformed to equivalent locations in other periods (and half-periods if stellarator symmetry is assumed).

class simsopt.geo.PortSet(ports=None, file=None, port_type=None)

Bases: object

Handles sets of port-like objects.

__add__(p)

Allow new port sets to be built with the “+” operator

__getitem__(key)

Allow member ports to be accessed with “[]” directly from the class instance

__init__(ports=None, file=None, port_type=None)

Initializes a list of ports by loading parameters from a file and/or incorporating input port class instances.

Parameters:

• ports (list of port class instances (optional)) – Ports to be added to the set

• file (string or list of strings (optional)) – Name(s) of the file(s) from which to load port data of type specified by port_type

• port_type (string or list of strings (optional)) – Type(s) of ports to load from file(s) if file is specified. Supported inputs currently include ‘circular’, ‘rectangular’, and the corresponding Python type instances.

add_ports(ports)

Add ports to the set

Parameters:

ports (list of port class instances)

collides(x, y, z, gap=0.0)

Determines if the user input point(s) collide with any port in the set.

Parameters:

• x (float arrays) – Cartesian x, y, and z coordinates of the point(s) to be assessed. Dimensions must be the same.

• y (float arrays) – Cartesian x, y, and z coordinates of the point(s) to be assessed. Dimensions must be the same.

• z (float arrays) – Cartesian x, y, and z coordinates of the point(s) to be assessed. Dimensions must be the same.

• gap (float or array-like (optional)) – Minimum gap spacing(s) required between the point and the exteriors of each port. If scalar, the same gap will be applied to all ports; if a list or array, the n-th element will be assigned to the n-th port in the set. Default is a scalar 0.

Returns:

colliding – Array of logical values indicating whether each of the input points collide with any of the ports in the set.

Return type:

logical array

load_circular_ports_from_file(file)

Loads a set of circular ports from a file. The file must have the CSV (comma-separated values) format and exactly one header line (the first line in the file will be ignored).

The columns must have the following parameters in order (see docstring for CircularPort for definitions): ox, oy, oz, ax, ay, az, ir, thick, l0, l1

Parameters:

file (string) – Name of the file from which to load the port parameters

load_rectangular_ports_from_file(file)

Loads a set of rectangular ports from a file. The file must have the CSV (comma-separated values) format and exactly one header line (the first line in the file will be ignored).

The columns must have the following parameters in order (see docstring for RectangularPort for definitions): ox, oy, oz, ax, ay, az, wx, wy, wz, iw, ih, thick, l0, l1

Parameters:

file (string) – Name of the file from which to load the port parameters

plot(n_edges=100, **kwargs)

Places representations of the ports on a 3D plot. Currently only works with the mayavi plotting package.

Parameters:

• n_edges (integer (optional)) – Number of edges for the polygons used to approximate the cross-section for any ports that are circular. Default is 100.

• **kwargs – Keyword arguments to be passed to the mayavi surface module for each port

Returns:

surfs – References to the plotted ports

Return type:

list of mayavi.modules.surface.Surface class instance

repeat_via_symmetries(nfp, stell_sym)

Creates a new set that contains the ports in the initial set plus additional equivalent ports the remaining field periods/half-periods.

Parameters:

• nfp (integer) – Number of toroidal field periods

• stell_sym (logical) – If true, stellarator symmetry will be assumed and equivalent ports will be created for every half-period

Returns:

ports_out – A set of all the symmetric ports, including the ones represented by the calling PortSet class instance

Return type:

PortSet class instance

save_ports_to_file(fname)

Save data on ports to csv-formatted files. Separate files will be created for each port type. Circular ports will be saved to a file ending in “_circ.csv”; rectangular ports will be saved to a file ending in “_rect.csv”.

Parameters:

fname (string) – Name of the file to save, not including the suffix

to_vtk(filename, n_edges=100)

Export a mesh representation of the port set to a VTK file, which can be read with ParaView.

In the VTK file, each mesh point will be associated with an integer value “index” corresponding to the index of its respective port within the port set.

Note: This function requires the pyevtk python package, which can be installed using pip install pyevtk.

Parameters:

• filename (string) – Name of the VTK file, without extension, to create.

• n_edges (integer (optional)) – Number of edges for the polygon used to approximate the circular cross-section of any circular ports. Default is 100.

class simsopt.geo.PrincipalCurvature(surface, kappamax1=1, kappamax2=1, weight1=0.05, weight2=0.05)

Bases: Optimizable

Given a Surface, evaluates a metric based on the principal curvatures, \(\kappa_1\) and \(\kappa_2\), where \(\kappa_1>\kappa_2\). This metric is designed to penalize \(\kappa_1 > \kappa_{\max,1}\) and \(-\kappa_2 > \kappa_{\max,2}\).

\[\begin{split}J &= \int d^2 x \exp \left(- ( \kappa_1 - \kappa_{\max,1})/w_1) \right) \\ &+ \int d^2 x \exp \left(- (-\kappa_2 - \kappa_{\max,2})/w_2) \right).\end{split}\]

This metric can be used as a regularization within fixed-boundary optimization to prevent, for example, surfaces with concave regions (large values of \(|\kappa_2|\)) or surfaces with large elongation (large values of \(\kappa_1\)).

J()

dJ(*args, partials=False, **kwargs)

class simsopt.geo.QfmResidual(surface, biotsavart)

Bases: Optimizable

For a given surface \(S\), this class computes the residual

\[f(S) = \frac{\int_{S} d^2 x \, (\textbf{B} \cdot \hat{\textbf{n}})^2}{\int_{S} d^2 x \, B^2}\]

where \(\textbf{B}\) is the magnetic field from biotsavart, \(\hat{\textbf{n}}\) is the unit normal on a given surface, and the integration is performed over the surface. Derivatives are computed wrt the surface dofs.

J()

dJ_by_dsurfacecoefficients()

Calculate the derivatives with respect to the surface coefficients

invalidate_cache()

recompute_bell(parent=None)

Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.

This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.

Need to be implemented by classes that provide a dof_setter for external handling of DOFs.

class simsopt.geo.QfmSurface(biotsavart, surface, label, targetlabel)

Bases: GSONable

QfmSurface is used to compute a quadratic-flux minimizing surface, defined as the minimum of the objective function,

\[f(S) = \frac{\int_{S} d^2 x \, \left(\textbf{B} \cdot \hat{\textbf{n}}\right)^2}{\int_{S} d^2 x \, B^2}\]

subject to a constraint on the surface label, such as the area, volume, or toroidal flux. This objective, \(f\), is computed in QfmResidual (see surfaceobjectives.py). If magnetic surfaces exists, \(f\) will be zero. If not, QFM surfaces approximate flux surfaces.

The label constraint can be enforced with a penalty formulation, defined in qfm_penalty_constraints(), and whose minimium is computed with LBFGS-B by minimize_qfm_penalty_constraints_LBFGS().

Alternatively, the label constraint can be enforced with a constrained optimization algorithm. This constrained optimization problem is solved with the SLSQP algorithm by minimize_qfm_exact_constraints_SLSQP().

minimize_qfm(tol=0.001, maxiter=1000, method='SLSQP', constraint_weight=1.0)

This function tries to find the surface that approximately solves

\[\min_{S} f(S)\]

subject to

\[\texttt{label} = \texttt{labeltarget}\]

where \(f(S)\) is the QFM residual. This is done using method, either ‘SLSQP’ or ‘LBFGS’. If ‘LBFGS’ is chosen, a penalized formulation is used defined by constraint_weight. If ‘SLSQP’ is chosen, constraint_weight is not used, and the constrained optimization problem is solved.

minimize_qfm_exact_constraints_SLSQP(tol=0.001, maxiter=1000)

This function tries to find the surface that approximately solves

\[\min_{S} f(S)\]

subject to

\[\texttt{label} = \texttt{labeltarget}\]

where \(f(S)\) is the QFM residual. This is done using SLSQP.

minimize_qfm_penalty_constraints_LBFGS(tol=0.001, maxiter=1000, constraint_weight=1.0)

This function tries to find the surface that approximately solves

\[\min_S f(S) + 0.5 \texttt{constraint_weight} (\texttt{label} - \texttt{labeltarget})^2\]

where \(f(S)\) is the QFM residual defined in surfaceobjectives.py. This is done using LBFGS.

qfm_label_constraint(x, derivatives=0)

Returns the residual

\[0.5 (\texttt{label} - \texttt{labeltarget})^2\]

and optionally the gradient wrt surface params.

qfm_objective(x, derivatives=0)

Returns the residual

\[f(S)\]

and optionally the gradient wrt surface params where \(f(S)`\) is the QFM residual defined in surfaceobjectives.py.

qfm_penalty_constraints(x, derivatives=0, constraint_weight=1)

Returns the residual

\[f(S) + 0.5 \texttt{constraint_weight} (\texttt{label} - \texttt{labeltarget})^2\]

and optionally the gradient and Hessian wrt surface params where \(f(S)\) is the QFM residual defined in surfaceobjectives.py.

class simsopt.geo.RectangularPort(ox=1.0, oy=0.0, oz=0.0, ax=1.0, ay=0.0, az=0.0, wx=0.0, wy=1.0, wz=0.0, iw=0.5, ih=0.5, thick=0.01, l0=0.0, l1=0.5)

Bases: Port

Class representing a port with a rectangular cross-section.

__init__(ox=1.0, oy=0.0, oz=0.0, ax=1.0, ay=0.0, az=0.0, wx=0.0, wy=1.0, wz=0.0, iw=0.5, ih=0.5, thick=0.01, l0=0.0, l1=0.5)

Initializes the RectangularPort class according to the geometric port parameters.

Parameters:

• ox (floats) – Cartesian x, y, and z coordinates of the origin point of the port in 3d space

• oy (floats) – Cartesian x, y, and z coordinates of the origin point of the port in 3d space

• oz (floats) – Cartesian x, y, and z coordinates of the origin point of the port in 3d space

• ax (floats) – Cartesian x, y, and z components of a vector in the direction of the port’s axis

• ay (floats) – Cartesian x, y, and z components of a vector in the direction of the port’s axis

• az (floats) – Cartesian x, y, and z components of a vector in the direction of the port’s axis

• wx (floats) – Cartesian x, y, and z components of a vector in the direction spanning the width of the cross-section, assumed perpendicular to the axis

• wy (floats) – Cartesian x, y, and z components of a vector in the direction spanning the width of the cross-section, assumed perpendicular to the axis

• wz (floats) – Cartesian x, y, and z components of a vector in the direction spanning the width of the cross-section, assumed perpendicular to the axis

• iw (float) – Inner width of the cross-section, i.e. the dimension spanned by the vector (wx, wy, wz)

• ih (float) – Inner height of the cross-section, i.e. the dimension spanned by the vector (hx, hy, hz)

• thick (float) – Thickness of the port wall

• l0 (double) – Locations of the ends of the port, specified by the signed distance along the port axis from the origin point

• l1 (double) – Locations of the ends of the port, specified by the signed distance along the port axis from the origin point

collides(x, y, z, gap=0.0)

Determines if the user input point(s) collide with the port

Parameters:

• x (float arrays) – Cartesian x, y, and z coordinates of the point(s) to be assessed. Dimensions must be the same.

• y (float arrays) – Cartesian x, y, and z coordinates of the point(s) to be assessed. Dimensions must be the same.

• z (float arrays) – Cartesian x, y, and z coordinates of the point(s) to be assessed. Dimensions must be the same.

• gap (float (optional)) – Minimum gap spacing required between the point and the exterior of the port. Default is 0.

Returns:

colliding – Array of logical values indicating whether each of the input points collide with the port.

Return type:

logical array

mesh_representation()

Constructs a triangular mesh representation of the port for plotting and visualization.

Returns:

• x, y, z (1D arrays) – Cartesian x, y, and z coordinates of the vertices of the mesh.

• triangles (2D array) – Indices of the vertices (as listed in the x, y, and z arrays) surrounding each triangular face within the mesh. Each row (dimension 1) represents a face.

plot(**kwargs)

Places a representation of the port on a 3D plot. Currently only works with the mayavi plotting package.

Parameters:

**kwargs – Keyword arguments to be passed to the mayavi surface module for the port

Returns:

surf – Reference to the plotted port

Return type:

mayavi.modules.surface.Surface class instance

repeat_via_symmetries(nfp, stell_sym, include_self=True)

Generates a set of equivalent ports to uphold toroidal and/or stellarator symmetry.

Parameters:

• nfp (integer) – Number of toroidal field periods

• stell_sym (logical) – If true, stellarator symmetry will be assumed and equivalent ports will be created for every half-period

• include_self (logical (optional)) – If true, the port instance calling this method will be included in the returned list of ports (see below). Default is True.

Returns:

ports – PortSet containing the equivalent ports, including the one represented by the calling instance if include_self == True.

Return type:

PortSet class instance

to_vtk(filename)

Export a mesh representation of the port to a VTK file, which can be read with ParaView.

Note: This function requires the pyevtk python package, which can be installed using pip install pyevtk.

Parameters:

filename (string) – Name of the VTK file, without extension, to create.

class simsopt.geo.RotatedCurve(curve, phi, flip)

Bases: Curve, Curve

RotatedCurve inherits from the Curve base class. It takes an input a Curve, rotates it about the z axis by a toroidal angle phi, and optionally completes a reflection when flip=True.

dgamma_by_dcoeff_impl(dgamma_by_dcoeff)

This function returns

\[\frac{\partial \Gamma}{\partial \mathbf c}\]

where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.

dgamma_by_dcoeff_vjp(v)

This function returns the vector Jacobian product

\[v^T \frac{\partial \Gamma}{\partial \mathbf c}\]

where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.

dgammadash_by_dcoeff_impl(dgammadash_by_dcoeff)

This function returns

\[\frac{\partial \Gamma'}{\partial \mathbf c}\]

where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.

dgammadash_by_dcoeff_vjp(v)

This function returns the vector Jacobian product

\[v^T \frac{\partial \Gamma'}{\partial \mathbf c}\]

where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.

dgammadashdash_by_dcoeff_impl(dgammadashdash_by_dcoeff)

This function returns

\[\frac{\partial \Gamma''}{\partial \mathbf c}\]

where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.

dgammadashdash_by_dcoeff_vjp(v)

This function returns the vector Jacobian product

\[v^T \frac{\partial \Gamma''}{\partial \mathbf c}\]

where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.

dgammadashdashdash_by_dcoeff_impl(dgammadashdashdash_by_dcoeff)

This function returns

\[\frac{\partial \Gamma'''}{\partial \mathbf c}\]

where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.

dgammadashdashdash_by_dcoeff_vjp(v)

This function returns the vector Jacobian product

\[v^T \frac{\partial \Gamma'''}{\partial \mathbf c}\]

where \(\mathbf{c}\) are the curve dofs, and \(\Gamma\) are the x, y, z coordinates of the curve.

property flip

gamma_impl(gamma, quadpoints)

This function returns the x,y,z coordinates of the curve, \(\Gamma\), where \(\Gamma\) are the x, y, z coordinates of the curve.

gammadash_impl(gammadash)

This function returns \(\Gamma'(\varphi)\), where \(\Gamma\) are the x, y, z coordinates of the curve.

gammadashdash_impl(gammadashdash)

This function returns \(\Gamma''(\varphi)\), where \(\Gamma\) are the x, y, z coordinates of the curve.

gammadashdashdash_impl(gammadashdashdash)

This function returns \(\Gamma'''(\varphi)\), where \(\Gamma\) are the x, y, z coordinates of the curve.

get_dofs()

RotatedCurve does not have any dofs of its own. This function returns null array

num_dofs()

This function returns the number of dofs associated to the curve.

set_dofs_impl(d)

RotatedCurve does not have any dofs of its own. This function does nothing.

class simsopt.geo.Surface(**kwargs)

Bases: Optimizable

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

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_FIELD_PERIOD = 'field period'

RANGE_FULL_TORUS = 'full torus'

RANGE_HALF_PERIOD = 'half period'

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

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.

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.

Parameters:

• 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.

d2aspect_ratio_by_dcoeff_dcoeff()

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

d2major_radius_by_dcoeff_dcoeff()

Return the second derivative of the major radius wrt surface coefficients

d2mean_cross_sectional_area_by_dcoeff_dcoeff()

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

d2minor_radius_by_dcoeff_dcoeff()

Return the second derivative of the minor radius wrt surface coefficients

daspect_ratio_by_dcoeff()

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

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

dmajor_radius_by_dcoeff()

Return the derivative of the major radius wrt surface coefficients

dmean_cross_sectional_area_by_dcoeff()

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

dminor_radius_by_dcoeff()

Return the derivative of the minor radius wrt surface coefficients

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.

Parameters:

• 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_phi_quadpoints(nphi=None, range=None, nfp=1)

Sets the phi grid points for Surface subclasses.

Parameters:

• 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\).

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.

Parameters:

• 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.

Parameters:

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

Returns:

List of grid points \(\theta_j\).

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:

2d array (numquadpoints_phi,numquadpoints_theta)

defining interpolated function on arclength angle along curve at constant phi

Return type:

function_interpolated

is_self_intersecting(angle=0.0, thetas=None)

This function computes a cross section of self at the input cylindrical angle. Then, approximating the cross section as a piecewise linear polyline, the Bentley-Ottmann algorithm is used to check for self-intersecting segments of the cross section. NOTE: if this function returns False, the surface may still be self-intersecting away from angle.

Parameters:

• angle (float) – the cylindrical angle at which we would like to check whether the surface is self-intersecting. Note that a surface might not be self-intersecting at a given angle, but may be self-intersecting elsewhere. To be certain that the surface is not self-intersecting, it is recommended to run this check at multiple angles. Also note that angle is assumed to be in radians, and not divided by 2*pi.

• thetas (int, array, optional) – can be (1) an integer indicating that the cross section should be calculated at the poloidal angles in the array np.linspace(0, 1, thetas, endpoint=False), or (2) an array containing the poloidal positions at which the cross section should be computed, between 0 and 1. If thetas is not provided, then the cross section will be calculated at the positions given by self.quadpoints_theta.

Returns:

True if surface is self-intersecting at angle, else False.

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.

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.

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.

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.

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.

Parameters:

• 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.

abstract to_RZFourier()

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

to_vtk(filename, extra_data=None)

Export the surface to a VTK format file, which can be read with Paraview. This function requires the pyevtk python package, which can be installed using pip install pyevtk.

Parameters:

• filename (str) – Name of the file to write

• extra_data (dict) – An optional data field (dictionary) on the surface, which can be associated with a colormap in Paraview.

class simsopt.geo.SurfaceClassifier(surface, p=1, h=0.05)

Bases: object

Takes in a toroidal surface and constructs an interpolant of the signed distance function \(f:R^3\to R\) that is positive inside the volume contained by the surface, (approximately) zero on the surface, and negative outisde the volume contained by the surface.

__init__(surface, p=1, h=0.05)

Parameters:

• surface – the surface to contruct the distance from.

• p – degree of the interpolant

• h – grid resolution of the interpolant

evaluate_rphiz(rphiz)

evaluate_xyz(xyz)

to_vtk(filename, h=0.01)

class simsopt.geo.SurfaceGarabedian(nfp=1, mmax=1, mmin=0, nmax=0, nmin=None, quadpoints_phi=None, quadpoints_theta=None, dofs=None)

Bases: Surface, Surface

SurfaceGarabedian represents a toroidal surface for which the shape is parameterized using Garabedian’s \(\Delta_{m,n}\) coefficients:

\[R + i Z = e^{i u} \sum_{m = m_\min}^{m_\max} \sum_{n = n_\min}^{n_\max} \Delta_{m,n} e^{-i m u + i n v}\]

where \(u = 2 \pi \theta\) is a poloidal angle on \([0, 2\pi]\), and \(v\) is the standard toroidal angle on \([0, 2\pi]\).

The present implementation assumes stellarator symmetry. Note that non-stellarator-symmetric surfaces require that the \(\Delta_{m,n}\) coefficients be imaginary.

For more information about the arguments quadpoints_phi, and quadpoints_theta, see the general documentation on Surfaces. Instead of supplying the quadrature point arrays along \(\phi\) and \(\theta\) directions, one could also specify the number of quadrature points for \(\phi\) and \(\theta\) using the class method from_nphi_ntheta().

Parameters:

• nfp – The number of field periods.

• mmin – Minimum poloidal mode number \(m\) included (usually 0 or negative).

• mmax – Maximum poloidal mode number \(m\) included.

• nmin – Minimum toroidal mode number \(n\) included (usually negative). If None, nmin = -nmax will be used.

• nmax – Maximum toroidal mode number \(n\) included.

• quadpoints_phi – Set this to a list or 1D array to set the \(\phi_j\) grid points directly.

• quadpoints_theta – Set this to a list or 1D array to set the \(\theta_j\) grid points directly.

property Delta

area()

Return the area of the surface.

area_volume()

Compute the surface area and the volume enclosed by the surface.

fix_range(mmin, mmax, nmin, nmax, fixed=True)

Fix the DOFs for a range of m and n values.

All modes with m in the interval [mmin, mmax] and n in the interval [nmin, nmax] will have their fixed property set to the value of the ‘fixed’ parameter. Note that mmax and nmax are included (unlike the upper bound in python’s range(min, max).)

classmethod from_RZFourier(surf)

Create a SurfaceGarabedian from a SurfaceRZFourier object of the identical shape.

For a derivation of the transformation here, see https://terpconnect.umd.edu/~mattland/assets/notes/toroidal_surface_parameterizations.pdf

get_Delta(m, n)

Return a particular \(\Delta_{m,n}\) coefficient.

get_dofs()

Return a 1D numpy array with all the degrees of freedom.

mmax

Validates that the input number is an integer. Optionally the user can specify the bounds for the number.

mmin

Validates that the input number is an integer. Optionally the user can specify the bounds for the number.

nfp

Validates that the input number is an integer. Optionally the user can specify the bounds for the number.

nmax

Validates that the input number is an integer. Optionally the user can specify the bounds for the number.

nmin

Validates that the input number is an integer. Optionally the user can specify the bounds for the number.

return_fn_map: Dict[str, Callable] = {'area': <function SurfaceGarabedian.area>, 'aspect-ratio': <function Surface.aspect_ratio>, 'volume': <function SurfaceGarabedian.volume>}

set_Delta(m, n, val)

Set a particular \(\Delta_{m,n}\) coefficient.

set_dofs(x)

Set the shape coefficients from a 1D list/array

to_RZFourier()

Return a SurfaceRZFourier object with the identical shape.

For a derivation of the transformation here, see https://terpconnect.umd.edu/~mattland/assets/notes/toroidal_surface_parameterizations.pdf

volume()

Return the volume of the surface.

class simsopt.geo.SurfaceHenneberg(nfp: int = 1, alpha_fac: int = 1, mmax: int = 1, nmax: int = 0, quadpoints_phi: Sequence[Real] | ndarray[tuple[Any, ...], dtype[float64]] = None, quadpoints_theta: Sequence[Real] | ndarray[tuple[Any, ...], dtype[float64]] = None, dofs: DOFs = None)

Bases: Surface, Surface

This class represents a toroidal surface using the parameterization in Henneberg, Helander, and Drevlak, Journal of Plasma Physics 87, 905870503 (2021). The main benefit of this representation is that there is no freedom in the poloidal angle, i.e. \(\theta\) is uniquely defined, in contrast to other parameterizations like SurfaceRZFourier. Stellarator symmetry is assumed.

In this representation by Henneberg et al, the cylindrical coordinates \((R,\phi,Z)\) are written in terms of a unique poloidal angle \(\theta\) as follows:

\[\begin{split}R(\theta,\phi) = R_0^H(\phi) + \rho(\theta,\phi) \cos(\alpha\phi) - \zeta(\theta,\phi) \sin(\alpha\phi), \\ Z(\theta,\phi) = Z_0^H(\phi) + \rho(\theta,\phi) \sin(\alpha\phi) + \zeta(\theta,\phi) \cos(\alpha\phi),\end{split}\]

where

\[\begin{split}R_0^H(\phi) &=& \sum_{n=0}^{nmax} R_{0,n}^H \cos(n_{fp} n \phi), \\ Z_0^H(\phi) &=& \sum_{n=1}^{nmax} Z_{0,n}^H \sin(n_{fp} n \phi), \\ \zeta(\theta,\phi) &=& \sum_{n=0}^{nmax} b_n \cos(n_{fp} n \phi) \sin(\theta - \alpha \phi), \\ \rho(\theta,\phi) &=& \sum_{n,m} \rho_{n,m} \cos(m \theta - n_{fp} n \phi - \alpha \phi).\end{split}\]

The continuous degrees of freedom are \(\{\rho_{m,n}, b_n, R_{0,n}^H, Z_{0,n}^H\}\). These variables correspond to the attributes rhomn, bn, R0nH, and Z0nH respectively, which are all numpy arrays. There is also a discrete degree of freedom \(\alpha\) which should be \(\pm n_{fp}/2\) where \(n_{fp}\) is the number of field periods. The attribute alpha_fac corresponds to \(2\alpha/n_{fp}\), so alpha_fac is either 1, 0, or -1. Using alpha_fac = 0 is appropriate for axisymmetry, while values of 1 or -1 are appropriate for a stellarator, depending on the handedness of the rotating elongation.

For \(R_{0,n}^H\) and \(b_n\), \(n\) is 0 or any positive integer up through nmax (inclusive). For \(Z_{0,n}^H\), \(n\) is any positive integer up through nmax. For \(\rho_{m,n}\), \(m\) is an integer from 0 through mmax (inclusive). For positive values of \(m\), \(n\) can be any integer from -nmax through nmax. For \(m=0\), \(n\) is restricted to integers from 1 through nmax. Note that we exclude the element of \(\rho_{m,n}\) with \(m=n=0\), because this degree of freedom is already represented in \(R_{0,0}^H\).

For the 2D array rhomn, functions set_rhomn() and get_rhomn() are provided for convenience so you can specify n, since the corresponding array index is shifted by nmax. There are no corresponding functions for the 1D arrays R0nH, Z0nH, and bn since these arrays all have a first index corresponding to n=0.

For more information about the arguments quadpoints_phi, and quadpoints_theta, see the general documentation on Surfaces. Instead of supplying the quadrature point arrays along \(\phi\) and \(\theta\) directions, one could also specify the number of quadrature points for \(\phi\) and \(\theta\) using the class method from_nphi_ntheta().

Parameters:

• nfp – The number of field periods.

• alpha_fac – Should be +1 or -1 for a stellarator, depending on the handedness by which the elongation rotates, or 0 for axisymmetry.

• mmax – Maximum poloidal mode number included.

• nmax – Maximum toroidal mode number included, divided by nfp.

• quadpoints_phi – Set this to a list or 1D array to set the \(\phi_j\) grid points directly.

• quadpoints_theta – Set this to a list or 1D array to set the \(\theta_j\) grid points directly.

_validate_mn(m, n)

Check whether given (m, n) values are allowed for \(\rho_{m,n}\).

allocate()

Create the arrays for the continuous degrees of freedom. Also set the names of the dofs.

alpha_fac

Validates that the input is among the specified integers.

Parameters:

*args – All the acceptable values for the integer

fixed_range(mmax, nmax, fixed=True)

Set the fixed property for a range of m and n values.

All modes with m <= mmax and |n| <= nmax will have their fixed property set to the value of the fixed parameter. Note that mmax and nmax are included.

Both mmax and nmax must be >= 0.

For any value of mmax, the fixed properties of R0nH, Z0nH, and rhomn are set. The fixed properties of bn are set only if mmax > 0. In other words, the bn modes are treated as having m=1.

classmethod from_RZFourier(surf, alpha_fac: int, mmax: int | None = None, nmax: int | None = None, ntheta: int | None = None, nphi: int | None = None)

Convert a SurfaceRZFourier surface to a SurfaceHenneberg surface.

Parameters:

• surf – The SurfaceRZFourier object to convert.

• mmax – Maximum poloidal mode number to include in the new surface. If None, the value mpol from the old surface will be used.

• nmax – Maximum toroidal mode number to include in the new surface. If None, the value ntor from the old surface will be used.

• ntheta – Number of grid points in the poloidal angle used for the transformation. If None, the value 3 * ntheta will be used.

• nphi – Number of grid points in the toroidal angle used for the transformation. If None, the value 3 * nphi will be used.

gamma_impl(data, quadpoints_phi, quadpoints_theta)

Evaluate the position vector on the surface in Cartesian coordinates, for a tensor product grid of points in theta and phi.

gamma_lin(data, quadpoints_phi, quadpoints_theta)

Evaluate the position vector on the surface in Cartesian coordinates, for a list of (phi, theta) points.

gammadash1_impl(data)

Evaluate the derivative of the position vector with respect to the toroidal angle phi.

gammadash2_impl(data)

Evaluate the derivative of the position vector with respect to theta.

get_dofs()

Return a 1D numpy array with all the degrees of freedom.

get_rhomn(m, n)

Return a particular \(\rho_{m,n}\) coefficient.

mmax

Validates that the input number is an integer. Optionally the user can specify the bounds for the number.

nfp

Validates that the input number is an integer. Optionally the user can specify the bounds for the number.

nmax

Validates that the input number is a positive integer number.

num_dofs()

Return the number of degrees of freedom.

set_dofs(self: simsoptpp.Surface, arg0: list[float]) → None

set_dofs_impl(v)

Set the shape coefficients from a 1D list/array

set_rhomn(m, n, val)

Set a particular \(\rho_{m,n}\) coefficient.

to_RZFourier()

Return a SurfaceRZFourier object with the identical shape. This routine implements eq (4.5)-(4.6) in the Henneberg paper, plus m=0 terms for R0 and Z0.

class simsopt.geo.SurfaceRZFourier(nfp=1, stellsym=True, mpol=1, ntor=0, quadpoints_phi=None, quadpoints_theta=None, dofs=None)

Bases: SurfaceRZFourier, Surface

SurfaceRZFourier is a surface that is represented in cylindrical coordinates using the following Fourier series:

\[r(\theta, \phi) = \sum_{m=0}^{m_{\text{pol}}} \sum_{n=-n_{\text{tor}}}^{n_\text{tor}} [ r_{c,m,n} \cos(m \theta - n_{\text{fp}} n \phi) + r_{s,m,n} \sin(m \theta - n_{\text{fp}} n \phi) ]\]

and the same for \(z(\theta, \phi)\).

Here, \((r,\phi, z)\) are standard cylindrical coordinates, and theta is any poloidal angle.

Note that for \(m=0\) we skip the \(n<0\) term for the cos terms, and the \(n \leq 0\) for the sin terms.

In addition, in the stellsym=True case, we skip the sin terms for \(r\), and the cos terms for \(z\).

For more information about the arguments quadpoints_phi, and quadpoints_theta, see the general documentation on Surfaces. Instead of supplying the quadrature point arrays along \(\phi\) and \(\theta\) directions, one could also specify the number of quadrature points for \(\phi\) and \(\theta\) using the class method from_nphi_ntheta().

Parameters:

• nfp – The number of field periods.

• stellsym – Whether the surface is stellarator-symmetric, i.e. symmetry under rotation by \(\pi\) about the x-axis.

• mpol – Maximum poloidal mode number included.

• ntor – Maximum toroidal mode number included, divided by nfp.

• quadpoints_phi – Set this to a list or 1D array to set the \(\phi_j\) grid points directly.

• quadpoints_theta – Set this to a list or 1D array to set the \(\theta_j\) grid points directly.

_make_mn()

Make the list of m and n values.

_make_names()

Form a list of names of the rc, zs, rs, or zc array elements. The order of these four arrays here must match the order in set_dofs_impl() and get_dofs() in src/simsoptpp/surfacerzfourier.h.

_validate_mn(m, n)

Check whether m and n are in the allowed range.

change_resolution(mpol, ntor)

return a new surface with Fourier resolution mpol, ntor :param mpol: new poloidal mode number :param ntor: new toroidal mode number

Returns:

A new SurfaceRZFourier object with the specified resolution.

Return type:

surf

condense_spectrum(n_theta=None, n_phi=None, power=2, maxiter=None, method='SLSQP', epsilon=0.001, Fourier_continuation=True, verbose=True)

Apply spectral condensation so the Fourier amplitudes rc and zs decay more rapidly with m and n.

This function has similarities to the spectral condensation used in VMEC, see e.g. Hirshman and Meier, Physics of Fluids 28, 1387 (1985), although the algorithm used here is different. The idea is to reparameterize the poloidal angle so that the shape can be represented with as few Fourier modes as possible.

Let \(\theta_1\) be the original poloidal angle, and let \(\theta_2\) be a new poloidal angle with \(\theta_2 = \theta_1 + \lambda(\theta_1, \phi)\) for some function \(\lambda\). This routine optimizes \(\lambda\) to minimize the spectral width while preserving the surface shape in real space as well as possible. The objective function minimized is the same function shown in the documentation for spectral_width(). At each iteration, i.e. for any specific choice of \(\lambda\), updated values of rc and zs are computed by Fourier-transforming points on the original surface, with the results used to evaluate the objective. Gradients are evaluated using JAX. Optionally, a constrained optimization method can be used to ensure that the maximum change to the surface shape is below a specified fraction of the minor radius, given by the epsilon argument.

If n_theta and n_phi are not specified, they default to 3*mpol+1 and 3*ntor+1, respectively.

The optimization algorithm is specified by method. The recommended settings are either SLSQP or trust-constr. This argument accepts any method supported by scipy.optimize.minimize (L-BFGS-B, BFGS, etc.) or any of the methods supported by scipy.optimize.least_squares (trf, dogleg, or lm). The constrained algorithms SLSQP or trust-constr are recommended to ensure that the surface shape does not change too much. You can try the unconstrained algorithms, which may be faster and result in a smaller spectral width, but not always, and sometimes the surface shape will not be preserved well. Based on experience, the least-squares methods achieve slightly lower values of the spectral width, but take more time per iteration. Any of the algorithms trf, BFGS, and lm are likely to work.

If maxiter is not specified, it defaults to 25 for the least-squares methods, and 100 for other methods.

If Fourier_continuation is True, first the Fourier modes of \(\lambda\) with \(|m|\) and \(|n|\) up to 1 are optimized, then modes with \(|m|\) and \(|n|\) up to 2, and so on. If False, all modes of \(\lambda\) are optimized at once. Using Fourier continuation generally results in better spectral condensation, but it requires more time.

The degrees of freedom of the original surface are not modified. Instead, a copy of the surface is returned with the new optimized degrees of freedom.

Spectral condensation always results in some finite change to the surface shape. Ideally this change is as small as possible. To measure this change, the return value max_RZ_error gives the maximum absolute change to R and Z on the \(\theta_1\) grid points, normalized to the minor radius. If this value is approximately equal to epsilon for the constrained algorithms, or unacceptably large for the unconstrained algorithms, you can use the change_resolution() method to increase mpol and ntor of the surface before calling condense_spectrum, so both the Fourier and grid resolutions used in condense_spectrum are increased. You may then be able to call change_resolution() again to reduce mpol and ntor after the condensation.

The data dictionary that is returned contains information about the optimization. The most noteworthy entries are max_RZ_error and spectral_width_reduction. The former indicates how well the surface shape was preserved, and the latter indicates how much the spectral width was reduced by the optimization.

The function plot_spectral_condensation() is useful for visualizing the results of spectral condensation.

Parameters:

• n_theta (int or None, optional) – Number of theta grid points to use for fitting.

• n_phi (int or None, optional) – Number of phi grid points to use for fitting.

• power (float, optional) – Power to use in the spectral width objective function.

• maxiter (int or None, optional) – Maximum number of iterations for the optimization.

• method (str or None, optional) – Optimization method to use. See above for details.

• epsilon (float, optional) – Maximum allowed change to R and Z on the θ₁ grid points, normalized to the minor radius. Only matters if method=="SLSQP" or method=="trust-constr".

• Fourier_continuation (bool, optional) – If True, increase the Fourier resolution of λ in steps (slower). If False, optimize all Fourier modes of λ at once (faster).

• verbose (bool, optional) – Whether to print progress messages.

Returns:

• surface (SurfaceRZFourier) – A new SurfaceRZFourier object with the condensed spectrum.

• data (dict) – A dictionary containing the following data from the optimization:

  • ”method” (str):

    The optimization algorithm used (value of the method argument).

  • ”maxiter” (int):

    The maximum number of iterations for the optimizer.

  • ”n_theta” (int):

    Number of theta grid points used for fitting.

  • ”n_phi” (int):

    Number of phi grid points used for fitting.

  • ”power” (float):

    The exponent used in the definition of spectral width.

  • ”epsilon” (float):

    The constraint tolerance on pointwise R and Z changes, normalized to the minor radius (only relevant for constrained optimizers).

  • ”Fourier_continuation” (bool):

    Whether Fourier continuation (progressively increasing mode count) was used.

  • ”initial_objective” (float):

    Value of the scalar objective (spectral width) before optimization.

  • ”final_objective” (float):

    Value of the scalar objective after the final optimization step.

  • ”spectral_width_reduction” (float):

    The ratio final_objective / initial_objective

  • ”max_RZ_error” (float):

    The maximum absolute change in R or Z on the θ₁ grid points, normalized by the minor radius, for the optimized poloidal angle.

  • ”m” (ndarray of int):

    The array of Fourier poloidal mode numbers used for the angle difference λ.

  • ”n” (ndarray of int):

    The array of Fourier toroidal mode numbers used for the angle difference λ.

  • ”lambda_mn” (ndarray of float):

    The optimized Fourier coefficients for the reparameterization function λ (sin coefficients only for the stellarator-symmetric case implemented here). These coefficients are in the same ordering as “m” and “n”.

  • ”x_scale” (ndarray of float):

    The exponential spectral scaling applied to λ modes.

  • ”elapsed_time” (float):

    Wall-clock time in seconds spent in the condensation routine.

  • ”minor_radius” (float):

    The minor radius used to normalize R and Z errors and the objective.

  • ”theta1_1d” (ndarray, shape (n_theta*n_phi,)):

    The flattened original poloidal angles (θ₁) corresponding to the quadrature points on the surface. Values are in radians in the interval [0, 2π).

  • ”phi_1d” (ndarray, shape (n_theta*n_phi,)):

    The flattened toroidal angles corresponding to the quadrature points on the surface. Values are in radians.

  • ”theta1_2d” (ndarray, shape (n_phi, n_theta)):

    The 2-D grid of original poloidal angles used to evaluate the surface (meshgrid of quadpoints_theta scaled by 2π).

  • ”phi_2d” (ndarray, shape (n_phi, n_theta)):

    The 2-D grid of toroidal angles used to evaluate the surface (meshgrid of quadpoints_phi scaled by 2π).

  • ”theta_optimized” (ndarray, shape (n_theta*n_phi,)):

    The optimized poloidal angles θ₂ = θ₁ + λ(θ₁, φ) evaluated at the quadrature points (flattened).

  • ”original_x” (ndarray):

    The Fourier amplitudes of the original surface prior to condensation.

  • ”condensed_x” (ndarray):

    The Fourier amplitudes of the surface after condensation.

copy(**kwargs)

Return a copy of the SurfaceRZFourier object. A range of relevant parameters of the surface can be passed to this function as keyword arguments in order to modify the properties of the returned copy. attributes changed. Keyword arguments accepted:

Kwargs:

ntheta (int): number of quadrature points in the theta direction nphi (int): number of quadrature points in the phi direction mpol (int): number of poloidal Fourier modes for the surface ntor (int): number of toroidal Fourier modes for the surface nfp (int): number of field periods stellsym (bool): whether the surface is stellarator-symmetric quadpoints_theta (NdArray[float]): theta grid points quadpoints_phi (NdArray[float]): phi grid points range (str): range of the gridpoints either ‘full torus’, ‘field period’ or ‘half period’. Ignored if quadponts are provided.

Returns:

A new SurfaceRZFourier object, with properties specified by kwargs changed.

Return type:

surf

darea()

Derivative of the area with respect to the surface Fourier coefficients. Short hand for Surface.darea_by_dcoeff()

dvolume()

Derivative of the volume with respect to the surface Fourier coefficients. Short hand for Surface.dvolume_by_dcoeff()

extend_via_normal(distance)

Extend the surface in the normal direction by a uniform distance.

NOTE this modifies the surface in place. use the surface copy method if you want to keep the original surface.

Parameters:

distance – The distance to extend the surface.

fixed_range(mmin, mmax, nmin, nmax, fixed=True)

Set the ‘fixed’ property for a range of m and n values.

All modes with m in the interval [mmin, mmax] and n in the interval [nmin, nmax] will have their fixed property set to the value of the fixed parameter. Note that mmax and nmax are included (unlike the upper bound in python’s range(min, max).)

flip_phi()

Flip the sign of the toroidal angle ϕ, i.e. mirror-reflect the surface about the x-z plane. This will reverse the sign of the rotational transform of a plasma bounded by this surface, without reversing the direction in which θ increases. This is the best way to flip the sign of the rotational transform for a vmec calculation.

flip_theta()

Flip the direction in which the poloidal angle θ increases. The physical shape of the surface in 3D will not change, only its parameterization. Note that vmec requires θ to increase as you move from the outboard to inboard side over the top of the surface. This transformation will reverse that direction.

flip_z()

Flip the sign of the z coordinate. This will flip the sign of the rotational transform of a plasma bounded by this surface. Note that vmec requires θ to increase as you move from the outboard to inboard side over the top of the surface. This z-flip transformation will reverse that direction.

fourier_transform_scalar(scalar, mpol=None, ntor=None, normalization=None, **kwargs)

Compute the Fourier components of a scalar on the surface. The scalar is evaluated at the quadrature points on the surface. The Fourier uses the conventions of the SurfaceRZFourier series, with npol going from -ntor to ntor and mpol from 0 to mpol i.e.:

\[f(\theta, \phi) = \sum_{m=0}^{mpol} \sum_{n=-npol}^{npol} A^{mn}_s \sin(m\theta - n N_{fp} \phi) + A^{mn}_c \cos(m\theta - n N_{fp} \phi)\]

Where the cosine series is only evaluated if the surface is not stellarator symmetric (if the scalar does not adhere to the symmetry of the surface, request the cosine series by setting the kwarg stellsym=False) By default, the poloidal and toroidal resolution are the same as those of the surface, but different quantities can be specified in the kwargs.

Parameters:

• scalar – 2D array of shape (numquadpoints_phi, numquadpoints_theta).

• mpol – maximum poloidal mode number of the transform, if None, the mpol attribute of the surface is used.

• ntor – maximum toroidal mode number of the transform if None, the ntor attribute of the surface is used.

• normalization – (optional) Fourier transform normalization. Can be: None: forward and back transform are not normalized. float: forward transform is divided by this number.

• stellsym – (optional) boolean to override the stellsym attribute of the surface if you want to force the calculation of the cosine series

Returns:

2-element tuple (A_mns, A_mnc), where A_mns is a 2D array of shape (mpol+1, 2*ntor+1) containing the sine coefficients, and A_mnc is a 2D array of shape (mpol+1, 2*ntor+1) containing the cosine coefficients (these are zero if the surface is stellarator symmetric).

classmethod from_focus(filename, **kwargs)

Read in a surface from a FOCUS-format file.

classmethod from_nescoil_input(filename, which_surf, **kwargs)

Read in a surface from a NESCOIL input file, as generated by the BNORM code.

Parameters:

• filename – Name of the nescin.* file to read.

• which_surf – either plasma or current, will select whether to import the plasma boundary or the current surface (i.e. winding surface) from the file

• kwargs – Any other arguments to pass to the SurfaceRZFourier constructor. You can specify quadpoints_theta and quadpoints_phi here.

classmethod from_pyQSC(stel: None, r: float = 0.1, ntheta=20, mpol=10, ntor=20, **kwargs)

Initialize the surface from a pyQSC object. This creates a surface from a near-axis equilibrium with a specified minor radius r (in meters).

Parameters:

• stel – Qsc object with a near-axis equilibrium.

• r – the near-axis coordinate radius (in meters).

• ntheta – number of points in the theta direction for the Fourier transform.

• mpol – number of poloidal Fourier modes for the surface.

• ntor – number of toroidal Fourier modes for the surface.

• kwargs – Any other arguments to pass to the SurfaceRZFourier constructor. You can specify quadpoints_theta and quadpoints_phi here.

classmethod from_vmec_input(filename: str, **kwargs)

Read in a surface from a VMEC input file. The INDATA namelist of this file will be read using f90nml. Note that this function does not require the VMEC python module.

Parameters:

• filename – Name of the input.* file to read.

• kwargs – Any other arguments to pass to the SurfaceRZFourier constructor. You can specify quadpoints_theta and quadpoints_phi here.

classmethod from_wout(filename: str, s: float = 1.0, interp_kind: str = 'linear', **kwargs)

Read in a surface from a VMEC wout output file. Note that this function does not require the VMEC python module.

Parameters:

• filename – Name of the wout_*.nc file to read.

• s – Value of normalized toroidal flux to use for the surface. The default value of 1.0 corresponds to the VMEC plasma boundary. Must lie in the interval [0, 1].

• interp_kind – Interpolation method in s. The available options correspond to the kind argument of scipy.interpolate.interp1d.

• kwargs – Any other arguments to pass to the SurfaceRZFourier constructor. You can specify quadpoints_theta and quadpoints_phi here.

get_dofs()

Return the dofs associated to this surface.

get_nml()

Generates a fortran namelist file containing the RBC/RBS/ZBC/ZBS coefficients, in the form used in VMEC and SPEC input files. The result will be returned as a string. For saving a file, see the write_nml() function.

get_rc(m, n)

Return a particular rc Parameter.

get_rs(m, n)

Return a particular rs Parameter.

get_zc(m, n)

Return a particular zc Parameter.

get_zs(m, n)

Return a particular zs Parameter.

inverse_fourier_transform_scalar(A_mns, A_mnc, normalization=None, **kwargs)

Compute the inverse Fourier transform of a scalar on the surface, specified by the Fourier coefficients. The quantity must be is evaluated at the quadrature points on the surface. The Fourier transform is defined as \(f(\theta, \phi) = \Sum_{m=0}^{mpol} \Sum_{n=-npol}^{npol} A^{mn}_s \sin(m\theta - n*Nfp*\phi) + A^{mn}_c \cos(m\theta - n*Nfp*\phi)\) Where the cosine series is only evaluated if the surface is not stellarator symmetric. Arguments:

• A_mns: 2D array of shape (mpol+1, 2*ntor+1) containing the sine coefficients

• A_mnc: 2D array of shape (mpol+1, 2*ntor+1) containing the cosine coefficients

  (these are zero if the surface is stellarator symmetric)

Optional keyword arguments:

• normalization: Fourier transform normalization. Can be:

  None: forward and back transform are not normalized float: inverse transform is multiplied by this number

• stellsym: boolean to override the stellsym attribute of the surface

make_rotating_ellipse(major_radius, minor_radius, elongation, torsion=0)

Set the surface shape to be a rotating ellipse with the given parameters.

Values of elongation larger than 1 will result in the elliptical cross-section at \(\phi=0\) being taller than it is wide. Values of elongation less than 1 will result in the elliptical cross-section at \(\phi=0\) being wider than it is tall.

The sign convention is such that both the rotating elongation and positive torsion will contribute positively to iota according to VMEC’s sign convention.

Parameters:

• major_radius – Average major radius of the surface.

• minor_radius – Average minor radius of the surface.

• elongation – Elongation of the elliptical cross-section.

• torsion – Value to use for the (m,n)=(0,1) mode of RC and -ZS, which controls the torsion of the magnetic axis.

recompute_bell(parent=None)

Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.

This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.

Need to be implemented by classes that provide a dof_setter for external handling of DOFs.

return_fn_map: Dict[str, Callable] = {'area': <instancemethod area>, 'aspect-ratio': <function Surface.aspect_ratio>, 'volume': <instancemethod volume>}

rotate_half_field_period()

Rotate the surface toroidally by half a field period.

This operation is useful when you have a surface with the bean cross-section at ϕ = π / nfp, and you want to rotate it so that the bean is at ϕ = 0.

set_dofs(self: simsoptpp.Surface, arg0: list[float]) → None

set_rc(m, n, val)

Set a particular rc Parameter.

set_rs(m, n, val)

Set a particular rs Parameter.

set_zc(m, n, val)

Set a particular zc Parameter.

set_zs(m, n, val)

Set a particular zs Parameter.

shift_theta_by_half()

Shift the origin of the poloidal angle θ by 1/2.

This operation is useful when you have a surface with θ=0 at the inboard side instead of the usual outboard side.

spectral_width(power=2)

Compute the spectral width of the surface Fourier coefficients.

For this function, the spectral width is defined as:

\[W = \frac{1}{2} \sum_{m,n} \left( \frac{ (r_{m,n}^c)^2 + (r_{m,n}^s)^2 + (z_{m,n}^c)^2 + (z_{m,n}^s)^2 }{ a^2 } \right) (m^2 + n^2)^{p}\]

where \(r_{m,n}^c\), \(r_{m,n}^s\), \(z_{m,n}^c\), and \(z_{m,n}^s\) are the Fourier coefficients of the surface shape, \(a\) is the minor radius of the surface, and \(p\) is a constant corresponding to the power argument.

This quantity is similar to the spectral width discussed in various papers related to VMEC, such as Hirshman and Meier, Physics of Fluids 28, 1387 (1985). However the definition here is not exactly the same, due to the 1/2 factor, different normalization, and inclusion of \(n\) in the weighting factor.

See also the method condense_spectrum(), which minimizes this quantity by adjusting the poloidal angle.

Parameters:

power – The power to which to raise the Fourier coefficients when computing the spectral width. Default is 2.

Returns:

The spectral width as a float.

to_RZFourier()

No conversion necessary.

write_nml(filename: str)

Writes a fortran namelist file containing the RBC/RBS/ZBC/ZBS coefficients, in the form used in VMEC and SPEC input files. To just generate the namelist as a string without saving a file, see the get_nml() function.

Parameters:

filename – Name of the file to write.

class simsopt.geo.SurfaceRZPseudospectral(mpol, ntor, nfp, r_shift=1.0, a_scale=1.0, **kwargs)

Bases: Optimizable

This class is used to replace the Fourier-space dofs of SurfaceRZFourier with real-space dofs, corresponding to the position of the surface on grid points. The advantage of having the dofs in real-space is that they are all of the same magnitude, so it is easier to know what reasonable box constraints are. This class may therefore be useful for stage-1 optimization using algorithms that require box constraints.

Presently, SurfaceRZPseudospectral assumes stellarator symmetry.

In this class, the position vector on the surface is specified on a tensor product grid of ntheta * nphi points per half period, where ntheta and nphi are both odd, phi is the standard toroidal angle, and theta is any poloidal angle. The maximum Fourier mode numbers that can be represented by this grid are mpol in the poloidal angle and ntor * nfp in the toroidal angle, where ntheta = 1 + 2 * mpol and nphi = 1 + 2 * ntor. However, due to stellarator symmetry, roughly half of the grid points are redundant. Therefore the dofs only correspond to the non-redundant points, and the remaining points are computed from the dofs using symmetry.

A SurfaceRZPseudospectral object with resolution parameters mpol and ntor has exactly the same number of dofs as a SurfaceRZFourier object with the same mpol and ntor. Specifically,

ndofs = 1 + 2 * (mpol + ntor + 2 * mpol * ntor)

This class also allows the coordinates r and z to be shifted and scaled, which may help to keep the dofs all of order 1. Letting r_dofs and z_dofs denote the dofs in this class, the actual r and z coordinates are determined via

r = r_dofs * a_scale + r_shift
z = z_dofs * a_scale

where r_shift and a_scale are optional arguments to the constructor, which would be set to roughly the major and minor radius.

Typical usage:

vmec = Vmec("input.your_filename_here")
vmec.boundary = SurfaceRZPseudospectral.from_RZFourier(vmec.boundary)

The dofs in this class are named r(jphi,jtheta) and z(jphi,jtheta), where jphi and jtheta are integer indices into the phi and theta grids.

This class does not presently implement the simsopt.geo.surface.Surface interface, e.g. there is not a gamma() function.

Parameters:

• mpol – Maximum poloidal Fourier mode number represented.

• ntor – The maximum toroidal Fourier mode number represented, divided by nfp.

• nfp – Number of field periods.

• r_shift – Constant added to the r(jphi,jtheta) dofs to get the actual major radius.

• a_scale – Dofs are multiplied by this factor to get the actual cylindrical coordinates.

_complete_grid()

Using stellarator symmetry, copy the real-space dofs to cover a full 2d (theta, phi) grid.

_make_names()

Create the list of names for the dofs.

change_resolution(mpol, ntor)

Increase or decrease the number of degrees of freedom. The new real-space dofs are obtained using Fourier interpolation. This function is useful for increasing the size of the parameter space during stage-1 optimization. If mpol and ntor are increased or unchanged, there is no loss of information. If mpol or ntor are decreased, information is lost.

Parameters:

• mpol – The new maximum poloidal mode number.

• ntor – The new maximum toroidal mode number, divided by nfp.

classmethod from_RZFourier(surff, **kwargs)

Convert a SurfaceRZFourier object to a SurfaceRZPseudospectral object.

Parameters:

• surff – The SurfaceRZFourier object to convert.

• kwargs – You can optionally provide the r_shift or a_scale arguments to the SurfaceRZPseudospectral constructor here.

to_RZFourier(**kwargs)

Convert to a SurfaceRZFourier describing the same shape.

Parameters:

kwargs – You can optionally provide the range, nphi, and/or ntheta arguments to the SurfaceRZFourier constructor here.

class simsopt.geo.SurfaceScaled(surf, scale_factors)

Bases: Optimizable

Allows you to take any Surface class and scale the dofs. This is useful for stage-1 optimization.

as_dict(serial_objs_dict) → dict

A JSON serializable dict representation of an object.

recompute_bell(parent=None)

Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.

This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.

Need to be implemented by classes that provide a dof_setter for external handling of DOFs.

to_RZFourier()

update_fixed()

Copy the fixed status from self.surf to self.

class simsopt.geo.SurfaceXYZFourier(nfp=1, stellsym=True, mpol=1, ntor=0, quadpoints_phi=None, quadpoints_theta=None, dofs=None)

Bases: SurfaceXYZFourier, Surface

SurfaceXYZFourier is a surface that is represented in Cartesian coordinates using the following Fourier series:

\[\begin{split}\hat x(\phi,\theta) &= \sum_{m=0}^{m_\text{pol}} \sum_{n=-n_{\text{tor}}}^{n_{tor}} [ x_{c,m,n} \cos(m \theta - n_\text{ fp } n \phi) + x_{s,m,n} \sin(m \theta - n_\text{ fp } n \phi)]\\ \hat y(\phi,\theta) &= \sum_{m=0}^{m_\text{pol}} \sum_{n=-n_\text{tor}}^{n_\text{tor}} [ y_{c,m,n} \cos(m \theta - n_\text{fp} n \phi) + y_{s,m,n} \sin(m \theta - n_\text{fp} n \phi)]\\ z(\phi,\theta) &= \sum_{m=0}^{m_\text{pol}} \sum_{n=-n_\text{tor}}^{n_\text{tor}} [ z_{c,m,n} \cos(m \theta - n_\text{fp}n \phi) + z_{s,m,n} \sin(m \theta - n_\text{fp}n \phi)]\end{split}\]

where

\[\begin{split}x &= \hat x \cos(\phi) - \hat y \sin(\phi)\\ y &= \hat x \sin(\phi) + \hat y \cos(\phi)\end{split}\]

Note that for \(m=0\) we skip the \(n<0\) term for the cos terms, and the \(n \leq 0\) for the sin terms.

When enforcing stellarator symmetry, we set the

\[x_{s,*,*}, ~y_{c,*,*}, \text{and} ~z_{c,*,*}\]

terms to zero.

For more information about the arguments quadpoints_phi`, and quadpoints_theta, see the general documentation on Surfaces. Instead of supplying the quadrature point arrays along \(\phi\) and \(\theta\) directions, one could also specify the number of quadrature points for \(\phi\) and \(\theta\) using the class method from_nphi_ntheta().

Parameters:

• nfp – The number of field periods.

• stellsym – Whether the surface is stellarator-symmetric, i.e. symmetry under rotation by \(\pi\) about the x-axis.

• mpol – Maximum poloidal mode number included.

• ntor – Maximum toroidal mode number included, divided by nfp.

• quadpoints_phi – Set this to a list or 1D array to set the \(\phi_j\) grid points directly.

• quadpoints_theta – Set this to a list or 1D array to set the \(\theta_j\) grid points directly.

_make_names()

Form a list of names of the xc, ys, zs, xs, yc, or zc array elements. The order of these four arrays here must match the order in set_dofs_impl() and get_dofs() in src/simsoptpp/surfacexyzfourier.h.

_make_names_helper(prefix, include0)

Helper function for _make_names method. Forms array of coefficients for :math:’m = [0, m_{pol}]’ and :math:’n = [-n_{tor}, n_{tor}]’. If :math:’m = 0’, only positive values of :math:’n’ are used. If it is a cosine term, the :math:’(0,0)’ term is included.

Parameters:

• prefix – The prefix for the name of the coefficients.

• include0 – Whether to include the (0,0) term.

extend_via_normal(distance)

Extend the surface in the normal direction by a uniform distance.

Parameters:

distance – The distance to extend the surface.

get_dofs()

Return the dofs associated to this surface.

recompute_bell(parent=None)

Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.

This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.

Need to be implemented by classes that provide a dof_setter for external handling of DOFs.

return_fn_map: Dict[str, Callable] = {'area': <instancemethod area>, 'aspect-ratio': <function Surface.aspect_ratio>, 'volume': <instancemethod volume>}

set_dofs(dofs)

Set the dofs associated to this surface.

to_RZFourier()

Return a SurfaceRZFourier instance corresponding to the shape of this surface.

class simsopt.geo.SurfaceXYZTensorFourier(nfp=1, stellsym=True, mpol=1, ntor=1, clamped_dims=[False, False, False], quadpoints_phi=None, quadpoints_theta=None, dofs=None)

Bases: SurfaceXYZTensorFourier, Surface

SurfaceXYZTensorFourier is a surface that is represented in cartesian coordinates using the following Fourier series:

\[\begin{split}\hat x(\theta, \phi) &= \sum_{i=0}^{2m_\text{pol}} \sum_{j=0}^{2n_\text{tor}} x_{ij} w_i(\theta)v_j(\phi)\\ \hat y(\theta, \phi) &= \sum_{i=0}^{2m_\text{pol}} \sum_{j=0}^{2n_\text{tor}} y_{ij} w_i(\theta)v_j(\phi)\\ x(\phi, \theta) &= \hat x(\phi, \theta) \cos(\phi) - \hat y(\phi, \theta) \sin(\phi)\\ y(\phi, \theta) &= \hat x(\phi, \theta) \sin(\phi) + \hat y(\phi, \theta) \cos(\phi)\\ z(\theta, \phi) &= \sum_{i=0}^{2m_\text{pol}} \sum_{j=0}^{2n_\text{tor}} z_{ij} w_i(\theta)v_j(\phi)\end{split}\]

where the basis functions \({v_j}\) are given by

\[\{1, \cos(1\,\mathrm{nfp}\,\phi), \ldots, \cos(n_\text{tor}\,\mathrm{nfp}\,\phi), \sin(1\,\mathrm{nfp}\,\phi), \ldots, \sin(n_\text{tor}\,\mathrm{nfp}\,\phi)\}\]

and \({w_i}\) are given by

\[\{1, \cos(1\theta), \ldots, \cos(m_\text{pol}\theta), \sin(1\theta), \ldots, \sin(m_\text{pol}\theta)\}\]

When stellsym=True the sums above change as follows:

\[\begin{split}\hat x(\theta, \phi) &= \sum_{i=0}^{m_\text{pol}} \sum_{j=0}^{n_\text{tor}} x_{ij} w_i(\theta)v_j(\phi) + \sum_{i=m_\text{pol}+1}^{2m_\text{pol}} \sum_{j=n_\text{tor}+1}^{2n_\text{tor}} x_{ij} w_i(\theta)v_j(\phi)\\ \hat y(\theta, \phi) &= \sum_{i=0}^{m_\text{pol}} \sum_{j=n_\text{tor}+1}^{2n_\text{tor}} y_{ij} w_i(\theta)v_j(\phi) + \sum_{i=m_\text{pol}+1}^{2m_\text{pol}} \sum_{j=0}^{n_\text{tor}} y_{ij} w_i(\theta)v_j(\phi)\\\\ z(\theta, \phi) &= \sum_{i=0}^{m_\text{pol}} \sum_{j=n_\text{tor}+1}^{2n_\text{tor}} z_{ij} w_i(\theta)v_j(\phi) + \sum_{i=m_\text{pol}+1}^{2m_\text{pol}} \sum_{j=0}^{n_\text{tor}} z_{ij} w_i(\theta)v_j(\phi)\end{split}\]

For more information about the arguments quadpoints_phi, and quadpoints_theta, see the general documentation on Surfaces. Instead of supplying the quadrature point arrays along \(\phi\) and \(\theta\) directions, one could also specify the number of quadrature points for \(\phi\) and \(\theta\) using the class method from_nphi_ntheta().

Parameters:

• nfp – The number of field periods.

• stellsym – Whether the surface is stellarator-symmetric, i.e. symmetry under rotation by \(\pi\) about the x-axis.

• mpol – Maximum poloidal mode number included.

• ntor – Maximum toroidal mode number included, divided by nfp.

• clamped_dims –

  ???

• quadpoints_phi – Set this to a list or 1D array to set the \(\phi_j\) grid points directly.

• quadpoints_theta – Set this to a list or 1D array to set the \(\theta_j\) grid points directly.

_make_names()

Form a list of names of the x, y, and z, array elements. The order of these four arrays here must match the order in set_dofs_impl() and get_dofs() in src/simsoptpp/surfacexyztensorfourier.h.

_make_names_helper(coord)

Helper function for _make_names method. Forms array of coefficients for :math:’i = [0, 2m_{pol}]’ and :math:’j = [0, 2n_{tor}]’. If stellarator symmetric, only some of the coefficients are used (see class docstring for more info)

Parameters:

coord – coordinate to make array for (‘x’, ‘y’, or z’).

extend_via_normal(distance)

Extend the surface in the normal direction by a uniform distance.

Parameters:

distance – The distance to extend the surface.

get_dofs()

Return the dofs associated to this surface.

get_stellsym_mask()

In the case of stellarator symmetry, some of the information is redundant, since the coordinates at (-phi, -theta) are the same (up to sign changes) to those at (phi, theta). The point of this function is to identify those angles phi and theta that we can ignore. This is difficult to do in general, so we focus on the following three common cases below:

phis = np.linspace(0, 1/self.nfp, 2*ntor+1, endpoint=False) thetas = np.linspace(0, 1, 2*mpol+1, endpoint=False)

or

phis = np.linspace(0, 1/self.nfp, 2*ntor+1, endpoint=False) thetas = np.linspace(0, 0.5, 2*mpol, endpoint=False)

or

phis = np.linspace(0, 1/(2*self.nfp), ntor+1, endpoint=False) thetas = np.linspace(0, 1, 2*mpol+1, endpoint=False)

This function could be extended to be aware of rotational symmetry as well. So far we assume that that redundancy was removed already (hence the phis only go to 1/nfp or 1/(2*nfp)).

recompute_bell(parent=None)

Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.

This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.

Need to be implemented by classes that provide a dof_setter for external handling of DOFs.

set_dofs(dofs)

Set the dofs associated to this surface.

skip(coord, i, j)

This function skips m,n values for SurfaceXYZTensorFourier that are not included in a stellarator symmetric surface. This is a Python version of skip() in surfacexyztensorfourier.h.

Parameters:

• coord – current coordinate (‘x’, ‘y’, or z’).

• i – poloidal mode index.

• j – toroidal mode index.

to_RZFourier()

Return a SurfaceRZFourier instance corresponding to the shape of this surface.

class simsopt.geo.ToroidalFlux(surface, biotsavart, idx=0, range=None, nphi=None, ntheta=None)

Bases: Optimizable

Given a surface and Biot Savart kernel, this objective calculates

\[\begin{split}J &= \int_{S_{\varphi}} \mathbf{B} \cdot \mathbf{n} ~ds, \\ &= \int_{S_{\varphi}} \text{curl} \mathbf{A} \cdot \mathbf{n} ~ds, \\ &= \int_{\partial S_{\varphi}} \mathbf{A} \cdot \mathbf{t}~dl,\end{split}\]

where \(S_{\varphi}\) is a surface of constant \(\varphi\), and \(\mathbf A\) is the magnetic vector potential.

J()

Compute the toroidal flux on the surface where \(\varphi = \texttt{quadpoints_varphi}[\texttt{idx}]\).

d2J_by_dsurfacecoefficientsdsurfacecoefficients()

Calculate the second partial derivatives with respect to the surface coefficients.

dJ(*args, partials=False, **kwargs)

dJ_by_dcoils()

Calculate the partial derivatives with respect to the coil coefficients.

dJ_by_dsurfacecoefficients()

Calculate the partial derivatives with respect to the surface coefficients.

invalidate_cache()

recompute_bell(parent=None)

Function to be called whenever new DOFs input is given or if the parent Optimizable’s data changed, so the output from the current Optimizable object is invalid.

This method gets called by various DOF setters. If only the local DOFs of an object are being set, the recompute_bell method is called in that object and also in the descendent objects that have a dependency on the object, whose local DOFs are being changed. If gloabl DOFs of an object are being set, the recompute_bell method is called in the object, ancestors of the object, as well as the descendents of the object.

Need to be implemented by classes that provide a dof_setter for external handling of DOFs.

class simsopt.geo.ToroidalWireframe(surface, n_phi, n_theta, constraint_atol=1e-10, constraint_rtol=1e-10)

Bases: object

ToroidalWireframe is a wireframe grid whose nodes are placed on a toroidal surface as a 2D grid with regular spacing in the poloidal and toroidal dimensions.

Currently only supports surfaces that exhibit stellarator symmetry.

Parameters:

• surface (SurfaceRZFourier class instance) – Toroidal surface on which on which the nodes will be placed.

• n_phi (integer) – Number of wireframe nodes per half-period in the toroidal dimension. Must be even; if an odd number is provided, it will be incremented by one.

• n_theta (integer) – Number of wireframe nodes in the poloidal dimension. Must be even; if an odd number is provided, it will be incremented by one.

• constraint_atol (double (optional)) – Absolute and relative tolerances against which constraint equations are to be evaluated (see docstring for method check_constraints for more details). Default for each is 1e-10.

• constraint_rtol (double (optional)) – Absolute and relative tolerances against which constraint equations are to be evaluated (see docstring for method check_constraints for more details). Default for each is 1e-10.

add_constraint(name, constraint_type, matrix_row, constant)

Add a linear equality constraint on the currents in the segments in the wireframe of the form

matrix_row * x = constant,

where x is the array of currents in each segment, matrix_row is a 1d array of coefficients for each segment, and constant is the constant appearing on the right-hand side

Parameters:

• name (string) – Unique name for the constraint

• constraint_type (string) – Type of constraint

• matrix_row (1d double array) – Array of coefficients as described above

• constant (double) – Constant on the right-hand side of the equation above

add_continuity_constraint(node_ind, ind_tor_in, ind_pol_in, ind_tor_out, ind_pol_out)

add_continuity_constraints()

Add constraints to ensure current continuity at each node. This is called automatically on initialization and doesn’t normally need to be called by the user.

add_poloidal_current_constraint(current)

Add constraint to require the total poloidal current through the inboard midplane to be a certain value (effectively sets the toroidal magnetic field).

Parameters:

current (double) – Total poloidal current; i.e. the sum of the currents in all poloidal segments passing through the inboard midplane. A positive poloidal current thereby creates a toroidal field in the negative toroidal direction (clockwise when viewed from above).

add_segment_constraints(segments, implicit=False)

Adds a constraint or constraints requiring the current to be zero in one or more given segments.

Parameters:

• segments (integer or array/list of integers) – Index of the segmenet or segments to be constrained

• implicit (boolean (optional)) – If true, the constraints will be marked as implicit (i.e. implied by other constraints rather than set by the user). Default is false. This option is meant primarily for internal use and should not normally be invoked by the user.

add_tfcoil_currents(n_tf, current_per_coil, constraint_atol=None, constraint_rtol=None)

Adds current to certain poloidal segments in order to form a set of planar TF coils. The loops will be spaced as evenly as possible within the wireframe.

Parameters:

• n_tf (integer) – Number of TF coils per half-period

• current_per_coil (double) – Current carried by each of the coils; positive current flows in the positive toroidal direction, thereby creating a negative toroidal magnetic field

• constraint_atol (double (optional)) – Absolute and relative tolerances to apply when checking to ensure that the added currents do not violate any existing constraints. Default is to use the value already stored within the class instance.

• constraint_rtol (double (optional)) – Absolute and relative tolerances to apply when checking to ensure that the added currents do not violate any existing constraints. Default is to use the value already stored within the class instance.

add_toroidal_current_constraint(current)

Add constraint to require the total toroidal current through a poloidal cross-section to be a certain value (effectively requires a helical current distribution when combined with a poloidal current constraint).

Parameters:

current (double) – Total toroidal current; i.e. the sum of the currents in all toroidal segments passing through a symmetry plane. A positive toroidal current thereby creates a dipole moment in the positive “z” direction.

check_constraints(currents=None, atol=None, rtol=None)

Verify that every constraint is satisfied by the present values of the segment currents. Specifically, for each constraint equation, confirm that:

\[|B * x - d| < atol + mean(|x|) * rtol,\]

where B is a vector of coefficients, x is the array of currents in each segment, d is a constant, atol is an absolute tolerance, and rtol is a relative tolerance.

Parameters:

• currents (double (optional)) – Array of segment currents to check against the constraints. If none is supplied, the internal currents array of the class instance will be checked against the constraints.

• atol (double (optional)) – Absolute and relative tolerances to apply when checking to ensure that the added currents do not violate any existing constraints. Default is to use the values already stored within the class instance.

• rtol (double (optional)) – Absolute and relative tolerances to apply when checking to ensure that the added currents do not violate any existing constraints. Default is to use the values already stored within the class instance.

Returns:

constraints_met – True if all constraints are met to within the tolerance; otherwise false

Return type:

boolean

constrain_colliding_segments(coll_func, pts_per_seg=10, **kwargs)

Constains segments found to be colliding with external objects or other spatial constraints according to a user-provided function.

Parameters:

• coll_func (function) – Function with the following interface: colliding = coll_func(x, y, z, **kwargs). Here, x, y, and z are arrays the Cartesian x, y, and z coordinates of a set of test points. The function returns a logical array (colliding) with the same dimensions of x, y, and z in which elements are True if the corresponding input point violates the spatial constraint and False otherwise.

• pts_per_seg (integer (optional)) – Number of spatial points to test along each segment for collisions. Must be at least 2; endpoints (i.e. nodes) are always included. Default is 10.

• **kwargs – Any keyword arguments to be supplied to coll_func.

constrained_segments(include='all', update=True)

Returns the IDs of the segments that are currently constrained (explicitly or implicitly) to have zero current.

Parameters:

• include (string (optional)) – ‘all’: (default) returns IDs of both explicitly and implicitly constrained segments. ‘explicit’: returns IDs only of explicitly constrained segments. ‘implicit’: returns IDs only of implicitly constrained segments.

• update (boolean (optional)) – Ensure that the implicit segments are updated before returning. Default is True. NOTE: this option is primarily for internal use and should not normally be changed by the user.

Returns:

segment_ids – IDs of the constrained segments.

Return type:

list of integers

constraint_matrices(remove_redundancies=True, remove_constrained_segments=False, assume_no_crossings=False)

Return the matrices for the system of equations that define the linear equality constraints for the wireframe segment currents. The equations have the form C * x = d, where x is a column vector with the segment currents, C is a matrix of coefficients for the segment currents in each equation, and d is a column vector of constant terms in each equation.

The matrices are initially constructed to have one row (equation) per constraint. However, if some of the constraints may be redundant. By default, this function will check the constraints for certain redundancies and remove rows from the matrix that are found to be redundant. This is necessary, e.g. for some constrained linear least-squares solvers in which the constraint matrix must be full-rank. However, the user may also opt out of row reduction such that the output matrices contain every constraint equation explicitly.

Note: the function does not put the output matrix into reduced row- echelon form or otherwise guarantee that it will be full-rank. Rather, it only checks for cases in which all four segments connected to a node are constrained to have zero current, in which case the node’s continuity constraint is redundant. If there are other redundancies, e.g. due to arbitrary constraints introduced by the user that aren’t of the usual constraint types, these may not be removed.

Parameters:

• remove_redundancies (boolean (optional)) – If true (the default option), rows in the matrix found to be redundant will be removed. If false, no checks for redundancy will be performed and all constraints will be represented in the output matrices.

• remove_constrained_segments (boolean (optional)) – If true (default is false), columns in the constraint matrix corresponding to constrained segments will be removed, and segment constraints will not appear explicitly in the constraints matrix.

• assume_no_crossings (boolean (optional)) – If true, will apply the assumption that the wireframe contains enough segment constraints such that the free segments form single-track loops with no forks or crossings. In this case, enough constraints will be removed to allow for one degree of freedom per loop.

Returns:

• constraints_C (2d double array) – The matrix B in the constraint equation

• constraints_d (1d double array (column vector)) – The column vector on the right-hand side of the constraint equation

determine_connected_segments()

Determine which segments are connected to each node.

find_inactive_nodes(assume_no_crossings=False)

Determines which nodes have no current flowing through them according to existing segment constraints (i.e. constraints that require individual segments to have zero current).

Additionally, if no crossings are assumed, identifies one node per loop of current for which the associated continuity constraint should be marked as redundant.

free_all_segments()

Remove any existing constraints that restrict individual segments to carry zero current.

get_cell_key()

Returns a matrix of the segments that border every rectangular cell in the wireframe. There is one row for every cell. The columns are defined as follows:

Column  Short name  Description
------  ----------  ----------------------------------------------------
1       ind_tor1    ID of toroidal segment with lower poloidal angle
2       ind_pol2    ID of poloidal segment with higher toroidal angle
3       ind_tor3    ID of toroidal segment with higher poloidal angle
4       ind_pol4    ID of poloidal segment with lower toroidal angle

get_cell_neighbors()

Returns a matrix of the indices of the four adjacent cells to each cell in the wireframe. There is one row for every cell. The columns are defined as follows:

Column  Short name  Description
------  ----------  ----------------------------------------------------
1       nbr_npol    ID of neighboring cell, negative poloidal direction
2       nbr_ptor    ID of neighboring cell, positive toroidal direction
3       nbr_ppol    ID of neighboring cell, positive poloidal direction
4       nbr_ntor    ID of neighboring cell, negative toroidal direction

get_free_cells(form='logical')

Returns the indices of the cells that are free; i.e. they do not border any constrained segments.

Parameters:

form (string (optional)) – If ‘indices’, will return an array giving the row indices of the free cells as they appear in the output of get_cell_key(). If ‘logical’, will return a logical array with one element per cell in which free cells are coded as true. Default is ‘logical’.

initialize_constraints()

make_plot_2d(extent='half period', quantity='currents', active_tol=1e-12, ax=None, add_colorbar=True, coordinates='indices', **kwargs)

Make a 2d plot of the segments in the wireframe grid.

Parameters:

• extent (string (optional)) – Portion of the torus to be plotted. Options are ‘half period’, ‘field period’ (default), and ‘full torus’.

• quantity (string (optional)) – Quantity to be represented in the color of each segment. Options are ‘currents’ (default), ‘nonzero currents’ (i.e. show only segments with nonzero current), and ‘constrained segments’.

• active_tol (float (optional)) – Minimum magnitude of current carried by a given segment for it to be plotted if quantity is set to ‘nonzero currents’. Default is 1e-12.

• ax (instance of the matplotlib.pyplot.Axis class (optional)) – Axis on which to generate the plot. If None, a new plot will be created.

• add_colorbar (logical (optional)) – If true, a colorbar will be added to accompany the 2d plot. Default is True.

• coordinates (string (optional)) – Coordinates to plot for the nodes of each segment. If ‘indices’, the coordinates will be the indices of the node in each dimension (toroidal and poloidal). If ‘radians’, the coordinates will be the angles in radians. If ‘degrees’, the coordinates will be the angles in degrees.

• kwargs (optional keyword arguments) – Keyword arguments to be passed to the LineCollection instance used to plot the segments

Returns:

• ax (instance of the matplotlib.pyplot.Axis class) – Axis instance on which the plot was created.

• lc (instance of the matplotlib.collections.LineCollection class) – LineCollection instance of the plotted segments.

• cb (instance of the matplotlib.colorbar class) – Colorbar instance to go with the plot on ax. If add_colorbar is False, will return None.

make_plot_3d(ax=None, engine='mayavi', extent='full torus', to_show='all', active_tol=1e-12, tube_radius=0.01, colormap=None, **kwargs)

Make a 3d plot of the wireframe grid.

Parameters:

• engine (string (optional)) – Plotting package to use, either ‘mayavi’ or ‘matplotlib’. Default is ‘mayavi’. ‘matplotlib’ is not recommended.

• extent (string (optional)) – Portion of the torus to be plotted. Options are ‘half period’, ‘field period’ (default), and ‘full torus’.

• to_show (string (optional)) – If ‘all’, will plot all segments, including those that carry no current. If ‘active’ will only show segments that carry a current of magnitude greater than active_tol. Default is ‘all’.

• active_tol (float (optional)) – Minimum magnitude of current carried by a given segment for it to be plotted if to_show is set to ‘active’. Default is 1e-12.

• tube_radius (float (optional)) – Radius of the tubes used to represent each segment in the mayavi rendering. Only used if engine is ‘mayavi’. Default is 0.01.

• colormap (2d array (optional)) – 4-column array of red, green, blue, and alpha (opacity) values on a scale of 0 to 255 applied to the coil currents

• ax (matplotlib.pyplot.axes class instance (optional)) – Axes on which to make the plot. Only used if engine is ‘matplotlib’. If not provided, a new set of axes will be generated.

• **kwargs – Additional keyword arguments to be passed to the plotting engine (mayavi only)

Returns:

ax – Axes on which the plot is generated. Only returned if engine is ‘matplotlib’.

Return type:

matplotlib.pyplot.axes class instance

plot_cells_2d(cell_values, value_label=None, ax=None)

Generates a 2d plot of the cells in one half-period of a wireframe, color-coded according to an input array of values.

Parameters:

• cell_values (1d array) – Values that determine the color of each cell.

• value_label (string (optional)) – Label for the colorbar

• ax (instance of the matplotlib.pyplot.Axis class (optional)) – Axis on which to generate the plot. If None, a new plot will be created.

Returns:

ax – Axis instance on which the plot was created.

Return type:

instance of the matplotlib.pyplot.Axis class

remove_constraint(names)

Remove a constraint from the wireframe’s set of constraints.

Parameters:

names (string or list of strings) – Name(s) of the constraints to be removed

remove_poloidal_current_constraint()

remove_segment_constraints(segments, implicit=False)

Removes constraints restricting the currents in given segment(s) to be zero.

Parameters:

• segments (integer or array/list of integers) – Index of the segmenet or segments for which constraints are to be removed

• implicit (boolean (optional)) – If true, the constraints will be marked as implicit (i.e. implied by other constraints rather than set by the user). Default is false.

remove_toroidal_current_constraint()

set_poloidal_current(current)

Set the constraint requiring the total poloidal current through the inboard midplane to be a certain value (effectively sets the toroidal magnetic field).

This method will replace an existing poloidal current constraint and create one if one does not exist.

Parameters:

current (double) – Total poloidal current; i.e. the sum of the currents in all poloidal segments passing through the inboard midplane. A positive poloidal current thereby creates a toroidal field in the negative toroidal direction (clockwise when viewed from above).

set_segments_constrained(segments, implicit=False)

Ensures that one or more given segments are constrained to have zero current.

Parameters:

• segments (integer or array/list of integers) – Index of the segmenet or segments to be constrained

• implicit (boolean (optional)) – If true, the constraints will be marked as implicit (i.e. implied by other constraints rather than set by the user). Default is false.

set_segments_free(segments)

Ensures that one or more given segments are unconstrained.

Parameters:

segments (integer or array/list of integers) – Index of the segmenet or segments to be unconstrained

set_toroidal_breaks(n_breaks, width, allow_pol_current=False)

Imposes segment constraints to prevent current from flowing toroidally through planes positioned at a given number of toroidal angles.

The spacing between the breaks will be as even as possible for the given grid dimension. For perfectly even spacing, the number of breaks must be an integer factor of the wireframe grid’s toroidal dimension (n_phi).

Parameters:

• n_breaks (integer) – Number of breaks to be inserted

• width (integer) – Toroidal extent of each break, in terms of the number of constrained toroidal segments (n_breaks*width must be less than n_phi for the wireframe)

• allow_pol_current (boolean (optional)) – If True, poloidal current will be permitted to flow in regions where toroidal current is forbidden. Otherwise, all poloidal segments within toroidal breaks will be constrained. Default is False.

set_toroidal_current(current)

Set the constraint requiring the total toroidal current through a poloidal cross-section to be a certain value (effectively requires a helical current distribution when combined with a poloidal current constraint).

This method will replace an existing toroidal current constraint and create one if one does not exist.

Parameters:

current (double) – Total toroidal current; i.e. the sum of the currents in all toroidal segments passing through a symmetry plane. A positive toroidal current thereby creates a dipole moment in the positive “z” direction.

set_up_cell_key()

Set up a matrix giving the indices of the segments forming each cell/loop in the wireframe. Also populate an array giving the indices of the four adjacent cells to each cell.

to_vtk(filename, extent='torus', extra_node_data=None, extra_segment_data=None)

Export the wireframe data to a VTK file, which can be read with ParaView.

Note: This function requires the pyevtk python package, which can be installed using pip install pyevtk.

The following datasets will be stored in the file:

current: Current [A] in each segment. Positive currents flow in the positive toroidal/poloidal direction.

constrained: 1 if the segment is constrained to carry no current; 0 otherwise

constrained_exp: 1 if the segment is explicitly constrained, i.e. set by the user to carry no current.

constrained_imp: 1 if the segment is implicitly constrained, i.e. it can carry no current due to neighboring segments being constrained.

Parameters:

• filename (string) – Name of the VTK file, without extension, to create.

• extent (string (optional)) – Portion of the torus to be represented in the file. Options are ‘half period’, ‘field period’, and ‘torus’ (default).

• extra_node_data (dictionary (optional)) – Data values to be associated with each node.

• extra_segment_data (dictionary (optional)) – Data values to be associated with each segment.

unconstrained_segments(update=True)

Returns the IDs of the segments that are unconstrained, explicitly or implicitly.

Parameters:

update (boolean (optional)) – Ensure that the implicit segments are updated before returning. Default is True. NOTE: this option is primarily for internal use and should not normally be changed by the user.

update_implicit_constraints()

Determines which segments are are implicitly constrained due to all connected segments on one or both ends being constrained.

class simsopt.geo.Volume(surface, range=None, nphi=None, ntheta=None)

Bases: Optimizable

Wrapper class for volume label.

J()

Compute the volume enclosed by the surface.

d2J_by_dsurfacecoefficientsdsurfacecoefficients()

Calculate the second derivatives with respect to the surface coefficients.

dJ(*args, partials=False, **kwargs)

dJ_by_dsurfacecoefficients()

Calculate the derivatives with respect to the surface coefficients.

class simsopt.geo.ZeroRotation(quadpoints)

Bases: Optimizable

__init__(quadpoints)

Dummy class that just returns zero for the rotation angle. Equivalent to using

rot = FrameRotation(...)
rot.fix_all()

alpha(quadpoints)

alphadash(quadpoints)

dalpha_by_dcoeff_vjp(quadpoints, v)

dalphadash_by_dcoeff_vjp(quadpoints, v)

simsopt.geo.best_nphi_over_ntheta(surf)

Given a surface, estimate the ratio of nphi / ntheta that minimizes the mesh anisotropy. This is useful for improving speed and accuracy of the virtual casing calculation. The result refers to the number of grid points in phi covering the full torus, not just one field period or half a field period. The input surface need not have range=="full torus" however; any range will work.

The result of this function will depend somewhat on the quadrature points of the input surface, but the dependence should be weak.

Parameters:

surf – A surface object.

Returns:

float with the best ratio nphi / ntheta.

simsopt.geo.boozer_surface_residual(surface, iota, G, biotsavart, derivatives=0, weight_inv_modB=False)

For a given surface, this function computes the residual

\[G\mathbf B_\text{BS}(\mathbf x) - \|\mathbf B_\text{BS}(\mathbf x)\|^2 (\mathbf x_\varphi + \iota \mathbf x_\theta)\]

as well as the derivatives of this residual with respect to surface dofs, iota, and G. In the above, \(\mathbf x\) are points on the surface, \(\iota\) is the rotational transform on that surface, and \(\mathbf B_{\text{BS}}\) is the magnetic field computed using the Biot-Savart law.

\(G\) is known for exact boozer surfaces, so if G=None is passed, then that value is used instead.

Parameters:

• surface – The surface to use for the computation

• iota – the surface rotational transform

• G – a constant that is a function of the coil currents in vacuum field

• biotsavart – the Biot-Savart magnetic field

• derivatives – how many spatial derivatives of the residual to compute

• weight_inv_modB – whether or not to weight the residual by \(1/\|\mathbf B\|\). This is useful to activate so that the residual does not scale with the coil currents.

Returns:

the residual at the surface quadrature points and optionally the spatial derivatives of the residual.

simsopt.geo.create_equally_spaced_curves(ncurves, nfp, stellsym, R0=1.0, R1=0.5, order=6, numquadpoints=None, use_jax_curve=False)

Create ncurves curves of type CurveXYZFourier of order order that will result in circular equally spaced coils (major radius R0 and minor radius R1) after applying coils_via_symmetries.

Usage example: create 4 base curves, which are then rotated 3 times and flipped for stellarator symmetry:

base_curves = create_equally_spaced_curves(4, 3, stellsym=True)
base_currents = [Current(1e5) for c in base_curves]
coils = coils_via_symmetries(base_curves, base_currents, 3, stellsym=True)

Parameters:

• ncurves – int Number of curves to create.

• nfp – int Field period symmetry of the plasma.

• stellsym – bool Whether the plasma has stellarator symmetry.

• R0 – float, optional, default=1.0 Major radius of the coils.

• R1 – float, optional, default=0.5 Minor radius of the coils.

• order – int, optional, default=6 Order of the Fourier series in the planar curve representation.

• numquadpoints – int, optional, default=None Number of quadrature points to use.

• use_jax_curve – bool, optional, default=False Whether to use JaxCurvePlanarFourier instead of CurvePlanarFourier.

Returns:

list

List of CurvePlanarFourier or JaxCurvePlanarFourier objects.

Return type:

curves

simsopt.geo.create_equally_spaced_planar_curves(ncurves, nfp, stellsym, R0=1.0, R1=0.5, order=6, numquadpoints=None, use_jax_curve=False)

Create ncurves curves of type CurvePlanarFourier of order order that will result in circular equally spaced coils (major radius R0 and minor radius R1) after applying coils_via_symmetries.

Parameters:

• ncurves – int Number of curves to create.

• nfp – int Field period symmetry of the plasma.

• stellsym – bool Whether the plasma has stellarator symmetry.

• R0 – float, optional, default=1.0 Major radius of the coils.

• R1 – float, optional, default=0.5 Minor radius of the coils.

• order – int, optional, default=6 Order of the Fourier series in the planar curve representation.

• numquadpoints – int, optional, default=None Number of quadrature points to use.

• use_jax_curve – bool, optional, default=False Whether to use JaxCurvePlanarFourier instead of CurvePlanarFourier.

Returns:

list

List of CurvePlanarFourier or JaxCurvePlanarFourier objects.

Return type:

curves

simsopt.geo.create_multifilament_grid(curve, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b, rotation_order=None, rotation_scaling=None, frame='centroid')

Create a regular grid of numfilaments_n * numfilaments_b many filaments to approximate a finite-build coil.

Note that “normal” and “binormal” in the function arguments here refer to either the Frenet frame or the “coil centroid frame” defined by Singh et al., before rotation.

Parameters:

• curve – The underlying curve.

• numfilaments_n – number of filaments in normal direction.

• numfilaments_b – number of filaments in bi-normal direction.

• gapsize_n – gap between filaments in normal direction.

• gapsize_b – gap between filaments in bi-normal direction.

• rotation_order – Fourier order (maximum mode number) to use in the expression for the rotation of the filament pack. None means that the rotation is not optimized.

• rotation_scaling – scaling for the rotation degrees of freedom. good scaling improves the convergence of first order optimization algorithms. If None, then the default of 1 / max(gapsize_n, gapsize_b) is used.

• frame – orthonormal frame to define normal and binormal before rotation (either ‘centroid’ or ‘frenet’)

simsopt.geo.create_planar_curves_between_two_toroidal_surfaces(s, s_inner, s_outer, Nx=10, Ny=10, Nz=10, order=1, use_jax_curve=False, numquadpoints=None, Nmin_factor=2.01)

Create a list of planar curves between two toroidal surfaces. The curves are initialized as circular coils of radius R and then the coils are rotated and flipped to satisfy stellarator symmetry. They are originally initialized on a uniform Cartesian grid and then filtered to only include points that are between the two toroidal surfaces.

Parameters:

• s – Surface The plasma surface object (used for nfp and geometry info).

• s_inner – Surface The inner toroidal surface (for grid bounding box). Typically generated from s.extend_via_normal().

• s_outer – Surface The outer toroidal surface (for grid bounding box). Typically generated from s.extend_via_normal() and should extend out further than s_inner.

• Nx – int, optional Number of uniform grid points in the x direction to initialize a uniform grid.

• Ny – int, optional Number of uniform grid points in the y direction to initialize a uniform grid.

• Nz – int, optional Number of uniform grid points in the z direction to initialize a uniform grid.

• order – int, optional Order of the Fourier series in the planar curve representation.

• use_jax_curve – bool, optional Whether to use JaxCurvePlanarFourier instead of CurvePlanarFourier.

• numquadpoints – int, optional Number of quadrature points to use.

• Nmin_factor – float, optional Factor to set minimum coil spacing (default: 2.01). The coil radius is set to Nmin / Nmin_factor, where Nmin is the minimum grid spacing. So as long as Nmin_factor > 2, then the coils (which are initialized as circles of radius R) will not overlap.

Returns:

list

List of CurvePlanarFourier or JaxCurvePlanarFourier objects.

all_curves : list

Return type:

curves

simsopt.geo.curves_to_vtk(curves, filename, close=False, extra_data=None)

Export a list of Curve objects in VTK format, so they can be viewed using Paraview. This function requires the python package pyevtk, which can be installed using pip install pyevtk.

Parameters:

• curves – A python list of Curve objects.

• filename – Name of the file to write.

• close – Whether to draw the segment from the last quadrature point back to the first.

simsopt.geo.fix_matplotlib_3d(ax)

Make axes of 3D plot have equal scale so that spheres appear as spheres, cubes as cubes, etc.. This is one possible solution to Matplotlib’s ax.set_aspect('equal') and ax.axis('equal') not working for 3D. This function is to be called after objects have been plotted.

This function was taken from https://stackoverflow.com/questions/13685386/matplotlib-equal-unit-length-with-equal-aspect-ratio-z-axis-is-not-equal-to

Parameters:

ax – a matplotlib axis, e.g., as output from plt.gca().

simsopt.geo.plot(items, ax=None, show=True, **kwargs)

Plot multiple Coil, Curve, and/or Surface objects together on the same axes. Any keyword arguments other than the ones listed here are passed to the plot() method of each item. A particularly useful argument is engine, which can be set to "matplotlib", "mayavi", or "plotly".

Parameters:

• items – A list of objects to plot, consisting of Coil, Curve, and Surface objects.

• ax – The axis object on which to plot. 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.

Returns:

The axis object used.

simsopt.geo.plot_spectral_condensation(surf1, surf2, data, show=True)

Plot results from condense_spectrum().

Three figures are generated.

Figure 1 shows λ as a function of θ₁ at various ϕ values (left), and the Fourier amplitudes of λ are plotted on the right.

Figure 2 shows the m- and n-dependence of rc and zs with respect to θ₁ and θ₂ in two ways: heatmaps on the left, and a scatter plot on the right.

Figure 3 shows eight cross-sections of the surface, including points that are uniformly spaced with respect to θ₁ and θ₂.

Parameters:

• surf1 (SurfaceRZFourier) – The original surface before spectral condensation.

• surf2 (SurfaceRZFourier) – The condensed surface after spectral condensation.

• data (dict) – The data dictionary returned by condense_spectrum().

• show (bool, optional) – Whether to call matplotlib’s show() function. Default is True.

Returns:

Matplotlib handles for the three figures

Return type:

fig1, fig2, fig3

simsopt.geo.signed_distance_from_surface(xyz, surface)

Compute the signed distances from points xyz to a surface. The sign is positive for points inside the volume surrounded by the surface.

simsopt.geo.windowpane_wireframe(surface, n_coils_tor, n_coils_pol, size_tor, size_pol, gap_tor, gap_pol, constraint_atol=1e-10, constraint_rtol=1e-10)

Create a ToroidalWireframe class instance with the current constrained to flow only within regularly spaced rectangular loops in the grid, i.e. windowpane coils.

NOTE: the assume_no_crossings parameter must always be set to True where relevant (e.g. for obtaining constraint matrices or for RCLS solves)

Parameters:

• surface (SurfaceRZFourier class instance) – Toroidal surface on which on which the nodes will be placed.

• n_coils_tor (integers) – Number of windowpane coils in the toroidal and poloidal dimensions

• n_coils_pol (integers) – Number of windowpane coils in the toroidal and poloidal dimensions

• size_tor (integers) – Number of wireframe grid cells per coil in the toroidal and poloidal dimensions; must be even

• size_pol (integers) – Number of wireframe grid cells per coil in the toroidal and poloidal dimensions; must be even

• gap_tor (integers) – Number of wireframe grid cells between adjacent windowpane coils in the toroidal and poloidal dimensions; must be even

• gap_pol (integers) – Number of wireframe grid cells between adjacent windowpane coils in the toroidal and poloidal dimensions; must be even

• constraint_atol (double (optional)) – Absolute and relative tolerances against which constraint equations are to be evaluated (see docstring for method check_constraints for more details). Default for each is 1e-10.

• constraint_rtol (double (optional)) – Absolute and relative tolerances against which constraint equations are to be evaluated (see docstring for method check_constraints for more details). Default for each is 1e-10.

Returns:

wframe – ToroidalWireframe instance with all segments constrained except for those that form the windowpane coils

Return type:

ToroidalWireframe class instance