Magnetic fields
Simsopt provides several representations of magnetic fields, including
the Biot-Savart field associated with coils, as well as others.
Magnetic fields are represented as subclasses of the base class
simsopt.field.magneticfield.MagneticField. The various field
types can be scaled and summed together; for instance you can add a
purely toroidal field to the field from coils. Each field object has
an associated set of evaluation points. At these points, you can
query the field B, the vector potential A, and their derivatives with
respect to position or the field parameters.
Field types
Coils and BiotSavart
In simsopt, a filamentary coil is represented by the class
simsopt.field.coil.Coil. A Coil is a pairing of a
Curve with a current magnitude. The latter
is represented by a simsopt.field.coil.Current object. For
information about Curve objects see the page on geometric
objects. If you wish for several filamentary curves to carry
the same current, you can reuse a Current
object in multiple Coil objects.
Once Coil objects are defined, the magnetic
field they produce can be evaluated using the class
simsopt.field.biotsavart.BiotSavart. This class provides an
implementation of the Biot-Savart law
where \(\mu_0=4\pi \times 10^{-7}\) is the vacuum permitivity, \(\Gamma_k\) is the position vector of coil \(k\), and \(I_k\) indicates the electric currents.
Example:
import numpy as np
from simsopt.geo.curvexyzfourier import CurveXYZFourier
from simsopt.field.coil import Current, Coil
from simsopt.field.biotsavart import BiotSavart
curve = CurveXYZFourier(100, 1) # 100 = Number of quadrature points, 1 = max Fourier mode number
curve.x = [0, 0, 1., 0., 1., 0., 0., 0., 0.] # Set Fourier amplitudes
coil = Coil(curve, Current(1.0e4)) # 10 kAmpere-turns
field = BiotSavart([coil]) # Multiple coils can be included in the list
field.set_points(np.array([[0.5, 0.5, 0.1], [0.1, 0.1, -0.3]]))
print(field.B())
For a more complex example of a
BiotSavart object used in coil
optimization, see
examples/2_Intermediate/stage_two_optimization.py.
ToroidalField
The simsopt.field.magneticfieldclasses.ToroidalField class
represents a purely toroidal magnetic field. The field is given by
\(\mathbf B = B_0 \frac{R_0}{R} \mathbf e_\phi\), where
\(R_0\) and \(B_0\) are input scalar quantities, with
\(R_0\) representing a reference major radius, \(B_0\) the
magnetic field at \(R_0\), and \(R\) is the radial coordinate
of the cylindrical coordinate system \((R,Z,\phi)\). The vector
\(\mathbf e_\phi\) is a unit vector pointing in the direction of
increasing \(\phi\), with \(\phi\) the standard azimuthal
angle. Given Cartesian coordinates \((x,y,z)\), \(\mathbf e_\phi\)
is calculated as \(\mathbf e_\phi=-\sin \phi \mathbf e_x+\cos
\phi \mathbf e_y\) with \(\phi=\arctan(y/x)\).
PoloidalField
The simsopt.field.magneticfieldclasses.PoloidalField class
represents a poloidal magnetic field acording to the formula
\(\mathbf B = B_0 \frac{r}{q R_0} \mathbf e_\theta\), where
\(R_0, q\) and \(B_0\) are input scalar
quantities. \(R_0\) represents the major radius of the magnetic
axis, \(B_0\) the magnetic field at \(r=R_0 q\) and \(q\)
the safety factor associated with the sum of a poloidal magnetic field
and a toroidal magnetic field with major radius \(R_0\) and
magnetic field on-axis \(B_0\). \(r\) is the radial coordinate
of the simple toroidal coordinate system
\((r,\phi,\theta)\). Given a set of points \((x,y,z)\),
\(r\) is calculated as
\(r=\sqrt{(\sqrt{x^2+y^2}-R_0)^2+z^2}\). The vector \(\mathbf
e_\theta\) is a unit vector pointing in the direction of increasing
\(\theta\), with \(\theta\) the poloidal angle in the simple
toroidal coordinate system \((r,\phi,\theta)\). Given a set of
points \((x,y,z)\), \(\mathbf e_\theta\) is calculated as
\(\mathbf e_\theta=-\sin \theta \cos \phi \mathbf e_x+\sin \theta
\sin \phi \mathbf e_y+\cos \theta \mathbf e_z\) with
\(\phi=\arctan(y/x)\) and
\(\theta=\arctan(z/(\sqrt{x^2+y^2}-R_0))\).
ScalarPotentialRZMagneticField
The
simsopt.field.magneticfieldclasses.ScalarPotentialRZMagneticField
class initializes a vacuum magnetic field \(\mathbf B = \nabla
\Phi\) defined via a scalar potential \(\Phi\) in cylindrical
coordinates \((R,Z,\phi)\). The field \(\Phi\) is specified as
an analytical expression via a string argument. Simsopt performs the
necessary partial derivatives in order find \(\mathbf B\) and its
derivatives. For example, the function
ScalarPotentialRZMagneticField("2*phi") represents a toroidal
magnetic field \(\mathbf B = \nabla (2\phi)=2/R \mathbf e_\phi\).
Note: this functions needs the library sympy for the analytical
derivatives.
CircularCoil
The simsopt.field.magneticfieldclasses.CircularCoil class
represents a magnetic field created by a single circular coil. It
takes four input quantities: \(a\), the radius of the coil,
\(\mathbf c=[c_x,c_y,c_z]\), the center of the coil, \(I\),
the current flowing through the coil and \(\mathbf n\), the normal
vector to the plane of the coil centered at the coil radius, which
could be specified either with its three Cartesian components
\(\mathbf n=[n_x,n_y,n_z]\) or as \(\mathbf n=[\theta,\phi]\)
with the spherical angles \(\theta\) and \(\phi\).
The magnetic field is calculated analitically using the following expressions (reference)
\(B_x=\frac{\mu_0 I}{2\pi}\frac{x z}{\alpha^2 \beta \rho^2}\left[(a^2+r^2)E(k^2)-\alpha^2 K(k^2)\right]\)
\(B_y=\frac{y}{x}B_x\)
\(B_z=\frac{\mu_0 I}{2\pi \alpha^2 \beta}\left[(a^2-r^2)E(k^2)+\alpha^2 K(k^2)\right]\)
where \(\rho^2=x^2+y^2\), \(r^2=x^2+y^2+z^2\), \(\alpha^2=a^2+r^2-2a\rho\), \(\beta^2=a^2+r^2+2 a \rho\), \(k^2=1-\alpha^2/\beta^2\).
Dommaschk
The simsopt.field.magneticfieldclasses.Dommaschk class
represents a vacuum magnetic field \(\mathbf B = \nabla \Phi\)
with basis functions for the scalar potential \(\Phi\) described
in W. Dommaschk (1986), Computer Physics Communications 40, 203-218. This
representation provides explicit analytic formulae for vacuum fields
with a mixture of flux surfaces, islands, and chaos. Following the
original reference, a toroidal field with \(B_0=R_0=1\) is already
included in the definition. As input parameters, it takes two arrays:
The first array is an \(N\times2\) array \([(m_1,n_1),(m_2,n_2),...]\) specifying which harmonic coefficients \(m\) and \(n\) are non-zero.
The second array is an \(N\times2\) array \([(b_1,c_1),(b_2,c_2),...]\) with \(b_i=b_{m_i,n_i}\) and \(c_i=c_{m_i,n_i}\) the coefficients used in the Dommaschk representation.
Reiman
The simsopt.field.magneticfieldclasses.Reiman provides the
magnetic field model in section 5 of Reiman and Greenside, Computer
Physics Communications 43 (1986) 157—167.
It is an analytical magnetic field representation that allows the
explicit calculation of the width of the magnetic field islands.
InterpolatedField
The simsopt.field.magneticfieldclasses.InterpolatedField
function takes an existing field and interpolates it on a regular grid
in \(r,\phi,z\). This resulting interpolant can then be evaluated
very quickly. This is useful for efficiently tracing field lines and
particle trajectories.
Scaling and summing fields
Magnetic field objects can be added together, either by using the
+ operator, or by creating an instance of the class
simsopt.field.magneticfield.MagneticFieldSum. (The +
operator creates the latter.)
Magnetic fields can also be scaled by a constant. This can be accomplished either using the * operator,
or by creating an instance of the class
simsopt.field.magneticfield.MagneticFieldMultiply. (The *
operator creates the latter.)
Example:
from simsopt.field.magneticfieldclasses import ToroidalField, CircularCoil
field1 = CircularCoil(I=1.e7, r0=1.)
field2 = ToroidalField(R0=1., B0=1.)
total_field = field1 + 2.5 * field2
Common operations
Magnetic field objects have a large number of functions available. Before evaluating the field, you must
set the evaluation points. This can be done using either Cartesian or cylindrical coordinates.
Let m be a MagneticField object, and suppose there are n points
at which you wish to evaluate the field.
m.set_points_cart()takes a numpy array of size(n, 3)with the Cartesian coordinates(x, y, z)of the points.m.set_points_cyl()takes a numpy array of size(n, 3)with the cylindrical coordinates(r, phi, z)of the points.m.set_points()is shorthand form.set_points_cart().m.get_points_cart()returns a numpy array of size(n, 3)with the Cartesian coordinates(x, y, z)of the points.m.get_points_cyl()returns a numpy array of size(n, 3)with the cylindrical coordinates(r, phi, z)of the points.
A variety of functions are available to return the magnetic field \(B\), vector potential \(A\), and their gradients. The most commonly used ones are the following:
m.B()returns an array of size(n, 3)with the Cartesian coordinates of \(B\).m.B_cyl()returns an array of size(n, 3)with the cylindrical(r, phi, z)coordinates of \(B\).m.A()returns an array of size(n, 3)with the Cartesian coordinates of \(A\).m.AbsB()returns an array of size(n, 1)with the field magnitude \(|B|\).m.dB_by_dX()returns an array of size(n, 3, 3)with the Cartesian coordinates of \(\nabla B\). Denoting the indices by \((i,j,l)\), the result contains \(\partial_j B_l(x_i)\).m.d2B_by_dXdX()returns an array of size(n, 3, 3, 3)with the Cartesian coordinates of \(\nabla\nabla B\). Denoting the indices by \((i,j,k,l)\), the result contains \(\partial_k \partial_j B_l(x_i)\).m.dA_by_dX()returns an array of size(n, 3, 3)with the Cartesian coordinates of \(\nabla A\). Denoting the indices by \((i,j,l)\), the result contains \(\partial_j A_l(x_i)\).m.d2A_by_dXdX()returns an array of size(n, 3, 3, 3)with the Cartesian coordinates of \(\nabla\nabla A\). Denoting the indices by \((i,j,k,l)\), the result contains \(\partial_k \partial_j A_l(x_i)\).m.GradAbsB()returns an array of size(n, 3)with the Cartesian components of \(\nabla |B|\).
Example:
import numpy as np
from simsopt.field.magneticfieldclasses import CircularCoil
field = CircularCoil(I=1.e7, r0=1.)
points = np.array([[0.5, 0.5, 0.1], [0.1, 0.1, -0.3]])
field.set_points(points)
print(field.B())
print(field.dB_by_dX())