simsopt.field package

Submodules

simsopt.field.biotsavart module

class simsopt.field.biotsavart.BiotSavart(coils)

Bases: BiotSavart, MagneticField

Computes the MagneticField induced by a list of closed curves \(\Gamma_k\) with electric currents \(I_k\). The field is given by

\[B(\mathbf{x}) = \frac{\mu_0}{4\pi} \sum_{k=1}^{n_\mathrm{coils}} I_k \int_0^1 \frac{(\Gamma_k(\phi)-\mathbf{x})\times \Gamma_k'(\phi)}{\|\Gamma_k(\phi)-\mathbf{x}\|^3} d\phi\]

where \(\mu_0=4\pi 10^{-7}\) is the magnetic constant.

Parameters:

coils – A list of simsopt.field.coil.Coil objects.

A_and_dA_vjp(v, vgrad)

Same as simsopt.geo.biotsavart.BiotSavart.A_vjp but returns the vector Jacobian product for \(A\) and \(\nabla A\), i.e. it returns

\[\{ \sum_{i=1}^{n} \mathbf{v}_i \cdot \partial_{\mathbf{c}_k} \mathbf{A}_i \}_k, \{ \sum_{i=1}^{n} {\mathbf{v}_\mathrm{grad}}_i \cdot \partial_{\mathbf{c}_k} \nabla \mathbf{A}_i \}_k.\]

A_vjp(v)

Assume the field was evaluated at points \(\mathbf{x}_i, i\in \{1, \ldots, n\}\) and denote the value of the field at those points by \(\{\mathbf{A}_i\}_{i=1}^n\). These values depend on the shape of the coils, i.e. on the dofs \(\mathbf{c}_k\) of each coil. This function returns the vector Jacobian product of this dependency, i.e.

\[\{ \sum_{i=1}^{n} \mathbf{v}_i \cdot \partial_{\mathbf{c}_k} \mathbf{A}_i \}_k.\]

B_and_dB_vjp(v, vgrad)

Same as simsopt.geo.biotsavart.BiotSavart.B_vjp but returns the vector Jacobian product for \(B\) and \(\nabla B\), i.e. it returns

\[\{ \sum_{i=1}^{n} \mathbf{v}_i \cdot \partial_{\mathbf{c}_k} \mathbf{B}_i \}_k, \{ \sum_{i=1}^{n} {\mathbf{v}_\mathrm{grad}}_i \cdot \partial_{\mathbf{c}_k} \nabla \mathbf{B}_i \}_k.\]

B_vjp(v)

Assume the field was evaluated at points \(\mathbf{x}_i, i\in \{1, \ldots, n\}\) and denote the value of the field at those points by \(\{\mathbf{B}_i\}_{i=1}^n\). These values depend on the shape of the coils, i.e. on the dofs \(\mathbf{c}_k\) of each coil. This function returns the vector Jacobian product of this dependency, i.e.

\[\{ \sum_{i=1}^{n} \mathbf{v}_i \cdot \partial_{\mathbf{c}_k} \mathbf{B}_i \}_k.\]

as_dict(serial_objs_dict) → dict

A JSON serializable dict representation of an object.

d2A_by_dXdcoilcurrents(compute_derivatives=1)

d2B_by_dXdcoilcurrents(compute_derivatives=1)

d3A_by_dXdXdcoilcurrents(compute_derivatives=2)

d3B_by_dXdXdcoilcurrents(compute_derivatives=2)

dA_by_dcoilcurrents(compute_derivatives=0)

dB_by_dcoilcurrents(compute_derivatives=0)

classmethod from_dict(d, serial_objs_dict, recon_objs)

Parameters:

d – Dict representation.

Returns:

GSONable class.

simsopt.field.boozermagneticfield module

class simsopt.field.boozermagneticfield.BoozerAnalytic(etabar, B0, N, G0, psi0, iota0, Bbar=1.0, I0=0.0, G1=0.0, I1=0.0, K1=0.0)

Bases: BoozerMagneticField

Computes a BoozerMagneticField based on a first-order expansion in distance from the magnetic axis (Landreman & Sengupta, Journal of Plasma Physics 2018). Here the magnetic field strength is expressed as,

\[B(s,\theta,\zeta) = B_0 \left(1 + \overline{\eta} \sqrt{2s\psi_0/\overline{B}}\cos(\theta - N \zeta)\right),\]

the covariant components are,

\[ \begin{align}\begin{aligned}G(s) = G_0 + \sqrt{2s\psi_0/\overline{B}} G_1\\I(s) = I_0 + \sqrt{2s\psi_0/\overline{B}} I_1\\K(s,\theta,\zeta) = \sqrt{2s\psi_0/\overline{B}} K_1 \sin(\theta - N \zeta),\end{aligned}\end{align} \]

and the rotational transform is,

\[\iota(s) = \iota_0.\]

While formally \(I_0 = I_1 = G_1 = K_1 = 0\), these terms have been included in order to test the guiding center equations at finite beta.

Parameters:

• etabar – magnitude of first order correction to magnetic field strength

• B0 – magnetic field strength on the axis

• N – helicity of symmetry (integer)

• G0 – lowest order toroidal covariant component

• psi0 – (toroidal flux)/ (2*pi) on the boundary

• iota0 – lowest order rotational transform

• Bbar – normalizing magnetic field strength (defaults to 1)

• I0 – lowest order poloidal covariant component (defaults to 0)

• G1 – first order correction to toroidal covariant component (defaults to 0)

• I1 – first order correction to poloidal covariant component (defaults to 0)

• K1 – first order correction to radial covariant component (defaults to 0)

set_B0(B0)

set_Bbar(Bbar)

set_G0(G0)

set_G1(G1)

set_I0(I0)

set_I1(I1)

set_K1(K1)

set_N(N)

set_etabar(etabar)

set_iota0(iota0)

set_psi0(psi0)

class simsopt.field.boozermagneticfield.BoozerMagneticField(psi0)

Bases: BoozerMagneticField

Generic class that represents a magnetic field in Boozer coordinates \((s,\theta,\zeta)\). Here \(s = \psi/\psi_0\) is the normalized toroidal flux where \(2\pi\psi_0\) is the toroidal flux at the boundary. The magnetic field in the covariant form is,

\[\textbf B(s,\theta,\zeta) = G(s) \nabla \zeta + I(s) \nabla \theta + K(s,\theta,\zeta) \nabla \psi,\]

and the contravariant form is,

\[\textbf B(s,\theta,\zeta) = \frac{1}{\sqrt{g}} \left(\frac{\partial \mathbf r}{\partial \zeta} + \iota(s)\frac{\partial \mathbf r}{\partial \theta}\right),\]

where,

\[\sqrt{g}(s,\theta,\zeta) = \frac{G(s) + \iota(s)I(s)}{B^2}.\]

Here \(\iota(s) = \psi_P'(\psi)\) where \(2\pi\psi_P\) is the poloidal flux and \(2\pi\psi\) is the toroidal flux. Each subclass of BoozerMagneticField implements functions to compute \(B\), \(G\), \(I\), \(\iota\), \(\psi_P\), and their derivatives. The cylindrical coordinates \(R(s,\theta,\zeta)\) and \(Z(s,\theta,\zeta)\) in addition to \(K(s,\theta,\zeta)\) and \(\nu\) where \(\zeta = \phi + \nu(s,\theta,\zeta)\) and \(\phi\) is the cylindrical azimuthal angle are also implemented by BoozerRadialInterpolant and InterpolatedBoozerField. The usage is similar to the MagneticField class.

The usage of BoozerMagneticField` is as follows:

booz = BoozerAnalytic(etabar,B0,N,G0,psi0,iota0) # An instance of BoozerMagneticField
points = ... # points is a (n, 3) numpy array defining :math:`(s,\theta,\zeta)`
booz.set_points(points)
modB = bfield.modB() # returns the magnetic field strength at `points`

BoozerMagneticField has a cache to avoid repeated calculations. To clear this cache manually, call the clear_cached_properties() function. The cache is automatically cleared when set_points() is called or one of the dependencies changes.

clear_cached_properties()

Clear the cache.

recompute_bell(parent=None)

class simsopt.field.boozermagneticfield.BoozerRadialInterpolant(equil, order, mpol=32, ntor=32, N=None, enforce_vacuum=False, rescale=False, ns_delete=0, no_K=False)

Bases: BoozerMagneticField

Given a Vmec instance, performs a Boozer coordinate transformation using BOOZXFORM. The magnetic field can be computed at any point in Boozer coordinates using radial spline interpolation (scipy.interpolate.InterpolatedUnivariateSpline) and an inverse Fourier transform in the two angles. Throughout stellarator symmetry is assumed.

Parameters:

• equil – instance of simsopt.mhd.vmec.Vmec or simsopt.mhd.boozer.Boozer. If it is an instance of simsopt.mhd.boozer.Boozer, the compute_surfs needs to include all of the grid points in the half-radius grid of the corresponding Vmec equilibrium.

• order – (int) order for radial interpolation. Must satisfy 1 <= order <= 5.

• mpol – (int) number of poloidal mode numbers for BOOZXFORM (defaults to 32)

• ntor – (int) number of toroidal mode numbers for BOOZXFORM (defaults to 32)

• N – Helicity of quasisymmetry to enforce. If specified, then the non-symmetric Fourier harmonics of \(B\) and \(K\) are filtered out. Otherwise, all harmonics are kept. (defaults to None)

• enforce_vacuum – If True, a vacuum field is assumed, \(G\) is set to its mean value, \(I = 0\), and \(K = 0\).

• rescale –

  If True, use the interpolation method in the DELTA5D code. Here, a few of the first radial grid points or (bmnc, rmnc, zmns, numns, kmns) are deleted (determined by ns_delete). The Fourier harmonics are then rescaled as:

  bmnc(s)/s^(1/2) for m = 1

  bmnc(s)/s for m even and >= 2

  bmnc(s)/s^(3/2) for m odd and >=3

  before performing interpolation and spline differentiation to obtain dbmncds. If False, interpolation of the unscaled Fourier harmonics and its finite-difference derivative wrt s is performed instead (defaults to False)

• ns_delete – (see rescale) (defaults to 0)

compute_K()

init_splines()

class simsopt.field.boozermagneticfield.InterpolatedBoozerField(field, degree, srange, thetarange, zetarange, extrapolate=True, nfp=1, stellsym=True)

Bases: InterpolatedBoozerField, BoozerMagneticField

This field takes an existing BoozerMagneticField and interpolates it on a regular grid in \(s,\theta,\zeta\). This resulting interpolant can then be evaluated very quickly. This is modeled after InterpolatedField.

__init__(field, degree, srange, thetarange, zetarange, extrapolate=True, nfp=1, stellsym=True)

Parameters:

• field – the underlying simsopt.field.boozermagneticfield.BoozerMagneticField to be interpolated.

• degree – the degree of the piecewise polynomial interpolant.

• srange – a 3-tuple of the form (smin, smax, ns). This mean that the interval [smin, smax] is split into ns many subintervals.

• thetarange – a 3-tuple of the form (thetamin, thetamax, ntheta). thetamin must be >= 0, and thetamax must be <=2*pi.

• zetarange – a 3-tuple of the form (zetamin, zetamax, nzeta). zetamin must be >= 0, and thetamax must be <=2*pi.

• extrapolate – whether to extrapolate the field when evaluate outside the integration domain or to throw an error.

• nfp – Whether to exploit rotational symmetry. In this case any toroidal angle is always mapped into the interval \([0, 2\pi/\mathrm{nfp})\), hence it makes sense to use zetamin=0 and zetamax=2*np.pi/nfp.

• stellsym – Whether to exploit stellarator symmetry. In this case theta is always mapped to the interval \([0, \pi]\), hence it makes sense to use thetamin=0 and thetamax=np.pi.

simsopt.field.coil module

class simsopt.field.coil.CircularRegularizedCoil(curve, current, a)

Bases: RegularizedCoil

A coil with a circular cross section. The regularization parameter is computed from the radius during initialization.

Parameters:

• curve (simsopt.geo.curve.Curve) – The geometric curve describing the coil shape.

• current (Current) – The current object describing the electric current in the coil.

• a (float) – The radius of the circular cross-section.

class simsopt.field.coil.Coil(curve, current)

Bases: Coil, Optimizable

Represents a magnetic coil as a combination of a geometric curve and an electric current.

This class combines a Curve and a Current object, and is used as input for BiotSavart field calculations.

Parameters:

• curve (simsopt.geo.curve.Curve) – The geometric curve describing the coil shape.

• current (Current) – The current object describing the electric current in the coil.

plot(**kwargs)

Plot the coil’s curve. This method is just shorthand for calling the plot() function on the underlying Curve. All arguments are passed to simsopt.geo.curve.Curve.plot()

Parameters:

kwargs (dictionary) – Additional keyword arguments.

vjp(v_gamma, v_gammadash, v_current)

Compute the vector-Jacobian product,

\[\frac{\partial \mathbf{\gamma}}{\partial \mathbf{x}}^T \mathbf{v}_\gamma + \frac{\partial \mathbf{\gamma'}}{\partial \mathbf{x}}^T \mathbf{v}_{\gamma'} + \frac{\partial \mathbf{I}}{\partial \mathbf{x}}^T \mathbf{v}_I\]

where \(\mathbf{x}\) are the degrees of freedom of the coil.

Parameters:

• v_gamma (array, shape (n, 3)) – Vector w.r.t. \(\gamma\); same shape as curve gamma.

• v_gammadash (array, shape (n, 3)) – Vector w.r.t. \(\gamma'\); same shape as curve gammadash.

• v_current (array, shape (1,)) – Vector w.r.t. coil current (scalar).

Returns:

The vector-Jacobian product of the coil.

class simsopt.field.coil.Current(current, dofs=None, **kwargs)

Bases: Current, CurrentBase

An optimizable object that wraps around a single scalar degree of freedom representing an electric current.

This class is used for the current in a coil, or in a set of coils constrained to use the same current.

Parameters:

• current (float) – Initial value of the current.

• dofs (array-like or None, optional) – Degrees of freedom for optimization. If None, uses the current value.

• kwargs (dictionary) – Additional keyword arguments.

property current

Get the current value of the current object.

Returns:

The current value of the current object.

vjp(v_current)

Compute the Jacobian-vector product of the current function.

Parameters:

v_current (array, shape (1,)) – The vector to multiply the Jacobian with.

Returns:

The Jacobian-vector product of the current function.

class simsopt.field.coil.RectangularRegularizedCoil(curve, current, a, b)

Bases: RegularizedCoil

A coil with a rectangular cross section. The regularization parameter is computed from the width and height during initialization.

Parameters:

• curve (simsopt.geo.curve.Curve) – The geometric curve describing the coil shape.

• current (Current) – The current object describing the electric current in the coil.

• a (float) – The width of the rectangular cross-section.

• b (float) – The height of the rectangular cross-section.

class simsopt.field.coil.RegularizedCoil(curve, current, regularization)

Bases: Coil

A coil with a model for its cross section. This cross section is used to compute the forces and torques on the coil.

Parameters:

• curve (simsopt.geo.curve.Curve) – The geometric curve describing the coil shape.

• current (Current) – The current object describing the electric current in the coil.

• regularization (float) – The regularization parameter for the coil cross section.

B_regularized()

Calculate the regularized field on this coil following the Landreman and Hurwitz method.

Returns:

The regularized field on the coil.

Return type:

array (shape (n,3))

static _coil_force_pure(B, I, t)

Compute the pointwise Lorentz force per unit length on a coil with n quadrature points, in Newtons/meter.

\[dF/d\ell = I \vec{t} \times \vec{B}\]

where \(\vec{t}\) is the tangent vector to the coil curve, \(\vec{B}\) is the magnetic field at the quadrature points, \(I\) is the coil current.

Parameters:

• B (array, shape (n,3)) – Array of magnetic field.

• I (float) – Coil current.

• t (array, shape (n,3)) – Array of coil tangent vectors.

Returns:

Array of force per unit length.

Return type:

array (shape (n,3))

force(source_coils)

Compute the force per unit length on this coil from other coils, in Newtons/meter.

\[dF_i/d\ell = I_i \vec{t_i} \times (\vec{B_{self}} + \vec{B_{mutual}})\]

where \(\vec{t_i}\) is the tangent vector to the ith coil curve, \(\vec{B_{self}}\) is the self-field of the ith coil, \(\vec{B_{mutual}}\) is the mutual field from the other coils.

Parameters:

source_coils (list of Coil or RegularizedCoil, shape (m,)) – List of coils contributing forces on this coil. Can be a mix of Coil and RegularizedCoil objects.

Returns:

Array of forces per unit length along the coil curve.

Return type:

array (shape (n,3))

net_force(source_coils)

Compute the net forces on this coil from other coils, in Newtons. This is the integrated pointwise force per unit length dF_i/dell on the coil curve.

\[F_net = \int (dF_i/d\ell) d\ell\]

Parameters:

source_coils (list of Coil or RegularizedCoil, shape (m,)) – List of coils contributing forces on this coil. Can be a mix of Coil and RegularizedCoil objects.

Returns:

Array of net forces.

Return type:

np.array (shape (3,))

net_torque(source_coils)

Compute the net torques on this coil from other coils, in Newton-meters. This is the integrated pointwise torque per unit length on the coil curve.

\[T_net = \int dT_i/d\ell d\ell\]

Parameters:

source_coils (list of Coil or RegularizedCoil, shape (m,)) – List of coils contributing torques on this coil. Can be a mix of Coil and RegularizedCoil objects.

Returns:

Array of net torques.

Return type:

np.array (shape (3,))

self_force()

Compute the self-force per unit length of this coil, in Newtons/meter.

Returns:

Array of self-force per unit length.

Return type:

array (shape (n,3))

torque(source_coils)

Compute the torques per unit length on this coil from other coils in Newtons (note that the force is per unit length, so the force has units of Newtons/meter and the torques per unit length have units of Newtons).

\[dT_i/d\ell = (\gamma_i - c_i) \times (dF_i/d\ell)\]

where \(\gamma_i\) is the position vector of the ith coil curve, \(c_i\) is the centroid of the ith coil curve, \(dF_i/d\ell\) is the pointwise force per unit length on the ith coil curve.

Parameters:

source_coils (list of Coil or RegularizedCoil, shape (m,)) – List of coils contributing torques on this coil. Can be a mix of Coil and RegularizedCoil objects.

Returns:

Array of torques per unit length along the coil curve.

Return type:

np.array (shape (n,3))

simsopt.field.coil.apply_symmetries_to_currents(base_currents, nfp, stellsym)

Generate a list of currents by applying rotational and (optionally) stellarator symmetries.

Take a list of n Current`s and return ``n * nfp * (1+int(stellsym))` Current objects obtained by copying (for nfp rotations) and sign-flipping (optionally for stellarator symmetry).

Parameters:

• base_currents (list of Current) – List of base current objects to replicate.

• nfp (int) – Number of field periods (rotational symmetry).

• stellsym (bool) – Whether to apply stellarator symmetry (sign flip).

Returns:

List of current objects with symmetries applied.

Return type:

currents (list of Current)

simsopt.field.coil.apply_symmetries_to_curves(base_curves, nfp, stellsym)

Generate a list of curves by applying rotational and (optionally) stellarator symmetries.

Take a list of n simsopt.geo.curve.Curve`s and return ``n * nfp * (1+int(stellsym))` simsopt.geo.curve.Curve objects obtained by applying rotations and flipping corresponding to nfp fold rotational symmetry and optionally stellarator symmetry.

Parameters:

• base_curves (list) – List of base curves to replicate.

• nfp (int) – Number of field periods (rotational symmetry).

• stellsym (bool) – Whether to apply stellarator symmetry (flipping).

Returns:

List of curves with symmetries applied.

Return type:

curves (list)

simsopt.field.coil.coils_to_focus(filename, curves, currents, nfp=1, stellsym=False, Ifree=False, Lfree=False)

Export a list of CurveXYZFourier objects and currents to a FOCUS input file.

The output can be used by FOCUS. The format is described at https://princetonuniversity.github.io/FOCUS/rdcoils.pdf

Parameters:

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

• curves (list) – list of CurveXYZFourier objects.

• currents (list) – list of Current objects.

• nfp (int, optional) – Number of field periods (default: 1).

• stellsym (bool, optional) – Whether to apply stellarator symmetry (default: False).

• Ifree (bool, optional) – Whether the coil current is free (default: False).

• Lfree (bool, optional) – Whether the coil geometry is free (default: False).

simsopt.field.coil.coils_to_makegrid(filename, curves, currents, groups=None, nfp=1, stellsym=False)

Export a list of Curve objects and currents to a mgrid (MAKEGRID) input file.

The output can be used by MAKEGRID and FOCUS. The format is described at https://princetonuniversity.github.io/STELLOPT/MAKEGRID

Parameters:

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

• curves (list) – list of Curve objects.

• currents (list) – list of Current objects.

• groups (list or None, optional) – Coil current group. Coils in the same group are assembled together.

• nfp (int, optional) – Number of field periods (default: 1).

• stellsym (bool, optional) – Whether to apply stellarator symmetry (default: False).

simsopt.field.coil.coils_to_vtk(coils, filename, close=False, extra_data=None)

Export a list of Coil 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.

Saves coil currents, net forces, net torques, and pointwise forces and torques.

Parameters:

• coils (list) – A python list of Coil objects.

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

• close (bool) – Whether to draw the segment from the last quadrature point back to the first.

• extra_data (dict) – Additional data to save to the VTK file.

simsopt.field.coil.coils_via_symmetries(curves, currents, nfp, stellsym, regularizations=None)

Generate a list of Coil objects by applying rotational and (optionally) stellarator symmetries.

Take a list of n curves and return n * nfp * (1+int(stellsym)) Coil objects obtained by applying rotations and flipping corresponding to nfp fold rotational symmetry and optionally stellarator symmetry.

If regularizations are provided for the base curves, then RegularizedCoil objects are returned. This is a coil class that carries around a regularization for a finite cross section, which is used for computing e.g. forces and torques on the coil. Format is e.g. regularizations = [regularization_circ(0.05) for _ in range(ncoils)]

Parameters:

• curves (list, shape (n_coils,)) – List of base curves.

• currents (list, shape (n_coils,)) – List of base current objects.

• nfp (int) – Number of field periods (rotational symmetry).

• stellsym (bool) – Whether to apply stellarator symmetry.

• regularizations (np.array, shape (n_coils,), optional) – The regularization objects for the coils representing the finite coil cross section.

Returns:

List of Coil or RegularizedCoil objects with symmetries applied. If regularizations are provided, then RegularizedCoil objects are returned.

Return type:

coils (list)

simsopt.field.coil.load_coils_from_makegrid_file(filename, order, ppp=20, group_names=None)

Load coils from a mgrid input file, returning a list of Coil objects.

This function loads a file in MAKEGRID input format containing the Cartesian coordinates and the currents for several coils and returns an array with the corresponding coils. The format is described at https://princetonuniversity.github.io/STELLOPT/MAKEGRID

Parameters:

• filename (str) – Path to the MAKEGRID input file.

• order (int) – Maximum mode number in the Fourier expansion.

• ppp (int, optional) – Points per period for quadrature (default: 20).

• group_names (list of str or str or None, optional) – If provided, only load coils in these groups.

Returns:

List of Coil objects loaded from the file.

Return type:

coils (list)

simsopt.field.coilset module

class simsopt.field.coilset.CoilSet(base_coils=None, coils=None, surface=None)

Bases: Optimizable

A set of coils as a single optimizable object, and a surface on which it evaluates the field.

The surfaces will be cast as a SurfaceRZFourier, as it will evaluate a field on the surface using the RZFourier/VMEC convention. The surfaces range will be adapted to ‘half period’ or ‘field period’ respectively if it is stellarator symmetric or not.

The Coilset contains convenience methods to acces coil optimization penalties: flux_penalty, length_penalty, cc_distance_penalty, cs_distance_penalty, lp_curvature_penalty, meansquared_curvature_penalty, arc_length_variation_penalty and total_length.

These functions can then be used to define a FOCUS-like optimization problem (note that the FOCUS algorithm is not a least-squares algorithm) that can be passed to scipy.minimize, or the CoilSet can be used as a parent for a SPEC NormalField object.

static _circlecurves_around_surface(surf, coils_per_period=4, order=6, R0=None, R1=None, use_stellsym=None, factor=2.0)

return a set of base curves for a surface using the surfaces properties where possible

arc_length_variation_penalty()

Return a penalty function on the arc length variation of the coils

property base_coils

cc_distance_penalty(DISTANCE_THRESHOLD)

Return a penalty function for the distance between coils :param DISTANCE_THRESHOLD: The threshold distance below which the penalty is applied

cs_distance_penalty(DISTANCE_THRESHOLD)

Return a penalty function for the distance between coils and the surface :param DISTANCE_THRESHOLD: The threshold distance below which the penalty is applied

flux_penalty(target=None)

Return the penalty function for the quadratic flux penalty on the surface

classmethod for_surface(surf, coil_current=100000.0, coils_per_period=5, nfp=None, current_constraint='fix_all', **kwargs)

Create a CoilSet for a given surface. The coils are created using simsopt.geo.create_equally_spaced_curves() with the given parameters.

The keyword arguments coils_per_period, order, R0, R1, factor, and use_stellsym are passed to the create_equally_spaced_curves() function.

Parameters:

• surf – The surface for which to create the coils

• total_current – the total current in the CoilSet

• coils_per_period – the number of coils per field period

• nfp – The number of field periods.

• current_constraint – "fix_one" or "fix_all" or "free_all"

• coils_per_period – The number of coils per field period

• order – The order of the Fourier expansion

• R0 – major radius of a torus on which the initial coils are placed

• R1 – The radius of the coils

• factor – if R0 or R1 are None, they are factor times the first sine/cosine coefficient (factor > 1)

• use_stellsym – Whether to use stellarator symmetry

classmethod from_makegrid_file(makegrid_file, surface, order=10, ppp=10)

Create a CoilSet from a MAKEGRID file and a surface.

get_dof_orders()

get the Fourier order of the degrees of freedom corresponding to the Fourier coefficients describing the coils. Useful for reducing the trust region for the higher order coefficients. in optimization.

Returns:

An array with the fourier order of each degree of freedom.

length_penalty(TOTAL_LENGTH, f)

Return a QuadraticPenalty on the total length of the coils if it is larger than TARGET_LENGTH (do not penalize if shorter) :param TOTAL_LENGTH: The threshold length above which the penalty is applied :param f: type of penalty, “min”, “max” or “identity”

lp_curvature_penalty(CURVATURE_THRESHOLD, p=2)

Return a penalty function on the curvature of the coils that is the Lp norm of the curvature of the coils. Defaults to L2. :param CURVATURE_THRESHOLD: The threshold curvature above which the penalty is applied :param p: The p in the Lp norm for calculating the penalty

meansquared_curvature_penalty()

Return a penalty function on the mean squared curvature of the coils

meansquared_curvature_threshold(CURVATURE_THRESHOLD)

Return a penalty function on the mean squared curvature of the coils

plot(**kwargs)

Plot the coil’s curve. This method is just shorthand for calling the plot() function on the underlying Curve. All arguments are passed to simsopt.geo.Curve/Surface.plot()

reduce(target_function, nsv='nonzero')

Return a ReducedCoilSet based on this CoilSet using the SVD of the mapping given by target_function.

Parameters:

• target_function – a function that takes a CoilSet and maps to a different space that is relevant for the optimization problem. The SVD will be calculated on the Jacobian of f: coilDOFs -> target. For example a function that calculates the fourier coefficients of the normal field on the surface.

• nsv – The number of singular vectors to keep, integer or "nonzero" (for all nonzero singular values)

property surface

to_makegrid_file(filename)

Write the CoilSet to a MAKEGRID file

to_vtk(filename, add_biotsavart=True, close=False)

Write the CoilSet and its surface to a vtk file

property total_length

Return the total length of the coils

total_length_penalty()

Return the total length of the coils

class simsopt.field.coilset.ReducedCoilSet(coilset=None, nsv=None, s_diag=None, u_matrix=None, vh_matrix=None, function=None, **kwargs)

Bases: CoilSet

An Optimizable that replaces a CoilSet and has as degrees of freedom linear combinations of the CoilSet degrees of freedom

Single-stage optimization is often accompanied with an explosion in the numbers of degrees of freedom needed to accurately represent the stellarator, as each coilneeds 6 DOFs per coil per Fourier order. This makes solving the optimization problem harder in two ways: 1: the computational cost of finite-difference gradient evaluation is linear in the number of DOFs, and 2: the optimization problem can become rank-deficient: changes to your DOFs cancel each other, and the minmum is no longer unique.

This class solves these problems by allowing physics-based dimensionality reduction of the optimization space using singluar value-decomposition of a fast-evalable map.

A ReducedCoilSet consists of a Coilset, and the three elements of an SVD of the Jacobian defined by a mapping from the coils’ DOFs to a physically relevant, fast-evaluating target.

__init__(coilset=None, nsv=None, s_diag=None, u_matrix=None, vh_matrix=None, function=None, **kwargs)

create a ReducedCoilSet from a CoilSet and the three elements of an SVD of the Jacobian of a mapping.

property base_coils

property coils

property coilset

classmethod from_function(coilset, target_function, nsv='nonzero')

Reduce an existing CoilSet to a ReducedCoilSet using the SVD of the mapping given by target_function.

Parameters:

• coilset – The CoilSet to reduce

• target_function – a function that takes a CoilSet and maps to a different space that is relevant for the optimization problem. The SVD will be calculated on the Jacobian of f: coilDOFs -> target For example a function that calculates the fourier coefficients of the normal field on the surface.

• nsv – The number of singular values to keep, either an integer or "nonzero". defaults to "nonzero". If an integer is given, this number of largest singular values will be kept. If "nonzero", all nonzero singular values are kept.

get_dof_orders()

Not available for a ReducedCoilSet

property lsv

left-singular vectors of the svd

property nsv

number of singular values used by the ReducedCoilSet

plot_singular_vector(n, eps=0.0001, show_delta_B=True, engine='mayavi', show=False, **kwargs)

Plot the n-th singular vector. :param n: the index of the singular vector :param eps: the magnitude of the displacement. :param show_delta_B: whether to show the change in the magnetic field :param engine: the plotting engine to use. Defaults to ‘mayavi’

recalculate_reduced_basis(target_function=None)

Recalculate the reduced basis using a new target function.

property rsv

right-singular vectors of the svd

set_dofs(x)

property singular_values

singular values of the svd

property surface

simsopt.field.force module

Implements the force on a coil in its own magnetic field and the field of other coils.

class simsopt.field.force.B2Energy(target_coils, downsample=1)

Bases: Optimizable

Optimizable class for minimizing the total vacuum magnetic field energy from a set of m coils.

The function is

\[J = \frac{1}{2}\sum_{i,j}I_i L_{ij} I_j\]

where \(L_{ij}\) is the coil inductance matrix (positive definite), and \(I_i\) is the current in the ith coil. The units of the objective function are MJ (megajoules).

Parameters:

• target_coils (list of RegularizedCoil, shape (m,)) – List of coils contributing to the total energy.

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

J()

Evaluate the B^2 energy objective.

Returns:

The total vacuum magnetic field energy

\(J = \frac{1}{2}\sum_{i,j} I_i L_{ij} I_j\) in MJ.

Return type:

float

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

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

class simsopt.field.force.LpCurveForce(target_coils, source_coils_coarse, source_coils_fine=None, p: float = 2.0, threshold: float = 0.0, downsample: int = 1)

Bases: Optimizable

Optimizable class to minimize the total Lp-Lorentz force density (force per unit length) integrated. Force density on a coil is computed on each coil in a set of m target coils, using the self-force from the coil itself, the force from the other m - 1 target coils and the force from a set of m’ source coils. If source_coils_coarse and target_coils have coils in common, they are removed during initialization of this class, to avoid double counting forces. A typical use case has the target_coils as the unique base_coils in a stellarator optimization, and source_coils_coarse are all the coils after applying symmetries. Typical initialization is LpCurveForce(base_coils, coils).

The objective function is

\[J = \frac{1}{p}\sum_i\frac{1}{L_i}\left(\int \text{max}(|d\vec{F}/d\ell_i| - F_0 , 0)^p d\ell_i\right)\]

where \(\frac{d\vec{F}_i}{d\ell_i}\) is the Lorentz force per unit length, in units of MN/m, where MN = meganewtons. The units of the objective function are therefore (MN/m)^p. \(d\ell_i\) is the arclength along the ith coil, \(L_i\) is the total coil length, and :math:`F_0 ` is a threshold force at the ith coil.

This class assumes there are two (or three) distinct lists of coils, which may have different finite-build parameters and/or different numbers of quadrature points. In order to avoid buildup of optimizable dependencies, it directly computes the BiotSavart law terms, instead of relying on the existing C++ code that computes BiotSavart related terms. Within each list of coils, all coils must have the same number of quadrature points. The source_coils_coarse and source_coils_fine lists allows one to optimize e.g. the torque on target_coils from a set of dipole coils (with barely any quadrature points) and a set of TF coils (with many quadrature points).

Parameters:

• target_coils (list of RegularizedCoil, shape (m,)) – List of coils on which the LpCurveForce is computed.

• source_coils_coarse (list of Coil or RegularizedCoil, shape (m',)) – Coarse-resolution source coils that provide forces on the target_coils. Forces are not computed on the source_coils.

• source_coils_fine (list of Coil or RegularizedCoil, optional) – Fine-resolution source coils, used in addition to coarse. Default: []. This functionality is provided for when there are two sets of source coils with very different numbers of quadrature points. This occurs e.g. when optimizing TF coils and dipole coils.

• p (float) – Power of the objective function.

• threshold (float) – Threshold force per unit length in units of MN/m (meganewtons per meter).

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

J()

Evaluate the Lp curve force objective.

_J_args()

Build arguments for evaluation of J and dJ.

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

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

class simsopt.field.force.LpCurveTorque(target_coils, source_coils_coarse, source_coils_fine=None, p: float = 2.0, threshold: float = 0.0, downsample: int = 1)

Bases: Optimizable

Optimizable class to minimize the total Lp-Lorentz torque density (torque per unit length) integrated. Torque density on a coil is computed on each coil in a set of m target coils, using the self-force from the coil itself, the force from the other m - 1 target coils and the force from a set of m’ source coils. If source_coils and target_coils have coils in common, they are removed during initialization of this class, to avoid double counting forces. A typical use case has the target_coils as the unique base_coils in a stellarator optimization, and source_coils are all the coils after applying symmetries. Typical initialization is LpCurveTorque(base_coils, coils).

The objective function is

\[J = \frac{1}{p}\sum_i\frac{1}{L_i}\left(\int \text{max}(|d\vec{T}/d\ell_i| - T_0 , 0)^p d\ell_i\right)\]

where \(\frac{d\vec{T}_i}{d\ell_i}\) is the Lorentz torque per unit length, in units of MN, where MN = meganewtons. The units of the objective function are therefore (MN)^p. \(d\ell_i\) is the arclength along the ith coil, \(L_i\) is the total coil length, and :math:`T_0 ` is a threshold torque per unit length at the ith coil.

This class assumes there are two (or three) distinct lists of coils, which may have different finite-build parameters and/or different numbers of quadrature points. In order to avoid buildup of optimizable dependencies, it directly computes the BiotSavart law terms, instead of relying on the existing C++ code that computes BiotSavart related terms. Within each list of coils, all coils must have the same number of quadrature points. The source_coils_coarse and source_coils_fine lists allows one to optimize e.g. the torque on target_coils from a set of dipole coils (with barely any quadrature points) and a set of TF coils (with many quadrature points).

Parameters:

• target_coils (list of RegularizedCoil, shape (m,)) – List of coils to use for computing LpCurveTorque.

• source_coils_coarse (list of Coil or RegularizedCoil, shape (m',)) – Coarse-resolution source coils that provide torques on the target_coils. Torques are not computed on the source_coils.

• source_coils_fine (list of Coil or RegularizedCoil, optional) – Fine-resolution source coils, used in addition to coarse. Default: []. This functionality is provided for when there are two sets of source coils with very different numbers of quadrature points. This occurs e.g. when optimizing TF coils and dipole coils.

• p (float) – Power of the objective function.

• threshold (float) – Threshold torque per unit length in units of MN (meganewtons).

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

J()

Evaluate the Lp curve torque objective.

_J_args()

Build arguments for evaluation of J and dJ.

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

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

class simsopt.field.force.NetFluxes(target_coil, source_coils, downsample=1)

Bases: Optimizable

Optimizable class for minimizing the total net flux from m coils through a single coil with n quadrature points.

The function is

\[\Psi = \int A_{ext}\cdot d\ell / L\]

where \(A_{ext}\) is the vector potential of an external magnetic field, evaluated along the quadpoints along the curve, L is the total length of the coil, and \(\ell\) is arclength along the coil.

The units of the objective function are Weber.

Parameters:

• target_coil (Coil) – Coil whose net flux is being computed.

• source_coils (list of Coil, shape (m,)) – List of coils to use for computing the net flux.

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

J()

Evaluate the net flux objective.

Computes \(\Psi = \int A_{ext} \cdot d\ell / L\) using the BiotSavart vector potential from the source coils evaluated at the target coil quadrature points.

Returns:

Net magnetic flux through the target coil in Weber.

Return type:

float

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

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

class simsopt.field.force.SquaredMeanForce(target_coils, source_coils_coarse, source_coils_fine=None, downsample: int = 1)

Bases: Optimizable

Optimizable class to minimize the (net (integrated) Lorentz force per unit length)^2 on a set of m coils from themselves and another set of m’ coils.

The objective function is

where \(\frac{d\vec{F}_i}{d\ell_i}\) is the Lorentz force per unit length, in units of MN/m. The units of the squared mean force are therefore (MN/m)^2. \(L_i\) is the total coil length, and \(\ell_i\) is arclength along the ith coil. The units of the objective function are (MN/m)^2, where MN = meganewtons.

This class assumes there are two (or three) distinct lists of coils, which may have different finite-build parameters and/or different numbers of quadrature points. In order to avoid buildup of optimizable dependencies, it directly computes the BiotSavart law terms, instead of relying on the existing C++ code that computes BiotSavart related terms. This is also useful for optimizing passive coils, which require a modified Jacobian calculation. Within each list of coils, all coils must have the same number of quadrature points. The source_coils_coarse and source_coils_fine lists allows one to optimize e.g. the force on target_coils from a set of dipole coils (with barely any quadrature points) and a set of TF coils (with many quadrature points).

Parameters:

• target_coils (list of Coil or RegularizedCoil, shape (m,)) – List of coils to use for computing SquaredMeanForce.

• source_coils_coarse (list of Coil or RegularizedCoil, shape (m',)) – Coarse-resolution source coils that provide forces on the target_coils. Forces are not computed on the source_coils.

• source_coils_fine (list of Coil or RegularizedCoil, optional) – Fine-resolution source coils, used in addition to coarse. Default: []. This functionality is provided for when there are two sets of source coils with very different numbers of quadrature points. This occurs e.g. when optimizing TF coils and dipole coils.

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

J()

Evaluate the squared mean force objective.

_J_args()

Build arguments for evaluation of J and dJ.

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

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

class simsopt.field.force.SquaredMeanTorque(target_coils, source_coils_coarse, source_coils_fine=None, downsample: int = 1)

Bases: Optimizable

Optimizable class to minimize the (net (integrated) Lorentz torque per unit length)^2 summed over a set of m coils due to themselves and another set of m’ coils.

The objective function is

where \(\frac{d\vec{T}_i}{d\ell_i}\) is the Lorentz torque per unit length, in units of MN. The units of the squared mean torque are therefore (MN)^2. \(d\ell_i\) is the arclength along the ith coil, \(L_i\) is the total coil length.

The units of the objective function are (MN)^2, where MN = meganewtons.

This class assumes there are two (or three) distinct lists of coils, which may have different finite-build parameters and/or different numbers of quadrature points. In order to avoid buildup of optimizable dependencies, it directly computes the BiotSavart law terms, instead of relying on the existing C++ code that computes BiotSavart related terms. Within each list of coils, all coils must have the same number of quadrature points. The source_coils_coarse and source_coils_fine lists allows one to optimize e.g. the torque on target_coils from a set of dipole coils (with barely any quadrature points) and a set of TF coils (with many quadrature points).

Parameters:

• target_coils (list of Coil or RegularizedCoil, shape (m,)) – List of coils to use for computing SquaredMeanTorque.

• source_coils_coarse (list of Coil or RegularizedCoil, shape (m',)) – Coarse-resolution source coils that provide torques on the target_coils. Torques are not computed on the source_coils.

• source_coils_fine (list of Coil or RegularizedCoil, optional) – Fine-resolution source coils, used in addition to coarse. Default: []. This functionality is provided for when there are two sets of source coils with very different numbers of quadrature points. This occurs e.g. when optimizing TF coils and dipole coils.

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

Returns:

Value of the objective function.

Return type:

float

J()

Evaluate the squared mean torque objective.

_J_args()

Build arguments for evaluation of J and dJ.

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

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

simsopt.field.force._B_at_point_from_coil_set_pure(pt, gammas, gammadashs, currents, exclude_index, eps)

Compute the magnetic field at a single point due to a set of coils via the Biot-Savart law, optionally excluding one coil (e.g. to avoid self-contribution).

This is a pure JAX implementation of the Biot-Savart integral used by force and torque objectives in this module. We do not use the BiotSavart class here because constructing one BiotSavart per objective (e.g. per coil in LpCurveForce) leads to a large number of optimizable dependencies and weak references. That makes operations like Jf.x = dofs scale poorly with the number of coils (tens of millions of function calls and tens of seconds for ~64 coils). See GitHub issue #487.

\[B = \frac{\mu_0}{4\pi} \frac{1}{n_{pts}} \sum_{j \neq \mathrm{exclude}} I_j \int \frac{d\vec{\ell}_j \times (\vec{r} - \vec{r}_j)}{|\vec{r} - \vec{r}_j|^3}\]

Parameters:

• pt – Array of shape (3,); evaluation point.

• gammas – Array of shape (m, n, 3); positions for m coils with n quadrature points.

• gammadashs – Array of shape (m, n, 3); tangent vectors.

• currents – Array of shape (m,); coil currents.

• exclude_index – Index of coil to exclude from the sum (use -1 to include all).

• eps – Small constant added to distances to avoid division by zero.

Returns:

Array of shape (3,); magnetic field contribution (without mu_0/(4*pi)).

simsopt.field.force._check_downsample(coils, downsample, label='coils')

Check that downsample evenly divides the number of quadrature points.

Parameters:

• coils – list of coils to check (must be non-empty).

• downsample – downsampling factor.

• label – descriptive label for the coil group (used in error message).

Raises:

ValueError – if downsample does not evenly divide the number of quadrature points.

simsopt.field.force._check_quadpoints_consistency(coils, label='coils')

Check that all coils in a list have the same number of quadrature points.

Parameters:

• coils – list of coils to check.

• label – descriptive label for the coil group (used in error message).

Raises:

ValueError – if not all coils have the same number of quadrature points.

simsopt.field.force._coil_coil_inductances_inv_pure(gammas, gammadashs, downsample, regularizations)

Pure function for computing the inverse of the coil inductance matrix L. This matrix is symmetric positive definite by definition.

Performs a Cholesky decomposition of the coil inductance matrix L and then solves for the inverse. Note that inverse of a (nonsingular)lower triangular matrix C is upper triangular and vice versa.

\[L = (C C^T)\]

where \(C\) is a lower triangular matrix from the Cholesky decomposition of \(L\). Then, we solve two triangular systems of equations:

\[ \begin{align}\begin{aligned}C^{-1}C = I\\L^{-1} = (C C^T)^{-1} = (C^T)^{-1} C^{-1}\end{aligned}\end{align} \]

so that we can solve for \(L^{-1}\) by multiplying both sides by \(C^T\) and solving it as an upper triangular system,

\[C^TL^{-1} = C^{-1}\]

The units of the inverse of the coil inductance matrix are 1/H, where H = henries.

Parameters:

• gammas (array, shape (m,n,3)) – Array of coil positions for all m coils.

• gammadashs (array, shape (m,n,3)) – Array of coil tangent vectors for all m coils.

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

• regularizations (array, shape (m,)) – Array of regularizations coming from finite cross-section for all m coils. The choices for each coil are regularization_circ and regularization_rect, although each coil can have different size and shape cross-sections in this list of regularization terms.

Returns:

Array of inverse of the coil inductance matrix.

Return type:

array (shape (m,m))

simsopt.field.force._coil_coil_inductances_pure(gammas, gammadashs, downsample, regularizations, eps=1e-10)

Compute the full inductance matrix for a set of coils, including both mutual and self-inductances. All coils are assumed to have the same number of quadrature points, denoted n. The units of the inductance matrix are H, where H = henries.

The mutual inductance between two coils is computed as:

\[M = \frac{\mu_0}{4\pi} \iint \frac{d\vec{r}_A \cdot d\vec{r}_B}{|\vec{r}_A - \vec{r}_B|}\]

and self-inductance of a regularized coil is computed as:

\[L = \frac{\mu_0}{4\pi} \int_0^{2\pi} d\phi \int_0^{2\pi} d\tilde{\phi} \frac{\vec{r}_c' \cdot \tilde{\vec{r}}_c'}{\sqrt{|\vec{r}_c - \tilde{\vec{r}}_c|^2 + \delta a b}}\]

where $delta a b$ is a regularization parameter depending on the cross-section. The units of the inductance matrices are henries.

Parameters:

• gammas (array, shape (m,n,3)) – Array of coil positions for all m coils.

• gammadashs (array, shape (m,n,3)) – Array of coil tangent vectors for all m coils.

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

• regularizations (array, shape (m,)) – Array of regularizations coming from finite cross-section for all m coils. The choices for each coil are regularization_circ and regularization_rect, although each coil can have different size and shape cross-sections in this list of regularization terms.

• eps (float) – Small constant to avoid division by zero for mutual inductance between coil_i and itself.

Returns:

Full inductance matrix Lij.

Return type:

array (shape (m,m))

simsopt.field.force._induced_currents_pure(gammas_targets, gammadashs_targets, gammas_sources, gammadashs_sources, currents_sources, downsample, regularizations)

Pure function for computing the induced currents in a set of m passive coils with n quadrature points due to a set of m’ source coils with n’ quadrature points (and themselves).

\[I = -L^{-1} \Psi\]

where \(L\) is the coil inductance matrix, \(\Psi\) is the net flux through the passive coils due to the source coils, and \(I\) is the induced currents in the passive coils. The units of the induced currents are Amperes.

Parameters:

• gammas_targets (array, shape (m,n,3)) – Array of passive coil positions for all m coils.

• gammadashs_targets (array, shape (m,n,3)) – Array of passive coil tangent vectors for all m coils.

• gammas_sources (array, shape (m',n',3)) – Array of source coil positions for all m’ coils.

• gammadashs_sources (array, shape (m',n',3)) – Array of source coil tangent vectors for all m’ coils.

• currents_sources (array, shape (m',)) – Array of source coil currents.

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

• regularizations (array, shape (m,)) – Array of regularizations coming from finite cross-section for all m coils. The choices for each coil are regularization_circ and regularization_rect, although each coil can have different size and shape cross-sections in this list of regularization terms.

Returns:

Array of induced currents.

Return type:

array (shape (m,))

simsopt.field.force._lorentz_force_density_pure(tangents, current, magnetic_field)

Compute Lorentz force density I * (t x B).

simsopt.field.force._mutual_B_field_at_point_pure(i, pt, gammas_targets, gammadashs_targets, currents_targets, gammas_sources_coarse, gammadashs_sources_coarse, currents_sources_coarse, gammas_sources_fine, gammadashs_sources_fine, currents_sources_fine, eps)

Compute the mutual magnetic field at a point on target coil i from all target coils (excluding coil i) and all source coils (coarse and fine) in Tesla.

Used by squared_mean_force_pure(), lp_force_pure(), lp_torque_pure(), and squared_mean_torque(). See _B_at_point_from_coil_set_pure() for why Biot-Savart is reimplemented here instead of using BiotSavart.

Parameters:

• i – Index of target coil.

• pt – Array of shape (3,); evaluation point.

• gammas_targets – Array of shape (m, n, 3); positions for m target coils with n quadrature points.

• gammadashs_targets – Array of shape (m, n, 3); tangent vectors for m target coils with n quadrature points.

• currents_targets – Array of shape (m,); currents for m target coils.

• gammas_sources_coarse – Array of shape (m’, n’, 3); positions for m’ coarse source coils.

• gammadashs_sources_coarse – Array of shape (m’, n’, 3); tangent vectors for coarse source coils.

• currents_sources_coarse – Array of shape (m’,); currents for coarse source coils.

• gammas_sources_fine – Array of shape (m’’, n’’, 3); positions for m’’ fine source coils (may be empty).

• gammadashs_sources_fine – Tangent vectors for fine source coils.

• currents_sources_fine – Currents for fine source coils.

• eps – Small constant added to distances to avoid division by zero.

Returns:

Array of shape (3,); mutual magnetic field at point pt in Tesla.

simsopt.field.force._net_fluxes_pure(gammas_targets, gammadashs_targets, gammas_sources, gammadashs_sources, currents_sources, downsample)

This function computes the total magnetic flux passing through a set of coils due to the magnetic field generated by another set of coils. The flux is calculated using the line integral of the vector potential along the coil paths.

math::

Psi = sum_i int_{C_i} A_{ext}cdot dell_i / L_i

where \(A_{ext}\) is the vector potential of an external magnetic field, evaluated along the quadpoints along the curve, \(L_i\) is the total length of the ith coil, and \(\ell_i\) is arclength along the ith coil.

Note that the first set of coils is assumed to all have the same number of quadrature points for the purposes of jit speed. Same with the second set of coils, although the number of points does not have to be the same between the two sets.

The units of the objective function are Weber.

Parameters:

• gammas_targets (array, shape (m,n,3)) – Position vectors for the coils receiving flux.

• gammadashs_targets (array, shape (m,n,3)) – Tangent vectors for the coils receiving flux.

• gammas_sources (array, shape (m',n',3)) – Position vectors for the coils generating flux.

• gammadashs_sources (array, shape (m',n',3)) – Tangent vectors for the coils generating flux.

• currents_sources (array, shape (m',)) – Current values for the coils generating flux.

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

Returns:

Net magnetic flux through each coil in the first set.

Return type:

array (shape (m,))

simsopt.field.force._prepare_regularized_target_source_inputs_pure(gammas_targets, gammadashs_targets, gammadashdashs_targets, quadpoints, gammas_sources, gammadashs_sources, currents_targets, currents_sources, regularizations, downsample)

Downsample/convert inputs for regularized Lp force/torque objectives. Just a wrapper around _prepare_target_source_inputs_pure that also prepares additional inputs for regularized coils.

Parameters:

• gammas_targets – Array of shape (m, n, 3); positions for m target coils with n quadrature points.

• gammadashs_targets – Array of shape (m, n, 3); tangent vectors for m target coils with n quadrature points.

• gammadashdashs_targets – Array of shape (m, n, 3); second derivatives of tangent vectors for m target coils with n quadrature points.

• quadpoints – Array of shape (m, n); quadrature points for m target coils with n quadrature points.

• gammas_sources – Array of shape (m’, n, 3); positions for m’ source coils with n quadrature points.

• gammadashs_sources – Array of shape (m’, n, 3); tangent vectors for m’ source coils with n quadrature points.

• currents_targets – Array of shape (m,); currents for m target coils.

• currents_sources – Array of shape (m’,); currents for m’ source coils.

• regularizations – Array of shape (m,); regularizations for m target coils.

• downsample – Factor by which to downsample the quadrature points.

Returns:

(gammas_targets, gammadashs_targets, gammadashdashs_targets, quadpoints, gammas_sources, gammadashs_sources, currents_targets, currents_sources, regularizations).

Return type:

Tuple of arrays

simsopt.field.force._prepare_target_source_inputs_pure(gammas_targets, gammadashs_targets, gammas_sources, gammadashs_sources, currents_targets, currents_sources, downsample)

Downsample and convert shared target/source inputs used by force/torque objectives.

Parameters:

• gammas_targets – Array of shape (m, n, 3); positions for m target coils with n quadrature points.

• gammadashs_targets – Array of shape (m, n, 3); tangent vectors for m target coils with n quadrature points.

• gammas_sources – Array of shape (m’, n, 3); positions for m’ source coils with n quadrature points.

• gammadashs_sources – Array of shape (m’, n, 3); tangent vectors for m’ source coils with n quadrature points.

• currents_targets – Array of shape (m,); currents for m target coils.

• currents_sources – Array of shape (m’,); currents for m’ source coils.

• downsample – Factor by which to downsample the quadrature points.

Returns:

(gammas_targets, gammadashs_targets, gammas_sources, gammadashs_sources, currents_targets, currents_sources).

Return type:

Tuple of arrays

simsopt.field.magnetic_axis_helpers module

simsopt.field.magnetic_axis_helpers.compute_on_axis_iota(axis, magnetic_field)

Computes the rotational transform on the magnetic axis of a device using a method based on equation (13) of Greene, Journal of Mathematical Physics 20, 1183 (1979); doi: 10.1063/1.524170. This method was shared by Stuart Hudson and Matt Landreman. NOTE: this function does not check that the provided axis is actually a magnetic axis of the input magnetic_field.

Parameters:

• axis – a CurveRZFourier corresponding to the magnetic axis of the magnetic field.

• magnetic_field – the magnetic field in which the axis lies.

Returns:

the rotational transform on the given axis.

Return type:

iota

simsopt.field.magneticfield module

class simsopt.field.magneticfield.MagneticField(**kwargs)

Bases: MagneticField, Optimizable

Generic class that represents any magnetic field from which each magnetic field class inherits. The usage of MagneticField is as follows:

bfield = BiotSavart(coils) # An instance of a MagneticField
points = ... # points is a (n, 3) numpy array
bfield.set_points(points)
B = bfield.B() # returns the Magnetic field at `points`
dA = bfield.dA_by_dX() # returns the gradient of the potential of the field at `points`

MagneticField has a cache to avoid repeated calculations. To clear this cache manually, call the clear_cached_properties() function. The cache is automatically cleared when set_points is called or one of the dependencies changes.

__add__(other)

Add two magnetic fields.

__mul__(other)

Multiply a field with a scalar.

__rmul__(other)

Multiply a field with a scalar.

clear_cached_properties()

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

set_points(self: simsoptpp.MagneticField, arg0: numpy.ndarray[numpy.float64]) → simsoptpp.MagneticField

Shorthand for set_points_cart.

set_points_cart(self: simsoptpp.MagneticField, arg0: numpy.ndarray[numpy.float64]) → simsoptpp.MagneticField

Set the points where to evaluate the magnetic fields, in cartesian coordinates.

set_points_cyl(self: simsoptpp.MagneticField, arg0: numpy.ndarray[numpy.float64]) → simsoptpp.MagneticField

Set the points where to evaluate the magnetic fields, in cylindrical coordinates (the order is \((r, \phi, z)\)).

to_mgrid(filename, nr=10, nphi=4, nz=12, rmin=1.0, rmax=2.0, zmin=-0.5, zmax=0.5, nfp=1, include_potential=False)

Export the field to the mgrid (NetCDF) file format for free boundary calculations.

An mgrid file contains a grid in cylindrical coordinates (R, phi, z) upon which the magnetic field is evaluated. The grid is defined by the number of points in the radius (nr), the number of planes per field period in the toroidal angle (nphi), and the number of points in the z coordinate (nz). A good rule here is to choose at least 4 times as many toroidal planes as the maximum toroidal mode number (ex. if ntor=6 then you should have at least 24 toroidal planes).

The grid boundary is defined by the minimum and maximum values of the radial coordinate (rmin, rmax), the minimum and maximum values of the z-coordinate (zmin, zmax). Note the R and z grids include both end points while the phi dimension includes the start point and excludes the end point. For free boundary calculations, the grid should be large enough to contain the entire plasma. For example, rmin should be smaller than min R(theta, phi) where R is the radial coordinate of the plasma. The same applies for the z-coordinate.

The choice of the number or radial and vertical gridpoints is not as straightforward. The VMEC code uses these grids to ‘deform’ the plasma boundary in a iterative sense. The complication revolves around the spectral condensation VMEC performs in the poloidal direction. The code calculates the location of the poloidal grid points so that an optimized choice is made for the form of the plasma. The user should decide how accurately the wish to know the plasma boundary. If [cm] precision is required then the number of grid points chosen should provide at least this resolution. Remember that the more datapoints, the slower VMEC will run.

Currently, this function only supports one “current group”: a set of coils that are electrically connected in series. Using multiple current groups emmulates groups of coils that are connected to distinct power supplies, and allows for the current in each group can be controlled independently. For example, W7-X has 7 current groups, one for each base coil (5 nonplanar coils + 2 planar coils), which allows all planar coils to be turned on without changing the currents in the nonplanar coils. In free-boundary vmec, the currents in each current group are scaled using the “extcur” array in the input file. Since only one current group is supported by this function, for free-boundary vmec, the vmec.indata.extcur array should have a single nonzero element, set to 1.0.

Parameters:

• filename (str) – Name of the NetCDF file to save.

• nr (int) – Number of grid points in the radial direction.

• nphi (int) – Number of planes in the toroidal angle.

• nz (int) – Number of grid points in the z coordinate.

• rmin (float) – Minimum value of major radius coordinate of the grid.

• rmax (float) – Maximum value of major radius coordinate of the grid.

• zmin (float) – Minimum value of z-coordinate of the grid.

• zmax (float) – Maximum value of z-coordinate of the grid.

• nfp (int) – Number of field periods.

• include_potential (bool) – Boolean to include the vector potential A. Defaults to false.

to_vtk(filename, nr=10, nphi=10, nz=10, rmin=1.0, rmax=2.0, zmin=-0.5, zmax=0.5)

Export the field evaluated on a regular grid for visualisation with e.g. Paraview.

class simsopt.field.magneticfield.MagneticFieldMultiply(scalar, Bfield)

Bases: MagneticField

Class used to multiply a magnetic field by a scalar. It takes as input a MagneticField class and a scalar and multiplies B, A and their derivatives by that value.

as_dict(serial_objs_dict) → dict

A JSON serializable dict representation of an object.

classmethod from_dict(d, serial_objs_dict, recon_objs)

Parameters:

d – Dict representation.

Returns:

GSONable class.

class simsopt.field.magneticfield.MagneticFieldSum(Bfields)

Bases: MagneticField

Class used to sum two or more magnetic field together. It can either be called directly with a list of magnetic fields given as input and outputing another magnetic field with B, A and its derivatives added together or it can be called by summing magnetic fields classes as Bfield1 + Bfield1

B_vjp(v)

as_dict(serial_objs_dict) → dict

A JSON serializable dict representation of an object.

classmethod from_dict(d, serial_objs_dict, recon_objs)

Parameters:

d – Dict representation.

Returns:

GSONable class.

simsopt.field.magneticfieldclasses module

class simsopt.field.magneticfieldclasses.CircularCoil(r0=0.1, center=[0, 0, 0], I=159154.94309189534, normal=[0, 0])

Bases: MagneticField

Magnetic field created by a single circular coil evaluated using analytical functions, including complete elliptic integrals of the first and second kind. As inputs, it takes the radius of the coil (r0), its center, current (I) and its normal vector [either spherical angle components (normal=[theta,phi]) or (x,y,z) components of a vector (normal=[x,y,z])]). The (theta,phi) angles are related to the (x,y,z) components of the normal vector via theta = np.arctan2(normal[1], normal[0]) and phi = np.arctan2(np.sqrt(normal[0]**2+normal[1]**2), normal[2]). Sign convention: CircularCoil with a positive current produces a magnetic field vector in the same direction as the normal when evaluated at the center of the coil.a

Parameters:

• r0 – radius of the coil

• center – point at the coil center

• I – current of the coil in Ampere’s

• normal – if list with two values treats it as spherical angles theta and phi of the normal vector to the plane of the coil centered at the coil center, if list with three values treats it a vector

property I

as_dict(serial_objs_dict)

A JSON serializable dict representation of an object.

classmethod from_dict(d, serial_objs_dict, recon_objs)

Parameters:

d – Dict representation.

Returns:

GSONable class.

gamma(points=64)

Export points of the coil.

get_dofs()

num_dofs()

set_dofs_impl(dofs)

to_vtk(filename, close=False)

Export circular coil to VTK format

Parameters:

• filename – Name of the file to write.

• close – Whether to draw the segment from the last quadrature point back to the first.

class simsopt.field.magneticfieldclasses.DipoleField(dipole_grid, dipole_vectors, stellsym=True, nfp=1, coordinate_flag='cartesian', m_maxima=None, R0=1)

Bases: MagneticField

Computes the MagneticField induced by N dipoles. The field is given by

\[B(\mathbf{x}) = \frac{\mu_0}{4\pi} \sum_{i=1}^{N} \left(\frac{3\mathbf{r}_i\cdot \mathbf{m}_i}{|\mathbf{r}_i|^5}\mathbf{r}_i - \frac{\mathbf{m}_i}{|\mathbf{r}_i|^3}\right)\]

where \(\mu_0=4\pi\times 10^{-7}\;N/A^2\) is the permeability of free space and \(\mathbf{r_i} = \mathbf{x} - \mathbf{x}^{dipole}_i\) is the vector between the field evaluation point and the dipole \(i\) position.

Parameters:

• dipole_grid – 2D numpy array, shape (ndipoles, 3). A set of points corresponding to the locations of magnetic dipoles.

• dipole_vectors – 2D numpy array, shape (ndipoles, 3). The dipole vectors of each of the dipoles in the grid.

• stellsym – bool (default True). Whether or not the dipole grid is stellarator symmetric.

• nfp – int (default 1). The field-period symmetry of the dipole-grid.

• coordinate_flag – string (default “cartesian”). The global coordinate system that should be considered grid-aligned in the calculation. The options are “cartesian” (rectangular bricks), “cylindrical” (cylindrical bricks), and “toroidal” (uniform grid in simple toroidal coordinates). Note that this ASSUMES that the global coordinate system for the dipole locations is one of these three choices, so be careful if your permanent magnets are shaped/arranged differently!

• m_maxima – 1D numpy array, shape (ndipoles,). The maximum dipole strengths of each magnet in the grid. If not specified, defaults to using the largest dipole strength of the magnets in dipole_grid, and using this value for all the dipoles. Needed for plotting normalized dipole magnitudes in the vtk functionality.

• R0 – double. The value of the major radius of the stellarator needed only for simple toroidal coordinates.

_dipole_fields_from_symmetries(dipole_grid, dipole_vectors, stellsym=True, nfp=1, coordinate_flag='cartesian', m_maxima=None, R0=1)

Takes the dipoles and grid initialized in a PermanentMagnetOptimizer (for a half-period surface) and generates the full dipole manifold so that the call to B() (the magnetic field from the dipoles) correctly returns contributions from all the dipoles from symmetries.

_toVTK(vtkname)

Write dipole data into a VTK file (acknowledgements to Caoxiang’s CoilPy code).

Parameters:

vtkname (str) – VTK filename, will be appended with .vts or .vtu.

class simsopt.field.magneticfieldclasses.Dommaschk(mn=[[0, 0]], coeffs=[[0, 0]])

Bases: MagneticField

Vacuum magnetic field created by an explicit representation of the magnetic field scalar potential as proposed by W. Dommaschk (1986), Computer Physics Communications 40, 203-218. As inputs, it takes the arrays for the harmonics m, n and its corresponding coefficients.

Parameters:

• m – first harmonic array

• n – second harmonic array

• coeffs – coefficient for Vml for each of the ith index of the harmonics m and n

as_dict(serial_objs_dict) → dict

A JSON serializable dict representation of an object.

classmethod from_dict(d, serial_objs_dict, recon_objs)

Parameters:

d – Dict representation.

Returns:

GSONable class.

property mn

class simsopt.field.magneticfieldclasses.InterpolatedField(field, degree, rrange, phirange, zrange, extrapolate=True, nfp=1, stellsym=False, skip=None)

Bases: InterpolatedField, MagneticField

This field takes an existing field and interpolates it on a regular grid in \(r,\phi,z\). This resulting interpolant can then be evaluated very quickly.

__init__(field, degree, rrange, phirange, zrange, extrapolate=True, nfp=1, stellsym=False, skip=None)

Parameters:

• field – the underlying simsopt.field.magneticfield.MagneticField to be interpolated.

• degree – the degree of the piecewise polynomial interpolant.

• rrange – a 3-tuple of the form (rmin, rmax, nr). This mean that the interval \([rmin, rmax]\) is split into nr many subintervals.

• phirange – a 3-tuple of the form (phimin, phimax, nphi).

• zrange – a 3-tuple of the form (zmin, zmax, nz).

• extrapolate – whether to extrapolate the field when evaluate outside the integration domain or to throw an error.

• nfp – Whether to exploit rotational symmetry. In this case any angle is always mapped into the interval \([0, 2\pi/\mathrm{nfp})\), hence it makes sense to use phimin=0 and phimax=2*np.pi/nfp.

• stellsym – Whether to exploit stellarator symmetry. In this case z is always mapped to be positive, hence it makes sense to use zmin=0.

• skip –

  a function that takes in a point (in cylindrical (r,phi,z) coordinates) and returns whether to skip that location when building the interpolant or not. The signature should be

  def skip(r: double, phi: double, z: double) -> bool:
      ...
  

  See also here https://github.com/hiddenSymmetries/simsopt/pull/227 for a graphical illustration.

to_vtk(filename)

Export the field evaluated on a regular grid for visualisation with e.g. Paraview.

class simsopt.field.magneticfieldclasses.MirrorModel(B0=6.51292, gamma=0.124904, Z_m=0.98)

Bases: MagneticField

Model magnetic field employed in https://arxiv.org/abs/2305.06372 to study the magnetic mirror experiment WHAM. The magnetic field is given by \(\vec{B}=B_R \vec{e}_R + B_Z \vec{e}_Z\), where \(\vec{e}_R\) and \(\vec{e}_Z\) are the cylindrical radial and axial unit vectors, respectively, and \(B_R\) and \(B_Z\) are given by

\[B_R = -\frac{1}{R} \frac{\partial\psi}{\partial Z}, \; B_Z = \frac{1}{R}. \frac{\partial\psi}{\partial R}\]

In this model, the magnetic flux function \(\psi\) is written as a double Lorentzian function

\[\psi = \frac{R^2 \mathcal{B}}{2 \pi \gamma}\left(\left[1+\left(\frac{Z-Z_m}{\gamma}\right)^2\right]^{-1}+\left[1+\left(\frac{Z+Z_m}{\gamma}\right)^2\right]^{-1}\right).\]

Note that this field is neither a vacuum field nor a solution of MHD force balance. The input parameters are B0, gamma and Z_m with the standard values the ones used in https://arxiv.org/abs/2305.06372, that is, B0 = 6.51292, gamma = 0.124904, and Z_m = 0.98.

Parameters:

• B0 – parameter \(\mathcal{B}\) of the flux surface function

• gamma – parameter \(\gamma\) of the flux surface function

• Z_m – parameter \(Z_m\) of the flux surface function

as_dict(serial_objs_dict) → dict

A JSON serializable dict representation of an object.

classmethod from_dict(d, serial_objs_dict, recon_objs)

Parameters:

d – Dict representation.

Returns:

GSONable class.

class simsopt.field.magneticfieldclasses.PoloidalField(R0, B0, q)

Bases: MagneticField

Magnetic field purely in the poloidal direction, that is, in the theta direction of a poloidal-toroidal coordinate system. Its modulus is given by B = B0 * r / (R0 * q) so that, together with the toroidal field, it creates a safety factor equals to q

Parameters:

• B0 – modulus of the magnetic field at R0

• R0 – major radius of the magnetic axis

• q – safety factor/pitch angle of the magnetic field lines

as_dict(serial_objs_dict) → dict

A JSON serializable dict representation of an object.

classmethod from_dict(d, serial_objs_dict, recon_objs)

Parameters:

d – Dict representation.

Returns:

GSONable class.

class simsopt.field.magneticfieldclasses.Reiman(iota0=0.15, iota1=0.38, k=[6], epsilonk=[0.01], m0=1)

Bases: MagneticField

Magnetic field model in section 5 of Reiman and Greenside, Computer Physics Communications 43 (1986) 157—167. This field allows for an analytical expression of the magnetic island width that can be used for island optimization. However, the field is not completely physical as it does not have nested flux surfaces.

Parameters:

• iota0 – unperturbed rotational transform

• iota1 – unperturbed global magnetic shear

• k – integer array specifying the Fourier modes used

• epsilonk – coefficient of the Fourier modes

• m0 – toroidal symmetry parameter (normally m0=1)

as_dict(serial_objs_dict)

A JSON serializable dict representation of an object.

classmethod from_dict(d, serial_objs_dict, recon_objs)

Parameters:

d – Dict representation.

Returns:

GSONable class.

class simsopt.field.magneticfieldclasses.ScalarPotentialRZMagneticField(phi_str)

Bases: MagneticField

Vacuum magnetic field as a solution of B = grad(Phi) where Phi is the magnetic field scalar potential. It takes Phi as an input string, which should contain an expression involving the standard cylindrical coordinates (R, phi, Z) Example: ScalarPotentialRZMagneticField(“2*phi”) yields a magnetic field B = grad(2*phi) = (0,2/R,0). In order for the analytical derivatives to be performed by sympy, a term 1e-30*Phi*R*Z is added to every entry. Note: this function needs sympy.

Parameters:

phi_str – string containing vacuum scalar potential expression as a function of R, Z and phi

as_dict(serial_objs_dict) → dict

A JSON serializable dict representation of an object.

classmethod from_dict(d, serial_objs_dict, recon_objs)

Parameters:

d – Dict representation.

Returns:

GSONable class.

class simsopt.field.magneticfieldclasses.ToroidalField(R0, B0)

Bases: MagneticField

Magnetic field purely in the toroidal direction, that is, in the phi direction with (R,phi,Z) the standard cylindrical coordinates. Its modulus is given by B = B0*R0/R where R0 is the first input and B0 the second input to the function.

Parameters:

• B0 – modulus of the magnetic field at R0

• R0 – radius of normalization

as_dict(serial_objs_dict) → dict

A JSON serializable dict representation of an object.

classmethod from_dict(d, serial_objs_dict, recon_objs)

Parameters:

d – Dict representation.

Returns:

GSONable class.

simsopt.field.mgrid module

class simsopt.field.mgrid.MGrid(nr: int = 51, nz: int = 51, nphi: int = 24, nfp: int = 2, rmin: float = 0.2, rmax: float = 0.4, zmin: float = -0.1, zmax: float = 0.1)

Bases: object

This class reads and writes mgrid (NetCDF) files for use in free boundary VMEC and other codes.

An mgrid file contains a grid in cylindrical coordinates (R, phi, z) upon which the magnetic field is evaluated. The grid is defined by the number of points in the radius (nr), the number of planes per field period in the toroidal angle (nphi), and the number of points in the z-coordinate (nz). A good rule here is to choose at least 4 times as many toroidal planes as the maximum toroidal mode number (ex. if ntor=6 then you should have at least 24 toroidal planes).

The grid boundary is defined by the minimum and maximum values of the radial coordinate (rmin, rmax), the minimum and maximum values of the z-coordinate (zmin, zmax). Note the R and z grids include both end points while the phi dimension includes the start point and excludes the end point. For free boundary calculations, the grid should be large enough to contain the entire plasma. For example, rmin should be smaller than min R(theta, phi) where R is the radial coordinate of the plasma. The same applies for the z-coordinate.

The choice of the number or radial and vertical gridpoints is not as straightforward as choosing nphi. The VMEC code uses these grids to “deform” the plasma boundary in a iterative sense. The complication revolves around the spectral condensation VMEC performs in the poloidal direction. The code calculates the location of the poloidal grid points so that an optimized choice is made for the form of the plasma. The user should decide how accurately the wish to know the plasma boundary. If centimeter precision is required then the number of grid points chosen should provide at least this resolution. Remember that the more datapoints, the slower VMEC will run.

Parameters:

• nr (int) – Number of grid points in the radial direction. Default is 51.

• nz (int) – Number of grid points in the z coordinate. Default is 51.

• nphi (int) – Number of planes in the toroidal angle. Default is 24.

• nfp (int) – Number of field periods. Default is 2.

• rmin (float) – Minimum value of major radius coordinate of the grid. Default is 0.2.

• rmax (float) – Maximum value of major radius coordinate of the grid. Default is 0.4.

• zmin (float) – Minimum value of z-coordinate of the grid. Default is -0.1.

• zmax (float) – Maximum value of z-coordinate of the grid. Default is 0.1.

add_field_cylindrical(br, bp, bz, ar=None, ap=None, az=None, name=None)

This function saves the magnetic field \(B\), and (optionally) the vector potential \(A\), to the MGrid object. \(B\) and \(A\) are provided on a tensor product grid in cylindrical components \((R, \phi, z)\).

This function may be called once for each current group, to save groups of fields that can be scaled using the vmec.indata.extcur array. Current groups emmulate groups of coils that are connected to distinct power supplies, and allows for the current in each group can be controlled independently. For example, W7-X has 7 current groups, one for each base coil (5 nonplanar coils + 2 planar coils), which allows all planar coils to be turned on without changing the currents in the nonplanar coils. In free-boundary vmec, the currents in each current group are scaled using the “extcur” array in the input file.

Parameters:

• br (ndarray) – (nphi, nz, nr) array of the radial component of B-field.

• bp (ndarray) – (nphi, nz, nr) array of the azimuthal component of B-field.

• bz (ndarray) – (nphi, nz, nr) array of the z-component of B-field.

• ar (ndarray, Optional) – (nphi, nz, nr) array of the radial component of the vector potential A. Default is None.

• ap (ndarray, Optional) – (nphi, nz, nr) array of the azimuthal component of the vector potential A. Default is None.

• az (ndarray, Optional) – (nphi, nz, nr) array of the axial component of the vector potential A. Default is None.

• name (str, Optional) – Name of the coil group. Default is None.

classmethod from_file(filename)

This method reads MGrid data from file.

Parameters:

filename (str) – mgrid netCDF input file name.

Returns:

MGrid object with data from file.

Return type:

MGrid

plot(jphi=0, bscale=0, show=True)

Creates a plot of the mgrid data. Shows the three components (br,bphi,bz) in a 2D plot for a fixed toroidal plane.

Parameters:

• jphi (int) – integer index for a toroidal slice. Default is 0.

• bscale (float) – sets saturation scale for colorbar. (This is useful, because the mgrid domain often includes coil currents, and arbitrarily close to the current source the Bfield values become singular.). Default is 0.

• show (bool) – If True, will call matplotlib’s show() function. Default is True.

Returns:

2-element tuple (fig, axs) with the Matplotlib figure and axes handles.

write(filename)

Export class data as a netCDF binary.

Parameters:

filename (str) – output file name.

simsopt.field.mgrid._pad_string(string)

Pads a string with 30 underscores (for writing coil group names).

simsopt.field.mgrid._unpack(binary_array)

Decrypt binary char array into a string. This function is used for reading coil group names.

simsopt.field.normal_field module

class simsopt.field.normal_field.CoilNormalField(coilset: CoilSet = None)

Bases: NormalField

A SPEC NormalField generated by a CoilSet.

The CoilNormalField provides the same interface as the NormalField, but its degrees of freedom are inherited from its CoilSet parent.

Parameters:

coilset – The CoilSet object from which to inherit the degrees of freedom

Properties:

surface: The computational boundary of the SPEC simulation, that is managed by the CoilSet. vns/vnc: fourier harmonics of the normal field. This property is cached, and recomputed only when the parents’ DOFS (the coils) change.

change_resolution(*args, **kwargs)

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.

property coilset

fixed_range(*args, **kwargs)

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

In case of non stellarator symmetric field, both vns and vnc will be fixed (or unfixed)

classmethod from_spec(*args, **kwargs)

Initialize using the harmonics in SPEC input file

classmethod from_spec_object(*args, **kwargs)

Initialize using the simsopt SPEC object’s attributes

get_dofs(*args, **kwargs)

get DOFs from vns and vnc

get_index_in_array(m, n, mpol=None, ntor=None)

Returns position of mode (m,n) in vns or vnc array

Args: - m: poloidal mode number - n: toroidal mode number (normalized by Nfp) - mpol: resolution of dofs array. If None (by default), use self.mpol - ntor: resolution of dofs array. If None (by default), use self.ntor - even: set to True to get vnc. Default is False

get_index_in_dofs(*args, **kwargs)

Returns position of mode (m,n) in dofs array

Args: - m: poloidal mode number - n: toroidal mode number (normalized by Nfp) - mpol: resolution of dofs array. If None (by default), use self.mpol - ntor: resolution of dofs array. If None (by default), use self.ntor - even: set to True to get vnc. Default is False

get_vnc(m, n)

get_vnc_asarray()

Return the vnc as a single array

get_vns(m, n)

get_vns_asarray()

Return the vns as a single array

get_vns_vnc_asarray()

Return the vns and vnc as two arrays single array

optimize_coils(targetvns, targetvnc=None, TARGET_LENGTH=1000, MAXITER=300)

Simple convenience function to optimize the coils to match the target vns and vnc using a FOCUS-style algorithm.

Uses the simplest FOCUS optimization consisting of only the quadratic flux penalty and a length penalty.

Parameters:

• targetvns – The target odd fourier modes of \(\mathbf{B}\cdot\mathbf{\vec{n}}\). 2D array of size (mpol+1)x(2ntor+1).

• targetvnc – The target even fourier modes of \(\mathbf{B}\cdot\mathbf{\vec{n}}\). 2D array of size (mpol+1)x(2ntor+1). Ignored if stellsym if True.

• TARGET_LENGTH – The target length of the coils. Default is 1000.

• MAXITER – The maximum number of iterations. Default is 1000.

Returns:

the result object from scipy.optimize.minimize

Return type:

res

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.

reduce_coilset(nsv='nonzero')

Replace the coilset with a Re:w ducedCoilSet keeping the first nsv singular values.

Note: Should this be done by proc0 and broadcast? Is SVD deterministic?

set_vnc(*args, **kwargs)

set_vnc_asarray(*args, **kwargs)

Set the vnc from a single array

set_vns(*args, **kwargs)

set_vns_asarray(*args, **kwargs)

Set the vns from a single array

set_vns_vnc_asarray(*args, **kwargs)

Set the vns and vnc from two single arrays

property surface

property vnc

property vns

class simsopt.field.normal_field.NormalField(nfp=1, stellsym=True, mpol=1, ntor=0, vns=None, vnc=None, surface=None)

Bases: Optimizable

NormalField represents the magnetic field normal to a toroidal surface, for example the computational boundary of SPEC free-boundary.

It consists a surface (the computational boundary), and a set of Fourier harmonics that describe the normal component of an externally provided field.

The Fourier harmonics are the degees of freedom, the computational boundary is kept fixed. The Fourier harmonics are stored in the SPEC convention, i.e. it is the fourier components of B.(partialvec{r}/ partialtheta times partialvec{r}/ partialzeta) on the surface with a 1/(2pi)^2 Fourier normalization factor.

Parameters:

• nfp – The number of field period

• stellsym – Whether (=True) or not (=False) stellarator symmetry is enforced.

• mpol – Poloidal Fourier resolution

• ntor – Toroidal Fourier resolution

• vns –

  Odd fourier modes of \(\mathbf{B}\cdot\mathbf{\hat{n}}\). 2D array of size (mpol+1)x(2ntor+1). Set to None to fill with zeros

  vns( mm, self.ntor+nn ) is the mode (mm,nn)

• vnc –

  Even fourier modes of \(\mathbf{B}\cdot\mathbf{\hat{n}}\). 2D array of size (mpol+1)x(2ntor+1). Ignored if stellsym if True. Set to None to fill with zeros

  vnc( mm, self.ntor+nn ) is the mode (mm,nn)

_make_names()

Form a list of names of the vns, vnc

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.

check_mn(m, n)

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

In case of non stellarator symmetric field, both vns and vnc will be fixed (or unfixed)

classmethod from_spec(filename)

Initialize using the harmonics in SPEC input file

classmethod from_spec_object(spec)

Initialize using the simsopt SPEC object’s attributes

get_dofs()

get DOFs from vns and vnc

get_index_in_array(m, n, mpol=None, ntor=None)

Returns position of mode (m,n) in vns or vnc array

Args: - m: poloidal mode number - n: toroidal mode number (normalized by Nfp) - mpol: resolution of dofs array. If None (by default), use self.mpol - ntor: resolution of dofs array. If None (by default), use self.ntor - even: set to True to get vnc. Default is False

get_index_in_dofs(m, n, mpol=None, ntor=None, even=False)

Returns position of mode (m,n) in dofs array

Args: - m: poloidal mode number - n: toroidal mode number (normalized by Nfp) - mpol: resolution of dofs array. If None (by default), use self.mpol - ntor: resolution of dofs array. If None (by default), use self.ntor - even: set to True to get vnc. Default is False

get_real_space_field()

Fourier transform the field and get the real-space values of the normal component of the externally provided field on the computational boundary. The returned array will be of size specified by the surfaces’ quadpoints and located on the quadpoints.

get_vnc(m, n)

get_vnc_asarray(mpol=None, ntor=None)

Return the vnc as a single array

get_vns(m, n)

get_vns_asarray(mpol=None, ntor=None)

Return the vns as a single array

get_vns_vnc_asarray(mpol=None, ntor=None)

Return the vns and vnc as two arrays single array

set_vnc(m, n, value)

set_vnc_asarray(vnc, mpol=None, ntor=None)

Set the vnc from a single array

set_vns(m, n, value)

set_vns_asarray(vns, mpol=None, ntor=None)

Set the vns from a single array

set_vns_vnc_asarray(vns, vnc, mpol=None, ntor=None)

Set the vns and vnc from two single arrays

property vnc

property vns

simsopt.field.sampling module

simsopt.field.sampling.draw_uniform_on_curve(curve, nsamples, safetyfactor=10)

Uses rejection sampling to sample points on a curve. Warning: assumes that the underlying quadrature points on the Curve are uniformly distributed.

Parameters:

• curve – The simsopt.geo.curve.Curve to spawn the particles on.

• nsamples – number of samples.

• safetyfactor – how many more samples than nsamples to generate for rejection/acceptance.

simsopt.field.sampling.draw_uniform_on_surface(surface, nsamples, safetyfactor=10)

Uses rejection sampling to sample points on a surface. Warning: assumes that the underlying quadrature points on the surface are uniformly distributed.

Parameters:

• surface – The simsopt.geo.surface.Surface to spawn the particles on.

• nsamples – number of samples.

• safetyfactor – how many more samples than nsamples to generate for rejection/acceptance.

simsopt.field.selffield module

This module contains functions for computing the self-field of a coil using the methods from:

Hurwitz, Siena, Matt Landreman, and Thomas M. Antonsen. “Efficient calculation of the self magnetic field, self-force, and self-inductance for electromagnetic coils.” IEEE Transactions on Magnetics (2024).

Landreman, Matt, Siena Hurwitz, and Thomas M. Antonsen. “Efficient calculation of self magnetic field, self-force, and self-inductance for electromagnetic coils with rectangular cross-section.” Nuclear Fusion 65.3 (2025): 036008.

simsopt.field.selffield.B_regularized_pure(gamma, gammadash, gammadashdash, quadpoints, current, regularization)

Compute the regularized field on a coil following the Landreman and Hurwitz method

Parameters:

• gamma (array (shape (n,3))) – The curve of the coil.

• gammadash (array (shape (n,3))) – The first derivative of the curve.

• gammadashdash (array (shape (n,3))) – The second derivative of the curve.

• quadpoints (array (shape (n,))) – The quadrature points of the curve.

• current (float) – The current in the coil.

• regularization (float) – The regularization parameter.

• simsopt (The factors of 2π in the next few lines come from the fact that)

• 2π. (uses a curve parameter that goes up to 1 rather than)

Returns:

The regularized field on the coil.

Return type:

array (shape (n,3))

simsopt.field.selffield._rectangular_xsection_delta(a, b)

Auxiliary function for regularization in rectangular conductor.

\[\delta = \exp \left( - \frac{25}{6} + K \right)\]

where K is the auxiliary function defined above.

Parameters:

• a (float) – The width of the rectangular conductor.

• b (float) – The height of the rectangular conductor.

Returns:

The regularization parameter.

Return type:

float

simsopt.field.selffield._rectangular_xsection_k(a, b)

Auxiliary function for regularization in rectangular conductor.

\[k = \frac{4 b}{3 a} \arctan \left( \frac a b \right) + \frac{4 a}{3 b} \arctan \left( \frac b a \right) + \frac{b^2}{6 a^2} \log \left( \frac b a \right) + \frac{a^2}{6 b^2} \log \left( \frac a b \right) - \frac{a^4 - 6 a^2 b^2 + b^4}{6 a^2 b^2} \log \left( \frac a b + \frac b a \right)\]

where a is the width of the rectangular conductor and b is the height.

Parameters:

• a (float) – The width of the rectangular conductor.

• b (float) – The height of the rectangular conductor.

Returns:

The regularization parameter.

Return type:

float

simsopt.field.selffield.regularization_circ(a)

Regularization for a circular conductor.

\[\delta = a^2 / \sqrt{e}\]

where e = 2.718… is the base of the natural logarithm and a is the radius of the circular conductor.

Parameters:

a (float) – The radius of the circular conductor.

Returns:

The regularization parameter.

Return type:

float

simsopt.field.selffield.regularization_rect(a, b)

Regularization for a rectangular conductor.

\[\delta = a b \exp \left( - \frac{25}{6} + K \right)\]

where K is the auxiliary function defined above, a is the width of the rectangular conductor and b is the height.

Parameters:

• a (float) – The width of the rectangular conductor.

• b (float) – The height of the rectangular conductor.

Returns:

The regularization parameter.

Return type:

float

simsopt.field.tracing module

class simsopt.field.tracing.IterationStoppingCriterion

Bases: IterationStoppingCriterion

Stop the iteration once the maximum number of iterations is reached.

class simsopt.field.tracing.LevelsetStoppingCriterion(classifier)

Bases: LevelsetStoppingCriterion

Based on a scalar function \(f:R^3\to R\), this criterion checks whether \(f(x, y, z) < 0\) and stops the iteration once this is true.

The idea is to use this for example with signed distance functions to a surface.

class simsopt.field.tracing.MaxRStoppingCriterion

Bases: MaxRStoppingCriterion

Stop the iteration once a particle goes above a critical value of R, the radial cylindrical coordinate.

Usage:

stopping_criteria=[MaxRStopingCriterion(crit_r)]

where crit_r is the value of the critical coordinate.

class simsopt.field.tracing.MaxToroidalFluxStoppingCriterion

Bases: MaxToroidalFluxStoppingCriterion

Stop the iteration once a particle falls above a critical value of s, the normalized toroidal flux. This should only be used when tracing trajectories in a flux coordinate system (i.e., trace_particles_boozer).

Usage:

stopping_criteria=[MaxToroidalFluxStopingCriterion(s)]

where s is the value of the maximum normalized toroidal flux.

class simsopt.field.tracing.MaxZStoppingCriterion

Bases: MaxZStoppingCriterion

Stop the iteration once a particle gove above a critical value of Z, the cylindrical vertical coordinate.

Usage:

stopping_criteria=[MaxZStopingCriterion(crit_z)]

where crit_z is the value of the critical coordinate.

class simsopt.field.tracing.MinRStoppingCriterion

Bases: MinRStoppingCriterion

Stop the iteration once a particle falls below a critical value of R, the radial cylindrical coordinate.

Usage:

stopping_criteria=[MinRStopingCriterion(crit_r)]

where crit_r is the value of the critical coordinate.

class simsopt.field.tracing.MinToroidalFluxStoppingCriterion

Bases: MinToroidalFluxStoppingCriterion

Stop the iteration once a particle falls below a critical value of s, the normalized toroidal flux. This StoppingCriterion is important to use when tracing particles in flux coordinates, as the poloidal angle becomes ill-defined at the magnetic axis. This should only be used when tracing trajectories in a flux coordinate system (i.e., trace_particles_boozer).

Usage:

stopping_criteria=[MinToroidalFluxStopingCriterion(s)]

where s is the value of the minimum normalized toroidal flux.

class simsopt.field.tracing.MinZStoppingCriterion

Bases: MinZStoppingCriterion

Stop the iteration once a particle falls below a critical value of Z, the cylindrical vertical coordinate.

Usage:

stopping_criteria=[MinZStopingCriterion(crit_z)]

where crit_z is the value of the critical coordinate.

class simsopt.field.tracing.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.field.tracing.ToroidalTransitStoppingCriterion

Bases: ToroidalTransitStoppingCriterion

Stop the iteration once the maximum number of toroidal transits is reached.

Usage:

stopping_criteria=[ToroidalTransitStoppingCriterion(ntransits,flux)]

where ntransits is the maximum number of toroidal transits and flux is a boolean indicating whether tracing is being performed in a flux coordinate system.

simsopt.field.tracing.compute_fieldlines(field, R0, Z0, tmax=200, tol=1e-07, phis=[], stopping_criteria=[], comm=None)

Compute magnetic field lines by solving

\[[\dot x, \dot y, \dot z] = B(x, y, z)\]

Integration is initialized on the \(\phi = 0\) plane.

Parameters:

• field – the magnetic field \(B\)

• R0 – list of radial components of initial points

• Z0 – list of vertical components of initial points

• tmax – for how long to trace. will do roughly |B|*tmax/(2*pi*r0) revolutions of the device

• tol – tolerance for the adaptive ode solver

• phis – list of angles in [0, 2pi] for which intersection with the plane corresponding to that phi should be computed

• stopping_criteria – list of stopping criteria, mostly used in combination with the LevelsetStoppingCriterion accessed via simsopt.field.tracing.SurfaceClassifier.

Returns: 2 element tuple containing

• res_tys:

  A list of numpy arrays (one for each particle) describing the solution over time. The numpy array is of shape (ntimesteps, 4). Each row contains the time and the position, i.e.`[t, x, y, z]`.

• res_phi_hits:

  A list of numpy arrays (one for each particle) containing information on each time the particle hits one of the phi planes or one of the stopping criteria. Each row of the array contains [time, idx, x, y, z], where idx tells us which of the phis or stopping_criteria was hit. If idx>=0, then phis[int(idx)] was hit. If idx<0, then stopping_criteria[int(-idx)-1] was hit.

simsopt.field.tracing.compute_poloidal_transits(res_tys, ma=None, flux=True)

Computes the number of poloidal transits of an orbit. For the case of particles traced in a MagneticField (not a BoozerMagneticField), the poloidal angle is computed using the arctangent angle in the poloidal plane with respect to the coordinate axis, ma,

\[\theta = \tan^{-1} \left( \frac{R(\phi)-R_{\mathrm{ma}}(\phi)}{Z(\phi)-Z_{\mathrm{ma}}(\phi)} \right),\]

where \((R,\phi,Z)\) are the cylindrical coordinates of the trajectory and \((R_{\mathrm{ma}}(\phi),Z_{\mathrm{ma}(\phi)})\) is the position of the coordinate axis.

Parameters:

• res_tys – trajectory solution computed from trace_particles() or trace_particles_boozer() with forget_exact_path=False.

• ma – an instance of Curve representing the coordinate axis with respect to which the poloidal angle is computed. If orbit is computed in Boozer coordinates, ma should be None.

• flux – if True, res_tys represents the position in flux coordinates (should be True if computed from trace_particles_boozer()). If True, ma is not used.

Returns:

array with length len(res_tys). Each element contains the

number of poloidal transits of the orbit.

Return type:

ntransits

simsopt.field.tracing.compute_resonances(res_tys, res_phi_hits, ma=None, delta=0.01)

Computes resonant particle orbits given the output of either trace_particles() or trace_particles_boozer(), res_tys and res_phi_hits/res_zeta_hits, with forget_exact_path=False. Resonance indicates a trajectory which returns to the same position at the \(\zeta = 0\) plane after mpol poloidal turns and ntor toroidal turns. For the case of particles traced in a MagneticField (not a BoozerMagneticField), the poloidal angle is computed using the arctangent angle in the poloidal plane with respect to the coordinate axis, ma,

\[\theta = \tan^{-1} \left( \frac{R(\phi)-R_{\mathrm{ma}}(\phi)}{Z(\phi)-Z_{\mathrm{ma}}(\phi)} \right),\]

where \((R,\phi,Z)\) are the cylindrical coordinates of the trajectory and \((R_{\mathrm{ma}}(\phi),Z_{\mathrm{ma}(\phi)})\) is the position of the coordinate axis.

Parameters:

• res_tys – trajectory solution computed from trace_particles() or trace_particles_boozer() with forget_exact_path=False

• res_phi_hits – output of trace_particles() or trace_particles_boozer() with phis/zetas = [0]

• ma – an instance of Curve representing the coordinate axis with respect to which the poloidal angle is computed. If the orbit is computed in flux coordinates, ma should be None. (defaults to None)

• delta – the distance tolerance in the poloidal plane used to compute a resonant orbit. (defaults to 1e-2)

Returns:

list of 7d arrays containing resonant particle orbits. The

elements of each array is [s0, theta0, zeta0, vpar0, t, mpol, ntor] if ma=None, and [R0, Z0, phi0, vpar0, t, mpol, ntor] otherwise. Here (s0, theta0, zeta0, vpar0)/(R0, Z0, phi0, vpar0) indicates the initial position and parallel velocity of the particle, t indicates the time of the resonance, mpol is the number of poloidal turns of the orbit, and ntor is the number of toroidal turns.

Return type:

resonances

simsopt.field.tracing.compute_toroidal_transits(res_tys, flux=True)

Computes the number of toroidal transits of an orbit.

Parameters:

• res_tys – trajectory solution computed from trace_particles() or trace_particles_boozer() with forget_exact_path=False.

• flux – if True, res_tys represents the position in flux coordinates (should be True if computed from trace_particles_boozer())

Returns:

array with length len(res_tys). Each element contains the

number of toroidal transits of the orbit.

Return type:

ntransits

simsopt.field.tracing.particles_to_vtk(res_tys, filename)

Export particle tracing or field lines to a vtk file. Expects that the xyz positions can be obtained by xyz[:, 1:4].

simsopt.field.tracing.plot_poincare_data(fieldlines_phi_hits, phis, filename, mark_lost=False, aspect='equal', dpi=300, xlims=None, ylims=None, surf=None, s=2, marker='o')

Create a poincare plot. Usage:

phis = np.linspace(0, 2*np.pi/nfp, nphis, endpoint=False)
res_tys, res_phi_hits = compute_fieldlines(
    bsh, R0, Z0, tmax=1000, phis=phis, stopping_criteria=[])
plot_poincare_data(res_phi_hits, phis, '/tmp/fieldlines.png')

Requires matplotlib to be installed.

simsopt.field.tracing.trace_particles(field: MagneticField, xyz_inits: Sequence[Real] | ndarray[tuple[Any, ...], dtype[float64]], parallel_speeds: Sequence[Real] | ndarray[tuple[Any, ...], dtype[float64]], tmax=0.0001, mass=6.69509884346e-27, charge=3.204353268e-19, Ekin=5.639661751679999e-13, tol=1e-09, comm=None, phis=[], stopping_criteria=[], mode='gc_vac', forget_exact_path=False, phase_angle=0)

Follow particles in a magnetic field.

In the case of mod='full' we solve

\[[\ddot x, \ddot y, \ddot z] = \frac{q}{m} [\dot x, \dot y, \dot z] \times B\]

in the case of mod='gc_vac' we solve the guiding center equations under the assumption \(\nabla p=0\), that is

\[\begin{split}[\dot x, \dot y, \dot z] &= v_{||}\frac{B}{|B|} + \frac{m}{q|B|^3} (0.5v_\perp^2 + v_{||}^2) B\times \nabla(|B|)\\ \dot v_{||} &= -\mu (B \cdot \nabla(|B|))\end{split}\]

where \(v_\perp = 2\mu|B|\). See equations (12) and (13) of [Guiding Center Motion, H.J. de Blank, https://doi.org/10.13182/FST04-A468].

Parameters:

• field – The magnetic field \(B\).

• xyz_inits – A (nparticles, 3) array with the initial positions of the particles.

• parallel_speeds – A (nparticles, ) array containing the speed in direction of the B field for each particle.

• tmax – integration time

• mass – particle mass in kg, defaults to the mass of an alpha particle

• charge – charge in Coulomb, defaults to the charge of an alpha particle

• Ekin – kinetic energy in Joule, defaults to 3.52MeV

• tol – tolerance for the adaptive ode solver

• comm – MPI communicator to parallelize over

• phis – list of angles in [0, 2pi] for which intersection with the plane corresponding to that phi should be computed

• stopping_criteria – list of stopping criteria, mostly used in combination with the LevelsetStoppingCriterion accessed via simsopt.field.tracing.SurfaceClassifier.

• mode – how to trace the particles. options are gc: general guiding center equations, gc_vac: simplified guiding center equations for the case \(\nabla p=0\), full: full orbit calculation (slow!)

• forget_exact_path – return only the first and last position of each particle for the res_tys. To be used when only res_phi_hits is of interest or one wants to reduce memory usage.

• phase_angle – the phase angle to use in the case of full orbit calculations

Returns: 2 element tuple containing

• res_tys:

  A list of numpy arrays (one for each particle) describing the solution over time. The numpy array is of shape (ntimesteps, M) with M depending on the mode. Each row contains the time and the state. So for mode=’gc’ and mode=’gc_vac’ the state consists of the xyz position and the parallel speed, hence each row contains [t, x, y, z, v_par]. For mode=’full’, the state consists of position and velocity vector, i.e. each row contains [t, x, y, z, vx, vy, vz].

• res_phi_hits:

  A list of numpy arrays (one for each particle) containing information on each time the particle hits one of the phi planes or one of the stopping criteria. Each row of the array contains [time] + [idx] + state, where idx tells us which of the phis or stopping_criteria was hit. If idx>=0, then phis[int(idx)] was hit. If idx<0, then stopping_criteria[int(-idx)-1] was hit.

simsopt.field.tracing.trace_particles_boozer(field: BoozerMagneticField, stz_inits: Sequence[Real] | ndarray[tuple[Any, ...], dtype[float64]], parallel_speeds: Sequence[Real] | ndarray[tuple[Any, ...], dtype[float64]], tmax=0.0001, mass=6.69509884346e-27, charge=3.204353268e-19, Ekin=5.639661751679999e-13, tol=1e-09, comm=None, zetas=[], stopping_criteria=[], mode='gc_vac', forget_exact_path=False)

Follow particles in a BoozerMagneticField. This is modeled after trace_particles().

In the case of mod='gc_vac' we solve the guiding center equations under the vacuum assumption, i.e \(G =\) const. and \(I = 0\):

\[ \begin{align}\begin{aligned}\dot s = -|B|_{,\theta} m(v_{||}^2/|B| + \mu)/(q \psi_0)\\\dot \theta = |B|_{,s} m(v_{||}^2/|B| + \mu)/(q \psi_0) + \iota v_{||} |B|/G\\\dot \zeta = v_{||}|B|/G\\\dot v_{||} = -(\iota |B|_{,\theta} + |B|_{,\zeta})\mu |B|/G,\end{aligned}\end{align} \]

where \(q\) is the charge, \(m\) is the mass, and \(v_\perp^2 = 2\mu|B|\).

In the case of mode='gc' we solve the general guiding center equations for an MHD equilibrium:

\[ \begin{align}\begin{aligned}\dot s = (I |B|_{,\zeta} - G |B|_{,\theta})m(v_{||}^2/|B| + \mu)/(\iota D \psi_0)\\\dot \theta = ((G |B|_{,\psi} - K |B|_{,\zeta}) m(v_{||}^2/|B| + \mu) - C v_{||} |B|)/(\iota D)\\\dot \zeta = (F v_{||} |B| - (|B|_{,\psi} I - |B|_{,\theta} K) m(\rho_{||}^2 |B| + \mu) )/(\iota D)\\\dot v_{||} = (C|B|_{,\theta} - F|B|_{,\zeta})\mu |B|/(\iota D)\\C = - m v_{||} K_{,\zeta}/|B| - q \iota + m v_{||}G'/|B|\\F = - m v_{||} K_{,\theta}/|B| + q + m v_{||}I'/|B|\\D = (F G - C I))/\iota\end{aligned}\end{align} \]

where primes indicate differentiation wrt \(\psi\). In the case mod='gc_noK', the above equations are used with \(K=0\).

Parameters:

• field – The BoozerMagneticField instance

• stz_inits – A (nparticles, 3) array with the initial positions of the particles in Boozer coordinates \((s,\theta,\zeta)\).

• parallel_speeds – A (nparticles, ) array containing the speed in direction of the B field for each particle.

• tmax – integration time

• mass – particle mass in kg, defaults to the mass of an alpha particle

• charge – charge in Coulomb, defaults to the charge of an alpha particle

• Ekin – kinetic energy in Joule, defaults to 3.52MeV

• tol – tolerance for the adaptive ode solver

• comm – MPI communicator to parallelize over

• zetas – list of angles in [0, 2pi] for which intersection with the plane corresponding to that zeta should be computed

• stopping_criteria – list of stopping criteria, mostly used in combination with the LevelsetStoppingCriterion accessed via simsopt.field.tracing.SurfaceClassifier.

• mode – how to trace the particles. Options are gc: general guiding center equations. gc_vac: simplified guiding center equations for the case \(G\) = const., \(I = 0\), and \(K = 0\). gc_noK: simplified guiding center equations for the case \(K = 0\).

• forget_exact_path – return only the first and last position of each particle for the res_tys. To be used when only res_zeta_hits is of interest or one wants to reduce memory usage.

Returns: 2 element tuple containing

• res_tys:

  A list of numpy arrays (one for each particle) describing the solution over time. The numpy array is of shape (ntimesteps, M) with M depending on the mode. Each row contains the time and the state. So for mode=’gc’ and mode=’gc_vac’ the state consists of the \((s,\theta,\zeta)\) position and the parallel speed, hence each row contains [t, s, t, z, v_par].

• res_zeta_hits:

  A list of numpy arrays (one for each particle) containing information on each time the particle hits one of the zeta planes or one of the stopping criteria. Each row of the array contains [time] + [idx] + state, where idx tells us which of the zetas or stopping_criteria was hit. If idx>=0, then zetas[int(idx)] was hit. If idx<0, then stopping_criteria[int(-idx)-1] was hit.

simsopt.field.tracing.trace_particles_starting_on_curve(curve, field, nparticles, tmax=0.0001, mass=6.69509884346e-27, charge=3.204353268e-19, Ekin=5.639661751679999e-13, tol=1e-09, comm=None, seed=1, umin=-1, umax=1, phis=[], stopping_criteria=[], mode='gc_vac', forget_exact_path=False, phase_angle=0)

Follows particles spawned at random locations on the magnetic axis with random pitch angle. See simsopt.field.tracing.trace_particles for the governing equations.

Parameters:

• curve – The simsopt.geo.curve.Curve to spawn the particles on. Uses rejection sampling to sample points on the curve. Warning: assumes that the underlying quadrature points on the Curve are uniformly distributed.

• field – The magnetic field \(B\).

• nparticles – number of particles to follow.

• tmax – integration time

• mass – particle mass in kg, defaults to the mass of an alpha particle

• charge – charge in Coulomb, defaults to the charge of an alpha particle

• Ekin – kinetic energy in Joule, defaults to 3.52MeV

• tol – tolerance for the adaptive ode solver

• comm – MPI communicator to parallelize over

• seed – random seed

• umin – the parallel speed is defined as v_par = u * speed_total where u is drawn uniformly in [umin, umax]

• umax – see umin

• phis – list of angles in [0, 2pi] for which intersection with the plane corresponding to that phi should be computed

• stopping_criteria – list of stopping criteria, mostly used in combination with the LevelsetStoppingCriterion accessed via simsopt.field.tracing.SurfaceClassifier.

• mode – how to trace the particles. options are gc: general guiding center equations, gc_vac: simplified guiding center equations for the case \(\nabla p=0\), full: full orbit calculation (slow!)

• forget_exact_path – return only the first and last position of each particle for the res_tys. To be used when only res_phi_hits is of interest or one wants to reduce memory usage.

• phase_angle – the phase angle to use in the case of full orbit calculations

Returns: see simsopt.field.tracing.trace_particles

simsopt.field.tracing.trace_particles_starting_on_surface(surface, field, nparticles, tmax=0.0001, mass=6.69509884346e-27, charge=3.204353268e-19, Ekin=5.639661751679999e-13, tol=1e-09, comm=None, seed=1, umin=-1, umax=1, phis=[], stopping_criteria=[], mode='gc_vac', forget_exact_path=False, phase_angle=0)

Follows particles spawned at random locations on the magnetic axis with random pitch angle. See simsopt.field.tracing.trace_particles for the governing equations.

Parameters:

• surface – The simsopt.geo.surface.Surface to spawn the particles on. Uses rejection sampling to sample points on the curve. Warning: assumes that the underlying quadrature points on the Curve as uniformly distributed.

• field – The magnetic field \(B\).

• nparticles – number of particles to follow.

• tmax – integration time

• mass – particle mass in kg, defaults to the mass of an alpha particle

• charge – charge in Coulomb, defaults to the charge of an alpha particle

• Ekin – kinetic energy in Joule, defaults to 3.52MeV

• tol – tolerance for the adaptive ode solver

• comm – MPI communicator to parallelize over

• seed – random seed

• umin – the parallel speed is defined as v_par = u * speed_total where u is drawn uniformly in [umin, umax].

• umax – see umin

• phis – list of angles in [0, 2pi] for which intersection with the plane corresponding to that phi should be computed

• stopping_criteria – list of stopping criteria, mostly used in combination with the LevelsetStoppingCriterion accessed via simsopt.field.tracing.SurfaceClassifier.

• mode – how to trace the particles. options are gc: general guiding center equations, gc_vac: simplified guiding center equations for the case \(\nabla p=0\), full: full orbit calculation (slow!)

• forget_exact_path – return only the first and last position of each particle for the res_tys. To be used when only res_phi_hits is of interest or one wants to reduce memory usage.

• phase_angle – the phase angle to use in the case of full orbit calculations

Returns: see simsopt.field.tracing.trace_particles

simsopt.field.wireframefield module

class simsopt.field.wireframefield.WireframeField(wframe)

Bases: WireframeField, MagneticField

Computes the magnetic field generated by currents in a wireframe grid.

Parameters:

wframe (instance of a ToroidalWireframe class)

dB_by_dsegmentcurrents(compute_derivatives)

Calculates the derivative of the magnetic field or its spatial derivatives at each field reference point with respect to the current in each wireframe segment.

Parameters:

compute_derivatives (integer (must be 0 or 1)) – If zero, will provide derivatives of the magnetic field with respect to segment currents. If one, will provide derivatives of the first spatial derivatives of the magnetic field with respect to segment currents.

Returns:

dB_by_dsegmentcurrents – Arrays of the requested derivatives, including one array for each segment in the wireframe.

Return type:

list of arrays

dBnormal_by_dsegmentcurrents_matrix(surface, area_weighted=False)

Generates a matrix with derivatives of normal magnetic field on a surface of interest with respect to the current in each wireframe segment, useful for optimization of the segment currents.

Parameters:

• surface (Surface class instance) – Surface on which to calculate the magnetic field. It is assumed that the quadrature points are evenly spaced in the toroidal and poloidal angles, and that they cover a whole number of half-periods.

• area_weighted (logical (optional)) – If true, will multiply each matrix element by the square root of the surface area ascribed to the corresponding quadrature point. In this way, the expression (A*x)^2 gives the surface integral of the squared flux, where A is the matrix with weighted elements and x is a vector with the currents in each wireframe segment. Default is False.

simsopt.field.wireframefield.enclosed_current(curve, field, n_quadpoints, preserve_points=True)

Calculates the integral of a vector field along a curve in 3D space. Useful for tests to verify consistency of wireframe solutions with their constraints.

Parameters:

• curve (CurveXYZFourier class instance) – Curve along which the integral is to be taken. The more quadrature points the curve has, the greater precision the result will have.

• field (MagneticField class instance) – Magnetic field in which the curve is integrated.

• n_quadpoints (integer) – Number of quadrature points for the integral.

• preserve_points (boolean) – If true, the existing field points of field will be restored before the function returns. If false, the field points will be changed to the quadrature points for integration. Default is true. Setting to false may save time.

Returns:

integral – Integral of the magnetic field along the curve.

Return type:

double