simsopt.geo package

Submodules

simsopt.geo.boozersurface module

class simsopt.geo.boozersurface.BoozerSurface(biotsavart, surface, label, targetlabel)

Bases: Optimizable

BoozerSurface and its associated methods can be used to compute the Boozer angles on a surface. It takes a Surface representation (e.g. SurfaceXYZFourier, or SurfaceXYZTensorFourier), a magnetic field evaluator, surface label evaluator, and a target surface label.

The Boozer angles are computed by solving a constrained least squares problem. The least squares objective is given by \(J(x) = \frac{1}{2} \mathbf r^T(x) \mathbf r(x)\), where \(\mathbf r\) is a vector of residuals computed by boozer_surface_residual (see surfaceobjectives.py), and some constraints. This objective is zero when the surface corresponds to a magnetic surface of the field and \((\phi,\theta)\) that parametrize the surface correspond to Boozer angles, and the constraints are satisfied.

The surface label can be area, volume, or toroidal flux. The label on the computed surface will be equal or close to the user-provided targetlabel, depending on how the label constraint is imposed. This constrained least squares problem can be solved by scalarizing and adding the constraint as an additional penalty term to the objective. This is done in

minimize_boozer_penalty_constraints_LBFGS

minimize_boozer_penalty_constraints_newton

minimize_boozer_penalty_constraints_ls

where LBFGS, Newton, or scipy.optimize.least_squares optimizers are used, respectively. Alternatively, the exactly constrained least squares optimization problem can be solved. This is done in

minimize_boozer_exact_constraints_newton

where Newton is used to solve the first order necessary conditions for optimality.

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}\]

where \(l\) is the surface label and \(l_0\) is the target surface label, \(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\). We can also optionally return the first derivatives of these optimality conditions.

boozer_penalty_constraints(x, derivatives=0, constraint_weight=1.0, scalarize=True, optimize_G=False)

Define the residual

\[\mathbf r(x) = [r_1(x),...,r_n(x), \sqrt{w_c} (l-l_0), \sqrt{w_c} (z(\varphi=0, \theta=0) - 0)]\]

where \(w_c\) is the constraint weight, \(r_i\) are the Boozer residuals at quadrature points \(1,\dots,n\), \(l\) is the surface label, and \(l_0\) is the target surface label.

For scalarized=False, this function returns \(\mathbf r(x)\) and optionally the Jacobian of \(\mathbf r(x)\).

for scalarized=True, this function returns

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

i.e. the least squares residual and optionally the gradient and the Hessian of \(J(x)\).

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 Lagrange multipliers and 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.

minimize_boozer_penalty_constraints_LBFGS(tol=0.001, maxiter=1000, constraint_weight=1.0, iota=0.0, G=None)

This function tries 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\). This is done using LBFGS.

minimize_boozer_penalty_constraints_ls(tol=1e-12, maxiter=10, constraint_weight=1.0, iota=0.0, G=None, method='lm'): 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.

minimize_boozer_penalty_constraints_newton(tol=1e-12, maxiter=10, constraint_weight=1.0, iota=0.0, G=None, stab=0.0): This function does the same as minimize_boozer_penalty_constraints_LBFGS, but instead of LBFGS it uses Newton’s method.

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.

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

This function solves the Boozer Surface residual equation exactly. For this to work, we need the right balance of quadrature points, degrees of freedom and constraints. For this reason, this only works 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 stellsym 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.

simsopt.geo.config module

simsopt.geo.curve module

class simsopt.geo.curve.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.

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

simsopt.geo.curve.create_equally_spaced_curves(ncurves, nfp, stellsym, R0=1.0, R1=0.5, order=6, numquadpoints=None)

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)

simsopt.geo.curve.create_equally_spaced_planar_curves(ncurves, nfp, stellsym, R0=1.0, R1=0.5, order=6, numquadpoints=None): 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.

simsopt.geo.curve.curves_to_vtk(curves, filename, close=False)

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.curvehelical module

class simsopt.geo.curvehelical.CurveHelical(quadpoints, order, n0, l0, R0, r0, **kwargs)

Bases: JaxCurve

Curve representation of a helical coil. The helical coil positions are specified by a poloidal angle eta that is a function of the toroidal angle phi with Fourier coefficients \(A_k\) and \(B_k\). The poloidal angle is represented as

\[\eta = m_0 \phi/l_0 + \sum_k A_k \cos(n_0 \phi k/l_0) + B_k \sin(n_0 \phi k/l_0)\]

Parameters:

quadpoints – number of grid points/resolution along the curve;
order – number of fourier coefficients, i.e., the length of the array A (or B);
n0 – toroidal periodicity/number of field periods
l0 – number of \(2\pi\) turns in \(\phi\)
R0 – major radius
r0 – minor radius

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

num_dofs(self: simsoptpp.Curve) → int

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

simsopt.geo.curveobjectives module

class simsopt.geo.curveobjectives.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.curveobjectives.CurveCurveDistance(curves, minimum_distance, num_basecurves=None)

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.curveobjectives.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.curveobjectives.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.curveobjectives.LinkingNumber(curves)

Bases: Optimizable

J()

curves

Compute the Linking number of a set of curves (whether the curves are interlocked or not).

The value is 1 if the are interlocked, 0 if not.

\[Link(c1,c2) = \frac{1}{4\pi} \oint_{c1}\oint_{c2}\frac{\textbf{R1} - \textbf{R2}}{|\textbf{R1}-\textbf{R2}|^3} (d\textbf{R1} \times d\textbf{R2})\]

where \(c1\) is the first curve and \(c2\) is the second curve, \(\textbf{R1}\) is the radius vector of the first curve, and \(\textbf{R2}\) is the radius vector of the second curve

