simsopt.geo package

Submodules

simsopt.geo.boozersurface module

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

Bases: MSONable

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

  1. minimize_boozer_penalty_constraints_LBFGS

  2. minimize_boozer_penalty_constraints_newton

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

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

classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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.

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.JaxCurve(quadpoints, gamma_pure, **kwargs)

Bases: Curve, Curve

dgamma_by_dcoeff_impl(dgamma_by_dcoeff)

This function returns

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

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

dgamma_by_dcoeff_vjp_impl(v)

This function returns the vector Jacobian product

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

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

dgammadash_by_dcoeff_impl(dgammadash_by_dcoeff)

This function returns

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

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

dgammadash_by_dcoeff_vjp_impl(v)

This function returns

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

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

dgammadashdash_by_dcoeff_impl(dgammadashdash_by_dcoeff)

This function returns

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

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

dgammadashdash_by_dcoeff_vjp_impl(v)

This function returns the vector Jacobian product

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

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

dgammadashdashdash_by_dcoeff_impl(dgammadashdashdash_by_dcoeff)

This function returns

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

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

dgammadashdashdash_by_dcoeff_vjp_impl(v)

This function returns the vector Jacobian product

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

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

dkappa_by_dcoeff_vjp(v)

This function returns the vector Jacobian product

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

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

dtorsion_by_dcoeff_vjp(v)

This function returns the vector Jacobian product

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

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

gamma_impl(gamma, quadpoints)

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

gammadash_impl(gammadash)

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

gammadashdash_impl(gammadashdash)

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

gammadashdashdash_impl(gammadashdashdash)

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

set_dofs(self: simsoptpp.Curve, arg0: List[float]) None
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.

as_dict() dict

A JSON serializable dict representation of an object.

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.

classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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.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.curve.incremental_arclength_pure(d1gamma)

This function is used in a Python+Jax implementation of the curve arc length formula.

simsopt.geo.curve.incremental_arclength_vjp(d1gamma, v)
simsopt.geo.curve.kappa_pure(d1gamma, d2gamma)

This function is used in a Python+Jax implementation of formula for curvature.

simsopt.geo.curve.kappagrad0(d1gamma, d2gamma)
simsopt.geo.curve.kappagrad1(d1gamma, d2gamma)
simsopt.geo.curve.kappavjp0(d1gamma, d2gamma, v)
simsopt.geo.curve.kappavjp1(d1gamma, d2gamma, v)
simsopt.geo.curve.torsion_pure(d1gamma, d2gamma, d3gamma)

This function is used in a Python+Jax implementation of formula for torsion.

simsopt.geo.curve.torsionvjp0(d1gamma, d2gamma, d3gamma, v)
simsopt.geo.curve.torsionvjp1(d1gamma, d2gamma, d3gamma, v)
simsopt.geo.curve.torsionvjp2(d1gamma, d2gamma, d3gamma, v)

simsopt.geo.curvehelical module

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

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

as_dict() dict

A JSON serializable dict representation of an object.

classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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.curvehelical.jaxHelicalfouriercurve_pure(dofs, quadpoints, order, n0, l0, R0, r0)

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.

as_dict() dict

A JSON serializable dict representation of an object.

dJ(*args, partials=False, **kwargs)
classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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.

as_dict() dict

A JSON serializable dict representation of an object.

compute_candidates()
dJ(*args, partials=False, **kwargs)
classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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.

as_dict() dict

A JSON serializable dict representation of an object.

dJ(*args, partials=False, **kwargs)
classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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.

as_dict() dict

A JSON serializable dict representation of an object.

compute_candidates()
dJ(*args, partials=False, **kwargs)
classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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

as_dict() dict

A JSON serializable dict representation of an object.

dJ(*args, partials=False, **kwargs)
classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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.

as_dict() dict

A JSON serializable dict representation of an object.

dJ(*args, partials=False, **kwargs)
classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

return_fn_map: Dict[str, Callable] = {'J': <function LpCurveTorsion.J>, 'dJ': <function derivative_dec.<locals>._derivative_dec>}
simsopt.geo.curveobjectives.Lp_curvature_pure(kappa, gammadash, p, desired_kappa)

This function is used in a Python+Jax implementation of the curvature penalty term.

simsopt.geo.curveobjectives.Lp_torsion_pure(torsion, gammadash, p, threshold)

This function is used in a Python+Jax implementation of the formula for the torsion penalty term.

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.

