simsopt.field package
- class simsopt.field.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.
- class simsopt.field.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(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 byBoozerRadialInterpolant
andInterpolatedBoozerField
. The usage is similar to theMagneticField
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 theclear_cached_properties()
function. The cache is automatically cleared whenset_points()
is called or one of the dependencies changes.- clear_cached_properties()
Clear the cache.
- recompute_bell(parent=None)
- class simsopt.field.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 usingBOOZXFORM
. 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
orsimsopt.mhd.boozer.Boozer
. If it is an instance ofsimsopt.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 byns_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
. IfFalse
, interpolation of the unscaled Fourier harmonics and its finite-difference derivative wrts
is performed instead (defaults toFalse
)ns_delete – (see
rescale
) (defaults to 0)
- compute_K()
- init_splines()
- class simsopt.field.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.
- class simsopt.field.Coil(curve, current)
Bases:
Coil
,Optimizable
A
Coil
combines aCurve
and aCurrent
and is used as input for aBiotSavart
field.- 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 tosimsopt.geo.curve.Curve.plot()
- vjp(v_gamma, v_gammadash, v_current)
- class simsopt.field.Current(current, dofs=None, **kwargs)
Bases:
Current
,CurrentBase
An optimizable object that wraps around a single scalar degree of freedom. It represents the electric current in a coil, or in a set of coils that are constrained to use the same current.
- property current
- vjp(v_current)
- class simsopt.field.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.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 afterInterpolatedField
.- __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 intons
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
andzetamax=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 usethetamin=0
andthetamax=np.pi
.
- class simsopt.field.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 intonr
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
andphimax=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 usezmin=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.IterationStoppingCriterion
Bases:
IterationStoppingCriterion
Stop the iteration once the maximum number of iterations is reached.
- class simsopt.field.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.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 whenset_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_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.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.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.
- class simsopt.field.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.MinToroidalFluxStoppingCriterion
Bases:
MinToroidalFluxStoppingCriterion
Stop the iteration once a particle falls below a critical value of
s
, the normalized toroidal flux. ThisStoppingCriterion
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.NormalField(nfp=1, stellsym=True, mpol=1, ntor=0, vns=None, vnc=None)
Bases:
Optimizable
NormalField
represents the magnetic field normal to a toroidal surface, for example the computational boundary of SPEC free-boundary.- 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
- 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_vnc(m, n)
- get_vns(m, n)
- set_vnc(m, n, value)
- set_vns(m, n, value)
- class simsopt.field.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.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.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.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.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.
- class simsopt.field.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 andflux
is a boolean indicating whether tracing is being performed in a flux coordinate system.
- simsopt.field.apply_symmetries_to_currents(base_currents, nfp, stellsym)
Take a list of
n
Current`s and return ``n * nfp * (1+int(stellsym))`
Current
objects obtained by copying (fornfp
rotations) and sign-flipping (optionally for stellarator symmetry).
- simsopt.field.apply_symmetries_to_curves(base_curves, nfp, stellsym)
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 tonfp
fold rotational symmetry and optionally stellarator symmetry.
- simsopt.field.coils_to_focus(filename, curves, currents, nfp=1, stellsym=False, Ifree=False, Lfree=False)
Export a list of Curve objects together with currents in FOCUS format, so they can be used by FOCUS. The format is introduced at https://princetonuniversity.github.io/FOCUS/rdcoils.pdf This routine only works with curves of type CurveXYZFourier, not other curve types.
- Parameters:
filename – Name of the file to write.
curves – A python list of CurveXYZFourier objects.
currents – Coil current of each curve.
nfp – The number of field periodicity. Defaults to 1.
stellsym – Whether or not following stellarator symmetry. Defaults to False.
Ifree – Flag specifying whether the coil current is free. Defaults to False.
Lfree – Flag specifying whether the coil geometry is free. Defaults to False.
- simsopt.field.coils_to_makegrid(filename, curves, currents, groups=None, nfp=1, stellsym=False)
Export a list of Curve objects together with currents in MAKEGRID input format, so they can be used by MAKEGRID and FOCUS. The format is introduced at https://princetonuniversity.github.io/STELLOPT/MAKEGRID Note that this function does not generate files with MAKEGRID’s output format.
- Parameters:
filename – Name of the file to write.
curves – A python list of Curve objects.
currents – Coil current of each curve.
groups – Coil current group. Coils in the same group will be assembled together. Defaults to None.
nfp – The number of field periodicity. Defaults to 1.
stellsym – Whether or not following stellarator symmetry. Defaults to False.
- simsopt.field.coils_via_symmetries(curves, currents, nfp, stellsym)
Take a list of
n
curves and returnn * nfp * (1+int(stellsym))
Coil
objects obtained by applying rotations and flipping corresponding tonfp
fold rotational symmetry and optionally stellarator symmetry.
- simsopt.field.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)\]- 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 devicetol – 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 viasimsopt.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.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 aBoozerMagneticField
), 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()
ortrace_particles_boozer()
withforget_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 beNone
.flux – if
True
,res_tys
represents the position in flux coordinates (should beTrue
if computed fromtrace_particles_boozer()
). IfTrue
,ma
is not used.
- Returns:
- array with length
len(res_tys)
. Each element contains the number of poloidal transits of the orbit.
- array with length
- Return type:
ntransits
- simsopt.field.compute_resonances(res_tys, res_phi_hits, ma=None, delta=0.01)
Computes resonant particle orbits given the output of either
trace_particles()
ortrace_particles_boozer()
,res_tys
andres_phi_hits
/res_zeta_hits
, withforget_exact_path=False
. Resonance indicates a trajectory which returns to the same position at the \(\zeta = 0\) plane aftermpol
poloidal turns andntor
toroidal turns. For the case of particles traced in aMagneticField
(not aBoozerMagneticField
), 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()
ortrace_particles_boozer()
withforget_exact_path=False
res_phi_hits – output of
trace_particles()
ortrace_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 beNone
. (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]
ifma=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, andntor
is the number of toroidal turns.
- Return type:
resonances
- simsopt.field.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()
ortrace_particles_boozer()
withforget_exact_path=False
.flux – if
True
,res_tys
represents the position in flux coordinates (should beTrue
if computed fromtrace_particles_boozer()
)
- Returns:
- array with length
len(res_tys)
. Each element contains the number of toroidal transits of the orbit.
- array with length
- Return type:
ntransits
- simsopt.field.load_coils_from_makegrid_file(filename, order, ppp=20)
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 – file to load.
order – maximum mode number in the Fourier expansion.
ppp – points-per-period: number of quadrature points per period.
- Returns:
A list of
Coil
objects with the Fourier coefficients and currents given by the file.
- simsopt.field.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.plot_poincare_data(fieldlines_phi_hits, phis, filename, mark_lost=False, aspect='equal', dpi=300)
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.trace_particles(field: MagneticField, xyz_inits: Sequence[Real] | NDArray[Any, float64], parallel_speeds: Sequence[Real] | NDArray[Any, 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 viasimsopt.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.trace_particles_boozer(field: BoozerMagneticField, stz_inits: Sequence[Real] | NDArray[Any, float64], parallel_speeds: Sequence[Real] | NDArray[Any, 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 aftertrace_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
instancestz_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 viasimsopt.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.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
whereu
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 viasimsopt.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.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
whereu
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 viasimsopt.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