Parameters:: curves – the set of curves on which the linking number should be computed.

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

class simsopt.geo.curveobjectives.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.curveobjectives.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.curveobjectives.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)

simsopt.geo.curveperturbed module

class simsopt.geo.curveperturbed.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 randomgen import SeedSequence, PCG64
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 = np.random.Generator(PCG64(seeds_sys[j], inc=0))
    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.curveperturbed.GaussianSampler(points: Sequence[Real] | NDArray[Any, 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[Any, float64]

sigma: float

class simsopt.geo.curveperturbed.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()

simsopt.geo.curverzfourier module

class simsopt.geo.curverzfourier.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.

simsopt.geo.curveplanarfourier module

class simsopt.geo.curveplanarfourier.CurvePlanarFourier(quadpoints, order, nfp, stellsym, 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. The resulting planar curve is then rotated in three dimensions using a quaternion, and finally a translation is applied. 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} &= [q_0,q_i,q_j,q_k]\\&= [\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}}, q_0, q_i, q_j, q_k, x_{\text{center}}, y_{\text{center}}, z_{\text{center}}]\]

get_dofs(): This function returns the dofs associated to this object.

set_dofs(dofs): This function sets the dofs associated to this object.

simsopt.geo.curvexyzfourier module

class simsopt.geo.curvexyzfourier.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]\]

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)

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.

Returns:

A list of CurveXYZFourier objects.

set_dofs(dofs): This function sets the dofs associated to this object.

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

Bases: JaxCurve

A Python+Jax implementation of the CurveXYZFourier class. There is actually no reason why one should use this over the C++ implementation in simsoptpp, but 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.

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.

simsopt.geo.finitebuild module

class simsopt.geo.finitebuild.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.

simsopt.geo.finitebuild.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.framedcurve module

class simsopt.geo.framedcurve.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.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.framedcurve.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.framedcurve.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.framedcurve.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.jit module

simsopt.geo.jit.jit(fun)

simsopt.geo.permanent_magnet_grid module

class simsopt.geo.permanent_magnet_grid.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.

simsopt.geo.plotting module

simsopt.geo.plotting.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.plotting.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.qfmsurface module

class simsopt.geo.qfmsurface.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.

simsopt.geo.strain_optimization module

class simsopt.geo.strain_optimization.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.strain_optimization.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.strain_optimization.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>}

simsopt.geo.surface module

class simsopt.geo.surface.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 poloidal angle based on arclenth along magnetic surface at constant phi. The resulting angle is in the range [0,1]. This is required for evaluating the adjoint shape gradient for free-boundary calculations.

Returns:: 2d array of shape (numquadpoints_phi, numquadpoints_theta) containing the arclength 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)

This function takes in a cylindrical angle \(\phi\) and returns the cross section of the surface in that plane evaluated at thetas. This is done using the method of bisection. This function takes in a cylindrical angle \(\phi\) and returns the cross section of the surface in that plane evaluated at thetas. This is done using the method of bisection.

This function assumes that the surface intersection with the plane is a single curve.

daspect_ratio_by_dcoeff(): Return the derivative of the aspect ratio with respect to the surface coefficients

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.

get_phi_quadpoints(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\).

get_quadpoints(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\).

get_theta_quadpoints()

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.

Returns:

2d array (numquadpoints_phi,numquadpoints_theta): defining interpolated function on arclength angle along curve at constant phi

Return type:

function_interpolated

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.

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 – Name of the file to write
extra_data – An optional data field on the surface, which can be associated with a colormap in Paraview.

class simsopt.geo.surface.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.surface.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.

simsopt.geo.surface.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.surface.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.surfacegarabedian module

class simsopt.geo.surfacegarabedian.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.

simsopt.geo.surfacehenneberg module

class simsopt.geo.surfacehenneberg.SurfaceHenneberg(nfp: int = 1, alpha_fac: int = 1, mmax: int = 1, nmax: int = 0, quadpoints_phi: Sequence[Real] | NDArray[Any, float64] | None = None, quadpoints_theta: Sequence[Real] | NDArray[Any, float64] | None = None, dofs: DOFs | None = 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.

simsopt.geo.surfaceobjectives module

class simsopt.geo.surfaceobjectives.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.surfaceobjectives.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.surfaceobjectives.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.surfaceobjectives.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.surfaceobjectives.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.surfaceobjectives.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.surfaceobjectives.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.surfaceobjectives.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.

simsopt.geo.surfaceobjectives.boozer_surface_residual(surface, iota, G, biotsavart, derivatives=0)

For a given surface, this function computes the residual

\[G\mathbf B(\mathbf x) - \|\mathbf B(\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\) is the magnetic field.

\(G\) is known for exact boozer surfaces, so if G=None is passed, then that value is used instead.

simsopt.geo.surfacerzfourier module

class simsopt.geo.surfacerzfourier.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): Change the values of mpol and ntor. Any new Fourier amplitudes will have a magnitude of zero. Any previous nonzero Fourier amplitudes that are not within the new range will be discarded.

darea(): Short hand for Surface.darea_by_dcoeff()

dvolume(): Short hand for Surface.dvolume_by_dcoeff()

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

classmethod from_focus(filename, **kwargs): Read in a surface from a FOCUS-format file.

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.

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

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

simsopt.geo.surfacexyzfourier module

class simsopt.geo.surfacexyzfourier.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.

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.

simsopt.geo.surfacexyztensorfourier module

class simsopt.geo.surfacexyztensorfourier.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.

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.

to_RZFourier(): Return a SurfaceRZFourier instance corresponding to the shape of this surface.