as_dict() dict

A JSON serializable dict representation of an object.

dJ(*args, partials=False, **kwargs)
classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

class simsopt.geo.curveobjectives.MinimumDistance(*args, **kwargs)

Bases: CurveCurveDistance

simsopt.geo.curveobjectives.cc_distance_pure(gamma1, l1, gamma2, l2, minimum_distance)

This function is used in a Python+Jax implementation of the curve-curve distance formula.

simsopt.geo.curveobjectives.cs_distance_pure(gammac, lc, gammas, ns, minimum_distance)

This function is used in a Python+Jax implementation of the curve-surface distance formula.

simsopt.geo.curveobjectives.curve_arclengthvariation_pure(l, mat)

This function is used in a Python+Jax implementation of the curve arclength variation.

simsopt.geo.curveobjectives.curve_length_pure(l)

This function is used in a Python+Jax implementation of the curve length formula.

simsopt.geo.curveobjectives.curve_msc_pure(kappa, gammadash)

This function is used in a Python+Jax implementation of the mean squared curvature objective.

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)
as_dict() dict

A JSON serializable dict representation of an object.

dgamma_by_dcoeff_vjp(v)
dgammadash_by_dcoeff_vjp(v)
dgammadashdash_by_dcoeff_vjp(v)
dgammadashdashdash_by_dcoeff_vjp(v)
classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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, sigma, length_scale, n_derivs=1)

Bases: MSONable

__init__(points, sigma, length_scale, n_derivs=1)

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.

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

Bases: MSONable

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)

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}} x_{c,m}\cos(n_{\text{fp}} m \phi) + \sum_{m=1}^{\text{order}} x_{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}]\]
as_dict() dict

A JSON serializable dict representation of an object.

classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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)

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]\]
as_dict() dict

A JSON serializable dict representation of an object.

classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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, …

set_dofs(dofs)

This function sets the dofs associated to this object.

class simsopt.geo.curvexyzfourier.JaxCurveXYZFourier(quadpoints, order)

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.

as_dict() dict

A JSON serializable dict representation of an object.

classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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.curvexyzfourier.jaxfouriercurve_pure(dofs, quadpoints, order)

simsopt.geo.finitebuild module

class simsopt.geo.finitebuild.CurveFilament(curve, dn, db, rotation=None)

Bases: Curve, Curve

__init__(curve, dn, db, rotation=None)

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 normal and binormal vector and then creates a grid of curves by shifting along the normal and binormal vector. In addition, we specify an angle along the curve that allows us to optimise for the rotation of the winding pack.

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

Note that “normal” and “binormal” in the function arguments here refer not to the Frenet frame but rather to the “coil centroid frame” defined by Singh et al., before rotation.

Parameters
  • curve – the underlying curve

  • dn – how far to move in normal direction

  • db – how far to move in binormal direction

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

dgamma_by_dcoeff_vjp(v)
dgammadash_by_dcoeff_vjp(v)
gamma_impl(self: simsoptpp.Curve, arg0: numpy.ndarray[numpy.float64], arg1: numpy.ndarray[numpy.float64]) None
gammadash_impl(gammadash)
recompute_bell(parent=None)

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

class simsopt.geo.finitebuild.FilamentRotation(quadpoints, order, scale=1.0)

Bases: Optimizable

__init__(quadpoints, order, scale=1.0)

The rotation of the 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.finitebuild.ZeroRotation(quadpoints)

Bases: Optimizable

__init__(quadpoints)

Dummy class that just returns zero for the rotation angle. Equivalent to using

rot = FilamentRotation(...)
rot.fix_all()
alpha(quadpoints)
alphadash(quadpoints)
dalpha_by_dcoeff_vjp(quadpoints, v)
dalphadash_by_dcoeff_vjp(quadpoints, v)
simsopt.geo.finitebuild.create_multifilament_grid(curve, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b, rotation_order=None, rotation_scaling=None)

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 not to the Frenet frame but rather to 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.

simsopt.geo.finitebuild.jaxrotation_pure(dofs, points, order)
simsopt.geo.finitebuild.jaxrotationdash_pure(dofs, points, order)
simsopt.geo.finitebuild.rotated_centroid_frame(gamma, gammadash, alpha)
simsopt.geo.finitebuild.rotated_centroid_frame_dash(gamma, gammadash, gammadashdash, alpha, alphadash)
simsopt.geo.finitebuild.rotated_centroid_frame_dash_dcoeff_vjp0(gamma, gammadash, gammadashdash, alpha, alphadash, v)
simsopt.geo.finitebuild.rotated_centroid_frame_dash_dcoeff_vjp1(gamma, gammadash, gammadashdash, alpha, alphadash, v)
simsopt.geo.finitebuild.rotated_centroid_frame_dash_dcoeff_vjp2(gamma, gammadash, gammadashdash, alpha, alphadash, v)
simsopt.geo.finitebuild.rotated_centroid_frame_dash_dcoeff_vjp3(gamma, gammadash, gammadashdash, alpha, alphadash, v)
simsopt.geo.finitebuild.rotated_centroid_frame_dash_dcoeff_vjp4(gamma, gammadash, gammadashdash, alpha, alphadash, v)
simsopt.geo.finitebuild.rotated_centroid_frame_dcoeff_vjp0(gamma, gammadash, alpha, v)
simsopt.geo.finitebuild.rotated_centroid_frame_dcoeff_vjp1(gamma, gammadash, alpha, v)
simsopt.geo.finitebuild.rotated_centroid_frame_dcoeff_vjp2(gamma, gammadash, alpha, v)
simsopt.geo.finitebuild.rotation_dcoeff(points, order)
simsopt.geo.finitebuild.rotationdash_dcoeff(points, order)

simsopt.geo.jit module

simsopt.geo.jit.jit(fun)

simsopt.geo.plot module

simsopt.geo.plot.fix_matplotlib_3d(ax)

Make axes of 3D plot have equal scale so that spheres appear as spheres, cubes as cubes, etc.. This is one possible solution to Matplotlib’s ax.set_aspect('equal') and ax.axis('equal') not working for 3D. This function is to be called after objects have been plotted.

This function was taken from https://stackoverflow.com/questions/13685386/matplotlib-equal-unit-length-with-equal-aspect-ratio-z-axis-is-not-equal-to

Parameters

ax – a matplotlib axis, e.g., as output from plt.gca().

simsopt.geo.plot.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: MSONable

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

as_dict() dict

A JSON serializable dict representation of an object.

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. The main difficult part of this calculation is the mean cross sectional area. This 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.

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.

get_quadpoints(quadpoints_theta=None, range='full torus', nphi=None, ntheta=None, nfp=1)

This function is used to set the theta and phi grid points for Surface subclasses. It is typically called in the constructor of each Surface subclass.

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. If quadpoints_phi is specified, range is irrelevant.

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

Returns

Tuple containing

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

  • quadpoints_theta: 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

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() dict

A JSON serializable dict representation of an object.

classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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, nphi=None, ntheta=None, range='full torus', quadpoints_phi=None, quadpoints_theta=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 nphi, ntheta, range, quadpoints_phi, and quadpoints_theta, see the general documentation on Surfaces.

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.

  • 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 SurfaceGarabedian.RANGE_FULL_TORUS) to generate points up to 1 (with no point at 1). Set to "field period" (or equivalently SurfaceGarabedian.RANGE_FIELD_PERIOD) to generate points up to \(1/n_{fp}\) (with no point at \(1/n_{fp}\)). Set to "half period" (or equivalently SurfaceGarabedian.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. If quadpoints_phi is specified, range is irrelevant.

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

as_dict() dict

A JSON serializable dict representation of an object.

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

classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

get_Delta(m, n)

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

get_dofs()

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

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, nphi: Optional[int] = None, ntheta: Optional[int] = None, range: str = 'full torus', quadpoints_phi: Optional[Union[Sequence[Real], NDArray[Any, float64]]] = None, quadpoints_theta: Optional[Union[Sequence[Real], NDArray[Any, float64]]] = 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 nphi, ntheta, range, quadpoints_phi, and quadpoints_theta, see the general documentation on Surfaces.

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.

  • 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 SurfaceHenneberg.RANGE_FULL_TORUS) to generate points up to 1 (with no point at 1). Set to "field period" (or equivalently SurfaceHenneberg.RANGE_FIELD_PERIOD) to generate points up to \(1/n_{fp}\) (with no point at \(1/n_{fp}\)). Set to "half period" (or equivalently SurfaceHenneberg.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. If quadpoints_phi is specified, range is irrelevant.

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

as_dict() dict

A JSON serializable dict representation of an object.

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: Optional[int] = None, nmax: Optional[int] = None, ntheta: Optional[int] = None, nphi: Optional[int] = 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.

classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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.

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)

Bases: Optimizable

Wrapper class for surface area label.

J()

Compute the area of a surface.

d2J_by_dsurfacecoefficientsdsurfacecoefficients()

Calculate the second derivatives with respect to the surface coefficients.

dJ_by_dsurfacecoefficients()

Calculate the derivatives with respect to the surface coefficients.

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

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 derivatives with respect to the surface coefficients.

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.Volume(surface)

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_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_\text{BS}(\mathbf x) - ||\mathbf B_\text{BS}(\mathbf x)||^2 (\mathbf x_\varphi + \iota \mathbf x_\theta)\]

as well as the derivatives of this residual with respect to surface dofs, iota, and G. In the above, \(\mathbf x\) are points on the surface, \(\iota\) is the rotational transform on that surface, and \(\mathbf B_{\text{BS}}\) is the magnetic field computed using the Biot-Savart law.

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

simsopt.geo.surfaceobjectives.parameter_derivatives(surface: Surface, shape_gradient: NDArray[Any, float64]) NDArray[Any, float64]

Converts the shape gradient of a given figure of merit, \(f\), to derivatives with respect to parameters defining a surface. For a perturbation to the surface \(\delta \vec{x}\), the resulting perturbation to the objective function is

\[\delta f(\delta \vec{x}) = \int d^2 x \, G \delta \vec{x} \cdot \vec{n}\]

where \(G\) is the shape gradient and \(\vec{n}\) is the unit normal. Given \(G\), the parameter derivatives are then computed as

\[\frac{\partial f}{\partial \Omega} = \int d^2 x \, G \frac{\partial\vec{x}}{\partial \Omega} \cdot \vec{n},\]

where \(\Omega\) is any parameter of the surface.

Parameters
  • surface – The surface to use for the computation

  • shape_gradient – 2d array of size (numquadpoints_phi,numquadpoints_theta)

Returns

1d array of size (ndofs)

simsopt.geo.surfacerzfourier module

class simsopt.geo.surfacerzfourier.SurfaceRZFourier(nfp=1, stellsym=True, mpol=1, ntor=0, nphi=None, ntheta=None, range='full torus', quadpoints_phi=None, quadpoints_theta=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 nphi, ntheta, range, quadpoints_phi, and quadpoints_theta, see the general documentation on Surfaces.

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.

  • 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 SurfaceRZFourier.RANGE_FULL_TORUS) to generate 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. If quadpoints_phi is specified, range is irrelevant.

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

as_dict() dict

A JSON serializable dict representation of an object.

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_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

classmethod from_focus(filename, **kwargs)

Read in a surface from a FOCUS-format file.

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

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.

as_dict() dict

A JSON serializable dict representation of an object.

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.

classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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, nphi=None, ntheta=None, range='full torus', quadpoints_phi=None, quadpoints_theta=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 nphi, ntheta, range, quadpoints_phi, and quadpoints_theta, see the general documentation on Surfaces.

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.

  • 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 SurfaceXYZFourier.RANGE_FULL_TORUS) to generate points up to 1 (with no point at 1). Set to "field period" (or equivalently SurfaceXYZFourier.RANGE_FIELD_PERIOD) to generate points up to \(1/n_{fp}\) (with no point at \(1/n_{fp}\)). Set to "half period" (or equivalently SurfaceXYZFourier.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. If quadpoints_phi is specified, range is irrelevant.

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

as_dict() dict

A JSON serializable dict representation of an object.

classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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], nphi=None, ntheta=None, range='full torus', quadpoints_phi=None, quadpoints_theta=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 nphi, ntheta, range, quadpoints_phi, and quadpoints_theta, see the general documentation on Surfaces.

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

    ???

  • 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 SurfaceXYZTensorFourier.RANGE_FULL_TORUS) to generate points up to 1 (with no point at 1). Set to "field period" (or equivalently SurfaceXYZTensorFourier.RANGE_FIELD_PERIOD) to generate points up to \(1/n_{fp}\) (with no point at \(1/n_{fp}\)). Set to "half period" (or equivalently SurfaceXYZTensorFourier.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. If quadpoints_phi is specified, range is irrelevant.

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

as_dict() dict

A JSON serializable dict representation of an object.

classmethod from_dict(d)
Parameters

d – Dict representation.

Returns

MSONable class.

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.

Module contents